在做配送调度和路径优化时“带时间窗的车辆路径问题”VRPTW是一个绕不开的经典问题客户有固定的服务时间段车辆有载重上限目标是用尽量少的车、跑尽量短的路把需求全部覆盖。这名字听起来专业其实拆开就是“几辆车、各跑一条线、每个客户卡着时间点送完货”。但真要把这个过程讲清楚、算明白光靠直觉是不行的必须把它写成数学模型交给求解器。这篇文章会把我自己完整跑通的一套路程整理出来如何用Python Gurobi 从零建模并求解VRPTW包括问题定义、数学公式、完整代码、结果可视化和性能调优。适合数学建模竞赛选手、物流调度相关开发人员以及所有想知道“优化算法到底怎么落地”的读者。代码我已经整理成可以直接复制运行的版本环境配好就能跑。1. 先明确VRPTW在解决什么问题定义与数学模型1.1 从现实场景理解VRPTW先说一个最典型的场景。某城市有1个配送中心、若干个客户点每个客户在某个时间段内收货比如“早上9点到11点必须送到”公司有若干辆载重相同的货车。现在要安排每一辆车的行驶路线使得所有客户都被服务车辆不超载并且车辆到达客户的时间不能晚于客户最晚可服务时间。这个问题的难点在于时间和路线互相牵制。A客户要求上午10点前到B客户要求下午2点后才方便如果先送A再送B中间等太久可能浪费反过来路线又影响到达时间。更麻烦的是车辆数量也是成本少一辆车能省一大笔固定开销但往往会让每辆车绕路更多、行驶距离更远。运筹优化要做的就是在这堆矛盾里找出一个全局最优解而不是靠经验拍脑袋。VRPTW在学术上属于车辆路径问题VRP的核心变种在物流配送、外卖调度、维修工单派单、校车接送、移动充电车调度里都有应用。数学建模竞赛里也反复出现这类题目的影子所以我一直建议接触优化方向的同学把VRPTW当成一个标准入门案例去啃。模型拆明白了后面再遇到多车型、多配送中心、取送一体等变种都是在这个骨架上加变量加约束的事。1.2 建模五要素集合、参数、变量、目标、约束新手建模最容易犯的错是上来就写代码。正确的顺序一定是先把数学模型写清楚哪怕只是写在纸面上。VRPTW的标准建模分五个部分。集合与参数节点集合0表示配送中心depot1到n表示客户。车辆集合K {1, 2, ..., m}所有车辆同质最大载重为Q。每个客户i有需求量q_i服务时间s_i服务时间窗[e_i, l_i]e_i是最早允许开始服务的时间l_i是最晚允许开始服务的时间。节点i到节点j的距离d_ij行驶时间t_ij。本文假设速度恒定所以直接t_ij d_ij。决策变量0-1变量x_ij^k车辆k是否从节点i直接行驶到节点j。连续变量T_i^k车辆k到达节点i并开始服务的时间。注意这里说的是“开始服务的时间”如果车辆早到了可以等无所谓。目标函数VRPTW的目标通常有两种写法最小化总行驶距离或者先最小化使用车辆数再最小化总距离。实际业务中车辆固定成本往往远大于油耗和里程成本所以优先减少车辆数更合理。数学上可以用多层目标处理Gurobi的setObjectiveN就是干这个的。这篇文章的主代码先按“给定车辆数最小化总行驶距离”来写后面在第4部分讲怎么改成分层目标。约束条件约束是模型的核心。按类划分每个客户恰好被一辆车服务一次。流守恒车辆到达某个客户后必须离开不能停在客户那里。车辆从配送中心出发最终必须返回配送中心。载重约束一条路径上的总需求量不能超过Q。时间窗约束每条弧上满足时间推进关系每个客户的开始服务时间落在[e_i, l_i]内。这个模型听着不难但把时间窗约束写成线性表达式的时候很多第一次接触的人会卡住。1.3 时间窗约束的线性化与子回路处理先看这段逻辑。如果车辆k从i直接开到j那么到达j并开始服务的时间一定不能早于“到达i开始服务的时间 服务时长 i到j的行驶时间”。写成逻辑表达式如果 x_ij^k 1则 T_i^k s_i t_ij ≤ T_j^k但Gurobi解决的是数学规划问题不能直接处理“如果……则……”这种条件逻辑。标准做法是引入一个足够大的常数M把它改写成T_i^k s_i t_ij ≤ T_j^k M(1 - x_ij^k)这个技巧叫大M法。当x_ij^k 1时不等式右侧的M(1 - x_ij^k) 0条件生效当x_ij^k 0时右边被放得很大约束自然松弛不影响其他解。用生活化的例子理解这相当于给约束加了一个“开关”路线上用了这条弧开关合上时间必须满足先后关系没走这条弧开关断开随便T_i和T_j怎么取值都不受约束。M怎么选是关键。太小会误伤可行解太大又会让模型数值变差求解器在分支定界时容易产生一堆靠近M的冗余计算。我习惯取M max(late) - min(early) max(service) max(d_ij)也就是一条路径上可能出现的最长时间跨度。这样M刚好够用又不会大到失控。实在不放心就取max(late) 100反正要比所有时间窗上界大出一截。再讲一个很多人不知道的细节。如果是不带时间窗的普通VRP模型里必须额外加“子回路消除约束”否则求解器会给出类似“1号车从客户2开到客户3再回到客户2”的非法回路。但VRPTW里时间窗约束天然做了这道防线只要出现子回路T在这些节点上会一圈一圈地累加服务时间和行驶时间很快就突破l_i子回路自动不可行。所以我在代码里没有单独写子回路消除这也是VRPTW比CVRP建模省事的地方。2. 环境准备Python、Gurobi与License实战2.1 为什么选Gurobi优化求解器里开源的有SCIP、OR-Tools、HiGHS商业的有Gurobi、CPLEX、Xpress。我平时做MIP混合整数规划项目Gurobi用得最多。原因很直接求解性能顶尖尤其在大规模整数规划上碾压大部分开源求解器Python接口设计得干净几乎和数学公式一一对应调试起来省心对学术界有免费license学生和老师用校园邮箱就能申请不用花一分钱。数学建模比赛里Gurobi的出场率也非常高老牌参赛选手基本人手一份。如果你的电脑上已经装了基于Matlab的YALMIP或CPLEX也不用担心冲突。Gurobi、CPLEX这类求解器是相互独立的软件各自有独立的license机制和接口装在同一台机器上互不干扰。平时用Python的话import gurobipy用的就是Gurobi不需要卸载其他任何求解器。2.2 安装与License激活实操环境配置这块很多人卡住的不是求解器本身而是Python环境混乱。我建议用conda或venv单独建一个项目环境别把包装到系统Python里。下面是一套我验证过的流程。第一步确认Python版本。Gurobi对Python版本有要求目前主流版本支持3.8到3.12建议用3.10或3.11兼容性最稳。第二步安装gurobipy。新版本Gurobi的Python包已经包含了完整求解器二进制直接pip装即可pip install gurobipy如果是通过Gurobi官方安装包安装的完整版也可以让Python自动找到gurobipy动态库。但pip方式最简单推荐。第三步申请并激活学术License。打开Gurobi官网注册账号注册时选择Academic用户用学校邮箱接收验证然后在用户后台申请学术license会得到一个license key。在命令行执行grbgetkey xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx按提示选择license文件存放路径一般默认放在用户主目录下生成一个gurobi.lic。激活完成后可通过环境变量GRB_LICENSE_FILE指向这个文件也可以用默认路径自动识别。第四步验证是否安装成功import gurobipy as gp print(gp.gurobi.version()) m gp.Model(test) x m.addVar(vtypegp.GRB.BINARY, namex) m.setObjective(x, gp.GRB.MAXIMIZE) m.optimize() print(x.x)如果输出版本号并且模型求解显示OPTIMAL说明环境没问题。很多同学会在无license状态下直接跑模型。Gurobi对未激活license的模型规模有限制通常只能跑最多2000个变量或200个约束的小规模问题再大就会报license错误。好消息是本文这个9客户的VRPTW案例规模很小即使没激活license也能求解只是控制台会提示使用限制。想认真做竞赛或项目还是建议把学术license激活。3. 核心代码实现从数据构造到模型求解3.1 构造测试数据我设计了一个小规模案例1个配送中心加9个客户。数据格式是(x坐标, y坐标, 需求量, 最早开始服务时间, 最晚开始服务时间, 服务时长)单位统一。import math import gurobipy as gp from gurobipy import GRB # (x, y, demand, early, late, service) data [ (50, 50, 0, 0, 999, 0), # 0: depot (37, 52, 7, 40, 80, 10), # 1 (49, 49, 3, 0, 60, 10), # 2 (52, 64, 5, 20, 90, 10), # 3 (20, 26, 6, 30, 100, 10), # 4 (40, 30, 4, 10, 70, 10), # 5 (61, 42, 5, 50, 120, 10), # 6 (35, 18, 8, 20, 90, 10), # 7 (47, 66, 2, 60, 130, 10), # 8 (52, 35, 9, 30, 90, 10), # 9 ] n len(data) - 1 N list(range(n 1)) C list(range(1, n 1)) K list(range(3)) # 3辆车 Q 25 # 单车载重上限 pos [(p[0], p[1]) for p in data] demand [p[2] for p in data] early [p[3] for p in data] late [p[4] for p in data] service [p[5] for p in data] def get_dist(i, j): return math.hypot(pos[i][0] - pos[j][0], pos[i][1] - pos[j][1]) d_ij {(i, j): get_dist(i, j) for i in N for j in N if i ! j} M max(late) 100 # 大M参数坐标直接用平面欧氏距离等于行驶时间这个假设在入门案例里非常常见。9个客户总需求量为49单车载重25理论上3辆车完全可以覆盖模型搜索可行解的余地是有的。时间窗设置得比较宽松目的就是让最优解的路由结构有更多组合空间演示效果更好。3.2 模型构建变量、目标与约束这一部分是全文的核心。我先把模型对象、决策变量和目标函数立起来。model gp.Model(VRPTW) # 决策变量 x[i,j,k]车辆k是否走弧(i,j) arcs [(i, j, k) for i in N for j in N if i ! j for k in K] x model.addVars(arcs, vtypeGRB.BINARY, namex) # 决策变量 s[i,k]车辆k到达节点i并开始服务的时间 s model.addVars(N, K, vtypeGRB.CONTINUOUS, names) # 目标最小化总行驶距离 model.setObjective( gp.quicksum(d_ij[i, j] * x[i, j, k] for i, j, k in arcs), GRB.MINIMIZE )这里有个细节arcs列表直接排除了i j的情况也就是不允许车辆从某节点开回自身避免出现无意义的“原地移动”变量。这比生成变量后再手动加约束x[i,i,k] 0要干净。接下来是约束。我直接把全部约束按类写在一个代码块里注释标清楚。# 约束1每个客户恰好被一辆车服务一次 # 等价于客户i必须有一条离开弧 for i in C: model.addConstr( gp.quicksum(x[i, j, k] for j in N if j ! i for k in K) 1 ) # 约束2流守恒进入客户h的车辆必须从h离开 for h in C: for k in K: model.addConstr( gp.quicksum(x[i, h, k] for i in N if i ! h) gp.quicksum(x[h, j, k] for j in N if j ! h) ) # 约束3每辆车最多从depot出发一次且出车后必须返回depot for k in K: model.addConstr( gp.quicksum(x[0, j, k] for j in C) 1 ) model.addConstr( gp.quicksum(x[i, 0, k] for i in C) gp.quicksum(x[0, j, k] for j in C) ) # 约束4载重约束 for k in K: model.addConstr( gp.quicksum(demand[i] * gp.quicksum(x[i, j, k] for j in N if j ! i) for i in C) Q ) # 约束5时间窗约束大M法 for k in K: for i in N: for j in N: if i ! j: model.addConstr( s[i, k] service[i] d_ij[i, j] s[j, k] M * (1 - x[i, j, k]) ) # 约束6服务开始时间必须落在时间窗内 for i in N: for k in K: model.addConstr(s[i, k] early[i]) model.addConstr(s[i, k] late[i])约束1要稍微解释一下。我没有写成“每个客户必须有一条进入弧”而是写成“客户i必须有一条离开弧”。由于约束2流守恒保证“进来多少就得出去多少”进入客户i的弧也必然存在。两个条件合在一起就确保了客户被服务一次且路径连贯。约束3允许车辆不出车。如果车辆k没用那它不产生任何从depot出发的弧也不产生返回弧如果用了就必须有一条离开depot的弧和一条回到depot的弧。这个写法比强制“每辆车必须出车”更贴近真实业务因为很多时候可用的车队是有冗余的。时间窗约束里的s[i,k]需要注意车辆可能提前到达在客户门口等着所以s[i,k]表示的是“开始服务时间”不是“到达时间”。车辆8点到了客户要求9点开始服务那s[i,k] 9等待的1小时不影响模型。3.3 求解、输出与可视化模型构建完成后直接调用optimize然后解析结果。解析最核心的一步是从0号节点出发沿着x[i,j,k].x 1的弧一路找下去直到回到0号节点得到每辆车的完整路径。model.optimize() if model.status GRB.OPTIMAL: print(f最优总行驶距离: {model.objVal:.2f}) for k in K: route [0] cur 0 while True: nxt None for j in N: if j ! cur and x[cur, j, k].x 0.5: nxt j break if nxt is None or nxt 0: break route.append(nxt) cur nxt if len(route) 1: load sum(demand[i] for i in route if i ! 0) print(f车辆{k1}路径: - .join(str(i) for i in route) f装载 {load}/{Q}) for i in route: if i ! 0: print(f 客户{i}: 开始服务时间 {s[i, k].x:.1f})这个解析循环用了一个小技巧只要某条弧上的x值为1就说明车辆走这条弧。Gurobi在MIP最优解里BINARY变量会给出0或1的精确值所以用0.5作为判断阈值是安全且常用的。如果想把路线可视化方便检查路径是否合理可以借助matplotlib画出节点分布和路线import matplotlib.pyplot as plt plt.figure(figsize(8, 8)) for idx, (xi, yi) in enumerate(pos): if idx 0: plt.scatter(xi, yi, cred, markers, s200, zorder5, labelDepot) else: plt.scatter(xi, yi, cblue, markero, s80, zorder5) plt.text(xi 1, yi 1, str(idx), fontsize10) colors [green, orange, purple] for k_idx, k in enumerate(K): route [0] cur 0 while True: nxt None for j in N: if j ! cur and x[cur, j, k].x 0.5: nxt j break if nxt is None or nxt 0: break route.append(nxt) cur nxt if len(route) 1: for t in range(len(route) - 1): a, b route[t], route[t 1] plt.plot( [pos[a][0], pos[b][0]], [pos[a][1], pos[b][1]], colorcolors[k_idx % len(colors)], linewidth2, ) plt.legend() plt.title(VRPTW Solution Routes) plt.show()运行之后你会在控制台看到3条路径并且所有客户都被服务一次每辆车装载不超过25。如果你调整某个客户的时间窗比如把客户9的最晚时间从90改成50路径大概率会发生变化这就是时间窗对路由结构的直接影响。这个实验我强烈建议新手自己动手做一遍比单纯抄代码理解深很多。4. 结果验证与性能优化4.1 结果合理性验证的三种方式很多入门者拿到OPTIMAL就以为万事大吉了。但模型输出“最优”不代表代码一定写对了模型本身可能有设计问题导致最优解不是真实问题的合理答案。我每建一个模型都会从三个维度做合理性体检。第一看硬约束是否满足。把打印出来的路径挨个检查每个客户是否只出现一次每辆车装载是否小于等于Q每个客户的服务时间是否落在这个客户的时间窗内。Gurobi本身不会破坏约束但如果约束漏写或写错它就会在“错误模型”里求最优。比如我刚开始学VRPTW时漏了流守恒约束结果解出来一辆车“陨落”在半路另一辆车“分身”去接续这种解一眼就能看出不对劲。第二观察路线是否“像人排出来的”。如果最优解里出现两条路线交叉严重、车辆绕远路、或者明明一次能送完的客户被拆到3辆车基本上可以怀疑目标函数或约束有问题。优化结果虽然不等同于人工经验但好的模型给出的解通常符合基本直觉。第三做灵敏度测试。把某个客户的时间窗收窄或者把某辆车去掉重新求解观察目标值和路径结构如何变化。如果时间窗收窄后目标值不变、路径不变说明该约束可能没有真正生效如果明显变化说明时间窗在模型里起作用了。这个测试特别能暴露“约束写错但没报错”的隐性bug。4.2 性能优化从通用模型到高效求解VRPTW看似简单一到大几十个客户、几十辆车模型规模迅速膨胀直接求解会非常吃力。我总结了四个实战中最有效的优化方向。第一打破车辆对称性。K辆车是同质的对模型来说“车辆1跑路线A、车辆2跑路线B”和“车辆1跑路线B、车辆2跑路线A”是两种候选解但实际是同一个方案。这会成倍增加搜索空间。解决方案是给车辆编号加一个先后约束for k in range(len(K) - 1): model.addConstr( gp.quicksum(x[0, j, k 1] for j in C) gp.quicksum(x[0, j, k] for j in C) )意思是编号大的车辆出车编号小的车辆至少也要出车。也就是“车辆1优先用车辆2其次车辆3最后”一下子把等价解削掉一大半。这个技巧在我处理多车模型时必加。第二提供一个好初始解。Gurobi自带启发式但很多时候一个手工构造的可行解能让求解器节省大量的初始搜索时间。最简单的方法是贪心从depot出发每次选择满足载重且时间窗允许的最近客户。把构造的路径转成x变量通过Gurobi的Start属性传入模型求解器会把它当作MIP start。第三调MIP参数。遇到中小规模但迟迟不收敛的情况我常用以下三个参数组合model.Params.MIPFocus 2 # 侧重证明最优性适合已有可行解但gap降不下去的情况 model.Params.TimeLimit 60 # 限制求解时间为60秒 model.Params.Threads 4 # 多线程加速如果当前连可行解都难找到更适合把MIPFocus设为1让求解器更早发现高质量可行解。第四当车辆数也是优化目标时不要直接用单目标加权比如写“总距离 1000 * 车辆数”。这种做法加权系数非常难调量纲也不统一。更推荐用Gurobi的分层多目标model.setObjectiveN( gp.quicksum(x[0, j, k] for j in C for k in K), index0, priority2, weight1.0 ) model.setObjectiveN( gp.quicksum(d_ij[i, j] * x[i, j, k] for i, j, k in arcs), index1, priority1, weight1.0 )priority数值越大优先级越高。Gurobi会先优化车辆数目标然后在车辆数最优的前提下优化总距离。这种方案在实际物流项目里几乎是标配。5. 实战中遇到的坑排查经验手册5.1 安装与License报错速查我在培训和答疑中遇到过大量环境问题把最高频的几个整理成了速查表。现象常见原因解决方法ModuleNotFoundError: No module named gurobipy包没有安装到当前Python环境检查是否用了正确的解释器执行pip install gurobipyGurobiError: license not foundlicense未激活或路径不对执行grbgetkey激活检查环境变量GRB_LICENSE_FILEFile grbgetkey not found没有安装Gurobi完整套件直接pip install gurobipy新版自带grbgetkeylicense无法激活没有用校园邮箱注册重新注册账号提交学校邮箱运行提示limited size用的是无license模式激活学术license或商业license如果Python和Anaconda并存最典型的坑就是pip装到了系统Python而Jupyter里用的是conda Python。解决方案是统一环境创建虚拟环境然后在虚拟环境里pip install gurobipy。在VSCode里也要确认右下角Python解释器选的是同一环境。别问我为什么知道问就是帮人排查过太多次。5.2 模型不可行的排查思路模型跑出来显示INFEASIBLE这不代表问题没有解更可能代表模型约束太紧或者写错了。我之前带竞赛组时学生最容易在三个地方把模型整到不可行车辆数太少。时间窗和容量叠加在一起可能3辆车就是覆盖不了9个客户。某个客户的时间窗设置不合理例如配送中心工作时间和客户时间窗完全没有交集。需求量超过单车载重或者数据录入时单位和换算出了错。排查不可行我强烈推荐用Gurobi的computeIIS功能。if model.status GRB.INFEASIBLE: model.computeIIS() model.write(model.ilp)computeIIS会找出一组“最小不可行约束集”写到model.ilp文件里。打开这个文件你能看到哪些约束组合在一起导致无解。把这些约束逐个放松或注释就可以定位问题。这里的技巧是对整数变量模型computeIIS有时会给出一个很大的冲突集不一定是最小的但依然有极强的参考价值。实操中我通常会先怀疑时间窗约束把大M值调大或放宽某几个客户的e_i、l_i再跑一次往往很快就能定位。5.3 时间窗建模的细节坑时间窗是最容易踩细节坑的地方我单独列几条。第一s[i,k]表示的是“开始服务时间”不是“到达时间”。如果客户最早9点开门车辆8点就到了那s[i,k] 9车辆等待1小时。这意味着模型允许早到等待但没有成本。真实业务里如果等待成本很高可以加一个“等待时间”变量并计入目标让模型尽量避免太早到达。第二depot时间窗的含义。depot的[e_0, l_0]代表车辆最早可发车时间和最晚必须回场时间。很多新手把所有节点同一个循环里约束时间窗结果把depot的l_0也设成某个客户的最晚时间导致车辆早出发被限制。正确做法是depot的l_0设为一天结束的horizon比如999。第三大M的选择会影响数值稳定性。M太小某些可行解被误杀M太大LP松弛非常弱求解速度显著下降。不要图省事直接写1e6你应该按时间跨度去算。前面的代码中M max(late) 100 230这个量级对这个案例来说完全够用。第四如果你想把“服务时间”纳入目标比如某些大件货物需要更长装卸时间那么服务时间通常还是作为参数s_i存在不会变成变量。只有当服务时间本身受资源影响时比如人手安排不同服务时长可变才需要升级成决策变量。入门阶段先保持它是固定参数除非题目明确要求。6. 从VRPTW继续扩展下一步可以做什么6.1 更接近真实场景的变种学会了基本VRPTW建模你已经掌握了一套适用范围很广的建模框架。真实业务里这个框架最常被扩展成以下方向多配送中心把单个depot扩展成多个depot需要为每个客户和车辆指定归属配送中心模型增加一个配送中心分配层。多车型不同车型载重、油耗、固定成本都不同车辆集合K换成异质的约束和目标要按车型区分。软时间窗客户可以接受晚到但要有惩罚成本。核心变化是时间窗约束从“必须满足”变成“违反就惩罚”目标函数里加一个晚到或早到的惩罚项。带回程取货车辆不仅要送货还要从客户处带回一些货物每辆车的总容量既要算送货量也要算取货量路径上还得考虑装载量的动态变化。考虑实际路网把欧氏距离换成真实道路网络的最短路需要用地图API或路网文件预计算OD矩阵。这些变种在实现上不会比基础版复杂太多核心加变量的思路完全一致。我每次给团队讲课都强调先保证基础模型能正确求解再一层一层叠加约束不要一上来就写大而全的模型。6.2 大规模问题怎么办当客户数量从几十涨到几百上千直接用Gurobi跑MIP基本扛不住。这时候有几条路可以走把模型改成“分簇路径”两阶段先用聚类把客户分给不同车辆再单独优化每一辆车的TSP路径。这种方法简单有效在很多实际项目里精度够用。用启发式算法如模拟退火、遗传算法、禁忌搜索、自适应大邻域搜索ALNS求高质量的近似解然后再用Gurobi在局部范围内做精确优化。用列生成或分支定价把路径作为列变量动态生成适合求解大规模VRPTW的精确解。这个方向数学要求偏高属于进阶内容。使用Gurobi的lazy constraint和pool search针对特殊结构问题做定制化分支切割。一个我在工业项目里反复验证过的组合思路是先用ALNS或贪心构造初始解再用Gurobi做局部重优化比如固定车辆数后对每辆车跑TSP。这样既发挥启发式在大规模上的速度优势又用Gurobi拿到了“局部精确”的可靠性性价比最高。根据自己的实际经验我在做优化项目时有一种很深的体会VRPTW这类问题的真正难点从来不是“会调用求解器”而是“能把现实问题准确翻译成约束和目标”。业务方一句“尽量别让客户等太久”落到模型里可能是很多种惩罚函数一句“司机不能开太远”可能是最大路径时长约束也可能是绕路惩罚。代码谁都会写但能把模糊需求拆成清晰的数学条件才是运筹优化工程师的核心竞争力。这篇文章里的代码只是开始你可以在它的基础上不断往上叠需求、加约束慢慢就会找到那种“拿到任何问题心里都有建模框架”的感觉。