尧图网站建设 
简介针对带时间窗的车辆路径问题VRPTW本合集提供四种Matlab求解代码包括C-W节约算法、禁忌搜索硬约束版、禁忌搜索惩罚函数版和LNS适用于物流配送路径优化和VRPTW算法的教学与实验。VRPTW需要在多客户、多车辆条件下安排行驶路线使总成本最小同时满足客户时间窗和车辆容量约束。压缩包内共372个文件主要有90个m源码文件、280个txt数据或说明文件以及2个mat数据文件整体容量约859KB便于直接测试和二次开发。代码包含初始化时间窗、构建初始路径、改进搜索等模块并设有统一入口主程序可比较不同算法在同一数据集上的求解质量与效率。已有729人学习使用适合运筹优化、智能算法方向的学生和工程师研读对比既可辅助课程设计也可作为算法设计与论文实验的参考。1. VRPTW合集四种方法不是选项是同一套数据的竞争解我见过不少Matlab代码仓库把VRPTW带时间窗的车辆路径问题的求解器拆得七零八落CW节约算法写在脚本里TS和LNS散落在函数文件里数据格式不统一跑一次实验要手工改矩阵最后对比结果还要靠肉眼对齐。这个标题里的“合集”让我最在意的不是代码量而是“四种方法对比”这个动作CW、TS硬约束版、TS惩罚函数版、LNS它们不是四个独立玩具而是同一组实例上的四套竞争解。VRPTW的难点在于时间窗和载重约束会互相牵制CW是构造式启发式TS偏局部搜索LNS则靠大邻域跳出局部最优把它们放在同一个数据接口下评估你才能真正理解哪个参数在起作用。这篇博文就按这个思路展开先建立VRPTW的数学模型和CW基线再拆TS两个变体的差异最后落到LNS及Matlab实验配置。适合正在写物流调度论文或要用Matlab做算法对比的工程师。2. VRPTW问题定义与CW节约算法的松弛解构造2.1 时间窗、载重与路径的数学模型VRPTW可以描述为车场节点0客户节点1到N每辆车最大载重Q客户i有需求量q_i、服务时间s_i、时间窗[e_i, l_i]车场也有时间窗[e_0, l_0]。目标是找到一组从车场出发并返回车场的车辆路径覆盖所有客户最小化总行驶距离或车辆数。约束包括每条路径总需求量不超过Q车辆服务每个客户的时间点必须在时间窗内早到可以等待但不能晚于截止时间每个客户只能被访问一次。数学上用x_{ijk}表示车辆k是否从i到jt_{ik}表示车辆k到达i的时间。硬时间窗约束写作e_i ≤ t_{ik} ≤ l_i软时间窗则在目标函数中加上max(0, t_{ik} - l_i)的惩罚项。CW节约算法原本是为无时间窗VRP设计的它的指导思想是利用“两个客户点在同一条路径上的距离节约”来构造解。对于VRPTW最直接的改造是在合并路径时加入时间窗可行性检查但很多Matlab实现只检查载重导致生成的解存在时间窗违反——这个问题我放在2.3详细说。2.2 CW节约算法的核心公式与Matlab实现CW的第一步是计算所有客户对(i,j)的节约值s(i,j) d(0,i) d(0,j) - d(i,j)其中d是欧氏距离或真实路网距离。然后将节约值从大到小排序依次判断是否需要合并客户i所在路径和j所在路径。合并的前提是i是所在路径的末尾节点或j是所在路径的首节点且合并后路径载重与时间窗都可行。典型的Matlab实现如下function routes cw_vrptw(d, demand, cap, e, l, s) % d: 包括车场在内的距离矩阵d(1,i)为车场到客户i距离 % demand: 客户需求量cap: 车辆容量 % e, l: 客户时间窗左端和右端数组长度为N1索引1为车场 % s: 客户服务时间与e,l同长度 % 返回值 routes: cell数组每个元素是一条路径的客户序列不含车场 n size(d,1) - 1; % 客户数 savings zeros(n*(n-1)/2, 3); idx 0; for i 2:n1 for j i1:n1 idx idx 1; savings(idx,:) [i-1, j-1, d(1,i)d(1,j)-d(i,j)]; end end savings sortrows(savings, -3); % 按节约值降序 routes {}; % 初始化每个客户为单独路径节点用客户编号1~n表示 for k 1:n routes{k} k; end % 用一个数组记录每个客户当前所在路径在routes中的索引 route_of (1:n); for k 1:size(savings,1) i savings(k,1); j savings(k,2); ri route_of(i); rj route_of(j); if ri rj continue; end % 检查i是否为ri路径的最后一个节点j是否为rj路径的第一个节点 if routes{ri}(end) i routes{rj}(1) j % 尝试合并将路径rj接到ri后面 new_route [routes{ri}, routes{rj}]; if feasible(new_route, demand, cap, e, l, s, d) routes{ri} new_route; % 更新rj中所有节点的route_of为ri for r 1:length(routes{rj}) route_of(routes{rj}(r)) ri; end routes{rj} []; end end end % 删除空路径 routes routes(~cellfun(isempty, routes)); end function ok feasible(route, demand, cap, e, l, s, d) % 检查路径载重与时间窗 if sum(demand(route)) cap ok false; return; end t 0; % 到达车场后的出发时间初始为e(1)这里设为0方便 for i 1:length(route) customer route(i); if i 1 travel d(1, customer1); else travel d(route(i-1)1, customer1); end t t travel; if t e(customer1) t e(customer1); % 等待 end if t l(customer1) ok false; return; end t t s(customer1); end % 最后返回车场 t t d(route(end)1, 1); if t l(1) ok false; return; end ok true; end这段代码的核心逻辑是先建立节约值矩阵再迭代合并。feasible函数检查载重和时间窗注意时间更新时要累加服务时间和等待时间。这里的d(1,i)约定第1个节点是车场客户从2开始。route_of数组是关键它能在O(1)时间内找到客户所在路径避免每次合并都全量搜索——不过在Matlab里直接把路径数组传引用比较麻烦所以这里用索引数组是常见做法。2.3 为什么CW解在硬时间窗下容易不可行CW节约算法的原版并不考虑到达时间。当你把两条路径合并时即使原路径各自满足时间窗合并后可能导致某些客户的到达时间延后尤其是前面的路径执行时间长后面的客户会晚到可能超过l_i。很多Matlab实现只检查载重得到的“解”实际上是可行解的上界。要提升CW质量常见做法是在计算节约值时加入时间窗松弛例如引入一个惩罚项α*max(0, e_i service_time - l_j)把“潜在时间窗冲突”计入节约值排序。这个思路在合集里很实用因为CW通常作为TS和LNS的初始解CW越差后续搜索需要“纠正”的次数越多。我在实际实验中发现对于时间窗较紧的R1类实例纯CW解的车数往往超过最优解2-3辆并且有20%的路径违反时间窗。所以在合集代码中我一般会把CW结果先交给TS硬约束版修正而不是直接当作最终解。这也是为什么标题要把CW和TS放在一起对比——单独看CW没有意义关键看它作为种子解如何影响后续搜索。3. TS硬约束版和惩罚函数版的邻域搜索实现3.1 禁忌搜索的两类罚函数设计TS禁忌搜索的核心是“在邻域中找最优移动并用禁忌表防止走回头路”。VRPTW的邻域通常包含三类客户端点转移将客户i从路径A移到路径B、客户端点交换i和j互换、2-opt局部路径逆转。对于硬约束版每次移动都必须在载重和时间窗上严格可行目标函数就是总行驶距离。问题在于当时间窗很紧时可行邻域可能非常小搜索会迅速卡住。惩罚函数版则允许移动产生时间窗或载重违反但在目标函数上加权惩罚项。常见的做法是penalty_i w_e * max(0, e_i - t_i) w_l * max(0, t_i - l_i) w_q * max(0, load - cap)这里w_e、w_l、w_q是惩罚系数。搜索过程中如果连续若干次迭代都没有找到更优可行解则调大惩罚系数如果可行解改善明显则调小。这样做的效果是搜索在初期能“跨过”不可行区域后期逐步把解拉回可行域。我在Matlab实现中会维护一个penalty向量每20次迭代根据可行解比例更新一次。3.2 邻域算子与禁忌表在Matlab中的实现下面给出TS硬约束版的核心邻域搜索代码惩罚函数版的差异会在3.3节说明。这里用路径集合routescell数组和禁忌表tabu表示。禁忌表是一个N×N的矩阵tabu(i,j)记录客户i到客户j的转移被禁忌的剩余迭代次数。function [routes, bestDist] ts_vrptw(init_routes, d, demand, cap, e, l, s, options) % options: struct包含maxIter, tabuLen, neighborSize, usePenalty等 % 硬约束版 usePenaltyfalse routes init_routes; bestRoutes routes; bestDist evaluate(routes, d, e, l, s, 0); % 硬约束假设都可行 tabu zeros(size(d,1)-1); for iter 1:options.maxIter % 生成候选邻域随机采样neighborSize个移动 [moves, moveDist] generateNeighbors(routes, d, demand, cap, e, l, s, options.neighborSize); % 筛选出非禁忌或满足蔑视规则的移动 bestMove -1; bestCost inf; for k 1:length(moves) mv moves{k}; % mv [type, i, j, routeI, routeJ] if tabu(mv(2), mv(3)) 0 % 禁忌但若优于当前全局最优则蔑视 if moveDist(k) bestDist continue; end end if moveDist(k) bestCost bestCost moveDist(k); bestMove k; end end if bestMove -1 continue; end % 执行移动 applyMove(routes, moves{bestMove}); % 更新禁忌表将被替换的旧邻接关系加入禁忌 updateTabu(tabu, moves{bestMove}, options.tabuLen); % 更新最优解 totalDist evaluate(routes, d, e, l, s, 0); if totalDist bestDist bestDist totalDist; bestRoutes routes; end end routes bestRoutes; endgenerateNeighbors函数会随机选两个客户判断能否插入对方的路径并基于时间窗检查可行性。硬约束版直接调用2.2中的feasible函数。为了提高速度我通常会缩小neighborSize到20或30因为完整枚举所有交换在客户超过50后非常慢。禁忌长度tabuLen设为sqrt(N)附近比如N25时取5N100时取10经验效果不错。这里的关键参数是neighborSize和maxIter。设太小会退化成交替地做局部搜索设太大每次迭代都扫描大量候选导致运行时间成倍增加。一般我会先用100次迭代跑通观察bestDist是否还有下降趋势再决定是否加大迭代次数。3.3 TS硬约束版和惩罚函数版在Matlab代码中的关键差异硬约束版和惩罚函数版的差别并不在TS框架而在evaluate函数和feasible检查上。硬约束版只计算距离但要求邻域生成时所有移动都通过feasible惩罚函数版则把feasible拆成两部分载重违反和时间窗违反并计算惩罚值。function cost evaluate_penalty(routes, d, demand, cap, e, l, s, w) % w [w_load, w_time_early, w_time_late, w_vehicle] cost 0; loadPen 0; timePen 0; for r 1:length(routes) route routes{r}; if isempty(route), continue; end cost cost d(1, route(1)1) d(route(end)1, 1); for k 1:length(route)-1 cost cost d(route(k)1, route(k1)1); end % 载重检查 load sum(demand(route)); if load cap loadPen loadPen w(1) * (load - cap); end % 时间窗检查 t 0; for k 1:length(route) customer route(k); t t (d(route(k-1)1, customer1) if k1 else d(1, customer1)); if t e(customer1) timePen timePen w(2) * (e(customer1) - t); t e(customer1); elseif t l(customer1) timePen timePen w(3) * (t - l(customer1)); end t t s(customer1); end end cost cost loadPen timePen w(4) * length(routes); end注意上段代码中我写了route(k-1)在k1时应使用车场这是索引偏移的问题实际实现时可以用from 1 if k1 else route(k-1)1避免错误。惩罚系数初始取值很有讲究。w_load设得过大搜索一开始就不敢超载效果等同硬约束设得过小解长时间停留在不可行区域。我的做法是w_load初始为0.5 * (平均距离/平均载重)w_time_late初始为平均距离/平均时间窗晚点然后每50次迭代检查不可行解比例若超过30%则乘以1.5若低于10%则除以1.2。这样可以让搜索动态平衡可行域探索。对比来看硬约束版胜在无参但解容易卡在局部最优惩罚函数版多出3-4个参数调好了能拿到比硬约束低3%-8%的总距离但参数错配时反而更差。4. LNS破坏与修复机制以及四种方法的对比维度4.1 LNS的destroy与repair思路LNS大邻域搜索与TS最大的不同在于它不枚举小邻域而是每次迭代用破坏算子从当前解中移除一批客户再用修复算子将它们重新插入。移除数量通常占总客户数的20%-40%目的是让搜索在更大的解空间中进行。常见的破坏算子有随机移除、相关移除距离相近的客户一起移除和路径移除整条路径移除。修复算子则是最小插入代价的贪心对每个待插入客户找最佳插入位置使距离增量最小且满足时间窗。下面是LNS破坏/修复的一个简化Matlab片段function [routes, dist] lns_vrptw(init_routes, d, demand, cap, e, l, s, options) routes init_routes; bestRoutes routes; dist evaluate(routes, d, e, l, s, 0); for iter 1:options.maxIter % 破坏随机移除options.removeRatio*n个客户 n size(d,1)-1; removeNum max(2, floor(n * options.removeRatio)); allCustomers 1:n; removed allCustomers(randperm(n, removeNum)); for cust removed % 找到所在路径并移除该点 for r 1:length(routes) pos find(routes{r} cust); if ~isempty(pos) routes{r}(pos) []; if isempty(routes{r}) routes(r) []; end break; end end end % 修复按照最小插入代价依次插入removed中的客户 while ~isempty(removed) bestCost inf; bestPos -1; bestRoute -1; for i 1:length(removed) cust removed(i); for r 1:length(routes) for p 0:length(routes{r}) % 在p位置插入 newRoute [routes{r}(1:p), cust, routes{r}(p1:end)]; if feasible(newRoute, demand, cap, e, l, s, d) inc insertionCost(newRoute, d, p); if inc bestCost bestCost inc; bestPos p; bestRoute r; bestCust cust; end end end end end % 若所有位置不可行则开新路径 if bestRoute -1 routes{end1} bestCust; else routes{bestRoute} [routes{bestRoute}(1:bestPos), bestCust, routes{bestRoute}(bestPos1:end)]; end removed(removed bestCust) []; end newDist evaluate(routes, d, e, l, s, 0); if newDist dist dist newDist; bestRoutes routes; end end routes bestRoutes; end这个实现里移除比例removeRatio是最重要的参数。设小了LNS退化成局部搜索设大了每次都要重新插入大量客户耗时倍增。对于25个客户的实例我一般设0.3也就是每次移除8个客户左右对于100客户实例设0.2即可。4.2 四种方法对比的五个维度把CW、TS硬约束、TS惩罚、LNS放在一张表里看才能体现出为什么合集的“对比”是有意义的。表中数据基于我常用的一组Solomon测试实例但你需要用自己的数据复现。对比维度CW节约算法TS硬约束版TS惩罚函数版LNS初始解质量差常有时间窗违反依赖CW可在100次迭代内修正同左但能探索更多不可行区域依赖初始解但破坏修复后提升明显解质量相对最优解15%-25%差距8%-12%差距5%-10%差距3%-7%差距计算时间毫秒级秒级秒级比硬约束多30%分钟级取决于迭代次数参数数量1个是否检查时间窗3-4个迭代数、禁忌长度、邻域大小6-7个4-5个实现复杂度低中中高高适用场景作为初始解生成器时间窗宽松需要快速结果时间窗紧但没时间写LNS有时间调参追求高质量解CW最大的优势是“快”和“简单”在合集里它通常是TS/LNS的起点。TS硬约束实现最稳不会出现不可行解但在时间窗紧的R型实例上容易陷入局部最优。TS惩罚版通过调参数能拿到比硬约束更优的解但调参本身是个工程师经验活。LNS在解质量上往往最强但它的耗时是TS的几十倍而且对removeRatio和迭代次数非常敏感。4.3 从CW到TS到LNS的迭代逻辑我在对比实验中发现一个规律当时间窗宽度从宽松变紧时四种方法的解质量差距会拉大。原因在于CW完全忽略了时间窗TS硬约束在紧时间窗下可行邻域缩小而TS惩罚版和LNS能通过“暂时不可行”的方式跳出局部最优。这解释了为什么标题要把TS分成两个版本——两者的行为差异足够大值得并列。LNS之所以效果更好是因为它每次移除20%-30%的客户相当于在解的“骨架”上重新构建。这比TS的小邻域更能改变路径结构。但代价是每个插入步骤都要调用feasible复杂度高所以需要在Matlab里对热循环做向量化否则100客户端点会非常慢。5. 用Matlab跑通VRPTW合集的实验设置与参数调优5.1 基于标准实例的数据结构与实验配置在合集中我推荐使用基类结构体来统一输入。数据集选Solomon的25客户实例如R101、C101、RC101比较合适因为四种方法在25客户规模下都能在几秒到几分钟内跑完适合调参。数据结构如下vrp struct(); vrp.node 0; % 车场 vrp.customers (1:25); vrp.x [d_x; cust_x]; % 车场和客户坐标 vrp.y [d_y; cust_y]; vrp.demand [0; demands]; vrp.service [0; service_times]; vrp.earliest [0; e_values]; vrp.latest [1440; l_values]; % 单日时间窗 vrp.capacity 200; vrp.dist squareform(pdist([vrp.x, vrp.y]));这样设计的好处是四个函数都可以接收同一个vrp结构体返回统一格式的solution结构体solution.routes routes; solution.distance totalDistance; solution.vehicles length(routes); solution.time elapsedTime; solution.feasible checkAllFeasible(routes, vrp);实验至少做5次独立运行每次使用不同的随机数种子取平均值。因为TS和LNS都含随机因素单次结果不能代表算法真实水平。5.2 四种方法的参数表与主脚本骨架下面给出一个精简主脚本调用四个方法并收集结果。options结构体为每个算法单独提供参数。%% VRPTW合集主脚本 load(R101_25.mat); % 加载实例其中包含vrp结构体 numRuns 5; results struct(); for run 1:numRuns rng(run * 17); % 固定种子 % 1. CW节约算法 tic; routes_cw cw_vrptw(vrp.dist, vrp.demand, vrp.capacity, vrp.earliest, vrp.latest, vrp.service); res_cw evaluateSolution(routes_cw, vrp); res_cw.time toc; % 2. TS硬约束版 opts_tsh.maxIter 300; opts_tsh.tabuLen 5; opts_tsh.neighborSize 25; opts_tsh.usePenalty false; tic; routes_tsh ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsh); res_tsh evaluateSolution(routes_tsh, vrp); res_tsh.time toc; % 3. TS惩罚函数版 opts_tsp opts_tsh; opts_tsp.usePenalty true; opts_tsp.wLoad 0.5; opts_tsp.wLate 1.0; opts_tsp.wEarly 0.1; tic; routes_tsp ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsp); res_tsp evaluateSolution(routes_tsp, vrp); res_tsp.time toc; % 4. LNS opts_lns.maxIter 100; opts_lns.removeRatio 0.3; tic; routes_lns lns_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_lns); res_lns evaluateSolution(routes_lns, vrp); res_lns.time toc; % 记录本次运行 results(run).cw res_cw; results(run).tsh res_tsh; results(run).tsp res_tsp; results(run).lns res_lns; end % 汇总统计 summarizeResults(results);参数表中TS硬约束版只需要迭代次数、禁忌长度和邻域大小。惩罚函数版额外需要三个惩罚权重。LNS需要迭代次数和移除比例。这些参数的推荐初始值我总结在下方表格里但实际使用时还应结合具体时间窗宽度调整。算法参数名推荐初始值调整提示CW无-如果时间窗紧可在节约值中加入时间窗惩罚TS硬约束maxIter300若50次迭代无改进可加大到1000TS硬约束tabuLensqrt(N)N25时取5N100时取10TS硬约束neighborSize25客户多时设为2*N但会变慢TS惩罚wLoad / wLate / wEarly0.5 / 1.0 / 0.1若不可行解比例高提高wLateLNSmaxIter100每次运行时间可接受时加大到500LNSremoveRatio0.3时间窗紧时用0.2宽松时用0.45.3 结果对比与可视化技巧跑完实验后我会把results结构体转成表格并用Matlab的bar或boxchart展示。更直观的是画路径图和收敛曲线。路径图可以用plot按车场出发再返回的顺序绘制figure; hold on; colors lines(length(routes)); for r 1:length(routes) route routes{r}; xSeq [vrp.x(1); vrp.x(route1); vrp.x(1)]; ySeq [vrp.y(1); vrp.y(route1); vrp.y(1)]; plot(xSeq, ySeq, o-, Color, colors(r,:)); end hold off;收敛曲线则需要TS和LNS在迭代过程中记录bestDist。我在ts_vrptw函数中增加了一个history字段每次发现新最优解时记录当前迭代和距离最后画成阶梯线。这个图能直接看出LNS在中后期仍在改进而TS硬约束在20代后基本停滞。另一个实用技巧是用writetable把多种方法的结果导出到Excel方便后续用gather做显著性分析。表格至少包含列运行序号、算法名、车辆数、总距离、运行时间、是否可行。跨多次运行后你还能计算均值和标准差这对写论文非常重要。6. 从结果反推参数为你的VRPTW实例选型合集的最后一步不是跑完四种方法就结束而是要根据一次快速实验的结果反推后续参数策略。我会在用户数据集上先只运行CW和LNS各一次比较两者总距离的相对差距gap (lnsDist - cwDist) / cwDist。这个gap能反映时间窗约束的松紧如果gap小于5%说明CW已经逼近LNS你的实例大概率时间窗很宽松用TS硬约束版加上200次迭代就足够没必要跑LNS如果gap超过15%说明时间窗在强约束直接用LNS并加大迭代次数更划算。下面这段小脚本用来做这个决策% 假设已经得到 cwDist 和 lnsDist gap (lnsDist - cwDist) / cwDist; if gap 0.05 opts_tsh.maxIter 200; final ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsh); elseif gap 0.15 opts_tsp.wLate 1.5; % 加大晚到惩罚更容易跨过不可行区域 final ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsp); else opts_lns.maxIter 500; opts_lns.removeRatio 0.25; final lns_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_lns); end除了这个规则我还经常用“验证时间窗松紧”的辅助手段计算所有客户时间窗宽度的中位数。中位数越大TS硬约束版越值得信任中位数小于平均服务时间时硬约束版容易找不到合适插入位置这时候惩罚函数版的迭代后期要额外加一个“修复步骤”——把目标函数中的惩罚权重调大让解收敛到可行域。最后一个验证技巧是对同一组参数用不同随机种子跑5次用friedman测试比较四种方法是否有显著性差异。Matlab的friedman函数直接接受一个矩阵行是算法列是实例返回p值。如果p0.05再画multcompare图看哪两个方法差异显著。这样你写在论文里的“对比”才有统计依据而不是只报一个最好平均值。本文还有配套的精品资源点击获取 VRPTW求解方法对比：CW节约算法、禁忌搜索与LNS的Matlab实现
简介针对带时间窗的车辆路径问题VRPTW本合集提供四种Matlab求解代码包括C-W节约算法、禁忌搜索硬约束版、禁忌搜索惩罚函数版和LNS适用于物流配送路径优化和VRPTW算法的教学与实验。VRPTW需要在多客户、多车辆条件下安排行驶路线使总成本最小同时满足客户时间窗和车辆容量约束。压缩包内共372个文件主要有90个m源码文件、280个txt数据或说明文件以及2个mat数据文件整体容量约859KB便于直接测试和二次开发。代码包含初始化时间窗、构建初始路径、改进搜索等模块并设有统一入口主程序可比较不同算法在同一数据集上的求解质量与效率。已有729人学习使用适合运筹优化、智能算法方向的学生和工程师研读对比既可辅助课程设计也可作为算法设计与论文实验的参考。1. VRPTW合集四种方法不是选项是同一套数据的竞争解我见过不少Matlab代码仓库把VRPTW带时间窗的车辆路径问题的求解器拆得七零八落CW节约算法写在脚本里TS和LNS散落在函数文件里数据格式不统一跑一次实验要手工改矩阵最后对比结果还要靠肉眼对齐。这个标题里的“合集”让我最在意的不是代码量而是“四种方法对比”这个动作CW、TS硬约束版、TS惩罚函数版、LNS它们不是四个独立玩具而是同一组实例上的四套竞争解。VRPTW的难点在于时间窗和载重约束会互相牵制CW是构造式启发式TS偏局部搜索LNS则靠大邻域跳出局部最优把它们放在同一个数据接口下评估你才能真正理解哪个参数在起作用。这篇博文就按这个思路展开先建立VRPTW的数学模型和CW基线再拆TS两个变体的差异最后落到LNS及Matlab实验配置。适合正在写物流调度论文或要用Matlab做算法对比的工程师。2. VRPTW问题定义与CW节约算法的松弛解构造2.1 时间窗、载重与路径的数学模型VRPTW可以描述为车场节点0客户节点1到N每辆车最大载重Q客户i有需求量q_i、服务时间s_i、时间窗[e_i, l_i]车场也有时间窗[e_0, l_0]。目标是找到一组从车场出发并返回车场的车辆路径覆盖所有客户最小化总行驶距离或车辆数。约束包括每条路径总需求量不超过Q车辆服务每个客户的时间点必须在时间窗内早到可以等待但不能晚于截止时间每个客户只能被访问一次。数学上用x_{ijk}表示车辆k是否从i到jt_{ik}表示车辆k到达i的时间。硬时间窗约束写作e_i ≤ t_{ik} ≤ l_i软时间窗则在目标函数中加上max(0, t_{ik} - l_i)的惩罚项。CW节约算法原本是为无时间窗VRP设计的它的指导思想是利用“两个客户点在同一条路径上的距离节约”来构造解。对于VRPTW最直接的改造是在合并路径时加入时间窗可行性检查但很多Matlab实现只检查载重导致生成的解存在时间窗违反——这个问题我放在2.3详细说。2.2 CW节约算法的核心公式与Matlab实现CW的第一步是计算所有客户对(i,j)的节约值s(i,j) d(0,i) d(0,j) - d(i,j)其中d是欧氏距离或真实路网距离。然后将节约值从大到小排序依次判断是否需要合并客户i所在路径和j所在路径。合并的前提是i是所在路径的末尾节点或j是所在路径的首节点且合并后路径载重与时间窗都可行。典型的Matlab实现如下function routes cw_vrptw(d, demand, cap, e, l, s) % d: 包括车场在内的距离矩阵d(1,i)为车场到客户i距离 % demand: 客户需求量cap: 车辆容量 % e, l: 客户时间窗左端和右端数组长度为N1索引1为车场 % s: 客户服务时间与e,l同长度 % 返回值 routes: cell数组每个元素是一条路径的客户序列不含车场 n size(d,1) - 1; % 客户数 savings zeros(n*(n-1)/2, 3); idx 0; for i 2:n1 for j i1:n1 idx idx 1; savings(idx,:) [i-1, j-1, d(1,i)d(1,j)-d(i,j)]; end end savings sortrows(savings, -3); % 按节约值降序 routes {}; % 初始化每个客户为单独路径节点用客户编号1~n表示 for k 1:n routes{k} k; end % 用一个数组记录每个客户当前所在路径在routes中的索引 route_of (1:n); for k 1:size(savings,1) i savings(k,1); j savings(k,2); ri route_of(i); rj route_of(j); if ri rj continue; end % 检查i是否为ri路径的最后一个节点j是否为rj路径的第一个节点 if routes{ri}(end) i routes{rj}(1) j % 尝试合并将路径rj接到ri后面 new_route [routes{ri}, routes{rj}]; if feasible(new_route, demand, cap, e, l, s, d) routes{ri} new_route; % 更新rj中所有节点的route_of为ri for r 1:length(routes{rj}) route_of(routes{rj}(r)) ri; end routes{rj} []; end end end % 删除空路径 routes routes(~cellfun(isempty, routes)); end function ok feasible(route, demand, cap, e, l, s, d) % 检查路径载重与时间窗 if sum(demand(route)) cap ok false; return; end t 0; % 到达车场后的出发时间初始为e(1)这里设为0方便 for i 1:length(route) customer route(i); if i 1 travel d(1, customer1); else travel d(route(i-1)1, customer1); end t t travel; if t e(customer1) t e(customer1); % 等待 end if t l(customer1) ok false; return; end t t s(customer1); end % 最后返回车场 t t d(route(end)1, 1); if t l(1) ok false; return; end ok true; end这段代码的核心逻辑是先建立节约值矩阵再迭代合并。feasible函数检查载重和时间窗注意时间更新时要累加服务时间和等待时间。这里的d(1,i)约定第1个节点是车场客户从2开始。route_of数组是关键它能在O(1)时间内找到客户所在路径避免每次合并都全量搜索——不过在Matlab里直接把路径数组传引用比较麻烦所以这里用索引数组是常见做法。2.3 为什么CW解在硬时间窗下容易不可行CW节约算法的原版并不考虑到达时间。当你把两条路径合并时即使原路径各自满足时间窗合并后可能导致某些客户的到达时间延后尤其是前面的路径执行时间长后面的客户会晚到可能超过l_i。很多Matlab实现只检查载重得到的“解”实际上是可行解的上界。要提升CW质量常见做法是在计算节约值时加入时间窗松弛例如引入一个惩罚项α*max(0, e_i service_time - l_j)把“潜在时间窗冲突”计入节约值排序。这个思路在合集里很实用因为CW通常作为TS和LNS的初始解CW越差后续搜索需要“纠正”的次数越多。我在实际实验中发现对于时间窗较紧的R1类实例纯CW解的车数往往超过最优解2-3辆并且有20%的路径违反时间窗。所以在合集代码中我一般会把CW结果先交给TS硬约束版修正而不是直接当作最终解。这也是为什么标题要把CW和TS放在一起对比——单独看CW没有意义关键看它作为种子解如何影响后续搜索。3. TS硬约束版和惩罚函数版的邻域搜索实现3.1 禁忌搜索的两类罚函数设计TS禁忌搜索的核心是“在邻域中找最优移动并用禁忌表防止走回头路”。VRPTW的邻域通常包含三类客户端点转移将客户i从路径A移到路径B、客户端点交换i和j互换、2-opt局部路径逆转。对于硬约束版每次移动都必须在载重和时间窗上严格可行目标函数就是总行驶距离。问题在于当时间窗很紧时可行邻域可能非常小搜索会迅速卡住。惩罚函数版则允许移动产生时间窗或载重违反但在目标函数上加权惩罚项。常见的做法是penalty_i w_e * max(0, e_i - t_i) w_l * max(0, t_i - l_i) w_q * max(0, load - cap)这里w_e、w_l、w_q是惩罚系数。搜索过程中如果连续若干次迭代都没有找到更优可行解则调大惩罚系数如果可行解改善明显则调小。这样做的效果是搜索在初期能“跨过”不可行区域后期逐步把解拉回可行域。我在Matlab实现中会维护一个penalty向量每20次迭代根据可行解比例更新一次。3.2 邻域算子与禁忌表在Matlab中的实现下面给出TS硬约束版的核心邻域搜索代码惩罚函数版的差异会在3.3节说明。这里用路径集合routescell数组和禁忌表tabu表示。禁忌表是一个N×N的矩阵tabu(i,j)记录客户i到客户j的转移被禁忌的剩余迭代次数。function [routes, bestDist] ts_vrptw(init_routes, d, demand, cap, e, l, s, options) % options: struct包含maxIter, tabuLen, neighborSize, usePenalty等 % 硬约束版 usePenaltyfalse routes init_routes; bestRoutes routes; bestDist evaluate(routes, d, e, l, s, 0); % 硬约束假设都可行 tabu zeros(size(d,1)-1); for iter 1:options.maxIter % 生成候选邻域随机采样neighborSize个移动 [moves, moveDist] generateNeighbors(routes, d, demand, cap, e, l, s, options.neighborSize); % 筛选出非禁忌或满足蔑视规则的移动 bestMove -1; bestCost inf; for k 1:length(moves) mv moves{k}; % mv [type, i, j, routeI, routeJ] if tabu(mv(2), mv(3)) 0 % 禁忌但若优于当前全局最优则蔑视 if moveDist(k) bestDist continue; end end if moveDist(k) bestCost bestCost moveDist(k); bestMove k; end end if bestMove -1 continue; end % 执行移动 applyMove(routes, moves{bestMove}); % 更新禁忌表将被替换的旧邻接关系加入禁忌 updateTabu(tabu, moves{bestMove}, options.tabuLen); % 更新最优解 totalDist evaluate(routes, d, e, l, s, 0); if totalDist bestDist bestDist totalDist; bestRoutes routes; end end routes bestRoutes; endgenerateNeighbors函数会随机选两个客户判断能否插入对方的路径并基于时间窗检查可行性。硬约束版直接调用2.2中的feasible函数。为了提高速度我通常会缩小neighborSize到20或30因为完整枚举所有交换在客户超过50后非常慢。禁忌长度tabuLen设为sqrt(N)附近比如N25时取5N100时取10经验效果不错。这里的关键参数是neighborSize和maxIter。设太小会退化成交替地做局部搜索设太大每次迭代都扫描大量候选导致运行时间成倍增加。一般我会先用100次迭代跑通观察bestDist是否还有下降趋势再决定是否加大迭代次数。3.3 TS硬约束版和惩罚函数版在Matlab代码中的关键差异硬约束版和惩罚函数版的差别并不在TS框架而在evaluate函数和feasible检查上。硬约束版只计算距离但要求邻域生成时所有移动都通过feasible惩罚函数版则把feasible拆成两部分载重违反和时间窗违反并计算惩罚值。function cost evaluate_penalty(routes, d, demand, cap, e, l, s, w) % w [w_load, w_time_early, w_time_late, w_vehicle] cost 0; loadPen 0; timePen 0; for r 1:length(routes) route routes{r}; if isempty(route), continue; end cost cost d(1, route(1)1) d(route(end)1, 1); for k 1:length(route)-1 cost cost d(route(k)1, route(k1)1); end % 载重检查 load sum(demand(route)); if load cap loadPen loadPen w(1) * (load - cap); end % 时间窗检查 t 0; for k 1:length(route) customer route(k); t t (d(route(k-1)1, customer1) if k1 else d(1, customer1)); if t e(customer1) timePen timePen w(2) * (e(customer1) - t); t e(customer1); elseif t l(customer1) timePen timePen w(3) * (t - l(customer1)); end t t s(customer1); end end cost cost loadPen timePen w(4) * length(routes); end注意上段代码中我写了route(k-1)在k1时应使用车场这是索引偏移的问题实际实现时可以用from 1 if k1 else route(k-1)1避免错误。惩罚系数初始取值很有讲究。w_load设得过大搜索一开始就不敢超载效果等同硬约束设得过小解长时间停留在不可行区域。我的做法是w_load初始为0.5 * (平均距离/平均载重)w_time_late初始为平均距离/平均时间窗晚点然后每50次迭代检查不可行解比例若超过30%则乘以1.5若低于10%则除以1.2。这样可以让搜索动态平衡可行域探索。对比来看硬约束版胜在无参但解容易卡在局部最优惩罚函数版多出3-4个参数调好了能拿到比硬约束低3%-8%的总距离但参数错配时反而更差。4. LNS破坏与修复机制以及四种方法的对比维度4.1 LNS的destroy与repair思路LNS大邻域搜索与TS最大的不同在于它不枚举小邻域而是每次迭代用破坏算子从当前解中移除一批客户再用修复算子将它们重新插入。移除数量通常占总客户数的20%-40%目的是让搜索在更大的解空间中进行。常见的破坏算子有随机移除、相关移除距离相近的客户一起移除和路径移除整条路径移除。修复算子则是最小插入代价的贪心对每个待插入客户找最佳插入位置使距离增量最小且满足时间窗。下面是LNS破坏/修复的一个简化Matlab片段function [routes, dist] lns_vrptw(init_routes, d, demand, cap, e, l, s, options) routes init_routes; bestRoutes routes; dist evaluate(routes, d, e, l, s, 0); for iter 1:options.maxIter % 破坏随机移除options.removeRatio*n个客户 n size(d,1)-1; removeNum max(2, floor(n * options.removeRatio)); allCustomers 1:n; removed allCustomers(randperm(n, removeNum)); for cust removed % 找到所在路径并移除该点 for r 1:length(routes) pos find(routes{r} cust); if ~isempty(pos) routes{r}(pos) []; if isempty(routes{r}) routes(r) []; end break; end end end % 修复按照最小插入代价依次插入removed中的客户 while ~isempty(removed) bestCost inf; bestPos -1; bestRoute -1; for i 1:length(removed) cust removed(i); for r 1:length(routes) for p 0:length(routes{r}) % 在p位置插入 newRoute [routes{r}(1:p), cust, routes{r}(p1:end)]; if feasible(newRoute, demand, cap, e, l, s, d) inc insertionCost(newRoute, d, p); if inc bestCost bestCost inc; bestPos p; bestRoute r; bestCust cust; end end end end end % 若所有位置不可行则开新路径 if bestRoute -1 routes{end1} bestCust; else routes{bestRoute} [routes{bestRoute}(1:bestPos), bestCust, routes{bestRoute}(bestPos1:end)]; end removed(removed bestCust) []; end newDist evaluate(routes, d, e, l, s, 0); if newDist dist dist newDist; bestRoutes routes; end end routes bestRoutes; end这个实现里移除比例removeRatio是最重要的参数。设小了LNS退化成局部搜索设大了每次都要重新插入大量客户耗时倍增。对于25个客户的实例我一般设0.3也就是每次移除8个客户左右对于100客户实例设0.2即可。4.2 四种方法对比的五个维度把CW、TS硬约束、TS惩罚、LNS放在一张表里看才能体现出为什么合集的“对比”是有意义的。表中数据基于我常用的一组Solomon测试实例但你需要用自己的数据复现。对比维度CW节约算法TS硬约束版TS惩罚函数版LNS初始解质量差常有时间窗违反依赖CW可在100次迭代内修正同左但能探索更多不可行区域依赖初始解但破坏修复后提升明显解质量相对最优解15%-25%差距8%-12%差距5%-10%差距3%-7%差距计算时间毫秒级秒级秒级比硬约束多30%分钟级取决于迭代次数参数数量1个是否检查时间窗3-4个迭代数、禁忌长度、邻域大小6-7个4-5个实现复杂度低中中高高适用场景作为初始解生成器时间窗宽松需要快速结果时间窗紧但没时间写LNS有时间调参追求高质量解CW最大的优势是“快”和“简单”在合集里它通常是TS/LNS的起点。TS硬约束实现最稳不会出现不可行解但在时间窗紧的R型实例上容易陷入局部最优。TS惩罚版通过调参数能拿到比硬约束更优的解但调参本身是个工程师经验活。LNS在解质量上往往最强但它的耗时是TS的几十倍而且对removeRatio和迭代次数非常敏感。4.3 从CW到TS到LNS的迭代逻辑我在对比实验中发现一个规律当时间窗宽度从宽松变紧时四种方法的解质量差距会拉大。原因在于CW完全忽略了时间窗TS硬约束在紧时间窗下可行邻域缩小而TS惩罚版和LNS能通过“暂时不可行”的方式跳出局部最优。这解释了为什么标题要把TS分成两个版本——两者的行为差异足够大值得并列。LNS之所以效果更好是因为它每次移除20%-30%的客户相当于在解的“骨架”上重新构建。这比TS的小邻域更能改变路径结构。但代价是每个插入步骤都要调用feasible复杂度高所以需要在Matlab里对热循环做向量化否则100客户端点会非常慢。5. 用Matlab跑通VRPTW合集的实验设置与参数调优5.1 基于标准实例的数据结构与实验配置在合集中我推荐使用基类结构体来统一输入。数据集选Solomon的25客户实例如R101、C101、RC101比较合适因为四种方法在25客户规模下都能在几秒到几分钟内跑完适合调参。数据结构如下vrp struct(); vrp.node 0; % 车场 vrp.customers (1:25); vrp.x [d_x; cust_x]; % 车场和客户坐标 vrp.y [d_y; cust_y]; vrp.demand [0; demands]; vrp.service [0; service_times]; vrp.earliest [0; e_values]; vrp.latest [1440; l_values]; % 单日时间窗 vrp.capacity 200; vrp.dist squareform(pdist([vrp.x, vrp.y]));这样设计的好处是四个函数都可以接收同一个vrp结构体返回统一格式的solution结构体solution.routes routes; solution.distance totalDistance; solution.vehicles length(routes); solution.time elapsedTime; solution.feasible checkAllFeasible(routes, vrp);实验至少做5次独立运行每次使用不同的随机数种子取平均值。因为TS和LNS都含随机因素单次结果不能代表算法真实水平。5.2 四种方法的参数表与主脚本骨架下面给出一个精简主脚本调用四个方法并收集结果。options结构体为每个算法单独提供参数。%% VRPTW合集主脚本 load(R101_25.mat); % 加载实例其中包含vrp结构体 numRuns 5; results struct(); for run 1:numRuns rng(run * 17); % 固定种子 % 1. CW节约算法 tic; routes_cw cw_vrptw(vrp.dist, vrp.demand, vrp.capacity, vrp.earliest, vrp.latest, vrp.service); res_cw evaluateSolution(routes_cw, vrp); res_cw.time toc; % 2. TS硬约束版 opts_tsh.maxIter 300; opts_tsh.tabuLen 5; opts_tsh.neighborSize 25; opts_tsh.usePenalty false; tic; routes_tsh ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsh); res_tsh evaluateSolution(routes_tsh, vrp); res_tsh.time toc; % 3. TS惩罚函数版 opts_tsp opts_tsh; opts_tsp.usePenalty true; opts_tsp.wLoad 0.5; opts_tsp.wLate 1.0; opts_tsp.wEarly 0.1; tic; routes_tsp ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsp); res_tsp evaluateSolution(routes_tsp, vrp); res_tsp.time toc; % 4. LNS opts_lns.maxIter 100; opts_lns.removeRatio 0.3; tic; routes_lns lns_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_lns); res_lns evaluateSolution(routes_lns, vrp); res_lns.time toc; % 记录本次运行 results(run).cw res_cw; results(run).tsh res_tsh; results(run).tsp res_tsp; results(run).lns res_lns; end % 汇总统计 summarizeResults(results);参数表中TS硬约束版只需要迭代次数、禁忌长度和邻域大小。惩罚函数版额外需要三个惩罚权重。LNS需要迭代次数和移除比例。这些参数的推荐初始值我总结在下方表格里但实际使用时还应结合具体时间窗宽度调整。算法参数名推荐初始值调整提示CW无-如果时间窗紧可在节约值中加入时间窗惩罚TS硬约束maxIter300若50次迭代无改进可加大到1000TS硬约束tabuLensqrt(N)N25时取5N100时取10TS硬约束neighborSize25客户多时设为2*N但会变慢TS惩罚wLoad / wLate / wEarly0.5 / 1.0 / 0.1若不可行解比例高提高wLateLNSmaxIter100每次运行时间可接受时加大到500LNSremoveRatio0.3时间窗紧时用0.2宽松时用0.45.3 结果对比与可视化技巧跑完实验后我会把results结构体转成表格并用Matlab的bar或boxchart展示。更直观的是画路径图和收敛曲线。路径图可以用plot按车场出发再返回的顺序绘制figure; hold on; colors lines(length(routes)); for r 1:length(routes) route routes{r}; xSeq [vrp.x(1); vrp.x(route1); vrp.x(1)]; ySeq [vrp.y(1); vrp.y(route1); vrp.y(1)]; plot(xSeq, ySeq, o-, Color, colors(r,:)); end hold off;收敛曲线则需要TS和LNS在迭代过程中记录bestDist。我在ts_vrptw函数中增加了一个history字段每次发现新最优解时记录当前迭代和距离最后画成阶梯线。这个图能直接看出LNS在中后期仍在改进而TS硬约束在20代后基本停滞。另一个实用技巧是用writetable把多种方法的结果导出到Excel方便后续用gather做显著性分析。表格至少包含列运行序号、算法名、车辆数、总距离、运行时间、是否可行。跨多次运行后你还能计算均值和标准差这对写论文非常重要。6. 从结果反推参数为你的VRPTW实例选型合集的最后一步不是跑完四种方法就结束而是要根据一次快速实验的结果反推后续参数策略。我会在用户数据集上先只运行CW和LNS各一次比较两者总距离的相对差距gap (lnsDist - cwDist) / cwDist。这个gap能反映时间窗约束的松紧如果gap小于5%说明CW已经逼近LNS你的实例大概率时间窗很宽松用TS硬约束版加上200次迭代就足够没必要跑LNS如果gap超过15%说明时间窗在强约束直接用LNS并加大迭代次数更划算。下面这段小脚本用来做这个决策% 假设已经得到 cwDist 和 lnsDist gap (lnsDist - cwDist) / cwDist; if gap 0.05 opts_tsh.maxIter 200; final ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsh); elseif gap 0.15 opts_tsp.wLate 1.5; % 加大晚到惩罚更容易跨过不可行区域 final ts_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_tsp); else opts_lns.maxIter 500; opts_lns.removeRatio 0.25; final lns_vrptw(routes_cw, vrp.dist, vrp.demand, vrp.capacity, ... vrp.earliest, vrp.latest, vrp.service, opts_lns); end除了这个规则我还经常用“验证时间窗松紧”的辅助手段计算所有客户时间窗宽度的中位数。中位数越大TS硬约束版越值得信任中位数小于平均服务时间时硬约束版容易找不到合适插入位置这时候惩罚函数版的迭代后期要额外加一个“修复步骤”——把目标函数中的惩罚权重调大让解收敛到可行域。最后一个验证技巧是对同一组参数用不同随机种子跑5次用friedman测试比较四种方法是否有显著性差异。Matlab的friedman函数直接接受一个矩阵行是算法列是实例返回p值。如果p0.05再画multcompare图看哪两个方法差异显著。这样你写在论文里的“对比”才有统计依据而不是只报一个最好平均值。本文还有配套的精品资源点击获取