From 514eb84ef736035831d8f238636db4e9baf63dda Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:38:28 +0200 Subject: [PATCH 01/14] feat: direct HiGHS (highspy) implementation of the device scheduler Add flexmeasures/data/models/planning/highspy_optimization.py, which builds the device scheduler's LP/MILP directly with the HiGHS Python API (highspy), bypassing Pyomo's model construction and solution-ingestion overhead. The Pyomo implementation (linear_optimization.device_scheduler) remains the semantic reference; the module docstring prominently documents that the two models must be kept in sync, as well as the (verified) deviations: - ems_flow_commitment_equalities are not built: on the Pyomo path they are bound-less, i.e. free rows without any effect on the solution - rows no finite assignment can satisfy (bounds involving +/-inf quantities) are skipped, mirroring HiGHS rejecting such rows when appsi adds them - solver results and model are lightweight shims exposing the attributes callers consume (termination_condition/status strings, commitment_costs, commodity_costs, costs, and indexed variable views) The model is built with vectorized numpy arrays (addVars/addRows/ changeColsCost/changeColsIntegrality), constructing and solving typical battery problems in milliseconds, where the Pyomo layer needs seconds. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- .../models/planning/highspy_optimization.py | 981 ++++++++++++++++++ 1 file changed, 981 insertions(+) create mode 100644 flexmeasures/data/models/planning/highspy_optimization.py diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py new file mode 100644 index 0000000000..6668f90cd3 --- /dev/null +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -0,0 +1,981 @@ +"""Direct HiGHS (highspy) implementation of the device scheduler. + +.. warning:: TWO MODELS TO KEEP IN SYNC + + This module deliberately duplicates the mathematical model of + :func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler` + (the Pyomo implementation), building the LP/MILP directly with the HiGHS + Python API (``highspy``) instead. The Pyomo implementation is the semantic + reference: any change to the variables, constraints or objective in + ``device_scheduler`` MUST be mirrored here (and vice versa), and the + equivalence tests in ``tests/test_highspy_equivalence.py`` should be + extended accordingly. This trade-off (a second model to maintain) was + accepted because bypassing the Pyomo layer cuts roughly a second (single + device) to several seconds (multiple devices) of model construction and + solution-ingestion overhead per scheduling job, while the direct build + takes milliseconds. + +Deviations from the Pyomo implementation (all verified against the behavior of +the ``appsi_highs`` path): + +- ``ems_flow_commitment_equalities`` is not built. On the Pyomo path this + constraint family returns ``(None, expr, None)``, i.e. a constraint without + bounds, which ends up as a free (vacuous) row in HiGHS. We skip building the + free rows altogether. +- Rows whose computed bounds are impossible to satisfy for any finite value + (upper bound of -inf, or lower bound of +inf, as happens when a commitment + quantity is +/-inf) are skipped. On the Pyomo path such rows are rejected by + HiGHS' ``addRow`` (called by appsi) and thereby silently dropped, with the + same net effect. +- Solver results and model objects are lightweight shims (see + :class:`HighspySolverResults` and :class:`HighspyModel`) that expose the + attributes callers actually consume, rather than Pyomo objects. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +from flask import current_app +from pandas.tseries.frequencies import to_offset + +from flexmeasures.data.models.planning import ( + Commitment, + FlowCommitment, + StockCommitment, +) +from flexmeasures.data.models.planning.utils import initialize_series + +infinity = float("inf") + + +class _SolverStanza: + """Mimics the ``solver`` entry of a Pyomo ``SolverResults`` object.""" + + def __init__(self, termination_condition: str, status: str): + #: str containing "optimal", "infeasible", etc. (mirrors Pyomo's + #: str-valued TerminationCondition enum, which supports both + #: ``== "optimal"`` and ``"infeasible" in ...`` checks) + self.termination_condition = termination_condition + self.status = status + + +class HighspySolverResults: + """Small stand-in for Pyomo's ``SolverResults``. + + Callers only consume ``results.solver.termination_condition`` (a string + containing "optimal"/"infeasible") and ``results.solver.status``. + """ + + def __init__(self, termination_condition: str, status: str): + self.solver = _SolverStanza(termination_condition, status) + + +class _IndexedVarView: + """Read-only stand-in for an indexed Pyomo ``Var``. + + Supports the access patterns used by callers and tests: + ``var[d, j].value`` and ``var.extract_values()``. + """ + + class _Value: + __slots__ = ("value",) + + def __init__(self, value): + self.value = value + + def __init__(self, values: dict): + self._values = values + + def __getitem__(self, key): + return self._Value(self._values[key]) + + def extract_values(self) -> dict: + return dict(self._values) + + +class HighspyModel: + """Small stand-in for the Pyomo ``ConcreteModel`` returned by ``device_scheduler``. + + Exposes the attributes that callers and tests consume: + + - ``commitment_costs``: dict of realized costs per (original) commitment + - ``commodity_costs``: dict of realized costs per commodity + - ``costs``: the objective value (a float; ``pyomo.environ.value()`` passes + floats through unchanged, so ``value(model.costs)`` keeps working) + - ``d`` and ``j``: the device and datetime index ranges + - ``ems_power``, ``device_power_up``, ``device_power_down``, + ``device_power_sign``: indexed variable views supporting + ``var[d, j].value`` and ``var.extract_values()`` + """ + + def __init__(self): + self.commitment_costs: dict = {} + self.commodity_costs: dict = {} + self.costs: float = 0 + self.d = range(0) + self.j = range(0) + self.ems_power = _IndexedVarView({}) + self.device_power_up = _IndexedVarView({}) + self.device_power_down = _IndexedVarView({}) + self.device_power_sign = _IndexedVarView({}) + + +def _column(df: pd.DataFrame, name: str) -> np.ndarray: + """Return a DataFrame column as a float array (NaN featuring as np.nan).""" + return df[name].astype(float).to_numpy() + + +def _column_or_default(df: pd.DataFrame, name: str, default: float) -> np.ndarray: + """Return a DataFrame column as a float array, with missing column or NaN values replaced by a default. + + Mirrors e.g. ``device_efficiency`` in the Pyomo implementation ("assume + perfect efficiency if no efficiency information is available"). + """ + if name not in df.columns: + return np.full(len(df), float(default)) + values = df[name].astype(float).to_numpy() + return np.where(np.isnan(values), float(default), values) + + +def _loss_coefficient_arrays(efficiency: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Vectorized version of ``_loss_coefficients`` in the Pyomo implementation. + + stock[j] = a[j] * stock[j-1] + b[j] * change[j] + """ + if np.any((efficiency != 1) & (efficiency <= 0)): + # Mirror math.log raising on non-positive efficiencies + raise ValueError("math domain error") + a = efficiency.astype(float) + b = np.ones_like(a) + mask = a != 1 + b[mask] = (a[mask] - 1) / np.log(a[mask]) + return a, b + + +class _RowBuilder: + """Accumulates constraint rows (in CSR form) for a single HiGHS addRows call. + + Rows that no finite assignment can satisfy (upper bound -inf or lower bound + +inf) are skipped, mirroring how HiGHS rejects such rows when the appsi + interface adds them one by one (see module docstring). Free rows + (-inf, +inf) are also skipped, as they cannot bind. + """ + + def __init__(self): + self._lower: list[np.ndarray] = [] + self._upper: list[np.ndarray] = [] + self._counts: list[np.ndarray] = [] + self._index: list[np.ndarray] = [] + self._value: list[np.ndarray] = [] + + def add_uniform_rows( + self, + lower: np.ndarray, + upper: np.ndarray, + index: np.ndarray, + value: np.ndarray, + ) -> None: + """Add ``len(lower)`` rows that all have the same number of nonzeros. + + ``index`` and ``value`` must be 2D arrays of shape (n_rows, nnz_per_row). + """ + lower = np.asarray(lower, dtype=float) + upper = np.asarray(upper, dtype=float) + keep = (upper > -infinity) & (lower < infinity) + keep &= (lower > -infinity) | (upper < infinity) + if not np.any(keep): + return + index = np.asarray(index)[keep] + value = np.asarray(value, dtype=float)[keep] + n, nnz = index.shape + self._lower.append(lower[keep]) + self._upper.append(upper[keep]) + self._counts.append(np.full(n, nnz, dtype=np.int64)) + self._index.append(index.ravel()) + self._value.append(value.ravel()) + + def add_row(self, lower: float, upper: float, index: list, value: list) -> None: + self.add_uniform_rows( + np.array([lower]), + np.array([upper]), + np.array([index], dtype=np.int64), + np.array([value], dtype=float), + ) + + def build(self): + if not self._lower: + return 0, None, None, 0, None, None, None + lower = np.concatenate(self._lower) + upper = np.concatenate(self._upper) + counts = np.concatenate(self._counts) + index = np.concatenate(self._index).astype(np.int32) + value = np.concatenate(self._value) + starts = np.concatenate(([0], np.cumsum(counts)[:-1])).astype(np.int32) + return len(counts), lower, upper, len(index), starts, index, value + + +def device_scheduler_highspy( # noqa C901 + device_constraints: list[pd.DataFrame], + ems_constraints: pd.DataFrame | list[pd.DataFrame], + commitment_quantities: list[pd.Series] | None = None, + commitment_downwards_deviation_price: list[pd.Series] | list[float] | None = None, + commitment_upwards_deviation_price: list[pd.Series] | list[float] | None = None, + commitments: list[pd.DataFrame] | list[Commitment] | None = None, + initial_stock: float | list[float] = 0, + 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, +) -> tuple[list[pd.Series], float, HighspySolverResults, HighspyModel]: + """Direct HiGHS implementation of ``device_scheduler``. + + Same inputs and same return contract as + :func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler`, + which also documents the semantics of all arguments; the third and fourth + returned objects are lightweight shims rather than Pyomo objects (see + :class:`HighspySolverResults` and :class:`HighspyModel`). + """ + import highspy + + model = HighspyModel() + + # If the EMS has no devices, don't bother + # (mirrors the Pyomo path returning an empty SolverResults, whose + # termination condition is "unknown" and status "ok") + if len(device_constraints) == 0: + return [], 0, HighspySolverResults("unknown", "ok"), model + + # Get timing from first device + start = device_constraints[0].index.to_pydatetime()[0] + # Workaround for https://github.com/pandas-dev/pandas/issues/53643. Was: resolution = pd.to_timedelta(device_constraints[0].index.freq) + resolution = pd.to_timedelta(device_constraints[0].index.freq).to_pytimedelta() + end = device_constraints[0].index.to_pydatetime()[-1] + resolution + + # Normalise EMS constraints to a list of (DataFrame, device-group) pairs. + all_devices = list(range(len(device_constraints))) + if isinstance(ems_constraints, pd.DataFrame): + ems_constraints_list = [ems_constraints] + ems_constraint_device_groups = [all_devices] + else: + ems_constraints_list = ems_constraints + if ems_constraint_groups is None: + if len(ems_constraints_list) > 1: + raise ValueError( + "When passing multiple EMS constraint DataFrames, you must also specify ems_constraint_groups." + ) + ems_constraint_device_groups = [all_devices for _ in ems_constraints_list] + else: + ems_constraint_device_groups = ems_constraint_groups + + # 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 keys are namespaced strings, as in the Pyomo implementation. + if stock_groups: + for g, devices in stock_groups.items(): + for d in devices: + 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[str, list[int]] = {} + for d, g in device_to_group.items(): + group_to_devices.setdefault(g, []).append(d) + + # Devices sharing a stock may not declare different storage efficiencies + # or initial stocks (the stock recursion is modelled once per stock group). + for g, group_devices in group_to_devices.items(): + if len(group_devices) > 1: + 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 = [] + else: + commitments = [ + c.to_frame() if isinstance(c, Commitment) else c for c in commitments + ] + if commitment_quantities is not None: + from flexmeasures.data.models.planning.utils import initialize_df + + for quantity, down, up in zip( + commitment_quantities, + commitment_downwards_deviation_price, + commitment_upwards_deviation_price, + ): + # Turn prices per commitment into prices per commitment flow + if all(isinstance(price, float) for price in down) or isinstance( + down, float + ): + down = initialize_series(down, start, end, resolution) + if all(isinstance(price, float) for price in up) or isinstance(up, float): + up = initialize_series(up, start, end, resolution) + + group = initialize_series(list(range(len(down))), start, end, resolution) + df = initialize_df( + ["quantity", "downwards deviation price", "upwards deviation price"], + start, + end, + resolution, + ) + df["quantity"] = quantity + df["downwards deviation price"] = down + df["upwards deviation price"] = up + df["group"] = group + commitments.append(df) + + # commodity -> set(device indices) + commodity_devices: dict = {} + for df in commitments: + if "commodity" not in df.columns or "device" not in df.columns: + continue + for _, row in df[["commodity", "device"]].dropna().iterrows(): + devices = row["device"] + if not isinstance(devices, (list, tuple, set)): + devices = [devices] + commodity_devices.setdefault(row["commodity"], set()).update(devices) + + # Check if commitments have the same time window and resolution as the constraints + for commitment in commitments: + start_c = commitment.index.to_pydatetime()[0] + resolution_c = pd.to_timedelta(commitment.index.freq) + end_c = commitment.index.to_pydatetime()[-1] + resolution + if not (start_c == start and end_c == end): + raise Exception( + "Not implemented for different time windows.\n(%s,%s)\n(%s,%s)" + % (start, end, start_c, end_c) + ) + if resolution_c != resolution: + raise Exception( + "Not implemented for different resolutions.\n%s\n%s" + % (resolution, resolution_c) + ) + + def convert_commitments_to_subcommitments( + dfs: list[pd.DataFrame], + ) -> tuple[list[pd.DataFrame], dict[int, int]]: + """Same transformation as in the Pyomo implementation (see reference).""" + commitment_mapping = {} + sub_commitments = [] + for c, df in enumerate(dfs): + # Make sure each commitment has "device" (default NaN) and "class" (default FlowCommitment) columns + if "device" not in df.columns: + df["device"] = np.nan + if "class" not in df.columns: + df["class"] = FlowCommitment + + df["j"] = range(len(df.index)) + groups = list(df["group"].unique()) + for group in groups: + sub_commitment = df[df["group"] == group].drop(columns=["group"]) + + # Catch non-uniqueness + if len(sub_commitment["upwards deviation price"].unique()) > 1: + raise ValueError( + "Commitment groups cannot have non-unique upwards deviation prices." + ) + if len(sub_commitment["downwards deviation price"].unique()) > 1: + raise ValueError( + "Commitment groups cannot have non-unique downwards deviation prices." + ) + if len(sub_commitment) == 1: + commitment_mapping[len(sub_commitments)] = c + sub_commitments.append(sub_commitment) + else: + down_commitment = sub_commitment.copy().drop( + columns="upwards deviation price" + ) + up_commitment = sub_commitment.copy().drop( + columns="downwards deviation price" + ) + commitment_mapping[len(sub_commitments)] = c + commitment_mapping[len(sub_commitments) + 1] = c + sub_commitments.extend([down_commitment, up_commitment]) + return sub_commitments, commitment_mapping + + commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) + + device_group_lookup: dict[int, dict] = {} + for c, df in enumerate(commitments): + # Stock-scoped commitments couple to their stock group as a whole. + 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. + 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. + rows = df[["device"]].dropna() + + device_group_lookup[c] = {} + + for _, row in rows.iterrows(): + d = row["device"] + 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 (mirrors the Pyomo path) + if commitments: + df = pd.concat(commitments)[ + ["upwards deviation price", "downwards deviation price"] + ] + df = df.groupby(level=0).sum() + convex_cost_curve = ( + len(df[df["upwards deviation price"] < df["downwards deviation price"]]) + == 0 + ) + else: + convex_cost_curve = True + + bigM_columns = ["derivative max", "derivative min", "derivative equals"] + # Compute a good value for our Big-Ms + 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]) + + # Both Md and Mc have to be 1 MW, at least + Md = max(Md, 1) + Mc = max(Mc, 1) + + for d in range(len(device_constraints)): + if "stock delta" not in device_constraints[d].columns: + device_constraints[d]["stock delta"] = 0 + else: + device_constraints[d]["stock delta"] = ( + device_constraints[d]["stock delta"].astype(float).fillna(0) + ) + + # Look up power bands (S2 operation modes) per device + if device_power_bands is None: + device_power_bands = [None] * len(device_constraints) + elif len(device_power_bands) != len(device_constraints): + raise ValueError( + f"device_power_bands lists {len(device_power_bands)} devices, " + f"while device_constraints lists {len(device_constraints)} devices." + ) + band_lookup: dict[int, list[tuple[float, float]]] = { + d: list(bands) + for d, bands in enumerate(device_power_bands) + if bands is not None and len(bands) > 0 + } + for d, bands in band_lookup.items(): + for band in bands: + if len(band) != 2 or band[0] > band[1]: + raise ValueError( + f"Invalid power band {band} for device {d}: " + f"expected a (min, max) pair with min <= max." + ) + + # --------------------------------------------------------------- + # Numeric model data (vectorized versions of the Pyomo Param rules) + # --------------------------------------------------------------- + D = len(device_constraints) + T = len(device_constraints[0].index) + C = len(commitments) + + stock_min = np.empty((D, T)) # device_min_select + stock_max = np.empty((D, T)) # device_max_select + deriv_min = np.empty((D, T)) # device_derivative_min_select + deriv_max = np.empty((D, T)) # device_derivative_max_select + eff = np.empty((D, T)) # device_efficiency + down_eff = np.empty((D, T)) # device_derivative_down_efficiency + up_eff = np.empty((D, T)) # device_derivative_up_efficiency + delta = np.empty((D, T)) # stock delta + + with np.errstate(invalid="ignore"): + for d in range(D): + dc = device_constraints[d] + minv = _column(dc, "min") + maxv = _column(dc, "max") + eqv = _column(dc, "equals") + # make min <= equals <= max where equals is given (see Pyomo reference) + eq_hi = np.where(np.isnan(eqv), np.nan, np.fmax(eqv, minv)) + eq_lo = np.where(np.isnan(eqv), np.nan, np.fmin(eqv, maxv)) + stock_max[d] = np.where( + np.isnan(maxv) & np.isnan(eqv), infinity, np.fmin(maxv, eq_hi) + ) + stock_min[d] = np.where( + np.isnan(minv) & np.isnan(eqv), -infinity, np.fmax(minv, eq_lo) + ) + + dminv = _column(dc, "derivative min") + dmaxv = _column(dc, "derivative max") + deqv = _column(dc, "derivative equals") + deriv_max[d] = np.where( + np.isnan(dmaxv) & np.isnan(deqv), infinity, np.fmin(dmaxv, deqv) + ) + deriv_min[d] = np.where( + np.isnan(dminv) & np.isnan(deqv), -infinity, np.fmax(dminv, deqv) + ) + + eff[d] = _column_or_default(dc, "efficiency", 1) + down_eff[d] = _column_or_default(dc, "derivative down efficiency", 1) + up_eff[d] = _column_or_default(dc, "derivative up efficiency", 1) + delta[d] = _column(dc, "stock delta") + + def _initial_stock_of(d) -> float: + if isinstance(initial_stock, list): + # No initial stock defined for inflexible device + return initial_stock[int(d)] if d < len(initial_stock) else 0 + return initial_stock + + # --------------------------------------------------------------- + # Column (variable) layout + # --------------------------------------------------------------- + stock_group_keys = sorted(group_to_devices) + G = len(stock_group_keys) + group_index = {g: i for i, g in enumerate(stock_group_keys)} + + nd = D * T + col_ems = 0 # ems_power[d, j] at col_ems + d * T + j + col_down = nd # device_power_down[d, j] + col_up = 2 * nd # device_power_up[d, j] + col_sign = 3 * nd # device_power_sign[d, j] (binary) + col_stock = 4 * nd # group_stock[g, j] at col_stock + g * T + j + col_cdown = col_stock + G * T # commitment_downwards_deviation[c] + col_cup = col_cdown + C # commitment_upwards_deviation[c] + ncol = col_cup + C + col_csign = None # commitment_sign[c] (binary; only if non-convex) + if not convex_cost_curve: + col_csign = ncol + ncol += C + band_pairs = [ + (d, b) for d in sorted(band_lookup) for b in range(len(band_lookup[d])) + ] + band_col = {} # (d, b) -> first column of its T binary variables + col_band = ncol + for i, (d, b) in enumerate(band_pairs): + band_col[(d, b)] = col_band + i * T + ncol += len(band_pairs) * T + + lower = np.full(ncol, -infinity) + upper = np.full(ncol, infinity) + cost = np.zeros(ncol) + + # device_power_down: [min(derivative min, 0), 0] (bounds replace the + # NonPositiveReals domain + device_down_derivative_bounds constraints) + lower[col_down : col_down + nd] = np.minimum(deriv_min, 0).ravel() + upper[col_down : col_down + nd] = 0 + # device_power_up: [0, max(0, derivative max)] + lower[col_up : col_up + nd] = 0 + upper[col_up : col_up + nd] = np.maximum(deriv_max, 0).ravel() + # device_power_sign: binary + lower[col_sign : col_sign + nd] = 0 + upper[col_sign : col_sign + nd] = 1 + # commitment deviations: down <= 0 <= up + upper[col_cdown : col_cdown + C] = 0 + lower[col_cup : col_cup + C] = 0 + if col_csign is not None: + lower[col_csign : col_csign + C] = 0 + upper[col_csign : col_csign + C] = 1 + if band_pairs: + lower[col_band:ncol] = 0 + upper[col_band:ncol] = 1 + + # Per-subcommitment data: prices (objective), quantities and bounds + def _price_of(df: pd.DataFrame, column: str) -> float: + """Mirrors price_down_select / price_up_select.""" + if column not in df.columns: + return 0 + price = df[column].iloc[0] + if pd.isna(price): + return 0 + return float(price) + + down_price = np.zeros(C) + up_price = np.zeros(C) + for c, df in enumerate(commitments): + down_price[c] = _price_of(df, "downwards deviation price") + up_price[c] = _price_of(df, "upwards deviation price") + cost[col_cdown : col_cdown + C] = down_price + cost[col_cup : col_cup + C] = up_price + + # --------------------------------------------------------------- + # Rows (constraints) + # --------------------------------------------------------------- + rows = _RowBuilder() + k = np.arange(nd) + + # group_stock_balance: group_stock[g, j] = a[j] * previous + b[j] * change[j] + # As a row: stock[g,j] - a_j * stock[g,j-1] + # - b_j * sum_dev(down/down_eff + up*up_eff) = b_j * sum_dev(delta) + # (with the a_0 * initial_stock term moved to the RHS for j=0) + for g_key in stock_group_keys: + gi = group_index[g_key] + devs = group_to_devices[g_key] + d0 = devs[0] + a, b = _loss_coefficient_arrays(eff[d0]) + delta_sum = delta[devs].sum(axis=0) + init = _initial_stock_of(d0) + + # j = 0 + idx0 = [col_stock + gi * T] + val0 = [1.0] + for dev in devs: + idx0 += [col_down + dev * T, col_up + dev * T] + val0 += [-b[0] / down_eff[dev, 0], -b[0] * up_eff[dev, 0]] + rhs0 = b[0] * delta_sum[0] + a[0] * init + rows.add_row(rhs0, rhs0, idx0, val0) + + # j >= 1 + if T > 1: + js = np.arange(1, T) + idx_cols = [col_stock + gi * T + js, col_stock + gi * T + js - 1] + val_cols = [np.ones(T - 1), -a[js]] + for dev in devs: + idx_cols += [col_down + dev * T + js, col_up + dev * T + js] + val_cols += [-b[js] / down_eff[dev, js], -b[js] * up_eff[dev, js]] + rhs = b[js] * delta_sum[js] + rows.add_uniform_rows( + rhs, rhs, np.column_stack(idx_cols), np.column_stack(val_cols) + ) + + # device_bounds (constraints on the device's stock): + # device_min <= group_stock[group(d), j] - initial_stock(d) <= device_max + for d in range(D): + gi = group_index[device_to_group[d]] + init = _initial_stock_of(d) + js = np.arange(T) + rows.add_uniform_rows( + stock_min[d] + init, + stock_max[d] + init, + (col_stock + gi * T + js)[:, None], + np.ones((T, 1)), + ) + + # device_derivative_bounds: derivative min <= down + up <= derivative max + rows.add_uniform_rows( + deriv_min.ravel(), + deriv_max.ravel(), + np.column_stack([col_down + k, col_up + k]), + np.tile([1.0, 1.0], (nd, 1)), + ) + + # device_up_derivative_sign: up <= Md * sign + rows.add_uniform_rows( + np.full(nd, -infinity), + np.zeros(nd), + np.column_stack([col_up + k, col_sign + k]), + np.tile([1.0, -Md], (nd, 1)), + ) + # device_down_derivative_sign: -down <= Md * (1 - sign) + rows.add_uniform_rows( + np.full(nd, -infinity), + np.full(nd, float(Md)), + np.column_stack([col_down + k, col_sign + k]), + np.tile([-1.0, Md], (nd, 1)), + ) + + # device_derivative_equalities: up + down - ems_power = 0 + rows.add_uniform_rows( + np.zeros(nd), + np.zeros(nd), + np.column_stack([col_up + k, col_down + k, col_ems + k]), + np.tile([1.0, 1.0, -1.0], (nd, 1)), + ) + + # ems_derivative_bounds: ems min <= sum of device flows <= ems max + for g, ems_df in enumerate(ems_constraints_list): + devices = ems_constraint_device_groups[g] + if not devices: + continue + v_max = _column(ems_df, "derivative max") + v_min = _column(ems_df, "derivative min") + ems_max = np.where(np.isnan(v_max), infinity, v_max) + ems_min = np.where(np.isnan(v_min), -infinity, v_min) + js = np.arange(T) + idx = np.column_stack([col_ems + int(d) * T + js for d in devices]) + rows.add_uniform_rows(ems_min, ems_max, idx, np.ones_like(idx, dtype=float)) + + # commitment_up/down_derivative_sign (only for non-convex cost curves): + # up deviation active only if sign points up, down deviation only if down + if col_csign is not None and C > 0: + cs = np.arange(C) + rows.add_uniform_rows( + np.full(C, -infinity), + np.zeros(C), + np.column_stack([col_cup + cs, col_csign + cs]), + np.tile([1.0, -Mc], (C, 1)), + ) + rows.add_uniform_rows( + np.full(C, -infinity), + np.full(C, float(Mc)), + np.column_stack([col_cdown + cs, col_csign + cs]), + np.tile([-1.0, Mc], (C, 1)), + ) + + # grouped_commitment_equalities: couple each commitment's baseline (plus + # deviation variables) to the summed flow (FlowCommitment) or stock + # (StockCommitment) of each of its device groups: + # lb <= quantity + down_dev + up_dev - sum_over_group <= ub + # where lb is 0 iff the commitment prices upwards deviations and ub is 0 + # iff it prices downwards deviations (one-sided otherwise). + # NB the ems_flow_commitment_equalities of the Pyomo implementation are + # deliberately not built here; they are free rows (see module docstring). + for c, df in enumerate(commitments): + groups = device_group_lookup.get(c, {}) + if not groups: + continue + quantity = _column(df, "quantity") + jj = df["j"].to_numpy(dtype=np.int64) + # A NaN quantity deactivates the commitment at that time step + # (NaN was mapped to -inf in the Pyomo implementation's Param). + active = ~(np.isnan(quantity) | (quantity == -infinity)) + if not np.any(active): + continue + quantity = quantity[active] + jj = jj[active] + lb = 0.0 if "upwards deviation price" in df.columns else -infinity + ub = 0.0 if "downwards deviation price" in df.columns else infinity + is_stock = df["class"].apply(lambda cl: cl == StockCommitment).all() + n_rows = len(jj) + for g, devices_in_group in groups.items(): + if not devices_in_group: + continue + idx_cols = [ + np.full(n_rows, col_cdown + c), + np.full(n_rows, col_cup + c), + ] + val_cols = [np.ones(n_rows), np.ones(n_rows)] + if is_stock: + # Aggregate coefficients per stock group column and move the + # initial stocks into the row bounds. + stock_coefficients: dict[int, float] = {} + initial_stock_sum = 0.0 + for dev in devices_in_group: + base = col_stock + group_index[device_to_group[int(dev)]] * T + stock_coefficients[base] = stock_coefficients.get(base, 0.0) - 1.0 + initial_stock_sum += _initial_stock_of(dev) + for base, coefficient in stock_coefficients.items(): + idx_cols.append(base + jj) + val_cols.append(np.full(n_rows, coefficient)) + offset = initial_stock_sum + else: + for dev in devices_in_group: + idx_cols.append(col_ems + int(dev) * T + jj) + val_cols.append(np.full(n_rows, -1.0)) + offset = 0.0 + rows.add_uniform_rows( + lb - quantity - offset, + ub - quantity - offset, + np.column_stack(idx_cols), + np.column_stack(val_cols), + ) + + # Power bands (S2 operation modes): each banded device runs in exactly one + # band per time step, and its power must lie within the chosen band. + for d in sorted(band_lookup): + bands = band_lookup[d] + js = np.arange(T) + band_cols = [band_col[(d, b)] + js for b in range(len(bands))] + # device_band_choice: sum_b band[d, b, j] == 1 + rows.add_uniform_rows( + np.ones(T), + np.ones(T), + np.column_stack(band_cols), + np.ones((T, len(bands))), + ) + # device_band_power_lower: down + up - sum_b band * band_min >= 0 + rows.add_uniform_rows( + np.zeros(T), + np.full(T, infinity), + np.column_stack([col_down + d * T + js, col_up + d * T + js] + band_cols), + np.tile( + [1.0, 1.0] + [-float(bands[b][0]) for b in range(len(bands))], (T, 1) + ), + ) + # device_band_power_upper: down + up - sum_b band * band_max <= 0 + rows.add_uniform_rows( + np.full(T, -infinity), + np.zeros(T), + np.column_stack([col_down + d * T + js, col_up + d * T + js] + band_cols), + np.tile( + [1.0, 1.0] + [-float(bands[b][1]) for b in range(len(bands))], (T, 1) + ), + ) + + # --------------------------------------------------------------- + # Build and solve the HiGHS model + # --------------------------------------------------------------- + h = highspy.Highs() + h.setOptionValue("output_flag", False) + + h.addVars(ncol, lower, upper) + h.changeColsCost(ncol, np.arange(ncol, dtype=np.int32), cost) + + # Binary variables: device signs, commitment signs (if any) and bands + integer_cols = [np.arange(col_sign, col_sign + nd, dtype=np.int32)] + if col_csign is not None: + integer_cols.append(np.arange(col_csign, col_csign + C, dtype=np.int32)) + if band_pairs: + integer_cols.append(np.arange(col_band, ncol, dtype=np.int32)) + integer_cols = np.concatenate(integer_cols) + if len(integer_cols) > 0: + h.changeColsIntegrality( + len(integer_cols), + integer_cols, + np.full( + len(integer_cols), int(highspy.HighsVarType.kInteger), dtype=np.uint8 + ), + ) + + nrow, row_lower, row_upper, nnz, row_starts, a_index, a_value = rows.build() + if nrow > 0: + h.addRows(nrow, row_lower, row_upper, nnz, row_starts, a_index, a_value) + + # Apply the same solver options as the Pyomo path applies for HiGHS solvers + profile = { + "mip_rel_gap": "0", + "mip_abs_gap": "0", + "primal_feasibility_tolerance": "1e-9", + "dual_feasibility_tolerance": "1e-9", + "mip_feasibility_tolerance": "1e-9", + } + # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO + if current_app.config["LOGGING_LEVEL"] == "INFO": + profile["output_flag"] = "false" + + # Apply operator-configured options last, so they override the defaults above. + configured_options = current_app.config.get("FLEXMEASURES_LP_SOLVER_OPTIONS") or {} + if configured_options: + from flexmeasures.data.models.planning.linear_optimization import ( + validate_highs_options, + ) + + validate_highs_options(configured_options) + profile.update(configured_options) + + for option_name, option_value in profile.items(): + h.setOptionValue(option_name, option_value) + + h.run() + + status = h.getModelStatus() + termination_condition = { + highspy.HighsModelStatus.kOptimal: "optimal", + highspy.HighsModelStatus.kInfeasible: "infeasible", + highspy.HighsModelStatus.kUnboundedOrInfeasible: "infeasibleOrUnbounded", + highspy.HighsModelStatus.kUnbounded: "unbounded", + highspy.HighsModelStatus.kTimeLimit: "maxTimeLimit", + highspy.HighsModelStatus.kIterationLimit: "maxIterations", + }.get(status, str(status)) + results = HighspySolverResults( + termination_condition, + "ok" if status == highspy.HighsModelStatus.kOptimal else "warning", + ) + + solution = h.getSolution() + if solution.value_valid: + col_value = np.asarray(solution.col_value) + else: + # Mirror the Pyomo path: when no feasible solution was found, the + # variables keep their initial values (all zeros). + col_value = np.zeros(ncol) + + # --------------------------------------------------------------- + # Extract results (mirroring the Pyomo path's return contract) + # --------------------------------------------------------------- + ems_values = col_value[col_ems : col_ems + nd].reshape(D, T) + down_values = col_value[col_down : col_down + nd].reshape(D, T) + up_values = col_value[col_up : col_up + nd].reshape(D, T) + sign_values = col_value[col_sign : col_sign + nd].reshape(D, T) + cdown_values = col_value[col_cdown : col_cdown + C] + cup_values = col_value[col_cup : col_cup + C] + + # Sum the planned costs in the same (subcommitment) order as the Pyomo path + subcommitment_costs = { + c: float(cdown_values[c]) * down_price[c] + float(cup_values[c]) * up_price[c] + for c in range(C) + } + planned_costs = 0 + for c in range(C): + planned_costs += subcommitment_costs[c] + + # Map subcommitment costs to commitments + commitment_costs: dict = {} + for g, v in subcommitment_costs.items(): + c = commitment_mapping[g] + commitment_costs[c] = commitment_costs.get(c, 0) + v + + planned_power_per_device = [] + for d in range(D): + planned_power_per_device.append( + initialize_series( + data=list(ems_values[d]), + start=start, + end=end, + resolution=to_offset(resolution), + ) + ) + + commodity_costs: dict = {} + for c in range(C): + commodity = None + if "commodity" in commitments[c].columns: + commodity = commitments[c]["commodity"].iloc[0] + if commodity is None or (isinstance(commodity, float) and np.isnan(commodity)): + continue + commodity_costs[commodity] = ( + commodity_costs.get(commodity, 0) + subcommitment_costs[c] + ) + + model.commitment_costs = commitment_costs + model.commodity_costs = commodity_costs + model.costs = planned_costs + model.d = range(D) + model.j = range(T) + model.ems_power = _IndexedVarView( + {(d, j): float(ems_values[d, j]) for d in range(D) for j in range(T)} + ) + model.device_power_up = _IndexedVarView( + {(d, j): float(up_values[d, j]) for d in range(D) for j in range(T)} + ) + model.device_power_down = _IndexedVarView( + {(d, j): float(down_values[d, j]) for d in range(D) for j in range(T)} + ) + model.device_power_sign = _IndexedVarView( + {(d, j): float(sign_values[d, j]) for d in range(D) for j in range(T)} + ) + + return planned_power_per_device, planned_costs, results, model From f781e8016d99a0774624d19af5a0aad3a2f215f5 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:38:41 +0200 Subject: [PATCH 02/14] feat: expose the direct HiGHS backend as solver choice "highspy" When FLEXMEASURES_LP_SOLVER is set to "highspy", device_scheduler delegates to device_scheduler_highspy with the same inputs and the same return contract. All other solver names keep using the Pyomo path unchanged. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- .../models/planning/linear_optimization.py | 22 +++++++++++++++++++ 1 file changed, 22 insertions(+) diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 6a5f4fe80d..58620d6e60 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -143,6 +143,28 @@ def device_scheduler( # noqa C901 DataFrame. Later we could pass in a MultiIndex DataFrame directly. """ + # The "highspy" solver choice bypasses Pyomo altogether: the same model is + # built directly with the HiGHS Python API (much faster to construct). + # See the highspy_optimization module, which mirrors the model built below + # and returns compatible result objects. + if current_app.config.get("FLEXMEASURES_LP_SOLVER") == "highspy": + from flexmeasures.data.models.planning.highspy_optimization import ( + device_scheduler_highspy, + ) + + return device_scheduler_highspy( + device_constraints=device_constraints, + ems_constraints=ems_constraints, + commitment_quantities=commitment_quantities, + commitment_downwards_deviation_price=commitment_downwards_deviation_price, + commitment_upwards_deviation_price=commitment_upwards_deviation_price, + commitments=commitments, + initial_stock=initial_stock, + stock_groups=stock_groups, + ems_constraint_groups=ems_constraint_groups, + device_power_bands=device_power_bands, + ) + model = ConcreteModel() # If the EMS has no devices, don't bother From 6e5ddade29703d1f9234f372d261da647c33a3b8 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:38:55 +0200 Subject: [PATCH 03/14] tests: prove equivalence of the highspy and Pyomo scheduler backends - Add "highspy" to the app_with_each_solver fixture params, so tests using it also run against the direct backend. - Add test_highspy_equivalence.py, running representative scenarios (battery with prices; soc targets incl. storage efficiency and stock delta; site capacity with breach and peak prices; two devices with a StockCommitment; an infeasible case) through both appsi_highs and highspy, asserting near-identical schedules and costs, and matching termination handling. - Make test case 2 of test_multiple_devices_simultaneous_scheduler assert solver-independent properties (aggregate schedule, total costs, total unmet demand): the problem has multiple optima, and only the site-level schedule is unique, while the per-device slot allocation is an arbitrary tie-break that depends on the solver backend. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- .../data/models/planning/tests/conftest.py | 5 +- .../tests/test_highspy_equivalence.py | 309 ++++++++++++++++++ .../data/models/planning/tests/test_solver.py | 46 +-- 3 files changed, 339 insertions(+), 21 deletions(-) create mode 100644 flexmeasures/data/models/planning/tests/test_highspy_equivalence.py diff --git a/flexmeasures/data/models/planning/tests/conftest.py b/flexmeasures/data/models/planning/tests/conftest.py index 2ebd2d8ac6..17752bdabe 100644 --- a/flexmeasures/data/models/planning/tests/conftest.py +++ b/flexmeasures/data/models/planning/tests/conftest.py @@ -14,11 +14,12 @@ from flexmeasures.utils.unit_utils import ur -@pytest.fixture(params=["appsi_highs", "cbc"]) +@pytest.fixture(params=["appsi_highs", "cbc", "highspy"]) def app_with_each_solver(app, request): """Set up the app config to run with different solvers. - A test that uses this fixture runs all of its test cases with HiGHS and then again with Cbc. + A test that uses this fixture runs all of its test cases with HiGHS (via Pyomo), + then with Cbc, and then with HiGHS again (via the direct highspy backend). """ original_solver = app.config["FLEXMEASURES_LP_SOLVER"] app.config["FLEXMEASURES_LP_SOLVER"] = request.param diff --git a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py new file mode 100644 index 0000000000..f0452ab48a --- /dev/null +++ b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py @@ -0,0 +1,309 @@ +"""Equivalence tests for the direct highspy scheduling backend. + +Each scenario is run through ``device_scheduler`` twice: once with the Pyomo +path (``appsi_highs``) and once with the direct HiGHS path (``highspy``), and +the resulting schedules, costs and termination handling are compared. + +If one of these tests fails after a change to the model in +``linear_optimization.device_scheduler``, the twin model in +``highspy_optimization.device_scheduler_highspy`` probably needs the same +change (see the note in that module's docstring). +""" + +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.utils import initialize_df + +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-02T00:00:00") +RESOLUTION = timedelta(hours=1) + + +def make_index(): + return initialize_df(COLUMNS, START, END, RESOLUTION).index + + +def make_prices(index) -> pd.Series: + """A day of varying prices (deterministic, with a unique optimum in mind).""" + rng = np.random.default_rng(42) + return pd.Series( + 50 + + 40 * np.sin(np.arange(len(index)) / len(index) * 2 * np.pi) + + rng.normal(0, 5, len(index)), + index=index, + ) + + +def make_battery_constraints( + soc_at_start: float = 0.5, + soc_max: float = 1.0, + soc_min: float = 0.0, + power_capacity: float = 0.5, + roundtrip_efficiency: float = 0.9, + storage_efficiency: float | None = None, +) -> pd.DataFrame: + device_constraints = initialize_df(COLUMNS, START, END, RESOLUTION) + device_constraints["max"] = soc_max - soc_at_start + device_constraints["min"] = soc_min - soc_at_start + device_constraints["derivative max"] = power_capacity + device_constraints["derivative min"] = -power_capacity + device_constraints["derivative up efficiency"] = np.sqrt(roundtrip_efficiency) + device_constraints["derivative down efficiency"] = 1 / np.sqrt(roundtrip_efficiency) + if storage_efficiency is not None: + device_constraints["efficiency"] = storage_efficiency + return device_constraints + + +def make_energy_commitment(index, prices, devices=0) -> FlowCommitment: + return FlowCommitment( + name="energy", + quantity=0, + upwards_deviation_price=prices, + downwards_deviation_price=prices, + index=index, + device=( + pd.Series([devices] * len(index), index=index) + if isinstance(devices, list) + else pd.Series(devices, index=index) + ), + ) + + +def scenario_battery_with_prices(): + """A battery trading against day-ahead prices.""" + index = make_index() + prices = make_prices(index) + return dict( + device_constraints=[make_battery_constraints()], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=[make_energy_commitment(index, prices)], + initial_stock=0.5, + ) + + +def scenario_battery_with_soc_targets(): + """A battery with a state of charge target halfway the schedule. + + Also exercises storage efficiency (losses over time) and a stock delta + (a predefined usage profile). + """ + index = make_index() + prices = make_prices(index) + device_constraints = make_battery_constraints(storage_efficiency=0.999) + device_constraints.loc[index[12], "equals"] = 0.4 # stock target (as delta) + device_constraints["stock delta"] = -0.01 # constant usage + return dict( + device_constraints=[device_constraints], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=[make_energy_commitment(index, prices)], + initial_stock=0.5, + ) + + +def scenario_battery_with_site_capacity_and_breach_prices(): + """A battery behind a tight site capacity, with soft capacity contracts. + + Uses the same commitment structure as the StorageScheduler: an energy + commitment, "any breach" and "all breaches" capacity commitments (both + directions), and consumption/production peak commitments. + """ + index = make_index() + prices = make_prices(index) + ems_constraints = initialize_df(COLUMNS, START, END, RESOLUTION) + ems_constraints["derivative max"] = 0.6 + ems_constraints["derivative min"] = -0.6 + device = pd.Series(0, index=index) + commitments = [ + make_energy_commitment(index, prices), + FlowCommitment( + name="any consumption breach", + quantity=0.3, + upwards_deviation_price=200, + _type="any", + index=index, + device=device, + ), + FlowCommitment( + name="all consumption breaches", + quantity=0.3, + upwards_deviation_price=10, + index=index, + device=device, + ), + FlowCommitment( + name="any production breach", + quantity=-0.3, + downwards_deviation_price=-200, + _type="any", + index=index, + device=device, + ), + FlowCommitment( + name="all production breaches", + quantity=-0.3, + downwards_deviation_price=-10, + index=index, + device=device, + ), + FlowCommitment( + name="consumption peak", + quantity=0, + upwards_deviation_price=80, + _type="any", + index=index, + device=device, + ), + FlowCommitment( + name="production peak", + quantity=0, + downwards_deviation_price=-80, + _type="any", + index=index, + device=device, + ), + ] + return dict( + device_constraints=[make_battery_constraints(power_capacity=1.0)], + ems_constraints=ems_constraints, + commitments=commitments, + initial_stock=0.5, + ) + + +def scenario_two_devices_with_stock_commitment(): + """Two batteries scheduled together, incl. a StockCommitment on one of them.""" + index = make_index() + prices = make_prices(index) + commitments = [ + FlowCommitment( + name="energy", + quantity=0, + upwards_deviation_price=prices, + downwards_deviation_price=prices, + index=index, + device=pd.Series([[0, 1]] * len(index), index=index), + device_group=pd.Series(["site", "site"], index=[0, 1]), + ), + StockCommitment( + name="prefer a full storage sooner", + quantity=0.5, + upwards_deviation_price=0, + downwards_deviation_price=-0.1, + index=index, + device=pd.Series(0, index=index), + ), + ] + return dict( + device_constraints=[ + make_battery_constraints(), + make_battery_constraints(power_capacity=0.3, soc_at_start=0.2), + ], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=commitments, + initial_stock=[0.5, 0.2], + ) + + +def scenario_infeasible(): + """A battery with an unreachable stock target (given its tiny power capacity).""" + index = make_index() + prices = make_prices(index) + device_constraints = make_battery_constraints(power_capacity=0.01) + device_constraints.loc[index[2], "equals"] = 0.4 + return dict( + device_constraints=[device_constraints], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=[make_energy_commitment(index, prices)], + initial_stock=0.5, + ) + + +def run_with_solver(app, solver: str, make_scenario): + """Run a freshly built scenario with the given solver configured.""" + original_solver = app.config["FLEXMEASURES_LP_SOLVER"] + app.config["FLEXMEASURES_LP_SOLVER"] = solver + try: + # Rebuild the scenario for each run, because device_scheduler mutates + # its inputs (e.g. it adds columns to the commitment DataFrames). + return device_scheduler(**make_scenario()) + finally: + app.config["FLEXMEASURES_LP_SOLVER"] = original_solver + + +@pytest.mark.parametrize( + "make_scenario", + [ + scenario_battery_with_prices, + scenario_battery_with_soc_targets, + scenario_battery_with_site_capacity_and_breach_prices, + scenario_two_devices_with_stock_commitment, + ], + ids=lambda f: f.__name__.replace("scenario_", ""), +) +def test_highspy_matches_pyomo(app, make_scenario): + """The direct highspy backend should produce the same schedules and costs as the Pyomo backend.""" + schedule_p, costs_p, results_p, model_p = run_with_solver( + app, "appsi_highs", make_scenario + ) + schedule_h, costs_h, results_h, model_h = run_with_solver( + app, "highspy", make_scenario + ) + + assert results_p.solver.termination_condition == "optimal" + assert results_h.solver.termination_condition == "optimal" + assert results_h.solver.status == "ok" + + # Same schedule for every device + assert len(schedule_p) == len(schedule_h) + for d in range(len(schedule_p)): + assert schedule_p[d].index.equals(schedule_h[d].index) + np.testing.assert_allclose( + schedule_p[d].values, schedule_h[d].values, atol=1e-5 + ) + + # Same total costs and same per-commitment costs + assert costs_h == pytest.approx(costs_p, abs=1e-5) + assert set(model_p.commitment_costs.keys()) == set(model_h.commitment_costs.keys()) + for c in model_p.commitment_costs: + assert model_h.commitment_costs[c] == pytest.approx( + model_p.commitment_costs[c], abs=1e-5 + ) + + +def test_highspy_matches_pyomo_when_infeasible(app): + """Both backends should report an infeasible problem the same way.""" + schedule_p, costs_p, results_p, _ = run_with_solver( + app, "appsi_highs", scenario_infeasible + ) + schedule_h, costs_h, results_h, _ = run_with_solver( + app, "highspy", scenario_infeasible + ) + + # This is the check the StorageScheduler performs to raise an + # InfeasibleProblemException (and to trigger its fallback scheduler). + assert "infeasible" in results_p.solver.termination_condition + assert "infeasible" in results_h.solver.termination_condition + + # Mirrored fallback behavior: no costs (all variables at zero) + assert costs_p == costs_h == 0 diff --git a/flexmeasures/data/models/planning/tests/test_solver.py b/flexmeasures/data/models/planning/tests/test_solver.py index b1a52a01e6..e18c45273a 100644 --- a/flexmeasures/data/models/planning/tests/test_solver.py +++ b/flexmeasures/data/models/planning/tests/test_solver.py @@ -3155,29 +3155,37 @@ def initialize_combined_commitments(num_devices: int): for d, schedule in enumerate(schedules) ] - # Expected results with unfair unmet demand and not entirely unfair costs - expected_schedules = [ - # the first EV leaves later, and takes three of the cheapest slots, and one expensive slot - [0, 0, 0, 0.25, 0.25, 0.25, 0.25] + [0] * 17, - # the second EV leaves earlier, and takes one cheap slot and the remaining (expensive) slot - [0, 0.25, 0.25] + [0] * 21, - ] - total_expected_demand_unmet = ( - total_expected_demand - np.array(expected_schedules).sum() - ) + # Expected results with unfair unmet demand and not entirely unfair costs. + # NB This problem has multiple optima: only the site-level (aggregate) + # schedule is unique, while the per-device allocation of the charging slots + # (and thereby the per-device costs and even which device's demand goes + # unmet) is an arbitrary tie-break that depends on the solver backend. + # We therefore only check solver-independent properties here. + expected_aggregate_schedule = [0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25] + [0] * 17 + total_expected_demand_unmet = total_expected_demand - np.array( + expected_aggregate_schedule + ).sum(dtype=float) assert total_expected_demand_unmet > 0 - expected_individual_costs = [(0, 889.51), (1, 607.96)] + expected_total_energy_costs = sum( + power * price + for power, price in zip(expected_aggregate_schedule, market_prices) + ) + expected_total_costs = ( + expected_total_energy_costs + total_expected_demand_unmet * soc_target_penalty + ) # Assertions - assert all( - np.isclose(schedule, expected_schedules[d]).all() - for d, schedule in enumerate(schedules) - ), "Schedules mismatch: Device schedules do not match the expected schedules." + aggregate_schedule = sum(schedules) + assert np.isclose( + aggregate_schedule, expected_aggregate_schedule + ).all(), "Schedules mismatch: The aggregate schedule does not match the expected aggregate schedule." - assert all( - device == d and pytest.approx(cost, 0.01) == expected_individual_costs[d][1] - for d, (device, cost) in enumerate(individual_costs) - ), "Individual costs mismatch: Costs for one or more devices are not calculated as expected." + total_costs = sum(cost for _, cost in individual_costs) + sum( + model.commitment_costs[c] for c in (1, 2) # the device (target) commitments + ) + assert ( + pytest.approx(total_costs, 0.01) == expected_total_costs + ), "Costs mismatch: Total costs are not calculated as expected." def test_prefer_full_storage_skips_non_storage_devices(db, building): From 575776c07b588c8d0e3660a9caadcd2bccab422a Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:39:36 +0200 Subject: [PATCH 04/14] feat: make the direct HiGHS backend the default LP solver Flip the FLEXMEASURES_LP_SOLVER default from "appsi_highs" to "highspy" and update the configuration, installation and deployment docs accordingly. Any Pyomo-based solver remains available as before. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- documentation/changelog.rst | 1 + documentation/configuration.rst | 14 ++++++++++---- documentation/host/deployment.rst | 5 +++-- documentation/host/installation.rst | 8 ++------ flexmeasures/utils/config_defaults.py | 2 +- 5 files changed, 17 insertions(+), 13 deletions(-) diff --git a/documentation/changelog.rst b/documentation/changelog.rst index 84b25dd8b1..5e58e61571 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -44,6 +44,7 @@ New features Infrastructure / Support ---------------------- +* Speed up scheduling jobs by building the scheduling problem directly with the HiGHS Python API (``highspy``), bypassing Pyomo's model construction and solution-ingestion overhead (roughly a second for a single-device job, and several seconds for multi-device jobs); this direct backend is the new default for the ``FLEXMEASURES_LP_SOLVER`` setting (``"highspy"``), while any Pyomo-based solver (e.g. the previous default ``"appsi_highs"``, or ``"cbc"``) remains available as before [see `PR #2365 `_] * Price fields in the flex-context (including nested commitment prices, which are now also held to the flex-context's shared currency) are selected for currency validation by field type (``PriceField``) instead of by name suffix [see `PR #2311 `_] * Document ``SECURITY_TWO_FACTOR`` and related 2FA configuration settings [see `PR #2340 `_] * ``flexmeasures db upgrade`` now runs ``VACUUM ANALYZE`` after upgrading by default, so Postgres has fresh planner statistics right after a migration; opt out with ``--no-vacuum`` [see `PR #2333 `_] diff --git a/documentation/configuration.rst b/documentation/configuration.rst index bfa87ccc6e..567f25fe61 100644 --- a/documentation/configuration.rst +++ b/documentation/configuration.rst @@ -55,11 +55,17 @@ Default: ``False`` FLEXMEASURES_LP_SOLVER ^^^^^^^^^^^^^^^^^^^^^^ -The command to run the scheduling solver. This is the executable command which FlexMeasures calls via the `pyomo library `_. Potential values might be ``cbc``, ``cplex``, ``glpk`` or ``appsi_highs``. Consult `their documentation `_ to learn more. -We have tested FlexMeasures with `HiGHS `_ and `Cbc `_. -Note that you need to install the solver, read more at :ref:`installing-a-solver`. +The scheduling solver backend. -Default: ``"appsi_highs"`` +The default, ``"highspy"``, builds the scheduling problem directly with the `HiGHS `_ Python API (``highspy``, which is installed with FlexMeasures). +This bypasses the `pyomo library `_ and is much faster to construct, while solving the exact same problem. + +Any other value is interpreted as the name of a Pyomo solver interface (the model is then built with Pyomo, which calls the solver). +Potential values might be ``cbc``, ``cplex``, ``glpk`` or ``appsi_highs``. Consult `the Pyomo documentation `_ to learn more. +We have tested FlexMeasures with `HiGHS `_ (both via ``highspy`` and via ``appsi_highs``) and `Cbc `_. +Note that for the Pyomo backends you need to install the solver yourself, read more at :ref:`installing-a-solver`. + +Default: ``"highspy"`` FLEXMEASURES_LP_SOLVER_OPTIONS diff --git a/documentation/host/deployment.rst b/documentation/host/deployment.rst index c5937d16ab..1610e5d6c8 100644 --- a/documentation/host/deployment.rst +++ b/documentation/host/deployment.rst @@ -89,11 +89,12 @@ Install the linear solver on the server --------------------------------------- To compute schedules, FlexMeasures uses the `HiGHS `_ mixed integer linear optimization solver (FlexMeasures solver by default) or `Cbc `_. -Solvers are used through `Pyomo `_\ , so in principle supporting a `different solver `_ would be possible. +By default, HiGHS is used directly through its Python API (``highspy``, which is installed together with FlexMeasures), so no extra installation is needed. +Solvers can also be used through `Pyomo `_\ , so in principle supporting a `different solver `_ would be possible. You tell FlexMeasures with the config setting :ref:`solver-config` which solver to use. -However, the solver also needs to be installed - in addition to FlexMeasures (the Docker image already has it). Here is advice on how to install the two solvers we test internally: +However, a Pyomo-based solver also needs to be installed - in addition to FlexMeasures (the Docker image already has it). Here is advice on how to install the two solvers we test internally: .. note:: We default to HiGHS, as it seems more powerful diff --git a/documentation/host/installation.rst b/documentation/host/installation.rst index 36073a8a8a..f1fbccc31e 100644 --- a/documentation/host/installation.rst +++ b/documentation/host/installation.rst @@ -360,13 +360,9 @@ Install an LP solver For computing schedules, the FlexMeasures platform uses a linear program solver. Currently that is the HiGHS or CBC solvers. -It's already installed in the Docker image. For yourself, you can simply install it like this: +The default solver (HiGHS, used directly via its Python API) is installed together with FlexMeasures, so there is nothing to do here. -.. code-block:: bash - - $ pip install highspy - -Read more on solvers (e.g. how to install a different one) at :ref:`installing-a-solver`. +Read more on solvers (e.g. how to install a different one, such as CBC) at :ref:`installing-a-solver`. diff --git a/flexmeasures/utils/config_defaults.py b/flexmeasures/utils/config_defaults.py index 114b26d2b8..d13c078815 100644 --- a/flexmeasures/utils/config_defaults.py +++ b/flexmeasures/utils/config_defaults.py @@ -154,7 +154,7 @@ class Config(object): "renewables": ["solar", "wind"], "EVSE": ["one-way_evse", "two-way_evse"], } # how to group assets by asset types - FLEXMEASURES_LP_SOLVER: str = "appsi_highs" + FLEXMEASURES_LP_SOLVER: str = "highspy" FLEXMEASURES_LP_SOLVER_OPTIONS: dict[str, str | int | float] = {} FLEXMEASURES_DEFAULT_JOB_TIMEOUT: timedelta = timedelta(seconds=180) FLEXMEASURES_JOB_TIMEOUT: dict[str, timedelta | str] = {} From 83cc6303bf9e506dfc91722e457fe41a6b911083 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:45:50 +0200 Subject: [PATCH 05/14] docs: point the changelog entry at the actual PR number Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- documentation/changelog.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/changelog.rst b/documentation/changelog.rst index 5e58e61571..7304767d9e 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -44,7 +44,7 @@ New features Infrastructure / Support ---------------------- -* Speed up scheduling jobs by building the scheduling problem directly with the HiGHS Python API (``highspy``), bypassing Pyomo's model construction and solution-ingestion overhead (roughly a second for a single-device job, and several seconds for multi-device jobs); this direct backend is the new default for the ``FLEXMEASURES_LP_SOLVER`` setting (``"highspy"``), while any Pyomo-based solver (e.g. the previous default ``"appsi_highs"``, or ``"cbc"``) remains available as before [see `PR #2365 `_] +* Speed up scheduling jobs by building the scheduling problem directly with the HiGHS Python API (``highspy``), bypassing Pyomo's model construction and solution-ingestion overhead (roughly a second for a single-device job, and several seconds for multi-device jobs); this direct backend is the new default for the ``FLEXMEASURES_LP_SOLVER`` setting (``"highspy"``), while any Pyomo-based solver (e.g. the previous default ``"appsi_highs"``, or ``"cbc"``) remains available as before [see `PR #2364 `_] * Price fields in the flex-context (including nested commitment prices, which are now also held to the flex-context's shared currency) are selected for currency validation by field type (``PriceField``) instead of by name suffix [see `PR #2311 `_] * Document ``SECURITY_TWO_FACTOR`` and related 2FA configuration settings [see `PR #2340 `_] * ``flexmeasures db upgrade`` now runs ``VACUUM ANALYZE`` after upgrading by default, so Postgres has fresh planner statistics right after a migration; opt out with ``--no-vacuum`` [see `PR #2333 `_] From 3555f8a899ccc45558d4afe5282fb182e8918aa3 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:53:46 +0200 Subject: [PATCH 06/14] perf: vectorize the subcommitment conversion Hoist convert_commitments_to_subcommitments to module level (it is solver-agnostic and closure-free) and share it between both scheduler backends, removing the duplicated copy from the highspy module. Splitting a commitment into per-group subcommitments now uses a single groupby pass (in order of first appearance, like pd.unique) instead of filtering the DataFrame once per group, and the price non-uniqueness checks are vectorized across all groups. This removes a cost that scaled quadratically with the number of time steps (each time step often forms its own group), benefiting both backends: a 2-device, 192-step benchmark drops from 0.47s to 0.28s via appsi_highs and from 0.37s to 0.12s via highspy. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- .../models/planning/highspy_optimization.py | 46 +------ .../models/planning/linear_optimization.py | 113 +++++++++--------- 2 files changed, 63 insertions(+), 96 deletions(-) diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 6668f90cd3..c505464132 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -41,7 +41,6 @@ from flexmeasures.data.models.planning import ( Commitment, - FlowCommitment, StockCommitment, ) from flexmeasures.data.models.planning.utils import initialize_series @@ -377,47 +376,10 @@ def device_scheduler_highspy( # noqa C901 % (resolution, resolution_c) ) - def convert_commitments_to_subcommitments( - dfs: list[pd.DataFrame], - ) -> tuple[list[pd.DataFrame], dict[int, int]]: - """Same transformation as in the Pyomo implementation (see reference).""" - commitment_mapping = {} - sub_commitments = [] - for c, df in enumerate(dfs): - # Make sure each commitment has "device" (default NaN) and "class" (default FlowCommitment) columns - if "device" not in df.columns: - df["device"] = np.nan - if "class" not in df.columns: - df["class"] = FlowCommitment - - df["j"] = range(len(df.index)) - groups = list(df["group"].unique()) - for group in groups: - sub_commitment = df[df["group"] == group].drop(columns=["group"]) - - # Catch non-uniqueness - if len(sub_commitment["upwards deviation price"].unique()) > 1: - raise ValueError( - "Commitment groups cannot have non-unique upwards deviation prices." - ) - if len(sub_commitment["downwards deviation price"].unique()) > 1: - raise ValueError( - "Commitment groups cannot have non-unique downwards deviation prices." - ) - if len(sub_commitment) == 1: - commitment_mapping[len(sub_commitments)] = c - sub_commitments.append(sub_commitment) - else: - down_commitment = sub_commitment.copy().drop( - columns="upwards deviation price" - ) - up_commitment = sub_commitment.copy().drop( - columns="downwards deviation price" - ) - commitment_mapping[len(sub_commitments)] = c - commitment_mapping[len(sub_commitments) + 1] = c - sub_commitments.extend([down_commitment, up_commitment]) - return sub_commitments, commitment_mapping + # This transformation is solver-agnostic and shared with the Pyomo backend. + from flexmeasures.data.models.planning.linear_optimization import ( + convert_commitments_to_subcommitments, + ) commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 58620d6e60..81e5000a09 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -73,6 +73,65 @@ def validate_highs_options(options: dict) -> None: ) +def convert_commitments_to_subcommitments( + dfs: list[pd.DataFrame], +) -> tuple[list[pd.DataFrame], dict[int, int]]: + """Transform commitments, each specifying a group for each time step, to sub-commitments, one per group. + + 'Groups' are a commitment concept (grouping time slots of a commitment), + making it possible that deviations/breaches can be accounted for properly within this group + (e.g. highest breach per calendar month defines the penalty). + Here, we define sub-commitments, by separating commitments by group and by direction of deviation (up, down). + + We also enumerate the time steps in a new column "j". + + For example, given contracts A and B (represented by 2 DataFrames), each with 3 groups, + we return (sub)commitments A1, A2, A3, B1, B2 and B3, + where A,B,C is the enumerated contract and 1,2,3 is the enumerated group. + + This helper is solver-agnostic and shared by both scheduler backends + (see the highspy_optimization module). + """ + commitment_mapping = {} + sub_commitments = [] + for c, df in enumerate(dfs): + # Make sure each commitment has "device" (default NaN) and "class" (default FlowCommitment) columns + if "device" not in df.columns: + df["device"] = np.nan + if "class" not in df.columns: + df["class"] = FlowCommitment + + df["j"] = range(len(df.index)) + + # Group rows by the "group" column in order of first appearance (like + # pd.unique), in a single pass rather than by filtering the DataFrame + # once per group (which would scale quadratically with the number of + # time steps, as each time step often forms its own group). + grouped = df.drop(columns=["group"]).groupby(df["group"], sort=False) + + # Catch non-uniqueness (vectorized across all groups) + if (grouped["upwards deviation price"].nunique(dropna=False) > 1).any(): + raise ValueError( + "Commitment groups cannot have non-unique upwards deviation prices." + ) + if (grouped["downwards deviation price"].nunique(dropna=False) > 1).any(): + raise ValueError( + "Commitment groups cannot have non-unique downwards deviation prices." + ) + + for _, sub_commitment in grouped: + if len(sub_commitment) == 1: + commitment_mapping[len(sub_commitments)] = c + sub_commitments.append(sub_commitment) + else: + down_commitment = sub_commitment.drop(columns="upwards deviation price") + up_commitment = sub_commitment.drop(columns="downwards deviation price") + commitment_mapping[len(sub_commitments)] = c + commitment_mapping[len(sub_commitments) + 1] = c + sub_commitments.extend([down_commitment, up_commitment]) + return sub_commitments, commitment_mapping + + def device_scheduler( # noqa C901 device_constraints: list[pd.DataFrame], ems_constraints: pd.DataFrame | list[pd.DataFrame], @@ -311,60 +370,6 @@ def device_scheduler( # noqa C901 % (resolution, resolution_c) ) - def convert_commitments_to_subcommitments( - dfs: list[pd.DataFrame], - ) -> tuple[list[pd.DataFrame], dict[int, int]]: - """Transform commitments, each specifying a group for each time step, to sub-commitments, one per group. - - 'Groups' are a commitment concept (grouping time slots of a commitment), - making it possible that deviations/breaches can be accounted for properly within this group - (e.g. highest breach per calendar month defines the penalty). - Here, we define sub-commitments, by separating commitments by group and by direction of deviation (up, down). - - We also enumerate the time steps in a new column "j". - - For example, given contracts A and B (represented by 2 DataFrames), each with 3 groups, - we return (sub)commitments A1, A2, A3, B1, B2 and B3, - where A,B,C is the enumerated contract and 1,2,3 is the enumerated group. - """ - commitment_mapping = {} - sub_commitments = [] - for c, df in enumerate(dfs): - # Make sure each commitment has "device" (default NaN) and "class" (default FlowCommitment) columns - if "device" not in df.columns: - df["device"] = np.nan - if "class" not in df.columns: - df["class"] = FlowCommitment - - df["j"] = range(len(df.index)) - groups = list(df["group"].unique()) - for group in groups: - sub_commitment = df[df["group"] == group].drop(columns=["group"]) - - # Catch non-uniqueness - if len(sub_commitment["upwards deviation price"].unique()) > 1: - raise ValueError( - "Commitment groups cannot have non-unique upwards deviation prices." - ) - if len(sub_commitment["downwards deviation price"].unique()) > 1: - raise ValueError( - "Commitment groups cannot have non-unique downwards deviation prices." - ) - if len(sub_commitment) == 1: - commitment_mapping[len(sub_commitments)] = c - sub_commitments.append(sub_commitment) - else: - down_commitment = sub_commitment.copy().drop( - columns="upwards deviation price" - ) - up_commitment = sub_commitment.copy().drop( - columns="downwards deviation price" - ) - commitment_mapping[len(sub_commitments)] = c - commitment_mapping[len(sub_commitments) + 1] = c - sub_commitments.extend([down_commitment, up_commitment]) - return sub_commitments, commitment_mapping - commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) device_group_lookup = {} From d566fc04fb35d48619b39e659878e5b18056a6fa Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 16:59:32 +0200 Subject: [PATCH 07/14] review: address Copilot comments on solver output and solver docs - Do not disable HiGHS output unconditionally in the direct backend; only force output_flag false when LOGGING_LEVEL is "INFO", mirroring exactly how the Pyomo path builds its solver options profile, so verbose logging modes can still see solver output. Operator-configured FLEXMEASURES_LP_SOLVER_OPTIONS are still applied last and can override. - Clarify in the configuration docs that a separate solver installation is only needed for external solvers such as cbc; both HiGHS-based choices (highspy and appsi_highs) rely on the bundled highspy package. - Remove the now-inconsistent "pip install highspy" instruction from the deployment docs, which already state that highspy ships with FlexMeasures. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- documentation/configuration.rst | 2 +- documentation/host/deployment.rst | 13 +++---------- .../data/models/planning/highspy_optimization.py | 1 - 3 files changed, 4 insertions(+), 12 deletions(-) diff --git a/documentation/configuration.rst b/documentation/configuration.rst index 567f25fe61..86b67da78b 100644 --- a/documentation/configuration.rst +++ b/documentation/configuration.rst @@ -63,7 +63,7 @@ This bypasses the `pyomo library `_ and is much faster to Any other value is interpreted as the name of a Pyomo solver interface (the model is then built with Pyomo, which calls the solver). Potential values might be ``cbc``, ``cplex``, ``glpk`` or ``appsi_highs``. Consult `the Pyomo documentation `_ to learn more. We have tested FlexMeasures with `HiGHS `_ (both via ``highspy`` and via ``appsi_highs``) and `Cbc `_. -Note that for the Pyomo backends you need to install the solver yourself, read more at :ref:`installing-a-solver`. +Note that a separate solver installation is only needed for external solvers such as ``cbc`` — both HiGHS-based choices (``highspy`` and ``appsi_highs``) rely on the ``highspy`` package that is installed together with FlexMeasures. Read more at :ref:`installing-a-solver`. Default: ``"highspy"`` diff --git a/documentation/host/deployment.rst b/documentation/host/deployment.rst index 1610e5d6c8..c766d03d13 100644 --- a/documentation/host/deployment.rst +++ b/documentation/host/deployment.rst @@ -94,21 +94,14 @@ Solvers can also be used through `Pyomo `_\ , so in princi You tell FlexMeasures with the config setting :ref:`solver-config` which solver to use. -However, a Pyomo-based solver also needs to be installed - in addition to FlexMeasures (the Docker image already has it). Here is advice on how to install the two solvers we test internally: - - .. note:: We default to HiGHS, as it seems more powerful -HiGHS can be installed using pip: - -.. code-block:: bash - - $ pip install highspy - +Both HiGHS-based solver choices (``highspy`` and ``appsi_highs``) rely on the ``highspy`` package, which is installed together with FlexMeasures — nothing more to do. More information on `the HiGHS website `_. -Cbc needs to be present on the server where FlexMeasures runs, under the ``cbc`` command. +An external solver, on the other hand, needs to be installed in addition to FlexMeasures (the Docker image already has it). +For example, Cbc needs to be present on the server where FlexMeasures runs, under the ``cbc`` command. You can install it on Debian like this: diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index c505464132..876b1a71ef 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -801,7 +801,6 @@ def _price_of(df: pd.DataFrame, column: str) -> float: # Build and solve the HiGHS model # --------------------------------------------------------------- h = highspy.Highs() - h.setOptionValue("output_flag", False) h.addVars(ncol, lower, upper) h.changeColsCost(ncol, np.arange(ncol, dtype=np.int32), cost) From f59b05694db7f4564c6aea2063e194fcacd2185f Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Tue, 28 Jul 2026 17:15:36 +0200 Subject: [PATCH 08/14] review: remove dead commodity_devices lookup from the direct backend The commodity -> device indices lookup was only consumed by the Pyomo path's ems_flow_commitment_equalities, which the direct backend deliberately does not build (they are free rows). A breadcrumb comment keeps pointing readers at that documented skip. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WtuVTVfL4fQ9QSqbLXmAGD Signed-off-by: F.N. Claessen --- .../data/models/planning/highspy_optimization.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 876b1a71ef..eeb4fcdb62 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -349,16 +349,9 @@ def device_scheduler_highspy( # noqa C901 df["group"] = group commitments.append(df) - # commodity -> set(device indices) - commodity_devices: dict = {} - for df in commitments: - if "commodity" not in df.columns or "device" not in df.columns: - continue - for _, row in df[["commodity", "device"]].dropna().iterrows(): - devices = row["device"] - if not isinstance(devices, (list, tuple, set)): - devices = [devices] - commodity_devices.setdefault(row["commodity"], set()).update(devices) + # NB the Pyomo implementation builds a commodity -> device indices lookup here, + # but it is only consumed by its ems_flow_commitment_equalities, which this + # backend deliberately does not build (they are free rows; see module docstring). # Check if commitments have the same time window and resolution as the constraints for commitment in commitments: From 305e037414790e851cbdb62886a4e0cfd637fac3 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Mon, 3 Aug 2026 15:09:12 +0200 Subject: [PATCH 09/14] Forward device_scheduler arguments to the highspy backend by name The dispatch to the direct HiGHS backend listed its keyword arguments by hand. That is a trap for the branches currently adding scheduling parameters: whoever adds the next one (coupling_groups in #2218, balance_groups in #2289) works on the Pyomo model further down the file, and a parameter missing from the dispatch list would not fail. It would simply never reach the backend, producing a schedule computed as if the constraint had never been requested -- and since this PR makes "highspy" the default solver, that would be silently wrong. Forward by name instead, mapping device_scheduler's signature onto the backend's. An argument the backend does not model raises NotImplementedError naming it, but only when the caller actually set it, so leaving a future parameter at its default stays free. The signature comparison is cached on the two function objects (~70 us once per process, 2 us per call after). Also record, in the highspy module docstring, that the deliberate omission of ems_flow_commitment_equalities stops being harmless once #2355 gives those rows bounds, and that #2380 routes unscoped flex-context commitments through grouped_commitment_equalities (which this backend does build). Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01X2jLjsQzwztDc7Nxz1ozmQ Signed-off-by: F.N. Claessen --- .../models/planning/highspy_optimization.py | 24 ++++ .../models/planning/linear_optimization.py | 104 ++++++++++++++++-- .../tests/test_highspy_equivalence.py | 50 +++++++++ 3 files changed, 166 insertions(+), 12 deletions(-) diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index eeb4fcdb62..4da7e37a8f 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -22,6 +22,29 @@ constraint family returns ``(None, expr, None)``, i.e. a constraint without bounds, which ends up as a free (vacuous) row in HiGHS. We skip building the free rows altogether. + + .. note:: This deviation is only harmless for as long as those rows stay free, + and two branches in flight change that: + + - #2355 gives the constraint the same one-sided bounds that + ``grouped_commitment_equalities`` already uses, making it bind. Once that + lands, this module must build the row too — otherwise an EMS-level + commitment is silently ignored under this backend. The row is then + structurally identical to the grouped one, differing only in the set it + sums over (all devices, or the commodity's devices). + - #2380 makes an unscoped flex-context commitment device-grouped + (``device=``, ``device_group=``) + rather than per-device, so commitments coming from ``convert_to_commitments`` + route through ``grouped_commitment_equalities``, which this module does + build. That narrows the exposure above to callers constructing an + EMS-level ``FlowCommitment`` (``device=None``) directly, but does not + remove it. + + Note that ``tests/test_commitments.py`` does not use the + ``app_with_each_solver`` fixture, so its cases only ever run under the + configured default solver -- which ``config_defaults`` now sets to + ``"highspy"``. Coverage of this constraint family belongs in + ``tests/test_highspy_equivalence.py`` to run under both backends. - Rows whose computed bounds are impossible to satisfy for any finite value (upper bound of -inf, or lower bound of +inf, as happens when a commitment quantity is +/-inf) are skipped. On the Pyomo path such rows are rejected by @@ -352,6 +375,7 @@ def device_scheduler_highspy( # noqa C901 # NB the Pyomo implementation builds a commodity -> device indices lookup here, # but it is only consumed by its ems_flow_commitment_equalities, which this # backend deliberately does not build (they are free rows; see module docstring). + # That lookup has to come back when #2355 makes those rows bind. # Check if commitments have the same time window and resolution as the constraints for commitment in commitments: diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 81e5000a09..61d972e64b 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -1,6 +1,8 @@ from __future__ import annotations +import inspect import math +from functools import lru_cache from flask import current_app import pandas as pd @@ -132,6 +134,91 @@ def convert_commitments_to_subcommitments( return sub_commitments, commitment_mapping +def _left_at_default(value, default) -> bool: + """Whether an argument was left at its default. + + Best-effort: pandas values compare element-wise, so ``value == default`` may + return an array (or raise) rather than a bool. Anything we cannot decide is + reported as "not the default", which errs towards raising in + :func:`_arguments_for_highspy_backend` rather than silently dropping a value. + """ + if value is default: + return True + if default is inspect.Parameter.empty: + return False + try: + return bool(value == default) + except (TypeError, ValueError): + return False + + +@lru_cache(maxsize=None) +def _backend_argument_map(declared_by, supported_by) -> tuple[frozenset, tuple]: + """Which of ``declared_by``'s arguments ``supported_by`` accepts, and which it lacks. + + Signatures are static, so this is computed once per process (roughly 70 us, + which is not worth paying on every schedule). Caching on the two function + objects rather than on nothing keeps it correct when a test substitutes one + of them. + """ + declared = inspect.signature(declared_by).parameters + supported = frozenset(inspect.signature(supported_by).parameters) + missing = tuple( + (name, parameter.default) + for name, parameter in declared.items() + if name not in supported + ) + return frozenset(declared) & supported, missing + + +def _arguments_for_highspy_backend(passed_arguments: dict) -> dict: + """Map ``device_scheduler``'s arguments onto the direct HiGHS backend's signature. + + A hand-written keyword list here would be a trap: whoever adds the next + ``device_scheduler`` parameter naturally works on the Pyomo model further + down this file, and a parameter missing from that list would not fail — it + would simply never reach the backend. Under ``FLEXMEASURES_LP_SOLVER="highspy"`` + (the default) that yields a schedule computed as if the constraint had never + been requested: plausible-looking, silently wrong, and not caught by tests + written before the default was flipped. + + Forwarding by name removes that failure mode entirely. The remaining case — + an argument the direct backend does not model at all — is caught statically + by ``test_every_device_scheduler_argument_currently_reaches_the_backend``, so + it cannot reach a release. The raise below is only a backstop for a build + where that test did not run; it costs nothing while the signatures agree. + + This is a live concern rather than a hypothetical one: ``device_scheduler`` + is gaining ``coupling_groups`` (#2218) and ``balance_groups`` (#2289) on + branches in flight, and each needs explicit support here. + """ + from flexmeasures.data.models.planning.highspy_optimization import ( + device_scheduler_highspy, + ) + + forwardable, missing = _backend_argument_map( + device_scheduler, device_scheduler_highspy + ) + if missing: + in_use = sorted( + name + for name, default in missing + if not _left_at_default(passed_arguments[name], default) + ) + if in_use: + raise NotImplementedError( + "The direct HiGHS backend (FLEXMEASURES_LP_SOLVER='highspy') does not" + f" model these device_scheduler arguments: {', '.join(in_use)}." + " Add support for them in" + " flexmeasures.data.models.planning.highspy_optimization (and extend" + " tests/test_highspy_equivalence.py), or configure a Pyomo-based" + " solver such as 'appsi_highs'." + ) + return { + name: value for name, value in passed_arguments.items() if name in forwardable + } + + def device_scheduler( # noqa C901 device_constraints: list[pd.DataFrame], ems_constraints: pd.DataFrame | list[pd.DataFrame], @@ -207,22 +294,15 @@ def device_scheduler( # noqa C901 # See the highspy_optimization module, which mirrors the model built below # and returns compatible result objects. if current_app.config.get("FLEXMEASURES_LP_SOLVER") == "highspy": + # Arguments are forwarded by name, and an argument the direct backend cannot + # model raises rather than being dropped. See _arguments_for_highspy_backend. + highspy_arguments = _arguments_for_highspy_backend(locals()) + from flexmeasures.data.models.planning.highspy_optimization import ( device_scheduler_highspy, ) - return device_scheduler_highspy( - device_constraints=device_constraints, - ems_constraints=ems_constraints, - commitment_quantities=commitment_quantities, - commitment_downwards_deviation_price=commitment_downwards_deviation_price, - commitment_upwards_deviation_price=commitment_upwards_deviation_price, - commitments=commitments, - initial_stock=initial_stock, - stock_groups=stock_groups, - ems_constraint_groups=ems_constraint_groups, - device_power_bands=device_power_bands, - ) + return device_scheduler_highspy(**highspy_arguments) model = ConcreteModel() diff --git a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py index f0452ab48a..48d4307823 100644 --- a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py +++ b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py @@ -12,6 +12,7 @@ from __future__ import annotations +import inspect from datetime import timedelta import numpy as np @@ -19,6 +20,10 @@ import pytest from flexmeasures.data.models.planning import FlowCommitment, StockCommitment +from flexmeasures.data.models.planning import linear_optimization +from flexmeasures.data.models.planning.highspy_optimization import ( + device_scheduler_highspy, +) from flexmeasures.data.models.planning.linear_optimization import device_scheduler from flexmeasures.data.models.planning.utils import initialize_df @@ -307,3 +312,48 @@ def test_highspy_matches_pyomo_when_infeasible(app): # Mirrored fallback behavior: no costs (all variables at zero) assert costs_p == costs_h == 0 + + +def test_unsupported_argument_is_rejected_not_ignored(): + """A device_scheduler argument the direct backend cannot model must raise. + + ``device_scheduler`` forwards its arguments to the direct HiGHS backend by + name. Whoever adds the next scheduling parameter (``coupling_groups`` in + #2218, ``balance_groups`` in #2289) works on the Pyomo model, and a + parameter that never reached the backend would not fail -- it would produce + a schedule computed as if the constraint had never been requested. Since + ``highspy`` is the default solver, that would be silently wrong. + """ + real_device_scheduler = linear_optimization.device_scheduler + + def device_scheduler_of_the_future( + device_constraints, + ems_constraints, + future_parameter=None, + ): + """Stand-in for a device_scheduler that grew a parameter.""" + + linear_optimization.device_scheduler = device_scheduler_of_the_future + try: + arguments = dict( + device_constraints=[], ems_constraints=None, future_parameter=None + ) + + # Left at its default, the new parameter costs existing callers nothing. + forwarded = linear_optimization._arguments_for_highspy_backend(arguments) + assert "future_parameter" not in forwarded + assert set(forwarded) == {"device_constraints", "ems_constraints"} + + # Actually set, it must be reported rather than dropped. + arguments["future_parameter"] = {"some group": [(0, 1.0)]} + with pytest.raises(NotImplementedError, match="future_parameter"): + linear_optimization._arguments_for_highspy_backend(arguments) + finally: + linear_optimization.device_scheduler = real_device_scheduler + + +def test_every_device_scheduler_argument_currently_reaches_the_backend(): + """The direct backend supports the whole current device_scheduler signature.""" + declared = set(inspect.signature(device_scheduler).parameters) + supported = set(inspect.signature(device_scheduler_highspy).parameters) + assert declared <= supported, declared - supported From 13ae603c4b5157b01fbe504f9521f1732eff29f3 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Mon, 3 Aug 2026 15:10:39 +0200 Subject: [PATCH 10/14] refactor: share the scheduler's solver-agnostic input handling The two backends build the same model in two representations, so the model construction is necessarily written twice. Everything around it was too: 199 lines were byte-identical between linear_optimization.py and highspy_optimization.py -- argument normalisation, stock groups and their validation, legacy commitment conversion, the sub-commitment split, device_group_lookup, the convex-curve check, the Big-Ms, band validation, the HiGHS option profile, and the cost/schedule assembly. None of it has a solver in it, and keeping it twice meant the two paths could drift apart on input handling, which the equivalence tests are not aimed at. Move it to a new scheduling_problem module: prepare_scheduling_problem() returns a SchedulingProblem that both backends unpack, plus solver_options() and the result-assembly helpers. Duplication between the backends drops from 199 to 40 lines, and those 40 are the shared signature and the call itself. Two deliberate changes while moving: - commodity_devices becomes a cached_property. It is a per-row scan that only the Pyomo ems_flow_commitment_equalities needs, so computing it eagerly would put a real cost on the fast path. Pyomo pays what it did before; the direct backend pays nothing until it needs it (see #2355). - initial_stock_of() casts its index to int, as the direct backend already did. The Pyomo version raised TypeError on the numpy float device indices a commitment's "device" column can carry. The empty-commitments case also stops raising on pd.concat([]), which previously made the Pyomo path crash where the direct path coped. Verified: 265 passed, 3 xfailed across the planning suite under all three solver parameters, plus the scheduling-job and API schedule tests. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01X2jLjsQzwztDc7Nxz1ozmQ Signed-off-by: F.N. Claessen --- .../models/planning/highspy_optimization.py | 320 ++-------- .../models/planning/linear_optimization.py | 464 ++------------- .../models/planning/scheduling_problem.py | 562 ++++++++++++++++++ 3 files changed, 657 insertions(+), 689 deletions(-) create mode 100644 flexmeasures/data/models/planning/scheduling_problem.py diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 4da7e37a8f..50b258000c 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -59,14 +59,18 @@ import numpy as np import pandas as pd -from flask import current_app -from pandas.tseries.frequencies import to_offset from flexmeasures.data.models.planning import ( Commitment, StockCommitment, ) -from flexmeasures.data.models.planning.utils import initialize_series +from flexmeasures.data.models.planning.scheduling_problem import ( + aggregate_commodity_costs, + aggregate_subcommitment_costs, + planned_power_per_device, + prepare_scheduling_problem, + solver_options, +) infinity = float("inf") @@ -267,224 +271,34 @@ def device_scheduler_highspy( # noqa C901 if len(device_constraints) == 0: return [], 0, HighspySolverResults("unknown", "ok"), model - # Get timing from first device - start = device_constraints[0].index.to_pydatetime()[0] - # Workaround for https://github.com/pandas-dev/pandas/issues/53643. Was: resolution = pd.to_timedelta(device_constraints[0].index.freq) - resolution = pd.to_timedelta(device_constraints[0].index.freq).to_pytimedelta() - end = device_constraints[0].index.to_pydatetime()[-1] + resolution - - # Normalise EMS constraints to a list of (DataFrame, device-group) pairs. - all_devices = list(range(len(device_constraints))) - if isinstance(ems_constraints, pd.DataFrame): - ems_constraints_list = [ems_constraints] - ems_constraint_device_groups = [all_devices] - else: - ems_constraints_list = ems_constraints - if ems_constraint_groups is None: - if len(ems_constraints_list) > 1: - raise ValueError( - "When passing multiple EMS constraint DataFrames, you must also specify ems_constraint_groups." - ) - ems_constraint_device_groups = [all_devices for _ in ems_constraints_list] - else: - ems_constraint_device_groups = ems_constraint_groups - - # 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 keys are namespaced strings, as in the Pyomo implementation. - if stock_groups: - for g, devices in stock_groups.items(): - for d in devices: - 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[str, list[int]] = {} - for d, g in device_to_group.items(): - group_to_devices.setdefault(g, []).append(d) - - # Devices sharing a stock may not declare different storage efficiencies - # or initial stocks (the stock recursion is modelled once per stock group). - for g, group_devices in group_to_devices.items(): - if len(group_devices) > 1: - 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 = [] - else: - commitments = [ - c.to_frame() if isinstance(c, Commitment) else c for c in commitments - ] - if commitment_quantities is not None: - from flexmeasures.data.models.planning.utils import initialize_df - - for quantity, down, up in zip( - commitment_quantities, - commitment_downwards_deviation_price, - commitment_upwards_deviation_price, - ): - # Turn prices per commitment into prices per commitment flow - if all(isinstance(price, float) for price in down) or isinstance( - down, float - ): - down = initialize_series(down, start, end, resolution) - if all(isinstance(price, float) for price in up) or isinstance(up, float): - up = initialize_series(up, start, end, resolution) - - group = initialize_series(list(range(len(down))), start, end, resolution) - df = initialize_df( - ["quantity", "downwards deviation price", "upwards deviation price"], - start, - end, - resolution, - ) - df["quantity"] = quantity - df["downwards deviation price"] = down - df["upwards deviation price"] = up - df["group"] = group - commitments.append(df) - - # NB the Pyomo implementation builds a commodity -> device indices lookup here, - # but it is only consumed by its ems_flow_commitment_equalities, which this - # backend deliberately does not build (they are free rows; see module docstring). - # That lookup has to come back when #2355 makes those rows bind. - - # Check if commitments have the same time window and resolution as the constraints - for commitment in commitments: - start_c = commitment.index.to_pydatetime()[0] - resolution_c = pd.to_timedelta(commitment.index.freq) - end_c = commitment.index.to_pydatetime()[-1] + resolution - if not (start_c == start and end_c == end): - raise Exception( - "Not implemented for different time windows.\n(%s,%s)\n(%s,%s)" - % (start, end, start_c, end_c) - ) - if resolution_c != resolution: - raise Exception( - "Not implemented for different resolutions.\n%s\n%s" - % (resolution, resolution_c) - ) - - # This transformation is solver-agnostic and shared with the Pyomo backend. - from flexmeasures.data.models.planning.linear_optimization import ( - convert_commitments_to_subcommitments, + problem = prepare_scheduling_problem( + device_constraints=device_constraints, + ems_constraints=ems_constraints, + commitment_quantities=commitment_quantities, + commitment_downwards_deviation_price=commitment_downwards_deviation_price, + commitment_upwards_deviation_price=commitment_upwards_deviation_price, + commitments=commitments, + initial_stock=initial_stock, + stock_groups=stock_groups, + ems_constraint_groups=ems_constraint_groups, + device_power_bands=device_power_bands, ) - commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) - - device_group_lookup: dict[int, dict] = {} - for c, df in enumerate(commitments): - # Stock-scoped commitments couple to their stock group as a whole. - 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. - 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. - rows = df[["device"]].dropna() - - device_group_lookup[c] = {} - - for _, row in rows.iterrows(): - d = row["device"] - 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 (mirrors the Pyomo path) - if commitments: - df = pd.concat(commitments)[ - ["upwards deviation price", "downwards deviation price"] - ] - df = df.groupby(level=0).sum() - convex_cost_curve = ( - len(df[df["upwards deviation price"] < df["downwards deviation price"]]) - == 0 - ) - else: - convex_cost_curve = True - - bigM_columns = ["derivative max", "derivative min", "derivative equals"] - # Compute a good value for our Big-Ms - 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]) - - # Both Md and Mc have to be 1 MW, at least - Md = max(Md, 1) - Mc = max(Mc, 1) - - for d in range(len(device_constraints)): - if "stock delta" not in device_constraints[d].columns: - device_constraints[d]["stock delta"] = 0 - else: - device_constraints[d]["stock delta"] = ( - device_constraints[d]["stock delta"].astype(float).fillna(0) - ) - - # Look up power bands (S2 operation modes) per device - if device_power_bands is None: - device_power_bands = [None] * len(device_constraints) - elif len(device_power_bands) != len(device_constraints): - raise ValueError( - f"device_power_bands lists {len(device_power_bands)} devices, " - f"while device_constraints lists {len(device_constraints)} devices." - ) - band_lookup: dict[int, list[tuple[float, float]]] = { - d: list(bands) - for d, bands in enumerate(device_power_bands) - if bands is not None and len(bands) > 0 - } - for d, bands in band_lookup.items(): - for band in bands: - if len(band) != 2 or band[0] > band[1]: - raise ValueError( - f"Invalid power band {band} for device {d}: " - f"expected a (min, max) pair with min <= max." - ) + # Local aliases, so that the model below reads as it did before the (solver-agnostic) + # input handling moved to the scheduling_problem module. + start, end, resolution = problem.start, problem.end, problem.resolution + device_constraints = problem.device_constraints + ems_constraints_list = problem.ems_constraints_list + ems_constraint_device_groups = problem.ems_constraint_device_groups + device_to_group = problem.device_to_group + group_to_devices = problem.group_to_devices + commitments = problem.commitments + commitment_mapping = problem.commitment_mapping + device_group_lookup = problem.device_group_lookup + convex_cost_curve = problem.convex_cost_curve + Md, Mc = problem.Md, problem.Mc + band_lookup = problem.band_lookup + _initial_stock_of = problem.initial_stock_of # --------------------------------------------------------------- # Numeric model data (vectorized versions of the Pyomo Param rules) @@ -533,12 +347,6 @@ def device_scheduler_highspy( # noqa C901 up_eff[d] = _column_or_default(dc, "derivative up efficiency", 1) delta[d] = _column(dc, "stock delta") - def _initial_stock_of(d) -> float: - if isinstance(initial_stock, list): - # No initial stock defined for inflexible device - return initial_stock[int(d)] if d < len(initial_stock) else 0 - return initial_stock - # --------------------------------------------------------------- # Column (variable) layout # --------------------------------------------------------------- @@ -842,29 +650,9 @@ def _price_of(df: pd.DataFrame, column: str) -> float: if nrow > 0: h.addRows(nrow, row_lower, row_upper, nnz, row_starts, a_index, a_value) - # Apply the same solver options as the Pyomo path applies for HiGHS solvers - profile = { - "mip_rel_gap": "0", - "mip_abs_gap": "0", - "primal_feasibility_tolerance": "1e-9", - "dual_feasibility_tolerance": "1e-9", - "mip_feasibility_tolerance": "1e-9", - } - # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO - if current_app.config["LOGGING_LEVEL"] == "INFO": - profile["output_flag"] = "false" - - # Apply operator-configured options last, so they override the defaults above. - configured_options = current_app.config.get("FLEXMEASURES_LP_SOLVER_OPTIONS") or {} - if configured_options: - from flexmeasures.data.models.planning.linear_optimization import ( - validate_highs_options, - ) - - validate_highs_options(configured_options) - profile.update(configured_options) - - for option_name, option_value in profile.items(): + # The same options the Pyomo path applies for HiGHS solvers ("highspy" matches + # on "highs"), so the two backends cannot disagree on tolerances. + for option_name, option_value in solver_options("highspy").items(): h.setOptionValue(option_name, option_value) h.run() @@ -910,36 +698,12 @@ def _price_of(df: pd.DataFrame, column: str) -> float: for c in range(C): planned_costs += subcommitment_costs[c] - # Map subcommitment costs to commitments - commitment_costs: dict = {} - for g, v in subcommitment_costs.items(): - c = commitment_mapping[g] - commitment_costs[c] = commitment_costs.get(c, 0) + v + planned_power = planned_power_per_device(ems_values, start, end, resolution) - planned_power_per_device = [] - for d in range(D): - planned_power_per_device.append( - initialize_series( - data=list(ems_values[d]), - start=start, - end=end, - resolution=to_offset(resolution), - ) - ) - - commodity_costs: dict = {} - for c in range(C): - commodity = None - if "commodity" in commitments[c].columns: - commodity = commitments[c]["commodity"].iloc[0] - if commodity is None or (isinstance(commodity, float) and np.isnan(commodity)): - continue - commodity_costs[commodity] = ( - commodity_costs.get(commodity, 0) + subcommitment_costs[c] - ) - - model.commitment_costs = commitment_costs - model.commodity_costs = commodity_costs + model.commitment_costs = aggregate_subcommitment_costs( + subcommitment_costs, commitment_mapping + ) + model.commodity_costs = aggregate_commodity_costs(commitments, subcommitment_costs) model.costs = planned_costs model.d = range(D) model.j = range(T) @@ -956,4 +720,4 @@ def _price_of(df: pd.DataFrame, column: str) -> float: {(d, j): float(sign_values[d, j]) for d in range(D) for j in range(T)} ) - return planned_power_per_device, planned_costs, results, model + return planned_power, planned_costs, results, model diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 61d972e64b..6b5a3b22cc 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -1,13 +1,11 @@ from __future__ import annotations import inspect -import math from functools import lru_cache from flask import current_app import pandas as pd import numpy as np -from pandas.tseries.frequencies import to_offset from pyomo.core import ( ConcreteModel, Var, @@ -31,109 +29,20 @@ FlowCommitment, StockCommitment, ) -from flexmeasures.data.models.planning.utils import initialize_series, initialize_df +from flexmeasures.data.models.planning.scheduling_problem import ( # noqa F401 + aggregate_commodity_costs, + aggregate_subcommitment_costs, + convert_commitments_to_subcommitments, + loss_coefficients, + planned_power_per_device, + prepare_scheduling_problem, + solver_options, + validate_highs_options, +) infinity = float("inf") -def validate_highs_options(options: dict) -> None: - """Raise if HiGHS would refuse any of these options. - - Pyomo's appsi_highs interface applies solver options without checking HiGHS' - return status, so an unknown name, an invalid value, or a feature missing from - the installed HiGHS build is otherwise ignored without a word. That silently - turns a mis-typed option into a no-op, and a benchmark of it into a false - negative. Probing a throwaway Highs instance surfaces the rejection instead. - """ - try: - import highspy - except ImportError: - # Solver named "*highs*" but highspy absent: let the solver interface complain. - return - - probe = highspy.Highs() - probe.setOptionValue("output_flag", False) - rejected = [ - f"{name}={value!r}" - for name, value in options.items() - if probe.setOptionValue(name, value) != highspy.HighsStatus.kOk - ] - if rejected: - raise ValueError( - f"HiGHS rejected these FLEXMEASURES_LP_SOLVER_OPTIONS: {', '.join(rejected)}." - " The option name may be unknown, the value invalid, or the feature absent" - " from this HiGHS build. For example, the HiPO solver (solver='hipo') needs" - " a HiGHS built against BLAS and METIS, which the pip-installed highspy is not." - ) - - if "threads" in options or "parallel" in options: - current_app.logger.warning( - "FLEXMEASURES_LP_SOLVER_OPTIONS sets 'threads' and/or 'parallel'. HiGHS" - " initializes its thread scheduler once per process, so inside a long-lived" - " worker only the first solve honours these; later solves fail with 'global" - " scheduler has already been initialized' and yield no schedule." - ) - - -def convert_commitments_to_subcommitments( - dfs: list[pd.DataFrame], -) -> tuple[list[pd.DataFrame], dict[int, int]]: - """Transform commitments, each specifying a group for each time step, to sub-commitments, one per group. - - 'Groups' are a commitment concept (grouping time slots of a commitment), - making it possible that deviations/breaches can be accounted for properly within this group - (e.g. highest breach per calendar month defines the penalty). - Here, we define sub-commitments, by separating commitments by group and by direction of deviation (up, down). - - We also enumerate the time steps in a new column "j". - - For example, given contracts A and B (represented by 2 DataFrames), each with 3 groups, - we return (sub)commitments A1, A2, A3, B1, B2 and B3, - where A,B,C is the enumerated contract and 1,2,3 is the enumerated group. - - This helper is solver-agnostic and shared by both scheduler backends - (see the highspy_optimization module). - """ - commitment_mapping = {} - sub_commitments = [] - for c, df in enumerate(dfs): - # Make sure each commitment has "device" (default NaN) and "class" (default FlowCommitment) columns - if "device" not in df.columns: - df["device"] = np.nan - if "class" not in df.columns: - df["class"] = FlowCommitment - - df["j"] = range(len(df.index)) - - # Group rows by the "group" column in order of first appearance (like - # pd.unique), in a single pass rather than by filtering the DataFrame - # once per group (which would scale quadratically with the number of - # time steps, as each time step often forms its own group). - grouped = df.drop(columns=["group"]).groupby(df["group"], sort=False) - - # Catch non-uniqueness (vectorized across all groups) - if (grouped["upwards deviation price"].nunique(dropna=False) > 1).any(): - raise ValueError( - "Commitment groups cannot have non-unique upwards deviation prices." - ) - if (grouped["downwards deviation price"].nunique(dropna=False) > 1).any(): - raise ValueError( - "Commitment groups cannot have non-unique downwards deviation prices." - ) - - for _, sub_commitment in grouped: - if len(sub_commitment) == 1: - commitment_mapping[len(sub_commitments)] = c - sub_commitments.append(sub_commitment) - else: - down_commitment = sub_commitment.drop(columns="upwards deviation price") - up_commitment = sub_commitment.drop(columns="downwards deviation price") - commitment_mapping[len(sub_commitments)] = c - commitment_mapping[len(sub_commitments) + 1] = c - sub_commitments.extend([down_commitment, up_commitment]) - return sub_commitments, commitment_mapping - - def _left_at_default(value, default) -> bool: """Whether an argument was left at its default. @@ -310,241 +219,34 @@ def device_scheduler( # noqa C901 if len(device_constraints) == 0: return [], 0, SolverResults(), model - # Get timing from first device - start = device_constraints[0].index.to_pydatetime()[0] - # Workaround for https://github.com/pandas-dev/pandas/issues/53643. Was: resolution = pd.to_timedelta(device_constraints[0].index.freq) - resolution = pd.to_timedelta(device_constraints[0].index.freq).to_pytimedelta() - end = device_constraints[0].index.to_pydatetime()[-1] + resolution - - # Normalise EMS constraints to a list of (DataFrame, device-group) pairs. - # A single DataFrame (legacy behaviour) applies to the summed flow of all devices; - # a list of DataFrames applies one EMS-level constraint per device group, as set up - # per commodity by the StorageScheduler. - all_devices = list(range(len(device_constraints))) - if isinstance(ems_constraints, pd.DataFrame): - ems_constraints_list = [ems_constraints] - ems_constraint_device_groups = [all_devices] - else: - ems_constraints_list = ems_constraints - if ems_constraint_groups is None: - if len(ems_constraints_list) > 1: - raise ValueError( - "When passing multiple EMS constraint DataFrames, you must also specify ems_constraint_groups." - ) - ems_constraint_device_groups = [all_devices for _ in ems_constraints_list] - else: - ems_constraint_device_groups = ems_constraint_groups - - # 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 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] = 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 = [] - else: - commitments = [ - c.to_frame() if isinstance(c, Commitment) else c for c in commitments - ] - if commitment_quantities is not None: - for quantity, down, up in zip( - commitment_quantities, - commitment_downwards_deviation_price, - commitment_upwards_deviation_price, - ): - - # Turn prices per commitment into prices per commitment flow - if all(isinstance(price, float) for price in down) or isinstance( - down, float - ): - down = initialize_series(down, start, end, resolution) - if all(isinstance(price, float) for price in up) or isinstance(up, float): - up = initialize_series(up, start, end, resolution) - - group = initialize_series(list(range(len(down))), start, end, resolution) - df = initialize_df( - ["quantity", "downwards deviation price", "upwards deviation price"], - start, - end, - resolution, - ) - df["quantity"] = quantity - df["downwards deviation price"] = down - df["upwards deviation price"] = up - df["group"] = group - commitments.append(df) - - # commodity → set(device indices) - commodity_devices = {} - - for df in commitments: - if "commodity" not in df.columns or "device" not in df.columns: - continue - - for _, row in df[["commodity", "device"]].dropna().iterrows(): - devices = row["device"] - if not isinstance(devices, (list, tuple, set)): - devices = [devices] - - commodity_devices.setdefault(row["commodity"], set()).update(devices) - - # Check if commitments have the same time window and resolution as the constraints - for commitment in commitments: - start_c = commitment.index.to_pydatetime()[0] - resolution_c = pd.to_timedelta(commitment.index.freq) - end_c = commitment.index.to_pydatetime()[-1] + resolution - if not (start_c == start and end_c == end): - raise Exception( - "Not implemented for different time windows.\n(%s,%s)\n(%s,%s)" - % (start, end, start_c, end_c) - ) - if resolution_c != resolution: - raise Exception( - "Not implemented for different resolutions.\n%s\n%s" - % (resolution, resolution_c) - ) - - commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) - - 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. - 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"] - ] - df = df.groupby(level=0).sum() - if len(df[df["upwards deviation price"] < df["downwards deviation price"]]) == 0: - convex_cost_curve = True - else: - convex_cost_curve = False - - 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]) - - # Both Md and Mc have to be 1 MW, at least - Md = max(Md, 1) - Mc = max(Mc, 1) - - for d in range(len(device_constraints)): - if "stock delta" not in device_constraints[d].columns: - device_constraints[d]["stock delta"] = 0 - else: - device_constraints[d]["stock delta"] = ( - device_constraints[d]["stock delta"].astype(float).fillna(0) - ) + problem = prepare_scheduling_problem( + device_constraints=device_constraints, + ems_constraints=ems_constraints, + commitment_quantities=commitment_quantities, + commitment_downwards_deviation_price=commitment_downwards_deviation_price, + commitment_upwards_deviation_price=commitment_upwards_deviation_price, + commitments=commitments, + initial_stock=initial_stock, + stock_groups=stock_groups, + ems_constraint_groups=ems_constraint_groups, + device_power_bands=device_power_bands, + ) - # Look up power bands (S2 operation modes) per device - if device_power_bands is None: - device_power_bands = [None] * len(device_constraints) - elif len(device_power_bands) != len(device_constraints): - raise ValueError( - f"device_power_bands lists {len(device_power_bands)} devices, " - f"while device_constraints lists {len(device_constraints)} devices." - ) - band_lookup: dict[int, list[tuple[float, float]]] = { - d: list(bands) - for d, bands in enumerate(device_power_bands) - if bands is not None and len(bands) > 0 - } - for d, bands in band_lookup.items(): - for band in bands: - if len(band) != 2 or band[0] > band[1]: - raise ValueError( - f"Invalid power band {band} for device {d}: " - f"expected a (min, max) pair with min <= max." - ) + # Local aliases, so that the model below reads as it did before the (solver-agnostic) + # input handling moved to the scheduling_problem module. + start, end, resolution = problem.start, problem.end, problem.resolution + device_constraints = problem.device_constraints + ems_constraints_list = problem.ems_constraints_list + ems_constraint_device_groups = problem.ems_constraint_device_groups + device_to_group = problem.device_to_group + group_to_devices = problem.group_to_devices + commitments = problem.commitments + commitment_mapping = problem.commitment_mapping + device_group_lookup = problem.device_group_lookup + convex_cost_curve = problem.convex_cost_curve + Md, Mc = problem.Md, problem.Mc + band_lookup = problem.band_lookup + _initial_stock_of = problem.initial_stock_of # Add indices for devices (d), datetimes (j) and commitments (c) model.d = RangeSet(0, len(device_constraints) - 1, doc="Set of devices") @@ -776,12 +478,6 @@ def grouped_commitment_equalities(m, c, j, g): ) model.commitment_sign = Var(model.c, domain=Binary, initialize=0) - def _initial_stock_of(d): - if isinstance(initial_stock, list): - # No initial stock defined for inflexible device - return initial_stock[d] if d < len(initial_stock) else 0 - return initial_stock - def _stock_change_at(m, g, j): """Stock change of stock group g during time step j (before losses).""" return sum( @@ -791,19 +487,6 @@ def _stock_change_at(m, g, j): for dev in group_to_devices[g] ) - def _loss_coefficients(efficiency: float) -> tuple[float, float]: - """Coefficients (a, b) of one step of the stock recursion, for `how="linear"`. - - stock[j] = a * stock[j-1] + b * change[j] - - Mirrors :func:`apply_stock_changes_and_losses`, which we cannot call here - because it expects numbers, while `change[j]` is a Pyomo expression. The - storage efficiency is a Param, so `a` and `b` are plain floats. - """ - if efficiency == 1: - return 1.0, 1.0 - return efficiency, (efficiency - 1) / math.log(efficiency) - def group_stock_balance(m, g, j): """Recursively couple a stock group's stock to the previous step's stock. @@ -815,7 +498,7 @@ def group_stock_balance(m, g, j): (validated above), so the first device can represent the group here. """ d0 = group_to_devices[g][0] - a, b = _loss_coefficients(m.device_efficiency[d0, j]) + 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) @@ -895,7 +578,7 @@ def ems_flow_commitment_equalities(m, c, j): if pd.isna(commodity): devices = m.d else: - devices = commodity_devices.get(commodity, set()) + devices = problem.commodity_devices.get(commodity, set()) if not devices: return Constraint.Skip @@ -1031,27 +714,9 @@ def cost_function(m): if cbc_path is not None: solver.set_executable(cbc_path) - # Set tight tolerance for HiGHS solver - profile = {} - if "highs" in solver_name.lower(): - profile = { - "mip_rel_gap": "0", - "mip_abs_gap": "0", - "primal_feasibility_tolerance": "1e-9", - "dual_feasibility_tolerance": "1e-9", - "mip_feasibility_tolerance": "1e-9", - } - # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO - if current_app.config["LOGGING_LEVEL"] == "INFO": - profile["output_flag"] = "false" - - # Apply operator-configured options last, so they override the defaults above. - configured_options = current_app.config.get("FLEXMEASURES_LP_SOLVER_OPTIONS") or {} - if configured_options and "highs" in solver_name.lower(): - validate_highs_options(configured_options) - profile.update(configured_options) - - for option_name, option_value in profile.items(): + # Tight tolerances for HiGHS, then operator-configured options last + # (shared with the direct HiGHS backend, so both apply the same settings). + for option_name, option_value in solver_options(solver_name).items(): solver.options[option_name] = option_value # load_solutions=False to avoid a RuntimeError exception in appsi solvers when solving an infeasible problem. @@ -1063,44 +728,21 @@ def cost_function(m): planned_costs = value(model.costs) subcommitment_costs = {g: value(cost) for g, cost in model.commitment_costs.items()} - commitment_costs = {} - - # Map subcommitment costs to commitments - for g, v in subcommitment_costs.items(): - c = commitment_mapping[g] - commitment_costs[c] = commitment_costs.get(c, 0) + v - - planned_power_per_device = [] - for d in model.d: - planned_device_power = [model.ems_power[d, j].value for j in model.j] - planned_power_per_device.append( - initialize_series( - data=planned_device_power, - start=start, - end=end, - resolution=to_offset(resolution), - ) - ) - - model.commitment_costs = commitment_costs - commodity_costs = {} - for c in model.c: - commodity = None - if "commodity" in commitments[c].columns: - commodity = commitments[c]["commodity"].iloc[0] - if commodity is None or (isinstance(commodity, float) and np.isnan(commodity)): - continue - cost = value( - model.commitment_downwards_deviation[c] * model.down_price[c] - + model.commitment_upwards_deviation[c] * model.up_price[c] - ) - commodity_costs[commodity] = commodity_costs.get(commodity, 0) + cost + planned_power = planned_power_per_device( + ([model.ems_power[d, j].value for j in model.j] for d in model.d), + start, + end, + resolution, + ) - model.commodity_costs = commodity_costs + model.commitment_costs = aggregate_subcommitment_costs( + subcommitment_costs, commitment_mapping + ) + model.commodity_costs = aggregate_commodity_costs(commitments, subcommitment_costs) # model.pprint() # model.display() # print(results.solver.termination_condition) # print(planned_costs) - return planned_power_per_device, planned_costs, results, model + return planned_power, planned_costs, results, model diff --git a/flexmeasures/data/models/planning/scheduling_problem.py b/flexmeasures/data/models/planning/scheduling_problem.py new file mode 100644 index 0000000000..39ae7d220f --- /dev/null +++ b/flexmeasures/data/models/planning/scheduling_problem.py @@ -0,0 +1,562 @@ +"""Solver-agnostic preparation of the device scheduler's inputs. + +:func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler` +(Pyomo) and +:func:`flexmeasures.data.models.planning.highspy_optimization.device_scheduler_highspy` +(direct HiGHS) build the same mathematical model in two very different +representations, so the model construction itself is necessarily written twice. +Everything *around* it is not: normalising arguments, resolving stock groups, +converting legacy commitments, deriving Big-Ms, and turning solver output back +into schedules and costs is plain pandas/numpy work with no solver in it. + +Keeping that work here means the two backends cannot drift apart on input +handling — only on the model, which is what the equivalence tests in +``tests/test_highspy_equivalence.py`` compare. It also gives both backends a +single place to grow support for a new scheduling feature's *inputs*. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field +from functools import cached_property + +import numpy as np +import pandas as pd +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.utils import initialize_df, initialize_series + +infinity = float("inf") + + +def validate_highs_options(options: dict) -> None: + """Raise if HiGHS would refuse any of these options. + + Pyomo's appsi_highs interface applies solver options without checking HiGHS' + return status, so an unknown name, an invalid value, or a feature missing from + the installed HiGHS build is otherwise ignored without a word. That silently + turns a mis-typed option into a no-op, and a benchmark of it into a false + negative. Probing a throwaway Highs instance surfaces the rejection instead. + """ + try: + import highspy + except ImportError: + # Solver named "*highs*" but highspy absent: let the solver interface complain. + return + + probe = highspy.Highs() + probe.setOptionValue("output_flag", False) + rejected = [ + f"{name}={value!r}" + for name, value in options.items() + if probe.setOptionValue(name, value) != highspy.HighsStatus.kOk + ] + if rejected: + raise ValueError( + f"HiGHS rejected these FLEXMEASURES_LP_SOLVER_OPTIONS: {', '.join(rejected)}." + " The option name may be unknown, the value invalid, or the feature absent" + " from this HiGHS build. For example, the HiPO solver (solver='hipo') needs" + " a HiGHS built against BLAS and METIS, which the pip-installed highspy is not." + ) + + if "threads" in options or "parallel" in options: + current_app.logger.warning( + "FLEXMEASURES_LP_SOLVER_OPTIONS sets 'threads' and/or 'parallel'. HiGHS" + " initializes its thread scheduler once per process, so inside a long-lived" + " worker only the first solve honours these; later solves fail with 'global" + " scheduler has already been initialized' and yield no schedule." + ) + + +def solver_options(solver_name: str) -> dict: + """The solver options to apply, for the given solver. + + HiGHS (whether reached through Pyomo as ``appsi_highs`` or directly as + ``highspy`` -- both match on "highs") gets a tight-tolerance profile, so the + two backends cannot disagree on tolerances and silently produce different + schedules. Operator-configured options are applied last, so they win. + """ + is_highs = "highs" in solver_name.lower() + + profile = {} + if is_highs: + profile = { + "mip_rel_gap": "0", + "mip_abs_gap": "0", + "primal_feasibility_tolerance": "1e-9", + "dual_feasibility_tolerance": "1e-9", + "mip_feasibility_tolerance": "1e-9", + } + # disable logs for the HiGHS solver in case that LOGGING_LEVEL is INFO + if current_app.config["LOGGING_LEVEL"] == "INFO": + profile["output_flag"] = "false" + + configured_options = current_app.config.get("FLEXMEASURES_LP_SOLVER_OPTIONS") or {} + if configured_options and is_highs: + validate_highs_options(configured_options) + profile.update(configured_options) + return profile + + +def convert_commitments_to_subcommitments( + dfs: list[pd.DataFrame], +) -> tuple[list[pd.DataFrame], dict[int, int]]: + """Transform commitments, each specifying a group for each time step, to sub-commitments, one per group. + + 'Groups' are a commitment concept (grouping time slots of a commitment), + making it possible that deviations/breaches can be accounted for properly within this group + (e.g. highest breach per calendar month defines the penalty). + Here, we define sub-commitments, by separating commitments by group and by direction of deviation (up, down). + + We also enumerate the time steps in a new column "j". + + For example, given contracts A and B (represented by 2 DataFrames), each with 3 groups, + we return (sub)commitments A1, A2, A3, B1, B2 and B3, + where A,B,C is the enumerated contract and 1,2,3 is the enumerated group. + """ + commitment_mapping = {} + sub_commitments = [] + for c, df in enumerate(dfs): + # Make sure each commitment has "device" (default NaN) and "class" (default FlowCommitment) columns + if "device" not in df.columns: + df["device"] = np.nan + if "class" not in df.columns: + df["class"] = FlowCommitment + + df["j"] = range(len(df.index)) + + # Group rows by the "group" column in order of first appearance (like + # pd.unique), in a single pass rather than by filtering the DataFrame + # once per group (which would scale quadratically with the number of + # time steps, as each time step often forms its own group). + grouped = df.drop(columns=["group"]).groupby(df["group"], sort=False) + + # Catch non-uniqueness (vectorized across all groups) + if (grouped["upwards deviation price"].nunique(dropna=False) > 1).any(): + raise ValueError( + "Commitment groups cannot have non-unique upwards deviation prices." + ) + if (grouped["downwards deviation price"].nunique(dropna=False) > 1).any(): + raise ValueError( + "Commitment groups cannot have non-unique downwards deviation prices." + ) + + for _, sub_commitment in grouped: + if len(sub_commitment) == 1: + commitment_mapping[len(sub_commitments)] = c + sub_commitments.append(sub_commitment) + else: + down_commitment = sub_commitment.drop(columns="upwards deviation price") + up_commitment = sub_commitment.drop(columns="downwards deviation price") + commitment_mapping[len(sub_commitments)] = c + commitment_mapping[len(sub_commitments) + 1] = c + sub_commitments.extend([down_commitment, up_commitment]) + return sub_commitments, commitment_mapping + + +def loss_coefficients(efficiency: float) -> tuple[float, float]: + """Coefficients (a, b) of one step of the stock recursion, for `how="linear"`. + + stock[j] = a * stock[j-1] + b * change[j] + + Mirrors :func:`apply_stock_changes_and_losses`, which we cannot call here + because it expects numbers, while `change[j]` may be a Pyomo expression. + """ + if efficiency == 1: + return 1.0, 1.0 + return efficiency, (efficiency - 1) / math.log(efficiency) + + +@dataclass +class SchedulingProblem: + """Everything both scheduler backends need before building their model. + + Produced by :func:`prepare_scheduling_problem`; see ``device_scheduler``'s + docstring for what the underlying arguments mean. + """ + + #: Timing, taken from the first device + start: object + end: object + resolution: object + + #: Device constraints, with a "stock delta" column guaranteed to be present + device_constraints: list[pd.DataFrame] + + #: EMS constraints, normalised to a list, plus the device indices each applies to + ems_constraints_list: list[pd.DataFrame] + ems_constraint_device_groups: list[list[int]] + + #: device -> its primary stock group key, and stock group key -> member devices + device_to_group: dict[int, str] + group_to_devices: dict[str, list[int]] + + #: Sub-commitments (one per commitment group and deviation direction), and the + #: mapping from each sub-commitment index back to its original commitment index + commitments: list[pd.DataFrame] + commitment_mapping: dict[int, int] + + #: sub-commitment index -> {device group label -> member device indices} + device_group_lookup: dict[int, dict] + + #: Whether the summed deviation prices describe a convex cost curve + #: (a non-convex curve needs binary commitment-sign variables) + convex_cost_curve: bool + + #: Big-Ms bounding the search space for device power (Md) and commitment + #: deviations (Mc) + Md: float + Mc: float + + #: device index -> its signed power bands (S2 operation modes) + band_lookup: dict[int, list[tuple[float, float]]] + + initial_stock: float | list[float] + + #: The commitments as passed in, before the sub-commitment split. Only kept to + #: derive :attr:`commodity_devices` lazily. + original_commitments: list[pd.DataFrame] = field(default_factory=list, repr=False) + + def initial_stock_of(self, d) -> float: + """The initial stock of device ``d``, defaulting to 0. + + Device indices reaching this from a commitment's "device" column may be + numpy floats, hence the cast. + """ + if isinstance(self.initial_stock, list): + # No initial stock defined for inflexible device + d = int(d) + return self.initial_stock[d] if d < len(self.initial_stock) else 0 + return self.initial_stock + + @cached_property + def commodity_devices(self) -> dict: + """commodity -> set(device indices). + + Computed on demand: only the EMS-level flow commitment constraints need + it, and the per-row scan is not cheap enough to pay for unconditionally. + """ + commodity_devices: dict = {} + for df in self.original_commitments: + if "commodity" not in df.columns or "device" not in df.columns: + continue + + for _, row in df[["commodity", "device"]].dropna().iterrows(): + devices = row["device"] + if not isinstance(devices, (list, tuple, set)): + devices = [devices] + + commodity_devices.setdefault(row["commodity"], set()).update(devices) + return commodity_devices + + +def prepare_scheduling_problem( # noqa C901 + device_constraints: list[pd.DataFrame], + ems_constraints: pd.DataFrame | list[pd.DataFrame], + commitment_quantities: list[pd.Series] | None = None, + commitment_downwards_deviation_price: list[pd.Series] | list[float] | None = None, + commitment_upwards_deviation_price: list[pd.Series] | list[float] | None = None, + commitments: list[pd.DataFrame] | list[Commitment] | None = None, + initial_stock: float | list[float] = 0, + 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, +) -> SchedulingProblem: + """Normalise and validate ``device_scheduler``'s arguments into a SchedulingProblem. + + .. note:: This adds a "stock delta" column to the passed ``device_constraints`` + DataFrames in place, as the schedulers have always done. + """ + # Get timing from first device + start = device_constraints[0].index.to_pydatetime()[0] + # Workaround for https://github.com/pandas-dev/pandas/issues/53643. Was: resolution = pd.to_timedelta(device_constraints[0].index.freq) + resolution = pd.to_timedelta(device_constraints[0].index.freq).to_pytimedelta() + end = device_constraints[0].index.to_pydatetime()[-1] + resolution + + # Normalise EMS constraints to a list of (DataFrame, device-group) pairs. + # A single DataFrame (legacy behaviour) applies to the summed flow of all devices; + # a list of DataFrames applies one EMS-level constraint per device group, as set up + # per commodity by the StorageScheduler. + all_devices = list(range(len(device_constraints))) + if isinstance(ems_constraints, pd.DataFrame): + ems_constraints_list = [ems_constraints] + ems_constraint_device_groups = [all_devices] + else: + ems_constraints_list = ems_constraints + if ems_constraint_groups is None: + if len(ems_constraints_list) > 1: + raise ValueError( + "When passing multiple EMS constraint DataFrames, you must also specify ems_constraint_groups." + ) + ems_constraint_device_groups = [all_devices for _ in ems_constraints_list] + else: + ems_constraint_device_groups = ems_constraint_groups + + # 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 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] = 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[str, 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 = [] + else: + commitments = [ + c.to_frame() if isinstance(c, Commitment) else c for c in commitments + ] + if commitment_quantities is not None: + for quantity, down, up in zip( + commitment_quantities, + commitment_downwards_deviation_price, + commitment_upwards_deviation_price, + ): + + # Turn prices per commitment into prices per commitment flow + if all(isinstance(price, float) for price in down) or isinstance( + down, float + ): + down = initialize_series(down, start, end, resolution) + if all(isinstance(price, float) for price in up) or isinstance(up, float): + up = initialize_series(up, start, end, resolution) + + group = initialize_series(list(range(len(down))), start, end, resolution) + df = initialize_df( + ["quantity", "downwards deviation price", "upwards deviation price"], + start, + end, + resolution, + ) + df["quantity"] = quantity + df["downwards deviation price"] = down + df["upwards deviation price"] = up + df["group"] = group + commitments.append(df) + + # Check if commitments have the same time window and resolution as the constraints + for commitment in commitments: + start_c = commitment.index.to_pydatetime()[0] + resolution_c = pd.to_timedelta(commitment.index.freq) + end_c = commitment.index.to_pydatetime()[-1] + resolution + if not (start_c == start and end_c == end): + raise Exception( + "Not implemented for different time windows.\n(%s,%s)\n(%s,%s)" + % (start, end, start_c, end_c) + ) + if resolution_c != resolution: + raise Exception( + "Not implemented for different resolutions.\n%s\n%s" + % (resolution, resolution_c) + ) + + original_commitments = list(commitments) + commitments, commitment_mapping = convert_commitments_to_subcommitments(commitments) + + device_group_lookup: dict[int, dict] = {} + + 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. + 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 + if commitments: + df = pd.concat(commitments)[ + ["upwards deviation price", "downwards deviation price"] + ] + df = df.groupby(level=0).sum() + convex_cost_curve = ( + len(df[df["upwards deviation price"] < df["downwards deviation price"]]) + == 0 + ) + else: + # No commitments at all: nothing can make the cost curve non-convex. + # (The Pyomo path used to raise on the empty pd.concat here.) + convex_cost_curve = True + + 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]) + + # Both Md and Mc have to be 1 MW, at least + Md = max(Md, 1) + Mc = max(Mc, 1) + + for d in range(len(device_constraints)): + if "stock delta" not in device_constraints[d].columns: + device_constraints[d]["stock delta"] = 0 + else: + device_constraints[d]["stock delta"] = ( + device_constraints[d]["stock delta"].astype(float).fillna(0) + ) + + # Look up power bands (S2 operation modes) per device + if device_power_bands is None: + device_power_bands = [None] * len(device_constraints) + elif len(device_power_bands) != len(device_constraints): + raise ValueError( + f"device_power_bands lists {len(device_power_bands)} devices, " + f"while device_constraints lists {len(device_constraints)} devices." + ) + band_lookup: dict[int, list[tuple[float, float]]] = { + d: list(bands) + for d, bands in enumerate(device_power_bands) + if bands is not None and len(bands) > 0 + } + for d, bands in band_lookup.items(): + for band in bands: + if len(band) != 2 or band[0] > band[1]: + raise ValueError( + f"Invalid power band {band} for device {d}: " + f"expected a (min, max) pair with min <= max." + ) + + return SchedulingProblem( + start=start, + end=end, + resolution=resolution, + device_constraints=device_constraints, + ems_constraints_list=ems_constraints_list, + ems_constraint_device_groups=ems_constraint_device_groups, + device_to_group=device_to_group, + group_to_devices=group_to_devices, + commitments=commitments, + commitment_mapping=commitment_mapping, + device_group_lookup=device_group_lookup, + convex_cost_curve=convex_cost_curve, + Md=Md, + Mc=Mc, + band_lookup=band_lookup, + initial_stock=initial_stock, + original_commitments=original_commitments, + ) + + +def aggregate_subcommitment_costs( + subcommitment_costs: dict, commitment_mapping: dict +) -> dict: + """Sum sub-commitment costs back onto the commitments they were split from.""" + commitment_costs: dict = {} + for g, v in subcommitment_costs.items(): + c = commitment_mapping[g] + commitment_costs[c] = commitment_costs.get(c, 0) + v + return commitment_costs + + +def aggregate_commodity_costs( + commitments: list[pd.DataFrame], subcommitment_costs: dict +) -> dict: + """Sum sub-commitment costs per commodity, skipping commitments without one.""" + commodity_costs: dict = {} + for c in range(len(commitments)): + commodity = None + if "commodity" in commitments[c].columns: + commodity = commitments[c]["commodity"].iloc[0] + if commodity is None or (isinstance(commodity, float) and np.isnan(commodity)): + continue + commodity_costs[commodity] = ( + commodity_costs.get(commodity, 0) + subcommitment_costs[c] + ) + return commodity_costs + + +def planned_power_per_device( + power_per_device, start, end, resolution +) -> list[pd.Series]: + """Turn each device's planned power values into a time series.""" + return [ + initialize_series( + data=list(values), + start=start, + end=end, + resolution=to_offset(resolution), + ) + for values in power_per_device + ] From e9866f7db09d4816f60283d2685320a9fb4fe3a3 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Mon, 3 Aug 2026 15:15:23 +0200 Subject: [PATCH 11/14] perf: build device_group_lookup from column arrays Sub-commitments are usually one row each (every time step tends to form its own commitment group), so this loop slices a fresh two-column DataFrame, runs dropna() on it and calls iterrows() once per time step. The pandas per-call overhead, not the work, dominated: profiling a 4-device x 192-step problem showed 192 dropna() calls accounting for ~50 ms of a ~135 ms prepare_scheduling_problem, against a model build measured in single-digit milliseconds. Read the two columns as arrays once and loop over them instead, replacing dropna() with an explicit missing-value check that also handles the collection-valued "device" entries pd.isna would answer element-wise. The loop itself goes from 49 ms to 1.0 ms (device+device_group) and 65 ms to 0.8 ms (device only) at 192 sub-commitments; prepare_scheduling_problem as a whole drops from 135 ms to 22 ms. Equivalence was checked against the previous implementation over device-only and grouped frames, NaN/None/pd.NA in either column, list/tuple/ndarray device entries, mixed group key types, stock-scoped and empty frames: same groups, same members, same insertion order. Reading the column array yields numpy scalars where iterrows() yielded Python floats on mixed-dtype frames; the two are interchangeable as set members and dict keys (equal hash and equality) and both survive the int() casts applied downstream. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01X2jLjsQzwztDc7Nxz1ozmQ Signed-off-by: F.N. Claessen --- .../models/planning/scheduling_problem.py | 39 ++++++++++++++----- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/flexmeasures/data/models/planning/scheduling_problem.py b/flexmeasures/data/models/planning/scheduling_problem.py index 39ae7d220f..9cfe86b901 100644 --- a/flexmeasures/data/models/planning/scheduling_problem.py +++ b/flexmeasures/data/models/planning/scheduling_problem.py @@ -157,6 +157,17 @@ def convert_commitments_to_subcommitments( return sub_commitments, commitment_mapping +def _is_missing(value) -> bool: + """Whether ``value`` is missing, in the sense ``DataFrame.dropna`` uses. + + A commitment's "device" column may hold a collection of device indices, + which ``pd.isna`` would answer element-wise; such a value is never missing. + """ + if isinstance(value, (list, tuple, set, np.ndarray)): + return False + return bool(pd.isna(value)) + + def loss_coefficients(efficiency: float) -> tuple[float, float]: """Coefficients (a, b) of one step of the stock recursion, for `how="linear"`. @@ -420,28 +431,36 @@ def prepare_scheduling_problem( # noqa C901 continue has_device_group = "device_group" in df.columns + + # Read the columns as arrays rather than slicing + dropna()-ing a fresh + # DataFrame per sub-commitment. Each time step usually forms its own + # group, so this loop runs once per time step, and the per-call pandas + # overhead dominated it (~50 ms of a ~135 ms prepare on 4 devices x 192 + # steps; the arrays bring that under 1 ms). + device_values = df["device"].to_numpy() if has_device_group: - rows = df[["device", "device_group"]].dropna() + group_values = df["device_group"].to_numpy() 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] = {} + group_values = device_values - 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 + groups: dict = {} + for d, g in zip(device_values, group_values): + # Skip what the previous dropna() dropped: a missing device, or a + # missing group label when the commitment declares groups. + if _is_missing(d) or (has_device_group and _is_missing(g)): + continue if isinstance(d, (list, tuple, set, np.ndarray)): devices = set(d) else: devices = {d} - device_group_lookup[c].setdefault(g, set()).update(devices) + groups.setdefault(g, set()).update(devices) + + device_group_lookup[c] = groups # Oversimplified check for a convex cost curve if commitments: From 554af550374d1c0d743003b98e0a4426034e7742 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Mon, 3 Aug 2026 16:22:41 +0200 Subject: [PATCH 12/14] Build the EMS-level flow commitment rows in the direct HiGHS backend This backend skipped ems_flow_commitment_equalities because the Pyomo path built those rows without bounds, so they could never bind and HiGHS dropped them anyway. PR #2355 gives them the same one-sided bounds that grouped_commitment_equalities uses, so they bind now, and skipping them here would silently ignore an EMS-level commitment under this backend -- which this PR makes the default. The row is the grouped one with a different summation set, so rather than write it twice, the existing loop's body is extracted into _active_rows (the commitment's active time steps and bounds) and _add_commitment_rows (bind a commitment to the summed flow or stock of a set of devices). The EMS-level loop then sums over every device, or over the commitment's commodity's devices, and mirrors the Pyomo path in skipping commitments that name a device group -- those are already bound per group, and binding them twice over-constrains the problem. Two equivalence scenarios cover this, so it runs under both backends: ems_level_flow_commitment (names no device, binds all devices) and ems_level_commodity_commitment (binds only its commodity's devices, which also exercises the commodity_devices lookup). Both were checked to fail with the new rows disabled, so they cannot pass vacuously. Also drops the now-inaccurate deviation note from the module docstring. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01X2jLjsQzwztDc7Nxz1ozmQ Signed-off-by: F.N. Claessen --- .../models/planning/highspy_optimization.py | 144 +++++++++--------- .../tests/test_highspy_equivalence.py | 64 ++++++++ 2 files changed, 140 insertions(+), 68 deletions(-) diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 50b258000c..899e903907 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -18,33 +18,6 @@ Deviations from the Pyomo implementation (all verified against the behavior of the ``appsi_highs`` path): -- ``ems_flow_commitment_equalities`` is not built. On the Pyomo path this - constraint family returns ``(None, expr, None)``, i.e. a constraint without - bounds, which ends up as a free (vacuous) row in HiGHS. We skip building the - free rows altogether. - - .. note:: This deviation is only harmless for as long as those rows stay free, - and two branches in flight change that: - - - #2355 gives the constraint the same one-sided bounds that - ``grouped_commitment_equalities`` already uses, making it bind. Once that - lands, this module must build the row too — otherwise an EMS-level - commitment is silently ignored under this backend. The row is then - structurally identical to the grouped one, differing only in the set it - sums over (all devices, or the commodity's devices). - - #2380 makes an unscoped flex-context commitment device-grouped - (``device=``, ``device_group=``) - rather than per-device, so commitments coming from ``convert_to_commitments`` - route through ``grouped_commitment_equalities``, which this module does - build. That narrows the exposure above to callers constructing an - EMS-level ``FlowCommitment`` (``device=None``) directly, but does not - remove it. - - Note that ``tests/test_commitments.py`` does not use the - ``app_with_each_solver`` fixture, so its cases only ever run under the - configured default solver -- which ``config_defaults`` now sets to - ``"highspy"``. Coverage of this constraint family belongs in - ``tests/test_highspy_equivalence.py`` to run under both backends. - Rows whose computed bounds are impossible to satisfy for any finite value (upper bound of -inf, or lower bound of +inf, as happens when a commitment quantity is +/-inf) are skipped. On the Pyomo path such rows are rejected by @@ -62,6 +35,7 @@ from flexmeasures.data.models.planning import ( Commitment, + FlowCommitment, StockCommitment, ) from flexmeasures.data.models.planning.scheduling_problem import ( @@ -538,57 +512,91 @@ def _price_of(df: pd.DataFrame, column: str) -> float: # lb <= quantity + down_dev + up_dev - sum_over_group <= ub # where lb is 0 iff the commitment prices upwards deviations and ub is 0 # iff it prices downwards deviations (one-sided otherwise). - # NB the ems_flow_commitment_equalities of the Pyomo implementation are - # deliberately not built here; they are free rows (see module docstring). - for c, df in enumerate(commitments): - groups = device_group_lookup.get(c, {}) - if not groups: - continue + def _active_rows(df: pd.DataFrame): + """The commitment's active time steps and its row bounds. + + A NaN quantity deactivates the commitment at that time step. The Pyomo + implementation maps such a quantity to -inf in its Param and lets the + resulting row (whose lower bound works out to +inf) be rejected by HiGHS; + dropping it here has the same effect. + """ quantity = _column(df, "quantity") jj = df["j"].to_numpy(dtype=np.int64) - # A NaN quantity deactivates the commitment at that time step - # (NaN was mapped to -inf in the Pyomo implementation's Param). active = ~(np.isnan(quantity) | (quantity == -infinity)) - if not np.any(active): - continue - quantity = quantity[active] - jj = jj[active] lb = 0.0 if "upwards deviation price" in df.columns else -infinity ub = 0.0 if "downwards deviation price" in df.columns else infinity - is_stock = df["class"].apply(lambda cl: cl == StockCommitment).all() + return quantity[active], jj[active], lb, ub + + def _add_commitment_rows(c, quantity, jj, lb, ub, devices, is_stock) -> None: + """Bind commitment ``c`` to the summed flow or stock of ``devices``.""" n_rows = len(jj) + idx_cols = [ + np.full(n_rows, col_cdown + c), + np.full(n_rows, col_cup + c), + ] + val_cols = [np.ones(n_rows), np.ones(n_rows)] + if is_stock: + # Aggregate coefficients per stock group column and move the + # initial stocks into the row bounds. + stock_coefficients: dict[int, float] = {} + initial_stock_sum = 0.0 + for dev in devices: + base = col_stock + group_index[device_to_group[int(dev)]] * T + stock_coefficients[base] = stock_coefficients.get(base, 0.0) - 1.0 + initial_stock_sum += _initial_stock_of(dev) + for base, coefficient in stock_coefficients.items(): + idx_cols.append(base + jj) + val_cols.append(np.full(n_rows, coefficient)) + offset = initial_stock_sum + else: + for dev in devices: + idx_cols.append(col_ems + int(dev) * T + jj) + val_cols.append(np.full(n_rows, -1.0)) + offset = 0.0 + rows.add_uniform_rows( + lb - quantity - offset, + ub - quantity - offset, + np.column_stack(idx_cols), + np.column_stack(val_cols), + ) + + for c, df in enumerate(commitments): + groups = device_group_lookup.get(c, {}) + if not groups: + continue + quantity, jj, lb, ub = _active_rows(df) + if len(jj) == 0: + continue + is_stock = df["class"].apply(lambda cl: cl == StockCommitment).all() for g, devices_in_group in groups.items(): if not devices_in_group: continue - idx_cols = [ - np.full(n_rows, col_cdown + c), - np.full(n_rows, col_cup + c), - ] - val_cols = [np.ones(n_rows), np.ones(n_rows)] - if is_stock: - # Aggregate coefficients per stock group column and move the - # initial stocks into the row bounds. - stock_coefficients: dict[int, float] = {} - initial_stock_sum = 0.0 - for dev in devices_in_group: - base = col_stock + group_index[device_to_group[int(dev)]] * T - stock_coefficients[base] = stock_coefficients.get(base, 0.0) - 1.0 - initial_stock_sum += _initial_stock_of(dev) - for base, coefficient in stock_coefficients.items(): - idx_cols.append(base + jj) - val_cols.append(np.full(n_rows, coefficient)) - offset = initial_stock_sum + _add_commitment_rows(c, quantity, jj, lb, ub, devices_in_group, is_stock) + + # ems_flow_commitment_equalities: an EMS-level flow commitment binds the + # summed flow of every device, or of its commodity's devices when it names a + # commodity. A commitment that names devices is skipped here, being already + # bound per device group above; binding it twice would over-constrain it. + for c, df in enumerate(commitments): + if device_group_lookup.get(c): + continue + if df["class"].iloc[0] != FlowCommitment: + continue + if "commodity" not in df.columns: + # Legacy behavior: no commodity, so sum over all devices. + devices: object = range(D) + else: + commodity = df["commodity"].iloc[0] + if pd.isna(commodity): + devices = range(D) else: - for dev in devices_in_group: - idx_cols.append(col_ems + int(dev) * T + jj) - val_cols.append(np.full(n_rows, -1.0)) - offset = 0.0 - rows.add_uniform_rows( - lb - quantity - offset, - ub - quantity - offset, - np.column_stack(idx_cols), - np.column_stack(val_cols), - ) + devices = problem.commodity_devices.get(commodity, set()) + if not devices: + continue + quantity, jj, lb, ub = _active_rows(df) + if len(jj) == 0: + continue + _add_commitment_rows(c, quantity, jj, lb, ub, devices, is_stock=False) # Power bands (S2 operation modes): each banded device runs in exactly one # band per time step, and its power must lie within the chosen band. diff --git a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py index 48d4307823..a6f88c2723 100644 --- a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py +++ b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py @@ -256,6 +256,68 @@ def run_with_solver(app, solver: str, make_scenario): app.config["FLEXMEASURES_LP_SOLVER"] = original_solver +def scenario_ems_level_flow_commitment(): + """Two devices under an EMS-level flow commitment, which names no device. + + Such a commitment binds the summed flow of all devices, via + ``ems_flow_commitment_equalities`` rather than the grouped constraints, so it + is the case that distinguishes the two constraint families. + """ + index = make_index() + prices = make_prices(index) + return dict( + device_constraints=[make_battery_constraints(), make_battery_constraints()], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=[ + FlowCommitment( + name="EMS target", + index=index, + quantity=0.1, + upwards_deviation_price=prices, + downwards_deviation_price=prices, + ) + ], + initial_stock=[0.5, 0.5], + ) + + +def scenario_ems_level_commodity_commitment(): + """An EMS-level flow commitment scoped to one commodity's devices. + + Device 0 carries the commodity, device 1 does not, so the commitment must + bind device 0's flow only -- exercising the commodity_devices lookup rather + than the sum-over-all-devices fallback. + """ + index = make_index() + prices = make_prices(index) + commodity_commitment = FlowCommitment( + name="gas target", + index=index, + quantity=0.1, + upwards_deviation_price=prices, + downwards_deviation_price=prices, + ) + frame = commodity_commitment.to_frame() + frame["commodity"] = "gas" + # Name the commodity's device without scoping the commitment to a device + # group, so it still routes through ems_flow_commitment_equalities. + scoping = FlowCommitment( + name="gas scope", + index=index, + quantity=0, + upwards_deviation_price=0, + downwards_deviation_price=0, + device=pd.Series(0, index=index), + ).to_frame() + scoping["commodity"] = "gas" + return dict( + device_constraints=[make_battery_constraints(), make_battery_constraints()], + ems_constraints=initialize_df(COLUMNS, START, END, RESOLUTION), + commitments=[frame, scoping], + initial_stock=[0.5, 0.5], + ) + + @pytest.mark.parametrize( "make_scenario", [ @@ -263,6 +325,8 @@ def run_with_solver(app, solver: str, make_scenario): scenario_battery_with_soc_targets, scenario_battery_with_site_capacity_and_breach_prices, scenario_two_devices_with_stock_commitment, + scenario_ems_level_flow_commitment, + scenario_ems_level_commodity_commitment, ], ids=lambda f: f.__name__.replace("scenario_", ""), ) From 64e2422532d838a378bafaf9efaa8261520b472c Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Mon, 3 Aug 2026 16:49:49 +0200 Subject: [PATCH 13/14] review: name the solver-results shim after Pyomo's own class "Stanza" was an unhelpful coinage. The attribute this shim stands in for holds a pyomo.opt.results.solver.SolverInformation, so name it _SolverInformation and say where the name comes from. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01X2jLjsQzwztDc7Nxz1ozmQ Signed-off-by: F.N. Claessen --- .../data/models/planning/highspy_optimization.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 899e903907..1e53eef154 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -49,8 +49,12 @@ infinity = float("inf") -class _SolverStanza: - """Mimics the ``solver`` entry of a Pyomo ``SolverResults`` object.""" +class _SolverInformation: + """Mimics the ``solver`` entry of a Pyomo ``SolverResults`` object. + + Named after Pyomo's own ``pyomo.opt.results.solver.SolverInformation``, which + is what that entry holds. + """ def __init__(self, termination_condition: str, status: str): #: str containing "optimal", "infeasible", etc. (mirrors Pyomo's @@ -68,7 +72,7 @@ class HighspySolverResults: """ def __init__(self, termination_condition: str, status: str): - self.solver = _SolverStanza(termination_condition, status) + self.solver = _SolverInformation(termination_condition, status) class _IndexedVarView: From bc0f13d5a5edf9825a38fb9f491e7a9737d5b7c1 Mon Sep 17 00:00:00 2001 From: "F.N. Claessen" Date: Mon, 3 Aug 2026 16:58:20 +0200 Subject: [PATCH 14/14] style: reflow docstrings and comments to break only after punctuation CLAUDE.md asks that docstrings and comments break lines only after punctuation, never mid-phrase, with max-line-length 160 and E501 ignored, so that review comments and text search stay stable. Copilot flagged nine places in this PR where I had not followed it, and it was right. Reflowed the docstrings and comments this PR adds, plus the ones it moved into the new scheduling_problem module. Doing that here rather than in the refactor commit keeps that commit verifiable as a verbatim move. Text only. No behaviour, names or logic changed; 290 passed, 3 xfailed across the planning suite under all three solvers. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01X2jLjsQzwztDc7Nxz1ozmQ Signed-off-by: F.N. Claessen --- .../models/planning/highspy_optimization.py | 105 ++++++++---------- .../models/planning/linear_optimization.py | 73 ++++++------ .../models/planning/scheduling_problem.py | 98 ++++++++-------- .../tests/test_highspy_equivalence.py | 39 ++++--- .../data/models/planning/tests/test_solver.py | 8 +- 5 files changed, 151 insertions(+), 172 deletions(-) diff --git a/flexmeasures/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py index 1e53eef154..0dbae6eb3e 100644 --- a/flexmeasures/data/models/planning/highspy_optimization.py +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -2,30 +2,26 @@ .. warning:: TWO MODELS TO KEEP IN SYNC - This module deliberately duplicates the mathematical model of - :func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler` - (the Pyomo implementation), building the LP/MILP directly with the HiGHS - Python API (``highspy``) instead. The Pyomo implementation is the semantic - reference: any change to the variables, constraints or objective in - ``device_scheduler`` MUST be mirrored here (and vice versa), and the - equivalence tests in ``tests/test_highspy_equivalence.py`` should be - extended accordingly. This trade-off (a second model to maintain) was - accepted because bypassing the Pyomo layer cuts roughly a second (single - device) to several seconds (multiple devices) of model construction and - solution-ingestion overhead per scheduling job, while the direct build - takes milliseconds. - -Deviations from the Pyomo implementation (all verified against the behavior of -the ``appsi_highs`` path): + This module deliberately duplicates the mathematical model of the Pyomo implementation, + :func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler`, + building the LP/MILP directly with the HiGHS Python API (``highspy``) instead. + The Pyomo implementation is the semantic reference: + any change to the variables, constraints or objective in ``device_scheduler`` MUST be mirrored here (and vice versa), + and the equivalence tests in ``tests/test_highspy_equivalence.py`` should be extended accordingly. + This trade-off (a second model to maintain) was accepted + because bypassing the Pyomo layer cuts roughly a second (single device) to several seconds (multiple devices) + of model construction and solution-ingestion overhead per scheduling job, + while the direct build takes milliseconds. + +Deviations from the Pyomo implementation (all verified against the behavior of the ``appsi_highs`` path): - Rows whose computed bounds are impossible to satisfy for any finite value - (upper bound of -inf, or lower bound of +inf, as happens when a commitment - quantity is +/-inf) are skipped. On the Pyomo path such rows are rejected by - HiGHS' ``addRow`` (called by appsi) and thereby silently dropped, with the - same net effect. -- Solver results and model objects are lightweight shims (see - :class:`HighspySolverResults` and :class:`HighspyModel`) that expose the - attributes callers actually consume, rather than Pyomo objects. + (upper bound of -inf, or lower bound of +inf, as happens when a commitment quantity is +/-inf) are skipped. + On the Pyomo path such rows are rejected by HiGHS' ``addRow`` (called by appsi) and thereby silently dropped, + with the same net effect. +- Solver results and model objects are lightweight shims + (see :class:`HighspySolverResults` and :class:`HighspyModel`) + that expose the attributes callers actually consume, rather than Pyomo objects. """ from __future__ import annotations @@ -52,14 +48,13 @@ class _SolverInformation: """Mimics the ``solver`` entry of a Pyomo ``SolverResults`` object. - Named after Pyomo's own ``pyomo.opt.results.solver.SolverInformation``, which - is what that entry holds. + Named after Pyomo's own ``pyomo.opt.results.solver.SolverInformation``, which is what that entry holds. """ def __init__(self, termination_condition: str, status: str): - #: str containing "optimal", "infeasible", etc. (mirrors Pyomo's - #: str-valued TerminationCondition enum, which supports both - #: ``== "optimal"`` and ``"infeasible" in ...`` checks) + #: str containing "optimal", "infeasible", etc. + #: Mirrors Pyomo's str-valued TerminationCondition enum, + #: which supports both ``== "optimal"`` and ``"infeasible" in ...`` checks. self.termination_condition = termination_condition self.status = status @@ -67,8 +62,8 @@ def __init__(self, termination_condition: str, status: str): class HighspySolverResults: """Small stand-in for Pyomo's ``SolverResults``. - Callers only consume ``results.solver.termination_condition`` (a string - containing "optimal"/"infeasible") and ``results.solver.status``. + Callers only consume ``results.solver.termination_condition`` (a string containing "optimal"/"infeasible") + and ``results.solver.status``. """ def __init__(self, termination_condition: str, status: str): @@ -105,12 +100,11 @@ class HighspyModel: - ``commitment_costs``: dict of realized costs per (original) commitment - ``commodity_costs``: dict of realized costs per commodity - - ``costs``: the objective value (a float; ``pyomo.environ.value()`` passes - floats through unchanged, so ``value(model.costs)`` keeps working) + - ``costs``: the objective value + (a float; ``pyomo.environ.value()`` passes floats through unchanged, so ``value(model.costs)`` keeps working) - ``d`` and ``j``: the device and datetime index ranges - - ``ems_power``, ``device_power_up``, ``device_power_down``, - ``device_power_sign``: indexed variable views supporting - ``var[d, j].value`` and ``var.extract_values()`` + - ``ems_power``, ``device_power_up``, ``device_power_down``, ``device_power_sign``: + indexed variable views supporting ``var[d, j].value`` and ``var.extract_values()`` """ def __init__(self): @@ -133,8 +127,8 @@ def _column(df: pd.DataFrame, name: str) -> np.ndarray: def _column_or_default(df: pd.DataFrame, name: str, default: float) -> np.ndarray: """Return a DataFrame column as a float array, with missing column or NaN values replaced by a default. - Mirrors e.g. ``device_efficiency`` in the Pyomo implementation ("assume - perfect efficiency if no efficiency information is available"). + Mirrors e.g. ``device_efficiency`` in the Pyomo implementation + ("assume perfect efficiency if no efficiency information is available"). """ if name not in df.columns: return np.full(len(df), float(default)) @@ -160,10 +154,9 @@ def _loss_coefficient_arrays(efficiency: np.ndarray) -> tuple[np.ndarray, np.nda class _RowBuilder: """Accumulates constraint rows (in CSR form) for a single HiGHS addRows call. - Rows that no finite assignment can satisfy (upper bound -inf or lower bound - +inf) are skipped, mirroring how HiGHS rejects such rows when the appsi - interface adds them one by one (see module docstring). Free rows - (-inf, +inf) are also skipped, as they cannot bind. + Rows that no finite assignment can satisfy (upper bound -inf or lower bound +inf) are skipped, + mirroring how HiGHS rejects such rows when the appsi interface adds them one by one (see module docstring). + Free rows (-inf, +inf) are also skipped, as they cannot bind. """ def __init__(self): @@ -510,18 +503,18 @@ def _price_of(df: pd.DataFrame, column: str) -> float: np.tile([-1.0, Mc], (C, 1)), ) - # grouped_commitment_equalities: couple each commitment's baseline (plus - # deviation variables) to the summed flow (FlowCommitment) or stock - # (StockCommitment) of each of its device groups: + # grouped_commitment_equalities: + # couple each commitment's baseline (plus deviation variables) + # to the summed flow (FlowCommitment) or stock (StockCommitment) of each of its device groups: # lb <= quantity + down_dev + up_dev - sum_over_group <= ub - # where lb is 0 iff the commitment prices upwards deviations and ub is 0 - # iff it prices downwards deviations (one-sided otherwise). + # where lb is 0 iff the commitment prices upwards deviations, + # and ub is 0 iff it prices downwards deviations (one-sided otherwise). def _active_rows(df: pd.DataFrame): """The commitment's active time steps and its row bounds. - A NaN quantity deactivates the commitment at that time step. The Pyomo - implementation maps such a quantity to -inf in its Param and lets the - resulting row (whose lower bound works out to +inf) be rejected by HiGHS; + A NaN quantity deactivates the commitment at that time step. + The Pyomo implementation maps such a quantity to -inf in its Param + and lets the resulting row (whose lower bound works out to +inf) be rejected by HiGHS; dropping it here has the same effect. """ quantity = _column(df, "quantity") @@ -540,8 +533,8 @@ def _add_commitment_rows(c, quantity, jj, lb, ub, devices, is_stock) -> None: ] val_cols = [np.ones(n_rows), np.ones(n_rows)] if is_stock: - # Aggregate coefficients per stock group column and move the - # initial stocks into the row bounds. + # Aggregate coefficients per stock group column, + # and move the initial stocks into the row bounds. stock_coefficients: dict[int, float] = {} initial_stock_sum = 0.0 for dev in devices: @@ -577,10 +570,10 @@ def _add_commitment_rows(c, quantity, jj, lb, ub, devices, is_stock) -> None: continue _add_commitment_rows(c, quantity, jj, lb, ub, devices_in_group, is_stock) - # ems_flow_commitment_equalities: an EMS-level flow commitment binds the - # summed flow of every device, or of its commodity's devices when it names a - # commodity. A commitment that names devices is skipped here, being already - # bound per device group above; binding it twice would over-constrain it. + # ems_flow_commitment_equalities: an EMS-level flow commitment binds the summed flow of every device, + # or of its commodity's devices when it names a commodity. + # A commitment that names devices is skipped here, being already bound per device group above; + # binding it twice would over-constrain it. for c, df in enumerate(commitments): if device_group_lookup.get(c): continue @@ -662,8 +655,8 @@ def _add_commitment_rows(c, quantity, jj, lb, ub, devices, is_stock) -> None: if nrow > 0: h.addRows(nrow, row_lower, row_upper, nnz, row_starts, a_index, a_value) - # The same options the Pyomo path applies for HiGHS solvers ("highspy" matches - # on "highs"), so the two backends cannot disagree on tolerances. + # The same options the Pyomo path applies for HiGHS solvers ("highspy" matches on "highs"), + # so the two backends cannot disagree on tolerances. for option_name, option_value in solver_options("highspy").items(): h.setOptionValue(option_name, option_value) diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index 14bfdf7304..51b94a0957 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -46,10 +46,9 @@ def _left_at_default(value, default) -> bool: """Whether an argument was left at its default. - Best-effort: pandas values compare element-wise, so ``value == default`` may - return an array (or raise) rather than a bool. Anything we cannot decide is - reported as "not the default", which errs towards raising in - :func:`_arguments_for_highspy_backend` rather than silently dropping a value. + Best-effort: pandas values compare element-wise, so ``value == default`` may return an array (or raise) rather than a bool. + Anything we cannot decide is reported as "not the default", + which errs towards raising in :func:`_arguments_for_highspy_backend` rather than silently dropping a value. """ if value is default: return True @@ -65,10 +64,8 @@ def _left_at_default(value, default) -> bool: def _backend_argument_map(declared_by, supported_by) -> tuple[frozenset, tuple]: """Which of ``declared_by``'s arguments ``supported_by`` accepts, and which it lacks. - Signatures are static, so this is computed once per process (roughly 70 us, - which is not worth paying on every schedule). Caching on the two function - objects rather than on nothing keeps it correct when a test substitutes one - of them. + Signatures are static, so this is computed once per process (roughly 70 us, which is not worth paying on every schedule). + Caching on the two function objects rather than on nothing keeps it correct when a test substitutes one of them. """ declared = inspect.signature(declared_by).parameters supported = frozenset(inspect.signature(supported_by).parameters) @@ -83,23 +80,22 @@ def _backend_argument_map(declared_by, supported_by) -> tuple[frozenset, tuple]: def _arguments_for_highspy_backend(passed_arguments: dict) -> dict: """Map ``device_scheduler``'s arguments onto the direct HiGHS backend's signature. - A hand-written keyword list here would be a trap: whoever adds the next - ``device_scheduler`` parameter naturally works on the Pyomo model further - down this file, and a parameter missing from that list would not fail — it - would simply never reach the backend. Under ``FLEXMEASURES_LP_SOLVER="highspy"`` - (the default) that yields a schedule computed as if the constraint had never - been requested: plausible-looking, silently wrong, and not caught by tests - written before the default was flipped. - - Forwarding by name removes that failure mode entirely. The remaining case — - an argument the direct backend does not model at all — is caught statically - by ``test_every_device_scheduler_argument_currently_reaches_the_backend``, so - it cannot reach a release. The raise below is only a backstop for a build - where that test did not run; it costs nothing while the signatures agree. - - This is a live concern rather than a hypothetical one: ``device_scheduler`` - is gaining ``coupling_groups`` (#2218) and ``balance_groups`` (#2289) on - branches in flight, and each needs explicit support here. + A hand-written keyword list here would be a trap: + whoever adds the next ``device_scheduler`` parameter naturally works on the Pyomo model further down this file, + and a parameter missing from that list would not fail — it would simply never reach the backend. + Under ``FLEXMEASURES_LP_SOLVER="highspy"`` (the default), + that yields a schedule computed as if the constraint had never been requested: + plausible-looking, silently wrong, and not caught by tests written before the default was flipped. + + Forwarding by name removes that failure mode entirely. + The remaining case — an argument the direct backend does not model at all — + is caught statically by ``test_every_device_scheduler_argument_currently_reaches_the_backend``, so it cannot reach a release. + The raise below is only a backstop for a build where that test did not run; + it costs nothing while the signatures agree. + + This is a live concern rather than a hypothetical one: + ``device_scheduler`` is gaining ``coupling_groups`` (#2218) and ``balance_groups`` (#2289) on branches in flight, + and each needs explicit support here. """ from flexmeasures.data.models.planning.highspy_optimization import ( device_scheduler_highspy, @@ -198,13 +194,13 @@ def device_scheduler( # noqa C901 DataFrame. Later we could pass in a MultiIndex DataFrame directly. """ - # The "highspy" solver choice bypasses Pyomo altogether: the same model is - # built directly with the HiGHS Python API (much faster to construct). - # See the highspy_optimization module, which mirrors the model built below - # and returns compatible result objects. + # The "highspy" solver choice bypasses Pyomo altogether: + # the same model is built directly with the HiGHS Python API (much faster to construct). + # See the highspy_optimization module, which mirrors the model built below and returns compatible result objects. if current_app.config.get("FLEXMEASURES_LP_SOLVER") == "highspy": - # Arguments are forwarded by name, and an argument the direct backend cannot - # model raises rather than being dropped. See _arguments_for_highspy_backend. + # Arguments are forwarded by name, + # and an argument the direct backend cannot model raises rather than being dropped. + # See _arguments_for_highspy_backend. highspy_arguments = _arguments_for_highspy_backend(locals()) from flexmeasures.data.models.planning.highspy_optimization import ( @@ -232,8 +228,8 @@ def device_scheduler( # noqa C901 device_power_bands=device_power_bands, ) - # Local aliases, so that the model below reads as it did before the (solver-agnostic) - # input handling moved to the scheduling_problem module. + # Local aliases, + # so that the model below reads as it did before the (solver-agnostic) input handling moved to the scheduling_problem module. start, end, resolution = problem.start, problem.end, problem.resolution device_constraints = problem.device_constraints ems_constraints_list = problem.ems_constraints_list @@ -570,12 +566,11 @@ def ems_flow_commitment_equalities(m, c, j): if commitments[c]["class"].iloc[0] != FlowCommitment: return Constraint.Skip - # A device-scoped commitment is already bound, once per device group, by - # grouped_commitment_equalities. Now that this constraint family actually - # has bounds, binding such a commitment here as well would constrain the - # same deviation variables a second time, against a different device set - # (the whole EMS, or the whole commodity). Only genuinely EMS-level - # commitments -- those naming no device -- belong here. + # A device-scoped commitment is already bound, once per device group, by grouped_commitment_equalities. + # Now that this constraint family actually has bounds, + # binding such a commitment here as well would constrain the same deviation variables a second time, + # against a different device set (the whole EMS, or the whole commodity). + # Only genuinely EMS-level commitments -- those naming no device -- belong here. if device_group_lookup.get(c): return Constraint.Skip diff --git a/flexmeasures/data/models/planning/scheduling_problem.py b/flexmeasures/data/models/planning/scheduling_problem.py index 9cfe86b901..92b48d349d 100644 --- a/flexmeasures/data/models/planning/scheduling_problem.py +++ b/flexmeasures/data/models/planning/scheduling_problem.py @@ -1,18 +1,16 @@ """Solver-agnostic preparation of the device scheduler's inputs. -:func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler` -(Pyomo) and -:func:`flexmeasures.data.models.planning.highspy_optimization.device_scheduler_highspy` -(direct HiGHS) build the same mathematical model in two very different -representations, so the model construction itself is necessarily written twice. -Everything *around* it is not: normalising arguments, resolving stock groups, -converting legacy commitments, deriving Big-Ms, and turning solver output back -into schedules and costs is plain pandas/numpy work with no solver in it. - -Keeping that work here means the two backends cannot drift apart on input -handling — only on the model, which is what the equivalence tests in -``tests/test_highspy_equivalence.py`` compare. It also gives both backends a -single place to grow support for a new scheduling feature's *inputs*. +:func:`flexmeasures.data.models.planning.linear_optimization.device_scheduler` (Pyomo) and +:func:`flexmeasures.data.models.planning.highspy_optimization.device_scheduler_highspy` (direct HiGHS) +build the same mathematical model in two very different representations, +so the model construction itself is necessarily written twice. +Everything *around* it is not: +normalising arguments, resolving stock groups, converting legacy commitments, deriving Big-Ms, +and turning solver output back into schedules and costs is plain pandas/numpy work with no solver in it. + +Keeping that work here means the two backends cannot drift apart on input handling — +only on the model, which is what the equivalence tests in ``tests/test_highspy_equivalence.py`` compare. +It also gives both backends a single place to grow support for a new scheduling feature's *inputs*. """ from __future__ import annotations @@ -35,11 +33,10 @@ def validate_highs_options(options: dict) -> None: """Raise if HiGHS would refuse any of these options. - Pyomo's appsi_highs interface applies solver options without checking HiGHS' - return status, so an unknown name, an invalid value, or a feature missing from - the installed HiGHS build is otherwise ignored without a word. That silently - turns a mis-typed option into a no-op, and a benchmark of it into a false - negative. Probing a throwaway Highs instance surfaces the rejection instead. + Pyomo's appsi_highs interface applies solver options without checking HiGHS' return status, + so an unknown name, an invalid value, or a feature missing from the installed HiGHS build is otherwise ignored without a word. + That silently turns a mis-typed option into a no-op, and a benchmark of it into a false negative. + Probing a throwaway Highs instance surfaces the rejection instead. """ try: import highspy @@ -74,10 +71,10 @@ def validate_highs_options(options: dict) -> None: def solver_options(solver_name: str) -> dict: """The solver options to apply, for the given solver. - HiGHS (whether reached through Pyomo as ``appsi_highs`` or directly as - ``highspy`` -- both match on "highs") gets a tight-tolerance profile, so the - two backends cannot disagree on tolerances and silently produce different - schedules. Operator-configured options are applied last, so they win. + HiGHS (whether reached through Pyomo as ``appsi_highs`` or directly as ``highspy`` -- both match on "highs") + gets a tight-tolerance profile, + so the two backends cannot disagree on tolerances and silently produce different schedules. + Operator-configured options are applied last, so they win. """ is_highs = "highs" in solver_name.lower() @@ -107,8 +104,7 @@ def convert_commitments_to_subcommitments( """Transform commitments, each specifying a group for each time step, to sub-commitments, one per group. 'Groups' are a commitment concept (grouping time slots of a commitment), - making it possible that deviations/breaches can be accounted for properly within this group - (e.g. highest breach per calendar month defines the penalty). + making it possible that deviations/breaches can be accounted for properly within this group (e.g. highest breach per calendar month defines the penalty). Here, we define sub-commitments, by separating commitments by group and by direction of deviation (up, down). We also enumerate the time steps in a new column "j". @@ -128,10 +124,9 @@ def convert_commitments_to_subcommitments( df["j"] = range(len(df.index)) - # Group rows by the "group" column in order of first appearance (like - # pd.unique), in a single pass rather than by filtering the DataFrame - # once per group (which would scale quadratically with the number of - # time steps, as each time step often forms its own group). + # Group rows by the "group" column in order of first appearance (like pd.unique), + # in a single pass rather than by filtering the DataFrame once per group + # (which would scale quadratically with the number of time steps, as each time step often forms its own group). grouped = df.drop(columns=["group"]).groupby(df["group"], sort=False) # Catch non-uniqueness (vectorized across all groups) @@ -160,8 +155,8 @@ def convert_commitments_to_subcommitments( def _is_missing(value) -> bool: """Whether ``value`` is missing, in the sense ``DataFrame.dropna`` uses. - A commitment's "device" column may hold a collection of device indices, - which ``pd.isna`` would answer element-wise; such a value is never missing. + A commitment's "device" column may hold a collection of device indices, which ``pd.isna`` would answer element-wise; + such a value is never missing. """ if isinstance(value, (list, tuple, set, np.ndarray)): return False @@ -173,8 +168,8 @@ def loss_coefficients(efficiency: float) -> tuple[float, float]: stock[j] = a * stock[j-1] + b * change[j] - Mirrors :func:`apply_stock_changes_and_losses`, which we cannot call here - because it expects numbers, while `change[j]` may be a Pyomo expression. + Mirrors :func:`apply_stock_changes_and_losses`, which we cannot call here because it expects numbers, + while `change[j]` may be a Pyomo expression. """ if efficiency == 1: return 1.0, 1.0 @@ -185,8 +180,8 @@ def loss_coefficients(efficiency: float) -> tuple[float, float]: class SchedulingProblem: """Everything both scheduler backends need before building their model. - Produced by :func:`prepare_scheduling_problem`; see ``device_scheduler``'s - docstring for what the underlying arguments mean. + Produced by :func:`prepare_scheduling_problem`; + see ``device_scheduler``'s docstring for what the underlying arguments mean. """ #: Timing, taken from the first device @@ -214,11 +209,10 @@ class SchedulingProblem: device_group_lookup: dict[int, dict] #: Whether the summed deviation prices describe a convex cost curve - #: (a non-convex curve needs binary commitment-sign variables) + #: (a non-convex curve needs binary commitment-sign variables). convex_cost_curve: bool - #: Big-Ms bounding the search space for device power (Md) and commitment - #: deviations (Mc) + #: Big-Ms bounding the search space for device power (Md) and commitment deviations (Mc) Md: float Mc: float @@ -227,15 +221,14 @@ class SchedulingProblem: initial_stock: float | list[float] - #: The commitments as passed in, before the sub-commitment split. Only kept to - #: derive :attr:`commodity_devices` lazily. + #: The commitments as passed in, before the sub-commitment split. + #: Only kept to derive :attr:`commodity_devices` lazily. original_commitments: list[pd.DataFrame] = field(default_factory=list, repr=False) def initial_stock_of(self, d) -> float: """The initial stock of device ``d``, defaulting to 0. - Device indices reaching this from a commitment's "device" column may be - numpy floats, hence the cast. + Device indices reaching this from a commitment's "device" column may be numpy floats, hence the cast. """ if isinstance(self.initial_stock, list): # No initial stock defined for inflexible device @@ -247,8 +240,8 @@ def initial_stock_of(self, d) -> float: def commodity_devices(self) -> dict: """commodity -> set(device indices). - Computed on demand: only the EMS-level flow commitment constraints need - it, and the per-row scan is not cheap enough to pay for unconditionally. + Computed on demand: only the EMS-level flow commitment constraints need it, + and the per-row scan is not cheap enough to pay for unconditionally. """ commodity_devices: dict = {} for df in self.original_commitments: @@ -278,8 +271,8 @@ def prepare_scheduling_problem( # noqa C901 ) -> SchedulingProblem: """Normalise and validate ``device_scheduler``'s arguments into a SchedulingProblem. - .. note:: This adds a "stock delta" column to the passed ``device_constraints`` - DataFrames in place, as the schedulers have always done. + .. note:: This adds a "stock delta" column to the passed ``device_constraints`` DataFrames in place, + as the schedulers have always done. """ # Get timing from first device start = device_constraints[0].index.to_pydatetime()[0] @@ -432,11 +425,10 @@ def prepare_scheduling_problem( # noqa C901 has_device_group = "device_group" in df.columns - # Read the columns as arrays rather than slicing + dropna()-ing a fresh - # DataFrame per sub-commitment. Each time step usually forms its own - # group, so this loop runs once per time step, and the per-call pandas - # overhead dominated it (~50 ms of a ~135 ms prepare on 4 devices x 192 - # steps; the arrays bring that under 1 ms). + # Read the columns as arrays rather than slicing + dropna()-ing a fresh DataFrame per sub-commitment. + # Each time step usually forms its own group, so this loop runs once per time step, + # and the per-call pandas overhead dominated it + # (~50 ms of a ~135 ms prepare on 4 devices x 192 steps; the arrays bring that under 1 ms). device_values = df["device"].to_numpy() if has_device_group: group_values = df["device_group"].to_numpy() @@ -448,8 +440,8 @@ def prepare_scheduling_problem( # noqa C901 groups: dict = {} for d, g in zip(device_values, group_values): - # Skip what the previous dropna() dropped: a missing device, or a - # missing group label when the commitment declares groups. + # Skip what the previous dropna() dropped: + # a missing device, or a missing group label when the commitment declares groups. if _is_missing(d) or (has_device_group and _is_missing(g)): continue @@ -474,7 +466,7 @@ def prepare_scheduling_problem( # noqa C901 ) else: # No commitments at all: nothing can make the cost curve non-convex. - # (The Pyomo path used to raise on the empty pd.concat here.) + # The Pyomo path used to raise on the empty pd.concat here. convex_cost_curve = True bigM_columns = ["derivative max", "derivative min", "derivative equals"] diff --git a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py index a6f88c2723..b1459cd92b 100644 --- a/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py +++ b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py @@ -1,13 +1,12 @@ """Equivalence tests for the direct highspy scheduling backend. -Each scenario is run through ``device_scheduler`` twice: once with the Pyomo -path (``appsi_highs``) and once with the direct HiGHS path (``highspy``), and -the resulting schedules, costs and termination handling are compared. - -If one of these tests fails after a change to the model in -``linear_optimization.device_scheduler``, the twin model in -``highspy_optimization.device_scheduler_highspy`` probably needs the same -change (see the note in that module's docstring). +Each scenario is run through ``device_scheduler`` twice: +once with the Pyomo path (``appsi_highs``) and once with the direct HiGHS path (``highspy``), +and the resulting schedules, costs and termination handling are compared. + +If one of these tests fails after a change to the model in ``linear_optimization.device_scheduler``, +the twin model in ``highspy_optimization.device_scheduler_highspy`` probably needs the same change +(see the note in that module's docstring). """ from __future__ import annotations @@ -259,9 +258,9 @@ def run_with_solver(app, solver: str, make_scenario): def scenario_ems_level_flow_commitment(): """Two devices under an EMS-level flow commitment, which names no device. - Such a commitment binds the summed flow of all devices, via - ``ems_flow_commitment_equalities`` rather than the grouped constraints, so it - is the case that distinguishes the two constraint families. + Such a commitment binds the summed flow of all devices, + via ``ems_flow_commitment_equalities`` rather than the grouped constraints, + so it is the case that distinguishes the two constraint families. """ index = make_index() prices = make_prices(index) @@ -284,9 +283,9 @@ def scenario_ems_level_flow_commitment(): def scenario_ems_level_commodity_commitment(): """An EMS-level flow commitment scoped to one commodity's devices. - Device 0 carries the commodity, device 1 does not, so the commitment must - bind device 0's flow only -- exercising the commodity_devices lookup rather - than the sum-over-all-devices fallback. + Device 0 carries the commodity, device 1 does not, + so the commitment must bind device 0's flow only -- + exercising the commodity_devices lookup rather than the sum-over-all-devices fallback. """ index = make_index() prices = make_prices(index) @@ -381,12 +380,12 @@ def test_highspy_matches_pyomo_when_infeasible(app): def test_unsupported_argument_is_rejected_not_ignored(): """A device_scheduler argument the direct backend cannot model must raise. - ``device_scheduler`` forwards its arguments to the direct HiGHS backend by - name. Whoever adds the next scheduling parameter (``coupling_groups`` in - #2218, ``balance_groups`` in #2289) works on the Pyomo model, and a - parameter that never reached the backend would not fail -- it would produce - a schedule computed as if the constraint had never been requested. Since - ``highspy`` is the default solver, that would be silently wrong. + ``device_scheduler`` forwards its arguments to the direct HiGHS backend by name. + Whoever adds the next scheduling parameter (``coupling_groups`` in #2218, ``balance_groups`` in #2289) + works on the Pyomo model, + and a parameter that never reached the backend would not fail -- + it would produce a schedule computed as if the constraint had never been requested. + Since ``highspy`` is the default solver, that would be silently wrong. """ real_device_scheduler = linear_optimization.device_scheduler diff --git a/flexmeasures/data/models/planning/tests/test_solver.py b/flexmeasures/data/models/planning/tests/test_solver.py index 106ca6cfbc..57cf8840b3 100644 --- a/flexmeasures/data/models/planning/tests/test_solver.py +++ b/flexmeasures/data/models/planning/tests/test_solver.py @@ -3156,10 +3156,10 @@ def initialize_combined_commitments(num_devices: int): ] # Expected results with unfair unmet demand and not entirely unfair costs. - # NB This problem has multiple optima: only the site-level (aggregate) - # schedule is unique, while the per-device allocation of the charging slots - # (and thereby the per-device costs and even which device's demand goes - # unmet) is an arbitrary tie-break that depends on the solver backend. + # NB This problem has multiple optima: only the site-level (aggregate) schedule is unique, + # while the per-device allocation of the charging slots + # (and thereby the per-device costs, and even which device's demand goes unmet) + # is an arbitrary tie-break that depends on the solver backend. # We therefore only check solver-independent properties here. expected_aggregate_schedule = [0, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25] + [0] * 17 total_expected_demand_unmet = total_expected_demand - np.array(