無線傳感器網(wǎng)絡的leach算法及其改進

1、無線傳感器網(wǎng)絡LEACH算法的研究與改進楊林(華南理工大學自動化科學與工程學院, 廣東廣州)摘要：本文重點研究基于無線傳感器網(wǎng)絡的LEACH和SEP式路由協(xié)議，無線傳感器網(wǎng)絡節(jié)點數(shù)量龐大、單個節(jié)點資源有限，其路由協(xié)議設計的首要目標是提高能量有效性，延長網(wǎng)絡壽命。本文選擇了LEACH路由協(xié)議為研究重點，分析了該路由算法的具體實現(xiàn)，針對傳感器節(jié)點能量及網(wǎng)絡生存周期等特點，介紹了其改進型SEP路由算法。并采用MATLAB網(wǎng)絡仿真工具對該路由協(xié)議進行了整體仿真，并對其數(shù)據(jù)進行了分析。關鍵詞：無線傳感器網(wǎng)絡；LEACH路由協(xié)議；SEP路由協(xié)議；MATLAB仿真Improvement of leach a

2、lgorithm for wireless sensor networksYang Lin(College of Automation Science and Technology, South China University of Technology, Guangzhou Guangdong , China) Abstract: This paper focuses on LEACH and SEP routing protocol based on wireless sensor networks, wireless sensor network node a huge number,

3、 the limited resources of a single node, the routing protocol design goal is to improve energy efficiency, and prolong the network lifetime .This paper studies the LEACH routing protocol analysis with a focus on the concrete realization of the routing algorithm for sensor node energy and network lif

4、e cycle characteristics, its improved SEP routing algorithm.Key words: WSN; LEACH routing protocols; SEP routing protocols;MATLAB Simulation 1 引言無線傳感器網(wǎng)絡(Wireless Sensor Network, WSN)是目前國際上研究的熱點，它融合了計算、通信和傳感器這 3 項技術的交叉應用，具有十分廣闊的應用前景。在網(wǎng)絡中，數(shù)據(jù)的傳輸是靠路由協(xié)議控制管理的。因此，要設計一個優(yōu)良的網(wǎng)絡系統(tǒng)離不開路由協(xié)議的研究。無線傳感器網(wǎng)絡的特點有：(1)節(jié)點沒有全

5、球唯一的標識符，傳統(tǒng)的互聯(lián)網(wǎng)路由協(xié)議無法應用在傳感器網(wǎng)絡中。(2)傳感器網(wǎng)絡中的所有節(jié)點都是源節(jié)點，向唯一的目的節(jié)點 Sink 發(fā)送數(shù)據(jù)。(3)由于在被測對象內部或附近部署了大量的節(jié)點，它們采集到的數(shù)據(jù)是相同或相近的。這就需要路由協(xié)議具有數(shù)據(jù)融合能力。(4)節(jié)點具備的處理能力、電能和存儲能力很有限。開發(fā)針對WSN的路由協(xié)議成為目前無線傳感器網(wǎng)絡研究的熱點和難點。2 WSN簡介2.1 WSN概念無線傳感器網(wǎng)絡（wireless sensor network）簡稱WSN，是一種由大量小型傳感器所組成的網(wǎng)絡。這些小型傳感器一般稱作sensor node（傳感器節(jié)點）或者mote（灰塵）。此種網(wǎng)絡中一

6、般也有一個或幾個基站(稱作sink)用來集中從小型傳感器收集的數(shù)據(jù)。2.2 WSN的結構2.2.1節(jié)點結構在不同應用中，傳感器節(jié)點的結構不盡相同，但一般都由傳感器模塊、處理器模塊、無線通信模塊和能量供應模塊四部分組成，如圖2.1所示。圖2.1 傳感器網(wǎng)絡節(jié)點結構Figure 2.1 Node structure of sensor network傳感器模塊負責監(jiān)測區(qū)域內信息的采集和數(shù)據(jù)轉換，傳感器的類型是由被監(jiān)測物理信號的形式?jīng)Q定的，如用于溫度監(jiān)測的鉑電阻傳感器，用于壓力傳感的電容式傳感器等;處理器模塊負責控制整個傳感器節(jié)點的操作，存儲和處理本身采集的數(shù)據(jù)以及其他節(jié)點發(fā)送來的數(shù)據(jù);無線通信模塊

7、負責與其他傳感器節(jié)點進行無線通信，交換控制信息和收發(fā)采集數(shù)據(jù);能量供應模塊為傳感器節(jié)點提供運行所需的能量。2.2.2 網(wǎng)絡體系結構無線傳感器網(wǎng)絡典型的體系結構如圖2-2所示，包括大量隨機分布的傳感器節(jié)點(sensor node)、匯聚節(jié)點(sink node)、互聯(lián)網(wǎng)(internet)、用戶終端(user)等。散布在監(jiān)測環(huán)境中的傳感器節(jié)點通過網(wǎng)絡自組織多跳方式，將數(shù)據(jù)向匯聚節(jié)點發(fā)送，匯聚節(jié)點可以使用多種方式與外部網(wǎng)絡通信，如Internet衛(wèi)星或移動通信網(wǎng)絡等等，將整個區(qū)域內的數(shù)據(jù)傳送到用戶終端，用戶通過管理節(jié)點對無線傳感器網(wǎng)絡進行配置管理。 圖2.2 無線傳感器網(wǎng)絡體系結構Figure 2

8、.2 Structure of wireless sensor network3 LEACH 路由協(xié)議的分析與研究LEACH 是第一個在無線傳感器網(wǎng)絡中提出的層次式路由協(xié)議。其后的很多層次式路由協(xié)議都是在它基礎上提出來的。因此這里我們著重介紹LEACH 協(xié)議。3.1 LEACH協(xié)議模型LEACH是MIT的Chandrakasan等人為無線傳感器網(wǎng)絡設計的低功耗自適應聚類路由算法。與一般的平面多跳路由協(xié)議和靜態(tài)聚類算法相比,LEACH可以將網(wǎng)絡生命周期延長15%,主要通過隨機選擇聚類首領,平均分擔中繼通信業(yè)務來實現(xiàn)。LEACH定義了“輪”(round)的概念,一輪由初始化和穩(wěn)定工作兩個階段組成。

9、為了避免額外的處理開銷,穩(wěn)定態(tài)一般持續(xù)相對較長的時間。如圖3.1所示圖3.1 LEACH協(xié)議的時序圖Figure 3.1 The sequence diagram of LEACH protocol在初始化階段,聚類首領是通過下面的機制產生的。傳感器節(jié)點生成0,1之間的隨機數(shù),如果大于閾值T,則選該節(jié)點為聚類首領T的計算方法如下： （1）其中p為節(jié)點中成為聚類首領的百分數(shù),r是當前的輪數(shù)。當簇頭選定之后，簇頭節(jié)點主動向網(wǎng)絡中節(jié)點廣播自己成為簇頭的消息(ADV_CH)。接收到此消息的節(jié)點，依據(jù)接收信號的強度，選擇它所要加入的簇，并發(fā)消息通知相應的簇頭(JOIN_REQ)?；跁r分多址(Time

10、Division Multiple Address，簡稱TDMA)的方式，簇頭節(jié)點為其中的每個成員分配通信時隙，并以廣播的形式通知所有的簇內節(jié)點(ADVSCH)。這樣保證了簇內每個節(jié)點在指定的傳輸時隙進行數(shù)據(jù)傳輸，而在其他時間進入休眠狀態(tài)，減少了能量消耗。在穩(wěn)定工作階段，節(jié)點持續(xù)采集監(jiān)測數(shù)據(jù)，在自身傳輸時隙到來時把監(jiān)測數(shù)據(jù)傳給簇頭節(jié)點(DATA)，如圖4.2所示。簇頭節(jié)點對接收到數(shù)據(jù)進行融合處理之后，發(fā)送到Sink節(jié)點，這是一種減小通信業(yè)務量的合理工作模式。持續(xù)一段時間以后，整個網(wǎng)絡進入下一輪工作周期，重新選擇簇頭節(jié)點。圖3.2 LEACH協(xié)議過程Figure 3.2 Progress of

11、the LEACH protocalsLEACH協(xié)議采用動態(tài)轉換簇頭的方法來平均網(wǎng)絡節(jié)點的能量消耗，使因能量耗盡而失效的節(jié)點呈隨機分布狀態(tài)，因而與一般的多跳路由協(xié)議和靜態(tài)簇算法相比，LEACH可以將網(wǎng)絡生命周期延長15%。但是LEACH協(xié)議在每輪固定簇頭節(jié)點后在劃分簇的過程中，簇頭節(jié)點開銷較大。并且簇頭節(jié)點的選擇無法達到最優(yōu)，有可能簇頭節(jié)點位于網(wǎng)絡的邊緣或者幾個簇頭節(jié)點相鄰，某些節(jié)點不得不傳輸較遠的距離來與簇頭通信，這就導致了大量能量消耗。而且LEACH協(xié)議所有簇頭節(jié)點直接與Sink節(jié)點通信，采用連續(xù)數(shù)據(jù)發(fā)送模式和單跳路徑選擇模式，使得每輪中簇頭節(jié)點能耗巨大，因此不適合在大規(guī)模的傳感器網(wǎng)絡中應

12、用。3.2 LEACH協(xié)議的優(yōu)缺點LEACH算法是第一個針對無線傳感器網(wǎng)絡提出的層次型網(wǎng)絡拓撲組織算法，與節(jié)點直接向BS發(fā)送數(shù)據(jù)的方式相比，在性能上有很大提升，但是該算法在某些方面還是存在一定局限性的。(1)優(yōu)點LEACH協(xié)議是一種層次結構的路由協(xié)議，簇首形成高一層的網(wǎng)絡，這樣簇內成員的功能就變得相對簡單，并且不需要維護復雜的路由信息，大大減少路由控制信息的數(shù)量。另外，LEACH算法隨機選取節(jié)點作為簇首，并且簇首是輪換選舉的，這樣能量的高消耗就平均分配到網(wǎng)絡的所有節(jié)點上，以此來延長網(wǎng)絡生命周期。最后，LEACH協(xié)議中簇的組織形式使網(wǎng)絡具有很好的擴展性。(2)缺陷LEACH協(xié)議中簇首不僅要接收簇

13、內成員節(jié)點采集的數(shù)據(jù)并對其進行融合，而且還要將數(shù)據(jù)發(fā)送給BS，因此消耗的能量比普通的成員節(jié)點多，很容易失效，從而導致頻繁分簇，分簇過程中所消耗的能量對于整個網(wǎng)絡中的能耗是一種額外的頭開銷，如果頻繁分簇的話，所產生的頭開銷就會增加，所以網(wǎng)絡節(jié)點能量的利用率就會相應降低。另外，LEACH協(xié)議只適用于小規(guī)模的無線傳感器網(wǎng)絡。由于簇內成員將原始數(shù)據(jù)經(jīng)過單跳直接發(fā)送給簇首，簇首將融合后的數(shù)據(jù)經(jīng)過單跳再發(fā)送給BS，所以簇內成員與簇首、簇首與BS之間必須布置在通信可達的范圍內。3.3 LEACH的改進型SEP協(xié)議經(jīng)典的聚類協(xié)議假定所有的節(jié)點都配備了相同數(shù)量的能量，因此，他們不能充分利用節(jié)點異質性的存在。本文

14、推薦一種能夠異構感知的SEP協(xié)議，該協(xié)議可以延長第一個節(jié)點的死亡時間（也就是穩(wěn)定期）。穩(wěn)定期的延長對于許多應用程序是至關重要的，因為它們從傳感器網(wǎng)絡的反饋要求必須是可靠的。為了延長穩(wěn)定期，SEP協(xié)議試圖維持均衡的能源消耗。直觀上，SEP協(xié)議高級節(jié)點（初始能量高的節(jié)點）成為簇頭的概率大于普通的節(jié)點（初始能量低的節(jié)點），這就相當于對于能源消耗的約束。SEP協(xié)議假定每個節(jié)點知道網(wǎng)絡的總能量，然后根據(jù)節(jié)點的剩余能量計算出成為簇頭的最佳概率。開始時先給每個節(jié)點分配一個權重對應于最優(yōu)概率，權重必須等于每個節(jié)點的初始能量除以普通節(jié)點的初始能量。并且還定義為普通節(jié)點加權選舉的概率，為高級節(jié)點加權選舉的概率。

15、（2） （3）其中a為高級節(jié)點的初始能量是普通節(jié)點初始能量的倍數(shù)，m為高級節(jié)點在總節(jié)點數(shù)中所占比例。普通節(jié)點與高級節(jié)點成為簇頭的閥值分別為和，計算公式如下： （4） （5）其中r是當前輪數(shù)。從而可以保證剩余能量高的節(jié)點成為簇頭的概率大于剩余能量低的節(jié)點。與傳統(tǒng)的LEACH相比有如下優(yōu)勢：（1） 不同節(jié)點的初始能量的不同，SEP協(xié)議選擇能量大的為簇頭從而延長了第一節(jié)點的死亡時間，也就是延長了穩(wěn)定期。（2） 由于穩(wěn)定期的延長，SEP協(xié)議的吞吐量也高于其他聚類協(xié)議。（3） 在一個理想的階段SEP協(xié)議的能量均勻的分布在傳感器網(wǎng)絡中的所有節(jié)點上。4 LEACH與SEP協(xié)議MATLAB仿真比較4.1 LE

16、ACH與SEP仿真結果我們在MATLAB編程環(huán)境中首先產生一個100100的區(qū)域，并在其內部隨機生成一個含有100個節(jié)點（坐標不同）的連通圖。假設每個節(jié)點在網(wǎng)絡形成前已知需要傳感的數(shù)據(jù)類型，且能檢測到出現(xiàn)在其檢測范圍內的事件。Sink在網(wǎng)絡區(qū)域內，不考慮Sink的能耗。普通節(jié)點初始初始能量為40mJ。MATLAB中的實驗結果圖如圖4.1a4.1d所示。圖4.1a 節(jié)點分布圖Figure 4.1a Distribution of nodes圖4.1b 存活節(jié)點數(shù)量圖Figure 4.1b The amount of alive nodes圖4.1c 剩余能量圖Figure 4.1c Remain

17、ing energy圖4.1d 最終節(jié)點分布圖Figure 4.1d Final distribution of nodes為具體比較LEACH與SEP協(xié)議性能，我們在MATLAB編程環(huán)境中首先產生一個100100的區(qū)域，并在其內部隨機生成一個含有100個節(jié)點（坐標不同）的連通圖，而且隨機選擇100*m個節(jié)點作為高級節(jié)點。假設每個節(jié)點在網(wǎng)絡形成前已知需要傳感的數(shù)據(jù)類型，且能檢測到出現(xiàn)在其檢測范圍內的事件。Sink在網(wǎng)絡區(qū)域內，不考慮Sink的能耗。普通節(jié)點初始初始能量為40mJ,100*m個節(jié)點為高級節(jié)點（所含能量為普通節(jié)點的a倍）。MATLAB中的實驗結果圖如圖4.2所示。圖 4.2 LEA

18、CH和SEP協(xié)議節(jié)點存活數(shù)量圖Figure 4.2 The amount of alive node of LEACH protocols and SEP protocols由圖4.2可以看出LEACH協(xié)議第一節(jié)點死亡比SEP節(jié)點快，也就是SEP協(xié)議的穩(wěn)定期更長。這是由于SEP協(xié)議采用加權的簇頭選擇方法，初始能量高的高級節(jié)點成為簇頭的概率更高，從而死亡的時間將會更久。SEP協(xié)議相對于LEACH協(xié)議穩(wěn)定期延長了，從而吞吐量也增加了。5 結語本文首先對目前傳感器網(wǎng)絡路由方面的研究進行了分析，并對常見的路由算法進行分類，詳細地描述了現(xiàn)有的算法，并闡述了它們各自的優(yōu)缺點。針對經(jīng)典的層次路由算法LEAC

19、H進行了分析研究針對傳感器第一節(jié)點死亡時間等特點，提出了改進型SEP算法，SEP算法通過給節(jié)點設定不同的初始能量，然后加權選擇簇頭，剩余能量高的成為簇頭的概率就大的方法成功的延長了第一節(jié)點的死亡時間也就是延長了穩(wěn)定期，而且吞吐量也更大。參考文獻1 于洪斌.曾鵬.梁偉智能無線傳感器網(wǎng)絡系統(tǒng)科學出版社2 史永彬.葉湘濱.劉培亮無線傳感器網(wǎng)絡技術研究現(xiàn)狀國外電子測量技術, 2005,(11)3 宮曉淵.周興社.李志剛.劉剛無線傳感器網(wǎng)絡組織結構研究 微電子學與計算機 4 李方敏.劉新華.曠海蘭無線傳感器網(wǎng)絡中一種高能效低延時的泛洪算法研究通信學報 2007年,第08期5I.Akyildiz,W.Su

20、,Y. Sankarasubramaniam, and E. Cayirci, “A surveyon sensor networks,” IEEE Communications Magazine, vol. 40, no. 8,pp. 102114, August 2002.6 T. J. Shepard, “A channel access scheme for large dense packet radio networks,” in Proccedings of ACM SIGCOMM, September 1996, pp.219230.7 黃利. LEACH協(xié)議的研究及改進 D.

21、武漢：武漢理工大學碩士學位論文，2011.8 Zhang W, Zhao Q H , etc. Dynamic Perfo rmance Analysis on Discrete PSO Alg or ithm with Parameter Selection J . Journal of Sy st em Simulation. 2010, 22 ( 8) : 1899- 1904.9 孫利民,李建中,陳渝,朱紅松.無線傳感器網(wǎng)絡.北京:清華大學出版社,2005. 10 唐勇,周明天,張欣.無線傳感器網(wǎng)絡路由協(xié)議研究進展.軟件學報,2006,17(3):410421.附錄clear;xm=1

22、00;%設置區(qū)域為100*100ym=100;sink.x=0.5*xm;%sink（匯聚）節(jié)點坐標sink.y=0.5*ym;n=100 %區(qū)域內的節(jié)點數(shù)目p=0.1;% 節(jié)點成為簇頭的概率 Eo=0.04;%節(jié)點初始能量ETX=50*0.;%發(fā)射單位報文損耗能量 ERX=50*0.;%接收單位報文損耗能量Efs=10*0.1;%自由空間能量Emp=0.0013*0.1;%衰減空間能量EDA=5*0.;%聚集數(shù)據(jù)所要消耗的能量m=0.1;%成為高級節(jié)點比率a=1;%參數(shù)rmax=100%最大的輪數(shù)do=sqrt(Efs/Emp); %計算do 通信半徑。能量損耗的界限值：大于它則符合fris

23、s free space model，小于它則符合two-ray ground modelfigure(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;%？簇頭數(shù)目？ S(i).type=N;%節(jié)點類型為普通 S(i).E=Eo;%設置初始能量為E0 S(i).ENERGY=0;%普通節(jié)點標志 plot(S(i).xd,S(i).yd,o);%輸出節(jié)點，用o表示 h

24、old on;endS(n+1).xd=sink.x;%匯聚節(jié)點X軸坐標S(n+1).yd=sink.y;%匯聚節(jié)點Y軸坐標plot(S(n+1).xd,S(n+1).yd,x); %輸出匯聚節(jié)點，用x表示 %第一次迭代figure(1);cluster=1;flag_first_dead=0;%第一個節(jié)點死亡的標志變量for r=0:1:rmaxrif(mod(r, round(1/p) )=0)%如果所有節(jié)點都當過簇頭，則全部節(jié)點清零，回到最初狀態(tài) for i=1:1:n S(i).G=0;%簇頭數(shù)目 endenddead=0;figure(4);for i=1:1:n if (S(i).

25、E0)%節(jié)點能量大于0 S(i).type=N; plot(S(i).xd,S(i).yd,o); hold on; endendplot(S(n+1).xd,S(n+1).yd,x);%sinkSTATISTICS_leach(r+1).DEAD=dead;%r輪后死亡節(jié)點數(shù)DEAD_leach(r+1)=dead;%r輪后死亡節(jié)點數(shù)if (dead=1)%第一個節(jié)點死亡 if(flag_first_dead=0)%第一個節(jié)點死亡周期 first_dead=r%第一個節(jié)點死亡輪數(shù) flag_first_dead=1;%第一個死亡節(jié)點標志 endendcountCHs=0;%簇頭的個數(shù)clus

26、ter=1;%簇頭的標號，初始值為1for i=1:1:n%i為矩陣1到n，間距為1 if(S(i).E0)%節(jié)點剩余能量大于0 temp_rand=rand; if ( (S(i).G)=0)%沒有當過簇頭？ if(temp_randdo)%距離大于通信半徑 S(i).E=S(i).E- ( (ETX+EDA)*(4000) + Emp*4000*( distance*distance*distance*distance ); %能量消耗 end if (distance0 )%處理普通節(jié)點 if(cluster-1=1)%簇頭總數(shù)大于2個 min_dis=sqrt( (S(i).xd-S(

27、n+1).xd)2 + (S(i).yd-S(n+1).yd)2 );%普通節(jié)點到sink節(jié)點間最短距離 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 ) );%選取節(jié)點到簇頭距離和到sink節(jié)點距離之間最近的一個 if ( tempdo) S(i).E=S(i).E- ( ETX*(4000) + Emp*4000*( min_dis * min_dis * min_dis * min_dis); end i

28、f (min_dis0)%能量消散 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 endendSTATISTICS(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輪后節(jié)點剩余能量加上節(jié)點i的剩余能量 end

29、endhold off;endfor i=2:1:rmax%當前節(jié)點數(shù)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);%輸出圖形2hold on;%保持曲線plot(mylive,color,r);%用紅色輸出存活節(jié)點數(shù)xlabel(周期數(shù));ylabel(存活節(jié)點);title(存活節(jié)點圖);figure(3);%輸出圖形3hold on;%保持曲線plot(myener

30、gy,color,r);%用紅色輸出剩余能量xlabel(周期數(shù));ylabel(剩余能量節(jié)點);title(剩余能量圖);figure(6); for j=1:1:n S(j).xd=rand(1,1)*xm; XR(j)=S(j).xd; S(j).yd=rand(1,1)*ym; YR(j)=S(j).yd; S(j).G=0; S(j).type=N; temp_rnd0=j; if (temp_rnd0=m*n+1) %普通節(jié)點 S(j).E=Eo; S(j).ENERGY=0; hold on; end if (temp_rnd0m*n+1) %高級節(jié)點 S(j).E=Eo*(1+

31、a) S(j).ENERGY=1; hold on; endend S(n+1).xd=sink.x;S(n+1).yd=sink.y;figure(6);countCHs=0;rcountCHs=0;cluster=1; countCHs;rcountCHs=rcountCHs+countCHs;flag_first_dead=0; for r=0:1:rmax r pnrm=( p/ (1+a*m) ); %普通節(jié)點的選舉 padv= ( p*(1+a)/(1+a*m) ); %高級節(jié)點選舉 if(mod(r, round(1/pnrm) )=0) for j=1:1:n S(j).G=0

32、; S(j).cl=0; end end if(mod(r, round(1/padv) )=0) for j=1:1:n if(S(j).ENERGY=1) S(j).G=0; S(j).cl=0; end end endhold off; dead=0;%節(jié)點死亡數(shù)dead_a=0;%高級節(jié)點死亡數(shù)dead_n=0;%普通節(jié)點死亡數(shù)packets_TO_BS=0;packets_TO_CH=0;PACKETS_TO_CH(r+1)=0;PACKETS_TO_BS(r+1)=0; figure(6);for j=1:1:n if (S(j).E0 S(j).type=N; if (S(j).

33、ENERGY=0) plot(S(j).xd,S(j).yd,o); end if (S(j).ENERGY=1) plot(S(j).xd,S(j).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; if (dead=1) if(flag_first_dead=0) first_dead=r flag_first_dead=1; endend countCHs=0;cluster=1;for j=1:1:n if(S(j).E0) temp_rand=rand; if ( (S(j).G)=0) if( (