armax

Estimate parameters of ARMAX model using time-domain data

Syntax

• sys = armax(data,[na nb nc nk])
  example
• sys = armax(data,[na nb nc nk],Name,Value)
  example

• sys = armax(data,init_sys)
  example

• sys = armax(data,___,opt)
  example

Description

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 a discrete-time linear model init_sys to configure the initial parameterization.

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

Input Arguments

data — Time-domain estimation dataiddata object

Time-domain estimation data, specified as an iddata object. You cannot use frequency-domain data for estimating ARMAX models.

[na nb nc nk] — Polynomial orders1-by-4 vector of positive integers | 1-by-4 vector of matrices

Polynomial orders of an ARMAX Model, specified as a 1-by-4 vector, [na nb nc nk].

For a model with Ny outputs and Nu inputs:

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

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

• nc is the order of the polynomial C(q), specified as a column vector of nonnegative integers of length Ny.

• nk is the input-output delay expressed as fixed leading zeros of the B polynomial.

  Specify nk as an Ny-by-Nu matrix of nonnegative integers.

init_sys — System for configuring initial parametrizationdiscrete-time linear model

System for configuring initial parametrization of sys, specified as a discrete-time linear model. You obtain init_sys by either performing an estimation using measured data or by direct construction using commands such as idpoly and idss.

If init_sys is an ARMAX model, armax uses the parameter values of init_sys as the initial guess for estimation. To configure initial guesses and constraints for A(q), B(q), and C(q), use the Structure property of init_sys. For example:

• To specify an initial guess for the A(q) term of init_sys, set init_sys.Structure.a.Value as the initial guess.

• To specify constraints for the B(q) term of init_sys:

  • set init_sys.Structure.b.Minimum to the minimum B(q) coefficient values.

  • set init_sys.Structure.b.Maximum to the maximum B(q) coefficient values.

  • set init_sys.Structure.b.Free to indicate which B(q) coefficients are free for estimation.

If init_sys is not a polynomial model of ARMAX structure, the software first converts init_sys to an ARMAX model. armax uses the parameters of the resulting model as the initial guess for estimating sys.

If opt is not specified, and init_sys was obtained by estimation, then the estimation options from init_sys.Report.OptionsUsed are used.

opt — Estimation optionsarmaxOptions option set

Estimation options for ARMAX model identification, specified as an armaxOptions 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 single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

'InputDelay' — Input delays0 (default) | scalar | vector of positive integers

Input delays, specified as the comma-separated pair consisting of 'InputDelay' and one of the following:

• Nu-by-1 vector, where Nu is the number of inputs — Each entry is a numerical value representing the input delay for the corresponding input channel. Specify input delays as integer multiples of the sample time Ts. For example, InputDelay = 3 means a delay of three sampling periods.

• Scalar value — The same delay is applied to all input channels.

'ioDelay' — Transport delays0 (default) | scalar | matrix

Transport delays for each input/output pair, specified as the comma-separated pair consisting of 'ioDelay' and one of the following:

• Ny-by-Nu matrix, where Ny is the number of outputs and Nu is the number of inputs — Each entry is an integer value representing the transport delay for the corresponding input/output pair. Specify transport delays as integers denoting delay of a multiple of the sample time Ts.

• Scalar value — The same delay is applied to all input/output pairs.

'ioDelay' is useful as a replacement for the nk order. You can factor out max(nk-1,0) lags as the ioDelay value. For nk>1, armax(na,nb,nk) is equivalent to armax(na,nb,1,'ioDelay',nk-1).

'IntegrateNoise' — Addition of integrators in noise channelfalse(Ny,1), where Ny is the number of outputs (default) | logical vector

Addition of integrators in noise channel, specified as the comma-separated pair consisting of ‘IntegrateNoise' and 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(q)y(t)=B(q)u(t-nk)+\frac{C(q)}{1-q^{-1}}e(t)$$

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

Use IntegrateNoise to create ARIMA or ARIMAX models.

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

Output Arguments

sys — ARMAX modelidpoly object

ARMAX model that fits the given 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 a string with one of the following values: 'zero' — The initial conditions were set to zero. 'estimate' — The initial conditions were treated as independent estimation parameters. 'backcast' — The initial conditions were estimated using the best least squares fit. 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 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 a percentage. 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 armaxOptions 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 in the MATLAB® documentation.
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. 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'. LastImprovement Criterion improvement in the last iteration, expressed as a percentage. Omitted when the search method is 'lsqnonlin'. Algorithm Algorithm used by 'lsqnonlin' 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.

Examples

Estimate ARMAX Model Using Regularization

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

Load data.

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.

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

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 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)

Estimate ARMAX Model Using State-Space Model as Inital Guess

Load the estimation data.

load iddata2 z2

Estimate a state-space model of order 3 from the estimation data.

sys0 = n4sid(z2,3);

Estimate an ARMAX model using the previously estimated state-space model as an initial guess.

sys = armax(z2,sys0);

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.

armax supports only time-domain data. For frequency-domain data, use oe to estimate an Output-Error (OE) model.

More About

ARMAX Model

The ARMAX model structure is:

$$\begin{array}{l}y(t)+a_{1}y(t-1)+\dots +a_{n_a}y(t-n_a)=\\ b_{1}u(t-n_k)+\dots +b_{n_b}u(t-n_k-n_b+1)+\\ c_{1}e(t-1)+\dots +c_{n_c}e(t-n_c)+e(t)\end{array}$$

A more compact way to write the difference equation is:

$$A(q)y(t)=B(q)u(t-n_k)+C(q)e(t)$$

where,

• $$y(t)$$— 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(t-1)\dots y(t-n_a)$$ — Previous outputs on which the current output depends.

• $$u(t-n_k)\dots u(t-n_k-n_b+1)$$ — Previous and delayed inputs on which the current output depends.

• $$e(t-1)\dots e(t-n_c)$$ — 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(q)=1+a_1{q}^{-1}+\dots +a_{n_a}{q}^{-n_a}$$

$$B(q)=b_1+b_2{q}^{-1}+\dots +b_{n_b}{q}^{-n_b+1}$$

$$C(q)=1+c_1{q}^{-1}+\dots +c_{n_c}{q}^{-n_c}$$

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

$$A(q)y(t)=C(q)e(t)$$

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(q)y(t)=B(q)u(t-nk)+\frac{C(q)}{(1-q^{-1})}e(t)$$

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

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

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 when any of the following is true:

• Maximum number of iterations is reached.

• Expected improvement is less than the specified tolerance.

• 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(q)$$ and $$F(q)$$ (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.

• Loss Function and Model Quality Metrics
• Regularized Estimates of Model Parameters
• Estimation Report

References

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