图说蚁群算法(ACO)附源码_蚁群算法图解-CSDN博客

本文链接：https://blog.csdn.net/Flag_ing/article/details/107010604

PS：再过几天就可以返校收拾东西了，想想还有点小激动呐hhh~回想疫情宅家的这半年，真是一段充满了焦虑、惊喜、忙碌、充实又时而无聊的时光。。。返校只能待三天又让人有点小遗憾呐。就想着趁还在家这几天，把毕设用到的一些算法整理一下吧！

废话不多说了，这次要讲的是蚁群算法(下称ACO)。

一、初识ACO：这是个啥玩意?

二、ACO的数学原理

三、MATLAB代码

一、初识ACO：这是个啥玩意?

ACO是一种用来规划最优路径的概率型算法，其原理参考了蚁群在寻找食物过程中发现最短路径的行为。

这里我们先不谈高大上的数学原理，我们先定性地理解下蚁群算法的基本原理。既然ACO参考了蚁群的觅食行为，那我们就不妨先了解下蚁群是如何在觅食过程中找到最优路径的吧，直接上图解释吧。

蚂蚁种群可以在行走过程中释放与感知一种叫做信息素的物质，并通过这种能力进行种群内的信息传递。每只蚂蚁都会在行走的路径上释放这种信息素，另外蚂蚁在选择出行路径时会一定程度上参考不同路径上其他蚂蚁留下的信息素浓度，并趋向于沿着信息素浓度较高的路径行走(如图1中的线路3所示)，越是短的路径上残留的信息素就会越多，久而久之，经过一段时间的摸索后，整个蚁群就会沿着信息素浓度最高的路径即最优路径行走了。

蚁群算法就是参照了上述蚁群觅食的过程而提出来的一种寻找最优路径的方法，最早用以解决旅行商(TSP)问题。蚁群算法的基本流程如下所示：

初始化蚂蚁的数量(足够多)，初始化每条路径上的信息素为一相等的常数；
所有蚂蚁同时开始并行搜索路径，并在行进的路径上释放信息素。蚂蚁选择下一条要走路径的概率受到路径上残留信息素的量以及路径长短的综合影响；
待同一迭代批次(同一世代)的蚂蚁都完成对路径的搜索之后再更新所有路径上的信息素，包括新加的信息素与挥发的信息素(自然界中真实的信息素存在随时间蒸发的现象)；
更新当次迭代中的最优路径，重新初始化所有的蚂蚁，准备下一轮迭代；
重复步骤2~4，直到所有蚂蚁都选择同一条完整的路径或者达到最大迭代次数，算法结束，以当前解作为问题的最优解。

在以上算法描述中，需要遵循两个原则：一个是在同一世代中的蚂蚁之间互不干扰，即同一世代的蚂蚁在本次行走中释放的信息素对该世代中其他的蚂蚁无影响。另外一个就是每只蚂蚁在每条路径上只能走一次，为此可以设置禁忌表(Tabu)来控制。

二、ACO的数学原理

在一中，我们了解了自然界蚁群的觅食过程，而将这一过程进行数学化描述就是ACO算法了。在上面对蚁群觅食的描述中我们可以发现几个关键词：蚁群( ${\color{Red} m}$ )，信息素( ${\color{Red} \tau }$ )，信息素变化量( ${\color{Red} \Delta \tau }$ )，信息素挥发(挥发因子 ${\color{Red} \rho }$ )，路径选择(转移概率 ${\color{Red} p}$ )，久而久之(迭代次数 ${\color{Red} K}$ )，以及路径长度( ${\color{Red} L}$ )。如何对这些关键词数学化呐？且听我慢慢道来。

假设我们有图2这么一张栅格地图，起点终点图中已标出，假设蚂蚁只能从一个栅格中点到下一个相邻(八邻域)栅格中点这样一步一步运动，黑色块为假想的障碍物区域，很明显，这是一次静态的路径规划，我们的目的就是规划一条从起点到终点的最短无碰路径。

看着上面的栅格图，我们可以想象这样的场景：一批(迭代批次)的蚂蚁从起点开始无重复地试探爬行，并在爬过的路径上留下信息素，直到爬到终点为止，之后这批蚂蚁下岗养老，新的一批蚂蚁上岗并开始执行同样的操作。后面的蚂蚁在前蚁留下的信息素的基础上并结合八邻域点距离终点的距离来综合考虑(转移概率)下一步应该走往哪个栅格点。对比记录每个批次中的最短路径。这样周而复始。。。结果你懂吧？——不懂？——那继续往下看吧

1、首先是蚁群数量 ${\color{Red} m}$ ，以及迭代次数 ${\color{Red} K}$

数量都不要太少，得有个一二十个往上吧，但也别太多，太多电脑跑起来也吃不消呀。

2、转移概率 ${\color{Red} p}$ 的表示

转移概率，顾名思义就是蚂蚁从当前点转移到下一个点的概率。

其中，α表示信息素相对重要的程度，β表示启发因子的相对重要程度。 $J_{k}\left ( i \right )$ 表示蚂蚁k在点i时下一步允许选择的城市集合。 $\tau _{ij}\left ( n \right )$ 表示从第i个目标点到第j个目标点间路径的信息素浓度。 $\eta_{ij}\left ( n \right )$ 表示启发因子的大小，启发因子的值等于路径长度的倒数，这里等于该点到终点距离的倒数，如下式所示。

那么，转移概率如何用来选择下一个行走的目标点呐？就用轮盘赌算法呗，关于轮盘赌算法，很简单的，不懂的可以百度~算了，还是说一下吧。

2.5、插播轮盘赌方法选择下一步目标点

在所有可选目标点的转移概率计算完毕后，使用轮盘赌算法进行最终的选择。转转盘抽奖都玩过吧？对就是这个道理！首先计算每个个体的累积概率 $q_{j}$ ，如下式：

$q_{j}$ 相当于转盘上的跨度，跨度越大的区域越容易选到，m代表下一步可选路径的数量。

之后随机生成一个(0，1)的小数，比较所有 $q_{j}$ 与r的大小，选出大于r的最小的那个 $q_{j}$ ，该 $q_{j}$ 对应的索引j即为在第n次迭代中第k只蚂蚁在第i条路径时下一步要选择的目标点 $J_{k}\left ( i \right )$ 。

3、信息素 ${\color{Red} \tau }$ 的表示

当所有蚂蚁完成1个世代的迭代后，更新各条路径上的信息素。 $\tau _{i,j}$ 表示从第i个栅格点到第j个栅格点中间路径上的信息素。初始化时 $\tau _{i,j}\left ( 0 \right )=C$ ，C为一个正常数，0表示第0次迭代，即刚开始时初始化所有路径上的信息素为一常数。 ρ为信息素挥发因子，模拟真实蚂蚁信息素在自然环境下的挥发，因此其取值范围为(0，1)。Q为一个正常数。 $\lambda$ 表示第k只蚂蚁是否走过i到j这条路，若走过则 $\lambda$ 为1，否则为0。

4、根据前述的蚁群算法的流程进行迭代运算即可找到一条当前最优路径。

三、MATLAB代码

% ACO算法实现
function [Lnum_max, num] = p2p(Alpha, Beta)
    G=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 1 1 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0 0;
        0 1 1 0 0 0 1 1 1 0 0 0 0 1 1 1 0 0 0 0;
        0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0;
        0 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0 1 1 0 0;
        0 1 1 1 0 0 1 1 1 0 0 0 0 0 0 0 0 0 0 0;
        0 1 1 1 0 0 1 1 1 0 1 1 1 1 0 0 0 1 1 0;
        0 1 1 1 0 0 0 0 0 0 1 1 1 1 0 0 0 1 1 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0;
        0 0 0 0 0 0 0 1 1 1 0 1 1 1 0 0 0 0 0 0;
        0 0 0 0 0 0 0 1 1 1 0 1 1 1 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 1 1 1 0 1 1 1 1 0;
        0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0;
        1 1 1 1 0 0 0 0 0 0 0 1 1 1 0 1 1 1 1 0;
        1 1 1 1 0 0 1 1 0 1 1 1 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 1 1 0 1 1 1 0 0 0 0 0 1 1 0;
        0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1 0;
        0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0;];
    MM=size(G,1);                  	   % 获取边长，G 地形图为01矩阵，如果为1表示障碍物

    Tau=ones(MM*MM,MM*MM);        % Tau 初始信息素矩阵
    Tau=8.*Tau;
    K=100;                       	   %迭代次数（指蚂蚁出动多少波）
    M=50;                        	   %蚂蚁个数
%     Alpha=1;                      	   % Alpha 表征信息素重要程度的参数
%     Beta=7;                       	   % Beta 表征启发式因子重要程度的参数
    Rho=0.3 ;                      	   % Rho 信息素蒸发系数
    Q=1;                               % Q 信息素增加强度系数

    minkl=inf;   % 最小正数
    mink=0;
    minl=0;
    D=G2D(G);
    N=size(D,1);               %N表示问题的规模（像素个数）
    a=1;                     %小方格像素的边长

    S=1 ;                         	   %最短路径的起始点
    E=MM*MM;                        %最短路径的目的点
    Ex=a*(mod(E,MM)-0.5);    %终止点横坐标
    if Ex==-0.5
        Ex=MM-0.5;
    end
    Ey=a*(MM+0.5-ceil(E/MM)); %终止点纵坐标 ceil向上取整

    Eta=zeros(N);             %启发式信息，取为至目标点的直线距离的倒数  N=L^2
    %以下启发式信息矩阵,每个点到终点距离的倒数
    for i=1:N
        ix=a*(mod(i,MM)-0.5);
        if ix==-0.5
            ix=MM-0.5;
        end
        iy=a*(MM+0.5-ceil(i/MM));
        if i~=E
            Eta(i)=1/((ix-Ex)^2+(iy-Ey)^2)^0.5;
        else
            Eta(i)=100;
        end
    end
    ROUTES=cell(K,M);     %用元胞结构存储每一代的每一只蚂蚁的爬行路线，K=100 迭代次数 M=50蚂蚁数量 cellplot:利用图形方式显示内容（奈斯！）
    PL=zeros(K,M);         %用矩阵存储每一代的每一只蚂蚁的爬行路线长度
    %启动K轮蚂蚁觅食活动，每轮派出M只蚂蚁
    for k=1:K
        for m=1:M
            %状态初始化
            W=S;                  %当前节点初始化为起始点 S为起始点
            Path=S;                %爬行路线初始化
            PLkm=0;               %爬行路线长度初始化
            TABUkm=ones(N);       %禁忌表初始化
            TABUkm(S)=0;          %已经在初始点了，因此要排除
            DD=D;                 %邻接矩阵初始化
            %下一步可以前往的节点
            DW=DD(W,:);
            DW1=find(DW);  % 当前点8邻域的索引，find：寻找非零元素的索引和值，返回线性索引
            for j=1:length(DW1)
                if TABUkm(DW1(j))==0
                    DW(DW1(j))=0;
                end
            end
            LJD=find(DW);  % 下一步的可选节点
            Len_LJD=length(LJD);%可选节点的个数
            %蚂蚁未遇到食物或者陷入死胡同或者觅食停止
            while W~=E&&Len_LJD>=1
                %轮盘赌法选择下一步怎么走
                PP=zeros(Len_LJD);
                for i=1:Len_LJD
                    PP(i)=(Tau(W,LJD(i))^Alpha)*((Eta(LJD(i)))^Beta);
                end
                sumpp=sum(PP);
                PP=PP/sumpp;%建立概率分布
                Pcum(1)=PP(1);
                for i=2:Len_LJD
                    Pcum(i)=Pcum(i-1)+PP(i);
                end
                Select=find(Pcum>=rand);
                to_visit=LJD(Select(1));
                %状态更新和记录
                Path=[Path,to_visit];       		 %路径增加
                PLkm=PLkm+DD(W,to_visit);    %路径长度增加
                W=to_visit;                   %蚂蚁移到下一个节点
                for kk=1:N   % 更新禁忌表
                    if TABUkm(kk)==0  % 禁忌表
                        DD(W,kk)=0;  % 双向关系
                        DD(kk,W)=0;
                    end
                end
                TABUkm(W)=0;				%已访问过的节点从禁忌表中删除
                DW=DD(W,:);
                DW1=find(DW);
                for j=1:length(DW1)
                    if TABUkm(DW1(j))==0
                        DW(j)=0;
                    end
                end
                LJD=find(DW);
                Len_LJD=length(LJD);%可选节点的个数
            end
            %记下每一代每一只蚂蚁的觅食路线和路线长度
            ROUTES{k,m}=Path;
            if Path(end)==E
                PL(k,m)=PLkm;
                if PLkm<minkl
                    mink=k;
                    minl=m;
                    minkl=PLkm;
                end
            else
                PL(k,m)=0;
            end
        end
        %更新信息素
        Delta_Tau=zeros(N,N);%更新量初始化
        for m=1:M
            if PL(k,m)
                ROUT=ROUTES{k,m};
                TS=length(ROUT)-1;%跳数
                PL_km=PL(k,m);
                for s=1:TS
                    x=ROUT(s);
                    y=ROUT(s+1);
                    Delta_Tau(x,y)=Delta_Tau(x,y)+Q/PL_km;
                    Delta_Tau(y,x)=Delta_Tau(y,x)+Q/PL_km;
                end
            end
        end
        Tau=(1-Rho).*Tau+Delta_Tau;%信息素挥发一部分，新增加一部分
    end
    
    for i=1:K
        PLK=PL(i,:);  % 第i次迭代每只蚂蚁的路程
        Nonzero=find(PLK);
        PLKPLK=PLK(Nonzero);
        if length(PLKPLK) == 0  % 该世代的蚂蚁全部没找到终点
            minPL(i) = 10000;  % 此时赋一个很大的值，代表没最优路径
            continue;
        end
        minPL(i)=min(PLKPLK);  % 本次迭代的最短路程
    end
    
    %绘图
    plotif=1;%是否绘图的控制参数
    if plotif==1 %绘收敛曲线
        minPL=zeros(K);
        for i=1:K
            PLK=PL(i,:);
            Nonzero=find(PLK);
            PLKPLK=PLK(Nonzero);
            minPL(i)=min(PLKPLK);
        end
        figure(1)
        plot(minPL);
        hold on
        grid on
        title('收敛曲线变化趋势');
        xlabel('迭代次数');
        ylabel('最小路径长度'); %绘爬行图
        figure(2)
        axis([0,MM,0,MM])
        for i=1:MM
            for j=1:MM
                if G(i,j)==1
                    x1=j-1;y1=MM-i;
                    x2=j;y2=MM-i;
                    x3=j;y3=MM-i+1;
                    x4=j-1;y4=MM-i+1;
                    fill([x1,x2,x3,x4],[y1,y2,y3,y4],[0.2,0.2,0.2]);
                    hold on
                else
                    x1=j-1;y1=MM-i;
                    x2=j;y2=MM-i;
                    x3=j;y3=MM-i+1;
                    x4=j-1;y4=MM-i+1;
                    fill([x1,x2,x3,x4],[y1,y2,y3,y4],[1,1,1]);
                    hold on
                end
            end
        end
        hold on
        title('本次迭代的最优路径');
        xlabel('x');
        ylabel('y');
        ROUT=ROUTES{mink,minl};
        LENROUT=length(ROUT);
        Rx=ROUT;
        Ry=ROUT;
        for ii=1:LENROUT
            Rx(ii)=a*(mod(ROUT(ii),MM)-0.5);
            if Rx(ii)==-0.5
                Rx(ii)=MM-0.5;
            end
            Ry(ii)=a*(MM+0.5-ceil(ROUT(ii)/MM));
        end
        plot(Rx,Ry)
    end
end

% G2D函数实现
function D=G2D(G)
    L=size(G,1);
    D=zeros(L*L,L*L);
    for i=1:L
        for j=1:L
            if G(i,j)==0
                for m=1:L
                    for n=1:L
                        if G(m,n)==0
                            im=abs(i-m);
                            jn=abs(j-n);
                            if im+jn==1||(im==1&&jn==1)  % 8邻域内
                                D((i-1)*L+j,(m-1)*L+n)=(im+jn)^0.5;  % D矩阵内每一行存储是一个点的8邻域内的距离
                            end
                        end
                    end
                end
            end
        end
    end
end

四、仿真结果