PMU Data Event Detection User Guide

· 40 pages
Find parts for this model

Manual text content

NREL is a national laboratory of the U.S. Department of Energy Office of Energy Efficiency & Renewable Energy Operated by the Alliance for Sustainable Energy, LLC This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. Contract No. DE-AC36-08GO28308 PMU Data Event Detection: A User Guide for Power Engineers A. Allen, M. Singh, and E. Muljadi National Renewable Energy Laboratory S. Santoso The University of Texas at Austin Link to Accompanying Zipped Data File (6.3 MB) Technical Report NREL/TP-5D00-61664 October 2014 NREL is a national laboratory of the U.S. Department of Energy Office of Energy Efficiency & Renewable Energy Operated by the Alliance for Sustainable Energy, LLC This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. Contract No. DE-AC36-08GO28308 National Renewable Energy Laboratory 15013 Denver West Parkway Golden, CO 80401 303-275- 3000 • www.nrel.gov PMU Data Event Detection: A User Guide for Power Engineers A. Allen, M. Singh, and E. Muljadi National Renewable Energy Laboratory S. Santoso The University of Texas at Austin Prepared under Task No. WE14.9A01 Technical Report NREL/TP-5D00-61664 October 2014 NOTICE This report was prepared as an account of work sponsored by an agency of the United States government. Neither the United States government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States government or any agency thereof. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. Available electronically at http://www.osti.gov/scitech Available for a processing fee to U.S. Department of Energy and its contractors, in paper, from: U.S. Department of Energy Office of Scientific and Technical Information P.O. Box 62 Oak Ridge, TN 37831-0062 phone: 865.576.8401 fax: 865.576.5728 email: mailto:reports@adonis.osti.gov Available for sale to the public, in paper, from: U.S. Department of Commerce National Technical Information Service 5285 Port Royal Road Springfield, VA 22161 phone: 800.553.6847 fax: 703.605.6900 email: orders@ntis.fedworld.gov online ordering: http://www.ntis.gov/help/ordermethods.aspx Cover Photos: (left to right) photo by Pat Corkery, NREL 16416, photo from SunEdison, NREL 17423, photo by Pat Corkery, NREL 16560, photo by Dennis Schroeder, NREL 17613, photo by Dean Armstrong, NREL 17436, photo by Pat Corkery, NREL 17721. NREL prints on paper that contains recycled content. iii List of Acronyms CSV comma-separated value ERCOT Electric Reliability Council of Texas FFT fast Fourier transform MATLAB Matrix Laboratory PMU phasor measurement unit PSD power spectral density RPAD relative phase angle difference UT University of Texas This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. iv Executive Summary This user guide is intended to accompany a software package containing a Matrix Laboratory (MATLAB) script and related functions for processing phasor measurement unit (PMU) data. This package and guide have been developed by the National Renewable Energy Laboratory and the University of Texas at Austin. The objective of this data processing exercise is to discover events in the vast quantities of data collected by PMUs. This document attempts to cover some of the theory behind processing the data to isolate events as well as the functioning of the MATLAB scripts. This user guide is divided into two parts. The first part describes the algorithms and mathematical background that the accompanying MATLAB codes use to detect events in PMU data. The second part describes the inputs required from the user and the outputs generated by the scripts. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. v Table of Contents List of Figures ......................................................................................................................................... vii List of Tables .......................................................................................................................................... viii PART 1 ....................................................................................................................................................... 1 Motivation .................................................................................................................................................. 1 Introduction to Signal-Processing Methods for the Analysis ................................................................ 1 Fast Fourier Transform .................................................................................................................. 1 Matrix-Pencil ................................................................................................................................. 3 Yule-Walker Spectral ..................................................................................................................... 5 Min-Max ........................................................................................................................................ 6 Event Detection.............................................................................................................................. 7 PMU Data Formats .................................................................................................................................... 8 PMU Data Description ......................................................................................................................... 8 Present Data Format ............................................................................................................................. 8 Protocol to Read in the PMU Data ....................................................................................................... 9 Preprocessing the Data .......................................................................................................................... 10 Unwrapping Voltage Phase Angle ...................................................................................................... 10 Angle Differences .............................................................................................................................. 11 Detection of Dropped Data Points ...................................................................................................... 11 Data Window Length ......................................................................................................................... 12 Detection of Events ................................................................................................................................ 12 Algorithm Description ........................................................................................................................ 12 Categorization of Events ........................................................................................................................ 13 PART 2 ..................................................................................................................................................... 14 Program Inputs ....................................................................................................................................... 14 Program Steps and Outputs .................................................................................................................. 15 Conclusion .............................................................................................................................................. 21 References............................................................................................................................................... 21 Appendix: Help and Documentation for Functions ............................................................................. 23 fReadPMUFile ................................................................................................................................... 23 Syntax .......................................................................................................................................... 23 Description................................................................................................................................... 23 Related Functions ........................................................................................................................ 23 Example ....................................................................................................................................... 23 fPMU_Protocol .................................................................................................................................. 24 Syntax .......................................................................................................................................... 24 Description................................................................................................................................... 24 Related Functions ........................................................................................................................ 24 Example ....................................................................................................................................... 25 fpeakdet .............................................................................................................................................. 25 Example ....................................................................................................................................... 25 fMatrixPencil ...................................................................................................................................... 26 Syntax .......................................................................................................................................... 26 Description................................................................................................................................... 26 Examples ..................................................................................................................................... 27 fEventCount ....................................................................................................................................... 28 Syntax .......................................................................................................................................... 28 Description................................................................................................................................... 28 Input ............................................................................................................................................. 28 Output .......................................................................................................................................... 28 fPlotandAnalysis ................................................................................................................................ 28 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. vi Syntax .......................................................................................................................................... 28 Description................................................................................................................................... 29 Input ............................................................................................................................................. 29 Output .......................................................................................................................................... 29 Example ....................................................................................................................................... 29 fAngleProcess_singlephase_unwrap .................................................................................................. 30 Syntax .......................................................................................................................................... 30 Description:.................................................................................................................................. 30 Example ....................................................................................................................................... 30 Reference ........................................................................................................................................... 31 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. vii List of Figures Figure 1. Window of voltage relative phase angle difference (RPAD) .................................................. 2 Figure 2. Amplitude of spectrum of voltage RPAD ................................................................................ 2 Figure 3. In the FFT method, data windows during which the peak magnitude exceeds three times the standard deviation of peak magnitudes for the entire hour are marked with an X. ............... 3 Figure 4. An example of the matrix-pencil method used to estimate parameters of the de-noised signal ................................................................................................................................................... 5 Figure 5. In the matrix-pencil method, data windows during which the peak magnitude exceeds three times the standard deviation of peak magnitudes for the entire hour are tagged. ............. 5 Figure 6. Yule-Walker method power of frequency content in a PMU signal during a large disturbance ......................................................................................................................................... 6 Figure 7. In the Yule-Walker method, data windows during which the peak magnitude exceeds three times the standard deviation are tagged. ............................................................................... 6 Figure 8. In the min-max method, data windows during which the difference between the minimum and maximum values in a 10-second data window exceeds three times the standard deviation are tagged. .......................................................................................................................................... 7 Figure 9. Comparison of all of the screening algorithms ...................................................................... 7 Figure 10. An entire day of PMU frequency from UT-Austin ................................................................. 8 Figure 11. (Left) Equipment for each PMU station and (right) the phasor data concentrator and PC within the Independent Texas Synchrophasor Network ................................................................. 9 Figure 12. Raw voltage phase angle for two PMU stations ................................................................. 11 Figure 13. RPAD between the same two PMU stations corrected by individually unwrapping and taking the difference ........................................................................................................................ 11 Figure 14. (Left) An illustration of a moving window applied to PMU data and (right) a 10-second data window for analysis ................................................................................................................. 12 Figure 15 The output plot of the first signal listed in the variable angle_difference ........................ 16 Figure 16. One window of data (10 seconds/300 data points) ............................................................ 16 Figure 17. FFT of window with peak detected ...................................................................................... 17 Figure 18. Peak of FFT for every 10-second window in the data, with the (black line) mean and (red line) 3.5 standard deviations from the mean delineated ............................................................... 17 Figure 19. Screening of suspect events to find events detected by three or more methods. (Such events are marked with an X.) ......................................................................................................... 18 Figure 20. An example of an RPAD single (large) impulse event (RPAD between UT-Austin and McDonald) ......................................................................................................................................... 18 Figure 21. An example of an RPAD damped transient Event 1 (RPAD between UT-Austin and Waco)................................................................................................................................................. 19 Figure 22. An example of an RPAD sustained transient (RPAD between UT-Austin and McDonald)19 Figure 23. A contour plot of the PSD (low-frequency content over time) for a voltage RPAD sustained transient event (RPAD between UT-Austin and McDonald) ........................................ 20 Figure 24. An example of an RPAD sustained step change event (RPAD between McDonald and Waco)................................................................................................................................................. 20 Figure 25. An example of an RPAD momentary step change event (RPAD between UT-Austin and Waco)................................................................................................................................................. 21 Figure 26. Example plot from function fReadPMUFile ......................................................................... 24 Figure 27. Signal with peaks detected .................................................................................................. 26 Figure 28. Sigma values for the matrix-pencil method ........................................................................ 27 Figure 29. Comparison of original and reconstituted signals............................................................. 28 Figure 30. Example signal ...................................................................................................................... 29 Figure 31. Peak of FFT for every 10-second window, with the (black line) mean and (red line) 3.5 standard deviation delineated ......................................................................................................... 30 Figure 32. Raw angle data with jumps and discontinuities present ................................................... 31 Figure 33. Processed angle data with discontinuities removed ......................................................... 31 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. viii List of Tables Table 1. An example of Type 1 PMU Data in CSV Format for Two PMU Stations With the Magnitude, Angle, and Frequency for Each Station Grouped Together............................................................ 9 Table 2. An Example of Type 2 PMU Data in CSV Format for Three PMU Stations With the Frequencies for Each Station Listed in the Last Columns ........................................................... 10 Table 3. An Example of the Protocol Function Output ........................................................................ 10 Table 4. Summary of Characteristics From Events Detected in RPAD Data ..................................... 13 Table 5. Summary of Characteristics from Events Detected in Frequency Data .............................. 14 Table 6. Output Protocol ........................................................................................................................ 25 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 1 PART 1 Motivation This guide describes how to analyze data from phasor measurement units (PMUs). It is based on data from the Texas Synchrophasor Network. 1 The Texas Synchrophasor Network has been continuously generating data since it was set up in the fall of 2008. Because so much data are being generated, it is difficult to detect events of interest and analyze them for power system studies. Resources such as the Electric Reliability Council of Texas’s (ERCOT’s) Daily Grid Operations Reports provide information on sudden losses of power above 450 MW caused by generating unit trips and information on line contingencies (ERCOT 2004–2010). (Total peak generating capacity in ERCOT is above 65 GW.) In addition, it is straightforward to visually detect large, sudden imbalances in generation and load by monitoring frequency; however, information related to events such as transmission line reclosing and trips and other equipment trips is not readily available. These system events are visible in the PMU data and this guide provides instructions on how to automatically detect these events of interest. In this guide, the voltage phase angle difference between two PMU locations is used to detect power system events that are visible in the voltage phase angle signal. However, when a large number of PMUs are installed in a system, a common reference phase angle for all PMUs is desired. This reference phase angle can be determined by using a center-of-gravity concept applied to PMU-derived local frequencies. A common phase angle reference can be subsequently obtained from this center-of-gravity frequency. Any PMUs located within the interconnected system can be referenced to this center-of-gravity derived reference frequency or reference phase angle. Many power system measurements are available to estimate the center of gravity, as is well known in the literature (Kundur 1994) and will not be discussed here. Introduction to Signal-Processing Methods for the Analysis The following signal-processing methods are used for the analysis of the PMU data. Familiarity with the mathematical underpinnings of these methods is necessary before examining the codes that will be used for the analysis. These analyses are applied to one window of data at a time. (More details about the data window are provided later.) It should be noted that every sample in the data is analyzed. The window is not used to smooth the data; instead, it is used to section the data into manageable lengths. Fast Fourier Transform Before the method is applied, the data window is differentiated using the Matrix Laboratory (MATLAB) function diff to remove the direct current (dc) offset in the data to facilitate the detection of events in the low-frequency range below 1 Hz. Differentiating the data increases the noise in the data; however, it has been found that the magnitude of frequencies associated with the events in the data tends to be very high. Differentiation also does not alter the frequency content of the signal. Therefore, increased noise is not a concern for this analysis. A digital filter can be applied as an alternative by the user. Leakage is also not a concern because the frequency content itself is not of as much importance as the relative magnitudes of the frequencies. If leakage is a concern, the user can apply a Hanning window or other discrete window function. 1 http://web.ecs.baylor.edu/faculty/grady/Texas_Synchrophasor_Network.html This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 2 The fast Fourier transform (FFT) utilizes the MATLAB function fft to detect events in the PMU data. This function returns the discrete Fourier transform calculated using the FFT algorithm. The discrete Fourier transform is used to find the strongest frequency component within the data window. The maximum amplitude, | 𝑌 ( 𝑓 )| , is saved for frequencies in the range of 0 < 𝑓 ≤ 2 𝐻𝑧 . The method is applied to a data window of the voltage phase angle in Figure 1 containing oscillations induced by a large event. The discrete Fourier transform applied to the data window containing event information is shown in Figure 2 . The Fourier transform in Figure 2 shows a strong 0.64-Hz component. In the algorithm, the peak magnitude, | 𝑌 ( 𝑓 )| = 0.03 , is saved for this data window. Figure 1. Window of voltage relative phase angle difference (RPAD) Figure 2. Amplitude of spectrum of voltage RPAD The FFT method is applied to each 10-second data window for the entire hour of PMU data. The averages and standard deviations are calculated for all saved peak magnitudes. Any saved magnitudes above three standard deviations are tagged as possible events. A three-standard- deviations window was selected based on examples of detected events in the Texas Synchrophasor Network data. A plot of the peak magnitudes from each 10-second data window 1928 1930 1932 1934 1936 1938 1940 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 Data Window During a Large Disturbance Time, seconds Voltage Phase Angle, degrees 0 1 2 3 4 5 0 0.005 0.01 0.015 0.02 0.025 0.03 Single-Sided Amplitude Spectrum of y(t) Frequency (Hz) |Y(f)| This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 3 is shown for an hour of PMU frequency in Figure 3 . The averages and standard deviations of the peak magnitudes are indicated. The outliers, or possible events, are also marked in Figure 3 . Figure 3. In the FFT method, data windows during which the peak magnitude exceeds three times the standard deviation of peak magnitudes for the entire hour are marked with an X. Matrix-Pencil The matrix-pencil method (Liu, Quintero, and Venkatasubramanian 2007) fits a sum of damped sinusoids to evenly sampled PMU data. The amplitude, phase angle, frequency, and damping are the parameters of the damped sinusoids estimated to fit the PMU data. The mathematical representation of the sum of damped sinusoid estimates fitted to the sampled data is given in ( 1 ). ( ) ∑ = = n i k i i z R k y 1 (1) where the discrete signal ( ) k y in (3) is equal to the sum of the product of the residues or complex amplitudes R and the poles z . The parameter n is the number of sinusoids or modes to be estimated. The damping and frequency estimates of the signal can be extracted from each zi , as shown in (2). ( ) i i i i i j t z ω σ λ λ ± = ∆ = exp (2) This method also utilizes singular value decomposition to remove noise from the signal before the estimates are made. Because a matrix-pencil function is not available in MATLAB, a user- defined function is created to calculate the parameters using the matrix-pencil method. It is described briefly by the following equations. A matrix [Y] is created from the noisy PMU data ( ) k y , as shown in (3). 0 500 1000 1500 2000 -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 FFT Method Data Window Peak Magnitudes Time, seconds |Y(f)| Peak Magnitudes > 3*std Peak Magnitudes Peak Magnitudes Mean 3*std -3*std This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 4 [ ] ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )             − − − − + = 1 1 1 2 1 1 0 N y L N y L N y L y y y L y y y Y        (3) The parameter L is the pencil parameter, and N is the total number of data points in ( ) k y . Singular value decomposition is applied to [Y], as shown in (4). [ ] H V U Y Σ = (4) The matrices U and V are unitary matrices, and the operator H is the conjugate transpose. The singular values of [Y] are located along the diagonal of matrix Σ in descending order. A threshold value is selected so that the n singular values above the threshold are considered and used to create the following matrices. [ ] [ ] [ ] [ ] [ ] [ ] H H n V U Y V U Y v v v V ' ' ' ' , , , ' 2 2 1 1 2 1 Σ = Σ = =  (5) The vs are the column vectors of V from (4), matrix ' Σ is the first n columns of Σ , [ ] ' 1 V is the matrix V with the final column removed, and [ ] ' 2 V is the matrix V with the final row removed. Using the matrix-pencil method as defined in Liu, Quintero, and Venkatasubramanian (2007), the matrices from (5) can be used to create the following new matrix-pencil definition for noisy data. [ ] [ ] [ ] + + = H H V V Y Y ' ' 1 2 1 1 (6) The next step is to find the eigenvalues for [ ] + H H V V ' ' 1 2 . The eigenvalues are used to solve for R from ( 1 ) and is rewritten in the following format. ( ) ( ) ( )                         =             − − − − n N n N N n R R R z z z z z z N y y y          2 1 1 1 2 1 1 2 1 1 1 1 1 1 0 (7) The estimates for the damping and frequency are extracted from z , as shown in (2), and the amplitude is extracted from R . The matrix-pencil method is applied to a 10-second window of PMU data. Figure 4 shows the reconstructed signal (solid line) using the estimated parameters compared to the original signal (dotted line). Because noise was removed using singular value decomposition, parameters were not estimated to fit the noise but only the “true” modes present in the data. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 5 Figure 4. An example of the matrix-pencil method used to estimate parameters of the de-noised signal Figure 5. In the matrix-pencil method, data windows during which the peak magnitude exceeds three times the standard deviation of peak magnitudes for the entire hour are tagged. The values of the two highest estimated amplitudes are saved for each data window. The method is repeated for the entire hour of data. It is assumed that the max amplitude would be significantly larger if an event is present in the data than the max amplitude from a data window without any events. After the two highest magnitudes are saved for all data windows, the averages and standard deviations for the highest and second highest magnitudes are calculated separately. The method used to detect possible events is similar to the FFT method; any magnitudes above three standard deviations are tagged as possible events, as shown in Figure 5 . Yule-Walker Spectral The Yule-Walker spectral method is similar to the FFT method but utilizes the pyulear function from the MATLAB Signal Processing Toolbox (2011). The pyulear function calculates the power spectral density (PSD) using the autoregressive Yule-Walker method. The PSD for a data window containing an event is shown in Figure 6 . 0 2 4 6 8 10 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 Time, seconds Voltage Phase Angle (detrended), degrees Matrix-Pencil Original Signal Estimated Signal 0 500 1000 1500 2000 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 Matrix-Pencil Method Peak Amplitudes Time, seconds Amplitude, degrees Peak Amplitude > 3*std Peak Amplitude Peak Amplitude Mean 3*std -3*std This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 6 Figure 6. Yule-Walker method power of frequency content in a PMU signal during a large disturbance Figure 7. In the Yule-Walker method, data windows during which the peak magnitude exceeds three times the standard deviation are tagged. The peak value in the PSD in dB is marked. The peak values for each PSD of the data window are saved. This is repeated for the entire hour. The averages and standard deviations of the maximum magnitudes are calculated. Magnitudes above three standard deviations are tagged as possible events, as shown in Figure 7 . Min-Max For this method, the difference between maximum and minimum values within the data window is calculated. The difference for each data window for the entire hour is saved. This method is sensitive because it also detects gradual changes in the signal above the 10-second data window that are not caused by power system events. Possible events are detected by calculating the averages and standard deviations for all of the data window differences for the entire hour of PMU data. The differences that exceed three times the standard deviation are marked as possible events, as shown in Figure 8 . 0 0.5 1 1.5 -90 -80 -70 -60 -50 -40 Yule-Walker Method Frequency, Hz Power, dB 0 500 1000 1500 2000 -90 -80 -70 -60 -50 -40 Yule-Walker Method Peak Magnitudes Time, seconds Peak Magnitudes (dB) Peak Magnitudes > 3*std Peak Magnitudes Peak Magnitudes Mean 3*std -3*std This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 7 Figure 8. In the min-max method, data windows during which the difference between the minimum and maximum values in a 10-second data window exceeds three times the standard deviation are tagged. Event Detection After all methods have been used to screen for possible events in the PMU data, the time stamps of data windows marked as containing possible events are compared. Figure 9 shows an example of the results for possible events in the RPAD between the University of Texas at Austin (UT- Austin) and McDonald PMUs. If two or more methods detect a possible event in the same data window, that data window is marked as containing an event and the PMU data for that window is saved and plotted. In Figure 9 , an event is detected at 1,285 seconds by both the matrix-pencil and min-max methods. Another event is detected at 1,935 seconds by all four methods. Figure 9. Comparison of all of the screening algorithms 0 500 1000 1500 2000 -0.5 0 0.5 1 1.5 2 Min-Max Method Time, seconds Min-Max, degrees Min-Max > 3*std Min-Max Min-Max Mean 3*std -3*std 0 500 1000 1500 2000 2500 3000 3500 4000 Time, seconds Comparison of Event Screening Methods (RPAD between UT-Ausitn - McDonald) Yule-Walker FFT Min-Max Matrix-Pencil (Amp 1) Matrix-Pencil (Amp 2) Event This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 8 Figure 10. An entire day of PMU frequency from UT-Austin PMU Data Formats PMU Data Description This program was created to analyze data from customer-level (120-V) voltages. The data are taken from the Texas Synchrophasor Network. Because they are taken from the 120-V level, only voltage magnitude, angle, and frequency are available. Current and power data are not available and are not considered in the program. Present Data Format A diagram of equipment at each PMU location and central phasor data concentrator within the synchrophasor network is shown in Figure 11 . All equipment was donated by Schweitzer Engineering Laboratories, Inc. As shown in Figure 11 , each PMU station has an antenna to receive a global positioning system signal. This global positioning system signal is the common time reference required to make synchronous phasor measurements. The global positioning system signal is sent to the global positioning system receiver, where it is converted to Coordinated Universal Time and used by each PMU to time stamp the measured phasor quantity. The time-stamped phasor from each PMU station is sent either through public internet or serial cable to the phasor data concentrator. The phasor data concentrator waits for data from all PMU locations and sorts the PMU data after it is received. The collected data are sent to a dedicated desktop PC, where software provided by Schweitzer Engineering Laboratories, Inc., archives the PMU data and displays it in real time. 0 5 10 15 20 59.75 59.8 59.85 59.9 59.95 60 60.05 60.1 60.15 60.2 60.25 Time, hours Frequency, Hz One Day of PMU Frequency from UT-Austin This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 9 Figure 11. (Left) Equipment for each PMU station and (right) the phasor data concentrator and PC within the Independent Texas Synchrophasor Network The synchronized phasor measurements and calculated frequency are stored at a rate of 30 data points per second. Because the data has a 30-Hz sampling rate, only low-frequency oscillations (below 15 Hz) can be analyzed. Events such as capacitor bank switching, which induces oscillations on the order of a few hundreds of hertz to 1.5 kHz, are not visible in the PMU data. Even with a sampling rate of 30 Hz, the archived hourly reports are very large. For one PMU, an hourly report of the voltage phase angle data consists of 108,000 data points. Because there are three signals for each PMU station (voltage phase angle, voltage magnitude, and frequency), each PMU station generates 324,000 data points each hour. These data points are stored hourly in a comma-separated value (CSV) format. A network with five PMU stations generates hourly CSV reports that are 13.8 MB in size. The synchrophasor network is in operation at all times. Protocol to Read in the PMU Data New PMU stations are installed and removed from the network for various reasons, causing the format of the archived hourly PMU files to constantly change. The grouping of PMU quantities for each station changes as well, resulting in two types of data files. The two types of data files are shown in Table 1 and Table 2 . The voltage magnitudes are arbitrarily selected. All measurements are taken at a customer-level voltage of 120 V. Table 1. An example of Type 1 PMU Data in CSV Format for Two PMU Stations With the Magnitude, Angle, and Frequency for Each Station Grouped Together t MAG1 (V) ANG1 (deg) FREQ1 (Hz) MAG2 (V) ANG2 (deg) FREQ2 (Hz) … … … … … … … 24:26.7 79924.90 120.18 59.98 40.80 -10.12 60.02 24:26.8 79930.60 119.89 59.98 40.81 -9.91 60.02 24:26.8 79933.17 119.60 59.98 40.82 -9.70 60.02 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 10 Table 2. An Example of Type 2 PMU Data in CSV Format for Three PMU Stations With the Frequencies for Each Station Listed in the Last Columns t MAG1 (V) ANG1 (deg) MAG2 (V) ANG2 (deg) MAG3 (V) ANG3 (deg) FREQ1 (Hz) FREQ2 (Hz) FREQ3 (Hz) … … … … … … … … … … 08:38.6 41.192 -22.827 79166. 508 -45.7 73482 99.924 59.949 59.947 60.008 08:38.7 41.192 -23.457 79141. 992 -46.3 0 0 59.948 59.947 0 08:38.7 41.187 -24.068 79138. 844 -46.9 73484. 742 100.14 59.948 59.948 60.009 Because of the two file types and the large variety in the number and location of stations in operation, it is cumbersome to read PMU data and begin analysis of the data. Without an automated method, it is first necessary to open each PMU file and count the number of stations and the order of the stations to conduct an accurate analysis. In addition, dropped data points appear in the voltage magnitude and frequency as zeros and need to be removed before analysis can begin. A protocol was created to analyze the PMU data files and indicate the type of file format, the number of PMU stations, and the quality of the PMU data. The file format type indicates if the PMU frequency for each station is grouped with the other PMU station data or if all PMU frequencies are grouped together, as shown in Table 1 and Table 2 , respectively. The number of PMU stations provides the number of PMU stations in operation at the time the data was archived. The quality of the PMU data is “true” if the PMU station contains zero or a limited number of dropped data points and “false” if the data contains many dropped data points or if the PMU station is not sending data to the phasor data concentrator. The protocol function determines the PMU file type and the number of PMU stations. It also indicates the column locations for each PMU station’s magnitude, angle, and frequency data and if the quality of the PMU data is “false” or contains missing data. An example is shown in Table 3 . This information is used to determine how to label the data type and how to label the data by PMU location. Table 3. An Example of the Protocol Function Output Type Number of Measurement Points (NMP) NMP1 MAG ANG FREQ Col Index Fault Col Index Fault Col Index Fault 2 6 1 “true” 2 “true” 13 “true” Preprocessing the Data Unwrapping Voltage Phase Angle After the file format type is known, the columns of the PMU data are identified and the voltage phase angle data are modified before analysis begins. Because the rotation of the voltage phasors is not constant and not at nominal frequency (60 Hz), the voltage phase angle increases or decreases depending on whether the speed of rotation is greater or less than nominal. When the This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 11 voltage phasor is rotating at greater-than-nominal speed and crosses from 180 to -180 degrees, the voltage phase angle jumps, as shown in Figure 12 . The voltage phase angle at each PMU station is unwrapped before the difference in the voltage angle between the two stations is taken. The unwrapped RPAD is shown in Figure 13 and can now be screened for events. All possible RPAD combinations are considered for each PMU file. Figure 12. Raw voltage phase angle for two PMU stations Figure 13. RPAD between the same two PMU stations corrected by individually unwrapping and taking the difference Angle Differences The user needs to determine which angle differences are to be examined. Angle differences between all stations are calculated and available for use. The next part of this manual describes how to enter the user choices. Detection of Dropped Data Points In the Texas Synchrophasor Network, data points were dropped mostly because of loss of internet connection. A drop appears as zero magnitude, zero angle, and zero frequency for that 12 14 16 18 20 22 -200 -150 -100 -50 0 50 100 150 200 Time, minutes Voltage Phase Angle, degrees Raw Voltage Phase Angle UT-Austin PMU McDonald PMU 12 14 16 18 20 22 -0.76 -0.74 -0.72 -0.7 -0.68 -0.66 -0.64 -0.62 -0.6 -0.58 -0.56 Time, minutes Voltage Phase Angle, degrees 'Unwrapped' Relative Phase Angle Difference (RPAD) This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 12 particular time stamp in the data record. It is necessary to record this in order to preprocess the data so that erroneous detection of system events does not occur. Data Window Length To screen for events within each hourly PMU file, an overlapping, moving window is applied to the PMU data. The window size is 10 seconds long and overlaps half of the previously windowed data. The moving window method is illustrated in Figure 14 . tWindow ... Figure 14. (Left) An illustration of a moving window applied to PMU data and (right) a 10-second data window for analysis Each window of data is analyzed for possible events by applying the four different detection methods mentioned. All methods except for the min-max are based on the strength of the estimated frequency content. Only frequencies below 2 Hz are examined (the frequency range of inter- and intra-area oscillations) (Rogers 1999). Each method is described in detail below. Detection of Events Algorithm Description The term event in the report is used to describe large disturbances on the power system and sudden changes in power output at wind farms. Large disturbances include but are not limited to short circuits on transmission lines; protective relay actions required to clear a fault, such as transmission line trip and reclosing actions; or a large and sudden loss of generation or load. A sudden change in renewable energy includes but is not limited to large ramping up or down of wind farm power output. Large disturbances excite a single or multiple modes within a power system; these can be measured as oscillations in PMU voltage data. These oscillations typically consist of one strong mode in which the frequency of the oscillation is between 0.1 Hz and 2.0 Hz, the amplitude of the oscillation is typically higher than oscillations induced by random load switching, and the oscillation is quickly damped out if the system is stable. For example, a low- frequency oscillation of 0.3 Hz with a damping constant of 10% has a time constant of 5.3 seconds (the signal would drop to 37% of its original value after 5.3 seconds). This information is used to detect events in the PMU data. However, it is difficult to verify the cause of many of the detected events. The root cause of a few of the detected events can be identified; however, not all causes are known because information on power system disturbances is not always This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 13 available or shared. Most of the guessed causes of events are based on clues present in the PMU data. Also, wind power–induced events are difficult to describe because detailed measurements of renewable generation power output is necessary for a thorough analysis. The purpose of the PMU network is to conduct analysis on transmission system behavior, but because the measurements are taken at the customer-level voltage of 120 V, events on the distribution system may be measured and recorded as well. However, the voltage magnitude at the customer level is more heavily impacted by events that occur on the distribution system. Also, large events that occur on the transmission system are visible at multiple PMU locations instead of at a single PMU location, as is the case for distribution events. If it is unclear whether a detected event occurred on the distribution system or the transmission system, the data from multiple PMUs can be compared during the time of the event and the voltage magnitude can be checked as well. If the event appears at only one PMU location and the voltage magnitude drops significantly (to 95% or less of rated voltage), then the event should not be considered a transmission system event. This type of event should be attributed to distribution system changes and should not be considered for transmission system analysis. Therefore, by using the guide- lines above, it is important to distinguish between these two types of events so that distribution system events are not incorrectly identified and interpreted as transmission system events. Categorization of Events Events in the PMU RPAD and frequency signals detected by the screening algorithm are analyzed further. To extract common characteristics, the events detected by the screening algorithm are placed into event categories. Event categories are created for events detected in the RPAD signal and for events detected in the frequency signal. These categories were selected after examining many events that occur at different times of the day, seasons, and are measured at different locations within the PMU network. A variety of methods are applied to the data to extract different characteristics from each group. The numerical characteristics from each event category are described in in the following tables. The ranges for each type of event are described for RPAD and frequency events in Table 4 and Table 5 , respectively. Refer to Allen, Santoso, and Muljadi (2013) for more details on these categories. Table 4. Summary of Characteristics From Events Detected in RPAD Data Category Subcategory MAG (deg) Ramp Rate (deg/sec) Duration FREQ Step Change Impulse Single (small) 0.18–0.4 0.5–5.5 0.167– 0.5 - - Single (large) 1.4–1.7 7.0–16.0 0.167– 0.5 - - Transient Damped 0.15–1.5 - <9 0.4– 1.1 - Sustained 0.1–1.6 - >9 0.8– 3.0 - Phase Angle Change Momentary (small) - 0.75–3.25 0.5–1.5 - 0.2–0.6 Momentary (large) - 2–10 0.5–1.5 - 0.95–1.15 Sustained - 1–5 - - 0.15–0.75 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 14 Table 5. Summary of Characteristics from Events Detected in Frequency Data Category Subcategory MAG (mHz) MAG Peak to Peak (mHz) Ramp Rate (mHz/ sec) Ramp Rate Peak to Peak (mHz/ sec) Duration (sec) FREQ (Hz) Delta FREQ (Hz) Impulse Single 3–9 - 10–900 - 0.267– 0.5 - - Multiple 20–25, 3–10 40–50 25–250 250–500 0.3– 0.6333 - - Transient Damped 5–325 - - - 0.4–7 0.3– 5.0 - Sustained 0.1– 1.6 - - - >9 0.8– 3.0 - Rise or Drop in Frequency Rise - - 0.8–3 - - - 0.015– 0.04 Drop - - 2–8 - - - 0.02– 0.1 Unit Trip - - 15–30 - - - 0.1– 0.21 PART 2 The zip file package accompanying this document contains the following files: 1. mainProg.m —This is the main program that calls the data processing functions. 2. Functions (described in the appendix): A. fReadPMUFile.m B. fPMU_Protocol.m C. fAngleProcess_singlephase_unwrap.m D. fMatrixPencil.m E. fpeakdet.m F. fPlotandAnalysis.m G. fEventCount.m 3. A .csv file containing PMU data (SEL PMU format; other formats may not work) Note that the MATLAB Signal Processing Toolbox is required for this analysis, specifically for the Yule-Walker method. Program Inputs This section describes the required inputs to run the software and codes. It is assumed that the data are in one of the two formats described previously. All PMU files should be located in the same directory as the MATLAB code. The first step is to open the code ( mainProg.m ) in MATLAB. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 15 The first input is the file name. It is entered into the fileName variable by the user. Data files are assumed to be in the .csv format. The format of the file name consists of the year, month, and day (yy,mm,dd), followed by the time stamp (hr,min,sec,msec) of the first data-point in the file. It is assumed by the code that the .csv file contains one hour of data. The file name will also contain the name of station at which the data was taken, but the program does not use this information. An example of a file name is given as follows: From this variable input, the date and hour are extracted. The next step in the code is to open the file and extract the data. The protocol sub-function determines the format of the PMU file, the number of stations, and whether the data from each station is of good quality (i.e., data are not dropped). The typical output of the protocol sub- function is described in Table 3 . The function fReadPMUFile will read the file using the protocol and output the time stamp, station titles, station data (which is a three-dimensional array, with rows as data points and columns as voltage magnitude, frequency, and voltage angle for each station), the names of all possible angle differences between stations (e.g., “between Station 1 and 2,” “between Station 1 and 3,” and so on), and finally the actual angle difference. The next step is to decide the window size ( window_rng ). The window size is the product of the samples/second with the number of seconds to be included in the window. An example of a 10- second window for a PMU sample rate of 30 samples/second is shown below. The next step is to select which signal to analyze and the number of stations to analyze. The above code can be changed by the user to determine which signal to use and at which stations to look. Variable if_freq is set to 1 if the frequency signal will be analyzed and the stations listed in variable freq_station will be examined. The frequency station data are located in the stationData variable. Variable if_freq is set to 0 if the phase angle will be analyzed and the angle difference between stations listed in variable angle_station will be examined. The angle difference data are located in the angle_difference variable. After these choices have been made by the user, the program will run to detect events in the chosen signal in the selected hour of data. Program Steps and Outputs A FOR loop is employed to conduct the same analysis for each signal specified. fileName = '120103,010000000,UT,Austin,3378,Phasor.csv'; window_rng = 30*10; % number of data points if_freq = 0; freq_station = [1,2]; angle_station=[1,3,4]; This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 16 The first output is the plot of the first signal in the list. For the above code snippet, the first signal is the first angle difference in the variable angle_difference . An example is shown below in Figure 15 . Figure 15 The output plot of the first signal listed in the variable angle_difference For the analysis, the moving window as previously described is applied to the signal. An example window for the first 10 seconds (300 data points) is given below. This is the first 10 seconds from the signal plotted above. Figure 16. One window of data (10 seconds/300 data points) Each type of signal-processing method (FFT, etc.) is then applied to the window of data, and max values from each analysis are saved. An example of the maximum value from the results of the FFT method applied to the above window of data is shown in the figure below. 0 500 1000 1500 2000 2500 3000 3500 4000 -63 -62.8 -62.6 -62.4 -62.2 -62 -61.8 -61.6 -61.4 -61.2 Austin V1LPM Angle and WACO V1YPM Angleangle pair 1 0 50 100 150 200 250 300 -61.6 -61.55 -61.5 -61.45 This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 17 Figure 17. FFT of window with peak detected The max value indicates the strength of oscillations or changes in the window of data being analyzed. This step is repeated until the entire signal is analyzed. After the entire signal is analyzed, the saved max values from each method are evaluated to detect suspected events. The peak values taken from each window of data are shown in the figure below. Values that lie above 3.5 standard deviations from the mean (indicated by the red line) are considered suspect events. Figure 18. Peak of FFT for every 10-second window in the data, with the (black line) mean and (red line) 3.5 standard deviations from the mean delineated If the suspected event shows up in the results from multiple processing methods, it is classified as an event. An example of all suspected events from each of the methods is given in the figure below. When three or more methods detect a suspect event in the same data window, the 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 3 x 10 -3 Peak Value from the FFT Method after applied to a Window of Data Frequency, Hz Magnitude 0 500 1000 1500 2000 2500 3000 3500 1 2 3 4 5 6 7 8 x 10 -3 Peak Values for Each Data Window Window No. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 18 algorithm defines this window of data as containing an event. These are indicated by an X in the figure below. Figure 19. Screening of suspect events to find events detected by three or more methods. (Such events are marked with an X.) After the events are detected, the events in the original signal can be examined and plotted. Examples of events detected in the phase angle are shown below. Figure 20. An example of an RPAD single (large) impulse event (RPAD between UT-Austin and McDonald) 0 500 1000 1500 2000 2500 3000 3500 4000 Time, seconds Comparison of Event Screening Methods (RPAD between UT-Ausitn - McDonald) Yule-Walker FFT Min-Max Matrix-Pencil (Amp 1) Matrix-Pencil (Amp 2) Event 570 575 580 585 48.5 49 49.5 50 Time, seconds of hour Voltage Phase Anlge, degrees This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 19 Figure 21. An example of an RPAD damped transient Event 1 (RPAD between UT-Austin and Waco) Figure 22. An example of an RPAD sustained transient (RPAD between UT-Austin and McDonald) 645 650 655 660 -65.1 -65.05 -65 -64.95 -64.9 -64.85 -64.8 -64.75 -64.7 Time, sec of hour Voltage Phase Angle, degrees 0 100 200 300 400 500 600 700 -8 -7.5 -7 -6.5 -6 -5.5 -5 Time, sec of hour Angle, degrees This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 20 Figure 23. A contour plot of the PSD (low-frequency content over time) for a voltage RPAD sustained transient event (RPAD between UT-Austin and McDonald) Figure 24. An example of an RPAD sustained step change event (RPAD between McDonald and Waco) Time, sec of hour Frequency, Hz PSD for RPAD between UT-Austin and McDonald 0 100 200 300 400 500 600 700 0 1 2 3 4 5 -150 -100 -50 3575 3580 3585 3590 -57.1 -57 -56.9 -56.8 -56.7 -56.6 -56.5 -56.4 Time, sec of hour Voltage Phase Angle, degrees This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 21 Figure 25. An example of an RPAD momentary step change event (RPAD between UT-Austin and Waco) Frequency events will look different from the results for the phase angle shown above. Refer to Allen, Santoso, and Muljadi (2013) for examples of events in frequency data. Conclusion This user guide, when used in conjunction with the report presented in Allen, Santoso, and Muljadi (2013) and the accompanying software package provides a basic introduction to PMU data analysis for power engineers. Although the results shown here are based on data from the Texas Synchrophasor Network, these same analysis codes can be used with minor modifications on any available PMU data. Future work will address steps to take after the events are detected and classified. References Allen, A.; Santoso, S.; Muljadi, E. (2013). Algorithm for Screening Phasor Measurement Unit Data for Power System Events and Categories and Common Characteristics for Events Seen in Phasor Measurement Unit Relative Phase-Angle Differences and Frequency Signals . NREL/TP- 5500-58611. Golden, CO: National Renewable Energy Laboratory. Accessed March 24, 2014: http://www.nrel.gov/docs/fy13osti/58611.pdf Electric Reliability Council of Texas (ERCOT). (2004–2010). Daily Grid Operations Reports. Austin, TX. Accessed March 24, 2014: www.ercot.com/gridinfo/congestion/operations/ Kundur, P. (1994). Power System Stability and Control. New York: McGraw Hill. Liu, G.; Quintero, J.; Venkatasubramanian, V. (2007). “Oscillation Monitoring System Based on Wide-Area Synchrophasors in Power Systems.” Bulk Power System Dynamics and Control—VII. Revitalizing Operational Reliability, iREP Symposium Proceedings ; Aug. 19–24, Charleston, South Carolina. 3265 3270 3275 3280 -67.5 -67.45 -67.4 -67.35 -67.3 -67.25 -67.2 -67.15 Time, sec of hour Voltage Phase Angle, degrees This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 22 MATLAB Version 7.12.0. (2011). Natick, Massachusetts: The Mathworks, Inc. Rogers, G. (1999). Power System Oscillations . New York: Springer. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 23 Appendix: Help and Documentation for Functions fReadPMUFile fReadPMUFile —Reads data from the PMU data file (.csv). Syntax [ timeStamp , stationTitle , stationData , angle_diff_title , angle_difference ] = fReadPMUFile(fileName , plts ) Description The arguments for fReadPMUFile are: • fileName —The name of the .csv file containing the PMU data • plts —Set to 1 if plots are desired, 0 if not The function fReadPMUFile will read the .csv file and will output: • timeStamp —The PMU time stamp located in the first column of the .csv file • stationTitle —The names of each station located in the first row of the .csv report • stationData —A three-dimensional array containing voltage magnitude, frequency, and voltage angle for all PMU stations in the .csv file • angle_diff_title —The names of all possible angle differences between stations (e.g., “between Station 1 and 2,” “between Station 1 and 3,” and so on) • angle_difference —The calculated angle differences between stations. (All angles are unwrapped using fAngleProcess_singlephase_unwrap .) Related Functions fPMU_protocol , fAngleProcess_singlephase_unwrap Example When the plts input is set to 1, the function plots the signals within the CSV file as well as every calculated angle difference (unwrapped; see fAngleProcess_singlephase_unwrap ). An example plot is shown below. fileName =120103,010000000,UT,Austin,3378,Phasor.csv; plts = 1; [timeStamp, stationTitle, stationData, angle_diff_title, angle_difference] = fReadPMUFile(fileName, plts) This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 24 Figure 26. Example plot from function fReadPMUFile fPMU_Protocol fPMU_protocol —Finds the format of the PMU data file Syntax [ Protocol , Vmeters , title ] = fPMU_protocol(fileName) Description The argument for fPMU_protocol is: • fileName —The name of the .csv file containing the PMU data The function the fPMU_protocol will read the .csv file and will output: • Protocol —The table describing the format of the PMU data file • Vmeters —The PMU time stamp, voltage angle, magnitude, and frequency are stored in this matrix. • title —The names and labels of all of the stations and data The function fPMU_protocol indicates the PMU file format, the number of PMU stations, and the quality of the PMU data. This function is necessary because of the changing format of the PMU data within the Texas Synchrophasor Network. Related Functions fReadPMUFile 0 500 1000 1500 2000 2500 3000 3500 4000 -111.5 -111 -110.5 -110 -109.5 -109 HARRIS V1LPM Angle and UT Pan Am V1LPM Angle This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 25 Example The example shows that for the CSV input provided, the function outputs Protocol , Vmeters , and title . A portion of the output Protocol is provided in the table below. Table 6. Output Protocol CSV File Type Number of Measure- ment Points ‘Austin_V 1LPM_Ma gnitude’ ‘Austin_V 1LPM_Ang le’ ‘Z_UT_337 8_AO[079] _Value’ ‘type2’ 6 ‘1’ ‘true’ ‘2’ ‘true’ ‘13’ ‘true’ The output Vmeters is a matrix of all of the data contained in the original PMU CSV file, and the output title contains all of the headers for the data from the original PMU CSV file. fpeakdet 2 PEAKDET —Detects peaks in a vector [ MAXTAB , MINTAB ] = PEAKDET(V , DELTA) —Finds the local maxima and minima (“peaks”) in the vector V . MAXTAB and MINTAB consist of two columns. Column 1 contains indices in V, and Column 2 contains the found values. With [ MAXTAB , MINTAB ] = PEAKDET(V , DELTA , X) the indices in MAXTAB and MINTAB are replaced with the corresponding X-values. A point is considered a maximum peak if it has the maximal value and was preceded (to the left) by a value lower by DELTA . Example The figure below shows an example signal with the peaks detected based on the DELTA threshold set in the code below. 2 This function and help is written by Eli Billauer, available online fileName =120103,010000000,UT,Austin,3378,Phasor.csv; [Protocol, Vmeters, title] = fPMU_protocol(fileName) delta = 0.025; [maxtab, mintab] = fpeakdet(signal, delta); figure; plot(signal) title('Signal'), xlabel('Data Point No.') hold on plot(maxtab(:,1),maxtab(:,2),'ro') % plot only the max peaks legend('Signal', 'Detected Peaks') This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 26 Figure 27. Signal with peaks detected fMatrixPencil fMatrixPencil —Decomposes the input signal into the sum of complex exponentials 𝑦 � ( 𝑘 ) = � 𝑎𝑚𝑝 ( 𝑛 ) ∗ exp ��𝑎𝑙𝑝ℎ𝑎 ( 𝑛 ) + 𝑗 ∗ 2 𝜋 ∗ 𝑓𝑟𝑒𝑞 ( 𝑛 ) � ∗ 𝑘 + 𝑗 ∗ 𝑡ℎ𝑒𝑡𝑎 ( 𝑛 ) � ; 1 ≤ 𝑘 ≤ 𝑁 𝑁 𝑛=1 where N is the length of signal y_noise . Syntax [ amp , theta , freq , alpha_percent , y_hat ] = fMatrixPencil(y_noise , L , dt , tol) Description The arguments for fMatrixPencil are: • y_noise —Input signal containing noise • L —Matrix-pencil parameter (see reference) • dt —Sample time of y_noise • tol —Tolerance. It determines number of estimated complex exponentials (i.e., a lower value leads to a higher number of estimated complex exponentials). The function fMatrixPencil will output: • amp —The vector containing the amplitude of each complex exponential • theta —The vector containing the angle theta of each complex exponential • freq —The vector containing the frequency of each complex exponential 0 50 100 150 200 250 300 -0.025 -0.02 -0.015 -0.01 -0.005 0 0.005 0.01 0.015 0.02 0.025 Signal Data Point No. Signal Detected Peaks This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 27 • alpha_percent —The vector containing the damping coefficient for each complex exponential • y_hat —The estimated signal from the summation of the complex exponentials. (Noise may be removed depending on the tol setting.) This function represents the input signal as a sum of complex exponentials. The matrix-pencil method is a one-step linear technique to accomplish this decomposition that is robust to noise, unlike polynomial methods (such as Prony analysis). Examples Singular value decomposition (SVD, see MATLAB function svd.m ) is applied, and the resulting sigma values and the selected tol are used to determine the order of the system. An example of the sigmas for a signal y is given below. The lower tol is, the larger the system order will be. Figure 28. Sigma values for the matrix-pencil method An example of the matrix-pencil method as applied to a signal y is given below. The original and matrix-pencil de-noised signals are given. 0 50 100 150 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 σ values n σ L = length(y)/2; % pencil parameter dt = 1/30; % sampling rate is 30 Hz tol = 10e-2; [amp, theta, freq, alpha_percent, y_hat] = fMatrixPencil(y, L, dt, tol); figure; plot(t,y,t,y_hat) legend('Original Signal','De-noised Matrix Pencil Reconstruction') xlabel('Time, seconds') This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 28 Figure 29. Comparison of original and reconstituted signals fEventCount fEventCount —Counts the number of events that were detected and prevents an event from being counted multiple times as a result of overlapping windows Syntax [ t_events_new ] = fEventCount(t_events_old , t_new) Description If any new events are in array t_new that do not belong to t_events_old , these new events will be included; if the event is already in t_events_old , it will be ignored and the event count (second column of t_events_old ) will be updated. Input • t_events_old —An n by 2 matrix in which the first column of the matrix includes all of the time stamps for the detected events and the second column contains the number of times the event has been detected • t_new —Array of new time stamps from a different detection method Output • t_events_new —New array of time stamps for the detected events from the old array t_events_old and from the events detected by the new method in t_new fPlotandAnalysis fPlotandAnalysis —Flags outliers and plots events detected by the method given Syntax [ pt_event ] = fPlotandAnalysis(timeStamp , signal , t , xAnalysis , std_cuttoff , if_diff , if_plot , sig) 0 1 2 3 4 5 6 7 8 9 10 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 Time, seconds Original Signal De-noised Matrix Pencil Reconstruction This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 29 Description Using the data from the analysis, this function finds events that are three (or other) times the standard deviations above the mean. The function return pt_event is a placeholder of where the events are in xAnalysis data. Input • timeStamp —Time stamp of the original signal that is being analyzed. This time is used to plot the detected event in the signal. • signal —Original signal that is being analyzed. Portions of this signal where events occur will be plotted. • t —The time at which the analysis window begins • xAnalysis —The peak values of the method results for the window of data analyzed • std_cuttoff —User set value. If a window contains a peak that is this number of standard deviations above the mean, it will be classified as an event-containing window. • if_diff —1 if the function is being applied to method diff, 0 otherwise • if_plot —1 if the plots of the detected events are desired, 0 otherwise • sig —If if_plot is set to 1, sig is the title of the plot. Output • pt_event —The index of the found events in the xAnalysis array Example An example signal is plotted below. The input timeStamp and signal are used in this plot. Figure 30. Example signal The peak values for the method applied to a moving window across the data are saved and are plotted as shown in the figure below. The inputs t and xAnalysis are used for this plot. The figure shows the (black line) mean for these window max values and (red line) 3.5 standard 0 500 1000 1500 2000 2500 3000 3500 -63 -62.8 -62.6 -62.4 -62.2 -62 -61.8 -61.6 -61.4 -61.2 Signal Time, seconds This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 30 deviations above the mean. All max values that are above the red line are considered possible events. In this case, five windows contain possible event data visible in the original signal. Figure 31. Peak of FFT for every 10-second window, with the (black line) mean and (red line) 3.5 standard deviation delineated If if_plot is set to 1, the function will plot a figure similar to the one above. fAngleProcess_singlephase_unwrap fAngleProcess_singlephase_unwrap —Unwraps the single-phase voltage angle data Syntax [ ang ] = fAngleProcess_singlephase_unwrap(V1ang) Description: The argument for fAngleProcess_singlephase_unwrap is: • V1ang —Raw PMU phase angle data The function fAngleProcess_singlephase_unwrap will output: • ang —Unwrapped PMU phase angle data Unwrapping angle data refers to the process of removing jumps that occur when the angle recorded by the PMU crosses 0/360 degrees or 180 to -180 degrees. The unwrapped angle is monotonically increasing with no jumps or discontinuities. Example A plot of the raw angle data before the function is applied is shown below. 0 500 1000 1500 2000 2500 3000 3500 1 2 3 4 5 6 7 8 x 10 -3 Peak Values for Each Data Window Window No. This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications. 31 Figure 32. Raw angle data with jumps and discontinuities present The processed angle data are given below. Figure 33. Processed angle data with discontinuities removed Reference Liu, G.; Quintero, J.; Venkatasubramanian, V. (2007). “Oscillation Monitoring System Based on Wide Area Synchrophasors in Power Systems.” In Bulk Power System Dynamics and Control- VII. Revitalizing Operational Reliability, 2007 iREP Symposium Proceedings ; August 19–24, 2007, Charleston, South Carolina. IEEE; pp. 1–13. 0 2 4 6 8 10 12 x 10 4 -200 -150 -100 -50 0 50 100 150 200 Raw Angle Data Data Point No. Degrees 0 2 4 6 8 10 12 x 10 4 -1000 -500 0 500 1000 1500 2000 2500 3000 Processed Angle Data Data Point No. Degrees This report is available at no cost from the National Renewable Energy Laboratory (NREL) at www.nrel.gov/publications.