简介面向多目标优化问题这套MATLAB代码包以差分进化DE为内核实现了无偏好与带偏好两类共六种算法方便研究者对比不同决策模式下的求解差异。无偏好部分包括基于非支配排序的DEMO和基于指标评价的IBEA偏好部分则涵盖R-DEMO、PBEA以及提出的PAR-DEMO(nds)与PAR-DEMO(ε)分别通过参考点或ε指标引导搜索。压缩包共24个文件以22个.m程序文件为主体辅以README.md和PDF说明文档总大小只有307KB便于快速下载和阅读。代码支持在Octave及Matlab环境运行示例中整合了DTLZ测试问题读者可以直接修改目标函数、参数或偏好设置评估不同算法在收敛性和多样性上的表现也可将模块嵌入自己的实验框架。目前已有384人参与学习适合正在学习多目标进化计算、需要可运行代码进行对比实验的研究生或工程师。1. 多目标差分进化从单目标到帕累托前沿的工程落地实际工程优化里冲突目标才是常态。比如传感器布局既要覆盖面积大又希望成本低复合材料铺层既要刚度大又要求重量轻。这类问题不存在唯一最优解只有一个帕累托前沿。差分进化DE的全局搜索能力在单目标问题上已经够用但它的选择机制是“新解更好就替换”一旦目标变成两个就必须引入支配关系和非支配排序这就是多目标差分进化MODE的核心改造点。对于手里有MATLAB的工程师与其去翻论文复现公式不如直接拿一段能跑的MODE代码在ZDT标准测试集上看到帕累托前沿再按自己的目标函数替换适配。下面我按照自己调这类算法的方式从基础实现讲到变体选型最后给出一份可直接运行的完整脚本。2. 基础多目标差分进化MODE的MATLAB实现与参数解析多目标差分进化并不是一套全新的算法它保留了DE的变异和交叉算子只把选择部分换成“非支配排序 拥挤距离”。所以先要把这两个概念在代码里落地否则后面的变体都无从谈。2.1 支配关系与拥挤距离环境选择怎么做在多目标优化中解A支配解B表示A在所有目标上都不比B差并且至少在一个目标上严格优于B。所有不被任何其他解支配的个体构成第一非支配层去掉第一层后再找第二层以此类推。环境选择时我优先保留层数靠前的个体。但同层内部如何取舍这需要用拥挤距离衡量个体的稀疏程度。以二维目标为例对当前层个体按每个目标方向排序取相邻个体在各目标坐标归一化后的差值和两端个体距离设成无穷大。保种时层数相同就看拥挤距离距离大的留这样种群就不会挤在解空间某一段。2.1.1 一句话理解DE的算子不变只换选择变异仍然用V X_r1 F * (X_r2 - X_r3)交叉仍然用二项式交叉但生成子代后不再“一换一”地替换而是把父代和子代合并再通过上述环境选择从2NP个个体里找出NP个。这是绝大多数MODE变体的共同骨架。2.2 一个可直接运行的基础MODE代码骨架下面这段函数是我常用的模板保存为MODE_basic.m。变异用DE/rand/1选择使用非支配排序加拥挤距离。你只需要把自己的目标函数做成句柄fun让fun(x)返回一个行向量即可。function [pf, pfObj] MODE_basic(fun, D, lb, ub, numObj, NP, maxGen, F, CR) % 多目标差分进化基础版 % 选择机制非支配排序 拥挤距离 lb lb(:); ub ub(:); if length(lb) 1 lb repmat(lb, 1, D); end if length(ub) 1 ub repmat(ub, 1, D); end % 初始化种群 X repmat(lb, NP, 1) rand(NP, D) .* repmat(ub-lb, NP, 1); Obj zeros(NP, numObj); for i 1:NP Obj(i, :) fun(X(i, :)); end for gen 1:maxGen U zeros(NP, D); for i 1:NP % 选择三个互不相同的随机个体且不能是当前i r randperm(NP, 3); while ismember(i, r) r randperm(NP, 3); end % DE/rand/1 变异 donor X(r(1), :) F * (X(r(2), :) - X(r(3), :)); % 二项式交叉 jrand randi(D); mask rand(1, D) CR; mask(jrand) true; u donor .* mask X(i, :) .* (1-mask); % 边界重置为边界值 u(u lb) lb(u lb); u(u ub) ub(u ub); U(i, :) u; end % 评估子代 subObj zeros(NP, numObj); for i 1:NP subObj(i, :) fun(U(i, :)); end % 合并父代与子代环境选择出NP个个体 [X, Obj] nonDomSelection([X; U], [Obj; subObj], NP); end % 提取最终第一前沿 [pf, pfObj] extractParetoFront(X, Obj); end %% 非支配排序 拥挤距离选择 function [Xsel, ObjSel] nonDomSelection(X, Obj, NP) N size(Obj, 1); % 计算支配关系 nP zeros(N, 1); S cell(N, 1); fronts {}; F 1; for i 1:N S{i} []; for j 1:N if dominates(Obj(i, :), Obj(j, :)) S{i} [S{i}, j]; elseif dominates(Obj(j, :), Obj(i, :)) nP(i) nP(i) 1; end end if nP(i) 0 fronts{1} [fronts{1}, i]; end end % 逐层剥离获得完整非支配层 while ~isempty(fronts{F}) cur []; for i 1:length(fronts{F}) idx fronts{F}(i); for k 1:length(S{idx}) j S{idx}(k); nP(j) nP(j) - 1; if nP(j) 0 fronts{F1} [fronts{F1}, j]; cur [cur, j]; end end end F F 1; fronts{F} cur; end % 逐层选取最后一层用拥挤距离截断 selected []; for f 1:F if isempty(fronts{f}) continue; end idx fronts{f}; if length(selected) length(idx) NP selected [selected, idx]; else need NP - length(selected); dist crowdingDistance(Obj(idx, :)); [~, order] sort(dist, descend); selected [selected, idx(order(1:need))]; break; end end Xsel X(selected, :); ObjSel Obj(selected, :); end function d dominates(a, b) d all(a b) any(a b); end function dist crowdingDistance(Obj) n size(Obj, 1); numObj size(Obj, 2); dist zeros(n, 1); for m 1:numObj [~, order] sort(Obj(:, m)); dist(order(1)) inf; dist(order(end)) inf; if n 2 fmin Obj(order(1), m); fmax Obj(order(end), m); if fmax fmin continue; end for k 2:n-1 dist(order(k)) dist(order(k)) ... (Obj(order(k1), m) - Obj(order(k-1), m)) / (fmax - fmin); end end end end %% 提取第一前沿 function [pf, pfObj] extractParetoFront(X, Obj) N size(Obj, 1); isPareto true(N, 1); for i 1:N for j 1:N if i ~ j dominates(Obj(j, :), Obj(i, :)) isPareto(i) false; break; end end end pf X(isPareto, :); pfObj Obj(isPareto, :); end代码里nonDomSelection的复杂度是O(N^2)NP不要设太大。extractParetoFront在最终种群上做一次全对全比较足以提取第一前沿。所有目标值默认按最小化处理如果你的问题是最大化先在fun里取负。2.3 必调参数NP、F、CR与收敛性的权衡这套模板真正的开关只有几个NP、F、CR、maxGen。它们对结果的影响我总结成下表方便你直接按表试参。参数推荐范围对优化过程的影响我的习惯NP50200太小种群易早熟太大单代计算量显著增长取决策变量维数的10倍至少100F0.30.9控制变异步长大F利于探索小F利于局部精调先用0.5收敛慢再降CR0.60.95交叉概率高子代保留父代信息少多样性变强ZDT类问题固定0.9maxGen100500影响最终收敛程度和运行时间先跑200代看前沿形状F和CR有一个常见的误区不是F越大越容易跳出局部最优。在多目标版本里F过大反而容易让变异后的个体频繁突破边界虽然代码会把边界重置但大量个体堆在边界上会损失内部多样性。如果你发现前沿缺失中间段先检查F是否超过0.7。边界处理这里用的是“越界重置为边界值”实现简单。但我在处理有界工程问题时如果发现种群有聚集在边界上的趋势会改成“越界随机重生”即越界维度重新在lb(j)到ub(j)内随机赋值代价是多写一行效果却更稳。3. 多目标差分进化变体对比与基于分解的MODE实现思路基础MODE能解决多数双目标问题但在目标数增加或问题形态复杂时单靠一种选择机制不够。变体的差异主要体现在选择机制、分解策略和种群维护上。3.1 主流变体的算法结构与适用场景这里列几个MATLAB社区里常见的多目标差分进化变体。变体名称核心机制典型适用场景GDE3每次迭代先扩展种群再用非支配排序和剪枝删除多余个体需要同时处理约束和目标连续问题DEMO每个父代强制生成一个子代合并后用NSGA-II风格选择双目标通用实现简单MOEA/D-DE将多目标分解为多个单目标子问题用邻域进行DE变异目标数23已知偏好或权重NSDE把NSGA-II的非支配排序直接替换DE的选择步骤决策变量多前沿不连续要注意GDE3和DEMO在代码上差别不大DEMO更接近基础MODE而GDE3多了一个“剪枝”阶段。选型上我没有固定答案一般如果问题维度在30以下我会直接用第二节的基础MODE如果需要更均匀的前沿分布就换MOEA/D-DE。3.2 在MATLAB中实现基于分解的MOEA/D-DE变体MOEA/D-DE的思路是用一组权重向量把多目标拆成若干子问题每个子问题对应一个聚合函数值。以Tchebycheff聚合函数为例子问题k的聚合值由权重向量lambda_k和目标值obj计算。更新时只在邻域内选择父代做变异和交叉新解如果聚合值更小就替换邻域内的旧解。下面是核心更新片段完整代码在你把权重向量生成加进去后就能跑% lambda: NP×numObj 的权重向量矩阵每行一个子问题 % B: NP×T 邻域索引矩阵假设每个子问题有T个邻居 % X, Obj: 当前种群和对应目标值 for i 1:NP % 从邻域B(i,:)中随机选择两个不同个体作为r1, r2 nb B(i, :); idx randperm(length(nb), 2); r1 nb(idx(1)); r2 nb(idx(2)); % DE变异和交叉与基础MODE相同 donor X(i, :) F * (X(r1, :) - X(r2, :)); mask rand(1, D) CR; mask(randi(D)) true; u donor .* mask X(i, :) .* (1 - mask); u min(max(u, lb), ub); % 计算新解目标值 newObj fun(u); % 更新邻域内聚合值更差解 for k 1:length(nb) j nb(k); % Tchebycheff聚合值z为理想点即当前各目标最小值 g_old max(lambda(j, :) .* abs(Obj(j, :) - z)); g_new max(lambda(j, :) .* abs(newObj - z)); if g_new g_old X(j, :) u; Obj(j, :) newObj; end end % 更新理想点z z min(z, newObj); end这个片段里最关键的部分是“用聚合值g_new和g_old比较”它直接替代了非支配排序。lambda的分布会决定前沿均匀性我一般用MATLAB自带的lhsdesign生成均匀权重再用knnsearch找每个子问题的邻近索引。注意z是动态更新的理想点如果目标量纲差异大先做归一化再算聚合值。3.3 变体选型根据问题特征选择算法选型没有银弹但可以参考这几个经验。目标只有两个且计算目标函数一次很贵用基础MODEmaxGen控制在200以内不要换分解法。目标前沿不连续或呈碎片状优先试NSDE这类基于非支配的变体因为分解法在碎片前沿上容易丢失部分子问题。参考向量或偏好已知比如工厂里指定“成本不能超过某值”用MOEA/D-DE权重向量可以直接带偏好。决策变量超过100建议先把代码的变异循环向量化否则无论哪种变体都会卡在非支配排序上。我在实际项目里往往先用基础MODE跑通流程验证目标函数无误再用MOEA/D-DE或GDE3去细调。这种顺序能省不少排错时间。4. 在MATLAB中跑通ZDT测试集完整脚本、性能指标与排错有了一套基础MODE下一步是在标准测试函数上看到实际效果。ZDT1是双目标连续凸前沿适合做第一个验证。4.1 定义ZDT1目标函数与算法入口ZDT1的决策变量维度可以任意设但常用30为了快速验证我这里设成10。它有两个目标f1 x(1)g 1 9 * sum(x(2:D))/(D-1)f2 g * (1 - sqrt(f1/g))。决策变量都在[0,1]之间。function f zdt1(x) D length(x); f1 x(1); g 1 9 * sum(x(2:end)) / (D - 1); f2 g * (1 - sqrt(f1/g)); f [f1, f2]; end把这段保存为zdt1.m放在和MODE_basic.m同一目录下。注意fun返回行向量这是第二节代码里的约定。4.2 运行完整脚本从MODE_basic到帕累托前沿绘图下面这个主脚本可以直接执行clc; clear; D 10; lb zeros(1, D); ub ones(1, D); numObj 2; NP 100; maxGen 200; F 0.5; CR 0.9; tic; [pf, pfObj] MODE_basic(zdt1, D, lb, ub, numObj, NP, maxGen, F, CR); toc; % 绘制算法前沿与真实前沿 t linspace(0, 1, 500); trueF [t, 1 - sqrt(t)]; figure; plot(pfObj(:,1), pfObj(:,2), bo, MarkerSize, 4); hold on; plot(trueF(:,1), trueF(:,2), r-, LineWidth, 1.2); xlabel(f1); ylabel(f2); legend(MODE result, True PF); % 计算IGD与HV igd computeIGD(pfObj, trueF); ref [1.1, 1.1]; hv computeHV(pfObj, ref); fprintf(IGD %.4f, HV %.4f\n, igd, hv);computeIGD计算反向世代距离computeHV计算超体积。两者的实现如下可以直接放在主脚本末尾或保存为子函数。function igd computeIGD(pf, truePF) n size(truePF, 1); dist zeros(n, 1); for i 1:n d sqrt(sum((pf - truePF(i, :)).^2, 2)); dist(i) min(d); end igd mean(dist); end function hv computeHV(pf, ref) pf sortrows(pf, 1); hv 0; lastX 0; for i 1:size(pf, 1) x pf(i, 1); if x lastX continue; end width x - lastX; hv hv width * (ref(2) - pf(i, 2)); lastX x; end hv hv (ref(1) - lastX) * ref(2); endcomputeHV只适用于两个目标且前沿按f1单调递增的情况ZDT1恰好满足。如果你换到其他测试函数建议用paretoset或MATLAB File Exchange上的通用HV函数。4.3 计算IGD与HV指标评估收敛性和多样性IGD越小越好它衡量算法前沿与真实前沿的平均距离。HV越大越好参考点通常取略大于各目标最大值。运行上面脚本在F0.5, CR0.9, NP100, maxGen200时我得到IGD约0.01、HV约0.65。不同随机种子会有波动属正常。如果IGD偏大说明前沿离真实前沿远先加大maxGen。如果HV不高说明前沿覆盖不够广泛优先尝试调大NP或降低F。4.4 四个常见运行错误及其解决多目标优化代码最常见的错误几乎都集中在这四处。现象原因解决方向结果全是NaN目标函数中存在除零或log负数检查zdt1里的g当D1时除零前沿只有一个点种群早熟或拥挤距离失效调大NP到150降低F到0.3运行时间暴涨非支配排序O(N^2)NP超过300将NP控制在200以内前沿分布不均边界处理把所有个体压到边界把越界重置改成随机重生第一条里的坑最隐蔽ZDT1当D1时sum(x(2:end))为空g19*0/0就会NaN。我建议在zdt1.m开头加一个断言assert(length(x) 2, ZDT1需要至少2个决策变量);另外MATLAB版本差异也会影响randperm(NP,3)的行为R2023b里它默认返回行向量R2018b也一致。如果你的版本很老可以改成randperm(NP)再取前三个。5. 进阶技巧向量化加速与自适应参数让多目标差分进化更快当问题收敛到前沿形状稳定后瓶颈通常在运行速度。这里给两个我能直接落地的技巧。5.1 向量化变异与交叉去掉内层for循环第二节的代码在NP100时内层循环尚可但NP200时每一代要算200次变异和200次目标函数。目标函数是黑盒时无法加速但变异和交叉可以整段向量化。关键在于一次生成所有变异向量而不逐个for% r1, r2, r3 都是NP×1的随机索引 r1 randperm(NP); r2 randperm(NP); r3 randperm(NP); % 保证互不相同 for i 1:NP while r1(i)r2(i) || r1(i)r3(i) || r2(i)r3(i) r1(i) randi(NP); r2(i) randi(NP); r3(i) randi(NP); end end donors X(r1,:) F * (X(r2,:) - X(r3,:)); mask rand(NP, D) CR; jrand randi(D, NP, 1); mask(sub2ind([NP, D], (1:NP), jrand)) true; U donors .* mask X .* (1 - mask);注意randperm(NP)生成的索引天然无重复但如果NP很小直接对每一维重新随机更安全。向量化后目标函数评估仍然是逐行进行的如果fun支持批量输入矩阵每一行是一个个体还可以再进一步向量化。5.2 自适应F与CR用历史成功率动态调整固定参数能跑通但前沿形状差时往往需要调两次。最简单的一种自适应策略是“中位数更新”每20代统计一下当前代较父代有所改进的个体占比若改进比例高说明搜索方向有效就调大F和CR一点反之调小。if mod(gen, 20) 0 improvement sum(newObjBetter) / NP; if improvement 0.3 F min(F * 1.1, 0.9); CR min(CR * 1.05, 1.0); elseif improvement 0.1 F max(F * 0.9, 0.2); CR max(CR * 0.95, 0.5); end end这个技巧不严谨但在工程上足够它不会让参数突变只做微调。你可以在第三节的MOEA/D-DE变体里也用同样的方法不需要改动选择部分。5.3 如何扩展自己的变体把约束处理加入MODE工程问题几乎都带约束。最常见的做法是把约束违反量作为第三个目标或者用“约束支配”修改2.1节里的dominates函数。我推荐后者改动最小function d dominates(a, b) cv (x) max(x(1), 0); % 简化示意实际约束在外部 % 下面需要比较两个解的约束违反量和目标值 d all(a(1:numObj) b(1:numObj)) any(a(1:numObj) b(1:numObj)); end更好的做法是保存每个解的约束违反总量在dominates里先比较约束违反总量违反小的支配违反大的。这样基础MODE就变成了约束多目标差分进化而变体逻辑完全不变。这个思路在GDE3里被大量使用你只需在fun返回目标值时同时返回约束值再在环境选择里加一个判断即可。本文还有配套的精品资源点击获取 多目标差分进化算法MATLAB实战：帕累托前沿求解与变体选型
简介面向多目标优化问题这套MATLAB代码包以差分进化DE为内核实现了无偏好与带偏好两类共六种算法方便研究者对比不同决策模式下的求解差异。无偏好部分包括基于非支配排序的DEMO和基于指标评价的IBEA偏好部分则涵盖R-DEMO、PBEA以及提出的PAR-DEMO(nds)与PAR-DEMO(ε)分别通过参考点或ε指标引导搜索。压缩包共24个文件以22个.m程序文件为主体辅以README.md和PDF说明文档总大小只有307KB便于快速下载和阅读。代码支持在Octave及Matlab环境运行示例中整合了DTLZ测试问题读者可以直接修改目标函数、参数或偏好设置评估不同算法在收敛性和多样性上的表现也可将模块嵌入自己的实验框架。目前已有384人参与学习适合正在学习多目标进化计算、需要可运行代码进行对比实验的研究生或工程师。1. 多目标差分进化从单目标到帕累托前沿的工程落地实际工程优化里冲突目标才是常态。比如传感器布局既要覆盖面积大又希望成本低复合材料铺层既要刚度大又要求重量轻。这类问题不存在唯一最优解只有一个帕累托前沿。差分进化DE的全局搜索能力在单目标问题上已经够用但它的选择机制是“新解更好就替换”一旦目标变成两个就必须引入支配关系和非支配排序这就是多目标差分进化MODE的核心改造点。对于手里有MATLAB的工程师与其去翻论文复现公式不如直接拿一段能跑的MODE代码在ZDT标准测试集上看到帕累托前沿再按自己的目标函数替换适配。下面我按照自己调这类算法的方式从基础实现讲到变体选型最后给出一份可直接运行的完整脚本。2. 基础多目标差分进化MODE的MATLAB实现与参数解析多目标差分进化并不是一套全新的算法它保留了DE的变异和交叉算子只把选择部分换成“非支配排序 拥挤距离”。所以先要把这两个概念在代码里落地否则后面的变体都无从谈。2.1 支配关系与拥挤距离环境选择怎么做在多目标优化中解A支配解B表示A在所有目标上都不比B差并且至少在一个目标上严格优于B。所有不被任何其他解支配的个体构成第一非支配层去掉第一层后再找第二层以此类推。环境选择时我优先保留层数靠前的个体。但同层内部如何取舍这需要用拥挤距离衡量个体的稀疏程度。以二维目标为例对当前层个体按每个目标方向排序取相邻个体在各目标坐标归一化后的差值和两端个体距离设成无穷大。保种时层数相同就看拥挤距离距离大的留这样种群就不会挤在解空间某一段。2.1.1 一句话理解DE的算子不变只换选择变异仍然用V X_r1 F * (X_r2 - X_r3)交叉仍然用二项式交叉但生成子代后不再“一换一”地替换而是把父代和子代合并再通过上述环境选择从2NP个个体里找出NP个。这是绝大多数MODE变体的共同骨架。2.2 一个可直接运行的基础MODE代码骨架下面这段函数是我常用的模板保存为MODE_basic.m。变异用DE/rand/1选择使用非支配排序加拥挤距离。你只需要把自己的目标函数做成句柄fun让fun(x)返回一个行向量即可。function [pf, pfObj] MODE_basic(fun, D, lb, ub, numObj, NP, maxGen, F, CR) % 多目标差分进化基础版 % 选择机制非支配排序 拥挤距离 lb lb(:); ub ub(:); if length(lb) 1 lb repmat(lb, 1, D); end if length(ub) 1 ub repmat(ub, 1, D); end % 初始化种群 X repmat(lb, NP, 1) rand(NP, D) .* repmat(ub-lb, NP, 1); Obj zeros(NP, numObj); for i 1:NP Obj(i, :) fun(X(i, :)); end for gen 1:maxGen U zeros(NP, D); for i 1:NP % 选择三个互不相同的随机个体且不能是当前i r randperm(NP, 3); while ismember(i, r) r randperm(NP, 3); end % DE/rand/1 变异 donor X(r(1), :) F * (X(r(2), :) - X(r(3), :)); % 二项式交叉 jrand randi(D); mask rand(1, D) CR; mask(jrand) true; u donor .* mask X(i, :) .* (1-mask); % 边界重置为边界值 u(u lb) lb(u lb); u(u ub) ub(u ub); U(i, :) u; end % 评估子代 subObj zeros(NP, numObj); for i 1:NP subObj(i, :) fun(U(i, :)); end % 合并父代与子代环境选择出NP个个体 [X, Obj] nonDomSelection([X; U], [Obj; subObj], NP); end % 提取最终第一前沿 [pf, pfObj] extractParetoFront(X, Obj); end %% 非支配排序 拥挤距离选择 function [Xsel, ObjSel] nonDomSelection(X, Obj, NP) N size(Obj, 1); % 计算支配关系 nP zeros(N, 1); S cell(N, 1); fronts {}; F 1; for i 1:N S{i} []; for j 1:N if dominates(Obj(i, :), Obj(j, :)) S{i} [S{i}, j]; elseif dominates(Obj(j, :), Obj(i, :)) nP(i) nP(i) 1; end end if nP(i) 0 fronts{1} [fronts{1}, i]; end end % 逐层剥离获得完整非支配层 while ~isempty(fronts{F}) cur []; for i 1:length(fronts{F}) idx fronts{F}(i); for k 1:length(S{idx}) j S{idx}(k); nP(j) nP(j) - 1; if nP(j) 0 fronts{F1} [fronts{F1}, j]; cur [cur, j]; end end end F F 1; fronts{F} cur; end % 逐层选取最后一层用拥挤距离截断 selected []; for f 1:F if isempty(fronts{f}) continue; end idx fronts{f}; if length(selected) length(idx) NP selected [selected, idx]; else need NP - length(selected); dist crowdingDistance(Obj(idx, :)); [~, order] sort(dist, descend); selected [selected, idx(order(1:need))]; break; end end Xsel X(selected, :); ObjSel Obj(selected, :); end function d dominates(a, b) d all(a b) any(a b); end function dist crowdingDistance(Obj) n size(Obj, 1); numObj size(Obj, 2); dist zeros(n, 1); for m 1:numObj [~, order] sort(Obj(:, m)); dist(order(1)) inf; dist(order(end)) inf; if n 2 fmin Obj(order(1), m); fmax Obj(order(end), m); if fmax fmin continue; end for k 2:n-1 dist(order(k)) dist(order(k)) ... (Obj(order(k1), m) - Obj(order(k-1), m)) / (fmax - fmin); end end end end %% 提取第一前沿 function [pf, pfObj] extractParetoFront(X, Obj) N size(Obj, 1); isPareto true(N, 1); for i 1:N for j 1:N if i ~ j dominates(Obj(j, :), Obj(i, :)) isPareto(i) false; break; end end end pf X(isPareto, :); pfObj Obj(isPareto, :); end代码里nonDomSelection的复杂度是O(N^2)NP不要设太大。extractParetoFront在最终种群上做一次全对全比较足以提取第一前沿。所有目标值默认按最小化处理如果你的问题是最大化先在fun里取负。2.3 必调参数NP、F、CR与收敛性的权衡这套模板真正的开关只有几个NP、F、CR、maxGen。它们对结果的影响我总结成下表方便你直接按表试参。参数推荐范围对优化过程的影响我的习惯NP50200太小种群易早熟太大单代计算量显著增长取决策变量维数的10倍至少100F0.30.9控制变异步长大F利于探索小F利于局部精调先用0.5收敛慢再降CR0.60.95交叉概率高子代保留父代信息少多样性变强ZDT类问题固定0.9maxGen100500影响最终收敛程度和运行时间先跑200代看前沿形状F和CR有一个常见的误区不是F越大越容易跳出局部最优。在多目标版本里F过大反而容易让变异后的个体频繁突破边界虽然代码会把边界重置但大量个体堆在边界上会损失内部多样性。如果你发现前沿缺失中间段先检查F是否超过0.7。边界处理这里用的是“越界重置为边界值”实现简单。但我在处理有界工程问题时如果发现种群有聚集在边界上的趋势会改成“越界随机重生”即越界维度重新在lb(j)到ub(j)内随机赋值代价是多写一行效果却更稳。3. 多目标差分进化变体对比与基于分解的MODE实现思路基础MODE能解决多数双目标问题但在目标数增加或问题形态复杂时单靠一种选择机制不够。变体的差异主要体现在选择机制、分解策略和种群维护上。3.1 主流变体的算法结构与适用场景这里列几个MATLAB社区里常见的多目标差分进化变体。变体名称核心机制典型适用场景GDE3每次迭代先扩展种群再用非支配排序和剪枝删除多余个体需要同时处理约束和目标连续问题DEMO每个父代强制生成一个子代合并后用NSGA-II风格选择双目标通用实现简单MOEA/D-DE将多目标分解为多个单目标子问题用邻域进行DE变异目标数23已知偏好或权重NSDE把NSGA-II的非支配排序直接替换DE的选择步骤决策变量多前沿不连续要注意GDE3和DEMO在代码上差别不大DEMO更接近基础MODE而GDE3多了一个“剪枝”阶段。选型上我没有固定答案一般如果问题维度在30以下我会直接用第二节的基础MODE如果需要更均匀的前沿分布就换MOEA/D-DE。3.2 在MATLAB中实现基于分解的MOEA/D-DE变体MOEA/D-DE的思路是用一组权重向量把多目标拆成若干子问题每个子问题对应一个聚合函数值。以Tchebycheff聚合函数为例子问题k的聚合值由权重向量lambda_k和目标值obj计算。更新时只在邻域内选择父代做变异和交叉新解如果聚合值更小就替换邻域内的旧解。下面是核心更新片段完整代码在你把权重向量生成加进去后就能跑% lambda: NP×numObj 的权重向量矩阵每行一个子问题 % B: NP×T 邻域索引矩阵假设每个子问题有T个邻居 % X, Obj: 当前种群和对应目标值 for i 1:NP % 从邻域B(i,:)中随机选择两个不同个体作为r1, r2 nb B(i, :); idx randperm(length(nb), 2); r1 nb(idx(1)); r2 nb(idx(2)); % DE变异和交叉与基础MODE相同 donor X(i, :) F * (X(r1, :) - X(r2, :)); mask rand(1, D) CR; mask(randi(D)) true; u donor .* mask X(i, :) .* (1 - mask); u min(max(u, lb), ub); % 计算新解目标值 newObj fun(u); % 更新邻域内聚合值更差解 for k 1:length(nb) j nb(k); % Tchebycheff聚合值z为理想点即当前各目标最小值 g_old max(lambda(j, :) .* abs(Obj(j, :) - z)); g_new max(lambda(j, :) .* abs(newObj - z)); if g_new g_old X(j, :) u; Obj(j, :) newObj; end end % 更新理想点z z min(z, newObj); end这个片段里最关键的部分是“用聚合值g_new和g_old比较”它直接替代了非支配排序。lambda的分布会决定前沿均匀性我一般用MATLAB自带的lhsdesign生成均匀权重再用knnsearch找每个子问题的邻近索引。注意z是动态更新的理想点如果目标量纲差异大先做归一化再算聚合值。3.3 变体选型根据问题特征选择算法选型没有银弹但可以参考这几个经验。目标只有两个且计算目标函数一次很贵用基础MODEmaxGen控制在200以内不要换分解法。目标前沿不连续或呈碎片状优先试NSDE这类基于非支配的变体因为分解法在碎片前沿上容易丢失部分子问题。参考向量或偏好已知比如工厂里指定“成本不能超过某值”用MOEA/D-DE权重向量可以直接带偏好。决策变量超过100建议先把代码的变异循环向量化否则无论哪种变体都会卡在非支配排序上。我在实际项目里往往先用基础MODE跑通流程验证目标函数无误再用MOEA/D-DE或GDE3去细调。这种顺序能省不少排错时间。4. 在MATLAB中跑通ZDT测试集完整脚本、性能指标与排错有了一套基础MODE下一步是在标准测试函数上看到实际效果。ZDT1是双目标连续凸前沿适合做第一个验证。4.1 定义ZDT1目标函数与算法入口ZDT1的决策变量维度可以任意设但常用30为了快速验证我这里设成10。它有两个目标f1 x(1)g 1 9 * sum(x(2:D))/(D-1)f2 g * (1 - sqrt(f1/g))。决策变量都在[0,1]之间。function f zdt1(x) D length(x); f1 x(1); g 1 9 * sum(x(2:end)) / (D - 1); f2 g * (1 - sqrt(f1/g)); f [f1, f2]; end把这段保存为zdt1.m放在和MODE_basic.m同一目录下。注意fun返回行向量这是第二节代码里的约定。4.2 运行完整脚本从MODE_basic到帕累托前沿绘图下面这个主脚本可以直接执行clc; clear; D 10; lb zeros(1, D); ub ones(1, D); numObj 2; NP 100; maxGen 200; F 0.5; CR 0.9; tic; [pf, pfObj] MODE_basic(zdt1, D, lb, ub, numObj, NP, maxGen, F, CR); toc; % 绘制算法前沿与真实前沿 t linspace(0, 1, 500); trueF [t, 1 - sqrt(t)]; figure; plot(pfObj(:,1), pfObj(:,2), bo, MarkerSize, 4); hold on; plot(trueF(:,1), trueF(:,2), r-, LineWidth, 1.2); xlabel(f1); ylabel(f2); legend(MODE result, True PF); % 计算IGD与HV igd computeIGD(pfObj, trueF); ref [1.1, 1.1]; hv computeHV(pfObj, ref); fprintf(IGD %.4f, HV %.4f\n, igd, hv);computeIGD计算反向世代距离computeHV计算超体积。两者的实现如下可以直接放在主脚本末尾或保存为子函数。function igd computeIGD(pf, truePF) n size(truePF, 1); dist zeros(n, 1); for i 1:n d sqrt(sum((pf - truePF(i, :)).^2, 2)); dist(i) min(d); end igd mean(dist); end function hv computeHV(pf, ref) pf sortrows(pf, 1); hv 0; lastX 0; for i 1:size(pf, 1) x pf(i, 1); if x lastX continue; end width x - lastX; hv hv width * (ref(2) - pf(i, 2)); lastX x; end hv hv (ref(1) - lastX) * ref(2); endcomputeHV只适用于两个目标且前沿按f1单调递增的情况ZDT1恰好满足。如果你换到其他测试函数建议用paretoset或MATLAB File Exchange上的通用HV函数。4.3 计算IGD与HV指标评估收敛性和多样性IGD越小越好它衡量算法前沿与真实前沿的平均距离。HV越大越好参考点通常取略大于各目标最大值。运行上面脚本在F0.5, CR0.9, NP100, maxGen200时我得到IGD约0.01、HV约0.65。不同随机种子会有波动属正常。如果IGD偏大说明前沿离真实前沿远先加大maxGen。如果HV不高说明前沿覆盖不够广泛优先尝试调大NP或降低F。4.4 四个常见运行错误及其解决多目标优化代码最常见的错误几乎都集中在这四处。现象原因解决方向结果全是NaN目标函数中存在除零或log负数检查zdt1里的g当D1时除零前沿只有一个点种群早熟或拥挤距离失效调大NP到150降低F到0.3运行时间暴涨非支配排序O(N^2)NP超过300将NP控制在200以内前沿分布不均边界处理把所有个体压到边界把越界重置改成随机重生第一条里的坑最隐蔽ZDT1当D1时sum(x(2:end))为空g19*0/0就会NaN。我建议在zdt1.m开头加一个断言assert(length(x) 2, ZDT1需要至少2个决策变量);另外MATLAB版本差异也会影响randperm(NP,3)的行为R2023b里它默认返回行向量R2018b也一致。如果你的版本很老可以改成randperm(NP)再取前三个。5. 进阶技巧向量化加速与自适应参数让多目标差分进化更快当问题收敛到前沿形状稳定后瓶颈通常在运行速度。这里给两个我能直接落地的技巧。5.1 向量化变异与交叉去掉内层for循环第二节的代码在NP100时内层循环尚可但NP200时每一代要算200次变异和200次目标函数。目标函数是黑盒时无法加速但变异和交叉可以整段向量化。关键在于一次生成所有变异向量而不逐个for% r1, r2, r3 都是NP×1的随机索引 r1 randperm(NP); r2 randperm(NP); r3 randperm(NP); % 保证互不相同 for i 1:NP while r1(i)r2(i) || r1(i)r3(i) || r2(i)r3(i) r1(i) randi(NP); r2(i) randi(NP); r3(i) randi(NP); end end donors X(r1,:) F * (X(r2,:) - X(r3,:)); mask rand(NP, D) CR; jrand randi(D, NP, 1); mask(sub2ind([NP, D], (1:NP), jrand)) true; U donors .* mask X .* (1 - mask);注意randperm(NP)生成的索引天然无重复但如果NP很小直接对每一维重新随机更安全。向量化后目标函数评估仍然是逐行进行的如果fun支持批量输入矩阵每一行是一个个体还可以再进一步向量化。5.2 自适应F与CR用历史成功率动态调整固定参数能跑通但前沿形状差时往往需要调两次。最简单的一种自适应策略是“中位数更新”每20代统计一下当前代较父代有所改进的个体占比若改进比例高说明搜索方向有效就调大F和CR一点反之调小。if mod(gen, 20) 0 improvement sum(newObjBetter) / NP; if improvement 0.3 F min(F * 1.1, 0.9); CR min(CR * 1.05, 1.0); elseif improvement 0.1 F max(F * 0.9, 0.2); CR max(CR * 0.95, 0.5); end end这个技巧不严谨但在工程上足够它不会让参数突变只做微调。你可以在第三节的MOEA/D-DE变体里也用同样的方法不需要改动选择部分。5.3 如何扩展自己的变体把约束处理加入MODE工程问题几乎都带约束。最常见的做法是把约束违反量作为第三个目标或者用“约束支配”修改2.1节里的dominates函数。我推荐后者改动最小function d dominates(a, b) cv (x) max(x(1), 0); % 简化示意实际约束在外部 % 下面需要比较两个解的约束违反量和目标值 d all(a(1:numObj) b(1:numObj)) any(a(1:numObj) b(1:numObj)); end更好的做法是保存每个解的约束违反总量在dominates里先比较约束违反总量违反小的支配违反大的。这样基础MODE就变成了约束多目标差分进化而变体逻辑完全不变。这个思路在GDE3里被大量使用你只需在fun返回目标值时同时返回约束值再在环境选择里加一个判断即可。本文还有配套的精品资源点击获取