System Identification
System identification procedures for LTI systems generally accept a collection of set of input data, and output data, and finds four arrays. These output arrays are typically denoted by \(A\), \(B\), \(C\), and \(D\) in the literature and characterize a state-space system according to the following equations:
\[ \begin{gathered} \boldsymbol{x}_{k+1}=A \boldsymbol{x}_k + B \boldsymbol{u}_k \\ \boldsymbol{y}_k = C \boldsymbol{x}_k + D \boldsymbol{u}_k \end{gathered} \]
where \(\boldsymbol{x}_k \in \mathbb{R}^n\) is an internal state vector at time index \(k\), \(\boldsymbol{u}_k \in \mathbb{R}^r\) is an input vector corresponding to \(r\) inputs, and \(\boldsymbol{y}_k \in \mathbb{R}^m\) is an output vector associated with \(m\) sensor measurements.
SRIM
The following outline of the SRIM procedure is reproduced from @juang1997system :
Choose an integer \(p\) such that \(p \geq(n / m)+1\), where \(n\) is the desired order of the system and \(m\) is the number of outputs.
Compute correlation matrices \(\mathcal{R}_{y y}\) of dimension \(p m \times p m\), \(\mathcal{R}_{y u}\) of dimension pm \(\times p r\), and \(\mathcal{R}_{uu}\) of dimension pr \(\times p r\) as defined in Eq. (13) using the matrices \(Y_{p}(k)\) of dimension \(p m \times N\) and \(U_{p}(k)\) of dimension \(p r \times N\) defined in Eq. (12). The integer \(r\) is the number of inputs. The index \(k\) is the data point used as the starting point for system identification. The integer \(N\) must be chosen such that \(\ell-k-p+2 \geq N \gg \min (p m, p r)\), where \(\ell\) is the length of the data.
Calculate the correlation matrix \(\mathcal{R}_{h h}\) of dimension \(p m \times p m\) defined in Eq. (13), i.e., \(\mathcal{R}_{h h}=\mathcal{R}_{y y}-\mathcal{R}_{y u} \mathcal{R}_{u u}^{-1} \mathcal{R}_{y u}^{T} .\)
Use singular value decomposition to factor \(\mathcal{R}_{h h}\) for the full decomposition method [Eq. (25)] or a portion of \(\mathcal{R}_{h h}\) for the partial decomposition method [Eq. (30)].
Determine the order \(n\) of the system by examining the singular values of \(\mathcal{R}_{h h}\) and obtain \(\mathcal{U}_{n}\) of dimension \(p m \times n\left[\mathrm{Eq} .\right.\) (25) and \(\mathcal{U}_{0}\) of dimension \(p m \times n_{0}\), where \(n_{0}=p m-n\) is the number of truncated small singularvalues. The integer \(n_{0}\) must satisfy the condition \(p n_{0} \geq(m+n)\) for the full decomposition method. For the partial decomposition method, \(\mathcal{U}_{n}\) is replaced by \(\mathcal{U}_{n}^{\prime}[\mathrm{Eq} .(30)]\), and the integer \(n_{0}\) is the sum of \(m\) and the number of truncated singular values.
Let \(\mathcal{U}_{n}=\mathcal{O}_{p}\) or \(\mathcal{U}_{n}^{\prime}=\mathcal{O}_{p} .\) Use Eq. (7) to determine the state matrix \(A\). The output matrix \(C\) is the first \(m\) rows of \(\mathcal{U}_{n}\).
Compute \(\mathcal{U}_0 \mathcal{R}=\mathcal{U}_0^T \mathcal{R}_{yu} \mathcal{R}_{uu}^{-1}\) shown in Eq. (38) for the indirect method and construct \(\mathcal{U}_{0 n}\) and \(\mathcal{U}_{0} \mathcal{T}\) shown in Eqs. (37) and (39). Determine the input matrix \(B\) and the direct transmission matrix \(D\) from Eq. (40), i.e., the first \(m\) rows of \(\mathcal{U}_{0 n}^{\dagger} \mathcal{U}_{0} \mathcal{T}\) form \(D\) and the last \(n\) rows produce \(B\). For the direct method, construct \(\mathcal{O}_{p_{\Gamma}}\) and \(\mathcal{O}_{p A}\) from Eq. (52) and solve for \(B\) and \(D\) by computing \(\mathcal{O}_{p A}^{\dagger} \mathcal{O}_{p \Gamma}\). The first \(m\) rows of \(\mathcal{O}_{p A}^{\dagger} \mathcal{O}_{p \Gamma}\) form \(D\), and the last \(n\) rows produce \(B\). For the output-error minimization method, construct \(\boldsymbol{y}_{N}(0)\) and \(\Phi\) from Eq. (63) and solve for \(B\) and \(D\) by computing \(\Phi^{\dagger} \boldsymbol{y}_{N}(0)\). The first \(n\) elements of \(\Phi^{\dagger} \boldsymbol{y}_{N}(0)\) form the initial state vector \(\boldsymbol{x}(0)\), the
Notes on the efficient implementation of the SRIM algorithm
- Precompute and store \(CA^n\)
- Precompute and store \(\mathcal{U}\)
- Parallelize assembly of \(\Phi\)