An EKF‐SVM machine learning‐based approach for fault detection and classification in three‐phase power transformers

Correspondence Farshid Naseri, Department of Control and Power Engineering, Shiraz University, Shiraz, Iran. Email: f.naseri@shirazu.ac.ir Abstract In this paper, a hybrid approach for effective diagnosis of power transformers is proposed. In the proposed method, the extended Kalman filter is used for the estimation of three-phase currents in the primary windings of the transformer. Three residual signals are defined as the differences between the measured and estimated three-phase currents. When the transformer is healthy, the EKF perfectly estimates the primary currents and hence, the residual signals are nearly zero. However, when the transformer is faulty, the EKF cannot suitably estimate the currents due to the large model mismatch resulting from the transformer internal faults. Consequently, large residual signals are generated, which are used as the key signatures for discriminating the internal faults from energisation conditions. Besides, the proposed method uses the entries of the covariance matrix of the estimation error to locate and classify the internal faults. For these purposes, support vector machine classifiers are used. The effectiveness of the proposed method is demonstrated by a large number of simulation test cases obtained using PSCAD/EMTDC software. Also, hardware-in-the-loop experiments are conducted using dSPACE1104 development platform and real-time feasibility of the proposed method is authenticated.

Traditionally, harmonic-based algorithms such as harmonic sharing or cross-blocking have been used to resolve the abovementioned problem [2]. The harmonic-based methods use this fact that the second harmonic component in the inrush current is considerably large. However, harmonic restraint methods are not dependable for the protection of modern transformers, which possess high-quality material and experience low second harmonic inrush events [3]. A large number of approaches have been proposed to rectify the deficiencies of the harmonic restraint methods, which include the use of wavelet transform [4], S-transform [5], artificial neural networks (ANNs) [6], fuzzy logic [7], Park's transform [8], higher-order statistics [9], extended Kalman filter (EKF) [13] etc. [10][11][12]. Despite the high performance of these approaches, in none of them, the classification of electrical internal faults has been fulfilled. In this sense, the researches that have so far been conducted mostly concern the transformer mechanical defects. However, in addition to the discrimination of the internal faults, the classification of the electrical internal faults is also very important due to the following reasons: • The downtime of the transformer will be substantially reduced by expediting the maintenance and repair processes. • The reliability and flexibility of the protection scheme will be enhanced.
In this paper, the foregoing task is fulfilled using a hybrid approach. In the proposed method, the electrical internal faults are first discriminated using an EKF. The EKF-based algorithm has recently been proposed for single-phase transformers by the authors of this paper in [13]. In this paper, the EKF-based approach is extended to adapt to the three-phase transformer case. Next, an innovative feature vector that is obtained from the prior stage is introduced to a set of support vector machine (SVM) classifiers for determining valid information regarding the type, location, and severity of the electrical internal faults. The usefulness of the KF algorithm and machine learning techniques such as SVM and ANN classifiers have been previously proven in earlier power system protection studies with various great achievements [14,15]. In this paper, these tools are used in the diagnosis of power transformers.
The rest of the paper is organised as follows: In Section 2, the state-space model of the under-study three-phase power transformer and EKF formulation are developed. In Section 3, the working principles of the proposed method are presented. To demonstrate the validity of the proposed approach, extensive simulation results are presented and discussed in Section 4. In Section 5, a comparison between the proposed method and some other algorithms is provided. To evaluate the robustness of the proposed EKF-based approach against measurement errors, a robustness analysis is also conducted and relevant results are included in Section 6. The feasibility of the embedded implementation is demonstrated in Section 7. Eventually, in Section 8, the main results are concluded. An Appendix provides some details about the studied power transformer.

MODEL DERIVATION AND EKF FORMULATION
In this section, the nonlinear state-space model of the understudy transformer is developed. Then, the formulation of the EKF is presented. The detailed description in these regards is provided in the following subsections:

Nonlinear state-space equations of the transformer
In this paper, a three-phase five-limb transformer with a wyewye connection is adopted for the study. Since the five-limb topology is the general case, the proposed method can be easily applied to three-limb and single-phase transformer topologies by some modifications. The topology of the five-limb transformer is depicted in Figure 1(a). The equivalent electric circuit of the transformer is also shown in Figure 1(b). Likewise, the equivalent magnetic circuit of the transformer is shown in Figure 2. The magnetic and electric circuits are merged to obtain FIGURE 1 (a) Topology of the five-limb three-phase transformer studied in this paper. (b) Equivalent electric circuit of the transformer FIGURE 2 Equivalent magnetic circuit of the under-study three-phase transformer a unified electromagnetic model, which considers the magnetic saturation, hysteresis, and electrical connections effects [16]. From Figure 1(b), the derivative of flux linkages i can be obtained by applying Kirchhoff 's voltage law (KVL) in each winding loop as follows: (1) In (1), v 1 , v 2 , and v 3 are the three-phase voltages of the primary windings considered as the system inputs, i 1 , i 2 , and i 3 are the currents at the primary windings of the transformer considered as the system outputs, and j ( j = 1, 2, … , 6) are flux linkages, which are considered as the system states. Likewise, i 4 , i 5 , and i 6 are the currents at the secondary windings of the transformer, r L is the load resistance, and r 1 to r 6 are the specific resistances of the transformer windings. Equation (1) can be rewritten in a matrix form as follows: where ∈ matbbR 6×1 , L ∈ ℝ 6×6 and i ∈ ℝ 6×1 are the vector of flux linkages, inductance matrix, and current vector, respectively. Equations (1) and (2) should be manipulated to obtain the well-known form of the state-space equations as follows: where x is the state vector, y is the output vector, u is the input vector, and A, B, C, and D are the system matrices. Likewise, and are the process and measurement noises with covariance matrices Q and S, respectively. The state matrices A, B, C, and D need to be derived to form the state-space model. Assuming a three-phase balanced load with resistance r L being connected to the secondary windings, the state-space equations can be obtained as (4). (4B) In (4), Γ is inverse of L and is changing according to the nonlinear characteristics of the magnetic material. Likewise, I ∈ ℝ 3×3 and 0 ∈ ℝ 3×3 are identity and zero matrices, respectively. From (4), it can be deduced that A and C matrices are also changing. Therefore, based on the current operating point of the system, L must be updated. In this paper, the Newton technique with Numerical Differentiation procedure is used to calculate L in a recursive manner [16]. Considering the transformer being energised in the no-load condition, one can write i 4 = i 5 = i 6 = 0. Figure 2 shows the equivalent magnetic circuit of Figure 1(a). In Figure 2,  j ( j = 1, 2, 3) are the magnetomotive forces produced by the primary windings of the transformer and denotes the magnetic flux in different core sections. Moreover,  1 to  7 are nonlinear reluctances of different sections including windings branches, yoks, and outer limbs, which are changing according to the saturation characteristics of the magnetic material. These reluctances can be expressed as follows: where l j and A j are the length and cross-section area of the core, respectively. Moreover, j is the permeability of the j th branch of the magnetic core. Since L is changing, j must be updated in each iteration using the incoming sample points. To fulfil the foregoing task, the equivalent magnetic circuit of Figure 2 can be analysed. Accordingly, the following equations are derived: Solving (6) and using the expressions 1 = 2 − 1 , 2 = 3 − 2 , 3 = 4 − 3 , 4 = − 2 , 5 = 3 , 6 = − 1 , and 7 = 4 , the values of the flux linkages can be calculated. To update the permeability, the magnetic flux intensity H should also be calculated for each branch of the magnetic core using Ampere's law ∑ H j l j =  . The magnetic flux intensity in each winding branch H 1 , H 2 , and H 3 can be calculated using the following formula: The magnetic flux intensity of the yokes can be calculated as follows: The magnetic flux intensity of H 6 and H 7 for outer limbs can also be obtained as follows: After the calculation of H, the magnetic flux density B can be obtained using the B-H relationship of the transformer. Several approaches have been proposed for modelling the nonlinear behaviour of the B-H magnetisation curve such as using polynomial, tangent, and hyperbolic tangent functions. In this paper, the sum of the arctangent and tangent hyperbolic functions is fitted to the data of the magnetisation curve to fulfil the foregoing objective. The curve fitting is accomplished using nonlinear least-squares optimisation and the Levenberg-Marquardt algorithm. The root mean square of errors (RMSE) for the fitted curve is 0.1351. The real and estimated magnetisation curves of the transformer are shown in Figure 3. Finally, the permeability is updated as follows: The above procedure is repeated for each sample point. Therefore, the state-space model of the three-phase transformer is achieved using the explained iterative manner. The EKF algorithm is discussed in the next subsection.

Formulation of the EKF
The algorithm of the EKF is explained through the flowchart of Figure 4 [17][18][19][20], where K denotes the Kalman filter gain and P is the covariance matrix of the estimation error. In this flowchart, A, B, and C matrices can be obtained using (4). However, the partial derivative matrices G and M and the covariance matrices of the process and measurement noises Q ∈ ℝ 6×6 and S ∈ ℝ 3×3 need to be calculated. The covariance matrix of the measurement noise S is selected based on the typical accuracy class of the protective instrument transformers (based on standard IEC 61869-1). Likewise, the process noise is tuned using the trial and error approach, which is the method adopted by the majority of the existing works [3,13]. Accordingly, these matrices are selected as follows: where I ∈ ℝ 3×3 is the identity matrix. Assuming the model and observation (measurement) noises are applied to the state-space model in an additive form [21][22][23][24], G ∈ ℝ 6×6 and M ∈ ℝ 3×3 matrices can also be written as follows: Since it is assumed that the observation (measurement) noises are applied to the state-space model in an additive form, no ill-conditioning problem in the matrix inversion process related to the Kalman gain calculation will happen in this work.
In Figure 4, T is the vector of the three-phase input voltages that are measured at each sample. Moreover, y = [ i 1 i 2 i 3 ] T is the vector of the three-phase currents (measured by CTs) in the primary windings of the transformer. Therefore, the proposed protection framework requires measurement of the three-phase current and voltage waveforms at the primary side (energisation side) of the transformer. In the generator substations, the three-phase current and voltage measurements are always available at the primary side of the transformer (terminals of the generator), which means that the proposed method can be straightforwardly used for generation substations. Likewise, in transmission substations, the distance relays are often used to protect the transmission lines. In some countries, the distance relays are also used as back up relays for the protection of power transformers. The distance relay uses three-phase voltage and current measurements to calculate impedances. These current and voltage signals can be reused for the differential relay of the substation transformer, as well. Hence, in many situations such as generator substations and transmission substations with distance protection, the proposed method can be applied with no additional installations.

PRINCIPLES OF THE PROPOSED HYBRID TECHNIQUE
The principles of the proposed hybrid method for detecting and classifying the transformer internal faults are described in the following subsections.

Discriminating the internal faults from inrush current
Initially, three residual signals are defined as the differences between the measured and estimated high-voltage (HV) currents. The residual signals are normalised to the rated current of the HV windings. The values of the percentage normalised residual signals (PNRSs) are compared with a predetermined threshold value to identify the energisation conditions. The PNRSs are defined as follows: where i Rated is the rated current of the primary windings andî j is the estimated value (by EKF) of i j . If the PNRSs are lower than the selected threshold value, the transformer is healthy and the protection is bypassed. If one or more PNRSs exceed the threshold value, the transformer is faulty. Hence, the protection relay disconnects the transformer using the main circuit breaker. Therefore, the protection relay operates when the following condition is satisfied for at least one j: To select the optimal threshold value, the Otsu thresholding method is used [25]. The main idea of the Otsu algorithm is to minimise the overlap between two classes of data, namely inrush current and internal fault classes. Based on this algorithm, a threshold value of 10% is selected. Promising results regarding the algorithm accuracy and operating time are obtained with the selected threshold value.

Classification of the internal faults
In (4), it is observed that each system state has a meaningful relationship with the parameters of the corresponding winding and core characteristics. For instance, 1 is linked with the parameters r 1 , Γ 11 , and input voltage v 1 . Moreover, it is expected that an internal fault within a certain winding of the transformer mostly influences the relevant parameters in that winding. Therefore, the state estimation error of the faulty winding will be increased. Accordingly, the entries of the covariance matrix of the estimation error can be interpreted to obtain valuable information about the internal faults. In this paper, the SVM classifiers are used to achieve the foregoing objective. The normalised averaged diagonal entries (NADEs) of the covariance matrix of the estimation error are first defined as follows: where p kk denotes the value of p kk (diagonal entries of the covariance matrix of the estimation error) in the j th sample. Also, N j is the sample number at which the EKF estimation is terminated (the internal fault is distinguished). Then, the inputs of the SVMs can be expressed as follows: The spider plots of Figure 5 indicate the variations of P ii (i = 1,2,..,6) for different fault types, locations, and severities. Normally, it is difficult to interpret the relationship between the specifications of an internal fault and P ii . However, SVM classifiers can fulfil this objective, perfectly. The SVM finds an optimal separating hyperplane by maximising the margin between the separating hyperplane and the data [26,27]. Three multi-class classification problems, namely SVM 1 , SVM 2 , and SVM 3 are defined and solved to determine the type, location, and severity of the internal faults, respectively. Each classification problem is handled using a one-versus-rest (1VR) strategy [26]. SVM 1 determines the fault type among five different classes, namely line-to-ground, line-to-line, line-to-line-toground, primary-to-secondary, and turn-to-turn internal faults. SVM 2 identifies the location of the internal fault among twelve different categories, namely XP, XS, XY, XPS, where XP and XS indicate an internal fault in primary and secondary windings of Phase X, respectively. Likewise, XY denotes a line-to-line fault between Phases X and Y, and XPS denotes a primary-tosecondary fault in Phase X. Likewise, the fault severity is classified into three different levels using SVM 3 (see also Table 1) The inputs of the SVM classifiers are defined in (16). The radial basis function (RBF) kernels are adopted for the nonlinear multi-class classification problem. The RBF kernel is chosen for its popularity, perfect nonlinear mapping capability, and having a fewer number of free parameters [26]. With the RBF kernel, the decision function f (x) related to the SVMs can be represented as follows: where l is the number of support vectors, i is the Lagrangian multiplier, y i is the target value, x i is the support vector, is the RBF parameter, and b is the bias parameter. The SVMs are trained using the Classification Learner application from MATLAB software. The parameters of the SVM classifiers, that is, penalty factor and gamma, are refined using the five-fold cross-validation technique. A total of 1677 test cases which consider different fault levels, locations, and types are recorded using PSCAD/EMTDC software for training and testing of the SVMs. The collected dataset is separated into two different datasets with 1227 and 450 test cases for training and testing of the SVM classifiers, respectively. The datasets are recorded using the Multiple Run Feature of the PSCAD software, which is a great option to repeatedly run a simulation and automatically change the value of influencing parameters from one simulation run to the next. The trained decision functions of the SVMs are then used to solve the three classification problems.

Eliminating the effect of current transformer saturation
To eliminate the effect of current transformer (CT) saturation, an effective CT saturation detection and compensation algo-rithm is used in this paper [28]. Different stages of the utilised algorithm for detection and compensation of the CT saturation are briefly described in the following: 1. Detection of the CT saturation inception.

Extension of the unsaturated interval of the fault (inrush)
current for three additional samples after the beginning of the CT saturation. For this purpose, the leastsquares digital filter is used to fit the data points of the unsaturated region to a model of the fault (inrush) current. 3. Estimation of the CT magnetic flux using the difference between the real secondary current (measured) and the extended samples in part 2 to obtain the CT magnetising current.

RESULTS AND DISCUSSIONS
A total of 2800 test cases are obtained using simulations in PSCAD/EMTDC software. These test cases consider different factors including the switching angle, remnant flux, fault severity, fault type etc. In some test cases, the effect of the CT saturation is also considered utilising another CT with a lower rating. The methods followed to simulate turn-to-earth and turnto-turn internal faults are based on dividing the faulty winding into two and three different windings, respectively [29]. To achieve the foregoing objective, a transformer model based on mutually coupled resistive and inductive branches is obtained using BCTRAN routine of the EMTP simulation software. The obtained model is then modified and used for simulation of the internal faults using a simple approach explained in [29]. Figure 7 shows the simulated system for obtaining the data sets. In Figure 7, V s , X s , and R s are the source voltage, source reactance, and source resistance, respectively. In addition to the simulation data, a small experimental dataset which includes  Table 2. The specifications of the simulated transformer and system parameters can also be found in the APPENDIX. The EKF is used for estimation of the currents at the primary side of the transformer and some typical results are shown in Figures 8-10. In all the test cases, the transformer is energised at t = 0.3s. In Figure 8(a), the three-phase measured and estimated currents for a typical inrush event are depicted. As seen, the HV currents are precisely estimated when the transformer is healthy. The PNRSs for this sample test case are shown in Figure 8(b). It is evident that all the PNRSs are far lower than the selected threshold level and thus, the energisation condition is successfully discriminated. However, in the case of energising a faulty transformer, the EKF-based estimator fails to precisely estimate the HV currents. As seen in Figure 9, when the transformer with 10% turn to turn internal fault in Phase C is energised, the PNRSs are increased and exceed the threshold value of 10%. From Figure 9(b), it can be observed that the operating time of the proposed algorithm is less than a quarter of a cycle of power system frequency. Figure 10 shows the performance of the proposed method under CT saturation condition. For simplicity, only waveforms relevant to Phase A are depicted. Upon the CT saturation occurs, the saturated waveform is reconstructed using a simple algorithm [28]. When a healthy transformer is energised, the PNRS is lower than the threshold value. However, when a transformer with an internal fault is energised, the PNRS of Phase A exceeds the threshold value in about 4 ms. Investigating the whole 2800 simulation test cases reveals that the proposed EKF-based method gives an accuracy of 99.79% and hence, provides suitable dependability for discriminating the inrush current from internal faults. The accuracy is measured using the following index: where N is the total number of test cases and N Succeeded is the number of test cases that have been successfully distinguished by the proposed method. In all the test cases, the decision making is accomplished within less than a quarter of a cycle, which means the operating time of the relay is also promising using the proposed method. The mean operating time of the proposed method is calculated 4.6 ms using the following formula: where Δt j is the operating time for the jth test case. Three multiclass classification problems for determining the type, location, and severity of the internal faults are defined using SVM classifiers. Different categories for the fault type, location, and severity are described in Table 2. Moreover, the classification of the fault severity is described in Table 1. Implementation of the SVMs is fulfilled in MATLAB using the Classification Learner environment. A total of 1227 test cases are used for training the SVM classifiers. The performances of the trained SVMs are then evaluated using a separate dataset containing 450 test cases. The results are summarised in Table 3. The overall accuracies of 98.65%, 98.19%, and 98.42% are obtained for SVM 1 , SVM 2 , and SVM 3 , respectively. It is evident that in addition to  Turn to turn internal faults with less than 1% of the winding influenced+ LG, LL, LLG, and primary to secondary faults applied at less than 5% of the winding from the transformer terminal LOW Turn to turn internal faults with 1-3% of the winding influenced+ LG, LL, LLG, and primary to secondary faults applied between 5% and 30% of the winding from the transformer terminal MEDIUM Turn to turn internal faults with more than 3% of the winding influenced+ LG, LL, LLG, and primary to secondary faults applied at more than 30% of the winding from the transformer terminal HIGH the fast and dependable performance of the proposed algorithm for discriminating the internal faults from inrush current, this method classifies the fault type, location, and severity with high precision.

COMPARISON WITH OTHER ALGORITHMS
In this section, a comparison between the proposed method and some of the existing algorithms is fulfilled. A random dataset of 100 test cases is selected from the main dataset. Based on these test cases, the threshold values of the algorithms are obtained using the Otsu method. The algorithms are evaluated using the same test cases to have a fair comparison. The comparison is fulfilled in terms of the algorithm speed, accuracy, robustness against measurement noises, complexity, and ability to classify the internal faults. The results are listed in Table 4. It can be concluded that the proposed method has superiority over the conventional algorithms in most of the terms, though these advantages are achieved at the cost of higher computational complexity.

ROBUSTNESS AGAINST MEASUREMENT NOISES AND CTS AND PTS ERRORS
The measurement noises and CT and potential transformer (PT) measurement errors are inevitable in the real field. To evaluate the effects of such errors on the performance of the proposed method, a robustness analysis is fulfilled. To consider the effect of the measurement noise, band-limited random white noises with different signal-to-noise ratios (SNRs) are added to the three-phase measured current and voltage waveforms. The settings of the EKF-based estimator are kept unchanged. Sample results for the estimation of the noisy current waveform in Phase A related to a typical energisation test case are shown in Figure 11. It can be deduced that even when the noise level is considerably large, the PNRS is far lower than the selected threshold value of 10%. Similar results are obtained for estimated currents in Phases B and C. To evaluate the robustness of the proposed method against the measurement noises more effectively, a random dataset of 100 test cases are selected from the main dataset. The random white noises with SNRs of up to 30 dB are added to the voltage and current measurements. By evaluating the whole 100 test cases, an overall accuracy of 99% is obtained, which demonstrates the high immunity of the proposed approach against the measurement noises.
To test the robustness of the fault detection method against the CTs and PTs errors, a total of 100 random test cases which include different energisation and internal fault conditions are selected from the main dataset. The three-phase current and voltage measurements are changed to account for ratio and phase shift errors of the CTs and PTs, respectively. All the settings of the EKF algorithm are kept unchanged. The overall accuracy of the fault detection algorithm (with EKF) is calculated for different levels of CTs and PTs errors and the results are shown in Figure 12(a). It can be concluded that even if the errors of the CTs and PTs are too high, the energisation conditions are discriminated from the internal fault conditions with reasonable precision. To test the robustness of the fault classification algorithm against the measurement errors, another    Figure 12(b). On the whole, it can be deduced that the overall accuracy of the proposed hybrid approach will be higher than 97% in the worst case when it is considered that the CTs and PTs have error values as large as 4%. In practice, however, the maximum error of most instrument transformers used for protection purposes does not exceed a few percents. Thus, the proposed approach effectively deals with the problem of CTs and PTs errors, as well.

FEASIBILITY OF EMBEDDED IMPLEMENTATION
The hardware requirements of the proposed method are estimated to assess the feasibility of embedded implementation with digital protective relays. To fulfil the foregoing task, hardware-in-the-loop (HIL) simulation is performed. The three-phase voltage and current measurements from PSCAD software are first converted to analog signals through a digitalto-analog (DAC) converter. These analog signals are then passed through simple resistive-capacitive (RC) low-pass filters (LPFs), which simulate the low-pass antialiasing filtering process of the real digital relays. The obtained analog signals are discretised with a sampling frequency of 5 kHz once again. The obtained digital waveforms are then processed using the EKF algorithm, which is completely realised using the dSPACE1104 control board. The needed memory for the proposed approach is about 12 kB. The execution time of the algorithm is about 10,500 cycles of the master digital signal processor (DSP) of the dSPACE1104 card, which equals to an operating time of about 70 µs. The obtained operating time is lower than the duty cycle of the master DSP based on a 5 kHz switching frequency, that is, only 21% of the DSP resource is used for running the proposed EKF algorithm. Therefore, there is no bottleneck for embedded implementation of the proposed EKF-based method. The computational effort of the proposed EKF-based approach can further be reduced by optimising the programming codes and using efficient compilers.

CONCLUSION AND FUTURE WORK
In this paper, an effective method for detection and classification of power transformer internal faults is presented. The proposed method detects and classifies the internal faults using a hybrid approach. The fault detecting is fulfilled comparing the PNRSs with a predetermined threshold. Additionally, the fault classification part of the proposed method determines the fault type, location, and severity using SVMs. The performance of the proposed method in fault detection and classification is confirmed using a large number of simulation test cases, which are obtained in various conditions using the PSCAD/EMTDC software. The results confirm that the operating time of the proposed fault detection method is lower than a quarter of the cycle of power system frequency. The overall accuracy of about 99% and 98% is obtained for detecting and classifying the internal faults. It is shown that the proposed method has suitable immunity against the measurement noises and errors. Some HIL tests are also conducted to confirm its real-time feasibility.
For future work, it would also be worthwhile to assess the performance of the proposed hybrid approach in detection of multiple faults occurring in the power transformer.