Bound learner covariance to stop RLS windup

With forgetting 0.999 and no heat excitation (summer, long OFF periods),
the heat-coefficient variance grew by 1/0.999 per sample without limit:
about 8.7e12 after 104 simulated OFF days. The first ON samples could
then swing coefficients to their bounds, and accumulated asymmetry
approached the absolute 1e-8 restore tolerance, past which a restart
silently reset the whole model.

Symmetrize P and bound its diagonal by the prior after every update
(row/column scaling keeps it positive definite). Restore uses a relative
round-off tolerance and applies the same bound, so wound-up v0.2.1
models load instead of being discarded. Payload schema is unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

custom_components/heating_mpc/physics.py

@@ -12,18 +12,40 @@ from typing import Any
 
 STEP_MINUTES = 5
 HEAT_LAG = 1.0 - exp(-STEP_MINUTES / 20.0)
+# 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))
 
 
 def bound(value: float, lo: float, hi: float) -> float:
     return max(lo, min(hi, value))
 
 
+def _bound_covariance(p: list[list[float]]) -> None:
+    """Keep P symmetric and no less certain than the prior, in place.
+
+    With forgetting, unexcited directions (e.g. heat input during a summer with the
+    relay off) grow by 1/forgetting each step without limit. Scaling row and column i
+    together is a congruence transform, so positive definiteness is preserved.
+    """
+    for i in range(3):
+        for j in range(i):
+            p[i][j] = p[j][i] = (p[i][j] + p[j][i]) / 2
+    for i in range(3):
+        if p[i][i] > PRIOR_P[i][i]:
+            scale = (PRIOR_P[i][i] / p[i][i]) ** .5
+            for j in range(3):
+                p[i][j] *= scale
+            for j in range(3):
+                p[j][i] *= scale
+            p[i][i] = PRIOR_P[i][i]
+
+
 @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: [[0.0005, 0., 0.], [0., 0.8, 0.], [0., 0., 0.08]])
+    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
@@ -57,6 +79,7 @@ class Learner:
         for i in range(3):
             for j in range(3):
                 self.p[i][j] = (self.p[i][j] - k[i] * px[j]) / forgetting
+        _bound_covariance(self.p)
         # Physical constraints; avoids unstable temperature predictions.
         self.theta = [
             bound(self.theta[0], 0.0, 0.04),
@@ -94,12 +117,17 @@ class Learner:
                 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 (Cholesky check).
             lower = [[0.0] * 3 for _ in range(3)]
             for i in range(3):
                 for j in range(i + 1):
-                    if abs(p[i][j] - p[j][i]) > 1e-8:
-                        return obj
                     value = p[i][j] - sum(lower[i][k] * lower[j][k] for k in range(j))
                     if i == j:
                         if value <= 0 or not isfinite(value):