LEACH剩余能量图、存活节点图

clear;xm=100;%设置区域为100*100ym=100;sink.x=0.5*xm;%sink（汇聚）节点坐标sink.y=0.5*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=1500%最大的轮数do=sqrt(Efs/Emp);%计算do通信半径。figure(1);%输出图形fori=1:1:n%i为矩阵1到n，间距为1S(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;%设置初始能量为E0S(i).ENERGY=0;%普通节点标志plot(S(i).xd,S(i).yd,'o');%输出节点，用o表示holdon;endS(n+1).xd=sink.x;%汇聚节点X轴坐标S(n+1).yd=sink.y;%汇聚节点Y轴坐标plot(S(n+1).xd,S(n+1).yd,'x');%输出汇聚节点，用x表示%第一次迭代figure(1);cluster=1;flag_first_dead=0;%第一个节点死亡的标志变量forr=0:1:rmaxrif(mod(r,round(1/p))==0)%如果所有节点都当过簇头，则全部节点清零，回到最初状态fori=1:1:nS(i).G=0;%簇头数目endenddead=0;figure(4);fori=1:1:nif(S(i).E=0)%检查是否有节点死亡plot(S(i).xd,S(i).yd,'red.')%输出节点，用红.表示dead=dead+1;%节点死亡数+1holdon;endif(S(i).E0)%节点能量大于0S(i).type='N';plot(S(i).xd,S(i).yd,'o');holdon;endendplot(S(n+1).xd,S(n+1).yd,'x');%sinkSTATISTICS_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;%第一个死亡节点标志endendcountCHs=0;%簇头的个数cluster=1;%簇头的标号，初始值为1fori=1:1:n%i为矩阵1到n，间距为1if(S(i).E0)%节点剩余能量大于0temp_rand=rand;if((S(i).G)=0)%没有当过簇头？if(temp_rand=(p/(1-p*mod(r,round(1/p)))))%对簇头节点进行处理countCHs=countCHs+1;%簇头数+1S(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(distancedo)%距离大于通信半径S(i).E=S(i).E-((ETX+EDA)*(4000)+Emp*4000*(distance*distance*distance*distance));%能量消耗endif(distance=do)%距离小于通信半径S(i).E=S(i).E-((ETX+EDA)*(4000)+Efs*4000*(distance*distance));%能量消耗endendendendendSTATISTICS(r+1).CLUSTERHEADS=cluster-1;%r轮后簇头数CLUSTERHS(r+1)=cluster-1;%r轮后簇头数fori=1:1:nif(S(i).type=='N'&&S(i).E0)%处理普通节点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;%初始距离最近的簇头默认为1forc=1:1:cluster-1temp=min(min_dis,sqrt((S(i).xd-C(c).xd)^2+(S(i).yd-C(c).yd)^2));%选取节点到簇头距离和到sink节点距离之间最近的一个if(tempmin_dis)%即离某个簇头距离更近min_dis=temp;min_dis_cluster=c;endend%循环结束后，即找到该节点对应的最近簇头或者sinkmin_dis;if(min_disdo)S(i).E=S(i).E-(ETX*(4000)+Emp*4000*(min_dis*min_dis*min_dis*min_dis));endif(min_dis=do)S(i).E=S(i).E-(ETX*(4000)+Efs*4000*(min_dis*min_dis));endif(min_dis0)%能量消散S(C(min_dis_cluster).id).E=S(C(min_dis_cluster).id).E-((ERX+EDA)*4000);endS(i).min_dis=min_dis;S(i).min_dis_cluster=min_dis_cluster;endendendSTATISTICS(r+1).ENERGY=0;fori=1:1:nifS(i).E0STATISTICS(r+1).ENERGY=STATISTICS(r+1).ENERGY+S(i).E;%r轮后节点剩余能量加上节点i的剩余能量endendholdoff;endfori=2:1:rmax%当前节点数mylive(i)=n-STATISTICS_leach(i).DEAD;myenergy(i)=STATISTICS(i).ENERGY;%剩余能量endmylive(1)=100;myenergy(1)=S(1).E+(n-1)*Eo;figure(2);%输出图形2holdon;%保持曲线plot(mylive,'color','r');%用红色输出存活节点数xlabel('周期数');ylabel('存活节点');title('存活节点图');figure(3);%输出图形3holdon;%保持曲线plot(myenergy,'color','r');%用红色输出剩余能量