Matrix-Riccati-Gleichung

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 obere 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]).

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.

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.

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