1. 为什么多旅行商问题比TSP难一个量级1.1 从一个外卖配送场景说起一个人送30个外卖点和三组人各送一批外卖点听起来像是同一件事的简单升级但实际做起来你会发现前者只需要回答“按什么顺序走”后者还要回答“谁送哪些点、每组分多少、顺序又怎么排”而这几个问题是互相绑定的。这就是多旅行商问题MTSP和经典旅行商问题TSP之间真正的差距。这篇文章我会用灰狼优化算法GWO求解MTSP给出完整的Matlab实现思路并分享我实测过程中踩过的几个关键的坑。适合正在做路径规划、组合优化相关毕业设计或者想把群体智能算法落到实际问题中的读者。换到工程场景里MTSP太常见了。快递网点有若干个配送员仓库有一队AGV小车或者一架无人机要给多个区域投放物资背后都是同一个问题有一批目标任务点、若干台执行设备怎么安排设备和路线让整体成本最低、负载相对均衡。很多人一上来就说“这不就是每个设备各跑一个TSP吗”等真正建模的时候才发现分组方案和组内顺序是强耦合的先分组再独立排序出来的解往往离最优解差得很远。1.2 MTSP的数学建模分组与排序的耦合难题先把MTSP用数学语言写清楚。设0号节点为配送中心depot1到n为需要服务的城市旅行商数量为m。每个旅行商负责一条闭合路径起点和终点都是配送中心。需要同时决策两件事把n个城市划分成m个集合每个集合对应一个旅行商在每个集合内部确定访问顺序。目标函数通常是最小化所有路径长度之和min Z Σ L(R_k)k 1, 2, ..., m约束条件包括每个城市恰好属于一条路径每条路径从depot出发并回到depot每个旅行商至少负责一个城市不允许出现空路线。理解这个问题的难处关键是要看到分组和排序不是独立的两步。假设你已经把城市分成了三组那每组内部确实可以当成TSP来解但你凭什么说这个分组是好的换个分法哪怕组内顺序不是最优总里程也可能大幅下降。反过来如果组内排序靠前又能反过来让某个本来很差的划分变得可用。这种分组与排序的组合爆炸使得哪怕n只有30、m只有3暴力枚举的空间也大到完全没法直接搜。所以这类问题适合用元启发式算法去逼近最优解而不是指望穷举或者精确算法。实际工程里大家也普遍接受“足够好”的解因为路径规划本身的动态扰动、道路拥堵等因素已经让“数学最优”失去了绝对意义。1.3 灰狼优化算法凭什么适合这类问题近几年做组合优化的人经常在TSP、VRP相关文献里看到灰狼优化算法Grey Wolf OptimizerGWO它由Mirjalili等人于2014年提出属于群体智能算法家族。选它来解MTSP我个人有几个很实际的理由不需要梯度信息。GWO只要求“给一个解能算出它的适应度值”这对于离散组合优化太重要了因为路径长度这类目标往往不可导、不连续。参数极少。GA要操心交叉概率、变异概率、选择压力PSO要调惯性权重和学习因子GWO核心只需要控制一个递减参数a对新手极其友好。有明确的领导机制。狼群向头狼、副头狼、第三头狼学习这种“少量精英引导”的结构既保证了收敛速度又保留了搜索多样性。大量文献已经验证了GWO在TSP、VRP类问题上的可行性但多数公开代码只做到TSP层面。MTSP需要在编码上额外做文章正好是本文要展开的重点。当然GWO不是银弹它在离散问题上容易出现后期收敛停滞这一点我放在第五章专门讲。先理解算法本身的运行机制后面才能知道怎么改、怎么补。2. 灰狼优化算法的核心机制模拟狼群狩猎的背后逻辑2.1 种群等级与角色分工GWO的思想很直观把每个候选解当成一只狼整个种群里有严格的等级结构。适应度最好的个体称为alpha头狼它掌握当前最优路线信息适应度第二、第三的个体分别是beta和delta其余个体统称为omega。更新过程中每一只omega个体不是漫无目的地随机游走而是同时参考alpha、beta、delta三个领导者的位置来调整自己。这样做的好处是搜索初期三只“头狼”可能差距很大普通狼群分布范围广搜索多样性有保证随着迭代推进头狼们逐渐统一到相近区域整个种群也会跟着收拢收敛速度就有了保障。相比于只朝全局最优个体靠拢三头狼并行的设计能减少过早陷入单一局部最优的概率。2.2 包围、追捕与进攻的数学表达GWO的位置更新可以拆成三部分来看。第一部分是搜索距离的估计。对任意一只狼X它会计算自己和“头狼”Xp之间的距离D |C · Xp(t) - X(t)|这里的C 2 · r2r2是[0,1]之间的随机数相当于给距离判断加了一层随机抖动。引入C的意义在于如果所有狼都严格计算直线距离算法会太机械容易同步聚集到同一个点。加了抖动之后狼与狼之间天然有了差异化搜索范围更宽。第二部分是收敛因子A。A 2a·r1 - a其中a从2线性递减到0。当|A| 1时狼群会远离当前领导者这对应全局搜索阶段当|A| 1时狼群向领导者位置逼近这对应局部开发阶段。通俗理解A就像一个油门前期数值大、冲得远后期数值小、收得稳。第三部分是综合更新。每只狼分别按alpha、beta、delta三个领导者计算一次候选位置然后取平均X1 Xα - A1·|C1·Xα - X| X2 Xβ - A2·|C2·Xβ - X| X3 Xδ - A3·|C3·Xδ - X| X(t1) (X1 X2 X3) / 3三条路径既考虑了当前最优方向又通过A和C引入随机性本质上是对“跟随精英”和“自主探索”做了一次加权折中。很多初学GWO的人只盯着公式本身忽略了a的时序设计。a线性递减意味着算法前半程以探索为主后半程以开发为主这也是为什么它和TSP这类多峰问题能匹配起来。2.3 从连续搜索到离散排列编码是第一道关卡GWO的所有公式都是为连续变量设计的。狼的位置X是一个实数向量但在MTSP里解是一组离散的城市序列和分组方案不可能直接把城市编号拿去做加减乘除。比如城市编号5和7做差得到2这个2没有实际意义经过一轮位置更新后个体还会莫名其妙出现重复编号、漏掉某些城市的情况。所以用GWO解MTSP必须设计一层“编码-解码”映射GWO负责在连续实数空间里迭代解出来的每个实数向量通过解码函数转成一条合法的MTSP路径方案。通俗地说GWO只负责“出主意”解码函数负责把“主意”翻译成能算长度的路线。这层翻译设计得好不好直接决定了算法能不能收敛到好解。这里给一个常见反例有人为了让GWO适配TSP用城市编号直接作为个体位置更新后取整结果每迭代几代解就非法了因为出现了重复城市和丢失城市。本质上就是少了这层映射。MTSP比TSP更麻烦的地方在于它不仅要把n个城市排成序列还要在序列里插入分组边界所以编码设计上要有两个信息维度而不仅仅是城市顺序。3. Matlab实现细节编码设计、适应度函数与主循环3.1 优先级-分割双段编码设计我在这篇文章里采用一种直观且工程上可靠的编码方式个体位置是一个长度为nm-1的实数向量前半段和后半段各管一件事。前n维城市的访问优先级。每个城市对应一个实数数值越小表示越优先被访问排序后得到一条完整的城市访问顺序。后m-1维子路径的分配信息。它决定刚才那条长序列从哪里切开分成m条子路径每条子路径对应一个旅行商。用具体例子走一遍。假设n6个城市m2个旅行商个体位置前半段排序后得到的城市顺序是[4, 2, 3, 6, 1, 5]后半段只有1个数值假设是0.42。n-m4round(0.42×4)2于是第一条路径分到213个城市第二条路径分到剩余4-213个城市这里的基本逻辑是先把每个旅行商保底分配1个城市再按后半段数值分配剩余的城市配额。得到两条路径旅行商1depot → 4 → 2 → 3 → depot旅行商2depot → 6 → 1 → 5 → depot可以看到这个编码天然保证了每个旅行商至少有一个城市因为保底机制先铺好了。排序解决“顺序”后半段解决“分组”两个信息全部揉在一个连续实数向量里GWO可以直接在[0,1]范围内更新位置更新完之后再解码整个过程合法且高效。3.2 解码函数与距离计算有了编码思路Matlab代码就好写了。核心函数是decode它接收一个个体位置、城市数、旅行商数量和距离矩阵返回每条路径的城市顺序、总距离和最长子路径距离。function [routes, totalDist, maxDist] decode(pos, nCities, m, distMat) % 输入 % pos - 个体位置长度 nCities m - 1 % nCities - 需要服务的城市数量 % m - 旅行商数量 % distMat - 距离矩阵包含配送中心和所有城市 % 输出 % routes - cell数组每项是一条子路径含首尾depot % totalDist - 所有子路径长度之和 % maxDist - 最长子路径长度 nRemain nCities - m; % 1. 排序得到城市访问顺序城市编号1到nCities [~, order] sort(pos(1:nCities)); % 2. 根据后 m-1 维计算分割点 if m 1 cumY round(pos(nCities1:end) * nRemain); cumY sort(cumY); % 防止出现重复分割点导致某条子路径长度为0 for j 2:m-1 if cumY(j) cumY(j-1) cumY(j) cumY(j-1) 1; end end cumY min(cumY, nRemain); splits [0, cumY, nRemain]; else splits [0, nRemain]; end % 3. 每个旅行商分到的城市数 segLens diff(splits) 1; routes cell(m, 1); totalDist 0; maxDist 0; startIdx 1; % 4. 按分配长度切分城市序列并计算各子路径长度 for k 1:m endIdx startIdx segLens(k) - 1; citySeq order(startIdx:endIdx); % 配送中心编号1城市编号整体加1保证索引从1开始 route [1, citySeq 1, 1]; routes{k} route; legDist 0; for j 1:length(route)-1 legDist legDist distMat(route(j), route(j1)); end totalDist totalDist legDist; maxDist max(maxDist, legDist); startIdx endIdx 1; end end这里有几个细节值得说明。城市编号加1是为了符合Matlab从1开始索引的习惯配送中心固定为1。分割点cumY的作用范围限制在[0, nRemain]保证最多把“剩余城市配额”分完不会出现城市越界。去重循环是很多人容易漏掉的一步我在第五章会专门讲它引发的坑。3.3 GWO主循环的代码骨架解码函数准备好之后GWO主循环非常直白核心过程就是初始化、评估、选三头狼、循环更新位置。% 主参数设置 nCities 30; % 城市数 m 3; % 旅行商数 NPop 60; % 种群规模 MaxIt 500; % 最大迭代次数 lambda 1.0; % 最长子路径惩罚权重 % 随机生成城市坐标配送中心固定 rng(42); cityXY rand(nCities, 2) * 100; coords [50, 50; cityXY]; dx coords(:,1) - coords(:,1); dy coords(:,2) - coords(:,2); distMat sqrt(dx.^2 dy.^2); % 初始化种群 nVar nCities m - 1; lb zeros(1, nVar); ub ones(1, nVar); Positions rand(NPop, nVar); fitness zeros(NPop, 1); for i 1:NPop [~, totalDist, maxDist] decode(Positions(i,:), nCities, m, distMat); fitness(i) totalDist lambda * maxDist; end % 初始化三头狼 [~, idx] sort(fitness); Alpha Positions(idx(1), :); AlphaScore fitness(idx(1)); Beta Positions(idx(2), :); BetaScore fitness(idx(2)); Delta Positions(idx(3), :); DeltaScore fitness(idx(3)); bestCost inf; history zeros(MaxIt, 1); % 主循环 for it 1:MaxIt a 2 - it * (2 / MaxIt); for i 1:NPop for d 1:nVar % 向alpha学习 r1 rand(); r2 rand(); A1 2*a*r1 - a; C1 2*r2; Dalpha abs(C1 * Alpha(d) - Positions(i,d)); X1 Alpha(d) - A1 * Dalpha; % 向beta学习 r1 rand(); r2 rand(); A2 2*a*r1 - a; C2 2*r2; Dbeta abs(C2 * Beta(d) - Positions(i,d)); X2 Beta(d) - A2 * Dbeta; % 向delta学习 r1 rand(); r2 rand(); A3 2*a*r1 - a; C3 2*r2; Ddelta abs(C3 * Delta(d) - Positions(i,d)); X3 Delta(d) - A3 * Ddelta; NewPos(i, d) (X1 X2 X3) / 3; end end % 边界处理 NewPos max(NewPos, lb); NewPos min(NewPos, ub); Positions NewPos; % 重新评估 for i 1:NPop [~, totalDist, maxDist] decode(Positions(i,:), nCities, m, distMat); fitness(i) totalDist lambda * maxDist; end [~, idx] sort(fitness); if fitness(idx(1)) AlphaScore Alpha Positions(idx(1), :); AlphaScore fitness(idx(1)); end if fitness(idx(2)) BetaScore Beta Positions(idx(2), :); BetaScore fitness(idx(2)); end if fitness(idx(3)) DeltaScore Delta Positions(idx(3), :); DeltaScore fitness(idx(3)); end history(it) AlphaScore; end这段代码结构很标准每个维度单独处理随机数三头狼分别引导最后平均得到新位置。实测下来这个版本在没有额外优化的前提下跑30城3商人的规模几百次迭代轻松收敛速度也很快。3.4 目标函数里的“总里程vs均衡性”权衡在MTSP里目标函数不是一刀切的。如果你运营的是自有车队关注总里程更合理因为油耗、电耗和总里程正相关如果你是平台调度把任务分给多个独立承运人那还要考虑司机之间的工作量是否公平这时候“最长的子路径”会比“总路径”更能反映服务质量。我在代码里把两个指标都算了出来并用一个权重lambda做折中fitness totalDist lambda * maxDistlambda0时只优化总里程lambda越大系统越倾向于把所有子路径压得差不多长。这个设计非常实用因为实际需求往往不是纯单目标而是希望在总成本和负载均衡之间找一个可接受的点。你甚至可以把lambda设为连续变量多跑几次看趋势再拍板用哪组解。4. 实验设计与结果分析30个城市、3个旅行商的测试4.1 测试场景与参数设置为了验证代码效果我设计了下面这个测试场景配送中心在平面正中心30个城市点坐标在[0,100]区间内均匀随机生成3个旅行商负责全部服务。固定随机种子保证任何人都能复现同一组数据。参数设置值城市数量 n30配送中心坐标(50, 50)城市坐标范围[0, 100] × [0, 100]旅行商数量 m3种群规模 NPop60最大迭代次数 MaxIt500惩罚权重 lambda0 或 1对比两种模式如果用Matlab没有统计工具箱pdist2可能没法用所以我在代码里直接用广播方式计算欧氏距离矩阵一行dx、dy就搞定不依赖任何额外工具箱。4.2 收敛曲线与路径规划结果按照上面的设置跑完用下面的代码可以绘制收敛曲线和路径图。% 收敛曲线 figure; plot(history, LineWidth, 1.5); xlabel(迭代次数); ylabel(适应度值); title(灰狼优化算法求解MTSP收敛曲线); grid on; % 路径图 [routes, totalDist, maxDist] decode(Alpha, nCities, m, distMat); figure; hold on; colors lines(m); plot(coords(1,1), coords(1,2), kp, MarkerSize, 12, MarkerFaceColor, k); for k 1:m route routes{k}; plot(coords(route,1), coords(route,2), o-, Color, colors(k,:), LineWidth, 1.2); end legend(配送中心, 旅行商1, 旅行商2, 旅行商3); xlabel(X坐标); ylabel(Y坐标); title(多旅行商路径规划结果); hold off;收敛曲线的典型走势是前100代左右下降非常快曲线近乎垂直100到300代之间进入缓慢下降阶段300代之后基本平缓偶尔有小幅波动但不改大局。这正是GWO“前期探索、后期开发”特性的直观体现。路径图上三条颜色不同的闭合回路都从配送中心出发各自覆盖一片区域整体上不会出现大面积交叉。我在一次固定种子下的实测结果大致如下。目标模式总距离最长子路径最短子路径lambda0只优化总里程512.6247.8131.9lambda1加入均衡惩罚589.4210.3176.8这个结果很有代表性。lambda0时算法把所有力气都用在全盘总距离上确实省了里程但代价是某一位旅行商的路线几乎是另一位的一半负载极度不均。lambda1时总距离上升了大约15%但最长子路径从247.8降到210.3最短子路径从131.9升到176.8整体均衡性明显改善。具体选择哪个权重完全取决于实际调度场景。4.3 怎么判断算法有没有跑对新手最容易遇到的问题不是“代码报错”而是“代码不报错但结果明显不合理”。我建议用两个手段自查。第一个是退化测试。把m设成1此时MTSP退化成标准TSPdecode里后m-1维字符串为空算法实际上就在用GWO解TSP。如果你的编码和主循环是正确的退化后的结果应该和独立TSP求解结果处于同一数量级不会出现明显的路径乱飞。这个测试还能同时验证解码函数在不同m下是否健壮。第二个是小时算例手工核对。把城市数设为4、旅行商数设为2手动排列组合一下就能算出最优解是多少然后看GWO能不能稳定找到这个解。如果连这种小规模都找不到那大概率不是算法问题而是解码函数或者适应度计算写错了比如距离矩阵下标错位、路线起点终点没闭合这类低级错误。5. 实测踩坑记录从“算法能跑”到“结果能用”之间隔着什么5.1 分割点重复导致空路线问题第一次写完代码我把问题规模设成n20、m4跑起来结果某次运行后解码出来的routes里出现了一个只有depot、depot两个元素的路径也就是某个旅行商一个城市都没分到。我排查了很久最后发现根因在后半段分割点的生成上。设想n20、m4nRemain16。如果后3维映射后得到的cumY是[3, 3, 10]排序之后仍是[3, 3, 10]那么第一个和第二个分割点重合中间那条子路径分配到的城市数量是0。虽然城市优先级排序没问题但split直接把一个合法方案切成了一条非法方案。这种问题在迭代中会随机出现不做处理的话算法偶尔会“偷偷丢掉”一个旅行商。解决办法就是在decode里加一个严格递增检查遇到重复或倒退的分割点就顺延一位并最终用min限制上限。代码就是第三章里展示的那几行。不要小看这个小修复它直接决定了解空间里有没有合法解。5.2 只优化总距离的隐性陷阱很多人会把TSP的思维直接搬到MTSP觉得目标函数写总距离最短就行了。但MTSP场景里如果你只写totalDist算法会倾向于把尽可能多的城市塞给同一辆“车”因为少派一辆车往往意味着更少的总行驶距离。结果是某个旅行商的路径奇长无比另外两个旅行商几乎在“散步”。我建议的工程做法是在目标函数里始终保留maxDist项。哪怕你的最终目标就是总里程最小也建议先跑一个lambda1的版本观察两者差异再做权衡。调度系统如果完全不管负载均衡最终上线时会引发很现实的员工满意度问题这一点客户不会在需求文档里写出来但验收时会提。5.3 收敛过早问题的对策标准GWO的一个通病是后期收敛过快种群迅速集中到某个局部最优附近之后无论怎么迭代位置更新都只是在小范围内扰动解的质量提升很有限。我在30城3商人的测试里也遇到了这个问题基本上到250代左右曲线就“黏住”了。可以尝试把a的递减方式从线性改成非线性让早期探索阶段更长一些a 2 * (1 - (it / MaxIt)^2);这样a在前期下降得慢狼群保持更长时间的多样化搜索后期再快速收拢用于精细开发。实测这种方式在部分随机坐标集上能多挤出5%到10%的目标改进。另一个更简单粗暴但有效的策略是种群重置。如果连续50次迭代最优值都没有变化就把当前种群中除了三头狼之外的大部分个体随机重新初始化相当于让算法跳出局部最优再跑一轮。这种带重启机制的GWO在组合优化里的稳定性会明显好于原始版本。5.4 随机种子与多次运行的工程化习惯群体智能算法有个绕不开的事实单次运行的结果带有随机性你在自己的机器上跑出来的最优路径换一个随机种子可能就差不少。尤其是MTSP这种解空间特别大的问题单次运行很难保证你得到了一个稳定的好解。我的习惯是写一个外层循环连续跑20次每次都记录最优适应度和对应路径最后统一比较bestHistory zeros(20, 1); for run 1:20 % 重新初始化种群并运行GWO bestHistory(run) AlphaScore; end [bestValue, bestRun] min(bestHistory);然后取bestRun对应的路径作为最终结果。如果时间允许这个循环放在任何真实项目里都值得做。它不改变算法本身但能显著提升最终交付解的质量和可信度。6. 扩展思路多仓库、容量约束与算法混合6.1 多仓库MTSP与容量约束的编码扩展本文的代码针对的是单仓库场景所有旅行商从同一个depot出发。实际项目里更常见的可能是多仓库比如一个城市有多个配送分站每个分站有自己的车辆车辆完成任务后回到自己的分站。这个改动其实不用动GWO主循环只需要在解码阶段增加一段“depot分配”信息。个体编码可以扩展成nm-1m段前n段表示城市优先级中间m-1段表示分割最后m段表示每个旅行商从哪个仓库出发。decode时先读仓库分配再按原有逻辑生成路径距离矩阵也换成对应仓库到城市间的距离。容量约束CVRP的改法就更简单了。每个城市有需求量每辆车有容量上限在decode切分子路径时维护一个累计需求量一旦超过容量就强制结束当前子路径、开启下一条。这种“约束放在解码层”的写法保持了GWO搜索空间的连续性是工程上很实用的扩展思路。6.2 灰狼2-opt局部搜索的组合策略标准GWO在MTSP上能快速得到一个“还不错的解”但精细度有限原因是它没有一个针对局部路径结构的优化手段。2-opt是TSP里最经典的局部搜索操作在一段路径里选两条边如果反向连接能缩短路径就翻转中间的一段路网。function route apply2opt(route, distMat) improved true; while improved improved false; for i 2:length(route)-2 for k i1:length(route)-1 oldDist distMat(route(i-1), route(i)) distMat(route(k), route(k1)); newDist distMat(route(i-1), route(k)) distMat(route(i), route(k1)); if newDist oldDist - 1e-12 route(i:k) route(k:-1:i); improved true; end end end end end注意route的第一个和最后一个元素都是depot不能翻转它们所以索引从2开始、到倒数第二个结束。把2-opt嵌进GWO循环时我建议不用每代都做而是每隔20代对当前最优个体做一次局部优化这样计算代价可控解的质量提升非常明显。这个组合策略在文献里通常被称为GWO-2opt算是求解路径类问题的一个性价比很高的改进方向。6.3 同一套解码可以移植到GA、PSO最后说一个我实际做项目才体会到的点真正可复用的部分不是GWO本身而是decode这层。你用GA去求解核心要解决的是交叉和变异怎么作用在这个nm-1长度的编码上你用PSO求解核心也只是把位置更新公式换成速度-位移公式。一旦理解了编码和解码的逻辑整套框架迁移到其他群体智能算法大约只需要几十分钟。所以我一直觉得MTSP这类问题的学习路径应该是“先啃编码设计再玩算法”。很多同学一上来就钻研GWO的数学推导反而忽略了工程上最关键的映射层最后代码跑起来总觉得哪里不对劲。希望这篇文章能把这条弯路帮你绕过去。