Python实现Lora节点拓扑自动规划

NodeClass.py:

# 小记：
# 先不考虑中继节点在classB模式下发Beacon的能耗
# 能量效率的计算公式是否有问题？只考虑发包的能耗和时间
# energyConsumption需统一成自己在lifetime中的能量消耗推导，不再用论文中笼统的折线图值
# 注意区分self.sendTo.index和一些函数中的node.index的区别，前者使用的场景在确定节点拓扑之后，后者使用在确定节点拓扑之前的一些遍历情况

import math
import numpy as np


# python中的public、private和protected是通过变量名前的下划线标识的，保护类型一个下划线，私有类型两个下划线
class Node:
	def __init__(self, name, locX, locY, nodeNum, _payload, _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE, _dataGen, _cycleT,
				 _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, isGateway=0, whichSet=0, index=0): # 默认缺省值isGateway=0
		self.name = name
		self.payload = _payload
		self.bandW = _bandW
		self.c = _c
		self.pi = _pi
		self.beta = _beta
		self.whiteNoise = _whiteNoise
		self.dutyCycle = _dutyCycle
		self.batteryE = _batteryE
		# 每确定要多为另一个节点进行中继，dataGen需要更新添加帮忙转发的数据量
		self.dataGen = _dataGen
		self.cycleT = _cycleT
		self.transRadioE = _transRadioE
		self.recvRadioE = _recvRadioE
		self.spareRadioE = _spareRadioE
		self.offRadioE = _offRadioE
		self.onMcuE = _onMcuE
		self.offMcuE = _offMcuE
		# 注意和dataGen的区别，这个变量是节点发送给节点/网关的总数据量，带有中继数据
		self.dataSend = _dataGen
		self.lifeT = 0.0
		self.cycleE = 0.0
		self.symbleT = 0.0
		self.pktT = 0.0
		self.transT = 0.0
		self.recvT = 0.0
		self.spareT = 0.0
		self.airtime = 0.0
		self.goodput = 0.0
		# 用于MST中判断节点在哪个集合中,初始集合为0，MST集合为1
		self.whichSet = whichSet
		# 标识一个“节点”是否是网关，输入中第一个“点”为网关，坐标[0, 0]
		self.isGateway = isGateway
		self.locX = locX
		self.locY = locY
		# 给节点一个索引标识
		self.index = index
		self.relayNum = 0
		# 存放了此节点通往其他节点时使用的参数值
		self.spreadFactor = [7 for i in range(nodeNum)]
		self.transPower = [13 for i in range(nodeNum)]
		self.channel = [1 for i in range(nodeNum)]
		self.codingRate = [5 for i in range(nodeNum)]
		self.pdr = 0.0
		# 标识此节点的数据直接发送给哪个节点，为Node类型
		self.sendTo = self
		# 标识此节点帮助哪些节点进行数据转发
		self.relayFrom = []
		# 标识此节点工作在哪个class下，class会影响cycleE
		self.workClass = 'A'

	# self._prop = []

	def distNode(self, node):
		dist = math.pow((math.pow((self.locX - node.locX), 2) + math.pow((self.locY - node.locY), 2)), 0.5)
		# print(self.locX)
		# print(self.locY)
		# print(node.locX)
		# print(node.locY)
		# print('Distance of ' + self.name + ' to ' + node.name + ' = ' + str(dist))
		return dist

	# 8个信道，从902.1开始每次加0.2
	def freq(self, node):
		frequency = 902.1 + 0.2 * self.channel[node.index]
		frequency *= math.pow(10, 6)
		# print('Frequency of channel ' + str(self.channel[node.index]) + ' = ' + str(frequency))
		return frequency

	# Rayleigh fading channel
	def g(self):
		reyleigh = np.random.exponential(1.0, size=None)
		# print(self.name + ' return rayleigh fading channel g of: ' + str(reyleigh))
		return reyleigh


	# 此为论文中的折线图值观测值实现的能耗函数，现替换成计算式方式实现
	# energy consumption, unit is mW
	def enerConsump(self, node):
		energyPacket = [[40, 50], [50, 67], [65, 85], [95, 105], [100, 125], [110, 135]]
		clsP = 0
		if self.transPower[node.index] <= 15:
			cls_p = 0
		else:
			cls_p = 1
		clsS = self.spreadFactor[node.index] - 7
		# rst = energy_packet[cls_s][cls_p] * 0.001
		rst = energyPacket[clsS][clsP]
		# print("ec(line 70, unit is mW):", rst)
		name = self.name
		sf = str(self.spreadFactor[node.index])
		tp = str(self.transPower[node.index])  # 只是为了下一行不超长...
		# print('Energy Consumption of ' + name + ' with sf=' + sf + ', tp=' + tp + ' is: ' + str(rst))
		return rst

	# 公式计算方式求energy consumption
	# def enerConsump(self):



	def th(self, si):
		if si == 7:
			return -6
		elif si == 8:
			return -9
		elif si == 9:
			return -12
		elif si == 10:
			return -15
		elif si == 11:
			return -17.5
		elif si == 12:
			return -20

	def ss(self, si):
		if si == 7:
			return -123
		elif si == 8:
			return -126
		elif si == 9:
			return -129
		elif si == 10:
			return -132
		elif si == 11:
			return -134.5
		elif si == 12:
			return -137

	def tomWatt(self, dBm):
		mWatt = math.pow(10, (dBm / 10))
		# print(self.name + ' dBm to mWatt is: ' + str(mWatt))
		return mWatt

	# 如果有relay，EE计算为payload * pdr1 * pdr2 / (energyInTime1 + energyInTime2)
	def energyEffi(self, node):
		tempVar = (8 * self.payload - 4 * self.spreadFactor[node.index] + 28 + 16) / (4 * self.spreadFactor[node.index])
		time = (20.25 + max(math.ceil(tempVar) * self.channel[node.index], 0)) * 2 ** self.spreadFactor[node.index] / self.bandW
		energyInTime = time * self.enerConsump(node)
		energyEfficiency = self.payload * self.selfToNodePdr(node) / energyInTime
		# print('EnergyEfficiency of ' + self.name + ' is: ' + str(energyEfficiency))
		return energyEfficiency

	def energyIntime(self, node):
		tempVar = (8 * self.payload - 4 * self.spreadFactor[node.index] + 28 + 16) / (4 * self.spreadFactor[node.index])
		time = (20.25 + max(math.ceil(tempVar) * self.channel[node.index], 0)) * 2 ** self.spreadFactor[
			node.index] / self.bandW
		energyInTime = time * self.enerConsump(node)
		return energyInTime

	def energyEffiRelay(self, node1, node2):
		tempVar1 = (8 * self.payload - 4 * self.spreadFactor[node1.index] + 28 + 16) / (4 * self.spreadFactor[node1.index])
		time1 = (20.25 + max(math.ceil(tempVar1) * self.channel[node1.index], 0)) * 2 ** self.spreadFactor[node1.index] / self.bandW
		energyInTime1 = time1 * self.enerConsump(node1)
		tempVar2 = (8 * node1.payload - 4 * node1.spreadFactor[node2.index] + 28 + 16) / (
					4 * node1.spreadFactor[node2.index])
		time2 = (20.25 + max(math.ceil(tempVar2) * node1.channel[node2.index], 0)) * 2 ** node1.spreadFactor[
			node2.index] / node1.bandW
		energyInTime2 = time2 * node1.enerConsump(node2)
		energyEfficiency = self.payload * self.selfToNodePdr(node1) * node1.selfToNodePdr(node2) / (energyInTime1 + energyInTime2)
		return energyEfficiency

	def relayPdr(self):
		if(self.isGateway == 1):
			return 1
		else:
			temp = self.sendTo.relayPdr() * self.selfToNodePdr(self.sendTo)
		return temp

	def relayEnergyInTime(self):
		if(self.isGateway == 1):
			return 0
		else:
			temp = self.sendTo.relayEnergyInTime() + self.energyIntime(self.sendTo)
		return temp

	# 将设备重叠数考虑进来，重叠、channel会降低pdr
	def selfToNodePdr(self, node):
		# tempVal1是路径损失函数
		tempVar1 = (self.c / (4 * self.pi * self.freq(node) * self.distNode(node))) ** self.beta
		# print('tempVar1:', tempVar1)
		# whiteNoise对pdr的影响大
		tempVar2 = (self.tomWatt(self.th(self.spreadFactor[node.index])) * self.whiteNoise + self.tomWatt(
			self.ss(self.spreadFactor[node.index])))
		# print('tempVar2:', tempVar2)
		# 这地方还是有问题！！！！为什么是/1000/1000，对比注释行
		# pdr = math.exp(-math.pow(tempVar2 / (self.tomWatt(self.transPower) * tempVar1), 0.5))
		pdr = math.exp(-math.pow(tempVar2 / (self.transPower[node.index] * tempVar1), 0.5) / 1000 / 1000)
		# print('Pdr of ' + self.name + ' to ' + node.name + ' is: ' + str(pdr))
		return pdr

	# 主文件里先确定sendTo谁，再更新transT值
	def updateTransT(self):
		self.symbleT = math.pow(2, self.spreadFactor[self.sendTo.index]) / self.bandW
		tempVal1 = 8 * self.payload - 4 * self.spreadFactor[self.sendTo.index] + 28 + 16
		tempVal2 = 20.25 + max(math.ceil(tempVal1 / (4 * self.spreadFactor[self.sendTo.index])) * self.codingRate[self.sendTo.index], 0)
		self.pktT = self.symbleT * tempVal2
		# 计算transT
		# 此节点将数据（包括自己的和中转的）发给其将要发至的节点时的pdr
		pdr = self.selfToNodePdr(self.sendTo)
		if pdr != 0:
			if self.relayFrom:
				for child in self.relayFrom:
					# 根据中继了哪些节点来更新dataSend的数据量，dataSend是此节点总共发送的数据量
					# 若此节点非中继节点，则dataSend == dataGen
					self.dataSend += child.dataGen
			self.transT = self.pktT * self.cycleT / self.dataSend / pdr
		else:
			self.transT = 10000
		return self.transT

	# 确定在self.relayFrom中添加完所有的远端节点后再进行更新
	def updateRecvT(self):
		# 计算recvT，此时间包括自身发送时的class A模式的部分和中继时class B模式的recv time
		# 没考虑作为发送端时的recvT，因为发送时工作在classA，只在发送数据的时候顺便开一个窗口接收数据，暂时不知道窗口总时长		if self.relayFrom:
		if self.relayFrom:
			for child in self.relayFrom:
				self.recvT += child.transT
		print('Update transT of ' + self.name + ' to: ' + str(self.recvT))
		return self.recvT

	def lifetime(self):
		# 设置transRadioE级别（87和120）
		if self.transPower[self.sendTo.index] < 15:
			self.transRadioE = 3.3 * 87e-3
		self.spareT = self.dutyCycle * self.cycleT - self.transT - self.recvT
		radioOnE = self.transT * self.transRadioE + self.recvT * self.recvRadioE + self.spareT * self.spareRadioE
		radioOffE = (1 - self.dutyCycle) * self.cycleT * self.offRadioE
		mcuE = self.cycleT * self.dutyCycle * self.onMcuE + self.cycleT * (1 - self.dutyCycle) * self.offMcuE
		self.cycleE = radioOnE + radioOffE + mcuE
		self.lifeT = self.cycleT * self.batteryE / self.cycleE / 3600 / 24
		print('Lifetime of ' + self.name + ' is: ' + str(self.lifeT))
		return self.lifeT

	def goodput(self):
		# 注意：goodput中计算的数据量用的是dataGen，而不是dataSend
		goodput = self.cycleT / self.dataGen * self.payload / self.transT
		print('Goodput of ' + self.name + ' is: ' + goodput)
		return goodput

TopoWithMST.py:

from MST.Node.NodeClass import Node
import random
import matplotlib.pyplot as plt
import numpy as np
import math
import scipy as sy
from scipy import stats
import xlsxwriter

# 添加node要修改三个地方，node初始化，nodes列表添加该node，nodeNum值

#  常量，全局变量均以“_”开头命名
# _nodeNum = 6  # device number
_payload = 4  # payload size of a packet, need to be modified
_bandW = 125000
_c = 3 * 10 ** 8
_pi = 3.1415926
_beta = 3  # path loss exponent, maybe a list, need to be init
_whiteNoise = random.random()   # gauss white noise
_dutyCycle = 0.02
_batteryE = 3.7 * 2 * 3600
_dataGen = 2000

# Radio耗电常量
_cycleT = 3600 * 12
_transRadioE = 3.3 * 120e-3  # 29, 87, 120
_recvRadioE = 3.3 * 11.5e-3
_spareRadioE = 3.3 * 1.6e-3
_offRadioE = 3.3 * 1.5e-6

# Known and Unknown... P10  MCU耗电情况
_onMcuE = 23.48e-4
_offMcuE = 174.65e-7

sfRange = [7,8,9,10,11,12]
tpRange = [13,14,15,16,17,18,19,20]
chRange = [1,2,3,4,5,6,7,8]
crRange = [5,6,7,8]

nodeNum = 18


# _node0是网关，所以isGateway=1, whichSet=1
_node0 = Node('_node0', 0, 0, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE, _dataGen,
              _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, isGateway=1, whichSet=1)
_node1 = Node('_node1', -2800, -500, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=1)
_node2 = Node('_node2', 2500, 2700, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=2)
_node3 = Node('_node3', 1000, 800, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=3)
_node4 = Node('_node4', 2500, 0, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=4)
_node5 = Node('_node5', 1000, 0, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=5)
_node6 = Node('_node6', 1500, 2000, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=6)
_node7 = Node('_node7', -800, 700, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=7)
_node8 = Node('_node8', -600, -400, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=8)
_node9 = Node('_node9', -1000, -1200, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
              _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=9)
_node10 = Node('_node10', -1800, -2000, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=10)
_node11 = Node('_node11', -2000, 1800, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=11)
_node12 = Node('_node12', 1500, -2000, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=12)
_node13 = Node('_node13', 500, -500, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=13)
_node14 = Node('_node14', 1200, -1200, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=14)
_node15 = Node('_node15', 2000, -2500, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=15)
_node16 = Node('_node16', -2500, -2000, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=16)
_node17 = Node('_node17', 2500, 1000, nodeNum, _payload,  _bandW, _c, _pi, _beta, _whiteNoise, _dutyCycle, _batteryE,
               _dataGen, _cycleT, _transRadioE, _recvRadioE, _spareRadioE, _offRadioE, _onMcuE, _offMcuE, index=17)


# 先求所有的边权（取所有sf，tp，cr选择中代价最小的）存储，再用MST连线（每次连线需根据isGateway确定代价是否叠加）划分集合，直到原本集合为空
indexes= [i for i in range(nodeNum)]
# print(indexes)
nodes = [_node0, _node1, _node2, _node3, _node4, _node5, _node6, _node7, _node8, _node9, _node10, _node11, _node12,
         _node13, _node14, _node15, _node16, _node17]
# 初始化代价矩阵
costMatrix = [[0.0 for i in range(nodeNum)] for j in range(nodeNum)]
# 注意python中不可用 energyEffiMatrix = costMatrix形式给矩阵赋值
# 否则energyEffiMatrix的值会随costMatrix的变化而变化

energyEffiMatrix = [[0.0 for i in range(nodeNum)] for p in range(nodeNum)]


# 计算每条边的代价，存入代价矩阵（不过记得后面MST时代价会根据是否relay确定是否进行叠加）
# 求每条边的代价实际上是先做预处理求最小代价，即遍历了可选的sf、tp、cr、
# 需要考虑sf和ch相同导致的重叠,暂未添加
# zip的使用
for rowIndex, node in zip(indexes, nodes):
    for colIndex, coloum in zip(indexes, nodes):
        if node is not coloum:
            # 针对某一条边求cheapest代价
            maxEE = 0
            setSpreadFactor = 7
            setTransPower = 13
            setChannel = 1
            setCodingRate =5
            for sf in sfRange:
                for tp in tpRange:
                    for ch in chRange:
                        for cr in crRange:
                            node.spreadFactor[coloum.index] = sf
                            node.transPower[coloum.index] = tp
                            node.channel[coloum.index] = ch
                            node.codingRate[coloum.index] = cr
                            tempEE = node.energyEffi(coloum)
                            if(maxEE < tempEE):
                                maxEE = tempEE
                                setSpreadFactor = sf
                                setTransPower = tp
                                setChannel = ch
                                setCodingRate = cr
            node.spreadFactor[coloum.index] = setSpreadFactor
            node.transPower[coloum.index] = setTransPower
            node.channel[coloum.index] = setChannel
            node.codingRate[coloum.index] = setCodingRate
            costMatrix[rowIndex][colIndex] = maxEE
        else:
            costMatrix[rowIndex][colIndex] = 0.1

for row in range(nodeNum):
    for col in range(nodeNum):
        energyEffiMatrix[row][col] = costMatrix[row][col]

for n in range(0, nodeNum):
    print(costMatrix[n])


# 代价是能量效率的倒数
for row in range(nodeNum):
    for col in range(nodeNum):
        costMatrix[row][col] = 1 / costMatrix[row][col]

print('cost statistic:')
for n in range(0, nodeNum):
    print(costMatrix[n])
# print(_node1.spreadFactor)
# print(_node1.transPower)
# print(_node1.channel)
# print(_node1.codingRate)
# 此部分仍有问题，需核对pdr和代价（能量效率EE）计算中各个参数对energy efficiency的影响

toWhichIndexes= [i for i in range(nodeNum)]
# u集合中的节点数
uNodeNum = 1
while(uNodeNum != nodeNum):
    minCost = 9999
    addRelay = nodes[0]
    ifRelay = False
    nodeToUSet = nodes[0]
    backTrack = nodes[0]
    for nodeIndex in range(0, nodeNum):
        # minCost = 9999
        # toWhichIndex为当前查看的v集合中的节点的index，可通过nodes列表反索引到相应的node
        # thisLineCost为单条边的cost，我们需要的是实际的cost（即带relay）
        # relay是在MST的过程中产生的
        # 网关初始应在u集合，对应的whichSet=1

        # addRelay = nodes[0]
        # ifRelay = False
        # nodeToUSet = nodes[0]
        for toWhichIndex, thisLineCost in zip(toWhichIndexes, costMatrix[nodeIndex]):
            # 一个节点属于u集合，一个节点属于v集合，才进行处理，否则不符合MST要求，进行下轮循环
            # addRelay = nodes[0]
            # ifRelay = False
            # nodeToUSet = nodes[0]
            if (nodes[nodeIndex].whichSet == 1 and nodes[toWhichIndex].whichSet == 0):
            # if ((nodes[nodeIndex].whichSet == 1 and nodes[toWhichIndex].whichSet == 0) or
            #         (nodes[nodeIndex].whichSet == 0 and nodes[toWhichIndex].whichSet == 1)):
                # 其中有一个是网关的话，说明与网关直接相连，thisLineCost即为realCost
                if (nodes[nodeIndex].isGateway == 1 or nodes[toWhichIndex].isGateway == 1):
                    realCost = thisLineCost
                    if realCost < minCost:
                        print('node:', nodes[nodeIndex].name, ' nodeTo:', nodes[toWhichIndex].name)
                        backTrack = nodes[nodeIndex]
                        ifRelay = False
                        minCost = realCost
                        nodeToUSet = nodes[toWhichIndex]

                # 所检索的边不是与网关直接相连，需要中继
                else:
                    # print(nodes[nodeIndex].name, nodes[toWhichIndex].name, thisLineCost, costMatrix[toWhichIndex][0], minCost)
                    # 需要中继的，先把代价求倒数转换成EE，两段EE相加之后再求倒数得到整段的代价
                    # [0]代表中继节点到网关的代价，因为网关为矩阵的第一列

                    # cost计算方式1
                    # realCost = thisLineCost + costMatrix[nodeIndex][0]

                    # cost计算方式2，见NodeClass中EE计算函数注释
                    # 这里energyEffiiRelay的实现实际应该是有递归，比如如果有两次relay，则
                    # payload * pdr1 * pdr2 * pdr3 / (energyInTime1 + energyInTime2 + energyInTime3)，怎么实现?
                    # 分别实现一个pdr的递归和一个energyInTime的递归
                    energyEffiRelay = nodes[toWhichIndex].payload * nodes[toWhichIndex].selfToNodePdr(nodes[nodeIndex]) * nodes[nodeIndex].relayPdr() / (nodes[toWhichIndex].energyIntime(nodes[nodeIndex]) + nodes[nodeIndex].relayEnergyInTime())
                    realCost = 1 / energyEffiRelay
                    if realCost < minCost:
                        # 如果12-0代价大于12-14-0，修改cost矩阵，将costMatrix[12][0]修改为12-14-0（经14中继）的代价
                        costMatrix[toWhichIndex][0] = realCost
                        # print(minCost)
                        print('node:', nodes[nodeIndex].name, ' nodeTo:', nodes[toWhichIndex].name)
                        print(1 / thisLineCost)
                        print(1 / costMatrix[nodeIndex][0])
                        print(realCost)
                        backTrack = nodes[nodeIndex]
                        addRelay = nodes[toWhichIndex]  # node类型
                        ifRelay = True
                        minCost = realCost
                        nodeToUSet = nodes[toWhichIndex]

        # 如果有中继的拓扑选择，为中继点修改relayNum, relayFrom信息
    if (ifRelay == True and addRelay):
        # addRelay.relayNum += 1
        # addRelay.relayFrom.append(backTrack.name)
        backTrack.relayNum += 1
        backTrack.relayFrom.append(addRelay.name)
                # 把新选择的节点添加到U集合
    nodeToUSet.whichSet = 1
                # 注意生成树的方向
                # 到最后应该_node0(网关)的sendTo没什么用
    nodeToUSet.sendTo = backTrack
    uNodeNum += 1
    print('one loop.')

for n in range(0, nodeNum):
    print(energyEffiMatrix[n])
    print(costMatrix[n])

# 测试结果
for n in range(0, nodeNum):
    print(costMatrix[n])

# 检查拓扑
print('Topo statistic:')
for n in range(0, nodeNum):
    print(nodes[n].name + ' send to ' + nodes[n].sendTo.name)


print('nodes belong to which set:')
# 检查集合
setList = []
for n in range(0, nodeNum):
    # print(nodes[n].whichSet)
    setList.append(nodes[n].whichSet)
print(setList)

# 检查每个节点的relayNum
# relayNum技术规则还有问题
print(' relay number statistic: ')
for n in range(0, nodeNum):
    print(nodes[n].name + ' relay number: ' + str(nodes[n].relayNum))

print('params choise:')
for n in range(0, nodeNum):
    # 现有结果表明优先增大tp，增到最大才开始增加sf，即增大sf带来的代价远超增大tp
    print(nodes[n].name + ' sf: ' + str(nodes[n].spreadFactor) + ' tp:' + str(nodes[n].transPower) + ' cr:' + str(nodes[n].codingRate))

print('relay:')
for n in range(0, nodeNum):
    print(nodes[n].name , 'relay:', nodes[n].relayFrom)

# 检查EE
print('EE origin:')
for n in range(1, nodeNum):
    print(nodes[n].name, ' EE origin:', energyEffiMatrix[n][0])
print('EE final:')
for n in range(1, nodeNum):
    print(nodes[n].name, ' EE final:', 1 / costMatrix[n][0])


# 绘图
plt.subplot(121)
for n in range(0, nodeNum):
    thisNodeLoc = [nodes[n].locX, nodes[n].locY]
    sendToNodeLoc = [nodes[n].sendTo.locX, nodes[n].sendTo.locY]
    plt.scatter(thisNodeLoc[0], thisNodeLoc[1], s=10, color='b')
    plt.plot([thisNodeLoc[0], sendToNodeLoc[0]], [thisNodeLoc[1], sendToNodeLoc[1]])
    plt.annotate(thisNodeLoc, xy=(thisNodeLoc[0], thisNodeLoc[1]), xytext=(thisNodeLoc[0]+1, thisNodeLoc[1]+1))
plt.xlim(-2800,2800)
plt.ylim(-2800,2800)


# 绘EE柱状图
# subplot绘制子图
plt.subplot(122)
nodeNameList = []
nodeEEOriginList = []
nodeEEFinalList = []
for n in range(0, nodeNum):
    nodeNameList.append(str(n))
    nodeEEOriginList.append(energyEffiMatrix[n][0])
    nodeEEFinalList.append(1 / costMatrix[n][0])
total_width, n = 0.8, 2
width = total_width / n
plt.bar(indexes, nodeEEOriginList, width=width, label='1', fc='b')
for i in range(len(indexes)):
    indexes[i] += width
plt.bar(indexes, nodeEEFinalList, width=width, label='2', tick_label=nodeNameList, fc='g')
plt.legend()


plt.show()