Optimisation problem with constraint in R: how to solve problem with "complex" logarithmic benefit functions?

Question

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?

Answer 1

Below are two options to solve the constrained optimization problem.

---

If you work with base R

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

If you would like to use 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

Answer 2

I don't know the package you're using but I tried base R. optim() apparently only works with a single parameter it can be a vector of length > 1.

foo <- function(x) {
   big_number <- 10^12
  -(
    2*10^9 + 160*x[1] - 7*x[1]*log(x[1])  + 
    3*10^9 + 170*x[2] - 7.5*x[2]*log(x[2])
  ) + 
    if (sum(x) > 3*10^9) big_number else 0 
}
init <- c(1 * 10^9, 2 * 10^9)
optim(par = init, fn = foo)

# $par
# [1] 1610058872 1389941050
# 
# $value
# [1] -40508596398
# ...

Answer 3

Similar to @s_baldur's answer, but you can work the constraints with optim using, e.g., tanh to keep the values between 0 and 3e9 and max to keep the sum less than 3e9:

# function parameters
A <- c(2e9, 3e9)
B <- c(160, 170)
C <- c(-7, -7.5)

reparam <- function(xy) { # reparameterize (x, y)
  xy <- (tanh(xy) + 1)*15e8 # min, max constraint
  xy[2] <- xy[2] - max(0, sum(xy) - 3e9) # x + y <= 3*10^9 constraint
  xy
}

FF <- function(xy) {
  xy <- reparam(xy)
  -sum(A + B*xy + C*xy*log(xy))
}

sol <- optim(c(0, 0), FF)[1:2]
sol$par <- reparam(sol$par)
sol
#> $par
#> [1] 1609878971 1390121029
#> 
#> $value
#> [1] -40508597105