蚁群算法用于航路规划的matlab简单实现
提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档
前言
本篇博客的内容其实来自学校的一门课程，实现部分代码编写之后尝试把它写到博客上练练手。
提示：以下是本篇文章正文内容，下面案例可供参考
一、蚁群算法
蚁群算法是在20世纪90年代由澳大利亚学者Marco Dorigo等人通过观察蚁群觅食的过程，发现众多蚂蚁在寻找食物的过程中，总能找到一条从蚂蚁巢穴到食物源之间的最短路径。随后他们在蚂蚁巢穴到食物源之间设置了一个障碍，一段时间以后发现蚂蚁又重新走出了一条到食物源最短的路径。通过对这种现象的不断研究，最后提出了蚁群算法。蚁群算法在解决旅行商问题（即TSP问题）时，取得了比较理想的结果。
图：蚁群搜索路径示例
图：蚁群算法用于多旅行商问题
二、背景
现实世界中，（初始状态）每个蚂蚁将在一定范围内随机游走，并通过在所经过路径上选留信息素的方式，不断把相关的食物信息反馈给蚁群。如果其他蚂蚁发现了那条路径，那么这些妈蚁很可能不再保持原来的随机游走，而跟随这条路径**（一条路径上的信息素越多，其他蚂蚁选择这条路径的可能性越大）**。如果由该路径他们最终发现了食物，那么他们将在这条路径上继续增加信息素。由于这条路径上信息素的不断增加，最终所有蚂蚁将找到食物。
无人机的航路规划需要考虑到空间中的障碍和威胁，在规划的过程中，不一定能够绕过所有的威胁源，在此情况下，将障碍设定为禁忌移动方向并将航路的威胁作为代价对航路进行规划有着一定的意义。
三、算法原理
（1）蚁群算法试验中，在各个蚂蚁在没有事先告诉它们食物在什么地方的前提下开始寻找；
其次，要让蚂蚁找到食物，就需要让它们遍历空间上的所有点；
（2）再次，如果要让蚂蚁找到最短的路找食物。当一只找到食物以后，它会向环境释放一种信息素，吸引其他的蚂蚁过来，这样越来越多的蚂蚁会找到食物。
（3）但有些蚂蚁并没有像其它蚂蚁一样总重复同样的路，它们会另辟蹊径，**如果令开辟的道路比原来的其他道路更短，那么更多的蚂蚁被吸引到这条较短的路上来。**最后，经过一段时间运行，可能会出现一条最短的路径被大多数蚂蚁重复着。
从数学角度来看, 蚂蚁构建的这一连接巢穴和食物源的路网组成了最小树, 使得蚂蚁在道路上为搜寻食物和把食物搬运到巢穴所消耗的能量最小。图论的知识已为我们提供了计算最小树的算法, 但蚂蚁采用了散布生物信息激素的方法来协同搜索最小树 。这一全局最优的道路结构会随着不同个体在路网中散布生物信息激素并感知其它蚂蚁留下的气味而逐步浮出水面”
蚁群算法的结构一般分为（1）状态转移规则、（2）路径代价计算和（3）信息素和启发信息的更新。
（1）状态转移规则：
每只蚂蚁在经过一个节点时，除到达终点附近以外该蚂蚁会发现其他的可能探索方向，在每次选取时蚂蚁并不是以等概率随机选取的，而是遵循一定的规则，先排除掉禁忌的方向，后根据空间中留下的信息素为不同方向赋予不同的“权重”或者概率进行带有概率的随机选取。对于人工蚂蚁来说，由于具有路线记忆能力，还可以将近期最优的路径和历史最优的路径代价作为启发信息加入概率规则中。
式中:Pk( r, s) 表示蚂蚁 k 从节点r 转移到节点s 的概率; τ(r,s)表示蚂蚁存储在边 V( r, s)上的生物信息激素强度;η( r, s)表示节点 s 相对于节点r 的可见性, η( r, s) =1/ crs, crs表示边 V( r, s)的代价 ;J k( r)是第 k 个蚂蚁由节点 r 可以到达的所有可行节点集合, 这些节点均是节点 r 的相邻节点, 而且它们比节点 r 更接近目标点;α、β 为控制参数, 确定生物信息激素和可见性的相对重要性 。蚂蚁从状态 r 转移到状态s 所选可行节点的概率会随着生物信息激素强度的增大而增大, 随着通路代价的增大而减少。
本程序加入了一定的转弯角度限制，上一段和下一段之间方向的夹角绝对值不能超过45度。
（2）路径代价计算
式中:L 为航路的长度 ;W 为广义代价函数 ;wt 为航路的威胁代价 ;wf 为航路的油耗代价 ;系数 k 表示根据任务安排, 航路制订人员在制订航路过程中所做出的倾向性选择 。油耗代价是航程的函数 。威胁代价与无人机的可探测性指标相关联, 而可探测性指标是根据无人机的可探测概率计算的。wi 为第i 条边的广义代价;wt, i为第i 条边的威胁代价;wf, i为第i 条边的油耗代价，以距离估算。
本程序简单地认为威胁源存在一定的范围且威胁强度反比于无人机与威胁源之间的距离（和文献[4]中不一样，属于是编写的时候疏忽了），并在简化计算中将边的威胁代价简单地定义为路径段重点的威胁代价并赋予在终点的位置上。
（3）信息素和启发信息更新规则
一旦所有蚂蚁完成了各自候选航路的选择过程(找到一条航路规划问题的可行解) , 必须对各边上的生物信息激素作一次全面的修正, 修正规则如下 :
Wk是蚂蚁k选择的航路的广义代价，而We是当前最小的航路代价，m是蚂蚁数。生物信息激素修正的目的是分配更多的生物信息激素到具有更小威胁代价航路的边上, 这类似于增强学习格式中的算子，如遗传算法中的选择算子。这个修正规则不仅存储生物信息激素, 还适当地蒸发它们。修正规则不是由个别蚂蚁来实现, 而是通过图的边来存储, 起到了一个分布式长期记忆的效果。
本程序还增加了信息素扩散的规则，每次迭代更新后同一点的信息素除边缘点外由该点和周围邻近的8个点加权平均获得。
参考自：
[1]https://blog.csdn.net/weixin_50570276/article/details/109449635
[2]https://blog.csdn.net/m0_37570854/article/details/83715944
[3]https://blog.csdn.net/qq_33829154/article/details/85258615
[4]柳长安,李为吉,王和平.基于蚁群算法的无人机航路规划[J].空军工程大学学报(自然科学版),2004(02):9-12.
四、算法步骤
五、程序实现
主函数
代码如下（示例）：
clc;
clear all;
num_iteration=200;
time_up=20;
alpha=1;
beta=1.1;
e =2;
rou=0.5;
Q=15;
radius =1;
j_ant=10;
% 地图空间20*20
x_max=40;%最小默认为0
y_max=40;%最小默认为0
map=[];
for x=0:x_max
for y=0:y_max
map(:,size(map,2)+1)=[x;y];
end
end
% 蚁群初始位置
start1=[3/4;1/4]*x_max;
start2=[3/4+0.1;1/4+0.1]*x_max;
start3=[3/4+0.2;1/4]*x_max;
% 目标位置
x=[];y=[];
target1=[0.3;0.75]*x_max;
target2=[0.3;0.75]*x_max;
target3=[0.3;0.75]*x_max;
rectangle('Position',[0,0,x_max,y_max],'Linewidth',1,'LineStyle','-','EdgeColor','b')
scatter(target(1),target(2),'filled')
axis([0 x_max 0 y_max])
for theta=linspace(0,2*pi,50)
x(end+1)=target(1)+radius*cos(theta);
y(end+1)=target(2)+radius*sin(theta);
end
plot(x,y,'b');
grid on;
hold on;
%威胁空间(简化）
threat=[
0.2,0.4,0.5,0.5,5/8
0.7,0.6,0.8,5/8,5/8
]*x_max;
thr_space=zeros(x_max+1,y_max+1);
for i=1:size(threat,2)
for x=0:x_max
for y=0:y_max
if norm([x;y]-threat(:,i))<5
thr_space(x,y)=thr_space(x,y)+1./(0.1+norm([x;y]-threat(:,i)) );
end
end
end
x1=[];
y1=[];
for theta=linspace(0,2*pi,50)
x1(end+1)=threat(1,i)+2*cos(theta);
y1(end+1)=threat(2,i)+2*sin(theta);
end
plot(x1,y1,'r');
hold on;
axis([0 x_max 0 y_max]);
end
%初始提示，加快收敛
tao=zeros(x_max+1,y_max+1);%[0,20],21个
deltax=1;
deltay=1;
if start(1)>target(1)
deltax=-1;
end
if start(2)>target(2)
deltay=-1;
end
for i=start(1):deltax:target(1)
%disp(i)
for j= start(2):deltay:target(2)
% disp(j)
if i==start(1) || j==start(2) || i==target(1) || j==target(2)
tao(i,j)=0.02; % 给个提示加快速度
else
tao(i,j)=0.02;
end
end
end
% 蚂蚁位置更新
direc=[
-1 0 1 0 1 -1 1 -1
0 -1 0 1 1 -1 -1 1
];
eta=0.01*ones(x_max+1,y_max+1);
hist_best_score=[];
%---------------迭代计算开始-------------------
for num =1:num_iteration % 迭代次数限制
% 蚂蚁路径
dots={};
% route={};
for k=1:j_ant% 蚁群规模
dots{k}=start;
% route{k}=[];
end
%{
update_log-------------------------
取消了每次搜索的时间限制，
增加了初始提示路径、不走回走过的点，和消灭被困蚂蚁路径的机制
-----------------------------------
%}
%位置更新
for k=1:length(dots) % j_ant Ants
while norm(dots{k}(:,end)-target)>radius %距离小于半径为捕获目标
%路径禁忌判定
permission=[];
for di =1:8
if dots{k}(1,end)+direc(1,di)>=0 && dots{k}(2,end)+direc(2,di)>=0 && dots{k}(1,end)+direc(1,di)<=x_max && dots{k}(2,end)+direc(2,di)<=y_max;
%禁止超出范围
permission(:,end+1)=direc(:,di);
end
% 转弯角限制（余弦定理）
if size(dots{k},2)>=2
if direc(:,di)'*(dots{k}(:,end)-dots{k}(:,end-1))/(norm(direc(:,di))*norm(dots{k}(:,end)-dots{k}(:,end-1)))<1/sqrt(2)
if ~isempty(permission)
permission(:,end)=[];
end
end
end
% 这一段注注释的情况下是取消了每只蚂蚁走过同一点的限制，如果取消注释则每只蚂蚁不能再经过走过的点
% for j=1:size(dots{k},2) %不过同一点
% if dots{k}(1,end)+direc(1,di)==dots{k}(1,j) && dots{k}(2,end)+direc(2,di)==dots{k}(2,j)
% if ~isempty(permission)
% permission(:,end)=[];
% end
% end
% end
end
ii=size(permission,2);
if ii>=1
a=1:ii;%许可方向表的索引
%{
update log:
通过函数更新概率更新规则，应该提供已知当前点位下周围所有可行点的更新概率
索引顺序和许可方向表相同
%}
p=prob(num,permission,dots,k,tao,eta,alpha,beta);
s=randsrc(1,1,[a;p]);%概率生成
dots{k}(:,end+1)=dots{k}(:,end)+permission(:,s);
else
break
end
end
end
%计算路径代价
[W,cost_road]=cost_track(dots,thr_space,target,radius);
%优选路径对比历史最优
if num<=1
minimum=min(cost_road);
which=find(cost_road==minimum);
rec_best=dots{which};
hist_best=rec_best;%XXXXXXXX
hist_best_score(end+1)=minimum;
end
[rec_best,hist_best,hist_best_score,hist_best_step]=compete(cost_road,dots,hist_best,hist_best_score,thr_space,target,radius);
%信息素和启发信息更新
%每次迭代计算更新一下信息素
tao=TAO(tao,dots,hist_best,hist_best_step,cost_road,W,rou,e,Q,x_max,y_max);
%还有启发信息
eta=0.01*ones(x_max+1,y_max+1);
for xmap=1:x_max+1
for ymap=1:y_max+1
for k=1:length(dots)
if cost_road(k)~=inf
for j1=1:size( dots{k} ,2)-1
if xmap==dots{k}(1,j1) && ymap==dots{k}(2,j1)
%fprintf('%d %d %d %d\n',xmap,ymap,k,j1)
eta(xmap,ymap)=eta(xmap,ymap)+1./(W{k}(j1)+0.01);
break %经过好几次也不给加太多信息素
end
end
end
end
end
end
% 显示每次迭代的最佳路线（如何做？）
disp(num)
end
%------------------迭代计算结束-----------------------
% 结果可视化
for i =2:length(hist_best)
plot(hist_best(1,i-1:i),hist_best(2,i-1:i),'g')
axis([0 x_max 0 y_max])
pause(0.03)
end
hold off
pause(0.5);
%路径代价显示
figure(2)
i =1:length(hist_best_score);
plot(i,hist_best_score);
% 说明
% 依照概率选取随机数
% a=[1 2 3 4 5]
% Prob=[0.31 0.07 0.33 0.09 0.2]
% S=randsrc(1,1,[a;Prob])
以下为用到的定义的函数
1、伪随机规则
%点概率也可以空间分布prob(x,y)
function[p]=prob(~,permission,dots,k,tao,eta,alpha,beta)
%初次迭代只有信息素没有启发信息
%求当前点位周围所有可行点的更新概率
%索引顺序和许可方向表相同
pawn=[];
begin=dots{k}(:,end);
for i3=1:size(permission,2)
x=begin(1)+permission(1,i3);
y=begin(2)+permission(2,i3);
% fprintf('%d %d %d\n',x,y,i3)
pawn(i3)=tao(x+1,y+1).^alpha.*eta(x+1,y+1).^beta;
end
% 防止被困
if max(pawn)==0
pawn=ones(1,length(pawn));
end
p=pawn/(sum(pawn));
end
2、代价函数
% 代价函数
% 对于空间中每个走过的点，代价函数需要提供需求W{k}(j1)
% 自变量分别是上一轮蚂蚁编号和路径段编号
%{
——————————————————————
update log:
增加了输出可行段代价和整条路径代价的函数
W 以每条边的终点近似计算威胁代价，并结合边长，将代价更新在边的终点位置上
W和dots为相同规模和顺序的元胞数组，因此通过和dots相同的方式索引
——————————————————————
cost_road记录每一条路径的综合代价，顺序和dots的元胞索引相同
——————————————————————
%}
function [W,cost_road]=cost_track(dots,thr_space,target,radius)
k=0.95;% 威胁代价权重，另一边是油耗
W={};
cost_road=[];
for i1 =1:length(dots)% 此时是蚂蚁的总数
W{i1}=[];
cost_road(i1)=inf;
if norm(dots{i1}(:,end)-target)<=radius % 达到目的地的才能算每条痕迹
for j1=2:size(dots{i1},2) % 每只蚂蚁路径的段数
xx=dots{i1}(1,j1);
yy=dots{i1}(2,j1);
% fprintf('%d ',xx);
% fprintf('%d ',yy);
% fprintf('%d ',i1);
% fprintf('%d ',j1);
% fprintf('\n');
wf=norm( dots{i1}(:,j1)-dots{i1}(:,j1-1) );
wti=norm( dots{i1}(:,j1)-dots{i1}(:,j1-1) )*thr_space(xx+1,yy+1);%路径终点的威胁估算整条路径
W{i1}(j1-1)=k*wti+(1-k)*wf;
end
cost_road(i1)=sum( W{i1} );
else
end
end
end
3、启发信息更新（记录历史最优）
function [rec_best,hist_best,hist_best_score,hist_best_step]=compete(cost_road,dots,hist_best,hist_best_score,thr_space,target,radius)
% Judge & recorder
% rec_best记录当前最优点列、hist_best记录历史最优点列,hist_best_score是记录点列得分进化过程
minimum=min(cost_road);
which=find(cost_road==minimum);
rec_best=dots{which};
if minimum<=hist_best_score(end)
hist_best=dots{which};
hist_best_score(end+1)=minimum;
end
hist_best_step=cost_track({hist_best},thr_space,target,radius);
end
4、信息素更新
% 信息素更新，直接更新信息素空间，输出的tao通过tao(x,y)索引
function [tao]=TAO(tao,dots,hist_best,hist_best_step,cost_road,W,rou,e,Q,x_max,y_max)
%*****没有考虑到扩散******
%***在当前代码结构下扩散可以在本函数上添加***
%现有的信息素在扩散
%上回走通的路沿途各点都会留下信息素
%历次走过的最好的路各点依旧强化信息素
for xmap=1:x_max+1
for ymap=1:y_max+1
tao(xmap,ymap)=(1-rou)*tao(xmap,ymap);
% 最近有效航路
for k=1:length(dots)
if cost_road(k)~=inf
for j1=1:size( dots{k} ,2)-1
if xmap==dots{k}(1,j1) && ymap==dots{k}(2,j1)
%fprintf('%d %d %d %d\n',xmap,ymap,k,j1)
tao(xmap,ymap)=tao(xmap,ymap)+rou*Q./(W{k}(j1)+0.01);
break %经过好几次也不给加太多信息素
end
end
end
end
% 历史最佳航路
for j2=1:size(hist_best,2)-1
if xmap==hist_best(1,j2) && ymap==hist_best(2,j2)
%fprintf('%d %d %d\n',xmap,ymap,j2)
if size(hist_best_step{1},2)>=j2
tao(xmap,ymap)=tao(xmap,ymap)+rou*e*Q./(hist_best_step{1}(j2)+0.01);
break %经过好几次也不给加太多信息素
end
end
end
end
end
tao1=tao;
for xmap=2:x_max
for ymap=2:y_max
for i=-1:1
for j=-1:1
if~ i==0&&j==0
tao(xmap,ymap)=0.6*tao1(xmap,ymap)+0.4/8*(tao1(xmap+i,ymap+j));
end
end
end
end
end
end
运行结果
图：左侧为实现的路径生成，右侧为最优代价——更新次数曲线。
写在最后
本人纯小白，本科期间没有参加过什么涉及算法的竞赛，目前对于这一方面也还处在摸索阶段。第一次写博客，对算法的理解还不太深入，有很多疏漏的地方，还请大家多多包涵，此外如有改进建议还请多多指教。
更多推荐
7
0- 0
已为社区贡献1条内容
所有评论(0)