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
|
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA
// Copyright (C) ENPC - J-Ph. Chancelier
//
// 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 [arc,resid]=armax1(r,s,q,y,u,b0f)
// [arc, resid] = armax1(r, s, q, y, u, [b0f])
//
// Compute the coefficient of monodimensional ARMAX
// A(z^-1)y= B(z^-1)u + D(z^-1)sig*e(t)
// e(t) is a white noise of variance 1
// A(z)= 1+a1*z+...+a_r*z^r; ( r=0 => A(z)=1)
// B(z)= b0+b1*z+...+b_s z^s ( s=-1 => B(z)=0)
// D(z)= 1+d1*z+...+d_q*z^q ( q=0 => D(z)=1)
//
// Intput:
// y: output process
// u: input process
// r, s and q: auto-regression orders r >= 0, s >= 1 and moving average q.
// b0f is a optional parameter. By default, b0f is 0 and it means that
// this parameter must be identified. If b0f is 1, then b0f is supposed
// to be zero and is not identified.
//
// Output:
// arc is the tlist with fields
// a is the vector [1,a1,...,a_r]
// b is the vector [b0, ..., b_s]
// d is the vector [1, d1, ..., d_q]
// sig is the estimated standard deviation
// resid = [sig*echap(1), ..., ]
//
// Method:
// Cfre : Eykhoff (trends and progress in system identification) page 91
// Introducing
// z(t)=[y(t-1),..,y(t-r),u(t),...,u(t-s),e(t-1),...,e(t-q)] and
// coef= [-a1,..,-ar,b0,...,b_s,d1,...,d_q]', we obtain
// y(t)= coef'* z(t) + sig*e(t).
// We use the sequential version of the AR estimation where e(t-i) is
// replaced by its estimated (Method RLLS).
// If q=0, it is a sequential version of the least squares algorithm given
// in armax function
[lhs,rhs]=argn(0)
if rhs<=5,b0f=0;end
if s==-1,b0f=0;end // Seems not natural, but makes things work
u=matrix(u,1,-1);y=matrix(y,1,-1); //make u and y row vectors
[n1,n2]=size(y)
if size(y,"*")<>size(u,"*") then
error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same numbers of elements expected.\n"),..
"armax1",4,5));
end
//
t0=max(max(r,s+1),1)+1;
if r<>0;XTM1=y((t0-1):-1:(t0-r));else XTM1=[];end
if s<>-1;UTM1=u(t0-b0f:-1:(t0-s));else UTM1=[];end
if q<>0;ETM1=0*ones(1,q);else ETM1=[];end
npar=r+s+1-b0f+q
CTM1=0*ones(npar,1);
ZTM1=[XTM1,UTM1,ETM1]';
PTM1=10.0*eye(npar,npar);
for t=t0+1:n2,
if r<>0;XT=[ y(t-1), XTM1(1:(r-1))];else XT=[];end
if s<>-1;UT=[ u(t-b0f), UTM1(1:(s-b0f))];else UT=[];end
eeTM1=y(t-1)- CTM1'*ZTM1;
if q<>0;ET=[ eeTM1, ETM1(1:(q-1))];else ET=[];end
ZT=[XT,UT,ET]';
//
KT=PTM1*ZT*(1/(1+ ZT'*PTM1*ZT))
CT=CTM1+KT*(y(t)-ZT'*CTM1)
PT=PTM1-KT*ZT'*PTM1
XTM1=XT;UTM1=UT;CTM1=CT;ETM1=ET;ZTM1=ZT;PTM1=PT;
end
// The coefficient a, b and d are extracted
//
if r<>0;a=[1;-CT(1:r)]';else a=1;end
if s<>-1;
if b0f==1,b=[0;CT(r+1:(r+s+1-b0f))]';else
b=[CT(r+1:(r+s+1-b0f))]';end
if q<>0;d=[1;CT(r+s+2-b0f:(r+s+q+1-b0f))]';else d=[1];end
else
b=0;
if q<>0;d=[1;CT(r+s+1-b0f:(r+s+q-b0f))]';else d=[1];end
end
// Simulation to get the prediction error
//
[sig,resid]=epred(r,s,q,CTM1,y,u,b0f);
arc=armac(a,b,d,1,1,sig);
endfunction
function [sig,resid]=epred(r,s,q,coef,y,u,b0f)
//=============================================
// [sig,resid] = epred(r,s,q,coef,y,u,b0f)
// Used by armax1 function to compute the prediction error
// coef= [-a1,..,-ar,b0,...,b_s,d1,...,d_q]'
// or
// coef= [-a1,..,-ar,b1,...,b_s,d1,...,d_q]' si b0f=1
//!
[n1,n2]=size(y);
t0=max(max(r,s+1),1)+1;
if r<>0;XTM1=y((t0-1):-1:(t0-r));else XTM1=[];end
if s<>-1;UTM1=u(t0-b0f:-1:(t0-s));else UTM1=[];end
if q<>0;ETM1=0*ones(1,q);else ETM1=[];end
npar=r+s+1-b0f+q
ZTM1=[XTM1,UTM1,ETM1]';
resid=0*ones(1,n2);
for t=t0+1:n2,
if r<>0;XT=[ y(t-1), XTM1(1:(r-1))];else XT=[];end
if s<>-1;UT=[ u(t-b0f), UTM1(1:(s-b0f))];else UT=[];end
resid(t)=y(t-1)- coef'*ZTM1;
if q<>0;ET=[ resid(t), ETM1(1:(q-1))];else ET=[];end
ZT=[XT,UT,ET]';
XTM1=XT;UTM1=UT;ETM1=ET;ZTM1=ZT;
end
sig=sqrt(sum(resid.*resid)/size(resid,"*"))
endfunction
|
