function TSPga(tspdist,termops,popsize,pm,alpha)%全國計算機仿真大獎賽第三題程序(旅行商TSP問題)%tspdist:距離矩陣%termops:種群代數 ，推薦40-300%popsize:每代染色體的個數， 推薦60-450 %pm:變異概率， 推薦0.1-0.3 pm越大，參與均勻變異的染色體越多, 比如pm=0.2，則參與變異的約為20%。%alpha:評價函數eval(vi)=alpha*(1-alpha).^(i-1)， 推薦0.04-0.2， alpha越大，則越擇優復制進入下一代，alpha=0.1。%bestpop:返回的最優種群%recordcost:進化軌跡%使用例子:%          tspdist=[ 0    17    13    24    10;%                   10     0    20     9     6;%                   17    29     0    21    28;%                   12    10    22     0    19;%                   12    18    31    20     0 ];%          TSPga(tspdist,20,15,0.25,0.18);%返回:%   Program begins to solve your TSP.................%   OK. The best touring path is:5  1  3  4  2  5%   And the best total traveling expense amounts to:62%   It takes me  0.25  second.%   You can get more information in the evolution graph window% ===============程序運行環境==========================% Windows XP , Pentium4 CPU 1.6GHz , 256MB , Matlab 6.5% =====================================================% Written by LiWeichao, OuyangBincitynum=size(tspdist,2);n=nargin;if n<2   % 設置缺省參數    termops=60;    popsize=100;    pm=0.2;    alpha=0.1;elseif n<3    popsize=100;    pm=0.2;    alpha=0.1;elseif n<4    pm=0.2;    alpha=0.1;elseif n<5    alpha=0.1;enddisp('  Program begins to solve your TSP.................');gt=initializega(popsize,citynum);% 初始化染色體庫（也叫原始種群）tic;for i=1:termops     % i         % 顯示進化代數,如果想運行時觀察當前在進化哪一代，則請把i前的“%”去掉。    rawcost=fitness(tspdist,gt);  % (rawcost是行向量) 計算目標函數      [x,y]=find(rawcost==min(rawcost)); % 找最優     recordcost(i)=rawcost(y(1));    bestpop=gt(y(1),:); % 當代最優染色體    gt=select(gt,rawcost,alpha); % 按alpha選擇和復制染色體    for j=1:2:(popsize-1)        [gt(j,:),gt(j+1,:)]=liucrossover(tspdist,gt(j,:),gt(j+1,:)); % 劉海交叉法    end    gt=grefenstette(gt); % grefenstette編碼    gt=mutation(gt,pm); % 按變異概率pm均勻變異    gt=congrefenstette(gt);% grefenstette反編碼endruntime=toc;rawcost=fitness(tspdist,gt);  % (行向量) 計算目標函數  [x,y]=find(rawcost==min(rawcost)); % 找最優 recordcost(termops+1)=rawcost(y(1));bestpop=gt(y(1),:); % 當代最優染色體bestpop(citynum+1)=bestpop(1);  % 構成閉合路徑plot([0:termops],recordcost,'.-');   % 畫種群進化圖set(gca,'xlim',[0 termops]);set(gcf,'NumberTitle','off','Name','Evolution Graph');grid on;xlabel('染色體種群進化代數');ylabel('每一代種群的最優值');legend(strcat('種群代數:',num2str(termops)),strcat('染色體的個數:',num2str(popsize)),...    strcat('變異概率:',num2str(pm)),strcat('評價函數alpha:',num2str(alpha)),strcat('城市個數:',num2str(citynum)));title('染色體種群進化過程');disp(['  OK. The best touring path is:',num2str(bestpop)]);% 顯示結果disp(['  And the best total traveling expense amounts to:',num2str(recordcost(i))]);disp(['  It takes me ',num2str(runtime),' second.']);disp('  You can get more information in the evolution graph window');% =========================================================% 以下是子函數% =========================================================function [t]=initializega(popsize,citynum)%%%%%%初始化基因庫for i=1:popsize    t(i,:)=randperm(citynum);end% -----------------------------------------------------------function rawcost=fitness(tspdist,gt)%%%%%%%計算目標適應度函數[m,n]=size(gt);for k=1:m    for i=1:n-1        distan(k,i)=tspdist(gt(k,i),gt(k,i+1));    end    distan(k,n)=tspdist(gt(k,n),gt(k,1));    rawcost(k)=sum(distan(k,:));end% -----------------------------------------------------------function [gt]=select(gt,rawcost,alpha)%%%%%%%選擇和復制染色體[m,n]=size(gt);gtb=gt;[aftercost,aftersort]=sort(rawcost,2);for k=1:m,    gt(k,:)=gtb(aftersort(k),:);endgtb=gt;for i=1:m    evalv(i)=alpha*(1-alpha).^(i-1);endq=[-1 cumsum(evalv)];qmax=max(q);for k=1:m    r=qmax*rand(1);    for j=1:m        if r>q(j)&r<=q(j+1)            gt(k,:)=gtb(j,:);            break;        end    endend% -------------------------------------------function [chromc,chromd]=liucrossover(tspdist,chroma,chromb)%%%染色體交叉，利用劉海的方法。nn=length(chroma);chroma=circshift(chroma',ceil(rand(1)*4))';chromc=zeros(1,nn);chromc(1)=chroma(1);for i=2:nn,    k=1;    while k>0,        am=(find(chroma==chromc(i-1))+k);        bm=(find(chromb==chromc(i-1))+k);        if am>nn,            if mod(am,nn)~=0,                am=mod(am,nn);            else ,                am=nn;            end        end        if bm>nn,            if mod(bm,nn)~=0,                bm=mod(bm,nn);            else ,                bm=nn;            end        end        if tspdist(chromc(i-1),chroma(am))<=tspdist(chromc(i-1),chromb(bm)),            if isempty(find(chromc([1:i-1])==chroma(am))),                chromc(i)=chroma(am); k=0;            elseif isempty(find(chromc([1:i-1])==chromb(bm))),                chromc(i)=chromb(bm); k=0;            else ,                k=k+1;            end        else ,            if isempty(find(chromc([1:i-1])==chromb(bm))),                chromc(i)=chromb(bm); k=0;            elseif isempty(find(chromc([1:i-1])==chroma(am))),                chromc(i)=chroma(am); k=0;            else ,                k=k+1;            end        end    endendchromd=chroma;chroma=chromb;chromb=chromd;chromd(1)=chroma(1);for i=2:nn,    k=1;    while k>0,        am=(find(chroma==chromd(i-1))+k);        bm=(find(chromb==chromd(i-1))+k);        if am>nn,            if mod(am,nn)~=0,                am=mod(am,nn);            else ,                am=nn;            end        end        if bm>nn,            if mod(bm,nn)~=0,                bm=mod(bm,nn);            else ,                bm=nn;            end        end        if tspdist(chromd(i-1),chroma(am))<=tspdist(chromd(i-1),chromb(bm)),            if isempty(find(chromd([1:i-1])==chroma(am))),                chromd(i)=chroma(am); k=0;            elseif isempty(find(chromd([1:i-1])==chromb(bm))),                chromd(i)=chromb(bm); k=0;            else ,                k=k+1;            end        else ,            if isempty(find(chromd([1:i-1])==chromb(bm))),                chromd(i)=chromb(bm); k=0;            elseif isempty(find(chromd([1:i-1])==chroma(am))),                chromd(i)=chroma(am); k=0;            else ,                k=k+1;            end        end    endend% ---------------------------------------------------function [gt]=mutation(gt,pm) %均勻變異%%%%%%%每條染色體獲得的變異概率是均等的[m,n]=size(gt);mutpot=ceil(rand*n/5)+1;  % 隨機確定染色體會變異的基因個數ran=rand(1,m);          % 給每條染色體隨機分配一個0-1均勻分布的數r=rand(1,mutpot);rr=floor(n*rand(1,mutpot)+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;        gt(mu(k),rr(i))=umin(i)+floor((umax(i)-umin(i))*r(i));    endend% ---------------------------------------------------function gt=gexchange(gt)nn=ceil(rand(1,round(length(gt)/5))*length(gt));% ---------------------------------------------------function [g]=grefenstette(t)%%%%%%%grefenstette編碼[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    endend% ---------------------------------------------------function [t]=congrefenstette(g)%%%%%%%grefenstette反編碼[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))=[];    endend % ---------------------------------------------------function tspsol=solvetsp(cm)% SOLVETSP  Solve the traveling salesman problem.%     TSPSOL=SOLVETSP(CM) where CM is a square matrix listing the travel%     cost between one site and another site, obtains a structure array%     which contains the minimal-cost way and the minimal cost,using the%     global search method.However,only one solution is given eventually %     even if there would be more than one.%     Direction of the way? for example, way=[1 5 4 3 2 1] means a touring%     cycle from the first site to the 5th site,(1-->5-->4-->3-->2-->1).By %     the way,the site marking at the Nth column or the Nth row of the cost %     matrix is mapped as N,an integer such as 5 or others.%     Example:%         cm=[ 0    17    13    24    10;%             10     0    20     9     6;%             17    29     0    21    28;%             12    10    22     0    19;%             12    18    31    20     0 ];%         mysol=solvetsp(cm);%% Written by OuyangBin, 2004-12-27% Variable Statement:%    cm is the cost matrix or distance matrix.%    pm is the permutation matrix;%    lpm is the length of pm array;%    tc is the total cost;%    tct is the temporary total cost;%    tk is the tag of minimal-cost way in the permutation array;%    tspsol is the solution structre array obtained eventually.% ======================Running Envirnment==========================% Windows XP , Pentium4 1.6GHz , 256MB , Matlab 6.5% ==================================================================% pre-processing[row col]=size(cm);if row~=col,    disp('   Bad TSP Matrix. Actually, The TSP matrix must be a square matrix.');    tspsol='Null';    return;endif sum(abs(diag(cm,0)))~=0,    disp('   Bad TSP Matrix.Actually, The main diagonal elements of TSP matrix must be zeros.');    tspsol='Null';    return;endif row>10,    disp('   Sorry,The scale of cost matrix is too great such that computer will runs out of memory(3GB).');    tspsol='Null';    return;end% main body of the program.% generate the cycle matrix.pm=allperm([2:row]);lpm=length(pm);pm=[ones(lpm,1) pm ones(lpm,1)];% refine the cycletic;tc=inf;for i=1:lpm,    tct=0;    for j=1:row,        tct=tct+cm(pm(i,j),pm(i,j+1));    end    if tct<tc,        tc=tct;        tk=i;    endendtspsol.way=pm(tk,:);tspsol.cost=tc;disp(strcat('  best touring path:',num2str(tspsol.way)));disp(strcat('  minimal cost:',num2str(tspsol.cost)));disp(strcat('  running time(second):',num2str(toc)));% end of the main program% the subroutine function P = allperm(V)% all possible permutation of input vector.V = V(:).'; % Make sure V is a row vectorn = length(V);if n <= 1, P = V;     return; endq = allperm(1:n-1);  % recursive callsm = size(q,1);P = zeros(n*m,n);P(1:m,:) = [n * ones(m,1) q];for i = n-1:-1:1,    t = q;    t(t == i) = n;    P((n-i)*m+1:(n-i+1)*m,:) = [i*ones(m,1) t]; % assign the next m    % rows in P.endP = V(P);v