diff --git a/.changeset/shared-storage-direct-highs.md b/.changeset/shared-storage-direct-highs.md new file mode 100644 index 00000000..9de3f30d --- /dev/null +++ b/.changeset/shared-storage-direct-highs.md @@ -0,0 +1,13 @@ +--- +"ftw": patch +--- + +The optional optimizer now uses highspy instead of CVXPY by default for +eligible shared storage plans, sending a sparse linear model straight to +HiGHS. It keeps the same shared-action constraints, service-first solve, +scenario risk cost and replay checks as the CVXPY model. It checks meter flow +against the post-curtailment baseline and retries an exact mixed-integer HiGHS +model within the same time budget when a relaxed candidate crosses import and +export. Auto mode uses this path only for continuous, cycle-safe HIGHS requests +with storage alone; commercial constraints, flexible loads, thermal loads, +guarded tariffs and direct-solver failures stay on CVXPY. diff --git a/go/internal/mpc/external_optimizer_test.go b/go/internal/mpc/external_optimizer_test.go index 876958e0..a0d01093 100644 --- a/go/internal/mpc/external_optimizer_test.go +++ b/go/internal/mpc/external_optimizer_test.go @@ -298,7 +298,8 @@ func TestExternalOptimizerEndToEnd(t *testing.T) { if err != nil { t.Fatalf("Optimize: %v", err) } - if plan.Solver == nil || plan.Solver.Engine != "cvxpy" || plan.Solver.Backend != "highs" { + if plan.Solver == nil || plan.Solver.Engine != "highspy" || plan.Solver.Backend != "highs" || + plan.Solver.ScenarioPolicy != "shared" || plan.Solver.PolicyVersion != "shared-v1" { t.Fatalf("unexpected solver metadata: %+v", plan.Solver) } if plan.Actions[0].BatteryW <= 0 || plan.Actions[1].BatteryW >= 0 { diff --git a/optimizer/ftw_optimizer/direct_highs.py b/optimizer/ftw_optimizer/direct_highs.py index b96a8fb5..fbdacab4 100644 --- a/optimizer/ftw_optimizer/direct_highs.py +++ b/optimizer/ftw_optimizer/direct_highs.py @@ -20,6 +20,21 @@ class DirectHighsError(RuntimeError): pass +SIMULTANEOUS_STORAGE_CYCLE_ERROR = ( + "HiGHS returned simultaneous storage charge and discharge" +) +SHARED_BASELINE_REPLAY_ERROR = ( + "direct shared mode violates the post-curtailment baseline" +) + + +class SharedBaselineReplayError(DirectHighsError): + def __init__(self, build_ms: float, solver_ms: float) -> None: + super().__init__(SHARED_BASELINE_REPLAY_ERROR) + self.build_ms = build_ms + self.solver_ms = solver_ms + + @dataclass class DirectScenarioVars: charge: list[list[int]] @@ -30,16 +45,36 @@ class DirectScenarioVars: grid_export: list[int] +@dataclass +class DirectSharedStorageVars: + charge: list[list[int]] + discharge: list[list[int]] + energy: list[list[int]] + total_charge: list[list[int]] + total_discharge: list[list[int]] + service: dict[int, float] + economic: dict[int, float] + + class SparseModel: def __init__(self) -> None: self.lower: list[float] = [] self.upper: list[float] = [] self.rows: list[tuple[dict[int, float], float, float]] = [] + self.integer: list[int] = [] - def variable(self, lower: float = 0.0, upper: float = highspy.kHighsInf) -> int: + def variable( + self, + lower: float = 0.0, + upper: float = highspy.kHighsInf, + *, + integer: bool = False, + ) -> int: index = len(self.lower) self.lower.append(lower) self.upper.append(upper) + if integer: + self.integer.append(index) return index def row( @@ -52,17 +87,38 @@ def row( self.rows.append((coefficients, lower, upper)) return index - def build(self, costs: np.ndarray, settings: dict[str, Any]) -> highspy.Highs: + def build( + self, + costs: np.ndarray, + settings: dict[str, Any], + *, + time_limit_s: float, + ) -> highspy.Highs: highs = highspy.Highs() highs.setOptionValue("output_flag", False) options = _solver_options(settings, "HIGHS") - highs.setOptionValue("time_limit", float(options["time_limit"])) + highs.setOptionValue( + "time_limit", min(float(options["time_limit"]), time_limit_s) + ) highs.setOptionValue("mip_rel_gap", float(options["mip_rel_gap"])) lower = np.asarray(self.lower, dtype=np.float64) upper = np.asarray(self.upper, dtype=np.float64) _require_ok(highs.addVars(len(lower), lower, upper), "add variables") indices = np.arange(len(lower), dtype=np.int32) _require_ok(highs.changeColsCost(len(lower), indices, costs), "set objective") + if self.integer: + integer_indices = np.asarray(self.integer, dtype=np.int32) + integer_types = np.full( + len(integer_indices), + highspy.HighsVarType.kInteger, + dtype=np.uint8, + ) + _require_ok( + highs.changeColsIntegrality( + len(integer_indices), integer_indices, integer_types + ), + "set integer variables", + ) starts = np.zeros(len(self.rows) + 1, dtype=np.int32) row_indices: list[int] = [] @@ -97,6 +153,12 @@ def solve_direct_highs( started: float, prepare_ms: float, decomposition: str, + *, + shared: bool = False, + exact_shared_baseline: bool = False, + deadline: float | None = None, + prior_build_ms: float = 0.0, + prior_solver_ms: float = 0.0, ) -> dict[str, Any]: if _storage_starts_above_maximum(prepared.storages): raise DirectHighsError( @@ -104,6 +166,10 @@ def solve_direct_highs( ) if prepared.discrete or prepared.unsafe_cycle or prepared.unsafe_meter_split: raise DirectHighsError("direct HiGHS path requires a cycle-safe continuous tariff") + if deadline is None: + deadline = started + float( + _solver_options(prepared.settings, "HIGHS")["time_limit"] + ) build_started = time.perf_counter() model = SparseModel() m = len(prepared.scenario_set.scenarios) @@ -113,6 +179,7 @@ def solve_direct_highs( ) service_terms: list[dict[int, float]] = [] economic_terms: list[dict[int, float]] = [] + risk_terms: list[dict[int, float]] = [] scenario_vars: list[DirectScenarioVars] = [] block_start_at = np.zeros(n, dtype=np.int64) @@ -120,6 +187,7 @@ def solve_direct_highs( block_start_at[block_start:block_end] = block_start storage_actions: dict[tuple[int, int, int], tuple[int, int]] = {} + shared_storage: DirectSharedStorageVars | None = None curtail_upper: dict[tuple[int, int], float] = {} for si, scenario in enumerate(prepared.scenario_set.scenarios): pv_generation = np.maximum(0.0, -scenario.pv) @@ -141,15 +209,25 @@ def solve_direct_highs( pv_generation = np.maximum(0.0, -scenario.pv) pv_surplus = np.maximum(0.0, pv_generation - scenario.load) base_import = np.maximum(0.0, scenario.load - pv_generation) - charges: list[list[int]] = [] - discharges: list[list[int]] = [] - energies: list[list[int]] = [] - total_charge: list[list[int]] = [[] for _ in range(n)] - total_discharge: list[list[int]] = [[] for _ in range(n)] - service: dict[int, float] = {} - economic: dict[int, float] = {} - - for storage_index, spec in enumerate(prepared.storages): + charges = shared_storage.charge if shared_storage is not None else [] + discharges = shared_storage.discharge if shared_storage is not None else [] + energies = shared_storage.energy if shared_storage is not None else [] + total_charge = ( + shared_storage.total_charge + if shared_storage is not None + else [[] for _ in range(n)] + ) + total_discharge = ( + shared_storage.total_discharge + if shared_storage is not None + else [[] for _ in range(n)] + ) + service = dict(shared_storage.service) if shared_storage is not None else {} + economic = dict(shared_storage.economic) if shared_storage is not None else {} + # Shared charge and discharge make storage state deterministic across + # scenarios. Reuse its variables and rows; only meter flow varies. + storages_to_build = () if shared_storage is not None else prepared.storages + for storage_index, spec in enumerate(storages_to_build): capacity = float(spec["capacity_wh"]) minimum = float(spec.get("min_energy_wh", 0)) maximum = float(spec.get("max_energy_wh", capacity)) @@ -217,19 +295,35 @@ def solve_direct_highs( _add(service, upper_recovery[t], 1.0 / (capacity * n)) if spec.get("target_energy_wh") is not None: - deadline = int(spec.get("target_slot", n - 1)) + target_slot = int(spec.get("target_slot", n - 1)) shortfall = model.variable() target = float(spec["target_energy_wh"]) - model.row({energy[deadline + 1]: 1.0, shortfall: 1.0}, target) + model.row( + {energy[target_slot + 1]: 1.0, shortfall: 1.0}, + target, + ) _add(service, shortfall, 1.0 / capacity) cycle_coefficient = spread + max(0.0, float(spec.get("cycle_cost_ore_kwh", 0))) + throughput_coefficient = ( + max(0.0, float(spec.get("throughput_cost_ore_kwh", 0))) + if shared + else 0.0 + ) for t in range(n): _add( economic, discharge[t], cycle_coefficient * float(prepared.dt_h[t]) / 1000.0, ) + if throughput_coefficient > 0: + coefficient = ( + throughput_coefficient + * float(prepared.dt_h[t]) + / 1000.0 + ) + _add(economic, charge[t], coefficient) + _add(economic, discharge[t], coefficient) _add( economic, energy[-1], @@ -238,6 +332,17 @@ def solve_direct_highs( charges.append(charge) discharges.append(discharge) energies.append(energy) + if shared and shared_storage is None: + shared_storage = DirectSharedStorageVars( + charges, + discharges, + energies, + total_charge, + total_discharge, + dict(service), + dict(economic), + ) + scenario_grid_cost: dict[int, float] = {} curtail = [ curtail_actions[(int(prepared.tree.node_at[si, t]), t)] for t in range(n) @@ -245,10 +350,70 @@ def solve_direct_highs( grid_import: list[int] = [] grid_export: list[int] = [] for t in range(n): + net = float(scenario.load[t] - pv_generation[t]) + shared_baseline_mode = shared and prepared.mode in { + "self_consumption", + "cheap_charge", + "passive_arbitrage", + } + baseline_stays_import = False + baseline_crosses = False + baseline_import_index: int | None = None + baseline_export_index: int | None = None + if shared_baseline_mode: + curtail_ceiling = float( + curtail_upper[ + (int(prepared.tree.node_at[si, t]), t) + ] + ) + baseline_crosses = net < 0 < net + curtail_ceiling + if baseline_crosses and exact_shared_baseline: + baseline_import_max = net + curtail_ceiling + baseline_export_max = -net + baseline_import_index = model.variable( + 0.0, baseline_import_max + ) + baseline_export_index = model.variable( + 0.0, baseline_export_max + ) + baseline_import_mode = model.variable( + 0.0, 1.0, integer=True + ) + model.row( + { + baseline_import_index: 1.0, + baseline_export_index: -1.0, + curtail[t]: -1.0, + }, + net, + net, + ) + model.row( + { + baseline_import_index: 1.0, + baseline_import_mode: -baseline_import_max, + }, + upper=0.0, + ) + model.row( + { + baseline_export_index: 1.0, + baseline_import_mode: baseline_export_max, + }, + upper=baseline_export_max, + ) + baseline_stays_import = net >= 0 + import_upper = float(prepared.import_bound[t]) export_upper = float(prepared.export_bound[t]) if prepared.mode == "self_consumption": - import_upper = min(import_upper, float(base_import[t]) + 50.0) + if shared_baseline_mode: + if not baseline_crosses and not baseline_stays_import: + import_upper = min(import_upper, 50.0) + else: + import_upper = min( + import_upper, float(base_import[t]) + 50.0 + ) import_index = model.variable(0.0, import_upper) export_index = model.variable(0.0, export_upper) grid_import.append(import_index) @@ -262,9 +427,59 @@ def solve_direct_highs( _add(balance, index, -1.0) for index in total_discharge[t]: _add(balance, index, 1.0) - net = float(scenario.load[t] - pv_generation[t]) model.row(balance, net, net) - if prepared.mode == "self_consumption": + if shared_baseline_mode and prepared.mode == "self_consumption": + if baseline_crosses and exact_shared_baseline: + assert baseline_import_index is not None + assert baseline_export_index is not None + model.row( + {import_index: 1.0, baseline_import_index: -1.0}, + upper=50.0, + ) + model.row( + {export_index: 1.0, baseline_export_index: -1.0}, + upper=50.0, + ) + elif baseline_crosses: + model.row( + {import_index: 1.0, curtail[t]: -1.0}, + upper=50.0, + ) + model.row( + {export_index: 1.0}, + upper=-net + 50.0, + ) + elif baseline_stays_import: + model.row( + {import_index: 1.0, curtail[t]: -1.0}, + upper=net + 50.0, + ) + model.row({export_index: 1.0}, upper=50.0) + else: + model.row( + {export_index: 1.0, curtail[t]: 1.0}, + upper=-net + 50.0, + ) + elif shared_baseline_mode: + if baseline_crosses and exact_shared_baseline: + assert baseline_export_index is not None + model.row( + {export_index: 1.0, baseline_export_index: -1.0}, + upper=1e-6, + ) + elif baseline_crosses: + model.row( + {export_index: 1.0}, + upper=-net + 1e-6, + ) + elif baseline_stays_import: + model.row({export_index: 1.0}, upper=1e-6) + else: + model.row( + {export_index: 1.0, curtail[t]: 1.0}, + upper=-net + 1e-6, + ) + elif prepared.mode == "self_consumption": model.row( {export_index: 1.0, curtail[t]: 1.0}, upper=float(pv_surplus[t]) + 50.0, @@ -279,11 +494,31 @@ def solve_direct_highs( import_index, float(prepared.effective_import[t] * prepared.dt_h[t] / 1000.0), ) + if shared: + _add( + scenario_grid_cost, + import_index, + float( + prepared.effective_import[t] + * prepared.dt_h[t] + / 1000.0 + ), + ) _add( economic, export_index, -float(prepared.effective_export[t] * prepared.dt_h[t] / 1000.0), ) + if shared: + _add( + scenario_grid_cost, + export_index, + -float( + prepared.effective_export[t] + * prepared.dt_h[t] + / 1000.0 + ), + ) if prepared.mode in {"self_consumption", "passive_arbitrage"}: house_import = model.variable() @@ -306,6 +541,7 @@ def solve_direct_highs( service_terms.append(service) economic_terms.append(economic) + risk_terms.append(scenario_grid_cost if shared else economic) scenario_vars.append( DirectScenarioVars( charges, discharges, energies, curtail, grid_import, grid_export @@ -347,9 +583,9 @@ def solve_direct_highs( service_costs = np.pad(service_costs, (0, 1 + m)) economic_costs = np.pad(economic_costs, (0, 1 + m)) economic_costs[threshold] += prepared.economic_cvar_weight - for si, economic in enumerate(economic_terms): + for si, risk in enumerate(risk_terms): row = {excess[si]: 1.0, threshold: 1.0} - for index, value in economic.items(): + for index, value in risk.items(): _add(row, index, -value) model.row(row, 0.0) economic_costs[excess[si]] += ( @@ -362,7 +598,11 @@ def solve_direct_highs( if len(economic_costs) < len(model.lower): economic_costs = np.pad(economic_costs, (0, len(model.lower) - len(economic_costs))) - highs = model.build(service_costs, prepared.settings) + highs = model.build( + service_costs, + prepared.settings, + time_limit_s=_remaining_time_s(deadline), + ) build_ms = (time.perf_counter() - build_started) * 1000.0 solver_started = time.perf_counter() _run_optimal(highs, "service") @@ -381,11 +621,30 @@ def solve_direct_highs( highs.changeColsCost(len(model.lower), column_indices, economic_costs), "set economic objective", ) + _require_ok( + highs.setOptionValue("time_limit", _remaining_time_s(deadline)), + "set economic time limit", + ) _run_optimal(highs, "economic") solver_ms = (time.perf_counter() - solver_started) * 1000.0 + mip_gap = float(highs.getInfo().mip_gap) if model.integer else None solution = np.asarray(highs.getSolution().col_value, dtype=np.float64) if len(solution) != len(model.lower) or not np.all(np.isfinite(solution)): raise DirectHighsError("HiGHS returned a non-finite solution") + if shared: + try: + _validate_shared_baseline_solution( + prepared, scenario_vars, solution + ) + except DirectHighsError as exc: + if ( + not exact_shared_baseline + and str(exc) == SHARED_BASELINE_REPLAY_ERROR + ): + raise SharedBaselineReplayError( + build_ms, solver_ms + ) from exc + raise return _response( prepared, @@ -395,11 +654,14 @@ def solve_direct_highs( best_service, started, prepare_ms, - build_ms, - solver_ms, + prior_build_ms + build_ms, + prior_solver_ms + solver_ms, decomposition, len(model.lower), len(model.rows), + len(model.integer), + mip_gap, + shared, ) @@ -416,9 +678,10 @@ def _response( decomposition: str, variables: int, constraints: int, + integer_variables: int, + mip_gap: float | None, + shared: bool, ) -> dict[str, Any]: - from .multistage import policy_config - scenarios = prepared.scenario_set.scenarios base_index = next((i for i, scenario in enumerate(scenarios) if scenario.id == "base"), 0) base = scenarios[base_index] @@ -427,9 +690,7 @@ def _response( 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" - ) + raise DirectHighsError(SIMULTANEOUS_STORAGE_CYCLE_ERROR) 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]] = [] @@ -457,6 +718,7 @@ def _response( ) raw_total_cost += raw_cost curtailed_w = max(0.0, float(solution[base_vars.curtail[t]])) + pv_forecast = prepared.base_pv if shared else base.pv actions.append( { "slot_start_ms": int(slot.get("start_ms", 0)), @@ -465,7 +727,7 @@ def _response( "grid_w": grid_w, "soc_pct": stored_wh / total_capacity * 100.0, "cost_ore": raw_cost, - "pv_limit_w": max(0.0, -base.pv[t] - curtailed_w) + "pv_limit_w": max(0.0, -pv_forecast[t] - curtailed_w) if curtailed_w > 1e-5 else 0.0, "storage_power_w": storage_power, @@ -477,42 +739,117 @@ def _response( } ) solve_ms = (time.perf_counter() - started) * 1000.0 + solver: dict[str, Any] = { + "engine": "highspy", + "backend": "highs", + "status": "optimal", + "formulation": ( + "milp" + if shared and integer_variables + else "convex" + if shared + else "multistage-lp" + ), + "objective_ore": objective, + "service_slack": best_service, + "solve_ms": solve_ms, + "prepare_ms": prepare_ms, + "build_ms": build_ms, + "solver_ms": solver_ms, + "cache_hit": False, + "dpp": False, + "mip_gap": mip_gap, + "scenario_count": len(scenarios), + "scenario_policy": "shared" if shared else "multistage", + "policy_version": "shared-v1" if shared else "storage-multistage-v1", + "non_anticipative_slots": prepared.first_stage_slots, + "model_variables": variables, + "model_constraints": constraints, + } + if shared: + probabilities = np.asarray([scenario.probability for scenario in scenarios]) + energy_cost = 0.0 + for probability, scenario in zip(probabilities, scenario_vars): + for t in range(prepared.n): + energy_cost += float(probability) * ( + prepared.effective_import[t] + * prepared.dt_h[t] + * solution[scenario.grid_import[t]] + / 1000.0 + - prepared.effective_export[t] + * prepared.dt_h[t] + * solution[scenario.grid_export[t]] + / 1000.0 + ) + degradation = 0.0 + terminal_value = 0.0 + spread = max( + 0.0, + finite_number( + prepared.settings.get("min_arbitrage_spread_ore_kwh", 0), + "settings.min_arbitrage_spread_ore_kwh", + ), + ) + for index, spec in enumerate(prepared.storages): + discharge_rate = spread + max( + 0.0, float(spec.get("cycle_cost_ore_kwh", 0)) + ) + throughput_rate = max( + 0.0, float(spec.get("throughput_cost_ore_kwh", 0)) + ) + for t in range(prepared.n): + degradation += ( + prepared.dt_h[t] + * ( + discharge_rate * solution[base_vars.discharge[index][t]] + + throughput_rate + * ( + solution[base_vars.charge[index][t]] + + solution[base_vars.discharge[index][t]] + ) + ) + / 1000.0 + ) + terminal_value -= ( + float(spec.get("terminal_price_ore_kwh", 0)) + * solution[base_vars.energy[index][-1]] + / 1000.0 + ) + solver.update( + { + "cvar_weight": prepared.economic_cvar_weight, + "cvar_alpha": prepared.economic_cvar_alpha, + "objective_breakdown_ore": { + "energy": float(energy_cost), + "demand_charge_increment": 0.0, + "degradation": float(degradation), + "terminal_energy_value": float(terminal_value), + }, + } + ) + else: + from .multistage import policy_config + + solver.update( + { + "scenario_original_count": prepared.scenario_set.original_count, + "scenario_reduction_error": prepared.scenario_set.reduction_error, + "policy_config": policy_config(prepared), + "tree_nodes": prepared.tree.node_count, + "move_blocks": len(prepared.blocks), + "decomposition": f"direct-highs-{decomposition}", + "risk_model": "service-cvar-then-expected-cost", + "service_cvar_weight": prepared.service_cvar_weight, + "service_cvar_alpha": prepared.service_cvar_alpha, + "economic_cvar_weight": prepared.economic_cvar_weight, + "economic_cvar_alpha": prepared.economic_cvar_alpha, + } + ) return { "schema_version": SCHEMA_VERSION, "request_id": str(prepared.payload["request_id"]), "ok": True, - "solver": { - "engine": "highspy", - "backend": "highs", - "status": "optimal", - "formulation": "multistage-lp", - "objective_ore": objective, - "service_slack": best_service, - "solve_ms": solve_ms, - "prepare_ms": prepare_ms, - "build_ms": build_ms, - "solver_ms": solver_ms, - "cache_hit": False, - "dpp": False, - "mip_gap": None, - "scenario_count": len(scenarios), - "scenario_original_count": prepared.scenario_set.original_count, - "scenario_reduction_error": prepared.scenario_set.reduction_error, - "scenario_policy": "multistage", - "policy_version": "storage-multistage-v1", - "policy_config": policy_config(prepared), - "non_anticipative_slots": prepared.first_stage_slots, - "tree_nodes": prepared.tree.node_count, - "move_blocks": len(prepared.blocks), - "decomposition": f"direct-highs-{decomposition}", - "risk_model": "service-cvar-then-expected-cost", - "service_cvar_weight": prepared.service_cvar_weight, - "service_cvar_alpha": prepared.service_cvar_alpha, - "economic_cvar_weight": prepared.economic_cvar_weight, - "economic_cvar_alpha": prepared.economic_cvar_alpha, - "model_variables": variables, - "model_constraints": constraints, - }, + "solver": solver, "plan": { "mode": prepared.mode, "horizon_slots": prepared.n, @@ -528,6 +865,49 @@ def _add(coefficients: dict[int, float], index: int, value: float) -> None: coefficients[index] = coefficients.get(index, 0.0) + value +def _remaining_time_s(deadline: float) -> float: + remaining = deadline - time.perf_counter() + if remaining <= 0.0: + raise DirectHighsError("direct HiGHS time budget exhausted") + return remaining + + +def _validate_shared_baseline_solution( + prepared: "PreparedMultistage", + scenario_vars: list[DirectScenarioVars], + solution: np.ndarray, +) -> None: + if prepared.mode not in { + "self_consumption", + "cheap_charge", + "passive_arbitrage", + }: + return + tolerance = 1e-4 + for scenario, variables in zip( + prepared.scenario_set.scenarios, scenario_vars + ): + pv_generation = np.maximum(0.0, -scenario.pv) + for t in range(prepared.n): + curtailed = float(solution[variables.curtail[t]]) + baseline = float( + scenario.load[t] - pv_generation[t] + curtailed + ) + baseline_import = max(baseline, 0.0) + baseline_export = max(-baseline, 0.0) + grid_import = float(solution[variables.grid_import[t]]) + grid_export = float(solution[variables.grid_export[t]]) + if prepared.mode == "self_consumption": + valid = ( + grid_import <= baseline_import + 50.0 + tolerance + and grid_export <= baseline_export + 50.0 + tolerance + ) + else: + valid = grid_export <= baseline_export + 1e-6 + tolerance + if not valid: + raise DirectHighsError(SHARED_BASELINE_REPLAY_ERROR) + + def _accumulate(target: np.ndarray, terms: dict[int, float], weight: float) -> None: for index, value in terms.items(): target[index] += weight * value diff --git a/optimizer/ftw_optimizer/model.py b/optimizer/ftw_optimizer/model.py index 4a4c33cf..6cf57fcd 100644 --- a/optimizer/ftw_optimizer/model.py +++ b/optimizer/ftw_optimizer/model.py @@ -268,6 +268,42 @@ def solve(payload: dict[str, Any]) -> dict[str, Any]: raise ProtocolError("settings.scenario_policy must be shared, recourse, or multistage") started = time.perf_counter() + shared_backend = str(settings.get("shared_backend", "auto")) + if shared_backend not in {"auto", "highs", "cvxpy"}: + raise ProtocolError("settings.shared_backend must be auto, highs, or cvxpy") + direct_fallback_reason = "" + direct_storage_guard = False + if shared_backend in {"auto", "highs"}: + from .direct_highs import ( + DirectHighsError, + SIMULTANEOUS_STORAGE_CYCLE_ERROR, + ) + from .shared_highs import DirectSharedIneligible, solve_shared_highs + + try: + response = solve_shared_highs(payload, started) + _validate_storage_replay( + response["plan"]["actions"], + require_list(payload.get("slots", []), "slots"), + require_list(payload.get("storages", []), "storages"), + ) + return response + except DirectSharedIneligible as exc: + if shared_backend == "highs": + raise ProtocolError(str(exc)) from exc + except Exception as exc: + if shared_backend == "highs": + raise + # The direct path is optional in auto mode. Let the reference + # model validate the request again as it builds the fallback. + direct_fallback_reason = str(exc) or type(exc).__name__ + direct_storage_guard = isinstance( + exc, ReplayConsistencyError + ) or ( + isinstance(exc, DirectHighsError) + and str(exc) == SIMULTANEOUS_STORAGE_CYCLE_ERROR + ) + slots = [require_dict(v, f"slots[{i}]") for i, v in enumerate(require_list(payload["slots"], "slots"))] n = len(slots) mode = _mode(payload) @@ -457,7 +493,12 @@ def solve(payload: dict[str, Any]) -> dict[str, Any]: 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_above_maximum[i] - if force_milp or initial_above_max or (formulation == "auto" and unsafe_cycle): + if ( + force_milp + or direct_storage_guard + 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 @@ -960,7 +1001,7 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: break solve_ms = (time.perf_counter() - started) * 1000.0 _validate_storage_replay(actions, slots, [storage.spec for storage in storages]) - return { + response = { "schema_version": SCHEMA_VERSION, "request_id": str(payload["request_id"]), "ok": True, @@ -995,3 +1036,7 @@ def run_problem(problem: cp.Problem, solver_name: str) -> None: "actions": actions, }, } + if direct_fallback_reason: + response["solver"]["fallback"] = True + response["solver"]["fallback_reason"] = direct_fallback_reason + return response diff --git a/optimizer/ftw_optimizer/shared_highs.py b/optimizer/ftw_optimizer/shared_highs.py new file mode 100644 index 00000000..ea5c5fbc --- /dev/null +++ b/optimizer/ftw_optimizer/shared_highs.py @@ -0,0 +1,253 @@ +from __future__ import annotations + +import copy +import time +from dataclasses import replace +from typing import Any + +import numpy as np + +from .direct_highs import ( + SharedBaselineReplayError, + solve_direct_highs, +) +from .model import _STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, _solver_options +from .multistage import _prepare +from .protocol import ProtocolError, finite_number, require_dict, require_list +from .scenario_tree import ScenarioTree + + +class DirectSharedIneligible(RuntimeError): + pass + + +def solve_shared_highs(payload: dict[str, Any], started: float) -> dict[str, Any]: + """Solve shared storage through the sparse HiGHS builder.""" + prepared_started = time.perf_counter() + direct_payload, risk_alpha = _direct_payload(payload) + prepared = replace( + _prepare(direct_payload), economic_cvar_alpha=risk_alpha + ) + + solver = str(prepared.settings.get("solver", "HIGHS")).upper() + if solver not in {"HIGHS", "CLARABEL"}: + raise ProtocolError("settings.solver must be HIGHS or CLARABEL") + if solver != "HIGHS": + raise DirectSharedIneligible("direct shared backend requires solver HIGHS") + if prepared.formulation == "milp": + raise DirectSharedIneligible( + "direct shared backend requires a continuous formulation" + ) + if any( + bool(spec.get(_STORAGE_INITIAL_ABOVE_MAXIMUM_KEY, False)) + for spec in prepared.storages + ): + raise DirectSharedIneligible( + "direct shared backend requires storage starts at or below its operating maximum" + ) + if ( + prepared.discrete + or prepared.unsafe_cycle + or bool(np.any(prepared.effective_import < 0)) + or prepared.unsafe_meter_split + ): + raise DirectSharedIneligible( + "direct shared backend requires a cycle-safe continuous tariff" + ) + + shared_pv_generation = np.minimum.reduce( + [ + np.maximum(0.0, -scenario.pv) + for scenario in prepared.scenario_set.scenarios + ] + ) + scenario_count = len(prepared.scenario_set.scenarios) + shared_tree = ScenarioTree( + node_at=np.zeros((scenario_count, prepared.n), dtype=np.int64), + branch_slots=(), + node_count=1, + ) + blocks = tuple((slot, slot + 1) for slot in range(prepared.n)) + + # Match the shared champion's fallback site bound. Explicit slot limits + # still take precedence in both implementations. + max_site_power = max( + 1000.0, + float(np.max(prepared.base_load + shared_pv_generation)) + + sum( + float(spec.get("max_charge_w", 0)) + + float(spec.get("max_discharge_w", 0)) + for spec in prepared.storages + ), + ) + raw_import_limit = np.asarray( + [ + max( + 0.0, + finite_number( + slot.get("max_import_w", 0), + f"slots[{index}].max_import_w", + ), + ) + for index, slot in enumerate(prepared.slots) + ] + ) + raw_export_limit = np.asarray( + [ + max( + 0.0, + finite_number( + slot.get("max_export_w", 0), + f"slots[{index}].max_export_w", + ), + ) + for index, slot in enumerate(prepared.slots) + ] + ) + prepared = replace( + prepared, + tree=shared_tree, + blocks=blocks, + first_stage_slots=prepared.n, + service_cvar_weight=0.0, + max_site_power=max_site_power, + import_bound=np.where( + raw_import_limit > 0, raw_import_limit, max_site_power + ), + export_bound=np.where( + raw_export_limit > 0, raw_export_limit, max_site_power + ), + ) + prepare_ms = (time.perf_counter() - prepared_started) * 1000.0 + deadline = started + float( + _solver_options(prepared.settings, "HIGHS")["time_limit"] + ) + try: + return solve_direct_highs( + prepared, + started, + prepare_ms, + "shared", + shared=True, + deadline=deadline, + ) + except SharedBaselineReplayError as exc: + return solve_direct_highs( + prepared, + started, + prepare_ms, + "shared", + shared=True, + exact_shared_baseline=True, + deadline=deadline, + prior_build_ms=exc.build_ms, + prior_solver_ms=exc.solver_ms, + ) + + +def _direct_payload(payload: dict[str, Any]) -> tuple[dict[str, Any], float]: + if require_dict( + payload.get("commercial_constraints", {}), "commercial_constraints" + ): + raise DirectSharedIneligible( + "direct shared backend does not support commercial constraints" + ) + if require_list(payload.get("flex_loads", []), "flex_loads"): + raise DirectSharedIneligible( + "direct shared backend does not support flex loads" + ) + if require_list(payload.get("thermal_loads", []), "thermal_loads"): + raise DirectSharedIneligible( + "direct shared backend does not support thermal loads" + ) + storages = require_list(payload.get("storages", []), "storages") + if not storages: + raise DirectSharedIneligible( + "direct shared backend requires at least one storage" + ) + slots = require_list(payload.get("slots", []), "slots") + if not slots: + raise ProtocolError("slots must not be empty") + + direct_payload = copy.deepcopy(payload) + direct_storages = require_list(direct_payload["storages"], "storages") + for index, raw in enumerate(direct_storages): + spec = require_dict(raw, f"storages[{index}]") + finite_number( + spec.get("max_charge_w", 0), f"storages[{index}].max_charge_w" + ) + finite_number( + spec.get("max_discharge_w", 0), + f"storages[{index}].max_discharge_w", + ) + finite_number( + spec.get("cycle_cost_ore_kwh", 0), + f"storages[{index}].cycle_cost_ore_kwh", + ) + finite_number( + spec.get("throughput_cost_ore_kwh", 0), + f"storages[{index}].throughput_cost_ore_kwh", + ) + finite_number( + spec.get("terminal_price_ore_kwh", 0), + f"storages[{index}].terminal_price_ore_kwh", + ) + if spec.get("target_energy_wh") is not None: + finite_number( + spec["target_energy_wh"], + f"storages[{index}].target_energy_wh", + ) + spec["target_slot"] = min( + len(slots) - 1, + max(0, int(spec.get("target_slot", len(slots) - 1))), + ) + + settings = dict( + require_dict(direct_payload.get("settings", {}), "settings") + ) + direct_payload["settings"] = settings + raw_scenarios = require_list(direct_payload.get("scenarios", []), "scenarios") + scenario_count = max(1, len(raw_scenarios)) + seen_scenario_ids: set[str] = set() + for index, raw in enumerate(raw_scenarios): + scenario = require_dict(raw, f"scenarios[{index}]") + scenario_id = str(scenario.get("id", f"scenario-{index}")) + if scenario_id in seen_scenario_ids: + suffix = 1 + unique_id = f"{scenario_id}-{index}-{suffix}" + while unique_id in seen_scenario_ids: + suffix += 1 + unique_id = f"{scenario_id}-{index}-{suffix}" + scenario["id"] = unique_id + scenario_id = unique_id + seen_scenario_ids.add(scenario_id) + risk_weight = max( + 0.0, + finite_number(settings.get("cvar_weight", 0), "settings.cvar_weight"), + ) + risk_alpha = finite_number( + settings.get("cvar_alpha", 0.9), "settings.cvar_alpha" + ) + if risk_weight > 0 and scenario_count > 1 and not 0 < risk_alpha < 1: + raise ProtocolError("settings.cvar_alpha must be between 0 and 1") + + # _prepare also serves multistage models. Pin its policy-only settings, + # then replace the generated tree with the exact shared policy above. + settings.update( + { + "scenario_limit": scenario_count, + "non_anticipative_slots": len(slots), + "branch_interval_slots": 1, + "branch_horizon_slots": len(slots), + "max_branching": 2, + "near_horizon_slots": len(slots), + "mid_horizon_slots": len(slots), + "mid_block_slots": 1, + "far_block_slots": 1, + "service_cvar_weight": 0, + "service_cvar_alpha": 0.95, + "economic_cvar_weight": risk_weight, + "economic_cvar_alpha": 0.9, + } + ) + return direct_payload, risk_alpha diff --git a/optimizer/tests/test_model.py b/optimizer/tests/test_model.py index d95d9131..d337ba62 100644 --- a/optimizer/tests/test_model.py +++ b/optimizer/tests/test_model.py @@ -1,12 +1,15 @@ from __future__ import annotations import copy +import json import math import threading +import time import numpy as np import pytest +from ftw_optimizer.direct_highs import DirectHighsError, _remaining_time_s from ftw_optimizer.multistage import clear_multistage_cache from ftw_optimizer.model import _canonicalize_storage_payload from ftw_optimizer.scenario_tree import ( @@ -153,6 +156,16 @@ def base_request() -> dict: } +def test_direct_highs_accepts_a_positive_sub_50ms_budget() -> None: + remaining = _remaining_time_s(time.perf_counter() + 0.01) + assert 0.0 < remaining <= 0.01 + + +def test_direct_highs_rejects_an_exhausted_budget() -> None: + with pytest.raises(DirectHighsError, match="time budget exhausted"): + _remaining_time_s(time.perf_counter() - 0.001) + + def assert_storage_replays(request: dict, response: dict, tolerance_wh: float = 2.1) -> None: energies = { str(spec["id"]): float(spec["initial_energy_wh"]) @@ -175,6 +188,885 @@ def assert_storage_replays(request: dict, response: dict, tolerance_wh: float = energies[storage_id] = replayed +def assert_nested_close( + direct: object, + reference: object, + *, + abs_tol: float = 1e-3, + path: str = "value", +) -> None: + if isinstance(direct, dict) and isinstance(reference, dict): + assert direct.keys() == reference.keys(), path + for key in direct: + assert_nested_close( + direct[key], + reference[key], + abs_tol=abs_tol, + path=f"{path}.{key}", + ) + return + if isinstance(direct, list) and isinstance(reference, list): + assert len(direct) == len(reference), path + for index, (direct_item, reference_item) in enumerate( + zip(direct, reference) + ): + assert_nested_close( + direct_item, + reference_item, + abs_tol=abs_tol, + path=f"{path}[{index}]", + ) + return + if ( + isinstance(direct, (int, float)) + and not isinstance(direct, bool) + and isinstance(reference, (int, float)) + and not isinstance(reference, bool) + ): + assert math.isclose( + float(direct), float(reference), abs_tol=abs_tol + ), f"{path}: {direct} != {reference}" + return + assert direct == reference, path + + +def assert_shared_plan_parity( + direct: dict, + reference: dict, + *, + abs_tol: float = 1e-3, +) -> None: + direct_plan = direct["plan"] + reference_plan = reference["plan"] + for key in ( + "mode", + "horizon_slots", + "capacity_wh", + "initial_soc_pct", + "total_cost_ore", + ): + assert_nested_close( + direct_plan[key], + reference_plan[key], + abs_tol=abs_tol, + path=f"plan.{key}", + ) + assert_nested_close( + direct_plan["actions"], + reference_plan["actions"], + abs_tol=abs_tol, + path="plan.actions", + ) + + +def test_shared_direct_highs_matches_cvxpy_with_risk_targets_and_costs() -> None: + request = base_request() + request["request_id"] = "shared-direct-parity" + request["settings"].update( + { + "shared_backend": "highs", + "cvar_weight": 0.25, + "cvar_alpha": 0.75, + "min_arbitrage_spread_ore_kwh": 7, + } + ) + prices = [20, 35, 70, 260, 310, 120] + loads = [1200, 1800, 900, 2600, 3200, 1600] + pv = [0, -500, -2400, -800, 0, -300] + request["slots"] = [ + { + "start_ms": 1 + index * 900_000, + "len_min": 15, + "price_ore": price, + "spot_ore": 0, + "confidence": 0.85, + "pv_w": pv[index], + "load_w": loads[index], + "max_import_w": 10_000, + "max_export_w": 10_000, + } + for index, price in enumerate(prices) + ] + request["storages"][0].update( + { + "initial_energy_wh": 2500, + "target_energy_wh": 4500, + "target_slot": 2, + "throughput_cost_ore_kwh": 2, + } + ) + request["storages"].append( + { + "id": "garage", + "capacity_wh": 6000, + "initial_energy_wh": 3500, + "min_energy_wh": 800, + "max_energy_wh": 5600, + "max_charge_w": 2600, + "max_discharge_w": 2200, + "charge_efficiency": 0.92, + "discharge_efficiency": 0.9, + "terminal_price_ore_kwh": 14, + "cycle_cost_ore_kwh": 11, + "throughput_cost_ore_kwh": 1.5, + } + ) + request["scenarios"] = [ + { + "id": "base", + "probability": 0.5, + "load_w": loads, + "pv_w": pv, + }, + { + "id": "cloudy", + "probability": 0.3, + "load_w": [value * 1.15 for value in loads], + "pv_w": [value * 0.65 for value in pv], + }, + { + "id": "sunny", + "probability": 0.2, + "load_w": [value * 0.9 for value in loads], + "pv_w": [value * 1.2 for value in pv], + }, + ] + + reference_request = copy.deepcopy(request) + reference_request["request_id"] = "shared-cvxpy-parity" + reference_request["settings"]["shared_backend"] = "cvxpy" + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert direct["solver"]["engine"] == "highspy" + assert reference["solver"]["engine"] == "cvxpy" + assert direct["solver"]["scenario_policy"] == "shared" + assert direct["solver"]["formulation"] == "convex" + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-4, + ) + for key in ( + "energy", + "demand_charge_increment", + "degradation", + "terminal_energy_value", + ): + assert math.isclose( + direct["solver"]["objective_breakdown_ore"][key], + reference["solver"]["objective_breakdown_ore"][key], + abs_tol=1e-4, + ) + for direct_action, reference_action in zip( + direct["plan"]["actions"], reference["plan"]["actions"] + ): + assert math.isclose( + direct_action["battery_w"], reference_action["battery_w"], abs_tol=1e-3 + ) + assert math.isclose( + direct_action["grid_w"], reference_action["grid_w"], abs_tol=1e-3 + ) + assert_storage_replays(request, direct) + + +@pytest.mark.parametrize( + "mode", + ["arbitrage", "cheap_charge", "passive_arbitrage", "self_consumption"], +) +def test_shared_direct_highs_matches_cvxpy_modes(mode: str) -> None: + request = base_request() + request["settings"].update({"mode": mode, "shared_backend": "highs"}) + for slot in request["slots"]: + slot["spot_ore"] = 0 + reference_request = copy.deepcopy(request) + reference_request["request_id"] = f"shared-{mode}-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-3, + ) + assert_storage_replays(request, direct) + + +def test_shared_direct_highs_matches_cvxpy_below_minimum_recovery() -> None: + request = base_request() + request["request_id"] = "shared-below-minimum-direct" + request["settings"].update( + { + "mode": "arbitrage", + "formulation": "relaxed", + "shared_backend": "highs", + } + ) + request["slots"] = [dict(request["slots"][0]) for _ in range(4)] + for index, slot in enumerate(request["slots"]): + slot.update( + { + "start_ms": 1 + index * 900_000, + "len_min": 15, + "price_ore": 20 + index * 40, + "spot_ore": 0, + } + ) + request["storages"][0].update( + { + "initial_energy_wh": 500, + "max_charge_w": 1000, + "max_discharge_w": 1000, + "terminal_price_ore_kwh": 0, + } + ) + reference_request = copy.deepcopy(request) + reference_request["request_id"] = "shared-below-minimum-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert direct["solver"]["engine"] == "highspy" + assert reference["solver"]["engine"] == "cvxpy" + assert math.isclose( + direct["solver"]["service_slack"], + reference["solver"]["service_slack"], + abs_tol=1e-7, + ) + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-4, + ) + assert_shared_plan_parity(direct, reference) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) + + +def test_shared_direct_highs_matches_strict_pv_surplus_and_limit() -> None: + request = base_request() + request["request_id"] = "shared-pv-surplus-direct" + request["settings"].update( + { + "mode": "passive_arbitrage", + "formulation": "relaxed", + "shared_backend": "highs", + } + ) + request["slots"] = [ + { + "start_ms": 1, + "len_min": 60, + "price_ore": 100, + "spot_ore": 50, + "confidence": 1, + "pv_w": -4000, + "load_w": 500, + "max_import_w": 8000, + "max_export_w": 100, + } + ] + request["storages"][0].update( + { + "initial_energy_wh": 9500, + "max_charge_w": 0, + "max_discharge_w": 0, + "terminal_price_ore_kwh": 0, + } + ) + reference_request = copy.deepcopy(request) + reference_request["request_id"] = "shared-pv-surplus-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert_shared_plan_parity(direct, reference) + action = direct["plan"]["actions"][0] + assert math.isclose(action["battery_w"], 0, abs_tol=1e-6) + assert math.isclose(action["grid_w"], -100, abs_tol=1e-6) + assert math.isclose(action["pv_limit_w"], 600, abs_tol=1e-6) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) + + +def shared_curtailment_request(mode: str, base_load_w: float) -> dict: + request = base_request() + request["request_id"] = f"shared-curtail-{mode}-{base_load_w}" + request["settings"].update( + { + "mode": mode, + "formulation": "relaxed", + "shared_backend": "auto", + } + ) + request["slots"] = [ + { + "start_ms": 1, + "len_min": 60, + "price_ore": 100, + "spot_ore": 20, + "confidence": 1, + "pv_w": -1000, + "load_w": base_load_w, + "max_import_w": 8000, + "max_export_w": 100, + } + ] + request["scenarios"] = [ + { + "id": "base", + "probability": 0.5, + "pv_w": [-1000], + "load_w": [base_load_w], + }, + { + "id": "sunny", + "probability": 0.5, + "pv_w": [-1600], + "load_w": [500], + }, + ] + request["storages"][0].update( + { + "initial_energy_wh": 5000, + "max_charge_w": 0, + "max_discharge_w": 0, + "terminal_price_ore_kwh": 0, + } + ) + return request + + +@pytest.mark.parametrize( + "mode", + ["self_consumption", "cheap_charge", "passive_arbitrage"], +) +def test_shared_direct_highs_models_post_curtailment_baseline(mode: str) -> None: + request = shared_curtailment_request(mode, 2000) + request["request_id"] = f"shared-post-curtail-{mode}-direct" + request["settings"]["shared_backend"] = "highs" + reference_request = copy.deepcopy(request) + reference_request["request_id"] = f"shared-post-curtail-{mode}-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert direct["solver"]["engine"] == "highspy" + assert reference["solver"]["engine"] == "cvxpy" + assert_shared_plan_parity(direct, reference) + assert math.isclose(direct["plan"]["actions"][0]["grid_w"], 2000) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) + + +@pytest.mark.parametrize( + "mode", + ["self_consumption", "cheap_charge", "passive_arbitrage"], +) +def test_shared_curtailment_direction_change_matches_cvxpy(mode: str) -> None: + request = shared_curtailment_request(mode, 900) + request["request_id"] = f"shared-curtail-cross-{mode}-auto" + reference_request = copy.deepcopy(request) + reference_request["request_id"] = f"shared-curtail-cross-{mode}-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert direct["solver"]["engine"] == "highspy" + assert direct["solver"]["formulation"] == "convex" + assert direct["solver"]["mip_gap"] is None + assert reference["solver"]["engine"] == "cvxpy" + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-4, + ) + assert_shared_plan_parity(direct, reference) + assert math.isclose(direct["plan"]["actions"][0]["grid_w"], 900) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) + + +def test_shared_baseline_replay_retries_with_exact_highs() -> None: + request = shared_curtailment_request("self_consumption", 900) + request["request_id"] = "shared-curtail-exact-retry" + request["settings"]["shared_backend"] = "highs" + request["storages"][0].update( + { + "max_charge_w": 5000, + "terminal_price_ore_kwh": 1000, + } + ) + reference_request = copy.deepcopy(request) + reference_request["request_id"] = "shared-curtail-exact-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert direct["solver"]["engine"] == "highspy" + assert direct["solver"]["formulation"] == "milp" + assert direct["solver"]["mip_gap"] is not None + assert direct["solver"]["build_ms"] > 0 + assert direct["solver"]["solver_ms"] > 0 + assert reference["solver"]["engine"] == "cvxpy" + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-4, + ) + assert_shared_plan_parity(direct, reference) + assert math.isclose(direct["plan"]["actions"][0]["battery_w"], 50) + assert math.isclose(direct["plan"]["actions"][0]["grid_w"], 900) + + +def test_shared_auto_falls_back_at_each_direct_eligibility_boundary() -> None: + cases: list[tuple[str, dict]] = [] + + commercial = base_request() + commercial["commercial_constraints"] = {"version": "srcful-commercial-v1"} + cases.append(("commercial", commercial)) + + flex = base_request() + flex["flex_loads"] = [ + { + "id": "car", + "capacity_wh": 40_000, + "initial_energy_wh": 20_000, + "max_energy_wh": 40_000, + "target_energy_wh": 20_000, + "target_slot": 1, + "charge_efficiency": 0.9, + "allowed_steps_w": [0, 2000], + } + ] + cases.append(("flex", flex)) + + thermal = base_request() + thermal["thermal_loads"] = [ + { + "id": "heater", + "initial_temp_c": 20, + "min_temp_c": 18, + "max_temp_c": 24, + "outside_temp_c": [10, 10], + "allowed_steps_w": [0, 1000], + "gain_c_per_kwh": 1, + "loss_per_hour": 0.1, + } + ] + cases.append(("thermal", thermal)) + + no_storage = base_request() + no_storage["storages"] = [] + cases.append(("no-storage", no_storage)) + + clarabel = base_request() + clarabel["settings"].update({"solver": "CLARABEL", "formulation": "relaxed"}) + cases.append(("clarabel", clarabel)) + + milp = base_request() + milp["settings"]["formulation"] = "milp" + cases.append(("milp", milp)) + + negative_import = base_request() + negative_import["slots"][0]["price_ore"] = -10 + cases.append(("unsafe-cycle", negative_import)) + + pv_charge_bonus = base_request() + pv_charge_bonus["settings"]["pv_charge_bonus_ore_kwh"] = 1 + cases.append(("pv-charge-bonus", pv_charge_bonus)) + + meter_split = base_request() + meter_split["settings"]["export_ore_per_kwh"] = 400 + cases.append(("unsafe-meter-split", meter_split)) + + above_maximum = base_request() + above_maximum["storages"][0]["initial_energy_wh"] = 9800 + cases.append(("initial-above-maximum", above_maximum)) + + for name, request in cases: + request["request_id"] = f"shared-auto-boundary-{name}" + request["settings"]["shared_backend"] = "auto" + response = handle(request) + assert response["ok"], (name, response) + assert response["solver"]["engine"] == "cvxpy", (name, response) + assert response["solver"]["scenario_policy"] == "shared" + + +def test_shared_auto_preserves_duplicate_ids_and_first_base_output() -> None: + request = base_request() + request["settings"]["shared_backend"] = "auto" + base_load = [slot["load_w"] for slot in request["slots"]] + base_pv = [slot["pv_w"] for slot in request["slots"]] + request["scenarios"] = [ + {"id": "base", "probability": 0.5, "load_w": base_load, "pv_w": base_pv}, + { + "id": "base", + "probability": 0.5, + "load_w": [3500, 6000], + "pv_w": [-500, -1000], + }, + ] + reference_request = copy.deepcopy(request) + reference_request["request_id"] = "duplicate-base-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert direct["solver"]["engine"] == "highspy" + assert direct["solver"]["scenario_count"] == 2 + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-4, + ) + assert_shared_plan_parity(direct, reference) + for action, load_w, pv_w in zip( + direct["plan"]["actions"], base_load, base_pv + ): + assert math.isclose( + action["grid_w"], + load_w + pv_w + action["battery_w"], + abs_tol=1e-4, + ) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) + + +@pytest.mark.parametrize("target_slot", [-4, 99]) +def test_shared_direct_clamps_target_slot_like_cvxpy(target_slot: int) -> None: + request = base_request() + request["settings"]["shared_backend"] = "highs" + request["storages"][0].update( + {"target_energy_wh": 5000, "target_slot": target_slot} + ) + reference_request = copy.deepcopy(request) + reference_request["request_id"] = f"target-slot-{target_slot}-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + assert math.isclose( + direct["solver"]["objective_ore"], + reference["solver"]["objective_ore"], + abs_tol=1e-4, + ) + assert_storage_replays(request, direct) + + +def realistic_shared_request(scenario_count: int) -> dict: + request = base_request() + request["request_id"] = f"realistic-shared-{scenario_count}" + request["settings"].update( + { + "mode": "passive_arbitrage", + "solver": "HIGHS", + "formulation": "relaxed", + "time_limit_s": 8, + "shared_backend": "highs", + "cvar_weight": 0.15, + "cvar_alpha": 0.9, + } + ) + slots = [] + base_load = [] + base_pv = [] + for index in range(192): + hour = (index % 96) / 4.0 + price = 80 + 180 * math.exp(-0.5 * ((hour - 18) / 2) ** 2) + pv_w = ( + -7000 * math.exp(-0.5 * ((hour - 12.5) / 3) ** 2) + if 5 < hour < 21 + else 0 + ) + load_w = 500 + 1800 * math.exp(-0.5 * ((hour - 19) / 2) ** 2) + base_load.append(load_w) + base_pv.append(pv_w) + slots.append( + { + "start_ms": 1 + index * 900_000, + "len_min": 15, + "price_ore": price, + "spot_ore": price * 0.7, + "confidence": 1 if index < 96 else 0.6, + "pv_w": pv_w, + "load_w": load_w, + "max_import_w": 11_000, + "max_export_w": 11_000, + } + ) + request["slots"] = slots + request["storages"] = [ + { + "id": "home", + "capacity_wh": 15_000, + "initial_energy_wh": 7500, + "min_energy_wh": 1500, + "max_energy_wh": 14_250, + "max_charge_w": 5000, + "max_discharge_w": 5000, + "charge_efficiency": 0.95, + "discharge_efficiency": 0.95, + "terminal_price_ore_kwh": 150, + "cycle_cost_ore_kwh": 10, + "throughput_cost_ore_kwh": 1.5, + } + ] + scenarios = [] + for index in range(scenario_count): + offset = index - (scenario_count - 1) / 2 + scenarios.append( + { + "id": "base" if index == 0 else f"scenario-{index}", + "probability": 1 / scenario_count, + "load_w": [max(0, value + offset * 250) for value in base_load], + "pv_w": [min(0, value + offset * 150) for value in base_pv], + } + ) + scenarios[0]["load_w"] = base_load + scenarios[0]["pv_w"] = base_pv + request["scenarios"] = scenarios + return request + + +@pytest.mark.parametrize("scenario_count", [3, 12]) +def test_shared_direct_realistic_horizon_matches_cvxpy( + scenario_count: int, +) -> None: + request = realistic_shared_request(scenario_count) + reference_request = copy.deepcopy(request) + reference_request["request_id"] += "-cvxpy" + reference_request["settings"]["shared_backend"] = "cvxpy" + + direct = handle(request) + reference = handle(reference_request) + + assert direct["ok"], direct + assert reference["ok"], reference + direct_solver = direct["solver"] + reference_solver = reference["solver"] + assert direct_solver["engine"] == "highspy" + assert direct_solver["backend"] == "highs" + assert direct_solver["status"] == "optimal" + assert direct_solver["formulation"] == "convex" + assert direct_solver["mip_gap"] is None + assert direct_solver["dpp"] is False + assert direct_solver["cache_hit"] is False + assert direct_solver["model_variables"] > 0 + assert direct_solver["model_constraints"] > 0 + assert reference_solver["engine"] == "cvxpy" + assert reference_solver["backend"] == "highs" + assert reference_solver["formulation"] == "milp" + for key, expected in ( + ("scenario_count", scenario_count), + ("scenario_policy", "shared"), + ("policy_version", "shared-v1"), + ("non_anticipative_slots", 192), + ("cvar_weight", 0.15), + ("cvar_alpha", 0.9), + ): + assert direct_solver[key] == expected + assert reference_solver[key] == expected + assert math.isclose( + direct_solver["service_slack"], + reference_solver["service_slack"], + abs_tol=1e-7, + ) + assert math.isclose( + direct_solver["objective_ore"], + reference_solver["objective_ore"], + abs_tol=1e-3, + ) + assert_nested_close( + direct_solver["objective_breakdown_ore"], + reference_solver["objective_breakdown_ore"], + abs_tol=1e-3, + path="solver.objective_breakdown_ore", + ) + assert_shared_plan_parity(direct, reference, abs_tol=2e-3) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) + + +def test_shared_auto_retries_with_storage_guard_after_direct_failure(monkeypatch) -> None: + from ftw_optimizer import direct_highs, shared_highs + + request = base_request() + request["settings"].update( + { + "mode": "arbitrage", + "formulation": "relaxed", + "shared_backend": "auto", + } + ) + + def reject_direct(*args, **kwargs): + raise direct_highs.DirectHighsError( + "HiGHS returned simultaneous storage charge and discharge" + ) + + monkeypatch.setattr(shared_highs, "solve_direct_highs", reject_direct) + response = handle(request) + + assert response["ok"], response + assert response["solver"]["engine"] == "cvxpy" + assert response["solver"]["formulation"] == "milp" + assert response["solver"]["fallback"] is True + assert "simultaneous" in response["solver"]["fallback_reason"] + assert_storage_replays(request, response) + + +def test_shared_auto_retries_with_storage_guard_after_replay_failure( + monkeypatch, +) -> None: + from ftw_optimizer import shared_highs + from ftw_optimizer.model import ReplayConsistencyError + + request = base_request() + request["settings"].update( + { + "mode": "arbitrage", + "formulation": "relaxed", + "shared_backend": "auto", + } + ) + + def reject_direct(*args, **kwargs): + raise ReplayConsistencyError("direct storage replay failed") + + monkeypatch.setattr(shared_highs, "solve_direct_highs", reject_direct) + response = handle(request) + + assert response["ok"], response + assert response["solver"]["engine"] == "cvxpy" + assert response["solver"]["formulation"] == "milp" + assert response["solver"]["fallback"] is True + assert response["solver"]["fallback_reason"] == "direct storage replay failed" + assert_storage_replays(request, response) + + +def test_shared_auto_retries_generic_direct_failure_without_storage_guard( + monkeypatch, +) -> None: + from ftw_optimizer import shared_highs + + request = base_request() + request["settings"].update( + { + "mode": "arbitrage", + "formulation": "relaxed", + "shared_backend": "auto", + } + ) + + def reject_direct(*args, **kwargs): + raise RuntimeError("direct API failed") + + monkeypatch.setattr(shared_highs, "solve_direct_highs", reject_direct) + response = handle(request) + + assert response["ok"], response + assert response["solver"]["engine"] == "cvxpy" + assert response["solver"]["formulation"] == "convex" + assert response["solver"]["fallback"] is True + assert response["solver"]["fallback_reason"] == "direct API failed" + assert_storage_replays(request, response) + + +def test_shared_auto_retries_other_direct_highs_error_without_storage_guard( + monkeypatch, +) -> None: + from ftw_optimizer import direct_highs, shared_highs + + request = base_request() + request["settings"].update( + { + "mode": "arbitrage", + "formulation": "relaxed", + "shared_backend": "auto", + } + ) + + def reject_direct(*args, **kwargs): + raise direct_highs.DirectHighsError( + "HiGHS economic solve failed with status kTimeLimit" + ) + + monkeypatch.setattr(shared_highs, "solve_direct_highs", reject_direct) + response = handle(request) + + assert response["ok"], response + assert response["solver"]["engine"] == "cvxpy" + assert response["solver"]["formulation"] == "convex" + assert response["solver"]["fallback"] is True + assert "kTimeLimit" in response["solver"]["fallback_reason"] + assert_storage_replays(request, response) + + +def test_shared_auto_lets_cvxpy_reject_invalid_input_after_direct_error() -> None: + request = base_request() + request["settings"].update( + {"formulation": "relaxed", "shared_backend": "auto"} + ) + request["storages"][0]["throughput_cost_ore_kwh"] = "invalid" + + response = handle(request) + + assert not response["ok"] + assert response["error"]["code"] == "invalid_request" + assert "throughput_cost_ore_kwh must be a number" in response["error"]["message"] + + +def test_shared_backend_rejects_unknown_value() -> None: + request = base_request() + request["settings"]["shared_backend"] = "other" + + response = handle(request) + + assert not response["ok"] + assert response["error"]["code"] == "invalid_request" + assert "shared_backend" in response["error"]["message"] + + +def test_shared_backend_highs_rejects_an_ineligible_request() -> None: + request = base_request() + request["settings"]["shared_backend"] = "highs" + request["storages"] = [] + + response = handle(request) + + assert not response["ok"] + assert response["error"]["code"] == "invalid_request" + assert "requires at least one storage" in response["error"]["message"] + + def test_arbitrage_moves_energy_from_cheap_to_expensive_slot() -> None: response = handle(base_request()) assert response["ok"], response @@ -515,8 +1407,49 @@ def test_direct_highs_matches_cvxpy_multistage_reference() -> None: assert direct["ok"], direct assert reference["ok"], reference assert direct["solver"]["engine"] == "highspy" + assert direct["solver"]["backend"] == "highs" + assert direct["solver"]["status"] == "optimal" assert direct["solver"]["formulation"] == "multistage-lp" assert direct["solver"]["dpp"] is False + assert direct["solver"]["cache_hit"] is False + assert direct["solver"]["mip_gap"] is None + assert direct["solver"]["scenario_count"] == 2 + assert direct["solver"]["scenario_original_count"] == 2 + assert direct["solver"]["scenario_reduction_error"] == 0 + assert direct["solver"]["scenario_policy"] == "multistage" + assert direct["solver"]["policy_version"] == "storage-multistage-v1" + assert direct["solver"]["non_anticipative_slots"] == 1 + assert direct["solver"]["tree_nodes"] == 1 + assert direct["solver"]["move_blocks"] == 2 + assert direct["solver"]["decomposition"] == "direct-highs-extensive" + assert direct["solver"]["risk_model"] == "service-cvar-then-expected-cost" + assert direct["solver"]["service_cvar_weight"] == 1 + assert direct["solver"]["service_cvar_alpha"] == 0.95 + assert direct["solver"]["economic_cvar_weight"] == 0.25 + assert direct["solver"]["economic_cvar_alpha"] == 0.9 + assert direct["solver"]["model_variables"] > 0 + assert direct["solver"]["model_constraints"] > 0 + direct_policy = json.loads(direct["solver"]["policy_config"]) + reference_policy = json.loads(reference["solver"]["policy_config"]) + assert direct_policy.pop("backend") == "highs" + assert reference_policy.pop("backend") == "cvxpy" + assert direct_policy == reference_policy + for key in ( + "scenario_count", + "scenario_original_count", + "scenario_reduction_error", + "scenario_policy", + "policy_version", + "non_anticipative_slots", + "tree_nodes", + "move_blocks", + "risk_model", + "service_cvar_weight", + "service_cvar_alpha", + "economic_cvar_weight", + "economic_cvar_alpha", + ): + assert direct["solver"][key] == reference["solver"][key] assert math.isclose( direct["solver"]["objective_ore"], reference["solver"]["objective_ore"], @@ -527,6 +1460,8 @@ def test_direct_highs_matches_cvxpy_multistage_reference() -> None: reference["plan"]["actions"][0]["battery_w"], abs_tol=1e-3, ) + assert_storage_replays(request, direct) + assert_storage_replays(reference_request, reference) def test_multistage_auto_keeps_binary_guards_for_unsafe_incentives() -> None: