- Email: [email protected]
Genetic programming and floating boom performance
Electric Power Systems Research 127 (2015) 118–125
Contents lists available at ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
Online coordination of directional overcurrent relays using binary integer programming Rafael Corrêa a,∗ , Ghendy Cardoso Jr. b , Olinto C.B. de Araújo b , Lenois Mariotto b a b
Federal Institute of Education, Science and Technology from Rio Grande do Sul, Farroupilha, RS, Brazil Federal University of Santa Maria, Santa Maria, RS, Brazil
a r t i c l e
i n f o
Article history: Received 7 October 2014 Received in revised form 21 May 2015 Accepted 23 May 2015 Keywords: Binary integer programming Directional overcurrent relay Online coordination SCADA system Smart grid Optimization
a b s t r a c t This work presents a binary programming model for online coordination of directional overcurrent relay problems in interconnected power distribution and subtransmission systems. The proposed model considers discrete time dial and pickup current settings as in microprocessor-based relays. The coordination problem is solved for every topological change in the network and whose optimized settings are remotely adjusted on each relay through a SCADA system, enhancing the speed and sensitivity of the protection system. A pre-processing step is performed to ensure that the operation time of each relay be below the melting curve of the protected cable. Computational results for the binary programming model solved by CPLEX optimization package are compared to heuristic-based techniques. The simulation results indicates that our approach is suitable for online applications, since an optimal solution can be found in reduced times. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Overcurrent relays are used as primary protection for distribution and subtransmission lines as well as backup protection for transmission lines, power transformers and generators. In fault situations, when the primary relay does not operate, a remote backup relay must operate and the associated circuit breaker eliminates the fault. In order to guarantee that the backup relay does not operate before the circuit breaker associated to the primary relay, an intentional time delay is considered to adjust the backup. This process is called coordination. This protection philosophy is used in power distribution and subtransmission networks. In interconnected power systems and considering distributed generation (DG), a known adversity in coordinating the non-directional overcurrent protection requires the application of directional overcurrent relays (DOCR). In some cases, selecting the appropriate settings for DOCR to ensure coordination, speed and sensitivity for the protection system may be a challenge. For instance, it is common that minimum fault currents at the end of the backup zone be lower than the maximum load current. In this situation, when a pickup current above
∗ Corresponding author. Tel.: +55 54 32602400. E-mail addresses: [email protected] (R. Corrêa), [email protected] (G. Cardoso Jr.), [email protected] (O.C.B.d. Araújo), [email protected] (L. Mariotto). http://dx.doi.org/10.1016/j.epsr.2015.05.017 0378-7796/© 2015 Elsevier B.V. All rights reserved.
the maximum load current is chosen, the relay will not provide backup protection for such fault currents. In an offline coordination process, considering the different n − 1 contingency topologies, the obtained settings may lead to high operation times for some topologies. To overcome this issue [1,2] propose to re-coordinate all DOCR for every topological change in the network. This concept is also used in this study, considering the possibilities and the growth of smart grids. The coordination of DOCR has been formulated as an optimization problem and solved with several techniques. To determine continuous time dial settings (TDS), linear programming solvers have been used in [3–6] considering fixed pickup current settings (PCS). In [7,8], mixed integer non-linear programming and mixed integer programming solvers, respectively, have been used to determine discrete PCS and continuous TDS. Continuous TDS and PCS have been determined using a non-linear programming solver in [9,10]. Evolutionary techniques have been used, such as the genetic algorithm (GA) [11–15], nondominated sorting genetic algorithmII (NSGA-II) [16], particle swarm optimization (PSO) [17,18], seeker optimization algorithm (SOA) [19], ant colony algorithm [1], differential evolution (DE) algorithm [1], opposition based chaotic differential evolution (OCDE) algorithm [20], adaptive differential evolution (ADE) algorithm [21] and teaching learning-based optimization algorithm [22]. Hybrid techniques have also been used in [23,24], combining metaheuristics with deterministic solvers. The main drawback from these heuristic-based techniques is the
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
presence of parameters that must be correctly chosen in order to provide the convergence of the algorithm. The optimal parameters may not be the same for different instances, thus this choice is performed in a trial-and-error procedure. Also, given the non deterministic nature of such algorithms, more than one simulation may be necessary to achieve the optimal. Nonconventional, or non-standardized, characteristics for inverse DOCR have been considered in [25,26]. The GA and a nonlinear programming problem solver, respectively, have been used to optimize TDS, PCS and the shape – characteristic – of the curve of each DOCR, taking advantage of numerical relays. In general, when the coordination of DOCR is addressed using the techniques previously listed, the TDS are determined as a continuous variable. However, if the step of the TDS of a microprocessor-based relay is not small enough, such as 0.01, miscoordinations may occur whether we round the continuous solution to the nearest available settings [12]. To overcome the described issues, this paper presents a binary integer programming (BIP) model considering discrete TDS and PCS, within the available range. The online coordination of DOCR for three test systems in presence of DG is performed using the CPLEX optimization package. In order to show the effectiveness of the proposed approach, the results are compared with recently used heuristic-based techniques. The novelties of this paper are described below: • The coordination of DOCR problem is modeled using a BIP model and solved by a mixed integer programming solver running on its default configuration, i.e., the user does not need to set parameters. Also, a single simulation is necessary to achieve the optimal; • The TDS and PCS are determined in the discrete form, as available in electromechanical and microprocessor-based units with noncontinuous settings; • A pre-processing step is performed to ensure that the operating times of the DOCR be below the melting curve of the protected cables. This step reduces the simulation times, making the proposed approach suitable for online applications. The other sections of this paper are organized as follows: in Section 2, the conventional formulation of the DOCR problem is described. The proposed online coordination procedure is detailed in Section 3. In Section 4, the proposed approach is applied over three test systems and the simulation results are discussed. Finally, the conclusions are provided in Section 5.
119
The constraints of the problem are defined as follows. 2.1. Limits of the settings The bounds of the settings can be expressed by TDSmin ≤ TDSi ≤ TDSmax , i i PCSmin i
≤ PCSi ≤
PCSmax , i
i = 1, . . ., m i = 1, . . ., m
(3) (4)
where TDSi min and TDSi max are the lower and upper limits of TDS of the ith relay, respectively; PCSi min and PCSi max are the lower and upper limits of PCS of the ith relay, respectively. Moreover, for phase protection, the PCS must be greater than maximum load current and lower than minimum fault current seen by each relay, including a safety margin, which depends on the relay technology and current transformer (CT) errors. 2.2. Limits of relay operation time Minimum and maximum operating times can be established for each relay by min max Tik ≤ Tik ≤ Tik , min
i = 1, . . ., m
(5)
max
and Tik are the minimum and maximum operwhere Tik ation times of the ith relay for a fault at k, respectively. 2.3. Coordination of primary and backup relays The backup relays must operate with an intentional time delay, called coordination time interval (CTI), in relation to the primary relay. Thereby, the backup operates only when the primary protection fails. This time interval includes factors such as the circuit breaker operation time, CT errors, relay overtravel time and a safety margin. Its value is usually selected between 0.2 s and 0.5 s. Tjk − Tik − CTI ≥ 0,
i = 1, . . ., m
(6)
where Tjk is the operation time of the jth backup relay for a fault at k inside the zone protected by the ith primary relay. This constraint must be satisfied for all primary/backup pairs. Usually, near-end and far-end fault currents are used, providing coordination for most of the fault situations. 2.4. Topological changes of the network
2. Conventional problem formulation
The coordination constraint (6) is considered not only to the main topology of the network, but also to other possible configurations.
The objective of the DOCR coordination problem is to minimize the operating times of the relays operating as primary protection. The objective function (OF) is then expressed by
3. Proposed approach
min OF =
PCSi ,TDSi
m i=1
Tik
(1)
k
where m is the number of relays in the network; PCSi and TDSi are the pickup current and time dial settings of the ith relay, respectively; and Tik is the operation time of the ith relay for a fault at k, defined as follows.
Tik = TDSi ×
A Iik /PCSi
P
−1
+B
(2)
where Iik is the fault current seen by the ith relay for a fault at k; A, B and P are standardized coefficients. In this work, the Type A curve from the IEC 255-3 has been used, where A = 0.14, B = 0 and P = 0.02.
This section presents the proposed online coordination procedure. The main steps are the supervisory control and data acquisition (SCADA) system, the pre-processing step and the process to create and solve the BIP model. 3.1. SCADA system The main tasks of the SCADA system are to run load flow and short-circuit analysis for each topological change of the power network, to optimize the TDSs and PCSs, and to re-set each relay remotely. Although the communication between a remote device and a supervisory system is already a reality in the state-of-the-art of digital relays, a fully implemented system to execute such resetting in interconnected power distribution and subtransmission networks is not available. Thus, this re-setting system is supposed to have been implemented.
120
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
According to each topological change of the network, the SCADA system processes the new topology and the following steps are performed: • The coordination pairs of relays are determined and stored as data; • Load flow analysis is executed for the n − 1 contingency topologies. The resulting maximum load current for each relay is stored as data; • Short-circuit analysis is executed considering the new topology. The fault currents related to each coordination pair of relays are stored as data.
Afterwards, the maximum load and fault currents are used to construct the BIP model. This model is solved and the optimal settings are adjusted on each relay via communication hardware. By the proposed approach, relays are coordinated for the current topology of the network. Whether any change occurs, the online coordination process starts, and coordination may be lost if a fault occurs before the new settings are determined and sent to each relay. However, the DOCR do not operate if a fault does not occur, since the pickup currents are higher than maximum load currents considering the n − 1 criterion. It is important to note that many types of relays cannot be remotely adjusted, like electromechanical units. If they require their settings to be changed, they must be locally adjusted. For example, in an expansion of the network, maximum load current may overcome the pickup current, thus the PCS must be increased. Therefore, two solutions are proposed:
(1) If these units require their settings to be changed, the BIP model is constructed considering fault currents for all possible topologies of the network. These units are then locally adjusted and digital relays are remotely adjusted. (2) If these units do not require their settings to be changed, the BIP model is constructed with fixed settings for these units, reducing the number of variables of the problem.
3.2. Binary integer programming A BIP model has been used to deal with TDS and PCS as discrete variables. The BIP is an integer programming problem with binary variables. A multivariable BIP problem can be written in the matrix form as follows. min FO = cT x s.t. Ax ≤ b
(7)
x ∈ Bn where cT is the transposed cost vector; x is the variables vector; A is a matrix with constants; b is a vector with constants; and Bn represents the space of vectors with n binary components. The branch and bound and the branch and cut are well known algorithms used for solving BIP problems. Commercial solvers for BIP problems incorporate both classical and recently developed techniques. For instance, the MIPOPT solver from CPLEX offers a dynamic search algorithm, which consists in LP relaxation, branching, cuts, and heuristics. Since it is not possible to detail all techniques implemented inside the CPLEX, the reader may find the description of the branch and bound and the branch and cut algorithms, and how they work, in [27].
3.3. Representation of the variables To represent TDS and PCS in the discrete form and to linearize the constraints and the objective function of the model presented in Section 2, the following changes have been made in the conventional problem formulation. A constant cirsk has been introduced, considering discrete TDSs and PCSs. Its value is calculated as follows.
cirsk = TDSir ×
A Iik /PCSis
P
−1
+B
(8)
where TDSir is the rth element of the set vTDSi . This set contains all nTDSi values of TDS for the ith relay; PCSis is the sth element of the set vPCSi . This set contains all nPCSi values of PCS for the ith relay; nTDSi is the number of available values of TDS for the ith relay; nPCSi is the number of available values of PCS for the ith relay; Iik is the fault current seen by the ith relay for a fault at k; A, B and P are standardized coefficients of the selected curve. The relay operating time is then written using (8) and binary variables xirs in the form
nTDSi nPCSi
Tik =
cirsk xirs
(9)
r=1 s=1
where xirs is a binary variable: equal to 1 if the rth TDS value and the sth PCS value are assigned to the ith relay; and 0 otherwise. Each variable xirs represents a combination of a TDS and a PCS for the ith relay. To ensure a single combination of TDS and PCS has been assigned to each relay, the following constraints are included in the model.
nTDSi nPCSi
xirs = 1,
i = 1, . . ., m
(10)
r=1 s=1
xirs ∈
0, 1 ,
i = 1, . . ., m;
r = 1, . . ., nTDSi ;
s = 1, . . ., nPCSi (11)
In (9), the relay operating time is linear with respect to xirs . The objective function (1) as well as the constraints (5) and (6) are then written using (8) and (9), and the resulting problem is linear. The sets vTDSi and vPCSi contain all possible values of TDS and PCS for the ith relay, respectively. For example, if relay number 3 has TDS available between 0.1 and 0.4 in steps of 0.1, the set is vTDS3 = {0.1; 0.2; 0.3; 0.4}. Using this representation, it is possible to model electromechanical and microprocessor-based relays. In electromechanical units, usually the PCS do not have a constant step, thus it can be written as vPCSi = {0.5; 0.6; 0.8; 1.0; 1.5; 2.0; 2.5}. Also, a step of 0.01 can be used to represent continuous settings of digital relays, without prejudice in the speed and sensitivity. As mentioned previously, each variable xirs represents a combination of a TDS and a PCS for the ith relay. For example, if the problem has 10 relays and each one presents 50 available TDSs and 40 available PCSs, then the problem has 10 × 50 × 40 = 20,000 variables. Similarly, if each relay has 50 values for TDS between 0.1 and 1.1 and 40 values for PCS between 0.5 and 4.0, the variable x5,1,1 represents a TDS 0.1 and a PCS 0.5 for relay number 5; the variable x2,1,40 represents a TDS 0.1 and PCS 4.0 for relay number 2; and so on. 3.4. Pre-processing step In real DOCR, hundreds of settings are available and the coordination problem has a large number of variables, as in the proposed
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
121
Fig. 1. Technical constraints of the pre-processing step.
BIP model each variable represents a combination of TDS and PCS. Although such a problem can be solved in the actual state-of-theart of the mixed integer programming solvers, in some cases the simulation times may be prohibitive for online purposes. In order to reduce the simulation time by reducing the number of variables, a pre-processing step is proposed to eliminate variables that do not match the following technical constraints: • The available PCS must be within the range defined in Section 2.1; • The relay operating curve must be below the melting curve of the protected cable; • The relay operating time cannot be greater than 5 s for faults in the primary zone; • The relay operating time cannot be greater than 10 s for faults in the remote backup zone. These technical constraints are shown in Fig. 1. The choice of maximum values for the operating times as primary and remote backup protection is based on experience. Before constructing the BIP model, all combinations of TDS and PCS of each relay are tested considering the technical constraints. The pairs that do not match these constraints are dismissed and do not compose the BIP model. For instance, in Fig. 1 it can be observed that if the TDS is increased, from a certain value the operating curve of the relay overcomes the melting curve of the cable. Therefore, any TDS above this limit does not compose the final model for that particular PCS value. It is important to note that without this pre-processing the variables that would be eliminated at this step are included in the BIP model, which could lead to higher simulation times due to the larger number of variables. Also, this step works as a constraint related to the melting curve of the cables, as the variables that could violate this constraint are not included in the final model. The flowchart of the pre-processing step is shown in Fig. 2.
Fig. 2. Flowchart of pre-processing step.
4.1. IEEE 6-bus test system The IEEE 6-bus test system is shown in Fig. 3. Three-phase fault currents of the main topology shown in Fig. 3 and CT ratios can be found in [20]. The available PCS are from 1.25 to 1.50, in steps of 0.01, and the available TDS are from 0.05 to 1.10, in steps of 0.01. Also, a CTI 0.2 s has been used. In [20], simulation results for DE algorithm and an OCDE algorithm, named OCDE2, have been presented. To make a clear comparison, the objective function of the BIP approach has been defined as in [20], considering near-end and far-end fault currents. The results are shown in Table 1. The simulation time using the BIP approach was 13.08 s. From Table 1, it can be observed that the OF using the BIP approach is lower than the DE algorithm, even with a worse TDS and PCS resolution. The OF using OCDE2 algorithm is lower than the proposed technique, but the continuous settings may not be available in some relays. On the contrary, the optimized settings
4. Results and discussion The BIP approach was applied in three test systems and the results have been compared to heuristic-based techniques. The coordination was carried out for only phase relays, since the ground relay coordination uses the same procedure with ground fault currents. In addition, the instantaneous units were not considered. Algorithms in C++ have been implemented to run the CPLEX 12.4 mixed integer programming solver. They have been executed on a Core i3, 2.13 GHz PC with 4 GB RAM.
Fig. 3. Single-line diagram of the IEEE 6-bus test system.
122
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
Table 1 Results for IEEE 6-bus test system. Relay
DE [20]
Table 2 Results of Case 1 for 8-bus test system. OCDE2 [20]
BIP
Relay
TDS
PCS
TDS
PCS
TDS
PCS
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.1173 0.2082 0.0997 0.1125 0.0500 0.0580 0.0500 0.0500 0.0500 0.0719 0.0649 0.0617 0.0500 0.0856
1.2505 1.2500 1.2512 1.2515 1.2500 1.2500 1.2500 1.2500 1.2502 1.2502 1.4998 1.2575 1.4805 1.2557
0.1014 0.1862 0.0946 0.1006 0.0500 0.0500 0.0500 0.0500 0.0500 0.0562 0.0649 0.0508 0.0500 0.0708
1.5000 1.4999 1.2772 1.4995 1.2501 1.3801 1.2500 1.2500 1.2500 1.4988 1.4999 1.5000 1.4608 1.4995
0.11 0.19 0.10 0.11 0.05 0.05 0.05 0.05 0.05 0.07 0.07 0.06 0.05 0.08
1.38 1.46 1.25 1.31 1.25 1.41 1.25 1.25 1.25 1.25 1.44 1.29 1.44 1.38
OF (s)
10.6272
10.3286
10.5380
Hybrid GA-LP [23]
SOA [19]
TDS
PCS
TDS
PCS
BIP TDS
PCS
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.3043 0.2917 0.2543 0.1851 0.1700 0.2711 0.5316 0.2387 0.1865 0.1895 0.2014 0.2890 0.2207 0.5278
1.0 2.5 2.5 2.5 1.5 2.5 0.5 2.5 2.0 2.5 2.5 2.5 1.5 0.5
0.113 0.260 0.225 0.160 0.100 0.173 0.243 0.170 0.147 0.176 0.187 0.266 0.114 0.246
2.0 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.0 2.5
0.10 0.28 0.24 0.19 0.10 0.18 0.26 0.17 0.15 0.18 0.19 0.27 0.10 0.25
2.5 2.5 2.5 2.0 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
OF (s)
10.9499
8.4270
8.6944
using the BIP approach are discrete and within the available values for each relay.
currents, with emphasis between relays number 11 and 12 where the left side of the coordination constraint (6) is equal to −0.111.
4.2. 8-Bus test system
4.2.2. Case 2 In Case 2, the line where relays number 1 and 8 are located was removed and the coordination was performed. These devices do not compose the BIP model. This simulation case is to show the effectiveness of the BIP approach dealing with line outages and different topologies of the network. The simulation time of the proposed BIP model was 38.75 s. Table 3 shows the simulation results. It can be seen that relays number 1, 8 and 13 have the highest available value for TDS. These relays do not operate for any fault on this topology, thus any TDS can be used for these units. Once the PCS are optimized considering the maximum load currents, these relays do not trip after a change in topology if a fault does not occur. After this change, new settings must be determined to ensure the coordination in fault situations. Additionally, relay number 6 has the lowest available TDS, as it is the downstream device for coordination purposes in this topology.
The 8-bus test system is shown in Fig. 4. The data of the network can be found in [17]. A CTI 0.3 s and cables 1/0 CAA have been considered. The available TDSs are from 0.10 and 1.10, in steps of 0.01, and the available PCSs are 0.5, 0.6, 0.8, 1.0, 1.5, 2.0 and 2.5. Three cases have been simulated, as described below. 4.2.1. Case 1 Hybrid GA-LP [23] and SOA [19] considered near-end threephase fault currents of the main topology shown in Fig. 4. The same fault currents have been used in the BIP approach. The results of Case 1 are shown in Table 2. The simulation time using the BIP approach was 34.11 s. Without the pre-processing procedure, the time was 55.30 s. It is important to note that a single execution of the algorithm is required to reach an optimal solution, which is not always true considering heuristic techniques. From Table 2, it can be noted that the OF using the BIP approach is lower than the hybrid GA-LP, even with a worse TDS resolution. The OF using SOA is lower than the proposed technique, but the TDS are continuous. To illustrate the miscoordinations that may occur if we round these settings to the nearest available value, we considered the SIPROTEC® 7SJ46 numerical overcurrent relay, which has a step of 0.05. Nine miscoordinations occur for near-end three-phase fault
4.2.3. Case 3 An outage of generator G1 was considered in Case 3. In this situation, all relays compose the BIP model and fault current levels are lower than in Case 1. The simulation time of the BIP approach was 22.23 s. Table 3 shows the simulation results. It can be seen that lower TDS values have been determined for most relays in comparison with Case 1. Table 3 Results of Case 2 and Case 3 for 8-bus test system. Relay
Fig. 4. Single-line diagram of the 8-bus test system.
Case 2
Case 3
TDS
PCS
TDS
PCS
1 2 3 4 5 6 7 8 9 10 11 12 13 14
1.10 0.32 0.33 0.23 0.16 0.10 0.35 1.10 0.24 0.25 0.25 0.33 1.10 0.36
2.0 2.5 2.0 2.5 2.5 0.5 2.5 2.0 2.5 2.5 2.5 2.5 2.0 2.5
0.11 0.15 0.13 0.13 0.12 0.16 0.22 0.10 0.11 0.14 0.12 0.21 0.10 0.10
2.0 2.5 2.5 2.5 1.5 2.5 2.5 2.0 2.5 2.5 2.5 2.0 0.5 1.5
OF (s)
9.0747
7.3972
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
123
Table 4 RTC and maximum load current for IEEE 30-bus test system.
Fig. 5. Single-line diagram of the IEEE 30-bus test system.
Once the inverse characteristic curves have been considered, lower TDS are expected to ensure the speed of the protection system for lower fault currents.
Relay
RTC
max Iload (A)
Relay
RTC
max Iload (A)
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
150 150 100 100 100 100 100 150 100 100 100 50 150 100 100 100 100 50 100 50 100 100 150 100 150 150 150 150 150 100 100 150 50 100 100 250 150 150 150
637 637 470 0 0 361 431 519 0 110 483 0 561 0 0 202 209 0 66 144 225 73 461 0 452 345 235 235 244 11 0 539 0 460 0 1086 288 451 300
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
50 50 100 50 100 100 50 100 50 50 100 50 50 50 50 100 50 50 100 50 100 100 100 100 50 100 50 50 100 100 100 50 50 50 50 100 100 100 50
0 191 110 183 278 258 228 304 134 0 307 81 153 220 113 220 162 158 194 98 0 0 108 347 173 269 220 0 445 244 247 0 201 47 0 0 350 119 0
4.3. IEEE 30-bus test system The IEEE 30-bus test system, shown in Fig. 5, has 78 DOCRs and 206 coordination pairs in the main topology. The RTC and the maximum load current for each relay are given in Table 4.
Fig. 6. Convergence of the DE algorithm for the IEEE 30-bus test system.
124
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125 Table 5 (Continued)
Table 5 Results for IEEE 30-bus test system. Relay
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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
Relay
BIP
BIP
DE [1]
TDS
PCS
TDS
PCS
0.01 0.01 0.02 0.06 0.06 0.04 0.07 0.06 0.10 0.04 0.01 0.17 0.03 0.09 0.10 0.07 0.12 0.22 0.02 0.09 0.10 0.08 0.05 0.01 0.10 0.07 0.07 0.07 0.14 0.05 0.08 0.08 0.04 0.03 0.11 0.05 0.07 0.08 0.04 0.12 0.08 0.05 0.08 0.11 0.10 0.12 0.07 0.12 0.13 0.08 0.10 0.10 0.15 0.07 0.07 0.09 0.08 0.05 0.07 0.04 0.04 0.08 0.04 0.06 0.08 0.09 0.11 0.01 0.02 0.02 0.01 0.01 0.03 0.01
1.60 1.60 2.70 0.80 0.80 1.80 1.10 0.80 0.50 2.70 2.50 0.60 1.90 0.50 1.50 1.20 1.20 0.50 1.90 1.50 0.50 0.50 0.80 0.50 2.40 3.00 2.00 2.00 0.60 2.60 0.50 1.10 1.50 1.90 0.60 1.60 3.00 3.00 2.90 0.70 2.30 2.90 2.10 2.80 1.00 3.00 3.00 1.60 1.00 3.00 2.90 2.80 2.60 2.30 3.00 2.30 2.70 2.90 2.80 0.60 0.60 3.00 2.80 2.60 0.60 2.90 1.00 1.80 1.90 2.30 0.50 2.90 1.50 0.50
0.0296 0.0346 0.0665 0.0793 0.0460 0.1169 0.0923 0.0805 0.0386 0.1117 0.0284 0.1115 0.0770 0.0953 0.1190 0.1713 0.1596 0.2000 0.1248 0.1550 0.0972 0.0579 0.0451 0.0100 0.2052 0.0857 0.1512 0.1469 0.1354 0.1383 0.0584 0.1147 0.0542 0.0630 0.1268 0.0684 0.0812 0.1765 0.1013 0.1953 0.1861 0.1519 0.0987 0.2485 0.1131 0.2386 0.1355 0.1521 0.1686 0.1066 0.1802 0.2351 0.2247 0.2365 0.1516 0.1793 0.2290 0.1669 0.1665 0.0656 0.0573 0.2794 0.1803 0.1457 0.1046 0.2021 0.1056 0.0653 0.0346 0.0560 0.0100 0.0474 0.0327 0.0100
1.1954 1.1085 1.7014 0.8681 1.5224 0.8042 0.9314 0.7609 1.4716 1.2382 1.6166 1.4398 0.8632 0.5685 1.3892 0.5702 0.9721 0.7917 0.5376 0.8066 0.6047 0.6794 0.8390 0.5000 1.0010 2.9939 1.3237 1.3724 1.1047 0.5266 0.7977 0.8001 1.3768 1.3302 0.7478 1.3484 2.9993 1.6597 1.8395 0.5001 1.1787 1.0852 1.7317 1.1345 1.1963 1.1189 1.6324 1.4450 0.9347 2.9986 1.5700 1.1836 1.7144 0.5279 1.7245 1.3946 0.7000 1.2399 0.8556 0.5001 0.5416 0.6463 0.7053 1.4774 0.9436 1.4395 1.4336 0.9122 1.4313 1.2673 0.5000 0.9583 1.3372 0.5000
DE [1]
TDS
PCS
TDS
PCS
75 76 77 78
0.08 0.07 0.05 0.11
0.50 2.70 2.50 0.50
0.0697 0.1789 0.1367 0.1051
0.6281 1.4193 0.5001 0.8290
OF (s) Simul. time (s)
89.2199 284.06
112.7958 1041.49
A CTI 0.2 s, cables 336.4 CAA and a large number of settings have been considered: the TDS is from 0.01 to 0.50, in steps of 0.01, and the PCS is from 0.5 to 3.0, in steps of 0.1. Also, the maximum operating time of the primary protection was reduced to 2 s. The BIP approach was applied in the main topology of the network. The following fault currents were considered to compose the coordination constraints and the objective function: three-phase near-end; and far-end with the remote circuit breaker opened. The far-end fault currents with the remote circuit breaker closed were also considered to compose the objective function. The results for this case are shown in Table 5. Using the preprocessing step, the simulation time was 284.06 s and the number of variables was 46,598. This step takes only 0.234 s. Without the pre-processing, the simulation time was 5593.07 s and the number of variables was 79,550. The reduction of 41.4% in the number of variables is due to the maximum operating time of the primary relay (2 s) considered in the pre-processing. The maximum operating time of the primary relay of 5 s was also used: the number of variables was 66,663 (reduction of 16.2%) and the pre-processing time was 0.343 s. For comparison purposes, the continuous DE algorithm used in [1] were executed 10 times with a population size of 200 individuals and 10,000 maximum iterations. The objective function used in [1] is different from (1), so we considered (12) to make a clear comparison with the BIP results. min OF =
PCSi ,TDSi
i
k
Tik + ˇ ×
j
|ejk |
(12)
k
where Tik is the operation time of the ith relay for a fault at k; ejk is the CTI error between the jth coordination pair for a fault at k. For each coordination constraint violated, ejk = Tbackup,k − Tprimary,k − CTI; else ejk = 0; ˇ is a penalty factor. The value of 1000 was used. The best result for continuous DE algorithm is shown in Table 5. Table 6 shows statistical data for 10 runs of DE. Fig. 6 shows the evolution of the best OF value per iteration, considering the best 3 runs. From Table 5, it can be seen that the OF value using the BIP approach are lower than using the DE algorithm. To improve the result of the DE algorithm, the maximum iterations must be increased and the parameters of the algorithm must be correctly chosen by a trial-and-error procedure. On the other hand, the BIP has the ability to find the optimal discrete solution for this problem with a single execution using the default configuration of the CPLEX solver. Table 6 Statistical data of the results of DE algorithm. Parameter OF values
Value
112.7958; 113.6175; 113.6447; 116.2549; 115.2693; 115.9573; 115.3835; 116.9318; 114.8772; 120.5272 Mean value of OF 115.5259 2.1820 Standard deviation of OF Mean value of simulation time (s) 1054.49
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
Another situation was considered: the simulation time using the BIP approach and the DE algorithm with 200 individuals was limited to 180 s. The OF value of the BIP approach and the best of 10 runs using the DE algorithm were 116.9039 and 118.5659, respectively. The OF value of the BIP is the best solution found at the set time limit and it is better than the DE solution. The CPLEX gives a gap of 28.51% for this solution, which indicates the relative optimality. The gap can be seen as an estimate of the distance of the best feasible solution from the optimal. In CPLEX, it is calculated by [|best bound − best integer|/(1e−10 + |best integer|)] × 100, where best bound is the best objective function value achievable and best integer is the objective function value of the integer solution at the time of termination [28]. As we can see, the truncated solution of the BIP is better than the DE at the set time limit. This approach is useful when we are handling with large power systems, where there might be an increase in computational time required to prove the optimality. Therefore, we can consider truncated runs of CPLEX as a strategy to find good solutions within acceptable times for online purposes. 5. Conclusion The coordination of DOCR in interconnected power systems with distributed generation was formulated as a BIP problem and solved using the CPLEX package. The variables of the problem are discrete time dial and pickup current settings, thus electromechanical and microprocessor-based relays can be modeled, overcoming the disadvantage of techniques that determine continuous TDS that may not be used in some relays. The coordination is performed online on each change in the topology of the network through a SCADA system, and the new settings are remotely adjusted. A pre-processing step has been proposed to reduce the number of variables in the BIP model as well as to consider the melting curve of cables as a new constraint. Simulations were carried out in three test systems and the results were compared to heuristic-based techniques recently used, showing that the BIP approach gives better discrete solutions than continuous solutions presented in the literature. The deterministic solver running on its default configuration dismisses the adjustment of any algorithm parameters, which are common in heuristic-based techniques. Additionally, an optimal solution is found with a single execution of the algorithm in reduced times, as required for online coordination purposes. References ˜ On-line coordination of directional [1] M.Y. Shih, A.C. Enríquez, L.M.T. Trevino, overcurrent relays: performance evaluation among optimization algorithms, Electr. Power Syst. Res. 110 (2014) 122–132. [2] F. Coffele, C. Booth, A. Dysko, An adaptive overcurrent protection scheme for distribution networks, IEEE Trans. Power Delivery 30 (2) (2015) 561–568. [3] A.J. Urdaneta, R. Nadira, L.G.P. Jiménez, Optimal coordination of directional overcurrent relays in interconnected power systems, IEEE Trans. Power Delivery 3 (3) (1988) 903–910. [4] A.J. Urdaneta, H. Restrepo, S. Márquez, J. Sánchez, Coordination of directional lvercurrent relay timing using linear programming, IEEE Trans. Power Delivery 11 (1) (1996) 122–129.
125
[5] A.J. Urdaneta, L.G. Pérez, H. Restrepo, Optimal coordination of directional overcurrent relays considering dinamic changes in the network topology, IEEE Trans. Power Delivery 12 (4) (1997) 1458–1464. [6] A.Y. Abdelaziz, H.E.A. Talaat, A.I. Nosseir, A.A. Hajjar, An adaptive protection scheme for optimal coordination of overcurrent relays, Electr. Power Syst. Res. 61 (2002) 1–9. [7] H. Zeineldin, E.F. El-Saadany, M.A. Salama, A novel problem formulation for directional overcurrent relay coordination, in: 2004 Large Eng. Systems Conf. on Power Eng., 2004, pp. 48–52. [8] H. Zeineldin, E.F. El-Saadany, M.A. Salama, Optimal coordination of directional overcurrent relay coordination, in: 2005 Power Eng. Society General Meeting, vol. 2, 2005, pp. 1101–1106. [9] D. Birla, R.P. Maheshwari, H.O. Gupta, An approach to tackle the threat of sympathy trips in directional overcurrent relay coordination, IEEE Trans. Power Delivery 22 (2) (2007) 851–858. [10] D. Birla, R.P. Maheshwari, H.O. Gupta, A new nonlinear directional overcurrent relay coordination technique, and Banes and Boons of near-end faults based approach, IEEE Trans. Power Delivery 21 (3) (2006) 1176–1182. [11] C.W. So, K.K. Li, K.T. Lai, K.Y. Fung, Application of genetic algorithm for overcurrent relay coordination, in: International Conf. of Developments in Power Syst. Protection, 1997, pp. 66–69. [12] 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 (2008) 713–720. [13] R.M. Chabanloo, H.A. Abyaneh, S.S.H. Kamangar, F. Razavi, Optimal combined overcurrent and distance relays coordination incorporing intelligent overcurrent relays characteristics selection, IEEE Trans. Power Delivery 26 (3) (2011) 1381–1391. [14] R. Mohammadi, H.A. Abyaneh, H.M. Rudsari, S.H. Fathi, H. Rastegar, Overcurrent relays coordination considering the priority of constraints, IEEE Trans. Power Delivery 26 (3) (2011) 1927–1938. [15] D. Solati, M.R. Vatani, M.J. Sanjari, G.B. Gharehpetian, A.H. Yatim, Overcurrent relays coordination in interconnected networks using accurate analytical method and based on determination of fault critical point, IEEE Trans. Power Delivery 30 (2) (2015) 870–877. [16] Z. Moravej, F. Adelnia, F. Abbasi, Optimal coordination of directional overcurrent relays using NSGA-II, Electr. Power Syst. Res. 119 (2015) 228–236. [17] H.H. Zeineldin, E.F. El-Saadany, M.M.A. Salama, Optimal coordination of overcurrent relays using a modified particle swarm optimization, Electr. Power Syst. Res. 76 (2006) 988–995. [18] M.M. Mansour, S.F. Mekhamer, N.E. El-Kharbawe, A modified particle swarm optimizer for the coordination of directional overcurrent relays, IEEE Trans. Power Delivery 22 (3) (2007) 1400–1410. [19] T. Amraee, Coordination of directional overcurrent relays using seeker algorithm, IEEE Trans. Power Delivery 27 (3) (2012) 1415–1422. [20] T.R. Chelliah, R. Thangaraj, S. Allamsetty, M. Pant, Coordination of directional overcurrent relays using opposition based chaotic differential evolution algorithm, Electr. Power Energy Syst. 55 (2014) 341–350. [21] J. Moirangthem, K.R. Krishnanand, S.S. Dash, R. Ramaswami, Adaptive differential evolution algorithm for solving non-linear coordination problem of directional overcurrent relays, IET Gener. Transm. Distrib. 7 (4) (2013) 329–336. [22] M. Singh, B.K. Panigrahi, A.R. Abhyankar, Optimal coordination of directional over-current relays using teaching learning-based optimization (TLBO) algorithm, Electr. Power Energy Syst. 50 (2013) 33–41. [23] A.S. Noghabi, J. Sadeh, H.R. Mashhadi, Considering different network topologies in optimal overcurrent relay coordination using a hybrid GA, IEEE Trans. Power Delivery 24 (4) (2009) 1857–1863. [24] P.P. Bedekar, S.R. Bhide, Optimum coordination of directional overcurrent relays using the hybrid GA-NLP approach, IEEE Trans. Power Delivery 26 (1) (2011) 109–119. [25] C.A.C. Salazar, A.C. Enríquez, S.E. Schaeffer, Directional overcurrent relay coordination considering non-standardized time curves, Electr. Power Syst. Res. 122 (2015) 42–49. [26] H.M. Sharaf, H.H. Zeineldin, D.K. Ibrahim, E.E.A. El-Zahab, A proposed coordination strategy for meshed distribution systems with DG considering user-defined characteristics of directional inverse time overcurrent relays, Electr. Power Energy Syst. 65 (2015) 49–58. [27] F.S. Hillier, G.J. Lieberman, Introduction to Operations Research, seventh ed., McGraw-Hill, New York, NY, 2001. [28] IBM ILOG CPLEX, IBM ILOG CPLEX Optimization Studio V12.4 Documentation, 2015, Available: http://www-01.ibm.com/support/knowledgecenter/SSSA5P 12.4.0/maps/welcome.html (accessed: May, 2015).
Contents lists available at ScienceDirect
Electric Power Systems Research journal homepage: www.elsevier.com/locate/epsr
Online coordination of directional overcurrent relays using binary integer programming Rafael Corrêa a,∗ , Ghendy Cardoso Jr. b , Olinto C.B. de Araújo b , Lenois Mariotto b a b
Federal Institute of Education, Science and Technology from Rio Grande do Sul, Farroupilha, RS, Brazil Federal University of Santa Maria, Santa Maria, RS, Brazil
a r t i c l e
i n f o
Article history: Received 7 October 2014 Received in revised form 21 May 2015 Accepted 23 May 2015 Keywords: Binary integer programming Directional overcurrent relay Online coordination SCADA system Smart grid Optimization
a b s t r a c t This work presents a binary programming model for online coordination of directional overcurrent relay problems in interconnected power distribution and subtransmission systems. The proposed model considers discrete time dial and pickup current settings as in microprocessor-based relays. The coordination problem is solved for every topological change in the network and whose optimized settings are remotely adjusted on each relay through a SCADA system, enhancing the speed and sensitivity of the protection system. A pre-processing step is performed to ensure that the operation time of each relay be below the melting curve of the protected cable. Computational results for the binary programming model solved by CPLEX optimization package are compared to heuristic-based techniques. The simulation results indicates that our approach is suitable for online applications, since an optimal solution can be found in reduced times. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Overcurrent relays are used as primary protection for distribution and subtransmission lines as well as backup protection for transmission lines, power transformers and generators. In fault situations, when the primary relay does not operate, a remote backup relay must operate and the associated circuit breaker eliminates the fault. In order to guarantee that the backup relay does not operate before the circuit breaker associated to the primary relay, an intentional time delay is considered to adjust the backup. This process is called coordination. This protection philosophy is used in power distribution and subtransmission networks. In interconnected power systems and considering distributed generation (DG), a known adversity in coordinating the non-directional overcurrent protection requires the application of directional overcurrent relays (DOCR). In some cases, selecting the appropriate settings for DOCR to ensure coordination, speed and sensitivity for the protection system may be a challenge. For instance, it is common that minimum fault currents at the end of the backup zone be lower than the maximum load current. In this situation, when a pickup current above
∗ Corresponding author. Tel.: +55 54 32602400. E-mail addresses: [email protected] (R. Corrêa), [email protected] (G. Cardoso Jr.), [email protected] (O.C.B.d. Araújo), [email protected] (L. Mariotto). http://dx.doi.org/10.1016/j.epsr.2015.05.017 0378-7796/© 2015 Elsevier B.V. All rights reserved.
the maximum load current is chosen, the relay will not provide backup protection for such fault currents. In an offline coordination process, considering the different n − 1 contingency topologies, the obtained settings may lead to high operation times for some topologies. To overcome this issue [1,2] propose to re-coordinate all DOCR for every topological change in the network. This concept is also used in this study, considering the possibilities and the growth of smart grids. The coordination of DOCR has been formulated as an optimization problem and solved with several techniques. To determine continuous time dial settings (TDS), linear programming solvers have been used in [3–6] considering fixed pickup current settings (PCS). In [7,8], mixed integer non-linear programming and mixed integer programming solvers, respectively, have been used to determine discrete PCS and continuous TDS. Continuous TDS and PCS have been determined using a non-linear programming solver in [9,10]. Evolutionary techniques have been used, such as the genetic algorithm (GA) [11–15], nondominated sorting genetic algorithmII (NSGA-II) [16], particle swarm optimization (PSO) [17,18], seeker optimization algorithm (SOA) [19], ant colony algorithm [1], differential evolution (DE) algorithm [1], opposition based chaotic differential evolution (OCDE) algorithm [20], adaptive differential evolution (ADE) algorithm [21] and teaching learning-based optimization algorithm [22]. Hybrid techniques have also been used in [23,24], combining metaheuristics with deterministic solvers. The main drawback from these heuristic-based techniques is the
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
presence of parameters that must be correctly chosen in order to provide the convergence of the algorithm. The optimal parameters may not be the same for different instances, thus this choice is performed in a trial-and-error procedure. Also, given the non deterministic nature of such algorithms, more than one simulation may be necessary to achieve the optimal. Nonconventional, or non-standardized, characteristics for inverse DOCR have been considered in [25,26]. The GA and a nonlinear programming problem solver, respectively, have been used to optimize TDS, PCS and the shape – characteristic – of the curve of each DOCR, taking advantage of numerical relays. In general, when the coordination of DOCR is addressed using the techniques previously listed, the TDS are determined as a continuous variable. However, if the step of the TDS of a microprocessor-based relay is not small enough, such as 0.01, miscoordinations may occur whether we round the continuous solution to the nearest available settings [12]. To overcome the described issues, this paper presents a binary integer programming (BIP) model considering discrete TDS and PCS, within the available range. The online coordination of DOCR for three test systems in presence of DG is performed using the CPLEX optimization package. In order to show the effectiveness of the proposed approach, the results are compared with recently used heuristic-based techniques. The novelties of this paper are described below: • The coordination of DOCR problem is modeled using a BIP model and solved by a mixed integer programming solver running on its default configuration, i.e., the user does not need to set parameters. Also, a single simulation is necessary to achieve the optimal; • The TDS and PCS are determined in the discrete form, as available in electromechanical and microprocessor-based units with noncontinuous settings; • A pre-processing step is performed to ensure that the operating times of the DOCR be below the melting curve of the protected cables. This step reduces the simulation times, making the proposed approach suitable for online applications. The other sections of this paper are organized as follows: in Section 2, the conventional formulation of the DOCR problem is described. The proposed online coordination procedure is detailed in Section 3. In Section 4, the proposed approach is applied over three test systems and the simulation results are discussed. Finally, the conclusions are provided in Section 5.
119
The constraints of the problem are defined as follows. 2.1. Limits of the settings The bounds of the settings can be expressed by TDSmin ≤ TDSi ≤ TDSmax , i i PCSmin i
≤ PCSi ≤
PCSmax , i
i = 1, . . ., m i = 1, . . ., m
(3) (4)
where TDSi min and TDSi max are the lower and upper limits of TDS of the ith relay, respectively; PCSi min and PCSi max are the lower and upper limits of PCS of the ith relay, respectively. Moreover, for phase protection, the PCS must be greater than maximum load current and lower than minimum fault current seen by each relay, including a safety margin, which depends on the relay technology and current transformer (CT) errors. 2.2. Limits of relay operation time Minimum and maximum operating times can be established for each relay by min max Tik ≤ Tik ≤ Tik , min
i = 1, . . ., m
(5)
max
and Tik are the minimum and maximum operwhere Tik ation times of the ith relay for a fault at k, respectively. 2.3. Coordination of primary and backup relays The backup relays must operate with an intentional time delay, called coordination time interval (CTI), in relation to the primary relay. Thereby, the backup operates only when the primary protection fails. This time interval includes factors such as the circuit breaker operation time, CT errors, relay overtravel time and a safety margin. Its value is usually selected between 0.2 s and 0.5 s. Tjk − Tik − CTI ≥ 0,
i = 1, . . ., m
(6)
where Tjk is the operation time of the jth backup relay for a fault at k inside the zone protected by the ith primary relay. This constraint must be satisfied for all primary/backup pairs. Usually, near-end and far-end fault currents are used, providing coordination for most of the fault situations. 2.4. Topological changes of the network
2. Conventional problem formulation
The coordination constraint (6) is considered not only to the main topology of the network, but also to other possible configurations.
The objective of the DOCR coordination problem is to minimize the operating times of the relays operating as primary protection. The objective function (OF) is then expressed by
3. Proposed approach
min OF =
PCSi ,TDSi
m i=1
Tik
(1)
k
where m is the number of relays in the network; PCSi and TDSi are the pickup current and time dial settings of the ith relay, respectively; and Tik is the operation time of the ith relay for a fault at k, defined as follows.
Tik = TDSi ×
A Iik /PCSi
P
−1
+B
(2)
where Iik is the fault current seen by the ith relay for a fault at k; A, B and P are standardized coefficients. In this work, the Type A curve from the IEC 255-3 has been used, where A = 0.14, B = 0 and P = 0.02.
This section presents the proposed online coordination procedure. The main steps are the supervisory control and data acquisition (SCADA) system, the pre-processing step and the process to create and solve the BIP model. 3.1. SCADA system The main tasks of the SCADA system are to run load flow and short-circuit analysis for each topological change of the power network, to optimize the TDSs and PCSs, and to re-set each relay remotely. Although the communication between a remote device and a supervisory system is already a reality in the state-of-the-art of digital relays, a fully implemented system to execute such resetting in interconnected power distribution and subtransmission networks is not available. Thus, this re-setting system is supposed to have been implemented.
120
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
According to each topological change of the network, the SCADA system processes the new topology and the following steps are performed: • The coordination pairs of relays are determined and stored as data; • Load flow analysis is executed for the n − 1 contingency topologies. The resulting maximum load current for each relay is stored as data; • Short-circuit analysis is executed considering the new topology. The fault currents related to each coordination pair of relays are stored as data.
Afterwards, the maximum load and fault currents are used to construct the BIP model. This model is solved and the optimal settings are adjusted on each relay via communication hardware. By the proposed approach, relays are coordinated for the current topology of the network. Whether any change occurs, the online coordination process starts, and coordination may be lost if a fault occurs before the new settings are determined and sent to each relay. However, the DOCR do not operate if a fault does not occur, since the pickup currents are higher than maximum load currents considering the n − 1 criterion. It is important to note that many types of relays cannot be remotely adjusted, like electromechanical units. If they require their settings to be changed, they must be locally adjusted. For example, in an expansion of the network, maximum load current may overcome the pickup current, thus the PCS must be increased. Therefore, two solutions are proposed:
(1) If these units require their settings to be changed, the BIP model is constructed considering fault currents for all possible topologies of the network. These units are then locally adjusted and digital relays are remotely adjusted. (2) If these units do not require their settings to be changed, the BIP model is constructed with fixed settings for these units, reducing the number of variables of the problem.
3.2. Binary integer programming A BIP model has been used to deal with TDS and PCS as discrete variables. The BIP is an integer programming problem with binary variables. A multivariable BIP problem can be written in the matrix form as follows. min FO = cT x s.t. Ax ≤ b
(7)
x ∈ Bn where cT is the transposed cost vector; x is the variables vector; A is a matrix with constants; b is a vector with constants; and Bn represents the space of vectors with n binary components. The branch and bound and the branch and cut are well known algorithms used for solving BIP problems. Commercial solvers for BIP problems incorporate both classical and recently developed techniques. For instance, the MIPOPT solver from CPLEX offers a dynamic search algorithm, which consists in LP relaxation, branching, cuts, and heuristics. Since it is not possible to detail all techniques implemented inside the CPLEX, the reader may find the description of the branch and bound and the branch and cut algorithms, and how they work, in [27].
3.3. Representation of the variables To represent TDS and PCS in the discrete form and to linearize the constraints and the objective function of the model presented in Section 2, the following changes have been made in the conventional problem formulation. A constant cirsk has been introduced, considering discrete TDSs and PCSs. Its value is calculated as follows.
cirsk = TDSir ×
A Iik /PCSis
P
−1
+B
(8)
where TDSir is the rth element of the set vTDSi . This set contains all nTDSi values of TDS for the ith relay; PCSis is the sth element of the set vPCSi . This set contains all nPCSi values of PCS for the ith relay; nTDSi is the number of available values of TDS for the ith relay; nPCSi is the number of available values of PCS for the ith relay; Iik is the fault current seen by the ith relay for a fault at k; A, B and P are standardized coefficients of the selected curve. The relay operating time is then written using (8) and binary variables xirs in the form
nTDSi nPCSi
Tik =
cirsk xirs
(9)
r=1 s=1
where xirs is a binary variable: equal to 1 if the rth TDS value and the sth PCS value are assigned to the ith relay; and 0 otherwise. Each variable xirs represents a combination of a TDS and a PCS for the ith relay. To ensure a single combination of TDS and PCS has been assigned to each relay, the following constraints are included in the model.
nTDSi nPCSi
xirs = 1,
i = 1, . . ., m
(10)
r=1 s=1
xirs ∈
0, 1 ,
i = 1, . . ., m;
r = 1, . . ., nTDSi ;
s = 1, . . ., nPCSi (11)
In (9), the relay operating time is linear with respect to xirs . The objective function (1) as well as the constraints (5) and (6) are then written using (8) and (9), and the resulting problem is linear. The sets vTDSi and vPCSi contain all possible values of TDS and PCS for the ith relay, respectively. For example, if relay number 3 has TDS available between 0.1 and 0.4 in steps of 0.1, the set is vTDS3 = {0.1; 0.2; 0.3; 0.4}. Using this representation, it is possible to model electromechanical and microprocessor-based relays. In electromechanical units, usually the PCS do not have a constant step, thus it can be written as vPCSi = {0.5; 0.6; 0.8; 1.0; 1.5; 2.0; 2.5}. Also, a step of 0.01 can be used to represent continuous settings of digital relays, without prejudice in the speed and sensitivity. As mentioned previously, each variable xirs represents a combination of a TDS and a PCS for the ith relay. For example, if the problem has 10 relays and each one presents 50 available TDSs and 40 available PCSs, then the problem has 10 × 50 × 40 = 20,000 variables. Similarly, if each relay has 50 values for TDS between 0.1 and 1.1 and 40 values for PCS between 0.5 and 4.0, the variable x5,1,1 represents a TDS 0.1 and a PCS 0.5 for relay number 5; the variable x2,1,40 represents a TDS 0.1 and PCS 4.0 for relay number 2; and so on. 3.4. Pre-processing step In real DOCR, hundreds of settings are available and the coordination problem has a large number of variables, as in the proposed
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
121
Fig. 1. Technical constraints of the pre-processing step.
BIP model each variable represents a combination of TDS and PCS. Although such a problem can be solved in the actual state-of-theart of the mixed integer programming solvers, in some cases the simulation times may be prohibitive for online purposes. In order to reduce the simulation time by reducing the number of variables, a pre-processing step is proposed to eliminate variables that do not match the following technical constraints: • The available PCS must be within the range defined in Section 2.1; • The relay operating curve must be below the melting curve of the protected cable; • The relay operating time cannot be greater than 5 s for faults in the primary zone; • The relay operating time cannot be greater than 10 s for faults in the remote backup zone. These technical constraints are shown in Fig. 1. The choice of maximum values for the operating times as primary and remote backup protection is based on experience. Before constructing the BIP model, all combinations of TDS and PCS of each relay are tested considering the technical constraints. The pairs that do not match these constraints are dismissed and do not compose the BIP model. For instance, in Fig. 1 it can be observed that if the TDS is increased, from a certain value the operating curve of the relay overcomes the melting curve of the cable. Therefore, any TDS above this limit does not compose the final model for that particular PCS value. It is important to note that without this pre-processing the variables that would be eliminated at this step are included in the BIP model, which could lead to higher simulation times due to the larger number of variables. Also, this step works as a constraint related to the melting curve of the cables, as the variables that could violate this constraint are not included in the final model. The flowchart of the pre-processing step is shown in Fig. 2.
Fig. 2. Flowchart of pre-processing step.
4.1. IEEE 6-bus test system The IEEE 6-bus test system is shown in Fig. 3. Three-phase fault currents of the main topology shown in Fig. 3 and CT ratios can be found in [20]. The available PCS are from 1.25 to 1.50, in steps of 0.01, and the available TDS are from 0.05 to 1.10, in steps of 0.01. Also, a CTI 0.2 s has been used. In [20], simulation results for DE algorithm and an OCDE algorithm, named OCDE2, have been presented. To make a clear comparison, the objective function of the BIP approach has been defined as in [20], considering near-end and far-end fault currents. The results are shown in Table 1. The simulation time using the BIP approach was 13.08 s. From Table 1, it can be observed that the OF using the BIP approach is lower than the DE algorithm, even with a worse TDS and PCS resolution. The OF using OCDE2 algorithm is lower than the proposed technique, but the continuous settings may not be available in some relays. On the contrary, the optimized settings
4. Results and discussion The BIP approach was applied in three test systems and the results have been compared to heuristic-based techniques. The coordination was carried out for only phase relays, since the ground relay coordination uses the same procedure with ground fault currents. In addition, the instantaneous units were not considered. Algorithms in C++ have been implemented to run the CPLEX 12.4 mixed integer programming solver. They have been executed on a Core i3, 2.13 GHz PC with 4 GB RAM.
Fig. 3. Single-line diagram of the IEEE 6-bus test system.
122
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
Table 1 Results for IEEE 6-bus test system. Relay
DE [20]
Table 2 Results of Case 1 for 8-bus test system. OCDE2 [20]
BIP
Relay
TDS
PCS
TDS
PCS
TDS
PCS
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.1173 0.2082 0.0997 0.1125 0.0500 0.0580 0.0500 0.0500 0.0500 0.0719 0.0649 0.0617 0.0500 0.0856
1.2505 1.2500 1.2512 1.2515 1.2500 1.2500 1.2500 1.2500 1.2502 1.2502 1.4998 1.2575 1.4805 1.2557
0.1014 0.1862 0.0946 0.1006 0.0500 0.0500 0.0500 0.0500 0.0500 0.0562 0.0649 0.0508 0.0500 0.0708
1.5000 1.4999 1.2772 1.4995 1.2501 1.3801 1.2500 1.2500 1.2500 1.4988 1.4999 1.5000 1.4608 1.4995
0.11 0.19 0.10 0.11 0.05 0.05 0.05 0.05 0.05 0.07 0.07 0.06 0.05 0.08
1.38 1.46 1.25 1.31 1.25 1.41 1.25 1.25 1.25 1.25 1.44 1.29 1.44 1.38
OF (s)
10.6272
10.3286
10.5380
Hybrid GA-LP [23]
SOA [19]
TDS
PCS
TDS
PCS
BIP TDS
PCS
1 2 3 4 5 6 7 8 9 10 11 12 13 14
0.3043 0.2917 0.2543 0.1851 0.1700 0.2711 0.5316 0.2387 0.1865 0.1895 0.2014 0.2890 0.2207 0.5278
1.0 2.5 2.5 2.5 1.5 2.5 0.5 2.5 2.0 2.5 2.5 2.5 1.5 0.5
0.113 0.260 0.225 0.160 0.100 0.173 0.243 0.170 0.147 0.176 0.187 0.266 0.114 0.246
2.0 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.0 2.5
0.10 0.28 0.24 0.19 0.10 0.18 0.26 0.17 0.15 0.18 0.19 0.27 0.10 0.25
2.5 2.5 2.5 2.0 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5
OF (s)
10.9499
8.4270
8.6944
using the BIP approach are discrete and within the available values for each relay.
currents, with emphasis between relays number 11 and 12 where the left side of the coordination constraint (6) is equal to −0.111.
4.2. 8-Bus test system
4.2.2. Case 2 In Case 2, the line where relays number 1 and 8 are located was removed and the coordination was performed. These devices do not compose the BIP model. This simulation case is to show the effectiveness of the BIP approach dealing with line outages and different topologies of the network. The simulation time of the proposed BIP model was 38.75 s. Table 3 shows the simulation results. It can be seen that relays number 1, 8 and 13 have the highest available value for TDS. These relays do not operate for any fault on this topology, thus any TDS can be used for these units. Once the PCS are optimized considering the maximum load currents, these relays do not trip after a change in topology if a fault does not occur. After this change, new settings must be determined to ensure the coordination in fault situations. Additionally, relay number 6 has the lowest available TDS, as it is the downstream device for coordination purposes in this topology.
The 8-bus test system is shown in Fig. 4. The data of the network can be found in [17]. A CTI 0.3 s and cables 1/0 CAA have been considered. The available TDSs are from 0.10 and 1.10, in steps of 0.01, and the available PCSs are 0.5, 0.6, 0.8, 1.0, 1.5, 2.0 and 2.5. Three cases have been simulated, as described below. 4.2.1. Case 1 Hybrid GA-LP [23] and SOA [19] considered near-end threephase fault currents of the main topology shown in Fig. 4. The same fault currents have been used in the BIP approach. The results of Case 1 are shown in Table 2. The simulation time using the BIP approach was 34.11 s. Without the pre-processing procedure, the time was 55.30 s. It is important to note that a single execution of the algorithm is required to reach an optimal solution, which is not always true considering heuristic techniques. From Table 2, it can be noted that the OF using the BIP approach is lower than the hybrid GA-LP, even with a worse TDS resolution. The OF using SOA is lower than the proposed technique, but the TDS are continuous. To illustrate the miscoordinations that may occur if we round these settings to the nearest available value, we considered the SIPROTEC® 7SJ46 numerical overcurrent relay, which has a step of 0.05. Nine miscoordinations occur for near-end three-phase fault
4.2.3. Case 3 An outage of generator G1 was considered in Case 3. In this situation, all relays compose the BIP model and fault current levels are lower than in Case 1. The simulation time of the BIP approach was 22.23 s. Table 3 shows the simulation results. It can be seen that lower TDS values have been determined for most relays in comparison with Case 1. Table 3 Results of Case 2 and Case 3 for 8-bus test system. Relay
Fig. 4. Single-line diagram of the 8-bus test system.
Case 2
Case 3
TDS
PCS
TDS
PCS
1 2 3 4 5 6 7 8 9 10 11 12 13 14
1.10 0.32 0.33 0.23 0.16 0.10 0.35 1.10 0.24 0.25 0.25 0.33 1.10 0.36
2.0 2.5 2.0 2.5 2.5 0.5 2.5 2.0 2.5 2.5 2.5 2.5 2.0 2.5
0.11 0.15 0.13 0.13 0.12 0.16 0.22 0.10 0.11 0.14 0.12 0.21 0.10 0.10
2.0 2.5 2.5 2.5 1.5 2.5 2.5 2.0 2.5 2.5 2.5 2.0 0.5 1.5
OF (s)
9.0747
7.3972
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
123
Table 4 RTC and maximum load current for IEEE 30-bus test system.
Fig. 5. Single-line diagram of the IEEE 30-bus test system.
Once the inverse characteristic curves have been considered, lower TDS are expected to ensure the speed of the protection system for lower fault currents.
Relay
RTC
max Iload (A)
Relay
RTC
max Iload (A)
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
150 150 100 100 100 100 100 150 100 100 100 50 150 100 100 100 100 50 100 50 100 100 150 100 150 150 150 150 150 100 100 150 50 100 100 250 150 150 150
637 637 470 0 0 361 431 519 0 110 483 0 561 0 0 202 209 0 66 144 225 73 461 0 452 345 235 235 244 11 0 539 0 460 0 1086 288 451 300
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78
50 50 100 50 100 100 50 100 50 50 100 50 50 50 50 100 50 50 100 50 100 100 100 100 50 100 50 50 100 100 100 50 50 50 50 100 100 100 50
0 191 110 183 278 258 228 304 134 0 307 81 153 220 113 220 162 158 194 98 0 0 108 347 173 269 220 0 445 244 247 0 201 47 0 0 350 119 0
4.3. IEEE 30-bus test system The IEEE 30-bus test system, shown in Fig. 5, has 78 DOCRs and 206 coordination pairs in the main topology. The RTC and the maximum load current for each relay are given in Table 4.
Fig. 6. Convergence of the DE algorithm for the IEEE 30-bus test system.
124
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125 Table 5 (Continued)
Table 5 Results for IEEE 30-bus test system. Relay
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 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
Relay
BIP
BIP
DE [1]
TDS
PCS
TDS
PCS
0.01 0.01 0.02 0.06 0.06 0.04 0.07 0.06 0.10 0.04 0.01 0.17 0.03 0.09 0.10 0.07 0.12 0.22 0.02 0.09 0.10 0.08 0.05 0.01 0.10 0.07 0.07 0.07 0.14 0.05 0.08 0.08 0.04 0.03 0.11 0.05 0.07 0.08 0.04 0.12 0.08 0.05 0.08 0.11 0.10 0.12 0.07 0.12 0.13 0.08 0.10 0.10 0.15 0.07 0.07 0.09 0.08 0.05 0.07 0.04 0.04 0.08 0.04 0.06 0.08 0.09 0.11 0.01 0.02 0.02 0.01 0.01 0.03 0.01
1.60 1.60 2.70 0.80 0.80 1.80 1.10 0.80 0.50 2.70 2.50 0.60 1.90 0.50 1.50 1.20 1.20 0.50 1.90 1.50 0.50 0.50 0.80 0.50 2.40 3.00 2.00 2.00 0.60 2.60 0.50 1.10 1.50 1.90 0.60 1.60 3.00 3.00 2.90 0.70 2.30 2.90 2.10 2.80 1.00 3.00 3.00 1.60 1.00 3.00 2.90 2.80 2.60 2.30 3.00 2.30 2.70 2.90 2.80 0.60 0.60 3.00 2.80 2.60 0.60 2.90 1.00 1.80 1.90 2.30 0.50 2.90 1.50 0.50
0.0296 0.0346 0.0665 0.0793 0.0460 0.1169 0.0923 0.0805 0.0386 0.1117 0.0284 0.1115 0.0770 0.0953 0.1190 0.1713 0.1596 0.2000 0.1248 0.1550 0.0972 0.0579 0.0451 0.0100 0.2052 0.0857 0.1512 0.1469 0.1354 0.1383 0.0584 0.1147 0.0542 0.0630 0.1268 0.0684 0.0812 0.1765 0.1013 0.1953 0.1861 0.1519 0.0987 0.2485 0.1131 0.2386 0.1355 0.1521 0.1686 0.1066 0.1802 0.2351 0.2247 0.2365 0.1516 0.1793 0.2290 0.1669 0.1665 0.0656 0.0573 0.2794 0.1803 0.1457 0.1046 0.2021 0.1056 0.0653 0.0346 0.0560 0.0100 0.0474 0.0327 0.0100
1.1954 1.1085 1.7014 0.8681 1.5224 0.8042 0.9314 0.7609 1.4716 1.2382 1.6166 1.4398 0.8632 0.5685 1.3892 0.5702 0.9721 0.7917 0.5376 0.8066 0.6047 0.6794 0.8390 0.5000 1.0010 2.9939 1.3237 1.3724 1.1047 0.5266 0.7977 0.8001 1.3768 1.3302 0.7478 1.3484 2.9993 1.6597 1.8395 0.5001 1.1787 1.0852 1.7317 1.1345 1.1963 1.1189 1.6324 1.4450 0.9347 2.9986 1.5700 1.1836 1.7144 0.5279 1.7245 1.3946 0.7000 1.2399 0.8556 0.5001 0.5416 0.6463 0.7053 1.4774 0.9436 1.4395 1.4336 0.9122 1.4313 1.2673 0.5000 0.9583 1.3372 0.5000
DE [1]
TDS
PCS
TDS
PCS
75 76 77 78
0.08 0.07 0.05 0.11
0.50 2.70 2.50 0.50
0.0697 0.1789 0.1367 0.1051
0.6281 1.4193 0.5001 0.8290
OF (s) Simul. time (s)
89.2199 284.06
112.7958 1041.49
A CTI 0.2 s, cables 336.4 CAA and a large number of settings have been considered: the TDS is from 0.01 to 0.50, in steps of 0.01, and the PCS is from 0.5 to 3.0, in steps of 0.1. Also, the maximum operating time of the primary protection was reduced to 2 s. The BIP approach was applied in the main topology of the network. The following fault currents were considered to compose the coordination constraints and the objective function: three-phase near-end; and far-end with the remote circuit breaker opened. The far-end fault currents with the remote circuit breaker closed were also considered to compose the objective function. The results for this case are shown in Table 5. Using the preprocessing step, the simulation time was 284.06 s and the number of variables was 46,598. This step takes only 0.234 s. Without the pre-processing, the simulation time was 5593.07 s and the number of variables was 79,550. The reduction of 41.4% in the number of variables is due to the maximum operating time of the primary relay (2 s) considered in the pre-processing. The maximum operating time of the primary relay of 5 s was also used: the number of variables was 66,663 (reduction of 16.2%) and the pre-processing time was 0.343 s. For comparison purposes, the continuous DE algorithm used in [1] were executed 10 times with a population size of 200 individuals and 10,000 maximum iterations. The objective function used in [1] is different from (1), so we considered (12) to make a clear comparison with the BIP results. min OF =
PCSi ,TDSi
i
k
Tik + ˇ ×
j
|ejk |
(12)
k
where Tik is the operation time of the ith relay for a fault at k; ejk is the CTI error between the jth coordination pair for a fault at k. For each coordination constraint violated, ejk = Tbackup,k − Tprimary,k − CTI; else ejk = 0; ˇ is a penalty factor. The value of 1000 was used. The best result for continuous DE algorithm is shown in Table 5. Table 6 shows statistical data for 10 runs of DE. Fig. 6 shows the evolution of the best OF value per iteration, considering the best 3 runs. From Table 5, it can be seen that the OF value using the BIP approach are lower than using the DE algorithm. To improve the result of the DE algorithm, the maximum iterations must be increased and the parameters of the algorithm must be correctly chosen by a trial-and-error procedure. On the other hand, the BIP has the ability to find the optimal discrete solution for this problem with a single execution using the default configuration of the CPLEX solver. Table 6 Statistical data of the results of DE algorithm. Parameter OF values
Value
112.7958; 113.6175; 113.6447; 116.2549; 115.2693; 115.9573; 115.3835; 116.9318; 114.8772; 120.5272 Mean value of OF 115.5259 2.1820 Standard deviation of OF Mean value of simulation time (s) 1054.49
R. Corrêa et al. / Electric Power Systems Research 127 (2015) 118–125
Another situation was considered: the simulation time using the BIP approach and the DE algorithm with 200 individuals was limited to 180 s. The OF value of the BIP approach and the best of 10 runs using the DE algorithm were 116.9039 and 118.5659, respectively. The OF value of the BIP is the best solution found at the set time limit and it is better than the DE solution. The CPLEX gives a gap of 28.51% for this solution, which indicates the relative optimality. The gap can be seen as an estimate of the distance of the best feasible solution from the optimal. In CPLEX, it is calculated by [|best bound − best integer|/(1e−10 + |best integer|)] × 100, where best bound is the best objective function value achievable and best integer is the objective function value of the integer solution at the time of termination [28]. As we can see, the truncated solution of the BIP is better than the DE at the set time limit. This approach is useful when we are handling with large power systems, where there might be an increase in computational time required to prove the optimality. Therefore, we can consider truncated runs of CPLEX as a strategy to find good solutions within acceptable times for online purposes. 5. Conclusion The coordination of DOCR in interconnected power systems with distributed generation was formulated as a BIP problem and solved using the CPLEX package. The variables of the problem are discrete time dial and pickup current settings, thus electromechanical and microprocessor-based relays can be modeled, overcoming the disadvantage of techniques that determine continuous TDS that may not be used in some relays. The coordination is performed online on each change in the topology of the network through a SCADA system, and the new settings are remotely adjusted. A pre-processing step has been proposed to reduce the number of variables in the BIP model as well as to consider the melting curve of cables as a new constraint. Simulations were carried out in three test systems and the results were compared to heuristic-based techniques recently used, showing that the BIP approach gives better discrete solutions than continuous solutions presented in the literature. The deterministic solver running on its default configuration dismisses the adjustment of any algorithm parameters, which are common in heuristic-based techniques. Additionally, an optimal solution is found with a single execution of the algorithm in reduced times, as required for online coordination purposes. References ˜ On-line coordination of directional [1] M.Y. Shih, A.C. Enríquez, L.M.T. Trevino, overcurrent relays: performance evaluation among optimization algorithms, Electr. Power Syst. Res. 110 (2014) 122–132. [2] F. Coffele, C. Booth, A. Dysko, An adaptive overcurrent protection scheme for distribution networks, IEEE Trans. Power Delivery 30 (2) (2015) 561–568. [3] A.J. Urdaneta, R. Nadira, L.G.P. Jiménez, Optimal coordination of directional overcurrent relays in interconnected power systems, IEEE Trans. Power Delivery 3 (3) (1988) 903–910. [4] A.J. Urdaneta, H. Restrepo, S. Márquez, J. Sánchez, Coordination of directional lvercurrent relay timing using linear programming, IEEE Trans. Power Delivery 11 (1) (1996) 122–129.
125
[5] A.J. Urdaneta, L.G. Pérez, H. Restrepo, Optimal coordination of directional overcurrent relays considering dinamic changes in the network topology, IEEE Trans. Power Delivery 12 (4) (1997) 1458–1464. [6] A.Y. Abdelaziz, H.E.A. Talaat, A.I. Nosseir, A.A. Hajjar, An adaptive protection scheme for optimal coordination of overcurrent relays, Electr. Power Syst. Res. 61 (2002) 1–9. [7] H. Zeineldin, E.F. El-Saadany, M.A. Salama, A novel problem formulation for directional overcurrent relay coordination, in: 2004 Large Eng. Systems Conf. on Power Eng., 2004, pp. 48–52. [8] H. Zeineldin, E.F. El-Saadany, M.A. Salama, Optimal coordination of directional overcurrent relay coordination, in: 2005 Power Eng. Society General Meeting, vol. 2, 2005, pp. 1101–1106. [9] D. Birla, R.P. Maheshwari, H.O. Gupta, An approach to tackle the threat of sympathy trips in directional overcurrent relay coordination, IEEE Trans. Power Delivery 22 (2) (2007) 851–858. [10] D. Birla, R.P. Maheshwari, H.O. Gupta, A new nonlinear directional overcurrent relay coordination technique, and Banes and Boons of near-end faults based approach, IEEE Trans. Power Delivery 21 (3) (2006) 1176–1182. [11] C.W. So, K.K. Li, K.T. Lai, K.Y. Fung, Application of genetic algorithm for overcurrent relay coordination, in: International Conf. of Developments in Power Syst. Protection, 1997, pp. 66–69. [12] 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 (2008) 713–720. [13] R.M. Chabanloo, H.A. Abyaneh, S.S.H. Kamangar, F. Razavi, Optimal combined overcurrent and distance relays coordination incorporing intelligent overcurrent relays characteristics selection, IEEE Trans. Power Delivery 26 (3) (2011) 1381–1391. [14] R. Mohammadi, H.A. Abyaneh, H.M. Rudsari, S.H. Fathi, H. Rastegar, Overcurrent relays coordination considering the priority of constraints, IEEE Trans. Power Delivery 26 (3) (2011) 1927–1938. [15] D. Solati, M.R. Vatani, M.J. Sanjari, G.B. Gharehpetian, A.H. Yatim, Overcurrent relays coordination in interconnected networks using accurate analytical method and based on determination of fault critical point, IEEE Trans. Power Delivery 30 (2) (2015) 870–877. [16] Z. Moravej, F. Adelnia, F. Abbasi, Optimal coordination of directional overcurrent relays using NSGA-II, Electr. Power Syst. Res. 119 (2015) 228–236. [17] H.H. Zeineldin, E.F. El-Saadany, M.M.A. Salama, Optimal coordination of overcurrent relays using a modified particle swarm optimization, Electr. Power Syst. Res. 76 (2006) 988–995. [18] M.M. Mansour, S.F. Mekhamer, N.E. El-Kharbawe, A modified particle swarm optimizer for the coordination of directional overcurrent relays, IEEE Trans. Power Delivery 22 (3) (2007) 1400–1410. [19] T. Amraee, Coordination of directional overcurrent relays using seeker algorithm, IEEE Trans. Power Delivery 27 (3) (2012) 1415–1422. [20] T.R. Chelliah, R. Thangaraj, S. Allamsetty, M. Pant, Coordination of directional overcurrent relays using opposition based chaotic differential evolution algorithm, Electr. Power Energy Syst. 55 (2014) 341–350. [21] J. Moirangthem, K.R. Krishnanand, S.S. Dash, R. Ramaswami, Adaptive differential evolution algorithm for solving non-linear coordination problem of directional overcurrent relays, IET Gener. Transm. Distrib. 7 (4) (2013) 329–336. [22] M. Singh, B.K. Panigrahi, A.R. Abhyankar, Optimal coordination of directional over-current relays using teaching learning-based optimization (TLBO) algorithm, Electr. Power Energy Syst. 50 (2013) 33–41. [23] A.S. Noghabi, J. Sadeh, H.R. Mashhadi, Considering different network topologies in optimal overcurrent relay coordination using a hybrid GA, IEEE Trans. Power Delivery 24 (4) (2009) 1857–1863. [24] P.P. Bedekar, S.R. Bhide, Optimum coordination of directional overcurrent relays using the hybrid GA-NLP approach, IEEE Trans. Power Delivery 26 (1) (2011) 109–119. [25] C.A.C. Salazar, A.C. Enríquez, S.E. Schaeffer, Directional overcurrent relay coordination considering non-standardized time curves, Electr. Power Syst. Res. 122 (2015) 42–49. [26] H.M. Sharaf, H.H. Zeineldin, D.K. Ibrahim, E.E.A. El-Zahab, A proposed coordination strategy for meshed distribution systems with DG considering user-defined characteristics of directional inverse time overcurrent relays, Electr. Power Energy Syst. 65 (2015) 49–58. [27] F.S. Hillier, G.J. Lieberman, Introduction to Operations Research, seventh ed., McGraw-Hill, New York, NY, 2001. [28] IBM ILOG CPLEX, IBM ILOG CPLEX Optimization Studio V12.4 Documentation, 2015, Available: http://www-01.ibm.com/support/knowledgecenter/SSSA5P 12.4.0/maps/welcome.html (accessed: May, 2015).