Das Kronecker Produkt und die Lyapunov-Gleichung

In diesem Blogeintrag wollen wir einen (naiven) Weg zeigen, wie man die zeitdiskrete Lyapunov-Gleichung

\[ A^TPA-P+Q=0 \]

löst.

Kronecker Produkt

Kronecker Produkt

Ist \(A\) eine \(m \times n\) und B eine \(p \times r\) Matrix, so ist das Kronecker-Produkt \(C = A \otimes B\) definiert als

\[ C = A \otimes B = (a_{ij} . B) = \begin{bmatrix} a_{11}B & \cdots & a_{11}B \\ \vdots & \ddots & \vdots \\ a_{m1}B & \cdots & a_{mn}B \end{bmatrix}. \]

Das Ergebnis ist also eine Matrix mit \(np\) Zeilen und \(nr\) Spalten, also \(np \times nr\).

Vektorisierung und “Matrixierung”

Gemeinsam mit dem Kronecker Produkt sind die Operation vec und mat hilfreich.

Gegeben sein eine quadratische Matrix

\[ H = \begin{bmatrix} h_{11} & h_{12} \\ h_{21} & h_{22} \end{bmatrix} \]

welche durch die Operation vec

\[ h = vec(H) = vec\left( \begin{bmatrix} h_{11} & h_{12} \\ h_{21} & h_{22} \end{bmatrix} \right) = \begin{bmatrix} h_{11} \\ h_{21} \\ h_{12} \\ h_{22} \end{bmatrix} \]

in einem Spaltenvektor überführt wird. Es werden einfach alle Spaltenvektoren der Matrix untereinander angeordnet womit der Spaltenvektor ensteht. In Python können wir das mit der Numpy-Methode flatten() erreichen.

Die Umgekehrung der Abbildung wollen wir mit mat bezeichnen. Es gilt also

\[ H = mat(vec(H)) = mat\left( \begin{bmatrix} h_{11} \\ h_{21} \\ h_{12} \\ h_{22} \end{bmatrix} \right) = \begin{bmatrix} h_{11} & h_{12} \\ h_{21} & h_{22} \end{bmatrix}. \]

Für uns relevant ist nun die Beziehung

\[ AXB=C \rightleftarrows (B^T\otimes A)vec(X)=vec(C). \]

Lyapunov Funktion

Gegeben sei die Lyapunov-Gleichung für zeitdiskrete Systeme

\[ A^TPA-P+Q=0 \]

mit den Matrizen \(A,P,Q\) der Dimension \(n \times n\), welche nun mit Hilfe des Kroneker Produktes gelöst werden kann. Dazu stellen wir die Gleichung um

\[ P-A^TPA=Q \]

und wenden die obige Rechregel

\[ vec(A^TPA) = (A^T \otimes A^T) vec(P) \]

auf das Produkt \(A^TPA\) an. Nun erhalten wir die vektorisierte Lyapunov-Gleichung

\[ vec(P) - (A^T \otimes A^T) vec(P) = vec(Q) \]

und können diese nach \(P\) auflösen

\[ vec(P) = \begin{bmatrix} I-A^T \otimes A^T \end{bmatrix}^{-1} vec(Q). \]

Die Matrix \(\begin{bmatrix} I-A^T \otimes A^T \end{bmatrix}\) ist eine \(n^2 \times n^2\) Matrix und der Vektor \(vec(Q)\) hat die Dimension \(n^2\). Die Invertierung besitzt dann die Laufzeitkomplexität \(\mathcal{O}(n^6)\).

Implementierung in Python mit scipy

Wir implementieren das nun in Python und verwenden technische Vektoren, also Zeilen- und Spaltenvektoren welche definitionsgemäß Matrizen sind, und somit zwei Dimension (Achsen) besitzen.

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

def vec(H):
    return np.expand_dims(H.flatten(),axis=1)

def mat(w):
    n = np.ceil(np.sqrt(len(w))).astype(int)
    return w.reshape(n,n)

Die Stabilität der Dynamikmatrix \(A\) kontrollieren wir anhand der Eigenwerte.

A = np.array([[0.5, 0.1],[0, 0.5]])
np.abs(linalg.eigvals(A)) < 1. # check if magnitued of all eigenvalues smaller than one

array([ True,  True])

Auch mit Hilfe der Lyapunov-Gleichung können wir die Stabilität überprüfen. (In scipy lautet die Lyapunov-Gleichung \(APA^T-P+Q=0\), weshalb die Transportierte von \(A\) eingesetzt wird. Mit help(modul.method) steht, wenn notwendig, die Hilfe zur Verfügung.)

#help(linalg.solve_discrete_lyapunov)

Q = np.eye(2)
P = linalg.solve_discrete_lyapunov(A.T,Q)

print(P)
print(np.real(linalg.eigvals(P)) > 0) # check if real part of all eigenvalues greater than zero

[[1.33333333 0.08888889]
 [0.08888889 1.36296296]]
[ True  True]

Alternativ wird hier die Lyapunov Gleichung mit dem obigen Algorithmus gelöst.

P_by_inv = mat(linalg.inv(np.eye(4)-np.kron(A.T,A.T))@vec(Q)) # inverse matrix
P_by_solve = mat(linalg.solve(np.eye(4)-np.kron(A.T,A.T),vec(Q))) # matlab A\b-notation 

print(np.allclose(P,P_by_inv))
print(np.allclose(P,P_by_solve))

True
True

Alle drei Methoden liefen dasselbe Ergebnis.

Fazit

Das Kroneker-Produkt ist oft hilfreich, Gleichungen so umzuformen, dass diese einfach lösbar werden. Diese Lösungen sind aber meist nicht die effizientesten Methoden hinsichtlich der Laufzeit (hier \(\mathcal{O}(n^6)\)). Das effiziente Lösen von Riccati-, Lyapunov- und Sylvester-Gleichungen ist aber alles andere als trivial. Meist gibt es auch nicht den einen Algorithmus, sonderen je nach Besetzung der Matrizen können verschieden Methoden zum Einsatz kommen. Als Anwender sollte man, wenn möglich, vorhandene Pakete verwenden. In Python stehen einige Methoden mit dem Modul scipy zur Verfügung. Im Falle der diskreten Lyapunov-Gleichung ist jedoch der obige Algorithmus in scipy wie folgt implementiert:

def _solve_discrete_lyapunov_direct(a, q):
    """
    Solves the discrete Lyapunov equation directly.
    This function is called by the `solve_discrete_lyapunov` function with
    `method=direct`. It is not supposed to be called directly.
    """

    lhs = kron(a, a.conj())
    lhs = np.eye(lhs.shape[0]) - lhs
    x = solve(lhs, q.flatten())

    return np.reshape(x, q.shape)

Referenzen

• https://de.wikipedia.org/wiki/Kronecker-Produkt
• https://en.wikipedia.org/wiki/Kronecker_product
• https://en.wikipedia.org/wiki/Lyapunov_equation
• https://github.com/scipy/scipy/blob/v1.7.1/scipy/linalg/_solvers.py#L234-L322