Global Index (short | long) | Local contents | Local Index (short | long)
This case study is based on the same "hairdyer" data as study # 1. The process consists of air being fanned through a tube. The air is heated at the inlet of the tube, and the input is the voltage applied to the heater. The output is the temperature at the outlet of the tube.
| This script calls | This script is called by |
|---|---|
echo on
echo off
% Lennart Ljung
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 2.3 $ $Date: 1997/12/02 03:43:15 $
echo on
% First, load the data:
load dryer2
% Form a data set for estimation of the first half, and a reference
% set for validation purposes of the second half:
ze = [y2(1:500) u2(1:500)];
zr = [y2(501:1000) u2(501:1000)];
% Detrend each of the sets:
ze = dtrend(ze); zr = dtrend(zr);
% Take a look at the data:
pause, idplot(ze,200:350), pause % Press any key to continue
% We shall work with models of the following kind:
% y(t) + a1*y(t-1) + ... + ana*y(t-na) = b1*u(t-nk) + ... + bnb*u(t-nb-nk+1)
% First we try and determine the time delay (nk). We then
% select a second order model (na=nb=2), and try out every time delay
% between 1 and 10. The loss function for the different models
% are computed using the reference data set:
V = arxstruc(ze,zr,struc(2,2,1:10));
% We now select that delay that gives the best fit for the
% reference set:
[nn,Vm] = selstruc(V,0);
% The chosen structure was
nn
pause % Press any key to continue.
% We can also check how the fit depends on the delay. The logarithms
% of a quadratic loss function are given as the first row, while
% the indexes na, nb and nk are given as a column below the correspon-
% ding loss function.
pause % Press any key to continue.
Vm
% The choice of three delays is thus rather clear. Now test the
% orders. We check the fit for all 25 combinations of up to 5
% b-parameters and up to 5 a-parameters, all with delay 3:
V = arxstruc(ze,zr,struc(1:5,1:5,3));
pause % Press any key to continue.
% The best fit for the reference data set is obtained for
nn = selstruc(V,0)
pause % Press any key to continue.
% It may be advisable to check yourself how much the fit is
% improved for the higher order models. The following plot
% shows the fit as a function of the number of parameters used.
% (Note that several different model structures use the same
% number of parameters.) You will then be asked to enter the
% number of parameters you think gives a sufficiently good fit.
% The routine then selects that structure with this many parameters
% that gives the best fit.
pause, nns = selstruc(V) % Press a key to continue.
% The best fit is thus obtained for nn = [4 4 3], while we
% see that the improved fit compared to nn = [2 2 3] is rather
% marginal.
pause % Press any key to continue.
% Let's compute the fourth order model:
th4 = arx(ze,[4 4 3]);
% We now check the pole-zero configuration for the fourth order
% model. We also include confidence regions for the poles and zeros
% corresponding to 3 standard deviations.
pause, zpplot(th2zp(th4),3), pause % Press any key to for plot.
% It is thus quite clear that the two complex conjugated pole-zero
% pairs are likely to cancel, so that a second order model would be
% adequate. We thus compute
th2 = arx(ze,[2 2 3]);
% We can test how well this model is capable of reproducing the
% reference data set. To compare the simulated output from the
% model with the actual output (plotting the mid 200 data points)
% we invoke (press key for plot)
pause,compare(zr,th2,inf,150:350);pause
% We can also check the residuals ("leftovers") of this model,
% i.e., what is left unexplained by the model.
pause % Press any key to continue.
e = resid(ze,th2); pause
plot(e), title('The residuals'),pause
% We see that the residual are quite small compared to the signal
% level of the output, that they are reasonably (although not
% perfectly) well uncorrelated with the input and between themselves.
% We can thus be (provisionally) satisfied with the model th2.
pause % Press any key to return to main menu.
echo off
set(gcf,'NextPlot','replace');