文章目录

LEACH算法仿真实验
- 算法原理简介
- 代码分析
- - 初始变量声明
  - 初始化所有节点
  - LEACH分簇（簇的建立阶段）
  - 普通节点加入生成的簇头
  - 稳定的数据通信阶段
  - 计算死亡节点
- MATLAB画图
完整代码

LEACH算法仿真实验

本文使用MATLAB实现，因为MATLAB画图方便。

如果不想看思路分析，点这里→查看完整代码。

算法原理简介

基本思想是以循环的方式随机选择簇头节点，将整个网络的能量负载均匀分配到网络中的每个传感器节点，从而达到降低网络能耗，提高网络生存周期的目的。

详细来说，LEACH算法是一种自应分簇拓扑算法，它的执行过程是周期性的，每轮循环分为簇的建立阶段和稳定的数据通信阶段。

在簇的建立阶段，相邻节点动态地形成簇，随机产生簇头；
在数据通信阶段，簇内节点把数据发送给簇头，簇头进行数据融合并把结果发送给汇聚节点。

由于簇头需要完成数据融合、与汇聚节点通信等工作，所以能量消耗大。LEACH 算法能够保证各节点等概率地担任簇头，使得网络中的节点相对均衡地消耗能量。

代码分析

初始变量声明可以结合原理一起看，写在前面比较方便后面介绍。
这里→查看LEACH原理

初始变量声明

%-- 节点、簇头设置相关量 --%

本实验一共进行times次实验，每一次实验进行的轮数为 r轮，一共有n个节点，高级节点有m个。

  本次实验中：times = 10; % 实验次数n = 100; % 节点个数m=0.1; % 控制高级节点的占比r = 500; % 总轮数

%-- 控制坐标相关量 --%

每一个节点保存在一个结构体数组S中。S的属性有：

  1. xd : 节点横坐标2. yd : 节点纵坐标3. G ： 4. type ： 'C'为该节点当选过簇头节点，'N'为该节点没有当选过簇头节点5. E : 节点当前的能量值6. ENERGY : 1为高级节点，0为普通节点7. min_dis : 与 加入的簇头节点 的距离8. min_dis_cluster : 加入的簇头节点的编号S定义如下：S = repmat(struct('xd', 0, 'yd', 0, 'G', 0', 'type', 'N', 'E', 0, 'ENERGY', 0, 'min_dis', 0, 'min_dis_cluster', 0), n+1, 1); % 记录所有节点的坐标，及其相关属性

设置测试范围，所有的节点都在100x100的坐标系内，基站在坐标为(50, 50)的中心。
```
  xm=100;%x轴范围ym=100;%y轴范围sink.x=0.5*xm; %基站x轴sink.y=0.5*ym; %基站y轴
```

%-- 计算能量相关量 --%

高级节点的能量为普通节点能量的(1+a)倍。
```
  a=1;% 初始额外增加的能量占比
```
当节点能量占比低于do时，节点只会发送数据，不会参与其他活动。
```
  do=30;% 判断是否广播的能耗
```

剩下的就是各种能量消耗数值定义。

  Bit=4000;% 每一次传输的bit量Eo=0.02;%初始能量ETX=50*0.000000001;%传输能量，每bitERX=50*0.000000001;%接收能量，每bitEfs=10*0.000000001;%耗散能量，每bitEDA=5*0.000000001;%融合能耗，每bitEmp=5*0.000000001;% 广播能量

%-- 记录节点相关量 --%

记录节点死亡数据

  dead=0;% 死亡节点个数 dead_a=0;% 死亡的特殊节点dead_n=0;% 死亡的普通节点PACKETS_FIRST_DEAD_X=zeros(1,times);% 记录每次实验第一次死掉节点时的横坐标PACKETS_FIRST_DEAD_Y=zeros(1,times);% 记录每次实验第一次死掉节点时的轮数PACKETS_DEAD_X=zeros(1,r);% 该变量为“记录每轮死亡节点死亡数量的值”的变化量的横坐标PACKETS_DEAD_Y=zeros(1,r);% 该变量为“记录每轮死亡节点死亡数量的值”的变化量的纵坐标值

记录簇头节点相关数据。

  cluster=0;% 记录当前的簇头节点的个数countCHs=0;% 记录簇头节点的个数。（跟cluster作用几乎一样，只是起到记录作用，并未像cluster一样进行运算）PACKETS_TO_BS=zeros(1,r);% 记录每一轮的簇头个数

结构体C记录了簇头节点的数据。C的属性有：

  1. xd : 簇头节点的x坐标2. yd : 簇头节点的y坐标3. distance : 簇头节点距离基站的距离4. id : 簇头节点的编号C定义如下：C=repmat(struct('xd',0,'yd',0,'distance',0,'id',0),n*p,1);% 记录当前轮簇头的相关数据

记录普通节点相关

  PACKETS_TO_CH=zeros(1,r);% 记录每一轮的普通节点个数

初始化所有节点

初始变量声明

将显示节点的窗口命名为1。
```
  figure(1);
```

初始化节点，循环n次，画出n个节点。其中前面的n*m个是高级节点，用'+'标记，后面的n-n*m个是普通节点，用'0'标记。

  if(temp_rnd0<=m*n+1) %前面的n*m个是+节点S(i).E=Eo*(1+a);  %能量是初始能量的2倍，a=1 S(i).ENERGY=1;plot(S(i).xd,S(i).yd,'+');hold on;endif(temp_rnd0>=m*n+1)  %后面的n-n*m个是0节点S(i).E=Eo;% S(i).E 是该节点当前的能量值S(i).ENERGY=0;plot(S(i).xd,S(i).yd,'o');% plot函数，在界面上画点hold on;end

LEACH分簇（簇的建立阶段）

初始变量声明

簇头选举初始阶段，会先设置一个建议网络簇头的百分比（注意是建议，每一轮实际选择的簇头节点不一定符合该百分比）
簇头选择阶段，每个节点产生一个0-1的随机数字，如果该数字小于阈值，节点成为当前轮的簇头。

阈值公式为:
```
  T(n) = p / {1 - p * [(r - 1) mod (1 / p)]}其中，p为预期的簇头百分比，r为当前轮数。本次实验中：p = 0.1
```
为了保证节点相对均衡的消耗能量，存在一个循环选举周期。在一个循环选举周期内，每个簇头只能担任一次簇头节点。当一个循环选举周期结束后，所有节点可以重新被选择成为簇头，并且开启新的一个循环选举周期。
```
  这个周期持续的轮数为：n * p , 即 建议网络中簇头存在个数 的值。
```
当i节点在选择中被选为簇头后，将 S(i).type 设置为'C'，表示该节点已经当选过簇头。
当循环周期结束的时候，根据算法可以知道所有的存活的节点都已经担任过一次簇头了，为了使节点可以重新担任簇头，将所有的节点的type全部设置为'N'。
当一个节点当选为簇头后，用'k*'标记。
```
  plot(S(i).xd,S(i).yd,'k*');
```

普通节点加入生成的簇头

初始变量声明

对于每一个普通节点，在一轮的数据传输中，该普通节点将与其距离最近的簇头节点进行通信。所以需要求出其附近据它最近的簇头节点是哪一个。

  %如果i节点是普通节点，并且i节点的能量够，则判断是否将该节点加入簇头if(S(i).type=='N'&& S(i).E>0)if(cluster>=1)min_dis=sqrt((S(i).xd-C(1).xd)^2+(S(i).yd-C(1).yd)^2);% 计算该点距离第一个簇头的距离（默认记录第一个簇头作为初始点）min_dis_cluster=1;% 记录第一个簇头的标号，（默认记录第一个簇头作为初始点）for c=1:1:clustertemp=min(min_dis,sqrt((S(i).xd-C(c).xd)^2+(S(i).yd-C(c).yd)^2)); % 计算c节点距离i节点的距离if(temp<min_dis) % 如果i到c比i到中心点距离短，则选择将数据传送到c节点min_dis=temp;min_dis_cluster=c;endendendend

稳定的数据通信阶段

初始变量声明

簇头节点进行的能量消耗

  if(distance>do)  % 能量大于阈值，消耗能量为：传输能量(将结果发送给汇聚节点)+数据融合能量+与汇聚节点通信能量%S(i).E=S(i).E-((ETX+EDA)*(Bit)+Emp*Bit*(distance*distance*distance*distance));S(i).E=S(i).E-((ETX+EDA)*(Bit)+Emp*Bit);endif (distance<=do)% 能量小于等于阈值，消耗能量为：传输能量(将结果发送给汇聚节点)+数据融合能量+耗散能量%S(i).E=S(i).E-((ETX+EDA)*(Bit)+Efs*Bit*(distance*distance));S(i).E=S(i).E-((ETX+EDA)*(Bit)+Efs*Bit);end

加入簇头节点的普通节点的能量消耗

  if(min_dis>do) % 如果最近的点的剩余能量可以进行广播，则参与活动。% 消耗的能量为：传输能量(将数据发送给簇头)S(i).E=S(i).E-(ETX*Bit);endif(min_dis<=do) % 如果最近的点的剩余能量不足以进行广播，则待机% 消耗的能量为：耗散能量S(i).E=S(i).E-(Efs*Bit);end%计算簇头节点的能量消耗if(min_dis>0) % 如果簇头节点还未能量耗尽，则消耗能量为：接收能量+融合能量。S(C(min_dis_cluster).id).E=S(C(min_dis_cluster).id).E-((ERX+EDA)*Bit);end

当普通节点与簇头节点通信时

  S(i).min_dis=min_dis;% 记录距离其连接的簇头的距离S(i).min_dis_cluster=min_dis_cluster;% 记录其连接的簇头的编号

计算死亡节点

初始变量声明

当节点剩余能量S.E小于或等于0的时候，标记该节点死亡。并用'red.'标记。

  %判断死亡节点。当节点能量为0时，该节点为死亡节点。for i=1:1:nif(S(i).E<=0) plot(S(i).xd,S(i).yd,'red.');hold on;dead=dead+1;if(dead == 1 && isFirst == 1)isFirst = 0;fprintf("times=%d\tr=%d\n",tests,temp_r);%pause;PACKETS_FIRST_DEAD_X(tests)=tests;PACKETS_FIRST_DEAD_Y(tests)=temp_r;endif(S(i).ENERGY==1) % 死亡的特殊节点dead_a=dead_a+1; endif(S(i).ENERGY==0) % 死亡的普通节点dead_n=dead_n+1;endhold on;endend

当所有节点都死亡时，界面1如下图所示：

MATLAB画图

在界面3中画出每一次实验节点死亡的折线图：

  figure(3);plot(PACKETS_DEAD_X,PACKETS_DEAD_Y);xlabel('轮数');ylabel('死亡节点个数');title('节点死亡曲线图');

在界面2中画出当全部实验都结束后，统计每一次实验第一次出现死亡节点的时候，所执行的轮数。

  figure(2);plot(PACKETS_FIRST_DEAD_X,PACKETS_FIRST_DEAD_Y);% 画折线图xlabel('测试次数');ylabel('运行轮数');title('第一次出现死亡节点时的运行轮数');

完整代码

想要其他的效果可以通过更改代码中变量的参数得到。
代码已加入较为详细的注释。
MATLAB代码实现如下：

```MATLAB
%仿真实验环境设定
%本实验设定有10个能量是1个单位的高级节点，用+ 表示;
%有90个能量是0.5个单位的普通节点，用o表示。
%在区域中间(50,50)的坐标处，有一个基站，用x表示.
% (CH:Cluster Head,即簇头)
clear;%清除內存变量%-- 节点、簇头设置相关量 --%
n=100;%节点总数
m=0.1;% 控制高级节点的占比
p=0.1;%簇头概率%-- 控制坐标相关量 --%
xm=100;%x轴范围
ym=100;%y轴范围
sink.x=0.5*xm; %基站x轴
sink.y=0.5*ym; %基站y轴
S=repmat(struct('xd',0,'yd',0,'G',0','type','N','E',0,'ENERGY',0,'min_dis',0,'min_dis_cluster',0),n+1,1);% 记录所有节点的坐标，及其相关属性%-- 计算能量相关量 --%
Bit=4000;% 每一次传输的bit量
Eo=0.02;%初始能量
ETX=50*0.000000001;%传输能量，每bit
ERX=50*0.000000001;%接收能量，每bit
Efs=10*0.000000001;%耗散能量，每bit
EDA=5*0.000000001;%融合能耗，每bit
Emp=5*0.000000001;% 广播能量
a=1;% 初始额外增加的能量占比
do=30;% 判断是否广播的能耗%-- 控制轮数相关量 --%
times=10;% 实验次数
r=500;% 总轮数（在簇头选举初始阶段，每个节点产生一个0-1的随机数字，如果该数字小于阈值，节点成为当前轮的簇头。）%-- 记录节点相关量 --%
dead=0;% 死亡节点个数
dead_a=0;% 死亡的特殊节点
dead_n=0;% 死亡的普通节点
cluster=0;% 记录当前的簇头节点的个数
countCHs=0;% 记录簇头节点的个数。（跟cluster作用几乎一样，只是起到记录作用，并未像cluster一样进行运算）
PACKETS_TO_BS=zeros(1,r);% 记录每一轮的簇头个数
PACKETS_TO_CH=zeros(1,r);% 记录每一轮的普通节点个数
C=repmat(struct('xd',0,'yd',0,'distance',0,'id',0),n*p,1);% 记录当前轮簇头的相关数据
PACKETS_FIRST_DEAD_X=zeros(1,times);% 记录每次实验第一次死掉节点时的横坐标
PACKETS_FIRST_DEAD_Y=zeros(1,times);% 记录每次实验第一次死掉节点时的轮数PACKETS_DEAD_X=zeros(1,r);% 该变量为“记录每轮死亡节点死亡数量的值”的变化量的横坐标
PACKETS_DEAD_Y=zeros(1,r);% 该变量为“记录每轮死亡节点死亡数量的值”的变化量的纵坐标值
for i=1:1:rPACKETS_DEAD_X(i)=i;
endfor tests=1:1:timesfigure(1);% 初始化节点，循环n次，画出n个节点XR=zeros(1,n);% 临时变量，记录X坐标YR=zeros(1,n);% 临时变量，记录Y坐标for i=1:1:n% 新建一个随机的点S(i).xd=rand(1,1)*xm;XR(i)=S(i).xd;S(i).yd=rand(1,1)*ym;YR(i)=S(i).yd;S(i).G=0;  %每一次周期结束此变量为0S(i).type='N';temp_rnd0=i;if(temp_rnd0<=m*n+1) %前面的n*m个是+节点S(i).E=Eo*(1+a);  %能量是初始能量的2倍，a=1 S(i).ENERGY=1;plot(S(i).xd,S(i).yd,'+');hold on;endif(temp_rnd0>=m*n+1)  %后面的n-n*m个是0节点S(i).E=Eo;% S(i).E 是该节点当前的能量值S(i).ENERGY=0;plot(S(i).xd,S(i).yd,'o');% plot函数，在界面上画点hold on;endendS(n+1).xd=sink.x;   %最后一个节点画，表示基站节点S(n+1).yd=sink.y;plot(S(n+1).xd,S(n+1).yd,'x');hold on;% 开始仿真isFirst = 1;for temp_r=1:1:r % 开始新的一轮hold off;% 使新添加到坐标区中的绘图清除现有绘图并重置所有的坐标区属性。figure(1);% 将界面反复显示在标号为1的窗口中% 新的一轮，当前簇头节点个数刷新cluster = 0;countCHs = 0;% 新的一轮，重新记录节点死亡数dead_a = 0;dead_n = 0;dead = 0;% 画点for i=1:1:nif(i<=n*m) % 前面的n*m个是+节点plot(S(i).xd,S(i).yd,'+');endif(i>n*m) % 后面的n-n*m个是0节点plot(S(i).xd,S(i).yd,'o');endhold on;end%判断死亡节点。当节点能量为0时，该节点为死亡节点。for i=1:1:nif(S(i).E<=0) plot(S(i).xd,S(i).yd,'red.');hold on;dead=dead+1;if(dead == 1 && isFirst == 1)isFirst = 0;fprintf("times=%d\tr=%d\n",tests,temp_r);%pause;PACKETS_FIRST_DEAD_X(tests)=tests;PACKETS_FIRST_DEAD_Y(tests)=temp_r;endif(S(i).ENERGY==1) % 死亡的特殊节点dead_a=dead_a+1; endif(S(i).ENERGY==0) % 死亡的普通节点dead_n=dead_n+1;endhold on;endendPACKETS_DEAD_Y(temp_r)=dead;plot(S(n+1).xd,S(n+1).yd,'x');hold on;%选举簇头节点。按 LEACH算法的公式选举簇头节点。 X=zeros(1,n*p);% 临时变量，记录X坐标Y=zeros(1,n*p);% 临时变量，记录Y坐标for i=1:1:nif(S(i).E>0) temp_rand=rand; % 给temp_rand随机附一个值（rand 返回一个(0-1)之间的一个随机数）if((S(i).G)<=0)% 该节点在本轮中没有当过簇头节点的二重保证if(temp_rand<=(p/(1-p*mod ((temp_r-1),round(1/p))))) countCHs=countCHs+1;% packets_T0_BS=packets_T0_BS+1; % PACKETS_TO_BS(temp_r+1)=packets_T0_BS; S(i).type='C';S(i).G=round(1/p)-1; % (round四舍五入为最近的整数)cluster=cluster+1;%簇头數加1C(cluster).xd=S(i).xd;C(cluster).yd=S(i).yd;plot(S(i).xd,S(i).yd,'k*');hold on;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;X(cluster)=S(i).xd;Y(cluster)=S(i).yd; distance;if(distance>do)  % 能量大于阈值，消耗能量为：传输能量(将结果发送给汇聚节点)+数据融合能量+与汇聚节点通信能量%S(i).E=S(i).E-((ETX+EDA)*(Bit)+Emp*Bit*(distance*distance*distance*distance));S(i).E=S(i).E-((ETX+EDA)*(Bit)+Emp*Bit);endif (distance<=do)% 能量小于等于阈值，消耗能量为：传输能量(将结果发送给汇聚节点)+数据融合能量+耗散能量%S(i).E=S(i).E-((ETX+EDA)*(Bit)+Efs*Bit*(distance*distance));S(i).E=S(i).E-((ETX+EDA)*(Bit)+Efs*Bit);endendendendendPACKETS_TO_BS(temp_r)=countCHs;% 本轮产生的簇头节点数%普通节点加入簇中for i=1:1:n%如果i节点是普通节点，并且i节点的能量够，则判断是否将该节点加入簇头if(S(i).type=='N'&& S(i).E>0)if(cluster>=1)min_dis=sqrt((S(i).xd-C(1).xd)^2+(S(i).yd-C(1).yd)^2);% 计算该点距离第一个簇头的距离（默认记录第一个簇头作为初始点）min_dis_cluster=1;% 记录第一个簇头的标号，（默认记录第一个簇头作为初始点）for c=1:1:clustertemp=min(min_dis,sqrt((S(i).xd-C(c).xd)^2+(S(i).yd-C(c).yd)^2)); % 计算c节点距离i节点的距离if(temp<min_dis) % 如果i到c比i到中心点距离短，则选择将数据传送到c节点min_dis=temp;min_dis_cluster=c;endendif(min_dis>do) % 如果最近的点的剩余能量可以进行广播，则参与活动。% 消耗的能量为：传输能量(将数据发送给簇头)S(i).E=S(i).E-(ETX*Bit);endif(min_dis<=do) % 如果最近的点的剩余能量不足以进行广播，则待机% 消耗的能量为：耗散能量S(i).E=S(i).E-(Efs*Bit);end%计算簇头节点的能量消耗if(min_dis>0) % 如果簇头节点还未能量耗尽，则消耗能量为：接收能量+融合能量。S(C(min_dis_cluster).id).E=S(C(min_dis_cluster).id).E-((ERX+EDA)*Bit);endS(i).min_dis=min_dis;% 记录距离其连接的簇头的距离S(i).min_dis_cluster=min_dis_cluster;% 记录其连接的簇头的编号endendendPACKETS_TO_CH(temp_r)=n-dead-cluster;% 记录当前存活非簇头节点个数fprintf("%d %d\n",temp_r, countCHs);% 当所有的节点都当过一次簇头后，将所有节点都设置为未曾担任过簇头节点，以便开始新的簇头选举if(mod (temp_r,round(n*p)) == 0)for i=1:1:nS(i).G=0;  %每一次周期结束此变量为0S(i).type='N';endendendfigure(3);plot(PACKETS_DEAD_X,PACKETS_DEAD_Y);xlabel('轮数');ylabel('死亡节点个数');title('节点死亡曲线图');hold on;
endfigure(2);
plot(PACKETS_FIRST_DEAD_X,PACKETS_FIRST_DEAD_Y);% 画折线图
xlabel('测试次数');
ylabel('运行轮数');
title('第一次出现死亡节点时的运行轮数');```