diff --git a/ashvale/api.py b/ashvale/api.py index ef87318..678498d 100644 --- a/ashvale/api.py +++ b/ashvale/api.py @@ -136,6 +136,9 @@ class SettingsIn(BaseModel): altitude_m: Optional[float] = Field(None, ge=-430, le=9000) latitude: Optional[float] = Field(None, ge=-90, le=90) longitude: Optional[float] = Field(None, ge=-180, le=180) + heating: Optional[bool] = None + heating_setpoint_c: Optional[float] = Field(None, ge=5, le=35) + thermal_time_constant_h: Optional[float] = Field(None, ge=0.1, le=24) hum_psychrometric: Optional[bool] = None led_enabled: Optional[bool] = None led_fps: Optional[float] = Field(None, ge=4, le=30) @@ -647,6 +650,9 @@ def get_settings() -> Dict: "latitude": CONFIG.site.latitude, "longitude": CONFIG.site.longitude, "timezone": CONFIG.site.timezone, + "heating": CONFIG.site.heating, + "heating_setpoint_c": CONFIG.site.heating_setpoint_c, + "thermal_time_constant_h": CONFIG.site.thermal_time_constant_h, "name": CONFIG.site.name}, "sensor": {"hum_psychrometric": CONFIG.sensor.hum_psychrometric, "cpu_heat_k": round(st.tracker.compensator.k, 4), @@ -693,6 +699,24 @@ def post_settings(body: SettingsIn) -> Dict: # history is now inconsistent with the new value until re-derived. needs_recompute = needs_recompute or name == "altitude_m" + # Turning the thermostat model on or off changes which process the heads are + # fitting, so it is a regime change and gets the same treatment as a door. + if body.heating is not None and body.heating != CONFIG.site.heating: + CONFIG.site.heating = bool(body.heating) + patch.setdefault("site", {})["heating"] = bool(body.heating) + applied.append(f"heating {'on' if body.heating else 'off'}") + st.store.log_event("discontinuity", "warn", + f"heating {'on' if body.heating else 'off'}") + st.monitor.retrain_requested = True + + for name, value, label in ( + ("heating_setpoint_c", body.heating_setpoint_c, "setpoint"), + ("thermal_time_constant_h", body.thermal_time_constant_h, "time constant")): + if value is not None and value != getattr(CONFIG.site, name): + applied.append(f"{label} {getattr(CONFIG.site, name)} -> {value}") + setattr(CONFIG.site, name, float(value)) + patch.setdefault("site", {})[name] = float(value) + if body.hum_psychrometric is not None and \ body.hum_psychrometric != CONFIG.sensor.hum_psychrometric: CONFIG.sensor.hum_psychrometric = bool(body.hum_psychrometric) diff --git a/ashvale/config.py b/ashvale/config.py index 6d4c768..36f935c 100644 --- a/ashvale/config.py +++ b/ashvale/config.py @@ -57,7 +57,25 @@ class SiteConfig: # keep predicting the old regime for two days unless told. POST # /api/environment marks the moment and asks for a retrain. environment: str = "indoor" # indoor | sheltered | outdoor - enclosure: str = "closed" # closed | ventilated | open # honest flag, changes how forecasts are worded + enclosure: str = "closed" # closed | ventilated | open + + # Central heating or air conditioning holding the room at a setpoint. + # + # This is a genuine change of process, not a label. A free-running room + # follows outdoor forcing and drifts; a thermostatted one is a closed loop + # that pulls back toward heating_setpoint_c whenever it strays. Persistence + # ("tomorrow equals today") is the wrong baseline for a controlled system, + # because the truth is "it returns to the setpoint". + # + # thermal_time_constant_h is how fast that pull acts: the time to close + # about 63% of a gap. A small well-insulated flat with responsive heating is + # under an hour; a large draughty house with slow radiators is several. If + # you do not know it, leave it: the ensemble weights this member against the + # others from measured error, so a wrong constant costs accuracy, not + # correctness. + heating: bool = False + heating_setpoint_c: float = 21.0 + thermal_time_constant_h: float = 1.5 # honest flag, changes how forecasts are worded @dataclass diff --git a/ashvale/dashboard.py b/ashvale/dashboard.py index 6e4dd57..a87bd82 100644 --- a/ashvale/dashboard.py +++ b/ashvale/dashboard.py @@ -489,6 +489,22 @@ DASHBOARD_HTML = r"""

Closing a door changes how strongly the sensor couples to outside. The heads carry about 55 hours of memory, so tell them rather than waiting two days.

+
+
+
+ Heated or cooled to a setpoint +

A thermostat makes the room a closed loop: it returns to the setpoint instead of drifting. Persistence is the wrong baseline for that, so this adds a fourth ensemble member and lets the Hedge weights decide if it earns its place. Humidity follows at constant vapour pressure, which is why a heated house is dry.

+
+ +
+
+ + +
+
+
@@ -1111,6 +1127,10 @@ async function loadSettings() { el('s-lon').value = site.longitude; el('s-psy').innerText = sen.hum_psychrometric ? 'on' : 'off'; el('s-psy').className = sen.hum_psychrometric ? PILL_ON : PILL_OFF; + el('s-heat').innerText = site.heating ? 'on' : 'off'; + el('s-heat').className = site.heating ? PILL_ON : PILL_OFF; + el('s-setpoint').value = site.heating_setpoint_c; + el('s-tau').value = site.thermal_time_constant_h; } loaders.settings = loadSettings; @@ -1147,6 +1167,13 @@ el('s-led').addEventListener('click', () => postSettings({led_enabled: !(settingsDoc.server||{}).led_enabled}, 's-msg')); el('s-psy').addEventListener('click', () => postSettings({hum_psychrometric: !(settingsDoc.sensor||{}).hum_psychrometric}, 's-msg')); +el('s-heat').addEventListener('click', () => + postSettings({heating: !(settingsDoc.site||{}).heating}, 's-msg')); +for (const id of ['s-setpoint','s-tau']) { + el(id).addEventListener('change', () => postSettings({ + heating_setpoint_c: Number(el('s-setpoint').value), + thermal_time_constant_h: Number(el('s-tau').value)}, 's-msg')); +} el('s-fps').addEventListener('input', e => el('s-fps-val').innerText = e.target.value + ' fps'); el('s-fps').addEventListener('change', e => postSettings({led_fps: Number(e.target.value)}, 's-msg')); diff --git a/ashvale/methods.py b/ashvale/methods.py index 3df7352..5abe74a 100644 --- a/ashvale/methods.py +++ b/ashvale/methods.py @@ -178,6 +178,10 @@ def pipeline(cfg) -> List[Dict[str, Any]]: "the model then detonates at sunrise. The trace is capped. " "This is the most common way a field RLS deployment dies.", "math": [ + r"dT_{set}(h) = (T_{set} - T_{now})\left(1 - e^{-h/\tau}\right)" + r"\qquad\text{(thermostat member, first-order closed loop)}", + r"RH(h) = RH_{now}\,\frac{e_s(T_{now})}{e_s(T_{now} + dT_{set}(h))}" + r"\qquad\text{(heating adds no moisture, so dew point is conserved)}", r"\hat{\theta} = \arg\min_{\theta}\; \sum_{i=1}^{t}" r"\lambda^{\,t-i}\big(y_i - \theta^{\top}x_i\big)^{2}" r"\qquad\text{(exponentially weighted least squares)}", diff --git a/ashvale/models/nowcast.py b/ashvale/models/nowcast.py index ad5ea19..9593f33 100644 --- a/ashvale/models/nowcast.py +++ b/ashvale/models/nowcast.py @@ -43,7 +43,7 @@ import numpy as np from ..features import N_FEATURES, Standardiser, supervised_pairs from .rls import AdaptiveConformal, RecursiveLeastSquares -MEMBERS = ("persistence", "climatology", "learned") +MEMBERS = ("persistence", "climatology", "learned", "setpoint") class ForecastHead: @@ -65,9 +65,11 @@ class ForecastHead: # -------------------------------------------------------- prediction def predict(self, x: np.ndarray, anchor: float, - climatology_delta: float = 0.0) -> Dict[str, float]: + climatology_delta: float = 0.0, + setpoint_delta: float = 0.0) -> Dict[str, float]: learned_delta = self.model.predict(x) - deltas = np.array([0.0, float(climatology_delta), float(learned_delta)]) + deltas = np.array([0.0, float(climatology_delta), float(learned_delta), + float(setpoint_delta)]) blended = float(np.dot(self.weights, deltas)) mu = float(anchor + blended) sigma = self.model.predict_std(x, self.model.noise_var) @@ -85,10 +87,12 @@ class ForecastHead: # ---------------------------------------------------------- learning def learn(self, x: np.ndarray, anchor: float, truth: float, - climatology_delta: float = 0.0) -> float: + climatology_delta: float = 0.0, + setpoint_delta: float = 0.0) -> float: """One supervised step given a matured target.""" deltas = np.array([0.0, float(climatology_delta), - float(self.model.predict(x))]) + float(self.model.predict(x)), + float(setpoint_delta)]) member_pred = anchor + deltas losses = np.abs(member_pred - truth) @@ -119,9 +123,23 @@ class ForecastHead: h = cls(s["target"], s["horizon_s"]) h.model = RecursiveLeastSquares.from_dict(s["model"]) h.conformal = AdaptiveConformal.from_dict(s["conformal"]) - h.weights = np.array(s["weights"], dtype=float) + w = np.array(s["weights"], dtype=float) + if w.size != len(MEMBERS): + # A saved head from before the setpoint member existed. Reinitialise + # uniformly rather than guessing: the Hedge weights re-converge in + # about a day, which is far cheaper than silently mismatching a + # member to the wrong loss and corrupting every blend until someone + # notices. + w = np.ones(len(MEMBERS)) / len(MEMBERS) + h.weights = w h.eta = s["eta"] - h.member_mae = np.array(s["member_mae"], dtype=float) + mae = np.array(s["member_mae"], dtype=float) + # Same migration as the weights. Missing this one did not fail on load, + # it failed later inside learn() on a shape mismatch, which is a worse + # place to find out. + if mae.size != len(MEMBERS): + mae = np.zeros(len(MEMBERS)) + h.member_mae = mae h.n_scored = s.get("n_scored", 0) return h @@ -149,7 +167,7 @@ class NowcastEnsemble: def fit(self, X: np.ndarray, valid: np.ndarray, series: Dict[str, np.ndarray], climatology=None, grid_ts: Optional[np.ndarray] = None, - passes: int = 1, max_pairs: int = 2500) -> Dict[str, int]: + passes: int = 1, max_pairs: int = 2500, setpoint_fn=None) -> Dict[str, int]: """Batch-update every head from history. `max_pairs` bounds the work per head to the most recent samples. @@ -187,7 +205,8 @@ class NowcastEnsemble: clim[-mask_len:] = clim_fut - clim_now for _ in range(max(int(passes), 1)): for i in range(Xa.shape[0]): - head.learn(Xa[i], anchor[i], anchor[i] + dy[i], clim[i]) + head.learn(Xa[i], anchor[i], anchor[i] + dy[i], clim[i], + setpoint_fn(target, h, anchor[i]) if setpoint_fn else 0.0) counts[f"{target}@{h}"] = int(Xa.shape[0]) self.trained_rows = int(X.shape[0]) return counts @@ -195,7 +214,7 @@ class NowcastEnsemble: # --------------------------------------------------------- inference def forecast(self, x_raw: np.ndarray, anchors: Dict[str, float], now: float, - climatology=None) -> Dict[str, Dict[int, Dict[str, float]]]: + climatology=None, setpoint_fn=None) -> Dict[str, Dict[int, Dict[str, float]]]: x = self.scaler.transform(np.atleast_2d(x_raw))[0] out: Dict[str, Dict[int, Dict[str, float]]] = {} for target in self.targets: @@ -206,7 +225,8 @@ class NowcastEnsemble: if climatology is not None and climatology.ready: clim_delta = float(climatology.predict(target, np.array([now + h]))[0] - climatology.predict(target, np.array([now]))[0]) - out[target][h] = self.heads[(target, h)].predict(x, anchor, clim_delta) + sp = setpoint_fn(target, h, anchor) if setpoint_fn else 0.0 + out[target][h] = self.heads[(target, h)].predict(x, anchor, clim_delta, sp) return out def diagnostics(self) -> List[Dict]: diff --git a/ashvale/station.py b/ashvale/station.py index 60d06dd..923fc1f 100644 --- a/ashvale/station.py +++ b/ashvale/station.py @@ -34,6 +34,7 @@ from __future__ import annotations import asyncio import json +import math import time from pathlib import Path from typing import Any, Dict, List, Optional @@ -382,6 +383,55 @@ class Station: return {"rows": written, "seconds": round(secs, 2), "k": comp.k, "hum_offset": hcomp.offset} + def _setpoint_delta(self, target: str, horizon_s: int, anchor: float) -> float: + """Where a thermostatted room is heading, as a delta from now. + + A controlled room is first order: the heating closes the gap to the + setpoint exponentially, so after time h the remaining error is + exp(-h/tau) of what it was. The expected change is therefore + + dT(h) = (T_set - T_now) * (1 - exp(-h / tau)) + + which is zero at h=0 and asymptotes to the full correction. That is a + much better statement about a heated room than persistence, which claims + the room stays wherever it happens to be. + + Humidity follows for free and is the part people get wrong. Heating adds + no moisture, so vapour pressure is what is conserved, not relative + humidity. Warm the air and RH falls even though nothing was dried: + + RH(h) = RH_now * es(T_now) / es(T_now + dT(h)) + + This is why a heated house in winter is dry. Pressure is unaffected: a + thermostat cannot move the synoptic field, so that member stays at zero + and the ensemble will correctly ignore it. + + Returns 0.0 when heating is off, which makes this member identical to + persistence and therefore harmless. + """ + site = self.cfg.site + if not site.heating: + return 0.0 + tau_s = max(float(site.thermal_time_constant_h), 0.05) * 3600.0 + closed = 1.0 - math.exp(-float(horizon_s) / tau_s) + + temp_now = self.live.get("temp_smooth") + if temp_now is None: + return 0.0 + d_temp = (float(site.heating_setpoint_c) - float(temp_now)) * closed + + if target == "temperature": + return d_temp + if target == "humidity": + # Constant vapour pressure, so RH moves only because es(T) moved. + es_now = float(physics.saturation_vapour_pressure(temp_now)) + es_fut = float(physics.saturation_vapour_pressure(temp_now + d_temp)) + if es_fut <= 1e-9: + return 0.0 + rh_now = float(anchor) + return float(np.clip(rh_now * es_now / es_fut, 0.0, 100.0)) - rh_now + return 0.0 + def set_environment(self, environment: Optional[str] = None, enclosure: Optional[str] = None, note: str = "") -> Dict: @@ -467,7 +517,8 @@ class Station: grid_ts, cols, X, valid = built clim_scores = self.climatology.fit(grid_ts, cols, valid) - counts = self.nowcast.fit(X, valid, cols, self.climatology, grid_ts) + counts = self.nowcast.fit(X, valid, cols, self.climatology, grid_ts, + setpoint_fn=self._setpoint_delta) self.last_train = time.time() self.monitor.clear_retrain_flag() @@ -503,7 +554,7 @@ class Station: "humidity": float(self.live.get("hum_smooth", cols["humidity"][-1])), "pressure": float(self.live.get("press_slp", cols["pressure"][-1])), } - fc = self.nowcast.forecast(x_now, anchors, now, self.climatology) + fc = self.nowcast.forecast(x_now, anchors, now, self.climatology, setpoint_fn=self._setpoint_delta) bundle: Dict[str, Any] = {"issued_ts": now, "anchors": anchors, "targets": {}} for target, per_h in fc.items(): diff --git a/docs/DESIGN.md b/docs/DESIGN.md index f707243..6acf332 100644 --- a/docs/DESIGN.md +++ b/docs/DESIGN.md @@ -240,6 +240,45 @@ re-weights within about a day when the season turns. 15-minute pressure it typically parks most of its weight on persistence. That is correct behaviour surfaced honestly, not a defect to engineer away. +### The thermostat member + +A room held at a setpoint is not the same process as a room that is free to +drift. It is a closed loop, and persistence, the baseline everything here is +scored against, is simply the wrong statement about it: the truth is not "it +stays where it is", it is "it returns to the setpoint". + +So when `site.heating` is on, the ensemble gains a fourth member: + +``` +dT_set(h) = (T_set - T_now) * (1 - exp(-h / tau)) +``` + +First order, because that is what a controlled system is: `tau` is the time to +close about 63% of the gap. Zero at h = 0, asymptotic to the full correction. + +Humidity follows and is the part that is easy to get wrong. Heating adds no +moisture, so what is conserved is vapour pressure, not relative humidity: + +``` +RH(h) = RH_now * es(T_now) / es(T_now + dT_set(h)) +``` + +Warm the air and RH falls although nothing was dried. This is why a heated house +in winter is dry, and the test asserts the dew point is unchanged to 1e-6. + +Pressure gets zero: a thermostat cannot move the synoptic field. + +**It is offered, not imposed.** The Hedge weights score this member against the +others on realised error like any other, so a wrong `tau` or a setpoint you +forgot to update costs accuracy and gets down-weighted, rather than quietly +biasing every forecast. With heating off the member returns zero, which makes it +identical to persistence and therefore harmless. + +Adding it changed the member count from three to four, so `ForecastHead.from_dict` +reinitialises `weights` **and** `member_mae` when a saved head has the old +length. Missing the second one did not fail on load: it failed later inside +`learn()` on a broadcast error, which is a much worse place to find out. + ### Adaptive conformal intervals Split conformal is valid only under exchangeability, and weather is emphatically diff --git a/tests/test_estimation.py b/tests/test_estimation.py index 787a480..18bd68b 100644 --- a/tests/test_estimation.py +++ b/tests/test_estimation.py @@ -177,3 +177,74 @@ def test_kalman_state_round_trips_through_dict(): back = KalmanCV.from_dict(kf.to_dict()) assert back.level == pytest.approx(kf.level) assert back.rate == pytest.approx(kf.rate) + + +# ---------------------------------------------------------------- thermostat + +def test_thermostat_reversion_is_first_order_and_preserves_dew_point(): + """A heated room is a closed loop, and heating adds no moisture. + + Two properties, both easy to get wrong. The temperature must close the gap + to the setpoint exponentially rather than jumping or drifting, and the + implied humidity change must leave the dew point exactly where it was: RH + falls only because es(T) rose, which is why a heated house in winter is dry. + """ + import math + + from ashvale.config import load_config + from ashvale.physics import dew_point, saturation_vapour_pressure + from ashvale.station import Station + + cfg = load_config() + cfg.site.heating = True + cfg.site.heating_setpoint_c = 23.0 + cfg.site.thermal_time_constant_h = 1.5 + st = Station(cfg) + st.live = {"temp_smooth": 18.0} + + tau = 1.5 * 3600.0 + for h in (900, 3600, 10800, 86400): + expected = (23.0 - 18.0) * (1.0 - math.exp(-h / tau)) + assert st._setpoint_delta("temperature", h, 18.0) == pytest.approx(expected, rel=1e-9) + + # monotonic toward the setpoint, never past it + deltas = [st._setpoint_delta("temperature", h, 18.0) + for h in (900, 3600, 10800, 21600, 86400)] + assert all(a < b for a, b in zip(deltas, deltas[1:])) + assert deltas[-1] <= 5.0 + 1e-9 + + # dew point invariant + t0, rh0 = 18.0, 55.0 + d_t = st._setpoint_delta("temperature", 86400, t0) + d_rh = st._setpoint_delta("humidity", 86400, rh0) + assert float(dew_point(t0 + d_t, rh0 + d_rh)) == pytest.approx( + float(dew_point(t0, rh0)), abs=1e-6) + assert d_rh < 0.0, "warming a room at constant moisture must lower RH" + assert float(saturation_vapour_pressure(t0 + d_t)) > float( + saturation_vapour_pressure(t0)) + + # a thermostat cannot move the synoptic field + assert st._setpoint_delta("pressure", 86400, 1013.0) == 0.0 + + # and off, the member is exactly persistence + cfg.site.heating = False + assert st._setpoint_delta("temperature", 86400, 18.0) == 0.0 + assert st._setpoint_delta("humidity", 86400, 55.0) == 0.0 + + +def test_forecast_head_migrates_state_from_before_the_setpoint_member(): + """An old save has three weights where there are now four.""" + from ashvale.models.nowcast import MEMBERS, ForecastHead + + h = ForecastHead(target="temperature", horizon_s=900, n_features=4) + state = h.to_dict() + state["weights"] = [0.2, 0.3, 0.5] # a pre-setpoint save + state["member_mae"] = [0.4, 0.5, 0.6] + back = ForecastHead.from_dict(state) + assert back.weights.size == len(MEMBERS) + assert float(back.weights.sum()) == pytest.approx(1.0) + # member_mae must migrate too. Missing it did not fail on load, it failed + # later inside learn() on a broadcast error, which is a worse place to + # discover a migration bug. + assert back.member_mae.size == len(MEMBERS) + back.learn(np.zeros(4), 20.0, 20.5, 0.1, 0.2) # must not raise