dhb_tsp.m

段海滨老师的TSP类型

% function bestroute = fun_basic_TSP(aN,initao,alpha,beta,Q,rou,ncmax)
% 


cities=reshape(citymat,3,cN);  % get the city coordinates in TSP
x=cities(2,:);
y=cities(3,:);

for i=1:cN
    for j=1:cN
        dist(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
    end
end

for i=1:cN
    for j=1:cN
        if j~=i
            tao(i,j)=initao;
            yita(i,j)=1/dist(i,j);
        end
    end
end

initao=initao.*ones(cN,cN);
detatao=zeros(cN,cN);
bestsolu=1e13;

%% This is the phase in which ants build their tours

for nc=1:ncmax
    tabu=zeros(aN,cN);
    for k=1:aN
        tabu(k,g(nc,k))=1;
        maxp(1)=0;
        for j=1:cN
            if tabu(k,j)==0
                psum_medium0(1,j)=(tao(g(nc,k),j)^alpha).*(yita(g(nc,k),j)^beta);
            else
                psum_medium0(1,j)=0;
            end
        end
        psum_medium=psum_medium0.';
        psum(k,1)=sum(psum_medium(:,1));
        for j=1:cN
            if tabu(k,j)==0
                p(1,j)=(tao(g(nc,k),j)^alpha).*(yita(g(nc,k),j)^beta)/psum(k,1);
            else
                p(1,j)=0;       % maybe it should be p(1,j)=0
            end
            if maxp(1)<p(1,j)
                maxp(1)=p(1,j);
            end
        end
        lindex=find(maxp(1)==p(1,:));
        size1=[1,n];
        [Rind,Cind]=ind2sub(size1,lindex(1));
        solu_medium(k,1)=dist(g(nc,k),Cind(1));
        route(k,1)=Cind(1);
        tabu(k,Cind(1))=1;
        tao(g(nc,k),Cind(1))=(1-rou)*tao(g(nc,k),Cind(1))+rou*initao(g(nc,k),Cind(1));  % local updating for the first iteration
        tao(Cind(1),g(nc,k))=(1-rou)*tao(Cind(1),g(nc,k))+rou*initao(Cind(1),g(nc,k));
        solu(k)=solu_medium(k,1);
        
        for s=2:(cN-1)
            Cind(s)=0;
            Rind(s)=0;
            lindex(s)=0;
            maxp(s)=0;
            for j=1:cN
                if tabu(k,j)==0
                    psum_medium0(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta);
                else
                    psum_medium0(s,j)=0;
                end
            end
            psum_medium=psum_medium0.';
            psum(k,s)=sum(psum_medium(:,s));
            for j=1:cN
                if tabu(k,j)==0
                    p(s,j)=(tao(route(k,s-1),j)^alpha).*(yita(route(k,s-1),j)^beta)/psum(k,s);
                else
                    p(s:(cN-1),j)=0;      
                end
                if maxp(s)<p(s,j)
                    maxp(s)=p(s,j);
                end
            end
            lindex=find(maxp(s)==p);
            size2=[cN-1,cN];
            [Rind(s),Cind(s)]=ind2sub(size2,lindex(1));
            solu_medium(k,s)=dist(Cind(s-1),Cind(s));
            route(k,s)=Cind(s);
            tabu(k,Cind(s))=1;
            tao(Cind(s-1),Cind(s))=(1-rou)*tao(Cind(s-1),Cind(s))+rou*initao(Cind(s-1),Cind(s));  % local updating for the first iteration
            tao(Cind(s),Cind(s-1))=(1-rou)*tao(Cind(s),Cind(s-1))+rou*initao(Cind(s),Cind(s-1));
            solu(k)=solu(k)+solu_medium(k,s);
        end
        tao(Cind(cN-1),g(nc,k))=(1-rou)*tao(Cind(cN-1),g(nc,k))+rou.*initao(Cind(cN-1),g(nc,k));
        tao(g(nc,k),Cind(cN-1))=(1-rou)*tao(g(nc,k),Cind(cN-1))+rou.*initao(g(nc,k),Cind(cN-1));  % local updating for the first iteration
        solu_medium(k,cN)=dist(Cind(cN-1),g(nc,k));
        solu(k)=solu(k)+solu_medium(k,cN);
        soluNC(nc,k)=solu(k);
        bestsoluNC(nc)=soluNC(nc,1);
    end
    for k=1:m
        if bestsoluNC(nc)>soluNC(nc,k)
            bestsoluNC(nc)=soluNC(nc,k);
        end
    end
    
%% In this phase global updating occurs

    lindex1=find(bestsoluNC(nc)==soluNC(nc,:));
    size3=[nc,aN];
    [Rind1(nc),Cind1(nc)]=ind2sub(size3,lindex1(1));
    bestrouteNC(nc,:)=route(Cind1(nc),:);
    [m,n]=size(lindex1);
    for i=1:n
        [Rind1t(nc,i),Cind1t(nc,i)]=ind2sub(size3,lindex1(i));
    
        detatao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),1))=Q/solu(Cind1t(nc,i));
        detatao(route(Cind1t(nc,i),1),g(nc,Cind1t(nc,i)))=Q/solu(Cind1t(nc,i));
        tao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),1))=...
            (1-rou).*tao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),1))+...
            rou.*detatao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),1));
        tao(route(Cind1t(nc,i),1),g(nc,Cind1t(nc,i)))=...
            (1-rou).*tao(route(Cind1t(nc,i),1),g(nc,Cind1t(nc,i)))+...
            rou.*detatao(route(Cind1t(nc,i),1),g(nc,Cind1t(nc,i)));
    
        detatao(route(Cind1t(nc,i),n-1),g(nc,Cind1t(nc,i)))=Q/solu(Cind1t(nc,i));
        detatao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),n-1))=Q/solu(Cind1t(nc,i));
        tao(route(Cind1t(nc,i),n-1),g(nc,Cind1t(nc,i)))=...
            (1-rou).*tao(route(Cind1t(nc,i),n-1),g(nc,Cind1t(nc,i)))+...
            rou.*detatao(route(Cind1t(nc,i),n-1),g(nc,Cind1t(nc,i)));
        tao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),n-1))=...
            (1-rou).*tao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),n-1))+...
            rou.*detatao(g(nc,Cind1t(nc,i)),route(Cind1t(nc,i),n-1));
    
        for s=2:n-1
            detatao(route(Cind1t(nc,i),s),route(Cind1t(nc,i),s-1))=Q/solu(Cind1t(nc,i));
            detatao(route(Cind1t(nc,i),s-1),route(Cind1t(nc,i),s))=Q/solu(Cind1t(nc,i));
            tao(route(Cind1t(nc,i),s),route(Cind1t(nc,i),s-1))=...
                (1-rou).*tao(route(Cind1t(nc,i),s),route(Cind1t(nc,i),s-1))+...
                rou.*detatao(route(Cind1t(nc,i),s),route(Cind1t(nc,i),s-1));
            tao(route(Cind1t(nc,i),s-1),route(Cind1t(nc,i),n-1))=...
                (1-rou).*tao(route(Cind1t(nc,i),s-1),route(Cind1t(nc,i),s))+...
                rou.*detatao(route(Cind1t(nc,i),s-1),route(Cind1t(nc,i),s));
        end
    end
    if bestsolu>bestsoluNC(nc)
        bestsolu=bestsoluNC(nc);
    end
    
    lindex2=find(bestsolu==bestsoluNC);
    size4=[1,nc];
    [Rind2,Cind2]=ind2sub(size4,lindex2(1));
    bestroute(1,:)=bestrouteNC(Cind2,:);
    initao=tao;
end
    
%% Output the best result for TSP

for nc=1:ncmax
    bestrouteNC(nc,cN)=g(nc,Cind1(nc));
    t(nc)=nc;
end
bestroute(1,cN)=bestrouteNC(Cind2,cN);
plot(x(bestroute(cN)),y(bestroute(cN)),'*');
hold on;
for i=1:cN
    plot(x(i),y(i),'o');
    hold on;
end
hold on;
line([x(bestroute(cN)) x(bestroute(1))],[y(bestroute(cN)) y(bestroute(1))]);
hold on;
for j=1:cN-1
    line([x(bestroute(j)) x(bestroute(j+1))],[y(bestroute(j)) y(bestroute(j+1))]);
    hold on;
end
hold off;
xlabel('X coordinate');ylabel('Y coordinate');
title('Best tour in NCmax iterations');