diff --git a/custom_components/heating_mpc/coordinator.py b/custom_components/heating_mpc/coordinator.py index 359f447..6ca878e 100644 --- a/custom_components/heating_mpc/coordinator.py +++ b/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 diff --git a/custom_components/heating_mpc/evaluation.py b/custom_components/heating_mpc/evaluation.py index 16fffa3..2d2e89e 100644 --- a/custom_components/heating_mpc/evaluation.py +++ b/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, + } diff --git a/custom_components/heating_mpc/sensor.py b/custom_components/heating_mpc/sensor.py index 321bce3..5a83a93 100644 --- a/custom_components/heating_mpc/sensor.py +++ b/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 diff --git a/tests/ha/test_lifecycle.py b/tests/ha/test_lifecycle.py index d4edf6f..01fa9f7 100644 --- a/tests/ha/test_lifecycle.py +++ b/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): diff --git a/tests/test_coordinator.py b/tests/test_coordinator.py index 7b2525b..5886912 100644 --- a/tests/test_coordinator.py +++ b/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 diff --git a/tests/test_evaluation.py b/tests/test_evaluation.py index 0514748..2a1ede4 100644 --- a/tests/test_evaluation.py +++ b/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()