系统辨识是从输入输出数据中建立动态系统数学模型的科学与艺术。在实际控制工程中，精确的系统模型往往难以从第一性原理推导，或者系统参数会随时间变化。系统辨识技术使我们能够从实验数据中提取系统的动态特性，为控制器设计提供可靠的模型基础。本章将介绍主要的辨识方法，从经典的参数估计到现代的子空间方法，并特别关注实际应用中的挑战，如闭环辨识和非线性系统辨识。
系统辨识的核心任务是找到一个数学模型 $M$，使其能够最好地描述观测到的输入输出数据：
\[y(t) = G(q, \theta)u(t) + H(q, \theta)e(t)\]其中：
实验设计 → 数据采集 → 模型结构选择 → 参数估计 → 模型验证
↑ ↓
←────────────── 不满足要求则重新设计 ←─────────────
考虑ARX模型： \(A(q)y(t) = B(q)u(t) + e(t)\)
其中： \(A(q) = 1 + a_1q^{-1} + \cdots + a_{n_a}q^{-n_a}\) \(B(q) = b_1q^{-1} + \cdots + b_{n_b}q^{-n_b}\)
写成线性回归形式： \(y(t) = \phi^T(t)\theta + e(t)\)
其中回归向量： \(\phi(t) = [-y(t-1), \ldots, -y(t-n_a), u(t-1), \ldots, u(t-n_b)]^T\)
参数向量： \(\theta = [a_1, \ldots, a_{n_a}, b_1, \ldots, b_{n_b}]^T\)
最小二乘估计： \(\hat{\theta}_{LS} = \arg\min_\theta \sum_{t=1}^N [y(t) - \phi^T(t)\theta]^2\)
解析解： \(\hat{\theta}_{LS} = \left[\sum_{t=1}^N \phi(t)\phi^T(t)\right]^{-1} \sum_{t=1}^N \phi(t)y(t)\)
实时更新参数估计：
\[\hat{\theta}(t) = \hat{\theta}(t-1) + K(t)[y(t) - \phi^T(t)\hat{\theta}(t-1)]\]其中增益矩阵： \(K(t) = P(t-1)\phi(t)[1 + \phi^T(t)P(t-1)\phi(t)]^{-1}\)
协方差矩阵更新： \(P(t) = [I - K(t)\phi^T(t)]P(t-1)\)
遗忘因子：处理时变系统 \(P(t) = \frac{1}{\lambda}[I - K(t)\phi^T(t)]P(t-1)\)
其中 $0.95 \leq \lambda < 1$ 是遗忘因子。
假设噪声服从高斯分布 $e(t) \sim \mathcal{N}(0, \sigma^2)$，似然函数：
\[L(\theta) = \prod_{t=1}^N \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{[y(t) - \hat{y}(t|\theta)]^2}{2\sigma^2}\right)\]对数似然： \(\log L(\theta) = -\frac{N}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{t=1}^N [y(t) - \hat{y}(t|\theta)]^2\)
MLE估计： \(\hat{\theta}_{MLE} = \arg\max_\theta \log L(\theta)\)
一般化的辨识准则：
\[\hat{\theta} = \arg\min_\theta V_N(\theta) = \arg\min_\theta \frac{1}{N}\sum_{t=1}^N \ell(\varepsilon(t, \theta))\]其中：
对于不同模型结构，预测误差有不同形式：
通过正弦激励获取频率响应：
输入：$u(t) = A\sin(\omega t)$ 稳态输出：$y(t) = |G(j\omega)|A\sin(\omega t + \angle G(j\omega))$
频率响应函数（FRF）： \(\hat{G}(j\omega) = \frac{Y(j\omega)}{U(j\omega)} = \frac{\mathcal{F}\{y(t)\}}{\mathcal{F}\{u(t)\}}\)
经验传递函数估计（ETFE）： \(\hat{G}_{ETFE}(e^{j\omega}) = \frac{Y_N(e^{j\omega})}{U_N(e^{j\omega})}\)
其中 $Y_N$ 和 $U_N$ 是DFT变换。
谱分析估计： \(\hat{G}(e^{j\omega}) = \frac{\hat{\Phi}_{yu}(e^{j\omega})}{\hat{\Phi}_u(e^{j\omega})}\)
其中：
给定频率响应数据 ${G_k, \omega_k}_{k=1}^M$，拟合参数模型：
\[\hat{\theta} = \arg\min_\theta \sum_{k=1}^M W_k |G_k - G(e^{j\omega_k}, \theta)|^2\]其中 $W_k$ 是频率相关的权重。
子空间方法直接从输入输出数据辨识状态空间模型：
\[\begin{aligned} x_{k+1} &= Ax_k + Bu_k + w_k \\ y_k &= Cx_k + Du_k + v_k \end{aligned}\]无需迭代优化，通过奇异值分解（SVD）获得模型。
构造Hankel矩阵：
\[Y_{p} = \begin{bmatrix} y_0 & y_1 & \cdots & y_{j-1} \\ y_1 & y_2 & \cdots & y_j \\ \vdots & \vdots & \ddots & \vdots \\ y_{i-1} & y_i & \cdots & y_{i+j-2} \end{bmatrix}\]类似构造 $U_p$（过去输入）、$Y_f$（未来输出）、$U_f$（未来输入）。
构造数据矩阵并进行QR分解： \(\begin{bmatrix} U_p \\ U_f \\ Y_p \\ Y_f \end{bmatrix} = \begin{bmatrix} R_{11} & 0 & 0 & 0 \\ R_{21} & R_{22} & 0 & 0 \\ R_{31} & R_{32} & R_{33} & 0 \\ R_{41} & R_{42} & R_{43} & R_{44} \end{bmatrix} \begin{bmatrix} Q_1^T \\ Q_2^T \\ Q_3^T \\ Q_4^T \end{bmatrix}\)
对 $R_{32}$ 进行SVD： \(R_{32} = U\Sigma V^T = \begin{bmatrix} U_1 & U_2 \end{bmatrix} \begin{bmatrix} \Sigma_1 & 0 \\ 0 & \Sigma_2 \end{bmatrix} \begin{bmatrix} V_1^T \\ V_2^T \end{bmatrix}\)
确定系统阶次 $n$（通过奇异值的跳变）
提取系统矩阵：
N4SID（Numerical algorithms for Subspace State Space System IDentification）是另一种流行的子空间方法：
斜投影： \(Y_f / U_f \perp W_p\)
其中 $W_p = \begin{bmatrix} U_p \ Y_p \end{bmatrix}$
通过SVD估计状态序列
使用最小二乘估计系统矩阵
级联结构：静态非线性 + 线性动态 + 静态非线性
u → [静态非线性 f(·)] → v → [线性系统 G(q)] → w → [静态非线性 h(·)] → y
辨识策略：
非线性自回归移动平均模型：
\[y(t) = f(y(t-1), \ldots, y(t-n_y), u(t-1), \ldots, u(t-n_u), e(t-1), \ldots, e(t-n_e)) + e(t)\]常用非线性函数形式：
循环神经网络（RNN）模型： \(\begin{aligned} h_t &= \sigma(W_h h_{t-1} + W_u u_t + b_h) \\ y_t &= W_y h_t + b_y \end{aligned}\)
LSTM模型：处理长期依赖 \(\begin{aligned} f_t &= \sigma(W_f [h_{t-1}, u_t] + b_f) \\ i_t &= \sigma(W_i [h_{t-1}, u_t] + b_i) \\ \tilde{C}_t &= \tanh(W_C [h_{t-1}, u_t] + b_C) \\ C_t &= f_t \odot C_{t-1} + i_t \odot \tilde{C}_t \\ o_t &= \sigma(W_o [h_{t-1}, u_t] + b_o) \\ h_t &= o_t \odot \tanh(C_t) \end{aligned}\)
系统动力学的稀疏辨识：
\[\dot{x} = \Theta(x, u)\xi\]其中：
通过L1正则化促进稀疏性： \(\hat{\xi} = \arg\min_\xi ||\dot{X} - \Theta\xi||_2^2 + \lambda||\xi||_1\)
在许多实际系统中，出于安全或稳定性考虑，必须在闭环条件下进行辨识：
r(t) →⊕→ [控制器 C] → u(t) → [系统 G] → y(t)
↑ ↓
←────────────────────────────────────
主要问题：
忽略反馈回路，直接使用开环辨识方法：
\[\hat{G} = \arg\min_G \sum_{t=1}^N [y(t) - G(q)u(t)]^2\]适用条件：
同时辨识闭环传递函数，然后推导开环模型：
闭环传递函数： \(T(q) = \frac{GC}{1 + GC}, \quad S(q) = \frac{1}{1 + GC}\)
从测量数据辨识 $T$ 和 $S$，然后计算： \(G = \frac{T}{C \cdot S}\)
使用输入输出的联合信息构造无偏估计：
\[\hat{\theta} = \arg\min_\theta \left|\left|\begin{bmatrix} Y \\ U \end{bmatrix} - \begin{bmatrix} G(\theta) \\ I \end{bmatrix} \begin{bmatrix} U \\ E \end{bmatrix}\right|\right|_F^2\]第一级：辨识闭环灵敏度函数 \(\hat{S} = \arg\min_S ||Y - S \cdot R||^2\)
第二级：使用模拟开环数据 \(\tilde{u} = \hat{S} \cdot r, \quad \tilde{y} = \hat{S} \cdot y\)
然后用开环方法辨识： \(\hat{G} = \arg\min_G ||\tilde{y} - G \cdot \tilde{u}||^2\)
考虑闭环下的预测误差：
\[\varepsilon(t, \theta) = H^{-1}(\theta)[y(t) - G(\theta)u(t)]\]其中 $u(t) = C[r(t) - y(t)]$
优化准则： \(\hat{\theta} = \arg\min_\theta \sum_{t=1}^N \varepsilon^2(t, \theta)\)
激励信号选择：
采样率选择：
数据长度：
残差分析：
交叉验证：
FIT值： \(\text{FIT} = 100 \times \left(1 - \frac{||y - \hat{y}||}{||y - \bar{y}||}\right) \%\)
信息准则：
其中 $p$ 是参数数量，$\hat{\sigma}^2$ 是残差方差。
锂离子电池的精确建模对电池管理系统（BMS）至关重要。等效电路模型（ECM）因其简单性和实用性被广泛采用。我们需要从充放电数据中辨识模型参数，用于SOC（State of Charge）估计和健康状态监测。
R0
I →──/\/\/──┬──────────────────→ V
│
┌┴┐ R1 R2
│ │OCV ║ ║
└┬┘ ┌┴┐ ┌┴┐
│ │C1│ │C2│
│ └┬┘ └┬┘
│ │ │
└──────┴─────────┘
模型方程： \(V(t) = OCV(SOC) - I(t)R_0 - V_1(t) - V_2(t)\)
其中RC网络动态： \(\begin{aligned} \dot{V}_1 &= -\frac{1}{R_1C_1}V_1 + \frac{1}{C_1}I \\ \dot{V}_2 &= -\frac{1}{R_2C_2}V_2 + \frac{1}{C_2}I \end{aligned}\)
步骤1：OCV-SOC曲线辨识
通过低倍率充放电实验，获取准静态OCV-SOC关系： \(OCV(SOC) = a_0 + a_1 \cdot SOC + a_2 \cdot SOC^2 + \cdots + a_n \cdot SOC^n\)
使用最小二乘拟合多项式系数。
步骤2：内阻R0辨识
通过脉冲放电测试，测量瞬时电压跌落： \(R_0 = \frac{\Delta V_{instant}}{\Delta I}\)
步骤3：RC参数辨识
离散化RC网络方程（采样时间 $T_s$）： \(\begin{aligned} V_1(k+1) &= e^{-T_s/(R_1C_1)} V_1(k) + R_1(1 - e^{-T_s/(R_1C_1)}) I(k) \\ V_2(k+1) &= e^{-T_s/(R_2C_2)} V_2(k) + R_2(1 - e^{-T_s/(R_2C_2)}) I(k) \end{aligned}\)
定义参数： \(\begin{aligned} \alpha_1 &= e^{-T_s/(R_1C_1)}, \quad \beta_1 = R_1(1 - \alpha_1) \\ \alpha_2 &= e^{-T_s/(R_2C_2)}, \quad \beta_2 = R_2(1 - \alpha_2) \end{aligned}\)
构造回归模型： \(V(k) = OCV(SOC_k) - R_0 I(k) - V_1(k) - V_2(k)\)
使用递推最小二乘（RLS）在线辨识 $\alpha_i, \beta_i$。
HPPC测试（Hybrid Pulse Power Characterization）：
电流 I
│ 10s放电 40s休息 10s充电
│ ┌─────────┐ ┌─────────┐
│ │ │ │ │
──┼──┘ └────────────────────┘ └──
│
└─────────────────────时间 t
这种脉冲序列能够激发不同时间常数的动态响应。
实际应用中，参数是SOC和温度的函数： \(\theta = f(SOC, T)\)
建立查找表或拟合函数： \(R_0(SOC, T) = r_{00} + r_{01} \cdot SOC + r_{10} \cdot T + r_{11} \cdot SOC \cdot T\)
使用动态应力测试（DST）或联邦城市驾驶循环（FUDS）验证模型：
电压误差统计：
- 最大误差: < 50mV
- 均方根误差(RMSE): < 20mV
- 平均绝对误差(MAE): < 15mV
考虑电池老化，使用自适应算法在线更新参数：
\[\begin{aligned} \hat{\theta}(k) &= \hat{\theta}(k-1) + K(k) \varepsilon(k) \\ K(k) &= \frac{P(k-1)\phi(k)}{\lambda + \phi^T(k)P(k-1)\phi(k)} \\ P(k) &= \frac{1}{\lambda}[P(k-1) - K(k)\phi^T(k)P(k-1)] \end{aligned}\]其中 $\lambda = 0.98$ 是遗忘因子，适应缓慢的老化过程。
Lennart Ljung 是系统辨识理论体系化的关键人物。他在1987年出版的经典著作《System Identification: Theory for the User》成为该领域的权威参考。Ljung的主要贡献包括：
他的工作将系统辨识从分散的技术集合发展成系统的学科，极大推动了控制工程的实践应用。
传统辨识方法需要预先假定模型结构，而实际系统的真实结构往往未知。稀疏辨识通过数据驱动的方式自动发现系统的控制方程，特别适合：
带控制输入的SINDy（SINDyc）： \(\dot{x} = \Theta(x)\Xi_x + \Theta(u)\Xi_u\)
噪声鲁棒的SINDy： 使用积分形式减少噪声影响： \(x(t) - x(0) = \int_0^t \Theta(x(\tau), u(\tau))\xi d\tau\)
约束SINDy： 加入物理约束（如能量守恒）： \(\begin{aligned} \min_\xi &\quad ||\dot{X} - \Theta\xi||_2^2 + \lambda||\xi||_1 \\ \text{s.t.} &\quad C\xi = d \end{aligned}\)
遗传编程（GP）： 通过进化算法搜索函数空间：
深度符号回归： 结合神经网络和符号回归：
模型预测控制的符号化： 将复杂的MPC策略蒸馏为简单的符号表达式： \(u = f_{symbolic}(x, r)\)
优势：
未知动力学的自适应控制：
系统辨识是连接理论与实践的桥梁，本章介绍了从经典到现代的主要辨识方法：
核心概念：
关键公式：
实践要点：
习题12.1 最小二乘估计的性质 给定线性回归模型 $y = X\theta + e$，其中 $e \sim \mathcal{N}(0, \sigma^2I)$。证明最小二乘估计 $\hat{\theta}_{LS}$ 是无偏的，并求其协方差矩阵。
提示：利用 $E[e] = 0$ 和 $\text{Cov}(e) = \sigma^2I$
习题12.2 持续激励条件 考虑一阶系统 $y(k) = ay(k-1) + bu(k-1) + e(k)$。如果输入信号是常数 $u(k) = c$，能否唯一辨识参数 $a$ 和 $b$？说明原因。
提示：考虑信息矩阵的秩
习题12.3 模型阶次选择 给定一组数据，使用不同阶次的ARX模型拟合，得到：
数据长度 $N = 100$。使用AIC准则选择最佳阶次。
提示：$\text{AIC} = N\log(\hat{\sigma}^2) + 2p$
习题12.4 闭环可辨识性分析 系统 $G(s) = \frac{K}{s+a}$ 在单位负反馈下运行。参考输入 $r(t)$ 是白噪声。分析： (a) 能否同时辨识 $K$ 和 $a$？ (b) 如果控制器改为 $C(s) = K_p + \frac{K_i}{s}$，可辨识性如何变化？
提示：分析闭环传递函数的参数化
习题12.5 子空间方法的数值实现 给定输入输出数据，实现MOESP算法的关键步骤。设系统阶次未知，如何从奇异值确定阶次？给出判断准则。
提示：观察奇异值的”膝点”
习题12.6 电池模型的实时辨识 设计一个实时算法，在电动汽车运行中持续更新电池模型参数。要求： (a) 计算复杂度 $O(n^2)$ 以内 (b) 能处理参数突变（如温度变化） (c) 保证数值稳定性
提示：考虑带遗忘因子的RLS或滑动窗口方法
习题12.7 非线性系统的SINDy辨识 考虑Lorenz系统的部分观测（只测量 $x$ 和 $y$）。设计稀疏辨识方案恢复完整动力学。讨论库函数选择和正则化参数调节。
提示：利用时延嵌入重构状态空间
习题12.8 开放性问题：深度学习与经典辨识的结合 讨论如何结合深度学习的表达能力和经典辨识的理论保证。设计一个混合框架，用于辨识具有未知非线性但需要稳定性保证的控制系统。考虑：