I have the following optimisation problem:
Maximisise combined benefit functions
F(x) = 2*10^9 + 160*x - 7*x*ln(x)
F(y) = 3*10^9 + 170*y - 7.5*y*ln(y)
subject to
x + y <= 3*10^9
x,y >= 0
I set up a Lagrangian
L(x,y,λ) = F(x) + F(y) + λ(3*10^9-(x+y))
I then reduced it algebraically before solving it with a graphing calculator
y = 3*10^9 - x
e^(19/15)*x^(14/15) + x - 3*10^9 = 0
y^(15/14) + e^(19/14)*y - e^(19/14)*3*10^9 = 0
x = 1.61 * 10^9
y = 1.39 * 10^9
Now, I would like to be able to achieve this solution in R because the actual function parameters may change. I have tried adapting a few different solutions I found online, but none seem to work.
If I understand the problem correctly, I need a non-linear solver. I therefore set up the problem in Rsolnp
as such (inspired from this answer):
library(Rsolnp)
opt_func_log <- function(x) {
a <- x[1]
b <- x[2]
ben_func <- 2e9 + 160*a - 7*a*log(a) + 3e9 + 170*b - 7.5*b*log(b)
-ben_func #invert to find minimum
}
equal_const <- function(x) {
a <- x[1]
b <- x[2]
a + b # budget constraint formula
}
eps <- .Machine$double.eps*10^2 # low number, but not zero due to logs
x0 <- c(0.1, 0.1) # starting values
budget <- 3e9 # overall budget constraint value
opt_solution_log <- solnp(pars = x0,
fun = opt_func_log,
eqfun = equal_const,
eqB = budget,
LB = c(eps,eps))
Unfortunately I don't get a viable solution
Iter: 1 fn: -5032442923.2173 Pars: 213333.43335 213333.43335
solnp-->Redundant constraints were found. Poor
solnp-->intermediate results may result.Suggest that you
solnp-->remove redundant constraints and re-OPTIMIZE
Iter: 2 fn: -5032442923.2173 Pars: 213333.43335 213333.43335
solnp--> Solution not reliable....Problem Inverting Hessian.
What am I doing wrong? What constraint is redundant in this problem? Have I defined the problem wrongly or is it just not solvable in this way?
Below are two options to solve the constrained optimization problem.
You can use constrOptim
if you need to set up linear constraints, e.g.,
ui <- rbind(-c(1, 1), c(1, 0), c(0, 1))
ci <- c(-3e9, 0, 0)
theta <- runif(2)
constrOptim(theta,
f = opt_func_log,
grad = NULL,
ui = ui,
ci = ci,
control = list(reltol = .Machine$double.eps)
)
and you will obtain
$par
[1] 1609771281 1390228719
$value
[1] -40508597161
$counts
function gradient
754 NA
$convergence
[1] 0
$message
NULL
$outer.iterations
[1] 3
$barrier.value
[1] -6339423
fmincon
from package pracma
library(pracma)
fmincon(runif(2),
opt_func_log,
A = t(c(1, 1)),
b = 3e9,
lb = c(0, 0),
tol = .Machine$double.eps
)
which gives
$par
[1] 1609771176 1390228824
$value
[1] -40508597161
$convergence
[1] 0
$info
$info$lambda
$info$lambda$lower
[,1]
[1,] 0
[2,] 0
$info$lambda$upper
[,1]
[1,] 0
[2,] 0
$info$lambda$ineqlin
[1] 4.604494 0.000000 0.000000
$info$grad
[,1]
[1,] -4.604494
[2,] -4.604494
$info$hessian
[,1] [,2]
[1,] 1.778707e+00 5.176693e-09
[2,] 5.176693e-09 2.693166e-09