- 1 -
基于贪心环路的减量迭代 TSP优化新算法
李洪波 1,马文军 1,陈军 2
1.鲁东大学管理学院,山东烟台(264025)
2.鲁东大学数学与信息学院,山东烟台(264025 )
摘 要:借用 Kruskal最小生成树算法的思想,按照最短边优先的次序,首先依次生成 TSP
问题环路中的各条边,构造出贪心环路初始解。然后,通过点减量迭代的方式对贪心环路进
行第一次优化。接着,采用单边减量迭代的方式进行第二次优化。最后,采用多边减量迭代
的方式进行第三次优化。为便于操作,将每个点的状态描述成六元组类型,将边描述成三元
组类型,将 TSPLIB中的坐标数据和邻接矩阵数据改造成由三元组数据元素组成的上三角邻
接矩阵。通过对称 TSPLIB中 12种实例和 CHN144总计 13种实例数据进行检测,能把目前
已知最好的解中 eil101、eil76和 brazil58进一步削减优化,将 eil101由 629下降到 620,将
eil76由 538下降到 533,将 brazil58由 25395下降到 18456,并给出了相应解的环路路径解
序列。减量迭代 TSP优化算法的时间复杂度为 O( , 2 3
,
max(( ) , )i j e k
i j e k
sum sum sum n n
n
− − −+ +∆ ∆ ∆
)。
关键词:TSP贪心环路,减量迭代,点减量,单边减量,多边减量
1. 引言
TSP问题是指给定一个城市集合和城市间的距离,从城市集合中的任一城市出发,经过所
有其他城市一次恰好一次又回到出发点的最短回路, 即给定图G= (V , E ,W ) , 其中V 为
顶点集合, |V|= n, E为边集合,W 为边权函数, 求集合{V1, V2, …, Vn}的一个排列P,
使
1
P( i) P( i+ 1) P( n) P(1)
1
W (v , v ) W (v , v )
n
i
−
=
+∑
最小。 该问题是一个经典的组合优化问题, 在VL S I 芯片设计和车辆选路等许多领域有广
泛的应用。
TSP 问题是经典的N P-Hard 问题[1 ] ,在多项式时间内找不到一个算法能保证得到最优
解。因此,人们采用近似算法找近似解, 以在多项式时间内能结束。 TSP 问题的近似算法
分为两类[ 2 ]: 环路构造算法和环路改进算法。前者从某个初始解开始, 通过某种增广策略
逐步求解, 直到得到一个合法环路解为止,这类算法包括最近邻算法、贪心算法、Clarke2 W
righ t 算法[ 3 ] 和Ch ristof ides 算法[ 4 ] 等。环路改进算法则在给定初始的合法解后使用某种策
略来改进初始解。这些策略包括局部搜索、模拟退火[ 5 ]、遗传算法[ 6]等, 其中最简单有效
的方法为局部搜索, 如2-O PT [7 ] , 3-O PT [ 8 ] ,L K[ 9]和循环L K[ 10 ]等. 通常这两类算法
被分开使用, 用环路构造算法来构造初始解,而用环路改进算法改进这个初始解。
本文首先借用Kruskal构造最小生成树的思想构造贪心TSP问题环路初始解,然后,提出
减量迭代新算法,即初始环路被分别用点减量迭代算法、边减量迭代算法和多边减量迭代算
法三种减量迭代新算法优化,在用均匀TSPLIB的部分实例检测后,得到了一些更好的结果。
2. 数据表示和预处理
上三角三元组邻接矩阵和边的描述
为在 O(1)的时间内找到任意两点间的距离,也为了按边长进行排序后的邻接矩阵不丢
- 2 -
失邻接关系信息,不用行号和列号来对应邻接点信息,将对称 TSPLIB 中的邻接矩阵数据和
坐标数据统一转化成三元组邻接矩阵。
定义 1 设 G=<V, E ,W >为无向网,其中 V={ a1, a2, …, an },|V|= n, E 为
边集合,W 为边权函数,W(ai,aj)返回 ai 和 aj 两点间的距离。定义一个三元组邻接矩阵 M,
其中各元素 mi,j 为
, 1 1 , 1 , 1 i( , , ) ( =W(a ,a ))i j n j i n j i n j i n jm a a L L− − − − − − − −=
邻接矩阵的数据元素描述为 1 1 ,( , , )n j i n j ia a L− − − − ,其中 an-1-j 和 ai互为邻接点,而 Ln-1-j,i
为两点间 an-1-j 和 ai 两点间的距离。因为邻接矩阵对称,进一步压缩成上三角三元组邻接矩
阵,如图 1 所示。
0 1 …… n-2
0 (an-1, a0, L n-1,0) (an-2, a0, L n-2,0) …… (a1, a0, L 1,0)
1 (an-1, a1, L n-1,1) (an-2, a1, L n-2,1) …… (a2, a1, L 2,1)
… …… …… ……
n-3 (an-1, an-3, L n-1,n-3) (an-2, an-3, L n-2,n-3)
n-2 (an-1, an-2, L n-1,n-2)
图 1 上三角三元组邻接矩阵示意图
在图 1 中,第 i 列表示 ai到 aj(j=1,2, …, i-1)的距离,我们称之为三元组距离列向
量;而第 j 行表示 ai(i=1,2, …, j-1)到 aj的距离,我们称之为三元组距离行向量。这样
可以避免对角线三元组数据元素(ai, ai, 0)被表示。
邻接矩阵的三元组数据元素实质就是边的描述,因此边也用邻接矩阵的三元组数据元素
来表示。
在实现时,用 C 语言实现,边的描述和上三角邻接矩阵的描述如下:
#define PointNumbers 76
#define MAX (PointNumbers*(PointNumbers-1))/2
typedef struct Edge
{
long length;
long p1;
long p2;
}TspEdge;
TspEdge TSP[MAX]; //存储上三角三元组邻接矩阵,排序操作针对此表进行。
TspEdge TSPOld[MAX]; //存储上三角三元组邻接矩阵,不在此表上进行排序操作。
例如 bayg29 的部分数据如下所示。
TSPOld [406]={
{28,0,97}, {27,0,205}, {26,0,139}, {25,0,86}, {24,0,60}, {23,0,220}, {22,
0,65}, {21,0,111}, {20,0,115}, {19,0,227}, {18,0,95}, {17,0,82}, {16,0,
225}, {15,0,168}, {14,0,103}, {13,0,266}, {12,0,205}, {11,0,149}, {10,0,120},
{9,0,58}, {8,0,257}, {7,0,152}, {6,0,52}, {5,0,180}, {4,0,136}, {3,0,
82}, {2,0,34}, {1,0,145}, {28,1,129}, {27,1,103}, {26,1,71}, {25,1,105}, {24,
1,258}, {23,1,154}, {22,1,112}, {21,1,65}, {20,1,204}, {19,1,150}, { 18,1,
- 3 -
87}, {17,1,176}, {16,1,137}, {15,1,142}, {14,1,204}, {13,1,148}, {12,1,148},
{11,1,49}, {10,1,41}, {9,1,211}, {8,1,226}, {7,1,116}, {6,1,197}, {5,1,
89}, {4,1,153}, {3,1,124}, {2,1,74}, {28,2,219}, {27,2,125}, {26,2,175}, {25,
2,386}, ……, {4,2,220}, {3,2,60}, ……};
点的描述
在本优化算法中,从点是否加入 TSP 解的环路中以及加入后点的状态描述成一个五元
组。
定义 2 设 G=<V, E>为无向网,其中 V={ a1, a2, …, an },对于任意点 ai定义成
一个五元组:
ai=(PointStatus1i, Adj_Point0 i, Adj_Point1 i, Adj_Edge0 i, Adj_Edge1 i, Degree i)。
其中,PointStatus1 i取值为 NoProceed、EndPoint 和 MiddlePoint,分别表示未处理、端
点和中间点三种情况,未处理表示该点还未处理到,端点表示处理到且为 TSP 解的路径的
一个端点,中间点表示处理到了且是 TSP 解的路径的一个中间点,如图 2 所示。
图 2 点状态示意图
在图 2 中,表示边(a1,a2)、(a3,a4)和(a4,a5)是 TSP 解的回路中的边,因此 a1、a2、a3
和 a5点的状态为端点, a4为中间点,如果存在其它点则均为未处理到的点。伴随着点和边
不断加入 TSP 目标环路解中,点的状态变化有两种情况,一种是由未处理状态变为端点,
另一种是由端点变为中间点。例如,假定,又有(a2,a6)加入目标环路解中,则 a2 由端点变
为中间点,而 a6由未处理到的点变为端点。
Adj_Point0i和Adj_Point1i表示该点已经加入到TSP目标环路解中后对应的两个邻接点,
如果该点为端点,则其邻接点只有一个,为 Adj_Point0i;如果该点为中间点,则其邻接点有
两个,现有的分别为 Adj_Point0 和 Adj_Point1 i。Degree i指明第 i 点的度,分别为 0、1 或 2,
用于控制生成 TSP 环路,求解完毕的最后环路中每个顶点的度均为 2,不能出现度为 3 及以
上的顶点。例如,在图 2 所示的示意图中,端点 a2 只有一个邻接点 a1,因此 a2的 Adj_Point0i
值为 a1,Degreei的值为 1。当(a2,a6)加入目标环路解后,端点 a2的另一个邻接点 a6,此时
a2 的 Adj_Point1i值为 a6,Degree i 的值为 2。
值得一提的是,当 V 集合中的点全部处理到,则所有点的状态均为中间点。且每个点
的度均为 2。
typedef enum PointStatus1 {NoProceed, EndPoint, MiddlePoint};
typedef struct {
PointStatus1 Point_Status;
long Adj_Point[2];//相邻的两个邻接点
long Adj_Edge[2];//相邻的两个邻接边的边长
long Degree;
}PointStatus;
PointStatus point_status[PointNumbers];
long sum=0;//存储已经求得的 TSP 环路周长值。
- 4 -
数据格式转换
为实现数据格式的转换,设计两个函数 CoordToAdjency 和 Adjency,CoordToAdjency
用于实现将TSPLIB中的二维平面坐标信息转化成本算法的上三角三元组邻接矩阵,Adjency
用于实现将 TSPLIB 中的上三角邻接矩阵转化成本算法的上三角三元组邻接矩阵。
3. 基于贪心环路的减量迭代 TSP 优化算法
构造基于最短边优先的贪心环路
按照边长递增的思想从小到大依次生成 TSP 环路解的各条边,直至最后形成一个包含
V 中全部点的环路为止。
设 G=(V, E)为无向网,TM 是 N 上环路解中边的集合,同时设 U 为已经加入 TM 中
边依赖的点的集合。算法从 U={},TM ={}开始,重复进行如下操作:在所有的 u∈U,v∈V-U
的边(u,v) ∈E 中找一条代价最小的边(u0, v0),如果并入 TM 集合后能保证 U 集合中每一
个点的度至多为 2,并且保证不生成小于顶点个数| V |的小环路,则将其并入 TM 集合,同
时{u0, v0}并入 U,直到 U=V 并且 U 集合中每一个点的度均为 2 为止。此时,TM 必是必
有包含 V 中所有点的 n 条边,且 U 中每个顶点的度均为 2。
借用 Kruskal 生成最小生成树算法的思想,假设连通网 G=(V, E),令 TSP 环路解的初
始状态为 T=(V, {}),图中每个顶点自成一个连通分量。在 E 中选择代价最小的边,若该
边依附的顶点落在 T 中两个不同的连通分量上,或者 T 已经是连接 V 中全部顶点的一个连
通分量,并且保证并入该边后 V 集合中每一个点的度至多为 2,则将此边加入 T 中,否则
舍弃此边而选择下一条代价最小的边。依次类推,直至 T 中所有顶点的度均为 2 停止。
例如,图 3 表示基于最短边优先的贪心环路 TSP 算法的求解过程。
- 5 -
图 3-b 是每个顶点自成一个连通分量的状态。图 3-c 是加入第一条最短边(A3,A4)后的
状态。图 3-d 是加入第二条最短边(A5,A4)后的状态。在图 3-e 中加入的是权值为 8 的边(A3,
A2),而不是值为 7 的边(A1,A4),因为如果加入(A1,A4),则 A4 顶点的度为 3,故跳过。
在图 3-f 中边(A3,A1)和(A3,A5)的 A3 与 A4相同,故跳过。而边(A2,A5)未加入的原因是
提前生成了不包含 A1点的只有其他 4 个顶点的环路,故也跳过,因此加入了权值为 15 的边
(A5,A1)。最后随着(A1,A2)的加入,形成了所有顶点的度均为 2 的环路,此时求得的 TSP
问题环路序列为 A3A4A5A1 A2A3,周长 49。
为实现按边长从小到大依次处理,对三元组邻接矩阵 M 按照边长进行升序排序,排序
后的三元组邻接矩阵为 M’, 在排序后的 M’上从小到大依次处理。假定采用堆排序的方法
进行排序,则其时间复杂度为 2 logn n 。为防止提前生成环路,引入一个邻接矩阵
PointRelation,其元素 Pi,j 值为 1 时指明 Vi和 Vj 处于一个连通分量上,为 0 时表明不在一个
连通分量上,当新插入环路中的边的连接的两个点如果是一个连通分量中的两个端点且该连
通分量不包含 V 中全部顶点则舍弃该边,继续处理下一条边。当加入一条边后,修改点的
信息、PointRelation 邻接矩阵信息。
算法的非形式化描述如下:
⑴ 对上三角三元组邻接矩阵 TSP 依据边的权值从小大大排序。将所有点对的关系矩阵
PointRelation 的各个元素清 0。将所有点的状态置为 NoProceed。处理过的点数变量 j 清 0。
当前处理的边的数组下标 i 清 0。
⑵ 取上三角三元组邻接矩阵 TSP 的第 i 条边,如果该边连接的点符合下述情况之一,
则允许插入,否则不能插入。
情况①: 该边连接的两个点的状态值为均为 NoProceed。
情况②: 该边连接的两个点的状态值一个为 NoProceed,另一个为 EndPoint,或者一
个为 EndPoint,另一个为 NoProceed;
情况③: 该边连接的两个点的状态值两个点均为 EndPoint,但两个点在两个不同的连
- 6 -
通分量中,或者在一个连通分量中但包含了全部顶点。
如果允许插入,则将该边长值累加到 sum 中,修改该边连接的两点信息和两点
PointRelation 连通分量信息,将该边添加到 ProceedTsp[j]中,j++。
⑶ i=i+1,如果 i 小于 MAX,则转向执行⑵,否则结束。
算法 1 的时间复杂度为 O(n3)。因为上三角三元组邻接矩阵 TSP 数组的数据元素个数为
n(n-1)/2,最坏情况下所有边均被处理到。如果允许插入,修改一个连通分量的两个端点的
时间复杂度为 O(n),原因如下:
在情况①中,直接将两个状态为 NoProceed 的点关系修改为一个连通分量的两个端点,
由下面两条语句实现:
PointRelation[TSP[i].p1][TSP[i].p2]=1; //一个连通分量的两个端点
PointRelation[TSP[i].p2][TSP[i].p1]=1; //一个连通分量的两个端点
在情况②中,首先找到该边的状态为 EndPoint 点对应连通分量的另一个端点,然后修
改另一个端点与该边的状态为 NoProceed 点的 PointRelation 关系信息,由下面语句实现:
for(int a=0; a<PointNumbers; a++)
if(PointRelation[TSP[i].p1][a]==1){
PointRelation[a][TSP[i].p2]=1, PointRelation[TSP[i].p2][a]=1;
PointRelation[a][TSP[i].p1]=0, PointRelation[TSP[i].p1][a]=0;
break;
}
其中,TSP[i].p2 点的状态为 EndPoint,而 TSP[i].p1 点的状态为 NoProceed。情况②中
的另一种情况与此相似,只不过是 TSP[i].p2 和 TSP[i].p1 对调。
在情况③中,因为 TSP[i].p2 和 TSP[i].p1 分别是两个连通分量的一个,查找到该两个连
通分量的各自的另一个端点,将两个连通分量的各自的另一个端点的 PointRelation 修改,是
两个连通分量成为一个连通分量,由下面语句实现:
if( !PointRelation[TSP[i].p1][TSP[i].p2] ){
for(int a=0; a<PointNumbers; a++){ if(PointRelation[TSP[i].p1][a]==1) break; }
for(int b=0; b<PointNumbers; b++){ if(PointRelation[TSP[i].p2][b]==1) break; }
PointRelation[b][TSP[i].p2]=0; PointRelation[TSP[i].p2][b]=0;
PointRelation[a][TSP[i].p1]=0; PointRelation[TSP[i].p1][a]=0;
PointRelation[a][b]=1; PointRelation[b][a]=1;
}
算法的 C 语言实现如算法 1 所示,在算法 1 中,用 ProceedTsp 向量依次存放 TSP 环路
中的各条边。
TspEdge ProceedTsp[PointNumbers];//按照权值递增的次序存储 TSP 环路中的各条边
void TspShortestFirst( )
{
int i, j=0,middle_n=0,kkk1,kkk2;
sort(TSP, 0, MAX - 1);//对上三角三元组邻接矩阵进行排序
for(i=0; i<PointNumbers; i++)
for(j=0; j<PointNumbers; j++) PointRelation[i][j]=0;
for(j=0; j<PointNumbers; j++) point_status[j].Point_Status=NoProceed;
j=0;
for(i=0; i<MAX ; i++) {
kkk1=point_status[TSP[i].p1].Degree;
kkk2=point_status[TSP[i].p2].Degree;
if(point_status[TSP[i].p1].Point_Status==NoProceed && point_status[TSP[i].p2].Point_Status==NoProceed){
sum+=TSP[i].length, j++;
point_status[TSP[i].p1].Point_Status=EndPoint; point_status[TSP[i].p2].Point_Status=EndPoint;
point_status[TSP[i].p1].Adj_Point[kkk1]=TSP[i].p2; point_status[TSP[i].p1].Adj_Edge[kkk1]=TSP[i].length;
point_status[TSP[i].p2].Adj_Point[kkk2]=TSP[i].p1; point_status[TSP[i].p2].Adj_Edge[kkk2]=TSP[i].length;
- 7 -
point_status[TSP[i].p2].Degree++; point_status[TSP[i].p1].Degree++;
ProceedTsp[j-1]=TSP[i];
PointRelation[TSP[i].p1][TSP[i].p2]=1; PointRelation[TSP[i].p2][TSP[i].p1]=1; //一个连通分量的两个端点
}
else if(point_status[TSP[i].p1].Point_Status==EndPoint && point_status[TSP[i].p2].Point_Status==NoProceed){
sum+=TSP[i].length, j++;
point_status[TSP[i].p1].Point_Status=MiddlePoint; point_status[TSP[i].p2].Point_Status=EndPoint;
point_status[TSP[i].p1].Adj_Point[kkk1]=TSP[i].p2; point_status[TSP[i].p1].Adj_Edge[kkk1]=TSP[i].length;
point_status[TSP[i].p2].Adj_Point[kkk2]=TSP[i].p1; point_status[TSP[i].p2].Adj_Edge[kkk2]=TSP[i].length;
point_status[TSP[i].p2].Degree++; point_status[TSP[i].p1].Degree++;
middle_n+=1;ProceedTsp[j-1]=TSP[i];
for(int a=0; a<PointNumbers; a++)
if(PointRelation[TSP[i].p1][a]==1){
PointRelation[a][TSP[i].p2]=1, PointRelation[TSP[i].p2][a]=1;
PointRelation[a][TSP[i].p1]=0, PointRelation[TSP[i].p1][a]=0;
break;
}
}
else if(point_status[TSP[i].p2].Point_Status==EndPoint && point_status[TSP[i].p1].Point_Status==NoProceed){
sum+=TSP[i].length, j++;
point_status[TSP[i].p2].Point_Status=MiddlePoint; point_status[TSP[i].p1].Point_Status=EndPoint;
point_status[TSP[i].p1].Adj_Point[kkk1]=TSP[i].p2; point_status[TSP[i].p1].Adj_Edge[kkk1]=TSP[i].length;
point_status[TSP[i].p2].Adj_Point[kkk2]=TSP[i].p1; point_status[TSP[i].p2].Adj_Edge[kkk2]=TSP[i].length;
point_status[TSP[i].p2].Degree++; point_status[TSP[i].p1].Degree++;
middle_n+=1;ProceedTsp[j-1]=TSP[i];
for(int a=0; a<PointNumbers; a++)
if(PointRelation[TSP[i].p2][a]==1){
PointRelation[a][TSP[i].p1]=1, PointRelation[TSP[i].p1][a]=1;
PointRelation[TSP[i].p2][a]=0, PointRelation[a][TSP[i].p2]=0;
break;
}
}
else if(point_status[TSP[i].p2].Point_Status==EndPoint && point_status[TSP[i].p1].Point_Status==EndPoint)
if( !PointRelation[TSP[i].p1][TSP[i].p2] || middle_n==PointNumbers-2 /*保证一个环路*/ ){
sum+=TSP[i].length, j++;
point_status[TSP[i].p2].Point_Status=MiddlePoint; point_status[TSP[i].p1].Point_Status=MiddlePoint;
point_status[TSP[i].p1].Adj_Point[kkk1]=TSP[i]. p2;
point_status[TSP[i].p1].Adj_Edge[kkk1]=TSP[i].length;
point_status[TSP[i].p2].Adj_Point[kkk2]=TSP[i].p1;
point_status[TSP[i].p2].Adj_Edge[kkk2]=TSP[i].length;
point_status[TSP[i].p2].Degree++; point_status[TSP[i].p1].Degree++;
ProceedTsp[j-1]=TSP[i]; middle_n+=2;
if( !PointRelation[TSP[i].p1][TSP[i].p2] ){
for(int a=0; a<PointNumbers; a++){ if(PointRelation[TSP[i].p1][a]==1) break; }
for(int b=0; b<PointNumbers; b++){ if(PointRelation[TSP[i].p2][b]==1) break;}
PointRelation[b][TSP[i].p2]=0; PointRelation[TSP[i].p2][b]=0;
PointRelation[a][TSP[i].p1]=0; PointRelation[TSP[i].p1][a]=0;
PointRelation[a][b]=1; PointRelation[b][a]=1;
}
}
}
}
算法 1 基于最短边优先的贪心环路算法
在算法 1 中,TSP 环路解信息以两种形式存储,一种形式为 ProceedTsp 边信息表,以
边的形式记录;另一种是 point_status 点信息表,记录点的信息,对于其中的任意元素
point_status[i],记录了顶点 ai在环路中的两个邻接点和两个邻接边权值。例如,在图 3-g 中
的 A1 点信息存储为 point_status[0],其值为{ MiddlePoint, {2,5}, {15,16},2},同理
整个表的最后点信息值为point_status[5]={{MiddlePoint, {2,5}, {15,16},2},{MiddlePoint,
- 8 -
{1,3}, {16,8},2},{MiddlePoint, {2,4}, {8,4},2},{MiddlePoint, {3,5}, {4,
6},2}, { MiddlePoint, {1,4}, {15,6},2}}。ProceedTsp 边信息表便于按边长进行排
序,point_status 点信息表可以在 O(1)时间内找到 TSP 环路解中任意点连接的两个邻接点及
其对应边的权值。sum 存储环路周长。
显然,算法 1 的时间复杂度为 O(n3)。
点减量迭代
点减量迭代算法如算法 2 所示:
算法 2:
⑴ 从初始环路中选择最长的边(ai,aj), ksum =0。
⑵ 从 V-{ai,aj}中选择 ak,假设 ak 在 TSP 环路中的两个邻接点分别为 ak1 和 ak2,函数
W((ai,aj))用于求边(ai,aj)的权值,选中的 ak 须满足公式 1 所示的关系:
公式 1:
1 2 1 2max{ (( , )) (( , )) (( , )) (( , )) (( , )) (( , ))}
, { , }
k i j k k k k k k i k k jk
i j
W a a W a a W a a W a a W a a W a a
k V a a
= + + − − −
∈ −
�
我们把 k� 称之为点 ak 的点减量。
⑶ 如果 0k >� ,则从 TSP 环路中删除边(ai,aj)、(ak1,ak)和(ak,ak2),添加(ak1,ak2)、
(ai,ak) 和(ak,aj),也就是用(ak1,ak2)、(ai,ak) 和(ak,aj)三条边置换(ai,aj)、(ak1,ak)和(ak,
ak2)三条边使得替换后的 TSP 环路周长最小,sum=sum- k� , ksum = ksum + k�,。转向步骤⑴
进行下一次迭代。
否则,如果 k� ≤0,则:
(3-1) 如果存在下一条次长的边,则选择环路中除(ai,aj)外的新的最长边(ai1,aj1),ai= ai1,
aj=aj1,转向步骤⑵。
(3-2) 否则,不存在下一条次长边,即所有的边已处理完毕,结束本算法。
在实现时,利用点信息表 ProceedTsp 能够迅速找到任意点 ak 对应的两个邻接点 ak1 和
ak2 及其对应(ak1,ak)和(ak,ak2)的边长信息。问题是如何在 O(1)时间内找到(ak1,ak2)、(ai,
ak) 和(ak,aj)三条边的边长或权值信息。因为 TSP 上三角邻接矩阵表已按边长进行排序,无
法在 O(1)时间内定位任意两点间的权值,所以我们采取以空间换时间的策略,定义一个
TSPOld,用于存储未排序的上三角三元组邻接矩阵,这样,对于任意两点(am+1,an+1)(用 C
语言实现时,对应的下标分别为 m 和 n),则其在 TSPOld 表中对应的元素位置关系如下:
(2 int 1 ) / 2 int 1 ,
(2 int 1 ) / 2 int 1 ,
m Po erNumber m Po erNumber n m n
k
n Po erNumber n Po erNumber m m n
− − + − − ≤⎧ ⎫= ⎨ ⎬− − + − − >⎩ ⎭
假设算法平均的点减量为 k
−∆ ,则算法的时间复杂度为 O( 2max( , )k
k
sum n n−∆
)。因为一次
迭代时在选边的时间复杂度为 O(n),算法的迭代次数最大为 k
k
sum
−∆
+n,因此,算法的时间
复杂度为 O(( k
k
sum
−∆
+n)n),也即 2max( , )k
k
sum n n−∆
。
- 9 -
算法的 C 语言实现如下所示。
void TspProceedPoint( )
{
long i, j=0,l1,l11=0,l2,l22=0,lastedge;
long k1,k2,k3,add_edge[3], delete_edge,add, row1,row2,row3, column1, column2,column3;
long t,mini,lastmini=PointNumbers+1,maxt,c,oldmaxt,max=0,add0,add1,add2,oldk[3],lastj,lastpair[3][2];
sort(ProceedTsp,0,PointNumbers-1);
for(j=PointNumbers-1;j>=0; ) {
max=ProceedTsp[j].length, l1 = ProceedTsp[j].p1, l2 = ProceedTsp[j].p2, lastedge = ProceedTsp[j].length;
maxt=0;
for(i=0; i<PointNumbers ; i++){
if (i==l1 || i==l2) continue;
delete_edge = point_status[i].Adj_Edge[0] + point_status[i].Adj_Edge[1] + lastedge;
row1 = point_status[i].Adj_Point[0], column1 = point_status[i].Adj_Point[1];
if(row1>column1) {int tt; tt=row1, row1=column1, column1=tt; }
column1=PointNumbers-1-column1;
k1 = (row1*(2*PointNumbers-1-row1))/2+column1; add_edge[0] = TSPOld[k1].length;
row2 = i, column2 = l1;
if(row2>column2) {int tt; tt=row2, row2=column2, column2=tt; }
column2=PointNumbers-1-column2;
k2 = (row2*(2*PointNumbers-1-row2))/2+column2, add_edge[1]=TSPOld[k2].length;
row2 = i, column3 = l2;
if(row2>column3) {int tt; tt=row2, row2=column3, column3=tt; }
column3=PointNumbers-1-column3;
k3 = (row2*(2*PointNumbers-1-row2))/2+column3, add_edge[2]=TSPOld[k3].length;
add=add_edge[0]+add_edge[1]+add_edge[2]; t=delete_edge-add;
if(t > maxt) {
mini=i,oldmaxt=maxt,maxt=t,add0=add_edge[0],add1=add_edge[1],add2=add_edge[2];
oldk[0]=k1, oldk[1]=k2, oldk[2]=k3;
lastpair[0][0]=point_status[i].Adj_Point[0], lastpair[0][1]= point_status[i].Adj_Point[1];
lastpair[1][0]=i, lastpair[1][1]=l1; lastpair[2][0]=i, lastpair[2][1]=l2;
}
}
if(maxt>0 && mini!=lastmini && l1!=l11 && l2!=l22) {
sum=sum-maxt;
int mini_left=point_status[mini].Adj_Point[0]; int mini_right=point_status[mini].Adj_Point[1];
if(point_status[mini_left].Adj_Point[0]==mini)
point_status[mini_left].Adj_Point[0]= mini_right, point_status[mini_left].Adj_Edge[0]=add0;
else point_status[mini_left].Adj_Point[1]= mini_right, point_status[mini_left].Adj_Edge[1]=add0;
if(point_status[mini_right].Adj_Point[0]==mini)
point_status[mini_right].Adj_Point[0]= mini_left, point_status[mini_right].Adj_Edge[0]=add0;
else point_status[mini_right].Adj_Point[1]= mini_left, point_status[mini_right].Adj_Edge[1]=add0;
point_status[mini].Adj_Point[0]=l1, point_status[mini].Adj_Edge[0]=add1;
point_status[mini].Adj_Point[1]=l2, point_status[mini].Adj_Edge[1]=add2;
if(point_status[l1].Adj_Point[0]==l2)
point_status[l1].Adj_Point[0]=mini,point_status[l1].Adj_Edge[0]=add1;
else point_status[l1].Adj_Point[1]=mini,point_status[l1].Adj_Edge[1]=add1;
if(point_status[l2].Adj_Point[0]==l1)
point_status[l2].Adj_Point[0]=mini,point_status[l2].Adj_Edge[0]=add2;
else point_status[l2].Adj_Point[1]=mini,point_status[l2].Adj_Edge[1]=add2;
for(c=0;c<PointNumbers;c++)//改变边的信息,置换三条边
if(ProceedTsp[c].p1==mini && ProceedTsp[c].p2==mini_left || ProceedTsp[c].p2==mini &&
ProceedTsp[c].p1==mini_left )
ProceedTsp[c].p1=mini, ProceedTsp[c].p2=l1, ProceedTsp[c].length=add1;
else if(ProceedTsp[c].p1==mini && ProceedTsp[c].p2==mini_right || ProceedTsp[c].p2==mini &&
ProceedTsp[c].p1==mini_right )
ProceedTsp[c].p1=mini, ProceedTsp[c].p2=l2, ProceedTsp[c].length=add2;
else if(ProceedTsp[c].p1==l1 && ProceedTsp[c].p2==l2 || ProceedTsp[c].p2==l1 &&
ProceedTsp[c].p1==l2 )
- 10 -
ProceedTsp[c].p1=mini_left, ProceedTsp[c].p2=mini_right, ProceedTsp[c].length=add0;
sort(ProceedTsp,0,PointNumbers-1);//选
lastj=j; j=PointNumbers-1;
l11=l1, l22=l2, lastmini=mini;
}
else
j--;
}
}
算法 2-1 点减量迭代的 C 程序
值得注意的是,为简化程序的实现,本算法并没有直接求 TSP 环路解中的最长边,而
是采用了一个排序算法,复杂度比算法 2 略高,但不影响最后的分析和结果。
单边减量迭代
单边减量迭代的思想如算法 3 所示。
算法 3:
⑴ 把经过点减量迭代的 TSP 环路 ProceedTsp 按照边长升序排序,令 i=0, esum =0,即
i 指向未处理过的当前最短边。
⑵ 假设 i 指向的最短边为(al,ar),在环路中找到 al和 ar 的另外两个邻接点 alb 和 arb(lb≠l
且 lb≠r, rb≠l 且 rb≠r),即边(al,alb) 和(ar,arb)为 (al,ar)的两个邻接边。在环路中选择除(al,
ar)、(al,alb) 和(ar,arb)的边(al1,al2),用最短边(al,ar)把边(al1,al2)试探性断开,用边(alb,
arb)连接点 alb 和 arb,在边对(al,al1)和(ar,al2)与边对(al,al2)和(ar,al1)中选择边对之和最小
的边对将边(al,ar)连接环路中,重新维持试探性环路,计算边减量迭代的减量,边减量的定
义如公式 2 所示。
公式 2:
1 2 1 2 1 2 1 2max{ (( , )) (( , )) (( , )) (( , )) (( , )) min( ( , ) ( , ), ( , ) ( , ))}
,
e lb l rb r l l k k lb rb l l r l l r l le
W a a W a a W a a W a a W a a Wa a Wa a Wa a Wa a
e E
= + + − − − + +
∈
�
我们把 e�称为边减量,把 1 2( , ), ( , ), ( , )lb l rb r l la a a a a a 称为待删除边,把
1 2( , ), ( , ), ( , )lb rb l l r la a a a a a 或
, , , ,,
0 max( ) , { | 0 int 1}i j i j i j i ji jdel add del add i j x x Po erNumber− > ∩ − ∩ ∈ ≤ ≤ −
1 2( , ), ( , ), ( , )lb rb l r l la a a a a a 三条边称为待插入边。
⑶ 如果 0e>� ,删除全部待删除边,插入三条待插入边,sum=sum- e�,sumε = sumε + e�,
把 ProceedTsp 重新排序,i=0,转向⑵。否则,i++,如果 i<PoinerNumber,转向⑵,否则结
束。
单边减量迭代的程序实现与点减量迭代的程序实现除了减量的计算方式相关程序段不
同外,其他基本相同,不再累述。
假设算法平均的点减量为 e
−∆ ,如果排序算法采用堆排序,对排序的时间复杂度为
O(nlogn),则算法的时间复杂度为 O( 3max( log , )e
e
sum n n n−∆
)。因为一次迭代时在第⑵步的
语句执行频度为 1,在第⑶用 TSP 环路中的边进行堆排序,其时间复杂度为 O(nlogn),算法
- 11 -
的迭代次数最大为 e
e
sum
−∆
+n,因此,算法的时间复杂度为 O(( e
e
sum
−∆
+n)n2),也即
2 3max( , )e
e
sum n n−∆
。
多边减量迭代
假设 TSP 环路路径序列为 ax1ax2… axm ax(m+1)…axn,多边减量迭代的思想是在 TM 中选择
当前最大的边(axm,ax(m+1)),然后形成 TSP 路径中至少由三点构成的但不包含 ax(m-1) 、axm 、
ax(m+1) 和 ax(m+2) 4 点的所有路径序列 axiax(i+1)…axj,选择这样的序列,将该序列插入到最大边
中形成新的环路路径序列的周长小于插入之前的路径序列且新周长是最小的,则形成新的
TSP 环路路径序列,重复上述过程;如果 TM 中当前最大的边插入任意路径序列都没有新周
长小于插入之前周长的,则选择 TM 中的下一条次长边,重复上述过程,直到 TM 中的所有
便被处理完毕为止,多边减量迭代的非形式化算法如算法 4 所示。
在算法 4 中,用边向量 Circle 记录 TSP 环路解路径序列,其数据类型为 TspEdge,
Circle[i].p1=Circle[(i+1)%PointNumber].p2(i=0,1,2,…,PointNumber-1)。用二维数组
AdjencyM 记录 V- {ax(m-1) , axm , ax(m+1) ,ax(m+2) }集合中任意边对间的距离(两个边对间
可以有任意多条边),对于数据元素 AdjencyM[i][j]记录 axiax(i+1)…axj 序列从 axi 到 axj 点对间的
距离,其值为 0 表示序列不合法,包含了{ax(m-1) , axm , ax(m+1) ,ax(m+2) }集合中元素。flags
标志向量用来表示 Circle 中的边是否已经被选过作为最长边,为 1 表示已经被选过,为 0
表示未被选过。
算法 4:
⑴ tt=0; , 0i jsum = , flags 向量各元素清 0。
⑵ 如果 tt< PointNumber,则转⑶;否则,转向⑹。
⑶ 在 Circle 向量中找出权值最大且未处理过的数据元素(flags 分量为 0 表示未处理),
其指针为 imax。flags[imax]=1,生成 AdjencyM。
⑷ 对于任意边对 i 和 j(i 和 j 指 CircleTsp 向量元素指针),令
delij=Circle[(i+PointNumbers-1)%PointNumbers].length+Circle[(j+1)%PointNumbers].lengt
h+Circle[imax].length;
a = Circle[(i+PointNumbers-1)%PointNumbers].p1;
b = Circle[(j+1)%PointNumbers].p2;
a1= Circle[imax].p1, b1= Circle[i].p1, a2= Circle[imax].p2,b2= Circle[j].p2
addi,j=len(a, b)+min(len(a1, b1)+ len(a2, b2)), len(a1, b2)+ len(a2, b1)))
在 AdjencyM 中找到满足下列条件的边对:
公式 3:
, , ,,
max( )( , { | 0 int 1})i j i j i ji j del add i j x x Po erNumber∆ = − ∈ ≤ ≤ −
我们把 ,i j∆ 称为边减量,如果 ,i j∆ >0,则转向步骤⑸;否则,t++,转向步骤⑵。
⑸ sum=sum- ,i j∆ ,更新 Circle 向量,flag 转向步骤⑵。
⑹ 结束。
- 12 -
假设平均的多边减量为 ,i j
−∆ ,则算法的时间复杂度为 , 2 3
,
max( , )i j
i j
sum
n n−∆
。因为一次迭
代时在第⑶步生成任意两点的距离 AdjencyM 的时间复杂度为 O(n2),在第⑷步在 AdjencyM
中找任意点对的最大距离的时间复杂度也为 O(n2),算法的迭代次数最大为 ,
,
i j
i j
sum
−∆
+n,因此,
算法的时间复杂度为 O(( ,
,
i j
i j
sum
−∆
+n)n2),也即 , 2 3
,
max( , )i j
i j
sum
n n−∆
。
基于贪心环路的减量迭代 TSP 优化算法的时间复杂度
综合算法 1、算法 2、算法 3 和算法 4 的时间复杂度,基于贪心环路的减量迭代 TSP 优
化算法的时间复杂度为
O( , 2 3
,
max(( ) , )i j e k
i j e k
sum sum sum n n
n
− − −+ +∆ ∆ ∆
)。
4. 实例数据
本算法在对称 TSPLIB 中 12 种实例和 CHN144 总计 13 种实例数据进行检测,检测的结
果如表 1 所示。
表 1 TSP 算法实例运行结果一览表
TSP 实力名称 本算法求得的最短环路 目前已知最短环路
eil101 620 629
eil76 533 538
brazil58 18456 25395
bayg29 1652 1610
bier127 120311 118282
ch130 6284 6110
CHN144 30626 30353
d198 16129 15780
fl417 12108 11861
gil262 2403 2378
kroA1503 26613 26524
kroB150 26918 26130
linhp318 43093 41345
在表 1 中,eil101 、eil76 和 brazil58 三个实例的运行结果明显优于目前已知的最好结果,
它们对应的路径序列如表 2 所示。在表 2 中,TSP 环路路径序列一列的数据是前述的三元组
表示的边序列。值得注意的是 eil101 和 eil76 中在 TSPLIB 给出的数据是坐标值,序号从 1
开始,本算法三元组给出的从 0 开始的两点序号及其对应的距离,而 brazil58 在 TSPLIB 给
出的数据是基于任意两点间的距离邻接矩阵。
- 13 -
表 2 TSP 问题实例环路路径序列表
实例名称 TSP 环路路径序列
eil101
(97,36,1),
(36,99,2), (99,90,3), (90,43,6), (43,13,5), (13,37,11),
(37,85,13), (85,15,6), (15,60,4), (60,4,6), (4,83,4),
(83,16,6), (16,44,8), (44,45,10), (45,7,9), (7,82,7),
(82,59,4), (59,17,6), (17,88,7), (88,51,8), (51,68,10),
(68,0,4), (0,49,6), (49,75,5), (75,76,4), (76,2,2),
(2,78,3), (78,32,5), (32,80,2), (80,8,6), (8,50,6),
(50,29,9), (29,69,5), (69,30,7), (30,87,5), (87,6,6),
(6,81,5), (81,47,5), (47,46,6), (46,35,7), (35,48,8),
(48,63,12), (63,10,13), (10,18,7), (18,61,9), (61,9,6),
(9,62,9), (62,89,4), (89,31,4), (31,19,10), (19,65,8),
(65,64,13), (64,70,10), (70,34,6), (34,33,10), (33,77,5),
(77,28,10), (28,23,7), (23,67,9), (67,79,2), (79,11,6),
(11,25,7), (25,27,8), (27,26,6), (26,100,5), (100,52,4),
(52,57,5), (57,39,4), (39,20,7), (20,72,3), (72,71,3),
(71,3,8), (3,53,9), (53,54,8), (54,24,3), (24,38,9),
(38,66,9), (66,22,12), (22,55,7), (55,74,4), (74,73,3),
(73,21,3), (21,40,4), (40,1,12), (1,56,5), (56,14,7),
(14,42,7), (42,41,9), (41,86,7), (86,12,7), (12,5,7),
(5,93,3), (93,94,3), (94,96,3), (96,91,3), (91,58,2),
(58,95,3), (95,98,2), (98,92,2), (92,84,2), (84,97,3)
eil76
(45,33,2),
(33,51,4), (51,26,4), (26,12,7), (12,53,7), (53,56,14),
(56,14,4), (14,4,8), (4,36,7), (36,19,6), (19,69,6),
(69,59,4), (59,70,5), (70,68,9), (68,35,7), (35,46,6),
(46,20,6), (20,60,10), (60,21,12), (21,61,8), (61,27,6),
(27,73,6), (73,1,7), (1,29,7), (29,47,7), (47,28,6),
(28,44,4), (44,3,7), (3,66,7), (66,25,5), (25,75,6),
(75,74,3), (74,67,4), (67,5,5), (5,50,6), (50,16,6),
(16,39,6), (39,11,5), (11,57,10), (57,71,5), (71,38,5),
(38,8,4), (8,31,7), (31,43,5), (43,2,3), (2,15,9),
(15,62,6), (62,32,6), (32,72,5), (72,0,5), (0,42,6),
(42,41,5), (41,63,9), (63,40,12), (40,55,10), (55,22,5),
(22,48,10), (48,23,7), (23,17,13), (17,49,6), (49,24,8),
(24,54,9), (54,30,21), (30,9,13), (9,37,7), (37,64,5),
(64,10,7), (10,65,7), (65,58,15), (58,13,10), (13,52,7),
(52,6,8), (6,34,5), (34,18,7), (18,7,7), (7,45,5)
brazil58
(54,32,290), (32,17,114), (17,9,98), (9,3,947), (3,12,806),
(12,41,211), (41,34,107), (34,42,171), (42,26,184), (26,2,158),
(2,20,276), (20,31,742), (31,50,297), (50,57,102), (57,53,112),
(53,27,283), (27,49,359), (49,16,129), (16,23,388), (23,24,259),
(24,43,125), (43,11,257), (11,7,436), (7,44,209), (44,29,656),
(29,0,307), (0,15,367), (15,19,1062), (19,1,98), (1,6,137),
(6,4,524), (4,22,523), (22,38,600), (38,5,124), (5,45,138), (45,40,184),
(40,47,217), (47,52,239), (52,28,354), (28,51,274), (51,18,178),
(18,30,691), (30,10,455), (10,21,696), (21,25,359), (25,37,528),
(37,8,343), (8,35,228), (35,13,137), (13,39,176), (39,14,262), (14,36,528),
(36,33,300), (33,55,257), (55,48,199), (48,56,97), (56,46,86), (46,54,72)
为便于观察算法 2、算法 3 和算法 4 对算法 1 的解的优化情况,表 3 列出了表 2 经过算
- 14 -
法 1 到算法 4 解的质量的变化过程。
表 3 TSP 问题实例经过优化算法求解后解的质量情况一览表
实例名称 算法 1 算法 2 算法 3 算法 4
eil101 775 659 655 620
eil76 583 569 555 533
brazil58 25956 22526 20876 18456
5. 结束语
本算法采用不同的构造算法,即借用 Kruskal 最小生成树算法的思想,按照最短边优先
的次序,依次构造贪心环路中的各条边。为便于算法的实现,对点和边进行了适合本算法的
抽象描述。在此基础上,对初始贪心环路进行了减量迭代优化,使得贪心环路解的质量得到
最大限度的提高。实测表明,该算法得到部分实例的结果比目前已知的解要好得多。值得进
一步研究的问题是实例数据具有何种性质采用本算法的质量才能显著是我们下一步工作。
参考文献
1. M R Garey, D S Johnson. Computers and Intractability: A Guide to the Theory of NP-Completeness, San
Francisco: W H Freeman, 1979
2. D S Johnson, L A Mc Geoch. The traveling salesman problem: A case study in local optimization, In: E H L A
arts, J K Lenstra eds. Local Search in Optimization. New York: John Wiley and Sons, 1997. 215~ 310.
3. G Clarke, J W W right. Scheduling of vehicles from a central depot to a number of delivery points. Operations
Res, 1964, 12 (4) : 568~ 581
4. N Christofides Worst-case analysis of a new heuristic for the traveling salesman problem. Graduate Schoo l of
Industrial Administration, Carnegie-Mellon University, Pittsburgh, PA , Tech Rep: 388, 1976
5. S Kirkpatrick, C D Gelatt, M P Vecch i. Optimization by simulated annealing. Science, 1983, 220 (4598) :
671~ 680
6. J H Holland. Adaptation in Natural and Artificial System s. Ann Arbor: University of Michigan Press, 1975
7. G A Croes. A method for solving traveling salesman problem s. Operations Res, 1958, 6 (6) : 791~ 812
8. S Lin. Computer solutions of the traveling salesman problem. Bell Syst Tech J , 1965, 44: 2245~ 2269
9. S Lin, BW Kernighan. A n effective heuristic algorithm fo r the traveling2salesman problem. Operations Res,
1973, 21 ( 2) :498~ 516
10. D S Johnson. Local optimization and the traveling salesman problem. In: P roc of the 17th Colloq on Automata,
Language, and Programming, Lecture Notes in Computer Science 443. Berlin: Sp ringer-Verlag, 1990. 446~461.
11. 陈国良,谢幸,徐云,顾钧. 随机算法重启策略的构造及其在TSP中的应用. 计算机学报. 2002, 25(5):
514~519
12. 金炳尧,蔚承建,何振亚. 一个用于优化搜索的学习算法. 软件学报. 2001, 12(3):448~453
13. 邹鹏,周智,陈国良,顾钧. 求解TSP问题的多级归约算法. 软件学报. 2003, 14(1):35~42
14. 万颖瑜,周智,陈国良,顾钧. SizeScale: 求解旅行商问题(TSP) 的新算法. 计算机研究与发展. 2002,
39(10):1294~1302
15. 朱文兴, 傅清祥. 一个基于填充函数变换的对称TSP问题的局部搜索算法. 计算机学报.2002, 25(7):
701~707
16. 邹鹏,周智,江贺,陈国良,顾钧.求解旅行商问题的循环局部搜索算法的运行时间和性能分布分析. 计算机
学报. 2006, 29(1):92~99
17. 蔡之华,彭锦国,高伟,魏巍,康立山. 一种改进的求解TSP 问题的演化算法. 计算机学报. 2005,
28(5):824~828
18. 杨 辉 康立山 陈毓屏.一种基于构建基因库求解TSP问题的遗传算法.计算机学报. 2006,
26(12):1753~1758
19. 吴斌,史忠植. 一种基于蚁群算法的TSP问题分段求解算法. 计算机学报. 2001, 24(12): 1329~1333
- 15 -
A new optimization algorithm of iterative decrease on
Traveling Salesman Problem based on greed circuit
Li Hongbo,Ma Wenjun,Chen Jun
College of Management,LuDong University,Yantai,Shandong (264025)
College of Mathematics and Information Science,LuDong University,Yantai,Shandong (264025)
Abstract
Use for reference Kruskal algorithm on the minimum spanning tree to produce the edges in the solution
of Traveling Salesman Problem and build the initial solution of the greed circuit in order of precedence
of the shortest edge. Then, make the first optimization arrangement with the iterative decrease of point
and the twice optimization arrangement with the iterative decrease of single edge. Finally, continue the
third optimization arrangement with the iterative decrease of multiple edges. In order to convenience of
computation, status of each point is described as a six tuple, and each edge is characterized as a three
tuple, and then these coordinate data and adjacency matrix data in TSPLIB is reconstructed as
adjacency matrix of upper triangle. Tested on thirteen illustrations, one part of which is twelve
illustrations in TSPLIB and another is CHN144, the new algorithm of which is able to largely decrease
these known best solutions and eil101, eil76 and brazil58 are deceased from 629 to 620, from 538 to
533 and from 25395 to 18456 respectively, at the same time show the circuit path sequence of their
solution on Traveling Salesman Problem. The time complexity of the new algorithm is
O( , 2 3
,
max(( ) , )i j e k
i j e k
sum sum sum n n
n
− − −+ +∆ ∆ ∆
).
Keywords:TSP greed circuit;iterative decrease;decrease of point;decrease of single edge;decrease
of multiple edges
作者简介:李洪波(1969-),男,汉族,辽宁昌图人,讲师,硕士,主要研究方向为算法分
析与设计、智能信息管理。