paper code

RAMP_v02-pre/stochastic_process.py → RAMP v.0.1-pre.py

@@ -1,11 +1,146 @@
 # -*- coding: utf-8 -*-
+"""
+Created on Fri Jun 08 11:46:00 2018
+This is the code for the open-source stochastic model for the generation of 
+multi-energy load profiles in off-grid areas, called RAMP, v.0.1-pre.
 
+@authors:
+- Francesco Lombardi, Politecnico di Milano
+- Sergio Balderrama, Université de Liège
+- Sylvain Quoilin, KU Leuven
+- Emanuela Colombo, Politecnico di Milano
+
+Copyright 2019 RAMP, contributors listed above.
+Licensed under the European Union Public Licence (EUPL), Version 1.1;
+you may not use this file except in compliance with the License.
+
+Unless required by applicable law or agreed to in writing,
+software distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and limitations
+under the License.
+
+"""
 #%% Import required libraries
+import matplotlib.pyplot as plt
 import numpy as np
 import numpy.ma as ma
 import random 
 import math
-from initialise import User_list, peak_enlarg, mu_peak, s_peak, Profile, num_profiles
+#%% Definition of Python classes that constitute the model architecture
+'''
+The code is based on two concatenated python classes, namely 'User' and
+'Appliance', which are used to define at the outer level the User classes and 
+at the inner level all the available appliances within each user class, with 
+their own characteristics. Within the Appliance class, some other functions are
+created to define windows of use and, if needed, specific duty cycles
+'''
+
+#Define the outer python class that represents 'User classes'
+class User():
+
+    def __init__(self, name = "", n_users = 1, us_pref = 0):
+        self.user_name = name
+        self.num_users = n_users #specifies the number of users within the class
+        self.user_preference = us_pref #allows to check if random number coincides with user preference, to distinguish between various appliance_use options (e.g. different cooking options)
+        self.App_list = [] #each instance of User (i.e. each user class) has its own list of Appliances
+
+#Define the inner class for modelling user's appliances within the correspoding user class
+    class Appliance():
+
+        def __init__(self,user, n = 1, P = 0, w = 1, t = 0, r_t = 0, c = 0, fixed = 'no', fixed_cycle = 0, occasional_use = 1, flat = 'no', thermal_P_var = 0, pref_index = 0):
+            self.user = user #user to which the appliance is bounded
+            self.number = n #number of appliances of the specified kind
+            self.POWER = P #nominal power of appliances of the specified kind
+            self.num_windows = w #number of functioning windows to be considered
+            self.func_time = t #total time the appliance is on during the day
+            self.r_t = r_t #percentage of total time of use that is subject to random variability
+            self.func_cycle = c #minimum time the appliance is kept on after switch-on event 
+            self.fixed = fixed #if 'yes', all the 'n' appliances of this kind are always switched-on together
+            self.activate = fixed_cycle #if equal to 1,2 or 3, respectively 1,2 or 3 duty cycles can be modelled, for different periods of the day
+            self.occasional_use = occasional_use #probability that the appliance is always (i.e. everyday) included in the mix of appliances that the user actually switches-on during the day
+            self.flat = flat #allows to model appliances that are not subject to any kind of random variability, such as public lighting
+            self.Thermal_P_var = thermal_P_var #allows to randomly variate the App power within a range
+            self.Pref_index = pref_index
+
+        def windows(self, w1 = np.array([0,0]), w2 = np.array([0,0]),r_w = 0, w3 = np.array([0,0])):    
+            self.window_1 = w1 #array of start and ending time for window of use #1
+            self.window_2 = w2 #array of start and ending time for window of use #2
+            self.window_3 = w3 #array of start and ending time for window of use #3
+            self.random_var_w = r_w #percentage of variability in the start and ending times of the windows
+            self.daily_use = np.zeros(1440) #create an empty daily use profile
+            self.daily_use[w1[0]:(w1[1])] = np.full(np.diff(w1),0.001) #fills the daily use profile with infinitesimal values that are just used to identify the functioning windows
+            self.daily_use[w2[0]:(w2[1])] = np.full(np.diff(w2),0.001) #same as above for window2
+            self.daily_use[w3[0]:(w3[1])] = np.full(np.diff(w3),0.001) #same as above for window3
+            self.daily_use_masked = np.zeros_like(ma.masked_not_equal(self.daily_use,0.001)) #apply a python mask to the daily_use array to make only functioning windows 'visibile'
+            self.random_var_1 = int(r_w*np.diff(w1)) #calculate the random variability of window1, i.e. the maximum range of time they can be enlarged or shortened
+            self.random_var_2 = int(r_w*np.diff(w2)) #same as above
+            self.random_var_3 = int(r_w*np.diff(w3)) #same as above
+            self.user.App_list.append(self) #automatically appends the appliance to the user's appliance list
+
+            #if needed, specific duty cycles can be defined for each Appliance, for a maximum of three different ones
+        def specific_cycle_1(self, P_11 = 0, t_11 = 0, P_12 = 0, t_12 = 0, r_c1 = 0):
+            self.P_11 = P_11 #power absorbed during first part of the duty cycle
+            self.t_11 = t_11 #duration of first part of the duty cycle
+            self.P_12 = P_12 #power absorbed during second part of the duty cycle
+            self.t_12 = t_12 #duration of second part of the duty cycle
+            self.r_c1 = r_c1 #random variability of duty cycle segments duration
+            self.fixed_cycle1 = np.concatenate(((np.ones(t_11)*P_11),(np.ones(t_12)*P_12))) #create numpy array representing the duty cycle
+
+        def specific_cycle_2(self, P_21 = 0, t_21 = 0, P_22 = 0, t_22 = 0, r_c2 = 0):
+            self.P_21 = P_21 #same as for cycle1
+            self.t_21 = t_21
+            self.P_22 = P_22
+            self.t_22 = t_22
+            self.r_c2 = r_c2
+            self.fixed_cycle2 = np.concatenate(((np.ones(t_21)*P_21),(np.ones(t_22)*P_22)))
+
+        def specific_cycle_3(self, P_31 = 0, t_31 = 0, P_32 = 0, t_32 = 0, r_c3 = 0):
+            self.P_31 = P_31 #same as for cycle1
+            self.t_31 = t_31
+            self.P_32 = P_32
+            self.t_32 = t_32
+            self.r_c3 = r_c3
+            self.fixed_cycle3 = np.concatenate(((np.ones(t_31)*P_31),(np.ones(t_32)*P_32)))
+
+        #different time windows can be associated with different specific duty cycles
+        def cycle_behaviour(self, cw11 = np.array([0,0]), cw12 = np.array([0,0]), cw21 = np.array([0,0]), cw22 = np.array([0,0]), cw31 = np.array([0,0]), cw32 = np.array([0,0])):
+            self.cw11 = cw11 #first window associated with cycle1
+            self.cw12 = cw12 #second window associated with cycle1
+            self.cw21 = cw21 #same for cycle2
+            self.cw22 = cw22
+            self.cw31 = cw31 #same for cycle 3
+            self.cw32 = cw32
+
+#%% Initialisation of a model instance
+'''
+The model is ready to be initialised
+'''
+User_list = [] #creates an empty list to store all the needed User classes
+num_profiles = int(input("please indicate the number of profiles to be generated: ")) #asks the user how many profiles (i.e. code runs) he wants
+print('Please wait...') 
+Profile = [] #creates an empty list to store the results of each code run, i.e. each stochastically generated profile
+
+#%% Definition of the inputs
+'''
+Input data definition (this is planned to be externalised in a separate script)
+'''
+
+#User classes definition
+HI = User("high income",1)
+User_list.append(HI)
+
+HI_Phone_charger = HI.Appliance(HI,1,10,2,300,0.2,5, thermal_P_var=0.2)
+HI_Phone_charger.windows([1110,1440],[0,30],0.35)
+
+#%%
+'''
+Calibration parameters. These can be changed in case the user has some real data against which the model can be calibrated
+They regulate the probabilities defining the largeness of the peak window and the probability of coincident switch-on within the peak window
+'''
+peak_enlarg = 0 #percentage random enlargement or reduction of peak time range length
+mu_peak = 0.5 #median value of gaussian distribution [0,1] by which the number of coincident switch_ons is randomly selected
+s_peak = 1 #standard deviation (as percentage of the median value) of the gaussian distribution [0,1] above mentioned
 
 #%% Core model stochastic script
 '''
@@ -276,3 +411,44 @@
     Profile.append(Tot_Classes) #appends the total load to the list that will contain all the generated profiles
     print('Profile',prof_i+1,'/',num_profiles,'completed') #screen update about progress of computation
 
+#%% Post-processing
+'''
+Just some additional code lines to calculate useful indicators and generate plots
+'''
+
+Profile_avg = np.zeros(1440)
+for pr in Profile:
+    Profile_avg = Profile_avg + pr
+Profile_avg = Profile_avg/len(Profile)
+Profile_avg_kW = Profile_avg/1000
+
+Profile_kW = []
+for kW in Profile:
+    Profile_kW.append(kW/1000)
+
+Profile_series = np.array([])
+for iii in Profile:
+    Profile_series = np.append(Profile_series,iii)
+
+x = np.arange(0,1440,5)
+plt.figure(figsize=(10,5))
+for n in Profile:
+    plt.plot(np.arange(1440),n,'#b0c4de')
+    plt.xlabel('Time (hours)')
+    plt.ylabel('Power (W)')
+    plt.ylim(ymin=0)
+    #plt.ylim(ymax=5000)
+    plt.margins(x=0)
+    plt.margins(y=0)
+plt.plot(np.arange(1440),Profile_avg,'#4169e1')
+plt.xticks([0,240,480,(60*12),(60*16),(60*20),(60*24)],[0,4,8,12,16,20,24])
+#plt.savefig('profile.eps', format='eps', dpi=1000)
+plt.show()
+
+
+#%% Export individual profiles
+'''
+for i in range (len(Profile)):
+    np.save('p0%d.npy' % (i), Profile[i])
+'''
+