feat(ramping): hydro ramping

Author: MPavicevic

scripts/run_hydrochipps.py

@@ -23,6 +23,7 @@
     - HydroCHIPPs (DoE-sponsored HydroWIRES project)
 """
 import os
+import pandas as pd
 import numpy as np
 from datetime import datetime
 from pownet.core import (
@@ -114,9 +115,10 @@ def update_contract_costs(inputs, supplier, new_values, start_day):
     weekly_prices = np.round(np.random.uniform(5, 240, 168), 1)
     weekly_contract_costs.append(weekly_prices)
 
-#%%########################################
-#     Input Processing and Simulation Loop
-###########################################
+####################################
+# Input Processing and Simulation Loop
+#####################################
+
 # Optional: re-process input data from scratch
 if to_process_inputs:
     data_processor = DataProcessor(
@@ -140,6 +142,7 @@ def update_contract_costs(inputs, supplier, new_values, start_day):
     line_capacity_factor=1,
     line_loss_factor=0,
 )
+
 inputs.load_and_check_data()
 
 # Dictionaries to store results and models
@@ -154,27 +157,31 @@ def update_contract_costs(inputs, supplier, new_values, start_day):
     update_weekly_hydro_capacity(inputs.weekly_hydro_capacity, start_day, first_values)
 
     # === Apply updated weekly prices for this day ===
-    update_contract_costs(inputs, 'supplier', weekly_contract_costs[start_day-1], start_day)
-
-    # Optional: push updates into the model (uncomment if needed)
-    # inputs.update_capacity()
+    update_contract_costs(inputs, 'supplier', weekly_contract_costs[start_day - 1], start_day)
 
     # Initialize model and record-keeping objects
     model_builder = ModelBuilder(inputs)
     record = SystemRecord(inputs)
+
     build_times = []
     opt_times = []
     objvals = []
 
-    # Get initial condition state (e.g., generator status)
-    init_conditions = create_init_condition(inputs.thermal_units)
+    # ========================================================================
+    # CRITICAL FIX: Initialize with proper hydro initial conditions
+    # ========================================================================
+    init_conditions = create_init_condition(
+        thermal_units=inputs.thermal_units,
+        storage_units=inputs.storage_units,
+        ess_max_capacity=inputs.ess_max_capacity,
+        hydro_units=inputs.hydro_units,
+    )
 
     # === Inner loop: run multiple steps from this starting day ===
     if steps_to_run is None:
-        steps_to_run = 10  # 365 - (sim_horizon // 24 - 1)
+        steps_to_run = 10
 
-    #todo:select the start date
-    for step_k in range(start_day, start_day+steps_to_run-1):
+    for step_k in range(start_day, start_day + steps_to_run - 1):
         start_time = datetime.now()
 
         # Build or update the model for the current step
@@ -188,6 +195,7 @@ def update_contract_costs(inputs, supplier, new_values, start_day):
                 step_k=step_k,
                 init_conds=init_conditions,
             )
+
         build_times.append((datetime.now() - start_time).total_seconds())
 
         # Solve the model using HiGHS solver
@@ -206,9 +214,54 @@ def update_contract_costs(inputs, supplier, new_values, start_day):
             runtime=power_system_model[start_day].get_runtime(),
             objval=power_system_model[start_day].get_objval(),
             solution=power_system_model[start_day].get_solution(),
-            # lmp=power_system_model[start_day].solve_for_lmp(),
             step_k=step_k,
         )
+
+        # ========================================================================
+        # CRITICAL FIX: Extract hydro dispatch at t=24 for next iteration
+        # ========================================================================
+        if step_k == start_day:
+            # Get the solution DataFrame
+            first_solution = power_system_model[start_day].get_solution()
+
+            # Initialize dictionary to store hydro dispatch at t=24
+            initial_phydro = {}
+
+            # Iterate through all variables to find phydro at t=24
+            for idx, row in first_solution.iterrows():
+                varname = row['varname']
+
+                # Check if this is a phydro variable
+                if varname.startswith('phydro['):
+                    # Extract unit and timestep using string operations
+                    # Format: phydro[unit, timestep]
+                    try:
+                        # Remove 'phydro[' prefix and ']' suffix
+                        content = varname[7:-1]  # Skip 'phydro[' and ']'
+
+                        # Split by comma
+                        parts = content.split(', ')
+
+                        if len(parts) == 2:
+                            unit = parts[0].strip()
+                            timestep = int(parts[1].strip())
+
+                            # Only keep values at t=24
+                            if timestep == 24:
+                                initial_phydro[unit] = row['value']
+                    except (ValueError, IndexError) as e:
+                        print(f"⚠️  Warning: Could not parse varname '{varname}': {e}")
+                        continue
+
+            # Update init_conditions if we found any hydro dispatch values
+            if initial_phydro:
+                init_conditions['initial_phydro'] = initial_phydro
+                print(f"✓ Initialized hydro ramping for {len(initial_phydro)} units at step {step_k}")
+                print(f"  Sample values: {dict(list(initial_phydro.items())[:3])}")
+            else:
+                print(f"⚠️  Warning: No hydro dispatch found at t=24 for step {step_k}")
+
+        # Update initial conditions from record for subsequent iterations
         init_conditions = record.get_init_conds()
 
         #########################################
@@ -232,5 +285,128 @@ def update_contract_costs(inputs, supplier, new_values, start_day):
         os.makedirs(output_folder, exist_ok=True)
         export_results[start_day].to_csv(output_folder + f'results_day_{str(start_day)}.csv', index=False)
 
+####################################
+##### Quick hydro ramping check ####
+####################################
+
+import matplotlib.pyplot as plt
+
+# Extract hydro dispatch from last solved day
+final_day_results = export_results[start_day + steps_to_run - 2]
+hydro_data = final_day_results[final_day_results['vartype'] == 'phydro']
+
+if not hydro_data.empty:
+    # ========================================================================
+    # CRITICAL FIX: Convert 'hour' to numeric and sort before pivoting
+    # ========================================================================
+    hydro_data = hydro_data.copy()
+    hydro_data['hour'] = pd.to_numeric(hydro_data['hour'], errors='coerce')
+
+    # Drop any rows where hour couldn't be converted
+    hydro_data = hydro_data.dropna(subset=['hour'])
+
+    # Convert to integer
+    hydro_data['hour'] = hydro_data['hour'].astype(int)
+
+    # Sort by hour before pivoting
+    hydro_data = hydro_data.sort_values('hour')
+
+    # Now pivot
+    hydro_pivot = hydro_data.pivot(columns='unit', index='hour', values='value')
+
+    # ========================================================================
+    # VERIFY: Check if index is properly sorted
+    # ========================================================================
+    print(f"\n📊 Hydro dispatch index check:")
+    print(f"   First 5 hours: {hydro_pivot.index[:5].tolist()}")
+    print(f"   Last 5 hours: {hydro_pivot.index[-5:].tolist()}")
+    print(f"   Is sorted: {hydro_pivot.index.is_monotonic_increasing}")
+
+    for unit in hydro_pivot.columns:
+        dispatch = hydro_pivot[unit]
+
+        # ========================================================================
+        # CRITICAL: Ensure dispatch is sorted by index before diff()
+        # ========================================================================
+        dispatch = dispatch.sort_index()
+        ramp_changes = dispatch.diff()
+
+        # ========================================================================
+        # DIAGNOSTIC: Check for ramping violations
+        # ========================================================================
+        ramp_up_limit = inputs.hydro_RU.get(unit, 1e6)
+        ramp_down_limit = inputs.hydro_RD.get(unit, 1e6)
+
+        violations_up = (ramp_changes > ramp_up_limit).sum()
+        violations_down = (ramp_changes < -ramp_down_limit).sum()
+
+        if violations_up > 0 or violations_down > 0:
+            print(f"\n⚠️  RAMPING VIOLATIONS DETECTED for {unit}:")
+            print(f"   Ramp-up violations: {violations_up}")
+            print(f"   Ramp-down violations: {violations_down}")
+            print(f"   Max ramp change: {ramp_changes.max():.2f} MW/h (limit: {ramp_up_limit:.2f})")
+            print(f"   Min ramp change: {ramp_changes.min():.2f} MW/h (limit: -{ramp_down_limit:.2f})")
+
+            # Show the specific violations
+            violation_hours = ramp_changes[
+                (ramp_changes > ramp_up_limit) | (ramp_changes < -ramp_down_limit)
+                ]
+            print(f"   Violation details:\n{violation_hours}")
+
+        # ========================================================================
+        # PLOTTING
+        # ========================================================================
+        fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(14, 8), sharex=True)
+
+        # Dispatch
+        ax1.plot(dispatch.index, dispatch.values, 'b-', linewidth=2, marker='o', markersize=3)
+        ax1.set_ylabel('Power (MW)', fontsize=12)
+        ax1.set_title(f'{unit} - Dispatch and Ramping', fontsize=14, weight='bold')
+        ax1.grid(True, alpha=0.3)
+
+        # Ramp changes with proper coloring
+        colors = []
+        for c in ramp_changes.values[1:]:  # Skip first NaN
+            if abs(c) > ramp_up_limit:
+                colors.append('red')
+            else:
+                colors.append('green')
+
+        # Plot bars
+        ax2.bar(ramp_changes.index[1:], ramp_changes.values[1:],
+                color=colors, alpha=0.7, edgecolor='black', linewidth=0.5)
+
+        # Add limit lines
+        ax2.axhline(ramp_up_limit, color='darkred', linestyle='--', linewidth=2,
+                    label=f'Ramp-up limit ({ramp_up_limit:.0f} MW/h)')
+        ax2.axhline(-ramp_down_limit, color='darkred', linestyle='--', linewidth=2,
+                    label=f'Ramp-down limit ({ramp_down_limit:.0f} MW/h)')
+        ax2.axhline(0, color='black', linestyle='-', linewidth=0.8)
+
+        ax2.set_xlabel('Hour', fontsize=12)
+        ax2.set_ylabel('Ramp (MW/h)', fontsize=12)
+        ax2.set_title('Hour-to-Hour Changes', fontsize=12)
+        ax2.legend(loc='upper right')
+        ax2.grid(True, alpha=0.3, axis='y')
+
+        # Add text annotation for violations
+        if violations_up > 0 or violations_down > 0:
+            ax2.text(0.02, 0.98,
+                     f'⚠️ Violations: {violations_up + violations_down} ({violations_up} up, {violations_down} down)',
+                     transform=ax2.transAxes, fontsize=11, weight='bold',
+                     verticalalignment='top',
+                     bbox=dict(boxstyle='round', facecolor='yellow', alpha=0.8))
+        else:
+            ax2.text(0.02, 0.98,
+                     '✓ All ramps within limits',
+                     transform=ax2.transAxes, fontsize=11, weight='bold',
+                     verticalalignment='top',
+                     bbox=dict(boxstyle='round', facecolor='lightgreen', alpha=0.8))
+
+        plt.tight_layout()
+        plt.savefig(f'{output_folder}/{unit}_ramping_check.png',
+                    dpi=300, bbox_inches='tight')
+        plt.show()
+
 # for i in range(1,6):
 #     print(export_results[i][(export_results[i]['vartype'] == 'phydro') & (export_results[i]['unit'] == 'Barkley')]['value'].sum())

src/pownet/builder/hydro.py

@@ -64,6 +64,11 @@ def __init__(self, model: gp.Model, inputs: SystemInput):
 
         self.c_link_hydro_pu = gp.tupledict()
 
+        self.c_hydro_ramp_up_init = gp.tupledict()
+        self.c_hydro_ramp_up = gp.tupledict()
+        self.c_hydro_ramp_down_init = gp.tupledict()
+        self.c_hydro_ramp_down = gp.tupledict()
+
     def add_variables(self, step_k: int) -> None:
         """Add variables to the model for hydro units.
 
@@ -129,6 +134,7 @@ def add_variables(self, step_k: int) -> None:
                 name="uhydro",
             )
 
+#TODO add fixed ramping cost term to the objective
     def get_fixed_objective_terms(self) -> gp.LinExpr:
         """Hydropower units have no fixed objective terms."""
         return self.total_fixed_objective_expr
@@ -155,6 +161,7 @@ def get_variable_objective_terms(self, step_k: int) -> gp.LinExpr:
 
         return self.total_energy_cost_expr
 
+# TODO check if you can just add the ramp constr here
     def add_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
         # Hourly upper bound
         # Limited by contracted capacity
@@ -217,6 +224,39 @@ def add_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
                 contracted_capacity_dict=self.inputs.hydro_contracted_capacity,
             )
 
+        # Hydro ramping constraints
+        self.c_hydro_ramp_up_init = nondispatch_constr.add_c_hydro_ramp_up_init(
+            model=self.model,
+            phydro=self.phydro,
+            hydro_units=self.inputs.hydro_units,
+            initial_phydro=init_conds.get("initial_phydro", {}),
+            hydro_RU=self.inputs.hydro_RU,
+        )
+
+        self.c_hydro_ramp_up = nondispatch_constr.add_c_hydro_ramp_up(
+            model=self.model,
+            phydro=self.phydro,
+            sim_horizon=self.inputs.sim_horizon,
+            hydro_units=self.inputs.hydro_units,
+            hydro_RU=self.inputs.hydro_RU,
+        )
+
+        self.c_hydro_ramp_down_init = nondispatch_constr.add_c_hydro_ramp_down_init(
+            model=self.model,
+            phydro=self.phydro,
+            hydro_units=self.inputs.hydro_units,
+            initial_phydro=init_conds.get("initial_phydro", {}),
+            hydro_RD=self.inputs.hydro_RD,
+        )
+
+        self.c_hydro_ramp_down = nondispatch_constr.add_c_hydro_ramp_down(
+            model=self.model,
+            phydro=self.phydro,
+            sim_horizon=self.inputs.sim_horizon,
+            hydro_units=self.inputs.hydro_units,
+            hydro_RD=self.inputs.hydro_RD,
+        )
+
     def update_variables(self, step_k: int) -> None:
         "Some hydropower units have hourly upper bounds."
         # Update the time-dependent upper bound of the variable
@@ -226,6 +266,7 @@ def update_variables(self, step_k: int) -> None:
             capacity_df=self.inputs.hydro_capacity,
         )
 
+
     def update_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
         """Update constraints for hydro units.
         Args:
@@ -268,6 +309,27 @@ def update_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
                 hydro_capacity_min=self.inputs.hydro_min_capacity,
             )
 
+        # Remove old ramping constraints
+        self.model.remove(self.c_hydro_ramp_up_init)
+        self.model.remove(self.c_hydro_ramp_down_init)
+
+        # Re-add with new initial conditions
+        self.c_hydro_ramp_up_init = nondispatch_constr.add_c_hydro_ramp_up_init(
+            model=self.model,
+            phydro=self.phydro,
+            hydro_units=self.inputs.hydro_units,
+            initial_phydro=init_conds.get("initial_phydro", {}),
+            hydro_RU=self.inputs.hydro_RU,
+        )
+
+        self.c_hydro_ramp_down_init = nondispatch_constr.add_c_hydro_ramp_down_init(
+            model=self.model,
+            phydro=self.phydro,
+            hydro_units=self.inputs.hydro_units,
+            initial_phydro=init_conds.get("initial_phydro", {}),
+            hydro_RD=self.inputs.hydro_RD,
+        )
+
     def update_daily_hydropower_capacity(
         self, step_k: int, new_capacity: dict[tuple[str, int], float]
     ) -> None:

src/pownet/builder/thermal.py

@@ -88,6 +88,7 @@ def __init__(self, model: gp.Model, inputs: SystemInput) -> None:
         self.c_peak_up_bound = gp.tupledict()
         self.c_ramp_down_init = gp.tupledict()
         self.c_ramp_up_init = gp.tupledict()
+        #TODO take those and add to hydro builder
         self.c_ramp_down = gp.tupledict()
         self.c_ramp_up = gp.tupledict()
 
@@ -401,6 +402,7 @@ def update_variables(self, step_k: int) -> None:
         for var_dict in thermal_unit_vars:
             update_var_with_variable_ub(var_dict, step_k, self.thermal_derated_capacity)
 
+
     def update_constraints(self, step_k: int, init_conds: dict, **kwargs) -> None:
         """Update time-dependent constraints:
         - c_link_uvw_init: initial_u is from the previous iteration