From 1cb49e394bfebf7fde2c20c1e906ae4e82799d01 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Tue, 4 Aug 2026 09:20:09 +0200 Subject: [PATCH 1/4] fix(optimizer): guard storage starts above maximum --- .changeset/fix-optimizer-soc-above-max.md | 5 +++ optimizer/ftw_optimizer/direct_highs.py | 15 ++++++++ optimizer/ftw_optimizer/model.py | 3 +- optimizer/ftw_optimizer/multistage.py | 11 +++++- optimizer/ftw_optimizer/recourse.py | 3 +- optimizer/tests/test_model.py | 45 +++++++++++++++++++---- 6 files changed, 71 insertions(+), 11 deletions(-) create mode 100644 .changeset/fix-optimizer-soc-above-max.md diff --git a/.changeset/fix-optimizer-soc-above-max.md b/.changeset/fix-optimizer-soc-above-max.md new file mode 100644 index 00000000..5b40a9f4 --- /dev/null +++ b/.changeset/fix-optimizer-soc-above-max.md @@ -0,0 +1,5 @@ +--- +"ftw": patch +--- + +Keep optimizer storage plans replay-consistent when a battery starts above its configured maximum state of charge. diff --git a/optimizer/ftw_optimizer/direct_highs.py b/optimizer/ftw_optimizer/direct_highs.py index f8316a13..43ea1efc 100644 --- a/optimizer/ftw_optimizer/direct_highs.py +++ b/optimizer/ftw_optimizer/direct_highs.py @@ -98,6 +98,14 @@ def solve_direct_highs( prepare_ms: float, decomposition: str, ) -> dict[str, Any]: + if any( + float(spec["initial_energy_wh"]) + > float(spec.get("max_energy_wh", spec["capacity_wh"])) + 1e-6 + for spec in prepared.storages + ): + raise DirectHighsError( + "direct HiGHS path requires storage starts at or below the operating maximum" + ) if prepared.discrete or prepared.unsafe_cycle or prepared.unsafe_meter_split: raise DirectHighsError("direct HiGHS path requires a cycle-safe continuous tariff") build_started = time.perf_counter() @@ -419,6 +427,13 @@ def _response( base_index = next((i for i, scenario in enumerate(scenarios) if scenario.id == "base"), 0) base = scenarios[base_index] base_vars = scenario_vars[base_index] + for scenario in scenario_vars: + for charges, discharges in zip(scenario.charge, scenario.discharge): + for charge_index, discharge_index in zip(charges, discharges): + if min(solution[charge_index], solution[discharge_index]) > 1e-6: + raise DirectHighsError( + "HiGHS returned simultaneous storage charge and discharge" + ) total_capacity = sum(float(spec["capacity_wh"]) for spec in prepared.storages) initial_total = sum(float(spec["initial_energy_wh"]) for spec in prepared.storages) actions: list[dict[str, Any]] = [] diff --git a/optimizer/ftw_optimizer/model.py b/optimizer/ftw_optimizer/model.py index db419eb1..2319df22 100644 --- a/optimizer/ftw_optimizer/model.py +++ b/optimizer/ftw_optimizer/model.py @@ -315,7 +315,8 @@ def solve(payload: dict[str, Any]) -> dict[str, Any]: # starts have zero recovery allowance, preserving hard min/max bounds. service_slack += cp.sum(lower_recovery[1:] + upper_recovery[1:]) / (capacity * n) unsafe_cycle = bool(np.any(eff_import < 0)) or pv_charge_bonus_ore > 0 - if force_milp or (formulation == "auto" and unsafe_cycle): + initial_above_max = initial > max_energy + 1e-6 + if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): direction = cp.Variable(n, boolean=True, name=f"storage_{i}_charge_mode") constraints += [charge <= max_charge * direction, discharge <= max_discharge * (1 - direction)] discrete = True diff --git a/optimizer/ftw_optimizer/multistage.py b/optimizer/ftw_optimizer/multistage.py index 64df8ecb..cd9fcf72 100644 --- a/optimizer/ftw_optimizer/multistage.py +++ b/optimizer/ftw_optimizer/multistage.py @@ -353,7 +353,6 @@ def _prepare(payload: dict[str, Any]) -> PreparedMultistage: unsafe_meter_split = bool( np.any(effective_import < effective_export - 1e-9) ) - storage_discrete = force_milp or (formulation == "auto" and unsafe_cycle) meter_discrete = force_milp or ( formulation == "auto" and unsafe_meter_split ) @@ -410,6 +409,16 @@ def _prepare(payload: dict[str, Any]) -> PreparedMultistage: _validate_storages(storage_specs, n) if not storage_specs: raise ProtocolError("multistage shadow requires at least one storage") + storage_above_max = any( + float(spec["initial_energy_wh"]) + > float(spec.get("max_energy_wh", spec["capacity_wh"])) + 1e-6 + for spec in storage_specs + ) + storage_discrete = ( + force_milp + or storage_above_max + or (formulation == "auto" and unsafe_cycle) + ) max_site_power = max( 1000.0, diff --git a/optimizer/ftw_optimizer/recourse.py b/optimizer/ftw_optimizer/recourse.py index ef328b54..d2e732ed 100644 --- a/optimizer/ftw_optimizer/recourse.py +++ b/optimizer/ftw_optimizer/recourse.py @@ -177,7 +177,8 @@ def solve_storage_recourse(payload: dict[str, Any]) -> dict[str, Any]: upper_recovery[1:] <= upper_recovery[:-1], ] scenario_service += cp.sum(lower_recovery[1:] + upper_recovery[1:]) / (capacity * n) - if force_milp or (formulation == "auto" and unsafe_cycle): + initial_above_max = initial > max_energy + 1e-6 + if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): direction = cp.Variable(n, boolean=True, name=f"scenario_{si}_storage_{i}_charge_mode") constraints += [charge <= max_charge * direction, discharge <= max_discharge * (1 - direction)] discrete = True diff --git a/optimizer/tests/test_model.py b/optimizer/tests/test_model.py index 7a153409..8b53147c 100644 --- a/optimizer/tests/test_model.py +++ b/optimizer/tests/test_model.py @@ -722,11 +722,40 @@ def test_storage_below_minimum_recovers_without_worsening() -> None: assert energies[-1] >= 1000 - 0.01 -def test_storage_above_maximum_recovers_without_worsening() -> None: - request = base_request() - request["storages"][0]["initial_energy_wh"] = 9800 - response = handle(request) - assert response["ok"], response - energies = [action["storage_energy_wh"]["home"] for action in response["plan"]["actions"]] - assert energies[0] <= 9800 + 1e-5 - assert energies[-1] <= 9500 + 1e-5 +def test_storage_above_maximum_replays_without_simultaneous_energy_loss() -> None: + expected_formulations = { + "shared": "milp", + "recourse": "stochastic-recourse-milp", + "multistage": "multistage-milp", + } + for scenario_policy in ("shared", "recourse", "multistage"): + for formulation in ("auto", "relaxed"): + request = base_request() + request["request_id"] = f"soc-above-max-{scenario_policy}-{formulation}" + request["settings"]["formulation"] = formulation + if scenario_policy != "shared": + request["settings"]["scenario_policy"] = scenario_policy + request["storages"][0]["initial_energy_wh"] = 9800 + + response = handle(request) + + assert response["ok"], response + assert response["solver"]["formulation"] == expected_formulations[scenario_policy] + energy_wh = request["storages"][0]["initial_energy_wh"] + storage = request["storages"][0] + for slot, action in zip(request["slots"], response["plan"]["actions"]): + power_w = action["storage_power_w"]["home"] + previous_energy_wh = energy_wh + dt_h = slot["len_min"] / 60.0 + if power_w >= 0: + energy_wh += power_w * dt_h * storage["charge_efficiency"] + else: + energy_wh += power_w * dt_h / storage["discharge_efficiency"] + assert math.isclose( + action["storage_energy_wh"]["home"], + energy_wh, + abs_tol=0.1, + ) + if abs(power_w) <= 1e-6: + assert action["storage_energy_wh"]["home"] >= previous_energy_wh - 0.1 + assert energy_wh <= storage["max_energy_wh"] + 0.1 From db7a7afb69ce198e856929973fc39bb07e7282f9 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Tue, 4 Aug 2026 09:57:35 +0200 Subject: [PATCH 2/4] fix(optimizer): make storage replay guard backend independent --- optimizer/ftw_optimizer/direct_highs.py | 8 +- optimizer/ftw_optimizer/model.py | 76 ++++++++++- optimizer/ftw_optimizer/multistage.py | 160 +++++++++++++++--------- optimizer/ftw_optimizer/recourse.py | 13 +- optimizer/tests/test_model.py | 82 ++++++++++++ 5 files changed, 268 insertions(+), 71 deletions(-) diff --git a/optimizer/ftw_optimizer/direct_highs.py b/optimizer/ftw_optimizer/direct_highs.py index 43ea1efc..b96a8fb5 100644 --- a/optimizer/ftw_optimizer/direct_highs.py +++ b/optimizer/ftw_optimizer/direct_highs.py @@ -9,7 +9,7 @@ import numpy as np from . import SCHEMA_VERSION -from .model import _solver_options +from .model import _solver_options, _storage_starts_above_maximum from .protocol import finite_number if TYPE_CHECKING: @@ -98,11 +98,7 @@ def solve_direct_highs( prepare_ms: float, decomposition: str, ) -> dict[str, Any]: - if any( - float(spec["initial_energy_wh"]) - > float(spec.get("max_energy_wh", spec["capacity_wh"])) + 1e-6 - for spec in prepared.storages - ): + if _storage_starts_above_maximum(prepared.storages): raise DirectHighsError( "direct HiGHS path requires storage starts at or below the operating maximum" ) diff --git a/optimizer/ftw_optimizer/model.py b/optimizer/ftw_optimizer/model.py index 2319df22..9fa236d2 100644 --- a/optimizer/ftw_optimizer/model.py +++ b/optimizer/ftw_optimizer/model.py @@ -41,6 +41,79 @@ class ThermalVars: upper_slack: cp.Variable +class ReplayConsistencyError(RuntimeError): + """The reported storage state cannot be replayed from reported power.""" + + +_REPLAY_TOLERANCE_FRACTION = 0.0002 +_REPLAY_TOLERANCE_MIN_WH = 1.0 + + +def _storage_starts_above_maximum(storages: Any) -> bool: + return any( + float(spec["initial_energy_wh"]) + > float(spec.get("max_energy_wh", spec["capacity_wh"])) + for spec in storages + ) + + +def _storage_replay_tolerance_wh(spec: dict[str, Any]) -> float: + return max( + _REPLAY_TOLERANCE_MIN_WH, + float(spec["capacity_wh"]) * _REPLAY_TOLERANCE_FRACTION, + ) + + +def _validate_storage_replay( + actions: list[dict[str, Any]], + slots: list[dict[str, Any]] | tuple[dict[str, Any], ...], + storages: list[dict[str, Any]] | tuple[dict[str, Any], ...], +) -> None: + if len(actions) != len(slots): + raise ReplayConsistencyError( + f"action count {len(actions)} does not match slot count {len(slots)}" + ) + energy = { + str(spec["id"]): float(spec["initial_energy_wh"]) for spec in storages + } + for slot_index, (slot, action) in enumerate(zip(slots, actions)): + dt_h = float(slot["len_min"]) / 60.0 + storage_power = action.get("storage_power_w", {}) + storage_energy = action.get("storage_energy_wh", {}) + for spec in storages: + storage_id = str(spec["id"]) + if storage_id not in storage_power or storage_id not in storage_energy: + raise ReplayConsistencyError( + f"slot {slot_index} storage {storage_id} output is missing" + ) + power = float(storage_power[storage_id]) + reported = float(storage_energy[storage_id]) + if not math.isfinite(power) or not math.isfinite(reported): + raise ReplayConsistencyError( + f"slot {slot_index} storage {storage_id} output is non-finite" + ) + if power >= 0: + replayed = energy[storage_id] + power * dt_h * float( + spec.get("charge_efficiency", 0.95) + ) + else: + replayed = energy[storage_id] + power * dt_h / float( + spec.get("discharge_efficiency", 0.95) + ) + tolerance = _storage_replay_tolerance_wh(spec) + if ( + replayed < -tolerance + or replayed > float(spec["capacity_wh"]) + tolerance + or abs(reported - replayed) > tolerance + ): + raise ReplayConsistencyError( + f"slot {slot_index} storage {storage_id} energy " + f"{reported:.6f} is inconsistent with replay " + f"{replayed:.6f} (tolerance {tolerance:.6f})" + ) + energy[storage_id] = replayed + + def _vector(value: Any, n: int, field: str) -> np.ndarray: items = require_list(value, field) if len(items) != n: @@ -315,7 +388,7 @@ def solve(payload: dict[str, Any]) -> dict[str, Any]: # starts have zero recovery allowance, preserving hard min/max bounds. service_slack += cp.sum(lower_recovery[1:] + upper_recovery[1:]) / (capacity * n) unsafe_cycle = bool(np.any(eff_import < 0)) or pv_charge_bonus_ore > 0 - initial_above_max = initial > max_energy + 1e-6 + initial_above_max = _storage_starts_above_maximum([spec]) if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): direction = cp.Variable(n, boolean=True, name=f"storage_{i}_charge_mode") constraints += [charge <= max_charge * direction, discharge <= max_discharge * (1 - direction)] @@ -818,6 +891,7 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: mip_gap = float(value) break solve_ms = (time.perf_counter() - started) * 1000.0 + _validate_storage_replay(actions, slots, [storage.spec for storage in storages]) return { "schema_version": SCHEMA_VERSION, "request_id": str(payload["request_id"]), diff --git a/optimizer/ftw_optimizer/multistage.py b/optimizer/ftw_optimizer/multistage.py index cd9fcf72..b4bd6e18 100644 --- a/optimizer/ftw_optimizer/multistage.py +++ b/optimizer/ftw_optimizer/multistage.py @@ -4,14 +4,22 @@ import math import time from collections import OrderedDict -from dataclasses import dataclass +from dataclasses import dataclass, replace from typing import Any import cvxpy as cp import numpy as np from . import SCHEMA_VERSION -from .model import OPTIMAL_STATUSES, _export_price, _mode, _solver_options +from .model import ( + OPTIMAL_STATUSES, + ReplayConsistencyError, + _export_price, + _mode, + _solver_options, + _storage_starts_above_maximum, + _validate_storage_replay, +) from .protocol import ProtocolError, finite_number, positive_number, require_dict, require_list from .scenario_tree import ( ScenarioSet, @@ -191,11 +199,20 @@ def solve_storage_multistage(payload: dict[str, Any]) -> dict[str, Any]: eligible, reason = ph_eligible(prepared) if eligible: try: - return solve_progressive_hedging(prepared, started, prepare_ms) + response = solve_progressive_hedging(prepared, started, prepare_ms) + _validate_storage_replay( + response["plan"]["actions"], prepared.slots, prepared.storages + ) + return response except ProgressiveHedgingNotConverged: if decomposition_method == "progressive_hedging": raise decomposition = "ph-fallback-scenario-reduction-extensive-dpp" + except ReplayConsistencyError as exc: + if decomposition_method == "progressive_hedging": + raise + prepared = _with_storage_discrete(prepared) + decomposition = f"ph-fallback-storage-replay-{exc}" elif decomposition_method == "progressive_hedging": raise ProtocolError(f"progressive hedging is not eligible: {reason}") @@ -239,69 +256,86 @@ def solve_storage_multistage(payload: dict[str, Any]) -> dict[str, Any]: from .direct_highs import DirectHighsError, solve_direct_highs try: - return solve_direct_highs( + response = solve_direct_highs( prepared, started, prepare_ms, decomposition.replace("-dpp", "") ) - except DirectHighsError as exc: + _validate_storage_replay( + response["plan"]["actions"], prepared.slots, prepared.storages + ) + return response + except (DirectHighsError, ReplayConsistencyError) as exc: if multistage_backend == "highs": raise direct_fallback_reason = str(exc) decomposition = f"direct-highs-fallback-{decomposition}" + prepared = _with_storage_discrete(prepared) - key = _cache_key(prepared) - compiled = _MODEL_CACHE.get(key) - cache_hit = compiled is not None - if compiled is None: - compiled = _compile(prepared, key) - _MODEL_CACHE[key] = compiled - while len(_MODEL_CACHE) > _CACHE_LIMIT: - _MODEL_CACHE.popitem(last=False) - else: - _MODEL_CACHE.move_to_end(key) - compiled.assign(prepared) + while True: + key = _cache_key(prepared) + compiled = _MODEL_CACHE.get(key) + cache_hit = compiled is not None + if compiled is None: + compiled = _compile(prepared, key) + _MODEL_CACHE[key] = compiled + while len(_MODEL_CACHE) > _CACHE_LIMIT: + _MODEL_CACHE.popitem(last=False) + else: + _MODEL_CACHE.move_to_end(key) + compiled.assign(prepared) - if prepared.discrete and solver_name == "CLARABEL": - solver_name = "HIGHS" - solver_started = time.perf_counter() - try: - _run_problem(compiled.service_problem, prepared.settings, solver_name) - except cp.error.SolverError: - if prepared.discrete or solver_name == "CLARABEL": - raise - solver_name = "CLARABEL" - _run_problem(compiled.service_problem, prepared.settings, solver_name) - if compiled.service_problem.status not in OPTIMAL_STATUSES or compiled.service_problem.value is None: - raise RuntimeError( - f"multistage service-level solve failed with status {compiled.service_problem.status}" - ) - best_service = max(0.0, float(compiled.service_problem.value)) - compiled.service_cap.value = best_service + 1e-7 - try: - _run_problem(compiled.economic_problem, prepared.settings, solver_name) - except cp.error.SolverError: - if prepared.discrete or solver_name == "CLARABEL": - raise - solver_name = "CLARABEL" - _run_problem(compiled.economic_problem, prepared.settings, solver_name) - if compiled.economic_problem.status not in OPTIMAL_STATUSES or compiled.economic_problem.value is None: - raise RuntimeError( - f"multistage economic solve failed with status {compiled.economic_problem.status}" - ) - solver_ms = (time.perf_counter() - solver_started) * 1000.0 + if prepared.discrete and solver_name == "CLARABEL": + solver_name = "HIGHS" + solver_started = time.perf_counter() + try: + _run_problem(compiled.service_problem, prepared.settings, solver_name) + except cp.error.SolverError: + if prepared.discrete or solver_name == "CLARABEL": + raise + solver_name = "CLARABEL" + _run_problem(compiled.service_problem, prepared.settings, solver_name) + if compiled.service_problem.status not in OPTIMAL_STATUSES or compiled.service_problem.value is None: + raise RuntimeError( + f"multistage service-level solve failed with status {compiled.service_problem.status}" + ) + best_service = max(0.0, float(compiled.service_problem.value)) + compiled.service_cap.value = best_service + 1e-7 + try: + _run_problem(compiled.economic_problem, prepared.settings, solver_name) + except cp.error.SolverError: + if prepared.discrete or solver_name == "CLARABEL": + raise + solver_name = "CLARABEL" + _run_problem(compiled.economic_problem, prepared.settings, solver_name) + if compiled.economic_problem.status not in OPTIMAL_STATUSES or compiled.economic_problem.value is None: + raise RuntimeError( + f"multistage economic solve failed with status {compiled.economic_problem.status}" + ) + solver_ms = (time.perf_counter() - solver_started) * 1000.0 - response = _response( - prepared, - compiled, - best_service, - solver_name, - started, - prepare_ms, - solver_ms, - cache_hit, - decomposition, - direct_fallback_reason, - ) - return response + response = _response( + prepared, + compiled, + best_service, + solver_name, + started, + prepare_ms, + solver_ms, + cache_hit, + decomposition, + direct_fallback_reason, + ) + try: + _validate_storage_replay( + response["plan"]["actions"], prepared.slots, prepared.storages + ) + except ReplayConsistencyError as exc: + if prepared.storage_discrete: + raise + prepared = _with_storage_discrete(prepared) + direct_fallback_reason = str(exc) + decomposition = f"storage-replay-fallback-{decomposition}" + continue + return response def clear_multistage_cache() -> None: @@ -409,11 +443,7 @@ def _prepare(payload: dict[str, Any]) -> PreparedMultistage: _validate_storages(storage_specs, n) if not storage_specs: raise ProtocolError("multistage shadow requires at least one storage") - storage_above_max = any( - float(spec["initial_energy_wh"]) - > float(spec.get("max_energy_wh", spec["capacity_wh"])) + 1e-6 - for spec in storage_specs - ) + storage_above_max = _storage_starts_above_maximum(storage_specs) storage_discrete = ( force_milp or storage_above_max @@ -538,6 +568,12 @@ def _replace_scenarios(prepared: PreparedMultistage, scenario_set: ScenarioSet) return PreparedMultistage(**{**prepared.__dict__, "scenario_set": scenario_set, "tree": tree, "blocks": blocks}) +def _with_storage_discrete(prepared: PreparedMultistage) -> PreparedMultistage: + if prepared.storage_discrete: + return prepared + return replace(prepared, storage_discrete=True, discrete=True) + + def _validate_storages(storages: tuple[dict[str, Any], ...], n: int) -> None: ids: set[str] = set() for i, spec in enumerate(storages): diff --git a/optimizer/ftw_optimizer/recourse.py b/optimizer/ftw_optimizer/recourse.py index d2e732ed..8d87a37a 100644 --- a/optimizer/ftw_optimizer/recourse.py +++ b/optimizer/ftw_optimizer/recourse.py @@ -9,7 +9,15 @@ import numpy as np from . import SCHEMA_VERSION -from .model import OPTIMAL_STATUSES, _export_price, _mode, _solver_options, _vector +from .model import ( + OPTIMAL_STATUSES, + _export_price, + _mode, + _solver_options, + _storage_starts_above_maximum, + _validate_storage_replay, + _vector, +) from .protocol import ProtocolError, finite_number, positive_number, require_dict, require_list @@ -177,7 +185,7 @@ def solve_storage_recourse(payload: dict[str, Any]) -> dict[str, Any]: upper_recovery[1:] <= upper_recovery[:-1], ] scenario_service += cp.sum(lower_recovery[1:] + upper_recovery[1:]) / (capacity * n) - initial_above_max = initial > max_energy + 1e-6 + initial_above_max = _storage_starts_above_maximum([spec]) if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): direction = cp.Variable(n, boolean=True, name=f"scenario_{si}_storage_{i}_charge_mode") constraints += [charge <= max_charge * direction, discharge <= max_discharge * (1 - direction)] @@ -371,6 +379,7 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: mip_gap = float(value) break solve_ms = (time.perf_counter() - started) * 1000.0 + _validate_storage_replay(actions, slots, storage_specs) return { "schema_version": SCHEMA_VERSION, "request_id": str(payload["request_id"]), diff --git a/optimizer/tests/test_model.py b/optimizer/tests/test_model.py index 8b53147c..31b62984 100644 --- a/optimizer/tests/test_model.py +++ b/optimizer/tests/test_model.py @@ -150,6 +150,28 @@ def base_request() -> dict: } +def assert_storage_replays(request: dict, response: dict, tolerance_wh: float = 2.1) -> None: + energies = { + str(spec["id"]): float(spec["initial_energy_wh"]) + for spec in request["storages"] + } + for slot, action in zip(request["slots"], response["plan"]["actions"]): + dt_h = slot["len_min"] / 60.0 + for spec in request["storages"]: + storage_id = str(spec["id"]) + power = action["storage_power_w"][storage_id] + previous = energies[storage_id] + if power >= 0: + replayed = previous + power * dt_h * spec["charge_efficiency"] + else: + replayed = previous + power * dt_h / spec["discharge_efficiency"] + reported = action["storage_energy_wh"][storage_id] + assert math.isclose(reported, replayed, abs_tol=tolerance_wh) + if abs(power) <= 1e-3: + assert reported >= previous - tolerance_wh + energies[storage_id] = replayed + + def test_arbitrage_moves_energy_from_cheap_to_expensive_slot() -> None: response = handle(base_request()) assert response["ok"], response @@ -759,3 +781,63 @@ def test_storage_above_maximum_replays_without_simultaneous_energy_loss() -> Non if abs(power_w) <= 1e-6: assert action["storage_energy_wh"]["home"] >= previous_energy_wh - 0.1 assert energy_wh <= storage["max_energy_wh"] + 0.1 + + +def test_storage_just_above_maximum_uses_replay_safe_guard() -> None: + expected_formulations = { + "shared": "milp", + "recourse": "stochastic-recourse-milp", + "multistage": "multistage-milp", + } + for delta in (0.5e-6, 1e-6): + for scenario_policy in ("shared", "recourse", "multistage"): + request = base_request() + request["request_id"] = f"soc-just-above-max-{scenario_policy}-{delta}" + request["settings"].update( + {"mode": "cheap_charge", "formulation": "relaxed"} + ) + if scenario_policy != "shared": + request["settings"]["scenario_policy"] = scenario_policy + request["storages"][0]["initial_energy_wh"] = 9500 + delta + + response = handle(request) + + assert response["ok"], response + assert response["solver"]["formulation"] == expected_formulations[scenario_policy] + assert_storage_replays(request, response) + + +def test_multistage_auto_retries_with_storage_guard_after_direct_cycle(monkeypatch) -> None: + from ftw_optimizer import direct_highs + + clear_multistage_cache() + request = base_request() + request["settings"].update( + { + "mode": "passive_arbitrage", + "formulation": "relaxed", + "scenario_policy": "multistage", + "multistage_backend": "auto", + } + ) + request["storages"][0]["initial_energy_wh"] = 2000 + + def reject_simultaneous_cycle(*args, **kwargs): + raise direct_highs.DirectHighsError( + "HiGHS returned simultaneous storage charge and discharge" + ) + + monkeypatch.setattr(direct_highs, "solve_direct_highs", reject_simultaneous_cycle) + response = handle(request) + + assert response["ok"], response + assert response["solver"]["formulation"] == "multistage-milp" + assert "simultaneous" in response["solver"]["fallback_reason"] + assert_storage_replays(request, response) + + cvxpy_request = copy.deepcopy(request) + cvxpy_request["request_id"] = "multistage-cvxpy-replay-guard" + cvxpy_request["settings"]["multistage_backend"] = "cvxpy" + cvxpy_response = handle(cvxpy_request) + assert cvxpy_response["ok"], cvxpy_response + assert_storage_replays(cvxpy_request, cvxpy_response) From dcc877acd7bac6b74370d6d1024ac0b0f1987e03 Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Tue, 4 Aug 2026 10:43:03 +0200 Subject: [PATCH 3/4] fix(optimizer): normalize tiny SoC bound noise --- .changeset/fix-optimizer-soc-above-max.md | 2 +- optimizer/ftw_optimizer/model.py | 58 ++++++++++++++++- optimizer/ftw_optimizer/multistage.py | 9 ++- optimizer/ftw_optimizer/recourse.py | 11 ++-- optimizer/tests/test_model.py | 78 +++++++++++++++++++++++ 5 files changed, 143 insertions(+), 15 deletions(-) diff --git a/.changeset/fix-optimizer-soc-above-max.md b/.changeset/fix-optimizer-soc-above-max.md index 5b40a9f4..20e93c48 100644 --- a/.changeset/fix-optimizer-soc-above-max.md +++ b/.changeset/fix-optimizer-soc-above-max.md @@ -2,4 +2,4 @@ "ftw": patch --- -Keep optimizer storage plans replay-consistent when a battery starts above its configured maximum state of charge. +Keep optimizer storage plans replay-consistent when a battery starts above its configured maximum state of charge, including solver-scale numeric noise at that boundary. diff --git a/optimizer/ftw_optimizer/model.py b/optimizer/ftw_optimizer/model.py index 9fa236d2..493af93e 100644 --- a/optimizer/ftw_optimizer/model.py +++ b/optimizer/ftw_optimizer/model.py @@ -47,6 +47,7 @@ class ReplayConsistencyError(RuntimeError): _REPLAY_TOLERANCE_FRACTION = 0.0002 _REPLAY_TOLERANCE_MIN_WH = 1.0 +_STORAGE_NUMERIC_TOLERANCE_WH = 1e-6 def _storage_starts_above_maximum(storages: Any) -> bool: @@ -57,6 +58,55 @@ def _storage_starts_above_maximum(storages: Any) -> bool: ) +def _within_storage_numeric_tolerance(value: float, bound: float) -> bool: + """Accept one micro-Wh of decimal input plus float representation error.""" + return value - bound <= _STORAGE_NUMERIC_TOLERANCE_WH + max( + math.ulp(value), math.ulp(bound) + ) + + +def _normalize_storage_specs( + storages: Any, +) -> tuple[tuple[dict[str, Any], ...], tuple[bool, ...]]: + """Clamp only solver-scale bound noise and retain the original guard signal.""" + normalized: list[dict[str, Any]] = [] + starts_above_maximum: list[bool] = [] + for i, raw in enumerate(storages): + spec = require_dict(raw, f"storages[{i}]") + capacity = positive_number(spec.get("capacity_wh"), f"storages[{i}].capacity_wh") + minimum = finite_number(spec.get("min_energy_wh", 0), f"storages[{i}].min_energy_wh") + maximum = finite_number( + spec.get("max_energy_wh", capacity), f"storages[{i}].max_energy_wh" + ) + initial = finite_number(spec.get("initial_energy_wh"), f"storages[{i}].initial_energy_wh") + if not ( + 0 <= minimum <= maximum <= capacity + _STORAGE_NUMERIC_TOLERANCE_WH + and 0 <= initial <= capacity + _STORAGE_NUMERIC_TOLERANCE_WH + ): + raise ProtocolError(f"storages[{i}] energy bounds are inconsistent") + + model_maximum = maximum + if model_maximum > capacity and _within_storage_numeric_tolerance( + model_maximum, capacity + ): + model_maximum = capacity + initial_above_maximum = initial > model_maximum + model_initial = initial + if initial_above_maximum and _within_storage_numeric_tolerance( + initial, model_maximum + ): + model_initial = model_maximum + + model_spec = dict(spec) + if model_maximum != maximum: + model_spec["max_energy_wh"] = model_maximum + if model_initial != initial: + model_spec["initial_energy_wh"] = model_initial + normalized.append(model_spec) + starts_above_maximum.append(initial_above_maximum) + return tuple(normalized), tuple(starts_above_maximum) + + def _storage_replay_tolerance_wh(spec: dict[str, Any]) -> float: return max( _REPLAY_TOLERANCE_MIN_WH, @@ -337,13 +387,15 @@ def solve(payload: dict[str, Any]) -> dict[str, Any]: constraints: list[cp.Constraint] = [] discrete = False + storage_specs, storage_above_maximum = _normalize_storage_specs( + require_list(payload.get("storages", []), "storages") + ) storages: list[StorageVars] = [] asset_ids: set[str] = set() total_charge: cp.Expression = cp.Constant(np.zeros(n)) total_discharge: cp.Expression = cp.Constant(np.zeros(n)) service_slack: cp.Expression = cp.Constant(0.0) - for i, raw in enumerate(require_list(payload.get("storages", []), "storages")): - spec = require_dict(raw, f"storages[{i}]") + for i, spec in enumerate(storage_specs): asset_id = spec.get("id") if not isinstance(asset_id, str) or not asset_id or asset_id in asset_ids: raise ProtocolError(f"storages[{i}].id must be non-empty and unique") @@ -388,7 +440,7 @@ def solve(payload: dict[str, Any]) -> dict[str, Any]: # starts have zero recovery allowance, preserving hard min/max bounds. service_slack += cp.sum(lower_recovery[1:] + upper_recovery[1:]) / (capacity * n) unsafe_cycle = bool(np.any(eff_import < 0)) or pv_charge_bonus_ore > 0 - initial_above_max = _storage_starts_above_maximum([spec]) + initial_above_max = storage_above_maximum[i] if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): direction = cp.Variable(n, boolean=True, name=f"storage_{i}_charge_mode") constraints += [charge <= max_charge * direction, discharge <= max_discharge * (1 - direction)] diff --git a/optimizer/ftw_optimizer/multistage.py b/optimizer/ftw_optimizer/multistage.py index b4bd6e18..1d667b32 100644 --- a/optimizer/ftw_optimizer/multistage.py +++ b/optimizer/ftw_optimizer/multistage.py @@ -17,7 +17,7 @@ _export_price, _mode, _solver_options, - _storage_starts_above_maximum, + _normalize_storage_specs, _validate_storage_replay, ) from .protocol import ProtocolError, finite_number, positive_number, require_dict, require_list @@ -436,14 +436,13 @@ def _prepare(payload: dict[str, Any]) -> PreparedMultistage: tree.branch_slots, ) - storage_specs = tuple( - require_dict(raw, f"storages[{i}]") - for i, raw in enumerate(require_list(payload.get("storages", []), "storages")) + storage_specs, storage_above_maximum = _normalize_storage_specs( + require_list(payload.get("storages", []), "storages") ) _validate_storages(storage_specs, n) if not storage_specs: raise ProtocolError("multistage shadow requires at least one storage") - storage_above_max = _storage_starts_above_maximum(storage_specs) + storage_above_max = any(storage_above_maximum) storage_discrete = ( force_milp or storage_above_max diff --git a/optimizer/ftw_optimizer/recourse.py b/optimizer/ftw_optimizer/recourse.py index 8d87a37a..df33d9c7 100644 --- a/optimizer/ftw_optimizer/recourse.py +++ b/optimizer/ftw_optimizer/recourse.py @@ -14,7 +14,7 @@ _export_price, _mode, _solver_options, - _storage_starts_above_maximum, + _normalize_storage_specs, _validate_storage_replay, _vector, ) @@ -108,10 +108,9 @@ def solve_storage_recourse(payload: dict[str, Any]) -> dict[str, Any]: force_milp = formulation == "milp" constraints: list[cp.Constraint] = [] discrete = False - storage_specs = [ - require_dict(raw, f"storages[{i}]") - for i, raw in enumerate(require_list(payload.get("storages", []), "storages")) - ] + storage_specs, storage_above_maximum = _normalize_storage_specs( + require_list(payload.get("storages", []), "storages") + ) asset_ids: set[str] = set() for i, spec in enumerate(storage_specs): asset_id = spec.get("id") @@ -185,7 +184,7 @@ def solve_storage_recourse(payload: dict[str, Any]) -> dict[str, Any]: upper_recovery[1:] <= upper_recovery[:-1], ] scenario_service += cp.sum(lower_recovery[1:] + upper_recovery[1:]) / (capacity * n) - initial_above_max = _storage_starts_above_maximum([spec]) + initial_above_max = storage_above_maximum[i] if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): direction = cp.Variable(n, boolean=True, name=f"scenario_{si}_storage_{i}_charge_mode") constraints += [charge <= max_charge * direction, discharge <= max_discharge * (1 - direction)] diff --git a/optimizer/tests/test_model.py b/optimizer/tests/test_model.py index 31b62984..cca2bb40 100644 --- a/optimizer/tests/test_model.py +++ b/optimizer/tests/test_model.py @@ -4,6 +4,7 @@ import math import numpy as np +import pytest from ftw_optimizer.multistage import clear_multistage_cache from ftw_optimizer.scenario_tree import ( @@ -807,6 +808,83 @@ def test_storage_just_above_maximum_uses_replay_safe_guard() -> None: assert_storage_replays(request, response) +def _multistage_soc_boundary_request(initial_energy_wh: float, backend: str) -> dict: + request = base_request() + request["settings"].update( + { + "mode": "cheap_charge", + "formulation": "relaxed", + "scenario_policy": "multistage", + "multistage_backend": backend, + } + ) + request["slots"] = [ + { + "start_ms": 1, + "len_min": 60, + "price_ore": 20, + "spot_ore": 10, + "confidence": 1, + "pv_w": 0, + "load_w": 0, + }, + { + "start_ms": 3600001, + "len_min": 60, + "price_ore": 20, + "spot_ore": 10, + "confidence": 1, + "pv_w": 0, + "load_w": 0, + }, + ] + request["storages"][0].update( + { + "capacity_wh": 10000, + "min_energy_wh": 1000, + "max_energy_wh": 9500, + "initial_energy_wh": initial_energy_wh, + "charge_efficiency": 0.95, + "discharge_efficiency": 0.95, + "cycle_cost_ore_kwh": 0, + "terminal_price_ore_kwh": 0, + } + ) + return request + + +@pytest.mark.parametrize("backend", ["auto", "cvxpy"]) +@pytest.mark.parametrize( + "delta", + [0.1e-6, 0.25e-6, 0.5e-6, 0.99e-6, 1e-6], +) +def test_multistage_normalizes_tiny_initial_over_maximum( + backend: str, delta: float +) -> None: + request = _multistage_soc_boundary_request(9500 + delta, backend) + clear_multistage_cache() + + response = handle(request) + + assert response["ok"], response + assert response["solver"]["formulation"] == "multistage-milp" + assert_storage_replays(request, response) + + +@pytest.mark.parametrize("backend", ["auto", "cvxpy"]) +def test_multistage_keeps_discrete_guard_for_material_initial_over_maximum( + backend: str, +) -> None: + request = _multistage_soc_boundary_request(9800, backend) + clear_multistage_cache() + + response = handle(request) + + assert response["ok"], response + assert response["solver"]["formulation"] == "multistage-milp" + assert_storage_replays(request, response) + + def test_multistage_auto_retries_with_storage_guard_after_direct_cycle(monkeypatch) -> None: from ftw_optimizer import direct_highs From d4ca6084f8aa140bc0b1913e6b46cb827abad0ca Mon Sep 17 00:00:00 2001 From: Fredrik Ahlgren Date: Tue, 4 Aug 2026 11:14:35 +0200 Subject: [PATCH 4/4] fix(optimizer): canonicalize storage state before builders --- .changeset/fix-optimizer-soc-above-max.md | 2 +- optimizer/ftw_optimizer/model.py | 20 ++++++++- optimizer/ftw_optimizer/multistage.py | 5 +++ optimizer/ftw_optimizer/recourse.py | 2 + optimizer/tests/test_model.py | 53 ++++++++++++++++++++--- 5 files changed, 73 insertions(+), 9 deletions(-) diff --git a/.changeset/fix-optimizer-soc-above-max.md b/.changeset/fix-optimizer-soc-above-max.md index 20e93c48..3d46ad13 100644 --- a/.changeset/fix-optimizer-soc-above-max.md +++ b/.changeset/fix-optimizer-soc-above-max.md @@ -2,4 +2,4 @@ "ftw": patch --- -Keep optimizer storage plans replay-consistent when a battery starts above its configured maximum state of charge, including solver-scale numeric noise at that boundary. +Canonicalize solver-scale state-of-charge noise to the exact operating maximum before building optimizer models, while keeping storage plans replay-consistent when a battery starts above its configured maximum. diff --git a/optimizer/ftw_optimizer/model.py b/optimizer/ftw_optimizer/model.py index 493af93e..4a4c33cf 100644 --- a/optimizer/ftw_optimizer/model.py +++ b/optimizer/ftw_optimizer/model.py @@ -48,11 +48,13 @@ class ReplayConsistencyError(RuntimeError): _REPLAY_TOLERANCE_FRACTION = 0.0002 _REPLAY_TOLERANCE_MIN_WH = 1.0 _STORAGE_NUMERIC_TOLERANCE_WH = 1e-6 +_STORAGE_INITIAL_ABOVE_MAXIMUM_KEY = "_optimizer_initial_above_maximum" def _storage_starts_above_maximum(storages: Any) -> bool: return any( - float(spec["initial_energy_wh"]) + bool(spec.get(_STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, False)) + or float(spec["initial_energy_wh"]) > float(spec.get("max_energy_wh", spec["capacity_wh"])) for spec in storages ) @@ -90,7 +92,9 @@ def _normalize_storage_specs( model_maximum, capacity ): model_maximum = capacity - initial_above_maximum = initial > model_maximum + initial_above_maximum = bool( + spec.get(_STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, False) + ) or initial > model_maximum model_initial = initial if initial_above_maximum and _within_storage_numeric_tolerance( initial, model_maximum @@ -102,11 +106,22 @@ def _normalize_storage_specs( model_spec["max_energy_wh"] = model_maximum if model_initial != initial: model_spec["initial_energy_wh"] = model_initial + model_spec[_STORAGE_INITIAL_ABOVE_MAXIMUM_KEY] = initial_above_maximum normalized.append(model_spec) starts_above_maximum.append(initial_above_maximum) return tuple(normalized), tuple(starts_above_maximum) +def _canonicalize_storage_payload(payload: dict[str, Any]) -> dict[str, Any]: + """Normalize storage state once before any policy or scenario builder runs.""" + storage_specs, _ = _normalize_storage_specs( + require_list(payload.get("storages", []), "storages") + ) + canonical = dict(payload) + canonical["storages"] = [dict(spec) for spec in storage_specs] + return canonical + + def _storage_replay_tolerance_wh(spec: dict[str, Any]) -> float: return max( _REPLAY_TOLERANCE_MIN_WH, @@ -223,6 +238,7 @@ def _solver_options(settings: dict[str, Any], solver: str) -> dict[str, Any]: def solve(payload: dict[str, Any]) -> dict[str, Any]: + payload = _canonicalize_storage_payload(payload) settings = require_dict(payload.get("settings", {}), "settings") commercial = require_dict( payload.get("commercial_constraints", {}), diff --git a/optimizer/ftw_optimizer/multistage.py b/optimizer/ftw_optimizer/multistage.py index 1d667b32..7f6bf97f 100644 --- a/optimizer/ftw_optimizer/multistage.py +++ b/optimizer/ftw_optimizer/multistage.py @@ -14,6 +14,8 @@ from .model import ( OPTIMAL_STATUSES, ReplayConsistencyError, + _STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, + _canonicalize_storage_payload, _export_price, _mode, _solver_options, @@ -343,6 +345,7 @@ def clear_multistage_cache() -> None: def _prepare(payload: dict[str, Any]) -> PreparedMultistage: + payload = _canonicalize_storage_payload(payload) settings = require_dict(payload.get("settings", {}), "settings") if require_list(payload.get("flex_loads", []), "flex_loads"): raise ProtocolError("multistage shadow does not yet support flex_loads") @@ -603,6 +606,8 @@ def _cache_key(prepared: PreparedMultistage) -> tuple[Any, ...]: float(spec["capacity_wh"]), float(spec.get("min_energy_wh", 0)), float(spec.get("max_energy_wh", spec["capacity_wh"])), + float(spec["initial_energy_wh"]), + bool(spec.get(_STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, False)), float(spec.get("max_charge_w", 0)), float(spec.get("max_discharge_w", 0)), float(spec.get("charge_efficiency", 0.95)), diff --git a/optimizer/ftw_optimizer/recourse.py b/optimizer/ftw_optimizer/recourse.py index df33d9c7..b878e522 100644 --- a/optimizer/ftw_optimizer/recourse.py +++ b/optimizer/ftw_optimizer/recourse.py @@ -11,6 +11,7 @@ from . import SCHEMA_VERSION from .model import ( OPTIMAL_STATUSES, + _canonicalize_storage_payload, _export_price, _mode, _solver_options, @@ -38,6 +39,7 @@ def solve_storage_recourse(payload: dict[str, Any]) -> dict[str, Any]: its shared first-stage action is intended for execution before replanning. """ + payload = _canonicalize_storage_payload(payload) started = time.perf_counter() settings = require_dict(payload.get("settings", {}), "settings") if require_list(payload.get("flex_loads", []), "flex_loads"): diff --git a/optimizer/tests/test_model.py b/optimizer/tests/test_model.py index cca2bb40..d95d9131 100644 --- a/optimizer/tests/test_model.py +++ b/optimizer/tests/test_model.py @@ -2,11 +2,13 @@ import copy import math +import threading import numpy as np import pytest from ftw_optimizer.multistage import clear_multistage_cache +from ftw_optimizer.model import _canonicalize_storage_payload from ftw_optimizer.scenario_tree import ( Scenario, build_scenario_tree, @@ -853,22 +855,61 @@ def _multistage_soc_boundary_request(initial_energy_wh: float, backend: str) -> return request +_MULTISTAGE_SOC_BOUNDARY_DELTAS = (0.1e-6, 0.25e-6, 0.5e-6, 0.99e-6, 1e-6) + + +def _assert_multistage_soc_boundary_case(backend: str, delta: float) -> None: + request = _multistage_soc_boundary_request(9500 + delta, backend) + canonical = _canonicalize_storage_payload(request) + storage = canonical["storages"][0] + assert storage["initial_energy_wh"] == storage["max_energy_wh"] == 9500.0 + + response = handle(request) + + assert response["ok"], response + assert response["solver"]["formulation"] == "multistage-milp" + assert_storage_replays(request, response) + + @pytest.mark.parametrize("backend", ["auto", "cvxpy"]) @pytest.mark.parametrize( "delta", - [0.1e-6, 0.25e-6, 0.5e-6, 0.99e-6, 1e-6], + _MULTISTAGE_SOC_BOUNDARY_DELTAS, ) def test_multistage_normalizes_tiny_initial_over_maximum( backend: str, delta: float ) -> None: - request = _multistage_soc_boundary_request(9500 + delta, backend) clear_multistage_cache() + _assert_multistage_soc_boundary_case(backend, delta) - response = handle(request) - assert response["ok"], response - assert response["solver"]["formulation"] == "multistage-milp" - assert_storage_replays(request, response) +@pytest.mark.parametrize("backend", ["auto", "cvxpy"]) +def test_multistage_soc_boundary_grid_is_deterministic(backend: str) -> None: + clear_multistage_cache() + for _ in range(4): + for delta in _MULTISTAGE_SOC_BOUNDARY_DELTAS: + _assert_multistage_soc_boundary_case(backend, delta) + + +def _burn_cpu(stop: threading.Event) -> None: + while not stop.is_set(): + sum(value * value for value in range(20_000)) + + +def test_multistage_soc_boundary_grid_survives_cpu_load() -> None: + clear_multistage_cache() + stop = threading.Event() + worker = threading.Thread(target=_burn_cpu, args=(stop,)) + worker.start() + try: + for _ in range(2): + for backend in ("auto", "cvxpy"): + for delta in _MULTISTAGE_SOC_BOUNDARY_DELTAS: + _assert_multistage_soc_boundary_case(backend, delta) + finally: + stop.set() + worker.join(timeout=5) + assert not worker.is_alive() @pytest.mark.parametrize("backend", ["auto", "cvxpy"])