fortyeight_pso_aco.m

PSO-ACO-TSP算法

%File: pso_aco.m
clear all
x=[41 37 54 25 7 2 68 71 54 83 64 18 22 83 91 20 24 31 70 1 50 15 40 90 ...
   25 24 58 71 74 87 18 13 82 62 58 45 41 44 4 38 42 69 71 78 76 40 7 32];
y=[94 84 67 62 64 99 58 44 62 69 60 54 60 46 38 90 19 8 50 70 80 11 24 20 ...
   38 42 69 71 78 76 40 7 32 35 21 26 35 66 50 25 24 58 71 74 87 18 13 82];
n=48;%48 cities
c=100;
q=1000000;
NC=100;
r=0.9;
a=1.5;
b=2;
m=48;%the number of ants
for i=1:n
    for j=1:n
         dij(i,j)=sqrt((x(i)-x(j))^2+(y(i)-y(j))^2);
    end
end
for i=1:n
    dij(i,i)=0.01;
end
minl1=1000000;
t=ones(n)*c;
for i=1:100  %i is what??
    road(i,:)=randperm(n);
    ltspc(i)=ca_tsp(n,road(i,:),dij);
end
ltspsort=sort(ltspc);
l48=ltspsort(m);%选择m个初始解，m只蚂蚁
i1=0;
for i=1:100
     if ltspc(i)<=l48
           for k=1:n-1
                t(road(i,k),road(i,k+1))=t(road(i,k),road(i,k+1))+10;
                t(road(i,k+1),road(i,k))=t(road(i,k),road(i,k+1));
           end
           t(road(i,1),road(i,n))=t(road(i,1),road(i,n))+10;
           t(road(i,n),road(i,1))=t(road(i,1),road(i,n));
           i1=i1+1;
           pcbest(i1,:)=road(i,:);    %初始化pcbest,路径
           plbest(i1)=ltspc(i);       %初始化pcbest,路径长度
     end
end
[glbest,j]=min(plbest);         %初始化glbest,路径长度
gcbest=pcbest(j,:);             %初始化gcbest,路径
for nc=1:NC;
     tabu=ones(m,n);
     tabu(:,1)=0;
     path=ones(m,n);
     for k=1:m
          for step=1:n-1
               ta=t.^a;
               tb=dij.^(-b);
               td=ta.*tb;
               pd=tabu(k,:).*td(path(k,step),:);
               pk=pd/sum(pd);
               rk=rand;
               spk=0;
               j=1;
               while j<=n
                 if rk<spk+pk(j)
                      break;
                 else
                    spk=spk+pk(j);
                    j=j+1;
                 end
               end
               if j==n+1  %added at 2006-10-3-------unless if will be "Index exceeds matrix dimensions." 
                   j=n;
               end        %added at 2006-10-3
               tabu(k,j)=0;
               path(k,step+1)=j;
           end
      end
      for i=1:m
           ltsp(i)=ca_tsp(n,path(i,:),dij);
           %交叉操作
           %四种交叉程序分别为cross_tsp_1,cross_tsp_2,cross_tsp_3,cross_tsp_4
      %     path1(i,:)=cross_tsp_2(path(i,:),gcbest,n);
      %     path1(i,:)=cross_tsp_2(path(i,:),pcbest(i,:),n);                                    %%%%10-31 what?

           %变异操作
           %四种变异程序分别为muta_tsp_1,muta_tsp_2,muta_tsp_3,muta_tsp_4
      %     path1(i,:)=mutation_1(path1(i,:),n);
 
           %计算新路径长度之和
      %     ltsp1(i)=ca_tsp(n,path1(i,:),dij);

           %求个体极值
     %      if ltsp1(i)<ltsp(i)
     %            ltsp(i)=ltsp(i);
     %            path(i,:)=path1(i,:);
     %      end
    %       if ltsp(i)<plbest(i);
    %             plbest(i)=ltsp(i);
    %             pcbest(i,:)=path(i,:);
    %       end
     end
     %找全局极值
     [glbest,j]=min(plbest);
     gcbest=pcbest(j,:);
     dt=zeros(n);
     for i=1:m
          for k=1:n-1
               dt(path(i,k),path(i,k+1))=dt(path(i,k),path(i,k+1))+q/ltsp(i);
               dt(path(i,k+1),path(i,k))=dt(path(i,k),path(i,k+1));
          end
          dt(path(i,n),path(i,1))=dt(path(i,n),path(i,1))+q/ltsp(i);
          dt(path(i,1),path(i,n))=dt(path(i,n),path(i,1));
     end
     t=r*t+dt;
end
ta=t.^a;
tb=dij.^(-b);
tb=ta.*tb;
k=3;
ts(1)=1;
td(:,1)=0;
[my,i]=max(td(1,:));
ts(2)=i;
td(:,i)=0;
while k<=n
    [my,i]=max(td(i,:));
    ts(k)=i;
    td(:,i)=0;
    k=k+1;
end
ltsp0=ca_tsp(n,ts,dij);
if glbest<ltsp0
     ts=gcbest;
     ltsp0=glbest;
end
k=1;
while k<=n
    xl(k)=x(ts(k));
    yl(k)=y(ts(k));
    k=k+1;
end
xl(n+1)=xl(1);
yl(n+1)=yl(1);
line(xl,yl)
hold on
plot(x,y,'o');
ltsp0