Add multi‑commodity support to device scheduler

documentation/dev/setup-and-guidelines.rst

@@ -63,7 +63,7 @@ On Linux and Windows, everything will be installed using Python packages.
 On MacOS, this will install all test dependencies, and locally install the HiGHS solver.
 For this to work, make sure you have `Homebrew <https://brew.sh/>`_ installed.
 
-Besides the HiGHS solver (as the current default), the CBC solver is required for tests as well. See `The install instructions <https://github.com/coin-or/Cbc?tab=readme-ov-file#binaries`_ for more information.
+Besides the HiGHS solver (as the current default), the CBC solver is required for tests as well. See `The install instructions <https://github.com/coin-or/Cbc?tab=readme-ov-file#binaries>`_ for more information.
 
 Configuration
 ^^^^^^^^^^^^^

flexmeasures/data/models/planning/__init__.py

@@ -288,12 +288,25 @@ 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):
+        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)
+            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)
@@ -386,6 +399,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:
@@ -435,12 +452,25 @@ def to_frame(self) -> pd.DataFrame:
             ],
             axis=1,
         )
-        # map device → device_group
+        # device_group
         if self.device is not None:
             df["device_group"] = map_device_to_group(self.device, self.device_group)
         else:
             df["device_group"] = 0
 
+        # 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]
@@ -579,45 +593,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):
@@ -726,6 +725,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)