Matrix-Riccati-Gleichung

Als Matrix-Riccati-Gleichungen oder algebraische Riccati-Gleichungen wird ein Typ von nichtlinearen Gleichungen für Matrizen bezeichnet, die sich, grob gesagt, bei Dimension 1 auf eine algebraische, quadratische Gleichung zurückführen lassen. Daher kommt auch die Bezeichnung des Problems in Anlehnung an die entsprechende Riccati-Differentialgleichung. Bei allgemeinen Dimensionen $m,n\in \mathbb {N}$ ist in einer recht allgemeinen Form der Matrix-Riccati-Gleichung eine Matrix $X\in \mathbb {R} ^{m\times n}$ gesucht, welche die Gleichung

$XBX+XA-DX-C=0\in \mathbb {R} ^{m\times n}$

erfüllt. Die anderen, vorgegebenen Matrizen haben die dazu passenden Dimensionen $C,B^{T}\in \mathbb {R} ^{m\times n}$, $A\in \mathbb {R} ^{n\times n}$, $D\in \mathbb {R} ^{m\times m}$. Ein Spezialfall dieser Gleichung ist $X^{2}=C$, welche als Lösungen die Quadratwurzel einer Matrix $X=C^{1/2}$ hat, wenn solche existieren.

Bedeutung der Riccati-Gleichung

Außer bei der Quadratwurzel treten Matrix-Riccati-Gleichungen bei weiteren wichtigen Problemen auf.

Eigenwertproblem, invariante Unterräume

Soll die $(m+n)\times (m+1)$-Blockmatrix

$M:={\begin{pmatrix}A&B\\C&D\end{pmatrix}}{\mbox{ mit }}{\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}$

auf obere Block-Dreieckform transformiert werden, bekommt man

${\begin{pmatrix}I_{n}&0\\-X&I_{m}\end{pmatrix}}{\begin{pmatrix}A&B\\C&D\end{pmatrix}}{\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}={\begin{pmatrix}A+BX&B\\0&D-XB\end{pmatrix}}=:{\hat {M}},$

wenn $X$ Lösung der obigen Riccati-Gleichung ist, dann verschwindet der linke untere Block $C+DX-XA-XBX=0$ in der transformierten Matrix. Bei den beiden Einheitsmatrizen ist die Dimension als Index vermerkt, $I_{k}\in \mathbb {R} ^{k\times k}$. Die Multiplikation der 3 Matrizen stellt tatsächlich eine Ähnlichkeitstransformation dar, da der linke und der rechte Faktor zueinander invers sind. Daher ergeben sich die Eigenwerte der Gesamtmatrix $M$ aus der Vereinigung der Eigenwerte der beiden Hauptdiagonalblöcke $A+BX$ und $D-XB$ vom ${\hat {M}}$. Darüber hinaus bilden die ersten $n$ Spalten ${\begin{pmatrix}I_{n}\\X\end{pmatrix}}$ der Transformationsmatrix eine Basis für den zu $A+BX$ gehörigen invarianten Unterraum (Summe von Eigenräumen) von $M$, aus dem sich bei Bedarf die Eigenvektoren bestimmen lassen. Es gilt also

${\begin{pmatrix}A&B\\C&D\end{pmatrix}}{\begin{pmatrix}I_{n}\\X\end{pmatrix}}={\begin{pmatrix}I_{n}\\X\end{pmatrix}}(A+BX).$

Anwendung findet diese Eigenschaft z. B. bei der Nachbesserung von Eigenvektor-Basen: wenn $M$ durch Störungen aus einer Block-Dreieckmatrix hervorging, ist $C$ klein und unter geeigneten Voraussetzungen auch $X$. Dann kann die Block-Dreieckform in der angegebenen Weise wiederhergestellt werden ([Stewart]).

Kontinuierliche, optimale Steuerung

Bei einem linearen System von Differentialgleichungen $y'(t)=Ay(t)+Su(t)$ für einen Zustand $y(t)\in \mathbb {R} ^{n}$ mit konstanten Koeffizienten $A\in \mathbb {R} ^{n\times n}$, $S\in \mathbb {R} ^{n\times p}$ soll diejenige optimale Steuerung $u(t)\in \mathbb {R} ^{p}$ bestimmt werden, welche bei unendlichem Zeithorizont das Funktional

$\int _{0}^{\infty }{\big (}y(t)^{T}Qy(t)+u(t)^{T}Ru(t){\big )}dt$

minimiert. Darin ist $R\in \mathbb {R} ^{p\times p}$ symmetrisch und positiv definit, $Q\in \mathbb {R} ^{n\times n}$ symmetrisch und positiv semi-definit. Verwendet man eine Steuerung durch Rückkopplung $u(t)=-Ky(t)$, ist das Optimum bei unendlichem Zeithorizont gegeben durch $u(t)=-R^{-1}S^{T}Xy(t)$, wobei $X=X^{T}$ die (maximale) symmetrische Lösung der Riccati-Gleichung

$Q+XA+A^{T}X-XSR^{-1}S^{T}X=0$

ist, für welche die Matrix $A-SK=A-SR^{-1}S^{T}X$ asymptotisch stabil ist mit allen Eigenwerten in der linken komplexen Halbebene. Für mehr Hintergrund wird auf den Artikel LQ-Regler verwiesen. Diese Gleichung ist also ein Spezialfall der Gleichung aus der Einleitung mit $m=n$, $C=-Q$, $D=-A^{T}$, $B=-SR^{-1}S^{T}=B^{T}$. Die hierzu gehörige Blockmatrix

$L={\begin{pmatrix}A&-SR^{-1}S^{T}\\-Q&-A^{T}\end{pmatrix}}$

ist eine hamiltonsche Matrix, da $B$ und $C$ hier symmetrisch sind. Bei dieser Matrix $L$ tritt mit jedem Eigenwert $\lambda$ auch $-\lambda$ als Eigenwert auf.

Numerische Lösung von Riccati-Gleichungen

Newton-Verfahren

Da die Matrix-Riccati-Gleichung eine algebraische Gleichung vom Grad 2 für die $m\cdot n$ Unbekannten in der Matrix $X$ ist, kann zur Lösung natürlich auch das Newton-Verfahren eingesetzt werden. Die Ableitung der Abbildung $X\mapsto XBX+XA-DX-C$ an der Stelle $X\in \mathbb {R} ^{m\times n}$ ist die lineare Abbildung

$H\mapsto H(A+BX)+(XB-D)H{\text{ für }}H\in \mathbb {R} ^{m\times n}.$

Mit einer aktuellen Näherung $X_{k}\in \mathbb {R} ^{m\times n}$ bekommt man das Inkrement $H_{k}=X_{k+1}-X_{k}$ zu einer verbesserten Näherung also aus dem linearen Gleichungssystem

$H_{k}(A+BX_{k})+(X_{k}B-D)H_{k}=C+DX_{k}-X_{k}A-X_{k}BX_{k},\quad k\geq 0,$

wo auf der rechten Seite, wie gewohnt, das negative Residuum der Riccati-Gleichung steht. Das Ganze stellt eine Sylvester-Gleichung dar, im zugehörigen Artikel werden numerische Methoden zu ihrer Auflösung behandelt. Diese lineare Gleichung ist eindeutig lösbar, wenn die beiden Matrizen $A+BX_{k}$ und $D-X_{k}B$ keine gemeinsamen Eigenwerte besitzen, z. B. wenn die Realteile aller Eigenwerte von $A+BX_{k}$ oberhalb und die von $D-X_{k}B$ unterhalb eines geeigneten Wertes (etwa null) liegen.

Lösung mit der Signum-Iteration

Involutorische Matrizen $V\in \mathbb {R} ^{N\times N}$ sind Lösungen der einfachen Riccati-Gleichung $V^{2}=I$. Auch die Newton-Iteration für diese spezielle Gleichung ist sehr einfach,

$V_{k+1}:={\frac {1}{2}}(V_{k}+V_{k}^{-1}),\ k=0,1,\ldots ,$

und man kann zeigen, dass diese Signum-Iteration immer und quadratisch konvergiert, sofern die Startmatrix $V_{0}:=M\in \mathbb {R} ^{N\times N}$ keine rein imaginären Eigenwerte (einschließlich null) besitzt. Alle Matrizen $V_{k},\,k\geq 0,$ kommutieren miteinander und besitzen daher die gleiche Jordan-Basis, und dies gilt auch für die Grenzwert-Matrix $S(M):=\lim _{k\to \infty }V_{k}$. Die zugehörigen Eigenwerte der $V_{k}$ konvergieren gegen $1$ bzw. $-1$, wenn der Realteil im Eigenwert von $V_{0}=M$ positiv bzw. negativ war. Daher besitzt $S(M)$ nur die beiden Eigenwerte $\pm 1$ und wird als Signum-Funktion von $M$ bezeichnet, $S(M)$ ist also eine Involution mit $S^{2}=I$. Da die Eigenwerte von $S(M)$ bekannt sind, bekommt man Basen für die invarianten Unterräume zu $+1$ bzw. $-1$, indem man Basen für die Kerne von $S(M)-I$ bzw. $S(M)+I$ bestimmt, etwa mit der QR-Zerlegung. Diese sind dann auch Basen für die invarianten Unterräume der Ausgangsmatrix $M$ zu den Eigenwerten mit positivem bzw. negativem Realteil.

Diesen Hintergrund kann man mit $N=m+n$ zur Lösung der ursprünglichen Riccati-Gleichung verwenden, wenn aufgrund der Struktur von $M$ bzw. $H$ die Zahl der Eigenwerte mit positivem und negativem Realteil klar ist. Das gilt für die Quadratwurzel und das Steuerungsproblem.

Für die Quadratwurzel verschwinden die Hauptdiagonalblöcke von $M$ und auch bei den $V_{k}$ ist das so, sei also

$M={\begin{pmatrix}0&I\\C&0\end{pmatrix}}=V_{0},\quad V_{k}={\begin{pmatrix}0&R_{k}\\P_{k}&0\end{pmatrix}},\,k\geq 0,$

mit $N=2n$. Die Iteration für die $V_{k}$ lautet dann für die Einzelblöcke

$P_{k+1}={\frac {1}{2}}(P_{k}+R_{k}^{-1}),\quad R_{k+1}:={\frac {1}{2}}(R_{k}+P_{k}^{-1}),\ k=0,1,\ldots .$

Falls $C$ keine reellen und nicht-positiven Eigenwerte besitzt, konvergiert die Iteration gegen die eindeutige Wurzel $\lim _{k\to \infty }P_{k}=C^{1/2}$, deren Eigenwert-Realteile positiv sind.

Bei der allgemeinen Gleichung ist die Signum-Funktion mit $N=m+n$ einsetzbar, wenn die Riccati-Gleichung eine Lösung $X$ besitzt, für die $A+BX$ und $XB-D$ asymptotisch stabil sind, also beide nur Eigenwerte mit negativem Realteil besitzen. Unter dieser Voraussetzung ist $S(A+BX)=-I_{n}$ und $S(D-XB)=+I_{m}$ und für die Blockmatrix $M$ folgt, dass

$S(M){\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}={\begin{pmatrix}I_{n}&0\\X&I_{m}\end{pmatrix}}{\begin{pmatrix}-I_{n}&G\\0&I_{m}\end{pmatrix}}$

ist mit einer geeigneten Matrix $G$. Die ersten $n$ Spalten dieser Gleichung zeigen mit

${\big (}S(M)+I_{N}{\big )}{\begin{pmatrix}I_{n}\\X\end{pmatrix}}=0,$

dass die Matrix ${\begin{pmatrix}I_{n}\\X\end{pmatrix}}$ eine spezielle Basis des Kerns von $S(M)+I_{N}$ ist. Zur Lösung der Riccati-Gleichung sind also mit der Startmatrix $V_{0}:=M$, bzw. $V_{0}:=L$ bei der optimalen Steuerung, die Matrizen $V_{k}$ und ihr Grenzwert $S(M)$ zu berechnen. Danach bekommt man $X$ bei Aufteilung von $S(M)+I_{N}$ in Blöcke aus dem folgenden Gleichungssystem

${\begin{pmatrix}G_{12}\\G_{22}\end{pmatrix}}X=-{\begin{pmatrix}G_{11}\\G_{21}\end{pmatrix}}{\text{ mit }}S(M)+I_{N}={\begin{pmatrix}G_{11}&G_{12}\\G_{21}&G_{22}\end{pmatrix}}.$

Hier sind $G_{11}\in \mathbb {R} ^{n\times n}$, $G_{12},G_{21}^{T}\in \mathbb {R} ^{n\times m}$, $G_{22}\in \mathbb {R} ^{m\times m}$.

Beispiel

Bei der Anwendung zur optimalen Steuerung sei mit $n=2$ und $p=1$,

$A={\begin{pmatrix}-{\frac {2}{3}}&-2\\-1&-{\frac {8}{3}}\end{pmatrix}},\ S={\begin{pmatrix}1\\{\frac {1}{2}}\end{pmatrix}},\ Q={\begin{pmatrix}1&{\frac {5}{3}}\\{\frac {5}{3}}&{\frac {20}{3}}\end{pmatrix}},$

sowie $R=1\in \mathbb {R} ^{1\times 1}$. Von den Eigenwerten $-{\frac {5}{3}}\pm {\sqrt {3}}$ der Matrix $A$ ist einer positiv, das ungeregelte System mit $u\equiv 0$ ist also instabil. Als Blockmatrix $M$ tritt hier die speziellere Form

$L={\begin{pmatrix}-{\frac {2}{3}}&-2&-1&-{\frac {1}{2}}\\-1&-{\frac {8}{3}}&-{\frac {1}{2}}&-{\frac {1}{4}}\\-1&-{\frac {5}{3}}&{\frac {2}{3}}&1\\-{\frac {5}{3}}&-{\frac {20}{3}}&2&{\frac {8}{3}}\end{pmatrix}}$

auf, sie besitzt die 4 Eigenwerte $\pm {\frac {13}{6}}\pm {\frac {1}{2}}{\sqrt {13}}$, von denen, wie erwähnt, tatsächlich 2 positiv und 2 negativ sind. In diesem Beispiel lässt sich die Signum-Funktion von $L$ noch über deren Jordan-Normalform berechnen, das Ergebnis ist

$S(L)={\begin{pmatrix}{\frac {25}{338}}&-{\frac {135}{169}}&-{\frac {114}{169}}&{\frac {21}{338}}\\-{\frac {75}{338}}&-{\frac {115}{169}}&{\frac {21}{338}}&-{\frac {87}{676}}\\-{\frac {789}{676}}&{\frac {163}{338}}&-{\frac {25}{338}}&{\frac {75}{338}}\\{\frac {163}{338}}&-{\frac {433}{169}}&{\frac {135}{169}}&{\frac {115}{169}}\end{pmatrix}}.$

Tatsächlich kann man direkt verifizieren, dass $S(L)$ involutorisch ist, $S(L)^{2}=I$, und mit $L$ kommutiert, $L\,S(L)-S(L)\,L=0$. Eine Basismatrix $Y$ des Kerns von $S(L)+I_{4}$, also mit $(S(L)+I_{4})Y=0$ ist gegeben durch

$Y={\begin{pmatrix}1&0\\{\frac {3}{2}}&1\\0&-1\\2&2\end{pmatrix}}=\left({\begin{array}{cc}1&0\\0&1\\\hline {\frac {3}{2}}&-1\\-1&2\end{array}}\right){\begin{pmatrix}1&0\\{\frac {3}{2}}&1\end{pmatrix}}={\begin{pmatrix}I_{2}\\X\end{pmatrix}}{\begin{pmatrix}1&0\\{\frac {3}{2}}&1\end{pmatrix}}.$

Durch spaltenweise Elimination in den ersten beiden Zeilen von $Y$ wurde dort eine Einheitsmatrix erzeugt und man kann daher im unteren Block die Lösung $X$ der Riccati-Gleichung ablesen mit

$X={\begin{pmatrix}{\frac {3}{2}}&-1\\-1&2\end{pmatrix}},{\mbox{ und es gilt }}A+BX=A-SR^{-1}S^{T}X={\begin{pmatrix}-{\frac {5}{3}}&-2\\-{\frac {3}{2}}&-{\frac {8}{3}}\end{pmatrix}}.$

Die gesteuerte Systemmatrix $A+BX$ hat jetzt also 2 negative Eigenwerte und das System ist daher asymptotisch stabil.

Die Berechnung der Jordan-Normalform umgeht man mit der in Abschnitt 2.2 beschriebenen Signum-Iteration. Die Konvergenz $V_{k}\to S(L),\,k\geq 0,$ ist quadratisch, man kann dies direkt an den Eigenwerten der Matrizen $V_{k}$ ablesen. Diese lauten:

${\begin{array}{c|cc}k=&{\text{Eigenwerte }}V_{k}\\\hline 0&\pm 0.363891029&\pm 3.969442305\\1&\pm 1.555983235&\pm 2.110683431\\2&\pm 1.099331841&\pm 1.292231811\\3&\pm 1.004487641&\pm 1.033043387\\4&\pm 1.000010024&\pm 1.000528470\\5&\pm 1.000000000&\pm 1.000000140\\6&\pm 1.000000000&\pm 1.000000000\end{array}}$

Tatsächlich ist $\|V_{6}^{2}-I_{4}\|\approx 10^{-14}$. Setzt man zur Berechnung der Lösung $X$ an Stelle von $S(L)$ die Näherung $V_{6}$ ein und teilt $V_{6}+I_{4}$ auf wie beschrieben,

$V_{6}+I_{N}={\begin{pmatrix}G_{11}&G_{12}\\G_{21}&G_{22}\end{pmatrix}}$

bekommt man mit Hilfe der reduzierten QR-Zerlegung

${\begin{pmatrix}G_{12}\\G_{22}\end{pmatrix}}={\hat {Q}}\cdot {\hat {R}}\approx {\begin{pmatrix}-0.48245&0.43832\\0.04444&-0.13345\\0.66238&-0.36922\\0.57138&0.80854\end{pmatrix}}{\begin{pmatrix}1.39805&1.07147\\0&1.32121\end{pmatrix}}$

(Angabe aus Platzgründen mit geringer Genauigkeit) die Näherungslösung

${\tilde {X}}=-{\hat {R}}^{-1}{\hat {Q}}^{T}{\begin{pmatrix}G_{11}\\G_{21}\end{pmatrix}}={\begin{pmatrix}1.499999999&-0.9999999997\\-1.000000000&2.000000001\end{pmatrix}}.$

Diese Näherung ist offensichtlich auf ca. 9 Stellen genau.

Literatur

• G.W. Stewart, Error and perturbation bounds for subspaces associated with certain eigenvalue problems, SIAM Review 15, 727–764
• N.J. Higham, Functions of matrices: Theory and computation, SIAM, Philadelphia, 2008.
• J.D. Roberts, Linear model reduction and solution of the algebraic Riccati equation by use of the sign function, Intern. J. Control 32, 677–687