nlhw

Estimate Hammerstein-Wiener model

Description

example

sys = nlhw(Data,Orders) creates and estimates a Hammerstein-Wiener model using the estimation data, model orders and delays, and default piecewise linear functions as input and output nonlinearity estimators.

example

sys = nlhw(Data,Orders,InputNL,OutputNL) specifies InputNL and OutputNL as the input and output nonlinearity estimators, respectively.

example

sys = nlhw(Data,LinModel) uses a linear model to specify the model orders and delays, and default piecewise linear functions for the input and output nonlinearity estimators.

example

sys = nlhw(Data,LinModel,InputNL,OutputNL) specifies InputNL and OutputNL as the input and output nonlinearity estimators, respectively.

example

sys = nlhw(Data,sys0) refines or estimates the parameters of a Hammerstein-Wiener model, sys0, using the estimation data.

Use this syntax to:

• Update the parameters of a previously estimated model to improve the fit to the estimation data. In this case, the estimation algorithm uses the parameters of sys0 as initial guesses.

• Estimate the parameters of a model previously created using the idnlhw constructor. Prior to estimation, you can configure the model properties using dot notation.

example

sys = nlhw(___,Options) specifies additional model estimation options. Use Options with any of the previous syntaxes.

Examples

collapse all

m1 = nlhw(z3,[4 2 1]);

z = iddata(y,u,0.2,'Name','Two tank system');
z1 = z(1:1000);

Create a saturation object with lower limit of 0 and upper limit of 5.

InputNL = idSaturation('LinearInterval',[0 5]);

Estimate model with no output nonlinearity.

m = nlhw(z1,[2 3 0],InputNL,[]);

Generating a custom network nonlinearity requires the definition of a user-defined unit function.

Define the unit function and save it as gaussunit.m.

function [f,g,a] = gaussunit(x)
% Custom unit function nonlinearity.
%
% Copyright 2015 The MathWorks, Inc.
f = exp(-x.*x);
if nargout>1
g = -2*x.*f;
a = 0.2;
end

Create a custom network nonlinearity using the gaussunit function.

H = @gaussunit;
CNet = idCustomNetwork(H);

z = iddata(y,u,0.2,'Name','Two tank system');
z1 = z(1:1000);

Estimate a Hammerstein-Wiener model using the custom network.

m = nlhw(z1,[5 1 3],CNet,[]);

Estimate linear OE model.

Tr = getTrend(ThrottleData);
Tr.OutputOffset = 15;
DetrendedData = detrend(ThrottleData, Tr);
opt = oeOptions('Focus','simulation');
LinearModel = oe(DetrendedData,[1 2 1],opt);

Estimate Hammerstein-Wiener model using OE model as its linear component and saturation as its output nonlinearity.

sys = nlhw(ThrottleData,LinearModel,[],idSaturation);

Construct a Hammerstein-Wiener model using idnlhw to define the model properties B and F.

sys0 = idnlhw([2,2,0],[],'idWaveletNetwork');
sys0.B{1} = [0.8,1];
sys0.F{1} = [1,-1.2,0.5];

Estimate the model.

sys = nlhw(z1,sys0);

Estimate a Hammerstein-Wiener model using nlhw to define the model properties B and F.

sys2 = nlhw(z1,[2,2,0],[],'idWaveletNetwork','B',{[0.8,1]},'F',{[1,-1.2,0.5]});

Compare the two estimated models to see that they are equivalent.

compare(z1,sys,'g',sys2,'r--'); Estimate a Hammerstein-Wiener Model.

sys = nlhw(z3,[4 2 1],'idSigmoidNetwork','idWaveletNetwork');

Refine the model, sys.

sys = nlhw(z3,sys);

Create estimation option set for nlhw to view estimation progress, use the Levenberg-Marquardt search method, and set the maximum iteration steps to 50.

opt = nlhwOptions;
opt.Display = 'on';
opt.SearchMethod = 'lm';
opt.SearchOptions.MaxIterations = 50;

Load data and estimate the model.

sys = nlhw(z3,[4 2 1],idSigmoidNetwork,idPiecewiseLinear,opt);

Input Arguments

collapse all

Time-domain estimation data, specified as an iddata

Order and delays of the linear subsystem transfer function, specified as a [nb nf nk] vector.

Dimensions of Orders:

• For a SISO transfer function, Orders is a vector of positive integers.

nb is the number of zeros plus 1, nf is the number of poles, and nk is the input delay.

• For a MIMO transfer function with nu inputs and ny outputs, Orders is a vector of matrices.

nb, nf, and nk are ny-by-nu matrices whose i-jth entry specifies the orders and delay of the transfer function from the jth input to the ith output.

Input static nonlinearity estimator, specified as one of the following.

 'idPiecewiseLinear' or idPiecewiseLinear object (default) Piecewise linear function 'idSigmoidNetwork' or idSigmoidNetwork object Sigmoid network 'idWaveletNetwork' or idWaveletNetwork object Wavelet network 'idSaturation' or idSaturation object Saturation 'idDeadZone' or idDeadZone object Dead zone 'idPolynomial1D' or idPolynomial1D object One-dimensional polynomial 'idUnitGain' or [] oridUnitGain object Unit gain idCustomNetwork object Custom network — Similar to idSigmoidNetwork, but with a user-defined replacement for the sigmoid function.

Specifying a character vector, for example 'idSigmoidNetwork', creates a nonlinearity estimator object with default settings. Use object representation instead to configure the properties of a nonlinearity estimator.

InputNL = idWaveletNetwork;
InputNL.NumberOfUnits = 10;

Alternatively, use the associated input nonlinearity estimator function with Name-Value pair arguments.

InputNL = idWaveletNetwork('NumberOfUnits',10);

For nu input channels, you can specify nonlinear estimators individually for each input channel by setting InputNL to an nu-by-1 array of nonlinearity estimators.

To specify the same nonlinearity for all inputs, specify a single input nonlinearity estimator.

Output static nonlinearity estimator, specified as one of the following:

 'idPiecewiseLinear' or idPiecewiseLinear object (default) Piecewise linear function 'idSigmoidNetwork' or idSigmoidNetwork object Sigmoid network 'idWaveletNetwork' or idWaveletNetwork object Wavelet network 'idSaturation' or idSaturation object Saturation 'idDeadZone' or idDeadZone object Dead zone 'idPolynomial1D' or idPolynomial1D object One-dimensional polynomial 'idUnitGain' or [] oridUnitGain object Unit gain idCustomNetwork object Custom network — Similar to idSigmoidNetwork, but with a user-defined replacement for the sigmoid function.

Specifying a character vector creates a nonlinearity estimator object with default settings. Use object representation instead to configure the properties of a nonlinearity estimator.

OutputNL = idSigmoidNetwork;
OutputNL.NumberOfUnits = 10;

Alternatively, use the associated input nonlinearity estimator function with Name-Value pair arguments.

OutputNL = idSigmoidNetwork('NumberOfUnits',10);

For ny output channels, you can specify nonlinear estimators individually for each output channel by setting OutputNL to an ny-by-1 array of nonlinearity estimators. To specify the same nonlinearity for all outputs, specify a single output nonlinearity estimator.

Discrete-time linear model used to specify the linear subsystem, specified as one of the following:

• Input-output polynomial model of Output-Error (OE) structure (idpoly)

• State-space model with no disturbance component (idss with K = 0)

• Transfer function model (idtf)

Typically, you estimate the model using oe, n4sid, or tfest.

Hammerstein-Wiener model, specified as an idnlhw object. sys0 can be:

• A model previously created using idnlhw to specify model properties.

• A model previously estimated using nlhw, that you want to update using a new estimation data set.

You can also refine sys0 using the original estimation data set. If the previous estimation stopped when the numerical search was stuck at a local minima of the cost function, use init to first randomize the parameters of sys0. See sys0.Report.Termination for search stopping conditions. Using init does not guarantee a better solution on further refinement.

Estimation options for Hammerstein-Wiener model identification, specified as an nlhwOptions option set.

Output Arguments

collapse all

Estimated Hammerstein-Wiener model, returned as an idnlhw object. The model is estimated using the specified model orders and delays, input and output nonlinearity estimators, and estimation options.

Information about the estimation results and options used is stored in the Report property of the model. Report has the following fields:

Report FieldDescription
Status

Summary of the model status, which indicates whether the model was created by construction or obtained by estimation.

Method

Estimation command used.

Fit

Quantitative assessment of the estimation, returned as a structure. See Loss Function and Model Quality Metrics for more information on these quality metrics. The structure has the following fields:

FieldDescription
FitPercent

Normalized root mean squared error (NRMSE) measure of how well the response of the model fits the estimation data, expressed as the percentage fitpercent = 100(1-NRMSE).

LossFcn

Value of the loss function when the estimation completes.

MSE

Mean squared error (MSE) measure of how well the response of the model fits the estimation data.

FPE

Final prediction error for the model.

AIC

Raw Akaike Information Criteria (AIC) measure of model quality.

AICc

Small-sample-size corrected AIC.

nAIC

Normalized AIC.

BIC

Bayesian Information Criteria (BIC).

Parameters

Estimated values of model parameters.

OptionsUsed

Option set used for estimation. If no custom options were configured, this is a set of default options. See nlhwOptions for more information.

RandState

State of the random number stream at the start of estimation. Empty, [], if randomization was not used during estimation. For more information, see rng.

DataUsed

Attributes of the data used for estimation, returned as a structure with the following fields.

FieldDescription
Name

Name of the data set.

Type

Data type.

Length

Number of data samples.

Ts

Sample time.

InterSample

Input intersample behavior, returned as one of the following values:

• 'zoh' — Zero-order hold maintains a piecewise-constant input signal between samples.

• 'foh' — First-order hold maintains a piecewise-linear input signal between samples.

• 'bl' — Band-limited behavior specifies that the continuous-time input signal has zero power above the Nyquist frequency.

InputOffset

Offset removed from time-domain input data during estimation. For nonlinear models, it is [].

OutputOffset

Offset removed from time-domain output data during estimation. For nonlinear models, it is [].

Termination

Termination conditions for the iterative search used for prediction error minimization, returned as a structure with the following fields:

FieldDescription
WhyStop

Reason for terminating the numerical search.

Iterations

Number of search iterations performed by the estimation algorithm.

FirstOrderOptimality

$\infty$-norm of the gradient search vector when the search algorithm terminates.

FcnCount

Number of times the objective function was called.

UpdateNorm

Norm of the gradient search vector in the last iteration. Omitted when the search method is 'lsqnonlin' or 'fmincon'.

LastImprovement

Criterion improvement in the last iteration, expressed as a percentage. Omitted when the search method is 'lsqnonlin' or 'fmincon'.

Algorithm

Algorithm used by 'lsqnonlin' or 'fmincon' search method. Omitted when other search methods are used.

For estimation methods that do not require numerical search optimization, the Termination field is omitted.

Compatibility Considerations

expand all

Not recommended starting in R2021b

Extended Capabilities

Introduced in R2007a

System Identification Toolbox Documentation Get trial now