- 1 -
中国科技论文在线
求解 CVRP 问题的快速迭代局部搜索算法#
刘万峰,李霞**
基金项目:国家自然科学基金资助项目(61171124);高等学校博士点基金资助课题 (200805900001)
作者简介:刘万峰(1974-),男,博士研究生,主要研究方向:智能计算
通信联系人:李霞(1968-),女,教授,主要研究方向:智能计算、图像处理. E-
(深圳大学信息工程学院,深圳 518060)
5 摘要:本文提出了一种求解带有容量约束的车辆路径问题(Capacitated VRP,CVRP)的快速
迭代局部搜索算法(Fast iterated local search,FILS)。该算法通过引入“前载重”和“后载重”
的概念,减少了在局部搜索算法中计算相邻解适应度值的复杂度,从而提高算法的运行速度。
此外,算法将派车成本折算为运输成本,使得 CVRP 问题简化为一个单目标优化问题。实
验结果表明,与其他相关算法相比较,该算法能在更短时间内求得满意解,具有很强的实用10
性。
关键词:车辆路径问题;迭代局部搜索;启发式算法
中图分类号:
A fast iterated local searching algorithm for capacitated 15
vehicle routing problem
LIU Wanfeng, LI Xia
(Information Engineering School, Shenzhen University, Shenzhen 518060)
Abstract: This paper proposes a fast iterated local searching algorithm (FILS) to solve the
capacitated vehicle routing problem. The concepts of “pre-load” and “post-load” are defined 20
which can be skillfully used to compute fitness, so that the computation complexity in the local
search algorithm is reduced. Meanwhile, the CVRP can be viewed as a single-objective
optimization problem if we incorporate the "vehicle cost" smartly into the transportation cost.
Experimental results show that the proposed algorithm is feasible with great practicability.
Compared with other algorithms, it obtains better achievement. 25
Key words: vehicle routing problem (VRP); iterated local search (ILS); Heuristic algorithm
0 引言
车辆路径问题(vehicle routing problem ,VRP)于 1959 年由 G. Dantzing 和 J . Ramser[1]提
出。自提出以来,VRP 问题一直是运筹学与组合优化领域的热点问题。VRP 是一 NP 难问30
题[2],因此,设计能较短时间内获得较好解的算法是该问题的主要研究方向。
带有容量约束的车辆路径问题(Capacitated VRP,CVRP)是最基本的 VRP 问题,Potvin
在文献[3]中对用于 CVRP 问题的进化算法进行了综述。近几年, 研究者又提出一些新的算法,
骆等[2]提出一种将幂律极值动力学优化(Power Law Extremal Optimization,τ -EO)融合于混
合蛙跳算法(Shuffled Frog Leaping Algorithm,SFLA)的改进算法实现 CVRP 求解,该算35
法在加强算法的局部搜索能力的同时,保证了算法全局收敛性。陈等[4]提出了基于多邻域的
迭代局部搜索算法(HILS 算法),取得较好的结果。
本文在求解 CVRP 算法中提出“前载重”和“后载重”的概念,使得在局部搜索算法
中计算相邻解的适应度值的复杂度由 O(n)减小为 O(1),从而极大地加快了局部优化算法的
执行效率。通过将派车成本转化为运输成本的方法,将 CVRP 问题的两个优化目标:车辆40
数最少和运输成本最小,统一为一个优化目标,即运输成本最小,简化了 CVRP 求解算法
- 2 -
中国科技论文在线
设计。在本文算法中,利用“前载重”和“后载重”信息对求解 CVRP 问题中的 4 种局部
搜索算子(插入、交换、2-opt、2-opt*)进行改进,极大地提高了 4 种局部搜索算子的运算
速度。最后,将改进后的局部搜索算子和 ILS[4]算法相结合设计了一种快速迭代局部搜索算
法 FILS。通过在 benchmark 问题上与文献[3]算法相比较,该算法能在更短的时间内找到较好45
的解。
1 CVRP 数学模型
CVRP 研究目标是对一系列的客户(城市)设计适当的路线,使车辆有序地通过,在满
足容量的约束条件下,达到使用车辆数最少、运输成本最小的优化目标。
一般情况下,CVRP 问题可以这样描述:设某配送中心有 M 辆车,需要对 K 个客户(城50
市)进行运输配送,每个客户的货物需求量是 gi i=(1,2,…,K),每辆配送车的最大载重量 Q。
设 ci,j 表示客户 i 到客户 j 的运输成本,如时间、路程、花费等。取配送中心编号为 0,各客
户编号为 i=(1,2,…,K),定义变量如下:
否则
配运送的货物由车客户
否则
驶向由车
,0
,1
,0
,1 si
y
jis
x isijs
则 CVRP 的数学模型如下: 55
K
i
K
j
M
s
ijsij xcZ
0 0 0
min (1)
MsQyg
K
i
isi ,...,2,1
0
(2)
0,
,...,2,1,1
1 iM
Ki
y
M
s
is
(3)
MsKjyx js
K
i
ijs ,...,2,1;,...,2,1
0
(4)
MsKiyx is
K
j
ijs ,...,2,1;,...,1,0
1
(5) 60
上述模型中,式(2)为车的容量约束;式(3)保证了每个客户的运输任务仅由 1 辆车完成,
而所有运输任务则由 M 辆车协同完成;式(4)和式(5)限制了到达和离开某一客户的车辆有且
仅有 1 辆。
2 快速迭代局部搜索算法求解 CVRP 问题
CVRP 问题解的构造 65
本文中 CVRP 的解 X 采用可变长编码方式,对于有 K 个客户的 CVRP 问题定义解向量
X 的数据结构为:
1 2 1
0 1
( , , , , )
0,1,2, 2,3, , 1
i
N N
i
x i i N
X x x x x
x K i N
当 或 =
当
(6)
其中 0 可以在编码中出现多次,1,2, … ,K 在编码中出现且仅出现一次。解 X 中相邻的
两个 0 之间的部分表示一条子线路(即某一配送车的配送线路)。例如一个 8 客户的 CVRP70
问题的解 X=(0,3,4, 5,0, 1,8,2,0,6,7,0),表示该配送方案共使用了 3 辆
配送车,3 条子线路分别为:
- 3 -
中国科技论文在线
子线路 1:0,3,4,5,0
子线路 2:0,1,8,2,0
子线路 3:0,6,7,0 75
在此这种编码方式下,编码的长度决定了所使用的配送车辆数目,同时允许在算法运行
过程中对车辆数进行增减。例如在对一个不可行解进行可行化时,只需在其超载子线路的适
当位置插入 0(增加配送车辆使其变为可行解);同样在优化过程中将编码中出现的连 0(表
示该配送车无配送任务)用单个 0 替代可实现车辆数的减少。
快速的局部搜索算法 80
前载重与后载重
在求解优化问题中,计算解的适应度值 f 占用了大量时间。为了减少适应度值的计算时
间,在求解 TSP 问题的局部搜索算法中,通常只计算解 X 和其近邻解 X ′的差异△f =f (X ′)-f
(X),从而将适应度值的计算复杂度由 O(n)减小为 O(1)。然而在 VRP 问题中,由于配送车辆
的容量约束,根据某一局部搜索算法得到的近邻解 X ′ 可能是不可行解,仅仅计算适应度值85
的改变量△f 不能正确反映解 X ′的适应度值。为此本文提出“前载重”和“后载重”的定义。
定义 1:假设 X= ( x1 ,x2 , … xN-1,xN) 为 CVRP 问题的一个可行解,xi 由车辆 v 进行配送,
其所在子线路 Xv=( P1 , … Pj-1 , xi , Pj+1 , … Pn-1 , Pn) (P1=0,Pn=0,P2…Pn-1表示车辆 v 依次服
务的客户,Pj =xi ),则 xi的前载重和后载重定义为:
1
1
j
k
P
pre
x ki
gg (7) 90
n
jk
post
x
kPi
gg
1
(8)
计算可行解 X 中每个元素 xi 的前载重和后载重的目的是将 xi 所在子线路的车辆载重信
息体现到 xi上。
快速局部搜索算法的实现
为了扩大局部搜索范围,本文算法依次使用了 4 种局部搜索算法。在算法中,对可行解95
X=( x1 , … ,xN) 中的所有 xi 搜索最佳位置进行最佳变换,直至不能获得更好的近邻解。采用
的局部搜索(局部变换)算法为:插入、交换、2-opt、2-opt*,分别描述如下:
a) 插入:将解 X 中的某个元素 xi 从当前位置移到另一个元素 xj 之后
b) 交换:将解 X 中的两个元素 xi 和 xj进行位置互换
c) 2-opt:假设 xi 和 xj 是解 X 中互不相邻的两个元素, xi+1和 xj+1分别是 xi和 xj在路100
径中的直接后继节点。 2-opt 是指删去弧(xi , xi+1)、 (xj , xj+1),添加弧(xi , xj)、
(xi+1,xj+1) , 从而得到一条新路径(xi 和 xj 可属于同一子线路, 也可属于不同子线路)。
如图 1 所示。
105
110
a. 同一路径内 2-opt 交换
xi
xi+1
xj
xj+1
xi
xi+1
xj
xj+1
Depot Depot
- 4 -
中国科技论文在线
115
图 1 2-opt 交换示意图
2-opt interchange 120
d) 2-opt*:假设 xi 和 xj 在解 X 中分属两条不同子线路,xi+1和 xj+1分别是 xi和 xj在路
径中的直接后继节点。2-opt*是指删去弧(xi , xi+1)、 (xj , xj+1),添加弧(xi , xj+1)、 (xj ,
xi+1), 从而得到一条新路径。如图 2 所示。
125
图 2 2-opt*交换示意图 130
2-opt* interchange
为了实现快速搜索,在以上 4 种局部搜索算法中采用计算解 X 和其近邻解 X ′的差异△f
的方法来实现对解 X ′适应度值的快速计算。首先根据(7)和(8)式计算可行解 X 中每个元素 xi
的前载重和后载重;其次判断近邻解 X ′是否为可行解,采用 4 种不同局部变换方法(插入、
交换、2-opt、2-opt*)所得近邻解是否可行的判断方法分别如式(9)~(12)所示。 135
QggggQgg postxxx
pre
x
post
x
pre
x jjijii
)(&)(
11
(9)
QgggQggg postxx
pre
x
post
xx
pre
x jijiji
)(&)( (10)
QggggQgggg postxxx
post
x
pre
xxx
pre
x jjiijjii
)(&)(
1111
(11)
QggggQgggg postxxx
pre
x
post
xxx
pre
x iijjjjii
)(&)(
1111
(12)
以上各式等于 1 时表示 X′是可行解,否则为非可行解。如果 X′是可行解,则计算 X ′ 和140
解 X 的差异△f 。4 种不同局部变换方法(插入、交换、2-opt、2-opt*)所得解 X ′的△f 的
计算公式如式(13)~(16)所示。
)()(
111111 ,,,,,,
jjiiiijiijii xxxxxxxxxxxx
ccccccf (13)
)()(
11111111 ,,,,,,,,
jjjjiiiijiijijji xxxxxxxxxxxxxxxx
ccccccccf (14)
)()(
1111 ,,,,
jjiijiji xxxxxxxx
ccccf (15) 145
)()(
1111 ,,,,
jjiiijji xxxxxxxx
ccccf (16)
以上各式中 xi-1 、xi+1 、xj-1 、xj+1 分别为解 X 中元素 xi 和 x 的前驱节点和后继节点。可
以看出 (9)~(16)式的计算复杂度均为 O(1),说明采用以上方法后,计算相邻解 X ′的适应度
值的复杂度由 O(n)降低为 O(1)。
构造 CVRP ′问题 150
CVRP 问题有两个优化目标:车辆数最少和运输成本最小。由于车辆数最少是首要优化
b. 不同路径间 2-opt 交换
Depot
xi
xi+1 xj
xj+1 Depot
xi
xi+1 xj
xj+1
Depot
xi
xi+1
xj
xj+1
xi
xi+1
xj
xj+1 Depot
- 5 -
中国科技论文在线
目标,若能保证运输成本最小的解也是车辆数最少的解就可以将两个优化目标统一为一个优
化目标,即运输成本最小。为此本文构造了 CVRP 问题的偏移问题 CVRP′,在 CVRP′问题
中引入派车成本的概念。派车成本即增加一辆配送车需要付出的成本 C,C 为一常数,算法
中将其折算为运输成本。其折算方法为:将每个客户点到仓库的运输成本增加 C/2 (如图 3155
所示)。在图 3 中形象地将增加的运输成本 C/2 表示为仓库在垂直方向上从 Depot 位置搬到
Depot′(图中 Depot′仅和 Depot 有路径相连),这样任一个客户到新仓库 Depot′的运输成本
比到 Depot 的运输成本增加了 C/2。一辆配送车一来一回增加的运输成本即为 C(等于派车
成本)。
由 CVRP′问题的构造方法易知,CVRP′问题的解 X 的适应度值 f ′(X)=f(X)+MC ( M 为解160
X 使用的车辆数,f ′(X)、f(X)分别表示解 X 在 CVRP′和 CVRP 问题中的适应度值 ),所以解
X 也是原 CVRP 问题的解,其适应度值可直接计算。对 CVRP′问题有如下定理:
定理 1:如果 C>2KCmax ( Cmax=max(ci,j) i, j=0,1,2 … K ),对于 CVRP′问题的任意两个
解 X1 ,X2 (车辆数分别为 K1,K2 ),若 K1 <K2 则 f ′(X1)< f ′(X2)
165
170
图 到 CVRP′转换示意图
CVRP to CVRP′ 175
证明:由 CVRP′问题的构造方法可知:
f ′(X1) = f (X1)+ K1C f ′(X2) = f (X2)+ K2C (17)
0 ≤ f (X1) ≤2KCmax 0 ≤ f (X2) ≤2KCmax (18)
0 ≤︱f (X1)- f (X2)︱≤ 2KCmax (19)
f ′(X1)-f ′(X2) = f (X1)-f (X2) +( K1 -K2) C ≤ 2KCmax+( K1 -K2) C 180
由于 K1 < K2 ( K1,K2 为整数) , C>2KCmax
因此 f ′(x1)-f ′(x2) ≤ 2KCmax-C < 0 证毕。
(18)式中当每辆配送车服务一个客户且每个客户到仓库的距离均等于 Cmax 时取等号。
定理 1 说明如果取 C>2KCmax,在 CVRP′问题中运输成本最小的解 Xmin 同时还是车辆数
最少的解,这样在 CVRP′问题中就可以将两个优化目标简化为一个。由 CVRP 问题和 CVRP′185
问题之间的关系可知:解 Xmin 为原 CVRP 问题中车辆数最少的解,同时其还是所有车辆数
最少解中运输成本最小的解,通过求解 CVRP′问题就能完成在使用车辆数最少的要求下实现
运输成本最小的优化目标。
FILS 算法
ILS(iterated local search, ILS)算法由 Baum[5]在 1986 年提出,Johnson[6]在此基础上设计190
了 ILK(iterated Lin-Kernighan algorithm)算法用于求解 TSP 问题,ILK 算法是目前求解 TSP
问题的最有效的算法之一。本文将求解 CVRP 问题的快速局部搜索算法和 ILS 算法相结合
6 6
1
2 3
4
5
7
8
Depot
Depot’
C/2
1
2 3
4
5
7
8
Depot
- 6 -
中国科技论文在线
设计了快速迭代局部搜索(FILS)算法, FILS 算法的伪代码如图 4 所示。
195
200
205
图 4 FILS算法的伪代码
pseudo code of FILS
在 FILS 算法中解采用的扰动方法为随机交换解 X*中 3 对元素位置;解 Xl的接收准则采
用类似模拟退火[7]的方法,即若 f (Xl )< f ( X* )接受 Xl , 否则产生[0,1]之间的随机数 s,以概210
率 s <exp(( f ( X* )- f ( Xl ))/ f ( X* )/Ti ) 接受 Xl 。
3 实验仿真及分析
本文的实验平台为 Intel E7500 G CPU,操作系统为 Windows XP,编程语言采用
C++,共进行了两组实验。第一组实验用于测试求解 CVRP′问题是否能够实现对车辆数的优
化。实验测试数据来源于 Vehicle Routing Data Sets[8],对于该测试集中多数实例,本文算法215
可以求得和最优解[9]车辆数一样的解,而有些实例求得解的车辆数比最优解的车辆数多。对
那些不能求得和最优解车辆数一样的实例,通过求解它们的偏移问题 CVRP′(取 C>2KCmax)
得到了和最优解车辆数一样的解,实验结果如表 1 所示。
表 1 求解 CVRP′ 问题实验结果
The results of CVRP′ problems 220
实例
最优解 本文算法(求解 CVRP) 本文算法(求解 CVRP′)
车辆数 运输成本 车辆数 运输成本 车辆数 运输成本
A-n61-k9 9 1034 10 1035 9 1035
B-n51-k7 7 1032 8 1016 7 1032
B-n57-k7 7 1153 8 1140 7 1157
P-n22-k8 8 603 9 590 8 603
P-n50-k8 8 631 9 630 8 638
从表 1 可见,通过求解 CVRP′问题能够获得了车辆数更少的解,这说明在求解 CVRP′
问题的算法中,仅对运输成本进行优化就能实现对车辆数的优化,实现了构造 CVRP′问题的
目的。
第二组实验是对算法的整体性能进行测试,在实验中,模拟退火的参数设置为:初始温
度 T0=,每隔 30 代按 Ti+1 = 0. 9×Ti 更新 1 次;算法的停止准则为:当连续 800 代 f ( X
*
)225
的值没有提高,算法结束。实验中采用Christofides, Mingozzi和Toth[10]提出的 14个 benchmark
问题中的 7 个 CVRP 问题进行验证(其它 7 个问题为带有最大行驶时间约束和容量约束的
Begin
初始化:随机生成初始解 X0
进行局部搜索:X*=LocalSearch(X0)
While (停止准则未满足)
对 X*进行扰动:X ′=Perturbation(X*)
进行局部搜索:Xl =LocalSearch(X ′ )
If (Xl 满足接收准则) then
X
*=
Xl
end
end
End
- 7 -
中国科技论文在线
VRP 问题,非基本的 CVRP 问题)。实验结果与文献[4](其实验平台为 Pentium IV ,
采用 C++语言)进行了比较,如表 2 所示。
表 2 算法性能测试结果 230
The performance of algorithm
实
例
客户
数
已知最
优解
文献[3] HILS 算法 本文算法(FILS)
最好解 平均解 标准差 t(s) 最好解 平均解 标准差 t(s)
C1 50 6
C2 75
C3 100
C4 150
C5 199
C11 120
C12 100 0 0
表 2 中最好解、平均解、标准差是独立运行 20 次,记录每个实例的最好解,同时计算
得到平均解和标准差。程序运行时间 t 为算法运行一次的平均时间。从表 2 中可见本文算法
在 6 个实例中所取得了比文献[4]更好或相同的结果,同时该算法运行速度提高了 17~37 倍。
例如,对于规模大小为 199 的实例 C5, 本文算法能在 秒左右的时间内求得和最优解相235
差在 %左右的解,这表明本文算法能在更短时间内求得较好结果,具有很强的实用性。
4 结论
在局部搜索算法中,计算相邻解的适应度值会占用大量的 CPU 时间。减少相邻解适应
度值的计算复杂度是设计高效局部搜索算法的关键因素。本文通过提出“前载重”和“后载
重”的概念,降低了在局部搜索算法中计算相邻解适应度值的复杂度,从而提高了算法的执240
行速度。进一步引入“派车成本”,使得决策者有一个量化指标去均衡车辆数和运输成本两
个优化目标。具体算法实现中将“派车成本”转化为运输成本并构建 CVRP′ 问题,使得车
辆数和运输成本统一为一个优化目标,简化了算法的设计。实验结果表明,该算法可同时对
车辆数和运输成本进行优化,并能够在较短时间内求得较大规模 CVRP 问题的满意解,具
有很强的使用价值。 245
[参考文献] (References)
[1] G. B. Dantzig and J. H. Ramser,The truck dispatching problem[J].Management science,1959, 6(1):80-91.
[2] 骆剑平, 李霞,陈泯融. 基于改进混合蛙跳算法的 CVRP 求解[J].电子与信息学报, 2011,33(2):429-434.
[3] Potvin Jean-Yves. State-of-the art review evolutionary algorithms for vehicle routing[J].Informs Journal On 250
Computing, 2009,21(4):518-548.
[4] 陈萍,黄厚宽,董兴业.基于多邻域的车辆路径优化迭代局部搜索算法[J]. 北京交通大学学报:自然科学
版,2009,33(2):1-5.
[5] , J. K. Lenstra. Local search in combinatorial optimization[M]. Chichester:John Wiley & Sons,1997.
[6] E. B. practical 'neural'computation for combinatorial optimization problems [A]. . 255
Neural Networks for Computing [C].AIP conference proceedings,1986,151-53.
[7] S. Kirkpatrick, C. D. Gelatt, M. P. Vecchi. Optimization by simulated annealing[J].science, 1983, 220: 671.
[8]
[9] R. Baldacci, N. Christofides, A. exact algorithm for the vehicle routing problem based on the set
partitioning formulation with additional cuts[J]. Mathematical Programming, ,2008,115:351-385. 260
[10]