pymoo复现车辆与无人机协同配送路径优化：从建模到实验全流程解析

论文复现这活儿，难的不是读懂论文里的公式，而是把两页纸的数学模型变成能跑的代码，还要保证结果跟论文对得上。最近我花了不少时间啃一个很热门的方向——低空经济下的车辆与无人机协同配送路径优化，求解框架选的是Python生态里的pymoo，对应的求解模式是论文里强调的"集中式协同配送模式"。这篇博文就把整个复现链路拆开讲清楚：从问题背景、数学模型、pymoo选型，到编码设计、约束处理、实验解读，再到现在为止踩过的一些坑。如果你也在复现类似的VRP变体论文，或者想把车机协同配送引入自己的项目里，这篇应该能帮你少走不少弯路。

1. 先说清楚问题：车辆与无人机协同配送到底在优化什么

1.1 无人机单独配送的明显瓶颈，逼出了"车机协同"方案

看低空经济相关的配送研究，绕不开一个现实问题：无人机单独飞着去送货，看着很美，但航程短、载重小、电池续航撑不起大范围的末端配送。单架无人机飞一个来回，可能只覆盖几公里，遇到稍微偏一点的郊区或者村镇客户，电量就成了硬约束。

所以论文里普遍采用一个思路：让车辆和无人机组成"编队"。车辆沿着实际路网移动，充当无人机移动的起降平台和补给站；无人机从车辆停靠点起飞，去服务周围不适合车辆进入或者距离车辆路径较远的客户，然后返回车辆的下一个汇合点。术语上管这个叫Vehicle-UAV Cooperative Delivery，中文通常是"车辆与无人机协同配送"。

我用一个生活化的类比来理解这事：卡车像一个移动基站加移动仓库，无人机像一个短跑运动员。运动员不可能从基地一路跑到几十公里外去送东西，但可以跟着车到附近，再冲刺一小段完成投递。这样车辆弥补无人机的续航短板，无人机弥补车辆在灵活性上的不足，整体配送范围和单位成本都会更优。

1.2 集中式协同配送模式，和常见的"主从式"有什么本质区别

标题里"集中式协同配送模式"这几个字，论文里是重点，复现时更是决定了你整个代码架构的方向。

主从式模式的做法是两阶段决策：先用标准的车辆路径规划方法，算出一条相对合理的车辆主路径，然后再贪心地或者按规则把一部分客户分配给无人机，在车辆路径的某些节点上插入无人机的起飞和回收动作。这种方法计算快、解空间小，但存在一个结构性缺陷——车辆路径的规划过程中完全没有考虑无人机起降的便利性。可能车辆路径经过了一个非常不适合放飞无人机的位置，无人机绕了很远才完成任务，那整个方案的质量天花板就被定死了。

集中式协同配送模式则不一样，它是把车辆路径、客户分配（谁由车服务、谁由无人机服务）、无人机起降点、车辆等待时间、时间同步全部放进同一个数学模型里统一优化。也就是说，车辆路径和无人机调度是一起长出来的，不是接起来的。这样解的质量上限更高，代价是决策变量种类多、解空间爆炸、求解难度大幅上升。

我在复现前也对比过两种模式，用小规模算例实测下来，集中式在总成本上通常能压掉主从式8%到15%的成本（随算例特征浮动），代价是计算时间翻了几倍甚至十几倍。论文里喜欢强调"集中式"，卖点就是这个。

1.3 问题的本质：一个多目标优化的VRP变体

把这个问题往运筹优化的框架里一放，本质上是经典车辆路径问题（Vehicle Routing Problem, VRP）的无人机版本，可以叫VU-VRP（Vehicle-UAV Routing Problem）或者TSP-UAV这类变种。经典VRP的目标是最小化总行驶距离或者总成本，但车机协同配送天然带有多个互相冲突的目标：

• 总作业成本尽量小：车辆油耗、无人机能耗、时间成本都折算成钱。
• 总完成时间尽量小：从所有车辆和无人机出发到所有客户服务结束的最晚结束时间，也叫makespan。
• 无人机的投入使用量尽量小：无人机数量越多成本越高，空域协调越复杂。

这几个目标经常互相打架，比如你想缩短总完成时间，就要多用无人机并行作业，但成本会上升。所以论文里普遍用多目标进化算法来解，pymoo里的NSGA-II、NSGA-III这些算法就是干这个的。这个"多目标"属性，也是我后面决定用pymoo而不是自己硬写一个单目标求解器的关键原因。

2. 数学模型怎么建：目标函数、决策变量与关键约束

2.1 决策变量的表示方式：客户分配与车辆路径

复现的第一步，是把论文描述的模型转换成代码能处理的决策变量。我复现的这篇论文采用了一种常见表示方式：

第一组决策变量是客户分配，对每个客户点，决策它由车辆服务还是由无人机服务。第二组决策变量是车辆服务客户的访问顺序，本质是一个排列。第三组是无人机的起降点归属，也就是无人机从哪个车辆停靠点起飞、回收到哪个停靠点。

但pymoo的NSGA-II默认不能直接处理"排列+0-1混合"的决策变量，它内部主要依赖实数编码的交叉变异算子。所以我在实现时做了一个转换：每个个体仍然是一串实数向量，长度等于客户数的两倍左右，前半部分客户的取值表示服务方式偏好（比如大于0.5倾向无人机服务），后半部分的取值通过排序映射成一个访问顺序。这种做法在论文复现里很常见，牺牲了一点点编码直观性，换来了与pymoo原生算子兼容。

复现时一定要在代码注释里把这个映射关系写清楚，不然隔一周再回来跑代码，你绝对看不懂当初这个向量是怎么变成路径的。

2.2 目标函数的构成：行驶成本、无人机成本与等待惩罚

我的代码里把目标函数设计成两个输出，对应pymoo的n_obj=2：

第一个目标是总成本，包含三块。车辆行驶成本：单位距离成本乘以车辆路径总距离。无人机飞行成本：单位飞行成本乘以无人机总飞行距离，有些论文会用包含载重、速度的能耗模型来算，更贴近真实物理；我复现的论文简化成了线性折算，先跑通结果再谈精细化。时间相关惩罚成本：如果车辆需要在某个节点等待无人机返回，这段时间也要折算成机会成本，等待时间越长，整体效率越差。

第二个目标是最晚完成时间，也就是makespan，包括车辆行驶时间、无人机飞行时间和车辆等待无人机的时间总和的最大值。

这里有一个我在实验里感受很深的点：如果两个目标量纲差太多，比如总成本动不动上万，makespan才几十，那帕累托前沿很容易被成本目标支配，看起来几乎所有解都"差不多"。遇到这种情况，建议在目标函数里做归一化或者加权重，保持两个目标的量级接近。论文里可能不会提这种细节，但实际跑出来的前沿图好不好看、信息量大不大，往往就取决于这种小处理。

2.3 容易被忽略的约束：续航、载重、起降点同步

约束条件才是这个模型真正"劝退"人的地方，因为每类约束都对应代码里的一堆if判断。我把它们分成三类：

续航约束。无人机每趟飞行距离不能超过最大续航里程，也就是从起降点飞到客户再返回的距离要小于电池可支持的最大航程。这里注意论文里说的是"每趟"还是"累计航程"，两种读法差异很大，累计航程会严格得多。我复现的论文用的是单次最大航程，理由是无人机每次返回车辆后可以更换电池。

载重约束。客户的需求量不能超过无人机最大载重，否则只能分配给车辆服务，或者说需要拆分订单。拆分订单会极大增加模型复杂度，因为同一个客户被拆分后，无人机和车辆的路径耦合更紧，代码里推荐先不加拆分假设。

起降点时间同步约束。这是最麻烦、最容易被忽略的一类。无人机起飞后，车辆不能马上开走，否则无人机找不到回收点；反过来，车辆开到了汇合点，无人机还没到，车辆就要等。为了处理这个约束，我给车辆引入了等待时间变量，并且要求无人机回到车辆的时刻等于车辆到达汇合点时刻加上等待时间。这个约束直接拉长了计算时间，因为它让车辆路径和无人机路径强耦合，任何一方的调整都会影响另一方的可行性。

我拿一个小算例说明这份约束的实际效果：假设配送中心O出发，3个客户A、B、C，车辆速度40km/h，无人机速度60km/h，无人机最大续航8km。车辆被分配服务A和C，路径是O-A-C-O，无人机服务B，从车辆行驶到A之前的某个停靠点起飞。代码里判断的逻辑是：先计算车辆到这个停靠点的时刻，再计算无人机飞到B再返回停靠点的往返距离是否小于8km，以及返回时刻是否晚于车辆到达时刻。如果晚，车辆需要等待，这个等待时长会同时进入成本和makespan的计算。就这么一个看似简单的判断，嵌套进NSGA-II的几千次评估里，就把计算量拉上去了。

3. 为什么选pymoo：多目标求解器选型逻辑与实践验证

3.1 pymoo帮我省掉了什么工作量

多目标进化算法的原理，通俗说就是"模拟自然选择，在多目标空间里同时逼近多个互不支配的最优解"。原理不复杂，但真让你从零手写NSGA-II，快速非支配排序、拥挤度距离计算、锦标赛选择、模拟二进制交叉、多项式变异、精英保留，再加上约束处理机制，工程量相当可观，而且特别容易出隐蔽bug。

pymoo最大的价值在于把这些算法组件全部封装好了。我需要做的只是定义一个继承自pymoo.Problem的类，把目标函数写进_evaluate方法，然后指定算法（比如NSGA-II）、种群大小、终止条件，剩下的事情pymoo自己处理。

对于VU-VRP这种定制化程度非常高的模型，pymoo还能自定义算子，虽然这需要理解它的Operator抽象机制，但总比自己连排序算法都重写要轻松。

顺便说一句，pymoo自带帕累托前沿可视化、超体积指标计算、历史记录抓取，这些功能在做论文复现实验时几乎全都能用上。我后来画帕累托前沿图，直接用pymoo的接口就能输出，省了很多事。

3.2 求解器选型的横向对比，为什么不是DEAP也不是手写

我做一个选型对比之前，先列了四个候选方案：pymoo、DEAP、jMetal、完全手写。表格如下：

最终的结论是pymoo。理由也很简单：这里是多目标问题，NSGA-II直接可用，同时pymoo约束处理用的是约束违反量的方式，跟我的模型天然契合。DEAP不是不好，而是当我需要约束支配关系、HV指标、Termination机制时，DEAP提供的封装级别太低了，我得自己拼装，这在论文复现场景里不是性价比最高的选择。

3.3 pymoo版本与API演进，复现论文时最常见的兼容性坑

这里必须提醒一句：pymoo从2.x到3.x，API发生了非常大的变化，很多2022年以前的论文代码用的是老写法。比如老版本用add()方法添加惩罚函数、老式的Minimize入口，新版本全部改成了继承Problem类、用pymoo.optimize.minimize函数。如果你拿老代码直接换新版pymoo跑，大概率会报出一堆AttributeError或者TypeError。

我自己的处理方式：不看论文附件的老代码，直接按pymoo 3.x的语法重写调用层，核心模型逻辑保留。这样既省去来回改API的时间，也能确保和当前pymoo版本兼容。建议你也这么做，别在旧版本的泥潭里挣扎。

4. 核心代码实现：问题定义、编码设计与算子配套

4.1 继承Problem类的骨架：_evaluate方法怎么写

pymoo里定义一个优化问题的标准方式是继承Problem类，然后实现_evaluate方法。我的VUVRP类骨架长这样：

import numpy as np from pymoo.core.problem import Problem class VUVRP(Problem): def __init__(self, data): self.data = data self.n_customer = data["n_customer"] # 决策变量：客户数（分配偏好） + 客户数（车辆路径顺序映射） n_var = self.n_customer * 2 super().__init__( n_var=n_var, n_obj=2, n_ieq_constr=6, # 约束数量根据实际需要调整 xl=0.0, xu=1.0, ) def _evaluate(self, X, out, *args, **kwargs): F = np.zeros((X.shape[0], 2)) G = np.zeros((X.shape[0], self.n_ieq_constr)) for i, x in enumerate(X): allocation, vehicle_route = self._decode(x) cost, makespan, violations = self._evaluate_solution(allocation, vehicle_route) F[i, 0] = cost F[i, 1] = makespan G[i, :] = violations out["F"] = F out["G"] = G

这段代码看起来简单，真正复杂的是_decode和_evaluate_solution两个函数。它们把决策变量解码成实际路径，然后模拟整个配送过程，计算目标函数和约束违反量。

4.2 编码策略：从客户分配到车辆路径的转换函数

_decode函数是核心。我的实现思路是：

• 决策向量X的前n_customer个元素表示客户的分配偏好。取值大于0.5的客户倾向于由无人机服务，小于等于0.5的由车辆服务。
• 后n_customer个元素用来生成车辆服务的客户顺序。做法是对这一部分向量做argsort，得到一组客户编号的排列，然后过滤掉已经分配给无人机的客户，剩下的顺序就是车辆服务的顺序。
• 无人机的起降点通过一个邻域匹配矩阵确定，简单来说就是每个无人机服务客户，找一个离它最近的车辆停靠点作为起飞/回收点。

这里有一个实操经验：纯随机初始化解里，会有大量个体把很多客户分配给无人机，导致续航约束严重违反，初始种群基本全是不可行解。解决办法是写一个可行性恢复启发式：如果某个客户被分配给无人机但往返距离超过最大续航，就强制把它的分配偏好改成车辆服务。这个启发式不改变NSGA-II的搜索逻辑，只是把初始解空间拉到合理范围内，能明显加快收敛。

4.3 约束违反量的计算与pymoo的约束处理机制

pymoo处理约束的机制是：你需要在out["G"]里返回每个约束的违反量，违反量大于0代表不可行，等于或小于0代表可行。NSGA-II在比较两个解时，先看可行性，如果都可行再看支配关系。

所以我的_evaluate_solution函数里，每个约束写成"实际值减允许值"的形式：

violations = [] # 续航约束：航程不能超过最大续航 violations.append(drone_route_distance - max_flight_range) # 载重约束：无人机服务客户需求量不能超过载重 violations.append(drone_demand - drone_capacity) # 时间同步约束：车辆等待时间不能超过允许上限 violations.append(vehicle_wait_time - max_wait_time) # 其他约束类似...

这里有一个非常隐蔽的坑：浮点误差。比如某个解的实际航程是7.9999999，续航上限是8.0，理论上可行，但因为浮点精度被判成不可行。解决方案很简单，所有约束违反量在返回前统一减去一个很小的epsilon，比如1e-6，给边界解留一点余地。这种问题不实际跑代码根本发现不了，但一旦出现，会污染整个种群，让帕累托前沿缺一大块。

4.4 交叉、变异的配套设定和参数选择

模型定义好之后，算法配置方面我踩过不少坑。pymoo默认给实数决策变量用SBX交叉（模拟二进制交叉）和多项式变异，这在一开始跑VUVRP时效果还行，但明显对路径顺序部分的改进效率不高，因为顺序信息靠argsort还原，SBX在这个映射下破坏性比较强。

我最后的做法是混合算子设计：分配偏好部分用SBX，路径顺序部分我自定义了一个类似顺序交叉（Order Crossover）的算子，只对后半段向量作用。这样路径结构能更完整地遗传给后代，不会每次都被交叉打乱重排。

参数上，我并没有盲目信任论文里的默认值，而是跑了小规模预实验。小规模算例（10个客户）种群大小80，进化200代，NSGA-II的交叉概率0.9，变异概率0.1，邻居参数eta设为20。这个组合在我的数据上表现相对稳定。这里也提醒一下，参数不是通用的，换一批坐标分布更密的客户点，你可能需要重新标定。

5. 实验设置与结果解读：从帕累托前沿看方案的取舍

5.1 小规模算例的实验设计与参数表

为了验证复现代码是否正确，我先设计了3组小规模算例：10个客户、20个客户、50个客户。每组算例里车辆数都是1辆（对应论文场景），无人机数量取1架和2架两种情况，客户需求值从1到5随机生成，坐标在10km乘10km的区域内随机撒点。

表：小规模算例参数设置

每个算例独立重复运行10次，收集帕累托前沿、非支配解数量、超体积指标（HV）和平均计算时间。这样做是为了减少NSGA-II随机性对结论的影响，论文复现阶段更要对数字负责。

5.2 帕累托前沿图怎么读：两个目标之间的权衡

跑完10个客户的小算例后，我画了两目标的帕累托前沿图。横轴是总成本，纵轴是最晚完成时间。典型的前沿形状是一条往左下角凸出的曲线：最左边的点对应成本最低的方案，通常是把大量客户分配给车辆，无人机用得很少，车辆慢慢跑；最下边的点对应完成时间最短的方案，车辆某些路径会让位给无人机，无人机满负荷飞行。

一个有意思的现象是：当无人机数量从1架增加到2架时，帕累托前沿整体向坐标轴左下方向收缩。这说明增加无人机能在压缩完成时间的同时并不显著抬高总成本，在这个模型假设下是"帕累托改进"。但注意，这不等同于无人机越多越好，因为模型里没计购置成本和空域管理成本，真正做工程决策时，这些成本应该作为第三个目标或参数进入模型。

读这张图时，我还注意到一个细节：前沿最左端和最下端的解数量明显偏少，中间解聚集较多。这是NSGA-II常见的现象，极端解不容易被保留，因为中间解的拥挤度距离更小，在同等支配关系下更有生存优势。论文里如果展示这样的前沿，说明算法设置基本正常。

5.3 与主从式模式的对比：集中式到底赢在哪里

为了验证"集中式协同配送模式"的价值，我实现了一个主从式baseline：先用单目标VRP求解器求出车辆最短路径，再用贪心算法把可受益的客户分配给无人机。结果整理成下面这张表：

从表里能清楚看到，集中式在成本和完成时间两个维度都压过了主从式。原因也不难理解：主从式把车辆路径先定死，无人机只能被动适应；集中式让车辆路径主动为无人机的起降创造便利条件，比如车辆愿意稍微绕一点路去一个更好的起飞点，换来无人机来回程短，整体反而更优。

但我也必须诚实记录一个代价：集中式求解时间远大于主从式，10客户算例有时候要2到3分钟，而主从式几秒钟就出结果。这个性价比问题在工程落地时非常关键，论文里通常不会明说，但你如果要做实时调度，动辄几分钟的求解时间是不能接受的。

6. 复现中踩过的坑与下一步改进方向

6.1 踩坑实录：pymoo版本、浮点约束、初始种群、性能瓶颈

第一个坑是pymoo版本混乱。我一开始在conda环境里装到了pymoo 2.x，跑新版代码报错无数，问题都集中在Population、Termination这些接口上。后来干脆新建环境，直接装3.x，并锁定在requirements.txt里，再也不折腾。

第二个坑是浮点约束，也就是前文提到的1e-6问题。第一次跑出来的帕累托前沿只有寥寥几个点，原因是大量在续航边界附近的解被浮点误差误杀。加了epsilon之后，前沿才恢复正常密度。

第三个坑是初始种群质量。纯随机生成的解里有90%以上违反无人机续航约束，NSGA-II前期绝大部分计算都在"修复不可行解"而不是"优化可行解"。引入最近邻启发式初始化之后，前期收敛速度肉眼可见地变快了。

第四个坑是性能瓶颈。当客户数到50个的时候，每代评估500个解，每个解都要解码并计算距离，纯Python循环跑得非常酸爽，200代等了好几分钟。我优化了两步：第一步，用一个预计算的客户间距离矩阵，解码时不重复调用距离公式；第二步，尽量把不涉及分支判断的计算用NumPy批量完成，只有必须逐解判断的地方才留循环。优化之后时间降到原来的三分之一左右。

6.2 大规模算例的计算瓶颈与加速思路

当客户数超过50个之后，pymoo默认的NSGA-II开始吃力了。我试验了几种加速思路，最有效的是限制最大评估次数而不是固定代数，让算法在预算用完后停止，再配合U-NSGA-III这种更适合高维目标和小种群的算法。另一个思路是给算法装上2-opt局部搜索作为增强算子，每进化几代就对当前非支配解做一轮局部改进，相当于"算法内置爬山"，效果很好但实现成本也不低。

如果项目要处理几百个客户，我的建议是冷静一点，别指望单个NSGA-II硬扛。更现实的做法是先对客户做空间聚类，把大问题拆成几个子区域，分别求解后再拼接。论文里可能不做这一步，但工程化时必须考虑。

6.3 后续可以扩展的方向：动态需求、随机时间、空域约束与多层网络

最后聊聊这个模型往工程落地走的扩展空间。我接触到的一些项目资料里，已经把静态模型往动态方向改了：客户订单不是一次性全部知道，而是随时间逐步到达，这时候需要把确定性优化改成滚动时域优化，每过一段时间用新数据重新求解一次。

另一个方向是随机行驶时间的鲁棒优化。现实里车辆速度、无人机风速都不可控，论文里假设匀速就是理想情况了。要贴近实际，得把行驶时间建模成随机变量，用期望值或者最坏情况作为目标，或者加分布鲁棒约束。

低空经济真正落地时还绕不开空域约束，比如禁飞区、噪音敏感区，甚至多无人机同时作业的安全间隔，这些约束在模型里可以表示成无人机飞行路径禁止进入某些多边形区域，或者最小安全距离约束。这是很实际的方向，但模型复杂度会再上一个台阶。

还有一个我觉得很有潜力的方向是多级配送网络：无人机从区域中转站起飞送到社区驿站，车辆从城市大仓运到中转站，两级协同优化。这种网络比单一车机编队更贴近真实物流体系，也把低空配送的价值从"最后一公里"扩展到了"最后三公里"。把这篇复现的模型理解透了，再看这些方向，你会发现它们本质上都是"时间同步+容量约束+成本权衡"这套思路的变体。

我个人在实际操作中的一个体会是：复现这类论文，最大的收获不是跑通了别人给的代码，而是建立了一套"从问题建模到算法配置再到结果调优"的闭环能力。做车辆与无人机协同配送项目，别一上来就挑战上百个客户的大规模算例，那样你很难判断结果是算法的还是模型的。先跑通3到5个客户的微型算例，手工把最优解算出来，拿元启发式的结果去对标，确认模型本身没问题，再逐步扩大规模。这一步做扎实了，后面的一切都是水到渠成。