From 09d57a15286e5328e60c209a42cc7d79f6f90533 Mon Sep 17 00:00:00 2001 From: Julien Schueller Date: Mon, 11 Nov 2024 17:15:20 +0100 Subject: [PATCH] optimize: transform box constraints Closes #528 --- NEWS.md | 5 +++ doc/docs/NLopt_Algorithms.md | 2 +- src/api/optimize.c | 75 +++++++++++++++++++++++++++++++----- test/CMakeLists.txt | 3 +- test/t_fbound.cxx | 46 ++++++++++++++++++++++ 5 files changed, 119 insertions(+), 12 deletions(-) create mode 100644 test/t_fbound.cxx diff --git a/NEWS.md b/NEWS.md index a51fd810..78efaa38 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,9 @@ # NLopt Release Notes +in-progress + +* Fixed PRAXIS box constraints ([#528]) + ## NLopt 2.9 10 November 2024 @@ -517,6 +521,7 @@ [#504]: https://github.com/stevengj/nlopt/issues/504 [#509]: https://github.com/stevengj/nlopt/issues/509 [#510]: https://github.com/stevengj/nlopt/issues/510 +[#528]: https://github.com/stevengj/nlopt/issues/528 [#533]: https://github.com/stevengj/nlopt/issues/533 [#534]: https://github.com/stevengj/nlopt/issues/534 [#535]: https://github.com/stevengj/nlopt/issues/535 diff --git a/doc/docs/NLopt_Algorithms.md b/doc/docs/NLopt_Algorithms.md index 40b2032b..fc8eb246 100644 --- a/doc/docs/NLopt_Algorithms.md +++ b/doc/docs/NLopt_Algorithms.md @@ -258,7 +258,7 @@ The original Fortran code was written by Richard Brent and made available by the Specified in NLopt as `NLOPT_LN_PRAXIS` -This algorithm was originally designed for unconstrained optimization. In NLopt, bound constraints are "implemented" in PRAXIS by the simple expedient of returning infinity (Inf) when the constraints are violated (this is done automatically—you don't have to do this in your own function). This seems to work, more-or-less, but appears to slow convergence significantly. If you have bound constraints, you are probably better off using COBYLA or BOBYQA. +This algorithm was originally designed for unconstrained optimization. In NLopt, bound constraints are emulated by variable transformation into an unbound space. ### Nelder-Mead Simplex diff --git a/src/api/optimize.c b/src/api/optimize.c index 5f542b5a..09af1909 100644 --- a/src/api/optimize.c +++ b/src/api/optimize.c @@ -57,23 +57,68 @@ #include "slsqp.h" /*********************************************************************/ +/* to emulate box constraint support with variable transform */ -static double f_bound(int n, const double *x, void *data_) +static void x_bound_inv(int n, double *x, const double *lb, const double *ub) +{ + int i; + double mid, width, th; + if (!lb || !ub) + return; + for (i = 0; i < n; ++ i) + { + if (!nlopt_isinf(lb[i]) && !nlopt_isinf(ub[i])) + { + mid = (lb[i] + ub[i]) * 0.5; + width = (ub[i] - lb[i]) * 0.5; + th = tanh(x[i]); + x[i] = mid + th * width; + } + else if (!nlopt_isinf(lb[i])) + x[i] = sqrt(x[i] - lb[i]); + else if (!nlopt_isinf(ub[i])) + x[i] = sqrt(ub[i] - x[i]); + } +} + +static void x_bound(int n, double *x, const double *lb, const double *ub) { int i; + double mid, width, th; + if (!lb || !ub) + return; + for (i = 0; i < n; ++ i) + { + if (!nlopt_isinf(lb[i]) && !nlopt_isinf(ub[i])) + { + mid = (lb[i] + ub[i]) * 0.5; + width = (ub[i] - lb[i]) * 0.5; + th = tanh(x[i]); + x[i] = mid + th * width; + } + else if (!nlopt_isinf(lb[i])) + x[i] = lb[i] + x[i] * x[i]; + else if (!nlopt_isinf(ub[i])) + x[i] = ub[i] - x[i] * x[i]; + } +} + +static double f_bound(int n, const double *x, void *data_) +{ nlopt_opt data = (nlopt_opt) data_; double f; + double *x_tr = NULL; - /* some methods do not support bound constraints, but support - discontinuous objectives so we can just return Inf for invalid x */ - for (i = 0; i < n; ++i) - if (x[i] < data->lb[i] || x[i] > data->ub[i]) - return HUGE_VAL; - - f = data->f((unsigned) n, x, NULL, data->f_data); - return (nlopt_isnan(f) || nlopt_isinf(f) ? HUGE_VAL : f); + x_tr = (double *) malloc(n * sizeof(double)); + memcpy(x_tr, x, n * sizeof(double)); + x_bound(n, x_tr, data->lb, data->ub); + f = data->f((unsigned) n, x_tr, NULL, data->f_data); + free(x_tr); + return f; } +/*********************************************************************/ + static double f_noderiv(int n, const double *x, void *data_) { nlopt_opt data = (nlopt_opt) data_; @@ -649,9 +694,19 @@ static nlopt_result nlopt_optimize_(nlopt_opt opt, double *x, double *minf) case NLOPT_LN_PRAXIS: { double step; + nlopt_result result; + if (initial_step(opt, x, &step) != NLOPT_SUCCESS) return NLOPT_OUT_OF_MEMORY; - return praxis_(nlopt_get_param(opt, "t0_tol", 0.0), DBL_EPSILON, step, ni, x, f_bound, opt, &stop, minf); + + /* transform starting point into unbound space */ + x_bound_inv(ni, x, opt->lb, opt->ub); + + result = praxis_(nlopt_get_param(opt, "t0_tol", 0.0), DBL_EPSILON, step, ni, x, f_bound, opt, &stop, minf); + + /* transform back optimal point to bounded space */ + x_bound(ni, x, opt->lb, opt->ub); + return result; } case NLOPT_LD_LBFGS: diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e2fcb027..cd367e9f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -15,7 +15,8 @@ macro(NLOPT_add_cpp_test test_name) endmacro() NLOPT_add_cpp_test(t_tutorial 23 24 30 39) -NLOPT_add_cpp_test(cpp_functor) +NLOPT_add_cpp_test(cpp_functor 0) +NLOPT_add_cpp_test(t_fbound 0) NLOPT_add_cpp_test(t_bounded 0 1 2 3 4 5 6 7 8 18 34 41 42) if (NOT NLOPT_CXX) diff --git a/test/t_fbound.cxx b/test/t_fbound.cxx new file mode 100644 index 00000000..b2fb583e --- /dev/null +++ b/test/t_fbound.cxx @@ -0,0 +1,46 @@ +#include +#include +#include +#include +#include +#include + +double myvfunc(const std::vector &x, std::vector &grad, void *data) +{ + (void)data; + if (!grad.empty()) { + grad[0] = 2*x[0]; + grad[1] = 2*x[1]; + } + if (std::isnan(x[0])) + throw std::invalid_argument("nan"); + return -x[0]*x[0]-x[1]*x[1]; +} + + +int main() { + nlopt::srand(0); + for (int repeat = 0; repeat < 100; ++repeat) + { + std::cout << "repeat="< lb = {-1.0, -10.0}; + std::vector ub = {10.0, 1.0}; + opt.set_lower_bounds(lb); + opt.set_upper_bounds(ub); + opt.set_min_objective(myvfunc, NULL); + opt.set_ftol_rel(1e-6); + std::vector x = {0.5, 0.5}; + double minf; + try{ + opt.optimize(x, minf); + std::cout << "found optimum at f(" << x[0] << "," << x[1] << ") = " + << std::setprecision(10) << minf <