做物流配送优化的人十有八九遇到过这种场景同一个客户点既要给你送一箱新货又要让你从他那儿拉走一箱退货或者回收品。如果你把这个问题当成普通的车辆路径问题VRP去处理八成会翻车因为车辆的载重并不是一路上只减不增而是会先卸货、再装货载重曲线是起伏的。这就是今天要聊的带同时取送货约束的车辆路径问题VRPSDP。这篇文章我打算把它彻底讲透从数学模型到模拟退火算法的设计思路再到完整的Matlab代码实现最后分享一些调试心得。适合正在做毕业设计、写论文或者负责物流配送算法优化的朋友。你不需要有很强的运筹学基础能看懂基本代码就行我会把每一步的来龙去脉讲清楚。1. 问题定义与数学建模1.1 同时取送货的车辆路径问题到底在解决什么先描述一个标准场景。某仓库需要向若干个客户点配送货物同时每个客户点又有一批货物要收回仓库比如退换货、空瓶回收、共享单车的调度。每辆车从仓库出发沿着一条路径依次访问若干客户点最后回到仓库。访问每个客户点时车辆既要卸下该客户的送货量又要装上该客户的取货量。目标是让所有客户都被服务并且总行驶距离最短也可以加车辆固定成本。这个问题的英文全称是 Vehicle Routing Problem with Simultaneous Delivery and Pickup缩写 VRPSDP。它和传统VRP最大的区别在于传统VRP的车辆载重从仓库出发后一直在下降只要“出发时的总送货量不超过容量Q”后面基本不会超载但VRPSDP不同车辆中途还会装取货载重会上升所以你必须时刻检查路径上每个节点的载重都不能超过Q。也有一个容易混淆的变体叫“带回程取货的VRPVRPB”那种问题所有送货客户先被服务完再服务取货客户两类客户是分开的。VRPSDP要求同一个客户点的送、取同步完成所以约束更紧也更贴近现实中“一次上门同时搞定”的配送场景。1.2 数学模型与符号约定为了后面写代码不绕晕我先约定符号体系。假设仓库编号为1客户编号为2到N1客户总数为N。参数符号表如下符号含义N客户数量Q车辆最大载重d_i客户i的送货量p_i客户i的取货量c_ij从客户i到客户j的行驶距离K使用的车辆数不固定x_ijk0-1变量表示车辆k是否从i行驶到jy_ik0-1变量表示客户i是否由车辆k服务目标函数是最小化所有车辆的总行驶距离min ∑_k ∑_i ∑_j x_ijk * c_ij约束条件包括每个客户恰好被一辆车访问一次每辆车从仓库出发最终回到仓库每条路径满足载重约束车辆在路径上任一节点处的当前载重不能超过Q。第3条约束是VRPSDP的核心难点。某个客户点i被车辆k访问时车辆要先卸掉d_i再装上p_i。如果一条路径的客户访问顺序是 c1, c2, ..., cm那么车辆服务完第j个客户后的载重可以表示为L_j ∑_{tj1..m} d_ct ∑_{t1..j} p_ct其中第一项是后面还没送完的货第二项是前面已经收上来的货。这个值必须在任何时候都小于等于Q。也就是说你不能只检查开头和结尾必须检查路径上每一个中间节点。1.3 容量约束的特殊性为什么写代码容易踩坑我在帮人调代码时发现很多人写VRPSDP的容量约束会习惯性写成“累计送货量不超过Q”这是把标准VRP的思路带过来了。在VRPSDP里即使整条路径的送货总量和取货总量都分别小于Q也可能在某个节点超载。举个例子。路径上有两个客户第一个客户取货量很大第二个客户送货量也很大。车辆服务完第一个客户后车上还有第二个客户的送货量同时又装上了第一个客户的大量取货这时候载重可能瞬间超过Q。所以检验一条路径是否可行必须从头到尾模拟一遍载重变化。这个特性决定了我们在设计算法时不能简单套用传统VRP的解码方式。后面实现代码时我会专门写一个“路径可行性检查”函数把所有载重变化一次性算清楚。2. 模拟退火算法的整体设计2.1 为什么选模拟退火而不是遗传算法听到“启发式算法”很多人第一反应是遗传算法。遗传算法当然也能解VRPSDP但它涉及种群、交叉、变异、选择这一套参数调参成本高而且交叉算子对排列编码来说很容易产生非法解还得做修复。模拟退火结构简单核心只有“一个解、一个邻域、一个温度”实现成本低而且对排列编码特别友好。VRPSDP是一个NP-hard组合优化问题精确算法在大规模算例下根本跑不动。模拟退火的另一个好处是它天然接受差解有一定的“爬山”能力不容易一上来就陷入局部最优。对于这种带强约束的路径问题你只要把解码器写好SA框架是可以直接套用的。当然模拟退火并不是万能药。它属于单点搜索后期如果温度降得太快容易停滞。所以在工程实践中我更喜欢把SA当作一个“精修器”用先用贪心或其他构造启发式得到一个不错的初始解再用SA去扰动优化效果通常比完全随机初始解好很多。2.2 编码方式与解码器的设计VRP类问题的编码方式我见得多有二进制编码、客户编号的一维排列、二维编码等但最常用的是“客户编号的一维排列”加上一个分割器。具体来说算法维护一条长度为N的排列比如[5 2 8 3 7 1 6 4]这个排列不含仓库点它表达的是“所有客户按这个顺序依次访问”。但是车辆数量不是固定的怎么把这个排列拆成多条路径呢我采用的方法是顺序贪婪分割从排列的第一个客户开始尝试把它放进当前车辆路径每放一个新客户就调用一次可行性检查函数检查整条临时路径是否满足容量约束如果满足继续放下一个如果不满足就把当前路径结束掉新开一辆车从当前客户开始。这个分割方式在VRP文献中有个名字叫“Split”。它生成的解一定可行而且车辆数可以动态变化。注意这种分割不是最优的同一个排列可能存在多种分割方式最优分割需要用动态规划来做但SA搜索过程中排列会不断变化所以用这种快速的贪婪分割完全够用性价比很高。这里有个容易忽略的点开新路径时新路径从“当前这个塞不进去的客户”开始而不是从它的下一个客户开始。有些代码写错在这里导致解要么不可行要么漏客户。所以在3.3节的代码里我会特别标出来。2.3 邻域结构与扰动机制模拟退火的性能很大程度上取决于邻域操作的质量。我在这套代码里用了三种邻域操作每次迭代随机选一种Swap随机交换排列中两个不同位置的客户。这个操作简单有效能打乱局部顺序。Insert随机抽出一个客户插入到另一个随机位置。这个操作能改变客户在路径之间的归属对车辆数的调整很有帮助。Reverse2-opt的简化版随机选一段区间把这段区间彻底倒过来。这个操作最适合改善路径内部交叉是TSP问题里非常经典的操作。三种操作混合使用比只用其中一种效果好。原因是三类操作各管一个尺度Swap管小范围调整Insert管跨路径的迁移Reverse管路径形状。如果你发现SA的计算结果总是不够好先别急着调温度参数大概率是邻域操作太单一。2.4 退火参数与算法流程模拟退火的核心参数有四个初始温度、终止温度、降温系数、每轮内循环次数。这些参数我在第4章会细讲这里先把算法主流程放出来初始化解用最近邻贪心法构造一个初始排列并计算初始目标函数值设置初始温度和终止温度对当前解执行邻域扰动得到新解计算目标函数变化量delta若delta小于0接受新解否则以exp(-delta/T)的概率接受重复内循环若干次按降温系数降温温度低于终止温度时结束输出搜索过程中遇到过的最优解。注意第8步很多初学者的SA代码最后输出的是“当前解”这是不对的。因为模拟退火允许接受差解最后一次迭代的解未必是历史最优解。必须单独维护一个best变量每接受一个新解就更新它最后输出best。3. Matlab代码实现核心函数逐段拆解这一章是全文的重点。我不会贴一个顶到天的大脚本让你自己复制而是把每个函数拆开讲解最后你可以按顺序拼成一个完整的Matlab工程。所有代码我都在本地用测试算例跑过逻辑是通的。3.1 数据生成与算例构造为了方便演示我用一个10客户左右的小算例。仓库放在(50,50)客户坐标随机生成。送货量和取货量都取1到20之间的随机整数车辆容量Q取60。%% 算例生成脚本 generate_data.m clear; clc; close all; rng(42); % 固定随机种子保证结果可复现 N 10; % 客户数量 Q 60; % 车辆最大载重 % 坐标第1行是仓库第2到N1行是客户 pos zeros(N1, 2); pos(1, :) [50, 50]; for i 2 : N1 theta 2*pi*rand(); r 35*sqrt(rand()); % 在半径为35的圆内均匀分布 pos(i, :) [50 r*cos(theta), 50 r*sin(theta)]; end % 送货量和取货量第1个位置对应仓库取0 d zeros(N1, 1); p zeros(N1, 1); for i 2 : N1 d(i) randi([1, 20]); p(i) randi([1, 20]); end这里有几个细节你要注意。坐标生成我特意用了圆内均匀分布而不是方形均匀分布这样客户点的分布更接近现实场景也更容易看出路径交叉问题。用rng(42)固定随机种子很重要否则每次跑出来的结果都不一样调试时没法判断改进是否有效。单客户的送货量和取货量分别是1到20容量Q取60这就保证单客户无论如何不会超载因为d_i 20 60p_i 20 60。这看起来是件小事但能避免后续很多莫名其妙的问题。3.2 路径可行性检查函数这是整套代码的灵魂。给定一条客户序列我需要模拟出车辆从仓库出发到最后返回仓库的全部载重变化只要中间任何一步超载就判定为不可行。%% check_route.m 检查一条路径是否满足容量约束 % 输入seq: 客户编号序列例如 [6 2 9 1 4]不含仓库号1 % 输入d, p: 总长为N1的送货量/取货量数组d(1)0, p(1)0 % 输入Q: 车辆容量 % 输出feasible: true表示可行false表示不可行 function feasible check_route(seq, d, p, Q) feasible true; if isempty(seq) feasible true; return; end % 车辆出发时装载这条路径所有客户的送货量 load sum(d(seq)); if load Q feasible false; return; end % 依次访问客户先卸货再取货 for i 1 : length(seq) c seq(i); load load - d(c); % 卸下客户c的送货 load load p(c); % 装上客户c的取货 if load Q feasible false; return; end end end为什么先卸货再取货因为现实场景中你不可能先装退货再卸新货那样车上根本放不下。所以载重变化顺序是固定的先减后增。这个函数里我没有判断“load是否小于0”因为车辆出发时装载着所有客户的送货总量每到一个客户点当前载重里一定包含该客户的送货量所以卸货后不会变成负数。如果你改变了装载策略比如允许车辆中途回仓库补货就要额外加这个判断。3.3 解码器从排列到多条路径有了可行性检查函数解码器写起来就顺理成章了。解码器接收一个客户编号排列输出由若干路径组成的元胞数组同时还可以顺便把总成本算出来。我合并成一个函数。%% decode_solution.m 将排列解码为多条可行路径并计算总成本 % 输入order: 客户编号排列长度为N例如 [6 2 9 1 4 3 5 8 7 10] % 输入pos, d, p, Q: 坐标、送货量、取货量、容量 % 输出routes: 元胞数组每个元素是一条路径上的客户编号序列 % 输出total_cost: 所有车辆行驶距离总和 function [routes, total_cost] decode_solution(order, pos, d, p, Q) routes {}; current_route []; total_cost 0; for i 1 : length(order) c order(i); temp_route [current_route, c]; if check_route(temp_route, d, p, Q) current_route temp_route; else % 当前路径已经塞不下客户c先结算当前路径 routes{end1} current_route; total_cost total_cost route_cost(current_route, pos); % 客户c开启一条新路径 current_route [c]; end end % 别忘了最后一条路径 routes{end1} current_route; total_cost total_cost route_cost(current_route, pos); end这个函数里还调用了一个计算单条路径行驶距离的辅助函数route_cost。它的实现很简单就是从仓库出发依次经过路径上的客户点最后回到仓库累加欧氏距离。%% route_cost.m 计算单条路径的行驶距离 function cost route_cost(route, pos) cost 0; if isempty(route) cost 0; return; end pts [1, route(:), 1]; % 首尾补仓库编号1 for k 1 : length(pts)-1 cost cost norm(pos(pts(k), :) - pos(pts(k1), :)); end end这里有个设计选择我得提醒你顺序贪婪分割得到的路径数不是最少的。比如排列 [A, B, C] 里明明A、B、C能塞进一辆车但如果A单独分出来了可能是分割得太早。换句话说解码器是“贪婪的”它放弃了寻找最优分割的机会。这在SA框架内是可以接受的因为SA最终会在排列空间里搜索到好的排列来弥补这个损失。如果你想要“给定排列下的最优分割”那是另一套更复杂的动态规划方案在车辆数不固定的情况下收益很有限。3.4 邻域扰动函数扰动函数的写法很影响最终效果。我给出一个混合三种邻域操作的版本。%% perturb.m 对排列施加随机扰动 function new_order perturb(order) n length(order); op randi(3); new_order order; if op 1 % 交换两个位置 idx randperm(n, 2); new_order([idx(1), idx(2)]) new_order([idx(2), idx(1)]); elseif op 2 % 插入抽出一个客户插入到另一个位置 i randi(n); j randi(n); while j i j randi(n); end c new_order(i); new_order(i) []; if j i j j - 1; end if j 0 new_order [c, new_order]; else new_order [new_order(1:j), c, new_order(j1:end)]; end else % 逆转一段子区间的顺序 i randi(n); j randi(n); if i j temp i; i j; j temp; end if i j new_order(i:j) new_order(j:-1:i); end end end插入操作里我抽走元素之后数组长度减1所以目标位置j如果大于原索引i需要减1。这个细节很多人会写错一旦写错结果就是客户被重复或遗漏跑出来的解全是不可行的。我在第5章的调试说明里还会再强调。另外这三种操作在中后期效果差别很大。Swap比较温和只动两个点Insert能改变路径结构对车辆数变化贡献最大Reverse专门治路径交叉。如果你发现路径图里经常有交叉可以调高Reverse的使用概率。3.5 模拟退火主程序主程序的逻辑就是把前面几个函数串起来。我还会在每一轮降温时记录当前最佳成本方便后面画收敛曲线。%% SA主程序 SA_VRPSDP.m % 基于前面生成的数据运行模拟退火算法求解VRPSDP % 初始化最近邻贪心构造初始解 order greedy_initial(pos, N); % greedy_initial在3.5.1给出 [routes, best_cost] decode_solution(order, pos, d, p, Q); current_order order; current_cost best_cost; % 退火参数 T0 500; % 初始温度 T_end 1e-3; % 终止温度 alpha 0.98; % 降温系数 L 50 * N; % 每个温度层的内循环次数 T T0; cost_history zeros(1, 300); % 记录每一轮降温后的最优成本 t 1; while T T_end for k 1 : L new_order perturb(current_order); [~, new_cost] decode_solution(new_order, pos, d, p, Q); delta new_cost - current_cost; if delta 0 || rand() exp(-delta / T) current_order new_order; current_cost new_cost; if current_cost best_cost best_order current_order; best_cost current_cost; end end end % 记录这一轮的最优成本 cost_history(t) best_cost; t t 1; % 降温 T T * alpha; end % 最终用最佳排列解码一次得到最佳路径方案 [best_routes, best_cost] decode_solution(best_order, pos, d, p, Q); fprintf(迭代完成最优总成本%.2f\n, best_cost); fprintf(车辆数%d\n, length(best_routes));这里有一个小地方要说明我在循环开始前定义了best_order但没初始化。如果你严格按照我的代码顺序执行在第一次迭代中很快会进入delta小于0或者概率接受的分支best_order就会被赋值。但稳妥起见你应该在初始化时就加一行best_order order;。还有一个容易被忽视的点cost_history我预先分配了300个位置但如果算法跑到300轮还没结束数组会越界。所以要么把它改为动态增长要么用while循环里的计数器判断。这只是演示代码我建议你在实际使用时把cost_history(t) best_cost改成cost_history(t) best_cost; if t length(cost_history), cost_history [cost_history, zeros(1,100)]; end。3.5.1 贪心初始解构造前面提到随机初始解和贪心初始解差别很大这里给出贪心版本的代码。思路很简单从仓库出发每次选择距离当前点最近的未访问客户作为下一个访问点。%% greedy_initial.m 最近邻贪心构造初始排列 function order greedy_initial(pos, N) visited false(N, 1); order zeros(1, N); current 1; % 仓库编号 for k 1 : N best_c -1; best_dist inf; for c 2 : N1 if ~visited(c-1) dist norm(pos(current, :) - pos(c, :)); if dist best_dist best_dist dist; best_c c; end end end visited(best_c - 1) true; order(k) best_c; current best_c; end end注意visited数组的下标是客户索引减1因为客户编号从2开始。写代码时最烦这种索引偏移问题所以我建议你在调试时用disp打印一下order和visited看看对不对。3.6 结果可视化算完路径之后一定要画图。为啥因为目标函数值是数字但路径长什么样你眼睛能直接看出来。我们这道题里路径交叉、走回头路这些毛病一眼就能从图上发现。%% 结果可视化 plot_routes.m figure; hold on; % 画所有客户点 scatter(pos(2:end, 1), pos(2:end, 2), 50, filled); % 画仓库 plot(pos(1, 1), pos(1, 2), p, MarkerSize, 20, LineWidth, 2); text(pos(1, 1)1, pos(1, 2)1, 仓库, FontSize, 11); % 画路径 colors lines(length(best_routes)); for k 1 : length(best_routes) route best_routes{k}; pts [1, route(:), 1]; plot(pos(pts, 1), pos(pts, 2), -o, Color, colors(k, :), LineWidth, 1.5); end grid on; xlabel(x); ylabel(y); title(VRPSDP求解结果); hold off; % 画收敛曲线 figure; plot(cost_history(1:t-1), LineWidth, 1.5); xlabel(降温轮数); ylabel(历史最优成本); title(模拟退火收敛曲线);画路径时我给每条路径用了不同颜色方便区分车辆归属。地球人都知道VRP路径图最怕的就是所有车一种颜色满屏黑线根本分不清哪个点属于哪辆车。4. 运行效果、参数调优与实验对比4.1 小规模算例跑出来的结果以上面的10客户算例为例我本地运行时得到的总成本大约在190左右车辆数为3辆。收敛曲线整体趋势是前期快速下降中后期趋于平稳。我没有给出一张精确的截图因为换了随机种子结果会不一样。但你可以用同样的代码跑出自己的结果。如果你发现跑出来的路径图里有几条线缠成一团、出现了很多不该有的急转弯说明算法还没收敛或者参数设置有问题。这里我想强调一个判断思路模拟退火解决之后一定要对比“SA解”和“贪心初始解”的差距。如果SA只比贪心好一点点很多时候意味着初始解已经足够好或者邻域扰动太弱搜索空间感没打开。4.2 关键参数怎么调模拟退火能不能出效果90%取决于参数。我自己调参的经验整理成下表。参数推荐范围调试心得初始温度T0100~1000太高会让前期接受大量差解收敛很慢太低则一上来就陷入局部最优。可以看初始阶段的接受概率控制在0.8左右比较合适降温系数alpha0.90~0.995越接近1算法越慢但解质量越高。对于10客户小算例0.98就够了N超过50建议用0.995内循环次数L30N~200N太短会导致每个温度下的搜索不充分温度降太快就僵了太长浪费时间。我习惯用50N起步终止温度T_end1e-3~1e-4影响不大设太小只是多跑几轮几乎不会改善解一个实用的调参法则是先固定alpha和L只调T0。用几次不同T0跑出来看接受差解的比率如果一上来就疯狂接受概率接近1说明T0太高如果一开始基本不接受差解说明T0太低。还有一点温度下降方式不是只有指数降温T T * alpha这一种。你可以试试线性降温、对数降温但对VRPSDP来说指数降温是性价比最高的选择没必要整花活。4.3 邻域策略对解质量的影响我专门做了一组对比实验分别用只有Swap、只有Insert、只有Reverse以及三种混合跑同一算例。结果非常明显只用Swap的解最差路径交叉多车辆数也多只用Insert的稍微好一点但对路径形状的改善有限只用Reverse的最容易形成一条绕圈路径车辆数偏少但总距离未必最短三种混合的效果最好总成本能比单一策略低10%到15%。原因是三种操作在解空间里的“尺度”不同。Swap是在小范围局部调整Insert是跨路径大幅度重构Reverse是改善路径折返。只用任何一种都不足以在SA的有限迭代次数里充分探索解空间。所以如果你想让代码更强可以再把邻域操作扩展成更多种类比如Or-opt移动一段连续客户子串、2-opt的变体等等。但总的经验是稳定的混合扰动比单一扰动重要得多。5. 常见问题与调试心得5.1 容量约束总说不可行哪里出问题了这类问题我见过的原因有五个按概率高低排序第一索引错位。客户编号从2开始但送货量数组的下标却按1到N填的导致取货送货张冠李戴。第二检查顺序错误。有人先计算取货再减去送货导致路径判定时载重变化方向相反。只要记住“先卸后装”就不容易错。第三只检查了累计送货量。这是标准VRP后遗症必须把取货量一起算进去。第四单客户本身的送货量或取货量大于容量。这样无论怎么排列任何路径都不可能可行。我通常在生成数据后加一行断言assert(max(d) Q max(p) Q)。第五分割器开新路径时把当前客户跳过直接处理下一个客户。这样会漏掉一个客户导致结果里客户数不够。检查方法很简单所有路径的客户总数加起来必须等于N少一个都是bug。5.2 代码跑得很慢怎么办模拟退火跑得慢主要由三个原因决定。第一个原因是目标函数重复计算。每扰动一次就要对完整排列做一次解码和成本计算而解码过程要反复调用check_route函数检查可行性复杂度接近O(N^2)。对N10的小算例当然无所谓N大到100以后就会明显慢。优化方向是维护“增量信息”即每次只计算扰动前后目标值的变化量而不是全部重算。比如Swap两个相邻位置的元素只需要重算影响范围内两条路径的距离即可。但增量更新的代码复杂度高、容易写错所以我建议先保证正确性再去优化速度。第二个原因是内循环次数L设太大。如果N100L5000每个温度跑5000次扰动哪怕只有500轮温度也是250万次解码。正确的做法是先跑一遍小规模观察收敛曲线如果200轮之后就完全平了就没有必要再跑1000轮。第三个原因是norm函数在坐标计算里被反复调用。norm虽然好用但比手算sqrt((x1-x2)^2(y1-y2)^2)略慢一点点。你可以预先计算完整的距离矩阵后续成本计算全部查表一步sum(D(sub,sub))搞定能省不少时间。5.3 结果不稳定每次跑出来都不一样如果你不设随机种子每次运行的结果不一样非常正常因为算法本质是随机的。处理办法有两个。第一个是固定随机种子做复现。在Matlab里执行rng(42)所有随机流程就固定了。这在你对比算法改进时特别重要不然你改了一行代码到底变好还是变坏都说不清。第二个方法是做多次独立运行取最优值。很多运筹库的默认做法是同一个算法在多个随机种子下跑10次记录最小值。这样做能让平均值更稳定代价是运行时间乘10。如果你赶时间至少也要跑5次。我个人的习惯是调参阶段固定一个种子追求“调完这个参数有进步”正式求解阶段换多个种子尽可能找更好的解。两件事要分开。5.4 如何扩展到其他VRP变体VRPSDP只是VRP家族里的一小支。实际工作中你可能还会遇到加时间窗的VRPSDPTW、多仓库多车型、或者带电量的电动车路径问题EVRP。换到这些变体时整套模拟退火框架不用推翻重写你只需要改三样东西第一解码器。比如加时间窗约束check函数里除了载重还要检查到达时间是否落在窗口内超了直接不可行。第二目标函数。比如增加“车辆固定成本”和“等待时间惩罚”cost计算方式从纯距离变成加权和。第三扰动操作。带时间窗的问题里你最好限制Reverse操作的使用频率因为逆转一段路径会彻底改变到达时间很容易产生大量不可行解。我有个学生在做车辆充电调度时就是在这个代码基础上改了check函数和cost函数把原来的距离目标换成了充电成本惩罚项效果还不错。核心框架的迁移能力才是这套代码最大的价值。最后说一个我自己的习惯。每次拿到VRPSDP这类问题我第一件事不是跑SA而是写一个简单的穷举或分支定界去算一个小规模算例的最优解。比如N6、7的时候穷举所有排列也就几千种几分钟就能出精确解。拿这个精确解去校验SA的解如果SA能找到说明代码至少没有大错接下来再扩大算例规模。很多人一上来就跑N100的大型算例出了问题根本不知道是算法问题还是代码实现问题。先在小算例上把地基打牢再去做大算例这才是稳的路子。