Fault section estimation in electric power systems using an optimization immune algorithm

Electric Power Systems Research 80 (2010) 1341–1352

Contents lists available at ScienceDirect

Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr

Fault section estimation in electric power systems using an optimization immune algorithm Fábio Bertequini Leão ∗ , Rodrigo A.F. Pereira, José R.S. Mantovani Research Group in Electric Power System Planning, Department of Electrical Engineering, Engineering College of Ilha Solteira, São Paulo State University, UNESP, Post Office Box 31, 15385-000, Ilha Solteira, SP, Brazil

a r t i c l e

i n f o

Article history: Received 4 February 2010 Received in revised form 24 May 2010 Accepted 27 May 2010 Available online 20 June 2010 Keywords: Fault section estimation Immune algorithm Power system protection Protective relaying

a b s t r a c t This paper proposes a methodology based on the unconstrained binary programming (UBP) model and an optimization immune algorithm (IA) to estimate fault sections in electric power systems. The UBP model is formulated using the parsimonious set covering theory for associating the alarms of the protective relay functions informed by the SCADA (supervisory control and data acquisition) system and the expected states of the protective relay functions. The IA is developed to minimize the UBP model and to estimate the fault sections as quickly and reliably as possible. The proposed methodology is tested using part of the South-Brazilian electric power system. The control parameters of the IA are set to reach the maximum computational efficiency and reduction of the processing time. The results show the methodology’s potential to estimate fault sections in electric power system control centers in real-time. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Modern control centers (CCs) have centralized SCADA systems which perform remote monitoring, control and real-time data acquisition from a number of substations. The operation status of sections, such as busbars, transformers and transmission lines in each substation are supervised by intelligent electronic devices (IEDs), e.g., multifunctional digital relays. The remote management of power system operations is performed by the energy management system (EMS), where the information is concentrated, processed and displayed by a human machine interface (HMI) to be analyzed and interpreted by operators. When a fault occurs in any section of the electric power system, a large flow of information, due to alarms and events produced by IEDs, overwhelms the operator at the CC. Based on the operator’s previous knowledge about: (1) the electric power system topology; (2) the electric power system behavior under various operating conditions; and (3) the corresponding performance philosophy of the protection system, he must quickly estimate the fault section in order to restore the system to normal conditions as soon as possible. Developing efficient methodologies to assist operators to estimate faults quickly and accurately is very important and also challenging in terms of its complexity. These methodologies can be integrated to the EMS as an additional function, in order to minimize equip-

∗ Corresponding author. Tel.: +55 18 3743 1150; fax: +55 18 3743 1163. E-mail addresses: [email protected] (F.B. Leão), [email protected] (R.A.F. Pereira), [email protected] (J.R.S. Mantovani). 0378-7796/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.epsr.2010.05.010

ment damages and also to reduce the electric supply interruption time. Estimating fault sections in electric power systems using protective relay alarms and circuit breaker status is difficult for bulk power systems that keep many protective and control devices installed. Furthermore, there are many decision variables to consider and large dimension solution space due to the combinatorial nature of the problem. Since Wollenberg proposed in 1986 [1] the use of expert systems (ES) to manipulate alarms in CCs, many works have been published in this area to solve the fault section estimation problem [2–4]. ES use a rule-based knowledge system and an inference engine to estimate the sections under fault. Although widely used, a limitation of ES is their generalization inability. The validation and maintenance of large knowledge bases also impose various difficulties. Contrary to the static knowledge base of ES, artificial neural networks (ANNs) are a promising technique to solve the problem. Capable of generalization, ANNs are ideal to deal with partial or corrupted alarms. Most studies in this area use only ANN [5–7], whereas a few works use hybrid methods combining ANN with fuzzy logic (FL) [8]. Despite the advantages of ANN, there are some drawbacks when applied to large electric power systems, such as: slow convergence on the training process, and difficulties in determining the network parameters such as hidden units, quantity of layers and training rate. FL is adequate to model imprecise systems and presents great flexibility. To deal with the inherent uncertainties of the protection systems, most of the developed methods in the literature use FL combined with other techniques [9–12]. On the other hand, FL presents the disadvantage of choosing membership functions by means of historical data, experience or

1342

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

trial-and-error. Genetic algorithm (GA) has proved to be successful for solving the problem [13–16]. GA presents several advantages, such as robustness, easy to code and reduced processing time. However, GA has difficulties in finding multiple solutions or suboptimal solutions and, in some cases, presents premature convergence for local optimum. Although every cited reference has developed special strategies to overcome these difficulties, the applications of GA are limited to test systems. Petri Nets (PN) is an adequate technique to describe and analyze information flow. PN is suitable to modeling concurrent asynchronous systems as protection systems [17]. The main disadvantage is its generalization incapacity and the graphical representation of the protection system is difficult to manage when applied to large electric power systems. Tabu Search (TS) [18], Ant Colony Optimization (ACO) [19] and IA [20] are solution techniques also used for solving fault section estimation problems. However they have not been fully explored and require further investigation in order to apply them to large real-life systems. This work proposes a new methodology to solve the fault section estimation problem. The main objective of the methodology is to overcome all drawbacks of early works and provide assistance to the operator. The methodology is based on the UBP model and an optimization IA. The UBP model, established by using the parsimonious set covering theory [15,18,21] seeking for an association between the alarms of the protective relay functions informed by the SCADA system and the relays expected states if a fault occurs in its protective zone. The optimization IA is inspired on the clonal selection principle and developed to estimate sections under faults faster and more reliably by the solution of combinatorial optimization UBP model. The proposed methodology is tested using part of the South-Brazilian electric power system [6]. Control parameters of the IA are tuned to increase the computational efficiency, assuring robustness, precision, searching capability for suboptimal solutions and reducing the processing time of the algorithm. The performance of the proposed IA is compared to a dedicated genetic algorithm [16]. 2. Mathematical model for the fault section estimation problem This paper proposes a new methodology that combines the mathematical models of protection system performance and of UBP in order to adequately solve the fault section estimation problem. This methodology is described in detail below. 2.1. Mathematical model of protection system operations The mathematical model of protection system operations is a set of equations of the expected state of the protective relay functions. Each equation of this set mathematically models the operation logic of the protective relay functions (differential, distance and overcurrent) located at the power system substations. This equation set is generated from a unique generic equation considering: (1) the electric system topology; (2) the relay alarms and the state of the circuit breakers informed by the SCADA system and (3) the protection philosophy used by experts on the protection specification, selectivity and coordination [22]. 2.1.1. Operation logic and protection philosophy In this work the overcurrent function is employed for backup protection of busbars and transformers since the differential functions are the main protection of these sections. The distance functions are employed for main and backup protection of transmission lines and backup protection for busbars. Two zones are considered for the distance functions: zone 1 or MP (main protection function) that reaches 100% of the line length; and zone 2 or BP (backup protection function) that reaches 150–200% of the line

length where the relay is located. Zone 2 also protects the busbars on the opposite terminal of the relay allocation point. Furthermore, it is considered that the distance functions are directional and do not have reverse zone, in other words, they protect only the section of the circuit as well as the adjacent lines, from the busbar where the relay is located. Fig. 1 illustrates a generic electric system used to establish the protection philosophy which is used to develop the mathematical model. The system is composed of four substations a, b, c, and d. Each substation has a DC source supplying the auxiliary substation service (relays and command circuits of circuit breaker, among others). If a DC source fails, it is expected that any protective function of the relay does not operate although there is a fault. Furthermore, one or more functions of a relay cannot trip due to a failure of the device itself. The part of a circuit between two circuit breakers is considered a section. For example, section T1 is within circuit breakers CB1 A and CB1 B and busbar B within circuit breakers CB1 B and CB2 B . Other sections are analogously defined. For instance, when there is a fault on line L1 of the system in Fig. 1, the main protection must trip, i.e., the functions MP1 B and MP1 C must send a trip signal to CB2 B and CB1 C respectively, provided the corresponding DC sources of substations a and b have not failed. The backup protection of line L1 is composed of functions BP1 B , BP1 C , BP1 D , and BP2 E and these functions must trip the circuit breakers when the main protection fails, provided the corresponding DC sources of substations a, b, c and d have not failed. The functions BP1 B and/or BP1 C (L1 local backup) must trip when the main protections fail and the functions BP1 D and BP2 E (L1 remote backup) must trip only if the protections located on line L1 were not able to isolate the fault. 2.1.2. Expected states of the protection functions From the operation logic of the protective functions and protection philosophy, a generic equation which represents the expected state of each protective relay functions is mathematically written as: fi = MAX j ∈ Ji

⎧ ⎨

⎡

⎛

sj × ⎣1 − ORC ⎝

⎩



⎞⎤

csck ⎠⎦ , sj × (1 − fsct )

k ∈ ˝li,j

× (1 − bh )

t ∈ ˝ti,j

⎫ ⎬ ⎭

(1)

Ji = ˝mi ∪ ˝bi ; h ∈ ˝hi ; ˝ci ; i = 1, . . . , nf ˝dp ; p = 1, . . . , ncb where nf: total number of protective relay functions; ncb: total number of circuit breakers; fi : expected state of the protective relay function i (1: actuated, 0: non-actuated); sj : state of section j (busbars, lines and transformers) (1: fault, 0: normal); csck : state of the circuit breaker k informed by the SCADA system (1: open, 0: closed); fsct : state of the protective relay function t informed by the SCADA system (1: actuated, 0: non-actuated); bh : state of the DC source of substation h (1: failure, 0: normal); ORC(x): results in 1 if x is equal or greater than 1. Otherwise, results in 0; MAX{·}: extracts the maximum value of the considered elements; Ji : index set of protected section(s) by protective function i; ˝mi : index set of protected section(s) by the main protective function i; ˝bi : index set of section(s) protected by the backup protective function i; ˝hi : index set of DC source that supplies the relay with protective function i; ˝ci : index set of circuit breaker(s) tripped by protective function i; ˝dp : Index set of sections connected to circuit breaker p; ˝ti,j : alarm index set of the main protective function that trip the same circuit breaker and protects the same section j and the backup protective function i. This set is empty for main protective functions; ˝li,j : index set of circuit breaker(s) connected to the circuit path between the circuit breaker tripped

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

by the backup protective function i and the protected section j. The index of the circuit breaker tripped by function i must be included on the set when it is connected to section j. This set is empty for main protective functions. Eq. (1) is a generic equation that must be applied to every nf protective function of the electric system. The quantity of terms between brackets depends on the quantity of sections that each function i protects, i.e., the quantity of elements pertaining to each set Ji . Furthermore, the indexes h, k and t of the independent variables are obtained by the sets ˝hi , ˝li,j and ˝ti,j , respectively. The sets Ji , ˝mi , ˝bi , ˝hi , ˝ci and ˝dp are named basic sets because they are formed from known data that can be easily extracted from the electric system. For the system illustrated in Fig. 1, and considering DP1 A = f1 , OR1 A = f2 , BP1 B = f3 and CB1 A = CB1 , CB1 B = CB2 , CB2 B = CB3 , the sets aforementioned are formed as follows: J1 = {1}, J2 = {1,6}, J3 = {3,7,8,10}, ˝m1 = {1}, ˝m2 = ˝m3 = ø, ˝b1 = ø, ˝b2 = {1,6}, ˝b3 = {3,7,8,10}, ˝h1 = ˝h2 = ˝h3 = {1}, ˝c1 = ˝c2 = {1}, ˝c3 = {3}, ˝d1 = {1,6}, ˝d2 = {2,6} and ˝d3 = {2,7}. Forming the sets ˝li,j and ˝ti,j from their definitions is not trivial and is difficult for large electric power systems. A generic method that can be applied to any system to obtain ˝li,j and ˝ti,j consists of transforming the basic sets Ji , ˝mi , ˝bi , ˝ci and ˝dp in binary basic matrices 0–1, and uses the conjunction and disjunction operations of sets in the Boolean operation form AND and OR element-byelement of the matrices. The binary basic matrices present some advantages in relation to the sets since they store the system information in a simple and compacted way. Furthermore, the matrices are scalable and can be easily enlarged when there is an expansion of the system. Each line of the basic matrices represents an index set and each column represents one element that pertains (1) or not (0) to this set. Eq. (2) presents the formation of the binary basic matrix generated by the sets Ji .

(2)

M (nf,ns) = Mm(nf,ns) ⊕ Mb(nf,ns)

(3)

where ns: total number of sections of the electric system; M(nf,ns) : binary matrix generated by sets Ji ; Mm(nf,ns) : binary matrix generated by sets ˝mi ; Mb(nf,ns) : binary matrix generated by sets ˝bi ; ⊕: Boolean disjunction operator – logic OR. Eq. (3) means that matrix M can be separated in two matrices that allow distinguishing which are main functions and which are backup protection. Element (i,j) of the matrices M, Mm and Mb is 1 when the function i protects section j, otherwise it is 0. For example, considering the sets formed for the system in Fig. 1, the three first rows of matrix M are as follows:

The same reasoning is used to generate the binary matrices given by (5)–(6) from the sets ˝ci and ˝dp , respectively.

(5)

(6) where Mc(nf,ncb) : binary matrix generated by sets ˝ci ; Md(ncb,ns) : binary matrix generated by sets ˝dp . Matrix in (5) shows that function 1 trips the circuit breaker 1, function 2 trips the circuit breakers ncb-1 and ncb, etc. The matrix in (6) contains the information of which circuit breaker is connected in such section. For example, from (6), circuit breaker 1 connects sections 1 and 2 and circuit breaker 2, the sections ns-1 and ns, etc. Appendix A shows details about the rules for forming the sets ˝li,j and ˝ti,j from the binary basic matrices defined by (2), (3), (5) and (6). 2.2. Unconstrained binary programming model Solving the fault section estimation problem means to estimate the state of all sections in the system that can explain a set of informed alarms by the SCADA system. The challenges in this area are related to the mathematical model of the protection system and solution technique in such a way that the solution hypotheses can be feasibly generated. The parsimonious set covering theory [21] applied to solve the fault section estimation problem [15,18] states that a plausible hypothesis or explanation for a set of informed alarms is achieved when: (1) The solution (fault sections) is a cover of the informed alarms, i.e., the informed alarms can be accounted for by the fault sections in the solution. (2) The informed alarms are as consistent as possible with the expected states of the protective functions of the relays that were calculated for the solution. (3) The solution with a minimum number of fault sections is preferred to explain the fault occurrence. The three conditions established by the parsimonious set covering theory can be translated mathematically by minimizing the following objective function: Min E = k1 ×

 nf

 Coi

+ k2 ×

i=1

×

 ns i=1

(4)

1343

si +

nb

 bi

 nf  Ini

+ k3

i=1

(7)

i=1

where nb: total number of DC sources of the electric system. The first term on the right hand side in (7) shows whether a solution covers or not the informed alarms for a given fault situation. This term is zero when the solution totally covers the informed

1344

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

Fig. 1. Basic diagram exemplifying the adopted protection philosophy.

alarms. Otherwise it provides the proximity that a solution has in order to cover the informed alarm. The second term on the right hand side in (7) shows the inconsistency with the informed alarms and the expected states of the protective functions. Therefore Ini is calculated by: Ini = |fsci − fi |

i = 1, . . . , nf

(8)

An alarm is totally consistent with the expected state when Ini = 0. When all alarms of a given fault situation are consistent with the expected state, then the second term in (7) is zero. The closer this term is to zero the more consistent the alarms are with the expected states. For each one of the nf alarms, Coi and Ini are obtained using the rules presented in Table 1. The third term on the right hand side of (7) is included in the objective function to comply with the third condition established by the parsimonious set covering theory. Constants k1, k2 and k3 are integers, and the corresponding values must satisfy k1 > k2 > k3. Since this is a minimization problem, the values for k1, k2 and k3 are chosen in a way that any solution will be discarded if, first, it does not cover the informed alarms; or, second, it is inconsistent with the informed alarms [15,16]. 3. Immune algorithm The optimization IA proposed in this work to solve the UBP model is based on the clonal selection principle [23,24]. The relay and the state of the circuit breaker informed by the SCADA system for each fault situation are considered as genes of an antigen presented to the IA in each generation. Each B cell represents one solution for the problem and is formed by a binary sequence of genes that represents the state of each

section and the state of each DC source of the electric system. The codification used for the B cells is presented in Fig. 2. Minimizing (7) means to find the B cell that complies with the conditions (1), (2) and (3) established by the parsimonious set covering theory. The quality of each B cell is measured by the value of (7). The IA is initialized with a B cell population randomly generated and this population evolves according to the proposed immune operators. In each generation the best clone produced, i.e., the one that presents the minimum value of the objective function is chosen to compose the memory set. The memory set is composed by the best solutions for the problem and the quantity of memory B cells is maintained constant. Moreover, memory B cells are not submitted to the hypermutation and replacement operators. The IA is concluded after a maximum number of generations and the fault section estimation is given by the memory B cell with the minimum value of E. 3.1. Immune operators In order to efficiently solve the fault section estimation problem, the immune operators are proposed and described as follows. 3.1.1. Selection and cloning operators Selection operator consists in selecting n B cells that present the lowest objective function values for cloning (elitist selection). The objective functions referred to each B cell are classified in ascending order. After selection, the cloning operator [23] clones the n B

Table 1 Rule for obtaining Coi and Ini . fsci

fi

Coi

Ini

0 0 1 1

0 1 0 1

0 0 1 0

0 1 1 0

Fig. 2. Codification of B cells. Example of a hypermacromutation applied to interval [i,j].

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

1345

Fig. 3. Calculation of the relay expected states and objective function.

cells proportionally to their objective functions. The quantity of the clones generated is obtained by the following equation: nc =

n

 round

i=1

ˇ × na i

 (9)

where nc: total number of generated clones; na: total number of the B cells in the population; n: quantity of B cells selected for cloning; ˇ: multiplication factor; round (·): operator that approaches the argument to the nearest integer. The proportional characteristic of the cloning operator allows B cells with the lowest objective function values to survive for a long time on the current population, reducing dynamically the diversity of the population during the evolution process. A rigorous diversification strategy of the memory set and the replacement of the B cells in each generation of the algorithm are proposed to adequately solve this feature. The maximum diversity in the memory set is performed through the hypermacromutation operator proposed in [25]. It maintains the memory set with a high percentage of diversification preserving the quantity and quality of the cells. These strategies contribute to the convergence of the algorithm accurately and with less computational time. Furthermore these strategies explore the search space finding multiple quality solutions (suboptimal solutions), and avoid the premature convergence of the algorithm. 3.1.2. Somatic hypermutation operator This operator is applied to the affinity maturation. The hypermutation consists in the mutation or modification of 0–1 or vice versa of one bit of each generated clone. This process is performed using the mutation rate directly proportional to the value of the objective function of the B cell that generates the clone. The mutation rate is defined by: tmi = 1 −

1 e(Ei /(k2×nf +k3(ns+nb)))

i = 1, . . . , n

(10)

where : damping constant. The function of the somatic hypermutation process is to explore local optima by small genetic variations on B cells. These variations can result in B cells that present better quality solutions for the problem. 3.1.3. Hypermacromutation operator It consists in maintaining the maximum diversity among cells that compose the memory set. In each generation of the algorithm, every memory cell is compared to each other by calculating the Hamming distance between two cells. The Hamming distance measures the similarity between two cells, being zero when the cells are identical, and ns + nb when the cells are completely different. The Hamming distance between two B cells, i and j, can be described

as:

  gene k gene k  − B cellj B celli 

ns+nb

Di,j =

(11)

k=1

The identical B cells are submitted to the process called hypermacromutation [25]. This process is similar to the hypermutation; however it is random and independent of the objective function value. The hypermacromutation causes an aggressive genetic modification, searching for B cells that are different from those located in the memory set, but with sufficient quality to remain in this set. The hypermacromutation process consists in randomly choosing two integer numbers [i,j] such as (i + 1) ≤ j ≤ , where  = ns + nb. The memory B cell under hypermacromutation receives T = j − i + 1 tentative of mutation from i to j. This process is shown in Fig. 2. At least one memory B cell of each group of equal cells must be preserved in the maintenance of the diversity process. Therefore, the memory B cells that can potentially lead to optimal solutions are not lost. 3.1.4. Replacement operator This operator seeks to maintain the diversity of the population. This strategy consists in replacing the d B cells with the highest objective functions by randomly generated cells in each generation of the IA. This process is important in the algorithm evolution since each generation allows those solutions with inferior quality to be substituted by better quality solutions. 3.2. Computational procedures In each IA generation the objective function value must be calculated for each B cell of the population. From nf equations established by (1) in an off-line manner, the procedure to calculate the nf relay expected states and the objective function in each generation is as follows: (1) The states of the nf equations established are calculated for each B cell of the population inserting the sections and DC source states (sj and bh ), alarms of the protective relay function (fsct ) and circuit breaker alarms (csck ) into the nf equations generated by (1); (2) With every nf estimated equation the Ini are calculated by using (8), Coi from Table 1 and finally function E is evaluated by applying (7). Fig. 3 illustrates the aforementioned procedure for jth B cell and a set of SCADA alarms (fault situation). The IA proposed and implemented is detailed as follows:

1346

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

Fig. 4. (a) IA proposed and (b) GA proposed in [16].

(i) Randomly generate an initial population of na B cells; (ii) Calculate the expected states of the protective relay functions by using (1) for each B cell of the current population; (iii) Calculate the objective functions of every B cell of the population by applying (7) and classify them in ascending order; (iv) Select n B cells, operators (3.1.1), that correspond to the least objective functions to proliferate by the cloning process (9), generating a population of nc clones; (v) The nc clones are submitted to the affinity maturation process, immune operator (3.1.2), generating a new population of nc maturated clones; (vi) Calculate the expected states of the protective relay functions by using (1) for each of the maturated clones; (vii) Calculate the objective functions for the maturated clones by applying (7) and classify them in ascending order; (viii) From the population of nc maturated clones, reselect the one with a minimum objective function value to be a candidate clone to pertain to the set of m memory B cells. If the objective function value of this maturated clone is less than the objective function of the worst memory B cell, then the clone substitutes the memory B cell in the memory set; (ix) Eliminate the equal configurations among those memory B cells by applying the immune operator (3.1.3); (x) Replace d B cells of the current population with new cells using the immune operator (3.1.4); (xi) Increase the number of q generations. If q is greater than the maximum number of ng generations then stop. The fault section estimation is provided by the memory B cell with the minimum objective function among those of the memory set. Otherwise, go to step ii.

Fig. 4 illustrates the flow chart of the proposed IA as well as the GA proposed in [16]. The GA uses variables crossing over and mutation rates and tournament selection where there are performed np games, with np as the population size. In each game 2 different configurations are selected randomically and is chosen that with the smaller objective function. In addition GA uses diversification rate calculated in each GA generational cycle after the selection. If this rate is below a given pre-established value, the mutation rate value is increased in order to provide new populations without saturation, to maintain diversity and to explore new searching spaces. The variables of Fig. 4(b) are: tc(q) : variable crossing over rate, tv(q) : variable mutation rate, Div: percentage diversification rate, Ceq: maximum number of equal configurations, kc and kv are constants and control the crossing over rate decrement and mutation rate increment, respectively. 4. Results and discussion The proposed methodology to solve the fault section estimation problem is implemented in C++ programming language and the simulations performed using a PC with Intel Core 2, 1.18 GHz and 2 Gb RAM. Part of the South-Brazilian system [6], with nf = 99, ncb = 55, ns = 41 and nb = 5, illustrated in Fig. 5, is used to evaluate the performance of the IA and compare it with a GA (Fig. 4). The control parameters of the IA are calibrated in such a way that the methodology is able to correctly estimate the sections under fault for every simulation (first column) presented in Table 2. Thus, as with any evolutionary algorithm, the calibration of the control parameters of the IA is performed by exhaustive tests and the objective is to assure the robustness and computational efficiency of the

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

Fig. 5. Part of the 230 kV South-Brazilian electric power system.

1347

1348

Table 2 Simulations and solutions—part of the 230 kV South-Brazilian electric power system. Simulation

Input alarms from SCADA system Tripped protective relay functions

Opened circuit breakers

Fault section estimation (providing by proposed methodology) Solution

Busbar(s)

Line(s)

Transformer(s)

DC source(s)

CB(s) failure

CB(s)

DP3

6

1(*)

C

–

–

JV

3, 8, 9

2

DP3, BP11, BP12

6, 9, 11, 12

1(*)

C

L5

–

–

3

BP11, BP12

6, 9, 11, 12

1 (*) 2 (*) 3

C – –

L5 L3, L5 L4, L5

– – –

CT CT CT

4

DP5, DP13, BP8, BP19, OR2, OR5 MP8, MP11, BP6, BP16, BP19, OR1, OR2, OR3, OR4

1(*) 8, 12, 13, 19, 28, 30, 32, 33 6, 8, 11, 16, 19, 28, 29, 1(*) 2 30, 32 3 4 5 6 7 8 9

E

–

TF1, TF2

E – – – – – – – –

L4 L4, L5 L4, L5 L4, L5 L4, L6 L4, L6 L4, L7 L4, L7 L4, L7

6

DP8, MP14, BP16, BP19, BP21, BP25, OR1, OR2, OR3, OR9, OR10

16, 19, 21, 25, 28, 29, 30, 37, 38, 41, 42

1(*) 2 3 4 5 6

E, H E, H H H H H

7

DP6, DP11, DP18, DP20, MP20, MP27, BP14, BP20, BP21

1(*)

8

DP22, MP22, BP21, BP25, OR9, OR10, OR11, OR12

14, 16, 17, 21, 24, 26, 27, 37, 38, 40, 41, 43, 53 21, 22, 25, 37, 38, 41, 42, 49, 51

9

DP10, DP22, MP20, MP21, MP22, OR19 DP13, MP6, MP9, BP6, BP9, BP12, BP16, BP19, OR2, OR5

5

10

Time × 10−2 s IA

GA

2

9

14

3, 8

32

11

16

6, 9 6, 9 6, 9

13 13 23

9 9 6

16 16 –

–

11, 14, 29

13

11

14

TF1 TF1, TF3, TF4 TF1, TF2, TF4 TF1, TF3, TF5 TF1, TF3, TF5 TF1, TF2, TF5 TF1, TF2, TF4 TF1, TF3, TF4 TF1, TF2, TF5

– – – – – – – – –

– – – – – – – – –

23 85 85 85 85 85 85 85 85

11 8 11 6 6 6 6 11 11

16 – – – – – – – –

L7 L7 L7 L7 L7 L7

– TF5 TF1, TF2, TF4 TF1, TF3, TF4 TF1, TF3, TF5 TF1, TF2, TF5

CT CT CT CT CT CT

14, 20, 40 14, 20, 40 14, 20, 40 14, 20, 40 14, 20, 40 14, 20, 40

34 55 86 86 86 86

11 6 11 11 11 6

16 – – – – –

F, K

L9, L14

TF6, TF8

–

20, 45

26

11

11

1(*) 2(*)

F H

L10 L10

TF10 TF10

JV JV

– –

14 14

11 11

16 16

20, 52

1(*)

J

L9, L10

TF10

–

21, 22, 49, 50, 51

54

11

16

9, 16, 19, 28, 29, 32, 33

1(*)

–

L2, L5

TF1, TF2

–

6, 12

1054

11

16

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

Relay functions(s) 1

Emin

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

1349

transformers TF6 and TF8 are under fault conditions and circuit breakers CB20 and CB45 failed. 4.2. Case study: simulation 6

Fig. 6. Performance of the immune and genetic algorithms – simulation 3, solution 2.

algorithm. The control parameters of the IA are: na = 100, ng = 50, n = 5, ˇ = 0.88, m = 10, d = 20,  = 100, k1 = 1000, k2 = 10 and k3 = 1. The parameter ˇ is attributed in such a way that the total number of generated clones is twice the number of B cells in the population, i.e., nc = 200. The control parameters of the GA are: np = 200, ng = 100, kc = 0.30, tc(0) = 0.75, kv = 0.03, tv(0) = 0.001, tp = 0.5 and Divx = 25. The solutions (fourth column) with the symbol (*) in Table 2 are obtained by both IA and GA and are optimal solutions for each fault simulation or fault scenario. The solutions that are not indicated with (*) are found only by the IA and are suboptimal solutions, therefore feasible for the problem. The suboptimal solutions are found only by the IA due to the efficiency of the immune operators proposed, which better exploits the search space. Suboptimal solutions can be interpreted as solutions that contain the sections that present less probability, when compared to the optimal one, of being under fault for a given set of alarms. Considering only optimal solutions, Table 2 shows that the IA has less processing time than the GA to obtain the correct solution for every simulation. An exception is simulation 7, where the time is the same for both solution techniques. Fig. 6 illustrates the performance of the algorithms for simulation 3, solution 2. The IA converges to the optimal solution on generation 35, whereas the GA on generation 79. The diagnosis of the state of the circuit breakers is performed after implementing the algorithms and is based on the following rule: a failure only occurs in one circuit breaker if it does not open after receiving the trip signal of the relay or if the circuit breaker opens without the trip signal of the relay. 4.1. Case study: simulation 7 The analysis of the fault diagnosis for simulation 7 allows concluding it is correct, according to the following explanation. DP6 trip is due to a fault at busbar F. Also DP6 sends a trip to circuit breakers CB16, CB17, CB37, CB38 and CB40 in order to isolate the fault. DP11 trip is due to a fault at busbar K and circuit breakers CB24, CB26, CB53 and CB27 are tripped in order to isolate the fault. Furthermore, MP27 trip is due to a fault at line L14 and this protective function sends a trip to CB27 in order to isolate the fault. Differential functions DP18 and DP20 trip due to a fault at transformers TF6 and TF8, respectively. DP18 sends a trip to CB41 and CB43. DP20 sends a trip to CB37 and CB45. It is observed that CB45 does not open, thus it is expected that there is a failure in this circuit breaker. Distance function MP20 sends a trip to CB20 in order to isolate the fault at line L9. Since CB20 does not open due to failure, local backup protection BP20 and remote backup protection BP14 were tripped. BP14 sends a trip to CB14 in order to isolate the fault propagation to more adjacent lines. Also BP21 sends a trip to CB21 in order to isolate the fault at line L9. Therefore, based on this explanation, it can be concluded that busbars F and K, lines L9 and L14,

In this case there are six different solutions presented in Table 2: one optimal (minimal) solution 1 and five suboptimal solutions. Based on Fig. 5 and the protection philosophy adopted, solution 1 is the best solution because its analysis enables concluding that is the most reliable and safe one to represent the possible status for the alarms informed by the SCADA system. The analysis of fault diagnosis is as follows. A DP8 trip is due to a fault at busbar H. DP8 sends a trip to CB19, CB20, CB40, CB41 and CB42. Since the circuit breakers CB20 and CB40 did not open, there are failures in these equipments. Due to CB20 failure, the busbar H backup protection BP21 sends a trip to CB21. Also, functions OR9, OR10 and BP25 send a trip to CB37, CB38 and CB25, respectively, in order to isolate the fault at busbar H. Distance functions MP14 and BP16 trip are due to fault at line L7. Since CB14 does not open, there is a failure in this circuit breaker. Since protective functions OR1, OR2 and OR3 are the backup protection of busbar E, transformers TF1 (OR1), TF2–TF3 (OR2), TF4–TF5 (OR3), and the main protection of these sections (DP5, DP13, DP14, DP15, DP16, DP17) did not trip, it is hard to point out precisely which of these sections are under fault. From an operator’s point of view it is reasonable to assume that busbar E presents more probability of being under fault than TF1–TF5 due to overcurrent functions OR1–OR3 to protect busbar E and these functions tripped simultaneously. Thus, busbar E is diagnosed under fault in solutions 1 and 2. Distance function BP16 trips is associated to both busbar E (CB14 failure) and line L7 (CB14 failure and CB13 does not open). Because the circuit breakers CB11 and CB12 did not open the backup protections BP8 and BP6, which protect busbar E, respectively, they should have tripped at the substation CT. Since backup protections BP8 and BP6 did not trip and any other protection in substation CT did not, there is DC source failure at the substation CT. The solutions 2–6 result in other possible diagnoses but with small probabilities (suboptimal solutions). Analogous reasoning used to explain solution 1 can be used for other solutions. It is important to notice that the suboptimal solutions indicate to the operator that these sections under fault are difficult to occur but they should be checked. In a practical case the simulation 6 solutions would indicate to the operator to check busbars E, H, line L7, transformers TF1, TF2, TF3, TF4 and TF5, DC source located at substation CT and circuit breakers CB14, CB20 and CB40. 4.3. Impact of noisy data and effects of system uncertainties From the point of view of the fault section estimation using alarms from relays and circuit breakers, noisy data can be considered such as: redundant alarms, false alarms, lack of alarms as communication failure between SCADA system and protective devices or remote terminal units (RTU). Redundant alarms can be filtered using an alarm processor or other equivalent technique. In any way redundant alarms are not a serious problem as compared to false alarms and lack of alarms. False alarms occur when the protection system does not operate and the SCADA system (due to noisy or corrupted alarms) indicates that determined alarm occurs. In this case the fault section estimation methodology can provide an incorrect result to the operator. However in these cases the solution will have more sections than the correct solution. Therefore the operator can check all sections and scrap the non-defective sections. Lack of alarms can be a major problem because the system protection has been operated and the alarms are not transmitted to the methodology’s input. The diagnosis can be incorrect in the

1350

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

(a) The methodology does not use knowledge base; therefore it does not require the availability of an expert human in order to construct the knowledge base. The equations of protective relay functions are based on operation logic and protection philosophy and do not suffer modifications when the system topology changes. Also the binary basic matrices do not change with system topology. This feature clearly present benefits because the methodology is able to generalize and is independent of fault situation. In addition the basic matrices are scalable and can be enlarged when there is an expansion of the system. (b) The methodology does not require validation of the mathematical model. This is especially important because rare fault situation can occur and it is very difficult to validate the ES’s knowledge base in order to diagnose this kind of fault precisely. When compared to ANN techniques the proposed methodology presents the benefits:

Fig. 7. Flow chart of proposed methodology implementation.

sense that fewer sections are diagnosed under fault. Therefore it is important that the operator can access the raw data and the SCADA system can inform the operator if any alarms are lacked. Modern SCADA systems have techniques to minimize the hazard of data contamination, such as, cyclic redundancy check (CRC). It is important notice that the quality of diagnosis depends on the quality of the alarms independently of which methodology is implemented. As a basic idea Fig. 7 illustrates a possible implementation from methodology proposed. The protection system uncertainties can be enumerated such as: operation of relay along with circuit breaks failure, failure of relays, or failure of relays and circuit breakers simultaneously. The effects of protection system uncertainties reflect directly to the multimodal characteristic of (7) and thus the existence of several suboptimal solutions for certain simulations (Table 2). Circuit breaker failure along with the main protection failure and backup protection operations make the search process difficult because backup protection usually covers more than one section. In case of relay failure the search space is increased and other sections that do not pertain to the minimal solution can be estimated under fault, i.e., suboptimal solutions. These solutions result from inconsistency with informed alarms and expected states of the protective functions that increase the second term on the right hand side in (7). Thus the expected states of the protective function that protect these sections are calculated to be 1 because the protective function should be operated. However, due to relay failure the alarm is 0, i.e., non-actuated (second row of Table 1) and there will be an inconsistency. Moreover, the complexity of the fault scenarios depends on the number of protective devices alarms and system topology at fault moment. Many alarms can generate one solution with many sections under fault. On the other hand a number of alarms inadequate to deduce an accurate diagnosis can result in several suboptimal solutions for the fault scenario.

4.4. Benefits and contributions of the proposed methodology The main drawbacks of the methodologies employed up to now in the literature are overcome by the proposed methodology in this work. When compared to ES the methodology presents the benefits:

(a) It does not require the training process. Training process requires great processing time when ANN is developed to large electric power systems. In these cases, the quantity of alarms sample set in order to ANN diagnosis any fault situation precisely is difficult to determine. Therefore the proposed methodology does not depend on any training or any sample in order to diagnose faults. Instead the control parameters of the IA are necessary to tune. Furthermore the control parameters do not suffer any modification with changes in system topology. (b) ANN requires determining of dimension of the network such as hidden units and quantity of layers. The proposed methodology does not require any dimensioning. Large electric power systems require large ANNs. This fact becomes difficult to determine the optimal weights of the neural network. The mathematical model of protection system operations developed in this work models the operation logic of the protective relay functions in a mathematical manner. Different from FL the operation logic of the protection system is not modeled through membership functions. In order to obtain the membership functions data is necessary from electric power system, such as, historic data from relays and circuit breakers operations. Also membership functions can be obtained through operator or engineer experience, or trial-and-error. Considering that the protection system operation presents uncertainty it is very difficult develop robust membership functions that are able to model the protective devices operation properly. The proposed methodology does not require any data in order to model the protection system. The only data is the basic data stored in binary basic matrices. Basically PN uses graphical representation and when applied to large electric power system it is very difficult to develop and to manage. Otherwise the proposed methodology does not use any graphical representation and the equations are simple and can be implemented in any computational language. In addition the mathematical model demand soft computation and can be implemented using addition and subtractions operations. Furthermore the proposed mathematical model is independent of solution technique and can be solved by any suitable solution technique (flexibility). Analyzing the algorithm proposed in [26] the methodology proposed in this work has the follows contributions: (a) The binary basic matrices are developed and an algorithm to obtain the sets ˝ti,j and ˝li,j (Appendix A) is proposed; (b) The fitness function is improved and the results are better than [26]. (c) The IA proposed has better performance and is more robust than [26]. This superiority is due to the usage of the hypermacromutation immune operator. Also the mutation rate used in [26]

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

is directly proportional to the fitness function and the results show that Eq. (10) is more powerful to solve the fault section estimation problem. 5. Conclusions This paper proposes a new methodology based on an unconstrained binary programming model and optimization immune algorithm to solve the fault section estimation problem. A real-life system is used in order to demonstrate the methodology’s applicability. The results show that the proposed methodology presents efficiency, precision, robustness, computational speed, and is adequate to estimate faults from alarm reports of simple, complex or misoperation situations of the protective devices. The immune algorithm is superior to the genetic algorithm proposed earlier, in terms of the capability of obtaining suboptimal solutions as well as the processing time required for the fault section estimation. This superiority is due to the proposed immune operators and the adequate calibration of the control parameters of the algorithm. From the operator’s point of view, suboptimal solutions are important, because they provide a more complete and in depth fault estimation, in addition to enabling a better analysis of the fault situation. The results show that the proposed methodology presents the potential for real-time fault estimation in electric power system control centers, hence overcoming certain drawbacks of earlier proposed techniques. Acknowledgments The authors acknowledge the financial support from FAPESP under Grant 06/02569-7 and CNPq under Grant 301060/2006-1 and Dr. A.D.P. Lotufo (Department of Electrical Engineering – Engineering College of Ilha Solteira – São Paulo State University – UNESP) for her collaboration.

1351

established according to the following Boolean function: r2(im, ib, p) = mc(im,p) • mc(ib,p) im = 1, . . . , nf ; ib = 1, . . . , nf ; p = 1, . . . , ncb

(A.2)

If r2(im,ib,p) = 1 then the functions im and ib trip the same circuit breaker p. If r2(im,ib,p) = 0 then there are three situations: (1) mc(im,p) = 1 and mc(ib,p) = 0, that means that only the main protection function im trips the circuit breaker p; (2) mc(im,p) = 0 and mc(ib,p) = 1, that means that only the backup protection function ib trips the circuit breaker p; and (3) mc(im,p) = 0 and mc(ib,p) = 0, that means that the circuit breaker p is not tripped by functions im and ib. From the Boolean functions (A.1) and (A.2) and considering that the indexes im, ib, j and p have the same value for both rules, the sets ˝ti,j can be obtained considering the following Boolean function:



r1(im, ib, j) • r2(im, ib, p) =

1 ⇒ ˝tib,j = {im}. 0 ⇒ ˝tib,j = ∅ ⇒ fsct = 1.

(A.3)

Eq. (A.3) shows that the set ˝ti,j are not empty only if the conditions given by rules 1 and 2 are true, i.e., equal to 1 simultaneously. This statement is according to the definition of these sets. A.2. Set ˝li,j The formation of the set ˝li,j is based on four rules established by Boolean functions considering that for each set ˝lib,je the circuit path must be examined between the circuit breaker tripped by the backup function ib and the protected section je. For this, it is necessary to identify: (1) which section je is protected by function ib; (2) which circuit breaker is tripped by function ib and (3) which section (or sections) is connected to the tripped circuit breaker. Rule 1 verifies if the protection ib is the backup protection of section je. Therefore: r1(ib, je) = mb(ib,je) ib = 1, . . . , nf ; je = 1, . . . , ns

(A.4)

Appendix A.

If r1(ib,je) = 1 then the backup function ib protects section je and rule 2 must be checked to identify which circuit breaker function ib trips. Rule 2 can be established by the following Boolean function:

For the elements (i,j) of matrices M, Mm, Mb, Mc and Md given by m(i,j) , mm(i,j) , mb(i,j) , mc(i,j) and md(i,j) , respectively, the following rules for forming the sets ˝ti,j and ˝li,j are established:

r2(ib, ps) = mc(ib,ps) ib = 1, . . . , nf ; ps = 1, . . . , ncb

A.1. Set ˝ti,j The formation of the set ˝ti,j is based on two rules that are checked separately. Rule 1 checks if the main protection function im protects the same section j as the backup protection function ib. Mathematically, the following Boolean function can be written for rule 1: r1(im, ib, j) = mm(im,j) • mb(ib,j) im = 1, . . . , nf ; ib = 1, . . . , nf ; j = 1, . . . , ns

(A.1)

where •: Boolean conjunction operator – logic AND. When r1(im,ib,j) = 1, then both main im and backup ib functions protect the section j. Otherwise, if r1(im,ib,j) = 0 then there are three situations: (1) mm(im,j) = 0 and mb(ib,j) = 1, that means that function ib is the backup protection of section j; (2) mm(im,j) = 1 and mb(ib,j) = 0, that means that function im is the main protection of section j, and (3) mm(im,j) = 0 and mb(ib,j) = 0 and therefore m(im,j) = m(ib,j) = 0 that means that section j is not protected by the considered functions and that sets Jim and Jib do not consider the index of the section as an element. If r1(im,ib,j) = 1 for section j, then according to rule 2 it must be checked if the functions im and ib besides protecting the same section, also trip the same circuit breaker. Rule 2 is mathematically

(A.5)

If r2(ib,ps) = 1 then the function ib trips the circuit breaker ps. When the index ps is obtained, then rule 3 must verify which section js is connected to circuit breaker ps. Section js is the initial section of the circuit path, from the circuit breaker ps, that contains the protected section je. The Boolean function that describes rule 3 can be established as: r3(ib, ps, js) = mb(ib,js) • md(ps,js) ib = 1, . . . , nf ; ps = 1, . . . , ncb; js = 1, . . . , ns

(A.6)

If r3(ib,ps,js) = 1 then js is the initial section. If there is more than one index js where r3(ib,ps,js) = 1, then js = je, with je obtained by (A.4). From the Boolean functions (A.4)–(A.6) and considering that the indexes ib, je, ps and js have the same values for every rule, the sets ˝lib,je can be obtained considering the following Boolean function: r1(ib, je) • r2(ib, ps) • r3(ib, ps, js)



=

1 ⇒ ˝lib,je = {ps}, with je = js. 1 ⇒ ˝lib,je = {index(es)  → rule 4}, with je =/ js. 0 ⇒ ˝lib,je = ∅ ⇒ csck = 0.

(A.7)

Rule 4 establishes the procedures for obtaining the index(es) of the circuit breaker(s) connected between section js and section je. These procedures use the properties of matrix Md that are:

1352

F.B. Leão et al. / Electric Power Systems Research 80 (2010) 1341–1352

(1) Examining the lines of matrix Md and maintaining the same column (vertical displacement) means to examine the same section in the electric system, from a circuit breaker connected to one of the terminals to the other circuit breaker that is connected to the other section terminal. (2) Examining the columns of matrix Md and maintaining the same line (horizontal displacement) means to examine the same circuit breaker from one section connected to one of the terminals of the circuit breaker to the other section connected to the other terminal. With the indexes ib, je, ps and js estimated by (A.7), with je = / js, rule 4 can be mathematically established by the following procedures: i. Let: ps∗ = ps

(A.8)

js∗ = js

(A.9)

ii. Verify the following Boolean function (vertical displacement): md(ps∗,js∗) • md(p,js∗)

p = 1, . . . , ncb

(A.10)

Go to step (iii) only if (A.10) is equal to 1. iii. Let: ps∗ = p

(A.11)

iv. Verify the following Boolean function (horizontal displacement): md(ps∗,js∗) • md(ps∗,j)

j = 1, . . . , ns

(A.12)

Go to step (v) only if (A.12) is equal to 1. v. Let: js∗ = j

(A.13)

vi. If js* = je then stop. The set ˝lib,je is formed by all indexes p that satisfy (A.10) and (A.12). Otherwise return to step (ii). The symbol (*) contained in the indexes of Eqs. (A.10) and (A.12) means that these indexes are fixed and therefore their values are changed by (A.11) and (A.13). References [1] B.F. Wollenberg, Feasibility study for an energy management system intelligent alarm processor, IEEE Trans. Power Syst. 2 (1986) 241–246. [2] C. Fukui, J. Kawakami, An expert system for fault section estimation using information from protective relays and circuit breakers, IEEE Trans. Power Deliv. 4 (1986) 83–90.

[3] C.A. Protopapas, K.P. Psaltiras, A.V. Machias, An expert system for substation fault diagnosis and alarm processing, IEEE Trans. Power Deliv. 6 (1991) 648–655. [4] H.T. Yang, W.Y. Chang, C.L. Huang, On-line fault diagnosis of power substation using connectionist expert system, IEEE Trans. Power Syst. 10 (1995) 323–331. [5] H.T. Yang, W.Y. Chang, C.L. Huang, A new neural networks approach to online fault section estimation using information of protective relays and circuit breakers, IEEE Trans. Power Deliv. 9 (1994) 220–230. [6] G.C. Junior, J.G. Rolim, H.H. Zürn, Application of neural-network modules to electric power system fault section estimation, IEEE Trans. Power Deliv. 19 (2004) 1034–1041. [7] C. Rodríguez, S. Rementería, J.I. Martín, A. Lafuente, J. Muguerza, J. Pérez, A modular neural network approach to fault diagnosis, IEEE Trans. Neural Netw. 7 (1996) 326–340. [8] J.C.S. de Souza, E.M. Meza, M.Th. Schilling, M.B.C. Filho, Alarm processing in electrical power systems through a neuro-fuzzy approach, IEEE Trans. Power Deliv. 19 (2004) 537–544. [9] H. Monsef, A.M. Ranjbar, S. Jadid, Fuzzy rule based expert system for power system fault diagnosis, IEE Proc. Gener. Transm. Distrib. 144 (1997) 186–192. [10] C.S. Chang, J.M. Chen, D. Srinivasan, F.S. Wen, A.C. Liew, Fuzzy logic approach in power system fault section identification, IEE Proc. Gener. Transm. Distrib. 144 (1997) 406–414. [11] H.C. Chin, Fault section diagnosis of power system using fuzzy logic, IEEE Trans. Power Syst. 18 (2003) 245–250. [12] S.W. Min, J.M. Sohn, J.K. Park, K.H. Kim, Adaptive fault section estimation using matrix representation with fuzzy relations, IEEE Trans. Power Syst. 19 (2004) 842–848. [13] F. Wen, Z. Han, Fault section estimation in power systems using a genetic algorithm, Electr. Power Syst. Res. 34 (1995) 165–172. [14] F.S. Wen, C.S. Chang, Probabilistic approach for fault section estimation in power systems based on a refined genetic algorithm, IEE Proc. Gener. Transm. Distrib. 144 (1997) 160–168. [15] F. Wen, C.S. Chang, A new approach to fault diagnosis in electrical distribution networks using a genetic algorithm, Artif. Intell. Eng. 12 (1998) 69–80. [16] F.B. Leão, L.G.W. da Silva, J.R.S. Mantovani, Defects location in power systems components through a dedicated genetic algorithm, IEEE/PES Transm. Distrib. Conf. Expo. Lat. Am. (2004) 57–62. [17] K.L. Lo, H.S. Ng, J. Trecat, Power systems fault diagnosis using petri nets, IEE Proc. Gen. Transm. Distrib. 144 (1997) 231–236. [18] F. Wen, C.S. Chang, A tabu search approach to fault section estimation in power systems, Electr. Power Syst. Res. 40 (1997) 63–73. [19] C.S. Chang, L. Tian, F.S. Wen, A new approach to fault section estimation in power systems using ant system, Electr. Power Syst. Res. 49 (1999) 63–70. [20] S.J. Huang, Application of immune-based optimization method for faultsection estimation in a distribution system, IEEE Trans. Power Deliv. 17 (2002) 779–784. [21] Y. Peng, J.A. Reggia, Abductive Inference Models for Diagnostic Problem–Solving, Springer-Verlag, New York, 1990. [22] P.M. Anderson, Power System Protection, McGraw-Hill IEEE Press Ser. on Power Eng., New York, 1999. [23] L.N. de Castro, F.J.V. Zuben, Learning and optimization using the clonal selection principle, IEEE Trans. Evol. Comput. 6 (2002) 239–251. [24] D. Dasgupta, N.A. Okine, Immunity-based systems: a survey, in: IEEE Int. Conf. on Syst., Man, and Cybern. 1, 1997, pp. 369–374. [25] V. Cutello, G. Nicosia, M. Pavone, J. Timmis, An immune algorithm for protein structure prediction on lattice models, IEEE Trans. Evol. Comput. 11 (2007) 101–117. [26] F.B. Leão, R.A.F. Pereira, J.R.S. Mantovani, Fault section estimation in electric power systems using an artificial immune system algorithm, Int. J. Innov. Energy Syst. Power 4 (2009) 14–43.