Inter-area oscillation analysis with PMU data

In this post, I present a simple, but effective way of inter-area oscillation analysis with phaser measurement unit (PMU) data. I simply run a Fourier analysis, more specifically a Fast Fourier Transform (FFT), on the measured power system frequency to identify the oscillation modes. So, if you have access to PMU data, you can use the provided code to analyze the oscillation behavior of your power system.

Quick recap: What are inter-area oscillations?

Inter-area oscillations are oscillations where groups of generators in different areas of a power system oscillate against each other.

In the continental European power system, two main oscillation modes have been identified. Namely, a North-South and a West-East oscillation mode with respective frequencies of approximately 0.25 Hz and 0.15 Hz. There is an interesting ENTSO-E report from 2017 about inter-area oscillations in the European power system triggered by the unexpected opening of a line in the French system.

PMU data

The Distributed Electrical Systems Laboratory at EPFL installed five PMUs in their 20-kV medium voltage grid on campus. More info about this monitoring activity can be found here. The best thing is that they provide open records of their PMU measurements with a time resolution of 20 ms for download. So, big thanks for their effort.

For my analysis, I used the data from the 1st – 7th of September 2018. Below the code to automatically download the EPFL PMU data for a specific period. Important: Before downloading data, check if the data for the desired period is available.

import os
import requests

#%% main
#path to store the download files
path_results = r'C:\Your Folder' 

#define the time period for data download (first check online if data for this time period is available)
year = '2018'
month, monthnum = 'Sep', '09' #leading zero
day_from, day_to = 1, 7
PMUid = '2' # avail: 1 - 6, 10 (availability must be checked online, sometimes not all PMUs are available. Check leading zero: before 2020 no leading zero!)

#first part of the url
url_part1 = 'http://nanotera-stg2.epfl.ch/data/'

for ii in range(day_from, day_to + 1): #loop through days
    if ii<10:
        day = '0' + str(ii)
    else:
        day = str(ii)

    for i in range(24): #loop through hours
        if i<10:
            hour = '0' + str(i)
        else:
            hour = str(i)
        
        #put together the file name and url of the data
        name = year + '-' + monthnum + '-' + day  + '_' + hour + 'h_UTC_PMUID' + PMUid + '.txt'
        url = url_part1 + year + '/' + month + '/' + day + '/' + name
        
        #save data
        filesavename = os.path.join(path_results,name)
        r = requests.get(url, allow_redirects=True)
        open(filesavename, 'wb').write(r.content)

Data preprocessing

Before using the data, I preprocessed it. Actually, I just extracted the time and frequency from the measurement files and saved them as pickle file. As the measurement time is given in two columns, one with the second and the other with the millisecond, I combined them into one vector providing the measurement time in seconds. Below the code to do it.

import pandas as pd
import numpy as np
import os
import pickle

#%% main
#path to the downloaded files
path = r'C:\Your Folder'

#name of the result file which will be saved with the processed data
nameres = 'PMUdat_201809_07-01.pkl'

#finds all the files in the folder
files = [x for x in os.listdir(path) if x.endswith(".txt")]
files.sort()

log_time = []
log_freq = []
for i, name in enumerate(files):
    print(name)
    #import measurement file
    dat = pd.read_csv(os.path.join(path,name), sep='\t', decimal='.', header=2, index_col=False)#, nrows=179998)
    
    #only continue if longer than 100 lines (to avoid error if file not found etc.)
    if np.size(dat,0) > 100:
        #remove lines with nans from imported data
        idx = ~(np.isnan(np.array(dat['Frequency [Hz]'])))
        dat = dat.iloc[idx,:]
        
        #combine the time columns to one time vector (in seconds)
        time1 = np.array(dat['SOC [s]'].astype(float))
        time2 = np.array(dat['FRACSEC [ms]'].astype(float))
        time = time1 + time2/1000
        
        #get frequency measurement
        freq = np.array(dat['Frequency [Hz]'])
        
        log_time.append(time)
        log_freq.append(freq)
        
time = np.hstack(log_time)    
freq = np.hstack(log_freq)
#just shift time to start from zero
time = time - time[0]  
#%% save data to pickle in same folder
respathname = os.path.join(path,nameres)
with open(respathname, 'wb') as f:
    pickle.dump([time, freq], f) # time in seconds, frequency in Hz

Inter-area oscillation analysis

In this section, I provide the code and the results of the inter-area oscillation analysis using FFT of the PMU frequency measurements.

Code

Below the code to calculate the FFT from the preprocessed measurement data. Since the FFT requires equidistant measurement samples, I interpolated the gaps in the frequency data according to the time measurements.

import numpy as np
import matplotlib.pyplot as plt
import os
import pickle 

#%% main
#path to the resfile
path = r'C:\Your Folder'
#name of the resfile
nameres = 'PMUdat_201809_07-01.pkl'

#combine path and name
respathname = os.path.join(path,nameres)

#load the preprocessed data
with open(respathname, 'rb') as f:
    in_time, in_freq = pickle.load(f)

Ts = 0.02 # (s) time step of data (20 ms)

#interpolate gaps in the data
time = np.arange(0, max(in_time),Ts)
frequency = np.interp(time, in_time, in_freq)

#calculate fft (omitting negative frequencies)
fft = np.fft.rfft(frequency)
fft_freq = np.fft.rfftfreq(len(time),Ts)

#amplitude of FFT normalized
amp_fft = np.abs(fft) / np.max(np.abs(fft))

#%% plot frequency measurement in time domain
plt.figure()
plt.plot(time, frequency, ls='-', color='k', marker='', label = 'f')
plt.grid(linestyle='--')
plt.xlabel('time (s)')
plt.ylabel('frequency (Hz)')
plt.locator_params(axis='y', nbins=10)
plt.xticks(np.arange(0, 86400*7+1, step = 86400), rotation='vertical')
plt.xlim(0, 86400*7)
plt.title('Frequency measurement \n 1 week - resolution 20 ms')

#%% plot FFT results (0- 0.5 Hz)
plt.figure()
plt.plot(fft_freq, amp_fft, ls='-', marker='', label = 'f')
plt.grid(True)
plt.xticks(np.arange(0, 11, step = 0.05))
plt.xlim(-0.01, 0.5)
plt.ylim(0, 0.000003)
plt.xlabel('oscillation modes (Hz)')
plt.ylabel('amplitude (pu)')
plt.title('FFT (0 - 0.5 Hz)')

Results

According to Nyquist, we can analyze frequency components of a signal up to half of its sample rate. In our case, the time resolution is 20 ms, which corresponds to a sample rate of 50 Hz. Thus, the FFT can provide the frequency components up to 25 Hz. However, I only show the frequency components up to 5 Hz since we want to analyze the inter-area oscillations which are generally in the area of 0.1 – 0.9 Hz.

Let’s first look at the measurement data which is the input for the FFT. Below the time domain plot of the frequency over one week with a time resolution of 20 ms. There are no gaps in the data and the values are in a reasonable range.

Below the results of the FFT calculation up to 5 Hz. There is a lot going on in the frequency range between 0 – 1 Hz (zoomed results in the next plot). However, there are also some higher modes appearing between 1 and 5 Hz. Modes between 2 – 5 Hz are referred to as local modes which come from single or multiple generators oscillating against the system. At this point, I am not able to assign them to a specific source. If you have an idea, please let me know in the comments.

The plot below shows the zoomed results for the frequency range 0 – 0.5 Hz. I have marked in orange the modes, which correspond closest to the North-South (0.25 Hz) and West-East (0.15 Hz) modes mentioned in the ENTSO-E report. However, they do not exactly match those values, here they are approximately 0.14 and 0.28 Hz.

I see two possible reasons for that. Either the system dynamics have changed, or there are some effects related to the measurements and FFT that shift the frequencies slightly. However, I am inclined to say it is the changed system dynamics.

Additionally, there is one mode at around 0.42 Hz and I have absolutely no idea where its coming from.

I would be happy to hear your comments on that.

If you find this article insightful, you might also like the following post: Analysis of frequency measurements of the Nordic, Continental European, and Japanese-50-Hz power systems.

More insights within power systems are to follow in the next article.

Michael