Selective estimation of harmonic components in noisy electrical signals for protective relaying purposes

Electrical Power and Energy Systems 56 (2014) 140–146

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Selective estimation of harmonic components in noisy electrical signals for protective relaying purposes J. Lázaro a,⇑, J.F. Miñambres b, M.A. Zorrozua b a b

Department of Applied Mathematics, Faculty of Engineering (Bilbao), University of the Basque Country UPV/EHU, Spain Department of Electrical Engineering, Faculty of Engineering (Bilbao), University of the Basque Country UPV/EHU, Spain

a r t i c l e

i n f o

Article history: Received 22 October 2012 Received in revised form 17 October 2013 Accepted 6 November 2013

Keywords: Digital filtering Digital relaying Noisy signal Phasor estimation algorithm

a b s t r a c t Digital filters used in protection relays must comply with a series of requirements. Among its most important features are: low computational load, good behaviour in the presence of harmonic and decaying dc components, accuracy against noise in the input signal and good behaviour in the transition from steady state to fault period. Before implementing a digital filter in a protection relay, it is necessary to subject it to stringent tests simulating a large array of real fault conditions. Electrical signals undergo a series of transformations before reaching the digital filter. As a result of this, the analysed signal carries errors and noise with it. The filter will only be apt for use in a digital relay if it rapidly and precisely responds to these types of signals. The objective of this paper is to present a new digital filter for selective estimation of harmonic components in noisy electrical signals for protective relaying purposes. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction All power systems must be provided with a protection system that effectively contributes to the reliability of electrical energy supply. Different technologies have been developed throughout time with the aim of providing protection relays with the capacity to carry out two basic tasks. On the one hand, they must be able to detect any faults occurring in the power system, while on the other hand, they must be able to give the necessary orders to isolate the faulty circuit and avoid the fault propagation. Both tasks must be carried out in as short a time as possible to minimise the impact of the fault on the power system. Digital technologies along with the need to optimise resources have currently led to the development of interesting additional features, such as fault location and data recording. Due to this, digital relays have become a fundamental tool for protecting power systems, thanks to their multi-functional features. In a digital relay (Fig. 1), signal processing can be divided schematically into the three following stages:  Measurement and preconditioning of signals: The analogical waves, from the Current Transformer (CT) and Voltage Transformer (VT), pass through an antialiasing (low pass) filter that eliminates the

⇑ Corresponding author. Address: Faculty of Engineering, Alda. Urquijo, s/n, 48013 Bilbao, Spain. Tel.: +34 946014150. E-mail address: [email protected] (J. Lázaro). 0142-0615/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ijepes.2013.11.007

high frequency components. Then, an Analogic Digital Converter (ADC) provides the discrete samples that will be numerically processed in the following step.  Digital filtering of the sampled signals: This process is carried out by a Digital Signal Processor (DSP). The digital filter implemented in the DSP estimates the harmonic components (in most cases the fundamental component) that will be used by the protection functions.  Protection functions: A set of logical and mathematical criteria which, based on the output data of the digital filter, determine the correct decision to be taken by the relay. The protection relay must take the correct decision in the shortest possible time. Therefore, two fundamental features of the digital relay are precision and speed in its response, especially during the transient period of the input signal. The key element to achieve both features is the digital filter. Typical algorithms used in digital filters to carry out phasor estimation are based on the Discrete Fourier Transform (DFT) [1,2]. But signals containing aperiodic components produce significant oscillations in the output of the DFT. Due to the fact that decaying dc-offset components are very common in transient periods, several authors have proposed modifications in the DFT to improve its convergence [3–7]. These methods are based on estimating and eliminating the exponential component and then applying the DFT to the resulting signal. One important drawback of these methods is that they require knowing the fault instant. As a consequence, an additional algorithm is needed to detect the

141

J. Lázaro et al. / Electrical Power and Energy Systems 56 (2014) 140–146

Fig. 1. Diagram showing signal processing in a digital relay.

fault instant. By this reason, these methods have the disadvantage of significantly increasing the computational load [8,9]. Other methods for phasor estimation are based on the Least Error Square (LES) technique [10–13]. In these methods, a system of equations is defined according to the components considered a priori in the analysed signal. The greater the number of components considered, the greater is the precision but also the greater is the computational load. Different solutions have been proposed to limit this computational load. Some authors propose the use of different algorithms according to the status (steady state/transient) of the input signal. However, this option requires knowing the fault instant [10], which has the disadvantage of having to implement an algorithm for fault detection. Other authors propose an application of the LES method to a limited set of harmonic components [11], which means losing precision in the results when the input signal is more complex. This paper proposes a new digital filter that can be implemented in digital relays. Its basic characteristics are the following:  Precision: The accuracy provided in the output enables the use of the proposed filter in digital relaying.  Low computational load: The mathematical techniques involved are well below the limits imposed by the DSP hardware for real time applications.  Independence of the fault instant: The proposed digital filter does not require knowing the fault instant, as it is applied in the same manner in all possible conditions (steady state/transient, pre-fault/fault, all fault types). In addition, it is also necessary to guarantee the correct performance of the filter in the presence of distortions due to errors of discretization and electromagnetic noise. In the case of the ADC [14,15], the errors depend on the maximum desired range of measurement and the number of bits of the ADC. Moreover, there are additional errors such as: offset error of the analogical circuits, measurement error of the internal transformers of the relay, error of derivation by temperature, electromagnetic noise and others. These errors of internal origin are random and difficult to predict. Therefore, the existence of distortions in the signals processed is inherent to the operation of all the digital relays. Consequently, the response of a digital filter must be both fast and accurate in the presence of noisy signals and decaying dc-offset in order to be successfully used for protective relaying purposes.

2. The new filtering technique: s-CharmDF definition The proposed method is a CharmDF-based modified digital filter. In comparison to the CharmDF [16,17], the proposed filter improves precision of results due to its flexibility of application. This flexibility enables it to optimise the estimation of any harmonic component. To illustrate the definition of this new digital filter a generic signal y(t) is used. This signal has the typical characteristics of electrical signals incoming to digital relays. According to the existing technical literature on the subject [18–21], the following signal is considered:

yðtÞ ¼ C 0 þ Cet=s þ

Xn r¼1

Ar cosðxr t þ ar Þ þ zðtÞ

ð1Þ

where C0 is a constant offset, C and s are the amplitude and time constant of decaying dc offset, Ar and ar are the amplitude and phase of the fundamental (r = 1) and harmonic components and, finally, z(t) is the noise component of the signal. The proposed methodology is based on defining an auxiliary wave xs(t) that contains all the information regarding the periodic components of y(t). The auxiliary signal xs(t) has been denominated s-Characteristic HARMonic wave (s-Charm wave) due to the fact that it is defined using samples of y(t) separated by a sample slip of s samples.

xs ðtÞ ¼ fs ðtÞ þ

Xn r¼1

X r cosðxr t þ br Þ

ð2Þ

where fs(t) is a residual signal that includes the influence of aperiodic component and noise. Moreover, Xr and br are the amplitude and phase of the fundamental (r = 1) and harmonic components of the s-Charm wave. The definition of the s-Charm wave implies a biunivocal relationship between the harmonic spectra of signals xs ðtÞ and yðtÞ. For an original signal recorded with a sampling rate of N samples per cycle, the relationship between the respective amplitudes and phases can be expressed through the following correction factors:

Ksr ¼

Ar 1   ¼ X r 2sin 180sr N

  N  2sr hsr ¼ bsr  ar ¼ 90 N

ð3Þ

ð4Þ

142

J. Lázaro et al. / Electrical Power and Energy Systems 56 (2014) 140–146

The sample slip s allows to minimise the influence of residual signal fs(t) in the s-Charm wave xs ðtÞ. By this reason, the application of a phasor estimation algorithm to xs ðtÞ provides more accurate results than its application to yðtÞ. The results shown in this paper have been obtained using FCDFT. It is very important to take into account that the correction factors (3) and (4) only depends of N, s and r, and are independent of the characteristics of the original signal being analysed. Also, it is important to notice that Ksr and hsr are independent of time and so they have the same value in pre-fault, fault and post-fault periods.

Table 1 Slip ranges up to the theoretical optimum (sot).

r=1 . . .2 . . .3 . . .4 . . .5 . . .6 . . .7 . . .8

N = 16

N = 32

N = 64

1– 1– 1– 1– 1– (1) (1) (1)

1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1

(8) (4) (3) (2) (2)

– – – – – – – –

(16) (8) (5) (4) (3) (3) (2) (2)

– – – – – – – –

(32) (16) (11) (8) (6) (5) (5) (4)

3. Optimisation of the new method To define the s-Charm-wave it is necessary to select previously the value of the sample slip s. This section describes the optimisation of the method s-CharmDF for selective estimation of the r-th harmonic component in noisy fault signals. The goal is to determine an optimum sample slip, sot, which will minimise the impact of error in the response of the phasor estimation algorithm when estimating the r-th harmonic. To achieve this objective, it is necessary to reach the maximum possible amplitude for the r-th harmonic component in xs ðtÞ. Thus, from the Eq. (3):

Xr ¼

Ar ¼ 2 Ar sinð180 s r=NÞ Ksr

ð5Þ

The value of s for which the r-th harmonic of the s-Charm wave reaches its maximum amplitude is denominated theoretical optimal sample slip (sot). From (5), this value is obtained when sin(180sotr/N) = 1. Therefore:

sot ¼

N 2r

ð6Þ

Eq. (6) shows that the theoretical optimal slip only depends on two parameters: the sampling rate N used for recording the signal and the order r of the estimated harmonic component. Therefore, due to the fact that digital relays operate at a constant N, there is a specific theoretical optimal sample slip for each harmonic component. For this value of the slip, Eqs. (7) and (8) provide the constants Ksr and hsr (defined in Section 2) that biunivocally relate the harmonic spectra of both the original signal and the s-Charm wave.

Ksr ¼

Ar 1 1 ¼ ¼ X r 2 sinð180 sot r=NÞ 2

  N  2 sot r ¼0 hsr ¼ bsr  ar ¼ 90 N

ð7Þ

ð8Þ

Therefore, for the slip given by Eq. (6), the r-th harmonic of the s-Charm wave will have double the amplitude of the r-th harmonic of the original signal and will be in phase with it. For example, for r = 1, the Eq. (6) establishes that a separation of a half cycle between samples will generate the greatest amplitude possible for the fundamental component of the s-Charm wave. Then, the range of possible slip values is:

16s6

N 2r

ð9Þ

The use of a slip value outside the range indicated in (9) has no practical sense, as it would imply a greater response time and would not improve the precision of the results. Table 1 shows the possible slips for different values of sampling rate and harmonic components to be estimated (r = 1, . . . and 8). Fig. 2 shows the flowchart for the proposed filtering technique.

Fig. 2. Flowchart for estimating the amplitude and phase for the selected r-th harmonic component at the time corresponding to sample ‘j’.

4. Validation of the optimum sample slip (sot) The parameter sot obtained in the previous section is validated below. For this purpose, a process has been designed in which a series of numerical tests are carried out with Mathematica software using a theoretical signal h(t) and a group of noisy signals yrj ðtÞ.

hðtÞ ¼ C 0 þ Cet=s þ

n X r¼1

Ar cosðxr t þ ar Þ

ð10Þ

143

J. Lázaro et al. / Electrical Power and Energy Systems 56 (2014) 140–146 Table 2 Characteristics of h(t) and yrj ðtÞ signals. C0

C

s, ms

n

Ar

xr, rad/s

ar, rad

100

100

100

8

100/r

2p50r

Rand[p, p]

yrj ðtÞ ¼ hðtÞ þ zj ðtÞ;

j ¼ 1; 2; . . . ; nr

ð11Þ

The values of parameters in (10) and (11) are shown in Table 2. These values have been chosen with sufficient generalisation as to represent a typical fault signal with an important exponential component (C = A1). The time constant of the exponential component is in the middle of the usual range of values for electrical signals. The group of noisy signals (11) is formed by nr different signals, each of which has been obtained superimposing a different noise zj(t) on the theoretical signal (10). These random noises (zj) are within the range of ±1% of the amplitude of the fundamental component of the analysed signal. Each of the nr noisy signals (11) has been filtered using the range of slips between s = 1 and s = 2sot. In each case, the maximum output error (Rjk) is calculated. Two parameters are computed with these maximum errors for each of the sample slips (sk) used in the filtering.

Rmax ¼ MaxfR1k ; R2k ; . . . ; Rnr k g k nr X R þR þþRnr k ¼ n1r Rjk Rakv ¼ 1k 2knr

ð12Þ

j¼1

Figs. 3 and 4 display the evolution of both parameters (Rmax and k according to the sample slip used in the definition of the s-Charm wave for nr = 1000. Both figures confirm the value of sot Rakv )

defined in (6). Similar results have been obtained for the rest of harmonic components and also for other values of the parameters shown in Table 2 (including a time constant range from 20 to 200 ms). It is interesting to point out that Figs. 3 and 4 show the great stability of the optimisation process presented in this section. Slips close to sot provide a reduction of the error very close to the maximum reduction possible. This is due to the shape that the graphs present in their central area. This enables an important reduction of the output error with slip values much lower than sot. Same conclusions are obtained using other ranges of noise. Therefore, based on the results of the validation process carried out, a new practical optimum slip is defined as follows:

sop ¼

sot N ¼ 4r 2

ð13Þ

For this practical optimum slip, the error reduction is very close to its maximum, and in addition, the filter response is significatively faster. Nevertheless, the most adequate choice of slip in each case will always be a compromise between the speed and accuracy required for the protection relay. Figs. 3 and 4 show that this choice is not critical for the correct performance of the method. The sample slip is an additional feature that gives flexibility to the s-CharmDF and allows the method to suit the specific requirements of each relay. To prove the validity of all the above, similar tests have been carried out for the estimation of the rest of harmonic components and for other values of the parameters shown in Table 2. The results were in all cases similar to those shown for the fundamental component in this section. Logically, taking into account Eqs. (6) and (13), the difference between the optimum theoretical slip (sot) and the optimum practical slip (sop) is reduced gradually as the order (r) of the harmonic component increases. 5. Test results To illustrate the behaviour of the s-CharmDF, this section presents the results obtained with different types of electrical signals. Signals of the S1 and S2 types have been defined to test the s-CharmDF response to extreme conditions of the signal to be analysed. Signals of the S3 type have been defined to test its response in the transition from pre-fault to fault conditions. The signals S1 and S2 are defined as follows:

yðtÞ ¼ C 0 þ C 1 et=s1 þ C 2 et=s2 þ Fig. 3. Maximum and average errors (°) in the estimation of the angle of the fundamental component for each slip s.

n X Ar cosðxr t þ ar Þ þ zðtÞ

ð14Þ

r¼1

The values considered for the parameters of these signals are those shown in Table 3. The values of the parameters and the inclusion of two decaying dc-offsets aim to test the s-CharmDF under extreme conditions. The corresponding noisy signals have been obtained taking into consideration the maximum random noise levels for z(t) shown in Table 4. These values have been defined taking into account practical aspects. For example, the 1% value corresponds to the possible maximum internal noise of a relay and the 5% value corresponds to the maximum error in a 5P20 current transformer. The signals corresponding to the cases S1 and S2 are shown in Figs. 5 and 6 respectively. Each figure shows the noisy signal and

Table 3 Characteristics of S1 and S2 signals.

Fig. 4. Maximum and average errors (%) in the estimation of the module of the fundamental component for each slip s.

S1 S2

C0

C1

C2

s1 ms

s2 ms

n

Ar

xr, rad/s

ar, rad

100 200

50 100

50 300

100 200

20 100

8

100/r 100

2p50r

pr/18 Rand[p, p]

144

J. Lázaro et al. / Electrical Power and Energy Systems 56 (2014) 140–146

Table 4 Noise considered in the tests of the s-CharmDF filter. TEST

Maximum random noise level (% A1) for z(t)

S1

S2

1

5

Fig. 8. Module (DFT (dashed line) and s-CharmDF (solid line)) of the fundamental component (S2 signal).

Fig. 5. Noisy signal (dashed line) and s-Charm wave (solid line) for the S1 case.

Fig. 9. Angle (DFT (dashed line) and s-CharmDF (solid line)) of the fundamental component (S1 signal).

Fig. 6. Noisy signal (dashed line) and s-Charm wave (solid line) for the S2 case.

Fig. 10. Angle (DFT (dashed line) and s-CharmDF (solid line)) of the fundamental component (S2 signal).

Fig. 7. Module (DFT (dashed line) and s-CharmDF (solid line)) of the fundamental component (S1 signal).

the s-Charm wave corresponding to the optimum practical slip used for estimating the fundamental component. Figs. 7–10 show the outputs provided by the s-CharmDF method for the module and phase of the fundamental component.

Tables 5 and 6 with numerical indices are included for both Discrete Fourier Transform (DFT) and the s-CharmDF showing the improvement in the convergence process. In the case of modules, the Maximum Percentage Deviation (MPD) [4] values corresponding to the DFT and the s-CharmDF are included. This parameter measures the maximum deviation of the filter output with regard to the reference value to which it converges:

145

J. Lázaro et al. / Electrical Power and Energy Systems 56 (2014) 140–146 Table 5 Comparison (DFT and s-CharmDF) for the module and the angle of the fundamental component corresponding to a S1 signal. Module S1 (r = 1) MPDfou MPDsch

Angle, ° 30.46 4.17

DANGfou DANGsch

13.95 2.81

Table 6 Comparison (DFT and s-CharmDF) for the module and the angle of the fundamental component corresponding to a S2 signal. Module S2 (r = 1) MPDfou MPDsch

Angle, ° 24.30 1.18

DANGfou DANGsch

13.88 0.93

Fig. 13. Angle in degrees (DFT (dashed line) and s-CharmDF (solid line)) of the fundamental component (S3 signal).

Table 7 Comparison (DFT and s-CharmDF) for the module and the angle of the fundamental component corresponding to a S3 signal. Module S3 (r = 1) MPDfou MPDsch

Angle, ° 3.25 0.07

DANGfou DANGsch

1.78 0.04

On the other hand, the MPD index is not useful for angular magnitudes and therefore it is replaced by a more adequate parameter defined as follows:

DANG ¼ MaxjaðtÞ  aref j Fig. 11. Noisy signal in amperes (dashed line) and s-Charm wave in amperes (solid line) corresponding to the S3 case.

Fig. 12. Module in amperes (DFT (dashed line) and s-CharmDF (solid line)) of the fundamental component (S3 signal).

  maxðXðtÞ  X ref Þ   100 MPDð%Þ ¼   X ref

ð15Þ

where X(t) is the instant output of the digital filter and Xref is the actual amplitude of the harmonic component which is being estimated.

ð16Þ

where a(t) is the instant output of the digital filter and aref is the actual phase of the harmonic component which is being estimated. Subscripts fou and sch in Tables 5 and 6 are used to indicate the values corresponding to DFT and s-CharmDF respectively. Signals of the S3 type were generated with the Power System Blockset of MATLAB. As an example of the cases carried out, the results of a S3 signal that includes the steady state (pre-fault) and the transient period (post-fault) are shown in Figs. 11–13. The signal level corresponds to the secondary of a current transformer. The noise level of this signal has a maximum amplitude of 1%. In pre-fault conditions, this signal is composed of a fundamental component (50 Hz). In post-fault conditions, the signal contains a decaying dc-offset, fundamental component and several harmonics. Fig. 11 shows the noisy signal and the s-Charm wave corresponding to the optimum practical slip for estimating the fundamental component. The results obtained for the fundamental component are shown in Figs. 12 and 13. Table 7 shows numerically the great improvement in results provided by the s-CharmDF in the estimation of the fundamental component. Similar results have been obtained in the estimation of other harmonic components, as well as in the application of s-CharmDF to other signals. 6. Conclusions This paper presents the s-CharmDF method, for selective estimation of harmonics in noisy fault signals. A detailed study has been carried out on the behaviour of the new algorithm to guarantee its stability and accuracy in the case of noisy signals. In the different case studies it has been proven that:

146

J. Lázaro et al. / Electrical Power and Energy Systems 56 (2014) 140–146

 Sample slip (s) is a parameter that gives the s-CharmDF filter flexibility to adapt to the different requirements (accuracy/ speed) of the relay without endangering the necessary stability and accuracy of the method.  Optimum theoretical slip (sot), defined in (6), minimises the impact of errors contained in the analysed signal in the output of the s-CharmDF method.  Optimum practical theoretical slip (sop), defined in (13), gives a good compromise between response speed and both stability and accuracy of the s-CharmDF method in case of noisy signals. All the aforementioned features make s-CharmDF suitable to be applied in real time digital relaying. References [1] Xia J, Jiale S, Song G, Wang L, He S, Liu K. Transmission line individual phase impedance and related pilot protection. Int J Electric Power Energy Syst 2011;33:1563–71. [2] Ghorbani A, Mozafari B, Ranjbar AM. Digital distance protection of transmission lines in the presence of SSSC. Int J Electric Power Energy Syst 2012;43:712–9. [3] Gu JC, Yu SL. Removal of DC offset in current and voltage signals using a novel Fourier filter algorithm. IEEE Trans Power Deliv 2000;15(1):73–9. [4] Sidhu TS, Zhang X, Albasri F, Sachdev MS. Discrete-Fourier-transform-based technique for removal of decaying DC offset from phasor estimates. IEE Proc Generat Trans Distrib 2003;150(6):745–52. [5] Guo Y, Kezunovic M, Chen D. Simplified algorithms for removal of the effect of exponentially decaying DC-offset on the Fourier algorithm. IEEE Trans Power Deliv 2003;18(3):711–7. [6] Radojevic ZM. A new spectral domain approach to the distance protection, fault location and arcing faults recognition on transmission lines. Int J Electric Power Energy Syst 2007;29:183–90. [7] Venkatesh K, Swarup KS. Estimation and elimination of DC component in digital relaying. In: International conference on Power, Signals, Controls and Computation (EPSCICON). Thrissur, Kerala, January 2012.

[8] Darwish HT, Fikri M. Practical considerations for recursive DFT implementation in numerical relays. IEEE Trans Power Deliv 2007;22(1):42–9. [9] Dash PK, Krishnanand KR, Patnaik RK. Dynamic phasor and frequency estimation of time-varing power system signals. Int J Electric Power Energy Syst 2013;44:971–80. [10] Sidhu TS, Ghotra DS, Sachdev MS. An adaptive distance relay and its performance comparison with a fixed data window distance relay. IEEE Trans Power Deliv 2002;17(3):691–7. [11] Pan J, Vu K, Hu Y. An efficient compensation algorithm for current transformer saturation effects. IEEE Trans Power Deliv 2004;19(4):1623–8. [12] Soliman SA, Alammari RA, El-Hawary ME. A new digital transformation for harmonics and DC offset removal for the distance fault locator algorithm. Int J Electric Power Energy Syst 2004;26:398–4055. [13] Barbosa D, Monaro RM, Coury DV, Oleskovicz M. Digital frequency relaying based on the modified least mean square method. Int J Electric Power Energy Syst 2010;32:236–42. [14] Muginov GA, Venetsanopoulos AV. A new approach for estimating high-speed analog-digital converter error. IEEE Trans Instrument Measure 1997;46(4):980–5. [15] Kuhlmann V, Sinton A, Dewe M, Arnold C. Effects of sampling rate and ADC width on the accuracy of amplitude and phase measurements in powerquality monitoring. IEEE Trans Power Deliv 2007;22(2):758–64. [16] Lázaro J. Nuevo filtro digital multipropósito para la eliminación de la componente exponencial de señales eléctricas y su aplicación a la estimación de la componente fundamental y armónicos. PhD thesis, Basque Country University, 2008. [17] Lázaro J, Miñambres JF, Zorrozua MA, Larrea B, Sánchez M, Antiza I. New Quick-converge invariant digital filter for phasor estimation. Elect Power Syst Res 2009;79(5):705–13. [18] El-Amin IM. Saturation of current transformers and its impact on digital overcurrent relays. In: IEEE PES transmission and distribution conference and exposition. Caracas, Venezuela, August 2006. p. 1–6. [19] Angell D, Hou D. Input source error concerns for protective relays’, 60th Annual Conference for Protective Relay Engineers. Texas, March 2007. p. 63– 70. [20] He B, Yiguan L, Bo ZQ. Adaptive distance relay based on transient error estimation of CVT. IEEE Trans Power Deliv 2006;21(4):1856–61. [21] Dash PK, Krishnanand KR, Padhee M. Fast recursive Gauss-Newton adaptive filter for the estimation of power system frequency and harmonics in a noisy environment. IET Generate Trans Distrib 2011;5(12):1277–89.