一、VRP的介绍
1.1 What is VRP？
车辆路径问题（Vehicle Routing Problem,VRP）最早由Dantzig和Ramser于1959年首次提出[1]，可用语言描述为：设有一集散中心(central depot），配备具有相同特征的K辆车为n位客户（customers）提供送货服务 ，每位顾客有其特定的需求量D（demand）。车辆从场站出发对客户进行配送，最后返回场站，要求所有客户都被配送，每位客户一次配送完成并且其demand被满足。目的是使总成本最小，也就是所有车辆行驶路线的总距离最小。
车辆路径问题一直是网络优化问题中最基本的问题之一，一直受到国内外学者的关注。图1-1展示的便是一个VRP网络示意图。
图1-1 VRP示意图
1.2 VRP的常见类型
历经数十载发展，VRP延伸出了非常多的变体问题，针对不同的适应场景和现实问题，这些变体各自具有鲜明的特征。针对初学者，本文仅列举出VRP的5种常见类型，各类型之间的关系如图1-2所示。
（1）具有容量限制的VRP（Capacitated VRP,CVRP）：任意车辆路径的的总需求不能超过车辆的最大容量；
（2）带有距离限制的VRP（Distance-Constrained-VRP,DVRP）：每条路径的距离（或车辆的行驶时间）是有限的，问题的目标是最小化总路线长度或者最小化车辆的服务持续时间（duration）；
（3）带有时间窗的VRP（VRP with Time Windows,VRPTW）：在VRP的基础上添加了客户点被访问的时间窗约束，除了行驶成本之外, 成本函数还要包括由于早到某个客户点而引起的等待时间和客户需要的服务时间。需求点的时窗限制可以分为两种，一种是硬时间窗（Hard Time Window），硬时窗要求车辆必须要在时窗内到达，早到必须等待，而迟到则拒收；另一种是软时间窗（Soft Time Window），不一定要在时窗内到达，但是在时窗之外到达必须要处罚。
（4）带有随机需求的的VRP （VRP withwith Stochastic Demands,VRPSD）：一个场站（depot）内具有容量的车辆向若干个具有随机需求的客户点服务，客户点的数量及坐标固定，其需求满足一定的随机分布，要求决策派出的车辆数和行驶路线。例如，银行运钞车对城市里银行的取款服务、垃圾回收中心对垃圾点的垃圾回收等均可视作VRPSD[2] 。
（5）带有pickup和delivery的VRP（VRP with Pickup and Delivery,VRPPD）：
客户既有需求（配送）又有供应（集货）的VRP问题，VRPPD的一个重要特征是可以利用回程车辆的空闲容量（提高车辆的装载率）降低成本。问题的目标是在已知客户需求和供应的数量前提下，找到一系列配送路线满足每个客户需求，根据实际情况可以添加车辆容量约束、时间窗限制及里程和时间限制约束等，并实现如路程最短、费用最少和时间尽量少等目标。VRPPD问题总体可分为：先配送后集货[3]、混合式的配送和集货[4]以及同时配送和集货三类。
图1-2 VRP的分类及相互间的关系
据图可知，CVRP是上述所有问题的基础，其重要性不言而喻，学界围绕CVRP开展的研究也最为广泛和深入。因此，本文以CVRP问题为例，并使用爬山算法（HC）进行求解，故有必要在此对CVRP做进一步的阐述与解释。
1.2.1 What is CVRP？
在标准CVRP中，所有客户对应的交付和需求是确定的，所有车辆完全相同且配送任务始于同一个场站，仅对这些车辆引入容量约束。CVRP的目标是最小化车辆行驶路线成本。
1.2.2 CVRP的约束条件
CVRP致力于寻找一组路径，满足：
（1）每条路径的起点和终点均为场站；
（2）每个客户只通过一条路径访问一次；
（3）分配到每条路径的客户总需求不超过Q。
1.3 VRP与TSP的联系与区别
上一篇博客我们介绍了HC求解TSP问题（LocalSearch之爬山算法（1））,这一篇继续探究如何利用HC求解VRP问题，因此有必要弄清VRP与TSP问题的联系和区别。
联系：
（1）TSP是VRP的特例，由于Gaery已经证明TSP是NP难问题，因此VRP也属于NP难问题。
（2）简单的VRP可视作多个TSP的组合，VRP是TSP的拓展。
区别：
TSP在现实物流场景中等价于：1个物流配送公司，只用1辆车，欲将n个客户的订货沿1条最短路线全部送到；而VRP的客户群体庞大，且存在容量、时间、距离等多重限制，只有1辆车或1条路径通常无法满足客户需求，基本上是多辆交通运输工具的行车顺序及路径规划两个问题的求解。相对于TSP，VRP更复杂、求解更困难，但也更接近实际情况。
1.4 VRP问题的常用求解方法
目前，求解VRP的方法非常多，大致可分为精确算法和启发式算法两大类。其中，精确算法主要有分支定界法（Branch and Bound Approach）、割平面法（Cutting Planes Approach）、网络流算法（Network Flow Approach）、动态规划算法（Dynamic Programming Approach）几种。由于精确算法的求解时间将随着客户点数目的增加呈指数增长，因此迄今为止，为解决VRP而设计的最有效的精确算法最多也只能包含50个客户点，其应用范围十分有限。相比之下，能在可接受时间内生成可接受的解的启发式算法更受学界青睐，模拟退火算法（Simulated Annealing）、禁忌搜索算法（Tabu Search）、遗传算法（Genetic Algorithm）、蚁群算法（Ant Colony）等均已被证实可用于求解VRP，是当前研究的热点。
1.5 VRP解的编码形式
编码形式多样，此处介绍两种方式：
（1）基于Permutation的编码：以本实验为例，设城市数目为I，车辆数目为J，VRP的某解可表示为{1,…I+1,…,I+J-1,…I}，大于I的元素视为分隔符/车场。例如：（3,6,4,7,1,8,2,5）说明有6个城市（1、2、3、4、5、6），车1的配送路径为[3,6,4]、车2的配送路径为[1]、车3的配送路径为[2,5]。
（2）两段式编码：通过查阅文献，介绍另一种VRP解的编码形式——两段式染色体编码方案。该编码由前后两段组成，以整数形式进行编码，如图1-3所示。
图1-3 两段式染色体编码方式
第一段编码为无序排列，表示n个客户被访问的顺序；第二段编码总长度为m，表示m台服务车辆在第一段编码中所对应的访问客户点和访问的客户数，并且m个数的总和等于n。如上图，车辆2服务的客户点为第一段中的第5-8位，依序为6,11,3,5，车辆2在第二段编码中所对应的数字4，即表示车辆2在此段编码所要服务的客户数为4。需要注意的是，不同编码方式一般对应不同的LocalSearch算子。
二、用HC求解CVRP问题
2.1 算法设计
本文使用Matlab工具，以CMT1算例为例使用HC求解CVRP问题。Matlab具体现实代码见文末链接。
（1）CVRP解的编码
本文采用整数编码的方式，对各需求点从1开始进行编号。假设需求点数量为 I，车辆数为J ，那么一个可行解的编码长度为 I+J-1，对1到I+J-1的自然数随机排列即可得到一个解，其中大于 的数字充当物流中心用于划分路径。
例如，假设有6个需求点、3辆车，有一个可行解为（3，6，4，7，1，8，2，5），那么其中大于6的数字：7和8，充当物流中心用于划分路径。实际上这个解被划分为以下3条路径：
物流中心 -> 3 -> 6 -> 4 -> 物流中心
物流中心 -> 1 -> 物流中心
物流中心 -> 2 -> 5 -> 物流中心
初始解采用随机生成的方式。且保证可行。
（2）初始数据
给定需求点数量I 、车辆数 J、需求点的需求量、车辆的承载量、各需求点以及物流中心的坐标。通过所有点（包括物流中心）的坐标求出各个点之间的距离，生成距离矩阵。
（3）新解产生
邻域算子：依旧采用基于Permutation的算子：swap/insertion/2-opt，详见blog：LocalSearch之爬山算法（1）3.2节
解的选择策略：first Improve/best Improve
（4）目标函数
本题的目标函数是车辆行驶总距离的最小化。累加路径上的各点距离即可。
（5）约束条件
每条路线的需求量总和不大于车辆的承载量。
（6）停止条件
当前解到达局部最优，算法停止。
（7）惩罚函数
惩罚函数是求解约束优化问题的一种有效方法，其将约束优化问题中的约束违反度乘以惩罚项加到目标函数中，从而构造出带参数的增广目标函数。这种构造的主要思想是把系列的约束优化问题转化为无约束的优化问题进行求解。
本文对车辆的超载量进行惩罚，代码如图2-1：
图2-1 超载惩罚机制
可以看出，objective代表的是目标函数值，penalty Capacity表示的是惩罚系数。若一个解的若干条路线上存在需求量大于车载量的情况，就将所有路线上需求量超出车容量的部分相加乘以惩罚系数，得到惩罚项penalty。最终，一个解的适应度值Cost等于目标函数值加上惩罚项，用于解的评估。
接下来介绍惩罚系数的设定方式：
Vidal[7]文章中惩罚系数（penalty parameters）的初始值设定如下：
其中， 表示AVG-c客户点之间距离的平均值， AVG-q表示各客户点需求量的平均值。
Vidal[7]的惩罚系数是动态变化的，且初始值采用均值的方式计算，为简化计算，我们借鉴Vidal[7]中惩罚系数初始值的计算方式设置静态的惩罚系数，且选取最大值计算。在程序中用penalty Capacity表示，代码实现如图2-2所示：
图2-2 penalty Capacity的代码实现
先求需求点之间距离的最大值与需求量最大值的比值，与1000进行比较，选较小值，再与0.1比较，选较大值。这样的设计下，惩罚系数的范围是0.1到1000，需求点间距离越大、需求量越小，惩罚系数会越大；需求点之间距离越小、需求量越大，惩罚系数就越小。
值得注意的是： 由于本文实现算法不允许非可行解产生，即每条路线的需求量总和不能大于车辆的承载量，所以可行解的惩罚项penalty均为0。对于本实现算法，惩罚函数并没有实质性作用。后续会尝试设计允许不可行解的算法架构，届时惩罚成本将发挥作用。
2.2 解决思路
以first Improve策略为例，伪代码如下：
伪代码具体阐述如下：
步骤1：随机生成一个初始可行解i，
步骤2：从当前解i的邻域Si中选取一个新可行解j
步骤3：如果j的目标函数值小于i的目标函数值，那么j成为当前解
步骤4：如果对于所有的j∈Si，都有f(j)≥f(i)（邻域Si中所有解的目标函数值均不小于当前最优目标函数值，解没有得到改进），爬山结束，当前解i就是局部最优解
步骤5：返回步骤2
2.3 Matlab代码实现
Matlab代码见文末链接。
本文使用CMT1算例，需求点数量为50，车辆数为50，每辆车的承载量均为160。
若邻域选择策略采用best Improve，LS方式采用2opt，最终的实验结果如图2-3所示：
图2-3 实验结果
该结果表示共有6条路线，编号代表不同的城市。以第一行为例，第一条路线的行车轨迹是：物流中心 ->3 -> 20->35 ->36 ->28 ->31 ->22 ->1 ->物流中心，该路线总长为114，路线上各城市需求量总和均没有超出车载量。
在不同条件下，进行多次实验，结果如下：
用爬山算法解决TSP问题详见LocalSearch之爬山算法，使用Berlin52算例的结果如下：
CMT1算例的解决过程中，一个可行解的编码长度为99，而Berlin52算例中，一个可行解的编码长度为52，所以在用邻域算子产生新解的时候，CMT1 算例中解的邻域更大，局部搜索次数更多，耗时会更长。
综合爬山算法解决VRP问题和TSP问题来看，2-opt邻域算子相比于swap和insert更适合这两个问题的求解，2-opt的求解速度比较快，而且解的质量也能得到一定的保证。
三、FAQ
Q1：基于HC求解CVRP、TSP问题的异同？
A1：CVRP和TSP的目标是一致的，即最小化路线成本cost。然而，基于HC求解的CVRP、TSP在解的编码、惩罚值、可行的保证等方面存在着不容忽视的差异。
在TSP问题中，对若干城市的不同访问顺序即为不同解的表示。假设有n个城市需要旅行商访问，那么可以对这些城市从1到n进行编号，不同的排列顺序即为不同的解。例如，有8个城市需要访问，将这些城市从1到8进行编号，对编号数任意排列，即可生成一个解。比如（3，4，6，7，1，8，2，5）就是一个解的表示。而上文1-5中提到，VRP常采用向I个城市的排列顺序中插入J-1个点的方法将原有序列分为J段，每段路径由1辆车负责配送；通过插入J-1个车场标志将VRP问题的解表示成了TSP解的表示形式，使得解决TSP问题的算法能够很好的解决VRP问题，增强了算法应用求解的通用性。
下图3-1、3-2分别是算例CMT1在par_hc.maxIter = 1、par_hc.localsearchType = '2opt'、par_hc.improveType = 'firstImporve'、par_hc.ISDEBUG = 1条件下和算例berlin8在相同条件下基于HC的局部最优解输出结果。
图3-1 CVRP的解的表示
图3-2 TSP的解的表示
可见，TSP的解的表现形式比较简单，只输出了最佳移动路线Position、路程/成本Cost以及最终解和初始解相差的Gap值。CVRP的解则相对复杂，罗列出参与送货的所有车辆各自的配送路径，并计算了每条路径上车辆的配送距离以及完成所有客户点的配送任务后行驶的总距离，同时，利用惩罚值penalty强制确保解的可行，如图3-3所示。
图3-3 penalty的作用
penalty代表的时车辆的超载量，若penalty>0，则判定不可行，在初始解生成和邻域过程中将展开重新搜索，直至满足penalty=0方可继续进行。
Q2：本实验代码编写的不足与改进方向？
A2：不足：虽然计算了penalty，但penalty没有发挥实质的作用。penalty=0时解才可行，初始解和local search中只接受可行解。
改进方向：1）不强制保证解可行，使用penalty进行惩罚，通过在搜索迭代的过程中抛弃差的不可行解；2）若要保证初始解一定可行，除了用random多次随机生成只接受可行解还可以使用修复算子修复random生成的不可行解，保证可行；3)邻域算子的设计也很重要，尝试更多针对VRP问题特征的邻域算子，如Relocate, Swap, K-opt, Cross-Operator, 限于篇幅，此处不予展开，具体可参考Groer [8]论文中的Fig.1。
图3-4 针对VRP问题特征的邻域算子示例
Q3：local search的邻域动作是不是需要避免车车交换？
A3：根据VRP的问题特征决定。同质车辆交换无意义，可将“for jj=1:I”调整为“for jj=1:N”，同时在if循环中增设限制条件：“tour1(ii)>I && tour1(jj)>I”即在访问序列中搜到车就重新找ii、jj，允许城—城交换、城—车交换，杜绝车—车交换，保证搜索质量，如图3-5所示。
图3-5 避免swap车车交换的途径
若车辆异质，则车车交换有意义，不需要避免该情况的发生。
四、小结
本文首先介绍了VRP问题的分类、常见求解方式和编码方式，着重介绍了VRP问题与TSP问题的区别与联系；然后通过Matlab实现了用爬山算法求解CVRP问题，并对比了HC解决TSP问题与CVRP问题的不同；最后指出了本文代码编写存在的不足之处，列举并回答了一些学习过程中可能产生的疑问。
【囿于作者的学术水平，本文的诸多表述或许不够严谨，用途仅限于优化算法入门者的非正式学术交流，不足之处还望各位读者指正。】
参考文献：
[1] Dantzig, G. and Ramser. The Truck Dispatching Problem[J]. Management Science,1959,6:80-91.
[2] 陈宝文,宋申民,陈兴林. 随机需求车辆路径问题及其启发式算法[J].计算机工程与设计,2007,(01):138-141+148.
[3] Wassan NA.A reactive tabu search for the vehicle routingpproblems[J]. J Oper Res Soc, 2006, 57(1) : 111-116.
[4] Crispim, Brandao.Metaheristcs applied to mixed and simulta-neous extensions of vehicle routing problems with backhauls[J].Operation Research, 2005, 56: 1296-1302.
[5] GendreauM, Potvin J Y. Handbook of Metaheuristics[M]. Switzerland: SpringerNature,2019.
[6] 朱智鹏,石宇强,蔡跃坤,韩婷.动态需求下车间生产物流VRP优化[J].西南科技大学学报,2020,35(03):68-74+85.
[7] Vidal T , Crainic T G , Gendreau M , et al. A Hybrid Genetic Algorithm for Multidepot and Periodic Vehicle Routing Problems[J]. Operations Research, 2012, 60(3):611-624.
[8] Groer C , Wasil E , Golden B . A library of local search heuristics for the vehicle routing problem.
相关资料：
作者：宋博文 张雨沁
校对：朱燕燕
