%Present an ARX-model from experimental data

load c:\MATLABR11\work\Diplom\prbsinclslave.txt
%reading output and input from a loaded file,
%the name is the same as the name to the loaded file
y=prbsinclslave(:,4); %select 2nd coloumn which contains temp cold water outlet
u=prbsinclslave(:,5); %select 2nd coloumn which contains the control signal
z=[y u]; %set the output and input in a vector
%constants
namin=1; %minimum degree of A polynomial
namax=5; %maximum degree of A polynomial
nbmin=1; %minimum degree of B polynomial
nbmax=5; %maximum degree of B polynomial
nkmin=1; %%minimum degree of delay
nkmax=5; %maximum degree of delay

%Initialized model
na=1;  %initialize the degree of A polynomial
nb=1;  %initialize the degree of A polynomial
nk=1;  %initialize the degree of delay
nn=[na nb nk]; %put the degrees in a vector
th=arx(z,nn); %develop the arx-model
FPEint=th(2,1);%initialize Aikake's Final Prediction Error.
%See help th for structure of the th matrix

%Check all the combinations of model structure that gives the best describtion of the process
for na=namin:1:namax
   for nb=nbmin:1:nbmax
      for nk=nkmin:1:nkmax
         nn=[na nb nk]; 
			th=arx(z,nn);%develop the arx-model
      	FPE=th(2,1); %Displays the FPE in the th matrix
      		if FPEint>FPE %checks the FPE, less value for FPE gives better model
         		FPEint=FPE; %update the FPE value if its less than the past value
         		naopt=na; %save the optimal degree of A
         		nbopt=nb; %save the optimal degree of B
      	   	nkopt=nk; %save the optimal degree of delay
     		   	optparameters=[naopt nbopt nkopt FPE];
           end      
       end
    end
end
%Displays the optimal modelparameters na,nb,nk and FPE and the modelcoeff. A and B
display('The optimal model paramters')
th=arx(z,[naopt nbopt nkopt]); %Find the arx-model with the optimal paramters
%present(th) %displays the paramters
display(optparameters)
[A,B]=th2arx(th)
loss_funciton=th(1,1)  %See help theta for exact structure of th
e=resid(z,th); %compute the residuals and display correlation graphs
stdev_e=std(e) %finds the standard devation of the error
var_e=var(e) %finds the variance of the error

%Transform the arx model to a transfer function
%[num_tf,den_tf]=th2tf(th)

