initial commit v0.2.1

custom_components/heating_mpc/physics.py

@@ -0,0 +1,169 @@
+"""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
+
+STEP_MINUTES = 5
+HEAT_LAG = 1.0 - exp(-STEP_MINUTES / 20.0)
+
+
+def bound(value: float, lo: float, hi: float) -> float:
+    return max(lo, min(hi, value))
+
+
+@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]])
+    count: int = 0
+    error_sum: float = 0.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) -> None:
+        if not all(isfinite(x) for x in (temp_before, temp_after, outside, effective_heat)):
+            return
+        y = temp_after - temp_before
+        if abs(y) > 0.8:  # sensor glitch/door opened/restarted device
+            return
+        x = [outside - temp_before, effective_heat, 1.0]
+        predicted = sum(t * v for t, v in zip(self.theta, x))
+        residual = y - predicted
+        px = [sum(self.p[i][j] * x[j] for j in range(3)) for i in range(3)]
+        forgetting = 0.999
+        denominator = forgetting + sum(x[i] * px[i] for i in range(3))
+        if denominator <= 0:
+            return
+        k = [value / denominator for value in px]
+        self.theta = [t + k[i] * residual for i, t in enumerate(self.theta)]
+        for i in range(3):
+            for j in range(3):
+                self.p[i][j] = (self.p[i][j] - k[i] * px[j]) / forgetting
+        # 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)
+
+    @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}
+
+    @classmethod
+    def from_dict(cls, source: dict[str, Any]) -> "Learner":
+        obj = cls()
+        try:
+            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
+            obj.theta, obj.p = theta, p
+            obj.count = max(0, int(source.get("count", 0)))
+            obj.error_sum = max(0., float(source.get("error_sum", 0)))
+        except (KeyError, TypeError, ValueError):
+            return cls()
+        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,
+) -> 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.
+    """
+    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(float(weights.get(room, 0)) <= 0 for room in temps):
+        raise ValueError("All rooms need a positive weight")
+    total_weight = sum(float(weights[room]) for room in temps)
+    norm = {room: float(weights[room]) / total_weight for room in temps}
+    blocks = horizon_minutes // 15
+    nsteps_per_block = 3
+    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 = simulate_step(state, outside, lag, bit, learners)
+                mismatch = sum(
+                    norm[room] * max(0.0, abs(value - target) - 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 - 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), dict(state))
+    if result is None:
+        raise RuntimeError("No feasible schedule")
+    plan, temps_60m, temps_end = result
+    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()},
+        "cost": round(best, 4),
+    }