diff --git a/Utilities/LoadProfile.py b/Utilities/LoadProfile.py index cbd5f65..6cfd060 100644 --- a/Utilities/LoadProfile.py +++ b/Utilities/LoadProfile.py @@ -1,16 +1,40 @@ import numpy as np -from Utilities.Time import generate_timestrings, index_peak_times, index_operating_hours +from Utilities.Time import ( + generate_timestrings, + index_peak_times, + index_operating_hours, + check_is_weekday, +) from scipy.optimize import root_scalar import pandas as pd +import matplotlib.pyplot as pl -def get_no_of_peaks(peak_bounds): - peak_occurences = np.random.randint(peak_bounds["min"], peak_bounds["max"], 1) +def get_no_of_peaks(peak_bounds, is_weekday, peak_duration, peak_energy, site): + # calculate theoretically maximum number of peaks based on daily consumption (kWh) + max_occurences = np.floor( + peak_energy / site["maximum_demand_kW"] / (peak_duration["max"] / 60) + ) + + if is_weekday: + peak_occurences = np.random.randint( + peak_bounds["weekdays"]["min"], max_occurences, 1 + ) + else: + peak_occurences = np.random.randint( + peak_bounds["weekends"]["min"], max_occurences, 1 + ) return peak_occurences -def generate_peak_info(c, dt): - no_of_peaks = get_no_of_peaks(c["site_info"]["no_of_peaks"]) +def generate_peak_info(c, dt, is_weekday, peak_energy, site): + no_of_peaks = get_no_of_peaks( + c["site_info"]["no_of_peaks"], + is_weekday, + c["site_info"]["peak_duration"], + peak_energy, + site, + ) operating_hours = generate_timestrings( c["site_info"]["operating hours"]["start"], c["site_info"]["operating hours"]["end"], @@ -26,6 +50,17 @@ def generate_peak_info(c, dt): return peak_times, peak_durations +def generate_peak_profile(idx_peak, c, site): + # Generate a peak profile based on the peak indices and site information + peak_profile = np.zeros(len(idx_peak)) + for i in range(1, np.max(idx_peak) + 1): + peak_profile[idx_peak == i] = site["maximum_demand_kW"] * np.random.uniform( + 1 - c["noise"]["range"], + 1 + c["noise"]["range"], + ) + return peak_profile + + def generate_out_of_hours_consumption_ratio(c): # Generate a random ratio for out-of-hours consumption ratio = np.random.uniform( @@ -35,18 +70,22 @@ def generate_out_of_hours_consumption_ratio(c): return ratio -def recompute_load_profile( - load_profile, - offset, -): +def recompute_load_profile(load_profile, offset, noise, peak_profile): + # apply noise to the load profile, including max demand + load_profile = load_profile * noise # apply offset to the load profile load_profile += offset + # overwrite with peak profile + for i in range(len(peak_profile)): + if peak_profile[i] > 0: + load_profile[i] = peak_profile[i] + return load_profile -def get_load_profile(c, dt, batch_start_time, batch_process_duration): +def get_load_profiles(c, dt, batch_start_time, batch_process_duration): # Generate load profile for each site # c is the configuration dictionary @@ -56,14 +95,26 @@ def get_load_profile(c, dt, batch_start_time, batch_process_duration): # start with indexing all the peak occurences # generate timeseries from start to end time + hours2seconds = 3600 + # check day of the week + is_weekday = check_is_weekday(batch_start_time) start_time = batch_start_time end_time = start_time + batch_process_duration - batch_process_duration_hours = batch_process_duration / 3600 # convert to hours - timestamps = np.arange(start_time, end_time + 1, dt) + batch_process_duration_hours = ( + batch_process_duration / hours2seconds + ) # convert to hours + timestamps = np.arange(start_time, end_time, dt) idx_operating_hours = index_operating_hours( timestamps, c["site_info"]["operating hours"] ) - no_of_operating_hours = np.sum(idx_operating_hours > 0) + no_of_operating_hours = ( + np.sum(idx_operating_hours) + / len(idx_operating_hours) + * batch_process_duration_hours + ) + + # initialise load profiles DataFrame + load_profiles = pd.DataFrame(index=timestamps) # loop through each site in the configuration for site in c["site_info"]["sites"]: @@ -71,14 +122,13 @@ def get_load_profile(c, dt, batch_start_time, batch_process_duration): load_profile = np.zeros(len(timestamps)) # generate noise to make the profile more realistic - noise = np.random.normal( + noise = np.random.uniform( 1 - c["noise"]["range"], 1 + c["noise"]["range"], len(timestamps) ) - # Generate peak times and durations - peak_times, peak_durations = generate_peak_info(c, dt) - # Generate peak times and durations - idx_peak = index_peak_times(timestamps, peak_times, peak_durations) + # make every 2 seconds the same + for i in range(0, len(noise), 2): + noise[i : i + 2] = noise[i] # Generate out-of-hours consumption ratio # The % of energy used outside of the operating hours @@ -91,20 +141,66 @@ def get_load_profile(c, dt, batch_start_time, batch_process_duration): ) out_of_hours_consumption = site["daily_consumption_kWh"] * out_of_hours_ratio - avg_operating_hour_consumption = ( - operating_hour_consumption / no_of_operating_hours - ) - avg_out_of_hours_consumption = out_of_hours_consumption / ( + avg_operating_hour_power = operating_hour_consumption / no_of_operating_hours + avg_out_of_hours_power = out_of_hours_consumption / ( batch_process_duration_hours - no_of_operating_hours ) + # baseline operating hour power is 40% higher than out-of-hours power + gain = 1.4 + assumed_operating_baseline_power = avg_out_of_hours_power * gain + baseline_energy = avg_out_of_hours_power * ( + batch_process_duration_hours - no_of_operating_hours + ) + (assumed_operating_baseline_power * no_of_operating_hours) + + peak_energy = site["daily_consumption_kWh"] - baseline_energy + + # Generate peak times and durations + peak_times, peak_durations = generate_peak_info( + c, dt, is_weekday, peak_energy, site + ) + + # Generate peak times and durations + idx_peak = index_peak_times(timestamps, peak_times, peak_durations) + + # Generate peak profile + peak_profile = generate_peak_profile(idx_peak, c, site) + # assign base load profile - load_profile[idx_operating_hours > 0] = avg_operating_hour_consumption - load_profile[idx_operating_hours == 0] = avg_out_of_hours_consumption + load_profile[idx_operating_hours > 0] = avg_operating_hour_power + load_profile[idx_operating_hours == 0] = avg_out_of_hours_power - # apply peak loads - for i in range(1, np.max(idx_peak) + 1): - load_profile[idx_peak == i] = site["maximum_demand_kW"] + # smoothen out sharp edges + load_profile = np.convolve(load_profile, np.ones(40) / 40, mode="same") - # apply noise to the load profile, including max demand - load_profile = load_profile * noise + def objective(x): + # Objective function to minimize the difference between the load profile and the target profile + # x is the offset + adjusted_profile = recompute_load_profile( + load_profile, x, noise, peak_profile + ) + # get energy consumption in kWh + energy_consumption = np.sum(adjusted_profile) * dt / 3600 + target_consumption = site["daily_consumption_kWh"] + delta = energy_consumption - target_consumption + return delta + + # Use root_scalar to find the optimal offset + result = root_scalar( + objective, + bracket=[-site["maximum_demand_kW"] * 10, site["maximum_demand_kW"] * 10], + method="bisect", + ) + + if result.converged: + offset = result.root + else: + raise ValueError("Root finding did not converge") + + # Recompute the load profile with the optimal offset + load_profile = recompute_load_profile(load_profile, offset, noise, peak_profile) + + # Add the load profile to the DataFrame + load_profiles[site["name"]] = load_profile + + return load_profiles diff --git a/Utilities/Time.py b/Utilities/Time.py index a368cf4..e3412d3 100644 --- a/Utilities/Time.py +++ b/Utilities/Time.py @@ -74,3 +74,14 @@ def index_operating_hours(timestamps, operating_hours): operating_indices[i] = 1 # mark as operating hour return operating_indices + + +def check_is_weekday(batch_start_time): + """Checks if the batch start time is on a weekday.""" + # batch_start_time is in seconds since the epoch + start_time = datetime.fromtimestamp(batch_start_time) + if start_time.weekday() >= 5: # Saturday or Sunday + is_weekday = False + else: + is_weekday = True + return is_weekday diff --git a/YAMLs/config.yml b/YAMLs/config.yml index 283dfc3..339a45c 100644 --- a/YAMLs/config.yml +++ b/YAMLs/config.yml @@ -4,7 +4,7 @@ sim_time: duration_days: 60 noise: - range: 0.15 + range: 0.3 paths: site_info: YAMLs/site_info.yaml diff --git a/YAMLs/site_info.yaml b/YAMLs/site_info.yaml index bf748d5..a8e8bde 100644 --- a/YAMLs/site_info.yaml +++ b/YAMLs/site_info.yaml @@ -26,8 +26,10 @@ operating hours: end: "19:00" time zone: Asia/Kuala_Lumpur no_of_peaks: - min: 30 - max: 100 + weekdays: + min: 5 + weekends: + min: 1 peak_duration: unit: minutes min: 1 diff --git a/main.py b/main.py index e8641cc..6dffd38 100644 --- a/main.py +++ b/main.py @@ -1,6 +1,8 @@ import yaml from Utilities.Time import get_start_time -from Utilities.LoadProfile import get_load_profile +from Utilities.LoadProfile import get_load_profiles +import matplotlib.pyplot as pl +import pandas as pd # read config file c = yaml.safe_load(open("YAMLs/config.yml")) @@ -20,5 +22,17 @@ c["sim_time"]["batch_process_seconds"] = c["sim_time"]["batch_process_hours"] * # load site info c["site_info"] = yaml.safe_load(open(c["paths"]["site_info"])) -# generate load profiles -get_load_profile(c, dt, c["sim_start_time"], c["sim_time"]["batch_process_seconds"]) +cumulative_load_profiles = pd.DataFrame() + +# loop through timesteps +for i in range( + c["sim_start_time"], c["sim_end_time"], c["sim_time"]["batch_process_seconds"] +): + + # generate load profiles + load_profiles = get_load_profiles( + c, dt, c["sim_start_time"], c["sim_time"]["batch_process_seconds"] + ) + + # add to cumulative load profiles + cumulative_load_profiles = pd.concat([cumulative_load_profiles, load_profiles], axis=1