function [MinD,BestPath]=MainAneal(CityPosition,pn)
function [MinD,BestPath]=MainAneal2(CityPosition,pn)
%此题以中国31省会城市的最短旅行路径为例，给出TSP问题的模拟退火程序
%CityPosition_31=[1304 2312;3639 1315;4177 2244;3712 1399;3488 1535;3326 1556;...
% 3238 1229;4196 1044;4312 790;4386 570;3007 1970;2562 1756;...
31省新增24例输入% 2788 1491;2381 1676;1332 695;3715 1678;3918 2179;4061 2370;...
% 3780 2212;3676 2578;4029 2838;4263 2931;3429 1908;3507 2376;...
% 3394 2643;3439 3201;2935 3240;3140 3550;2545 2357;2778 2826;2370 2975];
%T0=clock
global path p2 D;
[m,n]=size(CityPosition);
%生成初始解空间，这样可以比逐步分配空间运行快一些
TracePath=zeros(1e3,m);
Distance=inf*zeros(1,1e3);
D = sqrt((CityPosition( :, ones(1,m)) - CityPosition( :, ones(1,m))').^2 +...
(CityPosition( : ,2*ones(1,m)) - CityPosition( :,2*ones(1,m))').^2 );
%将城市的坐标矩阵转换为邻接矩阵（城市间距离矩阵）
for i=1:pn
path(i,:)=randperm(m);%构造一个初始可行解
end
t=zeros(1,pn);
p2=zeros(1,m);
iter_max=100;%input('请输入固定温度下最大迭代次数iter_max=' );
m_max=5;%input('请输入固定温度下目标函数值允许的最大连续未改进次数m_nax=' ) ;
%如果考虑到降温初期新解被吸收概率较大，容易陷入局部最优
%而随着降温的进行新解被吸收的概率逐渐减少，又难以跳出局限
%人为的使初期 iter_max,m_max 较小，然后使之随温度降低而逐步增大，可能
%会收到到比较好的效果
T=1e5;
N=1;
tau=1e-5;%input('请输入最低温度tau=' );
%nn=ceil(log10(tau/T)/log10(0.9));
while T>=tau%&m_num<m_max
iter_num=1;%某固定温度下迭代计数器
m_num=1;%某固定温度下目标函数值连续未改进次数计算器
%iter_max=100;
%m_max=10;%ceil(10+0.5*nn-0.3*N);
while m_num<m_max&iter_num<iter_max
%MRRTT(Metropolis, Rosenbluth, Rosenbluth, Teller, Teller)过程:
%用任意启发式算法在path的领域N(path)中出新的更优解
for i=1:pn
Len1(i)=sum([D(path(i,1:m-1)+m*(path(i,2:m)-1)) D(path(i,m)+m*(path(i,1)-1))]);
%计算一次行遍所有城市的总路程
[path2(i,: )]=ChangePath2(path(i,: ),m);%更新路线
Len2(i)=sum([D(path2(i,1:m-1)+m*(path2(i,2:m)-1)) D(path2(i,m)+m*(path2(i,1)-1))]);
end
%Len1
%Len2
%if Len2-Len1<0|exp((Len1-Len2)/(T))>rand
R=rand(1,pn);
%Len2-Len1<t|exp((Len1-Len2)/(T))>R
if find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0)
path(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0), : )=path2(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0), : );
Len1(find((Len2-Len1<t|exp((Len1-Len2)/(T))>R)~=0))=Len2(find((Len2-Len1
<t|exp((Len1-Len2)/(T))>R)~=0));
[TempMinD,TempIndex]=min(Len1);
%TempMinD
TracePath(N,: )=path(TempIndex,: );
Distance(N,: )=TempMinD;
N=N+1;
%T=T*0.9
m_num=0;
else
m_num=m_num+1;
end
iter_num=iter_num+1;
end
T=T*0.9
%m_num,iter_num,N
end
[MinD,Index]=min(Distance);
BestPath=TracePath(Index,: );
disp(MinD)
%T1=clock
%更新路线子程序
function [p2]=ChangePath2(p1,CityNum)
global p2;
while(1)
R=unidrnd(CityNum,1,2);
if abs(R(1)-R(2))>1
break;
end
end
R=unidrnd(CityNum,1,2);
I=R(1);J=R(2);
%len1=D(p(I),p(J))+D(p(I+1),p(J+1));
%len2=D(p(I),p(I+1))+D(p(J),p(J+1));
if I<J
p2(1:I)=p1(1:I);
p2(I+1:J)=p1(J:-1:I+1);
p2(J+1:CityNum)=p1(J+1:CityNum);
else
p2(1:J)=p1(1:J);
p2(J+1:I)=p1(I:-1:J+1);
p2(I+1:CityNum)=p1(I+1:CityNum);
end
版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。
发表评论