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

Power Grid Storage

parent 1ffb2401
__pycache__
fig
temp
test.py
main_test.py
N=*
test*
Vari*
maint_test*
\ No newline at end of file
No preview for this file type
import numpy as np
import sys
# SAMPLING
t_start = 0.
......@@ -9,7 +11,7 @@ SimulationTimeslots = int(t_end / step)
time = np.arange(t_start,t_end,step)
# NUMBER NODES AND GATEWAYS
NumberNodes = 1000
NumberNodes = 2000
NumberGateways = NumberNodes / 20
# Problem data.
......@@ -17,7 +19,7 @@ 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
alfa = .5 # Processing Factor
#Channel
Cg = .1 # Gateway Maximum Throughput
......@@ -29,7 +31,7 @@ k1 = .1 * k # Saturated Slope
# QoS
epsilon = .1
epsilon_energy = .2
epsilon_energy = .5
beta = .1
delta = 1.
eta = .1
......@@ -37,7 +39,9 @@ eta = .1
#Energy
energy_model = 1
MinLevel = 0
gamma_per_node = .15
energy_max_per_node = .3
energy_max = NumberNodes * energy_max_per_node
gamma_per_node = .32
gamma = NumberNodes * gamma_per_node
energy_harvesting = 1
energy_tx = 1
......@@ -46,4 +50,15 @@ energy_pl = 1
if not energy_model:
energy_harvesting = 0
energy_tx = 0
energy_pl = 0
\ No newline at end of file
energy_pl = 0
def get_config_string(params):
modname = globals()['__name__']
modobj = sys.modules[modname]
s = "Params: {"
for name in params:
s += str(name) + '=' + str(getattr(modobj, name)) + ','
s = s[0:len(s)-1] + "}"
return s
\ No newline at end of file
......@@ -2,12 +2,13 @@ import numpy as np
from tqdm import tqdm
from time import sleep
import math
import node
import importlib
import config
importlib.reload(config)
import objects
importlib.reload(objects)
#importlib.reload(objects)
import utils
importlib.reload(utils)
......@@ -50,4 +51,43 @@ ntx = StatsObj.Ntx / NumberNodes
npl = StatsObj.Npl / NumberNodes
nidle = 1 - ntx -npl
utils.phase_diagram(npl,ntx,nidle)
\ No newline at end of file
utils.phase_diagram(npl,ntx,nidle,params_s=['NumberNodes','gamma_per_node','alfa','lda','mu','k','v'])
#matrix = np.zeros((2,SimulationTimeslots))
#matrix[0,:] = npl
#matrix[1,:] = ntx
#np.save('N='+str(config.NumberNodes)+'.npy',matrix)
nt_lin_eq = 1/(1+config.k*config.v/(1-config.alfa)*(1/config.lda+config.alfa/config.mu))
np_lin_eq = 1/(1+config.mu/config.alfa*(1/config.lda+(1-config.alfa)/(config.k*config.v)))
eta = (1/config.alfa)*(1-config.alfa)*config.mu/(config.k*config.v)
np_eq_energy = config.energy_max_per_node/(eta + 1)
nt_eq_energy = eta * config.energy_max_per_node / (eta + 1)
if nt_lin_eq + np_lin_eq > config.energy_max_per_node:
nt_eq = nt_eq_energy
np_eq = np_eq_energy
else:
nt_eq = nt_lin_eq
np_eq = np_lin_eq
N = 100
EnergyGridAvg = np.zeros(config.SimulationTimeslots)
for i in range(N,config.SimulationTimeslots):
EnergyGridAvg[i] = np.mean(StorageObj.EnergyGrid[i-N:i])
from matplotlib import pyplot as plt
plt.figure()
plt.plot(StorageObj.EnergyGrid)
plt.axhline(y = (nt_eq + np_eq)*config.NumberNodes*(1-math.sqrt(config.gamma_per_node/(nt_eq + np_eq))), color= 'red', label = 'Model (equilibrium)')
plt.plot(EnergyGridAvg, label = 'Moving Average (N = '+str(N)+')')
plt.ylabel('Energy purchased Power Grid')
plt.xlabel('Time (s)')
plt.legend()
plt.show()
......@@ -147,7 +147,7 @@ class IoTNode:
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):
if StorageObj.Status == 'Charged' and StatsObj.energy_region(self.Timestamp):
self.update_status('Processing')
StorageObj.storage_access(self.Timestamp,-self.EnergyProcessing)
StatsObj.add_processor(self.Timestamp)
......@@ -182,7 +182,7 @@ class IoTNode:
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):
if StorageObj.Status == 'Charged' and StatsObj.energy_region(self.Timestamp):
self.update_status('Transmitting')
self.BitRate[self.Timestamp] = 1 / (config.v * TxTime)
ChannelObj.channel_access(self.NodeId,self.Timestamp)
......
import channel
import storage
import stats
import importlib
import config
importlib.reload(config)
import channel
import storage
import stats
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)
StorageObj = storage.Storage(100, config.gamma, config.MinLevel, Grid = True)
StatsObj = stats.Stats()
HarvestingVector = np.random.poisson(config.gamma, config.SimulationTimeslots)
\ No newline at end of file
......@@ -15,8 +15,8 @@ class Stats:
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):
def energy_region(self,timestamp):
if self.Ntx[timestamp] + self.Npl[timestamp] < config.energy_max:
value = 1
else:
value = 0
......
......@@ -6,7 +6,9 @@ importlib.reload(config)
import numpy as np
class Storage:
def __init__(self, Capacity, HarvestingRate, MinLevel):
def __init__(self, Capacity, HarvestingRate, MinLevel, Grid):
self.Grid = Grid
self.EnergyGrid = np.zeros(config.SimulationTimeslots)
self.Capacity = Capacity
self.HarvestingRate = HarvestingRate
self.HarvestingAmount = 1
......@@ -24,6 +26,9 @@ class Storage:
if (Timestamp < config.SimulationTimeslots):
if self.Level + NetFlow >= self.MinLevel:
self.Level += NetFlow
else:
if self.Grid:
self.EnergyGrid[Timestamp] = self.EnergyGrid[Timestamp] - NetFlow
self.update_state()
def get_level(self):
......
......@@ -44,6 +44,33 @@ def ode(values, t, params):
else:
derivs = [alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) - k*Y*v]
if X+Y>= config.energy_max_per_node:
derivs = [-(1-config.alfa)*config.mu*X + config.alfa*k*v*Y, (1-config.alfa)*config.mu*X - config.alfa*k*v*Y]
else:
derivs = [alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) - k*Y*v]
return derivs
def ode2(values, t, params):
X, Y = values # unpack current values
alfa,lda,mu,c,k,k1,v,restrictions = params # unpack parameters
if X < 0:
X = 0
if Y < 0:
Y = 0
if X > 1:
X = 1
if Y > 1:
Y = 1
if np.random.random() <= config.gamma_per_node / (X+Y):
derivs = [alfa*lda*(1-X-Y)-X*mu,(1-alfa)*lda*(1-X-Y) - k*Y*v]
else:
derivs = [-X*mu,- k*Y*v]
return derivs
def ode3(values, t, params):
X, Y = values # unpack current values
alfa,lda,mu,c,k,k1,v,restrictions = params # unpack parameters
derivs = [alfa*lda*(1-X-Y)*config.gamma_per_node / (X+Y)-X*mu,(1-alfa)*lda*(1-X-Y)*config.gamma_per_node / (X+Y) - k*Y*v]
return derivs
......@@ -56,11 +83,11 @@ def solve_ode(t_start,t_end,y0,alfa,lda,mu,c,k,k1,v,restrictions = False):
return n_p,n_t,t
def phase_diagram(npl,ntx,nidle):
def phase_diagram(npl,ntx,nidle,params_s):
constante = gamma / (NumberNodes*(1 - epsilon_energy))
np_eq_energy = constante / (1 + (mu/(k*v))*(1/alfa-1))
nt_eq_energy = gamma / (NumberNodes*(1 - epsilon_energy)) - np_eq_energy
eta = (1/config.alfa)*(1-config.alfa)*config.mu/(config.k*config.v)
np_eq_energy = config.energy_max_per_node/(eta + 1)
nt_eq_energy = eta * config.energy_max_per_node / (eta + 1)
plt.figure(figsize = (12,10))
plt.plot(time,ntx,color='r',label='Transmitting',alpha=.5)
......@@ -69,7 +96,8 @@ def phase_diagram(npl,ntx,nidle):
plt.plot(time,ntx+npl+nidle,color='k',label='Total',alpha=.5)
plt.axhline(y = np_eq_energy, linestyle='dotted',color='green', label = 'Equilibrio en recta de energia (np)')
plt.axhline(y = nt_eq_energy, linestyle='dotted',color='red',label = 'Equilibrio en recta de energia (nt)')
plt.xlabel('Time (s)')
plt.xlabel('Time (s) \n' + config.get_config_string(params_s))
plt.legend()
plt.show()
......@@ -114,25 +142,26 @@ def phase_diagram(npl,ntx,nidle):
else:
dy[i][j] = (1-alfa)*lda*(1-X[i][j]-Y[i][j])-k*v*Y[i][j]
constante = gamma / (NumberNodes*(1 - epsilon_energy))
np_eq_energy = constante / (1 + (mu/(k*v))*(1/alfa-1))
#constante = gamma / (NumberNodes*(1 - epsilon_energy))
#np_eq_energy = constante / (1 + (mu/(k*v))*(1/alfa-1))
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.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_eq_energy,gamma / (NumberNodes*(1 - epsilon_energy)) - np_eq_energy,marker='*',c='green',s=150,label='Equilibrium Energy Boundary (Model)')
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.plot(np.linspace(0,1,100),config.energy_max_per_node - np.linspace(0,1,100),'-',label = 'Energetic Boundary')
plt.scatter(np_eq_energy,nt_eq_energy,marker='*',c='green',s=150,label='Equilibrium Energy Boundary (Model)')
#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.xlabel(config.get_config_string(params_s))
plt.savefig('Evolucion del sistema.pdf')
plt.show()
......@@ -149,8 +178,8 @@ def phase_diagram(npl,ntx,nidle):
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.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)
......@@ -158,8 +187,11 @@ def phase_diagram(npl,ntx,nidle):
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.xlabel('Time (s) \n' + config.get_config_string(params_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
plt.show()
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