Solving an equation in R with given constraints

Question

I want to maximize the following function subject to the given constraints.

-p₁log(p₁) - p₃log(p₃) - p₅log(p₅)

subject to p₁ + p₃ + p₅ = 1 and p₁ + 3p₃ + 5p₅ = 3.5

p₁ , p₃ and p₅ all lie between 0 and 1 [They are probabilities].

My question is how do I solve this in R? From what I saw, constrOptim() is one of the functions commonly used to solve these type of problems. However I could not figure it out.

Any help is appreciated.

Answer 1

Package Rsolnp uses Lagrange multipliers to solve non-linear problems with equality constraints. Here is how it would be setup. eps is meant not to have the logarithms produce NaN values.

library(Rsolnp)

f <-function(X) {
  x <- X[1]
  y <- X[2]
  z <- X[3]
  res <- -x*log(x) - y*log(y) - z*log(z)
  -res
}
eq_f <- function(X){
  x <- X[1]
  y <- X[2]
  z <- X[3]
  c(
    x + y + z,
    x + 3*y + 5*z
  )
}

eps <- .Machine$double.eps*10^2
X0 <- c(0.1, 0.1, 0.1)

sol <- solnp(
  pars = X0,
  fun = f,
  eqfun = eq_f,
  eqB = c(1, 3.5),
  LB = rep(eps, 3),
  UB = rep(1, 3)
)
#
#Iter: 1 fn: -1.0512     Pars:  0.21624 0.31752 0.46624
#Iter: 2 fn: -1.0512     Pars:  0.21624 0.31752 0.46624
#solnp--> Completed in 2 iterations

sol$convergence
#[1] 0
sol$pars
#[1] 0.2162396 0.3175208 0.4662396
sol$values
#[1]  0.000000 -1.051173 -1.051173

The last value of sol$values is the function value at the optimal parameters.
We can check that the constraints are met.

sum(sol$pars)
#[1] 1
sum(sol$pars*c(1, 3, 5))
#[1] 3.5