function tsp()
clc
N=30;
a=0.6;
t0=10;
tf=0.001;
xx(1)=87;yy(1)=7;
xx(2)=91;yy(2)=38;
xx(3)=83;yy(3)=46;
xx(4)=71;yy(4)=44;
xx(5)=64;yy(5)=60;
xx(6)=68;yy(6)=69;
xx(7)=83;yy(7)=69;
xx(8)=87;yy(8)=76;
xx(9)=74;yy(9)=78;
xx(10)=71;yy(10)=71;
xx(11)=58;yy(11)=69;
xx(12)=54;yy(12)=62;
xx(13)=51;yy(13)=67;
xx(14)=37;yy(14)=84;
xx(15)=41;yy(15)=94;
xx(16)=2;yy(16)=99;
xx(17)=7;yy(17)=64;
xx(18)=22;yy(18)=60;
xx(19)=25;yy(19)=62;
xx(20)=18;yy(20)=54;
xx(21)=4;yy(21)=50;
xx(22)=13;yy(22)=40;
xx(23)=18;yy(23)=40;
xx(24)=24;yy(24)=42;
xx(25)=25;yy(25)=38;
xx(26)=41;yy(26)=26;
xx(27)=45;yy(27)=21;
xx(28)=44;yy(28)=35;
xx(29)=58;yy(29)=35;
xx(30)=62;yy(30)=32;
figure(1);
plot(xx,yy,'b+');text(87,7,'\rightarrow1');
xx1=xx;xx1(31)=xx(1);
yy1=yy;yy1(31)=yy(1);
figure(2);
plot(xx1,yy1);
for i=1:N
    for j=1:N
        if i==j
            continue;
        end
        d(i,j)=sqrt((xx(i)-xx(j)).^2+(yy(i)-yy(j)).^2);
    end
end
gen=1;
time1=datestr(now)
while gen<=10
    [f,T]=trp(a,d,t0,tf);
    opf(gen)=f;
    path(gen,:)=T;
    gen=gen+1;
end 
time2=datestr(now)
disp('最小路徑和')
min(opf)
disp('最優路線')
T
for i=1:N
    xx2(i)=xx(T(i));
    yy2(i)=yy(T(i));
end
xx2(31)=xx(1);yy2(31)=yy(1);
figure(3);plot(xx2,yy2);
function [f,T]=trp(a,d,t0,tf) %f為目標函數最優值，T為最優路線，d為距離矩陣
                            %t0為初始溫度，tf為結束溫度
   [m,n]=size(d);
   L=1000*n;
   t=t0;
   pi0=1:n;
   min_f=0;
   for k=1:(n-1)
        min_f=min_f+d(pi0(k),pi0(k+1));
    end
    min_f=min_f+d(pi0(n),pi0(1));
    p_min=pi0;
    while t>tf
        for k=1:L
            kk=rand;
            [d_f,pi_1]=exchange_2(pi0,d);
            r_r=rand;
            if d_f<0
                pi0=pi_1;
            elseif exp(-d_f/t)>r_r
                pi0=pi_1;
            else
                pi0=pi0;
            end
        end
        f_temp=0;
        for k=1:(n-1)
            f_temp=f_temp+d(pi0(k),pi0(k+1));
        end
        f_temp=f_temp+d(pi0(n),pi0(1));
        if min_f>f_temp
              min_f=f_temp;
              p_min=pi0;
          end
          t=a*t;
      end
      f=min_f;
      T=p_min;
      %下面的函數產生新解
      function [d_f,pi_r]=exchange_2(pi0,d)
      [m,n]=size(d);
      clear m;
      u=rand;
      u=u*(n-2);
      u=round(u);
      
      if u<2
          u=2;
      end
      
      if u>n-2
          u=n-2;
      end
      v=rand;
      v=v*(n-u+1);
      v=round(v);
      
      if v<1
          v=1;
      end
      v=u+v;
      
      if v>n
          v=n;
      end
      pi_1(u)=pi0(v);
      pi_1(v)=pi0(u);
      
      if u>1
          for k=1:(u-1)
              pi_1(k)=pi0(k);
          end
      end
      
      if v>(u+1)
          for k=1:(v-u-1)
              pi_1(u+k)=pi0(v-k);
          end
      end
      
      if v<n
          for k=(v+1):n
              pi_1(k)=pi0(k);
          end
      end
      d_f=0;
      
      if v<n
          d_f=d(pi0(u-1),pi0(v))+d(pi0(u),pi0(v+1));
          for k=(u+1):n
              d_f=d_f+d(pi0(k),pi0(k-1));
          end
          d_f=d_f-d(pi0(u-1),pi0(u))-d(pi0(v),pi0(v+1));
          for k=(u+1):n
              d_f=d_f-d(pi0(k-1),pi0(k));
          end
      else
          d_f=d(pi0(u-1),pi0(v))+d(pi0(u),pi0(1))-d(pi0(u-1),pi0(u))-d(pi0(v),pi0(1));
          for k=(u+1):n
              d_f=d_f+d(pi0(k),pi0(k-1));
          end
          for k=(u+1):n
              d_f=d_f-d(pi0(k-1),pi0(k));
          end
      end
      pi_r=pi_1;