Solve non-linear optimization problem with ROI package in R

Question

I want solve a non-linear objective function with two non-linear constraints in R. I have succesfully solved the model in GAMS but because of my general data workflow, I would like to use R instead. The idea behind the model is that I would like to calibrate a rational function (four parameters) so that it is consistent with a given dataset (resulting in two constraints), while minimizing the difference between the slope of the function and a target value (objective function). Hence the outcome of the optimization problem are the four parameters that determine the shape of the rational function. I am using the R Optimization Infrastructure (ROI) to solve the model. As I want to solve the model for several different cases, I created a series of functions that accept variables so the objective function and constraints can easily be updated.

My model:

library(ROI)

# I installed several non-linear solvers that were suggested after running ROI_applicable_solvers() on the model
#install.packages("ROI.plugin.alabama")
#install.packages("ROI.plugin.nlminb")
#install.packages("ROI.plugin.nloptr")

# Empirical data for the example
elas_sl_t <- -0.006648
elas_by_t <- 0.2831
gdp_cap_rel_t <- 51.69
cal_cap_t <- 27.75

# Rational function for which parameters x[] need to be solved
elas_f <- function(x, gdp_cap) {
  (x[1]*gdp_cap + x[2])/((gdp_cap+x[3])*(gdp_cap+x[4]))
}

# Integral of elas_f
int_elas_f <- function(x, gdp_cap) {
  -(x[1] * x[4] - x[2]) * log(gdp_cap + x[4])/(x[4]^2 - x[3] * x[4]) - 
    (x[2] - x[1] * x[3]) * log(gdp_cap + x[3])/(x[3] * x[4] - x[3]^2) + x[2] * log(gdp_cap)/(x[3] * x[4])
}

# Derivative of elas_f
der_elas_f <- function(x, gdp_cap){
  - (x[1] * gdp_cap + x[2]) / ((gdp_cap + x[3])^2 * (gdp_cap + x[4])) -
    (x[1] * gdp_cap + x[2]) / ((gdp_cap + x[3]) * (gdp_cap + x[4])^2) + x[1] / ((gdp_cap + x[3]) * (gdp_cap + x[4]))
}

# Objective function (should be minimized)
obj_f <- function(x, gdp_cap = 1, elas_sl_by = elas_sl_t){
  (der_elas_f(x, gdp_cap) - elas_sl_by)^2
}

# 1st constraint (should be equal to zero)
con1_f <- function(x, gdp_cap = 1, elas_by = elas_by_t){
  elas_f(x, gdp_cap) - elas_by
}

# 2nd constraint (should be equal to zero)
con2_f <- function(x, gdp_cap_ty = gdp_cap_rel_t, gdp_cap_by = 1, cal_ty = cal_cap_t) {
  int_elas_f(x, gdp_cap = gdp_cap_ty) - int_elas_f(x, gdp_cap = gdp_cap_by) - log(cal_ty)
}

# The solution of the model is
param_t <- c(1245.685,  30.987, 883.645,  4.097)

# check the solution
con1_f(param_t) # close to zero, ok
con2_f(param_t) # close to zero, ok
obj_f(param_t) # 0.05155

# Solving the model with ROI, using the actual solution as starting values does not work.
nlp <- OP(F_objective(F = obj_f, n = 4), 
          F_constraint(F = list(con1_f, con2_f), dir = rep("==", 2), rhs = c(0,0)),
          maximum = FALSE)
nlp
sol <- ROI_solve(nlp, start = param_t, solver = "auto")
sol
solution(sol)

It seems that the model cannot be solved or perhaps I am doing something wrong. Is it possible to solve the above problem using R (with ROI or perhaps an other package)? Perhaps the available solvers are not good enough? I used the GAMS CONOPT solver to solve my problem, which is not available for ROI/R. I also managed the solve the problems with the GAMS IPOPT solver, which is available for R (https://github.com/coin-or/Ipopt) but unfortunately not as a plugin for ROI (only IPOP, which is a quadratic solver). I also tried the NEOS solver in ROI but this results in an error message as the server cannot be contacted. Any help is much appreciated.

Answer 1

I'm not sure about the "solution" you got from GAMS. Also, what is considered "close to zero"?

# Modified Integral of elas_f
int_elas_f <- function(x, gdp_cap) {
  if (any(gdp_cap + c(0, x[3:4]) < 0)) return(NA_real_)
  -(x[1] * x[4] - x[2]) * log(gdp_cap + x[4])/(x[4]^2 - x[3] * x[4]) - 
    (x[2] - x[1] * x[3]) * log(gdp_cap + x[3])/(x[3] * x[4] - x[3]^2) + x[2] * log(gdp_cap)/(x[3] * x[4])
}

with(
  optim(
    c(1e3, 100, 1e3, 10),
    function(x) {
      y <- obj_f(x)
      y + 1e3*(abs(y) + 1)*(abs(con1_f(x)) + abs(con2_f(x)))
    }
  ), list(par = par, constraints = c(con1_f(par), con2_f(par)), obj = obj_f(par), value = value)
)
#> $par
#> [1] 1361.293833  400.099856  880.053420    6.061782
#> 
#> $constraints
#> [1]  1.676616e-08 -2.059692e-08
#> 
#> $obj
#> [1] 0.0342367
#> 
#> $value
#> [1] 0.03427534

The results are quite sensitive to the initial guess. A more comprehensive exploration indicates the solution may be divergent:

library(data.table)

f <- function(x) {
    with(
      optim(
        x,
        function(x) {
          y <- obj_f(x)
          y + 1e3*(abs(y) + 1)*(abs(con1_f(x)) + abs(con2_f(x)))
        }, method = "Nelder-Mead"
      ),
      c(con1_f(par), con2_f(par), obj_f(par))
    )
  }

m <- RcppAlgos::permuteGeneral(2^(8:20), 4)
n <- nrow(m)
df <- data.frame(con1 = numeric(n), con2 = numeric(n), obj = numeric(n))

for (i in 1:n) df[i,] <- f(m[i,])

dt <- setorder(setDT(df)[,r := .I][abs(con1) < 1e-6 & abs(con2) < 1e-6], obj)
m[dt$r[1],]
#> [1] 1048576    8192  131072    1024

sol <- with(
  optim(
    m[dt$r[1],],
    function(x) {
      y <- obj_f(x)
      y + 1e3*(abs(y) + 1)*(abs(con1_f(x)) + abs(con2_f(x)))
    }
  ),
  list(
    par = par,
    constraints = c(con1_f(par), con2_f(par)),
    obj = obj_f(par),
    value = value
  )
)

with(
  optim(
    sol$par,
    function(x) {
      y <- obj_f(x)
      y + 1e5*(abs(y) + 1)*(abs(con1_f(x)) + abs(con2_f(x)))
    }
  ),
  list(
    par = par,
    constraints = c(con1_f(par), con2_f(par)),
    obj = obj_f(par),
    value = value
  )
)
#> $par
#> [1]  4868827.420 24021111.619    30753.161     3317.202
#> 
#> $constraints
#> [1] -3.325118e-14  2.353673e-14
#> 
#> $obj
#> [1] 0.002944623
#> 
#> $value
#> [1] 0.002944629

As I continue to increase the initial values in the search space (e.g., m <- RcppAlgos::permuteGeneral(2^(12:24), 4)), the parameter values get increasingly larger while the value of the objective function continues to shrink.

The suspicion of non-convergence is supported by ROI, which has the sense to stop by the time it gets to 3.795066e-03:

 library(ROI)
 
 param_t <- c(1245.685,  30.987, 883.645,  4.097)

 nlp <- OP(F_objective(F = obj_f, n = 4), 
           F_constraint(F = list(con1_f, con2_f), dir = rep("==", 2), rhs = c(0,0)),
           maximum = FALSE)
 sol <- ROI_solve(nlp, start = param_t, solver = "auto")
#> Warning in nlminb2(start = start, objective = objective(x), eqFun = eqfun, :
#> gradient not applicable for *constrained* NLPs for solver 'nlminb'.
 sol
#> No optimal solution found.
#> The solver message was: No solution.
#> The objective value is: 3.795066e-03
 sol$solution
#> [1] 3.262698e+06 1.350195e+07 1.333212e+02 4.250796e+05
 sol$status$msg$symbol
#> [1] "NON_CONVERGENCE"
 sol$message$message
#> [1] "false convergence (8)"