Fitting Autoregressive Model into the Experimental/Plant Data

In this tutorial, we will learn how we can fit an autoregress model to an experimental data. For our case we generate data from a sample plant with some transfer function and fit a first order ARMAX model to it and compare the results.
Experimental Plant
We will use a plant with transfer function
\[G(z)=\frac{0.333}{z+0.6666}\]
Write following code in MATLAB to generate the model.

g=tf(1/3,[1 .66],0.1)
disp(g);

MATLAB will make a discreet model with name g and display it as

Transfer function:
 0.3333
--------
z + 0.66

Sampling time: 0.1

Simulation of the Plant
We will feed step input to this plant and add some noise to the output to make it more real.
Write following code into MATLAB

%simulate the system
[y t]=step(g,10); 

%add some random noise
n=0.001*randn(size(t));
y1=y+n;

%plot the data with noisy datastairs(t,y);
hold on
stairs(t,y1,'r');

x=ones(size(t)); %input step

ARMAX model
Let us fit a first order model to this data-set [x,y] with the following form

y[n+1]=a y[n] + bx[n]

This equation can be modified as
\[ \left[ \begin{array}{c}y[2]\\y[3]\\ ... \\ y[n+1] \end{array}\right]=  \left[ \begin{array}{cc}y[1] & x[1] \\y[2] & x[2]\\ ... & ... \\ y[n] & x[n] \end{array}\right] \times \left[ \begin{array}{c}a \\b \end{array}\right]\]
or

Y[n+1]=[Y[n]   X[n]] x [a b]'

Y[n+1]=R[n] x Theta

This can be solved for Theta=[a b]'  easily using least square fit method by

Theta= Inv(Y[n+1])*R[n]

So write in MATLAB the following script

Rmat=[y1(1:end-1) x(1:end-1)];
Yout=[y1(2:end)];

P=Rmat\Yout;
disp(P);

You will see the Parameters P as

 -0.6606
  0.3336  

so the fitted model is

 
y[n+1]=-0.6606 y[n] + 0.3336 x[n]

Error:

Let us calculate the fitting error in the model. The error comes due to inserted noise in the output. In real system, this can be anything: process noise, observation noise. So error is nothing but the difference between actual output and output from the fit model using the calculated value of a and b.

%simulated output
Yout1=Rmat*P;

%Mean Square Error
e=sqrt(mean((Yout-Yout1).^2));
disp(e);

This error e is the root mean square error (RMSE) which comes out to be

0.0013572                   

 which is close to the variation of the random noise we add added.

Converting a continuous statespace model to ARMA model

In many situation we need to convert a continuous statespace model in to a linear difference equation model or ARMA/ARMAX discreet time model. This post describes this conversion with some simple examples. 

Example 1: Conversion of First Order Differential Equation to Linear Difference Equation


we will start with a simple example here.  
Problem:
Let us take a simple first order differential equation: 

   dx
  -----  =  - x
   dt

This state space model is with A=-1 and B=0,

we want to discreetize it with T=1s and get a linear differential equation model;

Solution:

First calculate A1 and B1 as following:

1. Using the properties of z=e^-(st), we get that

A1=exp(AT);

A1=1/e;

 or in other words, we can say that

  A1=Linv((sI-A)^-1)  with t=T; 
  A1=Linv(1/(s+1))     =exp(-t) at (t=T)   =e^-1

then similarly

  B1= A^-1 (A1-I)B;= -1(e-1)*0 = 0

 

so the new model is

x[n+1]=A1*x[n]+B1*u[n]

x[n+1]=1/e x[n]   

So the final model in linear differential equation is

  
           x[n+1]-1/e x[n]=0

is the ARM model you wanted

Note :
       exp used is matrix exponential
       exp(A)= I+A+A^2/2!+A^3/3!+... = Linv{(sI-A)^-1}

Example 2: Considering the Inputs

Problem:
Now consider this

dx               
---- =-x+u(t);
dt               

Solution:

So we have

 A1=exp(AT);

 A1=1/e;

 or in other words, we can say that

  Linv((sI-A)^-1)  with t=T; 
 A1=Linv(1/(s+1))     =exp(-t) at (t=T)   = e^-1

then

  B1= inv(A)(A1-I)B;=-1(1/e-1)*1=1-1/e;

so the new model is

x[n+1]=A1 x[n]+B1u[n]

       x[n+1]=1/e x[n]+(1-1/e)u[n];

        x[n+1]-1/e x[n]=(1-1/e)u[n];

x[n+1]-1/e x[n]=(1-1/e)u[n];

Another Method Using Transfer Functions

Let us learn one more method for this.

first calculate the transfer function of this state space model
which is

 G(s)= 1/(s+1);

now discreetize it by using

z=exp(sT)=exp(s) as T=1;

This pretty tricky so we will take the help of MATLAB here to calculate it.

On MATLAB, write the following :

 %make the transfer function using tf
 g=tf(1,[1 1]);

 %convert it to discreet
 g1=c2d(g,1);

you will get

g1=(1-1/e)/(z-1/e);

Now it is very easy to convert it to time domain

X(z)/U(z)=(1-1/e)/(z-1/e);

z X(z)-X(z)/e=(1-1/e)U(z);

The final model is

x[n+1]-1/e x[n]= (1-1/e)u[n]

which is same as earlier one.