ch4_1k.m

《MATLAB7.0控制系统应用实例》 作者：刘叔军

%  Select a demo number: 12   
% 1) A glass tube manufacturing process
   
%CASE STUDY NUMBER 1 : A GLASS TUBE DRAWING FACTORY
%
% In this Case Study we shall study data from a glass tube factory. The
% experiments and the data are discussed in
%
% V. Wertz, G. Bastin and M. Heet: Identification of a glass tube
% drawing bench. Proc. of the 10th IFAC Congress, Vol 10, pp 334-339
% Paper number 14.5-5-2. Munich August 1987.
%
% The the output of the process is the thickness and the diameter of
% the manufactured tube. The inputs are the air-pressure inside the
% tube and the drawing speed.
%
% We shall in this tutorial study the process from the input speed
% to the output thickness.
%
load thispe25.mat
%
% The data ar contained in the variable glass:
 glass  % Strike a key to continue

                                   

% Split it into two halves, one for
% estimation and one for validation:
echo off

ze = glass(1001:1500);
zv = glass(1501:2000,:);
figure, plot(ze)
% A close-up: 
figure, plot(ze(101:200))

% Let us remove the mean values:
ze = detrend(ze);
zv = detrend(zv);

% The sampling interval of the data is one second. We may
% detect some rather high frequencies in the output. Let us therefore
% first compute the input and output spectra:
sy = spa(ze(:,1,[]));
su = spa(ze(:,[],1));
figure,  bode(sy,su)% Strike key to advance plots

% Note that the input has very little relative energy above 1 rad/sec
% while a substantial part of the output's energy comes from fequencies above
% 1 rad/sec. There are thus some high frequency disturbances that may cause
% some problem for the model building.
% We first compute the impulse response, using part of the data:
figure, impulse(ze,'sd',3,'fill')

% We see a quite substantial delay of about 12 samples in the impulse response.
% We also as a preliminary test compute the spectral analysis estimate:
g = spa(ze);
figure, bode(g,'sd',1,'fill')

% We note, among other things, that the high frequency behavior is quite
% uncertain.
% Let us do a quick check if we can pick up good dynamics by just computing
% a fourth order ARX model using the estimation data and simulate that
% model using the validation data. We know that the delay is about 12:

m1 = arx(ze,[4 4 12]);
figure, compare(zv,m1);

% A close up:
figure, compare(zv,m1,inf,101:200)

%
% There are clear difficulties to deal with the high frequency components
% of the output. That in conjunction with the long delay suggests that we
% decimate the data  by four (i.e. low-pass filter it and pick every fourth
% value):
zd = idresamp(detrend(glass),4);
zde = zd(1:500);
zdv = zd(501:size(zd,'N'));


% Let us find a good structure for the decimated data.
% First compute the impulse response:
figure, impulse(zde,'sd',3,'fill')

% We see that the delay is about 3 samples(which is consistent with what we saw
% above.
% What does the  default model give?
figure, compare(zdv,pem(zde))

% Seems much better that for the undecimated data. Let's go on to find
% a good model structure. First we look for the delay:
V = arxstruc(zde,zdv,struc(2,2,1:30));
nn = selstruc(V,0)


% Then we test several different orders with and around this delay:
V = arxstruc(zde,zdv,struc(1:5,1:5,nn(3)-1:nn(3)+1));

% Make your own choice of orders based on the figure to follow.
% (The "model quality" comments to be made below refer to the default choice)
selstruc(V)

% Let's compute and check the model for whatever model you just preferred:
m2 = arx(zde,nn);
 
figure, compare(zdv,m2,inf,21:150);

% This seems reasonably good!
% We can also compare the bodeplots of the model m1 and m2:
figure, bode(m1,m2)

% We see that the two models agree well up to the Nyquist frequency of
% the slower sampled data.
% Let's test the residuals:
 figure,resid(zdv,m2);

% This seems OK.
% What does the pole-zero diagram tell us?
 figure, zpplot(m2,'sd',3)

% Some quite clear indication of pole-zero cancellations!
% This shows that we should be able to do well with lower order
% models. We shall return to that later.
% Let us however also check if we can do better by modeling the noise
% separately in a Box-Jenkins model:
mb=bj(zde,[nn(2) 2 2 nn(1) nn(3)])

% Is this better than the ARX-model in simulation ?
% (press key)
 figure,compare(zdv,m2,mb)
 
% The Box-Jenkins model is thus capable of describing the validation
% data set somewhat better. How about the Bode plots?
 figure, bode(m2,mb)

% The models agree quite well, apart from the funny twist in the BJ model.
% This twist, however, seems to have some physical relevance, since we
% get better simulations in the latter case.
% Finally we can compare the FPE:s of the two models:
[fpe(m2) fpe(mb)]

% To summarize: After decimation it is quite possible to build simple models
% that are capable of reproducing the validation set in a good manner.
% We can however do quite well with a much simpler model; a [1 1 3] ARX model:
m3=arx(zde,[1 1 3]);

% Simulation:
 figure, compare(zdv,m3,mb)

% 5-step ahead prediction:
 figure,compare(zdv,m3,mb,5)