R optim() constraint optimization does not find the first best

Question

my problem is summarized in finding a vector X with the best solution to the problem: L is the profits, R is the restrictions, P is a constraint parameters matrix,

max SUM_i (x_i * l_i)

or max(t(L)%*%X)

restriction
- SUM_i(x_i*p_ij)<=r_j

or P%*%X <= R.

I find a solution for X, but not the best, which would be fb = c(.217,0,0,23,2865,0,13,427). How do I find the best solution?

code:

X<-matrix(rep(1,6),6,1)
P<-matrix(c(
    1, 1, 1, 2, 0, 0,
    0, 1, 1, 2, 1, 1,
    99.4, 37.75, 19.75, 54.40, 74.75, 53,
    2.400, 1.540, 0, 0, 0, 0,
    2.400, 1.960, 0, 0, 0, 0,
    1.800, 3.300, 5.330, 0, 0, 0,
    0, 0, 2.070, 0, 8.700, 0,
    0, 0, .436, 0, 19.100, 12.363,
    0, 3.000, .364, 0, 9.100, 26.737  ), 
9,6,1)

L <- matrix(c(83.4, 72.35, 27.3, 72.05, 217.25, 455), 6,1)
R <- matrix(c(60,60,2000,351,448,479,338,424,359),9,1)

farm<- function(par, P,R,  L){
   trues<- P%*%par<=R
 if (min(trues)==1 && min(par)>=0) {
    return(-t(L)%*%par)
     }
 else{
     return(0)
 }
}

mtds = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B", "SANN","Brent")
out <- optim(par = X,  # initial guess
             fn  = farm,
             P   = P,
             R   = R ,
             L   = L,
             method  = mtds[5]) 


# my result
t(L)%*%out$par
#A matrix: 1 × 1 of type dbl
#7419.596

# the first best
fb<- matrix(c(.217,0,0,23.2865,0,13.427),6,1)
t(L)%*%fb
#A matrix: 1 × 1 of type dbl
#7805.175

Is it the same result for all the methods? It could be that the solver gets stuck at a 'local' minimum rather than a 'global' minimum - this is fairly common. You can try a few different initial guesses to get the global min. Alternatively, you can also look at implementing a simple Newton-Raphson method here instead of using a library. — Gautam
No, it is the best solution between the methods. I tested them all from "Nelder-Mead" to "SANN" — Andre Suriane

ThomasIsCoding ThomasIsCoding · Accepted Answer · 2020-04-26T06:06:27

I think you can try fmincon from package pracma

library(pracma)
objfun <- function(x) -t(L)%*%x
res <- fmincon(x0 = X,fn = objfun,A = P,b = R,lb = rep(0,length(X)))

and you will see that

> res$par
[1]  4.201711e-16 -1.239088e-15  1.863081e-17  2.310286e+01
[5]  5.566620e-01  1.323762e+01

> -res$value
         [,1]
[1,] 7808.615

R optim() constraint optimization does not find the first best

2 Answers