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
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
| function [th,ice,M]=ivgpmf(z,nn,Ts,param,n0,ngpmf,methody,methodu,ic,Tsy,Tsu);
%IVGPMF Computes the continuous-time MISO model parameters from sampled I/O data
% by the Auxiliary Model Instrumental-Variable Generalized Poisson Moment Functionals Approach.
%
% numi(s) -nki.Ts b0i.s^m + ... + bmi
% Gi(s) = ------- = e . ------------------- with a0=1
% den(s) a0.s^n + ... + an
%
% Y(s) = G1(s).U1(s) + ... + Gnu(s).Unu(s)
%
% [th,ice]=ivgpmf(z,nn,Ts,param,n0,ngpmf,'methody','methodu','ic',Tsy,Tsu) returns
% the estimated numerators and denominator under the theta form and the estimated
% initial conditions.
%
% num1 = [b01 ... bm1], ..., numi= [b0i ... bmi]
% den = [1 a1 ... an]
%
% z : the output-input data z=[y u], with y and u as column vectors.
% nn : nn = [n m nk], the orders of the above model.
% n : order of denominator
% m : order of numerator contains as many columns as number of inputs
% nk : delay of the model contains as many columns as number of inputs
% (integer number of sampling period Ts)
% Ts : sampling period (s) at which i/o data are acquired
% param : vector of the Poisson filter parameters [lambda beta]
%
% optional inputs:
% n0 : starting estimation index
% ngpmf : order of the GPMF (minimal-order GPMF then ngpmf=n)
% methody : method used for the simulation of the Poisson filter chain for the output signal
% 'lsim' : simulation using the lsim.m function
% 'tustin' : simulation using the bilinear approximation
% methodu : method used for the simulation of the Poisson filter chain for the input signal
% 'zoh' : simulation assuming a zero order hold on the I/O data
% 'lsim' : simulation using the lsim.m function
% 'tustin' : simulation using the bilinear approximation
% 'ic' : if mentioned, the initial conditions are estimated
% Tsy : Sampling period for the digital simulation of output signal Poisson filter Chain
% Tsy must be smaller than Ts then Spline interpolation is used to minimize numerical
% errors introduced in the simulation of the PFC
% Tsu : Sampling period for the digital simulation of input signal Poisson filter Chain
% Tsu must be smaller than Ts then Spline interpolation is used to minimize numerical
% errors introduced in the simulation of the PFC
% See for further explanations :
% Garnier H.,
% "Identification de modèles paramétriques continus par moments de Poisson",
% Thèse de Doctorat de l'Université Henri Poincaré, Nancy 1, 1995.
% authors : Hugues Garnier - Michel Mensler
% date : 29 Oct. 1996
% revision : 8 Dec. 1998
% name : ivgpmf.m
%
% CRAN - Research Center in Automatic Control of Nancy
% e-mail : <a href="mailto:hugues.garnier@cran.u-nancy.fr">hugues.garnier@cran.u-nancy.fr</a>
%*** Preliminary calculations ***
if nargin<4,error('Number of input arguments incorrect!'),return,end
[Ncap,nz]=size(z); [nr,nc]=size(nn);
ny=1;
nu=(nz-1);
if nz>Ncap,error('Data should be organized in column vectors!'),return,end
if nc<3, nn(3)=0; end
m=nn(2:nu+1);n=nn(1);nk=nn(nu+2:2*nu+1);
nkmax=max(nk);
for counter=1:length(nk),
u(:,counter)=z(nkmax-nk(counter)+1:Ncap-nk(counter),counter+1);
end
y=z(nkmax+1:Ncap,1);
N=size(u,1); % Number of data points
lambda=param(1);
beta=param(2);
if (nargin<11)
Tsu=[];
if (nargin<10)
Tsy=[];
if (nargin<9)
ic=' ';
ice='no estimation of the ic';
if nargin<8
methodu='lsim';
if nargin<7
methody='lsim';
if nargin<6
ngpmf=n;
if nargin<5
n0=1;
end
end
end
end
end
end
end
if ngpmf<n,error('The GPMF order should be greater or equal to the model order n!'),return,end
if Tsy>Ts,
error('The sampling period Tsy should be smaller than the original sampling period Ts of the data!'),
return,end
if Tsu>Ts,
error('The sampling period Tsu should be smaller than the original sampling period Ts of the data!'),
return,end
T=0:Ts:(N-1)*Ts; % time vector
Tma=T; % time vector for the simulation of the auxiliary model
npu=n+1;
% 1. Generation of the I/O data moment functionals
% 1.1 State-Space Representation of the filter chain
I=eye(ngpmf+1);
D=diag([beta*ones(1,ngpmf)],-1)-lambda*I;
Q=[beta zeros(1,ngpmf)]';
% 1.2 Simulation of the I/O data filter chain for the output signal
if (isempty(Tsy)) % Simulation of the PFC at the original sampling time
Tsy=Ts;
Ty=T;
Int_Ratio=1; % Ratio for the output signal interpolation
else % New sampling time for the numerical simulation of the I/O PFC
Int_Ratio=Ts/Tsy;
Ti = (0:Tsy:T(N))';
y = interp1(T,y,Ti,'spline'); % New interpolated output vector
Ty=Ti;
end
Du=zeros(ngpmf+1,1);
if methody(1)=='l' % lsim
Yg=lsim(D,Q,I,Du,y,Ty);
elseif methody(1)=='t' % tustin
[Df,Qf,Cf,Dff]=c2dm(D,Q,I,Du,Tsy,'tustin');
Yg=dlsim(Df,Qf,Cf,Dff,y);
end
Yg=Yg(1:Int_Ratio:length(Yg),ngpmf+1-n:ngpmf+1);
% 2. Generation of the Poisson pulse functions
if ic=='ic' % tests if the initial conditions are wanted
for i=1:ngpmf+1
for j=1:N
P(i,j)=((beta^i)*(((j-1)*Ts)^(i-1))*(exp(-lambda*(j-1)*Ts)))/prod([1:(i-1)]);
end
end
P=P(ngpmf+1-n:ngpmf+1,:);
end
% 3. Definition of the matrices LAM LAMu LAMy and LAMp which are involved in the
% computation of the measures vector
% 3.1 Definition of LAM
for i=1:npu
for j=1:npu
if i>j LAM(i,j)=0;
else LAM(i,j)=((-1)^(j-i))*(prod([1:npu-i])/(prod([1:npu-j])*...
prod([1:j-i])))*(lambda^(j-i))*((beta^(npu-j)));
end
end
end
% 3.3 Definition of LAMy = LAM without its first line
LAMy=LAM(2:npu,:);
% 3.4 Definition of LAMp = LAMy
LAMp=LAMy;
% 4. Generation of the measures vector M (M'.theta = Y)
My = -LAMy*Yg(n0:N,:)';
Mu=[];
Ym = (LAM(1,:)*Yg(n0:N,:)')';
%Simulation of the I/O data filter chain for the input signal
if (isempty(Tsu))
Tsu=Ts; % Simulation of the PFC at the original sampling time
Tu=T;
Int_Ratio=1; % Ratio for the output signal interpolation
uint=u;
else % sampling time for the numerical simulation of the I/O PFC
Int_Ratio=Ts/Tsu;
Ti = (0:Tsu:T(N))';
uint=interp1(T,u,Ti,'spline'); % New interpolated input vector
Tu=Ti;
end
for kk=1:nu,
if methodu(1)=='z' % zoh
[Df,Qf]=c2d(D,Q,Tsu);
Ug=dlsim(Df,Qf,I,Du,uint(:,kk));
elseif methodu(1)=='l' % lsim
Ug=lsim(D,Q,I,Du,uint(:,kk),Tu);
elseif methodu(1)=='t' % tustin
[Df,Qf,Cf,Dff]=c2dm(D,Q,I,Du,Tsu,'tustin');
Ug=dlsim(Df,Qf,Cf,Dff,uint);
end
Ug=Ug(1:Int_Ratio:length(Ug),ngpmf+1-n:ngpmf+1);
LAMu=LAM(npu-m(kk):npu,1:npu);
Mu = [Mu ; LAMu*Ug(n0:N,:)'];
end
% 4.1 with the initial conditions terms
if ic=='ic'
Mic = -LAMp*P(:,n0:N);
M = [My' Mu' Mic'];
else
% 4.2 without the initial conditions term
M = [My' Mu'];
end
[N_LS,n_par]=size(M); % N_LS : number of data used in the LS algorithm
% n_par : number of parameters to be estimated
% 5. Estimation by off-line LS
thmc=pinv(M)*Ym;
% 6. Simulation of the estimated system (yma:aux. model output)
[N_LS,n_par]=size(M); % N_LS : number of data used in the LS algorithm
% n_par : number of parameters to be estimated
II=[1 -Ts nu 0 m+1 0 0 n*ones(1,nu) nk];
nb=n+sum(m+1);
th=zeros(nb+3,max([length(II) nb 7]));
th(1,1:length(II))=II;
ti=fix(clock); ti(1)=ti(1)/100;
th(2,2:6)=ti(1:5);
th(2,7)=1;
thmc_i=[];
for counter=1:nu,
thmc_i=[thmc_i thmc(1:n)'];
end
th(3,1:nu*n+sum(m+1))=[thmc(n+1:n+sum(m+1))' thmc_i];
[Amc,Bmc,Cmc,Dmc]=th2ss(th);
yma=lsim(Amc,Bmc,Cmc,Dmc,u,Tma);
% 7. Generation of the moment functionals of the aux. model output
if (isempty(Tsy)) % Simulation of the PFC at the original sampling time
Tsy=Ts;
Ty=T;
Int_Ratio=1; % Ratio for the output signal interpolation
else % New sampling time for the numerical simulation of the I/O PFC
Int_Ratio=Ts/Tsy;
Ti = (0:Tsy:T(N))';
yma = interp1(T,yma,Ti,'spline'); % New interpolated output vector
Ty=Ti;
end
if methody(1)=='z'
Ygma=dlsim(Df,Qf,I,Du,yma);
elseif methody(1)=='l'
Ygma=lsim(D,Q,I,Du,yma,Ty);
elseif methody(1)=='t'
Ygma=dlsim(Df,Qf,Cf,Dff,yma);
end
Ygma=Ygma(1:Int_Ratio:length(Ygma),ngpmf+1-n:ngpmf+1);
% 8. Generation of the aux. model measures vector
Myma=-LAMy*Ygma(n0:N,:)';
if ic=='ic'
% 8.1 with the initial conditions terms
Mma=[Myma' Mu' Mic'];
else
% 8.2 without the initial conditions terms
Mma=[Myma' Mu'];
end
% 9. Estimation by aux. model IV
thiv=(Mma'*M)\Mma'*Ym;
I=[1 -Ts nu n m+1 0 0 zeros(1,nu) nk];
nb=n+sum(m+1);
th=zeros(nb+3,max([length(I) nb 7]));
th(1,1:length(I))=I;
ti=fix(clock); ti(1)=ti(1)/100;
th(2,2:6)=ti(1:5);
th(2,7)=2;
th(3,1:n+sum(m+1))=thiv(1:n+sum(m+1))';
if ic=='ic'
ice=thiv(npu+sum(m+1):length(thiv))';
end |
