Inspect a bounded thermal PID simulation
Simulate derivative-on-measurement control, actuator saturation and conditional integration, with an explicitly assumed plant.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the research question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
from IPython.display import display as _display
def display(value):
if isinstance(value,pd.DataFrame):
table=value.to_html(max_rows=12,escape=True)
style="<style>.vibeit-readable-table{max-width:100%;overflow-x:auto}.vibeit-readable-table table{width:max-content;max-width:none;border-collapse:collapse}.vibeit-readable-table td,.vibeit-readable-table th{white-space:nowrap!important;word-break:normal!important;overflow-wrap:normal!important}</style>"
return _display(HTML(style+'<div class="vibeit-readable-table">'+table+'</div>'))
return _display(value)
RECIPE_ID='en06-pid-control'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
Load the embedded snapshot
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'models': '48fd2e6e706650fe2ee0c0ad80b6d794a546a1dec78dd9d370ac09db3cc3adbf'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['models']
Read the method functions
def pid_thermal_simulation(Kp=4,Ki=.2,Kd=6,tau_s=60,gain_degC_per_percent=1,setpoint_degC=60,ambient_degC=20,dt_s=1,duration_s=600,antiwindup=True,maximum_percent=100,disturbance_start_s=300,disturbance_equilibrium_offset_degC=-10):
"""Declared thermal plant and conditional-integration PID; derivative on measurement."""
if min(tau_s,dt_s,duration_s,maximum_percent)<=0 or min(Kp,Ki,Kd)<0 or dt_s>=tau_s:
raise ValueError('Invalid controller parameters or overly large plant time step')
count=int(np.floor(duration_s/dt_s));time=np.arange(count+1)*dt_s
temperature=np.full(count+1,ambient_degC,dtype=float);control=np.zeros(count+1);integral_trace=np.zeros(count+1);unclamped=np.zeros(count+1)
integral=0;previous=ambient_degC
for j in range(count):
error=setpoint_degC-temperature[j];derivative=(temperature[j]-previous)/dt_s
request=Kp*error+integral-Kd*derivative;applied=float(np.clip(request,0,maximum_percent))
blocked=(request>=maximum_percent and error>0) or (request<=0 and error<0)
if not antiwindup or not blocked:integral+=Ki*error*dt_s
control[j]=applied;integral_trace[j]=integral;unclamped[j]=request
disturbance=disturbance_equilibrium_offset_degC if time[j]>=disturbance_start_s else 0
previous=temperature[j]
temperature[j+1]=temperature[j]+dt_s*(ambient_degC-temperature[j]+gain_degC_per_percent*applied+disturbance)/tau_s
control[-1]=control[-2];integral_trace[-1]=integral_trace[-2];unclamped[-1]=unclamped[-2]
return dict(t=time,temperature_degC=temperature,control_percent=control,integral_percent=integral_trace,unclamped_percent=unclamped,error_degC=setpoint_degC-temperature)
print("Method ready: pid_thermal_simulation")
Method ready: pid_thermal_simulation
Data and model boundary
The building record contains real ten-minute energy and environmental observations from one source building. Airport weather fields include upstream interpolation and are not colocated calibration references. Beam, circuit and PID examples are declared numerical design models; computed responses do not establish hardware performance or commissioning acceptance.
Steps
1. Declare the plant and controller units
The assumed thermal plant relaxes toward ambient+gain·heater duty plus a declared disturbance. Heater duty is limited to 0–100%. Kp is %/°C, Ki %/(°C·s), Kd %·s/°C. Derivative is taken on the measured temperature, not on a step change of reference.
p=snapshot("models")["pid"]
arguments=dict(Kp=4,Ki=.2,Kd=6,tau_s=p["time_constant_s"],gain_degC_per_percent=p["plant_gain_degC_per_percent"],setpoint_degC=p["setpoint_degC"],ambient_degC=p["ambient_degC"],dt_s=p["step_s"],duration_s=p["duration_s"])
controlled=pid_thermal_simulation(**arguments,antiwindup=True)
unprotected=pid_thermal_simulation(**arguments,antiwindup=False)
fig,ax=plt.subplots(figsize=(9,4));ax.plot(controlled["t"],controlled["temperature_degC"],label="Conditional integration")
ax.plot(unprotected["t"],unprotected["temperature_degC"],label="No anti-windup")
ax.axhline(p["setpoint_degC"],color="black",linestyle="--",label="Declared setpoint");ax.axvline(300,color="gray",linestyle=":")
ax.set(xlabel="Model time (s)",ylabel="Model temperature (°C)",title="Bounded thermal control; equilibrium disturbance at 300 s");ax.legend();plt.tight_layout();plt.show()
2. Compare requested and applied control
The unclamped controller request can exceed actuator limits, while applied duty is bounded. This is an intentionally simulated actuator, not an instruction to control real hardware. Saturation is part of the declared model rather than a value to hide.
fig,ax=plt.subplots(figsize=(9,4));ax.plot(controlled["t"],controlled["unclamped_percent"],label="Requested duty")
ax.plot(controlled["t"],controlled["control_percent"],label="Applied duty")
ax.axhline(100,color="gray",linestyle="--");ax.axhline(0,color="gray",linestyle="--")
ax.set(xlabel="Model time (s)",ylabel="Heater duty (%)",title="Explicit actuator saturation");ax.legend();plt.tight_layout();plt.show()
3. Inspect the integral state
Conditional integration holds the integrator when an already saturated request and error would push farther into saturation. The integral contribution is a controller state in duty units; it is not temperature or physical stored heat.
fig,ax=plt.subplots(figsize=(9,4));ax.plot(controlled["t"],controlled["integral_percent"],label="Conditional integration")
ax.plot(unprotected["t"],unprotected["integral_percent"],label="No anti-windup")
ax.set(xlabel="Model time (s)",ylabel="Integral contribution (%)",title="Integrator behavior under the same plant model");ax.legend();plt.tight_layout();plt.show()
4. Reconcile bounded outputs and export a model audit
Check finite states and actuator bounds, then report integral absolute error in °C·s. Model gains, numerical step and disturbance must accompany the trace. Passing this simulation is not hardware validation or a universal controller tuning rule.
assert np.isfinite(controlled["temperature_degC"]).all()
assert np.all((controlled["control_percent"]>=0)&(controlled["control_percent"]<=100))
error_integral=float(np.trapezoid(abs(controlled["error_degC"]),controlled["t"]))
assert np.max(controlled["unclamped_percent"])>100
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
pd.DataFrame({"time_s":controlled["t"],"temperature_degC":controlled["temperature_degC"],"applied_duty_percent":controlled["control_percent"],"requested_duty_percent":controlled["unclamped_percent"],"integral_percent":controlled["integral_percent"],"error_degC":controlled["error_degC"]}).to_csv(OUTPUT_DIR/"thermal-controller-simulation.csv",index=False)
(OUTPUT_DIR/"declared-control-model.json").write_text(json.dumps({**arguments,"antiwindup":"conditional integration","disturbance_start_s":300,"disturbance_equilibrium_offset_degC":-10,"integral_absolute_error_degC_s":error_integral,"hardware_validation":False},indent=2))
print("Bounded-control checks passed; simulated IAE (°C·s):",error_integral)
Bounded-control checks passed; simulated IAE (°C·s): 794.0295016691224
Exercises and reference explanation
- Explain derivative-on-measurement units.
- What does anti-windup prevent in this declared model?
Reference explanation. The derivative term uses temperature rate and Kd to produce duty units. Conditional integration prevents an integral state from pushing farther into actuator saturation; it does not prove hardware suitability.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Using NumPy/Pandas/Matplotlib, compare two declared Ki values with identical plant and actuator limits.
Prompt 2
With the same packages, reduce maximum duty and explain whether the disturbed setpoint is achievable under the model.
Human acceptance. Applied duty stays within limits, states remain finite and plant/controller/time-step assumptions accompany every result.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://docs.scipy.org/doc/scipy/tutorial/linalg.html',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
Sources, snapshots and next steps
IAPWS saturation reference release — IAPWS SR1-86(1992): unrestricted publication allowed in all countries (cover page); retrieved 2026-10-09. SHA-256
a10d38d73339b0713c86cfac88f9adcba137202167ccdfa1cd58238c438361d7. Reference equations and numerical coefficients transcribed from release pages 2–3; pedagogical model values are not new experiments.Authored numerical teaching model parameters — Original teaching parameters/code MIT; reference equation coefficients attributed to IAPWS SR1-86(1992); retrieved 2026-10-09. SHA-256
4350c4305a83c772fa295a8be0ef86a247f2bfaae930a081720d3cd0c753abe0. Declared parameters and seed; generated trajectories and calibration observations must be labelled simulated.
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.