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