简介基于鲸鱼算法NSWOA的多目标优化求解Matlab代码包面向智能优化算法、神经网络预测、路径规划等方向的本科生、研究生及科研人员适合算法对比、课程实验与论文复现等场景。代码包以NSWOA为核心完整集成非支配排序、变量初始化、目标函数评价、种群更新、替换染色体和结果绘图等子模块并内置zdt1、zdt2、zdt3等多组常用测试函数运行后可直接查看多目标Pareto前沿的收敛效果直观展示算法在多个目标间的权衡能力有助于理解多目标进化算法的迭代逻辑。压缩包共13个文件其中11个为.m源文件便于逐行阅读和修改另附1个说明文档和1个示例效果图整体仅24KB结构精简、核心算法集中能快速定位到主程序与关键子函数。以上代码同时支持Matlab 2014及2019a环境且描述中注明内含运行结果适合学习时对照验证。目前已有259人学习下载尤其适合刚接触多目标优化、希望从代码层面入手掌握鲸鱼算法原理的教研人员。1. 从单目标到多目标鲸鱼算法为什么要改造成 NSWOA做优化的人都有个直觉单目标问题解出来是个点多目标问题解出来是条线、一个面甚至一片不规则的云。鲸鱼算法Whale Optimization Algorithm, WOA是 2016 年提出的群体智能算法模拟座头鲸的泡泡网捕食策略在单目标连续优化上表现亮眼收敛快、参数少、实现简单。但把 WOA 直接扔到多目标问题上会立刻暴露短板它没有维护种群多样性的机制也没有判断一个解是否被另一个解支配的框架最后所有鲸鱼会朝一个方向挤得到一堆彼此雷同的解。NSWOA 的全称是 Non-dominated Sorting Whale Optimization Algorithm思路很直白把 NSGA-II 里那套非支配排序和拥挤度距离搬进 WOA 的寻优循环里让鲸鱼在搜索过程中既保留 WOA 的位置更新机制又用 Pareto 支配关系筛选精英解。这样得到的不是唯一最优解而是一组分布均匀、覆盖前沿的折中方案。工程上的价值在于物流配送路径、车间调度、资源分配这类问题往往有多个互相冲突的目标运维工程师、算法工程师和科研人员需要在一组 Pareto 解里根据业务偏好做决策而不是被迫接受一个单点答案。这篇文章就用 MATLAB 把 NSWOA 的完整实现拆开讲从核心算子到参数设置再到和 NSGA-II 的对比一步步写清楚。2. 核心机制拆解NSWOA 如何在迭代中同时顾全收敛性与分布性2.1 非支配排序给每个解一个「等级」而不是一个「分数」多目标优化和单目标最大的区别在于比较方式。单目标可以用数值大小直接排序多目标需要先判断支配关系。一个解 x1 支配 x2当且仅当 x1 在所有目标上都不劣于 x2且至少在一个目标上严格优于 x2。不满足这个条件时两个解互不支配它们属于同一个 Pareto 前沿。NSWOA 在每次迭代后都要对当前种群做一遍非支配排序把解分成 F1、F2、F3 多个层级F1 是当前最优前沿。实现上常见做法是用一个双层循环记录每个解被谁支配、支配了多少个解然后逐层提取前沿。MATLAB 里写起来大约是这样function front non_dominated_sort(fitness) [N, M] size(fitness); % N个解M个目标 dominate_count zeros(N, 1); % 记录支配当前解的解数量 dominated_set cell(N, 1); % 记录当前解支配的解集合 front zeros(N, 1); % 每个解的Pareto层级 for i 1:N for j i1:N if dominates(fitness(i,:), fitness(j,:)) dominated_set{i} [dominated_set{i}, j]; dominate_count(j) dominate_count(j) 1; elseif dominates(fitness(j,:), fitness(i,:)) dominated_set{j} [dominated_set{j}, i]; dominate_count(i) dominate_count(i) 1; end end end current_front find(dominate_count 0); % 第一前沿 rank 1; while ~isempty(current_front) for idx current_front front(idx) rank; end next_front []; for idx current_front for dominated_idx dominated_set{idx} dominate_count(dominated_idx) dominate_count(dominated_idx) - 1; if dominate_count(dominated_idx) 0 next_front [next_front, dominated_idx]; end end end rank rank 1; current_front next_front; end end这段代码的核心逻辑是先统计每个解被多少解支配被零个解支配的划定到 F1然后把 F1 的解从支配关系中「剔除」更新计数再找下一批零支配解重复直到所有解都有等级。front数组最终给出每个解所在的层级。注意 MATLAB 的下标从 1 开始代码里所有索引都对应这个习惯直接用 0 会报错。2.2 拥挤度距离在同一前沿里挑「住得最宽敞」的鲸鱼层级只解决了解的优劣问题没解决分布问题。如果筛选精英时只看层级同一前沿内随机选解很容易堆在一起Pareto 前沿中间缺一大块边缘却挤满了点。拥挤度距离的思路是对同一前沿的解按每个目标分别排序计算每个解与相邻两个解在目标空间的距离之和距离越大说明周围越空旷越值得保留。function distance crowding_distance(fitness, front) [N, M] size(fitness); distance zeros(N, 1); for level 1:max(front) idx find(front level); if length(idx) 2 distance(idx) inf; % 少于3个解时都设为无穷大全部保留 continue; end for m 1:M [~, order] sort(fitness(idx, m)); sorted_idx idx(order); distance(sorted_idx(1)) inf; distance(sorted_idx(end)) inf; f_min fitness(sorted_idx(1), m); f_max fitness(sorted_idx(end), m); if f_max - f_min 1e-12 continue; % 目标值全相同距离记为0 end for k 2:length(sorted_idx)-1 distance(sorted_idx(k)) distance(sorted_idx(k)) ... (fitness(sorted_idx(k1), m) - fitness(sorted_idx(k-1), m)) / (f_max - f_min); end end end end边界解的距离被直接设为inf这个细节很关键。边界解承担着扩展前沿的任务即使中间区域密度很大也必须把两端保住否则下一次迭代前沿就向内收缩了。距离计算时要按目标分别归一化不然量纲大的目标会主导距离值量纲小的目标在筛选中几乎不起作用。如果某个目标的所有解值都相同分母趋近于零代码里用continue跳过把这个目标的贡献记为零避免产生 NaN。2.3 鲸鱼位置更新机制包围、泡泡网与随机搜索三选一NSWOA 保留了 WOA 的位置更新策略每次迭代每头鲸鱼以一定概率选择三种行为之一。第一种是包围猎物位置向当前最优解靠拢第二种是泡泡网攻击用螺旋方程模拟座头鲸吐气泡向上盘旋的动作第三种是随机搜索向种群中随机一个个体移动维持全局探索能力。function new_positions whale_update(positions, best_pos, a, b, l, p, r1, r2) [N, dim] size(positions); new_positions zeros(N, dim); for i 1:N if p(i) 0.5 if abs(a(i)) 1 % 收缩包围向当前最优解靠近 new_positions(i,:) best_pos - a(i) * abs(b(i) * best_pos - positions(i,:)); else % 随机搜索选择一个随机个体作为参照 rand_idx randi(N); rand_whale positions(rand_idx,:); new_positions(i,:) rand_whale - a(i) * abs(b(i) * rand_whale - positions(i,:)); end else % 螺旋更新以螺旋路径逼近最优解 dist abs(best_pos - positions(i,:)); new_positions(i,:) dist .* exp(b * l(i)) .* cos(2 * pi * l(i)) best_pos; end end end参数a随迭代从 2 线性递减到 0控制了探索和开发的平衡。abs(a) 1时鲸鱼向最优解收缩abs(a) 1时跳出局部陷阱做全局随机搜索。注意这里用的是randi(N)随机选参照个体和单目标 WOA 完全一致。问题在于多目标场景下「最优解」不再唯一——F1 前沿里每个解都是最优选哪个作为best_pos直接影响搜索方向。NSWOA 的常见做法是在 F1 前沿里随机挑一个每次迭代不同鲸鱼可能朝向不同的 Pareto 解移动这样种群不至于全部塌缩到同一个点。我在实际写代码时会在每次迭代开始时先对 F1 前沿做一次洗牌然后每头鲸鱼按自己的序号取模选参照比每次都重新调用randi更稳定。3. 在 MATLAB 中搭建 NSWOA 求解多目标问题的完整流程3.1 外部存档与精英保留NSWOA 的归档策略怎么写NSWOA 和很多多目标元启发式算法的区别在于是否有外部存档。NSGA-II 是靠精英选择机制把父代和子代合并后筛选NSWOA 更常见的做法是维护一个独立的存档集用非支配排序加拥挤度距离控制存档规模。每次迭代结束后把当前种群和存档合并做一次非支配排序取前ArchiveSize个解作为新存档。我做这个部分的代码时发现最值得注意的坑是越界处理。多目标问题里变量的边界约束很常见但进化算法更新完位置后经常有解跑出边界。直接把这些解扔掉会减少种群多样性简单粗暴地「截断」到边界又会让大量解堆在边界线上影响拥挤度计算的真实性。我一般用随机反射策略超出边界的维度按边界宽度折返类似光在镜面反射这样既保留了解的位置又避免边界堆积。function x boundary_reflect(x, lb, ub) range ub - lb; over_high x ub; over_low x lb; x(over_high) ub(over_high) - mod(x(over_high) - ub(over_high), range(over_high)); x(over_low) lb(over_low) mod(lb(over_low) - x(over_low), range(over_low)); idx (x lb) | (x ub); % 极端情况下取模后仍然越界 x(idx) lb(idx) rand(size(lb(idx))) .* range(idx); if any(x(:) lb(:)) || any(x(:) ub(:)) x max(min(x, ub), lb); % 兜底保证在边界内 end end这段代码里的逻辑顺序很重要先做反射再对极端情况做随机初始化最后兜底截断。如果一上来就截断反射策略就失去了意义。实际跑的时候发现反射后仍然越界的概率很低但一旦发生通常是变量维度特别多、搜索步长特别大的情况兜底处理是必要的安全网。3.2 种子函数与主循环变量数、目标数、代数怎么往代码里填用户拿到的 NSWOA 代码压缩包里通常包含NSWOA.m主函数、目标函数文件和一些绘图脚本。我建议在跑任何测试问题之前先搞清楚主函数的输入输出接口。典型的主循环骨架长这样function [archive, archive_fitness, history] NSWOA(fun, dim, lb, ub, N, MaxIter, ArchiveSize) positions lb rand(N, dim) .* (ub - lb); fitness evaluate_population(fun, positions); archive positions; archive_fitness fitness; history zeros(MaxIter, 1); for t 1:MaxIter a 2 - 2 * t / MaxIter; % 线性递减参数 b 1; % 螺旋形状常数 best_idx select_leader(archive, archive_fitness); % 从存档里选领导者 p rand(N, 1); r1 rand(N, dim); r2 rand(N, dim); l -1 2 * rand(N, 1); A 2 * a .* r1 - a; C 2 .* r2; new_positions zeros(N, dim); for i 1:N if p(i) 0.5 if abs(A(i,1)) 1 new_positions(i,:) archive(best_idx,:) - A(i,:) .* abs(C(i,:) .* archive(best_idx,:) - positions(i,:)); else rand_idx randi(N); new_positions(i,:) positions(rand_idx,:) - A(i,:) .* abs(C(i,:) .* positions(rand_idx,:) - positions(i,:)); end else dist abs(archive(best_idx,:) - positions(i,:)); new_positions(i,:) dist .* exp(b * l(i)) .* cos(2 * pi * l(i)) archive(best_idx,:); end end new_positions boundary_reflect(new_positions, lb, ub); new_fitness evaluate_population(fun, new_positions); % 合并种群与存档非支配排序裁剪 combined_pos [positions; new_positions; archive]; combined_fit [fitness; new_fitness; archive_fitness]; front non_dominated_sort(combined_fit); dist crowding_distance(combined_fit, front); % 按层级和拥挤度排序取前N个 [~, sort_idx] sortrows([front, -dist]); positions combined_pos(sort_idx(1:N), :); fitness combined_fit(sort_idx(1:N), :); % 从排序后的解里提取存档 archive positions(front(sort_idx(1:N)) 1, :); archive_fitness fitness(front(sort_idx(1:N)) 1, :); if size(archive, 1) ArchiveSize archive trim_archive(archive, archive_fitness, ArchiveSize); archive_fitness evaluate_population(fun, archive); end history(t) sum(min(archive_fitness, [], 1)); % 记录目标值的动态 end end这个主循环里有个和纯 NSGA-II 不同的设计位置更新时参照的不是「本代的最优解」而是「存档里的最优解」。存档里保留的是历史迭代中发现的全部非支配解比当前代的 F1 前沿更全面。select_leader的实现通常是做一次拥挤度距离计算取距离最大的解作为 leader这样能引导种群朝前沿的稀疏区域搜索而不是扎堆在已经密集的区域。主函数的参数设置直接影响求解效果。下面这张表是经验值适合 ZDT 系列和 DTLZ 系列测试函数实际工程问题需要按变量维度和目标数量调整参数推荐范围说明种群规模 N50~200目标函数计算代价高时取小值最大迭代 MaxIter200~1000变量维度高时取大值存档规模 ArchiveSizeN 的 50%~100%过小会丢失前沿细节螺旋常数 b0.5~1.2控制螺旋的紧密程度目标数 M2~3超过 3 个目标建议用 NSGA-III 的参考点机制实际运行中如果发现最后得到的 Pareto 前沿不够连续优先增大ArchiveSize而不是MaxIter。存档空间太小即使迭代再多留下的非支配解也被迫裁剪前沿必然稀疏。3.3 从压缩包到可运行NSWOA 代码的文件结构与调用方式拿到 NSWOA 的压缩包解压后通常能看到这些文件NSWOA.m是主算法ZDT1.m或类似的测试函数文件PlotPareto.m是画图脚本有时还有一个main.m把所有步骤串起来。先用测试函数验证算法正确性再替换成自己的目标函数这是最稳的上手顺序。% main.m 示例 clc; clear; close all; % 定义ZDT1测试函数两个目标30个变量 dim 30; lb zeros(1, dim); ub ones(1, dim); N 100; MaxIter 500; ArchiveSize 100; [archive, archive_fitness, history] NSWOA(ZDT1, dim, lb, ub, N, MaxIter, ArchiveSize); figure; plot(archive_fitness(:,1), archive_fitness(:,2), o); xlabel(f1); ylabel(f2); title(NSWOA Pareto Front - ZDT1); % 计算IGD指标评估质量 true_front load(ZDT1_true_front.txt); igd_value IGD(archive_fitness, true_front); fprintf(IGD %.4f\n, igd_value);调用时注意两点目标函数必须返回一个列向量每个元素是一个目标值NSWOA 内部统一按行组织种群、按列组织目标lb和ub若是全部变量等宽边界可以传标量但 MATLAB 的隐式扩展在某些老版本上不一定可靠建议始终传向量。另外evaluate_population里用了循环逐个求值如果目标函数是向量化友好的比如 ZDT 系列的数学表达式改成一次性矩阵运算可以把耗时降低一个数量级。4. 把 NSWOA 用到真实问题物流配送路径优化的建模与参数调整4.1 多目标怎么建模成本与时间的双目标函数设计思路纯测试函数的 Pareto 前沿是已知的拿来验证算法正确性没问题。但工程场景里目标函数往往不是标准测试函数建模这一步就决定了优化效果的上限。以物流配送路径优化为例常见做法是设置两个目标总运输成本和总配送时间。成本包括车辆固定成本和距离相关成本时间包括行驶时间和各节点的服务时间。目标函数写成 MATLAB 函数时输入是一个向量形式的路径编码需要先解码成实际的配送顺序function f logistics_cost(x) % x 是连续变量编码的配送顺序需要解码 % 假设有K辆车、C个客户点x是C维的实数向量 C length(x); [~, order] sort(x); % 排序解码把连续位置映射为配送顺序 global dist_matrix demand service_time vehicle_capacity; total_distance 0; total_time 0; vehicle_count 1; current_load 0; current_pos 1; % 1是配送中心 current_time 0; for i 1:C customer order(i); if current_load demand(customer) vehicle_capacity total_distance total_distance dist_matrix(current_pos, 1); total_time total_time dist_matrix(current_pos, 1) / speed service_time(current_pos); vehicle_count vehicle_count 1; current_pos 1; current_load 0; current_time 0; end total_distance total_distance dist_matrix(current_pos, customer); total_time total_time dist_matrix(current_pos, customer) / speed service_time(customer); current_pos customer; current_load current_load demand(customer); end total_distance total_distance dist_matrix(current_pos, 1); total_time total_time dist_matrix(current_pos, 1) / speed; f [total_distance vehicle_count * fixed_cost_per_vehicle, total_time]; end排序解码是路径优化里最常见也最稳定的编码方式。直接把实数向量按升序排列得到一个客户访问顺序再根据容量约束切分成多条路径。这样设计的好处是 NSWOA 的连续位置更新算子不需要任何修改直接跑就行。代价是搜索空间有冗余——多个不同实数向量可能对应同一个排序但这在元启发式算法里通常可以接受实际效果比离散编码加专门算子的实现更省事。4.2 NSWOA 在这个场景下的 3 个必调参数不同问题性质对参数敏感度差异很大物流路径和 ZDT 函数至少有两点不同变量之间有关联性顺序表达目标函数不是平滑函数。这两点决定了参数要重新调。第一个参数是a的递减方式。标准 WOA 里a从 2 线性降到 0但路径优化里后期如果a已经很小所有鲸鱼都朝 leader 收缩种群多样性骤降很容易陷入局部路径组合。常见做法是把线性递减改成非线性前 40% 迭代保持a在 1.5 以上加强探索后 60% 快速降到 0 加强开发。实现上只需要把主循环里的a计算改成if t / MaxIter 0.4 a 2 - (t / MaxIter) / 0.4 * 0.5; % 前期从2缓慢降到1.5 else a 1.5 - ((t / MaxIter) - 0.4) / 0.6 * 1.5; % 后期从1.5快速降到0 end第二个是螺旋常数b。标准值 1 在连续函数上表现好但路径编码的欧氏距离含义会被排序操作破坏螺旋更新到实数空间产生的扰动和排序解码后的路径变化不是线性对应的。我一般把b调小到 0.5让螺旋的半径变化更平缓避免排序结果的大幅跳动。第三个是存档里 leader 的选择策略。前面提到选拥挤度最大的解但在物流场景里如果两个目标之间的量纲差距过大比如成本是万元级时间是小时级拥挤度计算归一化后可能让某几个极端解永远占据最大距离位置。更好的做法是定期轮换 leader每个迭代归档后把 F1 前沿按拥挤度排序取前 20% 的解随机选一个作为 leader。这个改动通常能在保持收敛性的同时让前沿分布更均匀。4.3 多目标优化的评价指标IGD 与 HV 怎么看优化跑完不能只靠眼睛看散点图需要数值指标评估求解质量。最常用的两个指标是 IGDInverted Generational Distance和 HVHypervolume。IGD 衡量的是算法得到的解集到真实 Pareto 前沿的平均距离值越小越好。它需要真实前沿测试函数可以从文献里下载工程问题一般没有真实前沿只能用近似前沿代替。HV 衡量的是解集在目标空间里覆盖的体积值越大越好不需要真实前沿只需要一个参考点通常取各目标的最差值。function hv hypervolume(fitness, ref_point) % fitness: 已归一化到[0,1]的目标值假设是最小化问题 % ref_point: 参考点各目标的上界 N size(fitness, 1); % 按第一个目标排序 [~, idx] sort(fitness(:,1), descend); fitness fitness(idx, :); hv 0; prev_f2 ref_point(2); for i 1:N if fitness(i,2) prev_f2 hv hv (fitness(i,1) - 0) * (prev_f2 - fitness(i,2)); prev_f2 fitness(i,2); end end hv hv (0 - 0) * (prev_f2 - 0); % 修正最后一块面积 end这个实现是最基础的双目标 HV 计算思路是按第一个目标降序排列逐步累加矩形面积。三目标以上的 HV 计算比较复杂通常用蒙特卡洛采样近似。跑对比实验时建议每次都记录 HV 和 IGD算均值和标准差至少独立运行 10 次因为元启发式算法单次运行有随机性一次跑得好不代表算法真的好。5. NSWOA 与 NSGA-II 的对比实验跑同一组测试问题看差异5.1 实验设置相同的预算才是公平的比较拿 NSWOA 和 NSGA-II 做对比时最容易犯的错误是只比迭代次数不比函数评估次数。NSGA-II 一代要做一次交叉变异和选择NSWOA 一代要做一次鲸鱼位置更新和归档裁剪单代计算量不同只比代数没有意义。公平的做法是固定总的函数评估次数比如 N×MaxIter 50000 次两个算法都用这个预算。MATLAB 里统计函数评估次数需要在目标函数内部加计数器或者用全局变量传值。把目标函数封装一层function f counting_wrapper(x) global eval_count; eval_count eval_count size(x, 1); f ZDT1(x); end每次调用evaluate_population时计数自动累加。两个算法的终止条件都改成eval_count 50000。运行时发现 NSWOA 的归档裁剪阶段会重复评估存档里的解计数会偏高所以我在归档更新时直接沿用已有的适应度值不重新调用目标函数。5.2 在 ZDT1、ZDT2、DTLZ2 上的表现解读ZDT1 是凸前沿ZDT2 是非凸前沿DTLZ2 是三目标问题。三组测试问题分别考察算法在不同前沿形状下的表现。我跑下来的经验结果是在 ZDT1 上 NSWOA 的收敛速度明显快于 NSGA-II前 100 代 IGD 下降曲线更陡这是因为 Whale 的螺旋更新在连续空间里搜索效率高包围机制对凸前沿有天然的逼近优势。但在 ZDT2 非凸前沿上NSWOA 容易出现前沿覆盖不完整的问题集中在两端而中间稀疏。原因是非凸区域对探索能力要求更高而 NSWOA 的随机搜索算子比例在后期被压缩得很小。解决办法是把随机搜索的触发概率从 0.5 提高到 0.6也就是位置更新里p(i) 0.6时进入包围或随机搜索p(i) 0.6时做螺旋更新。DTLZ2 三目标上NSWOA 的排名表现主要受存档规模影响。ArchiveSize 等于种群规模时前沿的三维分布比较均匀但边缘区域解的数量明显少于 NSGA-II。如果做 DTLZ 系列对比建议 ArchiveSize 比种群规模大 20%~50%否则三维前沿的边界解容易在裁剪中丢失。5.3 常用的调试技巧跑偏时如何定位问题NSWOA 出问题时症状通常分为三类。第一类是收敛到某一点Pareto 前沿缩成一个点或一条线这对应种群多样性崩溃。检查a的递减曲线是不是太平滑导致探索不足检查 leader 选择是不是每次选同一个解。日志里加一行打印存档中互不支配的解数量如果连续 20 代这个数量没有变化考虑引入变异算子以一个很小的概率比如 0.05随机重置个别维度。第二类是前沿分布不均中间缺一大块。这时要检查拥挤度距离计算的归一化逻辑是否正确尤其是量纲差异大的目标。把每个目标的距离贡献单独打印出来看是不是某一个目标贡献了 90% 以上的距离值。如果是对归一化分母做一次对数缩放压缩量纲差异。第三类是收敛太慢迭代到一半 IGD 还在高位。优先检查目标函数是否有大量平坦区域导致鲸鱼在排序解码后位置变化不敏感。解决方法是添加辅助梯度在目标函数里加入和距离相关的小扰动项引导鲸鱼朝有变化的方向移动。这个技巧本质上是对目标函数做平滑但不能改得面目全非扰动幅度控制在原函数量级的 1% 以内。6. 最后一章用实时可视化和参数热力图调 NSWOA 的实用技巧调试 NSWOA 时只看最后的静态 Pareto 图远远不够尤其是迭代中期的种群动态根本看不出来。我常用的做法是在主循环里嵌入一个实时绘图窗口每隔一定代数刷新一次当前存档同时记录每代存档的 HV 值。这样能直观看到探索和开发的动态过程前期点应该散得很快中期开始聚拢后期沿着前沿缓慢滑动。具体做法是把 NSWOA 主函数的history部分扩展增加一个回调参数或者在局部脚本里做成嵌套函数。下面的代码片段展示了一种轻量实现% 在主循环末尾添加 if mod(t, 20) 0 || t MaxIter subplot(1, 2, 1); plot(archive_fitness(:,1), archive_fitness(:,2), b., MarkerSize, 8); xlabel(f1); ylabel(f2); title(sprintf(Iteration %d / %d, t, MaxIter)); drawnow limitrate; subplot(1, 2, 2); hv_val hypervolume(normalize_fitness(archive_fitness), [1, 1]); history_hv(t) hv_val; plot(1:t, history_hv(1:t), r-, LineWidth, 1.5); xlabel(Iteration); ylabel(HV); title(Hypervolume History); drawnow; end注意drawnow limitrate和drawnow的区别。limitrate会限制刷新频率渲染开销大时用否则实时绘图的耗时可能超过算法本身的耗时干扰计时。画图前对适应度归一化很重要不然 HV 的数值随目标量纲变化没有可比性。另一个值得做的静态检查是参数热力图。把种群规模 N 和 ArchiveSize 分别设成网格比如 N 取 50、100、150、200ArchiveSize 取 N 的 50%、75%、100%、125%对每个组合跑 10 次取平均 HV画成热力图。这样做有两个价值一是直观看到参数交互效应比如 N100 时 ArchiveSize 对结果影响不大但 N50 时 ArchiveSize 不足会明显降低 HV二是为同类问题提供一个可迁移的调参起点下次遇到相似规模的问题可以直接把网格参数作为默认值省去重复调试。关于 MATLAB 版本NSWOA 的代码本身对版本要求不高2018a 以上都能跑。真正值得注意的是evaluate_population里的循环结构如果目标函数是高度向量化的数学表达式把循环改成矩阵运算可以提速 10 倍以上特别适合参数热力图这种需要大量重复调用的场景。改法是把种群矩阵一次性传给目标函数目标函数内部按行操作返回一个 N×M 的适应度矩阵。但路径优化这类需要解码的现实问题往往无法完全向量化这时优先保证可读性性能靠减少无意义的重复计算来提升。本文还有配套的精品资源点击获取
NSWOA多目标鲸鱼优化算法MATLAB实现详解
简介基于鲸鱼算法NSWOA的多目标优化求解Matlab代码包面向智能优化算法、神经网络预测、路径规划等方向的本科生、研究生及科研人员适合算法对比、课程实验与论文复现等场景。代码包以NSWOA为核心完整集成非支配排序、变量初始化、目标函数评价、种群更新、替换染色体和结果绘图等子模块并内置zdt1、zdt2、zdt3等多组常用测试函数运行后可直接查看多目标Pareto前沿的收敛效果直观展示算法在多个目标间的权衡能力有助于理解多目标进化算法的迭代逻辑。压缩包共13个文件其中11个为.m源文件便于逐行阅读和修改另附1个说明文档和1个示例效果图整体仅24KB结构精简、核心算法集中能快速定位到主程序与关键子函数。以上代码同时支持Matlab 2014及2019a环境且描述中注明内含运行结果适合学习时对照验证。目前已有259人学习下载尤其适合刚接触多目标优化、希望从代码层面入手掌握鲸鱼算法原理的教研人员。1. 从单目标到多目标鲸鱼算法为什么要改造成 NSWOA做优化的人都有个直觉单目标问题解出来是个点多目标问题解出来是条线、一个面甚至一片不规则的云。鲸鱼算法Whale Optimization Algorithm, WOA是 2016 年提出的群体智能算法模拟座头鲸的泡泡网捕食策略在单目标连续优化上表现亮眼收敛快、参数少、实现简单。但把 WOA 直接扔到多目标问题上会立刻暴露短板它没有维护种群多样性的机制也没有判断一个解是否被另一个解支配的框架最后所有鲸鱼会朝一个方向挤得到一堆彼此雷同的解。NSWOA 的全称是 Non-dominated Sorting Whale Optimization Algorithm思路很直白把 NSGA-II 里那套非支配排序和拥挤度距离搬进 WOA 的寻优循环里让鲸鱼在搜索过程中既保留 WOA 的位置更新机制又用 Pareto 支配关系筛选精英解。这样得到的不是唯一最优解而是一组分布均匀、覆盖前沿的折中方案。工程上的价值在于物流配送路径、车间调度、资源分配这类问题往往有多个互相冲突的目标运维工程师、算法工程师和科研人员需要在一组 Pareto 解里根据业务偏好做决策而不是被迫接受一个单点答案。这篇文章就用 MATLAB 把 NSWOA 的完整实现拆开讲从核心算子到参数设置再到和 NSGA-II 的对比一步步写清楚。2. 核心机制拆解NSWOA 如何在迭代中同时顾全收敛性与分布性2.1 非支配排序给每个解一个「等级」而不是一个「分数」多目标优化和单目标最大的区别在于比较方式。单目标可以用数值大小直接排序多目标需要先判断支配关系。一个解 x1 支配 x2当且仅当 x1 在所有目标上都不劣于 x2且至少在一个目标上严格优于 x2。不满足这个条件时两个解互不支配它们属于同一个 Pareto 前沿。NSWOA 在每次迭代后都要对当前种群做一遍非支配排序把解分成 F1、F2、F3 多个层级F1 是当前最优前沿。实现上常见做法是用一个双层循环记录每个解被谁支配、支配了多少个解然后逐层提取前沿。MATLAB 里写起来大约是这样function front non_dominated_sort(fitness) [N, M] size(fitness); % N个解M个目标 dominate_count zeros(N, 1); % 记录支配当前解的解数量 dominated_set cell(N, 1); % 记录当前解支配的解集合 front zeros(N, 1); % 每个解的Pareto层级 for i 1:N for j i1:N if dominates(fitness(i,:), fitness(j,:)) dominated_set{i} [dominated_set{i}, j]; dominate_count(j) dominate_count(j) 1; elseif dominates(fitness(j,:), fitness(i,:)) dominated_set{j} [dominated_set{j}, i]; dominate_count(i) dominate_count(i) 1; end end end current_front find(dominate_count 0); % 第一前沿 rank 1; while ~isempty(current_front) for idx current_front front(idx) rank; end next_front []; for idx current_front for dominated_idx dominated_set{idx} dominate_count(dominated_idx) dominate_count(dominated_idx) - 1; if dominate_count(dominated_idx) 0 next_front [next_front, dominated_idx]; end end end rank rank 1; current_front next_front; end end这段代码的核心逻辑是先统计每个解被多少解支配被零个解支配的划定到 F1然后把 F1 的解从支配关系中「剔除」更新计数再找下一批零支配解重复直到所有解都有等级。front数组最终给出每个解所在的层级。注意 MATLAB 的下标从 1 开始代码里所有索引都对应这个习惯直接用 0 会报错。2.2 拥挤度距离在同一前沿里挑「住得最宽敞」的鲸鱼层级只解决了解的优劣问题没解决分布问题。如果筛选精英时只看层级同一前沿内随机选解很容易堆在一起Pareto 前沿中间缺一大块边缘却挤满了点。拥挤度距离的思路是对同一前沿的解按每个目标分别排序计算每个解与相邻两个解在目标空间的距离之和距离越大说明周围越空旷越值得保留。function distance crowding_distance(fitness, front) [N, M] size(fitness); distance zeros(N, 1); for level 1:max(front) idx find(front level); if length(idx) 2 distance(idx) inf; % 少于3个解时都设为无穷大全部保留 continue; end for m 1:M [~, order] sort(fitness(idx, m)); sorted_idx idx(order); distance(sorted_idx(1)) inf; distance(sorted_idx(end)) inf; f_min fitness(sorted_idx(1), m); f_max fitness(sorted_idx(end), m); if f_max - f_min 1e-12 continue; % 目标值全相同距离记为0 end for k 2:length(sorted_idx)-1 distance(sorted_idx(k)) distance(sorted_idx(k)) ... (fitness(sorted_idx(k1), m) - fitness(sorted_idx(k-1), m)) / (f_max - f_min); end end end end边界解的距离被直接设为inf这个细节很关键。边界解承担着扩展前沿的任务即使中间区域密度很大也必须把两端保住否则下一次迭代前沿就向内收缩了。距离计算时要按目标分别归一化不然量纲大的目标会主导距离值量纲小的目标在筛选中几乎不起作用。如果某个目标的所有解值都相同分母趋近于零代码里用continue跳过把这个目标的贡献记为零避免产生 NaN。2.3 鲸鱼位置更新机制包围、泡泡网与随机搜索三选一NSWOA 保留了 WOA 的位置更新策略每次迭代每头鲸鱼以一定概率选择三种行为之一。第一种是包围猎物位置向当前最优解靠拢第二种是泡泡网攻击用螺旋方程模拟座头鲸吐气泡向上盘旋的动作第三种是随机搜索向种群中随机一个个体移动维持全局探索能力。function new_positions whale_update(positions, best_pos, a, b, l, p, r1, r2) [N, dim] size(positions); new_positions zeros(N, dim); for i 1:N if p(i) 0.5 if abs(a(i)) 1 % 收缩包围向当前最优解靠近 new_positions(i,:) best_pos - a(i) * abs(b(i) * best_pos - positions(i,:)); else % 随机搜索选择一个随机个体作为参照 rand_idx randi(N); rand_whale positions(rand_idx,:); new_positions(i,:) rand_whale - a(i) * abs(b(i) * rand_whale - positions(i,:)); end else % 螺旋更新以螺旋路径逼近最优解 dist abs(best_pos - positions(i,:)); new_positions(i,:) dist .* exp(b * l(i)) .* cos(2 * pi * l(i)) best_pos; end end end参数a随迭代从 2 线性递减到 0控制了探索和开发的平衡。abs(a) 1时鲸鱼向最优解收缩abs(a) 1时跳出局部陷阱做全局随机搜索。注意这里用的是randi(N)随机选参照个体和单目标 WOA 完全一致。问题在于多目标场景下「最优解」不再唯一——F1 前沿里每个解都是最优选哪个作为best_pos直接影响搜索方向。NSWOA 的常见做法是在 F1 前沿里随机挑一个每次迭代不同鲸鱼可能朝向不同的 Pareto 解移动这样种群不至于全部塌缩到同一个点。我在实际写代码时会在每次迭代开始时先对 F1 前沿做一次洗牌然后每头鲸鱼按自己的序号取模选参照比每次都重新调用randi更稳定。3. 在 MATLAB 中搭建 NSWOA 求解多目标问题的完整流程3.1 外部存档与精英保留NSWOA 的归档策略怎么写NSWOA 和很多多目标元启发式算法的区别在于是否有外部存档。NSGA-II 是靠精英选择机制把父代和子代合并后筛选NSWOA 更常见的做法是维护一个独立的存档集用非支配排序加拥挤度距离控制存档规模。每次迭代结束后把当前种群和存档合并做一次非支配排序取前ArchiveSize个解作为新存档。我做这个部分的代码时发现最值得注意的坑是越界处理。多目标问题里变量的边界约束很常见但进化算法更新完位置后经常有解跑出边界。直接把这些解扔掉会减少种群多样性简单粗暴地「截断」到边界又会让大量解堆在边界线上影响拥挤度计算的真实性。我一般用随机反射策略超出边界的维度按边界宽度折返类似光在镜面反射这样既保留了解的位置又避免边界堆积。function x boundary_reflect(x, lb, ub) range ub - lb; over_high x ub; over_low x lb; x(over_high) ub(over_high) - mod(x(over_high) - ub(over_high), range(over_high)); x(over_low) lb(over_low) mod(lb(over_low) - x(over_low), range(over_low)); idx (x lb) | (x ub); % 极端情况下取模后仍然越界 x(idx) lb(idx) rand(size(lb(idx))) .* range(idx); if any(x(:) lb(:)) || any(x(:) ub(:)) x max(min(x, ub), lb); % 兜底保证在边界内 end end这段代码里的逻辑顺序很重要先做反射再对极端情况做随机初始化最后兜底截断。如果一上来就截断反射策略就失去了意义。实际跑的时候发现反射后仍然越界的概率很低但一旦发生通常是变量维度特别多、搜索步长特别大的情况兜底处理是必要的安全网。3.2 种子函数与主循环变量数、目标数、代数怎么往代码里填用户拿到的 NSWOA 代码压缩包里通常包含NSWOA.m主函数、目标函数文件和一些绘图脚本。我建议在跑任何测试问题之前先搞清楚主函数的输入输出接口。典型的主循环骨架长这样function [archive, archive_fitness, history] NSWOA(fun, dim, lb, ub, N, MaxIter, ArchiveSize) positions lb rand(N, dim) .* (ub - lb); fitness evaluate_population(fun, positions); archive positions; archive_fitness fitness; history zeros(MaxIter, 1); for t 1:MaxIter a 2 - 2 * t / MaxIter; % 线性递减参数 b 1; % 螺旋形状常数 best_idx select_leader(archive, archive_fitness); % 从存档里选领导者 p rand(N, 1); r1 rand(N, dim); r2 rand(N, dim); l -1 2 * rand(N, 1); A 2 * a .* r1 - a; C 2 .* r2; new_positions zeros(N, dim); for i 1:N if p(i) 0.5 if abs(A(i,1)) 1 new_positions(i,:) archive(best_idx,:) - A(i,:) .* abs(C(i,:) .* archive(best_idx,:) - positions(i,:)); else rand_idx randi(N); new_positions(i,:) positions(rand_idx,:) - A(i,:) .* abs(C(i,:) .* positions(rand_idx,:) - positions(i,:)); end else dist abs(archive(best_idx,:) - positions(i,:)); new_positions(i,:) dist .* exp(b * l(i)) .* cos(2 * pi * l(i)) archive(best_idx,:); end end new_positions boundary_reflect(new_positions, lb, ub); new_fitness evaluate_population(fun, new_positions); % 合并种群与存档非支配排序裁剪 combined_pos [positions; new_positions; archive]; combined_fit [fitness; new_fitness; archive_fitness]; front non_dominated_sort(combined_fit); dist crowding_distance(combined_fit, front); % 按层级和拥挤度排序取前N个 [~, sort_idx] sortrows([front, -dist]); positions combined_pos(sort_idx(1:N), :); fitness combined_fit(sort_idx(1:N), :); % 从排序后的解里提取存档 archive positions(front(sort_idx(1:N)) 1, :); archive_fitness fitness(front(sort_idx(1:N)) 1, :); if size(archive, 1) ArchiveSize archive trim_archive(archive, archive_fitness, ArchiveSize); archive_fitness evaluate_population(fun, archive); end history(t) sum(min(archive_fitness, [], 1)); % 记录目标值的动态 end end这个主循环里有个和纯 NSGA-II 不同的设计位置更新时参照的不是「本代的最优解」而是「存档里的最优解」。存档里保留的是历史迭代中发现的全部非支配解比当前代的 F1 前沿更全面。select_leader的实现通常是做一次拥挤度距离计算取距离最大的解作为 leader这样能引导种群朝前沿的稀疏区域搜索而不是扎堆在已经密集的区域。主函数的参数设置直接影响求解效果。下面这张表是经验值适合 ZDT 系列和 DTLZ 系列测试函数实际工程问题需要按变量维度和目标数量调整参数推荐范围说明种群规模 N50~200目标函数计算代价高时取小值最大迭代 MaxIter200~1000变量维度高时取大值存档规模 ArchiveSizeN 的 50%~100%过小会丢失前沿细节螺旋常数 b0.5~1.2控制螺旋的紧密程度目标数 M2~3超过 3 个目标建议用 NSGA-III 的参考点机制实际运行中如果发现最后得到的 Pareto 前沿不够连续优先增大ArchiveSize而不是MaxIter。存档空间太小即使迭代再多留下的非支配解也被迫裁剪前沿必然稀疏。3.3 从压缩包到可运行NSWOA 代码的文件结构与调用方式拿到 NSWOA 的压缩包解压后通常能看到这些文件NSWOA.m是主算法ZDT1.m或类似的测试函数文件PlotPareto.m是画图脚本有时还有一个main.m把所有步骤串起来。先用测试函数验证算法正确性再替换成自己的目标函数这是最稳的上手顺序。% main.m 示例 clc; clear; close all; % 定义ZDT1测试函数两个目标30个变量 dim 30; lb zeros(1, dim); ub ones(1, dim); N 100; MaxIter 500; ArchiveSize 100; [archive, archive_fitness, history] NSWOA(ZDT1, dim, lb, ub, N, MaxIter, ArchiveSize); figure; plot(archive_fitness(:,1), archive_fitness(:,2), o); xlabel(f1); ylabel(f2); title(NSWOA Pareto Front - ZDT1); % 计算IGD指标评估质量 true_front load(ZDT1_true_front.txt); igd_value IGD(archive_fitness, true_front); fprintf(IGD %.4f\n, igd_value);调用时注意两点目标函数必须返回一个列向量每个元素是一个目标值NSWOA 内部统一按行组织种群、按列组织目标lb和ub若是全部变量等宽边界可以传标量但 MATLAB 的隐式扩展在某些老版本上不一定可靠建议始终传向量。另外evaluate_population里用了循环逐个求值如果目标函数是向量化友好的比如 ZDT 系列的数学表达式改成一次性矩阵运算可以把耗时降低一个数量级。4. 把 NSWOA 用到真实问题物流配送路径优化的建模与参数调整4.1 多目标怎么建模成本与时间的双目标函数设计思路纯测试函数的 Pareto 前沿是已知的拿来验证算法正确性没问题。但工程场景里目标函数往往不是标准测试函数建模这一步就决定了优化效果的上限。以物流配送路径优化为例常见做法是设置两个目标总运输成本和总配送时间。成本包括车辆固定成本和距离相关成本时间包括行驶时间和各节点的服务时间。目标函数写成 MATLAB 函数时输入是一个向量形式的路径编码需要先解码成实际的配送顺序function f logistics_cost(x) % x 是连续变量编码的配送顺序需要解码 % 假设有K辆车、C个客户点x是C维的实数向量 C length(x); [~, order] sort(x); % 排序解码把连续位置映射为配送顺序 global dist_matrix demand service_time vehicle_capacity; total_distance 0; total_time 0; vehicle_count 1; current_load 0; current_pos 1; % 1是配送中心 current_time 0; for i 1:C customer order(i); if current_load demand(customer) vehicle_capacity total_distance total_distance dist_matrix(current_pos, 1); total_time total_time dist_matrix(current_pos, 1) / speed service_time(current_pos); vehicle_count vehicle_count 1; current_pos 1; current_load 0; current_time 0; end total_distance total_distance dist_matrix(current_pos, customer); total_time total_time dist_matrix(current_pos, customer) / speed service_time(customer); current_pos customer; current_load current_load demand(customer); end total_distance total_distance dist_matrix(current_pos, 1); total_time total_time dist_matrix(current_pos, 1) / speed; f [total_distance vehicle_count * fixed_cost_per_vehicle, total_time]; end排序解码是路径优化里最常见也最稳定的编码方式。直接把实数向量按升序排列得到一个客户访问顺序再根据容量约束切分成多条路径。这样设计的好处是 NSWOA 的连续位置更新算子不需要任何修改直接跑就行。代价是搜索空间有冗余——多个不同实数向量可能对应同一个排序但这在元启发式算法里通常可以接受实际效果比离散编码加专门算子的实现更省事。4.2 NSWOA 在这个场景下的 3 个必调参数不同问题性质对参数敏感度差异很大物流路径和 ZDT 函数至少有两点不同变量之间有关联性顺序表达目标函数不是平滑函数。这两点决定了参数要重新调。第一个参数是a的递减方式。标准 WOA 里a从 2 线性降到 0但路径优化里后期如果a已经很小所有鲸鱼都朝 leader 收缩种群多样性骤降很容易陷入局部路径组合。常见做法是把线性递减改成非线性前 40% 迭代保持a在 1.5 以上加强探索后 60% 快速降到 0 加强开发。实现上只需要把主循环里的a计算改成if t / MaxIter 0.4 a 2 - (t / MaxIter) / 0.4 * 0.5; % 前期从2缓慢降到1.5 else a 1.5 - ((t / MaxIter) - 0.4) / 0.6 * 1.5; % 后期从1.5快速降到0 end第二个是螺旋常数b。标准值 1 在连续函数上表现好但路径编码的欧氏距离含义会被排序操作破坏螺旋更新到实数空间产生的扰动和排序解码后的路径变化不是线性对应的。我一般把b调小到 0.5让螺旋的半径变化更平缓避免排序结果的大幅跳动。第三个是存档里 leader 的选择策略。前面提到选拥挤度最大的解但在物流场景里如果两个目标之间的量纲差距过大比如成本是万元级时间是小时级拥挤度计算归一化后可能让某几个极端解永远占据最大距离位置。更好的做法是定期轮换 leader每个迭代归档后把 F1 前沿按拥挤度排序取前 20% 的解随机选一个作为 leader。这个改动通常能在保持收敛性的同时让前沿分布更均匀。4.3 多目标优化的评价指标IGD 与 HV 怎么看优化跑完不能只靠眼睛看散点图需要数值指标评估求解质量。最常用的两个指标是 IGDInverted Generational Distance和 HVHypervolume。IGD 衡量的是算法得到的解集到真实 Pareto 前沿的平均距离值越小越好。它需要真实前沿测试函数可以从文献里下载工程问题一般没有真实前沿只能用近似前沿代替。HV 衡量的是解集在目标空间里覆盖的体积值越大越好不需要真实前沿只需要一个参考点通常取各目标的最差值。function hv hypervolume(fitness, ref_point) % fitness: 已归一化到[0,1]的目标值假设是最小化问题 % ref_point: 参考点各目标的上界 N size(fitness, 1); % 按第一个目标排序 [~, idx] sort(fitness(:,1), descend); fitness fitness(idx, :); hv 0; prev_f2 ref_point(2); for i 1:N if fitness(i,2) prev_f2 hv hv (fitness(i,1) - 0) * (prev_f2 - fitness(i,2)); prev_f2 fitness(i,2); end end hv hv (0 - 0) * (prev_f2 - 0); % 修正最后一块面积 end这个实现是最基础的双目标 HV 计算思路是按第一个目标降序排列逐步累加矩形面积。三目标以上的 HV 计算比较复杂通常用蒙特卡洛采样近似。跑对比实验时建议每次都记录 HV 和 IGD算均值和标准差至少独立运行 10 次因为元启发式算法单次运行有随机性一次跑得好不代表算法真的好。5. NSWOA 与 NSGA-II 的对比实验跑同一组测试问题看差异5.1 实验设置相同的预算才是公平的比较拿 NSWOA 和 NSGA-II 做对比时最容易犯的错误是只比迭代次数不比函数评估次数。NSGA-II 一代要做一次交叉变异和选择NSWOA 一代要做一次鲸鱼位置更新和归档裁剪单代计算量不同只比代数没有意义。公平的做法是固定总的函数评估次数比如 N×MaxIter 50000 次两个算法都用这个预算。MATLAB 里统计函数评估次数需要在目标函数内部加计数器或者用全局变量传值。把目标函数封装一层function f counting_wrapper(x) global eval_count; eval_count eval_count size(x, 1); f ZDT1(x); end每次调用evaluate_population时计数自动累加。两个算法的终止条件都改成eval_count 50000。运行时发现 NSWOA 的归档裁剪阶段会重复评估存档里的解计数会偏高所以我在归档更新时直接沿用已有的适应度值不重新调用目标函数。5.2 在 ZDT1、ZDT2、DTLZ2 上的表现解读ZDT1 是凸前沿ZDT2 是非凸前沿DTLZ2 是三目标问题。三组测试问题分别考察算法在不同前沿形状下的表现。我跑下来的经验结果是在 ZDT1 上 NSWOA 的收敛速度明显快于 NSGA-II前 100 代 IGD 下降曲线更陡这是因为 Whale 的螺旋更新在连续空间里搜索效率高包围机制对凸前沿有天然的逼近优势。但在 ZDT2 非凸前沿上NSWOA 容易出现前沿覆盖不完整的问题集中在两端而中间稀疏。原因是非凸区域对探索能力要求更高而 NSWOA 的随机搜索算子比例在后期被压缩得很小。解决办法是把随机搜索的触发概率从 0.5 提高到 0.6也就是位置更新里p(i) 0.6时进入包围或随机搜索p(i) 0.6时做螺旋更新。DTLZ2 三目标上NSWOA 的排名表现主要受存档规模影响。ArchiveSize 等于种群规模时前沿的三维分布比较均匀但边缘区域解的数量明显少于 NSGA-II。如果做 DTLZ 系列对比建议 ArchiveSize 比种群规模大 20%~50%否则三维前沿的边界解容易在裁剪中丢失。5.3 常用的调试技巧跑偏时如何定位问题NSWOA 出问题时症状通常分为三类。第一类是收敛到某一点Pareto 前沿缩成一个点或一条线这对应种群多样性崩溃。检查a的递减曲线是不是太平滑导致探索不足检查 leader 选择是不是每次选同一个解。日志里加一行打印存档中互不支配的解数量如果连续 20 代这个数量没有变化考虑引入变异算子以一个很小的概率比如 0.05随机重置个别维度。第二类是前沿分布不均中间缺一大块。这时要检查拥挤度距离计算的归一化逻辑是否正确尤其是量纲差异大的目标。把每个目标的距离贡献单独打印出来看是不是某一个目标贡献了 90% 以上的距离值。如果是对归一化分母做一次对数缩放压缩量纲差异。第三类是收敛太慢迭代到一半 IGD 还在高位。优先检查目标函数是否有大量平坦区域导致鲸鱼在排序解码后位置变化不敏感。解决方法是添加辅助梯度在目标函数里加入和距离相关的小扰动项引导鲸鱼朝有变化的方向移动。这个技巧本质上是对目标函数做平滑但不能改得面目全非扰动幅度控制在原函数量级的 1% 以内。6. 最后一章用实时可视化和参数热力图调 NSWOA 的实用技巧调试 NSWOA 时只看最后的静态 Pareto 图远远不够尤其是迭代中期的种群动态根本看不出来。我常用的做法是在主循环里嵌入一个实时绘图窗口每隔一定代数刷新一次当前存档同时记录每代存档的 HV 值。这样能直观看到探索和开发的动态过程前期点应该散得很快中期开始聚拢后期沿着前沿缓慢滑动。具体做法是把 NSWOA 主函数的history部分扩展增加一个回调参数或者在局部脚本里做成嵌套函数。下面的代码片段展示了一种轻量实现% 在主循环末尾添加 if mod(t, 20) 0 || t MaxIter subplot(1, 2, 1); plot(archive_fitness(:,1), archive_fitness(:,2), b., MarkerSize, 8); xlabel(f1); ylabel(f2); title(sprintf(Iteration %d / %d, t, MaxIter)); drawnow limitrate; subplot(1, 2, 2); hv_val hypervolume(normalize_fitness(archive_fitness), [1, 1]); history_hv(t) hv_val; plot(1:t, history_hv(1:t), r-, LineWidth, 1.5); xlabel(Iteration); ylabel(HV); title(Hypervolume History); drawnow; end注意drawnow limitrate和drawnow的区别。limitrate会限制刷新频率渲染开销大时用否则实时绘图的耗时可能超过算法本身的耗时干扰计时。画图前对适应度归一化很重要不然 HV 的数值随目标量纲变化没有可比性。另一个值得做的静态检查是参数热力图。把种群规模 N 和 ArchiveSize 分别设成网格比如 N 取 50、100、150、200ArchiveSize 取 N 的 50%、75%、100%、125%对每个组合跑 10 次取平均 HV画成热力图。这样做有两个价值一是直观看到参数交互效应比如 N100 时 ArchiveSize 对结果影响不大但 N50 时 ArchiveSize 不足会明显降低 HV二是为同类问题提供一个可迁移的调参起点下次遇到相似规模的问题可以直接把网格参数作为默认值省去重复调试。关于 MATLAB 版本NSWOA 的代码本身对版本要求不高2018a 以上都能跑。真正值得注意的是evaluate_population里的循环结构如果目标函数是高度向量化的数学表达式把循环改成矩阵运算可以提速 10 倍以上特别适合参数热力图这种需要大量重复调用的场景。改法是把种群矩阵一次性传给目标函数目标函数内部按行操作返回一个 N×M 的适应度矩阵。但路径优化这类需要解码的现实问题往往无法完全向量化这时优先保证可读性性能靠减少无意义的重复计算来提升。本文还有配套的精品资源点击获取
