diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..13566b8 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,8 @@ +# Default ignored files +/shelf/ +/workspace.xml +# Editor-based HTTP Client requests +/httpRequests/ +# Datasource local storage ignored files +/dataSources/ +/dataSources.local.xml diff --git a/.idea/ProjectRL2026.iml b/.idea/ProjectRL2026.iml new file mode 100644 index 0000000..f571432 --- /dev/null +++ b/.idea/ProjectRL2026.iml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/inspectionProfiles/profiles_settings.xml b/.idea/inspectionProfiles/profiles_settings.xml new file mode 100644 index 0000000..105ce2d --- /dev/null +++ b/.idea/inspectionProfiles/profiles_settings.xml @@ -0,0 +1,6 @@ + + + + \ No newline at end of file diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..db8786c --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,7 @@ + + + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..ba13012 --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..35eb1dd --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,6 @@ + + + + + + \ No newline at end of file diff --git a/Dam.py b/Dam.py new file mode 100644 index 0000000..913142b --- /dev/null +++ b/Dam.py @@ -0,0 +1,15 @@ + +class dam(object): + def __init__(self, size, ouput_max, input_max): + self.size = size + self.ouput_max = ouput_max + self.input_max = input_max + + self + + def fill(self, flow_in, time = 1) + if flow_in > self.input_max: + flow_in = self.input_max + print("Input flow is to high. Adjusted to max flow") + + diff --git a/Dam2.py b/Dam2.py new file mode 100644 index 0000000..271789d --- /dev/null +++ b/Dam2.py @@ -0,0 +1,84 @@ +import gymnasium as gym +import numpy as np +import pandas as pd +class HydroElectric_Test(gym.Env): +def __init__(self, path_to_test_data:str): +# Define a discrete action space, -1 0 or 1 +self.discrete_action_space = gym.spaces.Discrete(3) +# Define a continuous action space, -1 to 1 +self.continuous_action_space = gym.spaces.Box(low=-1, high=1, shape=(1,), +dtype=np.float32) +# Define the test data +self.test_data = pd.read_excel(path_to_test_data) +self.price_values = self.test_data.iloc[:, 1:25].to_numpy() +self.timestamps = self.test_data['PRICES'] +self.counter = 0 +self.hour = 1 +self.day = 1 +self.state = np.empty(7) +self.max_volume = 100000 # m^3 +self.volume = self.max_volume/2 # m^3 +self.max_flow = 18000 # m^3/h +self.pump_efficiency = 0.8 # - +self.flow_efficiency = 0.9 # - +self.water_mass = 1000 # kg/m^3 +self.dam_height = 30 # m +self.gravity_constant = 9.81 # m/s^2 +self.volume_to_MWh = +(self.water_mass*self.gravity_constant*self.dam_height)*2.77778e-10 # m^3 to MWh +def step(self, action): +reward = 0 +action = np.squeeze(action) # Remove the extra dimension +# Calculate the costs and volume change when pumping water (action >0) +if (action >0) and (self.volume <= self.max_volume): +if (self.volume + action*self.max_flow) > self.max_volume: +action = (self.max_volume - self.volume)/self.max_flow +pumped_water_volume = action * self.max_flow +pumped_water_costs = (1 / 0.8) * pumped_water_volume * +self.volume_to_MWh * self.price_values[self.day-1][self.hour-1] +reward = -pumped_water_costs +self.volume += pumped_water_volume +# Calculate the profits and volume change when selling water (action <0) +elif (action < 0) and (self.volume >= 0): +if (self.volume + action*self.max_flow) < 0: +action = -self.volume/self.max_flow +sold_water_volume = action * self.max_flow +sold_water_profits = +0.9*sold_water_volume*self.volume_to_MWh*self.price_values[self.day-1][self.hour-1] +reward = abs(sold_water_profits) +self.volume -= abs(sold_water_volume) +# No action (action =0) +elif action ==0: +reward = 0 +#volume safeguard +self.volume = np.clip(self.volume, 0, self.max_volume) +self.counter += 1 # Increase the counter +self.hour += 1 # Increase the hour +if self.counter % 24 == 0: # If the counter is a multiple of 24, increase +the day, reset hour to first hour +self.day += 1 +self.hour = 1 +if self.counter == len(self.price_values.flatten())-1: # If the counter is +equal to the number of hours in the test data, terminate the episode +terminated = True +truncated = True +else: # If the counter is not equal to the number of hours in the test +data, continue the episode +terminated = False +truncated = False +info = {} # No info +self.state = self.observation() # Update the state +return self.state, reward, terminated, truncated, info +def observation(self): # Returns the current state +dam_level = self.volume +price = self.price_values[self.day -1][self.hour-1] +hour = self.hour +day_of_week = self.timestamps[self.day -1].dayofweek # Monday = 0, Sunday = +6 +day_of_year = self.timestamps[self.day -1].dayofyear # January 1st = 1, +December 31st = 365 +month = self.timestamps[self.day -1].month # January = 1, December = 12 +year = self.timestamps[self.day -1].year +self.state = np.array([dam_level, price, int(hour), int(day_of_week), +int(day_of_year), int(month), int(year)]) +return self.state \ No newline at end of file diff --git a/TestEnv.py b/TestEnv.py new file mode 100644 index 0000000..900b78a --- /dev/null +++ b/TestEnv.py @@ -0,0 +1,106 @@ +import gymnasium as gym +import numpy as np +import pandas as pd + +class HydroElectric_Test(gym.Env): + + + def __init__(self, path_to_test_data:str): + # Define a discrete action space, -1 0 or 1 + self.discrete_action_space = gym.spaces.Discrete(3) + # Define a continuous action space, -1 to 1 + self.continuous_action_space = gym.spaces.Box(low=-1, high=1, shape=(1,), dtype=np.float32) + # Define the test data + self.test_data = pd.read_excel(path_to_test_data) + self.price_values = self.test_data.iloc[:, 1:25].to_numpy() + self.timestamps = self.test_data['PRICES'] + self.counter = 0 + self.hour = 1 + self.day = 1 + self.state = np.empty(7) + self.max_volume = 100000 # m^3 + self.volume = self.max_volume/2 # m^3 + self.max_flow = 18000 # m^3/h + self.pump_efficiency = 0.8 # - + self.flow_efficiency = 0.9 # - + + self.water_mass = 1000 # kg/m^3 + self.dam_height = 30 # m + self.gravity_constant = 9.81 # m/s^2 + self.volume_to_MWh = (self.water_mass*self.gravity_constant*self.dam_height)*2.77778e-10 # m^3 to MWh + + self.action_space = self.continuous_action_space #PICK ONE + + self.observation_space = gym.spaces.Box( + low=np.array([0, -np.inf, 1, 0, 1, 1, 1900], dtype=np.float32), + high=np.array([self.max_volume, np.inf, 24, 6, 366, 12, 2200], dtype=np.float32), + dtype=np.float32 + ) + + def step(self, action): + reward = 0 + action = np.squeeze(action) # Remove the extra dimension + # Calculate the costs and volume change when pumping water (action >0) + if (action >0) and (self.volume <= self.max_volume): + if (self.volume + action*self.max_flow) > self.max_volume: + action = (self.max_volume - self.volume)/self.max_flow + pumped_water_volume = action * self.max_flow + pumped_water_costs = (1 / 0.8) * pumped_water_volume * self.volume_to_MWh * self.price_values[self.day-1][self.hour-1] + reward = -pumped_water_costs + self.volume += pumped_water_volume + + # Calculate the profits and volume change when selling water (action <0) + elif (action < 0) and (self.volume >= 0): + if (self.volume + action*self.max_flow) < 0: + action = -self.volume/self.max_flow + sold_water_volume = action * self.max_flow + sold_water_profits = 0.9*sold_water_volume*self.volume_to_MWh*self.price_values[self.day-1][self.hour-1] + reward = abs(sold_water_profits) + self.volume -= abs(sold_water_volume) + # No action (action =0) + elif action ==0: + reward = 0 + #volume safeguard + self.volume = np.clip(self.volume, 0, self.max_volume) + + self.counter += 1 # Increase the counter + self.hour += 1 # Increase the hour + if self.counter % 24 == 0: # If the counter is a multiple of 24, increase the day, reset hour to first hour + self.day += 1 + self.hour = 1 + if self.counter == len(self.price_values.flatten())-1: # If the counter is equal to the number of hours in the test data, terminate the episode + terminated = True + truncated = True + else: # If the counter is not equal to the number of hours in the test data, continue the episode + terminated = False + truncated = False + info = {} # No info + self.state = self.observation() # Update the state + + return self.state, reward, terminated, truncated, info + + def observation(self): # Returns the current state + dam_level = self.volume + price = self.price_values[self.day -1][self.hour-1] + hour = self.hour + day_of_week = self.timestamps[self.day -1].dayofweek # Monday = 0, Sunday = 6 + day_of_year = self.timestamps[self.day -1].dayofyear # January 1st = 1, December 31st = 365 + month = self.timestamps[self.day -1].month # January = 1, December = 12 + year = self.timestamps[self.day -1].year + self.state = np.array([dam_level, price, int(hour), int(day_of_week), int(day_of_year), int(month), int(year)]) + + return self.state + + def reset(self, seed=None, options=None): + super().reset(seed=seed) + + self.counter = 0 + self.hour = 1 + self.day = 1 + self.volume = self.max_volume / 2 + + self.state = self.observation() + info = {} + return self.state, info + + diff --git a/figures/01_daily_average_timeseries.png b/figures/01_daily_average_timeseries.png index d1a240c..80e3c96 100644 Binary files a/figures/01_daily_average_timeseries.png and b/figures/01_daily_average_timeseries.png differ diff --git a/figures/02_within_day_std_timeseries.png b/figures/02_within_day_std_timeseries.png index a9c182b..1ced413 100644 Binary files a/figures/02_within_day_std_timeseries.png and b/figures/02_within_day_std_timeseries.png differ diff --git a/figures/03_all_values_hist.png b/figures/03_all_values_hist.png index a9eaa57..0aef798 100644 Binary files a/figures/03_all_values_hist.png and b/figures/03_all_values_hist.png differ diff --git a/figures/04_mean_std_overlay.png b/figures/04_mean_std_overlay.png index 32a2a51..4942ffb 100644 Binary files a/figures/04_mean_std_overlay.png and b/figures/04_mean_std_overlay.png differ diff --git a/figures/05_boxplot_by_hour.png b/figures/05_boxplot_by_hour.png index 79b98c2..f28818f 100644 Binary files a/figures/05_boxplot_by_hour.png and b/figures/05_boxplot_by_hour.png differ diff --git a/figures/06_hourly_mean_curve.png b/figures/06_hourly_mean_curve.png index 07187ee..3bbff02 100644 Binary files a/figures/06_hourly_mean_curve.png and b/figures/06_hourly_mean_curve.png differ diff --git a/figures/07_hourly_mean_with_std_errorbars.png b/figures/07_hourly_mean_with_std_errorbars.png index 4f21d5f..875d0c9 100644 Binary files a/figures/07_hourly_mean_with_std_errorbars.png and b/figures/07_hourly_mean_with_std_errorbars.png differ diff --git a/figures/08_heatmap_date_vs_hour.png b/figures/08_heatmap_date_vs_hour.png index 270ee64..96cc320 100644 Binary files a/figures/08_heatmap_date_vs_hour.png and b/figures/08_heatmap_date_vs_hour.png differ diff --git a/figures/09_hourly_correlation_matrix.png b/figures/09_hourly_correlation_matrix.png index 3d08417..0a5d220 100644 Binary files a/figures/09_hourly_correlation_matrix.png and b/figures/09_hourly_correlation_matrix.png differ diff --git a/figures/10_rolling_mean_daily_avg.png b/figures/10_rolling_mean_daily_avg.png index 0d92926..8e036ab 100644 Binary files a/figures/10_rolling_mean_daily_avg.png and b/figures/10_rolling_mean_daily_avg.png differ diff --git a/figures/11_scatter_daily_avg_vs_within_day_std.png b/figures/11_scatter_daily_avg_vs_within_day_std.png index a762e09..1e4775c 100644 Binary files a/figures/11_scatter_daily_avg_vs_within_day_std.png and b/figures/11_scatter_daily_avg_vs_within_day_std.png differ diff --git a/main (2).py b/main (2).py new file mode 100644 index 0000000..b251bd6 --- /dev/null +++ b/main (2).py @@ -0,0 +1,160 @@ +from TestEnv import HydroElectric_Test +import argparse +import matplotlib.pyplot as plt + +parser = argparse.ArgumentParser() +parser.add_argument('--excel_file', type=str, default='validate.xlsx') # Path to the excel file with the test data +args = parser.parse_args() + +env = HydroElectric_Test(path_to_test_data=args.excel_file) +total_reward = [] +cumulative_reward = [] + +observation = env.observation() + +def heuristic_action1(observation): + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + + if 0 <= hour <= 6: + return 1.0 # pump + elif 17 <= hour <= 21: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_action2(observation): + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + + if 9 <= hour <= 20: + return -1.0 # pump + else: + return 1 + +def heuristic_action3(observation): + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + + if 3 <= hour <= 7: + return 1.0 # pump + elif 11 <= hour <= 14: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_action4(observation): + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + + if 3 <= hour <= 7: + return 1.0 # pump + elif 11 <= hour <= 14: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_mees(observation): + #Takes into account that there is a peak in the winter month later in the day + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 2 <= hour <= 7: + return 1.0 # pump + elif 10 <= hour <= 12: + return -1.0 # sell + elif 18 <= hour <= 21: + return -1.0 # sell + else: + return 0.0 # hold + else: + if 2 <= hour <= 7: + return 1.0 # pump + elif 9 <= hour <= 14: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_mees2(observation, avg): + #Takes into account that there is a peak in the winter month later in the day + #Also uses avrage price to be sure the price is good + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 2 <= hour <= 7 and price < avg: + return 1.0 # pump + elif 10 <= hour <= 12 and price > avg: + return -1.0 # sell + elif 18 <= hour <= 21 and price > avg: + return -1.0 # sell + else: + return 0.0 # hold + else: + if 2 <= hour <= 7 and price < avg: + return 1.0 # pump + elif 9 <= hour <= 14 and price > avg: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_mees3(observation, avg): + #Takes into account that there is a peak in the winter month later in the day + #Also uses avrage price to be sure the price is good + #Include treshholds + + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + treshhold = 0 + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 2 <= hour <= 7 and price < avg-treshhold: + return 1.0 # pump + elif 10 <= hour <= 12 and price > avg+treshhold: + return -1.0 # sell + elif 18 <= hour <= 21 and price > avg+treshhold: + return -1.0 # sell + else: + return 0.0 # hold + else: + if 2 <= hour <= 7 and price < avg-treshhold: + return 1.0 # pump + elif 9 <= hour <= 14 and price > avg+treshhold: + return -1.0 # sell + else: + return 0.0 # hold + +prices = [] + +for i in range(730*24 -1): # Loop through 2 years -> 730 days * 24 hours + # Choose a random action between -1 (full capacity sell) and 1 (full capacity pump) + volume, price, hour, dow, doy, month, year = observation + + + prices.append(price) + + last_prices = prices[-24:] + avg_24h_price = sum(last_prices) / len(last_prices) + + + # action = env.continuous_action_space.sample() + action = heuristic_mees3(observation, avg_24h_price) + + # Or choose an action based on the observation using your RL agent!: + # action = RL_agent.act(observation) + # The observation is the tuple: [volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + next_observation, reward, terminated, truncated, info = env.step(action) + total_reward.append(reward) + cumulative_reward.append(sum(total_reward)) + + done = terminated or truncated + observation = next_observation + + if done: + print('Total reward: ', sum(total_reward)) + # Plot the cumulative reward over time + plt.plot(cumulative_reward) + plt.xlabel('Time (Hours)') + plt.show() + + + + diff --git a/main (3).py b/main (3).py new file mode 100644 index 0000000..f538ef5 --- /dev/null +++ b/main (3).py @@ -0,0 +1,98 @@ +from TestEnv import HydroElectric_Test +import argparse +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt + +parser = argparse.ArgumentParser() +parser.add_argument('--excel_file', type=str, default='validate.xlsx') +args = parser.parse_args() + +env = HydroElectric_Test(path_to_test_data=args.excel_file) + +def heuristic_action3(observation): + volume, price, hour, dow, doy, month, year = observation + if 3 <= hour <= 7: + return np.array([1.0], dtype=np.float32) # pump + elif 11 <= hour <= 14: + return np.array([-1.0], dtype=np.float32) # generate/sell + else: + return np.array([0.0], dtype=np.float32) # hold + +def run_and_plot_average_day(): + obs, _ = env.reset() + + rows = [] + done = False + total_reward = 0.0 + + while not done: + action = heuristic_action3(obs) + next_obs, reward, terminated, truncated, info = env.step(action) + done = terminated or truncated + total_reward += float(reward) + + # Log: hour-of-day is obs[2] (current hour), volume after step is next_obs[0] + rows.append({ + "hour": int(obs[2]), # 1..24 + "price": float(obs[1]), + "action": float(np.squeeze(action)), # -1, 0, 1 + "volume_after": float(next_obs[0]), + "reward": float(reward), + }) + + obs = next_obs + + print("Total reward:", total_reward) + + df = pd.DataFrame(rows) + + # Build "average day" (mean per hour-of-day) + avg = df.groupby("hour").agg( + avg_price=("price", "mean"), + avg_volume=("volume_after", "mean"), + pump_rate=("action", lambda x: np.mean(x > 0)), + gen_rate=("action", lambda x: np.mean(x < 0)), + hold_rate=("action", lambda x: np.mean(x == 0)), + avg_reward=("reward", "mean"), + ).reset_index().sort_values("hour") + + # ---------- Plot 1: Avg volume + avg price ---------- + fig, ax1 = plt.subplots() + + ax1.plot(avg["hour"], avg["avg_volume"]) + ax1.set_xlabel("Hour of day") + ax1.set_ylabel("Average reservoir volume (m³)") + ax1.set_xticks(range(1, 25)) + + ax2 = ax1.twinx() + ax2.plot(avg["hour"], avg["avg_price"]) + ax2.set_ylabel("Average price") + + plt.title("Average day: reservoir volume and price") + plt.show() + + # ---------- Plot 2: Action frequency by hour ---------- + plt.figure() + plt.plot(avg["hour"], avg["pump_rate"], label="Pump frequency") + plt.plot(avg["hour"], avg["gen_rate"], label="Generate frequency") + plt.plot(avg["hour"], avg["hold_rate"], label="Hold frequency") + plt.xlabel("Hour of day") + plt.ylabel("Fraction of days") + plt.title("Average day: how often the heuristic pumps/generates") + plt.xticks(range(1, 25)) + plt.ylim(-0.05, 1.05) + plt.legend() + plt.show() + + # ---------- Plot 3 (optional): Avg reward per hour ---------- + plt.figure() + plt.plot(avg["hour"], avg["avg_reward"]) + plt.xlabel("Hour of day") + plt.ylabel("Average reward per hour") + plt.title("Average day: reward contribution by hour") + plt.xticks(range(1, 25)) + plt.show() + +if __name__ == "__main__": + run_and_plot_average_day() diff --git a/main.py b/main.py new file mode 100644 index 0000000..e0b541a --- /dev/null +++ b/main.py @@ -0,0 +1,51 @@ +import numpy as np + +# import your env class from the file where you defined it +# Example: if your environment code is in hydro_env.py +from hydro_env import HydroElectric_Test + + +def time_window_heuristic(obs): + """ + obs = [dam_level, price, hour, day_of_week, day_of_year, month, year] + Pump 08-10, Generate 15-17, else idle. + """ + hour = int(obs[2]) # 1..24 + + # Pump between 8 and 10 inclusive + if 8 <= hour <= 10: + return np.array([1.0], dtype=np.float32) + + # Generate between 15 and 17 inclusive + if 15 <= hour <= 17: + return np.array([-1.0], dtype=np.float32) + + return np.array([0.0], dtype=np.float32) + + +def run_sim(path_to_xlsx): + env = HydroElectric_Test(path_to_xlsx) + + obs, info = env.reset() + total_reward = 0.0 + + done = False + step_i = 0 + + while not done: + action = time_window_heuristic(obs) + obs, reward, terminated, truncated, info = env.step(action) + + total_reward += float(reward) + done = terminated or truncated + step_i += 1 + + print("Finished simulation.") + print("Steps:", step_i) + print("Total reward (profit):", total_reward) + print("Final reservoir volume (m^3):", env.volume) + + +if __name__ == "__main__": + # Change this to your file path: + run_sim("train.xlsx") diff --git a/make_plots.py b/make_plots.py new file mode 100644 index 0000000..0f210cb --- /dev/null +++ b/make_plots.py @@ -0,0 +1,68 @@ +import pandas as pd +import matplotlib.pyplot as plt + +# =============================== +# LOAD DATA +# =============================== +FILE_PATH = "train.xlsx" # change if needed +DATE_COL = "PRICES" # your date column + +df = pd.read_excel(FILE_PATH) + +# Convert date column to datetime +df[DATE_COL] = pd.to_datetime(df[DATE_COL]) + +# Hour columns (Hour 01 ... Hour 24) +hour_cols = [c for c in df.columns if isinstance(c, str) and c.startswith("Hour ")] +if len(hour_cols) != 24: + raise ValueError(f"Expected 24 hour columns, found {len(hour_cols)}: {hour_cols}") + +# Ensure numeric +df[hour_cols] = df[hour_cols].apply(pd.to_numeric, errors="coerce") + +# Calendar features +df["month"] = df[DATE_COL].dt.month +df["weekday"] = df[DATE_COL].dt.weekday # Mon=0 ... Sun=6 +df["day_type"] = df["weekday"].apply(lambda x: "Weekend" if x >= 5 else "Weekday") + +hours = list(range(1, 25)) # x-axis + +# =============================== +# FIGURE 1: +# Typical 24h profile per month (12 lines) +# =============================== +plt.figure(figsize=(12, 6)) + +for month in range(1, 13): + month_profile = df[df["month"] == month][hour_cols].mean(axis=0) # mean over days, per hour + plt.plot(hours, month_profile.values, + label=pd.to_datetime(str(month), format="%m").strftime("%B")) + +plt.xlabel("Hour of Day") +plt.ylabel("Average Price") +plt.title("Typical Day (24h Profile) by Month") +plt.xticks(hours) +plt.grid(True) +plt.legend(ncol=3) +plt.tight_layout() +plt.show() + +# =============================== +# FIGURE 2: +# Typical weekday vs weekend 24h profile (2 lines) +# =============================== +weekday_profile = df[df["day_type"] == "Weekday"][hour_cols].mean(axis=0) +weekend_profile = df[df["day_type"] == "Weekend"][hour_cols].mean(axis=0) + +plt.figure(figsize=(12, 6)) +plt.plot(hours, weekday_profile.values, label="Weekday (avg)") +plt.plot(hours, weekend_profile.values, label="Weekend (avg)") + +plt.xlabel("Hour of Day") +plt.ylabel("Average Price") +plt.title("Typical Day (24h Profile): Weekday vs Weekend") +plt.xticks(hours) +plt.grid(True) +plt.legend() +plt.tight_layout() +plt.show() diff --git a/plot3.py b/plot3.py new file mode 100644 index 0000000..e7d4a34 --- /dev/null +++ b/plot3.py @@ -0,0 +1,124 @@ +import pandas as pd +import matplotlib.pyplot as plt + +# =============================== +# LOAD DATA +# =============================== +FILE_PATH = "train.xlsx" # change if needed +DATE_COL = "PRICES" # your date column + +df = pd.read_excel(FILE_PATH) + +# Convert date column to datetime +df[DATE_COL] = pd.to_datetime(df[DATE_COL]) + +# Hour columns (Hour 01 ... Hour 24) +hour_cols = [c for c in df.columns if isinstance(c, str) and c.startswith("Hour ")] +if len(hour_cols) != 24: + raise ValueError(f"Expected 24 hour columns, found {len(hour_cols)}: {hour_cols}") + +# Ensure numeric +df[hour_cols] = df[hour_cols].apply(pd.to_numeric, errors="coerce") + +# Calendar features +df["month"] = df[DATE_COL].dt.month +df["weekday"] = df[DATE_COL].dt.weekday # Mon=0 ... Sun=6 +df["day_type"] = df["weekday"].apply(lambda x: "Weekend" if x >= 5 else "Weekday") + +hours = list(range(1, 25)) # x-axis + +# =============================== +# Helper: add buy/sell windows +# =============================== +WINTER_MONTHS = {10, 11, 12, 1, 2} + +def add_trade_windows(ax): + """ + Adds shaded regions for buy/sell windows. + - Buy: 2..7 (both seasons) + - Sell winter: 10..12 and 18..21 + - Sell summer: 9..14 + """ + # Common buy window + buy_window = [(2, 7)] + # Winter vs summer sell windows + winter_sell = [(10, 12), (18, 21)] + summer_sell = [(9, 14)] + + # Shaded bands + for (a, b) in buy_window: + ax.axvspan(a, b, alpha=0.12, label="Buy window" if a == 2 else None) + for (a, b) in summer_sell: + ax.axvspan(a, b, alpha=0.12, label="Sell window (summer)") + for (a, b) in winter_sell: + ax.axvspan(a, b, alpha=0.12, label="Sell window (winter)" if a == 10 else None) + + # Make the band colors explicit via edge colors using lines (keeps default line colors untouched) + # Draw bold boundary lines for readability + def bold_bounds(windows, linestyle="-"): + for (a, b) in windows: + ax.axvline(a, linewidth=3, linestyle=linestyle) + ax.axvline(b, linewidth=3, linestyle=linestyle) + + # Bold boundaries (optional but helps a lot) + bold_bounds(buy_window, linestyle="--") + bold_bounds(summer_sell, linestyle=":") + bold_bounds(winter_sell, linestyle="-.") + +# =============================== +# FIGURE 1: +# Typical 24h profile per month (12 lines) +# with CLEAR buy/sell windows +# =============================== + +fig, ax = plt.subplots(figsize=(13, 6)) + +# ---------- BUY WINDOW (both seasons) ---------- +ax.axvspan(2, 7, + color="green", alpha=0.18, + label="Buy window (all months)") +ax.axvline(2, color="green", linewidth=3, linestyle="--") +ax.axvline(7, color="green", linewidth=3, linestyle="--") + +# ---------- SELL WINDOW (SUMMER) ---------- +ax.axvspan(9, 14, + color="tab:blue", alpha=0.18, + label="Sell window (summer)") +ax.axvline(9, color="tab:blue", linewidth=3, linestyle=":") +ax.axvline(14, color="tab:blue", linewidth=3, linestyle=":") + +# ---------- SELL WINDOWS (WINTER) ---------- +ax.axvspan(10, 12, + color="red", alpha=0.18, + label="Sell window (winter – morning)") +ax.axvspan(18, 21, + color="red", alpha=0.18, + label="Sell window (winter – evening)") +ax.axvline(10, color="red", linewidth=3) +ax.axvline(12, color="red", linewidth=3) +ax.axvline(18, color="red", linewidth=3) +ax.axvline(21, color="red", linewidth=3) + +# ---------- PLOT MONTH PROFILES ---------- +for month in range(1, 13): + month_profile = df[df["month"] == month][hour_cols].mean(axis=0) + ax.plot( + hours, + month_profile.values, + label=pd.to_datetime(str(month), format="%m").strftime("%B"), + linewidth=1.8 + ) + +ax.set_xlabel("Hour of Day") +ax.set_ylabel("Average Price") +ax.set_title("Typical Day (24h Profile) by Month with Buy/Sell Windows") +ax.set_xticks(hours) +ax.grid(True, alpha=0.3) + +# ---------- CLEAN LEGEND (no duplicates) ---------- +handles, labels = ax.get_legend_handles_labels() +unique = dict(zip(labels, handles)) +ax.legend(unique.values(), unique.keys(), ncol=3, fontsize=9) + +plt.tight_layout() +plt.show() diff --git a/plot_behavior.py b/plot_behavior.py new file mode 100644 index 0000000..6a236fb --- /dev/null +++ b/plot_behavior.py @@ -0,0 +1,322 @@ +from TestEnv import HydroElectric_Test +import argparse +import matplotlib.pyplot as plt + +# ----------------------------- +# Heuristics +# ----------------------------- +def heuristic_action1(observation): + # obs =[volume, price, hour_of_day, day_of_week, day_of_year, month_of_year, year] + volume, price, hour, dow, doy, month, year = observation + if 0 <= hour <= 6: + return 1.0 # pump + elif 17 <= hour <= 21: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_action2(observation): + volume, price, hour, dow, doy, month, year = observation + if 9 <= hour <= 20: + return -1.0 + else: + return 1.0 + +def heuristic_action3(observation): + volume, price, hour, dow, doy, month, year = observation + if 3 <= hour <= 7: + return 1.0 # pump + elif 11 <= hour <= 14: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_action4(observation): + volume, price, hour, dow, doy, month, year = observation + if 3 <= hour <= 7: + return 1.0 # pump + elif 11 <= hour <= 14: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_mees(observation): + # Takes into account that there is a peak in the winter month later in the day + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 2 <= hour <= 7: + return 1.0 # pump + elif 10 <= hour <= 12: + return -1.0 # sell + elif 18 <= hour <= 21: + return -1.0 # sell + else: + return 0.0 # hold + else: + if 2 <= hour <= 7: + return 1.0 # pump + elif 9 <= hour <= 14: + return -1.0 # sell + else: + return 0.0 # hold + +def heuristic_mees2(observation, avg): + # Winter peak later in day + compare against 24h average price + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 2 <= hour <= 7 and price < avg: + return 1.0 + elif 10 <= hour <= 12 and price > avg: + return -1.0 + elif 18 <= hour <= 21 and price > avg: + return -1.0 + else: + return 0.0 + else: + if 2 <= hour <= 7 and price < avg: + return 1.0 + elif 9 <= hour <= 14 and price > avg: + return -1.0 + else: + return 0.0 + +def heuristic_mees3(observation, avg, threshold=0.0): + # Winter peak later in day + 24h avg + threshold + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 2 <= hour <= 7 and price < avg - threshold: + return 1.0 + elif 10 <= hour <= 12 and price > avg + threshold: + return -1.0 + elif 18 <= hour <= 21 and price > avg + threshold: + return -1.0 + else: + return 0.0 + else: + if 2 <= hour <= 7 and price < avg - threshold: + return 1.0 + elif 9 <= hour <= 14 and price > avg + threshold: + return -1.0 + else: + return 0.0 + +def heuristic_mees4(observation, avg, threshold=0.0): + # Winter peak later in day + 24h avg + threshold + volume, price, hour, dow, doy, month, year = observation + if month in [11, 10, 12, 1, 2]: + if 1 <= hour <= 6 and price < avg - threshold: + return 1.0 + elif 10 <= hour <= 12 and price > avg + threshold: + return -1.0 + elif 18 <= hour <= 21 and price > avg + threshold: + return -1.0 + else: + return 0.0 + else: + if 1 <= hour <= 6 and price < avg - threshold: + return 1.0 + elif 9 <= hour <= 14 and price > avg + threshold: + return -1.0 + else: + return 0.0 + + +# ----------------------------- +# Plotting helpers +# ----------------------------- +def plot_day_from_logs(logs, target_year, target_doy, title): + """ + Plot one day (up to 24 hours) of: + - price (line) + - action (step) + - volume (dashed) + selected by (year, day_of_year). + """ + years = logs["year"] + doys = logs["doy"] + hours = logs["hour"] + prices = logs["price"] + actions = logs["action"] + volumes = logs["volume"] + + # find indices matching the requested day + idx = [i for i in range(len(actions)) if years[i] == target_year and doys[i] == target_doy] + if not idx: + print(f"[plot_day] No data for year={target_year}, doy={target_doy}") + return + + # Prefer starting at hour==0 (start of day) if present + start_candidates = [i for i in idx if hours[i] == 0] + start = start_candidates[0] if start_candidates else idx[0] + day_idx = list(range(start, min(start + 24, len(actions)))) + + h = [hours[i] for i in day_idx] + p = [prices[i] for i in day_idx] + a = [actions[i] for i in day_idx] + v = [volumes[i] for i in day_idx] + + fig, ax1 = plt.subplots() + ax1.plot(h, p, label = "Price", color = "orange") + ax1.set_xlabel("Hour of day") + ax1.set_ylabel("Price") + + ax2 = ax1.twinx() + ax2.plot(h, v, linestyle="--", label = "Reservoir volume") + ax2.set_ylabel("Action (step) / Volume (dashed)") + + # Combine legends from both axes + lines_1, labels_1 = ax1.get_legend_handles_labels() + lines_2, labels_2 = ax2.get_legend_handles_labels() + ax1.legend(lines_1 + lines_2, labels_1 + labels_2, loc="best") + + plt.title(title) + plt.xticks(range(0, 24, 2)) + plt.show() + + +def auto_pick_day_by_month(logs, target_month): + """ + Pick the first day-of-year (doy) we see for a given month. + Returns (year, doy) or None if not found. + """ + for y, m, d in zip(logs["year"], logs["month"], logs["doy"]): + if m == target_month: + return (y, d) + return None + +def compute_daily_profit(logs, target_year, target_doy): + """ + Returns total profit (sum of rewards) for the given (year, doy). + """ + years = logs["year"] + doys = logs["doy"] + hours = logs["hour"] + rewards = logs["reward"] + + idx = [i for i in range(len(hours)) if years[i] == target_year and doys[i] == target_doy] + if not idx: + return None + + # Prefer a full day starting at hour 0 + start_candidates = [i for i in idx if hours[i] == 0] + start = start_candidates[0] if start_candidates else idx[0] + + day_idx = list(range(start, min(start + 24, len(rewards)))) + + return sum(rewards[i] for i in day_idx) + + + +# ----------------------------- +# Main +# ----------------------------- +def main(): + parser = argparse.ArgumentParser() + parser.add_argument('--excel_file', type=str, default='validate.xlsx', help="Path to excel file with test data") + parser.add_argument('--threshold', type=float, default=0.0, help="Price threshold for heuristic_mees3") + parser.add_argument('--winter_doy', type=int, default=40, help="Explicit winter day-of-year to plot (optional)") + parser.add_argument('--summer_doy', type=int, default=210, help="Explicit summer day-of-year to plot (optional)") + parser.add_argument('--winter_month', type=int, default=1, help="Winter month to auto-pick (default Jan=1)") + parser.add_argument('--summer_month', type=int, default=7, help="Summer month to auto-pick (default Jul=7)") + args = parser.parse_args() + + env = HydroElectric_Test(path_to_test_data=args.excel_file) + + total_reward = [] + cumulative_reward = [] + + observation = env.observation() + + # Logs for plotting behavior + logs = { + "volume": [], + "price": [], + "hour": [], + "dow": [], + "doy": [], + "month": [], + "year": [], + "action": [], + "reward": [], + "cum_reward": [] + } + + prices_window = [] # for 24h rolling average + + # Run full episode (2 years -> 730 days * 24 hours - 1) + for t in range(730 * 24 - 1): + volume, price, hour, dow, doy, month, year = observation + + # update rolling window + prices_window.append(price) + last_prices = prices_window[-24:] + avg_24h_price = sum(last_prices) / len(last_prices) + + # choose action (your algorithm) + action = heuristic_mees3(observation, avg_24h_price, threshold=args.threshold) + + next_observation, reward, terminated, truncated, info = env.step(action) + + total_reward.append(reward) + cumulative_reward.append(sum(total_reward)) + + # log everything + logs["volume"].append(volume) + logs["price"].append(price) + logs["hour"].append(hour) + logs["dow"].append(dow) + logs["doy"].append(doy) + logs["month"].append(month) + logs["year"].append(year) + logs["action"].append(action) + logs["reward"].append(reward) + logs["cum_reward"].append(cumulative_reward[-1]) + + done = terminated or truncated + observation = next_observation + + if done: + break + + print("Total reward:", sum(total_reward)) + + # ----------------------------- + # Choose winter & summer days to plot + # ----------------------------- + # If user provided doy, use those. Otherwise auto-pick a day from winter_month / summer_month. + if args.winter_doy is not None: + winter = (logs["year"][0], args.winter_doy) + else: + winter = auto_pick_day_by_month(logs, args.winter_month) + + if args.summer_doy is not None: + summer = (logs["year"][0], args.summer_doy) + else: + summer = auto_pick_day_by_month(logs, args.summer_month) + + if winter is None: + print(f"Could not find any data for winter_month={args.winter_month}") + else: + wy, wd = winter + plot_day_from_logs(logs, wy, wd, title=f"Behavior on winter day (year={int(wy)}, doy={wd})") + winter_profit = compute_daily_profit(logs, wy, wd) + print("winter profit", winter_profit) + if summer is None: + print(f"Could not find any data for summer_month={args.summer_month}") + else: + sy, sd = summer + plot_day_from_logs(logs, sy, sd, title=f"Behavior on summer day (year={int(sy)}, doy={sd})") + summer_profit = compute_daily_profit(logs, sy, sd) + print("summer profit: ", summer_profit) + + # Optional: plot cumulative reward for the whole run + plt.figure() + plt.plot(logs["cum_reward"]) + plt.xlabel("Time (Hours)") + plt.ylabel("Cumulative reward") + plt.title("Cumulative reward over time") + plt.show() + + +if __name__ == "__main__": + main() diff --git a/validate.xlsx b/validate.xlsx new file mode 100644 index 0000000..0baa422 Binary files /dev/null and b/validate.xlsx differ