1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA -
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution. The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt
function [x1,x2]=riccati(a,b,c,dom,typ)
//[x1,[x2]]=riccati(a,b,c,dom,[typ]) is a Riccati solver
// x=x1/x2 solves:
// a'*x+x*a-x*b*x+c=0 (continuous time case)
// a'*x*a-(a'*x*b1/(b2+b1'*x*b1))*(b1'*x*a)+c-x with b=b1/b2*b1'
// (discrete time case)
// If called with LHS=1 (one output argument) riccati returns x.
//
// -- a,b,c real matrices nxn, b and c symmetric.
// -- dom = 'c' or 'd' for the time domain (continuous or discrete)
// -- typ = 'eigen' --->block diagonalization
// typ = 'schur' --->schur method
// See also ric_desc
//!
[lhs,rhs]=argn(0),
if rhs==4 then typ="eigen",end,
ham=[a -b;-c -a'],
[n,n]=size(a),
if part(dom,1)=="c" then
select typ,
case "schur" then
[s,d]=schur(ham,"c"),
if d<>n then
error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",d, n))
end
s=s(:,1:n),
case "eigen" then
[hb,u1]=bdiag(ham),
[u2,d]=schur(hb,"c"),
u=u1*u2,
if d<>n then
error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",d, n))
end
s=u(:,1:n),
end,
else
aa=[eye(n,n) b;0*ones(n,n) a'],bb=[a 0*ones(n,n);-c eye(n,n)],
[bs,as,s,n1]=schur(bb,aa,"d");
if n1<>n then
error(msprintf(gettext("%s: Wrong dimension (%d) of stable subspace: %d expected.\n"),"riccati",n1, n))
end
s=s(:,1:n1);
end,
x1=s(n+1:2*n,:),x2=s(1:n,:),
if lhs==1 then x1=x1/x2,end
endfunction
|
