%Present an ARX-model from simulation data for a multivariable case

load c:\MATLABR11\work\Diplom\MatlabHysys\NewTests\prbs3x3new.txt
%reading output and input from a loaded file,
%the name is the same as the name to the loaded file
u1=prbs3x3new(:,8);
ua=normaliz(u1);
u2=prbs3x3new(:,9);
ub=normaliz(u2);
u3=prbs3x3new(:,10);
uc=normaliz(u3);
y1=prbs3x3new(:,5);
ya=normaliz(y1);
y2=prbs3x3new(:,6);
yb=normaliz(y2);
y3=prbs3x3new(:,7);
yc=normaliz(y3);

z=[ya yb yc ua ub uc];%set the output and input in a vector
ny=3; %number of outputs
nu=3; %number of inputs

%constants
namin=1*ones(ny); %minimum degree of A polynomial
namax=5*ones(ny); %maximum degree of A polynomial
nbmin=1*ones(ny,nu); %minimum degree of B polynomial
nbmax=5*ones(ny,nu); %maximum degree of B polynomial
nkmin=1*ones(ny,nu); %%minimum degree of delay
nkmax=5*ones(ny,nu); %maximum degree of delay

%Initialized model
na=namin;  %initialize the degree of A polynomial
nb=nbmin;  %initialize the degree of A polynomial
nk=nkmin;  %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 theta" for structure of the th matrix
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];
optFPE=FPEint;
optloss_function=th(1,1);
opt_th=th;
%Check all the combinations of model structure that gives the best
%describtion of the process
for i=namin(1,1):1:namax(1,1)
   for j=nbmin(1,1):1:nbmax(1,1)
      for k= nkmin(1,1):1:nkmax(1,1)
         na=i*ones(ny,ny);
         nb=j*ones(ny,nu);
         nk=k*ones(ny,nu);
         nn=[na nb nk]; 
         th=arx(z,nn); %develop the arx-model
         FPE=th(2,1); %Displays the FPE in the th matrix
         loss_funciton=th(1,1);
         if FPE<FPEint %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];
            optFPE=FPE;
            optloss_function=th(1,1);
            opt_th=th;
         end %if
      end %for
   end %for
end %for

%optimal parameters
e=resid(z,opt_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

%Display modelparameters and other keyparamters
[Aopt,Bopt]=th2arx(opt_th)
optFPE
optloss_function
optparameters


   