## Friday, 15 April 2011

### 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);

n=0.001*randn(size(t));
y1=y+n;

%plot the data with noisy data
stairs(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.