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

前言

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

最小二乘估计

如下的系统辨识模型： { 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​

参数的最小二乘估计可表示为： θ ^ L S = ( H t T H t ) − 1 H t T Y t \hat\theta_{LS}=(H_t^TH_t)^{-1}H_t^TY_t θ^LS​=(HtT​Ht​)−1HtT​Yt​

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

递归最小二乘估计 原理

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

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)

于是，系统参数的递归最小二乘可表示为： θ ^ ( t ) = θ ^ ( t − 1 ) + L ( t ) [ y ( t ) − φ T ( t ) θ ^ ( t − 1 ) ] ] L ( t ) = P ( t − 1 ) φ ( t ) 1 + φ T ( t ) P ( t − 1 ) φ ( t ) P ( t ) = [ I − L ( t ) φ T ( t ) ] P ( t − 1 ) , P ( 0 ) = p 0 I , p 0 = inf ⁡ φ ( t ) = [ . . . ] \hat\theta(t)=\hat\theta(t-1)+L(t)[y(t)-\varphi^T(t)\hat\theta(t-1)]] \\ \quad \\ L(t)=\frac{P(t-1)\varphi(t)}{1+\varphi^T(t)P(t-1)\varphi(t)} \\ \quad \\ P(t)=[I-L(t)\varphi^T(t)]P(t-1), P(0)=p_0I, p_0=\inf \\ \quad \\ \varphi(t)=[...] θ^(t)=θ^(t−1)+L(t)[y(t)−φT(t)θ^(t−1)]]L(t)=1+φT(t)P(t−1)φ(t)P(t−1)φ(t)​P(t)=[I−L(t)φT(t)]P(t−1),P(0)=p0​I,p0​=infφ(t)=[...]

新息，残差

新息 e ( t ) = y ( t ) − φ T ( t ) θ ^ ( t − 1 ) e(t)=y(t)-\varphi^T(t)\hat\theta(t-1) e(t)=y(t)−φT(t)θ^(t−1) 残差 ε ( t ) = y ( t ) − φ T ( t ) θ ^ ( t ) \varepsilon(t)=y(t)-\varphi^T(t)\hat\theta(t) ε(t)=y(t)−φT(t)θ^(t)

数据饱和 现象

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

原因

P ( t ) = [ P − 1 ( t − 1 ) + φ ( t ) φ T ( t ) ] − 1 = > P ( t ) ≤ P ( t − 1 ) ( α I + 1 / p 0 ) t ≤ P − 1 ( t ) = P − 1 ( 0 ) + ∑ i = 1 t φ ( i ) φ T ( i ) ≤ ( β I + 1 / p 0 ) t lim ⁡ t − > inf ⁡ P − 1 ( t ) = 0 P(t)=[P^{-1}(t-1)+\varphi(t)\varphi^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\varphi(i)\varphi^T(i) \le (\beta I+1/p_0)t \\ \lim_{t->\inf}P^{-1}(t) = 0 P(t)=[P−1(t−1)+φ(t)φT(t)]−1=>P(t)≤P(t−1)(αI+1/p0​)t≤P−1(t)=P−1(0)+i=1∑t​φ(i)φT(i)≤(βI+1/p0​)tt−>inflim​P−1(t)=0

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

θ ^ ( t ) = θ ^ ( t − 1 ) + L ( t ) [ y ( t ) − φ T ( t ) θ ^ ( t − 1 ) ] ] L ( t ) = P ( t − 1 ) φ ( t ) λ + φ T ( t ) P ( t − 1 ) φ ( t ) P ( t ) = 1 λ [ I − L ( t ) φ T ( t ) ] P ( t − 1 ) , P ( 0 ) = p 0 I , p 0 = inf ⁡ φ ( t ) = [ . . . ] \hat\theta(t)=\hat\theta(t-1)+L(t)[y(t)-\varphi^T(t)\hat\theta(t-1)]] \\ \quad \\ L(t)=\frac{P(t-1)\varphi(t)}{\lambda+\varphi^T(t)P(t-1)\varphi(t)} \\ \quad \\ P(t)=\frac{1}{\lambda}[I-L(t)\varphi^T(t)]P(t-1), P(0)=p_0I, p_0=\inf \\ \quad \\ \varphi(t)=[...] θ^(t)=θ^(t−1)+L(t)[y(t)−φT(t)θ^(t−1)]]L(t)=λ+φT(t)P(t−1)φ(t)P(t−1)φ(t)​P(t)=λ1​[I−L(t)φT(t)]P(t−1),P(0)=p0​I,p0​=infφ(t)=[...]

后记

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