[R] Optimization with constraint.

[Paul Smith]

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