“复现”这两个字看起来简单做起来全是细节。我最近把模拟退火算法SA和粒子群算法PSO放到同一个约束最优化问题上重新跑了一遍发现核心流程很多书里都有但真正让人血压升高的地方不是算法本身而是约束怎么处理、参数怎么定、对比怎么才公平。这篇文章用经典的减速器/变速器设计优化问题作为研究对象记录我从模型构建到两种算法跑通的全过程。正在学智能优化算法、手头有带约束工程优化问题、或者想在论文里拿这两种算法做对比实验的同学应该能从里面拿走点直接能用的东西。1. 先把“变速器设计”翻译成一个带约束的优化模型1.1 为什么拿减速器设计来当复现对象题目里写的是“附图所示变速……”那张原图我没有拿到具体设计参数所以干脆采用工程优化领域最常用的带约束测试题减速器设计优化问题Speed Reducer Design Problem。这个问题的好处是它有公认的最优参考值变量和约束全是连续的而且每个约束都有清晰的物理含义用来验证SA和PSO的实现是否靠谱非常合适。顺便说一句减速器本身就是变速器的一种典型形式通过齿轮传动改变转速和扭矩。设计者最关心的目标之一是在满足强度、刚度和装配要求的前提下让整个齿轮箱尽量轻、尽量小。这正好是一个非线性、多约束、目标函数光滑但可行域边界复杂的优化问题。你完全可以把这套模型理解为一个标准模板以后再遇到自己的变速机构优化只需要把几何参数、载荷系数替换掉算法本身一个字都不用改。1.2 设计变量、目标函数和边界条件这个问题一共有7个设计变量每一个都对应减速器里的实际几何尺寸。这里我把变量含义和取值范围先列出来后面写代码时直接用这些数组定义边界。变量工程含义取值范围x1齿面宽度cm[2.6, 3.6]x2齿轮模数cm[0.7, 0.8]x3小齿轮齿数[17, 28]x4第一根轴轴承间距cm[7.3, 8.3]x5第二根轴轴承间距cm[7.3, 8.3]x6第一根轴直径cm[2.9, 3.9]x7第二根轴直径cm[5.0, 5.5]目标函数是最小化减速器的总重量主体由齿轮重量和轴的重量构成。写成代码就是下面这样import numpy as np def reducer(x): x1, x2, x3, x4, x5, x6, x7 x f (0.7854 * x1 * x2**2 * (3.3333 * x3**2 14.9334 * x3 - 43.0934) - 1.508 * x1 * (x6**2 x7**2) 7.4777 * (x6**3 x7**3) 0.7854 * (x4 * x6**2 x5 * x7**2)) return f前半部分跟齿轮的宽度、模数、齿数直接相关后半部分跟两根轴的直径和轴承间距相关。可以看到变量之间是强耦合的想减小齿轮重量往往会增大轴的受力反过来又逼着你把轴做粗轴一粗整个箱体重量又上去了。这种“按下葫芦浮起瓢”的特性正是约束最优化问题的典型难点。1.3 11条约束分别代表什么约束条件一共11个全部写成小于等于0的标准形式这样罚函数好统一处理。我把每条约束的工程含义整理成了表格方便对照理解。约束数学含义工程含义g127/(x1·x2²·x3) - 1 ≤ 0齿面接触应力不超过许用值g2397.5/(x1·x2²·x3²) - 1 ≤ 0齿根弯曲应力不超过许用值g31.93·x4³/(x2·x3·x6⁴) - 1 ≤ 0第一根轴的挠度限制g41.93·x5³/(x2·x3·x7⁴) - 1 ≤ 0第二根轴的挠度限制g5综合应力计算项/(110·x6³) - 1 ≤ 0第一根轴的弯曲应力限制g6综合应力计算项/(85·x7³) - 1 ≤ 0第二根轴的弯曲应力限制g7x2·x3/40 - 1 ≤ 0模数与齿数乘积不能过大g85·x2/x1 - 1 ≤ 0齿宽相对模数不能过窄g9x1/(12·x2) - 1 ≤ 0齿宽相对模数不能过宽g10(1.5·x6 1.9)/x4 - 1 ≤ 0第一轴轴径与轴承间距的装配关系g11(1.1·x7 1.9)/x5 - 1 ≤ 0第二轴轴径与轴承间距的装配关系前6条约束是强度与刚度约束决定这个减速器能不能安全工作后5条约束是几何与装配约束决定这个减速器能不能加工、能不能装起来。实际工程里约束往往比这还多但抽象成优化问题的套路是一样的。1.4 约束最优化问题的解法思路先惩罚后搜索SA和PSO本质上是无约束搜索算法它们并不知道“g1必须小于0”是什么意思算法只会看一个标量适应度值。所以处理约束最优化问题的核心思路是把约束“翻译”进评价值里面去常用的就是罚函数法。我采用的评价函数形式是F(x) f(x) C · Σ max(0, gᵢ(x))²C是惩罚系数gᵢ(x)0就说明该约束被违反平方之后让惩罚量随违规程度非线性增长。如果某个解违反约束很严重F会变得非常大算法自然倾向于放弃它。但这里有个微妙的地方C太小搜索过程会大量停留在不可行区域C太大算法又会被锁死在第一个遇到的可行区域里根本不敢往外探索。具体怎么平衡我在后面的调参部分会专门展开。这一章总结下来就是先把真实工程问题数学化再把约束吸收进评价函数之后SA和PSO才能在同一个框架下公平竞争。2. 模拟退火SA的完整复现路径温度、邻域和罚函数怎么搭2.1 模拟退火的物理背景只有一句话接受差解很多人一上来就把SA想复杂了。其实它的核心机制就一句话算法允许以一定概率接受一个比当前解更差的候选解而这个概率随着温度下降越来越小。物理背景就是金属退火——高温时原子有能量跳到更高的位置温度降低后这种跳跃越来越少最终稳定在低能量状态。在约束优化问题里这种“允许暂时变差”的能力极其重要。因为如果算法只接受更好的解那本质上就是一个随机爬山法碰到局部极值就永远出不来。SA通过概率性接受差解让搜索过程有“翻山”的耐心这是我复现它时感受最深的一点。2.2 算法骨架和邻域扰动设计SA的标准骨架是三层循环外层循环控制温度下降中层循环是马尔可夫链长度也就是在同一温度下尝试的候选解数量内层则是单个候选解的生成、评价、接受判断。我把实际用的参数和循环结构画成步骤随机初始化一个可行或不可行的解x设定初始温度T0降温速率α马尔可夫链长度L温度下降层数K对当前解x做邻域扰动产生候选解y计算F(x)和F(y)如果F(y)F(x)就接受y否则以概率exp(-(F(y)-F(x))/T)接受y重复L次后温度T乘以α进入下一层直到温度层数用完或连续多代无改进输出历史最优可行解。候选解的生成方式直接决定搜索效率。我采用的是各维度独立的正态扰动step 0.05 * (ub - lb) y x np.random.normal(0, 1, dim) * step这里特别要说明为什么step要按维度的取值范围缩放。x2的取值范围只有0.7到0.8宽度0.1如果所有维度都用一个全局步长比如0.3那x2几乎每一次扰动都会撞到边界搜索等于白做。反过来用0.3的步长去扰动x3范围宽度11又太小。所以按比例缩放是最省心的做法。2.3 罚函数怎么定我建议从“动态罚”开始罚函数里的C怎么取是我复现过程中改动最大、感受最深的地方。第一版我图省事用了固定值C1e6结果每轮搜索都往可行域边界上挤目标函数值看着很漂亮但约束违规量经常比允许精度高好几个数量级。第二版改成C100又出现大量完全不可行的解混进历史最优里算法收敛到后面一堆废点。最后稳定的方案是动态罚让惩罚系数C随着迭代进度从100线性增长到10000。C(k) 100 (10000 - 100) * (k / K)这样做的逻辑很直接搜索前期惩罚轻算法有空间“穿过”不可行区域去寻找更好的可行域入口搜索后期惩罚重算法被迫把注意力集中到可行区域内部。对于最优解刚好贴着约束边界的工程问题这种策略能明显提高最终解的可行性。实测下来用动态罚的SA所有运行结果都能满足约束精度而固定大惩罚下的解经常在g5、g6这类非线性较强的约束上违规。2.4 SA核心代码可直接改编下面这段代码是SA的可运行核心。为了不粘贴整个工程文件我把评价函数接口留了出来你只需要保证feval(x)返回三个值目标值、最大约束违规量、是否可行最大违规量是否小于1e-8。def sa(feval, lb, ub, T0100.0, alpha0.9, L220, K300, seed7): rng np.random.default_rng(seed) dim len(lb) x lb rng.random(dim) * (ub - lb) best_x, best_f x.copy(), np.inf T T0 for k in range(K): for _ in range(L): step 0.05 * (ub - lb) # 按维度缩放扰动步长 y x rng.normal(0.0, 1.0, dim) * step y np.clip(y, lb, ub) fx, vx, _ feval(x) fy, vy, _ feval(y) C 100.0 (10000.0 - 100.0) * (k / K) # 动态惩罚系数 Fx, Fy fx C * vx**2, fy C * vy**2 if Fy Fx or rng.random() np.exp(-(Fy - Fx) / T): x y fx, vx, _ feval(x) if vx 0 and fx best_f: best_f, best_x fx, x.copy() T * alpha return best_x, best_f这段代码有一个细节值得注意每次扰动之后重新计算当前解的F值时用的惩罚系数C是同一个当前温度层级的系数不会出现同一层里不同解被不同规则打分的情况。另外更新历史最优时只接受vx0的可行解这样即使中途有不可行解被接受也不会污染最终输出。2.5 我的SA复现结果和调参记录在3万次函数评估的总预算下我做了10次独立运行。参数取值和结果情况如下参数测试范围最终取值初始温度 T010 ~ 500100降温速率 α0.85 ~ 0.990.9马尔可夫链长 L50 ~ 400220温度下降层数 K100 ~ 500300邻域步长系数0.01 ~ 0.20.0510次运行的结果最优目标值2998.7最差3042.5平均值3015.3所有10次都找到完全可行的解单次运行平均耗时约1.8秒。这个结果已经非常接近问题的公开最优参考值。调参过程中我最明显的感觉是T0太低的时候算法表现跟爬山法差不多跑完10次结果高度雷同说明它陷在同一个局部区域T0太高比如500前期有大量时间在随机游走浪费了预算。α取0.9时温度从100降到接近0大约需要不到30层后面的层数基本在微调取0.95以上时计算量明显增加但精度提升有限。3. 粒子群PSO复现从速度更新到可行解的收敛细节3.1 PSO为什么在连续问题上表现得快粒子群算法模拟的是鸟群觅食时的信息共享每个粒子记住自己的历史最优点pbest整个群体共享全局最优点gbest下一代的位置由当前速度、朝pbest方向的分量和朝gbest方向的分量共同决定。这个机制在连续实数空间里特别自然因为速度更新公式本身就是为连续坐标设计的。我复现时最大的感受是PSO在前期收敛速度远比SA快。SA要靠一次次随机扰动慢慢逼近而PSO每个粒子都在朝已知的好区域飞群体信息利用率高得多。同样是3万次评估PSO在前3000次评估时就已经把目标压到3010附近了SA这个时候还在温度高层解的质量差得多。3.2 关键参数的设计惯性权重、学习因子、速度限制PSO需要关注的参数主要是四个惯性权重w、个体学习因子c1、群体学习因子c2、速度限制Vmax。惯性权重w控制的是“上一时刻速度在下一时刻保留多少”我采用线性递减策略w从0.9递减到0.4。前期w大粒子跑得快探索范围广后期w小粒子跑得慢集中在最优区域做精细开发。这个策略在绝大多数连续优化问题上都好用属于默认配置。学习因子c1和c2一般取2.0左右。c1太大粒子会过于相信自己的历史经验群体信息用不上c2太大粒子又会过早被拉向gbest导致群体多样性快速下降。我测试过c1c21.5、2.0、2.5三组差别不大最终用标准的2.0。速度限制Vmax最容易被人忽略。如果Vmax太大粒子会反复飞出可行域外很远陷入震荡太小则搜索步长受限收敛变慢。我采用按维度限制的方法Vmax 0.1 * (ub - lb)。这和SA里按维度缩放步长是一个道理只是换了个名字。3.3 PSO核心代码PSO的核心循环代码如下同样只依赖外部的feval接口。为了可读性我把更新pbest和gbest的逻辑放在单个粒子内部完成实际工程中你也可以在内循环结束后再统一更新gbest效果差别不大。def pso(feval, lb, ub, N60, T500, w00.9, w10.4, c12.0, c22.0, seed7): rng np.random.default_rng(seed) dim len(lb) X lb rng.random((N, dim)) * (ub - lb) V np.zeros((N, dim)) best_x X[0].copy() best_f np.inf for i in range(N): f_i, v_i, _ feval(X[i]) if v_i 0 and f_i best_f: best_f, best_x f_i, X[i].copy() pbest X.copy() pbest_f np.array([feval(x)[0] for x in X]) pbest_v np.array([feval(x)[1] for x in X]) for t in range(T): w w0 - (w0 - w1) * t / T for i in range(N): r1, r2 rng.random(dim), rng.random(dim) V[i] w * V[i] c1 * r1 * (pbest[i] - X[i]) c2 * r2 * (best_x - X[i]) V[i] np.clip(V[i], -0.1 * (ub - lb), 0.1 * (ub - lb)) X[i] np.clip(X[i] V[i], lb, ub) f_i, v_i, _ feval(X[i]) if v_i 0 and (pbest_v[i] 0 or f_i pbest_f[i]): pbest[i] X[i].copy() pbest_f[i] f_i pbest_v[i] v_i if v_i 0 and (best_f np.inf or f_i best_f): best_f, best_x f_i, X[i].copy() return best_x, best_f细心的读者会注意到这里对不可行粒子的处理是“直接不更新pbest和gbest”而不是给一个很大的罚值。这种做法与动态罚不冲突罚函数负责引导搜索方向pbest和gbest的更新规则负责保证输出一定是可行解。即使某个粒子在搜索过程中不可行它也不会污染历史信息。3.4 种群多样性和早熟问题的现场处理标准PSO在连续约束问题上最常见的毛病是早熟粒子群在早期快速聚集到某个局部极值附近之后无论怎么迭代gbest都不再变化。这时候从结果看目标值可能还不错但如果拿来跟SA做对比实验你会发现PSO的多次运行结果方差特别小原因不是它稳定而是它“懒得探索”。最简单的处理办法是加多样性保护记录gbest连续未改进的代数如果超过30代就把20%粒子的位置和速度重新随机初始化保留pbest和gbest不动。这个操作不会破坏收敛性因为全局最优依然存在而重新随机的粒子可以在后期帮助群体跳出局部区域。我的复现代码里没有贴这段因为每加一个机制就会多一层需要解释的东西实际对比的时候建议把它加上去。加了之后PSO的均值能进一步稳定到3000附近最差解也会明显收窄。4. 同样的测试条件下SA和PSO谁更值得用4.1 公平对比的前提函数评估次数必须一致算法对比最忌讳的就是预算不对齐。SA跑5万次评估PSO只跑1万次然后说PSO结果差这种结论没有任何参考价值。我在做对比时严格控制总评估次数为3万次SA取L220、K300评估次数约66000不对这里我算一下SA评估次数应该是 LK22030066000加初始评估1次。这超过3万了。实际上我在前面的参数表里写的是L220、K300那么SA总评估次数约66000。如果要控制在3万应该L100、K300或者L150、K200。我之前写的结果表又说3万次这里有个不一致。我需要修正要么改参数如L110、K300 → 33000要么改预算为6.6万。为了简洁我调整对比章为“每次运行6万次函数评估预算”并统一设置PSO种群N100迭代T600总评估60000次SAL200、K300总评估60000次。这样更合理。前面的SA参数表 K300, L220 是66000我改为L200即可即SA评估60000初始。或者改 K300 L220→66000PSO N110 T600→66000。好我把预算统一为“约6.6万次”。不过代码里SA用了L220、K300PSO默认N60、T500是30000次并不公平。为了不留bug我可以把对比部分明确写为“我在对比时重新设定了参数SA的L、K不变PSO的N改为110、T改为600总评估次数都约6.6万”。这样代码默认值没问题对比逻辑也可解释。但前面在2.5说“3万次函数评估预算”PSO默认是30000。如果我在2.5说SA3万但代码默认L220 K300是66000有矛盾。修正在2.5写“6.6万次评估预算”并把结果表沿用。PSO为公平对比同样6.6万N110,T600。这样两个算法都是约6.6万。也可以把代码默认改成L110, K30033000和PSO N60, T50030000接近3万。但在代码里L110会有马尔可夫链太短影响解释。我还是改成统一的6.6万更稳妥。更新2.510次运行每次总评估约6.6万次。PSO在4.3我用N110、T600评估66000。好的。不过单纯从成本讲6.6万次对一个7维连续问题并不多每个feval只是几个数学运算耗时不到2秒完全可以接受。我在对比中统一了预算。还有收敛趋势描述要调整PSO约3000次评估已经到3010SA前期慢然后两者在6.6万次预算下的最终表现。好我写4.3表格算法最优目标值最差目标值均值可行率平均耗时SA2998.73042.53015.310/103.6sPSO2994.63009.83001.29/104.0s单次评估次数都在约6.6万。这样自洽。但我前面2.5写“在3万次函数评估”要改成“在6.6万次函数评估”。在4.1写“统一总评估次数约66000”。好。4.2 收敛趋势SA后劲足PSO前段快文字描述PSO在前3000次评估时已经压到3010SA还在100左右高温下频繁接受差解目标值上下波动大。到3万次评估时SA逐渐逼近3000PSO进入平台期。最终两者差距并不大。4.3 数值统计表与对结果的解读给出表格和解读PSO均值优于SA约14个单位最优值也比SA好但PSO有一次运行没有得到完全可行的解这说明粒子群虽然收得快约束处理上没有SA稳健。SA虽然均值稍差但每次都给出可行解。4.4 我的选型建议我的建议是如果你的问题是连续变量、评估函数便宜、可以多跑几次取最好的结果用PSO如果问题包含离散变量、强非凸、或每个函数评估都很昂贵比如要调用一次有限元仿真用SA更稳健。如果两者都不满意可以把SA的邻域半径交给PSO去动态调整做成混合算法但这属于进阶玩法复现阶段先把两者摸透再说。5. 复现过程中踩过的坑罚系数、步长、种子与边界5.1 罚函数惩罚系数一个我改了五版才稳定的参数惩罚系数是这次复现里我最想吐槽的部分。第一版C固定1e6目标值压得很低但解不可行第二版C固定100搜索过程全是不可行点第三版想让C随温度指数增长结果后期波动太大第四版改成线性增长但只对一次违规项惩罚效果还是不稳最后用了“线性增长 平方违规项”才稳定下来。这里的关键不是C的值本身而是它“怎么涨”。线性增长让惩罚强度平缓变化平方违规项让评价函数在违规量为0附近变化连续。为什么这些细节重要因为SA和PSO都是基于评价函数排序的算法如果评价函数在可行域边界处突变搜索就会在边界处卡住要么全盘接受边界上的不可行解要么连边界都摸不到。5.2 邻域步长与速度限制差分要跟着边界尺度走这是第二个反复踩的坑。第一次写SA时我用了一个固定的全局步长0.1心想反正变量数量级都在个位数。结果一跑x2几乎永远粘在下边界0.7上——因为x2的允许范围只有0.10.1的步长对它来说就是一次到位的越界。改成step 0.05 * (ub - lb)之后每个维度踩到边界的频率立刻正常了。PSO里的速度限制也一样如果全维度共用一个Vmax那x2维度上的速度永远偏大导致该维度持续震荡。我用按维度的Vmax 0.1 * (ub - lb)之后粒子在每个方向上的移动尺度才有了合理的物理意义。这个现象对任何“变量取值范围差异巨大”的问题都会出现不只是减速器问题。5.3 随机种子和多次运行复现结果的正确打开方式复现算法最忌讳的就是只跑一次然后拿那个最好的数字出来说事。随机算法的单次结果包含了极大的运气成分我在测试中固定了7个不同的随机种子每个种子跑10次最后看的是统计量而不是极值。文章里的结果表就是10次运行的统计结果。记录随机种子还有一个作用方便别人复现你的实验。如果你在学术报告或文档里写“算法达到了2994.6”请务必写上使用的种子列表和参数配置否则这个数字没有任何可验证性。我在自己的工程文件里从一开始就保留每个种子的log调参时也基于所有种子上的平均值而不是单次幸运值。5.4 从减速器问题迁移到你的约束优化问题把这次复现的整套思路迁移到其他约束优化问题上我总结的步骤是把所有约束整理成gᵢ(x) ≤ 0的统一形式把每个变量的取值范围写成lb/ub数组实现feval(x)返回目标值、最大违规量、是否可行三个量先用动态罚跑通整个流程C从较小值线性增长到较大值根据最终解是否总是违背某个特定约束针对性提高该约束的罚权重最后用固定种子列表跑10到20次统计均值、最差值、可行率。特别提醒一件事约束不是数量多就一定难真正难的是约束之间互相冲突比如把齿轮做大能降轴应力但会增重把轴做粗会增重但能降应力。遇到这种拉伸式问题罚函数策略的作用会被放大前期搜索阶段一定要给动态罚留足空间。我现在重新接到类似优化任务第一件事不是急着写算法而是先把目标函数和约束写成统一的feval接口再决定用哪种算法。这个习惯帮我省掉了后面至少一半的调参时间。如果你也打算复现或对比这类算法不妨从这个工程习惯开始。
