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
|
// 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 [archap,la,lb,sig,resid]=armax(r,s,y,u,b0f,prf)
// [la, lb, sig, resid] = armax(r, s, y, u, [b0f, prf])
// armax identification
// Compute the coefficients of a n-dimensional ARX
// A(z^-1)y= B(z^-1)u + sig*e(t)
// e(t) is a white noise of n-dimensional Id variance
// sig is a n by n matrix
// 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)
// Method:
// Cfre : Eykhoff (trends and progress in system identification) page 96
// Introducing z(t)=[y(t-1),..,y(t-r),u(t),...,u(t-s)] and
// oeff = [-a1,..,-ar,b0,...,b_s], we obtain
// y(t)= coef* z(t) + sig*e(t)
// and the algorithm is finding coef which minimizes
// sum_{t=1}^N ( [y(t)- coef'z(t)]^2)
//
// Input:
// y: output process y(ny, n); ny: dimension,
// n: sample size
// u: input process u(nu, n); nu: dimension
// n: sample size
// r and s: auto-regression orders r >=0 and s >=-1
// b0f is 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.
// prf is optional parameter for display control.
// if prf = 1, a display is done (the default value)
// Output:
// la is the list(a,a+eta,a-eta); eta: estimated standard deviation
// a=[Id,a1,a2,...,ar] where ai: ny-by-ny matrix
// lb is the list(b,b+etb,b-etb); etb: estimated standard deviation
// b=[b0,.....,b_s] where bi is nu-by-nu matrix
// sig is the estimated standard deviation of noise
// resid=[ sig*e(t0),....];
// t0=max(max(r,s)+1,1));
//
// Example:
// Enter the command [a,b,sig,resid]=armax();
// to see an example in dimension 1.
// Auteur: J-Ph. Chancelier ENPC Cergrene
//!
// Copyright INRIA
[lhs,rhs]=argn(0)
if rhs==0,write(%io(2),"/ / y_t = 0.2*u_{t-1}+0.01*e(t)");
write(%io(2)," rand(''normal''); u=rand(1,1000);");
write(%io(2)," y=arsimul([1],[0,0.2],1,0.01,u);");
write(%io(2)," [archap,a,b,sig,resid]=armax(0,1,y,u)");
write(%io(2),"/ / we must find a=1,b=[0,0.2]''");
rand("normal"),u=rand(1,1000);
y=arsimul([1],[0,0.2],1,0.01,u);
[archap,la,lb,sig,resid]=armax(0,1,y,u,1);
return
end
if rhs<=5,prf=1;end
if rhs<=4,b0f=0;end
[ny,n2]=size(y)
[nu,n2u]=size(u)
// Compute zz matrix as
// zz(:,j)=[ y(t-1),...,y(t-r),u(t),...,u(t-s)]', with t=t0-1+j
// zz can be computed from t = t0
t0=max(max(r,s)+1,1);
if r==0;if s==-1;error(msprintf(gettext("%s: Wrong value for input arguments: If %s and %s nothing to identify.\n"),"armax","r==0","s==-1"))
end;end
z=[];
if r>=1;for i=1:r,z=[z ; y(:,t0-i:(n2-(i)))];end;end
if s>=-1;for i=b0f:s,z=[z ; u(:,t0-i:(n2-(i)))];end;end
zz= z*z';
zy= z*y(:,t0:n2)';
// Rank test
[nzl,nzc]=size(zz);
k=rank(zz);
if k<>nzl then
warning(msprintf(gettext("%s: %s is numerically singular.\n"),"armax","z*z''"));
end;
pv=pinv(zz);
coef=(pv*zy)';
// The residual noise
resid=y(:,t0:n2) - coef*z;
// The variance of the residual noise
sig2= resid*resid'/(n2-t0+1)
// The standard deviation
sig=sqrtm(sig2);
a=[eye(ny,ny),-coef(:,1:r*ny)];
if b0f==0 then
b=coef(:,r*ny+1:(s+1)*nu+r*ny);
else
b=[0*ones(ny,nu),coef(:,r*ny+1:r*ny+s*nu)];
end
// For the SISO systems, the estimated standard deviation is added.
// It is to be done for the MIMO
if ny == 1,
dve=sqrt(diag(sig*pv,0))';
la=list(a,a+[0,dve(1:r)],a-[0,dve(1:r)]);
if b0f==0 then
lb=list(b,b+dve(r+1:r+s+1),b-dve(r+1:r+s+1)),
else
lb=list(b,b+[0,dve(r+1:r+s)],b-[0,dve(r+1:r+s)]);
end
else
la=a;lb=b;
end
// If prf = 1, the display is done
//si prf vaut 1 on donne un display
archap=armac(a,b,eye(ny,ny),ny,nu,sig);
if prf==1;
if ny==1,
[nla,nca]=size(la(2));
mprintf(gettext("%s: Standard deviation of the estimator %s:\n"),"armax","a");
form_a="(5x,"+string(nca)+"(f7.5,1x))";
write(%io(2),la(2)-a,form_a);
if nu<>0 then
mprintf(gettext("%s: Standard deviation of the estimator %s:\n"),"armax","b");
[nlb,ncb]=size(lb(2));
write(%io(2),lb(2)-b,"(5x,"+string(ncb)+"(f7.5,1x))");
end
end
end
endfunction
|
