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}