Electrical Power and Energy Systems 83 (2016) 1–6
Contents lists available at ScienceDirect
Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes
Wide area backup protection algorithm for transmission lines based on fault component complex power Seyed-Sattar Mirhosseini a,⇑, Mahdi Akhbari b a b
Department of Electrical Engineering, Shahed University, Tehran, Iran Faculty of Electrical Engineering, Shahed University, Khalij Fars, P.O. Box 33191-18651, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 26 September 2015 Received in revised form 26 March 2016 Accepted 29 March 2016
Keywords: Complex power Fault component Phasor measurement unit (PMU) Wide area backup protection (WABP)
a b s t r a c t In this paper a new wide area backup protection algorithm based on the fault component complex power is proposed to overcome the problem of maloperation of conventional backup protection in highly stressed conditions of power system operation. Firstly, suspected faulty lines are detected using measured values of fault component voltages known as faulted area identification (FAI) criteria. Then, the fault component voltages and currents provided by phasor measurement units (PMUs) are applied to calculate the injected complex power to both terminals of the suspected faulty lines. The ratio between sum and difference of injected complex power to the both terminals of the suspected faulty lines is used as faulted line identification (FLI) criteria. The simulation studies performed on the IEEE 10 generator 39-bus system verify effectiveness of the proposed algorithm under various conditions and fault types. Ó 2016 Elsevier Ltd. All rights reserved.
Introduction Power system and its protection system are planned in such a way that the system can operate through a sequence of credible contingencies without contributing to widespread outages. However, it has been observed that maloperations of conventional protection system in abnormal condition of power system operation have led to large blackouts [1]. In order to comprehensive understanding of issue it is necessary to review conventional protection of transmission line. Recommended practice in the National Grid Company is that transmission line is protected using two main protections including phase and earth fault differential protection and phase and earth fault distance protection with time delayed backup [2,3]. Backup protection includes overreaching zone 2 and zone 3 of distance protection and earth fault overcurrent protection. If fault had not been cleared by main busbar protection, zone 2 elements are expected to detect phase and earth faults on the busbar at theremote end of the line with a time delay of the order of 15–30 cycles. Zone 3 elements are required to detect phase and earth faults on any transmission line connected to the remote end of the main protected line with a typically time delay about 90 cycles. Earth fault overcurrent elements are required to detect resistive earth faults in operating time normally about 1–3 s. ⇑ Corresponding author. Tel.: +98 21 51212066; fax: +98 21 51213564. E-mail address:
[email protected] (S.-S. Mirhosseini). http://dx.doi.org/10.1016/j.ijepes.2016.03.056 0142-0615/Ó 2016 Elsevier Ltd. All rights reserved.
Recent blackouts have attracted attention to the role of zone 3 of distance protection schemes. It has been reported that unwanted trips of zone 3 caused by unexpected loading conditions such as overload and flow transfer have often contributed to the cascading outages eventually result in large blackouts [1]. Moreover, setting mistakes due to complex principle of conventional backup protection setting lead to hidden failures [4,5]. Different solutions have been proposed to deal with this issue ranged from suggestion to completely eliminating the zone3 from the protection system to using blocking type distance protection schemes. With recent developments in computer networks, communication, information technology and advent of synchronised phasor measurement unit (PMU) and wide area measurement system (WAMS), nowadays it is possible and necessary to have a wider view to the power system in order to eliminate many shortcomings of conventional protection system [6,7]. One of the most effective solutions under investigation, which is proposed to eliminate conventional backup protection, is using WAMS to provide all backup protection as wide area backup protection (WABP) [8]. Recently, some algorithms have been proposed for WABP. A protection scheme based on comparing magnitude of positive sequence voltages and difference between positive sequence current phase angles at both terminals of lines is proposed in [9]. This algorithm is unable to reflect high resistance faults. An algorithm based on fault component voltage distribution is proposed in [10]. In this algorithm the measured values of fault component voltage and current at one terminal of the line are used to estimate
2
S.-S. Mirhosseini, M. Akhbari / Electrical Power and Energy Systems 83 (2016) 1–6
Nomenclature FAI PMU FLI WABP WAMS PCR FCF U1 U2 U0 DU Duf
faulted area identification phasor measurement unit faulted line identification wide area backup protection wide area measurement system protection correlation region fault correlation factor positive sequence bus voltage negative sequence bus voltage zero sequence bus voltage fault component bus voltage fault component voltage source
the fault component voltage at the other terminal. Then the faulted element is identified by using the ratio between the measured and estimated voltages. Moreover, the faulted element identification algorithm is accelerated by applying a faulted area detection scheme. This algorithm is able to identify faulted line in complex conditions such as high resistance fault and flow transfer. It needs to line impedance to estimate voltage at the other terminal of line. An algorithm based on fault steady state component is proposed in [11] to identify faulted branch. In this algorithm in normal condition of power system operation, on the basis of the network topology and PMUs placement, buses are classified into subsets named protection correlation regions (PCRs). When a fault occurs, the fault correlation region is determined by analysing the fault steady state component of differential current injected into each PCR. Then a fault correlation factor (FCF), using prefault and fault component currents and voltages, is calculated to identify the faulted branch. Unlike the two above-mentioned algorithms, this algorithm does not require to install PMU at all buses. The algorithm identifies high resistance fault correctly but it uses line impedance and bus impedance matrix of each PCR in its calculation. An analytical approach using dispersed PMUs and bus impedance matrix is proposed in [12]. In this approach, fault zone is first detected by local PMUs, then the suspected faulty lines are diagnosed and finally the fault line is identified and the fault point is located. This approach is successful even in case of the fault line, which is not equipped with PMU on either side. This paper proposes a new wide area backup protection algorithm based on fault component sequence voltages distribution and fault component complex power. The algorithm comprises two parts including faulted area identification (FAI) and faulted line identification (FLI). Phasor of sequence currents and voltages of all buses of power system, provided by PMUs, are monitored in the centre of WABP. In order to reduce the calculations some lines are selected as suspected faulty lines using FAI criteria, which are based on fault component sequence voltages distribution in the power system. Then, FLI criteria, which are defined based on injected complex power to both terminals of line, are deployed to identify faulted line. The algorithm uses only measured voltages and currents to calculate the criteria and it does not need to system parameters such as line impedance and bus impedance matrix. It also identifies faulted line in different conditions such as high resistance fault and flow transfer.
z zf zi zj a Si Sj k k0 k1 k2
line impedance ground fault impedance network equivalent impedance at bus i network equivalent impedance at bus j percentage of line length from bus i at which fault has occurred injected complex power to bus i injected complex power to bus j faulted line identification criteria zero sequence faulted line identification criteria positive sequence faulted line identification criteria negative sequence faulted line identification criteria
the sequence voltages distribution along the faulted line and equivalent sources of other elements of the power system is shown in Fig. 1. The fault is a resistive fault. According to Fig. 1, it can be concluded that: (1) The positive sequence voltage amplitude is minimum at the fault point and it is maximum at the source points. (2) The zero and negative sequence voltages amplitudes are maximum at the fault point and it is minimum at the equivalent source points and (3) among the voltages of all buses of a power system, the positive sequence voltage amplitudes of the two terminals of the faulted line are minimum, and the zero and negative sequence voltage amplitudes of the two terminals are maximum. Considering above points, the positive sequence voltage amplitudes of two buses at both terminals of the faulted line are minimum and the zero and negative sequence voltages amplitudes of these two buses are maximum among the voltages of the other buses of power system. On the basis of aforementioned explanation, the bus with minimum positive sequence voltage or maximum zero and negative
Faulted area identification When a fault occurs on a transmission line, the sequence currents in power system are calculated using symmetrical components method. Considering the sequence currents and impedances,
Fig. 1. Distribution of sequence voltages in a faulted power system.
3
S.-S. Mirhosseini, M. Akhbari / Electrical Power and Energy Systems 83 (2016) 1–6
sequence voltages can be considered as the nearest bus to the faulted line. Using only the positive sequence voltage it is possible to detect the occurrence of all types of fault in power system. However, some factors like short-term existence of the positive sequence fault component voltage in comparison with the zero and negative sequence voltages [13] and effect of measurement errors, require that positive, negative and zero sequence voltages to be simultaneously deployed for determination of suspected faulty lines. The synchronised phasors U 1 , U 2 and U 0 are provided by PMUs installed at each bus. In order to identify faulted area U 1 is ranked in increasing order and U 2 and U 0 are ranked in decreasing order. To overcome measurement errors and increasing reliability of algorithm, the three buses corresponding to top three values of U 1 , U 2 and U 0 are selected as the nearest buses to the faulted line. After selection of the nearest buses to the faulted line, the lines connected to the selected buses constitute faulted area. Then, the faulted line identification criteria are calculated for each line in the faulted area and the faulted line is identified. Faulted line identification In the presence of a fault in the power system, superposition theorem [14] allows that the voltages and currents can be considered as prefault and fault components. The fault equivalent network of a two terminal power system is shown in Fig. 2. In this figure DUs, Duf show the fault component bus voltages and the fault component voltage source, respectively. zi and zj are the equivalent impedances at buses i and j, z is the line impedance, zf is the ground fault impedance and a 2 ½0; 1 shows the percentage of the line length from bus i at which the fault has occurred. In case of an external fault, considering the equivalent network, voltages and currents at two terminals of the unfaulted line are calculated as follows:
DIi ¼ DIj ¼
zi Duf zf ðzi þ zj þ zÞ þ zi ðz þ zj Þ
DU i ¼ ðz þ zj ÞDIi DU j ¼ zj DIj ¼ zj DIi
Si ¼
zjzi j2 jDuf j2 zf ðzi þ zj þ zÞ þ zi ðz þ zj Þ
ð6Þ
Si Sj ¼
ðz þ 2zj Þjzi j2 jDuf j2 zf ðzi þ zj þ zÞ þ zi ðz þ zj Þ
ð7Þ
Considering (6) and (7) it can be concluded that in an unfaulted line, absolute sum of injected complex power to terminals of the line is smaller than absolute difference of injected complex powers.
2 jSi þ Sj j z <1 ¼ jSi Sj j z þ 2zj
ð8Þ
In case of an internal fault, considering the equivalent network, voltages and currents at two terminals of the faulted line are calculated as follows:
DI i ¼
zj þ ð1 aÞz Duf zf ðzi þ zj þ zÞ þ ðzi þ azÞðzj þ ð1 aÞzÞ
ð9Þ
DI j ¼
ðzi þ azÞ Duf zf ðzi þ zj þ zÞ þ ðzi þ azÞðzj þ ð1 aÞzÞ
ð10Þ ð11Þ
DU j ¼ zj DIj
ð12Þ
zj þ ð1 aÞz 2 jDuf j2 Si ¼ zi M
ð13Þ
z þ az 2 i 2 Sj ¼ zj jDuf j M
ð14Þ
ð1Þ
M ¼ zf ðzi þ zj þ zÞ þ ðzi þ azÞðzj þ ð1 aÞzÞ
ð15Þ
ð2Þ
Sum and difference of injected complex power to terminals of the faulted line are calculated as:
ð3Þ
ðz þ zj Þjzi j2 ¼ jDuf j2 zf ðzi þ zj þ zÞ þ zi ðz þ zj Þ
ð4Þ
zj jzi j2 jDuf j2 zf ðzi þ zj þ zÞ þ zi ðz þ zj Þ
ð5Þ
Sj ¼ DU j DIj ¼
Si þ Sj ¼
DU i ¼ zi DIi
Using calculated voltages and currents in (1)–(3), injected complex powers to terminals of the line are calculated as follows:
DU i DIi
Sum and difference of injected complex power to terminals of the line are obtained as:
Si þ Sj ¼ Si Sj ¼
zi jzj þ ð1 aÞzj2 þ zj jzi þ azj2 jMj2 zj jzi þ azj2 jMj2
jDuf j2
jDuf j2
ð16Þ ð17Þ
It can be also seen from (16) and (17) that in a faulted line, absolute sum of injected complex power to terminals of the line is greater than absolute difference of injected complex powers.
jSi þ Sj j jzi jzj þ ð1 aÞzj2 þ zj jzi þ azj2 j ¼ >1 jSi Sj j jzj jzi þ azj2 j
ð18Þ
Similar equations will be obtained for the positive and zero sequence networks. According to aforementioned explanations it can be concluded that when a fault occurs on a line in power system: (1) In case of an unfaulted line:
k¼
jSi þ Sj j <1 jSi Sj j
ð19Þ
(2) In case of a faulted line:
k¼
Fig. 2. Fault component equivalent system.
jSi þ Sj j >1 jSi Sj j
ð20Þ
Eqs. (19) and (20) are deployed to define the faulted line identification criteria using the positive, negative and zero sequence fault components as follows:
4
S.-S. Mirhosseini, M. Akhbari / Electrical Power and Energy Systems 83 (2016) 1–6
jSi1 þ Sj1 j jSi1 Sj1 j jSi2 þ Sj2 j k2 ¼ jSi2 Sj2 j jSi0 þ Sj0 j k0 ¼ jSi0 Sj0 j
k1 ¼
ð21Þ ð22Þ ð23Þ
When a fault occurs depending on the fault type some or all of k1; k2; k0 take a value greater than 1 for the faulted line and a value less than 1 for all the other lines. These factors are introduced as faulted line identification criteria. It should be noted that there is no zero sequence fault component during phase to phase fault. Therefore, k0 is not defined for phase to phase fault. Because of similar reason k0 and k2 are not defined for three phase fault. There will be no concern about this issue because it is possible to detect the faulty line using only k1 as it is shown in simulation study. However, simultaneous using of k1; k2; k0 enhance the reliability of the algorithm for asymmetrical faults. FLI criteria analysis For analysing characteristics of the FLI criteria, factors which may affect the criteria should be considered. These factors are ground fault impedance, fault position, faulted line impedance and equivalent network impedances of both terminals of the faulted line.
Fig. 3. Effect of line and network impedances on k.
Effect of ground fault impedance As it is obvious from (18); the ground fault impedance zf has no effect on the criteria. The simulation results shown in Table 4 in part 4 approve this subject. Effect of network impedances and fault position In order to analyse the effect of fault point position, the line impedance and the equivalent network impedances, suppose that ¼ zz , y ¼ zz . Substituting x, y in (18), the FLI criteria can be expressed i
j
as:
k¼
yj1 þ ð1 aÞyj2 þ yj yx þ ayj2 xj yx þ ayj2
ð24Þ
Neglecting the system resistances and attributing typical real values to x and y, the distribution of k along the faulted line is shown in Fig. 3. It can be seen that, if the ratio between the line impedance and the equivalent network impedances at both terminals of the faulted line are equal, k is independent of the parameter ‘‘a” and has a constant value, i.e. 2. In case of x > 1 and for all y, k decreases by increasing in parameter ‘‘a”. In case of x < 1 and for all y, k increases by increasing in parameter ‘‘a”. However, in all cases k factor is greater than 1. To evaluate the effect of impedances, the parameter ‘‘a” is set as 0.5 and k is plotted versus x, y in Fig. 4. As it is obvious, k increases by increasing in x and decreases by increasing in y. But, its value is greater than 1 for all values of x, y. In actual high voltage and extra high voltage systems z is generally greater than zi and zj . Therefore, in practice there will be no condition in which k is close to 1. Simulation studies The IEEE 10 generator 39 bus test system [15], as shown in Fig. 5, is used to evaluate the performance of the proposed algorithm using DigSilent Power Factory and Matlab softwares. Since the PMU data communication delay is much less than time delay
Fig. 4. Distribution of k along the faulted line.
in backup protection (about 150 ms in comparison with 90 cycle in case of using fibre optic link) [16,17], communication time delay is ignored in simulations. FAI and FLI test In order to test the FAI and FLI parts of the proposed WABP algorithm, different types of fault at different locations are incepted on lines L(4-14), L(17-18), L(8-9). As mentioned before there is no zero sequence fault component during phase to phase fault therefore, k0 is not defined for this type of fault. Similarly k0 and k2 are not defined for three phase fault. The values of k1, k2 and k0 of suspected faulty lines identified by FAI in case of a single phase to ground fault on L(4-14) at a = 0.5 is shown in Fig. 6. It is obvious that FLI criteria are greater than 1 for faulted line and are much smaller than 1 for all the other lines. The results are shown in Tables 1–3. The results show that the FAI properly detects the nearest buses to the fault and suspected faulty lines. The positive, negative and zero sequence FLI criteria are also greater than 1 for all types of fault at different locations.
5
S.-S. Mirhosseini, M. Akhbari / Electrical Power and Energy Systems 83 (2016) 1–6 Table 2 Faults on L(17-18). a (%)
Fault type
FAI
20
AG BC BCG ABC
B3, B3, B3, B3,
50
AG BC BCG ABC
B17, B18, B27 B17, B18, B27 B3, B17, B18 B3, B17, B18
80
AG BC BCG ABC
B17, B17, B17, B17,
B17, B17, B17, B17,
B18 B18 B18 B18
B18, B18, B18, B18,
B27 B27 B27 B27
K1
K2
K0
23.070 23.071 23.071 22.155
13.543 13.543 13.543 –
19.181 – 19.181 –
4.2758 4.2759 4.2759 4.2769
4.0695 4.0695 4.0695 –
8.457 – 8.457 –
2.3923 2.3924 2.3924 2.3973
2.4286 2.4286 2.4286 –
3.8119 – 3.8119 –
Table 3 Faults on L(8-9). a (%)
Fault type
FAI
K1
K2
K0
20
AG BC BCG ABC
B7, B7, B7, B7,
B8, B8, B8, B8,
B9 B9 B9 B9
1.6471 1.6475 1.6475 1.6751
1.8578 1.8578 1.8578 –
1.6806 – 1.6806 –
50
AG BC BCG ABC
B5, B5, B5, B5,
B7, B7, B7, B7,
B8 B8 B8 B8
4.4207 4.4155 4.4152 4.3923
17.893 17.893 17.893 –
9.0526 – 9.0526 –
80
AG BC BCG ABC
B5, B5, B5, B5,
B7, B7, B7, B7,
B8 B8 B8 B8
1.3207 1.3207 1.3207 1.3159
1.889 1.889 1.889 –
1.6071 – 1.6071 –
Fig. 5. The IEEE 10 generator 39 bus test system.
10
k1
5
k2 k0
0
Table 4 Effect of ground fault impedance. Fig. 6. FLI criteria of suspected faulty lines in single phase to ground fault.
Table 1 Faults on L(4-14). a (%)
Fault type
FAI
20
AG BC BCG ABC
B4, B4, B4, B4,
B13, B13, B13, B13,
B14 B14 B14 B14
50
AG BC BCG ABC
B4, B4, B4, B4,
B12, B12, B12, B12,
B14 B14 B14 B14
AG BC BCG ABC
B4, B4, B4, B4,
B13, B13, B13, B13,
B14 B14 B14 B14
80
K1 3.3924 3.3923 3.3923 3.4242 19.635 19.635 19.635 16.425 2.7579 2.758 2.758 2.7435
K2
K0
3.816 3.816 3.816 –
8.8647 – 8.8647 –
15.589 15.589 15.589 –
9.3942 – 9.3942 –
2.8933 2.8933 2.8933 –
3.3068 – 3.3068 –
Effect of ground fault impedance For analysing the effect of ground fault impedance on the criteria, different single phase to ground faults have been incepted on line L(4-14) and the criteria are calculated. The simulation results are shown in Table 4. As results show, the ground fault impedance has no effect on the FLI criteria. Effect of flow transfer In order to prevent maloperation of the protection algorithm it is essential that the criteria distinguish flow transfer condition
a (%)
Fault resistance ðXÞ
20
0 50 300
50
0 50 300
80
0 50 300
K1 3.3924 3.3924 3.3924 19.635 19.635 19.635 2.7579 2.7579 2.7579
K2
K0
3.816 3.816 3.815
8.8647 8.8647 8.8646
15.589 15.589 15.589
9.3942 9.3942 9.3941
2.8933 2.8933 2.8933
3.3068 3.3068 3.3066
from internal fault. For analysing the effectiveness of the FLI criteria in flow transfer condition, two heavily loaded lines are selected in such a way that fault occurrence on one line results in a large flow transfer to the other. Implementing load flow calculations, L (23-24) and L(21-22) are selected as two heavily loaded lines. A three phase fault is occurred on L(21-22) at t = 0.5 s and is cleared at t = 0.6 s. The positive sequence currents carried by the two lines are shown in Fig. 7. As it is shown, a large part of carried power by L (21-22) is transferred on L(23-24). The criteria values of the two lines during the flow transfer condition are shown in Fig. 8. The criteria of L(23-24) is less than 1 before fault occurrence on the other line, at the fault clearing time k1 violates than 1 for a time about three cycles then decreases rapidly to a very small value. This is due to an error in phasor computation created by switching transient occurred by opening the breakers. Since the time delay of backup protection is much greater than three cycles, this subject has no effects on performance of the algorithm. This subject can also be solved by setting a suitable threshold. Therefore, the algorithm can be blocked successfully in flow transfer condition.
6
S.-S. Mirhosseini, M. Akhbari / Electrical Power and Energy Systems 83 (2016) 1–6
and voltages of all buses of the power system, provided by phasor measurement units, are monitored in the centre of WABP. In order to reduce the calculations some lines are selected as suspected faulty lines using FAI criteria, which are based on the fault component sequence voltages distribution in the power system. Then, FLI criteria, which are calculated using sum and difference of injected complex power to both terminals of the suspected faulty lines, are defined to identify the faulted line. The algorithm uses only voltages and currents phasor and it is independent from system parameters such as line impedance and bus impedance matrix. The simulation results on IEEE 39 bus test power system validate the proposed algorithm and show the effectiveness of the algorithm to identify faulted line in different conditions such as high resistance fault and flow transfer. References
Fig. 7. Positive sequence currents of L(21-22), L(23-24).
Fig. 8. k1 of L(23-24) under flow transfer condition.
Conclusion A new wide area backup protection algorithm based on the fault component complex power is proposed. The algorithm comprises two parts including the faulted area identification (FAI) and the faulted line identification (FLI). Phasor of the sequence currents
[1] Horowitz SH, Phadke AG. Third zone revisited. IEEE Trans Power Deliv 2006;21 (January):23–9. [2] Phadke AG, Thorp JS. Computer relaying for power systems. Taunton (UK): Research Studies Press; 1988. [3] Anderson PM. Power system protection. New York: IEEE Press; 1999. [4] Phadke AG. Synchronized phasor measurements and their applications. New York: Springer; 2008. p. 211. [5] Novosel D, Bartok G, Henneberg G, Mysore P, Tziouvaras D, Ward S. IEEE PSRC report on performance of relaying during wide-area stressed conditions. IEEE Trans Power Deliv Jan. 2010;25(1):3–16. [6] Phadke AG. Synchronized phasor measurements in power systems. IEEE Comput Appl Power 1993;6(April):10–5. [7] Suranyi A, Bertsch J, Reinhardt P. Use of wide area monitoring, protection and control systems to supervise and maintain power system stability. In: Proc 8th IEE int conf AC DC power transm; 2006. p. 200–3. [8] Phadke AG, Thorp JS. History and applications of phasor measurements. In: Proc 2006 IEEE PES power syst conf expo; 2006. p. 331–5. [9] Eissa MM, Masoud ME, Elanwar MMM. A novel back up wide area protection technique for power transmission grids using phasor measurement unit. IEEE Trans Power Deliv 2010;25(1):270–8. [10] He Zhiqin, Zhang Zhe, Chen Wei, Malik Om P, Yin Xianggen. Wide-area backup protection algorithm based on fault component voltage distribution. IEEE Trans Power Deliv 2011;26(4):2752–60. [11] Ma Jing, Li Jinlong, Thorp James S, Arana Andrew J, Yang Qixun, Phadke AG. A fault steady state component-based wide area backup protection algorithm. IEEE Trans Smart Grid 2011;2(3):468–75. [12] Salehi-Dobakhshari A, Ranjbar AM. Application of synchronised phasor measurements to wide-area fault diagnosis and location. Gener Trans Distribut IET 2014;8(4):716–29. [13] Gao H, Crosslay PA. Design and evaluation of a directional algorithm for transmission-line protection based on positive sequence fault components. Proc Inst Elect Eng Gen Transm Distrib 2006;153(6):711–8. [14] Gao H, Crossley PA. Directional relay for EHV transmission lines using positive sequence fault components. 2005 IEEE Russia Power Tech 2004:1–5. [15] Pai MA. Energy function analysis for power system stability. Boston (MA): Kluwer; 1989. p. 223–7. [16] Chenine M, Zun K, Nordstrom L. Survey on priorities and communication requirements for PMU-based applications in the Nordic Region. Power Tech 2009. [17] Kim M, Damborg MJ, Huang J, Venkata SS. Wide-area adaptive protection using distributed control and high-speed communication. In: Proc of 2002 the 14th power system computation conference, Sevilla, Spain; 2002.