Fáze 3.2: solver_v2 — čisté ekonomické jádro plánovače

services/planning/solver_v2.py: MILP s objective = reálné peníze (cash +
degradace − terminal SoC value z DB faktoru). Tvrdá pravidla: bilance, SoC
dynamika, breaker (tvrdý), curtail jen A, GEN cutoff binárka, neg-buy/neg-sell
export bloky, export z baterie ⇒ arb floor (p.19), zákaz současného imp+exp,
EV deadline (placený slack 50 Kč/kWh místo infeasibility), TUV look-ahead,
provozní režimy. SQL masky allow_* vědomě ignorovány (heuristika, ne fyzika).

solver_v2_eval.py: v2 vs v1 na golden fixtures (SoC-fér):
  v2 lepší na VŠECH 5 řešitelných (+231.5 Kč ≈ +22 %), extreme_neg_buy den
  v1=INFEASIBLE → v2 OK (−674.5 Kč). Časy 0.4–10 s (2× na time limitu — TODO).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
Dusan Vojacek
2026-06-11 14:19:32 +02:00
parent 368291e562
commit 90a85b2727
2 changed files with 500 additions and 0 deletions

View File

@@ -0,0 +1,400 @@
# backend/services/planning/solver_v2.py
#
# EMS plánovač v2 — ČISTÉ ekonomické jádro (Fáze 3).
#
# Filozofie: objective = reálné peníze (nákup prodej + degradace terminal
# hodnota energie). Žádné heuristické penalty z constants.py, žádné pre-solver
# fáze/okna/kotvy. Chování (neg-sell příprava, evening export, arbitráž) má
# VYPLYNOUT z cen a fyziky, ne z ručně laděných vah.
#
# Co zůstává (tvrdá pravidla — fyzika, HW, CLAUDE.md):
# - bilance sběrnice, SoC dynamika s účinnostmi, výkonové stropy
# - curtailment jen pole A (pravidlo 5); GEN cutoff binárka pole B (pravidlo 6)
# - block_export_on_negative_sell → ge == 0 při sell < 0 (pravidlo 6, KV1)
# - buy < 0 → ge == 0 (žádná pumpa importexport přes jeden elektroměr; import
# je omezen breakerem — pravidlo 7)
# - export z BATERIE ⇒ koncové SoC ≥ arb floor (pravidlo 19; PV export floor nevynucuje)
# - zákaz současného importu a exportu (binárka)
# - load-first Deye: bc_pv + ge_pv jen z PV přebytku nad zátěží
# - EV deadline, TUV look-ahead, provozní režimy (legitimní constraints)
#
# Vědomé odchylky od v1 (změří harness):
# - SQL masky allow_charge / allow_discharge_export se IGNORUJÍ (jsou to
# výstupy charge-slot-budget heuristik, ne fyzika)
# - EV náklady jen přes bilanci (v1 je účtuje navíc v objective — dvojí započtení)
# - import breaker je tvrdý strop (v1 měkký s 10 Kč/kWh)
# - nedodaná EV energie má explicitní cenu místo infeasibility
from __future__ import annotations
import logging
import time
from typing import Any, Optional
import pulp
from services.planning.constants import (
INTERVAL_H,
SOLVER_TIME_LIMIT,
)
from services.planning.types import (
DispatchResult,
PlanningSlot,
_prague_dow_hour,
)
from services.planning.heuristics import _dispatch_grid_setpoint_w
logger = logging.getLogger(__name__)
V2_BUILD_TAG = "v2-clean-2026-06-11"
# Cena za vypnutí GEN portu (mikroinvertory pole B): reálné riziko/opotřebení
# cyklování stykače — drobná, ale nenulová, aby cutoff platil jen při sell < 0.
V2_GEN_CUTOFF_CZK_KWH = 2.0
# SELF_SUSTAIN: export je nežádoucí, ale tvrdé ge=0 by s neřiditelným polem B
# a plnou baterií bylo infeasible — vysoká cena funguje jako ventil.
V2_SELF_SUSTAIN_EXPORT_CZK_KWH = 100.0
# Cena nedodané EV energie do deadline (Kč/kWh) — místo tvrdé infeasibility.
V2_EV_UNMET_CZK_KWH = 50.0
# Nepatrný tie-break proti zbytečnému curtailu při cenové indiferenci (Kč/kWh).
V2_CURTAIL_TIEBREAK_CZK_KWH = 0.001
def _terminal_value_czk_per_wh(slots: list[PlanningSlot], battery: Any) -> float:
"""Shadow cena zbytkové energie: průměrný buy prvních 24 h × DB faktor (pravidlo 16)."""
n24 = min(len(slots), int(24 / INTERVAL_H))
avg_buy = sum(float(s.buy_price) for s in slots[:n24]) / max(1, n24)
factor = float(getattr(battery, "planner_terminal_soc_value_factor", 1.0) or 1.0)
return max(0.0, avg_buy) * factor / 1000.0
def _arb_floor_wh(battery: Any) -> float:
"""Podlaha SoC pro export z baterie (pravidlo 19): ekonomická rezerva z DB."""
floor = getattr(battery, "arb_floor_wh", None)
if floor is None:
floor = getattr(battery, "reserve_soc_wh", None)
return max(float(floor or 0.0), float(battery.min_soc_wh))
def solve_dispatch_v2(
slots: list[PlanningSlot],
battery: Any,
heat_pump: Any,
grid: Any,
ev_sessions: list,
vehicles: list,
current_soc_wh: float,
current_tuv_temp_c: float,
*,
tuv_delta_stats: Optional[dict[tuple[int, int], float]] = None,
operating_mode: str = "AUTO",
planner_version: str | None = None,
) -> tuple[list[DispatchResult], int, dict[str, Any]]:
"""Čistý ekonomický MILP; rozhraní kompatibilní se solve_dispatch (v1)."""
if not slots:
raise RuntimeError("solve_dispatch_v2 requires at least one slot")
t0 = time.monotonic()
T = len(slots)
om = (operating_mode or "AUTO").upper()
EV = min(len(vehicles), 2)
max_imp = float(grid.max_import_power_w)
max_exp = float(grid.max_export_power_w)
max_chg = float(battery.max_charge_power_w)
max_dis = float(battery.max_discharge_power_w)
eff_c = float(battery.charge_efficiency)
eff_d = float(battery.discharge_efficiency)
deg = float(battery.degradation_cost_czk_kwh)
soc_min = float(battery.min_soc_wh)
soc_max = float(battery.soc_max_wh)
usable = float(battery.usable_capacity_wh)
arb_floor = _arb_floor_wh(battery)
terminal = _terminal_value_czk_per_wh(slots, battery)
block_neg_sell = bool(getattr(grid, "block_export_on_negative_sell", False))
gen_cutoff_avail = bool(getattr(grid, "deye_gen_microinverter_cutoff_enabled", False))
soc0 = min(max(float(current_soc_wh), soc_min), soc_max)
prob = pulp.LpProblem("dispatch_v2", pulp.LpMinimize)
gi = [pulp.LpVariable(f"gi_{t}", 0, max_imp) for t in range(T)]
ge_pv = [pulp.LpVariable(f"gepv_{t}", 0, max_exp) for t in range(T)]
ge_bat = [pulp.LpVariable(f"gebat_{t}", 0, max_exp) for t in range(T)]
bc_pv = [pulp.LpVariable(f"bcpv_{t}", 0, max_chg) for t in range(T)]
bc_gi = [pulp.LpVariable(f"bcgi_{t}", 0, max_chg) for t in range(T)]
bd = [pulp.LpVariable(f"bd_{t}", 0, max_dis) for t in range(T)]
ca = [pulp.LpVariable(f"ca_{t}", 0, max(0, int(slots[t].pv_a_forecast_w))) for t in range(T)]
soc = [pulp.LpVariable(f"soc_{t}", soc_min, soc_max) for t in range(T)]
hp = [pulp.LpVariable(f"hp_{t}", 0, float(heat_pump.rated_heating_power_w)) for t in range(T)]
y_imp = [pulp.LpVariable(f"yimp_{t}", cat=pulp.LpBinary) for t in range(T)]
z_exp = [pulp.LpVariable(f"zexp_{t}", cat=pulp.LpBinary) for t in range(T)]
z_gen = (
[pulp.LpVariable(f"zgen_{t}", cat=pulp.LpBinary) for t in range(T)]
if gen_cutoff_avail
else None
)
ev_direct = [
[
pulp.LpVariable(f"evd_{e}_{t}", 0, min(float(vehicles[e].max_charge_power_w), max_imp))
for t in range(T)
]
for e in range(EV)
]
ev_via_bat = [
[
pulp.LpVariable(f"evb_{e}_{t}", 0, float(vehicles[e].max_charge_power_w))
for t in range(T)
]
for e in range(EV)
]
ev_unmet: list = [] # slack Wh per session (cena V2_EV_UNMET_CZK_KWH)
def _connected(e: int, t: int) -> bool:
return bool(slots[t].ev1_connected if e == 0 else slots[t].ev2_connected)
for t in range(T):
s = slots[t]
pv_a = max(0.0, float(s.pv_a_forecast_w))
pv_b = max(0.0, float(s.pv_b_forecast_w))
pv_a_net = pv_a - ca[t]
pv_b_eff = pv_b - (pv_b * z_gen[t] if z_gen is not None else 0.0)
ev_total_t = pulp.lpSum(
ev_direct[e][t] + ev_via_bat[e][t] for e in range(EV)
)
load_site = float(s.load_baseline_w) + ev_total_t + hp[t]
# bilance sběrnice (W)
prob += (
pv_a_net + pv_b_eff + gi[t] + bd[t]
== load_site + bc_pv[t] + bc_gi[t] + ge_pv[t] + ge_bat[t]
), f"balance_{t}"
# SoC dynamika (Wh)
prev = soc0 if t == 0 else soc[t - 1]
prob += (
soc[t]
== prev
+ (bc_pv[t] + bc_gi[t]) * eff_c * INTERVAL_H
- bd[t] / eff_d * INTERVAL_H
), f"soc_{t}"
# výkonové stropy
prob += bc_pv[t] + bc_gi[t] <= max_chg, f"chg_cap_{t}"
prob += ge_pv[t] + ge_bat[t] <= max_exp, f"exp_cap_{t}"
# PV cesty omezené dostupnou výrobou (load-first vynucuje HW; bilance účtuje energii)
prob += bc_pv[t] + ge_pv[t] <= pv_a_net + pv_b_eff, f"pv_src_{t}"
# bc_gi jen ze sítě:
prob += bc_gi[t] <= gi[t], f"bcgi_src_{t}"
# vybíjení kryje dům + EV-via-bat + export z baterie
prob += ge_bat[t] + pulp.lpSum(ev_via_bat[e][t] for e in range(EV)) <= bd[t], f"bd_split_{t}"
# zákaz současného importu a exportu
prob += gi[t] <= max_imp * y_imp[t], f"imp_excl_{t}"
prob += ge_pv[t] + ge_bat[t] <= max_exp * (1 - y_imp[t]), f"exp_excl_{t}"
# pravidlo 19: export z baterie ⇒ SoC ≥ arb floor
prob += ge_bat[t] <= max_exp * z_exp[t], f"zexp_link_{t}"
prob += soc[t] >= arb_floor - (soc_max - soc_min) * (1 - z_exp[t]), f"zexp_floor_{t}"
# tvrdá cenová pravidla
if float(s.buy_price) < 0.0:
prob += ge_pv[t] + ge_bat[t] == 0, f"neg_buy_noexp_{t}"
if float(s.sell_price) < 0.0 and block_neg_sell:
prob += ge_pv[t] + ge_bat[t] == 0, f"neg_sell_block_{t}"
# EV dostupnost
for e in range(EV):
if not _connected(e, t):
prob += ev_direct[e][t] == 0
prob += ev_via_bat[e][t] == 0
else:
prob += ev_direct[e][t] + ev_via_bat[e][t] <= float(
vehicles[e].max_charge_power_w
)
# provozní režimy (tvrdé constraints dle operating-modes.md)
if om == "SELF_SUSTAIN":
prob += gi[t] <= float(s.load_baseline_w), f"ss_gi_{t}"
elif om == "PRESERVE":
prob += bc_pv[t] == 0
prob += bc_gi[t] == 0
prob += bd[t] == 0
elif om == "CHARGE_CHEAP":
prob += ge_pv[t] + ge_bat[t] == 0
prob += bd[t] == 0
# EV deadline (s placeným slackem místo infeasibility)
for e in range(EV):
sess = ev_sessions[e] if e < len(ev_sessions) else None
if sess is None or not getattr(sess, "energy_needed_wh", 0):
continue
t_dl = next(
(t for t in range(T) if slots[t].interval_start >= sess.target_deadline),
T - 1,
)
unmet = pulp.LpVariable(f"ev_unmet_{e}", 0, float(sess.energy_needed_wh))
ev_unmet.append(unmet)
prob += (
pulp.lpSum(
(ev_direct[e][t] + ev_via_bat[e][t]) * INTERVAL_H
for t in range(t_dl + 1)
if _connected(e, t)
)
+ unmet
>= float(sess.energy_needed_wh)
), f"ev_deadline_{e}"
# TUV look-ahead (převzato z v1 — komfortní constraint, ne heuristika)
rated_hp = float(heat_pump.rated_heating_power_w)
if tuv_delta_stats and rated_hp > 0 and getattr(heat_pump, "tuv_min_temp_c", None):
tuv_pred = float(current_tuv_temp_c)
tgt = float(getattr(heat_pump, "tuv_target_temp_c", 55.0) or 55.0)
thr = float(heat_pump.tuv_min_temp_c) + 5.0
for t in range(T):
dow, hour = _prague_dow_hour(slots[t].interval_start)
delta = tuv_delta_stats.get((dow, hour), -0.1)
tuv_pred += float(delta) * INTERVAL_H
if tuv_pred < thr:
prob += (
pulp.lpSum(hp[s_] for s_ in range(max(0, t - 8), t + 1))
>= rated_hp * 0.5
), f"tuv_heat_{t}"
tuv_pred = tgt
if float(current_tuv_temp_c) < float(heat_pump.tuv_min_temp_c):
prob += hp[0] >= rated_hp * 0.8, "tuv_emergency"
# ---------------- objective: jen reálné peníze ----------------
wh = INTERVAL_H / 1000.0 # W → kWh za slot
cash = pulp.lpSum(
gi[t] * float(slots[t].buy_price) * wh
- (ge_pv[t] + ge_bat[t]) * float(slots[t].sell_price) * wh
for t in range(T)
)
degradation = pulp.lpSum(
0.5 * (bc_pv[t] + bc_gi[t] + bd[t]) * deg * wh for t in range(T)
)
extras = pulp.lpSum(ca[t] * V2_CURTAIL_TIEBREAK_CZK_KWH * wh for t in range(T))
if z_gen is not None:
extras += pulp.lpSum(
max(0.0, float(slots[t].pv_b_forecast_w)) * z_gen[t] * V2_GEN_CUTOFF_CZK_KWH * wh
for t in range(T)
)
if om == "SELF_SUSTAIN":
extras += pulp.lpSum(
(ge_pv[t] + ge_bat[t]) * V2_SELF_SUSTAIN_EXPORT_CZK_KWH * wh for t in range(T)
)
if ev_unmet:
extras += pulp.lpSum(u * V2_EV_UNMET_CZK_KWH / 1000.0 for u in ev_unmet)
prob += cash + degradation + extras - terminal * soc[T - 1]
solver = (
pulp.HiGHS_CMD(msg=False, timeLimit=SOLVER_TIME_LIMIT)
if pulp.HiGHS_CMD().available()
else pulp.PULP_CBC_CMD(msg=False, timeLimit=SOLVER_TIME_LIMIT)
)
status = prob.solve(solver)
duration_ms = int((time.monotonic() - t0) * 1000)
status_str = pulp.LpStatus[status]
if status_str != "Optimal":
# v2 nemá relax řetězec — model je navržen tak, aby byl feasible
# (placené slacky místo tvrdých kotev). Ne-Optimal je skutečná chyba.
raise RuntimeError(f"solver_v2: {status_str}")
# ---------------- DispatchResult assembly (parita s v1) ----------------
def _val(var) -> float:
v = pulp.value(var)
return float(v) if v is not None else 0.0
results: list[DispatchResult] = []
for t in range(T):
s = slots[t]
bc_tot = _val(bc_pv[t]) + _val(bc_gi[t])
bd_v = _val(bd[t])
batt_w = round(bc_tot - bd_v)
ge_pv_w = round(_val(ge_pv[t]))
ge_bat_w = round(_val(ge_bat[t]))
gi_w = _val(gi[t])
ge_w = float(ge_pv_w + ge_bat_w)
grid_w, export_mode = _dispatch_grid_setpoint_w(
gi_w=gi_w,
ge_w=ge_w,
ge_bat_w=float(ge_bat_w),
ge_pv_w=float(ge_pv_w),
max_export_power_w=int(max_exp),
)
if batt_w < 0 and grid_w < 0:
deye_mode = "SELL"
elif batt_w > 0 and grid_w > 0:
deye_mode = "CHARGE"
else:
deye_mode = "PASSIVE"
gen_cut = bool(round(_val(z_gen[t]))) if z_gen is not None else None
hp_v = _val(hp[t])
hp_on = hp_v > rated_hp * 0.5 if rated_hp > 0 else False
cash_t = gi_w * float(s.buy_price) * wh - ge_w * float(s.sell_price) * wh
pen_t = 0.0
if gen_cut:
pen_t += max(0.0, float(s.pv_b_forecast_w)) * V2_GEN_CUTOFF_CZK_KWH * wh
results.append(
DispatchResult(
interval_start=s.interval_start,
battery_setpoint_w=batt_w,
battery_soc_target=round(_val(soc[t]) / usable * 100.0, 2),
grid_setpoint_w=grid_w,
export_limit_w=int(max_exp) if grid_w < 0 else 0,
export_mode=export_mode,
deye_physical_mode=deye_mode,
deye_gen_cutoff_enabled=gen_cut,
ev1_setpoint_w=(
round(_val(ev_direct[0][t]) + _val(ev_via_bat[0][t]))
if EV > 0 and s.ev1_connected
else None
),
ev2_setpoint_w=(
round(_val(ev_direct[1][t]) + _val(ev_via_bat[1][t]))
if EV > 1 and s.ev2_connected
else None
),
ev1_via_bat_w=round(_val(ev_via_bat[0][t])) if EV > 0 else 0,
ev2_via_bat_w=round(_val(ev_via_bat[1][t])) if EV > 1 else 0,
heat_pump_enabled=hp_on,
heat_pump_setpoint_w=int(rated_hp) if hp_on else 0,
pv_a_curtailed_w=round(_val(ca[t])),
expected_cost_czk=round(cash_t, 4),
effective_buy_price=float(s.buy_price),
effective_sell_price=float(s.sell_price),
is_predicted_price=bool(s.is_predicted_price),
cashflow_czk=round(cash_t, 4),
battery_arbitrage_czk=0.0,
penalty_czk=round(pen_t, 4),
green_bonus_czk=float(getattr(s, "green_bonus_czk_per_slot", 0.0) or 0.0),
)
)
snapshot: dict[str, Any] = {
"version": planner_version or "v2-clean",
"planner_build_tag": V2_BUILD_TAG,
"inputs": {
"operating_mode": om,
"current_soc_wh": soc0,
"terminal_czk_per_wh": round(terminal, 8),
"arb_floor_wh": arb_floor,
"block_export_on_negative_sell": block_neg_sell,
"gen_cutoff_available": gen_cutoff_avail,
"slot_count": T,
"ev_sessions": sum(1 for x in ev_sessions if x is not None),
"masks_ignored": True,
},
"objective_terms": {
"cash_czk": round(float(pulp.value(cash)), 3),
"degradation_czk": round(float(pulp.value(degradation)), 3),
"extras_czk": round(float(pulp.value(extras)), 3) if not isinstance(extras, float) else 0.0,
"terminal_value_czk": round(terminal * _val(soc[T - 1]), 3),
"ev_unmet_wh": [round(_val(u), 1) for u in ev_unmet],
},
"solver_duration_ms": duration_ms,
"solver_status": status_str,
}
return results, duration_ms, snapshot

View File

@@ -0,0 +1,100 @@
#!/usr/bin/env python3
"""
Fáze 3 vyhodnocení solver_v2 (čisté jádro) proti v1 na golden fixtures.
Replay STEJNOU cestou jako golden gate (_load_site_context + _load_slots nad
FixtureDB), ale přes services.planning.solver_v2.solve_dispatch_v2. Porovnání
s golden snapshoty v1 (SoC-fér: koncový SoC obou oceněn terminal cenou v2).
Spouštět z backend/: python3 ../scripts/harness/solver_v2_eval.py
"""
from __future__ import annotations
import asyncio
import importlib.util
import json
import sys
from datetime import datetime
from pathlib import Path
BACKEND = Path(__file__).resolve().parents[2] / "backend"
sys.path.insert(0, str(BACKEND))
from services import planning_engine as pe # noqa: E402
from services.planning import solver_v2 as v2 # noqa: E402
_spec = importlib.util.spec_from_file_location(
"golden_replay", BACKEND / "tests" / "test_golden_replay.py"
)
_golden = importlib.util.module_from_spec(_spec)
sys.modules["golden_replay"] = _golden
_spec.loader.exec_module(_golden)
FIXTURES = sorted((BACKEND / "tests" / "golden" / "fixtures").glob("*.json"))
SNAPSHOTS = BACKEND / "tests" / "golden" / "snapshots"
def _replay_v2(fixture: dict):
async def _run():
db = _golden._FixtureDB(fixture)
meta = fixture["meta"]
(battery, heat_pump, grid, vehicles, ev_sessions, soc_wh, tuv_temp,
operating_mode, tuv_stats) = await pe._load_site_context(int(meta["site_id"]), db)
slots = await pe._load_slots(
int(meta["site_id"]),
datetime.fromisoformat(meta["window_from"]),
datetime.fromisoformat(meta["window_to"]),
db,
soc_wh=soc_wh,
)
results, ms, snap = v2.solve_dispatch_v2(
slots, battery, heat_pump, grid, ev_sessions, vehicles,
soc_wh, tuv_temp,
tuv_delta_stats=tuv_stats,
operating_mode=operating_mode or "AUTO",
)
return results, ms, snap, battery
return asyncio.run(_run())
def main() -> None:
header = f"{'fixture':<42} {'v1':>9} {'v2':>9} {'Δ':>8} {'v2 ms':>6} pozn."
print("# solver_v2 vs v1 — modelovaný cashflow, SoC-fér (Kč/horizont; Δ<0 = v2 lepší)")
print()
print(header)
print("-" * len(header))
tot1 = tot2 = 0.0
solved_both = 0
for path in FIXTURES:
fixture = json.loads(path.read_text(encoding="utf-8"))
try:
results, ms, snap, battery = _replay_v2(fixture)
except Exception as exc:
print(f"{path.stem:<42} {'?':>9} {'CHYBA':>9} {'':>8} {exc}")
continue
usable = float(battery.usable_capacity_wh)
term = float(snap["inputs"]["terminal_czk_per_wh"])
v2_cash = sum(r.cashflow_czk for r in results)
v2_soc_end = results[-1].battery_soc_target / 100.0 * usable
v2_adj = v2_cash - v2_soc_end * term
snap1 = json.loads((SNAPSHOTS / path.name).read_text(encoding="utf-8"))
if "solver_error" in snap1:
print(f"{path.stem:<42} {'INFEAS':>9} {v2_adj:>9.1f} {'':>8} {ms:>6} v1 selhal, v2 OK")
continue
v1_cash = snap1["totals"]["cashflow_czk"]
v1_soc_end = snap1["slots"][-1]["battery_soc_target"] / 100.0 * usable
v1_adj = v1_cash - v1_soc_end * term
d = v2_adj - v1_adj
tot1 += v1_adj
tot2 += v2_adj
solved_both += 1
print(f"{path.stem:<42} {v1_adj:>9.1f} {v2_adj:>9.1f} {d:>8.1f} {ms:>6}")
print("-" * len(header))
if solved_both:
print(f"{'CELKEM (oba řešitelné)':<42} {tot1:>9.1f} {tot2:>9.1f} {tot2 - tot1:>8.1f}")
if __name__ == "__main__":
main()