失眠网 > 无线传感器网络：LEACH路由协议优化python仿真代码

无线传感器网络：LEACH路由协议优化python仿真代码

时间：2022-11-16 18:36:08

仿真结果

生成基站以及随机生成400个WSN节点

存活节点数量比较、剩余能量比较：

死亡节点比例比较

python代码

import numpy as npimport matplotlib.pyplot as pltimport random# 解决中文显示问题plt.rcParams['font.sans-serif']=['SimHei']plt.rcParams['axes.unicode_minus'] = Falseclass BaseStation: # 定义基站类x=0y=0class SensorNode: # 定义传感器节点xd=0yd=0d=0Rc = 0temp_rand=0type = 'N'selected = 'N'power = 0CH = 0flag = 1N = []Num_N = 0FN = []Num_FN = 0CN = []Num_CN = 0num_join = 0# 初始化参数xm = 100# x轴范围ym = 100# y轴范围sink = BaseStation()sink.x = 50 # 基站x轴sink.y = 125 # 基站y轴n = 400 # 节点总数p = 0.08# 簇头概率Eelec = 50*(10**(-9))Efs=10*(10**(-12))Emp=0.0013*(10**(-12))ED=5*(10**(-9))d0 = 87packetLength = 4000ctrPacketLength = 100rmax = 2000 # 迭代次数E0 = 0.5 # 初始能量Emin = 0.001 # 节点存活所需的最小能量Rmax = 15# 初始通信距离## 节点随机分布fig1 = plt.figure(dpi=80)plt.grid(linestyle="dotted")Node = []plt.scatter(sink.x, sink.y,marker='*',s=200)for i in range(n):node = SensorNode()node.xd = random.random()*xmnode.yd = random.random()*ym # 随机产生100个点node.d = ((node.xd-sink.x)**2+(node.yd-sink.y)**2)**0.5 # 节点距基站的距离node.Rc = Rmax # 节点的通信距离node.temp_rand = random.random()# rand为(0,1)的随机数node.type = 'N'# 进行选举簇头前先将所有节点设为普通节点node.selected = 'N' # ’O'：当选过簇头，N：没有node.power = E0# 初始能量node.CH = 0# 保存普通节点的簇头节点，-1代表自己是簇头node.flag = 1 # 1代表存活；0代表死亡node.N = [0 for _ in range(n) ] # 邻居节点集node.Num_N = 0# 邻居节点集个数node.FN = [0 for _ in range(n)] # 前邻节点集node.Num_FN = 0# 前邻节点集个数 = [0 for _ in range(n)] # 前簇头节点集node.Num_CN = 0# 前簇头节点集个数Node.append(node)plt.scatter(node.xd, node.yd,marker='o')plt.legend(['基站', '节点'])plt.xlabel('x', fontdict={"family": "Times New Roman", "size": 15})plt.ylabel('y', fontdict={"family": "Times New Roman", "size": 15})# save dataflag = 1################IMP_LEACH################### 迭代alive_ima_leach = np.zeros((rmax, 1)) # 每轮存活节点数re_ima_leach = np.zeros((rmax, 1))# 每轮节点总能量for r in range(rmax):final_CH=[]for i in range(n):if Node[i].flag != 0:re_ima_leach[r] = re_ima_leach[r]+Node[i].power # 更新总能量alive_ima_leach[r] = alive_ima_leach[r]+1 # 更新存活节点f = 0if alive_ima_leach[r] == 0:stop = rf = 1breakfor i in range(n):Node[i].type = 'N'# 进行选举簇头前先将所有节点设为普通节点Node[i].selected = 'N'Node[i].temp_rand = random.random() # 节点取一个(0,1)的随机值，与p比较Node[i].Rc = Rmax*Node[i].power/E0 # 节点的通信距离Node[i].CH = 0Node[i].N = np.zeros(n)# 邻居节点集Node[i].Num_N = 0 # 邻居节点集个数Node[i].FN = np.zeros(n) # 前邻节点集Node[i].Num_FN = 0# 前邻节点集个数Node[i].CN = np.zeros(n) # 前簇头节点集Node[i].Num_CN = 0# 前簇头节点集个数Node[i].num_join = 1 # 簇成员的个数# 簇头选举count = 0 # 簇头个数for i in range(n):if Node[i].selected == 'N' and Node[i].flag != 0:if Node[i].d > d0:alpha = 4 # 能量损失指数else:alpha = 2Eavg = 0 # 系统节点平均能量m = 0# 存活节点个数for j in range(n):if Node[i].flag != 0:Eavg = Eavg + Node[i].powerm = m + 1if m != 0:Eavg = Eavg / n # 计算系统节点平均能量else:breakif Node[i].temp_rand <= (p / (1 - p * (r % round(1 / p)) * (Node[i].power / Eavg) ** (1 / alpha))) and \Node[i].d > Node[i].Rc:Node[i].type = 'C'# 节点类型为簇头Node[i].selected = 'O' # 该节点标记'O'，说明当选过簇头Node[i].CH = -1count = count + 1final_CH.append(i)# 加入簇头节点集合# 广播自成为簇头distanceBroad = (Node[i].Rc ** 2 + Node[i].Rc ** 2) ** 0.5if (distanceBroad > d0):Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Emp * ctrPacketLength * (distanceBroad ** 4))else:Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Efs * ctrPacketLength * distanceBroad ** 2)else:Node[i].type = 'N' # 节点类型为普通# 计算邻居节点集合for i in range(n):cnt = 0for j in range(n):if i != j:dist = ((Node[i].xd-Node[j].xd)**2+(Node[i].yd-Node[j].yd)**2)**0.5if dist < Node[i].Rc:cnt = cnt + 1Node[i].N[cnt] = j# if len(Node[i].N)<cnt:#Node[i].N.append(j)# else:#Node[i].N[cnt-1] = jif j == n:Node[i].Num_N = cnt# 计算前邻节点集for i in range(n):cnt = 0for j in range(Node[i].Num_N):ne = Node[i].N[j]if Node[ne].d < Node[i].d:cnt = cnt + 1Node[i].FN[cnt] = neif j == Node[i].Num_N:Node[i].Num_FN = cnt# 计算前簇头节点集for i in range(count):cnt = 0for j in range(Node[i].Num_FN):fne = Node[final_CH[i]].FN[j]if fne != 0 and Node[fne].d < Node[final_CH[i]].d and Node[fne].CH == -1:cnt = cnt + 1Node[final_CH[i]].CN[cnt] = fneif j == Node[i].Num_FN:Node[final_CH[i]].Num_CN = cnt# 加入簇for i in range(n):if Node[i].type == 'N' and Node[i].power > 0:E = np.zeros(count)for j in range(count):dist = ((Node[i].xd-Node[final_CH[j]].xd)**2+(Node[i].yd-Node[final_CH[j]].yd)**2)**0.5if dist < Node[final_CH[j]].Rc: # 满足条件1E[j] = (Node[final_CH[j]].power-Emin)/Node[final_CH[j]].num_joinif len(E)>0:max_value, max_index = np.max(E),np.argmax(E)else:max_value, max_index = 0,0# 节点发送加入簇的消息if len(final_CH) != 0:dist = ((Node[i].xd - Node[final_CH[max_index]].xd) ** 2 + (Node[i].yd - Node[final_CH[max_index]].yd) ** 2) ** 0.5if dist > Node[final_CH[max_index]].Rc: # 不满足条件1，选择最近的簇头加入Length = np.zeros(count)for j in range(count):Length[j] = ((Node[i].xd - Node[final_CH[j]].xd) ** 2 + (Node[i].yd - Node[final_CH[j]].yd) ** 2) ** 0.5min_value, min_index = np.min(Length), np.argmin(Length)Node[i].CH = final_CH[min_index]# 节点发送加入簇的消息Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Efs * ctrPacketLength * (dist ** 2))# 簇头接收消息Node[final_CH[min_index]].power = Node[final_CH[min_index]].power - Eelec * ctrPacketLengthNode[final_CH[min_index]].num_join = Node[final_CH[min_index]].num_join + 1else:# 节点发送加入簇的消息Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Efs * ctrPacketLength * (dist ** 2))# 簇头接收消息Node[final_CH[max_index]].power = Node[final_CH[max_index]].power - Eelec * ctrPacketLengthNode[final_CH[max_index]].Rc = Rmax * Node[final_CH[max_index]].power / E0Node[i].CH = final_CH[max_index]Node[final_CH[max_index]].num_join = Node[final_CH[max_index]].num_join + 1# 能量模型# 发送数据for i in range(n):if Node[i].flag != 0:if Node[i].type == 'N' and Node[i].CH != 0: # 普通节点dist = ((Node[i].xd - Node[Node[i].CH].xd) ** 2 + (Node[i].yd - Node[Node[i].CH].yd) ** 2) ** 0.5if dist > d0:Node[i].power = Node[i].power - (Eelec * packetLength + Emp * packetLength * (dist ** 4))else:Node[i].power = Node[i].power - (Eelec * packetLength + Efs * packetLength * (dist ** 2))else: # 簇头节点Node[i].power = Node[i].power - (Eelec + ED) * packetLength # 簇头接收数据if Node[i].d <= Node[i].Rc:Node[i].power = Node[i].power - (Eelec * packetLength + Efs * packetLength * (Node[i].d ** 2))else:if Node[i].Num_CN == 0:if Node[i].d > d0:Node[i].power = Node[i].power - (Eelec * packetLength + Emp * packetLength * (Node[i].d ** 4))else:Node[i].power = Node[i].power - (Eelec * packetLength + Efs * packetLength * (Node[i].d ** 2))else:# 选择中继节点dis = np.zeros((Node[i].Num_CN, 1))# 计算前簇头节点距基站的距离for j in range(Node[i].Num_CN):dis[j] = Node[Node[i].CN[j]].dindex = np.argsort(dis)di = dis[index]# 中继转发for j in range(Node[i].Num_CN):Node[i].power = Node[i].power - di[j] / np.sum(di) * (Eelec * packetLength + Emp * packetLength * (di[Node[i].Num_CN + 1 - j] ** 2))for i in range(n):if Node[i].power < Emin:Node[i].flag = 0final_CH = []if f == 0:stop = rmax# load data.mat 节点复位for i in range(n):# Node[i].temp_rand = random.random()# rand为(0,1)的随机数Node[i].type = 'N'# 进行选举簇头前先将所有节点设为普通节点Node[i].selected = 'N' # ’O'：当选过簇头，N：没有Node[i].power = E0# 初始能量Node[i].CH = 0# 保存普通节点的簇头节点，-1代表自己是簇头Node[i].flag = 1 # 1代表存活；0代表死亡Node[i].N = [0 for _ in range(n) ] # 邻居节点集Node[i].Num_N = 0# 邻居节点集个数Node[i].FN = [0 for _ in range(n)] # 前邻节点集Node[i].Num_FN = 0# 前邻节点集个数Node[i].CN = [0 for _ in range(n)] # 前簇头节点集Node[i].Num_CN = 0# 前簇头节点集个数################LEACH##################alive_leach = np.zeros((rmax, 1)) # 每轮存活节点数re_leach = np.zeros((rmax, 1)) # 每轮节点总能量for r in range(rmax):for i in range(n):if Node[i].flag != 0:re_leach[r] = re_leach[r]+Node[i].poweralive_leach[r] = alive_leach[r]+1f = 0if alive_leach[r] == 0:stop = rf = 1breakfor i in range(n):Node[i].type = 'N' # 进行选举簇头前先将所有节点设为普通节点Node[i].selected = 'N'Node[i].temp_rand = random.random() #节点取一个(0,1)的随机值，与p比较for i in range(n):if Node[i].selected == 'N' and Node[i].flag != 0:# if Node[i].type=='N' #只对普通节点进行选举，即已经当选簇头的节点不进行再选举if Node[i].temp_rand <= (p/(1-p*(r%round(1/p)))): # 选取随机数小于等于阈值，则为簇头Node[i].type = 'C' # 节点类型为蔟头Node[i].selected = 'O' # 该节点标记'O'，说明当选过簇头Node[i].CH = -1# 广播自成为簇头distanceBroad = (xm * xm + ym * ym) ** 0.5if distanceBroad > d0:Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Emp * ctrPacketLength * (distanceBroad ** 4))else:Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Efs * ctrPacketLength * (distanceBroad ** 2))else:Node[i].type = 'N' # 节点类型为普通# 判断最近的簇头结点，如何去判断，采用距离矩阵yy = np.zeros((n,n))Length = np.zeros((n,n))for i in range(n):if Node[i].type == 'N' and Node[i].flag != 0:for j in range(n):if Node[j].type == 'C' and Node[j].flag != 0: # 计算普通节点到簇头的距离Length[i, j] = ((Node[i].xd - Node[j].xd) ** 2 + (Node[i].yd - Node[j].yd) ** 2) ** 0.5else:Length[i, j] = float('inf')dist, index = np.min(Length[i, :]), np.argmin(Length[i, :]) # 找到距离簇头最近的簇成员节点# 加入这个簇if Length[i, index] < d0:Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Efs * ctrPacketLength * (Length[i, index] ** 2))else:Node[i].power = Node[i].power - (Eelec * ctrPacketLength + Emp * ctrPacketLength * (Length[i, index] ** 4))Node[i].CH = index# 接收簇头发来的广播的消耗Node[i].power = Node[i].power - Eelec * ctrPacketLength# 对应簇头接收确认加入的消息Node[index].power = Node[index].power - Eelec * ctrPacketLengthyy[i, index] = 1else:Length[i, :] = float('inf')for i in range(n):if Node[i].flag != 0:if Node[i].type == 'C':number = np.sum(yy[:, i])# 统计簇头节点i的成员数量# 簇头接收普通节点发来的数据Node[i].power = Node[i].power-(Eelec+ED)*packetLength# 簇头节点向基站发送数据len = ((Node[i].xd-sink.x)**2+(Node[i].yd-sink.y)**2)**0.5if len < d0:Node[i].power = Node[i].power-((Eelec+ED)*packetLength+Efs*packetLength*(len**2))else:Node[i].power = Node[i].power-((Eelec+ED)*packetLength+Emp*packetLength*len**4)else:# 普通节点向簇头发数据len = Length[i, Node[i].CH]if len < d0:Node[i].power = Node[i].power - (Eelec * packetLength + Efs * packetLength * len ** 2)else:Node[i].power = Node[i].power - (Eelec * packetLength + Emp * packetLength * len ** 4)for i in range(n):if Node[i].power < 0:Node[i].flag = 0if f == 0:stop = rmax## 绘图显示fig4 = plt.figure(dpi=80)plt.plot(range(rmax), alive_ima_leach, c='r',linewidth=2)plt.plot(range(rmax), alive_leach, c='b', linewidth=2)plt.legend(['IMP_LEACH', 'LEACH'])plt.xlabel('轮数')plt.ylabel('存活节点数')fig5 = plt.figure(dpi=80)plt.plot(range(rmax), re_ima_leach, c='r',linewidth=2)plt.plot(range(rmax), re_leach, c='b', linewidth=2)plt.legend(['IMP_LEACH', 'LEACH'])plt.xlabel('轮数')plt.ylabel('系统总能量')fig6 = plt.figure(dpi=80)for r in range(rmax):if alive_ima_leach[r] >= n:a1 = rif alive_leach[r] >= n:a2 = rif alive_ima_leach[r] >= (1-0.1)*n:b1 = rif alive_leach[r] >= (1-0.1)*n:b2 = rif alive_ima_leach[r] >= (1-0.2)*n:c1 = rif alive_leach[r] >= (1-0.2)*n:c2 = rif alive_ima_leach[r] >= (1-0.5)*n:d1 = rif alive_leach[r] >= (1-0.5)*n:d2 = rif alive_ima_leach[r] > 0:e1 = rif alive_leach[r] > 0:e2 = ry=[[a1, a2],[b1, b2],[c1, c2],[d1, d2],[e1, e2]]y=np.array(y)x=np.arange(5)width = 0.36plt.bar(x-width/2,y[:,0],width=width)plt.bar(x+width/2,y[:,1],width=width)plt.xticks(range(5),['0', '10', '20', '50', '100'])plt.legend(['IMP_LEACH', 'LEACH'])plt.xlabel('死亡比例')plt.ylabel('循环轮数')plt.show()