iv4

ARX model estimation using four-stage instrumental variable method

Syntax

sys = iv4(tt,[na nb nk])

sys = iv4(u,y,[na nb nk])

sys = iv4(data,[na nb nk])

sys = iv4(___,Name,Value)

sys = iv4(___,opt)

[sys,ic] = iv4(___)

Description

Estimate ARX Polynomial Model

sys = iv4(tt,[na nb nk]) estimates an ARX polynomial model sys using the time-domain data in the timetable tt. [na nb nk] specifies the ARX structure orders of the A and B polynomials and the input-to-output delay.

The software estimates sys using the four-stage instrumental variable method. The estimation algorithm is insensitive to the color of the noise term.

sys is an ARX model, which has the following form.

$$A(q)y(t)=B(q)u(t-nk)+v(t)$$

For more details on the ARX model structure, see arx.

sys = iv4(u,y,[na nb nk]) uses the time-domain input and output signals in the comma-separated matrices u,y. The software assumes that the data sample time is one second. To change the sample time, set Ts using name-value syntax.

sys = iv4(data,[na nb nk]) uses the time-domain or frequency-domain data in the data object data. Use this syntax especially when you want to estimate a model using frequency-domain or frequency-response data, or when you want to take advantage of the additional information, such as data sample time or experiment labeling, that data objects provide.

Specify Additional Options

sys = iv4(___,Name,Value) specifies additional options using one or more name-value arguments.

You can use this syntax with any of the previous input-argument combinations.

sys = iv4(___,opt) uses the option set opt to configure the estimation behavior.

Return Estimated Initial Conditions

[sys,ic] = iv4(___) returns the estimated initial conditions as an initialCondition object. For more information on ic, see the ic argument description.

Use this syntax if you plan to simulate or predict the model response using the same estimation input data and then compare the response with the same estimation output data. Incorporating the initial conditions yields a better match during the first part of the simulation.

Examples

Estimate ARX Model Using Four-Stage Instrumental Variable Method

Load estimation data.

load sdata7.mat tt7

This data has two inputs, u1 and u2, and one output, y. Visualize y.

T = tt7.Properties.RowTimes;
plot(T,tt7.y)

The plot shows a disturbance at around 314 seconds that should not be included in the estimation data.

Specify the ARX model orders, using the same orders for both inputs.

na = 4;
nb = [2 2];

Specify a delay of 3 samples for input u1 and a delay of 1 sample for input u2.

nk = [3 1];

Estimate an ARX model using the four-stage instrumental variable method and the first 300 samples of data.

sys = iv4(tt7((1:300),:),[na nb nk])

sys =
Discrete-time ARX model: A(z)y(t) = B(z)u(t) + e(t)                
  A(z) = 1 - 0.1071 z^-1 - 0.9286 z^-2 + 0.07906 z^-3 + 0.4891 z^-4
                                                                   
  B1(z) = 1.193 z^-3 + 1.241 z^-4                                  
                                                                   
  B2(z) = 0.655 z^-1 - 0.7329 z^-2                                 
                                                                   
Sample time: 1 seconds
  
Parameterization:
   Polynomial orders:   na=4   nb=[2 2]   nk=[3 1]
   Number of free coefficients: 8
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using IV4 on time domain data.         
Fit to estimation data: 64.76% (prediction focus)
FPE: 1.33, MSE: 1.228                            
 

Input Arguments

tt — Timetable-based estimation data

timetable | cell array of timetables

Timetable-based estimation data, specified as a timetable that uses a regularly spaced time vector. tt contains variables representing input and output channels.

For multiexperiment data, specify tt as an Ne-by-1 cell array of timetables, where Ne is the number of experiments. The sample times of all the experiments must match.

To select individual input and output channels to use for estimation, use the InputName and OutputName name-value arguments.

For example, use the following command to select timetable variables "u1" and "u3" as inputs and the variables "y2" and "y4" as outputs.

sys = iv4(tt,__,'InputName',["u1" "u3"],'OutputName',["y2" "y4"])

For more information about working with estimation data types, see Data Domains and Data Types in System Identification Toolbox.

u, y — Matrix-based estimation data

matrices | cell array of matrices

Estimation data, specified for SISO systems as a comma-separated pair of N_s-element numeric column vectors that contain uniformly sampled input and output time-domain signal values. Here, N_s is the number of samples.

For MIMO systems, specify u,y as an input/output matrix pair with the following dimensions:

• u — N_s-by-N_u, where N_u is the number of inputs

• y — N_s-by-N_y, where N_y is the number of outputs

For multiexperiment data, specify u,y as a pair of 1-by-N_e cell arrays, where N_e is the number of experiments. The sample times of all the experiments must match.

For more information about working with estimation data types, see Data Domains and Data Types in System Identification Toolbox.

data — Estimation data object

iddata object | frd object | idfrd object

Estimation data, specified as an iddata object, an frd (Control System Toolbox) object, or an idfrd frequency-response object. For frequency-domain data, data must be discrete time, that is, Ts > 0, and there must be only one output.

For more information about working with estimation data types, see Data Domains and Data Types in System Identification Toolbox.

[na nb nk] — ARX polynomial orders and delay

row vector of nonnegative integers | row vector of nonnegative integer matrices

ARX polynomial orders and delay, specified as a row vector of integers or integer matrices.

• For SISO models, [na nb nk] is a 1-by-3 row vector of nonnegative integers that correspond to the orders of the A(q) and B(q) polynomials and the input-to-output delays, respectively. The polynomial order is equal to the number of coefficients to estimate in that polynomial.

• For a model with N_y outputs and N_u inputs:

  • na is the order of polynomial A(q), specified as an N_y-by-N_y matrix of nonnegative integers.

  • nb is the order of polynomial B(q) + 1, specified as an N_y-by-N_u matrix of nonnegative integers.

  • nk is the input-output delay, also known as the transport delay, specified as an N_y-by-N_u matrix of nonnegative integers.

For more details on the ARX model structure, see arx.

Example: iv4(data,[2 1 1]) estimates, from an iddata object, a second-order ARX model with one input and one sample delay.

opt — Estimation options

iv4Options option set

Estimation options for ARX model identification, specified as an iv4Options option set. Options specified by opt include the following:

• Initial condition handling

• Input and output data offsets

• Regularization

For more information, see iv4Options.

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'IntegrateNoise',true adds an integrator in the noise sources

InputDelay — Delay at inputs

0 (default) | scalar | vector

Delay at each input, specified as a scalar or a vector. For a system with Nu inputs, set InputDelay to an Nu-by-1 vector. Each entry of this vector is a numerical value that represents the input delay for the corresponding input channel. For continuous-time models, specify input delays in the time unit stored in the TimeUnit property of the model object. For discrete-time models, specify input delays in integer multiples of the sample time Ts. For example, InputDelay = 3 means a delay of three sample times.

Set InputDelay to a scalar value to apply the same delay to all channels.

IODelay — Transport Delays

0 (default) | positive integer | array of positive integers

Transport delays, specified as positive integer or an array of positive integers. IODelay specifies a separate transport delay for each input/output pair.

For continuous-time systems, specify transport delays in the time unit stored in the TimeUnit property. For discrete-time systems, specify transport delays in integer multiples of the sample time, Ts.

For a MIMO system with Ny outputs and Nu inputs, set IODelay to an Ny-by-Nu array. Each entry of this array is a numerical value that represents the transport delay for the corresponding input/output pair. You can also set IODelay to a scalar value to apply the same delay to all input/output pairs.

IntegrateNoise — Addition of integrators to noise channel

false (default) | logical vector

Option whether to add integrators to the noise channel, specified as a logical vector of length Ny, where Ny is the number of outputs.

Adding an integrator creates an ARIX model represented by

$$A(q)y(t)=B(q)u(t-nk)+\frac{1}{1-q^{-1}}e(t)$$

where $$\frac{1}{1-q^{-1}}$$ is the integrator in the noise channel, e(t).

InputName — Input channel names

string | character vector | string array | cell array of character vectors

Input channel names, specified as a string, character vector, string array, or cell array of character vectors.

If you are using a timetable for the data source, the names in InputName must be a subset of the timetable variables.

Example: sys = iv4(tt,__,'InputName',["u1" "u2"]) selects the variables u1 and u2 as the input channels from the timetable tt to use for the estimation.

OutputName — Output channel names

string | character vector | string array | cell array of character vectors

Output channel names, specified as a string, character vector, string array, or cell array of character vectors.

If you are using a timetable for the data source, the names in OutputName must be a subset of the timetable variables.

Example: sys = iv4(tt,__,'OutputName',["y1" "y3"]) selects the variables y1 and y3 as the output channels from the timetable tt to use for the estimation.

Output Arguments

sys — ARX model

idpoly object

ARX model that fits the estimation data, returned as a discrete-time idpoly object. This model is created using the specified model orders, delays, 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 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
InitialCondition	Handling of initial conditions during model estimation, returned as one of the following values: 'zero' — The initial conditions were set to zero. 'estimate' — The initial conditions were treated as independent estimation parameters. This field is especially useful to view how the initial conditions were handled when the InitialCondition option in the estimation option set is 'auto'.
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 these fields. 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 iv4Options 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. 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' — A zero-order hold maintains a piecewise-constant input signal between samples. 'foh' — A 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 [].

For more information on using Report, see Estimation Report.

ic — Initial conditions

initialCondition object | object array of initialCondition values

Estimated initial conditions, returned as an initialCondition object or an object array of initialCondition values.

• For a single-experiment data set, ic represents, in state-space form, the free response of the transfer function model (A and C matrices) to the estimated initial states (x₀).

• For a multiple-experiment data set with N_e experiments, ic is an object array of length N_e that contains one set of initialCondition values for each experiment.

For more information, see initialCondition.

Algorithms

Estimation is performed in 4 stages. The first stage uses the arx function. The resulting model generates the instruments for a second-stage IV estimate. The residuals obtained from this model are modeled as a high-order AR model. At the fourth stage, the input-output data is filtered through this AR model and then subjected to the IV function with the same instrument filters as in the second stage.

For the multiple-output case, optimal instruments are obtained only if the noise sources at the different outputs have the same color. The estimates obtained with the routine are reasonably accurate, however, even in other cases.

References

[1] Ljung, Lennart. System Identification: Theory for the User, equations (15.21) through (15.26). 2nd ed. Prentice Hall Information and System Sciences Series. Upper Saddle River, NJ: Prentice Hall PTR, 1999.

Version History

Introduced before R2006a