Scilab 6.1.0
Change language to:
English - Français - 日本語 - Português -

See the recommended documentation of this function

Справка 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```

 Version Description 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.