TSP问题的动态规划解法

1、TSP 问题的动态规划解法第十七组： 3103038028 郑少斌3103038029 王瑞锋3103038035 江飞鸿3103038043 韩鑫3103055004 唐万强1TSP 问题简介旅行商问题（Traveling Salesman Problem,简称 TSP, 亦称为货单郎问题 ）可以描述为：对于 N 个城市，它们之间的距离已知，有一旅行商要从某一城市走遍所有的城市， 且每一城市只能经过一次， 最后回到出发的城市， 问如何选择路线可使他走过的路径最短。 这是一个典型的组合优化问题。 它有很强的现实意义，可以应用于交通运输，物资调配，旅游线路设置。对于了解某个国家地理分布也有一定的

2、现实意义。 这个问题的解法有很多种，在这里我们尝试使用最优控制中的动态规划的相关知识来进行求解。2TSP 问题分析对于这个问题，我们首先想到的是应用穷举法进行解答，但是这个方法时间和空间的复杂度很高。从表面上看，TSP问题很简单，其实则不然。对于N 个城市的 TSP,存在的可能路径为（ N-1 ）!/2 条，当 N 较大时，其数量是惊人的。计算每条路经都需求出N 个距离之和， 这样各种路径及其距离之和的计算量正比与N!/2. 用搜索法要求就规模大的TSP 是不现实的。例如使用 1GFLOPs次的计算机搜索TSP 所需的时间如下表所示城市数7152050100200加法量2.51036.5101

3、11.210181.5106451015710374搜索时2.510 5 s51048 y10142 y10358y1.8h350y间由上可知，对于这个问题采用穷举法进行解答是不现实的，这就要求我们采用其他的方法进行解答。3其他求解 TSP 问题的方法* 贪心法a. 所谓贪心法，就是在组合算法中，将每一步都取局部最优的求解方法。b. 下表表示用贪心法求解 TSP 的过程。先将各城市间的距离用行列式形式表示，主对角线上用表示。我们可以从城市 C1 出发，依次在每一行或列中选取元素最小的路径， 且每个城市只能访问一次。c. 按贪心法从 C1 出发所挑选的路径为C1C7C6C2C3C4C5C1L27

4、34431033不难看出，这种从局部最优原则出发的方法所得的结果的好坏，与城市间的距离的具体情况和从那个城市开始有关。例如，从 C7 出发时，用贪心法所得的结果为C7C1C4C5C3C2C6C7其路线长度减为L253743731*Hopfield 神经元网络法a. 全互连型神经网络求解TSP 问题：设有 N 个城市：C=C1C 2 .Cn C 中任意两个城市的距离D( cic j )= dij现在要找到一个城市的排列cn (1)cn (2).cn ( N )使得闭合路径d cn(i )cn (i1) mod n为最小我们用换位矩阵来表征一条路径。对于N 个城市来说，换位矩阵有 N*N 个元素，

5、需要用N*N 个神经元来表征。ABCDEA00001B10000C00100D01000E00010根据下面四方面的要求， 即：1.换位矩阵每行只能有一个一；2.换位矩阵每列只能有一个一；3.换位矩阵中元素一之和应为N；4.所构造的函数的极小值对应于最短路径。我们构造出与 TSP 相对应的计算能量函数为aNN1Nb N N1 NcEvxi vxj2 i 1 xvxi v yi22 x 1 i 1 j i 11 y i 1NN2dvxiNx 1i 12NNNd xyvxi (vy ,i 1vy,i _ 1 )x 1 y 1 i 1式中前三项与条件的1，2，3 的要求对应，上式的第四项为目标项，它

6、的最小值就对应于最短路径长度。上式中 v 的数值为 0 或者 1，是由表正换位矩阵中神经元的输出来表示的。4TSP 问题的动态规划解法主要特点：将一个问题分为若干个相互联系的阶段，每个阶段进行决策优化。在这种多阶段决策优化过程中，无论其初始状态和初始决策如何，以后的最优策略只取决于由最初策略所形成的当前策略。5问题分析我们应用动态规划的解法来求解五个城市的TSP 问题。我们采用一个矩阵A 表示城市之间的距离。0a12a13a14a15a16a17a18a19a110a120a23a24a25a26a27a28a29a210a13a230a34a35a36a37a38a39a310a14a24a

7、340a45a46a47a48a49a410a15a25a35a450a56a57a58a59a510Aa26a36a46a560a67a68a69a610a16a17a27a37a47a57a670a78a79a710a18a28a38a48a58a68a780a89a810a19a29a39a49a59a69a79a890a910a110a210a310a410a510a610a710a810a9100其中， aij ，i0,12 . 910j0,12 .910 表示第 i 个城市和第j 个城市之间的距离。这是个对称阵，而且对角线的元素均为0。假定我们已经找到一个最短的路径，设它是先从 c

8、1 到 ck ，则从 ck 出发沿着这条路径回到c1 的路，必定是经过Cc1, ck中每个城市一次， 最后回到 c1 的路径中的最短的一条，也就是符合最优原理。设 f c1 , S 表示从城市 c1 出发，通过 s 集合中所有城市一次且仅一次，最后回到出发城市 c1 的最短路径的长度。由 f 的定义知，所求最短路径长度可表示为 f c1 ,C c1 。根据最优原理，应有f c1, Cc1 min d1k f ck , Cc1,ck2 k 10一般的有： f ci , S min dij fc j , C c j 。ci S根据以上分析，应用Matlab 编程如下：clearclcD =0146

9、.13356.43286.99349.83;140.950228.38456.76201.68;364.21233.230431.53198.65;287.89471.96418.440222.73;340.68191.78213.68232.220;n=4;c1=1;c2=n*nchoosek(n-1,n-1);c3=n*nchoosek(n-1,n-2);c 4=n*nchoosek(n-1,n-3);c5=n*nchoosek(n-1,n-4);layer1=zeros(1,c1);layer2=zeros(1,c2);layer3=zeros(1,c3);layer4=zeros(1,

10、c4);layer5=zeros(1,c5);path1=0;path2=zeros(1,c2);path3=zeros(1,c3);path4=zeros(1,c4);path5=zeros(1,c5);%#第五层for i=1:1:c5+1layer5(1,i)=D(i,1);end%#第四层layer4(1,1)=D(3,2)+layer5(1,2); %3-2-1layer4(1,2)=D(2,3)+layer5(1,3); %2-3-1layer4(1,3)=D(4,2)+layer5(1,2); %4-2-1layer4(1,4)=D(2,4)+layer5(1,3); %2-4-

11、1layer4(1,5)=D(5,2)+layer5(1,2); %5-2-1layer4(1,6)=D(2,5)+layer5(1,5); %2-5-1layer4(1,7)=D(4,3)+layer5(1,3); %4-3-1layer4(1,8)=D(3,4)+layer5(1,4); %3-4-1layer4(1,9)=D(5,3)+layer5(1,3); %5-3-1layer4(1,10)=D(3,5)+layer5(1,5);%3-5-1layer4(1,11)=D(5,4)+layer5(1,4);%5-4-1 layer4(1,12)=D(4,5)+layer5(1,5);

12、%4-5-1 %#第三层dxy,p=min(D(4,3)+layer4(1,1)D(4,2)+layer4(1,2);%4-(2 3)-1layer3(1,1)=dxy; path3(1,1)=p;dxy,p=min(D(5,3)+layer4(1,1)D(5,2)+layer4(1,2);%5-(2 3)-1layer3(1,2)=dxy; path3(1,2)=p;dxy,p=min(D(3,4)+layer4(1,3)D(3,2)+layer4(1,4);%3-(2 4)-1layer3(1,3)=dxy; path3(1,3)=p;dxy,p=min(D(5,4)+layer4(1,3

13、)D(5,2)+layer4(1,4);%5-(2 4)-1layer3(1,4)=dxy; path3(1,4)=p;dxy,p=min(D(3,5)+layer4(1,5)D(3,2)+layer4(1,6);%3-(2 5)-1layer3(1,5)=dxy; path3(1,5)=p;dxy,p=min(D(4,5)+layer4(1,5)D(4,2)+layer4(1,6);%4-(2 5)-1layer3(1,6)=dxy; path3(1,6)=p;dxy,p=min(D(2,4)+layer4(1,7)D(2,3)+layer4(1,8);%2-(3 4)-1layer3(1,

14、7)=dxy; path3(1,7)=p;dxy,p=min(D(5,4)+layer4(1,7)D(5,3)+layer4(1,8);%5-(3 4)-1layer3(1,8)=dxy; path3(1,8)=p;dxy,p=min(D(2,5)+layer4(1,9)D(2,3)+layer4(1,10);%2-(3 5)-1layer3(1,9)=dxy; path3(1,9)=p;dxy,p=min(D(4,5)+layer4(1,9)D(4,3)+layer4(1,10);%4-(3 5)-1layer3(1,10)=dxy; path3(1,10)=p;dxy,p=min(D(2,

15、5)+layer4(1,11)D(2,4)+layer4(1,12); %2-(4 5)-1layer3(1,11)=dxy; path3(1,11)=p;dxy,p=min(D(3,5)+layer4(1,11)D(3,4)+layer4(1,12); %3-(4 5)-1 layer3(1,12)=dxy; path3(1,12)=p; %#第二层dxy,p=min(D(2,3)+layer3(1,12)D(2,4)+layer3(1,10)D(2,5)+layer3(1,8); %2-(3 4 5)-1layer2(1,1)=dxy;path2(1,1)=p;dxy,p=min(D(3,

16、2)+layer3(1,11)D(3,4)+layer3(1,6)D(3,5)+layer3(1,4); %3-(2 4 5)-1layer2(1,2)=dxy;path2(1,2)=p;dxy,p=min(D(4,2)+layer3(1,9)D(4,3)+layer3(1,5)D(4,5)+layer3(1,2);%4-(2 3 5)-1layer2(1,3)=dxy;path2(1,3)=p;dxy,p=min(D(5,2)+layer3(1,7)D(5,3)+layer3(1,3)D(5,4)+layer3(1,1);%5-(2 3 4)-1layer2(1,4)=dxy;path2(1

17、,4)=p;%#第 一层dxy,p=min(D(1,2)+layer2(1,1)D(1,3)+layer2(1,2)D(1,4)+layer2(1,3) D(1,5)+layer2(1,4); %1-(2 3 4 5)-1layer1(1,1)=dxypath1=p;path2=path2;path3=path3;optimal_path=1 2 3 5 4 1运行结果：layer1 =1.0933e+003optimal_path = 1235416附录附贪心法及神经网络法的程序及相应的运行结果（均为用Matlab 编写）贪心法程序：clearclcD=zeros(10,10);L_min=

18、zeros(1,10);Path=zeros(1,10);Dmin=0;Min_D=zeros(1,10);optimal=0;for i=1:1:10for j=1:1:10D(i,j)=unidrnd(1000);endendfor i=1:1:10for j=1:1:10if i=jD(i,j)=1001;elseD(i,j)=D(j,i);endendendfor n=1:1:10Path(1,1)=n;for i=1:1:9L_min(1,i)=min(D(Path(1,i),:);for y=1:1:10if D(Path(1,i),y)=L_min(1,i)Path(1,i+1)

19、=y;break;endendfor j=1:1:10D(Path(1,i),j)=1001;D(j,Path(1,i)=1001;endD;endfor k=1:1:9Dmin=Dmin+L_min(1,k);endPathDminMin_D(1,n)=Dmin;Dmin=0;endoptimal=min(Min_D)神经网络解法程序：%人工神经网络求解TSP200459clearclcn=10;k1=500;k2=500;k3=500;k4=500;u0=0.02;k=0;%初始参数dxy=0;const=0;u_xi=0;du_xi=0;T=0.01;Dmin=0;flag=0;coun

20、t=0; row=0;column=0;row_s=0;column_s=0; %初始化变量U=zeros(n,n);V=zeros(n,n);D=zeros(n,n);Utemp=zeros(n,n);Vtemp=zeros(n,n);%初始化数组Uini=zeros(n,n);Vini=zeros(n,n);City_Path=zeros(1,n);%初始化数组for x=1:1:nfor y=1:1:nD(x,y)=unifrnd(0.01,1);endendfor x=1:1:nfor y=1:1:nif x=yD(x,y)=0;endD(x,y)=D(y,x);endendN=n*n

21、;uoo=-0.5*u0*log(n-1); % 计算初值 u00while flag=0for i=1:1:nfor j=1:1:nr=unifrnd(-0.1*u0, 0.1*u0);U(i,j)= uoo + r;%按均匀分布随机产生人工神经网络每个神经元初始输入ui，i=1,2.100V(i,j)=0.5 * (1+tanh( U(i,j)/u0 ) ); % 把 ui 代入 S 函数得到每个神经元的初始输出vi ，i=1,2.100endendUini=U;%保存初态uiVini=V;%保存初态vifor k=0:1:200%迭代次数for x=1:1:n%城市循环for i=1:1

22、:n%路径按步循环for y=1:1:n %求微分方程中的最后一个和式： dxy*(V(y,i-1)+V(y,i+1),x,y from 1 to n.if i=1dxy=dxy+D(x,y)* (0+V(y,i+1);%当 i=1 时,不存在左顺联 ,所以 V(y,i-1) 0.elseif i=ndxy=dxy+D(x,y)*(V(y,i-1)+0);%当 i=n 时,不存在顺联 ,所以 V(y,i+1) 0.elsedxy=dxy+D(x,y)*(V(y,i-1)+V(y,i+1); % 当 1in 时,用原公式 .endendconst=-k1*(sum(V(x,:)-V(x,i)-k2*(sum(V(:,i)-V(x,i)- k3*(sum(sum(V)-n) - k4*dxy;% 求微分方程的常数项du_xi=-U(x,i)+const;%求ui对时间t的导数dui/dtu_xi=U(x,i)+du_xi*T;%用Eula法求u(t+T)=u(t)+du/dt * Tv_xi=0.5*(1+tanh(u_xi/u0);%把u(t+T) 代入S函数求t+T时刻的V(t+T)