diff --git a/documentation/changelog.rst b/documentation/changelog.rst index eda7bcad36..86e3a06d59 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -35,6 +35,7 @@ New features * The flex-context can now define multiple commodities, each specifying their own prices and grid capacities [see `PR #1946 `_, `PR #2172 `_, `PR #2235 `_, `PR #2271 `_, `PR #2355 `_ and `PR #2380 `_] * In the UI, the flex-context editor supports editing commitments (name, baseline and deviation prices, each accepting a fixed value or a sensor), also within each commodity context; the commitment's commodity follows the commodity tab being edited, and new commitments start with zero-valued baseline and prices [see `PR #2287 `_] * A commitment in the flex-context now requires a ``baseline`` and at least one deviation price (``up-price`` and/or ``down-price``), as already documented; commitment costs are reported in the scheduling results under the user-given name (with a ``(custom)`` suffix in the rare case the name collides with a scheduler-internal commitment name) [see `PR #2287 `_] +* 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 `_] * 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 bfbec970d3..9db8fb7732 100644 --- a/documentation/features/scheduling.rst +++ b/documentation/features/scheduling.rst @@ -229,6 +229,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 092586c2bc..de2e97f9a0 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 @@ -32,7 +33,6 @@ from flexmeasures.data.models.time_series import Sensor from flexmeasures.data.models.generic_assets import GenericAsset as Asset - #: The flex-context keys (as stored/passed, i.e. with hyphens) that together define #: the inflexible devices. They form one field family: the deprecated key must not #: be mixed with the newer sign-explicit keys in a single flex-context, and a @@ -59,8 +59,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). @@ -97,6 +99,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 #: For inflexible devices: whether positive sensor values mean consumption, #: as determined by the flex-context field the device was listed under #: (``inflexible-consumption`` → True, ``inflexible-production`` → False). @@ -195,6 +203,48 @@ 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: an input device, internally positive coefficient. + - only a (non-zero) ``production_capacity`` flows: an 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 _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: @@ -541,6 +591,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) @@ -662,6 +714,27 @@ 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 + @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/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 0dbae6eb3e..7312f77ca7 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -223,6 +223,7 @@ def device_scheduler_highspy( # noqa C901 stock_groups: dict[int, list[int]] | None = None, ems_constraint_groups: list[list[int]] | None = None, device_power_bands: list[list[tuple[float, float]] | None] | None = None, + coupling_groups: dict[str, list[tuple[int, float]]] | None = None, ) -> tuple[list[pd.Series], float, HighspySolverResults, HighspyModel]: """Direct HiGHS implementation of ``device_scheduler``. @@ -253,6 +254,7 @@ def device_scheduler_highspy( # noqa C901 stock_groups=stock_groups, ems_constraint_groups=ems_constraint_groups, device_power_bands=device_power_bands, + coupling_groups=coupling_groups, ) # Local aliases, so that the model below reads as it did before the (solver-agnostic) @@ -269,6 +271,7 @@ def device_scheduler_highspy( # noqa C901 convex_cost_curve = problem.convex_cost_curve Md, Mc = problem.Md, problem.Mc band_lookup = problem.band_lookup + coupling_device_specs = problem.coupling_device_specs _initial_stock_of = problem.initial_stock_of # --------------------------------------------------------------- @@ -347,6 +350,11 @@ def device_scheduler_highspy( # noqa C901 band_col[(d, b)] = col_band + i * T ncol += len(band_pairs) * T + # coupling_alpha[g, j]: the group's common normalised flow (free variable). + n_coupling_groups = 1 + max((g for g, _, _ in coupling_device_specs), default=-1) + col_alpha = ncol + ncol += n_coupling_groups * T + lower = np.full(ncol, -infinity) upper = np.full(ncol, infinity) cost = np.zeros(ncol) @@ -627,6 +635,20 @@ def _add_commitment_rows(c, quantity, jj, lb, ub, devices, is_stock) -> None: ), ) + # flow_coupling_constraints: pin every coupled device's flow to the group's + # common normalised flow, scaled by its coefficient: + # ems_power[d, j] - coeff * alpha[g, j] == 0 + # The coefficient's sign gives the direction (positive for inputs, negative for outputs). + if coupling_device_specs: + js = np.arange(T) + for g, d, coeff in coupling_device_specs: + rows.add_uniform_rows( + np.zeros(T), + np.zeros(T), + np.column_stack([col_ems + int(d) * T + js, col_alpha + g * T + js]), + np.tile([1.0, -float(coeff)], (T, 1)), + ) + # --------------------------------------------------------------- # Build and solve the HiGHS model # --------------------------------------------------------------- diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 51b94a0957..f174b5ca3a 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -133,6 +133,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, device_power_bands: list[list[tuple[float, float]] | None] | None = None, ) -> tuple[list[pd.Series], float, SolverResults, ConcreteModel]: @@ -172,6 +173,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)]} :param device_power_bands: optional per-device list of signed power bands (min, max), in flow units (e.g. MW, positive for consumption). A device with bands must operate within one of its bands at every time step (see S2 operation modes); this introduces @@ -226,6 +236,7 @@ def device_scheduler( # noqa C901 stock_groups=stock_groups, ems_constraint_groups=ems_constraint_groups, device_power_bands=device_power_bands, + coupling_groups=coupling_groups, ) # Local aliases, @@ -243,6 +254,7 @@ def device_scheduler( # noqa C901 Md, Mc = problem.Md, problem.Mc band_lookup = problem.band_lookup _initial_stock_of = problem.initial_stock_of + coupling_device_specs = problem.coupling_device_specs # Add indices for devices (d), datetimes (j) and commitments (c) model.d = RangeSet(0, len(device_constraints) - 1, doc="Set of devices") @@ -640,6 +652,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 + ) + # Power bands (S2 operation modes): a banded device must operate within # exactly one of its declared signed power ranges at every time step. # Which band it runs in (per time step) is a free binary decision variable; diff --git a/flexmeasures/data/models/planning/scheduling_problem.py b/flexmeasures/data/models/planning/scheduling_problem.py index 92b48d349d..840a5b98a1 100644 --- a/flexmeasures/data/models/planning/scheduling_problem.py +++ b/flexmeasures/data/models/planning/scheduling_problem.py @@ -219,6 +219,9 @@ class SchedulingProblem: #: device index -> its signed power bands (S2 operation modes) band_lookup: dict[int, list[tuple[float, float]]] + #: (group index, device index, coefficient) triples for hard flow-coupling constraints + coupling_device_specs: list[tuple[int, int, float]] + initial_stock: float | list[float] #: The commitments as passed in, before the sub-commitment split. @@ -268,6 +271,7 @@ def prepare_scheduling_problem( # noqa C901 stock_groups: dict[int, list[int]] | None = None, ems_constraint_groups: list[list[int]] | None = None, device_power_bands: list[list[tuple[float, float]] | None] | None = None, + coupling_groups: dict[str, list[tuple[int, float]]] | None = None, ) -> SchedulingProblem: """Normalise and validate ``device_scheduler``'s arguments into a SchedulingProblem. @@ -302,22 +306,29 @@ def prepare_scheduling_problem( # 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, + # 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[str, 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. @@ -349,6 +360,15 @@ def prepare_scheduling_problem( # 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 = [] @@ -525,6 +545,7 @@ def prepare_scheduling_problem( # noqa C901 Md=Md, Mc=Mc, band_lookup=band_lookup, + coupling_device_specs=coupling_device_specs, initial_stock=initial_stock, original_commitments=original_commitments, ) diff --git a/flexmeasures/data/models/planning/storage.py b/flexmeasures/data/models/planning/storage.py index 70cb4a01d3..bf68ed9a99 100644 --- a/flexmeasures/data/models/planning/storage.py +++ b/flexmeasures/data/models/planning/storage.py @@ -64,7 +64,6 @@ from flexmeasures.utils.time_utils import determine_minimum_resampling_resolution from flexmeasures.utils.unit_utils import ur, convert_units, units_are_convertible - storage_asset_types = ["one-way_evse", "two-way_evse", "battery", "heat-storage"] @@ -176,6 +175,10 @@ 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 + # 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 @@ -3326,6 +3329,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, device_power_bands=[ dc.attrs.get("operation_modes") for dc in device_constraints ], diff --git a/flexmeasures/data/models/planning/tests/test_commitments.py b/flexmeasures/data/models/planning/tests/test_commitments.py index 90663dd847..e9d94de23e 100644 --- a/flexmeasures/data/models/planning/tests/test_commitments.py +++ b/flexmeasures/data/models/planning/tests/test_commitments.py @@ -1575,6 +1575,667 @@ 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 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 67eec9e8f2..574a6986ed 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: @@ -372,6 +374,44 @@ def test_stock_constraint_device(): 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 test_inflexible_field_enumeration_and_signs(): """Inflexible devices are enumerated per field in canonical order, carrying the sign convention of the field they were listed under (None = defer to attribute).""" diff --git a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py index b1459cd92b..90407d91ad 100644 --- a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py +++ b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py @@ -317,6 +317,40 @@ def scenario_ems_level_commodity_commitment(): ) +def make_flow_device(derivative_min: float, derivative_max: float) -> pd.DataFrame: + """A device with no stock limits, free to flow within the given band.""" + device_constraints = initialize_df(COLUMNS, START, END, RESOLUTION) + device_constraints["min"] = -1000 + device_constraints["max"] = 1000 + device_constraints["derivative min"] = derivative_min + device_constraints["derivative max"] = derivative_max + return device_constraints + + +def scenario_chp_coupling_groups(): + """A CHP whose gas input, heat output and power output are hard-coupled. + + Every device's flow is pinned to the group's common normalised flow, scaled by its coefficient, + so the direct backend has to model coupling_groups or the three devices move independently. + """ + index = make_index() + prices = make_prices(index) + return dict( + device_constraints=[ + make_flow_device(0, 1), # gas in + make_flow_device(-1, 0), # heat out + make_flow_device(-1, 0), # power out + ], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=[ + make_energy_commitment(index, prices, devices=2), + make_energy_commitment(index, -prices / 4, devices=0), + ], + initial_stock=[0.0, 0.0, 0.0], + coupling_groups={"chp": [(0, 1.0), (1, -0.5), (2, -0.3)]}, + ) + + @pytest.mark.parametrize( "make_scenario", [ @@ -326,6 +360,7 @@ def scenario_ems_level_commodity_commitment(): scenario_two_devices_with_stock_commitment, scenario_ems_level_flow_commitment, scenario_ems_level_commodity_commitment, + scenario_chp_coupling_groups, ], ids=lambda f: f.__name__.replace("scenario_", ""), ) diff --git a/flexmeasures/data/models/planning/tests/test_storage.py b/flexmeasures/data/models/planning/tests/test_storage.py index 7b9391839f..1589a47280 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 @@ -1976,6 +1977,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 ae640c2f94..60c63ebad4 100644 --- a/flexmeasures/data/schemas/scheduling/__init__.py +++ b/flexmeasures/data/schemas/scheduling/__init__.py @@ -1406,6 +1406,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"], + }, "inflexible-consumption": { "default": None, "description": rst_to_openapi(metadata.INFLEXIBLE_CONSUMPTION.description), diff --git a/flexmeasures/data/schemas/scheduling/metadata.py b/flexmeasures/data/schemas/scheduling/metadata.py index fdd9a8a1d4..c4db606f3e 100644 --- a/flexmeasures/data/schemas/scheduling/metadata.py +++ b/flexmeasures/data/schemas/scheduling/metadata.py @@ -255,6 +255,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 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, +) 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 601d90dfe0..87ddc83ca0 100644 --- a/flexmeasures/data/schemas/scheduling/storage.py +++ b/flexmeasures/data/schemas/scheduling/storage.py @@ -59,6 +59,19 @@ ) +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", + ) + + def validate_inflexible_flex_model_entry(data: dict, original_data: dict): """Validate a flex-model entry that declares an inflexible device. @@ -484,6 +497,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, @@ -616,6 +642,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.""" @@ -798,6 +869,21 @@ 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(), + ) + mapped_schema_keys: dict def __init__(self, *args, **kwargs): @@ -818,6 +904,10 @@ def validate_group(self, group: dict, **kwargs): def validate_inflexible_device(self, data: dict, original_data: dict, **kwargs): validate_inflexible_flex_model_entry(data, original_data) + @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 af4bdd8313..6192aaefc8 100644 --- a/flexmeasures/data/schemas/tests/test_scheduling.py +++ b/flexmeasures/data/schemas/tests/test_scheduling.py @@ -981,7 +981,9 @@ def test_flex_model_schemas( and the second entry represents the expectation for the DBStorageFlexModelSchema. """ schemas = [ - StorageFlexModelSchema(start=datetime(2026, 6, 1), sensor=None), + StorageFlexModelSchema( + start=datetime(2026, 6, 1, tzinfo=pytz.utc), sensor=None + ), DBStorageFlexModelSchema(), ] if not isinstance(fails, list): @@ -1547,6 +1549,101 @@ def test_asset_trigger_schema_rejects_malformed_flex_context(app): 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): + """A device with a `coupling` field must have an unambiguous flow direction. + + The direction is 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, tzinfo=pytz.utc), 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): + """The coupling-direction check only applies to devices that define a `coupling` field.""" + schema = StorageFlexModelSchema( + start=datetime(2026, 6, 1, tzinfo=pytz.utc), 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): + """A provided coupling name must contain at least one non-whitespace character. + + Otherwise unrelated devices could be silently coupled under an empty group key. + This holds for both the scheduling schema and the db-stored one. + """ + 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, tzinfo=pytz.utc), 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): + """A db-stored flex-model accepts `coupling`/`coupling-coefficient` and round-trips them. + + Such flex-models are validated via DBStorageFlexModelSchema, e.g. by patch_asset. + """ + 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}) + + # Note: AssetTriggerSchema itself no longer aliases legacy field names (e.g. # force_new_job_creation) -- that's v3_0-specific backward compatibility, # layered on top by AssetTriggerSchemaV3 in flexmeasures/api/v3_0/assets.py, diff --git a/flexmeasures/ui/static/openapi-specs.json b/flexmeasures/ui/static/openapi-specs.json index eb6b0cc638..ebe6f3a407 100644 --- a/flexmeasures/ui/static/openapi-specs.json +++ b/flexmeasures/ui/static/openapi-specs.json @@ -7147,6 +7147,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 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 + }, "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 }