feat(weekly_hydro): weekly hydro introduced

Author: MPavicevic

scripts/run_pownet.py

@@ -1,3 +1,4 @@
+import os
 from datetime import datetime
 from pownet.core import (
     DataProcessor,
@@ -9,29 +10,32 @@
 )
 from pownet.data_utils import create_init_condition
 
+scen_path = ".\\model_library"
 ##### User inputs #####
 to_process_inputs = True
-sim_horizon = 24
-model_name = "dummy"
+sim_horizon = 24*7
+model_name = "cumberland"
 steps_to_run = 3  # Default is None
 do_plot = True
 #######################
 
 if to_process_inputs:
     data_processor = DataProcessor(
-        input_folder="../model_library", model_name=model_name, year=2016, frequency=50
+        input_folder=scen_path, model_name=model_name, year=2016, frequency=50
     )
     data_processor.execute_data_pipeline()
 
 inputs = SystemInput(
-    input_folder="../model_library",
+    input_folder=scen_path,
     model_name=model_name,
     year=2016,
     sim_horizon=sim_horizon,
-    spin_reserve_factor=0.15,
+    spin_reserve_factor=0,
     load_shortfall_penalty_factor=1000,
     load_curtail_penalty_factor=1,
     spin_shortfall_penalty_factor=1000,
+    line_capacity_factor=1,
+    line_loss_factor=0
 )
 inputs.load_and_check_data()
 
@@ -42,14 +46,15 @@
 build_times = []
 opt_times = []
 objvals = []
-init_conditions = create_init_condition(inputs.thermal_units)
+init_conditions = create_init_condition(inputs.thermal_units) # init conditions should be empty dict if we dont have any thermal generators
 
 if steps_to_run is None:
     steps_to_run = 10  # 365 - (sim_horizon // 24 - 1)
 
-for step_k in range(1, steps_to_run):
+start_day = 1
+for step_k in range(start_day, steps_to_run):
     start_time = datetime.now()
-    if step_k == 1:
+    if step_k == start_day:
         power_system_model = model_builder.build(
             step_k=step_k,
             init_conds=init_conditions,
@@ -61,6 +66,7 @@
         )
     build_times.append((datetime.now() - start_time).total_seconds())
 
+    power_system_model.model.write(str(step_k) + '.lp')
     power_system_model.optimize(mipgap=0.001)
 
     # Raise an error if the model is not feasible
@@ -79,45 +85,65 @@
     )
     init_conditions = record.get_init_conds()
 
+import os
+output_folder = ".//outputs"
+if not os.path.exists(output_folder):
+    os.makedirs(output_folder)
+record.write_simulation_results(output_folder)
 
+
+
+# output_folder = ".//outputs"
+# if not os.path.exists(output_folder):
+#     os.makedirs(output_folder)
+# results.write_results(output_folder)
+
+import pandas as pd
 # Process the results
-output_processor = OutputProcessor(inputs)
+output_processor = OutputProcessor()
+output_processor.load(inputs=inputs)
+# output_processor = OutputProcessor(inputs)
 node_var_df = record.get_node_variables()
-output_processor.load_from_dataframe(node_var_df)
+output_processor.get_hourly_generation(node_var_df)
+# output_processor.load_from_dataframe(node_var_df)
 
+output_processor.get_thermal_unit_hourly_dispatch(node_var_df)
+a = output_processor.get_unit_hourly_generation(node_var_df)
 # Visualize the results
 if do_plot:
     visualizer = Visualizer(inputs.model_id)
-    if steps_to_run <= 3:
+    if steps_to_run <= 7:
         visualizer.plot_fuelmix_bar(
-            dispatch=output_processor.get_hourly_dispatch(),
-            demand=output_processor.get_hourly_demand(),
+            dispatch=output_processor.get_hourly_generation(node_var_df),
+            demand=output_processor.get_hourly_demand(inputs.demand),
         )
+        visualizer.plot_lmp(record.lmp_df)
     else:
         visualizer = Visualizer(inputs.model_id)
         visualizer.plot_fuelmix_area(
             dispatch=output_processor.get_daily_dispatch(),
             demand=output_processor.get_daily_demand(),
         )
 
-# # Save other modeling statistics
-# import os
-# from pownet.folder_utils import get_output_dir
-# import pandas as pd
-
-# folder_dir = get_output_dir()
-# prefix_name = "warmstart"
-# df = pd.DataFrame(
-#     {
-#         "build_time": build_times,
-#         "opt_time": opt_times,
-#         "objval": objvals,
-#     }
-# )
-# df.to_csv(
-#     os.path.join(
-#         folder_dir,
-#         f"branch{prefix_name}_{inputs.model_name}_{sim_horizon}_modelstats.csv",
-#     ),
-#     index=False,
-# )
+
+# Save other modeling statistics
+import os
+from pownet.folder_utils import get_output_dir
+import pandas as pd
+
+folder_dir = get_output_dir()
+prefix_name = "warmstart"
+df = pd.DataFrame(
+    {
+        "build_time": build_times,
+        "opt_time": opt_times,
+        "objval": objvals,
+    }
+)
+df.to_csv(
+    os.path.join(
+        folder_dir,
+        f"branch{prefix_name}_{inputs.model_name}_{sim_horizon}_modelstats.csv",
+    ),
+    index=False,
+)

scripts/run_quickstart.py

@@ -9,15 +9,15 @@ def main():
 
     # --------- User inputs
 
-    input_folder = "..//model_library"
-    output_folder = "..//outputs"
+    input_folder = ".//model_library"
+    output_folder = ".//outputs"
 
-    model_name = "dummy"
+    model_name = "cumberland"
     model_year = 2016
 
     # Simulation parameters
     sim_horizon = 24
-    steps_to_run = 2  # 2 Simulation days or 48 hours
+    steps_to_run = 30  # 2 Simulation days or 48 hours
     solver = "gurobi"  # or highs
 
     # --------- End of user inputs
@@ -44,6 +44,8 @@ def main():
     simulator.plot_fuelmix("bar", output_folder)
     simulator.plot_thermal_units(output_folder)
 
+    simulator.plot_thermal_units(output_folder)
+
 
 if __name__ == "__main__":
     main()

src/pownet/core/builder.py

@@ -122,6 +122,7 @@ def __init__(self, inputs: SystemInput) -> None:
 
         self.c_hydro_curtail_ess = gp.tupledict()
         self.c_daily_hydro_curtail_ess = gp.tupledict()
+        self.c_weekly_hydro_curtail_ess = gp.tupledict()
         self.c_solar_curtail_ess = gp.tupledict()
         self.c_wind_curtail_ess = gp.tupledict()
         self.c_import_curtail_ess = gp.tupledict()
@@ -130,11 +131,13 @@ def __init__(self, inputs: SystemInput) -> None:
         # and provides indicators if the unit is online
         self.c_link_hydro_pu = gp.tupledict()
         self.c_link_daily_hydro_pu = gp.tupledict()
+        self.c_link_weekly_hydro_pu = gp.tupledict()
         self.c_link_solar_pu = gp.tupledict()
         self.c_link_wind_pu = gp.tupledict()
         self.c_link_import_pu = gp.tupledict()
 
         self.c_hydro_limit_daily = gp.tupledict()
+        self.c_hydro_limit_weekly = gp.tupledict()
 
         # Energy storage constraints
         self.c_link_ess_charge = gp.tupledict()
@@ -380,6 +383,12 @@ def _build_rnw_import_storage_objfunc_terms(self, step_k: int) -> gp.LinExpr:
                 self.inputs.daily_hydro_must_take_units,
                 self.inputs.nondispatch_contracts,
             ),
+            # Weekly hydro
+            (
+                self.phydro_curtail,
+                self.inputs.weekly_hydro_must_take_units,
+                self.inputs.nondispatch_contracts,
+            ),
             (self.psolar, self.inputs.solar_units, self.inputs.nondispatch_contracts),
             (
                 self.psolar_curtail,
@@ -531,6 +540,7 @@ def _add_daily_hydropower_constraints(self, step_k: int) -> None:
             sim_horizon=self.inputs.sim_horizon,
             hydro_units=self.inputs.daily_hydro_unit_node.keys(),
             hydro_capacity=self.inputs.daily_hydro_capacity,
+            hydro_capacity_min=self.inputs.hydro_capacity_min,
         )
         # (2) The curtailment is enforced for the day and not the hour.
         self.c_daily_hydro_curtail_ess = modeling.add_c_unit_curtail_ess_daily(
@@ -546,10 +556,36 @@ def _add_daily_hydropower_constraints(self, step_k: int) -> None:
             ess_attached=self.inputs.ess_daily_hydro_units,
         )
 
+    def _add_weekly_hydropower_constraints(self, step_k: int) -> None:
+        # (1) Define the weekly upper bound.
+        self.c_hydro_limit_weekly = modeling.add_c_hydro_limit_weekly(
+            model=self.model,
+            phydro=self.phydro,
+            step_k=step_k,
+            sim_horizon=self.inputs.sim_horizon,
+            hydro_units=self.inputs.weekly_hydro_unit_node.keys(),
+            hydro_capacity=self.inputs.weekly_hydro_capacity,
+            hydro_capacity_min=self.inputs.hydro_capacity_min,
+        )
+        # (2) The curtailment is enforced for the day and not the hour.
+        # self.c_weekly_hydro_curtail_ess = modeling.add_c_unit_curtail_ess_weekly(
+        #     model=self.model,
+        #     pdispatch=self.phydro,
+        #     pcurtail=self.phydro_curtail,
+        #     pcharge=self.pcharge,
+        #     type="hydro",
+        #     sim_horizon=self.inputs.sim_horizon,
+        #     step_k=step_k,
+        #     units=self.inputs.weekly_hydro_unit_node.keys(),
+        #     capacity_df=self.inputs.weekly_hydro_capacity,
+        #     ess_attached=self.inputs.ess_weekly_hydro_units,
+        # )
+
     def _add_hydropower_constraints(self, step_k: int) -> None:
 
         self._add_hourly_hydropower_constraints(step_k=step_k)
         self._add_daily_hydropower_constraints(step_k=step_k)
+        self._add_weekly_hydropower_constraints(step_k=step_k)
 
         # Does not update every step_k for this constraint
         self.c_link_daily_hydro_pu = modeling.add_c_link_unit_pu_constant(
@@ -561,6 +597,16 @@ def _add_hydropower_constraints(self, step_k: int) -> None:
             contracted_capacity=self.inputs.hydro_contracted_capacity,
         )
 
+        # Does not update every step_k for this constraint
+        self.c_link_weekly_hydro_pu = modeling.add_c_link_unit_pu_constant(
+            model=self.model,
+            pdispatch=self.phydro,
+            u=self.uhydro,
+            timesteps=self.timesteps,
+            units=self.inputs.weekly_hydro_unit_node.keys(),
+            contracted_capacity=self.inputs.hydro_contracted_capacity,
+        )
+
     def _update_hydropower_constraints(self, step_k: int) -> None:
         # Hourly constraints
         self.model.remove(self.c_link_hydro_pu)
@@ -571,6 +617,10 @@ def _update_hydropower_constraints(self, step_k: int) -> None:
         self.model.remove(self.c_daily_hydro_curtail_ess)
         self._add_daily_hydropower_constraints(step_k=step_k)
 
+        self.model.remove(self.c_hydro_limit_weekly)
+        # self.model.remove(self.c_weekly_hydro_curtail_ess)
+        self._add_weekly_hydropower_constraints(step_k=step_k)
+
     def _add_unit_link_pu(self, step_k: int) -> None:
         # Define parameters
         unit_params = {
@@ -1127,10 +1177,13 @@ def _update_constraints(self, step_k, init_conds: dict) -> None:
         - c_link_import_pu: Import capacity is a timeseries
 
         - c_hydro_limit_daily: Hydropower capacity is a timeseries
+        - c_hydro_limit_weekly: Hydropower capacity is a timeseries
 
         *Curtailment*
         - c_hydro_curtail_ess: Hydropower capacity is a timeseries
         - c_daily_hydro_curtail_ess: Daily hydropower capacity is a timeseries
+        - c_weekly_hydro_curtail_ess: Weekly hydropower capacity is a timeseries
+
         - c_solar_curtail_ess: Solar capacity is a timeseries
         - c_wind_curtail_ess: Wind capacity is a timeseries
         - c_import_curtail_ess: Import capacity is a timeseries
@@ -1376,15 +1429,18 @@ def print_added_constraints(self):
             "c_reserve_req",
             "c_hydro_curtail_ess",
             "c_daily_hydro_curtail_ess",
+            "c_weekly_hydro_curtail_ess",
             "c_solar_curtail_ess",
             "c_wind_curtail_ess",
             "c_import_curtail_ess",
             "c_link_hydro_pu",
             "c_link_daily_hydro_pu",
+            "c_link_weekly_hydro_pu",
             "c_link_solar_pu",
             "c_link_wind_pu",
             "c_link_import_pu",
             "c_hydro_limit_daily",
+            "c_hydro_limit_weekly",
             "c_link_ess_charge",
             "c_link_discharge",
             "c_link_ess_state",