企业官网建设流程全解析

简介：这份基于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 或设为 1），T 是当前温度。这个式子有两个关键性质：第一，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 的起点上限，而且耗时是两段叠加。

选择哪种结合方式，取决于你的目标函数评估成本。函数评估很快（比如毫秒级），可以直接用第一种，简单直接；函数评估很慢（比如每次要跑几秒的仿真），则用第二种，在少数关键点上做局部退火更有性价比。下表给出三种方式的对比。

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 的速度更新会被大量级维度支配，小量级维度几乎没有探索力度。一个简单的做法是在外层套一个归一化 wrapper：

norm_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 之间的连续优化问题。

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。比如 T0=100、T_end=1e-6、alpha=0.9，温度降到 T_end 需要大约 log(1e-8)/log(0.9) ≈ 175 代；如果 alpha=0.95，则需要约 360 代。如果 max_iter=200，那 alpha=0.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 的价值。至少选择四个地形特征差异明显的标准测试函数，覆盖单峰、多峰强震荡、复杂山谷三种情况。

其中 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 或错误赋值。

本文还有配套的精品资源，点击获取