www.gusucode.com > ident_featured 案例代码 matlab源码程序 > ident_featured/iddemo2.m
%% Comparison of Various Model Identification Methods
% This example shows several identification methods available in System
% Identification Toolbox(TM). We begin by simulating experimental data and
% use several estimation techniques to estimate models from the data. The
% following estimation routines are illustrated in this example: |spa|,
% |ssest|, |tfest|, |arx|, |oe|, |armax| and |bj|.
% Copyright 1986-2015 The MathWorks, Inc.
%% System Description
% A noise corrupted linear system can be described by:
%
% |y = G u + H e|
%
% where |y| is the output and |u| is the input, while |e| denotes the
% unmeasured (white) noise source. |G| is the system's transfer function
% and |H| gives the description of the noise disturbance.
%
% There are many ways to estimate |G| and |H|. Essentially they correspond
% to different ways of parameterizing these functions.
%% Defining a Model
% System Identification Toolbox provides users with the option of
% simulating data as would have been obtained from a physical process. Let
% us simulate the data from an IDPOLY model with a given set of
% coefficients.
B = [0 1 0.5];
A = [1 -1.5 0.7];
m0 = idpoly(A,B,[1 -1 0.2],'Ts',0.25,'Variable','q^-1'); % The sample time is 0.25 s.
%%
% The third argument |[1 -1 0.2]| gives a characterization of the
% disturbances that act on the system. Execution of these commands produces
% the following discrete-time idpoly model:
m0
%%
% Here |q| denotes the shift operator so that A(q)y(t) is really short for
% y(t) - 1.5 y(t-1) + 0.7 y(t-2). The display of model |m0| shows that it
% is an ARMAX model.
%
% This particular model structure is known as an ARMAX model, where AR
% (Autoregressive) refers to the A-polynomial, MA (Moving average) to the
% noise C-polynomial and X to the "eXtra" input B(q)u(t).
%
% In terms of the general transfer functions |G| and |H|, this model
% corresponds to a parameterization:
%
% |G(q) = B(q)/A(q)| and |H(q) = C(q/A(q)|, with common denominators
%% Simulating the Model
% We generate an input signal |u| and simulate the response of the model to
% these inputs. The |idinput| command can be used to create an input signal
% to the model and the |iddata| command will package the signal into a
% suitable format. The |sim| command can then be used to simulate the
% output as shown below:
prevRng = rng(12,'v5normal');
u = idinput(350,'rbs'); %Generates a random binary signal of length 350
u = iddata([],u,0.25); %Creates an IDDATA object. The sample time is 0.25 sec.
y = sim(m0,u,'noise') %Simulates the model's response for this
%input with added white Gaussian noise e
%according to the model description m0
rng(prevRng);
%%
% One aspect to note is that the signals |u| and |y| are both IDDATA
% objects. Next we collect this input-output data to form a single iddata
% object.
%
z = [y,u];
%%
% To get information on the data object (which now incorporates both the
% input and output data samples), just enter its name at the MATLAB(R)
% command window:
z
%%
% To plot the first 100 values of the input |u| and output |y|, use the
% |plot| command on the iddata object:
h = gcf; set(h,'DefaultLegendLocation','best')
h.Position = [100 100 780 520];
plot(z(1:100));
%%
% It is good practice to use only a portion of the data for estimation
% purposes, |ze| and save another part to validate the estimated models
% later on:
ze = z(1:200);
zv = z(201:350);
%% Performing Spectral Analysis
%
% Now that we have obtained the simulated data, we can estimate models and
% make comparisons. Let us start with spectral analysis. We invoke the
% |spa| command which returns a spectral analysis
% estimate of the frequency function and the noise spectrum.
GS = spa(ze);
%%
% The result of spectral analysis is a complex-valued function of
% frequency: the frequency response. It is packaged into an object
% called |IDFRD| object (Identified Frequency Response Data):
GS
%%
% To plot the frequency response it is customary to use the Bode plot, with
% the |bode| or |bodeplot| command:
clf
h = bodeplot(GS); % bodeplot returns a plot handle, which bode does not
ax = axis; axis([0.1 10 ax(3:4)])
%%
% The estimate |GS| is uncertain, since it is formed from noise corrupted
% data. The spectral analysis method provides also its own assessment of
% this uncertainty. To display the uncertainty (say for example 3 standard
% deviations) we can use the |showConfidence| command on the plot handle
% |h| returned by the previous |bodeplot| command.
%
showConfidence(h,3)
%%
% The plot says that although the nominal estimate of the frequency
% response (blue line) is not necessarily accurate, there is a 99.7%
% probability (3 standard deviations of the normal distribution) that the
% true response is within the shaded (light-blue) region.
%% Estimating Parametric State Space Models
% Next let us estimate a default (without us specifying a particular model
% structure) linear model. It will be computed as a
% state-space model using a prediction error method. We use the
% |ssest| function in this case:
%
m = ssest(ze) % The order of the model will be chosen automatically
%%
% At this point the state-space model matrices do not tell us very much. We
% could compare the frequency response of the model |m| with that of the
% spectral analysis estimate |GS|. Note that |GS| is a discrete-time model
% while |m| is a continuous-time model. It is possible to use |ssest| to
% compute discrete time models too.
bodeplot(m,GS)
ax = axis; axis([0.1 10 ax(3:4)])
%%
% We note that the two frequency responses are close, despite
% the fact that they have been estimated in very different ways.
%% Estimating Simple Transfer Functions
% There is a wide variety of linear model structures, all corresponding to
% linear difference or differential equations describing the relation
% between u and y. The different structures all correspond to various ways
% of modeling the noise influence. The simplest of these models are the
% transfer function and ARX models. A transfer function takes the form:
%
% Y(s) = NUM(s)/DEN(s) U(s) + E(s)
%
% where NUM(s) and DEN(s) are the numerator and denominator polynomials,
% and Y(s), U(s) and E(s) are the Laplace Transforms of the output, input
% and error signals (y(t), u(t) and e(t)) respectively. NUM and DEN can be
% parameterized by their orders which represent the number of zeros and
% poles. For a given number of poles and zeros, we can estimate a transfer
% function using the |tfest| command:
mtf = tfest(ze, 2, 2) % transfer function with 2 zeros and 2 poles
%%
% AN ARX model is represented as:
% A(q)y(t) = B(q)u(t) + e(t),
%
% or in long-hand,
%
% y(t) + a_1 y(t-1) + ... +y_na y(t-na)
% = b_1 u(t-nk) + ...+ b_nb u(t-nk-nb-1) + e(t)
%
% This model is linear in coefficients. The coefficients a_1, a_2, ...,
% b_1, b_2, .. can be computed efficiently using a least squares
% estimation technique. An estimate of an ARX model with a prescribed
% structure - two poles, one zero and a single lag on the input is obtained
% as ([na nb nk] = [2 2 1]):
mx = arx(ze,[2 2 1])
%% Comparing the Performance of Esti Models
% Are the models (|mtf|, |mx|) better or worse than the default state-space
% model |m|? One way to find out is to simulate each model (noise free)
% with the input from the validation data |zv| and compare the
% corresponding simulated outputs with the measured output in the set |zv|:
compare(zv,m,mtf,mx)
%%
% The fit is the percent of the output variation that is explained by the
% model. Clearly |m| is the better model, although |mtf| comes close. A
% discrete-time transfer function can be estimate using either the |tfest|
% or the |oe| command. The former creates an IDTF model while the latter
% creates an IDPOLY model of Output-Error structure (y(t) = B(q)/F(q)u(t) +
% e(t)), but they are mathematically equivalent.
md1 = tfest(ze,2,2,'Ts',0.25) % two poles and 2 zeros (counted as roots of polynomials in z^-1)
md2 = oe(ze,[2 2 1])
compare(zv, md1, md2)
%%
% The models |md1| and |md2| deliver identical fits to the data.
%% Residual Analysis
%
% A further way to gain insight into the quality of a model is to compute
% the "residuals": i.e. that part |e| in |y = Gu + He| that could not be
% explained by the model. Ideally, these should be uncorrelated with the
% input and also mutually uncorrelated. The residuals are computed and
% their correlation properties are displayed by the |resid| command. This
% function can be used to evaluate the residues both in the time and
% frequency domains. First let us obtain the residuals for the Output-Error
% model in the time domain:
resid(zv,md2) % plots the result of residual analysis
%%
% We see that the cross correlation between residuals and input lies in the
% confidence region, indicating that there is no
% significant correlation. The estimate of the dynamics |G| should then be
% considered as adequate. However, the (auto) correlation of |e| is
% significant, so |e| cannot be seen as white noise. This means that the
% noise model |H| is not adequate.
%% Estimating ARMAX and Box-Jenkins Models
% Let us now compute a second order ARMAX model and a second order
% Box-Jenkins model. The ARMAX model has the same noise characteristics as
% the simulated model |m0| and the Box-Jenkins model allows a more general
% noise description.
%
am2 = armax(ze,[2 2 2 1]) % 2nd order ARMAX model
%%
bj2 = bj(ze,[2 2 2 2 1]) % 2nd order BOX-JENKINS model
%%
%
%% Comparing Estimated Models - Simulation and Prediction Behavior
% Now that we have estimated so many different models, let us make some
% model comparisons. This can be done by comparing the simulated outputs as
% before:
%
clf
compare(zv,am2,md2,bj2,mx,mtf,m)
%%
% We can also compare how well the various estimated models are able to
% predict the output, say, 1 step ahead:
%
compare(zv,am2,md2,bj2,mx,mtf,m,1)
%%
% The residual analysis of model |bj2| shows that the prediction error is
% devoid of any information - it is not auto-correlated and also
% uncorrelated with the input. This shows that the extra dynamic elements
% of a Box-Jenkins model (the C and D polymials) were able to capture the
% noise dynamics.
resid(zv,bj2)
%% Comparing Frequency Functions
%
% In order to compare the frequency functions for the
% generated models we use again the |bodeplot| command:
%
clf
opt = bodeoptions; opt.PhaseMatching = 'on';
w = linspace(0.1,4*pi,100);
bodeplot(GS,m,mx,mtf,md2,am2,bj2,w,opt);
legend({'SPA','SS','ARX','TF','OE','ARMAX','BJ'},'Location','West');
%%
% The noise spectra of the estimated models can also be analyzed. For
% example here we compare the noise spectra the ARMAX and the Box-Jenkins
% models with the state-space and the Spectral Analysis models. For this,
% we use the |spectrum| command:
spectrum(GS,m,mx,am2,bj2,w)
legend('SPA','SS','ARX','ARMAX','BJ');
%% Comparing Estimated Models with the True System
%
% Here we validate the estimated models against the true system |m0| that
% we used to simulate the input and output data. Let us compare the
% frequency functions of the ARMAX model with the true system.
%
bode(m0,am2,w)
%%
%
% The responses practically coincide. The noise spectra can also be
% compared.
%
spectrum(m0,am2,w)
%%
%
% Let us also examine the pole-zero plot:
%
h = iopzplot(m0,am2);
%%
%
% It can be seen that the poles and the zeros of the true system (blue)
% and the ARMAX model (green) are very close.
%%
% We can also evaluate the uncertainty of the zeros and poles. To plot
% confidence regions around the estimated poles and zeros corresponding to
% 3 standard deviations, use |showConfidence| or turn on the "Confidence
% Region" characteristic from the plot's context (right-click) menu.
%
showConfidence(h,3)
%
%%
% We see that the true, blue zeros and poles are well inside the green
% uncertainty regions.
%% Additional Information
%
% For more information on the different estimation techniques that are
% available with System Identification Toolbox(TM) visit the
% <http://www.mathworks.com/products/sysid/ System Identification Toolbox>
% product information page.