Start subset modification

docs/tables/baseline/second_stage/variables_docs.csv

@@ -4,11 +4,11 @@
 ":math:`V_\text{BAS, END}^{b}`",,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
 ":math:`dV_\text{+}^{t,~b}`",,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
 ":math:`dV_\text{-}^{t,~b}`",,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}`
-":math:`P^{t,~h}`",Electrical power produced by a hydropower plant,:math:`\mathbb{R}^{+}`,:math:`\mathrm{MW}`
+":math:`P^{t,~h}`",Electrical power produced (or consumed) by a hydropower plant,:math:`\mathbb{R}^{+}`,:math:`\mathrm{MW}`
 ":math:`Q^{t,~h}`",Amount of water that flows through a hydropower plant,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
 ":math:`S_\text{BAS}^{t,~b,~s_b}`","Binary variable that indicates the basin's volume state derived from its volume:  :math:`S_\text{BAS}^{t, b, s_b} = \begin{cases} 1 \text{ if } V_\text{BAS}^{t, b} \in  \left[ V_\text{BAS, MIN}^{b,~s_b} ;V_\text{BAS, MAX}^{b,~s_b} \right] \\ 0 \text{ otherwise} \end{cases}`",:math:`\mathbb{B}`,
-":math:`F_\text{BAS}^{t,~h,~f}`",Binary variable that indicates the hydo power plant's flow state (for turbine and pump),:math:`\mathbb{B}`,
-":math:`P_\text{CAL}^{t,~h}`",Electrical power produced by a pump,:math:`\mathbb{R}^{+}`,:math:`\mathrm{MW}`
-":math:`Q_\text{CAL}^{t,~h}`",Amount of water that flows through a pump,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
-":math:`P_\text{S}^{t,~h,~s_h}`",Amount of water that flows through a turbine for each hydro powerplant state,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
-":math:`Q_\text{S}^{t,~h,~s_h}`",Amount of water that flows through a pump for each hydro powerplant state,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
+":math:`S_\text{FLOW}^{t,~h,~s_q}`",Binary variable that indicates the hydo power plant's flow state (for turbine and pump),:math:`\mathbb{B}`,
+":math:`P_\text{CAL}^{t,~h,~s_q}`",Calculated electrical power produced (or consumed) by a hydropower plant for each flow state,:math:`\mathbb{R}^{+}`,:math:`\mathrm{MW}`
+":math:`Q_\text{CAL}^{t,~h,~s_q}`",Calculated  water that flows through a hydropower plant for each flow state,:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
+":math:`P_\text{S}^{t,~h,~s_q}`","Real electrical power produced (or consumed) by a hydropower plant for each flow state: :math:`P_\text{S}^{t,~h,~s_q} = \begin{cases} P_\text{CAL}^{t,~h,~s_q} \text{ if } S_\text{FLOW}^{t,~h,~s_q} = 1 \\ 0 \text{ otherwise} \end{cases}`",:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`
+":math:`Q_\text{S}^{t,~h,~s_q}`","Real amount of water that flows through a pump for each hydro powerplant state: :math:`Q_\text{S}^{t,~h,~s_q} = \begin{cases} Q_\text{CAL}^{t,~h,~s_q} \text{ if } S_\text{FLOW}^{t,~h,~s_q} = 1 \\ 0",:math:`\mathbb{R}^{+}`,:math:`\mathrm{m}^{3}/\mathrm{s}`

scripts/baseline_optimizazion.py

@@ -22,7 +22,7 @@
 SIMULATION_SETTING = {
     "1": {"quantile": 0, "buffer": 0.2, "powered_volume_enabled": True, "with_penalty": True},
     "2": {"quantile": 0.15, "buffer": 0.3, "powered_volume_enabled": True, "with_penalty": True},
-    "3": {"quantile": 0.15, "buffer": 0.3, "powered_volume_enabled": False, "with_penalty": True},
+    "3": {"quantile": 0.25, "buffer": 0.3, "powered_volume_enabled": False, "with_penalty": True},
     # "4": {"quantile": 0, "buffer": 0.2, "powered_volume_enabled": True, "with_penalty": False},
 }
 output_file_names: dict[str, str] = json.load(open(settings.OUTPUT_FILE_NAMES))

src/pyomo_models/baseline/second_stage/second_stage_pipeline.py

@@ -1,6 +1,7 @@
 import warnings
 from typing import Optional
 from datetime import datetime, timedelta, timezone
+from annotated_types import T
 from patito import col
 import polars as pl
 from polars import col as c
@@ -394,17 +395,18 @@ def solve_model(self):
                 self.log_mip_gap(solution)
                 self.extract_result()
             else:
-                self.infeasible_constraints = check_infeasible_constraints(model=self.model_instance)
-                self.calculated_feasibility()
-                break
+                solved = self.solve_changing_powered_volume_constraint()
+                if not solved:
+                    break
+
         logging.getLogger('pyomo.core').setLevel(logging.WARNING)        
         self.finalizes_results_processing()
 
     def calculated_feasibility(self):
         start_basin_volume = self.start_basin_volume.filter(c("sim_nb") == self.sim_nb)[["B", "start_basin_volume"]]
         remaining_volume = self.remaining_volume.filter(c("sim_nb") == self.sim_nb)[["H", "remaining_volume"]]
         powered_volume = self.powered_volume.filter(c("sim_nb") == self.sim_nb)[["H", "powered_volume"]]
-        volume_buffer = self.powered_volume.filter(c("sim_nb") == self.sim_nb)[["H", "volume_buffer"]]
+        volume_buffer = self.volume_buffer.filter(c("sim_nb") == self.sim_nb)[["H", "volume_buffer"]]
         basin_power_mapping =pl.from_dicts(self.power_performance_table).drop("power_performance").with_columns(
             pl.all().cast(pl.UInt32)
         )
@@ -416,15 +418,49 @@ def calculated_feasibility(self):
             .join(volume_buffer, on="H", how="inner")\
             .join(self.index["water_basin"][["B", "volume_max", "volume_min"]], on="B", how="inner")\
             .with_columns(
-                (c("start_basin_volume") - c("powered_volume") + c("volume_buffer") +
-                pl.when(c("remaining_volume") > 0).then(c("remaining_volume")).otherwise(0)
-                ).alias("pos_boundary"),
-                (c("start_basin_volume") - c("powered_volume") - c("volume_buffer") +
+                (c("start_basin_volume") - c("powered_volume") + c("volume_buffer") -
                 pl.when(c("remaining_volume") < 0).then(c("remaining_volume")).otherwise(0)
+                ).alias("pos_boundary"),
+                (c("start_basin_volume") - c("powered_volume") - c("volume_buffer") -
+                pl.when(c("remaining_volume") > 0).then(c("remaining_volume")).otherwise(0)
                 ).alias("neg_boundary"),
             ).to_dicts()[0]
         )
 
+    def solve_changing_powered_volume_constraint(self) -> bool:
+        if self.powered_volume_enabled == True:
+            self.powered_volume_enabled = False
+            self.generate_model_instance()   
+            solution = self.solver.solve(self.model_instance, tee=self.log_solver_info)
+            if solution["Solver"][0]["Status"] == "ok":
+                self.extract_result()
+            elif solution["Solver"][0]["Status"] == "aborted":
+                self.log_mip_gap(solution)
+                self.extract_result()
+            else:
+                self.infeasible_constraints = check_infeasible_constraints(model=self.model_instance)
+                self.calculated_feasibility()
+                log.warning(f"Second stage optimization problem is infeasible for sim_nb: {self.sim_nb}")
+                return False
+        else:
+            self.infeasible_constraints = check_infeasible_constraints(model=self.model_instance)
+            self.calculated_feasibility()
+            log.warning(f"Second stage optimization problem is infeasible for sim_nb: {self.sim_nb}")
+            return False
+        self.powered_volume_enabled = True
+        return True
+
+    def log_unbounded(self):
+        self.log_book = pl.concat([
+            self.log_book,
+            pl.DataFrame(
+                {
+                    "sim_nb": [self.sim_nb],
+                    "unbounded": [True]
+                }
+            )
+        ], how="diagonal_relaxed")
+
     def log_mip_gap(self, solution):
         self.log_book = pl.concat([
             self.log_book,

src/pyomo_models/baseline/second_stage/sets.py

@@ -1,4 +1,4 @@
-"""
+r"""
 The sets :math:`T`, :math:`S_H` and :math:`S_B` are arranged in a specific order such that indexing a variable with 
 the first or last element of a set (i.e :math:`t^{0}` and :math:`t^{\\text{END}} \in T`) corresponds respectively to 
 the lowest/first or highest/final element of that variable.
@@ -23,6 +23,8 @@
 
 :math:`S_{BH} \in \{b, h\, s_b, s_h\}=\{(1, ~1, ~1, ~1),~(1, ~1, ~2, ~2),~(1, ~1, ~3, ~3)\}` 
 
+
+
 """
 
 import pyomo.environ as pyo
@@ -35,11 +37,12 @@ def baseline_sets(model):
     # index gathering the state per basin and the hydro powerplants
     model.S_B = pyo.Set(model.B)
     model.S_H = pyo.Set(model.H)
-    model.F = pyo.Set(model.H)
     # index gathering the state of every basin and hydro powerplants
     model.BS = pyo.Set(dimen=2, initialize=lambda model: [(b, s_b) for b in model.B for s_b in model.S_B[b]])
     model.HS = pyo.Set(dimen=2, initialize=lambda model: [(h, s_h) for h in model.H for s_h in model.S_H[h]])
-    model.HF = pyo.Set(dimen=2, initialize=lambda model: [(h, f) for h in model.H for f in model.F[h]])
+
+    model.S_Q = pyo.Set(model.HS)
+    model.HF = pyo.Set(dimen=3, initialize=lambda model: [(i, j, f) for (i, j) in model.JI for f in model.F[i, j]])
 
     model.HSF = pyo.Set(
         dimen=3,