Combined Heat and Power (CHP)

documentation/features/scheduling.rst

@@ -5,8 +5,8 @@ Scheduling
 
 Scheduling is the main value-drive of FlexMeasures. We have two major types of schedulers built-in, for storage devices (usually batteries or hot water storage) and processes (usually in industry).
 
-FlexMeasures computes schedules for energy systems that consist of multiple devices that consume and/or produce electricity.
-We model a device as an asset with a power sensor, and compute schedules only for flexible devices, while taking into account inflexible devices.
+FlexMeasures computes schedules for energy systems that consist of multiple devices that consume and/or produce a commodity (e.g. electricity or gas).
+We model a device as an asset with a consumption/production sensor recording power values, and compute schedules only for flexible devices, while taking into account inflexible devices.
 
 .. contents::
     :local:
@@ -58,6 +58,9 @@ And if the asset belongs to a larger system (a hierarchy of assets), the schedul
    * - Field
      - Example value
      - Description
+   * - ``commodity``
+     - |COMMODITY_FLEX_CONTEXT.example|
+     - .. include:: ../_autodoc/COMMODITY_FLEX_CONTEXT.rst
    * - ``inflexible-device-sensors``
      - |INFLEXIBLE_DEVICE_SENSORS.example|
      - .. include:: ../_autodoc/INFLEXIBLE_DEVICE_SENSORS.rst
@@ -70,9 +73,6 @@ 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
@@ -187,8 +187,8 @@ For more details on the possible formats for field values, see :ref:`variable_qu
      - Example value
      - Description
    * - ``commodity``
-     - |COMMODITY.example|
-     - .. include:: ../_autodoc/COMMODITY.rst
+     - |COMMODITY_FLEX_MODEL.example|
+     - .. include:: ../_autodoc/COMMODITY_FLEX_MODEL.rst
    * - ``consumption``
      - |CONSUMPTION.example|
      - .. include:: ../_autodoc/CONSUMPTION.rst

flexmeasures/data/models/planning/__init__.py

@@ -100,6 +100,70 @@ def _build_stock_groups(flex_model: list[dict]) -> dict:
 
         return dict(groups)
 
+    @staticmethod
+    def _build_coupling_groups(
+        flex_model: list[dict],
+    ) -> dict[str, list[tuple[int, float]]]:
+        """Build coupling groups from the 'coupling' and 'coupling_coefficient' fields
+        of each device model.
+
+        Devices sharing the same coupling name form a coupling group.
+        The optimization model introduces a decision variable ``alpha`` per group per time
+        step, and constrains every device by ``P[d] == coeff_d * alpha``.
+
+        Coupling coefficients in flex-models are user-facing positive magnitudes.
+        The internal sign is inferred from directional capacities:
+
+        - ``consumption_capacity == 0`` -> output device -> internally negative coefficient
+        - ``production_capacity == 0`` -> input device -> internally positive coefficient
+
+        If neither direction is explicitly blocked, the coefficient stays positive.
+
+        Example — a CHP with 50% heat efficiency and 30% power efficiency:
+
+            [
+                {"coupling": "chp", "coupling_coefficient": 1.0},  # gas input (alpha = P_gas)
+                {"coupling": "chp", "coupling_coefficient": 0.5},  # heat output (50% of gas)
+                {"coupling": "chp", "coupling_coefficient": 0.3},  # power output (30% of gas)
+            ]
+
+        :param flex_model: List of deserialized device flex-model dicts.
+        :returns: Mapping from coupling-group name to a list of
+                  ``(device_index, internal_signed_coefficient)`` tuples suitable for
+                  passing to ``device_scheduler(coupling_groups=...)``. Returns an empty dict
+                  when no device defines a ``coupling`` field.
+        """
+
+        def _is_zero_capacity(value: Any) -> bool:
+            """Return True if the capacity value is numerically zero."""
+
+            if value is None:
+                return False
+
+            # Pint quantities expose ``magnitude``.
+            magnitude = getattr(value, "magnitude", value)
+            try:
+                return bool(np.isclose(float(magnitude), 0.0))
+            except (TypeError, ValueError):
+                return False
+
+        groups: dict[str, list[tuple[int, float]]] = defaultdict(list)
+        for d, fm in enumerate(flex_model):
+            coupling_name = fm.get("coupling")
+            if coupling_name is None:
+                continue
+            coefficient = abs(float(fm.get("coupling_coefficient", 1.0)))
+
+            is_output = _is_zero_capacity(fm.get("consumption_capacity"))
+            is_input = _is_zero_capacity(fm.get("production_capacity"))
+
+            if is_output and not is_input:
+                coefficient = -coefficient
+
+            groups[coupling_name].append((d, coefficient))
+
+        return dict(groups)
+
     def __init__(
         self,
         sensor: Sensor | None = None,  # deprecated

flexmeasures/data/models/planning/linear_optimization.py

@@ -42,6 +42,7 @@ def device_scheduler(  # noqa C901
     commitments: list[pd.DataFrame] | list[Commitment] | None = None,
     initial_stock: float | list[float] = 0,
     stock_groups: dict[int, list[int]] | None = None,
+    coupling_groups: dict[str, list[tuple[int, float]]] | None = None,
     ems_constraint_groups: list[list[int]] | None = None,
 ) -> tuple[list[pd.Series], float, SolverResults, ConcreteModel]:
     """This generic device scheduler is able to handle an EMS with multiple devices,
@@ -80,6 +81,15 @@ def device_scheduler(  # noqa C901
                                     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.
+    :param coupling_groups:     Hard flow-coupling constraints between devices. Each entry maps a group name to a list of
+                                ``(device_index, coefficient)`` tuples. A decision variable ``alpha`` is introduced per group
+                                per time step and every device ``d`` in the group is constrained by ``P[d, j] == coeff_d * alpha[group, j]``.
+                                Sign convention: positive coefficient for input devices (consuming, positive ``ems_power``),
+                                negative coefficient for output devices (producing, negative ``ems_power``).
+                                Example — a CHP with gas input (d=0, coeff 1.0), heat output (d=1, coeff −0.5) and
+                                power output (d=2, coeff −0.3)::
+
+                                    coupling_groups={"chp": [(0, 1.0), (1, -0.5), (2, -0.3)]}
 
     Potentially deprecated arguments:
         commitment_quantities: amounts of flow specified in commitments (both previously ordered and newly requested)
@@ -131,19 +141,36 @@ def device_scheduler(  # noqa C901
     # map device -> primary stock group (used for per-device stock bounds)
     # and map stock group -> all member devices (used for stock accumulation).
     device_to_group = {}
+    group_to_devices = {}
 
     if stock_groups:
         for g, devices in stock_groups.items():
+            group_to_devices[g] = list(devices)
             for d in devices:
-                device_to_group[d] = g
-        # For devices not in any stock group (e.g., inflexible devices),
-        # map them to themselves so they're treated as individual groups
+                # Keep first assignment as the primary group. A device can still
+                # participate in multiple groups via ``group_to_devices``.
+                if d not in device_to_group:
+                    device_to_group[d] = g
+        # Devices not in any stock group are treated as single-device groups.
         for d in range(len(device_constraints)):
             if d not in device_to_group:
-                device_to_group[d] = d
+                g = f"_device_{d}"
+                device_to_group[d] = g
+                group_to_devices[g] = [d]
     else:
         for d in range(len(device_constraints)):
-            device_to_group[d] = d
+            g = f"_device_{d}"
+            device_to_group[d] = g
+            group_to_devices[g] = [d]
+
+    # Collect (group_index, device_index, coefficient) triples for coupling constraints.
+    # Each device in each group will be constrained: P[d, j] == coeff * alpha[group, j]
+    # where alpha is a free variable representing the common normalised flow.
+    coupling_device_specs: list[tuple[int, int, float]] = []
+    if coupling_groups:
+        for g_idx, (_group_name, members) in enumerate(coupling_groups.items()):
+            for d_idx, coeff in members:
+                coupling_device_specs.append((g_idx, d_idx, coeff))
 
     # Move commitments from old structure to new
     if commitments is None:
@@ -585,7 +612,7 @@ def _get_stock_change(m, d, j):
         group = device_to_group[d]
 
         # all devices belonging to this stock
-        devices = [dev for dev, g in device_to_group.items() if g == group]
+        devices = group_to_devices[group]
 
         # initial stock
         if isinstance(initial_stock, list):
@@ -780,6 +807,29 @@ def device_derivative_equalities(m, d, j):
         model.d, model.j, rule=device_derivative_equalities
     )
 
+    if coupling_device_specs:
+        n_coupling_groups = len(coupling_groups)
+
+        # One free variable per group per time step: the common normalised flow.
+        model.coupling_group_range = RangeSet(0, n_coupling_groups - 1)
+        model.coupling_alpha = Var(model.coupling_group_range, model.j, domain=Reals)
+
+        model.coupling_device_range = RangeSet(0, len(coupling_device_specs) - 1)
+
+        def flow_coupling_rule(m, c, j):
+            """Enforce P[d, j] == coeff * alpha[group, j] for each coupled device.
+
+            This pins every device's flow to the same normalised level ``alpha``,
+            scaled by its coupling coefficient. The coefficient sign indicates direction:
+            positive for inputs (consuming), negative for outputs (producing).
+            """
+            g, d, coeff = coupling_device_specs[c]
+            return m.ems_power[d, j] == coeff * m.coupling_alpha[g, j]
+
+        model.flow_coupling_constraints = Constraint(
+            model.coupling_device_range, model.j, rule=flow_coupling_rule
+        )
+
     # Add objective
     def cost_function(m):
         costs = 0

flexmeasures/data/models/planning/storage.py

@@ -3,7 +3,8 @@
 import re
 import copy
 from datetime import datetime, timedelta
-from typing import Type
+from itertools import chain
+from typing import Any, Type
 
 import pandas as pd
 import numpy as np
@@ -210,6 +211,10 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
         # This ensures the mapping aligns with the device indices
         self.stock_groups = self._build_stock_groups(device_models)
 
+        # Build coupling_groups from the 'coupling' and 'coupling_coefficient' fields
+        # of each device model. Devices sharing the same coupling name form a group.
+        self.coupling_groups = self._build_coupling_groups(device_models)
+
         # List the asset(s) and sensor(s) being scheduled
         if self.asset is not None:
             if not isinstance(self.flex_model, list):
@@ -286,7 +291,7 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
         ]
 
         # Get info from flex-context
-        inflexible_device_sensors = self.flex_context.get(
+        inflexible_device_sensors: list[Sensor] = self.flex_context.get(
             "inflexible_device_sensors", []
         )
 
@@ -327,22 +332,54 @@ def device_list_series(
         ) -> pd.Series:
             return pd.Series([tuple(devices)] * len(index), index=index, name="device")
 
+        # Set up a mapping between commodities and a list of enumerated flexible and inflexible devices
+        # 1. The enumeration starts with all flexible devices in the order that they are given in the flex-model
+        # 2. The enumeration continues with the inflexible devices referenced in the top-level flex-context, in the order given
+        # 3. Then the enumeration goes through the commodities under the commodity_contexts field, in the order given,
+        #    extending the enumeration with the inflexible devices referenced in these commodity contexts.
         commodity_to_devices = {}
+        # Step 1: enumerate the flexible devices
         for d, flex_model_d in enumerate(flex_model):
             commodity = flex_model_d.get("commodity", "electricity")
             commodity_to_devices.setdefault(commodity, []).append(d)
 
-        # inflexible devices are electricity by default
-        number_flexible_devices = len(flex_model)
-        number_inflexible_devices = len(
-            self.flex_context.get("inflexible_device_sensors", [])
-        )
+        # Step 2: enumerate the top-level inflexible devices (electric for backwards compatibility)
         num_flexible_devices = len(flex_model)
         commodity_to_devices["electricity"] += list(
             range(
-                number_flexible_devices,
-                number_flexible_devices + number_inflexible_devices,
+                num_flexible_devices,
+                num_flexible_devices + len(inflexible_device_sensors),
+            )
+        )
+
+        # Step 3: enumerate the inflexible devices per commodity
+        def list_to_map(listing: list, key: Any) -> dict:
+            """Note: the key is retained in the map values."""
+            return {l[key]: l for l in listing}
+
+        # Move inflexible-device-sensors per commodity into the device listing for the commodity
+        commodity_mapping: dict[str, dict] = list_to_map(
+            self.flex_context.get("commodity_contexts", []), key="commodity"
+        )
+        inflexible_devices_per_commodity = {
+            com: con.get("inflexible_device_sensors", [])
+            for com, con in commodity_mapping.items()
+        }
+        num_devices = num_flexible_devices + len(inflexible_device_sensors)
+        for (
+            commodity,
+            commodity_inflexible_device_sensors,
+        ) in inflexible_devices_per_commodity.items():
+            commodity_to_devices[commodity] += list(
+                range(
+                    num_devices,
+                    num_devices + len(commodity_inflexible_device_sensors),
+                )
             )
+            num_devices = num_devices + len(commodity_inflexible_device_sensors)
+
+        inflexible_device_sensors = inflexible_device_sensors + list(
+            chain.from_iterable(inflexible_devices_per_commodity.values())
         )
 
         commodity_contexts = self._get_commodity_contexts()
@@ -365,10 +402,12 @@ def device_list_series(
             if production_price is None:
                 production_price = consumption_price
 
-            if consumption_price is None:
-                raise ValueError(
-                    f"Missing consumption price for commodity '{commodity}'."
-                )
+            # todo: log info statement if commodity has no associated prices
+            # todo: raise if none of the commodities (or maybe electricity specifically) has prices
+            # if consumption_price is None:
+            #     raise ValueError(
+            #         f"Missing consumption price for commodity '{commodity}'."
+            #     )
 
             # Energy prices for this commodity.
             up_deviation_prices = get_continuous_series_sensor_or_quantity(
@@ -379,7 +418,8 @@ def device_list_series(
                 beliefs_before=belief_time,
                 fill_sides=True,
             ).to_frame(name="event_value")
-            ensure_prices_are_not_empty(up_deviation_prices, consumption_price)
+            # todo: see above todo
+            # ensure_prices_are_not_empty(up_deviation_prices, consumption_price)
 
             down_deviation_prices = get_continuous_series_sensor_or_quantity(
                 variable_quantity=production_price,
@@ -389,7 +429,8 @@ def device_list_series(
                 beliefs_before=belief_time,
                 fill_sides=True,
             ).to_frame(name="event_value")
-            ensure_prices_are_not_empty(down_deviation_prices, production_price)
+            # todo: see above todo
+            # ensure_prices_are_not_empty(down_deviation_prices, production_price)
 
             price_frames_by_commodity[commodity] = up_deviation_prices
 
@@ -2118,10 +2159,11 @@ def compute(self, skip_validation: bool = False) -> SchedulerOutputType:
         ems_schedule, expected_costs, scheduler_results, model = device_scheduler(
             device_constraints=device_constraints,
             ems_constraints=ems_constraints,
-            ems_constraint_groups=self.ems_constraint_groups,
             commitments=commitments,
             initial_stock=initial_stock,
             stock_groups=self.stock_groups,
+            coupling_groups=self.coupling_groups if self.coupling_groups else None,
+            ems_constraint_groups=self.ems_constraint_groups,
         )
         if "infeasible" in (tc := scheduler_results.solver.termination_condition):
             raise InfeasibleProblemException(tc)