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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
|
function [x0,V,rcnd]=inistate(A,B,C,D,y,u,tol,printw)
x0=[];V=[];rcnd=[];
[nargout,nargin] = argn(0)
//INISTATE Estimates the initial state of a discrete-time system, given the
// (estimated) system matrices, and a set of input/output data.
//
// X0 = INISTATE(SYS,Y,U,TOL,PRINTW) estimates the initial state X0 of
// the discrete-time system SYS = (A,B,C,D), using the output data Y
// and the input data U. The model structure is :
//
// x(k+1) = Ax(k) + Bu(k), k >= 1,
// y(k) = Cx(k) + Du(k),
//
// The vectors y(k) and u(k) are transposes of the k-th rows of Y and U,
// respectively.
// Instead of the first input parameter SYS (an ss object), equivalent
// information may be specified using matrix parameters, for instance,
// X0 = INISTATE(A,B,C,Y,U); or X0 = INISTATE(A,C,Y);
//
// TOL is the tolerance used for estimating the rank of matrices.
// If TOL > 0, then the given value of TOL is used as a lower bound
// for the reciprocal condition number.
// Default: prod(size(matrix))*epsilon_machine where epsilon_machine
// is the relative machine precision.
//
// PRINTW is a select for printing the warning messages.
// PRINTW = 1: print warning messages;
// = 0: do not print warning messages.
// Default: PRINTW = 0.
//
// [x0,V,rcnd] = INISTATE(SYS,Y,U,TOL,PRINTW) returns, besides x0,
// the orthogonal matrix V which reduces the system state matrix A to
// a real Schur form, as well as an estimate of the reciprocal condition
// number of the coefficient matrix of the least squares problem solved.
//
// See also FINDBD, FINDX0BD
//
// V. Sima 13-05-2000.
//
// For efficiency, most errors are checked in the mexfile findBD. Also,
// except for scalars, the input parameters are not copied, but renamed.
//
// Revisions:
// V. Sima, July 2000.
//
ni = nargin;
if mtlb_isa(A,"lti") then
// Get the system matrices of the ss object, and the remaining parameters.
// General call x0 = inistate(A,B,C,D,y,u,tol,printw);
// Special call x0 = inistate(sys,y,u,tol,printw);
//
if A.dt=="c" then
error(msprintf(gettext("%s: Wrong values for input argument #%d: Discrete time system expected.\n"),"inistate",1))
end
if ni<2 then
error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"inistate",2));
end
[As,Bs,Cs,Ds] = abcd(A)
[ny,nu] = size(A);
[ty,p] = size(B);
if ni>2 then
[tu,m] = size(C);
if ~(((tu==ty|tu==0)&(m==nu))&(ty>1)) then
tol = C;
// Special call x0 = inistate(sys,y,tol,printw);
if ni>3 then
printw = D;
else
printw = 0;
end
// Below, B means y !
[x0,Vl,rcndl] = findBD(1,3,As,Cs,B,tol,printw);
else
if ni>3 then
// Special call x0 = inistate(sys,y,u,tol,printw);
tol = D;
if ni>4 then
printw = y;
else
printw = 0;
end
else
tol = 0;
printw = 0;
end
// Below, B means y, and C means u !
if norm(Ds,1)==0 then
[x0,Vl,rcndl] = findBD(1,2,1,As,Bs,Cs,B,C,tol,printw);
else
[x0,Vl,rcndl] = findBD(1,2,2,As,Bs,Cs,Ds,B,C,tol,printw);
end
end
else
// Special call x0 = inistate(sys,y);
// Below, B means y !
[x0,Vl,rcndl] = findBD(1,3,As,Cs,B);
end
//
else
// The system matrices are directly specified.
// General call x0 = inistate(A,B,C,D,y,u,tol,printw);
// Special calls x0 = inistate(A,B,C,y,u,tol,printw);
// x0 = inistate(A,C,y,tol,printw);
//
if ni<3 then
error(msprintf(gettext("%s: Wrong number of input arguments: At least %d expected.\n"),"inistate",3));
end
[m2,n2] = size(B);
[m3,n3] = size(C);
if ni>=4 then
[m4,n4] = size(D);
if ni>=5 then
[m5,n5] = size(y);
if ni>=6 then
[m6,n6] = size(u);
if ni>=7 then
if (ni==7)&(m6*n6>1) then
// Special call x0 = inistate(A,B,C,D,y,u,tol);
[x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y,u,tol);
elseif ni==7 then
// Special call x0 = inistate(A,B,C,y,u,tol,printw);
// Below, D means y and y means u !
printw = tol;
tol = u;
[x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y,tol,printw);
else
[x0,Vl,rcndl] = findBD(1,2,2,A,B,C,D,y,u,tol,printw);
end
else
if m6*n6>1 then
[x0,Vl,rcndl] = findBD(1,2,2,A,B,C,D,y,u);
else
// Special call x0 = inistate(A,B,C,y,u,tol);
// Below, y means U and D means y !
tol = u;
[x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y,tol);
end
end
else
// Special calls x0 = inistate(A,B,C,y,u);
// x0 = inistate(A,C,y,tol,printw);
if m5*n5>1 then
// Below, y means u and D means y !
[x0,Vl,rcndl] = findBD(1,2,1,A,B,C,D,y);
else
// Below, C means y and B means C !
tol = D;
printw = y;
[x0,Vl,rcndl] = findBD(1,3,A,B,C,tol,printw);
end
end
else
// Below, D means tol, C means y, and B means C !
[x0,Vl,rcndl] = findBD(1,3,A,B,C,D);
end
else
// Below, C means y, and B means C !
[x0,Vl,rcndl] = findBD(1,3,A,B,C);
end
end
//
if nargout>2 then
V = Vl;
rcnd = rcndl;
elseif nargout>1 then
V = Vl;
end
//
// end inistate
endfunction
|
