Documentation Center

  • Trial Software
  • Product Updates

Estimating Nonlinear Grey-Box Models

Specifying the Nonlinear Grey-Box Model Structure

You must represent your system as a set of first-order nonlinear difference or differential equations:

where for continuous-time representation and for discrete-time representation with Ts as the sampling interval. F and H are arbitrary linear or nonlinear functions with Nx and Ny components, respectively. Nx is the number of states and Ny is the number of outputs.

After you establish the equations for your system, create a function or MEX-file. MEX-files, which can be created in C or Fortran, are dynamically linked subroutines that can be loaded and executed by the MATLAB® interpreter. For more information about MEX-files, see MEX-File Creation API.

The purpose of the model file is to return the state derivatives and model outputs as a function of time, states, inputs, and model parameters, as follows:

[dx,y] = MODFILENAME(t,x,u,p1,p2, ...,pN,FileArgument)

    Tip   The template file for writing the C MEX-file, IDNLGREY_MODEL_TEMPLATE.c, is located in matlab/toolbox/ident/nlident.

The output variables are:

  • dx — Represents the right side(s) of the state-space equation(s). A column vector with Nx entries. For static models, dx=[].

    For discrete-time models. dx is the value of the states at the next time step x(t+Ts).

    For continuous-time models. dx is the state derivatives at time t, or .

  • y — Represents the right side(s) of the output equation(s). A column vector with Ny entries.

The file inputs are:

  • t — Current time.

  • x — State vector at time t. For static models, equals [].

  • u — Input vector at time t. For time-series models, equals [].

  • p1,p2, ...,pN — Parameters, which can be real scalars, column vectors or two-dimensional matrices. N is the number of parameter objects. For scalar parameters, N is the total number of parameter elements.

  • FileArgument — Contains auxiliary variables that might be required for updating the constants in the state equations.

    Tip   After creating a model file, call it directly from the MATLAB software with reasonable inputs and verify the output values.

For an example of creating grey-box model files and idnlgrey model object, see Creating idnlgrey Model Files.

For examples of code files and MEX-files that specify model structure, see the toolbox/ident/iddemos/examples folder. For example, the model of a DC motor is described in files dcmotor_m and dcmotor_c.

Constructing the idnlgrey Object

After you create the function or MEX-file with your model structure, you must define an idnlgrey object. This object shares many of the properties of the linear idgrey model object.

Use the following syntax to define the idnlgrey model object:

m = idnlgrey('filename',Order,Parameters,InitialStates)

The idnlgrey arguments are defined as follows:

  • 'filename' — Name of the function or MEX-file storing the model structure. This file must be on the MATLAB path when you use this model object for model estimation, prediction, or simulation.

  • Order — Vector with three entries [Ny Nu Nx], specifying the number of model outputs Ny, the number of inputs Nu, and the number of states Nx.

  • Parameters — Parameters, specified as struct arrays, cell arrays, or double arrays.

  • InitialStates — Specified in the same way as parameters. Must be the fourth input to the idnlgrey constructor.

For detailed information about this object and its properties, see the idnlgrey reference page.

Use pem to estimate your grey-box model.

Using pem to Estimate Nonlinear Grey-Box Models

You can use the pem command to estimate the unknown idnlgrey model parameters and initial states using measured data.

The input-output dimensions of the data must be compatible with the input and output orders you specified for the idnlgrey model.

Use the following general estimation syntax:

m = pem(data,m)

where data is the estimation data and m is the idnlgrey model object you constructed.

You can pass additional property-value pairs to pem to specify the properties of the model or the estimation algorithm. Assignable properties include the ones returned by the get(idnlgrey) command and the algorithm properties returned by the get(idnlgrey, 'Algorithm'), such as MaxIter and Tolerance. For detailed information about these model properties, see the idnlgrey reference page.

For more information about validating your models, see Model Validation.

Nonlinear Grey-Box Model Estimation Algorithm Options

The Algorithm property of the model specifies the estimation algorithm, which simulates the model several times by trying various parameter values to reduce the prediction error.

The following algorithm properties can affect the quality of the results:

For detailed information about these and other model properties, see the idnlgrey reference page.

Simulation Method

You can specify the simulation method using the SimulationOptions (struct) fields of the model Algorithm property.

System Identification Toolbox™ software provides several variable-step and fixed-step solvers for simulating idnlgrey models. To view a list of available solvers and their properties, type the following command at the prompt:

idprops idnlgrey algorithm.simulationoptions

For discrete-time systems, the default solver is 'FixedStepDiscrete'. For continuous-time systems, the default solver is 'ode45'.

By default, SimulationOptions.Solver is set to 'Auto', which automatically selects either 'ode45' or 'FixedStepDiscrete' during estimation and simulation—depending on whether the system is continuous or discrete in time.

Search Method

You can specify the search method for estimating model parameters using the SearchMethod field of the Algorithm property. Two categories of methods are available for nonlinear grey-box modeling.

One category of methods consists of the minimization schemes that are based on line-search methods, including Gauss-Newton type methods, steepest-descent methods, and Levenberg-Marquardt methods.

The Trust-Region Reflective Newton method of nonlinear least-squares (lsqnonlin), where the cost is the sum of squares of errors between the measured and simulated outputs, requires Optimization Toolbox™ software. When the parameter bounds differ from the default +/- Inf, this search method handles the bounds better than the schemes based on a line search. However, unlike the line-search-based methods, lsqnonlin only works with Criterion='Trace'.

By default, SearchMethod is set to Auto, which automatically selects a method from the available minimizers. If the Optimization Toolbox product is installed, SearchMethod is set to 'lsqnonlin'. Otherwise, SearchMethod is a combination of line-search based schemes.

Gradient Options

You can specify the method for calculating gradients using the GradientOptions field of the Algorithm property. Gradients are the derivatives of errors with respect to unknown parameters and initial states.

Gradients are calculated by numerically perturbing unknown quantities and measuring their effects on the simulation error.

Option for gradient computation include the choice of the differencing scheme (forward, backward or central), the size of minimum perturbation of the unknown quantities, and whether the gradients are calculated simultaneously or individually.

Example – Specifying Algorithm Properties

You can specify the Algorithm fields directly in the estimation syntax, as property-value pairs.

For example, you can specify the following properties as part of the pem syntax:

m = pem(data,init_model,'Search','gn',...

Represent Nonlinear Dynamics Using MATLAB File for Grey-Box Estimation



About the Model


Figure 1: Schematic diagram of a DC-motor.


   d         | 0      1   |        |   0   |
             | 1   0 |
      y(t) = |       | x(t)
             | 0   1 |


About the Input-Output Data

1. Load the DC motor data.


2. Represent the estimation data as an iddata object.



set(z, 'InputName', 'Voltage', 'InputUnit', 'V');
set(z, 'OutputName', {'Angular position', 'Angular velocity'});
set(z, 'OutputUnit', {'rad', 'rad/s'});
set(z, 'Tstart', 0, 'TimeUnit', 's');

4. Plot the data.

The data is shown in two plot windows.



Linear Modeling of the DC-Motor

1. Represent the DC motor structure in a function.

In this example, you use a MATLAB® file, but you can also use C MEX-files (to gain computational speed), P-files or function handles. For more information, see "Creating IDNLGREY Model Files".


  % Output equations.
       x(2)                          ... % Angular velocity.
  % State equations.
  dx = [x(2);                        ... % Angular position.

The file must always be structured to return the following:

Output arguments:

  • dx is the right-hand side(s) of the state-space equation(s)

  • y is the output equation(s)

Input arguments:



  • varargin for the auxiliary input arguments

2. Represent the DC motor dynamics using an idnlgrey object.


Ts            = 0;                 % Time-continuous system.
                'Name', 'DC-motor');


3. Specify input and output names, and units.

          'TimeUnit', 's');




5. View the initial model.

a. Get basic information about the model.



b. View the initial states and parameters.

Both the initial states and parameters are structure arrays. The fields specify the properties of an individual initial state or parameter. Type idprops idnlgrey InitialStates and idprops idnlgrey Parameters for more information.

ans = 

       Name: 'Angular position'
       Unit: 'rad'
      Value: 0
    Minimum: -Inf
    Maximum: Inf
      Fixed: 1

ans = 

       Name: 'Static gain'
       Unit: 'rad/(V*s)'
      Value: 0.2800
    Minimum: -Inf
    Maximum: Inf
      Fixed: 0



getinit(nlgr, 'Fixed')
getpar(nlgr, 'Min')
ans = 


ans = 


d. Obtain basic information about the object:

nlgr =




            FileName: 'dcmotor_m'
               Order: [1x1 struct]
          Parameters: [2x1 struct]
       InitialStates: [2x1 struct]
        FileArgument: {}
    CovarianceMatrix: 'estimate'
      EstimationInfo: [1x1 struct]
        TimeVariable: 't'
       NoiseVariance: [2x2 double]
           Algorithm: [1x1 struct]
                  Ts: 0
            TimeUnit: 'seconds'
           InputName: {'Voltage'}
           InputUnit: {'V'}
          InputGroup: [1x1 struct]
          OutputName: {2x1 cell}
          OutputUnit: {2x1 cell}
         OutputGroup: [1x1 struct]
                Name: 'DC-motor'
               Notes: {}
            UserData: []

Performance Evaluation of the Initial DC-Motor Model



nlgr.Algorithm.SimulationOptions.AbsTol = 1e-6;
nlgr.Algorithm.SimulationOptions.RelTol = 1e-5;

2. Compare the simulated output with the measured data.


The simulated and measured outputs are shown in a plot window.

compare(z, nlgr);


Parameter Estimation

Estimate the parameters and initial states using pem (Prediction-Error identification Method).

nlgr = setinit(nlgr, 'Fixed', {false false});   % Estimate the initial state.
nlgr = pem(z, nlgr, 'Display', 'Full');

Performance Evaluation of the Estimated DC-Motor Model

1. Review the information about the estimation process.

This information is stored in the EstimationInfo property of the idnlgrey object. The property also contains information about how the model was estimated, such as solver and search method, data set, and why the estimation was terminated.

ans = 

             Status: 'Estimated model (PEM)'
             Method: 'Solver: ode45; Search: lsqnonlin'
            LossFcn: 0.0011
                FPE: 0.0011
           DataName: 'DC-motor'
         DataLength: 400
             DataTs: {[0.1000]}
         DataDomain: 'Time'
    DataInterSample: {'zoh'}
         UpdateNorm: []
    LastImprovement: []
         Iterations: 4
       InitialGuess: [1x1 struct]
            Warning: ''
     EstimationTime: 5.6600



compare(z, nlgr);



na = [2 2; 2 2];
nb = [2; 2];
nk = [1; 1];
compare(z, nlgr, dcarx);


4. Check the prediction errors.


pe(z, nlgr);


5. Check the residuals ("leftovers").

Residuals indicate what is left unexplained by the model and are small for good model quality. Execute the following two lines of code to generate the residual plot. Press any key to advance from one plot to another.

  resid(z, nlgr);


6. Plot the step response.




7. Examine the model covariance.


ans =

   1.0e-04 *

    0.1521    0.0015
    0.0015    0.0007

ans =

    0.0099   -0.0004
   -0.0004    0.1094


nlgr =                                                                               
    u(1)  Voltage(t) [V]                                                             
 States:                                initial value                                
    x(1)  Angular position(t) [rad]     xinit@exp1   0.0302986   (est) in [-Inf, Inf]
    x(2)  Angular velocity(t) [rad/s]   xinit@exp1   -0.133728   (est) in [-Inf, Inf]
    y(1)  Angular position(t) [rad]                                                  
 Parameters:                       value      standard dev                           
    p1   Time-constant [s]         0.243646   0.00390033   (est) in [-Inf, Inf]      
    p2   Static gain [rad/(V*s)]   0.249645   0.00027217   (est) in [-Inf, Inf]      
The model was estimated from the data set 'DC-motor', which                          
contains 400 data samples.                                                           
Loss function 0.00107462 and Akaike's FPE 0.00108536                                 
Created:       14-Jan-2014 00:59:46                                                  
Last modified: 14-Jan-2014 01:00:04                                                  




  • Writing and using C MEX model-files.

  • Handling nonscalar parameters.

  • Impact of certain algorithm choices.


Was this topic helpful?