简介SPEA2算法MATLAB源程序是一套实现强度帕累托进化算法第二版的多目标优化代码主要面向需要同时优化多个冲突目标的科研人员、硕博研究生和工程开发者。压缩包内共十个文件其中九个为点M文件分别完成初始化、适应度计算、强度值评估、环境选择、交叉变异等核心流程另有一个文本说明文件整体仅六KB轻巧完整。目前已有九百四十五人学习关注源码附带了典型测试函数和绘图工具用户可运行后直接观察非支配解集在帕累托前沿上的分布便于理解算法搜索过程与收敛特性。用户可以按需调整目标函数、决策变量范围等参数快速验证不同多目标优化问题对于希望掌握SPEA2原理或进行二次开发的人群这套代码模块划分清晰、注释到位适合作为学习和实验的基础平台。1. SPEA2 算法 MATLAB 源程序真问题不是代码是权重怎么定找 SPEA2 算法 MATLAB 源程序的人多半已经在多目标优化里尝过手动定权重的苦把两三个目标压成一个加权和权重 0.5/0.5 跑出来的解换到 0.8/0.2 就完全是另一种分布调到最后全靠玄学。SPEA2强度帕累托进化算法第二代不需要你给权重它直接按 Pareto 支配关系保留一组互不支配的解再用近邻距离把拥挤区域里的冗余点挤掉最终输出一条分布均匀的 Pareto 前沿。这套算法的 MATLAB 实现核心只有三块适应度赋值、环境选择和进化算子并不神秘。如果你正在做多目标优化方向的课题或者想在模拟项目 X 里拿到一组均衡解而不是单个最优解这套源程序可以直接当底座使用关键是得知道参数怎么设、哪里容易翻车。2. SPEA2 的原理与 MATLAB 程序模块划分三个改进点为什么值得抄SPEA2 脱胎于 SPEA但两者的代码结构差异很大。这一章先讲清楚它改了什么、为什么这样改再落到 MATLAB 程序里模块怎么拆、数据怎么流。理解了这三件事第 3 章的代码你才看得下去。2.1 SPEA2 比 SPEA 强在哪适应度赋值、密度估计与截断策略早期 SPEA 有三个明显痛点一是适应度只由「它支配了多少个体」决定对大量互不支配的个体打分区分度不够二是完全没有密度信息种群容易扎进局部区域三是存档溢出时的截断策略过于粗暴经常把处于前端边缘的极端解删掉导致覆盖范围收缩。SPEA2 针对这三点做了对应修正这也是你拿到任何一份 SPEA2 源程序时最先应该找的三段代码对比点SPEA 的做法SPEA2 的做法适应度构成强度值 S(i)即 i 支配的个体数R(i) 支配 i 的所有个体的强度之和再叠加密度项密度估计无目标空间第 k 近邻距离的倒数k 通常取 sqrt(N)存档截断聚类或简单删末尾每轮删除最近邻距离最小的个体直到存档达到上限适应度的计算是 SPEA2 的灵魂。设种群与存档合并后共有 N 个个体每个个体 i 先算强度值 S(i) 被 i 支配的个体数然后算原始适应度 R(i) sum(S(j))其中 j 遍历所有支配 i 的个体。注意这里的措辞方向R 越小越好R 0 意味着没有任何个体支配它即它是非支配个体。接着用目标空间里的第 k 近邻距离 sigma_i 构造密度项 D(i) 1 / (sigma_i 2)sigma 越小说明周围越挤密度项越大适应度越差。最终 fitness R D。这个设计把「质量」和「多样性」压缩进了同一个标量环境选择一步到位不用像 NSGA-II 那样先分层再算拥挤距离。很多初学者把 SPEA2 当成黑匣子跑其实核心逻辑就是这么一句话越靠近前端、越不拥挤的个体适应度越小越该留下。2.2 MATLAB 源程序的模块拆分与数据流拿到一份 SPEA2 的 MATLAB 源程序先别急着跑先对照一下它的文件结构。我一般会把代码拆成四层参数脚本、主循环函数、算子函数、目标函数。参数脚本里放种群规模、存档规模、进化代数、交叉变异概率主循环函数负责合并种群与存档、调用适应度赋值、调用环境选择算子函数包括锦标赛选择、交叉、变异目标函数独立成一个 m 文件换测试问题只改这一处。整套程序的数据流是这样的初始化随机生成种群 pop外部存档 arch 置空每一代把 arch 和 pop 纵向拼接成 combined调用目标函数批量算出所有个体的目标值 fval用 fval 算每个个体的 SPEA2 适应度环境选择从 combined 中挑出 archSize 个适应度最好的个体组成下一代存档再从存档里做二元锦标赛选出父代经过交叉变异生成新的 pop进入下一代。这里和 NSGA-II 最大的区别是存档是个体跨代保留的唯一场所种群每一代都被后代替换而存档里的好个体能一直活到被更优个体取代为止。所以主循环里第一步一定是combined [arch; pop]这是 SPEA2 的框架特征如果哪份代码没做合并那它八成是简化版或者写错了。2.3 向量化还是逐个体循环MATLAB 实现的取舍MATLAB 写进化算法最大的坑是「把别的语言逐行翻译过来」。SPEA2 的计算热点有两个适应度赋值里的两两支配比较密度估计里的 N 乘 N 距离矩阵。这两处必须向量化否则慢到你怀疑人生而交叉变异反而适合写小循环逐行生成子代更直观、更容易调试。两两支配比较可以用bsxfun一行完成我把具体写法放到第 3 章。距离矩阵我习惯用三维广播把 fval 变成 N 乘 1 乘 m 和 1 乘 N 乘 m两者相减再平方求和开方一次得到整个距离矩阵。pdist2也能做这事但它依赖统计工具箱不是每个人都有用基本的三维广播最稳。有一个取舍要提醒存档规模 archSize 直接决定距离矩阵的大小。archSize 100 时100 乘 100 的距离矩阵几乎不占内存到了 500 就是一秒级别1000 以上每代耗时明显增长。所以别盲目把存档调大分布性不够时先检查 k 近邻参数和归一化而不是无脑加存档。3. SPEA2 算法 MATLAB 源程序从零写最小可运行版本这一章直接给代码。我没有用任何工具箱全部基于 MATLAB 基本函数。你把这几个函数存成独立的 m 文件再写一个参数脚本调用就能在标准测试问题上跑通。建议按我下面这个小节顺序抄每抄一段就知道它在算法里的位置。3.1 主循环与初始化先把框架搭起来主循环函数是整套程序的骨架它只做四件事合并、算适应度、环境选择、进化生成新种群。function arch spea2_run(prob, opt) % prob 包含问题定义objfun, nvar, lb, ub % opt 包含算法参数pop, arch, gen, pc, pm, etac, etam popSize opt.pop; archSize opt.arch; maxGen opt.gen; % 初始化种群每行一个个体均匀随机采样 pop repmat(prob.lb, popSize, 1) ... rand(popSize, prob.nvar) .* repmat(prob.ub - prob.lb, popSize, 1); arch []; % 外部存档初始为空这是 SPEA2 的正常起点 for t 1:maxGen combined [arch; pop]; % 存档与种群合并行号即个体编号 fval prob.objfun(combined); % 批量计算所有目标值 n size(combined, 1); fit spea2_fitness(fval); % 适应度赋值见 3.2 archIdx spea2_env_select(fval, fit, archSize); % 环境选择见 3.3 arch combined(archIdx, :); % 保留选中的个体进入下一代存档 archFval fval(archIdx, :); pop spea2_variation(arch, opt, prob.lb, prob.ub); % 交叉变异生成新种群 end arch spea2_nondominated(arch); % 最后再过滤一次只留非支配个体 end逻辑说明combined行号就是个体编号fval的行号和它一一对应这是整个程序不出错的关键约定。arch []不是 bugSPEA2 的存档是从进化里「生长」出来的第一代完全由 pop 的个体竞争产生后面才慢慢填满。spea2_nondominated是一个过滤函数把最终存档里被其他个体支配的个体删掉。如果你在跑测试算例建议保留这步否则最终画出来的前沿里可能混入几个明显偏下的点。目标函数按可批量计算的格式写以标准双目标测试算例 ZDT1 为例function y zdt1(x) % x 是 n 行 nvar 列矩阵返回 y 是 n 行 2 列 x1 x(:,1); g 1 9 * sum(x(:,2:end), 2) / (size(x,2) - 1); y [x1, g .* (1 - sqrt(x1 ./ g))]; end参数怎么给我常用的起点是 pop 100、arch 100、gen 250、pc 0.9、pm 1 / nvar、etac 20、etam 20。pm 设为 1/nvar 意味着每个个体平均有一个决策变量发生变异这是实数编码进化算法最常用的经验值。etac 和 etam 是交叉变异算子的分布指数越大生成的子代越接近父代先不要动它们。提示决策变量边界在初始化时用 repmat 展开了变异之后也要再 clip 回边界否则目标函数里一个越界值就可能让整代结果变成 NaN这个问题第 4 章会专门说。3.2 适应度赋值与密度估计SPEA2 的核心函数这一节是整份源程序里最不能抄错的部分。fitness 的方向搞反整个选择机制就变成「选差解」了。function fitness spea2_fitness(fval) % fval: n 行 m 列每行是一个个体的目标值最小化问题 n size(fval, 1); % 1. 强度值 S(i)i 支配了多少个个体 S zeros(n, 1); for i 1:n dominates all(bsxfun(le, fval(i,:), fval), 2) ... any(bsxfun(lt, fval(i,:), fval), 2); S(i) sum(dominates); end % 2. 原始适应度 R(i)所有支配 i 的个体的强度之和 R zeros(n, 1); for i 1:n dominatedBy all(bsxfun(ge, fval(i,:), fval), 2) ... any(bsxfun(gt, fval(i,:), fval), 2); R(i) sum(S(dominatedBy)); end % 3. 密度项 D(i)目标空间第 k 近邻距离的倒数 k max(1, round(sqrt(n))); diff reshape(fval, n, 1, []) - reshape(fval, 1, n, []); d sqrt(sum(diff.^2, 3)); d(1:n1:end) inf; % 自己到自己的距离置为 inf避免干扰 d sort(d, 2); sigma d(:, k); D 1 ./ (sigma 2); fitness R D; end逻辑说明第一个循环里fval(i,:)与每一行比较保证所有目标不差于对方保证至少一个目标严格优于对方两者同时成立就是 i 支配 j。R 的计算方向相反和判断的是「谁支配了 i」把支配者的 S 累加给 i。之所以这样设计是因为被一群「强个体」包围的个体即使它暂时是非支配的也不是值得保留的方向。密度项的 k 取round(sqrt(n))这是 SPEA2 原文的建议。sigma 越小说明该个体周围越拥挤D 越大适应度越差sigma 越大说明它越孤单D 趋近于 0.5适应度主要由 R 决定。这个互补关系让前端两端稀疏的极端解天然获得低适应度分布性就是这样保持的。注意如果两个目标量纲差距很大比如 f1 在 [0, 1]f2 在 [0, 10000]距离矩阵会被 f2 完全主导密度项形同虚设。我习惯在算 d 之前先把 fval 每一列按当前合并集合的 min-max 归一到 [0, 1] 再算距离。3.3 环境选择与存档截断保证分布性的关键环境选择决定哪些个体能进入下一代存档。SPEA2 的做法比 NSGA-II 的分层排序更直接先取所有非支配个体数量不够就用适应度从低到高补数量超出就反复删掉「最近邻距离最小」的个体直到达标。function idx spea2_env_select(fval, fitness, archSize) n size(fval, 1); % 第一步找出所有非支配个体 isND true(n, 1); for i 1:n for j 1:n if i ~ j all(fval(j,:) fval(i,:)) any(fval(j,:) fval(i,:)) isND(i) false; break; end end end idx find(isND); % 第二步非支配个体不足用适应度补满 if numel(idx) archSize rest find(~isND); [~, ord] sort(fitness(rest)); idx [idx; rest(ord(1 : archSize - numel(idx)))]; end % 第三步非支配个体超出上限逐个删除最近邻距离最小的个体 while numel(idx) archSize sub fval(idx, :); d dist_matrix(sub); d(1:size(d,1)1:end) inf; % 对角线置 inf [minDist, ~] min(d, [], 2); % 每个个体到最近邻的距离 [~, del] min(minDist); % 找到最拥挤的那个个体 idx(del) []; end end function d dist_matrix(x) n size(x, 1); diff reshape(x, n, 1, []) - reshape(x, 1, n, []); d sqrt(sum(diff.^2, 3)); end逻辑说明第一步的双层循环是纯 MATLAB 实现可读性好n 在 300 以内速度完全能接受如果 n 上千这段可以用向量化重写但我建议先跑通再优化。第二步的sort(fitness(rest))必须升序因为 fitness 越小越优。第三步的删除策略是 SPEA2 的典型实现每次找一个「周围最挤」的个体删掉而不是一次性算完所有距离再批量删后者容易把分布性毁掉。变异算子这里用多项式变异交叉算子可以用模拟二进制交叉 SBX它们加一起生成新一代种群function newPop spea2_variation(arch, opt, lb, ub) n size(arch, 1); newPop zeros(opt.pop, size(arch, 2)); for i 1:2:opt.pop p1 tournament(arch); % 二元锦标赛选中适应度好的父代 p2 tournament(arch); if rand opt.pc [c1, c2] sbx_crossover(p1, p2, opt.etac); else c1 p1; c2 p2; end newPop(i, :) poly_mutation(c1, lb, ub, opt.pm, opt.etam); newPop(i1, :) poly_mutation(c2, lb, ub, opt.pm, opt.etam); end end锦标赛选择实现很简单从存档里随机抽两个个体比 fitness小的那个胜出。SBX 和多项式变异的公式是教科书标准式复制时注意两点beta 的指数要除以etac 1而不是 etac变异后必须做mutant min(max(mutant, lb), ub)边界钳制。我早期就是忘了这行 clip导致跑 ZDT 系列时经常在 30 代左右整片 NaN。4. SPEA2 调试避坑常见问题与排查代码抄对了不代表能一次跑通。这一章写的是我在调试 SPEA2 类实现时真正踩过、也帮别人排查过的五类问题每条按现象、原因、解决三步写你遇到类似症状可以直接对号入座。4.1 存档越界报错合并矩阵维数对不上现象程序前十几代正常跑到某代突然报「索引超出数组边界」或「数组索引必须为正整数」断点停在 spea2_fitness 或 env_select 里。原因最常见的是在函数里把popSize写死当成合并后的个体数。实际上combined [arch; pop]之后总个体数是n size(combined, 1)它会随存档是否满员变化。任何用固定行数创建的零矩阵在存档数量变化后都会索引错位。解决所有和个体数相关的矩阵一律用n size(combined, 1)现场取不要用传入的opt.pop或opt.arch硬算。我自己的习惯是每个函数入口第一行就写n size(fval, 1)保证后面的 R、S、isND 全部跟着实际行数走。排查时先用dbstop if error看报错那行是哪个矩阵越界再往回找它是不是用错了维度。4.2 种群糊成一个点密度项没起作用现象跑完画 Pareto 前沿一百多个存档个体挤在一小片区域前端两端完全没有点甚至退化到近似单个点。原因这题有三个常见根源。第一密度项没做目标归一化f2 量纲大时距离矩阵被它主导f1 方向的分布完全失效第二k取太大当k接近 n 时第 k 近邻距离变成了全局最大间距所有个体的 sigma 都差不多D 全等于 0.5 附近密度项退化成常数第三密度项的权重被 R 压住在 R 都等于 0 的非支配区里 D 的绝对值太小区分度不够。解决先把 fval 的每一列归一化到 [0,1] 再算密度把k固定为max(1, round(sqrt(n)))别手动改成一个别的值如果非支配个体太多导致 R 全是 0可以在 fitness 里给 D 乘一个放大系数比如把fitness R D改成fitness R 10 * D让分布压力起来。改完看前端的两个端点有没有逐渐外扩这是密度项恢复工作的标志。4.3 适应度里出现 NaN 或 Inf越界传播现象某代之后 fitness 整列变成 NaN环境选择选不出个体程序卡死或报错。有时前端图上出现一个点飞到极端位置像是坐标轴都被拉崩了。原因决策变量越界进入不可行域目标函数里开方、对数、除法的输入变成负数或零产生 NaN还有一种情况是多项式变异里(2*(1-r))在 r 恰好等于 1 时变成 0再取负指数直接给 Inf。最麻烦的是 NaN 有传染性——支配判断里只要包含 NaNall和any的结果全是 false整支配关系失效后续全是糊涂账。解决三重防线。第一变异函数末尾强制钳制边界第二算完目标函数后立刻断言一次all(isfinite(fval(:)))一旦发现非有限值直接报错并输出当前代数方便定位是哪类算子在爆雷第三在支配判断前用isnan(fval)把 NaN 视作被支配宁可丢个体也不让 NaN 继续扩散。前两项就够日常用了第三项是我在目标函数本身有坑时加的兜底。4.4 两次运行结果差异大随机性没有被固定现象同一套代码、同一个测试函数昨天跑的前沿光滑均匀今天再跑缺了一角或者边界点直接消失。换一台机器结果又不一样。原因MATLAB 的全局随机流是共享的你只要在跑 SPEA2 之前执行过任何带随机数的代码种子就变了。更隐蔽的是锦标赛选择里连续取两个父代如果随机数不独立局部随机波动会被存档截断放大小种群时尤其明显。解决主脚本开头写rng(2027)固定种子任何实验都从这个种子出发保证可复现正式报告结果时不要只跑一次跑 20 次独立实验记录 HV 或 IGD 指标的均值和方差。我自己的习惯是写一个run_experiments.m循环改种子把每代前沿的指标增量存成矩阵最后画箱线图。所谓「这算法不稳定」的结论十有八九是种子没固定造成的统计假象。4.5 跑得很慢循环里复制矩阵现象前 50 代很快到 100 代后每代耗时明显变长总数 250 代的程序跑了 40 分钟。profile 一查热点在 env_select 和 density。原因两个典型的慢点。一是距离矩阵用三层 for 循环逐对计算复杂度 O(n^3)二是存档截断的 while 循环里每次删除一个元素MATLAB 都要复制一遍剩余数组当 archSize 超过 300 时开销陡增。还有一种隐蔽情况目标函数接收的是单行向量输入被迫在 for 循环里逐个算目标值。解决距离矩阵改用三维广播一次成型前面 3.2 的代码就是这么写的截断循环保留但不要在里面重算整个距离矩阵可以先把所有个体到最近邻的距离向量算出来再排序后从远到近决定删谁目标函数写成接收整个矩阵。最后如果个体总数长期超过 800把距离矩阵的 dtype 明确为 single能再省一半内存。5. 验证 SPEA2 实现正确性的三个技巧从乱点云到光滑前端代码跑通不等于算法对。我验证一套新的 SPEA2 实现是否正常从来不看最终那一张图而是用三个递进的做法。第一个技巧用有已知真实前沿的双目标测试算例做回归。比如 ZDT1它的真实 Pareto 前端是 f2 1 - sqrt(f1)画出来是一条单调凸曲线。每次迭代把当前存档画成点云叠在这条曲线上眼睛一眼就能判断收敛方向对不对。具体做法是在主循环里加三行plot(fval(archIdx,1), fval(archIdx,2), o); hold on; ezplot((x) 1 - sqrt(x), [0.01, 1]); % 真实前端参考线 drawnow;如果点云从右上角往左下角逐步贴向这条参考线说明环境选择的方向正确如果点云停在中间某个位置就散掉八成是密度项或截断逻辑出了问题。第二个技巧用超体积指标 HV 做数值验收。双目标最小化问题里取参考点略大于每个目标的最大值把非支配解按 f1 升序排列HV 就是每个解与参考点围成的矩形面积逐段累加。正确的实现 HV 应该是单调上升的最终稳定在一个平台如果 HV 曲线出现下降说明存档里混入了被支配解或者截断时删掉了不该删的边界点。这个指标还能帮你调 archSizeHV 平台平缓后继续加大存档收益会明显递减这时候就该停。第三个技巧看间距而不是只看密度。把最终存档按 f1 排序计算相邻两个解在目标空间里的欧氏距离求标准差。标准差越小说明分布越均匀。我自己用下来ZDT1 这类光滑前端在 100 个存档个体时间距标准差一般能压到平均间距的 20% 以下如果这个比例超过 50%说明密度项被 R 压得没有话语权按第 4.2 节的方法放大 D 的权重即可。最后说一个我自己的教训最早写 SPEA2 源程序时我把所有函数塞进一个 300 行的 m 文件里每次调参都要从头跑一遍定位问题靠到处加 disp改一次参数等五分钟。后来拆成六个小函数每个函数只做一件事参数脚本单独放调试时间从半小时缩到两分钟。那之后我再拿到任何进化算法的参考代码第一件事永远是按模块拆开再跑。这套 SPEA2 的 MATLAB 实现也一样先跑通 ZDT1再换 ZDT2、ZDT4逐步替换你自己的业务目标函数比上来就改算法结构稳妥得多。希望帮到你。本文还有配套的精品资源点击获取
