aca_ant_system.m

蚂蚁系统三种模型：蚁量模型、蚁密模型、蚁圈模型的源代码。

function [Alocation,Newbest,traceInfo]=aca_ant_system()
% 参考文献：Distributed Optimization by Ant Colonies三种模型
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%               蚁群算法初始化程序开始                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%% 问题定义：获取城市位置坐标、计算距离矩阵  %%%%%%%%%
InitOps=[];
[Location,DistMatrix,Ncities,Bestx,Lengx] = eil51init(InitOps); %China31cityinit Kroa100init
close all;
figure (5);
hold on;
minx=min(DistMatrix(:,1));
maxx=max(DistMatrix(:,1));
miny=min(DistMatrix(:,2));
maxy=max(DistMatrix(:,2));
minm=min(minx,miny);
maxm=max(maxx,maxy);
l=(maxm-minm)/100;
for i=1:Ncities
    plot(Location(i,1),Location(i,2),'*b');
    text (Location(i,1)+l,Location(i,2)+l,num2str(i));
end
for i=1:Ncities-1
    line([Location(Bestx(i),1),Location(Bestx(i+1),1)] , [Location(Bestx(i),2),Location(Bestx(i+1),2)]) ;
end
line([Location(Bestx(1),1),Location(Bestx(Ncities),1)] , [Location(Bestx(1),2),Location(Bestx(Ncities),2)]) ;
grid on,title(['全局最优解路线图-',num2str(Lengx)]),xlabel('横坐标'),ylabel('纵坐标');
legend('城市位置');
hold off ;


Flag=3;
%%%%%%%%%%               设定系统有关参数           %%%%%%%%%%
% M=round(Ncities/2);% 蚂蚁数量
M=Ncities;% 蚂蚁数量
MaxNc=500;% 最大代数
A=1;% 信息素因子
B1=5;% 启发信息因子
B2=2;% 启发信息因子
B3=5;% 启发信息因子
Q1=100;% 蚁密模型信息素更新因子
Q2=100;% 蚁量模型信息素更新因子
Q3=100;% 蚁周模型信息素更新因子
P1=0.9;% 蚁密模型挥发系数初值
P2=0.9;% 蚁量模型挥发系数初值
P3=0.7;% 蚁周模型挥发系数初值
switch Flag
       case 1
              P=P1;
              Q=Q1;
              B=B1;
       case 2
              P=P2;
              Q=Q2;
              B=B2;
       case 3
              P=P3;
              Q=Q3;
              B=B3;
       otherwise
              fprintf(1,'Error and Exit',gen,bval);  
end

%%%%%%%%%%   初始化信息素、启发信息矩阵、确定蚂蚁最初位置及允许矩阵    %%%%%%%%%%
Q0=100;%信息素初始值
Pheromone=Q0*ones(Ncities,Ncities);% 信息素初始矩阵;
Temp=Q0*ones(1,Ncities);
Pheromone=Pheromone-diag(Temp);
Temp=ones(1,Ncities);
Heuristic=1./(DistMatrix+diag(Temp));% 启发信息初始矩阵
Heuristic=Heuristic-diag(Temp);
%%%%%%%%%%                运行参数初始化             %%%%%%%%%%
Nc=1;% 第一代
Lbestdis=inf;%局部最优解
Cbestdis=inf;%全局最优解
Fnewbest=0;%当前最优解
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%               蚁群算法初始化程序结束                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  蚁群算法主循环开始                      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

while(Nc<=MaxNc)
   
      % 运行参数初始化
      RandL=round(rand(M,1)*Ncities+0.5);%蚂蚁最初位置
      Alocation0=zeros(M,Ncities+1);% 存放M+1个蚂蚁游历的路径及长度矩阵初始化
      Alocation0(:,1)=RandL;
      Allow0=repmat(1:Ncities,M,1);% 允许访问的城市矩阵初始化
      for Ak=1:M
          Allow0(Ak,RandL(Ak))=0;
      end
      Alocation=Alocation0;% 存放M个蚂蚁游历的路径及长度矩阵初始化
      Allow=Allow0; % 允许矩阵赋初值
     
      % M个蚂蚁选择Ncities个城市
      for Cityi=2:Ncities
          for Ak=1:M
              i=Alocation(Ak,Cityi-1);% 当前城市
              j=Select_for_as(Ak,i,Allow,A,B,Pheromone,Heuristic);% 依据Pij选择下一个城市j
              Alocation(Ak,Cityi)=j;
              Allow(Ak,j)=0;% 更新允许矩阵
              if Flag==1
                 Pheromone(i,j)=P*Pheromone(i,j)+Q;
                 Pheromone(j,i)=Pheromone(i,j);
              elseif Flag==2
                 Pheromone(i,j)=P*Pheromone(i,j)+Q./DistMatrix(i,j);
                 Pheromone(j,i)=Pheromone(i,j);
              end
          end
      end
      
      % 计算每个蚂蚁找到的路径的长度
      for Ak=1:M
          Alocation(Ak,:)=AcatspEval(Alocation(Ak,:),DistMatrix,Ncities);% 计算适应值即每个蚂蚁找到的路径长度
      end
      
      % 保存最优解;
      [Cbestdis,Am]=min(Alocation(1:M,Ncities+1));
      Cbest=Alocation(Am,:);
      
      % 给输出变量赋值
      traceInfo(Nc,1)=Nc; 		          %current generation
      traceInfo(Nc,2)= Cbest(Ncities+1);       %Best fittness
      traceInfo(Nc,3)=mean(Alocation(:,Ncities+1));     %Avg fittness
      traceInfo(Nc,4)=std(Alocation(:,Ncities+1)); %计算标准方差
      if (Cbestdis<Lbestdis)||(Nc==1)
         Fnewbest=Fnewbest+1;
         Newbest(Fnewbest,1)=Nc;
         Newbest(Fnewbest,2:Ncities+2)=Cbest;
         Lbest=Cbest;
         Lbestdis=Cbest(Ncities+1);
      end
      % 更新信息素矩阵
      if Flag==3
         Pheromone=P*Pheromone;% 挥发
         for Ak=1:M
             Iindex=Alocation(Ak,1:Ncities);
             Jindex=[Alocation(Ak,2:Ncities) Alocation(Ak,1)];
             for k=1:Ncities
                 Pheromone(Iindex(k),Jindex(k))=Pheromone(Iindex(k),Jindex(k))+Q./Alocation(Ak,Ncities+1);
                 Pheromone(Jindex(k),Iindex(k))=Pheromone(Iindex(k),Jindex(k));
             end
         end
      end
      Nc=Nc+1;
  end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                  蚁群算法主循环结束                      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                 蚁群算法输出模块开始                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 以图形方式显示运行结果开始 %%
figure(1);
hold on;
plot ( traceInfo (: , 1) , traceInfo (: , 2) ) ;
grid on,title('离线性能'),xlabel('代数'),ylabel('最优个体适应值');
legend('最优个体适应值');
hold off;

figure(2);
hold on;
plot( traceInfo (: , 1) , traceInfo (: , 3) ) ;
grid on,title('在线性能'),xlabel('代数'),ylabel('平均适应值');
legend('平均适应值');
hold off;

figure(3);
hold on;
plot( Newbest (: , 1) , Newbest (: , Ncities+2) ) ;
grid on,title('全局最优解变化'),xlabel('代数'),ylabel('全局最优解适应度');
legend('全局最优解适应度');
hold off;

figure (4);
hold on;
minx=min(DistMatrix(:,1));
maxx=max(DistMatrix(:,1));
miny=min(DistMatrix(:,2));
maxy=max(DistMatrix(:,2));
minm=min(minx,miny);
maxm=max(maxx,maxy);
l=(maxm-minm)/100;

for i=1:Ncities
    plot(Location(i,1),Location(i,2),'*b');
    text (Location(i,1)+l,Location(i,2)+l,num2str(i));
end
for i=2:Ncities
    line([Location(Newbest(Fnewbest,i),1),Location(Newbest(Fnewbest,(i+1)),1)] , [Location(Newbest(Fnewbest,i),2),Location(Newbest(Fnewbest,(i+1)),2)]) ;
end
line([Location(Newbest(Fnewbest,2),1),Location(Newbest(Fnewbest,Ncities+1),1)] , [Location(Newbest(Fnewbest,2),2),Location(Newbest(Fnewbest,Ncities+1),2)]) ;
grid on,title(['全局最优解路线图-',num2str(Lbestdis)]),xlabel('横坐标'),ylabel('纵坐标');
legend('城市位置');
hold off ;


%% 以图形方式显示运行结果开始 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%