Dingnuooo's Notes

Back
16 min read
note /
control-science

控制科学笔记 - 子空间辨识

《控制科学Ⅰ》期末速通笔记。课件太抽象了 咱来点接地气的
views | comments

本章解决这样一个问题：对于一个 LTI 离散状态空间方程模型（注意与系统建模中的区别）

{xk+1=Axk+Bukyk=Cxk+Duk\left\{ \begin{aligned} x_{k+1}&= Ax_k+Bu_k\\ y_k&= Cx_k+Du_k \end{aligned} \right.

给定系统若干个输入 u0,⋯ ,us−1u_0,\cdots,u_{s-1}{} 和对应的 ss 个输出 yky_k（输入为 mm 维列向量、输出为 ll 维列向量），求出：

  • 系统的阶数 nn
  • 四个矩阵 AA、BB、CC、DD（相似变换意义下等价即可）

子空间方程#

推导#

假定当前时刻 k=ik=i. 取 k=0∼(i−1)k=0\sim (i-1) 这 ii 个下标，也就是相对于当前时刻的所有“过去状态”：

∣ x1=Ax0+Bu0x2=Ax1+Bu1=A2x0+ABu0+Bu1x3=Ax2+Bu2=A3x0+A2Bu0+ABu1+Bu2⋯xi=Aix0+∑k=0(i−1)A(i−1)−kBuk∣ y0=Cx0+Du0y1=Cx1+Du1=CAx0+CBu0+Du1y2=Cx2+Du2=CA2x0+CABu0+CBu1+Du2⋯yi−1=CAi−1x0+∑k=0(i−2)CA(i−2)−kBuk+Dui−1\left|\ \begin{aligned} x_1&= Ax_0+Bu_0\\ x_2&= Ax_1+Bu_1\\ &= A^2x_0+ABu_0+Bu_1\\ x_3&= Ax_2+Bu_2\\ &= A^3x_0+A^2Bu_0+ABu_1+Bu_2\\ &\cdots\\ x_i&= A^ix_0+\sum\limits_{k=0}^{(i-1)}A^{(i-1)-k}Bu_k \end{aligned} \right. \quad \left|\ \begin{aligned} y_0 &= C x_0 + D u_0\\ y_1 &= Cx_1+Du_1\\ &= C A x_0 + C B u_0 + D u_1\\ y_2 &= Cx_2+Du_2\\ &= C A^2 x_0 + C A B u_0 + C B u_1 + D u_2\\ &\cdots\\ y_{i-1} &= CA^{i-1}x_0+\sum\limits_{k=0}^{(i-2)}CA^{(i-2)-k}Bu_k+Du_{i-1} \end{aligned} \right.

把 yy 的部分写成矩阵形式

[y0y1⋮yi−2yi−1]=[CCA⋮CAi−2CAi−1]x0+[D0⋯00CBD⋯00⋮⋮⋱⋮⋮CAi−3BCAi−4B⋯D0CAi−2BCAi−3B⋯CBD][u0u1⋮ui−2ui−1]\begin{bmatrix} y_0\\y_1\\\vdots\\y_{i-2}\\y_{i-1} \end{bmatrix} = \begin{bmatrix} C\\ CA\\ \vdots\\ CA^{i-2}\\ CA^{i-1} \end{bmatrix} x_0 + \begin{bmatrix} D & 0 & \cdots & 0 & 0\\ CB & D & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ CA^{i-3}B & CA^{i-4}B & \cdots & D & 0\\ CA^{i-2}B & CA^{i-3}B & \cdots & CB & D \end{bmatrix} \begin{bmatrix} u_0\\u_1\\\vdots\\u_{i-2}\\u_{i-1} \end{bmatrix}

记 Γi\Gamma_i 为 x0x_0 前面的那一坨大矩阵、记 HiH_i 为 uu 前面的那一坨系数矩阵. 于是可将表达式压缩为

[y0y1⋮yi−1]=Γix0+Hi[u0u1⋮ui−1]\begin{bmatrix} y_0\\y_1\\ \vdots\\ y_{i-1} \end{bmatrix} = \Gamma_ix_0 + H_i \begin{bmatrix} u_0\\u_1\\\vdots\\u_{i-1} \end{bmatrix}

刚才我们的 kk 是从 0∼(i−1)0\sim(i-1). 如果我们把下标整体往后移一位，即 kk 从 1∼i1\sim i，这样等号左边的 yy 矩阵就是 y1∼yiy_1\sim y_i，右边就是 x1x_1，uu 矩阵就是 u1∼uiu_1\sim u_i. 很巧妙的在于，Γi\Gamma_i 和 HiH_i 并不会变，因为这两个系数矩阵实际上表达的是“相对于第一项的关系”，而量的个数没有变（都是 ii），也就是说相对于 x1x_1，后面的 ii 个数都有相同的矩阵表示. 同理，整体移动两位、三位……，两个系数矩阵都是一样的. 因此，把多个 yy 矩阵横着拼到一起时，系数矩阵可以提取出来.

从 00 拼到 j−1j-1（jj 只是一个参数，表示我们将使用多少数据进行分析. 有点像以前的迭代次数）

[y0y1⋯yj−1y1y2⋯yj⋮⋮⋱⋮yi−1yi⋯yi+j−2]=Γi[x0x1⋯xj−1]+Hi[u0u1⋯uj−1u1u2⋯uj⋮⋮⋱⋮ui−1ui⋯ui+j−2]\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} = \Gamma_i \begin{bmatrix} x_0 & x_1 & \cdots & x_{j-1} \end{bmatrix} + H_i \begin{bmatrix} u_0 & u_1 & \cdots & u_{j-1} \\ u_1 & u_2 & \cdots & u_{j} \\ \vdots & \vdots & \ddots & \vdots \\ u_{i-1} & u_{i} & \cdots & u_{i+j-2} \end{bmatrix}

我们再次压缩表达式为

Yp=ΓiXp+HiUpY_p=\Gamma_iX_p+H_iU_p

这里的 pp 就是 past“过去状态”的意思。

现在考虑“未来状态”，特别地，是未来状态相对于当前状态（即 k=ik=i）的情况，因此下标又要整体移动 ii，也即从 ii 拼到 i+j−1i+j-1：

[yiyi+1⋯yi+j−1yi+1yi+2⋯yi+j⋮⋮⋱⋮y2i−1y2i⋯y2i+j−2]=Γi[xixi+1⋯xi+j−1]+Hi[uiui+1⋯ui+j−1ui+1ui+2⋯ui+j⋮⋮⋱⋮u2i−1u2i⋯u2i+j−2]\begin{bmatrix} y_i & y_{i+1} & \cdots & y_{i+j-1} \\ y_{i+1} & y_{i+2} & \cdots & y_{i+j} \\ \vdots & \vdots & \ddots & \vdots \\ y_{2i-1} & y_{2i} & \cdots & y_{2i+j-2} \end{bmatrix} = \Gamma_i \begin{bmatrix} x_i & x_{i+1} & \cdots & x_{i+j-1} \end{bmatrix} + H_i \begin{bmatrix} u_i & u_{i+1} & \cdots & u_{i+j-1} \\ u_{i+1} & u_{i+2} & \cdots & u_{i+j} \\ \vdots & \vdots & \ddots & \vdots \\ u_{2i-1} & u_{2i} & \cdots & u_{2i+j-2} \end{bmatrix}

压缩表达式为

Yf=ΓiXf+HiUfY_f=\Gamma_iX_f+H_iU_f

这里的 ff 就是 future 的意思。

过去和未来的情况都有了，我们使用一开始算出来的各个 xx 的表达式，在过去与未来之间建立联系. 将 ii 时刻和 00 时刻之间的关系 xi=Aix0+∑k=0(i−1)A(i−1)−kBukx_i= A^ix_0+\sum\limits_{k=0}^{(i-1)}A^{(i-1)-k}Bu_k 写成矩阵形式、并压缩表达式为：

xi=Aix0+[Ai−1BAi−2B⋯ABB][u0u1⋮ui−2ui−1]=Aix0+Δi[u0u1⋮ui−1]x_i= A^ix_0+\begin{bmatrix}A^{i-1}B&A^{i-2}B&\cdots&AB&B\end{bmatrix}\begin{bmatrix}u_0\\u_1\\\vdots\\u_{i-2}\\u_{i-1}\end{bmatrix}=A^ix_0+\Delta_i\begin{bmatrix}u_0\\u_1\\\vdots\\u_{i-1}\end{bmatrix}

同样地，依然可以把下标整体移动，而两个系数矩阵不变. 然后再把 00 到 j−1j-1 拼到一起：

[xixi+1⋯xi+j−1]=Ai[x0x1⋯xj−1]+Δi[u0u1⋯uj−1u1u2⋯uj⋮⋮⋱⋮ui−1ui⋯ui+j−2]\begin{bmatrix}x_i&x_{i+1}&\cdots&x_{i+j-1}\end{bmatrix}=A^i\begin{bmatrix}x_0&x_1&\cdots&x_{j-1}\end{bmatrix}+\Delta_i\begin{bmatrix}u_0&u_1&\cdots&u_{j-1}\\u_1&u_2&\cdots&u_j\\\vdots&\vdots&\ddots&\vdots\\u_{i-1}&u_i&\cdots&u_{i+j-2}\end{bmatrix}

自然地，过去与未来 xx 被分别打包在一起了，因此可以再次压缩表达式为

Xf=AiXp+ΔiUpX_f=A^iX_p+\Delta_iU_p

最后一步. 我们实际上只想知道输入 YY 和输出 UU，状态变量可以进一步消去

  • 第①式移项后代入第③式消 XpX_p，压缩表达式为
Xf=AiΓi†(Yp−HiUp)+ΔiUp=[AiΓi†,Δi−AiΓi†Hi][YpUp]Xf=LpWp\begin{aligned}X_f&= A^i\Gamma_i^\dagger(Y_p-H_iU_p)+\Delta_iU_p\\&= \begin{bmatrix}A^i\Gamma_i^\dagger&,&\Delta_i-A^i\Gamma_i^\dagger H_i\end{bmatrix}\begin{bmatrix}Y_p\\U_p\end{bmatrix}\\X_f&= L_pW_p\end{aligned}
  • 相应地第②式改为 Yf=ΓiXf+HiUf=ΓiLpWp+HiUfY_f=\Gamma_iX_f+H_iU_f=\Gamma_iL_pW_p+H_iU_f

小结#

对一个 mm 维输入、ll 维输出的状态空间方程模型

{xk+1=Axk+Bukyk=Cxk+Duk\left\{ \begin{aligned} x_{k+1}&= Ax_k+Bu_k\\ y_k&= Cx_k+Du_k \end{aligned} \right.

我们有三条子空间辨识方程：

其中

Xp=[x0x1⋯xj−1]Xf=[xixi+1⋯xi+j−1]\begin{aligned} X_p&= \begin{bmatrix}x_0&x_1&\cdots&x_{j-1}\end{bmatrix}\\ X_f&= \begin{bmatrix}x_i&x_{i+1}&\cdots&x_{i+j-1}\end{bmatrix} \end{aligned} Yp=[y0y1⋯yj−1y1y2⋯yj⋮⋮⋱⋮yi−1yi⋯yi+j−2]Yf=[yiyi+1⋯yi+j−1yi+1yi+2⋯yi+j⋮⋮⋱⋮y2i−1y2i⋯y2i+j−2]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}\quad Y_f=\begin{bmatrix} y_i & y_{i+1} & \cdots & y_{i+j-1} \\ y_{i+1} & y_{i+2} & \cdots & y_{i+j} \\ \vdots & \vdots & \ddots & \vdots \\ y_{2i-1} & y_{2i} & \cdots & y_{2i+j-2} \end{bmatrix} Up=[u0u1⋯uj−1u1u2⋯uj⋮⋮⋱⋮ui−1ui⋯ui+j−2]Uf=[uiui+1⋯ui+j−1ui+1ui+2⋯ui+j⋮⋮⋱⋮u2i−1u2i⋯u2i+j−2]U_p=\begin{bmatrix}u_0&u_1&\cdots&u_{j-1}\\u_1&u_2&\cdots&u_j\\\vdots&\vdots&\ddots&\vdots\\u_{i-1}&u_i&\cdots&u_{i+j-2}\end{bmatrix}\quad U_f=\begin{bmatrix} u_i & u_{i+1} & \cdots & u_{i+j-1} \\ u_{i+1} & u_{i+2} & \cdots & u_{i+j} \\ \vdots & \vdots & \ddots & \vdots \\ u_{2i-1} & u_{2i} & \cdots & u_{2i+j-2} \end{bmatrix} Γi=[CCA⋮CAi−2CAi−1]Hi=[D0⋯00CBD⋯00⋮⋮⋱⋮⋮CAi−3BCAi−4B⋯D0CAi−2BCAi−3B⋯CBD]\Gamma_i=\begin{bmatrix} C\\ CA\\ \vdots\\ CA^{i-2}\\ CA^{i-1} \end{bmatrix}\quad H_i=\begin{bmatrix} D & 0 & \cdots & 0 & 0\\ CB & D & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ CA^{i-3}B & CA^{i-4}B & \cdots & D & 0\\ CA^{i-2}B & CA^{i-3}B & \cdots & CB & D \end{bmatrix} Δi=[Ai−1BAi−2B⋯ABB]\Delta_i=\begin{bmatrix}A^{i-1}B&A^{i-2}B&\cdots&AB&B\end{bmatrix}

消去状态变量的形式：

其中

Lp=[AiΓi†,Δi−AiΓi†Hi]Wp=[YpUp]L_p=\begin{bmatrix}A^i\Gamma_i^\dagger&,&\Delta_i-A^i\Gamma_i^\dagger H_i\end{bmatrix}\quad\quad W_p=\begin{bmatrix}Y_p\\U_p\end{bmatrix}

三条假设与两条推论#

为了让辨识算法有效，需要假定系统足够好，即满足：

若系统满足上述三个条件，有两个推论

证：Xf=AiXp+ΔiUp=[AiΔi][XpUp]X_f=A^iX_p+\Delta_iU_p=\begin{bmatrix}A^i & \Delta_i\end{bmatrix}\begin{bmatrix}X_p \\U_p\end{bmatrix}
由第二条，r(Up)=im{\rm r}(U_p) = im 行满秩；由第三条，两行空间不交，故 [XpUp]\begin{bmatrix}X_p \\U_p\end{bmatrix} 亦行满秩，故

r(Xf)=r([AiΔi])≥r(Δi)=n{\rm r}(X_f) = {\rm r}(\begin{bmatrix}A^i&\Delta_i\end{bmatrix})\geq {\rm r}(\Delta_i)=n

最后一步是因为 Δi\Delta_i 行满秩
又 span(Xf){\rm span}(X_f) 是 n 维状态空间的子空间，于是 r(Xf)≤n{\rm r}(X_f) \leq n. 故 r(Xf)=n{\rm r}(X_f) = n

证：(1) 证明 span(Xf)⊆span(Wp)∩span(Wf)\text{span}(X_f) \subseteq \text{span}(W_p) \cap \text{span}(W_f)

由于系统是可达和可观测，系统的状态 XfX_f 可以由输入输出数据唯一确定。因此，XfX_f 必定是 WpW_p 和 WfW_f 张成空间中的一个元素，所以 Xf∈span(Wp)X_f \in \text{span}(W_p) 且 Xf∈span(Wf)X_f \in \text{span}(W_f)，也即

span(Xf)⊆span(Wp)∩span(Wf)\text{span}(X_f) \subseteq \text{span}(W_p) \cap \text{span}(W_f)

(2) 证明 span(Wp)∩span(Wf)⊆span(Xf)\text{span}(W_p) \cap \text{span}(W_f) \subseteq \text{span}(X_f)

∀ v∈span(Wp)∩span(Wf)\forall\,v \in \text{span}(W_p) \cap \text{span}(W_f)，∃ a, b  s.t.\exists\,a,\,b\,\ \text{s.t.}{}

v=Wp⋅a=Wf⋅bv = W_p\cdot a = W_f\cdot b

而 WpW_p 和 WfW_f 都是由 XfX_f 以及他们自身线性组合得到的，所以交集中的向量 vv 必定可表示为 Xf⋅cX_f\cdot c，所以 v∈span(Xf)v \in \text{span}(X_f)，也即 、

span(Wp)∩span(Wf)⊆span(Xf)\text{span}(W_p) \cap \text{span}(W_f) \subseteq \text{span}(X_f)

确定性子空间辨识#

N4SID#

Numerical algorithm for Subspace State Space System IDentification 子空间状态空间系统辨识的数值算法

使用消去状态变量的形式 Yf=ΓiLpWp+HiUfY_f=\Gamma_iL_pW_p+H_iU_f，做斜交投影：

Yf/ ⁣UfWp=ΓiLpWp/ ⁣UfWp+HiUf/ ⁣UfWp=ΓiLpWp=ΓiXf\begin{aligned} Y_f/\!_{U_f}W_p&= \Gamma_iL_pW_p/\!_{U_f}W_p+H_iU_f/\!_{U_f}W_p\\ &= \Gamma_iL_pW_p\\ &= \Gamma_iX_f \end{aligned}

对 O=Yf/ ⁣UfWp\mathcal{O}=Y_f/\!_{U_f}W_p 做 SVD

O=[U1U2][Σr000][V1TV2T]=U1ΣrV1T  ( =ΓiXf)\mathcal{O}=\left[ U_1 \quad U_2 \right] \begin{bmatrix} \Sigma_r & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} V_1^\mathrm{T} \\ V_2^\mathrm{T}\end{bmatrix}=U_1\Sigma_rV_1^\mathrm{T}\;(\,=\Gamma_iX_f)

为了得到 ΓiXf\Gamma_iX_f，从 Σr\Sigma_r 中间劈成两个矩阵相乘，劈开处乘上相似变换阵（非奇异即可）：

  • Γi=U1Σr1/2T\Gamma_i = U_1 \Sigma_r^{1/2}\color{red}T
  • Xf=T−1Σr1/2V1TX_f = {\color{red}T^{-1} }\Sigma_r^{1/2} V_1^\mathrm{T}{}

下面就可以进行辨识了。由于 TT、V1V_1 都满秩，得 n=r(Xf)=r(Σr)n=\mathrm{r}(X_f)=\mathrm{r}(\Sigma_r)，也即系统维数为 O\mathcal O 矩阵非零奇异值的个数。

然后就是求 ABCDABCD。我们知道 Xf=[xi, xi+1, ⋯ , xi+j−1]X_f=[x_{i},\,x_{i+1},\,\cdots,\,x_{i+j-1}]，把所有涉及这些 xx 的状态空间方程都写出来，写成矩阵形式并压缩表达式为

[Xi+1, j−1Yi, j−1]=[ABCD][Xi, j−1Ui, j−1]\begin{bmatrix} X_{i+1,\,j-1} \\ Y_{i,\,j-1} \end{bmatrix} = \begin{bmatrix} A & B \\ C & D \end{bmatrix} \begin{bmatrix} X_{i,\,j-1} \\ U_{i,\,j-1} \end{bmatrix}

其中大写的 XX、YY、UU 表示把小写打包，下角标两个数，第一个数表示打包的起点，第二个数表示打包的个数，即

Xi+1, j−1=[xi+1, ⋯ , xi+j−1]Yi, j−1=[yi, yi+1, ⋯ , yi+j−2]Xi, j−1=[xi, xi+1, ⋯ , xi+j−2]Ui, j−1=[ui, ui+1, ⋯ , ui+j−2]\begin{aligned} X_{i+1,\,j-1} &= [x_{i+1},\,\cdots,\,x_{i+j-1}]\\ Y_{i,\,j-1} &= [y_{i},\,y_{i+1},\,\cdots,\,y_{i+j-2}]\\ X_{i,\,j-1} &= [x_{i},\,x_{i+1},\,\cdots,\,x_{i+j-2}]\\ U_{i,\,j-1} &= [u_{i},\,u_{i+1},\,\cdots,\,u_{i+j-2}]\\ \end{aligned}

既然 XfX_f 都算出来了，那刚才定义的这四个东西就都是已知的，因此 ABCDABCD 可以用最小二乘法解出（注意这里是估系数矩阵，所以和传统的最小二乘在结论上有一些区别，但推导是一样的）

y=Θxy= \Theta x ∂∥y−Θx∥2∂Θ=2(y−Θx)xT=0Θ^=yxT(xxT)−1\begin{aligned} \dfrac{ {\partial}\|y-\Theta x\|^2}{ {\partial}\Theta}&= 2(y-\Theta x)x^\mathrm{T}=0\\ \hat\Theta&= yx^\mathrm{T}(xx^\mathrm{T})^{-1} \end{aligned}

==现在剩下一个没弄明白的点：为什么要在 Σr\Sigma_r 中间劈开？？==

MOESP#

Multivariable Output Error State sPace 多变量输出误差状态空间

做LQ分解

[UpYp]=[L110L21L22][Q1TQ2T]\begin{bmatrix}U_p\\Y_p\end{bmatrix}= \begin{bmatrix}L_{11}&0\\L_{21}&L_{22}\end{bmatrix}\begin{bmatrix}Q_1^{\mathrm T}\\Q_2^{\mathrm T}\end{bmatrix} ⇒ {Up=L11Q1TYp=L21Q1T+L22Q2T\Rightarrow\,\left\{ \begin{aligned} U_p&= L_{11}Q_1^{\mathrm T}\\ Y_p&= L_{21}Q_1^{\mathrm T}+L_{22}Q_2^{\mathrm T}\\ \end{aligned} \right.

由子空间方程之“过去状态拼一起”：Yp=ΓiXp+HiUpY_p=\Gamma_iX_p+H_iU_p，把 YpY_p 和 UpU_p 的 LQ 表达式代进去：

ΓiXp+HiL11Q1T=L21Q1T+L22Q2T(1)\boxed{\Gamma_iX_p+H_iL_{11}Q_1^{\mathrm T}= L_{21}Q_1^{\mathrm T}+L_{22}Q_2^{\mathrm T} }\tag{1}

两边右乘 Q2Q_2 得

ΓiXpQ2=L22\Gamma_iX_pQ_2=L_{22}

和 N4SID 类似，为了得到 ΓiXpQ2\Gamma_iX_pQ_2，只需对 L22L_{22}{} 做奇异值分解：

L22=[U1U2][Σr000][V1TV2T]=U1ΣrV1T ( =ΓiXpQ2)L_{22}=\left[ U_1 \quad U_2 \right] \begin{bmatrix} \Sigma_r & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} V_1^\mathrm{T} \\ V_2^\mathrm{T}\end{bmatrix}=U_1\Sigma_rV_1^{\mathrm T}\ (\,=\Gamma_iX_pQ_2)

然后从 Σr\Sigma_r 中间劈开得（==同样的，没懂，为啥这样劈开？？==）

  • Γi=U1Σr1/2\Gamma_i= U_1\Sigma_r^{1/2}{}
  • XpQ2=Σr1/2V1TX_pQ_2= \Sigma_r^{1/2}V_1^\mathrm{T}{}

然后就可以开始辨识了：n=r(Xp)=r(Σr)n=\mathrm{r}(X_p)={\mathrm r}(\Sigma_r)（因为 Q2Q_2 正交，即满秩）

而 Γi=[CCA⋮CAi−1]\Gamma_i=\begin{bmatrix}C\\CA\\\vdots\\CA^{i-1}\end{bmatrix}，于是可以解 AA 和 CC：

  • Γi\Gamma_i 的开头 ll 行就是矩阵 CC
  • (Γi(\Gamma_i 去掉开头 ll 行)=(Γi)=(\Gamma_i 去掉最后 ll 行)×A)\times A，解方程即得 AA

至于 BB 和 DD，我们再次回到 (1) 式，两边同时左乘 U2TU_2^\mathrm{T}{}；因为 Γi\Gamma_i 和 L22L_{22}{} 都是以 U1U_1 开头的，由 UU 正交性得，左乘 U2TU_2^\mathrm{T}{} 时整个项都变成 0，于是

U2TΓiXp+U2THiL11Q1T=U2TL21Q1T+U2TL22Q2T\begin{aligned} \bcancel{U_2^\mathrm{T}\Gamma_iX_p}+U_2^\mathrm{T}H_iL_{11}Q_1^{\mathrm T}= U_2^\mathrm{T}L_{21}Q_1^{\mathrm T}+\bcancel{U_2^\mathrm{T}L_{22}Q_2^{\mathrm T} } \end{aligned}

去掉 Q1TQ_1^\mathrm{T}{}，L11L_{11}{} 移到右边去，把 HiH_i 的具体表达式代入：

U2T[D0⋯00CBD⋯00⋮⋮⋱⋮⋮CAi−3BCAi−4B⋯D0CAi−2BCAi−3B⋯CBD]=U2TL21L11−1U_2^\mathrm{T}\begin{bmatrix} D & 0 & \cdots & 0 & 0\\ CB & D & \cdots & 0 & 0\\ \vdots & \vdots & \ddots & \vdots & \vdots \\ CA^{i-3}B & CA^{i-4}B & \cdots & D & 0\\ CA^{i-2}B & CA^{i-3}B & \cdots & CB & D \end{bmatrix}=U_2^\mathrm{T}L_{21}L_{11}^{-1}

中间这个大矩阵 HiH_i 就是待求的，其他东西都是已知的。我们把左边 U2TU_2^\mathrm{T}{} 每 ll 列分块、右边每 mm 列分块（以匹配 HiH_i 的分块模式），可以分成 ii 块

U2T=[L1, L2, ⋯ , Li]U2TL21L11−1=[M1, M2, ⋯ , Mi]\begin{aligned} U_2^\mathrm{T}&=[\mathcal{L}_1,\,\mathcal{L}_2,\,\cdots,\,\mathcal{L}_i]\\ U_2^\mathrm{T}L_{21}L_{11}^{-1}&=[\mathcal{M}_1,\,\mathcal{M}_2,\,\cdots,\,\mathcal{M}_i] \end{aligned}

分块之后，实际上是阶梯型，而且 BB 和 DD 可以彻底分开

M1=L1D+(L2C+⋯+Li−1CAi−3+LiCAi−2)BM2=L2D+(L3C+⋯+LiCAi−3)B  ⋮Mi−1=Li−1D+(LiC)BMi=LiD\begin{aligned} \mathcal{M}_1&=\mathcal{L}_1 D + (\mathcal{L}_2 C + \cdots + \mathcal{L}_{i-1} C A^{i-3} + \mathcal{L}_i C A^{i-2})B \\ \mathcal{M}_2&=\mathcal{L}_2 D + (\mathcal{L}_3 C + \cdots + \mathcal{L}_i C A^{i-3}) B \\ &\ \ \vdots \\ \mathcal{M}_{i-1}&= \mathcal{L}_{i-1} D + (\mathcal{L}_i C )B \\ \mathcal{M}_{i}&=\mathcal{L}_i D \end{aligned}

注意到括号里面很像 Γi\Gamma_i，于是记 Lkˉ=[Lk,Lk+1,⋯ ,Li]\bar{\mathcal{L}_k}=\left[\mathcal{L}_k, \mathcal{L}_{k+1}, \cdots, \mathcal{L}_i\right]，就可以写成矩阵表达式

[M1M2⋮Mi−1Mi]=[L1L2ˉΓi−1L2L3ˉΓi−2⋮⋮Li−1LiˉΓ1Li0][DB]\begin{bmatrix} \mathcal{M}_1 \\ \mathcal{M}_2 \\ \vdots \\ \mathcal{M}_{i-1} \\ \mathcal{M}_i \end{bmatrix} = \begin{bmatrix} \mathcal{L}_1 & \bar{\mathcal{L}_2}\Gamma_{i-1} \\ \mathcal{L}_2 & \bar{\mathcal{L}_3}\Gamma_{i-2} \\ \vdots & \vdots \\ \mathcal{L}_{i-1} & \bar{\mathcal{L}_i}\Gamma_1 \\ \mathcal{L}_i & 0 \end{bmatrix} \begin{bmatrix} D \\ B \end{bmatrix}

然后只需要使用最小二乘法（传统的就行了）就能解 BB 和 DD

==BD这一段的构造，逻辑还是不够丝滑==

随机子空间辨识#

问题描述#

在状态空间方程中引入一个随机偏移

{xk+1=Axk+Buk+wkyk=Cxk+Duk+vk\left\{ \begin{aligned} x_{k+1}&= Ax_k+Bu_k+w_k\\ y_k&= Cx_k+Du_k+v_k \end{aligned} \right.

其中 wkw_k 和 vkv_k 是均值为 0⃗\vec 0 的白噪向量，它们的协方差如下定义

Q=E[wkwkT]R=E[vkvkT]S=E[wkvkT]\begin{aligned} Q&= E[w_kw_k^\mathrm{T}]\\ R&= E[v_kv_k^\mathrm{T}]\\ S&= E[w_kv_k^\mathrm{T}] \end{aligned}

注意中括号里头，两个角标必须一致，角标不一致时协方差为零（因为白噪的意思就是不同时刻互相独立）。并且进一步假定 E[xkwkT]=E[xkvkT]=0E[x_kw_k^\mathrm{T}]=E[x_kv_k^\mathrm{T}]=0

考试只考没有 BB 和 DD 的情形，也即只考虑

{xk+1=Axk+wkyk=Cxk+vk\left\{ \begin{aligned} x_{k+1}&= Ax_k+w_k\\ y_k&= Cx_k+v_k \end{aligned} \right.

记号#

定义：

  1. 状态协方差 Σ=E[xkxkT]\Sigma=E[x_kx_k^\mathrm{T}]
  2. 输出协方差 Λi=E[yk+i ykT]\Lambda_i=E[y_{k+i}\,y_k^\mathrm{T}]（相差 ii 时间的）
  3. 状态输出协方差 G=E[xk+1ykT]G=E[x_{k+1}y_k^\mathrm{T}]

这些协方差应当与 kk 无关，表示系统足够稳定。于是：

Σ=E[xk+1xk+1T]=E[(Axk+wk)(Axk+wk)T]=AE[xkxkT]AT+E[wkwkT]=AΣAT+Q\begin{aligned} \Sigma&= E[x_{k+1}x_{k+1}^\mathrm{T}]\\ &= E\big[(Ax_k+w_k)(Ax_k+w_k)^\mathrm{T}\big]\\ &= AE[x_kx_k^\mathrm{T}]A^\mathrm{T}+E[w_kw_k^\mathrm{T}]\\ &= A\Sigma A^\mathrm{T}+Q \end{aligned} G=E[xk+1ykT]=E[(Axk+wk)(Cxk+vk)T]=AE[xkxkT]CT+E[wkvkT]=AΣCT+S\begin{aligned} G&= E[x_{k+1}y_{k}^\mathrm{T}]\\ &= E\big[(Ax_k+w_k)(Cx_k+v_k)^\mathrm{T}\big]\\ &= AE[x_kx_k^\mathrm{T}]C^\mathrm{T}+E[w_kv_k^\mathrm{T}]\\ &= A\Sigma C^\mathrm{T}+S \end{aligned} Λ0=E[ykykT]=E[(Cxk+vk)(Cxk+vk)T]=CE[xkxkT]CT+E[vkvkT]=CΣCT+R\begin{aligned} \Lambda_0&= E[y_{k}y_{k}^\mathrm{T}]\\ &= E\big[(Cx_k+v_k)(Cx_k+v_k)^\mathrm{T}\big]\\ &= CE[x_kx_k^\mathrm{T}]C^\mathrm{T}+E[v_kv_k^\mathrm{T}]\\ &= C\Sigma C^\mathrm{T}+R \end{aligned}

现在考虑 Λi\Lambda_i，先把前面一个 yy 用 Cx+vCx+v 转成 xx，再用 Ax+wAx+w 一直向前递归直到 xk+1x_{k+1}{}：

Λi=E[yk+iykT]=E[(Cxk+i+vk+i)ykT]=CE[xk+iykT]+0=CE[(Axk+i−1+wk+i−1)ykT]=CAE[xk+i−1ykT]+0=CAE[(Axk+i−2+wk+i−2)ykT]=CA2E[xk+i−2ykT]+0=  ⋯  =CAi−1E[xk+1ykT]+0=CAi−1G\begin{aligned} \Lambda_i&= E[y_{k+i}y_{k}^\mathrm{T}]\\ &= E\big[(Cx_{k+i}+v_{k+i})y_{k}^\mathrm{T}\big]=CE\big[x_{k+i}y_{k}^\mathrm{T}\big]+0\\ &= CE\big[(Ax_{k+i-1}+w_{k+i-1})y_{k}^\mathrm{T}\big]=CAE\big[x_{k+i-1}y_{k}^\mathrm{T}\big]+0\\ &= CAE\big[(Ax_{k+i-2}+w_{k+i-2})y_{k}^\mathrm{T}\big]=CA^2E\big[x_{k+i-2}y_{k}^\mathrm{T}\big]+0\\ &=\ \ \cdots\ \ =CA^{i-1}E\big[x_{k+1}y_{k}^\mathrm{T}\big]+0\\ &= CA^{i-1}G \end{aligned}

核心原理#

说穿了其实就是“用样本均值代替期望”。也即，只要观测的次数 jj 足够大，我就可以认为

E[X]=1j∑k=0j−1xkE[X]=\dfrac1j\sum\limits_{k=0}^{j-1}x_k

用这个思路，有如下观察：

Λi=E[yk+iykT]=1j∑k=0j−1yk+iykT=1j[yiyi+1⋯yi+j−1][y0Ty1T⋮yj−1T]=1jYi, jY0, jT\begin{aligned} \Lambda_i&= E[y_{k+i}y_k^\mathrm{T}]\\ &= \dfrac1j\sum\limits_{k=0}^{j-1}y_{k+i}y_k^\mathrm{T}\\ &= \dfrac1j\begin{bmatrix}y_i&y_{i+1}&\cdots&y_{i+j-1}\end{bmatrix}\begin{bmatrix}y_0^\mathrm{T}\\y_{1}^\mathrm{T}\\\vdots\\y_{j-1}^\mathrm{T}\end{bmatrix}\\ &= \dfrac1jY_{i,\,j}Y_{0,\,j}^\mathrm{T} \end{aligned}

其中大写 YY 的含义和 N4SID 中提到的是一样的，下角标第一个数表示打包的起点，第二个数表示打包的个数

此时我定义：

Φ[A, B]=1jABT\Phi_{[A,\,B]}=\dfrac1jAB^\mathrm{T}

那么

Λi=Φ[Yi, j, Y0, j]\Lambda_i=\Phi_{[Y_{i,\,j},\,Y_{0,\,j}]}
控制科学笔记 - 子空间辨识
https://dingnuooo.top/blog/control-science/subspace-identification
Author Dingnuooo
Published at June 11, 2025
Copyright CC BY-NC-SA 4.0
Buy me a cup of coffee ☕.

控制科学笔记 - 参数估计

控制科学笔记 - 空间投影
Comment seems to stuck. Try to refresh?✨