一、内容说明
本文所采用的基于Kmeans算法的分簇路由协议是本人所想的一个简单的改进(不是复现的具体哪一篇论文的内容)，然后与Leach算法进行了指标对比。
Kmeans算法及其复现可以参考以下内容：
https://blog.csdn.net/qq_45004419/article/details/130636784?spm=1001.2014.3001.5502
Leach算法及其复现可以参考以下内容：
https://blog.csdn.net/qq_45004419/article/details/129584285?spm=1001.2014.3001.5502
二、算法简述
先成簇，后选簇头。
成簇阶段采用Kmeans算法，聚类数就是簇头数，簇头数采用SEP协议中给的最佳簇数公式:
成簇后，选择簇内距离聚类中心距离最短的点为聚类中心（因为只根据这个因素选簇头，所以本算法比较简单），当然，聚类中心选取的指标还可以增加关于能量、度等等指标。
簇间簇内通信和Leach协议一样。
思考：Kmeans算法是一种无监督的聚类算法，而无线传感器网络成簇也是一种聚类的形式，所以Kmeans算法可以很好的应用于无线传感器网络成簇方面，但是它们又有不同点，Kmenas算法及其相关的改进算法，如FCM算法、AFKM算法、MKM算法等，其聚类考虑的因素都是1～2个，这很难达到无线传感器网络负载均衡的要求，无线传感器网络通常需要考虑能量、各种距离、通信代价等相对复杂的因素，所以仅仅只与Kmeans算法结合是达不到如今相关领域研究的要求的。
但是Kemans可以作为一个改进点，聚类数的选取也是可以进行改进的，Kmenas初始化中心也可以进行改进，常用的改进点就是与启发式算法结合。
三、算法复现
本文提出的算法
%% 清空环境变量
clear;clc;
close all;
%% 1、初始参数设定模块
%监测区域大小（单位：米）
xm = 100;
ym = 100;
%基站的x和y坐标
sink.x = 0.5 * xm;
sink.y = 0.5 * ym;
%场地中的节点数
n = 100;
%最佳簇头概率(根据最优簇头公式得来)
p=0.1;
%最大循环轮数
rmax = 2000;
%能量模型（所有值均以焦耳为单位）
Eo = 0.5;%初始能量
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 = 4000;%数据包长度
ctrPacketLength = 100;%控制包长度
%% 2、无线传感器网络模型产生模块
%构建无线传感器网络,在区域内均匀投放所有节点,并画出图形
figure;
for i = 1:n
S(i).xd = rand(1,1)*xm;%坐标
S(i).yd = rand(1,1)*ym;
S(i).id = i;%节点的ID
S(i).class = 0;%所属的聚类号
S(i).E = Eo;%初始化节点能量
S(i).cond = 1;% 节点的当前状态：存活 1；死亡 0
%最初没有簇头只有簇成员
S(i).type = 'N';%簇成员(普通节点)
plot(S(i).xd,S(i).yd,'o');
hold on;
end
%加入基站信息
S(n+1).xd = sink.x;
S(n+1).yd = sink.y;
plot(S(n+1).xd,S(n+1).yd,'p');
%% 3、网络运行模块
% 死亡节点数指标
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;%每轮节点剩余总能量
for i = 1:n
%检查有无死亡节点
if S(i).E <= 0
S(i).cond = 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
end
X = zeros(n-dead,2);%质心坐标矩阵
j = 0;
for i = 1:n
if S(i).E > 0
energy = energy + S(i).E;%统计每轮节点总能量和
S(i).type = 'N';
j = j+1;
X(j,:) = [S(i).xd,S(i).yd];
end
end
%节点死亡，退出循环，但是一般都是与其他算法比较，由于不同算法最终循环伦次不同
%所以节点死亡并不退出循环，而是在达到最大循环次数后才结束循环
% if flag_all_dead == 1
% break;
% end
%r从0开始，但是计算机中数组从1开始，所以采用r+1的下标
STATISTICS.DEAD(r+1) = dead;%每轮死亡节点数
STATISTICS.alive(r+1) = n-dead;%每轮存活节点数
STATISTICS.energy(r+1) = energy;%每轮剩余总能量
%运用kmeans算法
k = n-dead;%质心个数（存活的节点）
C_number = round(sqrt(k/(2*pi))*2/0.765);%聚类数，该公式来源于SEP协议中的最佳簇数
if C_number > 0
[idx,center] = kmeans(X,C_number);%center返回聚类中心坐标,idx返回对应聚类号
[class_number,~] = size(center);%最终的聚类数
else
idx = [];
center = [];
class_number = 0;
end
j = 0;
%给节点标记聚类号
for i = 1:n
if S(i).cond == 1
j = j+1;
S(i).class = idx(j);
end
end
if class_number >= 1
%选簇头
for k = 1:class_number
min_dist = 1e6;%最小距离
for i = 1:n
if S(i).cond == k
dist = sqrt( (S(i).xd - center(k,1))^2 + (S(i).yd - center(k,2))^2 );
if dist<min_dist
min_dist = dist;
id = i;
end
end
end
S(id).type = 'C';%簇头标记
C(k).xd = S(id).xd;
C(k).yd = S(id).yd;
C(k).id = id;
distance=sqrt( (S(id).xd-(S(n+1).xd))^2 + (S(id).yd-(S(n+1).yd))^2 );%距离基站的距离
C(k).distance = distance;
%簇头向基站发送数据包能耗
if (distance > do)
S(id).E = S(id).E- ( (ETX+EDA)*(packetLength) + Emp * packetLength*distance^4);
else
S(id).E = S(id).E- ( (ETX+EDA)*(packetLength) + Efs * packetLength*distance^2);
end
end
STATISTICS.COUNTCHS(r+1) = class_number;
packets_TO_BS = packets_TO_BS + class_number;
if energy > 0%网络还剩余能量，这一步在节点都死亡后才生效
for k = 1:class_number
for i = 1:n
if ( S(i).class == k && S(i).type=='N' && S(i).E>0) %普通节点
dist = sqrt( (S(i).xd - C(k).xd)^2 + (S(i).yd - C(k).yd)^2 );
%向簇头发送数据包
if (dist > do)
S(i).E = S(i).E - ( ETX*packetLength + Emp*packetLength*dist^4); %向簇发送头数据包
else
S(i).E = S(i).E - ( ETX*packetLength + Efs*packetLength*dist^2); %向簇头发送数据包
end
%簇头接收簇成员的数据包
if(dist > 0)
S(C(k).id).E = S(C(k).id).E - ( (ERX + EDA)*packetLength ); %接受簇成员发来的数据包
end
packets_TO_CH = packets_TO_CH + 1;
end
end
end
end
STATISTICS.PACKETS_TO_CH(r+1) = packets_TO_CH;
STATISTICS.PACKETS_TO_BS(r+1) = packets_TO_BS;
end
end
%% 绘图显示
figure;
plot(0:rmax, STATISTICS.DEAD, 'r', 'LineWidth', 2);
xlabel('轮数'); ylabel('死亡节点数');
figure;
plot(0:rmax, STATISTICS.alive, 'b', 'LineWidth', 2);
xlabel('轮数'); ylabel('活动节点数');
figure;
plot(0:length(STATISTICS.PACKETS_TO_BS)-1, STATISTICS.PACKETS_TO_BS, 'm', 'LineWidth', 2);
xlabel('轮数'); ylabel('基站收到数据包数');
figure;
plot(0:length(STATISTICS.COUNTCHS)-1, STATISTICS.COUNTCHS, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('簇头数');
figure;
plot(0:rmax, STATISTICS.energy, 'k', 'LineWidth', 2);
xlabel('轮数'); ylabel('剩余能量');
本文提出的算法与Leach算法对比
注：因为控制包相关的能耗小，实验没有考虑与控制包相关的功能。
%% 清空环境变量
clear;clc;
close all;
%% 1、初始参数设定模块
%监测区域大小（单位：米）
xm = 100;
ym = 100;
%基站的x和y坐标
sink.x = 0.5 * xm;
sink.y = 0.5 * ym;
%场地中的节点数
n = 100;
%最佳簇头概率(根据最优簇头公式得来)
p=0.1;
%最大循环轮数
rmax = 2000;
%能量模型（所有值均以焦耳为单位）
Eo = 0.5;%初始能量
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 = 4000;%数据包长度
ctrPacketLength = 100;%控制包长度
%% 2、无线传感器网络模型产生模块
%构建无线传感器网络,在区域内均匀投放所有节点,并画出图形
figure;
for i = 1:n
S(i).xd = rand(1,1)*xm;%坐标
S(i).yd = rand(1,1)*ym;
S(i).id = i;%节点的ID
S(i).class = 0;%所属的聚类号
S(i).E = Eo;%初始化节点能量
S(i).cond = 1;% 节点的当前状态：存活 1；死亡 0
%最初没有簇头只有簇成员
S(i).type = 'N';%簇成员(普通节点)
plot(S(i).xd,S(i).yd,'o');
hold on;
N(i).xd = rand(1,1)*xm;%坐标
N(i).yd = rand(1,1)*ym;
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、网络运行模块(Based_on_Kmeans)
% 死亡节点数指标
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;%每轮节点剩余总能量
for i = 1:n
%检查有无死亡节点
if S(i).E <= 0
S(i).cond = 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
end
X = zeros(n-dead,2);%质心坐标矩阵
j = 0;
for i = 1:n
if S(i).E > 0
energy = energy + S(i).E;%统计每轮节点总能量和
S(i).type = 'N';
j = j+1;
X(j,:) = [S(i).xd,S(i).yd];
end
end
%节点死亡，退出循环，但是一般都是与其他算法比较，由于不同算法最终循环伦次不同
%所以节点死亡并不退出循环，而是在达到最大循环次数后才结束循环
% if flag_all_dead == 1
% break;
% end
%r从0开始，但是计算机中数组从1开始，所以采用r+1的下标
STATISTICS.DEAD(r+1) = dead;%每轮死亡节点数
STATISTICS.alive(r+1) = n-dead;%每轮存活节点数
STATISTICS.energy(r+1) = energy;%每轮剩余总能量
%运用kmeans算法
k = n-dead;%质心个数（存活的节点）
C_number = round(sqrt(k/(2*pi))*2/0.765);%聚类数，该公式来源于SEP协议中的最佳簇数
if C_number > 0
[idx,center] = kmeans(X,C_number);%center返回聚类中心坐标,idx返回对应聚类号
[class_number,~] = size(center);%最终的聚类数
else
idx = [];
center = [];
class_number = 0;
end
j = 0;
%给节点标记聚类号
for i = 1:n
if S(i).cond == 1
j = j+1;
S(i).class = idx(j);
end
end
if class_number >= 1
%选簇头
for k = 1:class_number
min_dist = 1e6;%最小距离
for i = 1:n
if S(i).cond == k
dist = sqrt( (S(i).xd - center(k,1))^2 + (S(i).yd - center(k,2))^2 );
if dist<min_dist
min_dist = dist;
id = i;
end
end
end
S(id).type = 'C';%簇头标记
C(k).xd = S(id).xd;
C(k).yd = S(id).yd;
C(k).id = id;
distance=sqrt( (S(id).xd-(S(n+1).xd))^2 + (S(id).yd-(S(n+1).yd))^2 );%距离基站的距离
C(k).distance = distance;
%簇头向基站发送数据包能耗
if (distance > do)
S(id).E = S(id).E- ( (ETX+EDA)*(packetLength) + Emp * packetLength*distance^4);
else
S(id).E = S(id).E- ( (ETX+EDA)*(packetLength) + Efs * packetLength*distance^2);
end
end
STATISTICS.COUNTCHS(r+1) = class_number;
packets_TO_BS = packets_TO_BS + class_number;
if energy > 0%网络还剩余能量，这一步在节点都死亡后才生效
for k = 1:class_number
for i = 1:n
if ( S(i).class == k && S(i).type=='N' && S(i).E>0) %普通节点
dist = sqrt( (S(i).xd - C(k).xd)^2 + (S(i).yd - C(k).yd)^2 );
%向簇头发送数据包
if (dist > do)
S(i).E = S(i).E - ( ETX*packetLength + Emp*packetLength*dist^4); %向簇发送头数据包
else
S(i).E = S(i).E - ( ETX*packetLength + Efs*packetLength*dist^2); %向簇头发送数据包
end
%簇头接收簇成员的数据包
if(dist > 0)
S(C(k).id).E = S(C(k).id).E - ( (ERX + EDA)*packetLength ); %接受簇成员发来的数据包
end
packets_TO_CH = packets_TO_CH + 1;
end
end
end
end
STATISTICS.PACKETS_TO_CH(r+1) = packets_TO_CH;
STATISTICS.PACKETS_TO_BS(r+1) = packets_TO_BS;
end
end
%% 3、网络运行模块(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;
%簇头向基站发送数据包能耗
if (distance > do)
N(i).E = N(i).E- ( (ETX+EDA)*(packetLength) + Emp * packetLength*distance^4);
else
N(i).E = N(i).E- ( (ETX+EDA)*(packetLength) + Efs * packetLength*distance^2);
end
packets_TO_BS1 = packets_TO_BS1+1;
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 );
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
N(i).min_dis = min_dis;
N(i).min_dis_cluster = min_dis_cluster;
packets_TO_CH1 = packets_TO_CH1 + 1;
%如果本轮没有簇头，则普通节点直接将数据直接发送给基站
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
packets_TO_BS1 = packets_TO_BS1+1;
end
end
end
end
STATISTICS1.PACKETS_TO_CH(r+1) = packets_TO_CH1;
STATISTICS1.PACKETS_TO_BS(r+1) = packets_TO_BS1;
end
%% 绘图显示
figure;
plot(0:rmax, STATISTICS.DEAD, 'r',0:rmax, STATISTICS1.DEAD, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('死亡节点数');
legend('Based on Kmeans','Leach')
figure;
plot(0:rmax, STATISTICS.alive, 'r',0:rmax, STATISTICS1.alive, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('活动节点数');
legend('Based on Kmeans','Leach')
figure;
plot(0:rmax, STATISTICS.energy, 'r',0:rmax, STATISTICS1.energy, 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('剩余能量');
legend('Based on Kmeans','Leach')
% 网络在没有能量后，循环还在进行，代码优化了这一点，导致部分数据在网络死亡后的轮数不更新了，所以维度不是0:rmax
% 用numel函数计算了实际数组长度，所以下面的做图代码不同于上面的代码
figure;
plot(0:numel(STATISTICS.COUNTCHS)-1, STATISTICS.COUNTCHS, 'r',0:numel(STATISTICS.COUNTCHS)-1, STATISTICS1.COUNTCHS(1:numel(STATISTICS.COUNTCHS)), 'g', 'LineWidth', 2);
xlabel('轮数'); ylabel('簇头数');
legend('Based on Kmeans','Leach')
四、运行截图
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
