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.
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)
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 #%% 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.
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)')
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.