简介这份基于MATLAB的模拟退火算法优化粒子群SA-PSO代码包面向需要改进群智能优化算法或进行目标函数寻优的本科及以上学习者适用于函数极值求解、参数优化等场景。压缩包共4个文件均为m脚本包含模拟退火算法主程序、粒子群初始化、适应度计算及迭代调用等模块代码注释完整便于二次扩展。已有362人学习资源体量仅2KB轻量易读适合快速掌握SA-PSO混合策略的核心思想。通过拆解代码框架读者既能理解模拟退火机制如何增强粒子群跳出局部最优的能力也能直接替换目标函数展开自己的优化实验省去从零搭建算法的时间。1. 把退火塞进粒子群到底在解决哪个痛点连续优化问题里粒子群PSO的早熟收敛几乎是每个用 MATLAB 做过优化的人都会撞上的墙。跑几十个测试函数算法前期收敛飞快到中后期速度向量逼近零整个种群挤在某个局部极小值附近再也没人肯往远处飞一步。模拟退火算法SA恰好擅长在这个阶段打破僵局——它在退火过程中以一定概率接受更差的解这种概率性跳跃是 PSO 缺失的机制。SA-PSO 是把两类算法的行为按时间尺度拆开再拼起来前期靠 PSO 的群体信息快速逼近有希望的区域后期靠 SA 的温度控制和 Metropolis 准则做局部逃逸。适合的对象很明确手头有 MATLAB 环境、跑过 PSO 但结果不理想、想在不重写整套寻优框架的前提下把全局搜索能力提一档的工程师。这篇不讲玄学直接给出能落地的 MATLAB 实现思路和参数设置经验。2. 粒子群原理和 SA 的数学互补性SA-PSO 为什么值得做2.1 粒子群的两条更新公式里藏着早熟收敛的病根标准 PSO 的核心只有两条公式。第 i 个粒子在第 d 维上的速度和位置更新如下v_id w * v_id c1 * r1 * (pbest_id - x_id) c2 * r2 * (gbest_d - x_d) x_id x_id v_id公式里的 w 是惯性权重c1 和 c2 是学习因子r1 和 r2 是 [0,1] 均匀分布的随机数pbest 是这个粒子自己历史最优位置gbest 是群体全局最优位置。速度更新由三项构成惯性项保留原来的飞行趋势认知项将粒子拉向自己的最佳记忆社会项将粒子拉向群体的最佳发现。问题恰恰出在这个结构上——如果 gbest 恰好落在一个局部极小值所有粒子的社会项都会指向这个局部极小值。经过若干代迭代后粒子的个体记忆 pbest 也会逐渐向 gbest 靠拢多样性随之丧失。最终所有粒子在局部极小值附近“冻结”速度趋于零算法不再有探索新区域的能力这就是教科书里说的早熟收敛。从控制论的角度看标准 PSO 是一个正反馈系统好的解吸引更多粒子靠近更多粒子靠近又强化了这个吸引。正反馈能加速收敛但也让系统失稳于局部最优。MATLAB 里做数值实验时可以观察到一种典型现象同样的代码跑同一个测试函数每次运行结果差异极大有时能收敛到全局最优附近有时 gbest 停在某个远离理论最优值的点上不动。这种高方差正是粒子群缺少负反馈机制、无法主动“拒绝劣质吸引子”的表现。2.2 模拟退火的 Metropolis 准则给算法装上概率逃逸阀模拟退火算法的数学基础是 Metropolis 准则。设当前解为 x_old适应度为 f_old新解为 x_new适应度为 f_new。如果 f_new f_old新解被无条件接受否则算法以概率P exp(-(f_new - f_old) / (k * T))接受这个更差的解。其中 k 为玻尔兹曼常数在工程实现中通常并入 T 或设为 1T 是当前温度。这个式子有两个关键性质第一f_new 与 f_old 的差值越大接受概率越低第二温度 T 越高接受概率越高。算法从高温开始以固定的冷却率逐渐降温最终在低温阶段只接受改善解收敛到一个稳定的最优解附近。SA 的“爬山能力”与 PSO 的“群体协同”是行为互补的。PSO 在迭代后期缺少的正是“有控制地变差”这种能力。SA 的退火过程本质上是时间轴上的探索策略高温期大肆探索低温期精细挖掘。把这个机制叠加到 PSO 上等于给粒子的速度更新增加了一条随机扰动通道——即使所有粒子的 pbest 和 gbest 都指向同一区域粒子仍有机会跳到远处且这种跳跃是概率性的、可控的而不是盲目的混沌扰动。这比简单加大速度扰动系数聪明得多扰动幅度会随温度自动收缩后期不会破坏已经找到的好解。2.3 SA 与 PSO 的结合方式混合粒度决定代码复杂度SA-PSO 在文献和工程实践中有几种不同的结合粒度它们的代码复杂度和行为差异很大。我按实际使用频率排序说明。第一种是“温度筛选”式也是本文后面重点给的实现。PSO 的粒子按常规速度公式更新位置后比较新旧适应度变好就接受变差就按 Metropolis 概率接受。温度按预设速率衰减。这种方式实现最简单只需在 PSO 主循环里插入两三行跟随机数比较的逻辑与原来的 PSO 代码融合度最高。第二种是“局部退火”式。确定当前 gbest 后在其邻域内以随机步长生成候选解用 Metropolis 准则决定是否替换 gbest。可以理解为在 PSO 外挂了一个局部搜索器。这种方式对代码结构改动小但每代多一次额外评估计算开销略有增加而且邻域步长的设定对效果影响很大。第三种是“级联”式。先用 PSO 跑完全程把最终 gbest 作为 SA 的初始解再单独做一次完整的退火搜索。这种做法的好处是两个算法互不干扰、参数调整独立缺点是 PSO 阶段早熟收敛产生的 gbest 质量直接决定了 SA 的起点上限而且耗时是两段叠加。选择哪种结合方式取决于你的目标函数评估成本。函数评估很快比如毫秒级可以直接用第一种简单直接函数评估很慢比如每次要跑几秒的仿真则用第二种在少数关键点上做局部退火更有性价比。下表给出三种方式的对比。结合方式代码改动量额外计算量对早熟收敛的改善适合场景温度筛选小极小明显常规连续优化函数评估成本低局部退火中每代 1 次评估强函数评估昂贵gbest 已接近真实解级联小翻倍依赖 PSO 输出思路清晰可分别调两个算法3. MATLAB 实现 SA-PSO从粒子结构体到主循环代码3.1 先定义粒子数据结构、目标函数和参数容器用 MATLAB 做这类优化我不建议把所有粒子数据摊开成零散变量来管理而是用结构体数组保存每个粒子状态。这样后续做并行评估、粒子数增减或状态可视化都比较方便代码可读性也好。下面是一个常见的初始化模板。% 目标函数Rastrigin用于验证算法的多峰寻优能力 fun (x) sum(x.^2 - 10*cos(2*pi*x) 10, 2); % 参数定义 NP 30; % 粒子数 D 10; % 问题维度 lb -5.12 * ones(1, D); % 变量下界 ub 5.12 * ones(1, D); % 变量上界 max_iter 1000; % 最大迭代次数 % SA 参数 T0 100; % 初始温度 T_end 1e-6; % 终止温度 alpha 0.90; % 温度衰减率 % PSO 参数 w_max 0.9; % 最大惯性权重 w_min 0.4; % 最小惯性权重 c1 1.5; % 个体学习因子 c2 1.5; % 群体学习因子 % 初始化粒子结构体数组 particles struct(pos, [], vel, [], fit, [], pbest, [], pbest_fit, []); for i 1:NP particles(i).pos lb (ub - lb) .* rand(1, D); particles(i).vel -0.1 * (ub - lb) .* rand(1, D); particles(i).fit fun(particles(i).pos); particles(i).pbest particles(i).pos; particles(i).pbest_fit particles(i).fit; end % 全局最优 [gbest_fit, best_idx] min([particles.pbest_fit]); gbest particles(best_idx).pbest;初始化时有几个细节值得注意。粒子速度不建议设为零向量否则前几次迭代会完全被认知项和社会项支配缺少自身的探索动量初期容易挤向某个随机点。速度量级设为变量范围的 10% 左右是比较保险的经验值。结构体数组里的 fit 和 pbest_fit 还必须分开存fit 是当前实际位置适应度pbest_fit 是这个粒子历史最优位置的适应度。在之后模拟退火接受劣解时粒子当前 fit 会变差而 pbest 不变这两个量必须由两个字段分别保存否则历史记忆会被污染。3.2 主循环在粒子群迭代中嵌 Metropolis 接受准则核心循环代码如下。每一代先更新惯性权重 w再逐粒子做速度位移更新随后用 Metropolis 准则决定是否接受变差解最后衰减温度。T T0; for iter 1:max_iter % 惯性权重线性递减 w w_max - (w_max - w_min) * (iter / max_iter); for i 1:NP % 标准粒子群速度更新 r1 rand(1, D); r2 rand(1, D); particles(i).vel w * particles(i).vel ... c1 * r1 .* (particles(i).pbest - particles(i).pos) ... c2 * r2 .* (gbest - particles(i).pos); % 位置更新与边界钳制 new_pos particles(i).pos particles(i).vel; new_pos min(max(new_pos, lb), ub); new_fit fun(new_pos); % 模拟退火接受准则 if new_fit particles(i).fit % 变好无条件接受 particles(i).pos new_pos; particles(i).fit new_fit; if new_fit particles(i).pbest_fit particles(i).pbest new_pos; particles(i).pbest_fit new_fit; end else % 变差按概率接受 delta new_fit - particles(i).fit; if exp(-delta / T) rand particles(i).pos new_pos; particles(i).fit new_fit; % 注意 pbest 不更新 end % 不接受则保持原位置 end end % 更新全局最优只依据 pbest_fit [gbest_fit, best_idx] min([particles.pbest_fit]); gbest particles(best_idx).pbest; % 温度衰减 T alpha * T; if T T_end T T_end; end % 可选输出每代最优值 fprintf(iter%d, gbest_fit%.4e, T%.4f\n, iter, gbest_fit, T); end这段代码里有几个设计决策必须解释清楚它们直接影响收敛行为。第一Metropolis 准则作用在粒子的当前位置上而不是全局最优上。粒子每次通过 PSO 公式生成一个新位置如果新位置比当前位置差并不立刻舍弃而是以 exp(-delta/T) 的概率接受它。这样粒子在温度较高时被允许暂时飞向不太好的区域增加种群空间分布多样性温度降低后接受概率指数级下降粒子群进入精细局部寻优。这个设计避免了把劣解写入 pbest 的危险——pbest 是粒子记忆里真正的好位置污染它会破坏整个搜索历史。第二全局最优 gbest 的更新只依赖 pbest_fit。即使某粒子在退火中接受了差解且 fit 指标变差其 pbest 仍然指向更优的历史位置所以全局最优不会因此被回退。这是防止“算法自我倒退”的关键约束。第三温度下限 T_end 的作用是让退火阶段自动终止。温度衰减到极低值后exp(-delta/T) 趋于零整个退火机制等同于关闭算法退化为纯 PSO 做最终收敛。这样 SA-PSO 既具备 SA 的前期探索能力又不影响 PSO 末期的收敛精度。fprintf 输出便于观察温度降速和最优值变化之间的对应关系。3.3 把测试函数换成真实目标函数时需要改哪里把模板接到自己的问题上时最常见的错误是只改 fun 函数就运行。以下三个地方必须同时检查。第一个是变量上下界 lb 和 ub。真实工程问题往往不是每个维度都有边界但 PSO 位置更新会产生越界值必须给出物理可行的边界范围。没有明确约束时建议从样本数据中取 min 和 max 作为先验边界而不是随意设一个很大的数——边界过宽会导致粒子在高维空间大片无效游荡。第二个是目标函数的方向。MATLAB 优化工具箱里的 fmincon 默认求极小值但如果你的业务指标是“吞吐量越大越好”“收益率越高越好”需要把目标函数写成负值。我见过不少人在这一步栽跟头算法收敛出来的“最优值”恰是业务意义下的最差值。第三个是评估函数的输入空间尺度。目标函数内部如果涉及单位换算、对数缩放等操作建议在 fun 开头统一处理。特别是当不同维度的物理量级差异极大时比如一个维度是温度几百量级、另一个维度是压力几兆帕量级需要在调用算法前做无量纲化归一化否则 PSO 的速度更新会被大量级维度支配小量级维度几乎没有探索力度。一个简单的做法是在外层套一个归一化 wrappernorm_fun (x_norm) fun(scaling .* x_norm offset);其中 scaling 和 offset 由实际边界换算得到。这样算法内部搜索空间始终规范在 [0,1] 范围内保持各维度探索力度均衡。4. SA-PSO 的关键参数设置温度曲线、惯性权重和混合时机4.1 一套可直接起步的参数表和经验范围SA-PSO 的参数数量比单一 PSO 至少多出三个初始温度 T0、终止温度 T_end、温度衰减率 alpha。新增参数的作用尺度与适应度函数数值量级强耦合因此没有能通吃所有问题的固定值。下面给出我常用的一套起步值适合维度 5 至 30、适应度量级在 1e-3 到 1e3 之间的连续优化问题。参数建议值调整方向作用说明粒子数 NP20 ~ 40复杂度高时可增至 60粒子数过多后期多样性冗余、计算浪费惯性权重 w0.9 线性降至 0.4前期偏探索、后期偏开发控制速度继承比例学习因子 c1, c21.5 左右c1c2 提升个体独立性c2 过大会加速聚集到 gbest初始温度 T0适应度典型差值的 2~5 倍过低则退火无效决定早期接受差解的概率上限终止温度 T_end1e-6 或 T0*1e-4过大会提前关闭退火低于此值退火机制等同于关闭衰减率 alpha0.85 ~ 0.95越大退火越慢、迭代越多主导退火探索持续时长最大迭代数1000 ~ 3000依函数复杂度调整与 alpha 配合保证温度降到位4.2 初始温度 T0 怎么定看适应度差异量级而不是拍脑袋初始温度是整个 SA-PSO 里最容易设错的参数。如果 T0 设得过大exp(-delta/T) 在早期几乎恒等于 1所有差解都会被接受粒子行为接近随机搜索前期收敛被严重拖慢。如果 T0 设得过小接受概率从一开始就很低SA 机制形同虚设混合算法退化为纯 PSO。合理的做法是根据初始粒子群的适应度分布来标定 T0。常见做法是让初始种群在退火开始时有大约 50% 到 80% 的概率接受一个“典型差解”。这里的典型差解可以用初始群体适应度的标准差来估计。一个快速启动办法是跑一次初始化统计所有粒子的适应度标准差然后取 T0 为标准差的 2 到 5 倍% 初始化后计算适应度标准差 fit_all arrayfun((s) s.fit, particles); fit_std std(fit_all); T0 3 * fit_std; % 按三倍标准差起步这种标定方式的直觉是初始群体的适应度离散程度直接反映了搜索空间的粗糙度。空间越粗糙不同随机点之间的适应度差异越大需要的初始温度越高。如果目标函数非常平坦标准差趋近于零说明任意两个随机点的适应度都差不多此时退火的核心不再是适应度差异而是空间位置差异可以考虑用位置扰动量来标定温度但大多数工程测试函数不会出现这种情况。4.3 退火持续范围与迭代次数的匹配SA-PSO 的温度衰减曲线与 PSO 迭代过程需要时间对齐。总迭代数 max_iter 固定时alpha 决定了退火在什么时间点降到 T_end。比如 T0100、T_end1e-6、alpha0.9温度降到 T_end 需要大约 log(1e-8)/log(0.9) ≈ 175 代如果 alpha0.95则需要约 360 代。如果 max_iter200那 alpha0.95 意味着整个迭代过程温度还没降到止点就已经结束退火在高温期被强行截断后期全部是随机扰动占主导。我一般会这样分配第一步固定 T0 和 T_end第二步根据期望退火在总迭代次数的前 40% 到 60% 完成即温度降到 T_end反推 alpha% 希望退火在前 40% 的迭代内完成 anneal_iter round(0.4 * max_iter); alpha exp((log(T_end) - log(T0)) / anneal_iter);这样做的理由比较实际迭代后期温度已经很低粒子群应当进入精细局部挖掘阶段此时退火机制已经无意义。继续让退火通行的劣解接受逻辑处于“几乎永远不接受”的状态除了白白增加计算开销外没有意义。反推 alpha 的做法让退火的探索强度分布与粒子群收敛阶段自动对齐不需要手动试算多组 alpha。4.4 惯性权重 w 的递减节奏要与退火探索强度互补如果说温度曲线控制的是“纵向跳跃”惯性权重递减控制的则是“横向飞行距离”。w 高时粒子飞得远w 低时飞得近。SA 的接受概率在高温期很高粒子本来就容易飞到远处如果此时 w 也很大粒子可能在一次迭代中飞出可接受范围之外且后续难以拉回反而降低搜索效率。因此两者应当按互补节奏设计前期 w 高但 T 高粒子经常飞远但也要靠边界钳制限制住中后期 w 降低粒子的探索主要靠残余的退火概率而非速度惯性。实践中我经常把 w 的递减函数从线性改为非线性比如余弦递减或指数递减配合温度衰减能产生更平滑的探索强度过渡。下面给出一个余弦递减的实现% 余弦递减惯性权重 w w_min 0.5 * (w_max - w_min) * (1 cos(pi * iter / max_iter));这段代码在早期保持接近 w_max 的探索强度在末尾平滑汇至 w_min与温度衰减的形态相似。两者都具“前期探索、后期开发”的同步趋势。如果算出来效果不稳定优先怀疑的是 w 的递减曲线与温度衰减率是否匹配而不是一上来就调 c1 和 c2。5. 用测试函数验证 SA-PSO多峰测试基准、统计性实验和状态同步排查5.1 选四个测试函数覆盖不同地形特征验证一个混合优化算法不能只跑一个 Sphere 函数。Sphere 是单峰平滑地形任何收敛性尚可的算法都能轻松解决显示不出 SA-PSO 的价值。至少选择四个地形特征差异明显的标准测试函数覆盖单峰、多峰强震荡、复杂山谷三种情况。函数名数学表达式最优解地形特征Spheresum(x_i^2)0单峰、平滑算法收敛速度测试Rosenbrocksum(100*(x(i1)-x(i)^2)^2 (x(i)-1)^2)0山谷狭窄弯曲测试可达精度Rastriginsum(x_i^2 - 10cos(2pi*x_i) 10)0多峰强震荡测试全局搜索能力Griewank1 sum(x_i^2)/4000 - prod(cos(x_i)/sqrt(i))0多峰且峰间嵌套混合难度高其中 Rastrigin 是最能体现 SA-PSO 相对 PSO 优势的函数。它有无穷多个局部极小值点且每个局部极小值的适应度相差不大标准 PSO 很容易陷入某个非最优的局部极小值。SA 的概率接受机制在高温期有条件地跳过这些局部极小值之间的势垒这正是对比实验中最容易出现明显差异的测试对象。5.2 跑 20 次、统计均值和标准差不比单次最优值验证算法的正确方式不是跑一次取最优而是多次运行取统计量。因为 PSO 和 SA-PSO 都是随机算法单次运行结果完全可能因为随机种子好坏而误导判断。下面给出一个简明的对比实验脚本框架用固定的参数配置分别运行标准 PSO 和 SA-PSO各跑 20 次% 对比实验脚本片段 n_runs 20; results_pso zeros(n_runs, 1); results_saps zeros(n_runs, 1); for r 1:n_runs % 每次运行前重置随机数保证两次算法在相同初始种子下对比 rng(r); results_pso(r) run_standard_pso(fun, lb, ub); rng(r); results_saps(r) run_sa_pso(fun, lb, ub); % 第二章程序封装成的函数 end % 统计输出 fprintf(PSO: mean%.4e, std%.4e, min%.4e\n, ... mean(results_pso), std(results_pso), min(results_pso)); fprintf(SA-PSO: mean%.4e, std%.4e, min%.4e\n, ... mean(results_saps), std(results_saps), min(results_saps));使用相同的 rng(r) 种子保证两个算法在每轮对比中从完全一致的初始种群出发这样消掉了初始种群随机性带来的干扰比较的是两种算法机制本身的差异。观察指标有三个均值代表平均表现标准差代表稳定性最小值代表最好运气下的上限。SA-PSO 在 Rastrigin 函数上的典型表现是均值显著低于标准 PSO同时标准差更小——这直接说明退火机制在减少对初始种子的依赖。5.3 一个容易被忽略的坑接受劣解后 pbest 与 fit 的状态同步前面代码里有一处非常容易写错的状态更新逻辑值得单独提出来强调粒子接受劣解时pos 字段和 fit 字段被更新但 pbest 和 pbest_fit 不能跟着更新。常见错误是把 pbest 直接写成 new_pos然后发现算法在某个迭代点突然出现“gbest 回退”的怪异现象。原因很直接pbest 保证这个粒子在整个搜索历史中记住的最好位置这是 PSO 认知项存在的根基。一旦劣解被接受粒子当前位置变差是暂时的是为了探索新区域付出的代价如果这个暂时的差位置被写进历史最优记忆下次迭代的速度更新会产生类似于“自我否定”的行为——认知项会把粒子拉向他刚刚故意放弃的差位置。而 gbest 的回退则会让所有粒子的社会项指向一个更差的方向形成系统性的倒退。调试这类问题可以在主循环每次更新 gbest 后加一个保护断言assert(islocal_comparable(gbest_fit, last_gbest_fit) 0, ... gbest 出现回退检查 pbest 是否被劣解污染);如果怀疑状态同步有问题打印每个粒子的 pbest_fit 和 fit逐代观察两者的差值变化。正常情况下每个粒子的 fit 可以高于 pbest_fit当前解比历史最优差但全局 gbest_fit 必须单调不增。这是一条直观、初步的算法健康度检查指标实验中发现违反这条约束时先检查 pbest 的更新分支是否有漏写 else 或错误赋值。本文还有配套的精品资源点击获取