diff --git a/documentation/changelog.rst b/documentation/changelog.rst index ed5b1019f1..c8b30b6a3c 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -22,6 +22,8 @@ New features * Let storage scheduling infer missing ``power-capacity`` from directional device capacities before falling back to site capacity, and default the missing opposite capacity to zero when only a non-zero ``consumption-capacity`` or ``production-capacity`` is configured [see `PR #2222 `_] * Support multiple feeders to a shared storage [see `PR #2001 `_, `PR #2321 `_, `PR #2322 `_ and `PR #2325 `_] * The flex-context can now define multiple commodities, each specifying their own prices and grid capacities [see `PR #1946 `_, `PR #2172 `_, `PR #2235 `_ and `PR #2271 `_] +* Support commodity-converting devices, such as a CHP, e-boiler or heat pump, by describing each of the device's commodity ports as a flex-model entry sharing one ``coupling`` group with fixed flow ratios (``coupling-coefficient``) [see `PR #2218 `_] +* Support unit-commitment of a ``coupling`` group, so a commodity-converting device (such as a cogeneration unit) can be modelled with an on/off binary, a minimum level when on (``coupling-min`` on the reference port), and a per-port no-load base gated by that binary (``coupling-base``), making the affine ``P = coefficient · level + base`` relation exact including no-load fuel [see `PR #2336 `_] * Commodity contexts that omit grid-connection fields (prices and site capacities) now get smart defaults instead of failing or silently leaving the grid unconstrained — for instance, a bare ``{"commodity": "gas"}`` is treated as having no grid connection; see :ref:`commodity_context_defaults` for the full rules [see `PR #2272 `_] * CLI support for adding/editing account attributes [see `PR #2242 `_] * Improve chart axis domain for event values not around zero, with a per-sub-chart ``y-axis`` option in ``sensors_to_show`` (default ``zero``, which pads the axis out to include zero) that can be set to ``data`` to fit a sub-chart's y-axis to the values shown, to an explicit ``[min, max]`` domain that the axis will cover at least (expanding to fit data beyond it), or to a strict ``{"min": min, "max": max}`` domain that the axis will never exceed (clamping data beyond it, with a warning when that happens), editable from the graph editor [see `PR #2244 `_] diff --git a/documentation/features/scheduling.rst b/documentation/features/scheduling.rst index 08e62e9267..fd24c487b3 100644 --- a/documentation/features/scheduling.rst +++ b/documentation/features/scheduling.rst @@ -223,6 +223,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 05bb7f8a9b..cdfb1f86a6 100644 --- a/flexmeasures/data/models/planning/devices.py +++ b/flexmeasures/data/models/planning/devices.py @@ -22,6 +22,7 @@ from __future__ import annotations +import math from dataclasses import dataclass, field from enum import Enum from functools import cached_property @@ -36,8 +37,10 @@ class DeviceRole(Enum): """The role a flex-model (or flex-context) entry plays in the scheduling problem. - Extension point (not yet implemented): 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`. + GROUP entries constrain the aggregate power of a set of member devices. """ #: A schedulable flexible device (usually with a power sensor). @@ -72,6 +75,19 @@ 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 + #: Signed per-port no-load base (in MW), gated by the group's on/off binary when the group is unit-committed. + #: Follows the same sign convention as the coefficient; 0.0 when no ``coupling-base`` is given. + coupling_base: float = 0.0 + #: Group minimum marginal level (in MW), declared on the reference port (|coefficient| == 1). None when not set. + coupling_min: float | None = None + #: The reference port's power capacity (in MW), used as the group's maximum marginal level. None when not resolvable. + coupling_max: float | None = None @property def sensor_id(self) -> int | None: @@ -155,6 +171,83 @@ 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 which directional capacity allows flow + (mirroring how a missing directional site/device capacity defaults to zero): + + - only a (non-zero) ``consumption_capacity`` flows -> input device -> + internally positive coefficient + - only a (non-zero) ``production_capacity`` flows -> output device -> + internally negative coefficient + + The unspecified direction is assumed to be zero, so the user no longer needs + to set the opposite direction to a fixed 0 (though doing so still works). + """ + coefficient = abs(float(flex_model.get("coupling_coefficient", 1.0))) + consumption = flex_model.get("consumption_capacity") + production = flex_model.get("production_capacity") + consumption_flows = consumption is not None and not _is_zero_capacity(consumption) + production_flows = production is not None and not _is_zero_capacity(production) + consumption_blocked = _is_zero_capacity(consumption) + production_blocked = _is_zero_capacity(production) + # A direction is active if it flows itself, or if the opposite direction is + # explicitly pinned to zero (the legacy way of marking a direction). + consumption_active = consumption_flows or production_blocked + production_active = production_flows or consumption_blocked + if production_active and not consumption_active: + # Output (producing) device -> internally negative coefficient. + coefficient = -coefficient + return coefficient + + +def _quantity_to_mw(value: Any) -> float | None: + """Convert a fixed power quantity to a magnitude in MW; None if it is not a fixed scalar quantity. + + Sensor references and time series (used for time-varying capacities) return None: + a unit-committed coupling group needs scalar bounds, so those are handled by the caller. + """ + if value is None: + return None + if hasattr(value, "to") and hasattr(value, "magnitude"): + try: + return float(value.to("MW").magnitude) + except Exception: + return None + try: + return float(value) + except (TypeError, ValueError): + return None + + +def _resolve_coupling_base(base_value: Any, coefficient: float) -> float: + """Resolve a coupled port's signed no-load base (in MW) from its raw ``coupling-base`` value. + + The ``coupling-base`` field is a user-facing positive power magnitude; its internal + sign follows the port's flow direction, i.e. the sign of the port's coupling + coefficient (mirroring :func:`_resolve_coupling_coefficient`). Returns 0.0 when no + base is given (proportional, non-unit-committed behaviour). + """ + magnitude = _quantity_to_mw(base_value) + if magnitude is None: + return 0.0 + return math.copysign(abs(magnitude), coefficient) + + def _ref_id(value: Any) -> int | None: """Return the id of a sensor/asset reference, which may be a model object or a raw id.""" if value is None: @@ -341,6 +434,8 @@ def register_stock_params(stock_key: int, fm: dict) -> None: inventory.stock_entries[stock_key] = fm for fm in flex_model_list: + # Each flex-model entry is a deserialized dict. + assert isinstance(fm, dict) # Group entry (multi-device mode only): this entry's own sensor/asset is # the aggregate sensor/asset referenced by another entry's "group" field. if _classify_group_entry(inventory, fm): @@ -399,6 +494,15 @@ 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=( + coupling_coefficient := _resolve_coupling_coefficient(fm) + ), + coupling_base=_resolve_coupling_base( + fm.get("coupling_base"), coupling_coefficient + ), + coupling_min=_quantity_to_mw(fm.get("coupling_min")), + coupling_max=_quantity_to_mw(fm.get("power_capacity_in_mw")), ) inventory.entries.append(device) inventory.devices.append(device) @@ -486,6 +590,96 @@ def stock_constraint_device(self, stock_key: int) -> int | None: group_devices = self.stock_groups.get(stock_key) return group_devices[0] if group_devices else None + @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 + + def _coupling_reference_port(self, members: list[tuple[int, float]]) -> FlexDevice: + """Return a coupling group's reference port: the port with |coefficient| == 1. + + The reference port is the driving variable of the group; it carries the group's + ``coupling-min`` and the ``power-capacity`` that bounds the group's marginal level. + + :raises ValueError: When no member has a unit coefficient. + """ + for d_idx, coeff in members: + if math.isclose(abs(coeff), 1.0, abs_tol=1e-9): + return self.devices[d_idx] + raise ValueError( + "A unit-committed coupling group must have a reference port with" + " coupling-coefficient 1 (the driving variable). None was found." + ) + + @cached_property + def coupling_uc(self) -> dict[str, tuple[float, float]]: + """Map each unit-committed coupling group to its ``(min, max)`` marginal-level bounds (in MW). + + A coupling group is unit-committed when its reference port declares a + ``coupling-min``, or any of its ports declares a non-zero ``coupling-base``. + The minimum comes from the reference port's ``coupling-min`` (0 when only a base + is given); the maximum from the reference port's ``power-capacity``. + Empty for purely proportional coupling (leaving the problem an LP). + + :raises ValueError: When a unit-committed group lacks a resolvable power-capacity + on its reference port. + """ + result: dict[str, tuple[float, float]] = {} + for name, members in self.coupling_groups.items(): + group_devices = [self.devices[d_idx] for d_idx, _ in members] + declared_mins = [ + dev.coupling_min + for dev in group_devices + if dev.coupling_min is not None + ] + has_base = any(abs(dev.coupling_base) > 1e-12 for dev in group_devices) + if not declared_mins and not has_base: + continue + reference = self._coupling_reference_port(members) + min_level = declared_mins[0] if declared_mins else 0.0 + max_level = reference.coupling_max + if max_level is None: + raise ValueError( + f"Unit-committed coupling group '{name}' needs a fixed power-capacity" + " on its reference port (the port with coupling-coefficient 1) to bound" + " its marginal level." + ) + result[name] = (min_level, max_level) + return result + + @cached_property + def coupling_bases(self) -> dict[str, list[tuple[int, float]]]: + """Map each unit-committed coupling group to its ports' (device index, signed base) pairs (in MW). + + Restricted to the groups in :attr:`coupling_uc`; suitable for passing to + ``device_scheduler(coupling_bases=...)`` alongside ``coupling_groups``. + """ + uc_groups = self.coupling_uc + result: dict[str, list[tuple[int, float]]] = {} + for name, members in self.coupling_groups.items(): + if name not in uc_groups: + continue + result[name] = [ + (d_idx, self.devices[d_idx].coupling_base) for d_idx, _ in members + ] + return result + @cached_property def group_to_devices(self) -> dict[tuple[str, int], list[int]]: """Map each group key to the indices of the (leaf) member devices of that group. diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index a53bc296ff..0fdc237f7e 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -43,6 +43,9 @@ 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, + coupling_uc: dict[str, tuple[float, float]] | None = None, + coupling_bases: 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, @@ -81,6 +84,25 @@ 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)]} + + :param coupling_uc: Optional unit-commitment bounds per coupling group, marking that group unit-committed. + Each entry maps a group name to a ``(min, max)`` pair for the group's marginal level ``alpha``. + A per-time-step binary ``u`` is introduced and ``min * u <= alpha[group, j] <= max * u``, + so ``alpha`` is zero when the group is off and at least ``min`` when on. + :param coupling_bases: Optional per-port no-load base per coupling group. Each entry maps a group name to a list of + ``(device_index, base)`` tuples. For a unit-committed group, every port ``d`` is then + constrained by ``P[d, j] == coeff_d * alpha[group, j] + base_d * u[group, j]``, so the + (signed) no-load base is gated by the on/off binary (e.g. a cogeneration unit's no-load fuel). + Ignored for groups not present in ``coupling_uc``. Potentially deprecated arguments: commitment_quantities: amounts of flow specified in commitments (both previously ordered and newly requested) @@ -132,22 +154,29 @@ 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: dict[str, list[int]] = {} # Group keys are namespaced strings: a declared stock group's key (a state-of-charge # sensor id) could otherwise collide with the device index of an ungrouped device, # silently merging that device into the stock group. + # + # A device may belong to more than one stock group — a commodity converter (e.g. a + # steamer bridging a heat node and a steam node) participates in every node it + # touches, so ``group_to_devices`` keeps the full (possibly overlapping) membership. + # ``device_to_group`` records only the primary group (first assignment wins), used + # where a single owning group is needed (per-device stock bounds). if stock_groups: for g, devices in stock_groups.items(): + gkey = f"stock:{g}" + group_to_devices[gkey] = list(devices) for d in devices: - device_to_group[d] = f"stock:{g}" + device_to_group.setdefault(d, gkey) # Devices not in any stock group (e.g. inflexible devices) form individual groups. for d in range(len(device_constraints)): if d not in device_to_group: - device_to_group[d] = f"device:{d}" - - group_to_devices: dict[int, list[int]] = {} - for d, g in device_to_group.items(): - group_to_devices.setdefault(g, []).append(d) + gkey = f"device:{d}" + device_to_group[d] = gkey + group_to_devices[gkey] = [d] # The stock recursion is modelled once per stock group, using the group's shared # storage efficiency, so devices sharing a stock may not declare different ones. @@ -179,6 +208,15 @@ def device_scheduler( # noqa C901 " stock, so define it once per stock group." ) + # 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: commitments = [] @@ -780,6 +818,79 @@ def device_derivative_equalities(m, d, j): model.d, model.j, rule=device_derivative_equalities ) + if coupling_device_specs: + # coupling_device_specs is only populated when coupling_groups is given. + assert coupling_groups is not None + group_names = list(coupling_groups.keys()) + 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) + + # Unit-commitment: a group listed in ``coupling_uc`` gains a per-time-step + # on/off binary. Its marginal level ``alpha`` is then gated by that binary + # (min * u <= alpha <= max * u), and every port picks up a no-load base + # (base * u) on top of the proportional term. Groups without UC keep their + # purely proportional behaviour (the model stays an LP unless some group is + # unit-committed). + coupling_uc = coupling_uc or {} + coupling_bases = coupling_bases or {} + uc_bounds_by_gidx: dict[int, tuple[float, float]] = {} + base_by_group_device: dict[tuple[int, int], float] = {} + for g_idx, name in enumerate(group_names): + if name in coupling_uc: + uc_bounds_by_gidx[g_idx] = coupling_uc[name] + for d_idx, base in coupling_bases.get(name, []): + base_by_group_device[(g_idx, d_idx)] = base + + if uc_bounds_by_gidx: + model.coupling_uc_range = Set(initialize=sorted(uc_bounds_by_gidx)) + model.coupling_on = Var( + model.coupling_uc_range, model.j, domain=Binary, initialize=0 + ) + + def coupling_alpha_lower_rule(m, g, j): + """alpha is at least the group minimum when the group is on, else zero.""" + min_level, _ = uc_bounds_by_gidx[g] + return m.coupling_alpha[g, j] >= min_level * m.coupling_on[g, j] + + def coupling_alpha_upper_rule(m, g, j): + """alpha is at most the group maximum when on, and forced to zero when off.""" + _, max_level = uc_bounds_by_gidx[g] + return m.coupling_alpha[g, j] <= max_level * m.coupling_on[g, j] + + model.coupling_alpha_lower_bounds = Constraint( + model.coupling_uc_range, model.j, rule=coupling_alpha_lower_rule + ) + model.coupling_alpha_upper_bounds = Constraint( + model.coupling_uc_range, model.j, rule=coupling_alpha_upper_rule + ) + + def flow_coupling_rule(m, c, j): + """Enforce the (affine) coupling of each port to its group's level. + + For a purely proportional group: ``P[d, j] == coeff * alpha[group, j]``. + For a unit-committed group: ``P[d, j] == coeff * alpha[group, j] + base * u[group, j]``, + so the (signed) no-load base is gated by the on/off binary. + The coefficient sign indicates direction: positive for inputs (consuming), + negative for outputs (producing); the base follows the same sign convention. + """ + g, d, coeff = coupling_device_specs[c] + if g in uc_bounds_by_gidx: + base = base_by_group_device.get((g, d), 0.0) + return ( + m.ems_power[d, j] + == coeff * m.coupling_alpha[g, j] + base * m.coupling_on[g, j] + ) + 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 2c2211e439..093b5a2032 100644 --- a/flexmeasures/data/models/planning/storage.py +++ b/flexmeasures/data/models/planning/storage.py @@ -172,6 +172,15 @@ def _prepare(self, skip_validation: bool = False) -> tuple: # noqa: C901 for d in group_devices } + # 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 + # Unit-commitment of a coupling group: per-group (min, max) marginal-level bounds + # and per-port no-load bases, gated by a per-time-step on/off binary. Empty for + # purely proportional coupling (keeping the problem an LP). + self.coupling_uc = inventory.coupling_uc + self.coupling_bases = inventory.coupling_bases + # Group entries (intermediate power constraints on groups of devices, e.g. a # sub-EMS) come classified from the inventory, together with the resolved # (leaf) group membership. Accessing `group_to_devices` also detects cyclic @@ -3172,6 +3181,9 @@ 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, + coupling_uc=self.coupling_uc if self.coupling_uc else None, + coupling_bases=self.coupling_bases if self.coupling_bases 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 7c37334299..e04616b8b6 100644 --- a/flexmeasures/data/models/planning/tests/test_commitments.py +++ b/flexmeasures/data/models/planning/tests/test_commitments.py @@ -1575,6 +1575,859 @@ def test_simulation_with_dynamic_consumption_capacity(app, db): ), "Electric heater should charge the heat buffer at full capacity just 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 _run_unit_committed_cogeneration_scenario(elec_price: float, gas_price: float): + """Run a 1-step unit-committed cogeneration dispatch and return (schedules, cost). + + A single cogeneration unit is modelled as an affine, unit-committed coupling group + with the electrical output as its driving variable (the reference port, coefficient 1): + + P_power = -alpha (electricity produced, coeff -1, no base) + P_gas = SLOPE_GAS * alpha + BASE_GAS (gas consumed, coeff +2, base +3 kW) + P_heat = -SLOPE_HEAT * alpha - BASE_HEAT (heat produced, coeff -1.5, base -1 kW) + + with the unit either off (alpha = 0, u = 0, so every port is exactly 0 and no + no-load fuel is burned) or on (alpha = P in [PMIN, PMAX], u = 1). PMIN comes from + ``coupling_uc`` and PMAX from the reference port's power capacity. + + The economics: producing electricity earns ``elec_price`` per kW, consuming gas + costs ``gas_price`` per kW. Both are modelled as two-sided soft flow commitments at + quantity 0, so an upward deviation (consumption) costs and a downward deviation + (production) earns. + """ + PMIN, PMAX = 4.0, 10.0 # kW electrical output when running + SLOPE_GAS, BASE_GAS = 2.0, 3.0 # gas = 2 * P_elec + 3 (no-load fuel) + SLOPE_HEAT, BASE_HEAT = 1.5, 1.0 # heat = 1.5 * P_elec + 1 (no-load heat) + + start = pd.Timestamp("2026-01-01T00:00+01:00") + end = pd.Timestamp("2026-01-01T01:00+01:00") + resolution = pd.Timedelta("1h") + index = initialize_index(start=start, end=end, resolution=resolution) + + def _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, + "stock delta": 0.0, + } + defaults.update(kwargs) + return pd.DataFrame(defaults, index=index) + + # d=0 electrical output (reference port): production only, P in [-PMAX, 0]. + # d=1 gas input: consumption only, P in {0} U [SLOPE_GAS*PMIN+BASE_GAS, SLOPE_GAS*PMAX+BASE_GAS]. + # d=2 heat output: production only; wide enough to hold the affine value. + device_constraints = [ + _df(**{"derivative min": -PMAX, "derivative max": 0.0}), + _df(**{"derivative min": 0.0, "derivative max": SLOPE_GAS * PMAX + BASE_GAS}), + _df( + **{ + "derivative min": -(SLOPE_HEAT * PMAX + BASE_HEAT), + "derivative max": 0.0, + } + ), + ] + + ems_constraints = pd.DataFrame( + {"derivative min": -1000.0, "derivative max": 1000.0}, index=index + ) + + # Coupling: electrical output is the reference port (coeff -1, |coeff| == 1). + coupling_groups = { + "cogen": [(0, -1.0), (1, SLOPE_GAS), (2, -SLOPE_HEAT)], + } + # Unit commitment: marginal level (alpha = electrical output) bounded to [PMIN, PMAX] + # when on, forced to 0 when off. + coupling_uc = {"cogen": (PMIN, PMAX)} + # Signed no-load bases (same sign convention as the coefficients). + coupling_bases = { + "cogen": [(0, 0.0), (1, BASE_GAS), (2, -BASE_HEAT)], + } + + elec_p = pd.Series(elec_price, index=index) + gas_p = pd.Series(gas_price, index=index) + commitments = [ + FlowCommitment( + name="electricity", + index=index, + quantity=pd.Series(0.0, index=index), + upwards_deviation_price=elec_p, + downwards_deviation_price=elec_p, + device=pd.Series(0, index=index), + ), + FlowCommitment( + name="gas", + index=index, + quantity=pd.Series(0.0, index=index), + upwards_deviation_price=gas_p, + downwards_deviation_price=gas_p, + device=pd.Series(1, index=index), + ), + ] + + schedules, costs, results, _model = device_scheduler( + device_constraints=device_constraints, + ems_constraints=ems_constraints, + commitments=commitments, + coupling_groups=coupling_groups, + coupling_uc=coupling_uc, + coupling_bases=coupling_bases, + ) + assert results.solver.termination_condition == "optimal", ( + f"Solver did not find an optimal solution " + f"(elec_price={elec_price}, gas_price={gas_price})" + ) + return schedules, costs + + +def test_factory_unit_committed_cogeneration_dispatch(app, db): + """Factory: a unit-committed cogeneration unit idles or runs at/above its minimum. + + The unit's driving variable is its electrical output P, with P in {0} U [PMIN, PMAX]: + + P_power = -P, P_gas = 2*P + 3 (kW), P_heat = -(1.5*P + 1) (kW) + + PMIN = 4 kW, PMAX = 10 kW, no-load fuel = 3 kW, no-load heat = 1 kW. + + Profitable step (elec 50, gas 10) + --------------------------------- + Marginal profit per kW of electricity = elec - gas*slope_gas = 50 - 10*2 = 30 > 0, + so the unit runs at its maximum P = PMAX = 10 kW (higher output is always better, + and the no-load fuel is comfortably paid for): + + P_power = -10 kW, P_gas = 2*10 + 3 = 23 kW, P_heat = -(1.5*10 + 1) = -16 kW + + Cost = gas*P_gas - elec*|P_power| = 10*23 - 50*10 = 230 - 500 = -270 (a profit). + Running at PMIN instead would earn only 4*30 - 30 = 90, so PMAX is optimal. + + Unprofitable step (elec 10, gas 20) + ----------------------------------- + Marginal profit per kW = 10 - 20*2 = -30 < 0, so every kW loses money and there is + a fixed no-load fuel cost on top. Even at PMIN the unit would lose + 20*(2*4 + 3) - 10*4 = 220 - 40 = 180. Idling costs 0, so the unit fully idles: + ALL ports are exactly 0 and NO no-load fuel is burned (base is gated by u). + + The objective equals the hand-computed profit of the running step. + """ + # ----- unprofitable step: the unit fully idles ----- + idle_schedules, idle_cost = _run_unit_committed_cogeneration_scenario( + elec_price=10.0, gas_price=20.0 + ) + power_idle, gas_idle, heat_idle = idle_schedules + + # Every port is exactly 0 — in particular no no-load fuel is burned. + np.testing.assert_allclose( + power_idle.iloc[0], 0.0, atol=1e-6, err_msg="Idle: electrical output must be 0" + ) + np.testing.assert_allclose( + gas_idle.iloc[0], + 0.0, + atol=1e-6, + err_msg="Idle: gas input must be 0 (no-load fuel gated by the on/off binary)", + ) + np.testing.assert_allclose( + heat_idle.iloc[0], 0.0, atol=1e-6, err_msg="Idle: heat output must be 0" + ) + np.testing.assert_allclose( + idle_cost, 0.0, atol=1e-6, err_msg="Idle: objective must be 0" + ) + + # ----- profitable step: the unit runs at its maximum (>= PMIN) ----- + run_schedules, run_cost = _run_unit_committed_cogeneration_scenario( + elec_price=50.0, gas_price=10.0 + ) + power_run, gas_run, heat_run = run_schedules + + # P = PMAX = 10 kW; the affine relations (including the no-load bases) hold. + np.testing.assert_allclose( + power_run.iloc[0], + -10.0, + rtol=1e-4, + err_msg="Run: electrical output must be -10 kW (at PMAX, never between 0 and PMIN)", + ) + np.testing.assert_allclose( + gas_run.iloc[0], + 23.0, # 2 * 10 + 3 + rtol=1e-4, + err_msg="Run: gas input must be 2*P + 3 = 23 kW (affine with no-load base)", + ) + np.testing.assert_allclose( + heat_run.iloc[0], + -16.0, # -(1.5 * 10 + 1) + rtol=1e-4, + err_msg="Run: heat output must be -(1.5*P + 1) = -16 kW (affine with no-load base)", + ) + # Objective: 10*23 - 50*10 = -270. + np.testing.assert_allclose( + run_cost, -270.0, rtol=1e-4, err_msg="Run: objective must equal -270" + ) + + 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_device_inventory.py b/flexmeasures/data/models/planning/tests/test_device_inventory.py index 779b07fb73..77d6d23a36 100644 --- a/flexmeasures/data/models/planning/tests/test_device_inventory.py +++ b/flexmeasures/data/models/planning/tests/test_device_inventory.py @@ -12,7 +12,9 @@ from flexmeasures.data.models.planning.devices import ( DeviceInventory, DeviceRole, + _resolve_coupling_coefficient, ) +from flexmeasures.utils.unit_utils import ur def make_sensor(sensor_id: int, unit: str = "kW") -> Sensor: @@ -352,3 +354,109 @@ def test_stock_constraint_device(): device_c = inventory.devices[2] assert inventory.stock_constraint_device(device_c.stock_key) == 2 assert inventory.stock_constraint_device(999) is None + + +@pytest.mark.parametrize( + "capacities, expected_sign", + [ + # Smart default: only a consumption-capacity -> input -> positive. + ({"consumption_capacity": ur.Quantity("5 kW")}, 1), + # Smart default: only a production-capacity -> output -> negative. + ({"production_capacity": ur.Quantity("5 kW")}, -1), + # Explicit zero (back-compat): input device. + ( + { + "consumption_capacity": ur.Quantity("5 kW"), + "production_capacity": ur.Quantity("0 kW"), + }, + 1, + ), + # Explicit zero (back-compat): output device. + ( + { + "consumption_capacity": ur.Quantity("0 kW"), + "production_capacity": ur.Quantity("5 kW"), + }, + -1, + ), + # Only a fixed zero on the opposite side (legacy): production blocked -> input. + ({"production_capacity": ur.Quantity("0 kW")}, 1), + # Only a fixed zero on the opposite side (legacy): consumption blocked -> output. + ({"consumption_capacity": ur.Quantity("0 kW")}, -1), + ], +) +def test_resolve_coupling_coefficient_direction(capacities, expected_sign): + """The coupling coefficient's sign follows which directional capacity flows; + the unspecified direction defaults to zero (no explicit zero required).""" + coefficient = _resolve_coupling_coefficient( + {"coupling_coefficient": 0.5, **capacities} + ) + assert coefficient == expected_sign * 0.5 + + +def _cogeneration_flex_model(with_unit_commitment: bool) -> list[dict]: + """A cogeneration unit as three coupled ports (deserialized flex-model entries). + + The electrical output is the reference port (|coefficient| == 1, output -> -1); the + gas input and heat output are affine in the electrical output. When + ``with_unit_commitment`` is set, the reference port carries a ``coupling-min`` and + the gas/heat ports carry a signed ``coupling-base`` (their no-load offset). + """ + power_port = { + "sensor": make_sensor(1), + "commodity": "electricity", + "coupling": "cogen", + "coupling_coefficient": 1.0, + "production_capacity": ur.Quantity("10 kW"), # output -> coeff -1 + "power_capacity_in_mw": ur.Quantity("10 kW"), # bounds the group's max level + } + gas_port = { + "sensor": make_sensor(2), + "commodity": "gas", + "coupling": "cogen", + "coupling_coefficient": 2.0, + "consumption_capacity": ur.Quantity("30 kW"), # input -> coeff +2 + } + heat_port = { + "sensor": make_sensor(3), + "commodity": "heat", + "coupling": "cogen", + "coupling_coefficient": 1.5, + "production_capacity": ur.Quantity("20 kW"), # output -> coeff -1.5 + } + if with_unit_commitment: + power_port["coupling_min"] = ur.Quantity("4 kW") + gas_port["coupling_base"] = ur.Quantity("3 kW") + heat_port["coupling_base"] = ur.Quantity("1 kW") + return [power_port, gas_port, heat_port] + + +def test_coupling_unit_commitment_plumbing(): + """A unit-committed coupling group resolves to signed bases and (min, max) bounds in MW.""" + inventory = DeviceInventory.from_flex_config( + _cogeneration_flex_model(with_unit_commitment=True) + ) + # Proportional coefficients are unchanged (signed by flow direction). + assert inventory.coupling_groups == {"cogen": [(0, -1.0), (1, 2.0), (2, -1.5)]} + # Group bounds come from the reference port: min from coupling-min (4 kW), + # max from its power-capacity (10 kW), both converted to MW. + assert set(inventory.coupling_uc) == {"cogen"} + min_level, max_level = inventory.coupling_uc["cogen"] + assert min_level == pytest.approx(0.004) + assert max_level == pytest.approx(0.010) + # Per-port no-load bases are signed by flow direction and expressed in MW. + assert inventory.coupling_bases["cogen"] == [ + (0, pytest.approx(0.0)), # electrical output: no base + (1, pytest.approx(0.003)), # gas input: +3 kW no-load fuel + (2, pytest.approx(-0.001)), # heat output: -1 kW no-load heat + ] + + +def test_coupling_without_unit_commitment_stays_proportional(): + """Without coupling-min or coupling-base, a coupling group is purely proportional (no UC).""" + inventory = DeviceInventory.from_flex_config( + _cogeneration_flex_model(with_unit_commitment=False) + ) + assert inventory.coupling_groups == {"cogen": [(0, -1.0), (1, 2.0), (2, -1.5)]} + assert inventory.coupling_uc == {} + assert inventory.coupling_bases == {} diff --git a/flexmeasures/data/models/planning/tests/test_storage.py b/flexmeasures/data/models/planning/tests/test_storage.py index 9435e2c3db..92a3f69335 100644 --- a/flexmeasures/data/models/planning/tests/test_storage.py +++ b/flexmeasures/data/models/planning/tests/test_storage.py @@ -7,12 +7,12 @@ 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.soc_projection import ( project_off_tick_soc_at_start, project_off_tick_soc_constraints, ) -from flexmeasures.data.models.generic_assets import GenericAsset from flexmeasures.data.models.planning.storage import StorageScheduler from flexmeasures.data.models.planning.utils import initialize_index from flexmeasures.data.models.time_series import Sensor, TimedBelief @@ -21,6 +21,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 @@ -1965,6 +1966,190 @@ def test_resolve_soc_at_start_from_percent_sensor_uses_device_sensor_fallback( ) +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)", + ) + + def test_off_tick_soc_relaxation_covers_all_devices_of_a_shared_stock( add_battery_assets, db ): diff --git a/flexmeasures/data/schemas/scheduling/__init__.py b/flexmeasures/data/schemas/scheduling/__init__.py index fdfb71f151..eb2a42fd96 100644 --- a/flexmeasures/data/schemas/scheduling/__init__.py +++ b/flexmeasures/data/schemas/scheduling/__init__.py @@ -1223,6 +1223,24 @@ def check_prices(self, data: dict, original_data: dict, **kwargs): }, "example-units": EXAMPLE_UNIT_TYPES["commodity"], }, + "coupling": { + "default": None, + "description": rst_to_openapi(metadata.COUPLING.description), + "types": { + "backend": "typeOne", + "ui": "One fixed value only (a coupling-group name shared by a converter's commodity ports).", + }, + "example-units": ['a coupling-group name, e.g. "CHP"'], + }, + "coupling-coefficient": { + "default": 1.0, + "description": rst_to_openapi(metadata.COUPLING_COEFFICIENT.description), + "types": { + "backend": "typeOne", + "ui": "One fixed value only (a positive number).", + }, + "example-units": ["a number, e.g. 0.5"], + }, } diff --git a/flexmeasures/data/schemas/scheduling/metadata.py b/flexmeasures/data/schemas/scheduling/metadata.py index be4e9d5452..038f65fdfd 100644 --- a/flexmeasures/data/schemas/scheduling/metadata.py +++ b/flexmeasures/data/schemas/scheduling/metadata.py @@ -223,6 +223,42 @@ 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 which directional capacity is set: a device given only a ``production-capacity`` is an output (producing) device, and a device given only a ``consumption-capacity`` is an input (consuming) device. +The unspecified direction is assumed to be zero (mirroring how a missing directional site capacity defaults to zero), so there is no need to set the opposite direction to a fixed 0 (though setting it explicitly still works). +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, +) +COUPLING_BASE = MetaData( + description="""Per-port no-load base power for a unit-committed coupling group. +Give a positive power magnitude; its direction follows the port's flow direction (the same way ``coupling-coefficient`` does). +When the coupling group is unit-committed, each port's power equals its coupling coefficient times the group's marginal level *plus* this base, all gated by the group's on/off binary: ``P = coefficient * level + base`` when on, and ``0`` when off. +Use it to model a no-load offset, such as a cogeneration unit's no-load fuel consumption that is burned only while the unit runs. +Setting a non-zero base (on any port) makes the coupling group unit-committed. +Defaults to 0. +""", + example="3 kW", +) +COUPLING_MIN = MetaData( + description="""Minimum marginal level of a coupling group, declared on its reference port (the port with ``coupling-coefficient`` 1). +Setting it makes the coupling group unit-committed: a per-time-step on/off binary is introduced, and when the group runs its marginal level stays between this minimum and the reference port's ``power-capacity`` (the maximum); when it is off, all ports are exactly 0. +Use it to model a minimum load, such as a cogeneration unit that must run at or above a minimum output when on and be fully off otherwise. +""", + example="4 kW", +) 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 006a68ee28..2dd160bd66 100644 --- a/flexmeasures/data/schemas/scheduling/storage.py +++ b/flexmeasures/data/schemas/scheduling/storage.py @@ -44,6 +44,20 @@ def _validate_group_sensor_is_power_sensor(group: dict): ) +def _validate_coupling_name(coupling: str | None): + """Reject blank/whitespace-only coupling names. + + A blank coupling name would become a coupling-group key, silently coupling + unrelated devices under an empty group. When provided, the name must contain + at least one non-whitespace character. + """ + if coupling is not None and not coupling.strip(): + raise ValidationError( + "The `coupling` field, when provided, must be a non-empty (non-whitespace) name.", + field_name="coupling", + ) + + class GroupReferenceSchema(SharedSensorReferenceSchema): """Reference to a group of devices whose aggregate power is constrained. @@ -301,6 +315,33 @@ 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(), + ) + coupling_base = QuantityField( + "MW", + data_key="coupling-base", + required=False, + validate=validate.Range(min=ur.Quantity("0 MW")), + metadata=metadata.COUPLING_BASE.to_dict(), + ) + coupling_min = QuantityField( + "MW", + data_key="coupling-min", + required=False, + validate=validate.Range(min=ur.Quantity("0 MW")), + metadata=metadata.COUPLING_MIN.to_dict(), + ) def __init__( self, @@ -429,6 +470,51 @@ 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("coupling") + def validate_coupling(self, coupling: str | None, **kwargs): + _validate_coupling_name(coupling) + + @validates_schema + def validate_coupling_direction_is_unambiguous(self, data: dict, **kwargs): + """A coupled device must have an inferable flow direction. + + The flow direction is inferred from which directional capacity is given: + a device with (only) a consumption-capacity is an input (consuming) device, + and a device with (only) a production-capacity is an output (producing) + device. The unspecified direction is assumed to be zero, mirroring how a + missing directional site capacity defaults to zero, so the user does not + need to set the opposite direction to a fixed 0 (though doing so still works). + + The direction is ambiguous only when both directions are active (each side + either flows itself or is marked active by a fixed zero on the opposite side) + or when neither is (both missing); such flex-models are rejected. + """ + if data.get("coupling") is None: + return + + def _is_fixed_zero(value) -> bool: + return isinstance(value, ur.Quantity) and float(value.magnitude) == 0.0 + + def _flows(value) -> bool: + # A capacity flows when it is given and not a fixed zero. + # Sensor references cannot be checked statically, so they flow. + return value is not None and not _is_fixed_zero(value) + + consumption = data.get("consumption_capacity") + production = data.get("production_capacity") + # A direction is active if it flows itself, or if the opposite direction is + # explicitly pinned to zero (the legacy way of marking a direction). + consumption_active = _flows(consumption) or _is_fixed_zero(production) + production_active = _flows(production) or _is_fixed_zero(consumption) + if consumption_active == production_active: + raise ValidationError( + "A device with a 'coupling' field must have an unambiguous flow direction: " + "provide exactly one directional capacity, either a consumption-capacity " + "(for an input/consuming device) or a production-capacity (for an " + "output/producing device). The opposite direction defaults to zero.", + field_name="coupling", + ) + @post_load def post_load_sequence(self, data: dict, **kwargs) -> dict: """Perform some checks and corrections after we loaded.""" @@ -598,6 +684,37 @@ class DBStorageFlexModelSchema(Schema): metadata=dict(description="Commodity label for this device/asset."), ) + 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(), + ) + + coupling_base = QuantityField( + "MW", + data_key="coupling-base", + required=False, + validate=validate.Range(min=ur.Quantity("0 MW")), + metadata=metadata.COUPLING_BASE.to_dict(), + ) + + coupling_min = QuantityField( + "MW", + data_key="coupling-min", + required=False, + validate=validate.Range(min=ur.Quantity("0 MW")), + metadata=metadata.COUPLING_MIN.to_dict(), + ) + mapped_schema_keys: dict def __init__(self, *args, **kwargs): @@ -614,6 +731,10 @@ def __init__(self, *args, **kwargs): def validate_group(self, group: dict, **kwargs): _validate_group_sensor_is_power_sensor(group) + @validates("coupling") + def validate_coupling(self, coupling: str | None, **kwargs): + _validate_coupling_name(coupling) + @validates_schema def forbid_time_series_specs(self, data: dict, **kwargs): """Do not allow time series specs for the flex-model fields saved in the db.""" diff --git a/flexmeasures/data/schemas/tests/test_scheduling.py b/flexmeasures/data/schemas/tests/test_scheduling.py index d3db401fcd..dc6bcc2cbd 100644 --- a/flexmeasures/data/schemas/tests/test_scheduling.py +++ b/flexmeasures/data/schemas/tests/test_scheduling.py @@ -1455,3 +1455,91 @@ 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), + # Smart default: only a consumption-capacity given -> input device + # (production defaults to zero), no explicit zero needed. + ({"consumption-capacity": "5 kW"}, False), + # Smart default: only a production-capacity given -> output device + # (consumption defaults to zero), no explicit zero needed. + ({"production-capacity": "5 kW"}, False), + # Neither direction given: 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 an unambiguous flow direction, inferred from which directional capacity is given + (the opposite direction defaults 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"}) + + +@pytest.mark.parametrize("blank_name", ["", " ", "\t", " \n "]) +def test_blank_coupling_name_is_rejected(app, blank_name): + """test_blank_coupling_name_is_rejected: a provided coupling name must contain at least + one non-whitespace character, so unrelated devices cannot be silently coupled under an + empty group key. This holds for both the scheduling and the db-stored schema.""" + scheduling_flex_model = { + "power-capacity": "20 kW", + "production-capacity": "0 kW", + "coupling": blank_name, + } + with pytest.raises(ValidationError) as e_info: + StorageFlexModelSchema(start=datetime(2026, 6, 1), sensor=None).load( + scheduling_flex_model + ) + assert "non-empty" in str(e_info.value) + + with pytest.raises(ValidationError) as e_info: + DBStorageFlexModelSchema().load({"coupling": blank_name}) + assert "non-empty" in str(e_info.value) + + +def test_db_flex_model_coupling_round_trips(app): + """test_db_flex_model_coupling_round_trips: a db-stored flex-model (validated via + DBStorageFlexModelSchema, e.g. by patch_asset) accepts `coupling`/`coupling-coefficient` + and round-trips them.""" + schema = DBStorageFlexModelSchema() + flex_model = { + "coupling": "chp", + "coupling-coefficient": 0.5, + } + loaded = schema.load(flex_model) + assert loaded["coupling"] == "chp" + assert loaded["coupling_coefficient"] == 0.5 + # coupling-coefficient must be strictly positive + with pytest.raises(ValidationError): + schema.load({"coupling": "chp", "coupling-coefficient": 0}) diff --git a/flexmeasures/ui/static/openapi-specs.json b/flexmeasures/ui/static/openapi-specs.json index 854f57305b..92ce838906 100644 --- a/flexmeasures/ui/static/openapi-specs.json +++ b/flexmeasures/ui/static/openapi-specs.json @@ -6392,6 +6392,34 @@ ], "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 coupling-coefficient 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 which directional capacity is set: a device given only a production-capacity is an output (producing) device, and a device given only a consumption-capacity is an input (consuming) device.\nThe unspecified direction is assumed to be zero (mirroring how a missing directional site capacity defaults to zero), so there is no need to set the opposite direction to a fixed 0 (though setting it explicitly still works).\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 + }, + "coupling-base": { + "type": "string", + "x-minimum": "0 MW", + "description": "Per-port no-load base power for a unit-committed coupling group.\nGive a positive power magnitude; its direction follows the port's flow direction (the same way coupling-coefficient does).\nWhen the coupling group is unit-committed, each port's power equals its coupling coefficient times the group's marginal level plus this base, all gated by the group's on/off binary: P = coefficient * level + base when on, and 0 when off.\nUse it to model a no-load offset, such as a cogeneration unit's no-load fuel consumption that is burned only while the unit runs.\nSetting a non-zero base (on any port) makes the coupling group unit-committed.\nDefaults to 0.\n", + "example": "3 kW" + }, + "coupling-min": { + "type": "string", + "x-minimum": "0 MW", + "description": "Minimum marginal level of a coupling group, declared on its reference port (the port with coupling-coefficient 1).\nSetting it makes the coupling group unit-committed: a per-time-step on/off binary is introduced, and when the group runs its marginal level stays between this minimum and the reference port's power-capacity (the maximum); when it is off, all ports are exactly 0.\nUse it to model a minimum load, such as a cogeneration unit that must run at or above a minimum output when on and be fully off otherwise.\n", + "example": "4 kW" + }, "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 }