%%======== Reads PMU data in "IEEE .csv" format
% ======== and analyses data on possible deficiencies
% Author: Slava Maslennikov; smaslennikov@iso-ne.com;  March 2018
%
% Instuctions for use
% 1. Assign proper file name for input and output and run this code
% 2. See the results in output file. Pay attentions to messages strting
%    with "*** Error:" and "*** Warning:"
% 3. See generated graphs for sanity check. Big difference between frequency 
%    spectrum for different signals could be a in dicatior of data issue. 
%    Peaks of spectrum indicate the presence of oscillations at corresponding 
%    frequencies. The larger is the peak the larger is the energy of 
%    oscillations and the wider is the spread in the power system.
%       Note: frequency spectrum here could be havily impacted by "bad" PMU
%    data. This code minimizes the impact of "bad" PMU data only given as zero
%    values in input.
%
% This code performs the following tests and verifications:
% 1. Types and dimensions of PMU measurements listed in header.
% 2. Types of transmission elements with measurements.
% 3. Cosistency of SignalName (first raw in header) structure to
%    recommendations.
% 4. Consistency of time stamps and interval between measurements.
% 5. Basic sanity check on average magnitude values for frequency, current
%    and voltage relative to reasonability range.
% 6. Count the number of measurements by type.
% 7. Verification of a possibility to calculate MW and MVAR of all transmission
%    elements where current is given in the input set.
% 8. Verification on whether the input includes a comprehensive set of
%    measurements for all transmission elements having current measurements.
% 9. Plots graphs of magnitudes for voltages, currents and frequencies.
% 10. Unwraps current and voltage angles and plots relative values of angles
%    with corresponding average value per measurement frame as a reference.
% 11. Calculates and plots frequency spectra for frequency, voltage and
%    current magnitudes for normolized values of oscillation emergy. 
%    This plot has a cutoff frequency to remove DC bias.
%%
clear all

%------- Input and Output file names
InputFile='ISO-NE_case1.csv';
OutFile='Report.txt';

%--- Constants ---- do not modify them if not necessary
typeAllowed={'F','VM','VA','IM','IA'};  %!!! do not change the sequence
dimAllowed={'HZ','KV','DEG','AMP','DEG'}; %!!! do not change the seqience
typeAllowedE={'Bus','Ln','Gen','Load','HVDC','FACTS'};  % allowed types of transmission elements
separatorAllowed={':'};
firstColumn={'Time';'T';'sec'}; %!!! do not change the seqience
thresholdIm=[0 3000]; % min and max AMP threshold for current sanity check
thresholdVm=[0.2 900]; % min and max KV threshold for voltage sanity check
thresholdF=1.0; % HZ threshold for frequency sanity check
F_cutoff=0.1; % Hz  low-band cutoff frequency for spectra plotting
%--------------------------------------------------------

disp('Start of data analysis')
tic

fid = fopen(OutFile,'wt'); % Open outut file
if exist(InputFile, 'file')~=2
    % PMU file does not exist
    fprintf(fid, ' *** Error: PMU file does not exist:   %s \n',InputFile);
    fprintf(' *** Error: PMU file does not exist:   %s \n',InputFile);
    fclose(fid) ;    
end
% PMU file exists
fprintf(fid, '--- Loading data from file: %s  \n\n',InputFile);
 
%--- Load PMU data   
[num,txt] = xlsread(InputFile);
    
%--- If last column(s) of data contain all NaN then Matlab does not read them
%   need to add columns with NaNs
[~,col]=size(txt);
[samples,j]=size(num);
if col>j; num(:,j+1:col)=NaN; end
       
%--- Assign specific lines of title and time vector
LineFull=txt(1, 2:end)';
typeS=txt(2, 2:end)';
dimS=txt(3, 2:end)';
LineShort=txt(4, 2:end)';
Time=num(:,1); % time
col=col-1; % the number of PMU data columns

Tstart=txt{4,1}; % Start date/time of PMU data period

%--- Verify the syntax of the first column header
if ~cellfun(@strcmp, txt(1,1),firstColumn(1))
    fprintf(fid, '*** Warning: First comumn header is   %s   instread of   %s\n',txt{1,1},firstColumn{1});
end
if ~cellfun(@strcmp, txt(2,1),firstColumn(2))
    fprintf(fid, '*** Warning: First comumn header is   %s   instread of   %s\n',txt{2,1},firstColumn{2});
end
if ~cellfun(@strcmp, txt(3,1),firstColumn(3))
    fprintf(fid, '*** Warning: First comumn header is   %s   instread of   %s\n',txt{3,1},firstColumn{3});
end

% -- Verify the dimension of measurements in header
[~,temp]=ismember(dimS,dimAllowed);
for i=1:length(dimS)
    if temp(i)==0
        fprintf(fid, '*** Warning: not allowed dimension of measurements  %s  in column %4.0f \n',dimS{i},(i+1));
    end 
end

%--- Verify the type of measurements in header
[~,temp]=ismember(typeS,typeAllowed);
for i=1:length(typeS)
    if temp(i)==0
        fprintf(fid, '*** Warning: not allowed type of measurements  %s  in column %4.0f \n',typeS{i},(i+1));
    end 
end

% -- Verify the match of type and dimension of measurements
for i=1:col
    for j=1:length(dimAllowed)
      if cellfun(@strcmp, typeAllowed(j),typeS(i))
         if ~cellfun(@strcmp, dimAllowed(j),dimS(i))
           fprintf(fid, '*** Error: not allowed combination of type  %s  and dimension  %s of measurements in column %4.0f \n',typeS{i},dimS{i},(i+1));
         end 
      end
    end
end

%--- Verify time vector
for i=1:samples-1; dT(i)=Time(i+1)-Time(i); end % time increments between samples
fs=floor(1/median(dT)+0.4);
if fs <30
    fprintf(fid, '*** Error: sampling rate is less than 30 fps: %4.0f \n',fs);
end
for i=1:samples-1
    if dT(i)>1/fs*1.1  %   10% tolerance
        fprintf(fid, '*** Error: time interval at sample %5.0f  (%4.4f) exceeds time for sampling rate %3.4f\n',i,abs(dT(i)),1/fs);
    end
end
for i=1:samples-1
    if dT(i)<0 
        fprintf(fid, '*** Error: time interval at sample %5.0f is negative: %5.3f\n',i,dT(i));
    end
end

%--- Verify the structure of SignalName (first line of header) per the number of separators
for i=1:col
    ind = find(LineFull{i,1}==separatorAllowed{1,1});
    if length(ind)<3
        fprintf(fid, '*** Warning: not recommended SignalName structure in the first line of header: %s   column %5.0f \n',LineFull{i,1},i);
    else
        % Verify the type of transmission element
        if ~ismember(typeAllowedE,{LineFull{i,1}(ind(2)+1:ind(3)-1)})
           fprintf(fid, '*** Warning: non-standard type of transmission element measurements  %s  in column %4.0f \n',LineFull{i,1}(ind(2)+1:ind(3)-1),(i+1));
        end
    end
end
 
F=[]; Vm=[]; Va=[]; Im=[]; Ia=[];
indF=[]; indVm=[]; indVa=[]; indIm=[]; indIa=[];
%--- Sepatate signals into corresponding matrices
for i=1:col
   switch temp(i)
        case 1
          F=[F num(:,i+1)]; indF=[indF i+1];
        case 2
          Vm=[Vm num(:,i+1)]; indVm=[indVm i+1];
        case 3
          Va=[Va num(:,i+1)]; indVa=[indVa i+1];
        case 4
          Im=[Im num(:,i+1)]; indIm=[indIm i+1];
        case 5
          Ia=[Ia num(:,i+1)]; indIa=[indIa i+1];
   end
end
Nsig=zeros(1,5); % number of signals 1-F, 2-Vm, 3-Va, 4-Im,5-Ia
if ~isempty(F); Nsig(1)=length(indF); end
if ~isempty(Vm); Nsig(2)=length(indVm); end
if ~isempty(Va); Nsig(3)=length(indVa); end
if ~isempty(Im); Nsig(4)=length(indIm); end
if ~isempty(Ia); Nsig(5)=length(indIa); end
fprintf(fid, 'Input data set consists of: \n');
fprintf(fid, ' %5.0f  frequency signals \n',Nsig(1));
fprintf(fid, ' %5.0f  voltage magnitude signals \n',Nsig(2));
fprintf(fid, ' %5.0f  voltage angle signals \n',Nsig(3));
fprintf(fid, ' %5.0f  current magnitude signals \n',Nsig(4));
fprintf(fid, ' %5.0f  current angle signals \n',Nsig(5));
fprintf(fid, '\n');

%--- Check the consistency of I and V phasors by the number of signals
if Nsig(2)~= Nsig(3)
    fprintf(fid, '\n*** Warning: the number of voltage magnitude and phase signals is not the same\n');
end
if Nsig(4)~= Nsig(5) 
    fprintf(fid, '\n*** Warning: the number of current magnitude and phase signals is not the same\n');
end

%--- Verify an ability to calculate MW and Mvar flow for lines with given current
%=== Count lines with I and fill in LineID for thse lines and Substation IDs for all signals
if Nsig(4)==0
    fprintf(fid, '*** Warning: data set does not contain current signals\n');
else
   for i=1:col
       buf1=find(LineFull{i,1}==separatorAllowed{1,1});
       Sub{i,1}=LineFull{i,1}(buf1(1)+1:buf1(2)-1); % name of Substation 
   end
%=== Filll in column indices of input signal for lines with current
%    ind: 1-Im 2-Ia 3-f 4-Vm 5-Va
   ind=zeros(Nsig(4),5);
   ind(:,1)=find(temp==4); % indices of Im
   LineID(:,1)=LineFull(ind(:,1));

   for j=1:Nsig(4) % loop per the number of lines with current
      for i=1:col
        if cellfun(@strcmp, LineFull(i),LineID(j)) & cellfun(@strcmp, typeS(i),typeAllowed(5)); ind(j,2)=i; end % Ia
        if cellfun(@strcmp, LineFull(i),LineID(j)) & cellfun(@strcmp, typeS(i),typeAllowed(1)); ind(j,3)=i; end % F
        if cellfun(@strcmp, LineFull(i),LineID(j)) & cellfun(@strcmp, typeS(i),typeAllowed(2)); ind(j,4)=i; end % Vm
        if cellfun(@strcmp, LineFull(i),LineID(j)) & cellfun(@strcmp, typeS(i),typeAllowed(3)); ind(j,5)=i; end % Va
      end
   end
% For a current signal, fiil in indices for F, Vm, Va from the same substation if F or V is not in input
% Use the quantity at the same substation listed first in input
   for j=1:Nsig(4)
      if ind(j,3)==0 % Line with current in input data does not have related frequency
          temp1=find(ismember(Sub,Sub(ind(j,1)))& temp==1);
          if ~isempty(temp1)
            ind(j,3)=temp1(1);
            fprintf(fid, '> Frequency for line  %s  can be linked to substation  %s (column %4.0f) \n',LineID{j},Sub{ind(j,3)},(ind(j,3)+1));
          end
      end
      if ind(j,4)==0 % Line with current in input data does not have related voltage magnitude
          temp1=find(ismember(Sub,Sub(ind(j,1)))& temp==2);
          if ~isempty(temp1)
            ind(j,4)=temp1(1);
            fprintf(fid, '> Voltage magnitude for line  %s   can be linked to substation  %s (column %4.0f) \n',LineID{j},Sub{ind(j,4)},(ind(j,4)+1));
          end
      end
      if ind(j,5)==0 % Line with current in input data does not have related voltage angle
          temp1=find(ismember(Sub,Sub(ind(j,1)))& temp==3);
          if ~isempty(temp1)
            ind(j,5)=temp1(1);
            fprintf(fid, '> Voltage angle for line  %s   can be linked to substation  %s (column %4.0f) \n',LineID{j},Sub{ind(j,5)},(ind(j,5)+1));
          end
      end        
   end
   
   %--- Check the completness of the data set for each line with current
   %    A complete set includes Im Ia F Vm Va
   temp1=sum(ismember(ind,0)');
   Nll=sum(ismember(temp1,0)'); % number of lines with full set of data
   if Nll>0
       fprintf(fid, '\n--There are %4.0f transmission elements with full data set(f,Vm,Va,Im,Ia) \n',Nll);
       for i=1:Nsig(4)
           if temp1(i)==0;fprintf(fid, '%s \n',LineID{i}); end
       end
   else
       fprintf(fid, '\n*** Warning: there are no transmission elements with full data set(f,Vm,Va,Im,Ia) \n');
   end
   fprintf(fid, '\n');
   for i=1:Nsig(4)
      if ind(i,2)==0; fprintf(fid, ' *** Warning: current angle for line  %s  is not in data set \n',LineID{i});  end
      if ind(i,3)==0; fprintf(fid, ' *** Warning: frequency for line  %s  is not in data set \n',LineID{i});  end
      if ind(i,4)==0; fprintf(fid, ' *** Warning: voltage magnitude for line  %s  is not in data set \n',LineID{i}); end
      if ind(i,5)==0; fprintf(fid, ' *** Warning: voltage angle for line  %s  is not in data set \n',LineID{i}); end      
   end
end

%--- Estimate frequency of the system
if Nsig(1)>0
   Fmed= median(F,'omitnan');
   Fmm=median(Fmed,'omitnan');
   Fsys=0;
   if Fmm >50-thresholdF & Fmm <50+thresholdF;    Fsys=50.0; end
   if Fmm >60-thresholdF & Fmm <60+thresholdF;    Fsys=60.0; end
   if Fsys==0; fprintf(fid, '*** Error: system frequency is not 50 nor 60 Hz\n'); end

   %--- Check the median frequency  per threshold
   for i=1:Nsig(1)
      if Fmed(i)~=NaN & (Fmed(i)<Fsys-thresholdF | Fmed(i)>Fsys+thresholdF)
        fprintf(fid, '*** Warning: median frequency in column  %5.0f is out of range; median is: %4.2f\n',indF(i),Fmed(i));
      end
   end
end

%--- Check the median of current magnitude  per threshold
if Nsig(4)>0
   for i=1:length(Im(1,:))
       temp=mean(Im(:,i),'omitnan');
       if temp<thresholdIm(1) | temp>thresholdIm(2)
          fprintf(fid, '*** Warning: Avarage magnitude of current in column  %5.0f is out of range; average value is: %5.2f\n',indIm(i),temp);
       end
   end
end

%--- Check the median of voltage magnitude  per threshold
if Nsig(2)>0
   for i=1:length(Vm(1,:))
       temp=mean(Vm(:,i),'omitnan');
       if temp<thresholdVm(1) | temp>thresholdVm(2)
          fprintf(fid, '*** Warning: Avarage magnitude of voltage in column  %5.0f is out of range; average value is: %5.2f\n',indVm(i),temp);
       end
   end
end

%--- Unwrap angles and remove average per measurement frame
% Unwarap angles
if Nsig(3)>0
   Va_un = unwrap(Va(:,1:end)/180*pi)*180/pi; % unwrapping must be done in radians but not in degrees
   Va_relative=Va_un-mean(Va_un','omitnan')';
end
if Nsig(5)>0
   Ia_un = unwrap(Ia(:,1:end)/180*pi)*180/pi;
   Ia_relative=Ia_un-mean(Ia_un','omitnan')';
end

%--- Estimate spectrum of oscillations from Vm and Im
f = fs/2*linspace(0,1,samples/2+1)'; %frequency

if Nsig(4)>0
   xx=Im;
   for j=1:samples; if abs(xx(j))<0.01; xx(j)=NaN; end; end % replace zeros by NaN
   xx=xx - ones(samples,1)*mean(xx,'omitnan'); % remove bias
   for i=1:length(xx(1,:)) % replace NaN by 1e-9
      for j=1:samples; if isnan(xx(j,i));  xx(j,i)=1e-9; end;  end
   end
   if Nsig(4)==1
       sIm =(abs(2*fft(xx, samples)/samples)')'; % magnitude of Im spectrum; sum of all transmission elements
   else
     sIm = sum(abs(2*fft(xx, samples)/samples)')'; % magnitude of Im spectrum; sum of all transmission elements
   end
end
if Nsig(3)>0
   xx=Vm;
   for j=1:samples; if abs(xx(j))<0.01; xx(j)=NaN; end; end % replace zeros by NaN
   xx=xx - ones(samples,1)*mean(xx,'omitnan');
   for i=1:length(xx(1,:))
      for j=1:samples;  if isnan(xx(j,i));  xx(j,i)=1e-9; end;  end
   end
   if Nsig(3)==1
     sVm =(abs(2*fft(xx, samples)/samples)')'; % magnitude of Vm spectrum; sum of all transmission elements
   else
     sVm = sum(abs(2*fft(xx, samples)/samples)')'; % magnitude of Vm spectrum; sum of all transmission elements
   end
end

%--- Estimate spectrum of oscillations from F
if Nsig(1)>0
   xx=F;
   for j=1:samples; if abs(xx(j))<0.01; xx(j)=NaN; end; end % replace zeros by NaN
   xx=xx - ones(samples,1)*mean(xx,'omitnan'); % remove bias
   for i=1:length(xx(1,:)) % replace NaN by 1e-9
      for j=1:samples; if isnan(xx(j,i));  xx(j,i)=1e-9; end;  end
   end
   if Nsig(1)==1
       sF =(abs(2*fft(xx, samples)/samples)')'; % magnitude of F spectrum; sum of all transmission elements
   else
     sF = sum(abs(2*fft(xx, samples)/samples)')'; % magnitude of F spectrum; sum of all transmission elements
   end
end



%--- Plot: current, voltage, frequency, angles and spectrum
if Nsig(4)>0
  figure (1) % Current magnitude
  plot(Time,Im)
  xlabel('Time, s','FontSize',10)
  ylabel('AMP','FontSize',10)
  title('Current magnitudes','FontSize',11)
end

if Nsig(2)>0
  figure (2) % Voltage magnitude
  plot(Time,Vm)
  xlabel('Time, s','FontSize',10)
  ylabel('KV','FontSize',10)
  title('Voltage magnitudes; phase-to-ground','FontSize',11)
end

if Nsig(1)>0
  figure (3) % Frequency
  plot(Time,F)
  xlabel('Time, s','FontSize',10)
  ylabel('Hz','FontSize',10)
  title('Frequency','FontSize',11)
  ylim([(Fsys-thresholdF) (Fsys+thresholdF)])
end

if Nsig(5)>0
  figure (4) % Current angle
  plot(Time,Ia_relative)
  xlabel('Time, s','FontSize',10)
  ylabel('Degrees','FontSize',10)
  title('Angle of currents; relative to average in measurement frame','FontSize',11)
end

if Nsig(3)>0
  figure (5)  % Voltage angle
  plot(Time,Va_relative)
  xlabel('Time, s','FontSize',10)
  ylabel('Degrees','FontSize',10)
  title('Angle of voltages; relative to average in measurement frame','FontSize',11)
end

if Nsig(1)>0 | Nsig(2)>0 | Nsig(4)>0 
   [~,ind2]=sort(abs(f-F_cutoff),'ascend'); % ind2(1) index of F_cutoff
   if Nsig(4)>0; sIm=sIm/max(sIm(ind2(1):length(f))); end % normolize to have 1 p.u. at F_cutoff
   if Nsig(2)>0; sVm=sVm/max(sVm(ind2(1):length(f))); end
   if Nsig(1)>0; sF=sF/max(sF(ind2(1):length(f))); end

   figure(6)  % Spectra of Im, Vm, F
   le={};
   if Nsig(4)>0
     plot(f(ind2(1):end),sIm(ind2(1):length(f)),'b')
     le=[le 'Current'];
     hold on
   end
   if Nsig(2)>0
      plot(f(ind2(1):end),sVm(ind2(1):length(f)),'r')
      le=[le 'Voltage'];
      hold on
   end
   if Nsig(1)>0
      plot(f(ind2(1):end),sF(ind2(1):length(f)),'k')
      le=[le 'Frequency'];
      hold on
   end
   hold off
   xlabel('Frequency,Hz','FontSize',10)
   ylabel('PSD','FontSize',10)
   h=legend(le);set(h,'Location','NorthEast','fontsize',8);
   title('Spectrum of oscillations from magnitudes','FontSize',11)
end

disp('Competed data analysis. See Report file for details')
toc

fclose all;
