工业具身智能的第一步：数据建模、系统辨识与控制策略

0开篇必读：工业具身智能落地的三道现实坎

在进入技术细节之前，先把三件最容易被忽略、也最常让项目中途夭折的事摆在前面。工业具身智能——把具备感知、决策、执行能力的智能体真正装到产线上——最大的困难通常不在算法本身，而在下面三件事。本章不谈算法，只谈立项前必须准备什么；后面第 1~18 章的方法论，全部是在承认这三道坎的前提下设计的。

0.1 坎一：仿真到真实（sim2real）的差距——沙盒里考满分，真机上不一定及格

智能策略（无论是强化学习还是 MPC 的最优轨迹）都是在「辨识模型沙盒」里算出来的，而模型与真实装置之间永远有差距：

后果：沙盒里节能 8% 的策略，上真机可能只有 3%，且偶发在训练中从未出现过的越界动作。

工程对策（按重要性排序）

1. Domain Randomization（域随机化）：训练 / 求解时随机扰动模型参数（±10~20%）、噪声水平、滞后拍数，逼策略对「一族环境」鲁棒，而不是对单个模型过拟合；
2. 残差 / 不确定性在线监测：记录模型预测误差，实时估计可信度；可信度低时降权或冻结优化（与 §12 工况包络监测联动）；
3. 保守策略：目标函数中加入 CV 越界概率约束或风险厌恶项（CVaR）；强化学习侧直接采用离线 RL（只从已有历史数据学策略，不做在线试错）；
4. 小步长 + 安全屏蔽：动作限制在 CV 量程的 ±(2~5%)，任何动作都必须过安全层（RL 见 safe_action()，MPC 见 §13.2.3 约束软化）与 §15 四道闸门；
5. 低频在线更新：真机数据回流后按周 / 月级重辨识并复算策略，而不是每分钟在线学习；
6. 反馈兜底永不撤销：上层只输出设定值建议，底层 PID 闭环与安全联锁永远在位。

0.2 坎二：仿真平台并不是总存在的——这是比 sim2real 更前置的困难

学术论文常常默认「有一个高保真仿真器可直接用」。工业现场的现实恰恰相反：

• 没有现成数字孪生：老产线往往只有一张工艺流程图和几块电表，设备图纸、物性参数、阀门特性曲线缺失；
• 机理仿真成本极高：搭建可信的机理-数据混合模型需要工艺、设备、控制三方长期投入，往往比算法本身贵得多，且换一台设备就要重建；
• 每条产线都是定制：钢铁、水泥、化工各工序差异巨大，一个厂的孪生很难直接复用到下一个厂，边际成本高；
• 数据基础薄弱：没有历史数据，或数据质量差（量程错、单位错、时间戳错位、缺测严重）——连数据驱动辨识都做不了，更谈不上训练环境。

没有仿真条件时的渐进路线（按投入产出排序）

1. 先补测点与数据：按 §5 做零侵入测点规划，按 §6 做低扰动 PRBS 激励实验——这一步同时解决「没孪生」和「没数据」；
2. 先解静态子问题：负荷分配、台数启停、供水温度这类静态优化，用 §13.3.2 的 PSO 就足够，只需要静态效率曲线（短期实验即可测得），不需要完整动态仿真；
3. 先单点、后全局：先给最关键的一台设备建「单点沙盒」，验证闭环再向上扩展，避免一上来就要建全厂孪生；
4. 先用监督优化 / MPC，再考虑 RL：数据不足、模型不敢轻信的阶段，优先用 §13 方案 2/3（增益调度 PID、MPC），它们对模型误差容忍度更高、行为更易解释；等模型质量足以支撑大量离线回放，再引入 RL。

0.3 坎三：效果验证闭环——生产线上不允许「对照组」

学术界可以做 A/B 实验，工厂里通常不行：同一条产线不能一半时间开优化、一半时间不开（产量与安全责任不允许随意切换），而「如果不上系统本会用多少电」是测不到的。

对策就是第 17 章的 M&V 反事实基线：用投运前冻结的「工况 → 能耗」回归模型推算基准能耗，再与实测对比；配合影子模式（策略在线运行但只记录不下发）积累同工况下的对照数据，作为补充证据。

1项目概述：零侵入节能减碳改造

传统节能改造往往涉及硬件更换、管路改造、PLC 程序修改，存在停机风险、安全风险、验收周期长。

具体做法：在不侵入原有设备控制逻辑前提下，通过只读采集现有传感器数据，必要时少量外挂低成本传感器；采用低扰动开环辨识，建立被控对象输入-输出动态模型；基于辨识模型做上层设定值优化，寻找满足工艺约束前提下能耗最小的工作点；将最优设定值推荐下发到底层原有控制器，底层原有 PID 保持闭环、保留原联锁保护，随时可以切回原始模式。

两层控制架构（全文最重要的图）：

┌─────────────────────────────────────────────────────────────┐
│  上层优化层（本次新增 · 零侵入 · 可随时切除）                    │
│                                                              │
│   只读采集数据 → 辨识模型 → 滚动优化 → 输出「CV 设定值建议」     │
│        ▲                                      │              │
│        │ 只读                                 ▼ 只写设定值     │
├────────┼──────────────────────────────────────┼──────────────┤
│  底层控制层（原有 PLC/PID · 完全不改动）          │              │
│                                                              │
│   接收设定值 → PID 闭环运算 → 驱动执行器（阀/泵/电源）           │
│   安全联锁、报警、急停 —— 全部保留在原系统内                     │
├─────────────────────────────────────────────────────────────┤
│  现场设备层：传感器（温度/压力/流量） + 执行器                   │
└─────────────────────────────────────────────────────────────┘

为什么这样能节能？本章是全篇的理论基点——后面所有的辨识、建模、优化工作，都是为了回答这一个问题服务。

1.1 一句话答案

这 10%~30% 不是某个神奇算法凭空变出来的，它本来就摆在那里，来自四个相互独立的来源。判断是否值得做改造，就是看这四条在对象上命中了几条。

1.2 来源一：为了「不出事」而层层加码的冗余 ⭐

这是现场最常见、量级最大的一类浪费。为了避免出问题，企业的习惯做法是在每一个环节都多给一点：容量选型按最大工况再乘系数、并联设备常年全开互为备用、设定值留足温度余量、冷与热同时备着。冗余本身没有错——安全必须有冗余。真正的问题在于：冗余一旦加进去就变成「常态」，从此再没有人回头区分「备用能力」和「日常投放量」。

1.3 来源二：把「不等式约束」当成「等式」去跟踪

这是上一节在控制层面的机理。工艺的真实要求是一条不等式，而 PID 天生只能解等式：

工艺的真实要求（不等式，可行解有无穷多个）：
  y ≤ 30.0 ℃          ← 30.0、27.5、24.0、22.0 ... 全都合格

底层 PID 实际在解的（等式，解唯一）：
  y = r               ← r 是人工给定的常数，例如 24.0 ℃

差出来的那一段： r 与 30.0 之间的余量
              = 每天 24 小时、每年 8760 小时一直在支付的溢价

而 r 这个常数，是多层保守叠加的产物，每一层单独看都很合理：

关键在于：没有人定期复核这个数。它一旦写进设定值就永久固化，直到设备报废。上层优化做的正是把这条不等式重新请回来当作不等式处理——最优解永远落在边界上，而不是落在边界内侧几度的某个「保险」位置上。

1.4 来源三：驱动势的边际代价递增（物理层面的根）

再往下追一层，为什么「多留一点余量」会贵得不成比例？因为制冷、输送、换热这些环节的能耗，对推动过程是超线性的：

1.5 来源四：多设备耦合，单回路 PID 看不见

• 并联机组的出力组合：同样总流量，「两台 50%」与「一台 100%」效率差异显著；
• 跨设备负耦合：冷冻水供水温度 ↑ 省了冷机，却使末端流量必须 ↑ 费了泵——系统最优解不在任何单回路的最优解上；
• 内部抵消：制冷的同时再加热、加湿的同时除湿，两笔能耗互相报销；
• 首尾耦合：冷却塔出水温度与冷机冷凝温度互相制约。

这不是 PID 的缺陷，而是单回路结构的固有边界：每个回路的反馈信息只有自己的那一个偏差，它在原理上就看不到「自己多压这一度」与「隔壁设备多耗那几度电」之间的关联。

1.6 为什么底层 PID 永远不会主动去做

1. 设定值是人给的，而人的设定值就是 §1.2 那串保守叠加的结果——PID 对此无从质疑，它只会忠实执行；
2. PID 没有目标函数，只有误差。它的世界只有 e = r − y，不知道「这一度是花了多少电买来的」；
3. PID 没有约束处理能力。在多约束下求极值是 MPC / RTO 的专长，PID 只能靠输出限幅、抗积分饱和被动应对；
4. 积分项会把工作点牢牢钉死在设定值上——这本身就是「不愿放松」的机制。

1.7 上层优化层实际在做什么

上层优化层（分钟级 · 本次新增）                      底层控制层（毫秒~秒级 · 原有 PID）
  读当前工况 DV（气温 / 产量 / 来料）                   接收上层给出的 CV 设定值
  + 辨识模型 ĝ(u, d) 与能耗模型 P(u)                    闭环运算 → 驱动阀 / 泵 / 电源
  求解   min P(u)   s.t.  y ∈ [y_low + δ, y_high − δ]   抑制扰动、处理全部动态与安全联锁
    ↓                                                          ▲
  输出「最优 CV 设定值」 ────────── 只写设定值 ─────────────────┘

• 它输出的是设定值，不是阀门开度。这一点决定了它的安全性：底层设备、程序、联锁一行不改；上层失效或失准，只要把设定值撤回到原值，系统立刻回到原来的运行方式（§15 硬旁路）；
• 数学模型——在一个被工艺约束围成的可行域内寻找能耗最低点（完整形式见 §13.1）：
  min_u P(u)    s.t.    y_ss = g(u, d) ∈ [y_low + δ, y_high − δ] ,    u ∈ [u_min, u_max]
  由于可行域随环境温度、产量等 DV 实时移动，最优点也必须滚动地重算；
• 所需模型 g(u, d) 与 P(u) 从哪来？正是本文后半段的全部内容：低扰动 PRBS 激励 → Hammerstein (FOPDT) 辨识 → 残差验证，得到的就是 §13.1 能耗目标函数里那两个函数。没有模型，就只能写「温度 ≥ 50℃ 就开泵」这类靠人拍阈值的规则；有了模型，才能连续寻优。

1.8 节能率的量级、适用边界与核算陷阱

前面那个 10%~30% 就是下表几项的叠加（不等于简单相加，多项往往同源）：

什么时候不成立，必须如实告知客户：

• 设备已在额定高效区附近运行，或负荷率本就极低/极高，没有压缩空间；
• 负荷波动快于优化周期——上层跟不上反而引起设定值振荡，比不优化更费电；
• 工艺约束其实是等式（必须精确控温的工段），此时没有松弛可榨，只能靠更换高效设备本身；
• 计量陷阱：节能率必须建立在同工况、同产量、同气象的可比基线上，否则极易把「今年夏天凉快」记成「我优化了」。这部分单独成章，见 §17 节能量核算（M&V）。

2适用边界与约束条件

✅ 适合本方案

• 热工类对象（冷却、加热、恒温，空压机冷却水、工艺水冷、热风烘干等），存在惯性 + 纯滞后，执行器存在非线性；
• 工艺允许小幅度、短时间扰动（可做低幅值 PRBS 激励）；
• 有可观测被控变量（温度、压力等工艺指标）；
• 可采集操纵变量（阀位、频率、电流、功率），以及可测扰动（环境温度、进水温度、流量）；
• 存在明显「过供给」痕迹（投入产出比最高的一类对象，见 §1.2）：设定值层层加裕度、多台并联设备常年同开、主要靠阀门/挡板节流调节、存在再热与冷热抵消、低负荷或夜间不降载。冗余越多，可压缩空间越大。

❌ 不适用

• 工艺完全不允许任何扰动，无法开展激励实验，只能使用历史闭环数据（需要闭环辨识，本方案不作为首选）；
• 强耦合 MIMO、多变量强交互，单回路 FOPDT 假设失效；
• 频繁启停、间歇批次生产，工况持续跳变，系统强时变；
• 安全联锁和工艺不允许外部设定值写入。

3术语与缩写对照表（含大白话解释）

4系统辨识整体方法论框架（迭代闭环流程）

现场调研与变量筛选 → 测点规划（零侵入，只读采集）
→ 激励方案设计（低扰动 PRBS，控制扰动幅度，保障生产）
→ 时序数据采集 → 数据预处理 → PE 检验（不通过 → 返回激励设计）
→ 参数辨识 NLS，Hammerstein (FOPDT)
→ 模型验证（拟合指标，ACF 残差，CCF 输入-残差交叉相关，AIC/BIC 择优）
→ 模型诊断（参数共线性，置信区间，插值/外推边界评估，假设校验）
→ 控制策略设计：MPC / 增益调度 PID / 模型辅助 PID 整定
→ 离线仿真验证（数字孪生，先在模型上跑闭环，评估超调、约束、能耗）
→ 软投切：影子模式 → 手动推荐 → 半自动 → 自动设定值下发（保留随时切回）
→ 在线模型监测，模型漂移告警 → 必要时重新激励、重新辨识，迭代更新模型

5现场调研、变量筛选与测点规划（零侵入测点选型）

5.1 变量分类

1. 操纵变量 MV：上层优化可以推荐调整的变量，如制冷机电流、泵频率、阀位；
2. 被控变量 CV：工艺约束指标（温度、压力，必须落在工艺上下限）；
3. 可测扰动 DV：不受控制，但可观测，影响 CV；环境温度、进水温度、来料温度。

5.2 测点规范

• 时间同步：全部测点统一时间戳（NTP 对时，误差 < 1 个采样周期）；不同源的采集程序必须打「源时间戳」，防止用「到达时间」冒充「测量时间」——无线/总线传感器（如 LoRa）的真实测量时刻可能比平台收到时刻早十几秒，直接用接收时间等于给系统人为加了一段随机滞后，辨识出的动态必然失真（本文 TEC 案例就踩过此坑，见 §19.6）；
• 采样周期：推荐 1~10s，按对象带宽确定。工程经验：采样周期取主导时间常数的 1/4 ~ 1/10。例如冷却水系统时间常数约 10 min，采样 30~60s 足够；风口级快响应约 1 min，采样 5~15s。「高于系统带宽 2 倍」只是奈奎斯特不混叠的理论下限，按它采样只能保证「不走样」，远不足以「看清动态形状」；
• 采集模式：旁路采集，只读；原系统不受采集程序影响；
• 记录元数据：设备状态（运行/停机）、事件日志（清洗、检修、批次切换），用于后续标记异常段。

📌 输出成果：变量清单，测点映射，采样周期，工艺安全上下限。

6实验设计，激励信号设计（最小扰动激励，兼顾生产安全）

6.1 什么是 PRBS？为什么要用它？

激励信号设计要点（多电平 PRBS，不是传统 2 电平）：

1. 激励幅值：在当前稳态点 ±Δ 范围，Δ 由工艺工程师确认，保证不越工艺阈值。同时有下限：输出响应的峰峰值应 ≥ 3~5 倍传感器噪声标准差，否则响应淹没在噪声里，实验白做。设计方法：先做 1~2 次小幅「试探性阶跃」，测出响应幅度与噪声水平，再定正式幅值；
2. 电平数量：一般取 3~5 电平；电平越多激励越丰富，对非线性辨识越有利，但扰动变大；
3. 保持时间（驻留时间）：取主导时间常数的 1/3 ~ 1 倍。热工对象时间常数常达 10 min 以上，若要求「驻留大于时间常数」，一次实验要跑数天，且驻留过长会使激励能量全堆在低频、丧失中高频信息；驻留太短（远小于时间常数）则响应起不来。1/3~1 倍是兼顾两者的工程折中；
4. 实验总时长：≥ 10~20 倍主导时间常数。时间常数 10 min 的系统，单次实验至少 2~3 小时。数据不够是辨识失败最常见的隐性原因；
5. 激励独立性：激励序列预先生成，不根据 CV 实时反馈调整。

6.2 重要澄清：「开环辨识」在工业现场的准确含义

工业语境下常说「开环激励 = 底层 PID 仍维持闭环，设定值由外部激励序列给出」。严格地说，这属于「闭环操作下的设定点激励」，不是教科书意义的开环辨识——回路并没有断开，底层 PID 还在工作。必须澄清两点，否则评审与落地都会出问题：

1. 为什么这样做仍然有效：闭环辨识产生偏差的根源是「输入与噪声通过反馈相关联」。而当激励是外部注入的设定值序列、且与过程噪声统计独立时，只要激励功率足够盖过噪声经反馈引起的输入分量，预测误差类方法（NLS/ARX）的偏差可以忽略。前提：① 激励序列离线生成、与扰动独立；② 激励幅值满足 §6.1 的信噪比下限；
2. 必须用实测 MV，不能用激励设定序列：底层 PID 会「抵抗」设定值激励（它要维持原 CV，会反向调节 MV），导致实际 MV 曲线 ≠ 激励序列。辨识的输入必须取执行器上实际测到的 MV（电流/阀位/频率实测值），否则模型增益会系统性偏小。本文 TEC 案例即按此执行（实测电流作为 MV，见 §19）。

7时序数据采集、数据预处理与数据质量诊断

处理顺序（注意：先划分、后归一化）：

异常标记 → 坏点剔除 → 缺失处理 → 稳态/非稳态分段 → 数据集划分（先划分！）
→ 去趋势，构造偏差变量 → 归一化（只用训练集统计量！）→ 输出辨识数据包

1. 异常点：传感器尖峰、通信丢包、设备启停、人工干预事件；直接标记，优先剔除。连续多点缺失，直接舍弃片段，禁止长段插值（短缺口 ≤ 3 个采样点可线性插值并打标记）；
2. 稳态/非稳态分段：实验初始段若系统未达稳态必须剔除——初始段的大幅自然趋势会淹没激励响应、主导参数估计（TEC 案例实测：水杯初始 41.9℃ 自然降温的前 16 min 混入后，测试段 R² 从 0.95 掉到 0.91，见 §19.6）；
3. 数据集划分（时序切分，禁止打乱时序）——必须先于归一化：训练集 70%，验证集 20%，测试集 10%，仅作为插值初步评估；黄金验证：独立第二条激励轨迹作为外测试集，评估泛化能力；
4. 去趋势：把所有变量转为相对于稳态工作点的增量，消除稳态偏移，适配动态辨识；
5. 归一化：不同量纲（℃，A，kW，Hz）映射至统一区间，改善 NLS 优化收敛；预测输出后执行反归一化。

import numpy as np

def reject_outliers(x, k=5.0):
    """稳健剔尖峰：MAD 判据 |x-median| > k*1.4826*MAD → NaN。
    为什么不用 3σ：尖峰会同时抬高 σ 本身，导致漏检；MAD 对尖峰不敏感。"""
    x = np.asarray(x, float).copy()
    mad = np.nanmedian(np.abs(x - np.nanmedian(x)))
    if mad > 0:
        bad = np.abs(x - np.nanmedian(x)) > k * 1.4826 * mad
        x[bad] = np.nan                      # 统一用 NaN 标记，后续统一处理
    return x

def fill_short_gaps(x, max_gap=3):
    """短缺口线性插值（≤3 点）；长缺口保持 NaN——禁止长段插值伪造动态。"""
    x = np.asarray(x, float).copy()
    for i in range(len(x)):
        if np.isnan(x[i]):
            j = i
            while j < len(x) and np.isnan(x[j]):
                j += 1
            if 0 < i and j < len(x) and (j - i) <= max_gap:   # 两端有真值且缺口短
                x[i:j] = np.interp(np.arange(i, j), [i - 1, j], [x[i - 1], x[j]])
    return x

def resample_align(t_arrive, x, age, dt):
    """时间轴还原 + 等间隔重采样（无线/总线测点必做）。
    t_true = t_arrive - age（age = 读数年龄），禁止直接用到达时刻。"""
    t_true = np.asarray(t_arrive, float) - np.asarray(age, float)
    grid = np.arange(t_true[0], t_true[-1], dt)
    return grid, np.interp(grid, t_true, x)

def split_then_scale(data, ratios=(0.7, 0.2, 0.1)):
    """先切分、后归一化——归一化统计量只来自训练集（防数据泄漏）。"""
    n = len(data)
    n_tr, n_va = int(n * ratios[0]), int(n * ratios[1])
    train, val, test = data[:n_tr], data[n_tr:n_tr + n_va], data[n_tr + n_va:]
    mu, sd = train.mean(0), train.std(0) + 1e-12
    return [(d - mu) / sd for d in (train, val, test)], (mu, sd)

8持续激励 PE 检验（辨识准入关卡）

8.1 检验对象（Hammerstein 特别注意）

PE 检验不是只检验原始 MV；而是检验经过静态非线性映射后的隐变量 z。就算输入电流是 PRBS，如果一直落在 TEC 非线性曲线的平坦区域，隐变量 z 几乎不变，此时内部 FOPDT 子系统不可辨识，参数多解，拟合看着好看，但外推预测完全失效。

8.2 检验方法（具体构造）

1. 取滞后阶数 n（n ≥ 待辨识模型的滞后参数个数，工程上取 10 左右）；
2. 构造回归向量 φ(k) = [u(k-1), u(k-2), ..., u(k-n)]^T，去均值后计算信息矩阵：

3. 判定（两条同时满足）：满秩 rank(R) = n，说明输入是 n 阶持续激励；条件数有界 κ(R) = λmax/λmin 小于阈值。阈值是经验值：设计良好的 m 序列 κ 通常在 10 以内；κ > 10³ 判定病态、激励不足。注意：含常数项（不去均值）的条件数会被信号均值放大，判据必须用去均值后的条件数；
4. 对 Hammerstein：用预估计的静态非线性得到隐变量 ẑ，对 ẑ 重复上述检验。

参考实现（采集脚本内置的 PE 自检）：

def check_pe(u, order=10, cond_limit=1000.0):
    # 持续激励检验：回归向量 φ=[u(k-1)..u(k-order)] 去均值后的信息矩阵条件数
    n = len(u)
    rows = [[u[k - i] for i in range(1, order + 1)] for k in range(order, n)]
    mean = [sum(r[j] for r in rows) / len(rows) for j in range(order)]
    R = [[sum((r[i]-mean[i]) * (r[j]-mean[j]) for r in rows) / len(rows)
          for j in range(order)] for i in range(order)]
    return cond_number(R) < cond_limit   # cond_number：幂迭代求 λmax/λmin

9模型结构选择：Hammerstein (FOPDT) 模型定义、备选模型对比

9.1 模型定义

Hammerstein 结构（静态非线性块在前）：

f(・) 是静态非线性映射（多项式/样条，系数待辨识），输入操纵变量（可辅以功率），输出等效热驱动隐变量 z(k)。

FOPDT 连续传递函数：

K 增益，T 时间常数，τ 纯滞后；离散化到采样周期 Ts，叠加扰动项：

τ_d 是离散滞后采样点数；Δd_i(k) 为可测扰动增量。

结构示意：

 电流/阀位 u ──►┌──────────────┐  隐变量 z  ┌───────────────┐
                │ 静态非线性 f  │ ─────────► │ FOPDT 线性动态 │ ──┐
 功率 P ───────►│ (多项式/样条) │            │ K/(Ts+1)·e^-τs │   │
                └──────────────┘            └───────────────┘   ▼
 环境温度等 DV ─────────────────────────────────────────────► ( + ) ──► 输出 y（温度）

9.2 ⚠️ 增益不可辨识性（教科书级的坑）

Hammerstein 模型有一个必须处理的问题：静态块 f 乘以任意常数 α、同时线性块增益 K 除以 α，输入-输出行为完全不变。也就是说 f 与 K 的增益分配是任意的，NLS 会遇到无穷多组等价解，表现为：参数置信区间巨大、不同初值收敛到看似不同但预测相同的参数、优化病态。

9.3 备选模型对比

✅ 阶次扫描原则：如果 FOPDT 的残差 CCF 显著不为 0，说明一阶假设不足，可以升级到高阶线性模型；优先保持 FOPDT，一阶模型可解释性强，更适合工程落地。

10参数估计（NLS 非线性最小二乘）

损失函数：最小化模型自由仿真输出与实测被控变量的残差平方和：

θ 包含：静态非线性 f 的系数（已做增益归一化）、FOPDT 参数 K, T、扰动增益、离散滞后 τ_d。

10.1 ⚠️ 离散滞后的处理

「θ 包含离散滞后 τ_d，用 LM 求解」这句话有一个算法陷阱：τ_d 是整数参数，损失函数对它不连续、不可导，LM 等梯度类算法无法直接优化，会卡在最近邻的整数上。

正确做法——双层优化：

外层：枚举 τ_d ∈ {0, 1, 2, ..., τ_max}      ← 网格/枚举（τ_max 按物理先验，如 滞后上限/Ts）
  内层：固定 τ_d，用 LM 优化连续参数 {f 系数, K, T, 扰动增益}   ← 可导，正常 NLS
取外层损失最小的组合

参考骨架：

best = None
for tau_d in range(0, tau_max + 1):            # 外层：枚举整数滞后
    for seed in range(N_INIT):                 # 多初值防局部极小
        theta0 = random_init(physical_bounds)
        theta, loss = levenberg_marquardt(loss_fn, theta0, fixed_tau_d=tau_d)
        if best is None or loss < best.loss:
            best = (tau_d, theta, loss)

10.2 工程要点

1. 多初值随机搜索：NLS 是非凸优化，极易卡在局部极小；多组随机初值 + 物理边界约束，筛选最优解；
2. 参数硬约束：物理约束，例如时间常数 T > 0、滞后 ≥ 0、静态块单调性（TEC 电流增大冷量必须单调不减），拒绝无物理意义参数；
3. 优化器推荐：Levenberg–Marquardt（LM），中小规模问题收敛快；大规模可用信赖域（trust-region）；
4. 模型假设：NLS 默认输入测量无噪声，输出为加性白噪声。

11模型验证、残差检验、模型择优（AIC/BIC，ACF+CCF）

❗ 不能只看拟合曲线，必须全套验证：

1. 时域量化指标：NRMSE、R²（在独立测试段 / 黄金验证轨迹上计算）：

2. 残差检验（核心），残差 ε(k) = y(k) − ŷ(k)：

3. AIC/BIC 信息准则，多个候选模型择优：

n_p：参数数量；BIC 对参数数量惩罚更强，工程优先使用；AIC/BIC 数值越小越好。

✅ 判定原则：AIC/BIC 不能单独作为判定依据，必须配合 ACF+CCF 残差检验一起评估。

参考实现（numpy）：验证指标与三项残差检验

def nrmse(y, yhat):
    """归一化均方根误差：<5% 优，5~10% 可接受，>15% 不可用"""
    return 100.0 * np.sqrt(np.mean((y - yhat) ** 2)) / (y.max() - y.min())

def r2_score(y, yhat):
    return 1.0 - np.sum((y - yhat) ** 2) / np.sum((y - y.mean()) ** 2)

def acf(e, lags=20):
    """残差自相关（去均值后归一化）。95% 置信带 ≈ ±1.96/√N，
    显著超出说明还有未建模动态（如一阶假设不够、缺滞后）。"""
    e = np.asarray(e) - np.mean(e)
    denom = np.dot(e, e)
    return np.array([np.dot(e[:-l], e[l:]) / denom for l in range(1, lags + 1)])

def ccf(e, u, lags=20):
    """残差与历史输入互相关：显著不为 0 ⇒ 结构性偏差（不是噪声问题，需改模型结构）"""
    e = np.asarray(e) - np.mean(e); u = np.asarray(u) - np.mean(u)
    return np.array([np.dot(e[l:], u[:-l]) / np.sqrt(np.dot(e, e) * np.dot(u, u))
                     for l in range(lags)])

def ljung_box_q(e, lags=20):
    """Ljung-Box 白噪声整体检验统计量：Q = N(N+2)·Σ r_k²/(N-k)。
    自由度 = lags − 模型参数个数；Q 超过卡方临界值（scipy.stats.chi2.ppf(0.95, dof)）
    ⇒ 残差非白噪声，模型结构还需改进。"""
    n = len(e); r = acf(e, lags)
    return float(n * (n + 2) * np.sum(r ** 2 / (n - np.arange(1, lags + 1))))

def aic_bic(sse, n, n_p):
    """修正后的信息准则（含除以 N，见上文公式）"""
    aic = n * np.log(sse / n) + 2 * n_p
    bic = n * np.log(sse / n) + n_p * np.log(n)
    return aic, bic

12模型诊断、假设校验、不确定性评估与风险清单

1. 参数可辨识性 & 参数共线性：热工模型中时间常数 T、增益 K、滞后 τ 经常强耦合。如果目标只是预测 / 上层 MPC 优化：参数共线性可以容忍，只要预测不确定性低；如果目标是提取物理热阻热容：共线性是致命问题，需要 Fisher 信息矩阵计算 95% 置信区间；置信区间跨 0 → 参数不可辨识，删减模型项。
2. 插值验证 vs 外推验证：插值是测试点落在训练多维工况凸包之内；外推是工况超出训练包络，Hammerstein 非线性模型外推风险很高。部署时做工况包络监测，轻度越界先降权+告警（预测可信度标记为低、优化动作冻结），持续越界或越界幅度大才切回原有控制逻辑，并设迟滞区间防止边界上反复横跳；
3. LTI 线性时不变假设检验：模型假设参数固定；使用滑动窗口分段辨识，持续监测参数随时间漂移（设备结垢、老化）。区分两个现象——非线性：参数随工况点变化（Hammerstein 模型就是用来描述这类）；时变：参数随时间缓慢漂移，和当前工况无关。漂移显著 → 模型失效，触发告警，评估是否需要重新辨识；
4. 采样混叠、时变滞后：采样周期按 §5.2 的 1/4~1/10 时间常数经验选取；FOPDT 模型假设滞后固定；很多热工系统滞后会随工况改变（如流量变、传输延迟变），这是模型本身的固有局限，必要时分工况辨识多组参数。

13控制策略设计，方案对比与智能寻优（PSO / 强化学习）

方案 1：模型辅助 PID 整定（最轻量）

用辨识得到 Hammerstein (FOPDT) 模型，在稳态工作点局部线性化，得到 FOPDT 传递函数；基于 FOPDT，用 Ziegler-Nichols/Cohen-Coon 计算 PID 初始参数。

✅ 适用场景：产线不想新增复杂优化器；仅优化底层 PID 参数；❌ 局限：固定参数 PID 只在该工作点附近效果好；PID 只消除误差，不会自动寻找能耗最优工作点。

方案 2：增益调度 PID（折中方案，推荐中小规模）

1. 在多个稳态工作点，对 Hammerstein (FOPDT) 局部线性化，得到多组 FOPDT；
2. 每组 FOPDT 对应一组 PID 参数；
3. 在线运行时，基于当前工况（MV、DV）查表插值切换 PID 参数。

底层 PID 仍然是 PLC 原生 PID；上层只维护一张调度参数表，不做在线滚动优化；算力需求低，边缘设备轻松部署。❌ 局限：调度表是离线预先计算；只能逼近最优，不是显式能耗约束优化。

方案 3：非线性 MPC 上层监督优化（节能减碳首选，主推方案）

13.1 决策变量的澄清

零侵入约束下，上层根本没有权限直接写 MV（阀位/电流是底层 PID 的输出），所以上层的决策变量不能是 MV，而必须直接是 CV 设定值。因此实际架构是分层的稳态目标优化（SSTO）：

上层 SSTO（每 1~5 min 执行一次）：
  输入：当前 DV（环境温度、来料温度）、工艺约束、能耗模型（Hammerstein 静态块 + 功率模型）
  求解：min  稳态能耗 P(u)                     ← 目标：省电
        s.t. 稳态 CV = y_ss(u, d) ∈ [CV_low + δ, CV_high − δ]   ← 守住工艺红线（δ 为动态安全裕度）
             u ∈ [u_min, u_max]
  输出：最优稳态对应的 CV 设定值 → 下发给底层 PID
底层 PID（原有，毫秒~秒级）：跟踪设定值，处理全部动态与扰动

能耗目标函数的完整定义：

• 第一项：预测时域 H 内的累计能耗（功率 P 由静态非线性块给出，TEC 案例中 P = 1.05·I² + 4.25·I，见 §19.4）；
• 第二项：MV 变化惩罚（ρ > 0），防止设定值频繁抖动冲击工艺与设备寿命；
• 约束：MV 上下限、MV 变化速率限、CV 工艺约束（含安全裕度 δ）。

13.2 MPC 的原理、滚动时域逻辑与工程实现

§13.1 给出了「算什么」，§13.2 回答「怎么算、怎么调、算不出来怎么办」。下面七节是一份可直接照着实现的 MPC 工程手册，全部基于第 9~12 章辨识出来的 Hammerstein (FOPDT) 模型。

13.2.1 预测模型：Hammerstein (FOPDT) 的离散化与多步预测

MPC 的全部本事都建立在「能预测未来」上。把连续 FOPDT 按上层采样周期 T_s 离散化，得到工程上唯一需要递推的一条式子：

• a 是「记忆系数」：a 越接近 1（τ ≫ T_s），系统越慢，当前值对未来影响越大；
• 静态非线性只在输入端出现一次（Hammerstein 的定义）：先把 MV 过静态块得到隐变量 z，再进线性动态。TEC 案例中 z = f(u) = u/(1+0.12u)（降温饱和），功率 P(u) = 1.05u² + 4.25u（§19.4）——同一组静态块同时给出「能省多少」和「花多少电」；
• DV 前馈是白送的：可测扰动 d（气温、进水温度、产量）直接写进预测式，MPC 能在扰动真正影响到 CV 之前就提前动作；PID 只能在偏差出现之后被动纠正。这是 MPC 相对 PID 最实打实的一项优势；
• 多步预测 = 把上一步的预测值喂回自己（自由仿真，与 §10 同口径）：y_k+1 用实测 y[k]，y_k+2 用 ŷ_k+1，依此类推，得到 H 步轨迹；
• 纯滞后用「队列」实现：维护最近 n_k 拍已下发的输入，预测第 j 步时取队首（即 n_k 拍前那个）。θ/T_s 非整数时向上取整——宁可多算一拍滞后（保守），也不要低估滞后。

13.2.2 代价函数：把「省电」写成可执行的目标

13.2.3 约束：硬约束、软约束与安全裕度 δ（不可省）

• 为什么 CV 约束必须软化：模型必有误差、扰动必未测全、执行器可能暂时到不了位。若写成硬约束，一旦某拍「无可行解」，求解器要么报错停机、要么输出跳变的伪解——这两种行为上真机都是事故。软化后即使不可行，优化器也会给出「越界最小的那个折中解」，同时 ε > 0 本身就是一条现成的告警信号；
• δ（动态安全裕度）怎么取：δ ≥ 3σ_model，其中 σ_model = 模型自由仿真 NRMSE × CV 量程（§16 验收指标已给出 NRMSE）；再叠加底层 PID 可能的动态超调量。TEC 案例中量程约 2 ℃、NRMSE 6.5% → σ ≈ 0.13 ℃，取 δ = 0.3 ℃（约 2.3σ）并配合 §15 硬旁路；
• 零侵入特有的一条：上层下发的是设定值而非 MV，因此 CV 约束实际约束的是「给底层的设定值区间」，必须把底层 PID 跟踪过程中的超调量也算进 δ，否则底层一超调就顶到工艺红线。

13.2.4 求解：零侵入下的两层分解，与完整 MPC

形态 A｜SSTO（零侵入默认，推荐先上）
  每个优化周期：读 DV → 求 min P(u) s.t. y_ss(u,d) ∈ [y_low+δ, y_high−δ]
              → 得到 (u*, y*) → 把 CV 设定值 = y* 下发给底层 PID（按速率限做斜坡）
  动态交给底层 PID；上层只负责「把目标定在最省电的边界上」

形态 B｜完整 MPC（需 MV 写权限）
  每个优化周期：读当前 y → 用模型做 H 步自由预测 → 优化 M 步输入序列 U = {u_k..u_{k+M-1}}
              → 只下发 u_k → 下一周期重新测量、重新预测、重新优化（滚动时域）

求解器怎么选（按维数递增）

1. 一维 / 二维单调问题（最常见：单台设备的电流或供水温度）：直接用二分 / 黄金分割求「刚好守住约束」的最小 u——稳态模型单调时这是全局最优且毫秒级，不需要任何优化库（见 §13.2.6 的 ssto()）；
2. 多维凸、可求导：SQP（SciPy minimize / OSQP）或投影梯度；CV 约束写成惩罚项时，投影梯度十几行就够用；
3. 多维、非凸、含启停组合（导数不存在）：转 §13.3.2 的 PSO，或「PSO 定台数组合 + 连续量用梯度」的两层法；
4. 求解失败 / 超时：保持上一拍输出并告警——绝不空发、绝不发未收敛的中间解。这条要写进代码的第一行，而不是靠人记住。

13.2.5 参数整定清单（工程起点值）

13.2.6 参考实现（纯 numpy，可直接运行）

用 §19 辨识出的 TEC 模型（τ = 909 s、θ = 20 s、K = −0.98 ℃/A、静态块含饱和、功率 P = 1.05u² + 4.25u）跑一个 2 小时闭环：第 30 min 冷却水温度阶跃 +1 ℃，对比「操作员把水温死守在 26.0 ℃」与「MPC 顶住 26.2 ℃ 约束边界」两种做法。

import numpy as np

# ===== 辨识得到的模型（§19.4 静态块 + §19.5 慢通道）=====
TAU, THETA, TS = 909.0, 20.0, 60.0        # 时间常数 s、纯滞后 s、上层优化周期 s
K      = -0.98                            # 静态增益 ℃/A（电流↑ → 水温↓）
NK     = max(1, int(round(THETA / TS)))   # 离散滞后拍数（滞后量化到整数拍）
A      = float(np.exp(-TS / TAU))         # 离散 FOPDT 系数 a
D_GAIN, Y_FREE = 0.85, 27.2               # 扰动增益、完全不制冷时水温（DV=22 ℃ 实测）
U_MIN, U_MAX = 0.0, 5.0

static_block = lambda u: u / (1 + 0.12 * u)       # Hammerstein 静态非线性（饱和）
power        = lambda u: 1.05 * u**2 + 4.25 * u   # §19.4 实测功率模型 W
y_free       = lambda d: Y_FREE + D_GAIN * (d - 22.0)


def rollout(y0, u_past, u_plan, d_plan):
    """用 Hammerstein (FOPDT) 模型做 H 步自由预测（含纯滞后队列）"""
    H = len(d_plan)
    u_full = list(u_plan) + [u_plan[-1]] * (H - len(u_plan))   # 控制时域之后保持末值
    y, out = y0, []
    q = list(u_past[-NK:]) + u_full                   # 队首是 nk 拍前已下发的输入
    for k in range(H):
        y = A * y + (1 - A) * (y_free(d_plan[k]) + K * static_block(q[k]))
        out.append(y)
    return np.array(out)


def cost(u_plan, y0, u_past, d_plan, u_prev, y_lo, y_hi, delta, w_du=0.5, w_v=1e3):
    y = rollout(y0, u_past, u_plan, d_plan)
    e_hi = np.maximum(0.0, y - (y_hi - delta))        # 松弛变量：越界量（软约束）
    e_lo = np.maximum(0.0, (y_lo + delta) - y)
    du = np.diff(np.concatenate([[u_prev], u_plan]))  # 移动惩罚：防止设定值抖动
    return (np.sum(power(u_plan)) * TS / 3600.0       # 预测时域内累计能耗 Wh
            + w_du * float(np.sum(du ** 2))
            + w_v * float(np.sum(e_hi ** 2 + e_lo ** 2)))


def mpc_solve(y0, u_past, u_prev, d_plan, y_lo, y_hi, delta, M=5, iters=60, lr=0.15):
    """滚动时域优化：规划 H 步，只返回序列第一步"""
    u = np.clip(np.full(M, u_prev), U_MIN, U_MAX)
    for _ in range(iters):
        J0 = cost(u, y0, u_past, d_plan, u_prev, y_lo, y_hi, delta)
        g = np.zeros(M)
        for j in range(M):                            # 数值梯度（M 很小，代价可忽略）
            up = u.copy(); up[j] += 1e-4
            g[j] = (cost(up, y0, u_past, d_plan, u_prev, y_lo, y_hi, delta) - J0) / 1e-4
        u = np.clip(u - lr * g / (np.linalg.norm(g) + 1e-9), U_MIN, U_MAX)
    return float(np.clip(u[0], U_MIN, U_MAX))         # 只执行第一步


def ssto(d, y_hi, delta):
    """稳态目标优化 SSTO：二分求「刚好守住工艺红线」的最小电流"""
    lo, hi = U_MIN, U_MAX
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if y_free(d) + K * static_block(mid) <= y_hi - delta:
            hi = mid                                  # 满足约束 → 继续往下压
        else:
            lo = mid
    return hi


def hold_setpoint(sp, d):
    """底层 PID 的理想跟踪：求把 CV 稳到 sp 所需的电流（逆稳态模型）"""
    lo, hi = U_MIN, U_MAX
    for _ in range(60):
        mid = 0.5 * (lo + hi)
        if y_free(d) + K * static_block(mid) > sp:    # 还不够冷 → 加大电流
            lo = mid
        else:
            hi = mid
    return hi


def run(mode, n=120, y0=26.0):
    """2 小时闭环仿真：n 步 × 60 s"""
    d_seq = np.full(n, 22.0)
    d_seq[30:] += 1.0                                 # 第 30 min 冷却水温度阶跃 +1 ℃（DV）
    y, u_prev, u_hist, e_sum, over = y0, 0.0, [0.0] * NK, 0.0, -9.9
    Y_LO, Y_HI, DELTA, SP_FIX = 20.0, 26.5, 0.3, 26.0
    for k in range(n):
        d_plan = np.full(30, d_seq[k])                # 预测时域 H=30 步（30 min ≈ 2τ）
        if mode == "fix":                             # 现状：操作员把水温死守在 26.0 ℃
            u = hold_setpoint(SP_FIX, d_seq[k])
        else:                                         # MPC：滚动求最优，只执行第一步
            u = mpc_solve(y, u_hist, u_prev, d_plan, Y_LO, Y_HI, DELTA)
        y = A * y + (1 - A) * (y_free(d_seq[k]) + K * static_block(u_hist[-NK]))
        u_hist.append(u); u_prev = u
        e_sum += power(u) * TS / 3600.0
        over = max(over, y - Y_HI)                    # CV 相对工艺红线（26.5 ℃）的最大越界
    return e_sum, over, y


d0 = 22.0
u_opt = ssto(d0, 26.5, 0.3)
print("SSTO 最优工作点: 电流 %.3f A → 稳态水温 %.2f ℃（顶住 26.2 ℃ 约束边界）  功率 %.2f W"
      % (u_opt, y_free(d0) + K * static_block(u_opt), power(u_opt)))
u_fix = hold_setpoint(26.0, d0)
print("现状死守 26.0 ℃: 电流 %.3f A   功率 %.2f W  → 稳态节能率 %.1f%%"
      % (u_fix, power(u_fix), 100 * (power(u_fix) - power(u_opt)) / power(u_fix)))
for m in ("fix", "mpc"):
    wh, over, yend = run(m)
    print("%-4s 模式: 2h 累计能耗 %.2f Wh   末态水温 %.2f ℃   CV 相对工艺红线最大越界 %+.2f ℃"
          % (m, wh, yend, over))

实际运行结果（本仓库 numpy 环境实测）：

SSTO 最优工作点: 电流 1.163 A → 稳态水温 26.20 ℃（顶住 26.2 ℃ 约束边界）  功率 6.36 W
现状死守 26.0 ℃: 电流 1.435 A   功率 8.26 W  → 稳态节能率 23.0%
fix  模式: 2h 累计能耗 34.22 Wh   末态水温 26.00 ℃   CV 相对工艺红线最大越界 -0.42 ℃
mpc  模式: 2h 累计能耗 28.14 Wh   末态水温 26.20 ℃   CV 相对工艺红线最大越界 -0.25 ℃
→ 2 小时累计节能 17.8%，且 CV 全过程未触碰工艺红线（26.5 ℃）

13.2.7 常见失效模式与对策（上线前逐条自查）

13.3 模型之上的智能寻优：粒子群（PSO）与强化学习（RL）

13.3.1 为什么模型建好后还需要智能算法？

§13.1 的 SSTO/MPC 求解的是一个带约束的非线性规划。当目标具备下列特征时，梯度类求解器会吃力甚至失效：

1. 目标非凸、不连续：存在设备启停组合、档位切换、查表型工艺曲线（导数不存在）；
2. 变量多、耦合强：多台冷机负荷分配、多泵并联变频组合，可行组合数爆炸；
3. 含难以写成解析式的经验规则：老师傅的「这条线夏天不能低于 xx℃」这类约束。

此时有两条补充（或替代）路线：PSO 做慢周期静态寻优，强化学习做快周期序贯决策。两者的共同前提是——用辨识模型当沙盒评估，绝不允许在真机上「试着调一调看看」。

13.3.2 路线 A：粒子群优化（PSO）——黑箱目标也能搜

典型适用场景

标准 PSO 步骤：① 初始化 NP 个粒子（位置 = 待优化变量，速度 = 随机小值，均限制在上下限内）；② 用辨识模型仿真评估每个粒子的适应度 J(x)（含约束惩罚）；③ 更新个体最优 pbest 与全局最优 gbest；④ 更新速度与位置：

⑤ 迭代至收敛或达到上限，输出 gbest，经 §15 安全闸门后才允许下发。

工程推荐参数（默认起点，按现场微调）

• 粒子数 NP = 5~10 × 变量维度，典型 20~40；
• 惯性权重 w 由 0.9 线性递减至 0.4（先广搜后精搜）；学习因子 c₁ = c₂ ≈ 1.5~2.0；
• 迭代 50~200 代；收敛判据：gbest 连续 20 代改善 < 0.1% 即停止；
• 边界处理：截断 + 速度反向（比随机重置更稳定）；约束用惩罚项，CV 越界的惩罚系数要比能耗项高 1~2 个数量级（宁可不省电，也不能出工艺事故）；
• 每代评估要跑多工况 / 多扰动场景，避免「过拟合到单一工况的最优解」。

参考实现（纯 numpy，可直接运行）：两台冷机负荷分配——满足总冷量需求的前提下最小化总功率。

import numpy as np

# 两台冷机功率模型（凸二次，可由 §9 静态非线性块标定得到）
# P_i(u_i) = a_i*u_i^2 + b_i*u_i + c_i   u_i ∈ [u_min, 1]（负荷率）
A = np.array([120.0, 150.0])      # 二次项系数 kW
B = np.array([180.0, 160.0])      # 一次项系数 kW
C = np.array([ 15.0,  18.0])      # 空载功率 kW
CAP = np.array([500.0, 650.0])    # 满负荷冷量 kW

U_MIN, U_MAX = 0.2, 1.0           # 负荷率上下限（低于 0.2 效率骤降/喘振）
DEMAND = 800.0                    # 当前冷量需求 kW
LAMBDA = 1e5                      # 需求不满足的惩罚系数（远大于功率尺度）

def power(u):
    """总功率 kW"""
    return np.sum(A * u**2 + B * u + C)

def objective(u):
    """目标函数：功率 + 需求缺口惩罚（PSO 只需能算分数，不要求可导）"""
    gap = max(0.0, DEMAND - np.sum(CAP * u))
    return power(u) + LAMBDA * gap

def pso(obj, n_dim, lb, ub, np_=30, iters=120, w0=0.9, w1=0.4, c1=1.7, c2=1.7, seed=0):
    rng = np.random.default_rng(seed)
    lb, ub = np.full(n_dim, lb), np.full(n_dim, ub)
    x = rng.uniform(lb, ub, (np_, n_dim))
    v = rng.uniform(-(ub - lb) * 0.1, (ub - lb) * 0.1, (np_, n_dim))
    pbest_x, pbest_f = x.copy(), np.array([obj(p) for p in x])
    gi = int(np.argmin(pbest_f)); gbest_x, gbest_f = pbest_x[gi].copy(), pbest_f[gi]
    for k in range(iters):
        w = w0 - (w0 - w1) * k / iters                      # 惯性权重线性递减
        r1, r2 = rng.random((np_, n_dim)), rng.random((np_, n_dim))
        v = w * v + c1 * r1 * (pbest_x - x) + c2 * r2 * (gbest_x - x)
        x = np.clip(x + v, lb, ub)
        v = np.where((x <= lb) | (x >= ub), -0.5 * v, v)    # 触边界：截断 + 速度反向
        f = np.array([obj(p) for p in x])
        improved = f < pbest_f
        pbest_x[improved], pbest_f[improved] = x[improved], f[improved]
        gi = int(np.argmin(pbest_f))
        if pbest_f[gi] < gbest_f:
            gbest_x, gbest_f = pbest_x[gi].copy(), pbest_f[gi]
    return gbest_x, gbest_f

u_opt, j_opt = pso(objective, n_dim=2, lb=U_MIN, ub=U_MAX)
print("最优负荷率:", np.round(u_opt, 3), " 总功率 kW:", round(power(u_opt), 1),
      " 供冷 kW:", round(float(np.sum(CAP * u_opt)), 1))

# 对比：平均分配的能耗（工程上最常见的现状做法）
u_avg = np.full(2, DEMAND / CAP.sum())
save = 100 * (power(u_avg) - power(u_opt)) / power(u_avg)
print("平均分配:", np.round(u_avg, 3), " 总功率 kW:", round(power(u_avg), 1),
      " → 优化后节能率: %.2f%%" % save)

实际运行结果（本仓库 numpy 环境实测）：

最优负荷率: [0.545 0.811]  总功率 kW: 395.4  供冷 kW: 800.0
平均分配: [0.696 0.696]  总功率 kW: 400.2  → 优化后节能率: 1.20%

13.3.3 路线 B：强化学习（RL）——把序贯决策学成一张策略网

三要素设计（以冷冻站 / TEC 类对象为例）

其中：c_elec 为分时电价（元/kWh）；EF 为电力排放因子（见 §17.3）；λ 为碳价（元/kg CO₂），可先取 0 或按碳市场价做敏感性分析；越界惩罚必须「一票否决」级（CV 触碰工艺红线即给大额负奖励并终止回合）。

算法选型

训练环境 = 辨识模型沙盒

class SandboxEnv:
    """用 Hammerstein (FOPDT) 模型做 RL 训练环境（真机零风险）"""
    def __init__(self, model, dv_traj, dt=60.0, sp0=25.0, horizon=288):
        """
        model   : Hammerstein(FOPDT) 模型，predict(sp, dv) -> 下一时刻 CV
        dv_traj : 历史 DV 轨迹（环境温度/负荷），训练时回放真实工况波动
        horizon : 单回合步数（如 288 步 × 1 min）
        """
        self.model, self.dv, self.dt = model, dv_traj, dt
        self.sp0, self.horizon = sp0, horizon

    def reset(self):
        self.t, self.sp, self.cv = 0, self.sp0, self.model.steady_state(self.sp0)
        return self._obs(self.cv)

    def _obs(self, cv):
        return np.array([cv, self.sp, cv - self.sp, self.dv[self.t % len(self.dv)]])

    def step(self, sp_delta):
        self.sp = self.sp + sp_delta
        cv_next = self.model.predict(self.sp, dv=self.dv[self.t % len(self.dv)])
        r = self._reward(cv_next)
        self.cv, self.t = cv_next, self.t + 1
        done = self.t >= self.horizon or bool(self._violate(cv_next))
        return self._obs(cv_next), r, done

    def _reward(self, cv):
        elec = 0.68                       # 分时电价 元/kWh（示例）
        ef, carbon_price = 0.5366, 0.08   # 排放因子 kg/kWh、碳价 元/kg（见 §17.3）
        p = self.model.power(self.sp)     # 由静态非线性块给出功率 kW
        cost = (elec + carbon_price * ef) * p * self.dt / 3600.0
        jitter = 0.5 * (self.sp - self.sp0) ** 2
        return -(cost + 20.0 * self._violate(cv) + jitter)

    def _violate(self, cv):
        return 1.0 if (cv > self.cv_hi or cv < self.cv_lo) else 0.0   # 工艺红线

零侵入安全封装（上线强制）

def safe_action(sp_curr, a_raw, sp_range, max_step, predict_fn, cv_bounds, horizon=10):
    """
    动作安全层：任何 RL 输出都必须经此函数才允许下发（对应 §15 的安全闸门）
    ① 增量限幅 ② 速率限幅 ③ 模型前瞻屏蔽（预测越界则回退）
    """
    sp_new = np.clip(sp_curr + np.clip(a_raw, -max_step, max_step), *sp_range)
    # 前瞻校验：用模型预测未来 horizon 步，若触碰工艺边界则拒绝下发
    if predict_violates(predict_fn, sp_new, cv_bounds, horizon):
        return sp_curr, False          # 拒绝下发，返回 False 触发告警
    return sp_new, True

训练与部署流程（工程顺序不能颠倒）

1. 离线预训练：在模型沙盒中用历史 DV 轨迹回放训练，配合 Domain Randomization（辨识参数随机扰动 ±10~20%、加大噪声），逼策略对模型误差鲁棒；
2. 影子模式（Shadow Mode）：策略在线运行但只算不下发，与现行设定值并行对比 2~4 周，验证指标稳定优于现状；
3. 建议模式：推送给操作员确认后执行（建议先跑满一个完整工况周期，如冬夏各一轮）；
4. 自动模式：在 §15 四道安全闸门（旁路、速率限幅、工艺硬约束、看门狗）全部就位后才允许开启；在线更新频率受限（如每天一次），更新后先回影子模式复评；
5. 永不绕过：任何时候安全机制优先级都高于 RL 输出；通信中断 / 模型失配立即旁路回归原始设定。

13.3.4 两条路线怎么选？（含与 MPC 的关系）

13.3.5 算法之外的三道现实坎（见开篇 §0）

三种方案对比

14模型部署、软投切、在线监测、模型维护与迭代更新

阶段1 影子模式     阶段2 手动推荐      阶段3 半自动        阶段4 全自动监督
(只预测不下发) → (推荐给操作员) → (推荐+一键确认) → (自动下发+硬旁路)
  ≥1~2 周          ≥1 周              ≥2 周             长期，随时可切回

1. 阶段 1：影子模式（Shadow mode）——模型并行跑，只做预测，不输出设定值。对比模型预测值和真实 CV，持续计算预测误差，评估模型失配。建议运行 ≥ 1~2 周，覆盖足够工况变化；
2. 阶段 2：手动推荐模式——上层算出最优设定，展示给操作员，操作员确认后手动写入 PLC；系统不下发指令；
3. 阶段 3：半自动模式——自动生成推荐设定，下发前等待操作员一键确认；一旦超限，自动禁止下发；
4. 阶段 4：全自动监督下发——自动下发设定；硬旁路开关永久保留，一键切回原系统。

15工程落地风险与安全保障（零侵入重点章节）

1. 独立硬旁路：一键切回原始设定值，切换逻辑不依赖上层软件运行状态（物理开关或 PLC 内部独立回路实现，上层系统断电/死机也能切）；
2. 速率限幅：设定值单次最大变化量限制，防止设定值跳变冲击工艺。推荐值：单次变化 ≤ CV 量程的 2%~5%，且 ≤ 工艺工程师确认的单步安全步长；温度类对象典型值 0.5~1.0℃/次；
3. 工艺上下限硬约束：CV 触碰工艺边界（含预警裕度），立刻停止优化计算，切回原设定；禁止继续下发；
4. 通信看门狗：上层网关和 PLC 通信中断，看门狗超时自动退出优化，恢复原始设定；告警推送运维。推荐值：超时取设定值刷新周期的 3 倍（如刷新周期 60s，则看门狗 180s）；PLC 侧在超时后自动保持或回落到最后一个「原生安全设定值」，而非保持上层最后下发值不变；
5. 激励实验审批：激励实验前工艺、设备、安全三方确认，现场专人值守；激励程序自带紧急停止逻辑（一键急停 + 超限自停）；
6. 权限隔离：上层网关对 PLC 默认只读；写设定权限单独管控；所有写操作留存审计日志（时间、原值、新值、操作人/模块）；禁止修改 PLC 联锁、报警逻辑。

看门狗参考逻辑（伪代码）：

def watchdog_loop():
    while True:
        if now() - last_plc_heartbeat() > 3 * SETPOINT_PERIOD:
            # 通信中断：停发设定、通知 PLC 回落、告警
            plc.fallback_to_native_setpoint()
            alarm("上层-PLC 通信中断，已切回原始设定")
            require_manual_rearm()          # 恢复后须人工确认才能重新投运
        sleep(1)

16模型验收指标清单

验收在影子模式采集独立测试数据，全部指标达标，才允许进入半自动投运：

备注：指标为推荐工程阈值；可根据工艺要求协商调整。若部分指标不达标，返回激励/数据采集阶段重新开展实验。

17节能量与碳减排核算（M&V）

节能改造没有可审计的收益证明，就无法验收、无法结算、也无法说服管理层。本章给出核算框架（参考 IPMVP《国际节能效果测量与验证规程》的「选项 C：整体设施/反事实基线」思路）。

17.1 核心思想：反事实基线

• E_actual：投运后实测能耗（电表直读）；
• E_baseline：基线模型推算的「若不投运，在当前工况（产量、气温、来料）下应有的能耗」。

17.2 基线模型建立

1. 投运前基线窗口：采集 ≥ 2~4 周的原始运行数据（能耗 + 工况变量），覆盖典型工况波动；
2. 回归建模：建立「工况 → 能耗」的基线回归模型 E_baseline = g(产量, 气温, 来料温度, …)，可以用简单的多元线性回归/分箱平均，不必复杂——基线模型必须在投运前冻结存档，防止事后「调基线凑节能率」；
3. 非例行调整：产线改造、产品切换、设备大修后，基线需要重新标定并留档。

17.3 碳减排换算

• 电力排放因子以生态环境部最新发布的全国电力平均二氧化碳排放因子为准（2021 年度为 0.5568 kg CO₂/kWh，2022 年度为 0.5366 kg CO₂/kWh；区域电网因子另见当年发布值）；
• 天然气、蒸汽等其他能源按对应排放因子折算后加总。

参考实现（numpy）：基线回归 + 反事实推算 + 碳减排换算

import numpy as np

def fit_baseline(X, y):
    """基线回归：y ≈ X·θ（X 含常数列）。y=能耗 kWh，X 列=产量/气温/来料温度等工况。"""
    X = np.column_stack([np.ones(len(X)), X])           # 加截距
    theta, *_ = np.linalg.lstsq(X, y, rcond=None)
    resid = y - X @ theta
    r2 = 1 - resid @ resid / ((y - y.mean()) @ (y - y.mean()))
    return theta, r2                                     # ⚠️ 要求 R² ≥ 0.7，否则换变量/分箱

def mv_report(X_now, E_actual, theta, ef=0.5366, price=0.68, carbon_price=0.08):
    """M&V 核算：反事实基线 vs 实测能耗 → 节能量/节能率/碳减排/经济效益"""
    X_now = np.column_stack([np.ones(len(X_now)), X_now])
    E_base = X_now @ theta                               # 若不投运应有的能耗 kWh
    save_kwh = float(np.sum(E_base - E_actual))          # 节能量 kWh
    rate = 100.0 * save_kwh / float(np.sum(E_base))      # 节能率 %
    co2 = save_kwh * ef                                  # 碳减排 kg CO₂
    money = save_kwh * price + co2 * carbon_price        # 电费节省 + 碳资产收益 元
    return dict(E_baseline=float(np.sum(E_base)), E_actual=float(np.sum(E_actual)),
                save_kwh=save_kwh, save_rate=rate, co2_kg=co2, benefit_yuan=money)

17.4 核算报告要求

• 报告周期：月度/季度；
• 必须同时披露：基线模型及其冻结版本、工况数据完整性、天气/产量修正说明、异常时段剔除理由；
• 节能率口径必须注明是「单耗（单位产品能耗）」还是「总量」，避免产量变化掩盖或夸大节能效果。

18故障处置预案

1. 通信中断（上层网关 ↔ PLC）：看门狗超时，自动停止下发设定，PLC 侧回落到原生安全设定；产生告警，推送到运维。通信恢复后，默认停留在影子模式，需要人工确认再开启自动下发。
2. 模型漂移告警（预测误差超限）：维持当前运行，不立即切回；告警通知工程师；评估是否需要重新激励辨识；如果误差持续恶化，自动退出自动优化。
3. 被控变量 CV 触碰工艺边界：立刻停止优化计算，切回原设定；记录事件日志；禁止继续下发新设定，直到人工确认工艺恢复。
4. 传感器故障（测点卡死/跳变）：标记该测点无效；如果 DV/CV 为关键测点，直接退出自动优化，切回原始设定。
5. 激励实验紧急事件：程序内置紧急停止，一旦越限，立即终止激励序列，恢复实验前原始设定；保存本次采集数据（异常段数据对事后分析同样有价值）。

✅ 预案演练：投运前离线仿真演练全部故障场景，验证逻辑有效性。

19案例：Hammerstein (FOPDT) 模型小样验证 ——TEC 制冷系统（真实实验数据）

19.1 小样场景

固定风扇转速；通过调节 TEC 电流（0~5A）改变制冷量；采集测点：时间戳、设定电流、实测电流、实时功率、输出电压、冷却水温度、环境室温、冷风风口温度；被控输出：水杯水温。

• MV：TEC 实测电流（注意：用实测值，不是设定序列，见 §6.2）；
• DV：冷却水温度、室温、风口温度；
• CV：水杯水温；
• 额外观测：实时功率、输出电压。

🎯 建模目标：辨识 Hammerstein (FOPDT) 模型，找到满足水温设定值前提下最小 TEC 功率，作为算法验证载体。

激励方案：5 电平 PRBS（GF(5) m 序列，电平 0 / 1.25 / 2.5 / 3.75 / 5 A），采样周期 5s，每电平最少驻留 25s；开环激励（预先生成电流序列，水温不反馈修改输入）；输出保存 CSV 时序数据。

19.2 实验原始数据总览

读图要点：

• 上图：设定电流（蓝）按 5 电平 PRBS 切换，实测电流（橙）紧密跟随（恒流电源跟踪误差 < 1.5%，8s 内到位）；
• 下图：冷端风口温度（绿）对电流切换响应灵敏（电流大→风口冷，增益为负）；水杯水温（红）变化缓慢平滑——这就是典型的「大惯性 + 滞后」对象；
• 橙色阴影：初始 16.7 min 系统未达稳态（水杯从 41.9℃ 自然降温），这段数据必须剔除（见 19.6 教训一）。

19.3 数据质量诊断（实测）

19.4 静态非线性辨识（Hammerstein 前级）

• 左图：功率-电流关系，二次拟合 P = 1.05·I² + 4.25·I（W）。功率随电流超线性增长（5A 时 48 W，是 2.5A 时 16 W 的 3 倍）——这正是 TEC 内阻焦耳热的体现，也是「为什么不能只用线性 FOPDT」的直接证据；
• 右图：冷端降温幅度-电流关系。降温效果并非电流的线性放大（存在饱和趋势），静态非线性块必不可少。

把这两条曲线写成公式，就是 Hammerstein 模型的前级静态块（它同时给出「能省多少」与「花多少电」）：

其中 β 由右图的饱和程度回归得到（本例量级 β ≈ 0.12，见 §13.2.6 的 static_block）；K 是后级线性动态的静态增益，由 §19.5 辨识给出。静态块负责「非线性」，FOPDT 负责「动态」，两者串联即 §9 的 Hammerstein 结构。

参考实现（numpy）：静态块拟合 + 边际代价递增验证

import numpy as np

I = np.array([0.00, 1.25, 2.50, 3.75, 5.00])   # 实测电流 A（5 电平 PRBS）
P = np.array([ 4.6,  9.1, 16.2, 30.0, 48.1])   # 实测功率 W（§19.4 表）
P0 = P[0]                                      # 空载功率：风扇+控制板，与电流无关
A = np.column_stack([I ** 2, I])               # 过原点二次型：P − P0 = a·I² + b·I
(a, b), *_ = np.linalg.lstsq(A, P - P0, rcond=None)
fit = P0 + a * I ** 2 + b * I
r2 = 1 - np.sum((P - fit) ** 2) / np.sum((P - P.mean()) ** 2)
print("空载功率 %.2f W   二次项 a=%.3f  一次项 b=%.3f   R²=%.4f" % (P0, a, b, r2))
print("边际功率代价 dP/dI (W/A):", np.round(2 * a * I + b, 2))   # 每多要 1 A 多付多少瓦

运行结果：

空载功率 4.60 W   二次项 a=1.523  一次项 b=1.067   R²=0.9992
边际功率代价 dP/dI (W/A): [ 1.07  4.87  8.68  12.49  16.29]

19.5 动态辨识结果（FOPDT 线性子系统）

先说清「辨识到底在估什么」——估的是 Hammerstein 模型的后级线性动态（前级静态块已在 §19.4 标定）：电流 u 先过静态块得到隐变量 z = f(u)，再进一阶惯性 + 纯滞后环节。

工程上真正回归的只有三个未知数——y[k+1] = a·y[k] + b·u[k−nk] + c，一次最小二乘解出；再用上面三行反解出 τ / K / y_free。这就是下表那三个数字的全部来源（n_k 靠外层枚举 0~12 拍、取残差平方和最小者）。

把辨识结果代回去，两条通道可直接落地的模型就是下面这两行（T_s = 5 s）：

按 §10.1 双层优化（外层枚举滞后 0~12 拍，内层最小二乘）+ 训练/测试 8:2 时序切分，得到的辨识结果：

解读：

• 两条通道都是负增益（电流越大越冷），快通道 76 s（风道热容小）、慢通道 909 s（水杯热容大），与物理直觉一致；
• 增益的物理意义：每加 1 A 电流，稳态下风口再降约 1.1 ℃、水杯再降约 1.0 ℃；
• NRMSE 略高于 §16 的 5% 验收线，但注意分母效应：水杯全程只波动约 2 ℃，6.5% 对应的绝对误差仅 ±0.13 ℃，与温度传感器分辨率（0.1 ℃）同量级——模型已逼近测量极限。这也是 §11 提醒「NRMSE 要结合绝对误差解读」的实例。

参考实现（numpy，可脱离数据一键跑通）：PRBS 生成 → 滞后枚举辨识 → 一步/自由仿真验证

下面这份代码完整走一遍 §10.1 的「外层枚举滞后 + 内层最小二乘」，并把 §19.6 的两条教训做成可复算的对照实验。示例取快通道（电流 → 风口温度，τ = 76 s、θ = 20 s、K = −1.08 ℃/A）；把 synth() 换成读 CSV（u, y = np.loadtxt("prbs_run.csv", delimiter=",", unpack=True)）、再把通道参数换成慢通道，即复现 §19.5 表中的慢通道结果。

import numpy as np

TS = 5.0                       # 采样周期 s（§19.1：815 点 × 5 s）


def prbs5(n, c=(1, 2, 4), seed=(1, 0, 3)):
    """GF(5) m 序列：x[k] = (c0·x[k-1] + c1·x[k-2] + c2·x[k-3]) mod 5，周期 5³−1 = 124"""
    x = list(seed)
    while len(x) < n:
        x.append((c[0] * x[-1] + c[1] * x[-2] + c[2] * x[-3]) % 5)
    return np.array(x[:n])


def synth(n=815, tau=76.0, theta=20.0, K=-1.08, y_free=22.15, sigma=0.18, seed=7):
    """合成数据（无 CSV 时一键跑通）：5 电平 PRBS + 每电平驻留 25 s + FOPDT 快通道 + 噪声"""
    rng = np.random.default_rng(seed)
    u = np.repeat(1.25 * prbs5(n), 5)[:n]                  # 0 / 1.25 / 2.5 / 3.75 / 5 A
    nk = int(round(theta / TS))
    a, b = np.exp(-TS / tau), K * (1 - np.exp(-TS / tau))  # 离散 FOPDT 系数（§13.2.1）
    yc = np.zeros(n); yc[0] = y_free + K * u[:50].mean()   # 从稳态起步的「干净」轨迹
    for k in range(n - 1):
        u_eff = u[k - nk] if k >= nk else 0.0              # 纯滞后：取 nk 拍前的输入
        yc[k + 1] = a * yc[k] + (1 - a) * y_free + b * u_eff + rng.normal(0, sigma)
    # 开机未稳态：叠加一段「另一时间常数」(200 s) 的自然降温，起点 30.0 ℃
    y = yc + (30.0 - yc[0]) * np.exp(-np.arange(n) * TS / 200.0)
    return u, y + rng.normal(0, 0.05, n)                   # 叠加测量噪声


def fit_fopdt(u, y, nk_grid=range(0, 13)):
    """滞后枚举 + 内层最小二乘：y[k+1] = a·y[k] + b·u[k−nk] + c"""
    best = None
    for nk in nk_grid:
        Y = y[nk + 1:]
        X = np.column_stack([y[nk:-1], u[:len(Y)], np.ones(len(Y))])
        th, *_ = np.linalg.lstsq(X, Y, rcond=None)
        r = Y - X @ th
        if best is None or r @ r < best[0]:
            best = (float(r @ r), nk, th)
    sse, nk, (a, b, c) = best
    return dict(nk=nk, tau=-TS / np.log(a), K=b / (1 - a), steady=c / (1 - a), a=a, b=b, c=c)


def simulate(u, y, p, free=False, n0=0):
    """一步预测（喂真实 y） vs 自由仿真（模型输出回馈自身，对应部署形态）"""
    out = np.zeros(len(y)); out[n0] = y[n0]
    for k in range(n0, len(y) - 1):
        yk = out[k] if free else y[k]                      # ← 唯一的差别，结果差 3 倍
        u_eff = u[k - p["nk"]] if k >= p["nk"] else 0.0
        out[k + 1] = p["a"] * yk + p["b"] * u_eff + p["c"]
    return out


def score(y, yhat, n0):
    e = y[n0:] - yhat[n0:]
    rng = float(np.max(y[n0:]) - np.min(y[n0:]))
    rmse = float(np.sqrt(np.mean(e ** 2)))
    return rmse, 100 * rmse / rng, 1 - float(e @ e) / float(((y[n0:] - y[n0:].mean()) ** 2).sum())


u, y = synth()
WARM = int(16.7 * 60 / TS)                    # 待剔除的初始非稳态段（§19.6 教训一）
n_tr = WARM + int(0.8 * (len(y) - WARM))      # 训练 / 测试 8:2 时序切分

for tag, seg in (("含初始非稳态段", slice(0, n_tr)), ("剔除后（正确做法）", slice(WARM, n_tr))):
    p = fit_fopdt(u[seg], y[seg])
    o = score(y, simulate(u, y, p, n0=n_tr), n_tr)         # 一步预测（指标看不出差别！）
    print("%s → 滞后 %2d 拍  τ=%4.0f s  K=%.2f ℃/A | 一步预测 NRMSE %.1f%%  R² %.3f"
          % (tag, p["nk"], p["tau"], p["K"], o[1], o[2]))

p = fit_fopdt(u[WARM:n_tr], y[WARM:n_tr])
print("最终模型：滞后 %d 拍（%.0f s）  τ = %.0f s   K = %.2f ℃/A（真实值 76 s / −1.08 ℃/A）"
      % (p["nk"], p["nk"] * TS, p["tau"], p["K"]))
for tag, free in (("一步预测", False), ("自由仿真", True)):
    rmse, nrmse, r2 = score(y, simulate(u, y, p, free=free, n0=n_tr), n_tr)
    print("%s（测试段）：RMSE %.3f ℃   NRMSE %.1f%%   R² %.3f" % (tag, rmse, nrmse, r2))

运行结果（本仓库 numpy 环境实测）：

含初始非稳态段 → 滞后  4 拍  τ= 135 s  K=-1.78 ℃/A | 一步预测 NRMSE 4.9%  R² 0.955
剔除后（正确做法） → 滞后  4 拍  τ=  71 s  K=-0.92 ℃/A | 一步预测 NRMSE 4.9%  R² 0.955
最终模型：滞后 4 拍（20 s）  τ = 71 s   K = -0.92 ℃/A（真实值 76 s / −1.08 ℃/A）
一步预测（测试段）：RMSE 0.198 ℃   NRMSE 4.9%   R² 0.955
自由仿真（测试段）：RMSE 0.586 ℃   NRMSE 14.6%   R² 0.607

19.6 三个真实踩坑教训（本章最有价值的部分）

1. 教训一：初始非稳态段必须剔除。本次实验开跑时水杯 41.9℃（未预热/预冷到稳态），前 16.7 min 的自然降温与 PRBS 响应混在一起。对比试验：含初始段辨识 → 测试段 R² = 0.913；剔除后 → R² = 0.950。→ 对应 §7 第 2 条；采集脚本虽有 warmup 参数，但本次未用足，代价是 1/4 数据作废。
2. 教训二：一步预测好看 ≠ 自由仿真可用。同一组参数：一步预测 NRMSE 6.5%，自由仿真（模型输出回馈自身，不开任何「答案」）NRMSE 13.3%，直接翻倍。→ 一步拟合准则会让误差被每拍的真实测量「掩盖」，部署形态是长时域开环预测，正式模型必须用仿真误差准则（NLS）精修（§10）；线性最小二乘只适合扫滞后/阶次与产初值。
3. 教训三：无线传感器的随机测量滞后必须还原时间轴。四路温度走 LoRa（15s 更新），平台返回的是「最近一次上行值」，真实测量时刻比名义采样时刻早 0~15s 且随机抖动。直接当本拍值用，等于给输出通道加了随机纯滞后，辨识出的滞后参数必然有偏。→ 采集时已记录每拍 age，按「采样时刻 − age」还原真实时刻（§7）；age > 45s 的陈旧读数按缺测丢弃（绝不把几分钟前的旧值当新值写库）。

时间轴还原的参考实现（教训三，插在辨识之前）：

# t_arrive: 平台收到时刻；age_s: 读数自带的上行时延（LoRa 15 s 更新下 0~15 s 随机抖动）
t_true = t_arrive - age_s                          # 还原真实测量时刻
valid  = age_s <= 45.0                             # 陈旧读数按缺测丢弃，绝不把旧值当新值写库
t_uni  = np.arange(t_true[valid][0], t_true[valid][-1], TS)
y_uni  = np.interp(t_uni, t_true[valid], y[valid]) # 重采样到等间隔时间轴后再送进 fit_fopdt()

19.7 模型验证图

• 上图：实测水杯水温（红）vs 模型一步预测（蓝），测试段（虚线右侧）贴合良好；
• 中图：预测残差，零均值、无趋势，幅度 ±0.2 ℃ 以内（与传感器分辨率同量级）；
• 下图：残差 ACF，滞后 1~20 阶全部落在 95% 置信界内 → 残差接近白噪声，一阶结构对本对象充分，无需升阶（§9.3 阶次扫描原则的实例）。

19.8 案例小结

小样完整走通了「采集 → 预处理 → PE 检验 → 辨识 → 残差诊断 → 插值验证」全流程，得到两个可直接用于 MPC 能耗优化的模型（快/慢通道），功率静态模型 P = 1.05I² + 4.25I 即为 §13.1 能耗目标函数中 P(u) 的真实形态。算法验证完成后，迁移到工业流水线零侵入节能改造时，必须同时配套全部工业安全机制（§14~§18）。

20文档总结与参考文献

零侵入节能减碳改造，核心不是硬件改造，而是数据 + 辨识 + 上层监督优化。推荐模型结构：Hammerstein (FOPDT)，静态非线性块描述 TEC/执行器非线性特性，FOPDT 捕捉热惯性与传输滞后动态。整套工作流是迭代闭环：调研 → 激励 → 采集 → PE 准入 → 辨识 → 模型验证诊断 → 离线仿真 → 影子模式 → 分阶段软投切 → 在线监测迭代更新 → 节能量核算（M&V）闭环收益。

参考文献

1. Ljung L. System Identification: Theory for the User, 2nd ed. Prentice Hall, 1999.（系统辨识经典教材：PE 条件、预测误差方法、模型验证）
2. Nelles O. Nonlinear System Identification. Springer, 2001.（Hammerstein/Wiener 结构辨识、增益归一化、局部线性化）
3. Ogata K. Modern Control Engineering, 5th ed. Prentice Hall, 2010.（FOPDT、PID 整定基础）
4. Qin S J, Badgwell T A. "A survey of industrial model predictive control technology." Control Engineering Practice, 2003.（工业 MPC 综述，SSTO 分层结构）
5. EVO (Efficiency Valuation Organization). International Performance Measurement and Verification Protocol (IPMVP).（节能量测量与验证国际规程，第 17 章框架来源）
6. 生态环境部、国家统计局：《关于发布 2022 年电力二氧化碳排放因子的公告》，2024.（电力排放因子官方来源）