From 72c427b1bb95992d3ddbcef0f4889a624d956764 Mon Sep 17 00:00:00 2001 From: stevenj Date: Mon, 18 Aug 2008 16:47:22 -0400 Subject: [PATCH] added nlopt_minimize_c function for constrained minimization (only via MMA for now); MMA now passes very simple test cases for constrained optimization darcs-hash:20080818204722-c8de0-3500ba9e001ac8d8f6a4593385c258ee6ed6ce4f.gz --- api/nlopt.c | 30 ++++++++++-- api/nlopt.h | 13 +++++ mma/mma.c | 135 +++++++++++++++++++++++++--------------------------- mma/mma.h | 4 +- 4 files changed, 106 insertions(+), 76 deletions(-) diff --git a/api/nlopt.c b/api/nlopt.c index b8019bf..e10448f 100644 --- a/api/nlopt.c +++ b/api/nlopt.c @@ -190,6 +190,7 @@ void nlopt_set_local_search_algorithm(nlopt_algorithm deriv, static nlopt_result nlopt_minimize_( nlopt_algorithm algorithm, int n, nlopt_func f, void *f_data, + int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, /* out: minimum */ @@ -210,6 +211,9 @@ static nlopt_result nlopt_minimize_( else if (n < 0 || !lb || !ub || !x) return NLOPT_INVALID_ARGS; + /* nonlinear constraints are only supported with MMA */ + if (m != 0 && algorithm != NLOPT_LD_MMA) return NLOPT_INVALID_ARGS; + d.f = f; d.f_data = f_data; d.lb = lb; @@ -409,7 +413,7 @@ static nlopt_result nlopt_minimize_( algorithm >= NLOPT_GN_MLSL_LDS); case NLOPT_LD_MMA: - return mma_minimize(n, f, f_data, 0, NULL, NULL, 0, + return mma_minimize(n, f, f_data, m, fc, fc_data, fc_datum_size, lb, ub, x, minf, &stop, local_search_alg_deriv, 1e-8, 100000); @@ -420,9 +424,10 @@ static nlopt_result nlopt_minimize_( return NLOPT_SUCCESS; } -nlopt_result nlopt_minimize( +nlopt_result nlopt_minimize_c( nlopt_algorithm algorithm, int n, nlopt_func f, void *f_data, + int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, /* out: minimum */ @@ -432,7 +437,8 @@ nlopt_result nlopt_minimize( { nlopt_result ret; if (xtol_abs) - ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub, + ret = nlopt_minimize_(algorithm, n, f, f_data, + m, fc, fc_data, fc_datum_size, lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, xtol_rel, xtol_abs, maxeval, maxtime); else { @@ -440,10 +446,26 @@ nlopt_result nlopt_minimize( double *xtol = (double *) malloc(sizeof(double) * n); if (!xtol) return NLOPT_OUT_OF_MEMORY; for (i = 0; i < n; ++i) xtol[i] = -1; - ret = nlopt_minimize_(algorithm, n, f, f_data, lb, ub, + ret = nlopt_minimize_(algorithm, n, f, f_data, + m, fc, fc_data, fc_datum_size, lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, xtol_rel, xtol, maxeval, maxtime); free(xtol); } return ret; } + +nlopt_result nlopt_minimize( + nlopt_algorithm algorithm, + int n, nlopt_func f, void *f_data, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, + int maxeval, double maxtime) +{ + return nlopt_minimize_c(algorithm, n, f, f_data, 0, NULL, NULL, 0, + lb, ub, x, minf, minf_max, ftol_rel, ftol_abs, + xtol_rel, xtol_abs, maxeval, maxtime); +} diff --git a/api/nlopt.h b/api/nlopt.h index edd04c4..fa7629a 100644 --- a/api/nlopt.h +++ b/api/nlopt.h @@ -1,6 +1,8 @@ #ifndef NLOPT_H #define NLOPT_H +#include /* for ptrdiff_t */ + #ifdef __cplusplus extern "C" { @@ -89,6 +91,17 @@ extern nlopt_result nlopt_minimize( double xtol_rel, const double *xtol_abs, int maxeval, double maxtime); +extern nlopt_result nlopt_minimize_c( + nlopt_algorithm algorithm, + int n, nlopt_func f, void *f_data, + int m, nlopt_func fc, void *fc_data, ptrdiff_t fc_datum_size, + const double *lb, const double *ub, /* bounds */ + double *x, /* in: initial guess, out: minimizer */ + double *minf, /* out: minimum */ + double minf_max, double ftol_rel, double ftol_abs, + double xtol_rel, const double *xtol_abs, + int maxeval, double maxtime); + extern void nlopt_srand(unsigned long seed); extern void nlopt_srand_time(void); diff --git a/mma/mma.c b/mma/mma.c index 4c595c4..e91a979 100644 --- a/mma/mma.c +++ b/mma/mma.c @@ -30,6 +30,8 @@ typedef struct { double gval, wval, *gcval; /* output each time (array length m) */ } dual_data; +static double sqr(double x) { return x * x; } + static double dual_func(int m, const double *y, double *grad, void *d_) { dual_data *d = (dual_data *) d_; @@ -47,79 +49,56 @@ static double dual_func(int m, const double *y, double *grad, void *d_) val = d->gval = fval; d->wval = 0; for (i = 0; i < m; ++i) val += y[i] * (gcval[i] = fcval[i]); - if (grad) { for (i = 0; i < m; ++i) grad[i] = fcval[i]; } for (j = 0; j < n; ++j) { - int x_pinned; - double a, b, dx, denom, denominv, ca, cb, sigma2; + double u, v, dx, denominv, c, sigma2, dx2; + + /* first, compute xcur[j] for y. Because this objective is + separable, we can minimize over x analytically, and the minimum + dx is given by the solution of a quadratic equation: + u dx^2 + 2 v sigma^2 dx + u sigma^2 = 0 + where u and v are defined by the sums below. Because of + the definitions, it is guaranteed that |u/v| <= sigma, + and it follows that the only dx solution with |dx| <= sigma + is given by: + (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2)) + (which goes to zero as u -> 0). */ - /* first, compute xcur[j] for y */ - a = dfdx[j]; - b = fabs(dfdx[j]) * sigma[j] + 0.5 * rho; + u = dfdx[j]; + v = fabs(dfdx[j]) * sigma[j] + 0.5 * rho; for (i = 0; i < m; ++i) { - a += dfcdx[i*n + j] * y[i]; - b += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i]; + u += dfcdx[i*n + j] * y[i]; + v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i]; } - sigma2 = sigma[j] * sigma[j]; - dx = -a * sigma2 / b; + u *= (sigma2 = sqr(sigma[j])); + dx = u==0 ? 0 : (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j])))); xcur[j] = x[j] + dx; if (xcur[j] > ub[j]) xcur[j] = ub[j]; - if (xcur[j] < lb[j]) xcur[j] = lb[j]; - if (xcur[j] > x[j] + 0.9*sigma[j]) xcur[j] = x[j] + 0.9*sigma[j]; - if (xcur[j] < x[j] - 0.9*sigma[j]) xcur[j] = x[j] - 0.9*sigma[j]; - x_pinned = xcur[j] != x[j] + dx; + else if (xcur[j] < lb[j]) xcur[j] = lb[j]; + if (xcur[j] > x[j]+0.9*sigma[j]) xcur[j] = x[j]+0.9*sigma[j]; + else if (xcur[j] < x[j]-0.9*sigma[j]) xcur[j] = x[j]-0.9*sigma[j]; dx = xcur[j] - x[j]; /* function value: */ - denom = sigma2 - dx*dx; - denominv = 1.0 / denom; - ca = sigma2 * dx; - cb = dx*dx; - val += (a*ca + b*cb) * denominv; + dx2 = dx * dx; + denominv = 1.0 / (sigma2 - dx2); + val += (u * dx + v * dx2) * denominv; - /* update gval, wval, gcval */ - d->gval += (dfdx[j] * ca + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * cb) + /* update gval, wval, gcval (approximant functions) */ + c = sigma2 * dx; + d->gval += (dfdx[j] * c + (fabs(dfdx[j])*sigma[j] + 0.5*rho) * dx2) * denominv; - d->wval += 0.5 * cb * denominv; + d->wval += 0.5 * dx2 * denominv; for (i = 0; i < m; ++i) - gcval[i] += (dfcdx[i*n+j] * ca - + (fabs(dfcdx[i*n+j])*sigma[j] + 0.5*rhoc[j]) * cb) + gcval[i] += (dfcdx[i*n+j] * c + (fabs(dfcdx[i*n+j])*sigma[j] + + 0.5*rhoc[j]) * dx2) * denominv; - - /* gradient */ - if (grad) - for (i = 0; i < m; ++i) { - double dady, dbdy; - dady = dfcdx[i*n + j]; - dbdy = fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]; - grad[i] += (dady*ca + dbdy*cb) * denominv; - if (!x_pinned) { /* dx/dy nonzero */ - int k; - double dxdy; - dxdy = (a*dbdy - dady*b) * (sigma2 / (b*b)); - grad[i] += - ((sigma2 * dfdx[j] + 2 * (fabs(dfdx[j])*sigma[j] - + 0.5*rho) * dx) - * denom - + 2*dx * (dfdx[j]*ca + (fabs(dfdx[j])*sigma[j] - + 0.5*rho) * cb)) - * (denominv*denominv * dxdy); - for (k = 0; k < m; ++k) - grad[i] += - ((sigma2 * dfcdx[k*n + j] - + 2 * (fabs(dfcdx[k*n + j])*sigma[j] - + 0.5*rhoc[k]) * dx) - * denom - + 2*dx * (dfcdx[k*n + j]*ca - + (fabs(dfcdx[k*n + j])*sigma[j] - + 0.5*rhoc[k]) * cb)) - * (denominv*denominv * dxdy) * y[k]; - } - } } - /* negate because we are maximizing */ - if (grad) { for (i = 0; i < m; ++i) grad[i] = -grad[i]; } + /* gradient is easy to compute: since we are at a minimum x (dval/dx=0), + we only need the partial derivative with respect to y, and + we negate because we are maximizing: */ + if (grad) for (i = 0; i < m; ++i) grad[i] = -gcval[i]; return -val; } @@ -127,7 +106,7 @@ static double dual_func(int m, const double *y, double *grad, void *d_) nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, int m, nlopt_func fc, - void *fc_data_, ptrdiff_t fc_data_size, + void *fc_data_, ptrdiff_t fc_datum_size, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, @@ -188,13 +167,13 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, dd.fval = fcur = *minf = f(n, x, dfdx, f_data); stop->nevals++; + memcpy(xcur, x, sizeof(double) * n); + feasible = 1; for (i = 0; i < m; ++i) { - fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_data_size * i); + fcval[i] = fc(n, x, dfcdx + i*n, fc_data + fc_datum_size * i); feasible = feasible && (fcval[i] <= 0); } - memcpy(xcur, x, sizeof(double) * n); - if (!feasible) { ret = NLOPT_FAILURE; goto done; } /* TODO: handle this */ while (1) { /* outer iterations */ @@ -209,13 +188,21 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, while (1) { /* inner iterations */ double min_dual; int feasible_cur, inner_done; + nlopt_result reti; /* solve dual problem */ dd.rho = rho; - nlopt_minimize(dual_alg, m, dual_func, &dd, - dual_lb, dual_ub, y, &min_dual, - -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval, - stop->maxtime - (nlopt_seconds() - stop->start)); + reti = nlopt_minimize( + dual_alg, m, dual_func, &dd, + dual_lb, dual_ub, y, &min_dual, + -HUGE_VAL, dual_tolrel,0., 0.,NULL, dual_maxeval, + stop->maxtime - (nlopt_seconds() - stop->start)); + if (reti == NLOPT_FAILURE) reti = NLOPT_SUCCESS; /* hack */ + if (reti < 0 || reti == NLOPT_MAXTIME_REACHED) { + ret = reti; + goto done; + } + dual_func(m, y, NULL, &dd); /* evaluate final xcur etc. */ fcur = f(n, xcur, dfdx_cur, f_data); @@ -224,7 +211,7 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, feasible_cur = 1; for (i = 0; i < m; ++i) { fcval_cur[i] = fc(n, xcur, dfcdx_cur + i*n, - fc_data + fc_data_size * i); + fc_data + fc_datum_size * i); feasible_cur = feasible_cur && (fcval_cur[i] <= 0); inner_done = inner_done && (dd.gcval[i] >= fcval_cur[i]); } @@ -255,6 +242,8 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, if (mma_verbose) printf("MMA inner iteration: rho -> %g\n", rho); + for (i = 0; i < mma_verbose; ++i) + printf(" rhoc[%d] -> %g\n", i,rhoc[i]); } if (nlopt_stop_ftol(stop, fcur, fprev)) @@ -269,14 +258,22 @@ nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, printf("MMA outer iteration: rho -> %g\n", rho); for (i = 0; i < m; ++i) rhoc[i] = MAX(0.1 * rhoc[i], MMA_RHOMIN); - if (k > 1) + for (i = 0; i < mma_verbose; ++i) + printf(" rhoc[%d] -> %g\n", i, rhoc[i]); + if (k > 1) { for (j = 0; j < n; ++j) { double dx2 = (xcur[j]-xprev[j]) * (xprev[j]-xprevprev[j]); double gam = dx2 < 0 ? 0.7 : (dx2 > 0 ? 1.2 : 1); sigma[j] *= gam; - sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j])); - sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j])); + if (!nlopt_isinf(ub[j]) && !nlopt_isinf(lb[j])) { + sigma[j] = MIN(sigma[j], 10*(ub[j]-lb[j])); + sigma[j] = MAX(sigma[j], 0.01*(ub[j]-lb[j])); + } } + for (j = 0; j < mma_verbose; ++j) + printf(" sigma[%d] -> %g\n", + j, sigma[j]); + } } done: diff --git a/mma/mma.h b/mma/mma.h index 929a228..f6aac37 100644 --- a/mma/mma.h +++ b/mma/mma.h @@ -4,8 +4,6 @@ #include "nlopt.h" #include "nlopt-util.h" -#include /* for ptrdiff_t */ - #ifdef __cplusplus extern "C" { @@ -15,7 +13,7 @@ extern int mma_verbose; nlopt_result mma_minimize(int n, nlopt_func f, void *f_data, int m, nlopt_func fc, - void *fc_data_, ptrdiff_t fc_data_size, + void *fc_data_, ptrdiff_t fc_datum_size, const double *lb, const double *ub, /* bounds */ double *x, /* in: initial guess, out: minimizer */ double *minf, -- 2.30.2