%模擬退火算法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