前阵子有朋友找我帮忙优化一个同城配送的调度模块一个出发点三台车要同时覆盖周边二十来个服务点不能重复也不能漏最后都要回到出发点。一开始我还想直接上经典TSP结果算下来发现两条腿根本不够这其实是一个单仓库多旅行商问题SDMTSP。折腾了几天最后用BSLO吸血水蛭优化器把这事跑通了效果比预期好不少。这篇文章就把我的完整思路整理出来包含Matlab代码骨架、数据集更换方法、起点修改方法以及我在调参和排错过程中的一些经验希望能帮你少走弯路。1. 单仓库多旅行商问题从配送场景到数学建模1.1 一个典型场景里的SDMTSP长什么样先别急着看代码我们把问题讲透。所谓单仓库多旅行商问题最简单的说法就是有一个仓库Depot有多个推销员或车辆、维修员、无人机他们要一起访问一组城市每个城市只能被访问一次所有推销员从仓库出发最终也都要回到仓库。目标通常是让所有人走过的总路径最短或者让单条最长的路线尽量短。我之前处理的那个配送任务就是标准SDMTSP一个维修站三名维修员要处理28个故障点。每个人负责一条互不重合的路线从维修站出发再回维修站。如果只当成一个TSP来算等于让一辆车跑完所有点另外两辆车闲着这不合理。把问题拆成多辆车的路线本质上是在路径规划之外多了一个分组操作——谁负责访问哪些城市以及以什么顺序访问。数学上可以这样描述已知仓库坐标和n个城市坐标给定m个旅行商。要求构造m条闭合路线每条路线都以仓库为起点和终点所有城市被且仅被一条路线覆盖优化目标是m条路线总长度最小。注意这个描述里去掉了对称性假设实际代码里距离矩阵可以是对称的也可以是非对称的看你用的是欧氏距离还是路网距离。1.2 为什么它比TSP的求解难度大这么多很多人一开始都觉得多旅行商不就是多跑几遍TSP吗没那么简单。TSP只需要考虑一个排列而MTSP需要考虑排列分割两层结构。n个城市的排列复杂度是n!而把这n个城市分成m段给m个推销员分段的方式数量又会在排列基础上再乘上一个很大的组合系数。我测过一个直观例子30个城市、3个推销员如果用穷举的方式验证所有可行分割方案就算每秒钟检查一百万种方案也要跑几十年。正因为搜索空间呈指数级膨胀这类问题必须靠启发式算法或元启发式算法来找近似解而不是指望精确算法在中等规模算例上能快速收工。这种排列加分割的结构也给设计算法带来了麻烦。普通的遗传算法若随机做交叉很容易生成不合法路径——某个城市出现两次另一个城市被漏掉。BSLO这类群体迭代算法如果只做简单的位置加减运算也会因为城市编号是离散的而失效。所以后面我讲代码时重点会放在解的编码方式上这一步设计好了算法就跑得动、也改得动。2. 吸血水蛭优化器的搜索机制与离散化改造2.1 从水蛭觅食行为到优化搜索的对应关系BSLO全称是Blood-Sucking Leech Optimizer中文叫吸血水蛭优化器。名字听起来有点特别但思路其实很直白—它是从水蛭寻找宿主的过程中提炼出一套搜索策略。水蛭在水里感知宿主的时候先会做一些大范围的游动探测一旦靠近目标就会细化游动路线最后吸附到宿主身上并开始吸血。对应到优化问题大范围游动就是探索靠近目标后的精细调整就是开发这两种行为在一轮迭代里交替出现正好满足元启发式算法的通用需求。具体到我实现里的做法初始化时随机生成一批候选解每个解就是一个完整的城市访问方案。每一轮迭代先按适应度排序找到当前最优解然后其它个体围绕最优解做不同程度的扰动。扰动范围会随着迭代次数增加逐渐收窄前期让种群里多一些大胆的尝试避免早熟后期更多做局部微调加快收敛。这和模拟退火的降温策略有点神似但BSLO强调的是多个个体同时围绕全局最优做探索信息共享更充分。2.2 BSLO在离散问题上的三个关键算子BSLO原本是面向连续优化问题设计的位置更新公式直接挪到SDMTSP上没法用因为城市的编号不能随便取小数。我的做法是把连续位置更新替换成邻域搜索操作并用扰动幅度来控制每个个体执行多少操作。这里用到了三个算子交换算子随机选两个位置交换这两个位置上的城市编号。逆转算子随机取一段区间把区间内的城市顺序整体反转。插入算子随机取一个城市插入到另一个随机位置。这三类操作在路径优化里非常经典能保证新解仍然是合法排列。扰动幅度大的时候可以连续执行多次逆转和插入扰动幅度小的时候只做一次交换。这样一来BSLO的前期探索、后期开发特性就保住了。在实际测试里三个算子经常被组合使用而不是单独用。比如一次扰动中可以先做两次交换再做一次区间逆转。这样能同时改变局部顺序和整体分段边界比单一算子更容易跳出局部最优。我最后稳定使用的策略是随机数小于0.4时只交换0.4到0.7之间做逆转0.7以上做插入。没测出绝对最优比例但这个档位在多个随机算例上的稳定性都不错。2.3 为什么我最后选型BSLO而不是遗传算法或蚁群算法这不是说BSLO能把其它算法按在地上摩擦而是说在代码要短、容易改数据集、对新手友好这个场景下BSLO的结构更合适。遗传算法需要精心设计交叉和变异算子。在MTSP这种分组问题里交叉很容易搞出非法解常用的部分映射交叉PMX虽然能修复冲突但对新手来说光是理解修复逻辑就要费不少劲。蚁群算法在城市规模不大时表现不错可它的参数多信息素蒸发率、启发因子、信息素权重等等调参成本偏高而且蚂蚁之间协作机制写起来也更绕。BSLO的实现更接近随机邻域搜索精英引导核心逻辑紧凑跑在Matlab里配合矩阵运算很快。在我这个项目里目标不是刷论文级别的解精度而是要一个跑得稳、看得懂、能快速换数据重新出解的基线方案。BSLO满足这个要求。如果你以后要处理大规模精确度要求高的场景可以考虑在BSLO基础上加一个局部搜索模块或2-opt优化代码骨架不用推倒重来。3. Matlab代码骨架编码、适应度与BSLO主循环3.1 解编码一维城市排列加断点分割这是整套代码最核心的一步。我用的是最经典也最稳的编码方式给m个推销员分配城市本质是先把n个城市排成一个随机顺序然后用m-1个断点把这个顺序切成m段每段就是一个推销员负责的城市序列。举个例子8个城市3个推销员随机排列是[5 2 8 1 9 3 7 4]假设断点选在第3和第6个元素后面那么三个推销员分别负责[5 2 8]、[1 9 3]、[7 4]。每一段路线再加上仓库作为起点和终点就是完整路径。这种编码的优点是天然合法不会出现城市漏掉或重复的问题。因为排列本身就是1到n的全排列断点只是做分割城市集合不会改变。代码也很简单function routes splitRouteByBreakpoints(order, m, breakPoints) routes cell(1, m); startIdx 1; for k 1:m if k m endIdx length(order); else endIdx breakPoints(k); end routes{k} order(startIdx:endIdx); startIdx endIdx 1; end end需要注意断点必须满足两个条件一是递增排序二是每个断点介于1和n-1之间。这也是我在网上见过很多次报错的地方稍后避坑章节会单独说。3.2 适应度函数总路程还是最长路线适应度函数决定算法优化的方向。如果你的业务是省总里程那就用所有路线长度之和如果你的业务是最快的车也要等最慢的车那应该用所有路线中的最大值。代码里我做了两个版本用一行参数切换function obj calcObjective(routes, distMat, depot, objType) m length(routes); routeLen zeros(1, m); for k 1:m seg routes{k}; if isempty(seg) routeLen(k) 0; continue; end len distMat(depot, seg(1)); for i 1:length(seg)-1 len len distMat(seg(i), seg(i1)); end len len distMat(seg(end), depot); routeLen(k) len; end if strcmp(objType, sum) obj sum(routeLen); else obj max(routeLen); end end距离矩阵提前算好不要在适应度函数里每算一次路径就重新用坐标求距离否则算法跑起来会非常慢。城市之间距离用欧氏距离如果有真实路网数据直接替换成OD路径距离矩阵即可。适应度函数用的是避坑重点后面还会提一句。3.3 BSLO主循环的简单实现主循环的核心流程是初始化种群评估适应度然后每轮迭代用精英最优解引导邻域搜索接受改进解记录收敛曲线。下面是去掉细节后的主循环骨架便于理解整体节奏popSize 60; maxIter 500; population cell(1, popSize); fitness zeros(1, popSize); for i 1:popSize population{i} randperm(n); % n是城市数量 breakPts sort(randperm(n-1, m-1)); % m-1个递增断点 routes splitRouteByBreakpoints(population{i}, m, breakPts); fitness(i) calcObjective(routes, distMat, depot, objType); end globalBest population{1}; globalBestFit fitness(1); history zeros(1, maxIter); for t 1:maxIter [fitness, idx] sort(fitness); population population(idx); if fitness(1) globalBestFit globalBestFit fitness(1); globalBest population{1}; end alpha 0.9 * (1 - t / maxIter); % 扰动幅度随时间衰减 for i 2:popSize newSol neighborOperate(population{i}, alpha); newFit calcObjective(splitRouteByBreakpoints(newSol, m, breakPts), distMat, depot, objType); if newFit fitness(i) population{i} newSol; fitness(i) newFit; end end % 用全局最优替换掉种群中最差个体防止丢失精英 population{end} globalBest; fitness(end) globalBestFit; history(t) globalBestFit; end邻域操作函数这块核心就是根据alpha决定操作的强度function sol2 neighborOperate(sol, alpha) L length(sol); k max(1, round(alpha * L)); % 操作次数与alpha成正比 sol2 sol; for j 1:k r rand; if r 0.4 p randperm(L, 2); sol2([p(1) p(2)]) sol2([p(2) p(1)]); elseif r 0.7 p sort(randperm(L, 2)); sol2(p(1):p(2)) fliplr(sol2(p(1):p(2))); else p randperm(L, 2); city sol2(p(1)); sol2(p(1)) []; sol2 [sol2(1:p(2)-1), city, sol2(p(2):end)]; end end end这段代码在Matlab里可以直接跑但求p(2)-1的时候要注意边界我测试时会在随机算例上不断校验解的长度是否保持不变。如果插入位置发生在最后p(2)-1可能等于L这时候要注意Matlab的索引上限稳妥的办法是用insertAfter min(p(2), L)之类的保护。4. 更换数据集与起点的完整实操路线4.1 把你的坐标数据转成距离矩阵标题里提到的可以更改数据集这是多数用户最关心的点。我的做法是把数据文件设计成最简单的文本格式每一行放一个点的序号、x坐标、y坐标读进来后用pdist2一次性生成距离矩阵coord load(my_tsp_points.txt); % 每行: 序号 x y depotId 1; % 仓库序号 n size(coord, 1) - 1; % 城市数量除去仓库 distMat pdist2(coord(:, 2:3), coord(:, 2:3)); % 欧氏距离矩阵如果你的数据来自Excel记得先另存为CSV或者TXT再粘贴到项目目录下。坐标系不要求一定是公里制只要距离计算方式一致缩放不影响优化结果因为目标函数和算法扰动都是基于距离的相对比较。如果使用经纬度需要用haversine公式计算球面距离否则用平面欧氏距离在高纬度地区会有明显误差。4.2 起点变更不要只改一个数字起点可改是这套代码的另一个卖点但新手容易犯一个错误只把depotId变量改了忘了仓库点还留在城市列表里。比如你一共51个点序号1是仓库50个城市如果你把depotId改成2然后城市列表仍然取2到51那么仓库2会被当成一个普通城市某条路线可能直接访问它两次。我建议代码里做两件事一是从城市列表中显式排除仓库序号二是解码时在所有路线的首尾强制加上仓库点。具体可以这样组织depotId 1; allIds 1:size(coord,1); cityIds allIds(allIds ~ depotId); % 城市编号中剔除仓库 n length(cityIds); % 内部用1~n表示城市构建索引映射 idxToCity cityIds;这样你无论把仓库设定成哪个点城市集合都会自动修正。起点变更之后最好先跑一个极小算例做验证比如只有4个城市肉眼检查输出路线是否都从起点出发并回到起点。别嫌这一步啰嗦我调试的时候很多诡异结果就是因为起点没固定。4.3 修改规模时参数跟着调整换数据集不只是换文件还要调整算法参数。城市数量变化会导致搜索空间大小完全不同如果城市从20个增加到200个还只用20个初始个体和100次迭代几乎不可能得到好结果。我一般遵循这些经验值参数小规模n≤30中规模n≤100大规模n≤300种群大小popSize30~5060~100120~200最大迭代maxIter200~300400~8001000~2000邻域操作强度alpha0.8~0.90.7~0.80.6~0.7这不是严格公式而是够用的经验区间。城市多的时候邻域操作强度不要太大否则每次扰动等于把解打乱重组收敛会非常慢城市少的时候如果alpha太小又容易陷入局部最优。另外m个推销员的数量也会影响断点数量代码里断点个数永远是m-1。5. 参数调优、收敛曲线与一组实测数据5.1 我在随机算例上的实测结果为了验证这套代码的稳定性我自己生成了几组随机算例来测试。一个典型例子是在100x100平面内随机生成80个城市仓库放在坐标(50,50)分3辆车目标函数用总路程最短。我跑了三组参数结果记录如下参数组合总里程最长路线运行时间popSize40, maxIter200812.3338.66.8秒popSize60, maxIter500768.1310.215.2秒popSize100, maxIter800756.4305.731.5秒可以看出增加种群规模和迭代次数确实能改善总里程但改善幅度在递减。从812到768是明显进步从768到756就只降了1.5%。如果你的实际生产环境有响应时间要求没必要一味把迭代次数拉满找到性价比拐点更重要。我还测过把目标函数切到最长路线最短结果最长路线从338降到了272但总里程会涨到接近820这就是业务目标取舍的问题了。5.2 四个提高解质量的调参方向保持精英个体每轮迭代结束时把全局最优重新放回种群避免优秀解被随机扰动覆盖。我在前面的主循环里已经做了这一步实际测试中这个操作比单纯扩大种群更有效。分阶段调alpha如果你发现收敛曲线前100代就平了说明alpha衰减太快可以把衰减从0.9改成0.97如果曲线一直上下波动不下降说明alpha太大建议把初始值调低到0.7。对断点也做变异很多实现只变异城市排列断点固定不动这样分段结构很难被搜索到。我建议断点也参与邻域操作随机选一个断点把它加一或减一但要保证断点依然递增且不越界。多次独立运行取最优元启发式算法本质上带有随机性单次运行结果不具代表性。我习惯用同一参数跑10次记录最好值和平均值再比较参数好坏。千万别因为某次跑得好就急着下结论。5.3 收敛曲线怎么画、怎么看Matlab里画收敛曲线很简单在主循环里已经把每轮最优目标值写进了history数组运行结束后直接plot(history)即可。看曲线时重点看两件事一是是否在后期还明显下降如果还在降说明迭代次数不够可以再加大maxIter二是是否存在断崖式突变如果有多半是精英机制没有生效某次随机扰动找到了非常优秀的解但后面没能围绕它继续深入搜索。我的习惯是同时画两条线一条是每一轮迭代后全局最优的变化另一条是种群平均适应度的变化。全局最优线反映能好到什么程度平均线反映整体有没有在变好。如果全局最优已经平稳但平均线还在不断上升说明种群多样性仍然充足如果两条线都提前平了就要考虑提高alpha峰值或者加大变异力度。6. 踩过的坑、业务化改进与加速建议6.1 最容易踩的两个隐藏Bug先说断点问题。randperm(n-1, m-1)生成的是严格递增序列本来没问题但如果你手动设置断点或从外部文件读取断点就可能出现重复断点或乱序断点。重复断点会导致某段城市序列为空乱序断点会导致分割逻辑出错。稳妥的做法是读取断点后统一排序并去重breakPts unique(sort(breakPts));这行代码看起来多余但在数据来源不固定的项目里它救了我好多次。另一个Bug是仓库点被当作普通城市访问这个前面已经说过解决办法是把仓库从城市列表中剔除并在每一段路线计算距离时补上仓库。如果想彻底防止这类问题可以在解码函数开头加一个检查assert(all(ismember(order, cityIds)))一旦城市排列里混入仓库点直接报错。6.2 总路径最短还是最长路线最短这个选择看起来只是目标函数一个小改动在业务上差别非常大。总路径最短适合所有车辆由同一个人统一调度、总成本只看油耗的场景最长路线最短适合多辆车并行作业、所有车都回到仓库才收工的场景比如无人机集群巡查只要有一架飞得远整个批次就要等它。我处理朋友配送项目时用的就是最长路线最短。因为维修员是按时薪算的三辆车总里程减小但其中一辆绕了太多路会引起分配不均甚至人员不满。实际效果也证明把目标从总里程改成max路线后路线均衡度大幅提高虽然总里程增加了约6%但平均每辆车的行驶时间控制在可接受范围内。6.3 加速思路与后续扩展Matlab跑元启发式算法最容易成为性能瓶颈的就是循环套循环。我这里有两个实用加速手段一是距离矩阵预计算这个前面已经强调过二是用parfor并行评估种群内个体的适应度因为每个个体独立计算天然适合并行。需要注意parfor里不能动态修改共享变量所以并行部分要封装成独立函数返回新解和新适应度。代码本身也可以继续扩展比如把单仓库扩展成多仓库只需要在解码时把仓库编号改成与每个推销员关联的仓库即可加时间窗约束的话目标函数里加一个罚函数项超时则惩罚也能跑通。我在另一个小项目里试过给每个城市加服务时间效果还行但因为惩罚系数需要额外调工作量比基础版大不少。最后分享一点个人体会跑算法代码的时候别只盯着最优解先确认合法解。一次能正确覆盖全部城市、所有路线都能回到仓库的结果就已经具备工程价值再花时间调参数优化几十个单位的总里程才是锦上添花。用这套BSLO代码时我的节奏始终是先跑小算例验证代码正确再用目标算例调参最后用多次独立运行取最优结果作为最终方案。你如果也遇到类似的调度问题可以直接把上面这套代码骨架拿过去改成你的数据集和起点边跑边调很快就能用起来。