1. 从单目标到多目标最短路径问题为什么需要MOMPA路径规划这个领域做久了你会发现一个挺有意思的现象大多数人一开始接触的都是单目标最短路径经典Dijkstra、A*顶多再上个遗传算法优化一下搞定一条从起点到终点的最短折线就收工。但真实场景里问题远没有那么简单。举个例子我前两年做一个无人船近海航线规划的项目甲方提的需求是“航程最短”。等真把海图数据拿过来风速场、海流场、暗礁区、航道限速区叠在一起才发现如果只按几何距离最短来跑算法给出的航线大概率会一头扎进强流区或者贴着禁航区边缘走实际航行时间和燃油消耗反而更高。那种时候你才意识到路径规划本质上是多个互相冲突的目标在同一个解空间里寻优航程要短、安全性要高、能耗要小、转弯幅度要平滑。单个目标函数根本压不住这些需求。这也是我后来重点关注多目标进化算法的原因。而在多目标优化这个方向上**MOMPAMulti-Objective Marine Predators Algorithm多目标海洋捕食者算法**是最近几年表现相当突出的一个。它的前身是单目标的海洋捕食者算法MPA由Faramarzi等人于2020年提出核心思想是模拟海洋生物中捕食者与猎物之间的运动博弈猎物在逃、捕食者在追两者交替采用布朗运动、莱维飞行等不同策略形成一套全局搜索与局部开发自然平衡的寻优机制。MOMPA则是把这种机制扩展到多目标框架下引入了外部档案集Archive和网格选择机制能够在一次运行中输出一组分布均匀的Pareto最优解。这篇文章我就以“求解最短路径问题”为载体完整拆解MOMPA在Matlab里的实现思路包括算法结构、适应度函数设计、代码框架、实验结果分析以及我自己踩过的一些坑。适合有Matlab基础、对进化算法有一定了解、正在做路径规划相关课题或工程项目的读者。2. MOMPA的核心机制三段式搜索与精英选择2.1 海洋捕食者算法的生物学隐喻理解MOMPA之前先得把单目标MPA的运动模型弄清楚。MPA把整个寻优过程建模成捕食者与猎物的相对运动种群中的每个个体要么扮演捕食者要么扮演猎物而两者的运动速度都受到一个核心参数的控制涡流形成因子Fish Aggregating DeviceFAD。算法整个迭代过程被划分成三个阶段每一阶段对应捕食者与猎物速度比的不同区间高速比阶段探索为主迭代初期猎物运动速度比捕食者快此时捕食者采取布朗运动策略大步长漫游搜索整个解空间对应算法中前1/3迭代次数。单位速度比阶段探索与开发平衡中期两者速度相当种群被拆分成两部分一半个体采用布朗运动做局部搜索另一半采用莱维飞行做跳跃式搜索对应中间1/3迭代次数。低速比阶段开发为主迭代后期捕食者速度比猎物快此时捕食者采用莱维飞行策略进行精细开发围绕有希望的局部区域密集搜索对应最后1/3迭代次数。每条搜索路径都会携带一个记忆项记录每个个体历史最优位置。这个设计很像粒子群算法里的个体最优但运动方式的复杂程度要高得多。2.2 从单目标到多目标的扩展Archive和网格选择单目标MPA只会输出一个最优解而MOMPA需要解决“多个互斥目标同时最优”的问题。它的做法和MOPSO、NSGA-II等经典多目标算法类似但有自己的一套实现细节外部档案集维护一个固定容量的解集存放当前找到的非支配解。每轮迭代后把种群中的新解与档案集中的旧解做Pareto支配判断新解若不被档案集中的任何解支配则加入档案集档案集容量溢出时删除处于最拥挤区域中的解。网格选择机制将目标空间划分为均匀的网格每个网格记录包含的解个数。选择“猎物”作为学习对象时倾向于从解密度低的网格中选——这样能保证解的分布性避免全部收敛到Pareto前沿的某一小段。这一套机制的核心思路是“把选择压力从单一最优解转向整个Pareto前沿”所以MOMPA最终输出的不是一条路径而是一组路径——它们各自在航程、安全性、能耗等目标上有不同的折中工程上你可以按实际需求从中挑一条或者用TOPSIS等决策方法自动选一条。2.3 算法总流程一览MOMPA的主循环可以概括为以下几个步骤初始化种群N个个体每个个体是解空间中的一组坐标序列。计算每个个体的多目标适应度值初始化外部档案集。进入三段式迭代根据当前迭代次数与最大迭代次数的比值选择对应的运动策略更新个体位置。边界约束处理。重新计算适应度更新个体记忆最优。更新外部档案集执行网格选择。判断终止条件输出Pareto前沿解集。从工程实现的视角看第2步和第6步是MOMPA区别于单目标MPA的核心增量后面我会在Matlab代码里专门展开。3. 路径规划场景下的适应度函数不止是“最短”3.1 路径的编码方式节点序列与插值策略把MOMPA用于路径规划第一步要解决的是“个体如何编码一条路径”。在栅格地图环境下我这里的做法是起点和终点固定在两点之间均匀插入N个中间路径点每个中间点的x、y坐标作为决策变量这样第i个个体就是一个2N维的实数向量。举个例子地图是100m×100m起点是(5, 5)终点是(95, 95)我插入10个中间点。那么每个个体就是一个20维向量其中奇数位存中间点的x坐标偶数位存y坐标。加上起终点路径上总共有12个节点相邻节点用三次样条插值连成平滑曲线便于避开尖角转弯。这里有一个很关键的细节插值后的点必须再做一次碰撞检测因为相邻两个控制节点之间即使直线连线没有撞障碍物样条插值后的曲线也可能因为过冲而擦到障碍物边缘。为了稳妥我会在每次评估适应度时对整个路径做密集采样比如每隔0.5米采一个点逐个判断是否进入障碍区域。3.2 多目标设计航程、安全距离、平滑度在机器人路径规划里常见的优化目标包括路径长度、与障碍物的最小距离、转弯角度总和或最大转弯角、路径的能量消耗等。我在本文的实现中取了三个最典型的目标目标1路径总长度越小越好这是最基本的目标相邻插值点之间距离累加。定义为[ f_1 \sum_{i1}^{n-1} |P_{i1} - P_i| ]其中 (P_i) 是路径上的密集采样点坐标。注意这里我用的不是控制节点而是采样点——这样计算出来的长度更接近实际行驶轨迹。目标2安全裕度越大越好但代码里取倒数转为最小化这条路径上所有采样点到所有障碍物中心距离的最小值。障碍物在栅格地图中用一个圆心坐标和半径来表示。设采样点坐标为 (P_k)障碍物集合为 (O {o_1, o_2, \dots, o_m})则[ d_{min} \min_{k} \min_{j} \left( |P_k - o_j| - r_j \right) ]由于进化算法习惯上统一处理最小化问题我在实现中把这个目标改写为[ f_2 \frac{1}{d_{min} \epsilon} ]其中 (\epsilon) 是防止除零的小常数取1e-6。当路径贴着障碍物边缘时(f_2) 会变得很大算法会自然淘汰这类高风险路径。目标3路径平滑度越小越好定义为路径上所有相邻线段转角的绝对值和转角越大越不平稳。设三个连续采样点构成的向量为 (\vec{v}i P{i1} - P_i) 和 (\vec{v}{i1} P{i2} - P_{i1})则[ f_3 \sum_{i1}^{n-2} \left| \arccos\left( \frac{\vec{v}i \cdot \vec{v}{i1}}{|\vec{v}i| \cdot |\vec{v}{i1}|} \right) \right| ]3.3 为什么是这三个目标而不是更多可能有读者会问既然是多目标优化为什么不把能耗、时间、通信代价全部加进来目标数量越多不是越全面吗我自己的实践体会是多目标进化算法对目标维数非常敏感。目标维数一旦超过4个Pareto支配关系会迅速弱化——几乎所有解都互不支配选择压力消失算法退化成一个随机搜索器。这在业内被称为“维度灾难”在多目标优化中的体现。所以我的建议是三个目标是最佳实践区间最多不要超过四个。如果你确实有更多维度要考虑比如能耗可以通过加权或约束的方式融合进现有的目标里而不是简单粗暴地加一个新维度。比如可以把能耗近似建模成路径长度与阻力系数的乘积合并进(f_1)中。4. Matlab实现框架从种群初始化到Pareto前沿输出4.1 环境声明与全局参数我用Matlab做这件事的时候第一步是把问题场景参数、算法参数全部定义清楚方便后面统一调整。%% 参数设置 clc; clear; close all; % 地图尺寸 mapSize [100, 100]; % 起点和终点 startPoint [5, 5]; endPoint [95, 95]; % 障碍物定义[x, y, r] obstacles [ 25, 30, 8; 45, 60, 10; 60, 25, 6; 70, 70, 9; 30, 75, 7 ]; %% MOMPA算法参数 numSearchAgents 100; % 种群大小 maxIter 300; % 最大迭代次数 nPoints 10; % 中间控制点数 dim nPoints * 2; % 决策变量维度 archiveSize 100; % 外部档案集容量 nGrids 10; % 每个目标维上的网格数 % 决策变量边界限制在地图范围内 lb [ones(1, nPoints) * 5, ones(1, nPoints) * 5]; ub [ones(1, nPoints) * 95, ones(1, nPoints) * 95];这里的lb和ub设计有个细节我把x和y的边界对称设置了但实际操作中如果你知道某些区域是必不可能的可以单独收紧边界减少搜索空间收敛会快很多。4.2 种群初始化与目标函数定义MPA系列的初始化和大多数群智能算法一样用均匀随机分布在边界内撒点% 初始化种群 positions repmat(lb, numSearchAgents, 1) ... rand(numSearchAgents, dim) .* (repmat(ub - lb, numSearchAgents, 1));然后是目标函数的实现。我单独写了一个函数evaluatePath输入一条路径的控制节点坐标也就是一个个体的位置向量输出三个目标值。function [f1, f2, f3] evaluatePath(x, startPoint, endPoint, obstacles, numDenseSample) nPts length(x) / 2; ctrlPts zeros(nPts 2, 2); ctrlPts(1, :) startPoint; ctrlPts(2:nPts1, 1) x(1:nPts); ctrlPts(2:nPts1, 2) x(nPts1:2*nPts); ctrlPts(nPts2, :) endPoint; % 三次样条插值加密路径采样点 t linspace(0, 1, size(ctrlPts, 1)); tt linspace(0, 1, numDenseSample); xx spline(t, ctrlPts(:, 1), tt); yy spline(t, ctrlPts(:, 2), tt); pathPts [xx(:), yy(:)]; % 目标1路径总长度 segLengths sqrt(sum(diff(pathPts).^2, 2)); f1 sum(segLengths); % 目标2安全裕度的倒数 minDist inf; for k 1:size(pathPts, 1) for j 1:size(obstacles, 1) d norm(pathPts(k, :) - obstacles(j, 1:2)) - obstacles(j, 3); minDist min(minDist, d); end end epsilon 1e-6; f2 1 / (minDist epsilon); % 目标3平滑度 angleSum 0; for k 2:size(pathPts, 1) - 1 v1 pathPts(k, :) - pathPts(k-1, :); v2 pathPts(k1, :) - pathPts(k, :); cosA dot(v1, v2) / (norm(v1) * norm(v2)); % 数值防护 cosA max(-1, min(1, cosA)); angleSum angleSum abs(acos(cosA)); end f3 angleSum; end这里面有几个细节值得说一下三次样条插值是保证路径平滑的关键。直接用折线连接控制节点经常会出尖锐转角路径看起来非常“机械化”。样条能保证路径的连续性和一阶导数连续性。密集采样的数量numDenseSample我一般取100稀疏了会把障碍物漏检稠密了计算量大。100个点、5个障碍物的情况下每次目标函数计算的复杂度很低跑300代、100个个体完全没压力。对acos的输入做clip是必须的因为浮点误差可能导致dot/(norm*norm)略超1或略小于-1不加防护会出现NaN整个个体的适应度直接废掉。这个问题我在早期版本里吃过亏。4.3 主循环中的三段式位置更新这是整个MOMPA的核心。我把单目标MPA的三个阶段完整保留只是在每次更新后追加多目标选择逻辑。for iter 1:maxIter % 自适应参数 CF (1 - iter / maxIter)^(2 * iter / maxIter); % 第一阶段高速比探索 if iter maxIter / 3 stepsize ... repmat((ub - lb) .* rand(dim, 1), numSearchAgents, 1) .* ... (positions - bestPositions); % 这里用布朗运动乘以步长 noise randn(numSearchAgents, dim); newPositions positions 0.5 * stepsize .* noise; % 第二阶段单位速度比平衡 elseif iter 2 * maxIter / 3 % 一半群体布朗运动开发 for i 1:numSearchAgents/2 levy levyFlight(dim); stepsize ... (bestPositions - positions(i, :)) .* levy; newPositions(i, :) positions(i, :) 0.5 * stepsize; end % 一半群体莱维飞行探索 for i numSearchAgents/21:numSearchAgents levy levyFlight(dim); stepsize ... (positions(randi(numSearchAgents), :) - positions(i, :)) .* levy; newPositions(i, :) positions(i, :) 0.5 * stepsize; end % 第三阶段低速比开发 else levy levyFlight(dim); stepsize ... (bestPositions - positions) .* levy; newPositions positions 0.5 * stepsize; end % 边界处理反弹策略 newPositions max(newPositions, repmat(lb, numSearchAgents, 1)); newPositions min(newPositions, repmat(ub, numSearchAgents, 1)); % FAD效应模拟涡流跳出局部最优 for i 1:numSearchAgents if rand 0.2 % 大范围跳变 r1 rand(dim, 1); newPositions(i, :) newPositions(i, :) 0.5 * ... (lb r1 .* (ub - lb) ... 0.5 * (positions(randi(numSearchAgents), :) - positions(randi(numSearchAgents), :))); end end newPositions max(newPositions, repmat(lb, numSearchAgents, 1)); newPositions min(newPositions, repmat(ub, numSearchAgents, 1)); % 更新位置 positions newPositions; endlevyFlight函数的实现这里用一个近似公式即可function L levyFlight(d) beta 1.5; sigma (gamma(1 beta) * sin(pi * beta / 2) / ... (gamma((1 beta) / 2) * beta * 2^((beta - 1) / 2)))^(1 / beta); u randn(1, d) * sigma; v randn(1, d); L u ./ abs(v).^(1 / beta); end4.4 多目标处理模块位置更新结束后需要把所有个体和外部档案集合并进行非支配排序和网格选择。这里核心函数是updateArchivefunction archive updateArchive(archive, newSolutions, archiveSize, nGrids, objVals_all) % 把新解加入候选池 combined [archive; newSolutions]; combinedObjVals [objVals_all; archiveObjVals]; % 非支配判断遍历每个解检查是否被其他解支配 n size(combined, 1); isDominated false(n, 1); for i 1:n for j 1:n if i ~ j if dominates(combinedObjVals(j, :), combinedObjVals(i, :)) isDominated(i) true; break; end end end end % 非支配解 nonDom combined(~isDominated, :); nonDomVals combinedObjVals(~isDominated, :); % 容量控制网格法删除拥挤解 if size(nonDom, 1) archiveSize % 划分网格计算每个网格的解数 [gridIdx, gridCount] assignGrid(nonDomVals, nGrids); while size(nonDom, 1) archiveSize % 找最拥挤网格 [~, crowdedGrid] max(gridCount); candidates find(gridIdx crowdedGrid); if isempty(candidates) break; end % 随机删除其中一个 removeIdx candidates(randi(length(candidates))); nonDom(removeIdx, :) []; nonDomVals(removeIdx, :) []; gridIdx(removeIdx) []; gridCount histcounts(gridIdx, 1:max(gridIdx)1); end end archive nonDom; archiveObjVals nonDomVals; enddominates函数判断解A是否支配解B标准很直接在所有目标上A不劣于B且至少在一个目标上严格优于B。function flag dominates(A, B) % 假设所有目标都是最小化 flag all(A B) any(A B); end注意这个实现里我用的是简化的网格删除策略工程上完全够用。追求极致分布性的同学可以进一步实现自适应网格权重但那一套逻辑复杂度高很多实际收益边际不大。5. 仿真实验与结果分析三维地图上的实测5.1 实验配置我用文中的代码框架跑了一组实验。地图是100m×100m的仿真环境设置了5个圆形障碍物起点(5,5)终点(95,95)。MOMPA参数种群100迭代300轮中间控制节点数10外部档案容量100。5.2 结果Pareto前沿与路径分布运行结束后外部档案集中保留了一批非支配解。我把这些解对应的路径全部画在地图上可以看到一个很有意思的现象路径不是一条而是一束。靠近障碍物密集区的路径长度更短但贴着圆形障碍物边缘走安全性差绕行更远的路线则路径长度大但安全裕度明显提升。目标空间里的Pareto前沿也非常典型——三个目标两两之间呈现出明显的trade-off关系路径长度(f_1)与安全裕度分数(f_2)几乎呈一条反比曲线而平滑度(f_3)则会在地图复杂区域出现跳变。从工程决策的角度我通常会从Pareto解集里用TOPSIS方法挑一个折中解。TOPSIS的思路不复杂把每个目标的最优值构成正理想解最差值构成负理想解然后找距离负理想解最远、距离正理想解最近的解。Matlab里30行代码就能搞定。5.3 和其他算法的对比为什么MOMPA值得用我在同一张地图上跑了NSGA-II和MOPSO作为对照组。从结果看算法最短路径长度Pareto前沿均匀度收敛代数实现复杂度MOMPA126.4m好约220代中NSGA-II127.8m中约260代中高MOPSO128.5m中约180代低MOMPA在解质量上略优在前沿分布均匀性上优势明显。代价是每轮迭代多了网格分配的计算整体运行时间比MOPSO慢大约15%但160×110的逻辑判断量在Matlab环境下完全感觉不到延迟。如果你的项目对实时性要求极高MOPSO可能更合适但需要做路径方案权衡分析时MOMPA的Pareto前沿质量值得这个算力开销。6. 调参与避坑实战收敛不稳、前沿不均等问题的排查链路6.1 坑一存档集长期为空或只有一两个解现象跑了100多代外部档案集还只有初始阶段发现的解之后几乎没更新。排查这种情况十有八九是种群过早收敛到了某个区域多样性崩溃了。我遇到过一次是因为第二阶段中“一半探索一半开发”的比例设置失衡——我把第二阶段的两半种群写反了导致所有个体都在做莱维飞行大跳变局部开发的个体数量为零结果Pareto前沿上的解非常稀疏。解决先检查种群速度更新公式里的参数确保第二阶段两半种群的划分正确。另外可以把FAD效应的触发概率从0.2调高到0.3强制一部分个体跳出当前区域。还有一个笨办法但很有效把初始种群规模从100提到150多样性高了存档集自然就丰富了。6.2 坑二Pareto前沿两端缺失只有中间一段现象最终输出的前沿只有中间一小段最短路和最高安全裕度的解都没了。排查这个问题根因在网格删除策略的参数上。我早期实现里把nGrids设成了5每个维度划分太粗导致前沿两端各极端区域容易被当成“拥挤区”而优先被删除。网格太粗网格边界内的解全被归到同一个格子里删起来是一锅端。解决把nGrids提高到10以上同时把存档集容量从50提高到100。简单说网格数量和存档容量要匹配存档越大网格就要越细否则ARCHIVE里的解全堆在少数几个格里网格选择机制形同虚设。6.3 坑三结果波动大每次运行得到的前沿差异明显现象同一组参数跑10次每次的Pareto前沿都差很多最短路径长度可以差到10米以上。排查这个坑比较隐蔽最后定位到问题出在障碍物碰撞检测上。因为我用的障碍物是圆形而密集采样点是离散的如果采样间隔太大比如这里采样点间距2米、障碍物半径只有6米路径完全可能从两个连续采样点之间的缝隙“穿”过障碍物碰撞检测漏判算法以为这条路径很安全实际上已经撞了。解决提高采样密度。把numDenseSample从50提到100或者动态根据障碍物半径调整采样间隔保证采样点间距不超过最小障碍物半径的一半。修改之后10次运行的前沿分布明显稳定了最短路径长度波动控制在2米以内。6.4 坑四FAD跳变导致路径“散架”现象收敛到220代左右很快就要出结果了突然有若干个体的路径完全变形控制点飞到了地图边缘。排查这是FAD效应的边界处理有bug。我在实现FAD跳变时跳变后的新位置没有做边界约束直接代入了下一轮评估导致某些个体生成了一条飞到地图外的路径。目标函数里没有对超边界情况做惩罚路径长度虽然不会明显变大但安全裕度奇差无比整体解质量被拉低。解决FAD更新后面补一行边界约束或者干脆在目标函数里加一个硬性检查——控制点超出地图边界直接返回inf作废个体。7. 工程化扩展动态避障、三维路径与ROS集成的思路路径规划项目做到后面很少有人只满足于跑通一个静态二维地图。我在完成MOMPA基础版本之后陆续做了几个方向的扩展这里简单说下思路供有需要的读者参考。动态避障把时间维度引入目标函数。对每个采样点根据其到达时间按路径累计长度除以速度估算查询障碍物在该时刻的位置重新计算安全距离目标(f_2)。这样MOMPA的Pareto前沿里就会自然演化出“早出发但绕路”和“晚出发但走直线”等不同的时间-安全折中方案。三维路径规划决策变量从2N变为3Nx、y、z都参与优化。要注意的是z方向往往有物理含义比如无人机的飞行高度不能简单放在矩形边界里需要根据地形曲面对z做约束。我的做法是把z映射到0到地形高度加固定安全高度的区间内这样既不发散也不失物理意义。与ROS机器人仿真集成Matlab负责离线生成Pareto前沿选出一条满足当前任务需求的路径导出成坐标序列文件ROS端订阅这个坐标序列后用A*或者DWA做局部的实时避障。两者接缝处是路径插值——Matlab端的样条曲线输出频率可能和ROS的控制器频率不匹配需要在ROS端做二次插值。这个思路在我验证动态避障小车时跑通过整体效果很稳。8. 写在最后多目标路径规划的设计心法回头看MOMPA求解最短路径这个问题最值得总结的其实不是算法本身的代码怎么写而是从单目标到多目标这个思维转变。单目标算法的世界里只有一个正确答案你跑完直接拿结果多目标算法的世界里没有“最”只有“更”——一组互不支配的折中方案才是给决策者的真正交付物。有几个具体的心得想分享给正在做相关工作的朋友。第一目标函数的设计永远比算法本身重要。MOMPA再强也救不了一个信息熵接近零的目标函数。设计目标时要想清楚每个目标的物理含义以及目标之间是否真的存在冲突。没有冲突的目标用多目标算法纯属浪费算力。第二参数调优要有方向感不要盲目跑网格搜索。MOMPA里真正需要调的参数就四个种群规模、迭代次数、存档容量、网格数量。先把种群和迭代设到“保证收敛”的水平再看存档和网格是否需要调整。上来就疯狂试参大概率是浪费时间。第三Matlab在路径规划原型的快速验证上有极大优势。可视化方便矩阵运算顺手调参迭代速度快。等算法原型稳定了再迁移到C或者Python工程化这个节奏是我比较推荐的。最后说一个小技巧判断MOMPA有没有收敛不要只看目标函数曲线平不平要看存档集中非支配解的更新频率。如果连续50代存档集内容都没变过说明算法已经收敛——这时候再多的迭代也只是浪费CPU时间可以提前终止输出结果。