A comparative study of metaheuristic optimization approaches for directional overcurrent relays coordination

Electric Power Systems Research 128 (2015) 39–52

Contents lists available at ScienceDirect

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

A comparative study of metaheuristic optimization approaches for directional overcurrent relays coordination Mahamad Nabab Alam ∗ , Biswarup Das, Vinay Pant Department of Electrical Engineering, Indian Institute of Technology, Roorkee, 247667, India

a r t i c l e

i n f o

Article history: Received 5 March 2015 Received in revised form 27 May 2015 Accepted 22 June 2015 Keywords: Directional overcurrent relays Meta-heuristic optimization approach Power system protection Protective relaying

a b s t r a c t In this paper, a comparative study of different meta-heuristic optimization approaches, which had been proposed in the literature for directional overcurrent relay coordination (DOCRs), is presented. Towards this goal, five most effective meta-heuristic optimization approaches such as genetic algorithm (GA), particle swarm optimization (PSO), differential evolution (DE), harmony search (HS) and seeker optimization algorithm (SOA) have been considered. The performances of these optmization methods have been investigated on several power system networks of different sizes. The comparative performances of these methods have been studied by executing each method 100 times with the same initial conditions and based on the obtained results, the best meta-heuristic optimization method for solving the coordination problem of directional overcurrent relays is identified. © 2015 Elsevier B.V. All rights reserved.

1. Introduction A properly coordinated protection scheme is one of the inherent requirements to operate a power system with highest reliability. A good protection scheme removes only the least possible portion of the system whenever a fault occurs so as to maintain supply to the rest of the healthy system unaffected by the fault. Each equipment of a power system is protected with two lines of defence, which are known as primary protection and backup protection. For reliable operation of the system, primary protection must react for a fault as quickly as possible to isolate the faulty parts from the healthy parts, but if primary protection fails to operate, the backup protection should operate. This condition is the most desired feature of any protection scheme as primary protection removes only faulted part whereas, whenever backup protection operates, a larger portion of the system has to suffer from outage unnecessarily. For ensuring that only the faulted portion of the network is disconnected thereby reducing the possibility of unwanted power outage, proper co-ordination among the protective devices is necessary. An economic and effective protection scheme for meshed or multi-sourced power systems requires directional overcurrent relays (DOCRs). The operation of DOCRs depends on its two parameters, namely time multiplier setting (TMS) and plug setting (PS).

∗ Corresponding author. Tel.: +91 9897998728. E-mail address: [email protected] (M.N. Alam). http://dx.doi.org/10.1016/j.epsr.2015.06.018 0378-7796/© 2015 Elsevier B.V. All rights reserved.

Through co-ordination, the proper TMS and PS of the relays are determined such that any fault is cleared by the corresponding primary relay as soon as possible. Also, both these settings of any relay should be properly coordinated with the relays protecting the adjacent equipments which, in turn, makes the co-ordination problem quite complex. To solve this complex problem, several methods have been developed in the literature. These methods include curve intersection approach [1], application of a graphical selection procedure for selecting the settings of the relays [2], identification of minimum break point set (MBPS) using expert system [3], linear graph theory [4], simplex method [5,6], Gauss–Seidel iterative procedure [7], and sequential quadratic programming (SQP) [8], etc. Recently, metaheuristic optimization approaches have shown great potential to solve protection coordination problems. Application of various meta-heuristic optimization approaches such as particle swarm optimization (PSO) [9,10], genetic algorithm (GA) [11–13], differential evolution (DE) [14], harmony search (HS) [15], seeker optimization algorithm (SOA) [16], etc., have also been proposed in the literature. The applications of these algorithms have been demonstrated on several small to medium sized systems in all these works. Although these approaches are somewhat time-consuming but they provide quite high quality solutions. Now, because of the availability of several meta-heuristic optimization methods for co-ordination of DOCRs, it is a natural curiosity to find the most effective meta-heuristic optimization method for practical implementation. However, in the literature,

40

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

Nomenclature individuals of population i ∈ {1, 2, . . ., N} components of an individual j ∈ {1, 2, . . ., D} iteration counter (k ∈ {1, 2, . . ., Maxite}) crossover factor mutation factor for genetic algorithm crossover rate mutation factor for differential evolution population size component of each individual of population maximum number of iterations index of the best individual in population index of the worst individual in population inertia factor minimum value of inertia factor maximum value of inertia factor bandwidth for harmony search minimum value of bandwidth maximum value of bandwidth pitch adjustment rate of harmony minimum value of pitch adjustment rate maximum value of pitch adjustment rate harmony memory consideration rate minimum value of membership degree maximum value of membership degree acceleration factor bc ∈ {bc1 , bc2 , bc3 } is an index of best seeker of first, second and third 1/3rd of population abs(·) return an absolute value of an input number f(·) objective function to be evaluated signum function on each variable of the input vector sign(·) rand() uniformly generated random number in the range [0, 1] randj uniformly generated random numbers in [0, 1] for jth component of an individual randi(N) return a uniformly generated random integer such that randi(N) ∈ {1, 2, . . ., N} RAND(ki , 1) uniformly generate random number in the

i j k CF MF CR F N D Maxite b c w wmin wmax BW BWmin BWmax PAR PARmin PARmax HMCR min max c1 and c2 bc

range [ki , 1] population of N individuals each having D components (variables) initial velocity of N individuals each having D components ith individual of population X at iteration k, i.e., Xki = k , X k , . . ., X k ] [Xi,1 i,2 i,D

X V Xki k Xi,j k Vi,j

Fik TFi TXi,j Ui,j Pbestki k

Pbest i,j Gbestk

jth component of ith individual of population at iteration k velocity of jth component of ith individual of population at iteration k value of objective function for ith individual of population at iteration k value of objective function for trial solution i at iteration k jth component of ith individual of trial solution jth component of ith individual of trial solution obtained after crossover operation personal best of ith individual of population up to iteration k personal best jth component of ith individual of population up to iteration k the global best individual of population up to iteration k

k

jth component of the best individual of population up to iteration k jth component of new harmony vector minimum value of jth component of an individual of population maximum value of jth component of an individual of population step length of jth component of a seeker at iteration k membership degree of ith seeker at iteration k

Lbestk k Lbest j

the local best seeker of population at iteration k jth component of local best seeker of population at iteration k jth component of the best seeker of a subpopulation bc ∈ {bc1 , bc2 , bc3 } at iteration k jth component of egotistic direction of ith seeker jth component of global altruistic direction of ith seeker jth component of local altruistic direction of ith seeker jth component of proactiveness direction of ith seeker jth component of step size of ith seeker at iteration k

Gbest j NHVj Xmin,j Xmax,j ıkj ki

k

Xbest bc,j di,j,ego di,j,alt1 di,j,alt2 di,j,pro ˛ki,j

no such comprehensive comparative study is available. Motivated by this fact, this paper endeavours to conduct this study for finding out the most effective method. This paper is organized as follows. In Section 2, the co-ordination problem of the DOCRs is described. In Section 3, the detailed algorithms of different meta-heuristic methods investigated in this work are described. Lastly, Section 4 gives the main results of this work and recommends the most effective method.

2. Problem formulation of protection coordination The protection coordination problem can be formulated as an optimization problem where the objective is to minimize the sum of the operating times of all the numerical DOCRs for the near-end three phase fault current [5,17]. Therefore, the objective function (OF) is expressed as,

min

n 

top,l

(1)

l=1

In Eq. (1), n is the number of relays in the system and top,l is the operating time of the relay Rl . The operating times of the relays are obtained from their characteristic curves which are defined by IEC/IEEE [18] as, top =

 × TMS 

(IF /PS) − 1

+L

(2)

In Eq. (2), ,  and L are the characteristic constants of the relays while IF is the fault current through the relay operating coil. For standard inverse definite minimum time (IDMT) relays  = 0.14,  = 0.02, and L = 0 [18]. Whereas, for other types of relays like very inverse (VI) relays  = 13.5,  = 1, and L = 0 and for extremely inverse (EI) relays  = 80,  = 2, and L = 0 [18]. The objective function defined above is subjected to the following sets of constraints [5,6]:

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

(1) Limitations on TMS, PS and top of relays: The limits on TMS, PS and top of relays are expressed as TMS l,min ≤ TMS l ≤ TMS l,max

(3)

PS l,min ≤ PS l ≤ PS l,max

(4)

tl,min ≤ top,l ≤ tl,max

(5)

In Eq. (3), TMSl,min and TMSl,max , in Eq. (4), PSl,min and PSl,max and in Eq. (5) tl,min and tl,max are the minimum and maximum limits of TMS, PS and operating time (top ) of relay Rl , respectively. Now, from Eqs. (2)–(5), it may appear that if Eqs. (3) and (4) are satisfied, then Eq. (5) is redundant. However, this is not so, as the operating time of the relay depends not only on TMS and PS but also on the fault current (IF ) according to Eq. (2). Even if TMS and PS are in the correct range, the operating time of the relay may not be in the correct range. For example, let us suppose that TMS and PS of a relay is 0.6 and 2.0, respectively, which are well within their respective ranges (these ranges are given in Section 4). Now, let us suppose that a fault occurs near the relay and causes 1600 A current to flow to the fault point. If current transformer (CT) ratio is 500:1, then 3.2 A current will flow through the secondary of the transformer and hence through the relay. Now, operating time of the relay (according to Eq. (2) for IDMT type) in this situation will be 8.8942 s (0.14 × 0.6/((3.2/2)0.02 − 1)) which is very large and not in the range of [0.1, 4.0] s (as mentioned in Section 4). (2) Protection coordination criteria: If primary relay Rl has a backup relay Rm for a fault at any point k, then the corresponding coordination constraint can be expressed as tob,m − top,l ≥ CTI

(6)

In Eq. (6), top,l and tob,m are the operating time of primary relay Rl and its backup relay Rm , respectively, for the same fault and CTI is the coordination time interval required for proper operation of primary/backup relays. 3. Optimization algorithms to solve protection coordination problems Five highly efficient algorithms suitable to solve the protection coordination problems are discussed here concisely. These five algorithms are GA, PSO, DE, HS and SOA. All these methods start from an initial solution to attain the optimum point in the search space. The size of the population is considered as N and the dimension of each element of the population is considered as D, where D represents the total number of variables. Thus, the initial solution is denoted as X = [X1 ,X2 ,. . .,XN ]T , where ‘T’ denotes the transpose operator. Each individual Xi (i = 1, 2, . . ., N) is given as Xi =[Xi,1 , Xi,2 , . . ., Xi,D ]. A protection scheme with n DOCRs will have 2 × n number of variables (i.e., D = 2n). The first n variables represent TMS values, whereas the second n variables represent PS values of the relays in the system. The detailed algorithms of various methods are described below for completeness. In all these algorithms, the index i varies from 1 to N, whereas the index j varies from 1 to D. 3.1. Genetic algorithm Genetic algorithm (GA) is a nature inspired meta-heuristic optimization approach suitable to solve any optimization problem irrespective of the type of objective function. The algorithm starts with a randomly generated solution called chromosomes and endeavors to reach the optimum point by reproduction, crossover and mutation operations. The crossover factor (CF) is expressed as the probability of pairs of chromosomes to produce offsprings (crossover), whereas mutation factor (MF) represents the probability of changing the status of a randomly selected binary bit of a

41

chromosome from 0 to 1 and vice versa (mutation). The different steps of GA are as follows [13]: Set parameter CF and MF of GA Initialize population of chromosomes, i.e., the solution X. Set iteration k = 1 Calculate fitness of chromosomes Fik = f (Xki ), ∀i and find the index of the best chromosome b ∈ {1, 2, . . ., N} 5. Perform selection of chromosomes, crossover of parents and mutation of offsprings to form a new set of chromosomes , ∀i Xk+1 i

1. 2. 3. 4.

), ∀i and identify the best chromo6. Evaluate fitness Fik+1 = f (Xk+1 i some b1 k+1 7. If Fb1 < Fbk then b = b1 8. If k < Maxite then k = k + 1 and goto step 5 else goto step 9 9. Print optimum solution as Xkb 3.2. Particle swarm optimization Particle swarm optimization (PSO) is inspired by social and cooperative behavior displayed by various species to fill their needs in the search space. The algorithm is guided by personal experience (Pbest), overall experience (Gbest) and the present movement of the particles to decide their next positions in the search space. Further, the experiences are accelerated by two factors c1 and c2 , and two random numbers generated between [0, 1], whereas the present movement is multiplied by an inertia factor w varying between [wmin , wmax ]. The initial velocity of the population is denoted as V = [V1 ,V2 ,. . .,VN ]T , where ‘T’ denotes the transpose operator. Thus, the velocity of each particle Xi (i = 1, 2, . . ., N) is given as Vi =[Vi,1 , Vi,2 , . . ., Vi,D ]. The different steps of PSO are as follows [19]: 1. Set parameter wmin , wmax , c1 and c2 of PSO 2. Initialize population of particles having positions X and velocities V 3. Set iteration k = 1 4. Calculate fitness of particles Fik = f (Xki ), ∀i and find the index of the best particle b 5. Select Pbestki = Xki , ∀i and Gbestk = Xkb 6. w = wmax − k × (wmax − wmin )/Maxite 7. Update velocity and position of particles k

k+1 k k Vi,j = w × Vi,j + c1 × rand( ) × (Pbest i,j − Xi,j ) + c2 × rand( ) k

k × (Gbest j − Xi,j ); ∀j and ∀i

k+1 k+1 k Xi,j = Xi,j + Vi,j ; ∀j and ∀i

8. Evaluate fitness Fik+1 = f (Xk+1 ), ∀i and find the index of the best i particle b1 9. Update Pbest of population ∀i = Xk+1 else Pbestk+1 = Pbestki If Fik+1 < Fik then Pbestk+1 i i i 10. Update Gbest of population k+1 If Fb1 < Fbk then Gbestk+1 = Pbestk+1 b1

and set b = b1 else Gbestk+1 = Gbestk

42

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

11. If k < Maxite then k = k + 1 and goto step 6 else goto step 12 12. Print optimum solution as Gbestk

7. Generate a new harmony with randomly generated numbers r1 ∈ [0, 1] and rj ∈ [0, 1], ∀ j

3.3. Differential evolution Differential evolution (DE) is a parallel direct search method guided by crossover, mutation and selection operations on a population. For each element of the population a mutant element is generated by adding the difference between two mutually independent elements multiplied by a mutation factor F. Subsequently, crossover is performed to make trial elements of the same size as the population with a crossover rate (CR). After that, selection is performed to pickup better trial elements. The different steps of DE used in this study are as follows [20]: Set parameter F and CR of DE Initialize population of elements X Set iteration k = 1 Calculate fitness of elements Fik = f (Xki ), ∀i and find the index of the best element b 5. Perform mutation with three mutually different random indices a1 , a2 , a3 ∈ {1, 2, . . ., N} to generate trial solution

1. 2. 3. 4.

TX i,j = Xak

1 ,j

+ F × (Xak

2 ,j

− Xak ,j ) ∀j and ∀i 3

6. Perform crossover on trial solution with the help of the CR and a randomly generated index randb where randb = randi(D) ∈ {1, 2, . . ., D}

 Ui,j =

TX i,j

if rand ( ) ≤ CR or j = randb

k Xi,j

if rand ( ) > CR and j = / randb

7. Evaluate trial fitness TFi = f(Ui,1 , Ui,2 , . . ., Ui,D ), ∀i 8. Perform selection of trial solution ∀i

 k+1 Xi,j =

NHV j =

⎧ Xrandi(N),j + rand( ) × BW, if r1 < HMCR and rj < PAR ⎪ ⎪ ⎨ Xrandi(N),j , if r1 < HMCR and rj > PAR

⎪ ⎪ ⎩ Xmin,j + rand( ) × (Xmax,j − Xmin,j ), otherwise

8. Evaluate trial fitness TF = f(NHV1 , NHV2 , . . ., NHVD ) k = NHV , ∀j elseif TF < F k then X k = 9. If TF < Fbk then Xb,j j c c,j

NHV j , ∀j and find the index of the worst harmony c else goto step 7 10. If k < Maxite then k = k + 1 and goto step 5 else goto step 11 11. Print optimum solution as Xkb 3.5. Seeker optimization algorithm

The seeker optimization algorithm (SOA) is a computational search algorithm inspired by the behaviour of human memory consideration, experience gained, uncertainty reasoning, and social learning. A simple fuzzy rule is used to evaluate seekers step length. The different steps of SOA are as follows [22]: 1. 2. 3. 4.

Set parameter wmin , wmax , min and min of SOA Initialize the solution, i.e., population of seekers X Set iteration k = 1 Calculate fitness of seekers Fik = f (Xki ), ∀i and find the index of the best seeker b and the indices of locally best seekers bc1 , bc2 and bc3 among first, second and third 1/3rd of population, respectively. 5. Select Pbestki = Xki , ∀i, Gbestk = Xkb , Lbestk = Xkb and Xbestkbc = Xkbc , ∀bc where bc ∈ {bc1 , bc2 , bc3 } 6. Calculate search direction ∀j and ∀i k k ) di,j,ego = sign(Pbest i,j − Xi,j k

k ) di,j,alt1 = sign(Gbest j − Xi,j k

k ) di,j,alt2 = sign(Lbest j − Xi,j

Ui,j ∀j if TF i < Fik

di,j,pro =

k ∀j if TF ≥ F k Xi,j i i

where k1 , k2 ∈ {k, k − 1, k − 2} such that f (ik ) < f (ik )

9. Evaluate fitness Fik+1 = f (Xk+1 ), ∀i and find the index of the best i element b 10. If k < Maxite then k = k + 1 and goto step 5 else goto step 11 11. Print optimum solution as Xkb 3.4. Harmony search Harmony search (HS) algorithm is inspired by the creative process of music composition where music players improvise the pitches of their instruments to obtain better harmony. The algorithm is guided by harmony memory consideration rate (HMCR), pitch adjustment rate (PAR) and bandwidth (BW) to obtain improved harmony memory (HM), i.e., solution vectors after each iteration. The different steps of HS are as follows [21]: 1. Set parameter BWmin , BWmax , PARmin , PARmax and HMCR of HS 2. Initialize population of harmonies X 3. Calculate fitness of harmonies Fik = f (Xki ), ∀i and find the index of the best harmony b and the index of the worst harmony c 4. Set iteration k = 1 5. PAR = PARmin + k × (PARmax − PARmin )/Maxite 6. BW = BWmax − k × (BWmax − BWmin )/Maxite

k1 sign(Xi,j

k2 ) − Xi,j

sdi,j = {di,j,ego , di,j,alt1 , di,j,alt2 , di,j,pro }

⎧ ⎨ 0, k+1 di,j

=

⎩

1

2

if randj < (no.ofzerosin sdi,j )/4

+1, elseif randj < (no. of zeros + no. of ones in sdi,j )/4 −1,

else

7. Calculate step size ∀j and ∀i w = wmax − k × (wmax − wmin )/Maxite k+1 = max − (s − Iik+1 ) × (max − min ), ∀i i k+1 1,j

ık+1 = w × abs(Xbest bc j ˛k+1 = ık+1 i,j j



k+1 − Xbc ), ∀j 2,j

− log(RAND(k+1 , 1)); i

∀j and ∀i

where s = {1, 2, . . ., N}, Ii = index of sorted fitness function, bc 1 and bc 2 are two mutually different indices of locally best seekers bc ∈ {bc1 , bc2 , bc3 } and RAND(k+1 , 1) generates i

, 1] random number uniformly in [k+1 i 8. Update seeker position k+1 k + ˛k+1 × dk+1 ; ∀j and ∀i = Xi,j Xi,j i,j i,j

9. Evaluate fitness Fik+1 = f (Xk+1 ), ∀i and find the index of the best i seeker b1 and the indices of locally best seekers bc1 , bc2 , and bc3 of the first, second and third 1/3rd of the newly generated solution Xk+1 10. Update Pbest, Gbest, Lbest and Xbest

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

43

MF, CR, F and HMCR have been varied randomly by ±20% (vis-avis their standard values), (b) the values of interval factor (for PSO and SOA), bandwidth and pitch adjustment rate (for HS) have been chosen randomly within their respective specified limits (instead of following the equations given in Sections 3.2, 3.4 and 3.5) and (c) following the recommendation of [22], the values of (min and max ) are kept unchanged. For further reference these modified values are termed as “perturbed values”. 4. Results and discussion

Fig. 1. Flowchart of the basic procedure.

∀i

If Fik+1 < Fik then Pbestk+1 = Xk+1 , ∀i else Pbestk+1 = Pbestki , i i i k+1 < Fbk then Gbestk+1 = Pbestk+1 and b = b1 else If Fb1 b1

Gbestk+1 = Gbestk Lbestk+1 = Pbestk+1 b1

k+1 Xbestk+1 bc = Pbestbc 11. If k < Maxite then k = k + 1 and goto step 6 else goto step 12 12. Print optimum solution as Gbestk

In this work, for evaluating the reproducibility of results obtained by different methods, each method has been executed 100 times. Subsequently, with these results obtained from 100 runs, various statistical parameters such as mean, standard deviation, etc., have been calculated for the purpose of overall comparison. Fig. 1 shows the comprehensive flowchart for the overall procedure adopted in this work. The following parameters have been considered for executing various methods. The population size of all the methods has been taken as 60 except DE and HS for which it is taken as 30 and 10, respectively, which gives better result. The crossover factor (CF) and the mutation factor (MF) for GA have been considered as 0.8 and 0.01, respectively [13]. For PSO, the two acceleration coefficients (c1 and c2 ) and the minimum and maximum value of inertial factor (wmin and wmax ) have been taken as 2.025, 2.025 and 0.4, 0.9, respectively [19]. For DE, the crossover rate (CR) and the mutation factor (F) have been considered as 0.4 and 0.5, respectively [20]. For HS, the harmony memory consideration rate (HMCR) and minimum and maximum values of bandwidth (BWmin and BWmax ) and pitch adjustment rate (PARmin and PARmax ) are considered as 0.9 and 0.0001, 1.0 and 0.4, 0.7, respectively [21]. Finally, for SOA, the minimum and maximum values of inertia factor (wmin and wmax ) and membership degree (min and max ) have been considered as 0.1, 0.9 and 0.0111, 0.95, respectively [22]. As can be seen from the above discussion, all these values of the parameters have been reported in the literature. For further reference, these values are termed as “standard values”. Initially, the performances of all the algorithms have been investigated with these standard values. However, for further investigations of the performances of the algorithms, these parameters have been varied from their “standard values” as follows: (a) the values of CF,

The algorithms described in Section 3 have been applied to solve the protection coordination problem of three power distribution systems of different sizes. In all these systems, numerical relays with standard IDMT characteristics have been considered. However, in the third test system, two cases have been considered. In the first case IDMT characteristics of all the relays have been assumed, whereas the second case considers other characteristic curves of the relays in addition to IDMT characteristics. The minimum and the maximum limits on TMS have been considered as 0.1 and 1.1, respectively, whereas the minimum and the maximum limits on PS are taken as: PSmin = max{0.5, min{1.25 × ILmax , 1/3 × IFmin }; and PSmax = min{2.5, 2/3 × IFmin }, where ILmax and IFmin are the maximum load current and minimum fault current, respectively [13,16]. Further, the minimum CTI considered for all the three system is 0.2 s [13]. Also, the minimum and the maximum limits on top has been considered as 0.1 s and 4.0 s, respectively. All these meta-heuristic optimization techniques have been implemented in MATLAB environment [23]. To investigate the performances of the algorithms comprehensively, all the algorithms have been executed using two sets of parameter values: “standard values” and “perturbed values” (as described in the previous section). In the following subsections, the results obtained on these three systems are presented in detail. 4.1. System I: 9-bus test system Fig. 2 shows the single line diagram of the 9-bus single source ring main distribution system. This system is fed at bus 1. Further, the detailed information about the system is available in [13]. In this system, there are 12 lines (L1,L2,. . ., L12) and 24 relays (R1,R2,. . .,R24) having 32 combinations of primary-backup relationship among them. Table 1 shows these 32 combinations of primary-backup relationship among all the 24 relays. The fault current passing through the primary circuit of the primary and the backup relays for various near-end three phase faults and the maximum as well as minimum fault current through the relays are given in [13] and hence are not repeated here. The current transformer (CT) ratio required for each relay is also given in [13]. The optimum settings, i.e., TMS and PS of the relays obtained by different methods using standard parameter values are given in Table 2, whereas the coordination time interval (CTI) between the time of operations of the primary relay (top ) and that of the corresponding backup relays (tob ) are plotted in Fig. 3. It is to be noted that the best results (for the minimum value of objective function) obtained after 100 runs of each method corresponding to the standard parameter values are given in Table 2 and Fig. 3. The last row in Table 2 gives the values of the sum of the operating times of all the relays, i.e., objective function value (OFV) corresponding to the set of TMSs and PSs obtained by each method. Table 3 shows the summary of the results obtained after 100 runs by the different methods with standard parameter values, whereas Table 4 shows the corresponding values obtained with the perturbed parameter values. These tables summarize the best, mean and the worst values of the objective function along with its standard deviation and the average elapsed time (s) for executing each algorithm.

44

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

Fig. 2. Single-line diagram of 9-bus system.

Table 1 Primary/backup relay pairs of the 9-bus system. Fault zone L1 L2 L3 L4 L5 L6

Primary relay

Backup relay

1 2 3 4 5 6 7 8 9 10 11 12

15, 17 4 1 6 3 8, 23 5, 23 10 7 12 9 14, 21

Fault zone L7 L8 L9 L10 L11 L12

Primary relay

Backup relay

13 14 15 16 17 18 19 20 21 22 23 24

11, 21 16, 19 13, 19 2, 17 – 2, 15 – 13, 16 – 11, 14 – 5, 8

Table 2 Optimum settings of relays obtained for the 9-bus system obtained with standard parameter values. Relays

GA TMS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 OFV

0.4395 0.3086 0.317 0.2194 0.4358 0.2289 0.4339 0.2093 0.3897 0.374 0.3491 0.2015 0.2756 0.1965 0.3812 0.3034 0.298 0.1 0.1686 0.1 0.2587 0.1 0.2607 0.1 14.5426

PSO

DE

HS

SOA

PS

TMS

PS

TMS

PS

TMS

PS

TMS

0.5 0.5001 1.3484 1.5477 0.5111 1.6634 0.5 1.7504 0.5 0.5 0.5 1.909 1.1158 1.607 0.6744 0.5002 0.8949 0.5193 1.6592 0.5 0.7805 0.5 0.9776 0.5001

0.2806 0.1636 0.3878 0.2015 0.2126 0.4266 0.2045 0.193 0.3401 0.4048 0.2246 0.4084 0.3046 0.1661 0.1832 0.3109 0.1555 0.1 0.2265 0.1001 0.1286 0.1001 0.2967 0.1 13.9472

1.2478 1.8196 0.5 2.4998 1.8823 0.5 2.5 2.5 0.5 0.5 1.0882 0.5 0.5537 2.4561 2.0707 0.5 1.9859 0.5041 0.75 0.5031 2.5 0.504 0.7501 0.5041

0.1241 0.1 0.137 0.1089 0.1237 0.1277 0.1277 0.1237 0.1089 0.137 0.1 0.1241 0.1 0.109 0.109 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1

2.5 2.0899 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.0899 2.5 2.2969 2.5 2.5 2.2969 2.1606 0.5 1.6462 0.5 2.1606 0.5 1.9435 0.5

0.1447 0.1 0.1684 0.1138 0.1309 0.1384 0.1388 0.13 0.1212 0.1598 0.1 0.1393 0.1021 0.1141 0.1165 0.1183 0.1 0.1002 0.1292 0.1001 0.1 0.1002 0.1334 0.1002

2.1741 2.294 1.8739 2.4472 2.4175 2.2897 2.3249 2.4176 2.2509 2.0335 2.3288 2.259 2.3465 2.4932 2.4666 1.936 2.3568 0.6198 1.2409 0.7395 2.4717 0.7203 1.4234 0.5572

0.2662 0.2076 0.2928 0.3192 0.2879 0.3677 0.3006 0.2905 0.2476 0.248 0.2578 0.3665 0.2581 0.3117 0.2921 0.3633 0.256 0.1038 0.2589 0.1002 0.2758 0.101 0.1757 0.1014

8.6822

9.2339

PS 1.2732 1.52 1.1975 0.6701 1.0785 0.6311 0.9637 1.1393 1.1994 1.7451 0.8454 0.6461 0.9784 0.886 0.8993 0.5004 0.9197 0.5003 0.7629 0.5041 0.8902 0.5008 1.5724 0.5017 14.2338

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

45

Table 3 Summary of results obtained after 100 independent runs for the 9-bus system with standard parameter values. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

14.5426 13.7150 8.6822 9.2339 14.2338

15.6592 20.6491 8.7162 9.6180 16.5710

16.8083 31.0790 9.4859 10.1411 19.9820

Standard deviation

Average elapsed time (s)

0.4110 3.6823 0.1233 0.2152 1.2133

360.70 22.22 7.29 122.15 30.20

Table 4 Summary of results obtained after 100 independent runs for the 9-bus system with perturbed parameter values. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

14.5888 13.6084 8.6822 9.2098 14.2501

15.888 22.7949 8.7060 9.2205 16.3973

18.0465 31.8319 9.5716 9.4667 21.6199

From Table 3 it is observed that among all these five methods, the DE gives the least value of the objective function (8.6822 s) and the standard deviation (0.1233). As the standard deviation is quite low for DE as compared to the other methods, it can be inferred that DE gives the most predictable results when it is executed repeatedly as compared to the other four methods. Further, the results obtained by DE are of high quality as the mean and the worst values of the objective function are close to the best value. Also, from Fig. 3, it is observed that among these five methods, the DE gives the lowest values of CTI. Again, from Table 4 it is observed that DE gives the best result (the least value of the objective function) among all these five methods. It is also observed from Tables 3 and 4 that with DE, the best value of the objective function remains the same with standard parameter values as well as with perturbed parameter values. Therefore, the best value (of the objective function) obtained by DE can be considered to be immune to variation in the parameter values. However, for the other four methods, the best values (and also the other values) of the objective functions change with variation in the parameter values.

Fig. 3. CTI obtained by various methods for 9-bus system obtained with standard parameter values.

Standard deviation

Average elapsed time (s)

0.6271 4.8349 0.1212 0.0368 1.1969

314.44 3.97 15.06 155.56 33.74

4.2. System II: 15-bus test system Fig. 4 shows single line diagram of the 15-bus network fed by an external grid (EG) at bus 8 and six distributed generators (DGs) connected at buses 1, 3, 4, 6, 13 and 15. The detailed information about the system is available in [16]. In this system, there are 21 lines (L1,L2,. . .,L21) and 42 relays (R1,R2,. . .,R42) having 82 combinations of primary-backup relationship among them. Table 5 shows these 82 combinations of primary-backup relationship among all the 42 relays. The fault current passing through the primary circuit of the primary and the backup relays for various near-end three phase faults are given in [16] and hence are not repeated here. The CT ratio required for each relay is also given in [16]. The optimum settings of the relays obtained by different methods obtained with standard parameter values are given in Table 6, whereas the coordination time interval (CTI) between the time of operations of the primary relay (top ) and that of the corresponding backup relays (tob ) are plotted in Fig. 5. Again, in this case, the best results obtained after 100 runs of each method corresponding to the standard parameter values are given in Table 6 and Fig. 5. Table 7 shows the summary of the results obtained after 100 runs by the different methods with standard parameter values, whereas Table 8 shows the corresponding values obtained with the perturbed parameter values. As in the previous case, these summaries include the best, mean and the worst values of the objective function along with its standard deviation and the average elapsed time (s) for executing each algorithm. From Table 7 it is again observed that for this system also the DE based technique is the best among the five meta-heuristic methods in terms of the values of the objective function (11.7591 s) and standard deviation (0.0004). Therefore, for this system also, out of the five methods, DE gives the best and most predictable results. Moreover, it is to be noted that the best objective function value obtained for this system so far reported in the literature is 12.2325 s as reported in [16], which is higher than the value of the objective function obtained by DE. Further, from Fig. 5, it is observed that among these five methods, DE gives the lowest values of CTI. Again, from Table 8 it is observed that DE gives the best result (the least value of the objective function) among all these five methods. It is also observed from Tables 7 and 8 that with DE, the best value of the objective function remains the same with standard parameter values as well as with the perturbed parameter values, thereby again establishing the fact that the best value obtained by DE is immune to the variation of the parameter values. However, again for the other four methods, the best values (and also

46

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

Table 5 Primary/backup relay pairs of the 15-bus system. Fault zone L1 L2 L3 L4 L5 L6 L7 L8 L9 L10 L11

Primary relay

Backup relay

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

6 4, 16 1, 16 7, 12, 20 2 8, 10 5, 10 3, 12, 20 5, 8 14 3, 7, 20 13, 24 9 11, 24 1, 4 18, 26 15, 26 19, 22, 30 3, 7, 12 17, 22, 30 17, 19, 30 23, 34

Fault zone L12 L13 L14 L15 L16 L17 L18 L19 L20 L21 –

Primary relay

Backup relay

23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 – –

11, 13 21, 34 15, 18 28, 36 25, 36 29, 32 17, 19, 22 27, 32 27, 29 33, 42 21, 23 31, 42 25, 28 38 35 40 37 41 31, 33 39 – –

Table 6 Optimum settings of relays obtained for the 15-bus system obtained with standard parameter values. Relays

GA TMS

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 OFV

0.2217 0.2199 0.2134 0.2169 0.2956 0.1814 0.2777 0.1503 0.2361 0.2701 0.2358 0.1692 0.2566 0.1435 0.2047 0.2355 0.2459 0.1775 0.2276 0.2037 0.245 0.2556 0.2422 0.2221 0.19 0.1858 0.2698 0.2826 0.2476 0.1931 0.2118 0.2553 0.3207 0.2809 0.1449 0.193 0.2321 0.1644 0.2657 0.2318 0.3092 0.1571 18.9033

PSO

DE

HS

SOA

PS

TMS

PS

TMS

PS

TMS

PS

TMS

0.5001 0.5001 0.9256 0.5372 0.5 1.2591 0.5001 1.3329 0.7428 0.5001 0.5 0.9745 0.5109 1.1425 0.5004 0.5017 0.5194 0.6713 0.7222 0.6488 0.5002 0.5029 0.5042 0.5001 1.2263 1.0327 0.5009 0.5035 0.5285 0.9965 0.9393 0.5001 0.5001 0.6255 1.8638 1.0986 0.9507 1.9288 0.5221 0.9407 0.5001 1.3041

0.1 0.1001 0.2134 0.1541 0.1636 0.2949 0.3441 0.3677 0.3529 0.1581 0.1 0.4058 0.4223 0.1 0.2635 0.3067 0.3741 0.1284 0.2671 0.1445 0.1 0.391 0.1243 0.1409 0.1867 0.1378 0.3496 1.0999 0.26 0.2592 0.3694 0.4379 0.3022 0.4165 0.3862 0.1763 0.423 0.1634 0.2078 0.2483 0.2052 0.2322 26.8093

2.5 1.4853 2.436 1.3103 2.5 0.5 0.7587 1.2353 0.6833 2.0514 2.5 0.5 0.5 2.5 0.5001 0.5 0.5239 1.4146 1.2694 2.2975 2.3325 0.5216 2.5 2.2981 2.5 2.4997 1.2764 0.5 2.5 0.5 1.8484 1.4562 2.3931 0.5 0.5 2.3456 0.5 2.5 2.5 0.9811 1.4688 2.4999

0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1043 0.1 0.1 0.1 0.1 0.1 0.1066 0.1 0.1 0.1023 0.1042 0.1014 0.1031 0.1035 0.1

1.2464 0.9762 2.0644 1.2439 2.1758 2.1084 2.0689 1.5847 2.1181 1.7132 1.2818 1.4139 2.1779 1.1156 1.0275 1.4475 1.5796 1.0498 1.809 1.3105 1.1437 1.6929 1.1259 1.4548 2.0088 1.7277 2.0028 2.5 1.5738 1.7436 1.8766 1.6002 2.4993 2.5 2.0723 1.8503 2.5 2.5 2.5 2.5 2.5 1.578

0.1 0.1001 0.1001 0.1001 0.1114 0.1057 0.1112 0.1001 0.1001 0.1001 0.1 0.1125 0.1001 0.114 0.1001 0.1 0.1001 0.1002 0.1133 0.1001 0.1028 0.1001 0.1 0.1043 0.1086 0.1001 0.1 0.1069 0.1001 0.1 0.1001 0.1001 0.1397 0.1222 0.1009 0.1004 0.108 0.1198 0.115 0.1361 0.1164 0.1

1.5933 1.2598 2.299 1.4479 1.9875 2.09 1.9172 1.7086 2.3959 2.3148 1.5747 1.259 2.3661 1.0666 1.2831 1.7563 1.7009 1.3617 1.5598 1.5965 1.2364 2.2202 1.6651 1.4555 1.9376 2.2356 2.3822 2.4603 1.8345 1.9207 2.4235 1.8187 1.635 2.1549 2.1173 1.974 2.428 2.0702 2.1461 1.7633 2.3633 1.7788

0.1555 0.2241 0.289 0.1642 0.3541 0.2449 0.2026 0.2693 0.3097 0.1793 0.1772 0.268 0.241 0.1642 0.1974 0.2421 0.2359 0.1698 0.2842 0.2564 0.2022 0.214 0.1224 0.2704 0.1819 0.2921 0.2058 0.3516 0.2023 0.259 0.2158 0.1352 0.3247 0.3283 0.3162 0.1966 0.3008 0.2967 0.2776 0.2201 0.3534 0.1491

11.7591

12.6225

PS 1.2031 0.6255 0.5859 1.0809 0.5 0.6088 1.3521 0.5414 0.5182 1.1962 0.6782 0.5052 0.9259 0.7968 0.583 0.6353 0.6312 0.7761 0.6357 0.5014 0.8159 1.1688 1.5117 0.5298 1.5846 0.5003 1.197 0.5066 1.0813 0.5573 1.1079 1.9317 0.6672 0.5265 0.5019 1.1711 0.5912 0.5405 0.5006 1.2674 0.5323 1.5585 20.4068

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

47

Fig. 4. Single-line diagram of 15-bus system.

the other values) of the objective functions change with variation in the parameter values. 4.3. System III: IEEE 30-bus test system Fig. 6 shows the 33 kV section of the IEEE 30-bus system. The system is fed through three 50 MVA, 132/33 kV transformers connected at buses 1, 6, and 13. In addition to the above three supply points, two DGs connected at bus no. 10 and 15 are also supplying the system. The detailed information about the system is available in [24,25]. The system has 20 lines (L1, L2, . . ., L20) and is protected

with 39 directional OCRs (R1,R2,. . .,R39) having 64 primary-backup combinations among them. Table 9 shows all the possible 64 combinations of the primary-backup relationships among the 39 relays. The fault current passing through the primary circuit of the primary and the backup relays for various near-end three phase faults are given in [24] and hence are not repeated here. The CT ratio for each relay is considered to be 500:1. It is to be noted that for this system some of the primary-backup relationships have been ignored while solving the coordination problem. These primary-backup relationships are 1-4, 17-4, 194, 28-34, 30-33, 31-34, 32-33 and 37-33. The reason behind this

Table 7 Summary of results obtained after 100 independent runs for the 15-bus system with standard parameter values. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

18.9033 26.8093 11.7591 12.6225 20.4068

19.9538 35.0252 11.7594 12.9014 24.9451

21.2223 44.8147 11.7615 13.2637 32.5519

Standard deviation

Average elapsed time (s)

0.6395 5.3496 0.0004 0.1983 1.8947

694.35 20.95 167.46 291.76 72.59

Table 8 Summary of results obtained after 100 independent runs for the 15-bus system with perturbed parameter values. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

24.6598 28.1412 11.7591 12.4895 20.1276

31.9816 46.9385 11.7610 12.5288 23.6208

41.3916 66.9602 11.8045 13.0173 31.4109

Standard deviation

Average elapsed time (s)

4.5996 9.1491 0.0073 0.0901 2.3619

421.42 9.71 98.39 311.20 80.98

48

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

Table 9 Primary/backup relay pairs of the 30-bus system. Fault zone L1 L2 L3 L4 L5 L6 L7 L8 L9 L10

Primary relay

Backup relay

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

4, 18, 20 6 2, 18, 20 5, 8 1 3, 8 3, 5 10, 36 7, 36 12 9 – 11 15 11 13 2, 4, 20 24 2, 4, 18 22

lies in the fact that for these combination the fault current passing through the corresponding backup relays are small (less than two times of the maximum load current of the relay) causing larger operating time of the backup relays and therefore, the minimum CTI requirements are always maintained for these combinations. For this test system two cases have been considered. In the first case, all the relays are assumed to have standard IDMT characteristics, whereas in the second case different characteristics curves of the relays have been considered. 4.3.1. Case A: relays with standard IDMT characteristics For this case, the optimum settings of the relays obtained (corresponding to the best run out of 100 independent runs) by the various methods corresponding to the standard parameter values are given in Table 10, whereas the coordination time interval (CTI) between the time of operations of the primary relay (top ) and that of the corresponding backup relays (tob ) are plotted in Fig. 7. Table 11 shows the summary of the results obtained after 100 runs by the different methods with standard parameter values, whereas

Fig. 5. CTI obtained by various methods for 15-bus system obtained with standard parameter values.

Fault zone L11 L12 L13 L14 L15 L16 L17 L18 L19 L20

Primary relay

Backup relay

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 –

19 26 17 28 21 30 23 32, 34 25 31, 33, 38 27, 34 29, 33, 38 27, 32 29, 31, 38 7, 10 37 29, 31, 33 35 9, 12 –

Table 12 shows the corresponding values obtained with the perturbed parameter values. As in the previous cases, these tables show the best, mean and the worst values of the objective function along with its standard deviation and the average elapsed time (s) for executing each algorithm. From Table 11, again it is observed that among all these five methods, the DE gives the least value of the objective function (17.8122 s). Further, the standard deviation (0.7226) is also quite low, establishing the high predictability of the results obtained in different runs of DE. Further, the results obtained by the DE are of high quality as the mean and the worst values of the objective function are close to that of the best value. From Fig. 7, it is again observed that among these five methods, DE gives the lowest values of CTI. Again, from Table 12 it is observed that DE gives the best result (the least value of the objective function) among all these five methods. It is also observed from Tables 11 and 12 that with DE, the best value of the objective function remain the same with standard parameter values as well as with the perturbed parameter values, thereby again establishing the fact that the best value obtained by DE is immune to the variation of the parameter values. However, again for the other four methods, the best values (and also the other values) of the objective functions change with variation in the parameter values. 4.3.2. Case B: relays with different characteristics In this case relays with different characteristic curves have been considered. For this case, the relays which are connected to the external supply buses are considered to have very inverse (VI) characteristics, the relays which are connected to the DG buses are considered to have extremely inverse (EI) characteristics and the rest are considered to be of standard IDMT type relays. Table 13 shows the details of various relays along with the coefficients of their characteristic curves [18]. For this case, the optimum settings of the relays obtained (corresponding to the best run out of 100 independent runs) by the various methods corresponding to the standard parameter values are given in Table 14, whereas the coordination time interval (CTI) between the time of operations of the primary relay (top ) and that of the corresponding backup relays (tob ) are plotted in Fig. 8. Table 15 shows the summary of the results obtained after 100 runs by the different methods with standard parameter values, whereas Table 16 shows the corresponding values obtained with the perturbed parameter values. As in the previous cases, these tables

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

49

Table 10 Optimum settings of relays obtained for the 30-bus system obtained with IDMT characteristics and standard parameter values. Relays

GA

PSO

TMS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 OFV

0.348 0.1001 0.3798 0.3554 0.3299 0.4352 0.1809 0.4889 0.1936 0.5572 0.1215 0.1 0.1001 0.1001 0.3973 0.4765 0.3487 0.4408 0.476 0.5129 0.2186 0.3895 0.3605 0.6013 0.2874 0.4992 0.4259 0.3122 0.3118 0.2663 0.1001 0.1 0.4798 0.4032 0.2961 0.3305 0.2389 0.1 0.5123 28.0195

DE

HS

SOA

PS

TMS

PS

TMS

PS

TMS

PS

TMS

0.5416 0.5001 1.1261 0.7318 1.1605 0.5002 0.941 0.5 0.5473 0.5001 0.5388 0.5581 0.5801 0.5025 0.6834 0.7692 1.2317 0.6306 0.7478 0.7213 1.9997 1.1956 1.1951 0.5393 1.5387 0.5 0.5 1.8681 1.494 1.7012 0.5001 0.5024 0.5001 0.6577 1.2358 1.6162 1.1937 0.5062 0.5001

0.4339 0.1 0.4759 0.2506 0.6753 0.771 0.2404 0.3494 0.1219 0.4973 0.1 0.1001 0.1001 0.1003 0.5204 0.9705 0.696 0.3678 0.4302 0.5525 0.2516 0.7658 0.6002 0.6159 0.5342 0.6434 0.2614 0.315 0.2962 0.2875 0.1004 0.1 0.7428 0.4862 0.5306 0.6806 0.2153 0.1001 0.8827 39.1836

0.7311 0.7793 1.5442 2.2529 0.5 0.5 0.6273 2.4872 2.4353 1.6339 2.4992 1.7099 0.8006 0.5 0.5002 0.5 0.5 2.451 2.291 1.4055 2.1498 0.8065 1.0292 1.0636 0.5 0.8444 2.5 2.5 2.5 2.4979 0.5 0.5001 0.5 1.2949 0.7374 0.5001 2.0325 0.5 0.717

0.1 0.1 0.1709 0.137 0.1372 0.1455 0.1 0.1514 0.1 0.1825 0.1 0.1 0.1 0.1 0.1654 0.181 0.1571 0.1387 0.1901 0.1924 0.135 0.1729 0.1602 0.209 0.151 0.1425 0.1267 0.173 0.1578 0.1239 0.1 0.1 0.1553 0.1389 0.1264 0.18 0.1 0.1 0.1657

1.9683 0.5 2.5 2.5 2.5 2.5 2.163 2.5 2.154 2.5 0.7809 0.5 0.5544 0.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 0.5 0.5 2.5 2.5 2.5 2.5 2.2539 0.5 2.5

0.1001 0.1001 0.2008 0.1534 0.1568 0.1618 0.1093 0.1678 0.1001 0.2191 0.104 0.1001 0.1052 0.1001 0.1832 0.2513 0.1935 0.168 0.2282 0.226 0.1631 0.2064 0.1889 0.2606 0.1656 0.173 0.1574 0.1903 0.1867 0.1395 0.1 0.1 0.1783 0.1659 0.1417 0.1983 0.1 0.1001 0.1898

2.2823 0.5006 2.0845 2.1177 2.2603 2.1949 1.9959 2.3441 2.1813 1.9803 0.73 0.5614 0.5235 0.5003 2.361 1.5888 2.0471 2.089 2.0034 2.128 2.0815 2.189 2.1125 1.89 2.414 2.1852 1.9668 2.4121 2.1393 2.4161 0.6023 0.7242 2.2162 2.0189 2.3048 2.3581 2.4106 0.6104 2.2205

0.2456 0.1006 0.4053 0.3748 0.4308 0.4538 0.1838 0.4547 0.1816 0.6056 0.1156 0.1 0.1289 0.1008 0.4841 0.7017 0.4025 0.3972 0.5379 0.4411 0.5264 0.4336 0.507 0.4521 0.5957 0.2889 0.4197 0.5891 0.4805 0.3784 0.1161 0.1075 0.3787 0.4874 0.5641 0.5417 0.5053 0.1074 0.4268

17.8122

19.2133

PS 1.419 0.51 1.2032 1.6489 0.7438 0.964 0.9266 0.894 0.6927 0.5109 0.6022 0.795 0.5031 0.5007 1.0237 0.5098 1.4966 1.1425 1.2087 1.4649 0.6436 1.4988 1.2602 1.4362 0.5518 2.0495 1.0216 0.5366 0.8305 0.8871 0.5024 0.5009 1.3351 0.6885 0.5061 0.7404 0.5226 0.5069 0.9769 33.7734

Table 11 Summary of results obtained after 100 independent runs for the 30-bus system with standard parameter values and IDMT characteristics. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

28.0195 39.1836 17.8122 19.2133 33.7734

29.2239 46.0857 18.4427 19.7119 35.8641

30.3700 58.5692 21.4340 20.2305 40.7814

Standard deviation

Average elapsed time (s)

0.6717 6.1688 0.7226 0.3626 2.0771

655.66 13.37 123.22 341.41 60.37

Table 12 Summary of results obtained after 100 independent runs for the 30-bus system with perturbed parameter values and IDMT characteristics. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

33.2164 34.7322 17.8122 19.0463 31.3752

38.3621 52.0418 18.0856 19.5602 35.4262

45.4163 76.1520 21.4402 20.5018 60.1147

Standard deviation

Average elapsed time (s)

2.8409 10.1870 0.6362 0.2779 3.8634

396.45 7.92 73.76 313.10 68.85

Table 13 Types of relays and its characteristics coefficients for the 30-bus system. Relays types

EI VI IDMT

Various relays

22, 25, 36, 37 1, 3, 12, 13, 15, 17, 19, 28, 31, 33 2-11, 14, 16, 18, 20, 21, 23, 24, 26–30, 32, 34, 35, 38, 39

Characteristic coefficients 



80 13.5 0.14

2 1 0.02

L 0 0 0

50

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

Fig. 6. Single-line diagram of 30-bus system.

show the best, mean and the worst values of the objective function along with its standard deviation and the average elapsed time (s) for executing each algorithm. From Table 15, again it is observed that among all these five methods, the DE gives the least value of the objective function (11.2804 s). Further, the standard deviation (0.0088) is also quite low, establishing the high predictability of the results obtained in different runs of DE. Further, the results obtained by the DE are of high quality as the mean and the worst values of the objective function is close to that of the best value. From Fig. 8, it is again

observed that among these five methods, DE gives the lowest values of CTI. In this case of IEEE 30-bus system, it is observed that the sum of the operating times of all the relays is very low (11.2804,s) as compared to case A (17.8122 s) in which all the relays are considered to be of IDMT characteristics. This is because of the fact that VI and EI types of relays operate much faster than the corresponding IDMT types of relays. Further, it is observed that the average execution time of various algorithms are higher than that in the previous

Fig. 7. CTI obtained by various methods for 30-bus system considering IDMT characteristic curves of relays and standard parameter values.

Fig. 8. CTI obtained by various methods for 30-bus system considering different characteristics curves of the relays and standard parameter values.

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

51

Table 14 Optimum settings of relays obtained for the 30-bus system considering different characteristics curves and standard parameter values. Relays

GA

PSO

TMS 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 OFV

0.647 0.1824 0.3398 0.1006 0.4587 0.2236 0.3806 0.2104 0.1477 0.3811 0.1355 0.5993 0.2001 0.1001 0.5759 0.1 0.4405 0.2514 0.1929 0.2166 0.2637 0.4443 0.2203 0.2663 0.3915 0.2757 0.1759 0.7854 0.2145 0.2286 0.6503 0.171 0.6897 0.1001 0.2758 0.54 0.5759 0.1966 0.1001 14.5616

DE

HS

SOA

PS

TMS

PS

TMS

PS

TMS

PS

TMS

1.4683 0.5008 1.2115 0.5016 0.5064 0.9054 0.5132 1.445 0.8937 0.5365 0.5276 1.9032 1.3997 0.5005 0.8173 0.5043 1.7051 0.5041 1.83 0.5048 0.5005 2.0787 0.5626 0.74 1.9055 0.7377 0.5397 0.9435 0.6507 1.2619 1.1948 0.5108 0.7988 0.5258 1.0407 1.4168 0.9354 2.1341 0.5038

0.2793 0.186 0.1001 0.1098 0.6303 0.2761 0.4464 0.1767 0.2112 0.25 0.1 1.0998 0.1001 0.1 1.0345 0.1 0.7189 0.282 0.9854 0.2281 0.1369 1.0999 0.1 0.1913 0.5152 0.2643 0.1011 0.8668 0.2398 0.4149 0.3959 0.1 0.6425 0.1 0.4552 0.1365 0.1 0.4588 0.1 17.127

2.4998 0.5 2.497 0.5001 0.5215 0.5 0.5028 2.0551 0.5 2.0413 1.8406 2.5 2.4887 0.5 0.6088 0.5 1.2503 0.5 0.8031 0.6114 2.5 1.3749 2.4877 1.7264 1.6743 0.9374 2.2132 0.8643 0.5015 0.5 1.692 2.4999 0.8308 0.5112 0.5 2.5 2.4948 0.5087 0.5001

0.1694 0.1 0.7994 0.1 0.1581 0.1171 0.1227 0.1306 0.1 0.1374 0.1 0.417 0.758 0.1 0.39 0.1 0.1655 0.1 0.1 0.1 0.1191 0.2546 0.1 0.1124 0.4745 0.1391 0.1 0.1 0.1 0.1277 0.1274 0.1 0.5672 0.1 0.1325 0.1176 0.1001 0.1447 0.1

2.3635 1.2324 0.7886 0.5 2.5 2.4999 2.5 2.5 1.7143 2.5 1.1477 1.9251 0.7184 0.5 0.9903 0.5001 2.4998 2.4118 2.3236 1.9251 2.4994 2.4995 2.1324 2.4999 1.5596 2.5 1.5244 2.1752 2.0647 2.5 2.3779 1.1364 0.8782 0.5 2.5 2.4999 1.9243 2.5 0.5

0.18 0.1 0.972 0.1001 0.1707 0.1605 0.1461 0.1396 0.1001 0.1661 0.1368 1.0044 0.2195 0.1001 0.3849 0.1001 0.256 0.1081 0.5019 0.1025 0.1479 0.7544 0.1003 0.1166 0.3654 0.1654 0.1097 0.6654 0.1 0.136 0.1568 0.1 0.1871 0.1 0.1451 0.9769 0.5773 0.1556 0.1

2.348 1.4679 0.7154 0.6203 2.4009 1.45 2.0543 2.4224 1.9709 1.9697 0.5098 1.3181 1.3337 1.0169 1.0011 0.5831 2.0345 2.2048 1.104 1.8833 1.7292 1.4803 2.1811 2.4503 1.8084 1.8923 1.3109 0.9027 2.3372 2.3901 2.1845 1.2598 1.5258 0.7002 2.2941 0.9235 0.862 2.3332 0.7078

0.5029 0.2059 0.4044 0.1041 0.3336 0.2345 0.3544 0.2742 0.1453 0.3785 0.1755 0.53 0.4676 0.1047 0.2319 0.1299 0.4782 0.2291 0.4061 0.2068 0.2771 0.6251 0.1984 0.2584 0.6648 0.2522 0.1117 0.4642 0.2658 0.2304 0.5575 0.1818 0.2935 0.1001 0.2338 0.6707 0.3389 0.29 0.109

11.2444

case. The reason behind this lies in the increased complexity of the modelling because of different relay curves. Again, from Table 16 it is observed that DE gives the best result (the least value of the objective function) among all these five methods. Also, from Tables 15 and 16 it is observed that with DE, the best

11.883

PS 1.5987 0.5063 1.1081 0.5176 1.03 0.9836 0.5 0.922 1.2266 0.5018 0.5 1.9725 0.9143 0.5001 1.2838 0.5021 1.6091 0.6349 1.3043 0.5503 0.5264 1.7435 0.6357 0.7423 1.5422 0.9963 1.5258 1.1711 0.5 1.2368 1.2866 0.5 1.2248 0.5074 1.4744 1.21 1.1365 0.984 0.5 14.7367

value of the objective function obtained with the standard parameter values and the perturbed parameter values is very close to each other. Thus, the best value obtained by DE is virtually immune to the variation of the parameter values for all practical purposes. However, again for the other four methods, the best values (and also

Table 15 Summary of results obtained after 100 independent runs for the 30-bus system with different characteristics curves and standard parameter values. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

14.5616 17.1270 11.2444 11.8830 14.7367

15.4757 23.1064 11.2778 12.1572 16.1104

16.5568 33.2713 11.3076 12.5398 18.6680

Standard deviation

Average elapsed time (s)

0.4148 3.6491 0.0134 0.1467 0.7146

394.74 101.77 716.32 624.96 124.54

Table 16 Summary of results obtained after 100 independent runs for the 30-bus system with different characteristics curves and perturbed parameter values. Methods

GA PSO DE HS SOA

Objective function value Best

Mean

Worst

14.5503 14.4910 11.2499 11.8067 14.1985

16.8485 23.8191 11.3053 12.1687 15.6728

21.9571 36.8190 11.6960 12.6070 18.3079

Standard deviation

Average elapsed time (s)

1.4913 4.8560 0.0709 0.1620 0.8336

261.7567 125.1248 751.7813 650.1724 181.7318

52

M.N. Alam et al. / Electric Power Systems Research 128 (2015) 39–52

the other values) of the objective functions change with variation in the parameter values. From the above studies carried out on three test systems it is observed that DE is the best technique among the meta-heuristic optimization techniques considered in the literature for solving protection coordination problems of directional overcurrent relays. 5. Conclusion In this work, the performances of five different metaheuristic algorithms, namely, GA, PSO, DE, HS and SOA have been investigated for solving the protection coordination problem of DOCRs. Each algorithm has been executed 100 times with same initial condition as well as with standard and perturbed parameter values in three test systems of various sizes. Based on the studies carried out in this work, DE algorithm has been found to be the best one among the five algorithms studied in this work on following counts: (a) The best value obtained by DE is always the lowest among the best values obtained by the five algorithms considered in the work. (b) The best value obtained by DE is virtually immune to the variations of the parameters of DE. (c) Even with different characteristic curves of the relays, DE gives the lowest value of the objective function. (d) Because of low value of standard deviation, the values obtained by DE are quite predictable as compared to the other four methods. As a result, DE can be considered to be the most suitable algorithm for the coordination of DOCRs among the five algorithms studied in this work. Acknowledgements The authors would like to thank the reviewers whose thought provoking comments helped to enhance the quality of the paper significantly. References [1] Y. Lu, J.-L. Chung, Detecting and solving the coordination curve intersection problem of overcurrent relays in subtransmission systems with a new method, Electr. Power Syst. Res. 95 (2013) 19–27. [2] M. Ezzeddine, R. Kaczmarek, M. Iftikhar, Coordination of directional overcurrent relays using a novel method to select their settings, IET Gener. Transm. Distrib. 5 (7) (2011) 743–750.

[3] H. Sharifian, H. Askarian Abyaneh, S. Salman, R. Mohammadi, F. Razavi, Determination of the minimum break point set using expert system and genetic algorithm, IEEE Trans. Power Deliv. 25 (3) (2010) 1284–1295. [4] M. Damborg, R. Ramaswami, S. Venkata, J. Postforoosh, Computer aided transmission protection system design. Part I: Algorithms, IEEE Trans. Power Appar. Syst. PAS-103 (1) (1984) 51–59. [5] B. Chattopadhyay, M. Sachdev, T. Sidhu, An on-line relay coordination algorithm for adaptive protection using linear programming technique, IEEE Trans. Power Deliv. 11 (1) (1996) 165–173. [6] A. Urdaneta, H. Restrepo, S. Marquez, J. Sanchez, Coordination of directional overcurrent relay timing using linear programming, IEEE Trans. Power Deliv. 11 (1) (1996) 122–129. [7] A. Urdaneta, R. Nadira, L. Perez Jimenez, Optimal coordination of directional overcurrent relays in interconnected power systems, IEEE Trans. Power Deliv. 3 (3) (1988) 903–911. [8] D. Birla, R. Maheshwari, H. Gupta, A new nonlinear directional overcurrent relay coordination technique, and banes and boons of near-end faults based approach, IEEE Trans. Power Deliv. 21 (3) (2006) 1176–1182. [9] H. Zeineldin, E. El-Saadany, M. Salama, Optimal coordination of overcurrent relays using a modified particle swarm optimization, Electr. Power Syst. Res. 76 (11) (2006) 988–995. [10] M.M. Mansour, S. Mekhamer, N.-S. El-Kharbawe, A modified particle swarm optimizer for the coordination of directional overcurrent relays, IEEE Trans. Power Deliv. 22 (3) (2007) 1400–1410. [11] F. Razavi, H.A. Abyaneh, M. Al-Dabbagh, R. Mohammadi, H. Torkaman, A new comprehensive genetic algorithm method for optimal overcurrent relays coordination, Electr. Power Syst. Res. 78 (4) (2008) 713–720. [12] A. Noghabi, J. Sadeh, H. Mashhadi, Considering different network topologies in optimal overcurrent relay coordination using a hybrid GA, IEEE Trans. Power Deliv. 24 (4) (2009) 1857–1863. [13] P. Bedekar, S. Bhide, Optimum coordination of directional overcurrent relays using the hybrid GA-NLP approach, IEEE Trans. Power Deliv. 26 (1) (2011) 109–119. [14] J. Moirangthem, K. Krishnanand, S. Dash, R. Ramaswami, Adaptive differential evolution algorithm for solving non-linear coordination problem of directional overcurrent relays, IET Transm. Gener. Distrib. 7 (4) (2013) 329–336. [15] M. Barzegari, S. Bathaee, M. Alizadeh, Optimal coordination of directional overcurrent relays using harmony search algorithm, in: Environment and Electrical Engineering (EEEIC), 2010 9th International Conference on, 2010, pp. 321–324. [16] T. Amraee, Coordination of directional overcurrent relays using seeker algorithm, IEEE Trans. Power Deliv. 27 (3) (2012) 1415–1422. [17] A. Abdelaziz, H. Talaat, A. Nosseir, A.A. Hajjar, An adaptive protection scheme for optimal coordination of overcurrent relays, Electr. Power Syst. Res. 61 (1) (2002) 1–9. [18] H.E.J. Gers, M. Juan, Protection of Electricity Distribution Networks, third ed., Institution of Engineering and Technology, 2011. [19] R. Eberhart, Y. Shi, Comparing inertia weights and constriction factors in particle swarm optimization, in: Evolutionary Computation, 2000. Proceedings of the 2000 Congress on, vol. 1, 2000, pp. 84–88. [20] S. Das, P. Suganthan, Differential evolution: a survey of the state-of-the-art, IEEE Trans. Evolut. Comput. 15 (1) (2011) 4–31. [21] S. Das, A. Mukhopadhyay, A. Roy, A. Abraham, B. Panigrahi, Exploratory power of the harmony search algorithm: analysis and improvements for global numerical optimization, IEEE Trans. Syst. Man Cybern. B: Cybern. 41 (1) (2011) 89–106. [22] C. Dai, W. Chen, Y. Zhu, X. Zhang, Seeker optimization algorithm for optimal reactive power dispatch, IEEE Trans. Power Syst. 24 (3) (2009) 1218–1231. [23] Matlab, Mathworks Inc., Massachusetts, USA, Version r2012a. [24] R. Mohammadi, H. Abyaneh, H. Rudsari, S. Fathi, H. Rastegar, Overcurrent relays coordination considering the priority of constraints, IEEE Trans. Power Deliv. 26 (3) (2011) 1927–1938. [25] R. Christie, Power system test cases, Available [Online]: www.ee.washington. edu/research/pstca, August 1993.