尧图网站开发 
简介面向多目标优化研究者与 MATLAB 学习者的 NSGA-II 算法实现压缩包针对多个冲突目标同时优化问题提供一套完整可运行的遗传算法框架。压缩包共 13 个文件以 12 个 m 文件为主涵盖非支配排序、拥挤距离计算、锦标赛选择、交叉、变异及帕累托前沿绘图等核心模块另含 1 个来源链接文件整体大小仅 6KB结构紧凑、便于阅读。当前已有 989 人学习适合入门与进阶。通过这套 MATLAB 源码读者可以直观观察 NSGA-II 的迭代收敛过程理解快速非支配排序如何划分解的质量层次、拥挤距离如何维持种群多样性并可根据自身优化目标修改适应度函数、调整种群参数快速迁移到工程设计、投资决策等实际场景中。模块化的函数设计也便于逐块调试与教学演示作为多目标进化算法学习与二次开发的轻量工具。1. NSGA-II 多目标优化在 MATLAB 里跑通收敛性的关键从来不是把代码复制下来NSGA-IINon-dominated Sorting Genetic Algorithm II是过去二十年多目标优化领域被引用最多的算法骨架。它用非支配排序维持解的层级用拥挤度距离保证种群在 Pareto 前沿上铺得开再用精英保留策略让每一代的最优个体不会丢失。MATLAB 里实现 NSGA-II 的代码在网上有成百上千份但多数人下载后第一轮跑出来的结果不是前沿残缺就是收敛曲线震荡甚至种群直接塌缩到某个单目标最优附近——这恰恰说明问题不在算法本身而在运行参数和收敛性判断这两件事上没被认真对待。这篇文章围绕「多目标运行」和「收敛性」两个核心词展开先拆 NSGA-II 在 MATLAB 中的完整实现链路再给出可直接运行的代码和参数表最后落到用 Hypervolume、IGD 等指标量化验证收敛性的具体脚本。读者如果是刚接触多目标优化跟着步骤可以把算法跑起来如果有几年优化算法经验重点看参数边界和收敛性指标计算部分的处理方式这两块是最容易踩坑的地方。2. NSGA-II 在 MATLAB 中的实现非支配排序、拥挤度距离与精英保留2.1 非支配排序把种群按 Pareto 支配关系分层NSGA-II 的核心第一步是对当前种群做非支配排序。假设种群规模为 N每个个体有 m 个目标函数值个体 A 支配个体 B 的条件是A 在所有目标上都不劣于 B且至少在一个目标上严格优于 B。排序的目标是把种群拆成 F1、F2、F3 等若干个前沿层F1 是当前种群中的非支配解集F2 是去掉 F1 后剩余个体中的非支配解集以此类推。MATLAB 工程实现里大多数人用双重循环逐个比较个体两两之间的支配关系复杂度是 O(MN²)M 是目标个数。种群规模在 100300 之间时这个开销可以接受但如果 N 超过 1000建议改成按目标维度先排序再剪枝的做法。下面给出一个标准实现function [F, rank] non_dominated_sort(pop, obj) % pop: 决策变量矩阵, 每行一个个体 % obj: 目标函数值矩阵, 每行对应 pop 中个体的 m 个目标 N size(pop, 1); dom_count zeros(N, 1); % 支配当前个体的数量 dominate_set cell(N, 1); % 当前个体支配的个体集合 F {}; % 前沿层集合, 每层存索引 for i 1:N for j 1:N if i j, continue; end if dominates(obj(i,:), obj(j,:)) dominate_set{i}(end1) j; elseif dominates(obj(j,:), obj(i,:)) dom_count(i) dom_count(i) 1; end end if dom_count(i) 0 F{1}(end1) i; end end k 1; while ~isempty(F{k}) Q []; for i F{k} for j dominate_set{i} dom_count(j) dom_count(j) - 1; if dom_count(j) 0 Q(end1) j; end end end k k 1; F{k} Q; end F(k:end) []; % 去掉末尾空层 rank zeros(N, 1); for i 1:length(F) rank(F{i}) i; end end function d dominates(a, b) d all(a b) any(a b); end这段代码的两个关键点在dom_count和dominate_set的用法上dom_count(i)记录有多少个体支配 i它为 0 的个体属于 F1每处理完一层就把该层个体从其余个体的支配计数中减掉新的 0 计数个体构成下一层。rank向量记录每个个体所属前沿的编号后面选择操作会直接拿它做比较。dominates函数里all(a b) any(a b)是严格意义上的 Pareto 支配判断如果想改成弱支配允许完全相等把any(a b)去掉即可但通常不建议这么做弱支配会让同一目标值重复的个体大量堆积在前沿上。注意这里F{k}中可能混入完全相同的个体决策变量和目标值都相同这是 NSGA-II 的已知问题标准做法是在排序前用unique按目标值去重否则后续拥挤度距离计算会产生大量 0 距离个体导致种群多样性退化。2.2 拥挤度距离维护多样性的关键计算非支配排序解决了「谁更优」的问题但同一前沿层内部的个体选择必须借助拥挤度距离。拥挤度距离的直观含义是某个个体周围有多少空间空间越大说明它越能代表该区域的分布性越应该被保留。对每个目标维度先把该前沿的个体按目标值排序边界个体的距离设为无穷大中间个体的距离用相邻个体的目标差除以该维度目标范围来归一化。function crowd crowding_distance(obj, idx) % idx: 某一前沿层个体的索引 layer_obj obj(idx, :); [m, n_obj] size(layer_obj); crowd zeros(m, 1); for j 1:n_obj [~, order] sort(layer_obj(:, j)); crowd(order(1)) inf; crowd(order(end)) inf; f_min layer_obj(order(1), j); f_max layer_obj(order(end), j); if f_max - f_min 1e-12 continue; % 该目标上所有值相同, 距离贡献为 0 end for k 2:m-1 crowd(order(k)) crowd(order(k)) ... (layer_obj(order(k1), j) - layer_obj(order(k-1), j)) / (f_max - f_min); end end end这段代码里容易出错的有两处。第一处是归一化分母的零值判断当某目标维度的所有值相等时直接除会得到 NaN处理方式是用f_max - f_min 1e-12跳过该维度。这在 ZDT1 这类测试函数上不常出现但换到真实工程问题、目标函数存在离散取值时几乎必然发生。第二处是边界个体设为 inf这样做保证了边界点在选择时永远优先被保留前沿端点不丢失。距离计算完成后选择操作就非常简单了比较两个个体时先看rankrank 小的胜出rank 相同时拥挤度距离大的胜出。这就是 NSGA-II 选择压力的全部来源替换 NSGA 的共享函数方法也正是它能让种群前沿均匀分布的根本原因。2.3 精英保留父代与子代合并后的截断策略NSGA-II 区别于早期遗传算法的核心改进在于精英保留策略。每一代生成子代后将父代种群 P_t 与子代种群 Q_t 合并成规模为 2N 的集合 R_t对 R_t 做非支配排序然后按前沿层从 F1 开始依次填入下一代种群 P_{t1}直到填满 N 个个体。最后一个能装下的前沿层内部用拥挤度距离从大到小选择把多余个体截断掉。function next_pop elitism_select(pop, obj, N) [F, rank] non_dominated_sort(pop, obj); next_pop []; for i 1:length(F) if length(next_pop) length(F{i}) N next_pop [next_pop; F{i}]; else need N - length(next_pop); idx F{i}; crowd crowding_distance(obj, idx); [~, order] sort(crowd, descend); next_pop [next_pop; idx(order(1:need))]; break; end end end这个实现的巧妙之处在于合并后的 2N 个体中最好的 N 个一定不会丢失因为即使某个优秀个体在子代生成过程中没有被选中参与交叉变异它也还存在于父代种群中合并后依然有资格被选入下一代。这也是 NSGA-II 收敛速度比 NSGA 快一截的直接原因。代码里的sort(crowd, descend)是升序取反的替代写法MATLAB 的 sort 函数在降序时注意descend参数拼写拼成descending会直接报错。3. 收敛性怎么判断NSGA-II 运行中的指标、绘图与参数调整3.1 种群进化过程可视化用前沿动态图观察收敛行为把 NSGA-II 跑起来后第一件事不是看最终结果而是观察每一代 Pareto 前沿的动态变化。最常见的做法是在主循环中每隔若干代调用一次绘图函数把当前种群的目标值散点画出来。对双目标问题一个简单的二维散点图就能暴露大部分收敛性问题前沿是否在向真实 Pareto 前沿移动、是否出现断层、边界是否收缩。function plot_front(obj, gen, n_gen) figure(1); plot(obj(:,1), obj(:,2), bo, MarkerSize, 4); xlabel(f_1); ylabel(f_2); title(sprintf(Generation %d / %d, gen, n_gen)); axis tight; grid on; drawnow; if mod(gen, 20) 0 saveas(gcf, sprintf(front_%04d.png, gen)); end end绘图函数里值得注意的细节是drawnow没有它 MATLAB 的 figure 窗口会在循环结束后才刷新动态观察效果完全丢失。saveas在迭代次数多的情况下会有 IO 开销通常每 20 代存一张图就够用了。另一个实践证明有效的做法是把真实 Pareto 前沿如果测试函数已知解析表达式画成实线作为参照这样收敛过程是否逼近真实前沿一目了然而不是只能看出「点在动」。判断收敛性的经验标准供参考当前沿连续多代一般 30 代以上的移动距离小于某个阈值时可视作进入停滞期。移动距离可以用相邻两代前沿个体之间的平均欧氏距离来度量也可以用 3.2 节的代际距离指标。3.2 收敛性量化指标IGD 与 GD 的计算脚本可视化只能定性判断定量判断需要收敛性指标。最常用的两个是 GDGenerational Distance世代距离和 IGDInverted Generational Distance反向世代距离。GD 度量的是算法求得的解集到真实 Pareto 前沿的最短距离平均值值越小说明解越接近真实前沿IGD 度量真实前沿上的每个点到算法解集的最短距离平均值它同时反映收敛性和分布性。function [gd, igd] compute_indicators(pf_approx, pf_true) % pf_approx: 算法得到的 Pareto 前沿 (n x m) % pf_true: 真实 Pareto 前沿的密集采样 (N x m), N 通常取 500~1000 n size(pf_approx, 1); N size(pf_true, 1); % GD: 对每个近似解找最近的真实前沿点 dist_app pdist2(pf_approx, pf_true); dmin_approx min(dist_app, [], 2); gd sqrt(mean(dmin_approx.^2)); % IGD: 对每个真实前沿点找最近的近似解 dist_true pdist2(pf_true, pf_approx); dmin_true min(dist_true, [], 2); igd sqrt(mean(dmin_true.^2)); endpdist2是 MATLAB 统计工具箱里的函数会自动计算两两欧氏距离比自己写双层循环快一个数量级。如果机器上没有工具箱用sqrt(bsxfun(plus, sum(a.^2,2), sum(b.^2,2)) - 2*a*b)也能达到同样效果。IGD 的计算对 pf_true 的采样密度很敏感真实前沿上采样点越多IGD 值越能反映真实的分布质量一般取 5001000 个点足够。跑实验时建议把每一代的 GD 和 IGD 值存到数组里迭代结束后用semilogy画对数坐标曲线收敛趋势会看得更清楚。GD 和 IGD 出现上升拐点往往意味着种群多样性崩溃即拥挤度距离计算或交叉变异参数出了问题。3.3 三个必须盯住的参数种群规模、交叉概率与变异概率参数设置是 NSGA-II 运行效果的分水岭也是最容易被网上下载的代码掩盖的部分。三个核心参数的取值范围和经验设置如下参数推荐范围常见默认值对收敛性的影响种群规模 N50500100太小导致前沿分层不充分IGD 指标高太大会拖慢每代排序速度交叉概率 Pc0.61.00.9过低时解空间探索不足代际距离停滞变异概率 Pm1/D 0.11/DD 为决策变量维度过大会破坏已收敛的前沿过小会早熟这里特别强调变异概率的计算方式。标准做法是取1/D其中 D 是决策变量的维度。对 ZDT1 这样的 30 维测试函数Pm 是 1/30 ≈ 0.033而有些博客干脆写死0.01这在低维问题上勉强能用到了高维问题基本等于没有变异算子。另一个工程经验是连续变量用模拟二进制交叉SBX和多项式变异不要用经典遗传算法的单点交叉和位翻转变异这两种算子是为离散编码设计的用于实数编码时收敛速度会明显变慢而且容易让解集中在局部区域。判断参数是否合理的快速方法是看 3.1 节的动态图如果前沿在 50 代内就冻结不再变化而 IGD 值还很高优先降低交叉概率、提高变异概率如果前沿始终分散无法形成连续曲线优先提高交叉概率和种群规模。4. NSGA-II 实战在 ZDT1 测试函数上的完整 MATLAB 运行代码4.1 问题定义与 ZDT1 目标函数ZDT1 是双目标测试问题中最经典的一个真实 Pareto 前沿解析形式为 f2 1 - sqrt(f1)f1 取值范围 [0,1]。它的决策变量维度通常取 30目标函数定义如下function [f1, f2] zdt1(x) % x: 决策变量向量, 长度 D n length(x); g 1 9 * sum(x(2:n)) / (n - 1); f1 x(1); f2 g * (1 - sqrt(f1)); endZDT1 的 g 函数设计使它存在大量局部前沿只有当 x(2) 到 x(n) 全部接近 0 时 g 才接近 1真实前沿在 f1 任意取值时对应 g1。测试一个 NSGA-II 实现的优劣ZDT1 是最快的判据跑 200 代如果 IGD 不能降到 1e-3 以下说明实现或参数大概率有问题。ZDT1 的前沿是凸的对算法相对友好建议跑通之后再换 ZDT2非凸前沿和 ZDT4大量局部前沿验证算法的极限。4.2 主循环初始化、锦标赛选择、交叉变异与迭代完整的 NSGA-II 主循环代码如下。这里使用 SBX 交叉和多项式变异两者是连续多目标优化的标准算子%% NSGA-II 主程序 - ZDT1 测试函数 clc; clear; close all; %% 参数设置 D 30; % 决策变量维度 N 100; % 种群规模 n_gen 250; % 迭代代数 Pc 0.9; % 交叉概率 Pm 1/D; % 变异概率 eta_c 20; % SBX 分布指数 eta_m 20; % 多项式变异分布指数 bounds zeros(D, 2); bounds(:, 1) 0; % 下界 bounds(:, 2) 1; % 上界 %% 初始化种群 pop rand(N, D) .* (bounds(:,2) - bounds(:,1)) bounds(:,1); %% 真实 Pareto 前沿采样 (用于 IGD 计算) pf_true linspace(0, 1, 500); pf_true [pf_true, 1 - sqrt(pf_true)]; %% 进化主循环 for gen 1:n_gen % 计算目标值 obj zeros(N, 2); for i 1:N [obj(i,1), obj(i,2)] zdt1(pop(i, :)); end % 锦标赛选择生成 mating pool [F, rank] non_dominated_sort(pop, obj); crowd zeros(N, 1); for k 1:length(F) crowd(F{k}) crowding_distance(obj, F{k}); end pool zeros(N, D); for i 1:N idx1 randi(N); idx2 randi(N); if rank(idx1) rank(idx2) pool(i, :) pop(idx1, :); elseif rank(idx1) rank(idx2) pool(i, :) pop(idx2, :); else if crowd(idx1) crowd(idx2) pool(i, :) pop(idx1, :); else pool(i, :) pop(idx2, :); end end end % SBX 交叉和多项式变异生成子代 offspring zeros(N, D); for i 1:2:N p1 pool(i, :); p2 pool(i1, :); if rand Pc [c1, c2] sbx_crossover(p1, p2, eta_c, bounds); else c1 p1; c2 p2; end offspring(i, :) mutation_poly(c1, Pm, eta_m, bounds); offspring(i1, :) mutation_poly(c2, Pm, eta_m, bounds); end % 合并父代与子代, 精英保留 combined_pop [pop; offspring]; combined_obj zeros(2*N, 2); for i 1:2*N [combined_obj(i,1), combined_obj(i,2)] zdt1(combined_pop(i, :)); end sel_idx elitism_select(combined_pop, combined_obj, N); pop combined_pop(sel_idx, :); % 每 10 代输出一次指标 if mod(gen, 10) 0 obj_now zeros(N, 2); for i 1:N [obj_now(i,1), obj_now(i,2)] zdt1(pop(i, :)); end [gd, igd] compute_indicators(obj_now, pf_true); fprintf(Gen %3d | GD: %.4e | IGD: %.4e\n, gen, gd, igd); end end %% SBX 交叉算子 function [c1, c2] sbx_crossover(p1, p2, eta_c, bounds) c1 p1; c2 p2; for j 1:length(p1) if rand 0.5 if abs(p1(j) - p2(j)) 1e-14 if p1(j) p2(j) x1 p1(j); x2 p2(j); else x1 p2(j); x2 p1(j); end u rand; if u 0.5 beta (2*u)^(1/(eta_c1)); else beta (1/(2*(1-u)))^(1/(eta_c1)); end c1(j) 0.5 * ((1beta)*x1 (1-beta)*x2); c2(j) 0.5 * ((1-beta)*x1 (1beta)*x2); % 边界截断 c1(j) min(max(c1(j), bounds(j,1)), bounds(j,2)); c2(j) min(max(c2(j), bounds(j,1)), bounds(j,2)); end end end end %% 多项式变异算子 function c mutation_poly(p, Pm, eta_m, bounds) c p; for j 1:length(p) if rand Pm u rand; delta 0; if u 0.5 delta (2*u)^(1/(eta_m1)) - 1; else delta 1 - (2*(1-u))^(1/(eta_m1)); end c(j) p(j) delta * (bounds(j,2) - bounds(j,1)); c(j) min(max(c(j), bounds(j,1)), bounds(j,2)); end end end这段代码的主循环执行顺序值得细看先算目标值做非支配排序和拥挤度计算然后锦标赛选择生成交配池交配池中的个体两两配对执行 SBX 交叉生成子代再对子代做多项式变异。子代生成后被合并进父代种群调用 2.3 节的elitism_select完成精英保留截断后进入下一代。fprintf每 10 代输出一次 GD 和 IGD便于实时观察收敛进程。锦标赛选择里用了 rank 优先、拥挤度次之的比较策略这和选择压力直接挂钩。如果randi(N)选出的两个个体恰好是同一个索引代码也不会出错因为它会再次比较该个体与自身的 rank 和 crowd结果仍是自己这在种群规模较大时影响可以忽略。SBX 交叉的实现有两个关键参数eta_c控制子代偏离父代的平均程度取值越大子代越接近父代分布指数 20 是原论文中的默认值beta的计算分两段确保子代在父代两侧的分布是对称的。多项式变异里delta的计算本质是把随机数映射到决策空间上的扰动幅度eta_m控制扰动大小的衰减速度。两个算子的末尾都有边界截断原因是算子在变量接近边界时可能产生越界值如果不截断会导致 zdt1 函数里出现超出定义域的输入。4.3 实验结果解读跑 250 代后怎么看收敛曲线在 MATLAB R2023b 上这个配置跑 250 代大约耗时 2040 秒取决于机器性能。输出内容类似下面这样Gen 10 | GD: 8.23e-03 | IGD: 1.91e-02 Gen 50 | GD: 1.56e-03 | IGD: 3.87e-03 Gen 100 | GD: 4.12e-04 | IGD: 1.03e-03 Gen 150 | GD: 2.08e-04 | IGD: 5.21e-04 Gen 200 | GD: 1.45e-04 | IGD: 3.66e-04 Gen 250 | GD: 1.21e-04 | IGD: 3.02e-04关键看两个规律GD 和 IGD 是否单调下降下降速度是否越来越慢。前 50 代下降最快这是算法从随机分布到接近前沿的快速收敛期100 代以后进入精细调整期指标下降趋缓但不应反弹。如果发现 IGD 在第 80 代以后出现上升优先检查精英保留选择出的个体是否正确填充了下一代——一个常见错误是截断时按拥挤度升序排序把最有价值的边界个体删掉了。5. 收敛性验证技巧把超体积指标和参考前沿结合判断算法是否真的收敛5.1 用 Hypervolume 指标验证收敛性与多样性IGD 和 GD 都需要真实 Pareto 前沿作为参照但真实工程问题的前沿是未知的。这时超体积指标Hypervolume, HV是唯一不依赖参考前沿的收敛性指标它以参考点通常取每个目标的 worst 值或稍大值为顶点计算算法解集在目标空间中覆盖的体积。HV 越大说明解集既接近真实前沿覆盖范围广又分布均匀体积没有空洞。对双目标问题HV 就是解集点与参考点围成的面积总和。function hv hypervolume(obj, ref_point) % obj: 算法得到的 Pareto 前沿点 (n x m) % ref_point: 参考点向量, 长度 m, 取各目标上界 % 实现思路: 按第一个目标排序后, 用梯形积分近似计算 pts sortrows(obj, 1); pts unique(pts, rows); % 去重, 防止重复点影响面积 n size(pts, 1); hv 0; prev_f2 ref_point(2); for i n:-1:1 f1_i pts(i, 1); if pts(i, 2) prev_f2 f2_i pts(i, 2); else f2_i prev_f2; % 非支配点修正, 保证面积不重叠 end if i n hv hv (ref_point(1) - f1_i) * (ref_point(2) - f2_i); else hv hv (pts(i1, 1) - f1_i) * (ref_point(2) - f2_i); end prev_f2 f2_i; end end这个梯形积分实现的关键在else分支遍历时从最右侧的点开始向左推进每一步用当前点的 f2 与「当前已覆盖区域的上边线」的较小值作为高度防止解集中非支配关系不严格时产生的矩形重叠。unique(pts, rows)去掉重复点否则重复点会在同位置产生零宽度的矩形虽不影响最终值但浪费计算。参考点的选择对 HV 值影响显著工程上通常取各目标在初始种群中的最大值再乘 1.1 作为参考点既保证能够覆盖整个区域又不会因参考点过远而掩盖分布不均的问题。5.2 多轮独立运行与指标统计避免单次运行误导NSGA-II 是随机算法单次运行的指标没有任何说服力。正式评估一个实现或一组参数的收敛性至少要独立运行 20 次记录每次的 IGD 和 HV 终值输出均值和标准差。标准差大说明算法稳定性差多半是变异概率过低或种群规模不足均值差而标准差小则是算法本身收敛能力不够。n_runs 20; igd_all zeros(n_runs, 1); hv_all zeros(n_runs, 1); for r 1:n_runs % 重新初始化种群并运行主循环 (略) % 结束后计算当前种群的 IGD 和 HV igd_all(r) igd_final; hv_all(r) hv_final; end fprintf(IGD: mean%.4e std%.4e\n, mean(igd_all), std(igd_all)); fprintf(HV: mean%.4e std%.4e\n, mean(hv_all), std(hv_all));实际做实验时建议同时输出 5.1 节的 HV 和 3.2 节的 IGD两者配合使用HV 不需要真实前沿适合工程问题IGD 需要真实前沿适合算法对比和测试函数验证。如果 HV 曲线在第 100 代后仍然明显上升说明算法尚未收敛需要增加迭代次数或调整参数。反之如果 HV 在很早就进入平台期但数值明显低于同类实现的文献结果就要回头检查非支配排序和精英保留的代码逻辑这两个函数一旦写错后面的所有优化都白费。最终验证通过的配置把随机种子固定下来再跑一次用 5.1 节的绘图函数保存最后一幅前沿图作为后续调参的基线对照。本文还有配套的精品资源点击获取 
简介面向多目标优化研究者与 MATLAB 学习者的 NSGA-II 算法实现压缩包针对多个冲突目标同时优化问题提供一套完整可运行的遗传算法框架。压缩包共 13 个文件以 12 个 m 文件为主涵盖非支配排序、拥挤距离计算、锦标赛选择、交叉、变异及帕累托前沿绘图等核心模块另含 1 个来源链接文件整体大小仅 6KB结构紧凑、便于阅读。当前已有 989 人学习适合入门与进阶。通过这套 MATLAB 源码读者可以直观观察 NSGA-II 的迭代收敛过程理解快速非支配排序如何划分解的质量层次、拥挤距离如何维持种群多样性并可根据自身优化目标修改适应度函数、调整种群参数快速迁移到工程设计、投资决策等实际场景中。模块化的函数设计也便于逐块调试与教学演示作为多目标进化算法学习与二次开发的轻量工具。1. NSGA-II 多目标优化在 MATLAB 里跑通收敛性的关键从来不是把代码复制下来NSGA-IINon-dominated Sorting Genetic Algorithm II是过去二十年多目标优化领域被引用最多的算法骨架。它用非支配排序维持解的层级用拥挤度距离保证种群在 Pareto 前沿上铺得开再用精英保留策略让每一代的最优个体不会丢失。MATLAB 里实现 NSGA-II 的代码在网上有成百上千份但多数人下载后第一轮跑出来的结果不是前沿残缺就是收敛曲线震荡甚至种群直接塌缩到某个单目标最优附近——这恰恰说明问题不在算法本身而在运行参数和收敛性判断这两件事上没被认真对待。这篇文章围绕「多目标运行」和「收敛性」两个核心词展开先拆 NSGA-II 在 MATLAB 中的完整实现链路再给出可直接运行的代码和参数表最后落到用 Hypervolume、IGD 等指标量化验证收敛性的具体脚本。读者如果是刚接触多目标优化跟着步骤可以把算法跑起来如果有几年优化算法经验重点看参数边界和收敛性指标计算部分的处理方式这两块是最容易踩坑的地方。2. NSGA-II 在 MATLAB 中的实现非支配排序、拥挤度距离与精英保留2.1 非支配排序把种群按 Pareto 支配关系分层NSGA-II 的核心第一步是对当前种群做非支配排序。假设种群规模为 N每个个体有 m 个目标函数值个体 A 支配个体 B 的条件是A 在所有目标上都不劣于 B且至少在一个目标上严格优于 B。排序的目标是把种群拆成 F1、F2、F3 等若干个前沿层F1 是当前种群中的非支配解集F2 是去掉 F1 后剩余个体中的非支配解集以此类推。MATLAB 工程实现里大多数人用双重循环逐个比较个体两两之间的支配关系复杂度是 O(MN²)M 是目标个数。种群规模在 100300 之间时这个开销可以接受但如果 N 超过 1000建议改成按目标维度先排序再剪枝的做法。下面给出一个标准实现function [F, rank] non_dominated_sort(pop, obj) % pop: 决策变量矩阵, 每行一个个体 % obj: 目标函数值矩阵, 每行对应 pop 中个体的 m 个目标 N size(pop, 1); dom_count zeros(N, 1); % 支配当前个体的数量 dominate_set cell(N, 1); % 当前个体支配的个体集合 F {}; % 前沿层集合, 每层存索引 for i 1:N for j 1:N if i j, continue; end if dominates(obj(i,:), obj(j,:)) dominate_set{i}(end1) j; elseif dominates(obj(j,:), obj(i,:)) dom_count(i) dom_count(i) 1; end end if dom_count(i) 0 F{1}(end1) i; end end k 1; while ~isempty(F{k}) Q []; for i F{k} for j dominate_set{i} dom_count(j) dom_count(j) - 1; if dom_count(j) 0 Q(end1) j; end end end k k 1; F{k} Q; end F(k:end) []; % 去掉末尾空层 rank zeros(N, 1); for i 1:length(F) rank(F{i}) i; end end function d dominates(a, b) d all(a b) any(a b); end这段代码的两个关键点在dom_count和dominate_set的用法上dom_count(i)记录有多少个体支配 i它为 0 的个体属于 F1每处理完一层就把该层个体从其余个体的支配计数中减掉新的 0 计数个体构成下一层。rank向量记录每个个体所属前沿的编号后面选择操作会直接拿它做比较。dominates函数里all(a b) any(a b)是严格意义上的 Pareto 支配判断如果想改成弱支配允许完全相等把any(a b)去掉即可但通常不建议这么做弱支配会让同一目标值重复的个体大量堆积在前沿上。注意这里F{k}中可能混入完全相同的个体决策变量和目标值都相同这是 NSGA-II 的已知问题标准做法是在排序前用unique按目标值去重否则后续拥挤度距离计算会产生大量 0 距离个体导致种群多样性退化。2.2 拥挤度距离维护多样性的关键计算非支配排序解决了「谁更优」的问题但同一前沿层内部的个体选择必须借助拥挤度距离。拥挤度距离的直观含义是某个个体周围有多少空间空间越大说明它越能代表该区域的分布性越应该被保留。对每个目标维度先把该前沿的个体按目标值排序边界个体的距离设为无穷大中间个体的距离用相邻个体的目标差除以该维度目标范围来归一化。function crowd crowding_distance(obj, idx) % idx: 某一前沿层个体的索引 layer_obj obj(idx, :); [m, n_obj] size(layer_obj); crowd zeros(m, 1); for j 1:n_obj [~, order] sort(layer_obj(:, j)); crowd(order(1)) inf; crowd(order(end)) inf; f_min layer_obj(order(1), j); f_max layer_obj(order(end), j); if f_max - f_min 1e-12 continue; % 该目标上所有值相同, 距离贡献为 0 end for k 2:m-1 crowd(order(k)) crowd(order(k)) ... (layer_obj(order(k1), j) - layer_obj(order(k-1), j)) / (f_max - f_min); end end end这段代码里容易出错的有两处。第一处是归一化分母的零值判断当某目标维度的所有值相等时直接除会得到 NaN处理方式是用f_max - f_min 1e-12跳过该维度。这在 ZDT1 这类测试函数上不常出现但换到真实工程问题、目标函数存在离散取值时几乎必然发生。第二处是边界个体设为 inf这样做保证了边界点在选择时永远优先被保留前沿端点不丢失。距离计算完成后选择操作就非常简单了比较两个个体时先看rankrank 小的胜出rank 相同时拥挤度距离大的胜出。这就是 NSGA-II 选择压力的全部来源替换 NSGA 的共享函数方法也正是它能让种群前沿均匀分布的根本原因。2.3 精英保留父代与子代合并后的截断策略NSGA-II 区别于早期遗传算法的核心改进在于精英保留策略。每一代生成子代后将父代种群 P_t 与子代种群 Q_t 合并成规模为 2N 的集合 R_t对 R_t 做非支配排序然后按前沿层从 F1 开始依次填入下一代种群 P_{t1}直到填满 N 个个体。最后一个能装下的前沿层内部用拥挤度距离从大到小选择把多余个体截断掉。function next_pop elitism_select(pop, obj, N) [F, rank] non_dominated_sort(pop, obj); next_pop []; for i 1:length(F) if length(next_pop) length(F{i}) N next_pop [next_pop; F{i}]; else need N - length(next_pop); idx F{i}; crowd crowding_distance(obj, idx); [~, order] sort(crowd, descend); next_pop [next_pop; idx(order(1:need))]; break; end end end这个实现的巧妙之处在于合并后的 2N 个体中最好的 N 个一定不会丢失因为即使某个优秀个体在子代生成过程中没有被选中参与交叉变异它也还存在于父代种群中合并后依然有资格被选入下一代。这也是 NSGA-II 收敛速度比 NSGA 快一截的直接原因。代码里的sort(crowd, descend)是升序取反的替代写法MATLAB 的 sort 函数在降序时注意descend参数拼写拼成descending会直接报错。3. 收敛性怎么判断NSGA-II 运行中的指标、绘图与参数调整3.1 种群进化过程可视化用前沿动态图观察收敛行为把 NSGA-II 跑起来后第一件事不是看最终结果而是观察每一代 Pareto 前沿的动态变化。最常见的做法是在主循环中每隔若干代调用一次绘图函数把当前种群的目标值散点画出来。对双目标问题一个简单的二维散点图就能暴露大部分收敛性问题前沿是否在向真实 Pareto 前沿移动、是否出现断层、边界是否收缩。function plot_front(obj, gen, n_gen) figure(1); plot(obj(:,1), obj(:,2), bo, MarkerSize, 4); xlabel(f_1); ylabel(f_2); title(sprintf(Generation %d / %d, gen, n_gen)); axis tight; grid on; drawnow; if mod(gen, 20) 0 saveas(gcf, sprintf(front_%04d.png, gen)); end end绘图函数里值得注意的细节是drawnow没有它 MATLAB 的 figure 窗口会在循环结束后才刷新动态观察效果完全丢失。saveas在迭代次数多的情况下会有 IO 开销通常每 20 代存一张图就够用了。另一个实践证明有效的做法是把真实 Pareto 前沿如果测试函数已知解析表达式画成实线作为参照这样收敛过程是否逼近真实前沿一目了然而不是只能看出「点在动」。判断收敛性的经验标准供参考当前沿连续多代一般 30 代以上的移动距离小于某个阈值时可视作进入停滞期。移动距离可以用相邻两代前沿个体之间的平均欧氏距离来度量也可以用 3.2 节的代际距离指标。3.2 收敛性量化指标IGD 与 GD 的计算脚本可视化只能定性判断定量判断需要收敛性指标。最常用的两个是 GDGenerational Distance世代距离和 IGDInverted Generational Distance反向世代距离。GD 度量的是算法求得的解集到真实 Pareto 前沿的最短距离平均值值越小说明解越接近真实前沿IGD 度量真实前沿上的每个点到算法解集的最短距离平均值它同时反映收敛性和分布性。function [gd, igd] compute_indicators(pf_approx, pf_true) % pf_approx: 算法得到的 Pareto 前沿 (n x m) % pf_true: 真实 Pareto 前沿的密集采样 (N x m), N 通常取 500~1000 n size(pf_approx, 1); N size(pf_true, 1); % GD: 对每个近似解找最近的真实前沿点 dist_app pdist2(pf_approx, pf_true); dmin_approx min(dist_app, [], 2); gd sqrt(mean(dmin_approx.^2)); % IGD: 对每个真实前沿点找最近的近似解 dist_true pdist2(pf_true, pf_approx); dmin_true min(dist_true, [], 2); igd sqrt(mean(dmin_true.^2)); endpdist2是 MATLAB 统计工具箱里的函数会自动计算两两欧氏距离比自己写双层循环快一个数量级。如果机器上没有工具箱用sqrt(bsxfun(plus, sum(a.^2,2), sum(b.^2,2)) - 2*a*b)也能达到同样效果。IGD 的计算对 pf_true 的采样密度很敏感真实前沿上采样点越多IGD 值越能反映真实的分布质量一般取 5001000 个点足够。跑实验时建议把每一代的 GD 和 IGD 值存到数组里迭代结束后用semilogy画对数坐标曲线收敛趋势会看得更清楚。GD 和 IGD 出现上升拐点往往意味着种群多样性崩溃即拥挤度距离计算或交叉变异参数出了问题。3.3 三个必须盯住的参数种群规模、交叉概率与变异概率参数设置是 NSGA-II 运行效果的分水岭也是最容易被网上下载的代码掩盖的部分。三个核心参数的取值范围和经验设置如下参数推荐范围常见默认值对收敛性的影响种群规模 N50500100太小导致前沿分层不充分IGD 指标高太大会拖慢每代排序速度交叉概率 Pc0.61.00.9过低时解空间探索不足代际距离停滞变异概率 Pm1/D 0.11/DD 为决策变量维度过大会破坏已收敛的前沿过小会早熟这里特别强调变异概率的计算方式。标准做法是取1/D其中 D 是决策变量的维度。对 ZDT1 这样的 30 维测试函数Pm 是 1/30 ≈ 0.033而有些博客干脆写死0.01这在低维问题上勉强能用到了高维问题基本等于没有变异算子。另一个工程经验是连续变量用模拟二进制交叉SBX和多项式变异不要用经典遗传算法的单点交叉和位翻转变异这两种算子是为离散编码设计的用于实数编码时收敛速度会明显变慢而且容易让解集中在局部区域。判断参数是否合理的快速方法是看 3.1 节的动态图如果前沿在 50 代内就冻结不再变化而 IGD 值还很高优先降低交叉概率、提高变异概率如果前沿始终分散无法形成连续曲线优先提高交叉概率和种群规模。4. NSGA-II 实战在 ZDT1 测试函数上的完整 MATLAB 运行代码4.1 问题定义与 ZDT1 目标函数ZDT1 是双目标测试问题中最经典的一个真实 Pareto 前沿解析形式为 f2 1 - sqrt(f1)f1 取值范围 [0,1]。它的决策变量维度通常取 30目标函数定义如下function [f1, f2] zdt1(x) % x: 决策变量向量, 长度 D n length(x); g 1 9 * sum(x(2:n)) / (n - 1); f1 x(1); f2 g * (1 - sqrt(f1)); endZDT1 的 g 函数设计使它存在大量局部前沿只有当 x(2) 到 x(n) 全部接近 0 时 g 才接近 1真实前沿在 f1 任意取值时对应 g1。测试一个 NSGA-II 实现的优劣ZDT1 是最快的判据跑 200 代如果 IGD 不能降到 1e-3 以下说明实现或参数大概率有问题。ZDT1 的前沿是凸的对算法相对友好建议跑通之后再换 ZDT2非凸前沿和 ZDT4大量局部前沿验证算法的极限。4.2 主循环初始化、锦标赛选择、交叉变异与迭代完整的 NSGA-II 主循环代码如下。这里使用 SBX 交叉和多项式变异两者是连续多目标优化的标准算子%% NSGA-II 主程序 - ZDT1 测试函数 clc; clear; close all; %% 参数设置 D 30; % 决策变量维度 N 100; % 种群规模 n_gen 250; % 迭代代数 Pc 0.9; % 交叉概率 Pm 1/D; % 变异概率 eta_c 20; % SBX 分布指数 eta_m 20; % 多项式变异分布指数 bounds zeros(D, 2); bounds(:, 1) 0; % 下界 bounds(:, 2) 1; % 上界 %% 初始化种群 pop rand(N, D) .* (bounds(:,2) - bounds(:,1)) bounds(:,1); %% 真实 Pareto 前沿采样 (用于 IGD 计算) pf_true linspace(0, 1, 500); pf_true [pf_true, 1 - sqrt(pf_true)]; %% 进化主循环 for gen 1:n_gen % 计算目标值 obj zeros(N, 2); for i 1:N [obj(i,1), obj(i,2)] zdt1(pop(i, :)); end % 锦标赛选择生成 mating pool [F, rank] non_dominated_sort(pop, obj); crowd zeros(N, 1); for k 1:length(F) crowd(F{k}) crowding_distance(obj, F{k}); end pool zeros(N, D); for i 1:N idx1 randi(N); idx2 randi(N); if rank(idx1) rank(idx2) pool(i, :) pop(idx1, :); elseif rank(idx1) rank(idx2) pool(i, :) pop(idx2, :); else if crowd(idx1) crowd(idx2) pool(i, :) pop(idx1, :); else pool(i, :) pop(idx2, :); end end end % SBX 交叉和多项式变异生成子代 offspring zeros(N, D); for i 1:2:N p1 pool(i, :); p2 pool(i1, :); if rand Pc [c1, c2] sbx_crossover(p1, p2, eta_c, bounds); else c1 p1; c2 p2; end offspring(i, :) mutation_poly(c1, Pm, eta_m, bounds); offspring(i1, :) mutation_poly(c2, Pm, eta_m, bounds); end % 合并父代与子代, 精英保留 combined_pop [pop; offspring]; combined_obj zeros(2*N, 2); for i 1:2*N [combined_obj(i,1), combined_obj(i,2)] zdt1(combined_pop(i, :)); end sel_idx elitism_select(combined_pop, combined_obj, N); pop combined_pop(sel_idx, :); % 每 10 代输出一次指标 if mod(gen, 10) 0 obj_now zeros(N, 2); for i 1:N [obj_now(i,1), obj_now(i,2)] zdt1(pop(i, :)); end [gd, igd] compute_indicators(obj_now, pf_true); fprintf(Gen %3d | GD: %.4e | IGD: %.4e\n, gen, gd, igd); end end %% SBX 交叉算子 function [c1, c2] sbx_crossover(p1, p2, eta_c, bounds) c1 p1; c2 p2; for j 1:length(p1) if rand 0.5 if abs(p1(j) - p2(j)) 1e-14 if p1(j) p2(j) x1 p1(j); x2 p2(j); else x1 p2(j); x2 p1(j); end u rand; if u 0.5 beta (2*u)^(1/(eta_c1)); else beta (1/(2*(1-u)))^(1/(eta_c1)); end c1(j) 0.5 * ((1beta)*x1 (1-beta)*x2); c2(j) 0.5 * ((1-beta)*x1 (1beta)*x2); % 边界截断 c1(j) min(max(c1(j), bounds(j,1)), bounds(j,2)); c2(j) min(max(c2(j), bounds(j,1)), bounds(j,2)); end end end end %% 多项式变异算子 function c mutation_poly(p, Pm, eta_m, bounds) c p; for j 1:length(p) if rand Pm u rand; delta 0; if u 0.5 delta (2*u)^(1/(eta_m1)) - 1; else delta 1 - (2*(1-u))^(1/(eta_m1)); end c(j) p(j) delta * (bounds(j,2) - bounds(j,1)); c(j) min(max(c(j), bounds(j,1)), bounds(j,2)); end end end这段代码的主循环执行顺序值得细看先算目标值做非支配排序和拥挤度计算然后锦标赛选择生成交配池交配池中的个体两两配对执行 SBX 交叉生成子代再对子代做多项式变异。子代生成后被合并进父代种群调用 2.3 节的elitism_select完成精英保留截断后进入下一代。fprintf每 10 代输出一次 GD 和 IGD便于实时观察收敛进程。锦标赛选择里用了 rank 优先、拥挤度次之的比较策略这和选择压力直接挂钩。如果randi(N)选出的两个个体恰好是同一个索引代码也不会出错因为它会再次比较该个体与自身的 rank 和 crowd结果仍是自己这在种群规模较大时影响可以忽略。SBX 交叉的实现有两个关键参数eta_c控制子代偏离父代的平均程度取值越大子代越接近父代分布指数 20 是原论文中的默认值beta的计算分两段确保子代在父代两侧的分布是对称的。多项式变异里delta的计算本质是把随机数映射到决策空间上的扰动幅度eta_m控制扰动大小的衰减速度。两个算子的末尾都有边界截断原因是算子在变量接近边界时可能产生越界值如果不截断会导致 zdt1 函数里出现超出定义域的输入。4.3 实验结果解读跑 250 代后怎么看收敛曲线在 MATLAB R2023b 上这个配置跑 250 代大约耗时 2040 秒取决于机器性能。输出内容类似下面这样Gen 10 | GD: 8.23e-03 | IGD: 1.91e-02 Gen 50 | GD: 1.56e-03 | IGD: 3.87e-03 Gen 100 | GD: 4.12e-04 | IGD: 1.03e-03 Gen 150 | GD: 2.08e-04 | IGD: 5.21e-04 Gen 200 | GD: 1.45e-04 | IGD: 3.66e-04 Gen 250 | GD: 1.21e-04 | IGD: 3.02e-04关键看两个规律GD 和 IGD 是否单调下降下降速度是否越来越慢。前 50 代下降最快这是算法从随机分布到接近前沿的快速收敛期100 代以后进入精细调整期指标下降趋缓但不应反弹。如果发现 IGD 在第 80 代以后出现上升优先检查精英保留选择出的个体是否正确填充了下一代——一个常见错误是截断时按拥挤度升序排序把最有价值的边界个体删掉了。5. 收敛性验证技巧把超体积指标和参考前沿结合判断算法是否真的收敛5.1 用 Hypervolume 指标验证收敛性与多样性IGD 和 GD 都需要真实 Pareto 前沿作为参照但真实工程问题的前沿是未知的。这时超体积指标Hypervolume, HV是唯一不依赖参考前沿的收敛性指标它以参考点通常取每个目标的 worst 值或稍大值为顶点计算算法解集在目标空间中覆盖的体积。HV 越大说明解集既接近真实前沿覆盖范围广又分布均匀体积没有空洞。对双目标问题HV 就是解集点与参考点围成的面积总和。function hv hypervolume(obj, ref_point) % obj: 算法得到的 Pareto 前沿点 (n x m) % ref_point: 参考点向量, 长度 m, 取各目标上界 % 实现思路: 按第一个目标排序后, 用梯形积分近似计算 pts sortrows(obj, 1); pts unique(pts, rows); % 去重, 防止重复点影响面积 n size(pts, 1); hv 0; prev_f2 ref_point(2); for i n:-1:1 f1_i pts(i, 1); if pts(i, 2) prev_f2 f2_i pts(i, 2); else f2_i prev_f2; % 非支配点修正, 保证面积不重叠 end if i n hv hv (ref_point(1) - f1_i) * (ref_point(2) - f2_i); else hv hv (pts(i1, 1) - f1_i) * (ref_point(2) - f2_i); end prev_f2 f2_i; end end这个梯形积分实现的关键在else分支遍历时从最右侧的点开始向左推进每一步用当前点的 f2 与「当前已覆盖区域的上边线」的较小值作为高度防止解集中非支配关系不严格时产生的矩形重叠。unique(pts, rows)去掉重复点否则重复点会在同位置产生零宽度的矩形虽不影响最终值但浪费计算。参考点的选择对 HV 值影响显著工程上通常取各目标在初始种群中的最大值再乘 1.1 作为参考点既保证能够覆盖整个区域又不会因参考点过远而掩盖分布不均的问题。5.2 多轮独立运行与指标统计避免单次运行误导NSGA-II 是随机算法单次运行的指标没有任何说服力。正式评估一个实现或一组参数的收敛性至少要独立运行 20 次记录每次的 IGD 和 HV 终值输出均值和标准差。标准差大说明算法稳定性差多半是变异概率过低或种群规模不足均值差而标准差小则是算法本身收敛能力不够。n_runs 20; igd_all zeros(n_runs, 1); hv_all zeros(n_runs, 1); for r 1:n_runs % 重新初始化种群并运行主循环 (略) % 结束后计算当前种群的 IGD 和 HV igd_all(r) igd_final; hv_all(r) hv_final; end fprintf(IGD: mean%.4e std%.4e\n, mean(igd_all), std(igd_all)); fprintf(HV: mean%.4e std%.4e\n, mean(hv_all), std(hv_all));实际做实验时建议同时输出 5.1 节的 HV 和 3.2 节的 IGD两者配合使用HV 不需要真实前沿适合工程问题IGD 需要真实前沿适合算法对比和测试函数验证。如果 HV 曲线在第 100 代后仍然明显上升说明算法尚未收敛需要增加迭代次数或调整参数。反之如果 HV 在很早就进入平台期但数值明显低于同类实现的文献结果就要回头检查非支配排序和精英保留的代码逻辑这两个函数一旦写错后面的所有优化都白费。最终验证通过的配置把随机种子固定下来再跑一次用 5.1 节的绘图函数保存最后一幅前沿图作为后续调参的基线对照。本文还有配套的精品资源点击获取