# A Novel Proposed Algorithm to Enhance the Overcurrent Relays’ Performance in Active Distribution Networks

Authors: Sally El-Tawab; Hassan S. Mohamed; A. M. Abdel-Aziz
Journal: International Transactions on Electrical Energy Systems
 (2022)

---

## Abstract

This research addresses two topics: the first is the optimal coordination of overcurrent relays (OCRs), which has recently become a major challenge to ensure reliability and speed of relays’ operation, and the second is a modified circuit added to the OCR to vanquish the barriers caused by the ever-increasing participation of distributed energy resources (DERs), at a lower cost. For robust and reliable OCRs coordination, various types of optimization techniques are frequently used to find optimal OCRs settings in distribution networks (DNs). In this study, the improved wild horse optimization (IWHO) algorithm is used as a novel metaheuristic method for solving optimization issues in coordination of OCRs in a distribution network, for the first time. This technique is used to determine the optimal time multiplier setting (TMS) and optimal pickup setting (PS) of OCRs in order to minimize the overall operating time. The results of the proposed algorithm are compared with those obtained from other recent metaheuristic techniques. The obtained TMS and PS settings of each relay show the effectiveness of the proposed IWHO technique in terms of accuracy and speed in a medium scale radial network configuration. To replace the costly mitigation methods employed in the DN comprising DER, a simple circuit was proposed to be added to the OCR. Using the angle’s sign of positive sequence current component (PSCC) of fault current solely, this proposed circuit can detect the fault direction and assist the OCR in making the correct decision in fault detection and clearing with high penetration of DERs. IEEE-33 bus electrical power system has been used to validate the newly proposed protective solutions.


---

## Body

## 1. Introduction

Electric power distribution networks (EPDN) are generally radial in nature, and the current philosophy in planning, management, and control is based on the assumption of the existence of unidirectional power flows transmitted from the highest voltage levels to the consuming levels. These assumptions allow for the implementation of relatively simple and inexpensive protection schemes with which a selective operation of the protection system is achieved [1]. Overcurrent relays (OCRs) are devices extensively used in EPDN as primary and backup protection devices. Backup protective relays are only activated if the primary relay fails to operate within the time limit specified. The backup relay must follow the coordination time interval (CTI) after fault initiation in its protection zone [2]. To enhance the fast operation and reliability of the protection strategy, OCRs must be operated with minimum time and be able to coordinate with another relay to improve the system reliability.Coordination of the OCR has been proposed using a number of techniques. These techniques may be divided into three categories: trial and error, topological analysis method, and optimization method [3]. However, in the solutions found by the first two classes, optimum coordination of relays was done manually by setting the plug multiplier setting (PMS) or pickup setting (PS) and the time multiplier setting (TMS) with expert knowledge or calculating analytically. Nowadays, the EPDN has become so complex that the conventional method of relay coordination is not a recommended choice. Therefore, optimal coordination of relays is essential to adapt [4]. The main objective of optimization is to minimize undesirable things or maximize desirable things in its mathematical model, subject to some constraints. Optimization is a commonly encountered mathematical problem in all engineering disciplines. Optimization problems are wide ranging and numerous, hence the various optimization techniques solve these problems. Some of these optimization strategies may be found in [5–10].Numerous methods and techniques have been published to solve the optimal setting parameter of OCR in order to get better results. In some research studies, the PMS has been predetermined and curve setting (CS) of each relay have been set to specific characteristic curve, hence this problem has been formulated as a linear optimization. The minimal overall operational time of the primary relays is the objective function (OF) subjected to the coordination constraint. As a consequence, the output for the optimal result is the TMS of each relay. The linear programming optimization (LPO) [11], the dual simplex and genetic algorithm (DS&GA) approaches [12], mixed integer linear programming (MILP) [13], Big-M (BM) strategy and the reverse simplex method by [14], and artificial bee colony [15] have all been suggested as solutions to these problems. While some other research refines the issue to a nonlinear optimization by using the TMS and PS as a decision variable, while the CS stayed unchanged. In order to perform this method, several artificial intelligence techniques were proposed. For example, adaptive modified fly fire algorithm (AMFA) was introduced in [16]. The gravitational search algorithm (GSA) and simulated annealing-based symbiotic organism search (SASOS) were proposed in [17, 18], respectively. Afterward, the researchers recommended using the PS as a step/discrete variable rather than a continuous variable to depict the actual industrial OCR. This issue was solved using a variety of techniques such as the modified seeker algorithm (MSA) [19], ant lion optimization (ALO) [20], and hybrid whale optimization algorithm (HWOA) [21]. Furthermore, some researchers concentrated on maximum PMS as a decision variable that can minimize the total CTI and total operating time of the relay [22].With the installation of distributed energy resources (DER) at medium and low voltage levels, the power flows and short-circuit currents can have different directions and values. Consequently, nondirectional OCRs alone are prone to loss of efficiency in such networks as unwanted tripping (maloperation) is likely to occur [23]. To remedy this issue, several “mitigation methods” have been proposed. An abstract about these methods and their drawbacks is given in [24].One of the mitigating strategies is to employ directional OCR (DOCR) and optimally set at both ends of every line. In fact, the need for DOCRs in distribution networks (DNs) has been raised due to the emergence of DER units because of the ability of energy to flow in two directions. Since conventional directional relaying schemes utilize voltage (line-side voltage or bus-side voltage) as a reference quantity to detect the fault direction (FD), so they are not suitable for DNs due to the absence of potential transformers. Therefore, the trend now is toward using the current-based directional relaying in DNs; this will be accomplished through this paper.Recently, many research studies have concentrated on the application of various current-dependent only schemes for FD detection, and these schemes were applied not only to relays in active distribution networks but also to directional relays. Among these research studies, Hosseini et al. [25] offered a scheme to directional detection possibility relies on the current, which can detect the fault direction by using the postfault current. To accomplish this, the fault direction was detected using a new reference phasor that was created from postfault current using the S-Transform. A technique based on the angle between the fault and prefault positive sequence current phasor was given in [26]. This method is used for faults on series-compensated high-voltage transmission lines when single-pole tripping conditions are present. The technique given in [27] is based on fault and prefault current samples. The summation and direction of power flow under normal conditions establishes the fault direction. Based on an investigation of prefault and fault phase angles, the method described in [23] provides a fault directional identification methodology that depends on the positive sequence components of the current.In this study, the newly proposed scheme also depends on the current only to determine the FD, but by using the PSCC angle sign of fault current without resorting to comparison with the pre/postfault current, which increases the performance of the proposed scheme’s speed and accuracy.Finally, the main contributions of this study can be summarized as follows:(i)
A neoteric optimal coordination algorithm for OCRs is implemented. The TMS and PS of all relays are considered optimization parameters in order to reduce the operating time. The performance of the developed technique is compared with other recently techniques to discover its superiority for solving the problem described.(ii)
Enhancing the performance of OCR by adding a new circuit to it: this circuit uses the sign angle of positive sequence current component (PSCC) to detect the FD. The modified circuit is proposed since it has been noticed that when the direction of the current is changed, the PSCC angle sign changes.(iii)
The new suggested circuit uses only a faulted current signal from the current transformer (CT), which is the main signal to the OCR and where no other measuring devices should be added; the resulting added cost is very small compared to the other used solutions. In addition, this circuit is simple to be connected with the OCR.(iv)
The optimal coordination and the proposed circuit are cooperating to assist the relay in making the correct decision in fault detection and clearing in the right form with high penetration of DERs. The proposed OCRs enhancement techniques are applied on a medium scale radial network configuration (the IEEE 33-node test feeder system) modeled in MATLAB/Coding and Simulink.


## 2. OCR Performance

OCR performance is described as follows.


### 2.1. Relay Characteristics

According to IEC and IEEE standards and in practical terms, the DN protection schemes consider one type of characteristic curve only. This characteristic determines the relation between the input current and operation time, which can be represented using a characteristic curve. The characteristics of such curve are described in [28] as follows:(1)TOP=β∗TMSiIiF/CTRi∗PSiα−1,where TOP is the operating time of relay, β and α are the relay characteristic constants selected from the values given in the IEC and IEEE standards, IiF is the fault current measured at primary terminal of CT, CTRi is the ith relay CT ratio, and PSi is the ith relay pickup current.Considering the IEC standard, inverse definite minimum time relay (IDMT) is used as OCR because of its best characteristic [4], especially in DNs, so it was selected here (β=0.14 and α=0.02). A typical IDMT has two units: a definite time unit and an inverse overcurrent unit (time dependent). The definite unit operates with constant intentional time delay when the current is above a predefined threshold value. However, for current that is less than the definite current setting but exceeds the maximum load current owing to a fault, a time-delay unit is employed. This unit is activated when a fault occurs with an intended time delay [29]. The time-delay unit is associated with two settings, namely, PS or PMS and TMS.


### 2.2. OCR Coordination in a Radial DN

We consider a radial DN as in Figure1, containing three buses and without DER penetration. If a fault is initiated at point k near bus B3, the relay R3 closest to the fault section must operate to isolate this section. If the relay R3 does not operate for any reason, R2 will work as a backup relay, where the setting of relay R2 is higher than the relay R3. In general, the setting of relay R2 is the addition of the operating time of R3, the circuit breaker (CB) operating time at a bus, and the overshoot time of relay R2. Similarly, R1 will work as a backup for R2 and R3. This implies that, in order to plan a protection order, the relay closest to the fault (primary) should be the quickest to act to isolate it, and what follows toward the upstream is its backup.Figure 1
Radial DN with 3 buses.


### 2.3. Impact of DERs on Protection System Using OCR

An increased DER penetration leads to the change of fault current magnitudes and their distribution in the network, thereby jeopardizing the existing protection system. The extent of DERs impact depends greatly on type, size, and location of the DERs. For instance, the fault current contribution from the rotating machines-based DER is significantly large and sustained compared to a small fault current contribution from inverter interfaced DER. Moreover, the impact rises with increasing DERs size. Therefore, the major protection challenges can be summarized as follows.


#### 2.3.1. Protection Blinding

It is a case in which the fault current seen by a relay gets changed due to DER connections, thereby leading to maloperation of protective devices.


#### 2.3.2. False Tripping

It is a case that occurs when DER causes a relay connected to a feeder to operate for a fault in an adjacent feeder.


#### 2.3.3. Intentional/Unintentional Islanding

Islanding is basically a condition where DER continues to supply power for a particular area during the loss of supply from the utility. Typically, the short circuit current during islanded operation is very low compared to the grid-connected mode. This might lead to a condition where relays having the setting for the grid-connected mode might not discriminate fault during islanded operations and vice versa. For more details and examples, we review [30].


## 2.1. Relay Characteristics

According to IEC and IEEE standards and in practical terms, the DN protection schemes consider one type of characteristic curve only. This characteristic determines the relation between the input current and operation time, which can be represented using a characteristic curve. The characteristics of such curve are described in [28] as follows:(1)TOP=β∗TMSiIiF/CTRi∗PSiα−1,where TOP is the operating time of relay, β and α are the relay characteristic constants selected from the values given in the IEC and IEEE standards, IiF is the fault current measured at primary terminal of CT, CTRi is the ith relay CT ratio, and PSi is the ith relay pickup current.Considering the IEC standard, inverse definite minimum time relay (IDMT) is used as OCR because of its best characteristic [4], especially in DNs, so it was selected here (β=0.14 and α=0.02). A typical IDMT has two units: a definite time unit and an inverse overcurrent unit (time dependent). The definite unit operates with constant intentional time delay when the current is above a predefined threshold value. However, for current that is less than the definite current setting but exceeds the maximum load current owing to a fault, a time-delay unit is employed. This unit is activated when a fault occurs with an intended time delay [29]. The time-delay unit is associated with two settings, namely, PS or PMS and TMS.


## 2.2. OCR Coordination in a Radial DN

We consider a radial DN as in Figure1, containing three buses and without DER penetration. If a fault is initiated at point k near bus B3, the relay R3 closest to the fault section must operate to isolate this section. If the relay R3 does not operate for any reason, R2 will work as a backup relay, where the setting of relay R2 is higher than the relay R3. In general, the setting of relay R2 is the addition of the operating time of R3, the circuit breaker (CB) operating time at a bus, and the overshoot time of relay R2. Similarly, R1 will work as a backup for R2 and R3. This implies that, in order to plan a protection order, the relay closest to the fault (primary) should be the quickest to act to isolate it, and what follows toward the upstream is its backup.Figure 1
Radial DN with 3 buses.


## 2.3. Impact of DERs on Protection System Using OCR

An increased DER penetration leads to the change of fault current magnitudes and their distribution in the network, thereby jeopardizing the existing protection system. The extent of DERs impact depends greatly on type, size, and location of the DERs. For instance, the fault current contribution from the rotating machines-based DER is significantly large and sustained compared to a small fault current contribution from inverter interfaced DER. Moreover, the impact rises with increasing DERs size. Therefore, the major protection challenges can be summarized as follows.


### 2.3.1. Protection Blinding

It is a case in which the fault current seen by a relay gets changed due to DER connections, thereby leading to maloperation of protective devices.


### 2.3.2. False Tripping

It is a case that occurs when DER causes a relay connected to a feeder to operate for a fault in an adjacent feeder.


### 2.3.3. Intentional/Unintentional Islanding

Islanding is basically a condition where DER continues to supply power for a particular area during the loss of supply from the utility. Typically, the short circuit current during islanded operation is very low compared to the grid-connected mode. This might lead to a condition where relays having the setting for the grid-connected mode might not discriminate fault during islanded operations and vice versa. For more details and examples, we review [30].


## 2.3.1. Protection Blinding

It is a case in which the fault current seen by a relay gets changed due to DER connections, thereby leading to maloperation of protective devices.


## 2.3.2. False Tripping

It is a case that occurs when DER causes a relay connected to a feeder to operate for a fault in an adjacent feeder.


## 2.3.3. Intentional/Unintentional Islanding

Islanding is basically a condition where DER continues to supply power for a particular area during the loss of supply from the utility. Typically, the short circuit current during islanded operation is very low compared to the grid-connected mode. This might lead to a condition where relays having the setting for the grid-connected mode might not discriminate fault during islanded operations and vice versa. For more details and examples, we review [30].


## 3. Problem Formulation for OCR Coordination Optimization Problem

In this study, the coordination among OCRs of a radial DN is considered as an optimization problem. The objective of the optimization problem is to minimize the total operating time of all relays.


### 3.1. Constraints of Relay Coordination

#### 3.1.1. Operating Time Constraints

As much as minimum operation time is very much desired for a relay, there are some limitations within the system that needs to be considered while coordinating OCRs. A relay needs to have a minimum and maximum operation time so that the protection engineer can set up the values within that range. The relay operation time can be mathematically represented as follows:(2)Timin≤Tik≥Timax,where Timin and Timax are lower and upper boundary of ith relay operating time. Typical values for Tmin = 0.1 s and Tmax = 1.1 s [31]. Whereas, Tmin is determined by relay manufacture, and Tmax is based on critical clearing time.


#### 3.1.2. Coordination Time Interval Constraints

The backup relay must operate only when primary relay fails to work in its protective zone under fault condition. The backup relay operating time must satisfyCTI constraints. The CTI value is mostly dictated by the practical limitations, which consist of the relay over travel time, the CB operating time, and the safety margin for the relay error [32], and can be defined as follows:(3)Tbk−Tpk≥CTI,where Tpk is operating time of primary relay for fault at k and Tbk is the operating time of backup relay for the same fault at k. The discrimination time between the backup and main OCRs, Δtbp, is obtained by(4)Λtbp=Tbk−Tpk−CTI.


#### 3.1.3. The TMS and PS Constraints for Each Relay (Control Variables and Constraints)

The OCR has two settings which are important in determining the relay operating time, namely,TMS and PS, which should be bounded in a range such that fast operation can be achieved with ignoring false operation due to switching external load. They are expressed as follows:(5)TMSimin≤TMSi≤TMSimax,PSimin≤PSi≤PSimax,where TMSimin is the minimum value of TMS, TMSmax is the maximum value of TMS, PSmin is the minimum value of PS, and PSmax is the maximum value of PS.


#### 3.1.4. Objective Function for Overcurrent Relay Coordination Optimization Problem

The OCRs have two control variables: TMS, which sets the relay operating time, and PS, which sets the operational current value. The PS value is calculated based on the maximum load current and minimum fault current. For the optimum relay coordination problem, the minimization of overall operating time is considered as the OF. The most common OF that models the coordination problem is described as follows [33]:(6)OF=min∑i=1NR∑kTik,where i is the number of relays in the network and Tik is the operation time of the ith relay for a fault at k.In this study, the proposed modified objective function (MOF) in [34] is implemented,but with the replacement of the absolute value with the squared value, since the results are observed to be balanced and improved, with this substitution. This MOF improves CTI and overcomes the extremely delayed of backup relay. The MOF can be expressed as follows:(7)MOF=∝1∗∑i=1NTopiM+∝2∗∑Δtbp−β∗Δtbp−Δtbp2N,where ∝1, ∝2, and β are weighting coefficients, M and N are constant values, and Topi is the operating time of ith OCR when a fault occurs near the relay.


## 3.1. Constraints of Relay Coordination

### 3.1.1. Operating Time Constraints

As much as minimum operation time is very much desired for a relay, there are some limitations within the system that needs to be considered while coordinating OCRs. A relay needs to have a minimum and maximum operation time so that the protection engineer can set up the values within that range. The relay operation time can be mathematically represented as follows:(2)Timin≤Tik≥Timax,where Timin and Timax are lower and upper boundary of ith relay operating time. Typical values for Tmin = 0.1 s and Tmax = 1.1 s [31]. Whereas, Tmin is determined by relay manufacture, and Tmax is based on critical clearing time.


### 3.1.2. Coordination Time Interval Constraints

The backup relay must operate only when primary relay fails to work in its protective zone under fault condition. The backup relay operating time must satisfyCTI constraints. The CTI value is mostly dictated by the practical limitations, which consist of the relay over travel time, the CB operating time, and the safety margin for the relay error [32], and can be defined as follows:(3)Tbk−Tpk≥CTI,where Tpk is operating time of primary relay for fault at k and Tbk is the operating time of backup relay for the same fault at k. The discrimination time between the backup and main OCRs, Δtbp, is obtained by(4)Λtbp=Tbk−Tpk−CTI.


### 3.1.3. The TMS and PS Constraints for Each Relay (Control Variables and Constraints)

The OCR has two settings which are important in determining the relay operating time, namely,TMS and PS, which should be bounded in a range such that fast operation can be achieved with ignoring false operation due to switching external load. They are expressed as follows:(5)TMSimin≤TMSi≤TMSimax,PSimin≤PSi≤PSimax,where TMSimin is the minimum value of TMS, TMSmax is the maximum value of TMS, PSmin is the minimum value of PS, and PSmax is the maximum value of PS.


### 3.1.4. Objective Function for Overcurrent Relay Coordination Optimization Problem

The OCRs have two control variables: TMS, which sets the relay operating time, and PS, which sets the operational current value. The PS value is calculated based on the maximum load current and minimum fault current. For the optimum relay coordination problem, the minimization of overall operating time is considered as the OF. The most common OF that models the coordination problem is described as follows [33]:(6)OF=min∑i=1NR∑kTik,where i is the number of relays in the network and Tik is the operation time of the ith relay for a fault at k.In this study, the proposed modified objective function (MOF) in [34] is implemented,but with the replacement of the absolute value with the squared value, since the results are observed to be balanced and improved, with this substitution. This MOF improves CTI and overcomes the extremely delayed of backup relay. The MOF can be expressed as follows:(7)MOF=∝1∗∑i=1NTopiM+∝2∗∑Δtbp−β∗Δtbp−Δtbp2N,where ∝1, ∝2, and β are weighting coefficients, M and N are constant values, and Topi is the operating time of ith OCR when a fault occurs near the relay.


## 3.1.1. Operating Time Constraints

As much as minimum operation time is very much desired for a relay, there are some limitations within the system that needs to be considered while coordinating OCRs. A relay needs to have a minimum and maximum operation time so that the protection engineer can set up the values within that range. The relay operation time can be mathematically represented as follows:(2)Timin≤Tik≥Timax,where Timin and Timax are lower and upper boundary of ith relay operating time. Typical values for Tmin = 0.1 s and Tmax = 1.1 s [31]. Whereas, Tmin is determined by relay manufacture, and Tmax is based on critical clearing time.


## 3.1.2. Coordination Time Interval Constraints

The backup relay must operate only when primary relay fails to work in its protective zone under fault condition. The backup relay operating time must satisfyCTI constraints. The CTI value is mostly dictated by the practical limitations, which consist of the relay over travel time, the CB operating time, and the safety margin for the relay error [32], and can be defined as follows:(3)Tbk−Tpk≥CTI,where Tpk is operating time of primary relay for fault at k and Tbk is the operating time of backup relay for the same fault at k. The discrimination time between the backup and main OCRs, Δtbp, is obtained by(4)Λtbp=Tbk−Tpk−CTI.


## 3.1.3. The TMS and PS Constraints for Each Relay (Control Variables and Constraints)

The OCR has two settings which are important in determining the relay operating time, namely,TMS and PS, which should be bounded in a range such that fast operation can be achieved with ignoring false operation due to switching external load. They are expressed as follows:(5)TMSimin≤TMSi≤TMSimax,PSimin≤PSi≤PSimax,where TMSimin is the minimum value of TMS, TMSmax is the maximum value of TMS, PSmin is the minimum value of PS, and PSmax is the maximum value of PS.


## 3.1.4. Objective Function for Overcurrent Relay Coordination Optimization Problem

The OCRs have two control variables: TMS, which sets the relay operating time, and PS, which sets the operational current value. The PS value is calculated based on the maximum load current and minimum fault current. For the optimum relay coordination problem, the minimization of overall operating time is considered as the OF. The most common OF that models the coordination problem is described as follows [33]:(6)OF=min∑i=1NR∑kTik,where i is the number of relays in the network and Tik is the operation time of the ith relay for a fault at k.In this study, the proposed modified objective function (MOF) in [34] is implemented,but with the replacement of the absolute value with the squared value, since the results are observed to be balanced and improved, with this substitution. This MOF improves CTI and overcomes the extremely delayed of backup relay. The MOF can be expressed as follows:(7)MOF=∝1∗∑i=1NTopiM+∝2∗∑Δtbp−β∗Δtbp−Δtbp2N,where ∝1, ∝2, and β are weighting coefficients, M and N are constant values, and Topi is the operating time of ith OCR when a fault occurs near the relay.


## 4. Proposed Improved Wild Horse Optimization (IWHO) Algorithm Mathematical Model

Before approaching the IWHO, wild horse optimizer (WHO) should be studied.


### 4.1. Wild Horse Optimizer

Nonterrestrial horses are herds made up of stable family groups or harems consisting of a stallion, one or more mares, and offspring. There are also other groups, such as mature stallions and adolescent horses. To know a lot of the social life behavior of these herds, we refer to [35]. The WHO is a novel algorithm proposed to mimic the social life behavior of these wild horses (WHs) in the wild [35]. The WH engages in a variety of behaviors, such as mating, grazing, chasing, dominating, and commanding. Accordingly, five steps of the WHO algorithm can be considered.


#### 4.1.1. Creating an Initial Population and Forming Horse Groups and Selecting Leaders

First, the initial population is separated into many groups. If N is the members number of the population and PC is a control parameter for the proposed algorithm that represents the percentage of stallions in the entire population, then the groups number isG = ⌈N × PC⌉. Each group has a leader (stallion), so the number of stallions in the algorithm is G. Foals and mares, which make up the remaining population (N-G), are dispersed evenly among these groups. Figure 2 shows how the foals and stallions have been chosen from the initial population to create the different groups.Figure 2
Groups’ formation from the original population.


#### 4.1.2. Grazing Behavior

To model how group members, move, and search (razing behavior) around the leader with a variable radius, equation (8) is proposed(8)X¯i,Gj=2∗Zcos2πRZ∗Stallionj−Xi,Gj+Stallionj,where Xi,Gj is the present position of the group member (foal or mare), Stallionj is the stallion position, and R is a uniform random number in the range [−2, 2] that causes the horses grazing at different angles (360 degrees) of group leader. The cos function by combining π and R causes the movement in different radius, X¯i,Gj is the new position of the group member when grazing, and Z is the adaptive mechanism calculated from the following equation:(9)AP=R1⟶<TDR…IDX=AP==0,Z=R2ΘIDX+R3⟶Θ∼IDX,where AP is a vector consisting of 0 to 1, which is equal to dimensions of the problem, R1⟶ and R3⟶ are random vectors with uniform distribution in the range [0, 1], R2 is a random number with uniform distribution in the range [0, 1], and IDX is indexes of the random vector R1⟶ returns that satisfy the condition (AP == 0). TDR is an adaptive parameter that begins with 1 and declines until it reaches 0 at the end of the algorithm implementation according to the following equation:(10)TDR=1−IT∗1ITmax,where IT is the current iteration and ITmax is the maximum number of algorithm iterations.


#### 4.1.3. Horse Mating Behavior

One of the distinctive behaviors of horses is that, before puberty, the foals leave the group, with female foals joining another family group and male foals joining a group of solitary horses. By leaving, the father will not be able to mate with the daughter or siblings. In order to implement this behavior, we suppose that there are 3 groups, that is,i, j, and h. A foal departs group i and joins a temporary group; a foal departs group j and joins a tentative group. We presume these two foals are female and male, and since these two foals have no familial connection, they can mate once they reach puberty. The resulting child must depart the tentative group and join another group, such as h. To simulate the horses mating behavior, the Crossover operator of the mean type is suggested as follows:(11)XG,hp=CrossoverXG,iq,XG,jzi≠j≠h,p=q=end,Crossover=Mean,where XG,hp is the position of horse p from group h that departs the group and is replaced by a horse whose parents are horses that must leave groups i (horse q) and j (horse z) and have reached puberty. They have no familial ties and have mated and reproduced.


#### 4.1.4. Group Leadership

In the WHO algorithm, if the stallions (group leaders) lead the group to a watering hole, for example, leaders struggle over this water hole so that the dominance group may use it first, and other groups are barred from using it until the domination group leaves. Equation (12) is recommended to do this approach and distance:(12)stallionGi=2∗Zcos2πRZ∗WH−stallionGi+WH,ifR3>0.52∗Zcos2πRZ∗WH−stallionGi−WH,ifR3≤0.5,where stallionGi is the leader’s current position of the i group and stallionGi is the next position of the leader. WH is the water hole location.


#### 4.1.5. Leaders Selection and Interchange

According the nature of the algorithm, the leaders are first selected at random and then chosen depending on fitness in the latter algorithm stages. If one of the group members is best fitness than the current group leader, the leader position and the relevant member will alter based on the following equation :(13)stallionGi=XG,iifcostXG,i<coststallionGistallionGiifcostXG,i>coststallionGi.For the full flowchart of WHO, see [35].


### 4.2. Improved Wild Horse Optimizer

IWHO relies on the cuckoo search algorithm as described in [36]. During the iteration of the proposed algorithm, the new solution is generated using the “Levy flight” [30] as provided in the following equation:(14)Xi,G=Xi,G−γXi,G−Xg⊕lavy′λ=Xi,G+0.01uv1/λXi,G−Xg,where Xi,G is the ith position of the group member, γ denotes the step scaling size, Xg denotes the global best solution, ⊕ refers to the procedure of element-wise multiplications, and λ refers to the Levy flight exponent, while u and v are defined as follows:(15)u∼N0,σu2,v∼N0,σv2.The standard deviationsσu and σv are represented as(16)σu=sinπλ/2.Γ1+λ2λ−1.λΓ1+λ/21/λ,σv=1,where Γ is the Gamma function; the new candidate solution is generated, and equation (10) is applied. The main advantage of this improvement is the ability to strike a balance local exploitation with global exploration and rapid convergence. Figure 3 shows the flowchart of the IWHO algorithm.Figure 3
Flowchart of IWHO algorithm.


## 4.1. Wild Horse Optimizer

Nonterrestrial horses are herds made up of stable family groups or harems consisting of a stallion, one or more mares, and offspring. There are also other groups, such as mature stallions and adolescent horses. To know a lot of the social life behavior of these herds, we refer to [35]. The WHO is a novel algorithm proposed to mimic the social life behavior of these wild horses (WHs) in the wild [35]. The WH engages in a variety of behaviors, such as mating, grazing, chasing, dominating, and commanding. Accordingly, five steps of the WHO algorithm can be considered.


### 4.1.1. Creating an Initial Population and Forming Horse Groups and Selecting Leaders

First, the initial population is separated into many groups. If N is the members number of the population and PC is a control parameter for the proposed algorithm that represents the percentage of stallions in the entire population, then the groups number isG = ⌈N × PC⌉. Each group has a leader (stallion), so the number of stallions in the algorithm is G. Foals and mares, which make up the remaining population (N-G), are dispersed evenly among these groups. Figure 2 shows how the foals and stallions have been chosen from the initial population to create the different groups.Figure 2
Groups’ formation from the original population.


### 4.1.2. Grazing Behavior

To model how group members, move, and search (razing behavior) around the leader with a variable radius, equation (8) is proposed(8)X¯i,Gj=2∗Zcos2πRZ∗Stallionj−Xi,Gj+Stallionj,where Xi,Gj is the present position of the group member (foal or mare), Stallionj is the stallion position, and R is a uniform random number in the range [−2, 2] that causes the horses grazing at different angles (360 degrees) of group leader. The cos function by combining π and R causes the movement in different radius, X¯i,Gj is the new position of the group member when grazing, and Z is the adaptive mechanism calculated from the following equation:(9)AP=R1⟶<TDR…IDX=AP==0,Z=R2ΘIDX+R3⟶Θ∼IDX,where AP is a vector consisting of 0 to 1, which is equal to dimensions of the problem, R1⟶ and R3⟶ are random vectors with uniform distribution in the range [0, 1], R2 is a random number with uniform distribution in the range [0, 1], and IDX is indexes of the random vector R1⟶ returns that satisfy the condition (AP == 0). TDR is an adaptive parameter that begins with 1 and declines until it reaches 0 at the end of the algorithm implementation according to the following equation:(10)TDR=1−IT∗1ITmax,where IT is the current iteration and ITmax is the maximum number of algorithm iterations.


### 4.1.3. Horse Mating Behavior

One of the distinctive behaviors of horses is that, before puberty, the foals leave the group, with female foals joining another family group and male foals joining a group of solitary horses. By leaving, the father will not be able to mate with the daughter or siblings. In order to implement this behavior, we suppose that there are 3 groups, that is,i, j, and h. A foal departs group i and joins a temporary group; a foal departs group j and joins a tentative group. We presume these two foals are female and male, and since these two foals have no familial connection, they can mate once they reach puberty. The resulting child must depart the tentative group and join another group, such as h. To simulate the horses mating behavior, the Crossover operator of the mean type is suggested as follows:(11)XG,hp=CrossoverXG,iq,XG,jzi≠j≠h,p=q=end,Crossover=Mean,where XG,hp is the position of horse p from group h that departs the group and is replaced by a horse whose parents are horses that must leave groups i (horse q) and j (horse z) and have reached puberty. They have no familial ties and have mated and reproduced.


### 4.1.4. Group Leadership

In the WHO algorithm, if the stallions (group leaders) lead the group to a watering hole, for example, leaders struggle over this water hole so that the dominance group may use it first, and other groups are barred from using it until the domination group leaves. Equation (12) is recommended to do this approach and distance:(12)stallionGi=2∗Zcos2πRZ∗WH−stallionGi+WH,ifR3>0.52∗Zcos2πRZ∗WH−stallionGi−WH,ifR3≤0.5,where stallionGi is the leader’s current position of the i group and stallionGi is the next position of the leader. WH is the water hole location.


### 4.1.5. Leaders Selection and Interchange

According the nature of the algorithm, the leaders are first selected at random and then chosen depending on fitness in the latter algorithm stages. If one of the group members is best fitness than the current group leader, the leader position and the relevant member will alter based on the following equation :(13)stallionGi=XG,iifcostXG,i<coststallionGistallionGiifcostXG,i>coststallionGi.For the full flowchart of WHO, see [35].


## 4.1.1. Creating an Initial Population and Forming Horse Groups and Selecting Leaders

First, the initial population is separated into many groups. If N is the members number of the population and PC is a control parameter for the proposed algorithm that represents the percentage of stallions in the entire population, then the groups number isG = ⌈N × PC⌉. Each group has a leader (stallion), so the number of stallions in the algorithm is G. Foals and mares, which make up the remaining population (N-G), are dispersed evenly among these groups. Figure 2 shows how the foals and stallions have been chosen from the initial population to create the different groups.Figure 2
Groups’ formation from the original population.


## 4.1.2. Grazing Behavior

To model how group members, move, and search (razing behavior) around the leader with a variable radius, equation (8) is proposed(8)X¯i,Gj=2∗Zcos2πRZ∗Stallionj−Xi,Gj+Stallionj,where Xi,Gj is the present position of the group member (foal or mare), Stallionj is the stallion position, and R is a uniform random number in the range [−2, 2] that causes the horses grazing at different angles (360 degrees) of group leader. The cos function by combining π and R causes the movement in different radius, X¯i,Gj is the new position of the group member when grazing, and Z is the adaptive mechanism calculated from the following equation:(9)AP=R1⟶<TDR…IDX=AP==0,Z=R2ΘIDX+R3⟶Θ∼IDX,where AP is a vector consisting of 0 to 1, which is equal to dimensions of the problem, R1⟶ and R3⟶ are random vectors with uniform distribution in the range [0, 1], R2 is a random number with uniform distribution in the range [0, 1], and IDX is indexes of the random vector R1⟶ returns that satisfy the condition (AP == 0). TDR is an adaptive parameter that begins with 1 and declines until it reaches 0 at the end of the algorithm implementation according to the following equation:(10)TDR=1−IT∗1ITmax,where IT is the current iteration and ITmax is the maximum number of algorithm iterations.


## 4.1.3. Horse Mating Behavior

One of the distinctive behaviors of horses is that, before puberty, the foals leave the group, with female foals joining another family group and male foals joining a group of solitary horses. By leaving, the father will not be able to mate with the daughter or siblings. In order to implement this behavior, we suppose that there are 3 groups, that is,i, j, and h. A foal departs group i and joins a temporary group; a foal departs group j and joins a tentative group. We presume these two foals are female and male, and since these two foals have no familial connection, they can mate once they reach puberty. The resulting child must depart the tentative group and join another group, such as h. To simulate the horses mating behavior, the Crossover operator of the mean type is suggested as follows:(11)XG,hp=CrossoverXG,iq,XG,jzi≠j≠h,p=q=end,Crossover=Mean,where XG,hp is the position of horse p from group h that departs the group and is replaced by a horse whose parents are horses that must leave groups i (horse q) and j (horse z) and have reached puberty. They have no familial ties and have mated and reproduced.


## 4.1.4. Group Leadership

In the WHO algorithm, if the stallions (group leaders) lead the group to a watering hole, for example, leaders struggle over this water hole so that the dominance group may use it first, and other groups are barred from using it until the domination group leaves. Equation (12) is recommended to do this approach and distance:(12)stallionGi=2∗Zcos2πRZ∗WH−stallionGi+WH,ifR3>0.52∗Zcos2πRZ∗WH−stallionGi−WH,ifR3≤0.5,where stallionGi is the leader’s current position of the i group and stallionGi is the next position of the leader. WH is the water hole location.


## 4.1.5. Leaders Selection and Interchange

According the nature of the algorithm, the leaders are first selected at random and then chosen depending on fitness in the latter algorithm stages. If one of the group members is best fitness than the current group leader, the leader position and the relevant member will alter based on the following equation :(13)stallionGi=XG,iifcostXG,i<coststallionGistallionGiifcostXG,i>coststallionGi.For the full flowchart of WHO, see [35].


## 4.2. Improved Wild Horse Optimizer

IWHO relies on the cuckoo search algorithm as described in [36]. During the iteration of the proposed algorithm, the new solution is generated using the “Levy flight” [30] as provided in the following equation:(14)Xi,G=Xi,G−γXi,G−Xg⊕lavy′λ=Xi,G+0.01uv1/λXi,G−Xg,where Xi,G is the ith position of the group member, γ denotes the step scaling size, Xg denotes the global best solution, ⊕ refers to the procedure of element-wise multiplications, and λ refers to the Levy flight exponent, while u and v are defined as follows:(15)u∼N0,σu2,v∼N0,σv2.The standard deviationsσu and σv are represented as(16)σu=sinπλ/2.Γ1+λ2λ−1.λΓ1+λ/21/λ,σv=1,where Γ is the Gamma function; the new candidate solution is generated, and equation (10) is applied. The main advantage of this improvement is the ability to strike a balance local exploitation with global exploration and rapid convergence. Figure 3 shows the flowchart of the IWHO algorithm.Figure 3
Flowchart of IWHO algorithm.


## 5. Newly Proposed Protection Scheme

A newly proposed protection scheme is described in the following sections.


### 5.1. Description of the Scheme

As mentioned before, the outstanding feature of the proposed scheme is the detection of the FD using only the fault current. Therefore, a scheme has been added to OCRs to determine the direction of current flow. For this purpose, the fundamental frequency components of current (from CT) and its positive sequence component should be extracted. An advantage of the positive sequence component is its availability for all types of faults; in contrast, negative and zero sequences are present in the case of asymmetrical faults, which is one of the biggest obstacles to use some DOCRs classifiers in practical networks.The PSCC can be calculated as follows:(17)PSCC=13IA+a∗IB+a2∗IC,where IA,IB and IC are the current of phaseA, B,andC, respectively, and a is the unit phase shift operator and equal to 1∠120o.The following steps describe how to access and implement the newly proposed scheme.After setting up, the optimal relay setting (PS,TMS) resulted from the optimal coordination technique on studied system (a redial DN containing a large number of buses). During the simulation on this system, the PSCC angle (from polar form) was monitored at each relay in two cases:(a)
The network is without DER, and it has one current flow direction, which is from the substation to the loads (from upstream to downstream). It was found that the values of PSCC angle take different values, but keep the same positive sign.(b)
The DER in the network is merged, and the PSCC angle values were monitored. It has been found that these values take either a positive sign on all buses feeding from a previous direction (from upstream to downstream), which is the original direction, or a negative sign, which happened when the loads feed from the opposite direction to the original direction, i.e., they do not feed from the substation but from DERs in the opposite direction.Depending on this characteristic, a scheme based on the sign of this angle has been developed and incorporated for each OCR to detect the FD. So, if the PSCC angle is positive, the trip operating time depends on the relay settings that have been optimally tuned in advance for each relay. If the angle sign is negative, this means that the flow is in the opposite direction, so the relay does not trip until it has a signal from the reverse coordination or control circuit (automation system), and this will be the second part to supplement this research. The proposed protection scheme flowchart is shown in Figure4.Figure 4
Flowchart for the proposed protection scheme.


## 5.1. Description of the Scheme

As mentioned before, the outstanding feature of the proposed scheme is the detection of the FD using only the fault current. Therefore, a scheme has been added to OCRs to determine the direction of current flow. For this purpose, the fundamental frequency components of current (from CT) and its positive sequence component should be extracted. An advantage of the positive sequence component is its availability for all types of faults; in contrast, negative and zero sequences are present in the case of asymmetrical faults, which is one of the biggest obstacles to use some DOCRs classifiers in practical networks.The PSCC can be calculated as follows:(17)PSCC=13IA+a∗IB+a2∗IC,where IA,IB and IC are the current of phaseA, B,andC, respectively, and a is the unit phase shift operator and equal to 1∠120o.The following steps describe how to access and implement the newly proposed scheme.After setting up, the optimal relay setting (PS,TMS) resulted from the optimal coordination technique on studied system (a redial DN containing a large number of buses). During the simulation on this system, the PSCC angle (from polar form) was monitored at each relay in two cases:(a)
The network is without DER, and it has one current flow direction, which is from the substation to the loads (from upstream to downstream). It was found that the values of PSCC angle take different values, but keep the same positive sign.(b)
The DER in the network is merged, and the PSCC angle values were monitored. It has been found that these values take either a positive sign on all buses feeding from a previous direction (from upstream to downstream), which is the original direction, or a negative sign, which happened when the loads feed from the opposite direction to the original direction, i.e., they do not feed from the substation but from DERs in the opposite direction.Depending on this characteristic, a scheme based on the sign of this angle has been developed and incorporated for each OCR to detect the FD. So, if the PSCC angle is positive, the trip operating time depends on the relay settings that have been optimally tuned in advance for each relay. If the angle sign is negative, this means that the flow is in the opposite direction, so the relay does not trip until it has a signal from the reverse coordination or control circuit (automation system), and this will be the second part to supplement this research. The proposed protection scheme flowchart is shown in Figure4.Figure 4
Flowchart for the proposed protection scheme.


## 6. Simulation Results and Discussion

The proposed algorithm and scheme ware implemented and tested on IEEE 33 bus DN. The test system has 32 relays with IDMT characteristics. The studied system is implemented by MATLAB/Simulink environment. The base kV and MVA of the system are 12.66 kV and 100 MVA, respectively. The system’s active and reactive load power are 3.715 MW and 2.3 MVAR, respectively. Branch data and loads connected in the destination bus are given in [37]. The system under investigation contains two DERs that provide 60% of the total load, its capacity, and location as in [38]. The IEEE 33 bus DN integrated with DERs single line diagram is shown in Figure 5.Figure 5
IEEE 33 bus distribution network with DER.


### 6.1. Optimal Coordination

Firstly, the DERs, which are combined with the test system, are not energized. As a result, the power flow and the short circuit current flow are in the same direction. Current flows from the substation to the load point or fault location. Figure6 demonstrate the load current passing through each OCR in the operating mode.Figure 6
Load current at each bus.The measured three-phase fault current (symmetrical fault current) passing through each OCR at near end and far end for each primary relay is given in Table1. According to Table 1, the primary and backup relay can see approximately the same fault current magnitude for a one fault source because this operating condition has a radial topology. Also, the fault current level is increased when the fault point is near the source (substation) and decreases when fault point is far from the source.Table 1
Three-phase fault current at near end and far end for each primary relay.
Faulted sectionPrimary relayBackup relayFault currentNearFarS1R1—11201100S2R2R111001058S3R3R210581021S4R4R31021983S5R5R4983875S6R6R5875811S7R7R6811768S8R8R7768684S9R9R8684614S10R10R9614606S11R11R10606590S12R12R11590512S13R13R12512480S14R14R13480455S15R15R14455430S16R16R15430379S17R17R16379360S18R18R111001075S19R19R181075858S20R20R19858801S21R21R20801710S22R22R21021999S23R23R22999882S24R24R23882785S25R25R5875857S26R26R25857832S27R27R26832722S28R28R27722658S29R29R28658630S30R30R29630561S31R31R30561541S32R32R31541515To achieve the relay coordination optimization goal, the control variables (TMS,PS) of all relays are acquired in order to clear the near-end fault and far-end fault, which means there are 64 decision variables for coordination (i.e., (TMS1-TMS32) and (PS1-PS32)). The relays’ operating time range is considered between 0.04–1.1 sec. The range of TMS are limited from 0.05–1.1 sec., while the PS lies between 1.25–1.5 and the minimum value of CTI is considered 0.2 sec [39].To verify the uniqueness and efficacy of the IWHO technique, other metaheuristic techniques based on the same MOF have been used to calculate the coordination issue for the sake of comparison. These algorithms are WHO, a mayfly optimization algorithm (MOA), and GSA. The results of the optimal controlled variables and OF obtained using the IWHO technique in comparison with other techniques for the IEEE 33 bus test system are tabulated in Tables2. By analyzing Table 2, it is found that the minimum operating time for the system is improved to 19.57667 sec when using IWHO algorithm, so these results were selected to set up its control variables on the IDMT relays on MATLAB/Simulink. Figure 7 also depicts the novelty of the IWHO algorithm, which converges rapidly (around 2 iteration) before the other algorithms to attain the optimal fitness function.Table 2
Optimal setting variables’ results for the test system’s OCRs.
RelayMethodGSAAMWHOIWHOTMSPSTMSPSTMSPSTMSPSR11.11.4484481.11.51.0817421.3050541.0318411.447165R21.0999371.4864621.0999521.50.8831911.3273580.885031.271911R31.0716711.4091981.092481.50.8818121.250110.7804661.499995R41.0674671.2825651.0664561.3909990.8460711.3479170.7364141.259884R51.0341371.4565951.0662021.4987460.833961.2885970.7130741.317918R61.0129931.4584471.0449491.50.7219011.4996760.5506331.499361R70.9885111.3272381.0315641.50.6673331.250.4536211.457282R80.8551031.4024450.981041.431290.5345321.3202220.3398961.27922R90.713011.3918720.9069791.4163060.4857871.4980190.3023151.250533R100.5736421.4844570.8228261.3009260.4495191.3058250.2979991.250022R110.4577611.4431170.6991211.444370.3165691.2939690.2819771.252573R120.3541871.4232150.5967871.250.2397731.4518120.2364781.289341R130.276941.3885240.4747571.4101560.2135141.253950.2304031.251034R140.2018021.4133530.3871941.2817020.1746921.458590.2148611.262904R150.136071.2806950.2822281.50.1629541.2721850.1866811.250296R160.0776821.2835760.1984341.250.0883181.2508720.1583321.322459R170.0500041.3566870.051.3048790.0510021.2506180.0504471.290968R181.0999811.4287941.11.51.0987541.2500381.0995241.498037R190.7618661.3447030.7531261.4526970.7047081.4263690.730441.313271R200.4332431.383540.4133231.4043030.3802381.2629670.36821.45484R210.0504461.326130.0500031.2502190.0502571.2514130.0514951.295263R221.0487731.3026531.0828141.4461460.7726771.4212780.7720551.499657R230.5490781.4225140.6227681.4264660.4040251.3008980.4212991.250451R240.1046521.4420780.1895161.3262170.0501631.250.0514171.250342R251.0826191.3055541.11.4985170.6588431.3118740.5687551.489642R260.8524621.4155920.965161.3962690.4521841.499930.4008491.481262R270.7035531.3137790.8320991.50.4018791.2587370.2407881.368357R280.5396941.4349310.7018831.3996960.3108731.2530390.2261641.342916R290.4027291.4073460.554111.410890.2285181.2620810.1732031.499112R300.2762631.3431290.3922781.2506060.1600471.250010.1253881.252482R310.1819671.4024740.2559331.2505280.125771.2505340.0846381.271112R320.051.3625530.051.3490390.0500961.2517560.0500361.476527Total operating time (sec)30.4523135.1855722.1799919.57667Figure 7
Convergence characteristic of proposed algorithms.


### 6.2. The Proposed Protection Scheme

After completing the optimization work to achieve optimal coordination on MATLAB/coding and depending on the best algorithm (IWHO), its results (TMS, PS) were taken and set up on IDMT relays on the studied system. During the simulation on this system, the PSCC angle was monitored at each relay for a number of scenarios that fall under two cases.


#### 6.2.1. Case 1: Without DERs’ Penetration

(1) At Normal Operation. In this scenario, the network has one current flow direction, which is from upstream to downstream. By analyzing Figure 8, it is found that the values of PSCC angle take different values, but keep the same positive sign, and this will be the base and considered as forward direction.Figure 8
PCSS angle values at normal operation (case 1).(2) At Abnormal Operation. In this scenario, the system is exposed to different fault types, sections, and fault resistance (Rf). In all fault types, the PSCC angle takes the same phenomena as the previous result, i.e., PSCC angle has a different value (affected by fault section) and a positive sign, which indicates that the direction of the current flow is also forward.Figures9–13 show the PSCC angle values and its sign at single line to ground (LG) fault at S3 with Rf  = 0.02 Ω, LG fault at S16 with Rf = 1000 Ω, double line to ground (LLG) fault at S20 with Rf = 0.02 Ω, three line to ground (LLLG) fault at S24 with Rf = 0.02 Ω, and line to line (LL) fault at S29 with Rf = 0.02 Ω, respectively. Note that the PSCC value of the buses after the fault, i.e., which is on the downstream side, is almost zero, but the PSCC angle value is recorded and plotted as observed from the simulation.Figure 9
PCSS angle values when LG fault occurs on S3 (case 1).Figure 10
PCSS angle values when LG fault occurs on S16 (case 1).Figure 11
PCSS angle values when LLG fault occurs on S20 (case 1).Figure 12
PCSS angle values when LLLG fault occurs on S24 (case 1).Figure 13
PCSS angle values when LL fault occurs on S29 (case 1).


#### 6.2.2. Case 2: With DERs’ Penetration

The system under investigation contains two DERs. DERs allocation (place and capacity) is based on reducing losses using optimization techniques. The DERs places, capacities, and the reduction of losses are given in Table3 as stipulated in [38]. It is clear from Table 3 that the DERs will provide 60% of the total load.Table 3
DERs allocation and power losses reduction.
DER numberDER locationDER sizePower loss reduction (%)MWMVAR1130.84650.396186.42301.13921.0611(1) At Normal Operation. In this scenario, the DERs are integrated with the network, so the network does not have a single current flow direction as it was in the first case. By analyzing Table 4, it was found that the PSCC angle values take different values, and at some relays, it has a positive sign and on the other has a negative sign. The positive sign means that the current flow from upstream to downstream, similar to the first case is the forward direction. As for the negative sign, it means that the direction of the current passage is reversed due to the presence of the DERs.Table 4
PCSS angle for all relays at normal operation (case 2).
Relay numberPSCC angleFlow directionR1124.1610ForwardR2124.1679ForwardR3123.0042ForwardR4115.7648ForwardR5114.2543ForwardR6111.6414ForwardR798.21968ForwardR824.69574ForwardR9−29.6059ReverseR10−33.1408ReverseR11−36.4600ReverseR12−39.3560ReverseR13−41.4905ReverseR14131.8906ForwardR15130.2809ForwardR16129.6603ForwardR17128.2808ForwardR18128.1673ForwardR19128.1388ForwardR20128.1286ForwardR21128.1181ForwardR22127.4158ForwardR23127.1179ForwardR24127.0909ForwardR25122.6414ForwardR26120.1087ForwardR2782.25281ForwardR28−34.9689ReverseR29−49.7763ReverseR30−60.4110ReverseR31126.0625ForwardR32122.3451ForwardTo clarify more, the PSCC angle for relays R9, R10, R11, R12, and R13 has a negative sign, which means that the loads at buses 9, 10, 11, 12, and 13, which these relays are located at, are fed by the DER1 that is located at bus 13 and not by the main source (substation). Of course, it is certain that the loads at buses 14, 15, 16, 17, and 18 are fed by this DER1 as well, but the direction of the current is the same as the original direction (forward direction), so its PSCC angle will take the positive sign as it was. By analogy with that, the loads at buses 28, 29, and 30 are fed from the DER2 that is located at bus 30, so they take the negative sign; as well as, loads at buses 31 and 32 are fed from this DER2, but the direction of the current is the same as the original direction, so it takes the positive sign.(2) At Abnormal Operation. In this scenario, the system is exposed to different fault types, sections, and fault resistance (Rf) as it was for the status in the first case scenarios at an abnormal operation. Table 5 gives the PSCC angle values and its sign at LG fault at S3 with Rf = 0.02 Ω.Table 5
The PCSS angle for all relays at LG fault at S3 (case 2).
Relay numberPSCC angleFlow directionR172.9926ForwardR272.7197ForwardR371.4626ForwardR4−98.572ReverseR5−98.299ReverseR6−98.009ReverseR7−93.952ReverseR8−92.143ReverseR9−90.462ReverseR10−89.938ReverseR11−89.607ReverseR12−89.163ReverseR13−88.721ReverseR14126.150ForwardR15124.540ForwardR16123.919ForwardR17122.539ForwardR18125.552ForwardR19125.408ForwardR20125.397ForwardR21125.387ForwardR22126.912ForwardR23126.716ForwardR24126.690ForwardR25−101.24ReverseR26−100.30ReverseR27−99.777ReverseR28−99.236ReverseR29−98.234ReverseR30−95.968ReverseR31120.632ForwardR32116.911ForwardBy following the same procedures as normal operation scenario of case two, results given in Table5 can be analyzed as follows:(i)
During the fault, substation feeds the fault, so PSCC angle of R1, R2, and R3 have a positive sign (forward direction)(ii)
During the fault, ER1 and DER2 feeds the fault, so the PSCC angle sign is negative (reverse direction) for R4, R5, R6, R7, R8, R9, R10, R11, R12, R13, R25, R26, R27, R28, R29, and R30Figures14(a)–14(e) show the PSCC angle values and its sign at LG fault at S8 with Rf = 100 Ω, LG fault at S16 with Rf = 0.02 Ω, LLG fault at S20 with Rf = 0.02, LLLG fault at S24 with Rf = 0.02 Ω, and LL fault at S29 with Rf = 0.02 Ω, respectively. By following the above procedures, the direction of the faulty current is clearly shown to any faulted section.Figure 14
PSCC angle values and their sing at all relays for (a) LG fault at S8, (b) LG fault at S16, (c) LLG fault at S20, (d) LLLG fault at S24, and (e) LL fault at S29.
(a)(b)(c)(d)(e)According to the previous results, the PSCC angle is one of the most important factors in determining the fault current direction. Therefore, a simple circuit or scheme should be added to the OCR to help it in determining the FD and take the right decision, and this proposed scheme on MATLAB/Simulink is shown in Figure15.Figure 15
Proposed protection scheme on MATLAB/Simulink.In accordance with Figure15, if a fault occurs and the angle is positive, the CB trip signal shall come from the OCR which was setting by the optimal control variable, but if it is negative, the trip signal will come from the reverse coordination or the control circuit (automation system) to isolate faulty section, and this will be the future work.


## 6.1. Optimal Coordination

Firstly, the DERs, which are combined with the test system, are not energized. As a result, the power flow and the short circuit current flow are in the same direction. Current flows from the substation to the load point or fault location. Figure6 demonstrate the load current passing through each OCR in the operating mode.Figure 6
Load current at each bus.The measured three-phase fault current (symmetrical fault current) passing through each OCR at near end and far end for each primary relay is given in Table1. According to Table 1, the primary and backup relay can see approximately the same fault current magnitude for a one fault source because this operating condition has a radial topology. Also, the fault current level is increased when the fault point is near the source (substation) and decreases when fault point is far from the source.Table 1
Three-phase fault current at near end and far end for each primary relay.
Faulted sectionPrimary relayBackup relayFault currentNearFarS1R1—11201100S2R2R111001058S3R3R210581021S4R4R31021983S5R5R4983875S6R6R5875811S7R7R6811768S8R8R7768684S9R9R8684614S10R10R9614606S11R11R10606590S12R12R11590512S13R13R12512480S14R14R13480455S15R15R14455430S16R16R15430379S17R17R16379360S18R18R111001075S19R19R181075858S20R20R19858801S21R21R20801710S22R22R21021999S23R23R22999882S24R24R23882785S25R25R5875857S26R26R25857832S27R27R26832722S28R28R27722658S29R29R28658630S30R30R29630561S31R31R30561541S32R32R31541515To achieve the relay coordination optimization goal, the control variables (TMS,PS) of all relays are acquired in order to clear the near-end fault and far-end fault, which means there are 64 decision variables for coordination (i.e., (TMS1-TMS32) and (PS1-PS32)). The relays’ operating time range is considered between 0.04–1.1 sec. The range of TMS are limited from 0.05–1.1 sec., while the PS lies between 1.25–1.5 and the minimum value of CTI is considered 0.2 sec [39].To verify the uniqueness and efficacy of the IWHO technique, other metaheuristic techniques based on the same MOF have been used to calculate the coordination issue for the sake of comparison. These algorithms are WHO, a mayfly optimization algorithm (MOA), and GSA. The results of the optimal controlled variables and OF obtained using the IWHO technique in comparison with other techniques for the IEEE 33 bus test system are tabulated in Tables2. By analyzing Table 2, it is found that the minimum operating time for the system is improved to 19.57667 sec when using IWHO algorithm, so these results were selected to set up its control variables on the IDMT relays on MATLAB/Simulink. Figure 7 also depicts the novelty of the IWHO algorithm, which converges rapidly (around 2 iteration) before the other algorithms to attain the optimal fitness function.Table 2
Optimal setting variables’ results for the test system’s OCRs.
RelayMethodGSAAMWHOIWHOTMSPSTMSPSTMSPSTMSPSR11.11.4484481.11.51.0817421.3050541.0318411.447165R21.0999371.4864621.0999521.50.8831911.3273580.885031.271911R31.0716711.4091981.092481.50.8818121.250110.7804661.499995R41.0674671.2825651.0664561.3909990.8460711.3479170.7364141.259884R51.0341371.4565951.0662021.4987460.833961.2885970.7130741.317918R61.0129931.4584471.0449491.50.7219011.4996760.5506331.499361R70.9885111.3272381.0315641.50.6673331.250.4536211.457282R80.8551031.4024450.981041.431290.5345321.3202220.3398961.27922R90.713011.3918720.9069791.4163060.4857871.4980190.3023151.250533R100.5736421.4844570.8228261.3009260.4495191.3058250.2979991.250022R110.4577611.4431170.6991211.444370.3165691.2939690.2819771.252573R120.3541871.4232150.5967871.250.2397731.4518120.2364781.289341R130.276941.3885240.4747571.4101560.2135141.253950.2304031.251034R140.2018021.4133530.3871941.2817020.1746921.458590.2148611.262904R150.136071.2806950.2822281.50.1629541.2721850.1866811.250296R160.0776821.2835760.1984341.250.0883181.2508720.1583321.322459R170.0500041.3566870.051.3048790.0510021.2506180.0504471.290968R181.0999811.4287941.11.51.0987541.2500381.0995241.498037R190.7618661.3447030.7531261.4526970.7047081.4263690.730441.313271R200.4332431.383540.4133231.4043030.3802381.2629670.36821.45484R210.0504461.326130.0500031.2502190.0502571.2514130.0514951.295263R221.0487731.3026531.0828141.4461460.7726771.4212780.7720551.499657R230.5490781.4225140.6227681.4264660.4040251.3008980.4212991.250451R240.1046521.4420780.1895161.3262170.0501631.250.0514171.250342R251.0826191.3055541.11.4985170.6588431.3118740.5687551.489642R260.8524621.4155920.965161.3962690.4521841.499930.4008491.481262R270.7035531.3137790.8320991.50.4018791.2587370.2407881.368357R280.5396941.4349310.7018831.3996960.3108731.2530390.2261641.342916R290.4027291.4073460.554111.410890.2285181.2620810.1732031.499112R300.2762631.3431290.3922781.2506060.1600471.250010.1253881.252482R310.1819671.4024740.2559331.2505280.125771.2505340.0846381.271112R320.051.3625530.051.3490390.0500961.2517560.0500361.476527Total operating time (sec)30.4523135.1855722.1799919.57667Figure 7
Convergence characteristic of proposed algorithms.


## 6.2. The Proposed Protection Scheme

After completing the optimization work to achieve optimal coordination on MATLAB/coding and depending on the best algorithm (IWHO), its results (TMS, PS) were taken and set up on IDMT relays on the studied system. During the simulation on this system, the PSCC angle was monitored at each relay for a number of scenarios that fall under two cases.


### 6.2.1. Case 1: Without DERs’ Penetration

(1) At Normal Operation. In this scenario, the network has one current flow direction, which is from upstream to downstream. By analyzing Figure 8, it is found that the values of PSCC angle take different values, but keep the same positive sign, and this will be the base and considered as forward direction.Figure 8
PCSS angle values at normal operation (case 1).(2) At Abnormal Operation. In this scenario, the system is exposed to different fault types, sections, and fault resistance (Rf). In all fault types, the PSCC angle takes the same phenomena as the previous result, i.e., PSCC angle has a different value (affected by fault section) and a positive sign, which indicates that the direction of the current flow is also forward.Figures9–13 show the PSCC angle values and its sign at single line to ground (LG) fault at S3 with Rf  = 0.02 Ω, LG fault at S16 with Rf = 1000 Ω, double line to ground (LLG) fault at S20 with Rf = 0.02 Ω, three line to ground (LLLG) fault at S24 with Rf = 0.02 Ω, and line to line (LL) fault at S29 with Rf = 0.02 Ω, respectively. Note that the PSCC value of the buses after the fault, i.e., which is on the downstream side, is almost zero, but the PSCC angle value is recorded and plotted as observed from the simulation.Figure 9
PCSS angle values when LG fault occurs on S3 (case 1).Figure 10
PCSS angle values when LG fault occurs on S16 (case 1).Figure 11
PCSS angle values when LLG fault occurs on S20 (case 1).Figure 12
PCSS angle values when LLLG fault occurs on S24 (case 1).Figure 13
PCSS angle values when LL fault occurs on S29 (case 1).


### 6.2.2. Case 2: With DERs’ Penetration

The system under investigation contains two DERs. DERs allocation (place and capacity) is based on reducing losses using optimization techniques. The DERs places, capacities, and the reduction of losses are given in Table3 as stipulated in [38]. It is clear from Table 3 that the DERs will provide 60% of the total load.Table 3
DERs allocation and power losses reduction.
DER numberDER locationDER sizePower loss reduction (%)MWMVAR1130.84650.396186.42301.13921.0611(1) At Normal Operation. In this scenario, the DERs are integrated with the network, so the network does not have a single current flow direction as it was in the first case. By analyzing Table 4, it was found that the PSCC angle values take different values, and at some relays, it has a positive sign and on the other has a negative sign. The positive sign means that the current flow from upstream to downstream, similar to the first case is the forward direction. As for the negative sign, it means that the direction of the current passage is reversed due to the presence of the DERs.Table 4
PCSS angle for all relays at normal operation (case 2).
Relay numberPSCC angleFlow directionR1124.1610ForwardR2124.1679ForwardR3123.0042ForwardR4115.7648ForwardR5114.2543ForwardR6111.6414ForwardR798.21968ForwardR824.69574ForwardR9−29.6059ReverseR10−33.1408ReverseR11−36.4600ReverseR12−39.3560ReverseR13−41.4905ReverseR14131.8906ForwardR15130.2809ForwardR16129.6603ForwardR17128.2808ForwardR18128.1673ForwardR19128.1388ForwardR20128.1286ForwardR21128.1181ForwardR22127.4158ForwardR23127.1179ForwardR24127.0909ForwardR25122.6414ForwardR26120.1087ForwardR2782.25281ForwardR28−34.9689ReverseR29−49.7763ReverseR30−60.4110ReverseR31126.0625ForwardR32122.3451ForwardTo clarify more, the PSCC angle for relays R9, R10, R11, R12, and R13 has a negative sign, which means that the loads at buses 9, 10, 11, 12, and 13, which these relays are located at, are fed by the DER1 that is located at bus 13 and not by the main source (substation). Of course, it is certain that the loads at buses 14, 15, 16, 17, and 18 are fed by this DER1 as well, but the direction of the current is the same as the original direction (forward direction), so its PSCC angle will take the positive sign as it was. By analogy with that, the loads at buses 28, 29, and 30 are fed from the DER2 that is located at bus 30, so they take the negative sign; as well as, loads at buses 31 and 32 are fed from this DER2, but the direction of the current is the same as the original direction, so it takes the positive sign.(2) At Abnormal Operation. In this scenario, the system is exposed to different fault types, sections, and fault resistance (Rf) as it was for the status in the first case scenarios at an abnormal operation. Table 5 gives the PSCC angle values and its sign at LG fault at S3 with Rf = 0.02 Ω.Table 5
The PCSS angle for all relays at LG fault at S3 (case 2).
Relay numberPSCC angleFlow directionR172.9926ForwardR272.7197ForwardR371.4626ForwardR4−98.572ReverseR5−98.299ReverseR6−98.009ReverseR7−93.952ReverseR8−92.143ReverseR9−90.462ReverseR10−89.938ReverseR11−89.607ReverseR12−89.163ReverseR13−88.721ReverseR14126.150ForwardR15124.540ForwardR16123.919ForwardR17122.539ForwardR18125.552ForwardR19125.408ForwardR20125.397ForwardR21125.387ForwardR22126.912ForwardR23126.716ForwardR24126.690ForwardR25−101.24ReverseR26−100.30ReverseR27−99.777ReverseR28−99.236ReverseR29−98.234ReverseR30−95.968ReverseR31120.632ForwardR32116.911ForwardBy following the same procedures as normal operation scenario of case two, results given in Table5 can be analyzed as follows:(i)
During the fault, substation feeds the fault, so PSCC angle of R1, R2, and R3 have a positive sign (forward direction)(ii)
During the fault, ER1 and DER2 feeds the fault, so the PSCC angle sign is negative (reverse direction) for R4, R5, R6, R7, R8, R9, R10, R11, R12, R13, R25, R26, R27, R28, R29, and R30Figures14(a)–14(e) show the PSCC angle values and its sign at LG fault at S8 with Rf = 100 Ω, LG fault at S16 with Rf = 0.02 Ω, LLG fault at S20 with Rf = 0.02, LLLG fault at S24 with Rf = 0.02 Ω, and LL fault at S29 with Rf = 0.02 Ω, respectively. By following the above procedures, the direction of the faulty current is clearly shown to any faulted section.Figure 14
PSCC angle values and their sing at all relays for (a) LG fault at S8, (b) LG fault at S16, (c) LLG fault at S20, (d) LLLG fault at S24, and (e) LL fault at S29.
(a)(b)(c)(d)(e)According to the previous results, the PSCC angle is one of the most important factors in determining the fault current direction. Therefore, a simple circuit or scheme should be added to the OCR to help it in determining the FD and take the right decision, and this proposed scheme on MATLAB/Simulink is shown in Figure15.Figure 15
Proposed protection scheme on MATLAB/Simulink.In accordance with Figure15, if a fault occurs and the angle is positive, the CB trip signal shall come from the OCR which was setting by the optimal control variable, but if it is negative, the trip signal will come from the reverse coordination or the control circuit (automation system) to isolate faulty section, and this will be the future work.


## 6.2.1. Case 1: Without DERs’ Penetration

(1) At Normal Operation. In this scenario, the network has one current flow direction, which is from upstream to downstream. By analyzing Figure 8, it is found that the values of PSCC angle take different values, but keep the same positive sign, and this will be the base and considered as forward direction.Figure 8
PCSS angle values at normal operation (case 1).(2) At Abnormal Operation. In this scenario, the system is exposed to different fault types, sections, and fault resistance (Rf). In all fault types, the PSCC angle takes the same phenomena as the previous result, i.e., PSCC angle has a different value (affected by fault section) and a positive sign, which indicates that the direction of the current flow is also forward.Figures9–13 show the PSCC angle values and its sign at single line to ground (LG) fault at S3 with Rf  = 0.02 Ω, LG fault at S16 with Rf = 1000 Ω, double line to ground (LLG) fault at S20 with Rf = 0.02 Ω, three line to ground (LLLG) fault at S24 with Rf = 0.02 Ω, and line to line (LL) fault at S29 with Rf = 0.02 Ω, respectively. Note that the PSCC value of the buses after the fault, i.e., which is on the downstream side, is almost zero, but the PSCC angle value is recorded and plotted as observed from the simulation.Figure 9
PCSS angle values when LG fault occurs on S3 (case 1).Figure 10
PCSS angle values when LG fault occurs on S16 (case 1).Figure 11
PCSS angle values when LLG fault occurs on S20 (case 1).Figure 12
PCSS angle values when LLLG fault occurs on S24 (case 1).Figure 13
PCSS angle values when LL fault occurs on S29 (case 1).


## 6.2.2. Case 2: With DERs’ Penetration

The system under investigation contains two DERs. DERs allocation (place and capacity) is based on reducing losses using optimization techniques. The DERs places, capacities, and the reduction of losses are given in Table3 as stipulated in [38]. It is clear from Table 3 that the DERs will provide 60% of the total load.Table 3
DERs allocation and power losses reduction.
DER numberDER locationDER sizePower loss reduction (%)MWMVAR1130.84650.396186.42301.13921.0611(1) At Normal Operation. In this scenario, the DERs are integrated with the network, so the network does not have a single current flow direction as it was in the first case. By analyzing Table 4, it was found that the PSCC angle values take different values, and at some relays, it has a positive sign and on the other has a negative sign. The positive sign means that the current flow from upstream to downstream, similar to the first case is the forward direction. As for the negative sign, it means that the direction of the current passage is reversed due to the presence of the DERs.Table 4
PCSS angle for all relays at normal operation (case 2).
Relay numberPSCC angleFlow directionR1124.1610ForwardR2124.1679ForwardR3123.0042ForwardR4115.7648ForwardR5114.2543ForwardR6111.6414ForwardR798.21968ForwardR824.69574ForwardR9−29.6059ReverseR10−33.1408ReverseR11−36.4600ReverseR12−39.3560ReverseR13−41.4905ReverseR14131.8906ForwardR15130.2809ForwardR16129.6603ForwardR17128.2808ForwardR18128.1673ForwardR19128.1388ForwardR20128.1286ForwardR21128.1181ForwardR22127.4158ForwardR23127.1179ForwardR24127.0909ForwardR25122.6414ForwardR26120.1087ForwardR2782.25281ForwardR28−34.9689ReverseR29−49.7763ReverseR30−60.4110ReverseR31126.0625ForwardR32122.3451ForwardTo clarify more, the PSCC angle for relays R9, R10, R11, R12, and R13 has a negative sign, which means that the loads at buses 9, 10, 11, 12, and 13, which these relays are located at, are fed by the DER1 that is located at bus 13 and not by the main source (substation). Of course, it is certain that the loads at buses 14, 15, 16, 17, and 18 are fed by this DER1 as well, but the direction of the current is the same as the original direction (forward direction), so its PSCC angle will take the positive sign as it was. By analogy with that, the loads at buses 28, 29, and 30 are fed from the DER2 that is located at bus 30, so they take the negative sign; as well as, loads at buses 31 and 32 are fed from this DER2, but the direction of the current is the same as the original direction, so it takes the positive sign.(2) At Abnormal Operation. In this scenario, the system is exposed to different fault types, sections, and fault resistance (Rf) as it was for the status in the first case scenarios at an abnormal operation. Table 5 gives the PSCC angle values and its sign at LG fault at S3 with Rf = 0.02 Ω.Table 5
The PCSS angle for all relays at LG fault at S3 (case 2).
Relay numberPSCC angleFlow directionR172.9926ForwardR272.7197ForwardR371.4626ForwardR4−98.572ReverseR5−98.299ReverseR6−98.009ReverseR7−93.952ReverseR8−92.143ReverseR9−90.462ReverseR10−89.938ReverseR11−89.607ReverseR12−89.163ReverseR13−88.721ReverseR14126.150ForwardR15124.540ForwardR16123.919ForwardR17122.539ForwardR18125.552ForwardR19125.408ForwardR20125.397ForwardR21125.387ForwardR22126.912ForwardR23126.716ForwardR24126.690ForwardR25−101.24ReverseR26−100.30ReverseR27−99.777ReverseR28−99.236ReverseR29−98.234ReverseR30−95.968ReverseR31120.632ForwardR32116.911ForwardBy following the same procedures as normal operation scenario of case two, results given in Table5 can be analyzed as follows:(i)
During the fault, substation feeds the fault, so PSCC angle of R1, R2, and R3 have a positive sign (forward direction)(ii)
During the fault, ER1 and DER2 feeds the fault, so the PSCC angle sign is negative (reverse direction) for R4, R5, R6, R7, R8, R9, R10, R11, R12, R13, R25, R26, R27, R28, R29, and R30Figures14(a)–14(e) show the PSCC angle values and its sign at LG fault at S8 with Rf = 100 Ω, LG fault at S16 with Rf = 0.02 Ω, LLG fault at S20 with Rf = 0.02, LLLG fault at S24 with Rf = 0.02 Ω, and LL fault at S29 with Rf = 0.02 Ω, respectively. By following the above procedures, the direction of the faulty current is clearly shown to any faulted section.Figure 14
PSCC angle values and their sing at all relays for (a) LG fault at S8, (b) LG fault at S16, (c) LLG fault at S20, (d) LLLG fault at S24, and (e) LL fault at S29.
(a)(b)(c)(d)(e)According to the previous results, the PSCC angle is one of the most important factors in determining the fault current direction. Therefore, a simple circuit or scheme should be added to the OCR to help it in determining the FD and take the right decision, and this proposed scheme on MATLAB/Simulink is shown in Figure15.Figure 15
Proposed protection scheme on MATLAB/Simulink.In accordance with Figure15, if a fault occurs and the angle is positive, the CB trip signal shall come from the OCR which was setting by the optimal control variable, but if it is negative, the trip signal will come from the reverse coordination or the control circuit (automation system) to isolate faulty section, and this will be the future work.


## 7. Conclusion

In this study, protection systems for DNs have been addressed on two sides. The first side is the optimum coordination for DN protection devices (IDMT OCRs), where this DN does not contain DERs. A comparison was made between a number of metaheuristics techniques, which are GSA, MA, WHO, and IWHO, and the best of them were chosen in terms of the least operating time and the fastest converges, which is IWHO. TheTMS and PS outputs of the IWHO algorithm were used to set IDMT OCRs on the MATLAB/Simulink.DERs are optimally integrated into the network, having a share rate of 60% from the total load in preparation for the second side. In the second side, a novel scheme was suggested and presented to help the OCRs to make the right decision with the DERs penetration. It was proposed because PSCC angle sign is changed if the direction of fault current is reversed. To identify the FD based on this property, a technique based on the sign of this angle has been created and implemented for each OCR. In light of this, if the PSCC angle is positive, the trip operating time is determined by the relay settings that have been optimally tuned in advance for each relay. If the sign is negative, it indicates that the flow is going in the opposite direction; therefore, the relay does not trip until it receives a signal from the reverse coordination or control circuit (automation system), which will be the future task. The main benefit of the proposed protection scheme is that it relies only on the current measurements from one end and, at the same time, reduces the cost because there is no potential transformer or compensation technique required. The proposed protective solutions are evaluated and tested on IEEE 33 bus system. The results showed the effectiveness of the proposed IWHO technique in terms of accuracy and speed and belay that the proposed protection scheme is unaffected by the fault location, transition resistance, and different types of faults.



---

*Source: 7630011-2022-11-07.xml*