chiark / gitweb /
initial (undocumented) support for equality constraints via augmented Lagrangian...
authorstevenj <stevenj@alum.mit.edu>
Wed, 26 Nov 2008 04:23:30 +0000 (23:23 -0500)
committerstevenj <stevenj@alum.mit.edu>
Wed, 26 Nov 2008 04:23:30 +0000 (23:23 -0500)
darcs-hash:20081126042330-c8de0-c6074b36c1bf1a41005c1b1526ce801659db703a.gz

Makefile.am
api/Makefile.am
api/nlopt.c
api/nlopt.h
auglag/Makefile.am [new file with mode: 0644]
auglag/README [new file with mode: 0644]
auglag/auglag.c [new file with mode: 0644]
auglag/auglag.h [new file with mode: 0644]
configure.ac
util/nlopt-util.h

index 9ddb4709acb6dee856789ac404af24b5ded0145b..ca25601b06ebc97ad8596a01b280103e3b2037d2 100644 (file)
@@ -8,7 +8,7 @@ CXX_DIRS = stogo
 CXX_LIBS = stogo/libstogo.la
 endif
 
-SUBDIRS = util direct cdirect $(CXX_DIRS) praxis luksan crs mlsl mma cobyla newuoa lbfgs neldermead api . octave test
+SUBDIRS = util direct cdirect $(CXX_DIRS) praxis luksan crs mlsl mma cobyla newuoa lbfgs neldermead auglag api . octave test
 EXTRA_DIST = autogen.sh nlopt.pc.in m4
 
 if WITH_NOCEDAL
@@ -19,7 +19,7 @@ libnlopt@NLOPT_SUFFIX@_la_SOURCES =
 libnlopt@NLOPT_SUFFIX@_la_LIBADD = \
 direct/libdirect.la cdirect/libcdirect.la $(CXX_LIBS)                  \
 praxis/libpraxis.la $(NOCEDAL_LBFGS) luksan/libluksan.la crs/libcrs.la \
-mlsl/libmlsl.la mma/libmma.la cobyla/libcobyla.la newuoa/libnewuoa.la neldermead/libneldermead.la api/libapi.la util/libutil.la
+mlsl/libmlsl.la mma/libmma.la cobyla/libcobyla.la newuoa/libnewuoa.la neldermead/libneldermead.la auglag/libauglag.la api/libapi.la util/libutil.la
 
 libnlopt@NLOPT_SUFFIX@_la_LDFLAGS = -no-undefined -version-info @SHARED_VERSION_INFO@
 
index 2e1767883123c0a810f24cbba9d8e08b82c28ebb..51356ae36c5aa09001b86093f59ebf1968eedc32 100644 (file)
@@ -1,4 +1,4 @@
-AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/praxis -I$(top_srcdir)/lbfgs -I$(top_srcdir)/luksan -I$(top_srcdir)/crs -I$(top_srcdir)/mlsl -I$(top_srcdir)/mma -I$(top_srcdir)/cobyla -I$(top_srcdir)/newuoa -I$(top_srcdir)/neldermead -I$(top_srcdir)/util
+AM_CPPFLAGS = -I$(top_srcdir)/cdirect -I$(top_srcdir)/direct -I$(top_srcdir)/stogo -I$(top_srcdir)/praxis -I$(top_srcdir)/lbfgs -I$(top_srcdir)/luksan -I$(top_srcdir)/crs -I$(top_srcdir)/mlsl -I$(top_srcdir)/mma -I$(top_srcdir)/cobyla -I$(top_srcdir)/newuoa -I$(top_srcdir)/neldermead -I$(top_srcdir)/auglag -I$(top_srcdir)/util
 
 include_HEADERS = nlopt.h nlopt.f
 noinst_LTLIBRARIES = libapi.la
index 28fd621da03a44f7393fe8d6140715f644ae937e..9ef5c3bb0f5e631f0b1ae91c8d0b04c900f09d3c 100644 (file)
@@ -95,7 +95,11 @@ static const char nlopt_algorithm_names[NLOPT_NUM_ALGORITHMS][256] = {
      "NEWUOA unconstrained optimization via quadratic models (local, no-derivative)",
      "Bound-constrained optimization via NEWUOA-based quadratic models (local, no-derivative)",
      "Nelder-Mead simplex algorithm (local, no-derivative)",
-     "Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)"
+     "Sbplx variant of Nelder-Mead (re-implementation of Rowan's Subplex) (local, no-derivative)",
+     "Augmented Lagrangian method (local, no-derivative)",
+     "Augmented Lagrangian method (local, derivative)",
+     "Augmented Lagrangian method for equality constraints (local, no-derivative)",
+     "Augmented Lagrangian method for equality constraints (local, derivative)",
 };
 
 const char *nlopt_algorithm_name(nlopt_algorithm a)
@@ -217,6 +221,12 @@ static double f_direct(int n, const double *x, int *undefined, void *data_)
 #include "cobyla.h"
 #include "newuoa.h"
 #include "neldermead.h"
+#include "auglag.h"
+
+#define AUGLAG_ALG(a) ((a) == NLOPT_LN_AUGLAG ||       \
+                      (a) == NLOPT_LN_AUGLAG_EQ ||     \
+                      (a) == NLOPT_LD_AUGLAG ||        \
+                      (a) == NLOPT_LD_AUGLAG_EQ)
 
 /*************************************************************************/
 
@@ -264,6 +274,7 @@ static nlopt_result nlopt_minimize_(
      double *minf, /* out: minimum */
      double minf_max, double ftol_rel, double ftol_abs,
      double xtol_rel, const double *xtol_abs,
+     double htol_rel, double htol_abs,
      int maxeval, double maxtime)
 {
      int i;
@@ -280,12 +291,13 @@ static nlopt_result nlopt_minimize_(
      else if (n < 0 || !lb || !ub || !x)
          return NLOPT_INVALID_ARGS;
 
-     /* nonlinear constraints are only supported with MMA or COBYLA */
-     if (m != 0 && algorithm != NLOPT_LD_MMA && algorithm != NLOPT_LN_COBYLA) 
+     /* nonlinear constraints are only supported with some algorithms */
+     if (m != 0 && algorithm != NLOPT_LD_MMA && algorithm != NLOPT_LN_COBYLA
+         && !AUGLAG_ALG(algorithm)) 
          return NLOPT_INVALID_ARGS;
 
-     /* nonlinear equality constraints (h(x) = 0) not yet supported */
-     if (p != 0)
+     /* nonlinear equality constraints (h(x) = 0) only via some algorithms */
+     if (p != 0 && !AUGLAG_ALG(algorithm))
          return NLOPT_INVALID_ARGS;
 
      d.f = f;
@@ -310,6 +322,8 @@ static nlopt_result nlopt_minimize_(
      stop.ftol_abs = ftol_abs;
      stop.xtol_rel = xtol_rel;
      stop.xtol_abs = xtol_abs;
+     stop.htol_rel = htol_rel;
+     stop.htol_abs = htol_abs;
      stop.nevals = 0;
      stop.maxeval = maxeval;
      stop.maxtime = maxtime;
@@ -514,6 +528,24 @@ static nlopt_result nlopt_minimize_(
              return ret;
         }
 
+        case NLOPT_LN_AUGLAG:
+        case NLOPT_LN_AUGLAG_EQ:
+             return auglag_minimize(n, f, f_data, 
+                                    m, fc, fc_data, fc_datum_size,
+                                    p, h, h_data, h_datum_size,
+                                    lb, ub, x, minf, &stop,
+                                    local_search_alg_nonderiv,
+                                    algorithm == NLOPT_LN_AUGLAG_EQ);
+
+        case NLOPT_LD_AUGLAG:
+        case NLOPT_LD_AUGLAG_EQ:
+             return auglag_minimize(n, f, f_data, 
+                                    m, fc, fc_data, fc_datum_size,
+                                    p, h, h_data, h_datum_size,
+                                    lb, ub, x, minf, &stop,
+                                    local_search_alg_deriv,
+                                    algorithm == NLOPT_LD_AUGLAG_EQ);
+
         default:
              return NLOPT_INVALID_ARGS;
      }
@@ -531,6 +563,7 @@ nlopt_result nlopt_minimize_econstrained(
      double *minf, /* out: minimum */
      double minf_max, double ftol_rel, double ftol_abs,
      double xtol_rel, const double *xtol_abs,
+     double htol_rel, double htol_abs,
      int maxeval, double maxtime)
 {
      nlopt_result ret;
@@ -540,7 +573,8 @@ nlopt_result nlopt_minimize_econstrained(
                                p, h, h_data, h_datum_size,
                                lb, ub,
                                x, minf, minf_max, ftol_rel, ftol_abs,
-                               xtol_rel, xtol_abs, maxeval, maxtime);
+                               xtol_rel, xtol_abs, htol_rel, htol_abs,
+                               maxeval, maxtime);
      else {
          int i;
          double *xtol = (double *) malloc(sizeof(double) * n);
@@ -551,7 +585,8 @@ nlopt_result nlopt_minimize_econstrained(
                                p, h, h_data, h_datum_size,
                                lb, ub,
                                x, minf, minf_max, ftol_rel, ftol_abs,
-                               xtol_rel, xtol, maxeval, maxtime);
+                               xtol_rel, xtol, htol_rel, htol_abs,
+                               maxeval, maxtime);
          free(xtol);
      }
      return ret;
@@ -572,7 +607,7 @@ nlopt_result nlopt_minimize_constrained(
          algorithm, n, f, f_data, 
          m, fc, fc_data, fc_datum_size, 0, NULL, NULL, 0,
          lb, ub, x, minf, minf_max, ftol_rel, ftol_abs,
-         xtol_rel, xtol_abs, maxeval, maxtime);
+         xtol_rel, xtol_abs, ftol_rel, ftol_abs, maxeval, maxtime);
 }
 
 nlopt_result nlopt_minimize(
index e2362be2dad202f7b9c7e1433df7fd44e952dd93..2b2b0423310f2fd6bf210ed8851e1b78da735aeb 100644 (file)
@@ -91,6 +91,11 @@ typedef enum {
      NLOPT_LN_NELDERMEAD,
      NLOPT_LN_SBPLX,
 
+     NLOPT_LN_AUGLAG,
+     NLOPT_LD_AUGLAG,
+     NLOPT_LN_AUGLAG_EQ,
+     NLOPT_LD_AUGLAG_EQ,
+
      NLOPT_NUM_ALGORITHMS /* not an algorithm, just the number of them */
 } nlopt_algorithm;
 
@@ -140,6 +145,7 @@ extern nlopt_result nlopt_minimize_econstrained(
      double *minf, /* out: minimum */
      double minf_max, double ftol_rel, double ftol_abs,
      double xtol_rel, const double *xtol_abs,
+     double htol_rel, double htol_abs,
      int maxeval, double maxtime);
 
 extern void nlopt_srand(unsigned long seed);
diff --git a/auglag/Makefile.am b/auglag/Makefile.am
new file mode 100644 (file)
index 0000000..7795415
--- /dev/null
@@ -0,0 +1,6 @@
+AM_CPPFLAGS = -I$(top_srcdir)/util -I$(top_srcdir)/api
+
+noinst_LTLIBRARIES = libauglag.la
+libauglag_la_SOURCES = auglag.c auglag.h
+
+EXTRA_DIST = README
diff --git a/auglag/README b/auglag/README
new file mode 100644 (file)
index 0000000..e1f8c2f
--- /dev/null
@@ -0,0 +1,25 @@
+This directory contains my implementation of the "augmented Lagrangian"
+method to express constrained optimization problems (including equality
+constraints) in terms of unconstrained optimization problems (or just
+inequality constraints, or just box constraints).
+
+I used the algorithm description (but no source code) from the outline
+in the paper:
+
+       E. G. Birgin and J. M. Martinez, "Improving ultimate convergence
+       of an augmented Lagrangian method," Optimization Methods and
+       Software vol. 23, no. 2, p. 177-195 (2008). 
+
+       http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.72.6121
+
+(The authors of this paper have their own free-software implementation
+of this idea and its variations, in the TANGO project:
+http://www.ime.usp.br/~egbirgin/tango/ ....I did *not* use any code
+from TANGO here, or even look at the TANGO code.  TANGO is GPLed, and
+I'm trying to keep NLopt at most LGPLed.)
+
+The code in this directory is under the same MIT license as the rest
+of my code in NLopt (see ../COPYRIGHT).
+
+Steven G. Johnson
+November 2008
diff --git a/auglag/auglag.c b/auglag/auglag.c
new file mode 100644 (file)
index 0000000..ad47f98
--- /dev/null
@@ -0,0 +1,169 @@
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "auglag.h"
+
+int auglag_verbose = 1;
+
+#define MIN(a,b) ((a) < (b) ? (a) : (b))
+#define MAX(a,b) ((a) > (b) ? (a) : (b))
+
+/***************************************************************************/
+
+typedef struct {
+     nlopt_func f; void *f_data;
+     int m; nlopt_func fc; char *fc_data; ptrdiff_t fc_datum_size;
+     int p; nlopt_func h; char *h_data; ptrdiff_t h_datum_size;
+     double rho, *lambda, *mu;
+     double *gradtmp;
+     nlopt_stopping *stop;
+} auglag_data;
+
+/* the augmented lagrangian objective function */
+static double auglag(int n, const double *x, double *grad, void *data)
+{
+     auglag_data *d = (auglag_data *) data;
+     double *gradtmp = grad ? d->gradtmp : NULL;
+     double rho = d->rho;
+     const double *lambda = d->lambda, *mu = d->mu;
+     double L;
+     int i, j;
+
+     L = d->f(n, x, grad, d->f_data);
+
+     for (i = 0; i < d->p; ++i) {
+         double h;
+         h = d->h(n, x, gradtmp, d->h_data + d->h_datum_size * i)
+              + lambda[i] / rho;
+         L += 0.5 * rho * h*h;
+         if (grad) for (j = 0; j < n; ++j) grad[j] += (rho * h) * gradtmp[j];
+     }
+
+     for (i = 0; i < d->m; ++i) {
+         double fc;
+         fc = d->fc(n, x, gradtmp, d->fc_data + d->fc_datum_size * i)
+              + mu[i] / rho;
+         if (fc > 0) {
+              L += 0.5 * rho * fc*fc;
+              if (grad) for (j = 0; j < n; ++j) 
+                   grad[j] += (rho * fc) * gradtmp[j];
+         }
+     }
+     
+     d->stop->nevals++;
+
+     return L;
+}
+
+/***************************************************************************/
+
+nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
+                            int m, nlopt_func fc,
+                            void *fc_data, ptrdiff_t fc_datum_size,
+                            int p, nlopt_func h,
+                            void *h_data, ptrdiff_t h_datum_size,
+                            const double *lb, const double *ub, /* bounds */
+                            double *x, /* in: initial guess, out: minimizer */
+                            double *minf,
+                            nlopt_stopping *stop,
+                            nlopt_algorithm sub_alg, int sub_has_fc)
+{
+     auglag_data d;
+     nlopt_result ret = NLOPT_SUCCESS;
+     double ICM = HUGE_VAL;
+     int i;
+
+     /* magic parameters from Birgin & Martinez */
+     const double tau = 0.5, gam = 10;
+     const double lam_min = -1e20, lam_max = 1e20, mu_max = 1e20;
+
+     d.f = f; d.f_data = f_data;
+     d.m = m; d.fc = fc; d.fc_data = (char *) fc_data;
+     d.fc_datum_size = fc_datum_size;
+     d.p = p; d.h = h; d.h_data = (char *) h_data;
+     d.h_datum_size = h_datum_size;
+     d.stop = stop;
+
+     /* whether we handle inequality constraints via the augmented
+       Lagrangian penalty function, or directly in the sub-algorithm */
+     if (sub_has_fc)
+         d.m = 0;
+     else
+         m = 0;
+
+     d.gradtmp = (double *) malloc(sizeof(double) * (n + d.p + d.m));
+     if (!d.gradtmp) return NLOPT_OUT_OF_MEMORY;
+     memset(d.gradtmp, 0, sizeof(double) * (n + d.p + d.m));
+     d.lambda = d.gradtmp + n;
+     d.mu = d.lambda + d.p;
+
+     /* starting rho suggested by B & M */
+     {
+         double con2 = 0;
+         *minf = f(n, x, NULL, f_data);
+         for (i = 0; i < d.p; ++i) {
+               double hi = h(n, x, NULL, d.h_data + i*h_datum_size);
+              con2 += hi * hi;
+         }
+         for (i = 0; i < d.m; ++i) {
+               double fci = fc(n, x, NULL, d.fc_data + i*fc_datum_size);
+              if (fci > 0) con2 += fci * fci;
+         }
+         d.rho = MAX(1e-6, MIN(10, 2 * fabs(*minf) / con2));
+     }
+
+     *minf = HUGE_VAL;
+
+     do {
+         double prev_ICM = ICM;
+
+         ret = nlopt_minimize_constrained(sub_alg, n, auglag, &d,
+                                          m, fc, fc_data, fc_datum_size,
+                                          lb, ub, x, minf,
+                                          -HUGE_VAL, 
+                                          stop->ftol_rel, stop->ftol_abs,
+                                          stop->xtol_rel, stop->xtol_abs,
+                                          stop->maxeval - stop->nevals,
+                                          stop->maxtime - (nlopt_seconds() 
+                                                           - stop->start));
+         if (ret < 0) break;
+         if (nlopt_stop_evals(stop)) {ret = NLOPT_MAXEVAL_REACHED; break;}
+          if (nlopt_stop_time(stop)) {ret = NLOPT_MAXTIME_REACHED; break;}
+         
+         *minf = f(n, x, NULL, f_data);
+         
+         ICM = 0;
+         for (i = 0; i < d.p; ++i) {
+              double hi = h(n, x, NULL, d.h_data + i*h_datum_size);
+              double newlam = d.lambda[i] + d.rho * hi;
+              ICM = MAX(ICM, fabs(hi));
+              d.lambda[i] = MIN(MAX(lam_min, newlam), lam_max);
+         }
+         for (i = 0; i < d.m; ++i) {
+              double fci = fc(n, x, NULL, d.fc_data + i*fc_datum_size);
+              double newmu = d.mu[i] + d.rho * fci;
+              ICM = MAX(ICM, fabs(MAX(fci, -d.mu[i] / d.rho)));
+              d.mu[i] = MIN(MAX(0.0, newmu), mu_max);
+         }
+         if (ICM > tau * prev_ICM)
+              d.rho *= gam;
+
+         if (auglag_verbose) {
+              printf("auglag: ICM=%g, rho=%g\nauglag lambda=", ICM, d.rho);
+              for (i = 0; i < d.p; ++i) printf(" %g", d.lambda[i]);
+              printf("\nauglag mu = ");
+              for (i = 0; i < d.m; ++i) printf(" %g", d.mu[i]);
+              printf("\n");
+         }
+
+         if (ICM <= stop->htol_abs || ICM <= stop->htol_rel * fabs(*minf)) {
+              ret = NLOPT_FTOL_REACHED;
+              break;
+         }
+     } while (1);
+
+     free(d.gradtmp);
+     return ret;
+}
diff --git a/auglag/auglag.h b/auglag/auglag.h
new file mode 100644 (file)
index 0000000..4bbbcfa
--- /dev/null
@@ -0,0 +1,52 @@
+/* Copyright (c) 2007-2008 Massachusetts Institute of Technology
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining
+ * a copy of this software and associated documentation files (the
+ * "Software"), to deal in the Software without restriction, including
+ * without limitation the rights to use, copy, modify, merge, publish,
+ * distribute, sublicense, and/or sell copies of the Software, and to
+ * permit persons to whom the Software is furnished to do so, subject to
+ * the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be
+ * included in all copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+ * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ * MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+ * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
+ * LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
+ * OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
+ * WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
+ */
+
+#ifndef AUGLAG_H
+#define AUGLAG_H
+
+#include "nlopt.h"
+#include "nlopt-util.h"
+
+#ifdef __cplusplus
+extern "C"
+{
+#endif /* __cplusplus */
+
+extern int auglag_verbose;
+
+nlopt_result auglag_minimize(int n, nlopt_func f, void *f_data,
+                            int m, nlopt_func fc,
+                            void *fc_data, ptrdiff_t fc_datum_size,
+                            int p, nlopt_func h,
+                            void *h_data, ptrdiff_t h_datum_size,
+                            const double *lb, const double *ub, /* bounds */
+                            double *x, /* in: initial guess, out: minimizer */
+                            double *minf,
+                            nlopt_stopping *stop,
+                            nlopt_algorithm sub_alg, int sub_has_fc);
+
+#ifdef __cplusplus
+}  /* extern "C" */
+#endif /* __cplusplus */
+
+#endif
+
index 992732d154d698705e59b47e38e7a620d5f003b2..93e7094375d300fc96fc651ea8546c5a0ee7e089 100644 (file)
@@ -236,6 +236,7 @@ AC_CONFIG_FILES([
    cobyla/Makefile
    newuoa/Makefile
    neldermead/Makefile
+   auglag/Makefile
    test/Makefile
 ])
 
index 586f35e0c19bc9c1b28df0a9082b12007730bc06..7a1875b3695e125887914de85a72dc93364d82ea 100644 (file)
@@ -70,6 +70,7 @@ typedef struct {
      double ftol_abs;
      double xtol_rel;
      const double *xtol_abs;
+     double htol_rel, htol_abs;
      int nevals, maxeval;
      double maxtime, start;
 } nlopt_stopping;