考虑无人机辅助的卡车配送路径优化_彭勇.pdf

第 28 卷 第 2 期2023 年 4 月工业工程与管理Industrial Engineering and ManagementVol.28 No.2Apr.2023考虑无人机辅助的卡车配送路径优化彭勇，张永辉，黎元钧（重庆交通大学 交通运输学院，重庆 400074）摘要:无人机具有飞行不受地形限制等优点而被应用于提升城市物流配送“最后一公里”的效率实践，但无人机配送受到最大飞行时间及最大载重等限制。结合卡车与无人机配送特点，研究了一类考虑无人机辅助的卡车配送路径优化问题。考虑无人机最大飞行时间、最大载重和飞行速度等因素，以配送完成时间最短为目标建立混合整数规划模型；采用结合自适应K-means聚类搜索的混合变邻域搜索算法进行求解。基于Solomon Benchmark C101、R101、RC101构建算例，分析结果表明：嵌入简单启发式算法的混合变邻域搜索算法能够较好求解所提出的路径优化问题，并能提高物流配送服务时效性；无人机飞行速度受载重影响越小，飞行持续时间越长，越有利于缩短配送时间。关键词:无人机辅助；路径优化；K-means聚类搜索算法；变邻域搜索算法中图分类号:U 492.3 文献标识码:AOptimization of Truck Distribution Route Considering Drone-assistedPENG Yong，ZHANG Yonghui，LI Yuanjun（School of Traffic&Transportation，Chongqing Jiaotong University，Chongqing 400074，China）Abstract:Drones have been used to improve the efficiency of urban logistics and distribution last mile because they are not restricted by terrain.However，they were limited by the maximum flight time and maximum load.Combining with the characteristics of truck and drone delivery，a class of truck distribution route optimization problem with drone assistance was studied.Considering the factors such as the maximum flight time，maximum load and flight speed of drone，a mixed integer programming model was established to minimize the delivery completion time.A hybrid variable neighborhood search algorithm combined with adaptive K-means clustering search was designed.A calculation example was constructed based on Solomon Benchmark C101，R101，and RC101.The analysis results show that the hybrid variable neighborhood search algorithm embedded with a simple heuristic algorithm can solve the proposed routing optimization problem and improve the timeliness of logistics and distribution services.The less the flight speed is affected by the load and the longer the flight duration，the more conducive to shortening the delivery time.Key words:drone-assisted;route optimization;K-means clustering search algorithm;variable neighbourhood search algorithm文章编号：1007-5429（2023）02-0031-09DOI：10.19495/ki.1007-5429.2023.02.004收稿日期：2021-06-09基金项目：重庆市社会科学规划项目（2019YBGL049）;教育部人文社会科学研究规划基金项目（17YJA630079）作者简介：彭勇（1973），重庆巴南人，教授，博士，主要研究方向为交通运输规划与管理。E-mail：。-31第 28 卷 彭勇，等：考虑无人机辅助的卡车配送路径优化1 引言 电子商务迅猛发展对城市物流配送“最后一公里”的效率提出了更高的要求。无人机与传统配送车辆相比，具有不受道路中断、操作灵活、飞行速度快等优势，但由于飞行里程、飞行时间和最大载重量受限制，无法适应不同类型的客户需求。故在传统配送模式中引入无人机配送有利于充分发挥两者的优势，提高配送效率。因此，研究无人机辅助的“最后一公里”配送问题更具有现实应用意义。无人机辅助车辆配送问题1可看作传统的旅行商问题（traveling salesman problem，TSP）和车辆路径问题（vehicle routing problem，VRP）的扩展，目前根据卡车及无人机数量主要分为无人机旅行商问题（traveling salesman problem with a drone，TSP-D）和多无人机多车辆路径优化问题（vehicle routing problem with drones，VRP-D）。文献 2-5 主要研究了卡车携带单架无人机配送路径问题，且无人机单次服务一个客户点。MURRAY和CHU2研究了飞行 助 手 旅 行 商 问 题（flying sidekick traveling salesman problem，FSTSP），即卡车携带无人机从配送中心出发服务客户点，在这些客户点处，无人机从卡车起飞为其附近的其他客户配送包裹，与此同时，卡车继续前往下一个客户点，无人机也返回至该客户点继续送货。FERRANDEZ 等3研究了卡车无人机协同的 TSP-D 问题，首先构造卡车的TSP路径，然后优化无人机在卡车路径中的发射和回收点。AGATZ等4在FSTSP基础上，考虑了无人机可以在同一个卡车服务节点被发送和回收的TSP-D问题。郭秀萍等5-6提出一种卡车-无人机联合配送模式，首先确定卡车配送点，其次优化经过所有配送点的卡车路径，最后解决无人机路径优化问题。目前有学者也考虑卡车携带多架无人机配送问题。LUO等7研究了无人机一次可服务多个客户点的TSP-D问题，即卡车与无人机协同进行配送，无人机的发射和回收都需要在卡车停在客户点才能够进行。胡觉亮等8考虑了由一辆卡车搭载多架无人机为客户进行运输的配送模式。上述研究2-4，7-8主要考虑卡车无人机同步运行的TSP-D问题，但实际配送中存在服务客户多、客户点需求量超过一辆卡车最大装载量的现象，因此有学者研究了VRP-D问题。WANG 等9研究了无人机辅助卡车配送的VRP-D 问题，即从单一卡车扩展到多辆卡车携带无人机完成配送服务，卡车可以在配送中心或客户点实现无人机的起飞和降落，且无人机保证由同一卡车完成收回操作。SCHERMER等 10 研究了无人机不在同一点起飞和降落，而且无人机起飞后电池确保瞬时达到满电状态的问题。目前在数学模型研究中，较多学者建立了以时间最短或以成本最低为目标的数学模型，如KITJACHAROENCHAI 等11忽略卡车最大容量和客户需求约束，建立了交货时间最少的目标函数。POIKONEN等12考虑无人机电量约束，建立了以时间最小化为主、成本节约为辅的目标函数。TAMKE 等13针对无人机飞行距离和速度约束，建立了以不同时间导向为目标的数学模型。然而，在配送实践中，配送方案的制定应综合考虑这些因素，且应考虑无人机的最大飞行时间及电池损耗受无人机载重的影响。同时，也应考虑无人机载重量限制对无人机配送的制约及道路限制等情况对卡车配送的制约。综上所述，本文充分考虑无人机最大飞行时间、载重量、载重因子等对方案的影响，提出了无人机辅助卡车配送模式，针对不同客户类型以配送服务完成时间最短为目标，建立了无人机辅助的卡车配送优化模型，基于问题的NP-hard特性，设计了一种混合变邻域搜索算法进行分析求解。2 问题描述与建模 2.1问题描述本文研究的问题具体可描述为：卡车在无人机辅助下从配送中心出发为多个客户点提供配送服务，完成所有客户点的配送后返回配送中心，并且卡车和无人机在配送过程中属于一一对应关系。由于不同因素影响，客户点被分为三类：一类是由于需求量超过了无人机的最大载重量，只能由卡车进行服务（TC）；一类是道路通行受限，只能由无人机提供服务（UC）；一类是客户点卡车和无人机都可以提供服务（FC），卡车最多搭载一架无人机完成配送，并且配送过程在确保两者时间约束一致性的-32第 2期工 业 工 程 与 管 理前提下，允许在卡车配送路径上的客户点处完成无人机发射和回收。图1是无人机辅助的卡车配送过程示意图。问题假设如下：（1）客户点的坐标和需求已知；（2）无人机一次只服务一个客户点；（3）无人机在某个客户点发射后，卡车必须在客户点原位置或者服务路径后续客户点（配送中心）对无人机进行回收；（4）无人机只能当卡车在客户点（配送中心）处于停留状态时才会发射；（5）无人机有最大飞行时间限制；（6）当卡车在某个客户点回收无人机时，存在无人机等待卡车或者卡车等待无人机或者同时到达客户点的情况；（7）无人机不考虑货物装卸时间和电量不足时的换电池时间；（8）每个客户点都可以通过卡车或无人机进行配送，但有且只能服务一次。2.2数学模型模型符号定义如表1所示。根据假设条件，建立如下无人机辅助的卡车配送路径优化问题的数学模型：Min Z=()lLtl0（1）s.t.jC-xl0j=iC+xli01，lL（2）xl0，0=0，lL（3）iC+xlij=0，lL，j C3：ij（4）iC+ijxlij=kC-jkxljk，lL，jC（5）lLiC+ijxlij+lLiC+ij kC-kjylijk=1，jC（6）eli-elj+1M(1-xlij)，lL，(i，j)A（7）eljMiC+ijxlij，lL jC-（8）jCjikC-kjylijk1，lL，i C+：ik（9）iC+ijjCjkylijk1，lL，k C-：ik（10）图1考虑无人机辅助的卡车配送问题示意图表1 模型符号定义名称参数连续变量决策变量符号dijd?ijvv?qiMC1C2C3CLDC+NC-iAiCtruck eliLRtiltilxlijylijkzlij定义客户点i到客户点j的车行距离客户点i到客户点j的无人机行驶距离卡车行驶速度无人机最大行驶速度无人机最大飞行时间无人机飞行速度载重影响因子客户点i的需求量一个无穷大数卡车无人机均可服务客户点集合，C1=1,2,c1仅卡车服务客户点集合，C2=c1+1,c1+2,c1+c2仅无人机服务客户点集合，C3=c1+c2+1,c1+c2+2,c1+c2+c3客户点集合，C=C1C2C3卡车（无人机）的集合，L=1,2,l 配送中心节点，D=0,0包含起点的所有客户点集合，C+=0 C网络中所有的节点，N=CD包含终点的所有客户点集合，C-=0C卡车在客户点i的服务时间卡车可行弧集合，A=(i,j):iC+,jC-,ij 无人机在客户点i的服务时间卡车的最大容量客户点i在卡车l路径中的序号无人机的发射时间无人机的回收时间卡车l到达客户点i的时间无人机l到达客户点i的时间卡车l经过可行弧(i,j)A，xlij=1；否则为0无人机l从节点 i 起飞服务节点j之后返回到节点k,ylijk=1;否则为0卡车l 的行驶路径中客户节点i 在客户节点j之前被访问,zlij=1;否则为0-33第 28 卷 彭勇，等：考虑无人机辅助的卡车配送路径优化iC+jikC-kjylijk=0，lL，j C2：ik（11）2ylijkhC+hixlhi+mCmkxlmk，lL，i，j C：ij，k C-：ik（12）jCqj()kC-jkxljk+iC+ijkC-jkylijkCtruck，lL（13）elj-eliMzlij，lL，(i，j)A（14）elj-eliM(zlij-1)+1，lL，(i，j)A（15）tlh+dhk/v+h+LjCkjmC-ylkjm+RnC+nkqCylnqktlk+M(1-xlhk)，lL，(h，k)A（16）tli+d?ij/(v?-qj)+L-M()1-kC-kjylijktlj，lL，(i，j)A（17）tlj+d?jk/v?+j+R-M()1-iC+ijylijktkl，lL，jC，kC-：kj（18）tli-M()1-jCjikC-kjylijktli，lL，iC+：ik（19）tli+M()1-jCjikC-kjylijktli，lL，iC+：ik（20）tlk-M()1-iC+ijjCkjylijktkl，lL，kC-：ik（21）tlk+M()1-iC+ijjCjkylijktkl，lL，kC-：ki（22）tlk-tli+M()1-jCjkylijk，lL，iC+：ij，kC-：ki（23）tmltlk-M()3-zlim-jCjmylijk-nCninkpC-pipkylmnp，lL，(i，k)A，mC：mi，gk（24）t0l，tl0=0，lL（25）eli，til，tli0，lL，iN（26）目标函数式（1）表示所有卡车完成配送服务总时间最短。式（2）表示卡车l至多提供一次配送服务，从配送中心离开一次。式（3）表示卡车l不允许从配送中心出发直接返回配送中心。式（4）表示卡车l不允许配送仅由无人机服务的客户。式（5）表示卡车配送过程中的出入流量平衡。式（6）表示每个客户点会由卡车或无人机服务。式（7）和式（8）表示卡车路径消除子回路。式（9）表示无人机l若从i点发射，有且只能服务一个客户点后回到卡车上。式（10）表示无人机l若从k点回收，回收前有且只能服务一个客户点。式（11）表示无人机l不允许配送仅由卡车服务的客户。式（12）表示无人机l在客户点i发射服务客户点j后在客户点k回收，发射点i和回收点k必须在卡车路径上。式（13）表示卡车的额定载重容量约束。式（14）和式（15）假设卡车l服务客户点i后再服务客户点j，则zlij=1，否则为0。式（16）表示卡车到达客户点k的时刻不晚于到达客户点h完成服务及无人机发射回收操作后的时刻。式（17）表示卡车l到达客户点i完成服务及发射的时刻不晚于无人机l到达客户点j的时刻。式（18）表示无人机l在客户点j服务完返回卡车完成回收操作的时刻不晚于无人机l到达客户点k的时刻。式（19）和式（22）表示卡车和无人机到达发射点和返回回收点的时刻一致性。式（23）表示无人机l的最大飞行时间约束。式（24）表示无人机l正在配送过程中，卡车上不会出现再次发射无人机的情形。式（25）初始化卡车和无人机l的出发时刻。式（26）表示变量的取值范围。3 变邻域搜索算法设计 变邻域搜索算法（variable neighborhood search，VNS）在1997年由MLADENOVI和 HANSEN14提出，目前被广泛应用于求解车辆路径问题的各种变形问题中15。本文针对无人机辅助的卡车配送路径优化问题，提出了一种混合变邻域搜索算法对其进行求解。算法流程如图2所示。其中，rand代表随机产生在（0，1）的数。首先，生成完整的卡车-无人机路径；其次，对生成的路径进行邻域操作和抖动操作；最后，采用局部搜索机制，进一步扩大搜索空间，对路径方案不断优化。-34第 2期工 业 工 程 与 管 理3.1构造初始解本文研究的是在卡车容量限制下有无人机辅助的卡车配送路径优化问题，为加快算法寻优速度，采取自适应K-means聚类算法聚类，降维生成较好的初始方案。在 K-means 聚类算法基础上增加一个最大化判断因子（maximum，MAX），主要通过在相同数据类型下改变不同的车辆容量和初始分区聚类数，多次对比得到不同目标值的设定，MAX的值不超过卡车容量。通过判断因子MAX的值改进K-means聚类算法，进一步得到比较合适的聚类数。自适应 K-means 聚类搜索算法的流程参见文献 5。根据聚类算法的结果，首先采用启发式方法得到卡车路径方案，其次选择距离无人机最近点和次最近点作为无人机可以发射和回收的节点，从而得到初始解。3.2抖动操作与局部搜索抖动操作包括插入操作和交换操作。插入操作是将某条子路径的一个客户点移至另一条子路径上，如图3所示。在选取客户点时，若被选择的客户点为无人机的起降点，则将无人机的起降点移至该子路径上与该客户点相邻的客户点上。交换操作是将两条不同的子路径上符合条件（交换后满足需求不大于卡车最大载重）的路径段进行交换，如图4所示。如遇到路径段上存在无人机起降点，则将起降点移至该路径上与该客户点相邻的客户点上。局部搜索是在子路径上的局部寻优过程，可描述为随机选取某条子路径中的两个客户点进行交换，如图 5所示。若交换客户点为无人机起降点，则无人机起降点随之移动。若发生无人机起降方向与卡车行驶方向不一致，则调整无人机起降顺序。若新方案目标值大于原方案目标值，则更新方案。图2混合变邻域搜索算法基本流程图3插入操作图4交换操作图5子路径局部搜索-35第 28 卷 彭勇，等：考虑无人机辅助的卡车配送路径优化4 算例分析 设卡车最大容量为500，车速为40，无人机最大飞行速度为30，发射与回收耗时均为2，最大持续飞行时间为50，飞行速度载重影响因子为0.06。上述参数都为相对数据，量纲为 1。自适应 K-means聚类算法中初始聚类数量为4，局部搜索次数为50，每次操作都以0.3的概率接受原来的劣解，迭代次数为500。测试环境为 MATLAB 2016a，CPLEX12.63，Intel（R）Core（TM）i5-1035G1 CPU 1.00GHz 1.19 GHz处理器、4G内存和Windows 10操作系统。4.1算法与优化求解器对比为验证不同规模算例采用CPLEX求解器与混合 变 邻 域 搜 索 算 法 的 有 效 性，选 取 Solomon Benchmark 标准实例中的 R101 数据随机生成了 7组不同规模的算例，每一组算例运行时间均不超过1 800 s，求解结果如表2所示。其中，C代表总客户点数量，D代表配送中心数量，-代表没有在给定的时间内找到最优目标值，Gap代表启发式算法与精确算法求解的偏差。由表2可知，混合变邻域搜索算法能得到比较接近精确解的目标值。但当算例规模适当增大时，CPLEX优化求解器已经不能在规定的时间内给出精确解，而混合变邻域搜索算法仍能在较短时间里给出全局近似最优解。4.2大规模案例算法有效性分析为进一步说明混合变邻域搜索算法的有效性，选取 Solomon Benchmark 标准实例的 C101、R101、RC101三类数据并适当调整，构建含有1个配送中心、60 个 FC、20 个 TC、20 个 UC 的大规模案例。图6为三类测试数据节点分布图。4.2.1自适应K-means聚类算法分析在三种数据类型下，给定初始分区聚类数量与车辆容量，分析其对自适应K-means聚类算法的影响。表3为不同参数组合下，算法独立重复运行10次得到的目标函数均值。由 表 3 可 知，在 R101、C101 和 RC101 三 组实验中，目标函数均值波动主要受卡车容量和初始分区聚类数量的影响。在 R101实验中，卡车容量和初始分区聚类数量对目标函数均值波动有较大影响，而在C101和RC101 实验中则影响较小。进一步分析，在RC101实验中，当卡车容量、初始分区聚类数量较大时，目标函数值更优。4.2.2混合变邻域搜索算法分析为分析抖动操作对混合变邻域搜索算法求解得到的目标函数最优值的影响，针对三类数据类型，分别移除不同操作模块，取10次运算结果如表4所示，其中最后一行加粗数字为均值。由表4可知，对于不同数据类型，抖动操作算法表2CPLEX与混合变邻域算法结果对比算例规模C+D5+16+18+118+128+1CPLEX优化求解器目标值/min40.6244.9054.55-运行时间/s7.2247.96262.441 800.001 800.00混合变邻域搜索算法目标值/min46.0051.0062.00169.00281.00运行时间/s9.3511.9512.9632.1667.52Gap/%13.2313.5913.66-图6三类数据类型分布-36第 2期工 业 工 程 与 管 理优化性能影响很大，抖动操作和抖动插入操作对算法求解稳定性产生较大影响。图7为不同数据类型下最优目标的收敛情况。由图7可知，C101、R101、RC101三种数据类型分别在迭代150次、250次、200次左右收敛。结果表明，不同数据类型对算法迭代收敛次数有一定的影响，但算法均在较少迭代次数下收敛。图8为三种类型测试数据最优配送路径方案。由图8可知，不同类型的数据都有7条最优配送路径。由于其分布的随机性导致配送路径呈现出不同特点，根据最优路径图可以看出，结合自适应K-means聚类搜索的混合邻域搜索算法有较好的不同表3不同车辆容量与初始分区聚类数量的目标函数均值车辆容量500600700800900数据类型C101R101RC101C101R101RC101C101R101RC101C101R101RC101C101R101RC101初始分区聚类数量3761.4816.8872.6862.6865.8842.6805.8842.2791.0809.8823.2863.4816.4815.4825.64737.2772.6849.4839.4782.4889.0738.4809.4839.2774.6834.8871.8738.2776.5845.45796.8785.2851.6743.8792.6801.5754.6769.0811.4786.2736.4821.5773.4816.4801.86803.0784.6839.0796.6847.6759.6748.0792.4776.8752.4816.0763.4720.8830.2755.47744.6812.8750.6781.6805.4769.2791.2822.8741.6785.2806.5744.5747.2813.6734.6表4三种数据类型分布下不同模块的目标函数最优值正常C101758744773755741752769786763772761.3R101885863921765751866855937906876862.5RC101911946831887845932849826855834871.6缺抖动操作C1019831 0069989831 0231 1011 0479769801 0061 010.3R1011 0471 0651 0391 1011 1231 0141 0449981 0361 0471 051.4RC1011 0561 1241 1561 0019788991 0251 1091 0421 0631 045.3缺抖动插入操作C101945865978932888904931962911863917.9R1011 1121 0231 0561 0891 0011 0229989871 0561 0851 042.9RC1019991 0251 0261 1111 0789869781 0231 0251 0011 025.2缺抖动交叉操作C101920835906967821852832799815917866.4R1011 0321 0351 0211 0029421 0561 1421 0891 0141 0091 034.2RC1011 0251 0021 0091 1021 0361 0691 0451 0241 0521 0531 041.7图7三类数据分布下算法求解迭代过程-37第 28 卷 彭勇，等：考虑无人机辅助的卡车配送路径优化数据类型的适用性。4.3无人机性能参数对方案的影响分析性能参数设置：飞行速度载重影响因子分别取0.06、0.08、0.10、0.12、0.14，最大飞行时间为20、30、40、50、60，并将参数两两组合，进行25组实验，每种实验算法运行10次并取均值，结果如图9所示。由图9可知，无人机飞行速度载重影响因子和无人机最大飞行时间均对目标函数值有影响。纵向分析，其目标函数值随无人机飞行速度载重影响因子的增加而不断增大，随无人机最大飞行时间的增加而不断减小。横向对比，无人机飞行速度载重影响因子对目标函数值影响程度低于无人机最大飞行时间。5 结论 本文研究了考虑无人机辅助的卡车配送路径优化问题，以配送完成总时间最短为目标建立数学模型，并设计了混合变邻域搜索算法。算法流程为：首先，采取自适应K-means聚类算法聚类；其次，图8三类数据分布下算法最优路径图图9载重影响因子和最大飞行时间变化对优化目标的影响-38第 2期工 业 工 程 与 管 理运用迭代最近点算法给出初始路径；最后，采用变邻域搜索算法优化初始路径，最终获得配送总时间最短的路径方案。算法有效性分析证明本文所设计算法能有效求解所提问题，案例分析还发现无人机性能对最优方案产生影响。考虑无人机辅助的卡车配送新模式在面临城市配送实际问题时可给决策者提供可行的解决方法，未来可从构建更有效的启发式算法、增加与benchmark 方法或结果的对比实验等角度进行研究，进一步提高求解的精度。参考文献：1 彭勇，黎元钧.考虑疫情影响的卡车无人机协同配送路径优化 J.中国公路学报，2020，33（11）：73-82.2 MURRAY C C，CHU A G.The flying sidekick traveling salesman problem：Optimization of drone-assisted parcel deliveryJ.Transportation Research Part C：Emerging Technologies，2015，54：86-109.3 FERRANDEZ S M，HARBISON T，WEBER T，et al.Optimization of a truck-drone in tandem delivery network using k-means and genetic algorithmJ.Journal of Industrial Engineering and Management，2016，9（2）：374-388.4 AGATZ N，BOUMAN P，SCHMIDT M.Optimization approaches for the traveling salesman problem with drone J.Transportation Science，2018，52（4）：965-981.5 郭秀萍，胡运霞.卡车与无人机联合配送模式下物流调度的优化研究 J.工业工程与管理，2021，26（1）：1-8.6 杨双鹏，郭秀萍，高娇娇.无接触式“卡车+无人机”联合配送问题研究 J.工业工程与管理，2022，27（1）：184-194.7 LUO Z，LIU Z，SHI J.A two-echelon cooperated routing problem for a ground vehicle and its carried unmanned aerial vehicle J.Sensors，2017，17（5）：1144.8 胡觉亮，杨航，张田会，等.无人机与卡车协同配送优化研究J.浙江理工大学学报（社会科学版），2020，44（5）：489-497.9 WANG X，POIKONEN S，GOLDEN B.The vehicle routing problem with drones：Several worst-case results J.Optimization Letters，2017，11（4）：679-697.10 SCHERMER D，MOEINI M，WENDT O.A hybrid VNS/Tabu search algorithm for solving the vehicle routing problem with drones and en route operations J.Computers&Operations Research，2019，109：134-158.11 KITJACHAROENCHAI P，VENTRESCA M，MOSHREF-JAVADI M，et al.Multiple traveling salesman problem with drones：Mathematical model and heuristic approachJ.Computers&Industrial Engineering，2019，129：14-30.12 POIKONEN S，WANG X，GOLDEN B.The vehicle routing problem with drones：Extended models and connections J.Networks，2017，70（1）：34-43.13 TAMKE F，BUSCHER U.A branch-and-cut algorithm for the vehicle routing problem with drones J.Transportation Research Part B：Methodological，2021，144：174-203.14 MLADENOVI N，HANSEN P.Variable neighborhood search J.Computers&Operations Research，1997，24（11）：1097-1100.15 REZGUI D，SIALA J C，AGGOUNE-MTALAA W，et al.Application of a variable neighborhood search algorithm to a fleet size and mix vehicle routing problem with electric modular vehicles J.Computers&Industrial Engineering，2019，130：537-550.-39