Eigensystem Realisierungsalgorithmus

Der Eigensystem Realisierungsalgorithmus, wenn man ihn wörtlich aus dem Englischen von Eigensystem Realization Algorithm, kurz ERA, übersetzt, ist ein Verfahren der Systemidentifikation zur Bestimmung eines linearen, zeitdiskreten Systems auf Basis von Eingangs-Ausgangsmessdaten. Der Algorithmus bietet sich dann an, wenn kein analytisches Modell auf Grundlagen physikalischer Zusammenhänge formuliert werden kann. Der ERA wurde erstmals 1985 von J. N. Juang und R. S. Pappa vorgestellt.

Theoretischer Hintergrund

Ein allgemeines lineares, zeitdiskretes System

$\begin{array}{cc}{\underset{\_}{x}}_{k+1}& =A{\underset{\_}{x}}_k+B{\underset{\_}{u}}_k\\ {\underset{\_}{y}}_k& =C{\underset{\_}{x}}_k+D{\underset{\_}{u}}_k\end{array}$

mit den Zuständen $\underset{\_}{x}\in {ℝ}^n$, dem Eingang $\underset{\_}{u}\in {ℝ}^p$ und dem Ausgang $\underset{\_}{y}\in {ℝ}^q$, zeigt auf eine impulsförmige Anregung mit der Anfangsbedingung ${\underset{\_}{x}}_0=\underset{\_}{0}$, die folgende zeitliche Evolution

$\begin{array}{cccc}& \underset{\_}{u}& \underset{\_}{x}& \underset{\_}{y}\\ k=0& I& \underset{\_}{0}& D\\ k=1& 0& B& CB\\ k=2& 0& AB& CAB\\ k=3& 0& A^{2}B& C{A}^{2}B\\ ⋮& ⋮& ⋮& ⋮\\ k=i& 0& A^{(i-1)}B& C{A}^{(i-1)}B\\ ⋮& ⋮& ⋮& ⋮\\ k=(N-1)& 0& A^{(N-2)}B& C{A}^{(N-2)}B\end{array}$

wobei $I$ den Dirac-Impuls im jeweiligen Eingang symbolisiert, $k$ den jeweiligen Zeitschritt bedeutet und $N$ die Anzahl der Samples darstellt. In der Impulsantwort $\underset{\_}{y}$ sind inhärent die Informationen über die Systemmatrizen $A$, $B$, $C$ und $D$ enthalten, was die einzigen Daten sind, die nach der Messung am System vorliegen. Mit dem ERA können diese Matrizen aus den Messdaten extrahiert werden. Das damit realisierte lineare und zeitdiskrete Modell ist ähnlich der Methode der balancierten Modellreduktion (Balanced Truncation) balanciert. Das bedeutet, dass die Zustände des realisierten Systems die am besten und gleichermaßen beobachtbaren und steuerbaren Zustände des identifizierten Systems sind. Letztlich ist also das realisierte Modell eine balancierte Modellreduktion des physikalischen Systems an dem die Messungen gemacht wurden. An dieser Stelle sei noch darauf hingewiesen, dass es nur sehr schwer, bis unmöglich ist, den Dirac-Impuls physikalisch zu erzeugen, weil er eine mathematische Idealfunktion darstellt und bestenfalls approximiert werden kann. Es gibt jedoch Möglichkeiten aus zufälligen Ein-Ausgangsmessdaten die Impulsantwort herauszufiltern. Dies kann beispielsweise mit der Observer Kalman Filter Identification erreicht werden. Außerdem sind die realisierten Zustände in der Regel nicht physikalisch interpretierbar. Es besteht jedoch die Möglichkeiten, bei einer vollständigen Zustandsmessung des Systems, die Steuerbarkeitsmatrix zu bilden und damit auf eine Unterraumtransformation rückzuschließen, mit der schließlich die realisierten Zustände interpretiert werden können. Dies ist verwandt mit der Ordnungsreduktionsmethode Balanced Proper Orthogonal Decomposition.

Der Algorithmus

Der erste Schritt ist offensichtlich das Sammeln von Messdaten. Im Folgenden wird von einem MIMO-System ausgegangen, um keine Unterscheidung zum SISO-Fall machen zu müssen. Ein SISO-System stellt insofern den Sonderfall dar, dass $p=1$ und $q=1$. Für die Messung wird jeder Eingang $u_i$ $(i=1,2,...,p)$ impulsartig angeregt und jeder Ausgang $y_j$ $(j=1,2,...,q)$ gemessen. Daraus folgen $(p\times q)$ Ausgangsvektoren ${\underset{\_}{y}}^{ij}$ mit jeweils $N$ Samples $y_k^{ij}$ in der Zeit. Exemplarische erhält man die Messdaten von Eingang $u_2$ zum Ausgang $y_3$ zu ${\underset{\_}{y}}^{23}=\left[\begin{array}{cccc}{y}_0^{23}& y_1^{23}& \dots & y_{N-1}^{23}\end{array}\right].$

Der zweite Schritt ist das Arrangieren der Messdaten in einem Array $\underset{\_}{\mathbf{\text{Y}}}\in {ℝ}^{q\times p\times N}$. Dieses Array entspricht der eingangs eingeführten Evolution des Ausgangs $\underset{\_}{y}$. Die $N$ Blöcke ${\tilde{y}}_k\in {ℝ}^{q\times p}$werden als Markov Parameter bezeichnet. Es folgt unmittelbar, dass $D={\tilde{y}}_0$.

Der dritte Schritt ist das Stapeln der Markov Parameter in eine Block Hankelmatrix $\mathbf{\text{H}}\in {ℝ}^{[bq]\times [sp]}$

$\mathbf{\text{H}}(k=1)=\mathbf{\text{H}}=\left[\begin{array}{cccc}{\tilde{y}}_1& {\tilde{y}}_2& \cdots & {\tilde{y}}_s\\ {\tilde{y}}_2& {\tilde{y}}_3& \cdots & {\tilde{y}}_{s+1}\\ ⋮& ⋮& & ⋮\\ {\tilde{y}}_b& {\tilde{y}}_{b+1}& \cdots & {\tilde{y}}_{b+s-1}\end{array}\right]$

und einer verschobenen Block Hankelmatrix ${\mathbf{\text{H}}}^{\prime }\in {ℝ}^{[bq]\times [sp]}$

$\mathbf{\text{H}}(k=2)={\mathbf{\text{H}}}^{\prime }=\left[\begin{array}{cccc}{\tilde{y}}_2& {\tilde{y}}_3& \cdots & {\tilde{y}}_{s+1}\\ {\tilde{y}}_3& {\tilde{y}}_4& \cdots & {\tilde{y}}_{s+2}\\ ⋮& ⋮& & ⋮\\ {\tilde{y}}_{b+1}& {\tilde{y}}_{b+2}& \cdots & {\tilde{y}}_{b+s}\end{array}\right]$.

Die Bedeutung der Parameter $s,b\in \mathbb{N}\mathrm{\setminus }\{0\}$ wird im Folgenden noch geklärt. Die beiden Block Hankelmatrizen können nun durch die Beobachtbarkeitsmatrix $Q_B$ und die Steuerbarkeitsmatrix $Q_S$ausgedrückt werden. Es gilt schließlich für die Block Hankelmatrix und die verschobene Block Hankelmatrix

$Q_B{Q}_S=\left[\begin{array}{c}C\\ CA\\ ⋮\\ C{A}^{b-1}\end{array}\right]\left[\begin{array}{cccc}B& AB& \cdots & A^{s-1}B\end{array}\right]=\left[\begin{array}{cccc}CB& CAB& \cdots & C{A}^{s-1}B\\ CAB& C{A}^{2}B& \cdots & C{A}^{s}B\\ ⋮& ⋮& & ⋮\\ C{A}^{b-1}B& C{A}^{b}B& \cdots & C{A}^{s+b-2}B\end{array}\right]=\mathbf{\text{H}}$

und $Q_{B}A{Q}_S={\mathbf{\text{H}}}^{\prime }$respektive.

An dieser Stelle wird deutlich, dass $s$ und $b$ die Größe der Steuerbarkeitsmatrix und der Beobachtbarkeitsmatrix respektive festlegen und somit die maximal realisierbare Systemordnung $r$ des zu realisierende Systems begrenzen. Sie können frei gewählt werden und sind lediglich durch die Anzahl $N$ der Samples begrenzt, da gelten muss $b+s=N-1$. Eine mögliche Wahl ist

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

was der maximalen Größe für $s$ und $b$ entspricht und die Block Hankelmatrizen quadratisch macht. Mit anderen Worten begrenzt $N$ die Größe der Hankelmatrizen und somit letztendlich die Güte der Realisierung.

Der 4. Schritt ist die (ökonomische) Singulärwertzerlegung von $\mathbf{\text{H}}$

$\mathbf{\text{H}}=U\mathrm{\Sigma }{V}^{\mathrm{\top }}=\left[\begin{array}{cc}\tilde{U}& U_t\end{array}\right]\left[\begin{array}{cc}\tilde{\mathrm{\Sigma }}& \underset{\_}{0}\\ \underset{\_}{0}& {\mathrm{\Sigma }}_t\end{array}\right]\left[\begin{array}{c}{\tilde{V}}^{\mathrm{\top }}\\ V_t^{\mathrm{\top }}\end{array}\right]$,

mit $\tilde{U}\in {ℝ}^{[bq]\times r}$, $\tilde{\mathrm{\Sigma }}\in {ℝ}^{r\times r}$und ${\tilde{V}}^{\mathrm{\top }}\in {ℝ}^{r\times [sp]}$, wobei Tilde den Teil bezeichnet, der für die Realisierung verwendet wird und $t$ den Teil, der abgetrennt (truncated) wird. An dieser Stelle muss also eine Entscheidung über die Systemordnung $r$ des zu realisierenden Systems getroffen werden. Die Entscheidung wird auf Grundlage der Hankelsingulärwerte $\sigma$, die in $\mathrm{\Sigma }$ enthalten sind getroffen. Die Hankelsingulärwerte sind dabei ein Maß für die transportierte Energie von Eingang zu Ausgang des jeweiligen Zustands. Wie bereits erwähnt sind die Zustände balanciert. Dadurch gehören die größten $r$ Hankelsingulärtwerte ${\sigma }_i$ zu den am besten und in gleichen Maßen beobachtbaren und steuerbaren Zuständen. In der Regel lässt sich eine markante Position finden, ab der die Singulärwerte deutlich kleiner werden.

${\sigma }_1>{\sigma }_2>\dots >{\sigma }_r>>{\sigma }_{r+1}>\dots$

Ist dies nicht gegeben, kann die kumulierte Energie der Singulärwerte bestimmt werden und über eine Vorgabe der mindestens berücksichtigten Energie die Anzahl $r$ der dafür benötigten Zustände bestimmt werden.

Der fünfte Schritt ist die Bestimmung des realisierten Systems. Betrachtet man zunächst

$\begin{array}{cc}\mathbf{\text{H}}& =U\mathrm{\Sigma }{V}^{\mathrm{\top }}\\ & =U{\mathrm{\Sigma }}^{1/2}{\mathrm{\Sigma }}^{1/2}{V}^{\mathrm{\top }}\\ & =Q_B{Q}_S\end{array}$

folgt

$\begin{array}{cc}{Q}_B& =U{\mathrm{\Sigma }}^{1/2}\\ Q_S& ={\mathrm{\Sigma }}^{1/2}{V}^{\mathrm{\top }}.\end{array}$

Mit dem Zusammenhang

${\mathbf{\text{H}}}^{\prime }=Q_{B}A{Q}_S$

folgt schließlich für die Systemmatrix $A$

$\begin{array}{cc}A& =Q_B^{-1}{\mathbf{\text{H}}}^{\prime }{Q}_S^{-1}\\ & ={\mathrm{\Sigma }}^{-1/2}{U}^{\mathrm{\top }}{\mathbf{\text{H}}}^{\prime }V{\mathrm{\Sigma }}^{-1/2}.\end{array}$

Beachte, für die unitären Matrizen $U$ und $V$ gilt $U^{-1}=U^{\mathrm{\top }}$ und $V^{-1}=V^{\mathrm{\top }}$.

Außerdem folgen aus $\mathbf{\text{H}}=Q_B{Q}_S$ die Zusammenhänge

$\begin{array}{cc}{Q}_B& =\mathbf{\text{H}}{Q}_S^{-1}=\mathbf{\text{H}}V{\mathrm{\Sigma }}^{-1/2}=\left[\begin{array}{c}{𝐇}_{q,sb}\\ {𝐇}_{(b-1)q,sb}\end{array}\right]\left[\begin{array}{cc}\tilde{V}& V_t\end{array}\right]\left[\begin{array}{cc}{\tilde{\mathrm{\Sigma }}}^{-1/2}& \\ & {\mathrm{\Sigma }}_t^{-1/2}\end{array}\right]=\left[\begin{array}{c}C\\ CA\\ ⋮\\ C{A}^{b-1}\end{array}\right]\\ Q_S& =Q_B^{-1}\mathbf{\text{H}}={\mathrm{\Sigma }}^{-1/2}{U}^{\mathrm{\top }}\mathbf{\text{H}}=\left[\begin{array}{cc}{\tilde{\mathrm{\Sigma }}}^{-1/2}& \\ & {\mathrm{\Sigma }}_t^{-1/2}\end{array}\right]\left[\begin{array}{c}{\tilde{U}}^{\mathrm{\top }}\\ U_t^{\mathrm{\top }}\end{array}\right]\left[\begin{array}{cc}{𝐇}_{bq,p}& {𝐇}_{bq,(s-1)b}\end{array}\right]=\left[\begin{array}{cccc}B& AB& \dots & A^{s-1}B\end{array}\right],\end{array}$

in denen die Ausgangsmatrix $C$ und die Eingangsmatrix $B$ enthalten sind. An dieser Stelle wurde die Hankelmatrix entsprechend partitioniert, um später auf die Matrizen $B$ und $C$ zu kommen. Schließlich und endlich folgt das realisierte System,

$\begin{array}{cc}{\underset{\_}{\tilde{x}}}_{k+1}& =\tilde{A}{\underset{\_}{\tilde{x}}}_k+\tilde{B}{\underset{\_}{u}}_k\\ {\underset{\_}{\hat{y}}}_k& =\tilde{C}{\underset{\_}{\tilde{x}}}_k+\tilde{D}{\underset{\_}{u}}_k\end{array}$

mit den Zuständen $\underset{\_}{\tilde{x}}\in {ℝ}^r$, dem Eingang $\underset{\_}{u}\in {ℝ}^p$ und dem approximierten Ausgang $\underset{\_}{\hat{y}}\in {ℝ}^q$, sowie den Systemmatrizen

$\begin{array}{cc}\tilde{A}& ={\tilde{\mathrm{\Sigma }}}^{-1/2}{\tilde{U}}^{\mathrm{\top }}{\mathbf{\text{H}}}^{\prime }\tilde{V}{\tilde{\mathrm{\Sigma }}}^{-1/2}\\ \tilde{B}& ={\tilde{\mathrm{\Sigma }}}^{-1/2}{\tilde{U}}^{\mathrm{\top }}{\mathbf{\text{H}}}_{bq,p}\\ \tilde{C}& ={\mathbf{\text{H}}}_{q,sp}\tilde{V}{\tilde{\mathrm{\Sigma }}}^{-1/2}\\ \tilde{D}& ={\tilde{y}}_0,\end{array}$

mit $\tilde{A}\in {ℝ}^{r\times r}$, $\tilde{B}\in {ℝ}^{r\times p}$, $\tilde{C}\in {ℝ}^{q\times r}$ und $\tilde{D}\in {ℝ}^{q\times p}$.