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
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
|
function [R,n,sval,rcnd]=findR(s,y,u,meth,alg,jobd,tol,printw)
R=[];n=[];sval=[];rcnd=[];
[nargout,nargin] = argn(0)
//FINDR Preprocesses the input-output data for estimating the matrices
//of a linear time-invariant dynamical system, using Cholesky or
//(fast) QR factorization and subspace identification techniques
//(MOESP or N4SID), and estimates the system order.
//
//[R,N] = FINDR(S,Y,U,METH,ALG,JOBD,TOL,PRINTW) returns the processed
//upper triangular factor R of the concatenated block-Hankel matrices
//built from the input-output data, and the order N of a discrete-time
//realization. The model structure is:
//
// x(k+1) = Ax(k) + Bu(k) + w(k), k >= 1,
// y(k) = Cx(k) + Du(k) + e(k).
//
//The vectors y(k) and u(k) are transposes of the k-th rows of Y and U,
//respectively.
//
//S is the number of block rows in the block-Hankel matrices.
//
//METH is an option for the method to use:
//METH = 1 : MOESP method with past inputs and outputs;
// = 2 : N4SID method.
//Default: METH = 1.
//
//ALG is an option for the algorithm to compute the triangular factor of
//the concatenated block-Hankel matrices built from the input-output data:
//ALG = 1 : Cholesky algorithm on the correlation matrix;
// = 2 : fast QR algorithm;
// = 3 : standard QR algorithm.
//Default: ALG = 1.
//
//JOBD is an option to specify if the matrices B and D should later
//be computed using the MOESP approach:
//JOBD = 1 : the matrices B and D should later be computed using
// the MOESP approach;
// = 2 : the matrices B and D should not be computed using
// the MOESP approach.
//Default: JOBD = 2.
//This parameter is not relevant for METH = 2.
//
//TOL is a vector of length 2 containing tolerances:
//TOL(1) is the tolerance for estimating the rank of matrices.
//If TOL(1) > 0, the given value of TOL(1) is used as a
//lower bound for the reciprocal condition number.
//Default: TOL(1) = prod(size(matrix))*epsilon_machine where
// epsilon_machine is the relative machine precision.
//TOL(2) is the tolerance for estimating the system order.
//If TOL(2) >= 0, the estimate is indicated by the index of
//the last singular value greater than or equal to TOL(2).
//(Singular values less than TOL(2) are considered as zero.)
//When TOL(2) = 0, then S*epsilon_machine*sval(1) is used instead
//TOL(2), where sval(1) is the maximal singular value.
//When TOL(2) < 0, the estimate is indicated by the index of the
//singular value that has the largest logarithmic gap to its successor.
//Default: TOL(2) = -1.
//
//PRINTW is a select for printing the warning messages.
//PRINTW = 1: print warning messages;
// = 0: do not print warning messages.
//Default: PRINTW = 0.
//
//[R,N,SVAL,RCND] = FINDR(S,Y,U,METH,ALG,JOBD,TOL,PRINTW) also returns
//the singular values SVAL, used for estimating the order, as well as,
//if meth = 2, the vector RCND of length 2 containing the reciprocal
//condition numbers of the matrices involved in rank decisions or least
//squares solutions.
//
//[R,N] = FINDR(S,Y) assumes U = [] and default values for the
//remaining input arguments.
//
//See also FINDABCD, FINDAC, FINDBD, FINDBDK, ORDER, SIDENT
//
// V. Sima 18-01-2000.
//
// Revisions:
// V. Sima, July 2000.
//
nin = nargin;
//
// Assumes one batch only.
batch = 4;
conct = 2;
//
if nin<8 then
printw = 0;
end
if nin<7 then
tol(1:2) = [0,-1]
end
if nin<6 then jobd = 2; end
if jobd==[] then jobd = 2,end
if nin<5 then alg = 1;end
if alg==[] then alg = 1;end
if nin<4 then meth = 1;end
if meth==[] then meth = 1;end
if nin<3 then
u = [];
end
//
if meth==1 then
[R,n,sval] = sorder(meth,alg,jobd,batch,conct,s,y,u,tol,printw);
else
[R,n,sval,rcnd] = sorder(meth,alg,jobd,batch,conct,s,y,u,tol,printw);
end
//
// end findR
endfunction
|
