完整word版TSP问题的动态规划解法.docx
《完整word版TSP问题的动态规划解法.docx》由会员分享，可在线阅读，更多相关《完整word版TSP问题的动态规划解法.docx（16页珍藏版）》请在冰豆网上搜索。
完整word版TSP问题的动态规划解法
TSP问题的动态规划解法
第十七组：
3103038028郑少斌
3103038029王瑞锋
3103038035江飞鸿
3103038043韩鑫
3103055004唐万强
1．TSP问题简介
旅行商问题（TravelingSalesmanProblem,简称TSP,亦称为货单郎问题）可以描述为：
对于N个城市，它们之间的距离已知，有一旅行商要从某一城市走遍所有的城市，且每一城市只能经过一次，最后回到出发的城市，问如何选择路线可使他走过的路径最短。
这是一个典型的组合优化问题。
它有很强的现实意义，可以应用于交通运输，物资调配，旅游线路设置。
对于了解某个国家地理分布也有一定的现实意义。
这个问题的解法有很多种，在这里我们尝试使用最优控制中的动态规划的相关知识来进行求解。
2．TSP问题分析
对于这个问题，我们首先想到的是应用穷举法进行解答，但是这个方法时间和空间的复杂度很高。
从表面上看，TSP问题很简单，其实则不然。
对于N个城市的TSP,存在的可能路径为（N-1）!
/2条，当N较大时，其数量是惊人的。
计算每条路经都需求出N个距离之和，这样各种路径及其距离之和的计算量正比与N!
/2.用搜索法要求就规模大的TSP是不现实的。
例如使用1GFLOPs次的计算机搜索TSP所需的时间如下表所示
城市数
7
15
20
50
100
200
加法量
搜索时间
1.8h
350y
由上可知，对于这个问题采用穷举法进行解答是不现实的，这就要求我们采用其他的方法进行解答。
3．其他求解TSP问题的方法
*贪心法
a.所谓贪心法，就是在组合算法中，将每一步都取局部最优的求解方法。
b.下表表示用贪心法求解TSP的过程。
先将各城市间的距离用行列式形式表示，主对角线上用∞表示。
我们可以从城市C1出发，依次在每一行或列中选取元素最小的路径，且每个城市只能访问一次。
c.按贪心法从C1出发所挑选的路径为
L＝2＋7＋3＋4＋4＋3＋10＝33
不难看出，这种从局部最优原则出发的方法所得的结果的好坏，与城市间的距离的具体情况和从那个城市开始有关。
例如，从C7出发时，用贪心法所得的结果为
其路线长度减为
L＝2＋5＋3＋7＋4＋3＋7＝31
*Hopfield神经元网络法
a.全互连型神经网络求解TSP问题：
设有N个城市：
C={
}
C中任意两个城市的距离D（
）=
现在要找到一个城市的排列
使得闭合路径
为最小
我们用换位矩阵来表征一条路径。
对于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相对应的计算能量函数为
式中前三项与条件的1，2，3的要求对应，上式的第四项为目标项，它的最小值就对应于最短路径长度。
上式中v的数值为0或者1，是由表正换位矩阵中神经元的输出来表示的。
4．TSP问题的动态规划解法
主要特点：
将一个问题分为若干个相互联系的阶段，每个阶段进行决策优化。
在这种多阶段决策优化过程中，无论其初始状态和初始决策如何，以后的最优策略只取决于由最初策略所形成的当前策略。
5．问题分析
我们应用动态规划的解法来求解五个城市的TSP问题。
我们采用一个矩阵A表示城市之间的距离。
其中，
，
表示第
个城市和第
个城市之间的距离。
这是个对称阵，而且对角线的元素均为0。
假定我们已经找到一个最短的路径，设它是先从
到
，则从
出发沿着这条路径回到
的路，必定是经过
中每个城市一次，最后回到
的路径中的最短的一条，也就是符合最优原理。
设
表示从城市
出发，通过s集合中所有城市一次且仅一次，最后回到出发城市
的最短路径的长度。
由f的定义知，所求最短路径长度可表示为
。
根据最优原理，应有
一般的有：
。
根据以上分析，应用Matlab编程如下：
clear
clc
D=[0146.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）;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）=dxy
path1=p;
path2=path2;
path3=path3;
optimal_path=[123541]
运行结果：
layer1=1.0933e+003
optimal_path=123541
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;
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:
10
Path（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;
end
D;
end
fork=1:
1:
9
Dmin=Dmin+L_min（1,k）;
end
Path
Dmin
Min_D（1,n）=Dmin;
Dmin=0;
end
optimal=min（Min_D）
神经网络解法程序：
%人工神经网络求解TSP2004－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）;%初始化数组
forx=1:
1:
n
fory=1:
1:
n
D（x,y）=unifrnd（0.01,1）;
end
end
forx=1:
1:
n
fory=1:
1:
n
ifx==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
whileflag==0
fori=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
end
end
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
end
V;
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;
end
row_s=0;
column_s=0;
end
row_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.9
City_Path（column）=row;
end
end
end
fori=1:
1:
n-1
Dmin=Dmin+D（City_Path（i）,City_Path（i+1））;%计算最优路径的距离
end
%ifDmin<=20.0
Uini=Uini;%显示初始输入Uini
Vini=Vini;%显示初始输出Vini
U=U;%显示稳态输入U
V=V%显示稳态输出V
City_Path=City_Path%显示最优路径
Dmin=Dmin%显示最短距离
%break;
%else
flag=0;
%end
end
Dmin=0;
count=count+1%寻最优的次数
ifcount>=50
break;
end
end
%result:
升级会员 