§6 补充答疑
纯滞后：泵的驱动电流变化后，水流经管道到达水箱需要一段固定时间，这段时间内输出完全不响应。这就是纯时延 $k=24$ 步（24秒）。
惯性：水箱有容积，水位不能瞬间变化——就像电容充电，需要时间积累。这产生了"惯性"（低通滤波效应）。
系统的"阶次"对应传递函数分母多项式的次数，等于系统的独立储能元件数量：
| 阶次 | 物理含义 | 传递函数 | 阶跃响应 |
|---|---|---|---|
| 一阶 | 1个储能元件（水箱容积） | $G(s) = \frac{K}{\tau s + 1} e^{-\theta s}$ | 指数上升，无超调 |
| 二阶 | 2个储能元件（如串联两个水箱，或水箱+管道惯性） | $G(s) = \frac{K}{(\tau_1 s+1)(\tau_2 s+1)} e^{-\theta s}$ | S形上升，可能有超调 |
为什么 AIC 选出 na=4 而不是 na=1？
真实水箱并不是理想的一阶系统。管道、阀门、传感器滤波都会引入额外的动态，使得实际响应比纯一阶更复杂。AIC 选 na=4 是因为数据里确实有这些高阶动态。但对于控制器设计，na=1 已经捕捉了主要动态（主极点 0.9846），RLS 在线更新可以弥补模型误差。
本实验 $a_1 = 0.9846$（极点），$T=1$s，则 $\tau = -T/\ln(0.9846) \approx 64$ 步 = 64s。
这与水箱的物理时间常数（约 1~2 分钟）吻合。
AIC（Akaike Information Criterion，赤池信息准则）的核心思想是在拟合精度和模型复杂度之间取平衡：
$$ \text{AIC} = N \ln(\hat\sigma_e^2) + 2 \cdot n_\text{params} \tag{6.1}$$- $N \ln(\hat\sigma_e^2)$：残差方差的对数，越小说明拟合越好（精度项）
- $2 n_\text{params}$：参数数量的惩罚（复杂度项），防止过拟合
- $n_\text{params} = n_a + n_b + 1$（ARX 模型的参数个数）
AIC 来自信息论，本质是最大化模型的对数似然函数。对于高斯噪声假设下的 ARX 模型，对数似然正比于 $-N\ln(\hat\sigma_e^2)$，所以最小化 AIC 等价于在精度和复杂度之间找最优平衡。
直接用 RMSE 比较不同阶次的模型会偏向高阶（参数越多 RMSE 越小），AIC 的惩罚项纠正了这个偏差。
量化依据：AIC 差值 $\Delta\text{AIC} = \text{AIC}_i - \text{AIC}_\text{min}$
- $\Delta\text{AIC} < 2$：两个模型几乎等价，选简单的
- $\Delta\text{AIC} \in [2, 10]$：较弱的证据支持更复杂的模型
- $\Delta\text{AIC} > 10$：强证据支持更复杂的模型
本实验 na=1 vs na=4 的 $\Delta\text{AIC} = -21843 - (-21452) = -391$，差距极大，说明 na=4 在统计上显著优于 na=1。但我们选 na=1 是出于工程考虑（简单、便于调试），不是统计最优。
问题圈定：已知 $N$ 组输入输出数据 $\{u(t), y(t)\}$，假设系统满足 ARX 模型：
$$ y(t) = \varphi(t)^T \theta + e(t) \tag{6.2}$$其中 $\varphi(t) = [-y(t-1), u(t-k)]^T$（na=1, nb=0），$\theta = [a_1, b_0]^T$ 是待估参数，$e(t)$ 是白噪声。
矩阵形式：把所有时刻堆叠：
$$ \mathbf{Y} = \mathbf{\Phi}\,\theta + \mathbf{e}, \quad \mathbf{Y} \in \mathbb{R}^M,\; \mathbf{\Phi} \in \mathbb{R}^{M \times 2} \tag{6.3}$$最小二乘目标：最小化残差平方和：
$$ \min_\theta \|\mathbf{Y} - \mathbf{\Phi}\theta\|^2 \tag{6.4}$$求导令其为零（正规方程）：
$$ \frac{\partial}{\partial\theta}\|\mathbf{Y} - \mathbf{\Phi}\theta\|^2 = -2\mathbf{\Phi}^T(\mathbf{Y} - \mathbf{\Phi}\theta) = 0 \tag{6.5}$$ $$ \Rightarrow \mathbf{\Phi}^T\mathbf{\Phi}\,\hat\theta = \mathbf{\Phi}^T\mathbf{Y} \tag{6.6}$$ $$ \Rightarrow \hat\theta = (\mathbf{\Phi}^T\mathbf{\Phi})^{-1}\mathbf{\Phi}^T\mathbf{Y} \tag{6.7}$$设原始空间的模型为 $y = a_1 y_{-1} + b_0^\text{raw} u_{-k}$，归一化后 $y_n = y/Y_R$，$u_n = u/U_R$（简化，忽略偏移），则：
$$ y_n \cdot Y_R = a_1 \cdot y_{n,-1} \cdot Y_R + b_0^\text{raw} \cdot u_{n,-k} \cdot U_R$$ $$ y_n = a_1 \cdot y_{n,-1} + \underbrace{b_0^\text{raw} \cdot \frac{U_R}{Y_R}}_{b_0^\text{norm}} \cdot u_{n,-k} \tag{6.8}$$所以 $b_0^\text{norm} = b_0^\text{raw} \times U_R/Y_R = 0.0055 \times 16/3.634 \approx 0.0242$，与直接在归一化空间辨识的结果 0.02436 吻合。
$a_1$ 不变，因为它是 $y/y$ 的比值，量纲相消。
教材（第五章 §5.2）的 CARIMA 模型式(5.1)：
$$ A(q^{-1})\Delta y(t) = B(q^{-1})\Delta u(t-1) + C(q^{-1})\xi(t) \tag{教材5.1}$$丢番图恒等式式(5.3)：$1 = E_j(q^{-1})\tilde A(q^{-1}) + q^{-j}F_j(q^{-1})$，其中 $\tilde A = A\Delta$。
教材 §5.3 式(5.9)~(5.11) 给出的 $F_j$ 初值（$j=1$）：
$$ F_1(q^{-1}) = (1+a_1) + (a_1+a_2)q^{-1} + \cdots + (a_{n_a-1}+a_{n_a})q^{-(n_a-1)} + a_{n_a}q^{-n_a} \tag{教材5.10}$$这里 $F_1$ 有 $n_a+1$ 项（$q^0$ 到 $q^{-n_a}$），与我们代码一致。
教材式(5.10) 的系数写法是 $(a_{i-1}+a_i)$，而我们代码用的是 $(\tilde A_i - \tilde A_{i+1})$。两者等价，只是展开方式不同：
$$ \tilde A_i - \tilde A_{i+1} = (a_{i-1} - a_i) - (a_i - a_{i+1}) \quad \text{（对 } \tilde A = A\Delta \text{ 展开后）}$$教材式(5.18) 的递推：$f_j^i = f_{j-1}^{i+1} - e_{j-1}(a_{i+1}-a_i)$，与我们代码完全一致。
结论：代码与教材完全一致，只是初值的书写形式不同。
$G_j$ 的长度：教材式(5.16) $G_j = E_j(q^{-1})B(q^{-1})$，长度 $= (j-1) + (n_b+1) = n_b+j$。本实验 $n_b=0$，所以 $G_j$ 长度 $= j$，与教材一致。
| 量 | 教材定义 | 本实验（na=1, nb=0） |
|---|---|---|
| $\deg E_j$ | $j-1$ | $j-1$（代码不显式存储） |
| $F_j$ 项数 | $n_a+1$（对 $\tilde A$ 递推） | $n_a+1 = 2$（但存3项，含 $q^{-2}$ 项） |
| $G_j$ 长度 | $n_b+j$ | $j$（$n_b=0$） |
教材 §5.2 式(5.2) 明确写出：$\tilde A(q^{-1}) = A(q^{-1})\Delta(q^{-1})$，阶次 $= n_a + 1$。
逐步拆解（na=1）：
$A = 1 + a_1 q^{-1}$（一阶，$a_1 = -0.9846$），$\Delta = 1 - q^{-1}$（一阶），则：
$$ \tilde A = A \cdot \Delta = (1 + a_1 q^{-1})(1 - q^{-1}) = 1 + (a_1-1)q^{-1} - a_1 q^{-2} \tag{6.9}$$代入数值：$\tilde A = [1,\ -1.9846,\ 0.9846]$，阶次 $= 2$，有 3 个系数。
恒等式 $1 = E_j \tilde A + q^{-j} F_j$ 中，两边次数平衡要求：
- $\deg(E_j \tilde A) = (j-1) + (n_a+1) = j + n_a$
- $\deg(q^{-j} F_j) = j + \deg F_j$
- 两边最高次相等：$\deg F_j = n_a$
$na=1$ 时 $\deg F_j = 1$，即 $F_j$ 有 $q^0$ 和 $q^{-1}$ 两项，系数向量长度 $= n_a + 1 = 2$。
但代码里存了 3 项（$f_j^0, f_j^1, f_j^2$），这是因为递推公式式(5.18) 中 $f_j^{n_a} = e_{j-1} \cdot a_{n_a}$ 用到了 $\tilde A$ 的第 $n_a+1$ 个系数（即 $\tilde A_2 = a_1 = 0.9846$），所以实现时把 $F_j$ 扩展为与 $\tilde A$ 等长（3项），第3项 $f_j^2$ 对应 $q^{-2}$ 系数。
自由响应为什么需要 $y(t-2)$：
$Y_1[j] = F_j(q^{-1}) y(t) = f_j^0 y(t) + f_j^1 y(t-1) + f_j^2 y(t-2)$。$\tilde A$ 是二阶，系统的"记忆"深度是 2 步，所以自由响应需要 $y(t), y(t-1), y(t-2)$ 三个值。
规律：$F_j$ 存储长度 $= n_a + 2$（$= \deg\tilde A + 1$），自由响应需要 $n_a+2$ 个历史 $y$ 值。
教材 §5.1 图1 给出了自控正控制的基本框图：辨识器（在线估计模型参数）+ 设计器（根据当前参数设计控制律）+ 控制器（执行控制量）。RLS 就是辨识器的具体实现。
RLS 的本质：批量最小二乘 $\hat\theta = (\Phi^T\Phi)^{-1}\Phi^T Y$ 的递推版本。利用矩阵求逆引理，每来一个新数据点只做增量更新，不重新解整个方程组：
$$ \mathbf{K}(t) = \frac{\mathbf{P}(t-1)\varphi(t)}{\lambda + \varphi^T(t)\mathbf{P}(t-1)\varphi(t)} \tag{6.10}$$ $$ \hat\theta(t) = \hat\theta(t-1) + \mathbf{K}(t)\left[y(t) - \varphi^T(t)\hat\theta(t-1)\right] \tag{6.11}$$ $$ \mathbf{P}(t) = \frac{1}{\lambda}\left[\mathbf{P}(t-1) - \mathbf{K}(t)\varphi^T(t)\mathbf{P}(t-1)\right] \tag{6.12}$$其中 $\mathbf{P}(t)$ 是参数估计的协方差矩阵，$\lambda$ 是遗忘因子。
带遗忘因子的 RLS 等价于对历史数据指数加权：
$$ \min_\theta \sum_{i=1}^{t} \lambda^{t-i} \left[y(i) - \varphi^T(i)\theta\right]^2 \tag{6.13}$$$\lambda < 1$ 时旧数据权重指数衰减，有效窗口 $N_\text{eff} = 1/(1-\lambda)$。$\lambda=0.995$ 对应 200 步有效窗口，适合水箱这种慢时变系统。
RLS 对 GPC 推导链的影响：
教材 §5.2~§5.4 的 GPC 推导假设 $A, B$ 已知固定。引入 RLS 后，每步 $\hat\theta$ 更新，整个推导链都要重算：
| 推导步骤 | 固定参数 GPC | 自校正 GPC（+RLS） |
|---|---|---|
| ①构造 $\tilde A = A\Delta$ | 离线一次 | 每步用最新 $\hat a_1$ 重算 |
| ②丢番图递推（§5.3） | 离线预计算 $F_j, G_j$ | 每步在线递推 |
| ③自由响应 $Y^1$（§5.4） | 用固定 $F_j$ 计算 | 用当步 $F_j$ 计算 |
| ④控制律 $\delta$（§5.5） | 用固定 $G_2$, denom | 用当步 $G_2$, denom |
| ⑤RLS 更新 | 无 | 用 $y(t)$ 更新 $\hat\theta$，供下步用 |
控制律用的是"当前对系统的最佳估计"，然后施加控制量 $u(t)$，等到下一步观测到 $y(t+1)$ 后再更新参数。如果先更新参数再计算控制律，会用到还没发生的观测，逻辑上不对。
本实验的顺序：读 $y(t)$ → 用旧 $\hat\theta$ 算控制律 → 输出 $u(t)$ → 用 $y(t)$ 更新 $\hat\theta$ → 等待下一步。