Ant colony algorithm is an optimization algorithm to simulate the foraging behavior of ant colonies. The ants scatter pheromones throughout the foraging process, and the ants use the amount of pheromone they perceive to decide which grid to choose.

In the initial stage, since there are no pheromones on the ground, the ant colony’s walking path is random, and the ants will continuously release pheromones to mark their walking path. Over time, several ants find food, and there are several routes from the burrow to the food. Because the behavioral trajectories of ants are randomly distributed, the number of ants on the short path is denser than those on the long path per unit time, and the pheromone concentration left by the short path is higher. This provides a strong direction for the ants behind them, and more and more ants flock to the shortest route. For an individual ant, it’s not looking for the shortest path, it’s just choosing according to probability; For the whole ant colony system, they achieve the objective effect of finding the optimal path. Assuming that the total number of ants in the ant colony is M, each ant moves in the grid environment, and selects the next grid according to the state transition rule. Assuming that ant K is located in grid I at time T, the probability of ant K selecting the next grid J is:

                                                        

In Formula (1), V represents the set of the next grid that ant K can choose; Alpha is the heuristic factor of pheromone concentration. The larger Alpha is, the more ant K tends to choose the path traveled by most ants. Beta represents the expectation heuristic factor, reflecting the effect of visibility information on ant’s selection of next position. The larger the Beta value is, the more ant K tends to choose the grid close to the target point, and the more inclined it is to move towards the visibility range. Is the pheromone concentration on the path (I,j) at time t; Represents the heuristic information on the path (I,j) at time t, which is defined as:

      

The core part of ant colony algorithm is to simulate the behavior of ant colony transition probability selection and calculate the transition probability by using pheromone and heuristic function value. In the process of ant state transfer, the reciprocal of the distance between the node and the target point is used as heuristic information, which is not conducive to the pre-avoidance of obstacles. And path planning in complex environment, the ant colony algorithm search in a large space, at the beginning of the optimized path pheromone concentration is small, positive feedback is not obvious especially in the process of random solution of “blind search” to produce a large number of local cross paths, reduce the operating efficiency of the ant colony algorithm, and easy to fall into local optimum, the search is to a certain degree, It is easy to stagnate, and the solutions found by all individuals are exactly the same, so further search is not conducive to finding better solutions.

t=[4 5;16 25.8;10 45;20 55;30 65;35 55;29 31;37 26;47 27;30 31.3;31 17;14 7;35.6 13.8;26.7 22.5;21 39;38 42;5 26;28 53;20 13;10 60;26 31;54 38;7 58;12 36;30 2] %24个点,第25个点事origin
save t.mat t
load t.mat
value=[1 1 1 2 3 2 1 3 3 2 2 2 2 2 1 2 3 3 1 1 2 1 1 1]; %24个目标的价值
value=value/100;
time=zeros(1,25); %侦察UAV时间数组,里面放的是飞机走的航程,除以速度便是时间,设速度为‘1’
attacktime=zeros(1,25); %打击UAV时间数组
%把侦察过的任务对应无人机走过的航程存到该任务的一个矩阵里,当做时间,然后打击任务如果选定某任务,check一下时间是否合格,合格的话,可以打击,并存入禁忌表,不合格的话,选次概率的
 
%注:目的是把所有目标执行完所有任务,所以每次迭代最后所有无人机收获的总价值都一样,都是所有目标的价值之和;所以本程序考虑优先执行价值大的目标,防止无人机飞很久、打很久后攻打效率变低
%的情况出现
 
%%计算城市间相互距离
n=size(t,1);
D=zeros(n,n);
for i=1:n
    for j=1:n
        if i~=j
            D(i,j)=sqrt(sum((t(i,:)-t(j,:)).^2));
        else
            D(i,j)=1e-2;
        end
    end
end
%%初始化参数
m=10;         %蚂蚁个数
alpha=1;      %信息素重要程度因子
beta=1;       %启发函数重要程度因子
gama=2;
rho=0.3;      %信息素挥发因子
Q=1.0;          %总量
eta=1./D;     %启发函数
tau=ones(n,n)+7.1192e-005;%信息素矩阵
iter=1;       %迭代次数初始值
iter_max=80;  %迭代次数最大值
length_best=zeros(iter_max,1);%每次迭代最佳路径长度(应该是一次比一次小)
length_ave=zeros(iter_max,1); %每次迭代路径平均长度 
%%迭代寻找最佳路径
while iter<=iter_max
    whta=cell(8,1);
    lieend=zeros(8,1);
    for zu=1:8
    city_index=1:25;       %城市来标号 
    table=[];
    start=zeros(4,1);
        temp=randperm(24);
        for i=1:4
        start(i)=temp(i);
        end
    table(:,1)=start;
    j=2;
 while (j<=30)
        for i=1:4
            if i==1 %UAV1只负责“侦察”任务
                if table(1,(j-1))~=25
                    table1=table(1,:);
                    table1=[table1;table(3:4,:)];
                    tabu1=table1(:); %UAV1的禁忌表出来了 %25如果也在tabu1里的话,那么
            allow_index1=~ismember(city_index,tabu1);  %【走过的变成0,能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意:allow_index与city_index同维】
            allow1=city_index(allow_index1);  %把还能走的序号摘出来了(待访问的城市集合)
            P1=allow1;
            %计算城市的转移概率
            if numel(allow1)~=0
              for k=1:max(size(allow1))-1
               P1(k)=(tau(table(1,(j-1)),allow1(k))^alpha)*(eta(table(1,(j-1)),allow1(k))^beta)*10000+7.1192e-004;
              end
            P1(max(size(allow1)))=7.1192e-005;
            P1=P1/sum(P1);
            [d1,ind1]=sort(P1,2,'descend');%从大到小排序是d1,对应的原序号是ind1
            target1=allow1(ind1(1));
            %轮盘赌法选择下一个城市
            %pc1=cumsum(P1);  % (p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....)【p1<->allow(1)  p2<->allow(2) ...】
            %target_index1=find(pc1>=rand); 
            %target1=allow1(target_index1(1));  %这次返回的是allow数组中城市的真正序号
            table(1,j)=target1;  %把选好这个点放到路径表里面
            rr=D(25,table(1,1));
            time(table(1,1))=rr;
            if j>2
            for c=2:(j-1)
                rr=rr+D(table(1,c-1),table(1,c));
            end 
            end
            rrr=rr+D(table(1,j-1),target1);%rrr就是UAV1到该点时走过的航程
            time(target1)=rrr;
            else
                table(1,j)=25;
            end
                end
                 if table(1,(j-1))==25
                    table(1,j)=25;
                 end
            end          
            if i==2 %UAV2只负责“打击”任务
                if (table(2,(j-1))~=25)
                table(2,1)=table(1,1); %设定它第一次打击的是UAV1侦察过的目标
                ta2=table(1:(4*(j-1)+1)); %当前元素之前所有的元素
                tabu21=[];
                tabu22=[];
                tabu2=[];
                for y=1:24
                    if sum(ta2==y)==2
                        tabu21=[tabu21;y];
                    end
                end                       %出现过两次的放在tabu21里
                tabu22=setdiff(1:24,ta2); %一次都没出现的放在tabu22里
                tabu2=[tabu21',tabu22];   %tabu2出来了
            allow_index2=~ismember(city_index,tabu2);  %【走过的变成0,能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意:allow_index与city_index同维】
            allow2=city_index(allow_index2);  %把还能走的序号摘出来了(待访问的城市集合)
            P2=allow2;
            %计算城市的转移概率
           for k=1:(length(allow2)-1)
               P2(k)=tau(table(2,(j-1)),allow2(k))*eta(table(2,(j-1)),allow2(k))*value(allow2(k))*10000;
           end
           P2(max(size(allow2)))=7.1192e-005;
            P2=P2/sum(P2);
            [d2,ind2]=sort(P2,2,'descend');%从大到小排序是d1,对应的原序号是ind1
            target2=allow2(ind2(1)); %target2=d1(1);
            %轮盘赌法选择下一个城市
            %pc2=cumsum(P2);  % (p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....)【p1<->allow(1)  p2<->allow(2) ...】
            %target_index2=find(pc2>=rand);  %选中那个概率较大的选中的点,返回的是allow数组中的序号
            %target2=allow2(target_index2(1));  %这次返回的是allow数组中城市的真正序号
            %table(2,j)=target2;  %把选好这个点放到路径表里面
            oo=D(25,table(2,1));
            attacktime(table(2,1))=oo;
            if j>2
            for c=2:(j-1)
                oo=oo+D(table(2,c-1),table(2,c));
            end 
            end
            ooo=oo+D(table(2,j-1),target2);%ooo就是UAV2到该点时走过的航程
            if numel(d2)>5
            u=2;
            while (ooo>time(target2)+20 & u<6)
                 target2=allow2(ind2(u));
                 ooo=oo+D(table(2,(j-1)),target2);
                 u=u+1;
            end
            end
            table(2,j)=target2;
            attacktime(target2)=ooo;
                end
                if table(2,(j-1))==25
                    table(2,j)=25;
                end
            end
            if i==3 %UAV3是“察打”任务
                if table(3,(j-1))~=25
                    ta3=table(1:(4*(j-1)+2));
                    tabu3=[];
                    tabu3c=[];
                for y=1:24
                    if sum(ta3==y)==2
                        tabu3=[tabu3;y];
                    end
                end    %出现两次的放在tabu3里   
                  for y=1:24
                    if sum(ta3==y)==1
                        tabu3c=[tabu3c;y];
                    end
                end %tabu3c是待打的任务,已侦查完的任务
            allow_index3=~ismember(city_index,tabu3);  %【走过的变成0,能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意:allow_index与city_index同维】
            allow3=city_index(allow_index3);  %把还能走的序号摘出来了(待访问的城市集合)
            P3=allow3;
            %计算城市的转移概率
           for k=1:(length(allow3)-1)
               %if ismember(allow3(k),tabu3c)==1     
               h=table(3,(j-1))
               P3(k)=(tau(table(3,j-1),allow3(k))^alpha)*(eta(table(3,(j-1)),allow3(k))^beta)*value(allow3(k))*10000+7.1192e-005;%这是要打的,需要价值
               %else
               %P3(k)=(tau(table(3,(j-1)),allow3(k))^alpha)*(eta(table(3,(j-1)),allow3(k))^beta)*100+7.1192e-005;%这些是待侦察的,没有价值
               %end
           end
            P3(max(size(allow3)))=7.1192e-009;
            P3=P3/sum(P3);
            [d3,ind3]=sort(P3,2,'descend');%从大到小排序是d1,对应的原序号是ind1
            target3=allow3(ind3(1));
            %轮盘赌法选择下一个城市
            %pc3=cumsum(P3);  % (p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....)【p1<->allow(1)  p2<->allow(2) ...】
            %target_index3=find(pc3>=rand);  %选中那个概率较大的选中的点,返回的是allow数组中的序号
            %target3=allow3(target_index3(1));  %这次返回的是allow数组中城市的真正序号
            %table(3,j)=target3;  %把选好这个点放到路径表里面
            ww=D(25,table(3,1));
            time(table(3,1))=ww;
            if j>2
            for c=2:(j-1)
                ww=ww+D(table(3,c-1),table(3,c));
            end 
            end
            www=ww+D(table(3,j-1),target3);%www就是UAV3到该点时走过的航程
            if ismember(target3,tabu3c)==0 %侦察任务
                time(target3)=www;
                table(3,j)=target3;
            else %打击任务
                attacktime(target3)=www;
                if numel(d3)>5
                u=2;
            while (www>time(target3)+20 & u<6)
                 target3=allow3(ind3(u));
                 www=ww+D(table(3,(j-1)),target3);
                 u=u+1;
            end
                end
            attacktime(target3)=www;
            table(3,j)=target3;%www<time(target3)+10 说明此打击任务合理
            end 
                end
                if table(3,(j-1))==25
                    table(3,j)=25;
                end
            end
            if i==4 %UAV4是“察打”任务
                if table(4,(j-1))~=25
                    ta4=table(1:(4*(j-1)+3));
                    tabu4=[];
                    tabu4c=[];
                for y=1:24
                    if sum(ta4==y)==2
                        tabu4=[tabu4;y];
                    end
                end    %出现两次的放在tabu4里、、可以把已经侦察过的放在tabu4c中,(即出现过一次的),如果选到的是在tabu4'中的说明是要打击的然后算一下它的航程,再和侦察路径比较 
                for y=1:24
                    if sum(ta4==y)==1
                        tabu4c=[tabu4c;y];
                    end
                end 
            allow_index4=~ismember(city_index,tabu4);  %【走过的变成0,能走的为1】【若tabu=(1 4)则allow_index=(0 1 1 0 1 1 1...)】【注意:allow_index与city_index同维】
            allow4=city_index(allow_index4);  %把还能走的序号摘出来了(待访问的城市集合)
            P4=allow4;
            %计算城市的转移概率
           for k=1:(max(size(allow4))-1)
               %if ismember(allow4(k),tabu4c)==1
               sxx=table(4,(j-1))
               P4(k)=(tau(table(4,(j-1)),allow4(k))^alpha)*(eta(table(4,(j-1)),allow4(k))^beta)*value(allow4(k))*10000+7.1192e-005;
               %else
               %P4(k)=(tau(table(4,(j-1)),allow4(k))^alpha)*(eta(table(4,(j-1)),allow4(k))^beta)*100+7.1192e-005;
               %end
           end
           P4(max(size(allow4)))=7.1192e-009;
            P4=P4/sum(P4);
            [d4,ind4]=sort(P4,2,'descend');%从大到小排序是d1,对应的原序号是ind1
            target4=allow4(ind4(1));
            %轮盘赌法选择下一个城市
            %pc4=cumsum(P4);  % (p1 p1+p2 p1+p2+p3 p1+p2+p3+p4 ....)【p1<->allow(1)  p2<->allow(2) ...】
            %target_index4=find(pc4>=rand);  %选中那个概率较大的选中的点,返回的是allow数组中的序号
            %target4=allow4(target_index4(1));  %这次返回的是allow数组中城市的真正序号
            %table(4,j)=target4;  %把选好这个点放到路径表里面
            qq=D(25,table(4,1));
            time(table(4,1))=qq;
            if j>2
            for c=2:(j-1)
                qq=qq+D(table(4,c-1),table(4,c));
            end 
            end
            qqq=qq+D(table(4,j-1),target4);%www就是UAV3到该点时走过的航程
            if ismember(target4,tabu4c)==0 %侦察任务
                time(target4)=qqq;
                table(4,j)=target4;
            else %打击任务
                attacktime(target4)=qqq;
                if numel(d4)>5
                u=2;
            while (qqq>time(target4)+20 & u<6)
                 target4=allow4(ind4(u));
                 qqq=qq+D(table(4,j-1),target4);
                 u=u+1;
            end
                end
            attacktime(target4)=qqq;
            table(4,j)=target4;%www<time(target4)+10 说明此打击任务合理
            end                
                end
                if table(4,(j-1))==25
                    table(4,j)=25;
                end
              end
        end %一列结束
Copy the code