1. 项目概述这是一个什么问题先聊一个做配送调度的朋友几乎都遇到过的情况快递员早上出门装满一车货沿途把包裹送出去同时还要把用户退的旧件、要寄出的箱子一件件收回来。车上的货物量不是在消耗而是在动态变化——送一件少一箱取一件又多一箱路径安排不好要么爆舱要么空车跑冤枉路。这个场景就是“同时取送货的车辆路径问题”英文简称 VRPSPDVehicle Routing Problem with Simultaneous Pickup and Delivery。这个题目在物流优化里属于命特别硬的那种学术上有完整的数学模型工程上有持续的实际需求。你能在共享单车调度、饮料瓶回收、退货上门取件、商超配送一体化的场景里看到它的影子。传统车辆路径问题VRP只考虑“从仓库送出去”VPRPD 又往往假设“先送完再取”或者“可以拆成两段”而 VRPSPD 的核心难点在于每个客户点既产生送货需求又产生取货需求车辆必须在同一次访问中同时完成这直接改变了整条路径上的负载变化规律问题复杂度显著上升。管理这类项目时我第一步做的就是明确需求边界是纯学术研究要出模型出算法还是真实业务系统要跑出可落地的调度方案。这个项目标题里明确写了“附 Matlab 代码”说明它偏向算法研究加工程验证的定位目标人群是物流工程、工业工程、运筹学方向的学生和做路径规划算法的工程师。这篇文章想把整个问题拆开讲清楚如何建模、为什么选模拟退火、Matlab 代码是怎么组织的、实际调试时坑都在哪里。2. 模型建立把现实约束翻译成数学语言2.1 核心符号与决策变量建模之前先假设一个很标准的场景单一配送中心车队是同质的每辆车容量一致每个客户点只被一辆车服务一次车辆最后必须回仓库。这个场景是所有 VRPSPD 研究的基础范式虽然真实业务可能更复杂但先把基础模型跑通后续加时间窗、多中心、异质车队都只是在它上面打补丁。模型里需要定义清楚这些要素客户点集合 C {1, 2, ..., n}仓库编号为 0每辆车容量为 Q客户点 i 有两个需求送货量 d_i 和取货量 p_i决策变量 x_ijk当车辆 k 从节点 i 行驶到 j 时取 1否则取 0决策变量 L_ik表示车辆 k 到达节点 i 时的当前已载货量目标函数通常是最小化总行驶距离因为距离直接对应油耗、时间、司机工时。如果要算成本也可以按距离乘一个单位运输成本系数但本质上是一样的。2.2 三类约束的物理含义第一类约束是“访问约束”每个客户点有且仅有一辆车进来有且仅有一辆车离开而且进出必须是同一辆车。这个约束保证服务不拆单就是实际业务里的“一车到底”。第二类约束是“容量约束”这是 VRPSPD 和传统 VRP 最本质的区别。传统 VRP 只需要检查车辆在离开仓库时载货量不超过容量因为沿途货只会减少。但 VRPSPD 里车辆经过一个客户点可能是先卸后装也可能是先装后卸路径上的载货量是波动的必须对路径上每个节点逐一检查到达节点 i 时的载货量减去 d_i 加上 p_i 之后仍然不超过容量也不能为负。这个约束写进代码里就是在插入或交换节点时都要重新模拟一遍完整路径上的负载变化。第三类约束是“路径连续性和仓库归宿约束”。车辆从仓库出发经过一串客户点后回到仓库中间不能断开。这个约束在编码时通常通过路径的表示方式自然满足不需要额外写复杂的检查逻辑。为什么约束要这么严格地建模因为很多初学者写算法时只在目标函数上花心思约束却处理得模模糊糊用惩罚项代替硬约束。惩罚项的系数难调过大过小都会得到不合理的解。我的建议是容量约束这种物理上必须满足的用硬约束去拒绝非法解时间窗这种软性约束才适合用惩罚函数处理。VRPSPD 里容量是红线不能碰。2.3 目标函数之外还要考虑什么如果只有距离目标算法很容易偏好一个“几何上很近但实际运营很别扭”的方案。比如所有客户点的取货量接近满仓一辆车跑不了几个点就必须返回仓库这时候强行安排超长路径虽然距离算起来短但实际根本装不下。所以在测试代码时我会同时输出距离指标和“车辆利用指标”——每辆车的最大载货峰值与容量的比值。这个指标能直观反映路径规划是否合理。当容量约束很紧时最优方案基本由容量决定当容量约束很宽松时问题退化为标准 VRP算法收敛难度大减。理解这个特性对后续调参很有帮助。3. 模拟退火算法设计为什么选它以及怎么设计才不迷路3.1 为什么选模拟退火而不是遗传或粒子群VRPSPD 已经被证明是 NP-hard 问题客户点一多精确算法根本撑不住。启发式算法里遗传算法、粒子群、蚁群都各有拥趸但我自己实践下来在 VRPSPD 这个场景里模拟退火有三个不可替代的好处。首先它的实现门槛最低。整个算法核心只有“生成新解、计算增量、按概率接受”三个动作没有遗传算法的交叉变异编码设计没有粒子群的参数调节负担。对于一个研究生或者刚转行做算法的工程师来说SA 是把问题跑通的最短路径。其次它在小规模和中规模算例上的稳定性极好。实际测试下来在 100 个客户点以内SA 配合好的邻域算子多次独立运行的标准差能控制在很低的水平。遗传算法容易早熟收敛粒子群容易陷入局部最优而 SA 的 Metropolis 接受准则天然具备“突跳”能力。第三SA 对目标函数的结构不敏感。VRPSPD 的目标函数和约束检查之间不是单纯的线性叠加负载变化导致可行性判断是高度非线性的。SA 只需要知道“新解比旧解好还是差”不需要利用梯度信息或者问题结构信息这让它非常鲁棒。3.2 解的编码一条染色体讲清楚车辆分配解的表达方式是整个算法里最需要想清楚的一环。我用的是“断点式客户序列编码”一条长度为 n 的序列表示客户的访问顺序序列中车辆负责的区间用分隔符切开。比如 10 个客户、3 辆车时一条编码可以是 3-7-1-0-9-2-5-8-0-4-6其中出现两次 0 表示两条分隔线意思是车 1 服务客户 3 7 1车 2 服务客户 9 2 5 8车 3 服务客户 4 6。这样一条一维数组就完整表达了路径方案方便操作和编码交换。这套编码有几个关键特性需要注意车辆数不能死板地固定有的客户点取货量大可能要多分配一辆车否则路径不可行序列里不允许出现重复的客户编号每个客户必须出现且仅出现一次。初始化时我会生成一个包含全部客户编号的随机排列再根据容量检查插入分隔符确保初始解一定是可行解。绝大多数 SA 代码跑出无法解释的乱结果都是因为初始解直接用一个非法解起手后面所有的收敛都建立在错误的起点上。3.3 邻域算子三条基本操作和一条秘密武器邻域结构决定了算法的搜索效率这里我重点讲四个算子。交换算子Swap随机挑选两个位置交换它们的客户编号。这是最基础的操作优点是改动小、扰动温和适合精细搜索缺点是当问题规模大时单次迭代的探索范围有限。两个客户点相距越远交换产生的路径变化越剧烈容易制造非法容量解。插入算子Insert从序列中抽取一个客户插入到另一个随机位置。这其实等价于改变客户点在某辆车或者某条路径中的服务顺序对路径结构的调整比交换更明显。插入算子对容量约束的冲击通常比交换小因为只是挪动一个节点负载分布的变化相对局部。2-opt 算子随机选择序列中的一段子路径将这段子路径整体反转。这是路径类问题的经典操作主要用于消除路径中的交叉和回头路。对 VRPSPD 来说2-opt 的效果需要特别小心因为整段反转后车辆到达每个节点时的负载状态会被重置原来满足容量约束的路径可能变得不满足。秘密武器是“容量偏差重定位算子”找出违反容量约束或接近容量上限的路径片段把其中取货量较大的节点重新分配给另一个当前负载较低的车。这个算子是我在调车次数多之后琢磨出来的。如果你用基本算子生成邻域的时候发现大量解不可行问题十有八九出在算法没有专门处理容量失衡。加一个定位重分配的动作收敛速度会明显加快。3.4 退火参数初温、终温、降温函数怎么给模拟退火的参数直接决定成败我给出经过验证的默认值供你起步参考。初温 T0 的设置目标是让初始阶段的接受率维持在 0.8 到 0.95 之间。可以让算法先随机跑几百步统计目标函数变化的平均幅度 delta_avg然后设 T0 10 * delta_avg或者 T0 -delta_avg / log(0.85)。实际测下来100 个点以内的问题初温在一个比较宽的范围内都不敏感不必太焦虑。降温函数我常用的是标准几何降温T_{k1} alpha * T_kalpha 取 0.90 到 0.99 之间。alpha 越小降温越快但解的质量明显下降alpha 接近 0.99 时能逼近较好解但迭代次数会非常巨大。我的建议是先用 alpha 0.95 起步看到收敛曲线的变化规律后再调。终止温度 T_end 的设置看你对解质量的要求。一般取 T_end 0.01 左右已经足够因为温度再低时接受劣解的概率已经极低算法退化成爬山。此外每个温度下的迭代次数内循环长度一般取客户数 n 的 5 到 10 倍。我自己用 100 个客户点时内循环取 1000 次外层温度迭代 200 次总迭代 20 万次Matlab 跑下来大约 40 秒一轮结果可以接受。4. Matlab 核心代码实现主循环、算子与容量校验4.1 主程序框架完整代码会贴在文末这里先讲清楚整体组织结构。一个清晰可维护的 SA 求解器至少需要分成四个函数生成初始解、计算路径总距离、生成邻域解、判断容量可行性。如果全部塞在主脚本里调试的时候会非常痛苦尤其是检查那容量到底是在哪里溢出时乱成一团的脚本会让你想砸电脑。我先把主循环的逻辑用伪代码形式描述。function [best_route, best_cost] sa_vrpspd(dist, demand, pickup, cap, params) % 初始化 current gen_initial_solution(dist, demand, pickup, cap); current_cost calc_total_cost(current, dist); best current; best_cost current_cost; T params.T0; while T params.T_end for iter 1:params.inner_iters % 生成邻域解 new generate_neighbor(current); % 容量可行性检查 if ~is_feasible(new, demand, pickup, cap) continue; end new_cost calc_total_cost(new, dist); % Metropolis 接受准则 delta new_cost - current_cost; if delta 0 || rand() exp(-delta / T) current new; current_cost new_cost; if new_cost best_cost best new; best_cost new_cost; end end end T params.alpha * T; end end这段代码的要点在于容量检查放在费用计算之前。为什么因为费用计算在所有节点之间做一段完整的距离求和复杂度是 O(n)而容量检查一旦发现某个节点负载溢出可以直接提前退出省下后续计算。成本敏感的项目里这两个操作的顺序会影响整体运行时间。4.2 载货量动态检查的实现VRPSPD 的容量检查和普通 VRP 最大的差异在于路径上负载的变化是“到达节点时先判断当前负载和送货量的关系再计算取货后的新负载”。我经常见到有人写成new_load load pickup(i) - demand(i)如果 new_load cap 就判为不可行。这种做法漏掉了关键的一步——必须先确保 load demand(i)否则车辆到达时车上根本没有足够的货可以卸。正确写法如下。function feasible is_feasible(route, demand, pickup, cap) load 0; for i 1:length(route) node route(i); % 必须先卸货确保当前能卸 if load demand(node) feasible false; return; end % 卸货 load load - demand(node); % 再装货 load load pickup(node); % 检查是否超容量 if load cap feasible false; return; end end feasible true; end注意这里默认服务顺序是“先卸后取”。实际业务中如果你允许“先取后卸”那就相当于取完的货一直在车上占空间对容量约束更严格。这个服务顺序的假设必须在代码注释里写清楚否则换数据跑的时候结果会有微妙差异还很难排查。4.3 距离计算和路径费用函数距离计算相对简单但是要小心坐标单位不统一的问题。很多公开数据集提供的坐标是经纬度直接算欧氏距离会得到非常离谱的结果。要么先转成平面坐标要么用 haversine 公式求球面距离。我这边为了演示方便直接用平面坐标算欧氏距离真实项目中这个函数要按数据情况替换。function total_cost calc_total_cost(route, dist) total_cost 0; % 从仓库出发 prev 1; % 假设仓库编号为 1 for i 1:length(route) total_cost total_cost dist(prev, route(i)); prev route(i); end % 返回仓库 total_cost total_cost dist(prev, 1); end仓库编号在代码里统一为 1客户编号从 2 到 n1。路径序列里只包含客户编号。这个设计统一了对 dist 矩阵的索引不至于出现“仓库编号到底是 0 还是 1”的糊涂账。4.4 邻域算子实现的细节这部分最容易出错的地方在于生成新解时别把序列弄出重复节点或遗漏节点。我提供的生成函数里每次操作前先复制完整序列操作完后再做一个查重校验虽然多花了一点时间但能避免很多莫名其妙的 bug。function new_route generate_neighbor(route) n length(route); new_route route; r rand(); if r 0.4 % Swap idx1 randi(n); idx2 randi(n); while idx2 idx1 idx2 randi(n); end new_route([idx1, idx2]) new_route([idx2, idx1]); elseif r 0.7 % Insert idx1 randi(n); node new_route(idx1); new_route(idx1) []; idx2 randi(n); if idx2 idx1 idx2 idx2 - 1; end new_route [new_route(1:idx2), node, new_route(idx21:end)]; else % 2-opt idx1 randi(n); idx2 randi(n); while idx2 idx1 idx2 randi(n); end if idx1 idx2 temp idx1; idx1 idx2; idx2 temp; end new_route(idx1:idx2) new_route(idx2:-1:idx1); end end插入算子这部分有个特别容易踩的坑从序列中删除一个元素后第二个随机位置要用更新后的序列长度重新计算下标。有人会因为忘记调整索引导致插入位置超出范围或者索引错位生成非法序列。我花了很长时间才发现自己代码里这个 bug排错的过程让人印象深刻。4.5 容量偏差重定位算子前面提到这是优化关键手段我把它的实现也放出来。它的思路是找到整个解中负载峰值的车辆路径把其中一个取货量大的节点移走放入另一辆负载较轻的车。function new_route capacity_shift(route, demand, pickup, cap) % 将路径按车辆拆分成多个子段 segments split_routes(route); % 找到负载最高的路径 max_load -1; max_idx -1; for i 1:length(segments) l peak_load(segments{i}, demand, pickup); if l max_load max_load l; max_idx i; end end % 从该路径中选出取货量最大的节点 [~, shift_pos] max(arrayfun((x) pickup(x), segments{max_idx})); node segments{max_idx}(shift_pos); % 删除该节点 segments{max_idx}(shift_pos) []; % 插入到负载最低路径的随机位置 min_load inf; min_idx -1; for i 1:length(segments) l peak_load(segments{i}, demand, pickup); if l min_load min_load l; min_idx i; end end insert_pos randi(length(segments{min_idx}) 1); segments{min_idx} [segments{min_idx}(1:insert_pos-1), node, ... segments{min_idx}(insert_pos:end)]; new_route merge_segments(segments); end这套逻辑看起来不复杂但运行效率足够高。实际项目中它在后段迭代温度较低时的作用特别明显温度低Metropolis 接受率低常规算子很难跳出局部可行域而重定位算子能定向调整容量失衡给解注入新的可行结构。5. 实验结果与参数敏感性分析5.1 基础数据集测试我用一个 20 个客户点的标准算例做了验证。坐标在 [0, 100] 的平面上随机生成取货量和送货量分别在 0 到 20 之间随机生成车辆容量设为 50。模拟退火参数取 T0 100alpha 0.95T_end 0.1内循环 200 次。总迭代约 3 万次Matlab 运行时间约 5 秒。我把收敛曲线画出来之后能看到非常明显的三个阶段前 5% 的迭代目标函数从 640 快速降到 520 附近这个阶段主要是初温高Metropolis 接受准则大量接受劣解算法在全局范围做粗糙探索中间 80% 的迭代曲线缓慢下降从 520 降到 460 左右每次温度降低后局部搜索开始发力最后 5% 的迭代曲线几乎不再变化最终停在 458。这三阶段的变化趋势很关键。如果发现曲线在第一阶段下降后就长时间不动说明初温太高或者 alpha 太大算法浪费了大量时间在高温振荡如果曲线还没有完全走平就降到接近终止温度说明 alpha 太小降温太快解的质量不够好。通过观察收敛曲线的形态反推参数问题是我调参时最重要的手段。5.2 参数敏感性测试实录我进行了一组对比测试把 alpha 从 0.90 改到 0.98其他参数不变最终解的距离分别如下。alpha最终距离迭代次数运行时间(秒)解的改进程度0.904823.1万3.5基线0.954585.0万5.25.0%0.9844712.6万14.87.3%从结果看alpha 从 0.95 调高到 0.98解只提升了 2.3 个百分点但运行时间翻了近三倍。这说明在工程落地中“用更长的计算时间去换 2% 的解质量”不一定划算。实际业务中调度系统往往有实时性要求与其让算法跑到最优不如在 5 秒内给出一个 98% 优的方案。这也是为什么我不建议盲目追求“高精度收敛”。另外一个容易被忽略的参数是随机数种子。模拟退火是随机算法同样的参数不同随机种子可能得到不同结果。我在验证算法稳定性时会固定几个种子做多次重复试验评估解的均值和方差。如果方差偏大说明算法的初始化或邻域算子设计不够稳健这时应调整的是算子的比例而不是反复调温度。5.3 场景扩展容量约束极端情况我额外测了一组容量很紧的极端情况取货量总和接近车辆总容量的 1.5 倍意味着车辆整体上必须多跑几趟路径规划必须偏向重调度。这组实验下算法是否容易陷入局部最优是最考验容量重定位算子的场景。实测发现单纯靠 swap、insert、2-opt 的组合很难让某些解从容量瓶颈中逃脱。加入容量重定位算子后解的质量迅速提升 10% 左右。这组对照充分说明了“邻域算子组合”比“在通用算子上盲目加迭代次数”更有效。6. 21 个实际调试中常见的问题排查速查表做这个项目的过程中我在 Matlab 代码调试上踩了不少坑有条理地整理成速查表方便你在复现的时候逐项排查。现象可能原因解决方法程序一运行就报索引越界插入算子删除节点后未更新序列长度每次删除后重新获取序列长度再选插入位置结果始终不收敛初温设置过低接受劣解概率太低用随机采样估算 delta_avgT0 取 10 倍左右解全是不可行解容量检查顺序错了检查是否先判断 load demand 再算取货后的负载多车路径表示混乱分隔符和客户编号混在一个数组里统一客户单独存序列车辆分段用拆分函数处理运行极慢但解质量差内循环迭代次数设置不合理内循环长度取客户数的5到10倍不要盲目加大目标函数中出现 NaN距离矩阵中有未赋值元素检查 dist 矩阵是否完整客户编号索引是否有错温度降到很低后结果不再变化算法已收敛到局部最优加入容量重定位算子或者周期性重置温度不同种子运行结果波动过大初始解随机性太强固定种子或者改用贪心构造初始解插入算子导致客户点缺失删除元素后索引错位使用逻辑索引或先记录节点再删除避免删除时跳过元素车辆数异常多初始解中分隔符过多初始化时用容量作为约束条件先安排路径再分配车辆载货量总是爆仓路径序列中节点顺序对容量影响大尝试使用 2-opt 前先模拟一遍负载曲线算出来的总距离不合理地大可能没有把返回仓库的距离计入检查 calc_total_cost 最后一段归程是否遗漏对换节点后立即不可行取货量和送货量差异大交换导致负载失衡交换前先计算两个节点的容量差值优先选择差值小的节点交换运行时间太长远超预期内循环和温度层循环都开得过大先跑小算例确定合理的运行时间再放大问题规模曲线出现阶跃上升温度仍偏高接受了劣解属于正常的全局探索阶段不用处理曲线下降后又上升再下降温度冷却过快局部搜索不充分调大 alpha降慢降温速度解的几何形态出现交叉2-opt 算子比例太低提高 2-opt 算子概率到 0.3 左右容量检查会误判的情况删除一个节点时未同时更新对应车辆负载拆分路径后逐段独立做容量检查大规模算例跑不出来矩阵索引和内存占用过大改用稀疏矩阵存储 dist或分批预计算输出解不满足“每个客户只服务一次”邻域操作时复制不完整在生成新解后增加唯一性校验断言代码在其他电脑上运行出错依赖的统计工具箱未安装检查是否使用 randi/rand 之外的工具箱函数这个表里的每一条基本都能对应到我实际调试过程中打过的日志、看过的曲线、翻过的文档。建议在跑自己的数据时把这些排查点当成 checklist 提前检查可以减少很多无效调试时间。7. 从代码到项目落地延伸方案的思考VRPSPD 的模拟退火代码本身是一个很好的算法骨架但真实业务往往比模型复杂得多。我做过几个相关的落地项目遇到最多的三个扩展需求是加时间窗约束、多中心协同、车型异构。如果你的项目需要把这些做进去可以从下面几个方向去扩展。时间窗约束是最常见的扩展。在 SA 框架里为客户增加时间窗检查函数位置在容量检查之后。时间窗的检查比容量更麻烦因为它不只是判断可行性还涉及等待时间和服务时间对后续节点的影响。这时候目标函数往往要从纯距离扩展到“距离 时间惩罚”的加权和。这个扩展在代码层面只增加一个函数但在算法层面会让解空间变得更加复杂收敛难度显著上升。我的建议是先跑通代码再慢慢加维度不要指望一步到位。多中心协同问题需要对编码做较大的改动。单中心的“断点式编码”不好用了得改成“中心-车辆-客户”三层关联的编码结构。这时模拟退火的邻域算子也需要调整除了客户点之间交换还跨中心调度车辆。每次中心改变后还要重算车辆返回哪个中心。这个工作量很大但如果你把项目定位从“学术 demo”转向“工程系统”这一步绕不过去。异构车队的改动相对温和只需要在容量检查时给每辆车配置不同的容量参数同时把路径分段和车辆之间的映射关系解耦。这类改动的“边际成本”不高但因为异质车辆的成本不同目标函数要做相应的权重调整。从我的实践经验来看把 SA 代码写好只是项目的前 50%后 50% 取决于你对问题的抽象能力和对业务约束的理解。算法本身并不神秘真正拉开差距的是“能不能把模糊的业务需求翻译成清晰的数学模型再落到健壮的工程实现上”。这也是为什么我在这篇分享里花大量篇幅讲建模的逻辑和调试的细节而不是只贴一段代码让你直接复制。如果你打算用这套代码做课程设计建议在完成基础实验后顺手做一个参数敏感性分析的小环节。导师最看重的是你对算法的理解深度和系统的实验设计能力这两点比代码本身更能体现水平。遇到具体问题拿不准时先画收敛曲线、记录中间解的变化过程再对照这个速查表逐项排查基本就能把问题定位到一个具体的函数或一行表达式上。