From 30dcb46c2cc2d36795d89cb07f54f488f2223d4f Mon Sep 17 00:00:00 2001 From: stevenj Date: Tue, 25 Nov 2008 23:23:30 -0500 Subject: [PATCH] initial (undocumented) support for equality constraints via augmented Lagrangian method darcs-hash:20081126042330-c8de0-c6074b36c1bf1a41005c1b1526ce801659db703a.gz --- Makefile.am | 4 +- api/Makefile.am | 2 +- api/nlopt.c | 51 +++++++++++--- api/nlopt.h | 6 ++ auglag/Makefile.am | 6 ++ auglag/README | 25 +++++++ auglag/auglag.c | 169 +++++++++++++++++++++++++++++++++++++++++++++ auglag/auglag.h | 52 ++++++++++++++ configure.ac | 1 + util/nlopt-util.h | 1 + 10 files changed, 306 insertions(+), 11 deletions(-) create mode 100644 auglag/Makefile.am create mode 100644 auglag/README create mode 100644 auglag/auglag.c create mode 100644 auglag/auglag.h diff --git a/Makefile.am b/Makefile.am index 9ddb470..ca25601 100644 --- a/Makefile.am +++ b/Makefile.am @@ -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@ diff --git a/api/Makefile.am b/api/Makefile.am index 2e17678..51356ae 100644 --- a/api/Makefile.am +++ b/api/Makefile.am @@ -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 diff --git a/api/nlopt.c b/api/nlopt.c index 28fd621..9ef5c3b 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -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( diff --git a/api/nlopt.h b/api/nlopt.h index e2362be..2b2b042 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -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 index 0000000..7795415 --- /dev/null +++ b/auglag/Makefile.am @@ -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 index 0000000..e1f8c2f --- /dev/null +++ b/auglag/README @@ -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 index 0000000..ad47f98 --- /dev/null +++ b/auglag/auglag.c @@ -0,0 +1,169 @@ +#include +#include +#include +#include + +#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 index 0000000..4bbbcfa --- /dev/null +++ b/auglag/auglag.h @@ -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 + diff --git a/configure.ac b/configure.ac index 992732d..93e7094 100644 --- a/configure.ac +++ b/configure.ac @@ -236,6 +236,7 @@ AC_CONFIG_FILES([ cobyla/Makefile newuoa/Makefile neldermead/Makefile + auglag/Makefile test/Makefile ]) diff --git a/util/nlopt-util.h b/util/nlopt-util.h index 586f35e..7a1875b 100644 --- a/util/nlopt-util.h +++ b/util/nlopt-util.h @@ -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; -- 2.30.2