LEACH协议matlab程序
clear;%清除內存变量
%%
t1 = clock;
xm=100;%x轴范围
ym=100;%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=1000;%防止簇头过早死亡，需定期更换簇头，即为总轮数
CM=32;%控制信息大小
DM=4000;%要传输的信息大小，一个数据包
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;%轮数
%如果轮数正好是一个周期的整数倍，则设置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; %每轮死亡节点数
% 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;%初始化簇头包含的节点数，默认有一个节点，故后面注意减一
end
end
end
end
% CH_Num(r+1)=cluster; %每轮的簇头数
%%
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(min_dis_cluster)=packet_To_BS(min_dis_cluster)+1;%将此节点加入的簇凑节点，簇所包含的节点数1，每个节点包含一个数据包
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数组中
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
hold on;
end
%%
for k = 1:1:rmax
figure(2)
if(k == first_dead)
plot(k,Dead(k),’*y’),
else
plot(k,Dead(k),’-.r’),
end
axis([0,rmax,0,100]);
hold on;
end
