Magnetizing inrush current identification using wavelet based gaussian mixture models

Simulation Modelling Practice and Theory 17 (2009) 991–1010

Contents lists available at ScienceDirect

Simulation Modelling Practice and Theory journal homepage: www.elsevier.com/locate/simpat

Magnetizing inrush current identification using wavelet based gaussian mixture models S. Jazebi a, B. Vahidi a,*, S.H. Hosseinian a, J. Faiz b a b

Department of Electrical Engineering, Amirkabir University of Technology, 424 Hafez Ave., Tehran 1591634311, Iran School of Electrical and Computer Engineering, University of Tehran, Tehran, Iran

a r t i c l e

i n f o

Article history: Received 26 August 2008 Received in revised form 11 February 2009 Accepted 13 February 2009 Available online 27 February 2009

Keywords: Gaussian mixture models Differential protection Inrush current Transformer internal fault Pattern recognition

a b s t r a c t The paper presents a novel differential protection approach based on Gaussian Mixture Models (GMM). This method has the advantages of high accuracy and low computational burden. Two common transients such as magnetizing inrush current and internal fault which their mis-identification may lead to mal-operation of differential relays are considered. GMM, the powerful probabilistic pattern classifier is trained with the features extracted by discrete wavelet transform to reduce the computational time of training procedure and enhance the discrimination accuracy. Training data sets are achieved using kmeans clustering algorithm. Based on the proposed algorithm, a high speed differential relaying (a quarter of a cycle) without dependency on a specific threshold is performed. The suitable performance of this method is demonstrated by simulation of different faults and switching conditions on a power transformer in PSCAD/EMTDC software. Sympathetic and recovery inrush currents were also simulated and investigated. The proposed algorithm is also evaluated by the data collected from a prototype laboratory power transformer. It provides a high operating sensitivity for internal faults and remains stable for inrush currents even in noisy conditions. Since the discrimination method is done with stochastic characteristics of signals without application of any deterministic index, more reliable and accurate classification is achieved. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction Differential protection has been employed as the primary protection of power transformers for many years. The major drawback of the differential relays stem from its potential for false tripping caused by the transient magnetizing inrush current, which flows when the transformer is energized. Many researches have focused on implementing methods to restrain tripping command when an inrush current flows in the transformer windings. An appropriate scheme should be reliable in operating conditions, secure while sending tripping commands, and have fast fault clearing performance. The most common technique used to prevent false trips during the initial energization is harmonic restraint relay [1]. Harmonic restrain principle based techniques often encounter problems when second harmonic generates in other operating conditions and the magnitude exceeds the predefined threshold. This may be due to current transformer saturation, system resonant conditions, presence of a shunt capacitor, or presence of nonlinear loads or distribution capacitance in a long extra high voltage (EHV) transmission line connected to the protected transformer [2,3]. However, the 2nd order harmonic component may also be generated during internal faults in the power transformer. To avoid the tripping operation dependency on such a threshold many novel methods have been proposed. These approaches include: autoregressive process based on power spectrum of

* Corresponding author. Tel.: +98 21 64543330; fax: +98 21 66406469. E-mail addresses: [email protected] (S. Jazebi), [email protected] (B. Vahidi), [email protected] (S.H. Hosseinian), [email protected] (J. Faiz). 1569-190X/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.simpat.2009.02.004

992

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

harmonics [2], wavelet transform based schemes [3–9], artificial neural networks [8–10], leakage inductance based techniques [11,12], harmonic characteristics based schemes [13], induced voltage based method [14], correlation transform [15], similarity degree between voltage and current [16], combination of wavelet transform and fuzzy logic [17], Bayes probability theory [18], non-stationary signal analysis using S-transform [19], space vector analysis of differential signals [20], statistical signal processing using Support Vector Machine (SVM) [21], and Hidden Markov Models (HMM) [22,23]. These schemes have their advantages and disadvantages. For instance, the method presented in [2] needs a predefined exact threshold which may be changed for different conditions. Wavelet based methods have better ability of time-frequency analysis but they usually require long data window and are also sensitive to noise and unpredicted disturbances [15]. Artificial Neural Networks can classify patterns perfectly but the results do not determine the classification accuracy. Some of priori methods suffer from being dependant on parameters of protected transformer, protective algorithm calculative cost, and some from voltage transformers requirements. Some are affected in a noisy environment or current distortions. The GMM has been applied specially to speaker identification [24]. It is also used for multiple limb motion classification [25], image processing and classification [26,27], control engineering [28], simplification of controller design [29], tracking elliptical living objects in video frames [30], and rotating machinery fault diagnosis [31], but not yet used in power engineering fault diagnosis. The GMM ability to solve nonlinear and high dimensional pattern problems makes this algorithm a powerful choice in power system transient classifications. A suitable classifier must be as accurate as possible, in the sense that it must be able to generalize well to novel data. It must be capable of being optimized to suite the unique patterns generated by individual users. It is essential for that to be computationally efficient in the act of classifying novel patterns, as it must satisfy the real-time constraints of fault identification procedure. It must also be easy to be implemented. Although not essential, it should be capable of being trained in a reasonable amount of time. The GMM satisfies each of these criterias [25]. This paper has proposed a new combinatorial classification method using GMM and wavelet transform to discriminate between magnetizing inrush current and internal faults in power transformers. First, some basic characteristics of the signals are extracted by using discrete wavelet transform. As shown later, the number of input samples in the system decreases using DWT and leads to less training and testing computation burden of GMM. Furthermore, it increases the discrimination accuracy of GMM. Finally, GMMs are adopted to correctly determine the fault type by a probabilistic manner from extracted features. The suitable performance of this method is demonstrated by simulation of different faults and switching conditions on power transformers. A typical 500 MVA, 400/230 kV, transformer is simulated using PSCAD/EMTDC software. The proposed scheme also has been tested using data collected from experimental tests on a prototype laboratory three phase power transformer. Following section will describe the power transformer modeling. In Section 3, the simulated test system and data-set collection are briefly reviewed. Gaussian Mixture Model is well presented in Section 4. Section 5 includes the proposed method definitions. Section 6 contains test results and discussion. Finally, Section 7 presents conclusion remarks. 2. Power transformer model The classical model of transformer presented in [32] is used for simulations. In this model each phase of the transformer is represented by a single phase transformer without winding to winding coupling. 2.1. Magnetic core saturation modeling The magnetic core saturation is modeled by compensating current source across the winding wound closest to the core [33]. The magnetizing current is represented by a current source Is which is related to the linkage flux through /(s)–I(s) characteristic and can be obtained from voltage/current measurements during no-load test. In high values of linkage flux, the slope of this curve tends towards the saturated core inductance of the transformer winding. This characteristic is calculated in PSCAD/EMTDC platform based on magnetizing current at rated voltage, the air core saturated reactance of the winding, and the position of the knee point. The simulation model is illustrated in Fig. 1. 2.2. Remanence modeling During the normal operation the flux linkage will vary sinusoidally and with the same magnitude but with 120 deg displace (in time-phase) in different phases. At the time the transformer is de-energized and disconnected from the supply, the linkage flux will be frozen. Usually the maximum remanence that might exist in any leg of the core is about 80% of the peak flux linkage generated at rated volts. It is essentially important to simulate the residual flux in the core while the degree of magnetizing inrush during energizing is a direct function of the remanence existing in the legs of the transformer core. Residual fluxlinkage can be modeled by inserting a controlled DC current source in parallel with each transformer winding [32]. 2.3. Some issues which have effect on inrush current  The point on the voltage wave at the instant the transformer is energized (switching angle (hA)).  The impedance of the circuit supplying the transformer (stiffness/weakness of systems).

993

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

+

Measured Terminal Voltage

Primary

I (t )

V (t )

-

1 S

φ (t ) = v (t ).dt

I (t )

φ (t )

Fig. 1. Mutually coupled winding saturation model.

 The value and sign of the residual flux linkage in the transformer core (Br).  The nonlinear magnetic saturation characteristics of the transformer core.  Switching in the case of close or open secondary.

2.4. Some issues which have effect on fault current    

Fault type {Single-line-to-ground (SLG), Double-line-to-ground (DLG), Line-to-line (LL) and Three-phase (3PH)}. Load conditions (different active and reactive loads). Fault incipient time instant. Fault impedance.

3. Data-set collection As it is obvious from Fig. 2 which presents the main steps of algorithm implementation, the first step is training and testing sample data preparation. In this paper training and testing data set is obtained through excessive simulation executions for different operating conditions. The simulated system shown in Fig. 3, includes a typical 500 MVA, 400/230 kV, three phase power transformer. The transformer saturation model is well presented in previous section. The affect of long lines and CT’s on both side of transformer are also modeled. Distributed model for lines is used. A total of 19,000 cases for inrush current and 2960 cases for internal fault conditions are simulated considering those major parameters presented in Section 2. The input data set for GMM models is organized by implementing wavelet transform on a data window with a fixed length of a quarter of a cycle. Sampling frequency in this paper is taken to be 10 kHz. The system basic frequency is 50 Hz. A new technique is presented to enhance learning capabilities of GMM using k-means clustering algorithm [34]. The kmeans approach has been used to select training and testing data sets. Actually the algorithm is applied to the training dataset to reduce amount of samples to be presented to the GMM by automatically selecting an optional set of samples. It was just engaged in the pre-training procedure to make the training procedure simple and feasible. The clustering of data is a technique for grouping large sets of data into clusters of small sets. This technique attempts to find natural clusters of components (or data) based on some similarities, and also centroid of a group of data sets. The centroid of a cluster is a point which parameter values are the means of parameter values of all points in the cluster. The main goal is to define k centroids (k clusters), one for each cluster. For more information refer to [34]. This algorithm tries to minimize the criterion function named a squared error function defined as follows:

J¼

k X n  2 X  j  xi  cj  j¼1

ð1Þ

i¼1

where kxji  cj k2 is a chosen distance measure between a data point xji and cluster center cj. The k-means algorithm is applied to the achieved data obtained from simulations. For each class we have defined 350 clusters. The centroids of clusters have been selected as the dataset. This selection can help to have more efficient number of testing and training dataset without losing any part of space. Training and test sets are randomly selected by cross-validation approach (0.7 for training and 0.3 for testing).

994

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Start

Simulation of different kinds of fault in different operating conditions

Defining a set of centroids for each kind of fault (each class)

Clustering the data for each class separately by using k-means Training data-set selection

Testing data-set selection

Training the models

Models evaluation

End

Fig. 2. Proposed method implementation basic steps.

4. Gaussian mixture models 4.1. Definition Mixture models are types of probabilistic density models which comprise a number of component functions, usually Gaussian. These component functions are combined to provide a multimodal density. Mixture models are semi-parametric alternatives to non-parametric histograms (which can also be used as densities) and provide greater flexibility and precision in modeling the underlying statistics of sample data. These can be employed to model faults in order to perform tasks such as real-time fault identification. For an S-class pattern recognition system, we could have a set of GMMs fI1 ; I2 ; . . . ; IS g associated with S classes. A random variable x with D-dimensions is said to follow a Gaussian mixture model, if its probability density function can be formulated by (2) following the constraints presented in (3):

PðxjIÞ ¼

m X

ak Pðxjhk Þ

ð2Þ

k¼1 m X

ak ¼ 1; ak > 0

ð3Þ

k¼1

so the Gaussian mixture density is a weighted linear combination of m component Gaussian density functions, hk, P P k = 1,2, . . . m. In the previous equations ak is mixture weight for kth component, I ¼ fa1 ; . . . ; am ; l1 ; . . . ; lm ; 1 ; . . . m g is P the complete set of parameters to define the model and hk ¼ flk ; k g is the kth component mean and covariance, respectively. Each component density is a D-variant Gaussian function parameterized by a D  1 mean vector and D  D covariance matrix. The component’s density P(x|hk) is a normal probability distribution which is defined by (4):

Pðxjhk Þ ¼

1 ð2pÞ

D=2

j

P

1

j

1=2

e2ðxlk Þ

P1 k

ðxlk ÞT

ð4Þ

4.2. Training process For a GMM-based classification system, the goal of the model training is to estimate the parameters of the GMM, Ig so that the Gaussian mixture density can best match the distribution of the training feature vectors. For a set of n independent and identically distribution vectors X = {x1,x2,. . .,xn}, the likelihood corresponding to a mixture is:

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

995

Fig. 3. Simulated power system.

PðXjIÞ ¼

n Y i¼1

Pðxi jIÞ

ð5Þ

996

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

which tells the likelihood of the data X given the distribution parameters of I. Specifically, by having the distribution param^ that maximizes the likelihood: eters in I, the goal is to find I

^ ¼ arg max PðXjIÞ I

ð6Þ

I

Usually this function is not maximized directly, but the logarithm of the above probability can be calculated by:

log PðXjIÞ ¼

n X i¼1

log

m X

ak Pðxi jhk Þ

ð7Þ

k¼1

called the log-likelihood function, which is easier to calculate. The logarithm of maximum likelihood estimate of each model I cannot be solved analytically. Expectation Maximization (EM) algorithm [30] is widely used to estimate the parameters of GMM. This is an iterative algorithm that maximizes the likelihood probability generated by each GMM, PðXjIg Þ, given the data for that class. EM algorithm has two main steps: 4.2.1. E-step In this stage the posterior probability of sample xi in the tth step is computed through the following equation:

atl Pðxi jhðtÞ l Þ ðtÞ ðtÞ a Pðx i jhk Þ k¼1 k

Pðljxi ; IðtÞ Þ ¼ Pm

ð8Þ

4.2.2. M-step n o P ðtþ1Þ ðtþ1Þ Pðtþ1Þ ^ ðtþ1Þ ¼ aðtþ1Þ ; . . . ; aðtþ1Þ ; l1 ; . . . ; lm ; 1 ; . . . ðtþ1Þ This re-estimation process will continue byo replacing I instead m n m 1 P ðtÞ ðtÞ ðtÞ ðtÞ PðtÞ by using following equations: of IðtÞ ¼ a1 ; . . . ; am ; l1 ; . . . ; lm ; 1 ; . . . ðtÞ m n 1X Pðkjxi ; IðtÞ Þ n i¼1 Pn xi Pðkjxi; IðtÞ Þ lðtþ1Þ ¼ Pi¼1 k n ðtÞ i¼1 Pðkjxi; I Þ   T  Pn ðtþ1Þ ðtþ1Þ ðtÞ xi  l m xi  lm tþ1 i¼1 Pðkjxi; I Þ X ¼ Pn ðtÞ k i¼1 Pðkjxi; I Þ

akðtþ1Þ ¼

ð9Þ ð10Þ

ð11Þ

For more information refer to [35]. Note that, it is preferred to assume the S classes GMMs independent in calculations. In the case of independent classes, the estimation problem of S class pdfs can be divided into S separate estimation problems. 5. Proposed method The protective device logic is described in detail in this section. The entire system consists of data acquisition, preprocessing, feature extraction and the GMM classifier. In this method, it is required to build GMM models for each kind of fault separately. Note that, we have 2 classes (S = 2), one for internal fault and the other for inrush current. The proposed relay logic is depicted in Fig. 4. The parameters of each GMM block should be calculated off-line by training process presented in previous sections. A windowed current waveform of about quarter of a cycle (50 samples per window) is fed into wavelet transform module which finally extracts an 8 sample array. In the fault inception moment, the extracted feature signals will be used as GMM models inputs. The probability for each model is then calculated. A decision will be made by selecting the model with higher likelihood. Blocking command will be established where inrush current GMM block output is higher otherwise the signal will be identified as internal fault current and the relay trips as the result. 5.1. Wavelet feature extraction Some special features should be extracted to improve efficiency, accuracy and also reducing input samples of GMM. Also it is not practical to directly use the windowed signals as the input of GMM, because this will result in large number of samples and thereby causing difficulty in GMM training. The other point is that the fault or inrush current may occur only in one or two phases, as the result, it is essentially important to use three phase differential currents as the input of protection device. Since each GMM block can have only one input at the time, three dimension current signals should be converted to one dimension. The wavelet transform has been introduced recently as a relatively new and powerful tool in the analysis of transient phenomena in power systems. It has the ability to extract information from the transient signals simultaneously in both time and frequency domains. The ability of wavelet transform to focus on short time intervals for high frequency components and long intervals for low-frequency components improves the analysis of signals with localized impulses and

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

997

Three phase differential currents (I ad, Ibd, Icd) from CT’s

Selected window (5 ms) Wavelet Transform to the windowed signals (I ad, Ibd, Icd) Feature Extraction

Inrush current GMM model (I)

Restrain command

Internal fault GMM model (II)

Comparator (Maximum Probability)

Internal fault?

Yes Trip

No Fig. 4. Differential protection logic.

oscillations. For this reason wavelet decomposition is ideal for studying transient signals and obtaining a better current characterization and more reliable discrimination technique [5]. Discrete wavelet transform signal spectral energy has been used in this paper. In wavelet transform analysis, the choice of mother wavelet plays a significant role in detecting and localizing different kinds of fault transients. A specific mother wavelet is commonly best suited for a particular application. For this purpose, mother wavelet type and decomposition level have been chosen based on experience and try and error. A certain high frequency component can be located better in time than a low-frequency component [8]. In this study, since we are interested in these components which are located better in time, higher frequency details (1–4) are employed for analysis and feature extraction. Besides, the research includes detecting and analyzing low amplitude, short duration, fast decaying, and oscillating type of high frequency current signals. For this purpose, Daubichies’s mother wavelet seems to be an appropriate choice [8]. In comparison with Harr wavelet, Daubichies are best suited for feature extraction due to their low-pass and high-pass filters. On the other hand because of its orthogonality inherent, it satisfies Parsaval theorem, not like biorthogonal wavelets such as Coiflet and Meyer wavelets [5]. db7 mother wavelet over level d6 have been chosen because the maximum energy localization in details (1–4) was obtained using these parameters. The selected data window is divided into 2 equal time periods. Detail components are extracted by DWT decomposition for each half as it is clearly shown in Fig. 5. Then, the spectral energies (SE) of details (1–4) are calculated through Eq. (12). Hence, 8 samples are obtained from details 1–4 for each phase. Finally the feature array presented in Eq. (13) is considered as the main feature. The array consists of the average of three phase extracted features. The transformation can retain principle features of the signal.

SE ¼

F X

jdðiÞj2

ð12Þ

i¼1

  x ¼ e11 ; e12 ; e13 ; e14 ; e24 ; e23 ; e22 ; e21

ð13Þ

where F is the number of samples in a specific detail (d). In the component ehd , e is the average of energy over three phases, h determines the first or second half and d is the detail number. The motivation of selecting this feature is that, the amplitude of the wavelet details in d1–d4 is larger than other low-frequency components. It should be mentioned that in high frequency details the first half signal’s spectral energy for fault current cases is greater than the second half and visa versa for inrush currents. This trend can also be concluded from simple principles presented in [3]:  The differential current due to the fault begins in higher slope and then its slope decreases. But differential current due to inrush current begins with a low slope and then its slope is raised.  A larger slope in the time domain shows that there are higher frequencies in the frequency domain.

998

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Second half of signal

First half of signal

a11

a12

a13

a14

a15

a16

d 12

d31

d 14

d51

d61

d11

a22

a32

a42

a52

a62

a12

d12

d 22

d32

d 42

d52

d 62

Fig. 5. Wavelet decomposition over 6 levels.

5.2. GMM classification In this stage, the objective is to find the model which has the maximum probability Pmax for a given observation sequence. The test patterns are provided for each of the GMM models ðP I1 ; P I2 ; . . . ; PIS Þ, and S probabilities are computed and compared, i.e. Pi, i = 1,2,. . .,S. The class f associated with the GMM that has the highest Gaussian mixture likelihood probability Pmax is chosen as the predicted class for the given feature set. Assume that we have a feature vector x produced by one type of faults. The identity of the fault producing is determined by finding the fault model If from S models which maximizes the probabilities across the fault set Ig : g 2 ð1; SÞ:

If ¼ arg max PðIg jxÞ

ð14Þ

16g6S

Using the Bayes rule (14) can be expressed as:

If ¼ arg max 16g6S

PðxjIg ÞPðIg Þ PðxÞ

ð15Þ

Assuming each fault model is equally the same and noting that P(x) is the same for all models, the identification task can be summarized as finding the logarithm of the following equation:

If ¼ arg max PðxjIg Þ

ð16Þ

16g6S

For simplicity, we investigate the logarithm of the above equation:

" If ¼ arg max log 16g6S

m X

#

akg Pðxjhkg Þ

ð17Þ

k¼1

where akg and hkg are mixture weights, means, and covariance’s of the gth fault model. 6. Proposed method evaluation 6.1. GMM parameters In order to achieve better performance of Gaussian mixture model an effort to investigate an optimum configuration of algorithmic issues must considered [25]. These issues include convergence threshold, model order selection, and the form of the covariance matrix: 6.1.1. Convergence threshold of EM algorithm The convergence threshold of the EM algorithm is defined as the difference of the probability PðXjIg Þ between iterations. The convergence threshold and the maximum iteration number are two conditions for stopping the EM algorithm. Different

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

999

values of the convergence threshold were carefully examined and a value of 0.0005 was found to be the best for achieving good performance. 6.1.2. GMM model order The GMM order selection is the most crucial factor in the performance of this classifier [25]. The objective is to choose the number of mixture components that yields the highest discriminative accuracy for the test dataset. Theoretically, too few mixture components can yield a GMM which does not accurately model the distinguishing characteristics of classes. However, too many components can also reduce performance, while increase the computation burden. It also makes complexity in classification. This is especially important for situations with smaller amount of training data [25]. A larger model order should be considered when larger amounts of training data are available. Therefore, the optimal number of mixture components for the GMM must be carefully examined to achieve high classification performance. The optimal number of mixture is achieved by computation of classification rate presented in (18) for test data set and varying number of mixtures. As it is obvious from Fig. 6, a Gaussian model with 5 mixtures is best suited for our problem.

Classification rate ¼

Correct classified transient events Total transient events

ð18Þ

6.1.3. The form of the covariance matrix Another effective issue on the performance of the classifier is form of covariance matrix. The form of the covariance matrix may be full or diagonal. In this work, diagonal covariance matrices are selected. This choice is based on empirical evidence that diagonal matrices outperform full matrices. The diagonal matrix GMMs are more computationally efficient than full covariance GMMs for training since repeated inversions of a D  D matrix are not required [25]. 6.2. Simulation results This scheme is tested using relaying signals obtained from simulation of various operating conditions in Section 3. Some typical results are well presented in Tables 1 and 2. Table 1 shows the probabilities of two kinds of fault models in inrush current cases. Results show that the proposed method, can effectively distinct transient phenomena. Table 2 shows the output of GMM models in the case of internal fault. In addition to the fault and inrush currents cases, it is necessary to study the merit of the algorithm for simultaneous fault and switching condition. Sympathetic inrush current and recovery inrush currents are also simulated and evaluated in this section. The winding fault identification is also very important case for power transformers which is considered as the contribution of this article. 6.2.1. Case 1: Simultaneous inrush and internal fault Fig. 7 shows three phase current waveforms under the case of simultaneous fault and magnetizing inrush current at t = 0.1016 s. the residual fluxes are taken Bra = 60%, Brb = 60%, Brc = 80% and the voltage switching angle is hA = 30. The log-

Fig. 6. Classification rate versus number of mixtures.

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

1000

Table 1 Logarithm values of models probabilities for different inrush current cases. Br

Fault model

hA

Bra = 20%; Brb = 60%; Brc = 80%

0

Inrush Internal Inrush Internal Inrush Internal Inrush Internal Inrush Internal Inrush Internal Inrush Internal Inrush Internal Inrush Internal Inrush Internal

60 90 120 270 Bra = 80%; Brb = 0%; Brc = 80%

0 60 90 120 270

No-load

On-load

Stiff system

Weak system

Stiff system

Weak system

11.5 89.4 8.7 56 30.3 98.1 34.3 55.1 24.5 35.8 17.6 140 19.8 105.3 6.9 67.3 16.7 55.3 28 41.7

29 76.1 31.2 49.3 17.1 25.9 15.4 68 19.9 40.6 42 82.1 21.7 52.6 15.7 35.4 47.3 61.5 31.7 57.1

18.3 93.2 31.8 54.6 9.1 11.3 24.2 98.1 31.7 44.8 15.7 152 36 66.3 14 32.1 17.8 107.9 48.3 65.6

26.9 85.5 43 64.3 13.4 44.9 12.5 63.2 13.2 24.9 29.5 163 5.7 123.6 14.9 28.4 8.9 20.4 11.3 21.5

Table 2 Logarithm values of models probabilities for different internal fault cases. hA

GMM Model

0

Inrush Internal Inrush Internal Inrush Internal

60 90

No-load

On-load

a–g

a–b

a–b–g

a–b–c

a–g

a–b

a–b–g

a–b–c

130.4 21.2 158.7 10.5 195.8 17.6

155.3 45.1 203.3 20.6 125.7 12.7

91.6 13 139.8 31.1 115.8 25.8

181.5 25.7 148.3 19.3 167.9 11.6

226.3 29.8 95.3 14.7 228.9 18

119.9 17.3 130.1 8.8 241.5 36.4

212.38 12.5 117.9 15 93.46 28.1

110.8 22.7 99.9 34.1 89.6 13.7

arithm value of model probabilities for inrush and fault are 67.8 and 32.6, respectively. The comparison between GMM models outputs represents the fault occurrence. 6.2.2. Case 2: Sympathetic inrush current This case analyzes the effect of the inrush current on the transformers that are already in operation. In practice, transformers are normally energized in parallel with other transformers that are already in operation. Assume a system which includes a generator connected to a transformer bus through a transmission line having a resistance R and an inductance L. Two transformers are potentially installed to feed the loads. Assume that one of the transformers is loaded or switched; the other one will be switched in parallel, with a delay. As the breaker closes, an inrush current is established in the primary winding of the second transformer. The inrush current has a decay dc component which produces a voltage drop in the resistance of the transmission line. The dc voltage drop has a specific polarity. Since the generator output is purely ac, and cannot be affected by this voltage drop, the transformer busbar acquires a negative dc component. This results in a negative change in the flux linkages of the two transformer cores. If the second transformer has a saturating flux in the positive direction, the effect of flux changes takes it out of saturation, and causes a possible saturation of the first transformer in the negative direction. Consequently, the inrush current in this transformer increases [36]. This is called sympathetic inrush current. To simulate this case, two 100 MVA, 400/230 kV transformers have been installed in the simulated system. The first transformer is switched in t = 0.1016 s, hA = 30, consequently the second one is switched in t = 0.308 s, hA = 60. The remnant fluxes in transformers have been taken B1ra = 80, B1rb = 80, B1rc = 80 and B2ra = 80, B2rb = 80, B2rc = 80, respectively. Three phase differential current and the achieved feature signal are shown in Fig. 8. The differential currents magnified in time interval of 0.2 s to 0.4 s are also well presented in Fig. 9. The logarithm value of model probabilities for inrush and fault are 31.8 and 57.6, respectively. 6.2.3. Case 3: Recovery inrush current A magnetizing inrush current can also flow if a voltage dip is followed by recovery to normal voltage. Typically this occurs upon removal of an external fault. The magnetizing inrush is usually less severe in this case than in initial energization as the

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

1001

Fig. 7. Three phase differential currents and feature signal as a result of simultaneous BC-G fault and inrush current.

Fig. 8. Three phase differential currents Ia, Ib, Ic and the related feature signal in the case of sympathetic inrush current.

transformer was not totally de-energized prior to voltage recovery. To simulate this case, a three phase fault is occurred in the primary side of transformer but out of differential relay protection area. The fault is simulated in time interval between 0.1 s and 0.905 s. The results are presented in Fig. 10. The differential currents start to flow when the fault is cleared. The differential currents are also presented with magnification in Fig. 11. The relay also has performed well as expected. The logarithm value of model probabilities for inrush and fault are 17.3 and 86.9, respectively.

1002

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Fig. 9. Three phase differential currents in the case of sympathetic inrush current (magnified in time range of 0.2 s to 0.4 s).

Fig. 10. Three phase differential currents Ia, Ib, Ic and the related feature signal in the case of recovery inrush current.

6.2.4. Case 4: Inter-turn fault current To investigate the winding faults, detail modeling of transformer is essentially required. A 132 kV/11.5 kV, 20 MVA and Y– Y connection transformer (with neutral connection) and inhomogeneous HV winding is chosen. HV and LV windings of each phase of transformer have been simulated according to the models presented in [37,38]. All modeling and simulation have been carried out in Matlab/Simulink software. Fig. 12 shows equivalent circuit for one phase (HV + LV) of transformer [39]. In

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

1003

Fig. 11. Three phase differential currents in the case of recovery inrush current (magnified in time range of 0.65 s to 0.85 s).

Fig. 12. Equivalent circuit for detailed transformer model.

order to calculate the parameters of equivalent circuit the HV winding is divided to eight sections and LV winding to two sections. For more information please refer to [37–39]. In this case a short circuit has been simulated between Sections 5 and 6 of HV winding in phase A (at t = 0.1 s). The differential currents are well presented in Fig. 13. As expected, the comparison of output probabilities could well identify the fault in the protection zone (24.1 for internal and 95.3 for inrush current models).

1004

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Fig. 13. Three phase differential currents and the related feature signal in the case of inter-turn fault.

6.2.5. Case 5: Turn to ground fault current The presented model in case 4 has also been used to simulate the turn to ground fault. This fault takes place between Section 6 of B phase of HV winding and the transformer tank (25% of winding turns have been short circuited). The inception time is t = 0.1 s. Three phase differential currents are shown in Fig. 14. The output probabilities of inrush and internal fault models are 67.3 and 14.8, respectively. 6.3. Experimental test setup The proposed GMM-based method has been also tested on a set of experimental data in order to validate the merit of the algorithm. The test setup consists of a three phase, 5 kVA, 230/550–575–600, 60-Hz, core type, D–Y-connected laboratory prototype power transformer, three identical CT’s connected in Y in the primary side and another three identical CT’s connected in D on the secondary side. The experimental setup is demonstrated in Fig. 15 [7]. The differential current is measured at point Q using Tektronics current probe. The differential current is sampled at 10 kHz and stored using Tektronics 2212 digital storage oscilloscope [7]. Three identical Triac switches and a control circuit to trigger them are used to provide a short contact between the supply and the transformer. The configuration of the switches and the control circuit are well presented in [7]. The following four cases are considered to evaluate the proposed method performance: 6.3.1. Case 6: Initial or switching magnetizing inrush current The most common type of inrush current is initial inrush current. Differential magnetizing inrush current is collected while the transformer is unloaded and the transformer is switched to the primary side. The supply voltage line is to be 230 V. Fig. 16 shows the three phase differential current waveforms of Ia, Ib and Ic and the feature signal. As seen the feature signal has an increasing trend. The logarithm of the GMM’s output probabilities are 25.4 and 53 for inrush and fault, respectively. The comparison between GMM models outputs implies the magnetizing inrush current occurrence (there is no fault and the maltrip is not issued). 6.3.2. Case 7: Single-phase-to ground fault current In the case of single line-to-ground, the transformer is loaded with a balanced Y-connected load. It should be noted that the primary line-to-line voltage is set to 50 V in all investigated fault cases in order to avoid the saturation or damage in the tested equipment. Fault is occurred in the secondary side (phase C) of transformer. Fig. 17 shows differential currents and also the feature signal obtained from wavelet transform. The feature has a decreasing trend which shows the occurrence of fault. The logarithm of output probabilities are 110.8 and 19.6 for inrush and fault models, respectively.

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

1005

Fig. 14. Three phase differential currents and the related feature signal in the case of turn to ground fault.

Fig. 15. Experimental setup for differential current data collection.

6.3.3. Case 8: Line-to-line fault current Fig. 18 shows differential currents for the case of unloaded secondary and line-to-line fault in the primary side of transformer (phases B and C). The criterion is –61.7 for inrush and 28.4 for fault models. The method effectively detects fault condition while the outputs have significant difference in their values.

1006

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Fig. 16. Three phase currents and feature signal for the case of unloaded transformer switching.

Fig. 17. Three phase currents and feature signal for the case of secondary loaded single-phase-to-ground fault.

6.3.4. Case 9: Three-phase-to-ground fault current The three-phase fault has occurred on the primary side of the unload power transformer. Fig. 19 shows the three phase differential currents Ia, Ib, Ic and the feature signal. The logarithm of output probabilities are 192.7 and 15.2 for inrush and fault models, respectively. The comparison between GMM models outputs represents the fault occurrence.

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

1007

Fig. 18. Three phase currents and feature signal for the case of primary side phase-to-phase fault.

Fig. 19. Three phase currents and feature signal for the case of primary side three-phase to ground fault.

6.3.5. Case 10: Effect of noise To investigate the noise effect on discrimination accuracy, the measured differential currents from simulations have been polluted to white Gaussian noise with two different signal to noise ratio (SNR). Three cases have been considered (one inrush

1008

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Table 3 The effect of noise on proposed method functionality. SNR = 20 Inrush current model Inrush current, Bra = 20%, Brb = 60%, Brc = 80%, hA = 0, No-load, Weak Sys. Internal fault, ag, hA = 60, No-load Internal fault, abc, hA = 90, On-load

SNR = 10 Internal fault model

Inrush current model

Internal fault model

43.1

52.6

91.3

106.7

185.6 143.7

61.5 45.7

201.8 179.13

83.2 135.9

current and two fault cases). The simulation results are presented in Table. 3. All simulation results show the robustness of the wavelet based Gaussian mixture model classifier. Although the discrimination accuracy have been reduced a little but the proposed method could also remain stable in these harsh situations. 6.4. Discussion From the results we can judge that the combinatorial scheme performs well in both accuracy and speed of response points of view. It has the probabilistic advantage of GMM while using time-frequency dependent features extracted by wavelet. The use of DWT naturally emphasizes the difference between fault and inrush currents since their time-frequency features are quite different. The proposed GMM based method recognizes different operating conditions of power transformers accurately and within 5 ms. Since the discrimination method is performed with probabilistic characteristics of signals without application of any deterministic index, high reliable and accurate classification is achieved. It provides a high operating sensitivity for internal faults and remains stable for inrush currents of the power transformers. Many novel methods have been enhanced in recent years but they usually suffer from lack of two main advantages simultaneously (accuracy and speed of response). To validate the proposed method’s merit, it is compared in both accuracy and speed of response aspects with the methods presented in [8,40]. The proposed scheme in [8] uses a combinatorial approach based on Artificial Neural Network (ANN) and DWT. Classification rate is reported 99% in [8] whereas we have been achieved 100% recognition rate. The input data set for ANN in [8] is organized by a data window with a length of 10 ms while we have been used a 5 ms window (faster differential relaying). The presented method in [40] has reported two different neural network structures: (1) Feed forward back propagation neural network, (2) probabilistic neural network. The recognition rate related to the methods has been reported 97% and 100%, respectively. But the drawback of these methods in comparison with GMM based scheme stems from the identification speed (larger input data requirement of about 20 ms). Although ANN can classify patterns perfectly but the results do not determine the classification accuracy. In contrast, GMM classifies patterns by using two probabilities. The variance between two probabilities could be used as a criterion in which higher extent in difference shows higher accuracy. The large amount of difference between two probabilities, imply the superiority of GMM. With carefully trained GMM, it has some advantages such as:  No preset threshold is needed.  Simple identification criterion.  It is best suited for a large power system protection which is full of uncertainty parameters because of the probabilistic manner.  The possibility to detect faults in a quarter of a cycle.

7. Conclusion This paper has introduced and evaluated the use of GMM’s for discrimination between internal faults and inrush currents. The problem has been broken down into the tasks of data collection, feature extraction, classification, and performance evaluation. Required data set has been achieved through excessive simulations in PSCAD/EMTDC software. Wavelet transform has been applied to extract discriminant features. The k-means clustering approach was applied to dataset to reduce training computation burden and possibility. Discrimination was carried out by computation of GMM models output probabilities. Successful operation of the trained model by using an experimental setup, which is completely different from training setup, proves the proposed method independency on system configurations, operating conditions and transformers. The proposed algorithm can be applied for about a quarter of a cycle (fast discrimination for on-line applications). Results show the stability of protection technique in presence of the magnetizing inrush current with robust fault diagnosis even in presence of noise which may insure secure and optimal operation of protection system. The high accuracy and fast performance of this method suggests more investments in application of GMMs in power systems.

1009

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

Appendix A.

Power transformer magnetization curve:

Current (% of rated current)

Voltage (pu.)

0 0.17 0.48 0.98 2 3.09 6.52 20.35 60.21 124.38

0 0.32 0.61 0.82 1 1.08 1.17 1.26 1.36 1.49

Leakage reactance: 0.1 pu. Knee point voltage: 1.2 pu.

List of mathematical symbols: xji cj J S D Ii

A data point in the array Cluster center Criterion function of k-means algorithm Number of classes (models) The dimension of variable xi The ith model A set of data vectors xi

X

xi m n hk

ak l Pk

k

^ I PðXjIÞ Dt e P Ii Pmax

The input vector of GMM (feature vector) The number of mixtures The dimension of X (number of data vectors xi) Gaussian density function The mixture weights The mean of Gaussian density function The covariance of Gaussian density function The re-estimated model I The probability of X given model I An small time interval The average of spectral energy of Wavelet Transform The output probability of model i The maximum probability of models

References [1] M.A. Rahman, B. Jeyasurya, A state of the art review of transformer protection algorithms, IEEE Trans. Power Delivery 3 (1988) 534–544. [2] H. Zhang, P. Liu, O.P. Malik, A new scheme for inrush identification in transformer protection, EPSR 63 (2002) 81–86. [3] J. Faiz, S. Lotfi-Fard, A novel wavelet based algorithm for discrimination of internal faults from magnetizing inrush currents in power transformers, IEEE Trans. Power Delivery 21 (2006) 1989–1996. [4] A.R. Seddighi, M.R. Haghifam, Detection of inrush current in distribution transformer using wavelet transform, EPES 27 (2005) 361–370. [5] H. Monsef, S. Lotfifard, Internal fault current identification based on wavelet transform in power transformers, EPSR 77 (2007) 1637–1645. [6] O.A.S. Youssef, Discrimination between faults and magnetizing inrush currents in transformers based on wavelet transforms, EPSR 63 (2002) 87–94. [7] S.A. Saleh, M.A. Rahman, Modeling and protection of a three phase power transformer using wavelet packet transform, IEEE Trans. Power Delivery 20 (2005) 1273–1282. [8] P.L. Mao, R.K. Aggarwal, A novel approach to the classification of the transient phenomena in power transformers using combined wavelet transform and neural network, IEEE Trans. Power Delivery 16 (2001) 654–660. [9] S. Sudha, A.E. Jeyakumar, Wavelet and ANN based relaying for power transformer protection, JCS 3 (6) (2007) 454–460. [10] M. Geethanjali, S.M.R. Slochanal, R. Bhavani, PSO trained ANN-based differential protection scheme for power transformers, Neurocomputing 71 (2007) 904–918. [11] M. Jing, W. Zengping, A novel algorithm for discrimination between inrush currents and internal faults based on equivalent instantaneous leakage inductance, power engineering society general meeting, IEEE (2007) 1–8. [12] G. Baoming, A.T. Almeida, Z. Qionglin, W. Xiangheng, An equivalent instantaneous inductance-based technique for discrimination between inrush current and internal faults in power transformers, IEEE Trans. Power Delivery 20 (2005) 2473–2482. [13] M.E.H. Golshan, M. Saghaian-nejad, A. Saha, H. Samet, A new method for recognizing internal faults from inrush current conditions in digital differential protection of power transformers, EPSR 71 (2004) 61–71.

1010

S. Jazebi et al. / Simulation Modelling Practice and Theory 17 (2009) 991–1010

[14] Y.C. Kang, B.E. Lee, S.H. Kang, Transformer protection relay base on induced voltages, EPES 29 (2007) 281–289. [15] H. Zhang, J.F. Wen, P. Liu, O.P. Malik, Discrimination between fault and magnetizing inrush current in transformers using short time correlation transform, EPES 24 (2002) 557–562. [16] S. Guocai, Y. Dachuan, Identifying internal faults of transformers through the similarity degree between voltage and current, power engineering society winter meeting, IEEE 3 (2000) 1868–1872. [17] E. Sayed, M.T. Eldin, A novel approach for classifying transient phenomena in power transformers, IJEEPS 1 (2004) 1010. Article. [18] T. Babnik, F. Gubina, Fast power transformer fault classification methods based on protection signals, IEE Proc. Gener. Transm. Distrib. 1 (2003) 205– 210. [19] S.R. Samantaray, B.K. Panigrahi, P.K. Dash, G. Panda, Power transformer protection using S-transform with complex window and pattern recognition approach, IET Gener. Transm. Distrib. 1 (2007) 278–286. [20] G. Diaz, P. Arboleya, A new transformer differential protection approach on the basis of space-vectors examination, EEJ 87 (2005) 129–135. [21] Z. Moravej, Power transformer protection using support vector machine network, PES (2008). [22] X. Ma, J. Shi, A new method for discrimination between fault and magnetizing inrush current using HMM, EPSR 56 (2000) 43–49. [23] S. Jazebi, B. Vahidi, S.H. Hosseinian, A new stochastic method based on hidden markov models to transformer differential protection, IEEE Proc. OPTIM’08 1 (2008) 179–184. [24] B. Xiang, T. Berger, Efficient text-independent speaker verification with structural gaussian mixture models and neural network, IEEE Trans. Speech Audio Proc. 11 (2003) 447–456. [25] Y. Huang, K.B. Englehart, B. Hudgins, A.D.C. Chan, A Gaussian mixture model based classification scheme for myoelectric control of powered upper limb prostheses, IEEE Trans. Biomed. Eng. 52 (2005) 1801–1811. [26] X. Zhou, X. Wang, Optimisation of Gaussian mixture model for satellite image classification, IEE Proc. Vis. Image Signal Process. 153 (2006) 349–356. [27] M. Hamouz, J. Kittler, J.K. Kamarainen, P. Paalanen, H. Kälviäinen, J. Matas, Feature-based affine-invariant detection and localization of faces, IEEE Trans. Pattern Anal. Mach. Intell. 27 (2005) 1490–1495. [28] J.C. Principe, M. Motter, Identification and Control of Aircrafts using Multiple Models and Adaptive Critics, NASA Proposal No. 04-0617, 2007. [29] A. Ruano, Intelligent control systems using computational intelligence techniques, IEE Control Series (2005). [30] G. Xiong, C. Feng, L. Ji, Dynamical Gaussian mixture model for tracking elliptical living objects, Pattern Recogn. Lett. 27 (2006) 838–842. [31] T. Marwala, U. Mahola, F.V. Nelwamondo, Hidden Markov models and Gaussian mixture models for bearing fault detection using fractals, 2006 Int. Joint Conf. Neural Networks (2006) 3237–3242. [32] D. Woodford, Introduction to PSCAD/EMTDC V3 (Manitoba HVDC Research Center Inc. (2000). [33] H.W. Dommel, Transformer Models in the Simulation of Electromagnetic Transients, in: Proceedings of Fifth Power Systems Computing Conference, England, Paper 3.1/4, 1975. [34] K.M. Faraoun, A. Boukelif, Neural networks learning improvement using the k-means clustering algorithm to detect network intrusions, IJCI 5 (2006) 28–36. [35] P. Paalanen, J.K. Kamarainen, J. Ilonen, H. Kälviäinen, Feature representation and discrimination based on Gaussian mixture model probability densities — Practices and algorithms, Pattern Recogn. 39 (2006) 1346–1358. [36] M. Sengül, B. Alboyaci, S. Öztürk, H.B. Cetinkaya, Case study of sympathetic interaction between transformers caused by inrush transients, in: International Conference on Power Systems Transients, Paper No. IPST05-125, 2005. [37] G.B. Gharepetian, H. Mohseni, K. Moller, Hybrid modeling of inhomogenous transformer windings for very fast transient overvoltage studies, IEEE Trans. Power Delivery 13 (1998) 157–163. [38] S.M.H. Hosseini, M. Vakilian, G.B. Gharepetian, Comparison of transformer detailed models for fast and very fast transient studies, IEEE Trans. Power Delivery 23 (2008) 733–741. [39] A. Shirvani Broujani, Identify the mechanical displacement of HV winding of transformer with the aid of detailed model, M.Sc. Thesis, Amirkabir University of Thechnology, Tehran, Iran, 2007. [40] M. Tripathy, R.P. Maheshvari, H.K. Verma, Probabilistic neural-network-based protection of power transformer, IET Elect. Power Appl. 1 (5) (2007) 793–798.