)遺傳算法解決 TSP問題(附MATLAB源程序)2007-04-25 16:26已知n個(gè)城市之間的相互距離，現(xiàn)有一個(gè)推銷員必須遍訪這n個(gè)城市，并且每個(gè)城市
只能訪問一次，最后又必須返回出發(fā)城市。如何安排他對(duì)這些城市的訪問次序，可使其
旅行路線的總長度最短？
     用圖論的術(shù)語來說，假設(shè)有一個(gè)圖g=(v,e)，其中v是頂點(diǎn)集，e是邊集，設(shè)d=(dij)
是由頂點(diǎn)i和頂點(diǎn)j之間的距離所組成的距離矩陣，旅行商問題就是求出一條通過所有頂
點(diǎn)且每個(gè)頂點(diǎn)只通過一次的具有最短距離的回路。
     這個(gè)問題可分為對(duì)稱旅行商問題(dij=dji,,任意i,j=1,2,3，…,n)和非對(duì)稱旅行商
問題(dij≠dji,,任意i,j=1,2,3，…,n)。
     若對(duì)于城市v={v1,v2,v3,…,vn}的一個(gè)訪問順序?yàn)閠=(t1,t2,t3,…,ti,…,tn),其中
ti∈v(i=1,2,3,…,n)，且記tn+1= t1，則旅行商問題的數(shù)學(xué)模型為：
      min      l=σd(t(i),t(i+1))   （i=1,…,n）
     旅行商問題是一個(gè)典型的組合優(yōu)化問題，并且是一個(gè)np難問題，其可能的路徑數(shù)目
與城市數(shù)目n是成指數(shù)型增長的，所以一般很難精確地求出其最優(yōu)解，本文采用遺傳算法
求其近似解。
     遺傳算法：
初始化過程：用v1,v2,v3,…,vn代表所選n個(gè)城市。定義整數(shù)pop-size作為染色體的個(gè)數(shù)
，并且隨機(jī)產(chǎn)生pop-size個(gè)初始染色體，每個(gè)染色體為1到18的整數(shù)組成的隨機(jī)序列。
適應(yīng)度f的計(jì)算：對(duì)種群中的每個(gè)染色體vi，計(jì)算其適應(yīng)度，f=σd(t(i),t(i+1)).
評(píng)價(jià)函數(shù)eval(vi)：用來對(duì)種群中的每個(gè)染色體vi設(shè)定一個(gè)概率，以使該染色體被選中
的可能性與其種群中其它染色體的適應(yīng)性成比例，既通過輪盤賭，適應(yīng)性強(qiáng)的染色體被
選擇產(chǎn)生后臺(tái)的機(jī)會(huì)要大，設(shè)alpha∈(0,1)，本文定義基于序的評(píng)價(jià)函數(shù)為eval(vi)=al
pha*(1-alpha).^(i-1) 。[隨機(jī)規(guī)劃與模糊規(guī)劃]
選擇過程：選擇過程是以旋轉(zhuǎn)賭輪pop-size次為基礎(chǔ)，每次旋轉(zhuǎn)都為新的種群選擇一個(gè)
染色體。賭輪是按每個(gè)染色體的適應(yīng)度進(jìn)行選擇染色體的。
    step1 、對(duì)每個(gè)染色體vi,計(jì)算累計(jì)概率qi，q0=0;qi=σeval(vj)    j=1,…,i;i=1,
…pop-size.
    step2、從區(qū)間(0,pop-size)中產(chǎn)生一個(gè)隨機(jī)數(shù)r；
    step3、若qi-1    step4、重復(fù)step2和step3共pop-size次，這樣可以得到pop-size個(gè)復(fù)制的染色體。
grefenstette編碼：由于常規(guī)的交叉運(yùn)算和變異運(yùn)算會(huì)使種群中產(chǎn)生一些無實(shí)際意義的
染色體，本文采用grefenstette編碼《遺傳算法原理及應(yīng)用》可以避免這種情況的出現(xiàn)
。所謂的grefenstette編碼就是用所選隊(duì)員在未選（不含淘汰）隊(duì)員中的位置，如：
           8 15 2 16 10 7 4 3 11 14 6 12 9 5 18 13 17 1
           對(duì)應(yīng)：
           8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1。
交叉過程：本文采用常規(guī)單點(diǎn)交叉。為確定交叉操作的父代，從 到pop-size重復(fù)以下過
程：從[0，1]中產(chǎn)生一個(gè)隨機(jī)數(shù)r，如果r             將所選的父代兩兩組隊(duì)，隨機(jī)產(chǎn)生一個(gè)位置進(jìn)行交叉，如：
           8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1
           6 12 3 5 6 8 5 6 3 1 8 5 6 3 3 2 1 1
交叉后為：
          8 14 2 13 8 6 3 2 5 1 8 5 6 3 3 2 1 1
          6 12 3 5 6 8 5 6 3 7 3 4 3 2 4 2 2 1
變異過程：本文采用均勻多點(diǎn)變異。類似交叉操作中選擇父代的過程，在r 選擇多個(gè)染色體vi作為父代。對(duì)每一個(gè)選擇的父代，隨機(jī)選擇多個(gè)位置，使其在每位置
按均勻變異（該變異點(diǎn)xk的取值范圍為[ukmin,ukmax],產(chǎn)生一個(gè)[0，1]中隨機(jī)數(shù)r，該點(diǎn)
變異為x'k=ukmin+r(ukmax-ukmin)）操作。如：
          8 14 2 13 8 6 3 2 5 7 3 4 3 2 4 2 2 1
       變異后：
         8 14 2 13 10 6 3 2 2 7 3 4 5 2 4 1 2 1
反grefenstette編碼：交叉和變異都是在grefenstette編碼之后進(jìn)行的，為了循環(huán)操作
和返回最終結(jié)果，必須逆grefenstette編碼過程，將編碼恢復(fù)到自然編碼。
循環(huán)操作：判斷是否滿足設(shè)定的帶數(shù)xzome，否，則跳入適應(yīng)度f的計(jì)算；是，結(jié)束遺傳
操作，跳出。 



matlab 代碼



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;