检查有界热系统 PID 模拟

在明确假定的对象模型中，模拟测量微分控制、执行器饱和及条件积分。

工程 · 第 06 课🔋 离线核心分析⏱ 60 分钟数值教学模型📦 numpy · pandas · matplotlib · IPython

目标与运行方法

本课适合具备 Python 基础的本科生。先阅读每一步的方法与图注，再修改参数。AI 练习为可选部分，不需要 AI 账号即可完成核心课程。数据写在 notebook metadata 中，不需要另下载数据文件。请保持 notebook 已保存到当前工作文件夹。

学习任务是理解本课的研究问题、执行透明的计算、校验结果，并说明结论的边界。

环境与参数

In [1] · Python

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

读取内嵌快照

In [2] · Python

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']

方法函数

In [3] · Python

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

数据与模型边界

建筑记录包含来源中单栋建筑的真实十分钟能耗及环境观测。机场天气字段含来源插值，不是同位置校准参考。梁、电路和 PID 例子是明确的数值设计模型；计算响应不能证明硬件性能或实际验收。

分步分析

1. 声明对象及控制器单位

假定热对象向环境温度+增益·加热占空比+声明扰动松弛。加热限制在 0–100%。Kp 为 %/°C，Ki 为 %/(°C·s)，Kd 为 %·s/°C；微分作用于测量温度，不作用于参考阶跃。

In [4] · Python

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. 比较请求与实际施加控制

未限幅请求可超过执行器范围，而施加占空比有界。这是模拟执行器，不指令真实硬件。饱和属于声明模型，应展示而非隐藏。

In [5] · Python

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. 检查积分状态

当已饱和请求和误差会继续加深饱和时，条件积分保持积分器状态。积分贡献单位为占空比，是控制器状态，不是温度或储存热量。

In [6] · Python

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. 核对有界输出并导出模型审计

检查有限状态及执行器边界，并报告 °C·s 的绝对误差积分。模型增益、数值步长及扰动随轨迹保存。模拟通过不等于硬件验证，也不产生普适调参规则。

In [7] · Python

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

自主练习与参考答案

解释测量微分项的单位。
抗积分饱和在此声明模型中防止什么？

参考解释。 温度变化率乘 Kd 得到占空比单位。条件积分防止积分状态继续推向执行器饱和，不证明硬件适用。

先保留当前 notebook 副本，再修改参数。把新图与原图一起比较，并记录改变了哪些输入、哪些计算和哪些结论。

可选的 AI 编程练习

提示词 1

仅用 NumPy、Pandas、Matplotlib 在相同对象和限制下比较两个声明 Ki。

提示词 2

使用同样的包降低最大占空比，并解释扰动下设定值在模型中是否可达。

人工验收。 施加占空比有界，状态有限，每项结果随同对象/控制器/步长假设。

检查修改后的源代码，再手动运行受影响单元及后续校验。要求助手解释修改的目的、使用的包和数据来源，并用实际输出核对解释。

可选联网扩展

In [8] · Python

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

来源、快照与下一步

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.

下一步。 在保留本课的数据检查、参数记录和结果校验之后，将相同方法用于自己的公开或已获授权数据。记录研究设计与方法限制，不把示例结果当成新实验的验证。完整数据许可与生成说明见 cookbook 网页。