Reproduced: a window unknown for a whole 15-minute span was replayed as closed and scored as a Model Error sample (MAE 0.76 °C on one sample). - A zone's replay is complete only if every step's model inputs were observed: window openings (and a confirmed recovery), the wind, direction and other zones' windows that an open window's terms use, and, for the extended layer, electrical heat and neighbours, whose own replay must be complete too. - Incomplete spans are counted, not scored; the base model's error on them is kept separately (it needs only the relay and outdoor temperature). Spans that could not be replayed (uncertain radiator lag, no reading at the end) are counted as skipped. Model Error sensors show incomplete, skipped and incomplete_base_model_mae; stored data gains an "excluded" section and older data loads unchanged. - Window inputs no longer stand in for unknown values: the facing term is unknown while a window is open and the wind direction is not, and the cross-ventilation input is unknown while another zone's window is. With every window closed, opening terms are 0 and need no weather input. Learning refuses such samples; forecasts still assume calm air. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
302 lines
18 KiB
Python
302 lines
18 KiB
Python
"""Walk-forward forecast evaluation, without Home Assistant."""
|
|
import sys
|
|
import types
|
|
import unittest
|
|
from datetime import datetime, timedelta, timezone
|
|
from pathlib import Path
|
|
|
|
COMPONENTS = Path(__file__).resolve().parents[1] / 'custom_components' / 'heating_mpc'
|
|
# Fake only the package container, to avoid importing Home Assistant from __init__.
|
|
sys.modules.setdefault('heating_mpc', types.ModuleType('heating_mpc'))
|
|
sys.modules['heating_mpc'].__path__ = [str(COMPONENTS)]
|
|
from heating_mpc import evaluation as ev # noqa: E402
|
|
|
|
THETA = (0.004, 0.12, 0.0005)
|
|
START = datetime(2026, 10, 26, tzinfo=timezone.utc)
|
|
|
|
|
|
def plant(polls, theta=THETA, model_theta=None, outside=5.0, lag_at=lambda k: 1.0 if k % 30 < 12 else 0.0,
|
|
evaluator=None, uncertainty_at=lambda k: 0.0, start=START, step=300):
|
|
"""Feed polls of a room that follows the 5-minute model exactly with `theta`."""
|
|
a, b, c = theta
|
|
evaluator = evaluator or ev.ForecastEvaluator(0.1)
|
|
temp = 20.0
|
|
for k in range(polls):
|
|
lag = lag_at(k)
|
|
if k:
|
|
temp += a * (outside - temp) + b * lag + c
|
|
evaluator.observe(start + timedelta(seconds=k * step), {'sensor.room': temp}, outside, lag,
|
|
uncertainty_at(k), {'sensor.room': model_theta or theta})
|
|
return evaluator
|
|
|
|
|
|
class EvaluationTests(unittest.TestCase):
|
|
def test_exact_model_scores_zero_error_and_beats_persistence(self):
|
|
summary = plant(60).summary(['sensor.room'])
|
|
for horizon, expected_samples in ((15, 57), (30, 54), (60, 48), (120, 36)):
|
|
with self.subTest(horizon=horizon):
|
|
stats = summary[horizon]
|
|
self.assertEqual(stats['samples'], expected_samples)
|
|
self.assertEqual(stats['rooms']['sensor.room']['samples'], expected_samples)
|
|
self.assertAlmostEqual(stats['mae'], 0, places=6)
|
|
self.assertGreater(stats['persistence_mae'], 0.05)
|
|
self.assertEqual(stats['skill'], 1)
|
|
|
|
def test_parameters_from_the_forecast_start_are_used(self):
|
|
# A model without heat gain under-predicts warming: negative bias.
|
|
stats = plant(10, model_theta=(0.004, 0.0, 0.0005), lag_at=lambda k: 1.0).summary(['sensor.room'])
|
|
self.assertLess(stats[15]['bias'], -0.3)
|
|
self.assertEqual(stats[15]['bias'], -stats[15]['mae'])
|
|
|
|
def test_replay_matches_hand_calculation(self):
|
|
temp = 20.0
|
|
for _ in range(3):
|
|
temp += 0.004 * (5 - temp) + 0.12 * 0.5 + 0.0005
|
|
self.assertAlmostEqual(ev.replay(20.0, 5.0, THETA, [(0.5, 300)] * 3), temp)
|
|
# Jitter is normalized like learning: 330 s counts as 1.1 steps.
|
|
self.assertAlmostEqual(ev.replay(20.0, 5.0, (0, 0.1, 0), [(1.0, 330)]), 20.11)
|
|
|
|
def test_breaks_uncertain_lag_and_missing_rooms_are_not_scored(self):
|
|
evaluator = plant(4)
|
|
self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 1)
|
|
evaluator.reset_window()
|
|
plant(3, evaluator=evaluator, start=START + timedelta(minutes=20))
|
|
self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 1)
|
|
# Irregular spacing restarts the window as well.
|
|
jittery = plant(8, step=400)
|
|
self.assertEqual(jittery.summary(['sensor.room'])[15]['samples'], 0)
|
|
uncertain = plant(8, uncertainty_at=lambda k: 0.5 if k == 4 else 0.0)
|
|
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']
|
|
plant(10, evaluator=evaluator, start=START + timedelta(days=7))
|
|
self.assertEqual(evaluator.summary(['sensor.room'])[15]['samples'], 7)
|
|
self.assertEqual(before, 7)
|
|
|
|
def test_full_model_replays_actual_inputs_and_is_compared_with_base(self):
|
|
base = (0.004, 0.12, 0.0005)
|
|
full = {'base': base, 'uses_extended': True,
|
|
'extended': (['loss', 'heat', 'bias', 'internal'], base + (0.15,)),
|
|
'ventilation': (['vent_natural'], (0.01,))}
|
|
evaluator, temp = ev.ForecastEvaluator(0.1), 20.0
|
|
for k in range(40):
|
|
internal = 0.4 if (k // 8) % 2 else 0.0 # gaming on and off
|
|
opening = 1.0 if k % 10 in (3, 4) else 0.0 # short ventilation
|
|
if k:
|
|
temp += 0.004 * (5 - temp) + 0.12 * 0.5 + 0.0005 + 0.15 * internal + 0.01 * opening * (5 - temp)
|
|
evaluator.observe(START + timedelta(seconds=300 * k), {'r': temp}, 5.0, 0.5, 0.0, {'r': base},
|
|
models={'r': full},
|
|
inputs={'r': {'internal': internal, 'ventilation': {'opening': opening}}})
|
|
stats = evaluator.summary(['r'])[60]
|
|
self.assertAlmostEqual(stats['mae'], 0, places=6)
|
|
self.assertGreater(stats['base_model_mae'], 0.05)
|
|
self.assertEqual(stats['improvement_vs_base'], 1.0)
|
|
|
|
def test_saved_statistics_round_trip_and_reject_garbage(self):
|
|
evaluator = plant(10)
|
|
restored = ev.ForecastEvaluator(0.1)
|
|
restored.load(evaluator.as_dict())
|
|
self.assertEqual(restored.summary(['sensor.room']), evaluator.summary(['sensor.room']))
|
|
day = START.date().isoformat()
|
|
restored.load({'days': {
|
|
'not-a-day': {'15': {'sensor.room': [1, 0.0, 0.1, 0.1]}},
|
|
day: {'45': {'sensor.room': [1, 0.0, 0.1, 0.1]},
|
|
'15': {'sensor.a': [1, 0.0, -0.1, 0.1], 'sensor.b': [True, 0, 0, 0],
|
|
'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]}}},
|
|
'excluded': {day: {'15': {'sensor.d': [1, 2, 0.3], 'sensor.x': [1, -2, 0.3], 'sensor.y': [1, 2]},
|
|
'45': {'sensor.d': [1, 2, 0.3]}}}})
|
|
# 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.
|
|
# 0.4 data has no excluded-span counts.
|
|
self.assertEqual(restored.as_dict(), {'days': {day: {'15': {
|
|
'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]}}},
|
|
'excluded': {day: {'15': {'sensor.d': [1, 2, 0.3]}}}})
|
|
self.assertIsNone(restored.summary(['sensor.d'])[15]['rmse'])
|
|
self.assertEqual({k: restored.summary(['sensor.d'])[15][k] for k in ('incomplete', 'skipped',
|
|
'incomplete_base_model_mae')},
|
|
{'incomplete': 1, 'skipped': 2, 'incomplete_base_model_mae': 0.3})
|
|
for garbage in (None, [], {'days': []}):
|
|
restored.load(garbage)
|
|
self.assertEqual(restored.as_dict(), {'days': {}, 'excluded': {}})
|
|
|
|
|
|
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()
|
|
|
|
|
|
BASE = (0.004, 0.12, 0.0005)
|
|
|
|
|
|
def ventilated(observed, actual=None, wind=4.0, features=('vent_natural',), zones=('r',), internal=None,
|
|
extended=None, lag_uncertain_at=(), missing_end_at=()):
|
|
"""Zones cooled by real ventilation; `observed[k]` is what the window sensor said
|
|
for the interval ending at poll k (None: unknown). Returns the summary at 15 min."""
|
|
actual = actual if actual is not None else [1.0 if o is None else o for o in observed]
|
|
snapshot = {'base': BASE, 'uses_extended': extended is not None, 'extended': extended,
|
|
'ventilation': (list(features), (0.02,) * len(features))}
|
|
evaluator, temps = ev.ForecastEvaluator(0.1), {zone: 20.0 for zone in zones}
|
|
for k, seen in enumerate(observed):
|
|
if k:
|
|
temps = {zone: t + 0.004 * (5 - t) + 0.12 * 0.5 + 0.0005
|
|
+ 0.02 * actual[k] * (5 - t) * (1 + (wind or 0) * ('vent_wind' in features))
|
|
for zone, t in temps.items()}
|
|
vent = {'opening': seen, 'facing': 0.0, 'other_opening': 0.0, 'recovery': 0.0, 't_ref': None,
|
|
'wind': wind}
|
|
readings = {z: t for z, t in temps.items() if not (k in missing_end_at and z == zones[0])}
|
|
evaluator.observe(START + timedelta(seconds=300 * k), readings, 5.0, 0.5,
|
|
0.5 if k in lag_uncertain_at else 0.0, {z: BASE for z in zones},
|
|
models={z: snapshot for z in zones},
|
|
inputs={z: {'internal': internal, 'ventilation': vent if z == zones[0] else
|
|
{**vent, 'opening': actual[k]}} for z in zones})
|
|
return evaluator.summary(list(zones))[15]
|
|
|
|
|
|
class IncompleteInputTests(unittest.TestCase):
|
|
"""Unknown ventilation or extended inputs never become Model Error samples."""
|
|
|
|
def test_window_unknown_for_the_whole_span_is_incomplete(self):
|
|
stats = ventilated([0.0, None, None, None])
|
|
self.assertEqual((stats['samples'], stats['mae'], stats['incomplete']), (0, None, 1))
|
|
# The base model needs no window input: its error is still kept, apart.
|
|
self.assertGreater(stats['incomplete_base_model_mae'], 0.1)
|
|
|
|
def test_window_unknown_halfway_or_recovering_during_the_span(self):
|
|
self.assertEqual(ventilated([0.0, 1.0, None, 1.0])['incomplete'], 1) # unavailable halfway
|
|
stats = ventilated([0.0, None, 1.0, 1.0, 1.0, 1.0, 1.0]) # back after one poll
|
|
# Spans ending at polls 3-6: only the first contains the unknown interval.
|
|
self.assertEqual((stats['incomplete'], stats['samples']), (1, 3))
|
|
self.assertAlmostEqual(stats['mae'], 0, places=6) # the known spans are exact
|
|
|
|
def test_one_zone_unknown_leaves_the_other_zone_scored(self):
|
|
stats = ventilated([0.0, None, None, None], zones=('a', 'b'))
|
|
self.assertEqual((stats['rooms']['a']['incomplete'], stats['rooms']['a']['samples']), (1, 0))
|
|
self.assertEqual((stats['rooms']['b']['incomplete'], stats['rooms']['b']['samples']), (0, 1))
|
|
|
|
def test_missing_wind_matters_only_while_a_window_is_open(self):
|
|
features = ('vent_natural', 'vent_wind')
|
|
open_span = ventilated([0.0, 1.0, 1.0, 1.0], wind=None, features=features)
|
|
closed_span = ventilated([0.0, 0.0, 0.0, 0.0], wind=None, features=features)
|
|
self.assertEqual((open_span['incomplete'], open_span['samples']), (1, 0))
|
|
self.assertEqual((closed_span['incomplete'], closed_span['samples']), (0, 1))
|
|
self.assertEqual(ventilated([0.0, 1.0, 1.0, 1.0], features=features)['samples'], 1)
|
|
|
|
def test_unknown_electrical_heat_makes_the_extended_replay_incomplete(self):
|
|
extended = (['loss', 'heat', 'bias', 'internal'], BASE + (0.15,))
|
|
stats = ventilated([0.0] * 4, extended=extended, internal=None)
|
|
self.assertEqual((stats['incomplete'], stats['samples']), (1, 0))
|
|
self.assertIsNotNone(stats['incomplete_base_model_mae'])
|
|
self.assertEqual(ventilated([0.0] * 4, extended=extended, internal=0.0)['samples'], 1)
|
|
|
|
def test_a_neighbour_with_unknown_inputs_spoils_a_coupled_zone(self):
|
|
coupled = {'base': BASE, 'uses_extended': True, 'ventilation': None,
|
|
'extended': (['loss', 'heat', 'bias', 'couple:a'], BASE + (0.01,))}
|
|
plain = {'base': BASE, 'uses_extended': False, 'extended': None,
|
|
'ventilation': (['vent_natural'], (0.02,))}
|
|
origin = {'temps': {'a': 20.0, 'b': 20.0}, 'outside': 5.0, 'models': {'a': plain, 'b': coupled}}
|
|
path = [origin] + [{'time': START + timedelta(seconds=300 * k), 'lag': 0.5,
|
|
'inputs': {'a': {'ventilation': {'opening': None}}, 'b': {}}} for k in range(4)]
|
|
origin['time'] = START - timedelta(seconds=300)
|
|
_, incomplete = ev.replay_zones(origin, path)
|
|
self.assertEqual(incomplete, {'a', 'b'})
|
|
path = [origin] + [{**poll, 'inputs': {'a': {'ventilation': {'opening': 0.0}}, 'b': {}}} for poll in path[1:]]
|
|
self.assertEqual(ev.replay_zones(origin, path)[1], set())
|
|
|
|
def test_skipped_spans_are_counted_and_samples_stay_correct(self):
|
|
stats = ventilated([0.0] * 9, lag_uncertain_at=(4,))
|
|
# Spans ending at 3..8: those crossing poll 4 (ending 4, 5, 6) are skipped.
|
|
self.assertEqual((stats['samples'], stats['skipped'], stats['incomplete']), (3, 3, 0))
|
|
stats = ventilated([0.0] * 5, missing_end_at=(4,), zones=('a', 'b'))
|
|
self.assertEqual((stats['rooms']['a']['samples'], stats['rooms']['a']['skipped']), (1, 1))
|
|
self.assertEqual(stats['rooms']['b']['samples'], 2)
|