做运筹和物流优化的朋友对模拟退火算法应该都不陌生但一碰到“同时取送货车辆路径问题”想手写一套能跑的Matlab代码很多人会卡在容量约束那块。这篇就专门解决这个痛点用模拟退火算法求解同时取送货车辆路径问题并给出完整可运行的Matlab实现。我会先把问题模型讲清楚再拆解算法设计思路最后贴代码、跑案例、聊参数尽量让拿到代码的人能直接改数据往自己业务上套。背景上同时取送货车辆路径问题在现实中遍地都是社区团购配送既要送菜又要回收前一天的空筐饮料批发商送新桶的同时要拉回空桶快递驿站既要派件又要揽收退件。这些场景都可以归纳成同一个数学结构车辆从配送中心出发为若干客户服务每个客户有一个送货需求和一个取货需求车辆在途中既要卸货又要装货容量需要全程不超限目标是找出一组总成本最小的路线。对于刚入门的研究者或者工程师这个问题比纯送货车更烧脑原因在于容量约束是“动态”的车上的货物量会随着服务进度变化不是简单的一车送完就结束。如果直接套用普通VRP的代码忽略装载量变化算出来的路线十有八九在现场执行不了。因此本文从模型、算法到Matlab代码全部围绕“同时取送货容量实时更新”来写适合需要快速实现原型、验证算法效果的读者。1. 问题模型同时取送货车辆路径问题到底在解什么1.1 从快递员“既送又收”讲起想象一个快递员骑着三轮车车上装满了今天要派送的快件同时还要沿途收取退换货的包裹。三段式客服话术永远都是“您好有您的快递顺便问一下有没有要退的件”。如果这个快递员服务20个点位他的麻烦在于一开始车上的空间被派送件占满每送掉一件就腾出一点空间但每收一件又会占回一点空间。怎么安排先后顺序才能让全程既不爆仓、又不绕路这就是同时取送货车辆路径问题的直观画面。术语上通常叫VRPSDPVehicle Routing Problem with Simultaneous Delivery and Pickup与“先送后取”不同它允许同一个客户点同时产生送货和取货车辆在同一个点既卸又装。对于每个客户i给定坐标、送货量、取货量。配送中心也有坐标车辆容量为Q。问题的目标是最小化总行驶距离也可以扩展为时间或费用同时满足每个客户必须被访问且仅访问一次车辆不得超容车辆从中心出发并返回中心。1.2 数学模型与容量约束写成标准的数学形式一般用三下标变量x_{ijv}表示车辆v是否从i行驶到j。目标函数为所有车辆行驶距离之和最小。常见约束包括流量守恒约束、每个客户被访问一次约束以及最关键的行车中容量约束。容量约束不能写成简单的“累计装载量不超过容量”因为装载量在每一个客户点之后都会改变。用生活化的语言理解容量约束设一辆车服务客户序列S[s1, s2, ..., sk]出发时车上装载的货物必须等于该车所有客户的送货需求之和。之后每访问一个客户执行“先卸后装”先减去该客户的送货量再加上该客户的取货量。于是车辆在这条路径上任意时刻的装载量都是一条动态变化的曲线。标准做法是让装载量变化满足两个上限出发装载量不超过容量且路径中任意位置装载量不超过容量。实践中可以统一用最大装载量来检查。这也就是后续代码里评估函数的核心逻辑。1.3 为什么这条路不好走复杂度分析很多人第一次在Matlab里实现VRPSDP会觉得“不就多了一个取货量嘛把TSP改改就行”。实际没那么简单。纯TSP是若干点的排列组合问题而VRPSDP每多一个取货需求就多一层容量可行性判断。它的解空间复杂度从组合爆炸升到组合爆炸加约束灾难属于强NP-hard问题。精确算法例如分支定界在客户数超过30时往往要跑上小时甚至无法收敛。因此工程上更多依赖启发式算法和元启发式算法。在元启发式中遗传算法需要设计交叉变异算子参数一多就难以调稳禁忌搜索需要维护禁忌表实现细节多粒子群算法对离散排序问题的编码绕来绕去。相比之下模拟退火算法结构简单、参数解释性强、跳出局部最优的能力有理论保证非常适合作为VRPSDP的快速求解原型。这也是为什么我在这篇文章里选择模拟退火而不是其他更花哨的算法。2. 模拟退火算法为什么适合这个问题2.1 启发式方法对比我用一张表总结常见方法的取舍大家可以根据自己的场景对号入座。方法实现复杂度跳出局部最优能力参数敏感度适用规模精确算法分支定界高全局最优高小规模客户数30遗传算法中高较强依赖多样性高中大规模禁忌搜索中高强靠禁忌表中中大规模模拟退火低强靠温度退火低几十到几百客户从工程落地角度模拟退火最大的优势是“一套循环走天下”只要定好编码、目标函数、邻域操作和冷却计划前面模型改成纯送货、带时间窗、多车场内核都不需要大改。经验不足的人也可以快速跑通一条基线结果后面再跟其他算法对比。2.2 核心思想从高温乱试到低温精细模拟退火的名字来自金属热处理。高温下原子运动剧烈系统可以跳到高能状态随着温度降低原子逐渐趋于稳定排列。优化问题中“能量”就是目标函数值“状态”就是一条路径方案。算法每步随机生成一个邻域解如果目标值下降就接受如果上升则以概率pexp(-ΔE/T)接受。ΔE是目标增量T是当前温度。温度高时p接近1意味着算法敢接受差解从而跳出局部最优温度低时p趋于0算法逐渐收敛。关键要理解“接受差解不是为了胡闹而是为了换一条路翻过山峰”。比如一个解在局部谷底如果只接受更优解算法永远出不来而一个较差解也许能引导到另一个更优谷底。这个机制让模拟退火在VRPSDP这类高度非线性的搜索空间里比贪心策略稳健得多。2.3 编码、目标函数与邻域结构的设计思路VRPSDP的编码方式直接决定后续算法好不好写。我采用最简单的客户序列编码用一个1×N的行向量保存所有客户编号的排列例如[3 5 1 2 4]。这个序列不代表单车路线而是所有待服务客户的“大序列”。具体几辆车、每辆车访问哪些客户由容量约束动态分割产生。用生活化比喻这一系列客户像一条待切分的香肠评估函数在扫描时如果当前车辆塞不下后续客户就从那里切一刀新开一辆车。这种方式的好处是模拟退火的邻域操作只需要重排序列不需要额外处理车辆数车队规模自动随解变化。目标函数定义为obj totalDist vehiclePenalty × numVehicle。为什么要加车辆数惩罚因为单纯最小化距离算法会把很多客户塞到同一辆车里只要能塞下哪怕绕远路也愿意这会造出不现实的路线。加入车辆数项后可以平衡行驶里程和用车数量。车辆数惩罚的设置也有讲究这我在参数章节会细说。邻域操作我一般维护三种混合使用交换swap随机选两个位置互换客户插入insert随机选一个客户插入另一个随机位置2-opt反转随机选一段子序列反转其顺序。三种操作分别对应改变车辆间分配、改变路径内顺序、消除交叉绕路。混合使用比只用一种更不容易陷入固定套路。3. Matlab代码实现3.1 代码文件结构和运行流程我用的是Matlab R2020a纯脚本加函数实现不需要额外优化工具箱。整个项目建议分成四个文件main_sa_vrpspd.m主脚本负责生成数据、初始化参数、执行模拟退火、输出结果并绘图evaluateSolution.m求解评估函数输入一个客户序列输出总距离、使用车辆数和每个车辆的路径neighborSolution.m邻域生成函数随机执行交换、插入或2-opt反转plotSolution.m可选的绘制函数把最终路径画在坐标平面上。三个函数文件保存时文件名必须和函数名严格一致否则Matlab会报“未定义函数”。这也是新手踩得最多的坑。3.2 主脚本完整的模拟退火流程下面这段代码可以直接复制到一个新脚本里运行我故意没有用太多技巧尽量把每一步写清楚方便调试。% main_sa_vrpspd.m % 同时取送货车辆路径问题 - 模拟退火算法求解示例 clc; clear; close all; rng(42); % 固定随机种子便于复现 %% 1. 生成测试数据 N 20; % 客户数量 data.x rand(1, N) * 100; % 客户横坐标 data.y rand(1, N) * 100; % 客户纵坐标 data.deliver randi([1, 10], 1, N); % 送货需求 data.pickup randi([1, 8], 1, N); % 取货需求 data.capacity 30; % 车辆容量 data.depotX 50; data.depotY 50; %% 2. 模拟退火参数 T0 1000; % 初始温度 Tend 1e-3; % 终止温度 alpha 0.99; % 降温系数 L 200; % 每个温度下的迭代次数 vehiclePenalty 100; % 车辆数在目标函数中的惩罚系数 %% 3. 初始化 curOrder randperm(N); % 随机初始客户序列 [curCost, curVehicleNum, ~] evaluateSolution(curOrder, data); bestOrder curOrder; bestCost curCost; T T0; costHistory zeros(1, 0); % 记录每一轮外层结束后的目标值 %% 4. 模拟退火主循环 while T Tend for k 1:L newOrder neighborSolution(curOrder); newCost evaluateSolution(newOrder, data); delta newCost - curCost; if delta 0 || rand exp(-delta / T) curOrder newOrder; curCost newCost; if curCost bestCost bestOrder curOrder; bestCost curCost; end end end T alpha * T; costHistory(end 1) bestCost; end %% 5. 输出与绘图 [bestCost, bestVehicleNum, bestPaths] evaluateSolution(bestOrder, data); fprintf(最优总距离%.2f\n, bestCost - vehiclePenalty * bestVehicleNum); fprintf(使用车辆数%d\n, bestVehicleNum); for v 1:length(bestPaths) seq bestPaths{v}; fprintf(车辆%d路线0 - , v); fprintf(%d - , seq); fprintf(0\n); end figure; plotSolution(bestOrder, data, bestPaths); title([模拟退火结果目标, num2str(bestCost)]); figure; plot(costHistory); xlabel(外层迭代次数); ylabel(最优目标值); title(收敛曲线);3.3 关键子函数路径分割、容量检查与目标函数计算evaluateSolution是整个代码的核心新手往往在这里翻车。函数需要完成三件事按顺序尝试把客户塞给当前车辆同时维护动态装载量塞不下就开下一辆车完成分割后计算所有车辆的距离和返回目标值、车辆数和路径。function [totalCost, vehicleNum, paths] evaluateSolution(order, data) % 输入order客户访问序列1xN整数向量 % 输出totalCost目标函数值包含距离和车辆数惩罚 % 输出vehicleNum所需车辆数 % 输出paths元胞数组每个元胞是该车辆访问的客户序列 x data.x; y data.y; del data.deliver; pk data.pickup; cap data.capacity; depotX data.depotX; depotY data.depotY; paths {}; idx 1; N length(order); while idx N currentPath []; % 尝试把后续客户逐个加入当前路径 while idx N candidate [currentPath, order(idx)]; % 容量检查 maxLoad checkMaxLoad(candidate, del, pk); if maxLoad cap currentPath candidate; idx idx 1; else break; end end paths{end 1} currentPath; end vehicleNum length(paths); % 计算总距离 dist 0; for v 1:vehicleNum seq paths{v}; if isempty(seq) continue; end prevX depotX; prevY depotY; for i 1:length(seq) c seq(i); dist dist sqrt((x(c) - prevX)^2 (y(c) - prevY)^2); prevX x(c); prevY y(c); end dist dist sqrt((depotX - prevX)^2 (depotY - prevY)^2); end vehiclePenalty 100; % 与主脚本保持一致 totalCost dist vehiclePenalty * vehicleNum; end function maxLoad checkMaxLoad(seq, del, pk) % 检查给定路径在同时取送货条件下的最大装载量 load sum(del(seq)); % 出发装载量等于该车所有客户的送货需求之和 maxLoad load; for i 1:length(seq) c seq(i); load load - del(c) pk(c); if load maxLoad maxLoad load; end end end这里有个经验点最大装载量可能不是出现在路径开头而是出现在中途。原因很简单如果先访问送货量大的客户后访问取货量大的客户装载量就会先降后升峰值不一定在起点。所以容量检查函数必须逐点扫描。3.4 邻域操作函数交换、插入与2-opt反转这个函数负责生成新解。我采用一个随机数决定三种操作中的一种。注意邻域操作不检查可行性可行性统一交给评估函数判断这大大简化了编程。function newOrder neighborSolution(order) % 随机选择一种邻域操作 % 1: 交换两个位置 % 2: 插入一个客户到另一个位置 % 3: 2-opt反转一段子序列 newOrder order; n length(order); r randi([1, 3]); switch r case 1 i randi([1, n]); j randi([1, n]); while j i j randi([1, n]); end newOrder([i, j]) newOrder([j, i]); case 2 i randi([1, n]); j randi([1, n]); while j i j randi([1, n]); end c newOrder(i); newOrder(i) []; if j i j j - 1; % 删除后的位置调整 end newOrder [newOrder(1:j), c, newOrder(j1:end)]; case 3 i randi([1, n]); j randi([1, n]); while j i j randi([1, n]); end newOrder(i:j) fliplr(newOrder(i:j)); end end插入操作的索引调整是最容易写错的地方。如果先删除客户i再插入位置j时j需要根据原序列位置做减一处理。上面代码用if j i, j j - 1处理能避免边界错误。3.5 绘图函数绘图对验证算法效果很有用特别是检查车辆之间的路线是否交叉、是否出现不合理绕路。我提供一个简单版本只画路径和坐标点。function plotSolution(~, data, paths) % 绘制最终路线 figure; hold on; plot(data.depotX, data.depotY, kp, MarkerSize, 12, MarkerFaceColor, k); colors lines(length(paths)); for v 1:length(paths) seq paths{v}; xp [data.depotX, data.x(seq)]; yp [data.depotY, data.y(seq)]; plot(xp, yp, o-, Color, colors(v, :), LineWidth, 1.2); end plot(data.x, data.y, ro, MarkerFaceColor, r); xlabel(X坐标); ylabel(Y坐标); grid on; axis equal; end4. 参数调优与实验结果分析4.1 关键参数的经验取值模拟退火的性能强依赖冷却计划和迭代次数。我把常用参数的经验值整理成一张表但要注意这些值需要根据问题规模调整。参数作用常用范围我的建议T0初始温度使初始接受概率约0.8-0.95可取初始目标值的5-20倍Tend终止温度接受概率低于0.0011e-3到1e-5alpha降温系数0.9-0.999客户数少用0.95客户多用0.99L每个温度迭代次数N×10到N×5020个客户用200100个客户用1000-2000vehiclePenalty车辆数惩罚与平均单车行驶距离同数量级暂定100再观察车辆数变化设置T0比较讲究。如果初始温度太低算法刚开始就不接受差解等于白跑退火如果太高前几百次迭代都在随机乱跳浪费计算。一个简单方法是先随机生成几十个解计算目标值标准差把T0设成该标准差的5-10倍。我在示例里用1000是针对20个客户和车辆惩罚100的情况换数据时最好重新标定。车辆数惩罚系数的调整经验是这样跑第一遍时把惩罚设成0看最优目标下的车辆数和距离再把惩罚设成平均每辆车距离的0.5倍重新跑比较结果。示例中距离是一百到几百单位所以惩罚系数取100是比较合理的折中。4.2 一个20客户点的测试案例我在固定随机种子42下跑了上述代码初始解的目标值大约在800左右含车辆惩罚。模拟退火完成后结果大概会收敛到500附近使用4-5辆车纯行驶距离在100-300之间浮动。由于随机种子和邻域操作的关系每次运行结果会有波动但总体距离比随机构造要低20%-40%。测试时有一点值得一提不同邻域操作的随机组合决定算法搜索路径所以同一组参数多跑几次取最好结果比单次记录更有参考价值。这也是元启发式在工程里的常规操作。4.3 收敛曲线怎么看收敛曲线一般分成三个阶段。第一阶段是快速下降温度还很高算法大量接受差解目标值会抖动但整体趋势快速走低。这个阶段实际上是“勘探”在寻找有潜力的盆地。第二阶段是平稳下降温度变低抖动量减小目标值稳步下降这个阶段是在局部区域精修。第三阶段是平台期再降温已经很难改善曲线趋于水平。如果曲线到了很好看的“指数衰减”说明初始温度或降温系数设置合理如果曲线一直在高位抖动且降不下来可能是初始温度太低或每个温度迭代次数不足如果最终目标还是不理想优先检查邻域操作是否被某种模式锁死。5. 常见问题与排查技巧5.1 问题速查表我把实际中容易踩的坑整理成了一张表格都是程序跑起来后最容易报错或结果异常的情况。症状可能原因解决办法未定义函数或变量无法识别函数文件名与函数名不一致将文件保存为对应名称例如evaluateSolution.m初始就报容量不可行某个客户送货需求大于车辆容量检查数据减小需求或更换大容量车辆每次结果完全不同没有固定随机种子起点加rng(seed)便于复现和对比目标值下降缓慢温度下降太快或L太小增大alpha比如从0.99调到0.995路线出现明显交叉且无法优化邻域操作过于单一混合交换、插入和2-opt三种操作车辆数特别多车辆惩罚太低增大车辆数惩罚系数车辆数很少但距离特别远车辆惩罚太高减小惩罚系数重新权衡排查时建议把每轮外层迭代的最优目标值打印出来或者画收敛曲线。如果曲线没有下降趋势多半是参数没有匹配问题规模。5.2 提高求解质量的三个细节第一个细节是初始解不要用随机解硬抗。虽然模拟退火能修正差初始解但一个好的初始解能大幅减少搜索时间。可以用最近邻贪心法生成一个包含所有客户的大序列例如从配送中心出发每次选离上一个客户最近的未访问客户。这样初始距离就比随机解低很多。第二个细节是温度调度可以分段。先快降采样再慢降精修能兼顾收敛速度和精度。简单做法是把alpha设成0.995但前一半迭代用0.999后一半用0.995。第三个细节是容量检查的提前剪枝。在容量检查函数中如果某个路径的送货需求总和超过容量其实可以直接判定不可行因为出发装载量已经超了。这样能少算几次逐点扫描尤其当客户数上百时能省不少时间。我个人的体会是这类问题最难点往往不是算法本身而是把“动态容量检查”这个业务现实转换成代码逻辑。很多现场跑不通的路线都是因为只看了总送货量没注意中途取货把车厢塞满的情况。我自己在调参时会故意设计几个“先送大货、后收大货”的极端数据来测试容量检查是否可靠确认无误后再交给业务方跑真实数据。你如果也想把这套代码用在自己的数据上优先做两件事一是把需求数据格式对齐二是把车辆惩罚和初始温度重新标定一遍。这样从模拟退火算法到Matlab实现整个链路就闭合了。