Big M for load shifting implemented

assume/strategies/naive_strategies.py

@@ -150,7 +150,7 @@ def calculate_bids(
         **kwargs,
     ) -> Orderbook:
         # calculate the optimal operation of the unit
-        unit.determine_optimal_operation_with_flex()
+        unit.determine_optimal_operation_without_flex()
 
         bids = []
         for product in product_tuples:

assume/units/dsm_load_shift.py

@@ -22,6 +22,7 @@ class DSMFlex:
             model,
         ),
     }
+    big_M = 10000000
 
     def __init__(self, components, **kwargs):
         super().__init__(**kwargs)
@@ -74,8 +75,8 @@ def cost_based_flexibility(self, model):
         model.total_cost = pyo.Param(initialize=0.0, mutable=True)
 
         # Variables
-        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
-        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
+        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
+        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
         model.shift_indicator = pyo.Var(model.time_steps, within=pyo.Binary)
 
         @model.Constraint()
@@ -84,22 +85,26 @@ def total_cost_upper_limit(m):
                 m.variable_cost[t] for t in m.time_steps
             ) <= m.total_cost * (1 + (m.cost_tolerance / 100))
 
+        @model.Constraint(model.time_steps)
+        def flex_constraint_upper(m, t):
+            return m.load_shift_pos[t] <= (1 - m.shift_indicator[t]) * self.big_M
+
+        @model.Constraint(model.time_steps)
+        def flex_constraint_lower(m, t):
+            return m.load_shift_neg[t] <= m.shift_indicator[t] * self.big_M
+
         @model.Constraint(model.time_steps)
         def total_power_input_constraint_with_flex(m, t):
             if self.has_electrolyser:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * model.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - model.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == self.model.dsm_blocks["electrolyser"].power_in[t]
                     + self.model.dsm_blocks["eaf"].power_in[t]
                     + self.model.dsm_blocks["dri_plant"].power_in[t]
                 )
             else:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * model.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - model.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == self.model.dsm_blocks["eaf"].power_in[t]
                     + self.model.dsm_blocks["dri_plant"].power_in[t]
                 )
@@ -144,8 +149,8 @@ def get_congestion_indicator(congestion_signal, threshold):
         )
 
         # Variables
-        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
-        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
+        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
+        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
         model.shift_indicator = pyo.Var(model.time_steps, within=pyo.Binary)
 
         # Constraint to manage total cost upper limit with cost tolerance
@@ -155,23 +160,27 @@ def total_cost_upper_limit(m):
                 m.variable_cost[t] for t in m.time_steps
             ) <= m.total_cost * (1 + (m.cost_tolerance / 100))
 
+        @model.Constraint(model.time_steps)
+        def flex_constraint_upper(m, t):
+            return m.load_shift_pos[t] <= (1 - m.shift_indicator[t]) * self.big_M
+
+        @model.Constraint(model.time_steps)
+        def flex_constraint_lower(m, t):
+            return m.load_shift_neg[t] <= m.shift_indicator[t] * self.big_M
+
         # Power input constraint with flexibility based on congestion
         @model.Constraint(model.time_steps)
         def total_power_input_constraint_with_flex(m, t):
             if self.has_electrolyser:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * model.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - model.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == self.model.dsm_blocks["electrolyser"].power_in[t]
                     + self.model.dsm_blocks["eaf"].power_in[t]
                     + self.model.dsm_blocks["dri_plant"].power_in[t]
                 )
             else:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * model.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - model.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == self.model.dsm_blocks["eaf"].power_in[t]
                     + self.model.dsm_blocks["dri_plant"].power_in[t]
                 )
@@ -198,8 +207,8 @@ def peak_load_shifting_flexibility(self, model):
         model.total_cost = pyo.Param(initialize=0.0, mutable=True)
 
         # Variables for load shifting
-        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
-        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
+        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
+        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
         model.shift_indicator = pyo.Var(model.time_steps, within=pyo.Binary)
 
         peak_periods = {
@@ -218,23 +227,27 @@ def total_cost_upper_limit(m):
                 m.variable_cost[t] for t in m.time_steps
             ) <= m.total_cost * (1 + (m.cost_tolerance / 100))
 
+        @model.Constraint(model.time_steps)
+        def flex_constraint_upper(m, t):
+            return m.load_shift_pos[t] <= (1 - m.shift_indicator[t]) * self.big_M
+
+        @model.Constraint(model.time_steps)
+        def flex_constraint_lower(m, t):
+            return m.load_shift_neg[t] <= m.shift_indicator[t] * self.big_M
+
         # Power input constraint with flexibility based on congestion
         @model.Constraint(model.time_steps)
         def total_power_input_constraint_with_peak_shift(m, t):
             if self.has_electrolyser:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * m.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - m.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == m.dsm_blocks["electrolyser"].power_in[t]
                     + m.dsm_blocks["eaf"].power_in[t]
                     + m.dsm_blocks["dri_plant"].power_in[t]
                 )
             else:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * m.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - m.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == m.dsm_blocks["eaf"].power_in[t]
                     + m.dsm_blocks["dri_plant"].power_in[t]
                 )
@@ -286,8 +299,8 @@ def renewable_utilisation(self, model):
         model.total_cost = pyo.Param(initialize=0.0, mutable=True)
 
         # Variables for load flexibility based on renewable intensity
-        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
-        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeReals)
+        model.load_shift_pos = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
+        model.load_shift_neg = pyo.Var(model.time_steps, within=pyo.NonNegativeIntegers)
         model.shift_indicator = pyo.Var(model.time_steps, within=pyo.Binary)
 
         # Constraint to manage total cost upper limit with cost tolerance
@@ -297,23 +310,27 @@ def total_cost_upper_limit(m):
                 m.variable_cost[t] for t in m.time_steps
             ) <= m.total_cost * (1 + (m.cost_tolerance / 100))
 
+        @model.Constraint(model.time_steps)
+        def flex_constraint_upper(m, t):
+            return m.load_shift_pos[t] <= (1 - m.shift_indicator[t]) * self.big_M
+
+        @model.Constraint(model.time_steps)
+        def flex_constraint_lower(m, t):
+            return m.load_shift_neg[t] <= m.shift_indicator[t] * self.big_M
+
         # Power input constraint integrating flexibility
         @model.Constraint(model.time_steps)
         def total_power_input_constraint_flex(m, t):
             if self.has_electrolyser:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * m.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - m.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == m.dsm_blocks["electrolyser"].power_in[t]
                     + m.dsm_blocks["eaf"].power_in[t]
                     + m.dsm_blocks["dri_plant"].power_in[t]
                 )
             else:
                 return (
-                    m.total_power_input[t]
-                    + m.load_shift_pos[t] * m.shift_indicator[t]
-                    - m.load_shift_neg[t] * (1 - m.shift_indicator[t])
+                    m.total_power_input[t] + m.load_shift_pos[t] - m.load_shift_neg[t]
                     == m.dsm_blocks["eaf"].power_in[t]
                     + m.dsm_blocks["dri_plant"].power_in[t]
                 )

assume/units/steel_plant.py

@@ -18,7 +18,7 @@
 from assume.common.utils import get_products_index
 from assume.units.dsm_load_shift import DSMFlex
 
-SOLVERS = ["gurobi", "appsi_highs", "glpk", "cbc", "cplex"]
+SOLVERS = ["appsi_highs", "gurobi", "glpk", "cbc", "cplex"]
 
 logger = logging.getLogger(__name__)
 
@@ -391,8 +391,7 @@ def obj_rule_flex(m):
                 """
 
                 maximise_load_shift = pyo.quicksum(
-                    m.load_shift_neg[t] * (1 - m.shift_indicator[t])
-                    for t in m.time_steps
+                    m.load_shift_pos[t] for t in m.time_steps
                 )
 
                 return maximise_load_shift
@@ -405,9 +404,7 @@ def obj_rule_flex(m):
                 Maximizes the load shift over all time steps.
                 """
                 maximise_load_shift = pyo.quicksum(
-                    m.load_shift_neg[t]
-                    * (1 - m.shift_indicator[t])
-                    * m.congestion_indicator[t]
+                    m.load_shift_neg[t] * m.congestion_indicator[t]
                     for t in m.time_steps
                 )
 
@@ -421,10 +418,7 @@ def obj_rule_flex(m):
                 Maximizes the load shift over all time steps.
                 """
                 maximise_load_shift = pyo.quicksum(
-                    m.load_shift_neg[t]
-                    * (1 - m.shift_indicator[t])
-                    * m.peak_indicator[t]
-                    for t in m.time_steps
+                    m.load_shift_neg[t] * m.peak_indicator[t] for t in m.time_steps
                 )
 
                 return maximise_load_shift
@@ -437,8 +431,7 @@ def obj_rule_flex(m):
                 Maximizes the load increase over all time steps based on renewable surplus.
                 """
                 maximise_renewable_utilisation = pyo.quicksum(
-                    m.load_shift_pos[t] * m.shift_indicator[t] * m.renewable_signal[t]
-                    for t in m.time_steps
+                    m.load_shift_pos[t] * m.renewable_signal[t] for t in m.time_steps
                 )
 
                 return maximise_renewable_utilisation