A Multi-Objective Optimal Power Flow Control of Electrical Transmission Networks Using Intelligent Meta-Heuristic Optimization Techniques

: Optimal power ﬂow (OPF) is considered one of the most critical challenges that can substantially impact the sustainable performance of power systems. Solving the OPF problem reduces three essential items: operation costs, transmission losses, and voltage drops. An intelligent controller is needed to adjust the power system’s control parameters to solve this problem optimally. However, many constraints must be considered that make the design process of the OPF algorithm exceedingly tricky due to the increased number of limitations and control variables. This paper proposes a multi-objective intelligent control technique based on three different meta-heuristic optimization algorithms: multi-verse optimization (MVO), grasshopper optimization (GOA), and Harris hawks optimization (HHO) to solve the OPF problem. The proposed control techniques were validated by applying them to the IEEE-30 bus system under different operating conditions through MATLAB simulations. The proposed techniques were then compared with the particle swarm optimization (PSO) algorithm, which is very popular in the literature studying how to solving the OPF problem. The obtained results show that the proposed methods are more effective in solving the OPF problem when compared to the commonly used PSO algorithm. The proposed HHO, in particular, shows that it can form a reliable candidate in solving power systems’ optimization problems.


Introduction
Electric power networks are very sophisticated and complex systems. This complexity is simply due to the fact that electric power systems consist of three sub-systems, each containing different elements. The first sub-system is the generating stations which mainly consist of generators and their additional items. The second sub-system is the transmission networks which contain transformer substations, transmission lines, and reactive power compensation units. The third sub-system is the distribution networks which are responsible for delivering electric power to the consumers. All these elements have to be operated in a reliable as well as economical manner. Each element of the power system acts as a control variable, as shown in Table 1.
Although all control variables mentioned in Table 1 can be set individually, they influence the power flows of the network. It is essential to mention the power balance between the generation and load consumption when using different setting combinations for the control variables of the power settings. However, the solution of the power flow equation is not enough for the proper operation of the power system as few setting combinations of the power system can achieve the economical and reliable operation of the power system. The optimal power flow (OPF) is the process responsible for selecting the setting combination of control variables of the power system to achieve optimal operation. OPF is considered a development of the optimal dispatch concept [1,2]. The objective function of the OPF aims to reduce either the operation cost or system losses and voltage drops through the adjustments of the control variables of the power system. Sometimes, the goal of the objective function is to reduce the three parameters at the same time. This means that the OPF objective function is a non-linear function with many constraints, which makes it not easy to solve. Due to this fact, traditional methods like Newton methods may never find optimum global variables. Additionally, in some cases where the objective function is not available in the algebraic form, the traditional methods fail in solving this type of function [3][4][5][6][7][8][9]. Based on this fact, metaheuristic optimization algorithms are highly recommended for solving the OPF problem since they can have the high possibility of finding the globally optimum values for objective functions that have no pre-defined algebraic form.
One of the famous metaheuristic optimization techniques that have been used in the literature in solving the OPF problem is particle swarm optimization (PSO). The PSO, in its basic form, faces some challenges in finding the globally optimum values for the OPF as the initial random allocation of the particles has a significant impact on the success rate of the algorithm [10][11][12][13]. A modified PSO algorithm is proposed in [10] to solve this problem by updating the position of the particles from other individuals.
Many researchers have also proposed genetic algorithms (GAs) for solving the OPF problems, as they have better accuracy in reaching the global optimum solutions [14][15][16]. However, the GA suffers from a severe problem which is the high computational burden.
In this paper, different new meta-heuristic optimization algorithms are presented to solve the OPF problem. The objective function of the OPF is designed so that it has a multi-objective instead of having a single objective. The objectives intended in this paper are reducing fuel costs, transmission power losses, and voltage drops. The proposed methods are tested using MATLAB simulations over the IEEE 30-bus test system. The performance of the proposed methods was assessed and evaluated in terms of accuracy and computational burden. The rest of this paper is organized as follows: Section 2 shows the mathematical modeling of the OPF objective function, including its relative constraints. Section 3 presents the proposed methods for solving the OPF problem. Four case studies are illustrated in Section 4, including a detailed analysis.

The Mathematical Modeling of OPF Objective Function
The OPF equation consists of two main vectors: the first is the input vector, which contains the parameters that are tuned so that OPF is achieved. The second vector is the output, formed from all parameters dependent on the input vector. The input vector is expressed as where NG is the number of generators, NT is the number of transformers, and NC is the number of reactive power compensation units. P Gi is the setting of the output power of the generator connected to bus i. T i is the transformer tap setting of the transformer connected to bus i. Q Ci is the reactive power value of the reactive power compensation unit connected to bus i. It is essential to mention that in the case of the output power of the generator, P G1 is the output power of the slack bus, so it is not included in the input vector. Instead, since it is dependent on the input vector. This is why the number of the bus of the output power in the input vector starts from bus 2. Regarding the input vector, it is formed according to Equation (2): where NL is the number of load buses, NTL is the number of transmission lines, V Li is the load voltage of bus i. Q Gi is the reactive power of the generator connected to bus i and S li is the load flowing in transmission line i. To understand the power flow equations, firstly, it is required to build the admittance matrix as per Equation (3): where I i is the current flowing through bus I, V i is the voltage of bus i and Y ij is the admittance element of the i th row and j th column. According to the load flow equations, the total generated power must equal the total power of the load and the power losses as shown in Equations (4) and (5): where P Di and Q Di are the active and reactive power demand of the load connected to bus i, θ ij is the voltage angle between buses i and j. G ij and B ij represent the real and imaginary parts of the admittance matrix, respectively. Although the power flow equations shown in Equations (4) and (5) look simple, the optimization of these equations is sophisticated due to the increased number of constraints due to the operation limits for the bus voltages and active and reactive power limits. The minimum and maximum limits for control variables in Equation (1) and the corresponding dependent variables in Equation (2) can be expressed as follows: Now, the OPF equation can be formed as per Equation (13): where λ p , λ v , λ Q and λ s are the penalty factors applied to keep the active power of the slack bus, bus load voltages, reactive power generation, and branch complex powers within acceptable limits. P lim Gl , V lim Li , Q lim Gi and S lim li are the limit values for the active power of the slack bus, bus load voltages, reactive power generation, and branch complex powers. Equations (14)-(17) represent the mathematical formulation of these limits: As mentioned before, there are three main objectives for the optimization of OPF, as shown in the following list: 1. The reduction in fuel costs; 2. The reduction in the power transmission losses; 3. The reduction in voltage drops.
To minimize the consumption of generator fuel, the following objective function is applied: where a i , b i and c i are the cost factors of the generator cost function.
In the case of reducing the power losses, the objective function can be expressed as where g ij is the conductance value of the transmission line between buses i and j. V i and θ i are the voltage magnitude and angle of bus i while V j and θ j are the voltage magnitude and angle of bus j. If the purpose of the OPF optimization is to reduce the voltage drop at the transmission lines, then the objective function is set as per Equation (20): (20) where V re f i is the reference value of voltage magnitude at bus i. Since the three objectives mentioned above are all critical, combining two or more objectives in the same objective function is possible. In this paper, a multi-objective function for achieving the three objectives is used as per Equation (21):

The Proposed Meta-Heuristic Algorithms for Solving the OPF Problem
As mentioned before, in this paper, three meta-heuristic optimization algorithms are used to optimize the solution of the OPF problem, as shown in the following subsections.

The Multi-Verse Optimization Algorithm
There are three main elements behind the multi-verse theory, which create the very first idea of the MVO. The first element was unobservable until now, which is the white hole; it occurs at the creation of a universe or the collision of two neighboring universes. The following element is the opposite of the first one in its behavior and is called the black hole. We can always observe black holes, and they are characterized by the vast gravitational forces that make every ambient object attracted to them. The last element is called the wormhole; it has the authority to exchange objects either between different universes or between different parts of one universe [38,39].
Additionally, the idea of the universe expansion process is clarified by the multi-verse theory, which depends mainly on the inflation rate. The universe elements are formed, and they are controlled by that rate. We can achieve a stable phase between two parallel universes by separating the three elements, white, black, and wormholes. This process is exactly like the MVO search process.
The above method has been applied in significant optimization applications and the management of processes related to renewable energy and power systems [38,39]. Figure 1 shows the main idea of the MVO, we have n universes, and each of them reflects a solution. The wormhole is the path of the objects from a high inflation rate universe to other lower inflation rate universes. The optimum case is when the universe can receive objects from all other universes, which means that it has the lowest inflation rate.
Since the three objectives mentioned above are all critical, comb objectives in the same objective function is possible. In this paper, a m tion for achieving the three objectives is used as per Equation (21):

The Proposed Meta-Heuristic Algorithms for Solving the OPF Pro
As mentioned before, in this paper, three meta-heuristic optimiza used to optimize the solution of the OPF problem, as shown in the foll

The Multi-Verse Optimization Algorithm
There are three main elements behind the multi-verse theory, wh first idea of the MVO. The first element was unobservable until now, hole; it occurs at the creation of a universe or the collision of two neig The following element is the opposite of the first one in its behavior and hole. We can always observe black holes, and they are characterized b tional forces that make every ambient object attracted to them. The la the wormhole; it has the authority to exchange objects either between or between different parts of one universe [38,39].
Additionally, the idea of the universe expansion process is clar verse theory, which depends mainly on the inflation rate. The uni formed, and they are controlled by that rate. We can achieve a stable p parallel universes by separating the three elements, white, black, an process is exactly like the MVO search process.
The above method has been applied in significant optimization a management of processes related to renewable energy and power syst Figure 1 shows the main idea of the MVO, we have universes reflects a solution. The wormhole is the path of the objects from a high verse to other lower inflation rate universes. The optimum case is whe receive objects from all other universes, which means that it has the lo To mathematically model the MVO algorithm, in Equation (22), tion for the roulette wheel mechanisms that can randomly arrange all u  To mathematically model the MVO algorithm, in Equation (22), there is an illustration for the roulette wheel mechanisms that can randomly arrange all universes, assuming that d is the number of variables and n is the number of candidate solutions [40][41][42]: . .
where x j i is the j th parameter of the i th universe. Each parameter can be calculated from Equation (23): where x j k is the j th parameter of the k th universe, r 1 is a random binary number which can be either 0 or 1, and N I(U i ) is the inflation rate of the i th universe. As shown in Equation (23), white holes are formed with different inflation rates using the roulette wheel mechanism. As we said, lower inflation universes can receive more objects through white/black holes. However, another mechanism describes more improvements in the inflation rate made by each universe using wormholes. This process can be shown as per Equation (24): where X j is the j th element of the best solution (best universe created). lb j and ub j are the lower and upper bounds of this element. r 2 , r 3 and r 4 are binary numbers. WEP is the wormhole existence probability while the TDR is traveling wave distance. As we go for more iterations, we observe the linearity of the increasing WEP which confirms the progress of the optimization algorithm and to what extent it can achieve the best solution.
The WEP is updated based on the adaptive equation as described in Equation (25): where min and max are the boundaries for the WEP coefficient; l is the order of the iteration, while L is the maximum number of iterations. As per Equation (12), this adaptive formula is used to update TDR, which acts similarly to the WEP; when the number of iterations increases, the value of the TDR increases, which guarantees a more precise local search in the path to find the best solution: where p is the coefficient that controls the accuracy and speed of algorithm convergence, Figure 2 shows the flow chart of the MVO algorithm.

The Grasshopper Optimization Algorithm
Grasshoppers are harmful insects. They are known for their harmful effect of reducing agriculture production. Figure 3 shows the change that occurs when grasshoppers travel and join a big group, among other creatures, despite these usually being seen individually [40,41]. The group size is big enough to be terrifying to ranchers. This behavior is seen in both the nymph and the adulthood, which makes them extraordinary [40,41].

The Grasshopper Optimization Algorithm
Grasshoppers are harmful insects. They are known for their harmful effect of reducing agriculture production. Figure 3 shows the change that occurs when grasshoppers travel and join a big group, among other creatures, despite these usually being seen individually [40,41]. The group size is big enough to be terrifying to ranchers. This behavior is seen in both the nymph and the adulthood, which makes them extraordinary [40,41]. Nymph grasshoppers gather in massive numbers, bouncing and moving like trigger clambers, eating their way through the harvest. After this stage, when they have grown into adults, they build a multitude of structures which are noticeable all around. This phase is the evacuation process of the grasshoppers. Nymph grasshoppers gather in massive numbers, bouncing and moving like trigger clambers, eating their way through the harvest. After this stage, when they have grown into adults, they build a multitude of structures which are noticeable all around. This phase is the evacuation process of the grasshoppers. In the icteric stage, the main characteristic of the herd is the moderate activity and little strides of the grasshopper. On the other hand, investigation for a food source and sudden movement among extended zones is the herd's most crucial ability. The natural inspired techniques categorize two denominations. The first denomination is exploration; in this denomination, the observation of the movement of the search agents is an unexpected motion. Subsequently, we need a mathematical model for the swarm behavior to accomplish the design of the developed inspired technique. The mathematical model of the swarming behavior of the grasshopper is presented as follows [40,41]:

Eggs
where is the current placement of any of the grasshoppers, is the social interaction, is the gravity force, and shows the wind advection. An imprecise attitude of the herd can be evaluated from the following equation: where 1 , 2 , and 3 are random numbers from [0,1]. The social interaction ( ) is obtained through the following equation: where ( ) is the displacement between two nearby grasshoppers, and it is evaluated from = | − |; ( ) is the function that indicates the vigor of social forces and is calculated as follows: Additionally, the (^) is a vector in which its magnitude equals one displacement and two grasshoppers and be obtained as follows: In the icteric stage, the main characteristic of the herd is the moderate activity and little strides of the grasshopper. On the other hand, investigation for a food source and sudden movement among extended zones is the herd's most crucial ability. The natural inspired techniques categorize two denominations. The first denomination is exploration; in this denomination, the observation of the movement of the search agents is an unexpected motion. Subsequently, we need a mathematical model for the swarm behavior to accomplish the design of the developed inspired technique. The mathematical model of the swarming behavior of the grasshopper is presented as follows [40,41]: where X n is the current placement of any of the grasshoppers, S n is the social interaction, G n is the gravity force, and A n shows the wind advection. An imprecise attitude of the herd can be evaluated from the following equation: where r 1 , r 2 , and r 3 are random numbers from [0,1]. The social interaction (S n ) is obtained through the following equation: where (dnm) is the displacement between two nearby grasshoppers, and it is evaluated from dnm = |X m − X n |; (s) is the function that indicates the vigor of social forces and is calculated as follows: Additionally, the (dˆm n ) is a vector in which its magnitude equals one displacement and two grasshoppers and be obtained as follows: where (f) represents the attraction vigor while (l) is the level of the attractive longitude. Figure 4 illustrates the effect of the value of (s) on the grasshopper's social interaction.
Additionally, it is shown in Figure 4a that the repulsion between grasshoppers befalls an interval of 0-2.079 with displacement changing from 0 to 15. A comfort zone occurs when the distance between two grasshoppers is equal to 2.079 units, making the attraction and repulsion between them vanish. In Figure 5, the variation of the two factors (l, f) is considered in plotting the value of (s) in Equation (30). It is noticeable that for a few values of l and f, for example, 1, 0.5, respectively, both the attraction and repulsion zones are microscopic. In Figure 6, the relationship between grasshoppers' interaction and the comfort zone is expressed by the value of (s) [40,41]. There is a prominent issue in implementing the value of (s) when applying strong forces between individual grasshoppers. Despite being able to determine the attraction and repulsion zones, its value reaches zero with an extended displacement of more than 10. The force of gravity is evaluated from the following: where (g) is the gravitational constant and (eˆg) shows a unity factor towards the center of the earth. The (A) component in Equation (27) is calculated as follows: where (u) is a constant drift and (eˆw) is a unit vector in the direction of the wind. As illustrated in Figure 3, the wind direction has a significant effect on the nymph grasshoppers as they have no wings. By replacing S, G, and A in Equation (27), the equation can be expressed as follows: (34) where N is the number of grasshoppers.
inability 2021, 13, x FOR PEER REVIEW 9 of where (f) represents the attraction vigor while (l) is the level of the attractive longitud Figure 4 illustrates the effect of the value of (s) on the grasshopper's social interactio Additionally, it is shown in Figure 4a that the repulsion between grasshoppers befalls interval of 0-2.079 with displacement changing from 0 to 15. A comfort zone occurs wh the distance between two grasshoppers is equal to 2.079 units, making the attraction a repulsion between them vanish. In Figure 5, the variation of the two factors (l, f) is co sidered in plotting the value of (s) in Equation (30). It is noticeable that for a few values l and f, for example, 1, 0.5, respectively, both the attraction and repulsion zones are micr scopic. In Figure 6, the relationship between grasshoppers' interaction and the comfo zone is expressed by the value of (s) [40,41]. There is a prominent issue in implementi the value of (s) when applying strong forces between individual grasshoppers. Desp being able to determine the attraction and repulsion zones, its value reaches zero with extended displacement of more than 10. The force of gravity is evaluated from the follo ing: is the gravitational constant and (^) shows a unity factor towards the center the earth. The (A) component in Equation (27) is calculated as follows: where (u) is a constant drift and (^) is a unit vector in the direction of the wind. As illu trated in Figure 3, the wind direction has a significant effect on the nymph grasshoppe as they have no wings. By replacing S, G, and A in Equation (27), the equation can expressed as follows: where N is the number of grasshoppers.
(a) (b)    At the point at which there is a change in the program's execution, the minor female of the grasshoppers is not allowed to reach that location, although these grasshoppers can land on the ground. As a result, in this case, the equation of the entire simulation blocked the algorithm from both the exploration and exploitation of a search agent around the solution, so it is never used. In conclusion, the implemented scheme of the herd took place in the free space. According to Equation (34), the interaction between the grasshoppers and each of the others in the swarm is implemented. However, when it comes to solving optimization problems, this mathematical model cannot be directly used because the grasshoppers reaches the comfort zone quickly, and the herd does not focus on a particular nearby point. Therefore, Equation (34) is modified and proposed as follows to solve the optimization problems: where and are the upper and lower bounds in the D-th dimensions, respectively. ^ is the optimal solution so far, and (c) is a decreasing degree to contract the attraction zone, repulsion zone, and comfort zone. The (G) component is neglected, and the wind direction is assumed to be always towards the target ( ^) . As it appears in Equation  At the point at which there is a change in the program's execution, the minor fe of the grasshoppers is not allowed to reach that location, although these grasshoppe land on the ground. As a result, in this case, the equation of the entire simulation blo the algorithm from both the exploration and exploitation of a search agent aroun solution, so it is never used. In conclusion, the implemented scheme of the herd took in the free space. According to Equation (34), the interaction between the grassho and each of the others in the swarm is implemented. However, when it comes to so optimization problems, this mathematical model cannot be directly used becaus grasshoppers reaches the comfort zone quickly, and the herd does not focus on a pa lar nearby point. Therefore, Equation (34) is modified and proposed as follows to the optimization problems: where and are the upper and lower bounds in the D-th dimensions, re tively. ^ is the optimal solution so far, and (c) is a decreasing degree to contract t traction zone, repulsion zone, and comfort zone. The (G) component is neglected, an wind direction is assumed to be always towards the target (^). As it appears in Equ (35), the following displacement of a given grasshopper is determined by three fa At the point at which there is a change in the program's execution, the minor female of the grasshoppers is not allowed to reach that location, although these grasshoppers can land on the ground. As a result, in this case, the equation of the entire simulation blocked the algorithm from both the exploration and exploitation of a search agent around the solution, so it is never used. In conclusion, the implemented scheme of the herd took place in the free space. According to Equation (34), the interaction between the grasshoppers and each of the others in the swarm is implemented. However, when it comes to solving optimization problems, this mathematical model cannot be directly used because the grasshoppers reaches the comfort zone quickly, and the herd does not focus on a particular nearby point. Therefore, Equation (34) is modified and proposed as follows to solve the optimization problems: where ub d and lb d are the upper and lower bounds in the D-th dimensions, respectively. Td is the optimal solution so far, and (c) is a decreasing degree to contract the attraction zone, repulsion zone, and comfort zone. The (G) component is neglected, and the wind direction is assumed to be always towards the target (Td). As it appears in Equation (35), the following displacement of a given grasshopper is determined by three factors. These factors are the current displacement of that given grasshopper, its objective displacement, and the locations of the other grasshoppers. This technique is different from that of PSO, as we mentioned before in the literature. The position and the velocity vectors are two critical factors that are needed to define each particle in the particle swarm optimization (PSO), while there is only one vectorthat is required to define the search agent in the grasshopper optimization (GOA).
Other factors differentiate between both techniques in determining the displacement of particles. According to the PSO, the essential factors in locating the position of particles are the current displacement, the best solution obtained by an individual, and the best solution obtained by the swarm. At the same time, concerning the GOA, it is the current location, the superior solution gained by an individual, the best solution obtained by the swarm, and the locations of the other search agents. According to Equation (35), it is clear that the adaptive element (c) has repeated two times for the following reasons: 1.
The first (c) is nearly comparable to the inertial weight (w) in PSO, and it is responsible for the remission of the motion of grasshoppers towards its target, which occurs by managing both exploration and exploitation.

2.
The second (c) aims to reduce the attraction, repulsion, and comfort zones between the grasshoppers.
With respect to Equation (35), it is evident that the element (c) inside the equation is directly proportional to the number of iterations as it participates in reducing the attraction and repulsion between the grasshoppers. Additionally, the outer element (c) plays a role in decreasing the concourse towards the target by increasing the number of iterations.
Finally, according to Equation (35), the start of this equation represents the location of the other grasshoppers and simulates their interaction in nature. Additionally, the second part which is identified by (Td) simulates its motion capability towards the target.
Generally, as in the icteric phase, when grasshoppers have no wings, they tend to stir and look locally for their food; in the next stage, they learn to move freely in the air as they explore much larger level zones.
In stochastic optimization techniques, finding promising regions of the search space is essential, and thus exploration is essential. After finding these regions, exploitation search agents search locally to find the global optimum as an accurate approximation value. The coefficient (c) is calculated as follows: where cmax and cmin are the limitation values, (i) indicates the present iteration, while (I) is the ultimate number of iterations. In this work, cmax and cmin are set to be 1 and 0.00004, respectively. In reality, the global optimum solution is unknown, so there is no target to achieve it. Therefore, there must be a clear objective for grasshoppers in each step as to which is the best objective value. This will help GOA to keep the most objective value in the search space in each iteration and require grasshoppers to move towards it. The flowchart of the grasshopper optimization technique is expressed in Figure 7. The GOA starts the optimization by initializing the behavior parameters such as S n , G n , A n , cmax, cmin, etc., then generating random solutions. Additionally, the fitness function is evaluated, leading to updating the locations of search agents based on Equation (33). The best target position was obtained and updated in each iteration. After that, the number of iterations is compared with the population size, and if the number of iterations is greater than the population size, then the best position will be observed if it reaches the best position. If not, the fitness function will be re-evaluated. Therefore, if the best position is achieved, it will be assigned to the senior position, and if not, the fitness function will be evaluated. Finally, the location and the objective of the best target are returned as the best approximation for the global optimum.

The Harris Hawks Optimization Algorithm
Harris hawks optimization is a metaheuristic optimization proposed in [42] that is inspired by the cooperative behavior of Harris hawks in hunting, chasing, and besieging their victims. The HHO is based on population optimization without having any gradients, which gives it a competitive edge over other techniques in terms of conversion speed.
The HHO consists of two main phases: exploration and exploitation. Additionally, there is a transition phase through which the algorithm is switched from exploration to exploitation.
In the exploration phase, Harris hawks start to search randomly for victims as per the following equation: where ( + 1) is the location of the hawks in the iteration ( + 1); ( ) is the location of the rabbit (the victim); 1 to 4 and are random numbers that can vary between 0 and 1; ( ) represents a hawk which is chosen randomly; and is the average location of the current population of hawks which can be calculated from Equation (38): where ( ) indicates the position of each hawk at iteration t while N represents the total number of hawks. As mentioned above, after finishing the exploration stage, there is a transient stage before moving to the exploitation stage. At this transient stage, it is necessary to model the energy of the rabbit as per Equation (39):

The Harris Hawks Optimization Algorithm
Harris hawks optimization is a metaheuristic optimization proposed in [42] that is inspired by the cooperative behavior of Harris hawks in hunting, chasing, and besieging their victims. The HHO is based on population optimization without having any gradients, which gives it a competitive edge over other techniques in terms of conversion speed.
The HHO consists of two main phases: exploration and exploitation. Additionally, there is a transition phase through which the algorithm is switched from exploration to exploitation.
In the exploration phase, Harris hawks start to search randomly for victims as per the following equation: where X(t + 1) is the location of the hawks in the iteration (t + 1); X rabbit (t) is the location of the rabbit (the victim); r 1 to r 4 and q are random numbers that can vary between 0 and 1; X rand (t) represents a hawk which is chosen randomly; and X m is the average location of the current population of hawks which can be calculated from Equation (38): where X i (t) indicates the position of each hawk at iteration t while N represents the total number of hawks. As mentioned above, after finishing the exploration stage, there is a transient stage before moving to the exploitation stage. At this transient stage, it is necessary to model the energy of the rabbit as per Equation (39): where E is the escaping energy of the rabbit, T is the maximum number of iterations and E 0 is the initial state of the rabbit energy. The value of E 0 is varies between −1 and 1 based on the physical fitness of the victim. When E 0 goes towards −1, this means that the victim is losing its energy and vice versa.
According to the behavior of rabbits, the relation between the rabbit energy and the time is inversely proportional. This means that as long as t increases, the E is decreased. Additionally, based on E, Harris hawks decide to either search different areas to detect the location of the rabbit when |E| ≥ 1, or move forward to the exploitation phase.
In the exploitation phase, two behaviors need to be modeled. The first is the soft besiege in which the rabbit energy is still high and can run fast; in this condition, Harris hawks try to softly follow and put it under surveillance until it starts to get exhausted. The second is the hard besiege; the prey in this behavior is tired and does not have sufficient energy to escape. As a result, the Harris hawks in this mode form closed circles to make a sudden attack. Figure 8 shows Harris hawk attack patterns.
where E is the escaping energy of the rabbit, T is the maximum number of iterations and 0 is the initial state of the rabbit energy. The value of 0 is varies between −1 and 1 based on the physical fitness of the victim. When 0 goes towards −1, this means that the victim is losing its energy and vice versa.
According to the behavior of rabbits, the relation between the rabbit energy and the time is inversely proportional. This means that as long as increases, the is decreased. Additionally, based on , Harris hawks decide to either search different areas to detect the location of the rabbit when | | ≥ 1, or move forward to the exploitation phase.
In the exploitation phase, two behaviors need to be modeled. The first is the soft besiege in which the rabbit energy is still high and can run fast; in this condition, Harris hawks try to softly follow and put it under surveillance until it starts to get exhausted. The second is the hard besiege; the prey in this behavior is tired and does not have sufficient energy to escape. As a result, the Harris hawks in this mode form closed circles to make a sudden attack. Figure 8 shows Harris hawk attack patterns.  Figure 8. Harris hawk prey attack patterns.

S o f t B e s i e g e
To mathematically model the two behaviors, let be the percentage of the successful escape of the rabbit. If | | ≥ 0.5 and ≥ 0.5, this means that the rabbit has relatively high escaping energy, and at the same time, the chance of successful escape is higher than 50%. This means that the Harris hawks will perform a soft besiege and will update their location according to Equation (40): where ∆ ( ) is the position difference between the rabbit and the hawks-this value can be calculated as follows: To mathematically model the two behaviors, let r be the percentage of the successful escape of the rabbit. If |E| ≥ 0.5 and r ≥ 0.5, this means that the rabbit has relatively high escaping energy, and at the same time, the chance of successful escape is higher than 50%. This means that the Harris hawks will perform a soft besiege and will update their location according to Equation (40): where ∆X(t) is the position difference between the rabbit and the hawks-this value can be calculated as follows: Moreover, J is a random number that represents the jump strength that can get from Equation (42) as follows: where r 5 is a random number that varies between 0 and 1. If |E| ≥ 0.5 and r < 0.5, this means that the rabbit has high energy. However, the chance of successful escape is not significant. In this case, the harris hawks will perform a soft besiege but with progressive and rapid dives. The next movement of the hawks will be updated according to: The hawks then will compare the current position with the previous dive to evaluate which is better. If the previous dive is better, the hawks will use it. If not, the hawks will then apply a new dive using the levy flight (LF) equation: where D is the problem dimension, S is a random vector with a size of 1 × D. The LF function can be calculated according to Equation (45): where u and v are a random number that varies between 0 and 1. β is a constant value of 1.5. σ which is calculated using: The Hariss hawks will then evaluate positions Y and Z and select the best position based on Equation (47): If |E| < 0.5 and r ≥ 0.5, this means that the rabbit has relatively low energy, but it has a moderate chance of successful escape. In this condition, the hawks will perform a hard besiege and will update their equation based on Equation (48): If |E| < 0.5 and r < 0.5, this means that the victim has low energy and also has a low chance to escape. In this situation, the hawks will also perform a hard besiege but with progressive rapid dives at which the next position of the hawks will be updated using Equation (21). Z will be calculated from Equation (18), and Y will be calculated using Equation (49) as follows: Figure 9 shows a flowchart of the proposed HHO.

Study Cases and Results Analysis
The IEEE 30 bus test system shown in Figure 10 is used for evaluating the proposed OPF algorithms. This system consists of six generating units, four transformers with tap changing units, and nine reactive power compensation units-in addition to 41 transmission lines. The total active power demand of the network is 283.4 MW, while the total reactive power demand is 126.2 MVAR. Figure 10 shows the single line diagram of the IEEE 30 bus system, while Tables 2-4 show the data of its components.  In this section, four case studies were presented to assess the performance during different operating modes. In the first case, the objective of the OPF was to reduce the cost, while the second cases show the performance of the OPF algorithm when the objective was minimizing the transmission losses. Case 3 presents the behavior of the OPF proposed algorithm when it aims to reduce the voltage drop of transmission lines. The last case simulates the action of the OPF proposed algorithm when it has three objectives that minimize the cost, power losses, and voltage drops. The performance of the three proposed algorithms is compared with the particle swarm optimization (PSO) algorithm since it is a famous technique in solving the OPF problem. All simulations are performed using MATLAB.

Case 1: Reducing Operating Costs
In this case, the proposed algorithms aim to reduce the operating cost by minimizing the consumption of generator fuel, which means that Equation (18) is selected as the objective function for the OPF algorithm. Table 5 shows the optimal results of the three proposed algorithms. All the obtained results are compared with those received from the PSO algorithms. The results shown in Table 5 clearly show the effectiveness of all proposed techniques compared with the PSO. The MVO algorithm has succeeded in obtaining lower costs than the PSO by 2.1%, while the GOA has lowered costs by 2.3%. The HHO algorithm gave the lowest cost, which was lower than the PSO by 2.9%. As shown in Figure 11, the HHO algorithm has an outstanding performance in terms of convergence speed as it converged after only 12 iterations.  Figure 11. The convergence curve of the proposed algorithms in Case 1.

Case 2: Reducing Power Losses
The proposed algorithms are required to lower the power losses, which means that Equation (19) is selected as the objective function for the OPF algorithm. As shown in Table 6, all proposed algorithms have succeeded in reducing the power losses and gave better results than the PSO algorithm, which is usually used in the literature.
Both MVO and GOA have set the power system parameters to give lower losses than the PSO by 5.9% and 44.1%, respectively. Again, the HHO succeeded in giving the best results, which are 56% lower than the PSO. As indicated in Figure 12, the GOA was the fastest algorithm to reach the convergence, followed by the PSO, the HHO, and the MVO. As a general evaluation for the performance of the proposed algorithms in terms of accuracy and speed of convergence, the HHO is considered the best as it gave the lowest power losses after only 24 iterations. Both MVO and GOA have set the power system parameters to give lower losses than the PSO by 5.9% and 44.1%, respectively. Again, the HHO succeeded in giving the bes results, which are 56% lower than the PSO. As indicated in Figure 12, the GOA was the fastest algorithm to reach the convergence, followed by the PSO, the HHO, and the MVO As a general evaluation for the performance of the proposed algorithms in terms of accu racy and speed of convergence, the HHO is considered the best as it gave the lowest powe losses after only 24 iterations.

Case 3: Reducing Voltage Drops
In this case, the proposed algorithms aim to reduce the voltage drops at each bus. This means that Equation (20) is selected as the objective function for the OPF algorithm. Table 7 shows the optimal results of the three proposed algorithms. All the obtained results are compared with those received from the PSO algorithms. The results shown in Table 7 clearly show the effectiveness of all proposed techniques compared with the PSO. The MVO algorithm has succeeded in obtaining costs lower than the PSO by 30.2%, while the GOA has lowered the cost by 51.3%. The HHO algorithm gave the lowest cost, which was lower than the PSO by 51%. As shown in Figure 13, the HHO algorithm has an outstanding performance in terms of convergence speed as it converged after only one iteration. It is essential to mention that the HHO algorithm has been tested many times, and at each attempt, it gave the same response.

Case 4: Reducing the Operation Costs, the Power Losses, and the Transmission Voltage Drops
In previous cases, the proposed algorithms have worked to optimize a single objective: either the operating costs or the power losses, or even the voltage drops. However, this has affected the other parameters. For example, in Table 5, the GOA algorithm succeeded in lowering the cost. However, based on these parameters, the voltage will drop

Case 4: Reducing the Operation Costs, the Power Losses, and the Transmission Voltage Drops
In previous cases, the proposed algorithms have worked to optimize a single objective: either the operating costs or the power losses, or even the voltage drops. However, this has affected the other parameters. For example, in Table 5, the GOA algorithm succeeded in lowering the cost. However, based on these parameters, the voltage will drop to reach 0.71 p.u.-which is not acceptable. To solve this problem, the proposed algorithms in Case 4 use Equation (21) as a multi-objective function which is utilized to update the parameters to lower the operation costs, the power losses, and the voltage drops. As indicated in Table 8, the proposed algorithms have given better results compared to the PSO. This case clearly shows the limitations of each method of the proposed methods applied to a multi-objective optimization function. Although the performance of the MVO was satisfactory in previous cases (as shown in Figure 14), it failed to reach the global optimum solution. Instead, it reached the optimum local value. Additionally, in terms of speed of conversion, it took too many iterations to reach convergence. Regarding the GOA, it also does not reach the best solutions for the three objectives. However, the speed of convergence is much better than the MVO. The HHO is the only algorithm that best performs as it finds global optimum solutions for the three objectives after only five iterations. optimum solution. Instead, it reached the optimum local value. Additionally, in term speed of conversion, it took too many iterations to reach convergence. Regarding GOA, it also does not reach the best solutions for the three objectives. However, the s of convergence is much better than the MVO. The HHO is the only algorithm that performs as it finds global optimum solutions for the three objectives after only five ations.

Conclusions
In this paper, three naturally inspired meta-heuristic algorithms were propose solve the OPF problem. Detailed mathematical modeling for the OPF problem was sented, including the power flow's multi-objective optimization by reducing opera fuel costs, transmission power losses, and voltage drops. The proposed techniques simulated using MATLAB and applied to the IEEE 30 bus bench-mark system to show effectiveness of each algorithm. Four case studies were formulated to assess the pe mance of each algorithm. In each case, the results of the proposed algorithms were c pared with the PSO algorithm, which is commonly used in the literature to solve the problem. The HHO algorithm showed the best performance in achieving a minimum in the first scenario, where it saved USD 24/h compared to the PSO. However, in the ond and third scenarios, the proposed HHO algorithm successfully resulted in 6 MW power loss and 0.37 p.u. more minor voltage deviation when compared to PSO res Finally, even in the multi-objective scenario, the proposed HHO proved to be a rel algorithm compared to all other algorithms under investigation. Research findings s

Conclusions
In this paper, three naturally inspired meta-heuristic algorithms were proposed to solve the OPF problem. Detailed mathematical modeling for the OPF problem was presented, including the power flow's multi-objective optimization by reducing operation fuel costs, transmission power losses, and voltage drops. The proposed techniques were simulated using MATLAB and applied to the IEEE 30 bus bench-mark system to show the effectiveness of each algorithm. Four case studies were formulated to assess the performance of each algorithm. In each case, the results of the proposed algorithms were compared with the PSO algorithm, which is commonly used in the literature to solve the OPF problem. The HHO algorithm showed the best performance in achieving a minimum cost in the first scenario, where it saved USD 24/h compared to the PSO. However, in the second and third scenarios, the proposed HHO algorithm successfully resulted in 6 MW less power loss and 0.37 p.u. more minor voltage deviation when compared to PSO results. Finally, even in the multi-objective scenario, the proposed HHO proved to be a reliable algorithm compared to all other algorithms under investigation. Research findings show that the HHO algorithm may form a very competitive algorithm for power system optimization problems.

Data Availability Statement:
No new data were created or analyzed in this study. Data sharing is not applicable to this article.