Scilab Home page | Wiki | Bug tracker | Forge | Mailing list archives | ATOMS | File exchange
Scilab-Branch-6.1-GIT
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

History

 Версия Описание 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.