
简介这份资源是NSGA-II多目标优化算法的MATLAB实现代码面向正在学习进化算法、准备课程设计或科研复现的高年级本科生与研究生。代码围绕非支配排序、拥挤距离计算、精英策略保留、基因操作与函数值评估等核心模块展开并针对二维测试函数ZDT1至ZDT6、三维测试函数DTLZ1至DTLZ6提供了完整测试数据与仿真图像便于读者对照验证算法收敛性与分布性。压缩包共27个文件以m脚本、txt测试数据与fig仿真图为主整体约2.41MB目录结构清晰主程序与各功能模块分离方便按需调用与二次修改。目前已有3204人学习下载适合作为理解NSGA-II流程、调试参数与撰写实验报告的参考范例读者可借此快速搭建多目标优化实验环境掌握非支配排序与拥挤距离的代码实现思路。1. 从一堆互相打架的目标说起NSGA-II 到底在算什么如果你做过参数调优、结构选型或者资源分配大概率遇到过这种局面三个指标互相拉扯改好一个另外两个就崩。比如某图像处理 Demo 里既要识别精度高又要推理耗时低还想让模型体积小——你手动调了一周最后发现只是在三个目标之间来回妥协根本找不到“最优解”。NSGA-II非支配排序遗传算法第二代就是干这个的它不给你一个解而是给你一族互不支配的解让你自己按场景挑。我第一次在 MATLAB 里跑通它的时候最直观的感受是——原来“没有最好只有最合适”这句话在算法层面是可以被量化的。这篇笔记就围绕 NSGA-II 的 MATLAB 代码把种群初始化、非支配排序、拥挤度计算、遗传算子到结果可视化整条链路拆开让你能自己改目标函数、调参数、看懂输出。2. 把 NSGA-II 拆成四块可写的 MATLAB 函数2.1 为什么选 MATLAB 而不是 Python从矩阵化非支配排序说起NSGA-II 的核心计算量集中在非支配排序和拥挤度距离。MATLAB 的矩阵操作在这两块上有天然优势——非支配排序本质上是两两比较用bsxfun或者隐式扩展可以一次性算出整个种群的支配关系矩阵不用写双重循环。我试过用 Python 的 NumPy 做同样的事代码量差不多但 MATLAB 的pdist2和逻辑索引在调试时更直观尤其是当你想在排序过程中临时查看某个个体的支配计数时直接在工作区点开变量就行。另一个现实原因是很多工程团队的目标函数本身就是 Simulink 模型或者 MATLAB 脚本比如某实验室做控制器参数整定目标函数是仿真出来的超调量和调节时间。这种情况下用 MATLAB 写 NSGA-II 可以避免跨语言调用的开销。当然如果你的目标函数是深度学习训练那还是 Python 生态更顺但算法逻辑完全一样移植过去就行。选 MATLAB 的代价是并行化不如 Python 灵活parfor在遗传算法里加速比有限因为每代之间必须同步。不过对于种群规模 100、迭代 200 代以内的场景单机跑几分钟到十几分钟完全在可接受范围。2.2 主循环骨架种群初始化、非支配排序、拥挤度、选择交叉变异先给一个能跑的最小骨架。假设目标函数是evaluate_objectives输入是决策变量矩阵输出是目标值矩阵每行一个个体每列一个目标。下面这段代码是我常用的主循环结构去掉了花哨的终止条件只保留最核心的流程。% 参数设置 pop_size 100; % 种群规模 max_gen 200; % 最大迭代代数 n_var 10; % 决策变量个数 var_min -5; % 变量下界 var_max 5; % 变量上界 n_obj 3; % 目标个数 % 初始化种群 pop repmat(var_min, pop_size, n_var) ... rand(pop_size, n_var) .* repmat(var_max - var_min, pop_size, n_var); obj evaluate_objectives(pop); % 计算目标值 for gen 1:max_gen % 非支配排序与拥挤度计算 [fronts, crowd_dist] non_dominated_sorting(obj); % 二元锦标赛选择 parent_idx tournament_selection(fronts, crowd_dist, pop_size); parents pop(parent_idx, :); % 模拟二进制交叉SBX与多项式变异 offspring sbx_crossover(parents, var_min, var_max, 1.0, 20); offspring polynomial_mutation(offspring, var_min, var_max, 1/n_var, 20); % 合并父代与子代 combined_pop [pop; offspring]; combined_obj [obj; evaluate_objectives(offspring)]; % 环境选择保留前 pop_size 个个体 [pop, obj] environmental_selection(combined_pop, combined_obj, pop_size); % 每 20 代输出一次进度 if mod(gen, 20) 0 fprintf(Generation %d completed\n, gen); end end这段代码的逻辑说明non_dominated_sorting返回每个个体所属的 Pareto 前沿层级和拥挤度距离tournament_selection根据层级优先、拥挤度次优的规则选父代sbx_crossover和polynomial_mutation是 NSGA-II 默认的遗传算子分布指数分别控制交叉和变异的“集中程度”。参数方面交叉概率通常取 0.9 到 1.0变异概率取1/n_var左右分布指数取 20 是文献里的常用值但如果你希望搜索更分散可以降到 10 到 15。environmental_selection是 NSGA-II 区别于普通遗传算法的关键它先把父代和子代合并然后按 Pareto 层级从低到高填充最后一层用拥挤度距离截断。这样既保证了精英保留又维持了种群多样性。2.3 非支配排序的 MATLAB 实现用支配矩阵替代双重循环非支配排序的经典写法是双重循环但 MATLAB 里可以用矩阵运算加速。核心思路是对种群中任意两个个体 i 和 j如果 i 在所有目标上都不差于 j且至少有一个目标严格优于 j则 i 支配 j。用矩阵表示就是function [fronts, crowd_dist] non_dominated_sorting(obj) [n, m] size(obj); dominated false(n, n); % dominated(i,j)true 表示 i 支配 j for i 1:m % 对每个目标比较所有个体 less obj(:, i) obj(:, i); % 小于矩阵 leq obj(:, i) obj(:, i); % 小于等于矩阵 dominated dominated | (less leq); end % 注意上面写法有误正确逻辑需要同时满足所有目标 % 更稳妥的写法是逐对比较但用向量化可以这样 dominated true(n, n); for i 1:m dominated dominated (obj(:, i) obj(:, i)); end strictly_better false(n, n); for i 1:m strictly_better strictly_better | (obj(:, i) obj(:, i)); end dominated dominated strictly_better; dominated(1:n1:end) false; % 自己不算支配自己 % 计算每个个体被支配的次数 n_dominated sum(dominated, 1); % 列求和得到每个个体被多少人支配 % 注意dominated(i,j)true 表示 i 支配 j所以 j 被 i 支配 % 因此 n_dominated(j) 是 j 被支配的次数 % 分层 fronts {}; remaining true(1, n); while any(remaining) current_front find(n_dominated 0 remaining); fronts{end1} current_front; remaining(current_front) false; for i current_front % 找到被 i 支配的个体减少它们的被支配计数 dominated_by_i find(dominated(i, :)); n_dominated(dominated_by_i) n_dominated(dominated_by_i) - 1; end end % 拥挤度距离计算略见下一节 crowd_dist calculate_crowding_distance(obj, fronts); end上面代码里我故意保留了一段“错误写法”的注释因为这是很多人第一次写非支配排序时容易踩的坑不能简单地把每个目标的比较结果做逻辑或必须同时满足“所有目标不差”和“至少一个目标严格好”。正确的做法是先构造dominated矩阵表示“所有目标小于等于”再构造strictly_better表示“至少一个目标严格小于”两者相与。参数说明obj是 n×m 矩阵n 是个体数m 是目标数。n_dominated初始为每个个体被支配的次数第一层就是次数为 0 的个体。每移除一层就把该层个体所支配的个体的计数减一直到所有个体被分层。这个实现的时间复杂度是 O(m·n²)对于 n100、m3 的情况矩阵运算瞬间完成。如果 n 上千可以考虑用快速非支配排序的改进版本但一般工程问题 n 在几百以内足够。2.4 拥挤度距离让解在 Pareto 前沿上均匀铺开拥挤度距离的作用是衡量一个个体周围解的密集程度。距离越大说明周围越稀疏越应该保留。计算方法是对每个目标按该目标值排序边界个体的拥挤度设为无穷大中间个体的拥挤度等于相邻个体在该目标上的差值归一化后累加。function crowd_dist calculate_crowding_distance(obj, fronts) n size(obj, 1); crowd_dist zeros(n, 1); for f 1:length(fronts) idx fronts{f}; if length(idx) 2 crowd_dist(idx) Inf; continue; end front_obj obj(idx, :); [m, ~] size(front_obj); dist zeros(length(idx), 1); for i 1:m [sorted_vals, sort_order] sort(front_obj(:, i)); dist(sort_order(1)) Inf; dist(sort_order(end)) Inf; if sorted_vals(end) - sorted_vals(1) 0 continue; % 该目标上所有值相同跳过 end for j 2:length(idx)-1 dist(sort_order(j)) dist(sort_order(j)) ... (sorted_vals(j1) - sorted_vals(j-1)) / ... (sorted_vals(end) - sorted_vals(1)); end end crowd_dist(idx) dist; end end逻辑说明对每个前沿逐个目标排序边界点拥挤度设为Inf保证一定被选中。中间点的累加项是相邻两点在该目标上的差值除以该目标的范围这样不同量纲的目标可以公平比较。如果某个目标上所有值相同分母为零直接跳过该目标避免除零错误。参数注意拥挤度距离只在同一前沿内比较才有意义不同前沿之间先比层级层级相同才比拥挤度。这也是tournament_selection里的比较规则。3. 目标函数怎么写从数学表达式到 MATLAB 向量化3.1 三个经典测试函数ZDT1、ZDT2、DTLZ2 的 MATLAB 写法调试 NSGA-II 时别急着上真实工程问题先用标准测试函数验证代码正确性。ZDT1 和 ZDT2 是两目标DTLZ2 可以扩展到三目标。下面给出 ZDT1 的写法决策变量维度通常取 30。function obj evaluate_objectives(x) % ZDT1: 两目标x 为 n×30 矩阵每行一个个体 n size(x, 2); f1 x(:, 1); g 1 9 * mean(x(:, 2:end), 2); h 1 - sqrt(f1 ./ g); f2 g .* h; obj [f1, f2]; endZDT1 的 Pareto 前沿是 f2 1 - sqrt(f1)f1 在 [0,1] 之间。跑完 NSGA-II 后把得到的解画在目标空间应该紧贴这条曲线。如果发现解偏离很远先检查非支配排序是否正确再检查交叉变异算子是否产生了可行解。DTLZ2 的三目标版本function obj evaluate_objectives(x) % DTLZ2: 三目标x 为 n×12 矩阵 n size(x, 2); k n - 2; % 对于三目标k n - (m-1) n - 2 g sum((x(:, 3:end) - 0.5).^2, 2); f1 (1 g) .* prod(cos(x(:, 1:2) * pi / 2), 2); f2 (1 g) .* cos(x(:, 1) * pi / 2) .* sin(x(:, 2) * pi / 2); f3 (1 g) .* sin(x(:, 1) * pi / 2); obj [f1, f2, f3]; endDTLZ2 的 Pareto 前沿是单位球面第一象限部分。三目标可视化用三维散点图能直观看到解是否均匀分布在球面上。3.2 把工程约束塞进目标函数罚函数法的参数怎么定真实问题往往带约束比如某资源分配问题要求总成本不超过预算。NSGA-II 本身不直接处理约束常见做法是罚函数法把约束违反量乘以罚系数加到目标上。但罚系数太大会导致种群迅速挤到可行域边界太小又压不住不可行解。我一般用自适应罚系数初始罚系数设为 1每代根据可行解比例调整。如果可行解少于 20%罚系数乘以 1.5如果可行解多于 80%罚系数乘以 0.9。这样能在探索和可行之间找平衡。function obj evaluate_objectives_with_constraint(x, budget) % 假设两个目标成本 f1收益 f2负值表示最大化 cost sum(x.^2, 2); profit -sum(sqrt(x), 2); violation max(0, cost - budget); penalty 1e3 * violation; % 初始罚系数 obj [cost penalty, profit penalty]; end参数说明budget是约束上限violation是超出量。罚系数 1e3 是经验值实际要根据目标量级调整。如果成本本身在 1e6 量级罚系数也要相应放大。更稳妥的做法是先跑几代观察可行解比例再手动调一次。3.3 向量化评估一次算完整个种群而不是循环MATLAB 里最拖慢速度的就是在目标函数里写for循环逐个个体计算。上面的 ZDT1 和 DTLZ2 都是向量化写法输入整个种群矩阵输出目标矩阵。如果你的目标函数依赖外部仿真没法向量化那就用arrayfun或者parfor并行。但注意parfor在遗传算法里加速有限因为每代都要同步而且启动并行池本身有开销。对于种群 100、迭代 200 的情况向量化后单代耗时通常在 0.01 秒以内总时间几秒钟没必要并行。一个常见的翻车场景目标函数里用了rand或者tic/toc导致每次评估结果不一致算法无法收敛。记住目标函数必须是确定性的随机性只能出现在遗传算子中。4. 避坑与排查NSGA-II 在 MATLAB 里最容易翻车的五个地方4.1 现象种群早熟所有解挤在一起原因变异率太低或选择压力过大早熟是遗传算法的通病。NSGA-II 里如果变异概率设成 0.001而决策变量有 30 个那每代期望变异次数不到 0.03种群很快失去多样性。解决方法是把变异概率设为1/n_var即每个变量平均每代变异一次。另外检查锦标赛选择的规模如果每次选 10 个比 1 个选择压力太大改成二元锦标赛选 2 个比 1 个。4.2 现象Pareto 前沿断断续续某些区域没有解原因拥挤度计算时边界处理不当如果拥挤度距离没有把边界点设为Inf边界点可能因为距离小而被淘汰导致前沿两端缺失。检查calculate_crowding_distance里是否对排序后的第一个和最后一个个体赋了Inf。另外如果某个目标上所有值相同分母为零代码里要跳过该目标否则会产生NaN污染整个距离向量。4.3 现象算法跑完得到的解还不如手动调参原因目标函数量纲差异过大三个目标里一个在 0.001 量级一个在 1000 量级非支配排序会偏向大量级的目标小量级目标的变化被淹没。解决办法是先归一化目标值或者在拥挤度计算时用归一化后的差值。我一般会在evaluate_objectives返回前做一次 min-max 归一化但注意归一化参数要固定不能每代变否则目标函数就不确定了。4.4 现象迭代到后面几代目标值不再改善原因环境选择时精英保留策略写错NSGA-II 的环境选择是合并父代和子代后选前 N 个。如果写成只从子代里选或者合并后没有按层级和拥挤度排序就会丢失精英。检查environmental_selection函数先按层级从低到高取最后一层按拥挤度从大到小取直到凑够pop_size。一个容易忽略的细节是最后一层截断时被截断的个体不能简单丢弃而是要用拥挤度排序后取前几个。4.5 现象MATLAB 报错“矩阵维度不一致”原因交叉变异后变量越界未处理SBX 交叉和多项式变异可能产生超出[var_min, var_max]的解。如果不做边界处理目标函数可能报错或者返回Inf。常见做法是越界后随机重置到边界内或者用反射法折回。我一般用简单截断offspring min(max(offspring, var_min), var_max)。但截断会损失多样性更好的做法是越界时重新生成该变量。5. 进阶技巧用 MATLAB 的并行池和结果可视化把效率拉满5.1 并行评估目标函数parfor 的正确打开方式当目标函数是外部仿真比如 Simulink 模型时单次评估可能几秒钟种群 100 就是几分钟。这时候parfor能派上用场。但注意parfor要求循环体独立且不能有嵌套的parfor。写法如下function obj evaluate_objectives_parallel(pop) n size(pop, 1); obj zeros(n, 2); parfor i 1:n obj(i, :) simulate_model(pop(i, :)); % 外部仿真函数 end end启动并行池用parpool(local, 4)4 是核心数。注意parfor第一次运行会花几秒启动池如果单次评估很快比如 0.01 秒并行反而更慢。所以只对耗时超过 0.1 秒的目标函数用并行。5.2 结果可视化二维和三维 Pareto 前沿的绘制脚本跑完算法后把obj里的非支配解画出来。二维用scatter三维用scatter3并加上真实 Pareto 前沿作为参考。% 假设 obj 是最终种群的目标值fronts 是分层结果 front1 obj(fronts{1}, :); figure; if size(obj, 2) 2 scatter(front1(:,1), front1(:,2), 40, filled); hold on; f1 linspace(0, 1, 100); plot(f1, 1 - sqrt(f1), r--, LineWidth, 1.5); % ZDT1 真实前沿 xlabel(f1); ylabel(f2); legend(NSGA-II, True Pareto Front); elseif size(obj, 2) 3 scatter3(front1(:,1), front1(:,2), front1(:,3), 40, filled); xlabel(f1); ylabel(f2); zlabel(f3); grid on; end参数说明fronts{1}是第一层非支配解通常就是最终输出的 Pareto 前沿。如果第一层解太多可以只画前 100 个。三维图用scatter3后旋转视角能直观看到解在球面上的分布。5.3 一个我常用的收敛性检查习惯每次跑完我会把每代第一层解的数量记录下来画一条曲线。如果数量从几十降到个位数说明种群多样性丢失需要调大变异率或减小选择压力。如果数量一直波动但目标值不改善可能是目标函数本身有问题。这个习惯帮我省了很多次重新跑的时间。另外我习惯把最终种群和fronts保存成.mat文件方便下次直接加载分析不用重新跑算法。文件名带上日期和参数比如nsga2_zdt1_pop100_gen200.mat避免版本混乱。希望帮到你。本文还有配套的精品资源点击获取