Global Index (short | long) | Local contents | Local Index (short | long)
This case study concerns data collected from a laboratory scale "hairdryer". (Feedback's Process Trainer PT326; See also page 440 in Ljung,1987). The process works as follows: Air is fanned through a tube and heated at the inlet. The air temperature is measured by a thermocouple at the outlet. The input (u) is the voltage over the heating device, which is just a mesh of resistor wires. The output is the outlet air temperature (or rather the voltage from the thermocouple).
| This script calls | This script is called by |
|---|---|
echo on,
echo off
% Lennart Ljung
% Copyright (c) 1986-98 by The MathWorks, Inc.
% $Revision: 3.4 $ $Date: 1997/12/02 03:43:18 $
echo on
% First load the data:
load dryer2
% Vector y2, the output, now contains 1000 measurements of
% temperature in the outlet airstream. Vector u2 contains 1000
% input data points, consisting of the voltage applied to the
% heater. The input was generated as a binary random sequence
% that switches from one level to the other with probability 0.2.
% The sampling interval is 0.08 seconds.
pause % Press a key to continue
% Select the 300 first values for building a model:
z2 = [y2(1:300) u2(1:300)];
% Plot the interval from sample 200 to 300:
% (0.08 is the sampling interval -- for correct scales--)
pause, idplot(z2,200:300,0.08), pause % Press any key for plot.
% Remove the constant levels and make the data zero-mean:
z2 = dtrend(z2);
% Let us first estimate the impulse response of the system
% by correlation analysis to get some idea of time constants
% and the like:
pause, ir = cra(z2); pause % Press a key for plot
% We can form the corresponding step response by integrating
% the impulse response:
stepr = cumsum(ir);
pause, % press a key for plot
newplot; plot(stepr), title('Step response'), xlabel('Time (samples)'),pause
% Compute a model with two poles, one zero, and three delays:
th = arx(z2,[2 2 3]);
th = sett(th,0.08); % Set the correct sampling interval.
% Let's take a look at the model:
pause, present(th) % Press any key to see model.
pause % Press any key to continue.
% How good is this model? One way to find out is to simulate it
% and compare the model output with measured output. We then
% select a portion of the original data that was not used to build
% the model, viz from sample 800 to 900:
u = dtrend(u2(800:900)); y = dtrend(y2(800:900));
yh = idsim(u,th);
pause % Press any key for plot.
plot([yh y]),title('Simulated (blue/solid) and measured (green/dashed) output')
xlabel('Time (samples)'),pause
% The same result is obtained with the more efficient command compare:
compare([y u],th);pause
% The agreement is quite good. Compute the poles and zeros of
% the model:
zpth = th2zp(th);
pause, zpplot(zpth), pause % Press any key for plot.
% Now compute the transfer function of the model:
gth = th2ff(th);
pause, bodeplot(gth), pause % Press any key for Bode plot.
% We can compare this transfer function with what is obtained
% by a non-parametric spectral analysis method:
gs = spa(z2);
gs = sett(gs,0.08); %setting the sampling interval
% Press a key for a comparison with the transfer function from the
% parametric model.
pause, bodeplot([gs gth]), pause
% We can compute the step response of the model as follows.
step = ones(20,1);
mstepr = idsim(step,th);
% This step response can be compared with that from correlation
% analysis:
pause, % Press any key for plot.
plot([stepr, mstepr])
title('Solid/blue: Step response from CRA; Dashed/green: Step response from ARX-model'),xlabel('Time (samples)'),pause
% Finally we can plot the uncertainty of the step response of the model.
% The model comes with an estimate of its own uncertainty. In this case
% 10 different step responses are plotted, corresponding to models drawn
% from the distribution of the true system (according to our model).
% Press a key for the plot.
pause, idsimsd(step,th), pause % Press any key for plot.
echo off
set(gcf,'NextPlot','replace');