混合狼群退火算法求解动态TWVRP的Matlab实战

做配送路径优化的同行应该都有类似体验凌晨刚排好的线路上午十点客户一个电话进来说下午三点前必须送到整条计划全乱。这就是TWVRP带时间窗车辆路径问题叠加动态需求后的真实日常。我前两年接了个城市生鲜配送的小项目越做越觉得静态排线只是基本功真正难的是怎么在客户不断插单的情况下快速给出一版新路线。当时试过纯狼群算法也试过裸模拟退火结果都不理想——单跑狼群容易早熟单跑退火又慢得没法落地。后来狠下心把两个算法揉进一个Matlab框架里做混合求解才把动态TWVRP给啃下来。这篇文章不聊虚的直接讲清楚问题建模、混合算法怎么设计、Matlab代码核心模块怎么拆、参数怎么调以及我踩过的几个坑。1. 先想明白TWVRP在求什么难点卡在哪几个环节1.1 从业务描述到数学模型目标函数和约束到底长什么样TWVRP的标准描述是有一个配送中心、若干客户、若干辆同型车每个客户有需求量和指定服务时间窗比如门店必须在8:30到10:00之间收货车辆从配送中心出发、完成服务后返回。问怎么样安排车辆路径使总成本最小同时尽量满足时间窗要求。落到数学上这其实是一个混合整数规划问题。定义决策变量x_ijk车辆k从客户i直接开到客户j则为1否则为0y_ik车辆k服务客户i则为1否则为0t_i客户i的开始服务时间。目标函数通常写成min Z ΣΣΣ d_ij · x_ijk λ · Σ penalty_i其中d_ij是i到j的行驶距离penalty_i是客户i的时间窗违反量λ是惩罚系数。时间窗违反量可以硬性定义也可以软性定义。例如对客户i预约窗为[e_i, l_i]车辆到达时间为arr_i则penalty_i max(0, e_i - arr_i) max(0, arr_i - l_i)这里有一个细节我在做项目时反复确认过早到和晚到都要付出代价。早到意味着车辆在客户门口干等、占用司机工时晚到更严重客户可能拒收。所以不能只罚迟到早到的等待时间也必须算进成本里。约束条件方面核心有这么几条每个客户必须且只能被一辆车服务每辆车从配送中心出发并最终回到配送中心车辆载货总量不能超过容量Q时间窗约束要满足开始服务时间落在[e_i, l_i]内硬窗或尽量满足软窗。此外还有服务时间的递推关系t_j ≥ t_i s_i t_ij其中s_i是客户i的服务时长t_ij是车辆从i到j的行驶时间。为什么要强调这个递推关系因为路径问题的本质不只是排序而是排序 时间传播。一个客户迟到10分钟后面所有客户的服务时间全部顺延这条路径的可行性可能立刻崩溃。这就是TWVRP比普通VRP难一个档次的核心原因。1.2 惩罚函数还是严格约束搜索空间设计决定算法死活很多刚接触TWVRP的人第一反应是时间窗是硬约束违反就判不可行这样不是更严谨吗理论上是但实际跑起来会发现一个大问题——当把时间窗硬约束直接塞进适应度函数时可行解的占比会变得非常低。尤其在客户时间窗狭窄、车辆数量紧张的算例里随机生成的路径几乎是十不存一。一旦可行解太少元启发式算法的搜索就退化成在可行解边缘反复试探狼群算法的探狼游走也变成无效的随机抖动。所以我在这类问题上一律采用惩罚函数法先把所有约束都改成目标函数里的惩罚项让算法在完整解空间里跑然后通过调节惩罚系数λ逐步把解拉向可行区域。这样做搜索更顺畅而且能处理时间窗只是软约束的实际业务场景。惩罚函数法的代价是把约束是否满足的问题变成了惩罚系数怎么标定的问题。λ太小惩罚形同虚设λ太大目标函数被时间窗主导车辆会为了避开早晚高峰绕远路。这个坑我在第5章单独展开讲。1.3 标题里的动态规划到底指什么和经典DP不是一回事这里必须澄清一个容易混淆的点。标题写的是车辆路径动态规划问题但这里的动态指的是需求随事件动态到达也就是文献里常说的DVRPTWDynamic Vehicle Routing Problem with Time Windows不是指动态规划Dynamic Programming算法。我第一次看这个标题也差点想歪以为要用DP解TWVRP那样问题规模稍微一大就爆炸了。动态TWVRP和静态TWVRP的区别在于静态版本假设所有客户在排线前已知一次性求解动态版本则假设初始只有部分客户车辆发车后新的客户订单会随机到达需要实时决定插入哪辆车、插在哪条路径的哪个位置。比如生鲜配送场景早上只派一批固定门店的订单中途会不断有补货订单进来这时候你要根据车辆当前实际位置、剩余载货量、已服务客户列表和剩余时间窗快速生成一版新路线。实际项目里的主流做法是滚动时域重调度rolling horizon reoptimization把时间切成片段车辆行驶过程中持续监控事件事件到达时就把当前未服务客户 新到达客户 车辆当前位置构造成一个新的子问题调用元启发式算法在有限时间预算内求解。这也是我采用狼群算法加模拟退火的根本原因——动态环境下真正稀缺的不是最优性而是在几十秒内给出足够好的解。2. 为什么是狼群退火两种算法各自的脾气与互补逻辑2.1 狼群算法不神秘头狼、探狼、猛狼的分工关系狼群算法Wolf Pack Algorithm, WPA是典型的群智能算法核心思路模拟狼群捕猎时的分工机制。整个过程有三个关键角色头狼当前种群中适应度最好的个体相当于全局最优解候选引导整个狼群的搜索方向探狼分布在解空间各处做独立游走负责探索新区域如果发现比头狼更好的位置就取而代之猛狼收到头狼召唤后朝头狼方向快速奔袭完成局部包围配合头狼做精细搜索。算法每一代大致分这几步初始化种群选出头狼探狼在各自邻域内游走若干步挑出邻域最优猛狼朝头狼位置靠近如果猛狼找到更优解则更新头狼种群更新后淘汰适应度最差的个体用新随机个体补充进来保持种群活性。和遗传算法GA相比WPA没有交叉算子主要靠邻域操作和位置更新来完成搜索。这意味着邻域操作的设计质量对算法效果影响极大。我在TWVRP上用的邻域操作是三类2-opt反转路径片段、insert把一个客户插入另一位置、swap交换两个客户位置。探狼游走时随机选一种算子执行执行完保留更优解。2.2 模拟退火不是来降温的而是来给算法留后路的模拟退火Simulated Annealing, SA的核心机制是Metropolis准则当前解S产生邻域解S如果S更优就接受如果更差则以概率exp(-ΔE/T)接受。温度T越高接受差解的概率越大随着T逐渐降低算法从广撒网过渡到精确打磨。单看这个机制SA最大的问题是爬坡速度慢温度从几百降到零点几每一层温度还要做多次内循环算一次TWVRP的适应度又涉及解码和可行性校验整体耗时非常难看。但它有一个WPA不具备的优点在搜索后期不排斥差解天然有跳出局部最优的能力。我把它定位成后路而不是主搜索器。WPA的问题在于头狼一旦锁定猛狼与探狼都围绕着它转复杂度不够时很容易早熟而SA恰好能在WPA每个迭代的末段对当前最优解做一轮局部扰动以一定概率接受略有变差的解帮助跳出土坑。2.3 两种耦合方式我选的是WPA迭代内嵌SA局部精化混合算法有几种常见搭法我实际试过三种串行级联WPA先跑完再用SA精化最终解。结构简单但动态场景下每次事件到达都要跑两段完整搜索时间预算很容易超内嵌局部算子在WPA每轮迭代末对头狼调用一轮SA局部搜索搜索结果参与下一轮竞争完全融合SA每次产生的解直接参与头狼竞争相当于给WPA持续注入多样性。对比一下三种方式的实测表现耦合方式计算开销解质量动态响应能力适用场景串行级联高较高差小规模静态问题内嵌局部算子中高好动态重调度完全融合中高高中需要强多样性的复杂算例我最终采用的是第二种WPA负责全局搜索和种群进化SA作为局部精化算子只在头狼和少数优秀个体上执行每次SA内循环次数控制在30~80次。这样既不拖慢WPA主体迭代又能把头部解打磨得更精细。动态事件到达时这个结构还能配合时间预算灵活截断——如果只剩5秒就把SA循环次数减半WPA迭代次数减半优先保证给出一个可用解。3. 从数据到路径Matlab实现的核心模块怎么搭3.1 数据组织客户、时间窗、车辆参数的统一结构Matlab写TWVRP最忌讳把数据散落在几十个全局变量里。我习惯用一个struct统一管理模型参数代码清晰又不容易出错。下面是一个典型的数据定义方式% model 结构体 % locations: [x坐标, y坐标, 需求量, 最早到达时间e, 最晚到达时间l, 服务时长s] model.locations [ 50 50 0 0 230 0; % 配送中心 24 38 11 30 90 90; % 客户1 18 40 13 50 120 80; % 客户2 22 45 9 80 150 100; % 客户3 35 30 16 15 60 75; % 客户4 28 52 10 90 180 85; % 客户5 ]; model.vehicleNum 3; model.capacity 80; model.speed 40; % 单位km/h用于距离转时间 model.depot 1; % 配送中心编号为第1行这里要特别提醒一个单位问题。客户给的最晚送达时间通常精确到分钟但你算出来的t_ij如果是小时时间窗比较就会出错。我踩过一次距离算出来12公里车速40km/h行驶时间0.3小时直接拿0.3和客户的时间窗单位分钟做比较结果全乱了。所以我在做数据预处理时统一把时间单位转换成分钟距离转时间的公式是t_ij d_ij / speed * 60千万别省这一步。3.2 编码与解码一条染色体怎么变成一组车辆路径车辆路径问题最常用的编码方式是客户排列编码也就是用一条客户序列表示解。比如有8个客户染色体的样子就是[3 1 5 2 7 4 6 8]表示访问客户的顺序。解码时按照顺序依次把客户分配给车辆直到车辆超载或违反硬时间窗再切换到下一辆车。解码函数的核心逻辑如下function routes decodeTWVRP(chrom, model) distMat model.distMat; Q model.capacity; demand model.locations(:,3); e model.locations(:,4); l model.locations(:,5); service model.locations(:,6); routes {}; curRoute []; curLoad 0; curTime 0; % 从配送中心发车时间 curPos model.depot; for i 1:numel(chrom) c chrom(i); travel distMat(curPos, c) / model.speed * 60; arr curTime travel; startTime max(arr, e(c)); % 早到则等待 if curLoad demand(c) Q startTime l(c) curRoute [curRoute, c]; curLoad curLoad demand(c); curTime startTime service(c); curPos c; else routes{end1} curRoute; % 结算当前车辆 curRoute [c]; curLoad demand(c); curTime distMat(model.depot, c) / model.speed * 60; curPos c; end end routes{end1} curRoute; end这只是一种按顺序贪心分割的简化解码实际项目里还要考虑回程时间是否超过配送中心的关门时间、车辆数量是否超限等细节。但这个结构已经能对付多数TWVRP算例而且解码速度快适合在算法迭代里反复调用。3.3 狼群三大行为的Matlab实现要点WPA在Matlab里的实现关键是把游走奔袭围攻翻译成针对客户排列的邻域操作。探狼游走这一段我这样写function newSol exploreWolf(sol, model, stepNum) % 探狼游走在邻域内做stepNum次搜索保留最优 best sol; for k 1:stepNum candidate neighOperator(best, model); % 随机选2-opt/insert/swap if fitness(candidate, model) fitness(best, model) best candidate; end end newSol best; endneighOperator是我自己的一个包装函数内容是从2-opt、insert、swap三个算子中随机挑一个执行。这里有个经验游走步数stepNum不是越大越好。步数太大探狼会变成局部爬山器过早收敛到头狼附近的区域步数太小又是随机抖动。我一般设置在5~15之间。猛狼奔袭的实现则是让个别个体向头狼靠近具体做法是取头狼路径的一段连续客户序列替换掉当前个体中对应的一段其他位置保持原样。这本质上是一种定向的路径片段传递比完全随机扰动更有目的性。种群更新阶段要注意保持种群规模恒定。每轮淘汰最差的若干个个体同时随机生成新个体补位。这里的随机个体不是纯随机排列而是基于当前头狼做一次大扰动得到——保留头狼的主要结构随机打乱其中一两个片段既维持多样性又不至于丢过头狼的优质信息。3.4 模拟退火做局部精化Metropolis准则写法与参数配套SA在整个框架里的角色是局部精化器代码如下function solNew saLocalSearch(sol, model, T0, T_end, alpha, innerLoop) T T0; cur sol; while T T_end for j 1:innerLoop cand neighOperator(cur, model); deltaE fitness(cand, model) - fitness(cur, model); if deltaE 0 || rand exp(-deltaE / T) cur cand; end end T T * alpha; end solNew cur; end这里面每个参数都很关键。T0是初始温度alpha是降温系数innerLoop是每一温度层的内循环次数。我调参时先把T_end固定为1alpha固定为0.95只调T0和innerLoop。初温对WPASA的耦合效果影响非常大原因我放到第5章讲。这里先给一组我实验中比较稳的初始值T0500T_end1alpha0.95innerLoop20。对25客户的算例这样一次SA局部搜索大约要做上百次邻域评估耗时在0.5秒以内可接受。3.5 主循环静态排线、事件到达、局部重优化三步走完整的动态TWVRP求解流程在主循环里串起来逻辑分三段静态阶段初始客户集合出发跑WPASA得到基线方案把各车辆的初始路径和时间计划下发事件推进模拟时钟推进车辆执行计划。每到一个事件时刻判断是否有新客户到达。有新订单时先把已服务客户从问题中剔除把各车辆当前位置和剩余容量收集起来参考当前时间修正剩余时间窗重优化阶段把未服务客户 新客户构造成新的子问题调用WPASA求解。这里要设置最大重优化时间比如5秒若超出时限则中止搜索用当前最好解作为新方案若连可行解都没有用最近车辆顺路插入的方式紧急兜底。第二段收集车辆当前位置最容易写漏。车辆不能当作还在配送中心它可能正在去往某个客户的路上。我的处理方式是在模型里增加一个虚拟出发点把每辆车的当前位置当作该车的0号节点该节点的时间窗取当前时间距离矩阵里加入该虚拟节点到各客户的距离。这样重优化时每辆车的起点就不一样了解出来的路径才是真实可执行的。4. 实验跑了三轮参数怎么调、动态场景怎么处理、效果如何4.1 静态基准算例对比混合算法在解质量上的收益为了验证WPASA不是花架子我在修改过的Solomon C101算例基础上跑了一组静态对比实验客户数分别取15、25、50每个规模跑20次取平均。25客户的结果如下算法平均路径距离平均违反时间窗数平均耗时WPA单独971.40.2111.2sSA单独968.70.059.6sWPASA953.5014.3s单独SA的表现其实不差因为它本质上是强局部搜索对中小规模算例有天然优势但它的缺点是稳定性差跑10次有几次会陷进很差的局部极小。WPA单独跑则明显偏早熟25客户时就已经出现违反时间窗的情况。混合后的优势在于WPA负责把搜索范围打开SA在每轮帮头狼做精修最终解质量稳定且很少违反时间窗。到50客户规模时差距更明显WPA单独平均路径距离1842.6违反时间窗1.8次WPASA是1806.3违反0.3次耗时约38秒。这说明问题规模越大耦合算法的优势越突出。4.2 参数敏感性哪些参数决定生死参数调优是这类项目的重头戏。我把MATLAB主程序里的核心参数单独做了一个参数扫描结论如下参数推荐范围过量后果不足后果狼群规模pop20~50耗时平方级增长探索不足早熟探狼游走步数5~15探狼退化成局部搜索接近随机抖动温度初值T0300~800前期浪费大量迭代局部精化失效降温系数alpha0.90~0.98收敛过慢温度骤降解差SA内循环数15~40单次调用太耗时局部搜索不充分调参顺序比我预想的更重要。我最开始是同时调五个参数结果跑了两个通宵也说不清是谁在起作用。后来改成先固定SA侧参数只扫描狼群规模和游走步数定下来之后再固定WPA侧参数扫T0和alpha。这样每次只有一个变量在动结论才干净。另外WPA还有一个隐藏参数是每轮淘汰比例我固定在20%。淘汰比例太高种群多样性崩塌太低劣质个体一直占用种群空间。20%是实验里比较中庸的选择。4.3 动态事件处理重优化频率和局部兜底策略动态场景实验我模拟了一个25客户、5次随机新增订单的配送周期比较两种重调度策略策略A每次事件到达都对整个未服务客户集合做全量重优化策略B只对受影响区域内的未服务客户重优化超出限时就用顺路插入兜底。策略单次平均重优化耗时总行驶距离司机路线变动总次数全量重优化8.5s972.314局部重优化兜底1.6s988.17策略A路径距离略优但有一个致命问题频繁全量重优化让司机收到的路线每次都在大改现场执行根本跟不上。策略B虽然在总距离上差约1.6%但路线变动次数少了一半实际落地时司机的接受度反而高很多。这里我学到一个做动态优化的通用原则动态问题里解的稳定性和解的最优性同等重要很多时候稳定性比最优性更值钱。5. 复盘几个坑时间窗惩罚、初温设置、动态解稳定性5.1 惩罚系数λ怎么标定才不会出为了守时绕远路第1章说了惩罚函数法这里讲它最坑的地方。如果λ没有做量纲换算你会看到算法为了少付时间窗惩罚金让车辆绕一个巨大的圈避开高峰期到达总行驶距离反而暴涨。这不是算法蠢是你把惩罚的单位价格定错了。我的标定方法是先跑若干次完全不加惩罚的随机解统计平均距离和平均时间窗违反量算出两者量级比。假设平均距离在1000公里量级平均违反时间窗总量在500分钟量级那么1公里距离约等于0.5分钟λ的合理起点就让两者价值相当即λ2左右。然后按2倍步长上下扫描看解的变化趋势。更稳的做法是分段惩罚时间窗违反在30分钟内罚0.530~60分钟罚1超过60分钟罚2。这种分段方式比单一λ更有韧性能避免个别严重迟到的解挤掉整体质量好的解。缺点是多了一个参数需要看业务成本来定但一旦定好后续算例基本可以直接复用。5.2 初温设置的连锁反应选错初温整个混合算法都失灵我试过几次把T0直接设成1000或更高结果狼群每轮的SA局部搜索都变成随机游走因为温度太高时Metropolis几乎接受一切差解SA精化的意义被清零。反过来T050时SA只接受小幅改良不到两轮就退化成爬山算法混合效果聊胜于无。一个比我拍脑袋可靠得多的经验做法是采样估算初始温差。在正式运行前随机生成30个候选解两两比较目标值取所有目标差值的均值作为ΔE_avg然后按公式T0 -ΔE_avg / ln(p_accept)其中p_accept取0.8~0.9表示希望SA在运行初期以80%~90%的概率接受中等幅度的差解。这个公式能让你在不同规模算例上都得到合理初始温度不用靠感觉。25客户算例算出来T0通常在300~600之间和我的经验范围基本吻合。5.3 动态重调度最大的坑解的颠簸与时间窗平移最后说一个比较隐蔽的问题解颠簸。动态事件到达后重优化生成的新方案可能在数值上更优但和上一版方案差异极大——司机刚收到新路线还没回过神来又收到一条完全不同路线。这种时域上的不稳定会让实际调度成本远超路径长度本身。我在项目里做了两件事抑制颠簸第一在目标函数里加入方案变动惩罚。新方案中每辆车的客户序列如果和旧方案的差异达到一定比例就增加一个惩罚项。这个惩罚随差异程度线性增长迫使算法在更优和更稳之间取平衡。第二设置重调度触发阈值。新订单到达时并不一定立刻重调度。先评估把新客户直接插入现有路线的代价如果插入导致的额外距离增量小于成本阈值就原地插入不跑重优化只有增量超过阈值时才触发WPASA重调度。这个策略能把重调度频率降下来路线变动次数显著减少。这套机制跑完实际项目之后我的体会是动态TWVRP真正的难点不在于把某个静态算例调出多好的解而在于当车辆已经在路上时你拿到的解既要让总成本可控又要让司机不造反。狼群算法负责打开搜索面模拟退火负责精修再加上一层稳定性约束三者缺一不可。如果读者想在这个方向上继续深入我建议下一步做两件事一是把解码逻辑改成支持多配送中心和时间窗平移二是尝试用机器学习预估新订单的空间分布提前预留车辆弹性这会让动态重调度的压力小很多。