最小二乘系统辨识课 中篇：递归最小二乘

前言

本篇记录系统辨识课中的递归最小二乘原理。

最小二乘估计

如下的系统辨识模型：
{ y ( 1 ) = φ T ( 1 ) θ + v ( 1 ) y ( 2 ) = φ T ( 2 ) θ + v ( 2 ) . . . y ( t ) = φ T ( t ) θ + v ( t ) 设 ： Y t = [ y ( 1 ) y ( 2 ) . . . y ( t ) ] T H t = [ φ ( 1 ) φ ( 2 ) . . . φ ( t ) ] T V t = [ v ( 1 ) v ( 2 ) . . . v ( n ) ] T 有 ： Y t = H t θ + V t \left \{

$$\begin{aligned} y(1)&=\varphi^T(1)\theta+v(1) \\ y(2)&=\varphi^T(2)\theta+v(2) \\ &... \\ y(t)&=\varphi^T(t)\theta+v(t) \end{aligned}$$ \right . \\ \quad \\ $$\begin{aligned} 设： Y_t&=[y(1)\quad y(2)\quad ... \quad y(t)]^T \\ H_t&=[\varphi(1) \quad \varphi(2) \quad ... \quad \varphi(t)]^T \\ V_t&=[v(1) \quad v(2) \quad ... \quad v(n)]^T \end{aligned}$$ \\ \quad \\ 有： Y_t=H_t\theta +V_t ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧​y(1)y(2)y(t)​=φT(1)θ+v(1)=φT(2)θ+v(2)...=φT(t)θ+v(t)​设：Yt​Ht​Vt​​=[y(1)y(2)...y(t)]T=[φ(1)φ(2)...φ(t)]T=[v(1)v(2)...v(n)]T​有：Yt​=Ht​θ+Vt​

参数的最小二乘估计可表示为：
$${\hat{\theta }}_{LS}=(H_t^T{H}_t{)}^{-1}{H}_t^T{Y}_t$$

上面的方法在每个新时刻的输入输出数据更新时，都要重估计，计算负载很大。

递归最小二乘估计

原理

递归最小二乘估计采用了增量思路，推导如下：

P − 1 ( t ) = H t T H t = ∑ i = 1 t φ ( i ) φ T ( i ) = ∑ i = 1 t − 1 φ ( i ) φ T ( i ) + φ ( t ) φ T ( t ) = P − 1 ( t − 1 ) + φ ( t ) φ T ( t ) θ ^ ( t ) = ( H t T H t ) − 1 H t T Y t = P ( t ) H t T Y t = P ( t ) [ H t − 1 T φ ( t ) ] [ Y t − 1 T y ( t ) ] T = P ( t ) H t − 1 T Y t − 1 T + P ( t ) φ ( t ) y ( t ) = P ( t ) P − 1 ( t − 1 ) P ( t − 1 ) H t − 1 T Y t − 1 T + P ( t ) φ ( t ) y ( t ) = P ( t ) P − 1 ( t − 1 ) θ ^ ( t − 1 ) + P ( t ) φ ( t ) y ( t ) = P ( t ) [ P − 1 ( t ) − φ ( t ) φ T ( t ) ] θ ^ ( t − 1 ) + P ( t ) φ ( t ) y ( t ) = θ ^ ( t − 1 ) − P ( t ) φ ( t ) φ T ( t ) θ ^ ( t − 1 ) + P ( t ) φ ( t ) y ( t ) = θ ^ ( t − 1 ) + P ( t ) φ ( t ) [ y ( t ) − φ T ( t ) θ ^ ( t − 1 ) ] 引 理 ： 矩 阵 A , B , C ， 若 ( I + C A − 1 B ) 与 A 都 可 逆 ， 则 ： ( A + B C ) − 1 = A − 1 − A − 1 B ( I + C A − 1 B ) − 1 C A − 1 P − 1 ( t ) = P − 1 ( t − 1 ) + φ ( t ) φ T ( t ) 把 后 面 的 式 子 取 逆 ， 根 据 引 理 有 ： P ( t ) = P ( t − 1 ) − P ( t − 1 ) φ ( t ) φ T ( t ) P ( t − 1 ) 1 + φ T ( t ) P ( t − 1 ) φ ( t ) P ( t ) φ ( t ) = P ( t − 1 ) φ ( t ) 1 + φ T ( t ) P ( t − 1 ) φ ( t ) L ( t ) = P ( t ) φ ( t ) P^{-1}(t)=H_t^TH_t=\sum_{i=1}^t \varphi(i)\varphi^T(i)=\sum_{i=1}^{t-1} \varphi(i)\varphi^T(i) + \varphi(t)\varphi^T(t) = P^{-1}(t-1)+ \varphi(t)\varphi^T(t) \\ \quad \\

$$\begin{aligned} \hat\theta(t)&=(H_t^TH_t)^{-1}H_t^TY_t \\ &=P(t)H_t^TY_t \\ &=P(t)\begin{bmatrix} H_{t-1}^T & \varphi(t) \end{bmatrix} \begin{bmatrix} Y_{t-1}^T & y(t) \end{bmatrix}^T \\ &=P(t)H_{t-1}^TY_{t-1}^T + P(t)\varphi(t)y(t) \\ &=P(t)P^{-1}(t-1)P(t-1)H_{t-1}^TY_{t-1}^T + P(t)\varphi(t)y(t) \\ &=P(t)P^{-1}(t-1)\hat \theta(t-1)+P(t)\varphi(t)y(t) \\ &=P(t)[P^{-1}(t)- \varphi(t)\varphi^T(t) ]\hat\theta(t-1)+P(t)\varphi(t)y(t) \\ &=\hat\theta(t-1)-P(t)\varphi(t)\varphi^T(t)\hat\theta(t-1)+P(t)\varphi(t)y(t) \\ &=\hat\theta(t-1)+P(t)\varphi(t)[y(t)-\varphi^T(t)\hat\theta(t-1)] \end{aligned}$$ \\ \quad \\ 引理：矩阵A,B,C，若(I+CA^{-1}B)与A都可逆，则：\\ (A+BC)^{-1}=A^{-1}-A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1} \\ \quad \\ P^{-1}(t)= P^{-1}(t-1)+ \varphi(t)\varphi^T(t) \\ \quad \\ 把后面的式子取逆，根据引理有：\\ P(t)=P(t-1)-\frac{P(t-1)\varphi(t)\varphi^T(t)P(t-1)}{1+\varphi^T(t)P(t-1)\varphi(t)} \\ \quad \\ P(t)\varphi(t)= \frac{P(t-1)\varphi(t)}{1+\varphi^T(t)P(t-1)\varphi(t)} \\ \quad \\ L(t)=P(t)\varphi(t) P−1(t)=HtT​Ht​=i=1∑t​φ(i)φT(i)=i=1∑t−1​φ(i)φT(i)+φ(t)φT(t)=P−1(t−1)+φ(t)φT(t)θ^(t)​=(HtT​Ht​)−1HtT​Yt​=P(t)HtT​Yt​=P(t)[Ht−1T​​φ(t)​][Yt−1T​​y(t)​]T=P(t)Ht−1T​Yt−1T​+P(t)φ(t)y(t)=P(t)P−1(t−1)P(t−1)Ht−1T​Yt−1T​+P(t)φ(t)y(t)=P(t)P−1(t−1)θ^(t−1)+P(t)φ(t)y(t)=P(t)[P−1(t)−φ(t)φT(t)]θ^(t−1)+P(t)φ(t)y(t)=θ^(t−1)−P(t)φ(t)φT(t)θ^(t−1)+P(t)φ(t)y(t)=θ^(t−1)+P(t)φ(t)[y(t)−φT(t)θ^(t−1)]​引理：矩阵A,B,C，若(I+CA−1B)与A都可逆，则：(A+BC)−1=A−1−A−1B(I+CA−1B)−1CA−1P−1(t)=P−1(t−1)+φ(t)φT(t)把后面的式子取逆，根据引理有：P(t)=P(t−1)−1+φT(t)P(t−1)φ(t)P(t−1)φ(t)φT(t)P(t−1)​P(t)φ(t)=1+φT(t)P(t−1)φ(t)P(t−1)φ(t)​L(t)=P(t)φ(t)

于是，系统参数的递归最小二乘可表示为：
$$\begin{gathered} \hat{\theta }(t)=\hat{\theta }(t-1)+L(t)[y(t)-{\phi }^T(t)\hat{\theta }(t-1)]] \\ L(t)=\frac{P(t-1)\phi (t)}{1+{\phi }^T(t)P(t-1)\phi (t)} \\ P(t)=[I-L(t){\phi }^T(t)]P(t-1),P(0)=p_{0}I,p_0=\inf \\ \phi (t)=[...] \end{gathered}$$

新息，残差

新息$e(t)=y(t)-{\phi }^T(t)\hat{\theta }(t-1)$
残差$\epsilon (t)=y(t)-{\phi }^T(t)\hat{\theta }(t)$

数据饱和

现象

当输入输出时间序列足够长时，新的输入输出对递归最小二乘估计结果没有贡献的现象，称为数据饱和。

原因

$$\begin{gathered} P(t)=[P^{-1}(t-1)+\phi (t){\phi }^T(t){]}^{-1}=>P(t)\le P(t-1) \\ (\alpha I+1/p_0)t\le P^{-1}(t)=P^{-1}(0)+\sum _{i=1}^t\phi (i){\phi }^T(i)\le (\beta I+1/p_0)t \\ \underset{t->\inf }{\lim }{P}^{-1}(t)=0 \end{gathered}$$

解决方法：带遗忘因子的递归最小二乘

$$\begin{gathered} \hat{\theta }(t)=\hat{\theta }(t-1)+L(t)[y(t)-{\phi }^T(t)\hat{\theta }(t-1)]] \\ L(t)=\frac{P(t-1)\phi (t)}{\lambda +{\phi }^T(t)P(t-1)\phi (t)} \\ P(t)=\frac{1}{\lambda }[I-L(t){\phi }^T(t)]P(t-1),P(0)=p_{0}I,p_0=\inf \\ \phi (t)=[...] \end{gathered}$$

后记

本篇记录了递归最小二乘的推导，数据饱和现象以及解决方法。下篇将记录最小二乘中，${\phi }^T(t)$的获取方法。