《MATLAB智能优化算法》学习笔记二
一、TSP概述 旅行商问题,即TSP问题(Traveling Salesman Problem) 又译为旅行推销员问题、货郎担问题,是数学领域中著名问题之一。假设有一个旅行商人要拜访n个城市,他必须选择所要走的路径,路径的限制是每个城市只能拜访一次,而且最后要回到原来出发的城市。路径的选择目标是要求得的路径路程为所有路径之中的最小值(等价于求图的最短哈密尔顿回路问题)。
举个栗子,下表是四个城市的位置信息:
因为一共需要访问4个城市,所以全部 的行走路线方案的数目为$A_4^4=24$。因为数据规模较小,所以把所有的行走路线方案全部列出来。如下表所示,其中黄色底色的方案表示行走距离最短的行走路线方案。
进一步观察发现,方案1、10、17、和19其实是一种方案,方案6、8、15和24其实也是同一种方案。虽然看似方案的起点和终点不同,但是旅行商人的行走路线实际上是一条闭环回路。
因为一条行走路线的起点和终点是同一个城市,所以在表示该路线时为了表现得更加简洁,将终点城市从行走路线中省略。以12341这条行走路线为例,在省略终点城市1后,该条行走路线变为1234。后面部分所有关于路线的表现形式全部省略终点城市。
二、算法简介 就TSP而言,最终的目的是为旅行商人找到条最短路线。因为一开始最短路线是未知的,所以只能先尝试随便为旅行商人设计出一条一个城市只能访问一次的行走路线。假设这条路线不是最短的那条路线,那么这条路线一定还可以通过调整旅行商人访问城市的顺序,从而达到缩短行走总距离的目的。因此,调整旅行商人访问城市顺序的方法就显得尤为重要,不同的调整方法对于缩短行走总距离的效果是不同的。
以上节中1243这条路线为例,由表可知,这条路线的行走总距离为41.4km。因为在上节已知最短总距离为33.3km,所以1243这条路线一定还有调整的空间。那究竟如何调整1243路线呢?二个容易想到的策略是对调相邻城市的访问顺序。因此,采用该策略对1243路线调整后的路线分别如下。
调整路线1:2143(将1和2对调).此时总距离为33.3km。
调整路线2:1423(将2和4对调),此时总距离为33.9km。
调整路线3:1234(将4和3对调).此时总距离为33.3km。
调整路线4:3241(将3和1对调),此时总距离为33.9km。
这里需要注意的一点是,3和1其实也是相邻城市,因此也可以进行对调。接下来从上述4条路线中选择最短的一条路线,因为调整路线1和调整路线3总距离是相等的,这里在经过随机选择后,假设选择的是调整路线3。
因为此时不知道1234这条路线是否为最短路线,所以还会继续对这条路线进行调整。此时用一种新的策略调整1234这条路线,即对调中间间隔一个城市的两个城市的访问顺序。采用该策略对1234路线调整后的路线分别如下。
调整路线5:3214(将1和3对调),此时总距离为33.3km。
调整路线6:1432(将2和4对调).此时总距离为33.3km。
虽然这2条路线总距离都为33.3km,但是此时依然不确定是否找到最短路线。因此,还会从调整路线5和调整路线6这2条路线中随机选择一条路线,然后继续对该路线进行调整,一直达到初始设定的条件为止,才会停止继续调整与寻找更短的路线。假设最初设定只允许调整路线3次,那么在调整3次路线后,就会自动停止寻找更短路线,此时所找到的最短路线就记为最短路线。
综上所述,对调相邻城市的访问顺序的策略称为一个规则,记为规则1;对调中间间隔一个城市的两个城市的访问顺序也是一个规则,记为规则 。第1、2、3、4条调整路线实际上表示一个集合,记为集合1;第5、6条调整路线实际上也表示一个集合,记为集合2只有按照规则对原始路线进行操作后才能形成集合,形成集合的目的是从集合中找出比原始路线更短的路线。
至此,引出“变邻域搜索”的概念。规则对应邻域动作,集合对应邻域,不同的集合体现出变邻域中的变.从集合中找出比原始路线更短的路线对应搜索。变邻域搜索就是在不同的领域中不断地搜索更优的解。 VNS 算法求解组合优化问题(以求最小化问题为例)的常规流程如图所示。
三、算法设计 (1)构造初始路线
就 TSP 而言,一条好的初始路线能够节省VNS的求解时间。在这采用常规且应用较为广泛的构造 TSP 初始路线的方法一贪婪算法。
假设 TSP 中城市数目为 N ,贪婪算法构造 TSP 初始路线的步骤如下。
STEP1:初始化已被访问的城市集合 visited 为空,初始化未被访问的城市集合 unvisited ={1,2…, N}。
STEP2:将 N 行 N 列距离矩阵 dist 的主对角线上的0全部赋值为无穷大。
STEP3:从距离矩阵 dist 中找出最小距离对应的行序号 row 和列序号 col ,如果存在多个最小距离,则选择行序号 row 中的第一个数 row (1)作为初始路线的起点 first 。
STEP4:更新 visited =[ visited , first ],即将 first 添加到 visited 中,同时也更新 unvisited ( unvisited = first ) = [ ] 即将 first 从 unvisited 中删除,然后将起点 first 赋值给紧前点 pre_point 。
STEP5:首先在距离矩阵 dist 中找到紧前点 pre_point 对应的那一行距离 pre_dist ,即紧前点 pre_point 与其他城市之间的距离;其次将已被访问的城市排除在外,即在 pre_dist 中将 visited 对应的列全设为无穷大;然后找出 pre_dist 中最小值对应的列序号作为下一个紧前点 pre_point 。
STEP6:更新 visited =[ visited , pre_point ], unvisited ( unvisited = pre_point ) = [ ]。
STEP7:若 unvisited 非空,则转至STEP5,否则转至STEP8。
STEP8:将 visited 赋值给初始路线 init_route ,初始路线构造完毕。
(2)交换操作
交换前路线: 1 2 3 4 5 6
若选择的交换位置$i=2$和$j=5$,那么交换第二位和第五位的两座城市
交换后路线: 1 5 3 4 2 6
(3)逆转操作
交换前路线: 1 2 3 4 5 6
若选择的逆转位置$i=2$和$j=5$,那么将位置二和五之间的城市顺序逆转
交换后路线: 1 5 3 4 2 6
(4)插入操作
例1:
交换前路线: 1 2 3 4 5 6
若选择的插入位置$i=2$和$j=5$,那么将第二位的城市插到第五位城市的后面
交换后路线: 1 3 4 5 2 6
例2:
交换前路线: 1 2 3 4 5 6
若选择的插入位置$i=5$和$j=2$,那么将第五位的城市插到第二位城市的后面
交换后路线: 1 2 5 3 4 6
(5)扰动操作
当使用某个邻域操作(交换、逆转和插入)前就先对当前路线使用该邻域操作以获得一个“扰动解”,然后再在这个“扰动解”上进行后续一些列操作。
(6)邻域搜索策略
在对一个当前解$S_{curr}$使用某个邻域操作得到$S_{curr}$的邻域集合后,如何对所得到的邻域集合进行处理才能使$S_{curr}$向“更好”的方向变换?
既然已经得到$S_{curr}$的邻域集合,那么就可以首先求出这个邻域集合中所有解的总距离;其次找到总距离最短的那个解$S_{min}$,并替换当前解$S_{curr}$,即令$S_{curr}=S_{min}$;然后求出$S_{curr}$的邻域集合,以及这个邻域集合中所有解的总距离,同样再找到总距离最短的那个解$S_{min}$,并替换当前解$S_{curr}$。
一直按照上述方式迭代,直至迭代 M 次后,停止对$S_{curr}$在当前邻域的搜索。
(7)邻域变化策略
邻域搜索策略详细介绍了如何针对一个当前解$S_{curr}$在一个邻域中进行搜索,因为本节介绍了3个邻域动作,所以会有3个不同结构的邻域。现在将这3个邻域分别编号为 k = 1、k = 2和 k = 3,即第1个邻域为交换操作得到的邻域,第2个邻域为逆转操作得到的邻域,第3个邻域为插人操作得到的邻域。
那么在对一个邻域搜索结束后,如何能够在下一个邻域中继续搜索呢?假设此时当前解为$S_{curr}$, $S_{curr}$的总距离为为$L_{curr}$,令最优路线$S_{best}=S_{curr}$,最优路线总距离$L_{best}=L_{curr}$
具体步骤如下。
STEP1:设 k = 1。
STEP2:如果 k = 1,转至STEP3;如果 k = 2,转至STEP4;如果 k = 3,转至STEP5;否则,转至STEP7。
STEP3:对$S_{curr}$进行扰动操作得到$S_{swap}$ ,$S_{swap}$的总距离为$L_{swap}$,令$L_{curr}=L_{swap}$。如果 $L_{curr}<L_{best}$ 则 $S_{best}=S_{curr}=S_{swap}$, $L_{best}=L_{curr}$ , k = 0,转至STEP6。
STEP4:对$S_{curr}$进行扰动操作得到$S_{reversion}$, $S_{reversion}$的总距离为$L_{reversion}$,令$L_{curr}=L_{reversion}$ 如如果 $L_{reversion}$ < $L_{best}$ ,则$S_{best}=S_{curr}=S_{reversion}$ , $L_{best}=L_{curr}$ , k =0,转至STEP6。
STEP5:对$S_{curr}$进行扰动操作得到$S_{insertion}$ , $S_{insertion}$ 的总距离为$L_{insertion}$ ,令$L_{curr}=L_{insertion}$ 。如果$L_{curr}<L_{best}$ ,则$S_{best}=S_{curr}=S_{insertion}$ , $L_{best}=L_{curr}$ , k =0,转至STEP6。
STEP6:k = k +1,转至STEP2。
STEP7:终止循环,输出$S_{curr}、S_{best}、L_{curr}$ 和$L_{best}$
(8)VNS求解TSP流程
四、MATLAB程序实现 (1)构造初始路线函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 function [init_route,init_len] =construct_route (dist) N=size (dist,1 ); for i =1 :N for j =1 :N if i ==j dist(i ,j )=inf ; end end end unvisited=1 :N; visited=[]; min_dist=min (min (dist)); [row,col]=find (dist==min_dist); first=row(1 ); unvisited(unvisited==first)=[]; visited=[visited,first]; pre_point=first; while ~isempty (unvisited) pre_dist=dist(pre_point,:); pre_dist(visited)=inf ; [~,pre_point]=min (pre_dist); unvisited(unvisited==pre_point)=[]; visited=[visited,pre_point]; end init_route=visited; init_len=route_length(init_route,dist); end
(2)路线总距离计算函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 function len =route_length (route,dist) n=numel (route); route=[route route(1 )]; len=0 ; for k=1 :n i =route(k); j =route(k+1 ); len=len+dist(i ,j ); end end
(3)交换操作函数:
1 2 3 4 5 6 7 8 9 10 function route2 =swap (route1,i,j) route2=route1; route2([i j ])=route1([j i ]); end
(4)逆转操作函数:
1 2 3 4 5 6 7 8 9 10 11 12 function route2 =reversion (route1,i,j) i1=min ([i ,j ]); i2=max ([i ,j ]); route2=route1; route2(i1:i2)=route1(i2:-1 :i1); end
(5)插入操作函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 function route2 =insertion (route1,i,j) if i <j route2=route1([1 :i -1 i +1 :j i j +1 :end ]); else route2=route1([1 :j i j +1 :i -1 i +1 :end ]); end end
(6)计算距离差值函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 function delta1 =cal_delta1 (route,dist,i,j) N=numel (route); if (i ==1 )&&(j ==N) delta1=-(dist(route(i ),route(i +1 ))+dist(route(j -1 ),route(j )))+... (dist(route(j ),route(i +1 ))+dist(route(j -1 ),route(i ))); elseif (i ==1 )&&(j ==2 ) delta1=-(dist(route(N),route(i ))+dist(route(j ),route(j +1 )))+... (dist(route(N),route(j ))+dist(route(i ),route(j +1 ))); elseif i ==1 delta1=-(dist(route(N),route(i ))+dist(route(i ),route(i +1 ))+... dist(route(j -1 ),route(j ))+dist(route(j ),route(j +1 )))+... (dist(route(N),route(j ))+dist(route(j ),route(i +1 ))+... dist(route(j -1 ),route(i ))+dist(route(i ),route(j +1 ))); elseif (i ==N-1 )&&(j ==N) delta1=-(dist(route(i -1 ),route(i ))+dist(route(j ),route(1 )))+... (dist(route(i -1 ),route(j ))+dist(route(i ),route(1 ))); elseif j ==N delta1=-(dist(route(i -1 ),route(i ))+dist(route(i ),route(i +1 ))+... dist(route(j -1 ),route(j ))+dist(route(j ),route(1 )))+... (dist(route(i -1 ),route(j ))+dist(route(j ),route(i +1 ))+... dist(route(j -1 ),route(i ))+dist(route(i ),route(1 ))); elseif abs (i -j )==1 delta1=-(dist(route(i -1 ),route(i ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(j ))+dist(route(i ),route(j +1 ))); else delta1=-(dist(route(i -1 ),route(i ))+dist(route(i ),route(i +1 ))+... dist(route(j -1 ),route(j ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(j ))+dist(route(j ),route(i +1 ))+... dist(route(j -1 ),route(i ))+dist(route(i ),route(j +1 ))); end end
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 function delta2 =cal_delta2 (route,dist,i,j) N=numel (route); if i ==1 if j ==N delta2=0 ; else delta2=-dist(route(j ),route(j +1 ))-dist(route(N),route(i ))+... dist(route(i ),route(j +1 ))+dist(route(N),route(j )); end else if j ==N delta2=-dist(route(i -1 ),route(i ))-dist(route(1 ),route(j ))+... dist(route(i -1 ),route(j ))+dist(route(i ),route(1 )); else delta2=-dist(route(i -1 ),route(i ))-dist(route(j ),route(j +1 ))+... dist(route(i -1 ),route(j ))+dist(route(i ),route(j +1 )); end end end
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 function delta3 =cal_delta3 (route,dist,i,j) N=numel (route); if i <j if (i ==1 ) && (j ==N) delta3=0 ; elseif (i ==1 ) && (j ==2 ) delta3=-(dist(route(N),route(i ))+dist(route(j ),route(j +1 )))+... (dist(route(N),route(j ))+dist(route(i ),route(j +1 ))); elseif i ==1 delta3=-(dist(route(N),route(i ))+dist(route(i ),route(i +1 ))+dist(route(j ),route(j +1 )))+... (dist(route(N),route(i +1 ))+dist(route(j ),route(i ))+dist(route(i ),route(j +1 ))); elseif (i ==N-1 )&&(j ==N) delta3=-(dist(route(i -1 ),route(i ))+dist(route(j ),route(1 )))+... (dist(route(i -1 ),route(j ))+dist(route(i ),route(1 ))); elseif j ==N delta3=-(dist(route(i -1 ),route(i ))+dist(route(i ),route(i +1 ))+dist(route(j ),route(1 )))+... (dist(route(i -1 ),route(i +1 ))+dist(route(j ),route(i ))+dist(route(i ),route(1 ))); elseif (j -i )==1 delta3=-(dist(route(i -1 ),route(i ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(j ))+dist(route(i ),route(j +1 ))); else delta3=-(dist(route(i -1 ),route(i ))+dist(route(i ),route(i +1 ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(i +1 ))+dist(route(j ),route(i ))+dist(route(i ),route(j +1 ))); end else if (i ==N) && (j ==1 ) delta3=-(dist(route(i -1 ),route(i ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(j ))+dist(route(i ),route(j +1 ))); elseif (i -j )==1 delta3=0 ; elseif i ==N delta3=-(dist(route(i -1 ),route(i ))+dist(route(i ),route(1 ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(1 ))+dist(route(j ),route(i ))+dist(route(i ),route(j +1 ))); else delta3=-(dist(route(i -1 ),route(i ))+dist(route(i ),route(i +1 ))+dist(route(j ),route(j +1 )))+... (dist(route(i -1 ),route(i +1 ))+dist(route(j ),route(i ))+dist(route(i ),route(j +1 ))); end end end
(7)更新距离差值函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 function Delta1 =Update1 (route,dist,i,j) N=numel (route); route2=swap(route,i ,j ); Delta1=zeros (N,N); for i =1 :N-1 for j =i +1 :N Delta1(i ,j )=cal_delta1(route2,dist,i ,j ); end end
1 2 3 4 5 6 7 8 9 10 11 12 13 14 function Delta2 =Update2 (route,dist,i,j) N=numel (route); route2=reversion(route,i ,j ); Delta2=zeros (N,N); for i =1 :N-1 for j =i +1 :N Delta2(i ,j )=cal_delta2(route2,dist,i ,j ); end end
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 function Delta3 =Update3 (route,dist,i,j) N=numel (route); route2=insertion(route,i ,j ); Delta3=zeros (N,N); for i =1 :N for j =1 :N if i ~=j Delta3(i ,j )=cal_delta3(route2,dist,i ,j ); end end end
(8)扰动操作函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 function [route_shake,len_shake] =shaking (route,dist,k) N=numel (route); select_no=randi([1 ,N],1 ,2 ); i =select_no(1 );j =select_no(2 );if k==1 route_shake=swap(route,i ,j ); elseif k==2 route_shake=reversion(route,i ,j ); else route_shake=insertion(route,i ,j ); end len_shake=route_length(route_shake,dist); end
(9)交换邻域搜索函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 function [swap_route,swap_len] =swap_neighbor (route,dist,M) N=numel (route); Delta1=zeros (N,N); for i =1 :N-1 for j =i +1 :N Delta1(i ,j )=cal_delta1(route,dist,i ,j ); end end cur_route=route; m=1 ; while m<=M min_value=min (min (Delta1)); if min_value<0 [min_row,min_col]=find (Delta1==min_value); Delta1=Update1(cur_route,dist,min_row(1 ),min_col(1 )); cur_route=swap(cur_route,min_row(1 ),min_col(1 )); else break end m=m+1 ; end swap_route=cur_route; swap_len=route_length(swap_route,dist); end
(10)逆转邻域搜索函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 function [reversion_route,reversion_len] =reversion_neighbor (route,dist,M) N=numel (route); Delta2=zeros (N,N); for i =1 :N-1 for j =i +1 :N Delta2(i ,j )=cal_delta2(route,dist,i ,j ); end end cur_route=route; m=1 ; while m<=M min_value=min (min (Delta2)); if min_value<0 [min_row,min_col]=find (Delta2==min_value); Delta2=Update2(cur_route,dist,min_row(1 ),min_col(1 )); cur_route=reversion(cur_route,min_row(1 ),min_col(1 )); else break end m=m+1 ; end reversion_route=cur_route; reversion_len=route_length(reversion_route,dist); end
(11)插入邻域搜索函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 function [insertion_route,insertion_len] =insertion_neighbor (route,dist,M) N=numel (route); Delta3=zeros (N,N); for i =1 :N-1 for j =i +1 :N Delta3(i ,j )=cal_delta3(route,dist,i ,j ); end end cur_route=route; m=1 ; while m<=M min_value=min (min (Delta3)); if min_value<0 [min_row,min_col]=find (Delta3==min_value); Delta3=Update3(cur_route,dist,min_row(1 ),min_col(1 )); cur_route=insertion(cur_route,min_row(1 ),min_col(1 )); else break end m=m+1 ; end insertion_route=cur_route; insertion_len=route_length(insertion_route,dist); end
(12)TSP路线图函数:
1 2 3 4 5 6 7 8 9 10 function plot_route (route,x,y) figure route=[route route(1 )]; plot (x(route),y(route),'k-o' ,'MarkerSize' ,10 ,'MarkerFaceColor' ,'w' ,'LineWidth' ,1.5 );xlabel('x' ); ylabel('y' ); end
(13)主函数:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 tic clear clc dataset=importdata('input.txt' ); x=dataset(:,2 ); y=dataset(:,3 ); vertexes=dataset(:,2 :3 ); N=size (dataset,1 ); h=pdist(vertexes); dist=squareform(h); MAXGEN=50 ; M=50 ; n=3 ; [init_route,init_len]=construct_route(dist); disp (['初始路线总距离为' ,num2str(init_len)]);cur_route=init_route; best_route=cur_route; best_len=route_length(cur_route,dist); BestL=zeros (MAXGEN,1 ); gen=1 ; while gen<=MAXGEN k=1 ; while (1 ) switch k case 1 cur_route=shaking(cur_route,dist,k); [swap_route,swap_len]=swap_neighbor(cur_route,dist,M); cur_len=swap_len; if cur_len<best_len cur_route=swap_route; best_len=cur_len; best_route=swap_route; k=0 ; end case 2 cur_route=shaking(cur_route,dist,k); [reversion_route,reversion_len]=reversion_neighbor(cur_route,dist,M); cur_len=reversion_len; if cur_len<best_len cur_route=reversion_route; best_len=cur_len; best_route=reversion_route; k=0 ; end case 3 cur_route=shaking(cur_route,dist,k); [insertion_route,insertion_len]=insertion_neighbor(cur_route,dist,M); cur_len=insertion_len; if cur_len<best_len cur_route=insertion_route; best_len=cur_len; best_route=insertion_route; k=0 ; end otherwise break ; end k=k+1 ; end disp (['第' ,num2str(gen),'代最优路线总距离为' ,num2str(best_len)]); BestL(gen,1 )=best_len; gen=gen+1 ; end figure ;plot (BestL,'LineWidth' ,1 );title('优化过程' ) xlabel('迭代次数' ); ylabel('总距离' ); plot_route(best_route,x,y); toc
五、实验结果展示
1 . 《MATLAB智能优化算法》(曹旺著) ↩
赞赏
感谢鼓励