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

Question

我想在 R 中求解具有两个非线性约束的非线性目标函数。我已经成功地在 GAMS 中求解了模型，但由于我的一般数据工作流程，我想改用 R。该模型背后的想法是，我想校准一个有理函数（四个参数），使其与给定的数据集一致（导致两个约束），同时最小化函数斜率与目标值之间的差异（目标函数）。因此，优化问题的结果是确定有理函数形状的四个参数。我正在使用 R 优化基础设施 (ROI) 来求解模型。由于我想针对几种不同的情况求解模型，因此我创建了一系列接受变量的函数，以便可以轻松更新目标函数和约束。

我的型号：

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)

似乎模型无法求解，或者可能是我做错了什么。是否可以使用 R（使用 ROI 或其他包）解决上述问题？也许可用的求解器不够好？我使用 GAMS CONOPT 求解器来解决我的问题，该求解器不适用于 ROI/R。我还设法使用 GAMS IPOPT 求解器解决问题，该求解器可用于 R (https://github.com/coin-or/Ipopt)，但不幸的是不能作为 ROI 的插件（仅限 IPOP，这是一个二次求解器）。我还尝试了 ROI 中的 NEOS 求解器，但这会导致错误消息，因为无法联系服务器。非常感谢任何帮助。

Answer 1

我不确定您从 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