8 int mma_verbose = 0; /* > 0 for verbose output */
10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
13 /* magic minimum value for rho in MMA ... the 2002 paper says it should
14 be a "fixed, strictly positive `small' number, e.g. 1e-5"
15 ... grrr, I hate these magic numbers, which seem like they
16 should depend on the objective function in some way ... in particular,
17 note that rho is dimensionful (= dimensions of objective function) */
18 #define MMA_RHOMIN 1e-5
20 /***********************************************************************/
21 /* function for MMA's dual solution of the approximate problem */
24 int count; /* evaluation count, incremented each call */
25 int n; /* must be set on input to dimension of x */
26 const double *x, *lb, *ub, *sigma, *dfdx; /* arrays of length n */
27 const double *dfcdx; /* m-by-n array of fc gradients */
28 double fval, rho; /* must be set on input */
29 const double *fcval, *rhoc; /* arrays of length m */
30 double *xcur; /* array of length n, output each time */
31 double gval, wval, *gcval; /* output each time (array length m) */
34 static double sqr(double x) { return x * x; }
36 static double dual_func(int m, const double *y, double *grad, void *d_)
38 dual_data *d = (dual_data *) d_;
40 const double *x = d->x, *lb = d->lb, *ub = d->ub, *sigma = d->sigma,
42 const double *dfcdx = d->dfcdx;
43 double rho = d->rho, fval = d->fval;
44 const double *rhoc = d->rhoc, *fcval = d->fcval;
45 double *xcur = d->xcur;
46 double *gcval = d->gcval;
54 for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]);
56 for (j = 0; j < n; ++j) {
57 double u, v, dx, denominv, c, sigma2, dx2;
59 /* first, compute xcur[j] for y. Because this objective is
60 separable, we can minimize over x analytically, and the minimum
61 dx is given by the solution of a quadratic equation:
62 u dx^2 + 2 v sigma^2 dx + u sigma^2 = 0
63 where u and v are defined by the sums below. Because of
64 the definitions, it is guaranteed that |u/v| <= sigma,
65 and it follows that the only dx solution with |dx| <= sigma
67 (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2))
68 (which goes to zero as u -> 0). */
71 v = fabs(dfdx[j]) * sigma[j] + 0.5 * rho;
72 for (i = 0; i < m; ++i) {
73 u += dfcdx[i*n + j] * y[i];
74 v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
76 u *= (sigma2 = sqr(sigma[j]));
77 if (fabs(u) < 1e-3 * (v*sigma[j])) { /* Taylor exp. for small u */
78 double a = u / (v*sigma[j]);
79 dx = -sigma[j] * (0.5 * a + 0.125 * a*a*a);
82 dx = (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j]))));
84 if (xcur[j] > ub[j]) xcur[j] = ub[j];
85 else if (xcur[j] < lb[j]) xcur[j] = lb[j];
86 if (xcur[j] > x[j]+0.9*sigma[j]) xcur[j] = x[j]+0.9*sigma[j];
87 else if (xcur[j] < x[j]-0.9*sigma[j]) xcur[j] = x[j]-0.9*sigma[j];
92 denominv = 1.0 / (sigma2 - dx2);
93 val += (u * dx + v * dx2) * denominv;
95 /* update gval, wval, gcval (approximant functions) */
97 d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2)
99 d->wval += 0.5 * dx2 * denominv;
100 for (i = 0; i < m; ++i)
101 gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j]
102 + 0.5*rhoc[j]) * dx2)
106 /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
107 we only need the partial derivative with respect to y, and
108 we negate because we are maximizing: */
109 if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
113 /***********************************************************************/
115 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
116 int m, nlopt_func fc,
117 void *fc_data_, ptrdiff_t fc_datum_size,
118 const double *lb, const double *ub, /* bounds */
119 double *x, /* in: initial guess, out: minimizer */
121 nlopt_stopping *stop,
122 nlopt_algorithm dual_alg,
123 double dual_tolrel, int dual_maxeval)
125 nlopt_result ret = NLOPT_SUCCESS;
126 double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
127 double *dfcdx, *dfcdx_cur;
128 double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
130 char *fc_data = (char *) fc_data_;
134 sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
135 if (!sigma) return NLOPT_OUT_OF_MEMORY;
140 xprevprev = xprev + n;
141 fcval = xprevprev + n;
142 fcval_cur = fcval + m;
143 rhoc = fcval_cur + m;
146 dual_ub = dual_lb + m;
149 dfcdx_cur = dfcdx + m*n;
163 for (j = 0; j < n; ++j) {
164 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
165 sigma[j] = 1.0; /* arbitrary default */
167 sigma[j] = 0.5 * (ub[j] - lb[j]);
170 for (i = 0; i < m; ++i) {
172 dual_lb[i] = y[i] = 0.0;
173 dual_ub[i] = HUGE_VAL;
176 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
178 memcpy(xcur, x, sizeof(double) * n);
181 for (i = 0; i < m; ++i) {
182 fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
183 feasible = feasible && (fcval[i] <= 0);
185 if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */
187 while (1) { /* outer iterations */
189 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
190 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
191 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
192 if (ret != NLOPT_SUCCESS) goto done;
193 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
194 memcpy(xprev, xcur, sizeof(double) * n);
196 while (1) { /* inner iterations */
198 int feasible_cur, inner_done;
201 /* solve dual problem */
202 dd.rho = rho; dd.count = 0;
203 reti = nlopt_minimize(
204 dual_alg, m, dual_func, &dd,
205 dual_lb, dual_ub, y, &min_dual,
206 -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
207 stop->maxtime - (nlopt_seconds() - stop->start));
208 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
213 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
215 printf("MMA dual converged in %d iterations to g=%g:\n",
217 for (i = 0; i < MIN(mma_verbose, m); ++i)
218 printf(" MMA y[%d]=%g, gc[%d]=%g\n",
219 i, y[i], i, dd.gcval[i]);
222 fcur = f(n, xcur, dfdx_cur, f_data);
225 inner_done = dd.gval >= fcur;
226 for (i = 0; i < m; ++i) {
227 fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
228 fc_data + fc_datum_size * i);
229 feasible_cur = feasible_cur && (fcval[i] <= 0);
230 inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
233 /* once we have reached a feasible solution, the
234 algorithm should never make the solution infeasible
235 again, although the constraints may be violated
236 slightly by rounding errors etc. so we must be a
237 little careful about checking feasibility */
238 if (feasible_cur) feasible = 1;
241 dd.fval = *minf = fcur;
242 memcpy(fcval, fcval_cur, sizeof(double)*m);
243 memcpy(x, xcur, sizeof(double)*n);
244 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
245 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
247 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
248 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
249 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
250 if (ret != NLOPT_SUCCESS) goto done;
252 if (inner_done) break;
255 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
256 for (i = 0; i < m; ++i)
257 if (fcval_cur[i] > dd.gcval[i])
260 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
264 printf("MMA inner iteration: rho -> %g\n", rho);
265 for (i = 0; i < MIN(mma_verbose, m); ++i)
266 printf(" rhoc[%d] -> %g\n", i,rhoc[i]);
269 if (nlopt_stop_ftol(stop, fcur, fprev))
270 ret = NLOPT_FTOL_REACHED;
271 if (nlopt_stop_x(stop, xcur, xprev))
272 ret = NLOPT_XTOL_REACHED;
273 if (ret != NLOPT_SUCCESS) goto done;
275 /* update rho and sigma for iteration k+1 */
276 rho = MAX(0.1 * rho, MMA_RHOMIN);
278 printf("MMA outer iteration: rho -> %g\n", rho);
279 for (i = 0; i < m; ++i)
280 rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
281 for (i = 0; i < MIN(mma_verbose, m); ++i)
282 printf(" rhoc[%d] -> %g\n", i, rhoc[i]);
284 for (j = 0; j < n; ++j) {
285 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
286 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
288 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
289 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
290 sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
293 for (j = 0; j < MIN(mma_verbose, n); ++j)
294 printf(" sigma[%d] -> %g\n",