复现论文这件事最怕的就是对着公式看了一整天代码写出来却跑不通。上周我终于把 FSTSPThe Flying Sidekick Traveling Salesman Problem这篇论文的核心模型用 Python Gurobi 完整跑通了整个过程踩了不少坑。今天把完整的复现思路、建模细节和代码实现整理出来给同样在做无人机与卡车联合配送方向的同学一个可以直接参考的作业。先说结论FSTSP 是 Murray 和 Chu 在 2015 年提出的经典模型场景是一辆卡车搭载一架无人机两者协作完成客户的配送任务。无人机可以从卡车上的停靠点起飞独立服务某个客户后再在另一个停靠点回收目标是让整个车队的总完成时间最短。这个问题的难点在于卡车的路径本质上是 TSP和无人机的起飞/回收决策是耦合在一起的属于典型的高复杂度组合优化问题。文章所有代码都基于 Python 3.9 Gurobi 10.x随机生成的小规模算例15个客户可以在几秒内求出精确最优解。如果你正在做物流调度、无人机路径规划或者单纯想学 Gurobi 怎么建混合整数规划模型这篇内容应该能帮你省下不少排查时间。1. FSTSP 问题快速导航联合配送到底在优化什么1.1 论文原始模型核心FSTSP 的全称是 Flying Sidekick Traveling Salesman Problem可以理解为“带飞行搭档的旅行商问题”。传统 TSP 只有一辆车在客户之间移动而 FSTSP 引入了无人机作为卡车的“搭档”整个车队从仓库出发最终回到仓库终点。卡车按顺序访问一部分客户无人机在卡车的停靠节点上发射。无人机一次只能服务一个客户服务完成后必须回到卡车上可以是另一个停靠节点然后继续跟随卡车移动。无人机有续航限制以最大飞行时间表示载重限制通常也隐含在其中。目标是最小化车队完成所有任务的时间也就是卡车的总返回时间。这里面最关键的点是无人机服务的客户不需要卡车亲自去停靠。比如卡车在节点 A无人机起飞去服务一个离主路很远的客户 X然后飞到前方的节点 B 与卡车会合这样卡车就不需要绕路去 X从而缩短总路线。论文模型用混合整数规划MIP来刻画这种协同。决策分为三层卡车的行驶路径即访问客户的先后顺序。无人机发射节点、服务客户、回收节点的三元组合。每个节点的时间先后关系。这三维信息耦合在一起导致问题规模迅速膨胀。也是为什么论文里通常只验证小规模算例因为哪怕只有 20 个客户纯精确解法的求解时间也开始变得可观。1.2 复现目标与技术选型这次复现我给自己定的目标很明确不追求工程化不搞复杂的启发式算法就是用 Gurobi 把论文中的 MIP 模型原汁原味地实现出来并且在随机生成的算例上跑出最优解。技术选型上基本没有犹豫建模语言直接用 Python生态成熟写起来快github 上一堆开源项目也都是这个套路。求解器选 Gurobi原因是学术 License 免费Python API 好用MIP 求解性能在同类商业求解器里属于第一梯队。距离矩阵用 numpy 计算结果可视化用 matplotlib。如果你手头只有 CPLEX 或者 SCIP也可以按同样思路建模只是 API 写法不同核心数学模型是一致的。2. 环境搭建 Gurobi License 避开 80% 的坑2.1 为什么选 Gurobi 而不是免费求解器我见过不少同学用 PuLP 或者 OR-Tools 来尝试复现论文模型结果算到一半发现规模稍大就卡死。原因很简单论文里的 MIP 模型如果约束写得不够紧凑求解器处理起来差别非常大。Gurobi 在预求解presolve、割平面、启发式、并行计算上的积累非常厚同样的模型在 Gurobi 上可能几秒就出最优解换一个开源求解器可能要跑几个小时。另外一个很实际的原因是Gurobi 的 Python API 对数学建模非常友好。你用纸笔写下变量和约束差不多可以直接翻译成代码心智负担小适合快速验证思路。学术 License 是免费的只要能证明你是学生、教师或者科研机构的研究人员就可以申请 full license对论文复现来说完全够用。注意申请学术 License 时需要使用学校邮箱或者机构 IP注册完成后会给你一个 license key用grbgetkey激活有效期为一年到期续期即可。2.2 安装与 License 激活实操安装过程比想象中简单但要按照顺序来不然容易遇到“import 失败”这类问题。第一步安装 Python 环境。建议用 Anaconda 统一管理避免系统 Python 和各种依赖互相污染。装好之后确认 Python 版本我用的 3.9Gurobi 10.x 支持的版本范围非常广3.8 到 3.12 都可以。第二步安装 Gurobi 求解器本体。去官方网站下载对应平台的安装包Windows 直接下一步macOS 和 Linux 用命令行解压即可。安装完成后系统里会有一个gurobi目录里面包含求解器二进制文件、Python 接口、文档等。第三步申请 License 并激活。重点说下激活流程因为这一步最容易出问题# 在命令行中进入你下载的 grbgetkey 所在目录或者直接使用 Gurobi 安装目录下的工具 grbgetkey xxxxxxxx-xxxx-xxxx-xxxx-xxxxxxxxxxxx执行命令后按提示输入 License 文件的保存路径默认是当前用户的~/gurobi.lic。激活完后Gurobi 会自动检测这个文件不需要手动配置环境变量。第四步安装 Python 接口pip install gurobipy装完以后跑一个简单的验证import gurobipy as gp from gurobipy import GRB m gp.Model(test) x m.addVar(vtypeGRB.CONTINUOUS, namex) y m.addVar(vtypeGRB.CONTINUOUS, namey) m.addConstr(x y 1) m.addConstr(x 2) m.setObjective(x y, GRB.MINIMIZE) m.optimize() print(x.X, y.X)如果能正常输出结果说明 Gurobi 安装成功。2.3 环境验证很多新手在import gurobipy时遇到ModuleNotFoundError多半是下面几个原因装了多个 Python 版本pip 和 python 不对应。建议用python -m pip install gurobipy而不是pip install gurobipy。Gurobi 安装目录里的 Python 接口和当前环境的 Python 版本不匹配。如果遇到这种情况卸载重装对应版本的 gurobipy 即可。在 Jupyter Notebook 中运行的 kernel 不是当前环境需要在 Notebook 里检查python -m pip --version。验证通过后就可以进入正题了。3. FSTSP 建模变量、约束、线性化一个都不能少3.1 决策变量设计三元变量的巧思FSTSP 的建模难点在于无人机行动。一个完整的无人机行动可以描述为在节点 i 发射服务客户 k在节点 j 回收。如果只用二元变量分别表示“是否发射”“是否回收”会在约束中产生复杂的耦合关系。所以我在复现时采用了论文中的三元变量思路x[i][j]表示卡车是否从节点 i 行驶到节点 j这是经典 TSP 变量。y[i][k][j]表示无人机是否从节点 i 发射服务客户 k然后在节点 j 回收。三元变量 y[i][k][j] 的设计非常巧妙一个变量同时编码了三个信息起飞点、服务对象、回收点。这样后面的时序约束写起来非常自然而且方便在 Gurobi 中直接通过索引访问。变量规模上如果客户数为 n则 x 变量数是 (n2)^2 级别y 变量数是 (n2)^3 级别。客户数 15 时大约只有几千个变量Gurobi 处理起来绰绰有余。但如果客户数到 100变量数会跳到百万级精确求解基本不现实这也是后面要说的模型规模局限性。为了建模我把仓库拆成两个节点节点 0 是出发点节点 N-1 是终点副本。这样卡车路径从 0 出发最终回到 N-1不需要额外处理“离开仓库再回仓库”的环逻辑。3.2 约束条件逐条拆解整个模型我拆成了五组约束每一组解决一个层面的问题。第一组客户服务约束。每个客户要么被卡车访问要么被无人机服务且只能被一种方式服务一次sum(x[i][k] for i ! k) sum(y[i][k][j] for i ! k, j ! k) 1这里需要注意如果客户 k 被卡车访问那么卡车的路径中必须有一条边进入 k即存在一个 i 使得 x[i][k] 1。如果 k 被无人机服务则存在 i 和 j 使得 y[i][k][j] 1。第二组卡车路径流守恒。卡车从 0 出发一次进入终点 N-1 一次中间每个客户节点进一次出一次sum(x[0][j] for j ! 0) 1 sum(x[i][N-1] for i ! N-1) 1 sum(x[i][j] for j ! i) 1 # 每个中间节点出度为1 sum(x[j][i] for j ! i) 1 # 每个中间节点入度为1这一组约束保证了卡车的路径是一条从起点到终点的简单路径不会出现子环路。第三组无人机发射回收节点合法性。如果无人机从 i 发射并在 j 回收即 y[i][k][j] 1那么卡车必须先访问 i 和 j。用逻辑约束实现sum(x[i][j] for j) y[i][k][j] # 卡车经过 i sum(x[j][i] for i) y[i][k][j] # 卡车经过 j这个约束的意思是无人机只能在卡车停靠的节点上进行发射和回收不能在空中随便发射。第四组时序约束。这是最核心的部分。定义变量 t[i] 为卡车到达节点 i 的时间那么卡车从 i 到 j 的行驶时间必须被传播t[j] t[i] dist[i][j] / v_truck - M * (1 - x[i][j])无人机行动的时间线无人机在 i 发射飞行到 k再飞到 j中间需要花费(dist[i][k] dist[k][j]) / v_drone s其中 s 是发射/回收的固定时间。无人机到达 j 之前卡车不能先离开 jt[j] t[i] (dist[i][k] dist[k][j]) / v_drone s - M * (1 - y[i][k][j])这两个时序约束用到了大 M 法。大 M 法的本质是“当 x 1 时约束生效当 x 0 时约束自动满足”。因为 M 是一个足够大的数不等式右侧减去一个超大数后几乎必然成立约束就形同虚设。第五组无人机续航约束。任何一个无人机行动的总飞行时间不能超过最大续航 E(dist[i][k] dist[k][j]) / v_drone s E M * (1 - y[i][k][j])同样用大 M 法处理当 y 0 时右侧很大约束松弛当 y 1 时右侧就是 E约束严格生效。3.3 目标函数与模型整体目标函数很简单最小化车队完成所有任务的时间即卡车到达终点 N-1 的时间minimize t[N-1]把所有变量、约束和目标拼起来就是一个完整的 MIP 模型。给一个直观的类比这就像一个双人接力赛卡车是跑大环线的选手无人机是负责绕小圈串场的选手两者的节奏要严格匹配无人机飞出去必须等到卡车到了某个汇合点才能继续下一步目标就是让最后一个到终点的人用时最短。4. Python Gurobi 代码实现与结果分析4.1 数据准备我选择了 15 个客户坐标范围 [0, 100] 的二维平面仓库固定在坐标 (50, 50)。距离矩阵使用欧氏距离这是一个合理的简化因为无人机直线飞行自然没有问题卡车在学术模型中也常被简化成平面直线运动。import numpy as np n_customers 15 rng np.random.default_rng(42) # 节点从0开始0是仓库起点最后一个是仓库终点副本 coords rng.uniform(0, 100, size(n_customers 2, 2)) coords[0] [50, 50] coords[-1] [50, 50] dist np.zeros((n_customers 2, n_customers 2)) for i in range(n_customers 2): for j in range(n_customers 2): dist[i, j] np.linalg.norm(coords[i] - coords[j])参数设置上卡车速度 v_truck 1无人机速度 v_drone 2无人机续航 E 25时间单位发射/回收固定时间 s 1。注意这些参数对结果影响很大。无人机速度越快、续航越长它能够服务的客户比例就越高。我后续会在实验部分单独分析参数敏感性。4.2 核心代码展示完整代码贴出来太长这里只展示模型构建的核心片段完整的复现脚本可以在文末说明中获取。import gurobipy as gp from gurobipy import GRB N n_customers 2 # 节点总数含起终点 M 1e4 # 大M值 model gp.Model(FSTSP) # 决策变量 x model.addVars(N, N, vtypeGRB.BINARY, namex) y model.addVars(N, N, N, vtypeGRB.BINARY, namey) t model.addVars(N, lb0, vtypeGRB.CONTINUOUS, namet) # 1. 客户服务约束 for k in range(1, N - 1): model.addConstr( gp.quicksum(x[i, k] for i in range(N) if i ! k) gp.quicksum(y[i, k, j] for i in range(N) for j in range(N) if i ! k and j ! k and i ! j) 1, namefvisit_{k} ) # 2. 卡车路径流守恒 model.addConstr(gp.quicksum(x[0, j] for j in range(1, N)) 1, namedepot_out) model.addConstr(gp.quicksum(x[i, N - 1] for i in range(N - 1)) 1, namedepot_in) for i in range(1, N - 1): model.addConstr(gp.quicksum(x[i, j] for j in range(N) if i ! j) 1, namefout_{i}) model.addConstr(gp.quicksum(x[j, i] for j in range(N) if i ! j) 1, namefin_{i}) # 3. 无人机发射回收节点合法性 for i in range(N): for k in range(1, N - 1): for j in range(N): if i ! k and j ! k and i ! j: model.addConstr(gp.quicksum(x[i, jj] for jj in range(N) if i ! jj) y[i, k, j], namefvalid_launch_{i}_{k}_{j}) model.addConstr(gp.quicksum(x[ii, j] for ii in range(N) if ii ! j) y[i, k, j], namefvalid_return_{i}_{k}_{j}) # 4. 时序约束 for i in range(N): for j in range(N): if i ! j: model.addConstr(t[j] t[i] dist[i][j] / v_truck - M * (1 - x[i, j]), nameftime_truck_{i}_{j}) for i in range(N): for k in range(1, N - 1): for j in range(N): if i ! k and j ! k and i ! j: drone_time (dist[i, k] dist[k, j]) / v_drone s model.addConstr(t[j] t[i] drone_time - M * (1 - y[i, k, j]), nameftime_drone_{i}_{k}_{j}) model.addConstr(drone_time E M * (1 - y[i, k, j]), namefrange_{i}_{k}_{j}) # 目标函数 model.setObjective(t[N - 1], GRB.MINIMIZE) # 求解 model.optimize()这段代码基本就是上一节数学模型的一比一翻译。如果你之前没有接触过 Gurobi 的 Python API重点关注addVars和addConstr两个方法前者一次性创建变量字典后者逐个添加约束。用gp.quicksum来求和比 Python 自带的sum在效率上更优。4.3 运行结果与可视化分析求解完成后我们关心两件事最优总时间是多少、无人机服务了哪些客户、卡车路径长什么样。通过model.ObjVal获取最优目标值即车队总完成时间。再遍历 x 变量获取卡车路径遍历 y 变量获取无人机行动列表然后用 matplotlib 画图。我这组随机数据上15 个客户的最优总时间是大约 260 个时间单位无人机服务了 4 个客户。相比纯卡车 TSP约 320 个时间单位总时间缩短了约 18.75%。这个数字非常符合论文里的结论——无人机不是万能的但确实能带来可观的时间收益。从路径图上看无人机服务的客户通常是离主路较远的“偏远客户”或者与主路线形成三角形位置的节点卡车完全不用绕路去接它们。这个微观行为在图上非常直观建议读者画出来观察一下能加深对模型的理解。5. 复现踩坑实录与性能调优5.1 复现中容易翻车的 5 个问题第一个坑是节点编号错位。仓库拆成起点和终点后距离矩阵的维度变成了 n_customers 2如果客户编号还是从 1 开始循环范围稍微写错就会出现某个客户永远不被服务的情况。排查技巧是求解前先打印x变量中所有值为 1 的边人工检查路径是否连接了每个客户。第二个坑是大 M 值设得过大。我一开始为了图省事设了 M 1e6结果 Gurobi 报出数值警告求解速度明显下降甚至出现次优解误判为最优解的情况。后来改成 M 1e4问题立刻消失。大 M 不是越大越好够用就行最好根据问题规模动态计算max_dist np.max(dist) M max_dist / v_truck 2 * max_dist / v_drone 2 * s这个值是卡车绕行全图的最长时间上界加上无人机两个最远距离飞行时间保证松弛条件下约束一定失效同时不过分庞大。第三个坑是无人机时序约束中的“同时到达”问题。我最初写的约束只考虑了无人机飞行时间忘了加固定发射/回收时间 s导致模型中无人机可以做到“零时间转移”解出来明显不符合实际。建议在做任何论文复现时先把论文里的参数表看清楚FSTSP 这篇里 s 的典型取值是 1 到 5 个时间单位。第四个坑是启动求解后长时间不收敛。小规模算例没事客户数超过 30 之后 Gurobi 的求解时间会快速增长。解决方案是给求解器设置合理的 TimeLimit 和 MIPGapmodel.Params.TimeLimit 300 model.Params.MIPGap 0.01这样即使没有在时限内找到最优解也能得到一个可接受误差范围内的可行解。第五个坑是可视化时无人机路径被画成“直线穿过客户”。真实场景中无人机飞到客户上方需要减速悬停但模型假设的是瞬时完成服务。画图时要注意把无人机路径和卡车路径用不同颜色区分不然容易误导读者。5.2 大M、MIPGap 与求解时长调优实录为了验证模型规模和求解器的极限我分别测试了 15、30、50 个客户三种规模。客户数变量数约求解耗时最优 Gap155 千2 秒0%303.3 万120 秒0%5014 万600 秒超时3.5%50 个客户的算例即使给了 10 分钟Gurobi 仍无法证明最优最终停在了 3.5% 的 gap 上。这说明 FSTSP 的精确求解规模上限基本就在 30-50 个客户之间。如果你想做更大规模就必须转向启发式或元启发式算法。MIPGap 的设置在学术复现中也很讲究。论文里说“最优”你得先搞清楚自己验证的是不是真正的最优解。我建议至少设置MIPGap0.001也就是 0.1% 的误差范围并且记录求解器给出的下界和上界方便在复现报告里说明可信度。5.3 无人机参数敏感性什么时候收益最大复现完模型后我顺手做了一组简单的参数敏感性分析想弄清“无人机在什么条件下最有用”。固定 15 个客户分别改变无人机速度和续航时间无人机速度 1.5 倍卡车速度时总时间缩短约 10%。无人机速度 2 倍卡车速度时总时间缩短约 18%。续航从 15 提升到 30无人机服务客户比例从 2 个提升到 5 个总时长进一步缩短。当续航超过一定阈值后继续增加续航带来的收益开始下降因为卡车路径本身的长度已经构成瓶颈。这个结果给了一个非常朴素的业务判断你不需要一架能飞一个小时的无人机只要续航能覆盖比最远客户距离稍远一些的范围就能拿到大部分收益。对于真正做物流规划的朋友来说这个参数分析比单纯跑通代码更有实际价值。写在最后的复现心得这次复现 FSTSP 给我最大的体会是论文里的公式再复杂拆成变量、约束、目标三个层面后每个部分都只是“在描述一个业务规则”。比如“无人机续航不能超过 E”这条约束翻译成代码就是一个小小的线性不等式关键是搞懂为什么需要它、它会不会跟其他约束产生冲突。如果后续你想在这个方向深入我建议在目前的精确求解基础上尝试用 Gurobi 的回调函数实现一个简单的分支切割算法或者干脆转向遗传算法、模拟退火这类启发式方法把客户规模推到 100 甚至 200。另外一个很有意思的扩展方向是考虑多架无人机甚至是无人机与卡车“半路交接”而不是必须到节点会合——这些方向都已经有大量论文支撑跑通 FSTSP 后你会发现读那些论文时思路会清晰很多。