icare

Implicit solver for continuous-time algebraic Riccati equations

Syntax

[X,K,L] = icare(A,B,Q,R,S,E,G)

[X,K,L,info] = icare(___)

[___] = icare(___,'noscaling')

[___] = icare(___,'anti')

Description

[X,K,L] = icare(A,B,Q,R,S,E,G) computes the unique stabilizing solution X, state-feedback gain K, and the closed-loop eigenvalues L of the following continuous-time algebraic Riccati equation.

$$A^{T}XE+E^{T}XA+E^{T}XGXE-\left(E^{T}XB+S\right)R^{-1}\left(B^{T}XE+S^T\right)+Q=0$$

The stabilizing solution X puts all the eigenvalues L in the left half-plane.

Algebraic Riccati equations play a key role in LQR/LQG control, H2- and H-infinity control, Kalman filtering, and spectral or co-prime factorizations.

[X,K,L,info] = icare(___) also returns a structure info which contains additional information about the solution to the continuous-time algebraic Riccati equation.

[___] = icare(___,'noscaling') turns off the built-in scaling and sets all entries of the scaling vectors info.Sx and info.Sr to 1. Turning off scaling speeds up computation but can be detrimental to accuracy when A,B,Q,R,S,E,G are poorly scaled.

[___] = icare(___,'anti') computes the anti-stabilizing solution X that puts all eigenvalues L in the right half-plane.

Examples

Solve Continuous-Time Algebraic Riccati Equation

To solve the algebraic Riccati equation ${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-\mathrm{XB}{\mathit{B}}^{\mathit{T}}\mathit{X}+\mathit{C}{\mathit{C}}^{\mathit{T}}=0$, consider the following matrices:

$$A=\left[\begin{array}{ccc}1& -2& 3\\ -4& 5& 6\\ 7& 8& 9\end{array}\right]B=\left[\begin{array}{c}5\\ 6\\ -7\end{array}\right]C=\left[\begin{array}{ccc}7& -8& 9\end{array}\right].$$

The least effort approach is to use $\mathit{G}=-\mathit{B}{\mathit{B}}^{\mathit{T}}$ and $\mathit{Q}={\mathit{C}}^{\mathit{T}}\mathit{C}$, and then find the solution using icare.

A = [-1,2,3;4,5,-6;7,-8,9];
B = [5;6;-7];
C = [7,-8,9];
G = -B*B';
Q = C'*C;
[X1,K1,L1] = icare(A,[],Q,[],[],[],G)

X1 = 3×3

   15.3201    4.2369   17.0090
    4.2369    2.6252    4.4123
   17.0090    4.4123   19.0374

K1 =

  0x3 empty double matrix

L1 = 3×1

   -3.2139
  -10.1191
  -76.9693

The above approach may lead to numerical inaccuracies when matrices B and C have large entries, since they are squared-up to compute the G and Q matrices. Due to limited numerical range, the computation may be less accurate or even fail. For example, if norm(B) is 1e6, then norm(G) is 1e12, and any eigenvalue within 1e-4 of the imaginary axis may be diagnosed as 'imaginary' due to numerical errors.

For greater numerical accuracy, re-write the algebraic Riccati equation in the following way:

${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}–[\mathrm{XB},{\mathit{C}}^{\mathit{T}}]*[I,0;0,-I][{\mathit{B}}^{\mathit{T}}\mathit{X};C]=0$.

The above equation is the standard form of ${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-(\mathrm{XB}+S){\mathit{R}}^{-1}(B^{\mathit{T}}\mathit{X}+S^{\mathit{T}})=0$,

where $\mathit{B}=\left[\mathit{B},\mathrm{0}\right],\mathit{S}=\left[0,{\mathit{C}}^{\mathit{T}}\right],\mathrm{and}\mathit{R}=\left[\mathit{I},0;0,-\mathit{I}\right]$.

Compute the solution using icare with the above values.

n = size(A,1);
m = size(B,2);
p = size(C,1);
R = blkdiag(eye(m),-eye(p));
BB = [B,zeros(n,p)];
S = [zeros(n,m),C'];
[X2,K2,L2,info] = icare(A,BB,0,R,S,[],[])

X2 = 3×3

   15.3201    4.2369   17.0090
    4.2369    2.6252    4.4123
   17.0090    4.4123   19.0374

K2 = 2×3

  -17.0406    6.0501  -21.7435
   -7.0000    8.0000   -9.0000

L2 = 3×1

   -3.2139
  -10.1191
  -76.9693

info = struct with fields:
        Sx: [3x1 double]
        Sr: [2x1 double]
         U: [3x3 double]
         V: [3x3 double]
         W: [2x3 double]
    Report: 0

Here, X2 is the unique stabilizing solution, K2 contains the state-feedback gain and L2 contains the closed-loop eigenvalues.

Anti-Stabilizing Solution of Continuous-Time Algebraic Riccati Equation

To find the anti-stabilizing solution of the continuous-time algebraic Riccati equation ${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-\mathrm{XB}{\mathit{B}}^{\mathit{T}}\mathit{X}+\mathit{C}{\mathit{C}}^{\mathit{T}}=0$, consider the following matrices:

$$A=\left[\begin{array}{ccc}1& -2& 3\\ -4& 5& 6\\ 7& 8& 9\end{array}\right]B=\left[\begin{array}{c}5\\ 6\\ -7\end{array}\right]C=\left[\begin{array}{ccc}7& -8& 9\end{array}\right].$$

For greater numerical accuracy, re-write the algebraic Riccati equation in the following way:

${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}–[\mathrm{XB},{\mathit{C}}^{\mathit{T}}]*[I,0;0,-I][{\mathit{B}}^{\mathit{T}}\mathit{X};C]=0$.

The above equation is the standard form of ${\mathit{A}}^{\mathit{T}}\mathit{X}+\mathrm{XA}-(\mathrm{XB}+S){\mathit{R}}^{-1}(B^{\mathit{T}}\mathit{X}+S^{\mathit{T}})=0$,

where $\mathit{B}=\left[\mathit{B},\mathrm{0}\right],\mathit{S}=\left[0,{\mathit{C}}^{\mathit{T}}\right],\mathrm{and}\mathit{R}=\left[\mathit{I},0;0,-\mathit{I}\right]$.

Compute the anti-stabilizing solution using the 'anti' option.

A = [-1,2,3;4,5,-6;7,-8,9];
B = [5;6;-7];
C = [7,-8,9];
n = size(A,1);
m = size(B,2);
p = size(C,1);
R = blkdiag(eye(m),-eye(p));
BB = [B,zeros(n,p)];
S = [zeros(n,m),C'];
[X,K,L] = icare(A,BB,0,R,S,[],[],'anti')

X = 3×3

  -18.0978   10.9090   -1.8466
   10.9090   -6.7195    1.4354
   -1.8466    1.4354   -0.9426

K = 2×3

  -12.1085    4.1803    5.9774
   -7.0000    8.0000   -9.0000

L = 3×1

   76.9693
   10.1191
    3.2139

Here, X is the unique anti-stabilizing solution, K contains the state-feedback gain and L contains the closed-loop eigenvalues.

Input Arguments

A,B,Q,R,S,E,G — Input matrices

matrices

Input matrices, specified as matrices.

The matrices Q, R and G must be Hermitian. A square matrix is Hermitian if it is equal to its complex conjugate transpose, that is, $$a_{i,j}={\overline{a}}_{j,i}$$.

For more information on Hermitian matrices, see ishermitian.

Matrices R and E must be invertible.

When matrices B, R, S, E and G are omitted or set to [], icare uses the following default values:

• B = 0

• R = I

• S = 0

• E = I

• G = 0

If the inputs Q, R and G are scalar-valued, icare interprets them as multiples of the identity matrix.

'noscaling' — Option to turn off built-in scaling

'noscaling'

Option to turn off built-in scaling, specified as 'noscaling'. When you turn off the built-in scaling, icare sets all entries in the scaling vectors info.Sx and info.Sr to 1. Turning off scaling speeds up computation but can be detrimental to accuracy when A,B,Q,R,S,E,G are poorly scaled.

'anti' — Option to compute the anti-stabilizing solution

'anti'

Option to compute the anti-stabilizing solution, specified as 'anti'. When you enable this option, icare computes the anti-stabilizing solution X that puts all eigenvalues of (A+G*X*E-B*K,E) in the right half-plane.

The unique stabilizing and anti-stabilizing are both needed to know the complete phase portrait of the Riccati differential equations.

Output Arguments

X — Unique solution to the Riccati equation

matrix

Unique solution to the continuous-time algebraic Riccati equation, returned as a matrix.

By default, X is the stabilizing solution of the continuous-time algebraic Riccati equation. When the 'anti' option is used, X is the anti-stabilizing solution.

icare returns [] for X when the associated Hamiltonian matrix has eigenvalues on the imaginary axis.

icare returns [] when pencil is singular, that is, [B;S;R] is rank deficient. (since R2024a)

K — State-feedback gain

matrix

State-feedback gain, returned as a matrix.

The state-feedback gain K is computed as:

$$K=R^{-1}\left(B^{T}XE+S^T\right).$$

icare returns [] for K when the associated Hamiltonian matrix has eigenvalues on the imaginary axis.

icare returns [] when pencil is singular, that is, [B;S;R] is rank deficient. (since R2024a)

L — Closed-loop eigenvalues

matrix

Closed-loop eigenvalues, returned as a matrix.

The closed-loop eigenvalues L is computed as:

$$L=eig\left(A+GXE-BK,E\right).$$

icare returns [] for X and K when the associated Hamiltonian matrix has eigenvalues on the imaginary axis. In other words, L is non-empty even when X and K are empty matrices.

icare returns L = NaN(n,1) when pencil is singular, that is, [B;S;R] is rank deficient. (since R2024a)

info — Information about the unique solution

structure

Information about the unique solution, returned as a structure with the following fields:

• Sx — Vector of values used to scale the states.

• Sr — Vector of values used to scale the R matrix.

• U, V and W — Vectors of values representing the basis of the stable invariant subspace of the associated scaled matrix pencil. For more information, see Algorithms.

• Report — A scalar with one of the following values:

  • 0 — The unique solution is accurate.

  • 1 — The solution accuracy is poor.

  • 2 — The solution is not finite.

  • 3 — No solution found since the Hamiltonian spectrum, denoted by [L;-L], has eigenvalues on the imaginary axis.

  • 4 — Pencil is singular ([B;S;R] is rank deficient) (since R2024a)

Limitations

• (A-sE,B) must be stabilizable, and E and R must be invertible for a finite stabilizing solution X to exist and be finite. While these conditions are not sufficient in general, they become sufficient when the following conditions are met:

  • $$\left[\begin{array}{cc}\begin{array}{c}Q\\ S^T\end{array}& \begin{array}{c}S\\ R\end{array}\end{array}\right]\ge 0$$

  • $$\left[\begin{array}{cc}A-B{R}^{-1}{S}^T& Q-S{R}^{-1}{S}^T\end{array}\right]$$ is detectible

Algorithms

Basis of the invariant subspace

icare works with the following pencil, and computes a basis [U;V;W] of the invariant subspace associated with the stable or anti-stable finite eigenvalues of this pencil.

$$M-sN=\left[\begin{array}{ccc}A& G& B\\ -Q& -A^T& -S\\ S^T& B^T& R\end{array}\right]-s\left[\begin{array}{ccc}E& 0& 0\\ 0& E^T& 0\\ 0& 0& 0\end{array}\right]$$

The data is automatically scaled to reduce the sensitivity of eigenvalues near the imaginary axis and increase separation between the stable and anti-stable invariant subspaces.

Relationship between the solution, the state-feedback gain, and the scaling vectors

The solution X and state-feedback gain K are related to the scaling vectors, and U,V,W by the following set of equations:

$$\begin{array}{l}X={\text{ D}}_{x}V{U}^{-1}{\text{ D}}_x{E}^{-1}\text{,}\\ K=-D_r{\text{ WU}}^{-1}{\text{ D}}_x,\end{array}$$

where,

$$\begin{array}{l}{D}_x=\mathrm{diag}\left(S_x\right),\\ D_r=\mathrm{diag}\left(S_r\right).\end{array}$$

Version History

Introduced in R2019a

R2024a: Improved singularity detection

The function now returns X = [], K = [], L = NaN(n,1), and info.Report = 4 when pencil is singular, that is, [B;S;R] is rank deficient.