Matrix Riccati equation

As matrix Riccati-equations or algebraic Riccati equations is called a type of non-linear equations for matrices which, roughly speaking, traced in dimension 1 on an algebraic, quadratic equation can be. Hence the name of the problem based on the corresponding Riccati differential equation . In the case of general dimensions , a matrix is sought in a very general form of the Matrix-Riccati equation that contains the equation ${\displaystyle m,n\in \mathbb {N} }$ ${\displaystyle X\in \mathbb {R} ^{m\times n}}$

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

Fulfills. The other, given matrices have the matching dimensions , , . A special case of this equation is which solution has the square root of a matrix , if there are any. ${\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}}$ ${\displaystyle X^{2}=C}$ ${\displaystyle X=C^{1/2}}$

Meaning of the Riccati equation

In addition to the square root, Matrix Riccati equations occur with other important problems.

Eigenvalue problem, invariant subspaces

Should the block matrix ${\displaystyle (m+n)\times (m+1)}$

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

are transformed to the upper block triangular shape, one gets

${\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}},}$

if the solution of the above Riccati equation is, then the lower left block in the transformed matrix disappears . The two identity matrices, the dimension is indicated as an index . The multiplication of the 3 matrices actually represents a similarity transformation, since the left and right factors are inverse to each other. Therefore, the eigenvalues ​​of the overall matrix result from the union of the eigenvalues ​​of the two main diagonal blocks and of . In addition, the first columns of the transformation matrix form a basis for the associated invariant subspace (sum of eigenspaces ) of , from which the eigenvectors can be determined if necessary. So it applies ${\displaystyle X}$ ${\displaystyle C+DX-XA-XBX=0}$ ${\displaystyle I_{k}\in \mathbb {R} ^{k\times k}}$ ${\displaystyle M}$ ${\displaystyle A+BX}$ ${\displaystyle D-XB}$ ${\displaystyle {\hat {M}}}$ ${\displaystyle n}$ ${\displaystyle {\begin{pmatrix}I_{n}\\X\end{pmatrix}}}$ ${\displaystyle A+BX}$ ${\displaystyle M}$

${\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).}$

This property is used z. B. in the improvement of eigenvector bases: if a block triangle matrix emerged from interference, is small and, under suitable conditions, too . Then the block triangle shape can be restored in the manner indicated ([Stewart]). ${\displaystyle M}$ ${\displaystyle C}$ ${\displaystyle X}$

Continuous, optimal control

In the case of a linear system of differential equations for a state with constant coefficients , that optimal control should be determined which, given an infinite time horizon, the functional ${\displaystyle y'(t)=Ay(t)+Su(t)}$ ${\displaystyle y(t)\in \mathbb {R} ^{n}}$ ${\displaystyle A\in \mathbb {R} ^{n\times n}}$ ${\displaystyle S\in \mathbb {R} ^{n\times p}}$ ${\displaystyle u(t)\in \mathbb {R} ^{p}}$

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

minimized. In it, symmetric and positive is definite, symmetric and positive is semi-definite. Using feedback control , the optimum at an infinite time horizon is given by , where the ( maximum ) symmetrical solution of the Riccati equation ${\displaystyle R\in \mathbb {R} ^{p\times p}}$ ${\displaystyle Q\in \mathbb {R} ^{n\times n}}$ ${\displaystyle u(t)=-Ky(t)}$ ${\displaystyle u(t)=-R^{-1}S^{T}Xy(t)}$ ${\displaystyle X=X^{T}}$

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

is, for which the matrix is asymptotically stable with all eigenvalues in the left complex half-plane. For more background, please refer to the article QL controller . This equation is therefore a special case of the equation from the introduction of , , , . The associated block matrix ${\displaystyle A-SK=A-SR^{-1}S^{T}X}$ ${\displaystyle m=n}$ ${\displaystyle C=-Q}$ ${\displaystyle D=-A^{T}}$ ${\displaystyle B=-SR^{-1}S^{T}=B^{T}}$

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

is a Hamiltonian matrix , since and here are symmetrical. With this matrix , each eigenvalue also appears as an eigenvalue. ${\displaystyle B}$ ${\displaystyle C}$ ${\displaystyle L}$ ${\displaystyle \lambda }$ ${\displaystyle -\lambda }$

Numerical solution of Riccati equations

Newton's method

Since the Matrix-Riccati equation is an algebraic equation of degree 2 for the unknowns in the matrix , Newton's method can of course also be used for the solution . The derivation of the mapping at the point is the linear mapping ${\displaystyle m\cdot n}$ ${\displaystyle X}$ ${\displaystyle X\mapsto XBX+XA-DX-C}$ ${\displaystyle X\in \mathbb {R} ^{m\times n}}$

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

With a current approximation , you get the increment for an improved approximation from the linear system of equations ${\displaystyle X_{k}\in \mathbb {R} ^{m\times n}}$ ${\displaystyle H_{k}=X_{k+1}-X_{k}}$

${\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,}$

where on the right side, as usual, is the negative residual of the Riccati equation. The whole thing represents a Sylvester equation , in the accompanying article numerical methods for its solution are discussed. This linear equation can be solved uniquely if the two matrices and have no common eigenvalues, e.g. B. if the real parts of all eigenvalues ​​are from above and those from below a suitable value (approximately zero). ${\displaystyle A+BX_{k}}$ ${\displaystyle D-X_{k}B}$ ${\displaystyle A+BX_{k}}$ ${\displaystyle D-X_{k}B}$

Solution with the Signum iteration

Involutive matrices are solutions to the simple Riccati equation . The Newton iteration for this particular equation is also very simple, ${\displaystyle V\in \mathbb {R} ^{N\times N}}$ ${\displaystyle V^{2}=I}$

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

and one can show that this Signum iteration always and quadratically converges as long as the starting matrix has no purely imaginary eigenvalues ​​(including zero). All matrices commute with one another and therefore have the same Jordan basis, and this also applies to the limit value matrix . The associated eigenvalues ​​of the converge to or if the real part in the eigenvalue was positive or negative. Therefore has only the two eigenvalues and is called the sign function of , so it is an involution with . Since the eigenvalues ​​of are known, bases for the invariant subspaces are obtained or , by determining bases for the kernels of and , for example with the QR decomposition . These are then also bases for the invariant subspaces of the output matrix for the eigenvalues ​​with a positive or negative real part. ${\displaystyle V_{0}:=M\in \mathbb {R} ^{N\times N}}$ ${\displaystyle V_{k},\,k\geq 0,}$ ${\displaystyle S(M):=\lim _{k\to \infty }V_{k}}$ ${\displaystyle V_{k}}$ ${\displaystyle 1}$ ${\displaystyle -1}$ ${\displaystyle V_{0}=M}$ ${\displaystyle S(M)}$ ${\displaystyle \pm 1}$ ${\displaystyle M}$ ${\displaystyle S(M)}$ ${\displaystyle S^{2}=I}$ ${\displaystyle S(M)}$ ${\displaystyle +1}$ ${\displaystyle -1}$ ${\displaystyle S(M)-I}$ ${\displaystyle S(M)-I}$ ${\displaystyle M}$

This background can be used to solve the original Riccati equation, if the structure of or the number of eigenvalues ​​with positive and negative real parts is clear. That goes for the square root and the control problem. ${\displaystyle N=m+n}$ ${\displaystyle M}$ ${\displaystyle H}$

For the square root the main diagonal blocks vanish from and this is also the case with them , so be ${\displaystyle M}$ ${\displaystyle V_{k}}$

${\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,}$

with . The iteration for is then for the individual blocks ${\displaystyle N=2n}$ ${\displaystyle V_{k}}$

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

If it has no real and non-positive eigenvalues, the iteration converges to the unique root whose eigenvalue real parts are positive. ${\displaystyle C}$ ${\displaystyle \lim _{k\to \infty }P_{k}=C^{1/2}}$

The signum function can also be used in the general equation if the Riccati equation has a solution for which and are asymptotically stable, i.e. both only have eigenvalues ​​with a negative real part. Given this assumption, and and for the block matrix it follows that ${\displaystyle N=m+n}$ ${\displaystyle X}$ ${\displaystyle A+BX}$ ${\displaystyle XB-D}$ ${\displaystyle S(A+BX)=-I_{n}}$ ${\displaystyle S(D-XB)=+I_{m}}$ ${\displaystyle M}$

${\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}}}$

is with a suitable matrix . The first columns of this equation show with ${\displaystyle G}$ ${\displaystyle n}$

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

that the matrix is a special basis of the core of . To solve the Riccati equation , the matrices and their limit value must be calculated with the start matrix or with the optimal control . Then you get by dividing into blocks from the following system of equations ${\displaystyle {\begin{pmatrix}I_{n}\\X\end{pmatrix}}}$ ${\displaystyle S(M)+I_{N}}$ ${\displaystyle V_{0}:=M}$ ${\displaystyle V_{0}:=L}$ ${\displaystyle V_{k}}$ ${\displaystyle S(M)}$ ${\displaystyle X}$ ${\displaystyle S(M)+I_{N}}$

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

Here , , . ${\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}}$

example

In the application for optimal control, with and , ${\displaystyle n=2}$ ${\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}},}$

as well . One of the eigenvalues ​​of the matrix is positive, so the unregulated system is unstable. The more special form occurs here as the block matrix ${\displaystyle R=1\in \mathbb {R} ^{1\times 1}}$ ${\displaystyle -{\frac {5}{3}}\pm {\sqrt {3}}}$ ${\displaystyle A}$ ${\displaystyle u\equiv 0}$ ${\displaystyle M}$

${\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}}}$

on, it has the 4 eigenvalues , of which, as mentioned, 2 are actually positive and 2 are negative. In this example, the Signum function of can be calculated using its Jordan normal form , the result is ${\displaystyle \pm {\frac {13}{6}}\pm {\frac {1}{2}}{\sqrt {13}}}$ ${\displaystyle L}$

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

Indeed, one can verify directly that is in involution, and with commutated . A basis matrix of the kernel of , so with is given by ${\displaystyle S(L)}$ ${\displaystyle S(L)^{2}=I}$ ${\displaystyle L}$ ${\displaystyle L\,S(L)-S(L)\,L=0}$ ${\displaystyle Y}$ ${\displaystyle S(L)+I_{4}}$ ${\displaystyle (S(L)+I_{4})Y=0}$

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

By column-wise elimination in the first two lines of , an identity matrix was generated there and you can therefore read the solution to the Riccati equation in the lower block with ${\displaystyle Y}$ ${\displaystyle X}$

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

The controlled system matrix now has 2 negative eigenvalues ​​and the system is therefore asymptotically stable. ${\displaystyle A+BX}$

The calculation of the Jordan normal form can be avoided with the Signum iteration described in Section 2.2 . The convergence is quadratic; this can be read off directly from the eigenvalues ​​of the matrices . These are: ${\displaystyle V_{k}\to S(L),\,k\geq 0,}$ ${\displaystyle V_{k}}$

${\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}}}$

Indeed it is . If one uses the approximation instead of to calculate the solution and divides as described, ${\displaystyle \|V_{6}^{2}-I_{4}\|\approx 10^{-14}}$ ${\displaystyle X}$ ${\displaystyle S(L)}$ ${\displaystyle V_{6}}$ ${\displaystyle V_{6}+I_{4}}$

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

can be obtained with the help of the reduced QR decomposition

${\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}}}$

(Specification with low accuracy for reasons of space) the approximate solution

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

This approximation is obviously accurate to about 9 places.

literature

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