用粒子群算法求解TWVRP的MATLAB实现与参数调优指南

简介：基于Matlab粒子群算法求解带时间窗车辆路径规划问题（TWVRP）的完整源码包，面向物流调度、运筹优化与智能算法学习者。针对单仓库多客户、硬时间窗约束场景，提供从问题建模到PSO迭代求解的整套Matlab实现，涵盖速度-位置更新、适应度计算、时间窗校验等关键模块，并配有运行结果图与说明文档，便于对照复现与二次开发。压缩包共12个文件，包含5个.m源码文件、5张结果图、1个md说明及1个数据文件，整体仅229KB，轻量易用。目前已有421人学习/下载。通过该资源可快速理解粒子群算法处理带时间窗VRP的编码思路与参数设置，适合课程设计、毕业设计及物流路径规划入门实践。

1. 用粒子群算法求解TWVRP：为什么MATLAB是最合适的落地工具

带时间窗的车辆路径规划问题（TWVRP）和普通VRP只差一个约束：每个客户必须在指定时间段内被服务。正是这一个约束，使得经典的节约算法和最近邻贪心几乎全部失效，因为“先到优先生效”的假设不再成立。面对这种NP难问题，工程上常见的做法是放弃精确求解，转而使用元启发式算法，其中粒子群算法（PSO）因实现简单、不依赖梯度信息、参数少而成为入门首选。这篇文章以“单仓库、多客户、多车辆”作为默认场景，从数学建模、排序编码、适应度函数设计到MATLAB源码和调参技巧，完整走一遍用粒子群优化算法求解TWVRP的落地路径。适合需要快速产出可运行原型，又不想被商用求解器绑定在特定平台上的工程师。

2. TWVRP问题建模与粒子群算法原理

2.1 单仓库多客户TWVRP的数学模型与目标函数

TWVRP的标准定义建立在带权完全图上，顶点集合V由仓库节点0和N个客户节点构成，边权通常取欧氏距离。每辆车的路线都从仓库出发，按某种顺序访问分配给它的客户，最后回到仓库。与基础VRP的区别在于，每个客户i附加了一个时间窗口[ei, li]，车辆到达早于ei必须原地等待，晚于li则构成违约。如果所有时间窗都必须是硬约束，那问题就变成带紧约束的组合优化，搜索空间里可行域很小；如果是软时间窗，则把“违约时长”计入目标函数，让算法在距离最优和时间窗违禁之间做权衡。

目标函数设计在不同论文里差异很大，但核心都可以写成：

min z = totalDistance + alpha * timeWindowViolation + beta * vehiclesUsed

其中alpha是时间窗惩罚系数，beta是车辆固定成本折算系数。很多初学实现只算总距离，忽略车辆数，结果算法会把N个客户拆成N条单人路线，总距离反而最小——因为每辆车只服务一个客户时路径就是点到点直线，距离最少但成本最不合理。因此，工程上做TWVRP时要把车辆数纳入目标，另一个更隐蔽的点是：软时间窗场景下alpha不能设得过大，否则粒子全部被推向“宁可绕远路也要压线到达”的解，整体路径质量反而变差。建模阶段就明确目标函数里的三个组成部分，比跑完看结果再去猜惩罚系数要高效得多。

2.2 粒子群优化算法的速度-位置更新公式

粒子群算法原理的核心是模拟鸟群觅食时的信息共享机制。每个粒子代表一个候选解，在D维搜索空间里有一个位置向量x和一个速度向量v，每次迭代依据自身历史最优pbest和全局最优gbest来调整飞行方向。速度和位置的更新公式为：

v_i(t+1) = w * v_i(t) + c1 * r1 * (pbest_i - x_i(t)) + c2 * r2 * (gbest - x_i(t)) x_i(t+1) = x_i(t) + v_i(t+1)

w是惯性权重，控制粒子延续上一轮速度的程度；c1是自我认知系数，c2是社会认知系数；r1、r2是[0,1]区间均匀分布的随机数。从粒子群算法原理的角度看，w大则全局探索能力强，容易跳出局部最优；w小则局部搜索精细，收敛快但易早熟。常见做法是让w随迭代线性递减，例如从0.9降到0.4，在MATLAB主循环里用一行代码实现：

w = wMax - (wMax - wMin) * it / maxIt;

学习因子c1、c2一般取1.5左右，如果对收敛速度要求高，可以让c1略小于c2，即让粒子更信赖群体的判断。即使公式本身很简单，应用到TWVRP时最大的问题不是迭代公式，而是位置向量和客户路径之间的映射——这是下一小节的核心。

2.3 排序解码：连续位置向量到客户访问序列的映射

PSO天然工作在连续变量空间，而TWVRP的解是离散的客户排列，因此必须设计解码器。最常用的方法是排序解码法（random key encoding）：假设有N个客户，粒子的位置向量就是N个连续实数，例如[3.1, 0.7, 5.8, 2.2, 4.9]。把这五个数按升序排序，排序后的索引顺序就是客户访问顺序。在MATLAB里只需要两行：

[~, order] = sort(position); order = order - 1; % 如果客户编号从1开始，减1后变为客户编号

这样做有两个好处：一是连续变量的编码天然兼容PSO的速度和位置更新公式，不需要额外的离散化操作；二是排序解码是双射关系，任何一个连续向量都唯一对应一个排列，不会出现传统0/1编码那种大量重复解占用种群的问题。实际项目里，我还习惯在排序前对位置向量做一次归一化处理，避免数值过大或过小导致排序结果趋于随机、失去粒子之间的差异信息。注意一个细节：排序得到的是访问顺序，而车辆如何分派、时间窗如何校验，要全部放到适应度函数中统一计算。解码和成环两者必须严格对应，否则画出来的路径图和算出来的目标值会互相矛盾。

3. MATLAB源码实现：数据准备、适应度函数和PSO参数

3.1 随机实例数据生成：仓库、需求与时间窗数据结构

在写优化算法前，先要把TWVRP的测试数据组织好。开发阶段不建议直接使用标准算例，因为问题规模大、不好定位bug；正确做法是先程序化生成一个小规模随机实例，比如20个客户，把每一行数据都打印出来人工核对。下面这段MATLAB源码生成仓库加客户的完整数据结构：

%% 生成TWVRP测试实例：1个仓库 + 20个客户 rng(42); % 固定随机种子，保证结果可复现 N = 20; depot = [50, 50]; coords = [depot; rand(N,2)*100]; % 第1行是仓库，其余为客户 demand = [0; randi([2,8], N, 1)]; % 需求量 serviceTime = [0; randi([5,15], N, 1)]; % 服务时间 % 客户时间窗 [最早到达, 最晚到达] timeWin = zeros(N+1, 2); timeWin(1,:) = [0, 1e4]; % 仓库不限制 for i = 2:N+1 early = randi([20, 60]); late = early + randi([30, 80]); timeWin(i,:) = [early, late]; end % 车辆容量 Q = 30; % 速度换算系数：距离单位/时间单位 v = 1; % 距离矩阵 dist = pdist2(coords, coords);

这段代码的逻辑澄清：demand和serviceTime的首元素设为0，是为了让矩阵索引和前面coords一致，仓库节点在所有矩阵中都是第1行。时间窗里仓库用[0, 1e4]相当于无约束，客户时间窗长度设为30～80个时间单位，服务时间5～15，这组参数会让问题有一定拥挤度，方便后续观察算法行为。另外，pdist2属于MATLAB统计工具箱，如果只用基础版MATLAB且没有该工具箱，改用两层for循环配合norm函数计算距离矩阵即可，数据量在几百个客户内性能差距可以忽略。

3.2 适应度函数：时间窗惩罚、容量约束和车辆数折算

适应度函数是TWVRP求解的灵魂，所有约束都要在这里变成成本。输入是一个客户访问顺序order，内部根据容量和时间窗把顺序切分成若干条车辆路径。

function cost = twvrpFitness(order, coords, demand, serviceTime, timeWin, Q, v) dist = pdist2(coords, coords); n = length(order); totalDist = 0; % 总行驶距离 twViol = 0; % 时间窗违约量 load = 0; % 当前车辆载重 vehicles = 1; % 已派车辆数 curTime = 0; % 当前累积时间 curPos = 1; % 当前位置索引，1是仓库 for i = 1:n custIdx = order(i) + 1; % 客户编号转矩阵索引（+1偏移） arrive = curTime + dist(curPos, custIdx) / v; % 容量超限 或 到达时间迟于最晚时间窗 => 当前车回库，换新车 if load + demand(custIdx) > Q || arrive > timeWin(custIdx, 2) totalDist = totalDist + dist(curPos, 1); curTime = curTime + dist(curPos, 1) / v; load = 0; vehicles = vehicles + 1; curPos = 1; arrive = curTime + dist(curPos, custIdx) / v; % 重算到达时间 end % 早到等待 if arrive < timeWin(custIdx, 1) curTime = timeWin(custIdx, 1); else curTime = arrive; end % 软时间窗惩罚 if curTime > timeWin(custIdx, 2) twViol = twViol + (curTime - timeWin(custIdx, 2)); end load = load + demand(custIdx); totalDist = totalDist + dist(curPos, custIdx); curPos = custIdx; end % 最后一辆车返回仓库 totalDist = totalDist + dist(curPos, 1); % 目标函数：距离 + 时间窗惩罚 + 车辆数惩罚 cost = totalDist + 200 * twViol + 30 * vehicles; end

参数说明：200是时间窗惩罚系数，设得比车辆数惩罚大一个量级，确保算法优先压掉时间窗违约；30是车辆数惩罚，它在目标函数中的实际作用是把“多派一辆车”折算成相当于30个单位的行驶距离。当客户数增多、路线距离变长时，这两个系数需要同步放大，否则车辆数惩罚会变得无足轻重。还有一处需要注意：车载容量不足或时间窗违约时强制回库换车，这属于贪心式的路径切分，得到的车辆数不保证全局最优，但胜在计算简单且易于工程实现，这也是目前绝大多数TPSO论文使用的做法。

3.3 关键参数表和初始化代码

PSO求解TWVRP时的参数设置直接决定收敛质量，下面这张表格是经过多组实验验证的起始值，适合作为第一次运行的基线：

参数初始化代码：

nVar = N; nPop = 40; maxIt = 200; wMax = 0.9; wMin = 0.4; c1 = 1.5; c2 = 1.5; vMax = nVar / 2; X = rand(nPop, nVar) * nVar; % 位置向量 V = zeros(nPop, nVar); % 速度向量 pBest = X; % 个体历史最优 pBestCost = inf(nPop, 1);

位置向量初始化为[0, N]区间内的随机实数，覆盖整个排序空间。位置更新的边界处理用截断策略：超出上限拉到N，低于下限拉到0，不要用弹性边界或对称映射，否则粒子大量聚集在边界端点，排序时会出现很多“并列最小值”，导致群体多样性快速丧失。如果多次运行后发现某次收敛曲线的最终值明显差于其他次，多半就是初始化时随机种子落在了一个拥挤区域，把nPop调大或增加迭代次数是最直接的缓解手段。

4. 主循环迭代、路线构建与结果可视化

4.1 主循环代码与线性递减惯性权重策略

PSO求解TWVRP的主循环就是反复执行“评估适应度 → 更新pbest和gbest → 更新速度和位置”三个动作。下面的代码段把三个动作写在一段MATLAB源码中，便于整体阅读和修改：

bestHistory = zeros(maxIt, 1); % 记录每轮全局最优 for it = 1:maxIt w = wMax - (wMax - wMin) * it / maxIt; % 线性递减惯性权重 % 计算所有粒子适应度，更新个体最优 for k = 1:nPop [~, order] = sort(X(k, :)); c = twvrpFitness(order, coords, demand, serviceTime, timeWin, Q, v); if c < pBestCost(k) pBestCost(k) = c; pBest(k, :) = X(k, :); end end % 更新全局最优 [gbestCost, gidx] = min(pBestCost); gbest = pBest(gidx, :); bestHistory(it) = gbestCost; % 同步更新全体粒子速度与位置 for k = 1:nPop V(k, :) = w * V(k, :) ... + c1 * rand * (pBest(k, :) - X(k, :)) ... + c2 * rand * (gbest - X(k, :)); V(k, :) = max(min(V(k, :), vMax), -vMax); % 速度限幅 X(k, :) = X(k, :) + V(k, :); X(k, :) = max(min(X(k, :), nVar), 0); % 位置边界截断 end end

这里采用的是同步更新模式，即先完成所有粒子的适应度评估，再统一更新速度和位置。另一种做法是异步更新，每算完一个粒子就立刻刷新gbest，收敛速度稍快但实现复杂度高，可复现性反而差。速度限幅vMax = N/2这个值需要根据客户规模调整：如果vMax设得比N还大，粒子每一轮的位移会大范围跳变，sort解码出来的排列顺序频繁从头洗牌，早期的优秀局部结构很难保留；如果vMax小于N/4，粒子又容易被困在一个小区间内无法逃逸。线性递减的惯性权重从0.9起步，让前期粒子大范围探索，后期小步慢跑精细收敛，这是工程上最稳的选择。

4.2 路线分解函数与多车辆路径还原

主循环里twvrpFitness只返回一个成本值，拿到order后想画车辆路径图，还要有一个和适应度函数分车逻辑完全一致的路线分解函数。两者的判定条件必须在代码上逐字对应，否则画出的路线和算出的成本对不上。

function routes = buildRoutes(order, coords, demand, serviceTime, timeWin, Q, v) dist = pdist2(coords, coords); routes = {}; % 元胞数组，每项存储一条路径的客户编号 currentRoute = []; load = 0; curTime = 0; curPos = 1; for i = 1:length(order) custIdx = order(i) + 1; arrive = curTime + dist(curPos, custIdx) / v; % 和适应度函数保持一致：容量超限或迟到则回库换车 if load + demand(custIdx) > Q || arrive > timeWin(custIdx, 2) if ~isempty(currentRoute) routes{end+1} = currentRoute; end currentRoute = []; load = 0; curTime = curTime + dist(curPos, 1) / v; curPos = 1; arrive = curTime + dist(curPos, custIdx) / v; end if arrive < timeWin(custIdx, 1) curTime = timeWin(custIdx, 1); else curTime = arrive; end currentRoute = [currentRoute, order(i)]; load = load + demand(custIdx); curPos = custIdx; end if ~isempty(currentRoute) routes{end+1} = currentRoute; end end

这里面有一个效率隐患：currentRoute使用方括号拼接，每加一个元素都会重新分配内存，客户数量超过几百时速度明显下降。如果目标是一两千个客户，建议预先分配路由矩阵，记录每条路径的起点和终点，再在气泡外补齐客户。另外，返回的routes元胞数组里每一项都是一个客户编号向量，首尾都不含仓库，后续绘图时手动把仓库节点插到路径两端即可。建议把twvrpFitness中的分车逻辑抽成一个单独函数，在两个地方复用，避免出现“适应度函数换车条件改了一行，路线分解函数没同步改”这种隐秘bug。

4.3 收敛曲线、路径图与导出

跑完主循环后，先看收敛曲线，再画路径图，能快速判断算法是否进入了正常优化的轨道：

% 绘制收敛曲线 figure; plot(1:maxIt, bestHistory, 'b-', 'LineWidth', 1.5); xlabel('迭代次数'); ylabel('全局最优成本'); title('PSO求解TWVRP收敛曲线'); grid on;

收敛曲线从高位快速下降、中后段趋于平缓是正常表现。如果200次迭代末尾还在明显下降，说明maxIt不够，再追加100次迭代；如果前30次就不动了，说明初始化或参数有问题，优先检查vMax和惯性权重是否设置过小。

路径图绘制代码：

[~, bestOrder] = sort(gbest); routes = buildRoutes(bestOrder, coords, demand, serviceTime, timeWin, Q, v); figure; hold on; plot(coords(2:end, 1), coords(2:end, 2), 'ko', 'MarkerFaceColor', 'k'); plot(depot(1), depot(2), 'rs', 'MarkerSize', 12, 'LineWidth', 2); colors = lines(length(routes)); for r = 1:length(routes) nodes = [1, routes{r} + 1, 1]; % 首尾插入仓库索引1 plot(coords(nodes, 1), coords(nodes, 2), 'Color', colors(r,:), 'LineWidth', 1.2); end % 坐标轴等比例显示，避免路径变形 axis equal; % 导出高清图片，便于写报告 exportgraphics(gcf, 'twvrp_routes.png', 'Resolution', 300);

axis equal这行千万别删，去掉后X轴和Y轴单位长度不一致，本来笔直的线路会被拉伸成弧线，干扰路径交叉的判断。exportgraphics在R2020a及以后版本可用，老版本MATLAB改用saveas。路径图上如果有多条线路大量交叉且绕行明显，可以先做一步简单的2-opt局部搜索，把同一路径内两点之间的交叉边交换后再评估成本，往往能提升2%到5%的解质量。

5. 参数调优、易踩的坑和Solomon基准验证

5.1 三个优先调节的参数

TWVRP场景下最优先调节的三个参数分别是惯性权重w、时间窗惩罚系数alpha和车辆数惩罚beta。w决定粒子的探索和开发平衡，w过大时收敛曲线会被拉出明显台阶，w过小则前期快速陷入局部最优，观察整条曲线形态后按0.05的步长微调。alpha的调节要看惩罚后的解是否仍然存在时间窗违约：如果twViol始终不为零，按步长0.5增大；如果完全无违约但总距离偏高，则适当减小到原来的0.8倍。beta的调试有个规律——当车辆数分布极不均匀（有的车装了1个客户，有的装了满额），说明beta太大或太小，参照最优解的车辆数平均值来反向微调。记住一条原则：每次只动一个参数，四到五次实验后就能锁定合理区间。

5.2 最容易出问题的三处实现细节

第一处是索引偏移。MATLAB矩阵索引从1开始，客户编号通常从1到N，仓库在第1行，解码生成的order中对应客户编号，传入适应度函数后必须执行custIdx = order(i) + 1，否则仓库节点会被当成普通客户参与服务。第二处是时间单位不统一。如果客户坐标是经纬度，而时间窗按分钟定义，直接用欧氏距离除以车速会得出完全错误的服务时长，先把经纬度投影为平面坐标，再统一时间单位。第三处是粒子位置大量重叠。排序解码法下位置向量的绝对值没有含义，只有相对顺序起作用，如果vMax过小导致位置集中在极窄区间，排序产生的排列就极度相似，群体多样性崩塌。解决方法是定期检测种群位置的标准差，低于某个阈值就随机重置部分粒子的位置。

5.3 用Solomon基准算例验证求解质量

验证PSO求解TWVRP是否正确，单靠肉眼观察路径图是不够的。工程上最稳妥的做法是用Solomon基准实例的子集做对比实验，C101和R101的25客户版本是行业内最常用的入门验证题目，已知最优解的数值公开可查。跑10次取最优值、平均值和最差值，与已知最优解对比：差距在5%以内说明算法骨架正确，15%以上说明编码或适应度函数有结构性问题，优先回头检查索引偏移和分车逻辑。

提示：对比时要把“算法单次最强解”和“算法平均水准”分开记录，只报单次最好值是很多项目后期复现失败的根本原因。用相同的随机种子分布跑满10次，再记录标准差，才能说明算法的稳定性。

在MATLAB中实现标准实例读取时，把Solomon文件的客户坐标、需求、时间窗解析成与本文一致的数据结构，注意Solomon的仓库节点编号为0，读取后加1对齐。最终对比时把参考解的路径和PSO解叠加画在同一张图上，观察哪些路段的走向是稳定的公共边，哪些路段属于随机波动——这两类边分别对应问题骨架和算法噪声。后续想进一步优化，可以沿稳定公共边附近做精细搜索，对波动边做2-opt或交换重排，这种混合策略比单纯加迭代次数有效得多。

本文还有配套的精品资源，点击获取