Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Please login or create an account
Change language to: English - Français - Português - Русский
Scilabヘルプ >> CACSD > Matrix Computation > riccati

riccati

Solves the matricial Riccati equation (continuous | discrete time domain)

Syntax

X = riccati(H)
X = riccati(H, E)
X = riccati(A, B, C, dom)
X = riccati(A, B, C, dom, method)
[X1, X2, residual] = riccati(...)

Arguments

A, B, C
Square matrices of real numbers, of size n x n: Matricial coefficients of the equation.

H, E
Square matrices of real numbers, of size 2n x 2n: Hamiltonian matrices | pencil of the equation.

dom
Type / time domain of the Riccati equation: "c" or "continuous", or "d" or "discrete".

method
string : "eigen" for block diagonalization (default), or "schur" for Schur method.

X1, X2, X
square matrices of real numbers (X2 invertible). X is symetric when B and C are so.

residual
real number: norm1 of the actual Equation(X) residue.

Description

Continuous time domain

X = riccati(A, B, C, "c") solves the matricial Riccati equation:

A'*X + X*A - X*B*X + C = 0

X = riccati(H) where H = [A, -B ; -C, -A'] is the Hamiltonian matrix, does the same, using the 'eigen' method, but is more stable.

Discrete time domain

X = riccati(A, B, C, "d") solves the Riccati equation:

A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X = 0

with B = B1 / B2 * B1'.

X = riccati(H, E) does the same, where H = [A, zeros(n,n) ; -C, eye(n,n)] and E = [eye(n,n), B ; zeros(n,n), A'] define the Hamiltonian pencil (E,H).

Output options

[X1, X2] = riccati(..) provides X1 and X2, with X2 invertible, such that X = X1 / X2.

residual is the L1-norm of the actual equation's result. If X is an actual solution, residual should be 0. Most often its value is close to %eps*norm(X).

In the discrete case, sometimes B1 and B2 can't be retrieved from B. Then residual can't be assessed and is set to %nan.
When no solution is found, X=[] or X1=[] is returned, and a warning is yielded.

Examples

Continuous time domain

n = 3;
// [A, B, C] = (grand(n,n,"uin",-9,9), grand(n,n,"uin",-9,9), grand(n,n,"uin",-9,9))
A = [
  -62.   91.   57.
  -43.  -45.  -19.
   58.   83.  -62.
   ];
B = [
   75.  -31.  -10.
  -79.   70.   68.
  -72.  -5.    32.
   ];
C = [
  -56.   70.   58.
  -41.   54.   50.
   90.   2.   -40.
   ];

// With A, B, C
// ------------
X = riccati(A, B, C, "c")
clean(A'*X + X*A - X*B*X + C)
[x1, x2, res] = riccati(A,B,C, "c");
x = x1 / x2;
and(x==X)
res

// With the Hamiltonian
// --------------------
H = [A, -B; -C, -A'];
X = riccati(H)
clean(A'*X + X*A - X*B*X + C)
[x1, x2] = riccati(H);
x = x1 / x2;
and(x==X)
--> X = riccati(A, B, C, "c")
 X  =
  -0.1790367   0.4166084   0.2319152
  -0.4977093   0.7993495   0.3086213
   0.5595286   0.3202094  -0.103394

--> clean(A'*X + X*A - X*B*X + C)
 ans  =
   0.   0.   0.
   0.   0.   0.
   0.   0.   0.

--> [x1, x2, res] = riccati(A,B,C, "c");
--> x = x1 / x2;
--> and(x==X)
 ans  =
  T

--> res
 res  =
   3.340D-13


--> // With the Hamiltonian
--> // --------------------
--> H = [A, -B; -C, -A'];
--> X = riccati(H);
 X  =
  -0.1790367   0.4166084   0.2319152
  -0.4977093   0.7993495   0.3086213
   0.5595286   0.3202094  -0.103394

--> clean(A'*X + X*A - X*B*X + C)
 ans  =
   0.   0.   0.
   0.   0.   0.
   0.   0.   0.

--> [x1, x2] = riccati(H);
--> x = x1 / x2;
--> and(x==X)
 ans  =
  T

Discrete time domain

// Eq: A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X = 0
n = 4;
[A, B1, B2, C] = (rand(n,n), rand(n,n), rand(n,n), rand(n,n));
B  = B1 / B2 * B1';

X = riccati(A, B, C, 'd', 'schur')
if X <> [] then
    clean(A'*X*A - (A'*X*B1) / (B2 + B1'*X*B1) * (B1'*X*A) + C - X)
    [x1, x2, res] = riccati(A, B, C, 'd', 'schur');
    x = x1 / x2;
    and(x==X)
    res
else
    disp("Riccati: No solution found")
end

See also

  • ricc — リカッチ方程式
  • srfaur — 平方根アルゴリズム
  • lindquist — Lindquistのアルゴリズム

History

VersionDescription
6.1.0
  • riccati(H) and riccati(H,E) added.
  • residual output added.
  • When no solution is found, X=[] | X1=[] is now returned, without error.
Scilab Enterprises
Copyright (c) 2011-2017 (Scilab Enterprises)
Copyright (c) 1989-2012 (INRIA)
Copyright (c) 1989-2007 (ENPC)
with contributors
Last updated:
Tue Feb 25 08:53:20 CET 2020