Files
heating-mpc-ha/custom_components/heating_mpc/evaluation.py
T
buzeraandClaude Opus 5.5 26c415fa6c 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>
2026-10-27 03:47:16 +02:00

467 lines
23 KiB
Python

"""Walk-forward evaluation of the thermal model against what actually happened.
Pure Python; no Home Assistant imports.
At every poll the readings, the outdoor temperature, the radiator-lag estimate
and each room's model parameters are recorded. When 15/30/60/120 minutes have
passed, each room is replayed from that earlier poll with the parameters it had
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.
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, timezone
from math import isfinite, sqrt
from typing import Any
from .physics import Learner
from .windows import ventilation_vector_for_prediction
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],
steps: list[tuple[float, float]]) -> float:
"""Advance one room through (end-of-step lag, step seconds) with fixed parameters."""
model = Learner(theta=list(theta))
for lag, seconds in steps:
# Same 5-minute model and jitter normalization as learning.
temp += model.predict_delta(temp, outside, lag) * seconds / STEP_SECONDS
return temp
def _dot(theta: tuple[float, ...], x: list[float]) -> float:
return sum(t * v for t, v in zip(theta, x))
def model_snapshot(model: Any) -> dict[str, Any]:
"""Frozen parameters of a ZoneModel, so later learning never leaks into a replay."""
def layer(learner):
return (list(learner.features), tuple(learner.theta)) if learner is not None else None
return {"base": tuple(model.base.theta), "uses_extended": model.uses_extended,
"extended": layer(model.extended), "ventilation": layer(model.ventilation)}
def snapshot_delta(snapshot: dict[str, Any], temp: float, inputs: dict[str, Any]) -> float:
"""5-minute change predicted by a model snapshot (see zone_model.ZoneModel.delta)."""
main = None
if snapshot["uses_extended"] and snapshot["extended"]:
features, theta = snapshot["extended"]
x = feature_vector(features, temp, inputs)
main = _dot(theta, x) if x is not None else None
if main is None:
a, b, c = snapshot["base"]
main = a * (inputs["outside"] - temp) + b * inputs["lag"] + c
ventilation = inputs.get("ventilation")
if snapshot["ventilation"] and ventilation:
features, theta = snapshot["ventilation"]
x = ventilation_vector_for_prediction(features, temp, inputs["outside"], ventilation)
if x is not None:
main += _dot(theta, x)
return main
def replay_zones(origin: dict[str, Any], path: list[dict[str, Any]]) -> dict[str, float]:
"""All zones advanced together with the origin's models and the actual inputs."""
state = dict(origin["temps"])
for before, poll in zip(path, path[1:]):
seconds = (poll["time"] - before["time"]).total_seconds()
new = {}
for zone, temp in state.items():
actual = poll["inputs"].get(zone, {})
inputs = {"outside": origin["outside"], "lag": poll["lag"], "temps": state,
"internal": actual.get("internal"), "ventilation": actual.get("ventilation")}
new[zone] = temp + snapshot_delta(origin["models"][zone], temp, inputs) * seconds / STEP_SECONDS
state = new
return state
class ForecastEvaluator:
"""Score replayed forecasts per room and horizon over the last WINDOW_DAYS days."""
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|, 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:
"""Inputs were missing: no forecast may be scored across this poll."""
self._history.clear()
def observe(self, time: datetime, temps: dict[str, float], outside: float, lag: float,
lag_uncertainty: float, theta: dict[str, tuple[float, float, float]],
models: dict[str, dict[str, Any]] | None = None,
inputs: dict[str, dict[str, Any]] | None = None) -> None:
"""Record a poll, then score every forecast that has reached its horizon.
`models`: model_snapshot per zone (base-only when omitted); `inputs`: this
poll's actual electrical heat and ventilation inputs per zone.
"""
if self._history and not 270 <= (time - self._history[-1]["time"]).total_seconds() <= 330:
self._history.clear() # irregular poll spacing: replays would be unreliable
models = models or {}
snapshots = {zone: models.get(zone) or {"base": tuple(theta[zone]), "uses_extended": False,
"extended": None, "ventilation": None}
for zone in temps if zone in theta}
self._history.append({"time": time, "temps": dict(temps), "outside": outside, "lag": lag,
"lag_uncertainty": lag_uncertainty, "theta": dict(theta),
"models": snapshots, "inputs": dict(inputs or {})})
self._score()
self._prune(time)
def _score(self) -> None:
history = list(self._history)
latest = history[-1]
day = latest["time"].date().isoformat()
for horizon in HORIZONS:
steps = horizon * 60 // STEP_SECONDS
if len(history) < steps + 1:
continue
path = history[-(steps + 1):]
# The actual heat input must be known along the way.
if any(poll["lag_uncertainty"] > self._max_lag_uncertainty for poll in path[1:]):
continue
origin = path[0]
inputs = [(poll["lag"], (poll["time"] - before["time"]).total_seconds())
for before, poll in zip(path, path[1:])]
full = replay_zones(origin, path)
for room, actual in latest["temps"].items():
start, theta = origin["temps"].get(room), origin["theta"].get(room)
if start is None or theta is None:
continue
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, 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()
for day in [day for day in self._days if day < oldest]:
del self._days[day]
def summary(self, rooms: list[str]) -> dict[int, dict[str, Any]]:
"""Error statistics per horizon, overall and per configured room."""
result = {}
for horizon in HORIZONS:
per_room: dict[str, list[float]] = {}
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, 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(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:
"""Drop a zone's statistics, for example when its temperature is measured differently."""
for by_horizon in self._days.values():
for by_room in by_horizon.values():
by_room.pop(room, None)
for poll in self._history:
poll["temps"].pop(room, None)
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: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()}}
def load(self, data: Any) -> None:
"""Restore stored statistics, keeping only well-formed entries."""
self._days = {}
days = data.get("days") if isinstance(data, dict) else None
if not isinstance(days, dict):
return
horizons = {str(horizon) for horizon in HORIZONS}
for day, by_horizon in days.items():
try:
datetime.strptime(day, "%Y-%m-%d")
except (TypeError, ValueError):
continue
if not isinstance(by_horizon, dict):
continue
for horizon, by_room in by_horizon.items():
if horizon not in horizons or not isinstance(by_room, dict):
continue
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.
# 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),
*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, 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
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], 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, "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).
"skill": round(1 - absolute / persistence, 3) if persistence > 0 else None,
"base_model_mae": round(base / count, 3),
# Share of the base model's error removed by electrical heat, neighbours and
# 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,
}