fortyeight_new_pso_aco.m

PSO-ACO-TSP算法

%File: pso_aco.m
clear all
x=[6734 2233 5530 401 3082 7608 7573 7265 6898 1112 5468 5989 4706 4612 6347 6107 7611 7462 7732 5900 4483 6101 5199 1633 ...
   4307 675 7555 7541 3177 7352 7545 3245 6426 4608 23 7248 7762 7392 3484 6271 4985 1916 7280 7509 10 6807 5185 3023];
y=[1453 10 1424 841 1644 4458 3716 1268 1885 2049 2606 2873 2674 2035 2683 669 5184 3590 4723 3561 3369 1110 2182 2809 ...
   2322 1006 4819 3981 756 4506 2801 3305 3173 1198 2216 3779 4595 2244 2829 2135 140 1569 4899 3239 2676 2993 3258 1942];
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