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