Global Index (short | long) | Local contents | Local Index (short | long)
In this example we compare several different methods of identification.
| 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:09 $
echo on
% We start by forming some simulated data. Here are the
% numerator and denominator coefficients of the transfer
% function model of a system:
B = [0 1 0.5];
A = [1 -1.5 0.7];
% These can be put into the special "theta" structure used by
% the toolbox:
th0 = poly2th(A,B,[1 -1 0.2]);
% The third argument, [1 -1 0.2], gives a characterization
% of the disturbances that act on the system.
pause % Press any key to continue.
% Now we generate an input signal U, a disturbance signal E,
% and simulate the response of the model to these inputs:
u = sign(randn(350,1));
e = randn(350,1);
y = idsim([u e],th0);
z = [y u]; % A matrix with the output and inputs as column vectors.
% Plot first 100 values of the input U and the output Y.
pause, idplot(z,1:100), pause % Press any key for plot.
% Now that we have corrupted, simulated data, we can estimate models
% and make comparisons. Let's start with spectral analysis. We
% invoke SPA which returns a spectral analysis estimate of
% the transfer function GS and the noise spectrum NSS:
[GS,NSS] = spa(z);
pause, bodeplot(GS), pause % Press any key for plot.
% Now use the Least Squares method to find an ARX-model with
% two poles, one zero, and a single delay on the input:
a2 = arx(z,[2 2 1]);
pause, present(a2), pause % Press any key to present results.
% Now use the instrumental variable method (IV4) to find a model
% with two poles, one zero, and a single delay on the input:
i2 = iv4(z,[2 2 1]);
pause, present(i2), pause % Press any key to present results.
% Calculate transfer functions for the two models obtained by
% ARX and IV4. Show Bode plots comparing the Spectral Analysis
% estimate with the ARX and IV4 estimates:
[Ga2,NSa2] = th2ff(a2); % From ARX
Gi2 = th2ff(i2); % From IV4
pause, bodeplot([GS Ga2 Gi2]), pause % Press any key for plots.
% Calculate and plot residuals for the model obtained by IV4:
e = resid(z,i2); pause
plot(e), title('Residuals, IV4 method'), pause
% Now compute second order ARMAX model:
am2 = armax(z,[2 2 2 1]);
pause % Press any key to continue.
% Now compute second order BOX-JENKINS model:
bj2 = bj(z,[2 2 2 2 1]);
% Compute transfer functions and noise spectra for these models:
[Gam2,NSam2] = th2ff(am2);
[Gbj2,NSbj2] = th2ff(bj2);
% Compare these transfer functions with those obtained by the
% IV4 method and the spectral analysis method:
pause % Press any key for plots.
bodeplot([Gam2 Gbj2 Gi2 GS]), pause
bodeplot([NSam2 NSbj2 NSS]), pause,
% Plot residuals for the ARMAX and BOX-JENKINS models
e1 = resid(z,am2); pause
e2 = resid(z,bj2); pause
plot([e1 e2]), title('ARMAX and BOX-JENKINS residuals'), pause
% Simulate the ARMAX and BOX-JENKINS models with the real input
yham2 = idsim(u,am2);
yhbj2 = idsim(u,bj2);
pause % Press any key for plot.
plot([y yham2 yhbj2]), ...
title('Real output and simulated model outputs '), pause
% Now we'll compute the poles and zeros of the ARMAX and
% BOX-JENKINS models
zpam2 = th2zp(am2);
zpbj2 = th2zp(bj2);
% Press any key for plot.
pause % Press keys to advance poles and zeros on plot, from the ARMAX
% model to the Box-Jenkins model.
zpplot([zpam2 zpbj2]), pause,
% Clearly, AM2 and BJ2 both give good residuals and simulated outputs.
% They are also in good mutual agreement. Let's check which has the
% smaller AIC:
[am2(2,1) bj2(2,1)]
pause % Press any key to continue.
% Since we generated the data, we enjoy the luxury of comparing
% the model to the true system:
[G0,NS0] = th2ff(th0);
zp0 = th2zp(th0);
pause % Press any key for plots of comparisons.
bodeplot([G0 Gam2]), pause
bodeplot([NS0 NSam2]), pause
zpplot([zp0 zpam2]), pause
echo off
set(gcf,'NextPlot','replace');