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
|
// 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 [p,s,t,l,rt,tt]=srfaur(h,f,g,r0,n,p,s,t,l)
//square-root algorithm for the algebraic Riccati equation.
//
// h,f,g : convenient matrices of the state-space model.
// r0 : E(yk*yk').
// n : number of iterations.
//
// p : estimate of the solution after n iterations.
// s,t,l : intermediate matrices for
// : successive iterations;
// rt,tt : gain matrices of the filter model after n iterations.
//
// : p,s,t,l may be given as input if more than one recursion
// : is desired (evaluation of intermediate values of p, e.g.).
//!
[lhs,rhs]=argn(0);
[d,m]=size(h);
if rhs==5,
//initializations
r0=.5*(r0+r0');
s=sqrtm(r0);
t=(g/s)';s=s';
l=-%i*t;
p=l'*l;
else,
if rhs<>9 then
error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"srfaur",9));
end;
end;
//recursion
for j=1:n,
a=[s t;l*h' l*f'];
[q,r]=qr(a);
s=r(1:d,1:d);
t=r(1:d,d+1:d+m);
l=r(d+1:2*d,d+1:d+m);
p=p+l'*l;
end;
//gain matrices of the filter.
rt=r0-h*p*h';
tt=(g-f*p*h')*inv(rt);
endfunction
|
