auglag/auglag.c

   1 #include <stdlib.h>
   2 #include <math.h>
   3 #include <string.h>
   4 #include <stdio.h>
   5
   6 #include "auglag.h"
   7
   8 int auglag_verbose = 1;
   9
  10 #define MIN(a,b) ((a) < (b) ? (a) : (b))
  11 #define MAX(a,b) ((a) > (b) ? (a) : (b))
  12
  13 /***************************************************************************/
  14
  15 typedef struct {
  16      nlopt_func f; void *f_data;
  17      int m; nlopt_constraint *fc;
  18      int p; nlopt_constraint *h;
  19      double rho, *lambda, *mu;
  20      double *gradtmp;
  21      nlopt_stopping *stop;
  22 } auglag_data;
  23
  24 /* the augmented lagrangian objective function */
  25 static double auglag(unsigned n, const double *x, double *grad, void *data)
  26 {
  27      auglag_data *d = (auglag_data *) data;
  28      double *gradtmp = grad ? d->gradtmp : NULL;
  29      double rho = d->rho;
  30      const double *lambda = d->lambda, *mu = d->mu;
  31      double L;
  32      int i;
  33      unsigned j;
  34
  35      L = d->f(n, x, grad, d->f_data);
  36
  37      for (i = 0; i < d->p; ++i) {
  38           double h;
  39           h = d->h[i].f(n, x, gradtmp, d->h[i].f_data) + lambda[i] / rho;
  40           L += 0.5 * rho * h*h;
  41           if (grad) for (j = 0; j < n; ++j) grad[j] += (rho * h) * gradtmp[j];
  42      }
  43
  44      for (i = 0; i < d->m; ++i) {
  45           double fc;
  46           fc = d->fc[i].f(n, x, gradtmp, d->fc[i].f_data) + mu[i] / rho;
  47           if (fc > 0) {
  48                L += 0.5 * rho * fc*fc;
  49                if (grad) for (j = 0; j < n; ++j)
  50                     grad[j] += (rho * fc) * gradtmp[j];
  51           }
  52      }
  53
  54      d->stop->nevals++;
  55
  56      return L;
  57 }
  58
  59 /***************************************************************************/
  60
  61 nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
  62                              int m, nlopt_constraint *fc,
  63                              int p, nlopt_constraint *h,
  64                              const double *lb, const double *ub, /* bounds */
  65                              double *x, /* in: initial guess, out: minimizer */
  66                              double *minf,
  67                              nlopt_stopping *stop,
  68                              nlopt_opt sub_opt, int sub_has_fc)
  69 {
  70      auglag_data d;
  71      nlopt_result ret = NLOPT_SUCCESS;
  72      double ICM = HUGE_VAL, minf_penalty = HUGE_VAL, penalty;
  73      double *xcur = NULL, fcur;
  74      int i, feasible, minf_feasible = 0;
  75      int auglag_iters = 0;
  76
  77      /* magic parameters from Birgin & Martinez */
  78      const double tau = 0.5, gam = 10;
  79      const double lam_min = -1e20, lam_max = 1e20, mu_max = 1e20;
  80
  81      d.f = f; d.f_data = f_data;
  82      d.m = m; d.fc = fc;
  83      d.p = p; d.h = h;
  84      d.stop = stop;
  85
  86      /* whether we handle inequality constraints via the augmented
  87         Lagrangian penalty function, or directly in the sub-algorithm */
  88      if (sub_has_fc)
  89           d.m = 0;
  90      else
  91           m = 0;
  92
  93      ret = nlopt_set_min_objective(sub_opt, auglag, &d); if (ret<0) return ret;
  94      ret = nlopt_set_lower_bounds(sub_opt, lb); if (ret<0) return ret;
  95      ret = nlopt_set_upper_bounds(sub_opt, ub); if (ret<0) return ret;
  96      ret = nlopt_set_stopval(sub_opt, stop->minf_max); if (ret<0) return ret;
  97      ret = nlopt_remove_inequality_constraints(sub_opt); if (ret<0) return ret;
  98      ret = nlopt_remove_equality_constraints(sub_opt); if (ret<0) return ret;
  99      for (i = 0; i < m; ++i) {
 100           ret = nlopt_add_inequality_constraint(sub_opt, fc[i].f, fc[i].f_data,
 101                                                 fc[i].tol);
 102           if (ret < 0) return ret;
 103      }
 104
 105      xcur = (double *) malloc(sizeof(double) * (2*n + d.p + d.m));
 106      if (!xcur) return NLOPT_OUT_OF_MEMORY;
 107      memcpy(xcur, x, sizeof(double) * n);
 108
 109      d.gradtmp = xcur + n;
 110      memset(d.gradtmp, 0, sizeof(double) * (n + d.p + d.m));
 111      d.lambda = d.gradtmp + n;
 112      d.mu = d.lambda + d.p;
 113
 114      *minf = HUGE_VAL;
 115
 116      /* starting rho suggested by B & M */
 117      if (d.p > 0 || d.m > 0) {
 118           double con2 = 0;
 119           d.stop->nevals++;
 120           fcur = f(n, xcur, NULL, f_data);
 121           penalty = 0;
 122           feasible = 1;
 123           for (i = 0; i < d.p; ++i) {
 124                double hi = h[i].f(n, xcur, NULL, d.h[i].f_data);
 125                penalty += fabs(hi);
 126                feasible = feasible && fabs(hi) <= h[i].tol;
 127                con2 += hi * hi;
 128           }
 129           for (i = 0; i < d.m; ++i) {
 130                double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data);
 131                penalty += fci > 0 ? fci : 0;
 132                feasible = feasible && fci <= fc[i].tol;
 133                if (fci > 0) con2 += fci * fci;
 134           }
 135           *minf = fcur;
 136           minf_penalty = penalty;
 137           minf_feasible = feasible;
 138           d.rho = MAX(1e-6, MIN(10, 2 * fabs(*minf) / con2));
 139      }
 140      else
 141           d.rho = 1; /* whatever, doesn't matter */
 142
 143      if (auglag_verbose) {
 144           printf("auglag: initial rho=%g\nauglag initial lambda=", d.rho);
 145           for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]);
 146           printf("\nauglag initial mu = ");
 147           for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]);
 148           printf("\n");
 149      }
 150
 151      do {
 152           double prev_ICM = ICM;
 153
 154           ret = nlopt_optimize_limited(sub_opt, xcur, &fcur,
 155                                        stop->maxeval - stop->nevals,
 156                                        stop->maxtime - (nlopt_seconds()
 157                                                         - stop->start));
 158           if (ret < 0) break;
 159
 160           d.stop->nevals++;
 161           fcur = f(n, xcur, NULL, f_data);
 162
 163           ICM = 0;
 164           penalty = 0;
 165           feasible = 1;
 166           for (i = 0; i < d.p; ++i) {
 167                double hi = h[i].f(n, xcur, NULL, d.h[i].f_data);
 168                double newlam = d.lambda[i] + d.rho * hi;
 169                penalty += fabs(hi);
 170                feasible = feasible && fabs(hi) <= h[i].tol;
 171                ICM = MAX(ICM, fabs(hi));
 172                d.lambda[i] = MIN(MAX(lam_min, newlam), lam_max);
 173           }
 174           for (i = 0; i < d.m; ++i) {
 175                double fci = fc[i].f(n, xcur, NULL, d.fc[i].f_data);
 176                double newmu = d.mu[i] + d.rho * fci;
 177                penalty += fci > 0 ? fci : 0;
 178                feasible = feasible && fci <= fc[i].tol;
 179                ICM = MAX(ICM, fabs(MAX(fci, -d.mu[i] / d.rho)));
 180                d.mu[i] = MIN(MAX(0.0, newmu), mu_max);
 181           }
 182           if (ICM > tau * prev_ICM) {
 183                d.rho *= gam;
 184           }
 185
 186           auglag_iters++;
 187
 188           if (auglag_verbose) {
 189                printf("auglag %d: ICM=%g (%sfeasible), rho=%g\nauglag lambda=",
 190                       auglag_iters, ICM, feasible ? "" : "not ", d.rho);
 191                for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]);
 192                printf("\nauglag %d: mu = ", auglag_iters);
 193                for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]);
 194                printf("\n");
 195           }
 196
 197           if ((feasible && (!minf_feasible || penalty < minf_penalty
 198                             || fcur < *minf)) ||
 199               (!minf_feasible && penalty < minf_penalty)) {
 200                ret = NLOPT_SUCCESS;
 201                if (feasible) {
 202                     if (fcur < stop->minf_max)
 203                          ret = NLOPT_MINF_MAX_REACHED;
 204                     else if (nlopt_stop_ftol(stop, fcur, *minf))
 205                          ret = NLOPT_FTOL_REACHED;
 206                     else if (nlopt_stop_x(stop, xcur, x))
 207                          ret = NLOPT_XTOL_REACHED;
 208                }
 209                *minf = fcur;
 210                minf_penalty = penalty;
 211                minf_feasible = feasible;
 212                memcpy(x, xcur, sizeof(double) * n);
 213                if (ret != NLOPT_SUCCESS) break;
 214           }
 215
 216           if (nlopt_stop_forced(stop)) {ret = NLOPT_FORCED_STOP; break;}
 217           if (nlopt_stop_evals(stop)) {ret = NLOPT_MAXEVAL_REACHED; break;}
 218           if (nlopt_stop_time(stop)) {ret = NLOPT_MAXTIME_REACHED; break;}
 219
 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;
 228      } while (1);
 229
 230      free(xcur);
 231      return ret;
 232 }