简介PDPTW带时间窗的取送货问题是城市物流配送中典型的带约束组合优化难题取货与送货需在时间窗、车辆容量及先后次序下协同规划。这套Matlab工程包主要面向学习车辆路径优化、元启发式算法的本专科生及算法研发人员给出了模拟退火与遗传算法相结合的混合求解实现。包内共40个文件以26个m脚本为骨架覆盖初始化、距离计算、请求配对、交叉变异、成本评估、路径绘制等关键环节另有6个xlsx数据表、6张原解结果图与2个mat基础数据文件并附有坐标数据、请求配对表、基础实验数据等辅助内容压缩包整体仅169KB无需大型依赖即可本地运行。目前已有177人学习/下载。使用者可直接运行主程序查看不同迭代阶段的目标值与车辆路径变化借助原解图片与基础实验数据还能对比p100至500等参数设定下的优化效果便于理解PDPTW编码设计、遗传算子细节以及SA-GA混合策略的调参与改进思路。1. 为什么说PDPTW的难点在“先取后送”与时间窗叠加PDPTW 的全称是 Pickup and Delivery Problem with Time Windows也就是带时间窗的取送货问题。它和普通 VRP 最大的区别是请求成对出现每个请求由取货点pickup和送货点delivery组成同一个请求的两点必须由同一辆车服务并且必须先取后送。仅看“最短路径”会得到一条在配对关系、装卸顺序和时间窗上全部失真的路线。我拆这个PDPTW.rar时重点看的是checkParingPrecedence.m、exchange.m、crossover.m、myEvalue.m这几个文件它们分别承担优先级校验、邻域交换、遗传交叉和成本评价。整体是一个“遗传算法做广度搜索、模拟退火做局部改进”的混合结构适合想在 Matlab 中直接复现 PDPTW 求解器、又不想从零写底层约束的从业者。下面我会按建模、混合算法、约束校验、参数调优到验证的顺序把它们串起来。2. PDPTW建模先立约束再谈Matlab数据表2.1 节点、变量与三组核心约束模型层面理解到位代码才不会越写越乱。首先定义节点集合仓库节点占一个 IDP是取货点集合D是送货点集合总节点数为2n1其中n是请求数。一组典型决策变量是x_{ij}^k车辆k是否从节点i行驶到节点jt_i车辆到达节点i的时刻Q_k与q_i车辆容量与节点需求取货为正、送货为负。三组约束里最容易写漏的是成对约束。每个请求的两个点不仅要被访问还必须分配给同一辆车。用公式表达就是对于请求r如果车辆k服务了它的取货点p_r则它也必须服务送货点d_r。第二条是关键的顺序约束t_{d_r}必须大于等于t_{p_r}加上服务时间和两点间的行驶时间。第三条是容量约束车辆在任何时刻的载货量必须在0到Q_k之间这也是把取货点设成正数、送货点设成负数的原因。在实际项目里我一般不会一次性把所有约束写进非线性规划而是先保证染色体表示天然满足一部分约束其余约束放进适应度函数做惩罚。这样求解速度比直接调用 matlab优化工具箱的整数规划求解器快得多尤其是请求数超过 30 时MILP 模型变量数量膨胀内存和分支定界都会失控。2.2 coodinate.xlsx、requestparing.xlsx 与 Dmand.xlsx 怎么组织这个压缩包里数据是拆开存放的coodinate.xlsx放所有节点的二维坐标requestparing.xlsx记录请求与取货/送货点的配对关系Dmand.xlsx放节点需求。三者通过节点 ID 关联常用字段如下文件字段示例说明coodinate.xlsxId, X, Y节点 ID 与横纵坐标requestparing.xlsxRequestId, PickupId, DeliveryId请求 ID、取货点 ID、送货点 IDDmand.xlsxNodeId, Demand节点需求取货为正、送货为负读取时有一个易错点readtable会把 Excel 第一行自动视为表头如果原表没有表头导入后字段名会变成Var1、Var2。所以我在读完之后习惯用coord.Properties.VariableNames确认一次。下面这段代码把三张表读进来并生成距离矩阵coord readtable(coodinate.xlsx); paring readtable(requestparing.xlsx); demand readtable(Dmand.xlsx); nReq height(paring); pickupNode paring.PickupId; % 取货点节点ID deliveryNode paring.DeliveryId; % 送货点节点ID nNode height(coord); dist zeros(nNode, nNode); for i 1:nNode xi coord.X(i); yi coord.Y(i); for j 1:nNode dx xi - coord.X(j); dy yi - coord.Y(j); dist(i,j) sqrt(dx*dx dy*dy); end end注意nNode里包含仓库节点这里我默认仓库 ID 为 1客户节点从 2 开始编排这样节点 ID 直接对应 Matlab 行号不会出现下标越界。distance.m在原始项目里就是做这件事只是它通常只返回二维距离矩阵不给时间窗矩阵。时间窗得另建一个nNode x 2的数组第一列是最早开始服务时间第二列是最晚开始服务时间。如果网点允许车辆提前到达并等待只要把到达时间更新为max(arrival, timeWin(node,1))即可。2.3 inatialize.m用“先取后送骨架”生成初始路线初始化这里最容易犯的错是直接randperm生成一段所有节点均参与排列的序列。这种方法产生的解大概率出现“送货点排在取货点前”的非法个体后面再怎么优化都是补救。我在重写初始化逻辑时采用先把请求编号随机打乱再把每个请求的取货点和送货点紧挨着展开的办法function route inatialize(nReq, pickupNode, deliveryNode, depot) % depot 是仓库节点ID建议取 1避免 Matlab 下标从0开始的问题 pOrder randperm(nReq); % 随机打乱请求编号 route []; for r pOrder route [route, pickupNode(r), deliveryNode(r)]; % 取货点后紧跟送货点 end route [depot, route, depot]; % 首尾都是仓库 end为什么要这样做因为“取货点紧跟送货点”是满足先取后送的最强充分条件任何基于这条序列的交换和交叉天然不会产生送前取后的非法顺序。缺点是它会把每个请求的两个点绑得太紧时间窗上可能不够松弛。不过没关系后续的exchange.m变异操作会逐渐把送货点向后挪从而在不破坏优先顺序的前提下放宽时间窗。真正的耗时不在初始化而在后面每一代里的合法性检查。这一章的内容到这里已经把建模、数据表和初始化解全部打通。接下来进入算法核心。3. 混合策略遗传算法做广度模拟退火做局部收尾3.1 为什么要把遗传和退火放在一起只用遗传算法求解 PDPTW常见问题是种群在 20 代以后迅速同质化交叉算子产生的子代与父代几乎没有结构差异然后陷入局部最优。只用模拟退火又受限于单点搜索邻域结构稍微设计得复杂一些比如引入 2-opt 加 exchange 的组合变异就很容易在接近最优时反复接受劣解、收敛缓慢。PDPTW.rar这套方案选择了折中遗传算法负责在“请求序列”这个离散空间里做全局勘探模拟退火则在外循环的每一代里用 Metropolis 准则决定是否接受交叉和变异后的新个体。这样做还有一个额外好处模拟退火本身不需要额外实现“种群”的概念只要把种群中的每个个体都看成一条独立的退火链温度统一按代数下降就可以保持种群多样性。相当于用温度给整个种群加了一个全局选择压力避免过早收敛。3.2 适应度函数必须包含时间窗惩罚myEvalue.m与costtest.m承担了适应度计算。PDPTW 的成本不只是距离时间窗惩罚和载重惩罚都要折算进去。我一般把目标函数写成function cost myEvalue(route, timeWin, serviceTime, dist, demand, Q, depot) totalDist 0; totalPenalty 0; t 0; % 当前时间 load 0; % 当前载重 for i 1:length(route)-1 from route(i); to route(i1); totalDist totalDist dist(from, to); t t dist(from, to); if to ~ depot if t timeWin(to, 1) t timeWin(to, 1); % 早到需要等待 end if t timeWin(to, 2) totalPenalty totalPenalty (t - timeWin(to, 2)); end t t serviceTime(to); load load demand(to); if load Q || load 0 totalPenalty totalPenalty 1e4 * (abs(load) 1); end end end cost totalDist 1e3 * totalPenalty; end这里的核心参数是惩罚系数1e3。如果调得太小算法会倾向于生成长距离绕行来避免时间窗惩罚路线会变得怪异调得太大可行解与不可行解差距悬殊搜索会直接跳过所有轻微超时的解导致无法利用“接近可行”的中间状态。常见做法是先固定为1e3观察成本曲线变化最后再微调。Dmand.xlsx中的Demand是取货正、送货负因此load load demand(to)这一句天然计算了装卸后的即时载重。3.3 在“请求编号”上做交叉而不是在节点序列上这是整个遗传设计里最值得注意的细节。直接在节点序列上做两点交叉会频繁产生同一请求只保留取货点或只保留送货点的个体。即使对缺失点做修补修复逻辑也很容易破坏另一条染色体的优秀片段。crossover.m里我采用的做法是染色体只保存请求编号顺序交叉后再展开成节点序列。function child crossoverOX(parent1, parent2, nReq) % 顺序交叉parent1 和 parent2 都是请求编号排列 cut randi([1, nReq-1]); % 保留父本1在切点前的一段 seg parent1(1:cut); % 从父本2中剔除已经选中的请求保持父本2的相对顺序 rest parent2(~ismember(parent2, seg)); child [seg, rest]; end这种交叉叫做顺序交叉Order Crossover。它保证每个请求只出现一次且请求级顺序不会出现重复或缺失。得到child请求序列后再用inatialize里的展开逻辑把它转成一个满足先取后送的节点序列。这样做之后交叉产生的子代永远满足“先取后送”的强条件代码里甚至可以省掉一部分优先级校验需要校验的场景只剩变异算子。variation.m里的变异操作也必须放在请求编号层。常见的变异是随机交换两个请求的位置或者把一段请求子序列倒置。倒置本质上是 2-opt 在请求层的一种体现它可以翻转整段路线的服务顺序在时间窗调整上比单点交换更有效。我一般让交换变异概率取 0.1倒置变异概率取 0.05整体变异概率控制在 0.15 到 0.3 之间。3.4 exchange 与 group在街道层面做邻域搜索请求级交叉负责“粗调”exchange.m和group.m负责“细调”。exchange在展开后的节点序列上进行常见动作是交换两个取货点的位置或把一个送货点插入到另一个位置。这种操作一旦做了就必须重新检查顺序和时间窗所以PDPTW.rar里会反复调用check.m。我在实际改造时把exchange的邻域限制为“只交换同一时间段内开始服务的节点”因为时间窗相差很远的节点交换没有意义只会产生大量无效扰动。用一个简单判断如果两个节点的时间窗中心点相差超过 20% 的总时间窗跨度就跳过这次交换。这个策略能显著降低适应度函数的调用次数在 100 个请求规模下运行速度能快三分之一左右。好策略部分先到这。接下来看约束校验和参数标定。4. 约束校验、车辆分配与参数调优把“合法解”落到实处4.1 checkParingPrecedence.m为什么要单独抽一个优先级校验在混合算法里交叉和变异已经尽量保证了先取后送但 exchange 操作仍然可能制造出“送货点提前于取货点”的非法片段。因此每次邻域搜索后都应该执行一次独立优先级校验。checkParingPrecedence.m的核心逻辑是遍历路线时维护一个“是否已取货”的布尔数组function ok checkParingPrecedence(route, pickupNode, deliveryNode, nReq) ok true; hasPickup false(1, nReq); % route 第一个和最后一个节点都是仓库跳过 for node route(2:end-1) pi find(pickupNode node, 1); if ~isempty(pi) hasPickup(pi) true; continue; end di find(deliveryNode node, 1); if ~isempty(di) ~hasPickup(di) ok false; return; end end end这段代码的巧妙之处在于利用find(..., 1)直接把节点 ID 映射回请求 ID不需要额外查表。如果路线中先出现某个请求的送货点hasPickup还是false直接返回非法。实际应用时我在交叉算子后不调用这个函数因为请求级交叉天然合法但在exchange.m和variation.m的节点级邻域搜索后必须调用。出现非法个体时最简单的处理是拒绝接受而不是修复因为修复一个配送节点很可能产生新的时间窗超限。4.2 vehiclelocation 与 group把一条长染色体拆给多辆车初始化和交叉产生的路线只是“一辆车”的超长访问序列。真实 PDPTW 里车辆有容量限制所以需要把节点序列切分给多辆车。vehiclelocation.m负责确定每个节点由哪辆车访问group.m负责把节点聚成车辆组。常见做法有三种我整理在这张表里切分策略实现思路适用场景顺序切分按路线顺序累加载重超过容量就切开请求规模小时间窗宽松时间窗切分按最早开始时间排序后分配时间窗密集需要减少等待按请求整体切分永远把一个请求的取、送两点放进同一辆车容量与配对同时受限PDPTW.rar里的group.m更接近第三种。它把请求级染色体按顺序往车辆里填新请求的载重加上当前车辆已有载重不超过容量时放入当前车否则换下一辆车。因为请求不可拆分所以不会出现运送同一请求的取货点和送货点使用不同车辆的情况。vehiclelocation.m再按每辆车服务的节点集合生成实际的行车路线。这里有一个隐藏毛病顺序切分在时间窗上往往是次优的。比如第 3 个请求虽然容量上可以装进第一辆车但它要求上午 10 点送达而第一辆车已经排到中午这时把它分给后面空闲车辆反而更好。为此我在group.m里增加了一个时间窗预判如果加入该请求后预计到达时间超出最晚时间窗超过 20%就强制换下一辆车。这样牺牲了一点载重利用率但能减少后续惩罚项。4.3 参数标定温度、冷却系数、种群规模怎么设混合算法的参数比单一遗传算法多调参顺序也相对固定。下面是我在请求数为 50 到 200 时常用的默认范围参数常见范围建议值说明种群规模 populationSize50 ~ 200请求数接近 200 时取 120 即可初始温度 T0初始最差解与最优解成本差的 5~20 倍越大越容易接受差解冷却系数 alpha0.90 ~ 0.990.95 起步跑不动就改成 0.98变异概率0.05 ~ 0.3时间窗紧时增大扰动外循环代数 maxIter200 ~ 2000每增加 50 个请求代数加 200外循环的骨架通常是这样T T0; for iter 1:maxIter T alpha * T; % 指数降温 for p 1:populationSize newRoute mutateOrCrossover(pop(p), pop(randi(populationSize))); newCost myEvalue(newRoute, timeWin, serviceTime, dist, demand, Q, depot); if newCost popCost(p) pop(p) newRoute; popCost(p) newCost; elseif rand exp(-(newCost - popCost(p)) / T) pop(p) newRoute; popCost(p) newCost; end end end这一段就是遗传与模拟退火的接口变异或交叉产生子代Metropolis 准则决定是否替换。和纯遗传算法不同这里的劣解替换概率随着T的下降而减小所以前期可以充分探索后期收敛到稳定路线。PDP.m常被作为顶层调用Run_VRP.m则承担算法主循环先读数据、初始化种群、计算T0然后进入这个外循环。注意T0不应该拍脑袋定。我一般先随机生成一百条初始路线统计它们的成本方差把T0设为方差的 10 倍左右alpha则根据期望迭代次数反推如果定了maxIter1000希望末尾温度是初始温度的1/100那么alpha (0.01)^(1/1000) 0.9954。这套推导比反复试错快得多。约束和参数都对齐后剩下的问题就是“怎么确认解是对的”。下一章讲验证技巧。5. 验证方法把解绘出来并量化低成本5.1 用 DrawRoute 与 circleDrawRoute 检查路线合法性DrawRoute.m和circleDrawRoute.m在项目里承担可视化验证。前者把节点按顺序连线后者在每个节点周围画圆圆的半径或颜色深浅体现时间窗松紧。检查画面时我重点看三件事同一请求的取货点与送货点是否连在一条颜色线内每一条线路上取货点是否都在送货点之前是否存在明显交叉的绕行路线。figure; DrawRoute(bestRoute, coord, b-); circleDrawRoute(bestRoute, coord, timeWin, 5); title(PDPTW best route);circleDrawRoute的最后一个参数是圆半径缩放系数用来让时间窗可视化更明显。如果某节点圆半径特别大说明该节点时间窗宽如果特别小说明时间窗紧车辆在这里的到站时间必须非常准。这条信息通常能直接看出哪些请求约束了整条路线。5.2 用成本分解确认“假优化”只看总成本下降曲线很容易被误导。一种典型现象是总成本在前 50 代快速下降之后就停下来但路线其实违反了载重约束只是惩罚项占比太小没体现出来。所以我建议在每次最优解更新后打印三个分量[distPart, penaltyPart, loadPart] myRecord(bestRoute, dist, timeWin, demand, Q, depot); fprintf(距离 %.2f | 时间窗惩罚 %.2f | 载重惩罚 %.2f\n, distPart, penaltyPart, loadPart);myRecord.m在这里的作用是把myEvalue的总成本还原成分项。如果loadPart不为零说明当前解连基本载重约束都不满足应该立即调高惩罚系数或加强check.m的拒绝机制如果timeWin惩罚收敛到 0说明算法已经找到可行解剩下的优化才是真正压缩行驶距离。5.3 一个提速技巧超时阈值前置拦截最后说一个对大算例很管用的技巧不要等完整计算完成本再接受或拒绝新解。在exchange或变异产生新路线后先快速检查累计时间窗超时量penaltyPart如果超过预设阈值maxPenalty直接continue不再计算剩余节点。这个方法能让无效解消耗的时间大幅减少。if penaltyPart maxPenalty continue; % 跳过严重超时解不做完整适应度计算 endmaxPenalty的初始值可以设为当前最优解时间窗惩罚的 2 倍每 20 代更新一次。把maxPenalty从1e3开始试观察收敛曲线通常能在前 100 代区分出是惩罚系数不够还是时间窗数据本身排得过密。本文还有配套的精品资源点击获取