简介面向MATLAB开发者的一套模拟退火-粒子群混合优化SAPSO算法实现资源用于求解复杂目标函数的最小值尤其适合处理多峰测试函数中传统PSO易陷入局部最优的问题。资源共10个文件包含7个m源文件和3个jpg图片m文件涵盖主程序、迭代过程、多个目标函数定义以及粒子群初始化、退火参数设置等完整模块图片直观呈现了算法收敛曲线与粒子分布变化。压缩包仅94KB轻量易读已有926人学习。通过这份资源读者既能掌握SA与PSO融合的思路也能理解温度控制、个体最优和全局最优更新等关键细节并能在Shubert等多峰函数上直接验证算法效果为后续优化算法改进或学术实验提供可复现的样例。尤其适合需要设计全局优化算法或进行算法对比研究的学生、教师与工程师使用。1. 粒子群算法早熟收敛时模拟退火算法怎样把目标函数最小值求解从局部最优中拉出来跑过 Rastrigin、Ackley 这类多峰测试函数的人大概率见过粒子群算法的翻车现场前 100 代下降很快之后几百代几乎横盘最终停在某个并非全局最小的区域。问题不在迭代次数而是粒子群算法原理里最要命的一点——所有粒子共享同一份全局最优信息一旦这份信息来自局部最优种群多样性就迅速归零再跑多少代都只是原地抖动。模拟退火算法恰恰是用温度控制对劣化解的接受概率高温阶段允许往“更差”的方向探索低温阶段收紧到精细搜索。把 Metropolis 接受准则接到粒子群算法的全局最优更新上就构成模拟退火算法优化粒子群算法SA-PSO。这篇文章按原理、最小实现、调参、验证的顺序回答两个问题模拟退火算法优化粒子群算法求解目标函数最小值怎么做以及如何用测试函数证明这种优化确实有效。2. 模拟退火算法优化粒子群算法的核心原理Metropolis规律与温度衰减的接入点2.1 粒子群算法原理中的共享信息机制为什么在多峰测试函数上必然早熟先看粒子群算法原理里决定收敛节奏的三个部分惯性项w * v保持原运动方向个体认知项c1 * r1 * (pbest - x)把粒子拉向自己的历史最优社会认知项c2 * r2 * (gbest - x)把粒子拉向全局最优。前 100 代里pbest 分布在解空间的多个区域粒子被不同的历史最优拉扯种群覆盖范围还比较大适应度下降明显。随着迭代推进粒子的个体最优逐渐汇集到同一个盆地pbest 之间的差异越来越小速度更新项里两个差距向量的模长趋近于零粒子几乎失去移动能力。这时候如果把种群适应度的方差打出来看往往在第 200 代就降到接近机器精度。表面看是“收敛”实质是失去多样性换成Griewank、Rastrigin这类具有大量规则局部最小值的测试函数粒子群几乎必然停在某一个对称的坑里。想让粒子重新具备翻越能力最直接的做法是在全局最优更新环节引入一个可以接受劣化解的概率门控而这个门控正是模拟退火算法最成熟的 Metropolis 接受准则。2.2 Metropolis 接受准则接入粒子群 gbest 更新的两种常见做法模拟退火算法原理中当前状态s_i产生候选状态s_new后计算能量差delta E(s_new) - E(s_i)。若delta 0直接接受新状态若delta 0以概率exp(-delta / T)接受其中温度T随时间下降。把这个机制映射到粒子群算法上第一类做法是改 gbest 的更新逻辑把“当前代最好的粒子”当作候选状态按 Metropolis 规则决定是否替换 gbest。# Metropolis 接受准则的最小骨架应用于粒子群算法 gbest 更新 delta candidate_fitness - gbest_fitness if delta 0: gbest candidate else: # exp(-delta/T) 在 T 大时接近 1.0在 T 小时接近 0 threshold np.exp(np.clip(-delta / T, -745, 0)) if np.random.rand() threshold: gbest candidate这段代码里np.clip(-delta / T, -745, 0)是为了防止delta / T过大导致exp下溢告警当温度T 1000而delta 1时接受概率是千分之九百九十九几乎来者不拒当温度降到T 0.1时同样 delta 的接受概率只有e^{-10}基本退化成标准 PSO 的“只接受更好解”。这就是模拟退火算法优化粒子群算法最核心的平衡高温阶段用劣化解保多样性低温阶段把选择压力还给适应度本身。第二类做法是对 gbest 本身做邻域扰动。每一代给当前 gbest 叠加一个与温度相关的随机向量得到一个候选解再用 Metropolis 准则判断是否采用。若采用就把这个候选解放入种群替换最差粒子从而在不破坏全局最优记忆的前提下增加种群扰动。实践中我一般会用第一类做法保证全局最优更新的合理性用第二类做法作为低温段的局部精修两者互补而不冲突。2.3 温度衰减与接受概率的配套关系T0、alpha、T_end 的决定性影响把 Metropolis 接入粒子群后整个 SA-PSO 的状态变量比基础 PSO 多出三个核心量初始温度T0、温度衰减系数alpha、终止温度T_end。它们的作用可以先用一张表定住边界状态变量常用取值在混合算法中的职责T0初始温度100 到 5000决定早期对劣化解的包容程度越大探索越强alpha衰减系数0.90 到 0.99决定温度下降快慢越接近 1.0 退火越慢T_end终止温度1e-5 到 1e-2防止除零也作为算法进入纯精细搜索的开关w惯性权重0.4 到 0.9控制粒子继承上一代速度的比例c1/c2加速系数1.2 到 2.0控制个体经验与社会经验对速度的影响温度衰减通常采用几何退火T max(T_end, T * alpha)。几何退火的好处是参数少、行为稳定alpha 0.95时温度从 1000 降到 1 需要约 135 次迭代alpha 0.99则需要约 690 次。前者适合迭代预算有限的场景后者适合求解精度要求高且允许跑 2000 代以上的测试函数。要注意delta与T的量级必须匹配T0设为 1000 而目标函数值是1e-6量级时exp(-delta/T)几乎始终等于 1SA 阶段会退化成均匀随机游走反过来T0 0.1而函数值普遍在1e2量级Metropolis 规则永远不会接受任何劣化解混合算法就退化为普通 PSO。这个匹配问题是调参失败时首先要排查的对象。3. 用 Python 实现 SA-PSO 求解目标函数最小值的最小可运行代码3.1 三个测试函数的定义Sphere、Rastrigin、Ackley 的形态差异测试函数的选择决定验证结果的说服力。Sphere 是单峰函数用来确认算法没有把简单问题搞坏Rastrigin 引入了大量规则排列的局部极小值能暴露粒子的早熟收敛Ackley 在原点附近有一个窄的全局最小区外围是起伏的丘陵对跳出局部最优的能力非常敏感。下表直接给出三个函数在 30 维下的标准测试配置函数名表达式搜索范围全局最优值Spheref(x) sum(x_i^2)[-100, 100]^D0Rastriginf(x) 10D sum(x_i^2 - 10cos(2πx_i))[-5.12, 5.12]^D0Ackleyf(x) -20exp(-0.2*sqrt(sum(x^2)/D)) - exp(sum(cos(2πx))/D) 20 e[-32, 32]^D0对应代码可以直接放进同一个模块后面对照实验时通过参数传入目标函数。import numpy as np def sphere(x): return np.sum(x ** 2) def rastrigin(x): return 10 * len(x) np.sum(x ** 2 - 10 * np.cos(2 * np.pi * x)) def ackley(x): a, b, c 20, 0.2, 2 * np.pi n len(x) sum1 np.sum(x ** 2) / n sum2 np.sum(np.cos(c * x)) / n return -a * np.exp(-b * np.sqrt(sum1)) - np.exp(sum2) a np.e这段代码里np.e对应欧拉常数ackley函数末尾的 a np.e把整个函数的全局最小值规整到 0。三个函数都以0为全局最优值方便后期比较收敛曲线的绝对精度。3.2 SAPSO 类实现速度更新、位置更新、Metropolis 接受与退火调度下面这个类把粒子群算法和模拟退火算法写在一个循环里去掉绘图和日志保留求解目标函数最小值的最小逻辑。class SAPSO: def __init__(self, func, dim30, pop_size40, max_iter1000, lb-5.12, ub5.12, T01000, T_end1e-5, alpha0.985, w0.72, c11.5, c21.5, seed1): self.func func self.dim dim self.pop_size pop_size self.max_iter max_iter self.lb, self.ub lb, ub self.T0, self.T, self.T_end, self.alpha T0, T0, T_end, alpha self.w, self.c1, self.c2 w, c1, c2 self.rng np.random.default_rng(seed) v_max (ub - lb) * 0.2 self.v_max v_max self.positions self.rng.uniform(lb, ub, (pop_size, dim)) self.velocities self.rng.uniform(-v_max, v_max, (pop_size, dim)) self.pbest self.positions.copy() self.pbest_fitness np.array([func(p) for p in self.positions]) best_idx np.argmin(self.pbest_fitness) self.gbest self.pbest[best_idx].copy() self.gbest_fitness self.pbest_fitness[best_idx] # best_ever 记录历史最优gbest 是允许被 Metropolis 暂时拉高的工作解 self.best_ever_pos self.gbest.copy() self.best_ever_fitness self.gbest_fitness self.best_history [] def optimize(self): for _ in range(self.max_iter): r1 self.rng.random((self.pop_size, self.dim)) r2 self.rng.random((self.pop_size, self.dim)) self.velocities (self.w * self.velocities self.c1 * r1 * (self.pbest - self.positions) self.c2 * r2 * (self.gbest - self.positions)) self.velocities np.clip(self.velocities, -self.v_max, self.v_max) self.positions np.clip(self.positions self.velocities, self.lb, self.ub) fitness np.array([self.func(p) for p in self.positions]) better fitness self.pbest_fitness self.pbest_fitness[better] fitness[better] self.pbest[better] self.positions[better] cur_idx np.argmin(self.pbest_fitness) cand_fitness self.pbest_fitness[cur_idx] delta cand_fitness - self.gbest_fitness # SA 核心入口delta 大于 0 时按温度接受劣化解 if delta 0 or self.rng.random() np.exp( np.clip(-delta / self.T, -745, 0)): self.gbest self.pbest[cur_idx].copy() self.gbest_fitness cand_fitness if self.gbest_fitness self.best_ever_fitness: self.best_ever_pos self.gbest.copy() self.best_ever_fitness self.gbest_fitness self.best_history.append(self.best_ever_fitness) self.T max(self.T_end, self.T * self.alpha) return self.best_ever_pos, self.best_ever_fitness, self.best_historyoptimize方法里每一代的执行顺序固定为五步更新速度、更新位置、计算适应度并刷新 pbest、对候选 gbest 做 Metropolis 接受、温度衰减。速度更新公式与标准粒子群算法原理完全一致唯一的不同是把gbest换成了可能被劣化解临时占据的工作最优解。best_ever_fitness独立于gbest_fitness因此即使 Metropolis 接受了更差解返回的仍然是一整轮实验里真正出现过的历史最优值。调用方式如下opt SAPSO(rastrigin, dim30, lb-5.12, ub5.12, T01000, alpha0.985, max_iter1000, seed42) best_pos, best_val, history opt.optimize() print(best_val, best_pos[:3])如果使用plt.plot(history)绘制收敛曲线会看到与普通 PSO 明显不同的形态曲线在高温段有小幅上涨这是 Metropolis 接受劣化解造成的随后在低温段快速下降最终收敛到比纯 PSO 更低的精度。3.3 用同一套 seed 对比 PSO 与 SA-PSO排除随机性干扰验证混合优化是否有效时常见做法是固定随机种子把 SA-PSO 与去掉 Metropolis 接受的普通 PSO 各跑一遍。去掉 SA 最简单的方式是把T0设置成极小值比如1e-10此时exp(-delta/T)几乎为 0接受准则退化为“只接受更好解”。将T01e-10的 SAPSO 看作 PSO 的近似等价实现在 Rastrigin 30 维上各跑 20 次随机种子记录平均最优值。两者使用完全相同的初始位置和速度更新逻辑差异只来自 Metropolis 接受是否生效这样得到的对比结果才不会被随机种子差异污染。4. SA-PSO 求解测试函数的参数调优实战温度与粒子群超参的协同配置4.1 退火温度初值 T0、衰减系数 alpha 与终止温度 T_end 的选择边界在求解目标函数最小值时T0的量级取决于目标函数适应度的典型尺度。对 Rastrigin 30 维随机初始解的适应度通常在几百到上千T0设 1000 到 5000 比较合理对 Sphere 30 维初始适应度可高达几十万但局部最优之间的落差远大于 RastriginT0反而可以设得更高。经验法则是先随机采样 100 个解统计适应度标准差T0取该标准差的 1 到 5 倍。alpha的选择与max_iter强相关。若迭代上限是 1000 代alpha0.985可以让温度经历完整的从高温到低温的过程若alpha0.99温度下降到原温度的十分之一需要约 230 代后面大部分时间都处于低温精细搜索若alpha0.90温度几乎在第 70 代就降到原来的千分之一以下SA 的探索能力过早丧失。T_end不必设置过小1e-5已经足够。因为一旦delta/T超过 745exp(-delta/T)就是数值零继续降低温度只是徒增计算开销。4.2 惯性权重 w 与接受概率的联动高低温阶段的节奏控制固定w0.72是最常见的取值但在 SA-PSO 中w与温度下降曲线需要协同。高温阶段粒子的速度本身比较大如果w也设置得过高粒子会出现大幅震荡粒子群算法原本靠社会认知项收敛的节奏会被打乱低温阶段粒子需要精细化搜索较大的w又会让粒子难以停在窄小的最优盆地。我一般会用线性递减惯性权重把w从 0.9 降到 0.3与温度衰减保持同一个节奏前 200 代高温配合大w增强探索后 800 代低温配合小w把速度控制在较小范围。对应的修改是在optimize方法内把固定self.w改为w_current 0.9 - 0.6 * iter / self.max_iter其余逻辑不变。这个改动对 Rastrigin 的影响往往比调alpha更直接原因是 Rastrigin 局部最小值排列密集低温阶段失败只要粒子还在运动就很容易撞进旁边另一个坑。4.3 12 组 T0-alpha 对照实验脚本与结果判读方法固定的参数组合容易掩盖问题我用一个网格脚本一次性跑 12 组组合每组重复 20 次随机种子输出均值和中位数。scores {} configs [] for T0 in [100, 1000, 5000]: for alpha in [0.90, 0.97, 0.99]: vals [] for seed in range(20): opt SAPSO(rastrigin, dim30, T0T0, alphaalpha, seedseed, max_iter1000) _, best_val, _ opt.optimize() vals.append(best_val) mean_val np.mean(vals) median_val np.median(vals) scores[(T0, alpha)] (mean_val, median_val) print(fT0{T0:5d} alpha{alpha:.2f} fmean{mean_val:.4e} median{median_val:.4e})这个脚本的价值不在数值本身而在判读方式。先看中位数它比均值更抗重尾干扰若中位数接近 0 而均值明显偏大说明少数种子陷入极差的局部最优SA 大概率在高温阶段接受了太多劣化解。再看alpha0.90和alpha0.99的差距如果前者明显更差说明退火过快如果两者差异不大说明瓶颈在粒子群算法自身参数而不是温度调度。把 12 组结果按(median, mean, min)排序优先选 median 最小且 mean 与 median 差距小的一组而不是选单次运行最优值最好的一组。下表给出我按照调试经验整理的参数行为对照供先跑通再细调T0alpha典型行为特征优先适用的测试函数1000.90退火极快退化接近 PSO早熟明显Sphere1000.97前期探索不足收敛偏慢Ackley1000.99低温段过长后期几乎无 SA 作用Rastrigin10000.90前 100 代接受劣解过频稳定性差GridWank10000.97探索与开发均衡常作为默认起点Rastrigin10000.99适合 2000 代以上长迭代Ackley, Rastrigin50000.90高温段搜索粗放位置更新震荡大Sphere50000.97探索充分但收敛速度慢Ackley50000.99需要充足迭代预算低温段表现稳定Griewank实际调试时把 20 次重复改成 10 次快速跑肉眼观察曲线就能判断出问题是出在高温过跳还是低温困住不需要每次都跑完 1000 代。5. 验证 SA-PSO 在测试函数上的改进不是运气的三个统计技巧5.1 重复 30 次实验把单点值换成成对分布的秩检验SA-PSO 引入随机接受后单次运行结果方差通常大于标准 PSO只看一次最优值很容易得出错误结论。验证改进时我会用相同随机种子分别跑标准 PSO 和 SA-PSO 各 30 次然后用 Wilcoxon 符号秩检验判断两组结果是否来自同一分布。from scipy.stats import wilcoxon pso_vals [] sapso_vals [] for seed in range(30): opt_pso SAPSO(rastrigin, T01e-10, seedseed) _, v_pso, _ opt_pso.optimize() pso_vals.append(v_pso) opt_sa SAPSO(rastrigin, T01000, alpha0.985, seedseed) _, v_sa, _ opt_sa.optimize() sapso_vals.append(v_sa) stat, p_value wilcoxon(sapso_vals, pso_vals) print(p_value)注意pso_vals和sapso_vals必须来自相同 seed 下的同一初始化否则配对检验的前提不成立。若p_value 0.05可以认为 SA-PSO 的改进显著若p_value 0.2再调参不如先检查T0与函数值量级是否匹配。5.2 同时看收敛曲线的“中断面积”而不仅是终值终值只反映最后一代的结果无法区分算法是第 10 代就找到好解还是第 990 代才碰巧找到。把两条收敛曲线的纵轴取对数后在每代取较差值与较优值的差值再对所有代求和得到中断面积。标准 PSO 往往曲线平滑但居高不下SA-PSO 会出现短暂抬升后再下降抬升次数和幅度直接对应逃逸局部最优的次数。若抬升后很快回到更低的水平说明 Metropolis 接受正在有效工作若抬升后长时间不能恢复则说明接受概率过大需要降低T0或调大alpha。5.3 一个提高低温段精度的实用技巧让 gbest 扰动幅度随温度收缩在optimize循环末尾加入一段对 gbest 邻域的局部搜索可以让低温阶段不再空转。做法是用当前温度与初始温度的比值作为扰动半径的比例因子候选解从 gbest 附近的高斯邻域中生成if self.T self.T0 * 0.05: step (self.ub - self.lb) * 0.001 candidate self.gbest self.rng.normal(0, step, self.dim) f_cand self.func(candidate) if f_cand self.best_ever_fitness: self.best_ever_pos candidate.copy() self.best_ever_fitness f_cand这段代码放在温度衰减之前相当于在低温段用 SA 的邻域搜索补充 PSO 的粗粒度更新。扰动半径固定为搜索范围的千分之一不随温度继续缩小避免数值精度限制导致粒子长时间无位移。对于 Rastrigin 这类最优值周围地形陡峭的函数这个技巧通常能把最终结果提高一个数量级而且不需要额外增加迭代次数。本文还有配套的精品资源点击获取