An Effective Coordination Setting for Directional Overcurrent Relays Using Modiﬁed Harris Hawk Optimization

: The relay optimization expresses quite a challenge for smooth and optimal operation of power system networks. The relay optimization is formulated as a mixed integer non-linear problem and is highly constrained. Furthermore, a reliable relaying system must be able to detect and isolate the faulted portion in a timely manner. Therefore, it is necessary to ﬁnd optimal parameters for relay settings to be able to respond in a timely way to the encountered fault and at the same time keep in consideration the operational and coordination constraints. This paper proposes modiﬁed Harris hawk optimization (MHHO), which is based on the intelligent preying tactics of Harris hawks and the improvement of intended modiﬁcations, crowding distance and roulette wheel selection. The proposed algorithm has been tested on IEEE 8 and 15-bus systems, using MATLAB programming. The test systems are the distribution networks covering the medium level voltage for consideration. The simulation results veriﬁed the success of MHHO to ﬁnd optimal settings for the relays. For IEEE 8-bus system, MHHO was able to give 35.45% improvement in the results in comparison to other algorithms. Furthermore, for the IEEE 15-bus system, MHHO showed 24.09% improvement on average. The comparison of the results obtained by MHHO with the other state-of-the-art algorithms proved that it is the strong candidate for optimization of the relay coordination problem.


Introduction
Power system protection is critically important for uninterrupted supply to healthy portion systems and timely fault isolation. As unwanted operation of a protection system may become the reason for an imbalance in a power supply and financial loss [1]. Overcurrent relays (OCR) are widely used in protection schema as a main protection as well as backup protection. They are being used as primary protection in distribution networks while they are also used in transmission and distribution networks as backup devices. OCRs are best known for their reliability, sensitivity, accuracy and selectivity [2].
In the event of a fault, primary relays are responsible for isolation of faulty portion from healthy portion of the network. Backup relays are supposed to be functional after a time margin for operations if primary relays fail to achieve the tripping. This time margin in referred to as the coordination time interval 'CTI' [3]. To enhance the reliability and security concerns in an interconnected meshed network, incorporating distributed generation, OCRs are linked with directional elements, referred to as directional overcurrent relays (DOCR) [4]. DOCRs are dependent upon time multiplier settings (TMS) and plug settings (PS) or pickup current (Ip) for their fine tuning and optimization. The optimization of relay operation schema is carried out by ensuring proper coordination of relays while achieving selectivity as well as fast isolation while fulfilling the performance and coordination constraints [5,6].

•
Development of MHHO by modifications of crowding distance and roulette wheel selection, with aim to enhance the ability to converge towards better solution while maintaining diversity, to the Harris hawk optimization algorithm.

•
Accessing the performance of MHHO by implementing it to the standard IEEE 8-bus and 15-bus test systems.
• Verification of the good performance of proposed MHHO by finding optimal settings of TMS and PS, satisfying the allowed limits, to minimize the overall operating time of DOCRs.
The rest of the paper is organized as follows: Section 2 for problem formulation, Section 3 comprises of proposed algorithm and its usage guidelines for DOCR coordination problem. Section 4 shows the results and simulation part to prove the effectiveness of MHHO for DOCR coordination problem by comparison with other algorithms. Section 5 concludes the paper.

Objective Function
The coordination problem of directional overcurrent relays is stated as an optimization problem that is to be minimized. The objective function is stated as a sum of primary relays for a system. Based on the objective function, the optimal setting of relays is achieved by minimizing the operating times of primary relays while satisfying the coordination constraints: where, 'n' shows the number of primary relays. 'w i ' shows the assigned weight to the relays, which is assigned as 1 for all relays. 'T i ' is operating times of the primary relays (where i = 1, . . . , no. of relays). The operating times of the primary relays can be constituted on their characteristics curve, by using Equation (2) [43]: where, 'TMS i ' is the time multiplier settings for the relays. The values of α and β varies according to type of characteristics used for relays as stated in Table 1 [44]. 'I fi ' is the fault current value seen by that relay and 'I pi ' is the pickup current, given be Equation (3). 'PS i ' stands for plug settings and 'CT i ' indicates the current transformer ratio value. In this paper, standard inverse characteristics of overcurrent relay has been used. Figure 1 shows the DOCR optimization goals and constraints necessary to be taken care of. This figure states the relay constraints: coordination and operational. The coordination constraints are fulfilled when backup relays operate after a certain CTI, in case of a primary relays failure. In addition, operational constraints are fulfilled when relay parameters, i.e., TMS and PS, remain within the allowable limits of maximum and minimum values. Furthermore, for DOCR optimization TMS and PS have been considered as design variables. The further information of constraints is explained in following section:

Constraints
In the case of a failure of primary relays, a backup relay should operate to fulfill the protection criteria. This constraint can be regarded as the coordination constraint. Along with the coordination constraint, DOCR coordination is subject to the operational constraints also. The operational constraints are imposed on the design variables. Therefore, the DOCR relay optimization and coordination problem is subjected to following constraints: These constraints state that the TMS value of a relay should be between 0.1 and 1.1, which are the minimum and maximum values of TMS. Furthermore, PS should also remain in the bounds of minimum and maximum values of 0.5 to 2.5. According to Equation (6), the difference between operating times of primary relays is: Tpri and backup relays Tbackup should also obey the defined limit, 0.3 s. This is set so that no dis-coordination happens between primary and backup relays and is called a selectivity constraint or coordination constraint. The set value for CTI varies between 0.2-0.5 s and it depends upon the type of relay or choice made by the users. Furthermore, the primary operating time of a relay should also be bounded between the limits of 0.05 to 1.0 s.

Constraints Enforcement Function
To ensure that no relay go out of bounds and defined constraints, a penalty function is added to the objective function of the algorithm. For DOCR coordination, the penalty function is taken as a large value. In the case of violations, this value is added to an objective function [45].

Constraints
In the case of a failure of primary relays, a backup relay should operate to fulfill the protection criteria. This constraint can be regarded as the coordination constraint. Along with the coordination constraint, DOCR coordination is subject to the operational constraints also. The operational constraints are imposed on the design variables. Therefore, the DOCR relay optimization and coordination problem is subjected to following constraints: These constraints state that the TMS value of a relay should be between 0.1 and 1.1, which are the minimum and maximum values of TMS. Furthermore, PS should also remain in the bounds of minimum and maximum values of 0.5 to 2.5. According to Equation (6), the difference between operating times of primary relays is: T pri and backup relays T backup should also obey the defined limit, 0.3 s. This is set so that no dis-coordination happens between primary and backup relays and is called a selectivity constraint or coordination constraint. The set value for CTI varies between 0.2-0.5 s and it depends upon the type of relay or choice made by the users. Furthermore, the primary operating time of a relay should also be bounded between the limits of 0.05 to 1.0 s.

Constraints Enforcement Function
To ensure that no relay go out of bounds and defined constraints, a penalty function is added to the objective function of the algorithm. For DOCR coordination, the penalty function is taken as a large value. In the case of violations, this value is added to an objective function [45].
where the Penalty function values varies from 1 to 'k' entries, which indicates the relay pairs. In Equation (8), the Penalty function gives zero value when the condition shown by Equation (6) is fulfilled, and δ is the big value for solution with the violation of constraints.
The function gives zero value if bounds are obeyed. For optimal minimization, the penalty function value must be zero.

Modified Harris Hawk Optimization (MHHO)
Living things found in nature have been evolving and adapting themselves for survival so nature can be regarded as the best optimizer. Furthermore, Harris hawks change their foraging and besieging strategies as per the prey's behavior. Therefore, this unique and effective method of the hunt became the inspiration for Harris hawk optimization (HHO). HHO can be applied as a black box to the optimization problems provided the proper definition and formulation of objective function. Analogous to other nature-inspired algorithms, HHO contains the exploration and exploitation phases [46]. Nevertheless, the exploitation phases are modelled as per the prey behavior during a hunt for different scenarios, as explained in the following sections.

Exploration Mode
Generally, hawks can detect and track their prey in no time with their powerful eyes. But sometimes it is harder to find the prey and hawks have to perch at any place and then scan the area for hours. In HHO this behavior or process is referred as an exploration process. Hawks are designated as candidate solutions and anticipated prey as best candidate solution or nearly optimum. Hawks perch based on two strategies for searching the prey. For equal consideration of q, the hawks can perch either on a position close to other members (equated in Equation (9) for q < 0.5) or perch randomly on a random location inside the group's range (equated in Equation (9) for q ≥ 0.5).
where, for current iteration; X rand (t) is a hawk selected randomly, X(t) is the position vector of hawks, X rabbit (t) is the position of prey which is rabbit, X m (t) is average position of hawks for a current population, L.B. and U.B. are the lower and upper bounds respectively for the variables and r 1 , r 2 , r 3 and r 4 are randomly selected numbers between the interval (0, 1). Lastly, for next iteration, X(t + 1) is the position vector for hawks. As per this framework, random locations are generated in the specified range, as by lower and upper bounds. The above equation has two parts; the first one states the generation of solutions with respect to random location of hawks and the current position. The second part states the best solution found and difference of average position of hawks' team and at the same time r 3 and r 4 act as random scaling elements. This gives diversity to the search space for exploration of an optimal solution. Furthermore, the average position of hawks is equated as follows: where, N states the total number of hawks and X i (t) states the location of the individual hawk for the iteration t.

Variation from Exploration to Exploitation Mode
HHO shifts from exploration to exploitation and then as per the escaping energy of prey, which is rabbit, different exploitation modes are adopted. While the escaping energy of a rabbit decreases as time passes during the hunt. The escaping energy of a rabbit can be equated as follows: where, E states the escaping energy of the rabbit, E o states the starting state of energy (it changes randomly at each iteration within the interval of (−1, 1)), t is the current iteration and T indicates the maximum number of iterations.
When the value of E o is between −1 and 0, this indicates that the escaping energy of rabbit is towards decreasing side. The value from 0 to 1 shows the escaping energy value towards increasing side. The escaping energy of the rabbit decreases with the increasing iterations. The transition of HHO from exploration to exploitation depends upon the escaping energy of the rabbit. When |E| ≥ 1, HHO is in exploration mode, for |E| < 1, exploitation mode is initiated.

Exploitation Mode
When in exploitation mode, hawks make a surprise pounce on the targeted rabbit. Here, parameter r illustrates the chance of the prey escaping. The value r < 0.5 indicates the successful escape of the rabbit whereas r ≥ 0.5 indicates the unsuccessful escape attempt of the rabbit. Depending upon the energy of rabbit and escaping pattern hawks adopts the besiege strategy between hard and soft besieging. Eventually, hawks bring the rabbit towards exhaustion by repeated pounce attempts. Then hawks carry out a hard besiege to successfully catch the prey.
For HHO, both the parameters r and E are used for switching between soft and hard besiege modes. Depending upon the energy of rabbit and chance of escape from striking hawks, besiege and hunt process can be classified into four categories, which are explained as follows.

Soft Besiege Mode
For r ≥ 0.5 and |E| ≥ 0.5 rabbit tries jump randomly to escape the hunt as it has enough energy. In such a scenario, hawks cooperatively give it surprise pounces, keeping soft besiege, to exhaust the rabbit. The following equation models the rabbit behavior for such a scenario: where, ∆X(t) shows the difference of position vector of rabbit with respect to current location, J is the random jump attempt by a rabbit to deceive the hunt, r 5 is the random number between the interval (0, 1) to elaborate the random jumping nature; whereas, ∆X(t) is equated as follows:

Hard Besiege Mode
For r ≥ 0.5 and |E| < 0.5, the escaping energy of a rabbit is considerably decreased. In such a scenario, hawks besiege the rabbit fiercely and give it a strong surprise pounce to capture it. For HHO, the following equation states the hawks' hunting step for next iteration:

Soft Besiege Mode with Rapid Dives
For r < 0.5 and |E| ≥ 0.5, the rabbit has enough escaping energy so that it can deceive the hawks. In such situations, hawks still adopt the soft besiege mode for surprise pounces. To include a mathematical framework of hawks' leapfrog movements [47] and deceptive pattern of rabbit to escape, the Lévy flight (LF) approach is incorporated in HHO. The rabbit performs zigzag deceptive moves to deceive the hawks. But in response to this, hawks analyze their moves and also implement sporadic and quick dives. The LF-model best explains the searching strategy, adopted by foragers [48,49]. Furthermore, the LF-model is also being utilized by other animals like moneys and sharks during their chase [50][51][52][53].
As hawks are regarded as quite intelligent foragers, it is assumed that they can sense the situation and select the best possible dive to target the rabbit. If hawks decide to do the soft besiege of the rabbit, then the following equation demonstrates the next move performed by them.
Electronics 2021, 10, 3007 Subsequently, hawks rethink of their dives according to the deceptive moves of the rabbit. If rabbit is showing more deceptive movements hawks change their strategy and they perform irregular pattern dives, which will be based on LF.
where, S is a random vector of size 1 × D and LF represents the Lévy flight function, which is equated as follows: where, u and v are randomly selected values between the interval (0, 1) and β is the constant equals to 1.5. Eventually, the soft besiege mode adopted by hawks, in such scenario, can be modelled by the following equation:

Hard Besiege Mode with Rapid Dives
For r < 0.5 and |E| < 0.5, the escaping energy of a rabbit is decreased significantly. At this stage, hawks decrease their average distance from the rabbit and conduct hard besiege for the surprise pounce. Here, along with the hard besiege strategy, the LF model is also being used by hawks to successfully capture the rabbit. This behavior of hawks is modelled for HHO by the following equation: where, Y and Z are equated as follows: For Equation (21), X m (t) can be accessed from Equation (10). In this way depending upon the values of r and E, hawks change their strategy of preying and switch between different modes of besiege. Finally, after capture, groups of hawks enjoy the feast cooperatively. In addition, this is the plus point of hawks that they capture and hunt in such an intelligent way that their success is guaranteed.

Modifications
This section explains the modifications done to HHO to increase its ability to find the optimum results. The following two modifications were made to increase the diversity and optimal selection capability:

Crowding Distance Sorting
Crowding distance is the estimation of the density of the solutions with respect to the solution they surround [54]. The crowded tournament selection is based on ranking and distance. In other words, if a solution x i has a better rank than x j , we select x i . If the ranks are the same but D i > D j , x i is selected based on its crowding distance. By this estimation, along with the increased density, the chances of better solution are higher, which will help in the optimal solution. In Figure 2, the cuboid shown is the crowding distance of i.

Crowding Distance Sorting
Crowding distance is the estimation of the density of the solutions with respect to the solution they surround [54]. The crowded tournament selection is based on ranking and distance. In other words, if a solution xi has a better rank than xj, we select xi. If the ranks are the same but Di > Dj, xi is selected based on its crowding distance. By this estimation, along with the increased density, the chances of better solution are higher, which will help in the optimal solution. In Figure 2, the cuboid shown is the crowding distance of i. The crowding distance is calculated by sorting the solutions as per the objective function values in ascending order. The crowding distance value is the average distance between the particular solutions with respect to its two neighbors. In this process, the maximum and minimum objective function boundary solutions are assigned infinity values, this ensures that the boundary solutions are always selected. Other than boundary solutions, all solutions are allocated with the value corresponding to absolute normalized difference between two adjacent solutions. This procedure is carried out for each and every objective function value. The conclusive crowding distance is the sum of all crowding distance values for every objective function. For minimization problems, the maximum value for crowding distance is the worst value. The pseudocode of crowding distance-working schema is shown in Algorithm 1.

Algorithm 1 Pseudocode of Crowding Distance 1:
Input a solution set 2: Calculate the size of solution set as far 'nPop' and 'nVar' 3: Initialize the distance vector D 4: Give 'infinity' value to min and max boundary solutions 5: for (j = 1 → nPop) 6: for j = 2 → (nPop-1) 7: D = normalize (all positions (i)-all positions (j)) 8: End 9: Crowding Distance = sum (D) 10: End The crowding distance is calculated by sorting the solutions as per the objective function values in ascending order. The crowding distance value is the average distance between the particular solutions with respect to its two neighbors. In this process, the maximum and minimum objective function boundary solutions are assigned infinity values, this ensures that the boundary solutions are always selected. Other than boundary solutions, all solutions are allocated with the value corresponding to absolute normalized difference between two adjacent solutions. This procedure is carried out for each and every objective function value. The conclusive crowding distance is the sum of all crowding distance values for every objective function. For minimization problems, the maximum value for crowding distance is the worst value. The pseudocode of crowding distanceworking schema is shown in Algorithm 1.

Roulette Wheel Selection
This technique was adopted under the inspiration of wide usage in a genetic algorithm. In a population with n individuals, for each element x with a corresponding fitness value f x , the corresponding probability P x of selection can be computed as: This sorting of solutions by crowding distance and then selecting on the basis of roulette wheel, increases the chances of selection of best value; the best rabbit in this case, is regarded as the best solution. Therefore, as a result along with increasing the diversity of solution space, the selection chances of best solution is boosted by these modifications.

MHHO Algorithm Workflow
The following section explains the details of steps involved to carry out MHHO search for DOCR optimization. Furthermore, in Figure 3, a flowchart is also included to explain the working of MHHO.

Simulation and Results
The simulation was carried out on two test systems: the IEEE 8-Bus system [55] 15-Bus system [56]. The results were obtained using MATLAB R2020a. The simula parameters used for the algorithm are shown below in Table 2.  Figure 3. Flowchart for modified Harris hawk optimization (MHHO).

Workflow Steps and Flowchart
The steps of MHHO workflow are explained as follows: i. Define the objective function for relay coordination and operation, which is bound to obey the constraints and read the fault currents and current transformer ratio data. ii. Initialize the hawk search mechanism by defining the upper and lower bounds for design variables PS and TMS, MHHO parameters (search agents, population size and maximum iterations). iii. Start the iteration counter from 1. iv. For each hawk search, check the objective function status using Equation (7). v.
Designate the best objective function value as the best hawk. vi. For each hawk, update the value of energy, initial energy and random jump of rabbit by using Equations (11) and (13). vii. Check the value of energy, if E ≥ 1, assign the position to hawk according to Equation (9), according to perching probability q, or else jump to step (viii). viii. Consider E ≥ 0.5, jump to step (ix) or else to step (x).
If r ≥ 0.5, assign hawk position as per Equation (15) or else by Equation (20). Afterwards, go to step (xi). xi. Calculate the crowding distance and sort the population in ascending order. xii. Select the position that have the shortest crowding distance and select the remaining positions according to best position found by Harris Hawks. xiii. Select the rabbit position for the next iteration based on the roulette wheel, on the basis of Equation (22). xiv. Check if maximum iteration limit is reached, if yes then go to step (xv) otherwise increase the iteration counter and start over from step (iv). xv. Return the optimal parameters found for relay coordination and operation.

Simulation and Results
The simulation was carried out on two test systems: the IEEE 8-Bus system [55] and 15-Bus system [56]. The results were obtained using MATLAB R2020a. The simulation parameters used for the algorithm are shown below in Table 2.

IEEE 8-Bus Network
The network consists of eight buses, two generating units, one external grid connection, two transformers and 14 DOCRs. At Bus-4, the test system is connected to an external grid having short-circuit capacity of 400 MVA. Meanwhile, both the generator units have equal ratings of 150 MVA, 10 kV and reactance of x = 15%. Furthermore, transformer units have equal ratings of 150 MVA, 10 kV and reactance of x = 4%. The relay coordination scheme is tested for three-phase faults at different points. A single line diagram of the eight-bus test system is shown in Figure 4. While current transformer ratios, fault current values and relevant relays are shown in Table 3. The max iterations were taken as 400 and population size was set as 30. have equal ratings of 150 MVA, 10 kV and reactance of x = 4%. The relay coordination scheme is tested for three-phase faults at different points. A single line diagram of the eight-bus test system is shown in Figure 4. While current transformer ratios, fault current values and relevant relays are shown in Table 3. The max iterations were taken as 400 and population size was set as 30.   Figure 4. IEEE 8-Bus network. The optimized parameters for a relay coordination scheme are shown in Table 4. The results of MHHO are compared with baseline HHO [57] and other algorithms: the seeker algorithm (SA) [58], genetic algorithm (GA), genetic algorithm-linear programming (GA-LP) [59], biogeography-based optimization (BBO), biogeography-based optimization linear programming (BBO-LP) [56] and JAYA [45] algorithms. The objective function value of MHHO was compared with HHO and other algorithms. There is a noticeable improvement obtained by MHHO, as shown in Figure 5. Furthermore, Table 5 shows the percentage of improvement in results by MHHO as compared to algorithms mentioned above for the IEEE 8-bus system. The objective function value of MHHO was compared with HHO and other algorithms. There is a noticeable improvement obtained by MHHO, as shown in Figure 5. Furthermore, Table 5 shows the percentage of improvement in results by MHHO as compared to algorithms mentioned above for the IEEE 8-bus system.

IEEE 15-Bus Network
The test system comprising 15 buses and 21 branches is being protected by 42 DOCRs. The three-phase close-in faults are considered for testing the relay protection system. The system is equipped with high distributed generation penetration. All the six generators have ratings of 15 MVA, 20 kV and x = 15% synchronous reactance. Furthermore, all the lines have equal impedance value of Z = 0.19 + j0.46 Ω/km. At Bus-8, the system is connected with an external grid having short-circuit capacity of 200 MVA. A one-line diagram of 15-bus system is shown in Figure 6.

IEEE 15-Bus Network
The test system comprising 15 buses and 21 branches is being protected by 42 DOCRs. The three-phase close-in faults are considered for testing the relay protection system. The system is equipped with high distributed generation penetration. All the six generators have ratings of 15 MVA, 20 kV and x = 15% synchronous reactance. Furthermore, all the lines have equal impedance value of Z = 0.19 + j0.46 Ω/km. At Bus-8, the system is connected with an external grid having short-circuit capacity of 200 MVA. A one-line diagram of 15-bus system is shown in Figure 6. The current transformer ratio for the relays are shown in Table 6. Furthermore, primary and backup relays data of a test system with relevant fault current values can be accessed from the reference [60]. The population size is 60 and iteration limit is set as 800. The optimal results obtained by MHHO are shown in Table 7. The range for design variables TMS and PS varies between (0.01-1.1) and (0.5-2.5) respectively.  The current transformer ratio for the relays are shown in Table 6. Furthermore, primary and backup relays data of a test system with relevant fault current values can be accessed from the reference [60]. The population size is 60 and iteration limit is set as 800. The optimal results obtained by MHHO are shown in Table 7. The range for design variables TMS and PS varies between (0.01-1.1) and (0.5-2.5) respectively.  The comparison of MHHO with HHO shows the better and optimal results for MHHO. Figure 7 shows the objective function value comparison of MHHO with HHO and other state-of-the-art algorithms. The results of MHHO are compared with baseline HHO and other algorithms: SA [58], mixed integer non-linear programming (MINLP) [58], analytic approach (AA) [61], differential evolution (DE) [60], harmony search (HS) [60], backtracking search algorithm (BSA) [62], modified adaptive teaching learning-based optimization algorithm (MTLBO) [63], group search optimization (GSO) [64], improved group search optimization (IGSO) [64], and modified electromagnetic field optimization (MEFO) [65]. Table 8 shows the percentage of improvement in results by MHHO as compared to algorithms mentioned above for IEEE 15-bus system. HHO and other algorithms: SA [58], mixed integer non-linear programming (MINLP) [58], analytic approach (AA) [61], differential evolution (DE) [60], harmony search (HS) [60], backtracking search algorithm (BSA) [62], modified adaptive teaching learning-based optimization algorithm (MTLBO) [63], group search optimization (GSO) [64], improved group search optimization (IGSO) [64], and modified electromagnetic field optimization (MEFO) [65]. Table 8 shows the percentage of improvement in results by MHHO as compared to algorithms mentioned above for IEEE 15-bus system.

Conclusions
In this paper, MHHO is proposed to address the coordination and optimization problem of DOCRs. The modification of crowding distance and roulette wheel selection have been made to boost the optimization capability, which are able to increase the diversity of overall search space as well as increase the ability to find optimal minimized settings for the DOCR problem. Owing to MINLP, two decision variables are considered to find optimal minimized settings for DOCRs; TMS and PS. MHHO is able to find more improved results as compared to baseline HHO, owing to the modifications. MHHO is able to achieve 16.73% improved performance for the IEEE 8-bus network, as compared to HHO. For the 15-bus network, MHHO has achieved 9.0491% better results, as compared to HHO. The better results shown by MHHO have proved its superiority over HHO.
To further investigate the superiority of MHHO, it has been compared to various stateof-the-art algorithms. For both of the test systems, MHHO is able to show better results. For the 8-bus network, MHHO is able to exhibit 35.45% superior results. Furthermore, for 15-bus network, MHHO has demonstrated 24.09% improvement. The better results of MHHO indicates its superiority and reliability as a valuable optimization tool. As a result, the modifications of crowding distance and roulette wheel selection could be made to various population-based metaheuristic algorithms like PSO, grey wolf optimization, whale optimization etc., to check for better performance and promising results for future research.
In future work, the MHHO will be applied in optimization of relays system used in bigger and complex test systems like 30-bus system and 42-bus systems. Furthermore, MHHO will be considered for applications to the more diversified power system environments like AC microgrids.