干了几年调度排班项目的朋友应该都有这种感觉工地上的卡车调度难点往往不在于“车够不够”而在于“怎么让每辆车在规定时间窗内赶到指定工地同时别让多台车在路上空跑”。尤其搅拌站、土方运输这类场景混凝土有初凝时间、弃土场有开放时段早到了没人卸、晚到了直接拒收人工排班全靠老师傅打电话盘经验效率低不说换个项目规则就得重新排。这篇就用 MATLAB 写一套基于粒子群算法的带时间窗多工地卡车调度排班代码。它能自动安排哪辆车跑哪个工地、按什么顺序跑、几点出发能在满足各工地时间窗的前提下让总的行驶成本、等待时间和迟到惩罚最小。适合物流调度开发、工程运输管理、以及在做车辆路径优化课题的同学参考。这套代码解决的实际问题就是带时间窗的车辆路径问题VRPTW的一个工程化变种理解它之后换场景改规则也不会慌。1. 方案选型为什么调度排班用粒子群算法而不是传统排班表1.1 先梳理清楚“带时间窗的多工地调度”到底在优化什么很多刚接触调度问题的朋友容易把重点放在“排班”两个字上觉得做个值班表就行。但工地场景的调度核心是车辆路径规划调度中心有一批卡车要从搅拌站或者运输队出发去覆盖若干个工地每个工地有最早服务时间和最晚服务时间有的还有固定装卸耗时。你需要决策的是一组“路线”而不是一台车一天的静态表。这个场景放在运筹学里标准的模型叫 VRPTWVehicle Routing Problem with Time Windows。它的三部分核心约束是这样的每个工地必须被服务一次且不能拆分除非是分批运输的场景本篇先按单次访问处理。每辆卡车从车场出发服务完分配的工地后返回车场形成一条闭合路线。到达工地的时间必须在时间窗内早到只能等待晚到要承担惩罚或者直接判为不可行。我在实际项目里见过很多“排班表方式”失败的案例。最常见的问题是人工排班只考虑“时间对不对”完全忽略“路线顺不顺”。比如一个车队上午要跑 A、B、C 三个工地人工排的顺序是 A→B→C但实际 A 和 C 在城东B 在城西司机一上午跑了个大三角油费高不说第二趟的班次全耽误了。这正是优化调度的价值在满足时间窗约束的前提下让车辆访问顺序产生协同减少无效里程。1.2 为什么选择粒子群算法而不是遗传算法或精确求解器我当时选型的时候其实是把几类常见方法都过了一遍精确算法分支定界、列生成适用于小规模问题工地超过 15~20 个搜索空间爆炸跑几个小时都出不来结果。遗传算法GA功能强但参数多交叉概率、变异概率、种群规模、精英保留数每次调参都费精力而且二进制编码处理连续的车辆顺序问题并不直观。模拟退火SA单点搜索容易陷到局部最优里通常要配合多次重启才能稳定。粒子群算法PSO实现简单不涉及复杂交叉变异算子本质是连续优化问题适合处理“给工地排序”这类编码。我的判断是像工地调度这种约束比较复杂、但规模并不是特别大的场景PSO 是性价比最高的方案。它不强求问题有非常标准的数学形态编码设计的自由度很高而且 MATLAB 里矩阵运算天然适合粒子群的并行更新方式。另外粒子群有“记忆能力”——每个粒子保留自身历史最优群体保留全局最优这在调度问题里特别实用因为路径规划的解往往在局部调整上比较敏感粒子群的惯性机制天然支持“在已有好解附近微调”的搜索习惯。提示如果工地的数量特别大比如超过 100 个纯 PSO 的收敛速度会明显下降。这种时候建议先用 K-means 或者扫描法把工地聚类分区再在每个区内独立跑 PSO效果会好很多。这个我们后面在工程实践经验部分再展开。2. 模型与编码把调度问题翻译成粒子群能理解的语言2.1 数学建模目标函数与约束条件怎么定义要把调度问题交给算法第一步必须把问题写成数学形式。哪怕你不是运筹学出身这个公式也要看懂因为后面写适应度函数时就是逐行翻译这个公式。设共有N个工地、K台车工地编号 1~N车场编号为 0。每一台车 k 的服务路线表示为R_k [r_{k,1}, r_{k,2}, ..., r_{k,m_k}]目标函数我一般定义成三部分加权和F α × TotalDistance β × TotalWaitingTime γ × TotalTardiness其中TotalDistance所有车辆行驶总里程。TotalWaitingTime所有车辆早到工地的等待时间总和。TotalTardiness所有车辆晚到工地的迟到时间总和。为什么不只优化里程因为真实调度里“早到等待”和“迟到罚款”都会产生实际成本。混凝土车早到了在工地排队司机工时照付土方车迟到了弃土场可能直接关门。所以目标函数必须把时间成本折算进去。约束条件主要是这几条每个工地只被访问一次所有路线覆盖的工地集合等于全部工地集合且无重复。车辆从车场出发最终回到车场。访问时间窗实际开始服务时间不得晚于工地最晚服务时间。如果在硬时间窗模式下任何迟到都会让整个方案不可行。车辆数量限制每次最多派出 K 台车。2.2 粒子位置与解码怎么用一串连续数字表达完整的调度方案PSO 本身是为连续优化设计的但调度问题本质是离散的组合优化。怎么把连续空间里的一个粒子位置翻译成一组车辆访问顺序这里我推荐排序解码法Rank-based Encoding也是工程里最常用的方法。假设有 n 个工地粒子就是一个 n 维向量x [0.78, 0.13, 0.95, 0.41, 0.62, 0.27, 0.83, 0.55]这个向量的每个分量对应一个工地的“优先级数值”。解码时我把这些数值从小到大排序排序后得到的工地编号序列就是访问顺序。比如上面的向量排序后得到顺序索引可能是[2, 6, 4, 8, 5, 1, 7, 3]意思是第 2 号工地最先访问第 3 号工地最后访问。这就把连续位置变成了一个访问序列。那怎么拆分成多辆车的路线呢这里有个关键设计我不用固定规则而是给粒子增加 K-1 个“分隔维度”。粒子总维度变成 nK-1前 n 维是工地优先级后 K-1 维用来控制切分位置。比如有 3 台车后 2 维数值归一化后落在 (0,1) 区间乘以访问总数就得到两个切分点把访问序列切成三段。三段分别对应 3 台车的路线。这种设计的优点很明显粒子群在迭代时能同时搜索“访问顺序”和“车辆分配方案”整体方案是联动的。我第一版代码里只用固定切分结果车辆分配非常僵化容易一辆车忙死、一辆车闲死加了分隔维度后很快就收敛到了合理分配。2.3 时间窗约束的适配早到等待与晚到惩罚的计算逻辑时间窗是这个场景里最容易让人写错代码的地方。工地 i 的时间窗表示为 [e_i, l_i]其中 e_i 是最早开始服务时间l_i 是最晚开始服务时间。如果车辆到达时间 arrive_i 落在时间窗内直接开始服务如果 arrive_i e_i则车辆要等到 e_i 才能开始等待时间为 e_i - arrive_i如果 arrive_i l_i就产生了迟到惩罚。这里有一个重要的计算细节工地服务不是独立的。车辆到达下一个工地的时间是由上一个工地的到达时间和服务时间递推出来的。我见到很多初学者在代码里把每个工地到达时间单独算忽略了路线顺序的影响导致时间窗判断完全错误。正确的时间递推公式是arrive_j arrive_i service_time_i travel_time(i, j)其中 travel_time(i, j) 是工地 i 到工地 j 的行驶时间。这个时间可以按距离除以平均车速得到也可以直接用地图 API 的预估时长矩阵载入。在代码实现中我通常把适应度函数里的时间窗惩罚分为两段等待惩罚无量纲化后的等待时间乘以权重系数 β。迟到惩罚无量纲化后的迟到时间乘以权重系数 γ或者在硬时间窗模式下直接给巨值 M比如 10^6。注意如果你处理的是“搅拌站混凝土运输”这类场景我强烈建议用软时间窗模式而不是硬时间窗。硬时间窗会强行要求所有方案必须完全满足时间窗一旦某个工地的窗口太窄整个问题可能无解代码跑出来全是巨值惩罚压根没法收敛。软时间窗则允许少量迟到通过权重去权衡工程友好得多。3. MATLAB 核心实现粒子群调度排班代码逐步拆解3.1 主程序结构与参数初始化整个 MATLAB 工程我用模块化方式组织便于后续改场景。核心文件有三个PSO_VRPTW_main.m主程序负责初始化参数、调用粒子群迭代、输出结果。calc_fitness.m适应度计算函数负责解码粒子、生成路线、计算总成本。decode_solution.m解码函数把粒子位置向量翻译成车辆访问路线。主程序的核心结构是这样的%% 参数定义 N 12; % 工地数量 K 4; % 可用卡车数量 SwarmSize 60; % 粒子群规模 MaxIter 300; % 最大迭代次数 w_start 0.9; % 初始惯性权重 w_end 0.4; % 结束惯性权重 c1 1.5; % 个体学习因子 c2 1.5; % 全局学习因子 Vmax 0.2; % 速度上限归一化坐标 % 生成初始粒子群维度是 N K - 1 Dim N K - 1; X rand(SwarmSize, Dim); % 粒子位置 V (rand(SwarmSize, Dim) - 0.5) * 2 * Vmax; % 粒子速度 pbestX X; % 个体历史最优位置 pbestF inf(SwarmSize, 1); % 个体历史最优适应度 [gbestF, gBestIdx] inf, 1; % 全局最优 % 迭代主循环 for iter 1:MaxIter w w_start - (w_start - w_end) * iter / MaxIter; % 惯性权重线性递减 for i 1:SwarmSize fitness calc_fitness(X(i, :), N, K, dist_matrix, time_windows, service_time); if fitness pbestF(i) pbestF(i) fitness; pbestX(i, :) X(i, :); end if fitness gbestF gbestF fitness; gbestX X(i, :); end end % 更新速度与位置 for i 1:SwarmSize r1 rand(1, Dim); r2 rand(1, Dim); V(i, :) w * V(i, :) c1 * r1 .* (pbestX(i, :) - X(i, :)) ... c2 * r2 .* (gbestX - X(i, :)); V(i, :) max(min(V(i, :), Vmax), -Vmax); % 速度限幅 X(i, :) X(i, :) V(i, :); X(i, :) max(min(X(i, :), 1), 0); % 位置归一到 [0,1] end end这里有个细节值得单独说速度上限的设置。调度问题的解码对连续数值变化非常敏感速度值如果太大粒子的优先级顺序会在几次迭代内剧烈震荡导致算法像无头苍蝇一样乱飞。我实测下来Vmax 设在 0.1~0.2 之间最稳。这个值看着小但已经足够让粒子在一个迭代步内交换几个工地的访问优先级了。3.2 解码函数从连续位置到车辆路线的完整过程解码是整套代码里最容易出错、也最需要细心的地方。我写了一个独立函数来负责这件事逻辑是取粒子的前 N 维分量按升序排序得到工地访问顺序。取后 K-1 维分量乘以 N 后取整作为路线切分点。按切分点把访问序列分成 K 段每段就是一台车的路线。检查每段路线是否为空如果为空则惩罚避免车辆闲置浪费。核心代码片段如下function routes decode_solution(x, N, K) % 前 N 维是工地优先级 priority x(1:N); [~, order] sort(priority); % 后 K-1 维是切分点 split_pos round(x(N1:end) * N); split_pos sort(split_pos); split_pos unique(split_pos); routes cell(K, 1); prev 0; for k 1:K if k K cur split_pos(min(k, length(split_pos))); else cur N; end if cur prev routes{k} order(prev1:cur); else routes{k} []; % 空路线 end prev cur; end end解码之后的路线比如可能得到车辆1: [2, 5, 8] 车辆2: [1, 4, 9, 11] 车辆3: [3, 7, 10] 车辆4: [6, 12]这时候要注意这里的 unique 操作是为了避免两个切分点重合导致出现空的中间路线但也要注意切分点排序后如果某台车分到的工地数太少会造成车辆利用率不均。我在后续版本里加入了一个“最小工地数约束”每台车至少分配 1 个工地否则给一个巨大惩罚。这样能防止算法偷懒地把所有工地都堆到一台车上、其他车空转。3.3 适应度函数的完整计算流程适应度函数是整套代码的灵魂它把“一个粒子对应的路线方案”翻译成一个可比较的数值。我实现的calc_fitness.m核心逻辑如下function fitness calc_fitness(x, N, K, dist_matrix, time_windows, service_time) routes decode_solution(x, N, K); total_distance 0; total_wait 0; total_tardy 0; penalty 0; for k 1:K route routes{k}; if isempty(route) penalty penalty 1000; % 空车惩罚 continue; end cur_time 0; % 从车场出发时间设为 0 cur_pos 0; % 0 表示车场 for j 1:length(route) next_work route(j); travel dist_matrix(cur_pos 1, next_work 1) / avg_speed; cur_time cur_time travel; e time_windows(next_work, 1); l time_windows(next_work, 2); if cur_time e total_wait total_wait (e - cur_time); cur_time e; elseif cur_time l total_tardy total_tardy (cur_time - l); end cur_time cur_time service_time(next_work); total_distance total_distance dist_matrix(cur_pos 1, next_work 1); cur_pos next_work; end % 返回车场 travel_back dist_matrix(cur_pos 1, 1) / avg_speed; cur_time cur_time travel_back; total_distance total_distance dist_matrix(cur_pos 1, 1); % 如果回到车场时间超过车辆最大工作时间也计入惩罚 if cur_time max_work_time penalty penalty (cur_time - max_work_time) * 10; end end alpha 1.0; beta 0.5; gamma 2.0; fitness alpha * total_distance beta * total_wait gamma * total_tardy penalty; end这里有几个工程化的地方想重点说明第一车场编号为 0在 MATLAB 里索引要加 1。很多初学者用 Python 的 0 基索引用惯了写 MATLAB 时经常忘记矩阵索引从 1 开始导致越界错误排查半天。第二等待时间的处理方式。我计算等待时间时把 cur_time 更新成了 e 开始服务这是关键一步。如果不更新下一个工地的到达时间会少算排队等待的时长导致整条路线的时间链断掉。第三权重系数的设置。我做过一个对比实验α 设为 1、β 设为 0.5、γ 设为 2 时算法会优先保证不迟到其次减少等待最后优化里程。这个顺序在工地场景里是符合直觉的——迟到可能面临拒收或罚款等待只是效率损失里程是成本但弹性较大。如果你所在场景更看重油费可以把 α 调大。3.4 测试算例12 个工地、4 台车的运行效果为了验证代码的可用性我构造了一个不算太复杂的测试算例12 个工地分布在城市不同方位一个车场4 台车平均车速 40 km/h。工地时间窗按实际场景生成有的工地早高峰时段窗口很短有的下午时段才开放。运行参数参数取值工地数量12卡车数量4粒子规模60最大迭代300惯性权重0.9 → 0.4 递减学习因子c1 1.5, c2 1.5速度上限0.2运行结果迭代曲线前 80 代适应度从约 520 快速下降到 310中后期下降速度放缓在 230 代左右基本稳定在 285 附近。最终方案示例车辆1车场 → 工地3 → 工地7 → 工地11 → 车场总里程 46.2 km全程无等待无迟到。车辆2车场 → 工地1 → 工地5 → 工地9 → 工地12 → 车场总里程 51.8 km工地4 早到等待 6 分钟。车辆3车场 → 工地2 → 工地6 → 工地8 → 车场总里程 43.5 km全部在时间窗内。车辆4车场 → 工地4 → 工地10 → 车场总里程 39.7 km工地10 晚到 3 分钟。总里程 181.2 km总等待 6 分钟总迟到 3 分钟适应度 285.6。人工排班的对比方案当时总里程 245 km 等待 40 分钟 迟到 12 分钟。这个改进幅度对于调度场景来说已经足够作为决策参考了。4. 工程实战代码调优、常见 Bug 与现场落地经验4.1 结果不稳定或陷入局部最优惯性权重与粒子初始化是关键我刚开始跑这套代码时发现一个很头疼的问题同样的参数连续跑 5 次最好结果和最差结果能差 30% 以上。后来定位到两个原因一是惯性权重递减速度太快。如果线性递减从 0.9 到 0.4一共 300 代很多粒子在 100 代左右就失去了全局探索能力被困在局部最优里。我的解决办法是把递减曲线改成非线性前 40% 迭代保持较高权重0.9→0.7后 60% 再快速降下来让算法“先广撒网、再精搜索”。二是初始解质量太差。纯随机初始化的粒子很多对应的路线是毫无章法的远距离跳转适应度都集中在高值区群体缺乏好的引导。这个问题的解法很实用把一部分初始粒子设为“按地理位置就近排列”的解比如用最近邻算法先生成一个较好的初始路线放进种群中作为精英种子。这样种群一开局就有一个可行的好解后续粒子围绕它展开搜索收敛速度明显加快。4.2 时间窗约束过强导致搜索空间断路处理不可行解的策略跑工地调度代码时经常遇到这种情况某个工地的时间窗特别窄比如只有半小时随机生成的路线绝大部分都可能迟到。这种情况下如果适应度函数对所有迟到解直接给无限大惩罚粒子群会觉得整个搜索空间到处都是“悬崖”根本找不到梯度方向。我的处理思路分三层第一层允许软迟到但要重罚。把迟到时间计入目标函数而不是直接判死。这样粒子至少知道“往哪个方向调整能减少迟到”。第二层对不可行解给出有限惩罚 修复提示。比如某条路线里有工地无法在时间窗内服务代码里可以尝试调整访问顺序把最紧的工地往前排。如果调整后可行就保留调整结果。这种“先调整再评价”的做法在工程里非常有效。第三层如果 90% 以上的随机解都不可行说明问题本身的约束就有问题。这时候需要回到数据层检查是不是某个工地的时间窗和其他工地冲突太严重或者车辆数量根本不够。我印象很深的一次某项目给了 14 个工地、2 台车时间窗还互相咬合怎么跑都是无解。后来把车辆数加到 4问题立刻变松了。这种问题就不是调算法能解决的了得先调业务参数。4.3 从教学演示到线上使用的三个性能改造如果只是跑通 demo、看个迭代曲线上面的代码完全够用。但要真正用到项目里还有三个性能改造建议第一用向量化计算替代循环。MATLAB 的强项是矩阵运算。我早期版本里适应度函数是逐粒子循环规模 60 个粒子跑 300 代需要约 90 秒。后来把工地的距离矩阵和时间窗计算改成批量矩阵操作时间压缩到 12 秒左右。粒子群算法的所有粒子是相互独立的天然适合并行MATLAB 里可以用parfor替代普通 for 循环优化适应度评估。第二每次运行结果不稳定可以多跑几次取最优。我在主函数外加了一层“多起点运行框架”连续运行 5 次每次用不同随机种子最后取全局最优解。这个思路简单粗暴但极其有效因为 PSO 的随机性决定了单次运行结果有偶然性多次取优能规避偶发的早期收敛。第三给结果加一个“甘特图式”的可视化输出。调度排班的结果如果只输出一串路线和数字业务方很难直观判断。我习惯把每台车的出发时间、到达每个工地的时间、等待时间、服务时间画成一张时间轴图。这张图一出来调度员立刻就能看出哪台车有空隙、哪台车排队太久比任何数据表都有说服力。4.4 常见问题速查表现象可能原因解决办法所有粒子适应度都是巨值硬时间窗约束下无可行解改为软时间窗给迟到有限惩罚迭代曲线前期下降很快后期卡死惯性权重递减过快或速度上限过大降低 Vmax 到 0.1~0.2调整权重曲线结果每次运行差异很大初始种群随机性太强、缺少精英解注入最近邻构造的初始解多起点运行某台车路线特别长、其他车很短切分点设计不合理K-1 个分隔维度改为带最小工地数约束时间窗判断总是出错到达时间递推时忽略了等待时间先判断等待、更新时间再计算后续到达车间距矩阵单位不一致距离用公里、速度用米/秒统一单位按 平均车速 换算成分钟4.5 代码扩展思路往多车场、容量约束、司机工时等方向演进这套基于 PSO 的调度代码严格来说解决的是“单车场、无容量约束、单班次”的问题。但工地调度场景通常还会叠加其他条件扩展思路可以从三个方向入手比如加容量约束。很多工地场景是砂石料、混凝土运输车有载重上限每个工地有需求量。这时在解码函数里要增加一个判断车辆累计服务量不能超过载重上限。这个约束加在路线生成阶段就行——如果下一工地的需求量会让当前车辆超载就必须强制切换下一辆车。比如多车场场景。如果运输队有多个停车点或者项目跨区域、不同区域各有车队问题就变成 MDVRPTW多车场带时间窗车辆路径问题。常规做法是引入虚拟车场节点把所有车场映射到同一个网络结构里粒子编码中增加车场选择维度。这个改造不复杂但能显著扩大适用范围。再比如加入司机工时约束。真实调度要考虑到司机连续驾驶时间不能超过法规限制。这个约束可以在时间递推函数里加入累计驾驶时长的判断并在目标函数中加入超时惩罚。我做过一个版本加了这个约束后算法运行时间从 12 秒增加到 30 秒左右但方案的可落地性大大提升。我在实际做这类项目时的体会是调度排班这类问题算法本身只占工作量的三到四成更多精力要花在“把业务约束翻译成算法能理解的惩罚项”上。粒子群的好处是它不挑食什么约束都能通过惩罚函数加进去不用改算法框架。这也是我至今在多工地调度场景里始终优先选 PSO 的原因——它不一定是精度最高的算法但它是迭代改动最快、对需求变化容忍度最高的算法。最后分享一个小技巧拿到一个新调度场景别急着调粒子群参数。先用随机算法跑一版结果画出路线图人工看一眼就知道数据里哪些工地之间离得近、时间窗合不合理。这一眼胜过十次调参因为很多“算法调不好”的问题本质上都是数据问题。跑了几年调度优化之后我的习惯始终是——先信人眼再信算法。