A time domain differential equation approach using distributed parameter line model for transmission line fault location algorithm

Electric Power Systems Research 46 (1998) 1 – 10

A time domain differential equation approach using distributed parameter line model for transmission line fault location algorithm Bao Lian a, M.M.A. Salama a, A.Y. Chikhani b,* a

b

Electrical and Computer Engineering Department, Uni6ersity of Waterloo, Waterloo, Ontario, Canada Department of Electrical and Computing Engineering, Royal Military College of Canada, Kingston, Ontario, Canada Received 12 February 1997; accepted 18 November 1997

Abstract The objective of this work is to develop an algorithm which can provide both a fast relaying time (1/4–1/2 cycles) and an accurate fault location (error within 1% of the total line length). The proposed algorithm combines the strengths of the two distinct approaches, namely, the steady-state phasor approaches (SSPA) and the differential equation approaches (DEA). In the proposed algorithm, a distributed line model is used, which is suitable for the short, the long, the transposed and the untransposed lines. The DEA is used to calculate the fault location and is implemented using the Laplace and the Z-transform techniques. There are several advantages in this proposed algorithm over the existing ones: (a) the variations in the line parameters with the frequency are considered; (b) the fault location estimate is less sensitive to the selection of the sampling rate, the fault inception angle, the signal spectrum and the equivalent sources behind both bus bars; and (c) the relay performance is greatly improved offering both the desired high accuracy and the speed at the same time. © 1998 Published by Elsevier Science S.A. All rights reserved. Keywords: Algorithm; Steady-state phasor approach; Differential equation approach

1. Introduction The steady-state phasor approaches (SSPA) and the differential equation approaches (DEA) have been used widely in the fault location algorithm. In the SSPA, the sequence or the modal networks are used and the fault distance is obtained using the post-fault steady-state voltage and current phasors. A high accuracy can be obtained if a sufficient filtering time is allowed for the extraction of the phasors. In a one-end lumped line model approach, line transposition is usually assumed and shunt capacitance is neglected. The source impedance and the load at the remote end need to be estimated and the effect of the fault impedance and ground conditions has to be reduced. These aspects are treated differently in each algorithm as shown in Refs. [1 – 4]. A two-end lumped model approach does not need the source impedance and the fault resistance can easily be eliminated [5,6]. The accuracy can be further improved with a distributed line model since the shunt

capacitance is automatically considered; however, in addition to the normal compensation, either a heavy computation burden arises from the calculation of the hyperbolic equations [7] or certain limitations have to be imposed on the system in order to simplify the calculation [8]. A two-end algorithm is also proposed [9], where the phasors at the two terminals are simply related by the voltage phasor at the fault location. The distance is then solved. In general, the required filtering time presents a trade-off between the speed and the accuracy, especially in the case of the fast relays with the relaying time of less than half a cycle [10]. A highly selective filter will always require a longer time delay regardless of the techniques to be applied [11–14]. The DEA uses the transient signals including the dc decay, the fundamental and non-fundamental components which satisfy a differential equation to calculate the fault location. The advantage is that the burden in filtering is reduced, increasing the relay speed and possibly the accuracy as well; however, in general, the DEA is sensitive to the fault inception angle, the sampling rate, the signal spectrum and the parameter variations with frequencies, and works well in the short lines only;

* Corresponding author. Tel.: +1 613 5416000; fax: + 1 613 5428612. 0378-7796/98/$19.00 © 1998 Published by Elsevier Science S.A. All rights reserved. PII S0378-7796(98)01231-5

2

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

while in the long lines, poor accuracy is observed due to neglecting the shunt capacitance in the line model. Despite the frequency dependence of the line parameters [15,16], R and L are always taken to be constant values for all frequencies lower than a cutoff frequency. Several methods have been proposed. A typical method is to integrate the differential equation in two intervals and solve for R and L directly [17]. Less computation burden and fast relay speed are obtained. Other methods have been proposed to improve the performance of the above algorithm by either eliminating certain harmonics [18] or reducing error in the difference approximation to the integration for certain harmonics [19]. The latter is especially suitable when a lower sampling rate is used. It should be noted that eliminating harmonics improves the approximation and reduces the effect of the shunt capacitance as well. However, since the shunt charging current due to even a lower frequency such as the fundamental in a long line represents possibly a countable error source, these algorithms are the best candidates for the short lines. Our objective in this work is to develop the algorithm which should provide both the fast relaying speed and the accurate fault location. A fast relay enables a shorter fault clearing time for the faults within the protection zone, while a longer relay time is often required for fault near the zone boundary or the balance point due to the error in the fault distance estimate. Thus, a higher accuracy in the fault location is very much desired, resulting in a larger and faster zone coverage. To this end, a relay speed of 1/4–1/2 cycles and an error within 1% of the total line length are determined for our algorithm. When such requirements are to be met, the previously discussed approaches need to be improved. An approach with the following two features seems attractive. First, the faulty system can be accurately modeled with the inclusion of the shunt capacitance. Second, the transient signals can be used without jeopardizing the accuracy of the result due to the shortcomings of the DEA previously discussed. An algorithm is proposed in this paper, which combines the strengths of the two distinct approaches (the SSPA and the DEA). In the proposed algorithm, a distributed line model is used, which is suitable for the short, the long, the transposed and the untransposed lines. The DEA is used to calculate the fault location and is implemented using the Laplace and the Z-transform techniques. The following advantages are observed in this proposed algorithm over the existing algorithms: 1. the variations in the line parameters with the frequency are considered in the algorithm; 2. the fault location estimate calculated in this algorithm is less sensitive to the various selections of the sampling rate, the fault inception angle and the

signal spectrum and he equivalent sources behind both busbars; and 3. the relay performance is greatly improved, offering both the desired high accuracy and speed at the same time. The proposed algorithm requires a data communication link between the two line terminals. The transmission delay; type of synchronization required and the data transfer capacity of the link will depend on both the transmission line characteristics (the transmission line length, the configuration of the line, the number of parallel circuits, etc.) and the type and complexity of the link. The paper is organized as follows. Section 2 presents a flow chart of the proposed algorithm. Several key aspects are discussed: 1. the criteria for the selection of the sampling rate and the frequency range; 2. the fine tuning process for controlling the error due to variations in the line parameters with frequency; and 3. the solution process to reduce error occurred in the final solution. Section 3 presents the simulation results using the EMTP. Section 4 presents the conclusion. Appendix A presents the derivation for the relationship between the fault location and the phase variables. Appendix B lists the expressions for some algorithm parameters.

2. Proposed algorithm Considering the sample system in Fig. 1, the line has parameters uniformly distributed along the length and the parameters can be calculated for any given frequency. In Appendix A, the relationship between the fault distance x F and the phase variables at both terminals is obtained as: F

F

(e − gi (v)x + egi (v)x )IS P(v)V Sp(v) F

F

+Z 0i (v)(e − gi (v)x − egi (v)x )IQ P(v)I Sp(v) F

F

= (e − gi (v)(1 − x ) + egi (v)(1 − x ))IS P(v)V R p (v) F

F

+ Z 0i (v)(e − gi (v)(1 − x ) − egi (v)(1 − x ))IQ P(v)I R p (v) (1) The notations can be found in Appendix A. As discussed before, the SSPA will solve the above complex equation for x F directly in the frequency domain. The differential equation approach will have to find a time domain expression for Eq. (1) that must use less time in filtering yet be able to provide the accuracy required.

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

Fig. 1. A sample line system with a fault.

2.1. Construction of the algorithm

Step (d7) Step (d8)

The proposed algorithm has been developed in the following flow chart, which consists of three phases shown as follows (Fig. 2):

2.1.1. Phase one: design phase

Step (d1) Step (d2)

input the line constants for a given line system input the error tolerance E, EV and ET and the relay time T

Step (d9) Step (d10)

Steps (d3–d11) determine the feasible frequency range for the Taylor expansion Step (d3) Step (d4) Step (d5) Step (d6)

j= 1 choose arbitrarily an initial frequency vj evaluate gi (vj ) find n from Eq. (11) for E and vj

Step (d11)

Steps (d12–d19) determine the frequency window Step (d12) Step (d13) Step (d14) Step (d15)

Step (d16) Step (d17)

Step (d18) Step (d19)

Fig. 2. Flowchart for the execution of the algorithm.

is n\2 yes! reducing frequency if j = 1, vj+1 = 0.628 vj, j= j+1, go to step (d5) if j\1 and vjBvj−1, then vj+1 = 0.6286vj, j= j−1−1, go to step (d5) if j\1 and vj\vj−1, then vj+1 =vj− 0.628(v;−vj−1 ), j =j+1, go to step (d5) no! Is n= 2? no! Then nB2, increasing frequency if j= 1, vj+1 = vj /0.628, j =j+l, goto step (d5) if j\1 and vjBvj−1, then vj+1 =vj− 0.628(vj−vj−1 ), j =j+1, go to step (d5) if j\1 and vj\vj−1, then vj+1 =vj / 0.628, j= j+1, go to step (d5) n= 2, vmax = vj, Wf = [0, vmax] obtained

determines the harmonics in Wf vi = i2p60, i= 1, …, k i= 1 evaluate the line parameters at v i estimate the parameter variation DPj (v) and the time step Dtj for EV using Eq. (15) for each parameter find the frequency range around v i: Wj = [vlj vuj ] for each parameter find Dt ip = min(Dtj ) and the intersection W i = (SWj )SWf = [v ic1 v ic2] , calculate the bandwidth vb i of W i using v ic2/v ic1 is i= k? If no, i= i+1, go to step (d14); if yes, go to step (d19) choose the maximum vb i, say with i= i0, the frequency window is finally selected as W 0 = [v ic10 v ic20 ], record w0 = w i0, vc1 = v ic10 and vc2 = v ic20 and Dt Pmax = Dt iP0

3

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

4

Steps (d20–d22) determine the bandpass filter Step (d20)

Step (d21) Step (d22)

input the stop band cut of frequencies [v ic1 v ic2] along with other filter specifications find the filter with the minimum order Nf} obtain the filter coefficients

Steps (d23–d24) determine the time step or the sampling rate Step (d23) Step (d24)

the maximum time step is Dtmax = z=Dt Pmax determine Dt using Eq. (18) such that Dt5Dtmax and TBTmax

2.1.2. Phase two: relay function Step (r1) Step (r2) Step (r3) Step Step Step Step Step

(r4) (r5) (r6) (r7) (r8)

input v0 and Dt evaluate all parameters ai, bi, RZ 0i , IZ 0i , RIS P and IIQ P at v 0 expand the exponential terms using Eq. (12) substitute the expansion into Eq. (9) jv= s {s =2[1−z−1]}/{Dt[1+z−1]} z−NG(z) = g(t0−NDt) obtain Eq. (13)

2.1.3. Phase three: solution process Step (s1) Step (s2) Step (s3) Step (s4)

Step (s5) Step Step Step Step Step

(s6) (s7) (s8) (s9) (s10)

Step (s11)

input ET, Ns., flag= 0, j, k, m = 0 is j =Ns? No, j =j+1, go to step (s3); if yes, go to step (s12) input the phase waveform, samples after the low pass evaluate the spectral phase samples using the filter coefficients obtained in step (d22) evaluate a, b, c in Eq. (13) using Appendix B solve Eq. (13) for x F evaluate x Fj using Eq. (20) if j =1 go to step (s2) is s Fj −x Fj-1 if yes and flag=0, then set flag=1 and k = k+1, record solutionk =x Fj , m =m+1 go to step (s2) if yes and flag= 1, then m= m+1, go to step (s2) if no and flag=0, go to step (s2) if no and flag= 1, then set flag=0, store the pair (solutionk =x Fj , number of solutions=m), set m =0, go to step (s2)

Step (s12)

Step (s13)

choose the solution with the largest number of solutions as the final estimate x F output the solution

2.2. Discussion In this section, the above algorithm is discussed in detail.

2.2.1. Signal conditioning The broad range of frequencies in the transients introduces variations in the parameters, gi (v), Z 0i (v)IS P(v) and IQ P(v) for i= 0, 1, 2. The earth mode parameters (i= 0) have much larger variations than the aerial mode ones (i= 1, 2). If the parameters are to be approximated by the values at a certain frequency, then variations in a large frequency range will cause great inaccuracy in the solution. For this reason, two measures are taken: 1. the aerial mode systems are used in the algorithm; and 2. in order to limit the variation, the range of the frequencies has to be limited. A bandpass filter will be applied for this purpose. The frequency response of the bandpass filter is denoted as W(v). By selecting the two stop band cutoff frequencies and the window response, the desired spectral signals are passed to the relay with the selected bandwidth, limiting, therefore, the parameter variations. To do so, Eq. (1) is multiplied by W(v), the resultant equation has the same form as Eq. (1) except that the phase variables are replaced with their corresponding conditioned spectral signals. F

F

(e − gi (v)x + egi (v)x )IS P(v)V. Sp(v) F

F

+Z 0i (v)(e − gi (v)x − egi (v)x )IQ P(v)I. Sp(v) F

F

= (e − gi (v)(1 − x ) + egi (v)(1 − x ))IS P(v)V. R p (v) F

F

+ Z 0i (v)(e − gi (v)(1 − x ) − egi (v)(1 − x ))IQ P(v)I. R p (v) (2) with V. Sp(v)= W(v)V Sp(v), R V. R p (v)= W(v)V p (v),

I. Sp(v)= W(v)I Sp(v) R I. R p (v)= W(v)I p (v)

(3) (4)

the superscript ^ indicating the conditioned spectral phase signals. Note that the shape of the response within the window does not affect the result, indicating that any type of filter can be applied regardless of the amount of ripples in the amplitude response between the two stop band cut off frequencies. This is shown in step (d21) of the flow chart, where the order of the filter becomes the only criteria in selecting a filter.

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

5

Fig. 3. Order–frequency curve for given error tolerance E.

2.2.2. Approximation of aerial parameters in the window As the window function is applied in Eq. (1), the spectral components in Eq. (2) have very small amplitude outside of the window due to large attenuation in the stop band. Thus, the values of parameters evaluated outside of the window have very little weight in the equation. Let the window be vc1 Bv B vc2 corresponding to the two stop band cut off frequencies, then, evaluating all the aerial parameters (i= 1, 2) at certain point v0 within the window results in gi (v)= ai (v0)+jvbi (v0)

(5)

Z 0i (v)= RZ 0i (v0)+ jvIZ 0i (v0)

(6)

IS P(v)=RIS P(v0)+ jvIIS P(v0)

(7)

IQ P(v)= RIQ P(v0)+ jvIIQ P(v0)

(8)

The selection of v0 will be discussed later. For simplicity, notations for the parameters will denote the constant values evaluated at v0 if (v) or (v0) is omitted. For example, ai, denotes ai (v0). These constant parameters are substituted into Eq. (2) to yield: (e

− (ai + jvbi )xF 0 i

+e

(ai + jvbi )xF

)(RIS P +jvIIS P)V. Sp(v)

0 i

+ (RZ +jvIZ )(e

− (ai + jvbi )xF

(ai + jvbi )xF

−e

)

(RIQ P +jvIIQ P)I. Sp(v) F

F

= (e − (ai + jvbi )(1 − x ) +e (ai + jvbi )(1 − x ))(RIS P +jvIIS P) 0 0 V. R p (v)+ (RZ i +jvIZ i ) F

F

(e − (ai + jvbi )(1 − x ) −e(ai + jvbi )(1 − x ))(RIQ P +jvIIQ P) I. R p (v)

(9)

The variation in each parameter due to the approxima-

tion in Eqs. (5)–(8) will be further discussed when the window is selected.

2.2.3. Taylor expansion for exponential terms in Eq. (9 ) Considering the nth order Taylor series expansion of e 9 g(v)l at origin gi (v)l= 0. e 9 gi (v)l = 19gi (v)l+ 12 g 2i (v)l 2 9 + … 1 + (9 1)ng ni (v)l n n!

(10)

Note that, gi (v)l= ail+ jvbi l, l is the line length. ai and bi are assumed to be evaluated at some point in the window. To start, choose arbitrarily a frequency vj, and evaluate them at vj. The error in a truncated expansion is to be investigated for the different frequencies v and orders of the expansion. The error occurred in this approximation is equal to the residual in the order of g ni − 1(v)l n + 1. The error function is denoted by s: d=





gi (v)l 1 +… (9 1)n + 1g ni + 1 (v)l n + 1 19 n+ 2 (n+1)!

(11) Note that this infinite summation always converges as there is the factorial term in every denominator. However, its magnitude, denoted by E, can be very large resulting from either a small n or a large exponent gi(v)l. The effect of the order n and frequency f on error magnitude E is analyzed in this work, which leads to the relationship between f and n for any given E as shown in (Fig. 3). Each curve gives the minimum order n at a frequency f for a given error tolerance E. The area above each curve represents all the feasible choices of the order for that specific error tolerance. With the above analysis, the feasible range of the

6

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

frequency window can be obtained according to the procedures listed in steps (d3 – d11) in the flow chart. One selects at first arbitrarily an initial frequency for which the minimum order is calculated. If the minimum order is B2, the next choice of frequency is to be increased, otherwise it is to be reduced. The step size is chosen as 0.628. It is found with this number, a maximum feasible frequency range vf from 0 to vmax can easily be located in this process for an order of 2. The choice of the second order approximation is due to the simplicity in finding a solution to a second order polynomial equation. Thus, the error in the Taylor expansion approximation can be controlled by the given error tolerance. The second order expansion as indicated in step (r3) is as follows: e 9 gi (v)x = 19(ai x+ jvbi x) + 12 (ai x +jvbi x)2

(12)

In order to apply the expansion to Eq. (9), one will first substitute x F or 1 −x F for x in Eq. (12) and then substitute the expansion for the respective terms in Eq. (9) as indicated in step (r4). A second order polynomial in x F and in jv is obtained with v varying in the frequency window. All the other parameters in the polynomial are constants to be evaluated at certain point in the window. Now, a form is obtained with which the frequency domain relationship can be easily transformed to the time domain as shown in steps (r5 – r8). The sequence is as follows: 1. from the v-domain to the s-domain, using jv= s requires zero initial conditions on all the variables involved. The superimposed post fault transients are then used; 2. from the s-domain to the z-domain. The bi-linear transform s = {2[1− z − 1]}/{Dt[1 +z − 1]} is used because of its anti-aliasing feature. Dt is the sampling interval. 3. finally, from the z-domain to the time domain, using z − NG(z) = g(t0 −NDt). G(z) denotes a function in the z-domain and g(t) its corresponding function in the time domain, to the current sampling instant and g(t0 −Ndt). G(z) the Nth historical sample of g(t) back in time from the current sample. In the resultant polynomial, x F has a second order, while N ranges from 0 to 3, g(t) denotes the corresponding spectral phase signals 6ˆ Sp(t), i. Sp, 6ˆ R p (t) and i. R (t). Therefore, four consecutive samples for each p spectral phase signal are needed to yield one estimate for x F by solving an equation of the following form: a(x F)2 + bx F +c=0

(13)

where, a, b and c are shown in Appendix B. It should be noticed that a new estimate for x F is generated for each set of new samples.

With Eq. (13) available the remaining work in the flow chart can be done easily.

2.2.4. The selection of frequency window in steps (d12) Having obtained the feasible frequency range in step (3) as Wf = [0, vmax] given the order n=2, the frequency window is to be found. It should be noticed that if a higher order n can be used, the range will be larger and more flexibility of selecting a particular frequency window is rendered. Factors in general influencing the selection of the window within Wf are among the following: (a) The maximum variation of the model parameters allowed within the window. An error due to the variation DPj (v) in a parameter, denoted by, Pj (v) is defined by EV(DPj (v)) as the first order deviation of x F with respect to DPj (v) around a given frequency v0. From Eq. (13), the first order deviation of x F due to DPj (v) can be calculated as: (a F 2 (b F (c (x ) + x + DPj (v) (Pj (Pj (Pj EV(dPj (v))= (14) 2ax F + b





Examining the coefficients in Eq. (14), EV(DPj (v)) can be further expressed in the form of EV(DPj (v))=

q3Dt 3 + q2Dt 2 + q1Dt +q0 DPj (v) r3Dt 3 + r2Dt 2 + r1Dt +r0

(15)

where qi and ri for i= 0, 1, 2, 3 consists of the linear combinations of the spectral signal samples. The coefficients of the combinations are determined by the model parameters at v0 and the distance. Given an error tolerance EV, (DPj (v)) will be found for the maximum variation allowed in parameter Pj (v) and the corresponding frequency range around v0 will be located. Several observations can be made regarding this process. (i) The magnitude of EV(DPj (v)) is controlled overall by DPj (v) and Dt as well, if q0 = 0 and r0 "0. (ii) If for each parameter, Pj (v)), the Dt, and is to be chosen small enough such that: 3

% qi Dt ij  q0, i=1

3

% ri Dt ij  r0

(16)

i=1

then approximately EV(DPj (v))=

q0 DPj (v) r0

(17)

The value of q0 or r0 is determined by the samples and the coefficients, and it also changes for each new sample. However, their magnitude and hence the ratio between them may be approximately estimated. Thus, given EV for all parameters, DPj (v) may be estimated. Two frequencies vlj and vuj can be chosen such that for a maximum bandwidth vlj B vBvuj, one has EVDPj (v) B %EV. The intersection of the

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

7

Table 1 Simulation results, length (km), angle (°), time (ms) Fault type

Fault angle

Actual location

Calculated location

Relay time

Error (%)

a a a ac abc

90 0 0 90 90

56.30 144.80 104.60 56.30 56.30

56.21 144.36 104.60 55.86 55.91

4–8 4–8 4–8 8 8

0.15 0.31 0 0.77 0.69

bands obtained for all the parameters is found to be the feasible bandwidth around v0. (iii) Since q0 and r0 are changing for each new sample, EVDPj (v) will vary for certain new samples and may be very large due to a large ratio of q0 to r0. This suggests that certain solution processes be applied to all the estimates of the fault distance such that the error in the final solution can be reduced. The above analysis is implemented in steps (d12– d19). Several points need further explanation: (i) k frequencies v i, i=1, 2, …, k in Wf are used in Eq. (15) for each parameter Pj (v). A largest bandwidth around one of the frequencies is to be located. These points have been selected as the fundamental frequency and its harmonics. For each frequency v i, a window Wj and a time step Dtj used for the parameter Pj (v) are both obtained. When this process is done for all parameters, the smallest time step is recorded as Dt iP intersection of the feasible range Wf and all the windows at v i are recorded as W i: v ic1 Bv B v ic2 and vb i =v ic2/v ic1 is also calculated as the window bandwidth. (ii) the process is then repeated for the rest of the frequencies until i= k. When it is done, one window and one time step are recorded for each frequency. Choose the window with the largest vb i, say with i= i0. i0 Then the window finally selected is W 0 =[v i0 c1 v c2]. The i0 corresponding v is the v0 at which the aerial parameters are to be evaluated in the algorithm, and v i0 c1 and are v and v for the stop band lower and upper v i0 c2 cl c2 cut off frequencies of the filter, respectively. Also the P corresponding time step Dt i0 P is recorded as Dt max.

will be selected. For example, an elliptic AR filter will be chosen compared to a Butterworth filter. The coefficients of the filter are stored and used to obtain the spectral phase samples. This process is shown in steps (d20–d22).

2.2.5. The bandpass filter designed in steps (d20 – d22) The bandpass filter can be designed with the selected stop band cut off frequencies and the given specifications such as the stop band attenuation and the roll off slope. As discussed before, the ripple in the pass band does not affect the accuracy. The roll off slope does not affect the accuracy either, since the parameter variation is within its limit between the stop band cut off frequencies. However, it should be noted that large ripples and small rolls of slope should be avoided for sufficient signal:noise ratio. In general, when the filter specification is given, the type of filter with the minimum order

1 x Fj = (( j− 1)x Fj − 1 + x F(t0)) (20) j t0 indicating the most current sample. An error tolerance, ET, is set for the convergence. When the difference between each consecutive mean x Fj is within the ET, a counter starts to record the period (number of the estimates) during which the difference remains in ET and also records the mean estimate for this period. When the difference is greater than the ET, the previous mean estimate and the period are stored and the counter resets to zero and is ready to record the next stabilized solution. This process continues until the

2.2.6. The selection of time step in steps (d23 –d24) Dt Pmax is already obtained from the window selection as shown in step (d19). The actual time step is obtained as follows. Given the maximum relay time Tmax required, the actual relay time T is: T= NfDt +3Dt + NsDt

(18)

with Nf the order of the bandpass filter, Ns, the number of estimates needed in the solution process. Thus, NfDt at is the filter time delay and 3Dt is approximately the time delay for solving Eq. (13) and Nst is approximately the time delay in the solution process. The time step Dt is to be chosen such that Dt BDtmax and T BTmax. It has been found that the algorithm is not sensitive to the selection of Dt and Ns, due to the solution process proposed for this algorithm. In the numerical test, Dt between 50 and 500 ms has been tested, all with very stable results.

2.2.7. Solution process in steps (s1 –s13) The mean value of the estimates is calculated. x Fj =

1 0 % x F(t0 − iDt) j i=j−1

(19)

or in the recursive form:

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

8

specified maximum relay time is reached. The mean with the longest period is taken as the final fault distance estimate. The numerical simulation has been done to examine the algorithm performance regarding the selection of frequency window, time step and number of solutions in addition to its performance in different systems and faults. A particular design of the algorithm is presented in the next section.

3. Application of the proposed method The results in Table 1 are obtained by applying the proposed method to a 400 kV, 209 km single circuit, three-phase overhead transmission line with no transposition. This line system has a flat-spaced horizontal configuration for the phase conductors. Each phase consists of two bundled conductors; while a solid sky ground wire is located above the middle phase. The design presented here uses a feasible frequency range of [0 100] Hz for an error tolerance 0.01 in a 2nd order exponential expansion. The window is selected as 30B f B 100 Hz for a solution error within 1%. A bandpass elliptic filter of order 10 is used with the sampling interval 50 ms which gives  60 samples in a quarter of a cycle. As shown in Table 1, the simulation using the EMTP exhibits fast and accurate fault distance calculation. The relay time is within half a cycle and the error in the fault distance is within 1%. The fault resistance is taken to be 20 V. A variety of tests have also been conducted. Different line lengths are used. The source impedances are also changed. The X/R ratio varies from 40 to 10 in different cases. Due to the limited space and the length of this paper, they are not presented here. Nevertheless, consistent results are observed. Regarding the variables in the algorithm, in general, the frequency window is around the fundamental in order to maintain a second order Taylor series expansion. When the time step is as large as 500 ms, the order of the filter is reduced to 6. The estimates are very stable without any large deviation from each other even though only a small number of solutions are available. When a small time step is used, the oscillation in estimates is observed; however, a large number of solutions are provided and the new counting and comparing solution process is very effective in dealing with the oscillation using a large number of data. Thus, a high accuracy is obtained in both cases for a relay time within half a cycle. Overall, the test results have been very consistent and promising.

4. Conclusions An algorithm is proposed in this paper to meet the

objective of a fast relaying time (1/4–1/2 cycle) and an accurate fault location (error within 1% of the total line length). The algorithm can be used in the short, the long the transposed and untransposed lines, and can be used with either an on-line digital protective relay or an accurate off-line fault locator. This is achieved mainly due to the following features of this algorithm: (a) the combination of an accurate line model with a fast time domain solution method is made possible through the use of the Laplace and the Z-transforms; (b) the stable performance of the algorithm with respect to the various sampling rate, the fault inception angle, the signal spectrum and the equivalent source impedances behind the both busbars; (c) the frequency dependence of the line parameters is considered, allowing the appropriate selection of the transient spectral signals; (d) the better approximation of the Taylor expansion and the proper selections of the frequency window and the time step allow error to be well controlled in the final solution; (e) simple and fast band pass filters can be used; and (f) the new counting and comparing scheme suggested in this algorithm is suitable for the time domain approaches. This method deals effectively with the oscillation in the instant estimate. Bad data with a large deviation from the mean can be rejected while the most stabilized mean estimate is chosen.

Appendix A. Derivation of Eq. (1) With a three-phase distributed parameter line model, a set of differential equations in time domain for phase variables at a distance x i from the sending end can be obtained: (6 (x, t) (i (x, t) (A1) − p = Rip(x, t)+ L p (x (t −

(ip(x, t) (6 (x, t) = G6p(x, t)+ C p (x (t

(A2)

6p(x, t) and ip(x, t) are 3× 1 vectors. R, L, G and C are 3× 3 matrices of the line series resistance and inductance, and the shunt conductance and capacitance in per unit length, respectively. The frequency dependence of R and L can be dealt with by taking the Fourier transform of Eqs. (A1) and (A2) as follows: −

(Vp(x, w) =(R(v)+jvL(v))Ip(x, v) (x

(A3)

−

(Ip(x, w) =(G+ jvC)Vp(x, v) (x

(A4)

The upper case V and I indicate the frequency domain variables. Zp(v)= R(v)+ jwL(v) and Yp(v) =G+ jvC are the 3 × 3 phase series impedance matrix and

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

the shunt admittance matrix in per unit length, respectively. The following equations can be shown equivalent to Eqs. (A3) and (A4): ( 2Vp(x, v) = Zp(v)Yp(v)Vp(x, v) (x 2

(A5)

gi (v)= ai (v)+ jvbi (v)

(A13) 0 i

and the surge impedance for the mode i, Z (v), can be found as the elements of the diagonal matrix Z 0(v) whereas Z 0(v)=g − 1(v)S − 1(v)Zp(v)Q(v)

(A14)

The superscript − 1 indicates the matrix inverse. Given the line length l, Ai (v) and Bi (v) can be obtained from either the sending end modal voltage V Smi(v) and current I Smi(v) by setting x =0 or the R receiving end modal voltage V R mi(v) and current I mi(v) by setting x = 1 in Eqns. A11 and A12. The upper case superscript indicates the location on the line, S, R and F being the sending end, receiving end and the fault point, respectively. Correspondingly, A Si (v), and B Si (v) or B Si (v) and B R i (v) can be obtained. Refer to Fig. 1, the voltage at the fault location V Fmi(v) can then be

9

obtained from Eqn. A11 by setting x= x F and using R either A Si (v) and B Si (v) or A R i (v) and B i (v). That is, F V mi(v) will be expressed simultaneously in two equations. Equating the two resultant equations and substituting the modal variables with the phase variables using V Smi (v)= IS P(v)V Sp(v)

(A15)

I Smi (v)= IQ P(v)I Sp(v)

(A16)

and P R VR mi (v)= IS (v)V p (v) R mi

P

(A17)

R p

I (v)= IQ (v)I (v)

(A18)

with IS P(v) and IQ P(v) for p =0, 1, 2 indicating the rows of matrices S − 1(v) and S − 1(v), respectively, one finally obtains F

F

(e − gi (v)x + egi (v)x )IS P(v)V Sp(v) F

F

+Z 0i (v)(e − gi (v)x − egi (v)x )IQ P(v)I Sp(v) F

F

= (e − gi (v)(1 − x ) + egi (v)(1 − x ))IS P(v)V R p (v) F

F

+ Z 0i (v)(e − gi (v)(1 − x ) − egi (v)(1 − x ))IQ P(v)I R p (v) (A19)

Appendix B. Coefficients in Eq. (13)

!

"

The coefficients of Eq. (13) can be obtained as follows: a = RIS P

!

"

4b 2i S 8b 2i S S S S P . . . . [y (0) −y ( − 1) −y ( −2) +y ( − 3)] + IIS [y. p(0)− 3y. Sp(− 1)− 3Sp(− 2)+ y. Sp(− 3)] p p p p Dt 2 Dt 3

!

S R S R S R +RIS P a 2i [y. Sp(− 0) −y. R p ( −0) +3y. p( − 1) −3y. p (− 1)+ 3y. p(− 2)− 3y. p (− 2)+ y. p(− 3)− y. p (− 3)]

+

!

2a 2i S S R S R S R [y. p(0) −y. R p (0) +y. p( −1) −y. p ( −1)− y. p(− 2)+ y. p (− 2)− y. p(− 3)+ y. p (− 3)] Dt

+IIS P +

"

4aibi S S R S R S R [y. p(0)− y. R p (0) +y. p( − 1) −y. p ( − 1) −y. p(− 2)+ y. p (− 2)− y. p(− 3)+ y. p (− 3)] Dt

8aibi S S R S R S R [y. p(0)− y. R p (0) −y. p( − 1) +y. p ( −1) − y. p(− 2)+ y. p (− 2)+ y. p(− 3)− y. p (− 3)] Dt 2

!

"

(B1)

ˆ Sp( − 1) +3ıˆ R ˆ Sp(− 2)+ 3ıˆ R ˆ Sp(− 3)+ 3ıˆ R b= −RIQ P 2ai IRZ 0i [ıˆ Sp(0) +ıˆ R p (0) +3ı p (− 1)+ 3ı p (− 2)+ 3ı p (− 3)] + +



n

4bi l 4a l RZ 0i + i IZ 0i [ıˆ Sp(0) +ıˆ R ˆ Sp( − 1)+ ˆı R ı Sp(− 2)− ˆı R ı Sp(− 3)− ˆı R p (0) +ı p (− 1)− ˆ p (− 2)− ˆ p (− 3)] Dt Dt

8bi l 0 S IZ i [ıˆ p(0)+ ˆı R ˆ Sp( −1) −ıˆ R ı Sp(− 2)− ˆı R ı Sp(− 3)− ˆı R p (0) −ı p ( −1) − ˆ p (− 2)+ ˆ p (− 3)] Dt 2

!

− IIQ P +



4ai l RZ 0i [ıˆ Sp(0) +ıˆ R ˆ Sp( − 1) +i R ı Sp(− 2)− ˆı R ˆ Sp(− 3)− ˆı R p (0) +ı p ( −1)− ˆ p (− 2)ı p (− 3)] Dt

n

8bi l 8a l RZ 0i + i2 IZ 0i [ıˆ Sp(0) +ıˆ R ˆ Sp( − 1)+ ˆı R ı Sp(− 2)− ˆı R ı Sp(− 3)− ˆı R p (0) +ı p (− 1)− ˆ p (− 2)− ˆ p (− 3)] 2 Dt Dt

!

R R R +RIS P 2a 2i l[y. R p (0) +3y. p ( − 1) +3y. p ( − 2) +y. p (− 3)]+

"

"

8aibi l R R R [y. p (0)+ y. R p (− 1)+ y. p (− 2)+ y. p (− 3)] Dt

10

!

+ IIS P

B. Lian et al. / Electric Power Systems Research 46 (1998) 1–10

"

16aibi l R 4a 2i l R R R R R [y. p (0)− y. R [y. p (0)+ y. R p ( − 1) −y. p ( − 2) −y. p ( − 3)]+ p (− 1)− y. p (− 2)+ y. p (− 3)] Dt Dt 2

S R c= RIS p[y. Sp(0)−y. R p (0) +y p( − 1) −y. p ( − 1)] +

! ! !

2 S R IIS p[y. Sp(0)− y. R p (0)− y. p(− 1)+ y. p (− 1)] Dt

R R − RIS p a 2i l 2[y. R p (0) +2y. p ( − 1) +y. p ( − 2)] +

− IIS p

(B2)

4aib 2i l 2 R [y. p (0)− y. R p (− 2)] Dt

" "

2a 2i l 2 R 8aibi l 2 R R [y. p (0) −y. R [y. p (0)− 2y. R p ( − 2)] + p (− 1)+ y. p (− 2)] Dt Dt 2

− RIQ 2ai lRZ 0i [ıˆ R ˆR ˆR ˆR p (0) +3ı p ( −1) +3ı p ( − 2) +ı p (− 3)] + +



n

4ai l 0 4bi l ˆR ˆR ıR RZ 0i + IZ i · [ıˆ R p (0) +ı p ( − 1) −ı p (− 2)− ˆ p (− 3)] Dt Dt

+ IIQ p + +

"

8bi l 0 R IZ i [ıˆ p (0)−ıˆ R ˆR ˆR p ( − 1) −ı p ( − 2) +ı p ( − 3)] Dt 2



!

4ai l RZ 0i [ıˆ R ˆR ˆR ˆR p (0) +ı p ( − 1) −ı p ( − 2) −ı p (− 3)] Dt

n

8bi l 8a l RZ 0i + i2 IZ 0i · [ıˆ R ˆR ˆR ıR p (0) −ı p ( −1) −ı p (− 2)+ ˆ p (− 3)] 2 Dt Dt

l6bi l 0 R IZ i [ıˆ p (0) −3ıˆ R ˆR ˆR p ( − 1) +3ı p (−2) −ı p (−3)] Dt 3

"

(B3)

All the parameters are constant values at v0 and the coefficients of the spectral variables can be calculated off-line to form an efficient computation scheme. References [1] A. Wiszniewski, Accurate fault impedance locating algorithm, IEE Proc. 130 (6) (1983) 311–314. [2] L. Eriksson, M.M. Saha, G.D. Rockefeller, An accurate fault locator wioth compensation for apparent reactance in the fault resistance negative sequence currents, IEEE Trans. Power Deliv. 5 (1) (1990) 79 – 84. [3] A.G. Phadke, M. Ibrahim, T. Hlibka, Fundamental basis for distance relaying with symmetrical components, IEEE Trans PAS PAS-96 (2) (1977) 635–646. [4] R.J. Marttila, Effect of transmission line loading on the performance characteristics of polyphase distance relay elements, IEEE Trans. Power Deliv. 3 (4) (1988) 1466–1474. [5] A.A. Girgis, D.G. Hart, W.L. Peterson, A new fault location technique for two- and three-terminal lines, IEEE Trans. Power Deliv. 7 (1) (1992) 98–107. [6] M.S. Sachdev, R. Agarwal, A technique for estimating transmissio line fault location from digital impedance relay measurement, IEEE Trans. Power Deliv. 3 (1) (1988) 121–129. [7] S.E. Westlin, J.A. Bubenko, Newton–Raphson technique applied to the fault location problem, IEEE PES Summer Meeting, Portland, OR, 18 – 23 July, 1976, Paper No. A 76 334–3. [8] T. Takagi, Y. Yamakoshi, M. Yamaura, R. Kondow, T. Matsushima, Development of a new type fault locator using the one-terminal voltage and current data, IEEE Trans. PAS PAS101 (8) (1982) 2892 –2898. [9] A.T. Johns, S.M. Haden, New accurate transmission line fault location equipment, IKE 4th Int. Conf. on Developments in Power System Protection, University of Edinburgh, UK, 11 – 13 .

.

April 1989, Publication No. 302, pp. 1 – 5. [10] J.S. Thorp, A.G. Phadke, S.H. Horowitz, J.E. Beehle, Limits to impedance relaying, IEEE Trans. PAS PAS-98 (1979) 246 –260. [11] B.J. Mann, I.F. Morrison, Digital calculation of impedance for transmission line protection, IEEE Franz PAS PAS-90 (1) (1971) 270 – 279. [12] G.R. Slemon, S.D.T. Robertson, M. Ramamoorty, High speed protection of power systems based on improved power system models, Int. Conf. on Large Electric Systems at High Tension (GIGRE), 22nd session, 1968, Paper no. 31-09. [13] M.S. Sachdev, N. Nagpal, A recursive least error squares algorithm for power system relaying and measurement applications, IEEE Trans. Power Deliv. 6 (3) (1991) 1008 – 1015. [14] A.A. Girgis, R.G. Brown, Application of Kalman filtering in computer relaying, IEEE Trans. PAS PAS-100 (7) (1981) 3387– 3397. [15] J.R. Carson, Wave propagation in overhead wires and ground return, Bell Sys. Tech. J. 5 (1926) 539 – 554. [16] L.M. Wedepohl, Application of matrix methods to the solution of traveling-wave phenomena in polyphase systems, Proc. IKE 110 (12) (1963) 2200 – 2212. [17] R. Poncelet, The use of digital computers for network protection, CIGRE, 24th Session, Paris, 1972, Paper no. 32-08. [18] A.M. Ranjbar, B.J. Cory, An improved method for the digital protection of high voltage transmission lines, IEEE Trans. PAS PAS-94 (2) (1975) 544 – 550. [19] Y. Ohura, T. Matsuda, M. Suzuki, M. Yamaura, Y. Kurosawa, T. Yokoyama, Digital distance relay with improved characteristics against distorted transient waveforms, IEEE Trans Power Deliv. 4 (4) (1989) 2025 – 2031.