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.
3.1. 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
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.
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
is the
jth parameter of the
ith universe. Each parameter can be calculated from Equation (23):
where
is the
jth parameter of the
kth universe,
is a random binary number which can be either 0 or 1, and
is the inflation rate of the
ith 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
is the
jth element of the best solution (best universe created).
and
are the lower and upper bounds of this element.
,
and
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.
3.2. 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.
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
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
,
, and
are random numbers from [0,1].
The social interaction (
is obtained through the following equation:
) 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:
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 (
) 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 (
) 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:
where
N is the number of grasshoppers.
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 (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:
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.
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 () 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
,
,
, 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.
3.3. 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
is the location of the hawks in the iteration
;
is the location of the rabbit (the victim);
to
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):
where
E is the escaping energy of the rabbit,
T is the maximum number of iterations and
is the initial state of the rabbit energy. The value of
is varies between −1 and 1 based on the physical fitness of the victim. When
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 , 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.
To mathematically model the two behaviors, let
be the percentage of the successful escape of the rabbit. If
and
, 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:
Moreover,
is a random number that represents the jump strength that can get from Equation (42) as follows:
where
is a random number that varies between 0 and 1.
If
and
, 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
and
, 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
and
, 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.