S-Transform based fault detection algorithm for enhancing distance protection performance

This paper presents a new fault detection algorithm based on the Fast Discrete Stockwell Transform. The algo- rithm can improve the functionality of existing distance protection and resolve shortcomings identified during the fault detection process in case of fault occurrence in systems with a high penetration of power electronics- based generators. Reported results of the operation of commercial distance relays of four different vendors show that all relays experience difficulties during ungrounded faults. An RTDS testbed is developed for extensive hardware in the loop testing, comprising electromagnetic transient models of Type-3 and Type-4 wind generators. The proposed algorithm successfully overcomes the identified problems for cases where commercial relays maloperate. The threshold parameters for the fault detection are set by using the energy content attributed to the Fast Discrete Stockwell Transform time – frequency domain signal. Other distance protection modules such as the determination of directionality, phase selection and the computation of the impedance, which are necessary for the protection selectivity, are developed based on currently available solutions applied in commercial relays. The new algorithm has been extensively tested using RTDS for various fault conditions and the obtained results are reported in the paper.


Introduction
In the future, a large number of coal and nuclear power plants will be decommissioned and replaced by renewable energy sources (RES). It is known that low system inertia reduces the fault current level. Furthermore, bi-directional fault currents and the different composition of energy production and load make the operation, control and protection of such power systems increasingly complicated. Distance protection will also suffer from this as it is essentially designed for grids dominated mainly by the classical synchronous generation.
RESs such as photovoltaic systems, as well as Type-3 and Type-4 wind turbines (WTs) are commonly interconnected to the grid through partial or full-scale power electronics (PE) voltage source converters (VSC). This gives RES the capability of fast switching when necessary, in order to control their output power upon disturbances. The VSC control strategy and its electric configuration are key factors for the determination of the fault current contribution of a PE-based generator [1,2]. This is significantly different from the synchronous generators' (SGs) contribution in terms of fault duration, magnitude and shape of voltage and current waveforms.
The essentials of testing protection behaviour for power systems dominated by converter-based renewable generators (RWGs) are addressed in [3] and [4]. It is shown that the use of Type-3 and Type-4 WTs considerably increases nuisance tripping of the protection system. In [5], an intelligent algorithm for the protection of modern power systems is presented. The developed algorithm is suitable for distribution voltage levels. In [6] and [7], distance protection performance is evaluated for different faults on a line supplied by a RWG from one side and by synchronous generation from the other side. Some cases in which protection may fail to clear the fault are highlighted and discussed. In [8], an adaptive model based on provided local information and artificial neural networks is used to improve the distance protection and prevent the wrong relay operation during changes in wind farm conditions. The loss of protection sensitivity and selectivity near DFIG is also observed in [9] where an adaptive distance protection scheme is proposed. In [10], the authors proposed a combination of distance and differential protection to overcome distance protection issues.
In this paper, a significant number of viable scenarios are studied to determine cases where commercial distance relays experience fault detection difficulties. It is observed that in more than 50% of the cases, relays fail to detect phase-to-phase (LL) and three-phase (LLL) faults. Delayed tripping is also reported during phase-to-phase-to-ground (LLG) faults. Therefore, in this paper a new approach with an improved fault detection algorithm is proposed that can detect and identify all fault types. A realistic test system is developed for testing purposes, which includes detailed electromagnetic transient (EMT) models for Type-3 and Type-4 WTs. The test system is suitable to produce different penetration levels of renewable generation in combination with conventional generation. The RES control strategy can be easily switched from positive to positive and negative control strategies, facilitating the testing with a broad range of possible scenarios that may occur in practice.
The rest of the paper is organized as follows. Section 2 presents the results of the hardware-in-the-loop (HiL) testing of four commercial distance relays. Section 3 explains the implementation of the Fast Discrete Stockwell Transform (FDST) as a new solution for reliable fault detection. In Section 4, the performance of the FDST algorithm is evaluated. For this purpose, the FDST algorithm is incorporated together with other distance relay modules which are important to guaranteeing protection selectivity. Thereafter, also in Section 4, a comparative analysis is performed. Finally, the conclusions are summarized in Section 5.

HiL testing of commercial distance relays during faults supplied by RWG
Due to the complexity of RWG dominated transmission systems, distance protection of these systems is quite challenging. This paper uses detailed EMT models of Type-3 and Type-4 WTs equipped with appropriate positive and negative sequence control. The modelling is performed according to the approach explained in [11][12][13][14][15][16]; however, the modelling details were beyond the scope of this paper. Fig. 1 shows the single-line diagram of the test system developed. The Type-3 and Type-4 WT aggregated models are connected to the grid by lines 4-6 and 5-7, respectively. The distance relays are located at buses 6 and 7. Four different manufacturers' relays, denoted by A, B, C, and D, are simultaneously tested in a HiL real-time testbed.
A C-code script is developed and used in order that a large number of tests can be run automatically. The same relay settings are applied to all relays. In order to investigate the validity of the setting values, the generators in Fig. 1 are replaced by SGs. The results of extensive number of HiL tests perfromed validate the accuracy of the relay operation for 99.97% of all simulated cases. The cases comprise all types of faults, different fault locations (including Zones 1, 2 and 3 reverse) and infeed power from the grid. Furthermore, the relay's performance is recorded and presented statistically for nearly 5000 different cases when PE-based generators are connected to the system. Fig. 2 shows the results related to missed trips for different fault types. The performance of each relay differs as each manufacturer has its own algorithm used to detect and pick up a fault. For all investigated relays, the majority of reported undetected faults are related to LL and LLL faults. In conclusion, the impact of the RWG is such that faults remain undetected in some cases, thus this paper proposes a solution for this problem. Fig. 3 shows the proposed enhanced relay (ER), which comprises different protection modules by applying a novel fault detection algorithm based on the FDST [17,19]. The FDST energy (S-energy) per phase is used to produce a selective trip signal by the ER. The S-energy changes abruptly after the fault inception even when a PE is connected to the protected line, making it an efficient indicator for fault detection. The other three modules are the phase selection [20,21], the directionality [22] and the impedance modules [23,24]. It should be noted that the developed models which represent these functions are needed in order to design a full relay model. The change of the fault impedance due to the fluctuating power resulting from the variable wind speed is challenging even for the ER. This case was beyond the scope of this paper, and it will be a topic for future research.

The FDST computation
The FDST of a discrete-time series x[kT], k = 1,…,N corresponding to a current signal x(t) where the sampling interval is T, can be expressed as: where n = 1,…,N-1 and m = 1,…,(N/2) represent the time and the frequency point indices of a given cycle, also known as a particular voice [18]. The matrix H is a concatenate and rotated matrix resulting from the FFT applied to x [kT]. The matrix is divided into N localized vectors as: where M is equal to N/2. The W term in (1) is a modified twodimensional (2D) Gaussian window [19] acquiring localization in frequency and time domain. The 2D Gaussian window is defined as:  where In the above equation, F is the window factor, b is the scaling factor that controls the number of oscillations in the window, and a and c are positive constants. The value of parameter c, which varies between 0 and 1, contributes toward the capture of damped hidden frequencies. The N shifted Gaussian window is used as a filter to decrease the computational burden of the discrete S-transform by filtering out unwanted frequency information [25]. In order to acquire the windowed time-frequency information, the concatenated matrix H is multiplied by the 2D Gaussian window matrix element so that G is a Hadamard product of H and W matrices: The FDST defined in (1) can be seen as the Inverse Fast Fourier Transform (IFFT) of each time step for every row of G: The resulting S-matrix contains the instantaneous phasor values for each frequency [25]. In this work, an FFT radix 4 is applied. The total number of mathematical operations to compute the FDST is 2Nlog 4 (N) additions and N(2 +log 4 (N)) multiplications [19] (for instance, for N = 16, 64 additions and 64 multiplications).

Energy from the fast discrete stockwell transform
The S-energy can be computed by the FDST matrix from During steady-state conditions, the S-energy is nearly a flat signal; however, it abruptly increases during transient conditions. The S-energy sensitivity is a suitable indicator for different disturbances, including faults.
The window length and the number of samples should be appropriately set to determine the S-energy. In this work, the window is adjusted to F = 1.5, a = 3, b = 0.2 and c = 1.4. The number of cycles per window and the sampling frequency should be selected in a way to obtain a useful signal with a minimum number of samples. In this way, sufficient data will be provided to limit rounding errors and avoid signal oscillations. Three window sizes with a sampling rate of 250 Hz as seen in Fig. 4(a) and three sampling rates with a window size of 16 points as seen in Fig. 4(b) are analyzed in order to choose the suitable window size and sampling rate. Fig. 5(a) shows the S-energy behaviour for the three different windows with a sampling rate of 250 Hz. As it is expected, the S-energy amplitude is directly proportional to the amount of data in the window. Nevertheless, there is no considerable difference between the S-energies within the first milliseconds after the fault inception. For all window sizes, the S-energy starts increasing 0.04 ms after the fault   occurrence. In Fig. 5 (b), three different sampling rates for a three-cycle data window are shown. The S-energy accuracy depends on the sampling rate and increases for higher sampling frequencies. However, fast sampling rates result in an increase in the computational burden. It can be seen that after 0.008 ms (1/4 cycle), the computed 250 Hz S-energy amplitude changes approximately by 2 dB and 5 dB in comparison with the S-energy obtained for 500 Hz and 1000 Hz sampling rates. The analysis also shows that 1000 Hz sampling rate increases the speed of fault detection at the expense of more computational burden.
Extensive simulations indicate an optimal sampling rate of 250 Hz and an optimal data window of 3 cycles. The number of measurement points per window is rounded to 16 as a power of two is needed to apply FFT.
The S-transform generates a two-dimensional signal array, defined simultaneously for both time and frequency domains. The S-energy is the energy from the S-transform signal, and it does not only change its magnitude because of changes in time, but also because of its frequency domain. Thus, it increases abruptly and quickly during faults, and significantly more than the increase resulting from one-dimensional signals [26]. Consequently, the indicator becomes more sensitive during disturbances. The energy is computed by adding the square magnitudes in the time and frequency domains, just after the S-transform is computed, and without an extra algorithm. The S-transform advantages over existing techniques, such as Short Time Fourier Transform (STFT) and Wavelet Transform (WT), are listed in [27]. The data window adjustability, the ideal time and frequency resolution [28], makes the Senergy an excellent indicator for signals with small distortions.

Determining S-energy threshold
To employ the S-energy as a fault detection module, it is necessary to introduce a threshold as a setting parameter for this function. The Senergy amplitude depends on the data window length, the frequency spectrum, the sampling rate, and the signal itself in steady-state condition. Based on conducted analyses of the broad range of current/power spectrum, three effective approaches are proposed to determine the threshold value. These approaches are all based on the steady-state RMS value of the signal.

Interpolation
Linear interpolation is performed by using the RMS values of the load current Irms of the signal and its corresponding S-Energy values (SE) as shown in Table 1. This table contains the infeed power, the measured RMS current, the S-energy and the S-energy with a safety margin (STHsm). The S-energy threshold with a safety margin is expressed in dBs and it also contains the maximum total harmonic distortion allowed (THDSE) by the Standard STD 519 (5%) [29] and a safety margin (smS E ) of 20% the S-Energy.

Visual inspection
A visual inspection can be made to obtain the S-Energy threshold by following the blue line in Fig. 6 and the corresponding threshold shown in red. The signal threshold is set based on the S-energy and threshold obtained in Table 1.

Empirical calculation
By fitting the threshold signal with a safety margin (STHsm) as shown in Fig. 6, an empirical equation is obtained, which can be used as a threshold.
Equation (8) is represented by a rational approximation using a 4th order polynomial function and it uses the provided RMS current. The sensitivity analysis shows that a 4th order polynomial produces the most accurate approximation for the threshold with a safety margin as proposed in (7).
The S-transform is an exceptionally useful tool, which has been implemented in the past for solving different types of technical problems. Typical examples are presented in [30][31][32].

Performance evaluation
This section provides details about the performance of the new approach by comparing the results to those obtained from the HiL tests of commercial relays and the RTDS distance relay existing in the library model. The operation of distance protection is tested for faults within Zone 1 and Zone 2 as well as for backward faults within Zone 3.

Overall performance
The tested cases with the commercial relays and the ER for the system shown in Fig. 1    Generators are modeled in detail as single wind generators. Hence, they are scaled in its interconnection transformer to represent the entire wind farm without neglecting any dynamics from their electrical, electronics and control elements.
The generator power levels in the performed test cases are 40 MW and 200 MW. The scenario with more unexpected cases (no trips or delayed trips) is when the RWG infeed is 40 MW. Table 2 summarizes the trip times of relays (denoted by A, B, C, and D) as well as the RTDS library model and the ER for 20 cases with bolted faults on line 4-6, out of around 5000 tests performed. Faults are applied at 0 km, 50 km, 70 km, 90 km and also right behind the relay (0 km back). The model successfully trips in all cases with an average trip time of 33 ms for Zone 1, 430 ms for Zone 2, and Zone 3 is disabled in all the relays. The obtained results confirm that the S-energy module (which is incorporated in the new model) shows good performance with regard to the fault type detection at any fault distance. In general, based on the extensive number of tests conducted, relay C performs better compared to the other investigated relays. Thus, further results of the analysis of ER performance will be only compared to the performance of relay C. For all line-to-ground faults, the ER successfully eliminates the fault by tripping the correct phase. A detailed inspection, regarding time operation, for the modules used by relay C and ER is presented in Table 3.

Overall performance
The RSCAD multi-function distance relay [33] is tested along with the commercial relays and the ER. The relay is based on starting, phase selection, and Mho or Quadrilateral characteristic. The starting function is based on the current magnitude. For all the tested relays, the settings are chosen in a way to test the performance under the same circumstances. As shown in Table 2, results are very close to those of the commercial relays. The current threshold is not restricted to a minimum or maximum value. Commercial relays have margin values of 0.5 A or 0.25 A, depending on the vendor. In order to test the limits of the RSCAD relay, the current threshold is set to 0.1 A. Under this condition, the relay trips for LG, LLG, and LL faults. For LL and LLG faults with a power infeed of 40 MW, however, the current threshold must be set lower than 0.1 A. The RSCAD model does not have phase selection and all the functions are performed in parallel.

Bus 6 line-to-line fault
According to Section 2, most failures occur during ungrounded faults, namely LL and LLL faults when having low power infeed from the RWG. In this subsection, the ER is examined during an LL fault occurring on bus 6, which is one of the critical cases when protection fails as reported for all tested relays. Fig. 7 shows the recorded waveform of relay C and its performance compared to that of the ER. In this case, the Type-3 wind turbine at bus 6 produces 40 MW active power, and a bolted LL fault at 70% of line 4-6 with a 0 • inception angle is applied. Fig. 7(a) shows the relay secondary voltage of an ideal voltage transformer.
During the fault, voltages in the faulty phases drop in amplitude and overlap. Fig. 7(b) shows the recorded currents of relay C. As it can be seen, during the fault, the current amplitude exceeds twice the steadystate current amplitude. As it is highlighted, Relay C trips whilst the rest of the relays do not operate. After the fault inception, the PE control regulates the voltages and the currents according to the grid code and converter constraints [12,13,16]. The ER performance is checked for the same fault case, and the trip signal as well as the circuit breaker operation are highlighted in Fig. 7(c). Fig. 8 (a) and Fig. 8 (b) summarize the fault detection results for faults occurring at 70% and 90% of the line. After 9 ms from the fault inception, the computed S-energy of phases A and B exceeds the threshold and triggers the fault computation functions (see Fig. 8 (a)). The signal corresponding to the S-energy AB (once A and B S-energies are ON) picks up 10 ms after the starting and a trip command is generated around 30 ms after the fault inception. The same type of the fault is applied at 90% of the line (Zone 2).
The starting and the pickup time are similar to the previous case as the S-energy shows a similar behaviour (see Fig. 8 (b)). However, in this case the trip command is generated 0.4 s after the pickup as the fault lies within Zone 2. Fig. 8 (a) and (b) show that the S-energy is immune to the fault distance. The binary signals from the ER for an LL fault occurring at 70% of the line are shown in Fig. 9. At the bottom of Fig. 9, the fault inception (0.799 s) and its duration are shown in blue. Then, from the top to the bottom, the graph shows the instant when S-energy occurs (0.7992 s), which is 4 ms after the fault. Thereafter, the faulty phase identification detects an AB fault at 0.8128 s, and the fault direction is determined as a forward fault at 0.7997 s. Finally, the last module computes the fault impedance, which drops inside Zone 1 at 0.8119 s. As we observe an LL fault between phases A and B, the ER provides three phase trip signal at 0.8128 s, which is 32 ms after the fault occurrence. Besides, the reclosing time action is seen as a penultimate indicator, which occurs at 1.82 s after the fault is cleared. The impedance trajectory corresponding to the foregoing fault case is shown in Fig. 10. In both cases, for the faults at 70% and 90% of the line length, the impedance entered the protection zones. This is an expected behaviour that is also recorded in commercial relays.

Bus 7 three phase fault
In this subsection, the model is tested during an LLL fault at bus 7. It is noticed that the commercial relays do not pickup and/or succeed in   phase selection during the first five cycles. An LLL fault is simulated at 70% of the line length and behind the relay for a fault resistance of 1 Ω. The resulting S-energy for the forward and reverse fault is shown in Fig. 11. The fault in Fig. 11(a) is detected 6 ms after the fault inception. The trip command is provided at 0.807 sec and it is cleared at 0.857 s. It can be seen from Fig. 11(b), that the fault is   detected at 4 ms. However, the directional module prevents the generation of a trip command as the condition for the direction is not fulfilled. In this case, the fault current and the energy computed by the relay corresponds to the fault current provided by the grid (larger in magnitude than the energy resulting from the WT fault current). From these cases where the fault resistance is 1 Ω, it can be concluded that the Senergy is immune to the fault resistance and the inception angle. The angles that determine the fault direction are shown by unit vectors in Fig. 12. The angle is computed by using the fault current and the positive polarization voltage (VAB*, VBC*, VCA*) [22]. During steady-state condition, the direction of the vectors (blue arrows in Fig. 12 (a) and Fig. 12 (b)) remains in the load area denoted as nondirectional. For a forward fault, the direction vectors abruptly change and reach a position lower than 90 • (red arrows Fig. 12 (a)). For a reverse fault, the angle is greater than 120 • (red arrows Fig. 12 (a)). The fault direction is determined within milliseconds after the fault inception (black arrows in Fig. 12). The direction is obtained only for the AB sequence after the faulty phase is determined. Fig. 13 shows the impedance trajectories for both cases. For a fault at 70% of the line, the impedance drops inside Zone 1 and for a fault behind the relay, inside Zone 3. The trajectories of sequences AB, BC, and CA are drawn accordingly. In Fig. 13, the black line refers to the forward fault trajectory and red line refers to the reverse fault trajectory. In this work, the third zone relay is set to 0.6 sec. Table 3 compares the performance of the ER and relay C in its available functions. It can be seen that the ER equipped with S-transform is approximately 7.5 times faster than relay C, on average, in terms of pickup time for LG faults. For LL faults, the pick-up time of the relay C is longer than 80 ms, which results in delayed trips for faults inside Zone 1. For most faults occurring in Zone 1, the ER operates faster than relay C. It can also be seen that the model overreaches less than relay C for high impedance faults. In general, the ER outperforms the relay C because the fault identification is much faster and the sensitivity during high   impedance faults is higher.
Regarding the influence of the fault resistance, according to the results shown in Table 3, the performance of the ER model and relay C is basically the same for LG faults. Both relays operate with full success for LG faults. The ER shows better performance in LLG and LLLG. The directionality function in relay C is based on the impedance trajectory. The direction and the impedance are computed at the same time.

Capacitor switching and power swings
The ER is also tested during capacitor bank (CB) switching and power swings. Relays are expected not to operate during these phenomena.
The capacitor or load switching produces a fast transient that can be erroneously identified as a fault inside a protection zone and trip the line. As a study case to test the ER performance, a CB is connected to bus 6 and switched on and off within one second. The capacitor bank has a star grounded connection with a capacitance of 1.2835 μF. The operation times for each module are shown in Table 4. As expected, the fault detection, phase selection, and directionality are active, whilst the impedance determination is blocked. In this case, the relay does not trip because the impedance trajectory does not enter the protection zones.
Large disturbances in the power system may cause power swings. The impedance trajectory can be seen in and out the protection zones. A 500 MW SG is added at bus 7 to produce a power swing in the grid. The relay can generate wrong or undesired tripping once the impedance trajectory enters the protection zones. However, the S-energy used in the ER does not rise enough to pick up the model and disables the trip command. The ER modules times for these study cases are reported in Table 4.

Conclusion
In this paper, an advanced algorithm for fault detection, based on Fast Discrete Stockwell Transform (FDST), is presented. It can be used to enhance the performance of distance protection. In order to investigate the performance of the proposed algorithm, it was integrated into a full model of a distance relay developed in RTDS environment. This significantly improves the sensitivity and the speed of fault detection and overcomes the difficulties experienced by present commercial relays during the detection of ungrounded faults. In addition, the low computation burden of the FDST makes this algorithm suitable to be implemented in commercial distance relays. Overall, the S-energy module is sensitive to different fault types, fault inception angles and fault resistances. The FDST is an essential module in the proposed ER. It is implemented as a settingless protection only with the RMS current value measured from the bus and it discriminates faulty and non-faulty conditions. However, in order to accurately identify any fault on the protected line, the FDST should be combined with other modules such as phase selection, directionality and impedance computation. The directionality and the phase selection module used in the ER are sufficient to determine faults in the vicinity of PE-based generators. A thorough analysis has been conducted to investigate the performance under different fault conditions in this work. The obtained results show that the ER outperforms four different commercial relays. The model presents accurate results in 99.88% of all conducted cases.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.