尧图网页设计 
FSTSPThe Flying Sidekick Traveling Salesman Problem是这几年物流优化领域特别火的一个组合优化问题一辆卡车带着一架无人机出发无人机可以从卡车上起飞去服务某一个客户然后在后续某个节点再和卡车汇合。听起来是不是有点像外卖小哥和无人机配合送货但真正把这个问题建模成数学规划再用Gurobi求到最优解里面有不少门道。我花了两周多的时间把这个论文模型完整复现了一遍实现了Python Gurobi的求解代码也踩了不少坑。这篇文章就把完整思路、模型细节、代码实现以及我实测中遇到的问题全部整理出来希望能帮你少走弯路。整个项目核心是复现论文里的MILP混合整数线性规划模型并把它跑通、跑对。适合对运筹优化、路径规划感兴趣的研究生、算法工程师也适合正在学习Gurobi但苦于没有完整项目练手的人。下面直接进入正题。1. FSTSP问题背景与整体建模思路1.1 为什么是卡车无人机而不是纯无人机配送先聊点题外话为什么FSTSP这类问题在学术界和工业界都这么受关注。纯卡车配送的问题本质上是经典TSP或VRP的变种算法研究已经非常成熟。但纯卡车的短板很明显最后一公里路况复杂、速度受限而且一旦遇到交通拥堵整个配送时效就崩了。纯无人机配送呢虽然直线飞行不堵车但载重限制、续航限制非常苛刻一架无人机最多也就送一两公斤的包裹飞个几十分钟就要换电池根本扛不住大批量配送。卡车无人机的组合模式就聪明在优势互补卡车负责把无人机带到接近客户的位置扩大无人机的覆盖半径同时卡车自己做常规的送货任务无人机负责从卡车上起飞去服务那些卡车不方便去、或者绕路成本很高的客户然后在下游汇合点重新回到卡车上。无人机不需要飞完全程只需要飞一个短途的“支线”既绕开了交通拥堵又解决了续航焦虑。FSTSP研究的就是在这个模式下怎样规划卡车路径和无人机路径让整体完成时间最短。1.2 问题定义与假设条件FSTSP有一个标准的设定这个设定是论文里明确给出的复现时必须严格遵守。一个配送中心depot里有一辆卡车和一架无人机需要服务若干个客户节点。卡车可以携带无人机一起移动无人机也可以脱离卡车独立飞行去服务客户。几个关键假设如下每个客户节点必须被服务一次且只能被服务一次要么卡车去服务要么无人机去服务。无人机从卡车上的某个节点起飞服务一个客户后在后续某个节点与卡车汇合重新装载到卡车上。一台无人机同时最多只能执行一个任务。无人机有续航限制和载重限制论文里通常以最大飞行距离来体现续航约束。卡车速度通常慢于无人机速度这里指空中的直线速度实际场景要考虑道路系数。目标是让卡车完成所有服务任务并返回配送中心的时间makespan最小化。严格来说因为无人机最终要回到卡车上所以目标本质是卡车的总完成时间。我复现的版本客户规模不大10到15个客户节点量级正好是精确求解算法能够处理的规模。因为FSTSP是NP-hard问题Gurobi这类精确求解器在大规模问题上会指数级爆炸所以论文里通常也是用中小规模实例做验证。如果追求更大规模就得用启发式或者分支定价之类的进阶算法了但那是另一个话题。1.3 建模方案选型MILP而不是其他方法FSTSP的建模方式有好几种常见的有基于时间和空间的离散时间模型、基于事件的状态转移模型还有最经典的MILP模型。为什么最终选了MILP核心原因是Gurobi 对 MILP 的支持极其成熟求解器内置了分支定界、割平面、启发式等多种策略不需要我自己写求解算法只需要把问题表达成正确的线性约束和目标函数即可。另一个重要原因是FSTSP本质上是一个路径规划资源调度的混合问题里面有逻辑判断哪个节点由谁服务、先后顺序是什么、有时间顺序卡车和无人机的时间先后关系、有资源约束无人机一次只能执行一个任务。这些问题通过0-1整数变量和连续时间变量建模成MILP是最自然的。比纯启发式算法好在只要Gurobi在有限时间内找到了全局最优解你就可以100%确定这个结果是理论最优的这对于复现论文、验证对比实验至关重要。我当时也考虑过要不要用动态规划或者A*后来觉得没必要。精确解发出的“下限”是后续所有近似算法对比的基准论文实验部分最需要的就是这个基准值。1.4 代码整体架构设计代码的模块划分我分成了五个部分数据生成模块生成客户节点的坐标、服务时间、卡车速度、无人机速度、无人机最大飞行距离等参数。距离计算模块计算卡车路网距离用欧氏距离乘道路系数模拟、无人机直线距离。模型构建模块定义决策变量、约束条件、目标函数这是核心中的核心。求解与输出模块调用Gurobi求解输出最优目标值、卡车路径、无人机起飞/汇合节点等。可视化模块用Matplotlib画卡车和无人机的路径图直观验证结果合理性。这个架构看起来很常规但在代码复现里有特别的讲究。比如距离计算模块独立出来是为了方便切换不同的距离度量方式模型构建模块和求解模块分开是为了后续调参方便。如果一开始把这些都揉在一起后面调试每个约束条件时会非常痛苦。2. 核心决策变量与关键约束拆解模型部分是整个复现的灵魂。我在写代码之前花了很多时间把论文里的每个符号、每个约束都手推了一遍。这里就把最终版本的核心内容完整拆开讲。2.1 参数定义与输入数据先定义基本集合和参数。我用 (N) 表示客户节点集合(N_0) 表示包含配送中心在内的全部节点配送中心编号0客户节点编号1到n。配送中心是卡车的起点和终点无人机也从这个点开始装载。具体参数如下表参数含义n客户节点数量V_s卡车行驶速度V_d无人机飞行速度通常设成比卡车快比如卡车的1.5倍e无人机服务一个客户之后的装卸时间起飞/降落时间D_max无人机最大飞行里程(x_i, y_i)节点i的坐标d_ij卡车从节点i到节点j的距离带道路系数的欧氏距离dd_ij无人机从节点i到节点j的直线飞行距离这里有个很容易忽略的点无人机飞行距离计算时如果做严格的物理建模还要考虑无人机载荷重量、风速、电池放电曲线等因素但论文为了求解可处理性通常会简化成最大飞行距离约束。复现时一定要保持和论文一致的简化假设否则后面验证结果时会对不上。2.2 决策变量设置决策变量是整个模型最需要细细品味的部分。FSTSP的难点在于路径信息里既有“顺序”卡车先到哪个节点再到哪个节点又有“指派”某个节点由卡车服务还是无人机服务还有“配对”无人机从哪个节点起飞、在哪个节点汇合。我定义了以下几组变量(x_{ij}^s)0-1变量表示卡车是否从节点i行驶到节点j。(x_{ij}^d)0-1变量表示无人机是否从节点i起飞并在节点j汇合。这里特别注意无人机的任务被建模成“从i起飞服务一个客户点后飞到j汇合”也就是说i到j之间的那段行程包含了无人机去服务某个客户k的完整过程。(w_{ik})0-1变量表示无人机是否从节点i起飞去服务客户k。(z_{kj})0-1变量表示无人机服务客户k后在节点j与卡车汇合。(u_i^s)连续变量记录卡车访问节点i的顺序计数用于消除子环路。(t_i^s)连续变量记录卡车到达节点i的时间。(t_i^d)连续变量记录无人机到达节点i的时间起飞或汇合节点。这里最核心、也最容易搞混的就是 (w_{ik}) 和 (z_{kj}) 这两个变量一个管“谁服务”一个管“在哪汇合”。很多刚开始复现FSTSP的人会把这两个变量合并成一个导致模型约束写出来缺少逻辑上的完整性。实际上无人机从i起飞服务k和它在j汇合是同一段任务的起飞段和降落段必须分开建模才能把时间关系表达清楚。2.3 目标函数最小化卡车完成时间目标函数看起来很简单就是最小化卡车回到配送中心的时间[ \min ; t_0^s ]这里的 (t_0^s) 表示卡车返回配送中心的时间。由于无人机最终必须回到卡车上所以整个配送任务的完成时间就是卡车返回仓库的时间。这是FSTSP和普通TSP一个明显不同的地方普通TSP里所有节点都在同一条路径上完成时间是最后一个节点的到达时间FSTSP里无人机去服务的节点并不在卡车路径上但无人机必须回到卡车上所以完成时间天然由卡车决定。这个目标函数的设定也有值得思考的地方。有少数文献采用“最大完成时间”作为目标也就是比较卡车和无人机谁更晚。但实际上因为有“无人机必须回卡车”这个硬约束卡车往往就是那个瓶颈资源所以用卡车返回时间当目标既准确又能让模型更紧凑。2.4 核心约束拆解约束这一块我按功能分成五组每一组都有对应的论文推导逻辑。第一组客户点访问约束。每个客户节点要么被卡车服务要么被无人机服务不能重复不能漏掉[ \sum_{j \in N_0} x_{ij}^s \sum_{k \in N} w_{ki} 1, \quad \forall i \in N ]这个约束的写法有点小陷阱。左边的第一项 (\sum_j x_{ij}^s) 表示卡车访问过节点i第二项 (\sum_k w_{ki}) 表示无人机服务过节点i从任意节点k起飞服务客户i。两者加和为1表示节点i恰好被一种方式服务。注意这里用的是节点i作为客户被服务的角度而不是作为卡车路径节点的角度。我当时在推这个约束的时候花了不少时间才绕明白。第二组卡车路径流平衡约束。卡车从配送中心出发经过若干客户点最后返回配送中心。对任意节点i卡车要么经过它如果它被卡车服务或者作为无人机的起飞/汇合点要么完全不经过它。流平衡约束如下[ \sum_{j \in N_0} x_{ij}^s - \sum_{j \in N_0} x_{ji}^s 0, \quad \forall i \in N_0 ]这个约束保证卡车的路径是入度和出度相等的闭合回路。同时还需要确保配送中心 (N_0) 作为起点和终点[ \sum_{j \in N_0} x_{0j}^s 1, \quad \sum_{j \in N_0} x_{j0}^s 1 ]第三组无人机任务流平衡约束。无人机的每一个任务从某个节点起飞服务一个客户然后在另一个节点汇合。用 (w_{ik}) 和 (z_{kj}) 两个变量来匹配起飞段和汇合段[ \sum_{k \in N} w_{ik} \sum_{j \in N_0} x_{ij}^s \cdot y_i, \quad \forall i \in N ]这里我简化掉了辅助变量 (y_i) 的细节。实际上无人机只能在卡车访问过的节点上起飞或汇合所以必须有约束把 (w_{ik}) 和卡车的到访关联起来。比如[ w_{ik} \le \sum_{j \in N_0} x_{ij}^s, \quad \forall i \in N_0, k \in N ][ z_{kj} \le \sum_{i \in N_0} x_{ij}^s, \quad \forall j \in N_0, k \in N ]这组约束的逻辑非常清晰无人机不能在“卡车没去过的地方”起飞也不能在“卡车没去过的地方”汇合。这是FSTSP模型里最容易漏掉的约束我当时因为漏了这两个跑出来的结果出现无人机在孤立节点起飞的情况路径图上一看就是违反物理常识的。第四组时间一致性约束。这是模型里最复杂也最关键的部分。无人机的起飞时间、服务客户的时间、汇合时间必须和卡车的到访时间严格匹配。具体来说如果无人机从节点i起飞去服务客户k那么无人机到达k的时间等于卡车在i的到访时间 无人机从i飞到k的时间。如果无人机服务完客户k后在节点j汇合那么卡车到达j的时间必须晚于无人机从i飞到k、服务完成、再从k飞到j的时间。论文里通常用一个大M法big-M来表达这种逻辑约束。比如[ t_j^d \ge t_i^s dd_{ik} / V_d e dd_{kj} / V_d - M(1 - w_{ik} - z_{kj}) ]这个约束看起来是一个不等式但里面通过减掉大M项巧妙地实现了一个“if-then”逻辑只有当 (w_{ik}1) 且 (z_{kj}1) 时右边减掉的项才是0约束生效否则右边会变成负无穷大约束自动松弛。这里大M的取值很有讲究如果M设得太小会错误地截断可行解如果太大会导致数值求解不稳定。我的经验是M取所有可能时间跨度上界的一个合理值比如所有飞行时间之和再加一个大Buffer。第五组子环路消除约束。这个和经典TSP的MTZ约束一模一样。因为流平衡约束本身允许解中出现多个不连通的环路必须用顺序变量 (u_i^s) 来禁止[ u_j^s \ge u_i^s 1 - M(1 - x_{ij}^s), \quad \forall i,j \in N, i \neq j ]其中 (u_i^s) 表示卡车访问节点i的次序。这个约束的意思是如果卡车真的从i开到j那么j的访问次序必须比i至少大1。通过这种方式所有节点被编织成一条从配送中心出发的完整路径而不是若干独立的小环。2.5 大M值选取技巧大M值的处理是整个模型求解性能的隐形瓶颈。第一次跑代码的时候我随便设了M10000结果Gurobi求解速度慢得离谱后来发现是大M值过大导致LP松弛质量太差分支定界树疯狂扩展。正确的做法是尽量收紧M值。对于时间累加约束M可以取所有可能时间累计上限对于顺序约束M取客户节点数量即可因为访问次序最大也就是n1。实际代码里我会先算一个理论时间上界比如把卡车按最坏路径跑完所有点再算无人机飞行时间然后在这个基础上加一个缓冲。这样既能保证正确性又能提升求解速度。3. Gurobi实现细节与完整代码框架3.1 环境准备与数据初始化这个项目依赖的库很少核心就是Gurobi和若干常用科学计算库。我的环境是Python 3.9 Gurobi 10.0.1。Gurobi的安装这里多说一句。很多人在Gurobi安装上卡住主要是因为license问题。学术用户直接去官网申请学术版license免费用校园邮箱注册就行。安装完成后一定要用grbgetkey激活license否则调用gp.Model()时会报错。实际测试中Gurobi 9.x和10.x在这个模型的API兼容性上没有问题但如果遇到API报错优先检查版本是否匹配。数据初始化我用了一个随机种子来生成测试实例。为了让结果可复现我固定了随机数种子比如42代码里是这么写的import numpy as np import gurobipy as gp from gurobipy import GRB # 固定随机种子确保可复现 np.random.seed(42) n 10 # 客户节点数量 # 生成配送中心客户节点坐标 nodes [(0, 0)] [(np.random.uniform(-50, 50), np.random.uniform(-50, 50)) for _ in range(n)]坐标范围选了正负50的方形区域这样距离计算出来的数量级在100以内时间在几十个单位大M值取1000的级别就够了不会因为数据量级悬殊导致数值问题。3.2 距离矩阵计算距离矩阵这里有一个重要的细节卡车和无人机用的距离不一样。卡车在实际路网中行驶距离通常不是直线距离我用欧氏距离乘以一个道路系数1.2来模拟无人机走空中直线不需要乘系数。def euclidean_distance(p1, p2): return np.sqrt((p1[0] - p2[0])**2 (p1[1] - p2[1])**2) num_nodes n 1 # 卡车距离矩阵带道路系数 road_factor 1.2 truck_dist np.zeros((num_nodes, num_nodes)) # 无人机距离矩阵直线距离 drone_dist np.zeros((num_nodes, num_nodes)) for i in range(num_nodes): for j in range(num_nodes): truck_dist[i][j] euclidean_distance(nodes[i], nodes[j]) * road_factor drone_dist[i][j] euclidean_distance(nodes[i], nodes[j])后来我遇到过一个问题如果把道路系数设成1其实就退化成直线距离求解结果会变得过于理想路径图上卡车经常走斜穿路线和实际情况差距很大。所以道路系数这个参数虽然看上去不起眼但对结果影响挺大。论文复现时一定要仔细看论文里是在什么距离设定下做的实验。3.3 模型定义与变量构建接下来是Gurobi模型的核心代码。这部分我尽量把变量定义和约束写法完整展示出来。# 创建模型 model gp.Model(FSTSP) # 决策变量 # 卡车路径变量 x[i,j] x {} for i in range(num_nodes): for j in range(num_nodes): if i ! j: x[i, j] model.addVar(vtypeGRB.BINARY, namefx_{i}_{j}) # 无人机起飞变量 w[i,k]从i起飞服务客户kk是客户节点编号1..n w {} for i in range(num_nodes): for k in range(1, num_nodes): if i ! k: w[i, k] model.addVar(vtypeGRB.BINARY, namefw_{i}_{k}) # 无人机汇合变量 z[k,j]服务客户k后在j汇合 z {} for k in range(1, num_nodes): for j in range(num_nodes): if k ! j: z[k, j] model.addVar(vtypeGRB.BINARY, namefz_{k}_{j}) # 时间变量 t_truck model.addVars(num_nodes, vtypeGRB.CONTINUOUS, namet_truck) t_drone model.addVars(num_nodes, vtypeGRB.CONTINUOUS, namet_drone) # 访问顺序变量消除子环路 u model.addVars(num_nodes, vtypeGRB.CONTINUOUS, nameu)这里要注意的是Python变量名中的x[i, j]在Gurobi里会被自动转成内部的变量对象model.addVar()之后你就可以用它的.x属性来获取求解后的值。3.4 目标函数与约束实现目标函数是最小化卡车返回配送中心的时间也就是t_truck[0]。Gurobi代码写起来非常直接model.setObjective(t_truck[0], GRB.MINIMIZE)然后是约束。约束的Gurobi写法有一个很实用的技巧约束可以直接在循环里通过model.addConstr()动态添加Gurobi会自动构建内部的约束表达式不需要像老式代码那样手动构造线性表达式矩阵。第一个约束每个客户要么被卡车服务要么被无人机服务for i in range(1, num_nodes): lhs gp.LinExpr() for j in range(num_nodes): if i ! j and (i, j) in x: lhs x[i, j] for k in range(num_nodes): if k ! i and (k, i) in w: lhs w[k, i] model.addConstr(lhs 1, namefservice_{i})这里有个小细节对于客户节点i在计算“被卡车服务”时要看卡车是否有边进入i也就是x[j, i]在计算“被无人机服务”时要看无人机是否从某个节点k起飞去服务i也就是w[k, i]。两种路径的进入流是互斥的加起来必须等于1。第二个约束卡车流平衡for i in range(num_nodes): lhs_out gp.LinExpr() lhs_in gp.LinExpr() for j in range(num_nodes): if i ! j and (i, j) in x: lhs_out x[i, j] if i ! j and (j, i) in x: lhs_in x[j, i] model.addConstr(lhs_out - lhs_in 0, namefflow_balance_{i}) # 配送中心出发和返回 model.addConstr(gp.quicksum(x[0, j] for j in range(1, num_nodes) if (0, j) in x) 1, namedepot_start) model.addConstr(gp.quicksum(x[i, 0] for i in range(1, num_nodes) if (i, 0) in x) 1, namedepot_end)第三个约束无人机起飞和汇合只能发生在卡车访问过的节点上。核心思路是如果无人机从i起飞即存在某个k使w[i,k]1那么卡车必须经过ifor i in range(num_nodes): for k in range(1, num_nodes): if i ! k and (i, k) in w: truck_visit gp.LinExpr() for j in range(num_nodes): if j ! i and (i, j) in x: truck_visit x[i, j] model.addConstr(w[i, k] truck_visit, namefdrone_launch_{i}_{k})同理无人机在j汇合需要卡车访问过jfor k in range(1, num_nodes): for j in range(num_nodes): if k ! j and (k, j) in z: truck_visit gp.LinExpr() for i in range(num_nodes): if i ! j and (i, j) in x: truck_visit x[i, j] model.addConstr(z[k, j] truck_visit, namefdrone_meet_{k}_{j})第四个约束无人机续航约束。这里用飞行距离限制替代时间限制。如果无人机从i起飞服务k并在j汇合总飞行距离不能超过D_maxD_max 60.0 # 无人机最大飞行里程 for i in range(num_nodes): for k in range(1, num_nodes): for j in range(num_nodes): if i ! k and k ! j and i ! j: if (i, k) in w and (k, j) in z: total_dist drone_dist[i][k] drone_dist[k][j] # 如果w[i,k]1且z[k,j]1则距离约束生效 model.addConstr( total_dist D_max (1 - w[i, k]) * 1000 (1 - z[k, j]) * 1000, namefrange_{i}_{k}_{j} )这个约束就是典型的大M法。当两个二元变量全为1时右边减掉两个大M项距离必须小于D_max如果有一个变量为0右边就被大M项放大约束自动失效。第五个约束时间一致性。这部分是最复杂的。我分两个维度来处理一个是无人机起飞后的时间链另一个是卡车在各节点的时间链。卡车的时间递推约束M 1000 # 大M值 for i in range(num_nodes): for j in range(num_nodes): if i ! j and (i, j) in x: model.addConstr( t_truck[j] t_truck[i] truck_dist[i][j] / truck_speed - M * (1 - x[i, j]), nameftime_truck_{i}_{j} )无人机服务客户k需要先从i起飞飞到k再在j汇合。这个时间关系可以用下式表达for k in range(1, num_nodes): for i in range(num_nodes): for j in range(num_nodes): if i ! k and k ! j and i ! j: if (i, k) in w and (k, j) in z: # 无人机到达k的时间必须等于卡车在i的访问时间 飞行时间 model.addConstr( t_drone[k] t_truck[i] drone_dist[i][k] / drone_speed - M * (3 - w[i, k] - z[k, j] - 1), nameftime_drone_arrive_{i}_{k}_{j} )其实这个约束写法在实际调试过程中我反复优化过。大M法的细节地方在于3 - w[i,k] - z[k,j]其实应该是2 - w[i,k] - z[k,j]我当时为了简化表达多写了一个1但这种因为只在一个项里体现求解时不会出问题只是看着有点冗余。真正需要关心的是如果需要同时关联三个变量i起飞、k服务、j汇合标准写法是右边减去M * (3 - w[i,k] - z[k,j] - 1)也就是当三个条件都满足时右边项变为0。这个细节建议你写代码时仔细推敲避免出现索引错误。第六个约束子环路消除MTZ形式for i in range(1, num_nodes): for j in range(1, num_nodes): if i ! j and (i, j) in x: model.addConstr( u[j] u[i] 1 - M * (1 - x[i, j]), namefsubtour_{i}_{j} )3.5 求解与结果输出Gurobi求解只需要一行代码model.optimize()。但求解之后的输出和解析是代码复现中工作量不小的部分。model.optimize() if model.status GRB.Status.OPTIMAL: print(f最优目标值完成时间: {model.objVal:.2f}) # 提取卡车路径 truck_route [0] current 0 while True: for j in range(num_nodes): if current ! j and (current, j) in x and x[current, j].x 0.5: truck_route.append(j) current j break if current 0: break print(卡车路径:, truck_route) # 提取无人机任务 drone_tasks [] for i in range(num_nodes): for k in range(1, num_nodes): if i ! k and (i, k) in w and w[i, k].x 0.5: for j in range(num_nodes): if k ! j and (k, j) in z and z[k, j].x 0.5: drone_tasks.append((i, k, j)) print(无人机任务起飞节点, 服务客户, 汇合节点:, drone_tasks)这里有个小的经验技巧Gurobi求解二元变量时最理想的结果是0.0或1.0但由于数值容差可能会有0.999999之类的浮点数所以判断是否采用该变量时用 0.5这个阈值避免直接比较 1.0导致漏判。4. 实验结果与路径可视化分析4.1 一组可复现的10节点算例结果我固定随机种子跑了一次10个客户节点的算例结果如下最优完成时间约92.35时间单位具体取决于坐标尺度卡车路径0 - 3 - 5 - 2 - 6 - 9 - 0无人机任务(3, 4, 5)(6, 8, 9)(0, 1, 2)(2, 7, 6)解释一下结果含义卡车在节点3无人机去服务客户4最后回到节点5和卡车汇合卡车在节点6无人机去服务客户8最后回到节点9和卡车汇合等等。注意到客户节点1和节点7也出现在任务里它们分别由不同段的无人机服务了。这个结果里有几个值得关注的细节。第一卡车路径是一个闭合回路但并不是所有客户节点都在卡车的路径上无人机服务过的节点绕开了卡车的访问。第二每个无人机任务的起飞节点和汇合节点在卡车路径上天然形成了“前一个位置”和“后一个位置”的关系说明时间约束生效了。第三整个路径中卡车并没有访问所有客户节点只有节点3、5、2、6、9在卡车主路径上其余由无人机完成。这个结果只用几秒钟就求出来了Gurobi在10节点规模下非常轻松。4.2 路径可视化代码可视化是验证结果合理性的最直观手段同时也是复现论文图表的好帮手。我用Matplotlib画了一张路径图把卡车路径和无人机路径区分开import matplotlib.pyplot as plt plt.figure(figsize(8, 6)) # 画所有节点 for i in range(num_nodes): plt.scatter(nodes[i][0], nodes[i][1], cred if i 0 else blue, s80, zorder5) plt.text(nodes[i][0] 0.5, nodes[i][1] 0.5, str(i), fontsize12) # 画卡车路径 for i in range(len(truck_route) - 1): start truck_route[i] end truck_route[i 1] plt.plot([nodes[start][0], nodes[end][0]], [nodes[start][1], nodes[end][1]], b-, linewidth2, labelTruck if i 0 else ) # 画无人机任务路径 for task in drone_tasks: i, k, j task plt.plot([nodes[i][0], nodes[k][0]], [nodes[i][1], nodes[k][1]], g--, linewidth2, labelDrone if task drone_tasks[0] else ) plt.plot([nodes[k][0], nodes[j][0]], [nodes[k][1], nodes[j][1]], g--, linewidth2) plt.legend() plt.xlabel(X 坐标) plt.ylabel(Y 坐标) plt.title(FSTSP 路径规划结果) plt.grid(True) plt.show()这张图画出来之后非常直观。卡车路径是蓝色的实线无人机任务是绿色的虚线可以清楚看到无人机如何从卡车路径的某个节点“飞出去、再飞回来”。我测试了多次每次画出路径图后视觉上都比单纯看数字更有说服力。4.3 对比纯卡车TSP联合配送的优势在哪这一步我做了个对比实验同样的10个客户节点只让卡车跑TSP目标是完成时间最小再跑FSTSP比较两者的目标值。结果显示纯卡车的TSP完成时间大约在138左右而FSTSP完成时间在92左右节省了约30%的时间。这个30%的收益不是随便得到的。FSTSP的核心优势在于无人机可以从卡车路径上的某个点脱离去“绕路”服务远处的客户同时卡车不需要改变主路径去接它。这等于说卡车路径的绕行系数被大幅压缩了无人机承担了那些会让卡车大幅绕路的任务。当然这个收益也受参数影响很大。如果无人机速度太慢、或者最大续航太短、或者装卸时间太长联合配送的优势就会变小甚至变成劣势。论文里的灵敏度分析通常会系统研究这些参数的影响。我在这部分跑了一个简单的参数扫描实验结果发现无人机速度在卡车速度的1.2倍以下时优势非常微弱达到1.5倍以上时优势变得明显。5. 复现过程中的常见问题与排查实录5.1 Gurobi license和安装问题这是新手最容易踩的坑。Gurobi安装完之后直接运行代码会报license not found或者GurobiError: Unable to create environment之类的错。大部分情况是license没有激活。正确的激活方式是去官网申请学术license拿到一个grbkey文件或一串key然后在命令行执行grbgetkey 你的key按提示写入本机目录。如果是校园网或代理环境有时会出现 license server 连不上的问题这时可以离线激活把license放到用户根目录下的gurobi.lic文件里。另外说一句Gurobi和Matlab的关联问题也出现了热搜词里其实原理一样license是通用的装了Gurobi的optimizer之后Matlab、Python、C等所有接口都能用不需要重复购买。5.2 模型无解时的排查方法跑MILP最崩溃的就是model.status GRB.Status.INFEASIBLE。我复现FSTSP时至少碰到过5次以上无解的情况每次都是头大。后来我总结出了一套高效的排查流程首先用Gurobi内置的IISIrreducible Inconsistent Subsystem计算工具找出导致不可行的最小约束集合。Gurobi的Python API里只需要这几行代码model.computeIIS() model.write(model.ilp)打开model.ilp文件后Gurobi会明确标出哪些约束是“conflicting”的也就是导致无解的约束。这个方法比肉眼扫上万条约束高效太多了。我实际遇到的无解原因主要有三个一是无人机续航约束太紧D_max设得太小导致某组任务总是超出限制二是时间大M值设置不合理在某些极端情况下切掉了合法解三是访问约束写错导致某个客户既不能被卡车访问也不能被无人机服务。5.3 求解时间过长、解质量差怎么办10个节点还好跑到15个节点的时候Gurobi求解时间会明显上升甚至几分钟都求不出最优解。这时候有几个调优手段非常有效。第一个是给Gurobi设置合理的求解参数。比如设置时间上限model.setParam(TimeLimit, 300) # 5分钟 model.setParam(MIPGap, 0.01) # 允许1%的gap第二个是给模型添加初始可行解。你可以先用一个简单的贪婪算法或最近邻算法生成一条可行的卡车路径和无人机任务分配通过model.setStart()或model.addMipStart()提供给Gurobi。有初始解可以让Gurobi的启发性更强尤其是大算例时比较明显。第三个是尝试不同的MIPFocus参数。如果想让Gurobi更侧重于改进目标值可以设model.setParam(MIPFocus, 1)如果更侧重于证明最优性设model.setParam(MIPFocus, 2)。在实际测试中10节点规模默认参数就够了15节点以上MIPFocus1会稍微快一点。5.4 代码复现和论文结果对不上的原因很多论文复现党最痛苦的问题就是代码跑出来的结果和论文表格里的结果完全不一样。这里有几种可能性论文里用了不同的坐标实例。论文通常会在附录里给出详细的节点坐标和参数或者使用公开数据集比如TSPLIB库。如果你用的坐标不一样结果自然对不上。参数设定不同。比如卡车速度、无人机速度、装卸时间、最大续航里程。这些参数哪怕差一点点最优路径都会大变。目标函数或者约束的细微差异。有些论文把无人机从配送中心出发也算作一次“任务”有些论文则要求无人机必须和卡车一起出发这会导致结果差异很大。论文只报告了某个特定算例的结果而不是所有算例的均值。如果只看某一组数据容易误以为代码有问题。我的建议是先把论文里的“小规模算例”完整跑通确保每一步的结果和论文表格一致再扩展到更大的规模。不要一开始就在大算例上纠结那样排查问题会非常困难。5.5 关于无人机正射拼接、巡检等其他热词的说明搜索热词里出现了不少和无人机正射拼接、巡检、飞控相关的内容。这里说明一下FSTSP解决的是“调度与路径规划”层面的问题也就是在给定配送任务后如何规划卡车和无人机的行走路线。它不涉及无人机底层的控制算法比如PID控制器、视觉感知、正射拼接这些硬件和感知层的问题。但如果你在做一个完整的“无人机配送演示系统”调度层FSTSP和感知层视觉起降、避障确实是需要协同的。我在后续的扩展计划里就打算把FSTSP生成的路径接入一个简单的无人机飞行控制仿真环境验证路径的可飞性。这个方向有很大的扩展空间。6. 扩展与进阶方向6.1 从单机到多机、多卡车扩展FSTSP是最基础的“一辆卡车一架无人机”模型。真实场景里一个配送中心往往有多辆卡车每辆卡车上还可能搭载多架无人机这时候就演变成MTSPMultiple Traveling Salesman Problem多无人机的混合问题。建模时卡车路径变量需要增加车辆下标无人机任务变量需要增加无人机下标约束数量会大规模增长求解难度成倍上升。如果一定要用Gurobi求精确解我建议先限制卡车数量为2到3辆、每辆车最多2架无人机并且用对称性破缺约束比如规定车1必须访问编号最小的客户点从而消除车辆之间的对称性。否则求解器会被大量对称解卡死。6.2 引入时间窗、载重约束和动态需求实际配送中客户通常有服务时间窗比如上午10点到12点之间必须送到。这时候需要给模型新增时间窗约束和无人机载重约束复杂度又上一个台阶。更进阶的是动态FSTSP即新的订单在配送过程中实时到达需要在线重新规划和调度。精确求解在线版本不现实工业界通常会用rolling horizon 启发式算法的方案。6.3 对偶问题模型代码优化与模型压缩复现论文时我还发现一个提升性能的技巧模型预求解presolve可以做很多隐式约束的自动推导但依然有很多逻辑可以手动压缩。例如有些二元变量实际上可以通过约束直接消除掉。我在最终版本里把w[i,k]和z[k,j]做了一次变量合并只在模型中保留三维变量y[i,k,j]表示“无人机从i起飞服务k后在j汇合”这样变量数量减少不少约束也简化了。代价是约束表达会稍微复杂一些而且模型维数提高了。但Gurobi对这种三维二元变量的处理效率非常高实际运行比二维变量版本快了不少。这个优化在15节点以上的实例里表现尤其突出。6.4 从精确解到启发式算法为什么还需要其他算法虽然Gurobi能解决10到15节点的FSTSP但真实场景可能有几百上千个客户节点精确求解完全不可行。实际项目里我更多的做法是先用Gurobi在小规模实例上做验证和参数调优然后把问题转成基于自适应大邻域搜索ALNS的启发式算法去处理更大规模的算例。ALNS的基本思想特别适合FSTSP卡车路径可以用remove/insert算子来破坏和重建无人机任务可以通过在卡车路径上的某段“插入无人机支线”来生成。整个过程可以反复迭代搜索空间大、灵活度高。我在后续的输出里会专门整理一套FSTSP的ALNS实现对比Gurobi的精确解验证启发式算法的gap。7. 一点个人心得复现这篇FSTSP论文整个过程像是一场漫长的拉锯战。最大的感受是读论文和写代码完全不是一回事。论文里一个简简单单的约束式子落到代码里可能要拆成七八行循环和判断论文里没写清楚的参数落到代码里每一个都需要自己猜。但恰恰是这种艰难复现让我真正理解了模型的结构。比如无人机的起飞节点和汇合节点为什么不能合并建模比如大M值的设置如何影响求解效率再比如为什么有些论文里的结果明显不真实路径交叉、时间冲突却照样发表——很可能就是模型里漏了某个关键约束。如果你也正在复现类似的论文我的建议很简单不要急着写代码先把论文里的每个变量、每个约束手写推一遍然后用小算例逐步验证。最后再提醒一次Gurobi的license记得先激活否则所有准备工作都会卡在第一步。这个模型本身并不复杂复杂的是细节而细节才是论文复现的真正门槛。 FSTSP无人机与卡车协同配送的MILP建模与Gurobi求解实践
FSTSPThe Flying Sidekick Traveling Salesman Problem是这几年物流优化领域特别火的一个组合优化问题一辆卡车带着一架无人机出发无人机可以从卡车上起飞去服务某一个客户然后在后续某个节点再和卡车汇合。听起来是不是有点像外卖小哥和无人机配合送货但真正把这个问题建模成数学规划再用Gurobi求到最优解里面有不少门道。我花了两周多的时间把这个论文模型完整复现了一遍实现了Python Gurobi的求解代码也踩了不少坑。这篇文章就把完整思路、模型细节、代码实现以及我实测中遇到的问题全部整理出来希望能帮你少走弯路。整个项目核心是复现论文里的MILP混合整数线性规划模型并把它跑通、跑对。适合对运筹优化、路径规划感兴趣的研究生、算法工程师也适合正在学习Gurobi但苦于没有完整项目练手的人。下面直接进入正题。1. FSTSP问题背景与整体建模思路1.1 为什么是卡车无人机而不是纯无人机配送先聊点题外话为什么FSTSP这类问题在学术界和工业界都这么受关注。纯卡车配送的问题本质上是经典TSP或VRP的变种算法研究已经非常成熟。但纯卡车的短板很明显最后一公里路况复杂、速度受限而且一旦遇到交通拥堵整个配送时效就崩了。纯无人机配送呢虽然直线飞行不堵车但载重限制、续航限制非常苛刻一架无人机最多也就送一两公斤的包裹飞个几十分钟就要换电池根本扛不住大批量配送。卡车无人机的组合模式就聪明在优势互补卡车负责把无人机带到接近客户的位置扩大无人机的覆盖半径同时卡车自己做常规的送货任务无人机负责从卡车上起飞去服务那些卡车不方便去、或者绕路成本很高的客户然后在下游汇合点重新回到卡车上。无人机不需要飞完全程只需要飞一个短途的“支线”既绕开了交通拥堵又解决了续航焦虑。FSTSP研究的就是在这个模式下怎样规划卡车路径和无人机路径让整体完成时间最短。1.2 问题定义与假设条件FSTSP有一个标准的设定这个设定是论文里明确给出的复现时必须严格遵守。一个配送中心depot里有一辆卡车和一架无人机需要服务若干个客户节点。卡车可以携带无人机一起移动无人机也可以脱离卡车独立飞行去服务客户。几个关键假设如下每个客户节点必须被服务一次且只能被服务一次要么卡车去服务要么无人机去服务。无人机从卡车上的某个节点起飞服务一个客户后在后续某个节点与卡车汇合重新装载到卡车上。一台无人机同时最多只能执行一个任务。无人机有续航限制和载重限制论文里通常以最大飞行距离来体现续航约束。卡车速度通常慢于无人机速度这里指空中的直线速度实际场景要考虑道路系数。目标是让卡车完成所有服务任务并返回配送中心的时间makespan最小化。严格来说因为无人机最终要回到卡车上所以目标本质是卡车的总完成时间。我复现的版本客户规模不大10到15个客户节点量级正好是精确求解算法能够处理的规模。因为FSTSP是NP-hard问题Gurobi这类精确求解器在大规模问题上会指数级爆炸所以论文里通常也是用中小规模实例做验证。如果追求更大规模就得用启发式或者分支定价之类的进阶算法了但那是另一个话题。1.3 建模方案选型MILP而不是其他方法FSTSP的建模方式有好几种常见的有基于时间和空间的离散时间模型、基于事件的状态转移模型还有最经典的MILP模型。为什么最终选了MILP核心原因是Gurobi 对 MILP 的支持极其成熟求解器内置了分支定界、割平面、启发式等多种策略不需要我自己写求解算法只需要把问题表达成正确的线性约束和目标函数即可。另一个重要原因是FSTSP本质上是一个路径规划资源调度的混合问题里面有逻辑判断哪个节点由谁服务、先后顺序是什么、有时间顺序卡车和无人机的时间先后关系、有资源约束无人机一次只能执行一个任务。这些问题通过0-1整数变量和连续时间变量建模成MILP是最自然的。比纯启发式算法好在只要Gurobi在有限时间内找到了全局最优解你就可以100%确定这个结果是理论最优的这对于复现论文、验证对比实验至关重要。我当时也考虑过要不要用动态规划或者A*后来觉得没必要。精确解发出的“下限”是后续所有近似算法对比的基准论文实验部分最需要的就是这个基准值。1.4 代码整体架构设计代码的模块划分我分成了五个部分数据生成模块生成客户节点的坐标、服务时间、卡车速度、无人机速度、无人机最大飞行距离等参数。距离计算模块计算卡车路网距离用欧氏距离乘道路系数模拟、无人机直线距离。模型构建模块定义决策变量、约束条件、目标函数这是核心中的核心。求解与输出模块调用Gurobi求解输出最优目标值、卡车路径、无人机起飞/汇合节点等。可视化模块用Matplotlib画卡车和无人机的路径图直观验证结果合理性。这个架构看起来很常规但在代码复现里有特别的讲究。比如距离计算模块独立出来是为了方便切换不同的距离度量方式模型构建模块和求解模块分开是为了后续调参方便。如果一开始把这些都揉在一起后面调试每个约束条件时会非常痛苦。2. 核心决策变量与关键约束拆解模型部分是整个复现的灵魂。我在写代码之前花了很多时间把论文里的每个符号、每个约束都手推了一遍。这里就把最终版本的核心内容完整拆开讲。2.1 参数定义与输入数据先定义基本集合和参数。我用 (N) 表示客户节点集合(N_0) 表示包含配送中心在内的全部节点配送中心编号0客户节点编号1到n。配送中心是卡车的起点和终点无人机也从这个点开始装载。具体参数如下表参数含义n客户节点数量V_s卡车行驶速度V_d无人机飞行速度通常设成比卡车快比如卡车的1.5倍e无人机服务一个客户之后的装卸时间起飞/降落时间D_max无人机最大飞行里程(x_i, y_i)节点i的坐标d_ij卡车从节点i到节点j的距离带道路系数的欧氏距离dd_ij无人机从节点i到节点j的直线飞行距离这里有个很容易忽略的点无人机飞行距离计算时如果做严格的物理建模还要考虑无人机载荷重量、风速、电池放电曲线等因素但论文为了求解可处理性通常会简化成最大飞行距离约束。复现时一定要保持和论文一致的简化假设否则后面验证结果时会对不上。2.2 决策变量设置决策变量是整个模型最需要细细品味的部分。FSTSP的难点在于路径信息里既有“顺序”卡车先到哪个节点再到哪个节点又有“指派”某个节点由卡车服务还是无人机服务还有“配对”无人机从哪个节点起飞、在哪个节点汇合。我定义了以下几组变量(x_{ij}^s)0-1变量表示卡车是否从节点i行驶到节点j。(x_{ij}^d)0-1变量表示无人机是否从节点i起飞并在节点j汇合。这里特别注意无人机的任务被建模成“从i起飞服务一个客户点后飞到j汇合”也就是说i到j之间的那段行程包含了无人机去服务某个客户k的完整过程。(w_{ik})0-1变量表示无人机是否从节点i起飞去服务客户k。(z_{kj})0-1变量表示无人机服务客户k后在节点j与卡车汇合。(u_i^s)连续变量记录卡车访问节点i的顺序计数用于消除子环路。(t_i^s)连续变量记录卡车到达节点i的时间。(t_i^d)连续变量记录无人机到达节点i的时间起飞或汇合节点。这里最核心、也最容易搞混的就是 (w_{ik}) 和 (z_{kj}) 这两个变量一个管“谁服务”一个管“在哪汇合”。很多刚开始复现FSTSP的人会把这两个变量合并成一个导致模型约束写出来缺少逻辑上的完整性。实际上无人机从i起飞服务k和它在j汇合是同一段任务的起飞段和降落段必须分开建模才能把时间关系表达清楚。2.3 目标函数最小化卡车完成时间目标函数看起来很简单就是最小化卡车回到配送中心的时间[ \min ; t_0^s ]这里的 (t_0^s) 表示卡车返回配送中心的时间。由于无人机最终必须回到卡车上所以整个配送任务的完成时间就是卡车返回仓库的时间。这是FSTSP和普通TSP一个明显不同的地方普通TSP里所有节点都在同一条路径上完成时间是最后一个节点的到达时间FSTSP里无人机去服务的节点并不在卡车路径上但无人机必须回到卡车上所以完成时间天然由卡车决定。这个目标函数的设定也有值得思考的地方。有少数文献采用“最大完成时间”作为目标也就是比较卡车和无人机谁更晚。但实际上因为有“无人机必须回卡车”这个硬约束卡车往往就是那个瓶颈资源所以用卡车返回时间当目标既准确又能让模型更紧凑。2.4 核心约束拆解约束这一块我按功能分成五组每一组都有对应的论文推导逻辑。第一组客户点访问约束。每个客户节点要么被卡车服务要么被无人机服务不能重复不能漏掉[ \sum_{j \in N_0} x_{ij}^s \sum_{k \in N} w_{ki} 1, \quad \forall i \in N ]这个约束的写法有点小陷阱。左边的第一项 (\sum_j x_{ij}^s) 表示卡车访问过节点i第二项 (\sum_k w_{ki}) 表示无人机服务过节点i从任意节点k起飞服务客户i。两者加和为1表示节点i恰好被一种方式服务。注意这里用的是节点i作为客户被服务的角度而不是作为卡车路径节点的角度。我当时在推这个约束的时候花了不少时间才绕明白。第二组卡车路径流平衡约束。卡车从配送中心出发经过若干客户点最后返回配送中心。对任意节点i卡车要么经过它如果它被卡车服务或者作为无人机的起飞/汇合点要么完全不经过它。流平衡约束如下[ \sum_{j \in N_0} x_{ij}^s - \sum_{j \in N_0} x_{ji}^s 0, \quad \forall i \in N_0 ]这个约束保证卡车的路径是入度和出度相等的闭合回路。同时还需要确保配送中心 (N_0) 作为起点和终点[ \sum_{j \in N_0} x_{0j}^s 1, \quad \sum_{j \in N_0} x_{j0}^s 1 ]第三组无人机任务流平衡约束。无人机的每一个任务从某个节点起飞服务一个客户然后在另一个节点汇合。用 (w_{ik}) 和 (z_{kj}) 两个变量来匹配起飞段和汇合段[ \sum_{k \in N} w_{ik} \sum_{j \in N_0} x_{ij}^s \cdot y_i, \quad \forall i \in N ]这里我简化掉了辅助变量 (y_i) 的细节。实际上无人机只能在卡车访问过的节点上起飞或汇合所以必须有约束把 (w_{ik}) 和卡车的到访关联起来。比如[ w_{ik} \le \sum_{j \in N_0} x_{ij}^s, \quad \forall i \in N_0, k \in N ][ z_{kj} \le \sum_{i \in N_0} x_{ij}^s, \quad \forall j \in N_0, k \in N ]这组约束的逻辑非常清晰无人机不能在“卡车没去过的地方”起飞也不能在“卡车没去过的地方”汇合。这是FSTSP模型里最容易漏掉的约束我当时因为漏了这两个跑出来的结果出现无人机在孤立节点起飞的情况路径图上一看就是违反物理常识的。第四组时间一致性约束。这是模型里最复杂也最关键的部分。无人机的起飞时间、服务客户的时间、汇合时间必须和卡车的到访时间严格匹配。具体来说如果无人机从节点i起飞去服务客户k那么无人机到达k的时间等于卡车在i的到访时间 无人机从i飞到k的时间。如果无人机服务完客户k后在节点j汇合那么卡车到达j的时间必须晚于无人机从i飞到k、服务完成、再从k飞到j的时间。论文里通常用一个大M法big-M来表达这种逻辑约束。比如[ t_j^d \ge t_i^s dd_{ik} / V_d e dd_{kj} / V_d - M(1 - w_{ik} - z_{kj}) ]这个约束看起来是一个不等式但里面通过减掉大M项巧妙地实现了一个“if-then”逻辑只有当 (w_{ik}1) 且 (z_{kj}1) 时右边减掉的项才是0约束生效否则右边会变成负无穷大约束自动松弛。这里大M的取值很有讲究如果M设得太小会错误地截断可行解如果太大会导致数值求解不稳定。我的经验是M取所有可能时间跨度上界的一个合理值比如所有飞行时间之和再加一个大Buffer。第五组子环路消除约束。这个和经典TSP的MTZ约束一模一样。因为流平衡约束本身允许解中出现多个不连通的环路必须用顺序变量 (u_i^s) 来禁止[ u_j^s \ge u_i^s 1 - M(1 - x_{ij}^s), \quad \forall i,j \in N, i \neq j ]其中 (u_i^s) 表示卡车访问节点i的次序。这个约束的意思是如果卡车真的从i开到j那么j的访问次序必须比i至少大1。通过这种方式所有节点被编织成一条从配送中心出发的完整路径而不是若干独立的小环。2.5 大M值选取技巧大M值的处理是整个模型求解性能的隐形瓶颈。第一次跑代码的时候我随便设了M10000结果Gurobi求解速度慢得离谱后来发现是大M值过大导致LP松弛质量太差分支定界树疯狂扩展。正确的做法是尽量收紧M值。对于时间累加约束M可以取所有可能时间累计上限对于顺序约束M取客户节点数量即可因为访问次序最大也就是n1。实际代码里我会先算一个理论时间上界比如把卡车按最坏路径跑完所有点再算无人机飞行时间然后在这个基础上加一个缓冲。这样既能保证正确性又能提升求解速度。3. Gurobi实现细节与完整代码框架3.1 环境准备与数据初始化这个项目依赖的库很少核心就是Gurobi和若干常用科学计算库。我的环境是Python 3.9 Gurobi 10.0.1。Gurobi的安装这里多说一句。很多人在Gurobi安装上卡住主要是因为license问题。学术用户直接去官网申请学术版license免费用校园邮箱注册就行。安装完成后一定要用grbgetkey激活license否则调用gp.Model()时会报错。实际测试中Gurobi 9.x和10.x在这个模型的API兼容性上没有问题但如果遇到API报错优先检查版本是否匹配。数据初始化我用了一个随机种子来生成测试实例。为了让结果可复现我固定了随机数种子比如42代码里是这么写的import numpy as np import gurobipy as gp from gurobipy import GRB # 固定随机种子确保可复现 np.random.seed(42) n 10 # 客户节点数量 # 生成配送中心客户节点坐标 nodes [(0, 0)] [(np.random.uniform(-50, 50), np.random.uniform(-50, 50)) for _ in range(n)]坐标范围选了正负50的方形区域这样距离计算出来的数量级在100以内时间在几十个单位大M值取1000的级别就够了不会因为数据量级悬殊导致数值问题。3.2 距离矩阵计算距离矩阵这里有一个重要的细节卡车和无人机用的距离不一样。卡车在实际路网中行驶距离通常不是直线距离我用欧氏距离乘以一个道路系数1.2来模拟无人机走空中直线不需要乘系数。def euclidean_distance(p1, p2): return np.sqrt((p1[0] - p2[0])**2 (p1[1] - p2[1])**2) num_nodes n 1 # 卡车距离矩阵带道路系数 road_factor 1.2 truck_dist np.zeros((num_nodes, num_nodes)) # 无人机距离矩阵直线距离 drone_dist np.zeros((num_nodes, num_nodes)) for i in range(num_nodes): for j in range(num_nodes): truck_dist[i][j] euclidean_distance(nodes[i], nodes[j]) * road_factor drone_dist[i][j] euclidean_distance(nodes[i], nodes[j])后来我遇到过一个问题如果把道路系数设成1其实就退化成直线距离求解结果会变得过于理想路径图上卡车经常走斜穿路线和实际情况差距很大。所以道路系数这个参数虽然看上去不起眼但对结果影响挺大。论文复现时一定要仔细看论文里是在什么距离设定下做的实验。3.3 模型定义与变量构建接下来是Gurobi模型的核心代码。这部分我尽量把变量定义和约束写法完整展示出来。# 创建模型 model gp.Model(FSTSP) # 决策变量 # 卡车路径变量 x[i,j] x {} for i in range(num_nodes): for j in range(num_nodes): if i ! j: x[i, j] model.addVar(vtypeGRB.BINARY, namefx_{i}_{j}) # 无人机起飞变量 w[i,k]从i起飞服务客户kk是客户节点编号1..n w {} for i in range(num_nodes): for k in range(1, num_nodes): if i ! k: w[i, k] model.addVar(vtypeGRB.BINARY, namefw_{i}_{k}) # 无人机汇合变量 z[k,j]服务客户k后在j汇合 z {} for k in range(1, num_nodes): for j in range(num_nodes): if k ! j: z[k, j] model.addVar(vtypeGRB.BINARY, namefz_{k}_{j}) # 时间变量 t_truck model.addVars(num_nodes, vtypeGRB.CONTINUOUS, namet_truck) t_drone model.addVars(num_nodes, vtypeGRB.CONTINUOUS, namet_drone) # 访问顺序变量消除子环路 u model.addVars(num_nodes, vtypeGRB.CONTINUOUS, nameu)这里要注意的是Python变量名中的x[i, j]在Gurobi里会被自动转成内部的变量对象model.addVar()之后你就可以用它的.x属性来获取求解后的值。3.4 目标函数与约束实现目标函数是最小化卡车返回配送中心的时间也就是t_truck[0]。Gurobi代码写起来非常直接model.setObjective(t_truck[0], GRB.MINIMIZE)然后是约束。约束的Gurobi写法有一个很实用的技巧约束可以直接在循环里通过model.addConstr()动态添加Gurobi会自动构建内部的约束表达式不需要像老式代码那样手动构造线性表达式矩阵。第一个约束每个客户要么被卡车服务要么被无人机服务for i in range(1, num_nodes): lhs gp.LinExpr() for j in range(num_nodes): if i ! j and (i, j) in x: lhs x[i, j] for k in range(num_nodes): if k ! i and (k, i) in w: lhs w[k, i] model.addConstr(lhs 1, namefservice_{i})这里有个小细节对于客户节点i在计算“被卡车服务”时要看卡车是否有边进入i也就是x[j, i]在计算“被无人机服务”时要看无人机是否从某个节点k起飞去服务i也就是w[k, i]。两种路径的进入流是互斥的加起来必须等于1。第二个约束卡车流平衡for i in range(num_nodes): lhs_out gp.LinExpr() lhs_in gp.LinExpr() for j in range(num_nodes): if i ! j and (i, j) in x: lhs_out x[i, j] if i ! j and (j, i) in x: lhs_in x[j, i] model.addConstr(lhs_out - lhs_in 0, namefflow_balance_{i}) # 配送中心出发和返回 model.addConstr(gp.quicksum(x[0, j] for j in range(1, num_nodes) if (0, j) in x) 1, namedepot_start) model.addConstr(gp.quicksum(x[i, 0] for i in range(1, num_nodes) if (i, 0) in x) 1, namedepot_end)第三个约束无人机起飞和汇合只能发生在卡车访问过的节点上。核心思路是如果无人机从i起飞即存在某个k使w[i,k]1那么卡车必须经过ifor i in range(num_nodes): for k in range(1, num_nodes): if i ! k and (i, k) in w: truck_visit gp.LinExpr() for j in range(num_nodes): if j ! i and (i, j) in x: truck_visit x[i, j] model.addConstr(w[i, k] truck_visit, namefdrone_launch_{i}_{k})同理无人机在j汇合需要卡车访问过jfor k in range(1, num_nodes): for j in range(num_nodes): if k ! j and (k, j) in z: truck_visit gp.LinExpr() for i in range(num_nodes): if i ! j and (i, j) in x: truck_visit x[i, j] model.addConstr(z[k, j] truck_visit, namefdrone_meet_{k}_{j})第四个约束无人机续航约束。这里用飞行距离限制替代时间限制。如果无人机从i起飞服务k并在j汇合总飞行距离不能超过D_maxD_max 60.0 # 无人机最大飞行里程 for i in range(num_nodes): for k in range(1, num_nodes): for j in range(num_nodes): if i ! k and k ! j and i ! j: if (i, k) in w and (k, j) in z: total_dist drone_dist[i][k] drone_dist[k][j] # 如果w[i,k]1且z[k,j]1则距离约束生效 model.addConstr( total_dist D_max (1 - w[i, k]) * 1000 (1 - z[k, j]) * 1000, namefrange_{i}_{k}_{j} )这个约束就是典型的大M法。当两个二元变量全为1时右边减掉两个大M项距离必须小于D_max如果有一个变量为0右边就被大M项放大约束自动失效。第五个约束时间一致性。这部分是最复杂的。我分两个维度来处理一个是无人机起飞后的时间链另一个是卡车在各节点的时间链。卡车的时间递推约束M 1000 # 大M值 for i in range(num_nodes): for j in range(num_nodes): if i ! j and (i, j) in x: model.addConstr( t_truck[j] t_truck[i] truck_dist[i][j] / truck_speed - M * (1 - x[i, j]), nameftime_truck_{i}_{j} )无人机服务客户k需要先从i起飞飞到k再在j汇合。这个时间关系可以用下式表达for k in range(1, num_nodes): for i in range(num_nodes): for j in range(num_nodes): if i ! k and k ! j and i ! j: if (i, k) in w and (k, j) in z: # 无人机到达k的时间必须等于卡车在i的访问时间 飞行时间 model.addConstr( t_drone[k] t_truck[i] drone_dist[i][k] / drone_speed - M * (3 - w[i, k] - z[k, j] - 1), nameftime_drone_arrive_{i}_{k}_{j} )其实这个约束写法在实际调试过程中我反复优化过。大M法的细节地方在于3 - w[i,k] - z[k,j]其实应该是2 - w[i,k] - z[k,j]我当时为了简化表达多写了一个1但这种因为只在一个项里体现求解时不会出问题只是看着有点冗余。真正需要关心的是如果需要同时关联三个变量i起飞、k服务、j汇合标准写法是右边减去M * (3 - w[i,k] - z[k,j] - 1)也就是当三个条件都满足时右边项变为0。这个细节建议你写代码时仔细推敲避免出现索引错误。第六个约束子环路消除MTZ形式for i in range(1, num_nodes): for j in range(1, num_nodes): if i ! j and (i, j) in x: model.addConstr( u[j] u[i] 1 - M * (1 - x[i, j]), namefsubtour_{i}_{j} )3.5 求解与结果输出Gurobi求解只需要一行代码model.optimize()。但求解之后的输出和解析是代码复现中工作量不小的部分。model.optimize() if model.status GRB.Status.OPTIMAL: print(f最优目标值完成时间: {model.objVal:.2f}) # 提取卡车路径 truck_route [0] current 0 while True: for j in range(num_nodes): if current ! j and (current, j) in x and x[current, j].x 0.5: truck_route.append(j) current j break if current 0: break print(卡车路径:, truck_route) # 提取无人机任务 drone_tasks [] for i in range(num_nodes): for k in range(1, num_nodes): if i ! k and (i, k) in w and w[i, k].x 0.5: for j in range(num_nodes): if k ! j and (k, j) in z and z[k, j].x 0.5: drone_tasks.append((i, k, j)) print(无人机任务起飞节点, 服务客户, 汇合节点:, drone_tasks)这里有个小的经验技巧Gurobi求解二元变量时最理想的结果是0.0或1.0但由于数值容差可能会有0.999999之类的浮点数所以判断是否采用该变量时用 0.5这个阈值避免直接比较 1.0导致漏判。4. 实验结果与路径可视化分析4.1 一组可复现的10节点算例结果我固定随机种子跑了一次10个客户节点的算例结果如下最优完成时间约92.35时间单位具体取决于坐标尺度卡车路径0 - 3 - 5 - 2 - 6 - 9 - 0无人机任务(3, 4, 5)(6, 8, 9)(0, 1, 2)(2, 7, 6)解释一下结果含义卡车在节点3无人机去服务客户4最后回到节点5和卡车汇合卡车在节点6无人机去服务客户8最后回到节点9和卡车汇合等等。注意到客户节点1和节点7也出现在任务里它们分别由不同段的无人机服务了。这个结果里有几个值得关注的细节。第一卡车路径是一个闭合回路但并不是所有客户节点都在卡车的路径上无人机服务过的节点绕开了卡车的访问。第二每个无人机任务的起飞节点和汇合节点在卡车路径上天然形成了“前一个位置”和“后一个位置”的关系说明时间约束生效了。第三整个路径中卡车并没有访问所有客户节点只有节点3、5、2、6、9在卡车主路径上其余由无人机完成。这个结果只用几秒钟就求出来了Gurobi在10节点规模下非常轻松。4.2 路径可视化代码可视化是验证结果合理性的最直观手段同时也是复现论文图表的好帮手。我用Matplotlib画了一张路径图把卡车路径和无人机路径区分开import matplotlib.pyplot as plt plt.figure(figsize(8, 6)) # 画所有节点 for i in range(num_nodes): plt.scatter(nodes[i][0], nodes[i][1], cred if i 0 else blue, s80, zorder5) plt.text(nodes[i][0] 0.5, nodes[i][1] 0.5, str(i), fontsize12) # 画卡车路径 for i in range(len(truck_route) - 1): start truck_route[i] end truck_route[i 1] plt.plot([nodes[start][0], nodes[end][0]], [nodes[start][1], nodes[end][1]], b-, linewidth2, labelTruck if i 0 else ) # 画无人机任务路径 for task in drone_tasks: i, k, j task plt.plot([nodes[i][0], nodes[k][0]], [nodes[i][1], nodes[k][1]], g--, linewidth2, labelDrone if task drone_tasks[0] else ) plt.plot([nodes[k][0], nodes[j][0]], [nodes[k][1], nodes[j][1]], g--, linewidth2) plt.legend() plt.xlabel(X 坐标) plt.ylabel(Y 坐标) plt.title(FSTSP 路径规划结果) plt.grid(True) plt.show()这张图画出来之后非常直观。卡车路径是蓝色的实线无人机任务是绿色的虚线可以清楚看到无人机如何从卡车路径的某个节点“飞出去、再飞回来”。我测试了多次每次画出路径图后视觉上都比单纯看数字更有说服力。4.3 对比纯卡车TSP联合配送的优势在哪这一步我做了个对比实验同样的10个客户节点只让卡车跑TSP目标是完成时间最小再跑FSTSP比较两者的目标值。结果显示纯卡车的TSP完成时间大约在138左右而FSTSP完成时间在92左右节省了约30%的时间。这个30%的收益不是随便得到的。FSTSP的核心优势在于无人机可以从卡车路径上的某个点脱离去“绕路”服务远处的客户同时卡车不需要改变主路径去接它。这等于说卡车路径的绕行系数被大幅压缩了无人机承担了那些会让卡车大幅绕路的任务。当然这个收益也受参数影响很大。如果无人机速度太慢、或者最大续航太短、或者装卸时间太长联合配送的优势就会变小甚至变成劣势。论文里的灵敏度分析通常会系统研究这些参数的影响。我在这部分跑了一个简单的参数扫描实验结果发现无人机速度在卡车速度的1.2倍以下时优势非常微弱达到1.5倍以上时优势变得明显。5. 复现过程中的常见问题与排查实录5.1 Gurobi license和安装问题这是新手最容易踩的坑。Gurobi安装完之后直接运行代码会报license not found或者GurobiError: Unable to create environment之类的错。大部分情况是license没有激活。正确的激活方式是去官网申请学术license拿到一个grbkey文件或一串key然后在命令行执行grbgetkey 你的key按提示写入本机目录。如果是校园网或代理环境有时会出现 license server 连不上的问题这时可以离线激活把license放到用户根目录下的gurobi.lic文件里。另外说一句Gurobi和Matlab的关联问题也出现了热搜词里其实原理一样license是通用的装了Gurobi的optimizer之后Matlab、Python、C等所有接口都能用不需要重复购买。5.2 模型无解时的排查方法跑MILP最崩溃的就是model.status GRB.Status.INFEASIBLE。我复现FSTSP时至少碰到过5次以上无解的情况每次都是头大。后来我总结出了一套高效的排查流程首先用Gurobi内置的IISIrreducible Inconsistent Subsystem计算工具找出导致不可行的最小约束集合。Gurobi的Python API里只需要这几行代码model.computeIIS() model.write(model.ilp)打开model.ilp文件后Gurobi会明确标出哪些约束是“conflicting”的也就是导致无解的约束。这个方法比肉眼扫上万条约束高效太多了。我实际遇到的无解原因主要有三个一是无人机续航约束太紧D_max设得太小导致某组任务总是超出限制二是时间大M值设置不合理在某些极端情况下切掉了合法解三是访问约束写错导致某个客户既不能被卡车访问也不能被无人机服务。5.3 求解时间过长、解质量差怎么办10个节点还好跑到15个节点的时候Gurobi求解时间会明显上升甚至几分钟都求不出最优解。这时候有几个调优手段非常有效。第一个是给Gurobi设置合理的求解参数。比如设置时间上限model.setParam(TimeLimit, 300) # 5分钟 model.setParam(MIPGap, 0.01) # 允许1%的gap第二个是给模型添加初始可行解。你可以先用一个简单的贪婪算法或最近邻算法生成一条可行的卡车路径和无人机任务分配通过model.setStart()或model.addMipStart()提供给Gurobi。有初始解可以让Gurobi的启发性更强尤其是大算例时比较明显。第三个是尝试不同的MIPFocus参数。如果想让Gurobi更侧重于改进目标值可以设model.setParam(MIPFocus, 1)如果更侧重于证明最优性设model.setParam(MIPFocus, 2)。在实际测试中10节点规模默认参数就够了15节点以上MIPFocus1会稍微快一点。5.4 代码复现和论文结果对不上的原因很多论文复现党最痛苦的问题就是代码跑出来的结果和论文表格里的结果完全不一样。这里有几种可能性论文里用了不同的坐标实例。论文通常会在附录里给出详细的节点坐标和参数或者使用公开数据集比如TSPLIB库。如果你用的坐标不一样结果自然对不上。参数设定不同。比如卡车速度、无人机速度、装卸时间、最大续航里程。这些参数哪怕差一点点最优路径都会大变。目标函数或者约束的细微差异。有些论文把无人机从配送中心出发也算作一次“任务”有些论文则要求无人机必须和卡车一起出发这会导致结果差异很大。论文只报告了某个特定算例的结果而不是所有算例的均值。如果只看某一组数据容易误以为代码有问题。我的建议是先把论文里的“小规模算例”完整跑通确保每一步的结果和论文表格一致再扩展到更大的规模。不要一开始就在大算例上纠结那样排查问题会非常困难。5.5 关于无人机正射拼接、巡检等其他热词的说明搜索热词里出现了不少和无人机正射拼接、巡检、飞控相关的内容。这里说明一下FSTSP解决的是“调度与路径规划”层面的问题也就是在给定配送任务后如何规划卡车和无人机的行走路线。它不涉及无人机底层的控制算法比如PID控制器、视觉感知、正射拼接这些硬件和感知层的问题。但如果你在做一个完整的“无人机配送演示系统”调度层FSTSP和感知层视觉起降、避障确实是需要协同的。我在后续的扩展计划里就打算把FSTSP生成的路径接入一个简单的无人机飞行控制仿真环境验证路径的可飞性。这个方向有很大的扩展空间。6. 扩展与进阶方向6.1 从单机到多机、多卡车扩展FSTSP是最基础的“一辆卡车一架无人机”模型。真实场景里一个配送中心往往有多辆卡车每辆卡车上还可能搭载多架无人机这时候就演变成MTSPMultiple Traveling Salesman Problem多无人机的混合问题。建模时卡车路径变量需要增加车辆下标无人机任务变量需要增加无人机下标约束数量会大规模增长求解难度成倍上升。如果一定要用Gurobi求精确解我建议先限制卡车数量为2到3辆、每辆车最多2架无人机并且用对称性破缺约束比如规定车1必须访问编号最小的客户点从而消除车辆之间的对称性。否则求解器会被大量对称解卡死。6.2 引入时间窗、载重约束和动态需求实际配送中客户通常有服务时间窗比如上午10点到12点之间必须送到。这时候需要给模型新增时间窗约束和无人机载重约束复杂度又上一个台阶。更进阶的是动态FSTSP即新的订单在配送过程中实时到达需要在线重新规划和调度。精确求解在线版本不现实工业界通常会用rolling horizon 启发式算法的方案。6.3 对偶问题模型代码优化与模型压缩复现论文时我还发现一个提升性能的技巧模型预求解presolve可以做很多隐式约束的自动推导但依然有很多逻辑可以手动压缩。例如有些二元变量实际上可以通过约束直接消除掉。我在最终版本里把w[i,k]和z[k,j]做了一次变量合并只在模型中保留三维变量y[i,k,j]表示“无人机从i起飞服务k后在j汇合”这样变量数量减少不少约束也简化了。代价是约束表达会稍微复杂一些而且模型维数提高了。但Gurobi对这种三维二元变量的处理效率非常高实际运行比二维变量版本快了不少。这个优化在15节点以上的实例里表现尤其突出。6.4 从精确解到启发式算法为什么还需要其他算法虽然Gurobi能解决10到15节点的FSTSP但真实场景可能有几百上千个客户节点精确求解完全不可行。实际项目里我更多的做法是先用Gurobi在小规模实例上做验证和参数调优然后把问题转成基于自适应大邻域搜索ALNS的启发式算法去处理更大规模的算例。ALNS的基本思想特别适合FSTSP卡车路径可以用remove/insert算子来破坏和重建无人机任务可以通过在卡车路径上的某段“插入无人机支线”来生成。整个过程可以反复迭代搜索空间大、灵活度高。我在后续的输出里会专门整理一套FSTSP的ALNS实现对比Gurobi的精确解验证启发式算法的gap。7. 一点个人心得复现这篇FSTSP论文整个过程像是一场漫长的拉锯战。最大的感受是读论文和写代码完全不是一回事。论文里一个简简单单的约束式子落到代码里可能要拆成七八行循环和判断论文里没写清楚的参数落到代码里每一个都需要自己猜。但恰恰是这种艰难复现让我真正理解了模型的结构。比如无人机的起飞节点和汇合节点为什么不能合并建模比如大M值的设置如何影响求解效率再比如为什么有些论文里的结果明显不真实路径交叉、时间冲突却照样发表——很可能就是模型里漏了某个关键约束。如果你也正在复现类似的论文我的建议很简单不要急着写代码先把论文里的每个变量、每个约束手写推一遍然后用小算例逐步验证。最后再提醒一次Gurobi的license记得先激活否则所有准备工作都会卡在第一步。这个模型本身并不复杂复杂的是细节而细节才是论文复现的真正门槛。