You can estimate linear discretetime and continuoustime greybox models for arbitrary ordinary differential or difference equations using singleoutput and multipleoutput timedomain data, or timeseries data (outputonly).
You must represent your system equations in statespace form. Statespace models use state variables x(t) to describe a system as a set of firstorder differential equations, rather than by one or more nthorder differential equations.
The first step in greybox modeling is to write a function that returns statespace matrices as a function of userdefined parameters and information about the model.
Use the following format to implement the linear greybox model in the file:
[A,B,C,D] = myfunc(par1,par2,...,parN,Ts,aux1,aux2,...)
myfunc
is
the name of the file. par1,par2,...,parN
are the N parameters
of the model. Each entry may be a scalar, vector or matrix.Ts
is
the sample time. aux1,aux2,...
are the optional
input arguments that myfunc
uses to compute the
statespace matrices in addition to the parameters and sample time. 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 file.
You can write the contents of myfunc
to parameterize
either a continuoustime, or a discretetime statespace model, or
both. When you create the idgrey
model using myfunc
,
you can specify the nature of the output arguments of myfunc
.
The continuoustime statespace model has the form:
In continuoustime, the statespace description has the following form:
$$\begin{array}{l}\dot{x}(t)=Ax(t)+Bu(t)+Ke(t)\\ y(t)=Cx(t)+Du(t)+e(t)\\ x(0)=x0\end{array}$$
Where, A,B,C and D are
matrices that are parameterized by the parameters par1,par2,...,parN
.
The noise matrix K and initial state vector, x0,
are not parameterized by myfunc
. They can still
be estimated as additional quantities along with the model parameters.
To configure the handling of initial states, x0,
and the disturbance component, K, during estimation,
use the greyestOptions
option
set.
In discretetime, the statespace description has a similar form:
$$\begin{array}{l}x(k+1)=Ax(k)+Bu(k)+Ke(t)\\ y(k)=Cx(k)+Du(k)+e(t)\\ x(0)=x0\end{array}$$
Where, A,B,C and D are
now the discretetime matrices that are parameterized by the parameters par1,par2,...,parN
. K and x0 are
not directly parameterized, but can be estimated if required by configuring
the corresponding estimation options.
In some applications, you may want to express K and x0 as
quantities that are parameterized by chosen parameters, just as the A, B, C and D matrices.
To handle such cases, you can write the ODE file, myfunc
,
to return K and x0 as additional
output arguments:
[A,B,C,D,K,x0] = myfunc(par1,par2,...,parN,Ts,aux1,aux2,...)
K and x0 are thus treated
in the same way as the A, B, C and D matrices.
They are all functions of the parameters par1,par2,...,parN
.
After creating the function or MEXfile with your model structure,
you must define as idgrey
object. For information
regarding creating this, see idgrey
.
This example shows how to represent the structure of the following continuoustime model:
$$\begin{array}{l}\dot{x}(t)=\left[\begin{array}{cc}0& 1\\ 0& {\theta}_{1}\end{array}\right]x(t)+\left[\begin{array}{c}0\\ {\theta}_{2}\end{array}\right]u(t)\\ y(t)=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]x(t)+e(t)\\ x(0)=\left[\begin{array}{c}{\theta}_{3}\\ 0\end{array}\right]\end{array}$$
This equation represents an electrical motor, where $${y}_{1}(t)={x}_{1}(t)$$ is the angular position of the motor shaft, and $${y}_{2}(t)={x}_{2}(t)$$ is the angular velocity. The parameter $${\theta}_{1}$$ is the inverse time constant of the motor, and $${\scriptscriptstyle \raisebox{1ex}{${\theta}_{2}$}\!\left/ \!\raisebox{1ex}{${\theta}_{1}$}\right.}$$ is the static gain from the input to the angular velocity.
The motor is at rest at t=0, but its angular position $${\theta}_{3}$$ is unknown. Suppose that the approximate nominal values of the unknown parameters are $${\theta}_{1}=1$$,$${\theta}_{2}=0.25$$ and $${\theta}_{3}=0$$. For more information about this example, see the section on statespace models in System Identification: Theory for the User, Second Edition, by Lennart Ljung, Prentice Hall PTR, 1999.
The continuoustime statespace model structure is defined by the following equation:
$$\begin{array}{l}\dot{x}(t)=Fx(t)+Gu(t)+\tilde{K}w(t)\\ y(t)=Hx(t)+Du(t)+w(t)\\ x(0)=x0\end{array}$$
To prepare this model for identification:
Create the following file to represent the model structure in this example:
function [A,B,C,D,K,x0] = myfunc(par,T)
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];
Use the following syntax to define
an idgrey
model object based
on the myfunc
file:
par = [1; 0.25; 0]; aux = {}; T = 0; m = idgrey('myfunc',par,'c',aux,T);
where par
represents
a vector of all the userdefined parameters and contains their nominal
(initial) values. In this example, all the scalarvalued parameters
are grouped in the par
vector. The scalarvalued
parameters could also have been treated as independent input arguments
to the ODE function myfunc
.'c'
specifies
that the underlying parameterization is in continuous time. aux
represents
optional arguments. As myfunc
does not have any
optional arguments, use aux = {}
. T
specifies
the sample interval; T = 0
indicates a continuoustime
model.
Use greyest
to
estimate the greybox parameter values:
m_est = greyest(data,m)
where data
is
the estimation data and m
is an estimation initialization idgrey
model. m_est
is
the estimated idgrey
model.
Note: Compare this example to Estimate Structured ContinuousTime StateSpace Models, where the same problem is solved using a structured statespace representation. 
This example shows how to estimate the heat conductivity and the heattransfer coefficient of a continuoustime greybox model for a heatedrod system.
This system consists of a wellinsulated metal rod of length L and a heatdiffusion coefficient $$\kappa $$. 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 heatdiffusion equation—which is a partial differential equation in space and time.
$$\frac{\partial x(t,\xi )}{\partial t}=\kappa \frac{{\partial}^{2}x(t,\xi )}{\partial {\xi}^{2}}$$
To get a continuoustime statespace model, you can represent the secondderivative using the following difference approximation:
$$\begin{array}{l}\frac{{\partial}^{2}x(t,\xi )}{\partial {\xi}^{2}}=\frac{x\left(t,\xi +\Delta L\right)2x(t,\xi )+x\left(t,\xi \Delta L\right)}{{\left(\Delta L\right)}^{2}}\\ \text{where}\xi =k\cdot \Delta L\end{array}$$
This transformation produces a statespace model of order $$n={\scriptscriptstyle \frac{L}{\Delta L}}$$, where the state variables $$x(t,k\cdot \Delta L)$$ are lumped representations for $$x(t,\xi )$$ for the following range of values:
$$k\cdot \Delta L\le \xi <\left(k+1\right)\Delta L$$
The dimension of x depends on the spatial grid size $$\Delta L$$ in the approximation.
The heatdiffusion equation is mapped to the following continuoustime statespace model structure to identify the statespace matrices:
$$\begin{array}{l}\dot{x}(t)=Fx(t)+Gu(t)+\tilde{K}w(t)\\ y(t)=Hx(t)+Du(t)+w(t)\\ x(0)=x0\end{array}$$
The statespace matrices are parameterized by the heat diffusion coefficient κ and the heat transfer coefficient at the far end of the rod h_{tf}. The expressions also depend upon the grid size, Ngrid, and the length of the rod L. The initial conditions x0 are a function of the initial room temperature, treated as a known quantity in this example.
Create a MATLAB^{®} file.
The following code describes the statespace equation for this
model. The parameters are κ and h_{tf} while
the auxiliary variables are Ngrid, L and
initial room temperature temp
. The grid size is
supplied as an auxiliary variable so that the ODE function can be
easily adapted for various grid sizes.
function [A,B,C,D,K,x0] = heatd(kappa,htf,T,Ngrid,L,temp) % ODE file parameterizing the heat diffusion model % kappa (first parameter)  heat diffusion coefficient % htf (second parameter)  heat transfer coefficient % at the far end of rod % Auxiliary variables for computing statespace matrices: % Ngrid: Number of points in the spacediscretization % L: Length of the rod % temp: Initial room temperature (uniform) % Compute space interval deltaL = L/Ngrid; % A matrix A = zeros(Ngrid,Ngrid); for kk = 2:Ngrid1 A(kk,kk1) = 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,Ngrid1) = 1; A(Ngrid,Ngrid) = 1; A = A*kappa/deltaL/deltaL; % B matrix B = zeros(Ngrid,1); B(Ngrid,1) = htf/deltaL; % C matrix C = zeros(1,Ngrid); C(1,1) = 1; % D matrix (fixed to zero) D = 0; % K matrix: fixed to zero K = zeros(Ngrid,1); % Initial states: fixed to room temperature x0 = temp*ones(Ngrid,1);
Use the following syntax to define an idgrey
model object based on the heatd
code
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 (grid size) 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 greyest
to estimate the greybox parameter
values:
me = greyest(data,m)
The following command shows how you can specify to estimate a new model with different auxiliary variables:
m.Structure.ExtraArgs = {20,1,22}; me = greyest(data,m);
This syntax uses the ExtraArgs
model structure
attribute to specify a finer grid using a larger value for Ngrid
.
For more information about linear greybox model properties, see the idgrey
reference page.
This example shows how to create a singleinput and singleoutput
greybox 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.
This example is based on a discrete, singleinput and singleoutput (SISO) system represented by the following statespace equations:
$$\begin{array}{l}x(kT+T)=\left[\begin{array}{cc}par1& par2\\ 1& 0\end{array}\right]x(kT)+\left[\begin{array}{c}1\\ 0\end{array}\right]u(kT)+w(kT)\\ y(kT)=\left[\begin{array}{cc}par3& par4\end{array}\right]x(kT)+e(kT)\\ x(0)=x0\end{array}$$
where w and e are independent whitenoise terms with covariance matrices R1 and R2, respectively. R1=E{ww'} is a 2–by2 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.
You can represent the system described in Description of the SISO System as an idgrey
(greybox)
model using a function. Then, you can use this file and the greyest
command
to estimate the model parameters based on initial parameter guesses.
To run this example, you must load an inputoutput data set
and represent it as an iddata
or idfrd
object
called data
. For more information about this operation,
see Representing Time and FrequencyDomain Data Using iddata
Objects or Representing FrequencyResponse Data Using idfrd Objects.
To estimate the parameters of a greybox model:
Create the file mynoise
that
computes the statespace matrices as a function of the five unknown
parameters and the auxiliary variable that represents the known variance R2
.
The initial conditions are not parameterized; they are assumed to
be zero during this estimation.
Note:

function [A,B,C,D,K] = 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]; [~,K] = kalman(ss(A,eye(2),C,0,T),R1,R2); % Uses Control System Toolbox product % u=[]
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
Construct an idgrey
model using
the mynoise
file:
Minit = idgrey('mynoise',Pvec,'d',auxVal);
The third input argument 'd'
specifies a
discretetime system.
Estimate the model parameter values from data:
opt = greyestOptions; opt.InitialState = 'zero'; opt.Display = 'full'; Model = greyest(data,Minit,opt)