基于密度峰值聚类的VRPTW问题研究

吴 斌，宋 琰，程 晶，董 敏

(南京工业大学 经济与管理学院，江苏 南京 211816)

带时间窗的车辆路径问题(vehicle routing problems with time windows, VRPTW)是车辆调度问题的拓展[1]，广泛存在于快递、商超和外卖配送中，对其研究具有重要的理论意义和现实价值。目前，VRPTW已经被运筹学、计算机、交通运输等领域研究人员关注[2]。VRPTW的优化算法有精确算法和启发式算法。精确算法可以得到精确解，但计算量随着问题规模增大呈指数增长。启发式算法拥有强大的全局搜索能力，因而对于大型复杂问题，只能使用启发式算法[3-5]。目前相关的文献已经有很多。文献[6]提出了混合粒度禁忌搜索算法求解VRPTW问题，该方法虽然局部搜索能力强，但时间复杂度高且结果严重依赖初始解。文献[7]克服了蚁群算法在最后阶段容易产生不可行解的弊端，建立了良好的积极反馈机制，但算法耗时且容易停滞。文献[8]利用人工蜂群算法解决VRPTW问题，并且算法收敛速度很快。文献[9]使用改进的遗传算法解决具有强制回程特征的双目标车辆路径问题，结果显示算法具有良好的全局搜索能力，计算速度快，但难以获得全局最优解。文献[10]研究DVRP (dynamic vehicle routing problem)问题，提出混合粒子群和变邻域搜索的求解方法，在DVRP场景里(维修服务、快递)能够得到有优势的解。Jiang 等[11]使用最佳客户插入原则PFIH(push forward insertion heuristic)构建初始解，再基于SA (simulated annealing)和LNS (large neighborhood search)的混合策略来优化初始解, 并且用回归迭代策略计算每辆车的最佳离开时间。Xu等[12]应用基于迭代次数的线性递减函数为粒子群算法提供全局和局部搜索能力之间的平衡，并且结合遗传算法的交叉操作避免算法早熟收敛。

上述研究侧重算法设计，缺乏针对问题特征的求解。在现实世界的配送中，如外卖、快递、商超配送等很多场景都是客户聚集的。如何通过聚类方法，减小问题规模，提高计算速度，获得更强的全局寻优能力，是亟待解决的问题。

2014年Rodriguez等[13]在提出了一种新型聚类算法——DPC (density peak clustering)。该算法能快速发现任意形状数据集的密度峰值点 (即聚类中心)，并高效进行样本点分配和离群点剔除。DPC算法在聚类过程中较少需要人为干预，能大幅提高聚类结果准确性。目前DPC已在一些领域成功应用。邹旭华等[14]将DPC算法运用于彩色图像分割，并且具有较好的分割效果。Chen等[15]应用DPC确定面部图像中年龄族群的密度峰值，再根据面部照片到密度峰值的距离来预测照片主人的年龄。文献[16]用DPC结合MDNS (multi-document news summarization)模型生成具有较小冗余的新闻摘要，解决了BOW(bag of words)模型存储量高，计算复杂的问题。DPC已经在计算机视觉、自然语言处理和社交网络[17-18]等领域成功应用，但在物流领域的应用还非常少。

综上所述，针对目前算法在处理大规模VRPTW问题时易陷入局部最优、计算时间长等弊端，本文提出一种基于密度峰值聚类 (DPC)与遗传算法(Genetic algorithm, GA)相结合的混合优化算法DGA (density peak clustering with genetic algorithm)。通过对客户进行聚类，将大规模问题划分为小规模的子问题，然后用遗传算法对每个子问题求解。

1 数学模型

1.1 问题描述

假设物流配送网络由一个配送中心(用0表示)和N个节点组成，其中N∈Z+。车辆均停靠在配送中心，配送中心提供货物并对车辆进行统一调度规划，车辆根据规划线路将货物运送至各个节点。对VRPTW的配送流程作如下定义。

1) 各节点的需求量已知，且每个节点的需求量mi小 于可用车型的最大容量qk。

2) 所有车辆从配送中心出发，完成配送任务后，需要返回配送中心，且车辆不能超载。

3) 配送中心拥有足够的货物来满足客户需求。

4) 所有的节点只能被某一辆车服务一次。

1.2 VRPTW模型[19]

决策变量如下。

ti为到达节点i的时间。

wi为在节点i的等待时间。

xijk∈0,1： 车辆k从i到j为1，否则为0。ij;i,j∈0,1,2,···,N。

参数如下。

K为配送车辆总数量；

N为节点总数量；

yi为任意实数；

dij为i,j节点之间的欧氏距离；

cij为i,j节点之间的配送成本；

mi为i点的需求；

qk为车辆k的总容量；

ei为 节点i的最早到达时间；

li为节点i的最晚到达时间；

fi为 在节点i的服务时间；

rk为车辆k的最大行驶时间。

目标函数为

约束条件为

模型中，0点为配送中心，车辆从配送中心出发。式(1)为目标函数，目的是最小化总配送成本。式(2)～(10)为模型需要满足的约束条件，其中，式(2)表示车数量约束；式(3)确保了每辆车从配送中心出发并返回配送中心；式(4)～(5)保证每个节点只能被一辆车服务，且只能服务一次；式(6)为车容量约束；式(7)为车辆k的最大行驶时间约束；式(8)~(10)为节点的时间窗约束。

2 密度峰值聚类结合遗传算法求解VRPTW问题

2.1 密度峰值聚类算法

DPC算法能够发现任意形状的簇，样本点无需迭代，调节参数少，思路直观明快，易于理解。该聚类算法的核心思想在于对聚类中心的刻画。Rodriguez等[13]认为聚类中心同时具有以下2个特点：1) 聚类中心点的自身密度很大；2) 聚类中心与比其密度大的数据点的距离相对较远。DPC主要有2个需要计算的量，局部密度ρi和相对距离δi。

定义1[13]当数据点离散时，

当x＜0时 ，χ(x)=1，否则，χ(x)=0。当数据点连续时，

其中，dij表示i,j两点之间的欧氏距离，dc为截断距离。

定义2[13]相对距离δi，

对于局部密度ρi最大的样本i，δi=maxjdij。

根据 ρi和δi选择聚类中心点，计算 γi=ρiδi。选取 γi值大的点作为聚类中心，这些点具有较大的局部密度和相对距离，具有作为聚类中心的特征。

DPC的不足在于聚类效果依赖截断距离dc的设定。文献[20]提出基于基尼系数自适应调整截断距离dc。Wang等[21]结合数据场理论，使用原始数据集的熵获取截断距离，提出的计算点势能的方法如式(14)所示。在一个数据集中，数据密集区域的点势能较大，稀疏区域的点势能较小。这表明数据域的势能与点的局部密度具有类似的效果。对于数据集{x1,x2,···,xn}, 每个点的势能计算公式为

式(14)类似Gaussian kernel的计算，其中∥x−xi∥代表欧氏距离，σ是需要确认的变量值。对于势能集{φ1,φ2,···,φn}，定义数据域的熵值H为

图1 σ值优化图Figure 1 σ-value optimization graph

从图1可以看出，随着σ 不断曾大，熵值先减小再增大，故可以取得熵最小时的σ 值。根据Gaussian分布3B原则[22]，每个点的影响半径为作为聚类算法，一个点只能影响位于其影响半径内的点。故取σ 为5.4437时的截断距离计算后再根据 γ−Point决策图选取聚类中心，选取过程更加简单直观。

2.2 VRPTW的遗传算法优化方法

2.2.1 编码和解码

针对本文问题，采用简单直观的序列编码方法。将客户点编号为1 ,2,···,N，然后根据此编号设计编码，随机生成长度为N的P个染色体，染色体中的每个基因为编号中的某个编号，代表一个客户，每一个染色体对应一种配送方案。解码按染色体顺序进行，从第1个基因开始分配车辆，直到该车满载或不能满足下一个客户的需求为止，然后安排下一辆车，最终该染色体所对应客户全部得到车辆分配。

2.2.2 遗传操作

交叉和变异是遗传算法操作算子的基本形式，通过交叉和变异操作，可以不断搜索靠近最优解。本文采用适用于车辆路径问题的部分匹配交叉(partially matching crossover, PMX)[23]法，具体交叉过程如图2所示。

图2 PMX交叉过程Figure 2 PMX crossover process

Step1在父代染色体中随机选取两个交叉点；

Step2交换两个交叉点间的基因，图中为基因3、4、5、6与6、9、2、1交换；

Step3从父代中继承未选定交叉的基因7、8；

Step4对剩余位置做映射及冲突检测。根据交换的2组基因建立一个映射关系，形成1 ↔6↔ 3、2↔ 5 、9↔ 4三组映射关系。父代1中基因1经过映射变为基因3(基因6存在冲突)，基因2变为基因5，基因9变为基因4。用相同的方法对父代2做映射，最终形成图2中的2个子代。

由于本文为整数编码，所以不同于二进制编码的变异方式。变异步骤如图3所示。首先选择一个父代染色体随机选择2个点作为起止点s和e，然后将起止点之间的基因反转，生成新的子代染色体。

图3 变异过程Figure 3 Mutation process

2.2.3 遗传算法流程及约束处理

Step1根据序列编码规则随机产生P个染色体以满足约束(4)~(5)，构成初始种群，设置迭代次数L=1。

Step2对染色体解码，把对不满足车容量、车辆最大行驶时间以及时间窗约束的染色体删除。评估第L代的种群，设计的适应度评价函数为

找到最优解及适应值，转Step 3。

Step3判断迭代次数L是否到达最大迭代次数Lmax，是则转Step 6。

Step4依据轮盘赌原则从种群中选择染色体作为父代染色体。

Step5采取复制，交叉，变异操作直至产生P个染色体，形成新的种群，转Step 2。

Step6获取最优解方案，输出结果。

2.3 DGA求解VRPTW问题的流程

1) 利用密度峰值聚类对客户点进行聚类，将位置相近的客户安排在同一个类簇中；2) 对每一个类簇应用遗传算法求解，求得的解根据车容量和客户需求约束来对车辆进行分配。这样不仅可以最大限度的减少车辆的行驶里程，还能够保证车辆能够严格地满足客户的时间窗要求。DGA算法流程如图4所示。

3 仿真实验

3.1 DPC算法进行客户聚类

选择Solomon benchmark中的C1数据集进行仿真实验。该数据集中每个实例都包含100个客户、25辆车，每辆车的容量为200。本文算法运行环境为64位Windows 10系统，处理器为Intel(R)Core(TM)i7-7700HQ 2.80 GHz，内存为8 G，编程语言为Python。1) 用DPC对C101数据集的100个客户的位置坐标进行聚类，让地理位置相似度高的客户安排在一起配送，简化配送流程。采用熵估计的方法计算C101数据集的聚类参数——截断距离dc，如图5所示，求得使熵H最小的σ 值 3.8429，与此对应的最优dc≈8.1519。2) 进行聚类中心的选择，作决策图，得到聚类中心点为10个，如图6所示。3) 对非聚类中心点进行归类，结果如图7所示，客户点被聚类到各个类簇，其中蓝色的点为配送中心，坐标值为(40, 50)，黑色点为各个类簇的聚类中心。

图4 DGA算法流程Figure 4 DGA algorithm flow

图5 熵优化Figure 5 Entropy optimization

图6 决策图 (dc=8.1519)Figure 6 Decision graph

图7 客户聚类结果(C101)Figure 7 Customer clustering results (C101)

对10个客户类簇进行可视化(图8)，图中横坐标为点索引，纵坐标为该点到对应聚类中心的距离。可以看到，簇中点到聚类中心的距离相对平均，没有距离差很大的客户点需要去安排配送。

图8 各类簇中样本点到聚类中心的距离Figure 8 The distance from the sample points in each cluster to the cluster center

3.2 遗传算法优化配送线路

Solomon数据集中假设车速为单位距离，这样可以使得各节点之间旅行时间tij和节点之间的欧氏距离dij相等。设置遗传算法交叉率为0.8，变异率为0.1，种群规模为100，迭代500代。以C101数据集为例，将遗传算法运行10次获得最优值，具体路径及车辆到达每个客户点的时间见表1(车辆列括号内为车装货量)，算法所求得最优解的车辆路径图如图9所示，图中每种颜色代表一辆车。最优值为828.94，车辆平均满载率为91%。

表1 DGA计算C101各车辆行驶路径及到达各客户点时间Table 1 Driving path and arrival time of C101 vehicles calculated by DGA

图9 各类簇车辆路径Figure 9 Vehicle paths of various clusters

3.3 算法比较

运用DGA对其他C1类型数据实例进行仿真实验，并与其他算法进行对比。各实例计算结果如表2所示，可视化如图10所示。从平均值角度分析，DGA的解优于OV[24]、SA[19]和Tabu[19]算法得出的解，并且接近于最优解。

文献[25]中的最优解是在禁忌搜索中引入精确算法得到的。将禁忌搜索算法找到的优秀解(350~450个)保存在集合T中，提出集合划分问题(该问题也是NP-难问题)，然后通过调用精确求解器Cplex来求解，得到了目前已知的最优解。该方法虽然找到了最优解，但算法运行时间较长。本文提出的方法与最优解之间有偏差，主要是在C103和C104这2个问题上未找到最优解。通过比较最优解和DGA的优化结果，发现每条线路中的客户是相同的，仅客户的配送路径不一致。这是由于C103和C104是时间窗十分宽松的VRPTW问题。此类问题配送线路的可行解空间更大，需要较强的路径优化算法。为突出聚类对算法的影响，本文仅采用基本遗传算法进行路径优化，没有对遗传算法进行改进或者引入启发式算法进行后优化，造成线路的优化结果有偏差。

图10 各算法计算结果可视化Figure 10 Visualize the calculation results of each algorithm

表2 各算法结果对比Table 2 Comparison of results of various algorithms

从表2所求得的最优值来看，OV、SA、Tabu算法在C102、C103和C104数据集中均会陷入局部最优，但本文算法在这3个数据集上仍能得到最接近最优值的解，最大误差为0.91%，说明本文提出的DGA算法具有较好的鲁棒性。聚类可以将原始数据集拆分成多个小数据集，通过减小问题规模来提高解的质量。表3记录了C102、C103 和C104三个实例聚类后每个类簇所要配送的客户数量。通过对比发现，聚类使原本的100个客户的优化问题，变成每个类簇进行10个客户左右的线路优化，大大减少了计算量和计算难度。因此，DGA在求解大规模问题时具有较大参考意义。聚类保留了客户的地理位置聚集特征，不仅能得到比GA[19]更快更好的解，而且得到的解在实际中更加方便物流配送。

表3 C102、C103、C104聚类计算结果Table 3 C102, C103, C104 cluster calculation results

4 结束语

本文提出密度峰值聚类与遗传算法相结合的混合算法求解VRPTW问题。通过聚类缩小遗传算法求解问题的规模，提高了结果的准确性，同时还加快了求解速度。该方法克服了传统遗传算法容易陷入局部最优的缺点。通过实例验证，并与其他算法进行对比，解的质量最多提升了26.4%，证明了DGA在处理VRPTW问题时具有一定的优越性。DGA在求解大规模问题时可以加快问题的求解过程，扩大了VRPTW问题的研究规模。但聚类算法在Solomon随机数据集上的效果不太好，因为客户点位置随机，DPC不能发挥最大效用。因此，如何针对随机客户数据集改进聚类算法是下一步需要研究的问题。