pronycoed.rar_Prony code_pmu_pmu data_prony

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