buzeraandClaude Opus 5.5
099c63cd02
Confirmed: learned models were matched only by zone key, so replacing a zone's sensors, changing sensor weights, moving or replacing a heat source or turning a window kept coefficients learned for a different configuration, silently. compatibility.py fingerprints exactly what each learned component depends on: the zone measurement (sensors and relative weights; equal scaling changes nothing), each zone's heat sources, the neighbour's measurement, the zone's windows, window orientation (direction term only), the recovery time (recovery only) and other zones' openings (cross term only). Comfort weight, target offset and names are not included. Heat-source baselines have their own fingerprint (sensor and included sources). On restore, a changed measurement recalibrates the zone: base and extended loss/heat/bias keep their values as a warm start, but their uncertainty, counts and readiness start over, other terms return to the prior, and the zone's forecast statistics are dropped. Otherwise only the changed terms are reset; everything else, including other zones, is kept. Each recalibration is logged (Status model_changes, info log) and saved at once. Data from before fingerprints (0.4.0) is assumed to match its current configuration and is fingerprinted at the first start, so the next edit is detected even before a learned sample. Tests cover comfort changes, sensor replacement and weights, window orientation and size, moved and replaced heat sources, re-added zones, 0.4.0 data and a second restart, plus an HA Core test through real subentry updates. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
297 lines
13 KiB
Python
297 lines
13 KiB
Python
"""Pure-Python small grey-box thermal model and binary model predictive control.
|
|
|
|
No Home Assistant imports; this module can be tested in isolation. Step=5 minutes.
|
|
This is an experimental thermal model, not a certified boiler controller.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
from dataclasses import dataclass, field
|
|
from itertools import product
|
|
from math import exp, isfinite
|
|
from typing import Any, Callable
|
|
|
|
# (zone temperatures, radiator lag at the end of the step, step index) -> new temperatures
|
|
StepModel = Callable[[dict[str, float], float, int], dict[str, float]]
|
|
|
|
from .learning import bound_covariance, positive_definite, rls_update
|
|
|
|
STEP_MINUTES = 5
|
|
# First-order radiator lag between the relay command and the heat reaching rooms.
|
|
LAG_TAU_SECONDS = 1200.0
|
|
HEAT_LAG = 1.0 - exp(-STEP_MINUTES * 60 / LAG_TAU_SECONDS)
|
|
# Learn only when the lag proxy is known within this half-width.
|
|
LAG_LEARN_MAX_UNCERTAINTY = 0.1
|
|
# Prior covariance; also the upper bound that stops RLS windup (see _bound_covariance).
|
|
PRIOR_P = ((0.0005, 0.0, 0.0), (0.0, 0.8, 0.0), (0.0, 0.0, 0.08))
|
|
PRIOR_VARIANCE = [PRIOR_P[i][i] for i in range(3)]
|
|
|
|
|
|
def bound(value: float, lo: float, hi: float) -> float:
|
|
return max(lo, min(hi, value))
|
|
|
|
|
|
def advance_lag(lag: float, command: bool, seconds: float) -> float:
|
|
"""Move the lag proxy towards the relay command over `seconds`."""
|
|
return lag + (1.0 - exp(-seconds / LAG_TAU_SECONDS)) * (float(command) - lag)
|
|
|
|
|
|
def relay_segments(
|
|
start_on: bool, start: float, changes: list[tuple[float, str]], end: float,
|
|
) -> list[tuple[bool, float]] | None:
|
|
"""Split [start, end] into (command, seconds) segments from ordered relay changes.
|
|
|
|
`changes` holds (timestamp, new_state) pairs. Return None when the command is
|
|
unknown for part of the interval: a state other than on/off, a change to the
|
|
state already in force (so transitions were missed), or out-of-order times.
|
|
"""
|
|
segments: list[tuple[bool, float]] = []
|
|
command, since = start_on, start
|
|
for when, state in changes:
|
|
if state not in ("on", "off") or not since <= when <= end or (state == "on") == command:
|
|
return None
|
|
segments.append((command, when - since))
|
|
command, since = state == "on", when
|
|
segments.append((command, end - since))
|
|
return segments
|
|
|
|
|
|
def lag_from_time_in_state(command: bool, seconds_in_state: float) -> tuple[float, float]:
|
|
"""Estimate the lag when only the time since the last relay change is known.
|
|
|
|
Starting from any unknown lag in [0, 1], `seconds_in_state` of a constant command
|
|
leaves the lag inside an interval of width exp(-t/tau). Return its midpoint and
|
|
half-width, so the error bound is explicit rather than assuming a steady state.
|
|
"""
|
|
width = exp(-max(0.0, seconds_in_state) / LAG_TAU_SECONDS)
|
|
return (1.0 - width / 2) if command else width / 2, width / 2
|
|
|
|
|
|
def _bound_covariance(p: list[list[float]]) -> None:
|
|
"""Covariance bound for the 3-parameter base model (see learning.bound_covariance)."""
|
|
bound_covariance(p, PRIOR_VARIANCE)
|
|
|
|
|
|
@dataclass
|
|
class Learner:
|
|
"""Recursive least squares for dT=a*(Tout-T)+b*effective_heat+c."""
|
|
|
|
theta: list[float] = field(default_factory=lambda: [0.002, 0.045, 0.0])
|
|
p: list[list[float]] = field(default_factory=lambda: [list(row) for row in PRIOR_P])
|
|
count: int = 0
|
|
error_sum: float = 0.0
|
|
on_count: int = 0
|
|
off_count: int = 0
|
|
# Samples offered but not fitted (outliers, implausible inputs).
|
|
rejected_count: int = 0
|
|
|
|
def predict_delta(self, temp: float, outside: float, heat: float) -> float:
|
|
a, b, c = self.theta
|
|
return a * (outside - temp) + b * heat + c
|
|
|
|
def learn(self, temp_before: float, temp_after: float, outside: float, effective_heat: float,
|
|
duration_seconds: float = 300, command: bool | None = None) -> bool:
|
|
if (not all(isfinite(x) for x in (temp_before, temp_after, outside, effective_heat))
|
|
or not 0 <= effective_heat <= 1
|
|
or not isfinite(duration_seconds) or not 270 <= duration_seconds <= 330):
|
|
return self._reject()
|
|
y = (temp_after - temp_before) * 300 / duration_seconds
|
|
if abs(y) > 0.8: # sensor glitch/door opened/restarted device
|
|
return self._reject()
|
|
x = [outside - temp_before, effective_heat, 1.0]
|
|
result = rls_update(self.theta, self.p, x, y)
|
|
if result is None:
|
|
return self._reject()
|
|
self.theta, residual = result
|
|
_bound_covariance(self.p)
|
|
# Physical constraints; avoids unstable temperature predictions.
|
|
self.theta = [
|
|
bound(self.theta[0], 0.0, 0.04),
|
|
bound(self.theta[1], 0.0, 0.25),
|
|
bound(self.theta[2], -0.06, 0.06),
|
|
]
|
|
self.count += 1
|
|
self.error_sum += abs(residual)
|
|
if command is not None:
|
|
if command:
|
|
self.on_count += 1
|
|
else:
|
|
self.off_count += 1
|
|
return True
|
|
|
|
def _reject(self) -> bool:
|
|
self.rejected_count += 1
|
|
return False
|
|
|
|
def recalibrate(self) -> None:
|
|
"""Re-learn after the zone's measurement changed: the coefficients stay as a
|
|
starting point, but their uncertainty and every count go back to the start,
|
|
so readiness has to be earned again."""
|
|
self.p = [list(row) for row in PRIOR_P]
|
|
self.count = self.on_count = self.off_count = self.rejected_count = 0
|
|
self.error_sum = 0.0
|
|
|
|
@property
|
|
def mae(self) -> float | None:
|
|
return self.error_sum / self.count if self.count else None
|
|
|
|
def as_dict(self) -> dict[str, Any]:
|
|
return {"theta": self.theta, "p": self.p, "count": self.count, "error_sum": self.error_sum,
|
|
"on_count": self.on_count, "off_count": self.off_count,
|
|
"rejected_count": self.rejected_count}
|
|
|
|
@classmethod
|
|
def from_dict(cls, source: dict[str, Any]) -> "Learner":
|
|
obj = cls()
|
|
try:
|
|
if not isinstance(source, dict):
|
|
return obj
|
|
theta = [float(v) for v in source["theta"]]
|
|
p = [[float(v) for v in r] for r in source["p"]]
|
|
if len(theta) != 3 or len(p) != 3 or any(len(r) != 3 for r in p):
|
|
return obj
|
|
if not all(isfinite(v) for v in theta + [v for row in p for v in row]):
|
|
return obj
|
|
if not (0 <= theta[0] <= .04 and 0 <= theta[1] <= .25 and -.06 <= theta[2] <= .06):
|
|
return obj
|
|
# Tolerate float round-off relative to magnitude (older versions could
|
|
# accumulate it during windup), then apply the current covariance bound.
|
|
for i in range(3):
|
|
for j in range(i):
|
|
if abs(p[i][j] - p[j][i]) > 1e-6 * max(1.0, abs(p[i][j]), abs(p[j][i])):
|
|
return obj
|
|
_bound_covariance(p)
|
|
# A covariance must be symmetric positive definite.
|
|
if not positive_definite(p):
|
|
return obj
|
|
counts = [source.get(key, 0)
|
|
for key in ("count", "on_count", "off_count", "rejected_count")]
|
|
if any(type(v) is not int or v < 0 for v in counts):
|
|
return obj
|
|
count, on_count, off_count, rejected_count = counts
|
|
error_sum = float(source.get("error_sum", 0))
|
|
if not isfinite(error_sum) or error_sum < 0 or on_count + off_count > count:
|
|
return obj
|
|
return cls(theta, p, count, error_sum, on_count, off_count, rejected_count)
|
|
except (KeyError, TypeError, ValueError, OverflowError):
|
|
return obj
|
|
|
|
|
|
def simulate_step(temps: dict[str, float], outside: float, lag: float,
|
|
command: bool, learners: dict[str, Learner]) -> tuple[dict[str, float], float]:
|
|
next_lag = lag + HEAT_LAG * (float(command) - lag)
|
|
new = {name: value + learners[name].predict_delta(value, outside, next_lag)
|
|
for name, value in temps.items()}
|
|
return new, next_lag
|
|
|
|
|
|
def optimize(
|
|
temps: dict[str, float], outside: float, effective_heat: float,
|
|
target: float, weights: dict[str, float], learners: dict[str, Learner],
|
|
is_on: bool, seconds_since_change: float, horizon_minutes: int = 120,
|
|
minimum_on_minutes: int = 15, minimum_off_minutes: int = 10,
|
|
switching_cost: float = .15, heating_cost: float = .025,
|
|
comfort_band: float = .15, target_offsets: dict[str, float] | None = None,
|
|
step_model: StepModel | None = None,
|
|
) -> dict[str, Any]:
|
|
"""Enumerate 15min ON/OFF plans; honor minimum ON/OFF durations.
|
|
|
|
120min = 256 plans (fast on Raspberry Pi 5). All rooms share one burner demand.
|
|
`target_offsets` shifts individual zones' targets (°C) relative to `target`.
|
|
`step_model` advances all zones together (electrical heat, ventilation,
|
|
neighbours); by default each zone uses its learner's base model.
|
|
"""
|
|
if not temps or set(temps) != set(learners):
|
|
raise ValueError("Missing room temperatures or learner")
|
|
if not 30 <= horizon_minutes <= 180 or horizon_minutes % 15:
|
|
raise ValueError("Horizon must be 30..180 minutes in 15-minute steps")
|
|
if not isfinite(target) or not isfinite(outside):
|
|
raise ValueError("Invalid temperature")
|
|
if any(not isfinite(t) for t in temps.values()):
|
|
raise ValueError("Invalid room temperature")
|
|
if any(not isfinite(float(weights.get(room, 0))) or float(weights.get(room, 0)) <= 0 for room in temps):
|
|
raise ValueError("All rooms need a positive weight")
|
|
if not isfinite(effective_heat) or not 0 <= effective_heat <= 1:
|
|
raise ValueError("Invalid heating lag")
|
|
if any(not isfinite(v) or v < 0 for v in (
|
|
seconds_since_change, minimum_on_minutes, minimum_off_minutes,
|
|
switching_cost, heating_cost, comfort_band,
|
|
)):
|
|
raise ValueError("Invalid optimizer parameter")
|
|
total_weight = sum(float(weights[room]) for room in temps)
|
|
if not isfinite(total_weight):
|
|
raise ValueError("Invalid total weight")
|
|
norm = {room: float(weights[room]) / total_weight for room in temps}
|
|
offsets = target_offsets or {}
|
|
if any(not isfinite(float(offsets.get(room, 0.0))) for room in temps):
|
|
raise ValueError("Invalid target offset")
|
|
targets = {room: target + float(offsets.get(room, 0.0)) for room in temps}
|
|
mean_target = sum(norm[room] * targets[room] for room in temps)
|
|
blocks = horizon_minutes // 15
|
|
nsteps_per_block = 3
|
|
|
|
def advance(state: dict[str, float], lag: float, bit: bool, k: int) -> tuple[dict[str, float], float]:
|
|
if step_model is None:
|
|
return simulate_step(state, outside, lag, bit, learners)
|
|
next_lag = lag + HEAT_LAG * (float(bit) - lag)
|
|
return step_model(state, next_lag, k), next_lag
|
|
best = float("inf")
|
|
result = None
|
|
min_on_s = minimum_on_minutes * 60
|
|
min_off_s = minimum_off_minutes * 60
|
|
elapsed_initial = max(0.0, seconds_since_change)
|
|
# More sequences than 4096 are deliberately excluded by horizon validation.
|
|
for bits in product((False, True), repeat=blocks):
|
|
current = bool(is_on)
|
|
elapsed_s = elapsed_initial
|
|
state = dict(temps)
|
|
lag = effective_heat
|
|
cost = 0.0
|
|
at_60 = None
|
|
feasible = True
|
|
for idx, bit in enumerate(bits):
|
|
if bit != current:
|
|
minimum = min_on_s if current else min_off_s
|
|
if elapsed_s < minimum:
|
|
feasible = False
|
|
break
|
|
cost += switching_cost
|
|
elapsed_s = 0.0
|
|
current = bit
|
|
for tick in range(nsteps_per_block):
|
|
state, lag = advance(state, lag, bit, idx * nsteps_per_block + tick)
|
|
if not all(isfinite(v) for v in state.values()):
|
|
raise ValueError("Non-finite simulated temperature")
|
|
mismatch = sum(
|
|
norm[room] * max(0.0, abs(value - targets[room]) - comfort_band) ** 2
|
|
for room, value in state.items()
|
|
)
|
|
# A lower bound on weighted average comfort also matters.
|
|
avg = sum(norm[r] * value for r, value in state.items())
|
|
mismatch += max(0.0, abs(avg - mean_target) - comfort_band) ** 2
|
|
cost += 8.0 * mismatch + (heating_cost if bit else 0.0)
|
|
if idx * 15 + (tick + 1) * STEP_MINUTES == 60:
|
|
at_60 = dict(state)
|
|
elapsed_s += 15 * 60
|
|
if feasible and cost < best:
|
|
best = cost
|
|
result = (bits, at_60 or {}, dict(state))
|
|
if result is None:
|
|
raise RuntimeError("No feasible schedule")
|
|
plan, temps_60m, temps_end = result
|
|
# Re-simulate the chosen plan for its trajectory, keeping the search loop lean.
|
|
trajectory: dict[str, list[float]] = {room: [] for room in temps}
|
|
state, lag = dict(temps), effective_heat
|
|
for idx, bit in enumerate(plan):
|
|
for tick in range(nsteps_per_block):
|
|
state, lag = advance(state, lag, bit, idx * nsteps_per_block + tick)
|
|
for room, value in state.items():
|
|
trajectory[room].append(round(value, 3))
|
|
return {
|
|
"proposal": bool(plan[0]),
|
|
"plan_15min": [int(v) for v in plan],
|
|
"predicted_60m": {r: round(t, 3) for r, t in temps_60m.items()},
|
|
"predicted_end": {r: round(t, 3) for r, t in temps_end.items()},
|
|
"trajectory_15min": trajectory,
|
|
"cost": round(best, 4),
|
|
}
|