#!/usr/bin/env python3 """OPT-11 P4: does the solar forecast pay? Replays recorded days through the REAL Controller with tools/opt11p3_replay.py's engine (house, batteries, prices, peak guard), plus a solar forecast sensor in Home Assistant. python3 tools/opt11p4_replay.py --variant p4_ls --fc recorded python3 tools/opt11p4_replay.py --variant p3_ls --prices real2026 Runs unchanged on release/1.0 7a00f12 (P2) and 5a4cdc2 (P3): an older tree ignores `pv_forecast_entity` and never reads the sensor. HOW THE FORECAST IS MADE (the site's own solar was not recorded for these days; the 2026-09-24 PV recording is only in the live Home Assistant): * The site's solar, hour by hour, is rebuilt from the recording: the house's net use (grid + battery) is house - solar, so solar = max(0, H - net), where H is the median net use of the night hours (22:00-06:00) over all days (`--fc perfect` feeds exactly this). * `recorded`: that solar times the day's error of the real Belgian day-ahead solar forecast (Elia, via energy-charts): forecast / actual, per day. * `national`: the Belgian day-ahead forecast itself, scaled to the site (one factor: the site's total over the national total, all days). It carries local weather differences too, so it is the harder test. * `zero`: a forecast of no sun at all (is the gain the forecast's, or only Leftover solar buying the evening's shortfall?). * `high`: `national` x 1.5. `missing`: the sensor is configured but reads `unavailable`. The sensor has Open-Meteo's shape (state = kWh left today, `wh_period` = Wh per hour); `--shape state` drops `wh_period` (Forecast.Solar's shape). MONEY: as opt11p3_replay.py (all-in import (wholesale + 0.12) x 1.06, export at wholesale), plus `stored_eur`: the change in usable energy (above the 10 % cut-out) from start to end, x 0.95, at the median all-in import price of the run. `net_after_eur` = net - stored_eur. """ import argparse import bisect import json import os import statistics import sys import time from datetime import datetime, timedelta, timezone HERE = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) sys.path.insert(0, HERE) sys.path.insert(0, os.path.join(HERE, "tools")) import opt11p3_replay as rp # noqa: E402 EC = os.path.join(HERE, "scenarios", "energy-charts-be-solar-price-2026-09-01_2026-09-17.json") FC_ENTITY = "sensor.energy_production_today_remaining" FC_TOMORROW = "sensor.energy_production_tomorrow" rp.VARIANTS.update({ "p4_ls": (True, False, {"a": (True, "solar_first", False), "b": (True, "leftover_solar", False)}), "p4_ch": (True, False, {"a": (True, "solar_first", False), "b": (True, "cheap_hours", False)}), }) LOCAL = timezone(timedelta(hours=2)) # CEST for every replayed day def ec(): return json.load(open(EC)) def real2026_points(): """The real day-ahead prices of the replayed days (15-minute).""" d = ec()["price"] return [(datetime.fromtimestamp(t, timezone.utc), 900, p / 1000.0) for t, p in zip(d["unix_seconds"], d["eur_mwh"]) if p is not None] def site_solar(house): """{hour start (UTC epoch): site solar kWh}, and H (W).""" by_h = {} for t, w in house: by_h.setdefault(int(t // 3600) * 3600, []).append(w) net = {h: statistics.fmean(v) for h, v in by_h.items() if len(v) > 20} night = [w for h, w in net.items() if datetime.fromtimestamp(h, LOCAL).hour in (22, 23, 0, 1, 2, 3, 4, 5)] base = statistics.median(night) return {h: max(0.0, base - w) / 1000.0 for h, w in net.items()}, base def national(): s = ec()["solar"] act = dict(zip(s["unix_seconds"], s["actual_mw"])) fc = dict(zip(s["dayahead_unix_seconds"], s["dayahead_mw"])) return act, fc def forecast_table(kind, house, days): """{hour start epoch: forecast kWh}, and a report of its accuracy.""" pv, base = site_solar(house) act, fc = national() day_of = {h: datetime.fromtimestamp(h, LOCAL).date() for h in pv} if kind in ("perfect", "missing"): out = dict(pv) elif kind == "zero": out = {h: 0.0 for h in pv} elif kind == "recorded": ratio = {} for d in set(day_of.values()): hs = [h for h in pv if day_of[h] == d] a = sum(act.get(h) or 0.0 for h in hs) f = sum(fc.get(h) or 0.0 for h in hs) ratio[d] = f / a if a > 0 else 1.0 out = {h: pv[h] * ratio[day_of[h]] for h in pv} else: k = sum(pv.values()) / max(1e-9, sum(act.get(h) or 0.0 for h in pv)) out = {h: (fc.get(h) or 0.0) * k * (1.5 if kind == "high" else 1.0) for h in pv} rep = [] for d in sorted(set(day_of.values())): a = sum(pv[h] for h in pv if day_of[h] == d) f = sum(out[h] for h in out if day_of[h] == d) rep.append((d.isoformat(), round(a, 1), round(f, 1))) return out, base, rep def install_sensor(table, kind, shape): """Give opt11p3_replay's Hass a `state()` for the forecast sensors.""" hours = sorted(table) t0 = rp_start[0] def day_bounds(now): loc = now.astimezone(LOCAL) d0 = loc.replace(hour=0, minute=0, second=0, microsecond=0) return d0, d0 + timedelta(days=1) def sensor(now, tomorrow): d0, d1 = day_bounds(now) if tomorrow: d0, d1 = d1, d1 + timedelta(days=1) a, b = d0.timestamp(), d1.timestamp() i, j = bisect.bisect_left(hours, a), bisect.bisect_left(hours, b) rows = [(h, table[h]) for h in hours[i:j]] if not rows: return None n = now.timestamp() left = 0.0 for h, kwh in rows: if tomorrow or h >= n: left += kwh elif h + 3600 > n: left += kwh * (h + 3600 - n) / 3600.0 attrs = {"unit_of_measurement": "kWh"} if shape == "periods": attrs["wh_period"] = { datetime.fromtimestamp(h, timezone.utc).isoformat(): round(kwh * 1000.0, 1) for h, kwh in rows} # The forecast refreshes every hour. rep = datetime.fromtimestamp(int(n // 3600) * 3600, timezone.utc) return {"state": f"{left:.3f}", "attributes": attrs, "last_updated": rep.isoformat(), "last_reported": rep.isoformat()} async def state(self, e): now = t0 + timedelta(seconds=self.t) if kind == "missing": return {"state": "unavailable", "attributes": {}, "last_updated": now.isoformat()} if e == FC_ENTITY: return sensor(now, False) if e == FC_TOMORROW: return sensor(now, True) return None rp.Hass.state = state rp_start = [None] def main(): ap = argparse.ArgumentParser() ap.add_argument("--variant", default="p4_ls") ap.add_argument("--fc", default="recorded", choices=("none", "recorded", "national", "perfect", "high", "missing", "zero")) ap.add_argument("--prices", default="epex2025", choices=("epex2025", "real2026")) ap.add_argument("--shape", default="periods", choices=("periods", "state")) ap.add_argument("--first", type=int, default=2) ap.add_argument("--days", type=int, default=16) ap.add_argument("--trust", type=int, default=None) ap.add_argument("--report", action="store_true", help="print the forecast's accuracy and exit") a = ap.parse_args() days = [f"2026-09-{d:02d}" for d in range(a.first, a.first + a.days)] house = rp.load_house(days) pts = rp.epex_points() if a.prices == "epex2025" else real2026_points() start = datetime.fromisoformat(days[0] + "T00:00:00+02:00") \ .astimezone(timezone.utc) rp_start[0] = start extra = {} if a.fc != "none": table, base, rep = forecast_table(a.fc, house, days) if a.report: print(json.dumps({"fc": a.fc, "night_base_w": round(base), "days": rep})) return install_sensor(table, a.fc, a.shape) extra = {"pv_forecast_entity": FC_ENTITY, "pv_forecast_tomorrow_entity": FC_TOMORROW} if a.trust is not None: extra["solar_trust_pct"] = a.trust socs = {} def check(c, hass, now, pr): if "start" not in socs: socs["start"] = {u: p.soc for u, p in hass.plants.items()} socs["end"] = {u: p.soc for u, p in hass.plants.items()} dp = getattr(c, "day_plan", None) if dp is not None: socs.setdefault("verdicts", {})[ now.astimezone(LOCAL).date().isoformat()] = dp.verdict t0 = time.time() stats, _c = rp.run(a.variant, house, pts, start, 86400.0 * len(days), extra=extra, check=check) out = rp.summary(stats) prices = sorted(rp.all_in(p) for _s, _n, p in pts if start <= _s < start + timedelta(days=len(days))) med = prices[len(prices) // 2] def usable(soc): return sum(max(0.0, v - 10.0) / 100.0 * rp.CAP_KWH for v in soc.values()) d_kwh = usable(socs["end"]) - usable(socs["start"]) out["stored_change_kwh"] = round(d_kwh, 2) out["stored_eur"] = round(d_kwh * 0.95 * med, 2) out["net_after_eur"] = round(out["net_eur"] - out["stored_eur"], 2) out["end_soc"] = {u: round(v, 1) for u, v in socs["end"].items()} out.update(variant=a.variant, fc=a.fc, prices=a.prices, shape=a.shape, first=days[0], days=len(days), verdicts=socs.get("verdicts", {}), wall_s=round(time.time() - t0, 1)) print(json.dumps(out)) if __name__ == "__main__": main()