diff --git a/contrib/backtest_forecast.py b/contrib/backtest_forecast.py new file mode 100644 index 0000000..61b4cf0 --- /dev/null +++ b/contrib/backtest_forecast.py @@ -0,0 +1,530 @@ +#!/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 (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. + +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)) +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 +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 + 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: + 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 + se_c: float | None = None # set when actual_c is a projection rather than an observation + + +# --------------------------------------------------------------------------- + + +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) + 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) + 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 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.""" + 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], + ) + 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: + 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 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 + 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, actual, predictions, se, + ) + + +# --------------------------------------------------------------------------- + + +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 _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, + 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(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), + 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") + 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() + + +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] = [] + trip_predictions: dict[tuple[int, int], dict[str, float]] = {} + 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): + prev = session_runs[index - 1] if index else None + ticks.extend(score_ticks(run, laws["n"])) + 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( + 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)} 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, trip_predictions) + 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..8574eb5 100644 --- a/docs/thermal-model.md +++ b/docs/thermal-model.md @@ -160,6 +160,71 @@ 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. + +**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 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]