http://www.ilovematlab.cn/thread-177006-1-1.html

LEACH協議
clear;%清除內存變量

xm=100;%x軸範圍
ym=100;%y軸範圍

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
    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;%第一個死亡節點的標誌變量

for r=1:1:rmax%開始每輪循環
  r+1%顯示輪數
    if(mod(r,round(1/p))==0)
       for i=1:1:n
           S(i).G=0;
       end
    end%如何輪數正好是一個週期的整數倍，則設置S(i).E為0

     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;%將能量小於等於0的節點繪製成紅色，並將死亡節點數增加1
         
         if(dead==1)
           if(flag_first_dead==0)
             first_dead=r %第一個節點的死亡輪數
             save ltest, first_dead;
             flag_first_dead=1;
           end
         end%將能量小於等於0的節點繪製成紅色，並將死亡節點數增加1
           
         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;%發送到基站的數據包數為1
            end
           end
          end
         end

        CH_Num(r+1)=cluster; %每輪的簇頭數
        save ltest,CH_Num(r+1);%保存每輪簇頭數到ltest
     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;
           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(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
  
    for i=1:1:n
    R(r+1,i)=S(i).E;  %每輪每節點的剩餘能量
    % save ltest,R(r+1,i);%保存此數據到ltest
    end

    hold on;

end