无线传感网络LEACH的MATLAB实现

xm=100;%设置区域为100*100
ym=100;

sink.x=0.5*xm;%sink（汇聚）节点坐标
sink.y=1.75*ym;

n=100 %区域内的节点数目
p=0.05;% 节点成为簇头的概率 
Eo=0.5;%节点初始能量
ETX=50*0.000000001;%发射单位报文损耗能量 
ERX=50*0.000000001;%接收单位报文损耗能量
Efs=10*0.000000000001;%自由空间能量
Emp=0.0013*0.000000000001;%衰减空间能量
EDA=5*0.000000001;%聚集数据所要消耗的能量
rmax=2000%最大的轮数
do=sqrt(Efs/Emp); %计算do 通信半径。

figure(1);%输出图形
for i=1:1:n %i为矩阵1到n，间距为1
    S(i).xd=rand(1,1)*xm;%1行1列矩阵
    XR(i)=S(i).xd;%随机生成的X轴
    S(i).yd=rand(1,1)*ym;
    YR(i)=S(i).yd;%随机生成的Y轴
    S(i).G=0;
    S(i).type='N';%节点类型为普通
    S(i).E=Eo;%设置初始能量为E0
    S(i).ENERGY=0;%普通节点标志
    plot(S(i).xd,S(i).yd,'o');%输出节点，用o表示
    hold on;
end

S(n+1).xd=sink.x;%汇聚节点X轴坐标
S(n+1).yd=sink.y;%汇聚节点Y轴坐标
plot(S(n+1).xd,S(n+1).yd,'x'); %输出汇聚节点，用x表示 
%第一次迭代
title('Wireless Sensor Network');
figure(1);
cluster=1;
flag_first_dead=0;%第一个节点死亡的标志变量

for r=0:1:rmax
r
if(mod(r, round(1/p) )==0)%如果所有节点都当过簇头，则全部节点清零，回到最初状态
    for i=1:1:n
        S(i).G=0;%簇头数目
    end
end
dead=0;

figure(4);

for i=1:1:n 
    if (S(i).E<=0)%检查是否有节点死亡
        plot(S(i).xd,S(i).yd,'red .')%输出节点，用红.表示
        dead=dead+1;%节点死亡数+1
        hold on; 
    end
    if (S(i).E>0)%节点能量大于0
        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');%sink


STATISTICS_leach(r+1).DEAD=dead;%r轮后死亡节点数
DEAD_leach(r+1)=dead;%r轮后死亡节点数
if (dead==1)%第一个节点死亡
    if(flag_first_dead==0)%第一个节点死亡周期
        first_dead=r%第一个节点死亡轮数
        flag_first_dead=1;%第一个死亡节点标志
    end
end
countCHs=0;%簇头的个数
cluster=1;%簇头的标号，初始值为1
for i=1:1:n%i为矩阵1到n，间距为1
    if(S(i).E>0)%节点剩余能量大于0
        temp_rand=rand; 
        if ( (S(i).G)<=0)%没有当过簇头？
            if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))%对簇头节点进行处理
                countCHs=countCHs+1;%簇头数+1
                S(i).type='C';%节点类型为簇头
                S(i).G=100;%S(i).G=100表示当过簇头
                C(cluster).xd=S(i).xd;%簇头X轴坐标 
                C(cluster).yd=S(i).yd;%簇头Y轴坐标
                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 );%相对应的簇头到sink的距离
                C(cluster).distance=distance;%距离
                C(cluster).id=i;%簇头对应的节点编号
                X(cluster)=S(i).xd;%X轴坐标
                Y(cluster)=S(i).yd;%Y轴坐标
                cluster=cluster+1;%簇头标号数+1！！
                
                distance;
                if (distance>do)%距离大于通信半径
                    S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance )); %能量消耗
                end
                if (distance<=do)%距离小于通信半径
                    S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Efs*4000*( distance * distance )); %能量消耗
                end
            end 
        end
    end 
end
STATISTICS(r+1).CLUSTERHEADS=cluster-1;%r轮后簇头数
CLUSTERHS(r+1)=cluster-1;%r轮后簇头数

for i=1:1:n
    if ( S(i).type=='N' && S(i).E>0 )%处理普通节点
        if(cluster-1>=1)%簇头总数大于2个
            min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );%普通节点到sink节点间最短距离
            min_dis_cluster=1;%初始距离最近的簇头默认为1
            for c=1:1:cluster-1
                temp=min(min_dis,sqrt( (S(i).xd-C(c).xd)^2 + (S(i).yd-C(c).yd)^2 ) );%选取节点到簇头距离和到sink节点距离之间最近的一个
                if ( temp<min_dis )%即离某个簇头距离更近
                    min_dis=temp;
                    min_dis_cluster=c;
                end
            end%循环结束后，即找到该节点对应的最近簇头或者sink
            min_dis;
            if (min_dis>do)
                S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis)); 
            end
            if (min_dis<=do)
                S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis)); 
            end

            if(min_dis>0)%能量消散
                S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (ERX + EDA)*4000 ); 
            end
            S(i).min_dis=min_dis;
            S(i).min_dis_cluster=min_dis_cluster;
        end
    end
end
STATISTICS(r+1).ENERGY=0;
for i=1:1:n
    if S(i).E > 0
        STATISTICS(r+1).ENERGY = STATISTICS(r+1).ENERGY +S(i).E;%r轮后节点剩余能量加上节点i的剩余能量
    end
end

hold off;
end
for i=2:1:rmax%当前节点数
mylive(i) = n - STATISTICS_leach(i).DEAD;
myenergy(i) = STATISTICS(i).ENERGY;%剩余能量
end
mylive(1)=100;
myenergy(1)=S(1).E+(n-1)*Eo;
figure(2);%输出图形2
hold on;%保持曲线
plot(mylive,'color','r');%用红色输出存活节点数
xlabel('Rounds');
ylabel('Operation Nodes');
title('LEACH Operating Nodes per Round');
figure(3);%输出图形3
hold on;%保持曲线
plot(myenergy,'color','r');%用红色输出剩余能量