Eigensystem implementation algorithm

The Eigensystem Realization Algorithm , if translated literally from the English Eigenensystem Realization Algorithm , ERA for short , is a system identification method for determining a linear, time-discrete system based on input-output measurement data . The algorithm is suitable when no analytical model can be formulated on the basis of physical relationships. The ERA was first introduced in 1985 by JN Juang and RS Pappa.

Theoretical background

A general linear, time-discrete system

${\displaystyle {\begin{aligned}{\underline {x}}_{k+1}&=A{\underline {x}}_{k}+B{\underline {u}}_{k}\\{\underline {y}}_{k}&=C{\underline {x}}_{k}+D{\underline {u}}_{k}\end{aligned}}}$

with the states , the input and the output , indicates a pulse-shaped excitation with the initial condition , the following evolution in time ${\displaystyle {\underline {x}}\in \mathbb {R} ^{n}}$ ${\displaystyle {\underline {u}}\in \mathbb {R} ^{p}}$ ${\displaystyle {\underline {y}}\in \mathbb {R} ^{q}}$ ${\displaystyle {\underline {x}}_{0}={\underline {0}}}$

${\displaystyle {\begin{array}{l l l l}&{\underline {u}}&{\underline {x}}&{\underline {y}}\\k=0&I&{\underline {0}}&D\\k=1&0&B&CB\\k=2&0&AB&CAB\\k=3&0&A^{2}B&CA^{2}B\\\vdots &\vdots &\vdots &\vdots \\k=i&0&A^{(i-1)}B&CA^{(i-1)}B\\\vdots &\vdots &\vdots &\vdots \\k=(N-1)&0&A^{(N-2)}B&CA^{(N-2)}B\\\end{array}}}$

where the Dirac impulse symbolizes in the respective input, means the respective time step and represents the number of samples. In the impulse response of the information inherent on the system matrices , , and contain what are the only data that are present after the measurement to the system. With the ERA, these matrices can be extracted from the measurement data. The linear and time-discrete model implemented in this way is similar to the method of balanced model reduction ( balanced truncation ). This means that the states of the implemented system are the best and equally observable and controllable states of the identified system. Ultimately, the realized model is a balanced model reduction of the physical system on which the measurements were made. At this point it should be pointed out that it is very difficult, if not impossible, to physically generate the Dirac impulse because it represents a mathematical ideal function and, at best, can be approximated. However, there are options to filter out the impulse response from random input / output measurement data. This can be achieved, for example, with the Observer Kalman Filter Identification . In addition, the realized states are usually not physically interpretable. However, there is the possibility, in the case of a complete state measurement of the system, to form the controllability matrix and thus to draw conclusions about a subspace transformation with which the realized states can finally be interpreted. This is related to the order reduction method Balanced Proper Orthogonal Decomposition . ${\displaystyle I}$ ${\displaystyle k}$ ${\displaystyle N}$ ${\displaystyle {\underline {y}}}$ ${\displaystyle A}$ ${\displaystyle B}$ ${\displaystyle C}$ ${\displaystyle D}$

The algorithm

The first step is obviously to collect measurement data. In the following, a MIMO system is assumed in order not to have to make a distinction to the SISO case. A SISO system represents the special case that and . For the measurement, each input is excited in a pulse-like manner and each output is measured. From this follow output vectors each with samples in time. Exemplary obtaining the measurement data from input to output to ${\displaystyle p=1}$ ${\displaystyle q=1}$ ${\displaystyle u_{i}}$ ${\displaystyle (i=1,2,...,p)}$ ${\displaystyle y_{j}}$ ${\displaystyle (j=1,2,...,q)}$ ${\displaystyle (p\times q)}$ ${\displaystyle {\underline {y}}^{ij}}$ ${\displaystyle N}$ ${\displaystyle y_{k}^{ij}}$ ${\displaystyle u_{2}}$ ${\displaystyle y_{3}}$ ${\displaystyle {\underline {y}}^{23}={\begin{bmatrix}y_{0}^{23}&y_{1}^{23}&\dots &y_{N-1}^{23}\end{bmatrix}}.}$

The second step is to arrange the measurement data in an array . This array corresponds to the evolution of the output introduced at the beginning . The blocks are called Markov parameters. It immediately follows that . ${\displaystyle {\underline {\textbf {Y}}}\in \mathbb {R} ^{q\times p\times N}}$ ${\displaystyle {\underline {y}}}$ ${\displaystyle N}$ ${\displaystyle {\tilde {y}}_{k}\in \mathbb {R} ^{q\times p}}$ ${\displaystyle D={\tilde {y}}_{0}}$

The third step is to stack the Markov parameters in a block Hankel matrix ${\displaystyle {\textbf {H}}\in \mathbb {R} ^{[bq]\times [sp]}}$

${\displaystyle {\textbf {H}}(k=1)={\textbf {H}}={\begin{bmatrix}{\tilde {y}}_{1}&{\tilde {y}}_{2}&\cdots &{\tilde {y}}_{s}\\{\tilde {y}}_{2}&{\tilde {y}}_{3}&\cdots &{\tilde {y}}_{s+1}\\\vdots &\vdots &&\vdots \\{\tilde {y}}_{b}&{\tilde {y}}_{b+1}&\cdots &{\tilde {y}}_{b+s-1}\\\end{bmatrix}}}$

and a shifted block Hankel matrix ${\displaystyle {\textbf {H}}'\in \mathbb {R} ^{[bq]\times [sp]}}$

${\displaystyle {\textbf {H}}(k=2)={\textbf {H}}'={\begin{bmatrix}{\tilde {y}}_{2}&{\tilde {y}}_{3}&\cdots &{\tilde {y}}_{s+1}\\{\tilde {y}}_{3}&{\tilde {y}}_{4}&\cdots &{\tilde {y}}_{s+2}\\\vdots &\vdots &&\vdots \\{\tilde {y}}_{b+1}&{\tilde {y}}_{b+2}&\cdots &{\tilde {y}}_{b+s}\\\end{bmatrix}}}$.

The meaning of the parameters will be clarified below. The two block Hankel matrices can now be expressed by the observability matrix and the controllability matrix . Finally, it applies to the block Hankel matrix and the shifted block Hankel matrix ${\displaystyle s,b\in \mathbb {N} \backslash \lbrace 0\rbrace }$ ${\displaystyle Q_{B}}$ ${\displaystyle Q_{S}}$

${\displaystyle Q_{B}Q_{S}={\begin{bmatrix}C\\CA\\\vdots \\CA^{b-1}\end{bmatrix}}{\begin{bmatrix}B&AB&\cdots &A^{s-1}B\end{bmatrix}}={\begin{bmatrix}CB&CAB&\cdots &CA^{s-1}B\\CAB&CA^{2}B&\cdots &CA^{s}B\\\vdots &\vdots &&\vdots \\CA^{b-1}B&CA^{b}B&\cdots &CA^{s+b-2}B\end{bmatrix}}={\textbf {H}}}$

and respectively. ${\displaystyle Q_{B}AQ_{S}={\textbf {H}}'}$

At this point it is clear that and respectively set the size of Steuerbarkeitsmatrix and the observability and thus the maximum achievable system order limit of the system to be implemented. They can be chosen freely and are only limited by the number of samples, as must apply . One possible choice is ${\displaystyle s}$ ${\displaystyle b}$ ${\displaystyle r}$ ${\displaystyle N}$ ${\displaystyle b+s=N-1}$

${\displaystyle s=b=\left\lfloor {\frac {N-1}{2}}\right\rfloor }$,

which corresponds to the maximum size for and and makes the block Hankel matrices square. In other words, the size of the Hankel matrices and thus ultimately the quality of the implementation limit. ${\displaystyle s}$ ${\displaystyle b}$ ${\displaystyle N}$

The fourth step is the (economic) singular value decomposition of ${\displaystyle {\textbf {H}}}$

${\displaystyle {\textbf {H}}=U\Sigma V^{\top }={\begin{bmatrix}{\tilde {U}}&U_{t}\end{bmatrix}}{\begin{bmatrix}{\tilde {\Sigma }}&{\underline {0}}\\{\underline {0}}&\Sigma _{t}\end{bmatrix}}{\begin{bmatrix}{\tilde {V}}^{\top }\\V_{t}^{\top }\end{bmatrix}}}$,

with , and , where tilde denotes the part that is used for the implementation and the part that is truncated . At this point, a decision about the system order of the system to be implemented must be made. The decision is based on the Hankelsingulärwerte that in hit are included. The Hankelsingular values ​​are a measure for the transported energy from input to output of the respective state. As already mentioned, the states are balanced. As a result, the largest Hankel singular values ​​are among the states that can best be observed and controlled to the same extent. As a rule, a prominent position can be found from which the singular values ​​become significantly smaller. ${\displaystyle {\tilde {U}}\in \mathbb {R} ^{[bq]\times r}}$ ${\displaystyle {\tilde {\Sigma }}\in \mathbb {R} ^{r\times r}}$ ${\displaystyle {\tilde {V}}^{\top }\in \mathbb {R} ^{r\times [sp]}}$ ${\displaystyle t}$ ${\displaystyle r}$ ${\displaystyle \sigma }$ ${\displaystyle \Sigma }$ ${\displaystyle r}$ ${\displaystyle \sigma _{i}}$

${\displaystyle \sigma _{1}>\sigma _{2}>\dots >\sigma _{r}>>\sigma _{r+1}>\dots }$

If this is not the case, the cumulative energy of the singular values ​​can be determined and the number of states required for this can be determined by specifying the minimum energy taken into account . ${\displaystyle r}$

The fifth step is the determination of the implemented system. If you look first

${\displaystyle {\begin{aligned}{\textbf {H}}&=U\Sigma V^{\top }\\&=U\Sigma ^{1/2}\Sigma ^{1/2}V^{\top }\\&=Q_{B}Q_{S}\end{aligned}}}$

follows

${\displaystyle {\begin{aligned}Q_{B}&=U\Sigma ^{1/2}\\Q_{S}&=\Sigma ^{1/2}V^{\top }.\end{aligned}}}$

With the context

${\displaystyle {\textbf {H}}'=Q_{B}AQ_{S}}$

finally follows for the system matrix ${\displaystyle A}$

${\displaystyle {\begin{aligned}A&=Q_{B}^{-1}{\textbf {H}}'Q_{S}^{-1}\\&=\Sigma ^{-1/2}U^{\top }{\textbf {H}}'V\Sigma ^{-1/2}.\end{aligned}}}$

Note that for the unitary matrices and we have and . ${\displaystyle U}$ ${\displaystyle V}$ ${\displaystyle U^{-1}=U^{\top }}$ ${\displaystyle V^{-1}=V^{\top }}$

In addition, the connections follow ${\displaystyle {\textbf {H}}=Q_{B}Q_{S}}$

${\displaystyle {\begin{aligned}Q_{B}&={\textbf {H}}Q_{S}^{-1}={\textbf {H}}V\Sigma ^{-1/2}={\begin{bmatrix}\mathbf {H} _{q,sb}\\\mathbf {H} _{(b-1)q,sb}\end{bmatrix}}{\begin{bmatrix}{\tilde {V}}&V_{t}\end{bmatrix}}{\begin{bmatrix}{\tilde {\Sigma }}^{-1/2}&\\&\Sigma _{t}^{-1/2}\end{bmatrix}}={\begin{bmatrix}C\\CA\\\vdots \\CA^{b-1}\end{bmatrix}}\\Q_{S}&=Q_{B}^{-1}{\textbf {H}}=\Sigma ^{-1/2}U^{\top }{\textbf {H}}={\begin{bmatrix}{\tilde {\Sigma }}^{-1/2}&\\&\Sigma _{t}^{-1/2}\end{bmatrix}}{\begin{bmatrix}{\tilde {U}}^{\top }\\U_{t}^{\top }\end{bmatrix}}{\begin{bmatrix}\mathbf {H} _{bq,p}&\mathbf {H} _{bq,(s-1)b}\end{bmatrix}}={\begin{bmatrix}B&AB&\dots &A^{s-1}B\end{bmatrix}},\end{aligned}}}$

which contain the output matrix and the input matrix. At this point the Hankel matrix was partitioned accordingly in order to come to the matrices and later . Finally, the implemented system follows, ${\displaystyle C}$ ${\displaystyle B}$ ${\displaystyle B}$ ${\displaystyle C}$

${\displaystyle {\begin{aligned}{\underline {\tilde {x}}}_{k+1}&={\tilde {A}}{\underline {\tilde {x}}}_{k}+{\tilde {B}}{\underline {u}}_{k}\\{\underline {\hat {y}}}_{k}&={\tilde {C}}{\underline {\tilde {x}}}_{k}+{\tilde {D}}{\underline {u}}_{k}\end{aligned}}}$

with the states , the input and the approximated output , as well as the system matrices ${\displaystyle {\underline {\tilde {x}}}\in \mathbb {R} ^{r}}$ ${\displaystyle {\underline {u}}\in \mathbb {R} ^{p}}$ ${\displaystyle {\underline {\hat {y}}}\in \mathbb {R} ^{q}}$

${\displaystyle {\begin{aligned}{\tilde {A}}&={\tilde {\Sigma }}^{-1/2}{\tilde {U}}^{\top }{\textbf {H}}'{\tilde {V}}{\tilde {\Sigma }}^{-1/2}\\{\tilde {B}}&={\tilde {\Sigma }}^{-1/2}{\tilde {U}}^{\top }{\textbf {H}}_{bq,p}\\{\tilde {C}}&={\textbf {H}}_{q,sp}{\tilde {V}}{\tilde {\Sigma }}^{-1/2}\\{\tilde {D}}&={\tilde {y}}_{0},\end{aligned}}}$

with , , and . ${\displaystyle {\tilde {A}}\in \mathbb {R} ^{r\times r}}$ ${\displaystyle {\tilde {B}}\in \mathbb {R} ^{r\times p}}$ ${\displaystyle {\tilde {C}}\in \mathbb {R} ^{q\times r}}$ ${\displaystyle {\tilde {D}}\in \mathbb {R} ^{q\times p}}$

Individual evidence

1. ↑ JN Juang, RS Pappa: An Eigenensystem Realization Algorithm (ERA) for modal parameter identification and model reduction . Ed .: NASA Langley Research Center. Hampton, VA 1985 ( aiaa.org ).
2. ^ Otto Föllinger: Control engineering . 12th edition. VDE Verlag, Berlin 2016, ISBN 978-3-8007-4201-1 .