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.

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 `
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:

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_heat.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:

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:

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:

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.