Distributed Generators Optimization Based on Multi-Objective Functions Using Manta Rays Foraging Optimization Algorithm (MRFO)

: Manta Ray Foraging Optimization Algorithm (MRFO) is a new bio-inspired, meta-heuristic algorithm. MRFO algorithm has been used for the ﬁrst time to optimize a multi-objective problem. The best size and location of distributed generations (DG) units have been determined to optimize three di ﬀ erent objective functions. Minimization of active power loss, minimization of voltage deviation, and maximization of voltage stability index has been achieved through optimizing DG units under di ﬀ erent power factor values, unity, 0.95, 0.866, and optimum value. MRFO has been applied to optimize DGs integrated with two well-known radial distribution power systems: IEEE 33-bus and 69-bus systems. The simulation results have been compared to di ﬀ erent optimization algorithms in di ﬀ erent cases. The results provide clear evidence of the superiority of MRFO that deﬁnd before (Manta Ray Foraging Optimization Algorithm. Quasi-Oppositional Di ﬀ erential Evolution L é vy Flights Algorithm (QODELFA), Stochastic Fractal Search Algorithm (SFSA), Genetics Algorithm (GA), Comprehensive Teaching Learning-Based Optimization (CTLBO), Comprehensive Teaching Learning-Based Optimization (CTLBO ( ε constraint)), Multi-Objective Harris Hawks Optimization (MOHHO), Multi-Objective Improved Harris Hawks Optimization (MOIHHO), Multi-Objective Particle Swarm Optimization (MOPSO), and Multi-Objective Particle Swarm Optimization (MOWOA) in terms of power loss, Voltage Stability Index (VSI), and voltage deviation for a wide range of operating conditions. It is clear that voltage buses are improved; and power losses are decreased in both IEEE 33-bus and IEEE 69-bus system for all studied cases. MRFO algorithm gives good results with a smaller number of iterations, which means saving the time required for solving the problem and saving energy. Using the new MRFO technique has a promising future in optimizing di ﬀ erent power system problems.


Introduction
Distributed generation (DG) represents the production of electrical energy from distributed units near the consumers with a small capacity. DGs may be renewable or nonrenewable resources such as: wind turbine, solar Photo Voltaic (PV) geothermal, hydro, and diesel generators [1]. In the (CSFS) [32]. Differential evolution (DE) and moth swarm algorithms (MSA) have been successfully used to optimize distributed energy resource allocation to reduce power loss and voltage deviation [33].
Multi-objective non-dominating solution based algorithms have been used to solve optimization problems as an effective strategy. Different objective functions have been studied, many types of optimization techniques have been used, and various kinds of constraints and limitations have been applied in the optimization of DG connected to the grid. Power losses and cost functions have been simultaneously minimized by finding the optimum size and location of PV-and WT-based DG through applying the non-dominated sorting genetic algorithm II (NSGA-II) [34]. Sizing and allocation of multiple DGs with unity power factor and capacitor bank have been adjusted to optimize different objective functions. Three objective functions are optimized in parallel: total power loss, voltage deviation, and load-balancing using multi-objective (hybrid weight improved particle swarm optimization (WIPSO)-gravitational search algorithm (GSA)) or (WIPSO-GSA) [35]. The multi-objective particle swarm optimization (MOPSO) algorithm has been applied to optimize two different objective functions simultaneously: cost function, and the total power losses during failures [36]. Grid-based multi-objective harmony search algorithm (GrMHSA) has been used to determine the optimal location and capacity of DG. This is to minimize multi-objectives simultaneously, total voltage deviation, active and reactive power loss [37]. Active power loss and Fast Voltage Stability Index (FVSI)) have been simultaneously minimized using multi-objective chaotic mutation immune evolutionary programming (MOCMIEP) through determining the best size and location of distributed generation photovoltaic (DGPV) [38].
The voltage deviation at any bus is the difference between the source voltage (reference voltage) and node voltage at each bus. It would be preferred that it be minimized. The voltage deviation index is expected to be one of the most important indices for safe operation, power quality, and increasing demand; it may result gradually in deteriorating voltage stability; the voltage may collapse, and increasing node voltages may reduce reactive power losses in the system [39,40].
Voltage instability has a damaging impact on the power system which in some cases may lead to partial or total blackout. Voltage Stability Index (VSI)is an indication of the voltage level at any bus, that level must be maintained within acceptable limits. It must be always greater than zero. VSI is mainly used to determine the most sensitive buses which have the lowest VSI. The maximizing voltage stability index has been introduced as an objective function in optimizing the size and location of DG in many pieces of research. DG must be allocated at buses with the lowest VSI (most sensitive) in order to increase system stability at that point in the system [41,42].
Due to the importance of power system efficiency and fixed voltage level, minimizing power loss and improving system voltage is expected to be a major objective function. Many pieces of research added some objective functions besides power loss and voltage profile. The most known objective functions besides power loss and voltage profile are reactive power loss, voltage stability index, investment cost, operational cost, total harmonic distortion, and system load ability, which have been introduced as a weighted sum multi-objective function to determine best site and size of DG integrated with DN [43][44][45][46][47][48][49][50].
Not all optimization algorithms are suitable for solving multi-objective problems. The weight sum method is suitable to modify multi-objective functions into a single objective one. It makes matters easy for different types of algorithms to cope with multi-objective problems. This method provides a dominated solution based on pre-determined weight factors. The decision-maker should accurately determine weight factors to achieve good results. Many optimization techniques have been used to optimize size and location in order to improve voltage profile and minimize power loss. Applying three different objective functions such as minimizing power loss, maximizing voltage stability index, and minimizing voltage deviation provide good results in improving power quality and enhancing system efficiency. Many optimization techniques have been applied based on these three objective functions. These techniques are the quasi-oppositional differential evolution Lévy flights algorithm (QODELFA), a novel stochastic fractal search algorithm (SFSA), genetics algorithm (GA), Energies 2020, 13, 3847 5 of 34 a comprehensive teaching learning-based optimization (CTLBO) technique, ant-lion optimization algorithm (ALOA), improved Harris Hawks optimizer (IHHO), and multi-objective improved Harris Hawks (MOIHHO) [51][52][53][54][55][56][57][58].
Bio-inspired optimization algorithms are artificial intelligence algorithms based on the characteristics or behaviors of living creatures. They are mostly inspired by the behaviors of insects, birds, fish, or any group-living animals. Many complicated industrial and engineering problems have been solved using bio-inspired optimization techniques. They also have many features that can be extensively used. For instance, they are flexible, simple, easy to understand and implement; they have the ability to avoid trapping in local optima. Swarm intelligence (SI)-based algorithms are types of bio-inspired meta-heuristic algorithms that consider that a group of individuals (SWARM) provides a better solution compared to that of separate individuals. Compared to the evolutionary type algorithm, swarm intelligence has fewer operators, fewer parameters, memorization of best results, and ease of implementation. The manta ray foraging optimizer (MRFO) is one of the swarm intelligence (SI) based algorithms that imitate foraging behaviors -chain, cyclone, and somersault of manta rays [59].
MRFOA is a new powerful optimization algorithm. In this paper, the MRFO algorithm has been applied for the first time to optimize multi-objective problems. The multi-objective problem is converted into a single objective one using the weight sum method. Optimizing power loss, voltage deviation, and voltage stability index have been achieved through optimizing the size and location of DG units connected to the radial power system. Two radial distributed power systems are integrated with three DG units, 33-bus, and 69 bus systems. The efficiency of the proposed technique is investigated through optimizing DG units with different power factors. The rest of the paper is organized as follows: Section 2, the foraging strategies of manta ray, especially those simulated in the MRFO algorithm, are explained briefly. Then, the problem is formulated in Section 3, i.e., power flow analysis, objective function; and constraints are explained in Section 3. The mathematical model of MRFO is simulated in Section 4. The simulation results, the parameters of power systems, and the parameters of the optimization algorithm are described in Section 5. Finally, the conclusions and future work are listed in Section 6.

Manta Ray Foraging Strategies
Manta rays are one of the biggest deep-sea creatures. Their body is flat with two pectoral fins. They swim smoothly like birds in the sky. They have a huge terminal mouth behind two extended vertical lobes. Their main food is plankton which does not require sharp teeth. Manta rays and the main parts of their body are shown in Figure 1 [60]. Feeding strategies of the manta ray can be classified into eight types according to the number of rays and their swimming behavior. These strategies are: straight, surface, chain, piggy-back, somersault, cyclone, sideways, and bottom. They can also be divided into group and individual feeding. Individual feeding strategies are straight, surface, somersault, sideways, and bottom. The group foraging strategies are chain, cyclone, and piggy-back. Chain, somersault, and straight foraging are the dominant foraging behaviors at 40.7%, 29.47%, and 24.21% respectively. Piggyback foraging comes next with 2.46%, followed by cyclone foraging 2.11%. Finally, surface and sideways foraging have the lowest percentages at 0.7% and 0.35% [61]. The MRFO algorithm imitates the feeding attitude of manta rays that include chain, cyclone, and somersault foraging [60]. Only chain, cyclone, and somersault foraging strategies are going to be explained as a preamble to the MRFO algorithm. Piggyback foraging comes next with 2.46%, followed by cyclone foraging 2.11%. Finally, surface and sideways foraging have the lowest percentages at 0.7% and 0.35% [61]. The MRFO algorithm imitates the feeding attitude of manta rays that include chain, cyclone, and somersault foraging [60]. Only chain, cyclone, and somersault foraging strategies are going to be explained as a preamble to the MRFO algorithm.

Chain Feeding
A group of manta rays forms a line head-to-tail moving horizontally. They open their fins in front of their mouth and move forward and backward through the same area. Feeding runs sometimes extend to several hundred meters based on prey's concentration and distribution. At the end of feeding runs, they keep their line formed around prey and each of them updates its position slightly above or below the one in front. The group may expand to form a line of over 40 individuals at large feeding events. Six samples of manta rays' chain foraging are shown in Figure 2

Chain Feeding
A group of manta rays forms a line head-to-tail moving horizontally. They open their fins in front of their mouth and move forward and backward through the same area. Feeding runs sometimes extend to several hundred meters based on prey's concentration and distribution. At the end of feeding runs, they keep their line formed around prey and each of them updates its position slightly above or below the one in front. The group may expand to form a line of over 40 individuals at large feeding events. Six samples of manta rays' chain foraging are shown in Figure 2 [61].
Piggyback foraging comes next with 2.46%, followed by cyclone foraging 2.11%. Finally, surface and sideways foraging have the lowest percentages at 0.7% and 0.35% [61]. The MRFO algorithm imitates the feeding attitude of manta rays that include chain, cyclone, and somersault foraging [60]. Only chain, cyclone, and somersault foraging strategies are going to be explained as a preamble to the MRFO algorithm.

Chain Feeding
A group of manta rays forms a line head-to-tail moving horizontally. They open their fins in front of their mouth and move forward and backward through the same area. Feeding runs sometimes extend to several hundred meters based on prey's concentration and distribution. At the end of feeding runs, they keep their line formed around prey and each of them updates its position slightly above or below the one in front. The group may expand to form a line of over 40 individuals at large feeding events. Six samples of manta rays' chain foraging are shown in Figure 2

Cyclone Feeding
This foraging type is used when the prey is extremely concentrated in a limited area (plankton-rich water). Each individual in the mantas' line (feeding chain) circulates around itself until dragging prey into a large feeding circle. This circle's diameter is proportional to the number of animals joining the circle, approximately 15-20 m. The manta rays spiral motion lasts for a few minutes on average. The cyclone always rotates anticlockwise. The mantas' movement creates a vortex and the rotating movement of prey creates a current that draws prey outside the feeding circle towards them. The cyclone foraging technique is introduced in Figure 3 [61]. water). Each individual in the mantas' line (feeding chain) circulates around itself until dragging prey into a large feeding circle. This circle's diameter is proportional to the number of animals joining the circle, approximately 15-20 m. The manta rays spiral motion lasts for a few minutes on average. The cyclone always rotates anticlockwise. The mantas' movement creates a vortex and the rotating movement of prey creates a current that draws prey outside the feeding circle towards them. The cyclone foraging technique is introduced in Figure 3 [61].

Somersault Feeding
Somersault feeding is always backward. When the concentration of prey is high, the manta ray performs a backward feeding somersault completing a loop that has a diameter less than its body width. This type of feeding is usually performed when the prey is concentrated near the surface to restrict the prey's movement and increase feeding efficiency. During somersault feeding, they open their mouths widely and position their fins in front of their lower jaws. Figure 4 introduces different pictures that demonstrate the somersault foraging strategy [61].

Somersault Feeding
Somersault feeding is always backward. When the concentration of prey is high, the manta ray performs a backward feeding somersault completing a loop that has a diameter less than its body width. This type of feeding is usually performed when the prey is concentrated near the surface to restrict the prey's movement and increase feeding efficiency. During somersault feeding, they open their mouths widely and position their fins in front of their lower jaws. Figure 4 introduces different pictures that demonstrate the somersault foraging strategy [61].

Problem Formulation
Load flow analysis plays an important role in finding an accurate solution for system parameters. This will definitely affect finding the optimal solution for the DG installation problem [62]. The radial distribution network (RDN) has a radial structure that consists of a large number of buses and the ratio of R/X is high, which leads to a considerable voltage drop. Conventional methods such as Gauss Siedel, Newton Raphson, and fast decoupled are not efficient or fast enough. They need a large number of iterations to provide good results [63,64]. Backward/Forward Sweep BFS load flow analysis is considered to be the most efficient method used to analyze RDN. This method is based on Kirchhoff's voltage law (KVL) and Kirchhoff's current law (KCL) [65]. Applying the BFS method includes three main steps. These three steps are the backward sweep, forward sweep, and nodal current analysis. This method is based on

Problem Formulation
Load flow analysis plays an important role in finding an accurate solution for system parameters. This will definitely affect finding the optimal solution for the DG installation problem [62]. The radial distribution network (RDN) has a radial structure that consists of a large number of buses and the ratio of R/X is high, which leads to a considerable voltage drop. Conventional methods such as Gauss Siedel, Newton Raphson, and fast decoupled are not efficient or fast enough. They need a large number of iterations to provide good results [63,64]. Backward/Forward Sweep BFS load flow analysis is considered to be the most efficient method used to analyze RDN. This method is based on Kirchhoff's voltage law (KVL) and Kirchhoff's current law (KCL) [65]. Applying the BFS method includes three main steps. These three steps are the backward sweep, forward sweep, and nodal current analysis. This method is based on achieving convergence when the voltages maximum deviation is less than the tolerance error (0.000001) [66].

Backward Forward Sweep Power Flow Method
To calculate system parameters using the BFS method, three steps are used. Firstly, identify the different layers of the system as shown in Figure 5. Secondly, calculate the injected current in each phase. Then, apply backward sweep. Finally, apply forward sweep. The final step is repeated until achieving convergence and providing good results [67]. The single line diagram of radial DN is represented in Figure 6. The mathematical formulation is derived as follows: Energies 2020, 13, x FOR PEER REVIEW 9 of 37 , represents the load currents of each bus, , represents system apparent power, Vi represents bus voltage, Iz is the total branch currents, and Ip is the sub-lines current.

Power Loss Calculation
The magnitudes of active and reactive powers are calculated based on Equations (4) and (5) respectively. To estimate the value of transmission line voltages, the mathematical formulation of Equation (6) can be used. Both active and reactive power losses are calculated for any line l th using Equations (7) and (8) respectively. The total active power losses can be calculated from Equations (9). To simulate active and reactive power injection in the system at bus k + 1, Equations (10) and (11) are used , represents the load currents of each bus, , represents system apparent power, Vi represents bus voltage, Iz is the total branch currents, and Ip is the sub-lines current.

Power Loss Calculation
The magnitudes of active and reactive powers are calculated based on Equations (4) and (5) respectively. To estimate the value of transmission line voltages, the mathematical formulation of Equation (6) can be used. Both active and reactive power losses are calculated for any line l th using Equations (7) and (8) respectively. The total active power losses can be calculated from Equations (9). To simulate active and reactive power injection in the system at bus k + 1, Equations (10) and (11) are used [68]. Step 1 define system 'layers' and calculate the load current of each bus I in the phase format depending on active and reactive power and the initial bus voltage using Equation (1).
Step 2 (backward sweep), the total branch currents are calculated starting from sub-lines and moving towards the main feeder which is represented by node #1, as given in Equation (2).
Step 3 (forward sweep), start to update bus' voltage from the main feeder at node #1 and move toward the last node using the mathematical formula shown in Equation (3) [65].
Energies 2020, 13, 3847 I L,k represents the load currents of each bus, S L,k represents system apparent power, V i represents bus voltage, I z is the total branch currents, and I p is the sub-lines current.

Power Loss Calculation
The magnitudes of active and reactive powers are calculated based on Equations (4) and (5) respectively. To estimate the value of transmission line voltages, the mathematical formulation of Equation (6) can be used. Both active and reactive power losses are calculated for any line l th using Equations (7) and (8) respectively. The total active power losses can be calculated from Equations (9). To simulate active and reactive power injection in the system at bus k + 1, Equations (10) and (11) are used [68].

Objective Function
In order to optimize the size and location of DG, we have three different objective functions: minimization of power loss, minimization of voltage deviation, and maximization of voltage stability index.

Minimization of Total Active Power Loss
The total active power loss could be optimized using the following formula in Equation (12). PL is the total active power loss described by Equation (9), PL b is the base power system losses [51].

Minimization of Voltage Deviation
Based on the following mathematical formula in Equation (13), voltage deviation is minimized. The voltage deviation denoted by VD is calculated using Equation (14). Base voltage deviation is denoted by VD b , whereas V i is the bus number, and V r is the rated voltage (1.0 p.u.) [51].

Maximization of the Voltage Stability Index
The voltage stability index is calculated using Equation (15). Then, the maximization of the voltage stability index (VSI) is turned into minimization using Equation (16). Finally, the objective function could be formulated using Equation (17) whereas VSI b is the voltage stability index of the base case [51].

Operational Constraints
In order to simulate the system accurately, some constraints must be applied to the system. Three constraints must be achieved while optimizing our objective functions in this work described as follows:

Equality Constraints
The active and reactive power flow Equations (19) and (20) are used as equality constraints [51].
P 0 and Q 0 stand for the base active and reactive powers of the system respectively. P Li (y) and Q Li (y) are active and reactive powers supplied from the reference bus. NDG is the number of distributed integrated units.  [51].
• DG sizing limits The distributed generators' active and reactive powers must be within limits according to Equations (22) and (23) [51].

Mathematical Model of MRFO
The MRFO algorithm starts to position individuals randomly using the following Equation (24): The position of manta rays is expressed by X i,j (:); also, lower and higher boundaries are denoted by Lb i,j and Ub i,j respectively; while N pop and N var represent the number of populations and the number of variables respectively [69]. Manta rays have no sharp teeth; their main foodstuff is plankton which is a microscopic animal living in the water [60]. Manta rays have many foraging techniques. According to their swimming position, these techniques can be categorized into eight types. These types are: straight, surface, chain, piggy-back, somersault, cyclone, sideways, and bottom [61]. Only three strategies are simulated using the MRFO algorithm, which are chain, cyclone, and somersaulting foraging [60,69,70]. Chain and cyclone foraging are group-based techniques, but somersault foraging is an individual-based technique [61].

Chain Foraging
Plankton is not concentrated in one area. Manta rays look for their prey, i.e., plankton, and after locating its location, they swim directly towards it. Of course, the best position is that which comprises the highest plankton concentrations. Manta rays line up head-to-tail forming a chain. According to the prey's best position so far, manta rays change their position. Manta rays update their position following the mathematical formula of the next Equation (25) [60].
The updated position is expressed by x d i−1 (t). Iteration number and dimension are expressed by t and d respectively. x d i (t) is the current position of i th individual; r is a random vector extended in the range of [0-1]. The weight coefficient is denoted by α; the best position is expressed by x d Best (t).

Cyclone Foraging
The movement of each manta ray is restricted by two factors: the position of food, and the position of the one in front. Manta rays move towards plankton taking a spiral motion. This spiral-shaped movement can be simulated following a mathematical formula represented in Equation (27) [60].
where ω represents a random number extended along the range of [0-1]. This motion is extended to an n-D space formula which simulates cyclone foraging. After this modification, the cyclone foraging can be expressed using the mathematical formula in Equation (28) [60].
where β stands for a weight coefficient; T max is the maximum number of iterations; r 1 is the random number in the interval [0-1] The cyclone foraging improves both exploration and exploitation because each individual updates its position relying on the food position which varies randomly. A new random position is introduced as a new reference. As a result, the global best position is improved. This part of the optimization technique can be expressed using Equation (31), where the random position is represented by x d rand [60].

Somersault Foraging
Somersault is an individual feeding strategy. Each individual moves toward the planktons' position then somersaults to a new position. Manta rays always swim around a higher food concentration area and update its position accordingly. Somersault foraging strategy can be derived from mathematical Equation (32) [60].
where, the somersault factor is expressed by S which decides the somersault range. Both r 2 and r 3 are random numbers extended in the range of [0-1]. Each individual is free to swim between its current position and the global best position determined so far. After some iteration, all individuals come closer to the optimum solution in the search space. As a result, increasing the number of iterations will decrease the range of somersault (inversely proportional) [60].

General Procedures of the MRFO Approach
The MRFO flow chart is demonstrated in Figure 7 and can be explained through the following steps [60,69]. • If t/ T max is < Rand, then update location using Equations (

Results and Discussion
The MRFO algorithm has been applied to optimize three DG units integrated with two test systems: the IEEE-33 bus, and 69-bus systems. The power flow analysis is achieved using the backward-forward sweep method (BFS). The two systems are optimized based on three objective functions, minimizing APL, VD, and VSI −1 . DG units could be categorized based on the capability of injecting active and reactive

Results and Discussion
The MRFO algorithm has been applied to optimize three DG units integrated with two test systems: the IEEE-33 bus, and 69-bus systems. The power flow analysis is achieved using the backward-forward sweep method (BFS). The two systems are optimized based on three objective functions, minimizing APL, VD, and VSI −1 . DG units could be categorized based on the capability of injecting active and reactive powers, in other words, (operating power factor). Some types could provide only an active power such as photovoltaic and fuel cells; that means that they are operating at unity power factor. DG operating based on synchronous machines delivers both active and reactive power, and they are operated at non-unity power factor [51].
To investigate the performance of the proposed algorithm, four different cases have been studied for each system based on the operating power factor of the DG units. These four cases are unity, 0.95, 0.866, and optimum power factors operating conditions. For each case, the resulting power loss at each branch, the voltage magnitude at each bus m, and the voltage stability index (VSI) at each branch have been compared to that of the base case (without connecting DGs). The objective function is simulated using weight sum method with the following factors ω 1 = 1, ω 2 = 0.65, ω 3 = 0.35. The parameters of MRFO are defined as follows: the MRFO algorithm was simulated at 50 search agents and 50 maximum iteration numbers and the somersault factor (S) = 2. The system was simulated using MATLAB 2014b running on a system core i7, 2.7GHZ, 4.0 GB ram. The studied cases are divided into two cases based on the power system and each case divided into four sub-cases based on different power factor operating conditions.

IEEE 33-Bus System
IEEE 33 bus-system is a standard electrical power system. It has a total load demand of 3.715 MW and 2.300 MVAR at 12.6 KV [54]. In this case, we have five subcases: The simulation results of power loss, voltage magnitudes, and VSI for the whole system are compared to the results of the base case and represented by Figures 8-10, respectively. Looking at Figure 8, the power loss is decreased significantly around every branch. The maximum power loss recorded at bus 2 equals 52.0726 and is minimized to be 14.0135; the minimum power loss recorded at branch 2 equals 0.0132 and is minimized to be 0.0115. The voltage at bus 18 was the minimum value of 0.90378 and is maximized to be 0.98; the maximum voltage at bus 2 is maximized to be 0.99909 as shown in Figure 9. The VSI is simulated as shown in Figure 10. The minimum value recorded at branch 17 equals 0.6672 and is enhanced to be 0.9182; the maximum VSI at branch 1 equals 0.9881 and is improved to be 0.9963.           Table 1 shows the results obtained by the MRFO algorithm to allocating and sizing three units of DGs operating at the unity power factor. Additionally, a comparison is carried out between those results and the results obtained by nine different recent optimization techniques [51][52][53][54][55][56][57]. Three objective functions are considered: power loss minimization, voltage deviation minimization, and voltage stability index maximization. The weight factors are 1, 0.65, and 0.35 respectively [51,52]. The MRFO algorithm improves the overall system performance and provides better results compared to those obtained by other techniques. The optimum bus locations for DG units obtained by the MRFO algorithm are 30, 24, and 13; the same results are obtained by QODELFA, SFSA, and CTLBO (ε constraint) [51,52,54]. The active power loss is minimized to 77.3793 kW with a power loss reduction equal to 63.32%. The obtained result for power reduction is better than that obtained from all other techniques. In addition, the voltage deviation is minimized to be 0.0063 more than the lowest value obtained by the GA with 0.0056, which is a very small value [53]. The voltage stability index (VSI) is maximized to 0.9182 which is the same for QODELFA and SFSA, but it is a lower value compared to other techniques. The highest VSI obtained by the GA and CTLBO (ε constraint) in [53,54], but with the highest power loss values of 95.8 and 96.17 respectively. Three basic characteristics, namely power loss, bus voltage, and VSI, are simulated as shown in Figures 11-13 as a comparison between DGs operating at 0.95 power factor and the base case. This is to investigate the effect of installing DG units and the validity of the MRFO algorithm. Figure 11 represents the system power loss. The maximum power loss recorded at bus 2 is reduced to be 3.1636 kW; the minimum power loss at branch 32 is decreased to be 0.0113. The minimum voltage magnitude evaluated at bus 18 is maximized to be 0.9943 as shown in Figure 12. The minimum VSI recorded at bus 17 is maximized to 0.9773 as shown in Figure 13. The provided results are better compared to the previous case.

Three DG Units at 0.95 Power Factor
Three basic characteristics, namely power loss, bus voltage, and VSI, are simulated as shown in Figures 11-13 as a comparison between DGs operating at 0.95 power factor and the base case. This is to investigate the effect of installing DG units and the validity of the MRFO algorithm. Figure 11 represents the system power loss. The maximum power loss recorded at bus 2 is reduced to be 3.1636 kW; the minimum power loss at branch 32 is decreased to be 0.0113. The minimum voltage magnitude evaluated at bus 18 is maximized to be 0.9943 as shown in Figure 12. The minimum VSI recorded at bus 17 is maximized to 0.9773 as shown in Figure 13. The provided results are better compared to the previous case.   Three basic characteristics, namely power loss, bus voltage, and VSI, are simulated as shown in Figures 11-13 as a comparison between DGs operating at 0.95 power factor and the base case. This is to investigate the effect of installing DG units and the validity of the MRFO algorithm. Figure 11 represents the system power loss. The maximum power loss recorded at bus 2 is reduced to be 3.1636 kW; the minimum power loss at branch 32 is decreased to be 0.0113. The minimum voltage magnitude evaluated at bus 18 is maximized to be 0.9943 as shown in Figure 12. The minimum VSI recorded at bus 17 is maximized to 0.9773 as shown in Figure 13. The provided results are better compared to the previous case.     Optimization of three DG units operating at 0.95 lagging power factor based on the MRFO algorithm; and the data obtained by different techniques are shown in Table 2. The MRFO algorithm determines the best location for DG units, namely 30, 24, and 13, which are similar to results obtained by QODELFA, SFSA, and MOIHHO techniques [51,52,56]. The MRFO algorithm minimizes the total power loss to 29.13 with a total loss reduction of 86.19% which is the best value compared to all other techniques. The VSI provided by the proposed techniques is 0.966; it is lower than the highest value obtained by MOIHHO [56] with approximately 0.0128 which is a very small value. The voltage deviation optimized based on MRFO is minimized to 0.0008, which is higher than the best deviation obtained by MOIHHO with 0.0004. The overall performance is better than the results provided by the previous case.  Figures 14-16. The maximum power loss at branch 2 is minimized to 0.4296 kW as represented in Figure 14. The minimum voltage recorded at bus 18 is maximized to 0.9964 as shown in Figure 15. The minimum voltage stability index recorded at bus 17 is improved to 0.9856 as described in Figure 16. This case provides results better than in the previous two cases.
Energies 2020, 13, x FOR PEER REVIEW 20 of 37 improved to 0.9856 as described in Figure 16. This case provides results better than in the previous two cases.          Table 3 lists the data obtained for optimizing three DG units operating at 0.866 lagging power factor. MRFO algorithm is compared to the QODELFA algorithm [51]. MRFO determines buses 13, 24, and 30 as the best location similar to results obtained by QODELFA. MRFO could not provide a significant improvement. But MRFO only takes 50 iterations to achieve optimum results, which is one-fourth of the number of iterations taken by QODELFA. The obtained results are almost the same for the two algorithms. The overall performance is improved compared to the former cases. Unlike the fixed power factor operation in the previous cases, the three DGs units optimized are operating at an optimum power factor which is not fixed. Compared to the former cases, this one gives the best results. As shown in Figure 17, the maximum power loss at bus 2 is minimized to a very low  Figure 18. Finally, the VSI which records its minimum value at bus 17, is maximized to a higher value of 0.9839 as shown in Figure 19. Unlike the fixed power factor operation in the previous cases, the three DGs units optimized are operating at an optimum power factor which is not fixed. Compared to the former cases, this one gives the best results. As shown in Figure 17, the maximum power loss at bus 2 is minimized to a very low value of 0.3347 kW. The minimum voltage recorded at bus 18 is maximized to 0.9960 as demonstrated in Figure 18. Finally, the VSI which records its minimum value at bus 17, is maximized to a higher value of 0.9839 as shown in Figure 19.   The MRFO algorithm is used to optimize DG units' size and location based on three different objective functions, namely APL, VD, and VSI; the obtained results are compared to those of the other  The MRFO algorithm is used to optimize DG units' size and location based on three different objective functions, namely APL, VD, and VSI; the obtained results are compared to those of the other three techniques, namely SFSA, MOHHO, and MOIHHO [52,56] as listed in Table 4. Buses 13, 30, and 24 were determined to be the best locations to install DG units operating at optimum power factor. The same Figure 19. Characteristics of VSI, 33-bus system at optimum pf. The MRFO algorithm is used to optimize DG units' size and location based on three different objective functions, namely APL, VD, and VSI; the obtained results are compared to those of the other three techniques, namely SFSA, MOHHO, and MOIHHO [52,56] as listed in Table 4. Buses 13, 30, and 24 were determined to be the best locations to install DG units operating at optimum power factor. The same locations were obtained by SFSA [36]. The optimum power factors obtained by the MRFO algorithm are 0.8916, 0.7258, and 0.8963. The reduced power factor values mean that much more reactive power is injected into the system. As a result, the simulated characteristics are highly improved. Compared to the previous three cases, this one gives the best performance. The results obtained by MRFO and SFSA are almost the same and provide better results compared to the other two techniques, namely MOHHO and MOIHHO [56]. Although SFSA provides better loss reduction compared to our MRFO-almost 0.5%, MRFO takes only 50 iterations to get optimum results, which is half the number of iterations taken by SFSA.  Table 5 is considered to be a collection of the previous four cases to investigate the effect of operating DG units at various power factor values. For the first objective function, the total active power loss reaches its minimum value of 11.918 with a total reduction of 93.86% when the DG units were operated at an optimum power factor. In addition, the second objective function, the voltage deviation, was minimized to 0.000338 as the lowest value for all four cases when the DG units were set to an optimum power factor. Finally, the voltage stability index was maximized to 0.976 as the highest value obtained by the optimum and 0.866 pf. SFSA [52] achieves results better than the MRFO algorithm; the total loss reduction is about 0.49% more than the loss reduction achieved by MRFO. Decreasing the operating power factor from unity to the optimum value improves the performance of the power system significantly due to increasing the injected reactive power suitably with the system conditions and the objective functions needed to be optimized. According to the objective functions, three characteristics of the power system are described, namely power loss, voltage profile, and voltage stability index (VSI). The MRFO algorithm optimizes the system to improve system performance. Power loss is shown in Figure 20; the maximum power loss recorded at branch 56 equals 49.6782 kW and is minimized to 14.8851 kW. The voltage characteristic is shown in Figure 21; the minimum voltage at bus 65 equals 0.90919 and is maximized to 0.98785. The minimum voltage stability index recorded at bus 64 equals 0.6833 and is maximized to 0.9522 as shown in Figure 22.

IEEE 69-Bus System
IEEE 69 bus-system is an electrical standard power system. It has a total load demand of 3.8 MW and 2.69 MVAR at 12.6 KV [54]. In this case, we have five subcases: According to the objective functions, three characteristics of the power system are described, namely power loss, voltage profile, and voltage stability index (VSI). The MRFO algorithm optimizes the system to improve system performance. Power loss is shown in Figure 20; the maximum power loss recorded at branch 56 equals 49.6782 kW and is minimized to 14.8851 kW. The voltage characteristic is shown in Figure 21; the minimum voltage at bus 65 equals 0.90919 and is maximized to 0.98785. The minimum voltage stability index recorded at bus 64 equals 0.6833 and is maximized to 0.9522 as shown in Figure 22.   MRFOA and eight other algorithms are compared to each other in optimizing three DG units based on the pre-determined objective functions, APL, VD, and VSI for the IEEE 69-bus system as shown in Table 6. The MRFO algorithm determines the optimum locations to be 19, 11, and 61 respectively; the same as SFSA [52]. The total active power loss is minimized to 71.02 with loss reduction equal to 68.427%, which is better than the results of all other techniques. Although the voltage deviation is minimized to 0.0022, it is not the minimum deviation value. It is 0.0019 higher than the minimum value reported by CTLPO (ε constraint), [54] but the second method estimated the total capacity for DG units as 3329.4 KW which is higher than that provided by MRFO, 2923.7 with approximately 405.7 kW. The voltage stability index is maximized to 0.94 lower than the best values obtained by CTLPO and CTLPO (ε constraint) 0.977, but those techniques provide higher power loss,76 kW, 37 kW, and 79.66 kW respectively, and that makes sense.  MRFOA and eight other algorithms are compared to each other in optimizing three DG units based on the pre-determined objective functions, APL, VD, and VSI for the IEEE 69-bus system as shown in Table 6. The MRFO algorithm determines the optimum locations to be 19, 11, and 61 respectively; the same as SFSA [52]. The total active power loss is minimized to 71.02 with loss reduction equal to 68.427%, which is better than the results of all other techniques. Although the voltage deviation is minimized to 0.0022, it is not the minimum deviation value. It is 0.0019 higher than the minimum value reported by CTLPO (ε constraint), [54] but the second method estimated the total capacity for DG units as 3329.4 KW which is higher than that provided by MRFO, 2923.7 with approximately 405.7 kW. The voltage stability index is maximized to 0.94 lower than the best values obtained by CTLPO and CTLPO (ε constraint) 0.977, but those techniques provide higher power loss,76 kW, 37 kW, and 79.66 kW respectively, and that makes sense. MRFOA and eight other algorithms are compared to each other in optimizing three DG units based on the pre-determined objective functions, APL, VD, and VSI for the IEEE 69-bus system as shown in Table 6. The MRFO algorithm determines the optimum locations to be 19, 11, and 61 respectively; the same as SFSA [52]. The total active power loss is minimized to 71.02 with loss reduction equal to 68.427%, which is better than the results of all other techniques. Although the voltage deviation is minimized to 0.0022, it is not the minimum deviation value. It is 0.0019 higher than the minimum value reported by CTLPO (ε constraint), [54] but the second method estimated the total capacity for DG units as 3329.4 KW which is higher than that provided by MRFO, 2923.7 with approximately 405.7 kW. The voltage stability index is maximized to 0.94 lower than the best values obtained by CTLPO and CTLPO (ε constraint) 0.977, but those techniques provide higher power loss,76 kW, 37 kW, and 79.66 kW respectively, and that makes sense. MRFO algorithm is applied to optimize three DG units operating at 0.95 power factor to optimize the proposed objective functions. The power loss is optimized as shown in Figure 23, i.e., the maximum power loss at branch 56 is reduced to 3.7085 kW. The minimum voltage magnitude recorded at bus 65 is maximized to 0.9963 as shown in Figure 24. The minimum value of VSI at branch 64 is maximized to 0.9855 as demonstrated by Figure 25. The obtained results are better compared to the previous case.        The IEEE 69-bus system is integrated with three units of DG. These DG units are optimized to improve system performance using the MRFO algorithm. The obtained results are compared to that of four other techniques and the data is recorded in Table 7. The selected buses to install DG units are 11, 18, and 61; this is similar to results obtained by QODELFA [51]. The results obtained by the MRFO algorithm and QODELFA are almost the same. The total power loss is minimized to 20.7702 kW which is the minimum value compared to other techniques. The total installed capacity 2919 kW and 959.65 KVAR close to QODELFA, 2915 kW, and 958 KVAR [51]. The voltage deviation is minimized to a value slightly higher than that for QODELFA with about 0.00001. The voltage stability index is maximized to 0.9772 which is the best value of all. The MRFO algorithm could not provide a significant improvement compared to the QODELF technique. However, MRFO is better as it takes one-fourth of the number of iterations taken by QODELF to achieve optimum results.
The obtained results are better compared to the previous case (unify power factor). In this case, the injected reactive power is higher compared to previous cases. The obtained results are improved significantly. The power loss characteristics shown in Figure 26 show an improved performance. The power loss at branch 56 is minimized to 0.0091 kW. The minimum voltage recorded at bus 65 is maximized to be 0.9976 as shown in Figure 27. Finally, the VSI simulated in Figure 28 demonstrates that the minimum VSI at branch 64 is maximized to 0.9905. This case provides better performance of the power system due to the increment of reactive power injected into the system.       . Table 8 shows the results of optimizing DG units based on MRFO and QODELFA techniques. The given results show the effectiveness of the proposed algorithm and the other ones. The selected buses to install DG units are 11, 61, and 18. Power loss is minimized to be 4.2952, with a total loss reduction of 98.09%. The voltage deviation is minimized to its minimum value, namely 0.0001. Finally, the voltage stability index is enhanced to reach 0.9773 slightly higher than that for QODELFA. The MRFO algorithm provides almost the same results obtained by QODELFA with no superiority. QODELFA seems to be as good as MRFO, but it takes four times the number of iterations taken by MRFO. Increasing the number of iterations means increasing the total time of simulations, which decreases the effectiveness of the algorithm. The main power system characteristics, power loss, voltage magnitude, and voltage stability index optimized by the MRFO algorithm are compared to the base case. The obtained results are demonstrating that MRFO is a powerful algorithm. The power loss for each branch is shown in Figure 29. The maximum power loss recorded at branch 56 is maximized to be 0.0056 kW. The new maximum power loss recorded at branch 48 with value 1.6312 kW. The minimum voltage magnitude recorded at bus 65 is enhanced to 0.9975 as shown in Figure 30. The minimum VSI recorded at branch 64 is maximized to 0.99. The new minimum VSI recorded at branch 49 has a value of 0.9773 as shown in Figure 31. This case gives the best solution compared to the former cases.
good as MRFO, but it takes four times the number of iterations taken by MRFO. Increasing the number of iterations means increasing the total time of simulations, which decreases the effectiveness of the algorithm. The main power system characteristics, power loss, voltage magnitude, and voltage stability index optimized by the MRFO algorithm are compared to the base case. The obtained results are demonstrating that MRFO is a powerful algorithm. The power loss for each branch is shown in Figure 29. The maximum power loss recorded at branch 56 is maximized to be 0.0056 kW. The new maximum power loss recorded at branch 48 with value 1.6312 kW. The minimum voltage magnitude recorded at bus 65 is enhanced to 0.9975 as shown in Figure 30. The minimum VSI recorded at branch 64 is maximized to 0.99. The new minimum VSI recorded at branch 49 has a value of 0.9773 as shown in Figure 31. This case gives the best solution compared to the former cases.    The final sub-case in the IEEE 69-bus system is optimizing the size and location of three DG units operating at an optimum power factor. The results obtained by MRFO and three other techniques are listed in Table 9. The best locations determined using MRFO are 17, 11, and 61 which are completely different from the results of other techniques. The maximum power loss is minimized to 4.2775 kW with a total reduction of 98.0984% which is the best-obtained value. The voltage deviation also provided by the MRFO algorithm is the minimum value of 0.0001. Finally, the voltage stability index obtained by the MRFO algorithm is almost the same as the results obtained by SFSA. MRFO and SFSA techniques provide the best results compared to others. The MRFO algorithm takes a smaller number of iterations, namely 50, to achieve the optimum solution, unlike SFSA that takes double the number of iterations, 100. That gives superiority to the MRFO algorithm. The final sub-case in the IEEE 69-bus system is optimizing the size and location of three DG units operating at an optimum power factor. The results obtained by MRFO and three other techniques are listed in Table 9. The best locations determined using MRFO are 17, 11, and 61 which are completely different from the results of other techniques. The maximum power loss is minimized to 4.2775 kW with a total reduction of 98.0984% which is the best-obtained value. The voltage deviation also provided by the MRFO algorithm is the minimum value of 0.0001. Finally, the voltage stability index obtained by the MRFO algorithm is almost the same as the results obtained by SFSA. MRFO and SFSA techniques provide the best results compared to others. The MRFO algorithm takes a smaller number of iterations, namely 50, to achieve the optimum solution, unlike SFSA that takes double the number of iterations, 100. That gives superiority to the MRFO algorithm. The former four cases are compared to each other to investigate the effect of varying power factor on the performance of the system. Table 10 provides results obtained by the MRFO algorithm and represents the operating DG units at different power factors. The overall performance of the system is enhanced when moving from the unity power factor towards the optimum values. When all the DG units operate at an optimum power factor, the power loss was at its minimum of 4.2775 kW. In addition, the voltage deviation was at its minimum value of 0.0001. Finally, the voltage stability index was maximized to its higher value of 0.9773, the same as for 0.82 power factor operating conditions.

Conclusions and Future Work
A new meta-heuristic algorithm inspired by nature, namely the manta ray foraging optimization algorithm (MRFO) is applied to optimize the size and location of three DG units. Two standard power systems, 33-bus and 69-bus, are optimized based on three objective functions: power loss, voltage deviation, and voltage stability index. The proposed MRFO was tested successfully on the IEEE 33-bus system, and the IEEE 69-bus system to improve the voltage stability index, reduce voltage deviation, and decrease power losses. To investigate the performance of the proposed technique, four cases are applied to each power system categorized based on the operating power factor of the DG units. The MRFO algorithm has been compared to other techniques for the same operating system and the same operating conditions. It can be concluded that the MRFO algorithm is an effective and powerful technique. Compared to other techniques, such as: QODELFA, SFSA, CTLPO, CTLPO (ε constraint), MOHHO, MOIHHO, MOPSO, and MOWOA, it provides superiority in two cases: unity and 0.95 lagging power factor. In addition, the MRFO algorithm provides almost the same results for lagging power factor and optimum power factor for the 69 bus system. Like all other techniques, the MRFO algorithm could not provide superiority in all cases; SFSA provides better results in optimum power factor for the 33 bus system. Taking into consideration the number of iterations for all optimization algorithms, superiority would be assigned to MRFO. MRFO takes only 50 iterations to achieve an optimum value which is half and one-fourth the number of iterations taken by SFSA and QODELFA respectively. Operating DG units at unity power factor gives the worst results of all cases due to the lack of injected reactive power. Decreasing power factor means increasing injected reactive power which improves the system performance. But decreasing power factor must not exceed the optimum values of power factor to avoid over-voltage. When the DG units operate at an optimum power factor, they give the best results of all. Although the weight sum is an effective method to deal with multi-objective problems, it must be constrained to bigger objective functions such as cost function, in order to be easily compared with different references. The MRFO algorithm is a powerful technique that must be evaluated for much more complicated power systems especially a real type one.