Fix merging HydroCHIPPs

Author: phumthep

src/pownet/builder/hydro.py

@@ -15,6 +15,7 @@ class HydroUnitBuilder(ComponentBuilder):
     Variables
     ===========================
     - `phydro`: Hydropower output. Unit: MW.
+    - `uhydro`: Hydropower unit status. Unit: binary (0 or 1).
 
     Fixed objective terms
     ===========================
@@ -46,6 +47,10 @@ def __init__(self, model: gp.Model, inputs: SystemInput):
 
         # Constraints
         self.c_link_hydro_pu = gp.tupledict()
+        self.c_link_weekly_hydro_pu = gp.tupledict()
+
+        self.c_hydro_limit_daily = gp.tupledict()
+        self.c_hydro_limit_weekly = gp.tupledict()
 
     def add_variables(self, step_k: int) -> None:
         """Add variables to the model for hydro units.
@@ -112,6 +117,17 @@ def add_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
             units=self.inputs.hydro_unit_node.keys(),
             capacity_df=self.inputs.hydro_capacity,
         )
+
+        # Weekly upper bound
+        self.c_link_weekly_hydro_pu = nondispatch_constr.add_c_link_unit_pu(
+            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,
+        )
+
         # Daily upper bound
         self.c_hydro_limit_daily = nondispatch_constr.add_c_hydro_limit_daily(
             model=self.model,
@@ -122,6 +138,17 @@ def add_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
             hydro_capacity=self.inputs.daily_hydro_capacity,
         )
 
+        # Weekly lower and upper bounds
+        self.c_hydro_limit_weekly = nondispatch_constr.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_min_capacity,
+        )
+
     def update_variables(self, step_k: int) -> None:
         "Hydropower variables do not have time-dependent lower/upper bounds."
         return
@@ -157,6 +184,17 @@ def update_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
             hydro_capacity=self.inputs.daily_hydro_capacity,
         )
 
+        self.model.remove(self.c_hydro_limit_weekly)
+        self.c_hydro_limit_weekly = nondispatch_constr.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_min_capacity,
+        )
+
     def update_daily_hydropower_capacity(
         self, step_k: int, new_capacity: dict[(str, int), float]
     ) -> None:

src/pownet/builder/system.py

@@ -108,6 +108,7 @@ def __init__(self, model: gp.Model, inputs: SystemInput):
         self.c_import_curtail_ess = gp.tupledict()
 
         self.c_daily_hydro_curtail_ess = gp.tupledict()
+        self.c_weekly_hydro_curtail_ess = gp.tupledict()
 
     def add_variables(self, step_k: int) -> None:
 
@@ -349,6 +350,19 @@ def get_variable_objective_terms(
             * self.c_daily_hydro_curtail_ess.prod(daily_hydro_coeffs)
         )
 
+        # Weekly hydro curtailment penalty
+        weekly_hydro_coeffs = get_marginal_cost_coeff(
+            step_k=step_k,
+            timesteps=self.timesteps,
+            units=self.inputs.weekly_hydro_must_take_units,
+            nondispatch_contracts=self.inputs.nondispatch_contracts,
+            contract_costs=self.inputs.contract_costs,
+        )
+        self.must_take_curtail_penalty_expr += (
+            curtail_cost_factor
+            * self.c_weekly_hydro_curtail_ess.prod(weekly_hydro_coeffs)
+        )
+
         # Solar, wind, and import curtailment penalties
         nondispatch_tuples = [
             (

src/pownet/core/output.py

@@ -71,6 +71,21 @@ def get_hourly_curtailment(
             node_variables["vartype"] == unit_type_map[unit_type]
         ].pivot(columns="node", index="hour", values="value")
 
+    def get_unit_hourly_generation(self, node_variables: pd.DataFrame) -> pd.DataFrame:
+        power_variables = self._get_power_variables(node_variables)
+        hourly_generation = (
+            power_variables[["unit", "value", "hour"]].groupby(["node", "hour"]).sum()
+        )
+        hourly_generation = hourly_generation.reset_index()
+        hourly_generation = hourly_generation.pivot(
+            columns=["hour"], index=["node"]
+        ).T.reset_index(drop=True)
+        # PowNet indexing starts at 1
+        hourly_generation.index += 1
+        hourly_generation.index.name = "Hour"
+
+        return hourly_generation
+
     def get_hourly_generation(self, node_variables: pd.DataFrame) -> pd.DataFrame:
         power_variables = self._get_power_variables(node_variables)
         hourly_generation = (

src/pownet/core/visualizer.py

@@ -42,7 +42,7 @@ def plot_fuelmix_bar(
         fig, ax = plt.subplots(figsize=(8, 5))
 
         dispatch.plot.bar(
-            stacked=True, ax=ax, linewidth=0, color=self.fuel_color_map, legend=True
+            stacked=True, ax=ax, linewidth=0, color=self.fuel_color_map, legend=False
         )
         ax.plot(
             range(0, total_timesteps),
@@ -72,7 +72,6 @@ def plot_fuelmix_bar(
                 bbox_inches="tight",
                 dpi=350,
             )
-        plt.tight_layout()
         plt.show()
 
     def plot_fuelmix_area(

src/pownet/input.py

@@ -499,8 +499,10 @@ def _load_hydropower(self) -> None:
                 )
             )
 
-        if os.path.exists(os.path.join(self.model_dir, "hydro_capacity_min.csv")):
-            self.hydro_capacity_min = pd.read_csv(os.path.join(self.model_dir, "hydro_capacity_min.csv"))
+        if os.path.exists(os.path.join(self.model_dir, "hydro_min_capacity.csv")):
+            self.hydro_min_capacity = pd.read_csv(
+                os.path.join(self.model_dir, "hydro_min_capacity.csv")
+            )
 
         # Check that the names do not repeat across different types
         repeated_units = set(self.hydro_unit_node.keys()).intersection(
@@ -521,9 +523,9 @@ def _load_hydropower(self) -> None:
             )
 
         # Check that the names do not repeat across different types
-        repeated_units_daily_weekly = set(self.daily_hydro_unit_node.keys()).intersection(
-            self.weekly_hydro_unit_node.keys()
-        )
+        repeated_units_daily_weekly = set(
+            self.daily_hydro_unit_node.keys()
+        ).intersection(self.weekly_hydro_unit_node.keys())
         if repeated_units_daily_weekly:
             raise ValueError(
                 f"PowNet: Found hydropower units to formulate with both daily and weekly formulations: {repeated_units_daily_weekly}"
@@ -783,9 +785,11 @@ def load_data(self):
         # List of units
         #################
         self.thermal_units = list(self.thermal_unit_node.keys())
-        self.hydro_units = (list(self.hydro_unit_node.keys())
-                            + list(self.daily_hydro_unit_node.keys())
-                            + list(self.weekly_hydro_unit_node.keys()))
+        self.hydro_units = (
+            list(self.hydro_unit_node.keys())
+            + list(self.daily_hydro_unit_node.keys())
+            + list(self.weekly_hydro_unit_node.keys())
+        )
         self.solar_units = list(self.solar_unit_node.keys())
         self.wind_units = list(self.wind_unit_node.keys())
         self.import_units = list(self.import_unit_node.keys())

src/pownet/optim_model/constraints/nondispatch_constr.py

@@ -142,3 +142,69 @@ def add_c_hydro_limit_daily_dict(
                 name=cname,
             )
     return constraints
+
+
+def add_c_hydro_limit_weekly(
+    model: gp.Model,
+    phydro: gp.tupledict,
+    step_k: int,
+    sim_horizon: int,
+    hydro_units: list,
+    hydro_capacity: pd.DataFrame,
+    hydro_capacity_min: pd.DataFrame,
+) -> gp.tupledict:
+    """
+    Defines the weekly limit of hydro generation. Assumes that a certain amount of water is
+    available for hydropower generation each day. In this case, the dataframe hydro_capacity has a length of 365 days
+    instead of 8760 hours.
+
+    Args:
+        model (gp.Model): The optimization model
+        phydro (gp.tupledict): The power output of hydro units
+        step_k (int): The current iteration
+        sim_horizon (int): The simulation horizon
+        hydro_units (list): The list of hydro units
+        hydro_capacity (pd.DataFrame): The capacity of the hydro unit
+
+    Returns:
+        gp.tupledict: The constraints for the weekly hydro limit
+
+    Raises:
+        ValueError: If the simulation horizon is not divisible by 24
+
+    """
+    # When formulating with weekly hydropower, sim_horizon must be divisible
+    # by 168 because the hydro_capacity is weekly.
+
+    if sim_horizon % 168 != 0:
+        raise ValueError(
+            "The simulation horizon must be divisible by 168 when using weekly hydropower capacity."
+        )
+    constraints = gp.tupledict()
+    max_week = sim_horizon // 168
+    for week in range(step_k, step_k + max_week):
+        for hydro_unit in hydro_units:
+            cname = f"hydro_limit_weekly_ub[{hydro_unit},{week}]"
+            cname_min = f"hydro_limit_weekly_lb[{hydro_unit},{week}]"
+            current_week = week - step_k + 1
+
+            # Upper bound constraint
+            constraints[cname] = model.addConstr(
+                gp.quicksum(
+                    phydro[hydro_unit, t]
+                    for t in range(1 + (current_week - 1) * 168, current_week * 168 + 1)
+                )
+                <= hydro_capacity.loc[week, hydro_unit],
+                name=cname,
+            )
+
+            # Lower bound constraint
+            constraints[cname_min] = model.addConstr(
+                gp.quicksum(
+                    phydro[hydro_unit, t]
+                    for t in range(1 + (current_week - 1) * 168, current_week * 168 + 1)
+                )
+                >= hydro_capacity_min.loc[week, hydro_unit],
+                name=cname_min,
+            )
+    return constraints