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

    %% Identifying Nonlinear ARX and Hammerstein-Wiener Models Using Measured Data 
% This example shows how to identify single-input-single-output (SISO)
% nonlinear black box models using measured input-output data. The example
% uses measured data from a two-tank system to explore various model
% structures and identification choices.

% Copyright 2005-2015 The MathWorks, Inc.

%% The Input-Output Data Set
% In this example, the data variables stored in the |twotankdata.mat| file
% will be used. It contains 3000 input-output data samples of a two-tank
% system generated using a sampling interval of 0.2 seconds. The input u(t)
% is the voltage [V] applied to a pump, which generates an inflow to the
% upper tank. A rather small hole at the bottom of this upper tank yields
% an outflow that goes into the lower tank, and the output y(t) of the two
% tank system is then the liquid level [m] of the lower tank. We create an
% IDDATA object |z| to encapsulate the loaded data:
%
% This data set is also used in the nonlinear grey box example "Two Tank
% System: C MEX-File Modeling of Time-Continuous SISO System" where
% physical equations describing the system are assumed. This example is
% about black box models, thus no physical equation is assumed.

load twotankdata
z = iddata(y, u, 0.2, 'Name', 'Two tank system');  
     
%% Plotting the Data
% Let us plot the whole 3000 samples of data, for the time ranging from 0
% to 600 seconds.
%
% The input is a multi-step signal with various levels. 

plot(z)

%% Splitting the Data
% Split this data set into 3 subsets of equal size, each containing 1000
% samples.
%
% According to the data plot, the input and output signals of z2 (the
% middle third of the plot) is similar to the signals of z1, but those of 
% z3 are quite different, slightly beyond the input and output ranges of
% z1.
%
% We shall estimate models with z1, and evaluate them on z2 and z3.
% In some sense, z2 allows to evaluate the interpolation ability of the
% models, and z3 allows to evaluate their extrapolation ability to some
% extent.
z1 = z(1:1000); z2 = z(1001:2000); z3 = z(2001:3000);

%% Model Order Selection
% The first step in estimating nonlinear black box models is to choose model 
% orders. For an IDNLARX model, Orders = [na nb nk], defining the numbers
% of past outputs, past inputs and the input delay used to predict the
% current output.
%
% Usually model orders are chosen by trials and errors. Sometimes
% it is useful to try linear ARX models first in order to get hints about
% model orders. So let us first try the tool for linear ARX model orders
% selection.
%
V = arxstruc(z1,z2,struc(1:5, 1:5,1:5));
nn = selstruc(V,'aic') % selection of nn=[na nb nk] by Akaike's information criterion (AIC)

%%
% The AIC criterion has selected nn = [na nb nk] = [5 1 3], meaning that,
% in the selected ARX model structure, y(t) is predicted by the 6 regressors 
% y(t-1),y(t-1),...,y(t-5) and u(t-3). We will try to use these values in 
% nonlinear models.

%% Nonlinear ARX (IDNLARX) Models
% In an IDNLARX model, the model output is a nonlinear function of 
% regressors, like model_output(t) = F(y(t-1),y(t-1),...,y(t-5), u(t-3)).
%
% IDNLARX models are typically estimated with the syntax:
%
%   M = NLARX(Data, Orders, Nonlinearity).
%
% It is similar to the linear ARX command, with the additional argument 
% Nonlinearity specifying the type of nonlinearity estimator.

%%
% To estimate an IDNLARX model, after the choice of model orders, we 
% should choose the nonlinearity estimator to be used. Let us first try the
% WAVENET estimator.
%
mw1 = nlarx(z1,[5 1 3], wavenet);

%% 
% The number of units (wavenet functions) in the WAVENET estimator has been
% automatically chosen. It can be accessed as shown below.
mw1.Nonlinearity.NumberOfUnits

%%
% Examine the result by comparing the output simulated with the estimated
% model and the output in the estimation data z1:
%
compare(z1,mw1); 

%% Playing with Model Properties
%
% The first model was not quite satisfactory, so let us try to modify some
% properties of the model. Instead of letting the estimation algorithm
% automatically choose the number of units in the WAVENET estimator, this
% number can be explicitly specified.
%
mw2 = nlarx(z1,[5 1 3], wavenet('NumberOfUnits',8));

%%
% Default InputName and OutputName have been used.
%
mw2.InputName
mw2.OutputName

%%
% The regressors of the IDNLARX model can be viewed with the command getreg.
% The regressors are made from |mw2.InputName|, |mw2.OutputName| and time delays.
%
getreg(mw2)

%% Nonlinear Regressors of IDNLARX Model
% An important property of IDNLARX models is NonlinearRegressors. The output
% of an IDNLARX model is a function of regressors. Usually this function has a 
% linear block and a nonlinear block. The model output is the sum of the
% outputs of the two blocks. In order to reduce model complexity, a subset of
% regressors can be chosen to enter the nonlinear block. These regressors
% are referred to as "nonlinear regressors".
%
% Nonlinear regressors can be specified as a vector of regressor indices 
% in the order of the regressors returned by GETREG. For the example of mw2,
% the entire list of regressors is displayed below, and
% NonlinearRegressors=[1 2 6] corresponds to y1(t-1), y1(t-2) and u1(t-3)
% in the following list.
%
getreg(mw2)

%%
% Nonlinear regressors can be indicated by the Property-Value pair as in
% the following example. 
% Notice the short-hand 'nlreg'='NonlinearRegressors'.
%
%mw3 = nlarx(z1,[5 1 3], wavenet, 'NonlinearRegressors', [1 2 6]);
mw3 = nlarx(z1,[5 1 3], wavenet, 'nlreg', [1 2 6]);  % using short-hand notation
getreg(mw3, 'nonlinear')  % get nonlinear regressors


%% 
% Instead of regressor indices, some keywords can also be used to
% indicate nonlinear regressors: 'input' for regressors consisting of
% delayed inputs, and 'output' for regressors consisting of delayed 
% outputs. In the following example, 'input' indicates the only input 
% regressor u1(t-3).
%
mw4 = nlarx(z1,[5 1 3],wavenet,'nlreg','input');
mw4.nlreg    % 'nlreg' is the short-hand of 'NonlinearRegressors' 
getreg(mw4,'nonlinear')  % get nonlinear regressor

%% Automatic Search for Nonlinear Regressors
% To automatically try all the possible subsets of regressors as nonlinear 
% regressors and to choose the one resulting in the lowest simulation error, 
% use the value NonlinearRegressors='search'. This exhaustive search usually 
% takes a long time. It is recommended to turn on the iteration display
% when using the search option.
%
% This example may take some time.
opt = nlarxOptions('Display','on');
mws = nlarx(z1,[5 1 3], wavenet, 'nlreg', 'search', opt);

%%
% The result of the search can be accessed in mws.nlreg
% The resulting value mws.nlreg=[4 6] means that the 4th and 6th 
% regressors are nonlinear in the following list of regressors.
% In other words, y1(t-4) and u1(t-3) are used as nonlinear
% regressors.
%
mws.nlreg
getreg(mws)

%% Evaluating Estimated Models
% Different models, both linear and nonlinear, can be compared together
% in the same COMPARE command.
%
mlin = arx(z1,[5 1 3]);

compare(z1,mlin,mw2,mws); 

%%
% Now let us see how mws behaves on the data set z2.
%
compare(z2,mws);

%%
% and on the data set z3.
%
compare(z3,mws); 

%%
% Function PLOT may be used to view the nonlinearity responses of various
% IDNLARX models.
plot(mw2,mws) % plot nonlinearity response as a function of selected regressors

%% Nonlinear ARX Model with SIGMOIDNET Nonlinearity Estimator
% At the place of WAVENET, other nonlinearity estimators can be used.
% Try the SIGMOIDNET estimator.
%
ms1 = nlarx(z1,[5 1 3], sigmoidnet('NumberOfUnits', 8));
compare(z1,ms1) 

%%
% Now evaluate the new model on the data set z2.
%
compare(z2,ms1);  

%%
% and on the data set z3.
%
compare(z3,ms1);  

%% Other Nonlinearity Estimators in IDNLARX Model
% CUSTOMNET is similar to SIGMOIDNET, but instead of the sigmoid unit
% function, the user can define other unit functions in MATLAB files. This 
% example uses the gaussunit function defined in the function gaussunit.m.
%
% TREEPARTITION is another nonlinearity estimator.
%
mc1 = nlarx(z1,[5 1 3], customnet(@gaussunit),'nlreg',[3 5 6]);
mt1 = nlarx(z1,[5 1 3], treepartition);

%% Using the Network Object from Neural Network Toolbox(TM) (NNTB)
% Neural networks created with the aid of the NNTB can also be used. This
% requires NNTB license.

%%
% Here, we will create a single-output network with an unknown number of
% inputs (denote by input size of zero in the network object). The number
% of inputs is set to the number of regressors (as determined by model
% orders) in the IDNLARX model during the estimation process
% automatically.

%%
% Compare responses of models |mc1| and |mt1| to data |z1|
compare(z1, mc1, mt1)

if exist('feedforwardnet','file')==2 % check if Neural Network Toolbox is available
    % This example is run only if the NNTB license is available.
    ff = feedforwardnet([5 20]);
    ff.layers{2}.transferFcn = 'radbas';
    ff.trainParam.epochs = 50;
    mn1 = nlarx(z1,[5 1 3], neuralnet(ff)); % estimation with neuralnet nonlinearity
    compare(z1, mn1) % compare fit to estimation data
end

%% Hammerstein-Wiener (IDNLHW) Models
% A Hammerstein-Wiener model is composed of up to 3 blocks: a linear
% transfer function block is preceded by a nonlinear static block and/or 
% followed by another nonlinear static block. It is called a Wiener model
% if the first nonlinear static block is absent, and a Hammerstein model 
% if the second nonlinear static block is absent.
%
% IDNLHW models are typically estimated with the syntax:
%
%   M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).
%
% where Orders = [nb bf nk] specifies the orders and delay of the linear
% transfer function, InputNonlinearity and OutputNonlinearity specify the
% nonlinearity estimators for the two nonlinear blocks. The linear output
% error (OE) model corresponds to the case of trivial identity
% nonlinearities.

%% Hammerstein-Wiener Model with the Piecewise Linear Nonlinearity Estimator
% The PWLINEAR (piecewise linear) nonlinearity estimator is useful in
% general IDNLHW models.
%
% Notice that, in Orders = [nb nf nk], nf specifies the number of poles of
% the linear transfer function, somewhat related to the |na| order of the
% IDNLARX model
mhw1 = nlhw(z1, [1 5 3], pwlinear, pwlinear);

%%
% The result can be evaluated with COMPARE as before.
compare(z1,mhw1)  

%%
% Model properties can be visualized by the PLOT command.
% 
% Click on the block-diagram to choose the view of the input nonlinearity
% (UNL), the linear block or the output nonlinearity (YNL). 
%
% For the linear block, it is possible to choose the type of plots within 
% Step response, Bode plot, Impulse response and Pole-zero map.
plot(mhw1)

%%
% The break points of the two piecewise linear estimators can be examined
% (notice that the short-hands 'u'='Input', 'y'='Output', 'nl'='Nonlinearity'
% can be used). 
% This is a two row matrix, the first row corresponds to the
% input and the second row to the output of the PWLINEAR estimator.
mhw1.unl.BreakPoints

%% Hammerstein-Wiener Model with Saturation and Dead Zone Nonlinearities
% The SATURATION and DEADZONE estimators can be used if such nonlinearities
% are supposed to be present in the system.
%
mhw2 = nlhw(z1, [1 5 3], deadzone, saturation);  
compare(z1,mhw2)  

%%
% The absence of a nonlinearity is indicated by the UNITGAIN object, or by
% the empty "[]" value. Abbreviated names of nonlinearity
% estimators can also be used.
%
mhw3 = nlhw(z1, [1 5 3], 'dead',unitgain);  
mhw4 = nlhw(z1, [1 5 3], [],'satur');  

%%
% The limit values of DEADZONE and SATURATION can be respectively examined.
%
mhw2.unl.ZeroInterval
mhw2.ynl.LinearInterval

%%
% The initial guess of ZeroInterval for DEADZONE or LinearInterval for 
% SATURATION can be indicated in the estimators:
mhw5 = nlhw(z1, [1 5 3], deadzone([0.1 0.2]), saturation([-1 1]));  


%% Hammerstein-Wiener Model - Specifying More Properties
% The estimation algorithm options can be specified using the NLHWOPTIONS
% command.
%  
opt = nlhwOptions();
opt.SearchMethod = 'gna';
opt.SearchOption.MaxIter = 3;
mhw6 = nlhw(z1, [1 5 3], deadzone, saturation, opt);

%%
% An already estimated model can be refined by more estimation iterations
%
% Evaluated on the estimation data z1, mhw7 should have a better fit than mhw6.
%
mhw7 = nlhw(z1, mhw6, opt);
compare(z1, mhw6, mhw7)

%% Hammerstein-Wiener Model - Use Other Nonlinearity Estimators
% The nonlinearity estimators PWLINEAR, SATURATION and DEADZONE have been
% mainly designed for use in IDNLHW models. They can only estimate single
% variable nonlinearities. The other more general nonlinearity estimators, 
% SIGMOIDNET, CUSTOMNET and WAVENET, can also be used in IDNLHW models.
% However, the non-differentiable estimators TREEPARTITION and NEURALNET
% are not applicable, because the estimation of IDNLHW models needs the
% derivatives of the nonlinearity estimators in iterative algorithms.
% Type "idprops idnlestimators" for a summary of nonlinearity estimators.