(完整word版)TSP问题的动态规划解法

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 所需的时间如下表所示 城市数

7

15

20

50

100

200

加法量 3105.2? 11105.6? 18102.1? 64105.1? 157105? 37410

搜索时间

s 5105.2-?

1.8h

350y

y 48105? y 14210

y 35810

由上可知，对于这个问题采用穷举法进行解答是不现实的，这就要求我们采用其他的方法进行解答。 3． 其他求解TSP 问题的方法

*贪心法

a. 所谓贪心法，就是在组合算法中，将每一步都取局部最优的求解方法。

b. 下表表示用贪心法求解TSP 的过程。先将各城市间的距离用行列式形式表示，主对角线上用∞表示。我们可以从城市C1出发，依次在每一行或列中选取元素最小的路径，且每个城市只能访问一次。

c. 按贪心法从C1出发所挑选的路径为 15432671C C C C C C C C →→→→→→→ L ＝2＋7＋3＋4＋4＋3＋10＝33

不难看出，这种从局部最优原则出发的方法所得的结果的好坏，与城市间的距离的具体情况和从那个城市开始有关。 例如，从C7出发时，用贪心法所得的结果为 76235417C C C C C C C C →→→→→→→ 其路线长度减为

L ＝2＋5＋3＋7＋4＋3＋7＝31

*Hopfield 神经元网络法

a. 全互连型神经网络求解TSP 问题： 设有N 个城市：

C={n C C C .....................21} C 中任意两个城市的距离 D(j i

c c )=ij d

现在要找到一个城市的排列{})(.....)2()1(N c c c n n n 使

得闭合路径

[]∑+n n n

i c i c d mod )1()

( 为最小

我们用换位矩阵来表征一条路径。对于N 个城市来说，换位矩阵有N*N 个元素，需要用N*N 个神经元来表征。

根据下面四方面的要求，即：1.换位矩阵每行只能有一个一； 2.换位矩阵每列只能有一个一；3.换位矩阵中元素一之和应为N ；4.所构造的函数的极小值对应于最短路径。 我们构造出与TSP 相对应的计算能量函数为

∑∑∑∑∑∑∑∑∑∑∑===+===-=+==-=+=++??

?

??-++=N x N y N i i y i y xi xy N x N i xi N i N x N i y yi xi N x N i N i j xj xi v v v d d N v c v v b v v a E 1111_,1,2

1111111111)

(2222 式中前三项与条件的1，2，3的要求对应，上式的第四项为目标 项，它的最小值就对应于最短路径长度。上式中v 的数值为0或者1，是由表正换位矩阵中神经元的输出来表示的。

4． TSP 问题的动态规划解法

主要特点：将一个问题分为若干个相互联系的阶段，每个阶段进行决策优化。在这种多阶段决策优化过程中，无论其初始状态和初始决策如何，以后的最优策略只取决于由最初策略所形成的当前策略。

5． 问题分析

我们应用动态规划的解法来求解五个城市的TSP 问题。我们采用一个矩阵A 表示城市之间的距离。

???

???

???

???

???

?????????

????????=0000000000910

810

710

610

510

410

310

210

110

910897969594939291981089786858483828187107978675747372717610696867564636261651059585756453525154104948474645342414310393837363534231321029282726252423121101918171615141312

a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a a A 其中，ij a ，109.......21,0=i 109.......21,0=j 表示第i 个城市和第j 个城市之间的距离。这是个对称阵，而且对角线的元素均为0。

假定我们已经找到一个最短的路径，设它是先从1c 到k c ，则从k c 出发沿着这条路径回到1c 的路，必定是经过{}k c c C ,1-中每个城市一次，最后回到1c 的路径中的最短的一条，也就是符合最优原理。

设()S c f ,1表示从城市1c 出发，通过s 集合中所有城市一次且仅一次，最后回到出发城市1c 的最短路径的长度。由f 的定义知，所求最短路径长度可表示为{}()11,c C c f -。

根据最优原理，应有

{}(){}(){}k k k k c c C c f d c C c f ,,,1110

211min -+=-≤≤

一般的有：(){}(){}j j ij S

c i c C c f

d S c f i

-+=?,min

,。

根据以上分析，应用Matlab编程如下：

clear

clc

D =[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);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,c4);layer5=zeros(1,c5);

path1=0;path2=zeros(1,c2);path3=zeros(1,c3);path4=zero s(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-1

layer4(1,2)=D(2,3)+layer5(1,3); %2-3-1

layer4(1,3)=D(4,2)+layer5(1,2); %4-2-1

layer4(1,4)=D(2,4)+layer5(1,3); %2-4-1

layer4(1,5)=D(5,2)+layer5(1,2); %5-2-1

layer4(1,6)=D(2,5)+layer5(1,5); %2-5-1

layer4(1,7)=D(4,3)+layer5(1,3); %4-3-1

layer4(1,8)=D(3,4)+layer5(1,4); %3-4-1

layer4(1,9)=D(5,3)+layer5(1,3); %5-3-1

layer4(1,10)=D(3,5)+layer5(1,5);%3-5-1

layer4(1,11)=D(5,4)+layer5(1,4);%5-4-1

layer4(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)-1

layer3(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)-1

layer3(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)-1 layer3(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)-1 layer3(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)-1 layer3(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)-1 layer3(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)-1 layer3(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)-1 layer3(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)-1

layer3(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)-1

layer3(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)-1

layer3(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)-1

layer2(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)-1

layer2(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)-1

layer2(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)-1

layer2(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)-1 layer1(1,1)=dxy

path1=p;

path2=path2;

path3=path3;

optimal_path=[1 2 3 5 4 1]

运行结果：

layer1 = 1.0933e+003

optimal_path = 1 2 3 5 4 1

6．附录

附贪心法及神经网络法的程序及相应的运行结果（均为用Matlab编写）

贪心法程序：

clear

clc

D=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);

end

end

for i=1:1:10

for j=1:1:10

if i==j

D(i,j)=1001;

else

D(i,j)=D(j,i);

end

end

for n=1:1:10

Path(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;

end

for k=1:1:9

Dmin=Dmin+L_min(1,k);

end

Path

Min_D(1,n)=Dmin;

Dmin=0;

end

optimal=min(Min_D)

神经网络解法程序：

%人工神经网络求解TSP 2004－5－9

clear

clc

n=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);

end

end

for x=1:1:n

for y=1:1:n

if x==y

D(x,y)=0;

end

D(x,y)=D(y,x);

end

end

N=n*n;

uoo=-0.5*u0*log(n-1); %计算初值u00

while 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)); %当1

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.1

V(i,j)=0;

elseif 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.1

% flag=0;

% break;

% end

% end

% if flag==0

% break;

% end

%end

if flag==1 %当标志flag=1时,说明网络收敛到稳态,同时找到最优路径.

for column=1:1:n %确定每一列1所在行的位置

for row=1:1:n

if V(row,column)>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)); %计