R中constrOptim在特定约束下得次优解的原因及解决方法

问题描述 #

我使用R语言的constrOptim函数在两种投入间分配预算以实现收益最大化，理论上最优分配应使两种投入的边际收益相等，但在预算设为2.5e9时，该结论不成立；当预算设为2.4e9或2.6e9时，边际收益则能相等。2.5e9预算下的运行结果显示两种投入的边际收益分别为7.16和7.33，怀疑是精度或容差问题但无法修正，需明确问题原因及调整方案。

示例代码 #

# 定义收益函数
benefit_c <- function(x){
  1832283199  + 162.1422 * x - 6.8976 * x * log(x)
}

benefit_m <- function(x){
  303009100 + 100.593 * x - 4.486285 * x * log(x)
}

# 对应的边际收益函数
marg_c <- function(x){
  (162.1422 - 6.8976)  - 6.8976 * log(x)
}

marg_m <- function(x) {
  (100.593 - 4.486285)  - 4.486285  * log(x)
}

# 总收益函数（取负值用于最小化求解）
benefit_function <- function(x) {
  c <- x[1] 
  m <- x[2] 
  
  ben_func <- - (benefit_c(c) + benefit_m(m))
  names(ben_func) <- "tot_benefit"
  
  return(ben_func)
}

# 设置参数
budget <- 2.5e9

# 初始值
eps <- .Machine$longdouble.eps # 避免对数取0
init_pars <- c(0.1, 0.1)

# 约束矩阵
constr_input_matrix <- matrix(c(1,0,0,1,-1,-1),
                              nrow = 3,ncol = 2,byrow = TRUE)
constr_output_matrix <- matrix(c(eps,eps,-budget),
                               nrow = 3,ncol = 1,byrow = TRUE)

# 运行求解器
log_solution <- constrOptim(theta = init_pars,
                            f = benefit_function,
                            grad = NULL,
                            ui = constr_input_matrix,
                            ci = constr_output_matrix,
                            control = list(reltol = eps)
)

# 查看边际收益
marg_c(log_solution$par[1])
marg_m(log_solution$par[2])

问题原因 #

1. 算法收敛性与数值精度限制
   constrOptim默认的Nelder-Mead无梯度算法，在目标函数平坦区域的收敛精度易受参数影响。当预算为2.5e9时，最优解附近的目标函数可能处于极平坦区域，算法达到reltol设定的容差后提前终止，得到的只是近似最优解，未满足边际收益严格相等的条件。

2. 数值尺度差异过大
   预算值（2.5e9）与初始参数（0.1）的尺度差异极大，会导致求解过程中数值稳定性下降，干扰算法对最优解的定位。

3. 容差参数设置不合理
   使用.Machine$longdouble.eps作为相对容差，该值过小，算法无法达到如此极端的精度要求，反而会提前终止迭代或陷入数值不稳定状态，接受接近但非最优的解。

调整方案 #

1. 优化初始值 #

放弃极小初始值，通过边际收益相等的理论条件计算接近最优解的初始值，帮助算法更快收敛到最优区域：

# 根据边际收益相等+预算约束计算近似初始值
approx_init <- function(budget) {
  f <- function(c) marg_c(c) - marg_m(budget - c)
  c_opt <- uniroot(f, interval = c(1e6, budget - 1e6))$root
  return(c(c_opt, budget - c_opt))
}
init_pars <- approx_init(budget)

2. 调整求解器控制参数 #

• 把reltol调整到合理范围（如1e-8），同时增加最大迭代次数，让算法有足够迭代次数收敛到更优解：

log_solution <- constrOptim(theta = init_pars,
                            f = benefit_function,
                            grad = NULL,
                            ui = constr_input_matrix,
                            ci = constr_output_matrix,
                            control = list(reltol = 1e-8, maxit = 10000)
)

• 利用已有的解析边际收益函数提供梯度，大幅提升求解精度与效率：

# 定义总收益的负梯度（适配最小化求解）
grad_function <- function(x) {
  c <- x[1]
  m <- x[2]
  return(c(-marg_c(c), -marg_m(m)))
}

log_solution <- constrOptim(theta = init_pars,
                            f = benefit_function,
                            grad = grad_function,
                            ui = constr_input_matrix,
                            ci = constr_output_matrix,
                            control = list(reltol = 1e-8)
)

3. 统一变量尺度 #

将投入量与预算除以大常数（如1e6）缩小数值尺度，提升求解过程的数值稳定性：

scale_factor <- 1e6
budget_scaled <- budget / scale_factor

# 重新定义缩放后的收益与边际收益函数
benefit_c_scaled <- function(x){
  x_original <- x * scale_factor
  1832283199  + 162.1422 * x_original - 6.8976 * x_original * log(x_original)
}

benefit_m_scaled <- function(x){
  x_original <- x * scale_factor
  303009100 + 100.593 * x_original - 4.486285 * x_original * log(x_original)
}

marg_c_scaled <- function(x){
  x_original <- x * scale_factor
  (162.1422 - 6.8976)  - 6.8976 * log(x_original)
}

marg_m_scaled <- function(x) {
  x_original <- x * scale_factor
  (100.593 - 4.486285)  - 4.486285  * log(x_original)
}

# 缩放后的总收益函数与约束
benefit_function_scaled <- function(x) {
  c <- x[1] 
  m <- x[2] 
  ben_func <- - (benefit_c_scaled(c) + benefit_m_scaled(m))
  names(ben_func) <- "tot_benefit"
  return(ben_func)
}

init_pars_scaled <- approx_init(budget) / scale_factor
constr_output_matrix_scaled <- matrix(c(eps/scale_factor, eps/scale_factor, -budget_scaled),
                                      nrow = 3,ncol = 1,byrow = TRUE)

# 求解缩放问题并还原原始尺度
log_solution_scaled <- constrOptim(theta = init_pars_scaled,
                                  f = benefit_function_scaled,
                                  grad = function(x) c(-marg_c_scaled(x[1]), -marg_m_scaled(x[2])),
                                  ui = constr_input_matrix,
                                  ci = constr_output_matrix_scaled,
                                  control = list(reltol = 1e-8)
)
log_solution$par <- log_solution_scaled$par * scale_factor

内容的提问来源于stack exchange，提问作者Modern_ibex