Separate model accuracy from forecast accuracy

The Forecast Error sensors replayed the model with the relay, electrical heat
and window history that actually followed, so they measured the model, not a
forecast anyone could have made. They are renamed Model Error (unique IDs
unchanged) and gain RMSE, effective sample count and a quality label; their
state stays unknown until there are at least 5 roughly independent spans.

A new Forecast Error sensor scores forecasts as published, with the plan and
assumptions of the moment they were made. A forecast is scored only when the
relay followed its plan; otherwise it is counted as counterfactual, not as
error. Statistics are saved with the model and dropped with a zone's
recalibration.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

custom_components/heating_mpc/coordinator.py

@@ -2,6 +2,8 @@
 from __future__ import annotations
 
 import asyncio
+import hashlib
+import json
 import logging
 from collections import deque
 from datetime import datetime, timedelta, timezone
@@ -23,7 +25,7 @@ from .const import (
 )
 from .compatibility import FINGERPRINT_VERSION, baseline_fingerprints, changed, reason, zone_fingerprints
 from .data_quality import ReportAgeHistory, freshness, inputs_differ
-from .evaluation import ForecastEvaluator, model_snapshot
+from .evaluation import ForecastEvaluator, OperationalForecastLedger, model_snapshot
 from .heat_sources import SUBENTRY_HEAT_SOURCE, HeatSource, HeatSourceTracker, build_sources, normalize_power
 from .learning import FeatureLearner
 from .physics import (
@@ -129,6 +131,9 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
         self._retired: dict[str, dict[str, Any]] = {}
         self._report_ages = ReportAgeHistory()
         self._evaluator = ForecastEvaluator(LAG_LEARN_MAX_UNCERTAINTY)
+        self._forecasts = OperationalForecastLedger()
+        # Seconds the relay was on during the last interval (None if unknown), and its length.
+        self._relay_interval: tuple[float | None, float | None] = (None, None)
         # How each poll was handled since start, to explain slow learning.
         self._sample_results: dict[str, int] = {}
         self._store = Store(hass, 1, f"{DOMAIN}.{entry.entry_id}")
@@ -167,6 +172,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             _LOGGER.warning("Heating MPC heater changed; starting a fresh thermal model")
             return
         self._evaluator.load(stored.get("evaluation"))
+        self._forecasts.load(stored.get("operational_evaluation"))
         changes = stored.get("model_changes")
         self._model_changes.extend(c for c in (changes if isinstance(changes, list) else []) if isinstance(c, dict))
         self._heat_sources = HeatSourceTracker(self.sources, stored.get("heat_sources"), self._baseline_prints)
@@ -227,6 +233,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
                 # The zone temperature is defined differently: re-learn it from a warm start.
                 self.learners[room].recalibrate()
                 self._evaluator.forget(room)
+                self._forecasts.forget(room)
                 resets.append("base")
             for layer in (self.extended.get(room), self.ventilation.get(room)):
                 if layer is None:
@@ -265,6 +272,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             "ventilation_models": {key: model.as_dict() for key, model in self.ventilation.items()},
             "ventilation": self._ventilation.as_dict(),
             "evaluation": self._evaluator.as_dict(),
+            "operational_evaluation": self._forecasts.as_dict(),
             "fingerprint_version": FINGERPRINT_VERSION,
             "fingerprints": self._fingerprints,
             "model_changes": list(self._model_changes),
@@ -316,6 +324,7 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
         tracker = self._lag
         steady = False
         lag = uncertainty = None
+        self._relay_interval = (None, None)  # unknown unless the relay history is known below
         if tracker is not None and now > tracker["time"]:
             start = tracker["time"]
             changes = [(when, state) for when, state in self._relay_events if start < when <= now]
@@ -326,10 +335,13 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             segments = None if self._relay_events_dropped else relay_segments(
                 tracker["heater_on"], start.timestamp(),
                 [(when.timestamp(), state) for when, state in changes], now.timestamp())
+            self._relay_interval = (None, (now - start).total_seconds())
             if segments is not None and segments[-1][0] == heater_on:
                 lag = tracker["lag"]
                 for command, seconds in segments:
                     lag = advance_lag(lag, command, seconds)
+                self._relay_interval = (sum(sec for command, sec in segments if command),
+                                        (now - start).total_seconds())
                 # The lag update is a contraction, so an initial error decays the same way.
                 uncertainty = tracker["uncertainty"] * exp(
                     -(now - start).total_seconds() / LAG_TAU_SECONDS)
@@ -343,6 +355,10 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
         self._relay_events_dropped = False
         return steady
 
+    def _revision(self) -> str:
+        """Short identity of the model configuration a forecast was made with."""
+        return hashlib.sha1(json.dumps(self._fingerprints, sort_keys=True).encode()).hexdigest()[:10]
+
     def _superseded(self, status: dict[str, Any]) -> dict[str, Any]:
         """Result of an update that a newer one (or unloading) overtook: keep what is shown."""
         current = getattr(self, "data", None)
@@ -612,6 +628,12 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             self._lag = None
             self._relay_events.clear()
             self._relay_events_dropped = False
+            self._relay_interval = (None, None)
+        # Published forecasts that are due are scored against these readings; nothing
+        # here is fed back into what those forecasts assumed.
+        self._forecasts.observe(now.timestamp(), {k: v for k, v in temps.items() if v is not None},
+                                *self._relay_interval)
+        status["operational_forecast"] = self._forecasts.summary(self.rooms)
         # Rooms without a usable reading are left out and the others carry on;
         # only the shared inputs, or losing every room, pause the integration.
         vals: dict[str, float] = {room: float(t) for room, t in temps.items() if t is not None}
@@ -782,12 +804,21 @@ class HeatingMPCCoordinator(DataUpdateCoordinator[dict[str, Any]]):
             _LOGGER.warning("MPC prediction failed: %s", err)
             status["status"] = "prediction_error"
             return status
+        self._forecasts.record(now.timestamp(), vals, output["trajectory_15min"], output["plan_15min"],
+                               self._revision(), {
+            "strategy": "mpc_plan", "outside": outside, "wind": weather.get("wind"),
+            "internal_half_life_minutes": float(self.config.get(CONF_INTERNAL_HALF_LIFE, 30)),
+            "window_open_minutes": float(self.config[CONF_WINDOW_OPEN_MINUTES]),
+            "internal_kw": {room: internal[room] for room in vals},
+            "windows_open": [w.name for w in self.windows
+                             if self._ventilation.window_state.get(w.key, {}).get("open")]})
         status.update({
             "status": "ready", "proposal": output["proposal"],
             "prediction_60m": output["predicted_60m"],
             "prediction_end": output["predicted_end"],
             "plan_15min": output["plan_15min"], "objective": output["cost"],
             "prediction_trajectory": output["trajectory_15min"],
-            "model_quality_note": "5-minute in-sample MAE; not a validated 60-minute forecast error",
+            "model_quality_note": ("mae_5m is the in-sample 5-minute error; Model Error sensors replay the "
+                                   "model with actual inputs; Forecast Error scores published forecasts"),
         })
         return status

custom_components/heating_mpc/evaluation.py

@@ -9,18 +9,28 @@ then (so later data is never used), the outdoor temperature held at its starting
 value, and the radiator lag that actually followed from the relay's real ON/OFF
 history. The replayed temperature is compared with the reading now.
 
-This scores the thermal model, not the MPC's proposed plan: Generic Thermostat
-controls the relay, so the plan is never executed. Two baselines are scored on
-the same spans: persistence (no change) and the base 3-parameter model. The full
-model replays all zones together with the electrical heat, neighbour exchange and
-window openings that actually happened, so its improvement over the base model
-shows whether those features earn their keep.
+Two kinds of accuracy are kept apart:
+
+* ForecastEvaluator, **model accuracy** (retrospective): the model is replayed
+  with the relay, electrical heat and window history that actually happened. It
+  answers "how well does the model reproduce the rooms under known inputs?" and
+  compares the full model with the base model and with persistence. It is not a
+  forecast accuracy: a real forecast does not know the future inputs.
+* OperationalForecastLedger, **forecast accuracy**: forecasts as published, with
+  the boiler plan, electrical heat and windows they assumed, compared with what
+  happened once their horizon has passed. A forecast assuming the MPC's plan can
+  only be checked when the relay actually followed that plan; otherwise it is
+  counted as counterfactual and not scored.
+
+Forecasts start every 5 minutes, so consecutive ones overlap: n forecasts at a
+horizon of h minutes carry roughly n x 5 / h independent samples. Metrics report
+that number and a quality label, and sensors show no value while it is too small.
 """
 from __future__ import annotations
 
 from collections import deque
-from datetime import datetime, timedelta
-from math import isfinite
+from datetime import datetime, timedelta, timezone
+from math import isfinite, sqrt
 from typing import Any
 
 from .physics import Learner
@@ -30,6 +40,20 @@ from .zone_model import feature_vector
 HORIZONS = (15, 30, 60, 120)  # minutes
 WINDOW_DAYS = 7
 STEP_SECONDS = 300
+# Independent samples (overlap-adjusted) needed before a metric is shown / trusted.
+MIN_EFFECTIVE_SAMPLES = 5
+GOOD_EFFECTIVE_SAMPLES = 30
+
+
+def effective_samples(count: int, horizon: int, zones: int = 1) -> float:
+    """Roughly independent samples among `count` forecasts started 5 minutes apart."""
+    return round(count / max(1, zones) * min(1.0, STEP_SECONDS / 60 / horizon), 1)
+
+
+def quality(effective: float) -> str:
+    if effective < MIN_EFFECTIVE_SAMPLES:
+        return "insufficient"
+    return "indicative" if effective < GOOD_EFFECTIVE_SAMPLES else "ok"
 
 
 def replay(temp: float, outside: float, theta: tuple[float, float, float],
@@ -94,8 +118,8 @@ class ForecastEvaluator:
     def __init__(self, max_lag_uncertainty: float) -> None:
         self._max_lag_uncertainty = max_lag_uncertainty
         self._history: deque[dict[str, Any]] = deque(maxlen=HORIZONS[-1] * 60 // STEP_SECONDS + 1)
-        # UTC day -> horizon -> room ->
-        #   [count, sum error, sum |error|, sum |persistence error|, sum |base model error|]
+        # UTC day -> horizon -> room -> [count, sum error, sum |error|, sum |persistence error|,
+        #   sum |base model error|, sum error^2, count with squares (0 for 0.3/0.4 data)]
         self._days: dict[str, dict[str, dict[str, list[float]]]] = {}
 
     def reset_window(self) -> None:
@@ -146,12 +170,14 @@ class ForecastEvaluator:
                 error = full[room] - actual
                 base_error = replay(start, origin["outside"], theta, inputs) - actual
                 bucket = self._days.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(
-                    room, [0, 0.0, 0.0, 0.0, 0.0])
+                    room, [0, 0.0, 0.0, 0.0, 0.0, 0.0, 0])
                 bucket[0] += 1
                 bucket[1] += error
                 bucket[2] += abs(error)
                 bucket[3] += abs(start - actual)
                 bucket[4] += abs(base_error)
+                bucket[5] += error * error
+                bucket[6] += 1
 
     def _prune(self, time: datetime) -> None:
         oldest = (time.date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
@@ -166,12 +192,12 @@ class ForecastEvaluator:
             for day in self._days.values():
                 for room, values in day.get(str(horizon), {}).items():
                     if room in rooms:
-                        totals = per_room.setdefault(room, [0, 0.0, 0.0, 0.0, 0.0])
+                        totals = per_room.setdefault(room, [0, 0.0, 0.0, 0.0, 0.0, 0.0, 0])
                         for i, value in enumerate(values):
                             totals[i] += value
-            overall = [sum(values[i] for values in per_room.values()) for i in range(5)]
-            result[horizon] = {**_stats(overall),
-                               "rooms": {room: _stats(values) for room, values in per_room.items()}}
+            overall = [sum(values[i] for values in per_room.values()) for i in range(7)]
+            result[horizon] = {**_stats(overall, horizon, len(per_room)),
+                               "rooms": {room: _stats(values, horizon) for room, values in per_room.items()}}
         return result
 
     def forget(self, room: str) -> None:
@@ -184,7 +210,7 @@ class ForecastEvaluator:
 
     def as_dict(self) -> dict[str, Any]:
         # Rounded copy: shorter JSON on every save; in-memory sums stay exact.
-        return {"days": {day: {horizon: {room: [v[0], *(round(x, 6) for x in v[1:])]
+        return {"days": {day: {horizon: {room: [v[0], *(round(x, 6) for x in v[1:6]), v[6]]
                                          for room, v in by_room.items()}
                                for horizon, by_room in by_horizon.items()}
                          for day, by_horizon in self._days.items()}}
@@ -209,28 +235,38 @@ class ForecastEvaluator:
                 for room, values in by_room.items():
                     if _valid_bucket(room, values):
                         # 0.3 buckets had no base column: the model then was the base model.
-                        base = values[4] if len(values) == 5 else values[2]
+                        # Older buckets have no squared errors: they do not count for RMSE.
+                        base = values[4] if len(values) >= 5 else values[2]
+                        squares = (float(values[5]), values[6]) if len(values) == 7 else (0.0, 0)
                         self._days.setdefault(day, {}).setdefault(horizon, {})[room] = [
-                            values[0], float(values[1]), float(values[2]), float(values[3]), float(base)]
+                            values[0], float(values[1]), float(values[2]), float(values[3]), float(base),
+                            *squares]
 
 
 def _valid_bucket(room: Any, values: Any) -> bool:
-    if not isinstance(room, str) or not isinstance(values, list) or len(values) not in (4, 5):
+    if not isinstance(room, str) or not isinstance(values, list) or len(values) not in (4, 5, 7):
+        return False
+    count, *sums = values[:6]
+    if len(values) == 7 and (type(values[6]) is not int or not 0 <= values[6] <= count):
         return False
-    count, *sums = values
     return (type(count) is int and count >= 0
             and all(type(v) in (int, float) and isfinite(v) for v in sums)
             and all(v >= 0 for v in sums[1:]))
 
 
-def _stats(values: list[float]) -> dict[str, Any]:
-    count, error, absolute, persistence, base = values
+def _stats(values: list[float], horizon: int, zones: int = 1) -> dict[str, Any]:
+    count, error, absolute, persistence, base, squares, square_count = values
+    effective = effective_samples(int(count), horizon, zones)
     if not count:
-        return {"samples": 0, "mae": None, "bias": None, "persistence_mae": None, "skill": None,
+        return {"samples": 0, "effective_samples": 0.0, "quality": "insufficient", "mae": None,
+                "rmse": None, "bias": None, "persistence_mae": None, "skill": None,
                 "base_model_mae": None, "improvement_vs_base": None}
     return {
         "samples": int(count),
+        "effective_samples": effective,
+        "quality": quality(effective),
         "mae": round(absolute / count, 3),
+        "rmse": round(sqrt(squares / square_count), 3) if square_count else None,
         "bias": round(error / count, 3),
         "persistence_mae": round(persistence / count, 3),
         # Share of the persistence error removed by the model (1 = perfect, < 0 = worse).
@@ -240,3 +276,191 @@ def _stats(values: list[float]) -> dict[str, Any]:
         # ventilation (0 = no gain, < 0 = the extra features make it worse).
         "improvement_vs_base": round(1 - absolute / base, 3) if base > 0 else None,
     }
+
+
+# A forecast's plan counts as executed when, in every 15-minute block, the share
+# of time the relay was on is within this of what the plan assumed (0 or 1).
+PLAN_TOLERANCE = 0.2
+DUE_TOLERANCE_SECONDS = 150  # half a poll: the poll closest to the forecast's due time
+VERDICTS = ("counterfactual", "unverifiable", "missing")
+
+
+class OperationalForecastLedger:
+    """Forecasts as published, scored against what happened once their horizon passed.
+
+    Nothing is replaced with hindsight: the start temperatures, trajectory, boiler
+    plan and assumptions are those of the moment the forecast was made. Because the
+    plan is the MPC's proposal while Generic Thermostat controls the relay, a forecast
+    is only scored when the relay actually followed the plan up to that horizon;
+    otherwise it is counted as counterfactual. Pending forecasts live in memory;
+    the statistics are saved.
+    """
+
+    def __init__(self, saved: Any = None) -> None:
+        self._pending: deque[dict[str, Any]] = deque(maxlen=60)
+        # (interval end, seconds the relay was on or None if unknown, interval seconds)
+        self._relay: deque[tuple[float, float | None, float]] = deque(maxlen=48)
+        # UTC day -> horizon -> zone -> [count, sum error, sum |error|, sum error^2, sum |persistence error|]
+        self._days: dict[str, dict[str, dict[str, list[float]]]] = {}
+        # UTC day -> horizon -> verdict -> count of forecasts not scored, and why
+        self._verdicts: dict[str, dict[str, dict[str, int]]] = {}
+        self.load(saved)
+
+    def record(self, time: float, start: dict[str, float], trajectory: dict[str, list[float]],
+               plan: list[int], revision: str, assumptions: dict[str, Any]) -> None:
+        self._pending.append({"time": time, "start": dict(start), "trajectory": dict(trajectory),
+                              "plan": list(plan), "revision": revision, "assumptions": assumptions,
+                              "done": set()})
+
+    def observe(self, now: float, actual: dict[str, float | None], on_seconds: float | None,
+                seconds: float | None) -> None:
+        """One poll: the relay's last interval, then score forecasts that are due."""
+        if seconds:
+            self._relay.append((now, on_seconds, seconds))
+        day = datetime.fromtimestamp(now, timezone.utc).date().isoformat()
+        for forecast in list(self._pending):
+            for horizon in HORIZONS:
+                index = horizon // 15 - 1
+                due = forecast["time"] + horizon * 60
+                if horizon in forecast["done"] or now < due - DUE_TOLERANCE_SECONDS:
+                    continue
+                forecast["done"].add(horizon)
+                if index >= len(forecast["plan"]):
+                    continue  # beyond this forecast's horizon
+                if now > due + DUE_TOLERANCE_SECONDS:
+                    self._count(day, horizon, "missing")  # no poll near the due time
+                    continue
+                verdict = self._verify(forecast, horizon)
+                if verdict:
+                    self._count(day, horizon, verdict)
+                    continue
+                for zone, path in forecast["trajectory"].items():
+                    measured, start = actual.get(zone), forecast["start"].get(zone)
+                    if measured is None or start is None or len(path) <= index:
+                        self._count(day, horizon, "missing")
+                        continue
+                    error = path[index] - measured
+                    bucket = self._days.setdefault(day, {}).setdefault(str(horizon), {}).setdefault(
+                        zone, [0, 0.0, 0.0, 0.0, 0.0])
+                    bucket[0] += 1
+                    bucket[1] += error
+                    bucket[2] += abs(error)
+                    bucket[3] += error * error
+                    bucket[4] += abs(start - measured)
+            if len(forecast["done"]) == len(HORIZONS) or now > forecast["time"] + (HORIZONS[-1] + 5) * 60:
+                self._pending.remove(forecast)
+        oldest = (datetime.fromtimestamp(now, timezone.utc).date() - timedelta(days=WINDOW_DAYS - 1)).isoformat()
+        for store in (self._days, self._verdicts):
+            for old in [d for d in store if d < oldest]:
+                del store[old]
+
+    def _verify(self, forecast: dict[str, Any], horizon: int) -> str | None:
+        """None when the relay followed the plan up to `horizon`, else the reason."""
+        for block in range(horizon // 15):
+            start = forecast["time"] + block * 900
+            end = start + 900
+            on = covered = 0.0
+            for interval_end, on_seconds, seconds in self._relay:
+                overlap = min(end, interval_end) - max(start, interval_end - seconds)
+                if overlap <= 0:
+                    continue
+                if on_seconds is None:
+                    return "unverifiable"
+                on += on_seconds * overlap / seconds
+                covered += overlap
+            if covered < 0.8 * 900:
+                return "unverifiable"
+            if abs(on / covered - forecast["plan"][block]) > PLAN_TOLERANCE:
+                return "counterfactual"
+        return None
+
+    def _count(self, day: str, horizon: int, verdict: str) -> None:
+        counts = self._verdicts.setdefault(day, {}).setdefault(str(horizon), {})
+        counts[verdict] = counts.get(verdict, 0) + 1
+
+    def forget(self, zone: str) -> None:
+        for by_horizon in self._days.values():
+            for by_zone in by_horizon.values():
+                by_zone.pop(zone, None)
+        for forecast in self._pending:
+            forecast["trajectory"].pop(zone, None)
+
+    @property
+    def pending(self) -> int:
+        return len(self._pending)
+
+    def summary(self, zones: list[str]) -> dict[int, dict[str, Any]]:
+        result = {}
+        for horizon in HORIZONS:
+            per_zone: dict[str, list[float]] = {}
+            verdicts = {verdict: 0 for verdict in VERDICTS}
+            for by_horizon in self._days.values():
+                for zone, values in by_horizon.get(str(horizon), {}).items():
+                    if zone in zones:
+                        totals = per_zone.setdefault(zone, [0, 0.0, 0.0, 0.0, 0.0])
+                        for i, value in enumerate(values):
+                            totals[i] += value
+            for by_horizon in self._verdicts.values():
+                for verdict, count in by_horizon.get(str(horizon), {}).items():
+                    verdicts[verdict] = verdicts.get(verdict, 0) + count
+            overall = [sum(values[i] for values in per_zone.values()) for i in range(5)]
+            result[horizon] = {**_operational_stats(overall, horizon, len(per_zone)), **verdicts,
+                               "zones": {zone: _operational_stats(values, horizon)
+                                         for zone, values in per_zone.items()}}
+        return result
+
+    def as_dict(self) -> dict[str, Any]:
+        return {"days": {day: {h: {zone: [v[0], *(round(x, 6) for x in v[1:])] for zone, v in by_zone.items()}
+                               for h, by_zone in by_horizon.items()} for day, by_horizon in self._days.items()},
+                "verdicts": self._verdicts}
+
+    def load(self, data: Any) -> None:
+        self._days, self._verdicts = {}, {}
+        if not isinstance(data, dict):
+            return
+        horizons = {str(h) for h in HORIZONS}
+        days = data.get("days") if isinstance(data.get("days"), dict) else {}
+        for day, by_horizon in days.items():
+            if not _valid_day(day) or not isinstance(by_horizon, dict):
+                continue
+            for horizon, by_zone in by_horizon.items():
+                if horizon not in horizons or not isinstance(by_zone, dict):
+                    continue
+                for zone, values in by_zone.items():
+                    if (isinstance(zone, str) and isinstance(values, list) and len(values) == 5
+                            and type(values[0]) is int and values[0] >= 0
+                            and all(type(v) in (int, float) and isfinite(v) for v in values[1:])
+                            and min(values[2:]) >= 0):
+                        self._days.setdefault(day, {}).setdefault(horizon, {})[zone] = [
+                            values[0], *(float(v) for v in values[1:])]
+        verdicts = data.get("verdicts") if isinstance(data.get("verdicts"), dict) else {}
+        for day, by_horizon in verdicts.items():
+            if not _valid_day(day) or not isinstance(by_horizon, dict):
+                continue
+            for horizon, counts in by_horizon.items():
+                if horizon in horizons and isinstance(counts, dict):
+                    clean = {k: v for k, v in counts.items() if k in VERDICTS and type(v) is int and v >= 0}
+                    if clean:
+                        self._verdicts.setdefault(day, {})[horizon] = clean
+
+
+def _valid_day(day: Any) -> bool:
+    try:
+        datetime.strptime(day, "%Y-%m-%d")
+        return True
+    except (TypeError, ValueError):
+        return False
+
+
+def _operational_stats(values: list[float], horizon: int, zones: int = 1) -> dict[str, Any]:
+    count, error, absolute, squares, persistence = values
+    effective = effective_samples(int(count), horizon, zones)
+    if not count:
+        return {"samples": 0, "effective_samples": 0.0, "quality": "insufficient", "mae": None, "rmse": None,
+                "bias": None, "persistence_mae": None, "improvement_vs_persistence": None}
+    return {
+        "samples": int(count), "effective_samples": effective, "quality": quality(effective),
+        "mae": round(absolute / count, 3), "rmse": round(sqrt(squares / count), 3),
+        "bias": round(error / count, 3), "persistence_mae": round(persistence / count, 3),
+        "improvement_vs_persistence": round(1 - absolute / persistence, 3) if persistence > 0 else None,
+    }

custom_components/heating_mpc/sensor.py

@@ -16,6 +16,7 @@ async def async_setup_entry(hass, entry, async_add_entities) -> None:
     entities = [MPCStatus(coordinator, entry.entry_id), MPCWeightedTemp(coordinator, entry.entry_id), MPCOutdoorTemperature(coordinator, entry.entry_id),
                 MPCLearningProgress(coordinator, entry.entry_id), MPCLastSample(coordinator, entry.entry_id)]
     entities += [MPCForecastError(coordinator, entry.entry_id, horizon) for horizon in HORIZONS]
+    entities.append(MPCOperationalForecastError(coordinator, entry.entry_id))
     # Zone forecasts belong to their zone subentry, so deleting a zone removes them.
     linked = {}
     for zone in coordinator.zones.values():
@@ -52,7 +53,7 @@ class MPCStatus(CoordinatorEntity, SensorEntity):
         # Forecast details live on their own entities.
         return {key: value for key, value in data.items() if key not in (
             "status", "prediction_60m", "prediction_end", "prediction_trajectory", "forecast_accuracy",
-            "zone_details")}
+            "zone_details", "operational_forecast")}
 
 
 class MPCWeightedTemp(CoordinatorEntity, SensorEntity):
@@ -133,7 +134,11 @@ class MPCOutdoorTemperature(CoordinatorEntity, SensorEntity):
 
 
 class MPCForecastError(CoordinatorEntity, SensorEntity):
-    """Mean absolute error of the thermal model at one horizon, replayed against reality."""
+    """Model accuracy at one horizon: the model replayed with the inputs that actually happened.
+
+    Not a forecast accuracy (a real forecast does not know the future relay, power or
+    windows); see MPCOperationalForecastError. The unique ID keeps its original name.
+    """
     _attr_has_entity_name = True
     _attr_device_class = SensorDeviceClass.TEMPERATURE_DELTA
     _attr_native_unit_of_measurement = UnitOfTemperature.CELSIUS
@@ -146,7 +151,7 @@ class MPCForecastError(CoordinatorEntity, SensorEntity):
         super().__init__(coordinator)
         self.horizon = horizon
         self._attr_unique_id = f"{entry_id}_forecast_error_{horizon}min"
-        self._attr_name = f"Heating MPC Forecast Error {horizon} min"
+        self._attr_name = f"Heating MPC Model Error {horizon} min"
 
     @property
     def _accuracy(self):
@@ -154,13 +159,20 @@ class MPCForecastError(CoordinatorEntity, SensorEntity):
 
     @property
     def native_value(self):
-        return self._accuracy.get("mae")
+        accuracy = self._accuracy
+        # No number from a handful of overlapping forecasts.
+        return accuracy.get("mae") if accuracy.get("quality", "insufficient") != "insufficient" else None
 
     @property
     def extra_state_attributes(self):
         accuracy = self._accuracy
         return {
+            "kind": "model accuracy: replayed with the actual relay, electrical heat and window history",
             "samples": accuracy.get("samples", 0),
+            "effective_samples": accuracy.get("effective_samples", 0),
+            "quality": accuracy.get("quality", "insufficient"),
+            "mae": accuracy.get("mae"),
+            "rmse": accuracy.get("rmse"),
             "bias": accuracy.get("bias"),
             "persistence_mae": accuracy.get("persistence_mae"),
             "skill": accuracy.get("skill"),
@@ -176,6 +188,47 @@ class MPCForecastError(CoordinatorEntity, SensorEntity):
         }
 
 
+class MPCOperationalForecastError(CoordinatorEntity, SensorEntity):
+    """Accuracy of forecasts as published, using only what was known when they were made.
+
+    A forecast assumes the MPC's plan; it is scored only when the relay actually
+    followed that plan, otherwise counted as counterfactual. State: the 60-minute
+    error; all horizons in the attributes.
+    """
+    _attr_has_entity_name = True
+    _attr_name = "Heating MPC Forecast Error"
+    _attr_device_class = SensorDeviceClass.TEMPERATURE_DELTA
+    _attr_native_unit_of_measurement = UnitOfTemperature.CELSIUS
+    _attr_state_class = SensorStateClass.MEASUREMENT
+    _attr_suggested_display_precision = 2
+    _attr_icon = "mdi:crystal-ball"
+    _unrecorded_attributes = frozenset({"horizons"})
+
+    def __init__(self, coordinator, entry_id):
+        super().__init__(coordinator)
+        self._attr_unique_id = f"{entry_id}_operational_forecast_error"
+
+    @property
+    def native_value(self):
+        hour = self.coordinator.data.get("operational_forecast", {}).get(60, {})
+        return hour.get("mae") if hour.get("quality", "insufficient") != "insufficient" else None
+
+    @property
+    def extra_state_attributes(self):
+        summary = self.coordinator.data.get("operational_forecast", {})
+        names = {key: zone.name for key, zone in self.coordinator.zones.items()}
+        horizons = {}
+        for horizon, stats in summary.items():
+            horizons[f"{horizon}_min"] = {**{k: v for k, v in stats.items() if k != "zones"},
+                                          "zones": {names.get(z, z): v for z, v in stats.get("zones", {}).items()}}
+        return {
+            "kind": "forecast accuracy: published forecasts, scored only where the relay followed the plan",
+            "horizons": horizons,
+            "window_days": WINDOW_DAYS,
+            "unit_of_attributes": "°C",
+        }
+
+
 class MPCLearningProgress(CoordinatorEntity, SensorEntity):
     """How close the least trained room is to being used for forecasts."""
     _attr_has_entity_name = True

tests/ha/test_lifecycle.py

@@ -39,7 +39,7 @@ async def test_setup_reload_unload(hass, entry, inputs):
     assert coordinator.config_entry is entry
     registry = er.async_get(hass)
     entities = er.async_entries_for_config_entry(registry, entry.entry_id)
-    assert len(entities) == 11  # 7 core + 4 forecast-error sensors
+    assert len(entities) == 12  # 7 core + 4 model-error sensors + 1 forecast-error sensor
     ids = {e.unique_id: e.entity_id for e in entities}
     assert coordinator.data["status"] == "learning"
     status = hass.states.get(ids[f"{entry.entry_id}_status"])
@@ -74,7 +74,7 @@ async def test_replacing_a_zone_removes_only_its_forecast(hass, entry, inputs):
     assert registry.async_get(old) is None
     new = registry.async_get_entity_id("sensor", DOMAIN, f"{entry.entry_id}_zone_new_predicted_60min")
     assert hass.states.get(new).attributes["friendly_name"] == "MPC New 60min Prediction"
-    assert len(er.async_entries_for_config_entry(registry, entry.entry_id)) == 11
+    assert len(er.async_entries_for_config_entry(registry, entry.entry_id)) == 12
 
 
 async def test_storage_restore_and_remove(hass, entry, inputs):
@@ -255,13 +255,36 @@ async def test_forecast_error_sensors_score_scheduled_polls(hass, entry, inputs,
         async_fire_time_changed(hass, future)
         await hass.async_block_till_done()
     state = hass.states.get(ids[15])
-    assert coordinator.data["forecast_accuracy"][15]["samples"] == 1
-    assert float(state.state) == coordinator.data["forecast_accuracy"][15]["mae"] > 0
+    accuracy = coordinator.data["forecast_accuracy"][15]
+    assert accuracy["samples"] == 1
+    # One replayed forecast is no accuracy figure: the state stays unknown, the numbers
+    # are in the attributes with their quality.
+    assert state.state == "unknown"
+    assert state.attributes["mae"] == accuracy["mae"] > 0
+    assert state.attributes["quality"] == "insufficient"
     assert state.attributes["samples"] == 1
     assert state.attributes["zones"]["room"]["samples"] == 1  # by zone name
     assert state.attributes["unit_of_measurement"] == "°C"
+    assert state.attributes["friendly_name"] == "Heating MPC Model Error 15 min"
     assert "zones" in state.state_info["unrecorded_attributes"]
     assert hass.states.get(ids[30]).state == "unknown"
+    for step in range(15):  # enough 15-minute spans for an indicative figure
+        future += timedelta(minutes=5, seconds=1)
+        freezer.move_to(future)
+        # Fresh readings, or the sensors turn stale after 30 minutes.
+        hass.states.async_set("sensor.room", str(20 - 0.01 * step), {"unit_of_measurement": "°C"})
+        hass.states.async_set("sensor.outside", "5", {"unit_of_measurement": "°C"}, force_update=True)
+        async_fire_time_changed(hass, future)
+        await hass.async_block_till_done()
+    accuracy = coordinator.data["forecast_accuracy"][15]
+    assert accuracy["quality"] == "indicative", accuracy
+    assert float(hass.states.get(ids[15]).state) == accuracy["mae"]
+    # Published forecasts are scored separately; nothing was published while learning.
+    operational = registry.async_get_entity_id("sensor", DOMAIN, f"{entry.entry_id}_operational_forecast_error")
+    state = hass.states.get(operational)
+    assert state.state == "unknown"
+    assert state.attributes["friendly_name"] == "Heating MPC Forecast Error"
+    assert "horizons" in state.state_info["unrecorded_attributes"]
 
 
 async def test_learning_progress_and_last_sample_sensors(hass, entry, inputs, freezer):

tests/test_coordinator.py

@@ -1117,8 +1117,8 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
         vent.count = 80
         vent.theta[vent.features.index('vent_natural')] = 0.01
         vent.theta[vent.features.index('vent_wind_dir')] = 0.0005
-        c._evaluator._days = {'2026-10-27': {'60': {'living': [10, 0.1, 0.2, 0.3, 0.4],
-                                                    'bed': [10, 0.1, 0.2, 0.3, 0.4]}}}
+        c._evaluator._days = {'2026-10-27': {'60': {'living': [10, 0.1, 0.2, 0.3, 0.4, 0.05, 10],
+                                                    'bed': [10, 0.1, 0.2, 0.3, 0.4, 0.05, 10]}}}
         c._heat_sources._baselines['ups'].samples.extend([200.0] * 24)
         return c._payload('sensor.outside')
 
@@ -1261,6 +1261,35 @@ class BoundaryTests(unittest.IsolatedAsyncioTestCase):
         self.assertEqual(newer['status'], 'ready')
         self.assertIs(result, newer)
 
+    async def test_published_forecasts_are_scored_only_when_the_relay_followed_the_plan(self):
+        # The relay stays off. A warm room's plan says off, so its forecast is scored;
+        # a cold room's plan says heat, so its forecast is counterfactual, never an error.
+        for room, follows in (('23.5', True), ('18', False)):
+            with self.subTest(room=room):
+                self.setUp()
+                self.ready_office()
+                self.set_state('sensor.room', room, unit_of_measurement='°C')
+                first = await self.sample()
+                self.assertEqual(first['status'], 'ready')
+                self.assertEqual(first['plan_15min'][0], 0 if follows else 1)
+                for _ in range(3):
+                    self.fresh('sensor.room', 'sensor.outside', 'sensor.ups_power')
+                    data = await self.sample()
+                fifteen = data['operational_forecast'][15]
+                if follows:
+                    self.assertEqual((fifteen['samples'], fifteen['counterfactual']), (1, 0))
+                    predicted = first['prediction_trajectory']['office'][0]
+                    self.assertAlmostEqual(fifteen['mae'], abs(predicted - 23.5), places=3)
+                    self.assertEqual(fifteen['persistence_mae'], 0.0)
+                    self.assertEqual(fifteen['quality'], 'insufficient')
+                else:
+                    self.assertEqual((fifteen['samples'], fifteen['counterfactual'], fifteen['mae']), (0, 1, None))
+                # Longer horizons are not due yet.
+                self.assertEqual(data['operational_forecast'][30]['samples'], 0)
+                stored = self.c._payload('sensor.outside')['operational_evaluation']
+                self.assertEqual(sum(sum(v.get('15', {}).values()) for v in stored['verdicts'].values()),
+                                 0 if follows else 1)
+
     async def test_unloading_during_save_or_optimization_publishes_nothing_new(self):
         self.ready_office()
         await self.sample()                       # a first sample to learn from

tests/test_evaluation.py

@@ -69,6 +69,16 @@ class EvaluationTests(unittest.TestCase):
         self.assertEqual(uncertain.summary(['sensor.room'])[15]['samples'], 2)  # k = 3 and 7 only
         self.assertEqual(plant(8).summary(['sensor.other'])[15]['samples'], 0)
 
+    def test_rmse_effective_samples_and_quality(self):
+        stats = plant(10, model_theta=(0.004, 0.0, 0.0005), lag_at=lambda k: 1.0).summary(['sensor.room'])
+        s15 = stats[15]
+        self.assertGreaterEqual(s15['rmse'], s15['mae'])
+        # 7 overlapping 15-minute forecasts started 5 minutes apart: ~2.3 independent ones.
+        self.assertEqual((s15['samples'], s15['effective_samples'], s15['quality']), (7, 2.3, 'insufficient'))
+        self.assertEqual(ev.quality(5), 'indicative')
+        self.assertEqual(ev.quality(30), 'ok')
+        self.assertEqual(ev.effective_samples(120, 60, zones=2), 5.0)
+
     def test_window_keeps_seven_days(self):
         evaluator = plant(10)
         before = evaluator.summary(['sensor.room'])[15]['samples']
@@ -108,12 +118,91 @@ class EvaluationTests(unittest.TestCase):
                          'sensor.c': [1, float('nan'), 0, 0], 'sensor.d': [2, 0.1, 0.2, 0.4],
                          'sensor.e': [2, 0.1, 0.2, 0.4, 0.3], 'sensor.f': [2, 0.1, 0.2, 0.4, -1]}}}})
         # A 0.3 bucket has no base-model column; its model was the base model then.
+        # Neither 0.3 nor 0.4 buckets have squared errors: they do not count for RMSE.
         self.assertEqual(restored.as_dict(), {'days': {day: {'15': {
-            'sensor.d': [2, 0.1, 0.2, 0.4, 0.2], 'sensor.e': [2, 0.1, 0.2, 0.4, 0.3]}}}})
+            'sensor.d': [2, 0.1, 0.2, 0.4, 0.2, 0.0, 0], 'sensor.e': [2, 0.1, 0.2, 0.4, 0.3, 0.0, 0]}}}})
+        self.assertIsNone(restored.summary(['sensor.d'])[15]['rmse'])
         for garbage in (None, [], {'days': []}):
             restored.load(garbage)
             self.assertEqual(restored.as_dict(), {'days': {}})
 
 
+T0 = START.timestamp()
+
+
+def run_ledger(plan, relay_on, temps, trajectory, start=20.0, polls=25, unknown_relay_at=None, skip=()):
+    """One forecast at T0, then a poll every 5 minutes with the given relay state per poll."""
+    ledger = ev.OperationalForecastLedger()
+    ledger.observe(T0, {'r': start}, None, None)
+    ledger.record(T0, {'r': start}, {'r': trajectory}, plan, 'rev1', {'strategy': 'mpc_plan'})
+    for k in range(1, polls):
+        if k in skip:
+            continue
+        on = None if k == unknown_relay_at else (300.0 if relay_on(k) else 0.0)
+        ledger.observe(T0 + 300 * k, {'r': temps(k)}, on, 300.0)
+    return ledger.summary(['r'])
+
+
+class OperationalForecastTests(unittest.TestCase):
+    PLAN_OFF = [0] * 8
+    TRAJECTORY = [19.9, 19.8, 19.7, 19.6, 19.5, 19.4, 19.3, 19.2]  # every 15 minutes
+
+    def test_followed_plan_is_scored_at_every_horizon_against_persistence(self):
+        stats = run_ledger(self.PLAN_OFF, lambda k: False, lambda k: 20.0 - k / 30, self.TRAJECTORY)
+        for horizon, actual in ((15, 19.9), (30, 19.8), (60, 19.6), (120, 19.2)):
+            with self.subTest(horizon=horizon):
+                s = stats[horizon]
+                self.assertEqual((s['samples'], s['counterfactual'], s['unverifiable']), (1, 0, 0))
+                self.assertAlmostEqual(s['mae'], 0, places=3)
+                self.assertAlmostEqual(s['persistence_mae'], 20.0 - actual, places=3)
+                self.assertEqual(s['improvement_vs_persistence'], 1.0)
+                self.assertEqual(s['quality'], 'insufficient')  # one forecast proves little
+
+    def test_plan_not_followed_is_counterfactual_not_scored(self):
+        stats = run_ledger([1] * 8, lambda k: False, lambda k: 20.0, self.TRAJECTORY)
+        for horizon in (15, 30, 60, 120):
+            self.assertEqual((stats[horizon]['samples'], stats[horizon]['counterfactual']), (0, 1))
+            self.assertIsNone(stats[horizon]['mae'])
+
+    def test_plan_followed_for_a_while(self):
+        # The thermostat switches the boiler on in the 4th block (after 45 minutes).
+        stats = run_ledger(self.PLAN_OFF, lambda k: k > 9, lambda k: 20.0, self.TRAJECTORY)
+        self.assertEqual([stats[h]['samples'] for h in (15, 30, 60, 120)], [1, 1, 0, 0])
+        self.assertEqual([stats[h]['counterfactual'] for h in (15, 30, 60, 120)], [0, 0, 1, 1])
+
+    def test_unknown_relay_history_or_missing_observations(self):
+        stats = run_ledger(self.PLAN_OFF, lambda k: False, lambda k: 20.0, self.TRAJECTORY, unknown_relay_at=2)
+        self.assertEqual((stats[15]['samples'], stats[15]['unverifiable'], stats[60]['unverifiable']), (0, 1, 1))
+        missing = run_ledger(self.PLAN_OFF, lambda k: False, lambda k: None if k == 12 else 20.0, self.TRAJECTORY)
+        self.assertEqual((missing[60]['samples'], missing[60]['missing'], missing[30]['samples']), (0, 1, 1))
+        # No poll near the due time (polls 11 to 13 missing): the 60-minute point cannot be scored.
+        gap = run_ledger(self.PLAN_OFF, lambda k: False, lambda k: 20.0, self.TRAJECTORY, skip=(11, 12, 13))
+        self.assertEqual((gap[60]['samples'], gap[60]['missing']), (0, 1))
+
+    def test_due_time_alignment(self):
+        ledger = ev.OperationalForecastLedger()
+        ledger.record(T0, {'r': 20.0}, {'r': self.TRAJECTORY}, self.PLAN_OFF, 'rev1', {})
+        for offset in (300, 600):
+            ledger.observe(T0 + offset, {'r': 20.0}, 0.0, 300.0)
+        ledger.observe(T0 + 900 - 200, {'r': 20.0}, 0.0, 100.0)  # too early for +15 min
+        self.assertEqual(ledger.summary(['r'])[15]['samples'], 0)
+        ledger.observe(T0 + 900 + 100, {'r': 19.9}, 0.0, 300.0)  # within half a poll of the due time
+        self.assertEqual(ledger.summary(['r'])[15]['samples'], 1)
+
+    def test_round_trip_forget_and_garbage(self):
+        stats_ledger = ev.OperationalForecastLedger()
+        stats_ledger.record(T0, {'r': 20.0}, {'r': self.TRAJECTORY}, [1] * 8, 'rev1', {})
+        stats_ledger.record(T0, {'r': 20.0}, {'r': self.TRAJECTORY}, self.PLAN_OFF, 'rev1', {})
+        for k in range(1, 4):
+            stats_ledger.observe(T0 + 300 * k, {'r': 19.9}, 0.0, 300.0)
+        restored = ev.OperationalForecastLedger(stats_ledger.as_dict())
+        self.assertEqual(restored.summary(['r']), stats_ledger.summary(['r']))
+        self.assertEqual(restored.summary(['r'])[15]['counterfactual'], 1)
+        restored.forget('r')
+        self.assertEqual(restored.summary(['r'])[15]['samples'], 0)
+        for garbage in (None, [], {'days': {'x': {}}, 'verdicts': {'2026-10-27': {'15': {'bogus': 3}}}}):
+            self.assertEqual(ev.OperationalForecastLedger(garbage).as_dict(), {'days': {}, 'verdicts': {}})
+
+
 if __name__ == '__main__':
     unittest.main()