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 {ℝ}^{m\times n}$ gesucht, welche die Gleichung

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

erfüllt. Die anderen, vorgegebenen Matrizen haben die dazu passenden Dimensionen $C,B^T\in {ℝ}^{m\times n}$, $A\in {ℝ}^{n\times n}$, $D\in {ℝ}^{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:=\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\text{ mit }\left(\begin{array}{cc}{I}_n& 0\\ X& I_m\end{array}\right)$

auf obere Block-Dreieckform transformiert werden, bekommt man

$\left(\begin{array}{cc}{I}_n& 0\\ -X& I_m\end{array}\right)\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{cc}{I}_n& 0\\ X& I_m\end{array}\right)=\left(\begin{array}{cc}A+BX& B\\ 0& D-XB\end{array}\right)=:\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 {ℝ}^{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 $\left(\begin{array}{c}{I}_n\\ X\end{array}\right)$ 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

$\left(\begin{array}{cc}A& B\\ C& D\end{array}\right)\left(\begin{array}{c}{I}_n\\ X\end{array}\right)=\left(\begin{array}{c}{I}_n\\ X\end{array}\right)(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^{\prime }(t)=Ay(t)+Su(t)$ für einen Zustand $y(t)\in {ℝ}^n$ mit konstanten Koeffizienten $A\in {ℝ}^{n\times n}$, $S\in {ℝ}^{n\times p}$ soll diejenige optimale Steuerung $u(t)\in {ℝ}^p$ bestimmt werden, welche bei unendlichem Zeithorizont das Funktional

$\underset{0}{\overset{\mathrm{\infty }}{\int }}(y(t{)}^{T}Qy(t)+u(t{)}^{T}Ru(t))dt$

minimiert. Darin ist $R\in {ℝ}^{p\times p}$ symmetrisch und positiv definit, $Q\in {ℝ}^{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-XS{R}^{-1}{S}^{T}X=0$

ist, für welche die Matrix $A-SK=A-S{R}^{-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=-S{R}^{-1}{S}^T=B^T$. Die hierzu gehörige Blockmatrix

$L=\left(\begin{array}{cc}A& -S{R}^{-1}{S}^T\\ -Q& -A^T\end{array}\right)$

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 {ℝ}^{m\times n}$ ist die lineare Abbildung

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

Mit einer aktuellen Näherung $X_k\in {ℝ}^{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+B{X}_k)+(X_{k}B-D)H_k=C+D{X}_k-X_{k}A-X_{k}B{X}_k,k\ge 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+B{X}_k$ und $D-X_{k}B$ keine gemeinsamen Eigenwerte besitzen, z. B. wenn die Realteile aller Eigenwerte von $A+B{X}_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 {ℝ}^{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,\dots ,$

und man kann zeigen, dass diese Signum-Iteration immer und quadratisch konvergiert, sofern die Startmatrix $V_0:=M\in {ℝ}^{N\times N}$ keine rein imaginären Eigenwerte (einschließlich null) besitzt. Alle Matrizen $V_k,k\ge 0,$ kommutieren miteinander und besitzen daher die gleiche Jordan-Basis, und dies gilt auch für die Grenzwert-Matrix $S(M):=\underset{k\to \mathrm{\infty }}{\lim }{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=\left(\begin{array}{cc}0& I\\ C& 0\end{array}\right)=V_0,V_k=\left(\begin{array}{cc}0& R_k\\ P_k& 0\end{array}\right),k\ge 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}),R_{k+1}:=\frac{1}{2}(R_k+P_k^{-1}),k=0,1,\dots .$

Falls $C$ keine reellen und nicht-positiven Eigenwerte besitzt, konvergiert die Iteration gegen die eindeutige Wurzel $\underset{k\to \mathrm{\infty }}{\lim }{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)\left(\begin{array}{cc}{I}_n& 0\\ X& I_m\end{array}\right)=\left(\begin{array}{cc}{I}_n& 0\\ X& I_m\end{array}\right)\left(\begin{array}{cc}-I_n& G\\ 0& I_m\end{array}\right)$

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

$(S(M)+I_N)\left(\begin{array}{c}{I}_n\\ X\end{array}\right)=0,$

dass die Matrix $\left(\begin{array}{c}{I}_n\\ X\end{array}\right)$ 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

$\left(\begin{array}{c}{G}_{12}\\ G_{22}\end{array}\right)X=-\left(\begin{array}{c}{G}_{11}\\ G_{21}\end{array}\right)\text{ mit }S(M)+I_N=\left(\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right).$

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

Beispiel

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

$A=\left(\begin{array}{cc}-\frac{2}{3}& -2\\ -1& -\frac{8}{3}\end{array}\right),S=\left(\begin{array}{c}1\\ \frac{1}{2}\end{array}\right),Q=\left(\begin{array}{cc}1& \frac{5}{3}\\ \frac{5}{3}& \frac{20}{3}\end{array}\right),$

sowie $R=1\in {ℝ}^{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=\left(\begin{array}{cccc}-\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{array}\right)$

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)=\left(\begin{array}{cccc}\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{array}\right).$

Tatsächlich kann man direkt verifizieren, dass $S(L)$ involutorisch ist, $S(L{)}^2=I$, und mit $L$ kommutiert, $LS(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=\left(\begin{array}{cc}1& 0\\ \frac{3}{2}& 1\\ 0& -1\\ 2& 2\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1\\ \frac{3}{2}& -1\\ -1& 2\end{array}\right)\left(\begin{array}{cc}1& 0\\ \frac{3}{2}& 1\end{array}\right)=\left(\begin{array}{c}{I}_2\\ X\end{array}\right)\left(\begin{array}{cc}1& 0\\ \frac{3}{2}& 1\end{array}\right).$

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=\left(\begin{array}{cc}\frac{3}{2}& -1\\ -1& 2\end{array}\right),\text{ und es gilt }A+BX=A-S{R}^{-1}{S}^{T}X=\left(\begin{array}{cc}-\frac{5}{3}& -2\\ -\frac{3}{2}& -\frac{8}{3}\end{array}\right).$

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\ge 0,$ ist quadratisch, man kann dies direkt an den Eigenwerten der Matrizen $V_k$ ablesen. Diese lauten:

$\begin{array}{cc}k=& \text{Eigenwerte }{V}_k\\ 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 $\Vert V_6^2-I_4\Vert \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=\left(\begin{array}{cc}{G}_{11}& G_{12}\\ G_{21}& G_{22}\end{array}\right)$

bekommt man mit Hilfe der reduzierten QR-Zerlegung

$\left(\begin{array}{c}{G}_{12}\\ G_{22}\end{array}\right)=\hat{Q}\cdot \hat{R}\approx \left(\begin{array}{cc}-0.48245& 0.43832\\ 0.04444& -0.13345\\ 0.66238& -0.36922\\ 0.57138& 0.80854\end{array}\right)\left(\begin{array}{cc}1.39805& 1.07147\\ 0& 1.32121\end{array}\right)$

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

$\tilde{X}=-{\hat{R}}^{-1}{\hat{Q}}^T\left(\begin{array}{c}{G}_{11}\\ G_{21}\end{array}\right)=\left(\begin{array}{cc}1.499999999& -0.9999999997\\ -1.000000000& 2.000000001\end{array}\right).$

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