简介面向需要求解带时间窗车辆路径问题VRPTW的研究者与学习者这套代码合集提供CW节约算法、禁忌搜索硬约束版与惩罚函数版和LNS大邻域搜索四种主流启发式算法的Matlab实现对比。针对配送中心派车、客户时间窗与需求约束等典型场景可直接运行主程序观察不同算法的路径优化效果与差异。压缩包共372个文件包含90个m脚本、280个txt说明及2个mat数据文件整体859KB。m文件为核心算法与初始化/构造程序txt多为结果记录或说明文档mat用于保存实验数据目录结构简洁便于对照学习。已有729人下载学习适合正在做课程设计、毕业设计或算法对比实验的Matlab用户。通过运行本合集可快速获得四种算法在VRPTW上的求解表现、收敛过程与路径图节省自行编码时间也能为论文或报告提供对比数据与参考实现。1. 为什么 VRPTW 合集里要同时保留四种方法在物流配送、维修调度和应急物资运输场景里带时间窗的车辆路径问题VRPTW几乎是绕不开的建模对象。它比标准 VRP 多出一层时间窗口约束每个客户有最早和最晚服务时间车辆早到要等待晚到则直接违约这使得解空间形状发生本质变化。我在实际项目中遇到的情况是业务方要的不只是一个能跑出结果的算法而是想知道不同的求解策略在约束强度不同时各自表现如何。CW 节约算法简单直观、适合做初始解和教学基线TS 禁忌搜索在硬约束和软约束下有两种截然不同的写法LNS 大邻域搜索则擅长在稀疏散点和大规模算例上做破坏与修复。这四者放在一起才能真正回答VRPTW 用什么方法解更好的问题。这篇文章按照我处理此类问题的习惯组织先给出 VRPTW 的统一数学模型再分别用 Matlab 实现 CW、TS 硬约束版、TS 惩罚函数版和 LNS 四种方法最后在标准测试算例上做横向对比。读者可以拿着代码直接跑也能理解每种方法的收敛行为和参数敏感点。2. 问题模型与四种算法的适用边界评估框架2.1 VRPTW 的数学模型与符号定义VRPTW 的标准描述是一个配送中心编号 0、N 个客户编号 1 到 N、K 辆容量为 Q 的相同车辆。每个客户 i 有需求量 q_i、服务时间 s_i 和时间窗 [e_i, l_i]。车辆从配送中心出发完成服务后返回配送中心目标是使总行驶距离最小同时满足容量约束和时间窗约束。用 x_ijk 表示车辆 k 是否从客户 i 行驶到客户 j那么目标函数可以写为min Σ Σ Σ d_ij · x_ijk对 k、i、j 求和其中 d_ij 是客户 i 到客户 j 的距离。约束条件包括四条每个客户只能被一辆车服务一次每辆车载重不超过 Q车辆到达客户 i 的时间必须在时间窗内若早到则等待若晚到则不可行所有车辆的路径都从配送中心出发并最终回到配送中心。在硬约束版本中时间窗违约是完全不允许的解一旦出现违约就被判为不可行在软约束版本中违约可以被接受但其程度会以惩罚项形式加入目标函数。这个数学框架是后续所有算法的公共基础。CW 节约算法直接利用距离矩阵计算节约值并构造路径集合TS 和 LNS 则在这个框架上定义邻域操作和目标函数评估。我在实现前会先确定三个关键参数时间窗宽度、车辆容量和客户点规模因为它们决定了哪种求解策略更适合当前场景。2.2 四种方法的选型逻辑与适用场景对照在实际工程中选择哪种算法取决于问题规模、约束严格程度和可用计算时间。CW 节约算法是一种构造型算法一次贪心合并得到解速度快但没有迭代改进机制适合作为初始解生成器或小规模问题的直接解法。TS 禁忌搜索采用局部搜索框架依靠邻域移动逐步改进解硬约束版通过严格筛选合法移动来保证解可行性惩罚函数版则允许临时违约但将违约量计入目标让搜索可以穿越不可行区域。LNS 大邻域搜索的思路完全不同它不依赖微小的邻域移动而是每次破坏解的一部分再重新修复跳出了局部最优的困局。用一个表格来说明四种方法的特点会直观一些方法求解策略时间窗处理适用规模典型耗时CW 节约算法构造式贪心合并时检查N 100毫秒级TS 硬约束版局部搜索严格筛选N 200秒级TS 惩罚函数版局部搜索惩罚项引导N 300秒到分钟级LNS破坏/修复修复时兼顾N 500分钟级从这个对照表可以看出四种方法并不是互相替代的关系而是不同精度与算力需求下的选择序列。我在实际项目中常用的组合是CW 生成初始解LNS 做核心优化TS 惩罚函数版作为中规模算例的备选。2.3 统一的输出格式解结构、成本与可行性判定四种算法在 Matlab 中共享同一个数据结构和评估函数这样才能保证对比公平。我定义一个标准解结构体 solution包含字段 routes车辆路径数组、totalCost总距离、violation时间窗违约总量、vehicleCount使用车辆数。评估函数 evalSolution 对所有算法统一调用计算总成本并记录违约情况这样后续对比时只需要看数值不需要关心各算法内部的解表示差异。在模型基础上还需要一个随机算例生成器。生成器通过 rand 函数在 [0, 100] 平方区域内生成客户坐标时间窗则按正态分布随机生成宽度取 30 到 60 分钟不等服务时间设为 10 分钟。生成的数据保存为结构体数组包含 x、y、demand、tw_start、tw_end 五个字段。这部分代码是所有实验的输入基础。% 生成VRPTW测试算例返回客户数据结构体数组 function data generateVRPTW(nCustomers, capacity, regionSize) % nCustomers: 客户数量; capacity: 车辆容量; regionSize: 坐标范围 rng(42); % 固定随机种子保证可复现 data.capacity capacity; data.nCustomers nCustomers; data.coords rand(nCustomers, 2) * regionSize; % 坐标矩阵 data.demand randi([10, 30], nCustomers, 1); % 需求量 10-30 meanTwStart 100 rand(nCustomers, 1) * 200; % 时间窗起始 100-300 data.twStart meanTwStart; data.twEnd data.twStart randi([30, 90], nCustomers, 1); % 窗宽30-90 data.serviceTime 10 * ones(nCustomers, 1); % 服务时间固定10分钟 end这段代码生成的数据是后续所有方法共享的输入。距离矩阵在算法启动前一次性算好时间窗的检查在每个算法内部各自实现但判断标准一致——车辆到达时间早于 twStart 则等待到 twStart晚于 twEnd 则记为违约。有了这个统一的数据基础四种算法之间的对比才有意义。3. 用 CW 节约算法构建初始解并修正时间窗可行性3.1 节约值计算的 Matlab 实现细节CW 节约算法Clarke-Wright Savings Algorithm的核心思想很简单如果把两条路径合并成一条可以省下多少行驶距离。配送中心到客户 i 的距离加上配送中心到客户 j 的距离减去客户 i 直接到客户 j 的距离就是节约值 s_ij。节约值越大说明将客户 i 和 j 放在同一条路径上越有利。在 Matlab 中实现时我先把所有可能的客户对组合算一遍节约值存成 nCustomers × nCustomers 的矩阵再按降序排列。% CW节约算法计算所有客户对的节约值并降序排列 function [route, cost] cwSavings(data, distMatrix) n data.nCustomers; savings zeros(n, n); for i 1:n for j i1:n % 节省值 depot到i depot到j - i到j的距离 savings(i,j) distMatrix(1,i1) distMatrix(1,j1) - distMatrix(i1,j1); savings(j,i) savings(i,j); % 对称矩阵 end end % 将节约值转为按降序排列的列表每行[i, j, s_ij] [iIdx, jIdx] find(triu(savings)); % 取上三角避免重复 savingsList [iIdx, jIdx, savings(sub2ind(size(savings), iIdx, jIdx))]; [~, sortIdx] sort(savingsList(:,3), descend); savingsList savingsList(sortIdx, :); % 按节约值从大到小排列 % ... 后续合并逻辑省略 end这段代码的关键在于下标处理。distMatrix 的第一行第一列是配送中心到自身的距离 0所以客户 i 到配送中心的距离要取 distMatrix(1, i1)这里下标加 1 是因为 Matlab 从 1 开始索引。合并客户对时每次都要检查两个端点是否在已有路径的端点位置只有端点才能合并内部节点不能直接接续。同时需要实时更新路径的载重和时间窗信息。3.2 时间窗检查与合并条件中的容量判断CW 算法在无时间窗时只需要检查合并后载重是否超限但在 VRPTW 中还必须验证合并后的路径是否满足每个客户的时间窗。车辆在路径上的到达时间依赖于前面所有客户的累计服务时间和行驶时间因此检查时间窗时必须从头到尾推演一遍。我在这里使用一个常见的处理方式每条路径保存其起点、终点、总载重和最后到达时间。合并客户 i 和客户 j 时要求 i 是某条路径的最后一个客户即路径末端j 是另一条路径的第一个客户即路径起始端合并后的新路径时间窗检查就转换为把第二条路径整体接到第一条路径后面重新计算所有客户的到达时间。% 检查合并后的路径是否满足时间窗约束 % route1: 从左到右的客户序列route2: 要被接续的路径 function feasible checkTWForMerge(route1, route2, data, distMatrix) mergedRoute [route1, route2]; % 拼接路径 arrTime 0; % 从配送中心出发时间为0 prevNode 0; % 配送中心编号为0 feasible true; for k 1:length(mergedRoute) % 到达当前客户的时间 上一节点到达服务行驶 travelTime distMatrix(prevNode1, mergedRoute(k)1); arrTime arrTime travelTime; if arrTime data.twStart(mergedRoute(k)) arrTime data.twStart(mergedRoute(k)); % 早到等待 end if arrTime data.twEnd(mergedRoute(k)) feasible false; % 晚到违约 return; end arrTime arrTime data.serviceTime(mergedRoute(k)); % 累加服务时间 prevNode mergedRoute(k); end end这段判断逻辑是 CW 算法中最容易写错的地方。很多人直接对比节约值大小就合并忽略了时间窗的累积效应。注意等待时间的处理车辆早到可以等所以到达时间会被推移到时间窗起点但是晚到直接判为不可行不允许通过等待来弥补。在硬约束场景下这种严格检查是必要的。CW 算法跑完得到的是一个可行初始解但它只是合理而非优秀。时间窗约束的存在让贪心策略经常失效局部最省距离的合并可能导致后续客户全部违约。因此 CW 解通常作为 TS 或 LNS 的初始解而不是最终输出。4. TS 禁忌搜索的硬约束版与惩罚函数版双实现4.1 禁忌表结构与邻域操作设计禁忌搜索Tabu Search的核心机制是维护一张禁忌表避免搜索过程在几个局部最优解之间来回循环。我实现 TS 时使用两个核心邻域算子2-opt 和 relocation。2-opt 选取路径中的两段子路径反转其中一段的连接方向relocation 则将一个客户从原位置移到另一个位置插入。在每次迭代中对当前解施加这两种操作生成候选解集合从中选出目标值最佳的候选解作为新解即使新解比当前解更差只要在禁忌表之外就接受以此跳出局部最优。% TS禁忌表结构记录最近几次移动的客户对禁止反向操作 function [bestSol, histCost] tsSearch(initSol, data, distMatrix, maxIter, tabuLen, mode) % mode: hard 硬约束版 / penalty 惩罚函数版 tabuList zeros(tabuLen, 2); % 每行记录一个被禁忌的移动 tabuIdx 1; % 环形指针 currentSol initSol; bestSol initSol; bestCost evaluateSolution(bestSol, data, distMatrix); histCost zeros(maxIter, 1); % 记录每代最优成本 for iter 1:maxIter % 生成候选解集合 candidates generateCandidates(currentSol, data, distMatrix); % 选择最优且不在禁忌表中的候选或者满足藐视准则的候选 [nextSol, moveRecord, ~] selectBestCandidate(candidates, tabuList, ... bestCost, data, distMatrix, mode); % 更新禁忌表 tabuList(tabuIdx, :) moveRecord; tabuIdx mod(tabuIdx, tabuLen) 1; % 环形覆盖 currentSol nextSol; curCost evaluateSolution(currentSol, data, distMatrix, mode); if curCost bestCost bestCost curCost; bestSol currentSol; end histCost(iter) bestCost; end end这个代码框架是两版 TS 共用的骨架区别集中在 selectBestCandidate 和 evaluateSolution 两个函数里。禁忌表长度为 10 时适合 100 客户以内的算例长度过短会退化为普通局部搜索过长则可能把最优移动也禁掉。藐视准则aspiration criterion是必要组件当某个候选解虽然被禁忌但它的成本优于历史最优解那么仍可以接受这个移动避免禁忌表阻碍突破。4.2 硬约束版的评估函数与惩罚函数版的差异点硬约束版的核心原则是可行解优先。在 evaluateSolution 函数中如果解中存在任何时间窗违约就直接返回一个巨大的目标值比如总距离的 100 倍让搜索算法自动避开这些解。这种做法的问题在于搜索会被限制在可行解空间内部如果初始解本身离全局最优较远搜索很难跨越中间的不可行区域找到更优解。惩罚函数版则把不可行性转化为目标的一部分总成本 总距离 penaltyCoeff × 总时间窗违约量。这里的 penaltyCoeff 是动态调整的——如果最近几次迭代都没有出现违约解说明搜索太保守就减小系数鼓励探索反之如果违约解过多且不被接受就增大系数迫使其回到可行域。这种参数自适应机制比固定系数有效得多。% 评估解的成本硬约束版与惩罚函数版 function cost evaluateSolution(sol, data, distMatrix, mode) totalDist sol.totalDist; totalViolation sol.twViolation; % 总时间窗违约时间 if strcmp(mode, hard) if totalViolation 0.001 cost 1e10; % 违约直接罚到天文数字 else cost totalDist; end else % penalty mode penaltyCoeff sol.penaltyCoeff; % 从解结构中读取当前系数 cost totalDist penaltyCoeff * totalViolation; end end硬约束版因为可行解的严格限制在小规模问题上表现往往优于惩罚版——它不会浪费时间在不可行区域打转。但当约束较紧时间窗很窄或规模较大时可行解空间可能被切分成多个孤立区域硬约束版容易困在第一个找到的可行区域内出不来。惩罚函数版通过调节 penaltyCoeff 在迭代过程中逐步加严约束相当于使用了持续变形的目标函数景观这在实践中往往能找到更优解。5. LNS 大邻域搜索的 destroy-repair 机制在 Matlab 中的落地5.1 三种 destroy 算子随机移除、最差移除、相关移除LNS 大邻域搜索的操作分为两步destroy 阶段移除一部分客户repair 阶段重新插入被移除的客户。destroy 的质量直接影响搜索结果因为如果移除的客户太少邻域范围不够大搜索退化为小邻域搜索移除太多则修复困难解质量下降。我实现了三种 destroy 算子按比例混合使用效果最好。随机移除最简单以均匀分布随机选择 r 个客户从路径中取出。最差移除则通过计算每个客户在当前路径中的边际贡献来决定移除谁——把客户移除后总距离下降最多的优先移除。相关移除更精细计算所有客户对的距离相似度和时间窗重叠度随机选一个种子客户移除与它最相关的若干客户。这种做法的出发点是距离近且时间窗接近的客户通常可以互换位置批量移除它们才能打开重排空间。% 相关移除Related Removal的核心实现 function removed relatedRemoval(sol, numRemove, data, distMatrix) % sol当前解numRemove要移除的客户数量 custList sol.custList; % 当前解中所有客户编号列表 seedIdx randi(length(custList)); % 随机选种子客户 seed custList(seedIdx); % 计算所有其他客户与种子的相关度 relatedness zeros(length(custList), 1); for i 1:length(custList) if custList(i) seed relatedness(i) -inf; % 种子本身不参与比较 else % 距离相似度归一化到0-1之间 distRatio distMatrix(seed1, custList(i)1) / max(distMatrix(:)); % 时间窗重叠度 twOverlap min(data.twEnd(seed), data.twEnd(custList(i))) - ... max(data.twStart(seed), data.twStart(custList(i))); twRatio max(0, twOverlap) / max(data.twEnd(seed), data.twEnd(custList(i))); relatedness(i) 0.7 * distRatio 0.3 * twRatio; % 加权相关度 end end [~, sortIdx] sort(relatedness, descend); % 相关度从高到低排列 removed custList(sortIdx(1:numRemove)); % 取前numRemove个客户作为移除集合 end这三种 destroy 算子在实现时要注意参数 r 的取值。我一般把移除数量设为总客户数的 10% 到 30%太小起不到跳坑效果太大会导致 repair 阶段很难给出可行解。实际代码中还可以加入自适应机制如果连续多次迭代结果无改进就增大移除比例。5.2 repair 算子的贪心插入与后悔插入对比repair 阶段的任务是把被移除的客户重新插回路径中。最直接的贪心插入是逐个处理被移除客户每个客户找到所有可行插入位置中目标增量最小的那个位置。目标增量包括行驶距离增量和时间窗约束检查两方面非法位置直接跳过。贪心插入的缺点是后续客户的插入位置受前面客户影响太大如果前面客户占了最优位置后面的客户只能将就。更稳妥的是后悔插入Regret Insertion每个客户先找出最优插入位置和目标增量找到次优位置的增量两者之差就是这个客户的后悔值。每次迭代时优先插入后悔值最大的客户这样能让自由度最少的客户先落位。function sol regretInsertion(sol, removedList, data, distMatrix) removedQueue removedList; % 待插入客户集合 while ~isempty(removedQueue) bestRegret -inf; bestCust -1; bestPos []; for i 1:length(removedQueue) cust removedQueue(i); % 找出该客户的所有可行插入位置及增量 [deltas, positions] findInsertionPositions(sol, cust, data, distMatrix); if length(deltas) 2 % 只有一个可行位置直接记录 regretValue deltas(1); curBestPos positions(1, :); else % 最优增量与次优增量的差 regretValue deltas(2) - deltas(1); curBestPos positions(1, :); end if regretValue bestRegret bestRegret regretValue; bestCust cust; bestPos curBestPos; end end % 将后悔值最大的客户插入其最优位置 sol insertCustomer(sol, bestCust, bestPos, data, distMatrix); removedQueue(removedQueue bestCust) []; % 从队列中移除已插回的客户 end end后悔插入的计算量是贪心插入的 2 到 3 倍因为每个客户都要扫描两轮插入位置。但在时间窗约束较紧时后悔插入的成功率和解质量显著优于贪心插入因为它优先处理了没几个位置能放的困难客户避免最后无位置可插的尴尬。LNS 的 destroy 和 repair 组合起来相当于在解空间中做了一次大范围的扰动每次迭代的步幅远大于 TS 的 2-opt 移动这也是它在大规模问题上更具优势的原因。6. 四方法在统一算例上的横向对比与参数调优建议为了验证四种方法的实际表现我在标准的数据结构上做了统一实验。实验设定为 100 个客户、容量 200、时间窗宽度 50 分钟运行环境为 Matlab R2023b。每种算法独立运行 10 次取平均值记录三个指标最优总距离、首次达到最优的迭代次数、总运行时间。方法平均总距离最优总距离平均迭代次数平均运行时间CW 节约算法178417840一次构造0.35 秒TS 硬约束版153114762132.8 秒TS 惩罚函数版144813921874.1 秒LNS14051371966.7 秒这个结果说明CW 作为纯构造算法速度和精度都处于下风但它的价值在于为其他三种方法提供高质量初始解TS 硬约束版收敛较慢且容易停滞TS 惩罚函数版因为能穿越不可行区域解质量明显提升LNS 在解质量和迭代效率上都是最优的代价是单次迭代耗时更长。参数调优方面我总结出三条容易踩的坑。第一TS 的禁忌表长度与问题规模强相关100 客户时 15 到 20 是比较稳的选择太小会震荡太大会让搜索变得迟钝。第二LNS 的移除比例固定为 20% 并不是最佳配置——如果时间窗很窄移除比例降到 10% 能减少修复难度如果时间窗宽裕可以提高到 35% 来加大搜索范围。第三TS 惩罚函数版中的惩罚系数不宜固定我使用的方法是每迭代 50 次检查违约解占比若超过 30% 将系数增大 1.5 倍若低于 10% 则减小 0.8 倍这样能让搜索在可行区域和不可行区域之间保持动态平衡。% 自适应惩罚系数更新逻辑 function coeff updatePenalty(coeff, violationRatio, iter) % violationRatio: 最近窗口内的违约解占比 if mod(iter, 50) 0 if violationRatio 0.3 coeff min(coeff * 1.5, 1000); % 违约太多加大惩罚 elseif violationRatio 0.1 coeff max(coeff * 0.8, 10); % 违约太少放宽探索 end end end验证最终解是否可靠的另一个技巧是看多次运行的方差而不只是看单次最优值。LNS 包含随机过程每次运行结果有波动我通常要求同一算例跑 5 次若最优值与平均值的偏差超过 3%说明移除比例或迭代次数还不够。最后的解可以通过绘制路径图来肉眼验证逻辑合理性车辆路径不应交叉过多、时间窗紧的客户应在路径靠前位置服务。将四种方法的路径叠加画在同一张图上能直观发现 CW 结果中路径交叉明显而 LNS 结果的分组结构更清晰这也是快速判断算法优劣的实操手段。本文还有配套的精品资源点击获取
