www.gusucode.com > ident_featured 案例代码 matlab源码程序 > ident_featured/iddemosl.m

    %% Estimating Continuous-Time Models using Simulink Data
%
% This example illustrates how models simulated in Simulink(R) can be 
% identified using System Identification Toolbox(TM). The example describes
% how to deal with continuous-time systems and delays, as well as the
% importance of the intersample behavior of the input.

% Copyright 1986-2016 The MathWorks, Inc. 

%%
if exist('start_simulink','file')~=2
    disp('This example requires Simulink.')
    return
end

%% Acquiring Simulation Data from a Simulink Model 
% Consider the system described by the following Simulink model:
open_system('iddemsl1')
set_param('iddemsl1/Random Number','seed','0')

%%
% The red part is the system, the blue part is the controller and
% the reference signal is a swept sinusoid (a chirp signal).
% The data sample time is set to 0.5 seconds.
%
% This system can be represented using an |idpoly| structure:
m0 = idpoly(1,0.1,1,1,[1 0.5],'Ts',0,'InputDelay',1,'NoiseVariance',0.01)

%%
% Let us simulate the model |iddemsl1| and save the data in an |iddata|
% object:
sim('iddemsl1')
dat1e = iddata(y,u,0.5); % The IDDATA object

%%
% Let us do a second simulation of the mode for validation purposes.
set_param('iddemsl1/Random Number','seed','13')
sim('iddemsl1')
dat1v = iddata(y,u,0.5);

%%
% Let us have a peek at the estimation data obtained during the first
% simulation:
plot(dat1e)

%% Estimating Discrete Models Using the Simulation Data
% Let us begin by evaluating a default-order discrete model to gain some
% preliminary insight into the data characteristics:
m1 = n4sid(dat1e, 'best')     % A default order model

%%
% Check how well the model reproduces the validation data
compare(dat1v,m1)

%%
% As observed, the validation data is predicted well by the model. To
% investigate more into the data characteristics, let us inspect the
% non-parametric impulse response computed using |dat1e| where the negative
% lags for anaysis are automatically determined:
ImpModel = impulseest(dat1e,[],'negative'); 
clf
h = impulseplot(ImpModel);
showConfidence(h,3)

%%
% |ImpModel| is an FIR model whose order (no. of coefficients) are
% automatically determined. We also choose to analyze feedback effects by
% computing impulse response for negative delays. Influences from negative
% lags are not all insignificant. This is due to the regulator (output
% feedback). This means that the impulse response estimate cannot be used
% to determine the time delay. Instead, build several low order ARX-models
% with different delays and find out the  best fit:
V = arxstruc(dat1e,dat1v,struc(1:2,1:2,1:10));
nn = selstruc(V,0) %delay is the third element of nn

%%
% The delay is determined to 3 lags. (This is correct: the deadtime of
% 1 second gives two lag-delays, and the ZOH-block another one.)
% The corresponding ARX-model can also be computed, as follows:
m2 = arx(dat1e,nn)
compare(dat1v,m1,m2);

%% Refining the Estimation
% The two models |m1| and |m2| behave similarly in simulation. Let us now
% try and fine-tune orders and delays. Fix the delay to 2 (which coupled
% with a lack of feedthrough gives a net delay of 3 samples) and find a
% default order state-space model with that delay:
m3 = n4sid(dat1e,'best','InputDelay',2,'Feedthrough',false);
% Refinement for prediction error minimization using pem (could also use
% |ssest|)
m3 = pem(dat1e, m3); 

%%
% Let as look at the estimated system matrix
m3.a  % the A-matrix of the resulting model

%%
% A third order dynamics is automatically chosen, which together with the
% 2 "extra" delays gives a 5th order state space model.
%
% It is always advisable not to blindly rely upon automatic order choices.
% They are influenced by random errors. A good way is to look at the model's
% zeros and poles, along with confidence regions:
clf
h = iopzplot(m3);
showConfidence(h,2) % Confidence region corresponding to 2 standard deviations

%%
% Clearly the two poles/zeros at the unit circle seem to cancel, indicating 
% that a first order dynamics might be sufficient. Using this information,
% let us do a  new first-order estimation:
m4 = ssest(dat1e,1,'Feedthrough',false,'InputDelay',2,'Ts',dat1e.Ts);
compare(dat1v,m4,m3,m1)

%%
% The |compare| plot shows that the simple first order model |m4| gives a
% very good description of the data. Thus we shall select this model as our
% final result.

%% Converting Discrete Model to Continuous-Time (LTI) 
% Convert this model to continuous time, and represent it in transfer
% function form:
mc = d2c(m4);
idtf(mc)

%% 
% A good description of the system has been obtained, as displayed above.

%% Estimating Continuous-Time Model Directly
% The continuous time model can also be estimated directly. The discrete
% model |m4| has 2 sample input delay which represents a 1 second delay. We
% use the |ssest| command for this estimation:
m5 = ssest(dat1e,1,'Feedthrough',false,'InputDelay',1);
present(m5)

%% Uncertainty Analysis
% The parameters of model |m5| exhibit high levels of uncertainty even
% though the model fits the data 87%. This is because the model uses more
% parameters than absolutely required leading to a loss of uniqueness in
% parameter estimates. To view the true effect of uncertainty in the model,
% there are two possible approaches:
%%
% 
% # View the uncertainty as confidence bounds on model's response rather
% than on the parameters. 
% # Estimate the model in canonical form.
% 
% Let use try both approaches. First we estimate the model in canonical
% form.
m5Canon = ssest(dat1e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical');
present(m5Canon)

%%
% |m5Canon| uses a canonical parameterization of the model. It fits the
% estimation data as good as the model |m5|. It shows small uncertainties
% in the values of its parameters giving an evidence of its reliability.
% However, as we saw for |m5|, a large uncertainty does not necessarily
% mean a "bad" model. To ascertain the quality of these models, let use
% view their responses in time and frequency domains with confidence
% regions corresponding to 3 standard deviations. We also plot the original
% system |m0| for comparison.

%%
% The bode plot.
clf
opt = bodeoptions;
opt.FreqScale = 'linear';
h = bodeplot(m0,m5,m5Canon,opt);
showConfidence(h,3)
legend show

%%
% The step plot.
clf
showConfidence(stepplot(m0,m5,m5Canon),3)
legend show

%%
% The uncertainty bounds for the two models are virtually identical. We
% can similarly generate pole-zero map (|iopzplot|) and Nyquist plot
% (|nyquistplot|) with confidence regions for these models.
idtf(m5)

%% Accounting for Intersample Behavior in Continuous-Time Estimation
% When comparing continuous time models computed from sampled data, it
% is important to consider the intersample behavior of the input signal.
% In the example so far, the input to the system was piece-wise constant, 
% due to the Zero-order-Hold (zoh) circuit in the controller. Now remove
% this circuit, and consider a truly continuous system. The input and output
% signals are still sampled a 2 Hz, and everything else is the same:
open_system('iddemsl3')
sim('iddemsl3')
dat2e = iddata(y,u,0.5);

%%
% Discrete time models will still do well on these data, since when they
% are adjusted to the measurements, they will incorporate the sampling
% properties, and intersample input behavior (for the current input).
% However, when building continuous time models, knowing the intersample 
% properties is essential. First build a model just as for the ZOH case:
m6 = ssest(dat2e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical'); 
idtf(m6)

%%
% Let us compare the estimated model (|m6|) against the true model (|m0|):
step(m6,m0)  % Compare with true system

%%
% The agreement is now not so good. We may, however, include in the
% data object information about the input. As an approximation, let it
% be described as piecewise linear (First-order-hold, FOH) between the
% sampling instants. This information is then used by the estimator for
% proper sampling:
dat2e.Intersample = 'foh';
m7 = ssest(dat2e,1,'Feedthrough',false,'InputDelay',1,'Form','canonical');  % new estimation with correct intersample behavior
idtf(m7)

%%
% Let us look at the step response comparison again:
step(m7,m0)  % Compare with true system

%%
% This model (|m7|) gives a much better result than |m6|. This concludes
% this example.
bdclose('iddemsl1');
bdclose('iddemsl3');