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.

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.

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

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

\[ \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

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

You will see the Parameters P as

so the fitted model is

**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.

MATLAB will make a discreet model with name g and display it asg=tf(1/3,[1 .66],0.1)

disp(g);

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

This equation can be modified as

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

\[ \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.

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

Yout1=Rmat*P;

%Mean Square Error

e=sqrt(mean((Yout-Yout1).^2));

disp(e);

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