linear_programming()

function linear_programming()
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;
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 m/s
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

Amat=A;
l_d=length(Amat(1,:));
y_out_sysID=[];
for ind=1:length(Amat)
    if (Amat(ind,l_d)<=w_h)&(Amat(ind,l_d)>=w_l)
        y_out_sysID=[y_out_sysID;Amat(ind,:)*parsol_H];%High speed range  Parameters for High speed data
    else
        y_out_sysID=[y_out_sysID;Amat(ind,:)*parsol_L];%Low speed range  Parameters for Low speed data
    end
end
figure
Yest1=y_out_sysID;
t=linspace(1,length(y_out_sysID),length(y_out_sysID));
vfy=y_out_vfy((mb+1):end);
plot(t,vfy,'b--',t,Yest1,'r');
legend('Actual','Estimate');
ylabel('Watts');
xlabel('Time minute');
title('Identification With ARX');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Linear Programming
Ain1=[];
for k=1:1:mb-na
    Ain1=[Ain1 -y_out((mb-k+1):(end-k))];
end

for p=1:1:mb
    k=p-1;
    Ain1=[Ain1 -u_in((mb-k+1):(end-k))];
end

n=length(Ain1);
one_I=ones(n,1);
fact=0.1;% Factor to scale error convergence limit t
Ain1=[Ain1 -fact*one_I];
Ain2=[];
for k=1:1:mb-na
    Ain2=[Ain2 y_out((mb-k+1):(end-k))];
end

for p=1:1:mb
    k=p-1;
    Ain2=[Ain2 u_in((mb-k+1):(end-k))];
end
Ain2=[Ain2 -fact*one_I];

Bin1=-y_out(mb+1:end);
Bin2=y_out(mb+1:end);
Ain=[Ain1;Ain2];
Bin=[Bin1;Bin2];
F=zeros(1,mb-na+mb);
F=[F fact];
parsol=linprog(F,Ain,Bin);
Yest=Amat*parsol(1:mb-na+mb);
t_vec=linspace(1,length(Yest),length(Yest));
figure
plot(t_vec,y_out(mb+1:end),'b--',t_vec,Yest,'r');
legend('Actual','Estimate');
title('Identification With Linear Programming');
ylabel('Watts');
xlabel('Time(minute)');
ind=[];
for i=1:10
    nr_Arx(i)=norm(Yest1-vfy,i);
    nr_lin(i)=norm(Yest-y_out(mb+1:end),i);
    ind=[ind i];
end
figure
plot(ind,nr_Arx,'r--',ind,nr_lin,'b')
legend('ARX','Linear Programming')
ylabel('\fontsize{16}(\Sigmax_i ^p) ^1^/^p')
xlabel('\fontsize{16}p')
title({'Comparision of P norms between error vectors'; 'identified by ARX and Linear Programming'})