Journal of Modern Power Systems and Clean Energy

ISSN 2196-5625 CN 32-1884/TK

网刊加载中。。。

使用Chrome浏览器效果最佳,继续浏览,你可能不会看到最佳的展示效果,

确定继续浏览么?

复制成功,请在其他浏览器进行阅读

Improved Leap-frog Method for Time-domain Fault Location  PDF

  • Džafić Izudin 1 (Senior Member, IEEE)
  • A. Jabr Rabih 2 (Fellow, IEEE)
1. Faculty of Electrical Engineering, University of Sarajevo, Sarajevo 71000, Bosnia and Herzegovina ; 2. Department of Electrical and Computer Engineering,American University of Beirut, Beirut 1107 2020, Lebanon

Updated:2024-03-26

DOI:10.35833/MPCE.2023.000175

  • Full Text
  • Figs & Tabs
  • References
  • Authors
  • About
CITE
OUTLINE

Abstract

The partial differential equation (PDE) solution of the telegrapher is a promising fault location method among time-domain and model-based techniques. Recent research works have shown that the leap-frog process is superior to other explicit methods for the PDE solution. However, its implementation is challenged by determining the initial conditions in time and the boundary conditions in space. This letter proposes two implicit solution methods for determining the initial conditions and an analytical way to obtain the boundary conditions founded on the signal decomposition. The results show that the proposal gives fault location accuracy superior to the existing leap-frog scheme, particularly in the presence of harmonics.

I. Introduction

FAULT location is an essential function in the operation of modern power systems for lessening the outage time following a fault event [

1], [2]. In addition, the accurate location of the fault contributes to speeding up the process of dispatching a specialized crew to attend to the fault restoration [3]. For long transmission lines, the time-domain and model-based techniques employ distributed parameters that necessitate a partial differential equation (PDE) solution for locating the fault, starting from measurements at both ends of the faulted line [4]-[6]. The digital computer solution is a finite-difference method in which discrete space-time replaces continuous space-time. In discrete space-time, the finite-difference equations substitute PDEs, giving rise to explicit solution methods where space-time stepping solution schemes can be implemented without solving large matrix equations. The explicit solution methods differ by how the partial derivatives are approximated in the discrete space-time and include the techniques such as the standard first-order difference form, the Lax-Wendreoff scheme, and the leap-frog (LF) method. Reference [4] has recently shown that the LF method is the best among other explicit techniques for fault location applications. The Courant-Friedrichs-Lewy criterion [7] governs its numerical stability.

This letter proposes improvements to the LF method applied to the fault location. The first one is calculating initial values along the x-axis (line length) using extracted phasors, which provides significant improvements, particularly when the system is contaminated with harmonics. The issue of harmonics has become critical due to the widespread use of power electronics. The second improvement is correctly estimating the initial voltage values for the first space-time stepping solution of LF. In particular, two new implicit-based schemes are introduced for calculating these initial values.

II. LF Method

The equations of telegrapher are a pair of coupled PDEs that govern the voltage u(x,t) and current i(x,t) on the transmission line in the function of distance x and time t:

-i(x,t)x=u(x,t)G'+C'u(x,t)t-u(x,t)x=R'i(x,t)+L'i(x,t)t (1)

The coefficients of the PDEs are in terms of distributed parameters: R', L', C', and G' are the the series resistance, series inductance, shunt capacitance, and shunt conductance (per unit length), respectively. By differentiating the set of PDEs in (1) with respect to t and x and eliminating the terms involving i(x,t), a second-order hyperbolic PDE can be derived as:

2u(x,t)x2-L'C'2u(x,t)t2=(R'C'+G'L')u(x,t)t+G'R'u(x,t) (2)

The LF method operates on (2) and employs second-order central differences for approximating the space and time derivatives. Reference [

4] reports that the voltage solution on the discrete space-time grid (as shown in Fig. 1) is obtained from the space-time stepping formula:

ux+1,t=-ux-1,t+(2-2v2b1+Δx2b3)ux,t+v2b1-vΔt2b2ux,t-1+v2b1+vΔt2b2ux,t+1 (3)

Fig. 1  Required data, calculation progress, and calculation space for LF method.

where b1=L'C'; b2=R'C'+L'G'; b3=R'G'; v=Δx/Δt; and the voltage subscripts x and t (without parenthesis) are integer values representing positions on the power grid with spacings Δx and Δt. As shown in (3), the LF calculation scheme requires the knowledge of values for x=0 and x=1 for all t. Therefore, the calculation in (3) begins with x=1 and t=1 and computes all values along the t-axis. After that, x is increased by one, and the calculation is repeated for all possible values along the t-axis, as shown in Fig. 1, where nodes enclosed with red ellipses are used to calculate the current calculated point, and the arrows indicate calculation steps.

The two-ended fault location technique solves PDEs twice using the LF method. The first solution starts at the sending end and the second at the receiving end. The fault location can be determined by comparing the voltage profiles computed from the sending and receiving ends; and the mismatch between these voltages is the least at the fault location. In the computation, the solution is found in the αβ0 domain, making the goal function (whose components are shown in Fig. 2) as (4). The voltage vector at point n (on x-axis) from the sending end is denoted as us,n(i), while the voltage vector at point N-n (on x-axis) from the receiving end is denoted as ur,N-n(i). Both vectors are computed for all available time instants while the fault occurs. represents the 𝓁1 vector norm.

f(n)=i{α,β,0}us,n(i)-ur,N-n(i)    n=0,1,,N (4)

Fig. 2  Goal functions |us-ur| in α, β, and 0 components for BCG fault (phases B and C to ground) on 240 km line, fault at x=85 km.

A. Initial Values Along x-axis

Although the initial voltage value for t=0 and x=0, i.e., u(0,0), is given by the intelligent electronic device (IED) at the sending end of the line, the values u(x,0) along the x-axis for t=0 are unavailable. Therefore, to compute u(x,0), the LF implementation in [

4] uses linear interpolated values of samples provided by IEDs on the line sending- and receiving-end terminals. Unfortunately, this method can give poor results whenever the voltage or current peak is at the midspan of the line or in the presence of significant harmonics.

Assuming that the system was in a steady state before the fault, one can consider the presence of a fundamental component voltage signal and its harmonics. This letter proposes an analytical procedure for calculating the pre-fault voltage and current profiles along the line length (x-axis) starting from the phasor estimation of the fundamental component and harmonics described in [

8]. The method in [8] employs orthogonal component (OC) phasor estimation with a Prony variation on sampled measured value (SMV) data to extract the initial (t=0) phasor values of the voltage and current at the sending end. The calculation process uses OC to compute the phasor of the fundamental component. Once obtained, the fundamental component is removed from the original signal, and the Prony method calculates the initial phasor values of harmonics. The initial values are then used as the inputs to the nonlinear least squares (NLS) technique [8], which has much better behavior in the presence of noise. Using the extracted phasors (u̲k and i̲k), the voltage values along the x-axis can be quickly obtained by (5) and transformed to the time-domain:

u̲x=k=1nh(u̲kcoshγkx-zki̲ksinhγkx) (5)

where zk=R'+ikωL'G'+ikωC', and ω is the fundamental system frequency; γk=(R'+ikωL')(G'+ikωC'); and nh is the number of harmonics.

The use of the wavelet transformation for double-ended fault location has been discussed in phasor extraction [

9]. Wavelets have been used in traveling-wave methods [10] to improve the performance when dealing with high harmonic content. The wavelet transform isolates the traveling waves from other signals such as harmonics on transmission lines.

B. Initial Values Along t-axis

Providing initial values for the LF method requires a solution of PDEs of the telegrapher (1) for a single step in the discrete space-time representation. The process is solved implicitly for x=1 (scheme 1 in Section II-B-1)) or x=-1 (scheme 2 in Section II-B-2)). The general algorithmic framework is as follows for sending-end calculations, and similar computations apply to the receiving end.

Step 1:   initialization. Collect the sampled data and convert the sampled voltage data, current data, and the line parameters to the αβ0 domain.

Step 2:   implicit formulation. Formulate the diagonally dominant sparse matrix A and the right-hand-side vector b of the PDE discrete space-time representation using (6), (7) for scheme 1 and (8), (9) for scheme 2.

Step 3:   iterative solution. Using an iterative solver, solve Ax=b, where x=[u1,1,u1,2,,u1,n,i1,1,i1,2,,i1,n]T for scheme 1 and x=[u-1,1,u-1,2,,u-1,n,i1,1,i1,2,,i1,n]T for scheme 2.

Step 4:   LF method. Apply the LF method using the measured voltage initial values umt=[u0,1,u0,2,,u0,n]T, the calculated initial values along the x-axis uct=[u1,0,u2,0,,um,0]T for scheme 1 and uct=[u-1,0,u1,0,,um,0]T for scheme 2, and the calculated initial voltage values uct=[u1,1,u1,2,,u1,n]T for scheme 1 and uct=[u-1,1,u-1,2,,u-1,n]T for scheme 2.

1) Scheme 1 (Forward Calculation)

To provide initial values required by the LF method, this scheme uses an implicit solution employing a centered derivative approximation with second-order accuracy along the t-axis and a forward approximation of first-order along the x-axis (represented by the green ellipses in Fig. 3) [

11]:

Fig. 3  Calculating initial values using unmodified space-time grid.

u(x,t)ti,j=u(xi,tj+1)-u(xi,tj-1)2Δt=u1,j+1-u1,j-12Δti(x,t)ti,j=i(xi,tj+1)-i(xi,tj-1)2Δt=i1,j+1-i1,j-12Δt (6)
u(x,t)xm,j=u(xm,tj)-u(xm-1,tj)Δx=u1,j-u0,jΔxi(x,t)xm,j=i(xm,tj)-i(xm-1,tj)Δx=i1,j-i0-1,jΔx (7)

The solution to Ax=b in Step 3 provides all values for x=1 (represented by the green dots in Fig. 3) after which the LF method can be started in Step 4. The LF method begins with the points enclosed in the red ellipses in Fig. 3 and calculates the red colored grid values.

2) Scheme 2 (Backward Calculation)

Scheme 2 is similar to the first but differs from it by using a backward approximation of first-order along the x-axis instead of a forward approximation [

11]:

u(x,t)ti,j=u(xi,tj+1)-u(xi,tj-1)2Δt=u-1,j+1-u-1,j-12Δti(x,t)ti,j=i(xi,tj+1)-i(xi,tj-1)2Δt=i-1,j+1-i-1,j-12Δt (8)
u(x,t)xm,j=u(xm,tj)-u(xm-1,tj)Δx=u0,j-u-1,jΔxi(x,t)xm,j=i(xm,tj)-i(xm-1,tj)Δx=i0,j-i-1,jΔx (9)

The solution to Ax=b provides all values for x=-1 (represented by the green dots in Fig. 4). As shown in scheme 1, the LF method begins with the points enclosed in the red ellipses in Fig. 4 and calculates the red colored grid values. Scheme 2 can be only applied with phasor extraction since it requires the values u-1,0 and i-1,0.

Fig. 4  Calculating initial values by extending space-time grid.

III. Numerical Results

The two-terminal fault location method is implemented by comparing the αβ0 components calculated using the samples from the sending- and receiving-end terminals. The calculation is done by solving the PDEs in the forward and backward directions, i.e., starting from the sending and receiving ends, using the LF method. The numerical methods were implemented in MATLAB using in-house developed solvers. Three initialization schemes were investigated, considering feasible combinations of two possibilities for the boundary condition (ux,0) and two for the initial voltage values required to initiate the LF computation.

1) Int+ refers to linear interpolation for the initial voltage profile ux,0 and estimation of the u1,t voltages using scheme 1. Int+ is closest to the initialization method in [

4]. Int- differs from Int+ by estimating the u-1,t voltages using scheme 2, so it is infeasible using linear interpolation. The reason is that the backward calculation requires the values u-1,0 and i-1,0, but the interpolation method cannot give these values.

2) Extr+ refers to phasor extraction followed by the analytical voltage solution (described in Section II-A) for the initial voltage profile ux,0 and estimation of the u1,t voltages using scheme 1. Extr- differs from Extr+ by estimating the u-1,t voltages using scheme 2.

The tests were carried out on a transmission line (as shown in Fig. 4) that is 240 km long and whose positive- and zero-sequence resistance, inductance, capacitance, and conductance are as follows: R1'=0.006365 Ω/km, R0'=0.01932 Ω/km, L1'=0.09337 mH/km, L0'=0.41264 mH/km, C1'=127.4 nF/km, C0'=127.4 nF/km, G1'=0 S/km, and G0'=0 S/km. The fault remains for 5 ms and has a resistance Rfalult=1 Ω. The sampling rate of IED is fs=20 kHz. The test setup of fault location is shown in Fig. 5.

Fig. 5  Test setup of fault location.

Table I shows the fault location error for single line-to-ground (AG) fault without harmonics. The first column shows the three actual fault locations (xf), the second column gives the fault inception angle (θs), and the last one gives the improvement factor (IF), i.e., the ratio of the worst error (Int+) relative to the best (Extr-). Table II is the counterpart of Table I, where harmonics giving a total harmonic distortion (THD) of 5% are simulated over 100 Monte-Carlo trials. The improvement, in this case, is quite remarkable and gives an average percentage error reaching around eight times lower than the standard Int+ initialization.

Table I  Fault Location Error for Single Line-to-ground (AG) Fault Without Harmonics
xf(km)θs (°)Error (%)IF
Int+Int-Extr+Extr-
20 0 0.716 × 0.688 0.627 1.141
45 1.147 × 1.064 1.008 1.138
90 1.322 × 1.193 1.163 1.137
65 0 1.163 × 1.086 1.034 1.124
45 1.205 × 1.107 1.079 1.117
90 1.970 × 1.876 1.789 1.101
120 0 0.683 × 0.610 0.607 1.126
45 0.831 × 0.764 0.731 1.137
90 0.998 × 0.951 0.899 1.110

Note:   × indicates that Int- cannot be applied.

Table II  Fault Location Error for AG Fault with Harmonics (5% THD Harmonics, 100 Monte-Carlo Trials)
xf (km)θs (°)Average error (%)IF
Int+Int-Extr+Extr-
20 0 15.075 × 3.765 3.033 4.971
45 17.089 × 4.325 4.182 4.086
90 20.711 × 5.802 4.959 4.177
65 0 13.801 × 3.774 2.965 4.655
45 15.943 × 3.765 3.127 5.098
90 17.805 × 4.097 3.436 5.182
120 0 11.066 × 3.876 2.139 5.173
45 13.721 × 3.096 1.996 6.875
90 15.097 × 2.973 1.890 7.987

Note:   × indicates that Int- cannot be applied.

Table III presents fault location error comparison between LF Extr- and the method of characteristics (MC) [

12]. The MC is a classical benchmark method that solves partial differential equations for model-based fault location. The results are based on different fault locations (xf) and fault inception angles (θs) with 5% THD harmonics. The LF Extr- method consistently yields more accurate results with less error (shown in both km and percentage). It also leads to an IF of up to 3.2.

Table III  Fault Location Error Comparison Between LF Extr- and MC for AG Fault with Harmonics (5% THD Harmonics, 100 Monte-Carlo Simulations)
xf (km)θs (°)Estimated (km)Error (km)Error (%)IF
MCLFMCLFMCLF
20 0 24.812 22.147 4.812 2.147 2.005 0.895 2.241
45 25.713 22.618 5.713 2.618 2.380 1.091 2.182
90 26.440 21.976 6.440 1.976 2.683 0.823 3.259
135 25.913 22.633 5.913 2.633 2.464 1.097 2.246
180 26.728 17.914 6.728 2.086 2.803 0.869 3.225
65 0 59.276 62.215 5.724 2.785 2.385 1.160 2.055
45 58.301 62.416 6.699 2.584 2.791 1.077 2.592
90 61.855 63.971 3.145 1.029 1.310 0.429 3.056
135 71.029 68.102 6.029 3.102 2.512 1.293 1.944
180 70.312 66.993 5.312 1.993 2.213 0.830 2.665
120 0 125.013 122.711 5.013 2.711 2.089 1.130 1.849
45 124.997 122.323 4.997 2.323 2.082 0.968 2.151
90 126.713 123.192 6.713 3.192 2.797 1.330 2.103
135 114.992 122.175 5.008 2.175 2.087 0.906 2.303
180 124.947 122.611 4.947 2.611 2.061 1.088 1.895

Table IV presents the computational time for LF fault location of different fault types, based on 10 Monte-Carlo trials. The results show that the calculation requires around 448 ms on a standard personal computer for an AG fault. Table IV also reports the average computational time using the LF method for the line-to-line (BC) fault, the double line-to-ground (BCG) fault, and the three-phase (ABC) fault. The simulations were conducted on a MacBook Pro, M1 max, 32 GB RAM, macOS Ventura 13.5. It is important to note that the computational time of the fault location method is not a valuable metric for comparing fault location methods, as the fault location application is not a real-time function [

8]. The most crucial factor is ensuring that the fault location is accurately determined. In practice, fault location starts once the fault occurrs, requiring a crew to be dispatched to fix the issue. It may take several hours for a utility crew to fix a problem on-site, mainly if it occurs in a remote mountain area without road access. The time required to fix a fault increases when the fault location is not accurately known. Therefore, pinpointing the fault location is crucial for prompt repair. The time it takes to locate faults is negligible compared with that for dispatching a team and repairing the faults.

Table IV  Computational Time for LF Fault Location of Different Fault Types Based on 10 Monte-Carlo Trials
Fault typeInvolved componentAverage computational time (ms)
αβ0
AG × 448.3
BC × × 262.1
BCG 603.7
ABC × 452.8

Note:   and × indicate that the component is involved and not involved, respectively.

IV. Conclusion

Among time-domain and model-based fault location techniques, the PDE solution of the telegrapher is the most effective. This letter builds upon recent research [

4], demonstrating that the LF process is superior to other explicit PDE solutions. Nevertheless, the LF implementation faces challenges in determining the initial conditions in time and the boundary conditions in space. The letter has presented a promising strategy for initializing the LF solution to PDEs involved in time-based fault location, thereby promoting its adoption in operational settings. In addition, when harmonics are present, the proposed initialization significantly improves the accuracy of fault location compared with [4] and [12].

Acknowledgements

The work of Izudin Džafić was supported by the Federal Ministry of Education and Science, Bosnia, through funding.

References

1

J. Hu, W. Hu, J. Chen et al., “Fault location and classification for distribution systems based on deep graph learning methods,” Journal of Modern Power Systems and Clean Energy, vol. 11, no. 1, pp. 35-51, Jan. 2023. [Baidu Scholar] 

2

Y. Zhu and H. Peng, “Multiple random forests based intelligent location of single-phase grounding fault in power lines of DFIG-based wind farm,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 5, pp. 1152-1163, Sept. 2022. [Baidu Scholar] 

3

R. A. Jabr and I. Džafić, “Distribution management systems for smart grid: architecture, work flows, and interoperability,” Journal of Modern Power Systems and Clean Energy, vol. 10, no. 2, pp. 300-308, Mar. 2022. [Baidu Scholar] 

4

D. Lu, Y. Liu, Q. Liao et al., “Time domain transmission line fault location method with full consideration of distributed parameters and line asymmetry,” IEEE Transactions on Power Delivery, vol. 35, no. 6, pp. 2651-2662, Dec. 2020. [Baidu Scholar] 

5

J. Suonan, G. Song, Q. Xu et al., “Time-domain fault location algorithm for parallel transmission lines using unsynchronized currents,” International Journal of Electrical Power & Energy Systems, vol. 28, no. 4, pp. 253-260, May 2006. [Baidu Scholar] 

6

A. Bendjabeur, A. Kouadri, and S. Mekhilef, “Transmission line fault location by solving line differential equations,” Electric Power Systems Research, vol. 192, p. 106912, Oct. 2021. [Baidu Scholar] 

7

R. Courant, K. Friedrichs, and H. Lewy, “On the partial difference equations of mathematical physics,” Mathematische Annalen, vol. 100, pp. 32-74, Jan. 1928. [Baidu Scholar] 

8

I. Džafić, R. A. Jabr, and N. Fetić, “Kung’s component extraction in power system fault location,” International Journal of Electrical Power & Energy Systems, vol. 119, p. 105888, Jul. 2020. [Baidu Scholar] 

9

T. Łobos, P. Kostyla, J. Pospieszna et al., “Location of faults on transmission lines using wavelet transforms,” in Proceedings of 2008 International Conference on High Voltage Engineering and Application, Chongqing, China, Nov. 2008, pp. 633-636. [Baidu Scholar] 

10

W. Fluty and Y. Liao, “Electric transmission fault location techniques using traveling wave method and discrete wavelet transform,” in Proceedings of 2020 Clemson University Power Systems Conference (PSC), Clemson, USA, Mar. 2020, pp. 1-8. [Baidu Scholar] 

11

B. Fornberg, “Generation of finite difference formulas on arbitrarily spaced grids,” Mathematics of Computation, vol. 51, no. 184, pp. 699-706, Oct. 1988. [Baidu Scholar] 

12

A. Gopalakrishnan, M. Kezunovic, S. McKenna et al., “Fault location using the distributed parameter transmission line model,” IEEE Transactions on Power Delivery, vol. 15, no. 4, pp. 1169-1174, Oct. 2000. [Baidu Scholar]