Share storage-efficiency per stock group and index solver stock by group (#2325)

* Speed up device_scheduler with a recursive stock balance

`_get_stock_change` expressed a device's stock at step j as a running sum
over all earlier steps, so the number of model nonzeros grew quadratically
with the scheduling horizon (e.g. ~12k nonzeros at 12h, ~349k at 3 days).
That dominated solve time on longer horizons.

Introduce an explicit `device_stock` variable coupled by a recursive
`device_stock_balance` constraint (stock[j] = a*stock[j-1] + b*change[j],
with a, b the per-step loss coefficients that apply_stock_changes_and_losses
computes). This is mathematically equivalent -- verified identical
objectives across horizons and against the planning test suite -- and keeps
the number of nonzeros linear in the horizon. Solve time drops ~10x at the
2-day default horizon and ~23x at 3 days.

Found while investigating SeitaBV/ems#172.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LeqGFrHfGrBHAJyjdAyr3y

* Reference PR #2282 in the changelog entry

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01LeqGFrHfGrBHAJyjdAyr3y

* Address review: apply changelog suggestion and cache stock-group device lists

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Fix degenerate optimum in alignment regression test after merging main

The recursion formulation lands on a cost-equal alternative vertex (A charges 3,
B discharges 1 for free under net site metering). Make device B lossy in both
directions so the optimum is unique again.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Share storage-efficiency per stock group and index solver stock by group (#2324)

- storage-efficiency may now be defined on the entry holding a shared stock's
  SoC parameters (or on exactly one member device) and applies to all members;
  conflicting definitions fail fast, mirroring how #2321 treats the other SoC
  parameters.
- device_scheduler now models one stock variable and one balance recursion per
  stock group instead of per device, dropping the redundant secondary
  recursions, and validates that grouped devices share their storage
  efficiency and initial stock.
- Fix a latent key collision in device_to_group: an ungrouped device whose
  index equalled a state-of-charge sensor id was silently merged into that
  stock group, feeding its flows into the wrong stock.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Add changelog entries for PR #2325

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Enable storage losses in the dynamic-capacity test and consolidate changelog

The buffer's storage-efficiency (previously ignored, with a 'does not work yet'
todo) is now honored, so enable it at 99% per 15 minutes and update the
expected schedules: the heater covers the loss-induced remainder beyond the
maxed-out boiler around the clock, and the boiler tops up the buffer in the
final hour of the cheap-electricity window.

Also fold the two changelog entries into the existing shared-storage feature
entry for v1.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Attach soft SoC commitments to their stock, not to a device index

Enabler for scoping #2194's SoC-constraint relaxation per stock:

- StockCommitment gains a 'stock' attribute holding the stock key; the
  StorageScheduler stamps it on all soft SoC commitments (soc-minima,
  soc-maxima, prefer-full).
- The solver couples stock-scoped commitments to the stock group as a whole
  (via the group's first device, whose stock is the group's stock), regardless
  of which device index the commitment names.
- DeviceInventory.stock_constraint_device(stock_key) exposes the device that
  applies a stock's SoC constraints, complementing stock_groups for
  deserialization-time scoping logic.
- The solver's grouped-device efficiency guard now tolerates a missing
  efficiency column (the default, lossless case), which direct callers of
  device_scheduler commonly omit.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

* Remove dead device_stock_commitment_equalities

It was never registered as a model constraint; all device-level commitments
route through grouped_commitment_equalities.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

---------

Co-authored-by: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

Author: Flix6x

documentation/changelog.rst

@@ -20,7 +20,7 @@ New features
 * Floor off-clock API datetimes to a non-instantaneous sensor's resolution by default when ingesting sensor data, uploading sensor data, and handling scheduler flex-model timed events; configurable with the ``floor_datetimes_to_resolution`` sensor attribute [see `PR #2146 <https://www.github.com/FlexMeasures/flexmeasures/pull/2146>`_]
 * Sensor references in flex-model and flex-context support various ways of filtering by source [see `PR #2209 <https://www.github.com/FlexMeasures/flexmeasures/pull/2209>`_]
 * Let storage scheduling infer missing ``power-capacity`` from directional device capacities before falling back to site capacity, and default the missing opposite capacity to zero when only a non-zero ``consumption-capacity`` or ``production-capacity`` is configured [see `PR #2222 <https://www.github.com/FlexMeasures/flexmeasures/pull/2222>`_]
-* Support multiple feeders to a shared storage [see `PR #2001 <https://www.github.com/FlexMeasures/flexmeasures/pull/2001>`_, `PR #2321 <https://www.github.com/FlexMeasures/flexmeasures/pull/2321>`_ and `PR #2322 <https://www.github.com/FlexMeasures/flexmeasures/pull/2322>`_]
+* Support multiple feeders to a shared storage [see `PR #2001 <https://www.github.com/FlexMeasures/flexmeasures/pull/2001>`_, `PR #2321 <https://www.github.com/FlexMeasures/flexmeasures/pull/2321>`_, `PR #2322 <https://www.github.com/FlexMeasures/flexmeasures/pull/2322>`_ and `PR #2325 <https://www.github.com/FlexMeasures/flexmeasures/pull/2325>`_]
 * The flex-context can now define multiple commodities, each specifying their own prices and grid capacities [see `PR #1946 <https://www.github.com/FlexMeasures/flexmeasures/pull/1946>`_, `PR #2172 <https://www.github.com/FlexMeasures/flexmeasures/pull/2172>`_, `PR #2235 <https://www.github.com/FlexMeasures/flexmeasures/pull/2235>`_ and `PR #2271 <https://www.github.com/FlexMeasures/flexmeasures/pull/2271>`_]
 * 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 <https://www.github.com/FlexMeasures/flexmeasures/pull/2272>`_]
 * CLI support for adding/editing account attributes [see `PR #2242 <https://www.github.com/FlexMeasures/flexmeasures/pull/2242>`_]

flexmeasures/data/models/planning/__init__.py

@@ -352,6 +352,10 @@ class Commitment:
     upwards_deviation_price: pd.Series = 0
     downwards_deviation_price: pd.Series = 0
     commodity: str | pd.Series | None = None
+    #: Stock key of the stock the commitment pertains to (StockCommitments only).
+    #: When set, the solver couples the commitment to the stock group as a whole,
+    #: rather than to the device index named by ``device``.
+    stock: int | None = None
 
     def __post_init__(self):
         # device_group is a device→label lookup table, not a time series;
@@ -542,6 +546,9 @@ def to_frame(self) -> pd.DataFrame:
             # scalar commodity
             df["commodity"] = self.commodity
 
+        # stock key (scalar; set on stock-scoped StockCommitments)
+        df["stock"] = self.stock
+
         return df
 
 

flexmeasures/data/models/planning/devices.py

@@ -344,6 +344,16 @@ def stock_params(self, stock_key: int) -> dict | None:
         """Return the flex-model entry holding the SoC parameters of the given stock."""
         return self.stock_entries.get(stock_key)
 
+    def stock_constraint_device(self, stock_key: int) -> int | None:
+        """Return the device that applies the given stock's SoC constraints.
+
+        Hard SoC constraints (and the device-indexed remains of softened ones) land on
+        the first device of the stock group, whose stock represents the group's stock
+        in the solver. Use :attr:`stock_groups` to look up all member devices instead.
+        """
+        group_devices = self.stock_groups.get(stock_key)
+        return group_devices[0] if group_devices else None
+
     @cached_property
     def commodity_to_devices(self) -> dict[str, list[int]]:
         """Map each commodity to its device indices, in canonical solver order."""

flexmeasures/data/models/planning/linear_optimization.py

@@ -133,23 +133,52 @@ def device_scheduler(  # noqa C901
     # and map stock group -> all member devices (used for stock accumulation).
     device_to_group = {}
 
+    # 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.
     if stock_groups:
         for g, devices in stock_groups.items():
             for d in devices:
-                device_to_group[d] = g
-        # For devices not in any stock group (e.g., inflexible devices),
-        # map them to themselves so they're treated as individual groups
-        for d in range(len(device_constraints)):
-            if d not in device_to_group:
-                device_to_group[d] = d
-    else:
-        for d in range(len(device_constraints)):
-            device_to_group[d] = d
+                device_to_group[d] = f"stock:{g}"
+    # Devices not in any stock group (e.g. inflexible devices) form individual groups.
+    for d in range(len(device_constraints)):
+        if d not in device_to_group:
+            device_to_group[d] = f"device:{d}"
 
     group_to_devices: dict[int, list[int]] = {}
     for d, g in device_to_group.items():
         group_to_devices.setdefault(g, []).append(d)
 
+    # 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.
+    for g, group_devices in group_to_devices.items():
+        if len(group_devices) > 1:
+            # A missing efficiency column means the default (no losses) applies.
+            group_efficiency = device_constraints[group_devices[0]].get("efficiency")
+            for d in group_devices[1:]:
+                efficiency = device_constraints[d].get("efficiency")
+                if (
+                    (efficiency is None) != (group_efficiency is None)
+                    or efficiency is not None
+                    and not efficiency.equals(group_efficiency)
+                ):
+                    raise ValueError(
+                        f"Devices {group_devices} share stock group {g} but have different"
+                        " storage efficiencies. The storage efficiency is a property of the"
+                        " shared stock, so define it once per stock group."
+                    )
+            if isinstance(initial_stock, list):
+                group_initial_stocks = {
+                    initial_stock[d] if d < len(initial_stock) else 0
+                    for d in group_devices
+                }
+                if len(group_initial_stocks) > 1:
+                    raise ValueError(
+                        f"Devices {group_devices} share stock group {g} but have different"
+                        " initial stocks. The initial stock is a property of the shared"
+                        " stock, so define it once per stock group."
+                    )
+
     # Move commitments from old structure to new
     if commitments is None:
         commitments = []
@@ -274,6 +303,18 @@ def convert_commitments_to_subcommitments(
     device_group_lookup = {}
 
     for c, df in enumerate(commitments):
+        # Stock-scoped commitments couple to their stock group as a whole, regardless
+        # of which device index they name: the group's first device carries the group's
+        # stock, so a single-member group suffices (also avoiding double-counting the
+        # shared stock when the commitment names multiple members).
+        if "stock" in df.columns and pd.notna(df["stock"].iloc[0]):
+            stock_group_key = f"stock:{int(df['stock'].iloc[0])}"
+            if stock_group_key in group_to_devices:
+                device_group_lookup[c] = {
+                    stock_group_key: {group_to_devices[stock_group_key][0]}
+                }
+                continue
+
         if "device" not in df.columns:
             # EMS-level commitment: no device grouping needed here;
             # handled by ems_flow_commitment_equalities.
@@ -541,10 +582,13 @@ def grouped_commitment_equalities(m, c, j, g):
     )
     model.device_power_up = Var(model.d, model.j, domain=NonNegativeReals, initialize=0)
     model.device_power_sign = Var(model.d, model.j, domain=Binary, initialize=0)
-    # Stock per device per time step, coupled recursively by device_stock_balance.
+    # Stock per stock group per time step, coupled recursively by group_stock_balance.
     # Having it as a variable (rather than a running sum expression) keeps the number
-    # of model nonzeros linear, rather than quadratic, in the scheduling horizon.
-    model.device_stock = Var(model.d, model.j, domain=Reals, initialize=0)
+    # of model nonzeros linear, rather than quadratic, in the scheduling horizon, and
+    # indexing it by stock group (rather than by device) avoids duplicating the
+    # recursion for each device sharing a stock.
+    model.sg = Set(initialize=sorted(group_to_devices), doc="Set of stock groups")
+    model.group_stock = Var(model.sg, model.j, domain=Reals, initialize=0)
     model.commitment_downwards_deviation = Var(
         model.c,
         domain=NonPositiveReals,
@@ -565,14 +609,13 @@ def _initial_stock_of(d):
             return initial_stock[d] if d < len(initial_stock) else 0
         return initial_stock
 
-    def _stock_change_at(m, d, j):
-        """Stock change of device d's stock group during time step j (before losses)."""
-        devices = group_to_devices[device_to_group[d]]
+    def _stock_change_at(m, g, j):
+        """Stock change of stock group g during time step j (before losses)."""
         return sum(
             m.device_power_down[dev, j] / m.device_derivative_down_efficiency[dev, j]
             + m.device_power_up[dev, j] * m.device_derivative_up_efficiency[dev, j]
             + m.stock_delta[dev, j]
-            for dev in devices
+            for dev in group_to_devices[g]
         )
 
     def _loss_coefficients(efficiency: float) -> tuple[float, float]:
@@ -588,20 +631,24 @@ def _loss_coefficients(efficiency: float) -> tuple[float, float]:
             return 1.0, 1.0
         return efficiency, (efficiency - 1) / math.log(efficiency)
 
-    def device_stock_balance(m, d, j):
-        """Recursively couple a device's stock to the previous step's stock.
+    def group_stock_balance(m, g, j):
+        """Recursively couple a stock group's stock to the previous step's stock.
 
         Expressing stock[j] as a running sum over all k <= j (as this once did) makes
         the number of nonzeros grow quadratically with the scheduling horizon. The
         recursion below is equivalent and keeps it linear.
+
+        The group's devices share their storage efficiency and initial stock
+        (validated above), so the first device can represent the group here.
         """
-        a, b = _loss_coefficients(m.device_efficiency[d, j])
-        previous = m.device_stock[d, j - 1] if j > 0 else _initial_stock_of(d)
-        return m.device_stock[d, j] == a * previous + b * _stock_change_at(m, d, j)
+        d0 = group_to_devices[g][0]
+        a, b = _loss_coefficients(m.device_efficiency[d0, j])
+        previous = m.group_stock[g, j - 1] if j > 0 else _initial_stock_of(d0)
+        return m.group_stock[g, j] == a * previous + b * _stock_change_at(m, g, j)
 
     def _get_stock_change(m, d, j):
         """Stock change of the stock group of device d, from the start until time j."""
-        return m.device_stock[d, j] - _initial_stock_of(d)
+        return m.group_stock[device_to_group[d], j] - _initial_stock_of(d)
 
     # Add constraints as a tuple of (lower bound, value, upper bound)
     def device_bounds(m, d, j):
@@ -661,34 +708,6 @@ def commitment_down_derivative_sign(m, c):
         """Down deviation active only if sign points down."""
         return -m.commitment_downwards_deviation[c] <= Mc * (1 - m.commitment_sign[c])
 
-    def device_stock_commitment_equalities(m, c, j, d):
-        """Couple device stocks to each commitment."""
-        if (
-            "device" not in commitments[c].columns
-            or (commitments[c]["device"] != d).all()
-            or m.commitment_quantity[c, j] == -infinity
-        ):
-            # Commitment c does not concern device d
-            return Constraint.Skip
-
-        # Determine center part of the lhs <= center part <= rhs constraint
-        center_part = (
-            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_part -= _get_stock_change(m, d, j)
-        elif commitments[c]["class"].apply(lambda cl: cl == FlowCommitment).all():
-            center_part -= m.ems_power[d, j]
-        else:
-            raise NotImplementedError("Unknown commitment class")
-        return (
-            0 if "upwards deviation price" in commitments[c].columns else None,
-            center_part,
-            0 if "downwards deviation price" in commitments[c].columns else None,
-        )
-
     def ems_flow_commitment_equalities(m, c, j):
         """Couple EMS flow commitments to device flows, optionally filtered by commodity."""
 
@@ -728,7 +747,7 @@ def device_derivative_equalities(m, d, j):
         model.cjg, rule=grouped_commitment_equalities
     )
 
-    model.device_stock_balance = Constraint(model.d, model.j, rule=device_stock_balance)
+    model.group_stock_balance = Constraint(model.sg, model.j, rule=group_stock_balance)
     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

flexmeasures/data/models/planning/storage.py

@@ -148,6 +148,13 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
         self.stock_models = inventory.stock_entries
         # The stock groups' device indices align with the device models
         self.stock_groups = inventory.stock_groups
+        # Soft SoC constraints are attached to their stock (see StockCommitment.stock),
+        # so the solver couples them to the stock group rather than to a device index.
+        device_stock_key = {
+            d: stock_key
+            for stock_key, group_devices in self.stock_groups.items()
+            for d in group_devices
+        }
 
         # List the asset(s) and sensor(s) being scheduled
         sensors: list[Sensor | None] = inventory.power_sensors
@@ -214,10 +221,37 @@ def _prepare(self, skip_validation: bool = False) -> tuple:  # noqa: C901
             prefer_charging_sooner[d0] = stock_model.get("prefer_charging_sooner")
             prefer_curtailing_later[d0] = stock_model.get("prefer_curtailing_later")
 
-        # todo: move storage-efficiency into a shared parameter for the first device belonging to a shared storage
         storage_efficiency = [
             flex_model_d.get("storage_efficiency") for flex_model_d in flex_model
         ]
+        # The storage efficiency is a property of the stock, not of a connected device:
+        # for shared stocks, it may be defined on the entry holding the stock's SoC
+        # parameters or on a single member device, and applies to all members.
+        for stock_id, stock_devices in self.stock_groups.items():
+            if len(stock_devices) <= 1:
+                continue
+            definitions = []
+            stock_model = self.stock_models.get(stock_id)
+            if (
+                stock_model is not None
+                and stock_model.get("storage_efficiency") is not None
+            ):
+                definitions.append(stock_model["storage_efficiency"])
+            definitions.extend(
+                storage_efficiency[d]
+                for d in stock_devices
+                if storage_efficiency[d] is not None
+            )
+            if len(set(map(id, definitions))) > 1:
+                raise ValueError(
+                    f"Multiple flex-model entries define a storage-efficiency for the same"
+                    f" stock (state-of-charge sensor {stock_id}). The storage efficiency"
+                    f" is a property of the shared stock, so please define it on a single"
+                    f" entry."
+                )
+            shared_efficiency = definitions[0] if definitions else None
+            for d in stock_devices:
+                storage_efficiency[d] = shared_efficiency
         consumption = [flex_model_d.get("consumption") for flex_model_d in flex_model]
         production = [flex_model_d.get("production") for flex_model_d in flex_model]
         consumption_capacity = [
@@ -650,6 +684,7 @@ def device_list_series(
                     downwards_deviation_price=-penalty,
                     index=index,
                     device=d,
+                    stock=device_stock_key.get(d),
                 )
                 commitments.append(commitment)
 
@@ -745,6 +780,7 @@ def device_list_series(
                     index=index,
                     _type="any",
                     device=d,
+                    stock=device_stock_key.get(d),
                 )
                 commitments.append(commitment)
 
@@ -755,6 +791,7 @@ def device_list_series(
                     downwards_deviation_price=-all_soc_minima_breach_price,
                     index=index,
                     device=d,
+                    stock=device_stock_key.get(d),
                 )
                 commitments.append(commitment)
 
@@ -819,6 +856,7 @@ def device_list_series(
                     index=index,
                     _type="any",
                     device=d,
+                    stock=device_stock_key.get(d),
                 )
                 commitments.append(commitment)
 
@@ -829,6 +867,7 @@ def device_list_series(
                     upwards_deviation_price=all_soc_maxima_breach_price,
                     index=index,
                     device=d,
+                    stock=device_stock_key.get(d),
                 )
                 commitments.append(commitment)