System Identification Toolbox

Estimation of Transfer Function Models

Here we give a very brief “cheat sheet” type explanation of how you estimate standard transfer function type models of the form

where

$G(q,theta) = q^{-k}frac{B(q,theta)}{A(q,theta)},quad H(q,theta) = frac{C(q,theta)}{D(q,theta)}$

with

$A(q,theta) = 1+a_1q^{-1}+a_2q^{-2}+cdots+a_{n_a}q^{-n_a},qquad B(q,theta) = b_0+b_1q^{-1}+b_2q^{-2}+cdots+b_{n_b}q^{-n_b}$
$D(q,theta) = 1+d_1q^{-1}+d_2q^{-2}+cdots+d_{n_d}q^{-n_d},qquad C(q,theta) = 1+c_1q^{-1}+c_2q^{-2}+cdots+c_{n_c}q^{-n_c}.$

First Fundamentals

Estimation of models via the command line involves executing the command

z = est(z,m,opt);

Here z specifies the data, m specifies the model structure, and opt determines optional aspects of how the estimate g is obtained. Because it is optional, it is not strictly necessary to use it, so g=est(z,m) is perfectly legal.

Example data

In order to try this command out on various transfer function model structures, we suggest you generate some test data. For example

u = randn(1,1000);             % Random white noise input
den = poly([0.9,0.4]);         % Specify tf for system
num = [1,0.7];
num = num*sum(den)/sum(num);   % Set unit dc gain
x=filter(num,den,u);           % Generate some data
y=x+sqrt(0.01)*randn(size(x)); % Add noise to make interesting

Specifying the model structure

One way this can be done is by specifying the model orders. In the toolbox, these correspond to n_a= m.A, n_b= m.B, n_c= m.C, n_d= m.D, =m.delay. Unspecified orders are replaced with default values.

Another way, is to specify the model type, such as m.type='arx’ and then default orders are chosen. The examples below hopefully will make this a bit clearer.

How to estimate various common structure types

Estimate FIR Model
z.y=y; z.u=u;   % Specify the input-output data
m.B=10;         % Specify 2nd order model for G
m.type='fir';   % All other order to suit FIR type (eg m_a=0)

g=est(z,m);     % Do the estimation
details(g);     % Present the results

Estimate ARX Model
z.y=y; z.u=u;   % Specify the input-output data
m.A=2;          % Specify 2nd order model for G
m.type='arx';   % All other order to suit ARX type (eg m_d=m_a)

g=est(z,m);     % Do the estimation
details(g);     % Present the results

Estimate ARMAX Model
z.y=y; z.u=u;   % Specify the input-output data
m.A=2;          % Specify 2nd order model for G
m.type='armax'; % All other order to suit ARX type (eg m_d=m_a)

g=est(z,m);     % Do the estimation
details(g);     % Present the results

Estimate Output Error Model
z.y=y; z.u=u;  % Specify the input-output data
m.A=2;         % Specify 2nd order model for G
m.type='oe';   % All other order to suit OE type (eg m_c=0)

g=est(z,m);    % Do the estimation
details(g);    % Present the results

Estimate Box-Jenkins Model
z.y=y; z.u=u;  % Specify the input-output data
m.A=2;         % Specify 2nd order model for G
m.C=1;         % Specify
m.type='bj';   % All other order to suit BJ type (eg m_c=m_d)

g=est(z,m);    % Do the estimation
details(g);    % Present the results

Time Series Modelling

Estimate ARMA Model
z.y=y;          % Specify the input-output data
m.A=4;          % Specify 2nd order model for G
m.type='ar';    % All other order to suit FIR type (eg m_a=0)

g=est(z,m);     % Do the estimation
details(g);     % Present the results

Estimate ARMA Model
z.y=y;          % Specify the input-output data
m.A=4;          % Specify 2nd order model for G
m.type='arma';  % All other order to suit FIR type (eg m_a=0)

g=est(z,m);     % Do the estimation
details(g);     % Present the results

Delta operator models

The toolbox supports estimation with respect to the Euler divided difference (delta) operator

$delta = frac{q-1}{T}$

where is the sampling period in seconds. To employ this, simply specify m.op='d’ in any of the above examples. To go back to the shift operator, specify m.op='q’.