TSP问题的动态规划解法.doc

TSP问题的动态规划解法 第十七组：3103038028 郑少斌 3103038029 王瑞锋 3103038035 江飞鸿 3103038043 韩鑫 3103055004 唐万强1 TSP问题简介 旅行商问题（Traveling Salesman Problem,简称TSP, 亦称为货单郎问题 ）可以描述为：对于N 个城市，它们之间的距离已知，有一旅行商要从某一城市走遍所有的城市，且每一城市只能经过一次，最后回到出发的城市，问如何选择路线可使他走过的路径最短。这是一个典型的组合优化问题。它有很强的现实意义，可以应用于交通运输，物资调配，旅游线路设置。对于了解某个国家地理分布也有一定的现实意义。这个问题的解法有很多种，在这里我们尝试使用最优控制中的动态规划的相关知识来进行求解。2 TSP问题分析对于这个问题，我们首先想到的是应用穷举法进行解答，但是这个方法时间和空间的复杂度很高。从表面上看，TSP 问题很简单，其实则不然。对于N 个城市的TSP,存在的可能路径为（N-1）!/2条，当N较大时，其数量是惊人的。计算每条路经都需求出N 个距离之和， 这样各种路径及其距离之和的计算量正比与N!/2.用搜索法要求就规模大的TSP是不现实的。例如使用1GFLOPs次的计算机搜索TSP所需的时间如下表所示城市数7152050100200加法量搜索时间1.8h350y 由上可知，对于这个问题采用穷举法进行解答是不现实的，这就要求我们采用其他的方法进行解答。3 其他求解TSP问题的方法*贪心法a. 所谓贪心法，就是在组合算法中，将每一步都取局部最优的求解方法。b. 下表表示用贪心法求解TSP的过程。先将各城市间的距离用行列式形式表示，主对角线上用表示。我们可以从城市C1出发，依次在每一行或列中选取元素最小的路径，且每个城市只能访问一次。 c. 按贪心法从C1出发所挑选的路径为 L2734431033不难看出，这种从局部最优原则出发的方法所得的结果的好坏，与城市间的距离的具体情况和从那个城市开始有关。例如，从C7出发时，用贪心法所得的结果为 其路线长度减为L253743731*Hopfield神经元网络法a. 全互连型神经网络求解TSP问题：设有N个城市： C= C中任意两个城市的距离 D()=现在要找到一个城市的排列使得闭合路径 为最小我们用换位矩阵来表征一条路径。对于N个城市来说，换位矩阵有N*N个元素，需要用N*N个神经元来表征。ABCDEA00001B10000C00100D01000E00010根据下面四方面的要求，即：1.换位矩阵每行只能有一个一；2.换位矩阵每列只能有一个一；3.换位矩阵中元素一之和应为N；4.所构造的函数的极小值对应于最短路径。我们构造出与TSP相对应的计算能量函数为 式中前三项与条件的1，2，3的要求对应，上式的第四项为目标 项，它的最小值就对应于最短路径长度。上式中v的数值为0或者1，是由表正换位矩阵中神经元的输出来表示的。4 TSP问题的动态规划解法主要特点：将一个问题分为若干个相互联系的阶段，每个阶段进行决策优化。在这种多阶段决策优化过程中，无论其初始状态和初始决策如何，以后的最优策略只取决于由最初策略所形成的当前策略。5 问题分析我们应用动态规划的解法来求解五个城市的TSP问题。我们采用一个矩阵A表示城市之间的距离。其中， 表示第个城市和第个城市之间的距离。这是个对称阵，而且对角线的元素均为0。 假定我们已经找到一个最短的路径，设它是先从到，则从出发沿着这条路径回到的路，必定是经过中每个城市一次，最后回到的路径中的最短的一条，也就是符合最优原理。 设表示从城市出发，通过s集合中所有城市一次且仅一次，最后回到出发城市的最短路径的长度。由f的定义知，所求最短路径长度可表示为。根据最优原理，应有一般的有：。根据以上分析，应用Matlab编程如下：clearclcD =0 146.13 356.43 286.99 349.83; 140.95 0 228.38 456.76 201.68; 364.21 233.23 0 431.53 198.65; 287.89 471.96 418.44 0 222.73; 340.68 191.78 213.68 232.22 0;n=4;c1=1;c2=n*nchoosek(n-1,n-1);c3=n*nchoosek(n-1,n-2);c4=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,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+1 layer5(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-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-1layer4(1,12)=D(4,5)+layer5(1,5);%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) 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,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,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)-1layer3(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,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,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 = 1 2 3 5 4 1 6 附录附贪心法及神经网络法的程序及相应的运行结果（均为用Matlab编写）贪心法程序：clearclcD=zeros(10,10);L_min=zeros(1,10);Path=zeros(1,10);Dmin=0;Min_D=zeros(1,10);optimal=0;for i=1:1:10 for j=1:1:10 D(i,j)=unidrnd(1000); endendfor i=1:1:10 for j=1:1:10 if i=j D(i,j)=1001; else D(i,j)=D(j,i); end endendfor n=1:1:10Path(1,1)=n;for i=1:1:9 L_min(1,i)=min(D(Path(1,i),:); for y=1:1:10 if D(Path(1,i),y)=L_min(1,i) Path(1,i+1)=y; break; end end for j=1:1:10 D(Path(1,i),j)=1001; D(j,Path(1,i)=1001; end D;endfor k=1:1:9 Dmin=Dmin+L_min(1,k);endPathDminMin_D(1,n)=Dmin;Dmin=0;endoptimal=min(Min_D) 神经网络解法程序：%人工神经网络求解TSP 200459clearclcn=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;count=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:n for y=1:1:n D(x,y)=unifrnd(0.01,1); endendfor x=1:1:n for y=1:1:n if x=y D(x,y)=0; end D(x,y)=D(y,x); endendN=n*n;uoo=-0.5*u0*log(n-1); %计算初值u00while flag=0 for i=1:1:n for j=1:1:n r=unifrnd(-0.1*u0, 0.1*u0); U(i,j)= uoo + r; %按均匀分布随机产生人工神经网络每个神经元初始输入ui，i=1,2.100 V(i,j)=0.5 * (1+tanh( U(i,j)/u0 ) ); %把ui代入 S 函数得到每个神经元的初始输出vi，i=1,2.100 end end Uini=U; %保存初态ui Vini=V; %保存初态vi for k=0:1:200 %迭代次数 for x=1:1:n %城市循环 for i=1:1:n %路径按步循环 for y=1:1:n %求微分方程中的最后一个和式：dxy*(V(y,i-1)+V(y,i+1),x,y from 1 to n. if i=1 dxy=dxy + D(x,y) * (0 + V(y,i+1); %当i=1时,不存在左顺联,所以V(y,i-1)0. elseif i=n dxy=dxy + D(x,y) * (V(y,i-1) + 0); %当i=n时,不存在顺联,所以V(y,i+1)0. else dxy=dxy + D(x,y) * (V(y,i-1) + V(y,i+1); %当1in时,用原公式. end end const=-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/dt u_xi=U(x,i)+du_xi*T; %用Eula法求u(t+T)=u(t)+du/dt * T v_xi=0.5*(1+tanh(u_xi/u0); %把u(t+T)代入S函数求t+T时刻的V(t+T) Utemp(x,i)=u_xi; %把u(t+T)存入转移数组Utemp Vtemp(x,i)=v_xi; %把V(t+T)存入转移数组Vtemp u_xi=0; %u_xi清零 dxy=0; %dxy清零 end end U=Utemp; %用t+T时刻的ui替换t时刻的ui V=Vtemp; %用t+T时刻的Vi替换t时刻的Vi end V; flag=1; for i=1:1:n for j=1:1:n if V(i,j)0.9 V(i,j)=1; else V(i,j)=2; end end end for x=1:1:n row_s=sum(V(x,:); %对V的每一行求和 column_s=sum(V(:,x); %对V的每一列求和 if column_s=1 %确保V的每一列只有一个1 flag=0; break; elseif row_s=1; %确保V的每一行只有一个1 flag=0; break; end row_s=0; column_s=0; end row_s=0; column_s=0; %for j=1:1:n %当某一列一个1也没有时,说明网络没有收敛到稳态,标志flag置零. % if City_Path(j) 0.1 % flag=0; % break; % end %end %for i=1:1:n %当某一行有多个1时,说明网络没有寻找到最优路径,标志flag置零. % for j=i:1:n-1 % if abs(City_Path(i)-City_Path(j+1)0.9 City_Path(column)=row; end end end for i=1:1:n-1 Dmin=Dmin+D(City_Path(i),City_Path(i+1); %计算最优路径的距离 end %if Dmin=50 break; endend% result: