matlabtsp.m

蚁群混合算法源代码 采用MATLAB

%模拟退火算法m源程序
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;...
%                 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