From 4f592777716581c78df36a8304de32e6e60c16df Mon Sep 17 00:00:00 2001 From: Fernando Gonzalez Date: Mon, 14 Sep 2026 20:59:11 -0400 Subject: [PATCH 1/2] feat: backtest the forecast against every recorded session Three model changes in one evening were each checked by replaying one or two hand-picked moments, and each found its counter-example within the hour. This scores every way the forecast can be made against what the handle actually did, over the whole history and per scenario, so a model change is judged on all of it at once. Every session is split into steady-current runs. A run that held its current for >= 3 tau yields a truth: an exponential fitted over the whole run with tau free, or (--truth last) the handle's own final three minutes - the two bracket the answer. Against it, two things are scored: the live forecast tick by tick as the current holds (trajectory basis vs model basis), and at every current change and session start the plateau that would have been predicted at the new current from each ambient on offer (sensor, previous run's trajectory, the warmer of the two, the idle handle) under each current law (I^2 prior vs fitted exponent), with the scored session left out of the fit. Errors are predicted minus actual: negative is optimistic, the direction that trips the charger. First run, 74 sessions, 40 scorable runs: the model-basis forecast is optimistic by 3-4 C on median and by more than 2 C in 79% of ticks; the mature trajectory by ~1 C; every cross-current method by 2.5-5 C at a step-down, under both truth definitions. The error grows with ambient - ~1 C on mild days, 3-5 C above 30 C - which is the ambient coefficient the drift regression had been reporting (0.35 C/C) and that its confidence interval alone did not make convincing. The fitted exponent beats I^2 everywhere scored. Whole-run tau runs 1-2 min longer than the 30 min fit windows' 11.25, which is where the trajectory's residual optimism comes from. No model change here. The next one gets scored by this before it ships. Co-Authored-By: Claude Fable 5.1 --- contrib/backtest_forecast.py | 423 ++++++++++++++++++++++++++++++++ docs/thermal-model.md | 33 +++ tests/test_backtest_forecast.py | 116 +++++++++ 3 files changed, 572 insertions(+) create mode 100644 contrib/backtest_forecast.py create mode 100644 tests/test_backtest_forecast.py diff --git a/contrib/backtest_forecast.py b/contrib/backtest_forecast.py new file mode 100644 index 0000000..8f5865e --- /dev/null +++ b/contrib/backtest_forecast.py @@ -0,0 +1,423 @@ +#!/usr/bin/env python3 +"""Score the thermal forecast against what the handle actually did, over +every recorded session. + +Three model changes in one evening were each checked by replaying one or +two hand-picked moments, and each found its counter-example within the +hour. This replays the whole history and scores every way the forecast can +be made against the plateau the handle actually reached, per scenario, so +a change to the model is judged on all of it at once. + +Two questions, scored separately: + +1. **In-run forecast.** While the current holds steady, how far off is the + projected plateau — the live forecast — as a function of how long it has + held? Trajectory basis (the run's own samples so far, the install's tau) + against model basis (an ambient plus the current law). + +2. **Cross-current prediction.** At every current change, and at every + session start, what would the plateau at the *new* current have been + predicted as from each ambient on offer — the LAN sensor, the ambient the + previous run's trajectory implies, the warmer of the two, the idle handle + before the session — under each current law: the I^2 prior and the + install's fitted exponent, fitted with the scored session left out. This + is the number a restore, or the end of a calibration probe, is decided + on. + +Ground truth is the plateau observed in a run that held its current for at +least --observe-tau time constants: an exponential fitted to the whole run +with tau free, the same fit the model's own parameters come from. Shorter +runs are still inputs (they imply an ambient) but never truth. + +Errors are predicted minus actual. Negative is optimistic — the handle ran +hotter than promised — which is the direction that trips the charger. + +Read-only. Point it at a copy of wallmonitor.db, or at the live file: SQLite +WAL readers do not block the writer. + +Example: + uv run python contrib/backtest_forecast.py --db ../wallmonitor.db + uv run python contrib/backtest_forecast.py --db ../wallmonitor.db --json out.json +""" + +from __future__ import annotations + +import argparse +import json +import sys +import time +from dataclasses import asdict, dataclass, field +from statistics import median + +from wallmonitor import thermal +from wallmonitor.db import Database + +MIN_RUN_S = 120.0 +TICK_S = 30.0 +FIRST_TICK_S = 120.0 +LAST_TICK_S = 1200.0 +OPTIMISTIC_C = 2.0 # a plateau under-read by this much is the dangerous kind of wrong +PROBE_CURRENT_FRAC = 0.75 +PROBE_MIN_S = 1800.0 +HOT_AMBIENT_C = 29.0 +WARM_START_C = 3.0 # handle this far above its idle level at a session's first run +BUCKETS_MIN = ((2, 5), (5, 10), (10, 20)) + + +@dataclass +class Run: + session_id: int + index: int + start_ts: float + end_ts: float + current_a: float + samples: list[tuple[float, float]] = field(repr=False) + handle_start_c: float = 0.0 + plateau_c: float | None = None # truth: see --truth + plateau_rmse_c: float | None = None + fit_c: float | None = None # exponential over the whole run, tau free + fit_tau_min: float | None = None + window_c: float | None = None # the same fit over the first PREFIX_SPAN_MIN_S — what the model's fitter sees + last_c: float | None = None # median handle over the run's last 3 min — model-free + max_c: float | None = None + sensor_ambient_c: float | None = None # LAN/car sensor at the run's start + idle_ambient_c: float | None = None # from the idle handle before the session (first run) + implied_ambient_c: float | None = None # this run's own trajectory, through the current law + kind: str = "" # cold_start | warm_start | step_down | step_up + probe: bool = False + hot: bool = False + + @property + def span_s(self) -> float: + return self.end_ts - self.start_ts + + +@dataclass +class TickScore: + session_id: int + kind: str + probe: bool + hot: bool + minutes: float + basis: str + error_c: float + + +@dataclass +class BoundaryScore: + session_id: int + kind: str + probe: bool + hot: bool + from_a: float | None + to_a: float + actual_c: float + predictions: dict[str, float] # "/" -> predicted plateau + + +# --------------------------------------------------------------------------- + + +def split_runs(rows: list[dict]) -> list[list[dict]]: + """Steady-current runs, the way the live forecast sees them: a run ends + when current leaves the band around its first sample, when the + contactor opens, or after a charging gap. Ramp-up samples fall into + runs too short to keep.""" + runs: list[list[dict]] = [] + current: list[dict] = [] + for row in rows: + amps = row.get("vehicle_current_a") or 0.0 + if not row.get("contactor_closed") or amps < 6.0 or row.get("handle_temp_c") is None: + if current: + runs.append(current) + current = [] + continue + if current: + ref = current[0]["vehicle_current_a"] + if abs(amps - ref) > max(2.0, 0.1 * ref) or row["ts"] - current[-1]["ts"] > thermal.SEGMENT_SPLIT_GAP_S: + runs.append(current) + current = [] + current.append(row) + if current: + runs.append(current) + return [run for run in runs if run[-1]["ts"] - run[0]["ts"] >= MIN_RUN_S] + + +def observe_plateau(run: Run, tau_min: float, observe_tau: float, truth: str) -> None: + """Fill in what the run actually reached. `truth` picks which reading + becomes plateau_c: "fit" — an exponential over the whole run with tau + free (needs >= observe_tau time constants; rejects a poor fit); "last" — + the handle's own median over the last three minutes, model-free (needs + >= observe_tau + 2 time constants, so the lag has worked itself out).""" + samples = run.samples + span = run.span_s + run.max_c = max(temp for _, temp in samples) + tail = [temp for ts, temp in samples if samples[-1][0] - ts <= 180.0] + run.last_c = median(tail) if tail else None + fit = thermal._fit_exponential(samples) + if fit is not None: + run.fit_tau_min, run.fit_c, run.plateau_rmse_c = fit[0] / 60.0, fit[1], fit[2] + head = [(ts, temp) for ts, temp in samples if ts - samples[0][0] <= thermal.PREFIX_SPAN_MIN_S] + if len(head) >= thermal.MIN_SEGMENT_SAMPLES: + window_fit = thermal._fit_exponential(head) + run.window_c = window_fit[1] if window_fit is not None else None + if truth == "fit": + if span >= observe_tau * tau_min * 60.0 and fit is not None and fit[2] <= thermal.MAX_FIT_RMSE_C: + run.plateau_c = fit[1] + elif span >= (observe_tau + 2.0) * tau_min * 60.0 and run.last_c is not None: + run.plateau_c = run.last_c + + +def params_without(fits: list[dict], sid: int, current_exp: float | None = None) -> thermal.ThermalParams: + """Model parameters from every fit but the scored session's own — + tau, the current exponent (or a forced one), and rise_ref under it.""" + others = [fit for fit in fits if fit["session_id"] != sid] + if current_exp is None: + current_exp, exp_fits, exp_se = thermal._fit_current_exponent(others) + else: + exp_fits, exp_se = 0, None + taus = [fit["tau_min"] for fit in others] + rises = [ + fit["rise_c"] * (thermal.REF_CURRENT_A / fit["current_a"]) ** current_exp + for fit in others + if fit.get("rise_c") is not None + ] + return thermal.ThermalParams( + tau_min=median(taus) if taus else thermal.DEFAULT_TAU_MIN, + rise_ref_c=median(rises) if rises else thermal.DEFAULT_RISE_REF_C, + tau_fits=len(taus), + rise_fits=len(rises), + current_exp=current_exp, + current_exp_fits=exp_fits, + current_exp_se=exp_se, + ) + + +def build_runs(db: Database, sess: dict, params: thermal.ThermalParams, observe_tau: float, + idle_model: thermal.IdleOffset, truth: str = "fit") -> list[Run]: + rows = db.vitals_range(sess["start_ts"] - 1, sess["end_ts"] + 1, 500_000) + runs: list[Run] = [] + for index, raw in enumerate(split_runs(rows)): + samples = [(row["ts"], row["handle_temp_c"]) for row in raw] + run = Run( + session_id=sess["id"], index=index, start_ts=samples[0][0], end_ts=samples[-1][0], + current_a=median(row["vehicle_current_a"] for row in raw), samples=samples, + handle_start_c=samples[0][1], + ) + observe_plateau(run, params.tau_min, observe_tau, truth) + measured = thermal._measured_ambient(db, run.start_ts - thermal.MEASURED_AMBIENT_WINDOW_S, run.start_ts + 60) + run.sensor_ambient_c = measured[0] if measured is not None else None + if index == 0: + run.idle_ambient_c = thermal._ambient_before(db, sess["start_ts"], idle_model) + if len(samples) >= thermal.TRAJECTORY_MIN_SAMPLES: + t_inf, _se = thermal._project_t_inf(samples, params.tau_min) + run.implied_ambient_c = t_inf - params.rise_at(run.current_a) + run.probe = run.current_a <= PROBE_CURRENT_FRAC * thermal.REF_CURRENT_A and run.span_s >= PROBE_MIN_S + ambient_for_hot = run.sensor_ambient_c if run.sensor_ambient_c is not None else run.idle_ambient_c + run.hot = ambient_for_hot is not None and ambient_for_hot >= HOT_AMBIENT_C + if index == 0: + baseline = ambient_for_hot + idle_level = thermal.idle_handle_c(baseline, idle_model) if baseline is not None else None + run.kind = ( + "warm_start" if idle_level is not None and run.handle_start_c > idle_level + WARM_START_C + else "cold_start" + ) + else: + run.kind = "step_down" if run.current_a < runs[-1].current_a else "step_up" + runs.append(run) + return runs + + +def score_ticks(run: Run, params: thermal.ThermalParams) -> list[TickScore]: + """The live forecast, tick by tick, against the plateau this run reached.""" + if run.plateau_c is None: + return [] + scores: list[TickScore] = [] + t0 = run.start_ts + model_sensor = ( + run.sensor_ambient_c + params.rise_at(run.current_a) if run.sensor_ambient_c is not None else None + ) + tick = FIRST_TICK_S + while tick <= min(LAST_TICK_S, run.span_s): + window = [(ts, temp) for ts, temp in run.samples if ts - t0 <= tick] + minutes = tick / 60.0 + if len(window) >= thermal.TRAJECTORY_MIN_SAMPLES: + t_inf, _se = thermal._project_t_inf(window, params.tau_min) + scores.append(TickScore(run.session_id, run.kind, run.probe, run.hot, minutes, "trajectory", t_inf - run.plateau_c)) + if model_sensor is not None: + scores.append(TickScore(run.session_id, run.kind, run.probe, run.hot, minutes, "model/sensor", model_sensor - run.plateau_c)) + tick += TICK_S + return scores + + +def score_boundary(run: Run, prev: Run | None, laws: dict[str, thermal.ThermalParams]) -> BoundaryScore | None: + """What each ambient, under each current law, would have predicted for + the plateau at this run's current, from what was known when it started.""" + if run.plateau_c is None: + return None + ambients: dict[str, float] = {} + if run.sensor_ambient_c is not None: + ambients["sensor"] = run.sensor_ambient_c + if prev is None: + if run.idle_ambient_c is not None: + ambients["idle"] = run.idle_ambient_c + else: + if prev.implied_ambient_c is not None: + ambients["implied"] = prev.implied_ambient_c + if run.sensor_ambient_c is not None: + ambients["warmer"] = max(run.sensor_ambient_c, prev.implied_ambient_c) + if not ambients: + return None + predictions = { + f"{name}/{law}": ambient + params.rise_at(run.current_a) + for name, ambient in ambients.items() + for law, params in laws.items() + } + return BoundaryScore( + run.session_id, run.kind, run.probe, run.hot, + prev.current_a if prev is not None else None, run.current_a, run.plateau_c, predictions, + ) + + +# --------------------------------------------------------------------------- + + +def _pct(values: list[float], p: float) -> float: + ordered = sorted(values) + return ordered[min(len(ordered) - 1, int(round(p * (len(ordered) - 1))))] + + +def _stats(errors: list[float]) -> str: + if not errors: + return f"{'-':>5} {'-':>6} {'-':>6} {'-':>6} {'-':>5}" + absolute = [abs(e) for e in errors] + optimistic = sum(1 for e in errors if e < -OPTIMISTIC_C) / len(errors) + return f"{len(errors):>5} {median(errors):>+6.2f} {median(absolute):>6.2f} {_pct(absolute, 0.9):>6.2f} {optimistic:>5.0%}" + + +STATS_HEADER = f"{'n':>5} {'bias':>6} {'|med|':>6} {'|p90|':>6} {'opt':>5}" + + +def report(runs: list[Run], ticks: list[TickScore], boundaries: list[BoundaryScore], + params: thermal.ThermalParams, observe_tau: float) -> None: + observed = [run for run in runs if run.plateau_c is not None] + print( + f"model: tau {params.tau_min:.2f} min, rise {params.rise_ref_c:.1f} C at {thermal.REF_CURRENT_A:g} A, " + f"n = {params.current_exp:.2f} ({params.current_exp_fits} fits)" + ) + print( + f"runs: {len(runs)} steady-current runs in {len({run.session_id for run in runs})} sessions; " + f"{len(observed)} held >= {observe_tau:g} tau and are scored as truth" + ) + kinds = sorted({run.kind for run in observed}) + print(" by kind:", ", ".join(f"{kind} {sum(1 for r in observed if r.kind == kind)}" for kind in kinds), + f"| probe {sum(1 for r in observed if r.probe)}, hot {sum(1 for r in observed if r.hot)}") + print() + print("== In-run forecast: projected plateau vs observed, by minutes at steady current") + print(" (bias = median signed error, predicted - actual; opt = share optimistic by > " + f"{OPTIMISTIC_C:g} C)") + bases = sorted({tick.basis for tick in ticks}) + print(f"{'basis':<14}" + "".join(f" {lo:>2}-{hi:<2} min {STATS_HEADER}" for lo, hi in BUCKETS_MIN)) + for basis in bases: + row = f"{basis:<14}" + for lo, hi in BUCKETS_MIN: + errors = [t.error_c for t in ticks if t.basis == basis and lo <= t.minutes < hi] + row += f" {'':>10}{_stats(errors)}" + print(row) + print() + print(" trajectory basis by scenario, 5-10 min:") + for kind in kinds: + errors = [t.error_c for t in ticks if t.basis == "trajectory" and t.kind == kind and 5 <= t.minutes < 10] + print(f" {kind:<12} {_stats(errors)}") + for flag in ("probe", "hot"): + errors = [t.error_c for t in ticks if t.basis == "trajectory" and getattr(t, flag) and 5 <= t.minutes < 10] + print(f" {flag:<12} {_stats(errors)}") + print() + print("== Cross-current prediction: plateau at the new current, from what was known at the change") + print(" ambient: sensor = LAN/car sensor at the change; implied = previous run's trajectory through the") + print(" current law; warmer = max(sensor, implied); idle = handle before the session (session start only)") + print(" law: I2 = the I^2 prior; n = the install's fitted exponent, scored session left out") + methods = sorted({m for b in boundaries for m in b.predictions}) + groups: list[tuple[str, list[BoundaryScore]]] = [(kind, [b for b in boundaries if b.kind == kind]) for kind in kinds] + groups += [("probe", [b for b in boundaries if b.probe]), ("hot", [b for b in boundaries if b.hot]), + ("all", boundaries)] + for label, group in groups: + if not group: + continue + print(f" {label} ({len(group)} changes) {STATS_HEADER}") + for method in methods: + errors = [b.predictions[method] - b.actual_c for b in group if method in b.predictions] + if errors: + print(f" {method:<14} {_stats(errors)}") + print() + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__.splitlines()[0]) + parser.add_argument("--db", required=True, help="path to a wallmonitor.db (read-only)") + parser.add_argument("--lookback-days", type=float, default=180.0) + parser.add_argument("--observe-tau", type=float, default=3.0, + help="a run must hold this many time constants for its plateau to count as truth (default %(default)s)") + parser.add_argument("--truth", choices=("fit", "last"), default="fit", + help="what counts as the plateau a run reached: an exponential fitted over the whole run " + "(fit), or the handle's own last three minutes, needing two more time constants (last)") + parser.add_argument("--json", help="also write every run, tick score and boundary score here") + parser.add_argument("--verbose", action="store_true", help="list every run") + args = parser.parse_args(argv) + + db = Database(args.db) + now = time.time() + fits = thermal.fit_sessions(db, now, lookback_days=args.lookback_days) + params = thermal.fit_history(db, now, fits=fits) + idle_model = thermal.load_idle_offset(db) + sessions = [ + sess for sess in db.sessions_range(now - args.lookback_days * 86400, now) + if sess.get("end_ts") and (sess.get("charging_s") or 0) >= thermal.MIN_SEGMENT_S + ] + sessions.sort(key=lambda sess: sess["start_ts"]) + + runs: list[Run] = [] + ticks: list[TickScore] = [] + boundaries: list[BoundaryScore] = [] + for sess in sessions: + laws = { + "n": params_without(fits, sess["id"]), + "I2": params_without(fits, sess["id"], current_exp=thermal.DEFAULT_CURRENT_EXP), + } + session_runs = build_runs(db, sess, laws["n"], args.observe_tau, idle_model, args.truth) + for index, run in enumerate(session_runs): + ticks.extend(score_ticks(run, laws["n"])) + boundary = score_boundary(run, session_runs[index - 1] if index else None, laws) + if boundary is not None: + boundaries.append(boundary) + if args.verbose: + fmt = lambda value: f"{value:5.1f}" if value is not None else " -" # noqa: E731 + print( + f"s{run.session_id:<4} run{run.index} {run.kind:<10} {run.current_a:5.1f}A " + f"{run.span_s / 60:5.0f}min handle {run.handle_start_c:5.1f} -> last {fmt(run.last_c)} " + f"max {fmt(run.max_c)} | fit {fmt(run.fit_c)} (tau {fmt(run.fit_tau_min)}) " + f"30min-fit {fmt(run.window_c)} | truth {fmt(run.plateau_c)} | " + f"sensor {fmt(run.sensor_ambient_c)} implied {fmt(run.implied_ambient_c)}" + f"{' probe' if run.probe else ''}{' hot' if run.hot else ''}" + ) + runs.extend(session_runs) + + report(runs, ticks, boundaries, params, args.observe_tau) + if args.json: + with open(args.json, "w") as fh: + json.dump( + { + "model": params.as_dict(), + "runs": [{k: v for k, v in asdict(run).items() if k != "samples"} for run in runs], + "ticks": [asdict(tick) for tick in ticks], + "boundaries": [asdict(b) for b in boundaries], + }, + fh, + indent=1, + ) + print(f"wrote {args.json}", file=sys.stderr) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/docs/thermal-model.md b/docs/thermal-model.md index f9d02d5..a13301d 100644 --- a/docs/thermal-model.md +++ b/docs/thermal-model.md @@ -160,6 +160,39 @@ alert-40 raise to within seconds. `/api/thermal` returns the fitted model, the live forecast, every per-segment fit, and the drift verdict. +## Measuring the forecast + +Every change to the model is judged against the whole recorded history, not +against the incident that prompted it: + +```bash +uv run python contrib/backtest_forecast.py --db /path/to/wallmonitor.db +``` + +The tool splits every session into steady-current runs, takes the plateau +each long-enough run actually reached as truth (an exponential fitted over +the whole run, or with `--truth last` the handle's own final minutes — the +two bracket the answer), and scores two things against it: the live +forecast tick by tick as the current holds, and — at every current change +and session start — the plateau that would have been predicted at the new +current from each ambient on offer under each current law, with the scored +session left out of the fit. Errors are predicted minus actual; negative is +optimistic, the direction that trips the charger. + +What it showed on the first run, over 74 sessions (40 runs long enough to +score): the model-basis forecast was optimistic by **3–4 °C** on median and +by more than 2 °C four times in five; the mature trajectory forecast by +about 1 °C; every cross-current method by 2.5–5 °C at a step-down. The +error grows with ambient — about 1 °C on mild days, 3–5 °C above 30 °C — +which is the same ambient coefficient the degradation watch's regression +had been reporting (0.35 °C/°C) and that the confidence interval alone did +not make convincing. The sensor reads the garage air; the handle's +environment runs hotter than the air by an amount that grows with the +heat. The fitted current exponent beat the I² prior everywhere it was +scored, and the whole-run τ ran 1–2 min longer than the 30-minute fit +windows' 11.25, which is where the trajectory's residual optimism comes +from. + ## Degradation watch The same per-segment fits feed a trend. Rising heat at unchanged current diff --git a/tests/test_backtest_forecast.py b/tests/test_backtest_forecast.py new file mode 100644 index 0000000..d8f7b84 --- /dev/null +++ b/tests/test_backtest_forecast.py @@ -0,0 +1,116 @@ +"""The forecast backtest on a synthetic history whose truth is known: a +cold-start run at full rate followed by a step-down, both long enough to +plateau. The tool must find both runs, read both plateaus, score the live +forecast against them, and score the cross-current prediction at the step +from what was known at the change.""" + +import importlib.util +import math +import pathlib +import sys +import time + +import pytest + +from wallmonitor import thermal +from wallmonitor.db import Database + +spec = importlib.util.spec_from_file_location( + "backtest_forecast", pathlib.Path(__file__).parent.parent / "contrib" / "backtest_forecast.py" +) +bt = importlib.util.module_from_spec(spec) +sys.modules[spec.name] = bt +spec.loader.exec_module(bt) + +TAU_S = 720.0 +RISE_REF = 36.0 + + +@pytest.fixture +def db(tmp_path): + database = Database(str(tmp_path / "test.db")) + yield database + database.close() + + +def _seed(db, start_ts, ambient_c, steps, dt=10.0): + """Idle lead-in, then charging segments [(amps, seconds), ...] that + follow the first-order model with the I^2 law, integrated so a + step-down cools toward its lower plateau. Sensor samples every minute.""" + ts = start_ts - 1800 + while ts < start_ts: + db.insert_vitals(ts, { + "vehicle_connected": 0, "contactor_closed": 0, "vehicle_current_a": 0.0, + "handle_temp_c": round(thermal.idle_handle_c(ambient_c), 2), "pcba_temp_c": 38.0, "mcu_temp_c": 46.0, + }, None, 0.0) + ts += dt + for minute in range(-30, int(sum(s for _, s in steps) / 60) + 1): + db.insert_ambient(start_ts + 60 * minute, ambient_c) + sid = db.start_session(start_ts) + temp = thermal.idle_handle_c(ambient_c) + ts = start_ts + for amps, seconds in steps: + end = ts + seconds + while ts < end: + db.insert_vitals(ts, { + "vehicle_connected": 1, "contactor_closed": 1, "vehicle_current_a": amps, + "handle_temp_c": round(temp, 3), "pcba_temp_c": 55.0, "mcu_temp_c": 50.0, + }, sid, amps * 233.0) + t_inf = ambient_c + RISE_REF * (amps / thermal.REF_CURRENT_A) ** 2 + temp += dt * (t_inf - temp) / TAU_S + ts += dt + db.close_session(sid, ts, "vehicle_disconnected") + return sid + + +def test_backtest_scores_a_cold_start_and_a_step_down(db, tmp_path, capsys): + now = time.time() + # History for the model to fit from, then the session under test. + for i in range(4): + _seed(db, now - (8 - i) * 7200, ambient_c=25.0, steps=[(48.6, 2400)]) + sid = _seed(db, now - 7200, ambient_c=25.0, steps=[(48.6, 2700), (40.0, 2700)]) + + out = tmp_path / "bt.json" + assert bt.main(["--db", str(tmp_path / "test.db"), "--json", str(out)]) == 0 + text = capsys.readouterr().out + assert "cold_start" in text and "step_down" in text + + import json + data = json.loads(out.read_text()) + runs = [r for r in data["runs"] if r["session_id"] == sid] + assert [r["kind"] for r in runs] == ["cold_start", "step_down"] + assert [round(r["current_a"], 1) for r in runs] == [48.6, 40.0] + # Both runs held >= 3 tau, so both plateaus are observed, and they are + # the seeded ones. + expected = [25.0 + RISE_REF * (a / 48.0) ** 2 for a in (48.6, 40.0)] + for run, plateau in zip(runs, expected): + assert run["plateau_c"] is not None and abs(run["plateau_c"] - plateau) < 1.0 + # The live forecast converges onto the plateau as the run matures. + late = [t["error_c"] for t in data["ticks"] + if t["session_id"] == sid and t["basis"] == "trajectory" and 10 <= t["minutes"] < 20] + assert late and max(abs(e) for e in late) < 1.5 + # At the step-down, every ambient on offer predicts the lower plateau — + # the history is I^2 and isothermal, so sensor, implied and warmer agree. + step = [b for b in data["boundaries"] if b["session_id"] == sid and b["kind"] == "step_down"] + assert len(step) == 1 and step[0]["from_a"] == pytest.approx(48.6, abs=0.1) + for method, predicted in step[0]["predictions"].items(): + assert abs(predicted - step[0]["actual_c"]) < 1.5, (method, predicted, step[0]["actual_c"]) + assert {m.split("/")[0] for m in step[0]["predictions"]} == {"sensor", "implied", "warmer"} + assert {m.split("/")[1] for m in step[0]["predictions"]} == {"I2", "n"} + # ...and the session start is scored from the sensor and the idle handle. + start = [b for b in data["boundaries"] if b["session_id"] == sid and b["kind"] == "cold_start"] + assert len(start) == 1 and start[0]["from_a"] is None + assert {m.split("/")[0] for m in start[0]["predictions"]} == {"sensor", "idle"} + + +def test_split_runs_drops_ramp_samples_and_splits_on_a_current_change(): + def row(i, amps, closed=1): + return {"ts": 1000.0 + 2.0 * i, "contactor_closed": closed, "vehicle_current_a": amps, "handle_temp_c": 30.0} + + rows = [row(i, a) for i, a in enumerate([6.1, 14.6, 23.9, 31.1, 32.6, 32.7, 32.9, 37.6, 43.3] + [32.5] * 200 + [44.0] * 100)] + runs = bt.split_runs(rows) + assert [len(r) for r in runs] == [200, 100] + assert runs[0][0]["vehicle_current_a"] == 32.5 and runs[1][0]["vehicle_current_a"] == 44.0 + # A contactor drop ends a run; a short remainder is dropped. + rows = [row(i, 48.0) for i in range(100)] + [row(100, 0.0, closed=0)] + [row(101 + i, 48.0) for i in range(10)] + assert [len(r) for r in bt.split_runs(rows)] == [100] From 9dfce2615f073d787c1288d366c7b11c08bed8bc Mon Sep 17 00:00:00 2001 From: Fernando Gonzalez Date: Mon, 14 Sep 2026 21:46:37 -0400 Subject: [PATCH 2/2] feat: backtest truths are only free plateaus; cut-short runs score against their own projection The history is not a neutral sample. The current and ambient every run happened at were chosen by the controller and the charger on the strength of this same model, which contaminates the truths three ways: a run whose current the charger trimmed has a plateau that is a setpoint, not an equilibrium; a run whose true plateau lay above the trip point tripped and ended before it could become truth, censoring exactly the optimistic errors the tool exists to find; and because the controller caps on hot days, hot and low-current arrive together, so an ambient effect and a current-law error wear the same signature. So: regulated runs (sag > 1.5%) are excluded from truth and counted; every run that tripped is listed with what each method predicted for it; and every free-running run cut short after >= 1 tau - capped, tripped, or ended - is scored against its own trajectory projection, SE alongside, as a proxy truth. The projection needs no ambient and no current law. On 74 sessions that reverses the first pass. Clean truths alone read the model as 3-4 C optimistic, worst on hot days, and an ambient term fitted to that looked like the fix. De-censored: at full-rate cold starts, hot days included, the deployed model is unbiased (+0.3 C from the sensor, |med| 0.7, 0% optimistic by > 2 C over 7 runs; +0.7 from the idle proxy over 18). The optimism is a history effect - a step-down right after a hot full-rate run reads 2.8-3.6 C optimistic from the sensor, and 0.7-1.4 C from the previous run's implied ambient, which carries the heat-soaked cable - and the deployed configuration (caps from the implied ambient, restores from the sensor) already sits on the right side of it. The ambient term is parked on wip/ambient-term with the evidence against it. Co-Authored-By: Claude Fable 5.1 --- contrib/backtest_forecast.py | 173 ++++++++++++++++++++++++++++------- docs/thermal-model.md | 56 +++++++++--- 2 files changed, 184 insertions(+), 45 deletions(-) diff --git a/contrib/backtest_forecast.py b/contrib/backtest_forecast.py index 8f5865e..61b4cf0 100644 --- a/contrib/backtest_forecast.py +++ b/contrib/backtest_forecast.py @@ -26,8 +26,26 @@ Ground truth is the plateau observed in a run that held its current for at least --observe-tau time constants: an exponential fitted to the whole run -with tau free, the same fit the model's own parameters come from. Shorter -runs are still inputs (they imply an ambient) but never truth. +with tau free, the same fit the model's own parameters come from (or, with +--truth last, the handle's own final minutes). Shorter runs are still +inputs (they imply an ambient) but never truth. + +The history is not a neutral sample, and the tool says so rather than +pretending. The charger trims current itself as the handle nears the trip +point and then holds it there, so a plateau from a run whose current +sagged is a setpoint, not an equilibrium: those runs are excluded from +truth and counted. A run whose true plateau lay above the trip point +tripped, folded back and ended, so it never lasted long enough to become +truth: that censors exactly the optimistic errors this tool exists to +find. Two things are done about it. Every run that tripped is listed with +what each method predicted for it — a prediction under the trip point is a +miss the tables cannot show. And every free-running run that was cut short +(capped, tripped, or simply ended) after at least one time constant is +scored against its own trajectory projection, with that projection's +standard error reported alongside, as a proxy truth: the projection needs +no ambient and no current law, so it is an independent reading of where +the run was heading. The hot-day full-rate population lives almost +entirely in that table. Errors are predicted minus actual. Negative is optimistic — the handle ran hotter than promised — which is the direction that trips the charger. @@ -62,6 +80,8 @@ HOT_AMBIENT_C = 29.0 WARM_START_C = 3.0 # handle this far above its idle level at a session's first run BUCKETS_MIN = ((2, 5), (5, 10), (10, 20)) +PROXY_MIN_TAU = 1.0 # a cut-short run is projected only once it has run this long +PROXY_MAX_SE_C = 1.5 # ...and only if the projection is this sure of itself @dataclass @@ -86,6 +106,11 @@ class Run: kind: str = "" # cold_start | warm_start | step_down | step_up probe: bool = False hot: bool = False + sag_a: float = 0.0 + free: bool = True # current held flat end to end: the plateau is the connector's own equilibrium + tripped: bool = False # alert 40 raised inside this run: its plateau was above the trip point + projected_c: float | None = None # proxy truth for a cut-short run: its own trajectory projection + projected_se_c: float | None = None @property def span_s(self) -> float: @@ -113,6 +138,7 @@ class BoundaryScore: to_a: float actual_c: float predictions: dict[str, float] # "/" -> predicted plateau + se_c: float | None = None # set when actual_c is a projection rather than an observation # --------------------------------------------------------------------------- @@ -152,6 +178,13 @@ def observe_plateau(run: Run, tau_min: float, observe_tau: float, truth: str) -> samples = run.samples span = run.span_s run.max_c = max(temp for _, temp in samples) + if not run.free or run.tripped: + # A regulated run's plateau is a setpoint the charger chose; a run + # that tripped has a plateau above the trip point that was never + # observed. Neither is a truth the forecast can be scored against — + # both are counted, because leaving them out silently censors + # exactly the optimistic errors this tool exists to find. + return tail = [temp for ts, temp in samples if samples[-1][0] - ts <= 180.0] run.last_c = median(tail) if tail else None fit = thermal._fit_exponential(samples) @@ -168,6 +201,23 @@ def observe_plateau(run: Run, tau_min: float, observe_tau: float, truth: str) -> run.plateau_c = run.last_c +def project_cut_short(run: Run, tau_min: float) -> None: + """Proxy truth for a run that ended before its plateau: where its own + trajectory was heading, with the projection's standard error. Only for + free-running runs — a regulated run's trajectory flattens because the + current fell — that ran at least PROXY_MIN_TAU and project with an SE + under PROXY_MAX_SE_C. A run that tripped qualifies: the projection is + the plateau it would have reached had the charger let it.""" + if run.plateau_c is not None or not run.free or run.span_s < PROXY_MIN_TAU * tau_min * 60.0: + return + if len(run.samples) < thermal.TRAJECTORY_MIN_SAMPLES: + return + t_inf, se = thermal._project_t_inf(run.samples, tau_min) + if se is None or se > PROXY_MAX_SE_C: + return + run.projected_c, run.projected_se_c = t_inf, se + + def params_without(fits: list[dict], sid: int, current_exp: float | None = None) -> thermal.ThermalParams: """Model parameters from every fit but the scored session's own — tau, the current exponent (or a forced one), and rise_ref under it.""" @@ -204,7 +254,14 @@ def build_runs(db: Database, sess: dict, params: thermal.ThermalParams, observe_ current_a=median(row["vehicle_current_a"] for row in raw), samples=samples, handle_start_c=samples[0][1], ) + run.sag_a = thermal._current_sag_a(raw) + run.free = thermal._free_plateau(raw, run.sag_a) + run.tripped = any( + alert.get("alert") == "40" and run.start_ts <= alert["first_ts"] <= run.end_ts + 60 + for alert in db.alerts_range(run.start_ts, run.end_ts + 60) + ) observe_plateau(run, params.tau_min, observe_tau, truth) + project_cut_short(run, params.tau_min) measured = thermal._measured_ambient(db, run.start_ts - thermal.MEASURED_AMBIENT_WINDOW_S, run.start_ts + 60) run.sensor_ambient_c = measured[0] if measured is not None else None if index == 0: @@ -250,32 +307,39 @@ def score_ticks(run: Run, params: thermal.ThermalParams) -> list[TickScore]: return scores -def score_boundary(run: Run, prev: Run | None, laws: dict[str, thermal.ThermalParams]) -> BoundaryScore | None: - """What each ambient, under each current law, would have predicted for - the plateau at this run's current, from what was known when it started.""" - if run.plateau_c is None: - return None +def predictions_for(run: Run, prev: Run | None, laws: dict[str, thermal.ThermalParams]) -> dict[str, float]: + """What each ambient, under each current law, predicts for the plateau + at this run's current, from what was known when it started.""" ambients: dict[str, float] = {} if run.sensor_ambient_c is not None: ambients["sensor"] = run.sensor_ambient_c if prev is None: if run.idle_ambient_c is not None: ambients["idle"] = run.idle_ambient_c - else: - if prev.implied_ambient_c is not None: - ambients["implied"] = prev.implied_ambient_c - if run.sensor_ambient_c is not None: - ambients["warmer"] = max(run.sensor_ambient_c, prev.implied_ambient_c) - if not ambients: - return None - predictions = { + elif prev.implied_ambient_c is not None: + ambients["implied"] = prev.implied_ambient_c + if run.sensor_ambient_c is not None: + ambients["warmer"] = max(run.sensor_ambient_c, prev.implied_ambient_c) + return { f"{name}/{law}": ambient + params.rise_at(run.current_a) for name, ambient in ambients.items() for law, params in laws.items() } + + +def score_boundary(run: Run, prev: Run | None, laws: dict[str, thermal.ThermalParams]) -> BoundaryScore | None: + """The cross-current prediction at this run's start against what the + run reached — its observed plateau, or for a cut-short run its own + projection (flagged by se_c).""" + actual, se = (run.plateau_c, None) if run.plateau_c is not None else (run.projected_c, run.projected_se_c) + if actual is None: + return None + predictions = predictions_for(run, prev, laws) + if not predictions: + return None return BoundaryScore( run.session_id, run.kind, run.probe, run.hot, - prev.current_a if prev is not None else None, run.current_a, run.plateau_c, predictions, + prev.current_a if prev is not None else None, run.current_a, actual, predictions, se, ) @@ -298,16 +362,46 @@ def _stats(errors: list[float]) -> str: STATS_HEADER = f"{'n':>5} {'bias':>6} {'|med|':>6} {'|p90|':>6} {'opt':>5}" +def _boundary_table(boundaries: list[BoundaryScore], kinds: list[str], with_se: bool = False) -> None: + methods = sorted({m for b in boundaries for m in b.predictions}) + groups: list[tuple[str, list[BoundaryScore]]] = [(kind, [b for b in boundaries if b.kind == kind]) for kind in kinds] + groups += [("probe", [b for b in boundaries if b.probe]), ("hot", [b for b in boundaries if b.hot]), + ("all", boundaries)] + for label, group in groups: + if not group: + continue + se_note = f", projection SE median {median(b.se_c for b in group):.2f} C" if with_se else "" + print(f" {label} ({len(group)} changes{se_note}) {STATS_HEADER}") + for method in methods: + errors = [b.predictions[method] - b.actual_c for b in group if method in b.predictions] + if errors: + print(f" {method:<14} {_stats(errors)}") + + def report(runs: list[Run], ticks: list[TickScore], boundaries: list[BoundaryScore], - params: thermal.ThermalParams, observe_tau: float) -> None: + params: thermal.ThermalParams, observe_tau: float, + trip_predictions: dict[tuple[int, int], dict[str, float]]) -> None: observed = [run for run in runs if run.plateau_c is not None] print( f"model: tau {params.tau_min:.2f} min, rise {params.rise_ref_c:.1f} C at {thermal.REF_CURRENT_A:g} A, " f"n = {params.current_exp:.2f} ({params.current_exp_fits} fits)" ) + long_enough = [run for run in runs if run.span_s >= observe_tau * params.tau_min * 60.0] + regulated = [run for run in long_enough if not run.free] + tripped = [run for run in runs if run.tripped] + projected = [run for run in runs if run.projected_c is not None] print( f"runs: {len(runs)} steady-current runs in {len({run.session_id for run in runs})} sessions; " - f"{len(observed)} held >= {observe_tau:g} tau and are scored as truth" + f"{len(long_enough)} held >= {observe_tau:g} tau, of which {len(observed)} are scored as truth" + ) + print( + f" not truth: {len(regulated)} regulated (current sagged > {thermal.FREE_CURRENT_SAG_FRAC:.1%} — the " + f"charger chose that plateau); {len(tripped)} runs tripped (plateau above {thermal.TRIP_HANDLE_C:g} C, " + f"never observed)" + ) + print( + f" proxy truth: {len(projected)} free-running runs cut short after >= {PROXY_MIN_TAU:g} tau, scored " + f"against their own projection (SE <= {PROXY_MAX_SE_C:g} C) in the last table" ) kinds = sorted({run.kind for run in observed}) print(" by kind:", ", ".join(f"{kind} {sum(1 for r in observed if r.kind == kind)}" for kind in kinds), @@ -337,18 +431,25 @@ def report(runs: list[Run], ticks: list[TickScore], boundaries: list[BoundarySco print(" ambient: sensor = LAN/car sensor at the change; implied = previous run's trajectory through the") print(" current law; warmer = max(sensor, implied); idle = handle before the session (session start only)") print(" law: I2 = the I^2 prior; n = the install's fitted exponent, scored session left out") - methods = sorted({m for b in boundaries for m in b.predictions}) - groups: list[tuple[str, list[BoundaryScore]]] = [(kind, [b for b in boundaries if b.kind == kind]) for kind in kinds] - groups += [("probe", [b for b in boundaries if b.probe]), ("hot", [b for b in boundaries if b.hot]), - ("all", boundaries)] - for label, group in groups: - if not group: - continue - print(f" {label} ({len(group)} changes) {STATS_HEADER}") - for method in methods: - errors = [b.predictions[method] - b.actual_c for b in group if method in b.predictions] - if errors: - print(f" {method:<14} {_stats(errors)}") + clean = [b for b in boundaries if b.se_c is None] + proxy = [b for b in boundaries if b.se_c is not None] + _boundary_table(clean, kinds) + if tripped: + print() + print("== Runs that tripped: what each method predicted for a plateau that was in fact above the trip point") + print(" (a prediction under 65 C is an optimistic miss the tables above cannot show)") + for run in tripped: + preds = trip_predictions.get((run.session_id, run.index), {}) + summary = ", ".join(f"{m} {v:.1f}" for m, v in sorted(preds.items())) + print(f" s{run.session_id:<4} run{run.index} {run.kind:<10} {run.current_a:5.1f}A " + f"{run.span_s / 60:4.0f}min max {run.max_c:.1f} | {summary or 'no ambient known'}") + if proxy: + print() + print("== Cut-short runs: the same prediction against the run's own trajectory projection (proxy truth)") + print(" These are the runs the controller capped or the charger tripped — the population the clean") + print(" tables cannot contain. The projection needs no ambient and no current law.") + proxy_kinds = sorted({b.kind for b in proxy}) + _boundary_table(proxy, proxy_kinds, with_se=True) print() @@ -379,6 +480,7 @@ def main(argv: list[str] | None = None) -> int: runs: list[Run] = [] ticks: list[TickScore] = [] boundaries: list[BoundaryScore] = [] + trip_predictions: dict[tuple[int, int], dict[str, float]] = {} for sess in sessions: laws = { "n": params_without(fits, sess["id"]), @@ -386,10 +488,13 @@ def main(argv: list[str] | None = None) -> int: } session_runs = build_runs(db, sess, laws["n"], args.observe_tau, idle_model, args.truth) for index, run in enumerate(session_runs): + prev = session_runs[index - 1] if index else None ticks.extend(score_ticks(run, laws["n"])) - boundary = score_boundary(run, session_runs[index - 1] if index else None, laws) + boundary = score_boundary(run, prev, laws) if boundary is not None: boundaries.append(boundary) + if run.tripped: + trip_predictions[(run.session_id, run.index)] = predictions_for(run, prev, laws) if args.verbose: fmt = lambda value: f"{value:5.1f}" if value is not None else " -" # noqa: E731 print( @@ -397,12 +502,14 @@ def main(argv: list[str] | None = None) -> int: f"{run.span_s / 60:5.0f}min handle {run.handle_start_c:5.1f} -> last {fmt(run.last_c)} " f"max {fmt(run.max_c)} | fit {fmt(run.fit_c)} (tau {fmt(run.fit_tau_min)}) " f"30min-fit {fmt(run.window_c)} | truth {fmt(run.plateau_c)} | " - f"sensor {fmt(run.sensor_ambient_c)} implied {fmt(run.implied_ambient_c)}" + f"sensor {fmt(run.sensor_ambient_c)} implied {fmt(run.implied_ambient_c)} sag {run.sag_a:+.1f}" f"{' probe' if run.probe else ''}{' hot' if run.hot else ''}" + f"{'' if run.free else ' REGULATED'}{' TRIPPED' if run.tripped else ''}" + f"{f' proj {run.projected_c:.1f}±{run.projected_se_c:.2f}' if run.projected_c is not None else ''}" ) runs.extend(session_runs) - report(runs, ticks, boundaries, params, args.observe_tau) + report(runs, ticks, boundaries, params, args.observe_tau, trip_predictions) if args.json: with open(args.json, "w") as fh: json.dump( diff --git a/docs/thermal-model.md b/docs/thermal-model.md index a13301d..8574eb5 100644 --- a/docs/thermal-model.md +++ b/docs/thermal-model.md @@ -179,18 +179,50 @@ current from each ambient on offer under each current law, with the scored session left out of the fit. Errors are predicted minus actual; negative is optimistic, the direction that trips the charger. -What it showed on the first run, over 74 sessions (40 runs long enough to -score): the model-basis forecast was optimistic by **3–4 °C** on median and -by more than 2 °C four times in five; the mature trajectory forecast by -about 1 °C; every cross-current method by 2.5–5 °C at a step-down. The -error grows with ambient — about 1 °C on mild days, 3–5 °C above 30 °C — -which is the same ambient coefficient the degradation watch's regression -had been reporting (0.35 °C/°C) and that the confidence interval alone did -not make convincing. The sensor reads the garage air; the handle's -environment runs hotter than the air by an amount that grows with the -heat. The fitted current exponent beat the I² prior everywhere it was -scored, and the whole-run τ ran 1–2 min longer than the 30-minute fit -windows' 11.25, which is where the trajectory's residual optimism comes +**The history is not a neutral sample, and the tool has to say so.** The +current and the ambient a run happened at were chosen by the controller +and by the charger, on the strength of this same model, and that +contaminates the truths three ways. The charger trims current as the +handle nears the trip point and then *holds* it there, so a run whose +current sagged has a plateau that is a setpoint, not an equilibrium — +those are excluded and counted (20 of the 52 long-enough runs on one +install). A run whose true plateau lay above the trip point tripped, +folded back and ended before it could become truth — so the clean tables +censor exactly the optimistic errors that matter, and every run that +tripped is listed instead with what each method predicted for it. And +because the controller caps on hot days, *hot* and *low current* arrive +together in the data, so an ambient effect and a current-law error wear +the same signature. + +The correction for the censoring is to score every free-running run that +was cut short — capped, tripped, or simply ended after at least one time +constant — against **its own trajectory projection**, with the +projection's standard error reported alongside. The projection needs no +ambient and no current law, so it is an independent reading of where the +run was heading, and the hot-day full-rate population lives almost +entirely in that table. + +What it showed, once corrected, on 74 sessions: the first pass — clean +truths only — read the model-basis forecast as optimistic by 3–4 °C, worst +on hot days, and an ambient term fitted to that would have "fixed" it. The +de-censored pass reversed the finding. At full-rate cold starts, hot days +included, the model is **unbiased** (+0.3 °C median from the sensor, +|median| 0.7, none optimistic by more than 2 °C; +0.7 from the idle +proxy over 18 runs). The optimism is a **history effect**: a step-down +right after a hot full-rate run reads 2.8–3.6 °C optimistic from the +sensor, because the cable and connector are still heat-soaked from the +run before, and the previous run's trajectory-implied ambient — which +carries that state — cuts it to 0.7–1.4 °C. That is why caps are worked +from the implied ambient and restores from the sensor, and why the +[degradation watch](#the-confounder-the-fits-cant-remove) reads an +"ambient coefficient" on an install whose hot days are also its +heaviest-charging days. A model with a second, slow time constant for the +cable would carry the effect properly; it needs designed data — probes at +one current with a cold cable and a warm one — not more of the history. + +The fitted current exponent beat the I² prior in every de-censored cell +below 48 A, and the whole-run τ ran 1–2 min longer than the 30-minute fit +windows', which is where the mature trajectory's residual 1 °C comes from. ## Degradation watch