Learn more about System Identification Toolbox

Estimating Continuous-Time Transfer Functions (Process Models)

On this page…
Specifying the Structure of the Process Model Viewing the Model Structure and Parameter Values Specifying Initial Guesses for Time Delays Estimating Model Parameters Using pem Validating the Process Model Estimating a Transfer Function with a Noise Model

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 n_a=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.

  • Delay of 10, corresponding to n_k=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'});

The argument of idproc is a cell array that contains two strings, where each string specifies 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).

Back to Top

Viewing the Model Structure and Parameter Values

To view the two resulting models, type the following command in the MATLAB Command Window:

midproc0

MATLAB displays the following model structure:

Process model with 2 inputs: y = G_1(s)u_1 + G_2(s)u_2
                                                      
where                                                 
                    1+Tz*s                            
G_1(s) = Kp * ---------------------- * exp(-Td*s)      
             1+2*Zeta*Tw*s+(Tw*s)^2                   
                                                      
                                                      
with    Kp = NaN                                       
       Tw = NaN                                       
     Zeta = NaN                                       
       Td = NaN                                       
       Tz = NaN                                       
                                                      
                                                      
           Kp                                          
G_2(s) = ---------- * exp(-Td*s)                      
          1+Tp1*s                                     
                                                      
with   Kp = NaN                                        
     Tp1 = NaN                                        
      Td = NaN                                        
                                                      
This model was not estimated from data.

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

Back to Top

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.Td.value = [2.5 5];
midproc0.Td.max = [3 7];

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

Back to Top

Estimating Model Parameters Using pem

pem is an iterative estimation command, 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, pem estimates different black-box polynomial models. You can use pem, for example, to estimate parameters for linear continuous-time transfer-function, state-space, or polynomial model structures. To use pem, 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 = pem(Ze1,midproc0);
present(midproc)

MATLAB displays the following model structure with estimated parameter values:

Process model with 2 inputs: y = G_1(s)u_1 + G_2(s)u_2

where                                                 
                    1+Tz*s                            
G_1(s) = Kp * ----------------------- * exp(-Td*s)    
               1+2*Zeta*Tw*s+(Tw*s)^2                 
                                                      
with    Kp = 0.11276+-7.3308                          
        Tw = 108.08+-196.4                            
      Zeta = 1.8735+-3.5244                           
        Td = 0+-1.5255                                
        Tz = 9031.6+-5.5446e+005                      
                                                      
                                                      
           Kp                                         
G_2(s) = ---------- * exp(-Td*s)                      
          1+Tp1*s                                     
                                                      
with  Kp = 10.068+-0.062327                           
     Tp1 = 2.1189+-0.072149                           
      Td = 4.8392+-0.046529                           
                                                      
Estimated using PEM using SearchMethod = Auto from    
data set Ze1 
Loss function 8.89822 and FPE 9.04059

Unlike discrete-time polynomial models, continuous-time models let you estimate the delays. In this case, the estimated delay values 0 and 4.8392 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 the corresponding topic in the System Identification Toolbox documentation.

Back to Top

Validating the Process 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,midproc)

The preceding plot shows that the model output reasonably agrees with the measured output: there is an agreement of 65.6% 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 first 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. However, the cross-correlation with each of the two inputs shows no significant correlation.

You can try to improve the estimation results by estimating the model again with more iterations and trying different search methods. For example, you can use the model midproc as an initial guess for the next iterations:

midproc1 = pem(Ze1,midproc,'SearchMethod','lm')
compare(Zv1,midproc,midproc1)
resid(Zv1, midproc1)

The resulting plots show an improved fit of the simulated model output to the validation data and a reduced correlation in the residuals. The SearchMethod='lm' setting uses the Levenberg-Marquardt method for iterative parameter estimation.

For more information about the SearchMethod field of the Algorithm model property of linear models, see the corresponding reference page.

Back to Top

Estimating a Transfer Function with a 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 command to specify a first-order ARMA noise:

midproc2 = pem(Ze1,midproc0,'DisturbanceModel','arma1')

Note   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 and perform residual analysis, as follows:

compare(Zv1,midproc2)
figure
resid(Zv1,midproc2)

The following plot shows that the model output maintains reasonable agreement with the validation-data output. 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.

Back to Top

Estimating Model Orders Using an ARX Model Structure

Estimating Black-Box Polynomial Models

© 1984-2009- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS