diff --git a/documentation/changelog.rst b/documentation/changelog.rst index bbe5a5ad80..eda7bcad36 100644 --- a/documentation/changelog.rst +++ b/documentation/changelog.rst @@ -51,6 +51,7 @@ Infrastructure / Support ---------------------- * The database migration for this release splits each stored flex-context's ``inflexible-device-sensors`` field into ``inflexible-consumption``/``inflexible-production`` sensor references, classifying each sensor by its ``consumption_is_positive`` attribute (behavior-preserving; sensor attributes themselves are kept). Downgrading merges them back into bare sensor IDs, dropping any source filters added in the meantime [see `PR #2358 `_] * Speed up listing assets: eager-load each asset's sensors instead of lazy-loading them one query per asset during serialization, and skip loading sensors entirely for field-filtered responses that do not include them [see `PR #2363 `_] +* 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 `_] diff --git a/documentation/configuration.rst b/documentation/configuration.rst index 2a76ca8570..91af658efb 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 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"`` FLEXMEASURES_LP_SOLVER_OPTIONS diff --git a/documentation/host/deployment.rst b/documentation/host/deployment.rst index c5937d16ab..c766d03d13 100644 --- a/documentation/host/deployment.rst +++ b/documentation/host/deployment.rst @@ -89,25 +89,19 @@ 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: - - .. 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/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/data/models/planning/highspy_optimization.py b/flexmeasures/data/models/planning/highspy_optimization.py new file mode 100644 index 0000000000..0dbae6eb3e --- /dev/null +++ b/flexmeasures/data/models/planning/highspy_optimization.py @@ -0,0 +1,728 @@ +"""Direct HiGHS (highspy) implementation of the device scheduler. + +.. warning:: TWO MODELS TO KEEP IN SYNC + + 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. +""" + +from __future__ import annotations + +import numpy as np +import pandas as pd + +from flexmeasures.data.models.planning import ( + Commitment, + FlowCommitment, + StockCommitment, +) +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") + + +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 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 = _SolverInformation(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 + + 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, + ) + + # 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) + # --------------------------------------------------------------- + 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") + + # --------------------------------------------------------------- + # 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). + 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) + active = ~(np.isnan(quantity) | (quantity == -infinity)) + lb = 0.0 if "upwards deviation price" in df.columns else -infinity + ub = 0.0 if "downwards deviation price" in df.columns else infinity + 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 + _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: + 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. + 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.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) + + # 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() + + 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] + + planned_power = planned_power_per_device(ems_values, start, end, resolution) + + 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) + 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, planned_costs, results, model diff --git a/flexmeasures/data/models/planning/linear_optimization.py b/flexmeasures/data/models/planning/linear_optimization.py index b42c9924e2..51b94a0957 100644 --- a/flexmeasures/data/models/planning/linear_optimization.py +++ b/flexmeasures/data/models/planning/linear_optimization.py @@ -1,11 +1,11 @@ from __future__ import annotations -import math +import inspect +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, @@ -29,48 +29,99 @@ 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. +def _left_at_default(value, default) -> bool: + """Whether an argument was left at its default. - 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. + 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: - 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." - ) + 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, + ) - 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." + 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 @@ -143,301 +194,55 @@ 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": + # 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(**highspy_arguments) + model = ConcreteModel() # If the EMS has no devices, don't bother 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) - ) - - 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 = {} - - 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") @@ -669,12 +474,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( @@ -684,19 +483,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. @@ -708,7 +494,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) @@ -780,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 @@ -797,7 +582,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 @@ -933,27 +718,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. @@ -965,44 +732,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..92b48d349d --- /dev/null +++ b/flexmeasures/data/models/planning/scheduling_problem.py @@ -0,0 +1,573 @@ +"""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 _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"`. + + 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 + + # 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() + 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). + group_values = device_values + + 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} + + groups.setdefault(g, set()).update(devices) + + device_group_lookup[c] = groups + + # 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 + ] 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..b1459cd92b --- /dev/null +++ b/flexmeasures/data/models/planning/tests/test_highspy_equivalence.py @@ -0,0 +1,422 @@ +"""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 + +import inspect +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 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 + +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 + + +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", + [ + scenario_battery_with_prices, + 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_", ""), +) +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 + + +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 diff --git a/flexmeasures/data/models/planning/tests/test_solver.py b/flexmeasures/data/models/planning/tests/test_solver.py index e1a4ecdd31..57cf8840b3 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): diff --git a/flexmeasures/utils/config_defaults.py b/flexmeasures/utils/config_defaults.py index f564e2db2b..f3828b87cd 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] = {}