在R中求解含对数收益函数的带约束非线性优化问题的报错排查

问题排查与核心原因分析 #

你遇到的报错其实是几个细节问题叠加导致的，并非约束真的冗余或者问题无解：

1. 初始值太离谱：你选的c(0.1, 0.1)和实际最优解的量级（1e9）差了9个数量级，优化器根本没法从这么小的起点收敛到正确区域；而且对数函数在接近0的地方导数极大，会严重干扰优化器的梯度计算逻辑。
2. 大数值的精度问题：变量x、y都是1e9级别，这么大的数值在计算导数、Hessian矩阵时容易出现数值精度损失，直接导致了“Hessian求逆失败”的报错。
3. 冗余约束提示是假阳性：这个提示不是说你的约束真的没用，而是优化器在初始值附近计算时，约束的梯度信息和目标函数梯度没形成有效约束方向，属于初始值不合理带来的误导。

修正后的解决方案 #

下面给你几个可行的修正方案，从简单到进阶：

方案1：调整初始值，继续用Rsolnp

先把初始值设成接近预算平分的合理值，让优化器从接近最优解的区域开始搜索，代码修改如下：

library(Rsolnp)

opt_func_log <- function(x) {
  a <- x[1]
  b <- x[2]
  # 计算总收益，取负转为最小化问题
  ben_func <- 2e9 + 160*a - 7*a*log(a) + 3e9 + 170*b - 7.5*b*log(b)
  -ben_func
}

equal_const <- function(x) {
  sum(x) # 直接求和更简洁
}

eps <- .Machine$double.eps*10^2 
# 初始值设为预算平分，贴近最优解区域
x0 <- c(1.5e9, 1.5e9) 
budget <- 3e9 

opt_solution_log <- solnp(pars = x0,
                          fun = opt_func_log,
                          eqfun = equal_const,
                          eqB = budget,
                          LB = c(eps, eps))

# 查看最优解
opt_solution_log$pars
# 计算对应的最大收益
-opt_func_log(opt_solution_log$pars)

运行后就能得到和你手动计算接近的结果：x≈1.61e9，y≈1.39e9，之前的报错也会消失。

方案2：变量缩放提升数值稳定性

为了彻底解决大数值的精度问题，你可以把x、y都除以1e9，转化为0-3之间的无量纲数，优化完成后再缩放回去，代码示例：

library(Rsolnp)

# 缩放因子
scale_factor <- 1e9

opt_func_scaled <- function(z) {
  a <- z[1] * scale_factor
  b <- z[2] * scale_factor
  ben_func <- 2e9 + 160*a - 7*a*log(a) + 3e9 + 170*b - 7.5*b*log(b)
  -ben_func
}

equal_const_scaled <- function(z) {
  sum(z) # 此时约束变为z1 + z2 = 3（对应原x+y=3e9）
}

x0_scaled <- c(1.5, 1.5) # 对应原初始值1.5e9
LB_scaled <- c(eps/scale_factor, eps/scale_factor)

opt_solution_scaled <- solnp(pars = x0_scaled,
                             fun = opt_func_scaled,
                             eqfun = equal_const_scaled,
                             eqB = 3,
                             LB = LB_scaled)

# 缩放回原变量量级
original_pars <- opt_solution_scaled$pars * scale_factor
original_pars

这种方法在参数变化、数值量级改变时，适配性会更强。

方案3：换用nloptr包（备选优化工具）

如果Rsolnp还是有奇怪的问题，你可以试试nloptr包，它的非线性优化接口更灵活，对大数值问题的鲁棒性也不错：

library(nloptr)

opt_func_nloptr <- function(x) {
  a <- x[1]
  b <- x[2]
  - (2e9 + 160*a - 7*a*log(a) + 3e9 + 170*b - 7.5*b*log(b))
}

# 等式约束：x1 + x2 = 3e9
constraints_eq <- function(x) {
  return(x[1] + x[2] - 3e9)
}

x0 <- c(1.5e9, 1.5e9)
lb <- c(eps, eps)

opts <- list("algorithm" = "NLOPT_LD_SLSQP",
             "xtol_rel" = 1e-8,
             "maxeval" = 1000)

result <- nloptr(x0 = x0,
                 eval_f = opt_func_nloptr,
                 eval_g_eq = constraints_eq,
                 lb = lb,
                 opts = opts)

# 查看最优解
result$solution

额外原理补充 #

你的原问题是凹函数的带约束最大化：收益函数F(x)的二阶导数是-7/x，F(y)是-7.5/y，都是负数，所以总收益函数是凹的，最优解必然出现在约束边界x+y=3e9上，所以用等式约束是完全正确的，不需要考虑不等式约束的情况。

备注：内容来源于stack exchange，提问作者Modern_ibex