| Title: | R Interface to NLopt |
|---|---|
| Description: | Solve optimization problems using an R interface to NLopt. NLopt is a free/open-source library for nonlinear optimization, providing a common interface for a number of different free optimization routines available online as well as original implementations of various other algorithms. See <https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/> for more information on the available algorithms. Building from included sources requires 'CMake'. On Linux and 'macOS', if a suitable system build of NLopt (2.7.0 or later) is found, it is used; otherwise, it is built from included sources via 'CMake'. On Windows, NLopt is obtained through 'rwinlib' for 'R <= 4.1.x' or grabbed from the appropriate toolchain for 'R >= 4.2.0'. |
| Authors: | Jelmer Ypma [aut], Steven G. Johnson [aut] (author of the NLopt
C library), Aymeric Stamm [ctb, cre]
 , Hans W. Borchers
[ctb], Dirk Eddelbuettel [ctb], Brian Ripley [ctb] (build
process on multiple OS), Kurt Hornik [ctb] (build process on
multiple OS), Julien Chiquet [ctb], Avraham Adler [ctb]
, Xiongtao Dai [ctb],
Jeroen Ooms [ctb] |
| Maintainer: | Aymeric Stamm <[email protected]> |
| License: | LGPL (>= 3) |
| Version: | 2.1.1 |
| Built: | 2024-11-03 06:40:17 UTC |
| Source: | CRAN |
The Augmented Lagrangian method adds additional terms to the unconstrained objective function, designed to emulate a Lagrangian multiplier.
auglag( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, heq = NULL, heqjac = NULL, localsolver = "COBYLA", localtol = 1e-06, ineq2local = FALSE, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )auglag( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, heq = NULL, heqjac = NULL, localsolver = "COBYLA", localtol = 1e-06, ineq2local = FALSE, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
gradient of the objective function; will be provided provided is
|
lower, upper
|
lower and upper bound constraints. |
hin, hinjac
|
defines the inequality constraints, |
heq, heqjac
|
defines the equality constraints, |
localsolver |
available local solvers: COBYLA, LBFGS, MMA, or SLSQP. |
localtol |
tolerance applied in the selected local solver. |
ineq2local |
logical; shall the inequality constraints be treated by the local solver?; not possible at the moment. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
deprecatedBehavior |
logical; if |
... |
additional arguments passed to the function. |
This method combines the objective function and the nonlinear inequality/equality constraints (if any) in to a single function: essentially, the objective plus a ‘penalty’ for any violated constraints.
This modified objective function is then passed to another optimization algorithm with no nonlinear constraints. If the constraints are violated by the solution of this sub-problem, then the size of the penalties is increased and the process is repeated; eventually, the process must converge to the desired solution (if it exists).
Since all of the actual optimization is performed in this subsidiary optimizer, the subsidiary algorithm that you specify determines whether the optimization is gradient-based or derivative-free.
The local solvers available at the moment are COBYLA'' (for the derivative-free approach) and LBFGS”, MMA'', or SLSQP” (for smooth
functions). The tolerance for the local solver has to be provided.
There is a variant that only uses penalty functions for equality constraints while inequality constraints are passed through to the subsidiary algorithm to be handled directly; in this case, the subsidiary algorithm must handle inequality constraints. (At the moment, this variant has been turned off because of problems with the NLOPT library.)
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
global_solver |
the global NLOPT solver used. |
local_solver |
the local NLOPT solver used, LBFGS or COBYLA. |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
Birgin and Martinez provide their own free implementation of the method as part of the TANGO project; other implementations can be found in semi-free packages like LANCELOT.
Hans W. Borchers
Andrew R. Conn, Nicholas I. M. Gould, and Philippe L. Toint, “A globally convergent augmented Lagrangian algorithm for optimization with general constraints and simple bounds,” SIAM J. Numer. Anal. vol. 28, no. 2, p. 545-572 (1991).
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).
alabama::auglag, Rsolnp::solnp
x0 <- c(1, 1) fn <- function(x) (x[1] - 2) ^ 2 + (x[2] - 1) ^ 2 hin <- function(x) 0.25 * x[1]^2 + x[2] ^ 2 - 1 # hin <= 0 heq <- function(x) x[1] - 2 * x[2] + 1 # heq = 0 gr <- function(x) nl.grad(x, fn) hinjac <- function(x) nl.jacobian(x, hin) heqjac <- function(x) nl.jacobian(x, heq) # with COBYLA auglag(x0, fn, gr = NULL, hin = hin, heq = heq, deprecatedBehavior = FALSE) # $par: 0.8228761 0.9114382 # $value: 1.393464 # $iter: 1001 auglag(x0, fn, gr = NULL, hin = hin, heq = heq, localsolver = "SLSQP", deprecatedBehavior = FALSE) # $par: 0.8228757 0.9114378 # $value: 1.393465 # $iter 184 ## Example from the alabama::auglag help page ## Parameters should be roughly (0, 0, 1) with an objective value of 1. fn <- function(x) (x[1] + 3 * x[2] + x[3]) ^ 2 + 4 * (x[1] - x[2]) ^ 2 heq <- function(x) x[1] + x[2] + x[3] - 1 # hin restated from alabama example to be <= 0. hin <- function(x) c(-6 * x[2] - 4 * x[3] + x[1] ^ 3 + 3, -x[1], -x[2], -x[3]) set.seed(12) auglag(runif(3), fn, hin = hin, heq = heq, localsolver= "lbfgs", deprecatedBehavior = FALSE) # $par: 4.861756e-08 4.732373e-08 9.999999e-01 # $value: 1 # $iter: 145 ## Powell problem from the Rsolnp::solnp help page ## Parameters should be roughly (-1.7171, 1.5957, 1.8272, -0.7636, -0.7636) ## with an objective value of 0.0539498478. x0 <- c(-2, 2, 2, -1, -1) fn1 <- function(x) exp(x[1] * x[2] * x[3] * x[4] * x[5]) eqn1 <-function(x) c(x[1] * x[1] + x[2] * x[2] + x[3] * x[3] + x[4] * x[4] + x[5] * x[5] - 10, x[2] * x[3] - 5 * x[4] * x[5], x[1] * x[1] * x[1] + x[2] * x[2] * x[2] + 1) auglag(x0, fn1, heq = eqn1, localsolver = "mma", deprecatedBehavior = FALSE) # $par: -1.7173645 1.5959655 1.8268352 -0.7636185 -0.7636185 # $value: 0.05394987 # $iter: 916x0 <- c(1, 1) fn <- function(x) (x[1] - 2) ^ 2 + (x[2] - 1) ^ 2 hin <- function(x) 0.25 * x[1]^2 + x[2] ^ 2 - 1 # hin <= 0 heq <- function(x) x[1] - 2 * x[2] + 1 # heq = 0 gr <- function(x) nl.grad(x, fn) hinjac <- function(x) nl.jacobian(x, hin) heqjac <- function(x) nl.jacobian(x, heq) # with COBYLA auglag(x0, fn, gr = NULL, hin = hin, heq = heq, deprecatedBehavior = FALSE) # $par: 0.8228761 0.9114382 # $value: 1.393464 # $iter: 1001 auglag(x0, fn, gr = NULL, hin = hin, heq = heq, localsolver = "SLSQP", deprecatedBehavior = FALSE) # $par: 0.8228757 0.9114378 # $value: 1.393465 # $iter 184 ## Example from the alabama::auglag help page ## Parameters should be roughly (0, 0, 1) with an objective value of 1. fn <- function(x) (x[1] + 3 * x[2] + x[3]) ^ 2 + 4 * (x[1] - x[2]) ^ 2 heq <- function(x) x[1] + x[2] + x[3] - 1 # hin restated from alabama example to be <= 0. hin <- function(x) c(-6 * x[2] - 4 * x[3] + x[1] ^ 3 + 3, -x[1], -x[2], -x[3]) set.seed(12) auglag(runif(3), fn, hin = hin, heq = heq, localsolver= "lbfgs", deprecatedBehavior = FALSE) # $par: 4.861756e-08 4.732373e-08 9.999999e-01 # $value: 1 # $iter: 145 ## Powell problem from the Rsolnp::solnp help page ## Parameters should be roughly (-1.7171, 1.5957, 1.8272, -0.7636, -0.7636) ## with an objective value of 0.0539498478. x0 <- c(-2, 2, 2, -1, -1) fn1 <- function(x) exp(x[1] * x[2] * x[3] * x[4] * x[5]) eqn1 <-function(x) c(x[1] * x[1] + x[2] * x[2] + x[3] * x[3] + x[4] * x[4] + x[5] * x[5] - 10, x[2] * x[3] - 5 * x[4] * x[5], x[1] * x[1] * x[1] + x[2] * x[2] * x[2] + 1) auglag(x0, fn1, heq = eqn1, localsolver = "mma", deprecatedBehavior = FALSE) # $par: -1.7173645 1.5959655 1.8268352 -0.7636185 -0.7636185 # $value: 0.05394987 # $iter: 916
BOBYQA performs derivative-free bound-constrained optimization using an iteratively constructed quadratic approximation for the objective function.
bobyqa( x0, fn, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )bobyqa( x0, fn, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
nl.info |
logical; shall the original NLopt info be shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
This is an algorithm derived from the BOBYQA Fortran subroutine of Powell, converted to C and modified for the NLopt stopping criteria.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
Because BOBYQA constructs a quadratic approximation of the objective, it may perform poorly for objective functions that are not twice-differentiable.
M. J. D. Powell. “The BOBYQA algorithm for bound constrained optimization without derivatives,” Department of Applied Mathematics and Theoretical Physics, Cambridge England, technical reportNA2009/06 (2009).
## Rosenbrock Banana function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} ## The function as written above has a minimum of 0 at (1, 1) S <- bobyqa(c(0, 0), rbf) S ## Rosenbrock Banana function with both parameters constrained to [0, 0.5] S <- bobyqa(c(0, 0), rbf, lower = c(0, 0), upper = c(0.5, 0.5)) S## Rosenbrock Banana function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} ## The function as written above has a minimum of 0 at (1, 1) S <- bobyqa(c(0, 0), rbf) S ## Rosenbrock Banana function with both parameters constrained to [0, 0.5] S <- bobyqa(c(0, 0), rbf, lower = c(0, 0), upper = c(0.5, 0.5)) S
This is a variant of CCSA ("conservative convex separable approximation") which, instead of constructing local MMA approximations, constructs simple quadratic approximations (or rather, affine approximations plus a quadratic penalty term to stay conservative)
ccsaq( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )ccsaq( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
gradient of function |
lower, upper
|
lower and upper bound constraints. |
hin |
function defining the inequality constraints, that is
|
hinjac |
Jacobian of function |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
deprecatedBehavior |
logical; if |
... |
additional arguments passed to the function. |
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 1) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
“Globally convergent” does not mean that this algorithm converges to the global optimum; it means that it is guaranteed to converge to some local minimum from any feasible starting point.
Krister Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM J. Optim. 12 (2), p. 555-573 (2002).
## Solve the Hock-Schittkowski problem no. 100 with analytic gradients ## See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } gr.hs100 <- function(x) { c( 2 * x[1] - 20, 10 * x[2] - 120, 4 * x[3] ^ 3, 6 * x[4] - 66, 60 * x[5] ^ 5, 14 * x[6] - 4 * x[7] - 10, 4 * x[7] ^ 3 - 4 * x[6] - 8) } hinjac.hs100 <- function(x) { matrix(c(4 * x[1], 12 * x[2] ^ 3, 1, 8 * x[4], 5, 0, 0, 7, 3, 20 * x[3], 1, -1, 0, 0, 23, 2 * x[2], 0, 0, 0, 12 * x[6], -8, 8 * x[1] - 3 * x[2], 2 * x[2] - 3 * x[1], 4 * x[3], 0, 0, 5, -11), nrow = 4, byrow = TRUE) } ## The optimum value of the objective function should be 680.6300573 ## A suitable parameter vector is roughly ## (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) # Results with exact Jacobian S <- ccsaq(x0.hs100, fn.hs100, gr = gr.hs100, hin = hin.hs100, hinjac = hinjac.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE) # Results without Jacobian S <- ccsaq(x0.hs100, fn.hs100, hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE)## Solve the Hock-Schittkowski problem no. 100 with analytic gradients ## See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } gr.hs100 <- function(x) { c( 2 * x[1] - 20, 10 * x[2] - 120, 4 * x[3] ^ 3, 6 * x[4] - 66, 60 * x[5] ^ 5, 14 * x[6] - 4 * x[7] - 10, 4 * x[7] ^ 3 - 4 * x[6] - 8) } hinjac.hs100 <- function(x) { matrix(c(4 * x[1], 12 * x[2] ^ 3, 1, 8 * x[4], 5, 0, 0, 7, 3, 20 * x[3], 1, -1, 0, 0, 23, 2 * x[2], 0, 0, 0, 12 * x[6], -8, 8 * x[1] - 3 * x[2], 2 * x[2] - 3 * x[1], 4 * x[3], 0, 0, 5, -11), nrow = 4, byrow = TRUE) } ## The optimum value of the objective function should be 680.6300573 ## A suitable parameter vector is roughly ## (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) # Results with exact Jacobian S <- ccsaq(x0.hs100, fn.hs100, gr = gr.hs100, hin = hin.hs100, hinjac = hinjac.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE) # Results without Jacobian S <- ccsaq(x0.hs100, fn.hs100, hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE)
This function compares the analytic gradients of a function with a finite difference approximation and prints the results of these checks.
check.derivatives( .x, func, func_grad, check_derivatives_tol = 1e-04, check_derivatives_print = "all", func_grad_name = "grad_f", ... )check.derivatives( .x, func, func_grad, check_derivatives_tol = 1e-04, check_derivatives_print = "all", func_grad_name = "grad_f", ... )
.x |
point at which the comparison is done. |
func |
function to be evaluated. |
func_grad |
function calculating the analytic gradients. |
check_derivatives_tol |
option determining when differences between the analytic gradient and its finite difference approximation are flagged as an error. |
check_derivatives_print |
option related to the amount of output. 'all' means that all comparisons are shown, 'errors' only shows comparisons that are flagged as an error, and 'none' shows the number of errors only. |
func_grad_name |
option to change the name of the gradient function that shows up in the output. |
... |
further arguments passed to the functions func and func_grad. |
The return value contains a list with the analytic gradient, its finite difference approximation, the relative errors, and vector comparing the relative errors to the tolerance.
Jelmer Ypma
library('nloptr') # example with correct gradient f <- function(x, a) sum((x - a) ^ 2) f_grad <- function(x, a) 2 * (x - a) check.derivatives(.x = 1:10, func = f, func_grad = f_grad, check_derivatives_print = 'none', a = runif(10)) # example with incorrect gradient f_grad <- function(x, a) 2 * (x - a) + c(0, 0.1, rep(0, 8)) check.derivatives(.x = 1:10, func = f, func_grad = f_grad, check_derivatives_print = 'errors', a = runif(10)) # example with incorrect gradient of vector-valued function g <- function(x, a) c(sum(x - a), sum((x - a) ^ 2)) g_grad <- function(x, a) { rbind(rep(1, length(x)) + c(0, 0.01, rep(0, 8)), 2 * (x - a) + c(0, 0.1, rep(0, 8))) } check.derivatives(.x = 1:10, func = g, func_grad = g_grad, check_derivatives_print = 'all', a = runif(10))library('nloptr') # example with correct gradient f <- function(x, a) sum((x - a) ^ 2) f_grad <- function(x, a) 2 * (x - a) check.derivatives(.x = 1:10, func = f, func_grad = f_grad, check_derivatives_print = 'none', a = runif(10)) # example with incorrect gradient f_grad <- function(x, a) 2 * (x - a) + c(0, 0.1, rep(0, 8)) check.derivatives(.x = 1:10, func = f, func_grad = f_grad, check_derivatives_print = 'errors', a = runif(10)) # example with incorrect gradient of vector-valued function g <- function(x, a) c(sum(x - a), sum((x - a) ^ 2)) g_grad <- function(x, a) { rbind(rep(1, length(x)) + c(0, 0.01, rep(0, 8)), 2 * (x - a) + c(0, 0.1, rep(0, 8))) } check.derivatives(.x = 1:10, func = g, func_grad = g_grad, check_derivatives_print = 'all', a = runif(10))
COBYLA is an algorithm for derivative-free optimization with nonlinear inequality and equality constraints (but see below).
cobyla( x0, fn, lower = NULL, upper = NULL, hin = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )cobyla( x0, fn, lower = NULL, upper = NULL, hin = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
hin |
function defining the inequality constraints, that is
|
nl.info |
logical; shall the original NLopt info be shown. |
control |
list of options, see |
deprecatedBehavior |
logical; if |
... |
additional arguments passed to the function. |
It constructs successive linear approximations of the objective function and
constraints via a simplex of points (in dimensions), and
optimizes these approximations in a trust region at each step.
COBYLA supports equality constraints by transforming them into two
inequality constraints. This functionality has not been added to the wrapper.
To use COBYLA with equality constraints, please use the full
nloptr invocation.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
The original code, written in Fortran by Powell, was converted in C for the SciPy project.
Hans W. Borchers
M. J. D. Powell, “A direct search optimization method that models the objective and constraint functions by linear interpolation,” in Advances in Optimization and Numerical Analysis, eds. S. Gomez and J.-P. Hennart (Kluwer Academic: Dordrecht, 1994), p. 51-67.
## Solve the Hock-Schittkowski problem no. 100 with analytic gradients ## See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } S <- cobyla(x0.hs100, fn.hs100, hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 2000), deprecatedBehavior = FALSE) ## The optimum value of the objective function should be 680.6300573 ## A suitable parameter vector is roughly ## (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) S## Solve the Hock-Schittkowski problem no. 100 with analytic gradients ## See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } S <- cobyla(x0.hs100, fn.hs100, hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 2000), deprecatedBehavior = FALSE) ## The optimum value of the objective function should be 680.6300573 ## A suitable parameter vector is roughly ## (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) S
The Controlled Random Search (CRS) algorithm (and in particular, the CRS2 variant) with the ‘local mutation’ modification.
crs2lm( x0, fn, lower, upper, maxeval = 10000, pop.size = 10 * (length(x0) + 1), ranseed = NULL, xtol_rel = 1e-06, nl.info = FALSE, ... )crs2lm( x0, fn, lower, upper, maxeval = 10000, pop.size = 10 * (length(x0) + 1), ranseed = NULL, xtol_rel = 1e-06, nl.info = FALSE, ... )
x0 |
initial point for searching the optimum. |
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
maxeval |
maximum number of function evaluations. |
pop.size |
population size. |
ranseed |
prescribe seed for random number generator. |
xtol_rel |
stopping criterion for relative change reached. |
nl.info |
logical; shall the original NLopt info be shown. |
... |
additional arguments passed to the function. |
The CRS algorithms are sometimes compared to genetic algorithms, in that they start with a random population of points, and randomly evolve these points by heuristic rules. In this case, the evolution somewhat resembles a randomized Nelder-Mead algorithm.
The published results for CRS seem to be largely empirical.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
The initial population size for CRS defaults to in
dimensions, but this can be changed. The initial population must be
at least .
W. L. Price, “Global optimization by controlled random search,” J. Optim. Theory Appl. 40 (3), p. 333-348 (1983).
P. Kaelo and M. M. Ali, “Some variants of the controlled random search algorithm for global optimization,” J. Optim. Theory Appl. 130 (2), 253-264 (2006).
## Minimize the Hartmann 6-Dimensional function ## See https://www.sfu.ca/~ssurjano/hart6.html a <- c(1.0, 1.2, 3.0, 3.2) A <- matrix(c(10, 0.05, 3, 17, 3, 10, 3.5, 8, 17, 17, 1.7, 0.05, 3.5, 0.1, 10, 10, 1.7, 8, 17, 0.1, 8, 14, 8, 14), nrow = 4) B <- matrix(c(.1312, .2329, .2348, .4047, .1696, .4135, .1451, .8828, .5569, .8307, .3522, .8732, .0124, .3736, .2883, .5743, .8283, .1004, .3047, .1091, .5886, .9991, .6650, .0381), nrow = 4) hartmann6 <- function(x, a, A, B) { fun <- 0 for (i in 1:4) { fun <- fun - a[i] * exp(-sum(A[i, ] * (x - B[i, ]) ^ 2)) } fun } ## The function has a global minimum of -3.32237 at ## (0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573) S <- crs2lm(x0 = rep(0, 6), hartmann6, lower = rep(0, 6), upper = rep(1, 6), ranseed = 10L, nl.info = TRUE, xtol_rel=1e-8, maxeval = 10000, a = a, A = A, B = B) S## Minimize the Hartmann 6-Dimensional function ## See https://www.sfu.ca/~ssurjano/hart6.html a <- c(1.0, 1.2, 3.0, 3.2) A <- matrix(c(10, 0.05, 3, 17, 3, 10, 3.5, 8, 17, 17, 1.7, 0.05, 3.5, 0.1, 10, 10, 1.7, 8, 17, 0.1, 8, 14, 8, 14), nrow = 4) B <- matrix(c(.1312, .2329, .2348, .4047, .1696, .4135, .1451, .8828, .5569, .8307, .3522, .8732, .0124, .3736, .2883, .5743, .8283, .1004, .3047, .1091, .5886, .9991, .6650, .0381), nrow = 4) hartmann6 <- function(x, a, A, B) { fun <- 0 for (i in 1:4) { fun <- fun - a[i] * exp(-sum(A[i, ] * (x - B[i, ]) ^ 2)) } fun } ## The function has a global minimum of -3.32237 at ## (0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573) S <- crs2lm(x0 = rep(0, 6), hartmann6, lower = rep(0, 6), upper = rep(1, 6), ranseed = 10L, nl.info = TRUE, xtol_rel=1e-8, maxeval = 10000, a = a, A = A, B = B) S
DIRECT is a deterministic search algorithm based on systematic division of the search domain into smaller and smaller hyperrectangles. The DIRECT_L makes the algorithm more biased towards local search (more efficient for functions without too many minima).
direct( fn, lower, upper, scaled = TRUE, original = FALSE, nl.info = FALSE, control = list(), ... ) directL( fn, lower, upper, randomized = FALSE, original = FALSE, nl.info = FALSE, control = list(), ... )direct( fn, lower, upper, scaled = TRUE, original = FALSE, nl.info = FALSE, control = list(), ... ) directL( fn, lower, upper, randomized = FALSE, original = FALSE, nl.info = FALSE, control = list(), ... )
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
scaled |
logical; shall the hypercube be scaled before starting. |
original |
logical; whether to use the original implementation by Gablonsky – the performance is mostly similar. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
randomized |
logical; shall some randomization be used to decide which dimension to halve next in the case of near-ties. |
The DIRECT and DIRECT-L algorithms start by rescaling the bound constraints to a hypercube, which gives all dimensions equal weight in the search procedure. If your dimensions do not have equal weight, e.g. if you have a “long and skinny” search space and your function varies at about the same speed in all directions, it may be better to use unscaled variant of the DIRECT algorithm.
The algorithms only handle finite bound constraints which must be provided. The original versions may include some support for arbitrary nonlinear inequality, but this has not been tested.
The original versions do not have randomized or unscaled variants, so these options will be disregarded for these versions.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
The DIRECT_L algorithm should be tried first.
Hans W. Borchers
D. R. Jones, C. D. Perttunen, and B. E. Stuckmann, “Lipschitzian optimization without the Lipschitz constant,” J. Optimization Theory and Applications, vol. 79, p. 157 (1993).
J. M. Gablonsky and C. T. Kelley, “A locally-biased form of the DIRECT algorithm," J. Global Optimization, vol. 21 (1), p. 27-37 (2001).
The dfoptim package will provide a pure R version of this
algorithm.
### Minimize the Hartmann6 function hartmann6 <- function(x) { a <- c(1.0, 1.2, 3.0, 3.2) A <- matrix(c(10.0, 0.05, 3.0, 17.0, 3.0, 10.0, 3.5, 8.0, 17.0, 17.0, 1.7, 0.05, 3.5, 0.1, 10.0, 10.0, 1.7, 8.0, 17.0, 0.1, 8.0, 14.0, 8.0, 14.0), nrow=4, ncol=6) B <- matrix(c(.1312,.2329,.2348,.4047, .1696,.4135,.1451,.8828, .5569,.8307,.3522,.8732, .0124,.3736,.2883,.5743, .8283,.1004,.3047,.1091, .5886,.9991,.6650,.0381), nrow=4, ncol=6) fun <- 0 for (i in 1:4) { fun <- fun - a[i] * exp(-sum(A[i,] * (x - B[i,]) ^ 2)) } fun } S <- directL(hartmann6, rep(0, 6), rep(1, 6), nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 1000)) ## Number of Iterations....: 1000 ## Termination conditions: stopval: -Inf ## xtol_rel: 1e-08, maxeval: 1000, ftol_rel: 0, ftol_abs: 0 ## Number of inequality constraints: 0 ## Number of equality constraints: 0 ## Current value of objective function: -3.32236800687327 ## Current value of controls: ## 0.2016884 0.1500025 0.4768667 0.2753391 0.311648 0.6572931### Minimize the Hartmann6 function hartmann6 <- function(x) { a <- c(1.0, 1.2, 3.0, 3.2) A <- matrix(c(10.0, 0.05, 3.0, 17.0, 3.0, 10.0, 3.5, 8.0, 17.0, 17.0, 1.7, 0.05, 3.5, 0.1, 10.0, 10.0, 1.7, 8.0, 17.0, 0.1, 8.0, 14.0, 8.0, 14.0), nrow=4, ncol=6) B <- matrix(c(.1312,.2329,.2348,.4047, .1696,.4135,.1451,.8828, .5569,.8307,.3522,.8732, .0124,.3736,.2883,.5743, .8283,.1004,.3047,.1091, .5886,.9991,.6650,.0381), nrow=4, ncol=6) fun <- 0 for (i in 1:4) { fun <- fun - a[i] * exp(-sum(A[i,] * (x - B[i,]) ^ 2)) } fun } S <- directL(hartmann6, rep(0, 6), rep(1, 6), nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 1000)) ## Number of Iterations....: 1000 ## Termination conditions: stopval: -Inf ## xtol_rel: 1e-08, maxeval: 1000, ftol_rel: 0, ftol_abs: 0 ## Number of inequality constraints: 0 ## Number of equality constraints: 0 ## Current value of objective function: -3.32236800687327 ## Current value of controls: ## 0.2016884 0.1500025 0.4768667 0.2753391 0.311648 0.6572931
is.nloptr preforms checks to see if a fully specified problem is supplied to nloptr. Mostly for internal use.
is.nloptr(x)is.nloptr(x)
x |
object to be tested. |
Logical. Return TRUE if all tests were passed, otherwise return FALSE or exit with Error.
Jelmer Ypma
The Improved Stochastic Ranking Evolution Strategy (ISRES) is an algorithm for nonlinearly constrained global optimization, or at least semi-global, although it has heuristics to escape local optima.
isres( x0, fn, lower, upper, hin = NULL, heq = NULL, maxeval = 10000, pop.size = 20 * (length(x0) + 1), xtol_rel = 1e-06, nl.info = FALSE, deprecatedBehavior = TRUE, ... )isres( x0, fn, lower, upper, hin = NULL, heq = NULL, maxeval = 10000, pop.size = 20 * (length(x0) + 1), xtol_rel = 1e-06, nl.info = FALSE, deprecatedBehavior = TRUE, ... )
x0 |
initial point for searching the optimum. |
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
hin |
function defining the inequality constraints, that is
|
heq |
function defining the equality constraints, that is |
maxeval |
maximum number of function evaluations. |
pop.size |
population size. |
xtol_rel |
stopping criterion for relative change reached. |
nl.info |
logical; shall the original NLopt info be shown. |
deprecatedBehavior |
logical; if |
... |
additional arguments passed to the function. |
The evolution strategy is based on a combination of a mutation rule—with a log-normal step-size update and exponential smoothing—and differential variation—a Nelder-Mead-like update rule). The fitness ranking is simply via the objective function for problems without nonlinear constraints, but when nonlinear constraints are included the stochastic ranking proposed by Runarsson and Yao is employed.
This method supports arbitrary nonlinear inequality and equality constraints in addition to the bounds constraints.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
The initial population size for CRS defaults to in
dimensions, but this can be changed. The initial population must be
at least .
Hans W. Borchers
Thomas Philip Runarsson and Xin Yao, “Search biases in constrained evolutionary optimization,” IEEE Trans. on Systems, Man, and Cybernetics Part C: Applications and Reviews, vol. 35 (no. 2), pp. 233-243 (2005).
## Rosenbrock Banana objective function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} x0 <- c(-1.2, 1) lb <- c(-3, -3) ub <- c(3, 3) ## The function as written above has a minimum of 0 at (1, 1) isres(x0 = x0, fn = rbf, lower = lb, upper = ub) ## Now subject to the inequality that x[1] + x[2] <= 1.5 hin <- function(x) {x[1] + x[2] - 1.5} S <- isres(x0 = x0, fn = rbf, hin = hin, lower = lb, upper = ub, maxeval = 2e5L, deprecatedBehavior = FALSE) S sum(S$par)## Rosenbrock Banana objective function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} x0 <- c(-1.2, 1) lb <- c(-3, -3) ub <- c(3, 3) ## The function as written above has a minimum of 0 at (1, 1) isres(x0 = x0, fn = rbf, lower = lb, upper = ub) ## Now subject to the inequality that x[1] + x[2] <= 1.5 hin <- function(x) {x[1] + x[2] - 1.5} S <- isres(x0 = x0, fn = rbf, hin = hin, lower = lb, upper = ub, maxeval = 2e5L, deprecatedBehavior = FALSE) S sum(S$par)
Low-storage version of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method.
lbfgs( x0, fn, gr = NULL, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )lbfgs( x0, fn, gr = NULL, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )
x0 |
initial point for searching the optimum. |
fn |
objective function to be minimized. |
gr |
gradient of function |
lower, upper
|
lower and upper bound constraints. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of control parameters, see |
... |
further arguments to be passed to the function. |
The low-storage (or limited-memory) algorithm is a member of the class of quasi-Newton optimization methods. It is well suited for optimization problems with a large number of variables.
One parameter of this algorithm is the number m of gradients to
remember from previous optimization steps. NLopt sets m to a
heuristic value by default. It can be changed by the NLopt function
set_vector_storage.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
Based on a Fortran implementation of the low-storage BFGS algorithm written by L. Luksan, and posted under the GNU LGPL license.
Hans W. Borchers
J. Nocedal, “Updating quasi-Newton matrices with limited storage,” Math. Comput. 35, 773-782 (1980).
D. C. Liu and J. Nocedal, “On the limited memory BFGS method for large scale optimization,” Math. Programming 45, p. 503-528 (1989).
flb <- function(x) { p <- length(x) sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) } # 25-dimensional box constrained: par[24] is *not* at the boundary S <- lbfgs(rep(3, 25), flb, lower=rep(2, 25), upper=rep(4, 25), nl.info = TRUE, control = list(xtol_rel=1e-8)) ## Optimal value of objective function: 368.105912874334 ## Optimal value of controls: 2 ... 2 2.109093 4flb <- function(x) { p <- length(x) sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) } # 25-dimensional box constrained: par[24] is *not* at the boundary S <- lbfgs(rep(3, 25), flb, lower=rep(2, 25), upper=rep(4, 25), nl.info = TRUE, control = list(xtol_rel=1e-8)) ## Optimal value of objective function: 368.105912874334 ## Optimal value of controls: 2 ... 2 2.109093 4
The “Multi-Level Single-Linkage” (MLSL) algorithm for global optimization searches by a sequence of local optimizations from random starting points. A modification of MLSL is included using a low-discrepancy sequence (LDS) instead of pseudorandom numbers.
mlsl( x0, fn, gr = NULL, lower, upper, local.method = "LBFGS", low.discrepancy = TRUE, nl.info = FALSE, control = list(), ... )mlsl( x0, fn, gr = NULL, lower, upper, local.method = "LBFGS", low.discrepancy = TRUE, nl.info = FALSE, control = list(), ... )
x0 |
initial point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
gradient of function |
lower, upper
|
lower and upper bound constraints. |
local.method |
only |
low.discrepancy |
logical; shall a low discrepancy variation be used. |
nl.info |
logical; shall the original NLopt info be shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
MLSL is a ‘multistart’ algorithm: it works by doing a sequence of local optimizations—using some other local optimization algorithm—from random or low-discrepancy starting points. MLSL is distinguished, however, by a ‘clustering’ heuristic that helps it to avoid repeated searches of the same local optima and also has some theoretical guarantees of finding all local optima in a finite number of local minimizations.
The local-search portion of MLSL can use any of the other algorithms in NLopt, and, in particular, can use either gradient-based or derivative-free algorithms. For this wrapper only gradient-based LBFGS is available as local method.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
If you don't set a stopping tolerance for your local-optimization
algorithm, MLSL defaults to ftol_rel = 1e-15 and
xtol_rel = 1e-7 for the local searches.
Hans W. Borchers
A. H. G. Rinnooy Kan and G. T. Timmer, “Stochastic global optimization methods” Mathematical Programming, vol. 39, p. 27-78 (1987).
Sergei Kucherenko and Yury Sytsko, “Application of deterministic low-discrepancy sequences in global optimization”, Computational Optimization and Applications, vol. 30, p. 297-318 (2005).
## Minimize the Hartmann 6-Dimensional function ## See https://www.sfu.ca/~ssurjano/hart6.html a <- c(1.0, 1.2, 3.0, 3.2) A <- matrix(c(10, 0.05, 3, 17, 3, 10, 3.5, 8, 17, 17, 1.7, 0.05, 3.5, 0.1, 10, 10, 1.7, 8, 17, 0.1, 8, 14, 8, 14), nrow = 4) B <- matrix(c(.1312, .2329, .2348, .4047, .1696, .4135, .1451, .8828, .5569, .8307, .3522, .8732, .0124, .3736, .2883, .5743, .8283, .1004, .3047, .1091, .5886, .9991, .6650, .0381), nrow = 4) hartmann6 <- function(x, a, A, B) { fun <- 0 for (i in 1:4) { fun <- fun - a[i] * exp(-sum(A[i, ] * (x - B[i, ]) ^ 2)) } fun } ## The function has a global minimum of -3.32237 at ## (0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573) S <- mlsl(x0 = rep(0, 6), hartmann6, lower = rep(0, 6), upper = rep(1, 6), nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 1000), a = a, A = A, B = B)## Minimize the Hartmann 6-Dimensional function ## See https://www.sfu.ca/~ssurjano/hart6.html a <- c(1.0, 1.2, 3.0, 3.2) A <- matrix(c(10, 0.05, 3, 17, 3, 10, 3.5, 8, 17, 17, 1.7, 0.05, 3.5, 0.1, 10, 10, 1.7, 8, 17, 0.1, 8, 14, 8, 14), nrow = 4) B <- matrix(c(.1312, .2329, .2348, .4047, .1696, .4135, .1451, .8828, .5569, .8307, .3522, .8732, .0124, .3736, .2883, .5743, .8283, .1004, .3047, .1091, .5886, .9991, .6650, .0381), nrow = 4) hartmann6 <- function(x, a, A, B) { fun <- 0 for (i in 1:4) { fun <- fun - a[i] * exp(-sum(A[i, ] * (x - B[i, ]) ^ 2)) } fun } ## The function has a global minimum of -3.32237 at ## (0.20169, 0.150011, 0.476874, 0.275332, 0.311652, 0.6573) S <- mlsl(x0 = rep(0, 6), hartmann6, lower = rep(0, 6), upper = rep(1, 6), nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 1000), a = a, A = A, B = B)
Globally-convergent method-of-moving-asymptotes (MMA) algorithm for gradient-based local optimization, including nonlinear inequality constraints (but not equality constraints).
mma( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )mma( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
gradient of function |
lower, upper
|
lower and upper bound constraints. |
hin |
function defining the inequality constraints, that is
|
hinjac |
Jacobian of function |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
deprecatedBehavior |
logical; if |
... |
additional arguments passed to the function. |
This is an improved CCSA ("conservative convex separable approximation") variant of the original MMA algorithm published by Svanberg in 1987, which has become popular for topology optimization.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 1) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
“Globally convergent” does not mean that this algorithm converges to the global optimum; rather, it means that the algorithm is guaranteed to converge to some local minimum from any feasible starting point.
Hans W. Borchers
Krister Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations”, SIAM J. Optim. 12 (2), p. 555-573 (2002).
# Solve the Hock-Schittkowski problem no. 100 with analytic gradients # See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } gr.hs100 <- function(x) { c( 2 * x[1] - 20, 10 * x[2] - 120, 4 * x[3] ^ 3, 6 * x[4] - 66, 60 * x[5] ^ 5, 14 * x[6] - 4 * x[7] - 10, 4 * x[7] ^ 3 - 4 * x[6] - 8) } hinjac.hs100 <- function(x) { matrix(c(4 * x[1], 12 * x[2] ^ 3, 1, 8 * x[4], 5, 0, 0, 7, 3, 20 * x[3], 1, -1, 0, 0, 23, 2 * x[2], 0, 0, 0, 12 * x[6], -8, 8 * x[1] - 3 * x[2], 2 * x[2] - 3 * x[1], 4 * x[3], 0, 0, 5, -11), nrow = 4, byrow = TRUE) } # The optimum value of the objective function should be 680.6300573 # A suitable parameter vector is roughly # (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) # Using analytic Jacobian S <- mma(x0.hs100, fn.hs100, gr = gr.hs100, hin = hin.hs100, hinjac = hinjac.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE) # Using computed Jacobian S <- mma(x0.hs100, fn.hs100, hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE)# Solve the Hock-Schittkowski problem no. 100 with analytic gradients # See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } gr.hs100 <- function(x) { c( 2 * x[1] - 20, 10 * x[2] - 120, 4 * x[3] ^ 3, 6 * x[4] - 66, 60 * x[5] ^ 5, 14 * x[6] - 4 * x[7] - 10, 4 * x[7] ^ 3 - 4 * x[6] - 8) } hinjac.hs100 <- function(x) { matrix(c(4 * x[1], 12 * x[2] ^ 3, 1, 8 * x[4], 5, 0, 0, 7, 3, 20 * x[3], 1, -1, 0, 0, 23, 2 * x[2], 0, 0, 0, 12 * x[6], -8, 8 * x[1] - 3 * x[2], 2 * x[2] - 3 * x[1], 4 * x[3], 0, 0, 5, -11), nrow = 4, byrow = TRUE) } # The optimum value of the objective function should be 680.6300573 # A suitable parameter vector is roughly # (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) # Using analytic Jacobian S <- mma(x0.hs100, fn.hs100, gr = gr.hs100, hin = hin.hs100, hinjac = hinjac.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE) # Using computed Jacobian S <- mma(x0.hs100, fn.hs100, hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8), deprecatedBehavior = FALSE)
An implementation of almost the original Nelder-Mead simplex algorithm.
neldermead( x0, fn, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )neldermead( x0, fn, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
Provides explicit support for bound constraints, using essentially the method proposed in Box. Whenever a new point would lie outside the bound constraints the point is moved back exactly onto the constraint.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
The author of NLopt would tend to recommend the Subplex method instead.
Hans W. Borchers
J. A. Nelder and R. Mead, “A simplex method for function minimization,” The Computer Journal 7, p. 308-313 (1965).
M. J. Box, “A new method of constrained optimization and a comparison with other methods,” Computer J. 8 (1), 42-52 (1965).
dfoptim::nmk
# Fletcher and Powell's helic valley fphv <- function(x) 100*(x[3] - 10*atan2(x[2], x[1])/(2*pi))^2 + (sqrt(x[1]^2 + x[2]^2) - 1)^2 +x[3]^2 x0 <- c(-1, 0, 0) neldermead(x0, fphv) # 1 0 0 # Powell's Singular Function (PSF) psf <- function(x) (x[1] + 10*x[2])^2 + 5*(x[3] - x[4])^2 + (x[2] - 2*x[3])^4 + 10*(x[1] - x[4])^4 x0 <- c(3, -1, 0, 1) neldermead(x0, psf) # 0 0 0 0, needs maximum number of function calls ## Not run: # Bounded version of Nelder-Mead rosenbrock <- function(x) { ## Rosenbrock Banana function 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2 + 100 * (x[3] - x[2]^2)^2 + (1 - x[2])^2 } lower <- c(-Inf, 0, 0) upper <- c( Inf, 0.5, 1) x0 <- c(0, 0.1, 0.1) S <- neldermead(c(0, 0.1, 0.1), rosenbrock, lower, upper, nl.info = TRUE) # $xmin = c(0.7085595, 0.5000000, 0.2500000) # $fmin = 0.3353605 ## End(Not run)# Fletcher and Powell's helic valley fphv <- function(x) 100*(x[3] - 10*atan2(x[2], x[1])/(2*pi))^2 + (sqrt(x[1]^2 + x[2]^2) - 1)^2 +x[3]^2 x0 <- c(-1, 0, 0) neldermead(x0, fphv) # 1 0 0 # Powell's Singular Function (PSF) psf <- function(x) (x[1] + 10*x[2])^2 + 5*(x[3] - x[4])^2 + (x[2] - 2*x[3])^4 + 10*(x[1] - x[4])^4 x0 <- c(3, -1, 0, 1) neldermead(x0, psf) # 0 0 0 0, needs maximum number of function calls ## Not run: # Bounded version of Nelder-Mead rosenbrock <- function(x) { ## Rosenbrock Banana function 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2 + 100 * (x[3] - x[2]^2)^2 + (1 - x[2])^2 } lower <- c(-Inf, 0, 0) upper <- c( Inf, 0.5, 1) x0 <- c(0, 0.1, 0.1) S <- neldermead(c(0, 0.1, 0.1), rosenbrock, lower, upper, nl.info = TRUE) # $xmin = c(0.7085595, 0.5000000, 0.2500000) # $fmin = 0.3353605 ## End(Not run)
NEWUOA solves quadratic subproblems in a spherical trust region via a truncated conjugate-gradient algorithm. For bound-constrained problems, BOBYQA should be used instead, as Powell developed it as an enhancement thereof for bound constraints.
newuoa(x0, fn, nl.info = FALSE, control = list(), ...)newuoa(x0, fn, nl.info = FALSE, control = list(), ...)
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
nl.info |
logical; shall the original NLopt info be shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
This is an algorithm derived from the NEWUOA Fortran subroutine of Powell, converted to C and modified for the NLopt stopping criteria.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
NEWUOA may be largely superseded by BOBYQA.
Hans W. Borchers
M. J. D. Powell. “The BOBYQA algorithm for bound constrained optimization without derivatives,” Department of Applied Mathematics and Theoretical Physics, Cambridge England, technical reportNA2009/06 (2009).
## Rosenbrock Banana function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} S <- newuoa(c(1, 2), rbf) ## The function as written above has a minimum of 0 at (1, 1) S## Rosenbrock Banana function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} S <- newuoa(c(1, 2), rbf) ## The function as written above has a minimum of 0 at (1, 1) S
Provides numerical gradients and Jacobians.
nl.grad(x0, fn, heps = .Machine$double.eps^(1/3), ...)nl.grad(x0, fn, heps = .Machine$double.eps^(1/3), ...)
x0 |
point as a vector where the gradient is to be calculated. |
fn |
scalar function of one or several variables. |
heps |
step size to be used. |
... |
additional arguments passed to the function. |
Both functions apply the “central difference formula” with step size as recommended in the literature.
grad returns the gradient as a vector; jacobian
returns the Jacobian as a matrix of usual dimensions.
Hans W. Borchers
fn1 <- function(x) sum(x ^ 2) nl.grad(seq(0, 1, by = 0.2), fn1) ## [1] 0.0 0.4 0.8 1.2 1.6 2.0 nl.grad(rep(1, 5), fn1) ## [1] 2 2 2 2 2 fn2 <- function(x) c(sin(x), cos(x)) x <- (0:1) * 2 * pi nl.jacobian(x, fn2) ## [,1] [,2] ## [1,] 1 0 ## [2,] 0 1 ## [3,] 0 0 ## [4,] 0 0fn1 <- function(x) sum(x ^ 2) nl.grad(seq(0, 1, by = 0.2), fn1) ## [1] 0.0 0.4 0.8 1.2 1.6 2.0 nl.grad(rep(1, 5), fn1) ## [1] 2 2 2 2 2 fn2 <- function(x) c(sin(x), cos(x)) x <- (0:1) * 2 * pi nl.jacobian(x, fn2) ## [,1] [,2] ## [1,] 1 0 ## [2,] 0 1 ## [3,] 0 0 ## [4,] 0 0
Sets and changes the NLOPT options.
nl.opts(optlist = NULL)nl.opts(optlist = NULL)
optlist |
list of options, see below. |
The following options can be set (here with default values):
stopval = -Inf, # stop minimization at this valuextol_rel =
1e-6, # stop on small optimization stepmaxeval = 1000, # stop on
this many function evaluationsftol_rel = 0.0, # stop on change
times function valueftol_abs = 0.0, # stop on small change of
function valuecheck_derivatives = FALSE
returns a list with default and changed options.
There are more options that can be set for solvers in NLOPT. These
cannot be set through their wrapper functions. To see the full list of
options and algorithms, type nloptr.print.options().
Hans W. Borchers
nl.opts(list(xtol_rel = 1e-8, maxeval = 2000))nl.opts(list(xtol_rel = 1e-8, maxeval = 2000))
nloptr is an R interface to NLopt, a free/open-source library for nonlinear optimization started by Steven G. Johnson, providing a common interface for a number of different free optimization routines available online as well as original implementations of various other algorithms. The NLopt library is available under the GNU Lesser General Public License (LGPL), and the copyrights are owned by a variety of authors. Most of the information here has been taken from the NLopt website, where more details are available.
nloptr( x0, eval_f, eval_grad_f = NULL, lb = NULL, ub = NULL, eval_g_ineq = NULL, eval_jac_g_ineq = NULL, eval_g_eq = NULL, eval_jac_g_eq = NULL, opts = list(), ... )nloptr( x0, eval_f, eval_grad_f = NULL, lb = NULL, ub = NULL, eval_g_ineq = NULL, eval_jac_g_ineq = NULL, eval_g_eq = NULL, eval_jac_g_eq = NULL, opts = list(), ... )
x0 |
vector with starting values for the optimization. |
||||||||
eval_f |
function that returns the value of the objective function. It can also return gradient information at the same time in a list with elements "objective" and "gradient" (see below for an example). |
||||||||
eval_grad_f |
function that returns the value of the gradient of the objective function. Not all of the algorithms require a gradient. |
||||||||
lb |
vector with lower bounds of the controls (use |
||||||||
ub |
vector with upper bounds of the controls (use |
||||||||
eval_g_ineq |
function to evaluate (non-)linear inequality constraints that should hold in the solution. It can also return gradient information at the same time in a list with elements "constraints" and "jacobian" (see below for an example). |
||||||||
eval_jac_g_ineq |
function to evaluate the Jacobian of the (non-)linear inequality constraints that should hold in the solution. |
||||||||
eval_g_eq |
function to evaluate (non-)linear equality constraints that should hold in the solution. It can also return gradient information at the same time in a list with elements "constraints" and "jacobian" (see below for an example). |
||||||||
eval_jac_g_eq |
function to evaluate the Jacobian of the (non-)linear equality constraints that should hold in the solution. |
||||||||
opts |
list with options. The option " Some algorithms with equality constraints require the option
The option
The option |
||||||||
... |
arguments that will be passed to the user-defined objective and constraints functions. |
NLopt addresses general nonlinear optimization problems of the form:
where is the objective function to be minimized and
represents the optimization parameters. This problem may optionally
be subject to the bound constraints (also called box constraints),
and . For partially or totally unconstrained problems the bounds can
take -Inf or Inf. One may also optionally have
nonlinear inequality constraints (sometimes called a nonlinear programming
problem), which can be specified in , and equality constraints that
can be specified in . Note that not all of the algorithms in NLopt
can handle constraints.
The return value contains a list with the inputs, and additional elements
call |
the call that was made to solve |
status |
integer value with the status of the optimization (0 is success) |
message |
more informative message with the status of the optimization |
iterations |
number of iterations that were executed |
objective |
value if the objective function in the solution |
solution |
optimal value of the controls |
version |
version of NLopt that was used |
See ?`nloptr-package` for an extended example.
Steven G. Johnson and others (C code)
Jelmer Ypma (R interface)
Steven G. Johnson, The NLopt nonlinear-optimization package, https://github.com/stevengj/nlopt
nloptr.print.options
check.derivatives
optim
nlm
nlminb
Rsolnp::Rsolnp
Rsolnp::solnp
library('nloptr') ## Rosenbrock Banana function and gradient in separate functions eval_f <- function(x) { return(100 * (x[2] - x[1] * x[1])^2 + (1 - x[1])^2) } eval_grad_f <- function(x) { return(c(-400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 * (x[2] - x[1] * x[1]))) } # initial values x0 <- c(-1.2, 1) opts <- list("algorithm"="NLOPT_LD_LBFGS", "xtol_rel"=1.0e-8) # solve Rosenbrock Banana function res <- nloptr(x0=x0, eval_f=eval_f, eval_grad_f=eval_grad_f, opts=opts) print(res) ## Rosenbrock Banana function and gradient in one function # this can be used to economize on calculations eval_f_list <- function(x) { return( list( "objective" = 100 * (x[2] - x[1] * x[1]) ^ 2 + (1 - x[1]) ^ 2, "gradient" = c(-400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 * (x[2] - x[1] * x[1])))) } # solve Rosenbrock Banana function using an objective function that # returns a list with the objective value and its gradient res <- nloptr(x0=x0, eval_f=eval_f_list, opts=opts) print(res) # Example showing how to solve the problem from the NLopt tutorial. # # min sqrt(x2) # s.t. x2 >= 0 # x2 >= (a1*x1 + b1)^3 # x2 >= (a2*x1 + b2)^3 # where # a1 = 2, b1 = 0, a2 = -1, b2 = 1 # # re-formulate constraints to be of form g(x) <= 0 # (a1*x1 + b1)^3 - x2 <= 0 # (a2*x1 + b2)^3 - x2 <= 0 library('nloptr') # objective function eval_f0 <- function(x, a, b) { return(sqrt(x[2])) } # constraint function eval_g0 <- function(x, a, b) { return((a*x[1] + b)^3 - x[2]) } # gradient of objective function eval_grad_f0 <- function(x, a, b) { return(c(0, .5/sqrt(x[2]))) } # Jacobian of constraint eval_jac_g0 <- function(x, a, b) { return(rbind(c(3*a[1]*(a[1]*x[1] + b[1])^2, -1.0), c(3*a[2]*(a[2]*x[1] + b[2])^2, -1.0))) } # functions with gradients in objective and constraint function # this can be useful if the same calculations are needed for # the function value and the gradient eval_f1 <- function(x, a, b) { return(list("objective"=sqrt(x[2]), "gradient"=c(0,.5/sqrt(x[2])))) } eval_g1 <- function(x, a, b) { return(list("constraints"=(a*x[1] + b)^3 - x[2], "jacobian"=rbind(c(3*a[1]*(a[1]*x[1] + b[1])^2, -1.0), c(3*a[2]*(a[2]*x[1] + b[2])^2, -1.0)))) } # define parameters a <- c(2,-1) b <- c(0, 1) # Solve using NLOPT_LD_MMA with gradient information supplied in separate # function. res0 <- nloptr(x0=c(1.234,5.678), eval_f=eval_f0, eval_grad_f=eval_grad_f0, lb = c(-Inf,0), ub = c(Inf,Inf), eval_g_ineq = eval_g0, eval_jac_g_ineq = eval_jac_g0, opts = list("algorithm"="NLOPT_LD_MMA"), a = a, b = b) print(res0) # Solve using NLOPT_LN_COBYLA without gradient information res1 <- nloptr(x0=c(1.234,5.678), eval_f=eval_f0, lb = c(-Inf, 0), ub = c(Inf, Inf), eval_g_ineq = eval_g0, opts = list("algorithm" = "NLOPT_LN_COBYLA"), a = a, b = b) print(res1) # Solve using NLOPT_LD_MMA with gradient information in objective function res2 <- nloptr(x0=c(1.234, 5.678), eval_f=eval_f1, lb = c(-Inf, 0), ub = c(Inf, Inf), eval_g_ineq = eval_g1, opts = list("algorithm"="NLOPT_LD_MMA", "check_derivatives" = TRUE), a = a, b = b) print(res2)library('nloptr') ## Rosenbrock Banana function and gradient in separate functions eval_f <- function(x) { return(100 * (x[2] - x[1] * x[1])^2 + (1 - x[1])^2) } eval_grad_f <- function(x) { return(c(-400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 * (x[2] - x[1] * x[1]))) } # initial values x0 <- c(-1.2, 1) opts <- list("algorithm"="NLOPT_LD_LBFGS", "xtol_rel"=1.0e-8) # solve Rosenbrock Banana function res <- nloptr(x0=x0, eval_f=eval_f, eval_grad_f=eval_grad_f, opts=opts) print(res) ## Rosenbrock Banana function and gradient in one function # this can be used to economize on calculations eval_f_list <- function(x) { return( list( "objective" = 100 * (x[2] - x[1] * x[1]) ^ 2 + (1 - x[1]) ^ 2, "gradient" = c(-400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 200 * (x[2] - x[1] * x[1])))) } # solve Rosenbrock Banana function using an objective function that # returns a list with the objective value and its gradient res <- nloptr(x0=x0, eval_f=eval_f_list, opts=opts) print(res) # Example showing how to solve the problem from the NLopt tutorial. # # min sqrt(x2) # s.t. x2 >= 0 # x2 >= (a1*x1 + b1)^3 # x2 >= (a2*x1 + b2)^3 # where # a1 = 2, b1 = 0, a2 = -1, b2 = 1 # # re-formulate constraints to be of form g(x) <= 0 # (a1*x1 + b1)^3 - x2 <= 0 # (a2*x1 + b2)^3 - x2 <= 0 library('nloptr') # objective function eval_f0 <- function(x, a, b) { return(sqrt(x[2])) } # constraint function eval_g0 <- function(x, a, b) { return((a*x[1] + b)^3 - x[2]) } # gradient of objective function eval_grad_f0 <- function(x, a, b) { return(c(0, .5/sqrt(x[2]))) } # Jacobian of constraint eval_jac_g0 <- function(x, a, b) { return(rbind(c(3*a[1]*(a[1]*x[1] + b[1])^2, -1.0), c(3*a[2]*(a[2]*x[1] + b[2])^2, -1.0))) } # functions with gradients in objective and constraint function # this can be useful if the same calculations are needed for # the function value and the gradient eval_f1 <- function(x, a, b) { return(list("objective"=sqrt(x[2]), "gradient"=c(0,.5/sqrt(x[2])))) } eval_g1 <- function(x, a, b) { return(list("constraints"=(a*x[1] + b)^3 - x[2], "jacobian"=rbind(c(3*a[1]*(a[1]*x[1] + b[1])^2, -1.0), c(3*a[2]*(a[2]*x[1] + b[2])^2, -1.0)))) } # define parameters a <- c(2,-1) b <- c(0, 1) # Solve using NLOPT_LD_MMA with gradient information supplied in separate # function. res0 <- nloptr(x0=c(1.234,5.678), eval_f=eval_f0, eval_grad_f=eval_grad_f0, lb = c(-Inf,0), ub = c(Inf,Inf), eval_g_ineq = eval_g0, eval_jac_g_ineq = eval_jac_g0, opts = list("algorithm"="NLOPT_LD_MMA"), a = a, b = b) print(res0) # Solve using NLOPT_LN_COBYLA without gradient information res1 <- nloptr(x0=c(1.234,5.678), eval_f=eval_f0, lb = c(-Inf, 0), ub = c(Inf, Inf), eval_g_ineq = eval_g0, opts = list("algorithm" = "NLOPT_LN_COBYLA"), a = a, b = b) print(res1) # Solve using NLOPT_LD_MMA with gradient information in objective function res2 <- nloptr(x0=c(1.234, 5.678), eval_f=eval_f1, lb = c(-Inf, 0), ub = c(Inf, Inf), eval_g_ineq = eval_g1, opts = list("algorithm"="NLOPT_LD_MMA", "check_derivatives" = TRUE), a = a, b = b) print(res2)
This function returns a data.frame with all the options that can be supplied
to nloptr. The data.frame contains the default
values of the options and an explanation. A user-friendly way to show these
options is by using the function
nloptr.print.options.
nloptr.get.default.options()nloptr.get.default.options()
The return value contains a data.frame with the following
elements
name |
name of the option |
type |
type (numeric, logical, integer, character) |
possible_values |
string explaining the values the option can take |
default |
default value of the option (as a string) |
is_termination_condition |
is this option part of the termination conditions? |
description |
description of the option (taken from NLopt website if it's an option that is passed on to NLopt). |
Jelmer Ypma
This function prints a list of all the options that can be set when solving
a minimization problem using nloptr.
nloptr.print.options(opts.show = NULL, opts.user = NULL)nloptr.print.options(opts.show = NULL, opts.user = NULL)
opts.show |
list or vector with names of options. A description will be shown for the options in this list. By default, a description of all options is shown. |
opts.user |
object containing user supplied options. This argument is
optional. It is used when |
Jelmer Ypma
library('nloptr') nloptr.print.options() nloptr.print.options(opts.show = c("algorithm", "check_derivatives")) opts <- list("algorithm"="NLOPT_LD_LBFGS", "xtol_rel"=1.0e-8) nloptr.print.options(opts.user = opts)library('nloptr') nloptr.print.options() nloptr.print.options(opts.show = c("algorithm", "check_derivatives")) opts <- list("algorithm"="NLOPT_LD_LBFGS", "xtol_rel"=1.0e-8) nloptr.print.options(opts.user = opts)
This function prints the nloptr object that holds the results from a
minimization using nloptr.
## S3 method for class 'nloptr' print(x, show.controls = TRUE, ...)## S3 method for class 'nloptr' print(x, show.controls = TRUE, ...)
x |
object containing result from minimization. |
show.controls |
Logical or vector with indices. Should we show the value
of the control variables in the solution? If |
... |
further arguments passed to or from other methods. |
Jelmer Ypma
Subplex is a variant of Nelder-Mead that uses Nelder-Mead on a sequence of subspaces.
sbplx( x0, fn, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )sbplx( x0, fn, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
lower, upper
|
lower and upper bound constraints. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
SUBPLEX is claimed to be much more efficient and robust than the original Nelder-Mead while retaining the latter's facility with discontinuous objectives.
This implementation has explicit support for bound constraints via the
method in the Box paper as described on the neldermead help page.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
It is the request of Tom Rowan that reimplementations of his algorithm shall not use the name ‘subplex’.
T. Rowan, “Functional Stability Analysis of Numerical Algorithms”, Ph.D. thesis, Department of Computer Sciences, University of Texas at Austin, 1990.
subplex::subplex
# Fletcher and Powell's helic valley fphv <- function(x) 100*(x[3] - 10*atan2(x[2], x[1])/(2*pi))^2 + (sqrt(x[1]^2 + x[2]^2) - 1)^2 +x[3]^2 x0 <- c(-1, 0, 0) sbplx(x0, fphv) # 1 0 0 # Powell's Singular Function (PSF) psf <- function(x) (x[1] + 10*x[2])^2 + 5*(x[3] - x[4])^2 + (x[2] - 2*x[3])^4 + 10*(x[1] - x[4])^4 x0 <- c(3, -1, 0, 1) sbplx(x0, psf, control = list(maxeval = Inf, ftol_rel = 1e-6)) # 0 0 0 0 (?)# Fletcher and Powell's helic valley fphv <- function(x) 100*(x[3] - 10*atan2(x[2], x[1])/(2*pi))^2 + (sqrt(x[1]^2 + x[2]^2) - 1)^2 +x[3]^2 x0 <- c(-1, 0, 0) sbplx(x0, fphv) # 1 0 0 # Powell's Singular Function (PSF) psf <- function(x) (x[1] + 10*x[2])^2 + 5*(x[3] - x[4])^2 + (x[2] - 2*x[3])^4 + 10*(x[1] - x[4])^4 x0 <- c(3, -1, 0, 1) sbplx(x0, psf, control = list(maxeval = Inf, ftol_rel = 1e-6)) # 0 0 0 0 (?)
Sequential (least-squares) quadratic programming (SQP) algorithm for nonlinearly constrained, gradient-based optimization, supporting both equality and inequality constraints.
slsqp( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, heq = NULL, heqjac = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )slsqp( x0, fn, gr = NULL, lower = NULL, upper = NULL, hin = NULL, hinjac = NULL, heq = NULL, heqjac = NULL, nl.info = FALSE, control = list(), deprecatedBehavior = TRUE, ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
gradient of function |
lower, upper
|
lower and upper bound constraints. |
hin |
function defining the inequality constraints, that is
|
hinjac |
Jacobian of function |
heq |
function defining the equality constraints, that is |
heqjac |
Jacobian of function |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
deprecatedBehavior |
logical; if |
... |
additional arguments passed to the function. |
The algorithm optimizes successive second-order (quadratic/least-squares) approximations of the objective function (via BFGS updates), with first-order (affine) approximations of the constraints.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 1) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
See more infos at https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/.
Hans W. Borchers
Dieter Kraft, “A software package for sequential quadratic programming”, Technical Report DFVLR-FB 88-28, Institut fuer Dynamik der Flugsysteme, Oberpfaffenhofen, July 1988.
alabama::auglag, Rsolnp::solnp,
Rdonlp2::donlp2
## Solve the Hock-Schittkowski problem no. 100 with analytic gradients ## See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } S <- slsqp(x0.hs100, fn = fn.hs100, # no gradients and jacobians provided hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8, check_derivatives = TRUE), deprecatedBehavior = FALSE) ## The optimum value of the objective function should be 680.6300573 ## A suitable parameter vector is roughly ## (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) S## Solve the Hock-Schittkowski problem no. 100 with analytic gradients ## See https://apmonitor.com/wiki/uploads/Apps/hs100.apm x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) fn.hs100 <- function(x) {(x[1] - 10) ^ 2 + 5 * (x[2] - 12) ^ 2 + x[3] ^ 4 + 3 * (x[4] - 11) ^ 2 + 10 * x[5] ^ 6 + 7 * x[6] ^ 2 + x[7] ^ 4 - 4 * x[6] * x[7] - 10 * x[6] - 8 * x[7]} hin.hs100 <- function(x) {c( 2 * x[1] ^ 2 + 3 * x[2] ^ 4 + x[3] + 4 * x[4] ^ 2 + 5 * x[5] - 127, 7 * x[1] + 3 * x[2] + 10 * x[3] ^ 2 + x[4] - x[5] - 282, 23 * x[1] + x[2] ^ 2 + 6 * x[6] ^ 2 - 8 * x[7] - 196, 4 * x[1] ^ 2 + x[2] ^ 2 - 3 * x[1] * x[2] + 2 * x[3] ^ 2 + 5 * x[6] - 11 * x[7]) } S <- slsqp(x0.hs100, fn = fn.hs100, # no gradients and jacobians provided hin = hin.hs100, nl.info = TRUE, control = list(xtol_rel = 1e-8, check_derivatives = TRUE), deprecatedBehavior = FALSE) ## The optimum value of the objective function should be 680.6300573 ## A suitable parameter vector is roughly ## (2.330, 1.9514, -0.4775, 4.3657, -0.6245, 1.0381, 1.5942) S
StoGO is a global optimization algorithm that works by systematically dividing the search space—which must be bound-constrained—into smaller hyper-rectangles via a branch-and-bound technique, and searching them using a gradient-based local-search algorithm (a BFGS variant), optionally including some randomness.
stogo( x0, fn, gr = NULL, lower = NULL, upper = NULL, maxeval = 10000, xtol_rel = 1e-06, randomized = FALSE, nl.info = FALSE, ... )stogo( x0, fn, gr = NULL, lower = NULL, upper = NULL, maxeval = 10000, xtol_rel = 1e-06, randomized = FALSE, nl.info = FALSE, ... )
x0 |
initial point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
optional gradient of the objective function. |
lower, upper
|
lower and upper bound constraints. |
maxeval |
maximum number of function evaluations. |
xtol_rel |
stopping criterion for relative change reached. |
randomized |
logical; shall a randomizing variant be used? |
nl.info |
logical; shall the original NLopt info be shown. |
... |
additional arguments passed to the function. |
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
Only bounds-constrained problems are supported by this algorithm.
Hans W. Borchers
S. Zertchaninov and K. Madsen, “A C++ Programme for Global Optimization,” IMM-REP-1998-04, Department of Mathematical Modelling, Technical University of Denmark.
## Rosenbrock Banana objective function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} x0 <- c(-1.2, 1) lb <- c(-3, -3) ub <- c(3, 3) ## The function as written above has a minimum of 0 at (1, 1) stogo(x0 = x0, fn = rbf, lower = lb, upper = ub)## Rosenbrock Banana objective function rbf <- function(x) {(1 - x[1]) ^ 2 + 100 * (x[2] - x[1] ^ 2) ^ 2} x0 <- c(-1.2, 1) lb <- c(-3, -3) ub <- c(3, 3) ## The function as written above has a minimum of 0 at (1, 1) stogo(x0 = x0, fn = rbf, lower = lb, upper = ub)
Truncated Newton methods, also called Newton-iterative methods, solve an approximating Newton system using a conjugate-gradient approach and are related to limited-memory BFGS.
tnewton( x0, fn, gr = NULL, lower = NULL, upper = NULL, precond = TRUE, restart = TRUE, nl.info = FALSE, control = list(), ... )tnewton( x0, fn, gr = NULL, lower = NULL, upper = NULL, precond = TRUE, restart = TRUE, nl.info = FALSE, control = list(), ... )
x0 |
starting point for searching the optimum. |
fn |
objective function that is to be minimized. |
gr |
gradient of function |
lower, upper
|
lower and upper bound constraints. |
precond |
logical; preset L-BFGS with steepest descent. |
restart |
logical; restarting L-BFGS with steepest descent. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of options, see |
... |
additional arguments passed to the function. |
Truncated Newton methods are based on approximating the objective with a quadratic function and applying an iterative scheme such as the linear conjugate-gradient algorithm.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 1) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
Less reliable than Newton's method, but can handle very large problems.
Hans W. Borchers
R. S. Dembo and T. Steihaug, “Truncated Newton algorithms for large-scale optimization,” Math. Programming 26, p. 190-212 (1982).
flb <- function(x) { p <- length(x) sum(c(1, rep(4, p - 1)) * (x - c(1, x[-p]) ^ 2) ^ 2) } # 25-dimensional box constrained: par[24] is *not* at boundary S <- tnewton(rep(3, 25L), flb, lower = rep(2, 25L), upper = rep(4, 25L), nl.info = TRUE, control = list(xtol_rel = 1e-8)) ## Optimal value of objective function: 368.105912874334 ## Optimal value of controls: 2 ... 2 2.109093 4flb <- function(x) { p <- length(x) sum(c(1, rep(4, p - 1)) * (x - c(1, x[-p]) ^ 2) ^ 2) } # 25-dimensional box constrained: par[24] is *not* at boundary S <- tnewton(rep(3, 25L), flb, lower = rep(2, 25L), upper = rep(4, 25L), nl.info = TRUE, control = list(xtol_rel = 1e-8)) ## Optimal value of objective function: 368.105912874334 ## Optimal value of controls: 2 ... 2 2.109093 4
Shifted limited-memory variable-metric algorithm.
varmetric( x0, fn, gr = NULL, rank2 = TRUE, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )varmetric( x0, fn, gr = NULL, rank2 = TRUE, lower = NULL, upper = NULL, nl.info = FALSE, control = list(), ... )
x0 |
initial point for searching the optimum. |
fn |
objective function to be minimized. |
gr |
gradient of function |
rank2 |
logical; if true uses a rank-2 update method, else rank-1. |
lower, upper
|
lower and upper bound constraints. |
nl.info |
logical; shall the original NLopt info been shown. |
control |
list of control parameters, see |
... |
further arguments to be passed to the function. |
Variable-metric methods are a variant of the quasi-Newton methods, especially adapted to large-scale unconstrained (or bound constrained) minimization.
List with components:
par |
the optimal solution found so far. |
value |
the function value corresponding to |
iter |
number of (outer) iterations, see |
convergence |
integer code indicating successful completion (> 0) or a possible error number (< 0). |
message |
character string produced by NLopt and giving additional information. |
Based on L. Luksan's Fortran implementation of a shifted limited-memory variable-metric algorithm.
Hans W. Borchers
J. Vlcek and L. Luksan, “Shifted limited-memory variable metric methods for large-scale unconstrained minimization,” J. Computational Appl. Math. 186, p. 365-390 (2006).
flb <- function(x) { p <- length(x) sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) } # 25-dimensional box constrained: par[24] is *not* at the boundary S <- varmetric(rep(3, 25), flb, lower=rep(2, 25), upper=rep(4, 25), nl.info = TRUE, control = list(xtol_rel=1e-8)) ## Optimal value of objective function: 368.105912874334 ## Optimal value of controls: 2 ... 2 2.109093 4flb <- function(x) { p <- length(x) sum(c(1, rep(4, p-1)) * (x - c(1, x[-p])^2)^2) } # 25-dimensional box constrained: par[24] is *not* at the boundary S <- varmetric(rep(3, 25), flb, lower=rep(2, 25), upper=rep(4, 25), nl.info = TRUE, control = list(xtol_rel=1e-8)) ## Optimal value of objective function: 368.105912874334 ## Optimal value of controls: 2 ... 2 2.109093 4