Documentation |
A polynomial model uses a generalized notion of transfer functions to express the relationship between the input, u(t), the output y(t), and the noise e(t) using the equation:
$$A(q)y(t)={\displaystyle \sum _{i=1}^{nu}\frac{{B}_{i}(q)}{{F}_{i}(q)}{u}_{i}\left(t-n{k}_{i}\right)}+\frac{C(q)}{D(q)}e(t)$$
The variables A, B, C, D, and F are polynomials expressed in the time-shift operator q^-1. u_{i} is the ith input, nu is the total number of inputs, and nk_{i} is the ith input delay that characterizes the transport delay. The variance of the white noise e(t) is assumed to be $$\lambda $$. For more information about the time-shift operator, see Understanding the Time-Shift Operator q.
In practice, not all the polynomials are simultaneously active. Often, simpler forms, such as ARX, ARMAX, Output-Error, and Box-Jenkins are employed. You also have the option of introducing an integrator in the noise source so that the general model takes the form:
$$A(q)y(t)={\displaystyle \sum _{i=1}^{nu}\frac{{B}_{i}(q)}{{F}_{i}(q)}{u}_{i}\left(t-n{k}_{i}\right)}+\frac{C(q)}{D(q)}\frac{1}{1-{q}^{-1}}e(t)$$
For more information, see Different Configurations of Polynomial Models.
You can estimate polynomial models using time or frequency domain data.
For estimation, you must specify the model order as a set of integers that represent the number of coefficients for each polynomial you include in your selected structure—na for A, nb for B, nc for C, nd for D, and nf for F. You must also specify the number of samples nk corresponding to the input delay—dead time—given by the number of samples before the output responds to the input.
The number of coefficients in denominator polynomials is equal to the number of poles, and the number of coefficients in the numerator polynomials is equal to the number of zeros plus 1. When the dynamics from u(t) to y(t) contain a delay of nk samples, then the first nk coefficients of B are zero.
For more information about the family of transfer-function models, see the corresponding section in System Identification: Theory for the User, Second Edition, by Lennart Ljung, Prentice Hall PTR, 1999.
The general polynomial equation is written in terms of the time-shift operator q^{–1}. To understand this time-shift operator, consider the following discrete-time difference equation:
$$\begin{array}{l}y(t)+{a}_{1}y(t-T)+{a}_{2}y(t-2T)=\\ \text{}{b}_{1}u(t-T)+{b}_{2}u(t-2T)\end{array}$$
where y(t) is the output, u(t) is the input, and T is the sampling interval. q^{-1} is a time-shift operator that compactly represents such difference equations using $${q}^{-1}u(t)=u(t-T)$$:
$$\begin{array}{l}y(t)+{a}_{1}{q}^{-1}y(t)+{a}_{2}{q}^{-2}y(t)=\\ {\text{b}}_{\text{1}}{q}^{-1}u(t)+{b}_{2}{q}^{-2}u(t)\\ \text{or}\\ A(q)y(t)=B(q)u(t)\end{array}$$
In this case, $$A(q)=1+{a}_{1}{q}^{-1}+{a}_{2}{q}^{-2}$$ and $$B(q)={b}_{1}{q}^{-1}+{b}_{2}{q}^{-2}$$.
These model structures are subsets of the following general polynomial equation:
$$A(q)y(t)={\displaystyle \sum _{i=1}^{nu}\frac{{B}_{i}(q)}{{F}_{i}(q)}{u}_{i}\left(t-n{k}_{i}\right)}+\frac{C(q)}{D(q)}e(t)$$
The model structures differ by how many of these polynomials are included in the structure. Thus, different model structures provide varying levels of flexibility for modeling the dynamics and noise characteristics.
The following table summarizes common linear polynomial model structures supported by the System Identification Toolbox™ product. If you have a specific structure in mind for your application, you can decide whether the dynamics and the noise have common or different poles. A(q) corresponds to poles that are common for the dynamic model and the noise model. Using common poles for dynamics and noise is useful when the disturbances enter the system at the input. F _{i} determines the poles unique to the system dynamics, and D determines the poles unique to the disturbances.
Model Structure | Equation | Description |
---|---|---|
ARX | $$A(q)y(t)={\displaystyle \sum _{i=1}^{nu}{B}_{i}(q){u}_{i}\left(t-n{k}_{i}\right)}+e(t)$$ | The noise model is $${\scriptscriptstyle \frac{1}{A}}$$ and the noise is coupled to the dynamics model. ARX does not let you model noise and dynamics independently. Estimate an ARX model to obtain a simple model at good signal-to-noise ratios. |
ARIX | $$Ay=Bu+\frac{1}{1-{q}^{-1}}e$$ | Extends the ARX structure by including an integrator in the noise source, e(t). This is useful in cases where the disturbance is not stationary. |
ARMAX | $$A(q)y(t)={\displaystyle \sum _{i=1}^{nu}{B}_{i}(q){u}_{i}\left(t-n{k}_{i}\right)}+C(q)e(t)$$ | Extends the ARX structure by providing more flexibility for modeling noise using the C parameters (a moving average of white noise). Use ARMAX when the dominating disturbances enter at the input. Such disturbances are called load disturbances. |
ARIMAX | $$Ay=Bu+C\frac{1}{1-{q}^{-1}}e$$ | Extends theARMAX structure by including an integrator in the noise source, e(t). This is useful in cases where the disturbance is not stationary. |
Box-Jenkins (BJ) | $$y(t)={\displaystyle \sum _{i=1}^{nu}\frac{{B}_{i}(q)}{{F}_{i}(q)}{u}_{i}\left(t-n{k}_{i}\right)}+\frac{C(q)}{D(q)}e(t)$$ | Provides completely independent parameterization for the dynamics
and the noise using rational polynomial functions. Use BJ models when the noise does not enter at the input, but is primary a measurement disturbance, This structure provides additional flexibility for modeling noise. |
Output-Error (OE) | $$y(t)={\displaystyle \sum _{i=1}^{nu}\frac{{B}_{i}(q)}{{F}_{i}(q)}{u}_{i}\left(t-n{k}_{i}\right)}+e(t)$$ | Use when you want to parameterize dynamics, but do not want to estimate a noise model. |
The polynomial models can contain one or more outputs and zero or more inputs.
The System Identification app supports direct estimation of ARX, ARMAX, OE and BJ models. You can add a noise integrator to the ARX, ARMAX and BJ forms. However, you can use polyest to estimate all five polynomial or any subset of polynomials in the general equation. For more information about working with pem, see Using polyest to Estimate Polynomial Models.
In continuous time, the general frequency-domain equation is written in terms of the Laplace transform variable s, which corresponds to a differentiation operation:
$$A(s)Y(s)=\frac{B(s)}{F(s)}U(s)+\frac{C(s)}{D(s)}E(s)$$
In the continuous-time case, the underlying time-domain model is a differential equation and the model order integers represent the number of estimated numerator and denominator coefficients. For example, n_{a}=3 and n_{b}=2 correspond to the following model:
$$\begin{array}{l}A(s)={s}^{4}+{a}_{1}{s}^{3}+{a}_{2}{s}^{2}+{a}_{3}\\ B(s)={b}_{1}s+{b}_{2}\end{array}$$
You can only estimate continuous-time polynomial models directly using continuous-time frequency-domain data. In this case, you must set the Ts data property to 0 to indicate that you have continuous-time frequency-domain data, and use the oe command to estimate an Output-Error polynomial model. Continuous-time models of other structures such as ARMAX or BJ cannot be estimated. You can obtain those forms only by direct construction (using idpoly), conversion from other model types, or by converting a discrete-time model into continuous-time (d2c). Note that the OE form represents a transfer function expressed as a ratio of numerator (B) and denominator (F) polynomials. For such forms consider using the transfer function models, represented by idtf models. You can estimate transfer function models using both time and frequency domain data. In addition to the numerator and denominator polynomials, you can also estimate transport delays. See idtf and tfest for more information.
You can create multi-output polynomial models by using the idpoly command or estimate them using ar, arx, bj, oe, armax, and polyest. In the app, you can estimate such models by choosing a multi-output data set and setting the orders appropriately in the Polynomial Models dialog box. For more details on the orders of multi-output models, see Polynomial Sizes and Orders of Multi-Output Polynomial Models.
You can estimate linear, black-box polynomial models from data with the following characteristics:
Time- or frequency-domain data (iddata or idfrd data objects).
To estimate polynomial models for time-series data, see Time-Series Model Identification.
Real data or complex data in any domain.
Single-output and multiple-output.
You must import your data into the MATLAB^{®} workspace, as described in Data Preparation.
To get a linear, continuous-time model of arbitrary structure for time-domain data, you can estimate a discrete-time model, and then use d2c to transform it to a continuous-time model.
For continuous-time frequency-domain data, you can estimate directly only Output-Error (OE) continuous-time models. Other structures include noise models, which is not supported for frequency-domain data.
Tip To denote continuous-time frequency-domain data, set the data sampling interval to 0. You can set the sampling interval when you import data into the app or set the Ts property of the data object at the command line. |
You can estimate arbitrary-order, linear state-space models for both time- or frequency-domain data.
Set the data property Ts to:
0, for frequency response data that is measured directly from an experiment.
Equal to the Ts of the original data, for frequency response data obtained by transforming time-domain iddata (using spa and etfe).
To estimate polynomial models, you must provide input delays and model orders. If you already have insight into the physics of your system, you can specify the number of poles and zeros.
In most cases, you do not know the model orders in advance. To get initial model orders and delays for your system, you can estimate several ARX models with a range of orders and delays and compare the performance of these models. You choose the model orders that correspond to the best model performance and use these orders as an initial guess for further modeling.
Because this estimation procedure uses the ARX model structure, which includes the A and B polynomials, you only get estimates for the na, nb, and nk parameters. However, you can use these results as initial guesses for the corresponding polynomial orders and input delays in other model structures, such as ARMAX, OE, and BJ.
If the estimated nk is too small, the leading nb coefficients are much smaller than their standard deviations. Conversely, if the estimated nk is too large, there is a significant correlation between the residuals and the input for lags that correspond to the missing B terms. For information about residual analysis plots, see Residual Analysis.
The following procedure assumes that you have already imported your data into the app and performed any necessary preprocessing operations. For more information, see Represent Data.
To estimate model orders and input delays in the System Identification app:
In the System Identification app, select Estimate > Polynomial Models to open the Polynomials Models dialog box.
The ARX model is already selected by default in the Structure list.
Edit the Orders field to specify a range of poles, zeros, and delays. For example, enter the following values for na, nb, and nk:
[1:10 1:10 1:10]
Click Estimate to open the ARX Model Structure Selection window, which displays the model performance for each combination of model parameters. The following figure shows an example plot.
Select a rectangle that represents the optimum parameter combination and click Insert to estimates a model with these parameters. For information about using this plot, see Selecting Model Orders from the Best ARX Structure.
This action adds a new model to the Model Board in the System Identification app. The default name of the parametric model contains the model type and the number of poles, zeros, and delays. For example, arx692 is an ARX model with n_{a}=6, n_{b}=9, and a delay of two samples.
Click Close to close the ARX Model Structure Selection window.
After estimating model orders and delays, use these values as initial guesses for estimating other model structures, as described in Estimate Polynomial Models in the App.
You can estimate model orders using the struc, arxstruc, and selstruc commands in combination.
If you are working with a multiple-output system, you must use the struc, arxstruc, and selstruc commands one output at a time. You must subreference the correct output channel in your estimation and validation data sets.
For each estimation, you use two independent data sets—an estimation data set and a validation data set. These independent data set can be from different experiments, or data subsets from a single experiment. For more information about subreferencing data, see Select Data Channels, I/O Data and Experiments in iddata Objects and Select I/O Channels and Data in idfrd Objects.
For an example of estimating model orders for a multiple-input system, see Estimating Delays in the Multiple-Input System in System Identification Toolbox Getting Started Guide.
struc. The struc command creates a matrix of possible model-order combinations for a specified range of n_{a}, n_{b}, and n_{k} values.
For example, the following command defines the range of model orders and delays na=2:5, nb=1:5, and nk=1:5:
NN = struc(2:5,1:5,1:5))
arxstruc. The arxstruc command takes the output from struc, estimates an ARX model for each model order, and compares the model output to the measured output. arxstruc returns the loss for each model, which is the normalized sum of squared prediction errors.
For example, the following command uses the range of specified orders NN to compute the loss function for single-input/single-output estimation data data_e and validation data data_v:
V = arxstruc(data_e,data_v,NN)
Each row in NN corresponds to one set of orders:
[na nb nk]
selstruc. The selstruc command takes the output from arxstruc and opens the ARX Model Structure Selection window to guide your choice of the model order with the best performance.
For example, to open the ARX Model Structure Selection window and interactively choose the optimum parameter combination, use the following command:
selstruc(V)
For more information about working with the ARX Model Structure Selection window, see Selecting Model Orders from the Best ARX Structure.
To find the structure that minimizes Akaike's Information Criterion, use the following command:
nn = selstruc(V,'AIC')
where nn contains the corresponding na, nb, and nk orders.
Similarly, to find the structure that minimizes the Rissanen's Minimum Description Length (MDL), use the following command:
nn = selstruc(V,'MDL')
To select the structure with the smallest loss function, use the following command:
nn = selstruc(V,0)
After estimating model orders and delays, use these values as initial guesses for estimating other model structures, as described in Using polyest to Estimate Polynomial Models.
The delayest command estimates the time delay in a dynamic system by estimating a low-order, discrete-time ARX model and treating the delay as an unknown parameter.
By default, delayest assumes that n_{a}=n_{b}=2 and that there is a good signal-to-noise ratio, and uses this information to estimate n_{k}.
To estimate the delay for a data set data, type the following at the prompt:
delayest(data)
If your data has a single input, MATLAB computes a scalar value for the input delay—equal to the number of data samples. If your data has multiple inputs, MATLAB returns a vector, where each value is the delay for the corresponding input signal.
To compute the actual delay time, you must multiply the input delay by the sampling interval of the data.
You can also use the ARX Model Structure Selection window to estimate input delays and model order together, as described in Estimating Model Orders at the Command Line.
You generate the ARX Model Structure Selection window for your data to select the best-fit model.
For a procedure on generating this plot in the System Identification app, see Estimating Orders and Delays in the App. To open this plot at the command line, see Estimating Model Orders at the Command Line.
The following figure shows a sample plot in the ARX Model Structure Selection window.
You use this plot to select the best-fit model.
The horizontal axis is the total number of parameters — n_{a} + n_{b}.
The vertical axis, called Unexplained output variance (in %), is the portion of the output not explained by the model—the ARX model prediction error for the number of parameters shown on the horizontal axis.
The prediction error is the sum of the squares of the differences between the validation data output and the model one-step-ahead predicted output.
n_{k} is the delay.
Three rectangles are highlighted on the plot in green, blue, and red. Each color indicates a type of best-fit criterion, as follows:
Red — Best fit minimizes the sum of the squares of the difference between the validation data output and the model output. This rectangle indicates the overall best fit.
Green — Best fit minimizes Rissanen MDL criterion.
Blue — Best fit minimizes Akaike AIC criterion.
In the ARX Model Structure Selection window, click any bar to view the orders that give the best fit. The area on the right is dynamically updated to show the orders and delays that give the best fit.
For more information about the AIC criterion, see Akaike's Criteria for Model Validation.
In the System Identification app, select Estimate > Polynomial Models to open the Polynomial Models dialog box.
For more information on the options in the dialog box, click Help.
In the Structure list, select the polynomial model structure you want to estimate from the following options:
ARX:[na nb nk]
ARMAX:[na nb nc nk]
OE:[nb nf nk]
BJ:[nb nc nd nf nk]
This action updates the options in the Polynomial Models dialog box to correspond with this model structure. For information about each model structure, see What Are Polynomial Models?.
Note: For time-series data, only AR and ARMA models are available. For more information about estimating time-series models, see Time-Series Model Identification. |
In the Orders field, specify the model orders and delays, as follows:
For single-output polynomial models. Enter the model orders and delays according to the sequence displayed in the Structure field. For multiple-input models, specify nb and nk as row vectors with as many elements as there are inputs. If you are estimating BJ and OE models, you must also specify nf as a vector.
For example, for a three-input system, nb can be [1 2 4], where each element corresponds to an input.
For multiple-output models. Enter the model orders, as described in Polynomial Sizes and Orders of Multi-Output Polynomial Models.
(ARX models only) Select the estimation Method as ARX or IV (instrumental variable method). For information about the algorithms, see Polynomial Model Estimation Algorithms.
(ARX, ARMAX, and BJ models only) Check the Add noise integration check box to add an integrator to the noise source, e.
Specify the delay using the Input delay edit box. The value must be a vector of length equal to the number of input channels in the data. For discrete-time estimations (any estimation using data with nonzero sample-time), the delay must be expressed in the number of lags. These delays are separate from the "in-model" delays specified by the nk order in the Orders edit box.
In the Name field, edit the name of the model or keep the default.
In the Focus list, select how to weigh the relative importance of the fit at different frequencies. For more information about each option, see Assigning Estimation Weightings.
In the Initial state list, specify how you want the algorithm to treat initial conditions. For more information about the available options, see Specifying Initial Conditions for Iterative Estimation Algorithms.
In the Covariance list, select Estimate if you want the algorithm to compute parameter uncertainties. Effects of such uncertainties are displayed on plots as model confidence regions.
To omit estimating uncertainty, select None. Skipping uncertainty computation for large, multiple-output models might reduce computation time.
Click Regularization to obtain regularized estimates of model parameters. Specify the regularization constants in the Regularization Options dialog box. To learn more, see Regularized Estimates of Model Parameters.
(ARMAX, OE, and BJ models only) To view the estimation progress in the MATLAB Command Window, select the Display progress check box. This launches a progress viewer window in which estimation progress is reported.
Click Estimate to add this model to the Model Board in the System Identification app.
(Prediction-error method only) To stop the search and save the results after the current iteration has been completed, click Stop Iterations. To continue iterations from the current model, click the Continue iter button to assign current parameter values as initial guesses for the next search.
Validate the model by selecting the appropriate check box in the Model Views area of the System Identification app. For more information about validating models, see Validating Models After Estimation.
Export the model to the MATLAB workspace for further analysis by dragging it to the To Workspace rectangle in the System Identification app.
Tip For ARX and OE models, you can use the exported model for initializing a nonlinear estimation at the command line. This initialization may improve the fit of the model. See Using Linear Model for Nonlinear ARX Estimation, and Using Linear Model for Hammerstein-Wiener Estimation. |
You can estimate single-output and multiple-output ARX models using the arx and iv4 commands. For information about the algorithms, see Polynomial Model Estimation Algorithms.
You can use the following general syntax to both configure and estimate ARX models:
% Using ARX method m = arx(data,[na nb nk],opt) % Using IV method m = iv4(data,[na nb nk],opt)
data is the estimation data and [na nb nk] specifies the model orders, as discussed in What Are Polynomial Models?.
The third input argument opt contains the options for configuring the estimation of the ARX model, such as handling of initial conditions and input offsets. You can create and configure the option set opt using the arxOptions and iv4Options commands. The three input arguments can also be followed by name and value pairs to specify optional model structure attributes such as InputDelay, ioDelay, and IntegrateNoise.
To get discrete-time models, use the time-domain data (iddata object).
For more information about validating you model, see Validating Models After Estimation.
You can use pem or polyest to refine parameter estimates of an existing polynomial model, as described in Refining Linear Parametric Models.
For detailed information about these commands, see the corresponding reference page.
Tip You can use the estimated ARX model for initializing a nonlinear estimation at the command line, which improves the fit of the model. See Using Linear Model for Nonlinear ARX Estimation. |
You can estimate any polynomial model using the iterative prediction-error estimation method polyest. For Gaussian disturbances of unknown variance, this method gives the maximum likelihood estimate. The resulting models are stored as idpoly model objects.
Use the following general syntax to both configure and estimate polynomial models:
m = polyest(data, [na nb nc nd nf nk], opt,Name,Value)
where data is the estimation data. na, nb, nc, nd, nf are integers that specify the model orders, and nk specifies the input delays for each input.For more information about model orders, see What Are Polynomial Models?.
If you want to estimate the coefficients of all five polynomials, A, B, C, D, and F, you must specify an integer order for each polynomial. However, if you want to specify an ARMAX model for example, which includes only the A, B, and C polynomials, you must set nd and nf to zero matrices of the appropriate size. For some simpler configurations, there are dedicated estimation commands such as arx, armax, bj, and oe, which deliver the required model by using just the required orders. For example, oe(data, [nb nf nk],opt) estimates an output-error structure polynomial model.
In addition to the polynomial models listed in What Are Polynomial Models?, you can use polyest to model the ARARX structure—called the generalized least-squares model—by setting nc=nf=0. You can also model the ARARMAX structure—called the extended matrix model—by setting nf=0.
The third input argument, opt, contains the options for configuring the estimation of the polynomial model, such as handling of initial conditions, input offsets and search algorithm. You can create and configure the option set opt using the polyestOptions command. The three input arguments can also be followed by name and value pairs to specify optional model structure attributes such as InputDelay, ioDelay, and IntegrateNoise.
For ARMAX, Box-Jenkins, and Output-Error models—which can only be estimated using the iterative prediction-error method—use the armax, bj, and oe estimation commands, respectively. These commands are versions of polyest with simplified syntax for these specific model structures, as follows:
m = armax(Data,[na nb nc nk]) m = oe(Data,[nb nf nk]) m = bj(Data,[nb nc nd nf nk])
Similar to polyest, you can specify as input arguments the option set configured using commands armaxOptions, oeOptions, and bjOptions for the estimators armax, oe, and bj respectively. You can also use name and value pairs to configure additional model structure attributes.
Tip If your data is sampled fast, it might help to apply a lowpass filter to the data before estimating the model, or specify a frequency range for the Focus property during estimation. For example, to model only data in the frequency range 0-10 rad/s, use the Focus property, as follows: opt = oeOptions('Focus',[0 10]) m = oe(Data, [nb nf nk], opt) |
For more information about validating your model, see Validating Models After Estimation.
You can use pem or polyest to refine parameter estimates of an existing polynomial model (of any configuration), as described in Refining Linear Parametric Models.
For a model with Ny (Ny > 1) outputs and Nu inputs, the polynomials A, B, C, D, and F are specified as cell arrays of row vectors. Each entry in the cell array contains the coefficients of a particular polynomial that relates input, output, and noise values. Orders are matrices of integers used as input arguments to the estimation commands.
Polynomial | Dimension | Relation Described | Orders |
---|---|---|---|
A | N_{y}-by-N_{y} array of row vectors | A{i,j} contains coefficients of relation between output y_{i} and output y_{j} | na: N_{y}-by-N_{y} matrix such that each entry contains the degree of the corresponding A polynomial. |
B | N_{y}-by-N_{u} array of row vectors | B{i,j} contain coefficients of relations between output y_{i} and input u_{j} | nk: N_{y}-by-N_{u} matrix such that each entry contains the number of leading fixed zeros of the corresponding B polynomial (input delay). nb: N_{y}-by-N_{u} matrix such nb(i,j) = length(B{i,j})- nk(i,j). |
C,D | N_{y}-by-1 array of row vectors | C{i} and D{i} contain coefficients of relations between output y_{i} and noise e_{i} | nc and nd are N_{y}-by-1 matrices such that each entry contains the degree of the corresponding C and D polynomial, respectively. |
F | N_{y}-by-N_{u} array of row vectors | F{i,j} contains coefficients of relations between output y_{i} and input u_{j} | nf: N_{y}-by-N_{u} matrix such that each entry contains the degree of the corresponding F polynomial. |
For more information, see idpoly.
For example, consider the ARMAX set of equations for a 2 output, 1 input model:
$$\begin{array}{l}{\text{y}}_{1}\text{(t)+0}{\text{.5y}}_{1}\text{(t-1)+0}{\text{.9y}}_{2}\text{(t-1)+0}{\text{.1y}}_{2}{\text{(t-2)=u(t)+5u(t-1)+2u(t-2)+e}}_{1}\text{(t)+0}{\text{.01e}}_{1}\text{(t-1)}\\ \text{}{\text{y}}_{2}\text{(t)+0}{\text{.05y}}_{2}\text{(t-1)+0}{\text{.3y}}_{2}{\text{(t-2)=10u(t-2)+e}}_{2}\text{(t)+0}{\text{.1e}}_{2}\text{(t-1)+0}{\text{.02e}}_{2}\text{(t-2)}\end{array}$$
y_{1} andy_{2} represent the two outputs and u represents the input variable. e_{1} and e_{2} represent the white noise disturbances on the outputs, y_{1} and y_{2}, respectively. To represent these equations as an ARMAX form polynomial using idpoly, configure the A, B, and C polynomials as follows:
A = cell(2,2); A{1,1} = [1 0.5]; A{1,2} = [0 0.9 0.1]; A{2,1} = [0]; A{2,2} = [1 0.05 0.3]; B = cell(2,1); B{1,1} = [1 5 2]; B{2,1} = [0 0 10]; C = cell(2,1); C{1} = [1 0.01]; C{2} = [1 0.1 0.02]; model = idpoly(A,B,C)
model = Discrete-time ARMAX model: Model for output number 1: A(z)y_1(t) = - A_i(z)y_i(t) + B(z)u(t) + C(z)e_1(t) A(z) = 1 + 0.5 z^-1 A_2(z) = 0.9 z^-1 + 0.1 z^-2 B(z) = 1 + 5 z^-1 + 2 z^-2 C(z) = 1 + 0.01 z^-1 Model for output number 2: A(z)y_2(t) = B(z)u(t) + C(z)e_2(t) A(z) = 1 + 0.05 z^-1 + 0.3 z^-2 B(z) = 10 z^-2 C(z) = 1 + 0.1 z^-1 + 0.02 z^-2 Sample time: unspecified Parameterization: Polynomial orders: na=[1 2;0 2] nb=[3;1] nc=[1;2] nk=[0;2] Number of free coefficients: 12 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Created by direct construction or transformation. Not estimated.
model is a discrete-time ARMAX model with unspecified sample-time. When estimating such models, you need to specify the orders of these polynomials as input arguments.
In the System Identification app. You can enter the matrices directly in the Orders field.
At the command line. Define variables that store the model order matrices and specify these variables in the model-estimation command.
You can specify how the estimation algorithm weighs the fit at various frequencies. This information supports the estimation procedures Estimate Polynomial Models in the App and Using polyest to Estimate Polynomial Models.
In the System Identification app. Set Focus to one of the following options:
Prediction — Uses the inverse of the noise model H to weigh the relative importance of how closely to fit the data in various frequency ranges. Corresponds to minimizing one-step-ahead prediction, which typically favors the fit over a short time interval. Optimized for output prediction applications.
Simulation — Uses the input spectrum to weigh the relative importance of the fit in a specific frequency range. Does not use the noise model to weigh the relative importance of how closely to fit the data in various frequency ranges. Optimized for output simulation applications.
Stability — Estimates the best stable model. For more information about model stability, see Unstable Models.
Filter — Specify a custom filter to open the Estimation Focus dialog box, where you can enter a filter, as described in Simple Passband Filter or Defining a Custom Filter. This prefiltering applies only for estimating the dynamics from input to output. The disturbance model is determined from the unfiltered estimation data.
At the command line. Specify the focus as an estimation option (created using polyestOptions, oeOptions etc.) using the same options as in the app. For example, use this command to estimate an ARX model and emphasize the frequency content related to the input spectrum only:
opt = arxOptions('Focus', 'simulation'); m = arx(data,[2 2 3],opt)
This Focus setting might produce more accurate simulation results, provided the orders picked are optimal for the given data..
When you use the pem or polyest to estimate ARMAX, Box-Jenkins (BJ), Output-Error (OE), you must specify how the algorithm treats initial conditions.
This information supports the estimation procedures Estimate Polynomial Models in the App and Using polyest to Estimate Polynomial Models.
In the System Identification app. For ARMAX, OE, and BJ models, set Initial state to one of the following options:
Auto — Automatically chooses Zero, Estimate, or Backcast based on the estimation data. If initial states have negligible effect on the prediction errors, the initial states are set to zero to optimize algorithm performance.
Zero — Sets all initial states to zero.
Estimate — Treats the initial states as an unknown vector of parameters and estimates these states from the data.
Backcast — Estimates initial states using a smoothing filter.
At the command line. Specify the initial conditions as an estimation option. Use polyestOptions to configure options for the polyest command, armaxOptions for the armax command etc. Set the InitialCondition option to the desired value in the option set. For example, use this command to estimate an ARMAX model and set the initial states to zero:
opt = armaxOptions('InitialCondition','zero') m = armax(data,[2 2 2 3],opt)
For a complete list of values for the InitialCondition estimation option, see the armaxOptions reference page.
For linear ARX and AR models, you can choose between the ARX and IV algorithms. ARX implements the least-squares estimation method that uses QR-factorization for overdetermined linear equations. IV is the instrument variable method. For more information about IV, see the section on variance-optimal instruments in System Identification: Theory for the User, Second Edition, by Lennart Ljung, Prentice Hall PTR, 1999.
The ARX and IV algorithms treat noise differently. ARX assumes white noise. However, the instrumental variable algorithm, IV, is not sensitive to noise color. Thus, use IV when the noise in your system is not completely white and it is incorrect to assume white noise. If the models you obtained using ARX are inaccurate, try using IV.
Note: AR models apply to time-series data, which has no input. For more information, see Time-Series Model Identification. For more information about working with AR and ARX models, see Identifying Input-Output Polynomial Models. |
This example shows how to estimate a linear, polynomial model with an ARMAX structure for a three-input and single-output (MISO) system using the iterative estimation method armax. For a summary of all available estimation commands in the toolbox, see Model Estimation Commands.
Load a sample data set z8 with three inputs and one output, measured at 1 -second intervals and containing 500 data samples.
load iddata8
Use armax to both construct the idpoly model object, and estimate the parameters:
Typically, you try different model orders and compare results, ultimately choosing the simplest model that best describes the system dynamics. The following command specifies the estimation data set, z8 , and the orders of the A , B , and C polynomials as na , nb , and nc, respectively. nk of [0 0 0] specifies that there is no input delay for all three input channels.
opt = armaxOptions;
opt.Focus = 'simulation';
opt.SearchOption.MaxIter = 50;
opt.SearchOption.Tolerance = 1e-5;
na = 4;
nb = [3 2 3];
nc = 4;
nk = [0 0 0];
m_armax = armax(z8, [na nb nc nk], opt);
Focus, Tolerance, and MaxIter are estimation options that configure the estimation objective function and the attributes of the search algorithm. The Focus option specifies whether the model is optimized for simulation or prediction applications. The Tolerance and MaxIter search options specify when to stop estimation. For more information about these properties, see the armaxOptions reference page.
armax is a version of polyest with simplified syntax for the ARMAX model structure. The armax method both constructs the idpoly model object and estimates its parameters.
View information about the resulting model object.
m_armax
m_armax = Discrete-time ARMAX model: A(z)y(t) = B(z)u(t) + C(z)e(t) A(z) = 1 - 1.284 z^-1 + 0.3048 z^-2 + 0.2648 z^-3 - 0.05708 z^-4 B1(z) = -0.07547 + 1.087 z^-1 + 0.7166 z^-2 B2(z) = 1.019 + 0.1142 z^-1 B3(z) = -0.06739 + 0.06828 z^-1 + 0.5509 z^-2 C(z) = 1 - 0.06096 z^-1 - 0.1296 z^-2 + 0.02489 z^-3 - 0.04699 z^-4 Sample time: 1 seconds Parameterization: Polynomial orders: na=4 nb=[3 2 3] nc=4 nk=[0 0 0] Number of free coefficients: 16 Use "polydata", "getpvec", "getcov" for parameters and their uncertainties. Status: Estimated using ARMAX on time domain data "z8". Fit to estimation data: 80.86% (simulation focus) FPE: 1.056, MSE: 0.9868
m_armax is an idpoly model object. The coefficients represent estimated parameters of this polynomial model. You can use present(m_armax) to show additional information about the model, including parameter uncertainties.
View all property values for this model.
get(m_armax)
a: [1 -1.2836 0.3048 0.2648 -0.0571] b: {[-0.0755 1.0870 0.7166] [1.0188 0.1142] [1x3 double]} c: [1 -0.0610 -0.1296 0.0249 -0.0470] d: 1 f: {[1] [1] [1]} IntegrateNoise: 0 Variable: 'z^-1' ioDelay: [0 0 0] Structure: [1x1 pmodel.polynomial] NoiseVariance: 0.9899 Report: [1x1 idresults.polyest] InputDelay: [3x1 double] OutputDelay: 0 Ts: 1 TimeUnit: 'seconds' InputName: {3x1 cell} InputUnit: {3x1 cell} InputGroup: [1x1 struct] OutputName: {'y1'} OutputUnit: {''} OutputGroup: [1x1 struct] Name: '' Notes: {} UserData: [] SamplingGrid: [1x1 struct]
The Report model property contains detailed information on the estimation results. To view the properties and values inside Report, use dot notation. For example:
m_armax.Report
Status: 'Estimated using POLYEST with Focus = "simulation"' Method: 'ARMAX' InitialCondition: 'zero' Fit: [1x1 struct] Parameters: [1x1 struct] OptionsUsed: [1x1 idoptions.polyest] RandState: [1x1 struct] DataUsed: [1x1 struct] Termination: [1x1 struct]
This action displays the contents of estimation report such as model quality measures (Fit), search termination criterion (Termination), and a record of estimation data (DataUsed) and options (OptionsUsed).