1# Copyright (C) 2014 Hans W. Borchers. All Rights Reserved. 2# This code is published under the L-GPL. 3# 4# File: cobyla.R 5# Author: Hans W. Borchers 6# Date: 27 January 2014 7# 8# Wrapper to solve optimization problem using COBYLA, BOBYQA, and NEWUOA. 9 10 11 12#' Constrained Optimization by Linear Approximations 13#' 14#' COBYLA is an algorithm for derivative-free optimization with nonlinear 15#' inequality and equality constraints (but see below). 16#' 17#' It constructs successive linear approximations of the objective function and 18#' constraints via a simplex of n+1 points (in n dimensions), and optimizes 19#' these approximations in a trust region at each step. 20#' 21#' COBYLA supports equality constraints by transforming them into two 22#' inequality constraints. As this does not give full satisfaction with the 23#' implementation in NLOPT, it has not been made available here. 24#' 25#' @param x0 starting point for searching the optimum. 26#' @param fn objective function that is to be minimized. 27#' @param lower,upper lower and upper bound constraints. 28#' @param hin function defining the inequality constraints, that is 29#' \code{hin>=0} for all components. 30#' @param nl.info logical; shall the original NLopt info been shown. 31#' @param control list of options, see \code{nl.opts} for help. 32#' @param ... additional arguments passed to the function. 33#' 34#' @return List with components: 35#' \item{par}{the optimal solution found so far.} 36#' \item{value}{the function value corresponding to \code{par}.} 37#' \item{iter}{number of (outer) iterations, see \code{maxeval}.} 38#' \item{convergence}{integer code indicating successful completion (> 0) 39#' or a possible error number (< 0).} 40#' \item{message}{character string produced by NLopt and giving additional 41#' information.} 42#' 43#' @author Hans W. Borchers 44#' 45#' @note The original code, written in Fortran by Powell, was converted in C 46#' for the Scipy project. 47#' 48#' @seealso \code{\link{bobyqa}}, \code{\link{newuoa}} 49#' 50#' @references M. J. D. Powell, ``A direct search optimization method that 51#' models the objective and constraint functions by linear interpolation,'' in 52#' Advances in Optimization and Numerical Analysis, eds. S. Gomez and J.-P. 53#' Hennart (Kluwer Academic: Dordrecht, 1994), p. 51-67. 54#' 55#' @examples 56#' 57#' ### Solve Hock-Schittkowski no. 100 58#' x0.hs100 <- c(1, 2, 0, 4, 0, 1, 1) 59#' fn.hs100 <- function(x) { 60#' (x[1]-10)^2 + 5*(x[2]-12)^2 + x[3]^4 + 3*(x[4]-11)^2 + 10*x[5]^6 + 61#' 7*x[6]^2 + x[7]^4 - 4*x[6]*x[7] - 10*x[6] - 8*x[7] 62#' } 63#' hin.hs100 <- function(x) { 64#' h <- numeric(4) 65#' h[1] <- 127 - 2*x[1]^2 - 3*x[2]^4 - x[3] - 4*x[4]^2 - 5*x[5] 66#' h[2] <- 282 - 7*x[1] - 3*x[2] - 10*x[3]^2 - x[4] + x[5] 67#' h[3] <- 196 - 23*x[1] - x[2]^2 - 6*x[6]^2 + 8*x[7] 68#' h[4] <- -4*x[1]^2 - x[2]^2 + 3*x[1]*x[2] -2*x[3]^2 - 5*x[6] +11*x[7] 69#' return(h) 70#' } 71#' 72#' S <- cobyla(x0.hs100, fn.hs100, hin = hin.hs100, 73#' nl.info = TRUE, control = list(xtol_rel = 1e-8, maxeval = 2000)) 74#' ## Optimal value of objective function: 680.630057374431 75#' 76#' @export cobyla 77cobyla <- 78function(x0, fn, lower = NULL, upper = NULL, hin = NULL, 79 nl.info = FALSE, control = list(), ...) 80{ 81 opts <- nl.opts(control) 82 opts["algorithm"] <- "NLOPT_LN_COBYLA" 83 84 f1 <- match.fun(fn) 85 fn <- function(x) f1(x, ...) 86 87 if (!is.null(hin)) { 88 if ( getOption('nloptr.show.inequality.warning') ) { 89 message('For consistency with the rest of the package the inequality sign may be switched from >= to <= in a future nloptr version.') 90 } 91 92 f2 <- match.fun(hin) 93 hin <- function(x) (-1)*f2(x, ...) # NLOPT expects hin <= 0 94 } 95 96 S0 <- nloptr(x0, 97 eval_f = fn, 98 lb = lower, 99 ub = upper, 100 eval_g_ineq = hin, 101 opts = opts) 102 103 if (nl.info) print(S0) 104 S1 <- list(par = S0$solution, value = S0$objective, iter = S0$iterations, 105 convergence = S0$status, message = S0$message) 106 return(S1) 107} 108 109 110 111 112#' Bound Optimization by Quadratic Approximation 113#' 114#' BOBYQA performs derivative-free bound-constrained optimization using an 115#' iteratively constructed quadratic approximation for the objective function. 116#' 117#' This is an algorithm derived from the BOBYQA Fortran subroutine of Powell, 118#' converted to C and modified for the NLOPT stopping criteria. 119#' 120#' @param x0 starting point for searching the optimum. 121#' @param fn objective function that is to be minimized. 122#' @param lower,upper lower and upper bound constraints. 123#' @param nl.info logical; shall the original NLopt info been shown. 124#' @param control list of options, see \code{nl.opts} for help. 125#' @param ... additional arguments passed to the function. 126#' 127#' @return List with components: 128#' \item{par}{the optimal solution found so far.} 129#' \item{value}{the function value corresponding to \code{par}.} 130#' \item{iter}{number of (outer) iterations, see \code{maxeval}.} 131#' \item{convergence}{integer code indicating successful completion (> 0) 132#' or a possible error number (< 0).} 133#' \item{message}{character string produced by NLopt and giving additional 134#' information.} 135#' 136#' @export bobyqa 137#' 138#' @note Because BOBYQA constructs a quadratic approximation of the objective, 139#' it may perform poorly for objective functions that are not 140#' twice-differentiable. 141#' 142#' @seealso \code{\link{cobyla}}, \code{\link{newuoa}} 143#' 144#' @references M. J. D. Powell. ``The BOBYQA algorithm for bound constrained 145#' optimization without derivatives,'' Department of Applied Mathematics and 146#' Theoretical Physics, Cambridge England, technical reportNA2009/06 (2009). 147#' 148#' @examples 149#' 150#' fr <- function(x) { ## Rosenbrock Banana function 151#' 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2 152#' } 153#' (S <- bobyqa(c(0, 0, 0), fr, lower = c(0, 0, 0), upper = c(0.5, 0.5, 0.5))) 154#' 155bobyqa <- 156function(x0, fn, lower = NULL, upper = NULL, 157 nl.info = FALSE, control = list(), ...) 158{ 159 opts <- nl.opts(control) 160 opts["algorithm"] <- "NLOPT_LN_BOBYQA" 161 162 fun <- match.fun(fn) 163 fn <- function(x) fun(x, ...) 164 165 S0 <- nloptr(x0, fn, lb = lower, ub = upper, 166 opts = opts) 167 168 if (nl.info) print(S0) 169 S1 <- list(par = S0$solution, value = S0$objective, iter = S0$iterations, 170 convergence = S0$status, message = S0$message) 171 return(S1) 172} 173 174 175 176 177#' New Unconstrained Optimization with quadratic Approximation 178#' 179#' NEWUOA solves quadratic subproblems in a spherical trust regionvia a 180#' truncated conjugate-gradient algorithm. For bound-constrained problems, 181#' BOBYQA shold be used instead, as Powell developed it as an enhancement 182#' thereof for bound constraints. 183#' 184#' This is an algorithm derived from the NEWUOA Fortran subroutine of Powell, 185#' converted to C and modified for the NLOPT stopping criteria. 186#' 187#' @param x0 starting point for searching the optimum. 188#' @param fn objective function that is to be minimized. 189#' @param nl.info logical; shall the original NLopt info been shown. 190#' @param control list of options, see \code{nl.opts} for help. 191#' @param ... additional arguments passed to the function. 192#' 193#' @return List with components: 194#' \item{par}{the optimal solution found so far.} 195#' \item{value}{the function value corresponding to \code{par}.} 196#' \item{iter}{number of (outer) iterations, see \code{maxeval}.} 197#' \item{convergence}{integer code indicating successful completion (> 0) 198#' or a possible error number (< 0).} 199#' \item{message}{character string produced by NLopt and giving additional 200#' information.} 201#' 202#' @export newuoa 203#' 204#' @author Hans W. Borchers 205#' 206#' @note NEWUOA may be largely superseded by BOBYQA. 207#' 208#' @seealso \code{\link{bobyqa}}, \code{\link{cobyla}} 209#' 210#' @references M. J. D. Powell. ``The BOBYQA algorithm for bound constrained 211#' optimization without derivatives,'' Department of Applied Mathematics and 212#' Theoretical Physics, Cambridge England, technical reportNA2009/06 (2009). 213#' 214#' @examples 215#' 216#' fr <- function(x) { ## Rosenbrock Banana function 217#' 100 * (x[2] - x[1]^2)^2 + (1 - x[1])^2 218#' } 219#' (S <- newuoa(c(1, 2), fr)) 220#' 221newuoa <- 222function(x0, fn, nl.info = FALSE, control = list(), ...) 223{ 224 opts <- nl.opts(control) 225 opts["algorithm"] <- "NLOPT_LN_NEWUOA" 226 227 fun <- match.fun(fn) 228 fn <- function(x) fun(x, ...) 229 230 S0 <- nloptr(x0, fn, opts = opts) 231 232 if (nl.info) print(S0) 233 S1 <- list(par = S0$solution, value = S0$objective, iter = S0$iterations, 234 convergence = S0$status, message = S0$message) 235 return(S1) 236} 237