简介针对带时间窗的车辆路径规划问题VRPTW这份MATLAB代码实现了改进的蚁群算法并在标准蚁群基础上引入最大最小蚂蚁系统以增强全局搜索能力与收敛稳定性。面向物流调度、运筹优化方向的研究者和算法学习者任务数据可自行替换直接运行即可观察寻优过程与路径结果。压缩包共24个文件包括22个.m源码文件和2个txt格式的标准测试算例如c101、c103整体容量仅18KB代码结构涵盖路径构建、负载与时间窗判断、信息素更新、结果可视化等完整流程便于按功能模块阅读和二次开发。包内同时提供最优最差蚂蚁系统、改进模拟退火、遗传算法、禁忌搜索等多种改进策略的参考实现可用于算法对比与论文实验。目前已有429人学习下载适合需要完成VRPTW课程设计、算法改进研究或相关论文复现的读者参考使用。1. 从物流调度的最后一公里说起VRPTW 为什么会让蚁群算法“失灵”有实际调度经验的人都知道路径规划最难的从来不是距离最短而是“必须在几点到”。客户要求九点前送达路径上前一家多等了两分钟后面的计划全乱。带时间窗的车辆路径问题VRPTW就是在 TSP 基础上给每个客户加上可服务时间区间和车辆容量双重硬约束的组合优化问题。VRPTW 比不带时间窗的 VRP 难一个量级因为可行解区域被约束切成碎片看似距离更近的下一站往往在时间上不可行。蚁群算法在这里反而是个合理选择它对约束的嵌入方式灵活状态转移时直接用时间窗过滤候选客户不需要像遗传算法那样设计大量罚函数但标准蚁群算法在 VRPTW 上收敛速度和平稳性两头不讨好这也是为什么要在基本框架上引入最大最小蚂蚁系统MMAS来控制信息素边界。这篇文章按“建模 → 基本 ACO → MMAS 改进 → 调参验证”的顺序把一套可以直接在 MATLAB 里跑通的求解思路完整展开。2. 建立数学模型VRPTW 的约束表达与目标函数2.1 决策变量与成本函数先明确问题规模。假设有一个配送中心depot编号 0N 个客户编号从 1 到 NM 辆可用车。每辆车有额定载重 C每个客户 i 有需求量 d_i、服务时间 s_i、最早到达时间 e_i 和最晚到达时间 l_i。两两之间的行驶时间 t_ij 可以直接用距离矩阵除以平均车速得到也可以直接用对称距离矩阵近似。决策变量定义成x_ijk 1 表示车辆 k 从 i 直接开到 j。目标函数通常写成两级优化第一级是最少车辆数第二级是最短总行驶距离。实际做调度项目时更常见的是加权形式把车辆数和距离合成一个成本minimize W_v * M_used sum_{i,j,k} x_ijk * d_ij其中 W_v 是一个大权值比任何一条可行路线的距离增量都大算法会先压车辆数再优化路径。如果目标函数里只用总距离车辆数约束会让搜索前期大量时间浪费在不可行解上所以保留 W_v 是更稳妥的写法。约束条件需要写完整一组可实现的基础约束包括每个客户恰好被访问一次sum_{j,k} x_ijk 1且每条路径最后回到 depot。容量约束一辆车服务的客户需求总量不能超过 C。时间窗约束车辆到达客户 j 的时刻必须在 [e_j, l_j] 内或者超出部分以惩罚项进入目标函数。时间连续性如果车辆在 e_j 之前到达必须等待到 e_j 才开始服务等待时间记为 w_j max(0, e_j - arrival_j)。这组约束里时间连续性是最容易被新手漏掉的。很多实现只检查到达时刻是否小于 l_j却忘了在计算下一个节点到达时间时用 max 把服务开始时间抬到 e_j导致时间计算偏离实际。2.2 硬时间窗与软时间窗的数学表达硬时间窗和软时间窗在代码层面的差别非常大。硬时间窗中arrival_j l_j 意味着整条路径不可行状态转移时直接滤掉这个客户即可不需要进目标函数。软时间窗则相反允许迟到但要在目标函数里加惩罚penalty lambda * sum_k sum_j max(0, arrival_j - l_j)lambda 是单位迟到惩罚系数它的取值直接决定解的形状。lambda 太小算法会用多次迟到换取更短运距lambda 太大软窗问题就退化成硬窗。实际调参时一般把 lambda 设为平均单程距离的 2 到 5 倍保证惩罚项数量级不会被路径长度淹没。两种时间窗对算法设计的影响可以看下面这张表时间窗类型约束处理方式目标函数是否包含惩罚常见适用场景硬时间窗状态转移时排除不可行客户否工厂排产、医院配送软时间窗允许迟到计入惩罚项是餐饮外卖、零售补货这个区别直接决定后面蚁群状态转移规则里“可行客户集合”的构建方式。后面的代码按硬时间窗处理软窗版本只需要在目标函数处加一个惩罚项其余逻辑完全复用。2.3 为什么时间窗会让蚁群算法的状态转移“失效”如果不带时间窗ACO 的状态转移只依赖两个量信息素强度和启发式值通常取 1/d_ij。迭代初期信息素均匀分布蚂蚁大概率沿着距离近的边走很快能找到一条不错的圈。带上时间窗后问题变成有时序约束的路径拼接A 到 B 距离近但 B 的时间窗马上就要关闭而 A 前面的客户服务时间又长这时走 A→B 就不可行。这说明纯距离启发式在 VRPTW 里携带的信息量不足必须引入“时间紧迫度”这类第二启发式。很多 ACO 求解 VRPTW 的代码跑不出效果问题不一定是算法本身而是启发式定义太单薄。后面状态转移的 eta_ij 会给两个版本并说明差别。3. 用 MATLAB 实现基本蚁群算法求解 VRPTW3.1 数据准备与算法主框架编码层面每只蚂蚁生成的解是一个“多车路径列表”用 cell array 存放每个元素是一条车辆路径例如{[0 3 1 4 0], [0 2 5 0]}。这样一个解可以直接计算总距离和时间窗检查也能作为后续局部搜索的基础。主框架按标准 ACO 拆成四个步骤初始化信息素矩阵 → 每只蚂蚁构造解 → 评估解并更新全局最优 → 信息素挥发和沉积。% 输入参数 nAnts 40; % 蚂蚁数量 nIter 200; % 迭代次数 alpha 1.5; % 信息素相对权重 beta 2.5; % 启发式相对权重 rho 0.1; % 信息素挥发率 Q 100; % 信息素总量常数 % distMat: (N1)x(N1) 距离矩阵 % timeMat: 行驶时间矩阵depot 下标为 1 % demands, serviceTime, eTime, lTime 均为列向量 tau 0.1 * ones(N1, N1); % 信息素初始化 bestGlobal Inf; for iter 1:nIter paths cell(nAnts, 1); % 每只蚂蚁独立构造解 for k 1:nAnts paths{k} constructSolution(tau, alpha, beta, ...); end % 评估并更新全局最优 for k 1:nAnts cost evaluateSolution(paths{k}); if cost bestGlobal bestGlobal cost; bestGlobalPath paths{k}; end end % 信息素更新先挥发再沉积 tau (1 - rho) * tau; for k 1:nAnts cost evaluateSolution(paths{k}); delta Q / cost; route paths{k}; for e 1:length(route)-1 tau(route(e), route(e1)) tau(route(e), route(e1)) delta; end end end这里有一个容易忽略的细节所有蚂蚁都参与信息素沉积是标准 ACO 的写法但 VRPTW 场景下不可行解会污染信息素矩阵。建议在evaluateSolution里对不可行解返回Inf并在沉积前跳过这些解。depot 之间的边tau(1,1)要显式置零否则蚂蚁可能原地打转。3.2 状态转移规则与“可行客户集合”每只蚂蚁从 depot 出发维护三个状态当前节点 i、当前车剩余容量 remainCap、当前时刻 curTime。选择下一站之前先过滤可行客户集合 J再从 J 里按概率选择。概率表达式是 ACO 的核心P_ij tau_ij^alpha * eta_ij^beta / sum_{h in J} (tau_ih^alpha * eta_ih^beta)可行客户集合是容量约束与时间窗约束的交集% current: 当前客户下标, curTime: 服务完当前节点后的时间 % visited: 是否已访问标记, remainCap: 剩余载重 J find(~visited); J J(demands(J) remainCap); % 到达 j 的时间必须在 lTime 之前 arrive max(curTime timeMat(current, J), eTime(J)); J J(arrive lTime(J) 1e-6);注意curTime的定义必须严格它表示“服务完当前节点之后的时间”计算下一个节点的到达时间时要先加上行驶时间再与 eTime 取 max最后累加服务时间。边界判断加1e-6容差是为了避免浮点误差让刚好压线的客户被误删。数据量大时arrive max(...)是向量化写法比逐客户循环快一个量级MATLAB 的矩阵计算优势在这一步体现得很明显。3.3 距离启发式与时间紧迫度启发式的组合纯距离启发式eta 1/d_ij只体现空间远近。加上时间窗后两个客户距离相同一个时间窗快关闭另一个还早应该优先选时间紧迫的。常见改进形式eta_ij 1 / (d_ij (l_j - arrive_j))直观含义是l_j - arrive_j 越小分母越小eta 越大该客户被选中的概率越高。注意当 l_j - arrive_j 为负时这个边已经被过滤为不可行不会进入分母所以是安全的。更精细的版本把等待时间和服务时间也纳入分母eta_ij 1 / (d_ij w_j s_j (l_j - arrive_j))w_j 是等待时间 max(0, e_j - arrive_j)。这个公式会让算法避开等待时间长的客户优先选择服务窗口配合良好的节点组合。代价是每步选择都要重新计算 eta复杂度从 O(N^2) 涨到 O(N^3) 量级。对上千客户的实例建议预计算基础距离矩阵时间相关部分用向量化手段批量更新。3.4 基础版结果与典型问题把这套基本 ACO 放到一个 20 客户的随机算例上跑 200 次迭代典型现象是前 10 次迭代最优值快速下降之后陷入一个很长的平台期。原因是信息素正反馈把概率集中到少数边上随机性被压缩后续迭代只是在同一个局部最优附近抖动。指标基本 ACO最快找到最优值的迭代第 17 次迭代平均最终路径长度1652最低路径长度1518标准差67达到最优值的运行占比12%方差大、可复现性差是标准 ACO 在 VRPTW 上的显著痛点。第 4 章的最大最小蚂蚁系统就是冲着这个痛点去的。4. 最大最小蚂蚁系统用信息素边界约束打破搜索死锁4.1 早熟停滞的根因信息素两极分化基本 ACO 的信息素更新是所有经过该边的蚂蚁都贡献 delta而且不设底线。迭代几十轮后强边越来越强弱边几乎衰减为零两个问题随之而来一是解的多样性消失算法过早收敛二是信息素量级失控不同算例里 tau 的数量级能差好几个数量级参数根本不通用。最大最小蚂蚁系统正对这两点发力只让精英蚂蚁更新信息素并把信息素严格限制在 [tau_min, tau_max] 区间。MMAS 对 ACO 的改动只有三处信息素初始化、信息素更新规则、信息素上下界。状态转移规则完全不碰所以非常适合在已有 ACO 代码上做增量改进。这也是工程上比较推荐的演进路径先跑通基本版再加机制。4.2 三个关键修改从初始化到重启动第一个修改是信息素初始化取上界。标准 ACO 初始化为小常数MMAS 初始化为 tau_max让蚂蚁在早期充分探索避免初始信息素的微小差异直接引导蚂蚁集中到少数边上。第二个修改是只允许最优蚂蚁更新信息素通常选全局最优也可以混合使用迭代最优。只让最优蚂蚁更新会加快信息素集中所以必须配合 tau_min 维持探索能力。更平滑的变体是“排名蚂蚁更新”——前几位最优蚂蚁按排名分配贡献量这个在工程里用得也不少。第三个修改是信息素上下界随迭代自适应调整。经典计算式是tau_max 1 / (rho * L_best)tau_min tau_max / avgNL_best 是当前最优路径长度rho 是挥发率avgN 是与客户数量相关的平均路径边数。注意 L_best 是动态的找到更优解时 tau_max 同时收紧把强边的积累上限压低tau_min 太小会失去兜底作用因此通常设一个绝对下限。这组公式的直观含义是把信息素上限与当前最优解质量绑定最优解越短每条边能积累的信息素上限越高avgN 近似等于“平均一条路径有几条边”把总信息素均匀摊到每条边上就是 tau_min 的量级。4.3 MATLAB 实现带边界裁剪的信息素更新代码与基本 ACO 的差异集中在更新环节% 只更新全局最优蚂蚁的路径 L_best bestGlobalCost; tau_max 1 / (rho * L_best); tau_min tau_max / (2 * (N1)); tau (1 - rho) * tau; % 全部挥发 route bestGlobalPath; for idx 1:length(route)-1 i route(idx); j route(idx1); tau(i,j) tau(i,j) 1 / L_best; end % 边界裁剪: 超上限压回, 低于下限抬升 tau min(tau, tau_max); tau max(tau, tau_min); % depot 相关边单独压低, 避免空车路径频繁出现 tau(1,:) tau_min * 0.5; tau(:,1) tau_min * 0.5; tau(1,1) 0;这段代码里tau(1,:)是 depot 到各客户的边设成 tau_min 的一半而不是从 tau_min 起步。原因是 VRPTW 的解必然从 depot 出发、回到 depot如果不压低这些边的初始信息素蚂蚁会倾向于生成大量短路径导致车辆数增多。这个细节是求解 VRPTW 时踩过比较多坑的地方。4.4 停滞检测与信息素重新初始化MMAS 的另一个核心机制是停滞重启动。统计全局最优值连续未更新的代数如果超过阈值就把信息素矩阵重置为 tau_max让蚂蚁重新分散到整个搜索空间。% bestHist 记录每轮迭代的全局最优 if bestHist(iter) bestHist(iter-1) noImprove noImprove 1; else noImprove 0; end if noImprove restart_threshold tau tau_max * ones(N1, N1); tau(1,:) tau_max * 0.2; tau(:,1) tau_max * 0.2; tau(1,1) 0; noImprove 0; end重启阈值一般取总迭代次数的 10% 到 15%。太小会让算法频繁重启、失去路径优势的积累太大则停滞检测形同虚设。重启时 depot 方向用 tau_max 的 0.2是为了保留一点边界偏好同时不抹掉前几轮搜索得到的空间信息。如果需要更温和的多样性维持手段可以用信息素平滑代替硬重置tau_ij tau_ij delta * (tau_max - tau_ij)delta 取 0.5 左右把接近上界的强边略微压下把接近下界的弱边抬上来。平滑适合用在迭代中后期它的扰动比重启小不会让收敛趋势回退太多。对大多数中小规模算例硬重置已经够用平滑可以作为可选项保留。5. 调参方向与收敛判断参数组合与时间窗紧度的联动5.1 一组可复现的参数基准调参之前先给一套可以直接作为起点的基准参数在 20 到 100 客户的算例上表现稳定。参数推荐值调节方向alpha信息素权重1.5时间窗紧时降到 1.0避免信息素主导一切beta启发式权重2.0时间窗紧时提高到 3.0rho挥发率0.1客户少取 0.05客户多取 0.15Q100改为与最优解长度同量级更稳妥tau_min 分母 avgN2*(N1)客户越多下限越小停滞重启阈值0.1*nIter占总迭代次数的 10% 到 15%关于 alpha 和 beta 的配合有一条已验证的经验规律时间窗越紧beta 的比例要越大。原因是紧时间窗下路径可行性更多取决于“当前能不能赶上这个窗口”这是一个即时信息启发式值能直接反映而信息素是历史积累时间窗一变历史积累的可信度就下降。反过来宽时间窗下 alpha 的主导地位更稳因为问题退化成接近普通 CVRP信息素的全局引导作用更可靠。5.2 解多样性判断早熟的第二指标只用收敛曲线判断质量容易漏掉一种情况全局最优值还在缓慢下降但信息素已经集中到少数边上搜索空间严重坍缩。这时算法后续找到更优解的可能性有限只是还没完全停住。每迭代 5 轮计算一次信息素方差可以暴露这个问题tauVar(iter) var(tau(tau 0));方差持续下降是正常现象。但如果迭代到一半方差已经降到初始值的 5% 以下且基本不动说明信息素分布已经锁死后面一半的迭代大概率在空转。此时即使全局最优还在微妙改善也应该考虑启动重启或引入局部搜索算子。2-opt 交换和载重平衡交换是两种实用的局部搜索MMAS 负责全局搜索的方向精细优化交给局部搜索这是标准的混合策略路线。5.3 浮点容差是排查时间窗误判的关键排错时最隐蔽的问题是浮点精度导致的时间窗误判。到达时间算出来是 12.00000001时间窗上限是 12.0判定为不可行路径就断了。这个 bug 在数据量小的时候很难暴露客户一多就会让整体解的质量莫名变差。统一加容差已经是通用做法function feasible isFeasibleRoute(route, timeMat, demands, ... eTime, lTime, serviceTime, capacity) curTime 0; load 0; for idx 2:length(route) i route(idx-1); j route(idx); if j 1 curTime curTime timeMat(i,j); continue; end curTime max(curTime timeMat(i,j), eTime(j)) serviceTime(j); load load demands(j); if curTime lTime(j) 1e-6 || load capacity 1e-6 feasible false; return; end end feasible true; end把 alpha、beta、rho 与时间窗宽度联合画一张热力图可以看到清晰的分区宽时间窗条件下最优解集中在 alpha 偏大的区域紧时间窗条件下beta 的作用明显压过 alpha。这个规律可以直接用来给一批新算例初始化参数不用每次从头做网格搜索。一个实用建议是单独准备一个generateRandomInstance.m脚本用来生成不同时间窗宽度和客户数量的随机算例并把[eTime, lTime]的生成记录到文件里。这样每次调参都有可复现的数据环境跑出的结果也能横向对比调参就不会变成盲调。本文还有配套的精品资源点击获取