brewpi/docs/temp_control_calibration.md

# Temperature controller calibration

## Controller overview

`components/pid/temp_controller_smith.py` implements a Smith predictor. Each tick it computes:

```
theta_ist = theta_ist_model + (theta_ist_plant − theta_ist_model_delay)
```

`theta_ist_model` comes from a fast internal `Pot` (transport delay Td = 0);
`theta_ist_model_delay` from a second `Pot` with the full configured Td.
The correction term `(theta_ist_plant − theta_ist_model_delay)` removes the dead
time from the feedback path.

The plant model (`components/plant/pot.py`) integrates:

```
temp += (delayed_power_in − L·M·(temp − theta_amb)) / (M·C) · dt
```

Four parameters must be calibrated:

| Parameter | Unit | Meaning |
|---|---|---|
| C | J/(kg·K) | Specific thermal capacity |
| M | kg | Mass of pot + water + grain |
| L | W/(kg·K) | Heat loss coefficient |
| Td | s | Transport (heater-to-sensor) delay |


## 1. Verify model-plant match: replay simulation

Feed the **recorded `power_eff`** from a log file into a fresh `Pot` instance
with the same parameters and compare the simulated temperature to `temp_ist`.

```python
from components.plant.pot import Pot

pot = Pot(dt=1.0)
pot.set_plant_params({'C': C, 'M': M, 'L': L, 'Td': Td})
pot.set_ambient_temperature(theta_amb)
pot.initial(samples[0]['temp_ist'])

sim_temp = []
for s in samples:
    pot.set_power(s['power_eff'])
    pot.process()
    sim_temp.append(pot.get_temperature())
```

Plot `sim_temp` (model) against `temp_ist` (real) and the residual
`temp_ist − sim_temp`. What divergence tells you:

| Symptom | Likely cause |
|---|---|
| Ramp slopes differ | C·M wrong |
| Phase shift between power step and temperature rise | Td wrong |
| Wrong equilibrium temperature during HOLD | L wrong |
| Residual grows over a long run | Model drift / L temperature-dependent |


## 2. Parameter identification from log data

### Td — transport delay

Cross-correlate `power_eff` with `temp_ist`. The lag at peak correlation is the
actual Td. Alternatively, find a sharp power step (start of a ramp) and measure
the visible delay before `temp_ist` begins rising.

### L — heat loss coefficient

At steady-state HOLD, `d(temp)/dt ≈ 0`, so all input power compensates losses:

```
P_hold = L · M · (T − theta_amb)
→  L = P_hold / (M · (T − theta_amb))
```

Repeat across hold phases at different temperatures and fit a line through
`(T − theta_amb)` vs `P_hold` for a more robust estimate.

### C — specific heat capacity

During a ramp where heat loss is small relative to input power:

```
C ≈ power_eff / (M · rate_ist_K_per_s)
```

`rate_ist` in the log is K/min; divide by 60 to get K/s. Use a mid-ramp window
where `power_eff` and `rate_ist` are both stable.


## 3. Controller performance metrics

These can be read directly from the log without re-simulation:

| Metric | How to compute |
|---|---|
| Overshoot | `max(temp_ist) − temp_soll` after each setpoint step |
| Settling time | First time `|temp_ist − temp_soll| < threshold` is sustained after a step |
| Steady-state error | Mean of `temp_soll − temp_ist` during HOLD phase |
| Rate tracking error | `rate_ist − rate_soll` during RAMP |

### Smith predictor correction signal

The term `theta_ist_plant − theta_ist_model_delay` is the Smith correction; a
growing correction signal indicates model drift. This is not currently written to
the log. To observe it, log `theta_ist_model` and `theta_ist_model_delay` from
`temp_controller_smith.py` alongside the existing fields.


## 4. Suggested addition to analyze_log.py

A replay panel added to `utils/analyze_log.py` would show:
- `temp_ist` (blue) vs `sim_temp` (red dashed) — model accuracy at a glance
- Residual `temp_ist − sim_temp` in a separate subplot — reveals systematic
  parameter errors vs noise floor


## 5. Heat rate noise from stirring

Stirring creates low-frequency temperature fluctuations at the sensor. The heat
rate is computed in both controllers (`temp_controller_smith.py:71-74`,
identical in `temp_controller.py`) as:

```python
heatrate = (self.theta_ist - self.last_theta_ist) / self.dt * 60
alpha = 0.1
self.heatrate_ist = (1-alpha) * self.heatrate_ist + alpha * heatrate
```

Two problems compound:

1. **Differentiation amplifies noise.** A temperature fluctuation δT at
   frequency f becomes a rate fluctuation of `2π·f·δT·60` K/min. Low-frequency
   stirrer noise (e.g. 0.05 Hz, δT = 0.05 °C) → ≈ 0.9 K/min rate noise, the
   same order as a typical `rate_soll`.

2. **The IIR filter comes after differentiation.** With α=0.1 and dt=1 s the
   time constant is τ = 9 s (cutoff ≈ 0.018 Hz). Filtering *after*
   differentiating is the least effective arrangement: noise has already been
   amplified before the filter sees it.

The noisy `heatrate_err` feeds into `pid_inner.process()` as the proportional
error. With `kp=0.08`, ±1 K/min rate ripple produces ±8% power output noise.

### Options

#### Option A — Pre-filter theta_ist before differentiating (best return for effort)

Apply a separate, heavier IIR to the raw temperature *before* the derivative:

```python
beta = 0.05  # smaller = more smoothing, more lag
self.theta_ist_filtered = (1-beta)*self.theta_ist_filtered + beta*self.theta_ist
heatrate = (self.theta_ist_filtered - self.last_theta_ist_filtered) / self.dt * 60
```

This is fundamentally better than post-filtering: the derivative of a smooth
signal stays smooth. `beta` can be tuned independently of `alpha`.

#### Option B — Sliding-window linear regression

Fit a line through the last N samples of `theta_ist` and use the slope:

```python
slope, _ = np.polyfit(t_window, theta_window, 1)
heatrate_ist = slope * 60  # K/s → K/min
```

Gives the minimum-variance rate estimate for additive white noise, with a clean
constant lag of N/2 seconds. N ≈ 20–30 s is a reasonable starting point.

#### Option C — Model-based rate (Smith predictor only)

`Pot.process()` already computes `p_pot = delayed_power_in − L·M·(temp − theta_amb)`
and integrates `d(temp)/dt = p_pot / (M·C)`. This can be used directly:

```python
rate_model = self.model.get_p_pot() / (self.model.M * self.model.C) * 60
```

(`get_p_pot()` is not yet public on `Pot` but trivially added.) The rate is
completely noise-free because it is derived from the heater power signal, not
from differentiating the sensor. Accuracy depends on model calibration — see
sections 1–2 above.

#### Option D — Reduce kp of the Heat/Cool PID

Halving `kp` from 0.08 to 0.04 halves the power noise at the cost of slower
rate convergence. Treats the symptom rather than the source.

### Recommendation

Start with **option A**: one extra IIR state, localised to the rate
computation, no other structural changes. Set `beta` so that
`τ = (1−beta)/beta · dt` is 2–3× the dominant stirrer-noise period (visible
as the ripple frequency in `rate_ist` in the log).

For the Smith predictor, **option C** is the most principled long-term
solution — it replaces a noisy numerical derivative with the model's own
analytical rate, which the Smith predictor is already tracking.