§4 RLS 在线辨识

递推最小二乘 · 遗忘因子 · 自校正 GPC 流程 · 符号约定 · 参数漂移修复

严格对照教材 §5.1 · 遗忘因子 λ=0.995，有效窗口 200 步 · 代码 gpc_na1.py

1. 为什么需要在线辨识

○ 背景

离线辨识的 ARX 参数是在特定工作点下得到的。实际水箱实验中：

水位变化范围大（1V → 3V），系统增益可能随工作点变化
泵的特性可能有漂移
离线辨识数据与实验当天的系统状态可能有偏差

RLS 在线辨识让控制器能自适应地更新模型参数，这就是"自校正 GPC"（Self-Tuning GPC）。

2. RLS 算法（教材§5.1）

★ 核心

ARX 模型（$n_a=1, n_b=0$）的仿真方程：$y(t) = -a_1 y(t-1) + b_0 u(t-k)$

RLS 辨识的参数向量和回归向量：

$$ \theta = [\theta_0,\ \theta_1]^T, \quad \varphi(t) = [y(t-1),\ u(t-k)]^T \tag{4.1}$$

符号约定（重要）

RLS 辨识的模型是 $y(t) = \theta_0 y(t-1) + \theta_1 u(t-k)$，所以 $\theta_0 = -a_1 > 0$（正值，如 0.9846）。
而 $A = [1, a_1]$ 中 $a_1 < 0$（如 $-0.9846$）。
代码中转换：A_cur = [1.0, -rls.theta[0]]，即 $a_1 = -\theta_0$。

RLS 更新公式（带遗忘因子 $\lambda$）：

$$ \mathbf{K}(t) = \frac{\mathbf{P}(t-1)\varphi(t)}{\lambda + \varphi^T(t)\mathbf{P}(t-1)\varphi(t)} \tag{4.2}$$ $$ \hat\theta(t) = \hat\theta(t-1) + \mathbf{K}(t)\left[y(t) - \varphi^T(t)\hat\theta(t-1)\right] \tag{4.3}$$ $$ \mathbf{P}(t) = \frac{1}{\lambda}\left[\mathbf{P}(t-1) - \mathbf{K}(t)\varphi^T(t)\mathbf{P}(t-1)\right] \tag{4.4}$$

初始值

$\hat\theta(0) = [0.984555,\ 0.024361]$（离线辨识结果），$\mathbf{P}(0) = 10^4 \mathbf{I}$（大初始不确定性，让 RLS 快速收敛到真实参数）

3. RLS 推导：从批量最小二乘到递推

★ 核心推导

第一步：批量最小二乘（LS）

已知 $t$ 时刻前的所有数据，最小化加权残差平方和：

$$ J_t(\theta) = \sum_{i=1}^{t} \lambda^{t-i}\left[y(i) - \varphi^T(i)\theta\right]^2 \tag{4.5}$$

令 $\partial J_t / \partial\theta = 0$，得正规方程：

$$ \underbrace{\left(\sum_{i=1}^{t} \lambda^{t-i}\varphi(i)\varphi^T(i)\right)}_{\mathbf{R}(t)}\,\hat\theta(t) = \sum_{i=1}^{t} \lambda^{t-i}\varphi(i)y(i) \tag{4.6}$$

定义 $\mathbf{P}(t) = \mathbf{R}^{-1}(t)$，则 $\hat\theta(t) = \mathbf{P}(t)\sum_{i=1}^{t}\lambda^{t-i}\varphi(i)y(i)$。

第二步：递推关系

注意到 $\mathbf{R}(t) = \lambda\mathbf{R}(t-1) + \varphi(t)\varphi^T(t)$，利用矩阵求逆引理（Sherman-Morrison 公式）：

$$ (A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^TA^{-1}}{1 + v^TA^{-1}u} \tag{4.7}$$

令 $A = \lambda\mathbf{R}(t-1)$，$u = \varphi(t)$，$v = \varphi(t)$，得：

$$ \mathbf{P}(t) = \frac{1}{\lambda}\left[\mathbf{P}(t-1) - \frac{\mathbf{P}(t-1)\varphi(t)\varphi^T(t)\mathbf{P}(t-1)}{\lambda + \varphi^T(t)\mathbf{P}(t-1)\varphi(t)}\right] \tag{4.8}$$

定义增益向量 $\mathbf{K}(t) = \mathbf{P}(t-1)\varphi(t) / (\lambda + \varphi^T(t)\mathbf{P}(t-1)\varphi(t))$，则式(4.8) 化简为教材式(4.4)。

第三步：参数更新

将 $\hat\theta(t) = \mathbf{P}(t)\sum_{i=1}^{t}\lambda^{t-i}\varphi(i)y(i)$ 展开，利用 $\mathbf{P}(t)$ 的递推关系，得：

$$ \hat\theta(t) = \hat\theta(t-1) + \mathbf{K}(t)\underbrace{\left[y(t) - \varphi^T(t)\hat\theta(t-1)\right]}_{\text{新息（预测误差）}} \tag{4.9}$$

这就是教材式(4.3)。每步只需做一次矩阵-向量乘法，计算量 $O(n^2)$，不需要重新解整个最小二乘。

$\mathbf{P}(t)$ 的物理含义

$\mathbf{P}(t) = \mathbf{R}^{-1}(t)$ 是参数估计的协方差矩阵（乘以噪声方差后）。$\mathbf{P}$ 越大，说明参数估计越不确定，新数据的权重越大（增益 $\mathbf{K}$ 越大）。初始 $\mathbf{P}(0) = 10^4\mathbf{I}$ 表示对初始参数完全不确定，让 RLS 快速收敛。

4. 遗忘因子的推导与适用范围

★ 核心推导

遗忘因子的来源

式(4.5) 中 $\lambda^{t-i}$ 对历史数据指数加权衰减。等效有效数据窗口长度：

$$ N_\text{eff} = \sum_{i=0}^{\infty} \lambda^i = \frac{1}{1-\lambda} \tag{4.10}$$

$\lambda=0.995$ 时 $N_\text{eff}=200$，即 RLS 相当于用最近 200 步数据做加权最小二乘。

稳态时的参数漂移问题

当系统处于稳态（$\Delta y \approx 0$，$\Delta u \approx 0$）时，回归向量 $\varphi(t) \approx \text{const}$，信息矩阵 $\mathbf{R}(t)$ 不再增长，但 $\mathbf{P}(t) = \mathbf{R}^{-1}(t)$ 会因遗忘因子而持续增大：

$$ \mathbf{P}(t) \approx \frac{1}{\lambda}\mathbf{P}(t-1) \quad \Rightarrow \quad \mathbf{P}(t) \to \infty \tag{4.11}$$

$\mathbf{P}$ 发散后，增益 $\mathbf{K}$ 变大，RLS 对噪声过度敏感，参数开始漂移。这就是 $\lambda=0.98$ 时 $b_0$ 从 0.024 漂到 0.011 的根本原因。

适用范围

条件	适用	说明
系统参数缓慢时变	✅	$\lambda$ 越小跟踪越快，但稳态漂移越严重
系统参数固定	⚠️	$\lambda=1$ 最优，但 $\lambda<1$ 会引入稳态漂移
持续激励（PE）条件满足	✅	$\varphi(t)$ 足够丰富，$\mathbf{R}(t)$ 不退化
稳态运行（无激励）	❌	$\mathbf{P}$ 发散，参数漂移，需要 $\lambda \geq 0.99$
参数突变	⚠️	需要 $\lambda$ 足够小（如 0.95），但会牺牲稳态稳定性

水箱系统的选择依据

水箱是慢时变系统（时间常数 $\tau \approx 64$ 步），参数变化缓慢。实验中系统会在多个工作点运行（1V→2V→3V），需要 RLS 能跟踪工作点变化，但不能对稳态噪声过度敏感。
$\lambda=0.995$（$N_\text{eff}=200$ 步）是合理折中：足够慢以抑制稳态漂移，足够快以跟踪工作点切换。

★ 核心

遗忘因子 $\lambda \in (0,1]$ 决定历史数据的权重衰减速度：

$$ N_\text{eff} = \frac{1}{1-\lambda} \quad \text{（有效数据窗口长度）} \tag{4.4}$$

$\lambda$	$N_\text{eff}$	适用场景
0.98	50 步	快速时变系统，但稳态时参数漂移严重 ❌
0.995	200 步	水箱系统（慢动态），参数稳定 ✅
0.999	1000 步	几乎不更新，接近离线辨识

第一次实验的教训

第一版用 $\lambda=0.98$，有效窗口只有 50 步。稳态时 RLS 持续用噪声数据更新，$b_0$ 从 0.024 漂到 0.011（减半），控制器认为系统增益变小，输出更大的控制量，导致超调加剧。改为 $\lambda=0.995$ 后参数稳定。

4. 自校正 GPC 流程

○ 每步执行顺序

读取传感器 $y(t)$，归一化得 $y_\text{dev} = (y-y_0)/Y_R$
从 RLS 取当前参数 $\hat\theta = [\theta_0, \theta_1]$，构造 $A = [1, -\theta_0]$
调用 diophantine(A, b0, p)，得 $F\_\text{list}[0..p-1]$、$G\_\text{list}[0..p-1]$
调用 free_response(F_list, G_list, y_buf, du_buf, p, k)，得 $Y^1 \in \mathbb{R}^p$
调用 soften_setpoint(y_buf[-1], SP, alpha, p)，得柔化轨迹 $W \in \mathbb{R}^p$
调用 stepped_gpc(G_list, Y1, W, p, pu, lam, beta)，得 $\delta$，限幅后输出 $u(t)$
更新 du_buf（追加 $\Delta u_t$），更新 y_buf、u_buf
RLS 更新：$\varphi = [y(t-1),\ u(t-k)]$，调用 rls.update(phi, y_t)

注意顺序

RLS 更新在控制律计算之后，用的是当前步的 $y(t)$ 作为观测值。这样 RLS 看到的是"施加了控制量 $u(t-k)$ 之后的真实输出"，回归向量和观测值时间对齐。

5. C++ 实现细节

○ 工程细节

RLS 的 $2\times2$ 协方差矩阵 $\mathbf{P}$ 展开为 4 个标量，避免动态内存：

static double rls_P00=1e4, rls_P01=0.0;
static double rls_P10=0.0, rls_P11=1e4;

// 更新：Pp = P*phi
double Pp0 = rls_P00*phi0 + rls_P01*phi1;
double Pp1 = rls_P10*phi0 + rls_P11*phi1;
double denom = RLS_LAM + phi0*Pp0 + phi1*Pp1;
double K0 = Pp0/denom,  K1 = Pp1/denom;
rls_a1 += K0*e;  rls_b0 += K1*e;
rls_P00 = (rls_P00 - K0*Pp0) / RLS_LAM;
// ... 其余3项类似

MSVC 旧版本 C2374 问题：所有 for (int i/j/k = ...) 循环用 { } 包裹，创建独立作用域。

6. 踩坑：参数漂移

⚠ 实际 bug（已修复）

现象：仿真中 $b_0$ 从初始值 0.024 持续下降到 0.011，控制量超调越来越大。

根本原因：$\lambda=0.98$ 有效窗口只有 50 步，稳态时系统几乎不动（$\Delta y \approx 0$），RLS 用噪声数据持续更新，参数向零漂移。

修复：$\lambda=0.995$，有效窗口 200 步，稳态时参数更新量极小。

	$\lambda=0.98$	$\lambda=0.995$
最终 $b_0$	0.011（漂移 55%）	0.0249（稳定）
RMSE	0.346 V	0.305 V