简介面向物流调度与优化算法学习者这份MATLAB源码包实现了粒子群算法PSO求解带时间窗的多客户单仓库车辆路径规划问题TWVRP。通过模拟鸟群觅食行为在距离最短与时间窗约束之间求得可行路径既适合作为智能优化算法的入门实践也可为实际配送调度提供参考。压缩包共12个文件包含5个m脚本算法主体与核心函数、5张运行结果截图、1份README说明和1个数据文件整体仅229KB便于快速下载与复现。目前已有421人学习浏览。读者可以从中了解PSO求解TWVRP的完整流程包括粒子初始化、适应度函数设计、速度位置更新以及时间窗约束处理等关键环节获得一个代码与结果并重、可直接运行的优化算法实验案例。1. 带时间窗的多客户单仓库问题为什么先看粒子群而不是扫描法做配送调度时最头疼的不是距离长短而是每个客户都有硬性的最早到达时间和最晚到达时间。带时间窗的多客户单仓库车辆路径规划问题TWVRP用粒子群算法求解时核心是把一条合法路线当作一个粒子通过速度和位置更新在解空间里翻找更优顺序。这类问题组合爆炸比普通VRP更猛扫描法、节约法这类确定性构造法在时间窗约束下很容易陷入局部选择。这套Matlab源码把 main.m、getCost.m、isValid.m、updateVelocity.m、adaMatrix.m 拆开覆盖解编码、成本计算、约束校验、速度更新和自适应参数调整五个环节适合用来跑课程设计、对比算法也适合第一次接触PSO的人拆开看粒子是怎么在离散路线问题上起作用的。2. 从 main.m 到 getCost.mTWVRP 的解编码与目标函数2.1 连续粒子位置如何对应一条离散路线标准粒子群算法的位置和速度都是连续实数向量它善于在连续域里做搜索。TWVRP 的候选解却是一组客户访问顺序属于离散排列。要让 PSO 跑起来必须建立一套“连续实数向量 → 客户排列”的映射。最常见的做法是实数排序编码粒子位置向量长度等于客户数第 i 维的数值大小决定客户 i 在访问顺序中的先后位置。数值最小的维度对应的客户排在最前数值最大的排在最后。例如一个有 5 个客户的粒子位置可能是x [0.81, 0.20, 0.55, 0.44, 0.72]按数值排序后得到顺序[2, 4, 3, 5, 1]这就是一条访问路线。这种编码的优点是位置和速度更新公式完全不需要改动解码时只做一次sort群体迭代的连续性不会被破坏。缺点是两个几何距离很远的位置向量可能解码出相近的路线反过来很小的位置扰动也可能让顺序发生大幅变化因此后期需要配合局部搜索。PSO 术语TWVRP 中的含义粒子位置 x一组客户顺序的连续编码粒子速度 v决策变量的扰动量和方向pbest当前轨迹上成本最低的路线gbest群体到目前为止发现的最优路线适应度函数getCost 返回的总成本2.2 多车辆时分清每辆车负责哪一段单仓库多客户意味着多辆车协同配送只有客户顺序还不够还要决定在哪个位置换下一辆车。源码中常见的处理方式是在位置向量里追加分隔维度维度数量为车辆数减一。解码时先对所有维度排序得到客户顺序后再用这些分隔维度切分路线。% 位置向量 x 前 nCustomer 维对应客户后 nVehicle-1 维对应分隔点 % 返回值 routes 是元胞数组routes{k} 表示第 k 辆车的客户序列 function routes decodeRoutes(x, nCustomer, nVehicle) seq zeros(1, nCustomer); % 客户顺序占位 [~, ord] sort(x(1:nCustomer)); seq ord; sep x(nCustomer1:end); [~, sepIdx] sort(sep); % 取分隔点在 seq 中的下标将 seq 切成多段 cutPos round(sepIdx / nVehicle * (nCustomer - 1)) 1; routes cell(1, nVehicle); prev 0; for k 1:nVehicle routes{k} seq(prev1:min(cutPos(k), nCustomer)); prev min(cutPos(k), nCustomer); end end这段代码里seq sort(x(1:nCustomer))得到客户访问编号序列sep是分隔点。实际工程中我更常把车辆数固定为已知值用容量约束去切而不是用分隔点因为后者容易出现空车或车辆数不匹配。cutPos计算可以改成直接cumsum按比例切分但无论哪种切法都要保证每辆车至少能服务一个客户否则 getCost 里会多出不必要的空驶成本。2.3 getCost 把距离、等待和超窗惩罚统一成一个可比较的数字getCost.m 是粒子群算法的适应度函数它的结果决定了 pbest 和 gbest 的更新。目标函数不能只写行驶总距离因为时间窗会让某些路段在物理上不可行。常见做法是构造一个加法模型总成本 行驶距离 等待时间权重 超窗惩罚。function totalCost getCost(x, data) routes decodeRoutes(x, data.nCustomer, data.nVehicle); distCost 0; waitCost 0; overCost 0; for k 1:length(routes) route [0, routes{k}, 0]; % 0 表示仓库便于统一计算 now data.leaveTime(k); % 每辆车有自己的出发时刻 for i 1:length(route)-1 from route(i) 1; % MATLAB 索引从 1 开始 to route(i1) 1; travel data.travelTime(from, to); now now travel; if route(i1) ~ 0 if now data.ready(route(i1)) waitCost waitCost (data.ready(route(i1)) - now); now data.ready(route(i1)); % 早到等待 end if now data.due(route(i1)) overCost overCost (now - data.due(route(i1))); end now now data.service(route(i1)); end end distCost distCost data.dist(from, 1); % 回仓库 end totalCost distCost 0.3 * waitCost 100 * overCost; endwaitCost是早到后的等待时长overCost是晚到总时长权重 0.3 和 100 不是固定值。这个系数就是源码里最值得调的地方如果超窗惩罚系数设得太大粒子群会优先保证时间窗完全不被破坏路线绕远也可能被接受设得太小最优解可能大面积超窗isValid 阶段又全部被过滤掉。我一般先把 100 设在比最大行驶距离还大一个数量级的位置再根据初次运行结果按 2 倍步长缩。注意waitCost和overCost的权重不是死的它们决定了算法最后输出的是“严格守时”还是“总成本最小”的路线。3. isValid 校验与时间窗约束硬性过滤还是惩罚项3.1 时间窗建模里的三个时间量TWVRP 的约束本质是每个客户都有一段可服务区间 [e_i, l_i]。车辆到达时间早于 e_i 时只能等待到达时间晚于 l_i 时判定超窗。除了这两个量时间窗问题还包含服务时间 s_i 和车辆在客户间的行驶时间 t_ij。三个时间量必须用同一把时间尺子仓库出发时间、路况折算、装卸时间都要加起来。常见错误是把行驶距离当成行驶时间使用。数据文件 data.txt 里面存的通常只有客户坐标和需求量因此源码里要按坐标算欧氏距离然后除以平均车速得到时间。如果直接用距离做时间窗判断早到晚到的误差会随着客户距离增大而放大。我们要在初始化 data 结构体时就预先算好data.travelTime矩阵而不是每次 getCost 都重复计算坐标距离。% 初始化时从坐标矩阵 xy 构造时间矩阵 n size(xy, 1); speed 50; % 假设平均车速 50 km/h travelTime zeros(n, n); for i 1:n for j 1:n travelTime(i, j) norm(xy(i, :) - xy(j, :)) / speed; end end上面代码的speed是全局参数50 km/h 只是初始值。实际项目中如果路网有单行线或者拥堵时段这个矩阵应该替换成外部 API 返回的耗时矩阵。对课程设计来说欧氏距离除以常数速度是最稳的方案既能保证逻辑正确也不会引入不必要的外部依赖。3.2 isValid 的递推判断从仓库出发到返回仓库isValid.m 的作用是判定一条路线是否完全满足时间窗。判断逻辑并不复杂从仓库出发时刻开始依次累加每段行驶时间、等待时间和服务时间每次到达客户点都检查当前时间是否落在 [e_i, l_i] 内。由于服务完成后车辆才会去下一个点所以时间轴只能沿路线前进方向推进。function ok isValid(route, e, l, s, travel, startTime) now startTime; % 仓库发车时刻 prev 1; % 仓库索引假设为 1 for i 1:length(route) now now travel(prev, route(i)); if now e(route(i)) now e(route(i)); % 早到则等待客户准备就绪 elseif now l(route(i)) ok false; % 晚到直接判为不可行 return; end now now s(route(i)); % 服务时间 prev route(i); end now now travel(prev, 1); % 返回仓库 ok true; end这段代码的核心在if now e(route(i))这一分支。早到时让当前时间跳到客户最早可服务时刻相当于车辆白白等待晚到则直接返回 false。注意仓库本身也有服务窗口源码里通常默认仓库窗口覆盖整个配送时段因此返回仓库不用检查超窗。实际项目中仓库也有最晚回场时间此时需要在函数末尾再加一次now l(1)判断。3.3 在 PSO 里用 isValid 还是用惩罚项两种做法边界isValid.m 在源码里并不是单独出现的它和 getCost.m 存在配合关系。如果每个粒子解码出来的路线都要先过 isValid再用 getCost 计算成本那么不合规的粒子会被直接丢弃这种策略叫做硬过滤。硬过滤的优点是最终输出的路线一定满足时间窗缺点是在迭代前期粒子群很难找到足够多的可行解群体容易陷入局部最优。另一种做法是只在 getCost 里对超窗量施加惩罚不调用 isValid迭代结束后再统一滤波。这种做法搜索空间更大群体不会被早期有限的可行解限制住但需要精心调整惩罚系数。源码目录里同时存在 isValid.m 和 getCost.m我倾向于把 isValid 放在最终输出前做确认而把 getCost 中的超窗惩罚作为引导搜索的主要信号。这样 pbest 和 gbest 可能在途中经历短暂超窗最终收敛的路线却一定是合法的。4. updateVelocity 与 adaMatrix粒子群迭代参数怎么联动4.1 标准速度更新公式和它的离散化陷阱粒子群算法原理其实只有两行更新式速度v w*v c1*r1*(pbest-x) c2*r2*(gbest-x)位置x x v。在 TWVRP 这种离散排列问题上位置向量本身是连续编码所以速度更新可以原样保留但解码时不能用四舍五入或者取整操作直接得到路线顺序而是要把位置向量在连续空间里移动然后统一按数值大小排序。function vNew updateVelocity(v, x, pbest, gbest, w, c1, c2) r1 rand(size(v)); r2 rand(size(v)); vNew w .* v c1 .* r1 .* (pbest - x) c2 .* r2 .* (gbest - x); vMax 0.15 * (max(x) - min(x)); vNew max(min(vNew, vMax), -vMax); endvNew限制在最大速度vMax内避免粒子一步跨过整个搜索区间。max(x)-min(x)在这里起到自适应界限的作用因为每次解码后位置向量的值域会随迭代变化。需要特别注意的是pbest - x和gbest - x得到的是向量差这里的向量维度是客户数加车辆分隔数不能和客户业务编号混用。4.2 adaMatrix 的自适应权重设计adaMatrix.m 这个名字直观理解是自适应矩阵。在多数 PSO 实现里惯性权重 w 是标量但更细的做法是让 w 随迭代次数和粒子状态变化。adaMatrix 的典型设计是迭代初期给较大权重保持全局搜索迭代后期逐渐变小侧重局部精化同时对适应度排名靠前的粒子使用更小权重让精英粒子在最优解附近精细搜索对排名靠后的粒子保留较大权重维持种群多样性。function wMat adaMatrix(x, costs, iter, maxIter) [~, order] sort(costs); % 适应度升序排列 nPop length(costs); wBase 0.9 - 0.5 * iter / maxIter; for i 1:nPop rank find(order i); % 当前粒子排名 wMat(i) wBase 0.1 * (rank / nPop - 0.5); end end上面这个矩阵按粒子维度展开每个粒子在下一代更新时使用不同的惯性权重。order的排名越靠前rank / nPop - 0.5越接近负值权重就越小。这样设计的好处是群体不会在迭代后期被 gbest 吸成一个密集的点因为最差粒子仍有较大的速度惯性可以跳出当前聚集区域再寻找新解。代价是每代都要做一次sort当种群规模达到 1000 以上时这部分开销会明显增加。4.3 参数组合与位置越界的处理把 updateVelocity 和 adaMatrix 放进 main.m 后还需要处理一组关键参数粒子数、迭代次数、学习因子、速度上限。这里给一张我在 TWVRP 这种规模下的默认参数表后续调参可以在这些数值的基础上按 1.5 倍或 0.7 倍缩放。参数推荐范围对 TWVRP 的影响粒子数 nPop60200太小容易早熟太大每代 getCost 计算量线性增加迭代次数 maxIter150400超过 300 代后收益递减可观察收敛曲线确认c1 / c21.52.5c2 偏大会让群体过快集中到 gbestw 初值0.850.95越大越利于前期探索离散排列空间w 终值0.20.4越小越利于后期精细逼近最优路线vMax0.10.25 倍位置范围过大导致频繁改变客户顺序过小导致收敛极慢位置越界处理有两种常见方案。一是直接把越界值裁剪回边界这种方案简单但在边界附近容易堆积粒子。二是把越界部分映射回区间内部比如x(xub) 2*ub - x(xub)让粒子穿过边界后反弹回搜索空间。源码中用哪种取决于 main.m 的赋值方式但从运行结果来看单纯的裁剪在客户数 50 以内表现稳定不需要过度设计。5. data.txt 到运行结果跑通 1407 期源码的完整流程5.1 data.txt 的字段组织与读取拿到源码后第一个动作不是运行 main.m而是先看 data.txt 的字段。常见的 TWVRP 数据文件格式是第一行放两个整数分别是客户数和车辆数从第二行开始每行依次保存客户编号、横坐标、纵坐标、需求量、最早服务时间、最晚服务时间、服务时间。坐标和服务时间单位保持一致需求量可以是件数也可以是重量关键是设好每辆车的容量上限。function data readTWVRP(filename) raw load(filename); data.nCustomer raw(1, 1); data.nVehicle raw(1, 2); info raw(2:end, :); data.xy info(:, 2:3); data.demand info(:, 4); data.ready info(:, 5); % 最早到达时间 data.due info(:, 6); % 最晚到达时间 data.service info(:, 7); % 服务时间 % 预计算距离与行驶时间矩阵供 getCost 和 isValid 复用 ... end这里load可以读入纯数字文本比fscanf更省事。需要留意raw中可能包含表头如果第一行有字符串标签load会直接报错此时要改成readmatrix或者importdata。文件里如果车辆数那一列后面还跟了容量、车场位置等字段也要同步调整raw(1,:)的取值。5.2 main.m 的主循环与停止条件main.m 是整套源码的入口。它按顺序完成五件事读取数据、初始化粒子群的群体矩阵、计算初始适应度、进入迭代更新 pbest 和 gbest、输出结果图片。下面是一个简化的主循环骨架和源码的结构一一对应。data readTWVRP(data.txt); nPop 100; maxIter 200; dim data.nCustomer data.nVehicle - 1; x rand(nPop, dim); % 随机初始化位置 v 0.1 * randn(nPop, dim); % 初始化速度 cost zeros(nPop, 1); pbestCost inf(nPop, 1); gbestCost inf; for iter 1:maxIter for i 1:nPop cost(i) getCost(x(i, :), data); end [~, gIdx] min(cost); for i 1:nPop if cost(i) pbestCost(i) pbestCost(i) cost(i); pbest(i, :) x(i, :); end end if cost(gIdx) gbestCost gbestCost cost(gIdx); gbest x(gIdx, :); end w adaMatrix(x, cost, iter, maxIter); for i 1:nPop v(i, :) updateVelocity(v(i, :), x(i, :), pbest(i, :), gbest, w(i), 2, 2); x(i, :) x(i, :) v(i, :); end end这个骨架里最容易被忽略的是pbestCost inf(nPop,1)和gbestCost inf的初始化。如果初始化为 0那么第一轮迭代的 pbest 和 gbest 永远无法被更新算法会停留在一个糟糕的粒子上。dim的设置必须与 decodeRoutes 的维度保持一致否则sort后切分时会出现数组越界。停止条件除了maxIter还可以加入“连续 20 代 gbest 不变则中断”这个条件很适合 TWVRP因为后期路线顺序很难发生明显推进。5.3 从运行结果图反查算法状态运行后的 5 张结果图通常包含路径图、迭代收敛曲线和车辆时间轴。路径图里每条线代表一辆车的行驶轨迹如果多条线路交叉明显且绕圈严重说明粒子群没有充分搜索排列空间需要增大粒子数或初始权重。收敛曲线要是下降过程突兀说明粒子群在某代突然跳变多半是 gbest 更新的判断写错了把最差解当成最优解。时间窗是否被满足可以从车辆时间轴看出来横轴是时刻纵轴是车辆编号每个客户点显示服务开始时刻。如果某个点的柱形超出了客户窗口区域说明目标函数的超窗惩罚权重偏低。调整时先固定其他参数只动惩罚系数观察整体超窗时长是否随系数上升而下降通常 3 次对比就能确定合适区间。6. 把 PSO-TWVRP 改到真实场景收敛验证与局部搜索6.1 用最近邻种子粒子替代全随机初始化PSO 对初始群体非常敏感。全随机初始化会让种群花掉几十代才凑出几条可行路线尤其时间窗窄的时候更明显。我惯用的做法是在初始化时生成 10% 的种子粒子用最近邻构造客户顺序。比如对仓库随机选一个未服务客户开始每次选择距离当前点最近且满足时间窗的客户加入路线直到车辆容量用完再开下一辆车。这样得到的种子粒子已经是一条合法路线能极大拉低第一代 gbestCost。% 最近邻构造种子路线其余粒子仍随机生成 seed nearestNeighborRoute(data); x(1, :) encodeRoute(seed, data.nCustomer); for i 2:nPop x(i, :) rand(1, dim); endencodeRoute把离散路线反向映射成连续位置向量最笨的方法是按客户出现顺序倒推出一个严格递增的数值序列再穿插分隔维度。对大多数规模在 30 个客户以内的 TWVRP10% 种子粒子不会引起早熟却能让收敛曲线从第一代就开始下降。6.2 用 2-opt 在 gbest 上做局部精化PSO 的全局搜索强但局部精化弱。gbest 可能已经很接近最优路线却因为一次 swap 操作能省掉几十公里而未被发现。常见做法是在每次迭代结束后对 gbest 解码出来的每一条子路线执行 2-opt反转子路线中的任一段如果总成本下降且时间窗合规则接受。function newRoute twoOpt(route, data) improved true; while improved improved false; for i 1:length(route)-1 for j i1:length(route) cand route; cand(i:j) cand(j:-1:i); if isValid(cand, data) routeCost(cand, data) routeCost(route, data) route cand; improved true; end end end end newRoute route; end这段 2-opt 的复杂度是 O(L^2)L 是单辆车客户数放在主循环里每代做一次会增加不少时间。实际使用时可以每隔 5 代才精化一次或者只在输出最终 gbest 前做一轮效果接近且时间开销小得多。6.3 验证收敛是否可靠固定随机种子做的三件事改完参数别急着验收先固定随机种子然后做三件事跑三遍同样参数对比 gbest 和平均成本把客户数减半做一次穷举或动态规划验证最优值的量级把时间窗全部放宽成 0 到当天结束检验 PSO 能否退化成普通 VRP 并得到接近最近邻的结果。这三件事能把参数误设和索引越界问题暴露出来。最后检查 getCost 里惩罚系数确保超窗成本在数值上不会大到掩盖距离项导致算法只保时间窗不保路径质量。本文还有配套的精品资源点击获取 
长河君 · 2026/9/10 11:34:40 · 阅读 1.2k 粒子群算法求解带时间窗车辆路径规划：Matlab源码解析与参数调优
简介面向物流调度与优化算法学习者这份MATLAB源码包实现了粒子群算法PSO求解带时间窗的多客户单仓库车辆路径规划问题TWVRP。通过模拟鸟群觅食行为在距离最短与时间窗约束之间求得可行路径既适合作为智能优化算法的入门实践也可为实际配送调度提供参考。压缩包共12个文件包含5个m脚本算法主体与核心函数、5张运行结果截图、1份README说明和1个数据文件整体仅229KB便于快速下载与复现。目前已有421人学习浏览。读者可以从中了解PSO求解TWVRP的完整流程包括粒子初始化、适应度函数设计、速度位置更新以及时间窗约束处理等关键环节获得一个代码与结果并重、可直接运行的优化算法实验案例。1. 带时间窗的多客户单仓库问题为什么先看粒子群而不是扫描法做配送调度时最头疼的不是距离长短而是每个客户都有硬性的最早到达时间和最晚到达时间。带时间窗的多客户单仓库车辆路径规划问题TWVRP用粒子群算法求解时核心是把一条合法路线当作一个粒子通过速度和位置更新在解空间里翻找更优顺序。这类问题组合爆炸比普通VRP更猛扫描法、节约法这类确定性构造法在时间窗约束下很容易陷入局部选择。这套Matlab源码把 main.m、getCost.m、isValid.m、updateVelocity.m、adaMatrix.m 拆开覆盖解编码、成本计算、约束校验、速度更新和自适应参数调整五个环节适合用来跑课程设计、对比算法也适合第一次接触PSO的人拆开看粒子是怎么在离散路线问题上起作用的。2. 从 main.m 到 getCost.mTWVRP 的解编码与目标函数2.1 连续粒子位置如何对应一条离散路线标准粒子群算法的位置和速度都是连续实数向量它善于在连续域里做搜索。TWVRP 的候选解却是一组客户访问顺序属于离散排列。要让 PSO 跑起来必须建立一套“连续实数向量 → 客户排列”的映射。最常见的做法是实数排序编码粒子位置向量长度等于客户数第 i 维的数值大小决定客户 i 在访问顺序中的先后位置。数值最小的维度对应的客户排在最前数值最大的排在最后。例如一个有 5 个客户的粒子位置可能是x [0.81, 0.20, 0.55, 0.44, 0.72]按数值排序后得到顺序[2, 4, 3, 5, 1]这就是一条访问路线。这种编码的优点是位置和速度更新公式完全不需要改动解码时只做一次sort群体迭代的连续性不会被破坏。缺点是两个几何距离很远的位置向量可能解码出相近的路线反过来很小的位置扰动也可能让顺序发生大幅变化因此后期需要配合局部搜索。PSO 术语TWVRP 中的含义粒子位置 x一组客户顺序的连续编码粒子速度 v决策变量的扰动量和方向pbest当前轨迹上成本最低的路线gbest群体到目前为止发现的最优路线适应度函数getCost 返回的总成本2.2 多车辆时分清每辆车负责哪一段单仓库多客户意味着多辆车协同配送只有客户顺序还不够还要决定在哪个位置换下一辆车。源码中常见的处理方式是在位置向量里追加分隔维度维度数量为车辆数减一。解码时先对所有维度排序得到客户顺序后再用这些分隔维度切分路线。% 位置向量 x 前 nCustomer 维对应客户后 nVehicle-1 维对应分隔点 % 返回值 routes 是元胞数组routes{k} 表示第 k 辆车的客户序列 function routes decodeRoutes(x, nCustomer, nVehicle) seq zeros(1, nCustomer); % 客户顺序占位 [~, ord] sort(x(1:nCustomer)); seq ord; sep x(nCustomer1:end); [~, sepIdx] sort(sep); % 取分隔点在 seq 中的下标将 seq 切成多段 cutPos round(sepIdx / nVehicle * (nCustomer - 1)) 1; routes cell(1, nVehicle); prev 0; for k 1:nVehicle routes{k} seq(prev1:min(cutPos(k), nCustomer)); prev min(cutPos(k), nCustomer); end end这段代码里seq sort(x(1:nCustomer))得到客户访问编号序列sep是分隔点。实际工程中我更常把车辆数固定为已知值用容量约束去切而不是用分隔点因为后者容易出现空车或车辆数不匹配。cutPos计算可以改成直接cumsum按比例切分但无论哪种切法都要保证每辆车至少能服务一个客户否则 getCost 里会多出不必要的空驶成本。2.3 getCost 把距离、等待和超窗惩罚统一成一个可比较的数字getCost.m 是粒子群算法的适应度函数它的结果决定了 pbest 和 gbest 的更新。目标函数不能只写行驶总距离因为时间窗会让某些路段在物理上不可行。常见做法是构造一个加法模型总成本 行驶距离 等待时间权重 超窗惩罚。function totalCost getCost(x, data) routes decodeRoutes(x, data.nCustomer, data.nVehicle); distCost 0; waitCost 0; overCost 0; for k 1:length(routes) route [0, routes{k}, 0]; % 0 表示仓库便于统一计算 now data.leaveTime(k); % 每辆车有自己的出发时刻 for i 1:length(route)-1 from route(i) 1; % MATLAB 索引从 1 开始 to route(i1) 1; travel data.travelTime(from, to); now now travel; if route(i1) ~ 0 if now data.ready(route(i1)) waitCost waitCost (data.ready(route(i1)) - now); now data.ready(route(i1)); % 早到等待 end if now data.due(route(i1)) overCost overCost (now - data.due(route(i1))); end now now data.service(route(i1)); end end distCost distCost data.dist(from, 1); % 回仓库 end totalCost distCost 0.3 * waitCost 100 * overCost; endwaitCost是早到后的等待时长overCost是晚到总时长权重 0.3 和 100 不是固定值。这个系数就是源码里最值得调的地方如果超窗惩罚系数设得太大粒子群会优先保证时间窗完全不被破坏路线绕远也可能被接受设得太小最优解可能大面积超窗isValid 阶段又全部被过滤掉。我一般先把 100 设在比最大行驶距离还大一个数量级的位置再根据初次运行结果按 2 倍步长缩。注意waitCost和overCost的权重不是死的它们决定了算法最后输出的是“严格守时”还是“总成本最小”的路线。3. isValid 校验与时间窗约束硬性过滤还是惩罚项3.1 时间窗建模里的三个时间量TWVRP 的约束本质是每个客户都有一段可服务区间 [e_i, l_i]。车辆到达时间早于 e_i 时只能等待到达时间晚于 l_i 时判定超窗。除了这两个量时间窗问题还包含服务时间 s_i 和车辆在客户间的行驶时间 t_ij。三个时间量必须用同一把时间尺子仓库出发时间、路况折算、装卸时间都要加起来。常见错误是把行驶距离当成行驶时间使用。数据文件 data.txt 里面存的通常只有客户坐标和需求量因此源码里要按坐标算欧氏距离然后除以平均车速得到时间。如果直接用距离做时间窗判断早到晚到的误差会随着客户距离增大而放大。我们要在初始化 data 结构体时就预先算好data.travelTime矩阵而不是每次 getCost 都重复计算坐标距离。% 初始化时从坐标矩阵 xy 构造时间矩阵 n size(xy, 1); speed 50; % 假设平均车速 50 km/h travelTime zeros(n, n); for i 1:n for j 1:n travelTime(i, j) norm(xy(i, :) - xy(j, :)) / speed; end end上面代码的speed是全局参数50 km/h 只是初始值。实际项目中如果路网有单行线或者拥堵时段这个矩阵应该替换成外部 API 返回的耗时矩阵。对课程设计来说欧氏距离除以常数速度是最稳的方案既能保证逻辑正确也不会引入不必要的外部依赖。3.2 isValid 的递推判断从仓库出发到返回仓库isValid.m 的作用是判定一条路线是否完全满足时间窗。判断逻辑并不复杂从仓库出发时刻开始依次累加每段行驶时间、等待时间和服务时间每次到达客户点都检查当前时间是否落在 [e_i, l_i] 内。由于服务完成后车辆才会去下一个点所以时间轴只能沿路线前进方向推进。function ok isValid(route, e, l, s, travel, startTime) now startTime; % 仓库发车时刻 prev 1; % 仓库索引假设为 1 for i 1:length(route) now now travel(prev, route(i)); if now e(route(i)) now e(route(i)); % 早到则等待客户准备就绪 elseif now l(route(i)) ok false; % 晚到直接判为不可行 return; end now now s(route(i)); % 服务时间 prev route(i); end now now travel(prev, 1); % 返回仓库 ok true; end这段代码的核心在if now e(route(i))这一分支。早到时让当前时间跳到客户最早可服务时刻相当于车辆白白等待晚到则直接返回 false。注意仓库本身也有服务窗口源码里通常默认仓库窗口覆盖整个配送时段因此返回仓库不用检查超窗。实际项目中仓库也有最晚回场时间此时需要在函数末尾再加一次now l(1)判断。3.3 在 PSO 里用 isValid 还是用惩罚项两种做法边界isValid.m 在源码里并不是单独出现的它和 getCost.m 存在配合关系。如果每个粒子解码出来的路线都要先过 isValid再用 getCost 计算成本那么不合规的粒子会被直接丢弃这种策略叫做硬过滤。硬过滤的优点是最终输出的路线一定满足时间窗缺点是在迭代前期粒子群很难找到足够多的可行解群体容易陷入局部最优。另一种做法是只在 getCost 里对超窗量施加惩罚不调用 isValid迭代结束后再统一滤波。这种做法搜索空间更大群体不会被早期有限的可行解限制住但需要精心调整惩罚系数。源码目录里同时存在 isValid.m 和 getCost.m我倾向于把 isValid 放在最终输出前做确认而把 getCost 中的超窗惩罚作为引导搜索的主要信号。这样 pbest 和 gbest 可能在途中经历短暂超窗最终收敛的路线却一定是合法的。4. updateVelocity 与 adaMatrix粒子群迭代参数怎么联动4.1 标准速度更新公式和它的离散化陷阱粒子群算法原理其实只有两行更新式速度v w*v c1*r1*(pbest-x) c2*r2*(gbest-x)位置x x v。在 TWVRP 这种离散排列问题上位置向量本身是连续编码所以速度更新可以原样保留但解码时不能用四舍五入或者取整操作直接得到路线顺序而是要把位置向量在连续空间里移动然后统一按数值大小排序。function vNew updateVelocity(v, x, pbest, gbest, w, c1, c2) r1 rand(size(v)); r2 rand(size(v)); vNew w .* v c1 .* r1 .* (pbest - x) c2 .* r2 .* (gbest - x); vMax 0.15 * (max(x) - min(x)); vNew max(min(vNew, vMax), -vMax); endvNew限制在最大速度vMax内避免粒子一步跨过整个搜索区间。max(x)-min(x)在这里起到自适应界限的作用因为每次解码后位置向量的值域会随迭代变化。需要特别注意的是pbest - x和gbest - x得到的是向量差这里的向量维度是客户数加车辆分隔数不能和客户业务编号混用。4.2 adaMatrix 的自适应权重设计adaMatrix.m 这个名字直观理解是自适应矩阵。在多数 PSO 实现里惯性权重 w 是标量但更细的做法是让 w 随迭代次数和粒子状态变化。adaMatrix 的典型设计是迭代初期给较大权重保持全局搜索迭代后期逐渐变小侧重局部精化同时对适应度排名靠前的粒子使用更小权重让精英粒子在最优解附近精细搜索对排名靠后的粒子保留较大权重维持种群多样性。function wMat adaMatrix(x, costs, iter, maxIter) [~, order] sort(costs); % 适应度升序排列 nPop length(costs); wBase 0.9 - 0.5 * iter / maxIter; for i 1:nPop rank find(order i); % 当前粒子排名 wMat(i) wBase 0.1 * (rank / nPop - 0.5); end end上面这个矩阵按粒子维度展开每个粒子在下一代更新时使用不同的惯性权重。order的排名越靠前rank / nPop - 0.5越接近负值权重就越小。这样设计的好处是群体不会在迭代后期被 gbest 吸成一个密集的点因为最差粒子仍有较大的速度惯性可以跳出当前聚集区域再寻找新解。代价是每代都要做一次sort当种群规模达到 1000 以上时这部分开销会明显增加。4.3 参数组合与位置越界的处理把 updateVelocity 和 adaMatrix 放进 main.m 后还需要处理一组关键参数粒子数、迭代次数、学习因子、速度上限。这里给一张我在 TWVRP 这种规模下的默认参数表后续调参可以在这些数值的基础上按 1.5 倍或 0.7 倍缩放。参数推荐范围对 TWVRP 的影响粒子数 nPop60200太小容易早熟太大每代 getCost 计算量线性增加迭代次数 maxIter150400超过 300 代后收益递减可观察收敛曲线确认c1 / c21.52.5c2 偏大会让群体过快集中到 gbestw 初值0.850.95越大越利于前期探索离散排列空间w 终值0.20.4越小越利于后期精细逼近最优路线vMax0.10.25 倍位置范围过大导致频繁改变客户顺序过小导致收敛极慢位置越界处理有两种常见方案。一是直接把越界值裁剪回边界这种方案简单但在边界附近容易堆积粒子。二是把越界部分映射回区间内部比如x(xub) 2*ub - x(xub)让粒子穿过边界后反弹回搜索空间。源码中用哪种取决于 main.m 的赋值方式但从运行结果来看单纯的裁剪在客户数 50 以内表现稳定不需要过度设计。5. data.txt 到运行结果跑通 1407 期源码的完整流程5.1 data.txt 的字段组织与读取拿到源码后第一个动作不是运行 main.m而是先看 data.txt 的字段。常见的 TWVRP 数据文件格式是第一行放两个整数分别是客户数和车辆数从第二行开始每行依次保存客户编号、横坐标、纵坐标、需求量、最早服务时间、最晚服务时间、服务时间。坐标和服务时间单位保持一致需求量可以是件数也可以是重量关键是设好每辆车的容量上限。function data readTWVRP(filename) raw load(filename); data.nCustomer raw(1, 1); data.nVehicle raw(1, 2); info raw(2:end, :); data.xy info(:, 2:3); data.demand info(:, 4); data.ready info(:, 5); % 最早到达时间 data.due info(:, 6); % 最晚到达时间 data.service info(:, 7); % 服务时间 % 预计算距离与行驶时间矩阵供 getCost 和 isValid 复用 ... end这里load可以读入纯数字文本比fscanf更省事。需要留意raw中可能包含表头如果第一行有字符串标签load会直接报错此时要改成readmatrix或者importdata。文件里如果车辆数那一列后面还跟了容量、车场位置等字段也要同步调整raw(1,:)的取值。5.2 main.m 的主循环与停止条件main.m 是整套源码的入口。它按顺序完成五件事读取数据、初始化粒子群的群体矩阵、计算初始适应度、进入迭代更新 pbest 和 gbest、输出结果图片。下面是一个简化的主循环骨架和源码的结构一一对应。data readTWVRP(data.txt); nPop 100; maxIter 200; dim data.nCustomer data.nVehicle - 1; x rand(nPop, dim); % 随机初始化位置 v 0.1 * randn(nPop, dim); % 初始化速度 cost zeros(nPop, 1); pbestCost inf(nPop, 1); gbestCost inf; for iter 1:maxIter for i 1:nPop cost(i) getCost(x(i, :), data); end [~, gIdx] min(cost); for i 1:nPop if cost(i) pbestCost(i) pbestCost(i) cost(i); pbest(i, :) x(i, :); end end if cost(gIdx) gbestCost gbestCost cost(gIdx); gbest x(gIdx, :); end w adaMatrix(x, cost, iter, maxIter); for i 1:nPop v(i, :) updateVelocity(v(i, :), x(i, :), pbest(i, :), gbest, w(i), 2, 2); x(i, :) x(i, :) v(i, :); end end这个骨架里最容易被忽略的是pbestCost inf(nPop,1)和gbestCost inf的初始化。如果初始化为 0那么第一轮迭代的 pbest 和 gbest 永远无法被更新算法会停留在一个糟糕的粒子上。dim的设置必须与 decodeRoutes 的维度保持一致否则sort后切分时会出现数组越界。停止条件除了maxIter还可以加入“连续 20 代 gbest 不变则中断”这个条件很适合 TWVRP因为后期路线顺序很难发生明显推进。5.3 从运行结果图反查算法状态运行后的 5 张结果图通常包含路径图、迭代收敛曲线和车辆时间轴。路径图里每条线代表一辆车的行驶轨迹如果多条线路交叉明显且绕圈严重说明粒子群没有充分搜索排列空间需要增大粒子数或初始权重。收敛曲线要是下降过程突兀说明粒子群在某代突然跳变多半是 gbest 更新的判断写错了把最差解当成最优解。时间窗是否被满足可以从车辆时间轴看出来横轴是时刻纵轴是车辆编号每个客户点显示服务开始时刻。如果某个点的柱形超出了客户窗口区域说明目标函数的超窗惩罚权重偏低。调整时先固定其他参数只动惩罚系数观察整体超窗时长是否随系数上升而下降通常 3 次对比就能确定合适区间。6. 把 PSO-TWVRP 改到真实场景收敛验证与局部搜索6.1 用最近邻种子粒子替代全随机初始化PSO 对初始群体非常敏感。全随机初始化会让种群花掉几十代才凑出几条可行路线尤其时间窗窄的时候更明显。我惯用的做法是在初始化时生成 10% 的种子粒子用最近邻构造客户顺序。比如对仓库随机选一个未服务客户开始每次选择距离当前点最近且满足时间窗的客户加入路线直到车辆容量用完再开下一辆车。这样得到的种子粒子已经是一条合法路线能极大拉低第一代 gbestCost。% 最近邻构造种子路线其余粒子仍随机生成 seed nearestNeighborRoute(data); x(1, :) encodeRoute(seed, data.nCustomer); for i 2:nPop x(i, :) rand(1, dim); endencodeRoute把离散路线反向映射成连续位置向量最笨的方法是按客户出现顺序倒推出一个严格递增的数值序列再穿插分隔维度。对大多数规模在 30 个客户以内的 TWVRP10% 种子粒子不会引起早熟却能让收敛曲线从第一代就开始下降。6.2 用 2-opt 在 gbest 上做局部精化PSO 的全局搜索强但局部精化弱。gbest 可能已经很接近最优路线却因为一次 swap 操作能省掉几十公里而未被发现。常见做法是在每次迭代结束后对 gbest 解码出来的每一条子路线执行 2-opt反转子路线中的任一段如果总成本下降且时间窗合规则接受。function newRoute twoOpt(route, data) improved true; while improved improved false; for i 1:length(route)-1 for j i1:length(route) cand route; cand(i:j) cand(j:-1:i); if isValid(cand, data) routeCost(cand, data) routeCost(route, data) route cand; improved true; end end end end newRoute route; end这段 2-opt 的复杂度是 O(L^2)L 是单辆车客户数放在主循环里每代做一次会增加不少时间。实际使用时可以每隔 5 代才精化一次或者只在输出最终 gbest 前做一轮效果接近且时间开销小得多。6.3 验证收敛是否可靠固定随机种子做的三件事改完参数别急着验收先固定随机种子然后做三件事跑三遍同样参数对比 gbest 和平均成本把客户数减半做一次穷举或动态规划验证最优值的量级把时间窗全部放宽成 0 到当天结束检验 PSO 能否退化成普通 VRP 并得到接近最近邻的结果。这三件事能把参数误设和索引越界问题暴露出来。最后检查 getCost 里惩罚系数确保超窗成本在数值上不会大到掩盖距离项导致算法只保时间窗不保路径质量。本文还有配套的精品资源点击获取