企业官网建设流程全解析

复现论文这件事，最怕的就是对着公式看了一整天，代码写出来却跑不通。上周我终于把 FSTSP（The 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.9，Gurobi 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(vtype=GRB.CONTINUOUS, name="x") y = m.addVar(vtype=GRB.CONTINUOUS, name="y") 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 变量数是 (n+2)^2 级别，y 变量数是 (n+2)^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 之前，卡车不能先离开 j：

t[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, vtype=GRB.BINARY, name="x") y = model.addVars(N, N, N, vtype=GRB.BINARY, name="y") t = model.addVars(N, lb=0, vtype=GRB.CONTINUOUS, name="t") # 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, name=f"visit_{k}" ) # 2. 卡车路径流守恒 model.addConstr(gp.quicksum(x[0, j] for j in range(1, N)) == 1, name="depot_out") model.addConstr(gp.quicksum(x[i, N - 1] for i in range(N - 1)) == 1, name="depot_in") for i in range(1, N - 1): model.addConstr(gp.quicksum(x[i, j] for j in range(N) if i != j) == 1, name=f"out_{i}") model.addConstr(gp.quicksum(x[j, i] for j in range(N) if i != j) == 1, name=f"in_{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], name=f"valid_launch_{i}_{k}_{j}") model.addConstr(gp.quicksum(x[ii, j] for ii in range(N) if ii != j) >= y[i, k, j], name=f"valid_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]), name=f"time_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]), name=f"time_drone_{i}_{k}_{j}") model.addConstr(drone_time <= E + M * (1 - y[i, k, j]), name=f"range_{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 和 MIPGap：

model.Params.TimeLimit = 300 model.Params.MIPGap = 0.01

这样即使没有在时限内找到最优解，也能得到一个可接受误差范围内的可行解。

第五个坑是可视化时无人机路径被画成“直线穿过客户”。真实场景中无人机飞到客户上方需要减速悬停，但模型假设的是瞬时完成服务。画图时要注意把无人机路径和卡车路径用不同颜色区分，不然容易误导读者。

5.2 大M、MIPGap 与求解时长调优实录

为了验证模型规模和求解器的极限，我分别测试了 15、30、50 个客户三种规模。

50 个客户的算例即使给了 10 分钟，Gurobi 仍无法证明最优，最终停在了 3.5% 的 gap 上。这说明 FSTSP 的精确求解规模上限基本就在 30-50 个客户之间。如果你想做更大规模，就必须转向启发式或元启发式算法。

MIPGap 的设置在学术复现中也很讲究。论文里说“最优”，你得先搞清楚自己验证的是不是真正的最优解。我建议至少设置MIPGap=0.001，也就是 0.1% 的误差范围，并且记录求解器给出的下界和上界，方便在复现报告里说明可信度。

5.3 无人机参数敏感性：什么时候收益最大

复现完模型后，我顺手做了一组简单的参数敏感性分析，想弄清“无人机在什么条件下最有用”。

固定 15 个客户，分别改变无人机速度和续航时间：

• 无人机速度 1.5 倍卡车速度时，总时间缩短约 10%。
• 无人机速度 2 倍卡车速度时，总时间缩短约 18%。
• 续航从 15 提升到 30，无人机服务客户比例从 2 个提升到 5 个，总时长进一步缩短。
• 当续航超过一定阈值后，继续增加续航带来的收益开始下降，因为卡车路径本身的长度已经构成瓶颈。

这个结果给了一个非常朴素的业务判断：你不需要一架能飞一个小时的无人机，只要续航能覆盖比最远客户距离稍远一些的范围，就能拿到大部分收益。对于真正做物流规划的朋友来说，这个参数分析比单纯跑通代码更有实际价值。

写在最后的复现心得

这次复现 FSTSP 给我最大的体会是：论文里的公式再复杂，拆成变量、约束、目标三个层面后，每个部分都只是“在描述一个业务规则”。比如“无人机续航不能超过 E”这条约束，翻译成代码就是一个小小的线性不等式，关键是搞懂为什么需要它、它会不会跟其他约束产生冲突。

如果后续你想在这个方向深入，我建议在目前的精确求解基础上，尝试用 Gurobi 的回调函数实现一个简单的分支切割算法，或者干脆转向遗传算法、模拟退火这类启发式方法，把客户规模推到 100 甚至 200。另外一个很有意思的扩展方向是考虑多架无人机，甚至是无人机与卡车“半路交接”而不是必须到节点会合——这些方向都已经有大量论文支撑，跑通 FSTSP 后，你会发现读那些论文时，思路会清晰很多。