Identify Linear Models Using the Command Line

Introduction

Objectives

Estimate and validate linear models from multiple-input/single-output (MISO) data to find the one that best describes the system dynamics.

After completing this tutorial, you will be able to accomplish the following tasks using the command line:

  • Create data objects to represent data.

  • Plot the data.

  • Process data by removing offsets from the input and output signals.

  • Estimate and validate linear models from the data.

  • Simulate and predict model output.

    Note:   This tutorial uses time-domain data to demonstrate how you can estimate linear models. The same workflow applies to fitting frequency-domain data.

Data Description

This tutorial uses the data file co2data.mat, which contains two experiments of two-input and single-output (MISO) time-domain data from a steady-state that the operator perturbed from equilibrium values.

In the first experiment, the operator introduced a pulse wave to both inputs. In the second experiment, the operator introduced a pulse wave to the first input and a step signal to the second input.

Preparing Data

Loading Data into the MATLAB Workspace

Load the data.

load co2data;

This command loads the following five variables into the MATLAB Workspace:

  • Input_exp1 and Output_exp1 are the input and output data from the first experiment, respectively.

  • Input_exp2 and Output_exp2 are the input and output data from the second experiment, respectively.

  • Time is the time vector from 0 to 1000 minutes, increasing in equal increments of 0.5 min.

For both experiments, the input data consists of two columns of values. The first column of values is the rate of chemical consumption (in kilograms per minute), and the second column of values is the current (in amperes). The output data is a single column of the rate of carbon-dioxide production (in milligrams per minute).

Plotting the Input/Output Data

Plot the input and output data from both experiments.

plot(Time,Input_exp1,Time,Output_exp1)
legend('Input 1','Input 2','Output 1')
figure
plot(Time,Input_exp2,Time,Output_exp2)
legend('Input 1','Input 2','Output 1')

The first plot shows the first experiment, where the operator applies a pulse wave to each input to perturb it from its steady-state equilibrium.

The second plot shows the second experiment, where the operator applies a pulse wave to the first input and a step signal to the second input.

Removing Equilibrium Values from the Data

Plotting the data, as described in Plotting the Input/Output Data, shows that the inputs and the outputs have nonzero equilibrium values. In this portion of the tutorial, you subtract equilibrium values from the data.

The reason you subtract the mean values from each signal is because, typically, you build linear models that describe the responses for deviations from a physical equilibrium. With steady-state data, it is reasonable to assume that the mean levels of the signals correspond to such an equilibrium. Thus, you can seek models around zero without modeling the absolute equilibrium levels in physical units.

Zoom in on the plots to see that the earliest moment when the operator applies a disturbance to the inputs occurs after 25 minutes of steady-state conditions (or after the first 50 samples). Thus, the average value of the first 50 samples represents the equilibrium conditions.

Use the following commands to remove the equilibrium values from inputs and outputs in both experiments:

Input_exp1 = Input_exp1-...
   ones(size(Input_exp1,1),1)*mean(Input_exp1(1:50,:));
Output_exp1 = Output_exp1-...
   mean(Output_exp1(1:50,:));
Input_exp2 = Input_exp2-...
   ones(size(Input_exp2,1),1)*mean(Input_exp2(1:50,:));
Output_exp2 = Output_exp2-...
   mean(Output_exp2(1:50,:));

Using Objects to Represent Data for System Identification

The System Identification Toolbox™ data objects, iddata and idfrd, encapsulate data values and data properties into a single entity. You can use the System Identification Toolbox commands to conveniently manipulate these data objects as single entities.

In this portion of the tutorial, you create two iddata objects, one for each of the two experiments. You use the data from Experiment 1 for model estimation, and the data from Experiment 2 for model validation. You work with two independent data sets because you use one data set for model estimation and the other for model validation.

    Note:   When two independent data sets are not available, you can split one data set into two parts, assuming that each part contains enough information to adequately represent the system dynamics.

Creating iddata Objects

You must have already loaded the sample data into the MATLAB® workspace, as described in Loading Data into the MATLAB Workspace.

Use these commands to create two data objects, ze and zv :

Ts = 0.5; % Sampling interval is 0.5 min
ze = iddata(Output_exp1,Input_exp1,Ts);
zv = iddata(Output_exp2,Input_exp2,Ts);

ze contains data from Experiment 1 and zv contains data from Experiment 2. Ts is the sampling interval.

The iddata constructor requires three arguments for time-domain data and has the following syntax:

data_obj = iddata(output,input,sampling_interval);

To view the properties of an iddata object, use the get command. For example, type this command to get the properties of the estimation data:

get(ze)
ans = 

              Domain: 'Time'
                Name: ''
          OutputData: [2001x1 double]
                   y: 'Same as OutputData'
          OutputName: {'y1'}
          OutputUnit: {''}
           InputData: [2001x2 double]
                   u: 'Same as InputData'
           InputName: {2x1 cell}
           InputUnit: {2x1 cell}
              Period: [2x1 double]
         InterSample: {2x1 cell}
                  Ts: 0.5000
              Tstart: []
    SamplingInstants: [2001x0 double]
            TimeUnit: 'seconds'
      ExperimentName: 'Exp1'
               Notes: {}
            UserData: []

To learn more about the properties of this data object, see the iddata reference page.

To modify data properties, you can use dot notation or the set command. For example, to assign channel names and units that label plot axes, type the following syntax in the MATLAB Command Window:

% Set time units to minutes
ze.TimeUnit = 'min';
% Set names of input channels
ze.InputName = {'ConsumptionRate','Current'};
% Set units for input variables
ze.InputUnit = {'kg/min','A'};
% Set name of output channel
ze.OutputName = 'Production';
% Set unit of output channel
ze.OutputUnit = 'mg/min';

% Set validation data properties
zv.TimeUnit = 'min';
zv.InputName = {'ConsumptionRate','Current'};
zv.InputUnit = {'kg/min','A'};
zv.OutputName = 'Production';
zv.OutputUnit = 'mg/min';

You can verify that the InputName property of ze is changed, or "index" into this property:

ze.inputname;

    Tip   Property names, such as InputUnit, are not case sensitive. You can also abbreviate property names that start with Input or Output by substituting u for Input and y for Output in the property name. For example, OutputUnit is equivalent to yunit.

Plotting the Data in a Data Object

You can plot iddata objects using the plot command.

Plot the estimation data.

plot(ze)

The bottom axes show inputs ConsumptionRate and Current, and the top axes show the output ProductionRate .

Plot the validation data in a new figure window.

figure   % Open a new MATLAB Figure window
plot(zv) % Plot the validation data

Zoom in on the plots to see that the experiment process amplifies the first input ( ConsumptionRate ) by a factor of 2, and amplifies the second input (Current ) by a factor of 10.

Selecting a Subset of the Data

Before you begin, create a new data set that contains only the first 1000 samples of the original estimation and validation data sets to speed up the calculations.

Ze1 = ze(1:1000);
Zv1 = zv(1:1000);

For more information about indexing into iddata objects, see the corresponding reference page.

Estimating Impulse Response Models

Why Estimate Step- and Frequency-Response Models?

Frequency-response and step-response are nonparametric models that can help you understand the dynamic characteristics of your system. These models are not represented by a compact mathematical formula with adjustable parameters. Instead, they consist of data tables.

In this portion of the tutorial, you estimate these models using the data set ze. You must have already created ze, as described in Creating iddata Objects.

The response plots from these models show the following characteristics of the system:

  • The response from the first input to the output might be a second-order function.

  • The response from the second input to the output might be a first-order or an overdamped function.

Estimating the Frequency Response

The System Identification Toolbox product provides three functions for estimating the frequency response:

  • etfe computes the empirical transfer function using Fourier analysis.

  • spa estimates the transfer function using spectral analysis for a fixed frequency resolution.

  • spafdr lets you specify a variable frequency resolution for estimating the frequency response.

Use spa to estimate the frequency response.

Ge=spa(ze);

Plot the frequency response as a Bode plot.

bode(Ge)

The amplitude peaks at the frequency of about 0.7 rad/s, which suggests a possible resonant behavior (complex poles) for the first input-to-output combination - ConsumptionRate to Production .

In both plots, the phase rolls off rapidly, which suggests a time delay for both input/output combinations.

Estimating the Empirical Step Response

To estimate the step response from the data, first estimate a non-parametric impulse response model (FIR filter) from data and then plot its step response.

% model estimation
Mimp = impulseest(Ze1,60);

% step response
step(Mimp)

The step response for the first input/output combination suggests an overshoot, which indicates the presence of an underdamped mode (complex poles) in the physical system.

The step response from the second input to the output shows no overshoot, which indicates either a first-order response or a higher-order response with real poles (overdamped response).

The step-response plot indicates a nonzero delay in the system, which is consistent with the rapid phase roll-off you got in the Bode plot you created in Estimating the Empirical Step Response.

Estimating Delays in the Multiple-Input System

Why Estimate Delays?

To identify parametric black-box models, you must specify the input/output delay as part of the model order.

If you do not know the input/output delays for your system from the experiment, you can use the System Identification Toolbox software to estimate the delay.

Estimating Delays Using the ARX Model Structure

In the case of single-input systems, you can read the delay on the impulse-response plot. However, in the case of multiple-input systems, such as the one in this tutorial, you might be unable to tell which input caused the initial change in the output and you can use the delayest command instead.

The command estimates the time delay in a dynamic system by estimating a low-order, discrete-time ARX model with a range of delays, and then choosing the delay that corresponding to the best fit.

The ARX model structure is one of the simplest black-box parametric structures. In discrete-time, the ARX structure is a difference equation with the following form:

y(t)+a1y(t1)++anay(tna)=         b1u(tnk)++bnbu(tnknb+1)+e(t)

y(t) represents the output at time t, u(t) represents the input at time t, na is the number of poles, nb is the number of b parameters (equal to the number of zeros plus 1), nk 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.

delayest assumes that na=nb=2 and that the noise e is white or insignificant, and estimates nk.

To estimate the delay in this system, type:

delayest(ze)
ans =

     5    10

This result includes two numbers because there are two inputs: the estimated delay for the first input is 5 data samples, and the estimated delay for the second input is 10 data samples. Because the sampling interval for the experiments is 0.5 min, this corresponds to a 2.5 -min delay before the first input affects the output, and a 5.0 -min delay before the second input affects the output.

Estimating Delays Using Alternative Methods

There are two alternative methods for estimating the time delay in the system:

  • Plot the time plot of the input and output data and read the time difference between the first change in the input and the first change in the output. This method is practical only for single-input/single-output system; in the case of multiple-input systems, you might be unable to tell which input caused the initial change in the output.

  • Plot the impulse response of the data with a 1-standard-deviation confidence region. You can estimate the time delay using the time when the impulse response is first outside the confidence region.

Estimating Model Orders Using an ARX Model Structure

Why Estimate Model Order?

Model order is one or more integers that define the complexity of the model. In general, model order is related to the number of poles, the number of zeros, and the response delay (time in terms of the number of samples before the output responds to the input). The specific meaning of model order depends on the model structure.

To compute parametric black-box models, you must provide the model order as an input. If you do not know the order of your system, you can estimate it.

After completing the steps in this section, you get the following results:

  • For the first input/output combination: na=2, nb=2, and the delay nk=5.

  • For the second input/output combination: na=1, nb=1, and the delay nk=10.

Later, you explore different model structures by specifying model-order values that are slight variations around these initial estimate.

Commands for Estimating the Model Order

In this portion of the tutorial, you use struc, arxstruc, and selstruc to estimate and compare simple polynomial (ARX) models for a range of model orders and delays, and select the best orders based on the quality of the model.

The following list describes the results of using each command:

  • struc creates a matrix of model-order combinations for a specified range of na, nb, and nk values.

  • arxstruc takes the output from struc, systematically estimates an ARX model for each model order, and compares the model output to the measured output. arxstruc returns the loss function for each model, which is the normalized sum of squared prediction errors.

  • selstruc takes the output from arxstruc and opens the ARX Model Structure Selection window, which resembles the following figure, to help you choose the model order.

    You use this plot to select the best-fit model.

    • The horizontal axis is the total number of parameters — na + nb.

    • 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.

    • nk 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 this tutorial, the Unexplained output variance (in %) value remains approximately constant for the combined number of parameters from 4 to 20. Such constancy indicates that model performance does not improve at higher orders. Thus, low-order models might fit the data equally well.

      Note:   When you use the same data set for estimation and validation, use the MDL and AIC criteria to select model orders. These criteria compensate for overfitting that results from using too many parameters. For more information about these criteria, see the selstruc reference page.

Model Order for the First Input-Output Combination

In this tutorial, there are two inputs to the system and one output and you estimate model orders for each input/output combination independently. You can either estimate the delays from the two inputs simultaneously or one input at a time.

It makes sense to try the following order combinations for the first input/output combination:

  • na=2:5

  • nb=1:5

  • nk=5

This is because the nonparametric models you estimated in Estimating Impulse Response Models show that the response for the first input/output combination might have a second-order response. Similarly, in Estimating Delays in the Multiple-Input System, the delay for this input/output combination was estimated to be 5.

To estimate model order for the first input/output combination:

  1. Use struc to create a matrix of possible model orders.

    NN1 = struc(2:5,1:5,5);
  2. Use selstruc to compute the loss functions for the ARX models with the orders in NN1.

    selstruc(arxstruc(ze(:,:,1),zv(:,:,1),NN1))

      Note:   (ze(:,:,1) selects the first input in the data.

    This command opens the interactive ARX Model Structure Selection window.

      Note:   The Rissanen MDL and Akaike AIC criteria produces equivalent results and are both indicated by a blue rectangle on the plot.

    The red rectangle represents the best overall fit, which occurs for na=2, nb=3, and nk=5. The height difference between the red and blue rectangles is insignificant. Therefore, you can choose the parameter combination that corresponds to the lowest model order and the simplest model.

  3. Click the blue rectangle, and then click Select to choose that combination of orders:

    na=2

    nb=2

    nk=5

  4. To continue, press any key while in the MATLAB Command Window.

Model Order for the Second Input-Output Combination

It makes sense to try the following order combinations for the second input/output combination:

  • na=1:3

  • nb=1:3

  • nk=10

This is because the nonparametric models you estimated in Estimating Impulse Response Models show that the response for the second input/output combination might have a first-order response. Similarly, in Estimating Delays in the Multiple-Input System, the delay for this input/output combination was estimated to be 10.

To estimate model order for the second input/output combination:

  1. Use struc to create a matrix of possible model orders.

    NN2 = struc(1:3,1:3,10);
  2. Use selstruc to compute the loss functions for the ARX models with the orders in NN2.

    selstruc(arxstruc(ze(:,:,2),zv(:,:,2),NN2))

    This command opens the interactive ARX Model Structure Selection window.

      Note:   The Akaike AIC and the overall best fit criteria produces equivalent results. Both are indicated by the same red rectangle on the plot.

    The height difference between all of the rectangles is insignificant and all of these model orders result in similar model performance. Therefore, you can choose the parameter combination that corresponds to the lowest model order and the simplest model.

  3. Click the yellow rectangle on the far left, and then click Select to choose the lowest order: na=1, nb=1, and nk=10.

  4. To continue, press any key while in the MATLAB Command Window.

Estimating Transfer Functions

Specifying the Structure of the Transfer Function

In this portion of the tutorial, you estimate a continuous-time transfer function. You must have already prepared your data, as described in Preparing Data. You can use the following results of estimated model orders to specify the orders of the model:

  • For the first input/output combination, use:

    • Two poles, corresponding to na=2 in the ARX model.

    • Delay of 5, corresponding to nk=5 samples (or 2.5 minutes) in the ARX model.

  • For the second input/output combination, use:

    • One pole, corresponding to na=1 in the ARX model

    • Delay of 10, corresponding to nk=10 samples (or 5 minutes) in the ARX model.

You can estimate a transfer function of these orders using the tfest command. For the estimation, you can also choose to view a progress report by setting the Display option to on in the option set created by the tfestOptions command.

Opt = tfestOptions('Display', 'on');

Collect the model orders and delays into variables to pass to tfest.

np = [2 1];
ioDelay = [2.5 5];

Estimate the transfer function.

mtf = tfest(Ze1, np, [], ioDelay, Opt);

View the model's coefficients.

mtf
mtf =
 
  From input "ConsumptionRate" to output "Production":
                   91.75 s - 1.272
  exp(-2.5*s) * ---------------------
                s^2 + 37.76 s + 1.009
 
  From input "Current" to output "Production":
                5.186
  exp(-5*s) * ----------
              s + 0.5066
 
Continuous-time identified transfer function.

Parameterization:
   Number of poles: [2 1]   Number of zeros: [1 0]
   Number of free coefficients: 6
   Use "tfdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                          
Estimated using TFEST on time domain data "Ze1". 
Fit to estimation data: 85.89% (simulation focus)
FPE: 6.424, MSE: 6.259                           

The model's display shows more than 85% fit to estimation data.

Validating the Model

In this portion of the tutorial, you create a plot that compares the actual output and the model output using the compare command.

compare(Zv1, mtf)

The comparison shows about 60% fit.

Residual Analysis

Use the resid command to evaluate the characteristics of the residuals.

resid(Zv1, mtf)

The residuals show high degree of auto-correlation. This is not unexpected since the model mtf does not have any components to describe the nature of the noise separately. To model both the measured input-output dynamics and the disturbance dynamics you will need to use a model structure that contains elements to describe the noise component. You can use bj, ssest and procest commands, which create models of polynomial, state-space and process structures respectively. These structures, among others, contain elements to capture the noise behavior.

Estimating Process Models

Specifying the Structure of the Process Model

In this portion of the tutorial, you estimate a low-order, continuous-time transfer function (process model). the System Identification Toolbox product supports continuous-time models with at most three poles (which might contain underdamped poles), one zero, a delay element, and an integrator.

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

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

  • For the first input/output combination, use:

    • Two poles, corresponding to na=2 in the ARX model.

    • Delay of 5, corresponding to nk=5 samples (or 2.5 minutes) in the ARX model.

  • For the second input/output combination, use:

    • One pole, corresponding to na=1 in the ARX model.

    • Delay of 10, corresponding to nk=10 samples (or 5 minutes) in the ARX model.

    Note:   Because there is no relationship between the number of zeros estimated by the discrete-time ARX model and its continuous-time counterpart, you do not have an estimate for the number of zeros. In this tutorial, you can specify one zero for the first input/output combination, and no zero for the second-output combination.

Use the idproc command to create two model structures, one for each input/output combination:

midproc0 = idproc({'P2ZUD','P1D'}, 'TimeUnit', 'minutes');

idproc accepts a cell array that contains two strings that specify the model structure for each input/output combination:

  • 'P2ZUD' represents a transfer function with two poles ( P2 ), one zero ( Z ), underdamped (complex-conjugate) poles ( U ) and a delay ( D ).

  • 'P1D' represents a transfer function with one pole ( P1 ) and a delay ( D ).

The example also uses the TimeUnit parameter to specify the unit of time used.

Viewing the Model Structure and Parameter Values

View the two resulting models.

midproc0
midproc0 =
Process model with 2 inputs: y = G11(s)u1 + G12(s)u2
  From input 1 to output 1:                         
                       1+Tz*s                       
  G11(s) = Kp * ---------------------- * exp(-Td*s) 
                1+2*Zeta*Tw*s+(Tw*s)^2              
                                                    
         Kp = NaN                                   
         Tw = NaN                                   
       Zeta = NaN                                   
         Td = NaN                                   
         Tz = NaN                                   
                                                    
  From input 2 to output 1:                         
             Kp                                     
  G12(s) = ---------- * exp(-Td*s)                  
            1+Tp1*s                                 
                                                    
        Kp = NaN                                    
       Tp1 = NaN                                    
        Td = NaN                                    
                                                    
Parameterization:
    'P2DUZ'    'P1D'
   Number of free coefficients: 8
   Use "getpvec", "getcov" for parameters and their uncertainties.

Status:                                                         
Created by direct construction or transformation. Not estimated.

The parameter values are set to NaN because they are not yet estimated.

Specifying Initial Guesses for Time Delays

Set the time delay property of the model object to 2.5 min and 5 min for each input/output pair as initial guesses. Also, set an upper limit on the delay because good initial guesses are available.

midproc0.Structure(1,1).Td.Value=2.5;
midproc0.Structure(1,2).Td.Value=5;
midproc0.Structure(1,1).Td.Maximum=3;
midproc0.Structure(1,2).Td.Maximum=7;

    Note:   When setting the delay constraints, you must specify the delays in terms of actual time units (minutes, in this case) and not the number of samples.

Estimating Model Parameters Using procest

procest is an iterative estimator of process models, which means that it uses an iterative nonlinear least-squares algorithm to minimize a cost function. The cost function is the weighted sum of the squares of the errors.

Depending on its arguments, procest estimates different black-box polynomial models. You can use procest, for example, to estimate parameters for linear continuous-time transfer-function, state-space, or polynomial model structures. To use procest, you must provide a model structure with unknown parameters and the estimation data as input arguments.

For this portion of the tutorial, you must have already defined the model structure, as described in Specifying the Structure of the Process Model. Use midproc0 as the model structure and Ze1 as the estimation data:

midproc = procest(Ze1, midproc0);
present(midproc)
                                                                  
midproc =                                                         
Process model with 2 inputs: y = G11(s)u1 + G12(s)u2              
  From input "ConsumptionRate" to output "Production":            
                       1+Tz*s                                     
  G11(s) = Kp * ---------------------- * exp(-Td*s)               
                1+2*Zeta*Tw*s+(Tw*s)^2                            
                                                                  
         Kp = -1.1804 +/- 0.29989                                 
         Tw = 1.6566 +/- 629.69                                   
       Zeta = 15.917 +/- 6039.2                                   
         Td = 2.425 +/- 56.972                                    
         Tz = -109.2 +/- 57.374                                   
                                                                  
  From input "Current" to output "Production":                    
             Kp                                                   
  G12(s) = ---------- * exp(-Td*s)                                
            1+Tp1*s                                               
                                                                  
        Kp = 10.264 +/- 0.048403                                  
       Tp1 = 2.0488 +/- 0.054899                                  
        Td = 4.9175 +/- 0.034373                                  
                                                                  
Parameterization:                                                 
    'P2DUZ'    'P1D'                                              
   Number of free coefficients: 8                                 
   Use "getpvec", "getcov" for parameters and their uncertainties.
                                                                  
Status:                                                           
Termination condition: Maximum number of iterations reached.      
Number of iterations: 20, Number of function evaluations: 276     
                                                                  
Estimated using PROCEST on time domain data "Ze1".                
Fit to estimation data: 86.2%                                     
FPE: 6.08, MSE: 5.984                                             
More information in model's "Report" property.                    

Unlike discrete-time polynomial models, continuous-time models let you estimate the delays. In this case, the estimated delay values are different from the initial guesses you specified of 2.5 and 5 , respectively. The large uncertainties in the estimated values of the parameters of G_1(s) indicate that the dynamics from input 1 to the output are not captured well.

To learn more about estimating models, see Process Models.

Validating the Model

In this portion of the tutorial, you create a plot that compares the actual output and the model output.

compare(Zv1,midproc)

The plot shows that the model output reasonably agrees with the measured output: there is an agreement of 60% between the model and the validation data.

Use resid to perform residual analysis:

resid(Zv1,midproc)

Because the sample system has two inputs, there are two cross-correlation plots of the residuals with each input, as shown in the following figure.

Autocorrelation and Cross-Correlations of Residuals with the First Input

The cross-correlation between the second input and residual errors is significant.

After MATLAB displays the first plot, press Enter to view the cross-correlation with the second input, as shown in the following figure.

Cross-Correlations of Residuals with the Second Input

In the preceding figure, the autocorrelation plot shows values outside the confidence region and indicates that the residuals are correlated.

Change the algorithm for iterative parameter estimation to Levenberg-Marquardt.

Opt = procestOptions;
Opt.SearchMethod = 'lm';
midproc1 = procest(Ze1, midproc, Opt);

Tweaking the algorithm properties or specifying initial parameter guesses and rerunning the estimation may improve the simulation results. Adding a noise model may improve prediction results but not necessarily the simulation results.

Estimating a Process Model with Noise Model

This portion of the tutorial shows how to estimate a process model and include a noise model in the estimation. Including a noise model typically improves model prediction results but not necessarily the simulation results.

Use the following commands to specify a first-order ARMA noise:

Opt = procestOptions;
Opt.DisturbanceModel = 'ARMA1';
midproc2 = procest(Ze1, midproc0, Opt);

You can type 'dist' instead of 'DisturbanceModel'. Property names are not case sensitive, and you only need to include the portion of the name that uniquely identifies the property.

Compare the resulting model to the measured data.

compare(Zv1,midproc2)

The plot shows that the model output maintains reasonable agreement with the validation-data output.

Perform residual analysis.

figure
resid(Zv1,midproc2)

Press Enter to view the cross-correlation of the residuals with the second input.

The next plot shows that adding a noise model produces uncorrelated residuals: the top set of axes show that the autocorrelation values are inside the confidence bounds. This indicates a more accurate model.

Estimating Black-Box Polynomial Models

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 na=2 in the ARX model.

    • One zero, corresponding to nb=2 in the ARX model.

    • Delay of 5, corresponding to nk=5 samples (or 2.5 minutes) in the ARX model.

  • For the second input/output combination, use:

    • One pole, corresponding to na=1 in the ARX model.

    • No zeros, corresponding to nb=1 in the ARX model.

    • Delay of 10, corresponding to nk=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)+a1y(t1)++anay(tna)=         b1u(tnk)++bnbu(tnknb+1)+e(t)

y(t) represents the output at time t, u(t) represents the input at time t, na is the number of poles, nb is the number of zeros plus 1, nk 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 na 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 a1an and b1bn 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
                                                                              
marx =                                                                        
Discrete-time ARX model:  A(z)y(t) = B(z)u(t) + e(t)                          
                                                                              
  A(z) = 1 - 1.027 (+/- 0.02917) z^-1 + 0.1675 (+/- 0.04214) z^-2             
                                              + 0.01307 (+/- 0.02591) z^-3    
                                                                              
  B1(z) = 1.86 (+/- 0.1896) z^-5 - 1.608 (+/- 0.1894) z^-6                    
                                                                              
  B2(z) = 1.612 (+/- 0.07417) z^-10                                           
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   na=3   nb=[2 1]   nk=[5 10]                           
   Number of free coefficients: 6                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Estimated using ARX on time domain data "Ze1".                                
Fit to estimation data: 90.7% (prediction focus)                              
FPE: 2.75, MSE: 2.719                                                         
More information in model's "Report" property.                                

MATLAB estimates the polynomials A , B1 , and B2. The uncertainty for each of the model parameters is computed to 1 standard deviation and appears in parentheses next to each parameter value.

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: {[0 0 0 0 0 1.8600 -1.6085]  [0 0 0 0 0 0 0 0 0 0 1.6117]}
                 c: 1
                 d: 1
                 f: {[1]  [1]}
    IntegrateNoise: 0
          Variable: 'z^-1'
           ioDelay: [0 0]
         Structure: [1x1 pmodel.polynomial]
     NoiseVariance: 2.7611
            Report: [1x1 idresults.arx]
        InputDelay: [2x1 double]
       OutputDelay: 0
                Ts: 0.5000
          TimeUnit: 'minutes'
         InputName: {2x1 cell}
         InputUnit: {2x1 cell}
        InputGroup: [1x1 struct]
        OutputName: {'Production'}
        OutputUnit: {'mg/min'}
       OutputGroup: [1x1 struct]
              Name: ''
             Notes: {}
          UserData: []
      SamplingGrid: [1x1 struct]

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 pole to compute the poles of a model directly.

Learn More.  To learn more about estimating polynomial models, see Input-Output Polynomial Models.

For more information about accessing model data, see Data Extraction.

Estimating State-Space Models

About State-Space Models.  The general state-space model structure is:

dxdt=Ax(t)+Bu(t)+Ke(t)y(t)=Cx(t)+Du(t)+e(t)

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,'InputDelay',[4 9]);

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=4, indicated by a red rectangle.

To select this model order, select 4 from the Model Order drop-down list and click Apply.

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' . You can also specify the input-to-output delay (in samples) using the 'InputDelay' property.

mCanonical = n4sid(Ze1, 3,'SSParameterization', 'canonical','InputDelay', [4 9]);
present(mCanonical);  %  Display model properties
                                                                              
mCanonical =                                                                  
  Discrete-time identified state-space model:                                 
    x(t+Ts) = A x(t) + B u(t) + K e(t)                                        
       y(t) = C x(t) + D u(t) + e(t)                                          
                                                                              
  A =                                                                         
                       x1                  x2                  x3             
   x1                   0                   1                   0             
   x2                   0                   0                   1             
   x3  0.0738 +/- 0.05922  -0.6083 +/- 0.1624    1.445 +/- 0.1284             
                                                                              
  B =                                                                         
            ConsumptionR            Current                                   
   x1   1.844 +/- 0.1756  0.5631 +/- 0.1225                                   
   x2   1.064 +/- 0.1679   2.309 +/- 0.1226                                   
   x3  0.2769 +/- 0.0966   1.878 +/- 0.1061                                   
                                                                              
  C =                                                                         
               x1  x2  x3                                                     
   Production   1   0   0                                                     
                                                                              
  D =                                                                         
               ConsumptionR       Current                                     
   Production             0             0                                     
                                                                              
  K =                                                                         
               Production                                                     
   x1  0.8676 +/- 0.03183                                                     
   x2  0.6864 +/- 0.04166                                                     
   x3  0.5113 +/- 0.04379                                                     
                                                                              
  Input delays (sampling periods): 4  9                                       
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   CANONICAL form with indices: 3.                                            
   Feedthrough: none                                                          
   Disturbance component: estimate                                            
   Number of free coefficients: 12                                            
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Estimated using N4SID on time domain data "Ze1".                              
Fit to estimation data: 91.37% (prediction focus)                             
FPE: 2.456, MSE: 2.339                                                        
More information in model's "Report" property.                                

    Note:   mn4sid and mCanonical are discrete-time models. To estimate a continuous-time model, set the 'Ts' property to 0 in the estimation command, or use the ssest command:

    mCT1 = n4sid(Ze1, 3, 'Ts', 0, 'InputDelay', [2.5 5])
    mCT2 = ssest(Ze1, 3,'InputDelay', [2.5 5])
    

Learn More.  To learn more about estimating state-space models, see State-Space Models.

Estimating a Box-Jenkins Model

About Box-Jenkins Models.  The general Box-Jenkins (BJ) structure is:

y(t)=i=1nuBi(q)Fi(q)ui(tnki)+C(q)D(q)e(t)

To estimate a BJ model, you need to specify the parameters nb, nf, nc, nd, and nk.

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)

To estimate the BJ model, type:

na = 0;
nb = [ 2 1 ];
nc = 1;
nd = 1;
nf = [ 2 1 ];
nk = [ 5 10];
mbj = polyest(Ze1,[na nb nc nd nf nk]);

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

Display the model information.

present(mbj)
                                                                              
mbj =                                                                         
Discrete-time BJ model:  y(t) = [B(z)/F(z)]u(t) + [C(z)/D(z)]e(t)             
  B1(z) = 1.823 (+/- 0.1792) z^-5 - 1.315 (+/- 0.2368) z^-6                   
                                                                              
  B2(z) = 1.791 (+/- 0.06435) z^-10                                           
                                                                              
  C(z) = 1 + 0.1066 (+/- 0.04015) z^-1                                        
                                                                              
  D(z) = 1 - 0.7453 (+/- 0.027) z^-1                                          
                                                                              
  F1(z) = 1 - 1.321 (+/- 0.06939) z^-1 + 0.5911 (+/- 0.05516) z^-2            
                                                                              
  F2(z) = 1 - 0.8314 (+/- 0.006445) z^-1                                      
                                                                              
Sample time: 0.5 minutes                                                      
                                                                              
Parameterization:                                                             
   Polynomial orders:   nb=[2 1]   nc=1   nd=1   nf=[2 1]                     
   nk=[5 10]                                                                  
   Number of free coefficients: 8                                             
   Use "polydata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Termination condition: Near (local) minimum, (norm(g) < tol).                 
Number of iterations: 6, Number of function evaluations: 13                   
                                                                              
Estimated using POLYEST on time domain data "Ze1".                            
Fit to estimation data: 90.75% (prediction focus)                             
FPE: 2.757, MSE: 2.689                                                        
More information in model's "Report" property.                                

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.

    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.

Learn More.  To learn more about identifying input-output polynomial models, such as BJ, see Input-Output Polynomial Models.

Comparing Model Output to Measured Output

Compare the output of the ARX, state-space, and Box-Jenkins models to the measured output.

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.

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.

Simulating and Predicting Model Output

Simulating the Model Output

In this portion of the tutorial, you simulate the model output. You must have already created the continuous-time model midproc2, as described in Estimating Process Models.

Simulating the model output requires the following information:

  • Input values to the model

  • Initial conditions for the simulation (also called initial states)

For example, the following commands use the iddata and idinput commands to construct an input data set, and use sim to simulate the model output:

% Create input for simulation
U = iddata([],idinput([200 2]),'Ts',0.5,'TimeUnit','min');
% Simulate the response setting initial conditions equal to zero
ysim_1 = sim(midproc2,U);

To maximize the fit between the simulated response of a model to the measured output for the same input, you can compute the initial conditions corresponding to the measured data. The best fitting initial conditions can be obtained by using findstates(idParametric) on the state-space version of the estimated model. The following commands estimate the initial states X0est from the data set Zv1:

% State-space version of the model needed for computing initial states
midproc2_ss = idss(midproc2);
X0est = findstates(midproc2_ss, Zv1);

Next, simulate the model using the initial states estimated from the data.

% Simulation input
Usim = Zv1(:,[],:);
Opt = simOptions('InitialCondition',X0est);
ysim_2 = sim(midproc2_ss, Usim, Opt);

Compare the simulated and the measured output on a plot.

figure
plot([ysim_2.y, Zv1.y])
legend({'model output','measured'})
xlabel('time'), ylabel('Output')

Predicting the Future Output

Many control-design applications require you to predict the future outputs of a dynamic system using the past input/output data.

For example, use predict to predict the model response five steps ahead:

predict(midproc2,Ze1,5)

Compare the predicted output values with the measured output values. The third argument of compare specifies a five-step-ahead prediction. When you do not specify a third argument, compare assumes an infinite prediction horizon and simulates the model output instead.

compare(Ze1,midproc2,5)

Use pe to compute the prediction error Err between the predicted output of midproc2 and the measured output. Then, plot the error spectrum using the spectrum command.

[Err] = pe(midproc2,Zv1);
spectrum(spa(Err,[],logspace(-2,2,200)))

The prediction errors are plotted with a 1-standard-deviation confidence interval. The errors are greater at high frequencies because of the high-frequency nature of the disturbance.

Was this topic helpful?