负载依赖型电动货运三轮车最后一千米同时取送货路径优化研究

Electric cargo bicycle’s last-mile simultaneous pickup and delivery routing problem with load-dependent travel time

摘要

电动货运三轮车在最后一千米配送活动中起着十分关键的作用.然而，在研究最后一千米配送相关问题时，极少有研究讨论电动货运三轮车行驶速度与其负载间的关系.本文研究了具有负载依赖行驶时间的电动货运三轮车最后一千米同时取送货问题（PDPLDTT），其中行驶速度取决于道路坡度和电动货运三轮车的负载.为了解决PDPLDTT，本研究首先构建了一个混合整数线性规划（MILP）模型，并使用了商业求解器CPLEX求解该模型的小规模案例；其次，提出了一种改进的混合蚁群优化（HACO）算法对该问题的大规模案例进行求解分析；最后，对比分析了HACO算法与对照算法的求解结果.结果表明：本文提出的MILP模型与HACO算法可以有效解决具有负载依赖行驶时间的电动货运三轮车最后一千米同时取送货问题，能为相关企业解决最后一千米取送货问题提供合理的决策建议.

Abstract

Electric cargo bicycles play a crucial role in the last mile delivery activity. However， few studies have discussed the relationship between the speed of electric cargo bicycles and their load when studying issues related to last-mile delivery. This paper investigates the electric cargo bicycles last-mile simultaneous pickup and delivery problem of with load-dependent travel time （PDPLDTT）， where the travel speed depends on the road slope and the load of the electric cargo bicycle. In order to solve the PDPLDTT， this study first developed a mixed-integer linear programming （MILP） model is proposed， and small-scale instances of the MILP model are solved by the commercial solver CPLEX； Then， an improved hybrid ant colony optimization （HACO） algorithm was proposed to solve large-scale instances of this problem； Finally， the solution results of HACO algorithm and other benchmark algorithms were compared and analyzed. The results show that the proposed MILP model and HACO algorithm in this paper can effectively solve the electric cargo bicycles last-mile simultaneous pickup and delivery problem with load-dependent travel time， and can provide reasonable decision-making suggestions for relevant enterprises to solve the last-mile pick and delivery problem.

Graphical abstract

关键词

最后一千米配送 / 负载依赖行驶时间 / 混合整数线性规划 / 混合蚁群优化算法

Key words

last-mile delivery / load-dependent travel time / mixed-integer linear programming / hybrid ant colony optimization algorithm

近年来，随着电子商务的蓬勃发展，网上购物已逐渐成为人们日常消费的重要途径.据国家统计局数据显示，2022年我国网上零售额达到13.79万亿元.庞大的网络购物市场使得网上订单包裹数量逐年递增，据《2021年中国快递发展指数报告》显示，2021年快递配送量达到1 083亿件，首次突破千亿件，日均快件处理量近3亿件.网上购物的迅猛发展给我国物流企业带来巨大的最后一千米配送压力^［1］.“最后一千米”虽然在整个物流环节中的距离最短，但是却消耗了总物流成本的30%左右^［2］.如何提升配送效率、降低配送成本、提高服务质量，是所有物流企业最后一千米配送时面临的重大难题.

传统的最后一千米配送通常由卡车或货车负责完成.然而，最后一千米配送大多发生在城市或者乡间小路，拥堵的交通将大幅度降低卡车的配送效率.同时，货运将增加交通拥堵程度，并导致大量尾气排放.因此，为了缓解交通拥堵、减少尾气排放并提高物流效率，包括DHL、顺丰快递和EMS在内的海内外快递服务商都使用了一种非常经济的模式，即使用电动货运三轮车进行最后一千米配送.在我国，电动货运三轮车是有效的最后一千米运输工具，但也有一些限制，例如不可载人、不可超载、不可超速.图1展示了中国邮政及顺丰快递在进行配送时使用的电动货运三轮车.电动货运三轮车货箱的装载能力可达200 kg以上^［3］，足以应对最后一千米配送所需的装载量.此外，由于电动货运三轮车不属于机动车，驾驶员可以经常选择非机动车道行驶以避免道路拥堵.

尽管交付方式将运输工具从卡车改变为电动货运三轮车，但物流公司仍需要为驾驶员设计从仓库访问客户的计划.这个问题属于车辆路径问题（vehicle routing problem，VRP）的变种问题^［4］.配送员驾驶电动货运三轮车给一组分布在指定区域内的客户进行包裹配送服务，由于电动三轮车的电机功率是固定的，其最大行驶速度受到三轮车负载、道路坡度和其他因素的影响.因此，使用电动货运三轮车配送的路径规划问题实际上属于具有负载依赖行驶时间的VRP问题（VRP with load-dependent travel time，VRPLDTT）^［5］.此外，最后一千米物流中的收取货情景也很重要.一些客户可能会在收到新商品时退回他们之前购买的一些商品；客户也可能因为不满意商品质量而拒收商品.因此，在实际的配送过程中，配送员不仅会送货，同时可能会收取货物.所以在VRPLDTT问题里考虑取货约束，可以使模型更加符合现实世界的问题情景.这也使得电动三轮车的实际负载更加具有随机性，从而导致行驶速度的变化更复杂.

在有关负载依赖的配送优化研究中，考虑负载依赖行驶时间的研究并不多，研究主要集中于负载依赖成本方面，其中最为经典的问题是由Bektacs等^［6］在2011年提出的污染路径问题（pollution routing problem， PRP）.该研究建立的混合整数线性规划模型（mixed-integer linear programming，MILP），通过使用速度水平来线性化处理速度变量，给很多研究带来了线性化处理变量的新思路.经过多年研究，目前已经开发出很多PRP及其衍生问题的求解算法.Demir等^［7］在2012年提出了ALNS算法来求解PRP问题，该算法可以在短时间内有效地解决PRP问题的大规模实例.Kramer等^［8］提出了基于局部搜索算法和数学求解器的数学启发式算法，实验结果证明了所提出方法的有效性.尽管这些研究在求解PRP问题方面具有一定的优点，但建模方法和求解算法难以有效解决VRPLDTT问题.

基于已有研究，本文研究了具有负载依赖行驶时间的电动货运三轮车最后一千米同时取送货问题（pickup and delivery problem with load-dependent travel time，PDPLDTT），其中行驶速度取决于道路坡度和电动货运三轮车的负载.为了求解PDPLDTT问题，本文构建了一个混合整数线性规划模型，并设计了一种改进的混合蚁群优化（hybrid ant colony optimization，HACO）算法以解决此问题的大规模实例.相比现有文献，本研究可能的边际贡献主要体现在以下方面：1） 首次研究了具有负载相关行驶时间的电动货运三轮车最后一千米同时取送货问题；2） 对此问题建立了数学规划模型，并使用了CPLEX验证模型的正确性；3） 提出了HACO算法来解决此问题的大规模案例，并使用了大量算例来验证HACO算法的有效性.

1 问题描述与数学模型

1.1 问题描述

本文研究的PDPLDTT问题里包含一个配送中心，配送中心里有若干待配送的包裹，每个包裹有一定的重量且重量已知，配送工作由若干同质的电动货运三轮车完成.具体来说，PDPLDTT问题可以描述为：假设有一个由一系列点N=\{\mathrm{0,1},\dots ,n,n+\mathrm{1}\}和相应的弧A组成的有向图G=\left(N,A\right)，其中\mathrm{0}和n+\mathrm{1}分别表示配送起点和终点，并指向同一个配送中心.N_{\mathrm{0}}=\{\mathrm{1,2},\dots ,n\}代表着一系列有配送需求的客户，其中，客户i所需要配送包裹的重量为q_i^d，所需要提取包裹的重量为q_i^p,所需要的服务时间为s_i.同时，客户i有服务时间窗\left[a_i,b_i\right]要求，即a_i是最早可以开始接受服务的时间，b_i是最晚可以开始接受服务的时间.对于点i\in \{\mathrm{0},n+\mathrm{1}\}，本文假设q_i^d=q_i^p=s_i=a_i=\mathrm{0}并且b_i=T，其中T表示电动货运三轮车最晚回到配送中心的时间.d_{ij}表示点i与点j之间的距离，t_{ij}表示点i与点j之间的行驶时间.Q代表电动货运三轮车的最大装载能力.

PDPLDTT问题是在一定的约束和目标下为配送人员设计合理的配送路线.通常来说，这类问题的目标是最大限度地减少总行驶时间或行驶距离.在本文中，电动货运三轮车在不同弧中的行驶速度不是一个固定常数，而是由负载决定的参数，进而电动货运三轮车在每段弧中的行驶时间也是由负载决定的参数.最小化总行程时间将比最小化行程距离更重要.因此，本文所研究的PDPLDTT问题考虑的目标函数是最小化总行驶时间.此外，本文做出如下假设：1） 在任意两个点（包括配送中心和所有客户）间，电动自行车将保持恒定的负载依赖速度行驶；2） 电动货运三轮车的最大行驶速度为25 km·h^-1；3） 电动货运三轮车将以最大额定功率进行行驶.

1.2 负载依赖的行驶时间

在PDPLDTT问题里，电动货运三轮车被用来对客户进行最后一千米配送服务，而电动货运三轮车主要是靠电机驱动车轮行驶.根据物理学的基本定理，电动三轮车在行驶过程中会克服四个方面的阻力^［9］，具体如下.

空气阻力F_D，其数学公式可以写为：

F_D=\frac{\rho C_{d}A}{\mathrm{2}}{v}^{\mathrm{2}},\left(\mathrm{1}\right)

其中：C_d表示空气阻力系数，\rho是空气密度，A表示电动自行车前部区域面积，v代表电动自行车的行驶速度.

滚动摩擦力F_R，其数学公式可以写为：

F_R=C_r\times m\times g\times \mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(h\right)\right),\left(\mathrm{2}\right)

其中：C_r表示滚动摩擦力系数，m指电动自行车自身及其负载的总质量，g是重力加速度，h表示道路的坡度.

重力阻力F_G，其数学公式可以写为：

F_G=m\times g\times \mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(h\right)\right).\left(\mathrm{3}\right)

能量损失F_F，能量损失是在链条和车轮轴承中能量传动时产生的.本文根据文献数据，设定能量损失为系统总功率的5%.

因此，根据物理学经典公式P=Fv可以得到：

P=\left(F_D+F_R+F_G+F_F\right)v,\left(\mathrm{4}\right)

其中，P表示电动货运三轮车的额定功率.

将四种阻力的公式代入上式，可以得到：

\begin{array}{l}P=\left(F_D+F_R+F_G\right)v+\mathrm{0.05}P=\frac{\left(F_D+F_R+F_G\right)v}{\mathrm{0.95}}=\\ \frac{\left(\begin{array}{c}\frac{\rho C_{d}A}{\mathrm{2}}{v}^{\mathrm{2}}+C_{r}mg\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(h\right)\right)+mg\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(h\right)\right)\end{array}\right)v}{\mathrm{0.95}}.(\mathrm{5})\end{array}

从公式（5）不难发现，针对一个特定的案例，额定功率P、空气密度\rho、空气阻力系数C_d、电动自行车前部区域面积A、滚动摩擦力系数C_r、重力加速度g以及坡度h均为特定的参数，仅有电动自行车自身及其负载的总质量m和电动自行车行驶速度v为变化的变量.此外，电动货运三轮车总质量由三部分组成，分别是电动货运三轮车空载质量、驾驶员重量以及货物载重.而电动货运三轮车的空载质量以及驾驶员质量可视为定值，仅有货物载重是随着配送过程进行变化的变量.根据此分析，可以得到如下的命题.

命题 1 对于一段指定的路程， 恒定功率的电动货运三轮车的最大行驶速度v仅仅只与电动货运三轮车的负载有关.

证明 公式（5）中，仅有m、h以及v三个变量.而对于指定路程而言，坡度h也不会变化.因此，速度v将仅仅只与电动货运三轮车的负载相关.得证.

根据命题1可知，电动货运三轮车在送货或者取货后由于货物载重的变化，其行驶速度也会发生变化，即负载依赖的行驶速度.因此，负载依赖的行驶时间模型的核心是公式（5），本文采纳了文献［5］中的参数设置进行了计算实验，其具体参数如表1所示.

将表1的参数设置代入式（5），可以得到：

\begin{array}{l}\mathrm{0.608}{v}^{\mathrm{3}}+[\mathrm{0.103}\mathrm{c}\mathrm{o}\mathrm{s}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(h\right)\right)+\\ \mathrm{10.326}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{a}\mathrm{n}\left(h\right)\right)]mv=\mathrm{350},\left(\mathrm{6}\right)\end{array}

从式（6）不难发现，在一段确定坡度的道路上，车辆行驶速度将随着负载的变化而变化.本文考虑了坡度分别为0、0.02、0.05以及0.1四种情景下，电动货运三轮车行驶速度在从货箱空载（140\mathrm{k}\mathrm{g}）到货箱满载（290\mathrm{k}\mathrm{g}）下的变化情况，如图2所示.可以看到，电动货运三轮车的行驶速度v随着货载总质量m的增加而降低.针对坡度h=\mathrm{0}情景，也满足行驶速度v随着货载总质量m的增加而降低这个结论.

由于电动货运三轮车的控制程序设置了速度上限25 km·h^-1，所以当行驶速度大于25 km·h^-1的时候，只能给车轮输出25的速度控制指令.

1.3 数学规划模型

本文通过构建MILP模型来解决所研究的PDPLDTT问题.该MILP模型涉及一个重要的二进制决策变量x_{ij}，其含义如下：

\begin{array}{c}\begin{array}{l}{x}_{ij}=\left\{\begin{array}{l}\mathrm{1},\mathrm{如果}\mathrm{配送}\mathrm{员先}\mathrm{服务}\mathrm{客户}i,\mathrm{再服}\mathrm{务客}\mathrm{户}j;\\ \mathrm{0},\mathrm{否则}.\end{array}\right.\\ \left(\mathrm{7}\right)\end{array}\end{array}

该模型涉及的另外两个重要变量分别是质量m和行驶速度v.由公式（6）可知，m与v之间存在着函数关系，本文假设v=f\left(m\right).由于行驶时间为距离除以速度，即t=\frac{d}{v}.综合这两个公式，可以得到：t=\frac{d}{f\left(m\right)}.公式中d为确定参数，而分母中的m为变量.因此，行驶时间t会呈现非线性变化，从而使得数学模型不是线性模型.为了线性化处理这一变量，本文借鉴了文献［6-7］中处理速度变量的方式，引入了一个新的参数与二进制决策变量，分别为质量等级参数l\in L=\{\mathrm{1,2},\dots ,l,\dots \}和二进制决策变量z_{ij}^l.

本文首先将电动货运三轮车的总负载能力Q划分为l\in L=\{\mathrm{1,2},\dots ,l,\dots \}个负载级别，即负载区间\left[\mathrm{0},Q\right]被分割成|L|个负载区间，分别为：\left[\mathrm{0},\frac{\mathrm{1}}{|L|}Q\right]，\left[\frac{\mathrm{1}}{|L|}Q,\frac{\mathrm{2}}{|L|}Q\right]，\dots，\left[\frac{|L|-\mathrm{1}}{|L|}Q,Q\right].为了简化处理，引入两个新的参数p以及r，并满足以下关系：

\begin{array}{r}{p}^l=\frac{l-\mathrm{1}}{\left|L\right|}Q,l\in L=\mathrm{1,2},\dots ,l,\dots ,\\ r^l=\frac{l}{\left|L\right|}Q,l\in L=\mathrm{1,2},\dots ,l,\dots .\end{array}\left(\mathrm{8}\right)

所以第l个负载等级对应的负载区间为\left[p^l,r^l\right],l\in L.本文使用负载区间中间值\frac{p^l+r^l}{\mathrm{2}}来代替落在这一负载区间的实际负载值，从而求得这一区间对应的行驶速度与行驶时间.即针对弧\left(i,j\right)，负载水平l对应的行驶速度v_{ij}^l和行驶时间t_{ij}^l可以根据公式（6）、弧\left(i,j\right)对应的坡度以及包含负载值\frac{p^l+r^l}{\mathrm{2}}在内的总质量求得.那么针对特定案例，v_{ij}^l和t_{ij}^l成为负载等级l对应的特定参数.故本文只需决定在弧\left(i,j\right)上的负载水平值.因此，本文引入了一个新的决策变量z_{ij}^l，其含义如下：

\begin{array}{c}\begin{array}{l}{z}_{ij}^l=\left\{\begin{array}{l}\mathrm{1},\mathrm{配送}\mathrm{员先}\mathrm{服务}i,\mathrm{再服}\mathrm{务}j,\mathrm{且负}\mathrm{载水}\mathrm{平为}l;\\ \mathrm{0},\mathrm{否则}.\end{array}\right.\\ \left(\mathrm{9}\right)\end{array}\end{array}

MILP模型中还涉及到一些其他的变量和参数，具体如表2所示.基于这些参数与变量，其具体的MILP模型如下：

Z=\mathrm{m}\mathrm{i}\mathrm{n}\sum _{\left(i,j\right)\in A}\sum _{l\in L}{t}_{ij}^l{z}_{ij}^l,\left(\mathrm{10}\right)

s.t.

\sum _{j\in N\setminus \{\mathrm{0}\}}{x}_{ij}=\mathrm{1},\forall i\in N_{\mathrm{0}},\left(\mathrm{11}\right)

\sum _{i\in N\setminus \{n+\mathrm{1}\}}{x}_{ij}=\mathrm{1},\forall j\in N_{\mathrm{0}},(\mathrm{12})

\sum _{j\in N}{f}_{ji}^d-\sum _{j\in N}{f}_{ij}^d=q_i^d,\forall i\in N_{\mathrm{0}},\left(\mathrm{13}\right)

\sum _{j\in N}{f}_{ij}^p-\sum _{j\in N}{f}_{ji}^p=q_i^p,\forall i\in N_{\mathrm{0},}\left(\mathrm{14}\right)

f_{ij}^d+f_{ij}^p\le Q{x}_{ij},\forall \left(i,j\right)\in A,\left(\mathrm{15}\right)

\sum _{l\in L}{p}^l{z}_{ij}^l\le f_{ij}^d+f_{ij}^p\le \sum _{l\in L}{r}^l{z}_{ij}^l,\forall \left(i,j\right)\in A,\left(\mathrm{16}\right)

\begin{array}{l}{y}_i-y_j+s_i+\sum _{l\in L}{t}_{ij}^l{z}_{ij}^l\le M_{ij}\left(\mathrm{1}-x_{ij}\right),\forall \left(i,j\right)\in A,\\ \left(\mathrm{17}\right)\end{array}

a_i\le y_i\le b_i,\forall i\in N_{\mathrm{0}},\left(\mathrm{18}\right)

\sum _{l\in L}{z}_{ij}^l=x_{ij},\forall \left(i,j\right)\in A,\left(\mathrm{19}\right)

x_{ij}\in \left\{\mathrm{0,1}\right\},\forall \left(i,j\right)\in A,\left(\mathrm{20}\right)

z_{ij}^l\in \left\{\mathrm{0,1}\right\},\forall \left(i,j\right)\in A,l\in L,\left(\mathrm{21}\right)

f_{ij}^d\ge \mathrm{0},f_{ij}^p\ge \mathrm{0},\forall \left(i,j\right)\in A,\left(\mathrm{22}\right)

y_i\ge \mathrm{0},\forall i\in N.\left(\mathrm{23}\right)

公式（10）为目标函数，意为最小化总行驶时间.公式（11）和（12）确保每个客户都被访问一次.公式（13）和（14）为负载流平衡约束，确保每个客户的配送与寄送需求被满足.公式（15）为负载约束，确保电动货运三轮车不会超载.公式（16）用来确定电动货运三轮车所处的负载区间，并决定决策变量z_{ij}^l的值.公式（17）和（18）组成了硬时间窗约束，其中M_{ij}为时间窗约束中的一个下界参数{}^{[\mathrm{6}]}，且M_{ij}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\mathrm{0},b_i-a_j+s_i+t_{ij}^{\mathrm{1}}\right\}.公式（19）为二进制决策变量间的关系，确保弧\left(i,j\right)上只有一个负载等级.公式（20）和（21）为二进制约束.公式（22）和（23）为非负约束.

2 混合蚁群优化算法

本文所研究的问题是VRP问题的衍生问题.VRP问题是典型的NP-hard问题，其求解难度会随着问题规模的增加呈指数级增长^［10］.由于负载依赖旅行时间约束的存在，本文所研究的PDPLDTT问题比传统VRP问题更为复杂，但同样是NP-hard问题.针对NP-hard问题，使用精确方法很难在短时间内求解此类问题的大规模案例.因此，学者们通常使用启发式算法来进行求解，如ALNS^［11］、蚁群优化算法（ant colony optimization，ACO）^［12］等.由于ACO算法在解决组合优化相关问题时有较好的求解效果^［13］，本文选择设计改进的ACO算法来求解PDPLDTT问题.ACO算法是意大利学者Dorigo在1996年受蚂蚁觅食的启发提出的一种全局优化算法^［14］.自提出后，ACO算法被广泛用于各类组合优化问题，如旅行商问题（traveling salesman problem，TSP）^［15］、VRP问题^［16］等.学者们在使用ACO算法时，通常会做出一定的改进以增强算法的求解效果.本文设计了一种混合蚁群优化（HACO）算法来求解PDPLDTT问题，具体改进方法是使用多种局部搜索算子与ACO算法进行结合，以提高解的质量.HACO算法的流程图如图3所示，算法具体过程在下文描述.

2.1 解的构造

在HACO算法解的构造方法中，解决方案是由蚂蚁构建的，并由一组路线组成.本研究使用实数编码方式对解进行构造.值得注意的是，解决方案是由一只蚂蚁而非多只蚂蚁构建.也就是说，蚂蚁构建了一条路线后，它会返回原始仓库并重新规划另一条路线，直到构建出完整的解决方案.在HACO算法中，蚂蚁访问完一个节点后，会根据概率规则选择下一个节点，具体的概率公式为：

\begin{array}{c}{P}_{ij}^k=\left\{\begin{array}{cc}\frac{{\left({\tau }_{ij}\right)}^{\alpha }{\left({\eta }_{ij}\right)}^{\beta }}{\sum _{l\in C_i^k}{\left({\tau }_{il}\right)}^{\alpha }{\left({\eta }_{il}\right)}^{\beta }},& \hspace{0.25em}\mathrm{如果}j\in C_i^k;\\ \mathrm{0},& \mathrm{否则},\end{array}\right.\left(\mathrm{24}\right)\end{array}

其中，{\tau }_{ij}表示节点i和节点j之间的信息素浓度，不同的蚂蚁通过更新信息素的浓度来获得周围环境的信息，并与其他蚂蚁进行交流；C_i^k表示蚂蚁k在访问节点i后可以访问的可行候选者的集合.\alpha和\beta是两个重要的权重参数.{\eta }_{ij}是可见性参数，其求解公式为：

{\eta }_{ij}=\frac{\mathrm{1}}{d_{ij}}.\left(\mathrm{25}\right)

HACO算法构造解的具体过程如下.

1） 蚂蚁从配送中心出发，准备服务客户集C中的客户（C表示该蚂蚁还未服务客户的集合，其初始状态为包含所有客户的集合）.

2） 该蚂蚁随机服务C中的客户i，并更新集合C.如果集合C不是空集，则根据约束条件计算在访问i后该蚂蚁可以访问的有效候选客户集合C\mathrm{\text{'}}.集合C\mathrm{\text{'}}是集合C的子集，即C\mathrm{\text{'}}\subseteq C.

3） 判断集合C\mathrm{\text{'}}是否为空集.如果集合C\mathrm{\text{'}}为空集，那么该蚂蚁将返回配送中心完成一条路线的规划，并从配送中心开始进行新的路线规划.

4） 如果集合C\mathrm{\text{'}}不是空集，那么蚂蚁根据公式（24）选择服务集合C\mathrm{\text{'}}中的客户j\in C\mathrm{\text{'}}，并更修集合C和C\mathrm{\text{'}}.

5） 如果集合C是空集，那么蚂蚁将回到配送中心并停止搜索.此时，蚂蚁已经构造出了该问题的可行解.

2.2 解的评估

在构造完一个可行解后，如何评估此可行解的质量是一个重要问题.本文所研究问题的目标是最小化行驶时间.因此，评估解的关键是求出任意两点间的行驶时间.那么根据上一节对负载依赖行驶时间的描述，在已知各个参数的数值以及两点间坡度大小后，可以算出不同负载级别下的行驶时间.即，先根据公式（6）计算负载依赖的行驶速度v_{ij}^l，然后得到行驶时间t_{ij}^l=\frac{d_{ij}}{v_{ij}^l}.

2.3 局部搜索算子

原始ACO算法在求解此类问题时易陷入局部最优，因此本文选择使用有较强搜索能力的局部搜索算子来提高HACO算法的全局优化能力^［17］.此种改进的有效性已经在很多文献中得到证明^［18-19］.本文设置了一个局部搜索判断阈值参数P_l，如果0至1之间的随机数大于阈值参数P_l，那么该蚂蚁将使用下面的局部搜索操作来提高解的质量.如果解的质量得到提高，那么接受局部搜索产生的新解；否则不接受该解.

本文使用了三种经典的局部搜索算子， 分别如下.

1） 2-opt操作.在一条完整路径中随机地选择两位客户，并反转两位客户间的路径（该路径包含这两位客户）.

2） 2-opt^*操作.在两条不同的路径中选择两个部分路径，并交换这两个部分路径的位置.

3） relocate操作.在一条路径中随机地选择一位客户，并将这位客户插入解中最合适的位置.

2.4 信息素更新和算法停止准则

HACO算法在完成局部搜索优化操作后，将进行信息素更新并判断算法是否达到终止条件.信息素更新包含两部分：全局衰减和局部增加.在完成局部搜索操作后，信息素将进行全局衰减操作，公式为：

{\tau }_{ij}=\left(\mathrm{1}-\rho \right){\tau }_{ij},\hspace{0.25em}\forall \left(i,j\right),\left(\mathrm{26}\right)

其中，\rho \in \left[\mathrm{0,1}\right]是信息素的调节参数.全局衰减后，信息素将进行局部增加操作，公式为：

{\tau }_{ij}={\tau }_{ij}+{\displaystyle \sum _{k=\mathrm{1}}^n}\Delta {\tau }_{ij},\left(\mathrm{27}\right)

其中，

\Delta {\tau }_{ij}=\left\{\begin{array}{ll}\frac{Q}{L_{An{t}_k}},& \mathrm{如果}\left(i,j\right)\mathrm{是}An{t}_k\mathrm{所构}\mathrm{造};\\ \mathrm{0},& \mathrm{否则}.\end{array}\right.

L_{An{t}_k}是An{t}_k所构造解的目标函数值，Q是一个常数.在完成信息素更新算子后，HACO算法将进行判断操作决定算法是否终止.HACO算法如果达到最大迭代次数\mathrm{M}\mathrm{a}\mathrm{x}\mathrm{I}\mathrm{t}，或者该算法持续\mathrm{M}\mathrm{a}\mathrm{x}\mathrm{C}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}代都没有搜索到更优的解决方案，算法将终止迭代.

3 计算实验

为测试所提出HACO算法和MILP模型的有效性，基于文献［5］中的实例生成数值实验所需算例.本文所提出的HACO算法使用c++编程语言编码，算法运行环境为：Intel（R）Core（TM）i7-9700 CPU @ 3.00 GHz\times8，16GB RAM.计算时间以秒为单位.

3.1 参数设置与测试实例

为了验证算法的有效性，将原始的ACO算法以及文献［5］中的ALNS算法作为对照算法，并进行大规模实例验证.通过参数调试实验，HACO算法、ACO算法以及文献［5］中的ALNS算法的参数设置如表3所示.为保证计算实验的公平性，本文进行的算例实验均采用此参数设置.

针对本文研究的PDPLDTT问题，基于文献［5］中的实际算例生成本文所需算例.文献［5］的算例是基于福冈、马德里、匹兹堡、西雅图和悉尼5个城市的真实地理数据生成的.基于每个城市生成了6个小规模实例和1个大规模实例，其中，小规模实例包含20位客户和1个配送中心，大规模实例包含100客户和1个配送中心.实例中还包含交付需求、时间窗和其他信息.电动货运三轮车空载质量与驾驶员质量之和设置为140\mathrm{k}\mathrm{g}，电动货运三轮车的载货容量设置为150\mathrm{k}\mathrm{g}，电动货运三轮车的额定功率为350\mathrm{W}.物流的规划时间为120 min.负载等级为10个级别.

由于文献［5］的算例没有取货约束，因此本文为每个实例随机生成了取货需求.生成取货需求的具体规则如下：对于某一个实例而言，为每一个客户点随机生成一个取货需求，其范围在该实例原有的最小的送货需求与最大的送货需求之间.

此外，本文还使用中大型实例对提出的模型和算法进行了验证.中大型实例包括50个客户和75个客户，分别是通过从大型实例中获取前50个和前75个客户而生成的.因此，算例实验总共有30个包含20个客户的小规模实例，5个包含50个客户的中等规模实例，5个含有75个客户的中大规模实例，还有5个包含100个客户的大规模实例，共计45个实例.对于小规模实例，本研究使用了HACO算法以及CPLEX求解器进行求解；对于大规实例，本研究使用了HACO算法、ACO算法以及ALNS算法进行求解.

3.2 算例实验

3.2.1 小规模实例分析

为了检验本文提出模型与算法的有效性，使用CPLEX以及HACO算法解决PDPLDTT问题的小规模实例，其计算结果如表4所示.本文给CPLEX设置了一个求解时间限制，为3 600 s.如果CPLEX在时间限制内解决了算例，那么上下界相等；反之，上界为CPLEX已找到的最好的可行解的目标值，下界为CPLEX在求解过程中通过松弛问题（通常是线性松弛）或其他方法计算出的目标值的下限，这个值表示在当前分支定界算法下，所有可能解的目标值中最低可能达到的值.通过比较上下界，可以衡量当前解接近最优解的程度.上下界之间的差距越小，当前解越接近最优解.本研究使用HACO算法求解PDPLDTT问题10次，以验证算法的稳定性.表4中的Gap为10次求解结果的最好值与CPLEX求解结果的上界之间的相对改进百分比，其计算公式如下：

Gap=\frac{\mathrm{最好}\mathrm{值}-\mathrm{上界}}{\mathrm{上界}}\times \mathrm{100}\mathrm{\%}.\left(\mathrm{28}\right)

从表4中不难看出，CPLEX在3 600 s时间限制下仅仅解决了6个小规模实例，而本文提出的HACO算法在较短时间内（平均求解时间3.4秒）解决了所有小规模实例，能够充分说明HACO算法的高效性.此外，HACO算法不仅具有高效的求解效率，在求解结果准确性上也有一定的保证.表4显示，所有的Gap值均是小于或等于0的，其意义是HACO算法的求解最好值均等于或优于CPLEX的上界.同时，表4显示HACO算法10次求解结果的标准差的平均值为0，说明HACO算法还具有较高的稳定性.

3.2.2 大规模实例分析

为进一步验证HACO算法的性能，本文使用该算法求解PDPLDTT问题的中大规模实例，并与算法ACO以及ALNS的求解结果进行比较，其求解结果如表5所示.表5中的Gap1和Gap2分别为不同算法10次求解结果的最好值之间的相对改进百分比，其计算公式如下：

Gap\mathrm{1}=\frac{\mathrm{H}\mathrm{A}\mathrm{C}\mathrm{O}.\mathrm{最好}\mathrm{值}-\mathrm{A}\mathrm{L}\mathrm{N}\mathrm{S}.\mathrm{最好}\mathrm{值}}{\mathrm{A}\mathrm{L}\mathrm{N}\mathrm{S}.\mathrm{最好}\mathrm{值}}\times \mathrm{100}\mathrm{\%},（29）

Gap\mathrm{2}=\frac{\mathrm{H}\mathrm{A}\mathrm{C}\mathrm{O}.\mathrm{最好}\mathrm{值}-\mathrm{A}\mathrm{C}\mathrm{O}.\mathrm{最好}\mathrm{值}}{\mathrm{A}\mathrm{C}\mathrm{O}.\mathrm{最好}\mathrm{值}}\times \mathrm{100}\mathrm{\%}.（30）

从Gap1可以看出，本文提出的HACO算法有比ALNS算法更好的求解效果，并且HACO算法的求解均值优于ALNS算法，其原因是基于种群的HACO算法有着比基于个体演化的ALNS算法更多的变化性，易于跳出局部收敛状态，因此求解结果的稳定性更高.从Gap2可以看出，HACO算法求解效果优于未经局部搜索优化的ACO算法，结果说明了局部搜索操作可以提高算法的求解能力.从标准差这一指标来看，HACO算法的求解稳定性优于ALNS算法与ACO算法.

本文设计的计算实验证明MILP模型与HACO算法可以有效解决具有负载依赖行驶时间的电动货运三轮车最后一千米同时取送货问题，其求解结果可以为物流企业设计合理的最后一千米物流配送路线，为相关企业提供合理的决策建议.如福冈_01_20实例，其配送方案如表6所示.可以看出该解决方案有三条路径安排，相关企业可以根据路径规划以及出发时间等信息，设计配送工作计划.根据配送计划、负载变化等情况，可绘制如图4的负载变化图，从而更好地展示负载变化的情况.本文提出的数学模型和求解方法可以有效提高决策者的工作效率，并在考虑负载变化的情况下更合理地安排配送规划信息，从而有效降低总体的配送时间.

3.3 HACO算法改进策略有效性分析

本文使用了2-opt、2-opt^*和relocate三种经典的局部搜索算子来提高HACO算法.三种提高策略在算法求解过程中是随机选取的，故很难判断哪种局部搜索算子的改进更为有效.因此，本节对使用的三种局部搜索算子进行有效性分析，使用的算例为福冈_50这一中等规模的算例.具体算法为分别仅仅使用1种局部搜索算子的HACO算法、ACO算法和完整的HACO算法，求解结果见表7.

通过对比发现，三种算子均可以找到最优解，但均容易陷入局部收敛，稳定性差于HACO算法.对三种算子进行比较发现，relocate算子的提升效果强于2-opt和2-opt^*算子，但从均值和标准差来看，此算子优化的HACO算法仍然易陷入局部收敛.

4 结束语

电动货运三轮车在最后一千米配送活动中起着十分关键的作用.然而，在研究最后一千米配送相关问题时，极少有研究讨论电动货运三轮车行驶速度与其负载间的关系.为此，本研究结合现实世界的问题情景，研究了具有负载依赖行驶时间的电动货运三轮车最后一千米同时取送货问题（PDPLDTT），其中行驶速度取决于道路坡度和电动货运三轮车的负载.为求解PDPLDTT问题，本研究首先分析了电动货运三轮车的行驶速度与其负载间的数学依赖关系；其次，对此问题进行了数学建模，构建了一个混合整数线性规划（MILP）模型，并使用了CPLEX求解该模型的小规模案例以验证MILP模型的正确性；再次，为了解此问题的大规模案例，本研究提出了一种使用局部搜索算子优化的混合蚁群优化（HACO）算法；最后，进行了大量的实例计算实验，以验证MILP模型和HACO算法的有效性.计算实验结果表明：本研究提出的MILP模型与HACO算法可以有效解决具有负载依赖行驶时间的电动货运三轮车最后一千米同时取送货问题，其求解结果可以给物流企业设计合理的最后一千米物流配送路线，从而有效提高决策者的工作效率.

本研究在求解PDPLDTT问题时，设计了HACO启发式算法.在后续工作中，可以设计一些基于分解技术（如Dantzig-Wolfe分解、Benders分解）的精确算法来求解负载依赖行驶时间的路径优化问题.此外，还可以进一步引入一些现实因素来丰富PDPLDTT问题.例如，可以考虑交通拥堵因素，将其置于时间依赖的路网环境；也可以考虑随机的取货需求、随机的服务时间等随机因素的影响.

参考文献

[1] 贺冰倩， 李昆鹏， 成幸幸. 快递企业“最后一公里”快件收派优化方案研究［J］. 运筹与管理， 2019， 28（1）： 27-34.

[2] HE B Q， LI K P， CHENG X X. The research on the optimization of courier companies “last-mile” express pickup and delivery process［J］. Operations Research and Management Science， 2019， 28（1）： 27-34.

[3] 王旭坪， 詹林敏， 张珺. 考虑碳税的电子商务物流最后一公里不同配送模式的成本研究［J］. 系统管理学报， 2018， 27（4）： 776-782.

[4] WANG X P， ZHAN L M， ZHANG J. Cost of different elogistics “last mile” delivery mode considering carbon tax［J］. Journal of Systems & Management， 2018， 27（4）： 776-782.

[5] ENTHOVEN D L， JARGALSAIKHAN B， ROODBERGEN K J， et al. The two-echelon vehicle routing problem with covering options： city logistics with cargo bikes and parcel lockers［J/OL］. Computers & Operations Research， 2020， 118［2024-08-23］.

[6] 葛显龙， 温鹏哲， 薛桂琴. 基于需求预测的两级动态配送路径优化研究［J］. 中国管理科学， 2022， 30（8）： 210-220.

[7] GE X L， WEN P Z， XUE G Q. Two-echelon dynamic vehicle routing problem with request forecasting［J］. Chinese Journal of Management Science， 2022， 30（8）： 210-220.

[8] FONTAINE P. The vehicle routing problem with load-dependent travel times for cargo bicycles［J］. European Journal of Operational Research， 2022， 300： 1005-1016.

[9] BEKTAS T， LAPORTE G. The pollution-routing problem［J］. Transportation Research Part B： Methodological， 2011， 45： 1232-1250.

[10] DEMIR E， BEKTAS T， LAPORTE G. An adaptive large neighborhood search heuristic for the pollution-routing problem［J］. European Journal of Operational Research， 2012， 223（2）： 346-359.

[11] KRAMER R， SUBRAMANIAN A， VIDAL T， et al. A matheuristic approach for the pollution-routing problem［J］. European Journal of Operational Research， 2015， 243： 523-539.

[12] HUNG N B， JAEWON S， LIM O. A study of the effects of input parameters on the dynamics and required power of an electric bicycle［J］. Applied Energy， 2017， 204： 1347-1362.

[13] 郭放， 杨珺， 杨超. 考虑充电策略与电池损耗的电动汽车路径优化问题研究［J］. 中国管理科学， 2018， 26（9）： 106-118.

[14] GUO F， YANG J， YANG C. Study on the electric vehicle routing problem in the present of charging strategy and battery consumption［J］. Chinese Journal of Management Science， 2018， 26（9）： 106-118.

[15] FRIEDRICH C， ELBERT R. Adaptive large neighborhood search for vehicle routing problems with transshipment facilities arising in city logistics［J/OL］. Computers & Operations Research， 2022， 137 ［2024-08-23］.

[16] HUANG S H， HUANG Y H， BLAZQUEZ C A， et al. Solving the vehicle routing problem with drone for delivery services using an ant colony optimization algorithm［J/OL］. Advanced Engineering Informatics， 2022， 51 ［2024-08-23］.

[17] COMERT S E， YAZGAN H R. A new approach based on hybrid ant colony optimization-artificial bee colony algorithm for multi-objective electric vehicle routing problems［J/OL］. Engineering Applications of Artificial Intelligence， 2023， 123 ［2024-08-23］.

[18] DORIGO M， MANIEZZO V， COLORNI A. Ant system： Optimization by a colony of cooperating agents［J］. IEEE Transactions on Systems， Man， and Cybernetics， Part B （Cybernetics）， 1996， 26（1）： 29-41.

[19] DORIGO M， BLUM C. Ant colony optimization theory： a survey［J］. Theoretical Computer Science， 2005， 344（2-3）： 243-278.

[20] ZHANG H Z， ZHANG Q W， MA L， et al. A hybrid ant colony optimization algorithm for a multi-objective vehicle routing problem with flexible time windows［J］. Information Sciences， 2019， 490： 166-190.

[21] 苏欣欣， 王红卫， 秦虎， 等. 混合启发式算法求解多配送人员车辆路径问题［J］. 运筹与管理， 2022， 31（2）： 42-47.

[22] SU X X， WANG H W， QIN H， et al. Hybrid heuristic algorithm for the vehicle routing problem with multiple deliverymen［J］. Operations Research and Management Science， 2022， 31（2）： 42-47.

[23] BALSEIRO S R， LOISEAU I， RAMONET J. An ant colony algorithm hybridized with insertion heuristics for the time dependent vehicle routing problem with time windows［J］. Computers & Operations Research， 2011， 38（6）： 954-966.

[24] MAVROVOUNIOTIS M， YANG S X. Ant colony optimization with immigrants schemes for the dynamic travelling salesman problem with traffic factors［J］. Applied Soft Computing， 2013， 13（10）： 4023-4037.

基金资助

国家自然科学基金项目(72404099)

中央高校基本科研业务费基金项目(CCNU23XJ015)

中央高校基本科研业务费基金项目(CCNU23XJ020)

中央高校基本科研业务费基金项目(CCNU24ZZ145)

中央高校基本科研业务费基金项目(CCNU24ZZ146)

中国博士后科学基金面上项目(2023M741302)

湖北省自然科学基金青年项目(2024AFB319)

华中师范大学湖北经济与社会发展研究院2022年立项课题(22HZJS10)