我入行做运筹优化这些年最常被问到的问题就是VRPTW到底怎么用代码写出来很多朋友看论文公式头头是道一打开Gurobi就不知道变量怎么定义、约束怎么加。这篇内容我就用Python加Gurobi的组合把带时间窗车辆路径规划问题的完整建模过程从头到尾捋一遍从问题定义、模型推导到代码实现、结果判读全部覆盖适合数学建模竞赛选手、物流调度方向的研究生、以及刚接触运筹优化想找一份可跑通Demo的开发者参考。我尽量不写废话每步操作都说明背后的逻辑。项目本身不复杂但很多细节会让初学者卡个半天这些我都会单独拎出来讲。VRPTW本质上解决的是有一群客户每个客户有货物需求和服务时间窗车队从仓库出发每个客户只能被访问一次车辆有容量限制目标是让总行驶距离或总车辆数最小。这类问题在物流配送、外卖调度、维修工单派发里非常常见建模实战的价值远不只是应付比赛而是真正能迁移到业务场景里用的抽象能力。我先从模型定义开始把VRPTW这个问题的边界彻底讲清楚。1. VRPTW到底是什么这个模型在解决哪一层的问题1.1 问题定义与适用场景边界带时间窗的车辆路径问题英文全称是Vehicle Routing Problem with Time Windows标准定义是在一个配送网络里有一座仓库 depot 和若干个客户点 customer 每辆运输车有载重上限每个客户有确定的货物需求量并且客户只在特定时间段内接受服务。你要规划一组车辆路线让所有客户都被服务到同时满足车辆容量约束、时间窗约束和道路连通性约束。这句话拆开看每一部分都是约束的来源。车辆容量约束意味着一条路线上所有客户的需求总和不能超过单车载重时间窗约束意味着车辆到达客户的时间不能早于最早允许服务时间也不能晚于最晚允许服务时间如果到早了必须等待道路连通性约束则保证每条路线是从仓库出发最终回到仓库的完整回路。不是所有配送问题都是VRPTW。很多新手会把VRPTW和TSP、VRP、MDVRP搞混。TSP是单个销售员走遍所有城市没有容量和时间窗概念VRP是VRPTW去掉时间窗的版本MDVRP是多个仓库的版本。VRPTW在最基础的VRP之上叠加了时间维度的约束这使得问题的可行域被大幅压缩求解难度也明显上升。你需要先判断自己的业务场景是否真的存在时间约束如果只是单纯安排路线且没有客户时间要求直接建模成VRP会更高效。1.2 惩罚成本模型和硬时间窗模型的取舍VRPTW有两种常见建模方式硬时间窗和软时间窗。硬时间窗要求车辆到达时间必须在时间窗内早到只能等待晚到直接不可行软时间窗允许违反时间窗但产生惩罚成本。从数学建模和Gurobi求解的角度看硬时间窗是纯MIP约束约束表达式为等待时间加行驶时间的递推关系软时间窗则需要引入非负偏差变量在目标函数中新增惩罚项。两者在Gurobi里的实现差别很大。我个人的建议是比赛场景优先用硬时间窗。原因有两点第一硬时间窗模型更干净不需要调惩罚系数惩罚系数一旦设置不当结果可能为了省成本大面积超时业务上不可接受第二Gurobi处理硬时间窗约束时线性松弛的质量通常更好分支定界效率更高。业务落地场景则用软时间窗更实际因为真实道路状况下硬性时间窗往往过于理想化但那是另一个话题本文先用硬时间窗。1.3 目标函数的设计逻辑VRPTW最常见的双目标是最小化车辆数和最小化总行驶距离。实际操作中大家通常采用字典序优化先把车辆数压到最小再在最少车辆的前提下缩短总距离。字典序优化的做法是给两个目标设置不同的权重系数比如最小化M * 车辆数 总距离其中M取一个非常大的数保证优先优化车辆数。这个M怎么取有讲究它必须大于所有可行方案下总距离的最大可能差值一般取总距离上界的估计值即可比如所有客户到仓库往返距离之和。如果M取太小模型可能为了省几公里路程多派一辆车结果得到的解与预期完全不符。2. 数学模型精读约束公式背后的物理含义与表达陷阱2.1 参数、集合与决策变量的完整定义在写代码之前必须先把数学模型用标准MIP形式表达清楚。数学符号虽然啰嗦但它是调试代码的唯一依据。集合定义节点集合N {0, 1, 2, ..., n}其中0表示仓库1到n表示客户车辆集合K {1, 2, ..., m}m是可用车辆数的上限参数定义q_i 表示客户i的货物需求量仓库的需求量 q_0 0Q 表示车辆载重上限t_ij 表示从节点i到节点j的行驶时间d_ij 表示从节点i到节点j的行驶距离e_i 表示客户i最早允许服务时间仓库的 e_0 0l_i 表示客户i最晚允许服务时间仓库的 l_0 大数s_i 表示客户i的服务时长仓库服务时长为0决策变量有两个x_ijk 是二进制变量等于1时表示车辆k从节点i直接行驶到节点jT_ik 表示车辆k到达节点i的时间。这里要特别注意到达时间和服务开始时间是两个概念。如果车辆早到了服务开始时间等于等待结束后的时间。硬时间窗约束实际约束的是服务开始时间不是到达时间。许多初学建模的朋友在这里犯错把到达时间直接限定在窗口内导致早到情况不可行但事实上早到是被允许的只是需要等待。2.2 目标函数与容量约束的代码化写法目标函数的表达式是minimize M * sum(k, x[0][j][k]) sum(i, j, k, d[i][j] * x[i][j][k])其中sum(k, x[0][j][k]) 对j求和对应实际使用了几辆车。严格来说车辆数目标应该用 y_k sum(j, x[0][j][k]) 这样的辅助变量更清晰但直接用出库弧的和也能达到相同效果因为一辆车要么从仓库出去一次要么不使用。后面代码里我会用辅助变量来表示方便阅读日志。容量约束的表达式是 sum(i, q[i] * x[i][j][k]) Q 对所有k成立表示任意车辆k服务的所有客户需求总和不超过载重上限。流守恒约束由三个表达式组成每个客户恰好被服务一次、车辆从仓库出发、车辆最终回到仓库。本质上是一个标准的不相交路径覆盖问题。2.3 时间窗约束的递推逻辑与消除子环的作用时间窗约束的核心公式是T_jk T_ik s_i t_ij - M * (1 - x_ijk)这个式子很多初学者看不懂整天纠结这个M怎么取。它的物理含义是如果车辆k从i跑到j那么到达j的时间至少等于到达i的时间加上服务时长和行驶时间如果车辆k没有从i跑到j那么x_ijk等于0右侧变成T_ik s_i t_ij - M当M足够大时约束自动松弛不产生任何限制。注意第一个约束只给了下界没给上界。车辆可以提前到达等待但不能晚于时间窗所以还需要一个窗口约束把T限制住每个窗口约束又分两种情况e_i T_ik 表示不早于最早服务时间如果早到就等T表示实际服务开始时间T_ik l_i 表示不晚于最晚服务时间这是硬性的。时间窗约束还有一个非常重要的附属作用消除子环。网络上流传的VRPTW模型通常还要单独加MTZ子环消除约束或DFJ子环消除约束但加入了严格时间窗之后子环往往被自然打破因为子环上的时间递推会产生矛盾。当节点数量适中且时间窗分布合理的时候仅靠时间窗约束就能避免子环出现。不过这个结论有前提时间窗约束必须是强约束车辆到达子环中某个节点后不可能再按原路回到仓库否则仍需额外加子环消除约束。3. 环境准备与测试数据构造没数据你没法调代码3.1 Gurobi安装与学术License配置要点Gurobi在学术和比赛场景下拥有完全免费的学术License学生或教师都可以申请。安装步骤网上有大量教程我简单说几个关键坑第一申请License时填写的邮箱必须是学校邮箱或教育机构认可的邮箱否则审批会被驳回。个人邮箱申请成功的概率很低如果你是参赛选手找指导老师要一个学校邮箱是最稳妥的办法。第二安装完Gurobi之后要运行grbgetkey命令把License文件下载到本地。这一步会在你的用户目录下生成一个gurobi.lic文件。很多人在这一步失败原因是环境变量没有配置好License文件找不到。你可以直接设置环境变量GRB_LICENSE_FILE指向license文件的路径绕开默认路径问题。第三Python调用Gurobi之前一定要确认版本匹配。最新版Gurobi对Python版本有要求如果你的Python是3.13旧版Gurobi根本不支持。建议直接安装最新版Gurobi同时Python不要追新3.10到3.12之间是最稳的选择。你还可以用pip install gurobipy安装Python接口但注意这只是API层底层求解器引擎仍然需要Gurobi安装包和License。不要觉得装个pip包就够了没有License你连一个小规模问题都解不了。3.2 测试数据设计原则与Solomon基准数据VRPTW建模的调试阶段非常推荐先跑通一份小规模随机数据再去挑战标准基准数据。业界公认的测试集是Solomon benchmark包含C系列聚类分布、R系列随机分布、RC系列混合分布规模从25个客户到1000个客户都有。自己构造数据时有几个参数必须相互协调仓库位置、客户坐标、需求量、服务时间、时间窗宽度、车辆载重。如果车辆载重远远大于总需求量容量约束退化成摆设问题实际上变成纯TSPTW这就不符合VRPTW的训练目的。时间窗宽度如果设置得过宽约束几乎不产生作用模型求解难度也大幅下降失去了测试意义。我建议按以下逻辑构造一个26个节点的测试实例1个仓库在坐标(50, 50)25个客户随机分布在100×100的方格内需求量在10到30之间车辆载重设为200时间窗按距离仓库远近来设计避免所有窗口重叠。这样既保证可行解存在又保证约束不是空气。3.3 Python数据结构的建模选择代码里我建议用面向对象或字典的方式来组织数据不建议用纯列表加下标因为VRPTW的约束表达式依赖大量的节点索引和属性查询字典的键值映射读起来清楚得多。常用设计是一个Node类包含id、x坐标、y坐标、需求、最早时间、最晚时间、服务时长这些字段然后节点列表的第一位固定为仓库客户索引从1开始。距离和行驶时间的计算用欧氏距离函数确保对称性。如果要用真实道路距离就把距离矩阵预先算好存成二维数组。数据准备这一步越细后面写约束就越省事。我的习惯是数据部分单独一个函数把所有的节点信息、参数配置集中管理主函数只负责建模型和求解。4. PythonGurobi核心代码逐段拆解4.1 建模型与变量定义从这一步开始进入核心代码。所有代码我都用完整可运行的方式呈现并加了注释。先做变量定义部分import gurobipy as gp from gurobipy import GRB def build_model(nodes, vehicle_num, capacity, dist_matrix, big_M): model gp.Model(VRPTW) n len(nodes) # x[i][j][k]: 车辆k是否从节点i行驶到节点j x {} for i in range(n): for j in range(n): if i ! j: for k in range(vehicle_num): x[i, j, k] model.addVar(vtypeGRB.BINARY, namefx_{i}_{j}_{k}) # T[i][k]: 车辆k到达节点i的时间 T {} for i in range(n): for k in range(vehicle_num): T[i, k] model.addVar(lb0, ubbig_M, vtypeGRB.CONTINUOUS, namefT_{i}_{k}) # y[k]: 车辆k是否被使用 y [] for k in range(vehicle_num): y.append(model.addVar(vtypeGRB.BINARY, namefy_{k})) model.update() return model, x, T, ybig_M这里取一个保守的大数。我通常直接用所有节点时间窗最大值的两倍加总服务时间如果不确定就取10000只要保证远大于任何可行解中的到达时间即可。注意T的ub也要设成big_M否则Gurobi的变量区间过大数值稳定性会下降。节点索引的对应关系要明确0号是仓库客户从1开始。x的索引三元组必须保证i不等于j否则会出现自环变量自环变量虽然约束下不会取1但白白增加求解空间。4.2 目标函数与车辆数约束def add_objective_and_constraints(model, x, T, y, nodes, vehicle_num, capacity, dist_matrix, big_M, big_M_coef): n len(nodes) # 目标函数先最小化车辆数再最小化总距离 total_distance gp.quicksum(dist_matrix[i][j] * x[i, j, k] for i in range(n) for j in range(n) if i ! j for k in range(vehicle_num)) total_vehicle gp.quicksum(y[k] for k in range(vehicle_num)) model.setObjective(big_M_coef * total_vehicle total_distance, GRB.MINIMIZE)big_M_coef就是我前面说的字典序权重系数。比赛场景大可直接取100000这个值远大于任何合理配送方案的总距离能保证优先压缩车辆数。在实际业务里如果你更看重距离和成本可以把权重调小甚至直接把目标函数改成总距离。接下来是车辆数量一致性约束保证车辆被使用的定义与出库弧一致# 车辆k被使用 车辆k至少从仓库出发一次 for k in range(vehicle_num): model.addConstr( gp.quicksum(x[0, j, k] for j in range(1, n)) y[k], namefvehicle_used_{k} )这个约束是整个模型里最容易被忽略的。没有它y变量和x变量之间没有耦合关系目标函数里的车辆数惩罚项就成了一个可以被任意置0的自由变量模型会优先让所有y为0车辆数目标失效。4.3 流守恒约束与容量约束流守恒约束是VRPTW的骨架分三部分每个客户必须被进入一次、每个客户必须被离开一次、仓库的出弧和入弧数量相同且等于车辆使用数。# 每个客户恰好被服务一次入弧为1 for j in range(1, n): model.addConstr( gp.quicksum(x[i, j, k] for i in range(n) if i ! j for k in range(vehicle_num)) 1, namefvisited_{j} ) # 每个客户恰好离开一次出弧为1 for i in range(1, n): model.addConstr( gp.quicksum(x[i, j, k] for j in range(n) if i ! j for k in range(vehicle_num)) 1, namefdeparture_{i} ) # 仓库的入弧数量等于出弧数量等于车辆使用数 model.addConstr( gp.quicksum(x[i, 0, k] for i in range(1, n) for k in range(vehicle_num)) gp.quicksum(x[0, j, k] for j in range(1, n) for k in range(vehicle_num)), namedepot_flow_balance )注意到仓库的流守恒约束我写的是入弧等于出弧没有显式让它等于总车辆数。因为前面vehicle_used约束已经把出库车辆数和y绑定了这里再强制等于总车辆数会造成冗余约束。Gurobi能处理冗余约束但去掉冗余能让求解日志更干净可读性更高。容量约束按车辆维度聚合一辆车在途中的所有需求总和不能超载# 容量约束 for k in range(vehicle_num): model.addConstr( gp.quicksum(nodes[i].demand * x[i, j, k] for i in range(1, n) for j in range(n) if i ! j for k_ in [k]) capacity, namefcapacity_{k} )这里用nodes[i].demand其中nodes是一个自定义Node对象的列表。代码里受限写法可能稍复杂但逻辑就是把客户i的送货需求乘上所有从i出发的弧变量之和一辆车服务的所有客户需求求和。4.4 时间窗约束最难也最容易出错的部分时间窗约束是整个模型里最微妙的部分。前面提到过核心递推公式是T_jk T_ik s_i t_ij - big_M * (1 - x_ijk)但实现时还有一个前置约束如果车辆k不经过节点i那么T[i][k]应该为0还是任意值这会影响时间递推的有效性。我的做法是加一个耦合约束强制T[i][k]只有在车辆访问i时才非零# T[i][k] 只有在车辆k访问i时才有意义 for i in range(1, n): for k in range(vehicle_num): visit gp.quicksum(x[h, i, k] for h in range(n) if h ! i) model.addConstr(T[i, k] big_M * visit, nameftime_link_{i}_{k})这个约束保证了T[i][k]在没有访问时为0。时间递推约束在此基础上可以放心写# 时间递推约束 for k in range(vehicle_num): for i in range(n): for j in range(n): if i ! j: model.addConstr( T[j, k] T[i, k] nodes[i].service_time travel_time[i][j] - big_M * (1 - x[i, j, k]), nameftime_{i}_{j}_{k} )时间窗约束分两步仓库的时间窗口要闭合客户的时间窗口是硬性上下界# 仓库时间窗所有车辆从仓库出发不早于e[0]回到仓库不晚于l[0] for k in range(vehicle_num): model.addConstr(T[0, k] nodes[0].early, namefdepot_start_{k}) model.addConstr(T[0, k] nodes[0].late, namefdepot_end_{k}) # 客户时间窗 for i in range(1, n): for k in range(vehicle_num): model.addConstr(T[i, k] nodes[i].early, namefwindow_early_{i}_{k}) model.addConstr(T[i, k] nodes[i].late, namefwindow_late_{i}_{k})这里时间变量的语义要统一T[i][k]是服务开始时间。如果车辆到达i的时间早于e_i模型会通过递推约束和窗口约束一起把T[i][k]推到e_i相当于强制等待。这是VRPTW建模的标准处理方式。4.5 求解与结果导出求解部分反而简单def solve_model(model, nodes, x, T, y, vehicle_num): model.optimize() if model.status ! GRB.OPTIMAL: print(未找到最优解状态:, model.status) return None routes [] for k in range(vehicle_num): if y[k].X 0.5: continue route [] current 0 while True: route.append(current) for j in range(len(nodes)): if j current: continue if x[current, j, k].X 0.5: current j break else: break if current 0: break routes.append(route) total_dist sum(dist_matrix[route[i - 1]][route[i]] for route in routes for i in range(1, len(route))) print(路线信息:) for idx, route in enumerate(routes): print(f车辆{idx 1}: {route}) print(总车辆数:, len(routes)) print(总行驶距离:, total_dist) return routes, total_dist路线提取的逻辑先把所有从仓库出发且y取出值为1的车辆标出来然后从仓库0开始沿着x取值为1的弧一步一步走直到回到仓库。这里的while循环每步找一个后继节点时间复杂度和节点数成线性关系规模在数千节点内都没问题。5. 完整的可复现代码一个可以直接跑通的最小demo下面的代码把前面所有片段整合成一个完整体数据用一段可复现的随机实例。你可以直接复制运行然后替换成自己的数据import math import random import numpy as np import gurobipy as gp from gurobipy import GRB class Node: def __init__(self, id, x, y, demand, early, late, service_time): self.id id self.x x self.y y self.demand demand self.early early self.late late self.service_time service_time def generate_demo_data(customer_num25): random.seed(42) nodes [] nodes.append(Node(0, 50, 50, 0, 0, 800, 0)) for i in range(1, customer_num 1): x random.randint(10, 90) y random.randint(10, 90) demand random.randint(10, 30) distance_to_depot math.hypot(x - 50, y - 50) early int(distance_to_depot * 0.8) late int(distance_to_depot * 0.8 120) service random.randint(10, 20) nodes.append(Node(i, x, y, demand, early, late, service)) n len(nodes) dist [[0.0] * n for _ in range(n)] travel [[0.0] * n for _ in range(n)] for i in range(n): for j in range(n): if i ! j: d math.hypot(nodes[i].x - nodes[j].x, nodes[i].y - nodes[j].y) dist[i][j] d travel[i][j] d return nodes, dist, travel def solve_vrptw(nodes, dist, travel, vehicle_num, capacity): n len(nodes) big_M 10000 model gp.Model(VRPTW) model.Params.TimeLimit 120 model.Params.MIPGap 0.01 x {} for i in range(n): for j in range(n): if i ! j: for k in range(vehicle_num): x[i, j, k] model.addVar(vtypeGRB.BINARY, namefx_{i}_{j}_{k}) T {} for i in range(n): for k in range(vehicle_num): T[i, k] model.addVar(lb0, ubbig_M, vtypeGRB.CONTINUOUS, namefT_{i}_{k}) y [] for k in range(vehicle_num): y.append(model.addVar(vtypeGRB.BINARY, namefy_{k})) model.update() total_distance gp.quicksum(dist[i][j] * x[i, j, k] for i in range(n) for j in range(n) if i ! j for k in range(vehicle_num)) total_vehicle gp.quicksum(y[k] for k in range(vehicle_num)) model.setObjective(100000 * total_vehicle total_distance, GRB.MINIMIZE) for k in range(vehicle_num): model.addConstr(gp.quicksum(x[0, j, k] for j in range(1, n)) y[k], namefvehicle_used_{k}) for j in range(1, n): model.addConstr(gp.quicksum(x[i, j, k] for i in range(n) if i ! j for k in range(vehicle_num)) 1, namefvisit_{j}) for i in range(1, n): model.addConstr(gp.quicksum(x[i, j, k] for j in range(n) if i ! j for k in range(vehicle_num)) 1, namefleave_{i}) model.addConstr( gp.quicksum(x[i, 0, k] for i in range(1, n) for k in range(vehicle_num)) gp.quicksum(x[0, j, k] for j in range(1, n) for k in range(vehicle_num)), namedepot_flow_balance ) for k in range(vehicle_num): model.addConstr( gp.quicksum(nodes[i].demand * x[i, j, k] for i in range(1, n) for j in range(n) if i ! j) capacity, namefcapacity_{k} ) for i in range(1, n): for k in range(vehicle_num): visit gp.quicksum(x[h, i, k] for h in range(n) if h ! i) model.addConstr(T[i, k] big_M * visit, nameftime_link_{i}_{k}) for k in range(vehicle_num): for i in range(n): for j in range(n): if i ! j: model.addConstr( T[j, k] T[i, k] nodes[i].service_time travel[i][j] - big_M * (1 - x[i, j, k]), nameftime_seq_{i}_{j}_{k} ) for k in range(vehicle_num): model.addConstr(T[0, k] nodes[0].early) model.addConstr(T[0, k] nodes[0].late) for i in range(1, n): for k in range(vehicle_num): model.addConstr(T[i, k] nodes[i].early) model.addConstr(T[i, k] nodes[i].late) model.optimize() if model.status not in (GRB.OPTIMAL, GRB.TIME_LIMIT, GRB.INTERRUPTED): print(求解失败状态码:, model.status) return None routes [] for k in range(vehicle_num): if y[k].X 0.5: continue route [] current 0 while True: route.append(current) for j in range(n): if j current: continue if x[current, j, k].X 0.5: current j break else: break if current 0: break routes.append(route) total_dist 0 for route in routes: for i in range(1, len(route)): total_dist dist[route[i - 1]][route[i]] print(求解状态:, model.status) print(Gurobi目标值:, model.ObjVal) print(车辆数:, len(routes)) print(总行驶距离:, total_dist) for idx, route in enumerate(routes): print(f车辆{idx 1}路线:, route) return routes, total_dist if __name__ __main__: nodes, dist, travel generate_demo_data(20) routes, distance solve_vrptw(nodes, dist, travel, vehicle_num8, capacity100)这份代码的设计意图可以直接说清楚generate_demo_data生成了20个客户和1个仓库时间窗宽120单位服务时间10到20单位车辆载重100。8辆车肯定够用但具体用几辆是求解器去决定的不需要你手动指定。5.1 运行结果的预期效果和日志判读正常求解完成后你会看到Gurobi的标准输出日志。日志里的关键指标有三个Incumbent当前找到的最好可行解BestBd当前松弛问题的理论下界Gap两者之间的相对差距Gap越小说明解越接近最优我的经验是20个客户规模的VRPTWGurobi默认参数下通常几秒到几十秒就能找到最优解。如果你看到Gap长时间不下降大概率是时间窗约束太紧或者big_M设置过大导致数值问题。此时先检查时间窗宽度把late下调或者early上调问题往往能缓解。日志下方生成的结果包含每辆车的访问序列。比如输出车辆1路线: [0, 5, 12, 3, 0]表示车辆1从仓库出发访问客户5、12、3最后回仓库。你可以用这份路线序列做后续可视化或业务系统对接。5.2 路线可视化用matplotlib验证结果合理性模型求解完最好画个图直观检查路线有没有交叉、是否明显不合理。我平时用matplotlib画节点和路线代码量很短import matplotlib.pyplot as plt def plot_routes(nodes, routes): fig, ax plt.subplots(figsize(8, 8)) depot nodes[0] ax.scatter(depot.x, depot.y, cred, s200, markers, labelDepot) for i in range(1, len(nodes)): ax.scatter(nodes[i].x, nodes[i].y, cblue, s80) ax.annotate(str(i), (nodes[i].x, nodes[i].y), textcoordsoffset points, xytext(5, 5)) colors plt.cm.tab10.colors for idx, route in enumerate(routes): xs [nodes[node].x for node in route] ys [nodes[node].y for node in route] ax.plot(xs, ys, colorcolors[idx % len(colors)], linewidth1.5, alpha0.8, labelfVehicle {idx1}) ax.set_xlabel(X coordinate) ax.set_ylabel(Y coordinate) ax.set_title(VRPTW Routes) ax.grid(True) ax.legend() plt.show()可视化不只是为了好看。我调试模型时经常发现Gurobi求出的路线在数值上可行但实际画出来会出现两条路线在中间交叉后各走各的这种方案往往是时间窗约束过松导致的次优解而不是模型错误。反过来如果客户坐标和时间窗分布有明显相关性路线交叉得非常厉害也要回头检查距离矩阵有没有算错。6. 建模过程中的典型坑这五个错误我踩过不止一次6.1 时间窗语义混淆到达时间还是服务开始时间新手在写时间窗约束时经常直接写到达时间 late这是错的。VRPTW标准定义中时间窗约束的是服务开始时间而车辆实际到达节点的时间可能早于early。如果没有等待机制模型会产生大量不可行解或者提前到达的车辆无法被分配。正确处理方式是T[i][k]定义为服务开始时间递推约束中T[j][k] max(T[i][k] s_i t_ij, e_j)但由于线性模型无法直接表达max就用递推加窗口约束的组合来实现。这种做法我在前面代码中已经呈现不要简化。6.2 big_M系数过大导致数值塌陷有些教材直接建议big_M取一个大到绝对够用的数比如10的8次方。这在理论上没错但在Gurobi里会引发严重的数值问题。求解器内部的容差和预处理机制面对大系数时可能导致约束被错误放松或紧化求解器在分支定界过程中不断产生无意义的分支。我实测的经验是big_M不要超过10000能用5000就不用10000。如果你的节点规模特别大时间跨度特别长可以考虑分段设置big_M比如针对每条弧单独算一个足够大的上界公式为节点i的最晚时间加上服务时间加上t_ij的最大值减去节点j的最早时间再加一个松弛量。6.3 车辆数固定与车辆数可变的建模差异有些业务场景要求必须固定用某几辆车有些场景允许车辆数自由变化。这两种建模方式在代码里差别巨大固定车辆数时目标函数不包含车辆数惩罚项y变量可以完全去掉可变车辆数时必须引入y变量并建立车辆使用与出库弧的一致性约束。最常见的问题是把固定车辆数方案里的约束直接搬到可变车辆数模型里导致车辆数目标失效最后解的路线数量远大于预期或者某些车辆被白开出去一趟只服务一个客户。修改的时候先确认你建的是哪个版本的模型。6.4 求解时间爆炸时的应对策略VRPTW是NP-hard问题客户数量超过50个后精确求解时间会急剧上升。如果比赛场景只有几十分钟求解时间建议做好三件事第一设置时间限制和Gurobi的MIPGap参数例如model.Params.TimeLimit 300、model.Params.MIPGap 0.02这样求解器会在2%的gap内返回一个可行解而不是无限跑下去。第二调整Gurobi的MIPFocus参数。MIPFocus1表示更努力找可行解MIPFocus2表示更努力证明最优性MIPFocus3表示两者平衡。求解卡住时把MIPFocus设成1通常能更快拿到一个可接受的解。第三手动给一个初始可行解。很多业务场景中你手头已经有一套人工规划路线可以把它作为MIP start喂给Gurobi帮助求解器更快剪枝。喂初始解的方法是用变量对象调用x[i, j, k].start 1或者直接构造一个Start对象一次性赋值。6.5 距离矩阵对称性假设的正确性很多VRPTW代码默认距离矩阵对称即dist[i][j] dist[j][i]在欧氏距离下成立。但真实业务中单向道路、单行道、收费不同的双向道路都会导致距离矩阵不对称。如果你在业务应用里直接套用对称假设求出的路线可能在实际行驶中无法执行或成本不一致。检查方法是在你计算距离矩阵的代码里加一个判断当dist[i][j] ! dist[j][i]时打印告警。我习惯于保留这个断言即使当时用的数据是对称的作为防止未来数据源更换导致模型静默出错的一道防线。7. 从这份代码出发你可以扩展的三条路线学会了VRPTW的Gurobi建模你实际上掌握了运筹优化里最核心的MIP建模技能。基于这份代码的框架你可以向三个方向继续深挖。第一是目标函数的业务化改造。把最小化距离改成最小化总成本包括固定出车成本、油耗成本、司机工时成本每辆车每个客户的服务时间都可以折算成成本系数。只需要调整目标函数表达式里的系数即可模型结构完全不用变。第二是加入更多业务约束。比如客户有服务时间优先级某些客户必须由特定类型的车辆服务车辆有最大连续行驶时间限制区域限行导致部分路段不可通行。这些约束大多数都可以通过新增0-1变量和约束表达式的组合来实现底层的流守恒和容量约束骨架保持稳定。第三是求解算法的替代。当客户规模超过200个Gurobi的精确求解会力不从心此时可以先跑一个相对短时间的Gurobi求解拿到一个上界再用这个上界作为初始解去跑ALNS或自适应大邻域搜索。Gurobi在VRPTW这类问题上真正的价值在于小规模示例的精确解验证和较大规模实例的可行解生成而不是终极武器。我建议想深入运筹优化的人一定要走完这条路径先亲手把VRPTW的MIP代码完整跑通理解每一行约束的含义再去看各种变种问题。这个基础打牢之后你会发现无论是CVRP还是带二维装箱约束的路径问题建模思路是相通的差别只在约束表达式的复杂度上。不要一上来就跳去学各种元启发式算法连精确模型都写不出来的话那些改进算子的设计理由你根本体会不到。我建这个demo的时候最深的体会是运筹建模的调试周期非常长很多时候不是语法错误而是模型语义和实际业务语义有微妙偏差。线性规划和整数规划的美妙之处在于约束即沟通语言你把物理世界的规则用数学公式说清楚求解器才能在这个框架里找到人类想不到的方案。这份能力是需要大量重复练习才能内化的希望这篇文章能让你少走几个弯路。
Python+Gurobi实现VRPTW：带时间窗车辆路径问题建模与代码全解析
我入行做运筹优化这些年最常被问到的问题就是VRPTW到底怎么用代码写出来很多朋友看论文公式头头是道一打开Gurobi就不知道变量怎么定义、约束怎么加。这篇内容我就用Python加Gurobi的组合把带时间窗车辆路径规划问题的完整建模过程从头到尾捋一遍从问题定义、模型推导到代码实现、结果判读全部覆盖适合数学建模竞赛选手、物流调度方向的研究生、以及刚接触运筹优化想找一份可跑通Demo的开发者参考。我尽量不写废话每步操作都说明背后的逻辑。项目本身不复杂但很多细节会让初学者卡个半天这些我都会单独拎出来讲。VRPTW本质上解决的是有一群客户每个客户有货物需求和服务时间窗车队从仓库出发每个客户只能被访问一次车辆有容量限制目标是让总行驶距离或总车辆数最小。这类问题在物流配送、外卖调度、维修工单派发里非常常见建模实战的价值远不只是应付比赛而是真正能迁移到业务场景里用的抽象能力。我先从模型定义开始把VRPTW这个问题的边界彻底讲清楚。1. VRPTW到底是什么这个模型在解决哪一层的问题1.1 问题定义与适用场景边界带时间窗的车辆路径问题英文全称是Vehicle Routing Problem with Time Windows标准定义是在一个配送网络里有一座仓库 depot 和若干个客户点 customer 每辆运输车有载重上限每个客户有确定的货物需求量并且客户只在特定时间段内接受服务。你要规划一组车辆路线让所有客户都被服务到同时满足车辆容量约束、时间窗约束和道路连通性约束。这句话拆开看每一部分都是约束的来源。车辆容量约束意味着一条路线上所有客户的需求总和不能超过单车载重时间窗约束意味着车辆到达客户的时间不能早于最早允许服务时间也不能晚于最晚允许服务时间如果到早了必须等待道路连通性约束则保证每条路线是从仓库出发最终回到仓库的完整回路。不是所有配送问题都是VRPTW。很多新手会把VRPTW和TSP、VRP、MDVRP搞混。TSP是单个销售员走遍所有城市没有容量和时间窗概念VRP是VRPTW去掉时间窗的版本MDVRP是多个仓库的版本。VRPTW在最基础的VRP之上叠加了时间维度的约束这使得问题的可行域被大幅压缩求解难度也明显上升。你需要先判断自己的业务场景是否真的存在时间约束如果只是单纯安排路线且没有客户时间要求直接建模成VRP会更高效。1.2 惩罚成本模型和硬时间窗模型的取舍VRPTW有两种常见建模方式硬时间窗和软时间窗。硬时间窗要求车辆到达时间必须在时间窗内早到只能等待晚到直接不可行软时间窗允许违反时间窗但产生惩罚成本。从数学建模和Gurobi求解的角度看硬时间窗是纯MIP约束约束表达式为等待时间加行驶时间的递推关系软时间窗则需要引入非负偏差变量在目标函数中新增惩罚项。两者在Gurobi里的实现差别很大。我个人的建议是比赛场景优先用硬时间窗。原因有两点第一硬时间窗模型更干净不需要调惩罚系数惩罚系数一旦设置不当结果可能为了省成本大面积超时业务上不可接受第二Gurobi处理硬时间窗约束时线性松弛的质量通常更好分支定界效率更高。业务落地场景则用软时间窗更实际因为真实道路状况下硬性时间窗往往过于理想化但那是另一个话题本文先用硬时间窗。1.3 目标函数的设计逻辑VRPTW最常见的双目标是最小化车辆数和最小化总行驶距离。实际操作中大家通常采用字典序优化先把车辆数压到最小再在最少车辆的前提下缩短总距离。字典序优化的做法是给两个目标设置不同的权重系数比如最小化M * 车辆数 总距离其中M取一个非常大的数保证优先优化车辆数。这个M怎么取有讲究它必须大于所有可行方案下总距离的最大可能差值一般取总距离上界的估计值即可比如所有客户到仓库往返距离之和。如果M取太小模型可能为了省几公里路程多派一辆车结果得到的解与预期完全不符。2. 数学模型精读约束公式背后的物理含义与表达陷阱2.1 参数、集合与决策变量的完整定义在写代码之前必须先把数学模型用标准MIP形式表达清楚。数学符号虽然啰嗦但它是调试代码的唯一依据。集合定义节点集合N {0, 1, 2, ..., n}其中0表示仓库1到n表示客户车辆集合K {1, 2, ..., m}m是可用车辆数的上限参数定义q_i 表示客户i的货物需求量仓库的需求量 q_0 0Q 表示车辆载重上限t_ij 表示从节点i到节点j的行驶时间d_ij 表示从节点i到节点j的行驶距离e_i 表示客户i最早允许服务时间仓库的 e_0 0l_i 表示客户i最晚允许服务时间仓库的 l_0 大数s_i 表示客户i的服务时长仓库服务时长为0决策变量有两个x_ijk 是二进制变量等于1时表示车辆k从节点i直接行驶到节点jT_ik 表示车辆k到达节点i的时间。这里要特别注意到达时间和服务开始时间是两个概念。如果车辆早到了服务开始时间等于等待结束后的时间。硬时间窗约束实际约束的是服务开始时间不是到达时间。许多初学建模的朋友在这里犯错把到达时间直接限定在窗口内导致早到情况不可行但事实上早到是被允许的只是需要等待。2.2 目标函数与容量约束的代码化写法目标函数的表达式是minimize M * sum(k, x[0][j][k]) sum(i, j, k, d[i][j] * x[i][j][k])其中sum(k, x[0][j][k]) 对j求和对应实际使用了几辆车。严格来说车辆数目标应该用 y_k sum(j, x[0][j][k]) 这样的辅助变量更清晰但直接用出库弧的和也能达到相同效果因为一辆车要么从仓库出去一次要么不使用。后面代码里我会用辅助变量来表示方便阅读日志。容量约束的表达式是 sum(i, q[i] * x[i][j][k]) Q 对所有k成立表示任意车辆k服务的所有客户需求总和不超过载重上限。流守恒约束由三个表达式组成每个客户恰好被服务一次、车辆从仓库出发、车辆最终回到仓库。本质上是一个标准的不相交路径覆盖问题。2.3 时间窗约束的递推逻辑与消除子环的作用时间窗约束的核心公式是T_jk T_ik s_i t_ij - M * (1 - x_ijk)这个式子很多初学者看不懂整天纠结这个M怎么取。它的物理含义是如果车辆k从i跑到j那么到达j的时间至少等于到达i的时间加上服务时长和行驶时间如果车辆k没有从i跑到j那么x_ijk等于0右侧变成T_ik s_i t_ij - M当M足够大时约束自动松弛不产生任何限制。注意第一个约束只给了下界没给上界。车辆可以提前到达等待但不能晚于时间窗所以还需要一个窗口约束把T限制住每个窗口约束又分两种情况e_i T_ik 表示不早于最早服务时间如果早到就等T表示实际服务开始时间T_ik l_i 表示不晚于最晚服务时间这是硬性的。时间窗约束还有一个非常重要的附属作用消除子环。网络上流传的VRPTW模型通常还要单独加MTZ子环消除约束或DFJ子环消除约束但加入了严格时间窗之后子环往往被自然打破因为子环上的时间递推会产生矛盾。当节点数量适中且时间窗分布合理的时候仅靠时间窗约束就能避免子环出现。不过这个结论有前提时间窗约束必须是强约束车辆到达子环中某个节点后不可能再按原路回到仓库否则仍需额外加子环消除约束。3. 环境准备与测试数据构造没数据你没法调代码3.1 Gurobi安装与学术License配置要点Gurobi在学术和比赛场景下拥有完全免费的学术License学生或教师都可以申请。安装步骤网上有大量教程我简单说几个关键坑第一申请License时填写的邮箱必须是学校邮箱或教育机构认可的邮箱否则审批会被驳回。个人邮箱申请成功的概率很低如果你是参赛选手找指导老师要一个学校邮箱是最稳妥的办法。第二安装完Gurobi之后要运行grbgetkey命令把License文件下载到本地。这一步会在你的用户目录下生成一个gurobi.lic文件。很多人在这一步失败原因是环境变量没有配置好License文件找不到。你可以直接设置环境变量GRB_LICENSE_FILE指向license文件的路径绕开默认路径问题。第三Python调用Gurobi之前一定要确认版本匹配。最新版Gurobi对Python版本有要求如果你的Python是3.13旧版Gurobi根本不支持。建议直接安装最新版Gurobi同时Python不要追新3.10到3.12之间是最稳的选择。你还可以用pip install gurobipy安装Python接口但注意这只是API层底层求解器引擎仍然需要Gurobi安装包和License。不要觉得装个pip包就够了没有License你连一个小规模问题都解不了。3.2 测试数据设计原则与Solomon基准数据VRPTW建模的调试阶段非常推荐先跑通一份小规模随机数据再去挑战标准基准数据。业界公认的测试集是Solomon benchmark包含C系列聚类分布、R系列随机分布、RC系列混合分布规模从25个客户到1000个客户都有。自己构造数据时有几个参数必须相互协调仓库位置、客户坐标、需求量、服务时间、时间窗宽度、车辆载重。如果车辆载重远远大于总需求量容量约束退化成摆设问题实际上变成纯TSPTW这就不符合VRPTW的训练目的。时间窗宽度如果设置得过宽约束几乎不产生作用模型求解难度也大幅下降失去了测试意义。我建议按以下逻辑构造一个26个节点的测试实例1个仓库在坐标(50, 50)25个客户随机分布在100×100的方格内需求量在10到30之间车辆载重设为200时间窗按距离仓库远近来设计避免所有窗口重叠。这样既保证可行解存在又保证约束不是空气。3.3 Python数据结构的建模选择代码里我建议用面向对象或字典的方式来组织数据不建议用纯列表加下标因为VRPTW的约束表达式依赖大量的节点索引和属性查询字典的键值映射读起来清楚得多。常用设计是一个Node类包含id、x坐标、y坐标、需求、最早时间、最晚时间、服务时长这些字段然后节点列表的第一位固定为仓库客户索引从1开始。距离和行驶时间的计算用欧氏距离函数确保对称性。如果要用真实道路距离就把距离矩阵预先算好存成二维数组。数据准备这一步越细后面写约束就越省事。我的习惯是数据部分单独一个函数把所有的节点信息、参数配置集中管理主函数只负责建模型和求解。4. PythonGurobi核心代码逐段拆解4.1 建模型与变量定义从这一步开始进入核心代码。所有代码我都用完整可运行的方式呈现并加了注释。先做变量定义部分import gurobipy as gp from gurobipy import GRB def build_model(nodes, vehicle_num, capacity, dist_matrix, big_M): model gp.Model(VRPTW) n len(nodes) # x[i][j][k]: 车辆k是否从节点i行驶到节点j x {} for i in range(n): for j in range(n): if i ! j: for k in range(vehicle_num): x[i, j, k] model.addVar(vtypeGRB.BINARY, namefx_{i}_{j}_{k}) # T[i][k]: 车辆k到达节点i的时间 T {} for i in range(n): for k in range(vehicle_num): T[i, k] model.addVar(lb0, ubbig_M, vtypeGRB.CONTINUOUS, namefT_{i}_{k}) # y[k]: 车辆k是否被使用 y [] for k in range(vehicle_num): y.append(model.addVar(vtypeGRB.BINARY, namefy_{k})) model.update() return model, x, T, ybig_M这里取一个保守的大数。我通常直接用所有节点时间窗最大值的两倍加总服务时间如果不确定就取10000只要保证远大于任何可行解中的到达时间即可。注意T的ub也要设成big_M否则Gurobi的变量区间过大数值稳定性会下降。节点索引的对应关系要明确0号是仓库客户从1开始。x的索引三元组必须保证i不等于j否则会出现自环变量自环变量虽然约束下不会取1但白白增加求解空间。4.2 目标函数与车辆数约束def add_objective_and_constraints(model, x, T, y, nodes, vehicle_num, capacity, dist_matrix, big_M, big_M_coef): n len(nodes) # 目标函数先最小化车辆数再最小化总距离 total_distance gp.quicksum(dist_matrix[i][j] * x[i, j, k] for i in range(n) for j in range(n) if i ! j for k in range(vehicle_num)) total_vehicle gp.quicksum(y[k] for k in range(vehicle_num)) model.setObjective(big_M_coef * total_vehicle total_distance, GRB.MINIMIZE)big_M_coef就是我前面说的字典序权重系数。比赛场景大可直接取100000这个值远大于任何合理配送方案的总距离能保证优先压缩车辆数。在实际业务里如果你更看重距离和成本可以把权重调小甚至直接把目标函数改成总距离。接下来是车辆数量一致性约束保证车辆被使用的定义与出库弧一致# 车辆k被使用 车辆k至少从仓库出发一次 for k in range(vehicle_num): model.addConstr( gp.quicksum(x[0, j, k] for j in range(1, n)) y[k], namefvehicle_used_{k} )这个约束是整个模型里最容易被忽略的。没有它y变量和x变量之间没有耦合关系目标函数里的车辆数惩罚项就成了一个可以被任意置0的自由变量模型会优先让所有y为0车辆数目标失效。4.3 流守恒约束与容量约束流守恒约束是VRPTW的骨架分三部分每个客户必须被进入一次、每个客户必须被离开一次、仓库的出弧和入弧数量相同且等于车辆使用数。# 每个客户恰好被服务一次入弧为1 for j in range(1, n): model.addConstr( gp.quicksum(x[i, j, k] for i in range(n) if i ! j for k in range(vehicle_num)) 1, namefvisited_{j} ) # 每个客户恰好离开一次出弧为1 for i in range(1, n): model.addConstr( gp.quicksum(x[i, j, k] for j in range(n) if i ! j for k in range(vehicle_num)) 1, namefdeparture_{i} ) # 仓库的入弧数量等于出弧数量等于车辆使用数 model.addConstr( gp.quicksum(x[i, 0, k] for i in range(1, n) for k in range(vehicle_num)) gp.quicksum(x[0, j, k] for j in range(1, n) for k in range(vehicle_num)), namedepot_flow_balance )注意到仓库的流守恒约束我写的是入弧等于出弧没有显式让它等于总车辆数。因为前面vehicle_used约束已经把出库车辆数和y绑定了这里再强制等于总车辆数会造成冗余约束。Gurobi能处理冗余约束但去掉冗余能让求解日志更干净可读性更高。容量约束按车辆维度聚合一辆车在途中的所有需求总和不能超载# 容量约束 for k in range(vehicle_num): model.addConstr( gp.quicksum(nodes[i].demand * x[i, j, k] for i in range(1, n) for j in range(n) if i ! j for k_ in [k]) capacity, namefcapacity_{k} )这里用nodes[i].demand其中nodes是一个自定义Node对象的列表。代码里受限写法可能稍复杂但逻辑就是把客户i的送货需求乘上所有从i出发的弧变量之和一辆车服务的所有客户需求求和。4.4 时间窗约束最难也最容易出错的部分时间窗约束是整个模型里最微妙的部分。前面提到过核心递推公式是T_jk T_ik s_i t_ij - big_M * (1 - x_ijk)但实现时还有一个前置约束如果车辆k不经过节点i那么T[i][k]应该为0还是任意值这会影响时间递推的有效性。我的做法是加一个耦合约束强制T[i][k]只有在车辆访问i时才非零# T[i][k] 只有在车辆k访问i时才有意义 for i in range(1, n): for k in range(vehicle_num): visit gp.quicksum(x[h, i, k] for h in range(n) if h ! i) model.addConstr(T[i, k] big_M * visit, nameftime_link_{i}_{k})这个约束保证了T[i][k]在没有访问时为0。时间递推约束在此基础上可以放心写# 时间递推约束 for k in range(vehicle_num): for i in range(n): for j in range(n): if i ! j: model.addConstr( T[j, k] T[i, k] nodes[i].service_time travel_time[i][j] - big_M * (1 - x[i, j, k]), nameftime_{i}_{j}_{k} )时间窗约束分两步仓库的时间窗口要闭合客户的时间窗口是硬性上下界# 仓库时间窗所有车辆从仓库出发不早于e[0]回到仓库不晚于l[0] for k in range(vehicle_num): model.addConstr(T[0, k] nodes[0].early, namefdepot_start_{k}) model.addConstr(T[0, k] nodes[0].late, namefdepot_end_{k}) # 客户时间窗 for i in range(1, n): for k in range(vehicle_num): model.addConstr(T[i, k] nodes[i].early, namefwindow_early_{i}_{k}) model.addConstr(T[i, k] nodes[i].late, namefwindow_late_{i}_{k})这里时间变量的语义要统一T[i][k]是服务开始时间。如果车辆到达i的时间早于e_i模型会通过递推约束和窗口约束一起把T[i][k]推到e_i相当于强制等待。这是VRPTW建模的标准处理方式。4.5 求解与结果导出求解部分反而简单def solve_model(model, nodes, x, T, y, vehicle_num): model.optimize() if model.status ! GRB.OPTIMAL: print(未找到最优解状态:, model.status) return None routes [] for k in range(vehicle_num): if y[k].X 0.5: continue route [] current 0 while True: route.append(current) for j in range(len(nodes)): if j current: continue if x[current, j, k].X 0.5: current j break else: break if current 0: break routes.append(route) total_dist sum(dist_matrix[route[i - 1]][route[i]] for route in routes for i in range(1, len(route))) print(路线信息:) for idx, route in enumerate(routes): print(f车辆{idx 1}: {route}) print(总车辆数:, len(routes)) print(总行驶距离:, total_dist) return routes, total_dist路线提取的逻辑先把所有从仓库出发且y取出值为1的车辆标出来然后从仓库0开始沿着x取值为1的弧一步一步走直到回到仓库。这里的while循环每步找一个后继节点时间复杂度和节点数成线性关系规模在数千节点内都没问题。5. 完整的可复现代码一个可以直接跑通的最小demo下面的代码把前面所有片段整合成一个完整体数据用一段可复现的随机实例。你可以直接复制运行然后替换成自己的数据import math import random import numpy as np import gurobipy as gp from gurobipy import GRB class Node: def __init__(self, id, x, y, demand, early, late, service_time): self.id id self.x x self.y y self.demand demand self.early early self.late late self.service_time service_time def generate_demo_data(customer_num25): random.seed(42) nodes [] nodes.append(Node(0, 50, 50, 0, 0, 800, 0)) for i in range(1, customer_num 1): x random.randint(10, 90) y random.randint(10, 90) demand random.randint(10, 30) distance_to_depot math.hypot(x - 50, y - 50) early int(distance_to_depot * 0.8) late int(distance_to_depot * 0.8 120) service random.randint(10, 20) nodes.append(Node(i, x, y, demand, early, late, service)) n len(nodes) dist [[0.0] * n for _ in range(n)] travel [[0.0] * n for _ in range(n)] for i in range(n): for j in range(n): if i ! j: d math.hypot(nodes[i].x - nodes[j].x, nodes[i].y - nodes[j].y) dist[i][j] d travel[i][j] d return nodes, dist, travel def solve_vrptw(nodes, dist, travel, vehicle_num, capacity): n len(nodes) big_M 10000 model gp.Model(VRPTW) model.Params.TimeLimit 120 model.Params.MIPGap 0.01 x {} for i in range(n): for j in range(n): if i ! j: for k in range(vehicle_num): x[i, j, k] model.addVar(vtypeGRB.BINARY, namefx_{i}_{j}_{k}) T {} for i in range(n): for k in range(vehicle_num): T[i, k] model.addVar(lb0, ubbig_M, vtypeGRB.CONTINUOUS, namefT_{i}_{k}) y [] for k in range(vehicle_num): y.append(model.addVar(vtypeGRB.BINARY, namefy_{k})) model.update() total_distance gp.quicksum(dist[i][j] * x[i, j, k] for i in range(n) for j in range(n) if i ! j for k in range(vehicle_num)) total_vehicle gp.quicksum(y[k] for k in range(vehicle_num)) model.setObjective(100000 * total_vehicle total_distance, GRB.MINIMIZE) for k in range(vehicle_num): model.addConstr(gp.quicksum(x[0, j, k] for j in range(1, n)) y[k], namefvehicle_used_{k}) for j in range(1, n): model.addConstr(gp.quicksum(x[i, j, k] for i in range(n) if i ! j for k in range(vehicle_num)) 1, namefvisit_{j}) for i in range(1, n): model.addConstr(gp.quicksum(x[i, j, k] for j in range(n) if i ! j for k in range(vehicle_num)) 1, namefleave_{i}) model.addConstr( gp.quicksum(x[i, 0, k] for i in range(1, n) for k in range(vehicle_num)) gp.quicksum(x[0, j, k] for j in range(1, n) for k in range(vehicle_num)), namedepot_flow_balance ) for k in range(vehicle_num): model.addConstr( gp.quicksum(nodes[i].demand * x[i, j, k] for i in range(1, n) for j in range(n) if i ! j) capacity, namefcapacity_{k} ) for i in range(1, n): for k in range(vehicle_num): visit gp.quicksum(x[h, i, k] for h in range(n) if h ! i) model.addConstr(T[i, k] big_M * visit, nameftime_link_{i}_{k}) for k in range(vehicle_num): for i in range(n): for j in range(n): if i ! j: model.addConstr( T[j, k] T[i, k] nodes[i].service_time travel[i][j] - big_M * (1 - x[i, j, k]), nameftime_seq_{i}_{j}_{k} ) for k in range(vehicle_num): model.addConstr(T[0, k] nodes[0].early) model.addConstr(T[0, k] nodes[0].late) for i in range(1, n): for k in range(vehicle_num): model.addConstr(T[i, k] nodes[i].early) model.addConstr(T[i, k] nodes[i].late) model.optimize() if model.status not in (GRB.OPTIMAL, GRB.TIME_LIMIT, GRB.INTERRUPTED): print(求解失败状态码:, model.status) return None routes [] for k in range(vehicle_num): if y[k].X 0.5: continue route [] current 0 while True: route.append(current) for j in range(n): if j current: continue if x[current, j, k].X 0.5: current j break else: break if current 0: break routes.append(route) total_dist 0 for route in routes: for i in range(1, len(route)): total_dist dist[route[i - 1]][route[i]] print(求解状态:, model.status) print(Gurobi目标值:, model.ObjVal) print(车辆数:, len(routes)) print(总行驶距离:, total_dist) for idx, route in enumerate(routes): print(f车辆{idx 1}路线:, route) return routes, total_dist if __name__ __main__: nodes, dist, travel generate_demo_data(20) routes, distance solve_vrptw(nodes, dist, travel, vehicle_num8, capacity100)这份代码的设计意图可以直接说清楚generate_demo_data生成了20个客户和1个仓库时间窗宽120单位服务时间10到20单位车辆载重100。8辆车肯定够用但具体用几辆是求解器去决定的不需要你手动指定。5.1 运行结果的预期效果和日志判读正常求解完成后你会看到Gurobi的标准输出日志。日志里的关键指标有三个Incumbent当前找到的最好可行解BestBd当前松弛问题的理论下界Gap两者之间的相对差距Gap越小说明解越接近最优我的经验是20个客户规模的VRPTWGurobi默认参数下通常几秒到几十秒就能找到最优解。如果你看到Gap长时间不下降大概率是时间窗约束太紧或者big_M设置过大导致数值问题。此时先检查时间窗宽度把late下调或者early上调问题往往能缓解。日志下方生成的结果包含每辆车的访问序列。比如输出车辆1路线: [0, 5, 12, 3, 0]表示车辆1从仓库出发访问客户5、12、3最后回仓库。你可以用这份路线序列做后续可视化或业务系统对接。5.2 路线可视化用matplotlib验证结果合理性模型求解完最好画个图直观检查路线有没有交叉、是否明显不合理。我平时用matplotlib画节点和路线代码量很短import matplotlib.pyplot as plt def plot_routes(nodes, routes): fig, ax plt.subplots(figsize(8, 8)) depot nodes[0] ax.scatter(depot.x, depot.y, cred, s200, markers, labelDepot) for i in range(1, len(nodes)): ax.scatter(nodes[i].x, nodes[i].y, cblue, s80) ax.annotate(str(i), (nodes[i].x, nodes[i].y), textcoordsoffset points, xytext(5, 5)) colors plt.cm.tab10.colors for idx, route in enumerate(routes): xs [nodes[node].x for node in route] ys [nodes[node].y for node in route] ax.plot(xs, ys, colorcolors[idx % len(colors)], linewidth1.5, alpha0.8, labelfVehicle {idx1}) ax.set_xlabel(X coordinate) ax.set_ylabel(Y coordinate) ax.set_title(VRPTW Routes) ax.grid(True) ax.legend() plt.show()可视化不只是为了好看。我调试模型时经常发现Gurobi求出的路线在数值上可行但实际画出来会出现两条路线在中间交叉后各走各的这种方案往往是时间窗约束过松导致的次优解而不是模型错误。反过来如果客户坐标和时间窗分布有明显相关性路线交叉得非常厉害也要回头检查距离矩阵有没有算错。6. 建模过程中的典型坑这五个错误我踩过不止一次6.1 时间窗语义混淆到达时间还是服务开始时间新手在写时间窗约束时经常直接写到达时间 late这是错的。VRPTW标准定义中时间窗约束的是服务开始时间而车辆实际到达节点的时间可能早于early。如果没有等待机制模型会产生大量不可行解或者提前到达的车辆无法被分配。正确处理方式是T[i][k]定义为服务开始时间递推约束中T[j][k] max(T[i][k] s_i t_ij, e_j)但由于线性模型无法直接表达max就用递推加窗口约束的组合来实现。这种做法我在前面代码中已经呈现不要简化。6.2 big_M系数过大导致数值塌陷有些教材直接建议big_M取一个大到绝对够用的数比如10的8次方。这在理论上没错但在Gurobi里会引发严重的数值问题。求解器内部的容差和预处理机制面对大系数时可能导致约束被错误放松或紧化求解器在分支定界过程中不断产生无意义的分支。我实测的经验是big_M不要超过10000能用5000就不用10000。如果你的节点规模特别大时间跨度特别长可以考虑分段设置big_M比如针对每条弧单独算一个足够大的上界公式为节点i的最晚时间加上服务时间加上t_ij的最大值减去节点j的最早时间再加一个松弛量。6.3 车辆数固定与车辆数可变的建模差异有些业务场景要求必须固定用某几辆车有些场景允许车辆数自由变化。这两种建模方式在代码里差别巨大固定车辆数时目标函数不包含车辆数惩罚项y变量可以完全去掉可变车辆数时必须引入y变量并建立车辆使用与出库弧的一致性约束。最常见的问题是把固定车辆数方案里的约束直接搬到可变车辆数模型里导致车辆数目标失效最后解的路线数量远大于预期或者某些车辆被白开出去一趟只服务一个客户。修改的时候先确认你建的是哪个版本的模型。6.4 求解时间爆炸时的应对策略VRPTW是NP-hard问题客户数量超过50个后精确求解时间会急剧上升。如果比赛场景只有几十分钟求解时间建议做好三件事第一设置时间限制和Gurobi的MIPGap参数例如model.Params.TimeLimit 300、model.Params.MIPGap 0.02这样求解器会在2%的gap内返回一个可行解而不是无限跑下去。第二调整Gurobi的MIPFocus参数。MIPFocus1表示更努力找可行解MIPFocus2表示更努力证明最优性MIPFocus3表示两者平衡。求解卡住时把MIPFocus设成1通常能更快拿到一个可接受的解。第三手动给一个初始可行解。很多业务场景中你手头已经有一套人工规划路线可以把它作为MIP start喂给Gurobi帮助求解器更快剪枝。喂初始解的方法是用变量对象调用x[i, j, k].start 1或者直接构造一个Start对象一次性赋值。6.5 距离矩阵对称性假设的正确性很多VRPTW代码默认距离矩阵对称即dist[i][j] dist[j][i]在欧氏距离下成立。但真实业务中单向道路、单行道、收费不同的双向道路都会导致距离矩阵不对称。如果你在业务应用里直接套用对称假设求出的路线可能在实际行驶中无法执行或成本不一致。检查方法是在你计算距离矩阵的代码里加一个判断当dist[i][j] ! dist[j][i]时打印告警。我习惯于保留这个断言即使当时用的数据是对称的作为防止未来数据源更换导致模型静默出错的一道防线。7. 从这份代码出发你可以扩展的三条路线学会了VRPTW的Gurobi建模你实际上掌握了运筹优化里最核心的MIP建模技能。基于这份代码的框架你可以向三个方向继续深挖。第一是目标函数的业务化改造。把最小化距离改成最小化总成本包括固定出车成本、油耗成本、司机工时成本每辆车每个客户的服务时间都可以折算成成本系数。只需要调整目标函数表达式里的系数即可模型结构完全不用变。第二是加入更多业务约束。比如客户有服务时间优先级某些客户必须由特定类型的车辆服务车辆有最大连续行驶时间限制区域限行导致部分路段不可通行。这些约束大多数都可以通过新增0-1变量和约束表达式的组合来实现底层的流守恒和容量约束骨架保持稳定。第三是求解算法的替代。当客户规模超过200个Gurobi的精确求解会力不从心此时可以先跑一个相对短时间的Gurobi求解拿到一个上界再用这个上界作为初始解去跑ALNS或自适应大邻域搜索。Gurobi在VRPTW这类问题上真正的价值在于小规模示例的精确解验证和较大规模实例的可行解生成而不是终极武器。我建议想深入运筹优化的人一定要走完这条路径先亲手把VRPTW的MIP代码完整跑通理解每一行约束的含义再去看各种变种问题。这个基础打牢之后你会发现无论是CVRP还是带二维装箱约束的路径问题建模思路是相通的差别只在约束表达式的复杂度上。不要一上来就跳去学各种元启发式算法连精确模型都写不出来的话那些改进算子的设计理由你根本体会不到。我建这个demo的时候最深的体会是运筹建模的调试周期非常长很多时候不是语法错误而是模型语义和实际业务语义有微妙偏差。线性规划和整数规划的美妙之处在于约束即沟通语言你把物理世界的规则用数学公式说清楚求解器才能在这个框架里找到人类想不到的方案。这份能力是需要大量重复练习才能内化的希望这篇文章能让你少走几个弯路。
 管道化去首元素变换)
