Accelerating the pace of engineering and science

# armax

Estimate parameters of ARMAX model using time-domain data

## Syntax

sys = armax(data,[na nb nc nk])
sys = armax(data,[na nb nc nk],Name,Value)
sys = armax(data,init_sys)
sys = armax(data,___,opt)

## Description

 Note:   armax supports only time-domain data. For frequency-domain data, use oe.

sys = armax(data,[na nb nc nk]) returns an idpoly model, sys, with estimated parameters and covariance (parameter uncertainties). Estimates the parameters using the prediction-error method and specified polynomial orders.

sys = armax(data,[na nb nc nk],Name,Value) returns an idpoly model, sys, with additional options specified by one or more Name,Value pair arguments.

sys = armax(data,init_sys) estimates a polynomial model using the ARMAX structure polynomial model init_sys to configure the initial parameterization.

sys = armax(data,___,opt) specifies estimation options using the option set opt.

## Input Arguments

 data Estimation data. Specify data as an iddata object containing the time-domain input-output data. You cannot use frequency-domain data for estimating ARMAX models. [na nb nc nk] Polynomial orders. [na nb nc nk] define the polynomial orders of an ARMAX Model. na — Order of the polynomial A(q). Specify na as an Ny-by-Ny matrix of nonnegative integers. Ny is the number of outputs.nb — Order of the polynomial B(q) + 1.nb is an Ny-by-Nu matrix of nonnegative integers. Ny is the number of outputs and Nu is the number of inputs.nc — Order of the polynomial C(q).nc is a column vector of nonnegative integers of length Ny. Ny is the number of outputs.nk — Input-output delay expressed as fixed leading zeros of the B polynomial. Specify nk as an Ny-by-Nu matrix of nonnegative integers. Ny is the number of outputs and Nu is the number of inputs. init_sys Linear polynomial model that configures the initial parameterization of sys. init_sys must be an ARMAX model. You may obtain init_sys by either performing an estimation using measured data, or by direct construction. Use the Structure property of init_sys to configure initial guesses and constraints for A(q), B(q), and C(q). To specify an initial guess for, say, the A(q) term of init_sys, set init_sys.Structure.a.Value as the initial guess. To specify constraints for, say, the B(q) term of init_sys: set init_sys.Structure.b.Minimum to the minimum B(q) coefficient valuesset init_sys.Structure.b.Maximum to the maximum B(q) coefficient valuesset init_sys.Structure.b.Free to indicate which B(q) coefficients are free for estimation You can similarly specify the initial guess and constraints for the other polynomials. If opt is not specified, and init_sys was created by estimation, then the estimation options from init_sys.Report.OptionsUsed are used. opt Estimation options. opt is an options set that specifies estimation options, including: estimation objectivehandling of initial conditionsnumerical search method to be used in estimation Use armaxOptions to create the options 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 single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

 'InputDelay' Input delays. InputDelay is a numeric vector specifying a time delay for each input channel. Specify input delays in integer multiples of the sampling period Ts. For example, InputDelay = 3 means a delay of three sampling periods. For a system with Nu inputs, set InputDelay to an Nu-by-1 vector, where each entry is a numerical value representing the input delay for the corresponding input channel. You can also set InputDelay to a scalar value to apply the same delay to all channels. Default: 0 for all input channels 'ioDelay' Transport delays. ioDelay is a numeric array specifying a separate transport delay for each input/output pair. Specify transport delays as integers denoting delay of a multiple of the sampling period Ts. For a MIMO system with Ny outputs and Nu inputs, set ioDelay to a Ny-by-Nu array, where each entry is a numerical value representing 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. Useful as a replacement for the nk order, you can factor out max(nk-1,0) lags as the ioDelay value. Default: 0 for all input/output pairs 'IntegrateNoise' Logical vector specifying integrators in the noise channel. IntegrateNoise is a logical vector of length Ny, where Ny is the number of outputs. Setting IntegrateNoise to true for a particular output results in the model: $A\left(q\right)y\left(t\right)=B\left(q\right)u\left(t-nk\right)+\frac{C\left(q\right)}{1-{q}^{-1}}e\left(t\right)$ Where, $\frac{1}{1-{q}^{-1}}$ is the integrator in the noise channel, e(t). Use IntegrateNoise to create an ARIMA model. For example, ```load iddata9 z9; z9.y = cumsum(z9.y); %integrated data sys = armax(z9,[4 1],'IntegrateNoise',true); compare(z9,sys,10) %10-step ahead prediction``` Default: false(Ny,1) (Ny is the number of outputs.)

## Output Arguments

 sys Identified ARMAX structure polynomial model. sys is a discrete-time idpoly model, which encapsulates the estimated A, B and C polynomials and the parameter covariance information.

## Examples

expand all

### Estimate ARMAX Model Using Regularization

Estimate a regularized ARMAX model by converting a regularized ARX model.

```load regularizationExampleData.mat m0simdata;
```

Estimate an unregularized ARMAX model of order 15.

```m1 = armax(m0simdata(1:150), [30 30 30 1]);
```

Estimate a regularized ARMAX model by determining Lambda value by trial and error.

```opt = armaxOptions;
opt.Regularization.Lambda = 1;
m2 = armax(m0simdata(1:150),  [30 30 30 1], opt);
```

Obtain a lower-order ARMAX model by converting a regularized ARX model followed by order reduction.

```opt1 = arxOptions;
[L,R] = arxRegul(m0simdata(1:150), [30 30 1]);
opt1.Regularization.Lambda = L;
opt1.Regularization.R = R;
m0 = arx(m0simdata(1:150),[30 30 1],opt1);
mr = idpoly(balred(idss(m0),7));
```

Compare the model outputs against data.

```compare(m0simdata(150:end), m1, m2, mr, compareOptions('InitialCondition','z'));
```

### Specify Estimation Options

Estimate an ARMAX model from measured data and specify the estimation options.

Estimate an ARMAX model with simulation focus, using 'lm' as the search method and maximum number of search iterations set to 10.

```load twotankdata;
z = iddata(y,u,0.2);
opt = armaxOptions;
opt.Focus = 'simulation';
opt.SearchMethod = 'lm';
opt.SearchOption.MaxIter = 10;
opt.Display = 'on';
sys = armax(z, [2 2 2 1], opt);
```

The termination conditions for measured component of the model shown in the progress viewer is that the maximum number of iterations were reached.

To improve results, re-estimate the model using a greater value for MaxIter, or continue iterations on the previously estimated model as follows:

```sys2 = armax(z, sys);
compare(z, sys, sys2)
```

where sys2 refines the parameters of sys to improve the fit to data.

### Estimate an ARIMA Model

Estimate a 4th order ARIMA model for univariate time-series data.

```load iddata9;
z9.y = cumsum(z9.y); % integrated data
model = armax(z9, [4 1], 'IntegrateNoise', true);
compare(z9, model, 10) % 10-step ahead prediction
```

### Estimate ARMAX Models Iteratively

Estimate ARMAX models of varying orders iteratively from measured data.

Estimate ARMAX models of orders varying between 1 and 4 for dryer data.

```load dryer2;
z = iddata(y2,u2,0.08,'Tstart',0);
na = 2:4; nc = 1:2; nk = 0:2;
models = cell(1,18);
ct = 1;
for i = 1:3
na_ = na(i);
nb_ = na_;
for j = 1:2
nc_ = nc(j);
for k = 1:3
nk_ = nk(k);
models{ct} = armax(z, [na_, nb_, nc_, nk_]);
ct = ct+1;
end
end
end
```

Stack the estimated models and compare their simulated responses to estimation data z.

```models = stack(1,models{:});
compare(z,models)
```

## Alternatives

armax does not support continuous-time model estimation. Use tfest to estimate a continuous-time transfer function model, or ssest to estimate a continuous-time state-space model.

expand all

### ARMAX Model

The ARMAX model structure is

A more compact way to write the difference equation is

$A\left(q\right)y\left(t\right)=B\left(q\right)u\left(t-{n}_{k}\right)+C\left(q\right)e\left(t\right)$

where

• $y\left(t\right)$— Output at time $t$.

• ${n}_{a}$ — Number of poles.

• ${n}_{b}$ — Number of zeroes plus 1.

• ${n}_{c}$ — Number of C coefficients.

• ${n}_{k}$ — Number of input samples that occur before the input affects the output, also called the dead time in the system.

• $y\left(t-1\right)\dots y\left(t-{n}_{a}\right)$ — Previous outputs on which the current output depends.

• $u\left(t-{n}_{k}\right)\dots u\left(t-{n}_{k}-{n}_{b}+1\right)$ — Previous and delayed inputs on which the current output depends.

• $e\left(t-1\right)\dots e\left(t-{n}_{c}\right)$ — White-noise disturbance value.

The parameters na, nb, and nc are the orders of the ARMAX model, and nk is the delay. q is the delay operator. Specifically,

$A\left(q\right)=1+{a}_{1}{q}^{-1}+\dots +{a}_{{n}_{a}}{q}^{-{n}_{a}}$

$B\left(q\right)={b}_{1}+{b}_{2}{q}^{-1}+\dots +{b}_{{n}_{b}}{q}^{-{n}_{b}+1}$

$C\left(q\right)=1+{c}_{1}{q}^{-1}+\dots +{c}_{{n}_{c}}{q}^{-{n}_{c}}$

If data is a time series, which has no input channels and one output channel, then armax calculates an ARMA model for the time series

$A\left(q\right)y\left(t\right)=e\left(t\right)$

In this case

```orders = [na nc]
```

### ARIMAX Model

An ARIMAX model structure is similar to ARMAX, except that it contains an integrator in the noise source e(t):

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

If there are no inputs, this reduces to an ARIMA model:

$A\left(q\right)y\left(t\right)=\frac{1}{\left(1-{q}^{-1}\right)}e\left(t\right)$

### Tips

• Use the IntegrateNoise property to add integrators to the noise source.

### Algorithms

An iterative search algorithm minimizes a robustified quadratic prediction error criterion. The iterations are terminated either when the maximum number of iterations is reached, or when the expected improvement is less than the specified tolerance, or when a lower value of the criterion cannot be found. You can get information about the stopping criteria using sys.Report.Termination.

Use the armaxOptions option set to create and configure options affecting the estimation results. In particular, set the search algorithm attributes, such as MaxIter and Tolerance, using the 'SearchOption' property.

When you do not specify initial parameter values for the iterative search as an initial model, they are constructed in a special four-stage LS-IV algorithm.

The cutoff value for the robustification is based on the Advanced.ErrorThreshold estimation option and on the estimated standard deviation of the residuals from the initial parameter estimate. It is not recalculated during the minimization. By default, no robustification is performed; the default value of ErrorThreshold option is 0.

To ensure that only models corresponding to stable predictors are tested, the algorithm performs a stability test of the predictor. Generally, both $C\left(q\right)$ and $F\left(q\right)$ (if applicable) must have all zeros inside the unit circle.

Minimization information is displayed on the screen when the estimation option 'Display' is 'On' or 'Full'. With 'Display' ='Full', both the current and the previous parameter estimates are displayed in column-vector form, listing parameters in alphabetical order. Also, the values of the criterion function (cost) are given and the Gauss-Newton vector and its norm are also displayed. With 'Display' = 'On' only the criterion values are displayed.

## References

Ljung, L. System Identification: Theory for the User, Upper Saddle River, NJ, Prentice-Hal PTR, 1999. See chapter about computing the estimate.