Documentation |
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 time-series data (output-only).
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.
The first step in grey-box modeling is to write a function 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 the file:
[A,B,C,D] = myfunc(par1,par2,...,parN,Ts,aux1,aux2,...)
where the output arguments are the state-space matrices and 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 state-space 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 continuous-time, or a discrete-time state-space model, or both. When you create the idgrey model using myfunc, you can specify the nature of the output arguments of myfunc. The continuous-time state-space model has the form:
In continuous-time, the state-space 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 discrete-time, the state-space 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 discrete-time 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 MEX-file 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 continuous-time 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 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:
$$\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 user-defined parameters and contains their nominal (initial) values. In this example, all the scalar-valued parameters are grouped in the par vector. The scalar-valued 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 continuous-time model.
Use greyest to estimate the grey-box 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 Continuous-Time State-Space Models, where the same problem is solved using a structured state-space representation. |
This example shows how to 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 $$\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 heat-diffusion 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 continuous-time state-space model, you can represent the second-derivative 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 state-space 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 heat-diffusion equation is mapped to the following continuous-time state-space model structure to identify the state-space 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 state-space 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 state-space 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 state-space matrices: % Ngrid: Number of points in the space-discretization % 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: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 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 grey-box 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 grey-box model properties, see the idgrey reference page.
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.
This example is based on a discrete, single-input and single-output (SISO) system represented by the following state-space 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 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.
You can represent the system described in Description of the SISO System as an idgrey (grey-box) 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 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:
Create the 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. The initial conditions are not parameterized; they are assumed to be zero during this estimation.
Note: R2 is treated as an auxiliary variable rather than assigned a value in the file to let you change this value directly at the command line and avoid editing the file. |
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 discrete-time system.
Estimate the model parameter values from data:
opt = greyestOptions; opt.InitialState = 'zero'; opt.Display = 'full'; Model = greyest(data,Minit,opt)