遗传算法改进leach协议(附MATLAB代码)

LEACH（低能量自适应聚簇层次）协议是无线传感器网络中最早的、最具有影响力的分簇协议之一。它通过随机选举一些传感器节点作为簇头，其他节点加入其中一个簇头组成一个簇，然后簇头负责汇聚数据并将数据发送给基站，从而达到减少能耗、提高网络寿命的目的。然而，LEACH协议的簇头选举过程随机性较强，可能导致簇头的能耗过快，进而影响网络性能。

为解决这个问题，可以采用遗传算法优化LEACH协议。遗传算法是一种模拟自然选择和遗传机制的全局搜索优化算法，具有优良的全局搜索能力和并行性。遗传算法可以用来优化簇头的选举和数据传输过程，以达到降低能耗、延长网络寿命的效果。

首先，构建初始种群，每个个体表示一种可能的簇头选举和数据传输方案。然后，对每个个体进行评价，得到适应度值，适应度函数可以设计为簇头的能耗和数据传输距离的函数。接着，根据适应度，采用遗传算子（选择、交叉、变异）产生新的种群。重复上述过程，种群将逐步向最优解区域进化。

通过遗传算法优化后的LEACH协议，可以有效地平衡簇头的能耗，保证簇内节点的能量消耗均衡，从而提高整个网络的能量效率，延长网络寿命。同时，遗传算法的并行性使得它可以在大规模网络中高效地寻找优化解，具有很高的实用价值和广阔的应用前景。

代码如下:

clear all;close all;clc;%清除内存变量
%% 遗传算法-多跳路由改进leach协议主程序

%% 输入参数列表
% BorderLength------------正方形区域的边长，单位：km
% NodeAmount--------------网络节点的个数
% Alpha-------------------网络特征参数，Alpha越大，短边相对长边的比例越大
% Beta--------------------网络特征参数，Beta越大，边的密度越大
% PlotIf------------------是否画网络拓扑图，如果为1则画图，否则不画图
% EdgeCostDUB-------------链路费用的控制参数，1×2，存储链路费用的下界和上界
% EdgeBandWideDUB---------链路带宽的控制参数，1×2，存储下界和上界
% VertexCostDUB-----------节点费用的控制参数，1×2，存储节点费用的下界和上界
% VertexDelayDUB----------节点时延的控制参数，1×2，存储节点时延的下界和上界
% VertexDelayJitterDUB----节点时延抖动的控制参数，1×2，存储节点时延抖动的下界和上界
% VertexPacketLossDUB-----节点丢包率的控制参数，1×2，存储节点丢包率的下界和上界
%% 输出参数列表
% Sxy---------------------3×N的矩阵，各列分别用于存储节点的序号，横坐标，纵坐标的矩阵
% AM----------------------01存储的邻接矩阵，AM(i,j)=1表示存在由i到j的有向边，N×N
% EdgeCost----------------链路费用矩阵，N×N
% EdgeDelay---------------链路时延矩阵，N×N
% EdgeBandWide------------链路带宽矩阵，N×N
% VertexCost--------------节点费用向量，1×N
% VertexDelay-------------节点时延向量，1×N
% VertexDelayJitter-------节点时延抖动向量，1×N
% VertexPacketLoss--------节点丢包率向量，1×N
% DisTabel----------------节点距离矩阵，N×N

BorderLength=1000; %正方形区域的边长，单位：km
NodeAmount=100; %网络节点的个数
Alpha=100000000; %网络特征参数，Alpha越大，短边相对长边的比例越大
Beta=200000000000; %网络特征参数，Beta越大，边的密度越大
PlotIf=0; %是否画网络拓扑图，如果为1则画图，否则不画图
EdgeCostDUB=[5,5]; %链路费用的下界和上界
EdgeBandWideDUB=[30,1000]; %链路带宽的下界和上界
VertexCostDUB=[3,3]; %节点费用的下界和上界
VertexDelayDUB=1e-4*[5,20]; %节点时延的下界和上界
VertexDelayJitterDUB=1e-4*[3,8]; %节点时延抖动的下界和上界
VertexPacketLossDUB=1e-4*[0,500]; %节点丢包率的下界和上界
% figure;
[Sxy,AM,EdgeCost,EdgeDelay,EdgeBandWide,VertexCost,VertexDelay,VertexDelayJitter,VertexPacketLoss,DisTabel]=...
NetCreateNew2(BorderLength,NodeAmount,Alpha,Beta,PlotIf,EdgeCostDUB,EdgeBandWideDUB,VertexCostDUB,VertexDelayDUB,VertexDelayJitterDUB,VertexPacketLossDUB);
% title('网络拓扑（满足带宽）');
xm=BorderLength;%x轴范围
ym=BorderLength;%y轴范围

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

n=NodeAmount;%节点总数

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=8;%总轮数

CM=32;%控制信息大小
DM=4000;%数据信息大小

figure(1);%显示图片

%% 设置数据
for i=1:1:n
S(i).xd=Sxy(2,i);
S(i).yd=Sxy(3,i);
S(i).G=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;%第一个死亡节点的标志变量

for r=1:1:rmax%开始每轮循环

if(mod(r,round(1/p))==0)
for i=1:1:n
S(i).G=0;
end
end%如何轮数正好是一个周期的整数倍，则设置S(i).G为0

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;%将能量小于等于0的节点绘制成红色，并将死亡节点数增加1

if(dead==1)
if(flag_first_dead==0)
first_dead=r; %第一个节点的死亡轮数
save ltest, first_dead;
flag_first_dead=1;
end
end%将能量小于等于0的节点绘制成红色，并将死亡节点数增加1

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; %每轮有死亡节点数
save ltest, Dead(r+1);%将此数据存入ltest文件
%% ---------------------选举簇头部分开始---------------
%不考虑簇头处理能力
%（1）计算簇头数cluster
nodeSET=[];
for i=1:1:n%对每个节点循环
if(S(i).E>0)%判断节点能量
nodeSET=[nodeSET i];%加入可用节点集
if(S(i).G<=0)
temp_rand=rand;%取一个随机数
if(temp_rand<=(p/(1-p*mod(r,round(1/p)))))%如果随机数小于等于
cluster=cluster+1;%簇头数加1
end
end
end
end%for i=1:1:n
%% 遗传算法
% （2）先找可用于选举的节点
popsize=20; %个体数目
MAXGEN=10;%最大遗传代数
GGAP=1;%代沟
crossp=0.9;%交叉率
mutionp=0.9;%变异率
trace1=[];
Chrom=genchrom01(nodeSET,cluster,popsize);%产生初始种群
gen=0;
bestchrom=cell(MAXGEN,1);
while gen<=MAXGEN
Value=DecodeLeachChrom(Chrom,S,nodeSET,cluster); %计算目标函数值
FitnV1=ranking(Value);%分配适应度值
Chrom=select('sus',Chrom,FitnV1,GGAP);%选择
Chrom=mutationLeach(Chrom,cluster,mutionp);%变异
Chrom=crossLeach(Chrom,crossp);%交叉
Value=DecodeLeachChrom(Chrom,S,nodeSET,cluster); %计算目标函数值
gen=gen+1;
[trace1(gen,1),index1]=min(Value);%记录最小值
bestchrom{gen,1}=Chrom(index1,:);
trace1(gen,2)=sum(Value)/popsize;%记录均值
end
temp01=bestchrom{gen,1};%最佳分簇方案
cluster_set=cell(cluster,1);
cluster_node=zeros(cluster,1);
packet_To_BS=zeros(cluster,1);
for k0001=1:cluster
[v01,index2]=find(temp01==k0001);
cluster_set{k0001,1}=nodeSET(index2);%分簇
%% ------------簇头选举开始：不考虑周期----------------------
temp02=cluster_set{k0001,1};
if ~isempty(temp02)
% 计算簇平均能量
E01=zeros(1,length(temp02));
for k0002=1:length(temp02)
E01(1,k0002)=S(temp02(k0002)).E;
end
distance01=calfitnessfun2(temp02,S);%到基站距离
%最终目标
OE1=1./E01.^2;
% fitness01=distance01.*OE1;
fitness01=distance01;
[v03,index3]=min(fitness01);
selectnode=temp02(index3);
cluster_node(k0001,1)=selectnode;
S(selectnode).type='C';%此节点为此轮簇头
S(selectnode).G=round(1/p)-1;%S(i).G设置为大于0，此周期不能再被选择为簇头
C(k0001).xd=S(selectnode).xd;
C(k0001).yd=S(selectnode).yd;%将此节点标志为簇头
plot(S(selectnode).xd,S(selectnode).yd,'k*');%绘制此簇头
distance=sqrt((S(selectnode).xd-(S(n+1).xd))^2+(S(selectnode).yd-(S(n+1).yd))^2);%簇头到基站的距离
C(k0001).distance=distance;%标志为此簇头的距离
C(k0001).id=selectnode; %此簇头的节点id
packet_To_BS(k0001)=1;%发送到基站的数据包数为1
end
%% ------------簇选举结束----------------------
end%for k0001=1:cluster
%% ---------------------选举簇头部分结束---------------
Delay01=0;%时延
DelayJitter01=0;%时延抖动
PacketLoss01=0;%丢包率
for c=1:1:cluster%各个簇头
cluster_set001=cluster_set{c,1};%节点集
[Ttable02,AM1]=topo_table23(AM,cluster_set001);
if ~isempty(cluster_set001)
number_dem=C(c).id;%得到簇头作为源节点
E=cluster_set001(cluster_set001~=number_dem);%目的节点
ChromRoute=zeros(length(E),length(cluster_set001));%路径记录
for w01=1:length(E)
number_source=E(w01);
%计算距离
distance_s2d=DisTabel(number_source,number_dem);%距离初始化
route=[number_source];%路径初始化
Ttable01=Ttable02;%临时拓扑赋值，防止搞坏原始拓扑集
%% GPRS主函数
state001=0;
while distance_s2d~=0
%% 判断死链
if isempty(route)
disp('无法产生路由路径');
break;
end
curr_node=route(end);%路径当前最后一个节点
curr_topolopy=Ttable01(curr_node).topolopy;%获取当前拓扑相邻集合
%% -------------判断终止开始---------------
if ~isempty(curr_topolopy)%拓扑相邻集合非空就找下
curr_distance1=DisTabel(curr_topolopy,number_dem);%获取当前拓扑相邻集合与目的地节点距离集合
distance01=DisTabel(number_dem,curr_node);%当前节点和目的节点距离
setindex01=curr_topolopy(curr_distance1<distance01);%比当前节点更靠近目的的邻近节点集
if isempty(setindex01)%就当前节点最近，那么直接结束
route=[route number_dem];
break;
else%如果不是空的，要继续判断
curr_distance3=DisTabel(setindex01,curr_node);%获取当前拓扑相邻集合与邻近节点距离集合
setindex02=setindex01(curr_distance3<distance01);
if isempty(setindex02)%即使距离目的近，但两者距离远过distance01，那么直接结束
route=[route number_dem];
break;
else%

end
end
end
%% -------------判断终止结束---------------
%% ------------------删除curr_topolopy中属于route的节点开始-----
n001=length(route);%计算route长度
for k=1:n001
v01=find(curr_topolopy==route(k));%route(k)是不是curr_topolopy的成员
if numel(v01)>0%是成员
curr_topolopy(v01)=[];%删除该元素
end
end
%% ------------------删除curr_topolopy中属于route的节点结束-----
if ~isempty(curr_topolopy)%拓扑相邻集合非空就找下
curr_distance1=DisTabel(curr_topolopy,number_dem);%获取当前拓扑相邻集合与目的地节点距离集合
distance01=DisTabel(number_dem,curr_node);%当前节点和目的节点距离
setindex01=curr_topolopy(curr_distance1<distance01);%比当前节点更靠近目的的邻近节点集
if isempty(setindex01)%就当前节点最近，那么直接结束
route=[route number_dem];
break;
else%如果不是空的，要继续判断
curr_distance3=DisTabel(setindex01,curr_node);%获取当前拓扑相邻集合与邻近节点距离集合
setindex02=setindex01(curr_distance3<distance01);
if isempty(setindex02)%即使距离目的近，但两者距离远过distance01，那么直接结束
route=[route number_dem];
break;
else%
[vt01,index01]=min(curr_distance3(curr_distance3<distance01));
route=[route setindex02(index01)];
end
end
end%if ~isempty(curr_topolopy)%拓扑相邻集合非空就找下
end%while distance_s2d~=0
ChromRoute(w01,1:length(route))=route;%输出路径
end%for w01=1:length(E)
% 画簇路由树
drawgoodroute23(Sxy,ChromRoute);

end%if ~isempty(cluster_set001)
% ChromRoute=zeros(length(E),length(cluster_set001));%路径记录
cset01=[];
als01=zeros(1,length(E));
for w009=1:length(E)
temp009=ChromRoute(w009,:);
temp009=temp009(temp009>0);
cset01=[cset01 temp009(end-1)];
als01(1,w009)=length(temp009);
end
[v02,index02]=sort(als01);
index02=fliplr(index02);
index03=index02(1:length(unique(cset01)));%前几个最长的
for w009=1:length(index03)
routenew01=ChromRoute(index03(w009),:);
routenew01=routenew01(routenew01>0);
Delay01=Delay01+sum(VertexDelay(routenew01));%节点时延
DelayJitter01=DelayJitter01+sum(VertexDelayJitter(routenew01));%节点时延抖动
PacketLoss01=PacketLoss01+sum(VertexPacketLoss(routenew01));%丢包率
for w010=1:length(routenew01)-1
% EdgeDelay,EdgeBandWide,VertexCost,VertexDelay,VertexDelayJitter,VertexPacketLoss
Delay01=Delay01+EdgeDelay(routenew01(w010),routenew01(w010+1));%路径时延
min_dis=norm([S(routenew01(w010)).xd,S(routenew01(w010)).yd]-[S(routenew01(w010+1)).xd,S(routenew01(w010+1)).yd]);
Er1=ERX*CM*(0+1);%此节点接收各个簇头的控制信息
%此节点加入的簇的簇头时隙控制信息的总接收能耗
Et1=ETX*(CM+DM)+Efs*(CM+DM)*min_dis*min_dis;%此节点发送加入信息和发送数据信息
%到簇头的能耗
S(routenew01(w010)).E=S(routenew01(w010)).E-Er1-Et1;%此轮后的剩余能量
end
end



packet_To_BS(c)=length(unique(cset01));%簇头需发送到基站的数据包个数
if packet_To_BS(c)~=0
CEr1=ERX*CM*(packet_To_BS(c)-1);%收到此簇各个节点加入信息的能耗
CEr2=ERX*DM*(packet_To_BS(c)-1);%收到此簇各个节点数据信息的能耗
CEt1=ETX*CM+Efs*CM*(sqrt(xm*ym))*(sqrt(xm*ym));%此簇头广播成簇信息的能耗
CEt2=(ETX+EDA)*DM*cc*packet_To_BS(c)+Efs*DM*cc*packet_To_BS(c)*C(c).distance*C(c).distance;%簇头将所以数据融合后发往基站的能耗
S(C(c).id).E=S(C(c).id).E-CEr1-CEr2-CEt1-CEt2;%此轮后簇头的剩余能量
end

end
CH_Num(r+1)=cluster; %每轮的簇头数
save ltest,CH_Num(r+1);%保存每轮簇头数到ltest
%保存数据
Delay02(r)=Delay01;%第r轮节点时延
DelayJitter02(r)=DelayJitter01;%第r轮时延抖动
PacketLoss02(r)=PacketLoss01;%第r轮丢包率

count002=0;
for i=1:1:n
R(r+1,i)=S(i).E; %每轮每节点的剩余能量
if S(i).E<=0
count002=count002+1;
end
end
if count002==n
break;
end
hold on;

end
title('遗传算法-多跳路由改进leach协议变化');
figure;
plot(Dead);
title('遗传算法-多跳路由改进leach协议死亡节点数');
figure;
plot(Delay02);
title('遗传算法-多跳路由改进leach协议每轮总时延');
figure;
plot(DelayJitter02);
title('遗传算法-多跳路由改进leach协议每轮总时延抖动');
figure;
plot(PacketLoss02);
title('遗传算法-多跳路由改进leach协议每轮总丢包率');

clear all;

[Sxy,AM,EdgeCost,EdgeDelay,EdgeBandWide,VertexCost,VertexDelay,VertexDelayJitter,VertexPacketLoss]=...
NetCreateNew2(BorderLength,NodeAmount,Alpha,Beta,PlotIf,EdgeCostDUB,EdgeBandWideDUB,VertexCostDUB,VertexDelayDUB,VertexDelayJitterDUB,VertexPacketLossDUB);

xm=BorderLength;%x轴范围
ym=BorderLength;%y轴范围

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

n=NodeAmount;%节点总数

p=0.1;%簇头概率

cc=0.6;%融合率

rmax=8;%总轮数

CM=32;%控制信息大小
DM=4000;%数据信息大小

figure(1);%显示图片

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;%第一个死亡节点的标志变量

for r=1:1:rmax%开始每轮循环
r+1;%显示轮数
if(mod(r,round(1/p))==0)
for i=1:1:n
S(i).G=0;
end
end%如何轮数正好是一个周期的整数倍，则设置S(i).G为0

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;%将能量小于等于0的节点绘制成红色，并将死亡节点数增加1

if(dead==1)
if(flag_first_dead==0)
first_dead=r %第一个节点的死亡轮数
save ltest, first_dead;
flag_first_dead=1;
end
end%将能量小于等于0的节点绘制成红色，并将死亡节点数增加1

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; %每轮有死亡节点数
save ltest, Dead(r+1);%将此数据存入ltest文件

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

packet_To_BS(cluster)=1;%发送到基站的数据包数为1
end
end
end
end

CH_Num(r+1)=cluster; %每轮的簇头数
save ltest,CH_Num(r+1);%保存每轮簇头数到ltest

Delay01=0;%时延
DelayJitter01=0;%时延抖动
PacketLoss01=0;%丢包率
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;
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(min_dis_cluster)=packet_To_BS(min_dis_cluster)+1;%将此节点加入的簇
%头节点数据包数加1
% EdgeDelay,EdgeBandWide,VertexCost,VertexDelay,VertexDelayJitter,VertexPacketLoss
Delay01=Delay01+VertexDelay(i)+VertexDelay(C(min_dis_cluster).id)+EdgeDelay(i,C(min_dis_cluster).id);%时延
DelayJitter01=DelayJitter01+VertexDelayJitter(i)+VertexDelayJitter(C(min_dis_cluster).id);%时延抖动
PacketLoss01=PacketLoss01+VertexPacketLoss(i)+VertexPacketLoss(C(min_dis_cluster).id);%丢包率
Er1=ERX*CM*(cluster+1);%此节点接收各个簇头的控制信息
%此节点加入的簇的簇头时隙控制信息的总接收能耗
Et1=ETX*(CM+DM)+Efs*(CM+DM)*min_dis*min_dis;%此节点发送加入信息和发送数据信息
%到簇头的能耗
S(i).E=S(i).E-Er1-Et1;%此轮后的剩余能量
end
end
for c=1:1:cluster%各个簇头
packet_To_BS(c);%簇头需发送到基站的数据包个数
CEr1=ERX*CM*(packet_To_BS(c)-1);%收到此簇各个节点加入信息的能耗
CEr2=ERX*DM*(packet_To_BS(c)-1);%收到此簇各个节点数据信息的能耗
CEt1=ETX*CM+Efs*CM*(sqrt(xm*ym))*(sqrt(xm*ym));%此簇头广播成簇信息的能耗
CEt2=(ETX+EDA)*DM*cc*packet_To_BS(c)+Efs*DM*cc*packet_To_BS(c)*C(c).distance*C(c).distance;%簇头将所以数据融合后发往基站的能耗
S(C(c).id).E=S(C(c).id).E-CEr1-CEr2-CEt1-CEt2;%此轮后簇头的剩余能量
end


%保存数据
Delay02(r)=Delay01;%第r轮节点时延
DelayJitter02(r)=DelayJitter01;%第r轮时延抖动
PacketLoss02(r)=PacketLoss01;%第r轮丢包率

for i=1:1:n
R(r+1,i)=S(i).E; %每轮每节点的剩余能量
% save ltest,R(r+1,i);%保存此数据到ltest
end

hold on;

end
title('基础leach协议');
figure;
plot(Dead);
title('基础leach协议死亡节点数');
figure;
plot(Delay02);
title('基础leach协议每轮总时延');
figure;
plot(DelayJitter02);
title('基础leach协议每轮总时延抖动');
figure;
plot(PacketLoss02);
title('基础leach协议每轮总丢包率');