Feature/stock commitments

documentation/features/scheduling.rst

@@ -105,6 +105,12 @@ For more details on the possible formats for field values, see :ref:`variable_qu
    * - ``site-peak-production-price``
      - ``"260 EUR/MWh"``
      - Production peaks above the ``site-peak-production`` are penalized against this per-kW price. [#penalty_field]_
+   * - ``soc-minima-breach-price``
+     - ``"6 EUR/kWh"``
+     - Penalty for not meeting ``soc-minima`` defined in the flex-model. [#penalty_field]_ [#soc_breach_prices]_
+   * - ``soc-maxima-breach-price``
+     - ``"6 EUR/kWh"``
+     - Penalty for not meeting ``soc-maxima`` defined in the flex-model. [#penalty_field]_ [#soc_breach_prices]_
 
 .. [#old_sensor_field] The old field only accepted an integer (sensor ID).
 
@@ -118,6 +124,8 @@ For more details on the possible formats for field values, see :ref:`variable_qu
 
 .. [#production] Example: with a connection capacity (``site-power-capacity``) of 1 MVA (apparent power) and a production capacity (``site-production-capacity``) of 400 kW (active power), the scheduler will make sure that the grid inflow doesn't exceed 400 kW.
 
+.. [#soc_breach_prices] The SoC breach prices (e.g. 6 EUR/kWh) to use for the schedule are applied over each time step equal to the sensor resolution. For example, a SoC breach price of 6 EUR/kWh per hour, for scheduling a 5-minute resolution sensor, should be passed as a SoC breach price of :math:`6*5/60 = 0.50` EUR/kWh.
+
 .. note:: If no (symmetric, consumption and production) site capacity is defined (also not as defaults), the scheduler will not enforce any bound on the site power.
           The flexible device can still have its own power limit defined in its flex-model.
 

flexmeasures/data/models/planning/__init__.py

@@ -214,6 +214,8 @@ class Commitment:
     ----------
     name:
         Name of the commitment.
+    device:
+        Device to which the commitment pertains. If None, the commitment pertains to the EMS.
     index:
         Pandas DatetimeIndex defining the time slots to which the commitment applies.
         The index is shared by the group, quantity. upwards_deviation_price and downwards_deviation_price Pandas Series.
@@ -234,6 +236,7 @@ class Commitment:
     """
 
     name: str
+    device: pd.Series = None
     index: pd.DatetimeIndex = field(repr=False, default=None)
     _type: str = field(repr=False, default="each")
     group: pd.Series = field(init=False)
@@ -262,6 +265,8 @@ def __post_init__(self):
                 )
 
         # Force type conversion of repr fields to pd.Series
+        if not isinstance(self.device, pd.Series):
+            self.device = pd.Series(self.device, index=self.index)
         if not isinstance(self.quantity, pd.Series):
             self.quantity = pd.Series(self.quantity, index=self.index)
         if not isinstance(self.upwards_deviation_price, pd.Series):
@@ -284,6 +289,7 @@ def __post_init__(self):
             raise ValueError('Commitment `_type` must be "any" or "each".')
 
         # Name the Series as expected by our device scheduler
+        self.device = self.device.rename("device")
         self.quantity = self.quantity.rename("quantity")
         self.upwards_deviation_price = self.upwards_deviation_price.rename(
             "upwards deviation price"
@@ -301,11 +307,20 @@ def to_frame(self) -> pd.DataFrame:
                 self.upwards_deviation_price,
                 self.downwards_deviation_price,
                 self.group,
+                pd.Series(self.__class__, index=self.index, name="class"),
             ],
             axis=1,
         )
 
 
+class FlowCommitment(Commitment):
+    pass
+
+
+class StockCommitment(Commitment):
+    pass
+
+
 """
 Deprecations
 """

flexmeasures/data/models/planning/linear_optimization.py

@@ -59,11 +59,12 @@ def device_scheduler(  # noqa C901
     :param ems_constraints:     EMS constraints are on an EMS level. Handled constraints (listed by column name):
                                     derivative max: maximum flow
                                     derivative min: minimum flow
-    :param commitments:         Commitments are on an EMS level. Handled parameters (listed by column name):
+    :param commitments:         Commitments are on an EMS level by default. Handled parameters (listed by column name):
                                     quantity:                   for example, 5.5
                                     downwards deviation price:  10.1
                                     upwards deviation price:    10.2
                                     group:                      1 (defaults to the enumerate time step j)
+                                    device:                     0 (corresponds to device d; if not set, commitment is on an EMS level)
     :param initial_stock:       initial stock for each device. Use a list with the same number of devices as device_constraints,
                                 or use a single value to set the initial stock to be the same for all devices.
 
@@ -444,12 +445,73 @@ def device_down_derivative_sign(m, d, j):
     def ems_derivative_bounds(m, j):
         return m.ems_derivative_min[j], sum(m.ems_power[:, j]), m.ems_derivative_max[j]
 
+    def device_stock_commitment_equalities(m, c, j, d):
+        """Couple device stocks to each commitment."""
+        if (
+            "device" not in commitments[c].columns
+            or (commitments[c]["device"] != d).all()
+        ):
+            # Commitment c does not concern device d
+            return Constraint.Skip
+        if not all(cl.__name__ == "StockCommitment" for cl in commitments[c]["class"]):
+            raise NotImplementedError(
+                "FlowCommitment on a device level has not been implemented. Please file a GitHub ticket explaining your use case."
+            )
+        if isinstance(initial_stock, list):
+            initial_stock_d = initial_stock[d]
+        else:
+            initial_stock_d = initial_stock
+
+        stock_changes = [
+            (
+                m.device_power_down[d, k] / m.device_derivative_down_efficiency[d, k]
+                + m.device_power_up[d, k] * m.device_derivative_up_efficiency[d, k]
+                + m.stock_delta[d, k]
+            )
+            for k in range(0, j + 1)
+        ]
+        efficiencies = [m.device_efficiency[d, k] for k in range(0, j + 1)]
+        return (
+            (
+                0
+                if len(commitments[c]) == 1
+                or "upwards deviation price" in commitments[c].columns
+                else None
+            ),
+            # 0 if "upwards deviation price" in commitments[c].columns else None,  # todo: possible simplification
+            m.commitment_quantity[c]
+            + m.commitment_downwards_deviation[c]
+            + m.commitment_upwards_deviation[c]
+            - [
+                stock - initial_stock_d
+                for stock in apply_stock_changes_and_losses(
+                    initial_stock_d, stock_changes, efficiencies
+                )
+            ][-1],
+            (
+                0
+                if len(commitments[c]) == 1
+                or "downwards deviation price" in commitments[c].columns
+                else None
+            ),
+            # 0 if "downwards deviation price" in commitments[c].columns else None,  # todo: possible simplification
+        )
+
     def ems_flow_commitment_equalities(m, c, j):
         """Couple EMS flows (sum over devices) to each commitment.
 
         - Creates an inequality for one-sided commitments.
         - Creates an equality for two-sided commitments and for groups of size 1.
         """
+        if "device" in commitments[c].columns:
+            # Commitment c does not concern EMS
+            return Constraint.Skip
+        if "class" in commitments[c].columns and not all(
+            cl.__name__ == "FlowCommitment" for cl in commitments[c]["class"]
+        ):
+            raise NotImplementedError(
+                "StockCommitment on an EMS level has not been implemented. Please file a GitHub ticket explaining your use case."
+            )
         return (
             (
                 0
@@ -499,6 +561,9 @@ def device_derivative_equalities(m, d, j):
     model.ems_power_commitment_equalities = Constraint(
         model.cj, rule=ems_flow_commitment_equalities
     )
+    model.device_energy_commitment_equalities = Constraint(
+        model.cj, model.d, rule=device_stock_commitment_equalities
+    )
     model.device_power_equalities = Constraint(
         model.d, model.j, rule=device_derivative_equalities
     )

flexmeasures/data/models/planning/storage.py

@@ -11,7 +11,12 @@
 
 
 from flexmeasures import Sensor
-from flexmeasures.data.models.planning import Commitment, Scheduler, SchedulerOutputType
+from flexmeasures.data.models.planning import (
+    FlowCommitment,
+    Scheduler,
+    SchedulerOutputType,
+    StockCommitment,
+)
 from flexmeasures.data.models.planning.linear_optimization import device_scheduler
 from flexmeasures.data.models.planning.utils import (
     get_prices,
@@ -271,7 +276,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
         )
 
         # Set up commitments DataFrame
-        commitment = Commitment(
+        commitment = FlowCommitment(
             name="energy",
             quantity=commitment_quantities,
             upwards_deviation_price=commitment_upwards_deviation_price,
@@ -316,7 +321,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             )
 
             # Set up commitments DataFrame
-            commitment = Commitment(
+            commitment = FlowCommitment(
                 name="consumption peak",
                 quantity=ems_peak_consumption,
                 # positive price because breaching in the upwards (consumption) direction is penalized
@@ -360,7 +365,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             )
 
             # Set up commitments DataFrame
-            commitment = Commitment(
+            commitment = FlowCommitment(
                 name="production peak",
                 quantity=-ems_peak_production,  # production is negative quantity
                 # negative price because peaking in the downwards (production) direction is penalized
@@ -405,7 +410,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             )
 
             # Set up commitments DataFrame to penalize any breach
-            commitment = Commitment(
+            commitment = FlowCommitment(
                 name="any consumption breach",
                 quantity=ems_consumption_capacity,
                 # positive price because breaching in the upwards (consumption) direction is penalized
@@ -416,7 +421,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             commitments.append(commitment)
 
             # Set up commitments DataFrame to penalize each breach
-            commitment = Commitment(
+            commitment = FlowCommitment(
                 name="all consumption breaches",
                 quantity=ems_consumption_capacity,
                 # positive price because breaching in the upwards (consumption) direction is penalized
@@ -454,7 +459,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             )
 
             # Set up commitments DataFrame to penalize any breach
-            commitment = Commitment(
+            commitment = FlowCommitment(
                 name="any production breach",
                 quantity=ems_production_capacity,
                 # positive price because breaching in the upwards (consumption) direction is penalized
@@ -465,7 +470,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             commitments.append(commitment)
 
             # Set up commitments DataFrame to penalize each breach
-            commitment = Commitment(
+            commitment = FlowCommitment(
                 name="all production breaches",
                 quantity=ems_production_capacity,
                 # positive price because breaching in the upwards (consumption) direction is penalized
@@ -511,6 +516,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
                     as_instantaneous_events=True,
                     resolve_overlaps="first",
                 )
+                # todo: check flex-model for soc_minima_breach_price and soc_maxima_breach_price fields; if these are defined, create a StockCommitment using both prices (if only 1 price is given, still create the commitment, but only penalize one direction)
             if isinstance(soc_minima[d], Sensor):
                 soc_minima[d] = get_continuous_series_sensor_or_quantity(
                     variable_quantity=soc_minima[d],
@@ -522,6 +528,43 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
                     as_instantaneous_events=True,
                     resolve_overlaps="max",
                 )
+            if self.flex_context.get("soc_minima_breach_price", None) is not None:
+                soc_minima_breach_price = self.flex_context.get(
+                    "soc_minima_breach_price"
+                )
+                soc_minima_breach_price = get_continuous_series_sensor_or_quantity(
+                    variable_quantity=soc_minima_breach_price,
+                    actuator=asset,
+                    unit=(
+                        soc_minima_breach_price.unit
+                        if isinstance(soc_minima_breach_price, Sensor)
+                        else (
+                            soc_minima_breach_price[0]["value"].units
+                            if isinstance(soc_minima_breach_price, list)
+                            else str(soc_minima_breach_price.units)
+                        )
+                    ),
+                    query_window=(start, end),
+                    resolution=resolution,
+                    beliefs_before=belief_time,
+                    fallback_attribute="soc-minima-breach-price",
+                    fill_sides=True,
+                )
+                # Set up commitments DataFrame
+                commitment = StockCommitment(
+                    name="soc minima",
+                    quantity=soc_minima[d],
+                    # negative price because breaching in the downwards (shortage) direction is penalized
+                    downwards_deviation_price=-soc_minima_breach_price,
+                    _type="any",
+                    index=index,
+                    device=0,
+                )
+                commitments.append(commitment)
+
+                # soc-minima will become a soft constraint (modelled as stock commitments), so remove hard constraint
+                soc_minima[d] = None
+
             if isinstance(soc_maxima[d], Sensor):
                 soc_maxima[d] = get_continuous_series_sensor_or_quantity(
                     variable_quantity=soc_maxima[d],
@@ -533,6 +576,42 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
                     as_instantaneous_events=True,
                     resolve_overlaps="min",
                 )
+            if self.flex_context.get("soc_maxima_breach_price", None) is not None:
+                soc_maxima_breach_price = self.flex_context.get(
+                    "soc_maxima_breach_price"
+                )
+                soc_maxima_breach_price = get_continuous_series_sensor_or_quantity(
+                    variable_quantity=soc_maxima_breach_price,
+                    actuator=asset,
+                    unit=(
+                        soc_maxima_breach_price.unit
+                        if isinstance(soc_maxima_breach_price, Sensor)
+                        else (
+                            soc_maxima_breach_price[0]["value"].units
+                            if isinstance(soc_maxima_breach_price, list)
+                            else str(soc_maxima_breach_price.units)
+                        )
+                    ),
+                    query_window=(start, end),
+                    resolution=resolution,
+                    beliefs_before=belief_time,
+                    fallback_attribute="soc-maxima-breach-price",
+                    fill_sides=True,
+                )
+                # Set up commitments DataFrame
+                commitment = StockCommitment(
+                    name="soc maxima",
+                    quantity=soc_maxima[d],
+                    # positive price because breaching in the upwards (surplus) direction is penalized
+                    upwards_deviation_price=soc_maxima_breach_price,
+                    _type="any",
+                    index=index,
+                    device=0,
+                )
+                commitments.append(commitment)
+
+                # soc-maxima will become a soft constraint (modelled as stock commitments), so remove hard constraint
+                soc_maxima[d] = None
 
             device_constraints[d] = add_storage_constraints(
                 start,
@@ -1207,7 +1286,7 @@ def add_storage_constraints(
 
     :param start:                       Start of the schedule.
     :param end:                         End of the schedule.
-    :param resolution:                  Timedelta used to resample the forecasts to the resolution of the schedule.
+    :param resolution:                  Timedelta used to resample the constraints to the resolution of the schedule.
     :param soc_at_start:                State of charge at the start time.
     :param soc_targets:                 Exact targets for the state of charge at each time.
     :param soc_maxima:                  Maximum state of charge at each time.

flexmeasures/data/models/planning/tests/test_solver.py

@@ -25,6 +25,7 @@
     get_sensors_from_db,
 )
 from flexmeasures.data.models.planning.utils import initialize_series, initialize_df
+from flexmeasures.data.models.planning import StockCommitment
 from flexmeasures.data.schemas.sensors import TimedEventSchema
 from flexmeasures.utils.calculations import (
     apply_stock_changes_and_losses,
@@ -2652,9 +2653,8 @@ def initialize_combined_commitments(num_devices: int):
             stock_commitment["downwards deviation price"] = -soc_target_penalty
             stock_commitment["upwards deviation price"] = soc_target_penalty
             stock_commitment["group"] = list(range(len(stock_commitment)))
-            # todo: amend test for https://github.com/FlexMeasures/flexmeasures/pull/1300
-            # stock_commitment["device"] = 0
-            # stock_commitment["class"] = StockCommitment
+            stock_commitment["device"] = 0
+            stock_commitment["class"] = StockCommitment
             commitments.append(stock_commitment)
 
         return commitments

flexmeasures/data/models/planning/utils.py

@@ -334,7 +334,7 @@ def idle_after_reaching_target(
 
 def get_quantity_from_attribute(
     entity: Asset | Sensor,
-    attribute: str,
+    attribute: str | None,
     unit: str | ur.Quantity,
 ) -> ur.Quantity:
     """Get the value (in the given unit) of a quantity stored as an entity attribute.
@@ -344,6 +344,9 @@ def get_quantity_from_attribute(
     :param unit:        The unit in which the value should be returned.
     :return:            The retrieved quantity or the provided default.
     """
+    if attribute is None:
+        return np.nan * ur.Quantity(unit)  # at least return result in the desired unit
+
     # Get the default value from the entity attribute
     value: str | float | int = entity.get_attribute(attribute, np.nan)