An Optimized Adaptive Protection Scheme for Numerical and Directional Overcurrent Relay Coordination Using Harris Hawk Optimization

: The relay coordination problem is of dire importance as it is critical to isolate the faulty portion in a timely way and thus ensure electrical network security and reliability. Meanwhile a relay protection optimization problem is highly constraint and complicated problem to be addressed. To fulﬁll this purpose, Harris Hawk Optimization (HHO) is adapted to solve the optimization problem for Directional Over-current Relays (DOCRs) and numerical relays. As it is inspired by the intelligent and collegial chasing and preying behavior of hawks for capturing the prey, it shows quite an impressive result for ﬁnding the global optimum values. Two decision variables; Time Dial Settings (TDS) and Plug Settings (PS) are chosen as the decision variables for minimization of overall operating time of relays. The proposed algorithm is implemented on three IEEE test systems. In comparison to other state-of-the-art nature inspired and traditional algorithms, the results demonstrate the superiority of HHO.


Introduction
The electrical power system is one of the most crucial systems running across the globe. For smooth operation of the electrical network, effective protection systems are a necessary requirement. In an electrical power system, the primary function of the protective system is to detect and isolate any failed or faulty components as fast as can be, so that the unfaulty portion continues to be operational. To ensure reliability and security, the components of the electrical network are protected by primary as well as backup protection. In an event of fault, the primary protection must operate in a timely way to confine the faulty part of the electrical network. In case of primary protection failure, backup protection must act to accomplish the protection task. This ought to be the favored situation of any protective system in light of the fact that the primary protection confines just the influenced region while the backup protection ensures that at whatever point, no more than a strictly necessary portion of the systems is exposed to the ill effects of blackouts. To ensure that only the affected segment of the system is confined thusly diminishing the likelihood of bothersome broader power outages, dependable and viable administration of protection equipment is required.
For a multi-loop network, a viable and productive protection scheme needs to consider Directional Over-current Relays (DOCRs). The functioning and configuration of a DOCR depend upon two parameters. The first one is the Time Dial Setting (TDS) and the second one is the Plug Setting (PS). The target of DOCR optimization is to provide the optimal settings for TDS and PS such that the primary relays respond promptly to any fault in their zone. Backup relays should also operate in a timely manner and not affect the unfaulty portions of the network. Therefore, minimization of operating times of primary relays as well as coordination of backup relays are the requirements of DOCRs, while obeying the constraints. To explain this confounded issue, various strategies have been proposed in the literature. The curve intersection method was utilized to take care of the coordination issue of overcurrent relays [1]. In [2], the graphical selection method was used to find the relay settings. The minimum break point set method has been used for an expert system in [3]. For DOCR coordination problems [4][5][6][7][8] and complex problems within the electrical power network, numerous techniques like the Nature-inspired optimization technique whale optimization [9] and the JAYA algorithm [10] have been developed to deal with the DOCR problem. With respect to the physical and logical change of arrangement of the system the optimal DOCR design was developed in [11,12]. The optimal TDS and PS settings of a DOCR were found using the hybrid genetic algorithm for the generation unit or line outage contingencies in [13,14]. The DOCRs problem in a multi loop transmission system due to single line outage contingencies has been resolved in [15,16]. A grid-connected photovoltaic system for regulating power quality issues using a supercapacitor-based STATCOM is used in [17]. The protection coordination problem is vital in the design of DOCRs, which sometime operate improperly due to alterations in the network configuration. In [18], fault current limiters were used to solve the protection coordination problem without considering line or DG outages. For line, substation and DG outages in microgrids, the protection settings were determined by using the new coordination parameters set [19]. Moreover, all these network arrangements makes them more complex, yielding discoordination of DOCRs, that constantly leads to contingencies. In order to resolve these issues and the power outages, optimal relay setting while considering the main system configuration has been recommended in [20][21][22][23]. The DOCRs problem was formulated as a mixed integer nonlinear programming (MINLP) problem and relay setting parameters were found using different populationbased optimization techniques in [24,25]. In [26,27] a few bio-motivated algorithms were developed to tackle the DOCR coordination issue by designing a linear formulation. In [28][29][30][31][32], a different version of particle swarm optimization (PSO) was used to determine the optimum values for DOCRs. A different version of the differential algorithm was reported in [33] to solve the DOCR coordination problem to point out the superiority of modified differential evolution algorithms. Many other Nature-inspired algorithms like the grey wolf optimizer (GWO), teaching learning-based optimization (TLBO), biography-based optimization (BBO), back-tracking algorithm, the improved firefly (IFA) metaheuristic and modified electromagnetic field optimization (MEFO) were used for DOCR coordination in [34][35][36][37][38][39][40]. A modified teaching-based optimization algorithm was implemented in [41]. An analytic approach to solve the DOCR coordination problem was utilized in [42]. In [43], for determining relay setting parameters, an improved group search algorithm was used. In [44], the comparison of several metaheuristic algorithms to solve DOCR problem was presented. In [45], in order to solve DOCR problem, multiple embedded crossover PSO algorithm was applied. For a multisource network, the DOCR problem can be regarded as an optimization problem. The downside of the past optimization procedures and that of the metaheuristic optimization, is that they possibly converge to settings which are not optimal and may trap into local optimum results. The mathematical optimization techniques work on the gradient based information about the used functions for finding solutions. For such type of algorithms, there is likelihood of converging to local optimum results. Further, convergence rate of such algorithms lowers with larger systems. These drawbacks are alleviated by nature-inspired algorithms. Such algorithms start optimization from random solutions, and this helps in avoidance of local optima. These algorithms also do not need the gradient-based information for the solution. The edge of Nature-inspired algorithms is that they can be incorporated into optimization solutions considering them as a black box. This helps in finding optimal solutions easily for problems with unknown search space [46,47]. To comprehend this issue, a HHO algorithm technique is examined in this work to determine the optimal DOCR parameters as compared to other state-of-the-art algorithms.
In this paper, the optimal settings for DOCRs were found by a novel populationbased, Nature-inspired optimization algorithm called the Harris Hawks Optimizer (HHO) deployed in a multi-loop power system. The primary motivation of HHO is the agreeable conduct and pursuing style of Harris hawks in Nature called surprise pounce. Furthermore, according to the no free lunch theorem, all proposed algorithms give almost equal results on average, when applied to all optimization problems [48]. In short, no single algorithm is universally good to solve all optimization problems, and this urges researchers to adapt and use more efficient optimization techniques. HHO is developed under the inspiration of the cooperative schema of hunting and chasing the prey by performing surprise pounces. The chasing style of Harris hawks changes dynamically according to the behavior of the prey. This paper proposes the HHO methodology that can accomplish the optimal DOCR coordination and operation. The implementation of the HHO is utilized for a test in the IEEE-9, 15 and 14-bus systems. The comparison of HHO with other algorithms proves its effectiveness. In this insightful procedure, a few Harris hawks cooperate to pounce on a prey from various bearings trying to surprise it. Harris hawks can reveal a variety of chasing patterns dependent on the dynamic nature of scenarios and the escape patterns of the prey. This work numerically impersonates such dynamic patterns and behavior to develop an optimization algorithm. The suggested HHO has more exploration capability in contrast to other metaheuristic algorithms. These qualities of HHO enhance the potential of search agents to seek optimal solutions. The important part of the proposed HHO is to decide the ideal estimations of TDS and PS to minimize the operational time of DOCRs regarding reinforcement and hand-off setting limitations.
The paper is organized as follows: Section 2 presents the problem formulation for DOCR. Section 3 contains an introduction to the proposed algorithm and its details. Section 4 explains about the results in the used test systems and a detailed discussion about the results. Finally, Section 5 concludes the paper.

DOCR Problem Formulation
The DOCR relay system is supposed to sense a fault event and isolate the faulty portion of a network in a timely manner, so that the healthy portion of the system does not experience the effects of the fault. Two design variables, TDS and PS, are selected to find the optimized and minimized operating time. The optimum values of TDS and PS will result in a minimized collective operating time of relays. The sum of operating times of the relays is referred to as an objective function as shown in Equation (1): where T ji indicates the operating times of relays for a fault in j zone and i varies from 1 to the number of relays in the system. The IEC normal inverse relays are selected for relay protection schema and operating time is given by Equation (2): where α and k are characteristic constants with values of 0.14 and 0.02, respectively [49]. Furthermore, TDS indicates time dial settings, IF shows the fault current value, PS shows plug settings value and CTR shows the currents transformer ratio value for the relays. T j stands for the operating time of the primary relay. The goal is to minimize the operating time of relays, while keeping constraints into consideration. The overall schema of DOCRs coordination and protection is shown in Figure 1. Which shows the goal of DOCRs protection keeping in view the backup protection as well as constraints.
plug settings value and CTR shows the currents transformer ratio value for the relays. Tj stands for the operating time of the primary relay. The goal is to minimize the operating time of relays, while keeping constraints into consideration. The overall schema of DOCRs coordination and protection is shown in Figure 1. Which shows the goal of DOCRs protection keeping in view the backup protection as well as constraints.

Coordination Criteria
For an effective protection scheme, backup protection is supposed to act as a second line of defense. It should be ensured that backup protection acts neither too fast nor too slow. Hence, backup protection is supposed to act after some coordination time interval (CTI). This value varies between 0.2-0.5 s, depending upon the choice of certain value or relay types used. The relation between backup and primary relays is shown by Equation (3): where Tb: the backup relay operating time; and Tj: the primary (or main) relay operating time.

Relay Setting Bounds
Along with coordination constraints, the relay protection scheme is subject to constraints imposed on the values of TDS and PS. Further information about minimum and maximum values of TDS and PS is mentioned along with the experimental data. The number of TDS and PS varies according to the number of primary relays. In following equations, i indicates the number of relays. The optimum values of TDS and PS are bound to obey these settings constraints of minimum and maximum values. The ultimate goal is the settings which will ensure constraints of coordination as well as settings constraints, indicated in this section, as shown by Equations (4) and (5):

Harris Hawks Optimization (HHO) Algorithm
In 2019, Heidari et al. presented a population-based Nature-inspired model-based algorithm called the Harris Hawk Optimization (HHO) algorithm to solve optimization problems [50]. The HHO was developed under the inspiration of one of highest IQ hunters, the Harris hawk. Harris Hawks sense and adapt their hunting strategies according to the dynamic moves and escape attempts of their prey. The hunting quest of several hawks

Coordination Criteria
For an effective protection scheme, backup protection is supposed to act as a second line of defense. It should be ensured that backup protection acts neither too fast nor too slow. Hence, backup protection is supposed to act after some coordination time interval (CTI). This value varies between 0.2-0.5 s, depending upon the choice of certain value or relay types used. The relation between backup and primary relays is shown by Equation (3): where T b : the backup relay operating time; and T j : the primary (or main) relay operating time.

Relay Setting Bounds
Along with coordination constraints, the relay protection scheme is subject to constraints imposed on the values of TDS and PS. Further information about minimum and maximum values of TDS and PS is mentioned along with the experimental data. The number of TDS and PS varies according to the number of primary relays. In following equations, i indicates the number of relays. The optimum values of TDS and PS are bound to obey these settings constraints of minimum and maximum values. The ultimate goal is the settings which will ensure constraints of coordination as well as settings constraints, indicated in this section, as shown by Equations (4) and (5):

Harris Hawks Optimization (HHO) Algorithm
In 2019, Heidari et al. presented a population-based Nature-inspired model-based algorithm called the Harris Hawk Optimization (HHO) algorithm to solve optimization problems [50]. The HHO was developed under the inspiration of one of highest IQ hunters, the Harris hawk. Harris Hawks sense and adapt their hunting strategies according to the dynamic moves and escape attempts of their prey. The hunting quest of several hawks starts with mutual encircling of the prey and then changing different hunting patterns according to the behavior of the prey. Unlike other algorithms, HHO adapts dynamically according to different patterns shown by the prey in the exploitation phase. In this section, the exploration and exploitation of HHO is modelled, elaborating mutual hunting strategy of surprise pounces and attacking by Harris hawks.

Exploration Phase
In this part, the exploration process by HHO is stated. Harris hawks scan the search area with their powerful eyes. They sit on some tall place to better cover a large search area and may spend several hours in locating and identifying a prey. It is important to note that in the HHO paradigm, Harris hawks are the candidate solutions and the prey is the optimal solution. There is an equal probability that Harris hawks will adapt to two different perching strategies. According to the first strategy, a Harris hawk may perch according to its hunting group, while, as a second strategy, Harris hawks may perch randomly on a tall tree, inside the hunting group range. Both of these perching strategies are summed up in Equation (6), indicating q ≥ 0.5 for first case and q < 0.5 for second case: Here, X(t + 1) is the random location of hawks for the iteration t + 1. X(t), given by Equation (7), is the hawks' position vector at current iteration t. Here, N D represents the number of decision variables and R N represents the number of relays being used in the system. This equation will produce the (N D × R N ) vector of uniformly distributed numbers in the interval (LB, UB). X rand (t) is a random hawk from the current population and X m (t) is the average position of hawks for current iteration. The position of the rabbit (prey) in the current iteration is given by X rabbit (t) . To depict the random nature, variables r 1 , r 2 , r 3 , r 4 and r 5 are included from interval (0, 1). To restrict the scenario to the valid search space, all variables are supposed to obey the upper and lower bounds, included as UB and LB in the Equation (6). According to the first part of Equation (6), the position of a hawk for the next iteration will be the result of the difference between the random hawk's position from the current iteration and the normalized difference of the position of the random hawk and other hawks in the current iteration. According to the second part of Equation (6), the position of hawks for the next iteration is based on the best location so far, the average position and a randomly scaled factor based on the group range. To impart randomness and explore search space more, a randomly scaled LB parameter is included. By using Equation (7), the average position of hawks for current iteration is calculated: Here, N indicates the number of hawks and X i (t) is the individual hawk location for the current iteration.

Transition from Exploration to Exploitation
The switching from exploration to exploitation and then different dynamic exploitation modes is triggered by the escaping energy of the prey. The escaping energy decreases during the hunt, whereas different escaping patterns can be adopted by prey to deceive the hawks. The escaping energy is of the prey is mathematically modelled in Equation (9): Here, E o is the energy of prey at the initial state, given by Equation (10). t is the current iteration and T indicates the total number of iterations and r 6 is a random number which varies between 0 and 1. The energy of the prey varies randomly between −1 and 1. For energies between 0 and −1, the prey starts weakening while for energies between 0 and 1 the prey starts strengthening and performs deceiving moves and tries its level best to get Energies 2021, 14, 5603 6 of 20 away from the hunting hawks. As iterations pass, the energy of the prey decreases. For escaping energies |E| ≥ 1, the hawks look through a different vicinity to locate the prey, indicating exploration, whereas for escaping energies |E| < 1 a detailed search is carried out, referred to as exploitation. To sum things up, for |E| ≥ 1, the HHO is in exploration phase, while for |E| < 1, HHO is in an exploitation phase. The escaping energy for two runs and 500 iterations is shown in Figure 2. Here, Eo is the energy of prey at the initial state, given by Equation (10). t is the current iteration and T indicates the total number of iterations and r6 is a random number which varies between 0 and 1. The energy of the prey varies randomly between −1 and 1. For energies between 0 and −1, the prey starts weakening while for energies between 0 and 1 the prey starts strengthening and performs deceiving moves and tries its level best to get away from the hunting hawks. As iterations pass, the energy of the prey decreases. For escaping energies |E| ≥ 1, the hawks look through a different vicinity to locate the prey indicating exploration, whereas for escaping energies |E| < 1 a detailed search is carried out, referred to as exploitation. To sum things up, for |E| ≥ 1, the HHO is in exploration phase, while for |E| < 1, HHO is in an exploitation phase. The escaping energy for two runs and 500 iterations is shown in Figure 2.

Exploitation Phase
Harris hawks implement a surprise pouncing strategy on the identified prey. As a resultant the prey becomes exhausted and is captured easily by the hawks, but during this process, the prey tries its best to deceive the hawks and save itself from danger. To mode this scenario, four different scenarios are presented in HHO. The prey strives hard to get away from hawks. Before the surprise pounce, the probability of a successful escape is characterized as (r < 0.5), while the probability of an unsuccessful escape is characterized as (r ≥ 0.5). Therefore, depending upon the parameters r and E, the prey's escape patterns can be characterized into four different scenarios, explained in detail in the following section. According to the different scenarios, Harris hawks change their cooperative hunting tactics. As a result, the prey gets exhausted and its energy is diminished. In this stage Harris hawks perform a hard or soft siege and capture the prey. In such a manner, the adaptation mode of a Harris hawk is classified into soft and hard siege, depending upon the condition |E| ≥ 0.5 and |E| < 0.5, respectively.

Exploitation Phase
Harris hawks implement a surprise pouncing strategy on the identified prey. As a resultant the prey becomes exhausted and is captured easily by the hawks, but during this process, the prey tries its best to deceive the hawks and save itself from danger. To model this scenario, four different scenarios are presented in HHO. The prey strives hard to get away from hawks. Before the surprise pounce, the probability of a successful escape is characterized as (r < 0.5), while the probability of an unsuccessful escape is characterized as (r ≥ 0.5). Therefore, depending upon the parameters r and E, the prey's escape patterns can be characterized into four different scenarios, explained in detail in the following section. According to the different scenarios, Harris hawks change their cooperative hunting tactics. As a result, the prey gets exhausted and its energy is diminished. In this stage, Harris hawks perform a hard or soft siege and capture the prey. In such a manner, the adaptation mode of a Harris hawk is classified into soft and hard siege, depending upon the condition |E| ≥ 0.5 and |E| < 0.5, respectively.

Soft Siege
For the scenario of unsuccessful escape the probability r ≥ 0.5 and the energy value |E| ≥ 0.5, the prey still has enough energy. The prey performs random jumps to mislead the Harris hawks. In such a scenario, Harris hawks adapt a soft siege mode to further exhaust the prey. This scenario is modelled in Equation (11): Here ∆X(t) in Equations (11) and (12) represents the difference between the rabbit position and hawk position for the current iteration: Furthermore, Equation (13) indicates the jump strength of the prey whilst escaping. Here, r 6 is random number from the interval (0, 1). This random number imparts a random nature to demonstrate the movement of the prey in every iteration

Hard Siege
For the scenario of unsuccessful escape probability r ≥ 0.5 and energy value |E|< 0.5, the prey is quite exhausted and has considerably low energy to escape the hunt. Harris hawks sense this situation and their siege tactics are intensified to capture the prey. This scenario is modelled in Equation (14): An example to demonstrate this scenario is shown in Figure 3.
For the scenario of unsuccessful escape the probability r ≥ 0.5 and the energy value |E| ≥ 0.5, the prey still has enough energy. The prey performs random jumps to mislead the Harris hawks. In such a scenario, Harris hawks adapt a soft siege mode to further ex haust the prey. This scenario is modelled in Equation (11): Here ∆X(t) in Equations (11) and (12) represents the difference between the rabbi position and hawk position for the current iteration: Furthermore, Equation (13) indicates the jump strength of the prey whilst escaping Here, r6 is random number from the interval (0, 1). This random number imparts a random nature to demonstrate the movement of the prey in every iteration

Hard Siege
For the scenario of unsuccessful escape probability r ≥ 0.5 and energy value |E|< 0.5 the prey is quite exhausted and has considerably low energy to escape the hunt. Harris hawks sense this situation and their siege tactics are intensified to capture the prey. Thi scenario is modelled in Equation (14): An example to demonstrate this scenario is shown in Figure 3.

Soft Siege with Progressive Rapid Dives
For the scenario of successful escape probability r < 0.5 and energy value |E|≥ 0.5 prey is energetic enough to deceive the hunting plans. Sensing this situation, Harris hawk switch their hunting schema to a more sophisticated and intelligent mode. The mode o hunting is soft siege but it combined with the technique of Levy flights (LF). In accordance to the deceiving motions and jumps of the prey, Harris hawks also start zigzag dives [51] In this way, Harris hawks try to match their hunting pattern with the random and deceiv ing jumps of the prey. LF is a quite favorable mode of capturing the prey for foragers and predators [52,53]. The LF technique is often adopted by predators, and can be seen in the chasing schemes used by sharks and monkeys [51][52][53][54][55][56][57]. The soft siege model used by hawk at this stage is shown in Equation (15)

Soft Siege with Progressive Rapid Dives
For the scenario of successful escape probability r < 0.5 and energy value |E|≥ 0.5, prey is energetic enough to deceive the hunting plans. Sensing this situation, Harris hawks switch their hunting schema to a more sophisticated and intelligent mode. The mode of hunting is soft siege but it combined with the technique of Levy flights (LF). In accordance to the deceiving motions and jumps of the prey, Harris hawks also start zigzag dives [51]. In this way, Harris hawks try to match their hunting pattern with the random and deceiving jumps of the prey. LF is a quite favorable mode of capturing the prey for foragers and predators [52,53]. The LF technique is often adopted by predators, and can be seen in the chasing schemes used by sharks and monkeys [51][52][53][54][55][56][57]. The soft siege model used by hawks at this stage is shown in Equation (15): The LF technique used by HHO is modelled on Equation (16). Harris hawks analyze their current move and the past moves and change their diving and hunting strategy accordingly. After comparison and realization, Harris hawks perform sporadic, unexpected and abrupt dives: Here, D indicates the dimension, S indicates a random vector of size (1 × D) and LF indicates the Levy flight function, which is modelled by using Equation (17): Here, u and v are the random values between the interval (0, 1). β is a constant with a value equal to 1.5. As a whole, the hunting and surprise pouncing pattern used by the hawks for this mode can be modelled by using Equation (18): Here, the parameters Y and Z refer to Equations (15) and (16) respectively. This step with one hawk is elaborated in Figure 4. The LF patterns shown by the hawk are also demonstrated for a course of a few iterations. The colored-dotted line shows the LF-based dives. For every next iteration, Y and Z are adopted to intelligently capture the prey. This behavior is adopted by the whole group of Harris hawks.
The LF technique used by HHO is modelled on Equation (16). Harris hawks analyze their current move and the past moves and change their diving and hunting strategy accordingly. After comparison and realization, Harris hawks perform sporadic, unexpected and abrupt dives: Here, D indicates the dimension, S indicates a random vector of size (1 × D) and LF indicates the Levy flight function, which is modelled by using Equation (17): Here, u and v are the random values between the interval (0, 1). β is a constant with a value equal to 1.5. As a whole, the hunting and surprise pouncing pattern used by the hawks for this mode can be modelled by using Equation (18): Here, the parameters Y and Z refer to Equations (15) and (16) respectively. This step with one hawk is elaborated in Figure 4. The LF patterns shown by the hawk are also demonstrated for a course of a few iterations. The colored-dotted line shows the LF-based dives. For every next iteration, Y and Z are adopted to intelligently capture the prey. This behavior is adopted by the whole group of Harris hawks.

Hard Siege with Progressive Rapid Dives
For the scenario of successful escape probability r < 0.5 and energy value |E|< 0.5, the prey is exhausted because of low energy. Now, Harris hawks decrease the distance between the prey and perform a hard siege. To tackle the deceiving jumps of the prey, the LF-based technique is also incorporated. This situation is modelled in Equation (19):

Hard Siege with Progressive Rapid Dives
For the scenario of successful escape probability r < 0.5 and energy value |E|< 0.5, the prey is exhausted because of low energy. Now, Harris hawks decrease the distance between the prey and perform a hard siege. To tackle the deceiving jumps of the prey, the LF-based technique is also incorporated. This situation is modelled in Equation (19): where Y and Z are obtained using the new rules in Equations (20) and (21).
Here, X m (t) can be accessed from Equation (7). This scenario is elaborated in Figure 5. The colored-dots indicate the LF patterns for the current iteration. For the next iteration, the decision between Y or Z will be made by the hawks, accordingly. where Y and Z are obtained using the new rules in Equations (20) and (21).
Here, Xm(t) can be accessed from Equation (7). This scenario is elaborated in Figure  5. The colored-dots indicate the LF patterns for the current iteration. For the next iteration the decision between Y or Z will be made by the hawks, accordingly.
(a) The process in 2D space (b) The process in 3D Space Figure 5. The position vectors for hard siege mode with progressive rapid dives in 2D and 3D space. [50]. Figure 5. The position vectors for hard siege mode with progressive rapid dives in 2D and 3D space [50].

Computational Complexity
The computational complexity is the result of three processes. The first one is initialization, the second one is fitness evaluation and the third one is updating of hawks. With N number of hawks, the initialization complexity is O(N). Due to updating, the complexity becomes O(T × N) + O(T × N × D). This includes the search for the best location and as well as updating of location for whole group of hawks. Here, T indicates the total number of iterations and D indicates the dimension of problems. Resultantly, the computational complexity can be given as O(N × (T + TD + 1)).

Results and Discussion
In this part, the HHO is successfully implemented to address the DOCR coordination problem and has been verified for three IEEE standard test systems: the IEEE-9, 14, and 15 bus systems. The results have been obtained by developing a simulation program using MATLAB software @ R2018b. The parameters used for HHO are listed in Table 1. Whereas, the parameters used for other compared algorithms are mentioned in Appendix A, Table A1.

IEEE Nine Bus System for DOCRs
The proposed algorithm is implemented on IEEE 9 bus system as shown in Figure 6. The system consists of nine buses and 12 lines. The system is powered by a generator located at bus 1. The primary and backup pairs of DOCRs and short circuit test is mentioned in reference [13] and not discussed further here. The current transformer ration is set at 500/1 for all DOCRs. The higher and lower values of TDS (1.2-0.1) and PS (2.5-0.5) are selected accordingly with a coordination interval of 0.2 s. The optimum results obtained for TDS and PS by the proposed HHO are shown in Table 2. Table 3 shows the net total operating time achieved by the HHO and other state of the art algorithms and it cn be realized that the HHO outperforms the other algorithms in minimizing the total operating time of DOCR for an IEEE 9 Bus system. By comparison, it is evident that HHO is giving desired optimal settings for DOCRs as compared to other algorithms. For example, the improvement in operating time achieved by using HHO with respect to particle swarm optimization (PSO) is 37.55%. Meanwhile, the improvements with respect to genetic algorithm (GA), non-linear programming (NLP), informative differential evolution (IDE) algorithm, harmony search (HS), biogeography-based optimization (BBO) and modified adaptive teaching learningbased optimization algorithm (MTLBO) are 73.23%, 55.03%, 85.37%, 11.28%, 69.73% and 79.17% respectively. Figure 7 depicts the convergence characteristics of HHO obtained during the course of the simulation. It can be noticed that HHO achieved a quite fast convergence rate to achieve the optimal values. At about 150 iterations, the HHO was able to attain the optimal values.    Table 3. Comparison of HHO with other algorithms.

IEEE 15 Bus System for DOCRs
The IEEE 15 bus system is a highly distributed generator (DG)-enlarged distribution system consisting of 21 lines and 42 relays as mentioned in Figure 8 and is has 82 constraints and 84 design variables. The short circuit test and primary and back up configuration are mentioned in reference [24]. Table 4 shows the current transformer ratio for the DOCR configuration. The higher and lower limits of TDS (1.2-0.1) and PS (2.5-0.5) are selected accordingly. A CTI of 0.3 s is selected.

IEEE 15 Bus System for DOCRs
The IEEE 15 bus system is a highly distributed generator (DG)-enlarged distribution system consisting of 21 lines and 42 relays as mentioned in Figure 8 and is has 82 constraints and 84 design variables. The short circuit test and primary and back up configuration are mentioned in reference [24]. Table 4   The optimum values obtained by the proposed HHO are listed in Table 5 which shows that HHO optimized the total operating time up to attain the minimum and optimum values. Table 6 shows the comparison of HHO with other algorithms applied for the  The optimum values obtained by the proposed HHO are listed in Table 5 which shows that HHO optimized the total operating time up to attain the minimum and optimum values. Table 6 shows the comparison of HHO with other algorithms applied for the same DOCR coordination problem. It proves that HHO is superior to other current algorithms in minimizing the total operating time up to a minimum value with a fast convergence rate as shown in Figure 9 and obtains the best value for the objective function in a smaller number of iterations. The comparison of optimal settings determined by HHO with seeker algorithm (SA), mixed integer non-linear programming (MINLP), analytic approach (AA), differential evolution (DE), HS, backtracking search algorithm (BSA), MTLBO, group search optimization (GSO), improved group search optimization (IGSO) and modified electromagnetic field optimization (MEFO) shows 5.64%, 24  same DOCR coordination problem. It proves that HHO is superior to other current algorithms in minimizing the total operating time up to a minimum value with a fast convergence rate as shown in Figure 9 and

Algorithm
Objective Function SA [24] 12.227 MINLP [24] 15.335 AA [42] 11.6618 DE [44] 11.7591 HS [44] 12.6225 BSA [37] 16.293 MTLBO [41] 52.5039 GSO [43] 13.6542 IGSO [43] 12.135 MEFO [40] 13.953 Proposed HHO 11.537 Total Operating Time (seconds) Iteration Figure 10. Convergence characteristics of HHO for a 15-bus system.  Figure 11 shows the single line diagram of the IEEE 14-bus system that consists of 14 buses and 40 relays. The listed relays' CT ratios (CTRs), P/B relay pairs, current transformer ration other information regarding this test system are shown in Table 7. The fault currents for the close-in 3φ faults and system information are mentioned in reference [58].  Table 8, which shows that the proposed algorithm optimized and minimized all the values up to optimum values. Figure 12 shows the convergence characteristics of the objective function value obtained as a result of the simulation, which show that the convergence is faster and achieved the optimum value in a smaller number of iterations. For the 14-bus network, HHO was able to attain optimal results with a fast convergence rate after about 115 iterations. Table 9 shows the comparison of the proposed algorithm with other state of the art algorithms, which confirms the superiority of the proposed HHO. The overall improvement in optimal settings determined by HHO as compared to the hybrid genetic algorithm linear programming (HGA-LP), mixed integer linear programming (MILP), multiple embedded crossover PSO (MECPSO) and modified adaptive PSO (MAPSO) algorithms is 4.18%, 1.63%, 0.06%, 8.48%, respectively. that the convergence is faster and achieved the optimum value in a smaller number of iterations. For the 14-bus network, HHO was able to attain optimal results with a fast convergence rate after about 115 iterations. Table 9 shows the comparison of the proposed algorithm with other state of the art algorithms, which confirms the superiority of the proposed HHO. The overall improvement in optimal settings determined by HHO as compared to the hybrid genetic algorithm linear programming (HGA-LP), mixed integer linear programming (MILP), multiple embedded crossover PSO (MECPSO) and modified adaptive PSO (MAPSO) algorithms is 4.18%, 1.63%, 0.06%, 8.48%, respectively.    Figure 12. Convergence characteristics of HHO for a 14-bus system.

Conclusions
In this paper, the network relay coordination and optimization problem has been solved by using HHO. Relay optimization is formulated as MINLP targeted to minimize Figure 12. Convergence characteristics of HHO for a 14-bus system.  Table 9. Comparison of HHO with other algorithms.

Conclusions
In this paper, the network relay coordination and optimization problem has been solved by using HHO. Relay optimization is formulated as MINLP targeted to minimize the overall operating time of relays by selecting as design parameters TDS and PS. For evaluation, three test systems with different scenarios have been considered. For DOCR optimization and coordination, a 9-bus system with a single generator and a 15-bus system with multiple DG penetration were tested, while for numerical relays, a 14-bus system having conventional as well as DGs is tested. The unique sieging and hunting capability of the HHO has been found affective in finding the global optimum values with robustness and better convergence as compared to other state of the art algorithms. The algorithm-wise comparison, for all three test systems, shows improved and optimum settings are found by HHO. The obtained results justify or claim that HHO is successful at finding better and optimum solutions for DOCRs and numerical relays proving it to be an effective tool for relay coordination and optimization.

Data Availability Statement:
The data used to support the finding of this study are included within the article.

Conflicts of Interest:
The authors declare no conflict of interest.