R中nloptr-COBYLA约束不满足问题求助及代码排查

问题梳理 #

你用nloptr的COBYLA做极小化建模，目标函数已经收敛，但解却不满足约束——尤其是你尝试通过两个不等式转换的等式约束完全没生效，这大概率是约束转换的写法出了核心问题，咱们一步步拆解修复：

核心问题与修复方案 #

1. 等式约束转换错误（最关键） #

COBYLA确实需要把等式约束g(x)=0拆成g(x) ≤ 0和-g(x) ≤ 0两个不等式，但你的代码里h[7]和h[8]是完全一样的表达式！这等于只加了一个重复的约束，根本没实现等式限制。

看你这段代码：

h[7]<-sum(wp$use)*x[1]+sum(wm$use)*x[2]+sum(wo$use)*x[3]+sum(sp$use)*x[4]+sum(sm$use)*x[5]+sum(so$use)*x[6]+sum(bp$use)*x[7]+sum(bm$use)*x[8]+sum(bo$use)*x[9]-1.1*cost+cost 
h[8]<-(sum(wp$use)*x[1]+sum(wm$use)*x[2]+sum(wo$use)*x[3]+sum(sp$use)*x[4]+sum(sm$use)*x[5]+sum(so$use)*x[6]+sum(bp$use)*x[7]+sum(bm$use)*x[8]+sum(bo$use)*x[9]-1.1*cost+cost) 

先简化你的表达式：-1.1*cost + cost = -0.1*cost，所以你要实现的等式是sum(...) = 0.1*cost，正确的转换应该是：

# 先把等式表达式提出来，更清晰
eq_val <- sum(wp$use*x[1] + wm$use*x[2] + wo$use*x[3] + 
              sp$use*x[4] + sm$use*x[5] + so$use*x[6] + 
              bp$use*x[7] + bm$use*x[8] + bo$use*x[9]) - 0.1*cost
h[7] <- eq_val       # 对应 eq_val ≤ 0 → sum(...) ≤ 0.1*cost
h[8] <- -eq_val      # 对应 -eq_val ≤ 0 → sum(...) ≥ 0.1*cost

这样两个约束结合起来，才是sum(...) = 0.1*cost的等式约束，你之前的写法等于只限制了sum(...) ≤ 0.1*cost，没有下限，解自然会偏离等式。

2. 数值容忍度设置 #

COBYLA是近似梯度算法，对约束满足有数值容忍度。你当前只设置了xtol_rel = 1e-8，可以在control里额外添加constrtol = 1e-8（或更小的数值），收紧约束的检查标准，避免因为数值误差导致约束轻微违反。

3. 初始点与迭代次数 #

你的初始值x0全设为100，可能离可行域较远；虽然迭代了4620次达到xtol_rel，但COBYLA在复杂约束下可能需要更多迭代次数。可以尝试：

• 先手动找一个满足所有约束的初始点（比如计算一组符合x2≥x1、x3≥x2等约束的x值），再启动优化
• 把maxeval增大到20000甚至更高，给算法更多时间收敛到可行解

4. 简化表达式减少误差 #

你的目标函数和约束里有大量重复的sum(...)计算，不仅低效，还可能引入数值误差。建议提前预计算每个x的系数：

# 预计算每个x变量的系数，后续直接调用
coeffs <- c(sum(wp$use), sum(wm$use), sum(wo$use),
            sum(sp$use), sum(sm$use), sum(so$use),
            sum(bp$use), sum(bm$use), sum(bo$use))
target_eq <- 0.1*cost  # 等式目标值

之后在约束里直接用sum(coeffs * x) - target_eq，既清晰又减少重复计算的误差。

修正后的代码片段 #

调整约束部分后的完整示例：

# 预计算系数，提升效率和准确性
coeffs <- c(sum(wp$use), sum(wm$use), sum(wo$use),
            sum(sp$use), sum(sm$use), sum(so$use),
            sum(bp$use), sum(bm$use), sum(bo$use))
target_eq <- 0.1*cost  # 简化后的等式目标

eval_f<-function(x){ 
 -(CCC-a1*x[1]^(0.9)-a2*x[2]^(0.9)-a3*x[3]^(0.9)-b1*x[4]^(0.9) 
 -b2*x[5]^(0.9)-b3*x[6]^(0.9)-c1*x[7]^(0.9)-c2*x[8]^(0.9)-c3*x[9]^(0.9)) 
}

eval_g_ineq<-function(x){ 
 h<-numeric(8) 
 h[1]<-x[2]-x[1] 
 h[2]<-x[3]-x[2] 
 h[3]<-x[5]-x[4] 
 h[4]<-x[6]-x[5] 
 h[5]<-x[8]-x[7] 
 h[6]<-x[9]-x[8] 
 # 正确转换等式约束
 eq_expr <- sum(coeffs * x) - target_eq
 h[7] <- eq_expr       
 h[8] <- -eq_expr      
 return(h) 
}

x0<-c(100,100,100,100,100,100,100,100,100) 
# 调整控制参数，收紧约束容忍度
res<-cobyla(x0, eval_f, lower = lb, upper = NULL, hin = eval_g_ineq,
            control = list(xtol_rel = 1e-8, maxeval = 20000, constrtol = 1e-8))
print(res)

额外检查项 #

• 确认lb（变量下界）是否合理，有没有和x2≥x1这类约束冲突的情况
• 计算求解结果的eq_expr值：如果是轻微的数值偏差（比如1e-7级别），调整constrtol即可；如果是明显偏离，那肯定是约束转换的问题（就是你之前重复表达式的错误）

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