nloptr算法xtol_rel终止条件的直觉解释问询及实例

一、先搞懂xtol_rel终止条件到底是什么 #

说白了，这个条件就是判断参数的相对变化量是否足够小：

• 公式拆解：|当前参数值 - 上一轮参数值| < xtol_rel * |上一轮参数值|，换算成比例就是|(当前值-上一轮值)/上一轮值| < xtol_rel。
• 拿你的实例举例子：比如x[1]初始是10000，你设的xtol_rel=1e-7，那允许的最大变化量就是10000*1e-7=0.001。如果迭代到某一步，x[1]从10000变成10000.0005，变化量0.0005小于0.001，那这个参数就达标了。
• 所有参数都满足这个条件时，算法就会终止。

二、为什么这个条件能判定局部最优？ #

核心逻辑就是：如果还能找到更优的参数组合，算法绝不会只迈这么小的步子。

1. 优化的本质是找目标函数的极值点：在局部最优附近，目标函数的“坡度”几乎为0——不管往哪个方向微调参数，目标函数值都不会有明显提升（你的例子是最大化目标，所以是不会再变大）。
2. 参数变化极小，说明算法已经“摸不到更好的方向”了：要是还有能让目标函数变好的路径，算法会直接调整参数到那个方向，步子不可能这么小。
3. 用相对变化而非绝对变化的原因：避免参数量级干扰。比如参数是10000时，变化1不算大；但参数是0.001时，变化1就完全是另一个量级了。用相对比例，能统一不同大小参数的判断标准，不会因为参数本身的数值大小误判。

你的优化实例代码 #

eval_f0 <- function(x){
  return(-((plogis(-2.15083 + 0.0000471*x[1])*exp(9.87243)*(x[1])^0.468)+
             (plogis(-2.15083 + 0.0000471*x[2])*exp(9.87243)*(x[2])^0.468)+
             (plogis(-2.15083 + 0.0000471*x[3])*exp(9.87243)*(x[3])^0.468)+
             (plogis(-2.15083 + 0.0000471*x[4])*exp(9.87243)*(x[4])^0.468)+
             (plogis(-1.86183 + 0.0000471*x[5])*exp(9.04243)*(x[5])^0.468)+
             (plogis(-1.86183 + 0.0000471*x[6])*exp(9.04243)*(x[6])^0.468)+
             (plogis(-1.86183 + 0.0000471*x[7])*exp(9.04243)*(x[7])^0.468)+
             (plogis(-1.86183 + 0.0000471*x[8])*exp(9.04243)*(x[8])^0.468)+
             (plogis(-1.06983 + 0.0000471*x[9])*exp(9.23543)*(x[9])^0.468)+
             (plogis(-1.06983 + 0.0000471*x[10])*exp(9.23543)*(x[10])^0.468)+
             (plogis(-1.06983 + 0.0000471*x[11])*exp(9.23543)*(x[11])^0.468)+
             (plogis(-1.06983 + 0.0000471*x[12])*exp(9.23543)*(x[12])^0.468)))
}

# 约束函数
eval_g0 <- function(x)
{
  return(x[1]+x[2]+x[3]+x[4]+x[5]+x[6]+x[7]+x[8]+x[9]+x[10]+x[11]+x[12]-50000)
}

# 设置优化选项
opts <- list( "algorithm"= "NLOPT_LN_COBYLA",
              "xtol_rel"= 1.0e-7,
              "maxeval"= 1000000,
              "print_level" = 3 )

# 求解
res2 <- nloptr( x0=c(10000,3666,2000,1000,10000,3666,2000,1000,10000,3666,2000,1000),
                eval_f=eval_f0,
                lb = c(0,0,0,0,0,0,0,0,0,0,0,0),
                ub = c(80000,80000,80000,80000,80000,80000,80000,80000,80000,80000,80000,80000),
                eval_g_ineq = eval_g0,
                opts = opts)

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