MATLAB实现NSGA-II多目标优化算法：非支配排序与拥挤距离详解

简介MATLAB环境下的NSGA-II算法实现面向需要解决多目标优化问题的科研与工程人员。资源以中文注释详细拆解了非支配排序遗传算法第二代的完整流程包括种群初始化、非支配排序、拥挤距离选择、精英保留策略、交叉变异等核心模块并附带自适应变异、正态分布标准差等改进对比实验便于理解算法原理与调参思路。资源共20个文件以11个m源码文件为主涵盖主程序、遗传算子、非支配排序、目标函数评估等关键函数另有7个bmp结果图直观展示收敛分布与改进效果1个xls保存实验数据1个txt提供运行说明。压缩包整体2.44MB轻量易部署。目前已有20880人学习下载口碑验证。通过运行代码可观察算法在特定问题上的帕累托前沿演化同时中文注释极大降低了上手门槛适合作为多目标优化入门及二次开发的参考脚本。 先小吐槽一下标题里的 MTALAB 我默认是按 MATLAB 来读。最近好几个人拿着“MTALAB NSGA2算法”这个关键词来问说网上讲 NSGA2 概念的 PPT 一堆但真正能在 MATLAB 里跑出 Pareto 前沿、能直接拿来改的代码特别少。这篇文章就按我自己的实现习惯把 NSGA-II 从原理到代码再到调参和踩坑完整讲一遍并以 ZDT1 测试函数做例子。适合课程作业、论文仿真和工程项目起步参考不需要额外工具箱基础版 MATLAB 就能跑。1. NSGA-II 的核心思路为什么值得在 MATLAB 里手写一遍1.1 多目标优化问题的本质没有唯一最优解单目标优化求的是一个最小值或最大值多目标优化不一样。比如设计一个结构想同时让重量最轻、刚度最大、成本最低这几个目标往往是冲突的重量轻了成本可能上去刚度大了重量也可能上去。最终得到的不该是一个点而是一组解这组解里没有一个解能在所有目标上同时碾压其他解。这里要引出支配关系如果一个解 A 在所有目标上都不比解 B 差且至少有一个目标严格优于 B就说 A 支配 B。所有不被其他解支配的解构成第一层非支配前沿也就是 Pareto 前沿。NSGA-II 的目标就是尽量得到一组分布均匀、覆盖面广的 Pareto 前沿解而不是只盯着某一个目标。1.2 两个核心机制非支配排序和拥挤距离NSGA-II 名字里的 Non-dominated Sorting Genetic Algorithm核心就是“非支配排序”。它把种群里的个体按支配关系分成好几层第一层是当前种群的非支配解第二层是去掉第一层后剩余个体里的非支配解以此类推。优化过程中越靠前的层意味着质量越高应该优先保留。但光分层还不够。同一层里可能有几十个个体到底留谁如果随便选最后结果容易挤在 Pareto 前沿的某一小段上。NSGA-II 用拥挤距离来解决把某个个体前后相邻个体在每个目标上的距离归一化后累加距离越大说明这个解周围越空越值得保留。简单说就是保证解的多样性别让它们全挤一块。和第一代 NSGA 比NSGA-II 最大的改进是去掉了共享参数引入了拥挤距离和精英保留机制。快速非支配排序的复杂度是 O(MN²)M 是目标数N 是种群规模在常见规模下完全能接受。也正因为这两个机制都很清楚用 MATLAB 从头实现一遍一点也不难。1.3 为什么选 MATLAB 实现MATLAB 做这个事最大的优势是矩阵思维。一个种群 N 个个体、每个个体有 dim 个决策变量直接用一个 N×dim 矩阵存目标值用一个 N×M 矩阵存处理起来非常顺手。排序、比较、索引都有一堆现成函数画图也方便算完一个 plot 就能看到 Pareto 前沿。还有一个现实好处不需要安装额外工具箱。网上关于 MATLAB 工具包的讨论很多什么下载安装、工具箱激活、版本密钥之类的搞得很复杂。但手写 NSGA-II 只需要基础 MATLAB 环境连 Optimization Toolbox 都不需要核心算法代码自己写了反而更可控。2. 环境准备与主循环框架2.1 基础环境要求建议用 R2016b 及以上版本因为从这一版开始脚本文件里可以写局部函数方便把所有代码放在一个文件里。如果版本更老就把每个函数单独存成 .m 文件名称和函数名保持一致。自定义函数我建议这样组织evaluate_zdt1.m计算目标函数non_dominated_sort.m快速非支配排序crowding_distance.m计算拥挤距离sbx_crossover.m模拟二进制交叉poly_mutation.m多项式变异主脚本负责初始化、循环和结果展示。这样结构清晰改起来也方便。2.2 NSGA-II 主循环逻辑主循环其实不复杂每次迭代做这几件事从当前种群中通过锦标赛选择选出父代对父代做交叉和变异生成子代种群把父代和子代合并成一个大小为 2N 的种群对合并后的种群做非支配排序计算每个个体的拥挤距离按“非支配层优先、同层拥挤距离大的优先”保留前 N 个个体这一步“合并父代和子代再筛选”就是精英保留机制。最优秀的解不会在进化中丢掉这是 NSGA-II 不容易退化的关键原因。3. 关键代码实现与核心环节拆解3.1 测试问题ZDT1讲代码前先选一个标准的测试函数。ZDT1 是双目标测试问题决策变量维度可以设成 30所有变量范围都在 [0,1] 之间。它的目标是f1(x) x1f2(x) g(x) * (1 - sqrt(x1 / g(x)))其中 g(x) 1 9 * sum(x2 到 xdim) / (dim - 1)ZDT1 的解析 Pareto 前沿是 f2 1 - sqrt(f1)适合用来验证算法写没写对。下面是目标函数代码function fit evaluate_zdt1(pop, dim) N size(pop, 1); fit zeros(N, 2); for i 1:N x pop(i, :); g 1 9 * sum(x(2:dim)) / (dim - 1); fit(i, 1) x(1); fit(i, 2) g * (1 - sqrt(x(1) / g)); end end3.2 非支配排序和拥挤距离的实现快速非支配排序我习惯用“支配计数 支配集合”的方式。对每个个体 i记录有多少个体支配它以及它支配了哪些个体。先找出所有计数为 0 的个体作为第一层然后遍历第一层里的每个个体把它支配集合里那些个体的计数减 1减到 0 就归入下一层。实现如下function [rank, F] non_dominated_sort(fit) N size(fit, 1); S cell(N, 1); n zeros(N, 1); for i 1:N for j 1:N if i j continue; end if dominates(fit(i, :), fit(j, :)) S{i}(end 1) j; elseif dominates(fit(j, :), fit(i, :)) n(i) n(i) 1; end end end rank zeros(N, 1); F {}; front []; for i 1:N if n(i) 0 rank(i) 1; front(end 1) i; end end F{1} front; k 1; while ~isempty(F{k}) Q []; for i F{k} for j S{i} n(j) n(j) - 1; if n(j) 0 rank(j) k 1; Q(end 1) j; end end end k k 1; F{k} Q; end F(end) []; end function isDom dominates(a, b) isDom all(a b) any(a b); end拥挤距离的计算要分目标做。对每一层里的个体按某个目标值排序排序后首尾两个个体的距离直接设为无穷大保证边界解一定能留下来。中间个体的距离用相邻两个个体的目标差除以该目标的最大最小值差再累加所有目标的贡献。function cd crowding_distance(fit, F) n size(fit, 1); cd zeros(n, 1); for i 1:numel(F) inds F{i}; m numel(inds); if m 2 cd(inds) inf; continue; end for obj 1:size(fit, 2) [~, order] sort(fit(inds, obj)); sorted inds(order); cd(sorted(1)) inf; cd(sorted(end)) inf; fmin fit(sorted(1), obj); fmax fit(sorted(end), obj); if fmax fmin continue; end for k 2:m - 1 cd(sorted(k)) cd(sorted(k)) ... (fit(sorted(k 1), obj) - fit(sorted(k - 1), obj)) / (fmax - fmin); end end end end3.3 主循环、交叉和变异主程序里我用了一个简化版锦标赛选择随机取两个个体谁支配谁就选谁如果互不支配就随便取第二个。正式项目中这一步建议改成“先看非支配层层数小的胜出同层再比拥挤距离”。教学版这样写已经能跑但别直接套到论文里。完整主脚本如下clear; clc; N 100; % 种群规模 maxGen 200; % 进化代数 dim 30; % 决策变量维数 lb zeros(1, dim); % 下界 ub ones(1, dim); % 上界 pc 0.9; % 交叉概率 pm 1 / dim; % 变异概率 etaC 15; % 交叉分布指数 etaM 20; % 变异分布指数 pop lb rand(N, dim) .* (ub - lb); fit evaluate_zdt1(pop, dim); for gen 1:maxGen offspring zeros(N, dim); for i 1:N p randi([1 N], 1, 2); if dominates(fit(p(1), :), fit(p(2), :)) p1 pop(p(1), :); else p1 pop(p(2), :); end p randi([1 N], 1, 2); if dominates(fit(p(1), :), fit(p(2), :)) p2 pop(p(1), :); else p2 pop(p(2), :); end if rand pc [offspring(i, :), ~] sbx_crossover(p1, p2, lb, ub, etaC); else offspring(i, :) p1; end offspring(i, :) poly_mutation(offspring(i, :), lb, ub, etaM); end allpop [pop; offspring]; allfit [fit; evaluate_zdt1(offspring, dim)]; [~, F] non_dominated_sort(allfit); cd crowding_distance(allfit, F); order []; for k 1:numel(F) [~, idx] sort(cd(F{k}), descend); order [order, F{k}(idx)]; end pop allpop(order(1:N), :); fit allfit(order(1:N), :); end [~, F] non_dominated_sort(fit); pareto fit(F{1}, :); figure; plot(pareto(:, 1), pareto(:, 2), o); xlabel(f1); ylabel(f2); title(NSGA-II on ZDT1); grid on;交叉算子用模拟二进制交叉 SBX它的思路是用一个随机数生成子代与父代之间的分布因子让子代尽量靠近父代同时保留一定探索能力。变异用多项式变异变异幅度由一个分布指数控制。function [c1, c2] sbx_crossover(p1, p2, lb, ub, etaC) u rand(size(p1)); beta zeros(size(p1)); same abs(p1 - p2) 1e-6; beta(same) 0; idx ~same u 0.5; beta(idx) (2 * u(idx)) .^ (1 / (etaC 1)); idx ~same u 0.5; beta(idx) (1 ./ (2 * (1 - u(idx)))) .^ (1 / (etaC 1)); c1 0.5 * ((1 beta) .* p1 (1 - beta) .* p2); c2 0.5 * ((1 - beta) .* p1 (1 beta) .* p2); c1 min(max(c1, lb), ub); c2 min(max(c2, lb), ub); endfunction c poly_mutation(p, lb, ub, etaM) r rand(size(p)); c p; idx r 0.5; c(idx) p(idx) (ub(idx) - lb(idx)) .* ... ((2 * r(idx)) .^ (1 / (etaM 1)) - 1); idx r 0.5; c(idx) p(idx) (ub(idx) - lb(idx)) .* ... (1 - (2 * (1 - r(idx))) .^ (1 / (etaM 1))); c min(max(c, lb), ub); end运行完主脚本ZDT1 的前沿点应该大致落在 f2 1 - sqrt(f1) 这条曲线上。如果差得很远优先检查非支配排序和拥挤距离两个函数有没有写错。4. 参数设置、常见问题与调试技巧4.1 参数怎么给不是越大越好很多新手一上来就把种群规模和迭代次数拉到很大结果跑半天不出结果。实际上 NSGA-II 的参数有常规范围够用就行。参数常用范围说明种群规模 N50 到 200太小容易早熟太大会让每次迭代都很慢迭代次数 maxGen100 到 500可以先跑 100 看趋势再适当增加交叉概率 pc0.8 到 1.0太高会影响收敛稳定性太低探索不够变异概率 pm1 / dim 左右维度越高变异概率越低交叉分布指数 etaC10 到 20控制子代离父代的远近变异分布指数 etaM20 到 100越大变异幅度越小交叉概率不一定要设成 1给一点概率保留父代个体反而能稳定搜索。变异概率如果设得太大算法很容易退化成随机搜索。4.2 常见问题速查我把自己跑 NSGA-II 时踩过的坑整理成了表格基本都是刚写完代码最容易碰到的问题。现象可能原因排查思路报错“数组索引超出范围”非支配排序返回的前沿层数不够或者 order 长度小于 N检查 non_dominated_sort 返回的 F 是否包含空层打印 numel(F) 和 numel(order)目标值出现 NaN 或 Inf目标函数里有除零、开根号负数检查变量是否越界比如 ZDT1 里 x1 不能为负变异后要 clamp 到边界Pareto 前沿点非常少种群太小或迭代没收敛增加 N 到 100 以上增加 maxGen也可能是拥挤距离没把边界个体设为 inf前沿点集中在一小段拥挤距离计算有误或者两个目标数量级差异太大对目标值做归一化后再算拥挤距离结果和解析前沿差很多交叉、变异参数不合适或选择压力不足先调 etaC 和 etaM再把简化版锦标赛替换成“先分层再比拥挤距离”的选择方式其中“数量级差异”是最容易被忽视的。比如目标一是成本数值在几万量级目标二是效率数值在 0 到 1 之间拥挤距离几乎会被成本目标主导效率目标对多样性的贡献被淹没。这种情况下建议先对每个目标做最小最大归一化再计算拥挤距离。4.3 结果验证和解析解对比ZDT1 的好处是知道理论前沿拿到结果后直接对比hold on; f1 linspace(0, 1, 100); plot(f1, 1 - sqrt(f1), k-, LineWidth, 1.5); legend(NSGA-II, Pareto optimal front);如果数值点能比较均匀地贴住这条曲线说明算法的排序、距离、选择逻辑基本没问题。想观察收敛过程可以在主循环里把每一代第一前沿存下来全部结束后用不同颜色画几代的前沿变化能直观看到解是怎么一步步往外扩的。这个习惯比只盯最终结果有用得多因为你能看到算法是否过早停滞。5. 工程落地中的几个额外建议如果是做工程项目不是跑个测试函数就结束后面还有几件事要做。约束处理可以先从罚函数开始。把违反约束的量加到目标函数里实现最快但罚因子不好调。更稳的做法是在个体比较时加入约束支配两个解都可行比较目标一个可行一个不可行可行解胜出两个都不可行违反程度小的胜出。这个思路和 NSGA-II 本身不冲突把支配判断加一段约束判断就行推荐优先尝试。性能优化方面如果种群规模和迭代次数都很高MATLAB 的循环会变慢。可以先学会用 profile 分析瓶颈。通常最后瓶颈在非支配排序的 O(MN²) 循环上。这个时候再用数组操作或者 parfor 并行代替两层 for 循环不要一开始就追求代码向量化先把逻辑调对。多个目标想可视化目标数超过 3 个就没有直观图形了。这时推荐用 Hypervolume 或者 IGD 这类指标做定量比较配合平行坐标图看每个解在不同目标上的数值分布。我在实际项目里常用的做法是跑 20 次独立实验记录每次的 Hypervolume 均值和方差再和别的算法对比比只贴一张前沿图有说服力得多。我实践下来最大的体会是NSGA-II 的框架不复杂真正影响结果的是非支配排序和拥挤距离这两个细节。边界个体忘记设成 inf、拥挤距离没做归一化、排序后索引弄错都会让结果看起来“好像差不多”但实际前沿质量很差。所以建议拿到任何一份参考代码都先用类似 ZDT1 这种有解析前沿的测试函数验证一遍再替换成自己的目标函数。这样即使后面改了问题也能快速定位是自己问题模型出错了还是算法本身的问题。本文还有配套的精品资源点击获取