This example shows how to estimate multi-input multi-output (MIMO) nonlinear black box models from data. Two types of nonlinear black box models are offered in the toolbox - Nonlinear ARX and Hammerstein-Wiener models.
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');
Let us first try nonlinear ARX models. Two important elements need to be chosen: the model orders and the nonlinearity estimators. To try simple things first, let us choose the orders [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. By using vector notations, the model can be written in the form of model_output(t) = F(y(t-1), u(t-1)). Let us choose Wavelet Network (wavenet) as nonlinearity estimator, thus the vectorial nonlinear function F will be estimated by two wavelet networks. The function nlarx facilitates estimation of Nonlinear ARX models.
mw1 = nlarx(z, [ones(2,2), ones(2,6), ones(2,6)], wavenet);
Examine the result by comparing the output simulated with the estimated model and the output in the measured data z:
Let us see 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 reliably calculated. 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)], wavenet); compare(z,mw2)
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)];
The numbers of units (wavelets) of the two WAVENET estimators have been automatically chosen by the estimation algorithm. These numbers are displayed below. Notice the abbreviations 'nl'='Nonlinearity' and 'num'='NumberOfUnits'.
mw2.Nonlinearity(1).NumberOfUnits %using full property names mw2.nl(2).num %using short-hand notations
ans = 27 ans = 25
The number of units in the WAVENET estimators can be explicitly specified instead of being automatically chosen by the estimation algorithm:
mw3 = nlarx(z, nanbnk, [wavenet('num',10); wavenet('num',5)]);
At the place of the WAVENET estimator, other nonlinearity estimators can also be used. Let us try the TREEPARTITION estimator.
mt1 = nlarx(z, nanbnk, 'treepartition');
In the above call, we have used a character vector ('treepartition') in place of an object to specify the nonlinearity estimator. However, this works only if the nonlinearity is to be used with default property values. In the following example, the treepartition object constructor is called to directly create a nonlinearity estimator object.
mt2 = nlarx(z, nanbnk, treepartition);
The SIGMOIDNET estimator can also be used. Estimation options such as maximum iterations (MaxIter) and iteration display can be specified using NLARXOPTIONS command.
opt = nlarxOptions('Display','on'); opt.SearchOption.MaxIter = 2; ms1 = nlarx(z, nanbnk, 'sigmoidnet', opt);
It is possible to use different nonlinearity estimators on different output channels in the same model. Suppose we want to use a tree partition nonlinearity estimator for prediction of first output and use a wavelet network for prediction of the second output. The model estimation is shown below. The third input argument (nonlinearity) is now an array of two different nonlinearity estimator objects.
mtw = nlarx(z, nanbnk, [treepartition; wavenet]);
The absence of nonlinearity in an output channel can be indicated by choosing a LINEAR estimator. 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 and a nonlinear component (SIGMOIDNET).
opt.Display = 'off'; % do not show estimation progress anymore mls = nlarx(z, nanbnk, [linear; sigmoidnet('NumberOfUnits',5)], opt); % no nonlinearity on first output
Various models can be compared in the same COMPARE command.
compare(z, mw2, ms1, mls)
Function PLOT may be used to view the nonlinearity responses of various models.
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.
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 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 (UNITGAIN) nonlinearities.
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 piecewise linear estimator. Notice that the entered single PWLINEAR object is automatically expanded to all the 6 input channels. UNITGAIN indicates absence of nonlinearity.
opt = nlhwOptions(); opt.SearchOption.MaxIter = 2; mhw1 = nlhw(z,[ones(2,6), ones(2,6), ones(2,6)], pwlinear, unitgain, opt);
Examine the result with COMPARE.
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.
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)];
Let us try a Wiener model. Notice that the absence of input nonlinearity can be indicated by "" instead of the UNITGAIN object.
opt.SearchOption.MaxIter = 4; mhw2 = nlhw(z, nbnfnk, , 'pwlinear', opt);
Indicate both input and output nonlinearities for a Hammerstein-Wiener model. As in the case of Nonlinear ARX models, we can use a character vector (rather than object) to specify the nonlinearity estimator.
mhw3 = nlhw(z, nbnfnk,'saturation', 'pwlinear', opt);
The limit values of the SATURATION estimators can be accessed. The short-hands 'u'='input', 'y'='output', and 'nl'='Nonlinearity' are available.
mhw3.InputNonlinearity(1).LinearInterval % view Linear Interval of SATURATION mhw3.unl(3).LinearInterval % access using using short-hand notation
ans = 0.0103 0.0562 ans = -0.0956 -0.0169
The break points of the PWLINEAR estimator can also be accessed.
ans = Columns 1 through 7 17.1246 34.2518 51.3765 68.5021 85.6296 102.7558 119.8835 2.6730 16.1133 45.5627 42.0046 62.3661 84.9454 112.3400 Columns 8 through 10 137.0097 154.1357 171.2604 135.5004 156.1533 173.3304
Different nonlinearity estimators 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 nonlinearity object estimators of chosen types and pass it to the estimation function NLHW as input nonlinearity.
inputNL = [pwlinear('num',3); saturation; deadzone; sigmoidnet; unitgain; sigmoidnet('num',5)]; opt.SearchOption.MaxIter = 2; mhw4 = nlhw(z, nbnfnk, inputNL, , opt); % "" for 4th input denotes no output nonlinearity
The initial guess for the linear interval of SATURATION and the zero interval of DEADZONE can be directly indicated when the estimator objects are created. 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.
mhw5 = idnlhw(nbnfnk, , [saturation([10 200]); deadzone([11 12])],'Ts',z.Ts); % Create nonlinearity estimator with initial guesses for properties. % 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 estimators can be accessed. They still have their initial values, because they are not yet estimated from data.
mhw5.outputNL(1).LinearInterval % show linear interval on saturation at first output channel after estimation mhw5.yNL(2).ZeroInterval % show 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.SearchOption.MaxIter = 5; mhw5 = nlhw(z, mhw5, opt); % estimate the model from data mhw5.outputNL(1).LinearInterval % show linear interval on saturation at first output channel after estimation mhw5.yNL(2).ZeroInterval % show zero interval on dead zone at second output channel after estimation
ans = 9.9997 199.9989 ans = 11.0001 12.0015
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
For more information on identification of dynamic systems with System Identification Toolbox visit the System Identification Toolbox product information page.