蚁群算法+LEACH 点滴（六） WSN簇头的路径规划（蚁群算法+LEACH）

2019年3月22日 191次阅读来源: 蚁群算法

今天BOSS来了，让他看了下实验，评价了下，说是思路不是很好，单纯的蚁群算法单在簇头路径上的应用效果肯定不是很好。想了想也是这样，虽然有点失望，但是只能就受了。决定还是将实验内容贴出来，这样对的起自己的努力吧~

实验效果只是在死亡节点的分布上取得了一些成果，死亡的节点分布比较开，也算上对自己的努力的一种肯定吧。

代码：

ANTLEACH1.m

clc;
clear;
%%
t1 = clock;
xm=400;%x轴范围
ym=400;%y轴范围
x = 0:180;%时间轴

sink.x=0.5*xm;%基站x轴
sink.y=0.5*ym;%基站y轴

n=100;%总的节点数
p=0.1;%总的簇头比例
E0=0.02;%初始能量

ETX=50*0.000000000001;%传输消耗能量，每bit
ERX=50*0.000000000001;%接收消耗能量，每bit
Efs=10*0.000000000001;%耗散能量，每bit
EDA=5*0.000000000001;%融合能耗，每bit

cc=0.6;%融合率
rmax=300;%防止簇头过早死亡，需定期更换簇头，即为总轮数
CM=32;%控制信息大小，包含自己簇内需要转发的总的数据包个数
DM=4000;%要传输的信息大小，一个数据包Bit数
figure(1);%显示图片

%%
for i=1:1:n
%随机产生节点
S(i).xd=rand(1,1)*xm;
S(i).yd=rand(1,1)*ym;

S(i).G=0;%若为0，标示当前节点可为簇头，每一周期全部置为0，重新指定全部簇头
S(i).E=E0;%设置初始化能量为E0
S(i).type=’N’;%节点类型为普通型

plot(S(i).xd,S(i).yd,’o’);
hold on;%保持所画的图像
end%为每个节点随机分配坐标，并设置初始化能量为E0，节点类型为普通型
%%
S(n+1).xd=sink.x;
S(n+1).yd=sink.y;
plot(S(n+1).xd,S(n+1).yd,’x’);%绘制基站节点
flag_first_dead=0;%第一个节点死亡的标志变量
syms count;
count = 0;

for r=1:1:rmax%开始每个循环
%%
r+1;%轮数
str1 = sprintf(‘第********** %d **********轮.’,r);
disp(str1);
%如果轮数正好是一个周期的整数倍，则设置S(i).E为0—此时所有节点均参见竞争
if(mod(r,round(1/p))==0)%round一个四舍五入函数
for i=1:1:n
S(i).G=0;
end
count = count + 1;
end
hold off;%每轮照片重新绘制
cluster=0;%初始化簇头为0
dead=0;%初始死亡节点数0
figure(1);
%重新绘制整体节点
%%
for i=1:1:n
if(S(i).E<=0)
plot(S(i).xd,S(i).yd,’red .’);
dead=dead+1;%节点死亡后变为红色，死亡节点数加一
if(dead==1&&flag_first_dead==0)
first_dead=r;%第一个节点的死亡轮数
% save ltest, first_dead;%将first_dead保存在ltest.mat中
flag_first_dead=1;
disp([‘第一个能量耗尽节点出现时间为第’,num2str(etime(clock,t1)),’秒’]);
pause
end
hold on;%新画图像之后不想覆盖原图像就要加上hold on
else
S(i).type=’N’;
plot(S(i).xd,S(i).yd,’o’);%绘制其他节点
hold on;
end
end
plot(S(n+1).xd,S(n+1).yd,’x’);%绘制基站
Dead(r+1)=dead; %每轮死亡节点数

%%
%获取簇头节点数组
for i=1:1:n
if(S(i).E>0)
if(S(i).G<=0)%不是簇头
temp_rand=rand;%取一个随机数
if(temp_rand<=(p/(1-p*mod(r,round(1/p)))))%如果随机数小于等于
S(i).type=’C’;%此节点为此轮簇头节点
S(i).G=round(1/p)-1;%S(i).G设置为大于0，此周期不能再被选为簇头节点
cluster=cluster+1;%簇头数加1
C(cluster).xd=S(i).xd;
C(cluster).yd=S(i).yd;%将此节点标志为簇头
plot(S(i).xd,S(i).yd,’k*’);%绘制此簇头，黑色

distance=sqrt((S(i).xd-(S(n+1).xd))^2+(S(i).yd-(S(n+1).yd))^2);%簇头到基站的距离
C(cluster).distance=distance;%标志位此簇头到基站的距离
C(cluster).id=i; %此簇头节点的id

%范围分类，最远为3的范围内
if(C(cluster).distance >=0 && C(cluster).distance <= xm/6)
C(cluster).group = 1.0;
elseif(C(cluster).distance >=xm/6 && C(cluster).distance <= (xm/6)*2)
C(cluster).group = 2.0;
else
C(cluster).group = 3.0;
end
packet_To_C(cluster)=1;%初始化簇头包含的节点数，默认有一个节点，故后面注意减一
end
end
end
end
str = sprintf(‘初始化生成 %d 个簇头.’,cluster);
disp(str);
% CH_Num(r+1)=cluster; %每轮的簇头数
disp(‘#################################################’);
[Route,D,NeedNodes,TiaoNum] = ANTLEACH(C,cluster,sink.x,sink.y,10);%NeedNodes范围3中节点数


%%
%普通节点的处理
for i=1:1:n

if(S(i).type==’N’&&S(i).E>0)%对每个能量大于0且非簇头节点

min_dis=sqrt((S(i).xd-(C(1).xd))^2+(S(i).yd-(C(1).yd))^2);%计算此节点到簇头1的距离
min_dis_cluster=1;%初始设置到簇头1的距离最近

for c=2:1:cluster
temp=sqrt((S(i).xd-(C(c).xd))^2+(S(i).yd-(C(c).yd))^2);
if(temp<min_dis)
min_dis=temp;
min_dis_cluster=c;
end
end%实现了选择此节点到哪个簇头节点的距离最短
%packet_To_BS为自己簇内需要转发的包数目
packet_To_C(min_dis_cluster)=packet_To_C(min_dis_cluster)+1;%该最近距离的簇头的非簇头节点数，一个非簇头节点包含一个数据包
Er1=ERX*CM*(cluster+1);%此节点接收各个簇头和BS（+1为BS的？）的控制信息
Et1=ETX*(CM+DM)+Efs*(CM+DM)*min_dis*min_dis;%此节点发送加入信息和发送数据信息到簇头的能耗，假定每一轮发送一次
S(i).E=S(i).E-Er1-Et1;%此轮的剩余能量

end

end



%%
%簇头节点的处理
if(TiaoNum == 3)%如果路径上除BS外存在两跳，即3个簇头
RouteDis=zeros(NeedNodes,3);%将路径最后一个节点到BS距离加上
for i=1:NeedNodes%所有蚂蚁的路径长度对应于L中
Dis=Route(i,:);
for j=1:2
RouteDis(i,j)=D(Dis(j),Dis(j+1)); %原距离加上第j个城市到第j+1个城市的距离
end
RouteDis(i,3)=D(Dis(end),cluster+1);

end

%%
%区域3与2中簇头，通过比较获取其到BS的最短距离

RouteRealtemp = cell(NeedNodes,1);%更新够最短路径表
RouteRealDistemp = cell(NeedNodes,1);%更新后最短距离表
RouteReal = cell(NeedNodes,1);%更新够最短路径表
% RouteRealDis = cell(NeedNodes,1);%更新后最短距离表
flag = -1;
shortdis = 10;
%
% Route

for i=1:NeedNodes%将BS添加到路径当中，通过条件判断添加的位置
%%%第一次裁剪，靠距离裁剪
%第一跳
if ((RouteDis(i,1)+shortdis) >= D(Route(i,1),cluster+1))
RouteRealtemp{i} = [Route(i,1) cluster+1];%cluster+1为BS编号
RouteRealDistemp{i} = D(Route(i,1),cluster+1);
flags = 1;
else
RouteRealtemp{i} = Route(i,:);
RouteRealDistemp{i} = RouteDis(i,:);
flags = -1;
end
%第二跳
if(flags == -1)%第二跳不为上面第一跳计算后的BS编号

if ((RouteDis(i,2)+shortdis) >= D(Route(i,2),cluster+1))
RouteRealtemp{i}(3) = (cluster+1);%距离不划算，将第三跳改为BS编号
RouteRealDistemp{i}(1) = RouteDis(i,1);
RouteRealDistemp{i}(2) = D(Route(i,2),cluster+1);

else
RouteRealtemp{i} = Route(i,:);
RouteRealtemp{i}(end + 1) = (cluster+1);%路径添加BS编号
RouteRealDistemp{i} = RouteDis(i,:);
end
end
% RouteRealtemp{i}
% RouteRealDistemp{i}

%%第二次裁剪,能量裁剪
if(RouteRealtemp{i}(2) ~= (cluster+1))
tempRoute = RouteRealtemp{i};
tempRouteDis = RouteRealDistemp{i};
templength = length(tempRouteDis);
Node_To_BS_Dis = D(tempRoute(1),cluster+1);%该节点到BS的距离
Node_Se_To_BS_E = 0;
Node_EX_To_BS_E = 0;
temp_Ex_Dis = 0;%转发的耗散距离
% packet = packet_To_C(i)%
% tR = tempRoute(1)
Node_Se_To_BS_E = Efs*DM*cc*( packet_To_C(tempRoute(1)) -1)*C(tempRoute(1)).distance*C(tempRoute(1)).distance;%该节点直接发送到BS的耗散能量
for j = 1:1:templength
temp_Ex_Dis = temp_Ex_Dis + RouteRealDistemp{i}(j);
end
%(templength – 1)为所经历的跳数
Node_EX_To_BS_E = (templength – 1)*DM*cc*(packet_To_C(tempRoute(1))-1)*(ETX+ERX) + Efs*DM*(packet_To_C(tempRoute(1))-1)*temp_Ex_Dis^2;%经转发的接收、转发与耗散的能量消耗
if(Node_EX_To_BS_E*9/10 >= Node_Se_To_BS_E)%转发能耗的2/3还大于直接发送，则没有转发的必要
RouteReal{i}(1) = tempRoute(1);
RouteReal{i}(2) = (cluster+1);
else
RouteReal{i} = tempRoute;
end
end
end
elseif(TiaoNum == 1)
RouteReal = cell(NeedNodes,1);%更新够最短路径表
% RouteRealDis = cell(NeedNodes,1);%更新后最短距离表
disp(‘TiaoNum == 1’);
for i=1:NeedNodes
RouteReal{i}(1) = Route(i,1);
RouteReal{i}(2) = (cluster+1);
% RouteRealDis{i} = D(Route(i,1),(cluster+1));
end
end

%%
%通过查找区域3中的簇头节点路径表，构成总的簇头路径路由表
RouteTotal = cell(cluster,1);%总的路由多跳表

for num = 1:cluster%生成所有的簇头路由多跳路由表
rowind = 0;
for i = 1:NeedNodes%区域3 簇头数

indtemp = find(RouteReal{i} == num);
if((indtemp ~= 0))
if(indtemp ~= length(RouteReal{i}))
rowind = i;
break;
end
end
end
if(rowind ~= 0)
columnind = min(indtemp);
RouteTotal{num} = RouteReal{rowind}(columnind:end);
else
RouteTotal{num} = [num,cluster+1];
end
end
%%
%生成簇头相互间直接通讯表，即跳表—最终需要的
DirectCon = zeros(cluster,2);
for i = 1:1:cluster%生成所有的簇头路由下一跳路由表
DirectCon(i,1) = RouteTotal{i}(1);
DirectCon(i,2) = RouteTotal{i}(2);
end



%%
%需要转发的数据包统计（不包含自己的）

Trans_Dis = zeros(cluster,1);%两个簇头节点之间的转发距离，即：该点到下一跳的距离，根据DirectCon得到
packet_To_BS_Transmit = zeros(cluster,1);%簇头需要转发其他簇头数据的总的数据包个数
for c = 1:1:cluster
% packet_To_BS_Transmit(c,1) = packet_To_BS(c);%初始化需要转发的数据量为本身的数据包
NodeNum = find(DirectCon(:,2) == c);
Num = length(NodeNum);
if(Num > 0)
for j = 1:1:Num
packet_To_BS_Transmit(c,1) = packet_To_BS_Transmit(c,1) + (packet_To_C(NodeNum(j))-1);
end
end
% disp(‘////////////////////////////////////////////’);
% cluster
% packet_To_BS_Transmit

Trans_Dis(c,1) = D(DirectCon(c,1),DirectCon(c,2));

end
% Trans_Dis
%%

for c=1:1:cluster%各个簇头，
% packet_To_BS(c);%簇头自己需发送到基站的数据包个数
% packet_To_BS_Transmit(c,1) = packet_To_BS(c);%初始化需要转发的数据量为本身的数据包
CEt1=ETX*CM+Efs*CM*(sqrt(xm*ym))*(sqrt(xm*ym));%此簇头广播成簇信息的能耗
CEr1=ERX*CM*(packet_To_C(c)-1+1);%簇头收到各个节点加入信息的能耗，初始化为一个节点，故要减一，加一为BS发送的跳表信息
CEr2=ERX*DM*(packet_To_C(c)-1);%簇内部接收数据包能耗

% CEt2=(ETX+EDA)*DM*cc*(packet_To_BS(c)-1)+Efs*DM*cc*(packet_To_BS(c)-1)*C(c).distance*C(c).distance;
CEbs = ERX*CM;%簇头接收BS发送的路由跳表的能耗
CEda = EDA*DM*cc*(packet_To_C(c)-1); %簇内部的数据融合能耗
CErc = ERX*DM*packet_To_BS_Transmit(c);%接收其它簇头的能耗
CEsc = ETX*DM*cc*(packet_To_C(c)-1) + ETX*packet_To_BS_Transmit(c);%总的数据包发送能耗（自己簇内融合+接收的）
CEfs = Efs*DM*(cc*(packet_To_C(c)-1) + packet_To_BS_Transmit(c))*(Trans_Dis(c))^2;%总的发送耗散能耗,自己+转发

S(C(c).id).E=S(C(c).id).E-CEr1-CEr2-CEt1-CEbs-CEda-CErc-CEsc-CEfs;%此轮后簇头的剩余能量
end
hold on;
end

%%
for k = 1:1:rmax
figure(2)
if(k == first_dead)
plot(k,Dead(k),’*b’),
else
plot(k,Dead(k),’-.r’),
end
axis([0,rmax,0,100]);
hold on;
end

ANTLEACH.m

function[Shortest_Route,D,cityS,RouteNodes] = ANTLEACH(G,n,mx,my,NC_max)%mx,my,基站坐标,n为节点规模,cityS表示范围3簇头数
%%
%%第一步：变量初始化
% 参考设置初始参数如下：
%NC_max=50
Alpha=1;Beta=5;Rho=0.1;Q=100;
CC = zeros(n,3);
for i = 1:n
CC(i,1) = G(i).xd;
CC(i,2) = G(i).yd;
CC(i,3) = G(i).group;
end

D=zeros(n,n+1);%D表示完全图的赋权邻接矩阵，n+1存储改点到基站的距离
m = fix(n/2);
cityS = 0;
cityNum = 0;
Eta = zeros(n,n+1);
%存放可访问城市
J=[]; %待访问的城市,定义为空矩阵
JS = [];
% Jc = 0;
digits(5);

for i=1:n
for j=1:n
if i~=j
if((CC(i,3) == 3&&CC(j,3) == 2)|| (CC(i,3) == 2&&CC(j,3) == 1))
D(i,j)=((CC(i,1)-CC(j,1))^2+(CC(i,2)-CC(j,2))^2)^0.5;
else
D(i,j)=100;
end
else
D(i,j)=100; %i=j时不计算，应该为0，但后面的启发因子要取倒数，用eps（浮点相对精度）表示

end
end
if(CC(i,3) == 3)

cityS = cityS+1;
JS(cityS)=i;%存储起始点
else

cityNum = cityNum+1;
J(cityNum)=i;%对可收索城市编号,也是在D中的编号

end
D(i,n+1) = ((CC(i,1)-mx)^2+(CC(i,2)-my)^2)^0.5;

end
Eta = 1./D;
D = vpa(D);
Eta = vpa(Eta);

%%
if (cityNum >= 3&&cityNum >= 0.5*cityS)%至少3个簇头节点可以参与转发，并且转发密度不太大

layers = 3;%路径长度
RouteNodes = layers;
Shortest_Route=zeros(cityS,layers); %存储并记录该节点最后的最佳路径
% Tabu = zeros(cityS,layers+1);
for i = 1:cityS
Shortest_Route(i,1) = JS(i);%城市标号
end

str = sprintf(‘共%d个远距离节点需要查找路径.’,cityS);
disp(str);

for st = 1:cityS%先寻找一个节点的路径
NC=1; %迭代计数器，记录迭代次数
Tau=ones(n,n); %Tau为信息素矩阵
L_best=inf.*ones(NC_max,1); %各代最佳路线的长度
R_best=zeros(NC_max,layers); %各迭代最佳路线
Delta_Tau=zeros(n,n); %开始时信息素为n*n的0矩阵
str = sprintf(‘第%d个节点路径查找中…’,st);
disp(str);
disp(‘Busy…’);

while NC<=NC_max %停止条件之一：达到最大迭代次数，停止

TabuTemp=zeros(m,layers); %存储并记录每只蚂蚁的所有路径上的城市编号
TabuTemp(:,1) =Shortest_Route(st,1);%从JS中编号st的城市开始访问￥
JJ = cell(m,1);
for i = 1:m
JJ{i} = J;%城市标号$$$$$$$$$;%JCA(i,:)
end

%%
%m只蚂蚁从第一个城市出发，发现最优路线，得到第一个城市最优路线
for j=2:layers %j表示当前收索到的最佳城市存储位置
for i=1:m%从第一个城市开始放m只蚂蚁
Jtemp = zeros(1,length(JJ{i}));
P = zeros(1,length(JJ{i}));
Jtemp = JJ{i};%第i只蚂蚁可以访问的集合

visited=TabuTemp(i,1:(j-1)); %记录已访问的城市，避免重复访问
P=Jtemp;

%下面计算待选城市的概率分布
for k=1:length(Jtemp)
P(k)=(Tau(visited(end),Jtemp(k))^Alpha)*(Eta(visited(end),Jtemp(k))^Beta);%end为矩阵末位索引
end

P=P/(sum(P));
PMax = max(P);
R = vpa(PMax);
rd = rand(1)*(0.001-R)+R;%0.001~R随机数,防止越界
select_index_list = find(P >=rd);%赌轮法
to_visit_city_index = Jtemp(select_index_list(1));
JJ{i}(select_index_list(1)) = [];
TabuTemp(i,j)=to_visit_city_index;%更新第i只蚂蚁的第j个城市

end
end

%%
%%记录本只蚂蚁迭代最佳路线
L=zeros(m,1); %开始距离为0，m*1的列向量
for i=1:m%所有蚂蚁的路径长度对应于L中
R=TabuTemp(i,:);
for j=1:(layers-1)
L(i)=L(i)+D(R(j),R(j+1)); %原距离加上第j个城市到第j+1个城市的距离
end
L(i)=L(i)+D(R(end),n+1); %路径上最后一个城市直接返回（这样合理？？）
end
L_best(NC)=min(L); %保存当前最佳距离
pos=find(L==L_best(NC));
R_best(NC,:)=TabuTemp(pos(1),:); %此轮迭代后的最佳路线

%%
%%第五步：更新信息素

for i=1:m
for j=1:(layers-1)
Delta_Tau(TabuTemp(i,j),TabuTemp(i,j+1))=Delta_Tau(TabuTemp(i,j),TabuTemp(i,j+1))+Q/L(i);
%此次循环在路径（i，j）上的信息素增量
end
Delta_Tau(TabuTemp(i,layers),TabuTemp(i,1))=Delta_Tau(TabuTemp(i,layers),TabuTemp(i,1))+Q/L(i);
%最后一次返回留下信息素，此次循环在整个路径上的信息素增量
end
Tau=(1-Rho).*Tau+Delta_Tau; %考虑信息素挥发，更新后的信息素
%%
%%第六步：禁忌表清零
TabuTemp = zeros(m,layers);
NC=NC+1; %迭代继续

end
%%
%m只蚂蚁路径最佳
Pos=find(L_best==min(L_best)); %找到最佳路径（非0为真）
Shortest_Route(st,:)=R_best(Pos(1),:); %最大迭代次数后最佳路径

end
else
Shortest_Route = zeros(cityS,1);
RouteNodes = 1;
for i = 1:1:cityS
Shortest_Route(i,1) = JS(i);
end
end

死亡节点分布：
《蚁群算法+LEACH 点滴（六） WSN簇头的路径规划（蚁群算法+LEACH）》

死亡节点数目：（红色为改进）

《蚁群算法+LEACH 点滴（六） WSN簇头的路径规划（蚁群算法+LEACH）》

    原文作者：蚁群算法
    原文地址: https://blog.csdn.net/hu1583141033/article/details/8891040
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。