goodwe-addon/goodwe_controller/tools/opt11p4_replay.py

#!/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()