%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                      %% SEP: A Stable Election Protocol for clustered                        %%      heterogeneous wireless sensor networks                          %%                                                                      %% (c) Georgios Smaragdakis                                             %% WING group, Computer Science Department, Boston University           %%                                                                      %% You can find full documentation and related information at:          %% http://csr.bu.edu/sep                                                %%                                                                      %  % To report your comment or any bug please send e-mail to:             %% gsmaragd@cs.bu.edu                                                   %%                                                                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                      %% This is the LEACH [1] code we have used.                             %% The same code can be used for FAIR if m=1                            %%                                                                      %% [1] W.R.Heinzelman, A.P.Chandrakasan and H.Balakrishnan,             %%     "An application-specific protocol architecture for wireless      % %      microsensor networks"                                           % %     IEEE Transactions on Wireless Communications, 1(4):660-670,2002  %%                                                                      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clear;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%Field Dimensions - x and y maximum (in meters)xm=100;ym=100;%x and y Coordinates of the Sinksink.x=0.5*xm;sink.y=0.5*ym;%Number of Nodes in the fieldn=100%Optimal Election Probability of a node%to become cluster headp=0.1;%Energy Model (all values in Joules)%Initial Energy Eo=0.5;%Eelec=Etx=ErxETX=50*0.000000001;ERX=50*0.000000001;%Transmit Amplifier typesEfs=10*0.000000000001;Emp=0.0013*0.000000000001;%Data Aggregation EnergyEDA=5*0.000000001;%Values for Hetereogeneity%Percentage of nodes than are advancedm=0.1;%\alphaa=1;%maximum number of roundsrmax=8Nmin=5%according Rc ,%%%%%%%%%%%%%%%%%%%%%%%%% END OF PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%Computation of dodo=sqrt(Efs/Emp);%Creation of the random Sensor Networkfigure(1);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;    %initially there are no cluster heads only nodes    S(i).type='N';       temp_rnd0=i;    %Random Election of Normal Nodes    if (temp_rnd0>=m*n+1)         S(i).E=Eo;        S(i).ENERGY=0;        plot(S(i).xd,S(i).yd,'o');        hold on;    end    %Random Election of Advanced Nodes    if (temp_rnd0<m*n+1)          S(i).E=Eo*(1+a)        S(i).ENERGY=1;        plot(S(i).xd,S(i).yd,'+');        hold on;    endendS(n+1).xd=sink.x;S(n+1).yd=sink.y;plot(S(n+1).xd,S(n+1).yd,'x');            %First Iterationfigure(1);%counter for CHscountCHs=0;%counter for CHs per roundrcountCHs=0;cluster=1;countCHs;rcountCHs=rcountCHs+countCHs;flag_first_dead=0;for r=0:1:rmax    r  %Operation for epoch  if(mod(r, round(1/p) )==0)    for i=1:1:n        S(i).G=0;        S(i).cl=0;    end  endhold off;%Number of dead nodesdead=0;%Number of dead Advanced Nodesdead_a=0;%Number of dead Normal Nodesdead_n=0;%counter for bit transmitted to Bases Station and to Cluster Headspackets_TO_BS=0;packets_TO_CH=0;%counter for bit transmitted to Bases Station and to Cluster Heads %per roundPACKETS_TO_CH(r+1)=0;PACKETS_TO_BS(r+1)=0;figure(1);for i=1:1:n    %checking if there is a dead node    if (S(i).E<=0)        plot(S(i).xd,S(i).yd,'red .');        dead=dead+1;        if(S(i).ENERGY==1)            dead_a=dead_a+1;        end        if(S(i).ENERGY==0)            dead_n=dead_n+1;        end        hold on;        end    if S(i).E>0        S(i).type='N';        if (S(i).ENERGY==0)          plot(S(i).xd,S(i).yd,'o');        end        if (S(i).ENERGY==1)          plot(S(i).xd,S(i).yd,'+');        end        hold on;    endendplot(S(n+1).xd,S(n+1).yd,'x');STATISTICS(r+1).DEAD=dead;DEAD(r+1)=dead;DEAD_N(r+1)=dead_n;DEAD_A(r+1)=dead_a;%When the first node diesif (dead==1)    if(flag_first_dead==0)        first_dead=r        flag_first_dead=1;    endend%Rc=35;%半徑countCHs=0;cluster=1;%lczC(1).xd=50;C(1).yd=50; C(1).distance=0;C(1).id=n+1;C(1).Ncount=0;%lczfor i=1:1:n   if(S(i).E>0)   temp_rand=rand;        if ( (S(i).G)<=0)min_temp1=36;% %Election of Cluster Heads if(temp_rand<= (p/(1-p*mod(r,round(1/p)))))        %均勻分布Cluster        for j=1:1:cluster          min_temp1           %temp=;          min_temp1=min(min_temp1,sqrt( (C(j).xd-S(i).xd )^2 + (C(j).yd- S(i).yd)^2))        end if(min_temp1>=Rc/4)        %if(min(c_distance)>=Rc)        %均勻分布Cluster            countCHs=countCHs+1;            packets_TO_BS=packets_TO_BS+1;            PACKETS_TO_BS(r+1)=packets_TO_BS;                        S(i).type='C';            S(i).G=round(1/p)-1;            C(cluster+1).xd=S(i).xd;            C(cluster+1).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;lcz  all next 4 line ,cluster            %convert cluster+1            C(cluster+1).distance=distance;            %Ncount will count the Associating Normal nodes for this            %Cluster            C(cluster+1).Ncount=0;            C(cluster+1).id=i;            X(cluster+1)=S(i).xd;            Y(cluster+1)=S(i).yd;            cluster=cluster+1;                        %add circle            alpha=0:pi/20:2*pi;%角度[0,2*pi]           % Rc=35;%半徑            Xc=S(i).xd+Rc*cos(alpha);            Yc=S(i).yd+Rc*sin(alpha);            plot(Xc,Yc)            axis equal                    %Calculation of Energy dissipated            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 endSTATISTICS(r+1).CLUSTERHEADS=cluster-1;CLUSTERHS(r+1)=cluster-1;%Election of Associated Cluster Head for Normal Nodesfor i=1:1:n   if ( S(i).type=='N' && S(i).E>0 )     if(cluster-1>=1)    %   if(cluster-1>=1)       min_dis=sqrt( (S(i).xd-S(n+1).xd)^2 + (S(i).yd-S(n+1).yd)^2 );       min_dis_temp=min_dis;       min_dis_cluster=1;       for c=1:1:cluster%cluster-1           temp=min(min_dis,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       C(min_dis_cluster).Ncount=C(min_dis_cluster).Ncount+1;       if (min_dis_temp==min_dis)         plot([S(i).xd,S(n+1).xd],[S(i).yd,S(n+1).yd])         else         plot([S(i).xd,C(min_dis_cluster).xd],[S(i).yd,C(min_dis_cluster).yd])       end       %Energy dissipated by associated Cluster Head  %          min_dis; %%          if (min_dis>do)    %            S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis));   %              lczflag=1;  %        end  %          if (min_dis<=do)  %              S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis));   %         end        %Energy dissipated   %     if(min_dis>0)   %       S(C(min_dis_cluster).id).E = S(C(min_dis_cluster).id).E- ( (ERX + EDA)*4000 );   %       PACKETS_TO_CH(r+1)=n-dead-cluster+1;   %      end    S(i).min_dis=min_dis;   S(i).min_dis_cluster=min_dis_cluster;           end endendhold on;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%sleep(r+1)=0;for i=1:1:n   if ( S(i).type=='N' && S(i).E>0 )       if(cluster-1>=1)       min_dis=S(i).min_dis          if (S(i).min_dis>Rc)% if the node out of round of Cluster,then is work node            if (S(i).min_dis>do)                S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis));                 lczflag=1;            end            if (S(i).min_dis<=do)                S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis));             end        else          temp_rand=rand;               if(temp_rand<= Nmin/C(S(i).min_dis_cluster).Ncount)                      if (S(i).min_dis>do)                S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis));                 lczflag=1;            end            if (S(i).min_dis<=do)                S(i).E=S(i).E- ( ETX*(4000) + Efs*4000*( min_dis * min_dis));             end          else            sleep(r+1)=sleep(r+1)+1;          end      end%according if the node out of range of Cluster,then is work node        %Energy dissipated        if(S(i).min_dis>0)          S(C(S(i).min_dis_cluster).id).E = S(C(S(i).min_dis_cluster).id).E- ( (ERX + EDA)*4000 );          PACKETS_TO_CH(r+1)=n-dead-cluster+1;         end              end               endend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%countCHs;rcountCHs=rcountCHs+countCHs;%Code for Voronoi Cells%Unfortynately if there is a small%number of cells, Matlab's voronoi%procedure has some problems%[vx,vy]=voronoi(X,Y);%plot(X,Y,'r*',vx,vy,'b-');% hold on;% voronoi(X,Y);% axis([0 xm 0 ym]);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   STATISTICS    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                                                                                     %%  DEAD  : a rmax x 1 array of number of dead nodes/round %  DEAD_A : a rmax x 1 array of number of dead Advanced nodes/round%  DEAD_N : a rmax x 1 array of number of dead Normal nodes/round%  CLUSTERHS : a rmax x 1 array of number of Cluster Heads/round%  PACKETS_TO_BS : a rmax x 1 array of number packets send to Base Station/round%  PACKETS_TO_CH : a rmax x 1 array of number of packets send to ClusterHeads/round%  first_dead: the round where the first node died                   %                                                                                     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%