0 引言

近年来, 随着互联网经济的蓬勃发展, 相关行业对快速、高效的物流配送服务的需求愈加迫切。我国物流行业的特点是起步晚、需求大、发展快, 在粗放式的发展过程中遗留下管理欠缺规范, 配送效率低下等问题[1]。利用计算机技术实现物流配送作业计划的合理制定, 是解决前述问题的有效手段。本文以单仓储中心的多车物流配送作为应用场景, 则多台配送车辆的路径规划问题可建模为图论中经典的车辆路径问题 (Vehicle Routing Problem, VRP) [2]。由于实际生活中配送任务往往带有一定的时间限制, 所以又出现了带时间窗的车辆路径问题 (Vehicle Routing Problem with Time Windows, VRPTW) 。VRPTW是在VRP的基础上增加客户要求访问的时间窗口, 且时间窗可分为软硬两种:硬时间窗VRP指每项任务必须在要求的时间范围内完成, 否则得到的解为不可行解, 而软时间窗VRP则是在超出要求时间范围时, 给予一定的惩罚[3]。VRP具有NP完全问题 (Non-deterministic Polynomial Complete, NPC) 求解复杂度[4], 因而不存在多项式时间算法, 一般以模拟退火[5]、禁忌搜索[6]等启发式算法或以遗传算法[7]、蚁群算法[8]等智能算法进行求解。

粒子群优化 (Particle Swarm Optimization, PSO) 算法是智能算法的优秀代表[9], 因其简洁的内在寻优机理与快速的收敛速度而被广泛应用于VRP以及VRPTW的求解, 并由此衍生出许多改进算法[10,11,12,13,14,15]。李宁等[3]将粒子种群划分为若干个部分重叠的子群, 各粒子根据所属子群内的最优解调整其位置。吴耀华等[10]将子群的概念同近邻机制[11]相结合提出局部近邻机制, 使粒子在更新速度时借鉴局部近邻粒子的优势, 提高算法收敛速度。马炫等[4]在粒子群算法的基础上引入粒子交换原理, 提出了一种整数粒子更新方法。马冬青等[12]为了增强算法的局部搜索能力, 在算法中增加了爬山操作。

以上研究大部分都是对算法中参数和调节因子的改进, 很少考虑粒子表达方式存在的问题。其中文献[3,4,10]均采用了双层的粒子模型, 按照先分组再排序的规则进行编码, 第一层表示各任务点对应的车辆, 第二层表示各任务点在对应的车辆运送路径中的优先顺序。由于分组与排列之间具有强耦合关系, 第一层编码的变化势必影响第二层编码更大的变动, 同时由于求解过程中直接套用公式, 易产生小数, 解码时必须进行取整和排序, 导致搜索效率低、不可行解数量多等问题。为解决上述编码存在的问题, 本文构造了一种适用于解决VRPTW的新型粒子表达方式———排位编码方式, 将每个站点的后继站点所组成的序列作为编码, 改变后继即可改变车辆路径;并依据编码特点设计相应的粒子更新方式, 避免因取整和排序造成的不可行解数量多、求解效率低等问题。

1 优化模型

本文仅考虑单仓储中心的多车物流配送场景, 对多车配送作业计划进行如下约定:

1) 各配送车辆以仓储中心为发车起点并返回;

2) 各配送点均须被访问且至多被访问一次;

3) 各车配送任务量不得超过车辆载重上限;

4) 各配送点间任意可达, 任意两点距离以欧氏距离计算。

在以上约定基础上, 若仓储中心有K辆车, L个配送点, 以站点1表示仓储中心, 客户站点2~L+1表示各配送点 (此后为便于区分将配送点称为客户站点) , 那么可将单仓储多车物流配送优化问题描述为:仓储中心用K辆车对L个客户站点进行配送, 各配送车辆的最大载重均为q;客户站点i的配送量为gi (i=2, …, L+1) 且max gi≤q;设完成客户站点i的配送任务需要的时间为Ti, 任务开始执行的时间范围为[ETi, LTi], 其中ETi为允许最早开始时间, LTi为允许最晚开始时间。若车辆到达站点i的时间早于ETi, 则需在i处等待;若车辆到达站点i的时间晚于LTi, 则任务需延迟执行[10]。求能够最小化车辆配送总成本且满足所有客户站点需求的配送作业计划。其中, 配送总成本可指代距离、费用、时间等一种或多种成本形式, 本文仅将其作为一种抽象的概念, 不作具体意义的解读, 但这并不影响问题的性质。

问题用数学语言描述如式 (1) ~ (5) :

上述模型中, cij表示车辆从点i行驶至j的行驶成本, 显然cij=cji。当车辆k完成从点i至j的配送时, xijk=1;否则, xijk=0。当点i的配送由车辆k完成时, yki=1;否则, yki=0;pE, pL分别为单位时间的等待成本和迟到罚金。式 (1) 中的目标函数min Z确保配送总成本最小。约束中, 式 (2) 为货车容量限制, 每辆车的配送总量不得超过货车载重量;式 (3) 保证每个客户站点只能由一辆车配送一次, 不能重复配送;式 (4) 和式 (5) 则保证所有客户站点都被配送到。

2 粒子群算法简介

PSO依赖于一个种群, 通过种群中各粒子在解空间中的协同飞行实现对问题潜在解的搜索。在一个n维搜索空间中, 粒子i的位置可表示为一个n维向量Xi, 其速度以向量Vi表示, 所经历的最优位置记为Pbest_i, 种群中所有粒子经历过的最优位置记为Gbest。粒子飞行过程中的速度和位置更新方法以式 (6) 、 (7) 描述:

其中:ω为惯性权重;c1和c2为加速因子, 亦称学习因子;r1, r2~U (0, 1) 。

粒子群算法的计算流程分为以下几步:

步骤1随机初始化各粒子的位置和速度, 将各粒子的初始位置作为Pbest_i, 以目标函数式 (1) 为依据, 对各粒子进行适应度的评价, 把种群最优个体作为Gbest;

步骤2根据式 (6) 、 (7) 分别更新各粒子的速度和位置;

步骤3计算更新后各粒子的适应度值并更新Pbest_i和Gbest;

步骤4检验是否满足算法终止条件:若满足, 终止迭代;反之, 返回步骤2。

3 求解VRPTW的粒子群算法

3.1 编码与解码

本文在构造编码方式时, 采用一种基于后继的排位编码方式, 先将所有站点按序号从小到大排列, 然后根据各站点实际配送顺序, 在每个站点下填入其后继站点。这样所有后继站点组成的编码就构成了一种粒子表达方式。

例如, 当使用3辆车为9个客户站点服务时, 车辆路径为:1-4-9-5-1;1-2-7-6-10-1;1-3-8-1。编码结果如图1。

图1 编码结果   

由于3辆车都从配送中心出发, 所以配送中心1的后继站点有3个, 而其余站点都只有一个后继站点。编码中1出现的次数等于车辆数, 表示所有车辆最后都回到配送中心, 2~10均只出现一次, 表示每个客户站点都被配送且只配送一次。

解码时, 以配送中心下的每个站点作为每辆车的第1个客户站点, 分别以它们为序号进行索引, 找到的对应元素就是需要配送的第2个站点, 再以第2个站点为序号继续索引, 直至最后所有车辆都回到配送中心。

这种编码方式有以下特点:首先选用整数编码, 与车辆路径问题的离散性相适应;其次从某个元素所对应的序号和该元素的大小可以断定它们代表的两站点之间存在一条车辆路径, 那么改变某个序号位置下的元素, 就改变了对应站点的后继站点, 所以可以通过变换编码, 实现车辆路径的改变。车辆路径的改变相当于粒子群算法中粒子的位置发生改变, 这一点又与算法的要求相吻合。

3.2 适应度函数和约束的处理

本算法要寻找的目标解是总成本最小的行驶方案, 而总成本中行驶成本占的比重最大, 所以在确定适应度函数时, 行驶成本是主要评判标准。其次考虑约束:针对时间窗约束, 采用增加等待成本和迟到罚金的手段;针对车辆的载重量约束, 引入外部罚函数, 赋予不可行解或非法解一个极大值, 将有约束最优化问题转化为无约束最优化问题。设某辆车分配的总任务量, 那么适应度函数写成表达式的形式如式 (8) :

式 (8) 中S是粒子当前位置的适应度值, 等号右侧前三项分别对应行驶成本, 等待成本和迟到罚金, 即目标函数, 最后一项为超载成本, p为单位超载的惩罚因子。

群体中粒子的位置每次发生变化后, 都要进行适应度的计算, 并将其与粒子历史个体最优进行比较。若当前位置的适应度小于历史个体最优, 则用当前位置替换历史个体最优, 否则保留。再通过比较所有粒子的历史个体最优, 找出适应度最小的替换种群的历史全局最优。

3.3 粒子更新

由于采用了全新的粒子编码形式, 粒子速度与位置的更新无法直接套用式 (6) 、 (7) , 因而需要在对粒子速度与位置更新公式理解的基础上, 重新设计相应的算子, 以实现粒子种群在解空间中的协同搜索。

对PSO粒子更新公式进行分析, 令G1=ω·Vi (t) 、G2=c1·r1· (Pbest_i (t) -Xi (t) ) 、G3=c2·r2· (Gbest (t) -Xi (t) ) , 则式 (6) 可简写为V=G1+G2+G3, 其中认知项G2和社会学习项G3决定了粒子对解空间的开采能力, 是粒子以向更优个体趋近的方式实现优化的一种学习机制;惯性项G1则确保了粒子对解空间未知区域的勘探能力。不失一般性, 一个群智能优化算法具备优秀寻优能力的前提是能够很好地平衡开采和勘探。只注重开采而缺少勘探, 算法表现为局部收敛;只注重勘探而缺少开采, 算法表现为盲目地随机搜索。

在对PSO粒子更新公式进行上述分析的基础上, 本文利用已有的编码特点, 设计了基于“学习”和“变异”的粒子更新方法, 具体操作步骤如下:

3.3.1 学习

如前所述, 开采本质上是粒子以向更优个体趋近的方式实现优化的一种学习机制。当前粒子向更优个体的趋近过程具体表现为两粒子编码差异的逐步缩小。据此, 学习机制中只考虑G2和G3项, 忽略G1项, 同时定义k1, k2为学习因子, 其作用同c1·r1, c2·r2, 那么式 (6) 、 (7) 可合并改写为式 (9) :

式 (9) 中涉及三种运算。本文对这三种运算定义如下:

1) 表示确定x中不同于y的元素及其位置信息;

2) 表示从x中随机选取部分元素, 且该部分元素占x中总元素个数的k倍;

3) 表示依据y中元素及其位置信息修正x的对应部分。

说明以上运算的具体实现步骤:

步骤1确定Pbest_i (t) 和Gbest (t) 中不同于Xi (t) 的部分, 即:和;

步骤2根据学习因子k1, 从中随机提取对应个元素, 完成的运算, 的运算同理;

步骤3根据定义3) , 以G2所得结果修正Xi (t) 中的对应部分, 得到新编码Xi (t) ', 再以G3所得结果修正Xi (t) ', 得到Xi (t+1) 。

需要强调的是, 若将G2所得结果中的站点称作目标点, Xi (t) 中对应位置的待修正站点称作初始点, 那么步骤3中“修正”具体指:将Xi (t) 中的初始点更换为G2所得结果中相应位置的目标点, 同时将Xi (t) 中原本等于目标点的元素更换为初始点, 实际操作时将Xi (t) 中等于目标点的元素与初始点交换位置即可。此外在执行步骤2时需检测是否有自环产生。所谓自环, 就是指几个客户站点独立成环, 没有遵循从配送中心出发并回到配送中心的原则。自环产生的原因在于目标点和其对应的初始点位于同一辆车的配送路线上。例如, 某辆车的配送路线为1-2-4-7-5-6-1, 当目标点为6, 初始点为4, 将4和6互换后, 这条线路被拆分为1-2-6-1和7-5-4-7-…, 这就导致了7, 5, 4这三个站点独立成环。为了避免自环, 在选定目标点时需检验其是否存在于初始点所在的线路中, 若存在, 则需重新选择目标点, 并再次检验。

举例当前粒子向个体最优粒子学习的操作:

Pbest＿i (t) 表示的车辆路径为:1-2-4-7-1;1-3-5-6-1。

Xi (t) 表示的车辆路径为:1-2-7-1;1-5-3-4-6-1。

首先找出两粒子编码中的不同之处, 共5处, 若k1=0.2, 则需要随机选取的元素个数为1。若选择第2站点的后继进行学习, 那么Xi (t) 中第2站点的后继7为初始点, Pbest＿i (t) 中第2站点的后继4为目标点。经检验在Xi (t) 中站点7和站点4不在同一条路径上, 所以符合交换条件。接下来将Xi (t) 中的7与4进行交换, 得到Xi (t) '的编码, 学习过程中粒子的编码变化如图2所示。

经过学习之后, Xi (t) '表示的车辆路径为:1-2-4-6-1;1-5-3-7-1。车辆路径从之前的配送完第2站点配送第7站点, 变为配送完第2站点配送第4站点, 这一路线的改变即是借鉴了Pbest＿i (t) 的结果。同理, 再以G3所得结果修正Xi (t) '即可得Xi (t+1) 。这种更新方式操作简洁, 且计算过程中不会产生小数, 解码便利, 同时避免了因取整和排序而出现不可行解的现象, 提高了算法的搜索效率。

图2 学习过程中粒子的编码变化 

3.3.2 变异

学习机制赋予粒子开采能力, 粒子之间通过互相学习使种群逐渐收敛, 虽然过程中会产生不同于初始解的新解, 但其搜索范围始终局限于初始解空间, 容易陷入局部最优解。为了扩大种群搜索范围, 需要粒子具备空间勘探能力, 即式 (6) 中的惯性项G1, 本文通过在迭代过程中引入变异机制来实现这一功能。变异机制中只考虑G1项, 忽略G2和G3项, 同时定义Pm为变异率, 那么式 (6) 、 (7) 可合并改写为式 (10) :

式 (10) 中的运算同学习机制中的定义, Vi (t) 为随机生成的粒子编码, ω为变异强度, rnd为区间 (0, 1) 上的随机数。

由于Vi (t) 是随机生成的, 所以在实际操作中借鉴学习机制, 随机交换某两个站点的后继站点, 即可完成变异操作。变异过程中同样需要避免自环的产生, 但这里使用和学习机制中不同的方法。因为在学习机制中, 更优粒子的编码是确定的, 当目标点选定时其对应的初始点也是确定的, 但在变异机制中, 初始点和目标点都是随机的, 所以可以先选定初始点, 然后把初始点所在线路上的所有客户站点从全部站点集合中剔去, 再从剩余站点中选取目标点直接进行交换, 从根源上杜绝自环的产生。

举例粒子的变异操作:

Xi (t) 表示的车辆路径为:1-2-3-7-1;1-4-8-1;1-6-5-1。

编码中共包含10个元素, 若ω=0.1, 则需要交换的元素对数为1。选择第3站点的后继站点7作为初始点, 将站点7所在线路上的所有客户站点即2, 3, 7从全部站点集合中剔去, 然后从剩余站点中随机挑选一个, 比如站点8作为目标站点, 进行交换, 得到Xi (t+1) 的编码。变异过程中粒子的编码变化如图3。

Xi (t+1) 表示的车辆路径为:1-2-3-8-1;1-4-7-1;1-6-5-1。

图3 变异过程中粒子的编码变化   

3.4算法流程

针对VRPTW优化模型的特点, 本文首先构造编码使问题更加适应算法, 然后根据实际问题的目标和约束设置适应度函数, 之后在编码方式的基础上设计了基于“学习”和“变异”的粒子更新方法, 在保障种群收敛的同时扩大搜索广度, 避免陷入局部最优, 最终找出全局最优解。粒子位置更新的完整数学表达如式 (11) 。本文所提算法可以图4所示的流程描述。

图4 求解VRPTW的粒子群算法流程  

4 仿真实验

为了验证算法的有效性并考察相关参数对算法性能的影响, 以Matlab对算法进行了编程实现。仿真平台为Intel Core i5-5200U CPU@2.2 GHz, 4 GB RAM。

4.1 仿真算例

采用文献[4]中20个任务点的算例, 对各站点进行重新编号, 配送中心用序号1表示, 20个客户站点用2到21表示, 各站点的要求如表1。车速为50 km/h, pE=pL=50, p=100, 单位运输成本为1。使用5辆车进行配送, 每辆车的载重量都为8 t。

表1 站点列表 

为了比较, 文中算法取与文献[4]相同的种群规模和迭代次数分别为100和1 000, 同样运行10次。本文算法参数设定如下:k1=k2=0.3, Pm=0.6, ω=0.15。对比结果如表2所示。本文最优解的车辆路径为:1-2-5-13-1;1-9-6-18-15-16-7-1;1-21-19-4-1;1-11-8-20-1;1-3-12-14-17-10-1, 如图5 (a) 所示。文献[4]最优解的车辆路径为1-21-13-1;1-2-3-5-6-4-1;1-9-12-14-10-7-1;1-19-17-18-15-16-1;1-11-8-20-1, 如图6 (b) 所示。

表2 本文和文献[4]算法实验结果对比  

图5 本文和文献[4]算法得到的车辆路径比较 

观察表2中的数据和图5中的车辆路径可以发现, 参考文献[4]中最优解的车辆行驶总路程较短, 所以行驶成本低于本文中的最优解, 但其等待成本过高, 也就是说文献[4]中的算法用过高的等待成本换取较低的行驶成本。而本文提出的算法则是在减小行驶成本的同时控制等待成本的增长, 所以得到的最优解优于文献[4]算法。

由此可见, 本文提出的编码方式与粒子群算法具备较好的寻优能力, 能够有效解决单仓储中心情景下的多车物流配送优化问题。

4.2 实验中参数的选取

4.2.1 学习因子的选取

学习因子在本算法中代表了粒子向个体最优和全局最优学习的强度, 用当前粒子向更优粒子学习的元素个数占两粒子海明距离的比例来表示。采用4.1节的算例, 保持算法中其他参数不变, 选取不同的学习因子进行实验, 为减小偶然因素导致的偏差, 每种情况运行程序20次并计算出平均值, 对算法收敛后的适应度值的分布进行三阶最小二乘拟合, 结果如图6 (a) 。

观察图形曲线可以发现, 当学习因子位于0.25~0.3时种群收敛后的适应度值达到最低, 所以, 学习因子在[0.25, 0.3]区间内取值时, 种群收敛于全局最优的概率最大。

4.2.2 变异率和变异强度的选取

变异率和变异强度的作用同标准粒子群算法中的惯性权重, 都是为了使算法在勘探和开采上取得平衡。在本问题中, 变异率指种群中发生变异的粒子数占粒子总数的比例, 变异强度指变异的粒子中发生交换的站点对数占总站点个数的比例。同样采用4.1节的算例, 实验过程和分析如下:

观察变异率对算法精度的影响时, 保持算法中其他参数不变, 对变异率取不同值, 分别运行程序20次并计算出平均值, 对算法收敛后的适应度值的分布进行四阶最小二乘拟合, 结果如图6 (b) 。观察变异强度对算法精度的影响时, 对变异强度取不同值, 结果如图6 (c) 。

图6 几种实验参数对算法精度的影响   

从图6 (b) 中可以看出, 变异率在[0.5, 0.7]区间内时, 种群收敛后的适应度值达到最小, 因此推断在该问题中将变异率设置在0.6附近时算法在搜索广度和深度上取得平衡, 求得最优解的概率较大。变异率的取值应遵循以下规律:如果种群较快收敛至局部最优, 则需要增大变异率, 使更多的粒子飞出当前聚集的区域, 去更远的解空间中寻找是否存在更优解。如果种群中的粒子运动过于分散, 长时间不收敛, 就需要减小变异率, 加强局部搜索, 保证种群在可接受的时间内收敛。

从图6 (c) 中可以看出, 当变异强度在0.2附近时, 算法求得最优解的概率较大。说明就当前种群的规模而言, 粒子的变异强度在[0.15, 0.25]区间内时, 取得比较适合的飞行步长, 太大则易导致开采不够细致, 太小则造成勘探能力不足。

总结实验结果, 算法中参数的选取具有一定的规律可寻, 当参数取值从小到大变化时, 算法收敛后的适应度值总是先减小后增大, 并且能够在一个较宽的区间内取得最优。说明算法收敛精度受参数影响相对较小, 鲁棒性较强。

5 结语

本文提出了一种适用于解决带时间窗车辆路径问题的粒子群优化算法。首先构造排位编码方式将算法与问题紧密联系起来, 其次依据PSO原理设计了基于“学习”和“变异”的粒子更新方法, 实现粒子种群对解空间的开采和勘探。通过学习使种群中的粒子在认知和社会部分的影响下逐渐收敛, 而变异则增强粒子在解空间中的勘探能力, 避免陷入局部最优解, 最终在可接受的时间内找出令人满意的解。

针对实际生产生活中的车辆路径问题可以在此算法的基础上根据不同要求进行优化, 例如进一步提高解的质量:可以将粒子群划分为若干个两两互相重叠的相邻子群, 从而进一步减小算法陷于局部最优解的可能;增加更多的约束, 使之更符合实际问题:考虑站点之间的实际距离和可达性;增加算法可解决问题的场景:考虑将算法应用到多仓储中心的物流配送优化中等。