X = Matrices.continuousRiccati(A, B, R, Q); (X, alphaReal, alphaImag) = Matrices.continuousRiccati(A, B, R, Q, true);
Function continuousRiccati computes the solution X of the continuous-time algebraic Riccati equation
A'*X + X*A - X*G*X + Q = 0
with G = B*inv(R)*B'
using the Schur vector approach proposed by Laub [1].
It is assumed that Q is symmetric and positive semidefinite and R is symmetric, nonsingular and positive definite, (A,B) is stabilizable and (A,Q) is detectable.
These assumptions are not checked in this function !!
The assumptions guarantee that the Hamiltonian matrix
H = [A, -G; -Q, -A']
has no pure imaginary eigenvalue and can be put to an ordered real Schur form
U'*H*U = S = [S11, S12; 0, S22]
with orthogonal similarity transformation U. S is ordered in such a way, that S11 contains the n stable eigenvalues of the closed loop system with system matrix A - B*inv(R)*B'*X. If U is partitioned to
U = [U11, U12; U21, U22]
with dimensions according to S, the solution X is calculated by
X*U11 = U21.
With optional input refinement=true
a subsequent iterative refinement based on Newton's method with exact line search is applied.
See continuousRiccatiIterative
for more information.
[1] Laub, A.J. A Schur Method for Solving Algebraic Riccati equations. IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
A = [0.0, 1.0; 0.0, 0.0]; B = [0.0; 1.0]; R = [1]; Q = [1.0, 0.0; 0.0, 2.0]; X = continuousRiccati(A, B, R, Q); results in: X = [2.0, 1.0; 1.0, 2.0];
Matrices.Utilities.continuousRiccatiIterative, Matrices.discreteRiccati
function continuousRiccati extends Modelica.Icons.Function; import Modelica.Math.Matrices; input Real A[:, size(A, 1)] "Square matrix A in CARE"; input Real B[size(A, 1), :] "Matrix B in CARE"; input Real R[size(B, 2), size(B, 2)] = identity(size(B, 2)) "Matrix R in CARE"; input Real Q[size(A, 1), size(A, 1)] = identity(size(A, 1)) "Matrix Q in CARE"; input Boolean refine = false "True for subsequent refinement"; output Real X[size(A, 1), size(A, 2)] "Stabilizing solution of CARE"; output Real alphaReal[2*size(A, 1)] "Real parts of eigenvalue=alphaReal+i*alphaImag"; output Real alphaImag[2*size(A, 1)] "Imaginary parts of eigenvalue=alphaReal+i*alphaImag"; end continuousRiccati;