
简介第三代非支配排序遗传算法NSGA-III的MATLAB完整实现代码源自开源项目Yarpiz主要面向需要处理高维、多目标复杂优化问题的科研人员、工程技术人员及相关专业学生。该算法借助非支配排序、拥挤距离与参考点机制来平衡收敛性和种群多样性尤其适合存在多个相互冲突目标的工程决策场景。压缩包共19个文件包含14个m格式源码文件完整算法逻辑、3个txt说明文档使用说明与版权信息以及2个原始项目链接整体大小仅15KB代码极为凝练。目前该资源已有2399人学习下载。实现代码完整覆盖主程序迭代流程、非支配排序、拥挤距离计算、参考点生成、交叉与变异等核心环节每个函数均配有细致注释可清晰理解NSGA-III的分层选择、精英保留及适应度分配过程。基于现有框架只需替换目标函数接口或调整种群规模、迭代次数等参数便能快速迁移到自己的多目标优化课题中既适合作为算法教学演示也可作为进一步扩展开发其他进化算法的轻量基础模板。1. 多目标优化与遗传算法NSGA-III 能在高维目标空间里做什么当目标数量超过三个的时候经典遗传算法NSGA-II里的拥挤度距离就会告诉你什么叫歇菜——高维解空间中几乎所有个体都被判定为不拥挤算法直接退化成随机搜索。这是我在做某工程调度优化时踩过的大坑。多目标优化这件事两个目标可以画Pareto前沿三个目标还能看三维图到了四个五个目标需要的是让算法在目标空间里均匀撒点保持多样性。NSGA-III就是为了解决这个问题出现的它用一组预定义的参考点替代拥挤度距离做多样性维持让解集在高维目标空间里依然分布均匀。这份MATLAB工程里参考点生成、非支配排序、归一化、关联选择、遗传算子全部带注释能直接跑通并迁移到你的业务优化问题上。2. 项目架构与模块设计从主循环到函数拆分的完整链路2.1 NSGA-III 与 NSGA-II 的本质区别拥挤度退场参考点登场NSGA-III 不是全新的算法它沿用了遗传算法的主框架初始化种群、交叉变异产生子代、合并父子种群、非支配排序、环境选择。和 NSGA-II 真正的分水岭在最后一步环境选择。NSGA-II 用拥挤度距离来保证解集多样性做法是把同一个非支配层里的个体按目标函数值排序计算每个个体两侧的矩形周长。这个机制在 2 目标、3 目标时非常好用但目标数量一涨每个个体的邻近个体变少拥挤度值迅速趋同选择压力被稀释种群漂移严重解集没法均匀覆盖整个前沿。NSGA-III 的思路是彻底换掉多样性度量在运行前预先在目标空间生成一组均匀分布的参考点环境选择时把每个个体关联到离它最近的参考点上然后按每个参考点附近应该保留多少个个体来做名额分配。这样解集的分布性是由参考点网眼决定的而不是由解之间的相对距离决定的目标数量增加时只需要加密参考点网格算法行为保持稳定。这也是为什么 NSGA-III 被称为基于参考点的非支配遗传算法。2.2 MATLAB 目录与函数职责这份代码的文件组织非常清爽每个.m文件只干一件事属于典型的模块化单文件设计。核心结构如下NSGA3_Code/ ├── main.m # 主入口参数定义、种群初始化、迭代主循环 ├── generate_reference_points.m # 参考点生成Das-Dennis 组合数枚举 ├── evaluate_objective.m # 目标函数默认 DTLZ2 三目标测试问题 ├── fast_non_dominated_sort.m # 快速非支配排序返回层级与前沿集合 ├── adaptive_normalization.m # 自适应归一化极值点、超平面、截距 ├── associate.m # 关联参考点垂直距离计算与最小距离匹配 ├── environmental_selection.m # 环境选择小生境计数与精英保留 ├── sbx_crossover.m # 模拟二进制交叉 SBX ├── polynomial_mutation.m # 多项式变异 └── genetic_operators.m # 调用交叉与变异生成完整子代种群每个函数的输入输出都是标准矩阵种群统一用 N×n_var 矩阵表示目标值统一用 N×M 矩阵表示参考点统一用 R×M 矩阵表示。模块之间没有全局变量传递所有数据通过函数参数显式流动。这种设计的好处是你可以单独抽出任意一个函数做单元测试比如只跑generate_reference_points.m观察参考点分布再决定要不要调整划分密度。2.3 主程序 main.m参数初始化与迭代主循环主循环的代码结构决定了整个算法的运行节奏先把骨架帖出来后面每一章再深入拆里面的模块。%% 参数设置 M 3; % 目标函数个数 p 12; % 参考点每个维度的划分数 R nchoosek(M p - 1, p); % 参考点数量Das-Dennis 组合数 N R; % 种群大小一般取参考点数或整数倍 max_gen 250; % 最大迭代代数 n_var 30; % 决策变量个数 lb zeros(1, n_var); % 决策变量下界 ub ones(1, n_var); % 决策变量上界 %% 初始化 population lb rand(N, n_var) .* (ub - lb); ref_points generate_reference_points(M, p); %% 迭代主循环 for gen 1:max_gen % 生成子代SBX 交叉 多项式变异 offspring genetic_operators(population, lb, ub); % 合并父子两代种群 combined [population; offspring]; % 计算全部个体的目标函数值 obj_vals evaluate_objective(combined); % 快速非支配排序得到 rank 向量 [rank, ~] fast_non_dominated_sort(obj_vals); % 自适应归一化更新理想点求截距得到 f_norm z_min min(obj_vals, [], 1); [f_norm, ~, ~] adaptive_normalization(obj_vals, z_min); % 每个个体关联到最近的参考点 [match_idx, match_dist] associate(f_norm, ref_points); % 环境选择保留 N 个精英个体进入下一代 selected_idx environmental_selection(combined, rank, match_idx, match_dist, ref_points, N); population combined(selected_idx, :); % 打印当前最优前沿信息 if mod(gen, 50) 0 fprintf(gen%d \\n, gen); end end这段代码有几个参数值得细说。p是参考点划分份数直接决定参考点数量 R当 M3、p12 时 R91种群一般就设成 91 或者 182 这样的整数倍。lb、ub是决策变量的上下界必须与你要优化的实际问题匹配。主循环里的关键步骤顺序是固定的先合并种群再做非支配排序然后归一化最后关联选择。前两步搞清个体优不优后两步解决多样性够不够缺一不可。3. 核心代码拆解参考点生成与自适应归一化3.1 Das-Dennis 参考点生成组合数学的递归实现参考点是 NSGA-III 的基石它决定了算法期望解集在目标空间里的分布形态。Das-Dennis 方法生成的是单纯形上的均匀网格点要求每个参考点坐标分量之和等于 1且每个分量是 p 等分中的整数倍。参考点数量由组合数C(Mp-1, p)决定比如 M3、p12 时就是 C(14,12)91 个点。function ref generate_reference_points(M, p) % M : 目标函数个数一般 3~15 % p : 每个维度上的划分份数正整数 % ref: 大小为 C(Mp-1, p) × M 的参考点矩阵每行和为 1 if M 1 % 单目标退化参考点就是 1 ref 1; return; end container []; % 递归枚举所有满足 x1 x2 ... xM p 的非负整数组合 for i 0:p sub generate_reference_points(M - 1, p - i); for j 1:size(sub, 1) container [container; i, sub(j, :)]; end end % 除以 p让每一行分量和为 1 ref container / p; end这段递归的实现思路是固定第一个维度的取值 i然后把剩余 p-i 的划分交给低一维的同样问题递归处理。当 M2 时递归到基准条件返回p即[p]除以 p 后变成 1回溯拼接后得到所有整数分区。最终除以 p 后每一行的分量都在 0 到 1 之间且总和为 1。需要说明的是参考点本身不是解只是目标空间里的虚拟锚点它不参与进化只是在环境选择时用来标定哪个区域该保人。这个函数在参数选择上的一个常见疑问是 p 取多少。p 越大参考点越密解集分布越细但计算量也会按组合数增长。工程上我一般这样定M3 时 p12 或 13 得到 91 或 105 个参考点M5 时 p6 得到 210 个点M8 时 p3 得到 120 个点。这个经验值能平衡覆盖密度和计算代价。3.2 理想点与极值点归一化前的两个关键输入归一化是整个 NSGA-III 中最容易翻车的一步。先说为什么要归一化多目标问题里不同目标的量纲往往不一致比如成本是元、时间是小时、风险是百分比直接用原始值算距离量级大的目标会完全压过量级小的目标。归一化的目的是把每个目标都变换到接近 0 到 1 的范围内让关联参考点的计算不偏向任何一维。归一化的第一步是平移用每个目标的历史最小值构成理想点 z_min从所有个体的目标值里减掉它让理想点落到原点。第二步是找极值点即每个目标轴方向上最突出的个体。这一步用的是加权切比雪夫函数 ASF它的核心思想是给目标 i 一个很大的权重 1其它目标给一个极小的权重比如 1e-6然后对每个个体计算 max(f_j / w_j)使这个值最小的个体就是沿目标 i 方向最极端的解。function [f_norm, z_min, intercepts] adaptive_normalization(f, z_min) % f : 种群目标值矩阵 N × M % z_min : 理想点各维目标的历史最小值 % f_norm : 归一化后的目标值 % intercepts: 超平面在每个目标轴上的截距 N size(f, 1); M size(f, 2); % 1. 平移让理想点落到原点 f_shifted f - repmat(z_min, N, 1); % 2. 用加权切比雪夫函数找每个目标轴的极值点 extreme zeros(M, M); for i 1:M w 1e-6 * ones(1, M); w(i) 1; % 计算 ASF 值取最小对应的个体作为极值点 asf max(f_shifted ./ repmat(w, N, 1), [], 2); [~, idx] min(asf); extreme(i, :) f_shifted(idx, :); end % 3. 构造超平面 extreme * a 1解出截距 if rcond(extreme) 1e-10 % 矩阵接近奇异改用伪逆避免数值爆炸 a pinv(extreme) * ones(M, 1); else a extreme \ ones(M, 1); end intercepts 1 ./ a; % 4. 基于截距做归一化 f_norm f_shifted ./ repmat(intercepts, N, 1); end代码里f_shifted ./ repmat(w, N, 1)中的w是行向量MATLAB 隐式扩展会用每个个体的目标值逐项除以权重向量。因为 w(i)1 而其它维度权重极小ASF 值主要由目标 i 分量主导所以min(asf)找出来的个体就是该目标轴方向上最极端的那个解。rcond是条件数评估函数值太小表示矩阵接近奇异此时普通除法会得到不可靠结果改用pinv是稳妥做法。实际应用中极值点矩阵奇异很常见比如两个目标的极值点是同一个个体几乎每次跑都会遇到所以这个判断不是防御性编程是必需品。3.3 超平面截距求解从线性方程到数值稳定性第三步的数学含义值得单独说清楚。拿到 M 个极值点后要构造一个穿过这些极值点的超平面然后求出这个超平面与每个目标轴的交点也就是截距。超平面方程可以写成extreme * a 1其中extreme是 M×M 矩阵每行是一个极值点的坐标a是待求系数向量截距就是1./a。实际运行中这个方程组经常会病态。常见情况是某些极值点坐标非常接近导致矩阵行向量近似线性相关rcond很小解出来的 a 有的分量接近零甚至为负。一旦截距为负归一化后该维目标值全部变负关联参考点算出来的距离全是错的迭代下去整个种群就废了。所以我在代码里做了两道防线第一道是用rcond判断是否改用伪逆第二道是后续在环境选择时对截距做合理性检查。经验来看如果目标函数本身存在强相关性或者种群早期极值点还没站稳这一步最容易出 NaN。4. 核心代码拆解非支配排序与精英选择机制4.1 快速非支配排序两两比较的支配链非支配排序和 NSGA-II 完全一致逻辑是如果一个个体 p 的所有目标值都不劣于另一个个体 q且至少一个目标严格优于 q那么 p 支配 q。把所有不被任何人支配的个体划为第一层去掉它们之后再找出新的不被支配个体划为第二层依此类推。层级越小个体越优秀。function [rank, fronts] fast_non_dominated_sort(obj_vals) % obj_vals : 种群目标值矩阵 N×M % rank : 每个个体的非支配层级1 为最优层 % fronts : 元胞数组每个元素是该层个体索引向量 N size(obj_vals, 1); dominated_set cell(N, 1); % p 所支配的个体集合 n_dominated zeros(N, 1); % 支配 p 的个体数量 for p 1:N for q 1:N if p q, continue; end if all(obj_vals(p, :) obj_vals(q, :)) any(obj_vals(p, :) obj_vals(q, :)) dominated_set{p}(end 1) q; elseif all(obj_vals(q, :) obj_vals(p, :)) any(obj_vals(q, :) obj_vals(p, :)) n_dominated(p) n_dominated(p) 1; end end end % 第一层不被任何个体支配的个体 fronts {}; current_front find(n_dominated 0); layer 1; % 循环去掉当前层后下一层的 n_dominated 减到 0 即为新层 while ~isempty(current_front) fronts{layer} current_front; rank(current_front) layer; next_front []; for i current_front for j dominated_set{i} n_dominated(j) n_dominated(j) - 1; if n_dominated(j) 0 next_front(end 1) j; end end end layer layer 1; current_front next_front; end end这段代码的时间复杂度是 O(M·N²)M 是目标数N 是种群大小。两层循环做两两比较之后每个个体的被支配计数逐层递减。当种群规模是 200、目标数是 3 时单次排序要做 4 万次比较在 MATLAB 里完全可接受。但如果目标数上到 10 且种群上千这个排序会成为瓶颈优化方向是把内层比较向量化用矩阵运算替代 for 循环。需要留意的是all(obj_vals(p,:) obj_vals(q,:)) any(...)的组合条件不能写错。只写all不判断any会把两个完全相同的个体也判成互相支配导致排序结果错乱。这是个非常容易忽略的细节我在写第一版时就踩过种群里重复个体一多所有层都乱套了。4.2 关联参考点垂直距离计算归一化完成后每个个体在目标空间里是一个点每个参考点是从原点出发的一条射线。关联操作要回答的问题是这个个体最靠近哪条射线。距离度量用的是点到射线的垂直距离不是点到点的欧氏距离。差别在于参考点是射线方向不是固定坐标点如果算欧氏距离所有个体都会涌向靠近原点的参考点分布性会彻底失效。function [match_idx, match_dist] associate(f_norm, ref) % f_norm : 归一化目标值矩阵 N×M % ref : 参考点矩阵 R×M % match_idx : 每个个体关联的参考点编号1~R % match_dist : 个体到参考点射线的垂直距离 N size(f_norm, 1); R size(ref, 1); match_idx zeros(N, 1); match_dist zeros(N, 1); for i 1:N f_i f_norm(i, :); best_d inf; best_r 0; for j 1:R w ref(j, :); w2 dot(w, w); if w2 1e-12 % 参考点为原点时垂直距离就是个体到原点的模长 d norm(f_i); else % 垂直距离公式: |f - dot(f,w)/dot(w,w) * w| proj dot(f_i, w) / w2 * w; d norm(f_i - proj); end if d best_d best_d d; best_r j; end end match_idx(i) best_r; match_dist(i) best_d; end endproj dot(f_i, w) / w2 * w是向量投影的标准公式把个体坐标向量投影到参考点方向向量上再用原始向量减去投影得到垂直分量它的模长就是垂直距离。当参考点落在原点即参考点向量模长为零时投影公式分母为零需要兜底处理。实际运行中这种极端情况很少发生但如果划分份数 p0 或者目标数为 1 时就会触发。我习惯在计算前先用norm(w)判断一下防止除零警告刷屏。这个嵌套循环在种群大、参考点多的时候是性能热点。优化手段通常是先用矩阵批量算出所有个体到所有参考点方向的投影参数再一次性求距离矩阵可以把耗时降一到两个数量级。代码里保留双循环是为了可读性性能敏感时值得重构。4.3 小生境计数与精英选择最后一道筛选环境选择的目标是从合并种群大小为 2N中留下 N 个精英。策略是逐层纳入非支配层如果加入某一层后总数正好等于 N直接结束如果会超过 N就需要对这一层的个体做小生境筛选。小生境的核心思想是先统计已选个体在各个参考点附近的密度然后优先从密度最小的参考点附近挑人。function selected environmental_selection(combined, rank, match_idx, match_dist, ref, N) % combined : 合并后的种群个体或索引 % rank : 每个个体的非支配层级 % match_idx : 每个个体关联的参考点编号 % match_dist : 每个个体到参考点的垂直距离 % ref : 参考点矩阵 % N : 需要保留的种群大小 selected []; layer 1; % 逐层纳入个体直到下一层会超过 N while length(selected) length(find(rank layer)) N selected [selected, find(rank layer)]; layer layer 1; end % 最后一层需要从中挑选 K 个 K N - length(selected); last_layer find(rank layer); % 统计已选个体在各参考点上的小生境计数 R size(ref, 1); rho zeros(1, R); for i 1:length(selected) rho(match_idx(selected(i))) rho(match_idx(selected(i))) 1; end chosen []; while length(chosen) K % 找小生境计数最小的参考点取第一个 [min_rho, j] min(rho); % 该参考点关联的、还在最后一层中的候选个体 linked last_layer(match_idx(last_layer) j); if isempty(linked) % 没有候选个体把该参考点的计数设为无穷大跳过 rho(j) inf; continue; end % 从候选个体中选垂直距离最小的 [~, idx] min(match_dist(linked)); chosen(end 1) linked(idx); % 移除已选个体更新小生境计数 last_layer(last_layer linked(idx)) []; rho(j) rho(j) 1; end selected [selected, chosen]; end这个选择逻辑体现了 NSGA-III 的核心策略它不按个体的绝对适应度来抢名额而是优先保证每个参考点附近都有人。rho记录的是已选个体在各参考点的分布密度min(rho)选出的参考点意味着该区域的人还很少那么从这一层里挑一个离它最近的个体补进去就有助于让解集均匀铺满前沿。实现细节上有个容易出问题的点当某参考点没有候选个体时如果直接continue会陷入死循环因为 rho 最小值一直不变。把该参考点的 rho 设为inf是常用的处理技巧相当于把这块区域从排队中移除。此外last_layer(last_layer linked(idx)) []删掉的是匹配索引的第一个匹配项如果候选集合里正好有两个相同索引的个体这种方式是安全的因为个体索引本身是唯一的。5. 跑通与避坑参数调优和常见问题排查5.1 换掉默认测试函数把DTLZ2改成你自己的问题默认的evaluate_objective.m实现的是 DTLZ2 三目标测试函数这段代码的价值在于它示范了多目标测试函数的写法。换成实际问题时只需要改这一个文件。首先定义决策变量维度和目标个数然后写清从决策变量到目标值的计算过程。如果目标是成本和时间这种最小化指标直接返回函数值即可如果是最大化收益取负值或在外部做转换。以 DTLZ2 为例三目标版本的实现是每个个体都要算一个距离函数 g然后根据角度变量投影到三个目标轴上function f evaluate_objective(x) % x : 决策变量矩阵每一行是一个个体 % f : 目标值矩阵每一行对应一个个体的 M 个目标 M 3; % 目标个数 n size(x, 1); f zeros(n, M); for i 1:n % 距离函数 g对下标 M 的决策变量求平方和 g sum((x(i, M:end) - 0.5).^2); % 角度投影 f(i, 1) (1 g) * cos(0.5 * pi * x(i, 1)) * cos(0.5 * pi * x(i, 2)); f(i, 2) (1 g) * cos(0.5 * pi * x(i, 1)) * sin(0.5 * pi * x(i, 2)); f(i, 3) (1 g) * sin(0.5 * pi * x(i, 1)); end end替换成实际问题时g的部分换成你的约束处理和损失计算x(i, M:end)这种取子集的操作换成你的变量映射。需要特别注意的是决策变量的范围如果实际问题里变量下限是 -5、上限是 5必须同时改main.m里的lb和ub否则种群初始化和遗传操作的边界裁剪会不一致子代很容易跑到不可行区域。5.2 参数匹配目标数、参考点数、种群规模的关系NSGA-III 的参数不是各调各的它们之间有强约束关系。参考点数量由目标数 M 和划分数 p 通过组合数决定种群大小一般直接取参考点数或它的整数倍。如果种群大小和参考点数差距太大会出现两个问题种群比参考点少很多时有些参考点附近永远没人解集覆盖不全种群比参考点多数倍时计算量浪费在过多的个体上收敛速度反而变慢。目标数 M划分数 p参考点数量 R推荐种群大小 N3129191 或 182330496496410286286 或 5725621021083120120 或 240max_gen的设定主要看问题复杂度。测试函数通常 250 到 500 代收敛工程问题建议先跑 500 代观察 Pareto 前沿是否还在移动如果最后几十代的解集变化很小说明可以提前停止如果变化明显继续加代数。迭代代数的判断没有玄学可言多跑几次对比 HV 指标或前沿分布就能找到合适的值。5.3 避坑记录5个真实踩过的坑第一个坑参考点数量与种群大小不匹配时环境选择死循环。现象是程序跑到选择阶段卡住不动CPU 占用 100%几十分钟不退出。原因是某参考点始终没有候选个体rho(j)一直是最小值代码不停重试同一个空集合。解决办法是给空的参考点设inf同时限制最大重试次数两种手段我都加上才彻底稳定。第二个坑归一化截距解出 NaN。现象是种群进化到某一代后目标值全是 NaN后续所有距离计算失效。原因是极值点矩阵接近奇异用普通除法得到了不可靠结果。我的解决方法是先判断rcond(extreme)小于阈值改用伪逆再对截距做一次合法性检查出现负值或 NaN 就回退为 1相当于不归一化该维度。第三个坑把垂直距离写成欧氏距离。现象是解集集中在少数几个参考点附近其它区域完全空白。网上不少简化实现确实把associate写成了sqrt(sum((f_i - w).^2))这算的是点到点的距离不是点到射线的距离。正确的做法是用投影公式让我吃了一次亏后现在每次检查别人的实现第一件事就是看associate函数里有没有除法dot(f_i,w)/dot(w,w)。第四个坑高维目标下参考点数量爆炸。现象是 M8、p10 时参考点数量 C(17,10)19448内存直接崩溃。这不是算法 bug是组合数增长的本质代价。解决思路是降低 p 或用两层参考点生成策略边界层用大 p、内部层用小 p总的参考点数能压到原来的十分之一。第五个坑目标量级悬殊时归一化失效。现象是两个目标一个范围在 0 到 1另一个在 0 到 1e6截距完全被大目标主导归一化后小目标的区分度接近零。解决方法是归一化前先对每个目标做独立缩放把量级拉到同一个水平线再进入自适应归一化流程。6. 二次开发把 NSGA-III 改造成自己的优化引擎6.1 从测试函数到实际工程问题约束处理与变量替换拿到代码后最典型的改造路径是把 DTLZ2 换成一个真实的调度问题或参数优化问题。第一步改evaluate_objective.m定义好输入输出的矩阵维度第二步改main.m里的lb、ub和n_var。如果你的问题带约束条件常见做法是在目标函数里加惩罚项罚函数系数从 1e3 开始调不要让惩罚项彻底压过真实目标否则早期种群全是不可行解算法探索效率极低。6.2 性能验证用超体积指标检查改造效果改造后不要只看前沿图超体积指标 HV 是更可靠的量化验证手段。HV 计算的是解集在目标空间中覆盖的体积值越大说明解集的收敛性和分布性越好。在 MATLAB 里可以用pareto函数画二维前沿三维用plot3四维以上就要靠 HV 了。每次改完代码我习惯固定相同随机种子跑 NSGA-II 和 NSGA-III 各 10 次对比平均 HV这样能量化出参考点带来的多样性收益。6.3 进阶改造方向并行化与热启动当评估函数计算成本很高时主循环里的evaluate_objective可以用parfor并行替换for在并行池开启的情况下性能提升接近核数倍率。热启动是另一个实用的改造技巧把上一轮进化得到的最终种群作为下一轮参数调整后的初始种群避免重复冷启动的收敛开销。从那以后我每次接手新的优化项目都会强制走一遍这套流程先跑通测试函数验证代码正确性再替换评估函数然后用 HV 量化基线最后才动手并行化和调参。希望帮到你。本文还有配套的精品资源点击获取