Cover the committed quantity in Mc, and tighten it

documentation/changelog.rst

@@ -101,6 +101,7 @@ Infrastructure / Support
 
 Bugfixes
 -----------
+* The big-M bounding commitment deviations now accounts for the committed quantity, so a committed quantity far beyond the devices' flow limits no longer renders the problem infeasible when a non-convex cost curve activates the commitment-sign constraints; the bound is also tighter (covering one time step rather than the whole horizon, for flow commitments), which improves MIP numerics [see `PR #2411 <https://www.github.com/FlexMeasures/flexmeasures/pull/2411>`_]
 * In a multi-device flex-model, a device without a stock (e.g. a converter port or curtailable generator) silently disabled constraint validation for all devices after it; validation now covers every device, and also newly checks that each device's power bounds do not contradict each other, so a contradictory hard bound fails with a clear per-time-step message instead of a bare solver infeasibility [see `PR #2252 <https://www.github.com/FlexMeasures/flexmeasures/pull/2252>`_]
 * The scheduler now rejects a commitment that no constraint would bind — a stock commitment naming no device or known stock group, or a commodity commitment for a commodity that no commitment maps devices to — instead of silently dropping it from the problem, or letting a favourably priced deviation make the problem unbounded; the error names the commitment [see `PR #2410 <https://www.github.com/FlexMeasures/flexmeasures/pull/2410>`_ and `PR #2413 <https://www.github.com/FlexMeasures/flexmeasures/pull/2413>`_]
 * Show icons for more asset types in the UI's asset structure view, which previously fell back to a question mark: the ``wind``, ``process`` and ``heat-storage`` types that FlexMeasures seeds by default, and EV infrastructure under its various names (such as ``one-way_evse``, ``two-way_evse``, ``evse``, ``charging_station`` and ``charging_hub``) and building services equipment (``hvac``, ``ahu``, ``dhw``, ``heatpump``, ``chiller``, ``lighting`` and ``other-loads``). Asset type names are now matched ignoring case and separators, so an asset type named ``charge-point`` gets the same icon as ``chargepoint`` [see `PR #2391 <https://www.github.com/FlexMeasures/flexmeasures/pull/2391>`_]

flexmeasures/data/models/planning/scheduling_problem.py

@@ -24,7 +24,11 @@
 from flask import current_app
 from pandas.tseries.frequencies import to_offset
 
-from flexmeasures.data.models.planning import Commitment, FlowCommitment
+from flexmeasures.data.models.planning import (
+    Commitment,
+    FlowCommitment,
+    StockCommitment,
+)
 from flexmeasures.data.models.planning.utils import initialize_df, initialize_series
 
 infinity = float("inf")
@@ -503,9 +507,61 @@ def prepare_scheduling_problem(  # noqa C901
     bigM_columns = ["derivative max", "derivative min", "derivative equals"]
     # Compute a good value for our Big-Ms
     # Md is used to constrain the search space for device power
-    # Mc is used to constrain the search space for commitment deviations
     Md = np.nanmax([np.nanmax(d[bigM_columns].abs()) for d in device_constraints])
-    Mc = np.nansum([np.nansum(d[bigM_columns].abs()) for d in device_constraints])
+
+    # Mc is used to constrain the search space for commitment deviations.
+    # A deviation makes up the gap between a committed quantity and an aggregated device flow (or stock change),
+    # so it is bounded by the largest absolute committed quantity plus the devices' summed flow limits:
+    # summed at any one time step for flow commitments,
+    # and summed over the whole horizon for stock commitments (a stock change accumulates flows since the start).
+    # A device without flow limits at some time step contributes nothing there (as it did before this bound was tightened),
+    # its power being unbounded anyway.
+    per_device_step_limits = [
+        d[bigM_columns].astype(float).abs().max(axis=1).fillna(0).to_numpy()
+        for d in device_constraints
+    ]
+    per_step_total = sum(per_device_step_limits)
+    Mc = float(np.max(per_step_total))
+    has_stock_commitment = any(
+        "class" in c.columns and (c["class"] == StockCommitment).any()
+        for c in commitments
+    )
+    if has_stock_commitment:
+        # A stock change is not a raw flow:
+        # it passes through the derivative efficiencies (P_up * eta_up + P_down / eta_down) and includes the explicit stock delta,
+        # so each device's flow limits are scaled by its worst-case conversion gain, and its stock deltas are added,
+        # before summing over the horizon.
+        horizon_stock_change_limit = 0.0
+        for d, flow_limits in zip(device_constraints, per_device_step_limits):
+            gain = np.ones(len(flow_limits))
+            if "derivative up efficiency" in d.columns:
+                gain = np.maximum(
+                    gain,
+                    d["derivative up efficiency"].astype(float).fillna(1).to_numpy(),
+                )
+            if "derivative down efficiency" in d.columns:
+                gain = np.maximum(
+                    gain,
+                    1
+                    / d["derivative down efficiency"]
+                    .astype(float)
+                    .fillna(1)
+                    .to_numpy(),
+                )
+            deltas = (
+                d["stock delta"].astype(float).fillna(0).abs().to_numpy()
+                if "stock delta" in d.columns
+                else 0.0
+            )
+            horizon_stock_change_limit += float(np.sum(flow_limits * gain + deltas))
+        Mc = max(Mc, horizon_stock_change_limit)
+    if commitments:
+        quantities = np.abs(
+            np.concatenate([c["quantity"].to_numpy(dtype=float) for c in commitments])
+        )
+        finite_quantities = quantities[np.isfinite(quantities)]
+        if len(finite_quantities):
+            Mc += float(np.max(finite_quantities))
 
     # Both Md and Mc have to be 1 MW, at least
     Md = max(Md, 1)

flexmeasures/data/models/planning/tests/test_big_m.py

@@ -0,0 +1,144 @@
+"""Tests for the big-M values bounding the scheduler's search space."""
+
+from __future__ import annotations
+
+from datetime import timedelta
+
+import numpy as np
+import pandas as pd
+import pytest
+
+from flexmeasures.data.models.planning import FlowCommitment, StockCommitment
+from flexmeasures.data.models.planning.linear_optimization import device_scheduler
+from flexmeasures.data.models.planning.scheduling_problem import (
+    prepare_scheduling_problem,
+)
+from flexmeasures.data.models.planning.utils import initialize_df
+
+#: Run every test in this module under both scheduler backends (see conftest).
+RUN_UNDER_EACH_SOLVER = True
+
+COLUMNS = [
+    "equals",
+    "max",
+    "min",
+    "efficiency",
+    "derivative equals",
+    "derivative max",
+    "derivative min",
+    "derivative down efficiency",
+    "derivative up efficiency",
+    "stock delta",
+]
+
+START = pd.Timestamp("2020-01-01T00:00:00")
+END = pd.Timestamp("2020-01-01T04:00:00")
+RESOLUTION = timedelta(hours=1)
+
+
+def make_device_constraints(power_capacity: float) -> pd.DataFrame:
+    device_constraints = initialize_df(COLUMNS, START, END, RESOLUTION)
+    device_constraints["derivative max"] = power_capacity
+    device_constraints["derivative min"] = -power_capacity
+    return device_constraints
+
+
+def make_index() -> pd.DatetimeIndex:
+    return initialize_df(COLUMNS, START, END, RESOLUTION).index
+
+
+def test_Mc_covers_flow_limits_per_step_plus_the_committed_quantity():
+    """For flow commitments, a deviation spans at most the committed quantity plus one time step's summed flow limits."""
+    commitment = FlowCommitment(
+        name="energy",
+        quantity=-100,
+        upwards_deviation_price=1,
+        downwards_deviation_price=-1,
+        index=make_index(),
+    )
+    problem = prepare_scheduling_problem(
+        device_constraints=[make_device_constraints(0.5), make_device_constraints(2)],
+        ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION),
+        commitments=[commitment],
+    )
+    assert problem.Md == 2
+    assert problem.Mc == 0.5 + 2 + 100
+
+
+def test_Mc_covers_the_horizon_for_stock_commitments():
+    """A stock commitment's deviation accumulates flows since the start, so Mc must cover the whole horizon."""
+    index = make_index()
+    commitment = StockCommitment(
+        name="soc",
+        quantity=0.5,
+        upwards_deviation_price=1,
+        downwards_deviation_price=-1,
+        device=pd.Series(0, index=index),
+        index=index,
+    )
+    problem = prepare_scheduling_problem(
+        device_constraints=[make_device_constraints(0.5), make_device_constraints(2)],
+        ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION),
+        commitments=[commitment],
+    )
+    # 4 time steps of 0.5 + 2 flow limits each, plus the committed quantity.
+    assert problem.Mc == 4 * 2.5 + 0.5
+
+
+def test_Mc_scales_stock_changes_by_conversion_efficiencies_and_stock_deltas():
+    """A stock change passes through the derivative efficiencies and includes the stock delta, unlike a raw flow."""
+    index = make_index()
+    commitment = StockCommitment(
+        name="soc",
+        quantity=0.5,
+        upwards_deviation_price=1,
+        downwards_deviation_price=-1,
+        device=pd.Series(0, index=index),
+        index=index,
+    )
+    device_constraints = make_device_constraints(0.5)
+    device_constraints["derivative up efficiency"] = 2
+    device_constraints["derivative down efficiency"] = 0.5
+    device_constraints["stock delta"] = 0.25
+    problem = prepare_scheduling_problem(
+        device_constraints=[device_constraints],
+        ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION),
+        commitments=[commitment],
+    )
+    # 4 time steps of 0.5 flow scaled by the worst-case conversion gain max(2, 1/0.5) plus a 0.25 stock delta,
+    # plus the committed quantity.
+    assert problem.Mc == 4 * (0.5 * 2 + 0.25) + 0.5
+
+
+def test_Mc_is_at_least_one():
+    problem = prepare_scheduling_problem(
+        device_constraints=[make_device_constraints(0.001)],
+        ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION),
+    )
+    assert problem.Mc == 1
+
+
+def test_large_committed_quantity_remains_feasible_under_a_non_convex_cost_curve():
+    """A committed quantity far beyond the devices' flow limits must not be cut off by Mc.
+
+    The non-convex prices (summed upwards price below summed downwards price) activate the commitment-sign constraints,
+    in which Mc caps the deviations.
+    Before Mc accounted for the committed quantity, the required deviation exceeded Mc and the problem was infeasible.
+    """
+    commitment = FlowCommitment(
+        name="energy",
+        quantity=-100,
+        upwards_deviation_price=-1,
+        downwards_deviation_price=1,
+        index=make_index(),
+    )
+    schedule, costs, results, model = device_scheduler(
+        device_constraints=[make_device_constraints(0.5)],
+        ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION),
+        commitments=[commitment],
+    )
+    assert results.solver.termination_condition == "optimal"
+    # The upwards deviation earns 1 per unit, so the device consumes at full power.
+    np.testing.assert_allclose(schedule[0].values, 0.5, atol=1e-6)
+    # Each of the 4 steps deviates upwards by 100.5 at price -1.
+    assert costs == pytest.approx(-4 * 100.5)

flexmeasures/data/tests/test_scheduling_simultaneous.py

@@ -134,13 +134,6 @@ def test_create_simultaneous_jobs(
 
     # Define expected costs based on resolution
     expected_total_cost = -3.2775
-    # Aggregate (unscoped) commitment semantics (issue #2379): the sample commitment
-    # rewarding supply binds the site's *aggregate* flow, so it stays inactive while
-    # the site is net-consuming and does not bias the per-device dispatch. Under the
-    # earlier per-device binding it wrongly rewarded the battery's supply, shifting the
-    # EV/battery split (EV costs were €2.3125); the total cost is unchanged either way.
-    expected_ev_costs = 2.2375
-    expected_battery_costs = expected_total_cost - expected_ev_costs
 
     # Check costs
     np.testing.assert_approx_equal(
@@ -149,6 +142,15 @@ def test_create_simultaneous_jobs(
         4,
         f"Total costs should be €{expected_total_cost}, got €{total_cost}",
     )
+    # Fairness benchmark: how the total cost is split between the devices.
+    # The optimum is degenerate (both devices face the same prices),
+    # so the split records which vertex the solver lands on, and model changes may move it.
+    # When one does, update the expectation deliberately, as a record of the change's fairness impact.
+    # PR #2411 (tightened Mc, changing the sign-constraint coefficients) moved the EV costs from €2.2375 to €2.3125,
+    # which happens to match the split seen under the old per-device commitment binding of issue #2379;
+    # that semantics is now pinned directly by the commitment-cost assertion below, rather than by the split.
+    expected_ev_costs = 2.3125
+    expected_battery_costs = expected_total_cost - expected_ev_costs
     np.testing.assert_approx_equal(
         ev_costs,
         expected_ev_costs,
@@ -161,6 +163,19 @@ def test_create_simultaneous_jobs(
         4,
         f"Battery costs should be €{expected_battery_costs}, got €{battery_costs}",
     )
+    # Aggregate (unscoped) commitment semantics (issue #2379):
+    # the sample commitment rewarding supply binds the site's *aggregate* flow,
+    # so it collects a reward only where the site as a whole net-produces.
+    # Under the earlier per-device binding it wrongly rewarded the battery's own supply,
+    # also while the site was net-consuming.
+    np.testing.assert_approx_equal(
+        job.meta["scheduler_info"]["commitment_costs"][
+            "a sample commitment rewarding supply"
+        ],
+        -0.038,
+        4,
+        "the aggregate flow commitment should be rewarded for the site's aggregate net supply only",
+    )
     np.testing.assert_approx_equal(
         job.meta["scheduler_info"]["commitment_costs"]["electricity net energy"],
         expected_total_cost,