H2-Norm und Kovarianz lineare Systeme

In Octave oder Matlab steht uns mit covar eine Methode zur Verfügung, mit welcher wir die stationäre Kovarianz des Zustandsvektors $x$ und die stationäre Kovarianz des Ausgangsvektors $y$ berechnen können. In Python existiert derzeit eine solche Methode weder in scipy noch in control.

Dabei ist covar wieder mit der ${\mathcal{H}}_2$-Norm verbunden.

Kovarianz eines zeitdiskretes Systems

Gegeben sei ein lineares System, das von Eingängen mit weißem Gauß’schen Rauschen gesteuert wird. Dieses Rauschen ist also kein Kontrolleingang sondern eine Störung. In Form eines Zustandsraums wird dieses System wie folgt dargestellt:

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

wobei $A\in {\mathbb{R}}^{n\times n}$, $B\in {\mathbb{R}}^{n\times p}$, $C\in {\mathbb{R}}^{q\times n}$ und $D\in {\mathbb{R}}^{q\times p}$. $n$, $m$ und $p$ beschreiben die Anzahl der Zustände, Eingänge bzw. Ausgänge des linearen Systems.

Bei den Eingangsstörungen $w$ handelt es sich um weißes Gauß’sches Rauschen mit der Spektraldichte ${\mathrm{\Sigma }}_w\in {\mathbb{R}}^{m\times m}$ mit Nullmittelwert. Anders ausgedückt: Der Erwartungswert von $w$ ist $\mathbb{E}(w_k)=0$ und seine Kovarianzmatrix ist $\mathbb{E}(w_k{w}_l^T)={\mathrm{\Sigma }}_w{\delta }_{kl}$, wobei das Kronecker-Symbol duch

$${\delta }_{kl}=\{\begin{array}{l}1\text{für }k=l\\ 0\text{für }k\ne l\end{array}$$

gegeben ist.

Analyse stochastischer linearer Systeme

Wir wollen hier nun stochastische zeitdiskrete zeitinvariante lineare Systeme analysieren. Dafür verwenden wir ein etwas andere Notation. Um den zeitlichen Zusammenhang stärker herauszustellen schreiben wir statt $x_k$, $x[k]$.

Gegeben sei das System

$$\begin{array}{rl}x[k+1]& =Ax[k]+Bw[k].\end{array}$$

Der Anfangszustand $x[0]$ ist ein Zufallsvektor, dessen Erwartungswert und Kovarianzmatrix wir kennen,

$$\begin{array}{r}\mathbb{E}\{x[0]\}={\mu }_x[0]=0\mathbb{E}\{[x[0]-{\mu }_x[0]][x[0]-{\mu }_x[0]{]}^T\}={\mathrm{\Sigma }}_x[0].\end{array}$$

Der Eingangsvektor $\left\{w_k\right\}$ ist ein (vom Anfangszustand $x[0]$ unabhäniges) weißes Rauschen, dessen momentanen Erwartungswert und Autokovarianzmatrix wir kennen,

$$\begin{array}{rlr}& \mathbb{E}\{w[k]\}={\mu }_w[k]=0& \text{ für }k\ge 0\\ & \mathbb{E}\{[w[k]-{\mu }_w[k]][w[l]-{\mu }_w[l]{]}^T\}={\mathrm{\Sigma }}_w{\delta }_{kl}& \text{ für }k,l\ge 0\text{ mit }{\mathrm{\Sigma }}_w={\mathrm{\Sigma }}_w^T\ge 0\\ & \mathbb{E}\{[w[k]-{\mu }_w[k]][x[0]-{\mu }_x[0]{]}^T\}=0.& \text{ für }k\ge 0\text{ (Unkorreliertheit)}.\end{array}$$

Der momentane Erwartungswert des Zustandsvektors kann offenbar aus der folgenden zeitdiskreten Bewegungsgleichung berechnet werden:

$$\begin{array}{rl}{\mu }_x[k+1]& =\mathbb{E}(A{x}_k+B{u}_k)=A{\mu }_x[k]+B{\mu }_w[k]=0\\ {\mu }_x[0]& =\mathbb{E}\{x[0]\}=0\end{array}$$

Da der Vektor-Zufallsprozeß $\left\{w_k\right\}$ ein weißes Rauschen und vom Anfangsvekor $x[0]$ unabhängig ist, kann die momentane Kovarianzmatrix ${\mathrm{\Sigma }}_k$ des Zustandsvektors $x_k$ ebenfalls rekursiv berechnet werden:

$$\begin{array}{rl}{\mathrm{\Sigma }}_x[k+1]& =\mathbb{E}\left\{[x[k+1]-{\mu }_x[k+1]]{[x[k+1]-{\mu }_x[k+1]]}^T\right\}\\ & =\mathbb{E}\left\{[A(x[k]-{\mu }_x[k])+B(w[k]-{\mu }_w[k])]{[A(x[k]-{\mu }_x[k])+B(w[k]-{\mu }_w[k])]}^T\right\}\\ & =A{\mathrm{\Sigma }}_x[k]A^T+B{\mathrm{\Sigma }}_w{B}^T\\ {\mathrm{\Sigma }}_x[0]& ={\mathrm{\Sigma }}_x[0]\end{array}$$

Für eine stabile Matrizen $A$ konvergiert diese Gleichung gegen die Lyapunov Gleichung

$${\mathrm{\Sigma }}_x=A{\mathrm{\Sigma }}_x{A}^T+B{\mathrm{\Sigma }}_w{B}^T=0.$$

Die Kovarianzmatrix des Ausgangs $y$

$$\begin{array}{rl}{\mathrm{\Sigma }}_y[k]& =\mathbb{E}\{y[k]y[k{]}^T\}\\ & =\mathbb{E}\{[Cx[k]+Dw[k]][Cx[k]+Dw[k]{]}^T\}\\ & =C{\mathrm{\Sigma }}_x[k]C^T+D{\mathrm{\Sigma }}_w{D}^T\end{array}$$

konvergiert für stabile System zu

$${\mathrm{\Sigma }}_y=C{\mathrm{\Sigma }}_x{C}^T+D{\mathrm{\Sigma }}_w{D}^T=0.$$

Für instabile System (Spektralradius $\rho (A)>1$) sind die Kovarinazmatrizen unendlich.

Die Verbindung zur ${\mathcal{H}}_2$-Norm ist nun mit

$$\Vert G{\Vert }_2^2=\underset{k\to \mathrm{\infty }}{\text{lim}}\mathbb{E}(y_k{y}_k^T)=tr(covar(G))=tr({\mathrm{\Sigma }}_y)$$

gegeben. Die ${\mathcal{H}}_2$-Norm ist ein Mass für die Energie am Ausgangs eines Systems, welches von einem weißen Gauß’schen Rauschen am Eingang angetrieben wird.

Tipp

Wir verwenden im Code statt $\mathrm{\Sigma }$ die Buchstaben $P$ und $W$. Ähnliche Notationen werden von vielen Programmpaketen verwendet:

$$\begin{gathered} {\mathrm{\Sigma }}_x=P_x \\ {\mathrm{\Sigma }}_y=P_y \\ {\mathrm{\Sigma }}_w=W \end{gathered}$$

Zusammengefasst lässt sich ein Zustandsraumverfahren wie folgt berechnen:

1. Berechne eine Zustandsrealisierung aus der Übertragungsfunktion

$$G(z):=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]$$

2. Löse die Lyapunov Gleichung, (dabei muss gelten $P_x\ge 0$)

$$\begin{array}{rl}{P}_x& =A{P}_x{A}^T{P}_x+BW{B}^T=0\\ P_y& =C{P}_x{C}^T+DW{D}^T\end{array}$$

Code

import numpy as np
import scipy.linalg as linalg
import scipy.signal as signal


import warnings
warnings.filterwarnings("error")

def covar_d(A,B,C,D,W):
    """ naive implementation of covar for linear time-discrete systems (d-LTI)
        (there might be better ways to do that)
    """
    
    if np.all(np.abs(linalg.eigvals(A)) <= 1):
        Px = linalg.solve_discrete_lyapunov(A.T, B@W@B.T)
        Py = C@Px@C.T + D@W@D.T   
    else:
        Px = np.inf
        Py = np.inf
    return [Px, Py]

def covar_d_data(tout,yout,xout):
    """ data based approximation of covar for linear time-discrete systems (d-LTI)
        (there might be better ways to do that)
    """
    
    # catch overflow
    try:
        Px = xout.T@xout/len(tout)
        Py = yout.T@yout/len(tout)
    except RuntimeWarning:
        Px = np.inf
        Py = np.inf
    return [Px, Py]

Stabiles System

Ad_stable = np.array([[0.5, 0.1],[0.1, 0.5]])
print(np.abs(linalg.eigvals(Ad_stable)) <= 1.)

Bd = np.array([[0.],[1.]])
Cd = np.eye(2)*0.5
Dd = np.zeros((2,1))
Wd = np.eye(1)*5

Px, Py = covar_d(Ad_stable,Bd,Cd,Dd,Wd)
print(Px)
print(Py)

[ True  True]
[[0.15174655 0.46502976]
 [0.46502976 6.73069392]]
[[0.03793664 0.11625744]
 [0.11625744 1.68267348]]

Wir wollen die modellbasierte Variante mit einer datenbasierten Variante vergleichen. Dazu erzeugen wir Eingangsdaten und Ausgangsdaten.

G_sys = signal.StateSpace(Ad_stable, Bd, Cd, Dd,  dt=1.)
w_in = np.sqrt(Wd)*np.random.randn(100000)
tout, yout, xout = signal.dlsim(G_sys, w_in.T)

Px_data, Py_data = covar_d_data(tout, yout, xout)
print(Px_data)
print(Py_data)

[[0.15053712 0.45885585]
 [0.45885585 6.70173513]]
[[0.03763428 0.11471396]
 [0.11471396 1.67543378]]

print(linalg.norm(Px - Px_data)/linalg.norm(Px))
print(linalg.norm(Py - Py_data)/linalg.norm(Py))

0.004474953153232987
0.004474953153232987

Beide Varianten stimmen sehr gut überein. Die Differenz erklärt sich einerseits durch die endlichen Daten und die numerische Ungenaugikeit. np.float128 könnte bessere Ergebnisse liefern.

Instabiles System

Ad_unstable = np.array([[1.5, 0.1], [0.1, 1.5]])
print(np.abs(linalg.eigvals(Ad_unstable)) < 1.)

Px, Py = covar_d(Ad_unstable,Bd,Cd,Dd,Wd)
print(Px)
print(Py)

[False False]
inf
inf

G_sys = signal.StateSpace(Ad_unstable, Bd, Cd, Dd, dt=1.)
w_in = np.sqrt(Wd)*np.random.randn(100000)
tout, yout, xout = signal.dlsim(G_sys, w_in.T)

Px_data, Py_data = covar_d_data(tout,yout,xout)
print(Px_data)
print(Py_data)

inf
inf

Die Kovarianmatrizen werden zu groß und erzeugen deshalb eine Laufzeitwarnung, welche abgefangen wird.

Fazit

Die Kovarianzmatrizen kommen in vielen wissenschaftlichen Aufsätzen zum Einsatz, vor allem in Aufsätzen in denen Schätzer analysiert werden (Kalmanfilter, Systemidentifikation, Adaptive Regelung, Bestärkendes Lernen, …).

Die Kovarianzmatrizen stehen mit der ${\mathcal{H}}_2$ in Verbindung, was für die Implementierungsebene von Bedeutung ist.

Referenzen

• Optimal Sampled-Data Control Systems (Tongwen Chen und Bruce Francis, 1994)
• Essentials of Robust Control (Kemin, Zhoun, John C. Doyle, 1997)
• Regelungstechnik, Mathematische Grundlagen, Entwurfsmethoden, Beispiele (Geering, 2004)