上一次学系统辨识已经两年多了，一直也都有在用，但是很多细节概念已经有些遗忘，现在再整体梳理一遍，方便自己也便于一些想入门系统辨识的小伙伴可以参考。
系统辨识三要素
- 数据
- 模型集（线性 or 非线性，连续 or 离散 等）：灰箱模型即结构已知（from physical law）参数带估计；黑箱模型即模型结构需要自己选择且参数没有物理意义
- 准则：用于确定什么样的模型才是best的
模型验证
- 根据先验知识判断是否符合系统的物理特性，比如增益方向，时间常数等；
- 根据未用于训练的实验数据进行验证，比较模型的仿真输出和实际观测输出。
线性系统
线性回归
系统以线性回归表示为
$$
y(k)=\phi^T(k)\theta
$$
其中 $\theta$ 为参数向量，$\phi(k)$ 为回归向量（通常包含过去时刻的输入输出值）。
等价于
$$
y(k)=G(q)u(k),G(q)=\frac{b_0+b_1q^{-1}+...+b_{nb}q^{-nb}}{1+a_1q^{-1}+...+a_{na}q^{-na}}
$$
非线性系统也可以表示成线性回归的形式，比如在回归向量里包含 $y^2(k-1);u(k)y(k-1)$ 等。
最小二乘估计
定义残差（residuals）为 $\epsilon(k,\theta)=y(k)-\hat{y}(k,\theta)=y(k)-\phi^T(k)\theta$
最小二乘估计即：
$$
\hat{\theta}_{LS}=\arg \min_{\theta} \sum_{k=1}^{N} \epsilon^2(k,\theta)=\arg \min_{\theta} ||Y-\Phi \theta||^2
$$
对 $\theta$ 求导，利用极值条件可以得到
$$
\hat{\theta}_{LS}=(\Phi^T\Phi)^{-1}Y=(\sum_{k=1}^N \phi(k)\phi^T(k))^{-1} \sum_{k=1}^N \phi(k)y(k)
$$
在 Matlab 里的实现就是
$$
\hat{\theta}{LS} = \Phi \backslash Y
$$
观察可以知道
$$
\Phi^T \Phi \hat{\theta}_{LS} = \Phi^T Y
$$
利用 Cholesky 分解 $\Phi^T \Phi = LL^T$，其中 $L$ 是下三角矩阵，matlab代码 $L=chol(\Phi^T\Phi,‘lower’)$
因此最小二乘估计的解为：
- 通过 forward substitution 求解 $Lz=\Phi^T Y$
- 通过 backward substitution 求解 $L^T\hat{\theta}_{LS}=z$
利用 QR 分解 $\Phi = QR=Q_1R_1$，其中 $Q=[Q_1,Q_2],R=[R_1,0]^T$，$Q_1^TQ=I$，$R_1$ 为上三角且当 $\Phi$ 满秩时 $R_1$ 对角线元素 $r_{ii}$ 大于0。matlab代码 $qr(\Phi)$。代入可得
$$
\hat{\theta}_{LS}=(\Phi^T\Phi)^{-1}\Phi^TY=((Q_1R_1)^T(Q_1R_1))^{-1}(R_1^TQ_1^TY)=R_1^{-1}Q_1^TY
$$
因此最小二乘估计的解为：
通过 backward substitution 求解
$R_1\hat{\theta}_{LS} = Q_1^TY$
Recursive 线性最小二乘
可以对 $k-1$ 时刻的参数估计 $\hat{\theta}_{LS}(k-1)$ 进行更新得到 $\hat{\theta}_{LS}(k)$。
定义 $P(k)=(\sum_{l=1}^k\phi(l)\phi^T(l))^{-1}$，满足 $P(k)^{-1}=P(k-1)^{-1}+\phi(k)\phi^T(k)$ 。
强推一波公式，记住思路就是利用误差进行修正：
\begin{aligned}
\hat{\theta}_{LS}(k)&=P(k)(\sum_{l=1}^{k-1}\phi(l)y(l) + \phi(k)y(k))\\
&=P(k)(P^{-1}(k-1)\hat{\theta}_{LS}(k-1)+\phi(k)y(k))\\
&=P(k)P^{-1}(k-1)\hat{\theta}_{LS}(k-1)+P(k)\phi(k)y(k)\\
&=P(k)(P(k)^{-1}-\phi(k)\phi^T(k))\hat{\theta}_{LS}(k-1)+P(k)\phi(k)y(k)\\
&=\hat{\theta}_{LS}(k-1) - P(k)\phi(k)\phi^T(k)\hat{\theta}_{LS}(k-1) +P(k)\phi(k)y(k)\\
&=\hat{\theta}_{LS}(k-1) + P(k)\phi(k)(y(k)-\phi^T(k)\hat{\theta}_{LS}(k-1))\\
&=\hat{\theta}_{LS}(k-1) + K(k)\epsilon(k)
\end{aligned}
其中 $K(k)=P(k)\phi(k)$ 可以看做是增益，$\epsilon(k)=y(k)-\phi^T(k)\hat{\theta}_{LS}(k-1)$ 看着是预测误差。
因此 Recursive 线性最小二乘总结为：
\begin{aligned}
\hat{\theta}_{LS}(k) &= \hat{\theta}_{LS}(k-1) + P(k)\phi(k)(y(k)-\phi^T(k)\hat{\theta}_{LS}(k-1))\\
P(k)^{-1}&=P(k-1)^{-1}+\phi(k)\phi^T(k)
\end{aligned}
由于矩阵涉及求逆费事，可以将 $P(k)^{-1}$ 转换为
$$
P(k)=P(k-1)+\frac{P(k-1)\phi(k)\phi^T(k)P(k-1)}{1+\phi^T(k)P(k-1)\phi(k)}
$$
带遗忘因子的 Recursive 线性最小二乘
在目标函数引入遗忘因子 $0< \lambda \le 1$，即
$$
\hat{\theta}=\arg \min_{\theta} \sum_{l=1}^k \lambda^{k-l}(y(l)-\phi^T(l)\theta)^2
$$
$\lambda$ 越小对过去信息忘得越多。同理可以得到递推公式：
\begin{aligned}
\hat{\theta}_{LS}(k) &= \hat{\theta}_{LS}(k-1) + P(k)\phi(k)(y(k)-\phi^T(k)\hat{\theta}_{LS}(k-1))\\
P(k)&=\frac{1}{\lambda} \left[ P(k-1)+\frac{P(k-1)\phi(k)\phi^T(k)P(k-1)}{\lambda+\phi^T(k)P(k-1)\phi(k)}\right]
\end{aligned}
有限脉冲响应模型估计
利用 FIR 模型来描述系统的动态：
$$
\hat{y}(k,g) = \sum_{l=0}^M g(l)u(k-l) = \phi^T(k)g
$$
其中 $\phi(k)=[u(k),u(k-1),…,u(k-M)]^T,g=[g(0),…,g(M)]^T$
- 如果 $M$ 足够大且系统 BIBO 稳定，FIR 模型能够很好近似原系统
- 需要数据量 $N >> M$
最小二乘估计为 $\hat{g}=(\Phi^T\Phi)^{-1}\Phi^TY$
假设真实系统可以被描述为
$$
y(k)=\sum_{l=0}^M g_0(l)u(k-l)+v_0(k)=\phi^T(k
)g_0+v_0(k)
$$
其中 $v_0(k)$ 是零均值且与输入 $u$ 独立的准平稳噪声。
一致性：
\begin{aligned}
\hat{g} &= \left(\frac{1}{N}\sum_{l=1}^N \phi(k)\phi^T(k)\right)^{-1}\frac{1}{N}\sum_{k=1}^N\phi(k)y(k)\\
&=g_0 + \left(\frac{1}{N}\sum_{l=1}^N \phi(k)\phi^T(k)\right)^{-1} \frac{1}{N}\sum_{l=1}^N \phi(k) v_0(k)\\
&=g_0 + R^*E[\phi(k)v_0(k)]\\
&=g_0
\end{aligned}
其中 $R$ 中每一项代表的是自相关或互相关函数的估计，$E[\phi(k)v_0(k)]=E[\phi(k)]\times 0=0$。
ARX 模型估计
模型结构为
$$
y(k)=\frac{B(q)}{A(q)}u(k)+\frac{1}{A(q)}e(k)
$$
输入输出关系为：
$$
y(k)=-a_1y(k-1)-...-a_{n_a}y(k-n_a)+b_1u(k-1)+...+b_{n_b}u(k-n_b) + e(k)
$$
同样可以写成线性回归的形式 $y(k)=\phi^T(k)\theta + e(k)$，同理得到最小二乘估计 $\hat{\theta}_{LS}=(\Phi^T\Phi)^{-1}\Phi^TY$。
假设真实系统为 $y(k)=\phi^T(k)\theta_0+e(k)$，一致性检验可以得到 $lim_{N \rightarrow \infty} \hat{\theta}_{LS}=\theta_0 $
线性参变 ARX 模型估计
模型描述
$$
y(k)=G(q^{-1},p(k))u(k)+v(k)
$$
等价于
$$
A(q^{-1},p(k))y(k)=B(q^{-1},p(k))u(k)+e(k)
$$
其中 $p$ 为调度变量。
通常对回归向量的系数进行相应的参数化，比如多项式
\begin{aligned}
A(q^{-1},p(k))=1+\sum_{i=1}^{n_a}a_i(p(k))q^{-i}\\
a_i(p(k))=a_{i,0}+\sum_{l=1}^{n_l}a_{i,l}p^{l}(k)
\end{aligned}
例子：
$$
y(k)=-[a_{1,0}+a_{1,1}p(k)]y(k-1)+[b_{1,0}+b_{1,1}p(k)]u(k-1)+e(k)
$$
此时有 $\phi^T(k)=[-y(k-1),-p(k)y(k-1),u(k-1),p(k)u(k-1)], \theta=[a{1,0},a{1,1},b{1,0},b{1,1}]^T$，同理还是利用最小二乘进行估计。当 $e(t)$ 为白噪声时，一致性能够得到保证。
输出误差模型的不一致性
模型结构
$$
y(k)=G(q)u(k)+H(q)e(k),H(q)=1
$$
对应的输入输出关系
$$
y_0(k)=-a_1y_0(k-1)-...-a_{n_a}y_0(k-n_a)+b_1u(k-1)+...+b_{n_b}u(k-n_b)\\
y(k)=y_0+e(k)
$$
OE 模型则表示为
\begin{aligned}
y(k)&=\phi^T(k)\theta + e(k)+a_1e(k-1)+...a_{n_a}e(k-n_a)\\
y(k)&=\phi^T(k)\theta + v(k)
\end{aligned}
注意！这里 $v(k)$ 不是白噪声了！
一致性：
\begin{aligned}
\hat{\theta}_{LS}&=(\Phi^T\Phi)^{-1}\Phi^TY=\left(\frac{1}{N}\sum_{k=1}^N\phi(k)\phi^T(k) \right)^{-1}\frac{1}{N}\sum_{k=1}^N \phi(k)y(k)\\
&=\left(\frac{1}{N}\sum_{k=1}^N\phi(k)\phi^T(k) \right)^{-1}\frac{1}{N}\sum_{k=1}^N \phi(k) (\phi^T(k)\theta_0 + v(k)) \\
&=\theta_0 + \left(\frac{1}{N}\sum_{k=1}^N\phi(k)\phi^T(k) \right)^{-1}\frac{1}{N}\sum_{k=1}^N \phi(k) v(k)
\end{aligned}
由于 $\phi(k)$ 和 $v(k)$ 相关，所以 $lim_{N \rightarrow \infty} \frac{1}{N} \sum_{k=1}^N \phi(k)v(k) \neq 0$。
辅助变量法
给定模型结构 $A(q^{-1}y(k))=B(q^{-1})u(k)+v(k)$，只有当 $v(k)$ 和回归向量不相关时，一致性才能得到保证。
选择与回归向量 $\phi(k)$ 维度相同的辅助变量 $z(k)$ 使得
$$
E[z(k)v(K)]=0
$$
修正最小二乘估计为
$$
\hat{\theta}_{IV}=(Z^T\Phi)^{-1}Z^TY=\left(\sum_{k=1}^N z(k)\phi^T(k) \right)^{-1}\sum_{k=1}^N z(k)y(k)
$$
一致性：
$$
\hat{\theta}_{IV}=\theta_0 + \left(\frac{1}{N}\sum_{k=1}^Nz(k)\phi^T(k) \right)^{-1}\frac{1}{N}\sum_{k=1}^N z(k) v(k)
$$
因此 $z(k)$ 应该和 $\phi(k)$ 相关否则 $R(N)$ 会收敛到零矩阵。
辅助变量的选择
$$
z(k)=[-y_0(k-1),...,-y_0(k-n_a),u(k-1),...,u(k-n_b)]^T
$$
也就是说，$z(k)$ 应该尽可能和无噪的回归向量相关，但是和 $v(k)$ 无关。
辅助变量法总结为：
- 最小二乘估计 $\hat{\theta}$ （biased）
- 开环仿真
$$
\hat{y}(k)=[-\hat{y}(k-1),…,-\hat{y}(k-n_a),u(k-1),…,u(k-n_b)]^T\hat{\theta}
$$
- 构造辅助变量
$$
z(k)=[-\hat{y}(k-1),…,-\hat{y}(k-n_a),u(k-1),…,u(k-n_b)]^T
$$
- 通过 IV 估计 $\hat{\theta}$ （consistent）
- 从第2步重复直到收敛
预报误差法 PEM
线性预报器
$$
\hat{y}(k|k-1;\theta) = L_y(q^{-1},\theta)y(k) + L_u(q^{-1},\theta) u(k); L_y(0,\theta)=0, L_u(0,\theta)=0
$$
只由过去的输入输出数据决定。
预报误差为
$$
\epsilon (k,\theta) = y(k)-\hat{y}(k|k-1,\theta)
$$
待估计的参数应该使得预报误差尽可能小。
- 模型结构，对 $G(q,\theta),H(q,\theta)$ 的参数化；
- 预报器，定义 $L_y,L_u$；
- 准则，预报误差的 scalar 函数 $V_N(\theta)$；
- 估计 $\hat{\theta}=\arg \min_{\theta} V_N(\theta)$
单变量模型结构
- ARX
$$y(k)=\frac{B(q,\theta)}{A(q,\theta)}u(k)+\frac{1}{A(q,\theta)}e(k)$$
- ARMAX
$$y(k)=\frac{B(q,\theta)}{A(q,\theta)}u(k)+\frac{C(q,\theta)}{A(q,\theta)}e(k)$$
- OE
$$y(k)=\frac{B(q,\theta)}{A(q,\theta)}u(k)+e(k)$$
- BJ
$$y(k)=\frac{B(q,\theta)}{A(q,\theta)}u(k)+\frac{C(q,\theta)}{D(q,\theta)}e(k)$$
预报滤波器
应该选择滤波器 $L_y,L_u$ 使得预报误差的方差最小。
\begin{aligned}
y(k)&=Gu(k)+He(k)=Gu(k)+(H-I)e(k)+e(k)\\
&=Gu(k)+(H-I)H^{-1}[y(k)-Gu(k)]+e(k)\\
&=[(I-H^{-1})y(k)+H^{-1}Gu(k)]+e(k)\\
&=z(k)+e(k)
\end{aligned}
注意 $z(k),e(k)$ 不相关。
假设 $y^*(k)$ 为 $y(k)$ 的任意预测，$z(k)$ 的最优性：
\begin{aligned}
E[(y-y^*)(y-y^*)^T]&=E[(z+e-y^*)(z+e-y^*)^T]\\
&=E[(z-y^*)(z-y^*)^T]+E[ee^T] \ge \Lambda_e
\end{aligned}
因此 $z(k)$ 是最优预报，$e(k)$ 是最优预报误差。
\begin{aligned}
\hat{y}(k|k-1,\theta)&=\left(I-H^{-1}(q^{-1},\theta)\right)y(k)+H^{-1}(q^{-1},\theta)G(q^{-1},\theta)u(k)\\
\epsilon (k,\theta)&=e(k)=H^{-1}(q^{-1},\theta)\left(y(k)-G(q^{-1},\theta)u(k)\right)
\end{aligned}
最小化准则
采样协方差矩阵
$$
R_N(\theta)=\frac{1}{N}\sum_{k=1}^N \epsilon(k,\theta)\epsilon^T(k,\theta)
$$
- 如果 $y$ 是标量，$R_N(\theta)$ 可以直接作为 $V_N(\theta)$ 用于最小化。
- 如果是多变量，我们可以最小化
$$
V_N(\theta)=h(R_N(\theta))
$$
其中，$h$ 为定义在正半定矩阵集合上的连续的单调递增函数满足 $h(Q+\Delta Q)\ge h(Q),\forall Q,\Delta Q \ge 0$，比如矩阵的迹 $V_N(\theta)=h(R_N(\theta))=tr(R_N(\theta))$。