Application of Wavelet Transform and Artificial Neural Network to Extract Power Quality Information from Voltage Oscillographic Signals in Electric Power Systems

Post-operation contingencies analysis in electrical power systems is of fundamental importance for the system secure operation, and also to maintain the quality of the electrical energy supplied to consumers. The electrical utilities use equipments as Digital Disturbance Registers (DDR), and Intelligent Electronics Devices (IED) for faults monitoring, and diagnosis about the electrical power systems operation and protection. In general, the DDR and IED are intended to monitor the protection system performance and detect failures in equipments and transmission lines, and also generate analog and digital oscillographic registers that better characterize the disturbing events.

One of the main difficulties in using measurements, obtained by DDR, in the evaluation of power quality as compared with those obtained by power quality monitors, is that many of the signal processing stages are not performed automatically by the first.For the oscillographic records to be useful as power quality indicators, it is first necessary to obtain certain parameters to classify the recorded signals according to the event type that has occurred.Considering the case of short duration voltage variations (SDVV), the parameters of interest are the event amplitude and time duration.Obtaining these parameters enables the application of statistical tools as presented in (Bollen, 2000), for results analysis and visualization, which allows having information about the electrical system behavior at certain time intervals, for example, months or years.
Another difficulty, perhaps the most critical, is the large volume of data obtained from oscillographic monitoring.Many of these recorded signals are due to switching maneuvers, or due to spurious signals or noise, without characterizing voltage changes in the electrical system.For this large amount of data to be evaluated, it is necessary that an automatic classification method be used so that only signals with the desired characteristics are used to determine the parameters of interest.This aspect is highlighted in several publications which present new methods for classification and characterization using digital signal processing and computational intelligence tools (Angrisani et al, 1998;Santoso et al, 2000a;Santoso et al, 2000b andHuang et al, 2002;Machado et al, 2009;Rodriguez et al, 2010 ).
The first use of wavelet transform in power systems is credited to (Ribeiro, 1994).In recent years, wavelet transform -WT, a powerful tool for digital signal processing, has been proposed as a new technique for monitoring and analysis of different disturbances types in power systems (Machado et al, 2009;Mokryani, 2010;A. Rodriguez et al, 2010;Gong Jing, 2010, 2011).Wavelets, along with computational intelligence techniques like artificial neural networks and fuzzy logic, have been used successfully in automatic classification of power quality problems.(Machado et al, 2009;Mokryani, 2010;Rodriguez et al, 2010) The present work aims to develop an automated system for classifying power quality problems with respect to the fault type that has occurred and the electric phase involved, and quantify SDVV in electrical power systems from the available oscillography in the electrical utilities post-operation centers, to form a parameter database characterizing power quality problems.The proposed methodology uses the wavelet transform to obtain a characteristic pattern to represent the phenomenon and a probabilistic neural network for classification.

Wavelet transform
Wavelets are functions that satisfy certain mathematical requirements.The wavelet name comes from the fact that they must be oscillatory (a wave), and be well placed, therefore exhibiting short time duration.There are several wavelet types, usually grouped into families, from which the Daubechies is one of the best known.
Wavelets are used to represent data or other functions in a similar way as the Fourier analysis uses sines and cosines.The signal analysis by wavelet transform has advantages over traditional methods using Fourier analysis when the signals have time discontinuities or present a non-stationary oscillatory behavior.
The mathematics main branch leading to wavelet analysis began with Joseph Fourier (1807) with his frequency analysis theory, known as Fourier analysis.The first wavelet mention appears in the appendix of A. Haar's thesis (1909).Paul Levy a 1930's physicist, investigating the Brownian motion, found that the Haar basis functions are superior to the Fourier basis functions for studying small and complicated details in the Brownian motion.In 1980, Grossman and Morlet, broadly defined wavelets in the context of quantum physics, providing a way of thinking about wavelets based on physical intuition.In 1985, Stephane Mallat gave wavelets an additional advance.Through his work in digital signal processing, he discovered some relationships among quadrature mirror filters -QMF, pyramidal algorithm, and orthogonal wavelet basis.Based partially on these results, Y. Meyer built the first non-trivial wavelets, which unlike the Haar wavelet, the Meyer wavelets are continuously differentiable, but do not have compact support.Years later, Ingrid Daubechies used Mallat's work to build a set of wavelets with orthogonal basis functions that have become the cornerstone of wavelet applications today.

Wavelet analysis
The wavelet transform is a technique similar to the windowed Fourier transform with the difference that the window width is variable.The wavelet analysis allows the use of large time intervals when it is desired to get low frequency information and shorter time intervals when the interest is to obtain high frequency information.Unlike Fourier analysis that uses sines and cosines, wavelet analysis uses wavelets.Figure 1 shows as an example, the Daubechies wavelet, db8.Wavelets sets are employed to approximate signals, and each set consists of scaled versions (compressed or expanded) and translated (time shifted) from a single wavelet, called mother wavelet.

Discrete wavelet transform
In the discrete wavelet transform the term "discrete" applies only to the parameters in the transformed domain, that is, scales and translations, and not to the independent variable time, of the function being transformed.The discrete wavelet transform provides a set of coefficients corresponding to points on a grid or two-dimensional lattice of discrete points in the time-scale domain.This grid is indexed by two integers, the first, denoted by m , corresponds to the discrete steps of the scale, while the second, denoted by n , corresponds to the discrete steps of translation (time displacement).The scale a becomes , where 0 a and 0 b are the discrete steps of the scale and translation, respectively (Young, 1995).Then the wavelet can be represented by: The discrete wavelet transform is given by: where, , mn Z  , and Z is the set of integer numbers.
The parameter m which is called level, determines the wavelet frequency, while the parameter n indicates its position.
The inverse discrete wavelet transform is given by: where k is a constant that depends on the redundancy of the combination of the lattice with the used mother wavelet (Young, 1995).
Along with the time-scale plane discretization, the independent variable (time) can also be discretized.The sequence of discrete points of the discretized signal can be represented by a discrete time wavelet series DTWS.The discrete time wavelet series is defined in relation to a discrete mother wavelet, () hk .The discrete wavelet time series maps a discrete finite energy sequence to a discrete grid of coefficients.The discrete time wavelet series is given by (Young, 1995). (4)

Multiresolution analysis
Multiresolution Analysis -MRA, aims to develop a signal () ft representation in terms of an orthogonal basis which is composed by the scale and wavelets functions.An efficient algorithm for this representation was developed in 1988 by Mallat (Mallat, 1989) considering a scale factor 0 2 a  and a translation factor 0 1 b  .This means that at each decomposition level m , scales are a power of 2 and translations are proportional to powers of 2. Scaling by powers of 2 can be easily implemented by decimation (sub-sampling) and over-sampling of a discrete signal by a factor of 2. Sub-sampling by a factor of 2, involves taking a signal sample from every two available ones, resulting in a signal with half the number of samples than the original one.Over-sampling by a factor of 2, consists of inserting zeros between each two samples resulting in a signal with twice the elements of the original one.

Analysis or decomposition
The structure of the multiresolution analysis is shown in Figure 2. The original signal passes through two filters, a low pass filter () gk , the function scale, and a high pass filter () hk , the mother wavelet.The impulse response of () hk is related to the impulse response of () gk by (Mallat, 1989): Filter () hk is the mirror of filter () gk and they are called quadrature mirror filters.
In the structure presented in Figure 2, the input signal is convolved with the impulse response of () hk , and () gk , obtaining two output signals.The low pass filter output represents the low frequency content of the input signal or an approximation of it.The high pass filter output represents the high frequency content of the input signal or a detail of it.It should be noted in Figure 2 that the output provided by the filters has together twice the number of samples of the original signal.This drawback is overcome by the process of decimation performed on each signal, thereby obtaining the signal cD , the wavelet coefficients that are the new signal representation in the wavelet domain, and the signal cA , the approximation coefficients which are used to feed the next stage of the decomposition process in an iterative manner resulting in a multi-level decomposition.

Fig. 2. Structure of the multiresolution analysis
The decomposition process in Figure 2 can be iterated with successive approximations being decomposed, then the signal being divided into several resolution levels.This scheme is called "wavelet decomposition tree" or "pyramidal structure" (Young, 1995 andMisit et al, 2000).Figure 3 shows the schematic representation of a signal being decomposed at multiple levels.Since the multiresolution analysis process is iterative, it can theoretically be continued indefinitely.In fact, the decomposition can proceed only up to 1 (one) detail, consisting of a single sample.The maximum number of decomposition levels for a signal having N samples is given by 2 log N .

Synthesis or reconstruction
The synthesis process or reconstruction is to obtain the original signal from the wavelet coefficients generated by the analysis or decomposition process.While the analysis process involves filtering and sub-sampling, the synthesis process performs a reverse sequence, over-sampling and filtering.The filters used in the synthesis process are called reconstruction filters, being () gk  the low pass filter, and () hk  the high pass filter.Figure 4 shows the reconstruction scheme from a single decomposition stage.It is observed from Figure 4 that to retrieve the original signal, it is necessary to reconstruct details and approximations.Details could be obtained with over-sampling of the cD coefficients, and a subsequent filtering with () hk  .Approximations are obtained with over- sampling of the coefficients cA , and a subsequent filtering with () gk  .The original signal is then obtained by: The scheme presented in Figure 4 can be extended to a multi-level decomposition.

Probabilistic neural network
The structure of a Probabilistic Neural Network (PNN) is similar to a feed forward network.
The main difference is that the activation function is no longer the sigmoid; it is replaced by a class of functions which includes, in particular, the exponential function.The main advantage of PNN is that it requires only one step for training and that the decision surfaces are close to the contours of the Bayes optimal decision when the number of training samples increases.Furthermore, the shape of the decision surface can be as complex as necessary, or as simple as desired (Specht, 1990).
The main drawback of PNN is that all samples used for the training process must be stored and used in the classification of new patterns.However, considering the use of high-density memories, problems with storage of training samples should not occur.In addition, the PNN processing speed in the classification of new patterns is quite satisfactory, and even several times faster than using back propagation algorithms as reported by (Maloney et al, 1989).

The Bayes strategy for pattern classification
One of the traditionally accepted strategies or decision rules used to patterns classification is that they minimize the "expected risk."Such strategies are called Bayes strategies, and can be applied to problems containing any number of categories (Specht, 1988).
To illustrate the Bayes decision rule formalism, it is considered the situation of two categories in which the state of known nature  , can be A  based on a measurements set represented by a n dimension vector x .Then the Bayes decision rule is given by: where ()  is given by: () () where: It should be noted that, in general, the decision surfaces of two categories defined by Eq. ( 8) can be arbitrarily complex, since there are no restrictions on the densities except for those conditions to which all probability density functions must satisfy, namely that they must be always non-negative, and integrable and their integrals over all space be equal to unity.
The ability to estimate the probability density functions, based on training patterns, is fundamental to the use of Eq. ( 8).Frequently, a priori probabilities can be known or estimated, and the loss functions require subjective evaluation.However, if the probability densities of the categories patterns to be separate are unknown, and all that is known is a set of training patterns, then, these patterns provide the only clue to the estimation of that unknown probability density.A particular estimator that can be used is (Specht, 1990): Where i is the pattern number, m is the total number of training patterns, ai x is the i-th training pattern of category A  , and  is the smoothing factor.It should be noted that () A fx is simply the sum of small Gaussian distributions centered at each training sample.

Structure of the Probabilistic Neural Network
The probabilistic neural network is basically a Bayesian classifier implemented in parallel.The PNN, as described by Specht (Specht, 1988)

Proposed procedure
The proposed procedure is shown schematically in Figure 6.The real data file contains phases A-B-C voltages and currents waveforms, as well as digital signals that indicate the statuses of protective devices, as relays and circuit breakers, acquired by DDR and IED installed in the electrical system substations.These raw data are coded in the COMTRADE format for power systems (IEEE Standard Common Format for Transient Data Exchange), (IEEE Std C37.111, 1999).So, to obtain the voltages and currents signals it is firstly necessary to decode the COMTRADE data, and select the desired waveforms to be analyzed.Before inputting the voltages waveforms to the processing stage, a pre-processing routine is accomplished to standardize the raw data due to the different voltage levels that are encountered in the power system topology.In the case study presented here, the power transmission system presents 230 kV and 500 kV voltage levels.The standardization is performed by converting the phase voltages to per unit (p.u.) values considering the voltage peak value as base voltage.

Processing stage
In the processing stage, the wavelet transform is applied to the voltage waveforms to obtaining signals patterns that characterize short duration voltage variations (SDVV) and transient variations (TV) due to system faults.These obtained patterns are used as inputs to two Probabilistic Neural Networks for SDVV classification (PNN1), as well as to classify the fault type that has occurred (PNN2).The classification results will form a database which can be used to evaluate power quality indices for the electrical system.

Input patterns
Power systems electromagnetic phenomena are characterized by categories according to their spectral content, magnitude and duration (IEEE Std 1250, 1995).These phenomena classification into categories requires an analysis methodology that very frequently must be individualized, which prevents this procedure applicability when the number of signals to be evaluated is very large.Then, procedures to extract signals relevant characteristics have been proposed, so that they can be automatically classified into a specific category.Obtaining parameters for characterizing a given signal usually requires a transformation from the time domain to another domain where the specific characteristics are highlighted.
The use of wavelet transform has proved adequate for obtaining electrical signals characteristics which can be used in classification processes.Studies such those presented in (Lee et al, 1997;Chan et al, 2000;Santoso et al, 2000c;Ramaswamy et al, 2003;Zwe-Lee et al, 2003;Zwe-Lee, 2004& Machado et al, 2009), use characteristic vectors based on the multiresolution analysis decomposition levels coefficients as input to computational intelligence-based systems to classify different power quality events.The characteristic vectors magnitudes depend on the number of decomposition levels used for the analysis, or the number of coefficients of a given decomposition level.The method proposed here uses the Daubechies wavelet, db4, and the voltage signals are decomposed into three levels.The first signal detail level is used to determine the time instant the disturbance has started and also to characterize the transients in the fault type identification, while the third signal approximation is used to characterize SDVV.The computational algorithms were implemented on MATLAB, and also coded in Java.
Figure 7(a) shows an original voltage waveform in p.u. obtained from a digital disturbance register (DDR) presenting a voltage sag.The original waveform is decomposed into three resolution levels.In Figures 7(b-d) the signal details from level 1 to level 3 are presented and in Figure 7(e) the signal approximation at level 3. Details retain the high-frequency information contained in the signal, divided into frequency bands which are function of the sampling rate used in the acquisition process.In case of Figure 7, the sampling rate is 96 samples per cycle of 60 Hz, or 5,760 samples per second.The wavelet transform performance to detect disturbances in electrical signals is substantially improved if a procedure for reducing noise level is applied to the decomposition level coefficients to be used in the detection process.This feature is highlighted in (Yang et al, 1999;2000& 2001).So, to better characterize the disturbance location in the signal, it is applied the following algorithm presented in (Misiti et al, 2000), to the previously selected decomposition level: The s  value used was 10% of the maximum absolute value of the decomposition level coefficients considered, as proposed in (Santoso et al, 1997).
A voltage waveform containing voltage sag is shown in Figure 8(a).In (b) it is presented the details level used to detect the disturbance beginning and (c) presents new details values after the noise reduction algorithm is applied.In (c) it can be observed smaller coefficients magnitudes over the entire signal which improves the algorithm performance used to detect the disturbance.The disturbance beginning point is found based on the following algorithm presented in (Gaouda et al, 2002) where: The algorithm (13) was originally proposed to find the disturbance start and end points.In this particular case, the interest is just the starting point, i p , which shall be considered as a reference for obtaining the phenomena pattern characterization in the classification stage.For this purpose the following algorithm is proposed:

End
Once the disturbance starting point is obtained, the next step is to determine the signal parameters to input the PNN in order to characterize SDVV and transients.

SDVV characterization parameters
As the signal magnitude and duration change during the SDVV occurrence, the norm value (Euclidian distance) will also change if the disturbed signal is considered.So, by monitoring changes in the norm of the third-level signal approximation (the level containing the fundamental frequency) and considering the signal disturbed portion, it can be obtained a standard value characterizing these signal changes.Figure 9 shows the signal norm variation as function of the signal magnitude for the third-level signal approximation of the multiresolution analysis.In this analysis, a 10 cycles signal window was considered and the disturbance magnitude ranging from zero to 1.8 p.u. So, the SDVV classification pattern is obtained by calculating the signal norm for 10 cycles counting from point i p , which represents the disturbance starting point.This procedure is applied to the voltage waveforms in phases A-B-C resulting a vector with three elements which is used as input to the PNN for classification purpose.

Transients characterization parameters
In the transient analysis case, a two cycles long window is selected from the disturbance starting point which, for real electrical systems, is a time interval within which most of the protective devices operate.This considered signal is then normalized based on the biggest magnitude coefficient, for creating a vector related to each fault type to be analyzed in the classification task.
In three-phase transmission lines, phases are mutually coupled and therefore the high frequency variations generated during a disturbance may also appear in non-faulted phases.Using a modal transformation allows the coupled three-phase system to be treated as a system with three independent single-phase circuits.Each phase values are transformed into three decoupled modes: mode 0 (zero), mode and mode , so the three phases are decomposed into nine modes, three for each phase.As mode 0 is the same for all phases, this mode can be calculated only once, reducing to seven the number of signals.Therefore, the three phase voltage signals are decomposed by the multiresolution analysis and the firstlevel detail 3-dimensioned array is used with the modal matrix to decoupling the original signals.
Mathematically the modal transformation consists of a matrix operation as follows: Where are the voltage wavelet coefficients corresponding to the coupled and decoupled phases respectively and W is the decoupling matrix.It is noteworthy that only the voltage signals can be decoupled by the method presented here and the operation described in Eq. ( 14) should be performed on each signal sample.The matrix W is described by (Silveira;et al, 1999): This way it is obtained a system that provides seven outputs, being mode and mode for each phase and a mode 0 which is common to the three phases.These modes contain the wavelet transform coefficients of the three-phase decoupled input signals.The linearity properties of the wavelet and modal transformations ensure that they can be carried out in a cascading way without causing problems to the classifier algorithm results.So, it is obtained a classification pattern that is represented by a matrix with seven columns and 192 rows.

Artificial neural networks structures
The ANN used for the SDVV classification, named PNN1, is composed of three classes, namely:  Class 1 -Voltage sags and interruptions, which are characterized by voltage magnitudes smaller than 0.9 p.u.


Class 2 -Adequate voltage, which is characterized by magnitudes between 0.9 p.u. and 1.1 p.u.;  Class 3 -Voltage swell, which is characterized by magnitudes between 1.1 p.u. and 1.8 p.u.
The training values of each class were obtained from points on the curve given in Figure 9

Results
In order to evaluate the performance of the proposed method in classifying SDVV, 311 voltage oscillographic signals obtained from a real power system were used.The oscillographic signals were numbered from 1 to 311 for the purpose of identification.The electrical power system is a 500 kV/230 kV transmission system connecting Tucuruí Hydroelectric Power Plant located in the south of the State of Pará-Brazil, to load centers in the northern region, which is operated by Eletronorte, a generation and transmission utility in the north of Brazil.The oscillography files used are from the 230 kV substation Guamá, located in Belém city, the capital of the state of Pará, and corresponds to a time period within 2004/2005.For the fault type classification and the faulted phase identification the same 230 kV/500 kV electrical power system was used in which short-circuits were simulated along the transmission lines by varying the incidence angle, and the short-circuit resistance to obtaining a set of voltage waveforms corresponding to the different simulated fault types, using the simulation software ATP.
The simulation studies included 1,029 single-phase to ground short-circuit; 2,058 two-phase and two-phase to ground short-circuits; and 686 three-phase and three-phase to ground short circuits.For the PNN2 training, seven case studies for each fault type as listed in Table 1 were used as input patterns, and the remaining cases were used for testing.Table 3 shows the classification results, noting that misclassification occurred for single-phase and twophase to ground short circuits, with 6% and 5.4% respectively.Also 58% of the three-phase short circuit were classified as three-phase to ground short circuits, but considering that these two fault types can be considered as a single class there would be no classification error in this case, as presented by the 100% result in Table 3

Conclusion
This work presented a methodology for automatic SDVV classification as well as fault type identification using digital signal processing and computational intelligence techniques.
Real power system data were used and satisfactory results for both SDVV and fault type classification were obtained.The implementation of the proposed methodology as part of a computational tool and its integration with the post-operational utility analysis routines will enable the automatic analysis of a larger number of signals waveforms, allowing the methodology proposed here to serve as a basis for future applications where automatic analysis procedures are needed.
One should also note that the wavelet used in this work was chosen due to its good performance in determining the disturbance location in the signal waveform.Various wavelets orders from db2 to db16 were tested and the db4 wavelet presented the best performance, and considering also the fact that it has filters with few coefficients, the processing time for the signals decomposition is greatly reduced, which is an important characteristic when a large number of signals are to be analyzed.

Fig. 3 .
Fig. 3. Schematic representation of a signal being decomposed at multiple levels.
probability density functions for categories A  and B  respectively, A l is the uncertainty function associated with the decision () l is the uncertainty function associated with the decision () , the boundary between the regions in which the Bayes decision () A dx   and () B dx  , is based on estimation of probability density functions for the various classes established by the training patterns.A schematic diagram for a PNN is shown in Figure 5.The input layer X is responsible for connecting the input pattern to the radial basis layer.

Fig. 8 .
Fig. 8. (a) Original voltage waveform with voltage sag, (b) second details level, and (c) second details level after noise reduction.

Fig. 9 .
Fig. 9. Signal norm variation as a function of the SDVV magnitude.
, resulting in 19 values stored in the PNN1.As each class covers a different magnitude range, it was established 9 values for class 1, 3 values for Class 2 and 7 values for Class 3. The weight matrix of the competitive layer is a 3x19 matrix, which corresponds to the 19 training values and the three classes considered.The input pattern to be classified consists of a three elements vector, each representing the characteristic of each phase voltage; and the PNN1 output consists of a three elements vector, each one indicating the classification corresponding to each phase.
( ) For transient analysis 11 classes were considered, which correspond to the short circuit types as listed in Table1.

Table 2
shows the results corresponding to the PNN1 output.The SDVV parameters represented in Table2are the time duration in cycles, and magnitude in p.u.As can be seen, 24 voltage signals were classified as having SDVV.According to data in Table2it may be seen that the PNN1 classification is consistent with the magnitude values calculated for the SDVV.It is observed that in most cases voltage sags were detected in all three phases (classification 1.1.1),and for signals 267 and 268 voltage sags were detected only in phase C, while phase A, and B exhibited adequate voltage magnitudes (classification 2.2.1).It is also worth noting that all these results were compared with the real original voltage waveforms, which proved the results correctness as obtained by the wavelet multiresolution analysis and by the PNN1 classification mechanism.

Table 2 .
SDVV classification and quantification results for three-phase voltage signals obtained from oscillographic records in a real electrical power system.

Table 3 .
. Results for fault type classification With the purpose of testing the performance of the proposed method in classifying real oscillographic signals, some Eletronorte operational reports in the period 2007/2008 were analyzed which contained 31 labeled transient occurrences, being 17 due to short circuits, and 14 due to lightning discharge.For considering lightning discharges (LD) a new class was added to PNN2, and 7 of the 14 signals were selected for training the PNN2 and the remaining signals were used for testing.The testing signals were applied to the trained PNN2 achieving 100% accuracy for short circuits and 85,7% for lightning discharges.