企业官网建设流程全解析

简介：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 源程序时最先应该找的三段代码：

适应度的计算是 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 = []不是 bug，SPEA2 的存档是从进化里「生长」出来的，第一代完全由 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:n+1: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(i+1, :) = 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 有传染性——支配判断里只要包含 NaN，all和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 曲线出现下降，说明存档里混入了被支配解，或者截断时删掉了不该删的边界点。这个指标还能帮你调 archSize：HV 平台平缓后继续加大存档，收益会明显递减，这时候就该停。

第三个技巧：看间距而不是只看密度。把最终存档按 f1 排序，计算相邻两个解在目标空间里的欧氏距离，求标准差。标准差越小说明分布越均匀。我自己用下来，ZDT1 这类光滑前端在 100 个存档个体时，间距标准差一般能压到平均间距的 20% 以下；如果这个比例超过 50%，说明密度项被 R 压得没有话语权，按第 4.2 节的方法放大 D 的权重即可。

最后说一个我自己的教训：最早写 SPEA2 源程序时，我把所有函数塞进一个 300 行的 m 文件里，每次调参都要从头跑一遍，定位问题靠到处加 disp，改一次参数等五分钟。后来拆成六个小函数，每个函数只做一件事，参数脚本单独放，调试时间从半小时缩到两分钟。那之后我再拿到任何进化算法的参考代码，第一件事永远是按模块拆开再跑。这套 SPEA2 的 MATLAB 实现也一样，先跑通 ZDT1，再换 ZDT2、ZDT4，逐步替换你自己的业务目标函数，比上来就改算法结构稳妥得多。希望帮到你。

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