TSP问题的动态规划解法.docx
《TSP问题的动态规划解法.docx》由会员分享,可在线阅读,更多相关《TSP问题的动态规划解法.docx(14页珍藏版)》请在冰豆网上搜索。
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));%当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;
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: