tsp.txt

matlab 源程序 说明见程序本身

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, OuyangBin

citynum=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;
end
disp('  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反编码
end
runtime=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),:);
end
gtb=gt;
for i=1:m
    evalv(i)=alpha*(1-alpha).^(i-1);
end
q=[-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
    end
end
% -------------------------------------------
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
    end
end
chromd=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
    end
end
% ---------------------------------------------------

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));
    end
end
% ---------------------------------------------------
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
    end
end

% ---------------------------------------------------
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))=[];
    end
end 
% ---------------------------------------------------






















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;
end
if sum(abs(diag(cm,0)))~=0,
    disp('   Bad TSP Matrix.Actually, The main diagonal elements of TSP matrix must be zeros.');
    tspsol='Null';
    return;
end
if 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 cycle
tic;
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;
    end
end
tspsol.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 vector
n = length(V);
if n <= 1, P = V; 
    return; 
end
q = allperm(1:n-1);  % recursive calls
m = 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.
end
P = V(P);




v