做物流调度的朋友应该都遇到过这个挺尴尬的现场一条线路上每个客户既要送货又要取货车辆容量一会儿腾空一块一会儿又被旧货填满。凭经验挑出来的“最短路径”可能跑到一半就因为装不下新收的货而卡住。这就是运筹学里常说的同时取送货车辆路径问题VRPSPD。我之前在做一个城市末端配送加回收的仿真项目时用模拟退火算法在Matlab里实现了完整求解器代码从数据组织到邻域搜索全部踩了一遍坑。这篇文章不适合翻课本而是直接把能跑的代码、参数设置逻辑、以及那些文档里不会写的坑位都摊开讲正在做物流优化算法、写学术论文或者做调度系统开发的同学可以直接拿去用。1. 同时取送货车辆路径问题到底是什么1.1 从经典VRP到VRPSPD多了一个动态容量经典的车容量受限车辆路径问题CVRP模型很清爽所有客户都只要求送货车从车场出发把货一件件卸掉车辆装载量只减不增所以容量约束是一条单调递减的曲线判断一条路线是否超载非常简单——从头加到尾任何时刻货量低于容量上限就行。VRPSPD把“取货”这个动作加了进来。每个客户身上挂两个需求一个是要送到他手里的货量delivery一个是要从他手里带走并运回车场的货量pickup。车辆到达一个客户后先卸下送货再装上取货。于是车辆装载量就不是单调的了卸货让负载下降装货又让负载上升可能卸得多、装得少也可能这个网点收上来的旧货比派下去的新货还多。这种“忽高忽低”的动态容量会让很多传统VRP的直观经验失效。举例来说一条路线 A-B-C 距离最短但B点取货量巨大车辆到B之后直接超载而绕一点走A-C-B总路程多几百米但装载曲线全程都在容量限制以内。传统VRP算法如果只盯着总距离很可能给出前者。VRPSPD的真正难点不只是搜索空间大更在于“路线距离”和“装载可行性”这两个目标或约束是强耦合的任何一个邻域改动都必须重新评估整条线路的装载过程。1.2 容量约束的数学表达假设车辆从车场出发负责访问客户序列 S (s_1, s_2, ..., s_m)车场编号为0。本辆车在其路径上负责配送的客户送货量总和记为Q0 Σ_{i1..m} delivery(s_i)车辆从车场出发时装载量就是Q0而不是全局所有客户的送货量总和——这一条特别容易写错。访问第k个客户后先卸货delivery再装货pickup车辆装载量为L_k Q0 - Σ_{i1..k} delivery(s_i) Σ_{i1..k} pickup(s_i)容量约束要求每一步都满足0 ≤ L_k ≤ Q注意L_k并不是单调序列所以不能用“只检查最大值”这种偷懒方式必须逐点检查。这也是VRPSPD评估函数比CVRP慢一点的根源不过对几十个客户规模来说逐点检查的开销完全可以忽略。1.3 现实中哪些场景在用VRPSPD在物流行业几乎无处不在我以前整理过一张表这些场景看起来千差万别底层模型完全一样应用场景送货方取货方典型行业快递末端派送新包裹送到客户客户退换货包裹收走电商物流、快递网点饮料分销新品饮料配送给门店空瓶和旧货回收快消品供应链共享单车调度把车辆运到热点区域把损坏车辆运回维修点出行运营制造业逆向物流新零部件配送旧件、包装载具回收汽车/电子制造图书流转系统新书配送到分馆读者归还的旧书集中回收图书馆、校园服务如果站在行业趋势上看电商退货率一直在往上走“送新收旧”一车跑完既省油费又能减少空驶这是逆向物流和循环经济背景下最典型的优化场景。所以不要觉得VRPSPD只是学术圈在玩的理论题它本质上是很多调度系统里“既要派件又要收件”那段业务逻辑的数学抽象。2. 为什么用模拟退火思路与参数选型2.1 邻域搜索与Metropolis准则VRPSPD是个NP-hard问题客户量稍微一大指望精确算法在可接受时间内求出最优解就不太现实。行业里常用做法是元启发式算法遗传算法GA、禁忌搜索TS、模拟退火SA、变邻域搜索VNS都能做。我个人偏向SA原因很朴素实现逻辑最短不需要维护禁忌表也没有遗传算法那群“种群、交叉、变异”参数要调它在局部搜索的基础上加了一个温度控制的概率接受机制差解在一定概率下会被接受从而跳出局部最优。Metropolis准则是整个算法的核心决策规则。假设当前解为x邻域扰动后得到新解x目标变化量为 Δ f(x) - f(x)如果 Δ 0说明新解更优无条件接受如果 Δ ≥ 0说明新解更差也并非直接拒绝而是以概率 exp(-Δ/T) 接受其中T是当前温度。温度T高的时候exp(-Δ/T)接近于1算法几乎不区分好坏解相当于在整个解空间里大范围乱逛温度慢慢降低后接受差解的概率越来越小算法逐步从“探索”过渡到“精细开采”最后稳定在某个优质解附近。这个“先高温撒网、后低温收网”的过程很像金属退火也是算法名字的由来。2.2 初始温度、降温系数、迭代次数怎么定参数选择是模拟退火最容易劝退新手的地方但本质上也就四个初始温度T0、降温系数alpha、每个温度下的迭代次数L、终止温度Tend。我给出一套实际项目中常用的经验值参数经验取值说明初始温度T0不手设抽样估计随机生成100个邻域解统计目标增量绝对值均值meanD取T0 5~10 * meanD降温系数alpha0.95 ~ 0.99越大搜索越充分但越慢常用0.98内循环次数L100 ~ 500或 n*20客户少取小值客户多取大值终止温度Tend0.01 ~ 1温度低于它是退出太小只增加后期无效迭代T0为什么用抽样而不是手设因为不同算例的目标函数量级差异很大坐标是100公里范围还是10公里范围总距离完全不一样。一个算例的“差解增量”可能是50另一个可能是5000直接用固定T0100根本不合适。我现在的习惯是任何一个新算例进来都先跑一遍这个抽样估计让初始接受率保持在0.7~0.9之间这能直接避免“初始温度太高浪费时间”或者“一上来就冻住”两个极端。2.3 多车辆编码怎么用一条序列表达整个车队VRPSPD默认有一个车队好几辆车同时工作。要在模拟退火中处理多车辆最方便的是采用“带0分隔符的访问序列”编码。举个例子[0, 5, 2, 0, 7, 1, 3, 0]这个序列表示两辆车第一辆访问客户5和2第二辆访问客户7、1和30是车场。车辆数上限由序列里0分隔符的个数决定。这种编码的好处是模拟退火的邻域操作可以同时在“车内的客户重排”和“车与车之间的客户交换”两个层面上进行不需要额外设计复杂的拆分合并算子。初始解可以直接用随机排列加随机切分先把所有客户随机打乱再随机选择切分点插入0。由于邻域算子只操作客户节点序列中的0分隔符数量通常保持不变这意味着车队规模基本稳定。当然如果翻转区间包含了0车辆分配结构会被重新洗牌这是允许的只要最终序列里没有语法错误即可。3. Matlab完整实现与核心代码拆解3.1 数据准备与距离矩阵代码里我统一用以下变量组织数据coord所有坐标点第一行是车场后面n行是客户demandD1×n向量每个客户的送货量demandP1×n向量每个客户的取货量vehicleCap单辆车容量numVehicle车队最大车辆数。把车场放在坐标第一行有个好处Matlab里pdist2计算距离矩阵后车场编号就是1客户编号从2到n1访问序列中的客户编号直接用2到n10在序列里作为分隔符存在评估时用dist(1, seg(1))表示从车场出发到第一个客户的距离清晰且不容易错。% 示例数据生成n30个客户 rng(2025); n 30; coord [50, 50; rand(n, 2) * 100]; % 车场在(50,50) demandD randi([1, 12], 1, n); % 每个客户送货量 demandP randi([1, 8], 1, n); % 每个客户取货量 vehicleCap 80; numVehicle 6; dist pdist2(coord, coord);3.2 核心评估函数分段计算装载量评估函数是整套算法的心脏。它接收一条带0分隔符的序列把序列拆成若干辆车访问的客户段对每一段独立计算行驶距离和容量惩罚。重点在于每辆车的初始装载量是本段客户送货量之和而不是全局送货总量。function totalCost evaluateSeq(seq, demandD, demandP, dist, vehicleCap) totalCost 0; i 1; seqLen length(seq); penaltyFactor 50 * mean(dist(:)); % 超载惩罚系数 while i seqLen if seq(i) 0 i i 1; continue; end % 收集从当前客户到下一个0之前的客户段 seg []; while i seqLen seq(i) ~ 0 seg(end 1) seq(i); i i 1; end if isempty(seg) continue; end % 车辆初始装载本段全部送货量之和 load sum(demandD(seg)); if load vehicleCap totalCost totalCost (load - vehicleCap) * penaltyFactor; end % 从车场到段内第一个客户 totalCost totalCost dist(1, seg(1)); load load - demandD(seg(1)) demandP(seg(1)); if load vehicleCap totalCost totalCost (load - vehicleCap) * penaltyFactor; end % 段内后续客户 for k 2:length(seg) totalCost totalCost dist(seg(k - 1), seg(k)); load load - demandD(seg(k)) demandP(seg(k)); if load vehicleCap totalCost totalCost (load - vehicleCap) * penaltyFactor; end end % 最后从段末客户回车场 totalCost totalCost dist(seg(end), 1); end end这段代码里有几个细节值得注意。第一个是load的计算更新时间点我把“减去送货量、加上取货量”放在距离累加之后表示车辆到达当前客户并完成服务后的装载状态即离开该客户时的实际载量。容量检查就在这个时间点做因为这是车辆在当前客户点上的最大载货时刻。第二个是空段处理如果序列里连续出现两个0说明这辆车没有被使用评估函数直接跳过不产生距离、也不产生惩罚这给了算法“关闭车辆”的能力。3.3 邻域算子交换、插入与2-optSA搜索质量高度依赖邻域算子。我这套实现里安排三个算子按概率混合使用交换算子Swap随机选两个客户交换它们在序列中的位置。这是最基础的扰动方式适合微调。插入算子Relocate把一个客户从当前位置取出插入到另一个客户之后。这个算子能直接改变一辆车的客户组成非常实用。2-opt反转随机选择序列中的一段区间整体反转。因为可能跨越0所以它经常造成“跨车辆重排”相当于一次较大的结构调整。三个算子的比例我习惯用0.4 / 0.4 / 0.2也就是交换和插入各占四成2-opt占两成。如果跑出来的结果老是陷入局部最优把2-opt比例调高到0.3左右会产生更剧烈的结构变化如果收敛一直不稳定优先提高交换和插入比例。function newSeq neighbor(seq) customerIdx find(seq 0); if length(customerIdx) 2 newSeq seq; return; end r rand(); if r 0.4 % 交换两个客户 iPos customerIdx(randi(length(customerIdx))); jPos customerIdx(randi(length(customerIdx))); while jPos iPos jPos customerIdx(randi(length(customerIdx))); end newSeq seq; newSeq([iPos, jPos]) newSeq([jPos, iPos]); elseif r 0.8 % 把客户i插入到客户j之后 iPos customerIdx(randi(length(customerIdx))); jPos customerIdx(randi(length(customerIdx))); while jPos iPos jPos customerIdx(randi(length(customerIdx))); end val seq(iPos); newSeq seq; newSeq(iPos) []; if jPos iPos insertAt jPos - 1; else insertAt jPos; end newSeq [newSeq(1:insertAt), val, newSeq(insertAt 1:end)]; else % 2-opt反转一段区间 a randi(length(seq)); b randi(length(seq)); if a b tmp a; a b; b tmp; end newSeq seq; newSeq(a:b) seq(b:-1:a); end end插入算子的索引处理容易写错。注意原序列里我先把iPos位置的值取出来删掉然后根据jPos是否在原序列iPos后面来重新计算插入位置。如果不处理这个索引偏移删除元素之后直接插到jPos大概率会插错位置导致客户重复或丢失。3.4 主循环与温度调度主循环结构不复杂但有一个很容易忽略的细节在进入降温循环前最好先做一次初始温度抽样估计。我给出的实现里如果调用者没有指定T0就从初始解出发生成100个随机邻域解统计目标函数增量均值的绝对值然后用10倍均值作为初始温度。这种做法比任何经验公式都更能适配具体算例。function [bestSeq, bestCost, rec] sa_vrpspd(demandD, demandP, coord, vehicleCap, numVehicle, params) n size(coord, 1) - 1; dist pdist2(coord, coord); if nargin 6 params struct(); end if ~isfield(params, alpha), params.alpha 0.98; end if ~isfield(params, Tend), params.Tend 0.01; end if ~isfield(params, L), params.L 300; end % 初始解随机客户排列 随机切分为numVehicle段 custOrder randperm(n); initSeq 0; if numVehicle 1 cutPositions sort(randperm(n - 1, numVehicle - 1)); for i 1:n initSeq [initSeq, custOrder(i)]; if ismember(i, cutPositions) initSeq [initSeq, 0]; end end else initSeq [0, custOrder]; end initSeq [initSeq, 0]; currentCost evaluateSeq(initSeq, demandD, demandP, dist, vehicleCap); currentSeq initSeq; bestSeq initSeq; bestCost currentCost; % 自动估计初始温度 if ~isfield(params, T0) sampleDeltas zeros(1, 100); baseSeq initSeq; for s 1:100 newSeq neighbor(baseSeq); sampleDeltas(s) abs(evaluateSeq(newSeq, demandD, demandP, dist, vehicleCap) - currentCost); end params.T0 10 * mean(sampleDeltas); end T params.T0; rec []; while T params.Tend for k 1:params.L newSeq neighbor(currentSeq); newCost evaluateSeq(newSeq, demandD, demandP, dist, vehicleCap); delta newCost - currentCost; if delta 0 || rand() exp(-delta / T) currentSeq newSeq; currentCost newCost; if currentCost bestCost bestSeq currentSeq; bestCost currentCost; end end end rec [rec, bestCost]; T T * params.alpha; end end注意我在每轮温度迭代结束后把当前最优成本bestCost记录到rec数组方便事后画收敛曲线。如果是正式项目rec里还可以记录当前温度值用来观察算法在哪个阶段完成了大部分改进。3.5 完整可运行算例把上面函数存成sa_vrpspd.m再新建一个脚本运行下面的代码就能看到结果和收敛曲线rng(2025); n 30; coord [50, 50; rand(n, 2) * 100]; demandD randi([1, 12], 1, n); demandP randi([1, 8], 1, n); vehicleCap 80; numVehicle 6; params.alpha 0.98; params.Tend 0.01; params.L 300; [bestSeq, bestCost, rec] sa_vrpspd(demandD, demandP, coord, vehicleCap, numVehicle, params); disp([最优总距离: , num2str(bestCost)]); disp([最优序列: , num2str(bestSeq)]); plot(rec, LineWidth, 1.5); xlabel(降温轮次); ylabel(当前最优总距离); title(模拟退火收敛曲线);我实测过这个数据组合一般运行几十秒内能收敛最优序列对应的总距离相比初始随机解通常能改进20%~30%。当然我的随机算例具体数字每次不完全一样但这套代码作为原型工具完全够用。4. 常见问题与排查技巧实录4.1 收敛太快结果明显不是最优如果算法在rec曲线上前几十轮就基本不动了大概率是初始温度太低导致算法根本没机会接受差解一上来就陷入某个局部区域“冻住”。解决办法是把T0调大。如果你不手设T0就检查一下抽样估计的代码是否漏了abs绝对值或者样本数太少导致均值失真。另外一个常见原因是alpha太小比如取到0.85甚至0.8温度下降呈指数级衰减没过几十轮就成了低温随机爬坡全局搜索能力被削减。我强烈建议alpha不要低于0.95宁可多跑几分钟也别让它提前“下班”。4.2 结果路线的超载一直消不掉这个现象出现时先别急着改算法重点检查惩罚因子。我的评估函数里罚因子取的是50倍平均距离如果你的算例平均距离特别大或重量单位很大这个系数可能不足以把不可行解“推开”。一个粗暴但有效的做法把penaltyFactor临时改成1e6如果超载解还是频繁出现说明问题根本不在罚函数而在于车辆数给得太少或者单个客户的取货量本身大于容量。还有一类隐蔽问题初始解里某些段的送货量总和已经超载。比如一辆车被随机分到好几个送货量很大的客户初始装载就是超的。SA在这种解的基础上要“绕很大一圈”才能把客户重新分配到其他车辆有时会转不出来。针对这种情况可以在初始解生成时做一个简单检查如果某段sum(demandD)vehicleCap就在切分点附近做一次局部重排。哪怕不做这种精修复也可以通过增加2-opt的占比来增强跨车辆调整能力。4.3 每次运行结果差异很大SA在本质上是随机搜索算法不同的随机种子对应不同结果这是正常现象。但如果你发现同一个算例跑十次最好最差差了20%以上那说明当前参数下算法稳定性不足。我的处理思路是固定随机种子保证项目交付时结果可复现多运行几次取总距离最小的解作为最终输出增加“低温精化”阶段在主循环结束后对bestSeq做一轮小步长局部搜索比如把每个客户的相邻位置都尝试一遍插入接受所有不劣化总距离的改进。这一步实现简单经常能再白捡2%~3%的优化效果。4.4 关于2-opt的一个隐蔽坑2-opt反转区间如果横跨了0分隔符会直接把一条路线里的客户段和分隔符顺序全部打乱。表面看序列仍然合法评估也能正常跑但它会引入一种“快速跨车辆重排”的副作用。这种剧烈变动在低温阶段出现太多会让收敛曲线剧烈抖动。如果发现低温阶段结果不稳定一个限制策略是把2-opt的区间限制在同一个0分隔符段内也就是只反转同一辆车内部的客户序列这样就不会把其他车辆的结构搅乱。具体实现时先找到区间内最接近的两个0把区间裁剪到那段客户范围内即可。这是我实际踩过坑之后的经验。5. 实测效果与进一步扩展5.1 小规模算例实测表现我用上面那组30客户数据做过一个简单对比随机初始解的总距离通常在450~550之间模拟退火跑完一遍最优总距离能稳定在350~400这个区间改进幅度在20%到30%之间。更关键的是最终解里没有出现任何超载段——每辆车每个时段的装载量都被压制在了容量以内。这说明罚函数Metropolis准则的组合在约束处理上足够有效。如果加大客户规模到50甚至100个运行时间会线性增长但总距离改进幅度不会明显缩小。SA的一个特点是问题规模越大初始解越差它能带来的相对改进越大。所以它特别适合作为快速原型帮你在不引入商业求解器的前提下先把业务方案的合理边界摸出来。5.2 可以继续升级的方向最直接的升级是把SA精化后的结果交给2-opt局部搜索做最后打磨。这基本不花什么开发时间却能在现有结果上再提高几个百分点。另一个方向是替换搜索策略比如用变邻域搜索VNS替代SA在多套邻域结构之间动态切换通常能在相同时间预算下找到更优解。如果项目对确定性有要求也可以用Gurobi等商业求解器求解小规模算例再用SA的解作为大规模算例的近似上界形成一个“精确启发式”混合框架。我在算法之外还有一个实际用的技巧不只用总距离作为目标而是在目标函数里叠加一个“车辆使用数量小惩罚”比如每启用一辆车加0.1倍平均单条线路距离。这样算法会在总距离相近的情况下倾向于压缩车辆数量从而降低人员成本。这个改进在物流调度项目里非常实用尤其当车队数量本身是可优化变量的时候。5.3 建议直接拿到项目里用的坑位提醒如果你要把这套代码接入实际调度系统请一定注意数据口径。VRPSPD里的delivery和pickup必须在同一个容量单位下要么都是件数要么都是重量不能一个车型用立方米、一个用吨就硬塞进同一个容量约束。这类问题在真实工程项目里遇到频率最高却最容易被代码演示忽略。另外真实路网里的距离通常不是欧氏距离建议用地图API先算好真实的OD距离矩阵把dist输入换成实际路网距离算法本身不需要任何修改只要距离矩阵合理结果就有参考价值。时间窗、装卸时长这类额外约束也可以继续在评估函数里以惩罚项方式扩展不需要推翻整体框架。我个人在实际操作中最深的体会是模拟退火这类元启发式算法真正的成败往往不在算子设计而在参数校准和惩罚函数量级匹配。同样的代码数据量级不同T0和penaltyFactor都可能要重调。所以拿到一个新算例我会先做一轮“不加约束、只看距离”的试跑把目标量级和邻域增量搞清楚再决定惩罚参数怎么设。这比一上来就调alpha、改算子有效率得多。希望这套带代码的完整实现能让你少走几步弯路。