On Fri, Mar 14, 2008 at 6:42 PM, Hans W Borchers <[EMAIL PROTECTED]> 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.5, 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.5, 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@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.
