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 ${\displaystyle m,n\in \mathbb {N} }$ ist in einer recht allgemeinen Form der Matrix-Riccati-Gleichung eine Matrix ${\displaystyle X\in \mathbb {R} ^{m\times n}}$ gesucht, welche die Gleichung

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

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

Bedeutung der Riccati-Gleichung

Außer b​ei der Quadratwurzel treten Matrix-Riccati-Gleichungen b​ei weiteren wichtigen Problemen auf.

Eigenwertproblem, invariante Unterräume

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

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

auf o​bere Block-Dreieckform transformiert werden, bekommt man

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

${\displaystyle {\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 ${\displaystyle M}$ durch Störungen aus einer Block-Dreieckmatrix hervorging, ist ${\displaystyle C}$ klein und unter geeigneten Voraussetzungen auch ${\displaystyle X}$. Dann kann die Block-Dreieckform in der angegebenen Weise wiederhergestellt werden ([Stewart]).

Kontinuierliche, optimale Steuerung

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

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

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

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

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

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

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

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

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

${\displaystyle 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 ${\displaystyle A+BX_{k}}$ und ${\displaystyle D-X_{k}B}$ keine gemeinsamen Eigenwerte besitzen, z. B. wenn die Realteile aller Eigenwerte von ${\displaystyle A+BX_{k}}$ oberhalb und die von ${\displaystyle D-X_{k}B}$ unterhalb eines geeigneten Wertes (etwa null) liegen.

Lösung mit der Signum-Iteration

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

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

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

${\displaystyle 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 ${\displaystyle N=2n}$. Die Iteration für die ${\displaystyle V_{k}}$ lautet dann für die Einzelblöcke

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

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

${\displaystyle 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 ${\displaystyle G}$. Die ersten ${\displaystyle n}$ Spalten dieser Gleichung zeigen mit

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

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

${\displaystyle {\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 ${\displaystyle G_{11}\in \mathbb {R} ^{n\times n}}$, ${\displaystyle G_{12},G_{21}^{T}\in \mathbb {R} ^{n\times m}}$, ${\displaystyle G_{22}\in \mathbb {R} ^{m\times m}}$.

Beispiel

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

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

${\displaystyle 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 ${\displaystyle \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 ${\displaystyle L}$ noch über deren Jordan-Normalform berechnen, das Ergebnis ist

${\displaystyle 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 ${\displaystyle S(L)}$ involutorisch ist, ${\displaystyle S(L)^{2}=I}$, und mit ${\displaystyle L}$ kommutiert, ${\displaystyle L\,S(L)-S(L)\,L=0}$. Eine Basismatrix ${\displaystyle Y}$ des Kerns von ${\displaystyle S(L)+I_{4}}$, also mit ${\displaystyle (S(L)+I_{4})Y=0}$ ist gegeben durch

${\displaystyle 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 ${\displaystyle Y}$ wurde dort eine Einheitsmatrix erzeugt und man kann daher im unteren Block die Lösung ${\displaystyle X}$ der Riccati-Gleichung ablesen mit

${\displaystyle 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 ${\displaystyle 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 ${\displaystyle V_{k}\to S(L),\,k\geq 0,}$ ist quadratisch, man kann dies direkt an den Eigenwerten der Matrizen ${\displaystyle V_{k}}$ ablesen. Diese lauten:

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

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

bekommt m​an mit Hilfe d​er reduzierten QR-Zerlegung

${\displaystyle {\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 a​us Platzgründen m​it geringer Genauigkeit) d​ie Näherungslösung

${\displaystyle {\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 i​st offensichtlich a​uf 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