%this program is written by 劉學智. Finished time is 05.1.23 16:03 
%utilizing it solving TSP problem by simulating stealing algorithm
%[fval,route]=sa_tsp(d,10,0.1,.87)
%d=[0	2	1	2	0	0	1	0	1	2	1	1	1	1
%2	0	1	4	1	0	1	1	1	3	1	0	2	1
%1	1	0	1	0	0	0	3	1	1	0	2	2	1
%2	4	1	0	1	1	2	1	0	2	1	0	1	1
%0	1	0	1	0	2	0	1	1	1	0	1	1	2
%0	0	0	1	2	0	1	2	1	1	1	2	1	2
%1	1	0	2	0	1	0	1	1	1	0	2	2	1
%0	1	3	1	1	2	1	0	1	2	1	4	2	2
%1	1	1	0	1	1	1	1	0	1	1	1	3	1
%2	3	1	2	1	1	1	2	1	0	1	0	0	3
%1	1	0	1	0	1	0	1	1	1	0	3	1	1
%1	0	2	0	1	2	2	4	1	0	3	0	1	0
%1	2	2	1	1	1	2	2	3	0	1	1	0	4
%1	1	1	1	2	2	1	2	1	3	1	0	4	0];

%the result is fval=2; route=14	 9	4	13	10	12	2	6	3	11	7	5	1	8

function [fval,route]=sa_tsp(d,t0,tf,alpha)
%d is the distance matrix;t0,tf is the initial and finil temperature;
%alpha is controling temperature coeffient
d=[0	2	1	2	0	0	1	0	1	2	1	1	1	1
2	0	1	4	1	0	1	1	1	3	1	0	2	1
1	1	0	1	0	0	0	3	1	1	0	2	2	1
2	4	1	0	1	1	2	1	0	2	1	0	1	1
0	1	0	1	0	2	0	1	1	1	0	1	1	2
0	0	0	1	2	0	1	2	1	1	1	2	1	2
1	1	0	2	0	1	0	1	1	1	0	2	2	1
0	1	3	1	1	2	1	0	1	2	1	4	2	2
1	1	1	0	1	1	1	1	0	1	1	1	3	1
2	3	1	2	1	1	1	2	1	0	1	0	0	3
1	1	0	1	0	1	0	1	1	1	0	3	1	1
1	0	2	0	1	2	2	4	1	0	3	0	1	0
1	2	2	1	1	1	2	2	3	0	1	1	0	4
1	1	1	1	2	2	1	2	1	3	1	0	4	0];
t0=10;
n=length(d);%the number of cities
L=100*n;%the length of Markov chain
route=randperm(n);%the initial traveling route
fval=value(route,d);%the initial goal value
t=t0;
tf=0.1;
alpha=.87;
tic
while t>tf
    for i=1:L
        [fval_after,route_after]=exchange(route,d);
        if fval_after<fval
           route=route_after;
           fval=fval_after;
        elseif exp((fval-fval_after)/t)>rand
               route=route_after;
               fval=fval_after;
        else route=route;
             fval=fval;
        end
        route
        fval
    end
    t=alpha*t;
end
toc
%----------------------------------------------------------------
function fval=value(route,d)%used for reckoning the goal value of the selected traveling route
n=length(d);
fval=0;
for i=1:n-1
    fval=fval+d(route(i),route(i+1));
end
%fval=fval+d(route(n),route(1));% if'%'is omited,it computes a circle,else
%a chain------------------------------------------------------------------
function [fval_after,route_after]=exchange(route,d)
%changing traveling route by inversing the sequence between two selected 2 locations 
n=length(d);
location1=ceil(n*rand);
location2=ceil(n*rand);%the location of two exchanged number
loc1=min(location1,location2);loc2=max(location1,location2);
middle_route=fliplr(route(loc1:loc2));%the part route which has been exchanged
route_after=[route(1:loc1-1) middle_route route(loc2+1:n)];%the after traveling route
fval_after=value(route_after,d);
%----------------------------------------------------------------