Project_Example_09.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sean Semper                                         System Identification
% Project Code 1                                                    MAE 566
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Last Modified: 04/26/09                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ARX Test Model Example from MATLAB
clc
clear all
close all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rpm	powero	powerr	windspeed
% load in_out_dat2.txt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% avg_rpm	avg_powero	avg_powerr	avg_windspeed	nettwatthrs        %
% load in_out_dat4.txt                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time%	nobs	avg_rpm 	avg_powero	avg_powerr	avg_windspeed      
load min_avg_3109.txt
data = min_avg_3109;
numHnkl = 100;  % Hankel Matrix with n, Markov Parameters
mm = 1;         % Number of States to Keep
in_num = 5;     % Input Index
out_num = 2;    % Output Index
in_out_ind = [in_num, out_num];
u = data(:,in_num);
y = data(:,out_num);
clear data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% time%	nobs	(7)avg_rpm 	(8)avg_powero	(9)avg_powerr
% (10)avg_windspeed     
% load Turbinedata_sky.mat -ascii   
%
% numHnkl = 40;   % Hankel Matrix with n, Markov Parameters
% mm = 3          % Number of States to Keep
% in_num = 10;     % Input Index
% out_num = 7;    % Output Index
% in_out_ind = [in_num, out_num];
% 
% u =Turbinedata_sky(1032:19260,in_num);
% y = Turbinedata_sky(1032:19260,out_num);
% pause
% clear Turbinedata_sky
% pause
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  break
% load Sheet1
% load Sheet2
% PreWindDAT = Sheet1;
% load in_out_data.mat

% v = min_avg(:,3);
% u = (0.5*1.1467*((3.72)^2)*((u*3.72)./v).*(v.^3))./u;

% break

L = length(u);
dfftu = fft(u)/L;
dffty = fft(y)/L;
dt = 1 ;%  1.91462 2.782647931
df = 1/(L*dt);
fmax = 1/(2*(dt));
omg = 0:df:fmax;

SYU = dffty.*conj(dfftu);
SUU = dfftu.*conj(dfftu);

% G = [zeros(1,5000),zeros(1,5000),zeros(1,5000),zeros(1,5000),zeros(1,5000)];
G = []

tic
for i = 1:length(SUU)
    G(i) = SYU(i)/SUU(i);
end
toc

invdfftG = ifft(G);


Yii=invdfftG(1:numHnkl);
display('Hankel Matrix of Markov Parameters Using the Frequency Response')
% Form the Hankel Matrix of Markov Parameters
H00 = hankel(Yii(1,1:numHnkl/2 ),Yii(1,numHnkl/2:numHnkl-1));
H11 = hankel(Yii(1,2:numHnkl/2+1 ),Yii(1,numHnkl/2+1:numHnkl));

clear R SS S P Q 
[R,SS,S] = svd(H00);

figure(3)
stem(diag(SS));
diag(SS)
ylabel('Eigenvalue ( \sigma_i ) Magnitude')
xlabel('\sigma_i')
title({'Using the Frequency Response';'(Considering only \sigma_1, & \sigma_2)'})


P = R*sqrtm(SS);
Q = sqrtm(SS)*S';

display('ERA Estimates for A, B, C, D = [0] ')
C_est = P(1,:)*1 %*2.782647931;
B_est = Q(:,1)*1 %*2.782647931;

A_est = inv(sqrtm(SS))*R'*(H11)*S*inv(sqrtm(SS));
A_hat = A_est(1:mm,1:mm)
B_hat = B_est(1:mm)
C_hat = C_est(1:mm)


save A_hat.mat A_hat
save B_hat.mat B_hat
save C_hat.mat C_hat
save in_out_ind.mat in_out_ind

Est_A_Eig = eig(A_est)



% Q_bar = sqrt(SS)*S';
% 
% [EigVec EigVal] = eig(A_est);
% 
% B_hat_m = inv(EigVec)*B_hat
 
% for n=1:numHnkl/2
% for ii = 1:numHnkl/2
% q_hat(n,ii) =  (EigVal(n)^(ii-1))*B_hat_m(n);

% end
% end
% 
% for kk=1:mm
% MAC(kk) = Q_bar(kk,:)*conj(q_hat(kk,:)')/sqrt((Q_bar(kk,:)*conj(Q_bar(kk,:)'))*(q_hat(kk,:)*conj(q_hat(kk,:)')))
% end
% 
% figure(4)
% plot(MAC)
% MAC = Q_bar*conj(q_hat')/sqrt((Q_bar*conj(Q_bar'))*(q_hat*conj(q_hat')))

% figure(1)
% subplot(211)
% semilogy(omg,abs(G(1:length(omg))))
% title('SUU(z_k)')
% ylabel('Magitude')
% subplot(212)
% plot(omg,angle(G(1:length(omg))))
% ylabel('Phase (rad)')
% xlabel('\omega (Hz)')
% 
% 
% figure(2)
% subplot(211)
% semilogy(omg,abs(SYU(1:length(omg))))
% title('SYU(z_k)')
% ylabel('Magitude')
% subplot(212)
% plot(omg,angle(SYU(1:length(omg))))
% ylabel('Phase (rad)')
% xlabel('\omega (Hz)')


% z = [y,u];
% m = arx(z,[4 4 0])
% m = arx(z,[6 1 0])

% u2 = Sheet2(:,1);
% y2 = Sheet2(:,2);
% yest = sim(m, [u2]);

% figure(1)
% plot(1:70,y2,'--',1:70,yest,'-r')

% load powero.txt
% load windspeed.txt


% Calculation of the SFT of input and output signals
% Len = 2^9
% y2 = powero(1:Len,1);
% u2 = windspeed(1:Len,1);

% % Model 
% A = [1  -1.5  0.7];
% B = [0 1 0.5];
% % Polynominal for System ID
% m0 = idpoly(A,B)
% % Package ID data
% u = iddata([],idinput(300,'rbs'));
% e = iddata([],randn(300,1));
% % Simulate Discrete System
% y = sim(m0, [u e]);