diff --git a/documentation/features/scheduling.rst b/documentation/features/scheduling.rst
index aec2ce3397..b3acbeefdd 100644
--- a/documentation/features/scheduling.rst
+++ b/documentation/features/scheduling.rst
@@ -196,6 +196,12 @@ For more details on the possible formats for field values, see :ref:`variable_qu
    * - ``commodity``
      - |COMMODITY_FLEX_MODEL.example|
      - .. include:: ../_autodoc/COMMODITY_FLEX_MODEL.rst
+   * - ``coupling``
+     - |COUPLING.example|
+     - .. include:: ../_autodoc/COUPLING.rst
+   * - ``coupling-coefficient``
+     - |COUPLING_COEFFICIENT.example|
+     - .. include:: ../_autodoc/COUPLING_COEFFICIENT.rst
    * - ``consumption``
      - |CONSUMPTION.example|
      - .. include:: ../_autodoc/CONSUMPTION.rst
diff --git a/flexmeasures/data/models/planning/devices.py b/flexmeasures/data/models/planning/devices.py
index b7eb4a98c2..ab70861b71 100644
--- a/flexmeasures/data/models/planning/devices.py
+++ b/flexmeasures/data/models/planning/devices.py
@@ -20,6 +20,7 @@
 
 from __future__ import annotations
 
+import math
 from dataclasses import dataclass, field
 from enum import Enum
 from functools import cached_property
@@ -34,9 +35,11 @@
 class DeviceRole(Enum):
     """The role a flex-model (or flex-context) entry plays in the scheduling problem.
 
-    Extension points (not yet implemented): GROUP (an entry constraining the aggregate
-    power of a set of member devices) and CONVERTER_PORT (a commodity port of a
-    multi-commodity converter).
+    Converter ports (the commodity ports of a multi-commodity converter,
+    such as a CHP unit) are DEVICE entries carrying a ``coupling`` field;
+    see :attr:`DeviceInventory.coupling_groups`.
+    Extension point (not yet implemented):
+    GROUP (an entry constraining the aggregate power of a set of member devices).
     """
 
     #: A schedulable flexible device (usually with a power sensor).
@@ -69,6 +72,12 @@ class FlexDevice:
     #: Key of the stock this device draws from: the id of its state-of-charge sensor, or a unique negative synthetic key for devices without one.
     #: None for inflexible devices.
     stock_key: int | None = None
+    #: Name of the coupling group this device belongs to (converter ports of one converter share a coupling name).
+    #: None for uncoupled devices.
+    coupling: str | None = None
+    #: Signed internal coupling coefficient: positive for input (consuming) ports, negative for output (producing) ports.
+    #: Meaningless (1.0) for uncoupled devices.
+    coupling_coefficient: float = 1.0
 
     @property
     def sensor_id(self) -> int | None:
@@ -147,6 +156,37 @@ def _resolve_stock_key(state_of_charge: Any) -> int | None:
     return key
 
 
+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 math.isclose(float(magnitude), 0.0, abs_tol=1e-08)
+    except (TypeError, ValueError):
+        return False
+
+
+def _resolve_coupling_coefficient(flex_model: dict) -> float:
+    """Resolve a coupled device's internal signed coupling coefficient.
+
+    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.
+    """
+    coefficient = abs(float(flex_model.get("coupling_coefficient", 1.0)))
+    is_output = _is_zero_capacity(flex_model.get("consumption_capacity"))
+    is_input = _is_zero_capacity(flex_model.get("production_capacity"))
+    if is_output and not is_input:
+        coefficient = -coefficient
+    return coefficient
+
+
 #: Flex-model fields that make a device entry (with a state-of-charge sensor)
 #: also carry the SoC parameters of its stock.
 SOC_PARAM_FIELDS = ("soc_at_start", "soc_min", "soc_max", "soc_targets")
@@ -267,6 +307,8 @@ def register_stock_params(stock_key: int, fm: dict) -> None:
                 ),
                 commodity=fm.get("commodity", "electricity"),
                 stock_key=stock_key,
+                coupling=fm.get("coupling"),
+                coupling_coefficient=_resolve_coupling_coefficient(fm),
             )
             inventory.entries.append(device)
             inventory.devices.append(device)
@@ -344,6 +386,27 @@ def stock_params(self, stock_key: int) -> dict | None:
         """Return the flex-model entry holding the SoC parameters of the given stock."""
         return self.stock_entries.get(stock_key)
 
+    @cached_property
+    def coupling_groups(self) -> dict[str, list[tuple[int, float]]]:
+        """Map each coupling-group name to its ports' (device index, signed coefficient) pairs.
+
+        Devices sharing a coupling name are the commodity ports of one converter (e.g. a CHP unit's gas input, heat output and electricity output).
+        The optimization model introduces a decision variable ``alpha`` per group per time step,
+        and constrains every port by ``P[d] == coeff_d * alpha``.
+        The coefficient signs follow the internal convention (see :func:`_resolve_coupling_coefficient`):
+        positive for inputs, negative for outputs.
+        The result is suitable for passing to ``device_scheduler(coupling_groups=...)``;
+        it is empty when no device defines a ``coupling`` field.
+        """
+        groups: dict[str, list[tuple[int, float]]] = {}
+        for device in self.devices:
+            if device.coupling is None:
+                continue
+            groups.setdefault(device.coupling, []).append(
+                (device.index, device.coupling_coefficient)
+            )
+        return groups
+
     @cached_property
     def commodity_to_devices(self) -> dict[str, list[int]]:
         """Map each commodity to its device indices, in canonical solver order."""
diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py
index 12b33fd70c..cff76d9416 100644
--- a/flexmeasures/data/models/planning/linear_optimization.py
+++ b/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:
@@ -550,6 +577,35 @@ def grouped_commitment_equalities(m, c, j, g):
     )
     model.commitment_sign = Var(model.c, domain=Binary, initialize=0)
 
+    # def _get_stock_change(m, d, j):
+    #     """Determine final stock change of device d until time j.
+    #
+    #     Apply conversion efficiencies to conversion from flow to stock change and vice versa,
+    #     and apply storage efficiencies to stock levels from one datetime to the next.
+    #     """
+    #     if isinstance(initial_stock, list):
+    #         # No initial stock defined for inflexible device
+    #         initial_stock_d = initial_stock[d] if d < len(initial_stock) else 0
+    #     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)]
+    #     final_stock_change = [
+    #         stock - initial_stock_d
+    #         for stock in apply_stock_changes_and_losses(
+    #             initial_stock_d, stock_changes, efficiencies
+    #         )
+    #     ][-1]
+    #     return final_stock_change
+
     def _get_stock_change(m, d, j):
         """Determine final stock change of the stock group of device d until time j.
 
@@ -583,7 +639,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):
@@ -778,6 +834,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
diff --git a/flexmeasures/data/models/planning/storage.py b/flexmeasures/data/models/planning/storage.py
index a6686748b8..993cd721cd 100644
--- a/flexmeasures/data/models/planning/storage.py
+++ b/flexmeasures/data/models/planning/storage.py
@@ -148,6 +148,10 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
         # The stock groups' device indices align with the device models
         self.stock_groups = inventory.stock_groups
 
+        # The coupling groups (converter ports sharing a coupling name) also derive from the inventory,
+        # with signed coefficients per canonical device index.
+        self.coupling_groups = inventory.coupling_groups
+
         # List the asset(s) and sensor(s) being scheduled
         sensors: list[Sensor | None] = inventory.power_sensors
         assets: list[Asset | None] = inventory.assets
@@ -2542,6 +2546,7 @@ def compute(self, skip_validation: bool = False) -> SchedulerOutputType:
             commitments=commitments,
             initial_stock=initial_stock,
             stock_groups=self.stock_groups,
+            coupling_groups=self.coupling_groups if self.coupling_groups else None,
         )
         if "infeasible" in (tc := scheduler_results.solver.termination_condition):
             raise InfeasibleProblemException(tc)
diff --git a/flexmeasures/data/models/planning/tests/test_commitments.py b/flexmeasures/data/models/planning/tests/test_commitments.py
index 4275892c53..5eb8620fc3 100644
--- a/flexmeasures/data/models/planning/tests/test_commitments.py
+++ b/flexmeasures/data/models/planning/tests/test_commitments.py
@@ -1556,6 +1556,667 @@ def test_simulation_with_dynamic_consumption_capacity(app, db):
     ), "Electric heater should have one expected partial 80 kW dispatch step before the first cheap-electricity window."
 
 
+def test_chp_coupling():
+    """Test that coupling_groups enforces fixed flow ratios between CHP devices.
+
+    Models a Combined Heat and Power unit with three pure flow devices:
+
+    - d=0  gas input:    can only consume gas          (derivative_min=0)
+    - d=1  heat output:  can only produce heat          (derivative_max=0)
+    - d=2  power output: can only produce electricity   (derivative_max=0)
+
+    The coupling group ``"chp"`` is specified with coefficients
+    ``[(0, 1.0), (1, -0.5), (2, -0.3)]``, introducing a decision variable ``alpha``
+    and enforcing ``P[d] == coeff * alpha`` for each device:
+
+        P_gas   =  1.0 * alpha   (input,  coeff =  1.0)
+        P_heat  = -0.5 * alpha   (output, coeff = -0.5, heat efficiency 50%)
+        P_power = -0.3 * alpha   (output, coeff = -0.3, power efficiency 30%)
+
+    Heat production is forced to exactly 10 kW via ``derivative equals = -10``
+    on device 1. Substituting ``P_heat = -10`` gives ``alpha = 20``, so:
+
+        P_gas   =  20 kW         (gas consumed)
+        P_heat  = -10 kW         (heat produced, forced)
+        P_power =  20 kW * -0.3
+                ≈  -6 kW         (electricity produced)
+
+    """
+    start = pd.Timestamp("2026-01-01T00:00+01:00")
+    end = pd.Timestamp("2026-01-01T04:00+01:00")
+    resolution = pd.Timedelta("1h")
+    index = initialize_index(start=start, end=end, resolution=resolution)
+
+    # d=0: gas input — can only consume (derivative_min=0), capacity 100 kW.
+    # NaN stock bounds mean no cumulative-stock constraint (pure flow device).
+    gas_constraints = pd.DataFrame(
+        {
+            "min": np.nan,
+            "max": np.nan,
+            "equals": np.nan,
+            "derivative min": 0.0,
+            "derivative max": 100.0,
+            "derivative equals": np.nan,
+            "derivative down efficiency": 1.0,
+            "derivative up efficiency": 1.0,
+        },
+        index=index,
+    )
+
+    # d=1: heat output — can only produce (derivative_max=0).
+    # Forced to exactly -10 kW via derivative equals.
+    heat_constraints = pd.DataFrame(
+        {
+            "min": np.nan,
+            "max": np.nan,
+            "equals": np.nan,
+            "derivative min": -100.0,
+            "derivative max": 0.0,
+            "derivative equals": -10.0,
+            "derivative down efficiency": 1.0,
+            "derivative up efficiency": 1.0,
+        },
+        index=index,
+    )
+
+    # d=2: power output — can only produce (derivative_max=0), capacity 100 kW.
+    # Flow is free; the coupling constraint will determine its value.
+    power_constraints = pd.DataFrame(
+        {
+            "min": np.nan,
+            "max": np.nan,
+            "equals": np.nan,
+            "derivative min": -100.0,
+            "derivative max": 0.0,
+            "derivative equals": np.nan,
+            "derivative down efficiency": 1.0,
+            "derivative up efficiency": 1.0,
+        },
+        index=index,
+    )
+
+    ems_constraints = pd.DataFrame(
+        {"derivative min": -200.0, "derivative max": 200.0},
+        index=index,
+    )
+
+    # Coupling group: one reference device (gas, coeff 1.0) and two coupled
+    # devices (heat with coeff -0.5, power with coeff -0.3).
+    coupling_groups = {"chp": [(0, 1.0), (1, -0.5), (2, -0.3)]}
+
+    # Gas-price commitment gives the objective a finite value and models the
+    # cost of consuming gas. With quantity=0 and both prices set the
+    # commitment acts as a two-sided soft equality: any upward deviation
+    # (gas consumption) incurs a cost of 1 EUR/kW.
+    gas_price_commitment = FlowCommitment(
+        name="gas cost",
+        index=index,
+        quantity=pd.Series(0.0, index=index),
+        upwards_deviation_price=pd.Series(1.0, index=index),
+        downwards_deviation_price=pd.Series(0.0, index=index),
+        device=pd.Series(0, index=index),
+    )
+
+    schedules, planned_costs, results, model = device_scheduler(
+        device_constraints=[gas_constraints, heat_constraints, power_constraints],
+        ems_constraints=ems_constraints,
+        commitments=[gas_price_commitment],
+        coupling_groups=coupling_groups,
+    )
+
+    assert (
+        results.solver.termination_condition == "optimal"
+    ), "Solver did not find an optimal solution."
+
+    # Heat is fixed to -10 kW by derivative_equals.
+    pd.testing.assert_series_equal(
+        schedules[1],
+        pd.Series(-10.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="heat output forced to -10 kW by derivative_equals",
+    )
+
+    # Coupling: P_gas / 1.0 == P_heat / -0.5  →  P_gas = -10 / -0.5 = 20 kW
+    pd.testing.assert_series_equal(
+        schedules[0],
+        pd.Series(20.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="gas consumption determined by coupling (20 kW from 10 kW heat at coeff -0.5)",
+    )
+
+    # Coupling: P_gas / 1.0 == P_power / -0.3  →  P_power = 20 / -0.3 = -6 kW
+    pd.testing.assert_series_equal(
+        schedules[2],
+        pd.Series(-6.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="power output determined by coupling (-0.3 * alpha = -0.3 * 20 = -6 kW)",
+    )
+
+
+def test_dual_fuel_chp_coupling():
+    """Test coupling_groups with two input devices (dual-fuel CHP).
+
+    Models a CHP unit that consumes equal parts natural gas and hydrogen,
+    producing heat and electricity:
+
+    - d=0  gas input:      can only consume gas           (derivative_min=0)
+    - d=1  hydrogen input: can only consume hydrogen      (derivative_min=0)
+    - d=2  heat output:    can only produce heat          (derivative_max=0)
+    - d=3  power output:   can only produce electricity   (derivative_max=0)
+
+    Coupling group ``"chp"`` with coefficients
+    ``[(0, 0.5), (1, 0.5), (2, -0.5), (3, -0.3)]`` introduces a free variable
+    ``alpha`` and enforces ``P[d] == coeff * alpha``:
+
+        P_gas      =  0.5 * alpha   (50% of total fuel from gas)
+        P_hydrogen =  0.5 * alpha   (50% of total fuel from hydrogen)
+        P_heat     = -0.5 * alpha   (heat efficiency 50% of total fuel)
+        P_power    = -0.3 * alpha   (power efficiency 30% of total fuel)
+
+    Because gas and hydrogen share the same coefficient the two fuel flows are
+    always equal, confirming that device order does not affect the result.
+
+    Heat production is forced to exactly 10 kW via ``derivative equals = -10`` on device 2.
+    Substituting ``P_heat = -10`` gives ``alpha = 20``, so:
+
+        P_gas      =  10 kW   (equal gas input)
+        P_hydrogen =  10 kW   (equal hydrogen input)
+        P_heat     = -10 kW   (heat produced, forced)
+        P_power    =  -6 kW   (electricity produced)
+    """
+    start = pd.Timestamp("2026-01-01T00:00+01:00")
+    end = pd.Timestamp("2026-01-01T04:00+01:00")
+    resolution = pd.Timedelta("1h")
+    index = initialize_index(start=start, end=end, resolution=resolution)
+
+    def _flow_df(**kwargs) -> pd.DataFrame:
+        defaults = {
+            "min": np.nan,
+            "max": np.nan,
+            "equals": np.nan,
+            "derivative min": 0.0,
+            "derivative max": 0.0,
+            "derivative equals": np.nan,
+            "derivative down efficiency": 1.0,
+            "derivative up efficiency": 1.0,
+        }
+        defaults.update(kwargs)
+        return pd.DataFrame(defaults, index=index)
+
+    # d=0: gas input — can only consume, capacity 100 kW
+    gas_constraints = _flow_df(**{"derivative max": 100.0})
+    # d=1: hydrogen input — can only consume, capacity 100 kW
+    hydrogen_constraints = _flow_df(**{"derivative max": 100.0})
+    # d=2: heat output — can only produce, forced to -10 kW
+    heat_constraints = _flow_df(
+        **{"derivative min": -100.0, "derivative equals": -10.0}
+    )
+    # d=3: power output — can only produce, free (coupling determines value)
+    power_constraints = _flow_df(**{"derivative min": -100.0})
+
+    ems_constraints = pd.DataFrame(
+        {"derivative min": -200.0, "derivative max": 200.0},
+        index=index,
+    )
+
+    # Both fuel inputs share coefficient 0.5, so they receive identical flows.
+    # Outputs have negative coefficients equal to their efficiency fractions.
+    coupling_groups = {"chp": [(0, 0.5), (1, 0.5), (2, -0.5), (3, -0.3)]}
+
+    # Gas-price commitment for device 0 just to give the objective a finite value
+    # Even though hydrogen is free, it will still be used because its consumption is coupled to gas.
+    fuel_cost_commitment = FlowCommitment(
+        name="fuel cost",
+        index=index,
+        quantity=pd.Series(0.0, index=index),
+        upwards_deviation_price=pd.Series(1.0, index=index),
+        downwards_deviation_price=pd.Series(0.0, index=index),
+        device=pd.Series(0, index=index),
+    )
+
+    schedules, _costs, results, _model = device_scheduler(
+        device_constraints=[
+            gas_constraints,
+            hydrogen_constraints,
+            heat_constraints,
+            power_constraints,
+        ],
+        ems_constraints=ems_constraints,
+        commitments=[fuel_cost_commitment],
+        coupling_groups=coupling_groups,
+    )
+
+    assert (
+        results.solver.termination_condition == "optimal"
+    ), "Solver did not find an optimal solution."
+
+    # Heat is fixed to -10 kW; alpha = -10 / -0.5 = 20.
+    pd.testing.assert_series_equal(
+        schedules[2],
+        pd.Series(-10.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="heat output forced to -10 kW by derivative_equals",
+    )
+
+    # Coupling: P_gas = 0.5 * alpha = 0.5 * 20 = 10 kW
+    pd.testing.assert_series_equal(
+        schedules[0],
+        pd.Series(10.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="gas input = 0.5 * alpha = 10 kW",
+    )
+
+    # Coupling: P_hydrogen = 0.5 * alpha = 10 kW (equal to gas)
+    pd.testing.assert_series_equal(
+        schedules[1],
+        pd.Series(10.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="hydrogen input = 0.5 * alpha = 10 kW (equal to gas input)",
+    )
+
+    # Coupling: P_power = -0.3 * alpha = -0.3 * 20 = -6 kW
+    pd.testing.assert_series_equal(
+        schedules[3],
+        pd.Series(-6.0, index=index),
+        check_names=False,
+        rtol=1e-4,
+        obj="power output = -0.3 * alpha = -6 kW",
+    )
+
+
+def _run_factory_scenario(
+    gas_price: float,
+    elec_price: float,
+) -> tuple:
+    """Run the simplified factory scenario and return the 7 device schedules.
+
+    Devices
+    ~~~~~~~
+    d=0  e-heater       electricity   → heat coupling (ems_power ≥ 0, i.e. consumes electricity)
+    d=1  gas boiler     gas           → heat coupling (ems_power ≥ 0, i.e. consumes gas)
+    d=2  steamer        heat coupling → steam         (ems_power ≤ 0, i.e. produces steam)
+    d=3  CHP gas input  gas           → chp coupling  (ems_power ≥ 0, i.e. consumes gas, coupling member = alpha)
+    d=4  CHP heat out   chp coupling  → steam         (ems_power ≤ 0, i.e. produces steam, coupling member = -0.5 alpha)
+    d=5  CHP power out  chp coupling  → electricity   (ems_power ≤ 0, i.e. produces electricity, coupling member = -0.3 alpha)
+    d=6  steam demand   steam         → fixed flow    (ems_power = 15, i.e. consumes steam)
+
+    CHP coupling coefficients
+    ~~~~~~~~~~~~~~~~~~~~~~~~~
+    The coupling constraint introduces a free variable ``alpha`` (the normalised gas flow)
+    and enforces ``P[d_i] == coeff_i * alpha`` for every device in the group.
+    Choosing thermal efficiency η_heat = 0.5 and power efficiency η_power = 0.3,
+    the coefficients simply become the signed efficiency fractions::
+
+        P_gas   =  1.0   * alpha   (input,  coeff = 1.0)
+        P_heat  = -0.5   * alpha   (output, coeff = η_heat  = -0.5)
+        P_power = -0.3   * alpha   (output, coeff = η_power = −0.3)
+
+    """
+    ETA_HEAT = 0.5  # fraction of CHP gas input that becomes heat
+    ETA_POWER = 0.3  # fraction of CHP gas input that becomes power
+    STEAM_DEMAND = 15.0  # kW, constant heat drain representing steam production
+    CHP_GAS_MAX = 20.0  # kW, maximum gas input to CHP
+    BOILER_GAS_MAX = 10.0  # kW, maximum gas input to gas boiler
+    HEATER_POWER_MAX = 100.0  # kW, maximum electricity input to e-heater
+
+    start = pd.Timestamp("2026-01-01T00:00+01:00")
+    end = pd.Timestamp("2026-01-01T04:00+01:00")
+    resolution = pd.Timedelta("1h")
+    index = initialize_index(start=start, end=end, resolution=resolution)
+
+    def _df(**kwargs) -> pd.DataFrame:
+        """Build a device-constraints DataFrame with defaults for unused columns."""
+        defaults = {
+            "min": np.nan,
+            "max": np.nan,
+            "equals": np.nan,
+            "derivative min": 0.0,
+            "derivative max": 0.0,
+            "derivative equals": np.nan,
+            "derivative down efficiency": 1.0,
+            "derivative up efficiency": 1.0,
+            "stock delta": 0.0,
+        }
+        defaults.update(kwargs)
+        return pd.DataFrame(defaults, index=index)
+
+    device_constraints = [
+        # d=0  e-heater: heat-node reference device. The min=max=0 forces the heat
+        #       node to balance at every step (zero-capacity flow node), making
+        #       the per-step dispatch deterministic despite flat prices.
+        _df(min=0.0, max=0.0, **{"derivative max": HEATER_POWER_MAX}),
+        # d=1  gas boiler: up to 100 kW gas → 100 kW heat (efficiency 1 for clean maths in test)
+        _df(**{"derivative max": BOILER_GAS_MAX, "commodity": "gas"}),
+        # d=2  steamer: can only produce steam (negative ems_power).
+        # The lower bound is finite to avoid unbounded model messages while still
+        # being looser than the upstream heat-supply limits.
+        _df(
+            **{
+                "derivative min": -(HEATER_POWER_MAX + BOILER_GAS_MAX),
+                "derivative max": 0.0,
+                "commodity": "steam",
+            }
+        ),
+        # d=3  CHP gas input: up to CHP_GAS_MAX kW gas
+        _df(**{"derivative max": CHP_GAS_MAX, "commodity": "gas"}),
+        # d=4  CHP heat output: positive ems_power adds heat to the steam node.
+        #      The min=max=0 forces the steam node to balance at every step.
+        _df(
+            min=0.0,
+            max=0.0,
+            **{
+                "derivative min": -CHP_GAS_MAX * ETA_HEAT,
+                "derivative max": 0.0,
+                "commodity": "steam",
+            },
+        ),
+        # d=5  CHP power output: negative ems_power only (production)
+        _df(**{"derivative min": -CHP_GAS_MAX * ETA_POWER, "derivative max": 0.0}),
+        # d=6  steam demand: fixed steam consumption at STEAM_DEMAND kW.
+        _df(
+            **{
+                "derivative min": STEAM_DEMAND,
+                "derivative max": STEAM_DEMAND,
+                "commodity": "steam",
+            }
+        ),
+    ]
+
+    ems_constraints = pd.DataFrame(
+        {"derivative min": -300.0, "derivative max": 300.0},
+        index=index,
+    )
+
+    # stock group: all heat-buffer devices share the same stock
+    # (key 0 is an arbitrary group id, not a device index)
+    heat_group_id = 0
+    steam_group_id = 1
+    stock_groups = {heat_group_id: [0, 1, 2], steam_group_id: [2, 4, 6]}
+
+    # CHP coupling: coefficients are signed efficiency fractions.
+    # coeff_heat  = -η_heat  = -0.5 →  P_heat = -0.5 * alpha = -0.5 * P_gas
+    # coeff_power = -η_power = -0.3 → P_power = -0.3 * alpha = -0.3 * P_gas
+    coupling_groups = {
+        "chp": [
+            (3, 1.0),
+            (4, -ETA_HEAT),  # = -0.5
+            (5, -ETA_POWER),  # = -0.3
+        ]
+    }
+
+    # --- energy-price commitments -------------------------------------------
+    # Gas price applies to gas boiler (d=1) and CHP gas input (d=3).
+    # Electricity price applies to e-heater (d=0) and CHP power output (d=5).
+    # Using both upwards and downwards prices makes each commitment a two-sided
+    # soft equality (quantity = 0):
+    #   • upward deviation  = consuming more than 0  → positive cost
+    #   • downward deviation = producing (negative flow) → negative cost (revenue)
+    gas_p = pd.Series(gas_price, index=index)
+    elec_p = pd.Series(elec_price, index=index)
+
+    commitments = []
+    for d, price in [(1, gas_p), (3, gas_p), (0, elec_p), (5, elec_p)]:
+        commitments.append(
+            FlowCommitment(
+                name="gas cost" if d in (1, 3) else "electricity cost",
+                index=index,
+                quantity=pd.Series(0.0, index=index),
+                upwards_deviation_price=price,
+                downwards_deviation_price=price,
+                device=pd.Series(d, index=index),
+            )
+        )
+
+    schedules, _costs, results, _model = device_scheduler(
+        device_constraints=device_constraints,
+        ems_constraints=ems_constraints,
+        commitments=commitments,
+        stock_groups=stock_groups,
+        coupling_groups=coupling_groups,
+    )
+
+    assert results.solver.termination_condition == "optimal", (
+        f"Solver did not find an optimal solution "
+        f"(gas_price={gas_price}, elec_price={elec_price})"
+    )
+    return tuple(schedules)
+
+
+def test_factory_chp_dispatch():
+    """Factory: CHP + gas boiler + e-heater competing to meet a fixed steam demand.
+
+    The shared heat buffer (modelled via ``stock_groups``) is drained at a
+    constant rate of 15 kW by the steam demand device. Two price scenarios
+    verify that the optimizer correctly chooses the cheapest heat source.
+
+    Scenario A — gas cheaper than electricity
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Prices: gas = 20 EUR/kW, electricity = 50 EUR/kW.
+
+    Effective cost per kW of heat delivered:
+    - CHP:       gas_cost − power_revenue  = (20·20 − 50·6) / 10 = 10 EUR/kW
+    - gas boiler: 20 EUR/kW  (efficiency = 1)
+    - e-heater:   50 EUR/kW  (efficiency = 1)
+
+    Merit order: CHP ≪ gas boiler ≪ e-heater.
+
+    With CHP at maximum (20 kW gas → 10 kW heat + 6 kW power):
+    - remaining heat demand = 15 − 10 = 5 kW → gas boiler
+    - e-heater not needed
+
+    Scenario B — electricity cheaper than gas
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Prices: gas = 100 EUR/kW, electricity = 10 EUR/kW.
+
+    Effective cost per kW of heat:
+    - CHP:       (100·20 − 10·6) / 10 = 194 EUR/kW
+    - gas boiler: 100 EUR/kW
+    - e-heater:   10 EUR/kW
+
+    Merit order: e-heater ≪ gas boiler ≪ CHP.
+
+    All 15 kW steam demand is met by the e-heater; CHP and gas boiler are off.
+
+    Scenario C — gas slightly cheaper
+    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    Prices: gas = 50 EUR/kW, electricity = 55 EUR/kW.
+
+    Effective cost per kW of heat delivered:
+    - CHP:       gas_cost − power_revenue  = (50·20 − 55·6) / 10 = 67 EUR/kW
+    - gas boiler: 50 EUR/kW
+    - e-heater:   55 EUR/kW
+
+    Merit order: gas boiler ≪ e-heater ≪ CHP.
+
+    With gas boiler at maximum (10 kW gas → 10 kW heat):
+    - remaining heat demand = 15 − 10 = 5 kW → e-heater
+    - CHP not needed
+    """
+    # ------------------------------------------------------------------ #
+    # Scenario A: gas cheaper — CHP at max, gas boiler fills the rest    #
+    # ------------------------------------------------------------------ #
+    (e_heater, gas_boiler, steamer, chp_gas, chp_heat, chp_power, demand) = (
+        _run_factory_scenario(gas_price=20.0, elec_price=50.0)
+    )
+
+    expected_chp_gas = pd.Series(20.0, index=e_heater.index)
+    expected_chp_heat = pd.Series(-10.0, index=e_heater.index)  # -0.5 * 20
+    expected_chp_power = pd.Series(-6.0, index=e_heater.index)  # -0.3 * 20
+    expected_boiler = pd.Series(5.0, index=e_heater.index)  # fills 15-10 kW gap
+    expected_steamer = pd.Series(-5.0, index=e_heater.index)
+    expected_demand = pd.Series(15.0, index=e_heater.index)
+    expected_eheater = pd.Series(0.0, index=e_heater.index)
+
+    pd.testing.assert_series_equal(
+        chp_gas,
+        expected_chp_gas,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario A: CHP gas input at maximum (20 kW)",
+    )
+    pd.testing.assert_series_equal(
+        chp_heat,
+        expected_chp_heat,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario A: CHP heat output = 0.5 × gas input (10 kW)",
+    )
+    pd.testing.assert_series_equal(
+        chp_power,
+        expected_chp_power,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario A: CHP power output = −0.3 × gas input (−6 kW)",
+    )
+    pd.testing.assert_series_equal(
+        gas_boiler,
+        expected_boiler,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario A: gas boiler fills remaining 5 kW heat demand",
+    )
+    pd.testing.assert_series_equal(
+        steamer,
+        expected_steamer,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario A: steamer supplies remaining 5 kW steam",
+    )
+    pd.testing.assert_series_equal(
+        demand,
+        expected_demand,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario A: steam demand fixed at 15 kW",
+    )
+    pd.testing.assert_series_equal(
+        e_heater,
+        expected_eheater,
+        check_names=False,
+        atol=1e-4,
+        obj="Scenario A: e-heater not used (gas is cheapest)",
+    )
+
+    # ------------------------------------------------------------------ #
+    # Scenario B: electricity cheaper — e-heater meets all demand        #
+    # ------------------------------------------------------------------ #
+    (e_heater, gas_boiler, steamer, chp_gas, chp_heat, chp_power, demand) = (
+        _run_factory_scenario(gas_price=100.0, elec_price=10.0)
+    )
+
+    expected_eheater_b = pd.Series(15.0, index=e_heater.index)
+    expected_zero = pd.Series(0.0, index=e_heater.index)
+    expected_steamer_b = pd.Series(-15.0, index=e_heater.index)
+    expected_demand_b = pd.Series(15.0, index=e_heater.index)
+
+    pd.testing.assert_series_equal(
+        e_heater,
+        expected_eheater_b,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario B: e-heater meets all 15 kW steam demand",
+    )
+    pd.testing.assert_series_equal(
+        chp_gas,
+        expected_zero,
+        check_names=False,
+        atol=1e-4,
+        obj="Scenario B: CHP not used (electricity is cheapest)",
+    )
+    pd.testing.assert_series_equal(
+        gas_boiler,
+        expected_zero,
+        check_names=False,
+        atol=1e-4,
+        obj="Scenario B: gas boiler not used (electricity is cheapest)",
+    )
+    pd.testing.assert_series_equal(
+        steamer,
+        expected_steamer_b,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario B: steamer supplies all 15 kW steam",
+    )
+    pd.testing.assert_series_equal(
+        demand,
+        expected_demand_b,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario B: steam demand fixed at 15 kW",
+    )
+
+    # --------------------------------------------------------------------------------- #
+    # Scenario C: gas slightly cheaper — gas boiler at max, e-heater fills the rest     #
+    # --------------------------------------------------------------------------------- #
+    (e_heater, gas_boiler, steamer, chp_gas, chp_heat, chp_power, demand) = (
+        _run_factory_scenario(gas_price=50.0, elec_price=55.0)
+    )
+
+    expected_chp_gas = pd.Series(0.0, index=e_heater.index)
+    expected_chp_heat = pd.Series(0.0, index=e_heater.index)
+    expected_chp_power = pd.Series(0.0, index=e_heater.index)
+    expected_boiler = pd.Series(10.0, index=e_heater.index)
+    expected_steamer = pd.Series(-15.0, index=e_heater.index)
+    expected_demand = pd.Series(15.0, index=e_heater.index)
+    expected_eheater = pd.Series(5.0, index=e_heater.index)  # fills 15-10 kW gap
+
+    pd.testing.assert_series_equal(
+        chp_gas,
+        expected_chp_gas,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario C: CHP not used",
+    )
+    pd.testing.assert_series_equal(
+        chp_heat,
+        expected_chp_heat,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario C: CHP not used",
+    )
+    pd.testing.assert_series_equal(
+        chp_power,
+        expected_chp_power,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario C: CHP not used",
+    )
+    pd.testing.assert_series_equal(
+        gas_boiler,
+        expected_boiler,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario C: gas boiler at maximum (10 kW)",
+    )
+    pd.testing.assert_series_equal(
+        steamer,
+        expected_steamer,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario C: steamer supplies all 15 kW steam",
+    )
+    pd.testing.assert_series_equal(
+        demand,
+        expected_demand,
+        check_names=False,
+        rtol=1e-4,
+        obj="Scenario C: steam demand fixed at 15 kW",
+    )
+    pd.testing.assert_series_equal(
+        e_heater,
+        expected_eheater,
+        check_names=False,
+        atol=1e-4,
+        obj="Scenario C: e-heater fills remaining 5 kW heat demand",
+    )
+
+
 def test_all_gas_flex_model_without_electricity_device(app, db):
     """test_all_gas_flex_model_without_electricity_device: a flex-model with only gas
     devices (no electricity device at all) should not raise a KeyError, now that
diff --git a/flexmeasures/data/models/planning/tests/test_storage.py b/flexmeasures/data/models/planning/tests/test_storage.py
index a03d09e05d..30314a3e24 100644
--- a/flexmeasures/data/models/planning/tests/test_storage.py
+++ b/flexmeasures/data/models/planning/tests/test_storage.py
@@ -7,6 +7,7 @@
 import numpy as np
 import pandas as pd
 
+from flexmeasures.data.models.generic_assets import GenericAsset, GenericAssetType
 from flexmeasures.data.models.planning import Scheduler
 from flexmeasures.data.models.planning.storage import StorageScheduler
 from flexmeasures.data.models.planning.utils import initialize_index
@@ -16,6 +17,7 @@
     get_sensors_from_db,
     series_to_ts_specs,
 )
+from flexmeasures.data.services.utils import get_or_create_model
 from flexmeasures.data.services.scheduling_result import SchedulingJobResult
 
 
@@ -1223,3 +1225,187 @@ def test_resolve_soc_at_start_from_percent_sensor_uses_device_sensor_fallback(
         )
         == 2.5
     )
+
+
+def test_storage_scheduler_chp_coupling(app, db):
+    """Test that the StorageScheduler enforces CHP coupling constraints between devices.
+
+    Models a Combined Heat and Power unit with three sensors.
+
+    In the flex-model, the coupling coefficients are entered as positive magnitudes::
+
+        gas input   -> 1.0
+        heat output -> 0.5
+        power output -> 0.3
+
+    Internally, the CHP is interpreted with the signed commodity-flow coefficients::
+
+        P_gas   ->  1.0
+        P_heat  -> -0.5
+        P_power -> -0.3
+
+    The returned storage schedule for the heat buffer is still positive, because this
+    test uses the storage sign convention for buffer charging.
+
+    - d=0  gas input:    CHP gas consumption
+    - d=1  heat output:  CHP heat -> heat buffer
+    - d=2  power output: CHP electricity production
+
+    The heat output is forced to exactly 5 kW per step by combining:
+    - ``production-capacity: "0 kW"``  (hard lower bound: derivative_min = 0)
+    - ``consumption-capacity: "5 kW"`` (hard upper bound: derivative_max = 0.005 MW)
+    - ``soc-targets`` requiring 20 kWh at the end of the 4-hour window
+
+    With soc_at_start = 0 and max 5 kW over 4 × 1-hour steps the only feasible
+    solution is P_heat = 5 kW every step. Substituting P_heat = 5 kW gives
+    alpha = 5 / 0.5 = 10 kW, so:
+
+        P_gas   =  1.0 × 10 kW = 10 kW
+        P_power = −0.3 × 10 kW = −3 kW
+    """
+    # ---- asset type + asset
+    chp_type = get_or_create_model(GenericAssetType, name="chp-plant")
+    chp = GenericAsset(name="CHP plant (coupling test)", generic_asset_type=chp_type)
+    db.session.add(chp)
+    db.session.flush()
+
+    # ---- schedule window
+    start = pd.Timestamp("2026-01-01T00:00:00+01:00")
+    end = pd.Timestamp("2026-01-01T04:00:00+01:00")
+    resolution = timedelta(hours=1)
+
+    # CHP efficiencies (same values as the factory scenario in test_commitments.py)
+    ETA_HEAT = 0.5  # fraction of gas input that becomes heat
+    ETA_POWER = 0.3  # fraction of gas input that becomes electricity
+
+    # ---- sensors
+    gas_input_sensor = Sensor(
+        name="CHP gas input (coupling test)",
+        generic_asset=chp,
+        unit="MW",
+        event_resolution=resolution,
+    )
+    heat_output_sensor = Sensor(
+        name="CHP heat output (coupling test)",
+        generic_asset=chp,
+        unit="MW",
+        event_resolution=resolution,
+    )
+    power_output_sensor = Sensor(
+        name="CHP power output (coupling test)",
+        generic_asset=chp,
+        unit="MW",
+        event_resolution=resolution,
+    )
+    db.session.add_all([gas_input_sensor, heat_output_sensor, power_output_sensor])
+    db.session.flush()
+
+    # ---- flex model
+    # Flex-model coupling-coefficients are user-facing positive magnitudes.
+    # The intended internal CHP coefficients are +1.0 for gas, -0.5 for heat,
+    # and -0.3 for power.
+    flex_model = [
+        {
+            # d=0: gas input — pure flow device (no SoC), can only consume gas.
+            "sensor": gas_input_sensor.id,
+            "power-capacity": "20 kW",
+            "production-capacity": "0 kW",  # derivative_min = 0
+            "coupling": "chp",
+            "coupling-coefficient": 1.0,
+        },
+        {
+            # d=1: heat output — tracks heat-buffer SoC, positive ems_power = heat
+            # added to buffer. The SoC target forces P_heat = 5 kW per step.
+            "sensor": heat_output_sensor.id,
+            "soc-at-start": "0 MWh",
+            "soc-min": "0 MWh",
+            "soc-max": "0.02 MWh",  # 20 kWh — matches the SoC target
+            "soc-targets": [
+                {
+                    # Single target at the schedule end: cumulative heat = 20 kWh.
+                    # With max 5 kW and 4 × 1 h steps the only feasible solution
+                    # is 5 kW every step.
+                    "start": "2026-01-01T04:00:00+01:00",
+                    "duration": "PT1H",
+                    "value": "0.02 MWh",
+                }
+            ],
+            "power-capacity": "5 kW",
+            "consumption-capacity": "5 kW",
+            "production-capacity": "0 kW",  # can only add heat, not extract
+            "prefer-charging-sooner": True,
+            "coupling": "chp",
+            "coupling-coefficient": ETA_HEAT,  # = 0.5
+        },
+        {
+            # d=2: power output — pure flow device (no SoC), can only produce
+            # electricity (negative ems_power).
+            "sensor": power_output_sensor.id,
+            "power-capacity": "6 kW",
+            "consumption-capacity": "0 kW",  # derivative_max = 0
+            "coupling": "chp",
+            "coupling-coefficient": ETA_POWER,  # = 0.3 (sign inferred from capacities)
+        },
+    ]
+
+    flex_context = {
+        "consumption-price": "50 EUR/MWh",
+        "production-price": "50 EUR/MWh",
+        "site-power-capacity": "1 MW",  # large enough to avoid EMS constraints
+    }
+
+    scheduler = StorageScheduler(
+        asset_or_sensor=chp,
+        start=start,
+        end=end,
+        resolution=resolution,
+        flex_model=flex_model,
+        flex_context=flex_context,
+        return_multiple=True,
+    )
+
+    results = scheduler.compute(skip_validation=True)
+
+    # ---- extract storage schedules per sensor
+    storage_schedules = {
+        r["sensor"]: r["data"] for r in results if r.get("name") == "storage_schedule"
+    }
+
+    assert gas_input_sensor in storage_schedules, "Gas input schedule missing"
+    assert heat_output_sensor in storage_schedules, "Heat output schedule missing"
+    assert power_output_sensor in storage_schedules, "Power output schedule missing"
+
+    gas_schedule = storage_schedules[gas_input_sensor]
+    heat_schedule = storage_schedules[heat_output_sensor]
+    power_schedule = storage_schedules[power_output_sensor]
+
+    # The SoC target of 20 kWh is met after 4 × 1-hour steps at 5 kW.
+    # The schedule index runs from ``start`` to ``end`` inclusive (5 time slots),
+    # so the last slot has no binding SoC constraint and the CHP is idle there.
+    # All assertions therefore apply to the first four active slots only.
+    active_steps = slice(None, -1)  # exclude the final trailing idle slot
+
+    # Heat output is forced to exactly 5 kW per step by the SoC target.
+    # alpha = P_heat / ETA_HEAT = 0.005 / 0.5 = 0.010 MW
+    np.testing.assert_allclose(
+        heat_schedule.iloc[active_steps],
+        0.005,  # 5 kW expressed in MW
+        rtol=1e-4,
+        err_msg="Heat output should be exactly 5 kW per step (forced by SoC target)",
+    )
+
+    # Coupling: P_gas = 1.0 * alpha = 0.010 MW = 10 kW
+    np.testing.assert_allclose(
+        gas_schedule.iloc[active_steps],
+        0.010,  # 10 kW expressed in MW
+        rtol=1e-4,
+        err_msg="Gas input must be 10 kW — determined by coupling (1.0 * alpha)",
+    )
+
+    # Coupling: P_power = -ETA_POWER * alpha = -0.3 * 0.010 MW = -0.003 MW = -3 kW
+    np.testing.assert_allclose(
+        power_schedule.iloc[active_steps],
+        -0.003,  # -3 kW expressed in MW
+        rtol=1e-4,
+        err_msg="Power output must be -3 kW — determined by coupling (-0.3 * alpha)",
+    )
diff --git a/flexmeasures/data/schemas/scheduling/metadata.py b/flexmeasures/data/schemas/scheduling/metadata.py
index 4f6f9d4298..e6b531b36e 100644
--- a/flexmeasures/data/schemas/scheduling/metadata.py
+++ b/flexmeasures/data/schemas/scheduling/metadata.py
@@ -223,6 +223,25 @@ def to_dict(self):
 """,
     examples=["electricity", "gas"],
 )
+COUPLING = MetaData(
+    description="""Name of the coupling group this device belongs to.
+Devices sharing the same coupling name are constrained to have proportionally related power flows, via a hard equality constraint.
+Use this to model a device that converts one commodity into another, by describing each of its commodity ports as a separate device.
+For example, a combined heat and power (CHP) unit is described as a gas input device, a heat output device and an electricity output device, all sharing one coupling name.
+Use together with ``coupling-coefficient`` to set the flow ratios.
+""",
+    example="chp",
+)
+COUPLING_COEFFICIENT = MetaData(
+    description="""Positive coupling magnitude for this device within its coupling group.
+The scheduler couples the power flows of all devices in the group: each device's power is its coupling coefficient times the group's common flow level.
+The flow direction of each device is inferred from its directional capacities: a device with ``consumption-capacity: "0 kW"`` is an output (producing) device, and a device with ``production-capacity: "0 kW"`` is an input (consuming) device.
+Exactly one of the two directional capacities must be set to a fixed zero for each coupled device.
+For example, a CHP unit with 50% thermal and 30% electrical efficiency uses a gas input device (coefficient 1), a heat output device (coefficient 0.5) and an electricity output device (coefficient 0.3).
+Defaults to 1.
+""",
+    example=0.5,
+)
 CONSUMPTION = MetaData(
     description="""Sensor used to record the scheduled power as seen from a consumption perspective.
 
diff --git a/flexmeasures/data/schemas/scheduling/storage.py b/flexmeasures/data/schemas/scheduling/storage.py
index 60d9da3448..06325375c8 100644
--- a/flexmeasures/data/schemas/scheduling/storage.py
+++ b/flexmeasures/data/schemas/scheduling/storage.py
@@ -245,6 +245,19 @@ class StorageFlexModelSchema(Schema):
         validate=validate.Length(min=1),
         metadata=metadata.SOC_USAGE.to_dict(),
     )
+    coupling = fields.Str(
+        data_key="coupling",
+        required=False,
+        load_default=None,
+        metadata=metadata.COUPLING.to_dict(),
+    )
+    coupling_coefficient = fields.Float(
+        data_key="coupling-coefficient",
+        required=False,
+        load_default=1.0,
+        validate=validate.Range(min=0, min_inclusive=False),
+        metadata=metadata.COUPLING_COEFFICIENT.to_dict(),
+    )
 
     def __init__(
         self,
@@ -394,6 +407,33 @@ def validate_commodity(self, commodity: str, **kwargs):
         if not isinstance(commodity, str) or not commodity.strip():
             raise ValidationError("commodity must be a non-empty string.")
 
+    @validates_schema
+    def validate_coupling_direction_is_unambiguous(self, data: dict, **kwargs):
+        """A coupled device must have an inferable flow direction.
+
+        The sign of the coupling coefficient is inferred from directional capacities:
+        a fixed zero consumption-capacity marks an output device, a fixed zero
+        production-capacity marks an input device. When both directional capacities
+        allow flow (or are given in a form whose value cannot be checked statically,
+        such as a sensor reference), the direction is ambiguous, so we reject the
+        flex-model rather than silently treating the device as an input.
+        """
+        if data.get("coupling") is None:
+            return
+
+        def _is_fixed_zero(value) -> bool:
+            return isinstance(value, ur.Quantity) and float(value.magnitude) == 0.0
+
+        consumption_blocked = _is_fixed_zero(data.get("consumption_capacity"))
+        production_blocked = _is_fixed_zero(data.get("production_capacity"))
+        if consumption_blocked == production_blocked:
+            raise ValidationError(
+                "A device with a 'coupling' field must have an unambiguous flow direction: "
+                "set either its consumption-capacity (for an output device) or its "
+                "production-capacity (for an input device) to a fixed 0, but not both.",
+                field_name="coupling",
+            )
+
     @post_load
     def post_load_sequence(self, data: dict, **kwargs) -> dict:
         """Perform some checks and corrections after we loaded."""
diff --git a/flexmeasures/data/schemas/tests/test_scheduling.py b/flexmeasures/data/schemas/tests/test_scheduling.py
index 90842d2cf9..7d557ee178 100644
--- a/flexmeasures/data/schemas/tests/test_scheduling.py
+++ b/flexmeasures/data/schemas/tests/test_scheduling.py
@@ -1197,3 +1197,46 @@ def test_asset_trigger_schema_rejects_malformed_flex_context(app):
     with pytest.raises(ValidationError) as e_info:
         schema.normalize_flex_context_format({"flex-context": "not-a-dict-or-list"})
     assert "flex-context" in str(e_info.value)
+
+
+@pytest.mark.parametrize(
+    "capacity_fields, fails",
+    [
+        # Input device: production blocked, direction is unambiguous
+        ({"production-capacity": "0 kW"}, False),
+        # Output device: consumption blocked, direction is unambiguous
+        ({"consumption-capacity": "0 kW"}, False),
+        # Output device with a bounded input side still has one blocked direction
+        ({"consumption-capacity": "5 kW", "production-capacity": "0 kW"}, False),
+        # Neither direction blocked: ambiguous
+        ({}, True),
+        # Both directions open: ambiguous
+        ({"consumption-capacity": "5 kW", "production-capacity": "5 kW"}, True),
+        # Both directions blocked: degenerate (device pinned to zero flow)
+        ({"consumption-capacity": "0 kW", "production-capacity": "0 kW"}, True),
+    ],
+)
+def test_coupling_direction_must_be_unambiguous(app, capacity_fields, fails):
+    """test_coupling_direction_must_be_unambiguous: a device with a `coupling` field must
+    have exactly one directional capacity fixed to zero, so the sign of its coupling
+    coefficient can be inferred."""
+    schema = StorageFlexModelSchema(start=datetime(2026, 6, 1), sensor=None)
+    flex_model = {
+        "power-capacity": "20 kW",
+        "coupling": "chp",
+        "coupling-coefficient": 0.5,
+        **capacity_fields,
+    }
+    if fails:
+        with pytest.raises(ValidationError) as e_info:
+            schema.load(flex_model)
+        assert "unambiguous flow direction" in str(e_info.value)
+    else:
+        schema.load(flex_model)
+
+
+def test_uncoupled_device_needs_no_directional_capacities(app):
+    """test_uncoupled_device_needs_no_directional_capacities: the coupling-direction check
+    only applies to devices that define a `coupling` field."""
+    schema = StorageFlexModelSchema(start=datetime(2026, 6, 1), sensor=None)
+    schema.load({"power-capacity": "20 kW"})
diff --git a/flexmeasures/ui/static/openapi-specs.json b/flexmeasures/ui/static/openapi-specs.json
index 0eab864eef..ec37ffbf09 100644
--- a/flexmeasures/ui/static/openapi-specs.json
+++ b/flexmeasures/ui/static/openapi-specs.json
@@ -6369,6 +6369,22 @@
             ],
             "items": {}
           },
+          "coupling": {
+            "type": [
+              "string",
+              "null"
+            ],
+            "default": null,
+            "description": "Name of the coupling group this device belongs to.\nDevices sharing the same coupling name are constrained to have proportionally related power flows, via a hard equality constraint.\nUse this to model a device that converts one commodity into another, by describing each of its commodity ports as a separate device.\nFor example, a combined heat and power (CHP) unit is described as a gas input device, a heat output device and an electricity output device, all sharing one coupling name.\nUse together with <code>coupling-coefficient</code> to set the flow ratios.\n",
+            "example": "chp"
+          },
+          "coupling-coefficient": {
+            "type": "number",
+            "default": 1.0,
+            "minimum": 0.0,
+            "description": "Positive coupling magnitude for this device within its coupling group.\nThe scheduler couples the power flows of all devices in the group: each device's power is its coupling coefficient times the group's common flow level.\nThe flow direction of each device is inferred from its directional capacities: a device with <code>consumption-capacity: \"0 kW\"</code> is an output (producing) device, and a device with <code>production-capacity: \"0 kW\"</code> is an input (consuming) device.\nExactly one of the two directional capacities must be set to a fixed zero for each coupled device.\nFor example, a CHP unit with 50% thermal and 30% electrical efficiency uses a gas input device (coefficient 1), a heat output device (coefficient 0.5) and an electricity output device (coefficient 0.3).\nDefaults to 1.\n",
+            "example": 0.5
+          },
           "sensor": {
             "type": "integer",
             "description": "ID of the device's power sensor."
diff --git a/tests/documentation/test_schemas.py b/tests/documentation/test_schemas.py
index fab2946647..ead6fb0a9e 100644
--- a/tests/documentation/test_schemas.py
+++ b/tests/documentation/test_schemas.py
@@ -10,11 +10,11 @@
 
 # Metadata constants that intentionally do not appear in the documentation
 EXCLUDED_METADATA = {
+    "COMMODITY_FLEX_CONTEXT",  # appears as `commodity` in the flex-context listing in scheduling.rst
+    "COMMODITY_FLEX_MODEL",  # appears as `commodity` in the flex-model listing in scheduling.rst
     "RELAX_CAPACITY_CONSTRAINTS",
     "RELAX_SITE_CAPACITY_CONSTRAINTS",
     "RELAX_SOC_CONSTRAINTS",
-    "COMMODITY_FLEX_CONTEXT",  # Documented as "commodity" in flex-context section
-    "COMMODITY_FLEX_MODEL",  # Documented as "commodity" in flex-model section
 }