简介：这份MATLAB版NSGA-III算法实现面向多目标进化算法学习者与研究者，聚焦两目标优化问题的代码理解与实验复现，也可作为课程设计、科研入门和算法对比的参考。压缩包共16个文件，以14个.m源码为主，另含1个原始链接文件和1个许可说明文本，整体约11KB，体积轻量，便于在MATLAB中直接查看与运行。已有5057人学习，说明其在NSGA-III入门实践中具有一定关注度。代码保留主流程、参考点生成、非支配排序、种群归一化、参考点关联、交叉变异及结果绘图等关键环节，并配有较完整的中文注释；作者对尚未完全弄懂之处留有空白或问号，适合读者边调试边讨论。对希望从源码层面理解NSGA-III选择机制、约束处理思路与两目标实验配置的读者，可据此梳理算法结构、对照原版实现并记录排错与改进方向。

1. 2 目标也用 NSGA-III？先看清这份 MATLAB 源码的边界

手里有一份 NSGA-III 算法的 MATLAB 版本，注释写满每一行，第一次打开时我第一反应是：2 目标问题用它是不是有点用力过猛。多目标进化算法里，2 目标常年是 NSGA-II 的地盘，NSGA-III 的名声更多挂在 3 目标及以上。但把这份源码跑完 ZDT1 之后，结论要修正：它不是不能打 2 目标，而是你得知道参考点机制在 2 目标下退化成了什么样子，否则种群规模、参考点层数两个参数一错，前沿分布会明显发飘。这份资源适合两类人：一是课程设计或论文复现里必须交出一版可读、可改的 NSGA-III 实现；二是手里有 3 目标以上问题、需要一套带中文注释的 MATLAB 骨架做二次开发。注释密度高，意味着你改目标函数、换约束处理时不用先花两天读黑匣子。

2. 参考点、归一化与小生境：NSGA-III 的三根支柱怎么在代码里落地

NSGA-III 和 NSGA-II 的差别不在交叉变异，而在环境选择：NSGA-II 用拥挤距离维持分布，NSGA-III 用一组预置参考点把种群往参考方向拽。这份 MATLAB 版本把这条链路拆成了独立函数，读的时候别一上来就啃主循环，先按“非支配排序 → 参考点生成 → 归一化与关联 → 小生境保留”的顺序看，每一段都能单独跑。下面三节按这个顺序拆，代码可以直接抄进自己的工程。

2.1 非支配排序：先分层，再谈分布

环境选择的第一步永远是把种群按 Pareto 支配关系分层。NSGA-III 的临界层处理比 NSGA-II 更依赖分层结果，因为只有临界层才进入小生境筛选，前面整层直接晋级。这份源码里的非支配排序用的是朴素双层循环，注释写得很直白，适合对照论文读，但规模上到几千个体就会明显拖慢。课程设计和小规模验证够用，真要做上万次评估的实验，建议换成快速非支配排序。

function [FrontNo, MaxFNo] = NonDominatedSorting(PopObj) % 输入 PopObj: N x M 目标值矩阵，每行一个解，已统一为最小化 % 输出 FrontNo: 1 x N 层级编号，1 为第一前沿 % 输出 MaxFNo: 最大层号 [N, M] = size(PopObj); FrontNo = inf(1, N); % 未分层前全部置 inf MaxFNo = 0; remaining = true(1, N); % 标记还没被分层的个体 while any(remaining) MaxFNo = MaxFNo + 1; current = find(remaining); dominated = false(1, numel(current)); for i = 1:numel(current) for j = 1:numel(current) if i == j continue; end % 支配判断：所有目标不差，且至少一个目标严格更好 if all(PopObj(current(i),:) <= PopObj(current(j),:)) && ... any(PopObj(current(i),:) < PopObj(current(j),:)) dominated(j) = true; end end end % 本轮没被任何个体支配的，就是当前前沿 FrontNo(current(~dominated)) = MaxFNo; remaining(current(~dominated)) = false; end end

代码的逻辑是每轮剥一层：dominated标记被当前层支配的个体，一轮结束后没被标记的就是当前前沿，从remaining里移除再进下一轮。FrontNo最终记录每个个体属于第几前沿，MaxFNo是总层数，环境选择里用它判断临界层。参数上唯一要留意的是PopObj必须已经转成最小化，如果目标函数是最大化，进这个函数之前要取负号，否则分层方向会反。这个实现的时间复杂度是每层 O(N²M)，层数多的时候更慢；我一般会在函数入口加一句if N > 2000, warning('建议换快速非支配排序'); end，算是给自己一个提醒。

2.2 Das-Dennis 参考点：p 一改，种群规模跟着变

参考点是 NSGA-III 的方向盘。Das-Dennis 方法把单位超平面按每维 p 段划分，生成一组均匀分布的点，每个点代表一个参考方向。M 个目标、p 段划分对应的参考点数是组合数C(M+p-1, p)，这个数直接决定种群规模该取多大。2 目标和 3 目标的差别在这里非常明显：2 目标 p=99 才 100 个参考点，3 目标 p=12 就有 91 个。很多人从 3 目标代码直接改成 2 目标，p 还留着 12，结果只有 13 个参考方向，种群规模却设成 100，环境选择阶段就会反复补随机解，前沿看起来又稀又散。

function Z = GenerateReferencePoints(M, p) % 生成 Das-Dennis 结构化参考点 % M: 目标个数；p: 每维划分段数 % Z: 参考点矩阵，每行一个点，满足 sum(Z,2) == 1 Z = []; Z = DasDennisRecursive(M, p, [], Z); Z = Z / p; % 从整数格点归一化到单位超平面 end function Z = DasDennisRecursive(M, p, prefix, Z) % 递归填充前 M-1 维，最后一维由 p - sum(prefix) 补齐 if M == 1 Z = [Z; prefix, p - sum(prefix)]; else for i = 0:(p - sum(prefix)) Z = DasDennisRecursive(M-1, p, [prefix, i], Z); end end end

递归的终止条件在M == 1，最后一维不再循环，直接用p - sum(prefix)补齐，保证每个参考点的坐标和为 p，除完 p 之后落在单位超平面上。prefix是行向量，拼接时注意方向，列向量拼进去会变成矩阵，后面关联距离计算会报维度错。参数p越大参考点越密，但点数按组合数增长，4 目标 p=10 就已经 286 个点，5 目标再往上要考虑计算量。我一般先按目标数查表定 p，再让种群规模N等于参考点数，这样环境选择里每个小生境至少能分到一个候选解，不会出现大面积空桶。

这张表不是让你死记，而是提醒一件事：改M或p之后，N要重新算，不能沿用上一组参数。源码里如果写了固定N = 100，切到 3 目标 p=12 时就会多出 9 个个体，临界层的小生境循环会多跑几轮，虽然不报错，但分布会偏。

2.3 归一化、关联与小生境：环境选择里最容易写错的三步

临界层进入小生境之前要先归一化，目的是把不同量纲的目标拉到同一尺度，否则参考线关联会被某个大数值目标主导。归一化的标准流程是：先减理想点平移，再求每个目标方向的极值点，用极值点构造超平面求截距，最后除以截距。这份源码把这一步单独写成NormalizePopulation，注释里标了“超平面构造可能奇异”，这是实话，极端点重复或共线时Extreme \ b会给出 Inf 或负截距，后面整层关联全乱。

function [Zn, intercepts] = NormalizePopulation(PopObj, ideal) % PopObj: N x M 目标值；ideal: 1 x M 理想点 % Zn: 归一化后的目标值；intercepts: 各维截距 F = PopObj - repmat(ideal, size(PopObj,1), 1); M = size(F, 2); % 构造 M 个 ASF 权重，略微偏离单位阵避免除零 W = eye(M) + 1e-6; ASF = zeros(size(F,1), M); for i = 1:M ASF(:,i) = max(F ./ repmat(W(i,:), size(F,1), 1), [], 2); end [~, idx] = min(ASF, [], 1); % 每个目标方向取 ASF 最小的个体 Extreme = F(idx, :); b = ones(M, 1); % 超平面可能奇异，先看条件数再决定用不用伪逆 if rcond(Extreme) < 1e-10 a = pinv(Extreme) * b; else a = Extreme \ b; end intercepts = 1 ./ a; % 截距出现非正数或 Inf 时，用该维最大值兜底 bad = ~isfinite(intercepts) | intercepts <= 0; if any(bad) intercepts(bad) = max(F(:,bad), [], 1) + 1e-6; end Zn = F ./ repmat(intercepts, size(F,1), 1); end

W加1e-6是为了避免 ASF 计算时除以零，代价是极值点会轻微偏移，但对关联结果影响很小。rcond判断条件数比直接try-catch更可控，条件数太小就用pinv，虽然慢一点但不会把整层带崩。截距兜底那三行是血泪经验：有一次跑 3 目标约束问题，两个极值点几乎重合，截距出现负值，归一化后所有解都关联到同一条参考线，前沿直接缩成一个点，排查了半天才定位到超平面。归一化之后是关联，每个解关联到距离最近的参考线，距离用点到参考线的垂直距离，不是欧氏距离。

function [pi, Dist] = AssociateToReferencePoint(Zn, Z) % Zn: 归一化目标值；Z: 参考点 % pi: 每个解关联的参考点索引；Dist: 对应的垂直距离 N = size(Zn, 1); K = size(Z, 1); Dist = zeros(N, K); for k = 1:K w = Z(k,:) / norm(Z(k,:)); % 参考线单位方向 proj = Zn * w'; % 投影长度 Dist(:,k) = sqrt(sum((Zn - proj * w).^2, 2)); end [~, pi] = min(Dist, [], 2); end

proj是解在参考线方向上的投影，Zn - proj * w是垂直分量，二范数就是垂直距离。这里要求Z的每一行不能全零，Das-Dennis 生成的参考点不会全零，但如果自己手工改参考点，记得检查。关联完进入小生境保留：统计每个参考点已经选了多少个解，优先选计数少的参考点，计数为 0 时选距离最近的，计数大于 0 时随机选，直到选够临界层该补的数量。

function Sel = Niching(Zn, Z, Last, K) % Zn: 临界层归一化目标值；Z: 参考点 % Last: 临界层在合并种群中的索引；K: 还需要选几个 % Sel: 1 x numel(Last) 逻辑向量 [pi, Dist] = AssociateToReferencePoint(Zn, Z); rho = histcounts(pi, 1:size(Z,1)+1); % 每个参考点已选计数 Sel = false(1, numel(Last)); while sum(Sel) < K minRho = min(rho); candidates = find(rho == minRho); k = candidates(randi(numel(candidates))); % 计数最小的参考点里随机挑一个 members = find(pi == k & ~Sel); if isempty(members) rho(k) = inf; % 该桶没有可用解，封掉避免死循环 else if rho(k) == 0 [~, idx] = min(Dist(members, k)); % 空桶选距离最近的 else idx = randi(numel(members)); % 非空桶随机选，保多样性 end Sel(members(idx)) = true; rho(k) = rho(k) + 1; end if all(isinf(rho)) Sel(find(~Sel, K - sum(Sel))) = true; % 兜底：参考点全封了就直接补 break; end end end

rho的初始值来自临界层之前的已选个体，不是从零开始，这样前面层已经占掉的参考方向会被降权。rho(k) = inf那行是防死循环的关键，没有它，某个参考点关联不到任何候选解时while会一直转。最后的兜底分支在参考点数少于K时触发，正常参数下不该走到，但 2 目标 p 设得太小、N又设得很大时就会遇到。

3. 从 ZDT1 跑通第一个 Pareto 前沿：参数表、主循环与出图

把源码解压后别急着改目标函数，先拿一个 2 目标基准问题把整条链路跑通，确认参考点数、种群规模、分层结果、出图都正常，再换自己的问题。这一章按入口脚本、参数配置、主循环、结果输出四段走，每一步都给可执行的代码和参数说明。

3.1 目录与入口脚本：先固定随机种子再谈复现

源码包的目录一般按“入口 + 核心 + 问题 + 工具”拆，常见结构是入口脚本放根目录，核心算子放NSGA3/，测试问题放Problems/，绘图和指标放Utils/。这种拆法不是为了好看，而是让你换目标函数时只动Problems/一个文件，核心选择逻辑不用碰。第一次运行前先确认 MATLAB 当前文件夹在根目录，否则addpath没加对，函数找不到会报Undefined function。

• run_nsga3.m：入口脚本，参数集中在这里
• NSGA3/EnvironmentalSelection.m：环境选择，含非支配排序、归一化、小生境
• NSGA3/GenerateReferencePoints.m：Das-Dennis 参考点
• NSGA3/GeneticOperator.m：SBX 交叉与多项式变异
• Problems/ZDT1.m：2 目标测试函数
• Utils/PlotFront.m：Pareto 前沿绘图

入口脚本第一行建议固定随机种子。进化算法本身是随机的，不固定种子，两次运行的 HV 能差出几个点，排查问题时根本分不清是参数改了还是运气变了。

%% run_nsga3.m clear; clc; rng(2024); % 固定随机种子，保证复现 addpath(genpath(pwd)); % 把子目录一次性加入路径 M = 2; % 目标个数 p = 99; % 2 目标取 99，生成 100 个参考点 Z = GenerateReferencePoints(M, p); N = size(Z, 1); % 种群规模与参考点数对齐 D = 10; % 决策变量维数 maxGen = 250; % 迭代代数 lb = zeros(1, D); % 下界 ub = ones(1, D); % 上界 Population = InitializePopulation(N, D, lb, ub); for gen = 1:maxGen MatingPool = TournamentSelection(Population, N); Offspring = GeneticOperator(MatingPool, lb, ub); Offspring = EvaluatePopulation(Offspring, @ZDT1); [Population, FrontNo] = EnvironmentalSelection([Population, Offspring], N, Z); for i = 1:numel(Population) Population(i).rank = FrontNo(i); % 回写层级，出图要用 end if mod(gen, 25) == 0 fprintf('第 %d 代，第一前沿个体数 %d\n', gen, sum(FrontNo == 1)); end end

rng(2024)的种子值不重要，重要的是固定下来之后每次跑同一组参数结果一致，调参时才有可比性。addpath(genpath(pwd))会把当前目录及所有子目录加进搜索路径，适合这种多文件夹结构；如果你的 MATLAB 当前目录不是源码根目录，换成绝对路径。N = size(Z,1)是这份代码能跑通的关键，2 目标 p=99 时 N=100，3 目标 p=12 时 N=91，换目标数不用手工改种群规模。

3.2 参数表：M、p、N、D、maxGen 怎么配

参数没有万能值，但有一组“不翻车”的起点。下面这张表是 2 目标和 3 目标各跑过几轮之后留下的常用组合，先按这个配，跑通再动。

2 目标用 p=99 看起来夸张，但 2 目标的参考点数就是 p+1，p=99 才 100 个点，计算量并不大。真正要小心的是 3 目标以上：p=12 出 91 个点，p=15 就出 136 个点，p=20 出 231 个点，每代环境选择都要做关联距离矩阵，点数翻倍时间大概翻三倍。我一般先按表跑，看第一前沿数量和 HV 再决定要不要加密。

3.3 主循环与遗传算子：SBX、多项式变异、合并选择

主循环的骨架是“选择 → 交叉变异 → 评估 → 合并环境选择”，这份源码把交叉和变异包在GeneticOperator里，内部用 SBX 和多项式变异。SBX 的核心是模拟二进制交叉的分布指数，eta_c越大子代越靠近父代，探索能力越弱；多项式变异同理。下面两个函数是源码里算子的简化版，注释保留了中文说明，可以直接替换进去。

function [c1, c2] = SBX(p1, p2, lb, ub, pc, eta) % 模拟二进制交叉 % pc: 交叉概率；eta: 分布指数，常用 15~20 D = numel(p1); c1 = p1; c2 = p2; if rand > pc return; end for j = 1:D if rand <= 0.5 && abs(p1(j) - p2(j)) > 1e-12 y1 = min(p1(j), p2(j)); y2 = max(p1(j), p2(j)); u = rand; beta = 1 + 2 * (y1 - lb(j)) / (y2 - y1); alpha = 2 - beta^(-(eta + 1)); if u <= 1 / alpha betaq = (u * alpha)^(1 / (eta + 1)); else betaq = (1 / (2 - u * alpha))^(1 / (eta + 1)); end c1(j) = 0.5 * ((y1 + y2) - betaq * (y2 - y1)); c2(j) = 0.5 * ((y1 + y2) + betaq * (y2 - y1)); c1(j) = min(max(c1(j), lb(j)), ub(j)); % 越界拉回 c2(j) = min(max(c2(j), lb(j)), ub(j)); end end end function y = PolynomialMutation(x, lb, ub, pm, eta) % 多项式变异 % pm: 每个变量的变异概率，常取 1/D；eta: 分布指数，常用 20 y = x; for j = 1:numel(x) if rand <= pm u = rand; if u < 0.5 delta = (2 * u)^(1 / (eta + 1)) - 1; else delta = 1 - (2 * (1 - u))^(1 / (eta + 1)); end y(j) = x(j) + delta * (ub(j) - lb(j)); y(j) = min(max(y(j), lb(j)), ub(j)); end end end

SBX里beta的分母是y2 - y1，两个父代在该维相等时跳过，避免除零；alpha和betaq是标准 SBX 公式，照着论文抄不会错。PolynomialMutation的pm常取1/D，意思是每个个体平均有一个变量发生变异，D 越大单变量变异概率越低，整体变异强度不变。这两个算子都不是 NSGA-III 特有，换 NSGA-II 也能用，所以改的时候不用动它们，把精力放在环境选择上。

3.4 出图与指标：Pareto 前沿、HV 粗算

跑完之后第一件事是把第一前沿画出来。2 目标问题直接散点图，横轴 f1 纵轴 f2，看前沿是否连续、是否被参考方向拽成均匀分布。源码里Population是结构体数组，目标值存在objs字段，层级存在rank字段，出图前先按rank == 1筛一遍。

Front = Population([Population.rank] == 1); figure; plot(Front.objs(:,1), Front.objs(:,2), 'o', 'MarkerSize', 4); xlabel('f_1'); ylabel('f_2'); title(sprintf('NSGA-III 第一前沿，个体数 %d', numel(Front))); grid on; saveas(gcf, 'pareto_front.png');

如果前沿出现明显的断点，先别怀疑算法，检查三件事：参考点数是否和 N 对齐、归一化截距是否有负值、临界层小生境是否走了兜底分支。这三件事在下一章展开。2 目标 ZDT1 的理论前沿是f2 = 1 - sqrt(f1)，可以把理论曲线叠上去对比，肉眼就能看出收敛和分布。3 目标以上散点图不够用，建议直接算 HV，第 5 章给了一个不依赖工具箱的 2 目标 HV 粗算函数。

4. 避坑记录：参考点数、矩阵奇异与 2 目标退化

这份源码注释多，但注释不会告诉你参数之间的耦合。下面五条是实际跑的时候最容易翻车的地方，每条按“现象 → 原因 → 解决”写，遇到问题先对照这里。

4.1 参考点数和种群规模不匹配，选择阶段选不满

现象：环境选择跑完，种群数量少于 N，或者临界层补进来的解明显偏多，前沿稀疏。原因：N没有跟着p和M变。2 目标 p=12 只有 13 个参考点，N还设 100，小生境循环里每个参考点最多关联几个解，补到后面rho全被置inf，兜底分支直接随机补，分布就散了。解决：入口脚本里用N = size(Z,1)，不要手写固定值；如果确实需要更大种群，取参考点数的整数倍，并在小生境里允许同一个参考点选多个解，同时把兜底分支的随机补改成按距离补。

4.2 超平面奇异，截距出现 Inf 或负数

现象：归一化之后Zn出现NaN或极大值，关联结果全部指向同一条参考线，第一前沿缩成一团。原因：ASF 选出的极值点重复或共线，Extreme \ b条件数极差，截距算出非正数。解决：NormalizePopulation里加rcond判断，条件数小于1e-10用pinv；截距算出后统一检查~isfinite或<= 0，用该维最大值加1e-6兜底。这个兜底不是万能药，如果频繁触发，说明种群多样性已经丢了，要回头检查交叉变异概率。

4.3 2 目标直接套 3 目标参数，前沿又稀又散

现象：从 3 目标代码切到 2 目标，p还是 12，N还是 100，跑出来前沿只有十几个点，分布还不如 NSGA-II。原因：2 目标参考点数是p+1，p=12 只有 13 个方向，参考方向退化成一条线段上的等分点，小生境分桶粒度太粗。解决：2 目标把 p 提到 49 或 99，让参考点数和 N 对齐；如果问题本身就是 2 目标、对分布要求高，直接换 NSGA-II 的拥挤距离更省事，NSGA-III 的优势在 3 目标以上才明显。

4.4 中文注释乱码，打开全是问号

现象：MATLAB 编辑器打开.m文件，中文注释显示为乱码或问号，代码本身能跑但没法读。原因：文件编码是 GBK，编辑器按 UTF-8 解码，或者反过来。解决：统一成 UTF-8，MATLAB 预设里把默认编码改成 UTF-8，保存时选 UTF-8 无 BOM；如果拿到的是 GBK 文件，用支持编码转换的编辑器另存为 UTF-8。改完编码后clear all再跑，避免旧缓存。

4.5 约束问题没接约束支配，种群全不可行

现象：带约束的问题跑完，第一前沿全是不可行解，目标值看起来很好但违反约束。原因：环境选择只按目标值非支配排序，没有把约束违反度纳入支配判断。解决：把支配函数换成约束支配：可行解优先于不可行解，两个可行解比目标，两个不可行解比违反度。源码里如果只有无约束版本，在NonDominatedSorting入口前把cv字段加进比较逻辑，改动不大但必须做。

5. 用 HV 和参考点自适应验收：把“能跑”变成“可信”

跑出前沿只是第一步，验收要看超体积 HV 是否随代数上升、是否稳定。2 目标 HV 不用工具箱也能算：先按 f1 升序排，维护当前 f2 的最低值，逐点累加矩形面积。下面这个函数只处理 2 目标，3 目标以上建议用专门的多目标指标库，不要自己硬写。

function hv = ComputeHV2D(F, ref) % F: N x 2 目标值，已转最小化；ref: 1 x 2 参考点，需被 F 支配 F = sortrows(F, 1); hv = 0; prev = ref(2); for i = 1:size(F,1) if F(i,2) < prev hv = hv + (ref(1) - F(i,1)) * (prev - F(i,2)); prev = F(i,2); end end end

ref取比前沿最差值略大的点，比如 ZDT1 取[1.1, 1.1]，每次迭代记录 HV，画一条随代数上升的曲线。如果 HV 中途掉下去，多半是归一化截距兜底触发或小生境选了劣解，回头查第 4 章。参考点自适应是另一个验收后的改造点：固定 p 在 3 目标以上容易在稀疏区域丢解，常见做法是先粗后细，前一半代数用较小 p 快速收敛，后一半在稀疏区域加密参考点，加密时注意同步更新N，否则环境选择数量对不上。

约束支配的改动也放在这一层做，不要动主循环：

function flag = ConstrainedDominates(x, y) % x.cv / y.cv 为约束违反度，0 表示可行 if x.cv == 0 && y.cv > 0 flag = true; elseif x.cv > 0 && y.cv == 0 flag = false; elseif x.cv > 0 && y.cv > 0 flag = x.cv < y.cv; else flag = all(x.obj <= y.obj) && any(x.obj < y.obj); end end

cv的计算放在目标函数里，等式约束取绝对值，不等式约束取max(0, g)，多个约束求和或取最大都行，关键是同一问题内保持一致。从那以后我每次跑 NSGA-III 都先打印size(Z,1)和N，再检查intercepts是否全为正，最后才看前沿图；这个习惯帮我省掉了至少一半的“玄学”排查。希望帮到你。

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