# Transmission Line Distance Protection Algorithm Based on Bayesian Inference

Authors: Rodrigo Rozenblit Tiferes; Giovanni Manassero Junior
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

This paper presents a distance protection method based on estimating the resistances and inductances seen by the relay using a least-squares (LS) approach combined with the application of Bayesian inference (BI) to determine fault probabilities in the protected zone over time. The use of BI in the presented algorithm increases the security and reliability of the LS-based line protection approach, which, although fast, may yield impedance estimates with numerical inaccuracies and therefore configures a solution capable of detecting faults internal to the protected zone rapidly and dependably. The authors tested the proposed algorithm using data sets obtained from transmission line fault simulations performed in Alternative Transients Program (ATP), which considered noisy measurements, current transformer (CT) saturation, capacitive-coupling voltage transformer (CCVT) transients, switching onto fault, close-in faults, frequency variations, uncertainties in the line’s parameters, and power swings. Furthermore, the authors tested the presented method against measurements recorded by field protection equipment during seventy-nine actual fault events on four existing 500 [kV] lines. The results indicate that the proposed algorithm is a fast, secure, and reliable distance protection solution and that using BI to the LS-based distance protection approach increases security while maintaining fast fault detection.


---

## Body

## 1. Introduction

Distance protection is a function widely used to protect power transmission lines, especially those with extra-high voltage (EHV) and extended lengths. Algorithms based on this function generally aim to determine the apparent impedance seen by the relay as the line impedance is proportional to its length. The modern microprocessor-based digital intelligent electronic devices (IEDs) combined with communication and signal processing technologies enable the development of various distance function solutions that seek to improve the protection’s speed and reliability [1].The authors of [2], for instance, propose in their research an algorithm based on solving an optimization problem that constructs an optimal zone for the protection of series-compensated lines. Reference [3], on the other hand, proposes an adaptive quadrilateral impedance zone whose boundaries are modified according to the results of a multioptimization problem solved through a genetic algorithm. In general, optimization-based protection algorithms constitute a trendy research line in the contemporary literature [4], and algorithms such as those proposed in the cited references may be advantageous as they allow the construction of adaptive protection zones rather than using predetermined fixed zones such as Mho or quadrilateral characteristics. Furthermore, adaptive zones can also be more effective at detecting high impedance faults than conventional zones. However, these methods generally require considerable computational effort to reach convergence and therefore need robust processing systems for real-time applications.Reference [5], in turn, presents an algorithm capable of compensating the zero-sequence couplings in double-circuit and in single-circuit lines, correcting the estimated impedance through a recursive approach. The main advantage of this method is that it compensates the zero-sequence magnetic couplings between circuits, unlike conventional double-circuit protection schemes. However, this solution is only applicable for double-circuit lines, and it demands a recursive approach with several operations, which requires a powerful real-time processing system.The algorithm proposed in [6] calculates the apparent impedance and fault distance by decomposing voltage and current phasors into sequence components. This method uses the sequence components theory to determine the fault’s distance in real-time and yields accurate results in most fault types and cases. Still, it may be affected by faults close to the relay bus and the fault distance’s estimate is adversely affected by the fault’s resistance. The authors of [7], in turn, present an effective solution to protect multiterminal lines that is not affected infeed currents produced from tapped lines or high fault resistances. Still, it demands synchronized measurements from all line terminals, which increases considerably the protection infrastructure’s complexity compared to solutions that use data from only one terminal. Reference [8], on the other hand, presents an algorithm that has an adaptive complex-plane characteristic that is not affected by the fault resistance, but like the other adaptive zone methods, it also needs robust real-time processing.The methods presented in [2–8] are examples of distance protection solutions that work with voltage and current phasors, estimated from the field measurements performed by the protection apparatus. Phasor-based protection solutions are the most common in commercial devices because of the low computational burden and effectiveness. In this context, the most consolidated phasor estimation technique is the discrete Fourier transform (DFT). However, the DFT is sensitive to the exponentially decaying direct current (dc) component present in the fault current signals and to variations in the system frequency, which may influence trip decisions and delay protection operation times. Besides, algorithms that use a full-cycle DFT (which, although slower, is more precise than the half-cycle DFT since the half-cycle DFT does not reject even harmonics) generally require a complete cycle of samples to estimate the apparent impedance accurately. Thus, it is common to find solutions that employ techniques other than Fourier analyses to estimate phasors.In [9], the authors present tests of a distance protection algorithm that estimates phasors through an LS formulation that yields the dc and harmonics components, initially proposed in [10]. The results indicate that the algorithm triggers the trip outputs in at least half a cycle, being faster than the methods based on the DFT, which usually present trip times close to one fundamental period. However, this algorithm may be affected by faults with fast decaying dc components and also by variations in the electric grid’s frequency. In [11], the authors use an adaptive Kalman filter to estimate phasors, obtaining accurate and precise results. Still, Kalman filters generally need high computational effort. In [12], the phasor calculation is based on phaselets, which are similar to the DFT but with fewer samples in the windows. Although phaselets estimate the faulted voltage and currents faster, they may be adversely affected by the dc and harmonic components. In [13], the estimation of phasors is performed by discrete wavelet transforms, which is also an approach that can estimate fault phasors rapidly and accurately. Still, the most appropriate wavelet transform may impact the phasor estimates, and this approach is also affected by the dc component and by frequency variations. In [14], in turn, the authors use Clarke’s transform to reduce the errors produced by frequency variations in the DFT, but using a Clarke-based routine to track the system’s exact frequency further delays the protection response and increases processing load. In addition, the trip times of all phasor-based protection methods are generally longer than half a cycle. Thus, to minimize trip times while maintaining reliability, researchers constantly propose non-phasor-based solutions.One of the main approaches used in current research that do not rely on phasor estimation is protection algorithms based on travelling waves (TW) theory. Reference [15], for example, presents a TW-based method that uses the polarities of the voltage and current TWs and arrival times to determine the fault direction. Generally, TW methods present considerably fast fault detection times. However, such solutions usually demand very high sampling frequencies ([MHz] magnitude) and robust processing systems. Besides, the wave-front detection may be affected by measurements with higher noise content.Another relevant approach is the knowledge-based protection solutions [4], which may use for instance artificial neural networks to detect and classify faults [16] or radial basis function neural networks [17] to prevent the relay to malfunction in situations concerning remote current infeed. This type of solutions may be effective in detecting all kinds of faults as long as they are properly trained. Still, knowledge-based algorithms require network training procedures concerning exclusively the line to be protected, which may lower these methods’ attractiveness compared to solutions that only need the real-time collected measurements.Finally, it is possible to mention the distance protection algorithms based on the line’s differential equation and model fitting, which are among the first proposed solutions that dispense voltages and currents phasors. Reference [18], for instance, presents an impedance estimation method that solves the line differential equation through interpolating integrals. This algorithm has the advantage of dismissing numerical derivative calculations. Nevertheless, solutions based on numerical integration favor the steady-state over high frequencies and, therefore, may be biased with error accumulation [19]. Thus, many authors have proposed solutions that calculate the current derivatives numerically and use such results in LS formulations to estimate the protected line’s impedance.In [20], the authors present an algorithm that calculates the resistance and inductance seen by the relay from the solution of an LS problem applied to the voltage and current samples collected in the relay bus. This method considers that the line can be modeled as a resistance in series with an inductance (lumped parameters) and that the current derivatives can be calculated from dit/dt=ik−ik−1/Ts, where Ts is the relay’s sampling period. This algorithm is considerably faster than those based on phasors, in addition to using fewer samples than one cycle. Furthermore, it is not influenced either by the dc component or by frequency variations and presents a similar computational burden compared to DFT-based solutions [19].In [21], the authors present an algorithm based on the LS technique applied to the line’s discrete-time model to estimate resistance and reactance values. However, the method is sensitive to frequency variations as it considers the reactance at the fundamental frequency. Reference [22], in turn, proposes an algorithm that detects faults by comparing the most recent samples with the corresponding samples from a cycle earlier and is based on a recursive LS solution to determine the fault distance and resistance. In [23], the fault distance and resistance are estimated through mathematical morphology and LS curve fitting. In [24], the LS-based proposed algorithm estimates the fault distance and arc voltage amplitude, assumed to have a square-shaped waveform.Distance protection methods based on estimating resistances and inductances/reactances using LS formulations [20–24] are advantageous over phasor-based solutions because they are considerably faster since they estimate impedance values within the zone more quickly. However, these algorithms may yield results that present considerable numerical oscillations and inaccuracies, especially in cases featuring high-noise content and high-frequency transient components in field measurements. In this way, the impedances estimated through phasor calculations are considerably more accurate than those estimated by solutions involving LS formulations, even if slower. Besides, LS-based line protection methods may be misled to protection failures (incorrect or nontripping) due to the numerical inaccuracies linked to the estimation of apparent impedances, which can compromise the security of solutions based on this approach.Therefore, the distance protection method presented in this paper combines the numerical estimation of the resistances and inductances seen by the relay through LS (based on [20] but adapted to consider mutual inductances) with the application of BI to determine the probability of an in-zone fault occurrence over time. The motivation for the development of this algorithm is to maintain the intrinsic speed of the estimation of apparent resistances and inductances by LS and increase the reliability of the trip decisions via the calculation of fault probabilities through BI, which are the basis to detect in-zone faults and trigger the trip outputs. The proposed algorithm was compared with the conventional DFT phasor-based solution since this is the most used approach by commercial distance relays and both methods are based on estimating impedance elements. Besides, it was compared with using only the LS impedance estimation, and the results show that the main advantage of using BI lies in cases featuring faults close to the zone boundary, where the LS-estimated resistances and inductances may fall in and out the zone repeatedly before settling. Furthermore, the authors compared the proposed solution with the current differential protection algorithm based on BI discussed in [25]. The main advantage of the distance function over the differential is that the distance relay does not need to exchange sampled values over a communication channel, as communication delays may slow the protection speed, especially on longer lines. However, the differential function generally does not struggle in cases of high-impedance faults. Therefore, when possible, the use of multifunctional relays can be highly beneficial to improve protection reliability.The authors used fault signals from ATP simulations, considering different lines and operating conditions, and data from seventy-nine actual faults recorded by IEDs of four different 500 [kV] lines from the Brazilian Interconnected Power System to test the proposed algorithm. The results indicate that it is fast and reliable, even in challenging situations featuring noisy voltage and current signals, CT saturation, CCVT transients, frequency variations, uncertainties in the line’s parameters, and power swings. Besides, the presented method does not require any improvement in the conventional distance protection infrastructure or high sampling rates.The upcoming sections of this paper are divided as follows: Section2 presents the distance protection topology used for implementing the proposed algorithm, the LS formulation used to estimate the resistances and inductances, the application of BI to calculate in-zone fault probabilities over time, the algorithm’s trip logic and settings, and the conventional algorithm based on DFT-calculated phasors that the authors developed for comparing with the proposed method. Section 3 presents the characteristics of the transmission lines and faults simulated in ATP and the algorithms’ responses. Section 4 shows the results of the sensitivity analyses, which considered conditions with noisy measurements, CT saturation, CCVT transients, close-in faults, line energizing onto faults, system frequency variations, uncertainties and errors in the line’s parameters, and power swings. Section 5 presents the test results of the proposed method against seventy-nine fault oscillographs recorded by the IEDs of four actual transmission lines from the Brazilian interconnected power system. Section 6, finally, concludes this article.


## 2. Methodology

### 2.1. Distance Protection Topology

Figure1 depicts a distance protection topology with IEDs in both terminals of the line: local L and remote R. In each IED, the respective phase voltages (measured by the voltage transformer (VT), which usually is a CCVT for high voltage (HV) and EHV lines) and currents (measured by the CT) are sampled at time intervals equal to Ts=1/fs, where fs is the sampling frequency.Figure 1
Two-terminal distance protection topology.Both IEDs execute the proposed algorithm in real time. The first IED that detects a fault sends a trip signal to its controlled circuit breaker (CB) and sends a transfer trip command to the other terminal through a dedicated communication channel. It is worth highlighting that the samples at each terminal do not need to be synchronized and the sampling frequency at each terminal does not necessarily need to be the same. Besides, a two-terminal protection infrastructure is not mandatory for the presented distance protection algorithm since it is possible to implement it at only one end of the line. However, a two-terminal solution increases reliability and fault detection speed.The authors set the primary zones (Zone 1⟶L and R) at 85 %of the lines’ positive-sequence resistance and inductance and considered a Mho characteristic [4] with the apparent inductance on the vertical axis and the apparent resistance on the horizontal axis for the application of the proposed method, in all cases (simulated and actual). Nevertheless, like phasor-based solutions, the presented method allows other zone configurations and characteristics.It is worth mentioning that, especially in one-terminal schemes and for high-impedance faults detection, secondary zones that overreach beyond the remote terminal and have intentional coordination delays can be used [28]. Still, as the evaluation of the proposed algorithm is based mainly on trip times, the authors considered in this research only the Mho characteristic’s first zone [4].


### 2.2. Estimation of Resistances and Inductances through LS

The proposed algorithm relies on six impedance elements: three line-to-ground (LG) elements (AG, BG, and CG) and three line-to-line (LL) elements (AB, BC, and CA). The LS formulations used to estimate the elements’ positive-sequence resistances and inductances (at every new sampling) are presented next. As mentioned earlier, the estimation performed by the presented method is based on [20]. Still, it is adapted to consider the mutual inductances between phases and the mutual inductance between circuits (in double-circuit lines).


#### 2.2.1. LG Elements

The algorithm considers that the voltage and current sample windows have a predetermined fixed sizeNLG to estimate the resistance and inductance of each LG element. At every new sampling, the most recent samples are incorporated into the respective windows and the oldest ones are discarded. The derivatives of the phase current signals at a generic discrete-time instant k, where k is a counter that grows from m=n−NLG+1 to n and indicates the position of each sample in the window (k=m⟶ oldest sample, k=n⟶ newest sample), can be calculated numerically as in(1)dit=k⋅Tsdt=Δik=ik−ik−dLGdLG⋅Ts,where dLG is a setting parameter (positive integer) and ik is the kth sample of the current data window. The current derivatives calculation in the presented method is more flexible compared to that presented in [20], in which dLG is always equal to one. Next, the formulations that involved the AG element are presented. The other LG elements are analogous.Single CircuitNeglecting the line capacitance and considering that it is perfectly transposed, the voltage in phase A to ground over timetvAt can be written as(2)vAt=RA⋅iAt+LA⋅diAtdt+LM⋅diBtdt+diCtdt,where RA and LA are, respectively, the resistance and inductance of phase A, LM is the mutual inductance between the circuit phases, and iA, iB, iC are the phase currents over time.Considering the discrete-time condition with the sampled values at a generic discrete-time instantk, (2) becomes(3)vAk=RA⋅iAk+LA⋅ΔiAk+LM⋅ΔiBk+ΔiCk.Since the voltage and current windows haveNLG positions for estimating the LG elements, the LS formulation for the AG element at every new sampling can be written from (3) as(4)IΔIAG⋅xAG=VA,where(5)IΔIAG=iAnΔiAnΔiBn+ΔiCn⋮⋮⋮iAkΔiAkΔiBk+ΔiCk⋮⋮⋮iAmΔiAmΔiBm+ΔiCm,(6)VA=vAn⋮vAk⋮vAm,xAG=RALALM.The LS solution of (4) yields vector xAG. Thus, the positive-sequence resistance of element AG RAG seen by the relay is equal to RA, while the inductance of element AG LAG seen by the relay is given by the difference between the inductance of phase ALA and the mutual inductance LM, as in(7)RAG=R,LAG=LA−LM.IfΔBk+ΔCk=−ΔAk for all instants k between m and n, the line is balanced and IΔIAG is singular since the third column is the opposite of the second. In this case, the third column of the matrix can be eliminated for the LS solution, which directly results in the AG element’s positive-sequence resistance and inductance (RAG and LAG). If the third column is not removed from the matrix, the estimated values of LA and LM respect the dependence relationship between the respective columns and have absolute values equal to the half of LAGLA=−LM=LAG/2. Thus, the verification of the mentioned singularity may be interesting from a computational point of view, since that if the second and third columns of IΔIAG are linearly dependent, the removal of the third column reduces the number of floating-point operations necessary to estimate the considered element.Double Circuit.In double-circuit lines, the mutual zero-sequence inductance between the circuitsLM0 has to be considered. Thus, (2) can be rewritten through(8)vAt=RA⋅iAt+LA⋅diAtdt+LM⋅diBtdt+diCtdt+3⋅LM0⋅di0ptdt,where(9)i0pt=iApt+iBpt+iCpt3,is the zero-sequence current and iAp,iBp,iCp are the phase currents of the parallel circuit over time. Therefore, for the inclusion of the zero-sequence magnetic coupling between circuits in discrete-time, (3) becomes(10)vAk=RA⋅iAk+LA⋅ΔiAk+LM⋅ΔiBk+ΔiCk+3⋅LM0⋅Δi0pk,where(11)Δi0pk=ΔiApk+ΔiBpk+ΔiCpk3,is the parallel circuit zero-sequence current derivative, calculated from the numerical derivatives of the respective phase currents. At every new sampling, the IED of each circuit (on the same bus) exchange the zero-sequence current derivative values so that the zero-sequence coupling is considered.For the LS formulation, the values of3⋅Δi0pk are added to matrix IΔIAG, forming a new column. Vector xAG now has four elements, the fourth being the estimated value of LM0. Vector VA, in turn, remains unchanged. After the LS solution, the calculation of the considered element’s resistance and inductance is the same as for the simple circuit Case (7).IfΔiAk+ΔiBk+ΔiCk=3⋅Δi0pk for all instants k between m and n, the derivatives of the zero-sequence currents in the two circuits are equal and IΔIAG is singular since the fourth column is the sum of the second with the third. If such singularity is verified, the algorithm removes the fourth column off and solves the LS in the same way as it does in the simple circuit case. Removing the fourth column decreases the operations necessary for estimating the considered element and directly results in the parameters of interest. If the fourth column is not removed when the singularity is verified, the resistance RAG is equal to the first element of xAG, while the inductance LAG equals the difference between the second and third elements of xAG.


#### 2.2.2. LL Elements

The algorithm considers that the voltage and current sample windows have a fixed sizeNLL to estimate the LL elements’ positive-sequence resistances and inductances. Therefore, the counter k goes from k=m=n−NLL+1 (oldest sample) to k=n (most recent sample) for the LL elements. At each new sampling, the algorithm incorporates the most recent samples into the windows and discards the oldest ones, as occurs in the sample windows of the LG elements. The calculation of current derivatives is the same as described by (1), with the exception that parameter dLG is replaced by dLL. It is noteworthy that the parameters NLL and dLL do not need to be equal to NLG and dLG.LL elements depend on voltage differences between two phases of the protected line. The LS formulation for estimating the positive-sequence resistance and inductance of element AB is presented below. The other LL elements are analogous.Starting from (2) and considering that the resistances and inductances in phases A and B are, respectively, equal to RAB and LAB, the voltage between phases A and B is given by(12)vABt=vAt−vBt=RAB⋅iAt−iBt+LAB⋅diAtdt−diBtdt.At the discrete-time instantk, (12) becomes(13)vABk=RAB⋅iABk+LAB⋅ΔiABk,where vABk is the difference between vAk and vBk, iABk is the difference between iAk and iBk, and ΔiABk is the difference between ΔiAk and ΔiBk.Therefore, from (13), the LS formulation for determining the resistance RAB and inductance LAB seen by the relay for element AB is given by(14)IΔIAB⋅xAB=VAB,where(15)IΔIAB=iABnΔiABn⋮⋮iABkΔiABk⋮⋮iABmΔiABm,VAB=vABn⋮vABk⋮vABm,xAB=RABLAB.The solution of (14) yields vector xAB. It is worth noting that the LL elements estimation is the same for single-circuit as for double-circuit transmission lines.


### 2.3. LS Solution Techniques

As the algorithm eliminates the linearly dependent columns when the LG elementsIΔI matrices are singular before solving the LS (to reduce computational burden), cases with rank-deficient matrices are not considered during the method’s execution. Thus, the LS can be solved, for example, through the Qℛ decomposition of IΔI, in which the matrix is decomposed into an orthogonal square matrix Q and an upper triangular matrix ℛ, or by determining the pseudoinverse of IΔI [29].The pseudoinverse can be obtained by the singular-value decomposition (SVD) of the original matrix or by the solution of the normal equation (29). The authors tested the proposed algorithm using the two approaches to obtain the pseudoinverse and Qℛ decomposition. The results obtained by Qℛ decomposition and calculating the pseudoinverse through SVD were the same. The results obtained by solving the normal equations, in turn, presented more significant numerical oscillations. As Qℛ decomposition requires fewer operations than SVD [29], it was the chosen LS solution technique.


### 2.4. Fault Probability Calculation through BI

The resistances and inductances estimated through LS usually enter the protection zone faster than the values calculated through voltage and current phasors during a fault. However, the LS estimates present more numerical oscillations than the impedance calculated by phasor-based algorithms, especially under conditions with noise and high-frequency transients. Thus, although the LS estimated values are generally faster in detecting in-zone faults, their accuracy and quality are inferior compared to phasor-calculated values.Therefore, the use of LS values for the distance function presents in average faster fault detection speeds than the phasor-based approaches. However, it is more susceptible to failures since the numerical fluctuations may cause the estimated values to fall in or out of the zone erroneously. Thus, to seek fast fault detection times and increase the trip decisions’ reliability, the presented algorithm uses the BI applied to the LS estimates to determine the in-zone fault occurrence probability over time. The calculation of fault probabilities presents a relatively low increase in the method’s computational burden in comparison with LS solving. Besides, it significantly increases the reliability of the trip decisions since the probability values oscillate considerably less than the LS estimated parameters.BI is based on Bayes’ theorem, which states that the occurrence probability of a given event changes as new evidence is collected [30, 31]. In the proposed algorithm, the event whose probability is to be determined is the occurrence of a fault in the protected zone, and the evidence is whether the LS estimates are inside or outside the zone. It is worth mentioning that each zone has its fault probability calculated individually in multiple zones applications.LetRk and Lk be the LS estimated positive-sequence resistance and inductance for a given element at the discrete-time instant k. From the setting of the protection zone in the R−L plane, defined before the method’s execution, the algorithm checks whether the point Pk, with coordinates Rk,Lk, is located inside or outside the zone. The logical value γk is true if Pk is inside the zone and false otherwise.The calculation of the in-zone fault probability considers a set ofNγ (setting parameter) logical values γk, which compose vector Γ. For the γk values, counter k points to the position of each value in Γ, where γk=m=n−Nγ+1 is the oldest logical value and γk=n is the newest. At every new sampling, the algorithm increments n by one unit and recalculates the fault probability. To do so, it incorporates the new γn value into Γ and discards the oldest one.The setting parametersα and β are, respectively, the probabilities of obtaining a true γk value when known that there is and when there is not a fault in the protected zone. Parameter ξ, in turn, is the initial probability that the in-zone fault exists if the method computes a γk value equal to one. During a fault in the protected zone ℱ, a given γk value presents a Bernoulli probability distribution pγk|ℱ depending on α, as in(16)pγk|ℱ=αγk⋅1−α1−γk.Therefore, the fault likelihood given the observations setΓ, ℒℱ|Γ, is given by the product of all Nγ probability distributions of elements γk, which results in(17)ℒℱ|Γ=α∑k=mnγk⋅1−αNγ−∑k=mnγk.Similarly, the no-faultNℱ likelihood depends on β, as in(18)ℒNℱ|Γ=β∑k=mnγk⋅1−βNγ−∑k=mnγk.Thus, with the faultℱ and no-fault Nℱ likelihoods and with parameter ξ, which is the initial in-zone fault probability for a true γk value (the initial no-fault probability is, therefore, equal to 1−ξ), Bayes’ theorem can be applied to determine the in-zone fault probability given the evidence contained in Γ, pℱ|Γ, as in(19)pℱ|Γ=ξ⋅ℒℱ|Γξ⋅ℒℱ|Γ+1−ξ⋅ℒNℱ|Γ.From (19), the fault probability pℱ|Γ can be rewritten as(20)pℱ|Γ=11+λ,where(21)λ=1−ξξ⋅β⋅1−αα⋅1−β∑k=mnγk⋅1−β1−αNγ.The method calculates the in-zone fault probabilitypℱ|Γ of each element through (20) and (21) at every new sampling. Setting parameters α, β, ξ, and Nγ are the same for all six elements. Each Γ vector, in contrast, is calculated by the algorithm from the respective element’s LS estimated resistances and inductances. The fault probability values pℱ|Γ are the basis for the algorithm’s trip logic, which is presented next.


### 2.5. Trip Logic

For the elaboration of the trip logic, the algorithm considers a minimum fault probability valuep, which is predetermined. If a pℱ|Γ calculated probability value exceeds p, there is a pick-up flag for the considered element. If and only if there is a pick-up flag for Npk consecutive samplings and pℱ|Γ calculations (criterion designed to avoid security failures by detecting the fault Npk consecutive times), the method concludes that there is an in-zone fault for the considered element and sends a trip command to the terminal’s CB and to the other terminal’s IED (in two-terminal schemes). There is one trip output per element, and the trip logic, presented in Figure 2, is the same for all six elements. The algorithm commands a line opening if at least one of the six trip outputs is triggered.Figure 2
Trip logic.


### 2.6. Algorithm Settings

As presented in the previous sections, the proposed algorithm has setting parameters that are fundamental to estimate the elements through LS, determine the in-zone fault probabilitiespℱ|Γ, and verify if these constitute a trip condition.Different values for the parameters may lead to different trip times. Therefore, through successive and extensive tests considering different parameter combinations, the authors determined default settings, which showed to be optimal considering the trip times and the algorithm’s security and dependability, and used them for all presented applications (including the real faults). The algorithm with the default settings obtained fast and secure results in all cases, which indicates that these settings present considerable applicability. Still, other settings may be tested and suggested if deemed necessary. It is noteworthy that the proposed algorithm’s setting parameters are independent from the IED sampling frequency.The authors determined the default values of parametersNLG, dLG, NLL, and dLL (for estimating the apparent resistances and inductances) as, respectively, equal to 8, 1, 9, and 3. With these settings, the values estimated through LS enter the zone during faults quickly, with reduced numerical oscillations compared to settings with smaller sample windows sizes and with few samplings compared to windows close to one full cycle (considering the typical sampling frequencies of modern commercial IEDs, which are generally at least on [kHz] magnitude).The default value of parameterα, which is the probability of obtaining a true γk value during a fault, is 0.95. Preferably, this value should be close to 1 as most values estimated by LS enter the zone during a fault. Still, it cannot be precisely equal to one as not all LS estimated values may enter the zone under fault conditions. Parameter β, which is the probability of obtaining a γk value equal to one in no-fault cases, has a default value of 0.05. Preferably, this value should be close to 0 as most LS-estimated values do not enter the zone during healthy line operation. Nevertheless, it cannot be precisely null, as true γk values can eventually be obtained even in normal operation. The default value of the initial fault probability if the method calculates a unitary γk value ξ, in turn, is 0.90. Even though this value should be preferably high, it cannot be precisely unitary since obtaining a single true γk logic result does not provide absolute certainty that there is a fault in the zone. The default number of γk values needed to calculate the probabilities, Nγ, is equal to 4. This setting proved itself to produce trustworthy and sensitive fault probability values pℱ|Γ, providing fast fault detection. The minimum fault probability p default value is 0.25, while the default value of Npk, which is the number of pick-up flags that must be verified consecutively to trigger the trip output, is also 4.Before the protected line’s energizing, the algorithm requires windows containingmaxNLG,NLL no-fault operation current and voltage values stored in the relay’s memory. The initial values are fundamental for determining the first numerical current derivatives since these calculations require past samples. The method discards the initial values as the protection apparatus performs new samplings and the data windows slide in time.The stored values increase the method’s security during the line’s energizing as they generate fault probability valuespℱ|Γ smaller than threshold p. These values can be collected and stored during real-time execution at standard line operation (to be used in future energizing events) or calculated from the line’s nominal loading conditions. The initial values influence the method’s response if a fault happens before maxNLG,NLL samplings are performed, especially in switching onto fault cases. Therefore, in Section 4, the authors analyzed the influence of the stored initial values in the method’s response to these cases.


### 2.7. Comparison with the DFT-Based Algorithm

The authors compared the proposed algorithm with a self-polarised distance method that calculates the LG and LL impedance elements through voltage and current phasors estimated by the full-cycle DFT [28].The AG element’s impedanceZAG in a simple-circuit line is(22)ZAG=VAGIA+z0−z1/z1⋅I0,where VAG is the voltage phasor between phase A and ground, IA is the current phasor in phase A, I0 is the zero-sequence current phasor of the circuit, and z0 and z1 are, respectively, the zero and positive-sequence impedance of the protected line, in Ω/km.In a double circuit, in turn, where the zero-sequence magnetic coupling between circuits has to be considered,ZAG is given by(23)ZAG=VAGIA+z0−z1/z1⋅I0+3⋅zM0/z1⋅I0p,where I\mbox0p is the zero-sequence current phasor of the parallel circuit and zM0 is the zero-sequence mutual impedance between circuits, in Ω/km. The calculation of the other LG elements (BG and CG) impedance values, both in single and double circuits, is analogous.For element AB, in both single- and double-circuit lines,ZAB is(24)ZAB=VAG−VBGIA−IB.For the other LL elements, BC and CA, the calculation is analogous.There is a pick-up flag for a given element whether the point whose coordinates are the resistanceR and inductance L, determined from the corresponding element’s impedance Z, as in(25)Z=R+j⋅2πf⋅L,where f is the fundamental frequency, is within the zone. The algorithm triggers the trip output of the considered element if it verifies Npk consecutive pick-up flags.The main comparison criterion used by the authors is the protection speed, characterized by the average trip times (interval between fault inception and the trigger of the trip outputs), of each algorithm. The authors also analyzed the stabilization times of each method, which are the intervals between fault inception and the instant when the trip magnitude, which is the probabilitypℱ|Γ in the proposed algorithm and the resistance and inductance in the conventional, stabilizes in a range of ± 5 % of the final measured value. The stabilization times are linked with the trip decision’s reliability since they indicate how long after the fault inception the numerical oscillations related to the trip magnitude are limited to variations of at most 5 %.


## 2.1. Distance Protection Topology

Figure1 depicts a distance protection topology with IEDs in both terminals of the line: local L and remote R. In each IED, the respective phase voltages (measured by the voltage transformer (VT), which usually is a CCVT for high voltage (HV) and EHV lines) and currents (measured by the CT) are sampled at time intervals equal to Ts=1/fs, where fs is the sampling frequency.Figure 1
Two-terminal distance protection topology.Both IEDs execute the proposed algorithm in real time. The first IED that detects a fault sends a trip signal to its controlled circuit breaker (CB) and sends a transfer trip command to the other terminal through a dedicated communication channel. It is worth highlighting that the samples at each terminal do not need to be synchronized and the sampling frequency at each terminal does not necessarily need to be the same. Besides, a two-terminal protection infrastructure is not mandatory for the presented distance protection algorithm since it is possible to implement it at only one end of the line. However, a two-terminal solution increases reliability and fault detection speed.The authors set the primary zones (Zone 1⟶L and R) at 85 %of the lines’ positive-sequence resistance and inductance and considered a Mho characteristic [4] with the apparent inductance on the vertical axis and the apparent resistance on the horizontal axis for the application of the proposed method, in all cases (simulated and actual). Nevertheless, like phasor-based solutions, the presented method allows other zone configurations and characteristics.It is worth mentioning that, especially in one-terminal schemes and for high-impedance faults detection, secondary zones that overreach beyond the remote terminal and have intentional coordination delays can be used [28]. Still, as the evaluation of the proposed algorithm is based mainly on trip times, the authors considered in this research only the Mho characteristic’s first zone [4].


## 2.2. Estimation of Resistances and Inductances through LS

The proposed algorithm relies on six impedance elements: three line-to-ground (LG) elements (AG, BG, and CG) and three line-to-line (LL) elements (AB, BC, and CA). The LS formulations used to estimate the elements’ positive-sequence resistances and inductances (at every new sampling) are presented next. As mentioned earlier, the estimation performed by the presented method is based on [20]. Still, it is adapted to consider the mutual inductances between phases and the mutual inductance between circuits (in double-circuit lines).


### 2.2.1. LG Elements

The algorithm considers that the voltage and current sample windows have a predetermined fixed sizeNLG to estimate the resistance and inductance of each LG element. At every new sampling, the most recent samples are incorporated into the respective windows and the oldest ones are discarded. The derivatives of the phase current signals at a generic discrete-time instant k, where k is a counter that grows from m=n−NLG+1 to n and indicates the position of each sample in the window (k=m⟶ oldest sample, k=n⟶ newest sample), can be calculated numerically as in(1)dit=k⋅Tsdt=Δik=ik−ik−dLGdLG⋅Ts,where dLG is a setting parameter (positive integer) and ik is the kth sample of the current data window. The current derivatives calculation in the presented method is more flexible compared to that presented in [20], in which dLG is always equal to one. Next, the formulations that involved the AG element are presented. The other LG elements are analogous.Single CircuitNeglecting the line capacitance and considering that it is perfectly transposed, the voltage in phase A to ground over timetvAt can be written as(2)vAt=RA⋅iAt+LA⋅diAtdt+LM⋅diBtdt+diCtdt,where RA and LA are, respectively, the resistance and inductance of phase A, LM is the mutual inductance between the circuit phases, and iA, iB, iC are the phase currents over time.Considering the discrete-time condition with the sampled values at a generic discrete-time instantk, (2) becomes(3)vAk=RA⋅iAk+LA⋅ΔiAk+LM⋅ΔiBk+ΔiCk.Since the voltage and current windows haveNLG positions for estimating the LG elements, the LS formulation for the AG element at every new sampling can be written from (3) as(4)IΔIAG⋅xAG=VA,where(5)IΔIAG=iAnΔiAnΔiBn+ΔiCn⋮⋮⋮iAkΔiAkΔiBk+ΔiCk⋮⋮⋮iAmΔiAmΔiBm+ΔiCm,(6)VA=vAn⋮vAk⋮vAm,xAG=RALALM.The LS solution of (4) yields vector xAG. Thus, the positive-sequence resistance of element AG RAG seen by the relay is equal to RA, while the inductance of element AG LAG seen by the relay is given by the difference between the inductance of phase ALA and the mutual inductance LM, as in(7)RAG=R,LAG=LA−LM.IfΔBk+ΔCk=−ΔAk for all instants k between m and n, the line is balanced and IΔIAG is singular since the third column is the opposite of the second. In this case, the third column of the matrix can be eliminated for the LS solution, which directly results in the AG element’s positive-sequence resistance and inductance (RAG and LAG). If the third column is not removed from the matrix, the estimated values of LA and LM respect the dependence relationship between the respective columns and have absolute values equal to the half of LAGLA=−LM=LAG/2. Thus, the verification of the mentioned singularity may be interesting from a computational point of view, since that if the second and third columns of IΔIAG are linearly dependent, the removal of the third column reduces the number of floating-point operations necessary to estimate the considered element.Double Circuit.In double-circuit lines, the mutual zero-sequence inductance between the circuitsLM0 has to be considered. Thus, (2) can be rewritten through(8)vAt=RA⋅iAt+LA⋅diAtdt+LM⋅diBtdt+diCtdt+3⋅LM0⋅di0ptdt,where(9)i0pt=iApt+iBpt+iCpt3,is the zero-sequence current and iAp,iBp,iCp are the phase currents of the parallel circuit over time. Therefore, for the inclusion of the zero-sequence magnetic coupling between circuits in discrete-time, (3) becomes(10)vAk=RA⋅iAk+LA⋅ΔiAk+LM⋅ΔiBk+ΔiCk+3⋅LM0⋅Δi0pk,where(11)Δi0pk=ΔiApk+ΔiBpk+ΔiCpk3,is the parallel circuit zero-sequence current derivative, calculated from the numerical derivatives of the respective phase currents. At every new sampling, the IED of each circuit (on the same bus) exchange the zero-sequence current derivative values so that the zero-sequence coupling is considered.For the LS formulation, the values of3⋅Δi0pk are added to matrix IΔIAG, forming a new column. Vector xAG now has four elements, the fourth being the estimated value of LM0. Vector VA, in turn, remains unchanged. After the LS solution, the calculation of the considered element’s resistance and inductance is the same as for the simple circuit Case (7).IfΔiAk+ΔiBk+ΔiCk=3⋅Δi0pk for all instants k between m and n, the derivatives of the zero-sequence currents in the two circuits are equal and IΔIAG is singular since the fourth column is the sum of the second with the third. If such singularity is verified, the algorithm removes the fourth column off and solves the LS in the same way as it does in the simple circuit case. Removing the fourth column decreases the operations necessary for estimating the considered element and directly results in the parameters of interest. If the fourth column is not removed when the singularity is verified, the resistance RAG is equal to the first element of xAG, while the inductance LAG equals the difference between the second and third elements of xAG.


### 2.2.2. LL Elements

The algorithm considers that the voltage and current sample windows have a fixed sizeNLL to estimate the LL elements’ positive-sequence resistances and inductances. Therefore, the counter k goes from k=m=n−NLL+1 (oldest sample) to k=n (most recent sample) for the LL elements. At each new sampling, the algorithm incorporates the most recent samples into the windows and discards the oldest ones, as occurs in the sample windows of the LG elements. The calculation of current derivatives is the same as described by (1), with the exception that parameter dLG is replaced by dLL. It is noteworthy that the parameters NLL and dLL do not need to be equal to NLG and dLG.LL elements depend on voltage differences between two phases of the protected line. The LS formulation for estimating the positive-sequence resistance and inductance of element AB is presented below. The other LL elements are analogous.Starting from (2) and considering that the resistances and inductances in phases A and B are, respectively, equal to RAB and LAB, the voltage between phases A and B is given by(12)vABt=vAt−vBt=RAB⋅iAt−iBt+LAB⋅diAtdt−diBtdt.At the discrete-time instantk, (12) becomes(13)vABk=RAB⋅iABk+LAB⋅ΔiABk,where vABk is the difference between vAk and vBk, iABk is the difference between iAk and iBk, and ΔiABk is the difference between ΔiAk and ΔiBk.Therefore, from (13), the LS formulation for determining the resistance RAB and inductance LAB seen by the relay for element AB is given by(14)IΔIAB⋅xAB=VAB,where(15)IΔIAB=iABnΔiABn⋮⋮iABkΔiABk⋮⋮iABmΔiABm,VAB=vABn⋮vABk⋮vABm,xAB=RABLAB.The solution of (14) yields vector xAB. It is worth noting that the LL elements estimation is the same for single-circuit as for double-circuit transmission lines.


## 2.2.1. LG Elements

The algorithm considers that the voltage and current sample windows have a predetermined fixed sizeNLG to estimate the resistance and inductance of each LG element. At every new sampling, the most recent samples are incorporated into the respective windows and the oldest ones are discarded. The derivatives of the phase current signals at a generic discrete-time instant k, where k is a counter that grows from m=n−NLG+1 to n and indicates the position of each sample in the window (k=m⟶ oldest sample, k=n⟶ newest sample), can be calculated numerically as in(1)dit=k⋅Tsdt=Δik=ik−ik−dLGdLG⋅Ts,where dLG is a setting parameter (positive integer) and ik is the kth sample of the current data window. The current derivatives calculation in the presented method is more flexible compared to that presented in [20], in which dLG is always equal to one. Next, the formulations that involved the AG element are presented. The other LG elements are analogous.Single CircuitNeglecting the line capacitance and considering that it is perfectly transposed, the voltage in phase A to ground over timetvAt can be written as(2)vAt=RA⋅iAt+LA⋅diAtdt+LM⋅diBtdt+diCtdt,where RA and LA are, respectively, the resistance and inductance of phase A, LM is the mutual inductance between the circuit phases, and iA, iB, iC are the phase currents over time.Considering the discrete-time condition with the sampled values at a generic discrete-time instantk, (2) becomes(3)vAk=RA⋅iAk+LA⋅ΔiAk+LM⋅ΔiBk+ΔiCk.Since the voltage and current windows haveNLG positions for estimating the LG elements, the LS formulation for the AG element at every new sampling can be written from (3) as(4)IΔIAG⋅xAG=VA,where(5)IΔIAG=iAnΔiAnΔiBn+ΔiCn⋮⋮⋮iAkΔiAkΔiBk+ΔiCk⋮⋮⋮iAmΔiAmΔiBm+ΔiCm,(6)VA=vAn⋮vAk⋮vAm,xAG=RALALM.The LS solution of (4) yields vector xAG. Thus, the positive-sequence resistance of element AG RAG seen by the relay is equal to RA, while the inductance of element AG LAG seen by the relay is given by the difference between the inductance of phase ALA and the mutual inductance LM, as in(7)RAG=R,LAG=LA−LM.IfΔBk+ΔCk=−ΔAk for all instants k between m and n, the line is balanced and IΔIAG is singular since the third column is the opposite of the second. In this case, the third column of the matrix can be eliminated for the LS solution, which directly results in the AG element’s positive-sequence resistance and inductance (RAG and LAG). If the third column is not removed from the matrix, the estimated values of LA and LM respect the dependence relationship between the respective columns and have absolute values equal to the half of LAGLA=−LM=LAG/2. Thus, the verification of the mentioned singularity may be interesting from a computational point of view, since that if the second and third columns of IΔIAG are linearly dependent, the removal of the third column reduces the number of floating-point operations necessary to estimate the considered element.Double Circuit.In double-circuit lines, the mutual zero-sequence inductance between the circuitsLM0 has to be considered. Thus, (2) can be rewritten through(8)vAt=RA⋅iAt+LA⋅diAtdt+LM⋅diBtdt+diCtdt+3⋅LM0⋅di0ptdt,where(9)i0pt=iApt+iBpt+iCpt3,is the zero-sequence current and iAp,iBp,iCp are the phase currents of the parallel circuit over time. Therefore, for the inclusion of the zero-sequence magnetic coupling between circuits in discrete-time, (3) becomes(10)vAk=RA⋅iAk+LA⋅ΔiAk+LM⋅ΔiBk+ΔiCk+3⋅LM0⋅Δi0pk,where(11)Δi0pk=ΔiApk+ΔiBpk+ΔiCpk3,is the parallel circuit zero-sequence current derivative, calculated from the numerical derivatives of the respective phase currents. At every new sampling, the IED of each circuit (on the same bus) exchange the zero-sequence current derivative values so that the zero-sequence coupling is considered.For the LS formulation, the values of3⋅Δi0pk are added to matrix IΔIAG, forming a new column. Vector xAG now has four elements, the fourth being the estimated value of LM0. Vector VA, in turn, remains unchanged. After the LS solution, the calculation of the considered element’s resistance and inductance is the same as for the simple circuit Case (7).IfΔiAk+ΔiBk+ΔiCk=3⋅Δi0pk for all instants k between m and n, the derivatives of the zero-sequence currents in the two circuits are equal and IΔIAG is singular since the fourth column is the sum of the second with the third. If such singularity is verified, the algorithm removes the fourth column off and solves the LS in the same way as it does in the simple circuit case. Removing the fourth column decreases the operations necessary for estimating the considered element and directly results in the parameters of interest. If the fourth column is not removed when the singularity is verified, the resistance RAG is equal to the first element of xAG, while the inductance LAG equals the difference between the second and third elements of xAG.


## 2.2.2. LL Elements

The algorithm considers that the voltage and current sample windows have a fixed sizeNLL to estimate the LL elements’ positive-sequence resistances and inductances. Therefore, the counter k goes from k=m=n−NLL+1 (oldest sample) to k=n (most recent sample) for the LL elements. At each new sampling, the algorithm incorporates the most recent samples into the windows and discards the oldest ones, as occurs in the sample windows of the LG elements. The calculation of current derivatives is the same as described by (1), with the exception that parameter dLG is replaced by dLL. It is noteworthy that the parameters NLL and dLL do not need to be equal to NLG and dLG.LL elements depend on voltage differences between two phases of the protected line. The LS formulation for estimating the positive-sequence resistance and inductance of element AB is presented below. The other LL elements are analogous.Starting from (2) and considering that the resistances and inductances in phases A and B are, respectively, equal to RAB and LAB, the voltage between phases A and B is given by(12)vABt=vAt−vBt=RAB⋅iAt−iBt+LAB⋅diAtdt−diBtdt.At the discrete-time instantk, (12) becomes(13)vABk=RAB⋅iABk+LAB⋅ΔiABk,where vABk is the difference between vAk and vBk, iABk is the difference between iAk and iBk, and ΔiABk is the difference between ΔiAk and ΔiBk.Therefore, from (13), the LS formulation for determining the resistance RAB and inductance LAB seen by the relay for element AB is given by(14)IΔIAB⋅xAB=VAB,where(15)IΔIAB=iABnΔiABn⋮⋮iABkΔiABk⋮⋮iABmΔiABm,VAB=vABn⋮vABk⋮vABm,xAB=RABLAB.The solution of (14) yields vector xAB. It is worth noting that the LL elements estimation is the same for single-circuit as for double-circuit transmission lines.


## 2.3. LS Solution Techniques

As the algorithm eliminates the linearly dependent columns when the LG elementsIΔI matrices are singular before solving the LS (to reduce computational burden), cases with rank-deficient matrices are not considered during the method’s execution. Thus, the LS can be solved, for example, through the Qℛ decomposition of IΔI, in which the matrix is decomposed into an orthogonal square matrix Q and an upper triangular matrix ℛ, or by determining the pseudoinverse of IΔI [29].The pseudoinverse can be obtained by the singular-value decomposition (SVD) of the original matrix or by the solution of the normal equation (29). The authors tested the proposed algorithm using the two approaches to obtain the pseudoinverse and Qℛ decomposition. The results obtained by Qℛ decomposition and calculating the pseudoinverse through SVD were the same. The results obtained by solving the normal equations, in turn, presented more significant numerical oscillations. As Qℛ decomposition requires fewer operations than SVD [29], it was the chosen LS solution technique.


## 2.4. Fault Probability Calculation through BI

The resistances and inductances estimated through LS usually enter the protection zone faster than the values calculated through voltage and current phasors during a fault. However, the LS estimates present more numerical oscillations than the impedance calculated by phasor-based algorithms, especially under conditions with noise and high-frequency transients. Thus, although the LS estimated values are generally faster in detecting in-zone faults, their accuracy and quality are inferior compared to phasor-calculated values.Therefore, the use of LS values for the distance function presents in average faster fault detection speeds than the phasor-based approaches. However, it is more susceptible to failures since the numerical fluctuations may cause the estimated values to fall in or out of the zone erroneously. Thus, to seek fast fault detection times and increase the trip decisions’ reliability, the presented algorithm uses the BI applied to the LS estimates to determine the in-zone fault occurrence probability over time. The calculation of fault probabilities presents a relatively low increase in the method’s computational burden in comparison with LS solving. Besides, it significantly increases the reliability of the trip decisions since the probability values oscillate considerably less than the LS estimated parameters.BI is based on Bayes’ theorem, which states that the occurrence probability of a given event changes as new evidence is collected [30, 31]. In the proposed algorithm, the event whose probability is to be determined is the occurrence of a fault in the protected zone, and the evidence is whether the LS estimates are inside or outside the zone. It is worth mentioning that each zone has its fault probability calculated individually in multiple zones applications.LetRk and Lk be the LS estimated positive-sequence resistance and inductance for a given element at the discrete-time instant k. From the setting of the protection zone in the R−L plane, defined before the method’s execution, the algorithm checks whether the point Pk, with coordinates Rk,Lk, is located inside or outside the zone. The logical value γk is true if Pk is inside the zone and false otherwise.The calculation of the in-zone fault probability considers a set ofNγ (setting parameter) logical values γk, which compose vector Γ. For the γk values, counter k points to the position of each value in Γ, where γk=m=n−Nγ+1 is the oldest logical value and γk=n is the newest. At every new sampling, the algorithm increments n by one unit and recalculates the fault probability. To do so, it incorporates the new γn value into Γ and discards the oldest one.The setting parametersα and β are, respectively, the probabilities of obtaining a true γk value when known that there is and when there is not a fault in the protected zone. Parameter ξ, in turn, is the initial probability that the in-zone fault exists if the method computes a γk value equal to one. During a fault in the protected zone ℱ, a given γk value presents a Bernoulli probability distribution pγk|ℱ depending on α, as in(16)pγk|ℱ=αγk⋅1−α1−γk.Therefore, the fault likelihood given the observations setΓ, ℒℱ|Γ, is given by the product of all Nγ probability distributions of elements γk, which results in(17)ℒℱ|Γ=α∑k=mnγk⋅1−αNγ−∑k=mnγk.Similarly, the no-faultNℱ likelihood depends on β, as in(18)ℒNℱ|Γ=β∑k=mnγk⋅1−βNγ−∑k=mnγk.Thus, with the faultℱ and no-fault Nℱ likelihoods and with parameter ξ, which is the initial in-zone fault probability for a true γk value (the initial no-fault probability is, therefore, equal to 1−ξ), Bayes’ theorem can be applied to determine the in-zone fault probability given the evidence contained in Γ, pℱ|Γ, as in(19)pℱ|Γ=ξ⋅ℒℱ|Γξ⋅ℒℱ|Γ+1−ξ⋅ℒNℱ|Γ.From (19), the fault probability pℱ|Γ can be rewritten as(20)pℱ|Γ=11+λ,where(21)λ=1−ξξ⋅β⋅1−αα⋅1−β∑k=mnγk⋅1−β1−αNγ.The method calculates the in-zone fault probabilitypℱ|Γ of each element through (20) and (21) at every new sampling. Setting parameters α, β, ξ, and Nγ are the same for all six elements. Each Γ vector, in contrast, is calculated by the algorithm from the respective element’s LS estimated resistances and inductances. The fault probability values pℱ|Γ are the basis for the algorithm’s trip logic, which is presented next.


## 2.5. Trip Logic

For the elaboration of the trip logic, the algorithm considers a minimum fault probability valuep, which is predetermined. If a pℱ|Γ calculated probability value exceeds p, there is a pick-up flag for the considered element. If and only if there is a pick-up flag for Npk consecutive samplings and pℱ|Γ calculations (criterion designed to avoid security failures by detecting the fault Npk consecutive times), the method concludes that there is an in-zone fault for the considered element and sends a trip command to the terminal’s CB and to the other terminal’s IED (in two-terminal schemes). There is one trip output per element, and the trip logic, presented in Figure 2, is the same for all six elements. The algorithm commands a line opening if at least one of the six trip outputs is triggered.Figure 2
Trip logic.


## 2.6. Algorithm Settings

As presented in the previous sections, the proposed algorithm has setting parameters that are fundamental to estimate the elements through LS, determine the in-zone fault probabilitiespℱ|Γ, and verify if these constitute a trip condition.Different values for the parameters may lead to different trip times. Therefore, through successive and extensive tests considering different parameter combinations, the authors determined default settings, which showed to be optimal considering the trip times and the algorithm’s security and dependability, and used them for all presented applications (including the real faults). The algorithm with the default settings obtained fast and secure results in all cases, which indicates that these settings present considerable applicability. Still, other settings may be tested and suggested if deemed necessary. It is noteworthy that the proposed algorithm’s setting parameters are independent from the IED sampling frequency.The authors determined the default values of parametersNLG, dLG, NLL, and dLL (for estimating the apparent resistances and inductances) as, respectively, equal to 8, 1, 9, and 3. With these settings, the values estimated through LS enter the zone during faults quickly, with reduced numerical oscillations compared to settings with smaller sample windows sizes and with few samplings compared to windows close to one full cycle (considering the typical sampling frequencies of modern commercial IEDs, which are generally at least on [kHz] magnitude).The default value of parameterα, which is the probability of obtaining a true γk value during a fault, is 0.95. Preferably, this value should be close to 1 as most values estimated by LS enter the zone during a fault. Still, it cannot be precisely equal to one as not all LS estimated values may enter the zone under fault conditions. Parameter β, which is the probability of obtaining a γk value equal to one in no-fault cases, has a default value of 0.05. Preferably, this value should be close to 0 as most LS-estimated values do not enter the zone during healthy line operation. Nevertheless, it cannot be precisely null, as true γk values can eventually be obtained even in normal operation. The default value of the initial fault probability if the method calculates a unitary γk value ξ, in turn, is 0.90. Even though this value should be preferably high, it cannot be precisely unitary since obtaining a single true γk logic result does not provide absolute certainty that there is a fault in the zone. The default number of γk values needed to calculate the probabilities, Nγ, is equal to 4. This setting proved itself to produce trustworthy and sensitive fault probability values pℱ|Γ, providing fast fault detection. The minimum fault probability p default value is 0.25, while the default value of Npk, which is the number of pick-up flags that must be verified consecutively to trigger the trip output, is also 4.Before the protected line’s energizing, the algorithm requires windows containingmaxNLG,NLL no-fault operation current and voltage values stored in the relay’s memory. The initial values are fundamental for determining the first numerical current derivatives since these calculations require past samples. The method discards the initial values as the protection apparatus performs new samplings and the data windows slide in time.The stored values increase the method’s security during the line’s energizing as they generate fault probability valuespℱ|Γ smaller than threshold p. These values can be collected and stored during real-time execution at standard line operation (to be used in future energizing events) or calculated from the line’s nominal loading conditions. The initial values influence the method’s response if a fault happens before maxNLG,NLL samplings are performed, especially in switching onto fault cases. Therefore, in Section 4, the authors analyzed the influence of the stored initial values in the method’s response to these cases.


## 2.7. Comparison with the DFT-Based Algorithm

The authors compared the proposed algorithm with a self-polarised distance method that calculates the LG and LL impedance elements through voltage and current phasors estimated by the full-cycle DFT [28].The AG element’s impedanceZAG in a simple-circuit line is(22)ZAG=VAGIA+z0−z1/z1⋅I0,where VAG is the voltage phasor between phase A and ground, IA is the current phasor in phase A, I0 is the zero-sequence current phasor of the circuit, and z0 and z1 are, respectively, the zero and positive-sequence impedance of the protected line, in Ω/km.In a double circuit, in turn, where the zero-sequence magnetic coupling between circuits has to be considered,ZAG is given by(23)ZAG=VAGIA+z0−z1/z1⋅I0+3⋅zM0/z1⋅I0p,where I\mbox0p is the zero-sequence current phasor of the parallel circuit and zM0 is the zero-sequence mutual impedance between circuits, in Ω/km. The calculation of the other LG elements (BG and CG) impedance values, both in single and double circuits, is analogous.For element AB, in both single- and double-circuit lines,ZAB is(24)ZAB=VAG−VBGIA−IB.For the other LL elements, BC and CA, the calculation is analogous.There is a pick-up flag for a given element whether the point whose coordinates are the resistanceR and inductance L, determined from the corresponding element’s impedance Z, as in(25)Z=R+j⋅2πf⋅L,where f is the fundamental frequency, is within the zone. The algorithm triggers the trip output of the considered element if it verifies Npk consecutive pick-up flags.The main comparison criterion used by the authors is the protection speed, characterized by the average trip times (interval between fault inception and the trigger of the trip outputs), of each algorithm. The authors also analyzed the stabilization times of each method, which are the intervals between fault inception and the instant when the trip magnitude, which is the probabilitypℱ|Γ in the proposed algorithm and the resistance and inductance in the conventional, stabilizes in a range of ± 5 % of the final measured value. The stabilization times are linked with the trip decision’s reliability since they indicate how long after the fault inception the numerical oscillations related to the trip magnitude are limited to variations of at most 5 %.


## 3. Algorithm Application and Results

### 3.1. Description of the Modeled Transmission Lines

The authors modeled five existing transmission lines in ATP for the fault simulations. Each line impedance was determined at 60 [Hz]. Table1 presents the characteristics of each line. The zero-sequence reactance between the two circuits in line 1 is XM0  = 0.39 [Ω/km].Table 1
Modeled lines characteristics.
IdV (kV)TypeL (km)Z1 (Ω/km)Z0 (Ω/km)1138Double740.21+j0.500.51+j1.742230Simple800.07+j0.490.40+j1.613345Simple1450.02+j0.350.39+j1.124500Simple1810.01+j0.310.38+j1.165765Simple2650.01+j0.330.36+j1.09The authors simulated single-phase (AG), double-phase (B-C), double-phase to ground (B-C-G), and three-phase faults (A-B-C) on all lines. Fault application points were between1/6 and 5/6 of each line’s length, fault inception instants were 0.017 and 0.020 [s], and fault resistances were 0.01, 1, 2, and 5 Ω for LL faults (B-C and A-B-C) and 0.01, 1, 2, 5, 20, and 50 Ω for LG faults (AG and B-C-G). The algorithms were applied in the two terminals of each line, according to Figure 1. Light and heavy loading conditions were considered for each line by varying the load angle of the remote terminal’s equivalent source. The authors resampled the voltage and current signals from the simulations to 1920 [Hz] (32 samples per cycle of 60 [Hz]). In total, 4080 cases were generated and simulated. Table 2 shows the line terminals’ three-phase (3-Ph) and single-phase (1-Ph) fault levels, in [GVA].Table 2
Transmission lines fault levels (GVA).
Id3-PhL3-PhR1-PhL1-PhR13.90∠86°1.99∠75°3.12∠ 88 °1.59∠ 81 °21.96∠81°4.82∠86°1.56∠ 85 °3.86∠ 88 °31.09∠83°1.87∠85°0.87∠ 85 °1.49∠ 89 °415.82∠87°14.82∠88°12.66∠ 89 °11.86∠ 89 °520.14∠87°22.34∠88°16.11∠ 89 °17.87∠ 89 °


### 3.2. Results

Table3 presents the average trip times of the algorithms in [ms] for both LG and LL elements. The proposed method’s average trip time (considering all 4080 cases) was equal to 7.46 [ms]. For the same cases, the average trip time of the algorithm with only the LS-estimated resistances and inductances (same trip logic of the DFT-conventional algorithm) was 7.77 [ms], while the average trip times of the DFT-conventional algorithm was 13.48 [ms]. The average trip times of the algorithms in faults close to the limit of the first zone (distance of 83.3 % 5/6 from the considered terminal) were respectively equal to 8.25 [ms] (LS + BI), 15.62 [ms] (LS only) and 17.71 [ms] (DFT-conventional). This result indicates that the use of BI significantly improves the identification of faults close to the first zone boundary, in which the estimated resistances and inductances enter and exit the zone more often before stabilizing.Table 3
Average trip times (ms).
IdTermLS + BI-LGDFT-LGLS + BI-LLDFT-LL1L7.5513.217.5113.73R6.8011.796.3811.802L6.9814.297.5514.53R6.4112.687.9812.603L9.6416.647.5817.02R7.2914.747.4014.574L7.4013.367.3212.87R6.6712.517.6412.475L7.9413.327.8312.71R7.5512.427.8212.23The average stabilization times were equal to 6.02 [ms] and 29.54 [ms] for the proposed (LS + BI) and the DFT-conventional algorithm, respectively. The average stabilization time of the algorithm with only the resistances and inductances estimated by LS was equal to 76.88 [ms]. This result indicates that the probabilitiespℱ|Γ calculated by BI were responsible for substantially reducing the numerical oscillations of the LS estimated values and thus improving the reliability of the trip decisions.


### 3.3. Comparison with BI-Based Differential Protection

#### 3.3.1. Trip times

For the same cases, the average trip time of the differential protection algorithm based on BI [25] was equal to 6.78 [ms]. Thus, both methods present answers with the same order of magnitude. However, the 6.78 [ms] average trip time did not consider communication delays so that in many practical situations, the proposed distance protection algorithm would be even faster than the compared differential scheme.


### 3.4. Regards concerning High-Impedance Faults

Additionally, the authors simulated both the proposed and the differential algorithms [25] in cases featuring high-impedance (500 and 1000 Ω) internal faults. In these cases, the proposed method only detects faults if zone 3 (see [26, 27]) is used, which yields trip times of at least 800 to 1000 [ms] [26]. Still, the differential algorithm detects high impedance faults with average trip times of 7.67 [ms].The proposed algorithm is capable of detecting faults with resistances of up to 150Ω in the first zone (average trip times of 18.87 [ms] with such resistance), which is a typical tower foot resistance value in the Brazilian Interconnected Grid.


#### 3.4.1. Advantages and Disadvantages of Each Approach

Therefore, the proposed distance scheme is advantageous over the differential protection approach as it dismisses current samples exchange between terminals, synchronization, and yields trip times with the same order of magnitude. In longer lines, the presented distance approach will certainly be faster due to the communication delays.Still, it is noteworthy that the differential approach, despite requiring more complex infrastructure and protection apparatus, can detect internal faults with very high resistance and does not need voltage sampled values nor the line’s parameters.


## 3.1. Description of the Modeled Transmission Lines

The authors modeled five existing transmission lines in ATP for the fault simulations. Each line impedance was determined at 60 [Hz]. Table1 presents the characteristics of each line. The zero-sequence reactance between the two circuits in line 1 is XM0  = 0.39 [Ω/km].Table 1
Modeled lines characteristics.
IdV (kV)TypeL (km)Z1 (Ω/km)Z0 (Ω/km)1138Double740.21+j0.500.51+j1.742230Simple800.07+j0.490.40+j1.613345Simple1450.02+j0.350.39+j1.124500Simple1810.01+j0.310.38+j1.165765Simple2650.01+j0.330.36+j1.09The authors simulated single-phase (AG), double-phase (B-C), double-phase to ground (B-C-G), and three-phase faults (A-B-C) on all lines. Fault application points were between1/6 and 5/6 of each line’s length, fault inception instants were 0.017 and 0.020 [s], and fault resistances were 0.01, 1, 2, and 5 Ω for LL faults (B-C and A-B-C) and 0.01, 1, 2, 5, 20, and 50 Ω for LG faults (AG and B-C-G). The algorithms were applied in the two terminals of each line, according to Figure 1. Light and heavy loading conditions were considered for each line by varying the load angle of the remote terminal’s equivalent source. The authors resampled the voltage and current signals from the simulations to 1920 [Hz] (32 samples per cycle of 60 [Hz]). In total, 4080 cases were generated and simulated. Table 2 shows the line terminals’ three-phase (3-Ph) and single-phase (1-Ph) fault levels, in [GVA].Table 2
Transmission lines fault levels (GVA).
Id3-PhL3-PhR1-PhL1-PhR13.90∠86°1.99∠75°3.12∠ 88 °1.59∠ 81 °21.96∠81°4.82∠86°1.56∠ 85 °3.86∠ 88 °31.09∠83°1.87∠85°0.87∠ 85 °1.49∠ 89 °415.82∠87°14.82∠88°12.66∠ 89 °11.86∠ 89 °520.14∠87°22.34∠88°16.11∠ 89 °17.87∠ 89 °


## 3.2. Results

Table3 presents the average trip times of the algorithms in [ms] for both LG and LL elements. The proposed method’s average trip time (considering all 4080 cases) was equal to 7.46 [ms]. For the same cases, the average trip time of the algorithm with only the LS-estimated resistances and inductances (same trip logic of the DFT-conventional algorithm) was 7.77 [ms], while the average trip times of the DFT-conventional algorithm was 13.48 [ms]. The average trip times of the algorithms in faults close to the limit of the first zone (distance of 83.3 % 5/6 from the considered terminal) were respectively equal to 8.25 [ms] (LS + BI), 15.62 [ms] (LS only) and 17.71 [ms] (DFT-conventional). This result indicates that the use of BI significantly improves the identification of faults close to the first zone boundary, in which the estimated resistances and inductances enter and exit the zone more often before stabilizing.Table 3
Average trip times (ms).
IdTermLS + BI-LGDFT-LGLS + BI-LLDFT-LL1L7.5513.217.5113.73R6.8011.796.3811.802L6.9814.297.5514.53R6.4112.687.9812.603L9.6416.647.5817.02R7.2914.747.4014.574L7.4013.367.3212.87R6.6712.517.6412.475L7.9413.327.8312.71R7.5512.427.8212.23The average stabilization times were equal to 6.02 [ms] and 29.54 [ms] for the proposed (LS + BI) and the DFT-conventional algorithm, respectively. The average stabilization time of the algorithm with only the resistances and inductances estimated by LS was equal to 76.88 [ms]. This result indicates that the probabilitiespℱ|Γ calculated by BI were responsible for substantially reducing the numerical oscillations of the LS estimated values and thus improving the reliability of the trip decisions.


## 3.3. Comparison with BI-Based Differential Protection

### 3.3.1. Trip times

For the same cases, the average trip time of the differential protection algorithm based on BI [25] was equal to 6.78 [ms]. Thus, both methods present answers with the same order of magnitude. However, the 6.78 [ms] average trip time did not consider communication delays so that in many practical situations, the proposed distance protection algorithm would be even faster than the compared differential scheme.


## 3.3.1. Trip times

For the same cases, the average trip time of the differential protection algorithm based on BI [25] was equal to 6.78 [ms]. Thus, both methods present answers with the same order of magnitude. However, the 6.78 [ms] average trip time did not consider communication delays so that in many practical situations, the proposed distance protection algorithm would be even faster than the compared differential scheme.


## 3.4. Regards concerning High-Impedance Faults

Additionally, the authors simulated both the proposed and the differential algorithms [25] in cases featuring high-impedance (500 and 1000 Ω) internal faults. In these cases, the proposed method only detects faults if zone 3 (see [26, 27]) is used, which yields trip times of at least 800 to 1000 [ms] [26]. Still, the differential algorithm detects high impedance faults with average trip times of 7.67 [ms].The proposed algorithm is capable of detecting faults with resistances of up to 150Ω in the first zone (average trip times of 18.87 [ms] with such resistance), which is a typical tower foot resistance value in the Brazilian Interconnected Grid.


### 3.4.1. Advantages and Disadvantages of Each Approach

Therefore, the proposed distance scheme is advantageous over the differential protection approach as it dismisses current samples exchange between terminals, synchronization, and yields trip times with the same order of magnitude. In longer lines, the presented distance approach will certainly be faster due to the communication delays.Still, it is noteworthy that the differential approach, despite requiring more complex infrastructure and protection apparatus, can detect internal faults with very high resistance and does not need voltage sampled values nor the line’s parameters.


## 3.4.1. Advantages and Disadvantages of Each Approach

Therefore, the proposed distance scheme is advantageous over the differential protection approach as it dismisses current samples exchange between terminals, synchronization, and yields trip times with the same order of magnitude. In longer lines, the presented distance approach will certainly be faster due to the communication delays.Still, it is noteworthy that the differential approach, despite requiring more complex infrastructure and protection apparatus, can detect internal faults with very high resistance and does not need voltage sampled values nor the line’s parameters.


## 4. Sensitivity Analysis

### 4.1. Noisy Measurements

The authors applied Gaussian noises with an average signal-to-noise ratio (SNR) of 15 [dB] to both voltage and current signals from lines 2 and 5 to verify noisy measurements’ effect on the algorithms’ responses. Figure3 presents the normal and noisy current signals for a three-phase fault online 2, and Figure 4 illustrates the probabilities pℱ|Γ for the CG element in both cases highlighting the methods’ trip instants. The average trip times in [ms] obtained in the L terminal of the lines in the noisy situation are presented in Table 4.Figure 3
Current signal with and without 15 [dB] Gaussian noise.Figure 4
Fault probabilities with and without noise over time.Table 4
Average trip timesL [ms]—Gaussian noise.
IdLS + BI-LGDFT-LGLS + BI-LLDFT-LL29.2014.297.9014.53510.3214.578.4114.38


### 4.2. CT Saturation and CCVT Transients

The authors modeled at the local terminal of line 2 a CCVT with the parameters obtained in [32] through an actual 230 [kV] equipment and a CT with 1200 : 5 transformation ratio and magnetisation characteristic based on [33] to verify the effects of CT saturation and CCVT transients on the trip times.Figure5 illustrates the saturated and normal currents for a three-phase fault, while Figure 6 presents the voltage with and without the CCVT for the same case. Figure 7, in turn, presents the AB element’s fault probability pℱ|Γ for this case with and without the CT and CCVT effects. Finally, Table 5 presents the average trip times with (1) and without (0) the instrumentation transformers.Figure 5
Normal and saturated current signal.Figure 6
Voltage signal with and without CCVT transients.Figure 7
Fault probabilities with and without the CT + CCVT effects.Table 5
Average trip timesL [ms]—instrumentation transformers.
CCVTCTLS + BI-LGDFT-LGLS + BI-LLDFT-LL107.1414.387.8114.55018.0525.8312.7827.22118.7626.2712.8727.50


### 4.3. Close-In Faults

For the algorithm’s evaluation in close-in fault events, the authors simulated the faults in line 2 at 2 [km] from the local terminal, which corresponds to 2.5 % of the line’s total length. Conditions with and without the instrumentation transformers were considered. Table6 presents the average trip times.Table 6
Average trip timesL [ms]—close-in faults.
CCVTCTLS + BI-LGDFT-LGLS + BI-LLDFT-LL004.749.254.698.85115.7814.204.7914.79


### 4.4. Line Energizing onto Faults

As mentioned in Section2.5, the initial voltage and current values, stored in the IED’s memory before line energizing, can be collected during standard line operation for future use or synthesized from the nominal loading conditions.The initial voltage and current stored values may influence the trip times, especially in switch onto fault conditions. Thus, to test their influence in these conditions, the authors simulated faults at line 2 energizing (by terminalL) and tested the algorithm with the following sets of maxNLG,NLL  = 9 (considering the default settings) initial values sets (1–3).(i)
1: collected in past standard operation(ii)
2: synthetic light loading values(iii)
3: synthetic heavy loading valuesTable7 presents the average trip times in [ms] for the three data sets. The conventional DFT-based algorithm does not need initial values. However, the results show that the proposed algorithm is faster with all considered data sets.Table 7
Average trip timesL [ms]—switch onto faults.
Data setLS + BI-LGDFT-LGLS + BI-LLDFT-LL110.3017.9510.3720.63210.4417.9510.3520.63310.9517.959.7620.63


### 4.5. System Frequency Variations

To evaluate the performance of the algorithms under operating conditions with frequencies other than the fundamental, which may occur due to variations in the system’s generation and load levels, the authors simulated the faults in line 2 considering the frequency range of 58–62 [Hz].The authors did not observe any variations in the trip times of the proposed algorithm and verified variations up to 1 [ms] in the DFT-conventional algorithm trip times. These variations happen because the algorithm based on phasors calculates reactance values at the fundamental frequency.


### 4.6. Uncertainties and Errors in Line Parameters

Generally, the protected line’s impedance values are determined through calculations in which the soil resistivity, towers and conductor heights, conductors’ transpositions, towers geometry, and spacing between consecutive towers are assumed as constants. However, these assumptions introduce uncertainties and errors to the actual line’s parameters in practical cases [36].Errors in the parameters of the protected line do not affect the proposed method’s calculation of resistances and inductances. However, these errors affect the region comprised by the zone since it depends on such parameters, and therefore may affect all impedance-based distance protection solutions.In this way, depending on the errors’ magnitude, the resistance and inductance values estimated by the algorithm may enter the zone a few sampling periods before or after the instant that they would in the ideal situation without uncertainties. Thus, the line parameters’ deviations may affect the fault probability.Therefore, to verify the effect of the parameters’ uncertainties and errors on the proposed method, the authors deliberately introduced errors/uncertainty levels of± 5%,± 10%, and ± 15% in each simulated line’s zero and positive-sequence impedances and reperformed all simulations. Table 8, then, presents the average, maximum, and minimum deviations in the trip times, considering all six impedance elements and all cases simulated, caused by each parameter uncertainty level, both for the proposed algorithm (LS + BI) as for the conventional phasor-based distance protection method (DFT).Table 8
Deviations on trip times [ms]—proposed (LS + BI) and conventional (DFT).
Uncertainty (%)Avg LS + BIMax LS + BIMin LS + BIAvg DFTMax DFTMin DFT+5%−0.09830−0.5201−0.13700−1.563−5%0.11530.520100.11661.5630+10%−0.18000−0.5201−0.24510−1.56310%0.18980.520100.22421.5630+15%−0.25230−0.5201−0.34000−1.563−15%0.21960.520100.30091.5630For the proposed method, the most significant deviations between trip times were of one sampling period,± 0.521 [ms]. For the conventional phasor-based approach, in turn, the most significant deviations between trip times were of three sampling periods, ± 1.563 [ms]. In addition, the average deviations were smaller for the proposed method, for all uncertainty levels. Furthermore, in all cases, with all considered errors, the presented algorithm still proved efficient in rapidly detecting all simulated faults. Therefore, the proposed solution is robust against the protected line’s parameters’ uncertainties. In addition, it is noteworthy that an increase in the line’s parameters may decrease trip times and vice versa. This happens because deviations in the parameters provoke changes in the zone boundaries, as mentioned earlier.Figure8, finally, displays the fault probabilities for the same fault case under the cases with no uncertainties and ± 15% uncertainties in the parameters, in which the error of + 15% in the line’s impedances decreased the trip time with no uncertainties by one sampling period while the error of −15% increased the trip time with no uncertainties by one sampling period.Figure 8
Fault probabilities and trip instants for AB element in the same fault case with 0 and± 15% uncertainties in line 2 impedance parameters.


### 4.7. Power Swings

To evaluate the proposed algorithm in a situation with power swings, usually challenging to distance protection solutions [4, 34, 35], the authors simulated the system illustrated in Figure 9, based on [37]. Protection was applied on bus 2 to protect both circuits of the line connecting buses 2 and 3. On the upper circuit (number 2), a three-phase fault located at 40 [km] from bus 2 starts at 0.453 [s] and is cleared in 0.5 [s]. After the fault clearance, the generator’s load angle δ begins to vary, causing the lower circuit (number 1) voltages and currents to oscillate. At 4 [s], the fault occurs in circuit 1.Figure 9
Power swing electrical system.The authors simulated three-phase symmetrical faults as these faults are usually the most challenging to be detected by the protection relays during power swings, precisely because they are balanced [34, 35].The line’s positive and zero-sequence impedance (at 50 [Hz]) are, respectively, equal to (0.07+j0.42) and (0.21+j0.94) [Ω/km]. The zero-sequence reactance between circuits is j0.43 [Ω/km]. Transformer and generator parameters are presented in [37]. The authors resampled the signals to 4 [kHz].The trip times of the LS + BI and DFT algorithms in circuit 2 were equal to 4.50 and 15.00 [ms], while in circuit 1, the trip times were 4.75 and 15.75 [ms]. Figure10 presents the CG element pℱ|Γ probability in circuit 2, BC element pℱ|Γ probability in circuit 1, and angle δ over time. In this simulated power swing critical scenario, the proposed algorithm did not present any protection failures.Figure 10
Fault probabilities and load angle.


## 4.1. Noisy Measurements

The authors applied Gaussian noises with an average signal-to-noise ratio (SNR) of 15 [dB] to both voltage and current signals from lines 2 and 5 to verify noisy measurements’ effect on the algorithms’ responses. Figure3 presents the normal and noisy current signals for a three-phase fault online 2, and Figure 4 illustrates the probabilities pℱ|Γ for the CG element in both cases highlighting the methods’ trip instants. The average trip times in [ms] obtained in the L terminal of the lines in the noisy situation are presented in Table 4.Figure 3
Current signal with and without 15 [dB] Gaussian noise.Figure 4
Fault probabilities with and without noise over time.Table 4
Average trip timesL [ms]—Gaussian noise.
IdLS + BI-LGDFT-LGLS + BI-LLDFT-LL29.2014.297.9014.53510.3214.578.4114.38


## 4.2. CT Saturation and CCVT Transients

The authors modeled at the local terminal of line 2 a CCVT with the parameters obtained in [32] through an actual 230 [kV] equipment and a CT with 1200 : 5 transformation ratio and magnetisation characteristic based on [33] to verify the effects of CT saturation and CCVT transients on the trip times.Figure5 illustrates the saturated and normal currents for a three-phase fault, while Figure 6 presents the voltage with and without the CCVT for the same case. Figure 7, in turn, presents the AB element’s fault probability pℱ|Γ for this case with and without the CT and CCVT effects. Finally, Table 5 presents the average trip times with (1) and without (0) the instrumentation transformers.Figure 5
Normal and saturated current signal.Figure 6
Voltage signal with and without CCVT transients.Figure 7
Fault probabilities with and without the CT + CCVT effects.Table 5
Average trip timesL [ms]—instrumentation transformers.
CCVTCTLS + BI-LGDFT-LGLS + BI-LLDFT-LL107.1414.387.8114.55018.0525.8312.7827.22118.7626.2712.8727.50


## 4.3. Close-In Faults

For the algorithm’s evaluation in close-in fault events, the authors simulated the faults in line 2 at 2 [km] from the local terminal, which corresponds to 2.5 % of the line’s total length. Conditions with and without the instrumentation transformers were considered. Table6 presents the average trip times.Table 6
Average trip timesL [ms]—close-in faults.
CCVTCTLS + BI-LGDFT-LGLS + BI-LLDFT-LL004.749.254.698.85115.7814.204.7914.79


## 4.4. Line Energizing onto Faults

As mentioned in Section2.5, the initial voltage and current values, stored in the IED’s memory before line energizing, can be collected during standard line operation for future use or synthesized from the nominal loading conditions.The initial voltage and current stored values may influence the trip times, especially in switch onto fault conditions. Thus, to test their influence in these conditions, the authors simulated faults at line 2 energizing (by terminalL) and tested the algorithm with the following sets of maxNLG,NLL  = 9 (considering the default settings) initial values sets (1–3).(i)
1: collected in past standard operation(ii)
2: synthetic light loading values(iii)
3: synthetic heavy loading valuesTable7 presents the average trip times in [ms] for the three data sets. The conventional DFT-based algorithm does not need initial values. However, the results show that the proposed algorithm is faster with all considered data sets.Table 7
Average trip timesL [ms]—switch onto faults.
Data setLS + BI-LGDFT-LGLS + BI-LLDFT-LL110.3017.9510.3720.63210.4417.9510.3520.63310.9517.959.7620.63


## 4.5. System Frequency Variations

To evaluate the performance of the algorithms under operating conditions with frequencies other than the fundamental, which may occur due to variations in the system’s generation and load levels, the authors simulated the faults in line 2 considering the frequency range of 58–62 [Hz].The authors did not observe any variations in the trip times of the proposed algorithm and verified variations up to 1 [ms] in the DFT-conventional algorithm trip times. These variations happen because the algorithm based on phasors calculates reactance values at the fundamental frequency.


## 4.6. Uncertainties and Errors in Line Parameters

Generally, the protected line’s impedance values are determined through calculations in which the soil resistivity, towers and conductor heights, conductors’ transpositions, towers geometry, and spacing between consecutive towers are assumed as constants. However, these assumptions introduce uncertainties and errors to the actual line’s parameters in practical cases [36].Errors in the parameters of the protected line do not affect the proposed method’s calculation of resistances and inductances. However, these errors affect the region comprised by the zone since it depends on such parameters, and therefore may affect all impedance-based distance protection solutions.In this way, depending on the errors’ magnitude, the resistance and inductance values estimated by the algorithm may enter the zone a few sampling periods before or after the instant that they would in the ideal situation without uncertainties. Thus, the line parameters’ deviations may affect the fault probability.Therefore, to verify the effect of the parameters’ uncertainties and errors on the proposed method, the authors deliberately introduced errors/uncertainty levels of± 5%,± 10%, and ± 15% in each simulated line’s zero and positive-sequence impedances and reperformed all simulations. Table 8, then, presents the average, maximum, and minimum deviations in the trip times, considering all six impedance elements and all cases simulated, caused by each parameter uncertainty level, both for the proposed algorithm (LS + BI) as for the conventional phasor-based distance protection method (DFT).Table 8
Deviations on trip times [ms]—proposed (LS + BI) and conventional (DFT).
Uncertainty (%)Avg LS + BIMax LS + BIMin LS + BIAvg DFTMax DFTMin DFT+5%−0.09830−0.5201−0.13700−1.563−5%0.11530.520100.11661.5630+10%−0.18000−0.5201−0.24510−1.56310%0.18980.520100.22421.5630+15%−0.25230−0.5201−0.34000−1.563−15%0.21960.520100.30091.5630For the proposed method, the most significant deviations between trip times were of one sampling period,± 0.521 [ms]. For the conventional phasor-based approach, in turn, the most significant deviations between trip times were of three sampling periods, ± 1.563 [ms]. In addition, the average deviations were smaller for the proposed method, for all uncertainty levels. Furthermore, in all cases, with all considered errors, the presented algorithm still proved efficient in rapidly detecting all simulated faults. Therefore, the proposed solution is robust against the protected line’s parameters’ uncertainties. In addition, it is noteworthy that an increase in the line’s parameters may decrease trip times and vice versa. This happens because deviations in the parameters provoke changes in the zone boundaries, as mentioned earlier.Figure8, finally, displays the fault probabilities for the same fault case under the cases with no uncertainties and ± 15% uncertainties in the parameters, in which the error of + 15% in the line’s impedances decreased the trip time with no uncertainties by one sampling period while the error of −15% increased the trip time with no uncertainties by one sampling period.Figure 8
Fault probabilities and trip instants for AB element in the same fault case with 0 and± 15% uncertainties in line 2 impedance parameters.


## 4.7. Power Swings

To evaluate the proposed algorithm in a situation with power swings, usually challenging to distance protection solutions [4, 34, 35], the authors simulated the system illustrated in Figure 9, based on [37]. Protection was applied on bus 2 to protect both circuits of the line connecting buses 2 and 3. On the upper circuit (number 2), a three-phase fault located at 40 [km] from bus 2 starts at 0.453 [s] and is cleared in 0.5 [s]. After the fault clearance, the generator’s load angle δ begins to vary, causing the lower circuit (number 1) voltages and currents to oscillate. At 4 [s], the fault occurs in circuit 1.Figure 9
Power swing electrical system.The authors simulated three-phase symmetrical faults as these faults are usually the most challenging to be detected by the protection relays during power swings, precisely because they are balanced [34, 35].The line’s positive and zero-sequence impedance (at 50 [Hz]) are, respectively, equal to (0.07+j0.42) and (0.21+j0.94) [Ω/km]. The zero-sequence reactance between circuits is j0.43 [Ω/km]. Transformer and generator parameters are presented in [37]. The authors resampled the signals to 4 [kHz].The trip times of the LS + BI and DFT algorithms in circuit 2 were equal to 4.50 and 15.00 [ms], while in circuit 1, the trip times were 4.75 and 15.75 [ms]. Figure10 presents the CG element pℱ|Γ probability in circuit 2, BC element pℱ|Γ probability in circuit 1, and angle δ over time. In this simulated power swing critical scenario, the proposed algorithm did not present any protection failures.Figure 10
Fault probabilities and load angle.


## 5. Actual Fault Records

The authors tested the presented algorithm with voltage and current oscillographs from seventy-nine faults recorded by the actual IEDs of four existing 500 [kV] and 60 [Hz] lines. The faults occurred between 2017 and 2020. Each IED has a distance function as primary protection and a directional overcurrent function as backup [38]. The authors collected the line’s impedance from the respective field IEDs to test the presented algorithm against the actual fault records, presented in Table 9.Table 9
Actual lines characteristics.
IdTypeL (km)Z1 [Ω/km]Z0 [Ω/km]6Single366.370.017+j0.2660.355+j1.2607Single290.930.016+j0.2620.322+j1.1628Single321.720.018+j0.2730.438+j1.5489Single248.570.015+j0.2630.364+j1.063The IEDs from lines 6,7, and 8 have a 2.0 [kHz] sampling rate, while the IED from line 9 has a 3.9 [kHz] sampling rate. Furthermore, it is noteworthy that lines 6, 7, and 9 present series compensation.The average trip times of the LS + BI algorithm (all elements), DFT-based (all elements), and real IEDs are presented in Table10. The trip logic and settings of the real IEDs are different from the used in the tested algorithms. Thus, the inclusion of the actual IEDs trip times has only an illustrative purpose.Table 10
Average trip times [ms]—actual faults.
IdNumber of faultsLS + BIDFTReal IED61912.8621.1724.307209.3819.4233.638208.0415.7631.769205.4810.9934.40Figure11 presents the current of a C-G fault that occurred online 6 at 12 : 32 pm (GMT-3) on September 24, 2018. This fault was caused by fire and was located at 200.9 [km] from the relay bus. For this fault, the trip time of the LS + BI method, DFT method, and the real IED was, respectively, equal to 5.88, 13.24, and 21.08 [ms]. Figure 12 presents the CG element’s probability pℱ|Γ over time.Figure 11
Phase C faulted phase current.Figure 12
CG element’s fault probability and trip instants.Figure13, in turn, illustrates the current of a A-G fault that occurred online 8 at 5 : 27 am (GMT-3) on February 8, 2020, caused by a lightning bolt that struck at 71.2 [km] from the relay bus. The prefault current has a total harmonic distortion (THD) of 17.71 %. For this fault, the trip time of the LS + BI method, DFT method, and the real IED was, respectively, equal to 5.39, 13.73, and 24.51 [ms]. Figure 14 presents the AG element’s probability pℱ|Γ over time.Figure 13
Phase A faulted phase current.Figure 14
AG element’s fault probability and trip instants.


## 6. Conclusion

This paper presented a distance protection algorithm based on the combination of estimating the resistances and inductances seen by the relay through LS formulations derived from the line differential equations in discrete time with applying BI to determine the in-zone fault probabilities over time. The probabilities are the basis for the method to trigger the trip outputs. The proposed algorithm differs from previous research as it introduces an application of BI to calculate, based on the estimates, the fault within the zone probability. The use of Bayesian probability increases the reliability of the classic LS method, as it supports a trip criterion that proved to be safe and reliable in all tested conditions.The authors tested the proposed method with fault simulations in ATP under different conditions (including noisy measurements, CT saturation, CCVT transients, close-in faults, line energizing onto faults, frequency variations, uncertainties and errors in the line’s parameters, and power swings) and real faults in actual transmission lines considering a single set of default settings, which proved to be reliable and with applicability in diverse electrical systems and situations. The obtained results indicate that the proposed distance protection solution presented is fast (capable of detecting faults in less than half a cycle in many cases), safe, and robust, configuring it as an interesting and potential approach for protecting transmission lines with the distance function.The use of BI to calculate fault probabilities increases the reliability of the LS-estimated resistances and inductances since it presents faster fault detection times, especially in faults near to the first zone boundary, and more rapid trip magnitude stabilization times. Therefore, the proposed method improves the methods based only on differential equations and LS-model fitting, which are not as used as DFT-based solutions in modern digital IEDs due to the numerical inaccuracies related to LS solving. It is also noteworthy that the proposed algorithm does not require any improvement in conventional protection infrastructures and high sampling frequencies.



---

*Source: 4912864-2022-06-01.xml*