简介面向多目标优化问题这份Matlab代码实现了经典的NSGA-II算法适合算法研究者、研究生及工程开发人员理解非支配排序遗传算法的核心机制。代码按功能拆分出拥挤距离计算、精英策略保留、遗传操作、非支配排序、目标函数定义等模块并基于ZDT1-6与DTLZ1-6标准测试函数给出了完整仿真附有测试数据和结果图像便于对照验证Pareto前沿求解效果。整个压缩包共27个文件其中10个m脚本对应各算法环节10个txt存放测试数据7个fig为生成的二维/三维前沿图像整体体积2.41MB结构清晰、便于按模块阅读和二次开发。目前已有3201人学习下载尤其适合初次接触多目标优化、希望快速跑通NSGA-II代码的读者。 写NSGA-II的Matlab代码这件事表面上是“照着论文撸一遍算法”实际上踩到的坑会比你预想的多得多。这篇内容是我自己把非支配排序遗传算法第二代从头实现一遍后的完整记录包含核心机制怎么理解、代码文件怎么组织、每个关键函数怎么实现以及我在调试、调参、扩展时遇到的问题和解决办法。不管你是刚接触多目标优化、想在毕业设计里用NSGA-II跑个对比实验还是想把算法改到像VRPTW这类实际问题上去这份记录都能给你省下不少搜代码、试错的时间。1. 先搞懂NSGA-II的三个支柱代码才有逻辑1.1 为什么多目标问题要Pareto而不是“加权求和”很多初学者一上来就问NSGA-II和普通遗传算法到底差在哪关键就一句话普通遗传算法处理单目标NSGA-II处理多个互相冲突的目标。比如我要买一台笔记本电脑既想要性能强又想要价格低——这两个目标打架没有绝对最优解只有“比某个解更好”的解。这种“别人没法全面打赢你”的解就叫Pareto最优解。NSGA-II要做的就是一次性找出一整组这样的折中方案供你根据偏好挑选。这里的核心机制有三块快速非支配排序Fast Non-Dominated Sort、拥挤度距离Crowding Distance、精英保留选择Elitism。业界常说NSGA-II解决了早期多目标算法的两大问题计算复杂度高和种群多样性差。前者通过O(MN²)的排序算法搞定后者靠拥挤度距离和锦标赛选择维持。后面写代码时会发现这三个机制环环相扣少一个整个算法都跑不动。1.2 主流程和MATLAB代码结构的对应关系NSGA-II的流程其实很固定论文里那张经典流程图看懂了代码结构也就定下来了初始化种群 → 评估目标值 → 非支配排序 → 拥挤度计算 → 锦标赛选择 → 交叉变异生成子代 → 合并父子代 → 环境选择重新排序修剪→ 循环直到达到最大迭代次数。我把这个流程映射到了MATLAB的工程文件上每个函数单独一个文件目录结构如下nsga2_project/ ├── nsga2_main.m % 主程序入口 ├── init_population.m % 初始化种群 ├── evaluate_objective.m % 目标函数改成自己的问题即可 ├── non_dominated_sort.m % 快速非支配排序 ├── crowding_distance.m % 拥挤度距离计算 ├── tournament_select.m % 锦标赛选择 ├── sbx_crossover.m % 模拟二进制交叉 ├── polynomial_mutation.m % 多项式变异 ├── environmental_select.m % 精英保留环境选择 └── plot_pareto.m % 画Pareto前沿这种“一函数一文件”的写法好处是在实验时能单独测试每个环节。比如你怀疑排序写错了直接在命令行调用non_dominated_sort输入几组目标向量看输出就完事不用跑整个算法。我自己最初是把所有代码塞在一个脚本里结果调试一次得跑十几分钟后来拆开才痛快。主程序里建议加一个“固定随机种子”的步骤% 固定随机种子保证实验结果可复现 rng(42);别小看这一行。多目标优化实验写论文时必须多次重复运行不固定种子你连自己都说不清楚结果是随机出来的还是算法真的有效。2. 初始化、参数设置与种群表示2.1 决策变量编码与边界约束NSGA-II默认采用实数编码每个个体是一个决策变量向量。对于连续优化问题这个编码方式很自然。比如ZDT1测试函数有n维决策变量每个变量范围是[0,1]初始化代码可以这样写function pop init_population(nPop, nVar, lb, ub) % nPop: 种群规模 % nVar: 决策变量维数 % lb, ub: 下界和上界向量 pop zeros(nPop, nVar); for i 1:nPop pop(i, :) lb (ub - lb) .* rand(1, nVar); end end如果你的实际问题里决策变量不是连续值比如VRPTW问题中每个客户点访问顺序是离散排列那这里的编码方式就要彻底换掉。很多初学者改不动标准NSGA-II代码多半是在“编码”这一层就卡住了。连续变量 → 实数编码离散顺序 → 排列编码这一点必须先想清楚再动手。2.2 一组能直接用起来的参数范围写NSGA-II之前先给自己定一套初始参数不用追求最优先让代码跑通。我常用的起步参数如下参数推荐取值作用种群规模 nPop100 ~ 200太小易早熟太慢则计算量大迭代次数 nGen200 ~ 500看问题和计算耗时交叉概率 pc0.8 ~ 0.9控制子代产生数量变异概率 pm1 / nVar一般取决策变量数的倒数交叉分布指数 eta_c15 ~ 20SBX交叉的分布程度变异分布指数 eta_m20 ~ 100多项式变异的分布程度这组参数不是拍脑袋来的。交叉概率太低子代多样性不足太高则近似随机搜索。变异概率取1/nVar是为了保证平均每个个体大约有一个决策变量发生变异这是遗传算法里常用的经验比例。eta_c、eta_m越大生成的后代越接近父代搜索更精细但可能陷入局部越小则后代偏离越远探索能力强但收敛慢。所以这两项我一般先取中间值观察收敛曲线后再调整。3. 核心函数实现与逐段代码讲解3.1 快速非支配排序怎么用O(MN²)完成前沿分级非支配排序是NSGA-II的灵魂。解释一下支配关系如果解A在所有目标上都不劣于解B且至少在一个目标上严格优于B那么A支配B。排序要做的是把所有个体划分到不同前沿第一前沿是所有不被任何其他解支配的解第二前沿是去掉第一前沿后剩下的不被支配的解以此类推。我用的是经典的两两比较逻辑很直白function [fronts, rank] non_dominated_sort(objValues) % objValues: nPop x nObj每一行是一个解的目标值 nPop size(objValues, 1); dominatedCount zeros(1, nPop); % 被多少个体支配 dominatedSet cell(1, nPop); % 它支配哪些个体 for i 1:nPop for j 1:nPop if i j continue; end if dominates(objValues(i, :), objValues(j, :)) dominatedSet{i} [dominatedSet{i}, j]; elseif dominates(objValues(j, :), objValues(i, :)) dominatedCount(i) dominatedCount(i) 1; end end end fronts {}; currentFront find(dominatedCount 0); while ~isempty(currentFront) fronts{end1} currentFront; nextFront []; for i currentFront for j dominatedSet{i} dominatedCount(j) dominatedCount(j) - 1; if dominatedCount(j) 0 nextFront [nextFront, j]; end end end currentFront nextFront; end rank zeros(nPop, 1); for k 1:length(fronts) rank(fronts{k}) k; end end function d dominates(x, y) d all(x y) any(x y); end这段代码的核心思想是“计数分层”先统计每个个体被谁支配、支配谁然后从“没人支配它”的个体开始逐层剥离。实际调试时我建议先用一个4个个体的小例子比如objValues [1,2; 2,1; 3,3; 1.5,1.5]手算一遍再来跑代码确认前沿划分正确后再接下游。注意这段代码是教学版当nPop达到几千时双重循环会变得很慢。工程上可以用排序优化但对常规实验这个版本完全够用。向量化技巧后面单独讲。3.2 拥挤度距离拿什么保证解的多样性光有排序还不够。同一前沿内的解有优劣之分吗从Pareto角度看它们都很优秀但我们需要选择一些“更有代表性”的解保留下来。NSGA-II的做法是计算每个解周围的拥挤程度——周围越空旷说明这个区域解越稀疏越值得保留。拥挤度距离的计算方式是对每个目标值排序边界个体最大值、最小值直接给一个无穷大的距离保证它们一定被选中内部个体的距离是相邻两个个体在该目标上归一化差值之和。function dist crowding_distance(objValues) nPop size(objValues, 1); nObj size(objValues, 2); dist zeros(1, nPop); for m 1:nObj [sortedValues, idx] sort(objValues(:, m)); fmin sortedValues(1); fmax sortedValues(end); dist(idx(1)) inf; dist(idx(end)) inf; for i 2:nPop-1 if fmax ~ fmin dist(idx(i)) dist(idx(i)) (sortedValues(i1) - sortedValues(i-1)) / (fmax - fmin); end end end end3.3 锦标赛选择压力和随机性的平衡有了排序等级rank和拥挤度距离dist选择父代就成了一个比较规则优先选择rank小的个体如果rank相同选择距离大的个体。这个规则在代码里体现为function parent tournament_select(pop, rank, dist, k) % pop: 决策变量种群 % k: 锦标赛规模一般取2 nPop size(pop, 1); parent zeros(size(pop, 1), size(pop, 2)); for i 1:nPop candidates randperm(nPop, k); best candidates(1); for j 2:k c candidates(j); if rank(c) rank(best) || (rank(c) rank(best) dist(c) dist(best)) best c; end end parent(i, :) pop(best, :); end end锦标赛选择是“有压力的随机抽样”规模k越大选择压力越大收敛快但更容易早熟。k2是标准设置先用它跑通整体流程。3.4 SBX交叉与多项式变异让种群“生”出新解SBX交叉模拟的是二进制编码中单点交叉的分布效果。给定父代p1、p2按概率生成子代c1、c2。核心是计算beta因子function [c1, c2] sbx_crossover(p1, p2, lb, ub, eta_c) % p1, p2是一个个体的决策变量向量 beta zeros(size(p1)); u rand(size(p1)); beta(u 0.5) (2 * u(u 0.5)).^(1 / (eta_c 1)); beta(u 0.5) (1 ./ (2 * (1 - u(u 0.5)))).^(1 / (eta_c 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); end多项式变异则是在当前解上叠加一个可控的扰动function child polynomial_mutation(x, lb, ub, eta_m) u rand(size(x)); delta zeros(size(x)); idx1 u 0.5; idx2 u 0.5; delta(idx1) (2 * u(idx1)).^(1 / (eta_m 1)) - 1; delta(idx2) 1 - (2 * (1 - u(idx2))).^(1 / (eta_m 1)); child x delta .* (ub - lb); child min(max(child, lb), ub); end注意一点SBX之后必须做变量边界截断。否则在边界附近的父代交叉后子代很容易越界后续计算目标时会报错或者导致结果无意义。这是我早期调试踩过最多的坑。4. 测试、调参和性能优化实录4.1 用ZDT系列测试函数验证算法对错拿到一套新写的优化算法第一步不是改参数而是用标准测试函数验证“它对不对”。ZDT1是最常用的两目标测试函数Pareto前沿是凸的f2 1 - sqrt(f1)非常适合验证function [f1, f2] zdt1(x) n length(x); f1 x(1); g 1 9 * sum(x(2:end)) / (n - 1); f2 g * (1 - sqrt(f1 / g)); end跑完程序后把每一代的最优前沿画出来和理论前沿叠在一起如果差距很大那基本可以断定代码有bug而不是参数问题。我个人的测试顺序是先用ZDT1凸前沿再用ZDT2凹前沿最后用ZDT4带局部最优陷阱的复杂地形。三个全过了代码才敢说“基本正确”。4.2 调参的实战经验先收敛后多样性调参这件事很多新手喜欢一上来就交叉、变异概率全面扫参结果跑了一宿也没搞清楚哪个参数影响大。我的经验是先固定种群规模和迭代次数单独调eta_c和eta_m如果最后一代的Pareto前沿明显偏离理论前沿优先增大eta_c让后代更接近父代收敛更细如果前沿“缺块”——有些区域没有解优先增大eta_m增强探索能力或者增大变异概率如果相同迭代次数下前沿越来越稳定但多样性变差试试增大种群规模。我用一个小建议跑一次实验后不只保存最终前沿还要保存每一代前沿保存到workspace。这样你可以看收敛动画直观判断是在“寻找新区域”还是“卡在某个区域细化”。4.3 向量化与并行评估让实验跑得快一倍NSGA-II的评估阶段是最大性能瓶颈。如果目标函数是仿真模型一次评估要几秒甚至几分钟你不可能等它慢慢跑完。Matlab有个很实用的手段是对多个个体做向量化评估——让目标函数能接收整个种群矩阵一次计算全部个体的目标值。如果目标函数无法向量化退而求其次用parforparfor i 1:nPop objValues(i, :) evaluate_objective(pop(i, :)); end注意parfor循环里使用的函数必须能被所有worker访问建议把问题定义函数写成一个独立的.m文件。提示工具箱方面Matlab自带gamultiobj也能求解多目标问题自己写NSGA-II的意义在于完全可定制比如改编码、加约束、自定义交叉算子。但如果你只是想快速得到一组结果先用gamultiobj对比一下结果再回来验证自己的实现也是一种高效路线。5. 常见问题与排查技巧实录5.1 错误排查速查表现象可能原因排查方法运行时提示“索引超出数组范围”非支配排序返回的前沿数量为0在non-dominated-sort后打印rank检查支配比较逻辑最后一轮前沿全是同一个点时变异概率过小、种群早熟增大pm或增大eta_m目标函数返回NaN或Inf决策变量越界或目标函数本身有问题检查交叉变异后是否做了边界截断运行很慢、每代都要十几秒评估目标函数没有向量化改用parfor或重构evaluate_objective多次运行结果差异巨大没有固定随机种子在main代码开头加rng(42)5.2 三个值得注意的工程细节第一个细节种群在matlab里的存储格式尽量用nPop × nVar的矩阵而不要用struct数组。矩阵运算可以直接用向量化struct数组取字段再操作很麻烦。只有当每个个体带约束值、ID等额外信息时才考虑struct。第二个细节环境选择合并父代与子代、再排序、再截断时最后一步“按front和拥挤度填满种群”要小心。排名靠前的前沿全要最后一个前沿只需填补剩余名额需要用拥挤度排序后取前面一部分而不是把这个前沿全塞进去。第三个细节调试时画图非常关键。我建议在main循环里加一个条件语句——每隔20代画一次当前最优前沿最终结果保存成动画格式。这能直观看到算法从“散点”到“收敛成线”的全过程一旦发现异常能立刻定位问题发生的大致迭代轮次。6. 从测试函数到实际问题VRPTW与更多扩展思路6.1 怎么把标准NSGA-II改造成解决VRPTW或调度问题很多人拿着标准测试函数版的NSGA-II想直接处理带时间窗的车辆路径问题VRPTW结果发现交叉、变异完全不起作用。原因很简单VRPTW的决策变量是一组车辆访问客户的排列序列标准实数编码的SBX在这里没有意义。需要做的改动是把每个个体编码成一串客户序列和车辆分配标志将SBX替换为顺序交叉Order CrossoverOX或部分映射交叉PMX将多项式变异替换为交换变异、插入变异或反转变异把时间窗超时量作为约束用罚函数加到目标中。这套改造思路本质上是“问题变了算子和表示跟着变”。理解这一点比自己硬套标准代码重要得多。改造我建议分三步走先不管约束跑通两个目标的优化再加上时间窗约束最后引入车辆数目标做成真正的“优化车辆数总路径惩罚项”的多目标。6.2 我自己的几个使用心得和避坑经验个人体会最深的一件事NSGA-II不是“调个参就能一劳永逸”的算法它对目标函数的尺度很敏感。如果两个目标值量级差太大比如一个在[0,1]区间一个在[1000,100000]区间归一化处理一定不能省否则拥挤度距离基本被大尺度目标主导小目标等于白算。这种情况在ZDT测试函数上看不出来到了实际工程一定会爆。另一个习惯是每跑完一组实验先把本次的种群、参数、随机种子存成.mat文件。后面整理数据、画图、复盘时这组变量就是你的原始证据。代码始终只是算法思想的刻板记录真正有价值的是你能解释清楚每一步为什么这么做。把上面这些函数跑通之后建议你换个测试函数、改改约束条件重新练一遍到那个时候才算真正吃透了NSGA-II。本文还有配套的精品资源点击获取