8 int auglag_verbose = 1;
10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
13 /***************************************************************************/
16 nlopt_func f; void *f_data;
17 int m; nlopt_constraint *fc;
18 int p; nlopt_constraint *h;
19 double rho, *lambda, *mu;
24 /* the augmented lagrangian objective function */
25 static double auglag(int n, const double *x, double *grad, void *data)
27 auglag_data *d = (auglag_data *) data;
28 double *gradtmp = grad ? d->gradtmp : NULL;
30 const double *lambda = d->lambda, *mu = d->mu;
34 L = d->f(n, x, grad, d->f_data);
36 for (i = 0; i < d->p; ++i) {
38 h = d->h[i].f(n, x, gradtmp, d->h[i].f_data) + lambda[i] / rho;
40 if (grad) for (j = 0; j < n; ++j) grad[j] += (rho * h) * gradtmp[j];
43 for (i = 0; i < d->m; ++i) {
45 fc = d->fc[i].f(n, x, gradtmp, d->fc[i].f_data) + mu[i] / rho;
47 L += 0.5 * rho * fc*fc;
48 if (grad) for (j = 0; j < n; ++j)
49 grad[j] += (rho * fc) * gradtmp[j];
58 /***************************************************************************/
60 nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
61 int m, nlopt_constraint *fc,
62 int p, nlopt_constraint *h,
63 const double *lb, const double *ub, /* bounds */
64 double *x, /* in: initial guess, out: minimizer */
67 nlopt_opt sub_opt, int sub_has_fc)
70 nlopt_result ret = NLOPT_SUCCESS;
71 double ICM = HUGE_VAL, minf_penalty = HUGE_VAL, penalty;
72 double *xcur = NULL, fcur;
73 int i, feasible, minf_feasible = 0;
75 /* magic parameters from Birgin & Martinez */
76 const double tau = 0.5, gam = 10;
77 const double lam_min = -1e20, lam_max = 1e20, mu_max = 1e20;
79 d.f = f; d.f_data = f_data;
84 /* whether we handle inequality constraints via the augmented
85 Lagrangian penalty function, or directly in the sub-algorithm */
91 ret = nlopt_set_min_objective(sub_opt, auglag, &d); if (ret<0) return ret;
92 ret = nlopt_set_lower_bounds(sub_opt, lb); if (ret<0) return ret;
93 ret = nlopt_set_upper_bounds(sub_opt, ub); if (ret<0) return ret;
94 ret = nlopt_remove_inequality_constraints(sub_opt); if (ret<0) return ret;
95 ret = nlopt_remove_equality_constraints(sub_opt); if (ret<0) return ret;
96 for (i = 0; i < m; ++i) {
97 ret = nlopt_add_inequality_constraint(sub_opt, fc[i].f, fc[i].f_data,
99 if (ret < 0) return ret;
102 xcur = (double *) malloc(sizeof(double) * (2*n + d.p + d.m));
103 if (!xcur) return NLOPT_OUT_OF_MEMORY;
104 memcpy(xcur, x, sizeof(double) * n);
106 d.gradtmp = xcur + n;
107 memset(d.gradtmp, 0, sizeof(double) * (n + d.p + d.m));
108 d.lambda = d.gradtmp + n;
109 d.mu = d.lambda + d.p;
113 /* starting rho suggested by B & M */
114 if (d.p > 0 || d.m > 0) {
117 fcur = f(n, xcur, NULL, f_data);
120 for (i = 0; i < d.p; ++i) {
121 double hi = h[i].f(n, xcur, NULL, d.h[i].f_data);
123 feasible = feasible && fabs(hi) <= h[i].tol;
126 for (i = 0; i < d.m; ++i) {
127 double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data);
128 penalty += fci > 0 ? fci : 0;
129 feasible = feasible && fci <= fc[i].tol;
130 if (fci > 0) con2 += fci * fci;
133 minf_penalty = penalty;
134 minf_feasible = feasible;
135 d.rho = MAX(1e-6, MIN(10, 2 * fabs(*minf) / con2));
138 d.rho = 1; /* whatever, doesn't matter */
140 if (auglag_verbose) {
141 printf("auglag: initial rho=%g\nauglag initial lambda=", d.rho);
142 for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]);
143 printf("\nauglag initial mu = ");
144 for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]);
149 double prev_ICM = ICM;
151 ret = nlopt_optimize_limited(sub_opt, xcur, &fcur,
152 stop->maxeval - stop->nevals,
153 stop->maxtime - (nlopt_seconds()
158 fcur = f(n, xcur, NULL, f_data);
163 for (i = 0; i < d.p; ++i) {
164 double hi = h[i].f(n, xcur, NULL, d.h[i].f_data);
165 double newlam = d.lambda[i] + d.rho * hi;
167 feasible = feasible && fabs(hi) <= h[i].tol;
168 ICM = MAX(ICM, fabs(hi));
169 d.lambda[i] = MIN(MAX(lam_min, newlam), lam_max);
171 for (i = 0; i < d.m; ++i) {
172 double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data);
173 double newmu = d.mu[i] + d.rho * fci;
174 penalty += fci > 0 ? fci : 0;
175 feasible = feasible && fci <= fc[i].tol;
176 ICM = MAX(ICM, fabs(MAX(fci, -d.mu[i] / d.rho)));
177 d.mu[i] = MIN(MAX(0.0, newmu), mu_max);
179 if (ICM > tau * prev_ICM) {
183 if (auglag_verbose) {
184 printf("auglag: ICM=%g, rho=%g\nauglag lambda=", ICM, d.rho);
185 for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]);
186 printf("\nauglag mu = ");
187 for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]);
191 if ((feasible && (!minf_feasible || penalty < minf_penalty
192 || fcur <= *minf)) ||
193 (!minf_feasible && penalty <= minf_penalty)) {
196 if (fcur < stop->minf_max)
197 ret = NLOPT_MINF_MAX_REACHED;
198 else if (nlopt_stop_ftol(stop, fcur, *minf))
199 ret = NLOPT_FTOL_REACHED;
200 else if (nlopt_stop_x(stop, xcur, x))
201 ret = NLOPT_XTOL_REACHED;
203 else { /* check if no progress towards feasibility */
204 if (nlopt_stop_ftol(stop, fcur, *minf)
205 && nlopt_stop_ftol(stop, penalty, minf_penalty))
206 ret = NLOPT_FTOL_REACHED;
207 else if (nlopt_stop_x(stop, xcur, x))
208 ret = NLOPT_XTOL_REACHED;
211 minf_penalty = penalty;
212 minf_feasible = feasible;
213 memcpy(x, xcur, sizeof(double) * n);
214 if (ret != NLOPT_SUCCESS) break;
217 if (nlopt_stop_evals(stop)) {ret = NLOPT_MAXEVAL_REACHED; break;}
218 if (nlopt_stop_time(stop)) {ret = NLOPT_MAXTIME_REACHED; break;}
220 /* TODO: use some other stopping criterion on ICM? */
221 /* The paper uses ICM <= epsilon and DFM <= epsilon, where
222 DFM is a measure of the size of the Lagrangian gradient.
223 Besides the fact that these kinds of absolute tolerances
224 (non-scale-invariant) are unsatisfying and it is not
225 clear how the user should specify it, the ICM <= epsilon
226 condition seems not too different from requiring feasibility. */
227 if (ICM == 0) return NLOPT_FTOL_REACHED;