H2-Norm linearer zeitdiskreter zeitinvarianter Systeme (d-LTI)

Die \(\mathcal{H}_2\)-Norm für lineare Systeme spielt eine große Rolle in der Systemanalyse, Reglersynthese und Modellreduktion. Auch im bestärkenden Lernen ist diese Norm von Bedeutung. Zum Beispiel werden Ansätze der \(\mathcal{H}_2\)-modellbasierten Reglersynthese, auf das \(\mathcal{H}_2\)-modellfreie bestärkende Lernen übertragen.

“In der Mathematik und Kontrolltheorie ist \(\mathcal{H}_2\) oder \(\mathcal{H}\)-Quadrat ein Hardy-Raum mit quadratischer Norm. Er ist ein Unterraum des \(\mathcal{L}_2\)-Raums (Lebesgue-Raum) und somit ein Hilbert-Raum. Insbesondere ist er ein reproduzierender Kern-Hilbertraum.” (wikipedia)

Regelungspakete in Octave oder Matlab bieten Methoden wie norm(sys, 2) um diese Norm zu berechnen. In Python steht derzeit ein solche Methode weder in scipy noch in control zur Verfügung. Wir werden hier eine naive Implementierung zeigen. Naiv bedeutet, dass es sich nicht zwingend um die beste Implementierung bezüglich Effizient, Exaktheit und Eleganz handelt. Um eine gute Implementierung sicherzustellen, müssen die wissenschaftlichen Aufsätze und Referenz-Implementierung studiert werden.

\(\mathcal{H}_2\)-Norm eines zeitdiskretes Systems

Die \(\mathcal{H}_2\)-Norm einen zeitdiskretes Systems mit der Übertragungsfunktion \(G(z)\) ist gegeben durch

\[ \Vert G \Vert_2 = \sqrt{tr\left[ \frac{1}{2 \pi} \int_{-\pi}^{\pi} G(e^{j \omega})^H G(e^{j \omega)} d \omega \right]} = \sqrt{\frac{1}{2 \pi} \int_{-\pi}^{\pi} tr[G(e^{j \omega})^H G(e^{j \omega)}] d \omega}. \]

Weiters gilt durch den Satz von Parseval

\[ \sum_{k=-\infty}^{\infty} g_k^{H}g_k = \frac{1}{2 \pi} \int_{-\pi}^{\pi} G(e^{j \omega})^H G(e^{j \omega)} d \omega \]

der Zusammenhang

\[ \Vert G \Vert_2 = \Vert g \Vert_2 = \Vert tr\left[ \sum_{k=-\infty}^{\infty} g_k^{H}g_k \right] \Vert_2 = \sqrt{tr\left[ g(0)^2+g(1)^2+\dots \right]} \]

wobei

\[ g = [D,CB,CAB,CA^2B,\dots] \]

die Impuls-Antwort des Systems darstellt.

Im Zeitbereich wird die \(\mathcal{H}_2\)-Norm einer Übertragungsfunktion unter der Annahme berechnet, dass die stabile Übertragungsfunktion \(G(z)\) eine Zustandsraumdarstellung hat

\[ \begin{split} x_{k+1} &= Ax_{k} + Bu_{k} \\ y_{k} &= Cx_{k} + Du_{k} \end{split} \]

so dass

\[ G(z) = \frac{Y(z)}{U(z)} = C(zI-A)^{-1}B+D := \left[ \begin{array}{c|c} A & B\\ \hline C & D \end{array} \right] \]

gilt.

Anhand der obigen Definitionen lässt sich zeigen, dass die \(\mathcal{H}_2\)-Norm eines zeitdiskreten Systems äquivalent ist zu

\[ \Vert G \Vert_2 = \sqrt{tr(CPC^T+DD^T)} \]

wobei die Matrix \(P \geq 0\) durch das Lösen

\[APA^T-P+BB^T = 0\]

einer diskreten Lyapunov-Gleichung berechnet wird. Wie wir diese Gleichung lösen können, wissen wir aus dem letzten Blogeintrag.

Alternativ lässt sich die Lösung auch mit der Hilfe einer Reihe veranschaulichen, wenn auch nicht in einer geschlossenen Form. Dazu stellen wir die Lyapunov-Gleichung

\[P = APA^T+BB^T\]

um und führen eine Rekursion durch:

\[ \begin{split} \rightarrow P(1) &= A(APA^T+BB^T)A^T+BB^T = A^2P(A^T)^2 + ABB^TA^T +BB^T \\ \rightarrow P(2) &= A^3P(A^T)^3+ A^2BB^T(A^T)^2 + ABB^TA^T +BB^T \\ \rightarrow P(3) &= \dots \\ & \vdots \end{split} \]

\(P\) ist also gleich der konvergenten Reihe von Matrizen

\[ P = BB^T + ABB^TA^T + A^2BB^T(A^T)^2 + \dots . \]

Betrachten wir nochmals \(\mathcal{H}_2\)-Norm unter der Verwendung der Impulsanwort

\[ \begin{split} \Vert G \Vert_2^2 &= \Vert g \Vert_2^2 \\ &= tr[D^2 + (CB)^2 + (CAB)^2 + (CA^2B)^2 + \dots] \\ &= tr[D^2 + (CBB^TC^T) + (CABB^TA^TC^T) + (CA^2BB^T(A^T)^2C^T) + \dots] \\ &= tr[D^2 + C(BB^T + BB^TA^T + A^2BB^T(A^T)^2 + \dots)C^T] \\ &= tr[D^2 + CPC^T] \end{split}. \]

Die Impulsantwort kann so umgeschrieben werden, dass die Reihe für die Matrix \(P\) auftaucht. Ersetzen wir diese Reihe in der Impulsantwort durch \(P\), ergibt sich die Berechnungsvorschrift für die \(\mathcal{H}_2\)-Norm.

Tipp

Zusammengefasst lässt sich ein Zustandsraumverfahren wie folgt berechnen:

Berechne eine Zustandsrealisierung aus der Übertrangungsfunktion

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

Löse die Lyapunov Gleichung, (dabei muss gelten \(P \geq 0\))

\[ APA^T-P+BB^T = 0 \]

Berechne die \(\mathcal{H}_2\)-Norm

\[ \Vert G \Vert_2 = \sqrt{tr(CPC^T+DD^T)} \]

import numpy as np
import scipy.linalg as linalg

import warnings
warnings.filterwarnings("error")

def h2norm_d(A,B,C,D):
    """ naive implementation of the h2 system norm of linear time-discrete time-invariant systems (d-LTI)
        (there might be better ways to do that, => check papers)
    """
    try:
        P = linalg.solve_discrete_lyapunov(A.T, B@B.T)
        if np.all(np.real(linalg.eigvals(P)) < 0):
            H2 = np.inf
        else:
            H2 = np.sqrt(np.trace(C@P@C.T+D@D.T))
    except RuntimeWarning:
        H2 = np.inf
    return H2

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)
Dd = np.zeros((2,1))

h2norm_d(Ad_stable,Bd,Cd,Dd)

[ True  True]

1.1732382943111324

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

h2norm_d(Ad_unstable,Bd,Cd,Dd)

[False False]

inf

Fazit

Für zeitdiskrete Systeme ist die \(\mathcal{H}_2\)-Norm endlich, wenn das LTI-System asymptotisch stabil ist. Daraus folgt, dass für instabile Systeme die \(\mathcal{H}_2\)-Norm unendlich ist.

Eine deterministische Interpretation der \(\mathcal{H}_2\)-Norm ist, dass sie die Energie der Impulsantwort des LTI-Systems misst.

Referenzen

Optimal Sampled-Data Control Systems (Tongwen Chen und Bruce Francis, 1994)
Essentials of Robust Control (Kemin, Zhoun, John C. Doyle, 1997)
Policy Optimization for H2 Linear Control with H∞ Robustness Guarantee: Implicit Regularization and Global Convergence (Kaiqing Zhang Bin Hu Tamer Basar, arxiv, 2019)