Mark Miller - 2 months ago 8
R Question

# constrained optimization R: another example

I am attempting to perform constrained optimization in R. I have looked at these posts and a couple of others:

constrained optimization in R

function constrained optimization in R

The first post above is quite helpful, but I am still not obtaining the correct answer to my problem.

My function is:

``````Fd <- 224 * d1 + 84 * d2 + d1 * d2 - 2 * d1^2 - d2^2
``````

and my constraint is:
`3 * d1 + d2 = 280`

First I find the correct answer using an unconstrained exhaustive search followed by a constrained exhaustive search:

``````my.data <- expand.grid(x1 = seq(0, 200, 1), x2 = seq(0, 200, 1))
dim(my.data)

d1     <- my.data[,1]
d2     <- my.data[,2]

Fd <- 224 * d1 + 84 * d2 + d1 * d2 - 2 * d1^2 - d2^2

new.data <- data.frame(Fd = Fd, d1 = d1, d2 = d2)

# identify values of d1 and d2 that maximize Fd without the constraint
new.data[new.data\$Fd == max(new.data\$Fd),]
# **This is the correct answer**
#         Fd d1 d2
# 6157 11872 76 80

# Impose constraint
new.data <- new.data[(3 * new.data\$d1 + new.data\$d2) == 280, ]

# identify values of d1 and d2 that maximize Fd with the constraint
new.data[new.data\$Fd == max(new.data\$Fd),]
# **This is the correct answer**
#          Fd d1 d2
# 14743 11774 69 73
``````

Now find unconstrained maxima using
`optim`
. This works.

``````    Fd <- function(betas) {

b1 = betas[1]
b2 = betas[2]

(224 * b1 + 84 * b2 + b1 * b2 - 2 * b1^2 - b2^2)

}

# unconstrained
optim(c(60, 100), Fd, control=list(fnscale=-1), method = "BFGS", hessian = TRUE)
# \$par
# [1] 75.99999 79.99995
``````

Now find constrained maxima using
`constrOptim`
. This does not work.

``````b1.lower.bound <- c(0, 280)
b1.upper.bound <- c(93.33333, 0)
b2.lower.bound <- c(93.33333, 0)
b2.upper.bound <- c(0, 280)

theta = c(60,100)                         # starting values
ui = rbind(c(280,0), c(0,93.33333))       # range of allowable values
theta %*% ui                              # obtain ci as -1 * theta %*% ui
#       [,1]     [,2]
# [1,] 16800 9333.333

constrOptim(c(60,100), Fd, NULL, ui = rbind(c(280,0), c(0,93.33333)), ci = c(-16800, -9333.333), control=list(fnscale=-1))
# \$par
# [1] 75.99951 80.00798
``````

I have tried playing around with
`ui`
and
`ci`
, but it seems like no matter what values I use for them I always get the same answer as with unconstrained
`optim`
.

`constrOptim()` uses linear inequality constraints and defines the feasible region by `ui %*% param - ci >= 0`. If the constraint is `3 * d1 + d2 <= 280`, `ui` is `c(-3, -1)` and `ci` is `-280`.

``````Fd <- function(betas) {
b1 = betas[1]
b2 = betas[2]
(224 * b1 + 84 * b2 + b1 * b2 - 2 * b1^2 - b2^2)
}

theta = c(59.999,100)    # because of needing `ui %*% inital_par - ci > 0`
ui = c(-3, -1)
ci = -280                # they mean -3*par[1] + -1*par[2] + 280 >= 0

constrOptim(theta, Fd, NULL, ui = ui, ci = ci, control=list(fnscale=-1))
# \$par
# [1] 69.00002 72.99993
``````
[Edited]

If you want not inequality but equality constraints, it would be better to use `Rsolnp` or `alabama` package.

``````library(Rsolnp); library(alabama);

Fd2 <- function(betas) {     #  -1 * Fd
b1 = betas[1]
b2 = betas[2]
-1 * (224 * b1 + 84 * b2 + b1 * b2 - 2 * b1^2 - b2^2)
}

eqFd <- function(betas) {  # the equality constraint
b1 = betas[1]
b2 = betas[2]
(3 * b1 + b2 -280)
}

solnp(pars = c(60, 100), fun = Fd2, eqfun = eqFd, eqB = 0)
auglag(par = c(60, 100), fn = Fd2, heq = eqFd)
``````
Source (Stackoverflow)