Code covered by the BSD License

# Developing Models from Experimental Data using System Identification Toolbox

24 Apr 2007 (Updated )

Describe linear and nonlinear data-centric modeling approaches.

webinar_walk_through.m
```%% Developing Models from Experimental Data using System Identification Toolbox
%
% On April 05, 2007 a webinar was presenting showcasing the capabilities
% of System Identification Toolbox 7.0 for linear and nonlinear model
% estimation. The title of the webinar was: "Developing Models from
% Experimental Data Using System Identification Toolbox".
%
% This file summarizes the examples used in this presentation.

% For accessing the webinar, visit the product page at:
% http://www.mathworks.com/products/sysid/
% or directly: http://www.mathworks.com/wbnr30317
%

%{
Author(s): Rajiv Singh
\$Date: 2007/04/23 \$
%}

%%
% Note: All the ideas presented in the webinar have been explained and
% illustrated in-depth in various demos of System Identification Toolbox.
% The demos cover several aspects of data handling, spectral, linear and
% nonlinear (parametric) model estimation and model usage. It is strongly
% advised that you look at these demos, which may be accessed from the
% product web page.

%% Contents
% 1. Linear Model Estimation: using time and frequency domain data
% 2. Nonlinear Black Box Model Estimation
% 3. Nonlinear Grey Box Model Estimation

%% Linear Model Estimation
% Consider a chemical production plant for which you are interested in
% modeling the dynamics between a certain voltage input and the rate of
% production of the chemical.
%
% Depending upon the type of experiment performed, we may have data
% available as either time signals or as frequency response.
%
% Load the hypothetical time and frequency domain data for the chemical
% plant.

%%
% Variable u represents the input (voltage) and y contains the output
% values (production rate). This is "time domain" data, which means that
% the measured signals represent values of variables as a function of time.
% The data is assumed to be sampled at 0.08 seconds.

%%
% Note that data processing and model estimation can be done in the GUI
% (System Identification Tool). However, these actions can also be
% performed in the MATLAB Command Window. In System Identification Toolbox,
% I/O data is stored in a "data object", called IDDATA object:

data = iddata(y, u, 0.08, 'InputName', 'Voltage', 'OutputName', 'Production');

%%
% 'data' represents time-domain data with sampling interval of 0.08 seconds
% and input name = 'Voltage' and output name = 'Production'. The object
% 'data' can now be used for pre-processing and estimation.

%%
% Detrend data:
datad = detrend(data,0); % remove mean

%%
% Select a portion of detrended data for estimation and another portion for
% validation (it is a good idea to use a data different from estimation
% data for validating the quality of the model):

eData = datad(1:500);   % select first 500 samples of datad for estimation
vData = datad(501:end); % select the second-half of data for validation

%%
% We can use eData to estimate various linear models such as ARX/ARMAX
% models, state-space models and continuous-time low order transfer
% functions ("process models").
%
% It is commonly advised to run the function ADVICE on data to get some
% information on its characteristics, such as inter-sampling behavior,
% indication of nonlinearity and feedback, and persistence of excitation:

%%
% For linear models, it is possible to inspect data to obtain suitable
% model orders and input delay. Use functions such as DELAYEST, ARXSTRUC
% and SELSTRUC for this purpose. You may also use the state-space
% estimator function N4SID to both determine the optimal number of states
% of the model as well as estimate the corresponding model. For example, we
% can search for the best model order in the range 1 to 10 by calling
% N4SID as follows:

model = n4sid(eData,1:10); % use a range of values for model order

%%
% This launches an interactive figure window that lets you choose the
% optimal model order (a 3rd order model is suggested in our case).

% For more details on model order determination, see the demo called
% "Model Structure Selection: Determining Model Order and Input Delay"
% (iddemo3.m)

%%
% Estimate a linear ARX model (Ay=Bu+e), with na = 3, nb=2, nk=1:
model = arx(eData, [3 2 1]);

%%
% If you did not know what orders to use, you could try several orders
% together and pick the best ones.

%%
% Estimate a linear Box-Jenkins model (y = B/Fu + C/De):
model = bj(eData,'nb',2,'nf',3,'nc',2,'nd',2);

%%
% The advantage of using BJ model is that it is able to capture non-white
% (colored) disturbances using a more complex noise model (=C/De).

%%
% You can separate out the measured and noise components of the model by
% sub-referencing as follows:

measurement_model = model('meas'); % model from measured input to output
noise_model = model('noise');      % model from input noise to output

%%
% You may also use other polynomial forms such as ARMAX and OE for
% estimation. Here, let us estimate a 3rd order state-space model:

model1 = n4sid(eData,3); % a subspace method of estimating state-space models

model2 = pem(eData,3); % model2 is usually more accurate than model1.

%%
% PEM is the main estimation function. Based on the inputs you provide,
% this function can estimate a variety of linear and nonlinear models. You
% can also use this function to refine the parameters of existing models
% using data. Refer to documentation on PEM for more details.

%%
% Now we estimate a continuous-time model ("Process Model") with 3 poles
% (containing a real pole and an underdamped pair), one zero and a delay
% element:

initial_model = idproc('p3uzd'); % IDPROC object represents a process model

model = pem(eData, initial_model); % estimate the parameters of initial_model and return the results in 'model'.

%%
% These estimations can be performed using the GUI (System Identification
% Tool) as well. Some of these estimations were shown during the webinar.

%%
% You may run PRESENT to view a description of the model (equation,
% parameter values, quality of model etc). PRESENT also shows the standard
% deviations on parameter values, whenever available:

present(model)  % view model equation with parameter uncertainty information

%%
% When using the GUI, right click on the model's icon in the model board to
% open a window that shows the model equation and reproduction steps. Press
% the "Present" button on this window to view the detailed information in
% the MATLAB Command Window. The most detailed information about the model
% is revealed by running GET on the model:
get(model)

%%
% GET provides a way of querying the model and obtain in-depth information
% on its parameter values, uncertainty values, covariance matrix,
% estimation related information (final prediction value, loss,
% optimization termination criteria, estimation data characteristics, etc),
% and algorithm settings.

%% Estimating Linear Models using Frequency Domain Data
% Estimation with frequency domain data uses the same syntax as that for
% estimation with time-domain data. The frequency responses data (FRF) is
% represented by an IDFRD object which encapsulates the complex frequency
% response and the corresponding frequency vector. The frequency response
% data is sometimes measured during actual experimentation using spectrum
% analyzers. If you have time domain data, you may use functions such as
% ETFE, SPA and SPAFDR to estimate the frequency response. For example, we
% may do:

G0 = spa(eData); % estimate non-parametric frequency response
get(G0) % inspect properties of G0 - an IDFRD object

bode(G0, 'sd', 1, 'fill') % plot response with 1 s.d. confidence region

%%
% For the above production process, the frequency response data was also
% collected as Amplitude and Phase (deg.) values for each frequency. We can
% package this data using an IDFRD object:

Phase = Phase*pi/180; % convert phase data into radians
Resp = Amp.*exp(i*Phase);

FRFdata =  idfrd(Resp, freq, 0.08,...
'InputName', 'Voltage', 'OutputName', 'Production');

%%
% Use BODE to view the frequency response data:
bode(FRFdata)

%%
% Now, FRFdata can be used in place of eData for estimation. Note that a
% noise model cannot be estimated using frequency response data. Hence
% functions ARMAX and BJ would not work. However, ARX, output-error (OE)
% and state-space models, as well as process models (IDPROC) can be
% estimated.

modelf1 = arx(FRFdata, [3 2 1]); % estimate ARX model
modelf2 = oe(FRFdata, [2 3 1]);  % estimate an OE model
modelf3 = pem(FRFdata, 4); %4th order state-space model
modelf4 = pem(FRFdata, idproc('p3dzu')); % process model with 3 poles, a zero and delay

%%
% When/Why do we use frequency domain data? Here are some possible reasons:
% 1. The obvious case is when the data was collected in frequency domain,
% such using spectrum analyzers.
%
% 2. Sometimes, we may convert time domain signals into FRF data to achieve
% data compression, such as convert large time-based data into a much smaller
% number of frequency response (FRF) values by using SPA or SPAFDR.
%
% 3. We may convert data from time domain into frequency domain to control
% resolution in frequency ranges of interest. For example, while time
% domain signals are required to be uniformly sampled, we may use a finer
% frequency grid on FRF data around the regions of interest such as
% resonant peaks and use a coarser grid elsewhere. Look up the function
% SPAFDR for achieving frequency-dependent resolution.
%
% 4. It is possible to obtain continuous-time models more easily using
% frequency domain data by setting the sampling interval to zero. For
% example, if we want a continuous-time state-space model, we may do:

FRFdatac = FRFdata;
FRFdatac.Ts = 0;
modelf4 = pem(FRFdatac, 4); %4th order CT state-space model

%% Inspecting Model Behavior and Validating Its Quality
% You can view model's responses in time (step, impulse, simulation using
% chosen input) and frequency (bode, nyquist) domains, and view the pole-zero
% maps. When viewing any of these responses, it is possible to view the
% corresponding response uncertainty regions too. For example, if we want
% to view the model's step response with 1 s.d. confidence region, we can
% do:

step(model, 'sd', 1)
step(modelf4, 'sd', 1)

%%
% Similarly, we may run BODE or PZMAP commands. Monte-Carlo simulations of
% the model may be performed using SIMSD. For example, we can generate 10
% responses consistent with the model's covariance information for a given
% sinusoidal input as:

u = sin(2*pi*5*[0:100]'/100); % 5 Hz signal
U = iddata([],u,0.1,'uname','Voltage'); % data packaged as iddata object to add sampling information (10 Hz).

simsd(model, U) % 10 random model simulations

%%
% Viewing the model responses along with the uncertainty information
% reveals some information about the quality of the model. You can
% ascertain if the uncertainty is too large (may imply bad data or too high
% a model order) or if the resonant frequencies (just as an example) are
% misplaced in the estimated model. Unnecessarily high orders may show up
% as pole-zero pairs with largely overlapping uncertainty regions in the
% pole-zero plot. Overlaps are indicative of an "effective" pole-zero
% cancellation.

% In addition to the visual inspection of response plots, the toolbox
% provides two dedicated tools for validation of estimated models - COMPARE
% and RESID. COMPARE facilitates the most obvious validation technique -
% comparing the model's response against the measured data. After estimation,
% it is advisable to view how good the model is in approximating the
% measured output data. To do so, run:

compare(eData, model) % plot shows 90% fit for IDPROC model

%%
% When comparing responses (simulated vs measured) it is advisable to use
% an independent ("validation") data set, that is different from the
% estimation data. Here, we use vData for validation:

compare(vData, model) % comparison against validation data set

%%
% In the GUI, the response comparisons are viewed by clicking on the Model
% output checkbox.

%%
% The second important function that facilitates validation is RESID. RESID
% plots the correlation of residual prediction error with itself and with
% the input signals. A good model should result in mostly "white"
% (uncorrelated) residuals. Practically speaking, the correlation of
% residuals should be "small", that is, within a 1 s.d. region around zero.
% This plot can be viewed as:

resid(vData, model) %plot residuals

%%
% What if there are significant values in the correlation plots? It may
% indicate that the model has not captured all the dynamics; we may need a
% higher order model. It may also denote uncaptured noise dynamics. For
% example, for the case of IDPROC models, we may also estimate a second
% order noise model by doing:

model2 = pem(eData, initial_model, 'disturbancemodel', 'arma2');

compare(vData, model, model2); % compare model and model2 simultaneously
resid(vData, model2) % view residuals corresponding to model2

%%
% The residuals plot shows that the prediction errors are mostly
% uncorrelated for model2.

%%
% The models obtained using frequency response data can be validated and
% used just as in the case of time-domain models.

compare(vData, modelf1, modelf3, modelf4)   % validate models using time-domain data
compare(FRFdata, modelf1, modelf3, modelf4) % validate models using frequency response data

%% Inserting Estimated Model into Simulink
% You can bring an estimated linear model into Simulink using the IDMODEL
% block from the toolbox's block library. To see the contents of the block
% library, type SLIDENT or access it from Simulink Library Browser:

slident

%%
open_system('ExampleModel') % an example model

sim('ExampleModel') % simulate model in Simulink

%%
% You can convert estimated models into LTI objects used by Control System
% Toolbox, Model Predictive Control Toolbox or Robust Control Toolbox:
sys1 = tf(model); % transfer function
sys2 = ss(model('measured')) %convert only the measured component of the model into State-space LTI

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Nonlinear Black Box Model Estimation

%%
% The above scheme of data estimation and comparison of models extends
% seamlessly to the nonlinear models too. However, nonlinear models require
% careful setup (good initial guesses, more data etc). Also, quite
% frequently, the models may not validate against data that have been
% obtained under different operating conditions (than the estimation data),
% or have a different range or amplitude. In spite of these constraints, the
% toolbox does let us evaluate models with various kinds of nonlinearities
% with minimal knowledge of the underlying dynamics. We can create
% nonlinear ARX models and models with static I/O nonlinearities
% (Hammerstein-Wiener models) and use a variety of non-linear
% approximators, such as wavelets, neural networks (requires Neural Network
% Toolbox), sigmoid and user-defined "custom" networks, piecewise-linear
% nonlinearity, saturation and deadzone elements etc.

%%
% Let us consider modeling the effect of air bypass valve's position on the
% idle speed deviation in the IC engine of a car. The air valve position is
% controlled by a voltage. Hence an experiment could be conceived of,
% wherein the voltage is changed (input) and the resulting change in
% idle speed (output) is recorded. The dynamics are likely to be quite
% nonlinear. In absence of a first-principles based model (which, in
% practice, may be hard to obtain or even unnecessary depending upon the
% application at hand), we will try to estimate "black box" models using
% data.

%%
% First, load the experimental data:
load icEngine % contains two variables - u and y

% (this data is at <matlabroot>/toolbox/ident/iddemos/icEngine.mat)

%%
% If working in the GUI, import u as the input and y as the output data and
% set the sampling interval to 0.04 seconds. If working in Command Window,
% package the data into an IDDATA object:
z = iddata(y,u, 0.04,'InputName','Voltage','OutputName',...
'Speed','InputUnit','V','OutputUnit','RPM/100');

%%
% We will split the data and use the first half for estimation and the
% latter half for validation.
eData = z(1:1000);     % estimation data
vData = z(1001:end);   % validation data

%% Nonlinear ARX models
% The nonlinear ARX models are represented by IDNLARX objects in the
% toolbox. They represent a flexible autoregressive structure capable of
% approximating arbitrary nonlinearities using various estimators such as
% wavelet network, tree partition etc.
%
% We need to make two main choices: the number of regressors to use and the
% type of nonlinearity to use. Regressors are specified by entering the
% "order" values, just as we do in the case of linear ARX model estimation.
% For the current example, we set na = 4 (four output regressors: y(t-1),
% y(t-2), .. y(t-4)), nk = 3 (three sample delay in input) and nb = 2 (two
% input regressors: u(t-3) and u(t-4)).

%%
% For nonlinearity, let us choose a wavelet network estimator. The
% advantage of this estimator is that it is quite versatile and flexible
% and thus suited for approximating arbitrary nonlinearities. Also,
% estimation of a wavelet network based model is non-iterative in nature.
% This makes the estimation fast and less prone to bad choices for initial
% values of its parameters. To estimate the model, we use the NLARX
% estimation function:

Order = [4 2 3];
Nonlinearity = wavenet; % number of wavelet units to use in the network is
% determined automatically, by default
M1 = nlarx(eData, Order, Nonlinearity);

%%
% If working in the GUI, open the "Nonlinear Models" window from the
% estimation popu menu. Then set the regressor orders and the nonlinearity
% type in the Regressors and the Model Properties tabs respectively.

%%
% Compare the model's (M1) response against both the estimation data and
% the validation data:

compare(eData,M1) % 81% fit to estimation data
compare(vData,M1) % 70% fit to estimation data

%%
% To see some details on the model, just type its name:

M1 % also, run GET(M1) to see detailed info on model

%%
% We can see what regressors were used and what the nonlinearity was. The
% display indicates that only 1 wavelet unit was used in the model M1.
% If we want to use more wavelet units in the network (say 5), we can set
% those explicitly:

Nonlinearity = wavenet('NumberOfUnits',5); % wavenet with 5 units
M2 = nlarx(eData, Order, Nonlinearity);

% Type M2.Nonlinearity (or just M2.nl) to inspect the properties of the
% wavelet network.

M2.nl

% We can view the parameter values for the wavelet network.

%%
% Note that nonlinear ARX models employ a parallel configuration - a linear
% block is put in parallel (in most cases) to a nonlinear block. A way to
% make sense of this configuration is to realize that the nonlinear block
% is capturing the elements of the overall response that could not be
% captured by the linear element alone. If we do not want the linear
% parallel block, we can remove it by using the LinearTerm property of the
% wavelet during estimation:

Nonlinearity = wavenet('NumberOfUnits',5,'LinearTerm','off'); % no linear block
M3 = nlarx(eData, Order, Nonlinearity);

%%
% We can similarly try other types of nonlinearities such as tree-partition
% or sigmoid networks.

M4 = nlarx(eData, Order, treepartition); % using tree-partition nonlinearity
M5 = nlarx(eData, Order, sigmoidnet, 'trace', 'on'); % using sigmoid network

%%
% Note that unlike the wavenet and treepartition estimation, estimation
% using sigmoidnet is an iterative process. We turned the trace on to view
% the display of the estimation process in the Command Window.

%%
% Any model can be iteratively refined using NLARX or PEM. In doing so, we
% may also set the focus of estimation to 'Simulation' rather than the
% default 'Prediction'. Note that with focus = 'prediction', the estimator
% tries to create the best results for 1-step ahead prediction of the
% output, while with focus = 'simulation', the aim of the estimator is to
% best fit the provided data over its measured time range.

% For example, if we want to refine the parameter values of a wavelet
% network based model (M2) we may do:

% (this estimation will go slower)
M2r = pem(eData, M2, 'focus', 'sim', 'maxiter', 5, 'trace', 'on'); % (or do: nlarx(eData, M2, ...) )

%%
% Note that an edge case of nonlinear ARX modeling is to not use the
% nonlinear block at all. This would make the model linear and similar in
% structure to a linear ARX model. This configuration is achieved by using
% LINEAR as a type of nonlinearity. It is quite useful to estimate such
% models for two main reasons: (1) This will give us a reference point to
% determine how nonlinear the model is, and (2) This would allow us to
% create linear models with customized regressors (discussed later), which,
% in some ways, are the simplest types of nonlinear models.

Mlin = nlarx(eData, Order, linear); % or, nlarx(eData, Order, []) -> "[]" denotes absence of nonlinearity

%% Creating Nonlinear ARX Models with Custom Regressors
% When we specify Order ([4 2 3] here), we are basically asking the
% estimator to automatically create simple regressors which are delayed
% version of the I/O variables. These regressors are referred to as
% "standard" regressors. In our example, they were: Speed(t-1),..,
% Speed(t-4), Voltage(t-3) and Voltage(t-4). However, frequently we want to
% employ inputs such as Speed(t-1)^4 or Speed(t-3)*Voltage(t-2). Such
% regressors can be created using CUSTOMREG function and can be injected
% into the nonlinear model as additional "custom" regressors. For simple
% formulas such as those specified above, we may specify them directly as
% strings during estimation.

Mc1 = nlarx(eData, Order, linear, ...
'customreg', {'Speed(t-4)^2', 'Speed(t-5)^2', 'Speed(t-3)*Voltage(t-2)'});

%%
% Display of this model shows both standard and custom regressors. We can
% remove standard regressors (set na = nb = 0) and use only custom regressors:

Mc2 = nlarx(eData, [0 0 0], linear, ...
'customreg', {'Speed(t-4)^2', 'Speed(t-5)^2', 'Speed(t-3)*Voltage(t-2)'});

%%
% Also, we may use custom regressors with models containing a nonlinearity
% (such as wavenet):

Mc3 = nlarx(eData, Order, wavenet, ...
'customreg', {'Speed(t-4)^2', 'Speed(t-5)^2', 'Speed(t-3)*Voltage(t-2)'});

%%
% Final note on custom regressors: the above mechanism provides a way for
% users to inject their knowledge about the linearity into the black box
% model of the system. For example, if a user knows that he/she has
% quadratic nonlinearity in his/her system, he/she may create regressors such
% as Speed(t-1)^2, Voltage(t-2)^2, Voltage(t-10)^2 etc and add them to the
% model. For more details, see demo idnlbbdemo_customreg, titled
% "Nonlinear ARX Models with Custom Regressors".

%% Inspecting Nonlinear ARX Models
% We can compare the fit obtained by these models against the estimation
% and validation data using COMPARE. We can also view their step responses.

compare(vData, M1, M2, Mlin, Mc1, Mc3) % just a subset of models; M1 appears to be the best

step(M1, M2, M3) %step responses of the models

%%
% We can inspect the shape and amplitude of the nonlinearity,  which is the
% output of the nonlinear block as a function of input regressors. For
% this, use the function PLOT. For example, if we want to plot and compare
% the nonlinearities of models M2 (contains wavenet), M5 (contains
% sigmoidnet) and Mc1 (linear with custom regressors), we may do:

plot(M2, M5, Mc1)

%%
% The plot, by default, shows amplitude of nonlinearities in the three models
% as a function of regressors Speed(t-1) and Speed(t-2). The other
% regressors are fixed to constant (mid-point) values. The plot is
% essentially "flat" for M2 and Mc1 indicating that the nonlinearity is
% quite weak in them. We can choose a different set of regressors to plot
% or change the cross-section (mid-point) values of the regressors not
% shown in the plot. Rotate the 3-D plot to inspect its shape.

%% Hammerstein-Wiener Models
% Hammerstein-Wiener models employ a series configuration wherein the
% nonlinear elements are placed in series with a linear block. Unlike the
% nonlinear ARX models, the nonlinearities in the Hammerstein-Wiener models
% are static (no time delay). We may think of these models as essentially
% linear models with static distortions (such as saturation or dead-zone)
% at their inputs and outputs. Sub-configurations, such as input
% nonlinearity only (Hammerstein models) or output nonlinearity only
% (Wiener models) are also possible. These models are represented by IDNLHW
% objects in the toolbox.

%%
% For estimation of Hammerstein-Wiener models, we need the following
% information:
%
% 1. Type of nonlinearity at each input and output channel. It is possible
% to specify no nonlinearities on one or more of the I/O channels.
%
% 2. Order of the linear block: We need to specify the number of numerator
% and denominator terms in the embedded linear model's transfer function.
% The linear model may be thought of as a output-error type (OE) linear
% model which is specified by coefficients nb, nf and nk to generate the
% transfer function: y = B/F u(t-nk).

%%
% Let us estimate a Hammerstein model using estimation function called
% NLHW. If estimating in the GUI, set the model structure to
% Hammerstein-Wiener in the Nonlinear Models window and use the I/O
% nonlinearities and Linear Model tabs to configure the model.

InputNL = pwlinear; %piecewise linear nonlinearity on input
OutputNL = unitgain; % no nonlinearity on output (o, we could use OutputNL = []
Order = [2 4 3];  % = [nb nf nk]

N1 = nlhw(eData, Order, InputNL, OutputNL); % Hammerstein model created as an IDNLHW object.

%%
% All Hammerstein-Wiener models are iterative in nature. We may turn trace
% display 'on' to view the progress of estimation. Note that in the above
% estimation, OutputNL was set to unit gain, which in effect, removes the
% nonlinearity from the output channel.

%%
% Similarly, we may estimate a Wiener model (no input nonlinearity) with
% output nonlinearity describe by a wavelet network using NLHW:

InputNL = []; %or, InputNL = unitgain
OutputNL = wavenet;
N2 = nlhw(eData, Order, InputNL, OutputNL, 'trace','on'); % Wiener model

%%
% The embedded linear block can be extracted out using the "LinearModel"
% property of the IDNLHW model:

linMod = N2.LinearModel

%%
% As before, we can invoke COMPARE to compare the models' responses against
% measured data, or use STEP/SIM to simulate the model.

compare(eData, N1, N2) % comparison against estimation data
compare(vData, N1, N2) % comparison against validation data

%%
% The plots show that the fit to estimation data is quite good, but the
% validation results are not that good. Let us estimate a full nonlinear
% model with nonlinearities at both input and output ports (both pwlinear)
% and also increase the order of the linear model:

InputNL = pwlinear;
OutputNL = pwlinear;
Order = [4 5 3];
N3 = nlhw(eData, Order, InputNL, OutputNL, 'trace','on', 'maxiter', 40); % Hammerstein-Wiener model

%%
% Hammerstein-Wiener models also support specification of saturation and
% deadzone as the types of I/O nonlinearities. For example, if we knew that
% our system is affected by voltage saturation, we could set the type of
% nonlinearity to saturation. Also, suppose we suspect presence of output
% nonlinearity, but we do not know its nature. Then, we can use a generic,
% flexible nonlinearity approximator such as pwlinear or wavenet to specify
% the output nonlinearity:

InputNL = saturation;
OutputNL = wavenet;
Order = [4 5 3];
N4 = nlhw(eData, Order, InputNL, OutputNL, 'trace','on', 'maxiter', 40);

%%
% For reference, let us also estimate a model with no nonlinearities:
Nlin =  nlhw(eData, [4 5 3], [] , []) % linear model

%%
% Compare the responses of a mix of nonlinear models:
compare(vData, M1, M2, Mlin, N1, N3, Nlin)

%% Inspecting Hammerstein-Wiener Models
% For closer visual inspection of Hammerstein-Wiener models, the function
% PLOT is available. Using this function, we can inspect the model on a
% block-by-block basis. We can view the shapes of the input and output
% nonlinearities and plot the time and frequency domain responses (step,
% impulse, bode and pole-zero map) of the embedded linear block.

plot(N1, N3, N4, Nlin) % inspect N1, N3, N4 and Nlin models

% In the resulting  figure, click on a block  (U_NL, Linear Block or Y_NL)
% to view the corresponding plot.

%%
% Some more examples of nonlinear ARX model and Hammerstein-Wiener model
% estimation are available in demos. See demos- "A Two Tank System -
% Single-Input Single-Output Nonlinear ARX and Hammerstein-Wiener Models"
% (idnlbbdemo_siso) and "A Motorized Camera - Multi-Input Multi-Output
% Nonlinear ARX and Hammerstein-Wiener Models" (idnlbbdemo_mimo).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimation of Nonlinear Grey Box Models
% Grey box models represent a class of those models whose structures are
% well known but the values of their parameters (various physical or
% empirical constants) may not be known exactly. It is not uncommon to
% use a spring-mass-damper combination to model (in lumped parameters) a
% dynamic system. In this model, the value of the damping coefficient and
% the stiffness constant may still need to be estimated from experiments.
% Grey box models, whether linear or nonlinear, (IDGREY and IDNLGREY models
% in System Identification Toolbox) support this activity.

%%
% In the webinar, only the IDNLGREY model was shown. The IDGREY (linear
% counterpart) usage scenario is simpler in many ways and uses a very
% similar way to handle the model file and the parameters. The example
% considered was that of the friction dynamics between two sliding bodies.
% The example was taken from the demo - "Dry Friction Between Two Bodies :
% Mex-File Modeling Using Multiple Experiments" (idnlgreydemo7.m).

% Model of dry friction between two bodies is defined by following parameters:
%{
m: Mass of the moving body [kg]
k: Sliding friction force coefficient [kg/s]
e: Distance-related dry friction parameter [m]
f: Force-related dry friction parameter [N]
%}

%%
% For the equations, see demo idnlgreydemo7. We
% basically have a set of nonlinear ordinary differential equations. Given
% these equations of motion, we can bring them in MATLAB using IDNLGREY
% object and set it up for simulation, estimation and analysis. For this to
% happen, the equations must first be written in "state-space form", that
% is, as a set of first-order nonlinear ODEs. We may write these equations
% in an M file or a MEX file (C/C++, Fortran). The resulting file is called
% an ODE file. The ODE file must return the values of state-derivatives and
% model output(s) as a function of time, inputs, states and parameter
% values. The model may be either continuous or discrete in time. If the
% model is discrete in time, the returned state "derivatives" are simply
% the values of the states at the next sampling instant (x(t+1)).

%%
% For our example, we can transform the equations of motion into
% state-space form by choosing and using the following three states:
% Position and Velocity of the moving body, and Dry Friction force between
% the bodies. Then, the resulting ODE file can be written as shown in
% files twobodies_m (M file implementation) or twobodies_c (MEX file
% implementation). In this example, we use the M-file version. Writing the
% ODE files in M is simpler but the MEX version may be faster.

%%
% Once we have the ODE file available, all we need to do is define an
% IDNLGREY model and call estimation and simulation commands on it. The
% IDNLGREY model is configured as shown below:

%%
% Step 1: designate the ODE file to be used
FileName = 'twobodies_m'; % ODE file

%%
% Step 2: define number of inputs, outputs and states
Order = [1 1 3]; % SISO system with three states

%%
% Step 3: Set initial values of Parameters
Parameters = [380; 2200; 0.00012; 1900];

%%
% Step 4: Specify sampling interval:
Ts = 0; % this is a continuous-time model

%%
% Create Model Object
% Now, construct the IDNLGREY model by calling its constructor:
FrictionModel = idnlgrey(FileName, Order, Parameters, [], Ts, ...
'Name', 'Two body system',                      ...
'InputName', 'Exogenous force',                 ...
'InputUnit', 'N',                               ...
'OutputName', 'Position of moving body',        ...
'OutputUnit', 'm',                              ...
'TimeUnit', 's')

%%
% The model  'FrictionModel' is ready to be used for simulation and analysis:

sim(FrictionModel,sin(2*pi*5*[0:200]'/100)) %simulate response to given input

step(FrictionModel,0.6); %step response for time 0 to 0.6 seconds

%%
% It is possible to configure which parameters and initial states to
% estimate and which to fix to their known values. By default, all
% parameters are treated as estimatable quantities and initial states are
% not estimated. You can change these default in the IDNLGREY model and
% configure each parameter and initial state separately.

%%
% Next, let us estimate parameters to fit experimental data. For this let
% us first load some data:

%%
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'twobodiesdata'));

%%
% Package data and its properties (sampling interval and I/O names) in a
% IDDATA object:
data = iddata(y1,u1,0.005,...
'InputName','Exogenous force','OutputName','Position of moving body');

%%
% Before estimation, let us compare the response of the original model
% against measured output data:
compare(data,FrictionModel) % fit is not good

%%
% Now, let us estimate parameters to fit data better. For this, we use the
% function PEM:
FrictionModel2 = pem(data, FrictionModel, 'Trace', 'Full');

present(FrictionModel2) % view model with parameter values and their std deviations

%%
% Finally, we compare the estimated model's response against measured data:
compare(data,FrictionModel2) % good fit

%%
% We can bring the estimated model into Simulink using IDNLGREY block that
% is a part of toolbox's block library.

inputData = data(:,[],:); %extract input signals for simulation of model in Simulink
open_system('Friction_Model')
sim('Friction_Model')

%% Conclusions
% In this walk through demo, we saw how to estimate a variety of linear and
% nonlinear black box models. The syntax for creation and comparison of
% these models are very similar. The toolbox offers a uniform and
% comprehensive way of creating and comparing linear and nonlinear models.
% Note that nonlinear block box models are inherently discrete in time. The
% nonlinear ARX models can be created for time-series also (no inputs). It
% is possible to linearize these models using the following functions:
%
% LINAPP - Linear approximation of nonlinear model for a given input range.
% LINEARIZE (used to be called LINTAN) - Small signal tangent linearization
% of nonlinear models.

%%
% We also saw how to estimate parameters of user-defined models using data.
% This is a useful way of estimating coefficients of ordinary differential
% equations. If you have legacy C or Fortran code describing the behavior
% of your system, you can bring this description directly into MATLAB by
% first creating a MEX wrapper for it (see the MEX template file:
% IDNLGREY_MODEL_TEMPLATE.c) and then creating an IDNLGREY model for it.
%
% For DAE systems or hybrid systems (mix of continuous and discrete time
% components), it is better to use Simulink. If creating the system in
% from data.

%% Using Identified Models
% The linear and nonlinear models created using data can
% be used for simulation and other analysis. You can bring them into
% Simulink using dedicated model blocks available as part of System
% Identification block library. Then, these models could be used as a
% component of larger systems. In most cases, the linear and nonlinear
% black box models can also be built using Real-Time Workshop(R) software
% which facilitates easy deployment of such models in real applications.

```