8. The ERA with Data Correlations (ERA/DC)
The Eigensystem Realization Algorithm with Data Correlations (ERA/DC) includes an additional fit to output correlations whereas the ERA is basically a least-square fit to the pulse response measurements only. The bias terms affecting the ERA when noise is present can, in principle, be omitted in the ERA/DC by properly tuning some of the parameters. The computational steps of the ERA/DC are outlined in this section.
8-1) Block Correlation Hankel Matrices
The ERA method with Data Correlation (ERA/DC) requires the definition of a square matrix of order \(\gamma = mp\),
\begin{align} \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k) &= \boldsymbol{H}(k)\boldsymbol{H}(0)^\intercal\\ &= \begin{bmatrix} h_{k+1} & h_{k+2} & \cdots & h_{k+q}\\ h_{k+2} & h_{k+3} & \cdots & h_{k+q+1}\\ \vdots & \vdots & \ddots & \vdots\\ h_{k+p} & h_{k+p+1} & \cdots & h_{k+p+q-1}\\ \end{bmatrix}\begin{bmatrix} h_{1} & h_{2} & \cdots & h_{q}\\ h_{2} & h_{3} & \cdots & h_{q+1}\\ \vdots & \vdots & \ddots & \vdots\\ h_{p} & h_{p+1} & \cdots & h_{p+q-1}\\ \end{bmatrix}^\intercal\\ &= \begin{bmatrix} \displaystyle\sum_{i=1}^qh_{k+i}h_i^\intercal & \displaystyle\sum_{i=1}^qh_{k+i}h_{i+1}^\intercal & \cdots & \displaystyle\sum_{i=1}^qh_{k+i}h_{p+i-1}^\intercal\\ \displaystyle\sum_{i=1}^qh_{k+i+1}h_i^\intercal & \displaystyle\sum_{i=1}^qh_{k+i+1}h_{i+1}^\intercal & \cdots & \displaystyle\sum_{i=1}^qh_{k+i+1}h_{p+i-1}^\intercal\\ \vdots & \vdots & \ddots & \vdots \\ \displaystyle\sum_{i=1}^qh_{k+p+i-1}h_i^\intercal & \displaystyle\sum_{i=1}^qh_{k+p+i-1}h_{i+1}^\intercal & \cdots & \displaystyle\sum_{i=1}^qh_{k+p+i-1}h_{p+i-1}^\intercal \end{bmatrix} \end{align}
Note that the data correlation matrix \(\boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k)\) can be smaller in size than the Hankel matrix \(\boldsymbol{H}(k)\) if \(qr\leq pm\).
For the case when \(k=0\), the correlation matrix \(\boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(0)\) becomes
\begin{align} \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(0) &= \boldsymbol{H}(0)\boldsymbol{H}(0)^\intercal\\ &= \begin{bmatrix} h_{1} & h_{2} & \cdots & h_{q}\\ h_{2} & h_{3} & \cdots & h_{q+1}\\ \vdots & \vdots & \ddots & \vdots\\ h_{p} & h_{p+1} & \cdots & h_{p+q-1}\\ \end{bmatrix}\begin{bmatrix} h_{1} & h_{2} & \cdots & h_{q}\\ h_{2} & h_{3} & \cdots & h_{q+1}\\ \vdots & \vdots & \ddots & \vdots\\ h_{p} & h_{p+1} & \cdots & h_{p+q-1}\\ \end{bmatrix}^\intercal\\ &= \begin{bmatrix} \displaystyle\sum_{i=1}^qh_{i}h_i^\intercal & \displaystyle\sum_{i=1}^qh_{i}h_{i+1}^\intercal & \cdots & \displaystyle\sum_{i=1}^qh_{i}h_{p+i-1}^\intercal\\ \displaystyle\sum_{i=1}^qh_{i+1}h_i^\intercal & \displaystyle\sum_{i=1}^qh_{i+1}h_{i+1}^\intercal & \cdots & \displaystyle\sum_{i=1}^qh_{i+1}h_{p+i-1}^\intercal\\ \vdots & \vdots & \ddots & \vdots \\ \displaystyle\sum_{i=1}^qh_{p+i-1}h_i^\intercal & \displaystyle\sum_{i=1}^qh_{p+i-1}h_{i+1}^\intercal & \cdots & \displaystyle\sum_{i=1}^qh_{p+i-1}h_{p+i-1}^\intercal \end{bmatrix} \end{align}
The matrix \(\boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(0)\) consists of auto-correlations of Markov parameters such as \(\displaystyle\sum_{i=1}^qh_{i}h_i^\intercal\) and cross-correlations between outputs such as \(\displaystyle\sum_{i=1}^qh_{i}h_{i+1}^\intercal\) at lag time values in the range \(\pm p\), summed over \(q\) values. If noises in the Markov parameters are not correlated, the correlation matrix \(\boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(0)\) will contain less noise than the Hankel matrix \(\boldsymbol{H}(0)\).
In terms of controllability and observability matrices, \(\boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k)\) can be written as
\begin{align} \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k) = \boldsymbol{O}_pA^{k}\boldsymbol{R}_q\boldsymbol{R}_q^\intercal \boldsymbol{O}_p^\intercal = \boldsymbol{O}_pA^{k}\boldsymbol{R}_\gamma, \end{align}
where \(\boldsymbol{R}_\gamma = \boldsymbol{R}_q\boldsymbol{R}_q^\intercal \boldsymbol{O}_p^\intercal\).
The data correlation matrix \(\boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k)\) can be used to form a block correlation Hankel matrix
\begin{align} \label{Eq: Block Correlation Hankel Matrix k} \boldsymbol{\mathcal{H}}(k) &= \begin{bmatrix} \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k) & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+\tau) & \cdots & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+\zeta\tau)\\ \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+\tau) & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+2\tau) & \cdots & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+(\zeta+1)\tau)\\ \vdots & \vdots & \ddots & \vdots\\ \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+\xi\tau) & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+(\xi+1)\tau) & \cdots & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(k+(\xi+\zeta)\tau) \end{bmatrix}\\ &=\begin{bmatrix} \boldsymbol{O}_p\\ \boldsymbol{O}_pA^\tau\\ \vdots\\ \boldsymbol{O}_pA^{\xi\tau} \end{bmatrix}A^k\begin{bmatrix} \boldsymbol{R}_\gamma & A^\tau\boldsymbol{R}_\gamma & \cdots & A^{\zeta\tau}\boldsymbol{R}_\gamma \end{bmatrix}\\ &=\boldsymbol{O}_\xi A^k\boldsymbol{R}_\zeta. \end{align}
For the case when \(k=0\),
\begin{align} \label{Eq: Block Correlation Hankel Matrix 0} \boldsymbol{\mathcal{H}}(0) &= \begin{bmatrix} \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(0) & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(\tau) & \cdots & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(\zeta\tau)\\ \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(\tau) & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(2\tau) & \cdots & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}((\zeta+1)\tau)\\ \vdots & \vdots & \ddots & \vdots\\ \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}(\xi\tau) & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}((\xi+1)\tau) & \cdots & \boldsymbol{\mathcal{R}}_{\boldsymbol{H}\boldsymbol{H}}((\xi+\zeta)\tau) \end{bmatrix}\\ &=\begin{bmatrix} \boldsymbol{O}_p\\ \boldsymbol{O}_pA^\tau\\ \vdots\\ \boldsymbol{O}_pA^{\xi\tau} \end{bmatrix}\begin{bmatrix} \boldsymbol{R}_\gamma & A^\tau\boldsymbol{R}_\gamma & \cdots & A^{\zeta\tau}\boldsymbol{R}_\gamma \end{bmatrix}\\ &=\boldsymbol{O}_\xi\boldsymbol{R}_\zeta. \end{align}
\(\tau\) is an integer chosen to chosen to prevent significant overlap of adjacent correlation blocks. The matrices \(\boldsymbol{R}_\zeta\) and \(\boldsymbol{O}_\xi\) are called the block correlation controllability and observability matrices.
8-2) Hankel Norm Approximation
Similarly to the ERA, the ERA/DC process continues with the factorization of the block correlation matrix \(\boldsymbol{\mathcal{H}}(0)\) using singular value decomposition so that
\begin{align} \boldsymbol{\mathcal{H}}(0) = \boldsymbol{\mathcal{U}}\boldsymbol{\Sigma}\boldsymbol{\mathcal{V}}^\intercal = \begin{bmatrix} \boldsymbol{\mathcal{U}}_n & \boldsymbol{\mathcal{U}}_{0} \end{bmatrix}\begin{bmatrix} \boldsymbol{\Sigma}_n & \boldsymbol{0}\\ \boldsymbol{0} & \boldsymbol{\Sigma}_{0} \end{bmatrix}\begin{bmatrix} \boldsymbol{\mathcal{V}}_n^\intercal\\ \boldsymbol{\mathcal{V}}_{0}^\intercal \end{bmatrix} = \boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n\boldsymbol{\mathcal{V}}_n^\intercal+\underbrace{\boldsymbol{\mathcal{U}}_{0}\boldsymbol{\Sigma}_{0}\boldsymbol{\mathcal{V}}_{0}^\intercal}_{\simeq\boldsymbol{0}} \simeq \boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n\boldsymbol{\mathcal{V}}_n^\intercal, \end{align}
and
\begin{align} \boldsymbol{\mathcal{H}}(0) = \boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n\boldsymbol{\mathcal{V}}_n^\intercal = \boldsymbol{O}_\xi\boldsymbol{R}_\zeta \Rightarrow \left\lbrace\begin{array}{ll} \boldsymbol{O}_\xi &\hspace{-0.7em}= \boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n^{1/2}\\ \boldsymbol{R}_\zeta &\hspace{-0.7em}= \boldsymbol{\Sigma}_n^{1/2}\boldsymbol{\mathcal{V}}_n^\intercal \end{array}\right.. \end{align}
Again, this choice makes both \(\boldsymbol{O}_\xi\) and \(\boldsymbol{R}_\zeta\) balanced.
8-3) Minimum Realization
From Eq. (\ref{Eq: Block Correlation Hankel Matrix 0}) we have directly
\begin{align} \boldsymbol{O}_p = \boldsymbol{E}_\gamma^\intercal\boldsymbol{O}_\xi = \boldsymbol{E}_\gamma^\intercal\boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n^{1/2}. \end{align}
From the classical definition of the controllability matrix, an expression of \(\boldsymbol{R}_q\) can be found
\begin{align} \boldsymbol{R}_q = \boldsymbol{O}_p^\dagger\boldsymbol{H}(0) = \left[\boldsymbol{E}_\gamma^\intercal\boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n^{1/2}\right]^\dagger\boldsymbol{H}(0), \end{align}
and a realization is shown to be
\begin{align} \widehat{A} &= \boldsymbol{O}_\xi^\dagger\boldsymbol{\mathcal{H}}(1)\boldsymbol{R}_\zeta^\dagger = \boldsymbol{\Sigma}_n^{-1/2}\boldsymbol{\mathcal{U}}_n^\intercal\boldsymbol{\mathcal{H}}(1)\boldsymbol{\mathcal{V}}_n\boldsymbol{\Sigma}_n^{-1/2},\\ \widehat{B} &= \boldsymbol{R}_q\boldsymbol{E}_r = \left[\boldsymbol{E}_\gamma^\intercal\boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n^{1/2}\right]^\dagger\boldsymbol{H}(0)\boldsymbol{E}_r,\\ \widehat{C} &= \boldsymbol{E}_m^\intercal\boldsymbol{O}_p = \boldsymbol{E}_m^\intercal\boldsymbol{E}_\gamma^\intercal\boldsymbol{\mathcal{U}}_n\boldsymbol{\Sigma}_n^{1/2},\\ \widehat{D} &= h_0. \end{align}