From: Paul Smith <phhs80_at_gmail.com>

Date: Sat, 15 Mar 2008 11:50:10 +0000

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 15 Mar 2008 - 11:52:48 GMT

Date: Sat, 15 Mar 2008 11:50:10 +0000

On Fri, Mar 14, 2008 at 6:42 PM, Hans W Borchers <hwborchers_at_gmail.com> wrote:

> > I have some problems, when I try to model an

*> > optimization problem with some constraints.
**> >
**> > The original problem cannot be solved analytically, so
**> > I have to use routines like "Simulated Annealing" or
**> > "Sequential Quadric Programming".
**> >
**> > But to see how all this works in R, I would like to
**> > start with some simple problem to get to know the
**> > basics:
**> >
**> > The Problem:
**> > min f(x1,x2)= (x1)^2 + (x2)^2
**> > s.t. x1 + x2 = 1
**> >
**> > The analytical solution:
**> > x1 = 0.5
**> > x2 = 0.5
**> >
**> > Does someone have some suggestions how to model it in
**> > R with the given functions optim or constrOptim with
**> > respect to the routines "SANN" or "SQP" to obtain the
**> > analytical solutions numerically?
**> >
**>
**> In optimization problems, very often you have to replace an equality by two
**> inequalities, that is you replace x1 + x2 = 1 with
**>
**> min f(x1,x2)= (x1)^2 + (x2)^2
**> s.t. x1 + x2 >= 1
**> x1 + x2 <= 1
**>
**> The problem with your example is that there is no 'interior' starting point for
**> this formulation while the documentation for constrOptim requests:
**>
**> The starting value must be in the interior of the feasible region,
**> but the minimum may be on the boundary.
**>
**> You can 'relax', e.g., the second inequality with x1 + x2 <= 1.0001 and use
**> (1.00005, 0.0) as starting point, and you will get a solution:
**>
**> >>> A <- matrix(c(1, 1, -1, -1), 2)
**> >>> b <- c(1, -1.0001)
**>
**> >>> fr <- function (x) { x1 <- x[1]; x2 <- x[2]; x1^2 + x2^2 }
**>
**> >>> constrOptim(c(1.00005, 0.0), fr, NULL, ui=t(A), ci=b)
**>
**> $par
**> [1] 0.5000232 0.4999768
**> $value
**> [1] 0.5
**> [...]
**> $barrier.value
**> [1] 9.21047e-08
**>
**> where the accuracy of the solution is certainly not excellent, but the solution
**> is correctly fulfilling x1 + x2 = 1.
*

One can get an accurate solution with optim via a penalty:

f <- function(x) {

x1 <- x[1]

x2 <- x[2]

return (-x1^2 - x2^2 + 1e10 * (x1 + x2 - 1)^2)

}

grad <- function(x) {

x1 <- x[1]

x2 <- x[2]

return (c(2e10*(x2+x1-1)-2*x1,2e10*(x2+x1-1)-2*x2))

}

optim(c(2,2),f,grad,control=list(maxit=1000),method="L-BFGS-B")

Paul

R-help_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sat 15 Mar 2008 - 11:52:48 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Sat 15 Mar 2008 - 12:30:21 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*