Products & Services Industries Academia Support User Community Company

Learn more about System Identification Toolbox   

Estimating Linear Grey-Box Models

Specifying the Linear Grey-Box Model Structure

You can estimate linear discrete-time and continuous-time grey-box models for arbitrary ordinary differential or difference equations using single-output and multiple-output time-domain data, or output-only time-series data.

You must represent your system equations in state-space form. State-space models use state variables x(t) to describe a system as a set of first-order differential equations, rather than by one or more nth-order differential equations.

In continuous-time, the state-space description has the following form:

The discrete-time state-space model structure is often written in the innovations form:

The first step in grey-box modeling is to write an M-file that returns state-space matrices as a function of user-defined parameters and information about the model.

Use the following format to implement the linear grey-box model in an M-file:

[A,B,C,D,K,x0] = myfunc(par,T,CDmfile,aux)

where the matrices A, B, C, D, K, and x0 represent both the continuous-time and discrete-time state-space description of the system, myfunc is the name of the M-file, par contains the parameters as a column vector, and T is the sampling interval. aux contains auxiliary variables in your system. You use auxiliary variables to vary system parameters at the input to the function, and avoid editing the M-file.

CDmfile is an optional argument that describes whether the resulting state-space matrices are in discrete time or continuous time. By default, CDmfile='cd', which means that the sampling interval property of the model Ts determines whether the model is continuous or discrete in time. For more information about these arguments, see the idgrey reference page.

Use pem to estimate your grey-box model.

Example – Representing a Grey-Box Model in an M-File

In this example, you represent the structure of the following continuous-time model:

This equation represents an electrical motor, where is the angular position of the motor shaft, and is the angular velocity. The parameter is the inverse time constant of the motor, and is the static gain from the input to the angular velocity.

The motor is at rest at t=0, but its angular position is unknown. Suppose that the approximate nominal values of the unknown parameters are and . For more information about this example, see the section on state-space models in System Identification: Theory for the User, Second Edition, by Lennart Ljung, Prentice Hall PTR, 1999.

The continuous-time state-space model structure is defined by the following equation:

To prepare this model for identification:

  1. Create the following M-file to represent the model structure in this example:

    function [A,B,C,D,K,x0] = myfunc(par,T,aux)
    A = [0 1; 0 par(1)]; 
    B = [0;par(2)];
    C = eye(2);
    D = zeros(2,1);
    K = zeros(2,2);
    x0 =[par(3);0];
    
  2. Use the following syntax to define an idgrey model object based on the myfunc M-file:

    m = idgrey('myfunc',par,'c',T,aux)

    where par represents user-defined parameters and contains their nominal (initial) values. 'c' specifies that the underlying parameterization is in continuous time. aux contains the values of the auxiliary parameters.

      Note   You must specify T and aux even if they are not used by the myfunc code.

Use pem to estimate the grey-box parameter values:

m = pem(data,m)

where data is the estimation data and m is the idgrey object with unknown parameters.

Example – Estimating a Continuous-Time Grey-Box Model for Heat Diffusion

In this example, you estimate the heat conductivity and the heat-transfer coefficient of a continuous-time grey-box model for a heated-rod system.

This system consists of a well-insulated metal rod of length L and a heat-diffusion coefficient . The input to the system is the heating power u(t) and the measured output y(t) is the temperature at the other end.

Under ideal conditions, this system is described by the heat-diffusion equation—which is a partial differential equation in space and time.

To get a continuous-time state-space model, you can represent the second-derivative using the following difference approximation:

This transformation produces a state-space model of order , where the state variables are lumped representations for for the following range of values:

The dimension of x depends on the spatial grid size in the approximation.

The heat-diffusion equation is mapped to the following continuous-time state-space model structure to identify the state-space matrices:

The following M-file describes the state-space equation for this model. In this case, the auxiliary variables specify grid-size variables, so that you can modify the grid size without the M-file.

function [A,B,C,D,K,x0] = heatd(pars,T,aux)
% Number of points in the space-discretization
Ngrid = aux(1);
% Length of the rod
L = aux(2);
% Initial rod temperature (uniform) 
temp = aux(3);
% Space interval
deltaL = L/Ngrid;
% Heat-diffusion coefficient
kappa = pars(1);
% Heat transfer coefficient at far end of rod
htf = pars(2);
A = zeros(Ngrid,Ngrid);
for kk = 2:Ngrid-1
    A(kk,kk-1) = 1;
    A(kk,kk) = -2;
    A(kk,kk+1) = 1;
end
% Boundary condition on insulated end
A(1,1) = -1; A(1,2) = 1; 
A(Ngrid,Ngrid-1) = 1;
A(Ngrid,Ngrid) = -1;
A = A*kappa/deltaL/deltaL;
B = zeros(Ngrid,1); 
B(Ngrid,1) = htf/deltaL;
C = zeros(1,Ngrid);
C(1,1) = 1;
D = 0;
K = zeros(Ngrid,1);
x0 = temp*ones(Ngrid,1);

Use the following syntax to define an idgrey model object based on the heatd M-file:

m = idgrey('heatd',[0.27 1],'c',[10,1,22])

This command specifies the auxiliary parameters as inputs to the function, include the model order 10, the rod length of 1 meter, and an initial temperature of 22 degrees Celsius. The command also specifies the initial values for heat conductivity as 0.27, and for the heat transfer coefficient as 1.

For given data, you can use pem to estimate the grey-box parameter values:

me = pem(data,m)

The following command shows how you can specify to estimate a new model with different auxiliary variables directly in the estimator command:

me = pem(data,m,'FileArgument',[20,1,22])

This syntax uses the FileArgument model property to specify a finer grid using a larger value for Ngrid. For more information about linear grey-box model properties, see the idgrey reference page.

Example – Estimating a Discrete-Time Grey-Box Model with Parameterized Disturbance

This example shows how to create a single-input and single-output grey-box model structure when you know the variance of the measurement noise. The code in this example uses the Control System Toolbox command kalman for computing the Kalman gain from the known and estimated noise variance.

Description of the SISO System

This example is based on a discrete, single-input and single-output (SISO) system represented by the following state-space equations:

where w and e are independent white-noise terms with covariance matrices R1 and R2, respectively. R1=E{ww'} is a 2–by-2 matrix and R2=E{ee'} is a scalar. par1, par2, par3, and par4 represent the unknown parameter values to be estimated.

Assume that you know the variance of the measurement noise R2 to be 1. R1(1,1) is unknown and is treated as an additional parameter par5. The remaining elements of R1 are known to be zero.

Estimating the Parameters of an idgrey Model

You can represent the system described in Description of the SISO System as an idgrey (grey-box) model using an M-file. Then, you can use this M-file and the pem command to estimate the model parameters based on initial parameter guesses.

To run this example, you must load an input-output data set and represent it as an iddata or idfrd object called data. For more information about this operation, see Representing Time- and Frequency-Domain Data Using iddata Objects or Representing Frequency-Response Data Using idfrd Objects.

To estimate the parameters of a grey-box model:

  1. Create the M-file mynoise that computes the state-space matrices as a function of the five unknown parameters and the auxiliary variable that represents the known variance R2.

      Note   R2 is treated as an auxiliary variable rather than assigned a value in the M-file to let you change this value directly at the command line and avoid editing the M-file.

    function [A,B,C,D,K,x0] = mynoise(par,T,aux)
    R2 = aux(1); % Known measurement noise variance
    A = [par(1) par(2);1 0];
    B = [1;0];
    C = [par(3) par(4)];
    D = 0;
    R1 = [par(5) 0;0 0];
    [est,K] = kalman(ss(A,eye(2),C,0,T),R1,R2); 
    					 % Uses Control System Toolbox product
    					 % u=[]
    x0 = [0;0];
  2. Specify initial guesses for the unknown parameter values and the auxiliary parameter value R2:

    par1 = 0.1; % Initial guess for A(1,1) 
    par2 = -2;  % Initial guess for A(1,2) 
    par3 = 1;   % Initial guess for C(1,1) 
    par4 = 3;   % Initial guess for C(1,2) 
    par5 = 0.2; % Initial guess for R1(1,1)
    Pvec = [par1; par2; par3; par4; par5]
    auxVal = 1; % R2=1
  3. Construct an idgrey model using the mynoise M-file:

    Minit = idgrey('mynoise',Pvec,'d',auxVal);

    The third input argument 'd' specifies a discrete-time system.

  4. Estimate the model parameter values from data:

    Model = pem(data,Minit)

  


Recommended Products

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

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