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 dx = u==0 ? 0 : (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j]))));
79 if (xcur[j] > ub[j]) xcur[j] = ub[j];
80 else if (xcur[j] < lb[j]) xcur[j] = lb[j];
81 if (xcur[j] > x[j]+0.9*sigma[j]) xcur[j] = x[j]+0.9*sigma[j];
82 else if (xcur[j] < x[j]-0.9*sigma[j]) xcur[j] = x[j]-0.9*sigma[j];
87 denominv = 1.0 / (sigma2 - dx2);
88 val += (u * dx + v * dx2) * denominv;
90 /* update gval, wval, gcval (approximant functions) */
92 d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2)
94 d->wval += 0.5 * dx2 * denominv;
95 for (i = 0; i < m; ++i)
96 gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j]
101 /* gradient is easy to compute: since we are at a minimum x (dval/dx=0),
102 we only need the partial derivative with respect to y, and
103 we negate because we are maximizing: */
104 if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i];
108 /***********************************************************************/
110 nlopt_result mma_minimize(int n, nlopt_func f, void *f_data,
111 int m, nlopt_func fc,
112 void *fc_data_, ptrdiff_t fc_datum_size,
113 const double *lb, const double *ub, /* bounds */
114 double *x, /* in: initial guess, out: minimizer */
116 nlopt_stopping *stop,
117 nlopt_algorithm dual_alg,
118 double dual_tolrel, int dual_maxeval)
120 nlopt_result ret = NLOPT_SUCCESS;
121 double *xcur, rho, *sigma, *dfdx, *dfdx_cur, *xprev, *xprevprev, fcur;
122 double *dfcdx, *dfcdx_cur;
123 double *fcval, *fcval_cur, *rhoc, *gcval, *y, *dual_lb, *dual_ub;
125 char *fc_data = (char *) fc_data_;
129 sigma = (double *) malloc(sizeof(double) * (6*n + 2*m*n + m*7));
130 if (!sigma) return NLOPT_OUT_OF_MEMORY;
135 xprevprev = xprev + n;
136 fcval = xprevprev + n;
137 fcval_cur = fcval + m;
138 rhoc = fcval_cur + m;
141 dual_ub = dual_lb + m;
144 dfcdx_cur = dfcdx + m*n;
158 for (j = 0; j < n; ++j) {
159 if (nlopt_isinf(ub[j]) || nlopt_isinf(lb[j]))
160 sigma[j] = 1.0; /* arbitrary default */
162 sigma[j] = 0.5 * (ub[j] - lb[j]);
165 for (i = 0; i < m; ++i) {
167 dual_lb[i] = y[i] = 0.0;
168 dual_ub[i] = HUGE_VAL;
171 dd.fval = fcur = *minf = f(n, x, dfdx, f_data);
173 memcpy(xcur, x, sizeof(double) * n);
176 for (i = 0; i < m; ++i) {
177 fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i);
178 feasible = feasible && (fcval[i] <= 0);
180 if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */
182 while (1) { /* outer iterations */
184 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
185 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
186 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
187 if (ret != NLOPT_SUCCESS) goto done;
188 if (++k > 1) memcpy(xprevprev, xprev, sizeof(double) * n);
189 memcpy(xprev, xcur, sizeof(double) * n);
191 while (1) { /* inner iterations */
193 int feasible_cur, inner_done;
196 /* solve dual problem */
197 dd.rho = rho; dd.count = 0;
198 reti = nlopt_minimize(
199 dual_alg, m, dual_func, &dd,
200 dual_lb, dual_ub, y, &min_dual,
201 -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval,
202 stop->maxtime - (nlopt_seconds() - stop->start));
203 if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) {
208 dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */
210 printf("MMA dual converged in %d iterations to g=%g:\n",
212 for (i = 0; i < mma_verbose; ++i)
213 printf(" MMA y[%d]=%g, gc[%d]=%g\n",
214 i, y[i], i, dd.gcval[i]);
217 fcur = f(n, xcur, dfdx_cur, f_data);
220 inner_done = dd.gval >= fcur;
221 for (i = 0; i < m; ++i) {
222 fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n,
223 fc_data + fc_datum_size * i);
224 feasible_cur = feasible_cur && (fcval[i] <= 0);
225 inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]);
228 /* once we have reached a feasible solution, the
229 algorithm should never make the solution infeasible
230 again, although the constraints may be violated
231 slightly by rounding errors etc. so we must be a
232 little careful about checking feasibility */
233 if (feasible_cur) feasible = 1;
236 dd.fval = *minf = fcur;
237 memcpy(fcval, fcval_cur, sizeof(double)*m);
238 memcpy(x, xcur, sizeof(double)*n);
239 memcpy(dfdx, dfdx_cur, sizeof(double)*n);
240 memcpy(dfcdx, dfcdx_cur, sizeof(double)*n*m);
242 if (nlopt_stop_evals(stop)) ret = NLOPT_MAXEVAL_REACHED;
243 else if (nlopt_stop_time(stop)) ret = NLOPT_MAXTIME_REACHED;
244 else if (*minf < stop->minf_max) ret = NLOPT_MINF_MAX_REACHED;
245 if (ret != NLOPT_SUCCESS) goto done;
247 if (inner_done) break;
250 rho = MIN(10*rho, 1.1 * (rho + (fcur-dd.gval) / dd.wval));
251 for (i = 0; i < m; ++i)
252 if (fcval_cur[i] > dd.gcval[i])
255 1.1 * (rhoc[i] + (fcval_cur[i]-dd.gcval[i])
259 printf("MMA inner iteration: rho -> %g\n", rho);
260 for (i = 0; i < mma_verbose; ++i)
261 printf(" rhoc[%d] -> %g\n", i,rhoc[i]);
264 if (nlopt_stop_ftol(stop, fcur, fprev))
265 ret = NLOPT_FTOL_REACHED;
266 if (nlopt_stop_x(stop, xcur, xprev))
267 ret = NLOPT_XTOL_REACHED;
268 if (ret != NLOPT_SUCCESS) goto done;
270 /* update rho and sigma for iteration k+1 */
271 rho = MAX(0.1 * rho, MMA_RHOMIN);
273 printf("MMA outer iteration: rho -> %g\n", rho);
274 for (i = 0; i < m; ++i)
275 rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN);
276 for (i = 0; i < mma_verbose; ++i)
277 printf(" rhoc[%d] -> %g\n", i, rhoc[i]);
279 for (j = 0; j < n; ++j) {
280 double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]);
281 double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1);
283 if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) {
284 sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j]));
285 sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j]));
288 for (j = 0; j < mma_verbose; ++j)
289 printf(" sigma[%d] -> %g\n",