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.
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)"