systemID_project_ARMAX()

function systemID_project_ARMAX()
clc;
clear all;
close all;
format long;
load 'Turbinedata_sky.m';
%%DATA Sequence
%%time(m) time(d) time(yyyy) time(h) time(mm) time(ss)	rpm	powero	powerr	windspeed	watthrs

Tdata=Turbinedata_sky(:,:);
N1=1032;
N2=2926;
count=400;% Number of minute data points to extract
w_l=3;% Low speed limit
w_h=20;% High Speed limit
w_filter1=w_l;
w_filter2=w_h;
% Function to extract minute data
[kwh wind]=minute_data(Tdata,count,N1,N2,w_filter1,w_filter2);
u_in=wind;
y_out=kwh;
mb=4;% Higest order of the system
na=1;
[parsol_H]=ARX_model(mb,na,u_in,y_out)%Parameters solved by wind speed 3 to 20
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Low wind speed Parameter calculation
count=600;
w_filter1=0;
w_filter2=w_l;
[kwh wind]=minute_data(Tdata,count,N1,N2,w_filter1,w_filter2);
u_in=wind;
y_out=kwh;
[parsol_L]=ARX_model(mb,na,u_in,y_out)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
count=200;%Extract Data points for verification
% Wind Speed 0 to 20
w_filter1=0;
w_filter2=20;
[kwh wind]=minute_data(Tdata,count,N1,N2,w_filter1,w_filter2);
u_in=wind;
y_out=kwh;
L=length(u_in);
u_in_ID=u_in;
u_in_vfy=u_in;

y_out_ID=y_out;
y_out_vfy=y_out;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n=length(y_out_vfy)-mb;
I=eye(n,1)+1;
I(1)=I(1)-1; 
A=I;
for k=1:1:mb-na
    A=[A y_out_vfy((mb-k+1):(end-k))];
end

for p=1:1:mb
    k=p-1;
    A=[A u_in_vfy((mb-k+1):(end-k))];
end
A(:,1)=[]; % Delete 1st column of Matrix A, which was used to Initialize it


l_d=length(A(1,:));
y_out_sysID=[];
for ind=1:length(A)
    if (A(ind,l_d)<=w_h)&(A(ind,l_d)>=w_l)
        y_out_sysID=[y_out_sysID;A(ind,:)*parsol_L];
    else
        y_out_sysID=[y_out_sysID;A(ind,:)*parsol_L];
    end
end
figure
t=linspace(1,length(y_out_sysID),length(y_out_sysID));
plot(t,y_out_sysID,'r','linewidth',1.5);
hold on;
vfy=y_out_vfy((mb+1):end);
plot(t,vfy,'g');
title('Predicton with only Low Wind Speed data: ARX');
legend('Identified Watts','Actual Watts');
ylabel('Watts');
xlabel('Time(min)');

l_d=length(A(1,:));
y_out_sysID=[];
for ind=1:length(A)
    if (A(ind,l_d)<=w_h)&(A(ind,l_d)>=w_l)
        y_out_sysID=[y_out_sysID;A(ind,:)*parsol_H];
    else
        y_out_sysID=[y_out_sysID;A(ind,:)*parsol_H];
    end
end
figure
t=linspace(1,length(y_out_sysID),length(y_out_sysID));
plot(t,y_out_sysID,'r','linewidth',1.5);
hold on;
vfy=y_out_vfy((mb+1):end);
plot(t,vfy,'g');
title('Predicton with only High Wind Speed data: ARX');
legend('Identified Watts','Actual Watts');
ylabel('Watts');
xlabel('Time(min)');

l_d=length(A(1,:));
y_out_sysID=[];
for ind=1:length(A)
    if (A(ind,l_d)<=w_h)&(A(ind,l_d)>=w_l)
        y_out_sysID=[y_out_sysID;A(ind,:)*parsol_H];
    else
        y_out_sysID=[y_out_sysID;A(ind,:)*parsol_L];
    end
end
figure
t=linspace(1,length(y_out_sysID),length(y_out_sysID));
plot(t,y_out_sysID,'r','linewidth',1.5);
hold on;
vfy=y_out_vfy((mb+1):end);
plot(t,vfy,'g');
title('Predicton with both High and Low Wind Speed data: ARX');
legend('Identified Watts','Actual Watts');
ylabel('Watts');
xlabel('Time(min)');
vfy_arx=vfy;
sysID_arx=y_out_sysID;
% %---------------------------------------------------------------------
%ARMAX MODEL
Y_actual=y_out_ID((mb+1):end); % Since 1st mb values are used in delay modeling in ARX
u=u_in_ID((mb+1):end);
Y_verify=y_out_vfy((mb+1):end);
U_verify=u_in_vfy((mb+1):end);
error =vfy- y_out_sysID;% Error vector

    n=length(Y_actual)-mb;
    I=eye(n,1)+1;
    I(1)=I(1)-1;
    A=I;% Initialize A
    
    Y=Y_actual((mb+1):end);
    for k=1:1:mb-na
        A=[A Y_actual((mb-k+1):(end-k))];
    end

    for p=1:1:mb
        k=p-1;
        A=[A u((mb-k+1):(end-k))];
    end

    for p=1:1:mb
        k=p;
        A=[A error((mb-k+1):(end-k))];
    end

    A(:,1)=[];% Delete 1st column of Matrix A, which was used to Initialize it
parsol=pinv(A)*Y
% Generate Identified Output vector based on previous
% outputs, previous errors, current and previous Inputs and Parameters solved
% by Least square method
n=length(Y_verify)-mb;
I=eye(n,1)+1;
I(1)=I(1)-1;
A=I;
Y=Y_verify((mb+1):end);
for k=1:1:mb-na
    A=[A Y_verify((mb-k+1):(end-k))];
end

for p=1:1:mb
    k=p-1;
    A=[A U_verify((mb-k+1):(end-k))];
end

for p=1:1:mb
    k=p;
    A=[A error((mb-k+1):(end-k))];
end
A(:,1)=[];
Y_sysID=A*parsol;

figure
t=linspace(1,length(Y_sysID),length(Y_sysID));
plot(t,Y_sysID,'r','linewidth',1.5);
hold on;
vfy=Y_verify((mb+1):end);
plot(t,vfy,'g');
title('Identification of Watts of Wind Turbine: ARMAX');
legend('Identified Watts','Actual Watts');
ylabel('Watts');
xlabel('Time(min)');
%%% Compare P norm of error vector from ARX and ARMAX
ind=[];
for i=1:10
    nr_Arx(i)=norm(vfy_arx-sysID_arx,i);
    nr_armax(i)=norm(vfy-Y_sysID,i);
    ind=[ind i];
end
figure
plot(ind,nr_Arx,'r--',ind,nr_armax,'b')
legend('ARX','ARMAX')
ylabel('\fontsize{16}(\Sigmax_i ^p) ^1^/^p')
xlabel('\fontsize{16}p')
title({'Comparision of P norms between error vectors'; 'identified by ARX andARMAX'})