??本文模擬仿真TL-EBC協(xié)議，將LEACH算法與TL-EBC對(duì)比。

一、算法簡介

1、PSO-C(2007)

??LEACH算法選簇頭是分布式隨機(jī)選取的，存在很多缺陷，集中式選取簇頭，一般是將能量不小于網(wǎng)絡(luò)平均能量的節(jié)點(diǎn)作為候選簇頭，然后再從這些候選簇頭中選出適應(yīng)度比較好的簇頭，這是一個(gè)NP問題，一般可采用元啟發(fā)算法解決。
??假設(shè)網(wǎng)絡(luò)包含N個(gè)節(jié)點(diǎn)，預(yù)先定義分為K個(gè)簇頭，候選簇頭數(shù)為M（一般情況M>>K），則可能的分簇方式有C(M,K)種，在其中確定最佳分簇方式，是一個(gè)最優(yōu)化問題。應(yīng)用PSO算法解決這個(gè)問題，使每個(gè)粒子代表一種可能的分簇方式，用目標(biāo)函數(shù)評(píng)價(jià)其性能，設(shè)置m個(gè)粒子組成群體在C(M,K)種可能的分簇方式種尋找最優(yōu)解，使目標(biāo)函數(shù)取得最小值，該目標(biāo)函數(shù)定義如下：

2、TL-EBC(2012)

??TL-EBC協(xié)議在PSO-C協(xié)議的基礎(chǔ)上加入了一個(gè)總簇頭的機(jī)制，將網(wǎng)絡(luò)設(shè)計(jì)成了兩層分簇的體系結(jié)構(gòu)。
??總簇頭的選擇綜合考慮節(jié)點(diǎn)能量、各簇頭節(jié)點(diǎn)距離和 與基站距離 3 個(gè)方面的因素。首先，在每個(gè)簇中選擇能量 最大的節(jié)點(diǎn)(不包括簇頭)為候選總簇頭。然后，計(jì)算每個(gè)候選總簇頭的評(píng)價(jià)函數(shù)，具有最小評(píng)價(jià)函數(shù)值的候選總簇 頭成為總簇頭。評(píng)價(jià)函數(shù)定義如下：

二、仿真與分析

1、參數(shù)設(shè)置

??將200個(gè)傳感器節(jié)點(diǎn)隨機(jī)部署在200m x 200m的范圍內(nèi)，基站坐標(biāo)(100,275)，位于網(wǎng)絡(luò)外部。數(shù)據(jù)包長度為2000bit；控制包長度50bit；簇頭的數(shù)量設(shè)定為節(jié)點(diǎn)總數(shù)的5%，則K=10；總簇頭的數(shù)據(jù)不進(jìn)行壓縮。PSO算法的參數(shù)設(shè)置為：種群規(guī)模Q=30，c1=c2=2,慣性權(quán)重w = 1（這里的設(shè)計(jì)不同于論文，本人實(shí)驗(yàn)過，在該模型種，固定w=1，比線性變化和非線性變化的w設(shè)定更加好，不論是在網(wǎng)絡(luò)循環(huán)中，還是在PSO循環(huán)中），最大循環(huán)次數(shù)MaxIter = 50（論文種的MaxIter = 30是不夠的）。評(píng)價(jià)因子權(quán)重系數(shù)$\beta$ = 0.5,$\alpha$ = 0.2;$\eta$ = 0.4;$\lambda$ = 0.4。網(wǎng)絡(luò)只進(jìn)行同構(gòu)網(wǎng)絡(luò)的分析。

2、MATLAB 程序?qū)崿F(xiàn)

實(shí)驗(yàn)1、不考慮控制包，對(duì)比LEACH與TL-EBC

%不考慮控制包
%% 清空環(huán)境變量
clear;clc;
close all;%% 1、初始參數(shù)設(shè)定模塊
%監(jiān)測(cè)區(qū)域大小（單位：米）
xm = 200;
ym = 200;%基站的x和y坐標(biāo)
sink.x = 0.5 * xm;
sink.y = 275;%場(chǎng)地中的節(jié)點(diǎn)數(shù)
n = 200;%簇頭概率
p=0.05;%最大循環(huán)輪數(shù)
rmax = 500;%能量模型(所有值均以焦耳為單位)
Eo = 0.1;%初始能量
ETX = 50*0.000000001;%傳輸能耗
ERX = 50*0.000000001;%接收能耗
EDA = 5*0.000000001;%數(shù)據(jù)聚合能量%發(fā)射放大器類型
Efs = 10*0.000000000001;%小于do
Emp = 0.0013*0.000000000001;%大于do
do = sqrt(Efs/Emp);%距離閾值%數(shù)據(jù)包和控制包長度
packetLength = 2000;%數(shù)據(jù)包長度 
ctrPacketLength = 50;%控制包長度%% 2、無線傳感器網(wǎng)絡(luò)模型產(chǎn)生模塊
%構(gòu)建無線傳感器網(wǎng)絡(luò),在區(qū)域內(nèi)均勻投放所有節(jié)點(diǎn),并畫出圖形
figure;
for i = 1:nS(i).xd = rand(1,1)*xm;%坐標(biāo)S(i).yd = rand(1,1)*ym;S(i).E = Eo;%初始化節(jié)點(diǎn)能量%最初沒有簇頭只有簇成員S(i).type = 'N';%簇成員(普通節(jié)點(diǎn)) plot(S(i).xd,S(i).yd,'o');hold on;N(i).xd = S(i).xd;%坐標(biāo)N(i).yd = S(i).yd;N(i).G = 0;%"G = 0"表示節(jié)點(diǎn)本周期內(nèi)未當(dāng)選過簇頭，當(dāng)選則改為"G = 1"N(i).E = Eo;%初始化節(jié)點(diǎn)能量%最初沒有簇頭只有簇成員N(i).type = 'N';%簇成員(普通節(jié)點(diǎn)) 
end
%加入基站信息
S(n+1).xd = sink.x;
S(n+1).yd = sink.y;
N(n+1).xd = sink.x;
N(n+1).yd = sink.y;
plot(S(n+1).xd,S(n+1).yd,'p');%% 3、PSO參數(shù)
MaxIter = 50;      % 最大迭代次數(shù)
nPop = 30;        % 種群規(guī)模
wmax = 0.9;
wmin = 0.4;
w = 1;
c1 = 2.0;         % 個(gè)人學(xué)習(xí)系數(shù)
c2 = 2.0;         % 全局學(xué)習(xí)系數(shù)
plimit = [0, 200];                % 位置參數(shù)限制
vlimit = [-20, 20];               % 速度限制%數(shù)據(jù)初始化
particle.Position = cell(1,nPop);
particle.Cost = zeros(1,nPop);
particle.Velocity = cell(1,nPop);
particle.Best.Position = cell(1,nPop);
particle.Best.Cost = ones(1,nPop)*inf;
GlobalBestCost = inf;
GlobalBestPosition = [];%% 4、網(wǎng)絡(luò)運(yùn)行模塊 
% 死亡節(jié)點(diǎn)數(shù)指標(biāo)
first_dead = 0;%第一個(gè)死亡節(jié)點(diǎn)
teenth_dead = 0;%10%的死亡節(jié)點(diǎn)
all_dead = 0;%節(jié)點(diǎn)都死亡%節(jié)點(diǎn)死亡數(shù)指標(biāo)對(duì)應(yīng)標(biāo)記
flag_first_dead = 0;
flag_teenth_dead = 0;
flag_all_dead = 0;%傳輸?shù)交竞痛仡^的比特計(jì)數(shù)器
packets_TO_BS = 0;
packets_TO_CH = 0;%主循環(huán)
for r = 0:rmax dead = 0;%計(jì)算死亡節(jié)點(diǎn)數(shù)energy = 0;%每輪節(jié)點(diǎn)剩余總能量Node = [];%存活節(jié)點(diǎn)集合C = [];%候選節(jié)點(diǎn)集合j = 0;for i = 1:n%檢查有無死亡節(jié)點(diǎn)if S(i).E <= 0dead = dead+1;%第一個(gè)節(jié)點(diǎn)死亡時(shí)間if dead == 1if flag_first_dead == 0first_dead = r+1;flag_first_dead = 1;endend% 10%的節(jié)點(diǎn)死亡時(shí)間if dead == 0.1*nif flag_teenth_dead ==0teenth_dead = r+1;flag_teenth_dead = 1;endendif dead == nif flag_all_dead == 0all_dead = r+1;flag_all_dead = 1;endendendif S(i).E > 0energy = energy + S(i).E;%統(tǒng)計(jì)每輪節(jié)點(diǎn)總能量和S(i).type = 'N';j = j+1;Node(j,:) = [S(i).xd S(i).yd];%存活節(jié)點(diǎn)集合endendSTATISTICS.DEAD(r+1) = dead;%每輪死亡節(jié)點(diǎn)數(shù)STATISTICS.alive(r+1) = n-dead;%每輪存活節(jié)點(diǎn)數(shù)STATISTICS.energy(r+1) = energy;%每輪剩余總能量if n-dead == 0STATISTICS.COUNTCHS(r+1) = 0;continue;end%% 4.1、應(yīng)用PSO選簇頭%候選簇頭集合j = 0;for i = 1:nif r ==0 %因?yàn)閙atlab最大精度16位，后面數(shù)據(jù)不準(zhǔn)，第一輪會(huì)出現(xiàn)這種情況，所以這樣編程j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候選簇頭坐標(biāo)C.E(j) = S(i).E;%候選簇頭剩余能量C.id(j) = i;%候選簇頭對(duì)應(yīng)的節(jié)點(diǎn)編號(hào)elseif S(i).E >= energy/(n-dead) %節(jié)點(diǎn)能量大于或等于總的平均能量j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候選簇頭坐標(biāo)C.E(j) = S(i).E;%候選簇頭剩余能量C.id(j) = i;%候選簇頭對(duì)應(yīng)的節(jié)點(diǎn)編號(hào)endendendM = size(C.d,1);%候選簇頭數(shù)K = ceil((n-dead)*p);% 簇頭數(shù)%初始化種群位置和速度for i = 1:nPop% 初始化位置C_k = randperm(M,K);particle.Position{i} = C.d(C_k,:);% 初始化速度v = zeros(K,2);particle.Velocity{i}=v;% 適應(yīng)度energe_Q = sum(C.E((C_k)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%傳感器簇成員節(jié)點(diǎn)矩陣particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = 1;%計(jì)數(shù)while Iter <= MaxIter%更新粒子個(gè)體最佳適應(yīng)度、位置   for i = 1:nPop      if particle.Cost(i) < particle.Best.Cost(i) particle.Best.Cost(i) = particle.Cost(i);%更新個(gè)體歷史最佳適應(yīng)度particle.Best.Position{i} = particle.Position{i};%更新個(gè)體歷史最佳位置end end%更新群體最佳適應(yīng)度、位置if min(particle.Best.Cost) < GlobalBestCost  [GlobalBestCost, index1] = min(particle.Best.Cost);%更新群體歷史最佳適應(yīng)度GlobalBestPosition = particle.Best.Position{index1};%更新群體歷史最佳位置end% 速度更新for i = 1:nPopparticle.Velocity{i} = particle.Velocity{i}.*w + c1*rand*(particle.Best.Position{i}-particle.Position{i}) ...+ c2*rand*(GlobalBestPosition-particle.Position{i});% 速度更新end% 邊界速度處理for i = 1:nPopparticle.Velocity{i}(particle.Velocity{i}>20) = 20;particle.Velocity{i}(particle.Velocity{i}<-20) = -20;end% 位置更新for i = 1:nPopparticle.Position{i} = particle.Position{i} + particle.Velocity{i};end% 邊界位置處理for i = 1:nPopparticle.Position{i}(particle.Position{i}>200)=200;particle.Position{i}(particle.Position{i}<0)=0;end%根據(jù)距離最近候選簇頭位置調(diào)整粒子位置,計(jì)算相應(yīng)適應(yīng)度for i = 1:nPopdistance = pdist2(C.d,particle.Position{i});[~, index2] = min(distance);particle.Position{i} = C.d(index2,:);% 適應(yīng)度energe_Q = sum(C.E((index2)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%傳感器簇成員節(jié)點(diǎn)矩陣particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = Iter+1;end%根據(jù)距離最近候選簇頭位置得到簇頭信息distance = pdist2(C.d,GlobalBestPosition);[~, index2] = min(distance);%如果簇頭重復(fù)，需要調(diào)整(這一步處理的不夠好)length_index2 = length(unique(index2));cha = K-length_index2;%重復(fù)個(gè)數(shù)while cha ~= 0for i = 1:K-1for j = i+1:Kif index2(i) == index2(j)distance_cut = distance(:,j);%取distance的第j列[~,Number] = sort(distance_cut);Num = 1;while length(unique(index2))~=K-(cha-1)  && Num<size(distance_cut,1)Num = Num+1;index2(j) = Number(Num);endcha = cha-1;endendendendfor i = 1:K%簇頭接受基站的廣播分簇信息S(C.id(index2(i))).E =  S(C.id(index2(i))).E - ERX * ctrPacketLength;end%重置群體和個(gè)體最佳適應(yīng)度、位置GlobalBestCost = inf;GlobalBestPosition = [];particle.Best.Cost = ones(1,nPop)*inf;particle.Best.Position = cell(1,nPop);%% 4.2選總簇頭%確定簇頭for i = 1:Kj = C.id(index2(i));S(j).type = 'C';endEach_cluster = zeros(K,2);%第一列放每個(gè)簇的能量最大的節(jié)點(diǎn)的能量信息，第二列放節(jié)點(diǎn)標(biāo)號(hào),行對(duì)應(yīng)簇號(hào)%預(yù)成簇for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通節(jié)點(diǎn)%加入最近的簇頭distance_to_cluster = pdist2(C.d(index2,:),[S(i).xd S(i).yd]);[min_dis, dis_index] = min(distance_to_cluster);%min_dis返回最短距離，dis_index返回最短距離對(duì)應(yīng)的簇號(hào)if S(i).E > Each_cluster(dis_index,1)Each_cluster(dis_index,1) = S(i).E;Each_cluster(dis_index,2) = i;end endend%若有空簇，將其候選總簇頭節(jié)點(diǎn)假定為1for i = 1:Kif Each_cluster(i,2) == 0Each_cluster(i,2) = 1;endend%每個(gè)候選總簇頭至所有簇頭距離和sum_distance_HCH_to_CH = sum(pdist2(C.d(index2,:),reshape([S(Each_cluster(:,2)).xd S(Each_cluster(:,2)).yd],[],2)));%每個(gè)候選總簇頭至基站距離distance_HCH_to_BS = pdist2([S(n+1).xd S(n+1).yd],C.d(index2,:));CostofHCH = zeros(1,K);for i = 1:KCostofHCH(i) = Cost_HCH(Each_cluster(:,1),sum_distance_HCH_to_CH,distance_HCH_to_BS,i);end[~,index3] = min(CostofHCH);S(Each_cluster(index3,2)).type = 'T';%總簇頭S(Each_cluster(index3,2)).E = S(Each_cluster(index3,2)).E - ERX * ctrPacketLength;%總簇頭接收基站發(fā)來的分簇信息
%% 4.3成簇，簇內(nèi)通信for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通節(jié)點(diǎn)%加入最近的簇頭distance_to_cluster = pdist2(C.d(index2,:),[S(i).xd S(i).yd]);[min_dis, dis_index] = min(distance_to_cluster);%min_dis返回最短距離，min_dis_cluster返回最短距離對(duì)應(yīng)的簇號(hào)%接收基站發(fā)來的廣播的消耗（收到基站發(fā)來的分簇信息）S(i).E = S(i).E - ERX * ctrPacketLength;%向簇頭發(fā)送數(shù)據(jù)包if (min_dis > do)S(i).E = S(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇發(fā)送頭數(shù)據(jù)包elseS(i).E = S(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇頭發(fā)送數(shù)據(jù)包end%簇頭接收簇成員的數(shù)據(jù)包if(min_dis > 0)S(C.id(dis_index)).E = S(C.id(dis_index)).E - ( (ERX + EDA)*packetLength ); %接受簇成員發(fā)來的數(shù)據(jù)包endendend
%% 4.4簇頭和總簇頭通信TCnumber = Each_cluster(index3,2);%總簇頭編號(hào)for i = 1:Kj = C.id(index2(i));S(j).type = 'C';%簇頭到總簇頭的距離distance = sqrt((S(j).xd-S(TCnumber).xd)^2 + (S(j).yd-(S(TCnumber).yd))^2 );%簇頭向總簇頭發(fā)送數(shù)據(jù)包能耗if (distance > do)S(j).E = S(j).E- ( ETX * (packetLength) + Emp * packetLength*distance^4);elseS(j).E = S(j).E- ( ETX * (packetLength) + Efs * packetLength*distance^2); end%總簇頭接收簇頭的的數(shù)據(jù)包(不融合數(shù)據(jù))if(distance > 0)S(TCnumber).E = S(TCnumber).E - ERX * packetLength; %接受簇頭發(fā)來的數(shù)據(jù)包endend
%% 4.5總簇頭和基站通信distance_to_BS = pdist2([S(n+1).xd S(n+1).yd],[S(TCnumber).xd S(TCnumber).yd]);if (distance_to_BS > do)S(TCnumber).E = S(TCnumber).E- (K * ETX * (packetLength) + Emp * packetLength*distance^4);elseS(TCnumber).E = S(TCnumber).E- (K * ETX * (packetLength) + Efs * packetLength*distance^2); endSTATISTICS.COUNTCHS(r+1) = K;
end%% 4、網(wǎng)絡(luò)運(yùn)行模塊(Leach)
% 死亡節(jié)點(diǎn)數(shù)指標(biāo)
first_dead1 = 0;%第一個(gè)死亡節(jié)點(diǎn)
teenth_dead1 = 0;%10%的死亡節(jié)點(diǎn)
all_dead1 = 0;%節(jié)點(diǎn)都死亡%節(jié)點(diǎn)死亡數(shù)指標(biāo)對(duì)應(yīng)標(biāo)記
flag_first_dead1 = 0;
flag_teenth_dead1 = 0;
flag_all_dead1 = 0;%傳輸?shù)交竞痛仡^的比特計(jì)數(shù)器
packets_TO_BS1 = 0;
packets_TO_CH1 = 0;%主循環(huán)
for r = 0:rmax %每過一個(gè)周期，重置Gif(mod(r, round(1/p)) == 0)for i = 1:nN(i).G = 0;endenddead1 = 0;%計(jì)算死亡節(jié)點(diǎn)數(shù)energy1 = 0;%每輪節(jié)點(diǎn)剩余總能量for i = 1:n%檢查有無死亡節(jié)點(diǎn)if N(i).E <= 0dead1 = dead1+1;%第一個(gè)節(jié)點(diǎn)死亡時(shí)間if dead1 == 1if flag_first_dead1 == 0first_dead1 = r+1;flag_first_dead1 = 1;endend% 10%的節(jié)點(diǎn)死亡時(shí)間if dead1 == 0.1*nif flag_teenth_dead1 ==0teenth_dead1 = r+1;flag_teenth_dead1 = 1;endendif dead1 == nif flag_all_dead1 == 0all_dead1 = r+1;flag_all_dead1 = 1;endendendif N(i).E > 0energy1 = energy1 + N(i).E;%統(tǒng)計(jì)每輪節(jié)點(diǎn)總能量和N(i).type = 'N';endend%節(jié)點(diǎn)死亡，退出循環(huán)，但是一般都是與其他算法比較，由于不同算法最終循環(huán)倫次不同%所以節(jié)點(diǎn)死亡并不退出循環(huán)，而是在達(dá)到最大循環(huán)次數(shù)后才結(jié)束循環(huán)
%     if flag_all_dead == 1
%         break;
%     end%r從0開始，但是計(jì)算機(jī)中數(shù)組從1開始，所以采用r+1的下標(biāo)STATISTICS1.DEAD(r+1) = dead1;%每輪死亡節(jié)點(diǎn)數(shù)STATISTICS1.alive(r+1) = n-dead1;%每輪存活節(jié)點(diǎn)數(shù)STATISTICS1.energy(r+1) = energy1;%每輪剩余總能量cluster=0;%計(jì)數(shù)簇頭數(shù)%次循環(huán)（選簇頭）if energy1 > 0%網(wǎng)絡(luò)還剩余能量，這一步在節(jié)點(diǎn)都死亡后才有用for i=1:nif(N(i).E>0)%節(jié)點(diǎn)必須有能量temp_rand=rand;     if ( (N(i).G)<=0) %該節(jié)點(diǎn)在候選集合中%選簇頭if( temp_rand <= (p/(1-p*mod(r,round(1/p)))))cluster = cluster+1;N(i).type = 'C'; %簇頭N(i).G = round(1/p)-1;%換成"S(i).G = 1"也行，只要不是0C(cluster).xd = N(i).xd;C(cluster).yd = N(i).yd;distance=sqrt( (N(i).xd-(N(n+1).xd))^2 + (N(i).yd-(N(n+1).yd))^2 );%距離基站的距離C(cluster).distance = distance;C(cluster).id = i;distanceBroad = sqrt(xm^2+ym^2);%全網(wǎng)廣播自成為簇頭,廣播TDMA時(shí)間表，距離是distanceBroadif (distanceBroad > do)N(i).E = N(i).E- (2* ETX * ctrPacketLength + Emp * ctrPacketLength*distanceBroad^4);elseN(i).E = N(i).E- (2* ETX * ctrPacketLength + Efs * ctrPacketLength*distanceBroad^2); end%簇頭向基站發(fā)送數(shù)據(jù)包能耗if (distance > do)N(i).E = N(i).E- ( ETX*(packetLength) + Emp * packetLength*distance^4);elseN(i).E = N(i).E- ( ETX*(packetLength) + Efs * packetLength*distance^2); endend     endend endendSTATISTICS1.COUNTCHS(r+1) = cluster;%次循環(huán)(成簇)if energy1 > 0%網(wǎng)絡(luò)還剩余能量，這一步在節(jié)點(diǎn)都死亡后才生效for i = 1:nif ( N(i).type=='N' && N(i).E>0 ) %普通節(jié)點(diǎn)%如果有簇頭存在if(cluster>= 1)length = zeros(cluster,1);%加入最近的簇頭for c = 1:clusterlength(c) = sqrt( (N(i).xd - C(c).xd)^2 + (N(i).yd - C(c).yd)^2 );%%接收簇頭發(fā)來的廣播消耗N(i).E = N(i).E - ERX * ctrPacketLength;end[min_dis, min_dis_cluster] = min(length);%min_dis返回最短距離，min_dis_cluster返回最短距離對(duì)應(yīng)的簇號(hào)%向簇頭發(fā)送數(shù)據(jù)包if (min_dis > do)N(i).E = N(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇發(fā)送頭數(shù)據(jù)包elseN(i).E = N(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇頭發(fā)送數(shù)據(jù)包end%簇頭接收簇成員的數(shù)據(jù)包if(min_dis > 0)N(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ( (ERX + EDA)*packetLength ); %接受簇成員發(fā)來的數(shù)據(jù)包end%如果本輪沒有簇頭，則普通節(jié)點(diǎn)直接將數(shù)據(jù)直接發(fā)送給基站elsemin_dis = sqrt((N(i).xd-N(n+1).xd)^2 + (N(i).yd-N(n+1).yd)^2);endendendend
end
%% 繪圖顯示
figure;
plot(0:rmax, STATISTICS.DEAD, 'r',0:rmax, STATISTICS1.DEAD, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('死亡節(jié)點(diǎn)數(shù)');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.alive, 'r',0:rmax, STATISTICS1.alive, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('活動(dòng)節(jié)點(diǎn)數(shù)');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.energy, 'r',0:rmax, STATISTICS1.energy, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('剩余能量');
legend('mypso','Leach')  
figure;
plot(0:rmax, STATISTICS.COUNTCHS, 'r',0:rmax, STATISTICS1.COUNTCHS, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('簇頭數(shù)');
legend('mypso','Leach')

實(shí)驗(yàn)2、慮控制包，對(duì)比LEACH與TL-EBC

%% 清空環(huán)境變量
clear;clc;
close all;%% 1、初始參數(shù)設(shè)定模塊
%監(jiān)測(cè)區(qū)域大小（單位：米）
xm = 200;
ym = 200;%基站的x和y坐標(biāo)
sink.x = 0.5 * xm;
sink.y = 275;%場(chǎng)地中的節(jié)點(diǎn)數(shù)
n = 200;%簇頭概率
p=0.05;%最大循環(huán)輪數(shù)
rmax = 500;%能量模型(所有值均以焦耳為單位)
Eo = 0.1;%初始能量
ETX = 50*0.000000001;%傳輸能耗
ERX = 50*0.000000001;%接收能耗
EDA = 5*0.000000001;%數(shù)據(jù)聚合能量%發(fā)射放大器類型
Efs = 10*0.000000000001;%小于do
Emp = 0.0013*0.000000000001;%大于do
do = sqrt(Efs/Emp);%距離閾值%數(shù)據(jù)包和控制包長度
packetLength = 2000;%數(shù)據(jù)包長度 
ctrPacketLength = 50;%控制包長度%% 2、無線傳感器網(wǎng)絡(luò)模型產(chǎn)生模塊
%構(gòu)建無線傳感器網(wǎng)絡(luò),在區(qū)域內(nèi)均勻投放所有節(jié)點(diǎn),并畫出圖形
figure;
for i = 1:nS(i).xd = rand(1,1)*xm;%坐標(biāo)S(i).yd = rand(1,1)*ym;S(i).E = Eo;%初始化節(jié)點(diǎn)能量%最初沒有簇頭只有簇成員S(i).type = 'N';%簇成員(普通節(jié)點(diǎn)) plot(S(i).xd,S(i).yd,'o');hold on;N(i).xd = S(i).xd;%坐標(biāo)N(i).yd = S(i).yd;N(i).G = 0;%"G = 0"表示節(jié)點(diǎn)本周期內(nèi)未當(dāng)選過簇頭，當(dāng)選則改為"G = 1"N(i).E = Eo;%初始化節(jié)點(diǎn)能量%最初沒有簇頭只有簇成員N(i).type = 'N';%簇成員(普通節(jié)點(diǎn)) 
end
%加入基站信息
S(n+1).xd = sink.x;
S(n+1).yd = sink.y;
N(n+1).xd = sink.x;
N(n+1).yd = sink.y;
plot(S(n+1).xd,S(n+1).yd,'p');%% 3、PSO參數(shù)
MaxIter = 50;      % 最大迭代次數(shù)
nPop = 30;        % 種群規(guī)模
wmax = 0.9;
wmin = 0.4;
w = 1;
c1 = 2.0;         % 個(gè)人學(xué)習(xí)系數(shù)
c2 = 2.0;         % 全局學(xué)習(xí)系數(shù)
plimit = [0, 200];                % 位置參數(shù)限制
vlimit = [-20, 20];               % 速度限制%數(shù)據(jù)初始化
particle.Position = cell(1,nPop);
particle.Cost = zeros(1,nPop);
particle.Velocity = cell(1,nPop);
particle.Best.Position = cell(1,nPop);
particle.Best.Cost = ones(1,nPop)*inf;
GlobalBestCost = inf;
GlobalBestPosition = [];%% 4、網(wǎng)絡(luò)運(yùn)行模塊 
% 死亡節(jié)點(diǎn)數(shù)指標(biāo)
first_dead = 0;%第一個(gè)死亡節(jié)點(diǎn)
teenth_dead = 0;%10%的死亡節(jié)點(diǎn)
all_dead = 0;%節(jié)點(diǎn)都死亡%節(jié)點(diǎn)死亡數(shù)指標(biāo)對(duì)應(yīng)標(biāo)記
flag_first_dead = 0;
flag_teenth_dead = 0;
flag_all_dead = 0;%傳輸?shù)交竞痛仡^的比特計(jì)數(shù)器
packets_TO_BS = 0;
packets_TO_CH = 0;%主循環(huán)
for r = 0:rmax dead = 0;%計(jì)算死亡節(jié)點(diǎn)數(shù)energy = 0;%每輪節(jié)點(diǎn)剩余總能量Node = [];%存活節(jié)點(diǎn)集合C = [];%候選節(jié)點(diǎn)集合j = 0;for i = 1:n%檢查有無死亡節(jié)點(diǎn)if S(i).E <= 0dead = dead+1;%第一個(gè)節(jié)點(diǎn)死亡時(shí)間if dead == 1if flag_first_dead == 0first_dead = r+1;flag_first_dead = 1;endend% 10%的節(jié)點(diǎn)死亡時(shí)間if dead == 0.1*nif flag_teenth_dead ==0teenth_dead = r+1;flag_teenth_dead = 1;endendif dead == nif flag_all_dead == 0all_dead = r+1;flag_all_dead = 1;endendendif S(i).E > 0energy = energy + S(i).E;%統(tǒng)計(jì)每輪節(jié)點(diǎn)總能量和S(i).type = 'N';j = j+1;Node(j,:) = [S(i).xd S(i).yd];%存活節(jié)點(diǎn)集合endendSTATISTICS.DEAD(r+1) = dead;%每輪死亡節(jié)點(diǎn)數(shù)STATISTICS.alive(r+1) = n-dead;%每輪存活節(jié)點(diǎn)數(shù)STATISTICS.energy(r+1) = energy;%每輪剩余總能量if n-dead == 0STATISTICS.COUNTCHS(r+1) = 0;continue;end%% 4.1、應(yīng)用PSO選簇頭%候選簇頭集合j = 0;for i = 1:nif r ==0 %因?yàn)閙atlab最大精度16位，后面數(shù)據(jù)不準(zhǔn)，第一輪會(huì)出現(xiàn)這種情況，所以這樣編程j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候選簇頭坐標(biāo)C.E(j) = S(i).E;%候選簇頭剩余能量C.id(j) = i;%候選簇頭對(duì)應(yīng)的節(jié)點(diǎn)編號(hào)elseif S(i).E >= energy/(n-dead) %節(jié)點(diǎn)能量大于或等于總的平均能量j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候選簇頭坐標(biāo)C.E(j) = S(i).E;%候選簇頭剩余能量C.id(j) = i;%候選簇頭對(duì)應(yīng)的節(jié)點(diǎn)編號(hào)endendendM = size(C.d,1);%候選簇頭數(shù)K = ceil((n-dead)*p);% 簇頭數(shù)
%     w = 0.9-0.5*(10-K);%慣性重量%初始化種群位置和速度for i = 1:nPop% 初始化位置C_k = randperm(M,K);particle.Position{i} = C.d(C_k,:);% 初始化速度v = zeros(K,2);particle.Velocity{i}=v;% 適應(yīng)度energe_Q = sum(C.E((C_k)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%傳感器簇成員節(jié)點(diǎn)矩陣particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = 1;%計(jì)數(shù)while Iter <= MaxIter%更新粒子個(gè)體最佳適應(yīng)度、位置   for i = 1:nPop      if particle.Cost(i) < particle.Best.Cost(i) particle.Best.Cost(i) = particle.Cost(i);%更新個(gè)體歷史最佳適應(yīng)度particle.Best.Position{i} = particle.Position{i};%更新個(gè)體歷史最佳位置end end%更新群體最佳適應(yīng)度、位置if min(particle.Best.Cost) < GlobalBestCost  [GlobalBestCost, index1] = min(particle.Best.Cost);%更新群體歷史最佳適應(yīng)度GlobalBestPosition = particle.Best.Position{index1};%更新群體歷史最佳位置end% 速度更新for i = 1:nPopparticle.Velocity{i} = particle.Velocity{i}.*w + c1*rand*(particle.Best.Position{i}-particle.Position{i}) ...+ c2*rand*(GlobalBestPosition-particle.Position{i});% 速度更新end% 邊界速度處理for i = 1:nPopparticle.Velocity{i}(particle.Velocity{i}>20) = 20;particle.Velocity{i}(particle.Velocity{i}<-20) = -20;end% 位置更新for i = 1:nPopparticle.Position{i} = particle.Position{i} + particle.Velocity{i};end% 邊界位置處理for i = 1:nPopparticle.Position{i}(particle.Position{i}>200)=200;particle.Position{i}(particle.Position{i}<0)=0;end%根據(jù)距離最近候選簇頭位置調(diào)整粒子位置,計(jì)算相應(yīng)適應(yīng)度for i = 1:nPopdistance = pdist2(C.d,particle.Position{i});[~, index2] = min(distance);particle.Position{i} = C.d(index2,:);% 適應(yīng)度energe_Q = sum(C.E((index2)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%傳感器簇成員節(jié)點(diǎn)矩陣particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = Iter+1;end%根據(jù)距離最近候選簇頭位置得到簇頭信息distance = pdist2(C.d,GlobalBestPosition);[~, index2] = min(distance);%如果簇頭重復(fù)，需要調(diào)整(這一步處理的不夠好)length_index2 = length(unique(index2));cha = K-length_index2;%重復(fù)個(gè)數(shù)while cha ~= 0for i = 1:K-1for j = i+1:Kif index2(i) == index2(j)distance_cut = distance(:,j);%取distance的第j列[~,Number] = sort(distance_cut);Num = 1;while length(unique(index2))~=K-(cha-1)  && Num<size(distance_cut,1)Num = Num+1;index2(j) = Number(Num);endcha = cha-1;endendendenddistanceBroad = sqrt(xm^2+ym^2)/3;%假定簇頭和簇成員最遠(yuǎn)距離是對(duì)角線的1/3(懶得去后面算這個(gè)實(shí)際距離了,在分簇均勻的情況下其實(shí)是對(duì)角線的1/6)                    for i = 1:K%簇頭接受基站的廣播分簇信息S(C.id(index2(i))).E =  S(C.id(index2(i))).E - ERX * ctrPacketLength;%簇頭廣播TDMA表if (distanceBroad > do)S(C.id(index2(i))).E = S(C.id(index2(i))).E - ( ETX*ctrPacketLength + Emp*packetLength*distanceBroad^4);elseS(C.id(index2(i))).E = S(C.id(index2(i))).E - ( ETX*ctrPacketLength + Efs*packetLength*distanceBroad^2);endend%重置群體和個(gè)體最佳適應(yīng)度、位置GlobalBestCost = inf;GlobalBestPosition = [];particle.Best.Cost = ones(1,nPop)*inf;particle.Best.Position = cell(1,nPop);%% 4.2選總簇頭%確定簇頭for i = 1:Kj = C.id(index2(i));S(j).type = 'C';endEach_cluster = zeros(K,2);%第一列放每個(gè)簇的能量最大的節(jié)點(diǎn)的能量信息，第二列放節(jié)點(diǎn)標(biāo)號(hào),行對(duì)應(yīng)簇號(hào)%預(yù)成簇for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通節(jié)點(diǎn)%加入最近的簇頭distance_to_cluster = pdist2(C.d(index2,:),[S(i).xd S(i).yd]);[min_dis, dis_index] = min(distance_to_cluster);%min_dis返回最短距離，dis_index返回最短距離對(duì)應(yīng)的簇號(hào)if S(i).E > Each_cluster(dis_index,1)Each_cluster(dis_index,1) = S(i).E;Each_cluster(dis_index,2) = i;end endend%若有空簇，將其候選總簇頭節(jié)點(diǎn)假定為1for i = 1:Kif Each_cluster(i,2) == 0Each_cluster(i,2) = 1;endend%每個(gè)候選總簇頭至所有簇頭距離和sum_distance_HCH_to_CH = sum(pdist2(C.d(index2,:),reshape([S(Each_cluster(:,2)).xd S(Each_cluster(:,2)).yd],[],2)));%每個(gè)候選總簇頭至基站距離distance_HCH_to_BS = pdist2([S(n+1).xd S(n+1).yd],C.d(index2,:));CostofHCH = zeros(1,K);for i = 1:KCostofHCH(i) = Cost_HCH(Each_cluster(:,1),sum_distance_HCH_to_CH,distance_HCH_to_BS,i);end[~,index3] = min(CostofHCH);S(Each_cluster(index3,2)).type = 'T';%總簇頭S(Each_cluster(index3,2)).E = S(Each_cluster(index3,2)).E - ERX * ctrPacketLength;%總簇頭接收基站發(fā)來的分簇信息
%% 4.3成簇，簇內(nèi)通信for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通節(jié)點(diǎn)%加入最近的簇頭distance_to_cluster = pdist2(C.d(index2,:),[S(i).xd S(i).yd]);[min_dis, dis_index] = min(distance_to_cluster);%min_dis返回最短距離，min_dis_cluster返回最短距離對(duì)應(yīng)的簇號(hào)%接收基站發(fā)來的廣播的消耗（收到基站發(fā)來的分簇信息）S(i).E = S(i).E - ERX * ctrPacketLength;%接收簇頭的廣播TDMA時(shí)間表信息S(i).E = S(i).E - ERX * ctrPacketLength;%向簇頭發(fā)送數(shù)據(jù)包if (min_dis > do)S(i).E = S(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇發(fā)送頭數(shù)據(jù)包elseS(i).E = S(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇頭發(fā)送數(shù)據(jù)包end%簇頭接收簇成員的數(shù)據(jù)包if(min_dis > 0)S(C.id(dis_index)).E = S(C.id(dis_index)).E - ( (ERX + EDA)*packetLength ); %接受簇成員發(fā)來的數(shù)據(jù)包endendend
%% 4.4簇頭和總簇頭通信TCnumber = Each_cluster(index3,2);%總簇頭編號(hào)for i = 1:Kj = C.id(index2(i));S(j).type = 'C';%簇頭到總簇頭的距離distance = sqrt((S(j).xd-S(TCnumber).xd)^2 + (S(j).yd-(S(TCnumber).yd))^2 );%接收總簇頭的廣播TDMA時(shí)間表信息S(j).E = S(i).E - ERX*ctrPacketLength;%簇頭向總簇頭發(fā)送數(shù)據(jù)包能耗if (distance > do)S(j).E = S(j).E- ( ETX * (packetLength) + Emp * packetLength*distance^4);elseS(j).E = S(j).E- ( ETX * (packetLength) + Efs * packetLength*distance^2); end%總簇頭接收簇頭的的數(shù)據(jù)包(不融合數(shù)據(jù))if(distance > 0)S(TCnumber).E = S(TCnumber).E - ERX * packetLength; %接受簇頭發(fā)來的數(shù)據(jù)包endend
%% 4.5總簇頭和基站通信%總簇頭廣播時(shí)隙表,距離是maxdistancemaxdistance = max(pdist2(C.d(index2,:),[S(TCnumber).xd S(TCnumber).yd]));%距離最遠(yuǎn)的簇頭的距離if (maxdistance > do)S(TCnumber).E = S(TCnumber).E- ( ETX * ctrPacketLength + Emp * ctrPacketLength*maxdistance^4);elseS(TCnumber).E = S(TCnumber).E- ( ETX * ctrPacketLength + Efs * ctrPacketLength*maxdistance^2); enddistance_to_BS = pdist2([S(n+1).xd S(n+1).yd],[S(TCnumber).xd S(TCnumber).yd]);if (distance_to_BS > do)S(TCnumber).E = S(TCnumber).E- (K * ETX * (packetLength) + Emp * packetLength*distance^4);elseS(TCnumber).E = S(TCnumber).E- (K * ETX * (packetLength) + Efs * packetLength*distance^2); endSTATISTICS.COUNTCHS(r+1) = K;
end%% 4、網(wǎng)絡(luò)運(yùn)行模塊(Leach)
% 死亡節(jié)點(diǎn)數(shù)指標(biāo)
first_dead1 = 0;%第一個(gè)死亡節(jié)點(diǎn)
teenth_dead1 = 0;%10%的死亡節(jié)點(diǎn)
all_dead1 = 0;%節(jié)點(diǎn)都死亡%節(jié)點(diǎn)死亡數(shù)指標(biāo)對(duì)應(yīng)標(biāo)記
flag_first_dead1 = 0;
flag_teenth_dead1 = 0;
flag_all_dead1 = 0;%傳輸?shù)交竞痛仡^的比特計(jì)數(shù)器
packets_TO_BS1 = 0;
packets_TO_CH1 = 0;%主循環(huán)
for r = 0:rmax %每過一個(gè)周期，重置Gif(mod(r, round(1/p)) == 0)for i = 1:nN(i).G = 0;endenddead1 = 0;%計(jì)算死亡節(jié)點(diǎn)數(shù)energy1 = 0;%每輪節(jié)點(diǎn)剩余總能量for i = 1:n%檢查有無死亡節(jié)點(diǎn)if N(i).E <= 0dead1 = dead1+1;%第一個(gè)節(jié)點(diǎn)死亡時(shí)間if dead1 == 1if flag_first_dead1 == 0first_dead1 = r+1;flag_first_dead1 = 1;endend% 10%的節(jié)點(diǎn)死亡時(shí)間if dead1 == 0.1*nif flag_teenth_dead1 ==0teenth_dead1 = r+1;flag_teenth_dead1 = 1;endendif dead1 == nif flag_all_dead1 == 0all_dead1 = r+1;flag_all_dead1 = 1;endendendif N(i).E > 0energy1 = energy1 + N(i).E;%統(tǒng)計(jì)每輪節(jié)點(diǎn)總能量和N(i).type = 'N';endend%節(jié)點(diǎn)死亡，退出循環(huán)，但是一般都是與其他算法比較，由于不同算法最終循環(huán)倫次不同%所以節(jié)點(diǎn)死亡并不退出循環(huán)，而是在達(dá)到最大循環(huán)次數(shù)后才結(jié)束循環(huán)
%     if flag_all_dead == 1
%         break;
%     end%r從0開始，但是計(jì)算機(jī)中數(shù)組從1開始，所以采用r+1的下標(biāo)STATISTICS1.DEAD(r+1) = dead1;%每輪死亡節(jié)點(diǎn)數(shù)STATISTICS1.alive(r+1) = n-dead1;%每輪存活節(jié)點(diǎn)數(shù)STATISTICS1.energy(r+1) = energy1;%每輪剩余總能量cluster=0;%計(jì)數(shù)簇頭數(shù)%次循環(huán)（選簇頭）if energy1 > 0%網(wǎng)絡(luò)還剩余能量，這一步在節(jié)點(diǎn)都死亡后才有用for i=1:nif(N(i).E>0)%節(jié)點(diǎn)必須有能量temp_rand=rand;     if ( (N(i).G)<=0) %該節(jié)點(diǎn)在候選集合中%選簇頭if( temp_rand <= (p/(1-p*mod(r,round(1/p)))))cluster = cluster+1;N(i).type = 'C'; %簇頭N(i).G = round(1/p)-1;%換成"S(i).G = 1"也行，只要不是0C(cluster).xd = N(i).xd;C(cluster).yd = N(i).yd;distance=sqrt( (N(i).xd-(N(n+1).xd))^2 + (N(i).yd-(N(n+1).yd))^2 );%距離基站的距離C(cluster).distance = distance;C(cluster).id = i;distanceBroad = sqrt(xm^2+ym^2);%全網(wǎng)廣播自成為簇頭,廣播TDMA時(shí)間表，距離是distanceBroadif (distanceBroad > do)N(i).E = N(i).E- (2* ETX * ctrPacketLength + Emp * ctrPacketLength*distanceBroad^4);elseN(i).E = N(i).E- (2* ETX * ctrPacketLength + Efs * ctrPacketLength*distanceBroad^2); end%簇頭向基站發(fā)送數(shù)據(jù)包能耗if (distance > do)N(i).E = N(i).E- ( ETX*(packetLength) + Emp * packetLength*distance^4);elseN(i).E = N(i).E- ( ETX*(packetLength) + Efs * packetLength*distance^2); endend     endend endendSTATISTICS1.COUNTCHS(r+1) = cluster;%次循環(huán)(成簇)if energy1 > 0%網(wǎng)絡(luò)還剩余能量，這一步在節(jié)點(diǎn)都死亡后才生效for i = 1:nif ( N(i).type=='N' && N(i).E>0 ) %普通節(jié)點(diǎn)%如果有簇頭存在if(cluster>= 1)length = zeros(cluster,1);%加入最近的簇頭for c = 1:clusterlength(c) = sqrt( (N(i).xd - C(c).xd)^2 + (N(i).yd - C(c).yd)^2 );%%接收簇頭發(fā)來的廣播和TDMA時(shí)間表的消耗N(i).E = N(i).E - 2*ERX * ctrPacketLength;end[min_dis, min_dis_cluster] = min(length);%min_dis返回最短距離，min_dis_cluster返回最短距離對(duì)應(yīng)的簇號(hào)%向簇頭發(fā)送加入消息、向簇頭發(fā)送數(shù)據(jù)包if (min_dis > do)N(i).E = N(i).E - ( ETX*ctrPacketLength + Emp * ctrPacketLength*min_dis^4); %向簇頭發(fā)送加入控制消息N(i).E = N(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇發(fā)送頭數(shù)據(jù)包elseN(i).E = N(i).E - ( ETX*ctrPacketLength + Efs*ctrPacketLength*min_dis^2); %向簇頭發(fā)送加入控制消息N(i).E = N(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇頭發(fā)送數(shù)據(jù)包end%接收簇頭確認(rèn)加入控制消息N(i).E = N(i).E - ERX*ctrPacketLength;%簇頭接收簇成員的加入消息和數(shù)據(jù)包，簇頭向簇成員發(fā)送確認(rèn)加入的信息if(min_dis > 0)N(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ERX *ctrPacketLength ;%接收加入消息N(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ( (ERX + EDA)*packetLength ); %接受簇成員發(fā)來的數(shù)據(jù)包%簇頭向簇成員發(fā)送確認(rèn)加入的信息if (min_dis > do)N(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ( ETX*ctrPacketLength + Emp * ctrPacketLength*min_dis^4);elseN(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ( ETX*ctrPacketLength + Efs * ctrPacketLength*min_dis^2);endend%如果本輪沒有簇頭，則普通節(jié)點(diǎn)直接將數(shù)據(jù)直接發(fā)送給基站elsemin_dis = sqrt((N(i).xd-N(n+1).xd)^2 + (N(i).yd-N(n+1).yd)^2);if min_dis > doN(i).E = N(i).E- (ETX*packetLength + Emp*packetLength*min_dis^4);elseN(i).E = N(i).E- (ETX*packetLength + Efs*packetLength*min_dis^2);endendendendend
end
%% 繪圖顯示
figure;
plot(0:rmax, STATISTICS.DEAD, 'r',0:rmax, STATISTICS1.DEAD, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('死亡節(jié)點(diǎn)數(shù)');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.alive, 'r',0:rmax, STATISTICS1.alive, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('活動(dòng)節(jié)點(diǎn)數(shù)');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.energy, 'r',0:rmax, STATISTICS1.energy, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('剩余能量');
legend('mypso','Leach')  
figure;
plot(0:rmax, STATISTICS.COUNTCHS, 'r',0:rmax, STATISTICS1.COUNTCHS, 'g', 'LineWidth', 2);
xlabel('輪數(shù)'); ylabel('簇頭數(shù)');
legend('mypso','Leach')

三、實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)1、不考慮控制包，對(duì)比LEACH與TL-EBC

實(shí)驗(yàn)2、考慮控制包，對(duì)比LEACH與TL-EBC

結(jié)果分析

??很明顯，這樣設(shè)計(jì)網(wǎng)絡(luò)性能是得到了優(yōu)化的，不論是否考慮控制包，對(duì)比LEACH算法，網(wǎng)絡(luò)性能都得到了很大的提升，一般網(wǎng)絡(luò)節(jié)點(diǎn)死亡超過30%，其網(wǎng)絡(luò)就認(rèn)定工作性能不達(dá)標(biāo)，TL-EBC算法實(shí)現(xiàn)的網(wǎng)絡(luò)節(jié)點(diǎn)幾乎同時(shí)死亡，使網(wǎng)絡(luò)有效運(yùn)轉(zhuǎn)壽命得到很大的提升。
??其實(shí)讀者會(huì)發(fā)現(xiàn)，是否加入控制包相關(guān)耗能，對(duì)LEACH算法的影響小，對(duì)TL-EBC算法的影響大，其實(shí)不是這樣的（這樣認(rèn)為說明讀者對(duì)該領(lǐng)域的了解太淺顯了），這是因?yàn)門L-EBC算法的網(wǎng)絡(luò)幾乎直到網(wǎng)絡(luò)死亡都是全部節(jié)點(diǎn)在工作，LEACH算法在大部分時(shí)間都是部分節(jié)點(diǎn)在工作，只有部分節(jié)點(diǎn)運(yùn)轉(zhuǎn)的網(wǎng)絡(luò)是有很多缺點(diǎn)的。
??本人對(duì)控制包相關(guān)的功能考慮的比較充分，下面我將論文的實(shí)驗(yàn)效果與我的實(shí)驗(yàn)效果做對(duì)比。
加入控制包相關(guān)能耗：

不加入控制包相關(guān)能耗：

??對(duì)比分析，本人沒加入控制包的圖比作者的好，加入了控制包的圖沒作者的好。不難發(fā)現(xiàn)，論文原作者是折中處理過了的，即加入了控制包，但控制包相關(guān)功能只考慮了一部分。

四、參考文獻(xiàn)

[1]蔣暢江,向敏,唐賢倫.基于PSO的無線傳感器網(wǎng)絡(luò)分簇路由協(xié)議[J].計(jì)算機(jī)工程,2012,38(17):59-62.
[2]Latiff N M A, Tsimenidis C C, Sharif B S. Energy-aware Clustering for Wireless Sensor Networks Using Particle Swarm Optimization[C]//Proc. of the 18th IEEE International Symposium on Personal, Indoor and Mobile Radio Communications. Athens, Greece: IEEE Press, 2007: 1-5.