# source:trunk/Ipopt/contrib/RInterface/inst/doc/ipoptr.Rnw@1861

Last change on this file since 1861 was 1861, checked in by andreasw, 3 years ago

moved Ipopt trunk from Common Public License to successor Eclispe Public License (see e.g. http://www.ibm.com/developerworks/library/os-cplfaq.html)

File size: 26.0 KB
Line
1\documentclass[a4paper]{article}
2\usepackage[english]{babel}
3\usepackage{apacite}
4\usepackage{graphicx}
5
6% \VignetteIndexEntry{Introduction to ipoptr: an R interface to Ipopt}
7% \VignetteKeyword{optimize}
8% \VignetteKeyword{interface}
9
10\SweaveOpts{keep.source=TRUE}
11\SweaveOpts{prefix.string = figs/plot, eps = FALSE, pdf = TRUE, tikz = FALSE}
12
13%\pgfrealjobname{ipoptr}
14
15\title{Introduction to \texttt{ipoptr}: an R interface to Ipopt
16\thanks{This package should be considered in beta and comments about any aspect of the package are welcome. Thanks to Alexios Ghalanos for comments. This document is an R vignette prepared with the aid of \texttt{Sweave}, Leisch(2002). Financial support of the UK Economic and Social Research Council through a grant (RES-589-28-0001) to the ESRC Centre for Microdata Methods and Practice (CeMMAP) is gratefully acknowledged.}}
17\author{Jelmer Ypma}
18\begin{document}
19\maketitle
20\nocite{Leisch2002}
21
22\DefineVerbatimEnvironment{Sinput}{Verbatim}{xleftmargin=2em}
23\DefineVerbatimEnvironment{Soutput}{Verbatim}{xleftmargin=2em}
24\DefineVerbatimEnvironment{Scode}{Verbatim}{xleftmargin=2em}
25\fvset{listparameters={\setlength{\topsep}{0pt}}}
26\renewenvironment{Schunk}{\vspace{\topsep}}{\vspace{\topsep}}
27
28
29<<setSweaveOptions,echo=FALSE>>=
30# have an (invisible) initialization noweb chunk
31# to remove the default continuation prompt '>'
32options(continue = " ")
33options(width = 60)
34
35# eliminate margin space above plots
36options(SweaveHooks=list(fig=function()
37    par(mar=c(5.1, 4.1, 1.1, 2.1))))
38@
39
40\begin{abstract}
41This document describes how to use \texttt{ipoptr}, which is an R interface to Ipopt (Interior Point Optimizer). Ipopt is an open source software package for large-scale nonlinear optimization \cite{WachterBiegler2006}. It can be used to solve general nonlinear programming problems with nonlinear constraints and lower and upper bounds for the controls. Ipopt is written in C++ and is released as open source code under the Eclipse Public License (EPL). It is available from the COIN-OR initiative. The code has been written by Carl Laird and Andreas W\"achter, who is the COIN project leader for Ipopt.
42\end{abstract}
43
44\section{Introduction}
45Ipopt is designed to find (local) solutions of mathematical optimization problems of the from
46\begin{eqnarray*}
47&&\min_{x \in R^n} f(x) \\
48&s.t.& g_L <= g(x) <= g_U \\
49&&     x_L <=  x   <= x_U
50\end{eqnarray*}
51where $f(x): R^n \rightarrow R$ is the objective function, and $g(x): R^n \rightarrow R^m$ are the constraint functions. The vectors $g_L$ and $g_U$ denote the lower and upper bounds on the constraints, and the vectors $x_L$ and $x_U$ are the bounds on the variables $x$. The functions $f(x)$ and $g(x)$ can be nonlinear and nonconvex, but should be twice continuously differentiable. Note that equality constraints can be formulated in the above formulation by setting the corresponding components of $g_L$ and $g_U$ to the same value.
52
53This vignette describes how to formulate minimization problems to be solved with the R interface to Ipopt. If you want to use the C++ interface directly or are interested in the Matlab interface, there are other sources of documentation avialable. Some of the information here is heavily based on the Ipopt Wiki\footnote{\texttt{https://projects.coin-or.org/Ipopt}} and generally that is a good source to find additional information, for instance on which options to use. All credit for implementing the C++ code for Ipopt should go to Andreas W\"achter and Carl Laird. Please show your appreciation by citing their paper.
54
55\section{Installation}
56Installing the \texttt{ipoptr} package is not as straightforward as most other R packages, because it depends on Ipopt. To install (and compile) Ipopt and the R interface a C/C++ compiler has to be available. On Windows I was succesful using MSYS to compile Ipopt and then use Rtools\footnote{\texttt{http://www.murdoch-sutherland.com/Rtools/}} to compile the R interface from source. On Ubuntu no additional tools were needed.
57
58Detailed installation instructions for Ipopt are available on \texttt{http://www.coin-or.org/Ipopt/documentation}. You should follow these first, before trying to install the R interface. Ipopt needs to be configured using the \texttt{-fPIC} flag for all GNU compilers. For 64bit Linux one needs to specify
59\begin{verbatim}
61\end{verbatim}
62
63For the installation of the R interface, I will assume that you have a working installation of Ipopt (i.e. \texttt{configure}, \texttt{make} and \texttt{make install} executed without problems).
64
65During the installation of Ipopt a file \texttt{Makevars} has been created in the source directory of the R interface, e.g. \texttt{\$IPOPTDIR/build/Ipopt/contrib/RInterface/src} if you used the same build directory as in the Ipopt installation notes, \texttt{\$IPOPTDIR/build}. The file \texttt{Makevars} in this directory has been configured for your system.
66
67To install the R interface, this file has to be copied to the Ipopt directory containing the source code, \texttt{\$IPOPTDIR/Ipopt/contrib/RInterface/src}. Notice that the path of this directory is different from the directory where you built Ipopt (\texttt{build} is not there). 68 69The source directory \texttt{\$IPOPTDIR/Ipopt/contrib/RInterface/src} contains four files, \texttt{ipoptr.cpp}, \texttt{IpoptRNLP.cpp}, \texttt{IpoptRNLP.hpp} and \texttt{Makevars.in}. Copy \texttt{Makevars} in this directory and, if you're on Windows, rename it to \texttt{Makevars.win}.
70
71You can then install the package from R with the command
72<<installIpoptPackage, eval=FALSE>>=
73install.packages('$IPOPTDIR/Ipopt/contrib/RInterface', repos=NULL, type='source') 74@ 75where the first argument specifies the directory where the source code for the R interface to Ipopt is located. You should now be able to load the R interface to Ipopt and read the help. 76<<testIpoptInstallation, eval=FALSE>>= 77library('ipoptr') 78?ipoptr 79@ 80 81\section{Minimizing the Rosenbrock Banana function} 82As a first example we will solve an unconstrained minimization problem. The function we look at is the Rosenbrock Banana function 83$84f( x ) = 100 \left( x_2 - x_1^2 \right)^2 + \left(1 - x_1 \right)^2, 85$ 86which is also used as an example in the documentation for the standard R optimizer \texttt{optim}. The gradient of the objective function is given by 87$88\nabla f( x ) = 89\left( \begin{array}[1]{c} 90-400 \cdot x_1 \cdot (x_2 - x_1^2) - 2 \cdot (1 - x_1) \\ 91 200 \cdot (x_2 - x_1^2) 92\end{array} \right). 93$ 94Ipopt always needs gradients to be supplied by the user. After loading the library 95<<loadLibrary>>= 96library(ipoptr) 97@ 98we start by specifying the objective function and its gradient 99<<defineRosenbrockBanana>>= 100## Rosenbrock Banana function 101eval_f <- function(x) { 102 return( 100 * (x[2] - x[1] * x[1])^2 + (1 - x[1])^2 ) 103} 104 105## Gradient of Rosenbrock Banana function 106eval_grad_f <- function(x) { 107 return( c( -400 * x[1] * (x[2] - x[1] * x[1]) - 2 * (1 - x[1]), 108 200 * (x[2] - x[1] * x[1]) ) ) 109} 110@ 111We define initial values 112<<setRosenbrockBananaInitialValues>>= 113# initial values 114x0 <- c( -1.2, 1 ) 115@ 116and then minimize the function using the \texttt{ipoptr} command. This command runs some checks on the supplied inputs and returns an object with the exit code of the solver, the optimal value of the objective function and the solution. The checks do not always return very informative messages, but usually there is something wrong with dimensions (e.g. \texttt{eval\_grad\_f} returns a vector that doesn't have the same size as \texttt{x0}). 117<<solveRosenbrockBanana>>= 118# solve Rosenbrock Banana function 119res <- ipoptr( x0=x0, 120 eval_f=eval_f, 121 eval_grad_f=eval_grad_f ) 122@ 123These are the minimal arguments that have to be supplied. If, like above, no Hessian is defined, Ipopt uses an approximation. We can see the results by printing the resulting object. 124<<printRosenbrockBanana>>= 125print( res ) 126@ 127It's advised to always check the exit code for convergence of the problem and in this case we can see that the algorithm terminated successfully. Ipopt used 47 iterations to find the solution and the optimal value of the objective function and the controls are given as well. 128 129If you do not want to, or cannot calculate the gradient analytically, you can supply a function \texttt{eval\_grad\_f} that approximates the gradient. However, this is not advisable and might result in convergence problems, for instance by not finding the minimum, or by finding the wrong minimum. We can see this from the following example where we approximate \texttt{eval\_grad\_f} using finite differences 130<<defineApproxGradF>>= 131# Approximate eval_f using finite differences 132# http://en.wikipedia.org/wiki/Numerical_differentiation 133approx_grad_f <- function( x ) { 134 minAbsValue <- 0 135 stepSize <- sqrt( .Machine$double.eps )
136
137    # if we evaluate at 0, we need a different step size
138    stepSizeVec     <- ifelse( abs(x) <= minAbsValue,
139                               stepSize^3,
140                               x * stepSize )
141
142    x_prime <- x
143    f       <- eval_f( x )
144    grad_f  <- rep( 0, length(x) )
145    for (i in 1:length(x)) {
146        x_prime[i]      <- x[i] + stepSizeVec[i]
147        stepSizeVec[i]  <- x_prime[i] - x[i]
148
149        f_prime         <- eval_f( x_prime )
150        grad_f[i]       <- ( f_prime - f )/stepSizeVec[i]
151        x_prime[i]      <- x[i]
152    }
153
155}
156@
157and using this approximation to minimize the same Rosenbrock Banana function.
158<<solveRosenbrockBananaApproximation1>>=
159# increase the maximum number of iterations
160opts <- list("tol"=1.0e-8, "max_iter"=5000)
161
162# solve Rosenbrock Banana function with approximated gradient
163print( ipoptr( x0=x0,
164               eval_f=eval_f,
166               opts=opts) )
167@
168In this case 5000 iterations are not enough to solve the minimization problem to the required tolerance. This has to do with the step size we choose to approximate the gradient
169<<showMachineEpsilon>>=
170sqrt( .Machine$double.eps ) 171@ 172which is of the same order of magnitude. If we decrease the tolerance, the algorithm converges, but the solution is less precise than if we supply gradients and it takes more iterations to get there. 173<<solveRosenbrockBananaApproximation2>>= 174# decrease the convergence criterium 175opts <- list("tol"=1.0e-7) 176 177# solve Rosenbrock Banana function with approximated gradient 178print( ipoptr( x0=x0, 179 eval_f=eval_f, 180 eval_grad_f=approx_grad_f, 181 opts=opts) ) 182@ 183 184\section{Sparse matrix structure} 185Ipopt can handle sparseness in the Jacobian of the constraints and the Hessian. The sparseness structure should be defined in advance and stay the same throughout the minimization procedure. A sparseness structure can be defined as a list of vectors, where each vector contains the indices of the non-zero elements of one row. E.g. the matrix 186$187\left( \begin{array}[4]{cccc} 188. & . & . & 1 \\ 1891 & 1 & . & . \\ 1901 & 1 & 1 & 1 191\end{array} \right) 192$ 193has a non-zero element in position 4 in the first row. In the second row it has non-zero elements in position 1 and 2, and the third row contains non-zero elements at every position. Its structure can be defined as 194<<sparseEx1>>= 195sparse_structure <- list( c( 4 ), c( 1, 2 ), c( 1, 2, 3, 4 ) ) 196@ 197The function \texttt{make.sparse} can simplify this procedure 198<<sparseEx2>>= 199make.sparse( rbind( c(0, 0, 0, 1), c( 1, 1, 0, 0 ), c( 1, 1, 1, 1 ) ) ) 200@ 201This function takes a matrix as argument. All non-zero elements in this matrix will be defined as non-zero in the sparseness structure, \texttt{NA} or \texttt{NaN} are not allowed. The function \texttt{print.sparseness} shows the non-zero elements 202<<sparseEx3>>= 203print.sparseness( sparse_structure ) 204@ 205By default \texttt{print.sparseness} shows the indices of the non-zero elements in the sparse matrix. Values for the non-zero elements of a sparse matrix have to be supplied in one vector, in the same order as the the non-zero elements occur in the structure. I.e. the order of the indices matters and the values of the following two matrices should be supplied in a different order 206<<sparseEx4>>= 207print.sparseness( list( c(1,3,6,8), c(2,5), c(3,7,9) ) ) 208print.sparseness( list( c(3,1,6,8), c(2,5), c(3,9,7) ) ) 209@ 210Since the sparseness structure defines the indices of non-zero elements by row, the order of the rows cannot be changed in the R implementation. In principle a more general order of the non-zero elements (independent of row or column) could be specified, which can be added as a feature on request. Below are two final examples on sparseness structure (see \texttt{?print.sparseness} for more options and examples) 211<<sparseEx5>>= 212# print lower-diagonal 5x5 matrix generated with make.sparse 213A_lower <- make.sparse( lower.tri( matrix(1, nrow=5, ncol=5), diag=TRUE ) ) 214print.sparseness( A_lower ) 215 216# print a diagonal 5x5 matrix without indices counts 217A_diag <- make.sparse( diag(5) > 0 ) 218print.sparseness( A_diag, indices=FALSE ) 219@ 220 221For larger matrices it is easier to plot them using the \texttt{plot.sparseness} command 222<<label=sparsefig,fig=TRUE,echo=TRUE,include=FALSE,results=hide,width=6,height=2.35>>= 223s <- do.call( "cbind", lapply( 1:5, function(i) { 224 diag(5) %x% matrix(1, nrow=5, ncol=20) 225 } ) ) 226s <- do.call( "rbind", lapply( 1:5, function(i) { s } ) ) 227s <- cbind( matrix( 1, nrow=nrow(s), ncol=40 ), s ) 228plot.sparseness( make.sparse( s ) ) 229@ 230The resulting sparse matrix structure from this code can be seen in figure \ref{fig:sparse}. All non-zero elements are shown as black dots by default. 231 232\begin{figure}[htbp] 233\begin{center} 234\includegraphics[width=12.0cm]{figs/plot-sparsefig.pdf} 235\caption{Plot of large sparseness structure} 236\label{fig:sparse} 237\end{center} 238\end{figure} 239 240\section{Supplying the Hessian} 241Now that we know how to define a sparseness structure we can supply the Hessian to the Rosenbrock Banana function from above. Its Hessian is given by 242$243\nabla^2 f( x ) = \left( \begin{array}[2]{rr} 2442 - 400 \cdot (x_2 - x_1^2) + 800 x_1^2 & -400 x_1 \\ 245-400 x_1 & 200 246\end{array} \right) 247$ 248Ipopt needs the Hessian of the Lagrangian in the following form 249$250\sigma_f \nabla^2 f(x) + \sum_{i=1}^m \lambda_i \nabla^2 g_i(x), 251$ 252where$g_i(x)$represents the$i$th of$m$constraints,$\lambda_i$are the multipliers of the constraints and$\sigma_f$is introduced so that Ipopt can ask for the Hessian of the objective or the constraints independently if required. 253 254In this case we don't have any constraints. The user-defined function \texttt{eval\_h} to define the Hessian takes three arguments. The first argument contains the value of the control variables,$x$, the second argument contains the multiplication factor of the Hessian of the objective function,$\sigma_f$, and the third argument contains a vector with the multipliers of the constraints,$\lambda$. We can define the structure of the Hessian and the function to evaluate the Hessian as follows 255<<defineRosenbrockBananaHessian>>= 256# The Hessian for this problem is actually dense, 257# This is a symmetric matrix, fill the lower left triangle only. 258eval_h_structure <- list( c(1), c(1,2) ) 259 260eval_h <- function( x, obj_factor, hessian_lambda ) { 261 return( obj_factor*c( 2 - 400*(x[2] - x[1]^2) + 800*x[1]^2, # 1,1 262 -400*x[1], # 2,1 263 200 ) ) # 2,2 264} 265@ 266Note that we only specify the lower half of the Hessian, since it is a symmetric matrix. Also, \texttt{eval\_h} returns a vector with all the non-zero elements of the Hessian in the same order as the non-zero indices in the sparseness structure. Then we minimize the function using the \texttt{ipoptr} command 267<<solveRosenbrockBananaHessian>>= 268opts <- list("print_level"=0, 269 "file_print_level"=12, 270 "output_file"="banana.out", 271 "tol"=1.0e-8) 272 273# solve Rosenbrock Banana function with analytic Hessian 274print( ipoptr( x0=x0, 275 eval_f=eval_f, 276 eval_grad_f=eval_grad_f, 277 eval_h=eval_h, 278 eval_h_structure=eval_h_structure, 279 opts=opts) ) 280@ 281Here we also supplied options to not print any intermediate information to the R screen (\texttt{print\_level=0}). Printing output to the screen directly from Ipopt does not work in all R terminals correctly, so it might be that even though you specify a positive number here, there will still be no output visible on the screen. If you want to print things to the screen, a workaround is to do this directly in the R functions you defined, such as \texttt{eval\_f}. 282 283Also, to inspect more details about the minimization we can write all the output to a file, which will be created in the current working directory. For larger problems, having a large number for \texttt{file\_print\_level} can easily generate very large files, which is probably not desirable. Many more options are available, and a full list of all the options can be found at the Ipopt website, \texttt{http://www.coin-or.org/Ipopt/documentation/node59.html\#app.options\_ref}. Options can also be supplied from an option file, which can be specified in \texttt{option\_file\_name}. 284 285\section{Adding constraints} 286To look at how we can add constraints to a problem, we take example problem number 71 from the Hock-Schittkowsky test suite, which is also used in the Ipopt C++ tutorial. The problem is 287\begin{eqnarray*} 288&&\min_{x} x_1 \cdot x_4 \cdot (x_1 + x_2 + x_3) + x_3 \\ 289&s.t.& \\ 290&& x_1 \cdot x_2 \cdot x_3 \cdot x_4 >= 25 \\ 291&& x_1^2 + x_2^2 + x_3^2 + x_4^2 = 40 \\ 292&& 1 <= x_1,x_2,x_3,x_4 <= 5, 293\end{eqnarray*} 294and we use$x = (1, 5, 5, 1)$as initial values. In this problem we have one inequality constraint, one equality constraint and upper and lower bounds for all the variables. The optimal solution is$(1.00000000, 4.74299963, 3.82114998, 1.37940829)$. First we define the objective function and its gradient 295<<defineHS071>>= 296eval_f <- function( x ) { 297 return( x[1]*x[4]*(x[1] + x[2] + x[3]) + x[3] ) 298} 299 300eval_grad_f <- function( x ) { 301 return( c( x[1] * x[4] + x[4] * (x[1] + x[2] + x[3]), 302 x[1] * x[4], 303 x[1] * x[4] + 1.0, 304 x[1] * (x[1] + x[2] + x[3]) ) ) 305} 306@ 307Then we define a function that returns the value of the two constraints. We define the bounds of the constraints (in this case the$g_L$and$g_U$are$25$and$40$) later. 308<<defineHS071constraints>>= 309# constraint functions 310eval_g <- function( x ) { 311 return( c( x[1] * x[2] * x[3] * x[4], 312 x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2 ) ) 313} 314@ 315Then we define the structure of the Jacobian, which is a dense matrix in this case, and function to evaluate it 316<<defineHS071jacobian>>= 317eval_jac_g_structure <- list( c(1,2,3,4), c(1,2,3,4) ) 318 319eval_jac_g <- function( x ) { 320 return( c ( x[2]*x[3]*x[4], 321 x[1]*x[3]*x[4], 322 x[1]*x[2]*x[4], 323 x[1]*x[2]*x[3], 324 2.0*x[1], 325 2.0*x[2], 326 2.0*x[3], 327 2.0*x[4] ) ) 328} 329@ 330The Hessian is also dense, but it looks slightly more complicated because we have to take into account the Hessian of the objective function and of the constraints at the same time, although you could write a function to calculate them both separately and then return the combined result in \texttt{eval\_h}. 331<<defineHS071hessian>>= 332# The Hessian for this problem is actually dense, 333# This is a symmetric matrix, fill the lower left triangle only. 334eval_h_structure <- list( c(1), c(1,2), c(1,2,3), c(1,2,3,4) ) 335 336eval_h <- function( x, obj_factor, hessian_lambda ) { 337 338 values <- numeric(10) 339 values[1] = obj_factor * (2*x[4]) # 1,1 340 341 values[2] = obj_factor * (x[4]) # 2,1 342 values[3] = 0 # 2,2 343 344 values[4] = obj_factor * (x[4]) # 3,1 345 values[5] = 0 # 4,2 346 values[6] = 0 # 3,3 347 348 values[7] = obj_factor * (2*x[1] + x[2] + x[3]) # 4,1 349 values[8] = obj_factor * (x[1]) # 4,2 350 values[9] = obj_factor * (x[1]) # 4,3 351 values[10] = 0 # 4,4 352 353 354 # add the portion for the first constraint 355 values[2] = values[2] + hessian_lambda[1] * (x[3] * x[4]) # 2,1 356 357 values[4] = values[4] + hessian_lambda[1] * (x[2] * x[4]) # 3,1 358 values[5] = values[5] + hessian_lambda[1] * (x[1] * x[4]) # 3,2 359 360 values[7] = values[7] + hessian_lambda[1] * (x[2] * x[3]) # 4,1 361 values[8] = values[8] + hessian_lambda[1] * (x[1] * x[3]) # 4,2 362 values[9] = values[9] + hessian_lambda[1] * (x[1] * x[2]) # 4,3 363 364 # add the portion for the second constraint 365 values[1] = values[1] + hessian_lambda[2] * 2 # 1,1 366 values[3] = values[3] + hessian_lambda[2] * 2 # 2,2 367 values[6] = values[6] + hessian_lambda[2] * 2 # 3,3 368 values[10] = values[10] + hessian_lambda[2] * 2 # 4,4 369 370 return ( values ) 371} 372@ 373After the hard part is done, we only have to define the initial values, the lower and upper bounds of the control variables, and the lower and upper bounds of the constraints. If a variable or a constraint does not have lower or upper bounds, the values \texttt{-Inf} or \texttt{Inf} can be used. If the upper and lower bounds of a constraint are equal, Ipopt recognizes this as an equality constraint and acts accordingly. 374<<solveHS071>>= 375# initial values 376x0 <- c( 1, 5, 5, 1 ) 377 378# lower and upper bounds of control 379lb <- c( 1, 1, 1, 1 ) 380ub <- c( 5, 5, 5, 5 ) 381 382# lower and upper bounds of constraints 383constraint_lb <- c( 25, 40 ) 384constraint_ub <- c( Inf, 40 ) 385 386opts <- list("print_level"=0, 387 "file_print_level"=12, 388 "output_file"="hs071_nlp.out") 389 390print( ipoptr( x0=x0, 391 eval_f=eval_f, 392 eval_grad_f=eval_grad_f, 393 lb=lb, 394 ub=ub, 395 eval_g=eval_g, 396 eval_jac_g=eval_jac_g, 397 constraint_lb=constraint_lb, 398 constraint_ub=constraint_ub, 399 eval_jac_g_structure=eval_jac_g_structure, 400 eval_h=eval_h, 401 eval_h_structure=eval_h_structure, 402 opts=opts) ) 403@ 404 405\section{Using data} 406The final subject we have to cover, is how to pass data to an objective function or the constraints. There are two ways to do this. The first is to supply additional arguments to the user defined functions and \texttt{ipoptr}. The second way is to define an environment that holds the data and pass this environment to \texttt{ipoptr}. Both methods are shown in \texttt{tests/parameters.R}. 407 408As a very simple example\footnote{A more interesting example is given in \texttt{tests/lasso.R}} suppose we want to find the minimum of 409$410f( x ) = a_1 x^2 + a_2 x + a_3 411$ 412for different values of the parameters$a_1$,$a_2$and$a_3$. 413 414First we define the objective function and its gradient using, assuming that there is some variable \texttt{params} that contains the values of the parameters. 415<<defineDataFunction1>>= 416eval_f_ex1 <- function(x, params) { 417 return( params[1]*x^2 + params[2]*x + params[3] ) 418} 419eval_grad_f_ex1 <- function(x, params) { 420 return( 2*params[1]*x + params[2] ) 421} 422@ 423Note that the first parameter should always be the control variable. All of the user-defined functions should contain the same set of additional parameters. You have to supply them as input argument to all functions, even if you're not using them in some of the functions. 424 425Then we can solve the problem for a specific set of parameters, in this case$a_1=1$,$a_2=2$and$a_3=3$, from initial value$x_0=0$, with the following command 426<<solveDataEnvironment1>>= 427# solve using ipoptr with additional parameters 428ipoptr( x0 = 0, 429 eval_f = eval_f_ex1, 430 eval_grad_f = eval_grad_f_ex1, 431 params = c(1,2,3) ) 432@ 433 434For the second method, we don't have to supply the parameters as additional arguments to the function. 435<<defineDataFunction2>>= 436eval_f_ex2 <- function(x) { 437 return( params[1]*x^2 + params[2]*x + params[3] ) 438} 439eval_grad_f_ex2 <- function(x) { 440 return( 2*params[1]*x + params[2] ) 441} 442@ 443Instead, we define an environment that contains specific values of \texttt{params} 444<<defineDataEnvironment2>>= 445# define a new environment that contains params 446auxdata <- new.env() 447auxdata$params <- c(1,2,3)
448@
449To solve this we supply \texttt{auxdata} as an argument to \texttt{ipoptr}, which will take care of evaluating the functions in the correct environment, so that auxiliary data is available.
450<<solveDataEnvironment2>>=
451# pass the environment that should be used to evaluate functions to ipoptr
452ipoptr( x0                 = 0,
453        eval_f             = eval_f_ex2,
455        ipoptr_environment = auxdata )
456@
457
458\section{Options}
459There are many options available, all of which are described on the Ipopt website. One of the options can test whether your derivatives are correct. This option is activated by setting \texttt{derivative\_test} to \texttt{first-order} or \texttt{second-order} if you want to test second derivatives as well. This process can take quite some time. To see all the output from this process you can set \texttt{derivative\_test\_print\_all} to \texttt{yes}, preferably when writing to a file, because of the problems with displaying on some terminals mentioned above. Without this optios the derivative checker only shows those lines where an error occurs if a high enough \texttt{print\_level} is supplied.
460
461\section{Remarks}
462If you run many large optimization problems in a row on Windows, at some point you'll get errors that Mumps is running out of memory and you won't get any solutions. On Linux this same problem hasn't occurred yet.
463
464The R terminal in Windows doesn't show any output. The Linux one does.
465
466\bibliographystyle{apacite}
467\bibliography{reflist}
468
469\end{document}
Note: See TracBrowser for help on using the repository browser.