X = Matrices.discreteRiccati(A, B, R, Q);
(X, alphaReal, alphaImag) = Matrices.discreteRiccati(A, B, R, Q, true);
Function discreteRiccati computes the solution X of the discrete-time algebraic Riccati equation
A'*X*A - X - A'*X*B*inv(R + B'*X*B)*B'*X*A + Q = 0
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. Using this method, A has also to be invertible.
These assumptions are not checked in this function !!!
The assumptions guarantee that the Hamiltonian matrix.
H = [A + G*T*Q, -G*T; -T*Q, T]
with
-T
T = A
and
-1
G = B*R *B'
has no eigenvalue on the unit circle 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
-1
A - B*(R + B'*X*B) *B'*X*A
If U is partitioned to
U = [U11, U12; U21, U22]
according to S, the solution X can be calculated by
X*U11 = U21.
[1] Laub, A.J.
A Schur Method for Solving Algebraic Riccati equations.
IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
A = [4.0 3.0]
-4.5, -3.5];
B = [ 1.0;
-1.0];
R = [1.0];
Q = [9.0, 6.0;
6.0, 4.0]
X = discreteRiccati(A, B, R, Q);
results in:
X = [14.5623, 9.7082;
9.7082, 6.4721];
function discreteRiccati extends Modelica.Icons.Function; import Modelica.Math.Matrices; input Real A[:, size(A, 1)] "Square matrix A in DARE"; input Real B[size(A, 1), :] "Matrix B in DARE"; input Real R[size(B, 2), size(B, 2)] = identity(size(B, 2)) "Matrix R in DARE"; input Real Q[size(A, 1), size(A, 1)] = identity(size(A, 1)) "Matrix Q in DARE"; input Boolean refine = false "True for subsequent refinement"; output Real X[size(A, 1), size(A, 2)] "orthogonal matrix of the Schur vectors associated to ordered rsf"; output Real alphaReal[2 * size(A, 1)] "Real part of eigenvalue=alphaReal+i*alphaImag"; output Real alphaImag[2 * size(A, 1)] "Imaginary part of eigenvalue=alphaReal+i*alphaImag"; end discreteRiccati;