备案 填写网站信息推广普通话宣传语100字

  本文模拟仿真TL-EBC协议，将LEACH算法与TL-EBC对比。

一、算法简介

1、PSO-C(2007)

  LEACH算法选簇头是分布式随机选取的，存在很多缺陷，集中式选取簇头，一般是将能量不小于网络平均能量的节点作为候选簇头，然后再从这些候选簇头中选出适应度比较好的簇头，这是一个NP问题，一般可采用元启发算法解决。
  假设网络包含N个节点，预先定义分为K个簇头，候选簇头数为M（一般情况M>>K），则可能的分簇方式有C(M,K)种，在其中确定最佳分簇方式，是一个最优化问题。应用PSO算法解决这个问题，使每个粒子代表一种可能的分簇方式，用目标函数评价其性能，设置m个粒子组成群体在C(M,K)种可能的分簇方式种寻找最优解，使目标函数取得最小值，该目标函数定义如下：

2、TL-EBC(2012)

  TL-EBC协议在PSO-C协议的基础上加入了一个总簇头的机制，将网络设计成了两层分簇的体系结构。
  总簇头的选择综合考虑节点能量、各簇头节点距离和 与基站距离 3 个方面的因素。首先，在每个簇中选择能量 最大的节点(不包括簇头)为候选总簇头。然后，计算每个候选总簇头的评价函数，具有最小评价函数值的候选总簇 头成为总簇头。评价函数定义如下：

二、仿真与分析

1、参数设置

  将200个传感器节点随机部署在200m x 200m的范围内，基站坐标(100,275)，位于网络外部。数据包长度为2000bit；控制包长度50bit；簇头的数量设定为节点总数的5%，则K=10；总簇头的数据不进行压缩。PSO算法的参数设置为：种群规模Q=30，c1=c2=2,惯性权重w = 1（这里的设计不同于论文，本人实验过，在该模型种，固定w=1，比线性变化和非线性变化的w设定更加好，不论是在网络循环中，还是在PSO循环中），最大循环次数MaxIter = 50（论文种的MaxIter = 30是不够的）。评价因子权重系数$\beta$ = 0.5,$\alpha$ = 0.2;$\eta$ = 0.4;$\lambda$ = 0.4。网络只进行同构网络的分析。

2、MATLAB 程序实现

实验1、不考虑控制包，对比LEACH与TL-EBC

%不考虑控制包
%% 清空环境变量
clear;clc;
close all;%% 1、初始参数设定模块
%监测区域大小（单位：米）
xm = 200;
ym = 200;%基站的x和y坐标
sink.x = 0.5 * xm;
sink.y = 275;%场地中的节点数
n = 200;%簇头概率
p=0.05;%最大循环轮数
rmax = 500;%能量模型(所有值均以焦耳为单位)
Eo = 0.1;%初始能量
ETX = 50*0.000000001;%传输能耗
ERX = 50*0.000000001;%接收能耗
EDA = 5*0.000000001;%数据聚合能量%发射放大器类型
Efs = 10*0.000000000001;%小于do
Emp = 0.0013*0.000000000001;%大于do
do = sqrt(Efs/Emp);%距离阈值%数据包和控制包长度
packetLength = 2000;%数据包长度 
ctrPacketLength = 50;%控制包长度%% 2、无线传感器网络模型产生模块
%构建无线传感器网络,在区域内均匀投放所有节点,并画出图形
figure;
for i = 1:nS(i).xd = rand(1,1)*xm;%坐标S(i).yd = rand(1,1)*ym;S(i).E = Eo;%初始化节点能量%最初没有簇头只有簇成员S(i).type = 'N';%簇成员(普通节点) plot(S(i).xd,S(i).yd,'o');hold on;N(i).xd = S(i).xd;%坐标N(i).yd = S(i).yd;N(i).G = 0;%"G = 0"表示节点本周期内未当选过簇头，当选则改为"G = 1"N(i).E = Eo;%初始化节点能量%最初没有簇头只有簇成员N(i).type = '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参数
MaxIter = 50;      % 最大迭代次数
nPop = 30;        % 种群规模
wmax = 0.9;
wmin = 0.4;
w = 1;
c1 = 2.0;         % 个人学习系数
c2 = 2.0;         % 全局学习系数
plimit = [0, 200];                % 位置参数限制
vlimit = [-20, 20];               % 速度限制%数据初始化
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、网络运行模块 
% 死亡节点数指标
first_dead = 0;%第一个死亡节点
teenth_dead = 0;%10%的死亡节点
all_dead = 0;%节点都死亡%节点死亡数指标对应标记
flag_first_dead = 0;
flag_teenth_dead = 0;
flag_all_dead = 0;%传输到基站和簇头的比特计数器
packets_TO_BS = 0;
packets_TO_CH = 0;%主循环
for r = 0:rmax dead = 0;%计算死亡节点数energy = 0;%每轮节点剩余总能量Node = [];%存活节点集合C = [];%候选节点集合j = 0;for i = 1:n%检查有无死亡节点if S(i).E <= 0dead = dead+1;%第一个节点死亡时间if dead == 1if flag_first_dead == 0first_dead = r+1;flag_first_dead = 1;endend% 10%的节点死亡时间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;%统计每轮节点总能量和S(i).type = 'N';j = j+1;Node(j,:) = [S(i).xd S(i).yd];%存活节点集合endendSTATISTICS.DEAD(r+1) = dead;%每轮死亡节点数STATISTICS.alive(r+1) = n-dead;%每轮存活节点数STATISTICS.energy(r+1) = energy;%每轮剩余总能量if n-dead == 0STATISTICS.COUNTCHS(r+1) = 0;continue;end%% 4.1、应用PSO选簇头%候选簇头集合j = 0;for i = 1:nif r ==0 %因为matlab最大精度16位，后面数据不准，第一轮会出现这种情况，所以这样编程j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候选簇头坐标C.E(j) = S(i).E;%候选簇头剩余能量C.id(j) = i;%候选簇头对应的节点编号elseif S(i).E >= energy/(n-dead) %节点能量大于或等于总的平均能量j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候选簇头坐标C.E(j) = S(i).E;%候选簇头剩余能量C.id(j) = i;%候选簇头对应的节点编号endendendM = size(C.d,1);%候选簇头数K = ceil((n-dead)*p);% 簇头数%初始化种群位置和速度for i = 1:nPop% 初始化位置C_k = randperm(M,K);particle.Position{i} = C.d(C_k,:);% 初始化速度v = zeros(K,2);particle.Velocity{i}=v;% 适应度energe_Q = sum(C.E((C_k)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%传感器簇成员节点矩阵particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = 1;%计数while Iter <= MaxIter%更新粒子个体最佳适应度、位置   for i = 1:nPop      if particle.Cost(i) < particle.Best.Cost(i) particle.Best.Cost(i) = particle.Cost(i);%更新个体历史最佳适应度particle.Best.Position{i} = particle.Position{i};%更新个体历史最佳位置end end%更新群体最佳适应度、位置if min(particle.Best.Cost) < GlobalBestCost  [GlobalBestCost, index1] = min(particle.Best.Cost);%更新群体历史最佳适应度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%根据距离最近候选簇头位置调整粒子位置,计算相应适应度for i = 1:nPopdistance = pdist2(C.d,particle.Position{i});[~, index2] = min(distance);particle.Position{i} = C.d(index2,:);% 适应度energe_Q = sum(C.E((index2)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%传感器簇成员节点矩阵particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = Iter+1;end%根据距离最近候选簇头位置得到簇头信息distance = pdist2(C.d,GlobalBestPosition);[~, index2] = min(distance);%如果簇头重复，需要调整(这一步处理的不够好)length_index2 = length(unique(index2));cha = K-length_index2;%重复个数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%重置群体和个体最佳适应度、位置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);%第一列放每个簇的能量最大的节点的能量信息，第二列放节点标号,行对应簇号%预成簇for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通节点%加入最近的簇头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返回最短距离对应的簇号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%若有空簇，将其候选总簇头节点假定为1for i = 1:Kif Each_cluster(i,2) == 0Each_cluster(i,2) = 1;endend%每个候选总簇头至所有簇头距离和sum_distance_HCH_to_CH = sum(pdist2(C.d(index2,:),reshape([S(Each_cluster(:,2)).xd S(Each_cluster(:,2)).yd],[],2)));%每个候选总簇头至基站距离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;%总簇头接收基站发来的分簇信息
%% 4.3成簇，簇内通信for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通节点%加入最近的簇头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返回最短距离对应的簇号%接收基站发来的广播的消耗（收到基站发来的分簇信息）S(i).E = S(i).E - ERX * ctrPacketLength;%向簇头发送数据包if (min_dis > do)S(i).E = S(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇发送头数据包elseS(i).E = S(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇头发送数据包end%簇头接收簇成员的数据包if(min_dis > 0)S(C.id(dis_index)).E = S(C.id(dis_index)).E - ( (ERX + EDA)*packetLength ); %接受簇成员发来的数据包endendend
%% 4.4簇头和总簇头通信TCnumber = Each_cluster(index3,2);%总簇头编号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 );%簇头向总簇头发送数据包能耗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%总簇头接收簇头的的数据包(不融合数据)if(distance > 0)S(TCnumber).E = S(TCnumber).E - ERX * packetLength; %接受簇头发来的数据包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、网络运行模块(Leach)
% 死亡节点数指标
first_dead1 = 0;%第一个死亡节点
teenth_dead1 = 0;%10%的死亡节点
all_dead1 = 0;%节点都死亡%节点死亡数指标对应标记
flag_first_dead1 = 0;
flag_teenth_dead1 = 0;
flag_all_dead1 = 0;%传输到基站和簇头的比特计数器
packets_TO_BS1 = 0;
packets_TO_CH1 = 0;%主循环
for r = 0:rmax %每过一个周期，重置Gif(mod(r, round(1/p)) == 0)for i = 1:nN(i).G = 0;endenddead1 = 0;%计算死亡节点数energy1 = 0;%每轮节点剩余总能量for i = 1:n%检查有无死亡节点if N(i).E <= 0dead1 = dead1+1;%第一个节点死亡时间if dead1 == 1if flag_first_dead1 == 0first_dead1 = r+1;flag_first_dead1 = 1;endend% 10%的节点死亡时间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;%统计每轮节点总能量和N(i).type = 'N';endend%节点死亡，退出循环，但是一般都是与其他算法比较，由于不同算法最终循环伦次不同%所以节点死亡并不退出循环，而是在达到最大循环次数后才结束循环
%     if flag_all_dead == 1
%         break;
%     end%r从0开始，但是计算机中数组从1开始，所以采用r+1的下标STATISTICS1.DEAD(r+1) = dead1;%每轮死亡节点数STATISTICS1.alive(r+1) = n-dead1;%每轮存活节点数STATISTICS1.energy(r+1) = energy1;%每轮剩余总能量cluster=0;%计数簇头数%次循环（选簇头）if energy1 > 0%网络还剩余能量，这一步在节点都死亡后才有用for i=1:nif(N(i).E>0)%节点必须有能量temp_rand=rand;     if ( (N(i).G)<=0) %该节点在候选集合中%选簇头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);%全网广播自成为簇头,广播TDMA时间表，距离是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%簇头向基站发送数据包能耗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;%次循环(成簇)if energy1 > 0%网络还剩余能量，这一步在节点都死亡后才生效for i = 1:nif ( N(i).type=='N' && N(i).E>0 ) %普通节点%如果有簇头存在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 );%%接收簇头发来的广播消耗N(i).E = N(i).E - ERX * ctrPacketLength;end[min_dis, min_dis_cluster] = min(length);%min_dis返回最短距离，min_dis_cluster返回最短距离对应的簇号%向簇头发送数据包if (min_dis > do)N(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); %向簇头发送数据包end%簇头接收簇成员的数据包if(min_dis > 0)N(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ( (ERX + EDA)*packetLength ); %接受簇成员发来的数据包end%如果本轮没有簇头，则普通节点直接将数据直接发送给基站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('轮数'); ylabel('死亡节点数');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.alive, 'r',0:rmax, STATISTICS1.alive, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('活动节点数');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.energy, 'r',0:rmax, STATISTICS1.energy, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('剩余能量');
legend('mypso','Leach')  
figure;
plot(0:rmax, STATISTICS.COUNTCHS, 'r',0:rmax, STATISTICS1.COUNTCHS, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('簇头数');
legend('mypso','Leach')

实验2、虑控制包，对比LEACH与TL-EBC

%% 清空环境变量
clear;clc;
close all;%% 1、初始参数设定模块
%监测区域大小（单位：米）
xm = 200;
ym = 200;%基站的x和y坐标
sink.x = 0.5 * xm;
sink.y = 275;%场地中的节点数
n = 200;%簇头概率
p=0.05;%最大循环轮数
rmax = 500;%能量模型(所有值均以焦耳为单位)
Eo = 0.1;%初始能量
ETX = 50*0.000000001;%传输能耗
ERX = 50*0.000000001;%接收能耗
EDA = 5*0.000000001;%数据聚合能量%发射放大器类型
Efs = 10*0.000000000001;%小于do
Emp = 0.0013*0.000000000001;%大于do
do = sqrt(Efs/Emp);%距离阈值%数据包和控制包长度
packetLength = 2000;%数据包长度 
ctrPacketLength = 50;%控制包长度%% 2、无线传感器网络模型产生模块
%构建无线传感器网络,在区域内均匀投放所有节点,并画出图形
figure;
for i = 1:nS(i).xd = rand(1,1)*xm;%坐标S(i).yd = rand(1,1)*ym;S(i).E = Eo;%初始化节点能量%最初没有簇头只有簇成员S(i).type = 'N';%簇成员(普通节点) plot(S(i).xd,S(i).yd,'o');hold on;N(i).xd = S(i).xd;%坐标N(i).yd = S(i).yd;N(i).G = 0;%"G = 0"表示节点本周期内未当选过簇头，当选则改为"G = 1"N(i).E = Eo;%初始化节点能量%最初没有簇头只有簇成员N(i).type = '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参数
MaxIter = 50;      % 最大迭代次数
nPop = 30;        % 种群规模
wmax = 0.9;
wmin = 0.4;
w = 1;
c1 = 2.0;         % 个人学习系数
c2 = 2.0;         % 全局学习系数
plimit = [0, 200];                % 位置参数限制
vlimit = [-20, 20];               % 速度限制%数据初始化
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、网络运行模块 
% 死亡节点数指标
first_dead = 0;%第一个死亡节点
teenth_dead = 0;%10%的死亡节点
all_dead = 0;%节点都死亡%节点死亡数指标对应标记
flag_first_dead = 0;
flag_teenth_dead = 0;
flag_all_dead = 0;%传输到基站和簇头的比特计数器
packets_TO_BS = 0;
packets_TO_CH = 0;%主循环
for r = 0:rmax dead = 0;%计算死亡节点数energy = 0;%每轮节点剩余总能量Node = [];%存活节点集合C = [];%候选节点集合j = 0;for i = 1:n%检查有无死亡节点if S(i).E <= 0dead = dead+1;%第一个节点死亡时间if dead == 1if flag_first_dead == 0first_dead = r+1;flag_first_dead = 1;endend% 10%的节点死亡时间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;%统计每轮节点总能量和S(i).type = 'N';j = j+1;Node(j,:) = [S(i).xd S(i).yd];%存活节点集合endendSTATISTICS.DEAD(r+1) = dead;%每轮死亡节点数STATISTICS.alive(r+1) = n-dead;%每轮存活节点数STATISTICS.energy(r+1) = energy;%每轮剩余总能量if n-dead == 0STATISTICS.COUNTCHS(r+1) = 0;continue;end%% 4.1、应用PSO选簇头%候选簇头集合j = 0;for i = 1:nif r ==0 %因为matlab最大精度16位，后面数据不准，第一轮会出现这种情况，所以这样编程j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候选簇头坐标C.E(j) = S(i).E;%候选簇头剩余能量C.id(j) = i;%候选簇头对应的节点编号elseif S(i).E >= energy/(n-dead) %节点能量大于或等于总的平均能量j = j+1;C.d(j,:) = [S(i).xd S(i).yd];%候选簇头坐标C.E(j) = S(i).E;%候选簇头剩余能量C.id(j) = i;%候选簇头对应的节点编号endendendM = size(C.d,1);%候选簇头数K = ceil((n-dead)*p);% 簇头数
%     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;% 适应度energe_Q = sum(C.E((C_k)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%传感器簇成员节点矩阵particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = 1;%计数while Iter <= MaxIter%更新粒子个体最佳适应度、位置   for i = 1:nPop      if particle.Cost(i) < particle.Best.Cost(i) particle.Best.Cost(i) = particle.Cost(i);%更新个体历史最佳适应度particle.Best.Position{i} = particle.Position{i};%更新个体历史最佳位置end end%更新群体最佳适应度、位置if min(particle.Best.Cost) < GlobalBestCost  [GlobalBestCost, index1] = min(particle.Best.Cost);%更新群体历史最佳适应度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%根据距离最近候选簇头位置调整粒子位置,计算相应适应度for i = 1:nPopdistance = pdist2(C.d,particle.Position{i});[~, index2] = min(distance);particle.Position{i} = C.d(index2,:);% 适应度energe_Q = sum(C.E((index2)));isB = ismember(Node,particle.Position{i});Cluster_member = reshape(Node(~isB),[],2);%传感器簇成员节点矩阵particle.Cost(i)=Cost(Cluster_member,particle.Position{i},K,energy,energe_Q);endIter = Iter+1;end%根据距离最近候选簇头位置得到簇头信息distance = pdist2(C.d,GlobalBestPosition);[~, index2] = min(distance);%如果簇头重复，需要调整(这一步处理的不够好)length_index2 = length(unique(index2));cha = K-length_index2;%重复个数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;%假定簇头和簇成员最远距离是对角线的1/3(懒得去后面算这个实际距离了,在分簇均匀的情况下其实是对角线的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%重置群体和个体最佳适应度、位置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);%第一列放每个簇的能量最大的节点的能量信息，第二列放节点标号,行对应簇号%预成簇for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通节点%加入最近的簇头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返回最短距离对应的簇号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%若有空簇，将其候选总簇头节点假定为1for i = 1:Kif Each_cluster(i,2) == 0Each_cluster(i,2) = 1;endend%每个候选总簇头至所有簇头距离和sum_distance_HCH_to_CH = sum(pdist2(C.d(index2,:),reshape([S(Each_cluster(:,2)).xd S(Each_cluster(:,2)).yd],[],2)));%每个候选总簇头至基站距离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;%总簇头接收基站发来的分簇信息
%% 4.3成簇，簇内通信for i = 1:nif ( S(i).type=='N' && S(i).E>0 ) %普通节点%加入最近的簇头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返回最短距离对应的簇号%接收基站发来的广播的消耗（收到基站发来的分簇信息）S(i).E = S(i).E - ERX * ctrPacketLength;%接收簇头的广播TDMA时间表信息S(i).E = S(i).E - ERX * ctrPacketLength;%向簇头发送数据包if (min_dis > do)S(i).E = S(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇发送头数据包elseS(i).E = S(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇头发送数据包end%簇头接收簇成员的数据包if(min_dis > 0)S(C.id(dis_index)).E = S(C.id(dis_index)).E - ( (ERX + EDA)*packetLength ); %接受簇成员发来的数据包endendend
%% 4.4簇头和总簇头通信TCnumber = Each_cluster(index3,2);%总簇头编号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时间表信息S(j).E = S(i).E - ERX*ctrPacketLength;%簇头向总簇头发送数据包能耗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%总簇头接收簇头的的数据包(不融合数据)if(distance > 0)S(TCnumber).E = S(TCnumber).E - ERX * packetLength; %接受簇头发来的数据包endend
%% 4.5总簇头和基站通信%总簇头广播时隙表,距离是maxdistancemaxdistance = max(pdist2(C.d(index2,:),[S(TCnumber).xd S(TCnumber).yd]));%距离最远的簇头的距离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、网络运行模块(Leach)
% 死亡节点数指标
first_dead1 = 0;%第一个死亡节点
teenth_dead1 = 0;%10%的死亡节点
all_dead1 = 0;%节点都死亡%节点死亡数指标对应标记
flag_first_dead1 = 0;
flag_teenth_dead1 = 0;
flag_all_dead1 = 0;%传输到基站和簇头的比特计数器
packets_TO_BS1 = 0;
packets_TO_CH1 = 0;%主循环
for r = 0:rmax %每过一个周期，重置Gif(mod(r, round(1/p)) == 0)for i = 1:nN(i).G = 0;endenddead1 = 0;%计算死亡节点数energy1 = 0;%每轮节点剩余总能量for i = 1:n%检查有无死亡节点if N(i).E <= 0dead1 = dead1+1;%第一个节点死亡时间if dead1 == 1if flag_first_dead1 == 0first_dead1 = r+1;flag_first_dead1 = 1;endend% 10%的节点死亡时间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;%统计每轮节点总能量和N(i).type = 'N';endend%节点死亡，退出循环，但是一般都是与其他算法比较，由于不同算法最终循环伦次不同%所以节点死亡并不退出循环，而是在达到最大循环次数后才结束循环
%     if flag_all_dead == 1
%         break;
%     end%r从0开始，但是计算机中数组从1开始，所以采用r+1的下标STATISTICS1.DEAD(r+1) = dead1;%每轮死亡节点数STATISTICS1.alive(r+1) = n-dead1;%每轮存活节点数STATISTICS1.energy(r+1) = energy1;%每轮剩余总能量cluster=0;%计数簇头数%次循环（选簇头）if energy1 > 0%网络还剩余能量，这一步在节点都死亡后才有用for i=1:nif(N(i).E>0)%节点必须有能量temp_rand=rand;     if ( (N(i).G)<=0) %该节点在候选集合中%选簇头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);%全网广播自成为簇头,广播TDMA时间表，距离是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%簇头向基站发送数据包能耗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;%次循环(成簇)if energy1 > 0%网络还剩余能量，这一步在节点都死亡后才生效for i = 1:nif ( N(i).type=='N' && N(i).E>0 ) %普通节点%如果有簇头存在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 );%%接收簇头发来的广播和TDMA时间表的消耗N(i).E = N(i).E - 2*ERX * ctrPacketLength;end[min_dis, min_dis_cluster] = min(length);%min_dis返回最短距离，min_dis_cluster返回最短距离对应的簇号%向簇头发送加入消息、向簇头发送数据包if (min_dis > do)N(i).E = N(i).E - ( ETX*ctrPacketLength + Emp * ctrPacketLength*min_dis^4); %向簇头发送加入控制消息N(i).E = N(i).E - ( ETX*packetLength + Emp*packetLength*min_dis^4); %向簇发送头数据包elseN(i).E = N(i).E - ( ETX*ctrPacketLength + Efs*ctrPacketLength*min_dis^2); %向簇头发送加入控制消息N(i).E = N(i).E - ( ETX*packetLength + Efs*packetLength*min_dis^2); %向簇头发送数据包end%接收簇头确认加入控制消息N(i).E = N(i).E - ERX*ctrPacketLength;%簇头接收簇成员的加入消息和数据包，簇头向簇成员发送确认加入的信息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 ); %接受簇成员发来的数据包%簇头向簇成员发送确认加入的信息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%如果本轮没有簇头，则普通节点直接将数据直接发送给基站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('轮数'); ylabel('死亡节点数');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.alive, 'r',0:rmax, STATISTICS1.alive, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('活动节点数');
legend('mypso','Leach')
figure;
plot(0:rmax, STATISTICS.energy, 'r',0:rmax, STATISTICS1.energy, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('剩余能量');
legend('mypso','Leach')  
figure;
plot(0:rmax, STATISTICS.COUNTCHS, 'r',0:rmax, STATISTICS1.COUNTCHS, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('簇头数');
legend('mypso','Leach')

三、实验结果

实验1、不考虑控制包，对比LEACH与TL-EBC

实验2、考虑控制包，对比LEACH与TL-EBC

结果分析

  很明显，这样设计网络性能是得到了优化的，不论是否考虑控制包，对比LEACH算法，网络性能都得到了很大的提升，一般网络节点死亡超过30%，其网络就认定工作性能不达标，TL-EBC算法实现的网络节点几乎同时死亡，使网络有效运转寿命得到很大的提升。
  其实读者会发现，是否加入控制包相关耗能，对LEACH算法的影响小，对TL-EBC算法的影响大，其实不是这样的（这样认为说明读者对该领域的了解太浅显了），这是因为TL-EBC算法的网络几乎直到网络死亡都是全部节点在工作，LEACH算法在大部分时间都是部分节点在工作，只有部分节点运转的网络是有很多缺点的。
  本人对控制包相关的功能考虑的比较充分，下面我将论文的实验效果与我的实验效果做对比。
加入控制包相关能耗：

不加入控制包相关能耗：

  对比分析，本人没加入控制包的图比作者的好，加入了控制包的图没作者的好。不难发现，论文原作者是折中处理过了的，即加入了控制包，但控制包相关功能只考虑了一部分。

四、参考文献

[1]蒋畅江,向敏,唐贤伦.基于PSO的无线传感器网络分簇路由协议[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.