最近在帮某课题组复现一篇关于“车辆与无人机协同配送路径优化”的论文,折腾了将近两周,把论文里那个“集中式协同配送模式优化问题”从公式一步步落到Python代码,最后用pymoo跑通了Pareto前沿。整个过程踩了不少坑,也积累了一些可以直接复用的代码框架、参数配置和调试思路,今天整理出来分享给同样在做运筹优化、物流配送路径规划,或者正在复现论文模型的读者。
这个项目的核心很好概括:在低空经济快速发展的背景下,末端配送不再只有卡车一种执行者,而是由卡车和无人机组成混合车队。无人机负责短距离、小批量、高时效的客户点,卡车承担长距离运输和中转。所谓“集中式协同配送模式”,就是所有卡车和无人机都由一个控制中心统一规划路线,目标是让整个车队完成配送任务的总成本最小。论文里通常把这个问题建模成一个混合整数规划模型,但用传统的精确求解器很难在大规模算例上收敛,所以需要引入启发式算法。我选用pymoo,是因为它内置了NSGA-II、NSGA-III等多个成熟的多目标进化算法框架,写自定义Problem、自定义解码器都很方便,适合做论文代码复现和二次开发。
如果你正在头疼“论文公式看得懂,但代码就是跑不出来”这个问题,这篇内容会重点拆解四大块:集中式协同配送的数学模型怎么建、pymoo怎么接入自定义优化问题、路径编码与解码怎么做、以及常见报错和调试技巧。我会尽量用白话把模型和代码讲透,同时给出可以直接改参数的Python代码。
1. 项目概述与背景:为什么值得复现这个模型
1.1 从标题拆出三个关键词
“低空经济”“车辆与无人机协同配送”“pymoo”,这三个词基本就能概括整篇论文代码复现的工作量。
低空经济指的是以无人机、低空飞行器为主要载体的新经济形态,落到物流配送上,最常见的形式就是无人机末端配送。前几年我们讨论无人机配送,更多是试飞、演示,真正进入学术研究和工程落地的,是它和传统地面车辆的协同问题。一辆卡车本身的续航长、载重大,但线路固定;无人机灵活、速度快,但航程短、载重小。两者结合,表面上是“物流效率提升”,本质上是一个组合优化问题:哪些客户点给卡车,哪些客户点给无人机,无人机从哪个点起飞、在哪回收,卡车路径怎么走。
“车辆与无人机协同配送路径优化”研究的是这类问题中的路径决策部分。论文题目里还有一个限定条件叫“集中式协同配送模式”,这意味着所有决策都由单一决策中心完成,不存在多主体博弈,是一个标准的全局优化问题。这对建模和求解来说反而是好事,复杂度可以完全集中在目标和约束建模上。
“pymoo求解”则是方法论层面的选择。pymoo是一个Python环境下的多目标优化库,提供NSGA-II、NSGA-III、U-NSGA-III、MOEA/D这类进化算法,使用者只需要把自己的优化问题定义为标准Problem类,就能直接调用内置算法求解。相比从零手写遗传算法,pymoo的代码可维护性和可扩展性都要好很多。
1.2 复现论文代码到底在复现什么
很多同学第一次接触论文复现,以为就是把论文里的公式变成代码,跑出一张差不多的图就完了。实际上论文代码复现包含四个层次的工作:
第一层是数学模型复现。把目标函数、约束条件、决策变量的定义用代码精确表达出来,任何一处符号理解错误,都会导致后续结果对不上。第二层是算法复现。论文如果用了改进的遗传算法、模拟退火或者混合算法,需要把算法的算子、参数、终止条件全部还原。第三层是数据层面的复现。使用论文提供的算例或者公开数据集,保证实验结果可以横向对比。第四层是结果的可视化和指标分析,用超体积、GD、IGD等指标评估算法的收敛性和分布性。
本次项目的难点主要集中在第一层和第三层。论文里的模型看着清楚,真正实现的时候,无人机续航约束、卡车容量约束、时间窗约束交织在一起,很容易出现索引混乱和约束漏写的问题。尤其是我在复现时还顺手扩展了模型,加入了无人机数量限制和发射/回收点约束,调试的复杂度进一步上升。
2. 数学模型与集中式协同配送模式拆解
2.1 集中式协同配送模式到底是什么
我用自己的话解释一下这个模式:有一个配送中心,中心里有若干辆卡车,每辆卡车可以搭载一架或多架无人机。车辆和无人机组成一个“车机编队”出发,卡车沿着一条主干路径前进,无人机在某些节点起飞,飞向附近的客户点完成配送,再飞回指定的回收点与卡车汇合。
这里的“集中式”强调的是:每辆卡车到底走哪条路、每架无人机服务哪些客户、在哪里起飞、在哪里回收,全部由统一计算得出,不存在实时讨价还价。这样做的好处是能拿到全局最优解,缺点是计算难度高。所以论文中常用启发式算法求解,而非分支定界这类精确算法。
我在复现时用了一个 5 个客户点的小规模案例来检验代码逻辑正确性,然后扩展到 30 个客户点的中等规模案例做性能测试。5 个客户点的场景下,车辆和无人机的分工比较直观,比如卡车从配送中心出发,路线是“中心–客户2–客户4–中心”,无人机在客户2起飞服务客户5,再到客户4汇合。这种手算能验证的结果,对排查代码bug特别重要,建议你也这样操作。
2.2 目标函数与约束条件怎么定义
论文中的集中式协同配送问题,通常可以建模为如下形式:
目标函数:
min f = Σ_{k} L_k (卡车路径总距离) + Σ_{u} L_u (无人机飞行总距离) + λ * U_used
其中 L_k 是卡车k的行驶距离,L_u 是无人机u的总飞行距离,U_used 是实际使用无人机架数,λ 是无人机启用成本系数。如果需要多目标优化,也可以把总距离和总完成时间拆成两个目标,此时就用多目标进化算法求解Pareto前沿。
约束条件我整理成表格:
| 约束类型 | 说明 | 建模建议 |
|---|---|---|
| 服务唯一性 | 每个客户点只能被卡车或无人机服务一次 | 分配变量二值化 |
| 路径连续性 | 卡车路径必须从配送中心出发并返回 | 流量守恒约束 |
| 无人机航程限制 | 无人机单次飞行距离不能超过最大航程E | 非线性约束,可在解码器内判断 |
| 无人机荷载限制 | 无人机单次携带包裹总重量不超过上限 | 需求求和约束 |
| 数量限制 | 实际使用的无人机数量不能超过可用数量 | 总架数约束 |
论文中常用所有客户点坐标已知、需求已知、卡车容量和无人机航程已知的前提。无人机采用“单次起飞服务多个客户之后回收”还是“每架次只服务一个客户”,决定了模型复杂度。我的复现采用后者,也就是发射–配送–回收的典型模式,这样便于和多数公开论文的结果对比。
2.3 决策变量与编码方式
数学模型中的决策变量一般有三组:卡车路径变量(卡车k是否从i行驶到j)、客户分配变量(客户i是否由无人机u服务)、时间变量(开始服务时间)。但在遗传算法的实现里,我不直接操作这些0-1变量,而是采用“随机密钥编码”。
随机密钥编码的核心思想是把一个复杂的组合决策压缩成一个固定长度的实数数组。数组的每个元素值在[0,1]区间内,解码器按照一定规则把它翻译成具体路径。这样做的好处很明显,适应度函数可以连续变化,方便遗传算子操作,避免大量非法解。缺点是解码器设计直接影响解的质量,编码到解码这一步是整个代码复现的关键难点。
2.4 为什么这道题不适合用精确求解器
刚接触这个问题时,我第一反应是试试商用求解器。模型虽然带有非线性项,但可以通过线性化处理转成MILP。但实测下来,当客户点数量超过20个时,分支定界的搜索空间会爆炸式增长,一台普通电脑上可能跑几个小时都没有可行解。而论文中通常要求30到100个客户点的算例,精确方法根本不现实。
遗传算法这类启发式方法虽然不保证全局最优,但能够在合理时间内找到质量不错的近似解,而且天然支持多目标优化。论文复现阶段,我们更关注的是“算法思路可复现”“参数可对比”,所以pymoo这套进化算法框架是最合理的选择。如果你想验证小规模案例,也可以用Gurobi或SCIP这类求解器去求精确解作为基准,但没必要作为主求解器。
3. 基于pymoo的求解框架搭建
3.1 为什么是pymoo而不是自己写遗传算法
我知道有些人习惯自己手写遗传算法,觉得可控性强。但手写GA有一个绕不开的问题:当你需要一个功能时,你会发现整个框架都得改。比如我最早加约束处理的时候,自己写的代码先判断合法性再给惩罚函数,逻辑是能跑,但换成另一个目标函数就要重写一大片。pymoo把这些问题封装好了,你只需要关注自己的问题定义。
pymoo的代码结构清晰,主要就几步:
from pymoo.core.problem import Problem from pymoo.algorithms.moo.nsga2 import NSGA2 from pymoo.operators.crossover.sbx import SBX from pymoo.operators.mutation.pm import PM from pymoo.optimize import minimizeProblem类负责定义优化问题,算法类负责进化搜索,minimize就是把两者组合起来跑。如果论文里用了自定义算子,pymoo也允许重载操作算子,不会限制扩展。
还有一点,pymoo内置了收敛性指标和Pareto前沿分析工具,像Hypervolume、IGD这些指标直接调用就行,不需要自己从头实现。
3.2 定义你的Problem子类
pymoo的使用核心是继承Problem类。你要做的就是在初始化函数里告诉框架这个问题的决策变量个数、目标函数个数、约束个数,然后在_evaluate函数里实现目标值和约束值的计算。
class CVDPProblem(Problem): def __init__(self, n_var, data): self.data = data super().__init__( n_var=n_var, n_obj=2, n_constr=2, xl=0.0, xu=1.0 ) def _evaluate(self, X, out, *args, **kwargs): # X shape为 (种群大小, 变量个数) f1 = np.zeros(X.shape[0]) f2 = np.zeros(X.shape[0]) g1 = np.zeros(X.shape[0]) g2 = np.zeros(X.shape[0]) for i in range(X.shape[0]): route_info = decode_individual(X[i], self.data) f1[i] = route_info["truck_distance"] f2[i] = route_info["drone_distance"] g1[i] = route_info["drone_range_violation"] g2[i] = route_info["capacity_violation"] out["F"] = np.column_stack([f1, f2]) out["G"] = np.column_stack([g1, g2])这里有个容易踩的坑:_evaluate里的X形状是二维数组,不是一维。我第一次写的时候习惯性当成单个个体处理,结果索引越界。pymoo的设计是一次性评估整个种群,这样能利用向量化运算提升速度。你在写代码时要注意循环内用row索引,避免修改全局变量。
3.3 算法与操作算子配置
我最终选择的是NSGA-II。原因很简单:问题规模适中,两个目标,非支配排序加上拥挤度距离已经够用。
algorithm = NSGA2( pop_size=120, sampling='real_random', crossover=SBX(prob=0.9, eta=15), mutation=PM(prob=0.2, eta=20), eliminate_duplicates=True ) res = minimize( problem, algorithm, termination=('n_gen', 300), seed=42, save_history=True, verbose=True )参数不是随便选的。种群大小120,对于30个客户点的算例来说不算大也不算小,太小容易早熟,太大单代计算太慢。模拟二进制交叉概率0.9,让大部分个体参与交叉,保持种群的多样性。多项式变异概率0.2,变异强度通过eta=20控制,变异过强会破坏已经收敛的结构,过弱会陷入局部最优。
这里建议你把终止代数设成300代以上。我一开始只跑50代,得到的前沿稀疏得一塌糊涂,后来发现论文里的收敛图都是400代左右才稳定。如果你用的是pymoo新版本,termination参数既可以传元组,也可以传Generation对象,注意区分。
4. 核心代码实现:从路径编码到目标评估
4.1 初始化数据与问题输入
为了便于复现,我用了一个伪随机生成算例。随机生成30个客户点坐标,生成方式比较简单,用numpy在二维平面上随机采样,单位距离设为公里。配送中心固定在坐标原点,需求在[1,5]公斤范围内随机生成。卡车容量50公斤,无人机最大航程10公里,无人机数量为3架。这些参数都可以在代码头部自定义。
import numpy as np import pandas as pd np.random.seed(42) n_customer = 30 coordinates = np.random.rand(n_customer, 2) * 20 demands = np.random.randint(1, 6, size=n_customer) depot = np.array([0.0, 0.0]) capacity = 50.0 drone_range = 10.0 drone_num = 3我用一个字典把这些参数存起来,传给CVDPProblem类。为什么不直接写在类里面?因为后续做灵敏度分析时,要反复改参数跑实验,数据驱动的方式更灵活。
4.2 随机密钥解码:从实数数组到可行路径
解码器是整个复现中最关键的一环。每个个体是一个长度为n_customer + drone_num的实数数组,前n_customer个值定义客户点的访问优先级,后drone_num个值定义无人机使用的倾向阈值。
解码逻辑是这样的:
- 根据前n_customer个值,按值大小排序得到客户访问序列。
- 逐个客户判断:如果某个客户的坐标点到当前卡车位置的距离超过无人机航程,或者客户需求超过无人机容量,则由卡车直接服务;否则,根据后段中的加权阈值决定是否分配给无人机。
- 被无人机服务的客户,从最近的卡车路径节点发射,无人机完成任务后返回最近的下一个卡车节点。
- 最终输出卡车路径序列和无人机发射回收位置表。
我把解码结果封装成一个字典,包含truck_route、drone_flights、truck_distance和drone_distance四个字段。这样目标函数评估时只需要读字典,逻辑非常清晰。
解码器还可以加约束判断:如果使用的无人机数量超过限制,直接在g1约束上记一个违反值。这个设计让惩罚不在目标函数里做,而是交给pymoo的约束支配机制处理,效果比简单的罚函数好很多。
4.3 目标函数与约束评估
在_evaluate函数里,我对每个个体跑一次解码器,得到卡车和无人机的距离后,把两个距离分别作为两个目标函数。为什么用两个目标而不是加权成一个?因为论文对比的通常是Pareto前沿,两个目标可以画出二维散点图,更直观地展示成本结构之间的权衡。
约束评估我做了两个:
第一个是无人机航程约束,检查每架无人机的飞行距离是否超过了drone_range,如果有超程则记录超出的最大值。这个值对应的g1为正表示违反,pymoo会用约束支配来处理。
第二个是卡车容量约束,检查每辆卡车的总载重是否超过capacity,同理记录超载量。这两个约束都是不等式约束,在pymoo里直接放在out["G"]里,注意G的每个元素表示违反程度,不违反时为0或负数。
4.4 求解循环与结果保存
minimize函数跑完之后,res.F里面保存的是Pareto前沿上每个解的目标值,res.X保存的是对应的决策变量。通常我会把结果存成CSV,然后用matplotlib画Pareto前沿图。
from pymoo.optimize import minimize res = minimize( problem, algorithm, termination=('n_gen', 400), seed=42, save_history=True, verbose=True ) pf = res.F pd.DataFrame(pf, columns=['truck_distance', 'drone_distance']).to_csv('pareto_front.csv', index=False) import matplotlib.pyplot as plt plt.scatter(pf[:, 0], pf[:, 1], alpha=0.7) plt.xlabel('Truck Distance (km)') plt.ylabel('Drone Distance (km)') plt.grid(True) plt.savefig('pareto_front.png', dpi=150)这里有个小细节,save_history=True能记录每一代的所有种群目标值,后面画收敛曲线时可以直接用history里的数据,不需要自己额外记录。history是个包含代数信息的对象列表,取用很方便。
5. 实验结果与灵敏度分析
5.1 小规模算例验证
我先把参数改成5个客户点、1辆无人机,手算一遍最优解。假设配送中心在原点,客户点坐标很简单,卡车绕一圈走三个客户,无人机负责另外两个。跑完代码以后,输出的truck_distance=12.3,drone_distance=4.1,和手算值完全一致,基本说明解码器和目标函数逻辑没有硬伤。
这一步非常关键。如果你一上来就跑30个客户点,错了根本不知道错在哪。用小规模案例做单元测试,是论文代码复现的标准做法。
5.2 算法性能指标怎么看
大规模案例下,我用HV指标评估Pareto前沿的质量。pymoo提供了计算超体积的工具:
from pymoo.indicators.hv import HV hv = HV(ref_point=np.array([60, 40])) value = hv(res.F)参考点要选得比前沿上所有点都大,否则HV值没有意义。在我的算例里,卡车距离通常在20到50之间,无人机距离在5到30之间,所以参考点设成了(60, 40)。HV值越高,说明前沿在目标空间中的覆盖面越大,解集质量越好。
除了HV,我还观察了收敛曲线。做法是把history里每一代种群的最小目标值取出来画折线图,可以看到大概在250代之后曲线变得平缓,说明算法已经收敛。
5.3 参数灵敏度实验
为了验证模型的响应是否符合直觉,我做了三组对比实验:
第一组改变无人机最大航程,从10公里提高到20公里,结果无人机目标距离明显上升,卡车距离下降。这说明航程放宽后,更多客户点可以转移给无人机服务,成本结构发生变化。第二组增加无人机数量,从1架增加到4架,总距离下降,但下降幅度递减,说明边际收益递减。第三组调整卡车容量,从50提高到100,对结果影响不大,说明在这个算例里容量不是瓶颈约束。
这类灵敏度测试在论文复现中很重要,它既能验证模型正确性,也能为后续算法改进提供方向。
6. 常见问题与排查技巧实录
6.1 常见报错与解决方案速查
复现过程中整理了一些高频问题,放在表里方便你自查:
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
| 报错X.shape不对 | 把单个个体传入_evaluate | 保持X为二维数组,循环处理每一行 |
| 目标值全是同一个数 | 解码器没有使用个体信息 | 检查随机密钥到路径的映射是否生效 |
| 约束违反值一直很大 | 初始化完全不满足可行域 | 用启发式方法生成初始可行解 |
| 跑几百代结果无变化 | 交叉概率过高导致结构破坏 | 适当降低变异概率,或提高eta |
| Pareto前沿点太少 | 种群太小或代数不够 | 增大pop_size和n_gen,至少300代 |
| 不同种子结果差异过大 | 问题多峰性强 | 固定seed做稳定性实验,取最优结果 |
第一次跑通后你会发现,大部分问题出在数据结构和索引上,真正复杂的算法逻辑反而好debug。所以写代码时多用字典保存中间结果,打印出解码后的路径,能省时间。
6.2 实操心得与避坑技巧
最后分享几条自己在项目中摸索出的经验。
第一,解码器一定要设计成“纯函数”。同样的数组输入必须得到同样的路径结构,不要在解码器内部修改外部状态。我刚开始在解码器里用了一个全局列表来记录已访问节点,结果多目标评估并行时相互污染,排查了一整天才发现问题。
第二,小算例验证时不要急着画Pareto前沿,先打印最优个体的卡车路径和无人机飞行列表,逐条核对路径上的约束。很多小错误在路径层面一眼就能看出来。
第三,做参数灵敏度分析时,一次只改一个参数,不要同时改两个。否则结果变化说不清楚是哪个参数引起的。
第四,pymoo的约束支配机制比你想象的更严格。如果你想让某个约束被“容忍”,可以给约束值乘以一个很小的系数,但不要用罚函数把约束揉进目标函数里,这样会破坏前沿的分布结构。
第五,如果打算复现论文的对比结果,最好保存每一代的种群目标值历史。我用save_history=True生成的记录,后来做多算法对比时直接就用上了,省了大量重跑的时间。
我在实际调试中比较大的体会是:模型公式里看似不起眼的一个下标,落到代码里可能就是一条完全不同的路径结构。复现论文的过程,本质上是在和论文作者“隔空对答案”。每次结果对上的时候,那种成就感比读十篇论文都强。希望你也能在自己的项目里跑出一个干净的Pareto前沿,有卡住的细节欢迎交流。