Matrix-Riccati-Gleichung

Als Matrix-Riccati-Gleichungen oder algebraische Riccati-Gleichungen wird ein Typ von nichtlinearen Gleichungen für Matrizen bezeichnet, die sich 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 <math>m,n\in\N</math> ist in allgemeiner Form der Matrix-Riccati-Gleichung eine Matrix <math>X\in\R^{m\times n}</math> gesucht, welche die Gleichung

<math> XBX+XA-DX-C=0 \in\R^{m\times n}

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

Bedeutung der Riccati-Gleichung

Außer bei der Quadratwurzel treten Matrix-Riccati-Gleichungen bei den folgenden Problemen auf.

Eigenwertproblem, invariante Unterräume

Soll die <math>(m+n)\times(m+1)</math>-Blockmatrix

<math> M:=\begin{pmatrix}

A& B\\ C&D

\end{pmatrix} \mbox{ mit } \begin{pmatrix}

I_n&0\\ X&I_m

\end{pmatrix} </math> auf obere Block-Dreieckform transformiert werden, bekommt man

<math> \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, </math> wenn <math>X</math> Lösung der obigen Riccati-Gleichung ist, dann verschwindet der linke untere Block <math>C+DX-XA-XBX=0</math> in der transformierten Matrix. Bei den beiden Einheitsmatrizen ist die Dimension als Index vermerkt, <math>I_k\in\R^{k\times k}</math>. 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 <math>M</math> aus der Vereinigung der Eigenwerte der beiden Hauptdiagonalblöcke <math>A+BX</math> und <math>D-XB</math> vom <math>\hat M</math>. Darüber hinaus bilden die ersten <math>n</math> Spalten <math>\begin{pmatrix} I_n\\ X\end{pmatrix}</math> der Transformationsmatrix eine Basis für den zu <math>A+BX</math> gehörigen invarianten Unterraum (Summe von Eigenräumen) von <math>M</math>, aus dem sich bei Bedarf die Eigenvektoren bestimmen lassen. Es gilt also

<math> \begin{pmatrix}

A& B\\ C&D

\end{pmatrix}\begin{pmatrix} I_n\\ X\end{pmatrix} =\begin{pmatrix} I_n\\ X\end{pmatrix}(A+BX). </math> Anwendung findet diese Eigenschaft z. B. bei der Nachbesserung von Eigenvektor-Basen: wenn <math>M</math> durch Störungen aus einer Block-Dreieckmatrix hervorging, ist <math>C</math> klein und unter geeigneten Voraussetzungen auch <math>X</math>. Dann kann die Block-Dreieckform in der angegebenen Weise wiederhergestellt werden ([Stewart]).

Kontinuierliche, optimale Steuerung

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

<math> \int_0^{\infty}\big(y(t)^T Qy(t)+u(t)^T Ru(t)\big)dt

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

<math> Q+XA+A^T X-XSR^{-1}S^T X=0

</math> ist, für welche die Matrix <math>A-SK=A-SR^{-1}S^T X</math> 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 <math>m=n</math>, <math>C=-Q</math>, <math>D=-A^T</math>, <math>B=-SR^{-1}S^T=B^T</math>. Die hierzu gehörige Blockmatrix

<math> L=\begin{pmatrix}

A &-SR^{-1}S^T \\ -Q&-A^T

\end{pmatrix} </math> ist eine hamiltonsche Matrix, da <math>B</math> und <math>C</math> hier symmetrisch sind. Bei dieser Matrix <math>L</math> tritt mit jedem Eigenwert <math>\lambda</math> auch <math>-\lambda</math> 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 <math>m\cdot n</math> Unbekannten in der Matrix <math>X</math> ist, kann zur Lösung das Newton-Verfahren eingesetzt werden. Die Ableitung der Abbildung <math>X\mapsto XBX+XA-DX-C</math> an der Stelle <math>X\in\R^{m\times n}</math> ist die lineare Abbildung

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

</math> Mit einer aktuellen Näherung <math>X_k\in\R^{m\times n}</math> bekommt man das Inkrement <math>H_k=X_{k+1}-X_k</math> zu einer verbesserten Näherung also aus dem linearen Gleichungssystem

<math> H_k(A+BX_k)+(X_kB-D)H_k=C+DX_k-X_kA-X_kBX_k,\quad k\ge0,

</math> 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 <math>A+BX_k</math> und <math>D-X_kB</math> keine gemeinsamen Eigenwerte besitzen, z. B. wenn die Realteile aller Eigenwerte von <math>A+BX_k</math> oberhalb und die von <math>D-X_kB</math> unterhalb eines geeigneten Wertes (etwa null) liegen.

Lösung mit der Signum-Iteration

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

<math> V_{k+1}:=\frac12(V_k+V_k^{-1}),\ k=0,1,\ldots,

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

Diesen Hintergrund kann man mit <math>N=m+n</math> zur Lösung der ursprünglichen Riccati-Gleichung verwenden, wenn aufgrund der Struktur von <math>M</math> bzw. <math>H</math> 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 <math>M</math> und auch bei den <math>V_k</math> ist das so, sei also

<math> 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\ge0,

</math> mit <math>N=2n</math>. Die Iteration für die <math>V_k</math> lautet dann für die Einzelblöcke

<math> P_{k+1}=\frac12(P_k+R_k^{-1}),\quad R_{k+1}:=\frac12(R_k+P_k^{-1}),\ k=0,1,\ldots.

</math> Falls <math>C</math> keine reellen und nicht-positiven Eigenwerte besitzt, konvergiert die Iteration gegen die eindeutige Wurzel <math>\lim_{k\to\infty}P_k=C^{1/2}</math>, deren Eigenwert-Realteile positiv sind.

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

<math> 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} </math> ist mit einer geeigneten Matrix <math>G</math>. Die ersten <math>n</math> Spalten dieser Gleichung zeigen mit

<math> \big(S(M)+I_N\big)\begin{pmatrix}I_n\\ X\end{pmatrix}=0,

</math> dass die Matrix <math>\begin{pmatrix}I_n\\ X\end{pmatrix}</math> eine spezielle Basis des Kerns von <math>S(M)+I_N</math> ist. Zur Lösung der Riccati-Gleichung sind also mit der Startmatrix <math>V_0:=M</math>, bzw. <math>V_0:=L</math> bei der optimalen Steuerung, die Matrizen <math>V_k</math> und ihr Grenzwert <math>S(M)</math> zu berechnen. Danach bekommt man <math>X</math> bei Aufteilung von <math>S(M)+I_N</math> in Blöcke aus dem folgenden Gleichungssystem

<math> \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}.

</math> Hier sind <math>G_{11}\in\R^{n\times n}</math>, <math>G_{12},G_{21}^T\in\R^{n\times m}</math>, <math>G_{22}\in\R^{m\times m}</math>.

Beispiel

Bei der Anwendung zur optimalen Steuerung sei mit <math>n=2</math> und <math>p=1</math>,

<math> A=\begin{pmatrix} -\frac{2}3&-2\\ -1&-\frac{8}3\end{pmatrix},

\ S=\begin{pmatrix}1\\\frac12 \end{pmatrix},
\ Q=\begin{pmatrix} 1&\frac53\\\frac53&\frac{20}3 \end{pmatrix},

</math> sowie <math>R=1\in\R^{1\times 1}</math>. Von den Eigenwerten <math>-\frac53\pm\sqrt3</math> der Matrix <math>A</math> ist einer positiv, das ungeregelte System mit <math>u\equiv0</math> ist also instabil. Als Blockmatrix <math>M</math> tritt hier die speziellere Form

<math> L=\begin{pmatrix}

-\frac23& -2& -1& -\frac12\\ -1& -\frac83& -\frac12& -\frac14\\ -1& -\frac53& \frac23& 1\\ -\frac53& -\frac{20}3& 2& \frac83 \end{pmatrix} </math> auf, sie besitzt die 4 Eigenwerte <math>\pm\frac{13}6\pm\frac12\sqrt{13}</math>, von denen, wie erwähnt, tatsächlich 2 positiv und 2 negativ sind. In diesem Beispiel lässt sich die Signum-Funktion von <math>L</math> noch über deren Jordan-Normalform berechnen, das Ergebnis ist

<math> 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}.

</math> Tatsächlich kann man direkt verifizieren, dass <math>S(L)</math> involutorisch ist, <math>S(L)^2=I</math>, und mit <math>L</math> kommutiert, <math>L\,S(L)-S(L)\, L=0</math>. Eine Basismatrix <math>Y</math> des Kerns von <math>S(L)+I_4</math>, also mit <math>(S(L)+I_4)Y=0</math> ist gegeben durch

<math> Y=\begin{pmatrix} 1&0\\\frac32&1\\0&-1\\2&2\end{pmatrix}

 =\left(\begin{array}{cc} 1&0\\0&1\\\hline
 \frac32&-1\\-1&2\end{array}\right)
 \begin{pmatrix} 1&0\\\frac32&1\end{pmatrix}
 =\begin{pmatrix} I_2\\X \end{pmatrix}
 \begin{pmatrix} 1&0\\\frac32&1\end{pmatrix}.

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

<math> X=\begin{pmatrix}\frac32&-1\\-1&2\end{pmatrix},

 \mbox{ und es gilt }
 A+BX=A-SR^{-1}S^TX=
 \begin{pmatrix} -\frac53& -2\\ -\frac32& -\frac83 \end{pmatrix}.

</math> Die gesteuerte Systemmatrix <math>A+BX</math> 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 <math>V_k\to S(L),\,k\ge0,</math> ist quadratisch, man kann dies direkt an den Eigenwerten der Matrizen <math>V_k</math> ablesen. Diese lauten:

<math> \begin{array}{c|cc}

k=&\text{Eigenwerte }V_k\\\hline
0&\pm0.363891029&\pm3.969442305\\
1&\pm1.555983235&\pm2.110683431\\
2&\pm1.099331841&\pm1.292231811\\
3&\pm1.004487641&\pm1.033043387\\
4&\pm1.000010024&\pm1.000528470\\
5&\pm1.000000000&\pm1.000000140\\
6&\pm1.000000000&\pm1.000000000
\end{array}

</math> Tatsächlich ist <math>\|V_6^2-I_4\|\approx10^{-14}</math>. Setzt man zur Berechnung der Lösung <math>X</math> an Stelle von <math>S(L)</math> die Näherung <math>V_6</math> ein und teilt <math>V_6+I_4</math> auf wie beschrieben,

<math> V_6+I_N=\begin{pmatrix} G_{11}& G_{12}\\G_{21}&G_{22}\end{pmatrix}

</math> bekommt man mit Hilfe der reduzierten QR-Zerlegung

<math> \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}

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

<math> \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}.

</math> 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