diff --git a/documentation/changelog.rst b/documentation/changelog.rst index 4074c3e9f4..aa1166e382 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -16,6 +16,7 @@ New features Infrastructure / Support ---------------------- * Upgraded dependencies [see `PR #1485 `_ and `PR #2215 `_] +* Prepare the ``device_scheduler`` to deal with commitments per device group [see `PR #1934 `_] Bugfixes ----------- diff --git a/flexmeasures/data/models/planning/__init__.py b/flexmeasures/data/models/planning/__init__.py index 46938c7eb2..997df93e4d 100644 --- a/flexmeasures/data/models/planning/__init__.py +++ b/flexmeasures/data/models/planning/__init__.py @@ -1,9 +1,12 @@ from __future__ import annotations +from collections.abc import Iterable from dataclasses import dataclass, field from datetime import datetime, timedelta +from tabulate import tabulate from typing import Any, Type +import numpy as np import pandas as pd from flask import current_app @@ -279,6 +282,7 @@ class Commitment: name: str device: pd.Series = None + device_group: pd.Series = None index: pd.DatetimeIndex = field(repr=False, default=None) _type: str = field(repr=False, default="each") group: pd.Series = field(init=False) @@ -287,10 +291,12 @@ class Commitment: downwards_deviation_price: pd.Series = 0 def __post_init__(self): + # device_group is a device→label lookup table, not a time series; + # exclude it from automatic time-series index coercion. series_attributes = [ attr for attr, _type in self.__annotations__.items() - if _type == "pd.Series" and hasattr(self, attr) + if _type == "pd.Series" and hasattr(self, attr) and attr != "device_group" ] for series_attr in series_attributes: val = getattr(self, series_attr) @@ -361,10 +367,67 @@ def __post_init__(self): "downwards deviation price" ) self.group = self.group.rename("group") + self._init_device_group() + + def _init_device_group(self): + # Extract device universe (empty for EMS-level commitments where device was None) + if isinstance(self.device, pd.Series): + devices = extract_devices(self.device) + else: + devices = [self.device] + + devices = list(devices) + + # EMS-level commitment: no device assignments, leave device_group empty. + if not devices: + self.device_group = pd.Series(dtype=object, name="device_group") + return + + # Default: one group per device (backwards compatible) + if self.device_group is None: + self.device_group = pd.Series( + range(len(devices)), index=devices, name="device_group" + ) + else: + # Validate custom grouping + missing = set(devices) - set(self.device_group.index) + if missing: + raise ValueError( + f"device_group missing assignments for devices: {missing}" + ) + self.device_group = self.device_group.loc[devices] + self.device_group.name = "device_group" + + def pretty_print(self): + """ + Pretty-print a list of Commitment objects as tabulated pandas DataFrames. + + For each Commitment, a DataFrame indexed by time is created containing + the commitment name, device values, group index, quantity, and any available + upward or downward deviation prices. Each commitment is printed separately + in a readable table format, making this function suitable for debugging, + logging, and interactive inspection. + """ + df = self.to_frame() + df = pd.DataFrame(index=df.device.index) + + df["commitment"] = self.name + df["device"] = self.device + df["group"] = self.group + df["quantity"] = self.quantity + + if hasattr(self, "upwards_deviation_price"): + df["up_price"] = self.upwards_deviation_price + + if hasattr(self, "downwards_deviation_price"): + df["down_price"] = self.downwards_deviation_price + + if not df.empty: + print(tabulate(df, headers=df.columns, tablefmt="fancy_grid")) def to_frame(self) -> pd.DataFrame: """Contains all info apart from the name.""" - return pd.concat( + df = pd.concat( [ self.device, self.quantity, @@ -375,6 +438,19 @@ def to_frame(self) -> pd.DataFrame: ], axis=1, ) + # Map device → device_group. + # For EMS-level commitments (device was None, now a NaN series) there are + # no device assignments, so device_group is an empty Series and we must not + # call map_device_to_group (it would KeyError on the NaN values). + devices = extract_devices(self.device) + if devices: + df["device_group"] = map_device_to_group(self.device, self.device_group) + else: + df["device_group"] = ( + np.nan + ) # EMS-level: handled by ems_flow_commitment_equalities + + return df class FlowCommitment(Commitment): @@ -396,3 +472,48 @@ class StockCommitment(Commitment): Scheduler.compute_schedule = deprecated(Scheduler.compute, "0.14")( Scheduler.compute_schedule ) + + +def extract_devices(device): + """ + Return a flat list of unique device identifiers from: + - scalar device + - Series of scalars + - Series of iterables (e.g. [0, 1]) + """ + if device is None: + return [] + + if isinstance(device, pd.Series): + values = device.dropna().values + else: + values = [device] + + devices = set() + for v in values: + if isinstance(v, Iterable) and not isinstance(v, (str, bytes)): + devices.update(v) + else: + devices.add(v) + + return list(devices) + + +def map_device_to_group(device_series, device_group_map): + """ + Map device identifiers to device_group. + + - scalar device → group label + - iterable of devices → group label (must be identical) + """ + + def resolve(v): + if isinstance(v, (list, tuple, set)): + groups = {device_group_map[d] for d in v} + if len(groups) != 1: + raise ValueError(f"Devices {v} map to multiple device groups: {groups}") + return groups.pop() + else: + return device_group_map[v] + + return device_series.apply(resolve) diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 93b919d6bb..1b85452f9d 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -207,6 +207,38 @@ def convert_commitments_to_subcommitments( commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) + device_group_lookup = {} + + for c, df in enumerate(commitments): + if "device" not in df.columns: + # EMS-level commitment: no device grouping needed here; + # handled by ems_flow_commitment_equalities. + continue + + has_device_group = "device_group" in df.columns + if has_device_group: + rows = df[["device", "device_group"]].dropna() + else: + # Backwards-compatible default: each device is its own group. + # This preserves the behaviour of old-style DataFrame commitments that + # pre-date the device_group feature (e.g. from initialize_device_commitment). + rows = df[["device"]].dropna() + + device_group_lookup[c] = {} + + for _, row in rows.iterrows(): + d = row["device"] + # When no device_group column is present, use the device id itself as + # the group label so that each device forms an independent group. + g = row["device_group"] if has_device_group else d + + if isinstance(d, (list, tuple, set, np.ndarray)): + devices = set(d) + else: + devices = {d} + + device_group_lookup[c].setdefault(g, set()).update(devices) + # Oversimplified check for a convex cost curve df = pd.concat(commitments)[ ["upwards deviation price", "downwards deviation price"] @@ -243,13 +275,24 @@ def convert_commitments_to_subcommitments( ) model.c = RangeSet(0, len(commitments) - 1, doc="Set of commitments") - # Add 2D indices for commitment datetimes (cj) + # Add 2D indices for commitment device groups (cg) + def commitment_device_groups_init(m): + return ((c, g) for c, groups in device_group_lookup.items() for g in groups) + + model.cg = Set(dimen=2, initialize=commitment_device_groups_init) + # Add 2D indices for commitment datetimes (cj) def commitments_init(m): return ((c, j) for c in m.c for j in commitments[c]["j"]) model.cj = Set(dimen=2, initialize=commitments_init) + # Add 3D indices for commitment datetime device groups (cjg) + def commitment_time_device_groups_init(m): + return ((c, j, g) for (c, j) in m.cj for (_, g) in m.cg if _ == c) + + model.cjg = Set(dimen=3, initialize=commitment_time_device_groups_init) + # Add parameters def price_down_select(m, c): if "downwards deviation price" not in commitments[c].columns: @@ -362,6 +405,40 @@ def device_derivative_up_efficiency(m, d, j): def device_stock_delta(m, d, j): return device_constraints[d]["stock delta"].iloc[j] + def grouped_commitment_equalities(m, c, j, g): + """ + Enforce a commitment deviation constraint on the aggregate of devices in a group. + + For commitment ``c`` at time index ``j``, this constraint couples the commitment + baseline (plus deviation variables) to the summed flow or stock of all devices + belonging to device group ``g``. StockCommitments aggregate device stocks, while + FlowCommitments aggregate device flows. Constraints are skipped if the commitment + is inactive at ``(c, j)`` or if the group contains no devices. + """ + if m.commitment_quantity[c, j] == -infinity: + return Constraint.Skip + + devices_in_group = device_group_lookup.get(c, {}).get(g, set()) + if not devices_in_group: + return Constraint.Skip + + center = ( + m.commitment_quantity[c, j] + + m.commitment_downwards_deviation[c] + + m.commitment_upwards_deviation[c] + ) + + if commitments[c]["class"].apply(lambda cl: cl == StockCommitment).all(): + center -= sum(_get_stock_change(m, d, j) for d in devices_in_group) + else: + center -= sum(m.ems_power[d, j] for d in devices_in_group) + + return ( + 0 if "upwards deviation price" in commitments[c].columns else None, + center, + 0 if "downwards deviation price" in commitments[c].columns else None, + ) + model.up_price = Param(model.c, initialize=price_up_select) model.down_price = Param(model.c, initialize=price_down_select) model.commitment_quantity = Param( @@ -565,6 +642,10 @@ def device_derivative_equalities(m, d, j): 0, ) + model.grouped_commitment_equalities = Constraint( + model.cjg, rule=grouped_commitment_equalities + ) + model.device_energy_bounds = Constraint(model.d, model.j, rule=device_bounds) model.device_power_bounds = Constraint( model.d, model.j, rule=device_derivative_bounds @@ -592,9 +673,7 @@ def device_derivative_equalities(m, d, j): model.ems_power_commitment_equalities = Constraint( model.cj, rule=ems_flow_commitment_equalities ) - model.device_energy_commitment_equalities = Constraint( - model.cj, model.d, rule=device_stock_commitment_equalities - ) + model.device_power_equalities = Constraint( model.d, model.j, rule=device_derivative_equalities ) diff --git a/flexmeasures/data/models/planning/tests/test_commitments.py b/flexmeasures/data/models/planning/tests/test_commitments.py index fa4d038ba6..195409afc5 100644 --- a/flexmeasures/data/models/planning/tests/test_commitments.py +++ b/flexmeasures/data/models/planning/tests/test_commitments.py @@ -1,6 +1,404 @@ +import pytest import pandas as pd +import numpy as np -from flexmeasures.data.models.planning import Commitment +from flexmeasures.data.models.planning import ( + Commitment, + StockCommitment, + FlowCommitment, +) +from flexmeasures.data.models.planning.utils import ( + initialize_index, +) +from flexmeasures.data.models.planning.linear_optimization import device_scheduler + + +def test_multi_feed_device_scheduler_shared_buffer(): + # ---- time setup + start = pd.Timestamp("2026-01-01T00:00+01") + end = pd.Timestamp("2026-01-02T00:00+01") + resolution = pd.Timedelta("PT1H") + index = initialize_index(start=start, end=end, resolution=resolution) + + # ---- three devices + devices = ["gas boiler", "heat pump power", "battery power"] + + # ---- device grouping + device_group = pd.Series( + { + 0: "shared thermal buffer", # gas boiler + 1: "shared thermal buffer", # "heat pump power" + 2: "battery SoC", # "battery power" + } + ) + device_commodity = pd.Series( + { + 0: "gas", # gas boiler + 1: "electricity", # "heat pump power" + 2: "electricity", # "battery power" + } + ) + equals = pd.Series(np.nan, index=index) + equals[-1] = 100 + device_constraints = [] + for d, device_name in enumerate(devices): + # 0 and 1 : derivative min 0 + # 2 : derivative min = - production capacity + + df = pd.DataFrame( + { + "min": 0, + "max": 100, + "equals": np.nan, + "derivative min": 0 if d in (0, 1) else -20, + "derivative max": 20, + "derivative equals": np.nan, + "derivative down efficiency": 0.9, + "derivative up efficiency": 0.9, + }, + index=index, + ) + device_constraints.append(df) + + ems_constraints = pd.DataFrame( + { + "derivative min": -40, + "derivative max": 40, + }, + index=index, + ) + + # ---- shared buffer max = 100 (soft) + max_soc = 100.0 + breach_price = 1_000.0 + min_soc = pd.Series(0, index=index) + min_soc[-1] = 100 + + # default commodity: electricity + # choice: electricity or gas + gas_price = pd.Series(300, index=index) + electricity_price = pd.Series(600, index=index, name="event_value") + electricity_price.iloc[12:14] = 200 + prices = {"gas": gas_price, "electricity": electricity_price} + + # Tie-breaking: prefer filling each device's storage as early as possible. + soc_max = 100.0 + penalty = 0.001 + + commitments = [] + + commitments.append( + StockCommitment( + name="buffer min", + index=index, + quantity=min_soc, + upwards_deviation_price=0, + downwards_deviation_price=-breach_price, + # instead of device=None, I considered to create a series for the devices that we need for this + # specific commitment. + device=pd.Series([[0, 1]] * len(index), index=index), + device_group=device_group, + ) + ) + for d, dev in enumerate(devices): + commitments.append( + StockCommitment( + name="buffer max", + index=index, + quantity=max_soc, + upwards_deviation_price=breach_price, + downwards_deviation_price=0, + device=pd.Series(d, index=index), + device_group=device_group, + ) + ) + commitments.append( + FlowCommitment( + name=device_commodity[d], + index=index, + quantity=0, + upwards_deviation_price=prices[device_commodity[d]], + downwards_deviation_price=prices[device_commodity[d]], + device=pd.Series(d, index=index), + device_group=device_commodity, + ) + ) + + commitments.append( + StockCommitment( + name=f"prefer a full storage {d} sooner", + index=index, + quantity=soc_max, + upwards_deviation_price=0, + downwards_deviation_price=-penalty, + device=pd.Series(d, index=index), + ) + ) + + # ---- run scheduler + planned_power, planned_costs, results, model = device_scheduler( + device_constraints=device_constraints, + ems_constraints=ems_constraints, + commitments=commitments, + initial_stock=0, + ) + + # ---- sanity: model solved optimally + assert results.solver.termination_condition == "optimal" + + # ---- key assertion: exactly TWO commitment groups + # - one for "shared thermal buffer" + # - one for "battery SoC" + # + # i.e. NOT three (which would indicate per-device baselines) + commitment_groups = set(commitments[0].device_group.values) + commodity_commitments = { + c.name + for c in commitments + if isinstance(c, FlowCommitment) and c.name in {"gas", "electricity"} + } + assert commodity_commitments == {"gas", "electricity"} + + # Sum per-commitment costs grouped by name so that duplicate names + # (e.g. "electricity" for both heat pump and battery) are accumulated. + electricity_cost = sum( + costs + for c, costs in zip(commitments, model.commitment_costs.values()) + if c.name == "electricity" + ) + gas_cost = sum( + costs + for c, costs in zip(commitments, model.commitment_costs.values()) + if c.name == "gas" + ) + commodity_costs = {"gas": gas_cost, "electricity": electricity_cost} + + assert set(commodity_costs.keys()) == {"gas", "electricity"} + + assert commitment_groups == {"shared thermal buffer"} + + # The shared buffer minimum (SoC ≥ 100 at the final step) must be met without any breach + buffer_min_cost = sum( + costs + for c, costs in zip(commitments, model.commitment_costs.values()) + if c.name == "buffer min" + ) + assert buffer_min_cost == 0, ( + f"Shared buffer target was breached (breach cost = {buffer_min_cost}). " + "This may indicate that baseline costs were incorrectly duplicated." + ) + + # At hours 12–13 electricity (200) is cheaper than gas (300), so the heat pump + # runs at full power: 2 h × 20 kW = 40 kWh flow → 40 × 0.9 = 36 kWh SoC + # contribution to the shared buffer. The gas boiler covers the remainder: + # (100 − 36) kWh SoC / 0.9 efficiency × 300 price. + assert gas_cost == pytest.approx((100 - 40 * 0.9) / 0.9 * 300, rel=1e-6) + + # Expect the total electricity costs to be: + # 2 * 20 for the heat pump + # 2 * 20 for the battery charging + # minus 2 * 20 * 0.9 * 0.9 for the battery discharging after roundtrip efficiency + assert electricity_cost == pytest.approx( + 200 * (40 + 40) - 600 * (40 * 0.9 * 0.9), rel=1e-6 + ) + + +def _run_hp_buffer_scenario(index, target_soc, shared: bool): + """ + Helper: run the two-heat-pump scheduler with either a shared or per-device buffer + commitment and return costs for assertions. + + Each heat pump (device 0 and 1) can supply at most + derivative_max × hours × efficiency = 20 kW × 24 h × 0.9 = 432 kWh. + + shared=True → one StockCommitment on the *combined* stock of devices 0+1 + (maximum reachable: 864 kWh). + shared=False → two separate StockCommitments, one per HP, each requiring + target_soc on its *own* stock (maximum per device: 432 kWh). + """ + n = len(index) + breach_price = 1_000.0 + energy_price = pd.Series(100, index=index) + + device_group = pd.Series( + { + 0: "shared thermal buffer", + 1: "shared thermal buffer", + 2: "battery SoC", + } + ) + + # ---- device constraints + # device 0: heat pump A (charge only) + # device 1: heat pump B (charge only) + # device 2: battery (charge and discharge) + device_constraints = [] + for d in range(3): + df = pd.DataFrame( + { + "min": 0, + "max": 500, + "equals": np.nan, + "derivative min": 0 if d in (0, 1) else -20, + "derivative max": 20, + "derivative equals": np.nan, + "derivative down efficiency": 0.9, + "derivative up efficiency": 0.9, + }, + index=index, + ) + device_constraints.append(df) + + ems_constraints = pd.DataFrame( + {"derivative min": -40, "derivative max": 40}, + index=index, + ) + + min_soc = pd.Series(0.0, index=index) + min_soc.iloc[-1] = target_soc + + # Tie-breaking: prefer filling each device's storage as early as possible. + soc_max = 500.0 # matches device_constraints["max"] + penalty = 0.001 + + commitments = [] + + if shared: + # One commitment covering the *combined* stock of both HPs. + commitments.append( + StockCommitment( + name="buffer min", + index=index, + quantity=min_soc, + upwards_deviation_price=0, + downwards_deviation_price=-breach_price, + device=pd.Series([[0, 1]] * n, index=index), + device_group=device_group, + ) + ) + else: + # Two separate commitments, one per HP — each must reach target_soc alone. + for d in range(2): + commitments.append( + StockCommitment( + name="buffer min", + index=index, + quantity=min_soc, + upwards_deviation_price=0, + downwards_deviation_price=-breach_price, + device=pd.Series(d, index=index), + device_group=device_group, + ) + ) + + # Individual upper bounds (soft) and energy price for all three devices. + for d in range(3): + commitments.append( + StockCommitment( + name=f"buffer max {d}", + index=index, + quantity=500.0, + upwards_deviation_price=breach_price, + downwards_deviation_price=0, + device=pd.Series(d, index=index), + device_group=device_group, + ) + ) + commitments.append( + FlowCommitment( + name="energy", + index=index, + quantity=0, + upwards_deviation_price=energy_price, + downwards_deviation_price=energy_price, + device=pd.Series(d, index=index), + device_group=device_group, + ) + ) + commitments.append( + StockCommitment( + name=f"prefer a full storage {d} sooner", + index=index, + quantity=soc_max, + upwards_deviation_price=0, + downwards_deviation_price=-penalty, + device=d, + ) + ) + + planned_power, planned_costs, results, model = device_scheduler( + device_constraints=device_constraints, + ems_constraints=ems_constraints, + commitments=commitments, + initial_stock=0, + ) + + assert results.solver.termination_condition == "optimal" + + buffer_min_cost = sum( + v + for c, v in zip(commitments, model.commitment_costs.values()) + if c.name == "buffer min" + ) + energy_cost = sum( + v + for c, v in zip(commitments, model.commitment_costs.values()) + if c.name == "energy" + ) + return {"buffer_min_cost": buffer_min_cost, "energy_cost": energy_cost} + + +def test_device_group_shared_buffer(): + """ + Two heat pumps (devices 0 and 1) charge a shared thermal buffer with a target of + 800 kWh by the last time slot. + + Each HP can supply at most 20 kW × 24 h × 0.9 = 432 kWh on its own, so neither + can reach 800 kWh individually. Together they can supply up to 864 kWh, so the + target is feasible when the commitment tracks their *combined* stock. + + This test verifies two contrasting scenarios: + + 1. Shared buffer (device_group): one StockCommitment on the combined stock of both + HPs. The optimizer fills 800 kWh across the two devices with zero breach cost. + + 2. Separate buffers (no device_group): one StockCommitment per HP, each requiring + 800 kWh on its own stock. Each HP falls short by ~368 kWh, so both commitments + incur a breach, and the total breach cost is positive. + """ + start = pd.Timestamp("2026-01-01T00:00+01") + end = pd.Timestamp("2026-01-02T00:00+01") + resolution = pd.Timedelta("PT1H") + index = initialize_index(start=start, end=end, resolution=resolution) + + # 800 kWh: above what one HP can reach (432 kWh), below what two can reach (864 kWh). + target_soc = 800.0 + + shared_result = _run_hp_buffer_scenario(index, target_soc, shared=True) + separate_result = _run_hp_buffer_scenario(index, target_soc, shared=False) + + shared_cost = shared_result["buffer_min_cost"] + separate_cost = separate_result["buffer_min_cost"] + shared_energy = shared_result["energy_cost"] + separate_energy = separate_result["energy_cost"] + + assert shared_cost == 0, ( + f"Shared buffer: both HPs together can reach {target_soc} kWh, " + f"so breach cost must be zero (got {shared_cost})" + ) + assert separate_cost > 0, ( + f"Separate buffers: each HP alone cannot reach {target_soc} kWh, " + f"so breach cost must be positive (got {separate_cost})" + ) + + # With shared buffer, the optimizer charges exactly 800 kWh combined. + # Total energy flow = 800 / 0.9 (accounting for charge efficiency). + assert shared_energy == pytest.approx(target_soc / 0.9 * 100, rel=1e-6) + + # With separate buffers, each HP charges at maximum power for all 24 hours + # since the 800 kWh individual target is unreachable. + # Total energy = 2 HPs × 20 kW × 24 h × 100 price. + assert separate_energy == pytest.approx(2 * 20 * 24 * 100, rel=1e-6) def make_index(n: int = 5) -> pd.DatetimeIndex: diff --git a/flexmeasures/data/tests/test_scheduling_simultaneous.py b/flexmeasures/data/tests/test_scheduling_simultaneous.py index c9a1fd8836..d58d6ab4c8 100644 --- a/flexmeasures/data/tests/test_scheduling_simultaneous.py +++ b/flexmeasures/data/tests/test_scheduling_simultaneous.py @@ -107,9 +107,9 @@ def test_create_simultaneous_jobs( total_cost = ev_costs + battery_costs # Define expected costs based on resolution - expected_ev_costs = 2.2375 - expected_battery_costs = -5.515 expected_total_cost = -3.2775 + expected_ev_costs = 2.3125 + expected_battery_costs = expected_total_cost - expected_ev_costs # Check costs assert (