的模拟退火求解与Matlab实现)
前阵子为了搞定一个带硬性交付时间的配送调度需求我把带时间窗的车辆路径规划问题VRPTW完整啃了一遍最终落地实现用的是模拟退火算法结合 Matlab。市面上的开源代码不少但真正能看懂、能改、能放进论文里的其实不多。这篇博文就把我的完整思路、核心代码结构、参数调优经验以及踩过的坑一次性讲明白适合刚接触 VRP 的学生、考研或竞赛党、以及想快速验证算法的从业者参考。模拟退火的优势在于实现门槛低、对初值不敏感、不容易写崩再加上问题规模控制在几百个节点以内时它拿到高质量可行解的概率非常高。配合时间窗的惩罚处理它能把“先满足时间约束、再优化总里程”的多目标压力转成一个可量化的目标函数这也是它能成为 VRPTW 经典入门算法的原因。这篇内容会从问题建模开始一直讲到 Matlab 代码怎么写、参数怎么定、结果怎么画、论文怎么描述最后附上几个我实际排查过的坑。别把模拟退火想得太玄本质就是一个“热的时候乱跳、冷的时候细调”的随机搜索难的是给这个随机搜索配一套不会跑飞的目标函数和解码逻辑。1. 先把问题本身拆干净什么是带时间窗的车辆路径规划1.1 问题定义和常见约束先看最基础的描述有一个配送中心有一组客户每个客户有坐标位置、货物需求量和允许服务的时间区间也就是时间窗。仓库需要安排若干台车辆从配送中心出发依次服务若干个客户最后回到配送中心。目标一般是最小化总行驶距离或总成本同时满足以下几个约束每辆车有最大载重客户需求量不能超过车辆剩余容量每个客户只能被一辆车服务一次且必须在规定的时间窗内开始服务车辆可以在时间窗开始前到达并等待但晚到通常视为严重违约除非模型允许软时间窗并附加惩罚每辆车从仓库出发后连续作业服务完成后回库。如果去掉时间窗这就是经典的容量约束车辆路径问题也就是 CVRP。而 VRPTW 在 CVRP 基础上多出来的核心难点是“时间窗”这种序列相关的约束。路径的位置顺序会直接影响到达时间晚到、等待、服务时长、行驶时间交错在一起导致解空间的可行域被切割得七零八落。实际物流场景里的时间窗通常分成两类。硬时间窗要求车辆必须在截止时间之前到达否则这条线路不能成立软时间窗则允许稍微迟到但会产生惩罚成本比如客户投诉折算、赔偿费用、下级转运延迟等待。算法上处理软时间窗更容易收敛因为可行域不要求严格断裂而是变成目标函数中的一个大惩罚项。1.2 数学建模和优化目标用比较简明的符号表示这个问题。假设有 n 个客户节点编号从 1 到 n配送中心记作 0用集合 N 表示所有节点。客户 i 的坐标、需求量、最早服务时间、最晚服务时间、服务时长分别记为 xi, yi, qi, ei, li, si。行驶时间矩阵为 t_{ij}车辆数量为 K载重为 Q。决策变量为 x_{ijk}表示车辆 k 是否从节点 i 行驶到节点 j另用变量 a_i 表示车辆到达节点 i 的时间。目标函数可以写成min Z ΣΣΣ t_{ij} * x_{ijk} penalty约束条件大体包括以下几点每个客户只能被访问一次Σ_j x_{ijk} 对于任意客户 i 和被选中的车辆 k 的求和等于1车辆进出守恒每辆车在每个节点的出度和入度相等载重约束每条路径上的累积需求不超过 Q时间窗约束ei ≤ a_i ≤ li如果提前到达则允许等待起终点约束所有车辆从仓库出发并回到仓库。实际编程时我不会字面意义上实现所有约束而是把约束拆成两类。硬性约束如“每个客户访问一次”和“载重不超”直接控制解的结构时间窗以及少量载重越界则放进目标函数里加惩罚。这样的处理方式让邻域搜索中的无效解比例大大降低退火过程也更平滑。1.3 为什么推荐模拟退火这个元启发式VRPTW 是典型的 NP-Hard 问题节点规模一大精确算法就跑不动了。工业场景和毕业论文里常用的方法是各类元启发式遗传算法、禁忌搜索、粒子群、蚁群以及模拟退火。我的经验是模拟退火在“代码可控性”这个维度上非常强几乎没有需要调节的复杂算子核心就是三件事目标函数、邻域动作、温度更新。相比遗传算法模拟退火不需要设计交叉和变异的比例不担心种群早熟也不会因为染色体长度不同出现解码困难。相比禁忌搜索模拟退火的邻域动作可以写得非常简单不需要维护禁忌表、特赦规则和评价函数。对时间窗这种强约束问题禁忌表设计不好反而会越搜越偏。但要注意模拟退火是单点搜索算法没有种群多样性理论上更容易陷入局部最优。因此工程上常常用“多起点”和“加大扰动”来对冲这个问题。后面我会说怎么在参数上做文章既保留退火的稳定性又能尽量接近全局最优。2. 算法设计与实现的关键点2.1 解的表达方式用客户序列生成路径VRPTW 的常见解编码方式有两种。一种是固定车辆数的分段编码比如每个客户直接指定给第几辆车再把每辆车的客户列表排序成路线另一种是客户全排列编码把所有客户排成一个序列然后按容量和时间窗约束切分成若干段。我强烈推荐第二种。原因在于在做邻域搜索时你只需要对一条字符串序列做交换、插入、反转不用时刻顾着车辆编号边界解码的时候再统一拆分这样做代码简单得多而且插入类算子本身就包含“重新分车”的语义。举个例子假设存在 8 个客户排列是 [3, 7, 1, 5, 2, 8, 6, 4]。解码时会从第一个客户开始累加需求和行驶时间如果某个客户加进来之后导致车辆超载或者到达时间已经无法满足时间窗那就在这个位置切开新开一条路线。每个客户会按排列顺序依次访问切割后的代码块各自形成一条从仓库出发并回归仓库的路径。这种编码的优点之一在于你不需要在每一步都直接记录完整路线只需要维护当前排列顺序路径切片完全由约束决定。这个设计大大降低了邻域算子的复杂度。2.2 初始解的重要性很多人直接随机生成一个客户排列当作初始解结果发现后面跑了很久效果都一般。原因不在于模拟退火本身而在于初始解距离可行域太远进入可行区域后温度已经降得很低几乎没有跳出局部极小的能力。我惯用的初始解构造方法是一个贪心插入策略从仓库出发每次都挑选一个“插入后增量最小”的客户。具体逻辑是先在已有路径的各个位置尝试插入当前候选客户计算行驶距离增加量、到达时间延迟、违反时间窗的惩罚增量然后选增量最小的位置插入。如果当前车辆无法承载就开一辆新车。这种初始解比纯随机解的质量高不少通常已经是一个相对合理的路线方案。即使后续退火过程不理想输出结果的下限也有保障。实测在中等规模算例里良好的初始解能帮助最终目标值提升约 10% 到 15%。2.3 邻域搜索算子怎么选模拟退火迭代的核心动作就是“在当前解附近生成一个新解”这个动作决定了算法的探索能力和搜索效率。我做 VRPTW 时最常用的邻域动作有三个。第一是插入算子从当前序列中随机抽出一个客户重新放到另一个随机位置。这个算子不仅是路径内部排序的变化还隐含着改变车辆分配的可能性。如果你只做段内的重新排序可能会把搜索范围限制在固定车辆路线里降低求解质量。第二是交换算子随机选择序列中的两个客户交换它们的位置。这个动作对路径长度和到达时间的影响比较随机升温时可用来做大幅度扰动低温时则会退化成为微调。第三是反转算子类似 2-opt 的思路取序列中随机一片区间把区间内客户顺序倒排。这个动作对长距离交叉线路特别有效可以瞬间消掉一条绕行的路段。实际执行时我会给三个算子分配不同概率通常插入占 0.4交换占 0.3反转占 0.3。每次迭代只执行其中一个算子生成的新序列与原序列只差一个局部结构目标函数变化也相对平滑有利于退火过程稳步收敛。2.4 目标函数与时间窗惩罚的设计目标函数怎么写直接影响搜索方向。如果只算总行驶距离忽略时间窗算法大概率倾向于把位置相近的客户放在一起但完全无视时间紧迫性最后产出大量不可行解。我采用的方式是加权求和总成本 总行驶距离 等待时间成本 时间窗延迟惩罚 载重超限惩罚。行驶距离是主要优化目标等待时间乘以一个较小系数时间窗延迟惩罚则用大权重控制。时间窗惩罚的设计尤其关键。硬时间窗场景下若晚到超过阈值就直接判为不可行那么目标值就会变得非常突兀搜索容易“碰壁”。更稳妥的做法是设置一个很大的惩罚系数 P并在退火早期允许少量违约使算法能穿过不可行区域在后期再逐渐强化时间窗约束。这实际上是一种变惩罚策略和退火过程天然契合。载重惩罚可以很简单车辆超载量和权重系数的成绩加到目标函数里。由于邻域动作大多是小幅度改动超载情况一般不会一直存在惩罚项更多起到引导角色。3. Matlab 工程实现细节与核心代码3.1 数据组织和距离矩阵在 Matlab 写 VRPTW第一步是把数据整理成标准化结构体。我经常这样组织数据nodesn×6 矩阵分别保存 x 坐标、y 坐标、需求量、最早开始时间、最晚开始时间、服务时长distMatrixn×n 对称矩阵元素是节点之间的欧氏距离或实际行驶距离travelTimeMatrixn×n 矩阵代表行驶耗时vehicleCapacity每辆车最大载重numCustomers客户总数。计算距离矩阵时要注意比例问题。如果坐标单位是公里行驶速度按 40km/h 算那么时间窗也要用相同的分钟或小时单位。否则可能出现“距离满足、时间完全不匹配”的荒谬局面。代码上通常直接这样生成coords nodes(:, 1:2); numNodes size(nodes, 1); distMatrix zeros(numNodes, numNodes); for i 1:numNodes for j 1:numNodes distMatrix(i, j) sqrt((coords(i,1)-coords(j,1))^2 (coords(i,2)-coords(j,2))^2); end end travelTimeMatrix distMatrix / speed;务必要在代码注释里标明单位换算。我见过不少案例坐标是千米、时间窗是分钟结果界面上下乱跳。后面这些误差极难排查因为它不是程序崩溃而是所有结果都很“合理但差”。3.2 解码与可行性检查函数带时间窗的解码听起来复杂但核心代码其实就是一条循环逐客户累加时间并判断新插入客户后是否会超载或超时。标准化解码函数的形式大致是function [routes, totalCost, feasible] decodeSequence(seq, nodes, distMatrix, travelTimeMatrix, cap, speed, timePenalty, capacityPenalty) % 这是一个示意结构实际运行时需要处理仓库节点索引 routes {}; totalCost 0; currentLoad 0; currentTime 0; currentRoute []; for k 1:length(seq) nodeIdx seq(k); if isempty(currentRoute) travelDist distMatrix(1, nodeIdx); travelTime travelTimeMatrix(1, nodeIdx); arrive currentTime travelTime; else lastIdx currentRoute(end); travelDist distMatrix(lastIdx, nodeIdx); travelTime travelTimeMatrix(lastIdx, nodeIdx); arrive currentTime travelTime; end wait max(0, nodes(nodeIdx,4) - arrive); serviceStart max(arrive, nodes(nodeIdx,4)); newLoad currentLoad nodes(nodeIdx,3); late max(0, serviceStart - nodes(nodeIdx,5)); % 如果超载或严重超时强制切分新路线 if newLoad cap || (late 1e-6 当前路线已有多个客户) % 先收尾当前路线再开新路线 routes{end1} currentRoute; currentRoute []; currentLoad 0; currentTime 0; continue; end currentRoute(end1) nodeIdx; currentLoad newLoad; currentTime serviceStart nodes(nodeIdx,6); totalCost totalCost travelDist wait * 0.1 late * timePenalty; end if ~isempty(currentRoute) routes{end1} currentRoute; end end上面的代码是示意真实使用还要把“回到仓库”的距离和耗时补上。这里的关键点是服务开始时间的计算如果车辆提前到达它要等到最早服务时间才开始如果晚到服务开始时间等于到达时间。arrive应该等于上一节点的服务结束时间加上行驶时间而不是上一节点到达时间直接加上行驶时间。还有一点很多新手会漏掉仓库节点本身没有时间窗但车辆从仓库出发的时间一般设定为 0。当一条新路线开始后currentTime要重新归零不能沿用上一条路线的返回时间。3.3 模拟退火主循环主循环的结构非常固定。先初始化温度和当前解然后在每个温度下进行固定次数的迭代。每次迭代随机挑一个邻域动作生成新序列解码出目标值计算增量 ΔE并按 Metropolis 准则决定是否接受。function [bestSequence, bestCost, costHistory] saSolve(nodes, distMatrix, travelTimeMatrix, vehicleCapacity, opts) seq initialSolution(nodes, distMatrix, travelTimeMatrix, vehicleCapacity); currentSeq seq; currentCost evaluateSequence(seq); bestSeq currentSeq; bestCost currentCost; T opts.T0; costHistory []; while T opts.Tend for inner 1:opts.innerIter newSeq neighborOperator(currentSeq); newCost evaluateSequence(newSeq); delta newCost - currentCost; if delta 0 || rand exp(-delta/T) currentSeq newSeq; currentCost newCost; end if currentCost bestCost bestSeq currentSeq; bestCost currentCost; end end T T * opts.alpha; costHistory(end1) bestCost; end end这里的关键参数有三个T0、Tend和alpha。innerIter是每个温度下的迭代次数一般设置成节点数量的线性或平方函数。如果问题有 100 个客户我会选innerIter在 500 到 1000 之间如果节点数到了 300这个数字可能要提升到 2000 以上。3.4 初始解生成与邻域算子实现初始解生成我用的是贪心插入。Matlab 里可以写成这样function route initialSolution(data, distMatrix, travelTimeMatrix, cap) customerList 1:data.numCustomers; route []; while ~isempty(customerList) bestNode -1; bestIncrease inf; bestPos -1; for i 1:length(customerList) cand customerList(i); % 计算插入到每个位置后的增量成本 for pos 1:length(route)1 inc computeInsertionPenalty(route, cand, pos, data, distMatrix, travelTimeMatrix, cap); if inc bestIncrease bestIncrease inc; bestNode cand; bestPos pos; end end end route [route(1:bestPos-1), bestNode, route(bestPos:end)]; customerList(customerList bestNode) []; end end邻域算子的实现就是在序列上做索引操作。插入算子抽出某个位置的元素再插入新位置交换算子直接newSeq(i) seq(j); newSeq(j) seq(i)反转算子则是对中间段使用fliplr。这里要记住一个细节生成新解后不要修改原序列而是复制一份再操作。如果直接原地修改后面计算 delta 用到的“原目标值”和“新目标值”就不是同一基线上的对比了。4. 参数调优与实验效果分析4.1 退火温度该怎么初始化和衰减初始温度 T0 直接影响算法在前面阶段的“容错能力”。如果 T0 太小算法一开始就不敢接受劣解搜索很快退化成局部贪婪如果 T0 太大前期大量接受劣解虽然探索范围大了但会把本来就优质的初始解折腾得面目全非。一个经验公式是先随机生成一组邻域解统计目标值变化量的标准差 σ然后取 T0 初始允许接受概率所需的倍数关系。更粗暴一点直接把 T0 设成初始目标值的 10% 到 30%。比如初始目标值是 800那 T0 可以设在 80 到 240 之间。后期温度降到 Tend 1 或者 0.1 就足够稳定因为此时每个劣解被接受的概率已经微乎其微。衰减系数 alpha 的选取与迭代次数密切相关。alpha 越小降温越快搜索越快但容易陷入局部最优alpha 越大降温越慢结果通常更好但时间更长。我常用的范围是 0.90 到 0.98 之间的一档。4.2 参数参考表和实测对比为了让大家好抄作业我把自己在 100 个客户规模算例上常用的参数整理成一个表。参数建议取值调高效果调低效果T050 ~ 200探索更广前期震荡大收敛快容易局部最优alpha0.90 ~ 0.98结果更稳耗时增加降温快结果波动大innerIter500 ~ 2000求解更彻底耗时线性上升迭代不足解质量下降插入算子概率0.3 ~ 0.4路径分配变化更灵活搜索偏局部交换算子概率0.2 ~ 0.4大范围扰动更多收敛相对平稳反转算子概率0.2 ~ 0.4消除交叉路径长距离绕行难改善时间惩罚系数100 ~ 1000约束硬可行性强允许晚到目标值偏低需要说明的是惩罚系数不是一个拍脑袋的数它必须明显大于“正常行驶一段路的成本”。如果时间窗违约一次的总成本是 50而车辆开一单位距离的成本是 1那么算法为了规避违约宁愿多绕 50 单位路程这是合理行为。4.3 收敛判断和终止条件判断收敛的方法有两种。一种直观的办法是记录目标值历史曲线如果在连续多个温度段内最优解的提升幅度小于 0.5%就可以提前终止。另一种办法是固定总运行时间比如限制整体循环 60 秒或 200 秒到时间直接输出当前最优解。我在实际求解时倾向使用“温度下限 无改进最大温度数”的组合终止条件。也就是即使 T 还没降到 Tend如果连续 30 个温度阶段都没有找到更优解就直接跳出。这样可以避免退火后期纯粹在最优解附近白耗时间对大批量算例尤其有效。不过论文实验里为了画平滑的收敛曲线我还是会保留完整退火过程并把 pre-termination 逻辑单独注释掉。4.4 多起点和重启策略模拟退火的单点性质决定了它容易陷入局部最优。最简单的补救办法就是跑多次随机设置不同初始解或者温度重启。比如每完成一轮完整退火后把当前最优解稍微扰动再把温度弹回 T0重新开始搜索这样做相当于做了一次粗糙的“多起点”搜索。实测结果表明对 100 个客户级别的问题跑 3 次独立退火取最优值通常比单次长迭代效果更稳定。这是因为每次的运行时长控制在可接受范围而不同起点能覆盖到不同的搜索区域组合起来比一条热链贯穿到底更稳妥。5. 常见问题与排查技巧实录5.1 解总是晚到但目标值却很漂亮这个问题我遇到太多次了。现象是输出的总距离明显低于合理值但拆开每条路线一看不少客户到达时间晚于最晚服务时间。原因几乎总是目标函数里时间窗惩罚权重开得太小等于给算法“吃了颗定心丸”让它觉得违约无所谓。排查步骤第一步检查解码函数中对late的统计第二步检查惩罚系数是否比“一单位总距离成本”高至少一个数量级第三步检查时间窗单位对不对。如果坐标单位是公里、速度是 40km/h那时间单位应该是小时而不是分钟。我曾经有过整套代码跑了半小时结果全是“晚到”最后查出来是速度写成了每小时公里数第一遍单位没对齐。5.2 邻域操作破坏了客户访问唯一性插入、交换、反转算子都是在客户序列上做操作理论上不会改变客户集合本身。但如果你在某一步用了带重复或缺失的随机索引或者不小心对序列seq直接赋值给newSeq就会出现“新解里某个客户出现两次、某客户消失”的情况。建议是在所有邻域操作里都写一个防呆检查assert(length(unique(newSeq)) length(newSeq), 邻域操作生成了非法序列);这个小断言在调试阶段非常有用。有时候明明目标值改进很多但跑着跑着就出现“客户数凭空少一个”的问题根源就是索引打架。5.3 退火过程快速收敛但结果很差如果目标值从第一个温度段开始就快速下降之后几乎不动大概率是 T0 设得太低或者衰减系数 alpha 太小。可以增加一个“扰动观察”在初始温度阶段随机接受多个劣解观察目标值是否能上升之后再下降。如果几次劣解直接全被拒绝说明初始温度不够需要把 T0 调大。另一种情况是邻域算子的“步长”太小每次只交换相邻客户导致搜索只能在原路径附近微调。这时可以临时把反转算子加大或者对序列中间段做更大幅度的重排。步长和温度是两个不同层面的控制维度二者不能用一个参数替代。5.4 Matlab 画图时路径乱连线画车辆路径图一般用 plot 函数把每辆车的路径节点按顺序连起来。如果解码后返回的 route 列表顺序不对画出来就会“飞线”。一个关键检查点是解码函数切分路线后每辆车的路线结尾是否补上了仓库节点索引 1。如果缺少回程节点图会缺一条边看起来像所有路径都直接悬空。画图代码示例figure; hold on; plot(nodes(2:end,1), nodes(2:end,2), bo, MarkerSize, 8); plot(nodes(1,1), nodes(1,2), rs, MarkerSize, 12); colormapLines lines(length(routes)); for i 1:length(routes) r [1, routes{i}, 1]; plot(nodes(r,1), nodes(r,2), -, Color, colormapLines(i,:), LineWidth, 1.2); end画完再看一眼如果某条边跨越了整个图面大概率是路线序列里少了或多了某个节点。5.5 数值稳定性问题时间窗服务时间和等待时间的小数点误差在某些判断里可能会造成问题。比如判断“到达时间是否晚于最晚服务时间”由于浮点数的表达误差一个精确位于边界上的数值可能被判为晚到从而影响目标函数值。处理办法很简单在比较时加上一个小阈值例如if arrive nodes(nodeIdx,5) 1e-6 % 计为晚到 end这个 1e-6 的松弛量不会产生实质影响但能避免很多“莫名其妙被切路线”的问题。6. 从代码到论文实验和结果呈现怎么写6.1 实验设计和算例准备论文里最害怕的问题是“只有一张结果表没有实验过程”。写 VRPTW 求解实验时至少要交代以下几个方面测试算例的规模节点数量、客户数量、时间窗宽度、车辆容量限制参数设定初始温度、终止温度、衰减系数、内层迭代次数、邻域算子概率对比基线是使用人工经验路径、贪心算法、其他启发式还是公开已知最优解评价指标总行驶里程、总耗时、平均等待时间、超时次数、车辆使用数。建议在论文中放一张“算例参数表”列出所有节点的坐标范围、需求量分布、时间窗分布。这样读者才能判断你的实验是否可信。6.2 结果可视化与收敛曲线收敛曲线图是模拟退火论文的标配。横轴可以选“迭代次数”或“温度下降阶段”纵轴选当前最优目标值。不要只画一条最终最优值的曲线那样看不出算法过程。常规做法是同时画两条线一条是当前解在每个温度段的目标值另一条是历史最优目标值。两条线之间的差距能直观体现退火搜索的波动性。路径规划图建议按不同车辆用不同颜色区分并标注仓库位置。如果时间窗比较紧还可以在节点旁边标注“最早开始时间/最晚开始时间”帮助审稿人理解算法为什么给某些客户安排了特殊顺序。6.3 对比实验怎么更有说服力单纯的模拟退火结果很难让审稿人眼前一亮。更稳妥的策略是设置三组对比第一组是纯贪心算法结果第二组是模拟退火使用随机初始解的运行结果第三组是模拟退火使用贪心初始解加调优参数的运行结果。通过三组数值对比可以证明算法的改进不只是来源于元启发式自身更来源于初始解和惩罚参数设计的组合优势。如果时间允许可以再加一组“不同惩罚系数下的结果对比”说明过大惩罚会让解偏保守过小惩罚会让时间窗失去意义以此论证你选择的惩罚系数有实验依据。不要只列最佳结果也要给出平均数、最差值、标准差这会显著提升论文的实验可信度。6.4 代码附录的组织很多论文要求附代码但很多人直接贴整个脚本阅读体验很差。更理想的做法是拆分成几个有注释的函数块数据导入模块、初始解生成模块、邻域搜索模块、解码评估模块、主循环模块、可视化模块。每个模块前用一段文字说明输入输出含义并把参数设置独立成一个小节方便读者调整。Matlab 代码里尽量加中文注释并写明运行环境版本。如果代码中使用了自定义矩阵或结构体字段最好在附录开头写一个字段说明表能省掉读者大量对字段名发愁的时间。一些实操层面的经验心得我自己的体会是模拟退火求解 VRPTW 的难点不在算法本身而在“时间窗检查”和“路径解码”这两个容易犯迷糊的地方。很多开源代码喜欢把所有逻辑都堆在主循环里看起来很短但改一个参数就要全局排查远不如我上面这种“解码函数 邻域算子 主循环”的结构清晰。最后再分享一个小技巧求解之前先把仓库坐标和数据规模在图上画一遍直观感受一下客户分布。如果客户太分散时间窗又窄那无论算法调得多好也不可能达到理想中的“几辆车全走完”。先通过画图确认问题本身是不是“有资格难”再决定要不要用更复杂的算法。别一上来就写模拟退火主循环跑完之后才发现数据存在严重单位或坐标错误那样一天时间就白费了。
前阵子为了搞定一个带硬性交付时间的配送调度需求我把带时间窗的车辆路径规划问题VRPTW完整啃了一遍最终落地实现用的是模拟退火算法结合 Matlab。市面上的开源代码不少但真正能看懂、能改、能放进论文里的其实不多。这篇博文就把我的完整思路、核心代码结构、参数调优经验以及踩过的坑一次性讲明白适合刚接触 VRP 的学生、考研或竞赛党、以及想快速验证算法的从业者参考。模拟退火的优势在于实现门槛低、对初值不敏感、不容易写崩再加上问题规模控制在几百个节点以内时它拿到高质量可行解的概率非常高。配合时间窗的惩罚处理它能把“先满足时间约束、再优化总里程”的多目标压力转成一个可量化的目标函数这也是它能成为 VRPTW 经典入门算法的原因。这篇内容会从问题建模开始一直讲到 Matlab 代码怎么写、参数怎么定、结果怎么画、论文怎么描述最后附上几个我实际排查过的坑。别把模拟退火想得太玄本质就是一个“热的时候乱跳、冷的时候细调”的随机搜索难的是给这个随机搜索配一套不会跑飞的目标函数和解码逻辑。1. 先把问题本身拆干净什么是带时间窗的车辆路径规划1.1 问题定义和常见约束先看最基础的描述有一个配送中心有一组客户每个客户有坐标位置、货物需求量和允许服务的时间区间也就是时间窗。仓库需要安排若干台车辆从配送中心出发依次服务若干个客户最后回到配送中心。目标一般是最小化总行驶距离或总成本同时满足以下几个约束每辆车有最大载重客户需求量不能超过车辆剩余容量每个客户只能被一辆车服务一次且必须在规定的时间窗内开始服务车辆可以在时间窗开始前到达并等待但晚到通常视为严重违约除非模型允许软时间窗并附加惩罚每辆车从仓库出发后连续作业服务完成后回库。如果去掉时间窗这就是经典的容量约束车辆路径问题也就是 CVRP。而 VRPTW 在 CVRP 基础上多出来的核心难点是“时间窗”这种序列相关的约束。路径的位置顺序会直接影响到达时间晚到、等待、服务时长、行驶时间交错在一起导致解空间的可行域被切割得七零八落。实际物流场景里的时间窗通常分成两类。硬时间窗要求车辆必须在截止时间之前到达否则这条线路不能成立软时间窗则允许稍微迟到但会产生惩罚成本比如客户投诉折算、赔偿费用、下级转运延迟等待。算法上处理软时间窗更容易收敛因为可行域不要求严格断裂而是变成目标函数中的一个大惩罚项。1.2 数学建模和优化目标用比较简明的符号表示这个问题。假设有 n 个客户节点编号从 1 到 n配送中心记作 0用集合 N 表示所有节点。客户 i 的坐标、需求量、最早服务时间、最晚服务时间、服务时长分别记为 xi, yi, qi, ei, li, si。行驶时间矩阵为 t_{ij}车辆数量为 K载重为 Q。决策变量为 x_{ijk}表示车辆 k 是否从节点 i 行驶到节点 j另用变量 a_i 表示车辆到达节点 i 的时间。目标函数可以写成min Z ΣΣΣ t_{ij} * x_{ijk} penalty约束条件大体包括以下几点每个客户只能被访问一次Σ_j x_{ijk} 对于任意客户 i 和被选中的车辆 k 的求和等于1车辆进出守恒每辆车在每个节点的出度和入度相等载重约束每条路径上的累积需求不超过 Q时间窗约束ei ≤ a_i ≤ li如果提前到达则允许等待起终点约束所有车辆从仓库出发并回到仓库。实际编程时我不会字面意义上实现所有约束而是把约束拆成两类。硬性约束如“每个客户访问一次”和“载重不超”直接控制解的结构时间窗以及少量载重越界则放进目标函数里加惩罚。这样的处理方式让邻域搜索中的无效解比例大大降低退火过程也更平滑。1.3 为什么推荐模拟退火这个元启发式VRPTW 是典型的 NP-Hard 问题节点规模一大精确算法就跑不动了。工业场景和毕业论文里常用的方法是各类元启发式遗传算法、禁忌搜索、粒子群、蚁群以及模拟退火。我的经验是模拟退火在“代码可控性”这个维度上非常强几乎没有需要调节的复杂算子核心就是三件事目标函数、邻域动作、温度更新。相比遗传算法模拟退火不需要设计交叉和变异的比例不担心种群早熟也不会因为染色体长度不同出现解码困难。相比禁忌搜索模拟退火的邻域动作可以写得非常简单不需要维护禁忌表、特赦规则和评价函数。对时间窗这种强约束问题禁忌表设计不好反而会越搜越偏。但要注意模拟退火是单点搜索算法没有种群多样性理论上更容易陷入局部最优。因此工程上常常用“多起点”和“加大扰动”来对冲这个问题。后面我会说怎么在参数上做文章既保留退火的稳定性又能尽量接近全局最优。2. 算法设计与实现的关键点2.1 解的表达方式用客户序列生成路径VRPTW 的常见解编码方式有两种。一种是固定车辆数的分段编码比如每个客户直接指定给第几辆车再把每辆车的客户列表排序成路线另一种是客户全排列编码把所有客户排成一个序列然后按容量和时间窗约束切分成若干段。我强烈推荐第二种。原因在于在做邻域搜索时你只需要对一条字符串序列做交换、插入、反转不用时刻顾着车辆编号边界解码的时候再统一拆分这样做代码简单得多而且插入类算子本身就包含“重新分车”的语义。举个例子假设存在 8 个客户排列是 [3, 7, 1, 5, 2, 8, 6, 4]。解码时会从第一个客户开始累加需求和行驶时间如果某个客户加进来之后导致车辆超载或者到达时间已经无法满足时间窗那就在这个位置切开新开一条路线。每个客户会按排列顺序依次访问切割后的代码块各自形成一条从仓库出发并回归仓库的路径。这种编码的优点之一在于你不需要在每一步都直接记录完整路线只需要维护当前排列顺序路径切片完全由约束决定。这个设计大大降低了邻域算子的复杂度。2.2 初始解的重要性很多人直接随机生成一个客户排列当作初始解结果发现后面跑了很久效果都一般。原因不在于模拟退火本身而在于初始解距离可行域太远进入可行区域后温度已经降得很低几乎没有跳出局部极小的能力。我惯用的初始解构造方法是一个贪心插入策略从仓库出发每次都挑选一个“插入后增量最小”的客户。具体逻辑是先在已有路径的各个位置尝试插入当前候选客户计算行驶距离增加量、到达时间延迟、违反时间窗的惩罚增量然后选增量最小的位置插入。如果当前车辆无法承载就开一辆新车。这种初始解比纯随机解的质量高不少通常已经是一个相对合理的路线方案。即使后续退火过程不理想输出结果的下限也有保障。实测在中等规模算例里良好的初始解能帮助最终目标值提升约 10% 到 15%。2.3 邻域搜索算子怎么选模拟退火迭代的核心动作就是“在当前解附近生成一个新解”这个动作决定了算法的探索能力和搜索效率。我做 VRPTW 时最常用的邻域动作有三个。第一是插入算子从当前序列中随机抽出一个客户重新放到另一个随机位置。这个算子不仅是路径内部排序的变化还隐含着改变车辆分配的可能性。如果你只做段内的重新排序可能会把搜索范围限制在固定车辆路线里降低求解质量。第二是交换算子随机选择序列中的两个客户交换它们的位置。这个动作对路径长度和到达时间的影响比较随机升温时可用来做大幅度扰动低温时则会退化成为微调。第三是反转算子类似 2-opt 的思路取序列中随机一片区间把区间内客户顺序倒排。这个动作对长距离交叉线路特别有效可以瞬间消掉一条绕行的路段。实际执行时我会给三个算子分配不同概率通常插入占 0.4交换占 0.3反转占 0.3。每次迭代只执行其中一个算子生成的新序列与原序列只差一个局部结构目标函数变化也相对平滑有利于退火过程稳步收敛。2.4 目标函数与时间窗惩罚的设计目标函数怎么写直接影响搜索方向。如果只算总行驶距离忽略时间窗算法大概率倾向于把位置相近的客户放在一起但完全无视时间紧迫性最后产出大量不可行解。我采用的方式是加权求和总成本 总行驶距离 等待时间成本 时间窗延迟惩罚 载重超限惩罚。行驶距离是主要优化目标等待时间乘以一个较小系数时间窗延迟惩罚则用大权重控制。时间窗惩罚的设计尤其关键。硬时间窗场景下若晚到超过阈值就直接判为不可行那么目标值就会变得非常突兀搜索容易“碰壁”。更稳妥的做法是设置一个很大的惩罚系数 P并在退火早期允许少量违约使算法能穿过不可行区域在后期再逐渐强化时间窗约束。这实际上是一种变惩罚策略和退火过程天然契合。载重惩罚可以很简单车辆超载量和权重系数的成绩加到目标函数里。由于邻域动作大多是小幅度改动超载情况一般不会一直存在惩罚项更多起到引导角色。3. Matlab 工程实现细节与核心代码3.1 数据组织和距离矩阵在 Matlab 写 VRPTW第一步是把数据整理成标准化结构体。我经常这样组织数据nodesn×6 矩阵分别保存 x 坐标、y 坐标、需求量、最早开始时间、最晚开始时间、服务时长distMatrixn×n 对称矩阵元素是节点之间的欧氏距离或实际行驶距离travelTimeMatrixn×n 矩阵代表行驶耗时vehicleCapacity每辆车最大载重numCustomers客户总数。计算距离矩阵时要注意比例问题。如果坐标单位是公里行驶速度按 40km/h 算那么时间窗也要用相同的分钟或小时单位。否则可能出现“距离满足、时间完全不匹配”的荒谬局面。代码上通常直接这样生成coords nodes(:, 1:2); numNodes size(nodes, 1); distMatrix zeros(numNodes, numNodes); for i 1:numNodes for j 1:numNodes distMatrix(i, j) sqrt((coords(i,1)-coords(j,1))^2 (coords(i,2)-coords(j,2))^2); end end travelTimeMatrix distMatrix / speed;务必要在代码注释里标明单位换算。我见过不少案例坐标是千米、时间窗是分钟结果界面上下乱跳。后面这些误差极难排查因为它不是程序崩溃而是所有结果都很“合理但差”。3.2 解码与可行性检查函数带时间窗的解码听起来复杂但核心代码其实就是一条循环逐客户累加时间并判断新插入客户后是否会超载或超时。标准化解码函数的形式大致是function [routes, totalCost, feasible] decodeSequence(seq, nodes, distMatrix, travelTimeMatrix, cap, speed, timePenalty, capacityPenalty) % 这是一个示意结构实际运行时需要处理仓库节点索引 routes {}; totalCost 0; currentLoad 0; currentTime 0; currentRoute []; for k 1:length(seq) nodeIdx seq(k); if isempty(currentRoute) travelDist distMatrix(1, nodeIdx); travelTime travelTimeMatrix(1, nodeIdx); arrive currentTime travelTime; else lastIdx currentRoute(end); travelDist distMatrix(lastIdx, nodeIdx); travelTime travelTimeMatrix(lastIdx, nodeIdx); arrive currentTime travelTime; end wait max(0, nodes(nodeIdx,4) - arrive); serviceStart max(arrive, nodes(nodeIdx,4)); newLoad currentLoad nodes(nodeIdx,3); late max(0, serviceStart - nodes(nodeIdx,5)); % 如果超载或严重超时强制切分新路线 if newLoad cap || (late 1e-6 当前路线已有多个客户) % 先收尾当前路线再开新路线 routes{end1} currentRoute; currentRoute []; currentLoad 0; currentTime 0; continue; end currentRoute(end1) nodeIdx; currentLoad newLoad; currentTime serviceStart nodes(nodeIdx,6); totalCost totalCost travelDist wait * 0.1 late * timePenalty; end if ~isempty(currentRoute) routes{end1} currentRoute; end end上面的代码是示意真实使用还要把“回到仓库”的距离和耗时补上。这里的关键点是服务开始时间的计算如果车辆提前到达它要等到最早服务时间才开始如果晚到服务开始时间等于到达时间。arrive应该等于上一节点的服务结束时间加上行驶时间而不是上一节点到达时间直接加上行驶时间。还有一点很多新手会漏掉仓库节点本身没有时间窗但车辆从仓库出发的时间一般设定为 0。当一条新路线开始后currentTime要重新归零不能沿用上一条路线的返回时间。3.3 模拟退火主循环主循环的结构非常固定。先初始化温度和当前解然后在每个温度下进行固定次数的迭代。每次迭代随机挑一个邻域动作生成新序列解码出目标值计算增量 ΔE并按 Metropolis 准则决定是否接受。function [bestSequence, bestCost, costHistory] saSolve(nodes, distMatrix, travelTimeMatrix, vehicleCapacity, opts) seq initialSolution(nodes, distMatrix, travelTimeMatrix, vehicleCapacity); currentSeq seq; currentCost evaluateSequence(seq); bestSeq currentSeq; bestCost currentCost; T opts.T0; costHistory []; while T opts.Tend for inner 1:opts.innerIter newSeq neighborOperator(currentSeq); newCost evaluateSequence(newSeq); delta newCost - currentCost; if delta 0 || rand exp(-delta/T) currentSeq newSeq; currentCost newCost; end if currentCost bestCost bestSeq currentSeq; bestCost currentCost; end end T T * opts.alpha; costHistory(end1) bestCost; end end这里的关键参数有三个T0、Tend和alpha。innerIter是每个温度下的迭代次数一般设置成节点数量的线性或平方函数。如果问题有 100 个客户我会选innerIter在 500 到 1000 之间如果节点数到了 300这个数字可能要提升到 2000 以上。3.4 初始解生成与邻域算子实现初始解生成我用的是贪心插入。Matlab 里可以写成这样function route initialSolution(data, distMatrix, travelTimeMatrix, cap) customerList 1:data.numCustomers; route []; while ~isempty(customerList) bestNode -1; bestIncrease inf; bestPos -1; for i 1:length(customerList) cand customerList(i); % 计算插入到每个位置后的增量成本 for pos 1:length(route)1 inc computeInsertionPenalty(route, cand, pos, data, distMatrix, travelTimeMatrix, cap); if inc bestIncrease bestIncrease inc; bestNode cand; bestPos pos; end end end route [route(1:bestPos-1), bestNode, route(bestPos:end)]; customerList(customerList bestNode) []; end end邻域算子的实现就是在序列上做索引操作。插入算子抽出某个位置的元素再插入新位置交换算子直接newSeq(i) seq(j); newSeq(j) seq(i)反转算子则是对中间段使用fliplr。这里要记住一个细节生成新解后不要修改原序列而是复制一份再操作。如果直接原地修改后面计算 delta 用到的“原目标值”和“新目标值”就不是同一基线上的对比了。4. 参数调优与实验效果分析4.1 退火温度该怎么初始化和衰减初始温度 T0 直接影响算法在前面阶段的“容错能力”。如果 T0 太小算法一开始就不敢接受劣解搜索很快退化成局部贪婪如果 T0 太大前期大量接受劣解虽然探索范围大了但会把本来就优质的初始解折腾得面目全非。一个经验公式是先随机生成一组邻域解统计目标值变化量的标准差 σ然后取 T0 初始允许接受概率所需的倍数关系。更粗暴一点直接把 T0 设成初始目标值的 10% 到 30%。比如初始目标值是 800那 T0 可以设在 80 到 240 之间。后期温度降到 Tend 1 或者 0.1 就足够稳定因为此时每个劣解被接受的概率已经微乎其微。衰减系数 alpha 的选取与迭代次数密切相关。alpha 越小降温越快搜索越快但容易陷入局部最优alpha 越大降温越慢结果通常更好但时间更长。我常用的范围是 0.90 到 0.98 之间的一档。4.2 参数参考表和实测对比为了让大家好抄作业我把自己在 100 个客户规模算例上常用的参数整理成一个表。参数建议取值调高效果调低效果T050 ~ 200探索更广前期震荡大收敛快容易局部最优alpha0.90 ~ 0.98结果更稳耗时增加降温快结果波动大innerIter500 ~ 2000求解更彻底耗时线性上升迭代不足解质量下降插入算子概率0.3 ~ 0.4路径分配变化更灵活搜索偏局部交换算子概率0.2 ~ 0.4大范围扰动更多收敛相对平稳反转算子概率0.2 ~ 0.4消除交叉路径长距离绕行难改善时间惩罚系数100 ~ 1000约束硬可行性强允许晚到目标值偏低需要说明的是惩罚系数不是一个拍脑袋的数它必须明显大于“正常行驶一段路的成本”。如果时间窗违约一次的总成本是 50而车辆开一单位距离的成本是 1那么算法为了规避违约宁愿多绕 50 单位路程这是合理行为。4.3 收敛判断和终止条件判断收敛的方法有两种。一种直观的办法是记录目标值历史曲线如果在连续多个温度段内最优解的提升幅度小于 0.5%就可以提前终止。另一种办法是固定总运行时间比如限制整体循环 60 秒或 200 秒到时间直接输出当前最优解。我在实际求解时倾向使用“温度下限 无改进最大温度数”的组合终止条件。也就是即使 T 还没降到 Tend如果连续 30 个温度阶段都没有找到更优解就直接跳出。这样可以避免退火后期纯粹在最优解附近白耗时间对大批量算例尤其有效。不过论文实验里为了画平滑的收敛曲线我还是会保留完整退火过程并把 pre-termination 逻辑单独注释掉。4.4 多起点和重启策略模拟退火的单点性质决定了它容易陷入局部最优。最简单的补救办法就是跑多次随机设置不同初始解或者温度重启。比如每完成一轮完整退火后把当前最优解稍微扰动再把温度弹回 T0重新开始搜索这样做相当于做了一次粗糙的“多起点”搜索。实测结果表明对 100 个客户级别的问题跑 3 次独立退火取最优值通常比单次长迭代效果更稳定。这是因为每次的运行时长控制在可接受范围而不同起点能覆盖到不同的搜索区域组合起来比一条热链贯穿到底更稳妥。5. 常见问题与排查技巧实录5.1 解总是晚到但目标值却很漂亮这个问题我遇到太多次了。现象是输出的总距离明显低于合理值但拆开每条路线一看不少客户到达时间晚于最晚服务时间。原因几乎总是目标函数里时间窗惩罚权重开得太小等于给算法“吃了颗定心丸”让它觉得违约无所谓。排查步骤第一步检查解码函数中对late的统计第二步检查惩罚系数是否比“一单位总距离成本”高至少一个数量级第三步检查时间窗单位对不对。如果坐标单位是公里、速度是 40km/h那时间单位应该是小时而不是分钟。我曾经有过整套代码跑了半小时结果全是“晚到”最后查出来是速度写成了每小时公里数第一遍单位没对齐。5.2 邻域操作破坏了客户访问唯一性插入、交换、反转算子都是在客户序列上做操作理论上不会改变客户集合本身。但如果你在某一步用了带重复或缺失的随机索引或者不小心对序列seq直接赋值给newSeq就会出现“新解里某个客户出现两次、某客户消失”的情况。建议是在所有邻域操作里都写一个防呆检查assert(length(unique(newSeq)) length(newSeq), 邻域操作生成了非法序列);这个小断言在调试阶段非常有用。有时候明明目标值改进很多但跑着跑着就出现“客户数凭空少一个”的问题根源就是索引打架。5.3 退火过程快速收敛但结果很差如果目标值从第一个温度段开始就快速下降之后几乎不动大概率是 T0 设得太低或者衰减系数 alpha 太小。可以增加一个“扰动观察”在初始温度阶段随机接受多个劣解观察目标值是否能上升之后再下降。如果几次劣解直接全被拒绝说明初始温度不够需要把 T0 调大。另一种情况是邻域算子的“步长”太小每次只交换相邻客户导致搜索只能在原路径附近微调。这时可以临时把反转算子加大或者对序列中间段做更大幅度的重排。步长和温度是两个不同层面的控制维度二者不能用一个参数替代。5.4 Matlab 画图时路径乱连线画车辆路径图一般用 plot 函数把每辆车的路径节点按顺序连起来。如果解码后返回的 route 列表顺序不对画出来就会“飞线”。一个关键检查点是解码函数切分路线后每辆车的路线结尾是否补上了仓库节点索引 1。如果缺少回程节点图会缺一条边看起来像所有路径都直接悬空。画图代码示例figure; hold on; plot(nodes(2:end,1), nodes(2:end,2), bo, MarkerSize, 8); plot(nodes(1,1), nodes(1,2), rs, MarkerSize, 12); colormapLines lines(length(routes)); for i 1:length(routes) r [1, routes{i}, 1]; plot(nodes(r,1), nodes(r,2), -, Color, colormapLines(i,:), LineWidth, 1.2); end画完再看一眼如果某条边跨越了整个图面大概率是路线序列里少了或多了某个节点。5.5 数值稳定性问题时间窗服务时间和等待时间的小数点误差在某些判断里可能会造成问题。比如判断“到达时间是否晚于最晚服务时间”由于浮点数的表达误差一个精确位于边界上的数值可能被判为晚到从而影响目标函数值。处理办法很简单在比较时加上一个小阈值例如if arrive nodes(nodeIdx,5) 1e-6 % 计为晚到 end这个 1e-6 的松弛量不会产生实质影响但能避免很多“莫名其妙被切路线”的问题。6. 从代码到论文实验和结果呈现怎么写6.1 实验设计和算例准备论文里最害怕的问题是“只有一张结果表没有实验过程”。写 VRPTW 求解实验时至少要交代以下几个方面测试算例的规模节点数量、客户数量、时间窗宽度、车辆容量限制参数设定初始温度、终止温度、衰减系数、内层迭代次数、邻域算子概率对比基线是使用人工经验路径、贪心算法、其他启发式还是公开已知最优解评价指标总行驶里程、总耗时、平均等待时间、超时次数、车辆使用数。建议在论文中放一张“算例参数表”列出所有节点的坐标范围、需求量分布、时间窗分布。这样读者才能判断你的实验是否可信。6.2 结果可视化与收敛曲线收敛曲线图是模拟退火论文的标配。横轴可以选“迭代次数”或“温度下降阶段”纵轴选当前最优目标值。不要只画一条最终最优值的曲线那样看不出算法过程。常规做法是同时画两条线一条是当前解在每个温度段的目标值另一条是历史最优目标值。两条线之间的差距能直观体现退火搜索的波动性。路径规划图建议按不同车辆用不同颜色区分并标注仓库位置。如果时间窗比较紧还可以在节点旁边标注“最早开始时间/最晚开始时间”帮助审稿人理解算法为什么给某些客户安排了特殊顺序。6.3 对比实验怎么更有说服力单纯的模拟退火结果很难让审稿人眼前一亮。更稳妥的策略是设置三组对比第一组是纯贪心算法结果第二组是模拟退火使用随机初始解的运行结果第三组是模拟退火使用贪心初始解加调优参数的运行结果。通过三组数值对比可以证明算法的改进不只是来源于元启发式自身更来源于初始解和惩罚参数设计的组合优势。如果时间允许可以再加一组“不同惩罚系数下的结果对比”说明过大惩罚会让解偏保守过小惩罚会让时间窗失去意义以此论证你选择的惩罚系数有实验依据。不要只列最佳结果也要给出平均数、最差值、标准差这会显著提升论文的实验可信度。6.4 代码附录的组织很多论文要求附代码但很多人直接贴整个脚本阅读体验很差。更理想的做法是拆分成几个有注释的函数块数据导入模块、初始解生成模块、邻域搜索模块、解码评估模块、主循环模块、可视化模块。每个模块前用一段文字说明输入输出含义并把参数设置独立成一个小节方便读者调整。Matlab 代码里尽量加中文注释并写明运行环境版本。如果代码中使用了自定义矩阵或结构体字段最好在附录开头写一个字段说明表能省掉读者大量对字段名发愁的时间。一些实操层面的经验心得我自己的体会是模拟退火求解 VRPTW 的难点不在算法本身而在“时间窗检查”和“路径解码”这两个容易犯迷糊的地方。很多开源代码喜欢把所有逻辑都堆在主循环里看起来很短但改一个参数就要全局排查远不如我上面这种“解码函数 邻域算子 主循环”的结构清晰。最后再分享一个小技巧求解之前先把仓库坐标和数据规模在图上画一遍直观感受一下客户分布。如果客户太分散时间窗又窄那无论算法调得多好也不可能达到理想中的“几辆车全走完”。先通过画图确认问题本身是不是“有资格难”再决定要不要用更复杂的算法。别一上来就写模拟退火主循环跑完之后才发现数据存在严重单位或坐标错误那样一天时间就白费了。
