distTSP.txt
0 6 18 4 8
7 0 17 3 7
4 4 0 4 5 
20 19 24 0 22
8 8 16 6 0 
%GATSP.m
function gatsp1()
clear;
load distTSP.txt;
distance=distTSP;
N=5;
ngen=100;
ngpool=10;
%ngen=input('# of generations to evolve = ');
%ngpool=input('# of chromosoms in the gene pool = '); % size of genepool
gpool=zeros(ngpool,N+1); % gene pool
for i=1:ngpool, % intialize gene pool 
gpool(i,:)=[1 randomize([2:N]')' 1]; 
for j=1:i-1
  while gpool(i,:)==gpool(j,:) 
       gpool(i,:)=[1 randomize([2:N]')' 1]; 
                end 
             end
          end 

  costmin=100000; 
    tourmin=zeros(1,N); 
      cost=zeros(1,ngpool); 
  increase=1;resultincrease=1;
      for i=1:ngpool, 
          cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))')));
     end
% record current best solution
[costmin,idx]=min(cost);
tourmin=gpool(idx,:);
disp([num2str(increase) 'minmum trip length = ' num2str(costmin)])

costminold2=200000;costminold1=150000;resultcost=100000;
tourminold2=zeros(1,N);
tourminold1=zeros(1,N);
resulttour=zeros(1,N);
while (abs(costminold2-costminold1) ;100)&(abs(costminold1-costmin) ;100)&(increase ;500) 

costminold2=costminold1; tourminold2=tourminold1; 
costminold1=costmin;tourminold1=tourmin; 
increase=increase+1; 
 if resultcost>costmin 
   resultcost=costmin; 
   resulttour=tourmin; 
   resultincrease=increase-1; 
         end 
 for i=1:ngpool, 
           cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))'))); 
end 
% record current best solution 
[costmin,idx]=min(cost); 
tourmin=gpool(idx,:); 
%============== 
% copy gens in th gpool according to the probility ratio  
% >1.1 copy twice 
%  >=0.9 copy once 
%  ;0.9 remove 
[csort,ridx]=sort(cost); 
% sort from small to big. 
csum=sum(csort); 
caverage=csum/ngpool; 
cprobilities=caverage./csort; 
copynumbers=0;removenumbers=0; 
for i=1:ngpool, 
    if cprobilities(i) >1.1 
             copynumbers=copynumbers+1; 
                    end 
           if cprobilities(i) <0.9 
                   removenumbers=removenumbers+1; 
                           end 
                end 
   copygpool=min(copynumbers,removenumbers); 
               for i=1:copygpool 
                  for j=ngpool:-1:2*i+2 gpool(j,:)=gpool(j-1,:); 
            end  
                   gpool(2*i+1,:)=gpool(i,:); 
          end 
                 if copygpool==0 
                       gpool(ngpool,:)=gpool(1,:); 
                  end 
%========= 
%when genaration is more than 50,or the patterns in a couple are too close,do mutation 
for i=1:ngpool/2 
        % 
sameidx=[gpool(2*i-1,:)==gpool(2*i,:)]; 
diffidx=find(sameidx==0); 
           if length(diffidx)<=2 
                gpool(2*i,:)=[1 randomize([2:12]')' 1]; 
                           end 
                               end 
%=========== 
%cross gens in couples 
           for i=1:ngpool/2 
                  [gpool(2*i-1,:),gpool(2*i,:)]=crossgens(gpool(2*i-1,:),gpool(2*i,:)); 
       end  
  
        for i=1:ngpool, 
              cost(i)=sum(diag(distance(gpool(i,:)',rshift(gpool(i,:))'))); 
       end 
% record current best solution 
  [costmin,idx]=min(cost); 
  tourmin=gpool(idx,:); 
  disp([num2str(increase) 'minmum trip length = ' num2str(costmin)])
end   
 
disp(['cost function evaluation: ' int2str(increase) ' times!'])
disp(['n:' int2str(resultincrease)])
disp(['minmum trip length = ' num2str(resultcost)])
disp('optimum tour = ')
disp(num2str(resulttour)) 
%====================================================
function B=randomize(A,rowcol)
% Usage: B=randomize(A,rowcol)
% randomize row orders or column orders of A matrix
% rowcol: if =0 or omitted, row order (default)
% if = 1, column order 

rand('state',sum(100*clock))
if nargin == 1, 
        rowcol=0;
end
         if rowcol==0,  
              [m,n]=size(A); 
              p=rand(m,1); 
              [p1,I]=sort(p); 
              B=A(I,:);
elseif rowcol==1,
          Ap=A'; 
          [m,n]=size(Ap); 
          p=rand(m,1); 
          [p1,I]=sort(p); 
          B=Ap(I,:)';
 end
%=====================================================
function y=rshift(x,dir)
% Usage: y=rshift(x,dir)
% rotate x vector to right (down) by 1 if dir = 0 (default)
% or rotate x to left (up) by 1 if dir = 1 
if nargin ;2, dir=0; end 
[m,n]=size(x); 
if m>1, 
if n == 1,  
    col=1;  
elseif n>1, 
    error('x must be a vector! break'); 
end % x is a column vectorelseif m == 1,  
if n == 1,  y=x; 
return 
  elseif n>1, 
     col=0; % x is a row vector endend 
if dir==1, % rotate left or up 
       if col==0, % row vector, rotate left 
             y = [x(2:n) x(1)]; 
       elseif col==1, 
             y = [x(2:n); x(1)]; % rotate up 
end
   elseif dir==0, % default rotate right or down  
              if col==0,  
                    y = [x(n) x(1:n-1)]; 
             elseif col==1 % column vector 
                       y = [x(n); x(1:n-1)]; 
                   end
             end
%==================================================
function [L1,L2]=crossgens(X1,X2)
% Usage:[L1,L2]=crossgens(X1,X2)
s=randomize([2:12]')';
n1=min(s(1),s(11));n2=max(s(1),s(11));
X3=X1;X4=X2;
for i=n1:n2, 
                for j=1:13, 
                     if X2(i)==X3(j), 
                          X3(j)=0; 
                             end 
                  if X1(i)==X4(j),                          X4(j)=0; 
               end 
           end
        end
   j=13;k=13;
    for i=12:-1:2, 
          if X3(i)~=0, 
               j=j-1; 
                 t=X3(j);X3(j)=X3(i);X3(i)=t; 
               end 
                    if X4(i)~=0, 
                           k=k-1; 
                      t=X4(k);X4(k)=X4(i);X4(i)=t; 
                   end
               end
           for i=n1:n2 
              X3(2+i-n1)=X2(i); 
              X4(2+i-n1)=X1(i);
           end
L1=X3;L2=X4;
%=======================