You can estimate Hammerstein-Wiener models after performing the following tasks:
Prepare your data, as described in Preparing Data for Nonlinear Identification.
(Optional) Choose a nonlinearity estimator in Available Nonlinearity Estimators for Hammerstein-Wiener Models.
(Optional) Estimate or construct an input-output polynomial
model of Output-Error (OE) structure (
or a state-space model with no disturbance component (
idss with K=0)
for initialization of Hammerstein-Wiener model. See Initialize Hammerstein-Wiener Estimation Using Linear Model.
Start with the simplest estimation using
m = nlhw(data,[nb nf nk]). For example:
load iddata3; % nb = nf = 2 and nk = 1 m = nlhw(z3,[2 2 1])
m = Hammerstein-Wiener model with 1 output and 1 input Linear transfer function corresponding to the orders nb = 2, nf = 2, nk = 1 Input nonlinearity: Piecewise linear with 10 break-points Output nonlinearity: Piecewise linear with 10 break-points Sample time: 1 seconds Status: Estimated using NLHW on time domain data "z3". Fit to estimation data: 75.31% FPE: 2.019, MSE: 1.472
The second input argument [
nk] sets the order of the linear transfer function, where
nb is the number of zeros plus 1,
nf is the number of poles, and
nk is the input delay. By default, both the input and
output nonlinearity estimators are piecewise linear functions (see the
idPiecewiseLinear reference page).
m is an
For MIMO systems,
nk are ny-by-nu matrices.
nlhw reference page
for more information about MIMO estimation.
You can specify a different nonlinearity estimator than the default piecewise linear estimators.
m = nlhw(data,[nb,nf,nk],InputNL,OutputNL)
OutputNL are nonlinearity estimator objects.
If your input signal is binary, set
To use nonlinearity estimators with default settings, specify
OutputNL using constructors with no input arguments or their names as character vectors (such as 'idWaveletNetwork' for wavelet network or 'idSigmoidNetwork' for sigmoid network).
load iddata3; m = nlhw(z3,[2 2 1],'idSigmoidNetwork','idDeadZone');
If you need to configure the properties of a nonlinearity estimator, use its object representation. For example, to estimate a Hammerstein-Wiener model that uses saturation as its input nonlinearity and one-dimensional polynomial of degree 3 as its output nonlinearity:
m = nlhw(z3,[2 2 1],'idSaturation',idPolynomial1D(3));
The third input
'idSaturation' specifies the saturation nonlinearity with default property values.
idPolynomial1D(3) creates a one-dimensional polynomial object of degree 3. Of course, you could have used the constructor idSaturation directly in place of the character vector 'idSaturation'.
For MIMO models, specify the nonlinearities using objects, or cell array of character vectors representing the nonlinearity names. If single name (character vector) used, the same value is applied to all the channels.
This table summarizes values that specify the nonlinearity estimators.
|Nonlinearity||Value (Default Nonlinearity Configuration)||Class|
|One layer sigmoid network|
Additional available nonlinearities include custom networks that you create. Specify a custom
network by defining a function called
gaussunit.m, as described
idCustomNetwork reference page.
Define the custom network object
CNetw = idCustomNetwork(@gaussunit); m = nlhw(z3,[2 2 1],'idSaturation',CNetw);
For more information, see Available Nonlinearity Estimators for Hammerstein-Wiener Models.
Exclude a nonlinearity for a specific channel by specifying the
idUnitGain value for the
If the input signal is binary, set
For a description of each nonlinearity estimator, see Available Nonlinearity Estimators for Hammerstein-Wiener Models.
Estimate a Hammerstein-Wiener model and then use
nlhw command to iteratively refine the model.
load iddata3; m1 = nlhw(z3,[2 2 1],idSigmoidNetwork, idWaveletNetwork); m2 = nlhw(z3,m1);
pem to refine the model.
m2 = pem(z3,m1);
Check the search termination criterion in
WhyStop indicates that the estimation reached the maximum number of iterations, try repeating the estimation and possibly specifying a larger value for the
Run 30 more iterations starting at model
opt = nlhwOptions; opt.SearchOptions.MaxIterations = 30; m2 = nlhw(z3,m1,opt);
m.Report.Termination.WhyStop value is
Near (local) minimum, (norm(g) < tol or
No improvement along the search direction with line search, validate your model to see if this model adequately fits the data. If not, the solution might be stuck in a local minimum of the cost-function surface. Try adjusting the
SearchOptions.Tolerance or the
SearchMethod option of the
nlhw option set, and repeat the estimation.
You can also try perturbing the parameters of the last model using
init, and then refine the model using
Randomly perturb parameters of original model
m1 about nominal values.
m1p = init(m1);
Estimate the parameters of perturbed model.
M2 = nlhw(z3,m1p);
Note that using
init does not guarantee a better solution on further refinement.
If your estimated Hammerstein-Wiener model provides a poor fit to measured data, you can repeat the estimation using the initial state values estimated from the data. By default, the initial states corresponding to the linear block of the Hammerstein-Wiener model are zero.
To specify estimating initial states during model estimation:
load iddata3; opt = nlhwOptions('InitialCondition', 'estimate'); m = nlhw(z3,[2 2 1],idSigmoidNetwork,,opt);
If you do not get a satisfactory model after many trials with various model structures and estimation options, it is possible that the data is poor. For example, your data might be missing important input or output variables and does not sufficiently cover all the operating points of the system.
Nonlinear black-box system identification usually requires more data than linear model identification to gain enough information about the system. See also Troubleshooting Model Estimation.
This example shows how to estimate and compare multiple Hammerstein-Wiener models using measured input-output data.
Load estimation and validation data.
load twotankdata z = iddata(y,u,0.2); ze = z(1:1000); zv = z(1001:3000);
Estimate several models using the estimation data
ze and different model orders, delays, and nonlinearity settings.
m1 = nlhw(ze,[2 3 1]); m2 = nlhw(ze,[2 2 3]); m3 = nlhw(ze,[2 2 3],idPiecewiseLinear('NumberofUnits',13),idPiecewiseLinear('NumberofUnits',10)); m4 = nlhw(ze,[2 2 3],idSigmoidNetwork('NumberofUnits',2),idPiecewiseLinear('NumberofUnits',10));
An alternative way to perform the estimation is to configure the model structure first using
idnlhw, and then estimate the model.
m5 = idnlhw([2 2 3],idDeadZone,idSaturation); m5 = nlhw(ze,m5);
Compare the resulting models by plotting the model outputs and the measured output in validation data
This example shows how to use the Hammerstein-Wiener model structure to improve a previously estimated linear model.
After estimating the linear model, insert it into the Hammerstein-Wiener structure that includes input or output nonlinearities.
Estimate a linear model.
load iddata1 LM = arx(z1,[2 2 1]);
Extract the transfer function coefficients from the linear model.
[Num,Den] = tfdata(LM);
Create a Hammerstein-Wiener model, where you initialize the linear block properties
nb = 1; % In general, nb = ones(ny,nu) % ny is number of outputs and nu is number of inputs nf = nb; nk = 0; % In general, nk = zeros(ny,nu) % ny is number of outputs and nu is number of inputs M = idnlhw([nb nf nk],,'idPiecewiseLinear'); M.B = Num; M.F = Den;
Estimate the model coefficients, which refines the linear model coefficients in
M = nlhw(z1,M);
Compare responses of linear and nonlinear model against measured data.