R quadratic multi dimensional optimization with constraints

Question

I am tryinig to solve the quadratic optimization of the following function:

b<-4.7e-09
a<-(-2e-05)
M<-100
beta<-0.5
min<-fuction(x){
    x1=x[1]
    x2=x[2]
    x3=x[3]
    E=a*x1+b*x1^2+a*x2+b*x2^2+a*x3+b*x3^2
    V=(M-x1)+(M-x1-x2)+(M-x1-x2-x3)
    return (E+beta*V)
}

And the constrains are

x1+x2+x3=M
x1>=0,x2>=0,x3>=0

Is there a way I can solve this problem using constrOptim or solve.QP?

Not necessarily, but integral optimization would be better.

I appreciate any comment.

Answer 1

[Note: In what follows I call your function f(...) to avoid confusion with the built-in R function min(...). Also, I'm assuming that x2=x[3] in your code is an error, and you want x2=x[2].]

First, before you resort to numerical optimization, you should do some basic math. If x_i ≥ 0 and sum(x) = M, then x_i ≤ M. So we are operating in a cube with sides (0,M). Further, if sum(x) = M then we really have only 2 independent variables (say x₁ and x₂) and x₃ = M - (x₁ + x₂). We can determine the minimum relatively easily this way:

x    <- seq(0,M,len=101)
df   <- expand.grid(x=x,y=x)
df$f <- mapply(function(x,y) f(c(x,y,M-(x+y))),df$x,df$y)
df$f <- ifelse(df$x+df$y>M,NA,df$f)
df[which.min(df$f),]
#       x y         f
# 101 100 0 -0.001953

So the minimum of f occurs at x₁=M, x₂ = x₃ = 0.

Since the function f(...), is a surface, we can plot this to confirm, as follows (it is always a good idea to plot the function if at all possible!!).

library(reshape2)     # for dcast(...)
library(rgl)          # for surface3d(...), etc.
z       <- dcast(df,x~y,value.var="f")[-1]
zlim    <- range(z[!is.na(zz)])
palette <- rev(heat.colors(10))
col     <- palette[9*(df$f-zlim[1])/diff(zlim) + 1] 
surface3d(x,x,as.matrix(zz),color=col)
axes3d()
title3d(xlab="X",ylab="Y",zlab="Z")

So the surface turns out to be a plane, and the minimum is indeed at (100,0,0).

Finally, we can of course use a numerical optimizer (which IMO is overkill for this problem - unless of course this is a homework assignment??). Here we use nloptr(...) from the package of the same name. f(...) is the function to be minimized, g(...) is the constraint expressed as an inequality, abs(sum(x)-M) <= 0. and lb is a vector of the lower bounds on x. You can also specify the constraint as an equality using eval_g_eq=.... Read the documentation for more details.

f <-function(x){                   # objective function
  x1=x[1]
  x2=x[2]
  x3=x[3]
  E=a*x1+b*x1^2+a*x2+b*x2^2+a*x3+b*x3^2
  V=(M-x1)+(M-x1-x2)+(M-x1-x2-x3)
  return (E+beta*V)
}
g <- function(x) {abs(sum(x)-M)}   # constraint function

library(nloptr)
result <-nloptr(c(0,0,0), f, lb=c(0,0,0), eval_g_ineq=g,
                opts = list(algorithm="NLOPT_LN_COBYLA"))
result$solution
# [1] 1.000000e+02 4.440892e-16 4.835780e-16