clc;
clear;
close all;
warning off;
addpath(genpath(pwd));
N = 100; %节点数
W = 200; %网络长度
L = 200; %网络宽度
Ei = 2; %每个节点的初始能量(焦耳)
CHpl = 3000; %每轮簇头的数据包大小(位)
p = 5/100; %期望的簇头百分比
R = 50; %群集的范围
pMin = 10^-4; %可能的最低CH_prop
num_rounds = 2000; %最大模拟轮数
NonCHpl = 200; %每轮正常节点的数据包大小(位)
Tsetup = 4; %安装阶段的平均时间(秒)
Tss = 10; %稳态阶段的平均时间(秒)
Etrans = 1.0000e-05; %传输一位的能量
Erec = 1.0000e-05; %接收一位的能量
Eagg = 1.0000e-07; %数据聚合能量
Efs = 0.3400e-9; %自由空间模型放大器的能量
% Position of sink
SX = W/2; SY = L/2;
net = [ones(1,N);rand([1,N])*W;rand([1,N])*L];
% Preallocation for energy calculations
E = Ei*ones(1,N); % Energy left in each node
EH = zeros(1,num_rounds);
% Preallocation for dead nodes calculations
Ecrit = 0; % Critical energy left in node to call it alive
Ncrit = fix((95/100)*N); % Critical number for dead nodes to stop simulation
Dead = false(1,N); % Status of all nodes 0:Alive 1:Dead
DeadH = zeros(1,num_rounds);
% Preallocation for Bits sent calculations
BitsH = zeros(1,num_rounds);
figure('Position',[34 30 792 613]);
for r=1:num_rounds % iterating on each round
[net(1,:),CH] = Leach_algo(net(1,:),Dead,p,r);
tmp = find(CH);
for i=1:N
if isempty(net(2,CH))
else
[~,aa]=min(sqrt((net(2,CH) - net(2,i)).^2 + (net(3,CH) - net(3,i)).^2));
net(1,i) = tmp(aa);
end
end
EH(r) = sum(E);
numClust = length(find(CH));
D = sqrt((net(2,CH) - SX).^2 + (net(3,CH) - SY).^2);
E(CH) = E(CH) - (((Etrans+Eagg)*CHpl)+(Efs*CHpl*(D.^ 2))+(NonCHpl*Erec*round(N/numClust)));
% second rest of nodes
rest = N-numClust-sum(double(Dead));
mD = zeros(1,rest); tmp = net(2:3,~CH&~Dead);
for i=1:rest, mD(i) = fun(tmp(1,i),tmp(2,i),net,CH,SX,SY); end
E(~CH&~Dead) = E(~CH&~Dead) - ((NonCHpl*Etrans) + (Efs*CHpl*(mD.^2)) + ((Erec+Eagg)*CHpl));
%finally updating alive status to all nodes
E(Dead) = 0;
Dead(E<=Ecrit) = true ; DeadH(r)=sum(double(Dead));
%%%% sent bits %%%%
BitsH(r+1) = BitsH(r) + numClust*CHpl + rest*NonCHpl;
net = DrawNet(net,N,CH,Dead,SX,SY,1);
title(['Normal nodes:Black ---- CH:Red ---- Dead:Empty circle --- round (',num2str(r),')']);
drawnow
if DeadH(r)>=Ncrit,break;end % Stop simulation when 5% or less is alive
end
close all
T_L = (Tsetup+Tss)*(0:r-1);
EH = EH(1:r); EHdis_L = (N*Ei)-EH;
DeadH = DeadH(1:r); AliveH_L = N-DeadH;
BitsH_L = BitsH(2:r+1);
net = [zeros(1,N);net(2:3,:)];
cost = zeros(1,N);
for i=1:N
Dist = sqrt(((net(2,:)-net(2,i)).^2) + ((net(3,:)-net(3,i)).^2));
Snbr = Dist <= R;
cost(i) = sum(Dist(Snbr))/(sum(Snbr)-1);
end
E = Ei*ones(1,N);
EH = zeros(1,num_rounds);
Ecrit = 0;
Ncrit = fix((98/100)*N);
Dead = false(1,N);
DeadH = zeros(1,num_rounds);
% Preallocation for Bits sent calculations
BitsH = zeros(1,num_rounds);
figure('Position',[34 30 792 613]);
% Simulating for each round
for r=1:num_rounds
[CH,net] = HEED_algo(R,Dead,p,pMin,E,Ei,net,cost);
EH(r) = sum(E); %get all energy left in all nodes
% first CH
numClust = length(find(CH));
D = sqrt((net(2,CH) - SX).^2 + (net(3,CH) - SY).^2);
E(CH) = E(CH) - (((Etrans+Eagg)*CHpl)+(Efs*CHpl*(D.^ 2))+(NonCHpl*Erec*round(N/numClust)));
% second rest of nodes
rest = N-numClust-sum(double(Dead));
mD = zeros(1,rest); tmp = net(2:3,~CH&~Dead);
for i=1:rest, mD(i) = funH(tmp(1,i),tmp(2,i),net,CH,SX,SY); end
E(~CH&~Dead) = E(~CH&~Dead) - ((NonCHpl*Etrans) + (Efs*CHpl*(mD.^2)) + ((Erec+Eagg)*CHpl));
%finally updating alive status to all nodes
E(Dead) = 0; CH(Dead) = false;
Dead(E<=Ecrit) = true ; DeadH(r)=sum(double(Dead));
BitsH(r+1) = BitsH(r) + numClust*CHpl + rest*NonCHpl;
net = DrawNet(net,N,CH,Dead,SX,SY,1);
title(['Normal nodes:Black ---- CH:Red ---- Dead:Empty circle --- round (',num2str(r),')']);
drawnow
if sum(Dead)>=Ncrit,break;end
end
T = (Tsetup+Tss)*(0:r-1);
EH = EH(1:r); EHdis = (N*Ei)-EH;
DeadH = DeadH(1:r); AliveH = N-DeadH;
BitsH = BitsH(2:r+1);
close all
figure('Position',[131 59 792 613]);
plot(T,EHdis,'-x')
plot(T,EHdis,'-x',T_L,EHdis_L,'-x');
xlabel('Time (s)'); ylabel('Energy (j)')
title('Total energy dissipated')
legend('HEED','LEACH')
figure('Position',[298 66 792 613]);
plot(T,AliveH,'-x')
plot(T,AliveH,'-x',T_L,AliveH_L,'-x');
xlabel('Time (s)'); ylabel('No of live Sensors nodes')
title('Life time of sensor nodes')
legend('HEED','LEACH')
figure('Position',[511 61 792 613]);
plot(T,BitsH,'-x')
plot(T,BitsH,'-x',T_L,BitsH_L,'-x');
xlabel('Time (s)'); ylabel('Throughput (no of packets)')
title(['Throughput (' num2str(N) ' Sensor nodes)'])
legend('HEED','LEACH')