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 f5ed8586a956679673c5a68661839f480a6c315b Mon Sep 17 00:00:00 2001 From: Fernando Gonzalez Date: Mon, 14 Sep 2026 21:13:27 -0400 Subject: [PATCH 2/2] feat: fit an ambient term into the rise law, and fit tau over 4 tau windows What the backtest found, fixed, and scored again before shipping. The sensor reads the garage air. The handle's environment runs hotter than the air by an amount that grows with the heat - sun on the wall, a heat-soaked structure and cable - and the model had no term for it. Over 74 sessions the model-basis forecast was optimistic by ~1 C on mild days and 3-5 C above 30 C, in 79% of ticks by more than 2 C, and every cross-current prediction (the number a restore is decided on) by 2.5-5 C. The degradation watch's regression had been reporting the same slope as an unexplained "ambient coefficient" (0.35 C/C) for weeks; its confidence interval alone had not made it convincing. Forty runs did. So rise(I, a) = rise_ref * (I/48)^n + k * (a - 25), with n and k fitted together from the free-running fits: n profiled over a grid, k and rise_ref by least squares at each n, each term only when the history spans enough to identify it (>= 6 A, >= 4 C). rise_ref_c now means the rise at 48 A and 25 C; every fit is re-normalized under the law; every site that turned an ambient into a plateau or back goes through plateau_at / ambient_from_plateau. On this install: n = 1.52 +/- 0.25, k = 0.25 +/- 0.13; a full-rate charge settles at 66 C at 28 C ambient and 75 C at 35 C, and the sustainable current reads 45 / 43 / 41 / 37 A at 28 / 30 / 32 / 35 C. The fit window also runs to 4 tau instead of 30 min. Whole-run fits read tau 1-2 min longer than the 30 min windows did, and that gap was the trajectory forecast's residual optimism: a window that ends at 2.7 tau has seen 93% of the rise and still trades a little tau for a little rise. tau moves 11.25 -> 12.0 min. Before -> after on production, as (whole-run-fit truth / model-free truth), bias in C and share optimistic by more than 2 C: model basis, in-run -3.3/-3.7 -> -2.1/-1.2 79%/80% -> 52%/20% step-down, from sensor -4.6/-4.0 -> -2.2/-1.3 92%/100% -> 64%/20% probe end, from sensor -3.8/-3.9 -> -0.6/-0.8 67%/100% -> 0%/0% trajectory at 10-20 min -1.0/-0.6 -> -0.6/+0.4 31%/31% -> 21%/18% The same run settled the restore-ambient question the anecdotes could not: at a step-down the sensor beats the ambient the previous run's trajectory implies, and "the warmer of the two" ties the sensor within 0.1 C. The sensor stays; nothing cleverer ships. Two fits that read free-running over 30 min windows read regulated over 4 tau ones (the sag shows), so the law is fitted from 12 free fits rather than 14, and the drift regression's residual scatter rises to 1.85 C - the three off-reference free fits still disagree with each other about the law, and the watch's threshold self-calibrates to 4.75 C until monthly probes pin it. Verdict unchanged: flat. Co-Authored-By: Claude Fable 5.1 --- contrib/backtest_forecast.py | 67 +++++----- docs/thermal-model.md | 60 ++++++--- tests/test_backtest_forecast.py | 2 +- tests/test_wallmonitor.py | 41 ++++++- wallmonitor/static/app.js | 6 +- wallmonitor/thermal.py | 210 ++++++++++++++++++++++++-------- 6 files changed, 271 insertions(+), 115 deletions(-) diff --git a/contrib/backtest_forecast.py b/contrib/backtest_forecast.py index 8f5865e..f14b8d1 100644 --- a/contrib/backtest_forecast.py +++ b/contrib/backtest_forecast.py @@ -19,10 +19,10 @@ 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. + before the session — under each current law: the I^2 prior, the + install's fitted exponent, and the fitted exponent with the ambient + term, each 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 @@ -168,29 +168,19 @@ def observe_plateau(run: Run, tau_min: float, observe_tau: float, truth: str) -> 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, - ) +LAWS = { + # name -> (exponent pinned?, ambient_coef pinned?) — None means fitted + "I2": dict(exponent=thermal.DEFAULT_CURRENT_EXP, ambient_coef=thermal.DEFAULT_AMBIENT_COEF), + "n": dict(ambient_coef=thermal.DEFAULT_AMBIENT_COEF), + "n+k": dict(), +} + + +def params_without(fits: list[dict], sid: int, **law) -> thermal.ThermalParams: + """Model parameters from every fit but the scored session's own, under a + current law with the given terms pinned (see LAWS).""" + others = [dict(fit) for fit in fits if fit["session_id"] != sid] + return thermal.params_from_fits(others, **law) def build_runs(db: Database, sess: dict, params: thermal.ThermalParams, observe_tau: float, @@ -211,7 +201,7 @@ def build_runs(db: Database, sess: dict, params: thermal.ThermalParams, observe_ 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.implied_ambient_c = params.ambient_from_plateau(t_inf, 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 @@ -235,7 +225,7 @@ def score_ticks(run: Run, params: thermal.ThermalParams) -> list[TickScore]: 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 + params.plateau_at(run.current_a, run.sensor_ambient_c) if run.sensor_ambient_c is not None else None ) tick = FIRST_TICK_S while tick <= min(LAST_TICK_S, run.span_s): @@ -269,7 +259,7 @@ def score_boundary(run: Run, prev: Run | None, laws: dict[str, thermal.ThermalPa if not ambients: return None predictions = { - f"{name}/{law}": ambient + params.rise_at(run.current_a) + f"{name}/{law}": params.plateau_at(run.current_a, ambient) for name, ambient in ambients.items() for law, params in laws.items() } @@ -302,8 +292,9 @@ def report(runs: list[Run], ticks: list[TickScore], boundaries: list[BoundarySco 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)" + f"model: tau {params.tau_min:.2f} min, rise {params.rise_ref_c:.1f} C at {thermal.REF_CURRENT_A:g} A " + f"and {thermal.AMBIENT_REF_C:g} C, n = {params.current_exp:.2f}, k = {params.ambient_coef:+.3f} C/C " + f"({params.current_exp_fits} fits)" ) print( f"runs: {len(runs)} steady-current runs in {len({run.session_id for run in runs})} sessions; " @@ -336,7 +327,8 @@ def report(runs: list[Run], ticks: list[TickScore], boundaries: list[BoundarySco 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") + print(" law: I2 = the I^2 prior; n = fitted exponent; n+k = fitted exponent and ambient term;") + print(" each fitted with the 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]), @@ -380,13 +372,10 @@ def main(argv: list[str] | None = None) -> int: 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) + laws = {name: params_without(fits, sess["id"], **pins) for name, pins in LAWS.items()} + session_runs = build_runs(db, sess, laws["n+k"], args.observe_tau, idle_model, args.truth) for index, run in enumerate(session_runs): - ticks.extend(score_ticks(run, laws["n"])) + ticks.extend(score_ticks(run, laws["n+k"])) boundary = score_boundary(run, session_runs[index - 1] if index else None, laws) if boundary is not None: boundaries.append(boundary) diff --git a/docs/thermal-model.md b/docs/thermal-model.md index a13301d..7c229aa 100644 --- a/docs/thermal-model.md +++ b/docs/thermal-model.md @@ -31,9 +31,22 @@ the trip happens. forecast at an off-reference current inherited the error — including the current the amp controller was told to restore to. Once the free-running fits span ≥ 6 A of current, the exponent *n* in - rise = rise₄₈ · (I/48)ⁿ is fitted by log-log regression (`current_exp` - in `/api/thermal`, with its standard error), every fit's `rise_ref_c` is - re-normalized with it, and the model note says so. + rise = rise₄₈ · (I/48)ⁿ is fitted (`current_exp` in `/api/thermal`, with + its standard error), every fit's `rise_ref_c` is re-normalized with it, + and the model note says so. +- **And so is how it scales with ambient.** The sensor reads the garage + air; the handle's environment runs hotter than the air by an amount that + grows with the heat — sun on the wall, a heat-soaked structure and cable. + [Backtested](#measuring-the-forecast) over 74 sessions, a model without + that term was optimistic by ~1 °C on mild days and 3–5 °C above 30 °C. + So the rise carries an ambient term, k · (ambient − 25 °C), fitted + jointly with *n* once the free-running fits span ≥ 4 °C of ambient + (`ambient_coef`, with its standard error; `rise_ref_c` is the rise at + 48 A *and* 25 °C). On the install above k = 0.25 °C/°C. The + [degradation watch](#the-confounder-the-fits-cant-remove) had been + reporting the same slope as an unexplained "ambient coefficient" for + weeks; with the term in the model, that coefficient reads ~0 and the + watch is checking the model rather than doing its job for it. - **The charger is its own thermometer.** Idle, the handle sits ~1–2 °C above ambient (an ambient-dependent offset), so ambient can be read without any extra sensor. The offset model ships as a seed from one install and is @@ -179,19 +192,34 @@ 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. +What it showed on its first run, over 74 sessions: the model-basis +forecast was optimistic by **3–4 °C** on median and by more than 2 °C four +times in five; every cross-current method by 2.5–5 °C at a step-down. The +error grew with ambient — about 1 °C on mild days, 3–5 °C above 30 °C — +which was the same ambient coefficient the degradation watch's regression +had been reporting (0.35 °C/°C) and that its confidence interval alone had +not made convincing. 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 came from. + +Two model changes followed, each scored by the tool before it shipped: the +ambient term above, and fit windows of 4 τ instead of 30 min. Before → +after, as (whole-run-fit truth / model-free truth): + +| | bias, °C | optimistic by > 2 °C | +|---|---|---| +| model basis, in-run | −3.3 / −3.7 → **−2.1 / −1.2** | 79 % / 80 % → **52 % / 20 %** | +| trajectory, 10–20 min in | −1.0 / −0.6 → −0.6 / +0.4 | 31 % / 31 % → 21 % / 18 % | +| step-down, from the sensor | −4.6 / −4.0 → **−2.2 / −1.3** | 92 % / 100 % → **64 % / 20 %** | +| probe end, from the sensor | −3.8 / −3.9 → **−0.6 / −0.8** | 67 % / 100 % → **0 % / 0 %** | + +The remaining optimism at a step-down (1–2 °C) is a history effect no +static ambient carries: the cable is still warm from the higher current +that preceded it. The `SUGGEST_MARGIN_C` of 2 °C covers the median of it. +The same run settled a rule the anecdotes could not: at a step-down the +sensor's ambient beats the one the previous run's trajectory implies, and +"the warmer of the two" ties the sensor within 0.1 °C — so the sensor is +used, and nothing cleverer. ## Degradation watch diff --git a/tests/test_backtest_forecast.py b/tests/test_backtest_forecast.py index d8f7b84..9d45cce 100644 --- a/tests/test_backtest_forecast.py +++ b/tests/test_backtest_forecast.py @@ -96,7 +96,7 @@ def test_backtest_scores_a_cold_start_and_a_step_down(db, tmp_path, capsys): 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"} + assert {m.split("/")[1] for m in step[0]["predictions"]} == {"I2", "n", "n+k"} # ...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 diff --git a/tests/test_wallmonitor.py b/tests/test_wallmonitor.py index 761c5b7..2bc249e 100644 --- a/tests/test_wallmonitor.py +++ b/tests/test_wallmonitor.py @@ -336,7 +336,7 @@ def _seed_idle(db, t_from, t_to, ambient_c, dt=10.0): def _seed_thermal_session(db, start_ts, ambient_c, tau_s=720.0, rise_ref_c=36.0, amps=48.6, charge_s=1500.0, dt=10.0, ambient_end_c=None, cooldown_s=0.0, sag_to_a=None, - current_exp=2.0): + current_exp=2.0, ambient_coef=0.0): """Idle lead-in plus a charging ramp that follows the first-order model. With ambient_end_c set, ambient drifts linearly across the charge (the @@ -354,7 +354,7 @@ def _seed_thermal_session(db, start_ts, ambient_c, tau_s=720.0, rise_ref_c=36.0, _seed_idle(db, start_ts - 1800, start_ts, ambient_c, dt) sid = db.start_session(start_ts) t0_temp = thermal.idle_handle_c(ambient_c) - rise_at = rise_ref_c * (amps / thermal.REF_CURRENT_A) ** current_exp + rise_at = rise_ref_c * (amps / thermal.REF_CURRENT_A) ** current_exp + ambient_coef * (ambient_c - 25.0) integrate = ambient_end_c is not None or sag_to_a is not None temp = t0_temp ts = start_ts @@ -371,7 +371,7 @@ def _seed_thermal_session(db, start_ts, ambient_c, tau_s=720.0, rise_ref_c=36.0, if integrate: ambient_now = ambient_c if ambient_end_c is None else \ ambient_c + (ambient_end_c - ambient_c) * elapsed - rise_now = rise_ref_c * (amps_now / thermal.REF_CURRENT_A) ** current_exp + rise_now = rise_ref_c * (amps_now / thermal.REF_CURRENT_A) ** current_exp + ambient_coef * (ambient_now - 25.0) temp += dt * ((ambient_now + rise_now - temp) / tau_s) ts += dt db.close_session(sid, start_ts + charge_s, "vehicle_disconnected") @@ -410,7 +410,8 @@ async def test_thermal_fits_the_current_exponent_from_a_current_spread(db): params = thermal.fit_history(db, now, fits=fits) assert params.current_exp_fits == 6 assert abs(params.current_exp - 1.5) < 0.15 - assert params.current_exp_se is not None and params.current_exp_se < 0.15 + assert params.current_exp_se is not None and params.current_exp_se < 0.2 + assert params.ambient_coef == 0.0 and params.ambient_coef_se is None # one ambient: not identifiable assert abs(params.rise_ref_c - 36.0) < 2.0 for fit in fits: assert abs(fit["rise_ref_c"] - 36.0) < 2.5, fit @@ -419,6 +420,38 @@ async def test_thermal_fits_the_current_exponent_from_a_current_spread(db): assert all(fit["rise_c"] * (48.0 / fit["current_a"]) ** 2 > 40.0 for fit in low) +async def test_thermal_fits_the_ambient_term_and_predicts_hot_days(db): + # Backtested over 74 real sessions, the forecast was ~1 C optimistic on + # mild days and 3-5 C above 30 C: the handle's environment runs hotter + # than the garage air by an amount that grows with the heat. Seeded here + # with k = 0.3 C/C across a spread of ambients and two currents; the + # fitter must recover k and n together, normalize every fit to 48 A and + # 25 C, and the forecast must then land on a hot day it has never seen. + now = time.time() + plan = [(22.0, 48.6), (25.0, 48.6), (30.0, 48.6), (34.0, 48.6), (27.0, 40.0), (32.0, 40.0), (24.0, 40.0)] + for i, (ambient, amps) in enumerate(plan): + _seed_thermal_session(db, now - (len(plan) + 1 - i) * 7200, ambient_c=ambient, amps=amps, + current_exp=1.5, ambient_coef=0.3) + fits = thermal.fit_sessions(db, now) + params = thermal.fit_history(db, now, fits=fits) + assert params.current_exp_fits == len(plan) + assert abs(params.ambient_coef - 0.3) < 0.08, params.ambient_coef + assert params.ambient_coef_se is not None and params.ambient_coef_se < 0.08 + assert abs(params.current_exp - 1.5) < 0.2 + assert abs(params.rise_ref_c - 36.0) < 2.0 + for fit in fits: + assert abs(fit["rise_ref_c"] - 36.0) < 2.5, fit + # A 36 C day at full rate: the seeded truth is 36 + 36.4 + 3.3 = 75.7. + truth = 36.0 + 36.0 * (48.6 / 48.0) ** 1.5 + 0.3 * 11.0 + assert abs(params.plateau_at(48.6, 36.0) - truth) < 1.5 + assert abs(params.ambient_from_plateau(truth, 48.6) - 36.0) < 1.0 + # ...and the sustainable current shrinks accordingly: without k, 36 C + # leaves 27 C of headroom; with k = 0.3 it leaves 23.7. + without = thermal.ThermalParams(rise_ref_c=params.rise_ref_c, current_exp=params.current_exp) + assert thermal.sustainable_max_current(36.0, params) < thermal.sustainable_max_current(36.0, without) + assert params.safe_ambient_max_c() < without.safe_ambient_max_c() + + def test_steady_prefix_restarts_after_a_ramp_up_overshoot(): # 2026-09-14, first live calibration probe: the car overshot toward 48 A # for two samples while the 32 A cap was taking effect, then held 32.5 A diff --git a/wallmonitor/static/app.js b/wallmonitor/static/app.js index a5e4bf6..1dd0618 100644 --- a/wallmonitor/static/app.js +++ b/wallmonitor/static/app.js @@ -1145,8 +1145,10 @@ async function viewLive(root) { // The current law is a prior (I²) until this install's own fits span // enough current to measure it; say which one every forecast rests on. const expNote = model.current_exp_fits > 0 - ? ` Heat rise scales as I^${fmtNum(model.current_exp, 2)} here (fitted across ${model.current_exp_fits} free-running ` + - "sessions; the I² prior over-reads how much lower currents cool the handle)." + ? ` Heat rise scales as I^${fmtNum(model.current_exp, 2)} here` + + (model.ambient_coef ? ` and grows ${fmtNum(model.ambient_coef, 2)} °C per °C of ambient above ${fmtNum(model.ambient_ref_c, 0)} °C` : "") + + ` (fitted across ${model.current_exp_fits} free-running sessions; the I² prior over-reads how much lower currents ` + + "cool the handle, and the garage air under-reads how hot the handle's surroundings run on a hot day)." : ""; const modelNote = `Model: τ ≈ ${fmtNum(model.tau_min, 1)} min, +${fmtNum(model.rise_ref_c, 0)} °C at ${fmtNum(model.ref_current_a, 0)} A — ` + (model.fitted ? `fitted from ${model.tau_fits} recorded session ramp${model.tau_fits === 1 ? "" : "s"}.` + priorNote + expNote diff --git a/wallmonitor/thermal.py b/wallmonitor/thermal.py index 09e7699..55bc9a5 100644 --- a/wallmonitor/thermal.py +++ b/wallmonitor/thermal.py @@ -199,6 +199,21 @@ def ambient_from_idle_handle(handle_c: float, model: IdleOffset = BUILTIN_IDLE_O CURRENT_EXP_MIN_FITS = 4 CURRENT_EXP_MIN_SPAN_A = 6.0 CURRENT_EXP_RANGE = (1.0, 2.5) +CURRENT_EXP_GRID_STEP = 0.02 + +# The rise also grows with ambient: rise(I, a) = rise_ref * (I/48)^n + +# k * (a - 25). The sensor reads the garage air, and the handle's +# environment runs hotter than the air by an amount that scales with the +# heat — sun on the wall, a heat-soaked structure and cable. Backtested over +# 74 sessions with k = 0, the model-basis forecast was optimistic by ~1 C on +# mild days and 3-5 C above 30 C, in 79% of ticks by more than 2 C, and the +# degradation watch's regression had been reporting the same slope +# (0.35 C/C) as an "ambient coefficient" for weeks. So it is a model term, +# fitted from the same free-running fits once they span enough ambient, +# and rise_ref_c means the rise at 48 A and 25 C. +AMBIENT_REF_C = 25.0 +DEFAULT_AMBIENT_COEF = 0.0 +AMBIENT_COEF_MIN_SPREAD_C = 4.0 # Fit acceptance gates: a segment must actually contain a thermal ramp and # the exponential must describe it well, or it teaches the model nothing. @@ -220,7 +235,12 @@ def ambient_from_idle_handle(handle_c: float, model: IdleOffset = BUILTIN_IDLE_O # scales with the same tau estimate (PREFIX_SPAN_TAU) so a slow-tau install # is not starved of fits by a fixed cap it can never clear. MIN_SPAN_TAU = 1.8 -PREFIX_SPAN_TAU = 2.5 +# The fit window runs to 4 tau (was 2.5, which the 30 min floor made moot at +# a typical tau). Whole-run fits over the recorded history read tau 1-2 min +# longer than the 30 min windows did, and that gap is where the trajectory +# forecast's residual optimism came from: a window that ends at 2.7 tau has +# seen 93% of the rise and still trades a little tau for a little rise. +PREFIX_SPAN_TAU = 4.0 PREFIX_SPAN_MIN_S = 1800.0 # A steady run that ends within its first minute did not end; the ramp-up # was still wobbling. Vehicles overshoot on the way to a cap: one session @@ -294,19 +314,36 @@ class ThermalParams: current_exp: float = DEFAULT_CURRENT_EXP current_exp_fits: int = 0 current_exp_se: float | None = None + ambient_coef: float = DEFAULT_AMBIENT_COEF + ambient_coef_se: float | None = None @property def fitted(self) -> bool: return self.tau_fits > 0 and self.rise_fits > 0 - def rise_at(self, current_a: float) -> float: - """Steady-state rise above ambient at a charge current.""" - return self.rise_ref_c * (current_a / REF_CURRENT_A) ** self.current_exp + def rise_at(self, current_a: float, ambient_c: float = AMBIENT_REF_C) -> float: + """Steady-state rise above ambient at a charge current and ambient.""" + return ( + self.rise_ref_c * (current_a / REF_CURRENT_A) ** self.current_exp + + self.ambient_coef * (ambient_c - AMBIENT_REF_C) + ) + + def plateau_at(self, current_a: float, ambient_c: float) -> float: + return ambient_c + self.rise_at(current_a, ambient_c) + + def ambient_from_plateau(self, plateau_c: float, current_a: float) -> float: + """plateau_at inverted: the ambient at which this current settles here.""" + base = self.rise_ref_c * (current_a / REF_CURRENT_A) ** self.current_exp + return (plateau_c - base + self.ambient_coef * AMBIENT_REF_C) / (1.0 + self.ambient_coef) def current_for_rise(self, rise_c: float) -> float: - """The charge current whose steady-state rise is rise_c — rise_at inverted.""" + """The charge current whose current-dependent rise is rise_c.""" return REF_CURRENT_A * (rise_c / self.rise_ref_c) ** (1.0 / self.current_exp) + def safe_ambient_max_c(self) -> float: + """The ambient above which a full-rate charge settles at the trip point.""" + return (TRIP_HANDLE_C - self.rise_ref_c + self.ambient_coef * AMBIENT_REF_C) / (1.0 + self.ambient_coef) + # How far a fitted value may sit from the default before the dashboard # says the priors were a poor fit for this install. A heuristic, not a # statistic: 30% is roughly where the default-driven forecast's plateau @@ -352,6 +389,9 @@ def as_dict(self) -> dict: "current_exp": round(self.current_exp, 2), "current_exp_fits": self.current_exp_fits, "current_exp_se": round(self.current_exp_se, 2) if self.current_exp_se is not None else None, + "ambient_coef": round(self.ambient_coef, 3), + "ambient_coef_se": round(self.ambient_coef_se, 3) if self.ambient_coef_se is not None else None, + "ambient_ref_c": AMBIENT_REF_C, } @@ -810,51 +850,105 @@ def fit_sessions(db: Database, now: float, lookback_days: float = 120.0) -> list return fits -def _fit_current_exponent(fits: list[dict]) -> tuple[float, int, float | None]: - """The exponent n in rise = rise_ref * (I/48)^n, from free-running fits: - (n, fits used, standard error). The prior with no fits used when the - history cannot identify it — too few free-running fits, or all at one - current, where any n explains the data equally well. - - Log-log least squares: ln(rise) = ln(rise_ref) + n * ln(I/48). Only - free-running windows count — a regulated window's plateau is lower for a - reason that has nothing to do with current, and those windows cluster at - high current, which would bend n downward for the wrong reason.""" +@dataclass(frozen=True) +class CurrentLaw: + """How rise depends on current and ambient: rise = rise_ref * (I/48)^n + + k * (ambient - 25), with how many fits identified it and the standard + error of each fitted term (None for a term held at its prior).""" + + exponent: float = DEFAULT_CURRENT_EXP + ambient_coef: float = DEFAULT_AMBIENT_COEF + fits: int = 0 + exponent_se: float | None = None + ambient_coef_se: float | None = None + + +def _fit_current_law(fits: list[dict], exponent: float | None = None, + ambient_coef: float | None = None) -> CurrentLaw: + """Fit n and k from free-running fits; either can be pinned instead. + + n is profiled over a grid: for each n the rise is linear in rise_ref + and k, so ordinary least squares gives both and the n with the smallest + residual wins. Its standard error comes from the curvature of the + residual sum of squares at the optimum. A term is only fitted when the + history moved enough to identify it — >= CURRENT_EXP_MIN_SPAN_A of + current for n, >= AMBIENT_COEF_MIN_SPREAD_C of ambient for k — and + stays at its prior otherwise. Only free-running windows count: a + regulated window's plateau is lower for a reason that has nothing to do + with current or ambient, and those windows cluster at high current on + hot days, which would bend both terms for the wrong reason.""" points = [ - (math.log(fit["current_a"] / REF_CURRENT_A), math.log(fit["rise_c"])) + (fit["current_a"], fit["ambient_c"], fit["rise_c"]) for fit in fits - if fit.get("rise_c") is not None and fit["rise_c"] > 0 - and fit.get("current_a") and fit.get("free_plateau", True) + if fit.get("rise_c") is not None and fit["rise_c"] > 0 and fit.get("current_a") + and fit.get("ambient_c") is not None and fit.get("free_plateau", True) ] + pinned = CurrentLaw( + exponent=DEFAULT_CURRENT_EXP if exponent is None else exponent, + ambient_coef=DEFAULT_AMBIENT_COEF if ambient_coef is None else ambient_coef, + ) if len(points) < CURRENT_EXP_MIN_FITS: - return DEFAULT_CURRENT_EXP, 0, None - currents = [REF_CURRENT_A * math.exp(x) for x, _ in points] - if max(currents) - min(currents) < CURRENT_EXP_MIN_SPAN_A: - return DEFAULT_CURRENT_EXP, 0, None - result = _ols([y for _, y in points], [[1.0, x] for x, _ in points]) - if result is None: - return DEFAULT_CURRENT_EXP, 0, None - beta, se, _resid, _dof = result - exponent = min(CURRENT_EXP_RANGE[1], max(CURRENT_EXP_RANGE[0], beta[1])) - return exponent, len(points), se[1] - - -def fit_history(db: Database, now: float, lookback_days: float = 120.0, - fits: list[dict] | None = None) -> ThermalParams: - """Aggregate per-session fits into model parameters; defaults where thin. - - Fits the install's current exponent first and re-normalizes every fit's - rise_ref_c with it, in place — so the fits handed on to the API and the - degradation watch, and the rise_ref_c median taken here, all share one - current law. With the I^2 prior, off-reference fits carry a bias that the - watch's current term then has to absorb; with a fitted exponent that - term has nothing left to explain.""" - if fits is None: - fits = fit_sessions(db, now, lookback_days) - exponent, exp_fits, exp_se = _fit_current_exponent(fits) + return pinned + currents = [current for current, _, _ in points] + ambients = [ambient for _, ambient, _ in points] + rises = [rise for _, _, rise in points] + fit_n = exponent is None and max(currents) - min(currents) >= CURRENT_EXP_MIN_SPAN_A + fit_k = ambient_coef is None and max(ambients) - min(ambients) >= AMBIENT_COEF_MIN_SPREAD_C + if not fit_n and not fit_k: + return pinned + k_pinned = pinned.ambient_coef + + def solve(n: float) -> tuple[list[float], list[float], float, int] | None: + # With k pinned its contribution is moved to the left-hand side. + y = rises if fit_k else [rise - k_pinned * (ambient - AMBIENT_REF_C) for ambient, rise in zip(ambients, rises)] + design = [ + [(current / REF_CURRENT_A) ** n] + ([ambient - AMBIENT_REF_C] if fit_k else []) + for current, ambient in zip(currents, ambients) + ] + return _ols(y, design) + + if not fit_n: + result = solve(pinned.exponent) + if result is None: + return pinned + beta, se, _resid, _dof = result + return CurrentLaw(pinned.exponent, beta[1], len(points), None, se[1]) + + steps = int(round((CURRENT_EXP_RANGE[1] - CURRENT_EXP_RANGE[0]) / CURRENT_EXP_GRID_STEP)) + grid = [CURRENT_EXP_RANGE[0] + i * CURRENT_EXP_GRID_STEP for i in range(steps + 1)] + solved = [(n, solve(n)) for n in grid] + candidates = [(n, result) for n, result in solved if result is not None] + if not candidates: + return pinned + best_index = min(range(len(candidates)), key=lambda i: candidates[i][1][2]) + n, (beta, se, resid_sd, dof) = candidates[best_index] + n_se = None + if 0 < best_index < len(candidates) - 1: + # SSE(n) is locally quadratic; Var(n) = 2 sigma^2 / SSE''(n). + sse = [candidates[i][1][2] ** 2 * dof for i in (best_index - 1, best_index, best_index + 1)] + curvature = (sse[0] - 2 * sse[1] + sse[2]) / CURRENT_EXP_GRID_STEP**2 + if curvature > 0: + n_se = math.sqrt(2 * resid_sd**2 / curvature) + return CurrentLaw(n, beta[1] if fit_k else k_pinned, len(points), n_se, se[1] if fit_k else None) + + +def params_from_fits(fits: list[dict], exponent: float | None = None, + ambient_coef: float | None = None) -> ThermalParams: + """Model parameters from per-segment fits; defaults where thin. + + Fits the current law first (either term can be pinned, for comparing + laws) and re-normalizes every fit's rise_ref_c under it, in place — so + the fits handed on to the API and the degradation watch, and the + rise_ref_c median taken here, all share one law. Under the I^2, k = 0 + prior, off-reference and hot-day fits carry a bias that the watch's + current and ambient terms then have to absorb; under the fitted law + those terms have nothing left to explain.""" + law = _fit_current_law(fits, exponent, ambient_coef) for fit in fits: if fit.get("rise_c") is not None and fit.get("current_a"): - fit["rise_ref_c"] = round(fit["rise_c"] * (REF_CURRENT_A / fit["current_a"]) ** exponent, 2) + ambient = fit.get("ambient_c") + adjusted = fit["rise_c"] - (law.ambient_coef * (ambient - AMBIENT_REF_C) if ambient is not None else 0.0) + fit["rise_ref_c"] = round(adjusted * (REF_CURRENT_A / fit["current_a"]) ** law.exponent, 2) taus = [fit["tau_min"] for fit in fits] rises = [fit["rise_ref_c"] for fit in fits if fit["rise_ref_c"] is not None] rmses = [fit["rmse_c"] for fit in fits] @@ -864,12 +958,22 @@ def fit_history(db: Database, now: float, lookback_days: float = 120.0, tau_fits=len(taus), rise_fits=len(rises), fit_rmse_c=median(rmses) if rmses else None, - current_exp=exponent, - current_exp_fits=exp_fits, - current_exp_se=exp_se, + current_exp=law.exponent, + current_exp_fits=law.fits, + current_exp_se=law.exponent_se, + ambient_coef=law.ambient_coef, + ambient_coef_se=law.ambient_coef_se, ) +def fit_history(db: Database, now: float, lookback_days: float = 120.0, + fits: list[dict] | None = None) -> ThermalParams: + """Aggregate per-session fits into model parameters; defaults where thin.""" + if fits is None: + fits = fit_sessions(db, now, lookback_days) + return params_from_fits(fits) + + # ---------------- degradation watch ---------------- # A loose lug or degrading contact shows up as extra resistance: more heat @@ -1255,7 +1359,7 @@ def sustainable_max_current(ambient_c: float, params: ThermalParams) -> float | trim the last amp or two. Worked from a measured ambient when a sensor is reporting — see predict() for why the trajectory-implied one is only the fallback.""" - headroom = TRIP_HANDLE_C - SUGGEST_MARGIN_C - ambient_c + headroom = TRIP_HANDLE_C - SUGGEST_MARGIN_C - ambient_c - params.ambient_coef * (ambient_c - AMBIENT_REF_C) if headroom <= 0: return None amps = math.floor(params.current_for_rise(headroom)) @@ -1346,7 +1450,7 @@ def _recent_steady_ambient(recent: list[dict], params: ThermalParams) -> float | window = [(sample["ts"], sample["handle_temp_c"]) for sample in run] t_inf, _se = _project_t_inf(window, params.tau_min) run_current = median(sample["vehicle_current_a"] for sample in run) - ambient = t_inf - params.rise_at(run_current) + ambient = params.ambient_from_plateau(t_inf, run_current) if -30.0 <= ambient <= TRIP_HANDLE_C: return ambient return None @@ -1427,7 +1531,7 @@ def predict(db: Database, now: float, params: ThermalParams) -> dict: gap_reason = "warming_up" out["forecast"] = {"basis": "insufficient", "will_trip": None, "reason": gap_reason} return out - t_inf = ambient + params.rise_at(current) + t_inf = params.plateau_at(current, ambient) forecast["basis"] = "model" forecast["ambient_source"] = source # A model-basis plateau is only as good as its ambient: a sensor @@ -1447,7 +1551,7 @@ def predict(db: Database, now: float, params: ThermalParams) -> dict: # reporting — interpolation inside the fitted data rather than a # round trip to 32 A and back (on one install, after a 32 A probe, # the implied route named 47 A; the measured one 44 A, which held). - implied_ambient = t_inf - params.rise_at(current) + implied_ambient = params.ambient_from_plateau(t_inf, current) measured = _latest_measured_ambient(db, now) sustain_ambient, sustain_source = measured if measured is not None else (implied_ambient, "implied") forecast.update( @@ -1485,7 +1589,7 @@ def predict(db: Database, now: float, params: ThermalParams) -> dict: out["ambient_se_c"] = 0.3 if measured is not None else idle_model.ambient_se_c out["ambient_stable"] = stable # Hypothetical: a full-rate session started right now. - t_inf = ambient + params.rise_ref_c + t_inf = params.plateau_at(REF_CURRENT_A, ambient) minutes = _minutes_to_trip(last["handle_temp_c"], t_inf, tau_min) out["forecast"] = { "basis": "hypothetical", @@ -1493,7 +1597,7 @@ def predict(db: Database, now: float, params: ThermalParams) -> dict: "will_trip": minutes is not None, "minutes_to_trip": round(minutes, 1) if minutes is not None else None, "trip_ts": None, - "safe_ambient_max_c": round(TRIP_HANDLE_C - params.rise_ref_c, 1), + "safe_ambient_max_c": round(params.safe_ambient_max_c(), 1), "suggested_max_a": suggest_max_current(ambient, params) if minutes is not None else None, "sustainable_max_a": sustainable_max_current(ambient, params), }