Commit 5aa04bdb authored by Gonzalo Belcredi's avatar Gonzalo Belcredi

Second commit

parent 39f3795c
import numpy as np
import config
class Channel:
def __init__(self, ChannelId, NumberNodes, Capacity, SlopeLin, SlopeSat):
self.Capacity = Capacity
self.SlopeLin = SlopeLin
self.SlopeSat = SlopeSat
self.ChannelId = ChannelId
self.Linear = np.ones(config.SimulationTimeslots)
self.NrUsers = np.zeros(config.SimulationTimeslots)
self.AccessMatrix = np.zeros((NumberNodes,config.SimulationTimeslots))
def channel_access(self,NodeId,Timestamp):
if (Timestamp < config.SimulationTimeslots):
self.AccessMatrix[NodeId,Timestamp] = 1
self.NrUsers[Timestamp]+=1
def update_state(self,Timestamp):
if (Timestamp < config.SimulationTimeslots):
if (self.NrUsers[Timestamp] < (self.Capacity / self.SlopeLin)):
self.Linear[Timestamp] = 1
else:
self.Linear[Timestamp] = 0
\ No newline at end of file
import numpy as np
# SAMPLING
t_start = 0.
t_end = 80. # units of time
step = .01
SimulationTimeslots = int(t_end / step)
time = np.arange(t_start,t_end,step)
# NUMBER NODES AND GATEWAYS
NumberNodes = 5
NumberGateways = NumberNodes / 20
# Problem data.
mu = .003 # Processing Rate (tasks per time step)
lda = .004 # Incoming Rate (tasks per time step)
RxRate = lda
v = 1. # Tasks per Bit
alfa = .1 # Processing Factor
#Channel
Cg = .1 # Gateway Maximum Throughput
M = 15 # Each Gateway can service M users in linear zone
c = Cg * NumberGateways / NumberNodes # Maximum Throughput per user
C = Cg * NumberGateways # Total Maximum Thoughput
k = Cg / M # Linear Slope
k1 = .1 * k # Saturated Slope
# QoS
epsilon = .1
epsilon_energy = .2
beta = .1
delta = 1.
eta = .1
#Energy
energy_model = 1
MinLevel = 0
gamma = NumberNodes * .15
energy_harvesting = 1
energy_tx = 1
energy_pl = 1
if not energy_model:
energy_harvesting = 0
energy_tx = 0
energy_pl = 0
\ No newline at end of file
import numpy as np
from tqdm import tqdm
from time import sleep
import objects
import node
import config
import utils
StorageObj = objects.StorageObj
ChannelObj = objects.ChannelObj
StatsObj = objects.StatsObj
HarvestingVector = np.random.poisson(config.gamma, config.SimulationTimeslots)
NumberNodes = config.NumberNodes
SimulationTimeslots = config.SimulationTimeslots
nodos = []
for i in range(NumberNodes):
np.random.seed(i)
nodos.append(node.IoTNode(config.SimulationTimeslots, config.mu, config.alfa, config.energy_tx, config.energy_pl,
config.energy_harvesting,NodeId=i, ReceivingRate = config.RxRate, HarvestingRate = config.gamma,
BatteryCapacity=10))
ntx = np.zeros(SimulationTimeslots)
npl = np.zeros(SimulationTimeslots)
nidle = np.zeros(SimulationTimeslots)
battery_level = np.ones(SimulationTimeslots)
consumption = np.zeros(SimulationTimeslots)
with tqdm(total=SimulationTimeslots) as pbar:
for timeslot in range(SimulationTimeslots):
StorageObj.Level += HarvestingVector[timeslot]
EnergyStart = StorageObj.Level
StorageObj.update_state()
for nodo in nodos:
ChannelObj.update_state(timeslot)
nodo.run()
battery_level[timeslot] = StorageObj.get_level()
consumption[timeslot] = battery_level[timeslot] - EnergyStart
pbar.update(1)
ntx = StatsObj.Ntx / NumberNodes
npl = StatsObj.Npl / NumberNodes
nidle = 1 - ntx -npl
utils.phase_diagram(npl,ntx,nidle)
\ No newline at end of file
import numpy as np
import config
import math
import objects
import packet
debug = 0
StorageObj = objects.StorageObj
ChannelObj = objects.ChannelObj
StatsObj = objects.StatsObj
class IoTNode:
def __init__(self,SimulationTimeslots, mu, alfa, energy_tx, energy_pl, energy_harvesting,NodeId,ReceivingRate,HarvestingRate, BatteryCapacity):
self.ReceivingRate = ReceivingRate
self.BatteryCapacity = BatteryCapacity
self.SimulationTimeslots = SimulationTimeslots
self.Status = 'Idle'
self.Timestamp = 0
self.NodeId = NodeId
self.ProcessingRate = mu
self.EnergyRadio = energy_tx
self.EnergyProcessing = energy_pl
self.StatusLog = np.zeros((3,SimulationTimeslots))
self.StatusLog[0,:] = np.ones(SimulationTimeslots)
self.EventsLog = np.zeros(SimulationTimeslots)
self.Availability = np.ones(SimulationTimeslots)
self.BitRate = np.zeros(SimulationTimeslots)
self.TxTime = 1
self.alpha = alfa
self.EventCounter = 0
self.events = ()
self.Loss = np.zeros(SimulationTimeslots)
self.update_status(self.Status)
self.CurrentEvent = {}
self.CurrentEventFinished = True
self.NumberEvents = 0
self.Terminate = 0
self.events = self.generate_packet_events()
self.harvesting_rate = HarvestingRate
self.harvesting_amount = energy_harvesting
self.harvesting_vector = np.zeros(SimulationTimeslots)
self.energy_consumption = np.zeros(SimulationTimeslots)
self.run()
def simulate_energy_harvesting(self):
last_timestamp = 0
np.random.seed(self.NodeId)
while last_timestamp <= self.SimulationTimeslots:
last_timestamp += self.nextTime(self.harvesting_rate)
if last_timestamp < self.SimulationTimeslots:
self.harvesting_vector[last_timestamp] += self.harvesting_amount
def run(self):
if self.Timestamp < self.SimulationTimeslots:
if len(self.events)>0:#Si tengo eventos en la lista
if self.NodeId == 1 and debug:
print("Nodo ", self.NodeId, " Timeslot: ", self.Timestamp, " | Status: ", self.Status)
print("Eventos pendientes: ", [ev[1] for ev in self.events] )
if self.CurrentEventFinished: # Si terminé de procesar el evento anterior
self.CurrentEvent = {}
if len(self.events)>0:
while (len(self.events)>0 and self.events[0][1] < self.Timestamp): # Para el caso de eventos repetidos en timeslot, esto se activa, el evento se saca de la lista y se computa una pérdida.
if len(self.events)>0:
event_tuple = self.events.pop(0)
if self.NodeId == 1 and debug:
print("Evento cancelado: ", event_tuple[1])
self.Loss[event_tuple[1]] =+ 1
if len(self.events)>0:
if self.NodeId == 1 and debug:
print("Próximo evento ", self.events[0][1])
if self.events[0][1] == self.Timestamp: # Si el primer evento en la lista está programado para el tiempo actual
event_tuple = self.events.pop(0) # Me quedo con el primer evento de la lista de eventos y lo saco de la lista
if StorageObj.Status == 'Charged':
if np.random.rand() < self.alpha:
self.Status = 'Processing'
#StorageObj.storage_access(self.Timestamp,-self.EnergyProcessing)
self.CurrentEvent = {'Event': event_tuple,'Status': self.Status}
self.CurrentEventFinished = self.packet_pl(self.CurrentEvent['Event'])
else:
self.Status = 'Transmitting'
#StorageObj.storage_access(self.Timestamp,-self.EnergyRadio)
self.CurrentEvent = {'Event': event_tuple, 'Status': self.Status}
self.CurrentEventFinished = self.packet_tx(self.CurrentEvent['Event'])
elif StorageObj.Status == 'Discharged':
self.Loss[self.Timestamp] += 1
self.CurrentEventFinished = 1
self.Status = 'Idle'
elif self.events[0][1] > self.Timestamp:
self.Status = 'Idle'
else:
if self.CurrentEvent['Status'] == 'Processing':
self.Status = 'Processing'
self.CurrentEventFinished = self.packet_pl(self.CurrentEvent['Event'])
else:
self.Status = 'Transmitting'
self.CurrentEventFinished = self.packet_tx(self.CurrentEvent['Event'])
self.update_status(self.Status)
self.update_time()
def nextTime(self,rateParameter):
# https://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/
return int(np.round(-math.log(1.0 - np.random.random_sample()) / rateParameter))
def generate_packet_events(self):
last_timestamp = 0
events_tuple = []
packet_received = packet.Packet(Type = 'Data', Payload ='23.1232')
np.random.seed(self.NodeId)
while last_timestamp <= self.SimulationTimeslots*1.5:
last_timestamp += self.nextTime(self.ReceivingRate)
if last_timestamp < self.SimulationTimeslots:
events_tuple.append(('packet_event',last_timestamp,packet_received))
self.NumberEvents +=1
return events_tuple
def update_status(self,NewStatus):
if (self.Timestamp < self.SimulationTimeslots):
if NewStatus == 'Processing' :
self.StatusLog[0,self.Timestamp] = 0
self.StatusLog[1,self.Timestamp] = 1
self.StatusLog[2,self.Timestamp] = 0
self.Availability[self.Timestamp] = 0
elif NewStatus == 'Transmitting':
self.StatusLog[0,self.Timestamp] = 0
self.StatusLog[1,self.Timestamp] = 0
self.StatusLog[2,self.Timestamp] = 1
self.Availability[self.Timestamp] = 0
def update_time(self):
self.Timestamp +=1
if (self.Timestamp < self.SimulationTimeslots):
self.Availability[self.Timestamp] = 0
def update_battery(self, NetFlow):
if self.Timestamp < self.SimulationTimeslots:
StorageObj.storage_access(self.Timestamp,NetFlow)
def packet_pl(self,event):
TimeStart = event[1]
TimeEnd = TimeStart + int(np.round(1./(self.ProcessingRate)))
if self.Timestamp <= TimeEnd:
finished = 0
if StorageObj.Status == 'Charged' and StatsObj.energy_region(self.Timestamp,config.epsilon_energy,config.gamma):
self.update_status('Processing')
StorageObj.storage_access(self.Timestamp,-self.EnergyProcessing)
StatsObj.add_processor(self.Timestamp)
else:
finished = 1
self.Loss[self.Timestamp] +=1
self.CurrentEventFinished = 1
self.status = 'Idle'
elif self.Timestamp > TimeEnd:
finished = 1
self.status = 'Idle'
return finished
def packet_tx(self,event):
finished = 0
TimeStart = event[1]
if self.Timestamp == TimeStart:
if (ChannelObj.NrUsers[self.Timestamp - 1 ] < (ChannelObj.Capacity / ChannelObj.SlopeLin)):
lineal = True
else:
lineal = False
if ChannelObj.NrUsers[self.Timestamp - 1] > 0:
contendants = ChannelObj.NrUsers[self.Timestamp - 1]
else:
contendants = 1
if lineal:
TxTime = 1 / (config.v * (1.0*ChannelObj.SlopeLin))
else:
TxTime = 1 / (config.v * ((ChannelObj.Capacity/contendants) *(1 + ChannelObj.SlopeSat/ChannelObj.SlopeLin)-ChannelObj.SlopeSat))
self.TxTime = TxTime
else:
TxTime = self.TxTime
TimeEnd = TimeStart + int(np.round(TxTime))
if self.Timestamp <= TimeEnd:
if StorageObj.Status == 'Charged' and StatsObj.energy_region(self.Timestamp,config.epsilon_energy,config.gamma):
self.update_status('Transmitting')
self.BitRate[self.Timestamp] = 1 / (config.v * TxTime)
ChannelObj.channel_access(self.NodeId,self.Timestamp)
StorageObj.storage_access(self.Timestamp,-self.EnergyRadio)
StatsObj.add_transmitter(self.Timestamp)
else:
finished = 1
self.Loss[self.Timestamp] +=1
self.CurrentEventFinished = 1
self.status = 'Idle'
elif self.Timestamp > TimeEnd:
finished = 1
self.status = 'Idle'
return finished
\ No newline at end of file
import channel
import storage
import stats
import config
import numpy as np
ChannelObj = channel.Channel(ChannelId=1, NumberNodes=config.NumberNodes,
Capacity=config.C, SlopeLin=config.k, SlopeSat=config.k1)
StorageObj = storage.Storage(100, config.gamma, config.MinLevel)
StatsObj = stats.Stats()
HarvestingVector = np.random.poisson(config.gamma, config.SimulationTimeslots)
\ No newline at end of file
class Packet:
def __init__(self, Type = 'Data', Payload ='23.1232'):
self.Type = Type
self.Length = len(Payload)
self.Payload = Payload
\ No newline at end of file
import config
import numpy as np
class Stats:
def __init__(self):
self.Ntx = np.zeros(config.SimulationTimeslots)
self.Npl = np.zeros(config.SimulationTimeslots)
def add_transmitter(self, timestamp):
self.Ntx[timestamp] += 1
def add_processor(self, timestamp):
self.Npl[timestamp] += 1
def energy_region(self,timestamp,epsilon,gamma):
if self.Ntx[timestamp] + self.Npl[timestamp] < config.gamma / (1 - config.epsilon_energy):
value = 1
else:
value = 0
return value
\ No newline at end of file
import random
import config
import numpy as np
class Storage:
def __init__(self, Capacity, HarvestingRate, MinLevel):
self.Capacity = Capacity
self.HarvestingRate = HarvestingRate
self.HarvestingAmount = 1
self.Level = 10
self.Status = 'Charged'
self.MinLevel = MinLevel
self.Leakage = 0
self.HarvestingVector = np.zeros(config.SimulationTimeslots)
def nextTime(self,rateParameter):
value = random.expovariate(rateParameter)
return value
def storage_access(self,Timestamp,NetFlow):
if (Timestamp < config.SimulationTimeslots):
if self.Level + NetFlow >= self.MinLevel:
self.Level += NetFlow
self.update_state()
def get_level(self):
return self.Level
def get_status(self):
return self.Status
def update_state(self):
if self.Level > self.MinLevel:
self.Status = 'Charged'
else:
self.Status = 'Discharged'
\ No newline at end of file
This diff is collapsed.
import config
import numpy as np
import math
import matplotlib.lines as lines
from matplotlib import pyplot as plt
from matplotlib.patches import Circle, Wedge, Polygon
from scipy.integrate import odeint
from matplotlib.collections import PatchCollection
time = config.time
t_start = config.t_start
t_end = config.t_end
alfa = config.alfa
mu = config.mu
c = config.c
lda= config.lda
k= config.k
k1= config.k1
v= config.v
step = config.step
gamma = config.gamma
NumberNodes = config.NumberNodes
epsilon_energy = config.epsilon_energy
SimulationTimeslots = config.SimulationTimeslots
def ode(values, t, params):
X, Y = values # unpack current values
alfa,lda,mu,c,k,k1,v,restrictions = params # unpack parameters
if Y >= c/k:
derivs = [alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) -(c*(1+k1/k) - k1 *Y)*v]
#if (c*(1+k1/k) - k1 *Y) < 0: #Capacidad del canal es positiva, de lo contrario hago que C(n_t)=0
# derivs = [alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y)]
if X+Y >= 1 and (restrictions==True):
proy_x = np.dot(np.array([alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) -(c*(1+k1/k) - k1 *Y)*v]),W[:,0])
proy_y = np.dot(np.array([alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) -(c*(1+k1/k) - k1 *Y)*v]),W[:,1])
derivs = [-X*mu,-(c*(1+k1/k) - k1 *Y)*v]
else:
derivs = [alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) - k*Y*v]
return derivs
def solve_ode(t_start,t_end,y0,alfa,lda,mu,c,k,k1,v,restrictions = False):
t = np.arange(t_start,int(t_end/step),1)
params = [alfa,lda,mu,c,k,k1,v,restrictions]
y = odeint(ode, y0, t, args=(params,))
n_p = y[:,0]
n_t = y[:,1]
return n_p,n_t,t
def phase_diagram(npl,ntx,nidle):
plt.figure(figsize = (12,10))
plt.plot(time,ntx,color='r',label='Transmitting',alpha=.5)
plt.plot(time,npl,color='g',label='Processing',alpha=.5)
plt.plot(time,nidle,color='b',label='Idle',alpha=.5)
plt.plot(time,ntx+npl+nidle,color='k',label='Total',alpha=.5)
plt.xlabel('Time (s)')
plt.legend()
plt.show()
# Initial values
X0 = 0.0 # initial n_p
Y0 = 0.0 # initial n_t
y0 = [X0,Y0]
#y = solve_lin(t_start,t_end,np.array([[X0],[Y0]]),alfa,lda,mu,c,k,k1,v)
npl_t,ntx_t,t = solve_ode(t_start,t_end,y0,alfa,lda,mu,c,k,k1,v)
''' ODE's SET '''
patches = []
polygon = Polygon([[0,c/k],[0,1],[1-c/k,c/k]], True)
patches.append(polygon)
polygon = Polygon([[0,0],[0,c/k],[1-c/k,c/k],[1,0]], True)
patches.append(polygon)
nt_sat_eq = (c*(1+k1/k)*v-mu*(1-alfa)*lda/(mu+alfa*lda)) / (k1*v-mu*(1-alfa)*lda/(mu+alfa*lda))
np_sat_eq = (alfa*lda/(mu+alfa*lda))*(k1*v-c*(1+k1/k)*v)/ (k1*v+(1-alfa)*lda*(alfa*lda/(mu+alfa*lda)-1))
nt_lin_eq = 1/(1+k*v/(1-alfa)*(1/lda+alfa/mu))
np_lin_eq = 1/(1+mu/alfa*(1/lda+(1-alfa)/(k*v)))
#alfa_tilde = (c*(1/k+v/lda)-1)/(c/k-c*v/mu-1)
#alfa_vector = [alfa_tilde]
#alfa = .6*(mu**2+lda**2+(k*v)**2+2*lda*(k*v-mu)-2*k*v*mu)/(4*lda*(k*v-mu))
#alfa = 1.2*alfa_tilde
''' GRID '''
nx, ny = .05, .05
x = np.arange(0, 1, nx)
y = np.arange(0, 1, ny)
X, Y = np.meshgrid(x, y)
dx = np.zeros((len(y),len(x)))
dy = np.zeros((len(y),len(x)))
dx = alfa*lda*(1-X-Y)-X*mu
for i in range(len(y)):
for j in range(len(x)) :
if Y[i][j] > c/k:
dy[i][j] = (1-alfa)*lda*(1-X[i][j]-Y[i][j])-(c*(1+k1/k)-k1*Y[i][j])*v
else:
dy[i][j] = (1-alfa)*lda*(1-X[i][j]-Y[i][j])-k*v*Y[i][j]
fig, ax = plt.subplots(figsize=(14,14))
p = PatchCollection(patches, alpha=0.1)
p.set_color(['red','blue'])
ax.add_collection(p)
plt.quiver(X, Y, dx, dy,np.hypot(dx, dy), headlength=9)
plt.xlim(-0, 1)
plt.ylim(-0, 1)
plt.xlabel('n_p')
plt.ylabel('n_t')
plt.plot(npl,ntx,'-',markersize=4,color='r',alpha=.5,label='Simulation Trajectory')
plt.plot(npl_t,ntx_t,linewidth=3,c='blue',label = 'Model Trajectory')
plt.plot(np.linspace(0,1,100),gamma / (NumberNodes*(1 - epsilon_energy)) - np.linspace(0,1,100),'-',label = 'Energetic Boundary')
plt.scatter(np_lin_eq,nt_lin_eq,marker='x',c='green',s=150,label='Equilibrium Lin Point (Model)')
plt.scatter(np_sat_eq,nt_sat_eq,marker='+',c='orange',s=150,label='Equilibrium Sat Point (Model)')
plt.legend()
plt.savefig('fig/Evolucion del sistema.pdf')
plt.show()
nidle_t = 1 - npl_t - ntx_t
ntotal = npl + ntx + nidle
ntotal_t = npl_t + ntx_t + nidle_t
if ntx_t[-1]<c/k:
np_eq = np_lin_eq
nt_eq = nt_lin_eq
else:
np_eq = np_sat_eq
nt_eq = nt_sat_eq
plt.figure(figsize = (14,10))
plt.grid('on')
plt.axhline(np_eq,linestyle='dotted',color='green',label='Equilibrium Processing (Model)')
plt.axhline(nt_eq,linestyle='dotted',color='red',label='Equilibrium Transmitting (Model)')
plt.plot(np.linspace(0,t_end,SimulationTimeslots),ntx,color='r',label='Transmitting',alpha=.2)
plt.plot(time,ntx_t,linestyle='dashed',linewidth=3,color='r',label='Transmitting (Model)')
plt.plot(np.linspace(0,t_end,SimulationTimeslots),npl,color='g',label='Processing',alpha=.2)
plt.plot(time,npl_t,linestyle='dashed',linewidth=3,color='g',label='Processing (Model)')
plt.plot(np.linspace(0,t_end,SimulationTimeslots),nidle,color='b',label='Idle',alpha=.2)
plt.plot(time,nidle_t,linestyle='dashed',linewidth=3,color='b',label='Idle (Model)')
plt.ylabel('Population Density')
plt.xlabel('Time (s)')
plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc='lower left',
ncol=2, mode="expand", borderaxespad=0.)
plt.savefig('Population Evolution.pdf')
plt.show()
\ No newline at end of file
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment