Estimating Black-Box Polynomial Models :: Tutorial - Identifying Linear Models Using the Command Line (System Identification Toolbox™)

Learn more about System Identification Toolbox

Estimating Black-Box Polynomial Models

On this page…
Model Orders for Estimating Polynomial Models Estimating a Linear ARX Model Estimating a State-Space Model Estimating a Box-Jenkins Model Comparing Model Output to Measured Output

On this page…

Model Orders for Estimating Polynomial Models

Estimating a Linear ARX Model

Estimating a State-Space Model

Estimating a Box-Jenkins Model

Comparing Model Output to Measured Output

Model Orders for Estimating Polynomial Models

In this portion of the tutorial, you estimate several different types of black-box, input-output polynomial models.

You must have already prepared your data, as described in Preparing Data.

You can use the following previous results of estimated model orders to specify the orders of the polynomial model:

For the first input/output combination, use:
- Two poles, corresponding to n_a=2 in the ARX model.
- One zero, corresponding to n_b=2 in the ARX model.
- Delay of 5, corresponding to n_k=5 samples (or 2.5 minutes) in the ARX model.
For the second input/output combination, use:
- One pole, corresponding to n_a=1 in the ARX model.
- No zeros, corresponding to n_b=1 in the ARX model.
- Delay of 10, corresponding to n_k=10 samples (or 5 minutes) in the ARX model.

Estimating a Linear ARX Model

About ARX Models

For a single-input/single-output system (SISO), the ARX model structure is:

y(t) represents the output at time t, u(t) represents the input at time t, n_a is the number of poles, n_b is the number of zeros plus 1, n_k is the number of samples before the input affects output of the system (called the delay or dead time of the model), and e(t) is the white-noise disturbance.

The ARX model structure does not distinguish between the poles for individual input/output paths: dividing the ARX equation by A, which contains the poles, shows that A appears in the denominator for both inputs. Therefore, you can set n_a to the sum of the poles for each input/output pair, which is equal to 3 in this case.

The System Identification Toolbox product estimates the parameters and using the data and the model orders you specify.

Estimating ARX Models Using arx

Use arx to compute the polynomial coefficients using the fast, noniterative method arx:

marx = arx(Ze1,'na',3,'nb',[2 1],'nk',[5 10]);
present(marx) % Displays model parameters
              % with uncertainty information

MATLAB estimates the polynomials A, B1, and B2:

Discrete-time IDPOLY model: A(q)y(t) = B(q)u(t) + e(t)
A(q) = 1 - 1.027 (+-0.02917) q^-1 
         + 0.1675 (+-0.04214) q^-2 
         + 0.01307 (+-0.02591) q^-3
B1(q) = 1.86 (+-0.1896) q^-5 - 1.608 (+-0.1894) q^-6
B2(q) = 1.612 (+-0.07417) q^-10

The uncertainty for each of the model parameters is computed to 1 standard deviation and appears in parentheses next to each parameter value.

Tip Alternatively, you can use the following shorthand syntax and specify model orders as a single vector:

marx = arx(Ze1,[3 2 1 5 10])

Accessing Model Data

The model you estimated, marx, is a discrete-time idpoly object. To get the properties of this model object, you can use the get function:

get(marx) 
         a: [1 -1.0266 0.1675 0.0131]
                   b: [2x11 double]
                   c: 1
                   d: 1
                   f: [2x1 double]
                  da: [0 0.0292 0.0421 0.0259]
                  db: [2x11 double]
                  dc: 0
                  dd: 0
                  df: [2x1 double]
                  na: 3
                  nb: [2 1]
                  nc: 0
                  nd: 0
                  nf: [0 0]
                  nk: [5 10]
        InitialState: 'Auto'
                Name: ''
                  Ts: 0.5000
           InputName: {2x1 cell}
           InputUnit: {2x1 cell}
          OutputName: {'ProductionRate'}
          OutputUnit: {'mg/min'}
            TimeUnit: 'min'
     ParameterVector: [6x1 double]
               PName: {}
    CovarianceMatrix: [6x6 double]
       NoiseVariance: 2.7611
          InputDelay: [2x1 double]
           Algorithm: [1x1 struct]
      EstimationInfo: [1x1 struct]
               Notes: {}
            UserData: []

You can access the information stored by these properties using dot notation. For example, you can compute the discrete poles of the model by finding the roots of the A polynomial:

marx_poles=roots(marx.a)

marx_poles =

    0.7953
    0.2883
   -0.0570

In this case, you access the A polynomial using marx.a.

The model marx describes system dynamics using three discrete poles.

Tip You can also use the zpkdata command to compute the poles of a model directly.

Learn More

To learn more about estimating polynomial models, see the corresponding sections in the System Identification Toolbox User's Guide.

For more information about accessing model data, see the topic on extracting numerical data from linear models in the System Identification Toolbox User's Guide.

Estimating a State-Space Model

About State-Space Models

The general state-space model structure is:

y(t) represents the output at time t, u(t) represents the input at time t, x(t) is the state vector at time t, and e(t) is the white-noise disturbance.

You must specify a single integer as the model order (dimension of the state vector) to estimate a state-space model. By default, the delay equals 1.

The System Identification Toolbox product estimates the state-space matrices A, B, C, D, and K using the model order and the data you specify.

The state-space model structure is a good choice for quick estimation because it contains only two parameters: n is the number of poles (the size of the A matrix) and nk is the delay.

Estimating State-Space Models Using n4sid

Use the n4sid command to specify a range of model orders and evaluate the performance of several state-space models (orders 2 to 8):

mn4sid = n4sid(Ze1,2:8,'nk',[5 10]);

This command uses the fast, noniterative (subspace) method and opens the following plot. You use this plot to decide which states provide a significant relative contribution to the input/output behavior, and which states provide the smallest contribution.

The vertical axis is a relative measure of how much each state contributes to the input/output behavior of the model (log of singular values of the covariance matrix). The horizontal axis corresponds to the model order n. This plot recommends n=3, indicated by a red rectangle.

To select this model order, type 3 in the MATLAB Command Window, and press Enter.

By default, n4sid uses a free parameterization of the state-space form. To estimate a canonical form instead, set the value of the SSParameterization property to 'Canonical':

mCanonical = n4sid(Ze1,3,'nk',[5 10],...
             'ssparameterization','canonical');
present(mCanonical)    % Displays model properties

Note When you examine the displayed properties, notice that the model order is high. This high order occurs because the model uses additional states to incorporate input delays.

Learn More

To learn more about estimating state-space models, see the corresponding topic in the System Identification Toolbox User's Guide.

Estimating a Box-Jenkins Model

About Box-Jenkins Models

The general Box-Jenkins (BJ) structure is:

To estimate a BJ model, you need to specify the parameters n_b, n_f, n_c, n_d, and n_k.

Whereas the ARX model structure does not distinguish between the poles for individual input/output paths, the BJ model provides more flexibility in modeling the poles and zeros of the disturbance separately from the poles and zeros of the system dynamics.

Estimating a BJ Model Using pem

You can use pem to estimate the BJ model. pem is an iterative method and has the following general syntax:

pem(data,'na',na,'nb',nb,'nc',nc,'nd',nd,'nf',nf,'nk',nk)

In this case, data is an iddata or idfrd object, and na, nb, nc, nd, nf, and nk specify the model order.

To estimate the BJ model, type the following command in the MATLAB Command Window:

mbj = pem(Ze1,'nf',[2 1],'nb',[2 1],'nc',1,'nd',1,'nk',[5 10]);
present(mbj)

This command specifies nf=2, nb=2, nk=5 for the first input, and nf=nb=1 and nk=10 for the second input.

Tip Alternatively, you can use the following shorthand syntax that specifies the orders as a single vector:

mbj = bj(Ze1,[2 1 1 1 2 1 5 10]);

bj is a version of pem that specifically estimates the BJ model structure.

MATLAB estimates the following polynomial coefficients:

Discrete-time IDPOLY model: 
y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)
B1(q) = 1.823(+-0.1792)q^-5 - 1.315(+-0.2368)q^-6          
B2(q) = 1.791(+-0.06435)q^-10                              
C(q) = 1 + 0.1066(+-0.04015)q^-1                           
D(q) = 1 - 0.7453(+-0.027)q^-1                             
F1(q) = 1 - 1.321(+-0.06939)q^-1 + 0.5911(+-0.05516)q^-2   
F2(q) = 1 - 0.8314(+-0.006445)q^-1

The uncertainty for each of the model parameters is computed to 1 standard deviation and appears in parentheses next to each parameter value.

The polynomials C and D give the numerator and the denominator of the noise model, respectively.

Learn More

To learn more about estimating BJ models, see the corresponding topic in the System Identification Toolbox User's Guide.

Comparing Model Output to Measured Output

To compare the output of the ARX, state-space, and Box-Jenkins models to the measured output, use the compare function:

compare(Zv1,marx,mn4sid,mbj)

compare plots the measured output in the validation data set against the simulated output from the models. The input data from the validation data set serves as input to the models.

Measured Output and Simulated Outputs

To perform residual analysis on the ARX model, type the following command:

resid(Zv1,marx)

Because the sample system has two inputs, there are two plots showing the cross-correlation of the residuals with each input. Press Enter to view the cross-correlation with the second input.

To perform residual analysis on the state-space model, type the following command:

resid(Zv1,mn4sid)

Finally, to perform residual analysis on the BJ model, type the following command:

resid(Zv1,mbj)

All three models simulate the output equally well and have uncorrelated residuals. Therefore, choose the ARX model because it is the simplest of the three input-output polynomial models and adequately captures the process dynamics.

Estimating Continuous-Time Transfer Functions (Process Models)

Simulating and Predicting Model Output

Recommended Products

Get MATLAB trial software

Includes the most popular MATLAB recorded presentations with Q&A sessions led by MATLAB experts.

Get the Interactive Kit