使用 R 中的 ROI 包解决非线性优化问题

我不确定你从GAMS得到的"解决方案".另外，什么被认为是"接近零"？

# 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

结果对最初的猜测相当敏感.一项更全面的探索表明，解决方案可能存在分歧:

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 <- Rcpp算法s::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

随着我继续增加搜索空间中的初始值(例如，m <- Rcpp算法s::permuteGeneral(2^(12:24), 4))，参数值变得越来越大，而目标函数的值继续缩小.

不收敛的怀疑得到了ROI的支持，当它达到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)"