本文模拟仿真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:n
S(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 <= 0
dead = dead+1;
%第一个节点死亡时间
if dead == 1
if flag_first_dead == 0
first_dead = r+1;
flag_first_dead = 1;
end
end
% 10%的节点死亡时间
if dead == 0.1*n
if flag_teenth_dead ==0
teenth_dead = r+1;
flag_teenth_dead = 1;
end
end
if dead == n
if flag_all_dead == 0
all_dead = r+1;
flag_all_dead = 1;
end
end
end
if S(i).E > 0
energy = energy + S(i).E;%统计每轮节点总能量和
S(i).type = 'N';
j = j+1;
Node(j,:) = [S(i).xd S(i).yd];%存活节点集合
end
end
STATISTICS.DEAD(r+1) = dead;%每轮死亡节点数
STATISTICS.alive(r+1) = n-dead;%每轮存活节点数
STATISTICS.energy(r+1) = energy;%每轮剩余总能量
if n-dead == 0
STATISTICS.COUNTCHS(r+1) = 0;
continue;
end
%% 4.1、应用PSO选簇头
%候选簇头集合
j = 0;
for i = 1:n
if 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;%候选簇头对应的节点编号
else
if 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;%候选簇头对应的节点编号
end
end
end
M = 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);
end
Iter = 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:nPop
particle.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:nPop
particle.Velocity{i}(particle.Velocity{i}>20) = 20;
particle.Velocity{i}(particle.Velocity{i}<-20) = -20;
end
% 位置更新
for i = 1:nPop
particle.Position{i} = particle.Position{i} + particle.Velocity{i};
end
% 边界位置处理
for i = 1:nPop
particle.Position{i}(particle.Position{i}>200)=200;
particle.Position{i}(particle.Position{i}<0)=0;
end
%根据距离最近候选簇头位置调整粒子位置,计算相应适应度
for i = 1:nPop
distance = 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);
end
Iter = Iter+1;
end
%根据距离最近候选簇头位置得到簇头信息
distance = pdist2(C.d,GlobalBestPosition);
[~, index2] = min(distance);
%如果簇头重复，需要调整(这一步处理的不够好)
length_index2 = length(unique(index2));
cha = K-length_index2;%重复个数
while cha ~= 0
for i = 1:K-1
for j = i+1:K
if 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);
end
cha = cha-1;
end
end
end
end
for 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:K
j = C.id(index2(i));
S(j).type = 'C';
end
Each_cluster = zeros(K,2);%第一列放每个簇的能量最大的节点的能量信息，第二列放节点标号,行对应簇号
%预成簇
for i = 1:n
if ( 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
end
end
%若有空簇，将其候选总簇头节点假定为1
for i = 1:K
if Each_cluster(i,2) == 0
Each_cluster(i,2) = 1;
end
end
%每个候选总簇头至所有簇头距离和
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:K
CostofHCH(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:n
if ( 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); %向簇发送头数据包
else
S(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 ); %接受簇成员发来的数据包
end
end
end
%% 4.4簇头和总簇头通信
TCnumber = Each_cluster(index3,2);%总簇头编号
for i = 1:K
j = 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);
else
S(j).E = S(j).E- ( ETX * (packetLength) + Efs * packetLength*distance^2);
end
%总簇头接收簇头的的数据包(不融合数据)
if(distance > 0)
S(TCnumber).E = S(TCnumber).E - ERX * packetLength; %接受簇头发来的数据包
end
end
%% 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);
else
S(TCnumber).E = S(TCnumber).E- (K * ETX * (packetLength) + Efs * packetLength*distance^2);
end
STATISTICS.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
%每过一个周期，重置G
if(mod(r, round(1/p)) == 0)
for i = 1:n
N(i).G = 0;
end
end
dead1 = 0;%计算死亡节点数
energy1 = 0;%每轮节点剩余总能量
for i = 1:n
%检查有无死亡节点
if N(i).E <= 0
dead1 = dead1+1;
%第一个节点死亡时间
if dead1 == 1
if flag_first_dead1 == 0
first_dead1 = r+1;
flag_first_dead1 = 1;
end
end
% 10%的节点死亡时间
if dead1 == 0.1*n
if flag_teenth_dead1 ==0
teenth_dead1 = r+1;
flag_teenth_dead1 = 1;
end
end
if dead1 == n
if flag_all_dead1 == 0
all_dead1 = r+1;
flag_all_dead1 = 1;
end
end
end
if N(i).E > 0
energy1 = energy1 + N(i).E;%统计每轮节点总能量和
N(i).type = 'N';
end
end
%节点死亡，退出循环，但是一般都是与其他算法比较，由于不同算法最终循环伦次不同
%所以节点死亡并不退出循环，而是在达到最大循环次数后才结束循环
% 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:n
if(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"也行，只要不是0
C(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时间表，距离是distanceBroad
if (distanceBroad > do)
N(i).E = N(i).E- (2* ETX * ctrPacketLength + Emp * ctrPacketLength*distanceBroad^4);
else
N(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);
else
N(i).E = N(i).E- ( ETX*(packetLength) + Efs * packetLength*distance^2);
end
end
end
end
end
end
STATISTICS1.COUNTCHS(r+1) = cluster;
%次循环(成簇)
if energy1 > 0%网络还剩余能量，这一步在节点都死亡后才生效
for i = 1:n
if ( N(i).type=='N' && N(i).E>0 ) %普通节点
%如果有簇头存在
if(cluster>= 1)
length = zeros(cluster,1);
%加入最近的簇头
for c = 1:cluster
length(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); %向簇发送头数据包
else
N(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
%如果本轮没有簇头，则普通节点直接将数据直接发送给基站
else
min_dis = sqrt((N(i).xd-N(n+1).xd)^2 + (N(i).yd-N(n+1).yd)^2);
end
end
end
end
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:n
S(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 <= 0
dead = dead+1;
%第一个节点死亡时间
if dead == 1
if flag_first_dead == 0
first_dead = r+1;
flag_first_dead = 1;
end
end
% 10%的节点死亡时间
if dead == 0.1*n
if flag_teenth_dead ==0
teenth_dead = r+1;
flag_teenth_dead = 1;
end
end
if dead == n
if flag_all_dead == 0
all_dead = r+1;
flag_all_dead = 1;
end
end
end
if S(i).E > 0
energy = energy + S(i).E;%统计每轮节点总能量和
S(i).type = 'N';
j = j+1;
Node(j,:) = [S(i).xd S(i).yd];%存活节点集合
end
end
STATISTICS.DEAD(r+1) = dead;%每轮死亡节点数
STATISTICS.alive(r+1) = n-dead;%每轮存活节点数
STATISTICS.energy(r+1) = energy;%每轮剩余总能量
if n-dead == 0
STATISTICS.COUNTCHS(r+1) = 0;
continue;
end
%% 4.1、应用PSO选簇头
%候选簇头集合
j = 0;
for i = 1:n
if 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;%候选簇头对应的节点编号
else
if 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;%候选簇头对应的节点编号
end
end
end
M = 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);
end
Iter = 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:nPop
particle.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:nPop
particle.Velocity{i}(particle.Velocity{i}>20) = 20;
particle.Velocity{i}(particle.Velocity{i}<-20) = -20;
end
% 位置更新
for i = 1:nPop
particle.Position{i} = particle.Position{i} + particle.Velocity{i};
end
% 边界位置处理
for i = 1:nPop
particle.Position{i}(particle.Position{i}>200)=200;
particle.Position{i}(particle.Position{i}<0)=0;
end
%根据距离最近候选簇头位置调整粒子位置,计算相应适应度
for i = 1:nPop
distance = 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);
end
Iter = Iter+1;
end
%根据距离最近候选簇头位置得到簇头信息
distance = pdist2(C.d,GlobalBestPosition);
[~, index2] = min(distance);
%如果簇头重复，需要调整(这一步处理的不够好)
length_index2 = length(unique(index2));
cha = K-length_index2;%重复个数
while cha ~= 0
for i = 1:K-1
for j = i+1:K
if 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);
end
cha = cha-1;
end
end
end
end
distanceBroad = 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);
else
S(C.id(index2(i))).E = S(C.id(index2(i))).E - ( ETX*ctrPacketLength + Efs*packetLength*distanceBroad^2);
end
end
%重置群体和个体最佳适应度、位置
GlobalBestCost = inf;
GlobalBestPosition = [];
particle.Best.Cost = ones(1,nPop)*inf;
particle.Best.Position = cell(1,nPop);
%% 4.2选总簇头
%确定簇头
for i = 1:K
j = C.id(index2(i));
S(j).type = 'C';
end
Each_cluster = zeros(K,2);%第一列放每个簇的能量最大的节点的能量信息，第二列放节点标号,行对应簇号
%预成簇
for i = 1:n
if ( 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
end
end
%若有空簇，将其候选总簇头节点假定为1
for i = 1:K
if Each_cluster(i,2) == 0
Each_cluster(i,2) = 1;
end
end
%每个候选总簇头至所有簇头距离和
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:K
CostofHCH(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:n
if ( 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); %向簇发送头数据包
else
S(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 ); %接受簇成员发来的数据包
end
end
end
%% 4.4簇头和总簇头通信
TCnumber = Each_cluster(index3,2);%总簇头编号
for i = 1:K
j = 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);
else
S(j).E = S(j).E- ( ETX * (packetLength) + Efs * packetLength*distance^2);
end
%总簇头接收簇头的的数据包(不融合数据)
if(distance > 0)
S(TCnumber).E = S(TCnumber).E - ERX * packetLength; %接受簇头发来的数据包
end
end
%% 4.5总簇头和基站通信
%总簇头广播时隙表,距离是maxdistance
maxdistance = 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);
else
S(TCnumber).E = S(TCnumber).E- ( ETX * ctrPacketLength + Efs * ctrPacketLength*maxdistance^2);
end
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);
else
S(TCnumber).E = S(TCnumber).E- (K * ETX * (packetLength) + Efs * packetLength*distance^2);
end
STATISTICS.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
%每过一个周期，重置G
if(mod(r, round(1/p)) == 0)
for i = 1:n
N(i).G = 0;
end
end
dead1 = 0;%计算死亡节点数
energy1 = 0;%每轮节点剩余总能量
for i = 1:n
%检查有无死亡节点
if N(i).E <= 0
dead1 = dead1+1;
%第一个节点死亡时间
if dead1 == 1
if flag_first_dead1 == 0
first_dead1 = r+1;
flag_first_dead1 = 1;
end
end
% 10%的节点死亡时间
if dead1 == 0.1*n
if flag_teenth_dead1 ==0
teenth_dead1 = r+1;
flag_teenth_dead1 = 1;
end
end
if dead1 == n
if flag_all_dead1 == 0
all_dead1 = r+1;
flag_all_dead1 = 1;
end
end
end
if N(i).E > 0
energy1 = energy1 + N(i).E;%统计每轮节点总能量和
N(i).type = 'N';
end
end
%节点死亡，退出循环，但是一般都是与其他算法比较，由于不同算法最终循环伦次不同
%所以节点死亡并不退出循环，而是在达到最大循环次数后才结束循环
% 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:n
if(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"也行，只要不是0
C(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时间表，距离是distanceBroad
if (distanceBroad > do)
N(i).E = N(i).E- (2* ETX * ctrPacketLength + Emp * ctrPacketLength*distanceBroad^4);
else
N(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);
else
N(i).E = N(i).E- ( ETX*(packetLength) + Efs * packetLength*distance^2);
end
end
end
end
end
end
STATISTICS1.COUNTCHS(r+1) = cluster;
%次循环(成簇)
if energy1 > 0%网络还剩余能量，这一步在节点都死亡后才生效
for i = 1:n
if ( N(i).type=='N' && N(i).E>0 ) %普通节点
%如果有簇头存在
if(cluster>= 1)
length = zeros(cluster,1);
%加入最近的簇头
for c = 1:cluster
length(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); %向簇发送头数据包
else
N(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);
else
N(C(min_dis_cluster).id).E = N(C(min_dis_cluster).id).E - ( ETX*ctrPacketLength + Efs * ctrPacketLength*min_dis^2);
end
end
%如果本轮没有簇头，则普通节点直接将数据直接发送给基站
else
min_dis = sqrt((N(i).xd-N(n+1).xd)^2 + (N(i).yd-N(n+1).yd)^2);
if min_dis > do
N(i).E = N(i).E- (ETX*packetLength + Emp*packetLength*min_dis^4);
else
N(i).E = N(i).E- (ETX*packetLength + Efs*packetLength*min_dis^2);
end
end
end
end
end
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.
9765
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
