TSP问题的动态规划解法可编辑修改word版.docx
《TSP问题的动态规划解法可编辑修改word版.docx》由会员分享，可在线阅读，更多相关《TSP问题的动态规划解法可编辑修改word版.docx（19页珍藏版）》请在冰点文库上搜索。
TSP问题的动态规划解法可编辑修改word版
TSP问题的动态规划解法
第十七组：
3103038028郑少斌
3103038029王瑞锋
3103038035江飞鸿
3103038043韩鑫
3103055004唐万强
1.TSP问题简介
旅行商问题（TravelingSalesmanProblem,简称TSP,亦称为货单郎问题）可以描述为：
对于N个城市，它们之间的距离已知，有一旅行商要从某一城市走遍所有的城市，且每一城市只能经过一次，最后回到出发的城市，问如何选择路线可使他走过的路径最短。
这是一个典型的组合优化问题。
它有很强的现实意义，可以应用于交通运输，物资调配，旅游线路设置。
对于了解某个国家地理分布也有一定的现实意义。
这个问题的解法有很多种，在这里我们尝试使用最优控制中的动态规划的相关知识来进行求解。
2.TSP问题分析
对于这个问题，我们首先想到的是应用穷举法进行解答，但是这个方法时间和空间的复杂度很高。
从表面上看，TSP问题很简单，其实则不然。
对于N个城市的TSP,存在的可能路径为（N-1）!
/2条，当N较大时，其数量是惊人的。
计算每条路经都需求出N个距离之和，这样各种路径及其距离之和的计算量正比与N!
/2.用搜索法要求就规模大的TSP是不现实的。
城市数
7
15
20
50
100
200
加法量
2.5⨯103
6.5⨯1011
1.2⨯1018
1.5⨯1064
5⨯10157
10374
搜索时
间
2.5⨯10-5s
1.8h
350y
5⨯1048y
10142y
10358y
例如使用1GFLOPs次的计算机搜索TSP所需的时间如下表所示
由上可知，对于这个问题采用穷举法进行解答是不现实的，这就要求我们采用其他的方法进行解答。
3.其他求解TSP问题的方法
*贪心法
a.所谓贪心法，就是在组合算法中，将每一步都取局部最优的求解方法。
b.
下表表示用贪心法求解TSP的过程。
先将各城市间的距离用行列式形式表示，主对角线上用∞表示。
我们可以从城市C1出发，依次在每一行或列中选取元素最小的路径，且每个城市只能访问一次。
c.按贪心法从C1出发所挑选的路径为
C1→C7→C6→C2→C3→C4→C5→C1
L＝2＋7＋3＋4＋4＋3＋10＝33
不难看出，这种从局部最优原则出发的方法所得的结果的好坏，与城市间的距离的具体情况和从那个城市开始有关。
例如，从C7出发时，用贪心法所得的结果为
C7→C1→C4→C5→C3→C2→C6→C7
其路线长度减为
L＝2＋5＋3＋7＋4＋3＋7＝31
*Hopfield神经元网络法
a.全互连型神经网络求解TSP问题：
设有N个城市：
C={C1C2......
.......
........Cn}
C中任意两个城市的距离D（cicj）=dij
现在要找到一个城市的排列{cn
（1）
cn
（2）
.....
cn（N）}使
得闭合路径
∑d[cn（i）
cn（i+1）modn]
为最小
我们用换位矩阵来表征一条路径。
对于N个城市来说，换位矩阵有N*N个元素，需要用N*N个神经元来表征。
A
B
C
D
E
A
0
0
0
0
1
B
1
0
0
0
0
C
0
0
1
0
0
D
0
1
0
0
0
E
0
0
0
1
0
根据下面四方面的要求，即：
1.换位矩阵每行只能有一个一；2.换位矩阵每列只能有一个一；3.换位矩阵中元素一之和应为N；4.所构造的函数的极小值对应于最短路径。
我们构造出与TSP相对应的计算能量函数为
aNN-1N
bNN-1N
c⎛NN
⎫2dNNN
E=2∑∑∑vxivxj+2∑∑∑vxivyi+2ç∑∑vxi-N⎪+2∑∑∑dxyvxi（vy,i+1+vy,i_1）
x=1i=1j=i+1
i=1x=1y=i+1
⎝x=1i=1
⎭x=1y=1i=1
式中前三项与条件的1，2，3的要求对应，上式的第四项为目标项，它的最小值就对应于最短路径长度。
上式中v的数值为0或者1，是由表正换位矩阵中神经元的输出来表示的。
4.TSP问题的动态规划解法
主要特点：
将一个问题分为若干个相互联系的阶段，每个阶段进行决策优化。
在这种多阶段决策优化过程中，无论其初始状态和初始决策如何，以后的最优策略只取决于由最初策略所形成的当前策略。
5.问题分析
我们应用动态规划的解法来求解五个城市的TSP问题。
我们采用一个矩阵A表示城市之间的距离。
⎡0a12
a13
a14
a15
a16
a17
a18
a19
a110⎤
⎢a0aa
aaa
aaa⎥
⎢12
2324
2526
272829
210⎥
⎢a13
a230
a34
a35a36
a37a38a39
a310⎥
⎢aaa0
aaa
aaa⎥
⎢14
2434
4546
474849
410⎥
⎢
A=⎢a15
⎢a16
a25a35
a26a36
a45a46
0a56
a560
a57a58a59
a67a68a69
a510⎥
⎥
a610⎥
⎢aaaa
aa0
aaa⎥
⎢17273747
5767
7879710⎥
⎢a18a28a38a48
a58a68
a780
a89a810⎥
⎢aaaaaaaa0a⎥
⎢19
293949
5969
7989
910⎥
⎢⎣a110
a210a310a410
a510a610
a710a810
a9100⎥⎦
其中，
aij，i=0,12
.......
910
j=0,12
.......
910
表示第i个城市和第j个城市之间的距离。
这是个对称阵，
而且对角线的元素均为0。
假定我们已经找到一个最短的路径，设它是先从c1到ck，
则从c
出发沿着这条路径回到c的路，必定是经过C-{cc}
k11,k
中每个城市一次，最后回到c1的路径中的最短的一条，也就
是符合最优原理。
设f（c1,S）表示从城市c1出发，通过s集合中所有城市一次且仅一次，最后回到出发城市c1的最短路径的长度。
由f的定义知，所求最短路径长度可表示为f（c1,C-{c1}）。
根据最优原理，应有
f（c1,C-{c1}）=min{d1k+f（ck,C-{c1,ck}）}
2≤k≤10
一般的有：
f（c,S）=min{d+f（c,C-{c}）}。
ici∉Sijjj
根据以上分析，应用Matlab编程如下：
clear
clc
D=[0
146.13
356.43
286.99349.83;
140.95
0
228.38
456.76201.68;
364.21
233.23
0
431.53198.65;
287.89
471.96
418.44
0222.73;
340.68
191.78
213.68
232.220];
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）;
%###############################第五层
fori=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-（23）-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-（23）-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-（24）-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-（24）-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-（25）-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-（25）-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-（34）-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-（34）-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-（35）-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-（35）-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-（45）-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-（45）-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-（345）-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-（245）-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-（235）-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-（234）-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-（2345）-1
layer1（1,1）=dxypath1=p;path2=path2;path3=path3;
optimal_path=[123541]
运行结果：
layer1=1.0933e+003
optimal_path=123541
6.附录
附贪心法及神经网络法的程序及相应的运行结果（均为用
Matlab编写）贪心法程序：
clear
clcD=zeros（10,10）;L_min=zeros（1,10）;Path=zeros（1,10）;Dmin=0;Min_D=zeros（1,10）;optimal=0;
fori=1:
1:
10
forj=1:
1:
10
D（i,j）=unidrnd（1000）;
end
end
fori=1:
1:
10
forj=1:
1:
10
ifi==j
D（i,j）=1001;
else
D（i,j）=D（j,i）;
end
end
end
forn=1:
1:
10Path（1,1）=n;fori=1:
1:
9
L_min（1,i）=min（D（Path（1,i）,:
））;fory=1:
1:
10
ifD（Path（1,i）,y）==L_min（1,i）Path（1,i+1）=y;
break;
end
end
forj=1:
1:
10
D（Path（1,i）,j）=1001;
D（j,Path（1,i））=1001;
endD;
end
fork=1:
1:
9
Dmin=Dmin+L_min（1,k）;
endPath
DminMin_D（1,n）=Dmin;Dmin=0;
endoptimal=min（Min_D）
神经网络解法程序：
%人工神经网络求解TSP2004－5－9
clearclc
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）;
%初始化数组forx=1:
1:
n
fory=1:
1:
n
D（x,y）=unifrnd（0.01,1）;
end
end
forx=1:
1:
nfory=1:
1:
n
ifx==y
D（x,y）=0;
endD（x,y）=D（y,x）;
endendN=n*n;
uoo=-0.5*u0*log（n-1）;%计算初值u00
whileflag==0fori=1:
1:
n
forj=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
endend
Uini=U;%保存初态ui
Vini=V;%保存初态vi
fork=0:
1:
200%迭代次数
forx=1:
1:
n%城市循环
fori=1:
1:
n%路径按步循环
fory=1:
1:
n%求微分方程中的最后一个和式：
∑dxy*（V（y,i-1）+V（y,i+1））,x,yfrom1ton.ifi==1
dxy=dxy+D（x,y）*（0+V（y,i+1））;
%当i=1时,不存在左顺联,所以V（y,i-1）＝0.
elseifi==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））;%当1end
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
endV;
flag=1;
fori=1:
1:
n
forj=1:
1:
n
ifV（i,j）<0.1
V（i,j）=0;
elseifV（i,j）>0.9
V（i,j）=1;
else
V（i,j）=2;
end
end
end
forx=1:
1:
n
row_s=sum（V（x,:
））;%对V的每一行求和column_s=sum（V（:
x））;%对V的每一列求和
ifcolumn_s~=1%确保V的每一列只有一个1
flag=0;break;
elseifrow_s~=1;%确保V的每一
行只有一个1
flag=0;break;
endrow_s=0;column_s=0;
endrow_s=0;column_s=0;
%forj=1:
1:
n%当某一列一个1
也没有时,说明网络没有收敛到稳态,标志flag置零.
%ifCity_Path（j）<0.1
%flag=0;
%break;
%end
%end
%fori=1:
1:
n%当某一行有多个1时,说明网络没有寻找到最优路径,标志flag置零.
%forj=i:
1:
n-1
%ifabs（City_Path（i）-City_Path（j+1））<0.1
%flag=0;
%break;
%end
%end
%ifflag==0
%break;
%end
%end
ifflag==1%当标志flag=1时,说明网络收敛到稳态,同时找到最优路径.
forcolumn=1:
1:
n%确定每一列1所在行的位置
forrow=1:
1:
n
ifV（row,column）>0.9City_Path（column）=row;
end
endend
fori=1:
1:
n-1
Dmin=Dmin+D（City_Path（i）,City_Path（i+1））;%计算最优路径的距离
end
%ifDmin<=20.0
Uini=Uini;%显示初始输入Uini
Vini=Vini;%显示初始输出ViniU=U;%显示稳态输入U
V=V%显示稳态输出VCity_Path=City_Path%显示最优路径Dmin=Dmin%显示最短距离
%break;
%elseflag=0;
%end
endDmin=0;
count=count+1%寻最优的次数
ifcount>=50break;
endend
%result:
升级会员 