clear all;clc;close all;
%% formation of intial data.
format long
digits(14);
data=xlsread('C:\Users\Unmesh\Documents\Project PMU\Codes\Prony Code, Reference paper and handout\Faulty240forprony.xlsx');
%data=xlsread('C:\Users\Unmesh\Documents\Project PMU\Codes\Prony Code, Reference paper and handout\240forprony.xlsx');
ang=data(:,3); % currents phase
mag=data(:,4); % currents magnitude
%ang=data(:,11); % voltage phase
%mag=data(:,12); % voltage magnitude
freq=data(:,5);
t1=data(:,1);
t2=data(:,2);
t=t1+t2;
omega_t=2*pi*times(freq,t);
f=times(mag,sin(omega_t+ang));
N=100;
n=40;
startPoint=5;
endPoint=n+startPoint;
Y=f(endPoint:N+startPoint);
k=0;
for i=1:1:(N-n+1)
for j=1:1:n
D(i,j)=f(endPoint+k-j);
end
k=k+1;
end
A=
EE=pinv(D)*Y;
eeee=eig(EE);
freqqqq=imag(eeee)./(2*pi*0.2);
dampppp=real(eeee);
[Ub,Sb,Vb]=svd(Ab);
[Ua,Sa,Va]=svd(Aa);
VV=pinv((Vb.')')*(Va.')';
vvvv=eig(VV);
freqqqq1=imag(vvvv);
VVVV=eeee-vvvv;
% a=B'*B;
% b=B'*A;
%c=(A'*A)\(A'*B);%applying LSM for finding cofficent of polynomial equatins
%c=mldivide(A'*A, A'*B);
c=pinv(A'*A)*A'*B;
d=[1;-c];%formation of polymonial equation
u=roots(d);%finding roots of polynomial equation
% formation of matrix for finding residue(Ci)
for i=1:1:N
E(i,:)=f(i);
end
% formation of matrix for finding residue(Ci)
for j=1:1:n
for i=1:1:(N)
F(i,j)=(u(j))^((i-1));
end
end
% ff=F'*F;
% e=F'*E;
C=(F'*F)\(F'*E);%residuse by applying LSM
fn=F*C;%reconstructed signal
yn=fn';
figure(1)
plot(t,f)
title('signal from pmu')
xlabel('time (seconds)')
ylabel('amplitude (per unit)')
%% Reconstructed data
figure(2)
plot(t,yn);%ploting reconstructed data
title('signal reconstructed')
xlabel('time (seconds)')
ylabel('amplitude (per unit)')
er=(fn-f);
%% Ploting error
figure(3)
plot(t,er);%ploting error
title('error');
xlabel('time (seconds)')
ylabel('amplitude (per unit)')
%% calculation of damping freqency phase
l=5*log(u);%continous pole
w=imag(l);
%% damping
damp=real(l);%damping
dampr=real(l)./abs(w);%damping ratio
%% freqeuncy
freq=imag(l)/(2*pi);%freqeuncy
figure(4)
scatter(dampr,freq);
title('frequency and damping of modes included')
xlabel('damping of modes')
ylabel('frequency of modes')
%% Phase
pha=angle(yn);%phase
figure(5)
scatter(t,pha);
title('phase of modes')
xlabel('time (seconds)')
ylabel('phase (radian)')
eng=0;
% U=F';
%% calculation of energy
for i=1:1:n
if angle(l(i,:))==0
m(i,:)=(real(C(i).*F(:,i))).^2;
else
m(i,:)=2.*((real(C(i).*F(:,i)))).^2;
end
end
h=m';
for i=1:1:N
Energy=m(:,i)+eng;
eng=Energy;
end
%% displaying energy,freqency,damping;
for i=1:4
if i==1
U(:,i)=Energy;
end
if i==2
U(:,i)=freq;
end
if i==3
U(:,i)=dampr;
end
end
%display(U)
%% sorting by energy;
sortedEnergy=sortrows(U,1);
display('sorted modes according to energy')
display(sortedEnergy)
% for i=1:1:n
% for j=1:1:N
% En(i)=(real(C(i)*F(j,i)))^2+eng;
% eng=En(i);
% end
% end
% sys=ss(c,0,C,0);
% [kest,L,P] = kalman(sys,1,0.01,0);
评论0