Project_Example_09.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sean Semper                                         System Identification
%/ Sumit Tripathi
% Project Code 1                                                    MAE 566
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Created Spring 2009                                                     %
% Last Modified: 05/7/09                                                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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                                                 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% To Pick Relationships Set the Following 
% in_num = #;               % Input Index
% out_num = #;              % Output Index
% load min_avg_####.txt     % Set the dataset March 2009, Oct and Sept 2008
% data = min_avg_####;
% Then run the code...And follow Prompts
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%time%	nobs(1)	avg_rpm(2) 	avg_powero(3)	avg_powerr(4)  avg_windspeed(5)      
load min_avg_0908.txt
data = min_avg_0908;
display('Loading')
numHnkl = 100;  % Hankel Matrix with n, Markov Parameters
mm = 1;         % Number of States to Keep
in_num = 2;     % Input Index
out_num = 3;    % Output Index
in_out_ind = [in_num, out_num];
u = data(:,in_num);
y = data(:,out_num);
clear data


L = length(u);
dfftu = fft(u)/L;
dffty = fft(y)/L;
dt = 1;
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 = [];
display('Computing')
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(1)
stem(diag(SS));
diag(SS)
ylabel('Eigenvalue ( \sigma_i ) Magnitude')
xlabel('\sigma_i')
title({'Using the Frequency Response';'(Considering only \sigma_1, to \sigma_n)'})


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

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

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);

[EigVec, EigA] = eig(A_est);
diag(EigA)

q_bar = inv(EigVec)*sqrt(SS)*S';

b_hat = inv(EigVec)*B_est;

for kp = 1:1:mm
kk=0;
for k=1:mm
q_hat(kp,k) = [EigA(kp,kp)^kk*b_hat(kp)] ;
kk=kk+1;
end
end

k=1;
for k=1:mm
MAC(k,1) = abs(q_hat(k,1:mm)*conj(q_bar(k,1:mm))')/ ...
           sqrt(abs((q_bar(k,1:mm)*conj(q_bar(k,1:mm))'))*abs((q_hat(k,1:mm)*conj(q_hat(k,1:mm))')));
end

figure(2)
plot(1:mm,MAC,'o-')
ylim([0 1.2])
title('Modal Amplitude Coherence (MAC) ')
ylabel('MAC Value')
xlabel('Eigenvalue Index \it i')
set(gca,'XTick',1:1:4)
set(gca,'XTickLabel',{'1','2','3','4'})


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Logic to Continue on to Simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display('The code is paused for output review...')
pause(1)
reply = input('Please hit Enter Key to continue...Or Type "n" and Hit "Enter" to break:-->   ', 's');
if isempty(reply)
    Project_Example2_09
else
    break    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%