简介这套基于Matlab的NSGA-III多目标优化代码面向需要解决多目标优化问题的科研人员、工程师和算法学习者能帮助快速搜索复杂问题的帕累托前沿。资源共含十个m文件压缩包仅十一KB结构清晰主程序负责整体进化流程算法主体实现非支配排序、参考点生成与环境选择等关键环节另有锦标赛选择、均匀点生成、测试函数等子模块所有文件均可按实际问题修改目标函数与参数直接运行主文件即可看到示例结果。目前已有1165人学习浏览适合作为学习多目标进化算法或开展工程优化的实践起点。运行整套代码后可以直观理解参考点机制、精英保留策略以及多样性与收敛性的平衡方法并借助逆向世代距离指标评估解集质量必要时迁移到资源配置、工程设计、投资组合优化等真实场景中。1. NSGA3 与 NSGA2 的分野为什么第三版不再靠拥挤距离你在三个或更多目标上跑过遗传算法大概率会遇到这样的场面NSGA2 跑出来的帕累托前沿在某一片区域扎堆边界目标始终推不出去。原因是 NSGA2 的拥挤距离在三维及以上目标空间里起不到足够的多样性压力——距离差不多的个体太多删谁留谁几乎随机。NSGA3非支配排序遗传算法第三版用一组均匀分布在超平面上的参考点替代拥挤距离把“均匀”从个体间距转化成对参考方向的覆盖这是它在多目标优化问题中被大量采用的核心原因。这里提到的多目标优化 NSGA3 代码是一份 MATLAB 工程包含 NSGAIII_main.m 入口、NSGAIII.m 主循环、GA.m 遗传算子、UniformPoint.m 参考点生成、NDSort.m 非支配排序、EnvironmentalSelection.m 环境选择以及 funfun.m 和 CalObj.m 构成的目标函数层。完整链路可以跑通 DTLZ、ZDT 等基准问题也方便替换成自己的模型。适合两类人一是刚接触 NSGA3需要边跑代码边理解非支配排序和参考点机制的读者二是已经用过 NSGA2想快速评估第三版算法在自己问题上的效果的工程师。2. 代码结构拆解UniformPoint、NDSort 与 EnvironmentalSelection 的协作压缩包里的文件不算多但每个文件的职责和调用顺序如果只看文件名容易搞反。我拿到的版本中NSGAIII_main.m 是入口它会先设置问题维度和种群规模然后调用 UniformPoint.m 生成参考点再把初始种群交给 NSGAIII.m 进入主循环。主循环里每一代的流程是TournamentSelection.m 从当前种群选父代GA.m 做交叉变异得到子代CalObj.m 计算子代目标值最后 EnvironmentalSelection.m 从父代和子代的合并集合里选出进入下一代的个体。NDSort.m 在这个流程里被调用两次一次在环境选择前做非支配分层一次在指标计算时取第一前沿。下表是我整理的文件功能和关键接口方便在 IDE 里按图索骥。文件名在算法中的角色主要输入输出NSGAIII_main.m入口脚本负责设置 M、D、N、T 等全局参数无返回值输出最终种群和 IGD 曲线NSGAIII.m主循环初始化、进化迭代、记录指标输入种群与参数输出最优前端GA.m遗传操作常见做法是 SBX 交叉加多项式变异输入父代决策变量输出子代决策变量TournamentSelection.m二元锦标赛选择根据排名或目标值挑选父代输入选择压力与适应值输出父代索引NDSort.m非支配排序返回每个个体的帕累托层号输入目标矩阵输出层号和总层数UniformPoint.m在超平面上生成均匀参考点输入目标数 M 与规模 N输出参考点 WEnvironmentalSelection.m基于参考点和小生境计数的环境选择输入合并种群、目标值、参考点输出下一代funfun.m / CalObj.m目标函数和批量计算封装输入决策变量矩阵输出目标矩阵IGD.m计算世代距离指标输入算法前沿与理论前沿输出标量调用链中最容易忽略的是 CalObj 这层封装。很多人会把 funfun 的循环直接写进主循环导致每次评估都要复制一遍代价函数。CalObj 把“决策变量矩阵到目标值矩阵”做成统一入口后续想改并行评估或 GPU 加速都不用动算法主结构。实际调试时我也是先单独跑一行PopObj CalObj(PopDec)确认输出维度是 N x M再去追算法内部效率比一头扎进 NSGAIII.m 里高得多。2.1 UniformPoint.m参考点为什么不能靠随机生成NSGA3 的参考点要求在标准单纯形上均匀分布也就是每个参考点满足所有目标方向的权重之和为 1同时点与点之间尽量等距。如果你用 rand 随机生成点会在单纯形上聚成一团最终环境选择时某些参考方向旁挤满个体另一些方向一个也没有。这就是为什么这套代码要单独用 UniformPoint.m 而不是在 EnvironmentalSelection 里直接随机拉权重。function [W, N] UniformPoint(N, M) % 生成 M 维目标空间中的均匀参考点 % N 是期望的点数M 是目标个数 % 返回 W 是每行一个参考点且满足 sum(W,2)1 H1 1; while nchoosek(H1M-1, M-1) N H1 H1 1; end H1 H1 - 1; W nchoosek(1:H1M-1, M-1) - repmat(0:M-2, nchoosek(H1M-1, M-1), 1); W [W, H1 - sum(W, 2)] / H1; N size(W, 1); end这段代码的要点是 H1 的选取。nchoosek(H1M-1, M-1) 是组合数表示在 M 个目标方向上把每维间隔分成 H1 份后得到的参考点数量。循环先找到第一个不小于期望 N 的 H1 再减 1是为了让实际生成点数量尽量贴近传入的种群规模。组合数结果是一个维度为“参考点数 x (M-1)”的坐标序列最后拼上 H1 - sum(W,2)保证每行求和为 1。实际运行中你会发现 N 不一定等于传入的 N比如三目标时指定 N100可能返回 105 或 91。主循环里应该以返回值 N 为准而不是用写死的 100否则环境选择中参考点与种群规模对不上矩阵维度就会直接报错。2.2 NDSort.m非支配排序在 NSGA3 中仍然不可省环境选择虽然不用拥挤距离但非支配排序仍然在每一代承担“分层剪枝”的工作。环境选择的第一步是把合并后的父代和子代按非支配层排序然后一层一层往下一代里放。当放到某一层时下一代容量可能不够这时才轮到参考点机制来决定这一层里谁留下、谁淘汰。如果没有 NDSort参考点会把搜索方向引向非支配层相同的一堆个体收敛方向会乱掉。function [rank, maxRank] NDSort(PopObj, nSort) % 简化但可读的非支配排序版本 % PopObj: N 行 M 列的目标矩阵假定全部为最小化目标 % nSort: 需要保留的最大层数不传则保留全部 % rank: 每个个体所属的帕累托层号 N size(PopObj, 1); rank ones(N, 1); for i 1:N for j 1:N if i j continue; end if dominates(PopObj(i, :), PopObj(j, :)) rank(j) rank(j) 1; end end end if nargin 1 rank(rank nSort) nSort 1; end maxRank max(rank); end function flag dominates(a, b) % a 支配 b 当且仅当 a 所有目标不差于 b 且有至少一个目标更优 flag all(a b) any(a b); end这段代码是 O(N^2) 的教学版实际 NDSort.m 会按被支配计数和支配集合并的方式把复杂度压到接近 O(N log N) 的数量级。需要留神的是MATLAB 里如果 PopObj 存在数值一样的个体dominates 的判断会变成 all(ab) 成立但 any(ab) 不成立两个完全相同的个体互不支配它们会落在同一层。这种等值个体在工程问题里很常见你可以在进入 NDSort 前给目标矩阵加一个机器精度级别的抖动或者修改比较条件为 all(a b 1e-10)。2.3 EnvironmentalSelection.m关联参考线与小生境计数的选择压环境选择是 NSGA3 和 NSGA2 差异最集中的地方。拿到临界层个体后算法先把所有目标值归一化到参考点所在的超平面上然后计算每个个体到每条参考线的垂直距离把个体关联到距离最近的参考点。接着统计每个参考点当前已关联的个体数量这个数量叫小生境计数用小写 rho 表示。选择时优先保留 rho 最小的参考点方向上的个体从而让种群从拥挤方向往稀疏方向迁移。% EnvironmentalSelection.m 里的关键步骤示意 PopObjNorm (PopObj - repmat(zmin, size(PopObj, 1), 1)) ... ./ repmat(zmax - zmin, size(PopObj, 1), 1); Distance pdist2(PopObjNorm, W, cosine); [~, assoc] min(Distance, [], 2); rho accumarray(assoc, 1, [size(W, 1), 1]); % 对不满足数量要求的方向重新选择 rho 最小的参考点这里用余弦距离还是欧氏距离并不是关键关键是要把目标点映射到以理想点 zmin 为原点的空间。如果 zmax 和 zmin 接近归一化会出现放大噪声可以在这行之前打印 zmax-zmin 看一下尺度。我曾经遇到一个实际模型的两个目标在某一代几乎不变zmax-zmin 掉到 1e-6结果 EnvironmentalSelection 的关联结果等于随机分配整代进化直接退化。后来在 CalObj 里给目标做了 log 或 sqrt 变换情况才稳定下来。另一个坑是 pdist2 需要 Statistics Toolbox如果你没有这个工具箱可以用两层循环加 sqrt(sum((A-B).^2)) 实现多目标规模下性能差距可以接受。3. 把 NSGAIII_main 跑起来参数设置、目标函数与第一张帕累托前沿第 2 章把文件关系理清之后下一步是让它在你的 MATLAB 环境里真的跑起来。这套代码不依赖优化工具箱只要 nchoosek、pdist2 这些基础能力可用就能运行。我建议先不要动 NSGAIII.m 内部先从入口脚本改起确定目标个数 M、决策变量维度 D、种群规模 N 和最大进化代数 T把 funfun.m 替换成要优化的函数然后运行 NSGAIII_main.m。3.1 修改 funfun 和 CalObj测试函数的最小改法funfun.m 的约定是接收一个 N 行 D 列的决策变量矩阵返回一个 N 行 M 列的目标值矩阵。很多从 NSGA2 代码转过来的用户喜欢在 funfun 里写 for 循环逐行计算这在 N 只有 50 时没问题但 NSGA3 的种群经常是 100 到 300还要跑 500 代循环代价会被放大。最好的做法是把计算写成向量形式至少要把标量计算改成用 x(:,k) 这样的列切片。% funfun.m - 以 ZDT1 为例的两个目标测试函数 function f funfun(x) % x: N x D 决策变量矩阵每一行为一个个体 % f: N x M 目标矩阵这里 M2全部为最小化 n size(x, 2); g 1 9 * sum(x(:, 2:end), 2) / (n - 1); f(:, 1) x(:, 1); f(:, 2) g .* (1 - sqrt(x(:, 1) ./ g)); end这段代码中 g 是 ZDT1 的距离函数x(:,1) 控制分布x(:,2:end) 控制收敛。sum(x(:,2:end),2) 是沿行求和结果为 N 维列向量9/(n-1) 是标量两者做乘除后 g 仍保持列向量。f 的第一列固定取第一维决策变量第二列用 g 做缩放保证真实前沿是凸的。替换为实际模型时最容易出的错是行方向没对齐比如直接用 x(:,1) 和某个标量内部变量相乘结果变成 Nx1 和 1xN 的隐式扩展最后 f 的维度变成 NxN。遇到维度报错第一步先 size 检查 x 和中间结果。CalObj.m 在这套代码里通常只做一件事调用 funfun 并处理极端值。如果你的实际问题有最大化的目标不要写成长率最大化然后指望算法自动处理NSGA3 默认按最小化设计统一在 funfun 里把最大化目标取负后面所有指标和绘图才不会出现反向结果。3.2 GA.m 中的交叉变异参数SBX 和多项式变异的典型取值遗传操作在 NSGA3 里一般沿用 NSGA2 的配置模拟二进制交叉 SBX 负责在父代附近生成子代多项式变异负责局部扰动。交叉概率一般取 1.0表示所有配对的父代都要参与交叉分布指数 disC 通常取 20数值越大子代越容易靠近父代变异概率 proM 取 1/D即每个决策变量平均有一个发生变异分布指数 disM 也取 20。% GA.m 中的参数设置和调用 proC 1; disC 20; proM 1 / size(ParentDec, 2); disM 20; [OffDec1, OffDec2] SBX(ParentDec(1:2:end, :), ParentDec(2:2:end, :), lower, upper, proC, disC); OffDec [OffDec1; OffDec2]; OffDec PolynomialMutation(OffDec, lower, upper, proM, disM);这里 ParentDec 是锦标赛选择后的父代按邻对方式配对。如果种群 N 是奇数2:2:end 最后会空出来主循环里最好先判断 N 是否为偶数或者在初始化时把 N 改成偶数。实际测试中我发现决策变量维度 D 很大时 proM1/D 会变得很小变异几乎不起作用。常见做法是把变异概率的下限设为 0.05或者对当前代数做一个自适应映射让算法在后期更容易跳出局部前沿。这些不属于 NSGA3 框架的强制要求但能明显改善真实工程的收敛速度。3.3 参数表与主循环顺序下面这张参数表可以作为第一批试验的基线。NSGA3 相比 NSGA2 多了一个参考点规模的问题因此 N 不再是一个随意指定的偶数而是由 UniformPoint 返回的点数决定。第一次跑先按表中取值确认 IGD 能下降后再动 N 和 T最后再调 disC 和 disM避免多个参数同时改导致不知道是谁在起作用。参数第一次跑的取值说明M 目标数2 或 3从 M2 开始排查代码再切三目标D 决策变量维数10 到 30ZDT1 常用 30DTLZ 系列常用 12 或 15N 种群规模由 UniformPoint(N,M) 返回的实际点数三目标时常用 91 或 105T 最大代数500benchmark 可跑到 1000工程问题看 IGD 平台proC / disC1.0 / 20交叉概率低时种群多样性不足proM / disM1/D / 20变异率太小容易早熟主循环的顺序不要随意颠倒先用 UniformPoint 拿到参考点 W再用 lower 和 upper 生成初始种群然后才进入 NSGAIII.m 的 for 循环。如果先初始化再调用 UniformPoint可能出现初始种群的规模与参考点数量不一致环境选择会报维度错误。% NSGAIII_main.m 的关键片段 M 2; D 30; N 100; T 500; lower zeros(1, D); upper ones(1, D); [W, N] UniformPoint(N, M); PopDec lower rand(N, D) .* (upper - lower); PopObj CalObj(PopDec); % 之后把 PopDec, PopObj, W 交给 NSGAIII.m 迭代生成初始种群时用 rand(N,D) 后乘上 upper-lower 再加 lower得到的是 [lower, upper) 区间内均匀分布的实数值。如果你的决策变量是整数或离散的直接对连续的进化结果做 round 会引入重复个体更好的做法是在 funfun 内部对决策变量取整这样遗传操作仍保持连续性但目标函数反映的是离散模型。这个技巧对排产、组合优化很实用。3.4 画帕累托前沿二维与三维跑完一代或全部代数后最直接的检查是画图。先取第一非支配层个体再根据 M 选择 plot 或 scatter3。这一步能很快暴露目标函数写错、方向反向、前沿边界缺失等问题。画图用的数据要直接从当前种群取不要用历史存档因为后者可能包含旧代数的被支配个体。% 取第一非支配层 [rank, ~] NDSort(PopObj); Front PopObj(rank 1, :); if M 2 plot(Front(:, 1), Front(:, 2), o, MarkerSize, 4); elseif M 3 scatter3(Front(:, 1), Front(:, 2), Front(:, 3), .); end画图时要注意 rank 变量没有排序意义rank1 的就是当前最优前沿。如果你看到前沿上某一段明显缺块先别急着调参考点去检查 funfun 是不是在局部区域返回了 NaN 或 Inf。NaN 进入 NDSort 后比较运算会全部返回 false导致这一层里大量个体无法排序。4. 把 funfun 换成自己的问题决策变量、约束与归一化跑通基准测试只是第一步实际工程问题里很少有现成的 ZDT 函数。替换自定义模型时主要的改动集中在 funfun.m、入口脚本的 lower/upper以及约束的处理。这一章我会用一个双目标示例走一遍替换流程并标出 NSGA3 里最容易因为归一化而翻车的位置。4.1 双目标真实模型的目标函数写法假设你要优化一个简单的生产方案决策变量有两个原材料配比 x1 和加工温度 x2。目标一是材料成本随配比上升而增加目标二是产出率随温度和配比有一个非线性关系。NSGA3 里全部按最小化处理产出率取负即可。function f funfun(x) % x: N x 2 cost 2 * x(:, 1) 5 * x(:, 2); yield -(x(:, 1).^1.5 .* exp(-x(:, 2) ./ 3)); f [cost, yield]; end这里 cost 越小越好yield 是更大产出更好所以负号让它在排序时与 cost 同向。变量范围在 main 脚本里定义lower[0, 0]upper[1, 100] 之类。要注意的是x 两维的量纲差异非常大0-1 和 0-100SBX 交叉和多项式变异在决策空间里是按同样步长生成的温度维的微小变化会被 x1 的尺度掩盖。常见做法是在 funfun 内部对 x 先做归一化再换算真实值或者把 upper 设成同等量纲优化完再映射回真实温度。我一般倾向于前者因为不会影响遗传算子的比例。4.2 约束处理罚函数还是可行优先标准 NSGA3 选择机制没有显式的约束处理罚函数是最容易落地的方式。实现时把约束违反量变成目标函数里的惩罚项让不可行解在非支配排序时被支配或落到较后的层。% 在 funfun 的末尾加约束惩罚 g1 x(:, 2) - 80; % 约束温度不能超过 80g10 viol max(0, g1); f f 1e6 * viol; % 对两个目标同时加惩罚这里的 1e6 不是固定值原则是让不可行解的任意一个目标值都显著劣于可行解。若两个目标的正常量级分别是 10 和 1001e6 会过大可能导致目标矩阵数值范围太宽归一化后可行区域被压缩到看不见。可以先跑一次不约束的版本统计目标正常范围的量级再取量级的 100 到 1000 倍作为惩罚系数。多约束时把每个约束的 viol 拼成一个矩阵最后取 max 或加权和均可。罚函数的一个弱点是会让大量不可行解进入初始种群环境选择时不保证这些不可行解一定消失。另一种做法是可行优先在种群初始化时逐个判断约束只保留可行个体直到填满 N。缺点是可行域很小时初始化几乎死循环。对于工程模型我建议先用罚函数跑通流程观察可行比例如果可行比例长时间低于 10%再考虑给算法外层包一个修复算子或专门处理约束的变体。4.3 归一化与 Nadir 点的调试为什么参考点总是偏向某些区域NSGA3 的归一化发生在 EnvironmentalSelection 中目的是让不同尺度的目标值都映射到参考点所在的超平面。大多数实现会先找每一维目标的理想点 zmin再用当前非支配层中每一维的最大值来估计 Nadir 点 zmax。这个估计并不总是准确尤其是前沿不完整时zmax 可能来自一个极端个体导致参考点整体被拉向某个方向。症状可能原因检查方法前沿偏向目标 1 很小的区域zmax 中 f2 维估计过大打印 rank1 个体的 max(PopObj) 和 min(PopObj)某代以后参考点关联几乎都指向同一个点zmax-zmin 接近 0在 EnvironmentalSelection 入口计算数值范围两个目标尺度差 1e6 时前沿散成一片归一化前没有做 log/sqrt 变换在 funfun 中对目标做单调变换% 调试归一化在 EnvironmentalSelection 之前打印尺度 zmin min(PopObj, [], 1); zmax max(PopObj, [], 1); fprintf(range: %e %e\n, zmax - zmin);如果看到 range 中有 1e-10 级别的分量说明这一维目标在所有个体上几乎没有区分度那不只是归一化问题而是目标函数设计问题。这时先去检查 funfun 的输入输出有没有写反列或者某一维决策变量根本没参与计算。把这一问题修掉后再回来调参考点数量才有意义。我之前在一个工程模型里花了一整个下午调 UniformPoint 的 H 值最后发现是 funfun 第二列目标永远等于常数所有参考点关联都在第一列方向这属于最容易踩但最不容易想到的坑。5. 收尾章用 IGD 验证参考点均匀性和选择压力5.1 计算 IGD 时容易忽略的两个细节IGD 是衡量算法得到的前沿与理论前沿平均距离的指标数值越小越好。IGD.m 的核心只有两行计算理论前沿到算法前沿的最小距离再取平均。虽然代码量少但两个使用细节会影响你判断的准确性尤其是前沿中混入被支配解时IGD 会被这些点拉得异常低。% IGD.m function score IGD(PF, PFtrue) if isempty(PF) score inf; return; end d pdist2(PFtrue, PF); score mean(min(d, [], 2)); end第一个细节是在调用前先把 PF 里的支配点去掉否则 pdist2 会对每个理论点找到最近的算法点但如果算法前沿里有大量被支配点堆在一个角落它们不仅不会改善 IGD还会掩盖真实前沿的不均匀。第二个细节是如果理论前沿 PFtrue 采样点太少mean(min(...)) 会低估误差DTLZ1 的理论前沿需要至少 10000 个采样点ZDT1 的前沿是连续曲线常用 500 到 1000 个点已够。5.2 用 IGD 曲线判断参考点分布是否合理在 NSGAIII_main.m 的循环里记录每一代 IGD最后用 semilogy 画出来比看最终数字有用得多。history zeros(1, T); for t 1:T % ... 一次进化 ... [rank, ~] NDSort(PopObj); Front PopObj(rank 1, :); history(t) IGD(Front, PFtrue); end semilogy(1:T, history); xlabel(Generation); ylabel(IGD);曲线如果快速下降后进入平台说明收敛和多样性已经稳定继续增加代数收益不大如果曲线长期不下调小 disM 或增大变异概率比增加代数更有效如果曲线先降后升说明种群多样性崩溃常见原因是参考点数量太少或归一化失败。一个具体技巧是同时运行 M2 和 M3 的两份配置用同一把随机种子观察三目标的平台是否明显高于双目标。若是优先怀疑参考点分布而不是遗传算法参数。5.3 一个验证 UniformPoint 的三行检查最后回到参考点本身。无论替换成什么模型先单独验证 UniformPoint 的输出[W, N] UniformPoint(100, 3); assert(all(abs(sum(W, 2) - 1) 1e-6)); scatter3(W(:, 1), W(:, 2), W(:, 3), filled);看到点均匀铺在一个三角形内说明参考点生成正常再回到 NSGAIII_main.m 调种群规模。若点出现聚集或者有一侧明显稀疏检查 nchoosek 输入参数是否在你的 MATLAB 版本中被当作 double 导致精度丢失必要时改用 nchoosek(int32(H1M-1), int32(M-1)) 重新生成。本文还有配套的精品资源点击获取
NSGA3多目标优化算法详解：基于参考点的MATLAB实现与代码解析
简介这套基于Matlab的NSGA-III多目标优化代码面向需要解决多目标优化问题的科研人员、工程师和算法学习者能帮助快速搜索复杂问题的帕累托前沿。资源共含十个m文件压缩包仅十一KB结构清晰主程序负责整体进化流程算法主体实现非支配排序、参考点生成与环境选择等关键环节另有锦标赛选择、均匀点生成、测试函数等子模块所有文件均可按实际问题修改目标函数与参数直接运行主文件即可看到示例结果。目前已有1165人学习浏览适合作为学习多目标进化算法或开展工程优化的实践起点。运行整套代码后可以直观理解参考点机制、精英保留策略以及多样性与收敛性的平衡方法并借助逆向世代距离指标评估解集质量必要时迁移到资源配置、工程设计、投资组合优化等真实场景中。1. NSGA3 与 NSGA2 的分野为什么第三版不再靠拥挤距离你在三个或更多目标上跑过遗传算法大概率会遇到这样的场面NSGA2 跑出来的帕累托前沿在某一片区域扎堆边界目标始终推不出去。原因是 NSGA2 的拥挤距离在三维及以上目标空间里起不到足够的多样性压力——距离差不多的个体太多删谁留谁几乎随机。NSGA3非支配排序遗传算法第三版用一组均匀分布在超平面上的参考点替代拥挤距离把“均匀”从个体间距转化成对参考方向的覆盖这是它在多目标优化问题中被大量采用的核心原因。这里提到的多目标优化 NSGA3 代码是一份 MATLAB 工程包含 NSGAIII_main.m 入口、NSGAIII.m 主循环、GA.m 遗传算子、UniformPoint.m 参考点生成、NDSort.m 非支配排序、EnvironmentalSelection.m 环境选择以及 funfun.m 和 CalObj.m 构成的目标函数层。完整链路可以跑通 DTLZ、ZDT 等基准问题也方便替换成自己的模型。适合两类人一是刚接触 NSGA3需要边跑代码边理解非支配排序和参考点机制的读者二是已经用过 NSGA2想快速评估第三版算法在自己问题上的效果的工程师。2. 代码结构拆解UniformPoint、NDSort 与 EnvironmentalSelection 的协作压缩包里的文件不算多但每个文件的职责和调用顺序如果只看文件名容易搞反。我拿到的版本中NSGAIII_main.m 是入口它会先设置问题维度和种群规模然后调用 UniformPoint.m 生成参考点再把初始种群交给 NSGAIII.m 进入主循环。主循环里每一代的流程是TournamentSelection.m 从当前种群选父代GA.m 做交叉变异得到子代CalObj.m 计算子代目标值最后 EnvironmentalSelection.m 从父代和子代的合并集合里选出进入下一代的个体。NDSort.m 在这个流程里被调用两次一次在环境选择前做非支配分层一次在指标计算时取第一前沿。下表是我整理的文件功能和关键接口方便在 IDE 里按图索骥。文件名在算法中的角色主要输入输出NSGAIII_main.m入口脚本负责设置 M、D、N、T 等全局参数无返回值输出最终种群和 IGD 曲线NSGAIII.m主循环初始化、进化迭代、记录指标输入种群与参数输出最优前端GA.m遗传操作常见做法是 SBX 交叉加多项式变异输入父代决策变量输出子代决策变量TournamentSelection.m二元锦标赛选择根据排名或目标值挑选父代输入选择压力与适应值输出父代索引NDSort.m非支配排序返回每个个体的帕累托层号输入目标矩阵输出层号和总层数UniformPoint.m在超平面上生成均匀参考点输入目标数 M 与规模 N输出参考点 WEnvironmentalSelection.m基于参考点和小生境计数的环境选择输入合并种群、目标值、参考点输出下一代funfun.m / CalObj.m目标函数和批量计算封装输入决策变量矩阵输出目标矩阵IGD.m计算世代距离指标输入算法前沿与理论前沿输出标量调用链中最容易忽略的是 CalObj 这层封装。很多人会把 funfun 的循环直接写进主循环导致每次评估都要复制一遍代价函数。CalObj 把“决策变量矩阵到目标值矩阵”做成统一入口后续想改并行评估或 GPU 加速都不用动算法主结构。实际调试时我也是先单独跑一行PopObj CalObj(PopDec)确认输出维度是 N x M再去追算法内部效率比一头扎进 NSGAIII.m 里高得多。2.1 UniformPoint.m参考点为什么不能靠随机生成NSGA3 的参考点要求在标准单纯形上均匀分布也就是每个参考点满足所有目标方向的权重之和为 1同时点与点之间尽量等距。如果你用 rand 随机生成点会在单纯形上聚成一团最终环境选择时某些参考方向旁挤满个体另一些方向一个也没有。这就是为什么这套代码要单独用 UniformPoint.m 而不是在 EnvironmentalSelection 里直接随机拉权重。function [W, N] UniformPoint(N, M) % 生成 M 维目标空间中的均匀参考点 % N 是期望的点数M 是目标个数 % 返回 W 是每行一个参考点且满足 sum(W,2)1 H1 1; while nchoosek(H1M-1, M-1) N H1 H1 1; end H1 H1 - 1; W nchoosek(1:H1M-1, M-1) - repmat(0:M-2, nchoosek(H1M-1, M-1), 1); W [W, H1 - sum(W, 2)] / H1; N size(W, 1); end这段代码的要点是 H1 的选取。nchoosek(H1M-1, M-1) 是组合数表示在 M 个目标方向上把每维间隔分成 H1 份后得到的参考点数量。循环先找到第一个不小于期望 N 的 H1 再减 1是为了让实际生成点数量尽量贴近传入的种群规模。组合数结果是一个维度为“参考点数 x (M-1)”的坐标序列最后拼上 H1 - sum(W,2)保证每行求和为 1。实际运行中你会发现 N 不一定等于传入的 N比如三目标时指定 N100可能返回 105 或 91。主循环里应该以返回值 N 为准而不是用写死的 100否则环境选择中参考点与种群规模对不上矩阵维度就会直接报错。2.2 NDSort.m非支配排序在 NSGA3 中仍然不可省环境选择虽然不用拥挤距离但非支配排序仍然在每一代承担“分层剪枝”的工作。环境选择的第一步是把合并后的父代和子代按非支配层排序然后一层一层往下一代里放。当放到某一层时下一代容量可能不够这时才轮到参考点机制来决定这一层里谁留下、谁淘汰。如果没有 NDSort参考点会把搜索方向引向非支配层相同的一堆个体收敛方向会乱掉。function [rank, maxRank] NDSort(PopObj, nSort) % 简化但可读的非支配排序版本 % PopObj: N 行 M 列的目标矩阵假定全部为最小化目标 % nSort: 需要保留的最大层数不传则保留全部 % rank: 每个个体所属的帕累托层号 N size(PopObj, 1); rank ones(N, 1); for i 1:N for j 1:N if i j continue; end if dominates(PopObj(i, :), PopObj(j, :)) rank(j) rank(j) 1; end end end if nargin 1 rank(rank nSort) nSort 1; end maxRank max(rank); end function flag dominates(a, b) % a 支配 b 当且仅当 a 所有目标不差于 b 且有至少一个目标更优 flag all(a b) any(a b); end这段代码是 O(N^2) 的教学版实际 NDSort.m 会按被支配计数和支配集合并的方式把复杂度压到接近 O(N log N) 的数量级。需要留神的是MATLAB 里如果 PopObj 存在数值一样的个体dominates 的判断会变成 all(ab) 成立但 any(ab) 不成立两个完全相同的个体互不支配它们会落在同一层。这种等值个体在工程问题里很常见你可以在进入 NDSort 前给目标矩阵加一个机器精度级别的抖动或者修改比较条件为 all(a b 1e-10)。2.3 EnvironmentalSelection.m关联参考线与小生境计数的选择压环境选择是 NSGA3 和 NSGA2 差异最集中的地方。拿到临界层个体后算法先把所有目标值归一化到参考点所在的超平面上然后计算每个个体到每条参考线的垂直距离把个体关联到距离最近的参考点。接着统计每个参考点当前已关联的个体数量这个数量叫小生境计数用小写 rho 表示。选择时优先保留 rho 最小的参考点方向上的个体从而让种群从拥挤方向往稀疏方向迁移。% EnvironmentalSelection.m 里的关键步骤示意 PopObjNorm (PopObj - repmat(zmin, size(PopObj, 1), 1)) ... ./ repmat(zmax - zmin, size(PopObj, 1), 1); Distance pdist2(PopObjNorm, W, cosine); [~, assoc] min(Distance, [], 2); rho accumarray(assoc, 1, [size(W, 1), 1]); % 对不满足数量要求的方向重新选择 rho 最小的参考点这里用余弦距离还是欧氏距离并不是关键关键是要把目标点映射到以理想点 zmin 为原点的空间。如果 zmax 和 zmin 接近归一化会出现放大噪声可以在这行之前打印 zmax-zmin 看一下尺度。我曾经遇到一个实际模型的两个目标在某一代几乎不变zmax-zmin 掉到 1e-6结果 EnvironmentalSelection 的关联结果等于随机分配整代进化直接退化。后来在 CalObj 里给目标做了 log 或 sqrt 变换情况才稳定下来。另一个坑是 pdist2 需要 Statistics Toolbox如果你没有这个工具箱可以用两层循环加 sqrt(sum((A-B).^2)) 实现多目标规模下性能差距可以接受。3. 把 NSGAIII_main 跑起来参数设置、目标函数与第一张帕累托前沿第 2 章把文件关系理清之后下一步是让它在你的 MATLAB 环境里真的跑起来。这套代码不依赖优化工具箱只要 nchoosek、pdist2 这些基础能力可用就能运行。我建议先不要动 NSGAIII.m 内部先从入口脚本改起确定目标个数 M、决策变量维度 D、种群规模 N 和最大进化代数 T把 funfun.m 替换成要优化的函数然后运行 NSGAIII_main.m。3.1 修改 funfun 和 CalObj测试函数的最小改法funfun.m 的约定是接收一个 N 行 D 列的决策变量矩阵返回一个 N 行 M 列的目标值矩阵。很多从 NSGA2 代码转过来的用户喜欢在 funfun 里写 for 循环逐行计算这在 N 只有 50 时没问题但 NSGA3 的种群经常是 100 到 300还要跑 500 代循环代价会被放大。最好的做法是把计算写成向量形式至少要把标量计算改成用 x(:,k) 这样的列切片。% funfun.m - 以 ZDT1 为例的两个目标测试函数 function f funfun(x) % x: N x D 决策变量矩阵每一行为一个个体 % f: N x M 目标矩阵这里 M2全部为最小化 n size(x, 2); g 1 9 * sum(x(:, 2:end), 2) / (n - 1); f(:, 1) x(:, 1); f(:, 2) g .* (1 - sqrt(x(:, 1) ./ g)); end这段代码中 g 是 ZDT1 的距离函数x(:,1) 控制分布x(:,2:end) 控制收敛。sum(x(:,2:end),2) 是沿行求和结果为 N 维列向量9/(n-1) 是标量两者做乘除后 g 仍保持列向量。f 的第一列固定取第一维决策变量第二列用 g 做缩放保证真实前沿是凸的。替换为实际模型时最容易出的错是行方向没对齐比如直接用 x(:,1) 和某个标量内部变量相乘结果变成 Nx1 和 1xN 的隐式扩展最后 f 的维度变成 NxN。遇到维度报错第一步先 size 检查 x 和中间结果。CalObj.m 在这套代码里通常只做一件事调用 funfun 并处理极端值。如果你的实际问题有最大化的目标不要写成长率最大化然后指望算法自动处理NSGA3 默认按最小化设计统一在 funfun 里把最大化目标取负后面所有指标和绘图才不会出现反向结果。3.2 GA.m 中的交叉变异参数SBX 和多项式变异的典型取值遗传操作在 NSGA3 里一般沿用 NSGA2 的配置模拟二进制交叉 SBX 负责在父代附近生成子代多项式变异负责局部扰动。交叉概率一般取 1.0表示所有配对的父代都要参与交叉分布指数 disC 通常取 20数值越大子代越容易靠近父代变异概率 proM 取 1/D即每个决策变量平均有一个发生变异分布指数 disM 也取 20。% GA.m 中的参数设置和调用 proC 1; disC 20; proM 1 / size(ParentDec, 2); disM 20; [OffDec1, OffDec2] SBX(ParentDec(1:2:end, :), ParentDec(2:2:end, :), lower, upper, proC, disC); OffDec [OffDec1; OffDec2]; OffDec PolynomialMutation(OffDec, lower, upper, proM, disM);这里 ParentDec 是锦标赛选择后的父代按邻对方式配对。如果种群 N 是奇数2:2:end 最后会空出来主循环里最好先判断 N 是否为偶数或者在初始化时把 N 改成偶数。实际测试中我发现决策变量维度 D 很大时 proM1/D 会变得很小变异几乎不起作用。常见做法是把变异概率的下限设为 0.05或者对当前代数做一个自适应映射让算法在后期更容易跳出局部前沿。这些不属于 NSGA3 框架的强制要求但能明显改善真实工程的收敛速度。3.3 参数表与主循环顺序下面这张参数表可以作为第一批试验的基线。NSGA3 相比 NSGA2 多了一个参考点规模的问题因此 N 不再是一个随意指定的偶数而是由 UniformPoint 返回的点数决定。第一次跑先按表中取值确认 IGD 能下降后再动 N 和 T最后再调 disC 和 disM避免多个参数同时改导致不知道是谁在起作用。参数第一次跑的取值说明M 目标数2 或 3从 M2 开始排查代码再切三目标D 决策变量维数10 到 30ZDT1 常用 30DTLZ 系列常用 12 或 15N 种群规模由 UniformPoint(N,M) 返回的实际点数三目标时常用 91 或 105T 最大代数500benchmark 可跑到 1000工程问题看 IGD 平台proC / disC1.0 / 20交叉概率低时种群多样性不足proM / disM1/D / 20变异率太小容易早熟主循环的顺序不要随意颠倒先用 UniformPoint 拿到参考点 W再用 lower 和 upper 生成初始种群然后才进入 NSGAIII.m 的 for 循环。如果先初始化再调用 UniformPoint可能出现初始种群的规模与参考点数量不一致环境选择会报维度错误。% NSGAIII_main.m 的关键片段 M 2; D 30; N 100; T 500; lower zeros(1, D); upper ones(1, D); [W, N] UniformPoint(N, M); PopDec lower rand(N, D) .* (upper - lower); PopObj CalObj(PopDec); % 之后把 PopDec, PopObj, W 交给 NSGAIII.m 迭代生成初始种群时用 rand(N,D) 后乘上 upper-lower 再加 lower得到的是 [lower, upper) 区间内均匀分布的实数值。如果你的决策变量是整数或离散的直接对连续的进化结果做 round 会引入重复个体更好的做法是在 funfun 内部对决策变量取整这样遗传操作仍保持连续性但目标函数反映的是离散模型。这个技巧对排产、组合优化很实用。3.4 画帕累托前沿二维与三维跑完一代或全部代数后最直接的检查是画图。先取第一非支配层个体再根据 M 选择 plot 或 scatter3。这一步能很快暴露目标函数写错、方向反向、前沿边界缺失等问题。画图用的数据要直接从当前种群取不要用历史存档因为后者可能包含旧代数的被支配个体。% 取第一非支配层 [rank, ~] NDSort(PopObj); Front PopObj(rank 1, :); if M 2 plot(Front(:, 1), Front(:, 2), o, MarkerSize, 4); elseif M 3 scatter3(Front(:, 1), Front(:, 2), Front(:, 3), .); end画图时要注意 rank 变量没有排序意义rank1 的就是当前最优前沿。如果你看到前沿上某一段明显缺块先别急着调参考点去检查 funfun 是不是在局部区域返回了 NaN 或 Inf。NaN 进入 NDSort 后比较运算会全部返回 false导致这一层里大量个体无法排序。4. 把 funfun 换成自己的问题决策变量、约束与归一化跑通基准测试只是第一步实际工程问题里很少有现成的 ZDT 函数。替换自定义模型时主要的改动集中在 funfun.m、入口脚本的 lower/upper以及约束的处理。这一章我会用一个双目标示例走一遍替换流程并标出 NSGA3 里最容易因为归一化而翻车的位置。4.1 双目标真实模型的目标函数写法假设你要优化一个简单的生产方案决策变量有两个原材料配比 x1 和加工温度 x2。目标一是材料成本随配比上升而增加目标二是产出率随温度和配比有一个非线性关系。NSGA3 里全部按最小化处理产出率取负即可。function f funfun(x) % x: N x 2 cost 2 * x(:, 1) 5 * x(:, 2); yield -(x(:, 1).^1.5 .* exp(-x(:, 2) ./ 3)); f [cost, yield]; end这里 cost 越小越好yield 是更大产出更好所以负号让它在排序时与 cost 同向。变量范围在 main 脚本里定义lower[0, 0]upper[1, 100] 之类。要注意的是x 两维的量纲差异非常大0-1 和 0-100SBX 交叉和多项式变异在决策空间里是按同样步长生成的温度维的微小变化会被 x1 的尺度掩盖。常见做法是在 funfun 内部对 x 先做归一化再换算真实值或者把 upper 设成同等量纲优化完再映射回真实温度。我一般倾向于前者因为不会影响遗传算子的比例。4.2 约束处理罚函数还是可行优先标准 NSGA3 选择机制没有显式的约束处理罚函数是最容易落地的方式。实现时把约束违反量变成目标函数里的惩罚项让不可行解在非支配排序时被支配或落到较后的层。% 在 funfun 的末尾加约束惩罚 g1 x(:, 2) - 80; % 约束温度不能超过 80g10 viol max(0, g1); f f 1e6 * viol; % 对两个目标同时加惩罚这里的 1e6 不是固定值原则是让不可行解的任意一个目标值都显著劣于可行解。若两个目标的正常量级分别是 10 和 1001e6 会过大可能导致目标矩阵数值范围太宽归一化后可行区域被压缩到看不见。可以先跑一次不约束的版本统计目标正常范围的量级再取量级的 100 到 1000 倍作为惩罚系数。多约束时把每个约束的 viol 拼成一个矩阵最后取 max 或加权和均可。罚函数的一个弱点是会让大量不可行解进入初始种群环境选择时不保证这些不可行解一定消失。另一种做法是可行优先在种群初始化时逐个判断约束只保留可行个体直到填满 N。缺点是可行域很小时初始化几乎死循环。对于工程模型我建议先用罚函数跑通流程观察可行比例如果可行比例长时间低于 10%再考虑给算法外层包一个修复算子或专门处理约束的变体。4.3 归一化与 Nadir 点的调试为什么参考点总是偏向某些区域NSGA3 的归一化发生在 EnvironmentalSelection 中目的是让不同尺度的目标值都映射到参考点所在的超平面。大多数实现会先找每一维目标的理想点 zmin再用当前非支配层中每一维的最大值来估计 Nadir 点 zmax。这个估计并不总是准确尤其是前沿不完整时zmax 可能来自一个极端个体导致参考点整体被拉向某个方向。症状可能原因检查方法前沿偏向目标 1 很小的区域zmax 中 f2 维估计过大打印 rank1 个体的 max(PopObj) 和 min(PopObj)某代以后参考点关联几乎都指向同一个点zmax-zmin 接近 0在 EnvironmentalSelection 入口计算数值范围两个目标尺度差 1e6 时前沿散成一片归一化前没有做 log/sqrt 变换在 funfun 中对目标做单调变换% 调试归一化在 EnvironmentalSelection 之前打印尺度 zmin min(PopObj, [], 1); zmax max(PopObj, [], 1); fprintf(range: %e %e\n, zmax - zmin);如果看到 range 中有 1e-10 级别的分量说明这一维目标在所有个体上几乎没有区分度那不只是归一化问题而是目标函数设计问题。这时先去检查 funfun 的输入输出有没有写反列或者某一维决策变量根本没参与计算。把这一问题修掉后再回来调参考点数量才有意义。我之前在一个工程模型里花了一整个下午调 UniformPoint 的 H 值最后发现是 funfun 第二列目标永远等于常数所有参考点关联都在第一列方向这属于最容易踩但最不容易想到的坑。5. 收尾章用 IGD 验证参考点均匀性和选择压力5.1 计算 IGD 时容易忽略的两个细节IGD 是衡量算法得到的前沿与理论前沿平均距离的指标数值越小越好。IGD.m 的核心只有两行计算理论前沿到算法前沿的最小距离再取平均。虽然代码量少但两个使用细节会影响你判断的准确性尤其是前沿中混入被支配解时IGD 会被这些点拉得异常低。% IGD.m function score IGD(PF, PFtrue) if isempty(PF) score inf; return; end d pdist2(PFtrue, PF); score mean(min(d, [], 2)); end第一个细节是在调用前先把 PF 里的支配点去掉否则 pdist2 会对每个理论点找到最近的算法点但如果算法前沿里有大量被支配点堆在一个角落它们不仅不会改善 IGD还会掩盖真实前沿的不均匀。第二个细节是如果理论前沿 PFtrue 采样点太少mean(min(...)) 会低估误差DTLZ1 的理论前沿需要至少 10000 个采样点ZDT1 的前沿是连续曲线常用 500 到 1000 个点已够。5.2 用 IGD 曲线判断参考点分布是否合理在 NSGAIII_main.m 的循环里记录每一代 IGD最后用 semilogy 画出来比看最终数字有用得多。history zeros(1, T); for t 1:T % ... 一次进化 ... [rank, ~] NDSort(PopObj); Front PopObj(rank 1, :); history(t) IGD(Front, PFtrue); end semilogy(1:T, history); xlabel(Generation); ylabel(IGD);曲线如果快速下降后进入平台说明收敛和多样性已经稳定继续增加代数收益不大如果曲线长期不下调小 disM 或增大变异概率比增加代数更有效如果曲线先降后升说明种群多样性崩溃常见原因是参考点数量太少或归一化失败。一个具体技巧是同时运行 M2 和 M3 的两份配置用同一把随机种子观察三目标的平台是否明显高于双目标。若是优先怀疑参考点分布而不是遗传算法参数。5.3 一个验证 UniformPoint 的三行检查最后回到参考点本身。无论替换成什么模型先单独验证 UniformPoint 的输出[W, N] UniformPoint(100, 3); assert(all(abs(sum(W, 2) - 1) 1e-6)); scatter3(W(:, 1), W(:, 2), W(:, 3), filled);看到点均匀铺在一个三角形内说明参考点生成正常再回到 NSGAIII_main.m 调种群规模。若点出现聚集或者有一侧明显稀疏检查 nchoosek 输入参数是否在你的 MATLAB 版本中被当作 double 导致精度丢失必要时改用 nchoosek(int32(H1M-1), int32(M-1)) 重新生成。本文还有配套的精品资源点击获取