SAPSO算法MATLAB实现：融合模拟退火解决粒子群局部最优

简介一套基于MATLAB实现的模拟退火算法优化粒子群算法SAPSO源码包面向需要求解多峰函数全局最小值的算法学习者与研究人员。资源以10个文件打包压缩包仅94KB主要包含7个.m源文件与3张结果图sapso.m为主入口iterateSAPSO.m提供SA-PSO迭代流程shubertfun.m、fun2.m、fun3.m等对应不同测试函数jpg图片直观展示目标函数变化曲线或粒子分布便于对照代码理解算法运行过程。已有926人学习下载。资料从SA与PSO的基本原理出发完整覆盖参数初始化、适应度计算、粒子与速度更新、模拟退火温度控制、全局最优解输出等关键模块尤其适合处理易陷入局部最优的复杂测试函数场景。通过阅读与运行这些代码可快速掌握SAPSO的耦合思路并在此基础上调整目标函数或参数迁移至实际工程优化问题。1. SAPSO为什么单纯PSO容易卡在局部最优跑过粒子群算法的人应该都有这种体验连续几个测试函数都能收敛到不错的解但一旦换成Shubert这类极值点分布密集的多峰函数PSO经常在迭代中期就“抱团”停在一个局部深谷里gbest不再变化粒子速度趋近于零。这是因为PSO的收敛机制本质上是向群体历史最优位置靠拢一旦gbest落在某个局部极小附近整个种群的搜索方向就被锁定。模拟退火的思路恰好补上了这块短板——它以一定概率接受比当前更差的解在温度较高时允许粒子“跳出”局部最优随着温度降低再逐步收敛。本文基于一套MATLAB实现的SAPSO源码包sapso.m、iterateSAPSO.m以及多个测试函数文件从算法机制、代码逐模块拆解到Shubert函数上的实测调参把这套组合优化策略讲清楚。适合正在做智能优化算法对比实验、需要求解多峰函数最小值或者想在PSO基础上加全局探索能力的读者。2. 模拟退火与粒子群融合的机制拆解2.1 PSO的位置-速度更新与早熟成因标准PSO的粒子i在第t次迭代时速度与位置更新公式如下v_i(t1) w * v_i(t) c1 * r1 * (pbest_i - x_i(t)) c2 * r2 * (gbest - x_i(t)); x_i(t1) x_i(t) v_i(t1);w是惯性权重c1和c2分别为个体学习因子和社会学习因子r1、r2是[0,1]上的均匀随机数。代码中会先评估粒子当前适应度然后更新pbest、gbest再进入下一轮迭代。参数比较常见的取法是w从0.9线性递减到0.4c1 c2 2.0种群规模取30到50。问题在于gbest一旦被某个较优的局部解占据其他粒子会被持续拉向该位置即使某几个粒子发现了更远的低谷方向也会因为飞行速度被w和c2压制而难以到达。Shubert函数在[-10,10]的定义域内有大量深度相近的局部极小点这种地形对PSO非常不友好。从群体行为角度看PSO缺乏“撤退”机制。每个粒子的运动只受个人经验和群体经验牵引没有显式的回退选项。模拟退火则赋予了搜索过程接受劣质解的能力这种能力在算法早期尤其重要因为它让种群保持探索多样性而不是过早收缩到某个吸引域。2.2 Metropolis准则与退火温度曲线模拟退火的基础是Metropolis接受准则。设当前解为x扰动后的新解为x_new对应的适应度值为f(x)和f(x_new)。对于最小化问题新解的接受概率为delta f(x_new) - f(x); if delta 0 x x_new; else p exp(-delta / T); if rand p x x_new; end end温度T是核心控制变量。T越大exp(-delta/T)越接近1算法越容忍差解T趋近于0时接受概率趋近于0算法退化为爬山法。常见的降温策略有三种线性降温T T0 * (1 - k / K)指数降温T T0 * 0.9^k以及自适应降温。指数降温在工程中用得最多因为参数少、实现简单而且前期降温快、后期温度变化放缓正好匹配“先探索后利用”的搜索节奏。2.3 SAPSO的融合方式与参数组合SA和PSO的融合主要有两种实现路线。第一种是将Metropolis准则嵌入PSO迭代粒子飞行后产生候选解用SA概率决定是否接受该位置相当于给每个粒子的移动加了一个“允许变差”的门槛。第二种是把SA作为对gbest的扰动机制在每轮迭代结束时对全局最优解做随机扰动如果扰动后的解更优则替换gbest。这套MATLAB代码从文件分布上看sapso.m是主程序iterateSAPSO.m负责迭代循环融合方式采用第一种路线即在粒子更新后对每个粒子执行SA接受判断。融合后需要配套调整的参数主要有参数典型取值范围作用初始温度T0100 ~ 1000控制初期对差解的容忍度降温系数alpha0.85 ~ 0.98控制温度下降速率温度迭代步数每代或每若干代降温一次决定SA节奏和PSO代数的配合比例PSO种群规模20 ~ 50兼顾客群多样性和计算开销SA扰动步长0.01 ~ 0.1倍定义域范围决定了新解的空间跨度初始温度如果设置过低Metropolis准则在迭代早期就几乎不接受差解SA的全局探索能力被浪费如果太高则前几十代粒子几乎没有收敛趋势浪费计算资源。实际调试时我一般先固定alpha 0.95再根据目标函数极值波动范围倒推T0让初始接受概率维持在0.5到0.8之间。3. MATLAB代码逐模块拆解sapso.m、iterateSAPSO.m与测试函数文件3.1 sapso.m主函数参数初始化与文件组织逻辑直接打开sapso.m会发现主程序的结构非常直观大致分四段定义目标函数句柄、初始化PSO参数、调用迭代函数、输出结果。下面是核心初始化代码段clear; clc; nPop 30; % 粒子群规模 maxIter 500; % 最大迭代次数 c1 2.0; c2 2.0; % 学习因子 wmax 0.9; wmin 0.4; % 惯性权重范围 T0 500; % 模拟退火初始温度 alpha 0.95; % 温度衰减系数 dim 2; % 目标函数维度 lb -10 * ones(1, dim); % 变量下界 ub 10 * ones(1, dim); % 变量上界 fun shubertfun; % 目标函数句柄可替换为fun2/fun3等 [x_best, f_min] iterateSAPSO(fun, nPop, maxIter, c1, c2, ... wmax, wmin, T0, alpha, lb, ub); fprintf(最优解: [%.6f, %.6f]\n, x_best(1), x_best(2)); fprintf(目标函数最小值: %.6f\n, f_min);注意这里的fun是函数句柄通过符号指向shubertfun.m或fun2.m等文件这样sapso.m不需要改动即可切换测试函数。lb和ub是决策变量的边界约束PSO粒子初始位置会在边界内随机分布。惯性权重w没有直接写成常数而是按迭代进度线性递减这是在主函数中通过wmax、wmin和当前迭代次数计算出来的。这种做法在单峰函数上加速收敛效果明显但在多峰地形上过快降到wmin会削弱种群多样性配合SA的接受机制可以适当缓解。3.2 iterateSAPSO.m迭代核心温度更新与新解接受逻辑iterateSAPSO.m是整套代码的发动机内部完成粒子初始化、适应度计算、速度与位置更新、SA接受判断、温度衰减、结果记录六个动作。核心迭代逻辑如下for iter 1:maxIter % 线性递减惯性权重 w wmax - (wmax - wmin) * iter / maxIter; for i 1:nPop % 更新速度与位置 v(i, :) w * v(i, :) c1 * rand(1, dim) .* (pbest(i, :) - x(i, :)) ... c2 * rand(1, dim) .* (gbest(1, :) - x(i, :)); x(i, :) x(i, :) v(i, :); % 边界约束处理 x(i, :) max(x(i, :), lb); x(i, :) min(x(i, :), ub); % 计算新适应度并应用SA准则 newFit fun(x(i, :)); delta newFit - fit(i); if delta 0 x(i, :) x(i, :); % 新解更优直接接受 fit(i) newFit; else pAccept exp(-delta / T); if rand pAccept x(i, :) x(i, :); fit(i) newFit; else % 拒绝差解粒子回退到上一个位置 x(i, :) x(i, :) - v(i, :); end end % 更新pbest与gbest if fit(i) pbest_fit(i) pbest_fit(i) fit(i); pbest(i, :) x(i, :); end if fit(i) f_min f_min fit(i); gbest(1, :) x(i, :); end end % 温度衰减 T alpha * T; record(iter) f_min; end这段代码有几个关键参数需要理解。delta是目标函数的变化量如果新解更优则delta为负无条件接受如果新解更差delta为正用exp(-delta/T)计算接受概率同时用rand生成均匀随机数来决策是否接受。注意温度下降放到了每代结束后所以粒子飞行和SA判断共用同一个T粒子的探索能力随代数和温度同步下降。边界约束处理用的max和min是MATLAB的逐元素比较函数防止粒子飞出定义域。拒绝差解时粒子通过减去当前速度v(i, :)回退到上一轮的位置。这里的回退操作与标准PSO不同标准PSO不做任何拒绝判断粒子差的位置也会保留并参与后续演化SAPSO保留了SA的判据被拒绝的粒子位置不发生改变相当于本次飞行被取消。3.3 funv.m、funx.m、fun2.m、fun3.m与shubertfun.m的职责划分代码包中多个fun开头的文件容易让人迷惑实际分为两类。funv.m和funx.m是PSO迭代过程中固定使用的辅助函数funv.m用于计算粒子在约束边界处的额外速度约束funx.m用于边界越界时的坐标映射而fun2.m、fun3.m、shubertfun.m则分别是不同的目标测试函数可以独立替换到sapso.m的fun句柄位置。文件类型作用funv.m辅助函数计算粒子速度越界时受到的边界反弹或截断逻辑funx.m辅助函数计算粒子位置越界后的投影修正值fun2.m测试函数通常为二维球形函数或Rosenbrock类单峰函数fun3.m测试函数通常为Rastrigin类周期多峰函数shubertfun.m测试函数Shubert函数测试全局优化能力的关键用例以shubertfun.m为例Shubert函数的表达式为f(x_1, x_2) Σ_{i1}^{5} i * cos((i1) * x_1 i) * Σ_{i1}^{5} i * cos((i1) * x_2 i)该函数在[-10, 10] × [-10, 10]内有大量局部极小点和多个全局极小点是检验全局优化算法跳出局部最优能力的经典标准测试函数。把测试函数从shubertfun换成fun2或fun3只需要改动一行函数句柄整套迭代逻辑不用动这也是这套代码设计上比较实用的地方。4. 在Shubert多峰函数上验证SAPSO的有效性4.1 Shubert函数的极值分布与测试意义Shubert函数之所以难优化是因为它的极小值分布非常密集。在定义域内存在许多深度相近的局部极小点且这些极小点之间的距离小于PSO粒子速度的常规尺度导致粒子群很容易在某个局部极小附近“误以为”已经找到了全局最优。在二维Shubert函数中全局最小值的近似值为-186.731对应的最优点有多组。这意味着算法只给出目标函数最小值还不够还要看是否稳定到达-186.7附近。如果运行SAPSO后f_min落在-140到-180区间说明算法找到了某个较优的局部极值但未命中全局只有多次运行都能稳定收敛到-186.73才算真正验证了SA对PSO的优化效果。单次实验不构成结论建议每组参数至少跑20次并统计最小值、平均值和方差。4.2 运行方式与对比实验设计直接运行sapso.m即可看到输出结果但要想验证SA的贡献需要同时跑一个标准PSO作为对照。可以在MATLAB命令行中运行以下脚本对比两种算法在相同预算下的表现% 对照组标准PSO rng(1); [pso_best, pso_record] iteratePSO(shubertfun, 30, 500, 2, 2, 0.9, 0.4, -10, 10); fprintf(标准PSO最小值: %.6f\n, pso_best); % 实验组SAPSO rng(1); [sapso_best, sapso_record] iterateSAPSO(shubertfun, 30, 500, 2, 2, 0.9, 0.4, 500, 0.95, -10, 10); fprintf(SAPSO最小值: %.6f\n, sapso_best);rng(1)确保了两个算法在相同的随机数序列下初始化这样对比的是纯粹由算法机制差异导致的效果区别而不是随机运气。iteratePSO是标准PSO迭代函数可以自己按第2章的公式写一个结构与iterateSAPSO.m基本一致只是去掉SA接受判断和温度衰减两个环节。判断有效性的主要视角有两点。第一是收敛精度SAPSO最终稳定结果应明显优于标准PSO。第二是收敛路径画出每次迭代的record曲线SAPSO的曲线通常呈现“阶梯式下降”即在某一代温度还比较高时接受差解后跳出当前区域后续迭代中找到更好的局部极值。4.3 参数调优路径与常见陷阱调参顺序我一般按照T0、alpha、nPop三步走。先固定alpha 0.95、nPop 30跑一组实验观察前50代内接受较差解的比例。如果粒子从未离开过初始区域说明T0太低或降温太快将T0调大或alpha调高到0.98如果一直到迭代后期还能看到明显的位置跳跃说明温度还没有降到足够低此时应降低T0或降低alpha。完成T0和alpha的匹配后再适当调整nPop。现象可能原因调整方向前几十代就完全收敛目标函数值很低但未到-186alpha过小温度下降过快alpha增大到0.95~0.98大量粒子长期逃离较优区域收敛速度慢T0过高差解接受概率太大T0降低到100~300实验结果方差极大部分运行陷入-120附近种群规模过小探索不充分nPop增到40~50或增加maxIter曲线几乎不下降粒子一直随机飞行w衰减过快与T0过低双重抑制将wmin调高至0.5或放慢w衰减一个很常见的误用是把SA温度判断加在了gbest上而不是每个粒子上。一旦只对gbest做SA扰动扰动步长很难控制步长太大等于随机重启步长太小则完全无效。把这套代码中的SA判断放在每个粒子位置更新之后的好处是较差粒子有机会反复探索不同区域而gbest的更新策略没变整体收敛性不会受到破坏。在测试fun3这类Rastrigin函数时如果发现SAPSO比标准PSO的方差还大建议先检查边界约束用的是截断还是投影funx.m的修正逻辑对多峰函数影响显著。另外值得注意的一点是T0和fitness的尺度需要匹配。Shubert函数的fitness范围在-186到几百之间delta的量级通常是几十T0取500时exp(-delta/T0)约为0.9接受概率比较高如果换成fun2这类球形函数fitness动辄上千同样的T0会让接受概率趋近于1SA退化为随机搜索。切换测试函数时务必根据新函数的适应度量级重新校准T0这也是这类SAPSO代码在复用中最容易忽略的问题。5. 再进一步把SAPSO改造成通用求解器的四个细节5.1 加入退火重启机制标准SAPSO在温度降到很低之后粒子的探索能力基本消失如果此时gbest仍然不是全局最优算法已经没有自救能力。做法是在温度下降到初始温度的5%以下时如果gbest连续50代没有更新则触发重启保留当前gbest重新随机初始化一半粒子的位置和速度同时将温度回升到T0的30%给予新一轮跳出机会。这个策略不改变算法整体框架只需要在iterateSAPSO.m的温度衰减后加一个条件判断。重启会把收敛曲线拉出平台期但代价是增加迭代轮次实际使用时要设置一个最大重启次数比如3次防止无限循环。5.2 用历史gbest变化率做早停maxIter不总是越多越好。在工程场景中目标函数可能是耗时较长的仿真程序每一轮fitness评估都耗费真实时间。可以在iterateSAPSO.m中记录每次迭代f_min的变化量如果连续20代变化量小于1e-6就直接跳出主循环并返回当前gbest。需要注意的是早停阈值要结合目标函数数值尺度设置Shubert函数可以用1e-4而量级更大的函数需要按比例放大。这个改动能把计算成本降低20%到40%在批量跑参数实验时尤其划算。5.3 与函数评估预算结合与其固定maxIter不如把终止条件改成“总函数评估次数不超过N”。每次迭代包含nPop次fitness评估所以总评估次数等于nPop × maxIter。这样做的好处是在对比不同算法时保持计算预算一致而不是迭代次数一致因为不同算法的每代计算量可能有差异。修改方式是在主循环内部维护一个evalCounter变量每调用一次fun就加1超过上限直接终止并输出当前最优解。5.4 轨迹可视化的做法如果想观察SAPSO的搜索轨迹尤其是温度下降过程中粒子跳出局部极值的行为可以在iterateSAPSO.m的每个粒子更新位置后记录其坐标叠加在Shubert函数的等高线图上。绘制方法如下figure; [x1, x2] meshgrid(-10:0.1:10, -10:0.1:10); z arrayfun((a, b) shubertfun([a, b]), x1, x2); contour(x1, x2, z, 40); hold on; plot(history(1:100:end, 1), history(1:100:end, 2), r.-);history数组在每次迭代结束后记录gbest的坐标这样能直观看到gbest在不同局部极小点之间的跃迁路径。当温度下降后轨迹会逐渐稳定并最终停在-186.73附近——这是判断SA是否发挥作用的最后一道视觉验证。本文还有配套的精品资源点击获取