From f77a56755724e4c161d86a7055a8cda2bf9d4fe8 Mon Sep 17 00:00:00 2001 From: Fernando Gonzalez Date: Sat, 3 Oct 2026 17:50:38 -0400 Subject: [PATCH] feat: a settle-on-target amp controller and a closed-loop backtest to judge it (experimental, not wired in) 2026-10-02 (session 101): the amp controller moved 49 -> 40 -> 46 -> 44 -> 46 -> 48 -> 45 -> 41 A in eighteen minutes, then the quick-reversal back-off pinned the session at 41 A when ~44 A was sustainable. Step-ups are decided on the plateau at the present current, never the target's; caps trust a two-minute-old trajectory fit; and a 2-3 A step stays inside predict()'s 10% window band, so the next fit mixes two currents. contrib/settle_amp_control.py is a candidate replacement for the derate decision. It tracks one number per session - the effective ambient, plateau minus the model's rise - fused from the ambient prior and one reading per steady-current stretch, and targets the largest whole amp whose plateau stays a soak margin under the trip. It front-loads at full rate while the handle is cold, makes one planned drop, and afterwards steps down on whichever is higher of the fused estimate and the live stretch's own reading, with a handle-driven cut near the trip that ignores the estimate. Not wired into the daemon. contrib/backtest_settle_amp_control.py closes the loop the existing backtest leaves open: it replays a session's real samples through the real predict() until the simulated command diverges from what the car drew, then hands over to a two-time-constant plant fitted to that session, and reruns each session against plants 1-2 C hotter than the fit. Whole history, 82 sessions with a usable plant (15 short, 4 poor fits skipped), fitted plant, deployed vs settle at a 1.5 C soak margin: trips 0 / 0, worst peak 64.5 / 63.7 C, amp changes 432 / 288, 3542 / 3561 Ah. On the 54 sessions where full rate is not sustainable: median changes 7 / 5, mean reversals 3.1 / 0.9, sessions with no reversal 3 / 31. No trips for either on any of the four hotter plants. A 2.0 C margin delivers less than the deployed controller and is no safer. Two drafts failed this backtest before this one passed, both pinned by tests: counting every tick of a long stretch as a new reading (held 48 A into a 64.3 C peak on 09-30), and averaging the live reading with earlier, cooler ones (seven trips on hot days across the history). Known limits: the plant is fitted to the day it is judged on; the forecast still under-reads the plateau until the cable has soaked, which is why hot days take several 1 A steps down rather than one move; 09-06 is worse than the deployed controller (12 changes, 6 reversals against 10 and 4). Co-Authored-By: Claude Opus 5.5 --- contrib/backtest_settle_amp_control.py | 501 +++++++++++++++++++++++++ contrib/settle_amp_control.py | 249 ++++++++++++ tests/test_settle_amp_control.py | 134 +++++++ 3 files changed, 884 insertions(+) create mode 100644 contrib/backtest_settle_amp_control.py create mode 100644 contrib/settle_amp_control.py create mode 100644 tests/test_settle_amp_control.py diff --git a/contrib/backtest_settle_amp_control.py b/contrib/backtest_settle_amp_control.py new file mode 100644 index 0000000..8b6bc86 --- /dev/null +++ b/contrib/backtest_settle_amp_control.py @@ -0,0 +1,501 @@ +#!/usr/bin/env python3 +"""Closed-loop backtest of an amp controller: what would the handle have +done, session by session, had a different controller been setting the +current? + +backtest_derate_amp_control.py replays recorded temperatures, so it can say +when a controller *would have acted* but not what would have happened next +— the recorded handle followed the amps the deployed controller chose. This +script closes the loop: + +1. **Replay until divergence.** The session's real samples are fed to the + real ``thermal.predict()`` for as long as the simulated controller's + command matches what the car really drew. The first minutes — the + S-shaped start and the noise that fool a young trajectory fit — are the + recorded ones, not a model's idea of them. +2. **Plant after.** From the first tick the commands differ, a per-session + plant takes over: two time constants (handle, and the slower cable soak) + plus a short sensor lag, fitted to that session's own recorded current + and handle temperature. It has been run alongside the real currents up to + that point, so its state is warm, and the real-minus-plant offset at the + handover decays over the fast time constant. +3. **Model as of the session's start.** ``fit_history`` is frozen just + before the session began, as in backtest_derate_amp_control.py, so a + session never leaks into its own forecast. + +A plant fitted to the day it is judged on flatters every controller, so +each session is also run against **hotter plants** (--stress): the ambient +1 and 2 C above the fitted one, and the cable soaking 1 and 2 C more at +full rate. The perturbation rides on the replayed samples too. + +Scored per controller: charger trips (must be zero), peak handle +temperature, amp changes, amp-hours delivered in the session's real +charging time. Sessions whose plant fits poorly are counted and left out +of the totals rather than trusted. + +Read-only: never writes to --db. Point it at a copy. + +Example: + ./backtest_settle_amp_control.py --db ../wallmonitor-copy.db --stress +""" + +from __future__ import annotations + +import argparse +import datetime +import math +import sqlite3 +import sys +from dataclasses import dataclass, field, replace +from multiprocessing import Pool + +import derate_amp_control +import settle_amp_control + +from wallmonitor import thermal +from wallmonitor.db import Database + +SAMPLE_S = 2.0 # the poller's cadence +TICK_S = 30.0 # the daemon's +FULL_A = 48.7 # what the car draws at an uncapped 48 A +CMD_OFFSET_A = 0.6 # a 45 A cap draws ~45.6 A +DERATE_A = 24.0 # the charger's 50% fold-back, held for the rest of the session +POOR_FIT_RMSE_C = 0.6 +STRESS_CASES = (("ambient", 1.0), ("ambient", 2.0), ("soak", 1.0), ("soak", 2.0)) + + +# --------------------------------------------------------------------------- +# plant + + +@dataclass +class Plant: + """T = a + y1 + y2, where x1 and x2 relax toward g*r(I) and h*r(I) with + time constants tau1 and tau2 (minutes), r(I) = (I/48)^n, and y trails x + by `lag` minutes (the reading trails the heat).""" + + a: float + g: float + h: float + tau1: float + tau2: float + n: float + lag: float + rmse: float + x1: float = 0.0 + x2: float = 0.0 + y1: float = 0.0 + y2: float = 0.0 + + def r(self, amps: float) -> float: + return (max(amps, 0.0) / thermal.REF_CURRENT_A) ** self.n + + def temp(self) -> float: + return self.a + (self.y1 + self.y2 if self.lag else self.x1 + self.x2) + + def step(self, amps: float, dt_s: float) -> float: + dt = dt_s / 60.0 + drive = self.r(amps) + self.x1 += (self.g * drive - self.x1) * (1 - math.exp(-dt / self.tau1)) + if self.h: + self.x2 += (self.h * drive - self.x2) * (1 - math.exp(-dt / self.tau2)) + if self.lag: + blend = 1 - math.exp(-dt / self.lag) + self.y1 += (self.x1 - self.y1) * blend + self.y2 += (self.x2 - self.y2) * blend + return self.temp() + + def steady(self, amps: float) -> float: + return self.a + (self.g + self.h) * self.r(amps) + + +def charging_run(db_path: str, session_id: int, max_gap_s: float = 60.0) -> list[tuple[float, float, float]]: + """(ts, amps, handle_c) over the session's first unbroken charging run.""" + con = sqlite3.connect(db_path) + rows = con.execute( + "SELECT ts, vehicle_current_a, handle_temp_c FROM vitals_samples " + "WHERE session_id = ? AND contactor_closed = 1 AND vehicle_current_a >= 6 " + "AND handle_temp_c IS NOT NULL AND handle_temp_c < 255 ORDER BY ts", + (session_id,), + ).fetchall() + con.close() + run: list[tuple[float, float, float]] = [] + for row in rows: + if run and row[0] - run[-1][0] > max_gap_s: + break + run.append(row) + return run + + +def _unit_responses(run, tau1: float, tau2: float | None, n: float, lag: float) -> list[list[float]]: + """The design matrix rows [1, fast response, slow response] to the run's + own current, for unit gains — the plant is linear in (a, g, h).""" + x1 = x2 = y1 = y2 = 0.0 + prev = run[0][0] + rows = [] + for ts, amps, _temp in run: + dt = (ts - prev) / 60.0 + prev = ts + drive = (amps / thermal.REF_CURRENT_A) ** n + x1 += (drive - x1) * (1 - math.exp(-dt / tau1)) + if tau2: + x2 += (drive - x2) * (1 - math.exp(-dt / tau2)) + if lag: + blend = 1 - math.exp(-dt / lag) if dt > 0 else 0.0 + y1 += (x1 - y1) * blend + y2 += (x2 - y2) * blend + fast, slow = (y1, y2) if lag else (x1, x2) + rows.append([1.0, fast, slow] if tau2 else [1.0, fast]) + return rows + + +def fit_plant(run, n: float, two_tau: bool = True) -> Plant | None: + """Grid over the time constants and lag, OLS for the gains inside.""" + temps = [temp for _, _, temp in run] + best: tuple | None = None + for tau1 in [x / 2 for x in range(6, 41)]: # 3..20 min + for tau2 in ((12, 15, 20, 25, 30, 40, 50, 60, 90) if two_tau else (None,)): + if tau2 is not None and tau2 < 1.5 * tau1: + continue + for lag in (0.0, 0.25, 0.5, 1.0): + solved = thermal._ols(temps, _unit_responses(run, tau1, tau2, n, lag)) + if solved is None: + continue + coef, _se, rmse, _dof = solved + if any(gain < 0 for gain in coef[1:]): + continue + if best is None or rmse < best[0]: + best = (rmse, tau1, tau2, lag, coef) + if best is None: + return None + rmse, tau1, tau2, lag, coef = best + return Plant(a=coef[0], g=coef[1], h=coef[2] if tau2 else 0.0, tau1=tau1, tau2=tau2 or 1.0, n=n, lag=lag, rmse=rmse) + + +def fit_session_plant(run, n: float) -> tuple[Plant | None, str]: + """Two time constants when the session's data can tell them apart; + one when it cannot (a single current held throughout leaves the slow + gain free to run away); None when neither fits.""" + plant = fit_plant(run, n, two_tau=True) + if plant is not None and plant.rmse <= POOR_FIT_RMSE_C and plant.steady(FULL_A) - plant.a <= 60.0: + return plant, "two-tau" + plant = fit_plant(run, n, two_tau=False) + if plant is not None and plant.rmse <= POOR_FIT_RMSE_C and plant.steady(FULL_A) - plant.a <= 60.0: + return plant, "one-tau" + return None, "poor fit" + + +def perturbed(plant: Plant, kind: str, delta_c: float) -> Plant: + if kind == "ambient": # the day runs hotter than the fit (and the sensor) says + return replace(plant, a=plant.a + delta_c) + if kind == "soak": # the cable soaks delta_c more at full rate + return replace(plant, h=plant.h + delta_c / plant.r(FULL_A), tau2=plant.tau2 if plant.h else 20.0) + return replace(plant) + + +# --------------------------------------------------------------------------- +# closed loop + + +class SimDB(Database): + """The real DB, except vitals from sim_start on are the simulated ones.""" + + def __init__(self, path: str, sim_start: float): + super().__init__(path) + self.sim_start = sim_start + self.synth: list[dict] = [] + + def vitals_range(self, t_from, t_to, max_points=1500): + real = [] + if t_from < self.sim_start: + real = [ + row + for row in super().vitals_range(t_from, min(t_to, self.sim_start - 0.001), max_points) + if row["ts"] < self.sim_start + ] + return real + [row for row in self.synth if t_from <= row["ts"] <= t_to] + + +def drawn(cmd: float | None) -> float: + return FULL_A if cmd is None else min(FULL_A, cmd + CMD_OFFSET_A) + + +def settle_forecast(result: dict, db: SimDB, now: float, params: thermal.ThermalParams, session_start: float) -> dict: + """The forecast fields settle_amp_control wants and predict() does not + serve yet: the prior ambient, the trajectory window's current, and the + plateau projected over a tight current band (3% rather than 10%), so a + 2-3 A step starts a fresh window instead of mixing two currents.""" + forecast = result.get("forecast") + if result.get("state") != "charging" or not forecast: + return result + idle_model = thermal.load_idle_offset(db) + measured = thermal._latest_measured_ambient(db, now) + if measured is not None: + forecast["prior_ambient_c"], forecast["prior_ambient_sd_c"] = measured[0], 0.3 + else: + forecast["prior_ambient_c"] = thermal._ambient_before(db, session_start, idle_model) + forecast["prior_ambient_sd_c"] = idle_model.ambient_se_c + recent = [row for row in db.vitals_range(now - 900, now, 2000) if row.get("handle_temp_c") is not None] + last = recent[-1] + current = last["vehicle_current_a"] + band = max(1.0, 0.03 * current) + window = [] + for sample in reversed(recent): + if not sample.get("contactor_closed") or last["ts"] - sample["ts"] > 360: + break + if abs(sample["vehicle_current_a"] - current) > band: + break + window.append(sample) + window.reverse() + forecast["window_current_a"] = sorted(s["vehicle_current_a"] for s in window)[len(window) // 2] + points = [(s["ts"], s["handle_temp_c"]) for s in window] + if len(points) >= thermal.TRAJECTORY_MIN_SAMPLES and points[-1][0] - points[0][0] >= thermal.TRAJECTORY_MIN_SPAN_S: + t_inf, se = thermal._project_t_inf(points, params.tau_min) + forecast["basis"] = "trajectory" + forecast["steady_state_c"] = round(t_inf, 1) + forecast["steady_state_se_c"] = round(max(se, 0.1), 2) if se is not None else None + elif forecast.get("basis") == "trajectory": + forecast["basis"] = "model" # the tight band leaves too short a window for a slope + return result + + +@dataclass +class Controller: + name: str + decide: object + new_state: object + cfg: object + settle_fields: bool + + +def controllers(soak_margins: list[float], unguarded: bool = False) -> list[Controller]: + out = [ + Controller( + "deployed", + derate_amp_control.decide, + lambda: derate_amp_control.State(last_session_state="idle"), + derate_amp_control.Config(forecast_confidence_k=4.0, probe_amps=()), + False, + ) + ] + for margin in soak_margins: + out.append( + Controller( + f"settle m{margin:g}", + settle_amp_control.decide, + lambda: settle_amp_control.State(last_session_state="idle"), + settle_amp_control.Config(soak_margin_c=margin), + True, + ) + ) + if unguarded: + # As first reviewed: no young-stretch guard, "rising" on one tick. + out.append( + Controller( + "settle raw", + settle_amp_control.decide, + lambda: settle_amp_control.State(last_session_state="idle"), + settle_amp_control.Config(young_guard_tau=0.0, rising_window_s=0.0, rising_min_c=0.05), + True, + ) + ) + return out + + +@dataclass +class Result: + trip_min: float | None = None + peak_c: float = 0.0 + amp_hours: float = 0.0 + final_a: float = 0.0 + moves: list[tuple[float, float, str]] = field(default_factory=list) # (minute, to amps, why) + + +def simulate(db_path, session_id, run, plant: Plant, base: Plant, params, ctl: Controller, session_start: float) -> Result: + """One session under one controller. `plant` may be a perturbed copy of + `base` (the fitted one); while replaying, the perturbation rides on the + real samples as (plant - base).""" + t0, t_end = run[0][0], run[-1][0] + db = SimDB(db_path, t0) + plant, base = replace(plant), replace(base) + state = ctl.new_state() + cmd: float | None = None + replaying, derated = True, False + offset = 0.0 + index = 0 + out = Result() + t = prev_t = t0 + next_tick = t0 + TICK_S + while t <= t_end: + if replaying and index < len(run): + ts, real_a, real_temp = run[index] + settled = abs(run[index - 1][1] - real_a) < 1.0 + if ts - t0 > 30 and settled and abs(drawn(cmd) - real_a) > 1.5: + replaying = False + offset = real_temp - base.temp() + if replaying and index < len(run): + t, amps, temp = run[index] + index += 1 + plant.step(amps, t - prev_t) + base.step(amps, t - prev_t) + temp += plant.temp() - base.temp() + else: + replaying = False + amps = DERATE_A if derated else drawn(cmd) + offset *= math.exp(-(t - prev_t) / (plant.tau1 * 60.0)) + temp = plant.step(amps, t - prev_t) + offset + if temp >= thermal.TRIP_HANDLE_C and not derated: + derated = True + out.trip_min = (t - t0) / 60.0 + db.synth.append( + {"ts": t, "session_id": session_id, "vehicle_current_a": round(amps, 1), "handle_temp_c": round(temp, 1), + "contactor_closed": 1, "vehicle_connected": 1} + ) + out.peak_c = max(out.peak_c, temp) + out.amp_hours += amps * (t - prev_t) / 3600.0 + out.final_a = amps + prev_t = t + if t >= next_tick: + next_tick += TICK_S + result = thermal.predict(db, t, params) + if ctl.settle_fields: + settle_forecast(result, db, t, params, session_start) + action, state, reason = ctl.decide(result, state, ctl.cfg) + new_cmd = float(action.value) if action.kind == "cap" else None if action.kind == "restore" else cmd + if new_cmd != cmd: + out.moves.append(((t - t0) / 60.0, new_cmd if new_cmd is not None else 48.0, reason)) + cmd = new_cmd + t = run[index][0] if replaying and index < len(run) else t + SAMPLE_S + db.close() + return out + + +def best_constant(plant: Plant, duration_s: float) -> int: + """Highest whole-amp command held from the start that never trips.""" + best = 6 + for cmd in range(6, 49): + trial = replace(plant) + amps = drawn(None if cmd >= 48 else cmd) + elapsed = 0.0 + while elapsed <= duration_s and trial.step(amps, SAMPLE_S) < thermal.TRIP_HANDLE_C: + elapsed += SAMPLE_S + if elapsed > duration_s: + best = cmd + return best + + +# --------------------------------------------------------------------------- +# driver + + +def run_session(job: tuple) -> dict: + db_path, session_id, soak_margins, stress, min_run_min, unguarded = job + db = Database(db_path) + session = db.session(session_id) + run = charging_run(db_path, session_id) + info = {"id": session_id, "start": session["start_ts"], "cases": {}} + if len(run) < 10 or run[-1][0] - run[0][0] < min_run_min * 60.0: + db.close() + return {**info, "skip": "short"} + start = session["start_ts"] + params = thermal.fit_history(db, start - 1, fits=thermal.fit_sessions(db, start - 1)) + ambient = thermal._latest_measured_ambient(db, run[0][0]) + real_trip = any(str(a["alert"]) == "40" for a in db.alerts_range(run[0][0], run[-1][0])) + real_moves = len(db.events_range(run[0][0], run[-1][0], kinds=["amp_capped", "amp_restored"])) + db.close() + plant, fit_kind = fit_session_plant(run, params.current_exp) + info.update( + minutes=(run[-1][0] - run[0][0]) / 60.0, ambient=ambient[0] if ambient else None, + real_peak=max(temp for *_, temp in run), real_trip=real_trip, real_moves=real_moves, fit=fit_kind, + ) + if plant is None: + return {**info, "skip": "poor plant fit"} + info.update(rmse=plant.rmse, steady_full=plant.steady(FULL_A)) + cases = [("nominal", 0.0)] + (list(STRESS_CASES) if stress else []) + for kind, delta in cases: + hot = perturbed(plant, kind, delta) + entry = {"best": best_constant(hot, run[-1][0] - run[0][0]), "ctl": {}} + for ctl in controllers(soak_margins, unguarded): + entry["ctl"][ctl.name] = simulate(db_path, session_id, run, hot, plant, params, ctl, start) + info["cases"][f"{kind}+{delta:g}"] = entry + return info + + +def main(argv: list[str] | None = None) -> int: + parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0]) + parser.add_argument("--db", required=True, help="path to a wallmonitor.db (read-only; use a copy)") + parser.add_argument("--sessions", help="comma-separated session ids (default: every session)") + parser.add_argument("--soak-margin", default="1.5", help="settle soak margins to score, comma-separated") + parser.add_argument("--stress", action="store_true", help="also run each session against hotter plants") + parser.add_argument("--unguarded", action="store_true", help="also score settle without its review guards") + parser.add_argument("--min-run-min", type=float, default=15.0, help="skip shorter charging runs") + parser.add_argument("--jobs", type=int, default=8) + parser.add_argument("-v", "--verbose", action="store_true", help="print every amp change") + args = parser.parse_args(argv) + + margins = [float(x) for x in args.soak_margin.split(",")] + if args.sessions: + ids = [int(x) for x in args.sessions.split(",")] + else: + con = sqlite3.connect(args.db) + ids = [row[0] for row in con.execute("SELECT id FROM sessions ORDER BY id")] + con.close() + jobs = [(args.db, sid, margins, args.stress, args.min_run_min, args.unguarded) for sid in ids] + with Pool(args.jobs) as pool: + results = pool.map(run_session, jobs) + + totals: dict[tuple[str, str, str], dict] = {} + skipped: dict[str, int] = {} + for info in results: + day = datetime.datetime.fromtimestamp(info["start"]).strftime("%a %m-%d") + if "skip" in info: + skipped[info["skip"]] = skipped.get(info["skip"], 0) + 1 + if info["skip"] != "short": + print(f"session {info['id']} {day}: skipped ({info['skip']})") + continue + ambient = f"{info['ambient']:.1f}C" if info["ambient"] is not None else "no sensor" + print( + f"session {info['id']} {day}: {info['minutes']:.0f} min, ambient {ambient}, real peak " + f"{info['real_peak']:.1f}C, real amp changes {info['real_moves']}" + f"{', REAL TRIP' if info['real_trip'] else ''}; plant {info['fit']} rmse {info['rmse']:.2f}C, " + f"plateau at full rate {info['steady_full']:.1f}C" + ) + for case, entry in info["cases"].items(): + print(f" {case}: best constant {entry['best']} A") + for name, res in entry["ctl"].items(): + trip = f"TRIP at {res.trip_min:.0f}min" if res.trip_min is not None else "no trip" + path = " ".join(f"{minute:.0f}m:{amps:g}" for minute, amps, _ in res.moves) + print( + f" {name:<12} {trip}, peak {res.peak_c:.1f}C, {len(res.moves)} changes, " + f"{res.amp_hours:.1f} Ah, ends {res.final_a:.1f}A [{path}]" + ) + if args.verbose: + for minute, amps, why in res.moves: + print(f" {minute:5.1f}m -> {amps:g}A {why}") + # Only a session whose full rate is not sustainable can tell + # controllers apart; the rest check that nothing caps needlessly. + limited = info["cases"]["nominal+0"]["best"] < 48 + for group in ("all", "limited" if limited else "full-rate ok"): + agg = totals.setdefault( + (group, case, name), {"n": 0, "trips": 0, "moves": 0, "ah": 0.0, "peak": 0.0, "moved": 0} + ) + agg["n"] += 1 + agg["trips"] += res.trip_min is not None + agg["moves"] += len(res.moves) + agg["moved"] += bool(res.moves) + agg["ah"] += res.amp_hours + agg["peak"] = max(agg["peak"], res.peak_c) + + print("\nTOTALS (sessions with a usable plant)") + for (group, case, name), agg in sorted(totals.items(), key=lambda item: item[0][0]): + print( + f" {group:<12} {case:<10} {name:<12} sessions {agg['n']}, trips {agg['trips']}, amp changes {agg['moves']} " + f"(in {agg['moved']} sessions), {agg['ah']:.1f} Ah, worst peak {agg['peak']:.1f}C" + ) + if skipped: + print("skipped: " + ", ".join(f"{count} {why}" for why, count in skipped.items())) + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/contrib/settle_amp_control.py b/contrib/settle_amp_control.py new file mode 100644 index 0000000..692c532 --- /dev/null +++ b/contrib/settle_amp_control.py @@ -0,0 +1,249 @@ +#!/usr/bin/env python3 +"""Candidate replacement for derate_amp_control's derate decision: settle on +the highest sustainable charge current in one or two moves instead of +hunting for it. + +**Experimental — not wired into the daemon.** It exists so the closed-loop +backtest (backtest_settle_amp_control.py) can score it against the deployed +``decide()`` on the whole history before anything is deployed. + +Why: derate_amp_control's restore path decides to step up from the plateau +at the *present* current plus the handle's instantaneous margin, never from +the plateau at the *target* current, and its cap path trusts a +two-minute-old trajectory fit. 2026-10-02 (session 101): 49 -> 40 -> 46 -> +44 -> 46 -> 48 -> 45 -> 41 A in eighteen minutes, then the quick-reversal +back-off pinned the session at 41 A when ~44 A was sustainable. + +The idea is one number tracked per session — the **effective ambient** +``theta = plateau - rise(I)``, what the handle would settle at with no +current — and one target derived from it: + + I* = largest whole amp with theta + k * sd + rise(I) <= trip - soak_margin + +``theta`` is an inverse-variance fusion of + +- a prior: the ambient sensor (or the pre-session idle-handle proxy), and +- **one reading per steady-current stretch**: that stretch's projected + plateau minus the model's rise at its current. A stretch is a single + reading however long it runs — consecutive ticks are the same data, and + counting them as independent lets stale early readings outvote the live + one (that bug, in a draft of this file, held 48 A into a 64.3 C peak). A + young stretch is distrusted explicitly: its sd is floored at + ``floor_sd_c + young_sd_c * exp(-age / tau)``, because the projection's + own standard error measures residual noise, not how wrong a single + exponential is two minutes after a current change. + +The moves: + +- **Front-load.** The first minutes at full rate are thermally free (the + handle is 20-30 C under the trip) and are the best reading of the day + there is, taken at the reference current with no exponent extrapolation. + Hold full rate until the handle is within ``frontload_c`` of the target + plateau, then make one planned drop to I* (at ``k_down``: this is a first + landing, not a step-up, and landing 1 A high is corrected by an ordinary + 1 A step down). +- **Down** whenever I* (at ``k_down``) is at least 1 A under the present + current. Not rate-limited, and judged on the higher of the fused theta + and the present stretch's own reading: the cable soaks through a hot + session, so earlier readings are of a cooler cable. +- **Up** only when I* (at ``k_up``) is at least 1 A over, the present + stretch is at least a time constant old, the handle has margin, and a + time constant has passed since the last move. A step-up undone within a + time constant widens ``k_up`` for the rest of the session — it never + gives up and pins the session low. +- **Emergency:** handle at ``emergency_c`` and still rising (against itself + two minutes ago, not the last tick) — cut 2 A below the present current, + or below I* if that is lower, whatever the estimate says. + +``decide()`` takes the same ``/api/thermal`` dict as the daemon's, plus +three forecast fields the server does not serve yet (the backtest computes +them): ``window_current_a``, ``prior_ambient_c``/``prior_ambient_sd_c``, and +a ``steady_state_c`` projected over a tight current band (3% of the +current, not 10%) so a 2-3 A step starts a fresh window instead of mixing +two currents in one fit. +""" + +from __future__ import annotations + +import math +from dataclasses import dataclass, field, replace + +from derate_amp_control import TRIP_HANDLE_C, Action + + +@dataclass +class Config: + normal_amps: float = 48.0 + min_amps: float = 16.0 + k_up: float = 1.5 + k_down: float = 0.5 + soak_margin_c: float = 1.5 # the cable keeps soaking after the handle looks settled + prior_sd_c: float = 2.5 # floor on the prior's sd, whatever its source claims + floor_sd_c: float = 1.0 # sd of a mature stretch's reading + young_sd_c: float = 4.0 # extra sd of a brand-new stretch, decaying over tau + drift_c_per_min: float = 0.025 # a finished stretch's reading ages as the cable soaks + frontload_c: float | None = 4.0 # None: land on I* at the first tick instead + up_handle_max_c: float = 63.5 + emergency_c: float = 63.5 + reversal_k_step: float = 0.5 + # A stretch younger than this many time constants reads low after a + # step down (the fixed-tau projection on a still-cooling cable: session + # 101 read 58.4 C at 24 min against 60.8 C settled), so it may raise the + # estimate but not pull it under the last long stretch's. 0 disables. + young_guard_tau: float = 2.0 + # "Rising" is the handle against itself this long ago, by at least + # rising_min_c: one tick of a 0.1 C sensor is a coin flip at the plateau. + live_down_tau: float = 0.5 # stretch age (in tau) from which its own reading can force a step down + rising_window_s: float = 120.0 + rising_min_c: float = 0.3 + + +@dataclass +class State: + last_session_state: str | None = None + cap: float | None = None # None: uncapped, drawing normal_amps + frontloading: bool = False + stretch_start_ts: float | None = None + stretch_current_a: float | None = None + live: tuple[float, float] | None = None # (theta, sd) of the present stretch + done: list[tuple[float, float, float]] = field(default_factory=list) # finished stretches: (ts, theta, sd) + last_move_ts: float | None = None + last_up_ts: float | None = None + k_extra: float = 0.0 + handle_log: list[tuple[float, float]] = field(default_factory=list) # recent (ts, handle_c) + anchor_theta_c: float | None = None # theta of the last stretch that ran long enough to trust + theta_c: float | None = None # the fused estimate, kept for the log + + +def _rise(model: dict, amps: float) -> float: + return model["rise_ref_c"] * (amps / model["ref_current_a"]) ** model["current_exp"] + + +def _amps_for(theta: float, headroom_c: float, model: dict, cfg: Config) -> float: + """Largest whole amp whose plateau sits headroom_c + soak_margin_c under the trip.""" + budget = model.get("trip_c", TRIP_HANDLE_C) - cfg.soak_margin_c - theta - headroom_c + if budget <= 0: + return cfg.min_amps + amps = model["ref_current_a"] * (budget / model["rise_ref_c"]) ** (1.0 / model["current_exp"]) + return float(max(cfg.min_amps, min(cfg.normal_amps, math.floor(amps)))) + + +def decide(thermal: dict, state: State, cfg: Config) -> tuple[Action, State, str]: + """What to do, the state to persist, and why.""" + forecast = thermal.get("forecast") or {} + model = thermal.get("model") or {} + now = thermal.get("ts") + handle = thermal.get("handle_c") + + if thermal.get("state") != "charging": + new = State(last_session_state=thermal.get("state")) + if state.cap is not None: + return Action("restore", cfg.normal_amps), new, "session ended: restoring normal rate" + return Action("none"), new, "idle: nothing to do" + + session_started = state.last_session_state != "charging" + st = State(last_session_state="charging") if session_started else replace(state) + tau_s = model["tau_min"] * 60.0 + trip_c = model.get("trip_c", TRIP_HANDLE_C) + + # One steady-current stretch is one reading. Tracked here rather than + # read off the server's sliding sample buffer, whose window start moves + # every tick once a stretch outlasts it. + current = forecast.get("window_current_a") or thermal.get("current_a") + if st.stretch_current_a is None or abs(current - st.stretch_current_a) > max(1.0, 0.03 * current): + if st.live is not None: + st.done = st.done + [(now, *st.live)] + if cfg.young_guard_tau and now - st.stretch_start_ts >= cfg.young_guard_tau * tau_s: + st.anchor_theta_c = st.live[0] + st.stretch_start_ts, st.stretch_current_a, st.live = now, current, None + age_s = now - st.stretch_start_ts + if forecast.get("basis") == "trajectory" and isinstance(forecast.get("steady_state_c"), (int, float)): + se = forecast.get("steady_state_se_c") or 0.1 + sd = max(se, cfg.floor_sd_c + cfg.young_sd_c * math.exp(-age_s / tau_s)) + st.live = (forecast["steady_state_c"] - _rise(model, current), sd) + + readings: list[tuple[float, float]] = [] + prior = forecast.get("prior_ambient_c") + if isinstance(prior, (int, float)): + readings.append((prior, max(cfg.prior_sd_c, forecast.get("prior_ambient_sd_c") or 0.0))) + for ts, theta, sd in st.done: + readings.append((theta, math.hypot(sd, cfg.drift_c_per_min * (now - ts) / 60.0))) + if st.live is not None: + readings.append(st.live) + if not readings: + st.frontloading = st.frontloading or session_started + return Action("none"), st, "no ambient and no trajectory yet: holding" + weight = sum(1.0 / sd**2 for _, sd in readings) + theta = sum(value / sd**2 for value, sd in readings) / weight + sd = math.sqrt(1.0 / weight) + if st.anchor_theta_c is not None and age_s < cfg.young_guard_tau * tau_s: + theta = max(theta, st.anchor_theta_c) + st.theta_c = theta + + # The fusion treats theta as a constant; it is not — the cable soaks + # through a hot session and every earlier reading is of a cooler one. + # Averaging them in is optimistic exactly when it matters, so stepping + # *down* is decided on whichever is higher: the fused estimate or the + # present stretch's own reading once it is half a time constant old. + # (Whole-history backtest: without this, hot days landed 3-4 A high and + # walked down 1 A at a time into the trip — seven sessions.) + theta_down = theta + if st.live is not None and age_s >= cfg.live_down_tau * tau_s: + theta_down = max(theta, st.live[0]) + + cur = st.cap if st.cap is not None else cfg.normal_amps + k_up = cfg.k_up + st.k_extra + i_down = _amps_for(theta_down, cfg.k_down * sd, model, cfg) + i_up = _amps_for(theta, k_up * sd, model, cfg) + rising = False + if isinstance(handle, (int, float)): + log = [(ts, value) for ts, value in st.handle_log if now - ts <= cfg.rising_window_s + 45.0] + [(now, handle)] + st.handle_log = log + earlier = [value for ts, value in log if now - ts >= cfg.rising_window_s] + rising = bool(earlier) and handle - earlier[-1] >= cfg.rising_min_c + tag = f"theta {theta:.1f}+-{sd:.1f}C from {len(readings)} readings, I* {i_down:g}A down / {i_up:g}A up" + + def move(target: float, why: str) -> tuple[Action, State, str]: + going_up = target > cur + undone = not going_up and st.last_up_ts is not None and now - st.last_up_ts <= tau_s + new = replace( + st, + cap=None if target >= cfg.normal_amps else float(target), + last_move_ts=now, + last_up_ts=now if going_up else st.last_up_ts, + k_extra=st.k_extra + (cfg.reversal_k_step if undone else 0.0), + ) + if target >= cfg.normal_amps: + return Action("restore", cfg.normal_amps), new, f"{why}: {cur:g} -> {cfg.normal_amps:g}A ({tag})" + return Action("cap", float(target)), new, f"{why}: {cur:g} -> {target:g}A ({tag})" + + if session_started: + if cfg.frontload_c is not None: + st.frontloading = True + return Action("none"), st, f"front-loading at full rate ({tag})" + if i_up < cur: + return move(i_up, "landing") + return Action("none"), st, f"landing at full rate ({tag})" + + if st.frontloading: + drop_at = trip_c - cfg.soak_margin_c - k_up * sd - cfg.frontload_c + if isinstance(handle, (int, float)) and handle >= drop_at: + st.frontloading = False + if i_down < cur: + return move(i_down, f"front-load done at {handle:.1f}C") + return Action("none"), st, f"front-load done at {handle:.1f}C, full rate holds ({tag})" + return Action("none"), st, f"front-loading at full rate ({tag})" + + # The handle itself outranks every estimate: this close to the trip and + # still climbing, cut 2 A below wherever we are (or below I*, if lower). + if isinstance(handle, (int, float)) and handle >= cfg.emergency_c and rising: + return move(max(cfg.min_amps, min(cur, i_down) - 2.0), f"handle {handle:.1f}C and rising") + if i_down <= cur - 1.0: + return move(i_down, "stepping down") + if i_up >= cur + 1.0 and cur < cfg.normal_amps: + mature = st.live is not None and (age_s >= tau_s or st.live[1] <= 1.5) + rested = st.last_move_ts is None or now - st.last_move_ts >= tau_s + cool = isinstance(handle, (int, float)) and handle <= cfg.up_handle_max_c + if mature and rested and cool: + return move(i_up, "stepping up") + return Action("none"), st, f"holding {cur:g}A ({tag})" diff --git a/tests/test_settle_amp_control.py b/tests/test_settle_amp_control.py new file mode 100644 index 0000000..3803d4e --- /dev/null +++ b/tests/test_settle_amp_control.py @@ -0,0 +1,134 @@ +"""settle_amp_control.decide() is pure, like the daemon's. The behaviors +pinned here are the ones the closed-loop backtest found the hard way: one +reading per steady stretch (a draft counted every tick and held 48 A into a +64.3 C peak), stepping down on the live stretch's own reading when it runs +hotter than the session's average (the whole-history run: seven trips +without it), and the handle outranking the estimate near the trip.""" + +import importlib.util +import pathlib +import sys + +CONTRIB = pathlib.Path(__file__).parent.parent / "contrib" + + +def _load(name): + if name in sys.modules: + return sys.modules[name] + spec = importlib.util.spec_from_file_location(name, CONTRIB / f"{name}.py") + module = importlib.util.module_from_spec(spec) + sys.modules[name] = module # dataclass annotation resolution needs this + spec.loader.exec_module(module) + return module + + +_load("derate_amp_control") # settle imports Action and the trip point from it +sac = _load("settle_amp_control") + +MODEL = {"tau_min": 11.0, "rise_ref_c": 36.0, "ref_current_a": 48.0, "current_exp": 1.5, "trip_c": 65.0} +TAU_S = 11.0 * 60.0 + + +def _thermal(ts, handle, current=48.7, steady=None, prior=25.0, state="charging"): + forecast = {"basis": "trajectory" if steady is not None else "model", "window_current_a": current, + "prior_ambient_c": prior, "prior_ambient_sd_c": 0.3} + if steady is not None: + forecast.update(steady_state_c=steady, steady_state_se_c=0.1) + return {"state": state, "ts": ts, "handle_c": handle, "current_a": current, "model": MODEL, "forecast": forecast} + + +def _run(ticks, cfg=None, state=None): + """Feed (ts, handle, current, steady) ticks; return the actions and final state.""" + cfg = cfg or sac.Config() + state = state or sac.State(last_session_state="idle") + actions = [] + for ts, handle, current, steady in ticks: + action, state, _why = sac.decide(_thermal(ts, handle, current, steady), state, cfg) + actions.append(action) + return actions, state + + +def test_idle_does_nothing_and_session_end_restores(): + action, state, _ = sac.decide(_thermal(0, 25.0, state="idle"), sac.State(), sac.Config()) + assert action.kind == "none" + capped = sac.State(last_session_state="charging", cap=42.0) + action, state, _ = sac.decide(_thermal(0, 40.0, state="idle"), capped, sac.Config()) + assert (action.kind, action.value) == ("restore", 48.0) + assert state.cap is None + + +def test_front_load_holds_full_rate_while_the_handle_is_cold(): + # A plateau far over the trip, but the handle is 30 C under it: no move yet. + ticks = [(30.0 * i, 30.0 + i, 48.7, 72.0 if i >= 4 else None) for i in range(1, 12)] + actions, state = _run(ticks) + assert all(action.kind == "none" for action in actions) + assert state.frontloading + + +def test_front_load_drops_once_to_the_target_when_the_handle_nears_the_plateau(): + # Plateau 70 C at 48.7 A => theta ~33; the handle crosses the drop point. + ticks = [(30.0 * i, 30.0 + 1.2 * i, 48.7, 70.0 if i >= 4 else None) for i in range(1, 40)] + actions, state = _run(ticks) + caps = [action for action in actions if action.kind == "cap"] + assert len(caps) >= 1 + assert not state.frontloading + # The first move is a real drop, several amps, not a 1 A nibble. + assert caps[0].value <= 45.0 + + +def test_a_long_stretch_is_one_reading_not_one_per_tick(): + # 40 minutes at one current: the stretch must never be banked as finished. + ticks = [(30.0 * i, 50.0, 48.7, 60.0) for i in range(1, 81)] + _actions, state = _run(ticks) + assert state.done == [] + assert state.live is not None + + +def test_a_current_change_banks_the_stretch_and_starts_a_new_one(): + ticks = [(30.0 * i, 50.0, 48.7, 60.0) for i in range(1, 30)] + ticks += [(30.0 * i, 50.0, 44.6, 57.0) for i in range(30, 40)] + _actions, state = _run(ticks) + assert len(state.done) == 1 + assert state.stretch_current_a == 44.6 + + +def test_steps_down_on_the_live_reading_when_it_runs_hotter_than_the_average(): + # Capped at 44 A with a cool prior and a cool earlier stretch; the present + # stretch, well past half a time constant, projects over the trip. The + # fused average alone would hold; the live reading must force the step. + cfg = sac.Config(frontload_c=None) + state = sac.State(last_session_state="charging", cap=44.0, stretch_start_ts=0.0, stretch_current_a=44.6, + done=[(0.0, 26.0, 1.0)], last_move_ts=0.0) + action, state, why = sac.decide(_thermal(TAU_S, 60.0, 44.6, steady=66.0, prior=25.0), state, cfg) + assert action.kind == "cap" and action.value < 44.0, why + + +def test_emergency_cut_ignores_the_estimate(): + # The estimate says 44 A is fine; the handle is at 63.6 C and has climbed + # 0.5 C in two minutes. Cut 2 A regardless. + cfg = sac.Config(frontload_c=None) + state = sac.State(last_session_state="charging", cap=44.0, stretch_start_ts=0.0, stretch_current_a=44.6, + handle_log=[(TAU_S - 120.0, 63.1)], last_move_ts=0.0) + action, _state, why = sac.decide(_thermal(TAU_S, 63.6, 44.6, steady=62.0, prior=25.0), state, cfg) + assert (action.kind, action.value) == ("cap", 42.0), why + + +def test_one_tick_of_sensor_noise_is_not_rising(): + cfg = sac.Config(frontload_c=None) + state = sac.State(last_session_state="charging", cap=44.0, stretch_start_ts=0.0, stretch_current_a=44.6, + handle_log=[(TAU_S - 120.0, 63.5), (TAU_S - 30.0, 63.5)], last_move_ts=0.0) + action, _state, _why = sac.decide(_thermal(TAU_S, 63.6, 44.6, steady=62.0, prior=25.0), state, cfg) + assert action.kind == "none" + + +def test_step_up_waits_for_a_mature_stretch_and_a_rested_controller(): + cfg = sac.Config(frontload_c=None) + base = dict(last_session_state="charging", cap=40.0, stretch_current_a=40.6) + # Young stretch, just moved: hold even though the estimate allows more. + young = sac.State(**base, stretch_start_ts=0.0, last_move_ts=0.0) + action, _s, _w = sac.decide(_thermal(120.0, 50.0, 40.6, steady=52.0, prior=22.0), young, cfg) + assert action.kind == "none" + # A time constant on, same reading: now it may go up, and by more than 1 A. + mature = sac.State(**base, stretch_start_ts=0.0, last_move_ts=0.0) + action, _s, why = sac.decide(_thermal(TAU_S + 30.0, 50.0, 40.6, steady=52.0, prior=22.0), mature, cfg) + assert action.kind in ("cap", "restore") and action.value > 41.0, why