基于蚁群算法的路径规划_蚁群算法路径规划

蚁群算法（Ant Colony Algorithm, ACO）于1991年首次提出，该算法模拟了自然界中蚂蚁的觅食行为。蚂蚁在寻找食物源时，会在其经过的路径上释放一种信息素，并能够感知其它蚂蚁释放的信息素。信息素浓度的大小表征路径的远近，信息素浓度越高，表示对应的路径距离越短。通常，蚂蚁会以较大的概率优先选择信息素浓度较高的路径，并释放一定量的信息素，以增强该条路径上的信息素浓度，这样，会形成一个正反馈。最终，蚂蚁能够找到一条从巢穴到食物源的最佳路径，即距离最短。

用蚂蚁的行走路径表示待优化问题的可行解， 整个蚂蚁群体的所有路径构成待优化问题的解空间。

路径较短的蚂蚁释放的信息素量较多， 随着时间的推进， 较短的路径上累积的信息素浓度逐渐增高， 选择该路径的蚂蚁个数也愈来愈多。

最终， 整个蚂蚁会在正反馈的作用下集中到最佳的路径上， 此时对应的便是待优化问题的最优解。

具体步骤如下：

设整个蚂蚊群体中蚂蚊的数量为m, 路径节点的数量为n, 节点i与节点 j 之间的相互距离为dij ( i , j = 1 , 2 , … , n ) , t时刻节点i与节点j的连接路径上的信息素浓度为τij (t)，初始时刻, 各个节点间连接路径上的信息素浓度相同。

蚂蚁k( k =1,2,…,m)根据各个节点间连接路径上的信息素浓度决定其下一个访问节点, 设Pijk (t)表示t时刻蚂蚊 k从节点 i转移到节点 j的概率, 其计算公式如下:

ηij(t)为启发函数, ηij(t)=1/dij，表示蚂蚊从节点i转移到节点j的期望程度。

allowk(k=1,2,…,m)为蚂蚁k待访问节点的集合。开始时,allowk中有n-1个元素,即包括除了蚂蚁k出发节点的其它所有节点。随着时间的推进, allowk中的元素不断减少, 直至为空, 即表示所有的节点均访问完毕。

α为信息素重要程度因子,其值越大, 蚂蚁选择之前走过的路径可能性就越大,搜索路径的随机性减弱,其值越小,蚁群搜索范围就会减少,容易陷入局部最优。一般取值范围为[0,5]。

β为启发函数重要程度因子,其值越大,表示启发函数在转移中的作用越大, 即蚂蚊会以较大的摡率转移到距离短的节点，蚁群就越容易选择局部较短路径,这时算法的收敛速度是加快了，但是随机性却不高，容易得到局部的相对最优。一般取值范围为[0,5]。

计算完节点间的转移概率后，采用轮盘赌方法选择下一个待访问的节点。依据轮盘赌法来选择下一个待访问的节点，而不是直接按概率大小选择，是因为这样可以扩大搜索范围，进而寻找全局最优，避免陷入局部最优。首先计算每个个体的累积概率qj ，之后随机生成一个(0，1)的小数r，比较所有qj与r的大小，选出大于r的最小的那个qj，该qj对应的索引j即为第k只蚂蚁在第i条路径时下一步要选择的目标点。

蚂蚁信息素更新的模型包括蚁周模型（Ant-Cycle模型）、蚁量模型（Ant-Quantity模型）、蚁密模型（Ant-Density模型）等。

区别：

蚁周模型利用的是全局信息，即蚂蚁完成一个循环后更新所有路径上的信息素。

蚁量和蚁密模型利用的是局部信息，即蚂蚁完成一步后更新路径上的信息素。

在蚂蚁释放信息素的同时，各个节点间连接路径上的信息素逐渐消失,设参数ρ(0<ρ<1，一般取值为0.1～0.99 )表示信息素的挥发程度。当所有的蚂蚁完成一次循环后，各个节点间链接路径上的信息素浓度需进行更新，计算公式为

τijt+1=1-ρτijt+Δτij

Δτij=k=1nΔτijk​

其中，Δτijk表示第k只蚂蚁在节点i与节点j连接路径上释放的信息素浓度；Δτij表示所有蚂蚁在节点i与节点j连接路径上释放的信息素浓度之和。

算法流程图如下：

在实验中，设定的蚂蚁数量为50，迭代200次。在100*100的地图上随机生成地图（通过MATLAB的binornd函数实现）。从左上角顶点处（起点）走到右下角顶点处（终点），结果如下。

实验代码如下

clc

clear

G=binornd(1,0.05,100,100);

G(1,1)=0;

G(100,100)=0;

MM=size(G,1);                 %G地形图为01矩阵，如果为1表示障碍物

Tau=ones(MM*MM,MM*MM);        %Tau初始信息素矩阵 大小要与Alpha系数匹配

%Tau=8.*Tau;

K=200;                        %迭代次数（指蚂蚁出动多少波）

M=50;                         %蚂蚁个数

S=1;                          %最短路径的起始点

E=MM*MM;                      %最短路径的目的点

Alpha=1;                      %Alpha表征信息素重要程度的参数

Beta=7;                       %Beta表征启发式因子重要程度的参数

Rho=0.4;                      %Rho信息素蒸发系数

Q=1;                         %Q信息素增加强度系数

minkl=inf;                    %minkl表示当前最短路径长度

mink=0;                       %当前完成最短路径为第几次迭代

minl=0;                       %当前完成最短路径为第只蚂蚁

D=G2D(G);                     %每个栅格至各自邻域无障碍栅格的代价值

N=size(D,1);                  %N表示问题的规模（像素个数）

a=1;                          %小方格像素的边长

Ex=a*(mod(E,MM)-0.5);         %目的点横坐标

if Ex==-0.5

    Ex=MM-0.5;

end

Ey=a*(MM+0.5-ceil(E/MM));     %目的点纵坐标

Eta=zeros(1,N);               %启发式信息矩阵，记录启发式信息

%以下开始建立启发式信息矩阵

for i=1:N                     %每个栅格的索引号，一共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次迭代中的第M只蚂蚁的爬行路线

PL=zeros(K,M);        %用矩阵存储第K次迭代中的第M只蚂蚁爬行路线的长度

%开始K次迭代，每轮派出M只蚂蚁，开始寻找路径

for k=1:K               %第K次迭代

    for m=1:M           %第M只蚂蚁

        %状态初始化

        W=S;                  %将当前节点初始化为起始点

        Path=S;               %爬行路线初始化

        PLkm=0;               %爬行路线长度初始化

        TABUkm=ones(N);       %禁忌表初始化为1，禁忌表记录走过的位置，将走过的位置由1变0

        TABUkm(S)=0;          %将禁忌表起点位置置为0

        DD=D;                 %邻栅格点初始化

        %下一步可以前往的栅格点

        DW=DD(W,:);   %取G2D矩阵中以当前点为局部起始点的一行

        DW1=find(DW); %DW1矩阵存储该行中所有无障碍相邻栅格点（元素不为0）的索引位置

        for j=1:length(DW1)

            if TABUkm(DW1(j))==0  %判断TABUkm禁忌表中，该位置是否为之前走过的点

                DW(DW1(j))=0;     %删除DW中所有之前已经走过的相邻栅格点

            end

        end %现在DW中为当前栅格点可以选择的所有相邻栅格点了

        %计算各可选择邻节点的选择概率

        LJD=find(DW);       %LJD记录未走过的点的索引号，即下一步可选的点

        Len_LJD=length(LJD);%可选点的个数

        %蚂蚁未遇到食物或者陷入死胡同或者觅食停止

        while W~=E&&Len_LJD>=1 %起始点不等于终止点，且可选节点个数大于等于1

            %转轮赌法选择下一步怎么走

            PP=zeros(Len_LJD);

            for i=1:Len_LJD

                PP(i)=(Tau(W,LJD(i))^Alpha)*((Eta(LJD(i)))^Beta);

                %PP(i)=(Tau(W,LJD(i))^Alpha)*((D(W,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;                  %蚂蚁移到下一个点

            %对应禁忌表更新D中可选择的相邻栅格点

            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(DW1(j))=0;  

                end

            end

            LJD=find(DW);

            Len_LJD=length(LJD);%可选节点的个数

        end                 %本次迭代的当前蚂蚁寻路完毕

        %记下每一代每一只蚂蚁的觅食路线和路线长度

        ROUTES{k,m}=Path; %记录本次迭代中的当前蚂蚁的行走路线

        if Path(end)==E   %判断本只蚂蚁寻找路径的最后一个节点是否为终点

            PL(k,m)=PLkm; %若该蚂蚁到达终点，将本次路线长度放到PL的第k行m列

            if PLkm<minkl %若本次路径长度<当前已知的最短路径长度

                mink=k;   %记录完成本次最短路径的迭代次数

                minl=m;   %记录完成本次最短路程的哪只蚂蚁

                minkl=PLkm; %记录本次最短路线的长度

            end

        else

            PL(k,m)=0; %若该蚂蚁没有到达终点，则本次路径长度为0

        end

    end %返回进行下一只蚂蚁的寻路

    %更新信息素

    Delta_Tau=zeros(N,N);%初始化信息素增量

    for m=1:M

        if PL(k,m) 

            ROUT=ROUTES{k,m}; %ROUT取本次迭代的所有蚂蚁的行走路线

            TS=length(ROUT)-1;

            PL_km=PL(k,m); %PL_km取本次迭代当前蚂蚁的路径长度

            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

%绘图

plotif=1;%是否绘图的控制参数

if plotif==1 %绘收敛曲线

    minPL=zeros(K);

    for i=1:K %选出每次迭代的最短路径

        PLK=PL(i,:); %取第i次迭代的所有路径长度

        Nonzero=find(PLK); %提取本次迭代路径长度不为0的索引号存储至Nonzero

        if isempty(Nonzero)

            minPL(i)=0;

        else

            PLKPLK=PLK(Nonzero);

            minPL(i)=min(PLKPLK); %选出本次迭代最短路径存储至minPL的相应迭代位置

        end

    end

    %绘制“收敛曲线变化趋势”图，将每次迭代的最短路径放在图中表示

    FFF=[];

    WWJ=minPL(:,1);

    for i=1:length(WWJ)

        if WWJ(i)~=0

            FFF=[FFF;[i,WWJ(i)]];

        end

    end

    figure(1)

    plot(FFF(:,1),FFF(:,2));

    axis([0 200 0 1000]);

    hold on

    grid on

    title('收敛曲线变化趋势');

    xlabel('迭代次数');

    ylabel('最小路径长度');

    %绘爬行图

    figure(2)

    axis([0,MM,0,MM]) %设置图的横纵坐标，MM为地图矩阵的行数或列数

    for i=1:MM

        for j=1:MM

            if G(i,j)==1 %1是黑色代表障碍，0为白色无障碍

                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]); %将1234点所围成的图形进行黑色填充

                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]); %将1234点所围成的图形进行白色填充

                hold on

                set(gca,'YTickLabel',[100 90 80 70 60 50 40 30 20 10 0]); %使得地图的矩阵的行列与正常坐标轴的行列一致

                %set(gca,'YTickLabel');

            end

        end

    end

    hold on

    title('最优路径');

    xlabel('坐标x');

    ylabel('坐标y');

    ROUT=ROUTES{mink,minl};

    %ROUT取最短行走路线，mink为该路线的迭代号,minl为该路线的蚂蚁号

    LENROUT=length(ROUT);

    %Rx与Ry中分别存储具体的该路线

    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