TSP問題的遺傳算法+轉載

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


function [bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%
%————————————————————————
%[bestpop,trace]=ga(d,termops,num,pc,cxops,pm,alpha)
%d:距離矩陣
%termops:種群帶數
%num:每帶染色體的個數
%pc:交叉概率
%cxops:由于本程序采用單點交叉，交叉點的設置在本程序中沒有很好的解決，所以本文了采用定點，即第cxops，可以隨機產生。
%pm:變異概率
%alpha:評價函數eval(vi)=alpha*(1-alpha).^(i-1).
%bestpop:返回的最優種群
%trace:進化軌跡
%------------------------------------------------
%####@@@##版權所有！歡迎廣大網友改正，改進！##@@@####
%e-mail:tobysidney33@sohu.com
%####################################################
%
citynum=size(d,2);
n=nargin;
if n<2
    disp('缺少變量！！')
    disp('^_^開個玩笑^_^')
end
if n<2
    termops=500;
    num=50;
    pc=0.25;
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<3
    num=50;
    pc=0.25;
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<4
    pc=0.25;
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<5
    cxops=3;
    pm=0.30;
    alpha=0.10;
end
if n<6
    pm=0.30;
    alpha=0.10;
end
if n<7
    alpha=0.10;
end
if isempty(cxops)
    cxops=3;
end

[t]=initializega(num,citynum);
for i=1:termops
[l]=f(d,t);
[x,y]=find(l==max(l));
trace(i)=-l(y(1));
bestpop=t(y(1),:);
[t]=select(t,l,alpha);
[g]=grefenstette(t);
[g1]=crossover(g,pc,cxops);
[g]=mutation(g1,pm);  %均勻變異
[t]=congrefenstette(g);
end

---------------------------------------------------------
function [t]=initializega(num,citynum)
for i=1:num
    t(i,:)=randperm(citynum);
end
-----------------------------------------------------------
function [l]=f(d,t)
[m,n]=size(t);
for k=1:m
    for i=1:n-1
      l(k,i)=d(t(k,i),t(k,i+1));
    end
      l(k,n)=d(t(k,n),t(k,1));
      l(k)=-sum(l(k,:));
end
-----------------------------------------------------------
function [t]=select(t,l,alpha)
[m,n]=size(l);
t1=t;
[beforesort,aftersort1]=sort(l,2);%fsort from l to u
for i=1:n
    aftersort(i)=aftersort1(n+1-i);      %change 
end
for k=1:n;
    t(k,:)=t1(aftersort(k),:);
    l1(k)=l(aftersort(k));
end
t1=t;
l=l1;
for i=1:size(aftersort,2)
    evalv(i)=alpha*(1-alpha).^(i-1);
end
m=size(t,1);
q=cumsum(evalv);
qmax=max(q);
for k=1:m
    r=qmax*rand(1);
    for j=1:m
        if j==1&r<=q(1)
            t(k,:)=t1(1,:);
        elseif j~=1&r>q(j-1)&r<=q(j)
            t(k,:)=t1(j,:);
        end
    end
end
--------------------------------------------------
function [g]=grefenstette(t)
[m,n]=size(t);
for k=1:m
    t0=1:n;
   for i=1:n
       for j=1:length(t0)
           if t(k,i)==t0(j)
              g(k,i)=j;
              t0(j)=[];
              break
           end
       end
    end
end
-------------------------------------------
function [g]=crossover(g,pc,cxops)
[m,n]=size(g);
ran=rand(1,m);
r=cxops;
[x,ru]=find(ran<pc);
if ru>=2
    for k=1:2:length(ru)-1
       g1(ru(k),:)=[g(ru(k),[1:r]),g(ru(k+1),[(r+1):n])];
       g(ru(k+1),:)=[g(ru(k+1),[1:r]),g(ru(k),[(r+1):n])];
       g(ru(k),:)=g1(ru(k),:);
    end
end
--------------------------------------------
function [g]=mutation(g,pm)    %均勻變異
[m,n]=size(g);
ran=rand(1,m);
r=rand(1,3);      %dai gai jin
rr=floor(n*rand(1,3)+1);
[x,mu]=find(ran<pm);
for k=1:length(mu)
    for i=1:length(r)
        umax(i)=n+1-rr(i);
        umin(i)=1;
        g(mu(k),rr(i))=umin(i)+floor((umax(i)-umin(i))*r(i));
    end
end
---------------------------------------------------
function [t]=congrefenstette(g)
[m,n]=size(g);
for k=1:m
    t0=1:n;
   for i=1:n
      t(k,i)=t0(g(k,i));
      t0(g(k,i))=[];
   end
end