- Email: [email protected]
Distance protection using a wavelet-based filtering algorithm
Electric Power Systems Research 80 (2010) 84–90
Contents lists available at ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
Distance protection using a wavelet-based filtering algorithm Kleber M. Silva a,∗ , Washington L.A. Neves b , Benemar A. Souza b a b
University of Brasília, Electrical Engineering Department, 70919-970 Brasília, Brazil Federal University of Campina Grande, Paraiba, Brazil, Electrical Engineering Department, 882 Aprigio Veloso Ave, Bodocongo, 58109-970 Campina Grande, Brazil
a r t i c l e
i n f o
Article history: Received 2 December 2008 Received in revised form 1 July 2009 Accepted 21 August 2009 Available online 27 September 2009 Keywords: High-speed distance protection Wavelet transform Phasor estimation
a b s t r a c t This paper presents a novel filter design technique based on the maximal overlap discrete wavelet transform (MODWT). Time and frequency response of the designed filters are compared to the ones of traditional discrete Fourier transform (DFT) based filters. The obtained results show improvements on response speed when comparing to full cycle traditional filters and improvements on frequency response when compared to half cycle traditional filters. According to these characteristics, the designed filter may be used in secure high-speed distance protection, since it provides fast relay operating times, without compromising their frequency response. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Due to the growth of the complexity and interconnectivity of the HV and EHV power systems, modern protection schemes must be designed to clear faults rapidly in order to maintain the system stability and enhance the system security. Thus, it is necessary to detect a fault within sub-cycles of the power system frequency [1]. Distance protection relays are widely used to protect transmission lines. Their operation depends on both voltage and current phasors at the fundamental frequency. In fact, the accuracy, speed and security of these relays are directly related to the operation of the digital filtering algorithms used to estimate phasors [2], which must have certain characteristics, such as: bandpass response about the system frequency, decaying DC rejection, harmonic attenuation or rejection and good transient behavior [3]. The most popular phasor estimation algorithms are the full cycle DFT (FCDFT) and the full cycle least error square (LES) algorithms [4]. However, they are not suitable for high-speed distance protection, because relays are required to operate in less than one cycle. In this way, many short window algorithms have been reported [5,6]. Although these algorithms are very fast, they are very inaccurate because they do not eliminate harmonics. Another important fast algorithm is the half cycle DFT (HCDFT) [4], but it also suffers from inaccuracies because it neither eliminates even harmonics nor rejects decaying DC components. Many improved half cycle algorithms have been reported [7,8], but they do not eliminate even harmonics. To overcome these drawbacks, some adaptive algorithms have been developed, which varies the
window length after detecting the fault [9,10]. However, they are slower than the HCDFT algorithm. Alternatively, the discrete wavelet transform (DWT) has been also used to develop innovative phasor estimation algorithms [11,12]. In [11], the authors use the orthogonality property of the DWT to compute synchronized phasors. This method was used in [13] to propose a new DWT-based transmission line distance protection algorithm. However, this algorithm cannot be used in high-speed distance protection, since it does not provide a fast relay operating time and is heavily affected by harmonic and decaying DC components. An algorithm based on the redundant version of the DWT, named MODWT, was proposed in [12]. It uses the good signal approximation characteristic provided by the MODWT, computing the phasor by means of a two-samples algorithm applied to the smooth version of the original voltage and current signals. The reported results show that this algorithm is faster than FCDFT. Unfortunately, the authors did not analyze the frequency response of the overall filtering system. This paper presents a novel filter design technique based on the MODWT. This technique may be understood as a generalization of the phasor computation algorithm proposed in [12]. Time and frequency responses of proposed filters are compared to the ones of traditional DFT-based filters. The obtained results shown that the proposed and HCDFT filters have similar time responses, with the advantage that the frequency responses of the proposed filters are close to the ones of FCDFT filters. It indicates that the filters developed here are suitable for secure high-speed distance protection.
2. MODWT fundamentals ∗ Corresponding author. Tel.: +55 61 3307 2308; fax: +55 61 3274 6651. E-mail address: [email protected] (K.M. Silva). 0378-7796/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2009.08.012
The MODWT is similar to the DWT in that both produce a set of time-dependent wavelet and scaling coefficients, but the MODWT
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
85
is a highly redundant nonorthogonal transform, retaining downsampled values at each level of the decomposition that would be otherwise discarded by the DWT. In addition, the DWT cause more distortions in signal reconstruction and requires that the total number of samples must be a power of two [12,14]. Using wavelet transform a signal is decomposed into smoothed and detailed versions for different levels of resolution. For instance, a discrete signal X may be written using its smoothed S1 and detailed D1 versions at the 1st MODWT level [14]: X = S1 + D1
(1)
and using its smoothed S2 and detailed D2 versions at the 2nd MODWT level: X = S2 + D2 + D1
(2)
and so on. In this decomposition the signals Sj and Dj are computed as: Sj = Mj X
and
Dj = Nj X,
(3)
where the matrices Mj and Nj are obtained for each mother wavelet as shown in Appendix A. This analysis may also be understood as a filtering process, in which each MODWT level of decomposition divides the bandwidth of X into octave bands, as shown in Fig. 1 that depicts the frequency response of the 16th rows of the matrices N1 , N2 , N3 , M1 , M2 and M3 computed using the Daubechies 8 mother wavelet [14], for a generic sampling frequency fs . The rows of the matrix Nj work as low-pass filters, whereas the rows of the matrix Mj work as bandpass filters. Since the focus are the phasors at the fundamental frequency, only matrices Mj are analyzed. A discrete signal X is considered to illustrate the MODWT multiresolution analysis: xk = sin ωtk +
1 1 sin 3ωtk + sin 5ωtk , 3 5
Fig. 2. Original signal x and its smoothed version S3 at 3rd MODWT level, obtained using the mother wavelet Daubechies 8.
Fig. 3. Frequency response of the 16th, 24th and 32nd rows of M3 using Daubechies 8.
(4)
where xk is the k th sample of the signal X, ω = 2f is the fundamental angular frequency, f the fundamental frequency and tk is the k th sampling interval. Fig. 2 illustrates one cycle of the original signal X and its smoothed version S3 , computed using the mother wavelet Daubechies 8 and assuming the sampling rate of 32 samples per cycle. It is observed that the smoothed version S3 contains only the fundamental frequency component. In fact, the 3rd and 5th harmonics of X are eliminated by the rows of M3 , as shown in Fig. 3, which depicts the frequency response of the 16th, 24th and 32nd rows of M3 . In the same way, the rows of the matrices M2 and M1 eliminate some harmonics too.
From the above discussion, it is clear that the rows of the matrices Mj eliminate some harmonics. This filtering characteristic is quite suitable for phasor estimation and is taken into account in the proposed filter design technique. Further details are discussed next. 3. Proposed filter design technique In order to present the proposed filter design technique, consider one cycle of a sinusoidal signal X, which has N samples within a cycle. Its smoothed version Sj at the j th MODWT level may be computed by Eq. (3), in which the order of the matrix Mj is N × N. Thus, the signal Sj have also N samples, which may be written as: sj,k =
H
(Ych sin hωtk + Ysh cos hωtk ) + εk
(5)
h=1
where sj,k is the k th sample of the signal Sj ; H is the total of specified harmonics; Ych and Ysh are, respectively, the real and imaginary parts of the h th harmonic and εk is the estimation error of the k th sample. The Eq. (5) may be applied to M consecutive samples of the signal Sj , with M ≤ N, resulting in a linear system of equations: AY + E = B,
(6)
where A is a M × 2H matrix:
⎡
Fig. 1. Frequency response of the 16th row of the matrices N1 , N2 , N3 , M1 , M2 and M3 computed to Daubechies 8.
sin ωt0 ⎢ sin ωt1 A = ⎢. ⎣ .. sin ωtM−1
cos ωt0 cos ωt1 .. . cos ωtM−1
··· ··· .. . ···
sin Hωt0 sin Hωt1 .. . sin HωtM−1
⎤
cos Hωt0 cos Hωt1 ⎥ ⎥ .. ⎦ . cos HωtM−1
86
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
Y is a 2H × 1 vector of variables:
4.3. On choosing the Mj matrix
Y = Yc1
and E is a M × 1 vector of estimation errors:
Once the mother wavelet and the sampling rate are defined, the total number of MODWT levels is computed for one cycle of the signal X (Appendix A). The Mj matrix to be used is the last one. For example, considering a sampling rate of 32 samples per cycle and the mother wavelet Daubechies 16, it is possible to compute two MODWT levels. Thus, M2 is the chosen matrix.
E = ε0
4.4. Harmonic filtering
Ys1
· · · YcH
YsH
T
,
B is a vector with M consecutive samples of the signal Sj :
B = sj,0
sj,1
ε1
sj,2
ε2
· · · sj,M−2
· · · εM−2
sj,M−1
εM−1
T
T
.
The solution of Eq. (6) that minimizes the sum ET E, is obtained using the left pseudo-inverse A+ matrix [4]: Y = A+ B,
(7)
As discussed before, the rows of the matrices Mj eliminate some harmonics. For example, the 16th row of the matrix M2 in Fig. 1 eliminates all harmonics beyond the 6th harmonic, thus it is necessary to consider in the model only the 1st up to the 5th harmonics. The frequency response of the Mj rows defines which harmonics must be considered in the signal estimation model.
(8)
4.5. On choosing the M rows of the matrix Mj
−1
where A+ = AT A AT . From Eqs. (3) and (7), one can obtain:
Y = A+ Mj
M
X = GX,
G
where [Mj ]M is formed by M consecutive rows of the matrix Mj , in such way that G is a 2H × N matrix. The matrix G in Eq. (8) is the base of the proposed filter design technique. It is computed in advance in an off-line mode. The real part hc and imaginary part hs filters proposed to estimate the fundamental frequency phasors are the 1st and 2nd rows of G, respectively. The frequency response of the Mj rows defines the designed filters. It should be noticed that these filters are full cycle filters. Experience has shown that the more samples within one cycle are used, the more harmonics are filtered out, that is the reason why all samples within one cycle were used. The proposed technique may be understood as a generalization of the phasor computation algorithm proposed in [12]. This algorithm computes fundamental phasors using equations, formulated specifically to two consecutive samples of the smoothed signals at the 3rd MOWDT level, using the mother wavelet Daubechies 8. In fact, one can see that these equations are a particular case of Eq. (8). 4. Filter design considerations This section discusses general design considerations that must be taken into account to design filters that present good time and frequency responses. 4.1. Sampling rate As it can be seen in Figs. 1 and 3, the frequency responses of the rows of the Mj matrix have zeros. An appropriate choice of the sampling rate can put those zeros to work to notch out the harmonics. This can be achieved by choosing a sampling rate that is multiple of the fundamental frequency.
It is desirable that the chosen consecutive M rows of matrix Mj lead to designed filters hc and hs with good frequency response and small time delays. For instance, in Fig. 3, the 32nd row of M3 leads to the smallest time delay and to the poorest frequency response. On the other hand, the 16th row of M3 leads to a smoothest frequency response but to a larger time delay. In this case, the 24th row would be chosen, since it leads to good frequency response and a small time delay. M is directly related to the total number of harmonics to be reject in the signal estimation model, in such way that the more harmonics to be reject, the more consecutive rows must be taken into account. 4.6. Decaying DC component Traditionally, in LES algorithm applications, the decaying DC component is estimated by including the first terms of its Taylor series expansion in Eq. (5). However, in the cases discussed here, it leads to time delay in phasor computation. In this paper, the traditional digital mimic filter, which is a first order high-pass filter with unitary gain at the fundamental frequency, is used to eliminate the decaying DC component [15]. The final characteristics of the proposed filters are the combination of the one of the mimic filter with those of hc and hs filters. 5. Designed filters evaluation This paper presents the analysis of three filter pairs designed using the proposed technique. Table 1 shows the chosen design parameters: the mother wavelet, the MODWT level, the rows of Mj and the harmonics included in Eq. (5). In special, the filter pair FP1 is a particular case leading to the phasor computation algorithm proposed in [12]. The frequency and time responses of these filters were evaluated and the obtained results are discussed next. 5.1. Frequency response analysis
4.2. Mother wavelet The proposed filter design technique can be applied to any mother wavelet. In fact, different mother wavelets (Daubechies, Symlets, Coiflets, etc., [14]) with the same order leads to designed filters with nearly the same frequency responses. For example, filters designed using mother wavelets Daubechies 8 and Symlet 8 are quite similar, but those designed using mother wavelets Daubechies 8 and 4 are quite different. In addition, it was observed that the larger the order of the mother wavelet is, the better frequency responses the designed filters hc and hs have.
In this section, the frequency response of designed filter pairs are compared to the ones of traditional Fourier HCDFT and FCDFT filters. In order to do that, the designed FP1, FP2 and FP3 and Fourier filters use the same 2 cycle digital mimic filter [15]. The sampling rate is 32 samples per cycle. The frequency responses of FCDFT and HCDFT filters combined with mimic filter are shown in Fig. 4[4]. It is well known that FCDFT filters eliminate all harmonics, whereas HCDFT filters do not eliminate even harmonics. In addition, HCDFT filters are more affected by decaying DC and off nominal frequency components [3].
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
87
Table 1 Parameters used to design the evaluated filter pairs. Filter pair
Mother wavelet
MODWT level j
Rows of Mj
Harmonics in Eq. (5)
FP1 FP2 FP3
Daubechies 8 Daubechies 16 Daubechies 16
3 2 2
24–25 16–31 15–32
1 1, 3, 4, 5 1, 2, 3, 4, 5
Fig. 4. Combined frequency responses of the 2 cycle digital mimic filter and: (a) FCDFT filters and (b) HCDFT filters.
Fig. 5 shows the frequency responses of filters FP1, FP2 and FP3 combined with the digital mimic filter. It can be seen that both FP1 and FP2 filters do not eliminate the 2nd harmonic, because it neither was filtered out in wavelet filtering nor was included in the signal estimation model (Eq. (5)). However, FP2 filters are less affected by off nominal frequencies and reject harmonics better than FP1 filters. In fact, FP1 filters do not completely eliminate the 3rd, 6th and 10th harmonics. On the other hand, FP3 filters eliminate all harmonics. The proposed filters have better frequency responses than HCDFT filters, eliminating more harmonics and attenuating better off nominal frequencies. The decaying DC component effect on filters frequency responses can be observed at low frequencies, where it is directly affected by the filter gain. It can be seen that its effect on FP1, FP2
Fig. 5. Combined frequency responses of the 2 cycle digital mimic filter and: (a) FP1 filters (b) FP2 filters and (c) FP3 filters.
Fig. 6. 230 kV power system single line diagram [17].
88
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90 Table 2 Fault data of additional test cases. Case
Fault data Type
Location (miles)
Incidence angle (◦ )
Resistance ()
I II III
AG BG ABG
7.5 30.0 15.0
30 60 45
15 10 20
IV V VI
BCG AB BC
7.5 7.5 15.0
90 30 90
10 1 1
VII VIII IX
CA ABC ABC
22.5 7.5 15.0
60 45 60
5 1 1
and HCDFT filters is quite similar, but FP3 filters are more sensitive to decaying DC component. 5.2. Time response analysis In this section, the focus is on time domain responses analysis of the designed filter pairs. Like frequency response analysis, the time responses of the designed filters are compared to the ones of HCDFT and FCDFT filters. The basic system model shown in Fig. 6 was used to generate fault data to test the performance of the designed filters. It is a 230 kV ATP (Alternative Transient Program [16]) power system reference model proposed by IEEE Power System Relaying Committee. In this model, there are two ideal sources S1 and S3 and a synchronous machine S2, which is connected to bus 4 through a grounded Y- two-winding transformer. The transmission lines consist of one pair of mutually coupled lines (between buses 1 and 2), out of which one is a three terminal line connected to bus 3, and a single circuit line to bus 4. The instrument transformers are also included in this system model. A detailed description of this reference system and the atp files is found in [17]. For the sake of simplicity, the time responses of the designed filters were evaluated considering faults at line 4 of the system model. In fact, these filters were used in a distance relay at the line end connected to bus 2 to estimate voltage and current fundamental phasors. It was considered a mho auto-polarized characteristic, and the 1st and 2nd zone of this relay was set to cover 85% and 125% of the line, respectively. In addition, the relay characteristic angle was set to 60◦ [1]. Fig. 7 shows the magnitude of voltage and current fundamental phasors and the apparent impedance seen by the ZAG unit of the relay for an AG fault (7.5 miles from the bus 2, with incidence angle of 20◦ and fault resistance of 5 ). In respect to time response of phasor computation, it can be seen that both FP1 and FP2 filters have response speed and convergence quite similar to HCDFT filters, even being full cycle filters. On the other hand, FP3 filters have performance in between HCDFT and FCDFT filters, because they are more affected by decaying DC component. In this example, using the FP1 and FP2 filters, the distance relay detects the fault within its 1st zone at 4.69 ms, and using the HCDFT filters at 4.17 ms, whereas using FP3 and FCDFT filters the fault was detected at 5.31 and 9.38 ms, respectively. Table 2 shows the data of additional test cases used to evaluate the performance of distance relay for each filter pair. The results of each additional test case is shown in Table 3. This table compares the fault detection time of the 1st and 2nd zones, considering FP1, FP2, FP3, HCDFT and FCDFT filters. By the analysis of these results, one can see that FP1 and FP2 filters always provide faster fault detection than FCDFT filters in both 1st and 2nd zones, and have response speed and convergence quite similar to HCDFT filters. In
Fig. 7. (a) Voltage phasor magnitude, (b) current phasor magnitude and (c) apparent impedance seen by the ZAG unit of the relay for an AG fault at line 4 (location: 7.5 miles; incidence angle: 20◦ ; resistance: 5 ).
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
89
Table 3 Summary of fault detection time for the additional test cases. Case
Zone
Fault detection time (ms)
I
1st 2nd
6.77 4.69
6.25 4.69
6.25 4.69
5.73 4.17
11.46 9.90
II
1st 2nd
9.38 8.33
10.42 8.33
12.50 10.94
9.38 8.33
17.19 11.46
III
1st 2nd
8.33 6.25
8.33 5.73
10.94 5.21
8.33 5.21
11.46 9.90
IV
1st 2nd
7.29 4.69
6.77 4.69
9.46 4.69
6.25 4.17
11.98 9.90
V
1st 2nd
9.38 6.25
9.38 6.25
11.54 6.25
8.85 5.73
14.06 11.98
VI
1st 2nd
8.85 5.73
9.38 5.73
11.46 5.73
8.85 5.21
13.02 11.46
VII
1st 2nd
10.42 7.81
11.98 7.29
13.58 6.77
9.90 6.77
17.19 12.50
VIII
1st 2nd
7.29 5.21
6.77 5.21
6.77 5.21
6.77 4.69
12.50 10.42
IX
1st 2nd
9.38 8.33
9.90 8.33
12.50 10.94
9.38 7.81
15.63 11.46
FP1
FP2
addition, the FP3 filters provide fault detection times in between the ones provided by HCDFT and FCDFT filters.
FP3
HCDFT
FCDFT
In accordance with the results discussed in this paper, filters designed by the proposed filter design technique show great promise to be used in high-speed distance protection.
6. Conclusions The paper describes a novel MODWT-based technique to design digital filters used to estimate voltage and current fundamental phasors. This technique was used to design three different filter pairs, whose time and frequency responses were compared to the ones of traditional Fourier HCDFT and FCDFT filters. The obtained results show that designed filters have better frequency responses than HCDFT filters, eliminating more harmonics. In addition, although the designed filters are full cycle filters, they always provide faster fault detection than FCDFT filters, and have time responses quite similar to the ones of HCDFT filters, specially the filters FP1 and FP2. This technique may be understood as a generalization of the phasor computation algorithm proposed in [12]. In fact, the FP1 fil-
⎡
h0 ⎢0 B1 = ⎢ . ⎣ .. h1
h1 h0 .. . h2
h2 h1 .. . h3
h3 h2 .. . ···
ters designed here are the particular case of the use of the proposed technique, leading to the phasor computation algorithm reported in [12]. It can be seen that although they provide faster fault detection than FP2 and FP3 filters for some test cases, they have worse frequency responses, do not completely eliminating the 3rd, 6th and 10th harmonics. Comparing FP2 and FP3 filters performances, it is observed that to eliminate the 2nd harmonic the latter suffers more from inaccuracies because of the decaying DC component, delaying the fault detection. Thus, if the 2nd harmonic content is not significant, the FP2 filters may be used instead of FP3 filters, providing faster fault detection.
Appendix A. Mathematical MODWT background The 1st level coefficients of the MODWT pyramid algorithm can be computed by [14]:
W1 V1
=
B1 A1
X,
(A.1)
where X is a N × 1 vector of samples of a discrete signal; W1 and V1 are N × 1 vectors of wavelet and scaling coefficients, respectively; B1 and A1 are N × N matrices. Consider a even-length wavelet filter denoted by {hl : l = 0, 1, . . . , L − 1}, where L is the length of the filter. The matrix B1 takes the form: ··· h3 .. . hL−1
hL−1 ··· .. . 0
0 hL−1 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 ··· 0 ··· .. . . . . 0 ···
⎤
0 0 0 0 ⎥ .. .. ⎥ ⎦ . . 0 h0
(A.2)
It can be seen that the rows of the matrix B1 are circular shifted version of its 1st row, considering a one-sample shift between consecutive rows. The matrix A1 has a similar structure to B1 , with gl s replacing hl s, where {gl : l = 0, 1, . . . , L − 1} is a even-length scaling filter. In fact, the filters {hl √ } and {gl } are rescaled versions of the DWT filters dividing them by 2. In the 2nd MODWT level, V1 is treated in exactly the same manner as X is treated in the 1st level. Thus, the 2nd level wavelet and scaling coefficients can be computed by:
W2 V2
=
B2 A2
V1 =
B2 A1 A2 A1
X,
(A.3)
90
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
where B2 and A2 are N × N matrices. The matrix B2 takes the form:
⎡
h0 ⎢0 B2 = ⎢ . ⎣ .. 0
0 h0 .. . h1
h1 0 .. . 0
h2 0 .. . 0
0 h1 .. . h2
0 h2 .. . h3
h3 0 .. . ···
··· h3 .. . 0
0 ··· .. . hL−1
hL−1 0 .. . 0
0 hL−1 .. . 0
0 0 .. . 0
0 0 .. . 0
One can see that the rows of the matrix B2 are upsampling version of the ones of B1 , putting one zero between the filter coefficients. The matrix A2 has a similar structure to B2 , with gl s replacing hl s. Thus, these filters has width L2 equal to 2L − 1. The j th MODWT level coefficients are computed by:
Wj Vj
=
Bj Aj−1 · · ·A1
Aj Aj−1 · · ·A1
X=
Wj Vj
X,
(A.5)
where Wj = Bj Aj−1 · · ·A1 and Vj = Aj Aj−1 · · ·A1 are N × N matrices. The rows of the matrices Bj and Aj are upsampling version of the ones of matrices B1 and A1 , respectively, putting 2j−1 − 1 zeros between the filter coefficients. Thus, these filters have width Lj equal to 2j−1 (L − 1) + 1. For instance, the 1st row of Bj is: h0 0· · ·0 h1 · · ·hL−2 0· · ·0 hL−1 0· · ·0
2j−1 −1
2j−1 −1
(A.6)
N−Lj
According to MODWT multiresolution analysis, the original signal X can be reconstructed using its smoothed S1 and detailed D1 versions at the 1st MODWT level: X = S1 + D1 = M1 X + N1 X ,
(A.7)
S1
D1
where M1 = VT1 V1 and N1 = WT1 W1 . In a similar way, S1 may be written using the smoothed S2 and detailed D2 versions of the original signal X at the 2nd MODWT level: X = S2 + D2 + D1
S1
= M2 X + N2 X + N1 X ,
(A.8)
S2
D2
D1
VT2 V2
where M2 = and N2 = WT2 W2 . One can see by induction that for a j th MODWT level, the original signal X can be reconstructed by: X = Sj +
j
Dk
k=1
= Mj X +
Sj
j k=1
(A.9) N X ,
Dk
where Mj = VTj Vj and Nj = WTj Wj .
0 ··· 0 ··· .. . . . . 0 0
0 0 .. . ···
⎤
0 0 ⎥ .. ⎥ ⎦ . h0
(A.4)
References [1] G. Ziegler, Numerical Distance Protection: Principles and Applications. Siemens AG, Berlin, 2nd ed., 2006. [2] M.G. Adamiak, G.E. Alexander, W. Premerlani, Advancements in Adaptive Algorithms for Secure High-Speed Distance Protection, in: Western Protective Relay Conf, Spokane, Washington, October, 1996. [3] E.O. Schweitzer, D. Hou III, Filtering for Protective Relays, in: 47th Anual Georgia Tech Protective Relaying Conference, Atlanta, Georgia, 1993. [4] A.T. Johns, S.K. Salman, Digital Protection for Power Systems, IEE Power Series 15, Peter Peregrinus Ltd., United Kingdom, 1995. [5] S. Jiale, S. Guobing, Z. Jiankang, S. Yonghua, A Novel Short-Window Algorithm of Phasor Computation for EHV Line Protection. Electric Power System Research, no. 73, September 2004. [6] B.J. Mann, I.F. Morrison, Digital calculation of impedance for transmission line protection, IEEE Trans Power Appl. Syst. PAS 90 (1) (1971) 270–279. [7] T.S. Sidhu, X. Zhang, V. Balamourougan, A new half-cycle phasor estimation algorithm, IEEE Trans. Power Deliv. 20 (April (2)) (2005) 1299– 1305. [8] J.-C. Gu, S.-L. Yu, Removal of DC Offset in Current and Voltage signals using a novel Fourier algorithm, IEEE Trans. Power Deliv. 15 (1) (2000) 73–79. [9] C.-S. Chen, C.-W. Liu, J.A. Jiang, Application of combined adaptive Fourier filtering techinique and fault detector to fast distance protection, IEEE Trans. Power Deliv. 21 (April (2)) (2006) 619–626. [10] T.S. Sidhu, D.S. Ghotra, M.S. Sachdev, An adaptive distance relay and its performance comparison with a fixed data window distance relay, IEEE Trans. Power Deliv. 17 (July (3)) (2002) 691–697. [11] C.-K. Wong, I.-T. Leong, C.-S. Lei, J.-T. Wu, Y.-D. Han, A novel algorithm for phasor calculation based on wavelet analysis, in: IEEE PES Summer Meeting, Vancouver, BC, Canada, 2001. [12] F. Liang, B. Jeyasurya, Transmission line distance protection using wavelet transform algorithm, IEEE Trans. Power Deliv. 19 (April (2)) (2004) 545– 553. [13] A.H. Osman, O.P. Malik, Transmission line distance protection based on wavelet transform, IEEE Trans. Power Deliv. 19 (April (2)) (2004) 515–523. [14] D.B. Percival, A.T. Walden, Wavelet Methods for Time Series Analysis, Cambridge University Press, New York, 2000. [15] G. Benmouyal, Removal of DC-offset in current waveforms using digital mimic filtering, IEEE Trans. Power Deliv. 10 (April (2)) (1995) 621–630. [16] CanAm EMTP User Group, Alternative Transient Program (ATP) Rule Book, Portland, USA, 2001. [17] IEEE PES/PSRC Special Publication, EMTP Reference Models for Transmission Line Relay Testing, 2004. Available: http://www.pespsrc.org/Reports/D10/ Report rev0.zip.
Contents lists available at ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
Distance protection using a wavelet-based filtering algorithm Kleber M. Silva a,∗ , Washington L.A. Neves b , Benemar A. Souza b a b
University of Brasília, Electrical Engineering Department, 70919-970 Brasília, Brazil Federal University of Campina Grande, Paraiba, Brazil, Electrical Engineering Department, 882 Aprigio Veloso Ave, Bodocongo, 58109-970 Campina Grande, Brazil
a r t i c l e
i n f o
Article history: Received 2 December 2008 Received in revised form 1 July 2009 Accepted 21 August 2009 Available online 27 September 2009 Keywords: High-speed distance protection Wavelet transform Phasor estimation
a b s t r a c t This paper presents a novel filter design technique based on the maximal overlap discrete wavelet transform (MODWT). Time and frequency response of the designed filters are compared to the ones of traditional discrete Fourier transform (DFT) based filters. The obtained results show improvements on response speed when comparing to full cycle traditional filters and improvements on frequency response when compared to half cycle traditional filters. According to these characteristics, the designed filter may be used in secure high-speed distance protection, since it provides fast relay operating times, without compromising their frequency response. © 2009 Elsevier B.V. All rights reserved.
1. Introduction Due to the growth of the complexity and interconnectivity of the HV and EHV power systems, modern protection schemes must be designed to clear faults rapidly in order to maintain the system stability and enhance the system security. Thus, it is necessary to detect a fault within sub-cycles of the power system frequency [1]. Distance protection relays are widely used to protect transmission lines. Their operation depends on both voltage and current phasors at the fundamental frequency. In fact, the accuracy, speed and security of these relays are directly related to the operation of the digital filtering algorithms used to estimate phasors [2], which must have certain characteristics, such as: bandpass response about the system frequency, decaying DC rejection, harmonic attenuation or rejection and good transient behavior [3]. The most popular phasor estimation algorithms are the full cycle DFT (FCDFT) and the full cycle least error square (LES) algorithms [4]. However, they are not suitable for high-speed distance protection, because relays are required to operate in less than one cycle. In this way, many short window algorithms have been reported [5,6]. Although these algorithms are very fast, they are very inaccurate because they do not eliminate harmonics. Another important fast algorithm is the half cycle DFT (HCDFT) [4], but it also suffers from inaccuracies because it neither eliminates even harmonics nor rejects decaying DC components. Many improved half cycle algorithms have been reported [7,8], but they do not eliminate even harmonics. To overcome these drawbacks, some adaptive algorithms have been developed, which varies the
window length after detecting the fault [9,10]. However, they are slower than the HCDFT algorithm. Alternatively, the discrete wavelet transform (DWT) has been also used to develop innovative phasor estimation algorithms [11,12]. In [11], the authors use the orthogonality property of the DWT to compute synchronized phasors. This method was used in [13] to propose a new DWT-based transmission line distance protection algorithm. However, this algorithm cannot be used in high-speed distance protection, since it does not provide a fast relay operating time and is heavily affected by harmonic and decaying DC components. An algorithm based on the redundant version of the DWT, named MODWT, was proposed in [12]. It uses the good signal approximation characteristic provided by the MODWT, computing the phasor by means of a two-samples algorithm applied to the smooth version of the original voltage and current signals. The reported results show that this algorithm is faster than FCDFT. Unfortunately, the authors did not analyze the frequency response of the overall filtering system. This paper presents a novel filter design technique based on the MODWT. This technique may be understood as a generalization of the phasor computation algorithm proposed in [12]. Time and frequency responses of proposed filters are compared to the ones of traditional DFT-based filters. The obtained results shown that the proposed and HCDFT filters have similar time responses, with the advantage that the frequency responses of the proposed filters are close to the ones of FCDFT filters. It indicates that the filters developed here are suitable for secure high-speed distance protection.
2. MODWT fundamentals ∗ Corresponding author. Tel.: +55 61 3307 2308; fax: +55 61 3274 6651. E-mail address: [email protected] (K.M. Silva). 0378-7796/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2009.08.012
The MODWT is similar to the DWT in that both produce a set of time-dependent wavelet and scaling coefficients, but the MODWT
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
85
is a highly redundant nonorthogonal transform, retaining downsampled values at each level of the decomposition that would be otherwise discarded by the DWT. In addition, the DWT cause more distortions in signal reconstruction and requires that the total number of samples must be a power of two [12,14]. Using wavelet transform a signal is decomposed into smoothed and detailed versions for different levels of resolution. For instance, a discrete signal X may be written using its smoothed S1 and detailed D1 versions at the 1st MODWT level [14]: X = S1 + D1
(1)
and using its smoothed S2 and detailed D2 versions at the 2nd MODWT level: X = S2 + D2 + D1
(2)
and so on. In this decomposition the signals Sj and Dj are computed as: Sj = Mj X
and
Dj = Nj X,
(3)
where the matrices Mj and Nj are obtained for each mother wavelet as shown in Appendix A. This analysis may also be understood as a filtering process, in which each MODWT level of decomposition divides the bandwidth of X into octave bands, as shown in Fig. 1 that depicts the frequency response of the 16th rows of the matrices N1 , N2 , N3 , M1 , M2 and M3 computed using the Daubechies 8 mother wavelet [14], for a generic sampling frequency fs . The rows of the matrix Nj work as low-pass filters, whereas the rows of the matrix Mj work as bandpass filters. Since the focus are the phasors at the fundamental frequency, only matrices Mj are analyzed. A discrete signal X is considered to illustrate the MODWT multiresolution analysis: xk = sin ωtk +
1 1 sin 3ωtk + sin 5ωtk , 3 5
Fig. 2. Original signal x and its smoothed version S3 at 3rd MODWT level, obtained using the mother wavelet Daubechies 8.
Fig. 3. Frequency response of the 16th, 24th and 32nd rows of M3 using Daubechies 8.
(4)
where xk is the k th sample of the signal X, ω = 2f is the fundamental angular frequency, f the fundamental frequency and tk is the k th sampling interval. Fig. 2 illustrates one cycle of the original signal X and its smoothed version S3 , computed using the mother wavelet Daubechies 8 and assuming the sampling rate of 32 samples per cycle. It is observed that the smoothed version S3 contains only the fundamental frequency component. In fact, the 3rd and 5th harmonics of X are eliminated by the rows of M3 , as shown in Fig. 3, which depicts the frequency response of the 16th, 24th and 32nd rows of M3 . In the same way, the rows of the matrices M2 and M1 eliminate some harmonics too.
From the above discussion, it is clear that the rows of the matrices Mj eliminate some harmonics. This filtering characteristic is quite suitable for phasor estimation and is taken into account in the proposed filter design technique. Further details are discussed next. 3. Proposed filter design technique In order to present the proposed filter design technique, consider one cycle of a sinusoidal signal X, which has N samples within a cycle. Its smoothed version Sj at the j th MODWT level may be computed by Eq. (3), in which the order of the matrix Mj is N × N. Thus, the signal Sj have also N samples, which may be written as: sj,k =
H
(Ych sin hωtk + Ysh cos hωtk ) + εk
(5)
h=1
where sj,k is the k th sample of the signal Sj ; H is the total of specified harmonics; Ych and Ysh are, respectively, the real and imaginary parts of the h th harmonic and εk is the estimation error of the k th sample. The Eq. (5) may be applied to M consecutive samples of the signal Sj , with M ≤ N, resulting in a linear system of equations: AY + E = B,
(6)
where A is a M × 2H matrix:
⎡
Fig. 1. Frequency response of the 16th row of the matrices N1 , N2 , N3 , M1 , M2 and M3 computed to Daubechies 8.
sin ωt0 ⎢ sin ωt1 A = ⎢. ⎣ .. sin ωtM−1
cos ωt0 cos ωt1 .. . cos ωtM−1
··· ··· .. . ···
sin Hωt0 sin Hωt1 .. . sin HωtM−1
⎤
cos Hωt0 cos Hωt1 ⎥ ⎥ .. ⎦ . cos HωtM−1
86
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
Y is a 2H × 1 vector of variables:
4.3. On choosing the Mj matrix
Y = Yc1
and E is a M × 1 vector of estimation errors:
Once the mother wavelet and the sampling rate are defined, the total number of MODWT levels is computed for one cycle of the signal X (Appendix A). The Mj matrix to be used is the last one. For example, considering a sampling rate of 32 samples per cycle and the mother wavelet Daubechies 16, it is possible to compute two MODWT levels. Thus, M2 is the chosen matrix.
E = ε0
4.4. Harmonic filtering
Ys1
· · · YcH
YsH
T
,
B is a vector with M consecutive samples of the signal Sj :
B = sj,0
sj,1
ε1
sj,2
ε2
· · · sj,M−2
· · · εM−2
sj,M−1
εM−1
T
T
.
The solution of Eq. (6) that minimizes the sum ET E, is obtained using the left pseudo-inverse A+ matrix [4]: Y = A+ B,
(7)
As discussed before, the rows of the matrices Mj eliminate some harmonics. For example, the 16th row of the matrix M2 in Fig. 1 eliminates all harmonics beyond the 6th harmonic, thus it is necessary to consider in the model only the 1st up to the 5th harmonics. The frequency response of the Mj rows defines which harmonics must be considered in the signal estimation model.
(8)
4.5. On choosing the M rows of the matrix Mj
−1
where A+ = AT A AT . From Eqs. (3) and (7), one can obtain:
Y = A+ Mj
M
X = GX,
G
where [Mj ]M is formed by M consecutive rows of the matrix Mj , in such way that G is a 2H × N matrix. The matrix G in Eq. (8) is the base of the proposed filter design technique. It is computed in advance in an off-line mode. The real part hc and imaginary part hs filters proposed to estimate the fundamental frequency phasors are the 1st and 2nd rows of G, respectively. The frequency response of the Mj rows defines the designed filters. It should be noticed that these filters are full cycle filters. Experience has shown that the more samples within one cycle are used, the more harmonics are filtered out, that is the reason why all samples within one cycle were used. The proposed technique may be understood as a generalization of the phasor computation algorithm proposed in [12]. This algorithm computes fundamental phasors using equations, formulated specifically to two consecutive samples of the smoothed signals at the 3rd MOWDT level, using the mother wavelet Daubechies 8. In fact, one can see that these equations are a particular case of Eq. (8). 4. Filter design considerations This section discusses general design considerations that must be taken into account to design filters that present good time and frequency responses. 4.1. Sampling rate As it can be seen in Figs. 1 and 3, the frequency responses of the rows of the Mj matrix have zeros. An appropriate choice of the sampling rate can put those zeros to work to notch out the harmonics. This can be achieved by choosing a sampling rate that is multiple of the fundamental frequency.
It is desirable that the chosen consecutive M rows of matrix Mj lead to designed filters hc and hs with good frequency response and small time delays. For instance, in Fig. 3, the 32nd row of M3 leads to the smallest time delay and to the poorest frequency response. On the other hand, the 16th row of M3 leads to a smoothest frequency response but to a larger time delay. In this case, the 24th row would be chosen, since it leads to good frequency response and a small time delay. M is directly related to the total number of harmonics to be reject in the signal estimation model, in such way that the more harmonics to be reject, the more consecutive rows must be taken into account. 4.6. Decaying DC component Traditionally, in LES algorithm applications, the decaying DC component is estimated by including the first terms of its Taylor series expansion in Eq. (5). However, in the cases discussed here, it leads to time delay in phasor computation. In this paper, the traditional digital mimic filter, which is a first order high-pass filter with unitary gain at the fundamental frequency, is used to eliminate the decaying DC component [15]. The final characteristics of the proposed filters are the combination of the one of the mimic filter with those of hc and hs filters. 5. Designed filters evaluation This paper presents the analysis of three filter pairs designed using the proposed technique. Table 1 shows the chosen design parameters: the mother wavelet, the MODWT level, the rows of Mj and the harmonics included in Eq. (5). In special, the filter pair FP1 is a particular case leading to the phasor computation algorithm proposed in [12]. The frequency and time responses of these filters were evaluated and the obtained results are discussed next. 5.1. Frequency response analysis
4.2. Mother wavelet The proposed filter design technique can be applied to any mother wavelet. In fact, different mother wavelets (Daubechies, Symlets, Coiflets, etc., [14]) with the same order leads to designed filters with nearly the same frequency responses. For example, filters designed using mother wavelets Daubechies 8 and Symlet 8 are quite similar, but those designed using mother wavelets Daubechies 8 and 4 are quite different. In addition, it was observed that the larger the order of the mother wavelet is, the better frequency responses the designed filters hc and hs have.
In this section, the frequency response of designed filter pairs are compared to the ones of traditional Fourier HCDFT and FCDFT filters. In order to do that, the designed FP1, FP2 and FP3 and Fourier filters use the same 2 cycle digital mimic filter [15]. The sampling rate is 32 samples per cycle. The frequency responses of FCDFT and HCDFT filters combined with mimic filter are shown in Fig. 4[4]. It is well known that FCDFT filters eliminate all harmonics, whereas HCDFT filters do not eliminate even harmonics. In addition, HCDFT filters are more affected by decaying DC and off nominal frequency components [3].
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
87
Table 1 Parameters used to design the evaluated filter pairs. Filter pair
Mother wavelet
MODWT level j
Rows of Mj
Harmonics in Eq. (5)
FP1 FP2 FP3
Daubechies 8 Daubechies 16 Daubechies 16
3 2 2
24–25 16–31 15–32
1 1, 3, 4, 5 1, 2, 3, 4, 5
Fig. 4. Combined frequency responses of the 2 cycle digital mimic filter and: (a) FCDFT filters and (b) HCDFT filters.
Fig. 5 shows the frequency responses of filters FP1, FP2 and FP3 combined with the digital mimic filter. It can be seen that both FP1 and FP2 filters do not eliminate the 2nd harmonic, because it neither was filtered out in wavelet filtering nor was included in the signal estimation model (Eq. (5)). However, FP2 filters are less affected by off nominal frequencies and reject harmonics better than FP1 filters. In fact, FP1 filters do not completely eliminate the 3rd, 6th and 10th harmonics. On the other hand, FP3 filters eliminate all harmonics. The proposed filters have better frequency responses than HCDFT filters, eliminating more harmonics and attenuating better off nominal frequencies. The decaying DC component effect on filters frequency responses can be observed at low frequencies, where it is directly affected by the filter gain. It can be seen that its effect on FP1, FP2
Fig. 5. Combined frequency responses of the 2 cycle digital mimic filter and: (a) FP1 filters (b) FP2 filters and (c) FP3 filters.
Fig. 6. 230 kV power system single line diagram [17].
88
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90 Table 2 Fault data of additional test cases. Case
Fault data Type
Location (miles)
Incidence angle (◦ )
Resistance ()
I II III
AG BG ABG
7.5 30.0 15.0
30 60 45
15 10 20
IV V VI
BCG AB BC
7.5 7.5 15.0
90 30 90
10 1 1
VII VIII IX
CA ABC ABC
22.5 7.5 15.0
60 45 60
5 1 1
and HCDFT filters is quite similar, but FP3 filters are more sensitive to decaying DC component. 5.2. Time response analysis In this section, the focus is on time domain responses analysis of the designed filter pairs. Like frequency response analysis, the time responses of the designed filters are compared to the ones of HCDFT and FCDFT filters. The basic system model shown in Fig. 6 was used to generate fault data to test the performance of the designed filters. It is a 230 kV ATP (Alternative Transient Program [16]) power system reference model proposed by IEEE Power System Relaying Committee. In this model, there are two ideal sources S1 and S3 and a synchronous machine S2, which is connected to bus 4 through a grounded Y- two-winding transformer. The transmission lines consist of one pair of mutually coupled lines (between buses 1 and 2), out of which one is a three terminal line connected to bus 3, and a single circuit line to bus 4. The instrument transformers are also included in this system model. A detailed description of this reference system and the atp files is found in [17]. For the sake of simplicity, the time responses of the designed filters were evaluated considering faults at line 4 of the system model. In fact, these filters were used in a distance relay at the line end connected to bus 2 to estimate voltage and current fundamental phasors. It was considered a mho auto-polarized characteristic, and the 1st and 2nd zone of this relay was set to cover 85% and 125% of the line, respectively. In addition, the relay characteristic angle was set to 60◦ [1]. Fig. 7 shows the magnitude of voltage and current fundamental phasors and the apparent impedance seen by the ZAG unit of the relay for an AG fault (7.5 miles from the bus 2, with incidence angle of 20◦ and fault resistance of 5 ). In respect to time response of phasor computation, it can be seen that both FP1 and FP2 filters have response speed and convergence quite similar to HCDFT filters, even being full cycle filters. On the other hand, FP3 filters have performance in between HCDFT and FCDFT filters, because they are more affected by decaying DC component. In this example, using the FP1 and FP2 filters, the distance relay detects the fault within its 1st zone at 4.69 ms, and using the HCDFT filters at 4.17 ms, whereas using FP3 and FCDFT filters the fault was detected at 5.31 and 9.38 ms, respectively. Table 2 shows the data of additional test cases used to evaluate the performance of distance relay for each filter pair. The results of each additional test case is shown in Table 3. This table compares the fault detection time of the 1st and 2nd zones, considering FP1, FP2, FP3, HCDFT and FCDFT filters. By the analysis of these results, one can see that FP1 and FP2 filters always provide faster fault detection than FCDFT filters in both 1st and 2nd zones, and have response speed and convergence quite similar to HCDFT filters. In
Fig. 7. (a) Voltage phasor magnitude, (b) current phasor magnitude and (c) apparent impedance seen by the ZAG unit of the relay for an AG fault at line 4 (location: 7.5 miles; incidence angle: 20◦ ; resistance: 5 ).
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
89
Table 3 Summary of fault detection time for the additional test cases. Case
Zone
Fault detection time (ms)
I
1st 2nd
6.77 4.69
6.25 4.69
6.25 4.69
5.73 4.17
11.46 9.90
II
1st 2nd
9.38 8.33
10.42 8.33
12.50 10.94
9.38 8.33
17.19 11.46
III
1st 2nd
8.33 6.25
8.33 5.73
10.94 5.21
8.33 5.21
11.46 9.90
IV
1st 2nd
7.29 4.69
6.77 4.69
9.46 4.69
6.25 4.17
11.98 9.90
V
1st 2nd
9.38 6.25
9.38 6.25
11.54 6.25
8.85 5.73
14.06 11.98
VI
1st 2nd
8.85 5.73
9.38 5.73
11.46 5.73
8.85 5.21
13.02 11.46
VII
1st 2nd
10.42 7.81
11.98 7.29
13.58 6.77
9.90 6.77
17.19 12.50
VIII
1st 2nd
7.29 5.21
6.77 5.21
6.77 5.21
6.77 4.69
12.50 10.42
IX
1st 2nd
9.38 8.33
9.90 8.33
12.50 10.94
9.38 7.81
15.63 11.46
FP1
FP2
addition, the FP3 filters provide fault detection times in between the ones provided by HCDFT and FCDFT filters.
FP3
HCDFT
FCDFT
In accordance with the results discussed in this paper, filters designed by the proposed filter design technique show great promise to be used in high-speed distance protection.
6. Conclusions The paper describes a novel MODWT-based technique to design digital filters used to estimate voltage and current fundamental phasors. This technique was used to design three different filter pairs, whose time and frequency responses were compared to the ones of traditional Fourier HCDFT and FCDFT filters. The obtained results show that designed filters have better frequency responses than HCDFT filters, eliminating more harmonics. In addition, although the designed filters are full cycle filters, they always provide faster fault detection than FCDFT filters, and have time responses quite similar to the ones of HCDFT filters, specially the filters FP1 and FP2. This technique may be understood as a generalization of the phasor computation algorithm proposed in [12]. In fact, the FP1 fil-
⎡
h0 ⎢0 B1 = ⎢ . ⎣ .. h1
h1 h0 .. . h2
h2 h1 .. . h3
h3 h2 .. . ···
ters designed here are the particular case of the use of the proposed technique, leading to the phasor computation algorithm reported in [12]. It can be seen that although they provide faster fault detection than FP2 and FP3 filters for some test cases, they have worse frequency responses, do not completely eliminating the 3rd, 6th and 10th harmonics. Comparing FP2 and FP3 filters performances, it is observed that to eliminate the 2nd harmonic the latter suffers more from inaccuracies because of the decaying DC component, delaying the fault detection. Thus, if the 2nd harmonic content is not significant, the FP2 filters may be used instead of FP3 filters, providing faster fault detection.
Appendix A. Mathematical MODWT background The 1st level coefficients of the MODWT pyramid algorithm can be computed by [14]:
W1 V1
=
B1 A1
X,
(A.1)
where X is a N × 1 vector of samples of a discrete signal; W1 and V1 are N × 1 vectors of wavelet and scaling coefficients, respectively; B1 and A1 are N × N matrices. Consider a even-length wavelet filter denoted by {hl : l = 0, 1, . . . , L − 1}, where L is the length of the filter. The matrix B1 takes the form: ··· h3 .. . hL−1
hL−1 ··· .. . 0
0 hL−1 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 0 .. . 0
0 ··· 0 ··· .. . . . . 0 ···
⎤
0 0 0 0 ⎥ .. .. ⎥ ⎦ . . 0 h0
(A.2)
It can be seen that the rows of the matrix B1 are circular shifted version of its 1st row, considering a one-sample shift between consecutive rows. The matrix A1 has a similar structure to B1 , with gl s replacing hl s, where {gl : l = 0, 1, . . . , L − 1} is a even-length scaling filter. In fact, the filters {hl √ } and {gl } are rescaled versions of the DWT filters dividing them by 2. In the 2nd MODWT level, V1 is treated in exactly the same manner as X is treated in the 1st level. Thus, the 2nd level wavelet and scaling coefficients can be computed by:
W2 V2
=
B2 A2
V1 =
B2 A1 A2 A1
X,
(A.3)
90
K.M. Silva et al. / Electric Power Systems Research 80 (2010) 84–90
where B2 and A2 are N × N matrices. The matrix B2 takes the form:
⎡
h0 ⎢0 B2 = ⎢ . ⎣ .. 0
0 h0 .. . h1
h1 0 .. . 0
h2 0 .. . 0
0 h1 .. . h2
0 h2 .. . h3
h3 0 .. . ···
··· h3 .. . 0
0 ··· .. . hL−1
hL−1 0 .. . 0
0 hL−1 .. . 0
0 0 .. . 0
0 0 .. . 0
One can see that the rows of the matrix B2 are upsampling version of the ones of B1 , putting one zero between the filter coefficients. The matrix A2 has a similar structure to B2 , with gl s replacing hl s. Thus, these filters has width L2 equal to 2L − 1. The j th MODWT level coefficients are computed by:
Wj Vj
=
Bj Aj−1 · · ·A1
Aj Aj−1 · · ·A1
X=
Wj Vj
X,
(A.5)
where Wj = Bj Aj−1 · · ·A1 and Vj = Aj Aj−1 · · ·A1 are N × N matrices. The rows of the matrices Bj and Aj are upsampling version of the ones of matrices B1 and A1 , respectively, putting 2j−1 − 1 zeros between the filter coefficients. Thus, these filters have width Lj equal to 2j−1 (L − 1) + 1. For instance, the 1st row of Bj is: h0 0· · ·0 h1 · · ·hL−2 0· · ·0 hL−1 0· · ·0
2j−1 −1
2j−1 −1
(A.6)
N−Lj
According to MODWT multiresolution analysis, the original signal X can be reconstructed using its smoothed S1 and detailed D1 versions at the 1st MODWT level: X = S1 + D1 = M1 X + N1 X ,
(A.7)
S1
D1
where M1 = VT1 V1 and N1 = WT1 W1 . In a similar way, S1 may be written using the smoothed S2 and detailed D2 versions of the original signal X at the 2nd MODWT level: X = S2 + D2 + D1
S1
= M2 X + N2 X + N1 X ,
(A.8)
S2
D2
D1
VT2 V2
where M2 = and N2 = WT2 W2 . One can see by induction that for a j th MODWT level, the original signal X can be reconstructed by: X = Sj +
j
Dk
k=1
= Mj X +
Sj
j k=1
(A.9) N X ,
Dk
where Mj = VTj Vj and Nj = WTj Wj .
0 ··· 0 ··· .. . . . . 0 0
0 0 .. . ···
⎤
0 0 ⎥ .. ⎥ ⎦ . h0
(A.4)
References [1] G. Ziegler, Numerical Distance Protection: Principles and Applications. Siemens AG, Berlin, 2nd ed., 2006. [2] M.G. Adamiak, G.E. Alexander, W. Premerlani, Advancements in Adaptive Algorithms for Secure High-Speed Distance Protection, in: Western Protective Relay Conf, Spokane, Washington, October, 1996. [3] E.O. Schweitzer, D. Hou III, Filtering for Protective Relays, in: 47th Anual Georgia Tech Protective Relaying Conference, Atlanta, Georgia, 1993. [4] A.T. Johns, S.K. Salman, Digital Protection for Power Systems, IEE Power Series 15, Peter Peregrinus Ltd., United Kingdom, 1995. [5] S. Jiale, S. Guobing, Z. Jiankang, S. Yonghua, A Novel Short-Window Algorithm of Phasor Computation for EHV Line Protection. Electric Power System Research, no. 73, September 2004. [6] B.J. Mann, I.F. Morrison, Digital calculation of impedance for transmission line protection, IEEE Trans Power Appl. Syst. PAS 90 (1) (1971) 270–279. [7] T.S. Sidhu, X. Zhang, V. Balamourougan, A new half-cycle phasor estimation algorithm, IEEE Trans. Power Deliv. 20 (April (2)) (2005) 1299– 1305. [8] J.-C. Gu, S.-L. Yu, Removal of DC Offset in Current and Voltage signals using a novel Fourier algorithm, IEEE Trans. Power Deliv. 15 (1) (2000) 73–79. [9] C.-S. Chen, C.-W. Liu, J.A. Jiang, Application of combined adaptive Fourier filtering techinique and fault detector to fast distance protection, IEEE Trans. Power Deliv. 21 (April (2)) (2006) 619–626. [10] T.S. Sidhu, D.S. Ghotra, M.S. Sachdev, An adaptive distance relay and its performance comparison with a fixed data window distance relay, IEEE Trans. Power Deliv. 17 (July (3)) (2002) 691–697. [11] C.-K. Wong, I.-T. Leong, C.-S. Lei, J.-T. Wu, Y.-D. Han, A novel algorithm for phasor calculation based on wavelet analysis, in: IEEE PES Summer Meeting, Vancouver, BC, Canada, 2001. [12] F. Liang, B. Jeyasurya, Transmission line distance protection using wavelet transform algorithm, IEEE Trans. Power Deliv. 19 (April (2)) (2004) 545– 553. [13] A.H. Osman, O.P. Malik, Transmission line distance protection based on wavelet transform, IEEE Trans. Power Deliv. 19 (April (2)) (2004) 515–523. [14] D.B. Percival, A.T. Walden, Wavelet Methods for Time Series Analysis, Cambridge University Press, New York, 2000. [15] G. Benmouyal, Removal of DC-offset in current waveforms using digital mimic filtering, IEEE Trans. Power Deliv. 10 (April (2)) (1995) 621–630. [16] CanAm EMTP User Group, Alternative Transient Program (ATP) Rule Book, Portland, USA, 2001. [17] IEEE PES/PSRC Special Publication, EMTP Reference Models for Transmission Line Relay Testing, 2004. Available: http://www.pespsrc.org/Reports/D10/ Report rev0.zip.