nlarx

Estimate parameters of nonlinear ARX model

Syntax

sys = nlarx(Data,Orders)

sys = nlarx(Data,Orders,Nonlinearity)

sys = nlarx(Data,Orders,Nonlinearity,Name,Value)

sys = nlarx(Data,LinModel)

sys = nlarx(Data,LinModel,Nonlinearity)

sys = nlarx(Data,LinModel,Nonlinearity,Name,Value)

sys = nlarx(Data,sys0)

sys = nlarx(___,Options)

Description

sys = nlarx(Data,Orders) estimates a nonlinear ARX model to fit the given estimation data using the specified orders and a default wavelet network nonlinearity estimator.

sys = nlarx(Data,Orders,Nonlinearity) specifies the nonlinearity to use for model estimation.

sys = nlarx(Data,Orders,Nonlinearity,Name,Value) specifies additional attributes of the estimated model using one or more Name,Value pair arguments. These attributes include the nonlinear and custom regressor structure, and the data properties of the idnlarx model.

sys = nlarx(Data,LinModel) uses a linear ARX model, LinModel, to specify the model orders and the initial values of the linear coefficients of the model.

sys = nlarx(Data,LinModel,Nonlinearity) specifies the nonlinearity to use for model estimation.

sys = nlarx(Data,LinModel,Nonlinearity,Name,Value) specifies additional attributes of the idnlarx model structure using one or more Name,Value pair arguments.

sys = nlarx(Data,sys0) refines the parameters of the nonlinear ARX model, sys0.

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 idnlarx constructor. Prior to estimation, you can configure the model properties using dot notation.

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

Examples

Estimate Nonlinear ARX Model with Default Settings

Load the estimation data.

load twotankdata;

Create an iddata object from the estimation data with a sample time of 0.2 min.

Ts = 0.2;
z = iddata(y,u,Ts);

Estimate the nonlinear ARX model.

sys = nlarx(z,[4 4 1]);

Estimate Nonlinear ARX Model from Time Series Data

Create time and data arrays.

dt = 0.01;
t = 0:dt:10;
y = 10*sin(2*pi*t)+rand(size(t));

Create an iddata object with no input signal specified.

z = iddata(y',[],dt);

Estimate the nonlinear ARX model.

sys = nlarx(z,2);

Estimate Nonlinear ARX Model with Specific Nonlinearity

Load the estimation data.

load twotankdata;

Create an iddata object from the estimation data.

z = iddata(y,u,0.2);

Create a wavelet network nonlinearity estimator with 5 units.

NL = wavenet('NumberOfUnits',5);

Estimate the nonlinear ARX model.

sys = nlarx(z,[4 4 1],NL);

Estimate Nonlinear ARX Model Using Custom Network Nonlinearity

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

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

% Copyright 2015 The MathWorks, Inc.

function [f, g, a] = gaussunit(x)
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 = customnet(H);

Load the estimation data.

load iddata1;

Estimate a nonlinear ARX model using the custom network.

sys = nlarx(z1,[1 2 1],CNet);

Estimate MIMO Nonlinear ARX Model

Load the estimation data.

load motorizedcamera;

Create an iddata object.

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

z is an iddata object with 6 inputs and 2 outputs.

Specify the model orders.

Orders = [ones(2,2),2*ones(2,6),ones(2,6)];

Specify different nonlinearity estimators for each output channel.

NL = [wavenet('NumberOfUnits',2),linear];

Estimate the nonlinear ARX model.

sys = nlarx(z,Orders,NL);

Estimate MIMO Nonlinear ARX Model with Same Nonlinearity for All Outputs

Load the estimation data and create an iddata object.

load motorizedcamera;
z = iddata(y,u,0.02);

Specify the model orders.

Orders = [ones(2,2),2*ones(2,6),ones(2,6)];

Estimate a nonlinear ARX model using a sigmoidnet nonlinearity with 4 units for all output channels.

m = nlarx(z,Orders,sigmoidnet('numberOfUnits',4));

Estimate Nonlinear ARX Model with Custom Regressors

Load the estimation data.

load iddata1;

Create a cell array with two custom regressors.

C = {'y1(t-1)^2','y1(t-2)*u1(t-3)'};

Estimate a nonlinear ARX model with custom regressors and no standard regressors.

sys = nlarx(z1,[0 0 0],'linear','CustomRegressors',C);

Estimate Nonlinear ARX Model with Custom Regressor Objects

Load the estimation data.

load iddata1;

Define a custom regressor object for y1(t-1)^2.

C1 = customreg(@(x)x^2,{'y1'},[1]);

Define a custom regressor object for y1(t-2)*u1(t-3).

C2 = customreg(@(x,y)x*y,{'y1','u1'},[2 3]);

Create a custom regressor object array.

C = [C1,C2];

Estimate a nonlinear ARX model with custom regressors.

sys = nlarx(z1,[0 0 0],'linear','CustomRegressors',C);

List the model regressors.

getreg(sys);

Regressors:
    y1(t-1)^2
    y1(t-2)*u1(t-3)

Estimate Nonlinear ARX Model Searching for Optimum Nonlinear Regressors

Load the estimation data.

load iddata1;

Estimate a Nonlinear ARX model using the 'search' option.

sys = nlarx(z1,[4 4 1],'sigmoidnet','NonlinearRegressors','search');

List the model nonlinear regressor indices.

sys.NonlinearRegressors

ans = 1×4

     3     5     6     7

List all of the model regressors.

getreg(sys)

Regressors:
    y1(t-1)
    y1(t-2)
    y1(t-3)
    y1(t-4)
    u1(t-1)
    u1(t-2)
    u1(t-3)
    u1(t-4)

The optimum set of nonlinear regressors for this model includes y1(t-3), u1(t-1), u1(t-2), and u1(t-3).

Estimate Nonlinear ARX Model with No Linear Term in Nonlinearity Estimator

Load the estimation data.

load iddata1;

Create a sigmoid network nonlinearity estimator with no linear term.

SNL = sigmoidnet('LinearTerm','off');

Estimate the nonlinear ARX model.

sys = nlarx(z1,[2 2 1],SNL);

Specify Nonlinear ARX Orders and Linear Parameters Using Linear ARX Model

Load the estimation data.

load throttledata;

Detrend the data.

Tr = getTrend(ThrottleData);
Tr.OutputOffset = 15;
DetrendedData = detrend(ThrottleData,Tr);

Estimate the linear ARX model.

LinearModel = arx(DetrendedData,[2 1 1]);

Estimate the nonlinear ARX model using the linear model. The model orders, delays, and linear parameters of NonlinearModel are derived from LinearModel.

NonlinearModel = nlarx(ThrottleData,LinearModel);

Estimate Nonlinear ARX Model Using Constructed idnlarx Object

Load the estimation data.

load iddata1;

Create an idnlarx model.

sys = idnlarx([2 2 1]);

Configure the model using dot notation to set the following parameters:

• Use a sigmoid network nonlinearity

• Search for an optimum nonlinear regressor subset

sys.Nonlinearity = 'sigmoidnet';
sys.NonlinearRegressors = 'search';

Estimate a nonlinear ARX model with the structure and properties specified in the idnlarx object.

sys = nlarx(z1,sys);

Estimate Nonlinear ARX Model and Avoid Local Minima

If an estimation stops at a local minimum, you can perturb the model using init and reestimate the model.

Load the estimation data.

load iddata1;

Estimate the initial nonlinear model using specific nonlinear regressors.

sys1 = nlarx(z1,[4 2 1],'sigmoidnet','NonlinearRegressors',[1:3]);

Randomly perturb the model parameters to avoid local minima.

sys2 = init(sys1);

Estimate the new nonlinear model with the perturbed values.

sys2 = nlarx(z1,sys2);

Estimate Nonlinear ARX Model Using Specific Options

Load the estimation data.

load twotankdata;

Create an iddata object from the estimation data.

z = iddata(y,u,0.2);

Create an nlarxOptions option set specifying a simulation error minimization objective and a maximum of 50 estimation iterations.

opt = nlarxOptions;
opt.Focus = 'simulation';
opt.SearchOptions.MaxIterations = 50;

Estimate the nonlinear ARX model.

sys = nlarx(z,[4 4 1],'sigmoidnet',opt);

Estimate Regularized Nonlinear ARX Model with Large Number of Units

Load the regularization example data.

load regularizationExampleData.mat nldata;

Create a sigmoidnet nonlinearity with 30 units, and specify the model orders.

NL = sigmoidnet('NumberOfUnits',30);
Orders = [1 2 1];

Create an estimation option set and set the estimation search method to lm.

opt = nlarxOptions('SearchMethod','lm');

Estimate an unregularized model.

sys = nlarx(nldata,Orders,NL,opt);

Configure the regularization Lambda parameter.

opt.Regularization.Lambda = 1e-8;

Estimate a regularized model.

sysR = nlarx(nldata,Orders,NL,opt);

Compare the two models.

compare(nldata,sys,sysR)

The large negative fit result for the unregularized model indicates a poor fit to the data. Estimating a regularized model produces a significantly better result.

Input Arguments

Data — Time-domain estimation data

iddata object

Time-domain estimation data, specified as an iddata object. Data can have one or more output channels and zero or more input channels. Data must be uniformly sampled and cannot contain missing (NaN) samples.

Orders — Model orders and delays

1-by-3 vector of positive integers | 1-by-3 vector of matrices

Model orders and delays for defining the regressor configuration, specified as a 1-by-3 vector, [na nb nk].

For a model with n_y output channels and n_u input channels:

• na is an n_y-by-n_y matrix, where na(i,j) specifies the number of regressors from the jth output used to predict the ith output.

• nb is an n_y-by-n_u matrix, where nb(i,j) specifies the number of regressors from the jth input used to predict the ith output.

• nk is an n_y-by-n_u matrix, where nk(i,j) specifies the lag in the jth input used to predict the ith output.

na = [1 2; 2 3]
nb = [1 2 3; 2 3 1];
nk = [2 0 3; 1 0 5];

The estimation data for this system has three inputs (u1, u2, u3) and two outputs (y1, y2). Consider the regressors used to predict output, y2(t):

• Since na(2,:) is [2 3], the contributing regressors from the outputs are:

  • y1(t-1) and y1(t-2)

  • y2(t-1), y2(t-2), and y2(t-3)

• Since nb(2,:) is [2 3 1] and nk(2,:) is [1 0 5], the contributing regressors from the inputs are:

  • u1(t-1) and u1(t-2)

  • u2(t), u2(t-1), and u2(t-2)

  • u3(t-5)

Note

The minimum lag for regressors based on output variables is always 1, while the minimum lag for regressors based on input variables is dictated by nk. Use getreg to view the complete set of regressors used by the nonlinear ARX model.

Nonlinearity — Nonlinearity estimator

'wavenet' (default) | 'sigmoidnet' | 'treepartition' | 'linear' | nonlinearity estimator object | array of nonlinearity estimator objects

Nonlinearity estimator, specified as one of the following:

'wavenet' or wavenet object	Wavelet network
'sigmoidnet' or sigmoidnet object	Sigmoid network
'treepartition' or treepartition object	Binary-tree
'linear' or [] or linear object	Linear function
neuralnet object	Neural network — Requires Deep Learning Toolbox™.
customnet object	Custom network — Similar to sigmoidnet, but with a user-defined replacement for the sigmoid function.

For more information, see Available Nonlinearity Estimators for Nonlinear ARX Models.

Specifying a character vector, for example 'sigmoidnet', creates a nonlinearity estimator object with default settings. Alternatively, you can specify nonlinearity estimator settings in two ways:

• Use the associated nonlinearity estimator function with Name-Value pair arguments.

  NL = sigmoidnet('NumberOfUnits',10);

• Create and modify a default nonlinearity estimator object.

  NL = sigmoidnet;
  NL.NumberOfUnits = 10;

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

Example: 'sigmoidnet' specifies a sigmoid network nonlinearity with a default configuration.

Example: treepartition('NumberOfUnits',5) specifies a binary-tree nonlinearity with 5 terms in the binary tree expansion.

Example: [wavenet('NumberOfUnits',10);sigmoidnet] specifies different nonlinearity estimators for two output channels.

LinModel — Discrete time input-output polynomial model of ARX structure

idpoly model

Discrete time input-output polynomial model of ARX structure, specified as an idpoly model. Create this object using the idpoly constructor or estimate it using the arx command.

sys0 — Nonlinear ARX model

idnlarx model

Nonlinear ARX model, specified as an idnlarx model. sys0 can be:

• A model previously estimated using nlarx. The estimation algorithm uses the parameters of sys0 as initial guesses. In this case, use init to slightly perturb the model properties to avoid being trapped in local minima.

  sys = init(sys);
  sys = nlarx(data,sys);

• A model previously created using idnlarx and with properties set using dot notation. Use this method to avoid complicated Name-Value pair syntax when configuring multiple model properties. For example, use

  sys1 = idnlarx([4 3 1]);
  sys1.Nonlinearity = 'treepartition';
  sys1.CustomRegressors = {'sin(u1(t-1))'};
  sys1.NonlinearRegressors = 'search';
  sys2 = nlarx(data,sys1);

  in place of the equivalent

  sys2 = nlarx(data,[4,3,1],'treepartition','CustomRegressors',...
      {'sin(u1(t-1))'},'NonlinearRegressors','search');

Options — Estimation options

nlarxOptions option set

Estimation options for nonlinear ARX model identification, specified as an nlarxOptions option set.

Name-Value Pair Arguments

Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside quotes. You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

Example: 'NonlinearRegressors','output' specifies that only the regressors containing output variables are used as inputs to the nonlinear block of the model.

'TimeVariable' — Independent variable name

't' (default) | character vector

Independent variable name, specified as the comma-separated pair consisting of 'TimeVariable' and a character vector. For example, 't'.

'CustomRegressors' — Regressors constructed from combinations of inputs and outputs

{} (default) | cell array of character vectors | array of customreg objects

Regressors constructed from combinations of inputs and outputs, specified as the comma-separated pair consisting of 'CustomRegressors' and one of the following for single-output systems:

• Cell array of character vectors. For example:

  • {'y1(t-3)^3','y2(t-1)*u1(t-3)','sin(u3(t-2))'}

  Each character vector must represent a valid formula for a regressor contributing towards the prediction of the model output. The formula must be written using the input and output names and the time variable name as variables.

• Array of custom regressor objects, created using customreg or polyreg.

For a model with n_y outputs, specify an n_y-by-1 cell array of customreg object arrays or character arrays.

These regressors are in addition to the standard regressors based on Orders.

Example: 'CustomRegressors',{'y1(t-3)^3','y2(t-1)*u1(t-3)'}

Example: 'CustomRegressors',{'sin(u3(t-2))'}

'NonlinearRegressors' — Subset of regressors that enter as inputs to the nonlinear block of the model

'all' (default) | 'output' | 'input' | 'standard' | 'custom' | 'search' | vector of positive integers | [] | cell array

Subset of regressors that enter as inputs to the nonlinear block of the model, specified as the comma-separated pair consisting of 'NonlinearRegressors' and one of the following values:

• 'all' — All regressors

• 'output' — Regressors containing output variables

• 'input' — Regressors containing input variables

• 'standard' — Standard regressors

• 'custom' — Custom regressors

• 'search' — The estimation algorithm performs a search for the best regressor subset. This is useful when you want to reduce a large number of regressors entering the nonlinear function block of the nonlinearity estimator. This option must be applied to all output models simultaneously.

• [] — No regressors. This creates a linear-in-regressor model.

• Vector of regressor indices. To determine the number and order of regressors, use getreg.

For a model with multiple outputs, specify a cell array of n_y elements, where n_y is the number of output channels. For each output, specify one of the preceding options. Alternatively, to apply the same regressor subset to all model outputs, specify [] or any of the character vector options alone, for example 'standard'.

Example: 'NonlinearRegressors','search' performs a best regressor search for the only output of a single output model, or all of the outputs of a multiple output model.

Example: 'NonlinearReg','input' applies only input regressors to the inputs of the nonlinear function.

Example: 'NonlinearRegressors',{'input','output'} applies input regressors to the first output, and output regressors to the second output of a model with two outputs.

Output Arguments

sys — Nonlinear ARX model

idnlarx object

Nonlinear ARX model that fits the given estimation data, returned as an idnlarx object. This model is created using the specified model orders, nonlinearity estimator, and estimation options.

Information about the estimation results and options used is stored in the Report property of the model. The contents of Report depend upon the choice of nonlinearity and estimation focus you specified for nlarx. Report has the following fields:

Report Field	Description
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: Field Description FitPercent Normalized root mean squared error (NRMSE) measure of how well the response of the model fits the estimation data, expressed as the percentage fit = 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 nlarxOptions 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: Field Description 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: Field Description 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.

For more information on using Report, see Estimation Report.

Algorithms

Nonlinear ARX Model Structure

A nonlinear ARX model consists of model regressors and a nonlinearity estimator. The nonlinearity estimator comprises both linear and nonlinear functions that act on the model regressors to give the model output. This block diagram represents the structure of a nonlinear ARX model in a simulation scenario.

The software computes the nonlinear ARX model output y in two stages:

1. It computes regressor values from the current and past input values and past output data.

   In the simplest case, regressors are delayed inputs and outputs, such as u(t-1) and y(t-3). These kind of regressors are called standard regressors. You specify the standard regressors using the model orders and delay. For more information, see Nonlinear ARX Model Orders and Delay. You can also specify custom regressors, which are nonlinear functions of delayed inputs and outputs. For example, u(t-1)*y(t-3). To create a set of polynomial type regressors, use polyreg.

   By default, all regressors are inputs to both the linear and the nonlinear function blocks of the nonlinearity estimator. You can choose a subset of regressors as inputs to the nonlinear function block.

2. It maps the regressors to the model output using the nonlinearity estimator block. The nonlinearity estimator block can include linear and nonlinear blocks in parallel. For example:

   $$F(x)=L^T(x-r)+d+g\left(Q(x-r)\right)$$

   Here, x is a vector of the regressors, and r is the mean of the regressors x. $$L^T(x)+d$$ is the output of the linear function block and is affine when d ≠ 0. d is a scalar offset. $$g\left(Q(x-r)\right)$$ represents the output of the nonlinear function block. Q is a projection matrix that makes the calculations well conditioned. The exact form of F(x) depends on your choice of the nonlinearity estimator. You can select from available nonlinearity estimators, such as tree-partition networks, wavelet networks, and multilayer neural networks. You can also exclude either the linear or the nonlinear function block from the nonlinearity estimator.

   When estimating a nonlinear ARX model, the software computes the model parameter values, such as L, r, d, Q, and other parameters specifying g.

Resulting nonlinear ARX models are idnlarx objects that store all model data, including model regressors and parameters of the nonlinearity estimator. For more information about these objects, see Nonlinear Model Structures.

Extended Capabilities

Automatic Parallel Support

Accelerate code by automatically running computation in parallel using Parallel Computing Toolbox™.

Parallel computing support is available for estimation using the lsqnonlin search method (requires Optimization Toolbox™). To enable parallel computing, use nlarxOptions, set SearchMethod to 'lsqnonlin', and set SearchOptions.Advanced.UseParallel to true.

For example:

opt = nlarxOptions;
opt.SearchMethod = 'lsqnonlin';
opt.SearchOptions.Advanced.UseParallel = true;