Add multi‑commodity support to device scheduler

documentation/features/scheduling.rst

@@ -70,6 +70,9 @@ And if the asset belongs to a larger system (a hierarchy of assets), the schedul
    * - ``production-price``
      - |PRODUCTION_PRICE.example|
      - .. include:: ../_autodoc/PRODUCTION_PRICE.rst
+   * - ``gas-price``
+     - |GAS_PRICE.example|
+     - .. include:: ../_autodoc/GAS_PRICE.rst
    * - ``site-power-capacity``
      - |SITE_POWER_CAPACITY.example|
      - .. include:: ../_autodoc/SITE_POWER_CAPACITY.rst
@@ -183,6 +186,9 @@ For more details on the possible formats for field values, see :ref:`variable_qu
    * - Field
      - Example value
      - Description
+   * - ``commodity``
+     - |COMMODITY.example|
+     - .. include:: ../_autodoc/COMMODITY.rst
    * - ``consumption``
      - |CONSUMPTION.example|
      - .. include:: ../_autodoc/CONSUMPTION.rst

flexmeasures/data/models/planning/__init__.py

@@ -289,14 +289,27 @@ class Commitment:
     quantity: pd.Series = 0
     upwards_deviation_price: pd.Series = 0
     downwards_deviation_price: pd.Series = 0
+    commodity: str | pd.Series | None = None
 
     def __post_init__(self):
         # device_group is a device→label lookup table, not a time series;
         # exclude it from automatic time-series index coercion.
+        if (
+            isinstance(self, FlowCommitment)
+            and isinstance(self.commodity, pd.Series)
+            and self.device is not None
+        ):
+            devices = extract_devices(self.device)
+            missing = set(devices) - set(self.commodity.index)
+            if missing:
+                raise ValueError(f"commodity mapping missing for devices: {missing}")
+
         series_attributes = [
             attr
             for attr, _type in self.__annotations__.items()
-            if _type == "pd.Series" and hasattr(self, attr) and attr != "device_group"
+            if _type == "pd.Series"
+            and hasattr(self, attr)
+            and attr not in ("device_group", "commodity")
         ]
         for series_attr in series_attributes:
             val = getattr(self, series_attr)
@@ -389,6 +402,10 @@ def _init_device_group(self):
                 range(len(devices)), index=devices, name="device_group"
             )
         else:
+            if not isinstance(self.device_group, pd.Series):
+                self.device_group = pd.Series(
+                    self.device_group, index=devices, name="device_group"
+                )
             # Validate custom grouping
             missing = set(devices) - set(self.device_group.index)
             if missing:
@@ -450,6 +467,19 @@ def to_frame(self) -> pd.DataFrame:
                 np.nan
             )  # EMS-level: handled by ems_flow_commitment_equalities
 
+        # commodity
+        if getattr(self, "commodity", None) is None:
+            df["commodity"] = None
+        elif isinstance(self.commodity, pd.Series):
+            # commodity is a device→commodity mapping, like device_group
+            if self.device is None:
+                df["commodity"] = None
+            else:
+                df["commodity"] = map_device_to_group(self.device, self.commodity)
+        else:
+            # scalar commodity
+            df["commodity"] = self.commodity
+
         return df
 
 

flexmeasures/data/models/planning/linear_optimization.py

@@ -135,6 +135,20 @@ def device_scheduler(  # noqa C901
             df["group"] = group
             commitments.append(df)
 
+    # commodity → set(device indices)
+    commodity_devices = {}
+
+    for df in commitments:
+        if "commodity" not in df.columns or "device" not in df.columns:
+            continue
+
+        for _, row in df[["commodity", "device"]].dropna().iterrows():
+            devices = row["device"]
+            if not isinstance(devices, (list, tuple, set)):
+                devices = [devices]
+
+            commodity_devices.setdefault(row["commodity"], set()).update(devices)
+
     # Check if commitments have the same time window and resolution as the constraints
     for commitment in commitments:
         start_c = commitment.index.to_pydatetime()[0]
@@ -593,45 +607,30 @@ def device_stock_commitment_equalities(m, c, j, d):
         )
 
     def ems_flow_commitment_equalities(m, c, j):
-        """Couple EMS flows (sum over devices) to each commitment.
+        """Couple EMS flow commitments to device flows, optionally filtered by commodity."""
 
-        - 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
-            and not pd.isnull(commitments[c]["device"]).all()
-        ) or m.commitment_quantity[c, j] == -infinity:
-            # Commitment c does not concern EMS
+        if commitments[c]["class"].iloc[0] != FlowCommitment:
             return Constraint.Skip
-        if (
-            "class" in commitments[c].columns
-            and not (
-                commitments[c]["class"].apply(lambda cl: cl == FlowCommitment)
-            ).all()
-        ):
-            raise NotImplementedError(
-                "StockCommitment on an EMS level has not been implemented. Please file a GitHub ticket explaining your use case."
-            )
+
+        # Legacy behavior: no commodity → sum over all devices
+        if "commodity" not in commitments[c].columns:
+            devices = m.d
+        else:
+            commodity = commitments[c]["commodity"].iloc[0]
+            if pd.isna(commodity):
+                devices = m.d
+            else:
+                devices = commodity_devices.get(commodity, set())
+                if not devices:
+                    return Constraint.Skip
+
         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
+            None,
             m.commitment_quantity[c, j]
             + m.commitment_downwards_deviation[c]
             + m.commitment_upwards_deviation[c]
-            - sum(m.ems_power[:, j]),
-            (
-                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
+            - sum(m.ems_power[d, j] for d in devices),
+            None,
         )
 
     def device_derivative_equalities(m, d, j):
@@ -751,6 +750,22 @@ def cost_function(m):
         )
 
     model.commitment_costs = commitment_costs
+    commodity_costs = {}
+    for c in model.c:
+        commodity = None
+        if "commodity" in commitments[c].columns:
+            commodity = commitments[c]["commodity"].iloc[0]
+        if commodity is None or (isinstance(commodity, float) and np.isnan(commodity)):
+            continue
+
+        cost = value(
+            model.commitment_downwards_deviation[c] * model.down_price[c]
+            + model.commitment_upwards_deviation[c] * model.up_price[c]
+        )
+        commodity_costs[commodity] = commodity_costs.get(commodity, 0) + cost
+
+    model.commodity_costs = commodity_costs
+
     # model.pprint()
     # model.display()
     # print(results.solver.termination_condition)