Masse Feder System, Experimentelle Modellbildung (ERA)

In diesem Notebook wollen wir einen anderen Weg beschreiten. Statt mittels Lagrange II ein physikalisches Modell aufzustellen und die Parameter an die Messdaten anzupassen, wollen wir hier das Modell direkt aus Daten erstellen.

Die Experimentelle Modellbildung (aka Black-Box Systemidentifikation*) hat innerhalb der Ingenieurwissenschaften und der Regelungstheorie eine lange Tradition. Erste wissenschaftliche Aufsätze gehen auf die späten 1950er und frühen 1960er Jahre zurück.

Ähnliche Techniken sind derzeit innerhalb der künstlichen Intelligenz Forschung sehr beliebt. Im Englischen spricht man dann oft von Machine Learning oder Data Driven Modelling.

HoKalman / ERA

Ho und Kalman haben 1966 ihren Algorithmus vorgestellt. Die HoKalman Methode erlaubt das Lernen/Identifizieren von linearen zeitinvarianten Systemen (engl. LTI) direkt aus Messdaten. Kung hat 1978 diese Methode für den verrauschten Fall angepasst. Die moderne Version wurde 1985 von J.N. Juang und R. S. Pappa am NASA Langley Research Center vorgestellt. Ähnliche Versionen werden z.B von Luft- und Raumfahrtagenturen oder Autohersteller eingesetzt. Noch heute wird diese Methode innerhalb der Kontrolltheorie und dem maschinellen Lernen in wissenschaftlichen Aufsätzen untersucht.

Betrachten wir ein lineares zeitdiskretes System

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

mit den Signaldimensionen

${\mathbf{x}}_k\in {\mathbb{R}}^n,{\mathbf{u}}_k\in {\mathbb{R}}^p,{\mathbf{y}}_k\in {\mathbb{R}}^q$

und den Matrizen mit den Dimensionen

$\mathbf{A}\in {\mathbb{R}}^{n\times n},\mathbf{B}\in {\mathbb{R}}^{n\times p},\mathbf{C}\in {\mathbb{R}}^{q\times n},\mathbf{D}\in {\mathbb{R}}^{q\times p}.$

Für die ERA/HoKalman Methode benötigen wir Messdaten welche von einer Impulsantwort stammen. Den Eingangsimpulse wollen wir mit

$${\mathbf{u}}_k^{\delta }:=\mathbf{u}(k\Delta t)=\{\begin{array}{ll}\mathbf{I},& k=0\\ \mathbf{0},& k=1,2,3\end{array}$$

kennzeichnen und die Impulsantwort mit

$${\mathbf{y}}_k^{\delta }:=\mathbf{y}(k\Delta t)=\{\begin{array}{ll}\mathbf{D},& k=0\\ \mathbf{C}{\mathbf{A}}^{k-1}\mathbf{B},& k=1,2,3\end{array}.$$

Markov-Parameter / Impulsantwort

Als Nächstes wollen wir die Markov-Parameter bzw. die Impulsantwort verstehen.

Dazu setzen wir in die Ausgangsgleichung

$$\begin{array}{ccc}{\mathbf{y}}_k& =& \mathbf{C}{\mathbf{x}}_k+\mathbf{D}{\mathbf{u}}_k\end{array}$$

die zeitverschobene Zustandsgleichung

$$\begin{array}{ccc}{\mathbf{x}}_k& =& \mathbf{A}{\mathbf{x}}_{k-1}+\mathbf{B}{\mathbf{u}}_{k-1}\end{array}$$

ein. Damit erhalten wir eine Ausgangsgleichung welche statt ${\mathbf{x}}_k$ und ${\mathbf{u}}_k$

$$\begin{array}{c}{\mathbf{y}}_k=\mathbf{C}{\mathbf{x}}_k+\mathbf{D}{\mathbf{u}}_k=\mathbf{C}\mathbf{A}{\mathbf{x}}_{k-1}+\mathbf{C}\mathbf{B}{\mathbf{u}}_{k-1}+\mathbf{D}{\mathbf{u}}_k\end{array}$$

von ${\mathbf{x}}_{k-1}$, ${\mathbf{u}}_{k-1}$ und ${\mathbf{u}}_k$ abhängig ist. Wiederholtes anwenden führt auf eine Ausgangsgleichung welche nur noch von ${\mathbf{x}}_0$ und ${\mathbf{u}}_k,{\mathbf{u}}_{k-1},\dots ,{\mathbf{u}}_0$ abhängt.

Um diese Erkenntnis noch mal zu unterstreichen wiederholen wir dieses Vorgehen

$$\begin{array}{ccccc}{\mathbf{y}}_0& =& \mathbf{C}{\mathbf{x}}_0+\mathbf{D}{\mathbf{u}}_0& =& \mathbf{D}{\mathbf{u}}_0\\ {\mathbf{y}}_1& =& \mathbf{C}\mathbf{A}{\mathbf{x}}_0+\mathbf{C}\mathbf{B}{\mathbf{u}}_0+\mathbf{C}\mathbf{B}{\mathbf{u}}_1& =& \mathbf{C}\mathbf{B}{\mathbf{u}}_0+\mathbf{C}\mathbf{B}{\mathbf{u}}_1\\ {\mathbf{y}}_2& =& \mathbf{C}\mathbf{A}{\mathbf{x}}_1+\mathbf{C}\mathbf{B}{\mathbf{u}}_1+\mathbf{C}\mathbf{B}{\mathbf{u}}_2& =& \mathbf{C}\mathbf{A}\mathbf{B}{\mathbf{u}}_0+\mathbf{C}\mathbf{B}{\mathbf{u}}_1+\mathbf{C}\mathbf{B}{\mathbf{u}}_2\\ {\mathbf{y}}_3& =& \mathbf{C}\mathbf{A}{\mathbf{x}}_2+\mathbf{C}\mathbf{B}{\mathbf{u}}_2+\mathbf{C}\mathbf{B}{\mathbf{u}}_3& =& \mathbf{C}{\mathbf{A}}^2\mathbf{B}{\mathbf{u}}_0+\mathbf{C}\mathbf{A}\mathbf{B}{\mathbf{u}}_1+\mathbf{C}\mathbf{B}{\mathbf{u}}_2+\mathbf{D}{\mathbf{u}}_3\end{array}$$

wobei wir ${\mathbf{x}}_0=\mathbf{0}$ angenommen haben.

Die Ausgangsgleichung

$$\begin{array}{ccc}{\mathbf{y}}_k& =& \underset{\text{Markov-Parameter}}{\underbrace{\left[\begin{array}{ccccccc}\mathbf{D}& \mathbf{C}\mathbf{B}& \mathbf{C}\mathbf{A}\mathbf{B}& \mathbf{C}{\mathbf{A}}^2\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{(i-2)}\mathbf{B}& \mathbf{C}{\mathbf{A}}^{(i-1)}\mathbf{B}\end{array}\right]}}\left[\begin{array}{c}{\mathbf{u}}_k\\ {\mathbf{u}}_{k-1}\\ {\mathbf{u}}_{k-2}\\ {\mathbf{u}}_{k-3}\\ ⋮\\ {\mathbf{u}}_1\\ {\mathbf{u}}_0\end{array}\right]\end{array}$$

ergibt sich also durch eine Faltung der Markov-Parameter / (theoretische) Impulsantwort mit dem Eingang. Ist die Impulsantwort durch Daten gegeben wollen wir folgende Notation verwenden:

$$\begin{array}{ccc}{\mathbf{y}}_k& =& \underset{\text{Markov-Parameter}}{\underbrace{\left[\begin{array}{ccccccc}{\mathbf{y}}_0^{\delta }& {\mathbf{y}}_1^{\delta }& {\mathbf{y}}_2^{\delta }& {\mathbf{y}}_3^{\delta }& \cdots & {\mathbf{y}}_{k-1}^{\delta }& {\mathbf{y}}_k^{\delta }\end{array}\right]}}\left[\begin{array}{c}{\mathbf{u}}_k\\ {\mathbf{u}}_{k-1}\\ {\mathbf{u}}_{k-2}\\ {\mathbf{u}}_{k-3}\\ ⋮\\ {\mathbf{u}}_1\\ {\mathbf{u}}_0\end{array}\right]\end{array}$$

Diese Faltung können wir auch kurz mit

$$\begin{array}{c}{\mathbf{y}}_k=\sum _{i=0}^k{\mathbf{y}}_i^{\delta }{\mathbf{u}}_{k-i}\end{array}$$

anschreiben.

Steuerbarkeit und Beobachtbarkeit

Für die ERA/HoKalman benötigen wir die Konzepte der erweiterten Steuerbarkeits- und Beobachtbarkeitsmatrizen. Aus der linearen Kontrolltheorie kennen wir die Steuerbarkeitsmatrix

$${\mathcal{C}}_n=\left[\begin{array}{cccc}\mathbf{B}& \mathbf{A}\mathbf{B}& \cdots & {\mathbf{A}}^{n-1}\mathbf{B}\end{array}\right].$$

Die erweiterte Steuerbarkeitsmatrix

$${\mathcal{C}}_{m_c}=\left[\begin{array}{cccc}\mathbf{B}& \mathbf{A}\mathbf{B}& \cdots & {\mathbf{A}}^{m_c-1}\mathbf{B}\end{array}\right]$$

ist prinzipiell ähnlich aufgebaut, jedoch gilt immer $m_c>n$. In der Praxis besitzt die erweiterte Steuerbarkeitsmatrix eine deutlich größere Dimension $m_c>>n$ als die Steuerbarkeitsmatrix.

Die gleichen Überlegungen gelten für die Beobachtbarkeitsmatrix

$${\mathcal{O}}_n=\left[\begin{array}{c}\mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \cdots \\ \mathbf{C}{\mathbf{A}}^{n-1}\end{array}\right]$$

und für die erweiterte Beobachtbarkeitsmatrix

$${\mathcal{O}}_{m_o}=\left[\begin{array}{c}\mathbf{C}\\ \mathbf{C}\mathbf{A}\\ \cdots \\ \mathbf{C}{\mathbf{A}}^{m_o-1}\end{array}\right].$$

wobei wieder $m_o>>n$ gilt.

Diese Matrizen sind in vielen Methoden Systemidentifikation anzutreffen.

Block Hankel Matrizen

Zentral für viele Systemidentifikationsmethoden ist die sogenannte Block Hankel Matrix. Diese Matrix

$${\mathbf{H}}_{\mathrm{\infty },\mathrm{\infty }}=\left[\begin{array}{ccccc}{\mathbf{y}}_1^{\delta }& {\mathbf{y}}_2^{\delta }& {\mathbf{y}}_3^{\delta }& {\mathbf{y}}_4^{\delta }& \cdots \\ {\mathbf{y}}_2^{\delta }& {\mathbf{y}}_3^{\delta }& {\mathbf{y}}_4^{\delta }& {\mathbf{y}}_5^{\delta }& \cdots \\ {\mathbf{y}}_3^{\delta }& {\mathbf{y}}_4^{\delta }& {\mathbf{y}}_5^{\delta }& {\mathbf{y}}_6^{\delta }& \cdots \\ {\mathbf{y}}_4^{\delta }& {\mathbf{y}}_5^{\delta }& {\mathbf{y}}_6^{\delta }& {\mathbf{y}}_7^{\delta }& \cdots \\ ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right]$$

wird durch die Impulsantwort konstruiert. Die Impulsantwort wird in die erste Zeile geschrieben und dann um einen Zeitschritt versetzt in die Zweite, und so weiter und so fort.

Weiters gelten die Beziehungen

$$\begin{array}{c}{\mathbf{H}}_{\mathrm{\infty },\mathrm{\infty }}={\mathcal{O}}_{\mathrm{\infty }}{\mathcal{C}}_{\mathrm{\infty }}={\mathcal{O}}_{\mathrm{\infty }}\mathbf{T}{\mathbf{T}}^{-1}{\mathcal{C}}_{\mathrm{\infty }}\end{array}$$

und

$${\mathbf{y}}_{+}=\mathbf{H}{\mathbf{u}}_{-}$$

wobei die Vektoren wie folgt konstruiert

$${\mathbf{y}}_{+}={\left[\begin{array}{ccc}{\mathbf{y}}_0& {\mathbf{y}}_1& \cdots \end{array}\right]}^T;{\mathbf{u}}_{-}={\left[\begin{array}{ccc}{\mathbf{u}}_{-1}& {\mathbf{u}}_{-2}& \cdots \end{array}\right]}^T$$

sind. Die Hankel Matrix verbindet alle vergangenen Eingänge mit allen zukünftigen Ausgängen. Die Matrix $\mathbf{T}$ ist eine beliebige unitäre Transformationsmatrix.

In der Praxis operieren wir immer mit einer endlichen Anzahl von Messdaten

$$\begin{array}{c}{\mathbf{H}}_{m_o,m_c}={\mathcal{O}}_{m_o}{\mathcal{C}}_{m_c}={\mathcal{O}}_{m_o}\mathbf{T}{\mathbf{T}}^{-1}{\mathcal{C}}_{m_c}\end{array}.$$

Durch das Einsetzen der Markov-Parameter in die Hankel Matrix

$$\begin{array}{rl}{\mathbf{H}}_{m_o,m_c}& =\left[\begin{array}{cccc}{\mathbf{y}}_1^{\delta }& {\mathbf{y}}_3^{\delta }& \cdots & {\mathbf{y}}_{m_c}^{\delta }\\ {\mathbf{y}}_2^{\delta }& {\mathbf{y}}_3^{\delta }& \cdots & {\mathbf{y}}_{m_c+1}^{\delta }\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{y}}_{m_o}^{\delta }& {\mathbf{y}}_{m_o+1}^{\delta }& \cdots & {\mathbf{y}}_{m_o+m_c-1}^{\delta }\end{array}\right]\\ & =\left[\begin{array}{cccc}\mathbf{C}\mathbf{B}& \mathbf{C}\mathbf{A}\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{m_c-1}\mathbf{B}\\ \mathbf{C}\mathbf{A}\mathbf{B}& \mathbf{C}{\mathbf{A}}^2\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{m_c}\mathbf{B}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{C}{\mathbf{A}}^{m_o-1}\mathbf{B}& \mathbf{C}{\mathbf{A}}^{m_o}\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{m_o+m_c-2}\mathbf{B}\end{array}\right]\\ & =\left[\begin{array}{c}\mathbf{C}\\ \mathbf{C}\mathbf{A}\\ ⋮\\ \mathbf{C}{\mathbf{A}}^{m_o-1}\end{array}\right]\left[\begin{array}{cccc}\mathbf{B}& \mathbf{A}\mathbf{B}& \cdots & {\mathbf{A}}^{m_c-1}\mathbf{B}\end{array}\right]\\ & ={\mathcal{O}}_{m_o}{\mathcal{C}}_{m_c}\end{array}$$

können wir die Beziehung ${\mathbf{H}}_{m_o,m_c}={\mathcal{O}}_{m_o}{\mathcal{C}}_{m_c}$ leicht erkennen. Diese Zerlegung erfolgt mittels einer Singulärwertzerlegung ${\mathcal{O}}_{m_o}{\mathcal{C}}_{m_c}=SVD({\mathbf{H}}_{m_o,m_c})$ der Hankel Matrix.

Sind die beiden Matrizen ${\mathcal{O}}_{m_o}$ und ${\mathcal{C}}_{m_c}$ vorhanden so können die Matrizen $\mathbf{C}$ und $\mathbf{B}$ leicht selektiert werden, weil die Eingangs- und Ausgangsdimensionen bekannt sind.

Für die Rekonstruktion $\mathbf{A}$ von benötigen wir eine weiter Hankel Matrix

$$\begin{array}{rl}{\mathbf{H}}_{m_o,m_c}^{{}^{\prime }}& =\left[\begin{array}{cccc}{\mathbf{y}}_2^{\delta }& {\mathbf{y}}_3^{\delta }& \cdots & {\mathbf{y}}_{m_c+1}^{\delta }\\ {\mathbf{y}}_3^{\delta }& {\mathbf{y}}_4^{\delta }& \cdots & {\mathbf{y}}_{m_c+2}^{\delta }\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{y}}_{m_o+1}^{\delta }& {\mathbf{y}}_{m_o+2}^{\delta }& \cdots & {\mathbf{y}}_{m_o+m_c}^{\delta }\end{array}\right]\\ & =\left[\begin{array}{cccc}\mathbf{C}\mathbf{A}\mathbf{B}& \mathbf{C}{\mathbf{A}}^2\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{m_c}\mathbf{B}\\ \mathbf{C}{\mathbf{A}}^2\mathbf{B}& \mathbf{C}{\mathbf{A}}^3\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{m_c+1}\mathbf{B}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{C}{\mathbf{A}}^{m_o}\mathbf{B}& \mathbf{C}{\mathbf{A}}^{m_o+1}\mathbf{B}& \cdots & \mathbf{C}{\mathbf{A}}^{m_o+m_c-1}\mathbf{B}\end{array}\right]\\ & =\left[\begin{array}{c}\mathbf{C}\\ \mathbf{C}\mathbf{A}\\ ⋮\\ \mathbf{C}{\mathbf{A}}^{m_o-1}\end{array}\right]\mathbf{A}\left[\begin{array}{cccc}\mathbf{B}& \mathbf{A}\mathbf{B}& \cdots & {\mathbf{A}}^{m_c-1}\mathbf{B}\end{array}\right]\\ & ={\mathcal{O}}_{m_o}\mathbf{A}{\mathcal{C}}_{m_c}\end{array}$$

welche jedoch um einen Zeitschritt verschoben ist. Dadurch gewinnen wir die Beziehung ${\mathbf{H}}_{m_o,m_c}^{{}^{\prime }}={\mathcal{O}}_{m_o}\mathbf{A}{\mathcal{C}}_{m_c}$.

Durch die Verwendung der Moore-Penrose-Inverse erfolgt ein Umstellen der obigen Gleichung auf $\mathbf{A}={\mathcal{O}}_{m_o}^{†}{\mathbf{H}}_{m_o,m_c}^{{}^{\prime }}{\mathcal{C}}_{m_c}^{†}$ womit die Rekonstruktion von $\mathbf{A}$ gelungen ist.

Die Beziehung $\mathbf{D}={\mathbf{y}}_0^{\delta }$ ist aus der Definition der Markov Parameter klar ersichtlich.

Hinweis

Eigensystem Realization Algorithm:

• Konstruieren der Hankel-Matrizen ${\mathbf{H}}_{m_o,m_c}$ und ${\mathbf{H}}_{m_o,m_c}^{\prime }$
• Durchführung der SVD der Hankel-Matrix

$${\mathbf{H}}_{m_o,m_c}=\left[\begin{array}{cc}{\mathbf{U}}_s& {\mathbf{U}}_n\end{array}\right]\left[\begin{array}{cc}{\mathbf{\Sigma }}_s& \mathbf{0}\\ \mathbf{0}& {\mathbf{\Sigma }}_n\end{array}\right]\left[\begin{array}{c}{\mathbf{V}}_s^T\\ {\mathbf{V}}_n^T\end{array}\right]$$

• Schätzung der Matrizen aus den Daten

  $$\begin{array}{rl}& {\mathcal{O}}_{m_o}={\mathbf{U}}_s{\mathbf{\Sigma }}_s^{1/2}\mathbf{T}\\ & {\mathcal{C}}_{m_c}={\mathbf{T}}^{-1}{\mathbf{\Sigma }}_s^{1/2}{\mathbf{V}}_s^T\end{array}$$

• Schätzungen der Zustandsraummatrizen anhand obiger Gleichungen

  $$\begin{array}{rl}& \mathbf{A}={\mathbf{\Sigma }}_s^{-1/2}{\mathbf{U}}_s^T{\mathbf{H}}_{m_o,m_c}^{\prime }{\mathbf{V}}_s{\mathbf{\Sigma }}_s^{-1/2}\\ & \mathbf{B}={\mathcal{C}}_{m_c}(:,1:q)\\ & \mathbf{C}={\mathcal{O}}_{m_o}(1:p,:)\\ & \mathbf{D}={\mathbf{y}}_0^{\delta }\end{array}$$

Code Beispiel

In diesen Abschnitt wollen wir den Algorithmus in Python implementieren und an dem Beispiel eines Masse Feder Systems testen.

Masse Feder System (2-FG)

Gegeben ist wieder das bekannte System wobei hier nur die Positionen der Massen gemessen werden können. D.h. es kann nur ein Teil der Zustände gemessen werden. (Partially Observable Markov Decision Process (POMDP), Teilweise beobachtbarer Markov-Entscheidungsprozess)

$$A=\left[\begin{array}{cccc}0& 0& 1& 0\\ 0& 0& 0& 1\\ -\frac{c_1+c_2}{m_1}& \frac{c_2}{m_1}& -\frac{d_1+d_2}{m_1}& \frac{d_2}{m_1}\\ \frac{c_2}{m_2}& \frac{-c_2-c_3}{m_2}& \frac{d_2}{m_2}& \frac{-d_2-d_3}{m_2}\end{array}\right]$$

$$B=\left[\begin{array}{cc}0& 0\\ 0& 0\\ \frac{1}{m_1}& 0\\ 0& \frac{1}{m_2}\end{array}\right]$$

$$C=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\end{array}\right]$$

$$D=\left[\begin{array}{cc}0& 0\\ 0& 0\end{array}\right]$$

Erstelle das System

Wir wollen ein zeitdiskretes System verwenden welches aus dem zeitkontinuierlichen System durch die scipy Methode to_discrete erstellt wird.

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt

m1 = 1
d1 = 0.15
c1 = 0.3

m2 = 1
d2 = 0.15
c2 = 0.3

d3 = 0.15
c3 = 0.3

A = np.array([[0., 0., 1., 0.],
              [0., 0., 0., 1.], 
              [-(c1+c2)/m1, c2/m1, -(d1+d2)/m1, d2/m1],
              [c2/m2, -(c2+c3)/m2, (d2)/m2, -(d2+d3)/m2]])
B = np.array([[0., 0.], [0., 0.], [1./m1, 0.], [0., 1/m2]])
C = np.array([[1., 0., 0., 0], [0., 1., 0., 0.]])
D = np.array([[0., 0.],[0., 0.]])

C

array([[1., 0., 0., 0.],
       [0., 1., 0., 0.]])

sys_c = signal.StateSpace(A, B, C, D)
print(sys_c)

StateSpaceContinuous(
array([[ 0.  ,  0.  ,  1.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  1.  ],
       [-0.6 ,  0.3 , -0.3 ,  0.15],
       [ 0.3 , -0.6 ,  0.15, -0.3 ]]),
array([[0., 0.],
       [0., 0.],
       [1., 0.],
       [0., 1.]]),
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.]]),
array([[0., 0.],
       [0., 0.]]),
dt: None
)

eig_values, _ = np.linalg.eig(sys_c.A)
eig_values

array([-0.225+0.92161543j, -0.225-0.92161543j, -0.075+0.54256336j,
       -0.075-0.54256336j])

dt = 0.1
sys_d = sys_c.to_discrete(dt)
print(sys_d)

StateSpaceDiscrete(
array([[ 9.97038953e-01,  1.46889167e-03,  9.84204432e-02,
         7.83673580e-04],
       [ 1.46889167e-03,  9.97038953e-01,  7.83673580e-04,
         9.84204432e-02],
       [-5.88171638e-02,  2.90559288e-02,  9.67630371e-01,
         1.59968561e-02],
       [ 2.90559288e-02, -5.88171638e-02,  1.59968561e-02,
         9.67630371e-01]]),
array([[4.94800256e-03, 2.58485012e-05],
       [2.58485012e-05, 4.94800256e-03],
       [9.84204432e-02, 7.83673580e-04],
       [7.83673580e-04, 9.84204432e-02]]),
array([[1., 0., 0., 0.],
       [0., 1., 0., 0.]]),
array([[0., 0.],
       [0., 0.]]),
dt: 0.1
)

eig_values, _ = np.linalg.eig(sys_d.A)
eig_values

array([0.97360179+0.08998355j, 0.97360179-0.08998355j,
       0.99106754+0.05382452j, 0.99106754-0.05382452j])

Erzeuge Messdaten für ERA

Wie theoretisch besprochen, benötigen wir eine Impulsantwort des Systems. In der Praxis werden die Impulsantworten für jeden Eingang separat erfasst und dann gestapelt.

Weiters überlagern wir die Impulsantwort mit einem Messrauschen.

t = np.arange(0,50,dt)
print(t.shape)

(500,)

u1 = np.zeros((len(t),2))
print(u1.shape)
u1[0,0] = 1
print(u1[:6,:])

u2 = np.zeros((len(t),2))
print(u2.shape)
u2[0,1] = 1
print(u2[:6,:])

(500, 2)
[[1. 0.]
 [0. 0.]
 [0. 0.]
 [0. 0.]
 [0. 0.]
 [0. 0.]]
(500, 2)
[[0. 1.]
 [0. 0.]
 [0. 0.]
 [0. 0.]
 [0. 0.]
 [0. 0.]]

tout1, yout1_, xout1 = signal.dlsim(sys_d,u1)
print(yout1_.shape)

tout2, yout2_, xout2 = signal.dlsim(sys_d,u2)
print(yout2_.shape)

(500, 2)
(500, 2)

yout1 = yout1_ + np.random.randn(len(yout1_),2)*0.01
yout2 = yout2_ + np.random.randn(len(yout2_),2)*0.01

fig = plt.figure(figsize=(10, 5))

plt.subplot(2,2,1)
plt.plot(tout1,u1)
plt.grid()

plt.subplot(2,2,2)
plt.plot(tout1,yout1)
plt.grid()

plt.subplot(2,2,3)
plt.plot(tout1,u2)
plt.grid()

plt.subplot(2,2,4)
plt.plot(tout2,yout2)
plt.grid()

yout = np.stack((yout1,yout2))
yout.shape

(2, 500, 2)

YY = yout.transpose(2,0,1)
YY.shape

(2, 2, 500)

Systemidentifizierung mittels ERA

def BlockHankel(Y, mo, mc):
    
    nr, nc, L = Y.shape
    
    H = np.zeros((nr*mo,nc*mc))
    
    for r in range(mo):
        for c in range(mc):
            rowBeg = nr*r
            rowEnd = nr*(r+1)
            colBeg = nc*c
            colEnd = nc*(c+1)
            
            #print(rowBeg,rowEnd)
            #print(colBeg,colEnd)
            if (r+c<L):
                H[rowBeg:rowEnd, colBeg:colEnd] = Y[:,:,r+c]

    return H

def ERA(Y, mo, mc, r):
    
    nr, nc, L = Y.shape
    
    Hr = BlockHankel(Y[:,:,1:], mo, mc)
    Hr_p = BlockHankel(Y[:,:,2:], mo, mc)
    
    nx = np.linalg.matrix_rank(Hr)
    
    U,S,Vh = np.linalg.svd(Hr, True)
    S_mat = np.diag(S)
    Sigma = S_mat[0:r,0:r]
    Ur =U[:,0:r]
    Vhr =Vh[0:r,:]
    
    Sigma_inv = np.linalg.inv(Sigma)
    
    Ar = Sigma_inv**(0.5) @ Ur.T @ Hr_p @ Vhr.T @ Sigma_inv**(0.5)
    Br = Sigma_inv**(0.5) @ Ur.T @ Hr[:,0:nr]
    Cr = Hr[0:nc,:] @ Vhr.T @ Sigma_inv**(0.5)
    
    Dr = Y[:,:,0]
    
    HSVs = S
        
    return Ar, Br, Cr, Dr, HSVs, nx

r = 4
Ar, Br, Cr, Dr, HSVs, nx = ERA(YY,mo=100,mc=100,r=r)
nx

200

sysr = signal.StateSpace(Ar,Br,Cr,Dr,dt=dt)
print(sysr)

StateSpaceDiscrete(
array([[ 9.90596372e-01,  5.41136505e-02,  8.04812446e-04,
         6.29245067e-04],
       [-5.41254292e-02,  9.91365664e-01, -1.18051982e-03,
        -1.22550639e-03],
       [-1.26784150e-03, -1.12264851e-03,  9.79708374e-01,
         8.95006660e-02],
       [-3.07530892e-04,  7.77785959e-04, -8.95483693e-02,
         9.66642986e-01]]),
array([[ 0.21456143,  0.22054022],
       [ 0.21562074,  0.2121428 ],
       [-0.15103836,  0.1746997 ],
       [-0.15850599,  0.15819699]]),
array([[ 0.21590535, -0.20778459, -0.17197816,  0.15145799],
       [ 0.21912918, -0.21986231,  0.15552216, -0.16481435]]),
array([[0.01162235, 0.01012659],
       [0.00740089, 0.01118258]]),
dt: 0.1
)

tout1_simr,yout1_simr, xout1_simr = signal.dlsim(sysr,u1)
tout2_simr,yout2_simr, xout2_simr = signal.dlsim(sysr,u2)

fig = plt.figure(figsize=(10, 8))

plt.subplot(2,2,1)
plt.plot(tout1,yout1,label='Orginal System')
plt.plot(tout1_simr,yout1_simr,linestyle='--',label='Sysid ERA')
plt.grid()
plt.legend()

plt.subplot(2,2,2)
plt.plot(tout2,yout2,label='Orginal System')
plt.plot(tout2_simr,yout2_simr,linestyle='--',label='Sysid ERA')
plt.grid()
plt.legend()

Validierung

Eine in der Praxis übliche Methode der Validierung verlangt neue Messdaten zu sammeln, welche nicht für das Identifizieren des Systems benutzt wurden. Oft genügen 2 neue Signale. Hier verwenden wir einen sin und einen step. Der Vorteil dieser Methode ist, dass die identifizierten Modelle sehr gut generalisieren müssen um diesen Test zu bestehen.

Sowohl das Feld der Systemidentifikation als auch das Feld des Maschinellen Lernens beinhalten weitere Theorien und Methoden wie man die Güte eines gelernten Systems validiert.

Validierung 1

• u = sin(t)

t_sim = np.arange(0,50,dt)
print(t_sim.shape)

(500,)

u1 = np.zeros((len(t_sim),2))
u1[:,0] = np.sin(t_sim)
tout1, yout1, xout1 =signal.dlsim(sys_d,u1)
tout1_simr,yout1_simr, xout1_simr = signal.dlsim(sysr,u1)

plt.subplots(1, 1, figsize = (10, 5))
plt.subplot(1,2,1)
plt.plot(tout1,yout1[:,0],label='True')
plt.plot(tout1_simr,yout1_simr[:,0],linestyle='--',label='Sysid ERA')
plt.grid()
plt.legend()
plt.subplot(1,2,2)
plt.plot(tout1,yout1[:,1],label='True')
plt.plot(tout1_simr,yout1_simr[:,1],linestyle='--',label='Sysid ERA')
plt.grid()
plt.legend()

Validierung 2

• u = step(t)

u1 = np.ones((len(t_sim),2))
u1.shape

(500, 2)

# Input sequence
u1 = np.ones((len(t_sim),2))
u1[:10,:] = 0
tout1, yout1, xout1 =signal.dlsim(sys_d,u1)
tout1_simr,yout1_simr, xout1_simr = signal.dlsim(sysr,u1)

plt.subplots(1, 1, figsize = (10, 5))
plt.subplot(1,2,1)
plt.plot(tout1,yout1[:,0],label='True')
plt.plot(tout1_simr,yout1_simr[:,0],linestyle='--',label='Sysid ERA')
plt.grid()
plt.legend()
plt.subplot(1,2,2)
plt.plot(tout1,yout1[:,1],label='True')
plt.plot(tout1_simr,yout1_simr[:,1],linestyle='--',label='Sysid ERA')
plt.grid()
plt.legend()

Fazit

Wir haben eine experimentelle Modellbildung mit der Hilfe der HoKalman/ERA Methode durchgeführt. Dadurch ist es uns gelungen ein Modell $\left(\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}\right)$ direkt aus Messdaten zu identifizieren.

Referenzen

• Wikipedia: Eigensystem Realisierungsalgorithmus
• Wikipedia: Systemidentifikation
• The rise of data-driven modelling, Nature Reviews Physics, June 2021
• Samet Oymak and Necmiye Ozay, “Non-asymptotic identification of lti systems from a single trajectory”, 2019 American control conference (ACC), 2019.