Motorized Camera - Multi-Input Multi-Output Nonlinear ARX and Hammerstein-Wiener Models

The Measured Data Set

The data saved in the file motorizedcamera.mat will be used in this example. It contains 188 data samples, collected from a motorized camera with a sample time of 0.02 seconds. The input vector u(t) is composed of 6 variables: the 3 translation velocity components in the orthogonal X-Y-Z coordinate system fixed to the camera [m/s], and the 3 rotation velocity components around the X-Y-Z axes [rad/s]. The output vector y(t) contains 2 variables: the position (in pixels) of a point which is the image taken by the camera of a fixed point in the 3D space. We create an IDDATA object z to hold the loaded data:

load motorizedcamera
z = iddata(y, u, 0.02, 'Name', 'Motorized Camera', 'TimeUnit', 's');

Nonlinear ARX (IDNLARX) Model - Picking Regressors

Let us first try nonlinear ARX models. Two important elements need to be chosen: the model regressors and the nonlinear mapping functions.

The regressors are simple formulas based on time-delayed I/O variables, the simplest case being the variables lagged by a small set of consecutive values. For example, if "u" in the name of the input variable, and "y" the name of the output variables, then an example regressor set can be { y(t-1), y(t-2), u(t), u(t-1), u(t-2) }, where "t" denotes the time variable. Another example involving polynomial terms could be {y(t-2)^4, y(t-2)*u(t-1), u(t-4)^2}. More complex, arbitrary formulas in the delayed variables are also possible.

The easiest way of specifying linear regressors is by using an "orders matrix". This matrix takes the form NN = [na nb nk] and it specifies by how many lags each output (na) and each input variable (nb, nk) are delayed. This is the same idea that is used when estimating the linear ARX models (see ARX, IDPOLY). For example, NN = [2 3 4] implies the regressor set {y(t-2), u(t-4), u(t-5), u(t-6)}. In the general case of models with NY outputs and NU inputs, NN is a matrix with NY rows and NY+2*NU rows.

Nonlinear ARX (IDNLARX) Model - Preliminary Estimation Using Wavelet Network

To start, let us choose the orders NN = [na nb nk] = [ones(2,2), ones(2,6), ones(2,6)]. It means that each output variable is predicted by all the output and input variables, each being delayed by 1 sample. The model equation can be written as y_i(t) = F_i(y1(t-1), y2(t-1), u1(t-1), u2(t-1), u3(t-1)), i = 1, 2. Let us choose a Wavelet Network (idWaveletNetwork object) as the nonlinear mapping function for both outputs. The function NLARX estimates the parameters of the Nonlinear ARX model.

NN = [ones(2,2), ones(2,6), ones(2,6)]; % the orders
mw1 = nlarx(z, NN, idWaveletNetwork);

Examine the result by comparing the output simulated with the estimated model and the output in the measured data z:

compare(z,mw1)

Nonlinear ARX Model - Trying Higher Orders

Let us check if we can do better with higher orders. Note that when identifying models using basis function expansions to express the nonlinearity, the number of model parameters can exceed the number of data samples. In such cases, some estimation metrics such as Noise Variance and Final Prediction Error (FPE) cannot be determined reliably. For the current example, we turn off the warning informing us about this limitation.

ws = warning('off','Ident:estimation:NparGTNsamp');
%
mw2 = nlarx(z, [ones(2,2), 2*ones(2,6), ones(2,6)], idWaveletNetwork)
compare(z,mw2)

mw2 =

<strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong>
  Inputs: u1, u2, u3, u4, u5, u6
  Outputs: y1, y2

Regressors:
  Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6

Output functions:
  Output 1:  Wavelet network with 27 units
  Output 2:  Wavelet network with 25 units

Sample time: 0.02 seconds

Status:                                                      
Estimated using NLARX on time domain data "Motorized Camera".
Fit to estimation data: [99.22;99.15]% (prediction focus)    
FPE: 0.1262, MSE: 0.4453

The second model mw2 is pretty good. So let us keep this choice of model orders in the following examples.

nanbnk = [ones(2,2), 2*ones(2,6), ones(2,6)]; % final orders

To view the regressor set implied by this choice of orders, use the GETREG command:

getreg(mw2)

ans =

  14x1 cell array

    {'y1(t-1)'}
    {'y2(t-1)'}
    {'u1(t-1)'}
    {'u1(t-2)'}
    {'u2(t-1)'}
    {'u2(t-2)'}
    {'u3(t-1)'}
    {'u3(t-2)'}
    {'u4(t-1)'}
    {'u4(t-2)'}
    {'u5(t-1)'}
    {'u5(t-2)'}
    {'u6(t-1)'}
    {'u6(t-2)'}

The numbers of units ("wavelets") of the two idWaveletNetwork functions have been automatically chosen by the estimation algorithm. These numbers are displayed below.

mw2.OutputFcn(1).NonlinearFcn.NumberOfUnits
mw2.OutputFcn(2).NonlinearFcn.NumberOfUnits

ans =

    27


ans =

    25

Nonlinear ARX Model - Adjusting Number of Units of Nonlinear Functions

The number of units in the idWaveletNetwork function can be explicitly specified instead of being automatically chosen by the estimation algorithm. Suppose we want to use 10 units for the first nonlinear mapping function, and 5 for the second one (note that the model for each output uses its own mapping function; the array of all mapping functions is stored in the model's "OutputFcn" property).

Fcn1 = idWaveletNetwork(10); % output function for the first output
Fcn2 = idWaveletNetwork(5);  % output function for the second output
mw3 = nlarx(z, nanbnk, [Fcn1; Fcn2])

mw3 =

<strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong>
  Inputs: u1, u2, u3, u4, u5, u6
  Outputs: y1, y2

Regressors:
  Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6

Output functions:
  Output 1:  Wavelet network with 10 units
  Output 2:  Wavelet network with 5 units

Sample time: 0.02 seconds

Status:                                                      
Estimated using NLARX on time domain data "Motorized Camera".
Fit to estimation data: [99.01;98.89]% (prediction focus)    
FPE: 0.2273, MSE: 0.7443

Nonlinear ARX Model - Trying Other Nonlinear Mapping Functions

In place of the idWaveletNetwork function, other nonlinear functions can also be used. Let us try the idTreePartition for both outputs.

mt1 = nlarx(z, nanbnk, idTreePartition);

In the above call, a tree partition mapping function with default configuration is used for both the outputs.

Similarly, we can use a Sigmoid Network (idSigmoidNetwork object) as the mapping function. Estimation options such as maximum iterations (MaxIterations) and iteration display can be specified using the nlarxOptions command.

opt = nlarxOptions('Display','on');
opt.SearchOptions.MaxIterations = 5;
ms1 = nlarx(z, nanbnk, idSigmoidNetwork, opt);

This calling syntax for NLARX is very similar to the ones used before - specifying the data, the orders and the nonlinear mapping functions as their primary input arguments. However, in order to modify the estimation options from their default values, we constructed an option set, opt, using the nlarxOptions command and passed it to the NLARX command as an additional input argument.

Nonlinear ARX Model with Mixed Nonlinear Functions

It is possible to use different nonlinear functions on different output channels in the same model. Suppose we want to use a tree partition function to describe the first output and use a wavelet network for the second output. The model estimation is shown below. The third input argument (the nonlinear mapping function) is now an array of two different functions.

Fcn1 = idTreePartition;
Fcn2 = idWaveletNetwork;
mtw = nlarx(z, nanbnk, [Fcn1; Fcn2]);

Sometimes the function that maps the regressors to the model output need not be nonlinear; it could simply be a weighted sum of the regressor vectors, plus an optional offset. This is similar to the linear ARX models (except for the offset term). The absence of nonlinearity in an output channel can be indicated by choosing a idLinear function. The following example means that, in model_output(t) = F(y(t-1), u(t-1), u(t-2)), the function F is composed of a linear component for the first output, and a nonlinear component (idSigmoidNetwork) for the second output.

Fcn1 = idLinear;
Fcn2 = idSigmoidNetwork(2);
opt.Display = 'off'; % do not show estimation progress anymore
mls = nlarx(z, nanbnk, [Fcn1; Fcn2], opt)

mls =

<strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong>
  Inputs: u1, u2, u3, u4, u5, u6
  Outputs: y1, y2

Regressors:
  Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6

Output functions:
  Output 1:  Linear with offset
  Output 2:  Sigmoid network with 2 units

Sample time: 0.02 seconds

Status:                                                      
Estimated using NLARX on time domain data "Motorized Camera".
Fit to estimation data: [98.72;98.8]% (prediction focus)     
FPE: 0.5481, MSE: 1.041

There is no nonlinearity in the model for the first output. It is not, however, purely linear because of the presence of an offset term.

disp(mls.OutputFcn(1))

<strong>Linear Function</strong>
Inputs: y1(t-1), y2(t-1), u1(t-1), u1(t-2), u2(t-1), u2(t-2), u3(t-1), u3(t-2), u4(t-1), u4(t-2), u5(t-1), u5(t-2), u6(t-1), u6(t-2)
Output: y1(t)

 Linear Function: initialized to [48.3 -32.2 -0.229 -0.0706 -0.113 -0.0516 -0.182 0.297 0.199 -0.133 -0.337 0.583 -0.448 0.167]
 Output Offset: initialized to -9.03e-06

    LinearFcn: '<Linear function parameters>'
       Offset: '<Offset parameters>'

We have the option of making it purely linear by choosing to not use the offset term. This choice can be made by fixing the offset to a zero value before estimation.

Fcn1.Offset.Value = 0;
Fcn1.Offset.Free = false;
mlsNoOffset = nlarx(z, nanbnk, [Fcn1; Fcn2], opt);

Using Regression Algorithms from Statistics and Machine Learning Toolbox

If you have access to the Statistics and Machine Learning Toolbox™ you can also use Gaussian process (GP) functions and bagged or boosted tree ensembles to create your Nonlinear ARX models. The GPs are represented by the idGaussianProcess object, while regression tree ensembles are represented by the idTreeEnsemble object. GPs use squared-exponential kernel by default. Other kernel types may be specified when creating the object. In the example here, we use a 'ARDMatern52' kernel based GP for the first output, and a bagged tree ensemble for the second output.

if exist('fitrgp','file')==2
    warning('off','stats:lasso:MaxIterReached');
    NL1 = idGaussianProcess('ardmatern52');
    NL2 = idTreeEnsemble;
    % estimate model
    mML = nlarx(z, nanbnk, [NL1; NL2])
end

mML =

<strong>Nonlinear ARX model with 2 outputs and 6 inputs</strong>
  Inputs: u1, u2, u3, u4, u5, u6
  Outputs: y1, y2

Regressors:
  Linear regressors in variables y1, y2, u1, u2, u3, u4, u5, u6

Output functions:
  Output 1:  Gaussian process function using a ARDMatern52 kernel
  Output 2:  Bagged Regression Tree Ensemble

Sample time: 0.02 seconds

Status:                                                      
Estimated using NLARX on time domain data "Motorized Camera".
Fit to estimation data: [99.29;98.44]% (prediction focus)    
MSE: 0.968

Inspection of Estimation Results

Various models can be compared in the same COMPARE command. For example, compare the simulated outputs of the models mw2, ms1, mls, and mlsNoOffset to the measured output in dataset "z".

compare(z, mw2, ms1, mls, mlsNoOffset)

Function plot may be used to view the nonlinear mapping function "shape" for various models.

plot(mt1,mtw,ms1,mls)

Note that the control panel on the right hand side of the plot allows for regressor selection and configuration.

Other functions such as resid, predict and pe may be used on the estimated models in the same way as in the case of linear models.

Hammerstein-Wiener (IDNLHW) Model - Preliminary Estimation

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 using the syntax:

M = NLHW(Data, Orders, InputNonlinearity, OutputNonlinearity).

where Orders = [nb nf nk] specifies the orders and delay of the linear transfer function. InputNonlinearity and OutputNonlinearity specify the nonlinear functions for the two nonlinear blocks. The linear output error (OE) model corresponds to the case of trivial (unit gain mapping) nonlinearities.

Estimation of Hammerstein Model (No Output Nonlinearity)

Let us choose Orders = [nb bf nk] = [ones(2,6), ones(2,6), ones(2,6)]. It means that, in the linear block, each output is the sum of 6 first order transfer functions driven by the 6 inputs. Try a Hammerstein model (no output nonlinearity) with the input nonlinearity described by a piecewise linear function. Notice that the entered single idPiecewiseLinear object is automatically expanded to all the 6 input channels. idUnitGain indicates absence of nonlinearity.

opt = nlhwOptions();
opt.SearchOptions.MaxIterations = 15;
NN = [ones(2,6), ones(2,6), ones(2,6)];
InputNL  = idPiecewiseLinear;  % input nonlinearity
OutputNL = idUnitGain;  % output nonlinearity
mhw1 = nlhw(z, NN, InputNL, OutputNL, opt)

mhw1 =

<strong>Hammerstein-Wiener model with 2 outputs and 6 inputs</strong>

Linear transfer function matrix corresponding to the orders:
   nb = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nf = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nk = [1 1 1 1 1 1; 1 1 1 1 1 1]

Input nonlinearities:
  Input 1:  Piecewise linear with 10 break-points
  Input 2:  Piecewise linear with 10 break-points
  Input 3:  Piecewise linear with 10 break-points
  Input 4:  Piecewise linear with 10 break-points
  Input 5:  Piecewise linear with 10 break-points
  Input 6:  Piecewise linear with 10 break-points

Output nonlinearities:
  Output 1:  None
  Output 2:  None

Sample time: 0.02 seconds

Status:                                                     
Estimated using NLHW on time domain data "Motorized Camera".
Fit to estimation data: [98.46;97.93]%                      
FPE: 7.931, MSE: 2.217

Examine the result with compare.

compare(z, 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).

When the linear block view is chosen, by default all the 12 channels are shown in very reduced sizes. Choose one of the channels to have a better view. It is possible to choose the type of plots within Step response, Bode plot, Impulse response and Pole-zero map.

plot(mhw1)

The result shown by compare was quite good, so let us keep the same model orders in the following trials.

nbnfnk = [ones(2,6), ones(2,6), ones(2,6)];

Estimation of Wiener Model (No Input Nonlinearity)

Let us try a Wiener model. Notice that the absence of input nonlinearity can be indicated by "[]" instead of the idUnitGain object.

opt.SearchOptions.MaxIterations = 10;
mhw2 = nlhw(z, nbnfnk, [], idPiecewiseLinear, opt)

mhw2 =

<strong>Hammerstein-Wiener model with 2 outputs and 6 inputs</strong>

Linear transfer function matrix corresponding to the orders:
   nb = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nf = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nk = [1 1 1 1 1 1; 1 1 1 1 1 1]

Input nonlinearities:
  Input 1:  None
  Input 2:  None
  Input 3:  None
  Input 4:  None
  Input 5:  None
  Input 6:  None

Output nonlinearities:
  Output 1:  Piecewise linear with 10 break-points
  Output 2:  Piecewise linear with 10 break-points

Sample time: 0.02 seconds

Status:                                                     
Estimated using NLHW on time domain data "Motorized Camera".
Fit to estimation data: [73.85;71.36]%                      
FPE: 1.314e+05, MSE: 503.8

Estimation of Hammerstein-Wiener Model (Both Input and Output Nonlinearities)

Indicate both input and output nonlinearities for a Hammerstein-Wiener model.

mhw3 = nlhw(z, nbnfnk,idSaturation, idPiecewiseLinear, opt)

mhw3 =

<strong>Hammerstein-Wiener model with 2 outputs and 6 inputs</strong>

Linear transfer function matrix corresponding to the orders:
   nb = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nf = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nk = [1 1 1 1 1 1; 1 1 1 1 1 1]

Input nonlinearities:
  Input 1:  Saturation with linear interval: [0.0103 0.0562]
  Input 2:  Saturation with linear interval: [-0.00143 0.000909]
  Input 3:  Saturation with linear interval: [-0.0947 -0.0185]
  Input 4:  Saturation with linear interval: [-0.00385 0.00527]
  Input 5:  Saturation with linear interval: [0.0195 0.13]
  Input 6:  Saturation with linear interval: [-0.00302 0.000387]

Output nonlinearities:
  Output 1:  Piecewise linear with 10 break-points
  Output 2:  Piecewise linear with 10 break-points

Sample time: 0.02 seconds

Status:                                                     
Estimated using NLHW on time domain data "Motorized Camera".
Fit to estimation data: [86.88;84.55]%                      
FPE: 1.111e+04, MSE: 137.3

The limit values of the idSaturation function can be accessed as follows:

mhw3.InputNonlinearity(1).LinearInterval % view Linear Interval of saturation function

ans =

    0.0103    0.0562

Similarly, the break points of the idPiecewiseLinear function can be accessed as follows:

mhw3.OutputNonlinearity(1).BreakPoints

ans =

  Columns 1 through 7

   17.1233   34.2491   51.3726   68.4968   85.6230  102.7478  119.8742
    2.6184   16.0645   45.5178   41.9621   62.3246   84.9038  112.2970

  Columns 8 through 10

  136.9991  154.1238  171.2472
  135.4543  156.1016  173.2701

Hammerstein-Wiener Model - Using Mixed Nonlinear Functions

Different nonlinear functions can be mixed in a same model. Suppose we want a model with: - No nonlinearity on any output channels ("Hammerstein Model") - Piecewise linear nonlinearity on input channel #1 with 3 units - Saturation nonlinearity on input channel #2 - Dead Zone nonlinearity on input channel #3 - Sigmoid Network nonlinearity on input channel #4 - No nonlinearity (specified by a unitgain object) on input channel #5 - Sigmoid Network nonlinearity on input channel #6 with 5 units

We can create an array of nonlinear mapping function objects of chosen types and pass it to the estimation function NLHW as input nonlinearity.

inputNL = [idPiecewiseLinear; ...
    idSaturation; ...
    idDeadZone; ...
    idSigmoidNetwork; ...
    idUnitGain; ...
    idSigmoidNetwork(5)];
inputNL(1).NumberOfUnits = 3;
opt.SearchOptions.MaxIterations = 25;
mhw4 = nlhw(z, nbnfnk, inputNL, [], opt); % "[]" as the 4th input argument  denotes no output nonlinearity

Hammerstein-Wiener Model - Specifying Initial Guess for SATURATION and DEADZONE

The initial guess for the linear interval of saturation and the zero interval of dead zone can be directly indicated when these objects are created; you can also specify constraints on these values such as whether they are fixed to their specified values (by setting Free attribute to false), or if their estimations are subject to minimum/maximum bounds (using the Minimum and Maximum attributes).

Suppose we want to set the saturation's linear interval to [10 200] and the deadzone's zero interval to [11 12] as initial guesses. Furthermore, we want the upper limit of the saturation to remain fixed. We can achieve this configuration as follows.

% Create  nonlinear functions with initial guesses for properties.
OutputNL1 = idSaturation([10 200]);
OutputNL1.Free(2) = false; % the upper limit is a fixed value
OutputNL2 = idDeadZone([11 12]);

mhw5 = idnlhw(nbnfnk, [], [OutputNL1; OutputNL2], 'Ts',z.Ts);

Notice the use of the IDNLHW model object constructor idnlhw, not the estimator nlhw. The resulting model object mhw5 is not yet estimated from data. The limit values of the saturation and deadzone functions can be accessed. They still have their initial values, because they are not yet estimated from data.

mhw5.OutputNonlinearity(1).LinearInterval % view the linear interval on saturation at first output channel after estimation
mhw5.OutputNonlinearity(2).ZeroInterval   % view the zero interval on dead zone at second output channel after estimation

ans =

    10   200


ans =

    11    12

What are the values of these limits after estimation?

opt.SearchOptions.MaxIterations = 15;
mhw5 = nlhw(z, mhw5, opt)  % estimate the model from data
mhw5.OutputNonlinearity(1).LinearInterval % show linear interval on saturation at first output channel after estimation
mhw5.OutputNonlinearity(2).ZeroInterval    % show zero interval on dead zone at second output channel after estimation

mhw5 =

<strong>Hammerstein-Wiener model with 2 outputs and 6 inputs</strong>

Linear transfer function matrix corresponding to the orders:
   nb = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nf = [1 1 1 1 1 1; 1 1 1 1 1 1]
   nk = [1 1 1 1 1 1; 1 1 1 1 1 1]

Input nonlinearities:
  Input 1:  None
  Input 2:  None
  Input 3:  None
  Input 4:  None
  Input 5:  None
  Input 6:  None

Output nonlinearities:
  Output 1:  Saturation with linear interval: [10 200]
  Output 2:  Dead Zone with zero interval: [11 12]

Sample time: 0.02 seconds

Status:                                                     
Estimated using NLHW on time domain data "Motorized Camera".
Fit to estimation data: [27.12;6.857]%                      
FPE: 3.373e+06, MSE: 4666

ans =

    9.9974  200.0000


ans =

   11.0020   12.0011

Post Estimation Analysis - Comparing Different Models

Models of different nature (IDNLARX and IDNLHW) can be compared in the same compare command.

compare(z,mw2,mhw1)

warning(ws) % reset the warning state