ch4_1b.m

清华大学《matlab 控制系统应用与实例》第一部分M文件的源码

%   In this example we compare several different methods of
%   identification.
echo off

%       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 model object used by
%       the toolbox:
m0 = idpoly(A,B,[1 -1 0.2]);

%       The third argument, [1 -1 0.2], gives a characterization
%       of the disturbances that act on the system.

%       Now we generate an input signal U, a disturbance signal E,
%       and simulate the response of the model to these inputs:
u = idinput(350,'rbs'); % A random binary signal of length 350
u = iddata([],u);
e = iddata([],randn(350,1)); % White Gaussian noise with unit variance
y = sim(m0,[u e]) % Both u and y are now IDDATA objects

 %       Collect the input/output data as an IDDATA data object: 
z = [y,u];

 z  % Press any key to get information about the data.

%       Plot first 100 values of the input U and the output Y.
figure, plot(z(1:100))

%       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 frequency function and the noise spectrum:
GS = spa(z);

figure, bode(GS)

figure, bode(GS,'sd',3,'fill') % Showing the uncertainty (3 standard deviations)

%       We start with a default model (state-space model computed
%       by a prediction error method:
m = pem(z);

m  % Press any key to present results.

%       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]);

 a2  % 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]);

i2    % Press any key to present results.

%       Show Bode plots comparing the Spectral Analysis
%       estimate with the ARX and IV4 estimates:
figure, bode(GS,m,a2,i2) % Press any key for plots.

%       To set the line markers use
figure, bode(GS,'b',m,'m',a2,'g',i2,'k')  % Press any key for plots.

%       Calculate and plot residuals for the model obtained by IV4:
resid(z,i2);  %Plots the result of residual analyssi
 resid(z,i2,'fr') % presents the result in the frequncy domain.
e = resid(z,i2);  % Just computes the residuals
figure,plot(e(:,1,[])), % e(:,1,[]) selects all samples, output #1 and no inputs
title('Residuals, IV4 method')

%       Now compute second order ARMAX model:
am2 = armax(z,[2 2 2 1]);

%       Now compute second order BOX-JENKINS model:
bj2 = bj(z,[2 2 2 2 1]);

%       Compare the frequency functions with those obtained by the
%       PEM method and the spectral analysis method:


figure,bode(am2,bj2,m,GS)

%       % To compare the noise spectra:
figure,bode(am2('noise'),bj2('n'),m('n'),GS('n'))

%       Plot residuals for the ARMAX and BOX-JENKINS models
resid(z,am2); 
resid(z,bj2); 

%       Compare the BJ model with the ARMAX model and with the 
%       default PEM model:
figure,compare(z,bj2,am2,m);

%       Now we'll compute the poles and zeros of the ARMAX and
%       BOX-JENKINS models
figure, zpplot(am2,bj2),  % Press a key for plot

% 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 (Akaike's  AIC-criterion):
[aic(am2), aic(bj2)]




%       Since we generated the data, we enjoy the luxury of comparing
%       the model to the true system:
 
  % Press any key for plots of comparisons.

figure,bode(m0,am2)
%       % The responses practically coincide
figure,bode(m0('noise'),am2('noise'))
figure,zpplot(m0,am2)