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

Abstract

Manta Ray Foraging Optimization Algorithm (MRFO) is a new bio-inspired, meta-heuristic algorithm. MRFO algorithm has been used for the first time to optimize a multi-objective problem. The best size and location of distributed generations (DG) units have been determined to optimize three different 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 different 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 different optimization algorithms in different cases. The results provide clear evidence of the superiority of MRFO that defind before (Manta Ray Foraging Optimization Algorithm. Quasi-Oppositional Differential 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 different power system problems.

1. 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 future, the shape of electrical power systems is expected to be one of four types. These types are centralized integrated with distributed generation, centralized with increased decentralization, partially decentralized, and fully decentralized. This shape will be formed based on many factors that may be summarized into three main categories. Firstly, high level of uncertainty that includes social factors, demand response, regulatory, policy and political factors; secondly, medium level of uncertainty that includes system challenges and technology requirements, and; finally, low level of uncertainty, which includes geographical, climatic considerations, availability of natural resources, population density and distribution, and existing infrastructure [2].

Attention towards DGs is increasing rapidly due to their effectiveness in solving overloading capacity problems and reducing the capital cost of the system. Furthermore, they require a small area, they are easy to construct within a small time and reduce greenhouse gases. Recently, DGs have gained much more interest due to the developed technologies of DG, restrictions applied for installing new transmission lines, increasing demand, liberalization of the electricity market, and concerns about greenhouse emissions [3]. Conventional energy resources such as thermal, hydro-power plants, have many limitations such as high investment costs for construction and transmission over a long distance; they are complex and inefficient. Distributed generations (DGs), which are decentralized power plants, are used to eliminate these limitations [4]. DGs can be categorized into two classifications: dispatchable DG, which includes all controllable generation types, such as gas turbine, combined heat and power, hydro, and biomass; and non-dispatchable DG, which includes generators characterized by uncertainty and whose generated power cannot be accurately predicted, such as wind turbines and photovoltaic cells [5].

Sizing and allocation of DG is an important issue to improve system efficiency, reliability, and power quality [6]. Allocating DG inappropriately increases system losses, and, consequently, increases total cost [7]. On the other hand, a properly sized and located DG provides many technological and economic benefits. For technological advantages, it could minimize line losses, improve power, minimize emissions of CO₂, increase security, provide better power standard, and prevent transmission and distribution congestion. For economic advantages, it saves the investment required to upgrade the services, minimizes maintenance and operation costs, reduces expenditure required for curing environmental damages, reduces fuel rate, and security costs for important loads. The use of an unsuitably placed or improperly sized DG has an adverse effect on system performance. Adding two or more DGs provides more advantages both economically and environmentally [8]. Renewable energy (RE) based DGs have been optimized for many objective functions such as: minimization of power losses, maximization of voltage stability, keeping stability margin, minimizing voltage deviation, minimizing total harmonic distortion, minimization of total costs, maximization of system reliability (minimizing system failures), maximizing RE penetration, and minimizing CO₂ emission. In order to assure system limits, some constraints are used. For instance, bus voltage limits, power flow equality constraint, overloading constraint, and DG capacity constraint [9].

A renewable energy (RE) based DG has been studied extensively during the last two decades. The effect of the Wind Turbine (WT)generator’s uncertainty on DG has been studied using different algorithms [10]. The DG-based PV cell has been overviewed taking into account optimization techniques, constraints, and optimization techniques [11]. Hybrid combinations of WT and PV cells have been studied in order to minimize power loss and total cost and to reduce CO₂ emissions [12,13]. WT and PV generators have also been combined with other types of renewable and nonrenewable DGs such as: micro turbine, fuel cell, biomass, battery storage system, and capacitor bank in order to enhance power quality and to increase benefits of these hybrid combinations [14,15,16,17]. In real life, the multi-objective optimization problem has been applied to residential homes to maximize profits and minimize electricity charges by utilizing an energy storage system (ESS), heat pump, and an electric vehicle [18]. On the microgrid level, also, hybridization between PV, WT, ESS, and a fuel cell has been applied on a remote island in Japan in order to overcome the intermittent nature of a renewable energy system, reduce energy costs, and to improve system performance [19]. Due to load variance and the need to meet load demand, an ESS has been introduced as an effective solution to the optimal unit commitment problem to minimize operating costs, maximize profits, and to assure system reliability [20].

Optimization problems can be solved using single or multi-objective functions. Multi-objective problems may be solved using one of two methods. The first is to convert the multi-objective problem into a single-objective problem using weight factors which gives a single dependent solution. Any change in the parameter’s weight, even if it is very small, causes a significant change in the prediction of the DG’s size and location. The second method is to develop meta-heuristic algorithms in order to generate a set of Pareto-optimal solutions which are non-dominated solutions [21,22,23].

Many optimization techniques have been used to optimize the DG site and size. These techniques are categorized as follows: analytical approach, classical (non-heuristic) approach, meta-heuristic optimization approach, hybrid approach, and some other assorted approaches. The performance of the analytical and classical (non-heuristic) approach is effective for simple systems, but not for complex systems. Meta-heuristic algorithms offer fast convergence with high accuracy and suitability for complex and large systems. Hybrid optimization techniques provide better performance with higher accuracy and are more suited to complex multi-objective problems. Distribution Network (DN) contains uncertainties resulting from the variable output of DG, load variation, and energy cost. This produces a complex system that is not easy to optimize. Meta-heuristic and hybrid techniques provide a better solution for optimizing the allocation and capacity of DGs [6]. Meta-heuristic algorithms are a powerful tool in optimizing a wide range of problems. They are problem-independent algorithms. Meta-heuristic algorithms are categorized into two divisions: evolutionary algorithm (EA), and swarm intelligence algorithms (SI). The evolutionary algorithm contains many types of algorithms such as: genetic algorithm (GA), evolutionary strategy (ES), evolution programming immune algorithm (AIA), and Jaya algorithm (JA). Swarm intelligence algorithms represent the biggest category that contains particle swarm optimization (PSO), shuffled frog leaping optimization (SFLO), ant lion optimization algorithm (ALOA), ant colony optimization (ACO), firefly algorithm (FFA), grey wolf optimizer (GWO), dragon algorithm (DA), whale optimization algorithm (WOA), and social learning PSO (SLPSO) [24]. Meta-heuristic algorithms can also be classified into four main sections: methods based on an evolutionary algorithm, methods based on physics, methods based on swarm algorithms, and methods based on human activities [22].

Based on the type of objective function (single or multi-objective) and optimization techniques, DG integrated with DN has been extensively studied. Minimizing power losses as a single objective function has been studied widely through different optimization techniques and different constraints. Teaching-learning-based (TLBO), WOA, DA, and Moth-Flame Optimization Algorithm (MFO) optimization algorithms have been applied to allocating and sizing DG to minimize power loss [25,26]. Determination of the size and location of the DG and capacitor bank has been studied to minimize total power loss using a hybrid approach based on GA and Artificial Bee Colony Algorithm (ABC) [27]. The coyote algorithm (COA) has effectively minimized power losses through determining the size and location of DGs, the results have also been compared to that of many types of algorithms (Genetic Algorithm (GA), Particle Swarm Optimization (PSO), Cuckoo Search (CSA), Fireworks Algorithm (FWA), Stochastic Fractal Search (SFS), Harmony Search Algorithm (HSA), and Salp Swarm Algorithm (SSA)) [28]. The size and location of DG integrated with a 15 and 33 bus system have been optimized using the whale optimization algorithm (WOA) to enhance power loss reduction [29]. The Loss Sensitivity Factor (LSF)) and PSO have been used to minimize power loss using the exact loss’ formula as a result of optimizing the size and location of the DG [30]. Compared to three types of optimization techniques (PSO, Improved Analytical (IA) method, hybrid approach), the genetic algorithm (GA) minimizes total power losses in 33 and 69 bus systems through optimizing the site and size of DG integrated with the two systems [31]. Allocation and sizing of DG connected to 33, 69, and 118-bus systems have been investigated to minimize power losses using chaotic stochastic fractal search (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), 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.

2. 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.

2.1. 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].

2.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].

2.3. 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].

3. 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].

3.1. 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:

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].

$$\overline{I_{L,k}}={\left(\frac{\overline{S_{L,k}}}{\overline{V_k}}\right)}^{*}$$

(1)

$$\overline{I_L}=\overline{I_{L,k+1}}+{\displaystyle \sum }_{p\in M}\overline{I_{L+1}}$$

(2)

$$\overline{V_{k+1}}=\overline{V_k}-\left(R_l+j{X}_l\right)\overline{I_L}$$

(3)

$\overline{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.

3.2. 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].

$$P_{l+1}=P_l-P_{L,k+1}-R_l\left(\frac{P_l^2+Q_l^2}{{\left|V_k\right|}^2}\right)$$

(4)

$$Q_{l+1}=Q_l-Q_{L,k+1}-X_l\left(\frac{P_l^2+j{Q}_l^2}{{\left|V_k\right|}^2}\right)$$

(5)

$$V_{k+1}^2=V_k^2-2\left(R_l{P}_l+X_l{Q}_l\right)+\left(R_l^2+X_l^2\right)\left(\frac{P_l^2+Q_l^2}{{\left|V_k\right|}^2}\right)$$

(6)

$$P_{loss\left(k,k+1\right)}=R_l\left(\frac{P_l^2+Q_l^2}{{\left|V_k\right|}^2}\right)$$

(7)

$$Q_{loss\left(k,k+1\right)}=X_l\left(\frac{P_l^2+Q_l^2}{{\left|V_k\right|}^2}\right)$$

(8)

$$P_{Tloss}={\displaystyle \sum }_{k=1}^{n-1}{P}_{loss\left(k,k+1\right)}$$

(9)

$$P_k=P_{k+1}+P_{L,k+1}+r_k\left(\frac{P_k^2+j{Q}_k^2}{{\left|V_k\right|}^2}\right)-P_{DG}$$

(10)

$$Q_k=Q_{k+1}+Q_{L,k+1}+x_k\left(\frac{P_k^2+j{Q}_k^2}{{\left|V_k\right|}^2}\right)-Q_{DG}$$

(11)

3.3. 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.

3.3.1. 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].

$$O{F}_1:\mathrm{Min}\left(\frac{PL}{P{L}_b}\right)$$

(12)

3.3.2. 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].

$$O{F}_2:\mathrm{Min}\left(\frac{VD}{V{D}_b}\right)$$

(13)

$$\mathrm{VD}={\displaystyle \sum }_{i=1}^n{\left(V_i-V_r\right)}^2$$

(14)

3.3.3. 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 $VS{I}_b$ is the voltage stability index of the base case [51].

$$VS{I}_j={\left|V_i\right|}^4-4\times {\left(P_j{X}_k-Q_j{R}_k\right)}^2-4\times \left(P_j{R}_k+Q_j{X}_k\right)\times {\left|V_i\right|}^2$$

(15)

$$\mathrm{Max}\left(VS{I}_j\right)=\mathrm{Min}\left(\frac{1}{VS{I}_j}\right)$$

(16)

$$O{F}_3:\mathrm{Min}\left(\frac{VS{I}_b}{VS{I}_j}\right)$$

(17)

3.3.4. The Overall Objective Function

In order to use the MRFO algorithm to solve multi-objective problems, the single-objective functions are combined together using the weight sum method as shown in Equation (18) [51]. The weighting factors are taken as ω₁ = 1, ω₂ = 0.65, ω₃ = 0.35 [51,52,53].

$$F=\mathrm{Min}\left({\omega }_1\times O{F}_1+{\omega }_2\times O{F}_2+{\omega }_3\times O{F}_3\right)$$

(18)

3.4. 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:

3.4.1. Equality Constraints

The active and reactive power flow Equations (19) and (20) are used as equality constraints [51].

$$P_0+{\displaystyle \sum }_{i=1}^{NDG}{P}_{DG}\left(i\right)={\displaystyle \sum }_{i=1}^k{P}_{Li}\left(y\right)+{\displaystyle \sum }_{j=1}^{nb}{P}_{loss}\left(j\right)$$

(19)

$$Q_0+{\displaystyle \sum }_{i=1}^{NDG}{Q}_{DG}\left(i\right)={\displaystyle \sum }_{i=1}^k{Q}_{Li}\left(y\right)+{\displaystyle \sum }_{j=1}^{nb}{Q}_{loss}\left(j\right)$$

(20)

$P_0$ and $Q_0$ stand for the base active and reactive powers of the system respectively. $P_{Li}\left(y\right)$ and $Q_{Li}\left(y\right)$ are active and reactive powers supplied from the reference bus. $NDG$ is the number of distributed integrated units.

3.4.2. Inequality Constraints

• Bus voltage constraints
  The magnitudes of voltage for all buses must be restricted in the range (0.95: 1.05)per unit, according to Equation (21) [51].

  $$V_{min}\le V_i\le V_{max}$$

  (21)
• DG sizing limits
  The distributed generators’ active and reactive powers must be within limits according to Equations (22) and (23) [51].

  $$P_{DG,min}\le P_{DG,i}\le P_{DG,max}$$

  (22)

  $$P{F}_{DG,min}\le P{F}_{DG,i}\le P{F}_{DG,max}$$

  (23)

4. Mathematical Model of MRFO

The MRFO algorithm starts to position individuals randomly using the following Equation (24):

$$X_{i,j}(:)=L{b}_{i,j}+r(:).\left(U{b}_{i,j}-L{b}_{i,j}\right)\forall i\in N_{pop}\&\&j\in N_{var}$$

(24)

The position of manta rays is expressed by $X_{i,j}(:)$; also, lower and higher boundaries are denoted by $L{b}_{i,j}$ and $U{b}_{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].

4.1. 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].

$$\begin{array}{l}{x}_i^d\left(t+1\right)\\ =\{\begin{array}{cc}{x}_i^d\left(t\right)+r.\left(x_{Best}^d\left(t\right)-x_i^d\left(t\right)\right)+\alpha .\left(x_{Best}^d\left(t\right)-x_i^d\left(t\right)\right)\hfill & \hfill i=1\\ x_i^d\left(t\right)+r.\left(x^d{}_{i-1}\left(t\right)-x_i^d\left(t\right)\right)+\alpha .\left(x_{Best}^d\left(t\right)-x_i^d\left(t\right)\right)\hfill & \hfill i=2,\dots \dots ..,N\end{array}\end{array}$$

(25)

$$\alpha =2.r.\sqrt{\left|\log \left(r\right)\right|}$$

(26)

The updated position is expressed by $x^d{}_{i-1}\left(t\right)$. Iteration number and dimension are expressed by t and d respectively. $x_i^d\left(t\right)$ 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 $\alpha$; the best position is expressed by $x_{Best}^d\left(t\right)$.

4.2. 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].

$$\{\begin{array}{c}{X}_i\left(t+1\right)=X_{best}+r\times \left(X_{i-1}\left(t\right)-X_i\left(t\right)\right)+e^{b\omega }\times \cos \left(2\pi \omega \right)\times (X_{best}-X_i\left(t\right))\\ Y_i\left(t+1\right)=Y_{best}+r\times \left(Y_{i-1}\left(t\right)-Y_i\left(t\right)\right)+e^{b\omega }\times \sin \left(2\pi \omega \right)\times (Y_{best}-Y_i\left(t\right))\end{array}$$

(27)

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].

$$\begin{array}{l}{x}_i^d\left(t+1\right)\\ =\{\begin{array}{cc}{x}_{best}^d+r\times \left(x_{best}^d\left(t\right)-x_i^d\left(t\right)\right)+\beta \times \left(x_{best}^d\left(t\right)-x_i^d\left(t\right)\right)\hfill & \hfill i=1\\ x_{best}^d+r\times \left(x^d{}_{i-1}\left(t\right)-x_i^d\left(t\right)\right)+\beta \times \left(x_{best}^d\left(t\right)-x_i^d\left(t\right)\right)\hfill & \hfill i=2,\dots \dots ..,N\end{array}\end{array}$$

(28)

$$\beta =2e^{r_1\frac{Tmax-t+1}{T}}\times \sin \left(2\pi r_1\right)$$

(29)

where β stands for a weight coefficient; T_max is the maximum number of iterations; r₁ 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_{rand}^d$ [60].

$$x_{rand}^d=L{b}^d+r\times \left(U{b}^d-L{b}^d\right)$$

(30)

$$\begin{array}{l}{x}_i^d\left(t+1\right)\\ =\{\begin{array}{cc}{x}_{rand}^d+r\times \left(x_{rand}^d\left(t\right)-x_i^d\left(t\right)\right)+\beta \times \left(x_{rand}^d\left(t\right)-x_i^d\left(t\right)\right)\hfill & \hfill i=1\\ x_{rand}^d+r\times \left(x^d{}_{i-1}\left(t\right)-x_i^d\left(t\right)\right)+\beta \times \left(x_{rand}^d\left(t\right)-x_i^d\left(t\right)\right)\hfill & \hfill i=2,\dots \dots ..,N\end{array}\end{array}$$

(31)

4.3. 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].

$$x_i^d\left(t+1\right)=x_i^d\left(t\right)+S\times \left(r_2\times x_{best}^d-r_3\times x_i^d\left(t\right)\right),i=1,\dots \dots \dots \dots N$$

(32)

where, the somersault factor is expressed by S which decides the somersault range. Both r₂ and r₃ 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].

4.4. 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].

1 -

Formulating the optimization problem and determining boundary limits.

2 -

Inserting control parameters, number of iterations (T_max), number of populations (N_pop), and somersault factor (S).

3 -

Initially positioning individuals randomly and calculating the fitness of each to determine the best solution so far.

4 -

Starting the main loop for i = 1: N_pop, If the stop criteria is not satisfied.

5 -

If Rand is >0.5, then apply cyclone foraging.

• If t/ T_max is < Rand, then update location using Equations (31)
• Else update location using Equation (28)
• End if

6 -

Else (if Rand is <0.5) apply chain foraging.

• Update location using Equation (25)
• End if

7 -

Evaluating the fitness value of each individual and updating position according to the best position.

8 -

Then, update location using somersault foraging Equation (32)

5. 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⁻¹. 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, ω₂ = 0.65, ω₃ = 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.

5.1. 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:

• Three DG units operating at a unity power factor
• Three DG units operating at 0.95 power factor
• Three DG units operating at 0.866 power factor
• Three DG units operating at an optimum power factor
• Comparing results for different power factors

5.1.1. Three DG Units at Unity Power Factor

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 Figure 8, Figure 9 and Figure 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. 33 bus system, three DG units, at unity Power factor (pf)

Reference	Method	Location	Size kW	APL kW	VD	VSI	LR (%)
Proposed	MRFO	30	1302.5	77.3793	0.0063	0.9182	63.3239
		24	1136.4
		13	962.292
[51]	QODELFA	13	964.7	77.408	0.00621	0.9182	63.31
		24	1133.4
		30	1301.7
[52]	SFSA	13	964.7	77.410	0.006232	0.9182	63.31
		24	1133.7
		30	1301.8
[53]	GA	25	909.0	95.8	0.0007	0.9701	54.6
		30	1684.0
		13	1658.0
[54]	CTLBO	13	1036.4	85.9595	0.0026	0.9481	59.26
		25	1163.1
		30	1521.7
[54]	CTLBO ε constraint	13	1.1926	96.1732	0.0009	0.9638	54.4
		24	0.8706
		30	1.6296
[56]	MOHHO	13	1207	92.95	0.002	0.9654	55.94
		25	763
		31	1400
[56]	MOIHHO	14	1223	92.25	0.0019	0.9580	56.27
		24	1144
		31	1290
[57]	MOPSO	12	1200	83.99	0.0053	0.919	60.19
		25	949.8
		33	1142.7
[57]	MOWOA	14	1021.6	79.72	0.0045	0.9249	62.214
		24	1200
		31	1200

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.

5.1.2. Three DG Units at 0.95 Power Factor

Three basic characteristics, namely power loss, bus voltage, and VSI, are simulated as shown in Figure 11, Figure 12 and Figure 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. 33-bus, three DG units at 0.95 lag pf.

Reference	Method	Location	Size		APL (kW)	VD	VSI	LR (%)
			(kW)	(KVAR)
proposed	MRFO	30	1297.20	426.36	29.1317	0.0008	0.9662	86.1922
		24	1098.00	360.88
		13	911.31	299.53
[51]	QODELFA	13	916.9	301.3	29.386	0.0007	0.9698	86.07
		24	1146.6	376.8
		30	1316.7	432.7
[52]	SFSA	13	917.4	301.5	29.383	0.000673	0.9697	86.07
		24	1146.3	376.8
		30	1315.7	432.4
[56]	MOHHO	13	1008.0	331.0	31.4	0.0005	0.976	85.1171
		25	910.0	299.0
		30	1334.0	439.0
[56]	MOIHHO	13	924	304	30.6	0.0004	0.979	85.4963
		24	1312	431
		30	1356	446

. 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.

5.1.3. Three DG Units at 0.866 Power Factor

Compared to the base case, three DG units are optimized and the simulation results of power loss, voltage magnitude, and VSI are shown in Figure 14, Figure 15 and Figure 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.

Table 3. 33-bus, three distributed generator (DG) units, 0.866 pf.

Reference	Method	Location	Size		APL (kW)	VD	VSI	LR (%)
			(kW)	(KVAR)
proposed	MRFO	13	792.710	457.72	15.4956	0.00035	0.9763	92.6554
		24	1039.7	600.35
		30	1239.6	715.75
[51]	QODELFA	13	791.1	456.7	15.498	0.0003	0.9764	92.65
		24	1041.1	599.1
		30	1243.1	717.8

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.

5.1.4. Three DG Units at Optimum Power Factor

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 three techniques, namely SFSA, MOHHO, and MOIHHO [52,56] as listed in

Table 4. 33-bus, three DG units at optimum pf.

Reference	Method	Location	DG Capacity and PF			APL (kW)	VD	VSI	LR (%)
			(kW)	(KVAR)	(Pf)
Proposed	MRFO	13	809.9112	411.27	0.8916	11.918	0.000338	0.9760	93.8559
		30	1072.4	1016.2	0.7258
		24	1079.0	533.9	0.8963
[52]	SFSA	13	834.0	391.6	0.905	11.911	0.000334	0.9763	94.35
		24	1064.8	531.4	0.895
		30	1059.2	1025.9	0.718
[56]	MOHHO	12	951	516	0.88	18.8	0.0005	0.978	91.0892
		25	786	436	0.87
		30	1381	809	0.86
[56]	MOIHHO	12	916	576	0.85	15.0	0.0003	0.978	92.8903
		24	1088	386	0.94
		30	1171	830	0.82

. 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.

5.1.5. Comparing Results Obtained for Different Power Factors

Table 5. Comparing results for different power factors.

Case	Location	DG Capacity and PF			APL (kW)	VD	VSI	LR (%)
		(kW)	(KVAR)	(Pf)
Optimum pf	13	809.9112	411.27	0.8916	11.918	0.000338	0.976	93.8559
	30	1072.4	1016.2	0.7258
	24	1079.0	533.9	0.8963
0.866 pf	13	792.710	457.72	0.866	15.4956	0.00035	0.976	92.6554
	24	1039.7	600.35	0.866
	30	1239.6	715.75	0.866
0.95 pf	30	1297.20	426.36	0.95	29.1317	0.0008	0.9662	86.1922
	24	1098.00	360.88	0.95
	13	911.31	299.53	0.95
Unity pf	30	1302.5	0	1	77.3793	0.0063	0.9182	63.3239
	24	1136.4	0	1
	13	962.292	0	1

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.

5.2. 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:

• Three DG units operating at a unity power factor
• Three DG units operating at 0.95 power factor
• Three DG units operating at 0.82 power factor
• Three DG units operating at an optimum power factor
• Comparing results for different power factors

5.2.1. Three DG Units at Unity Power Factor

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. 69-bus system, three DG units at unity Power Factor (PF)

Reference	Method	Location	Size (kW)	APL (kW)	VD	VSI	LR (%)
Proposed	MRFO	19	473.1375	71.0207	0.0022	0.9402	68.4268
		11	591.3010
		61	1859.3
[51]	QODELFA	11	629.4	72.295	0.00150	0.9525	67.87
		20	438.6
		61	1953.7
[52]	SFSA	11	570.3	72.445	0.001434	0.9537	67.80
		19	466.1
		61	1967.4
[54]	CTLPO	11	560.3	76.372	0.0008	0.9770	66.0478
		18	427.4
		61	2153.4
[54]	CTLPO ε constraint	12	965.8	79.66	0.0003	0.9770	64.5861
		25	230.7
		61	2133.6
[56]	MOHHO	20	643.6	81.0	0.0008	0.9720	63.9904
		60	971.4
		61	1328.2
[56]	MOIHHO	18	796.2	80.8	0.0007	0.978	64.0793
		61	1447.1
		64	707.5
[57]	MOPSO	21	383.6	82.79	0.0015	0.9455	63.1946
		61	1770.8
		66	1450.2
[57]	MOWOA	20	558.8	75.56	0.00083	0.9643	66.4088
		53	784.9
		61	1970.1

. 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.

5.2.2. Three DG Units at 0.95 Power Factor

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. 69-bus, three DG units at 0.95 lag pf.

Reference	Method	Location	Size		APL (kW)	VD	VSI	LR (%)
			(kW)	(KVAR)
proposed	MRFO	11	598.0106	196.5566	20.7702	0.00016	0.9772	90.77
		18	425.9067	139.9888
		61	1895.7	623.1016
[51]	QODELFA	11	579.7	190.5	20.774	0.00015	0.9770	90.77
		18	434.0	142.6
		61	1901.3	624.9
[52]	SFSA	11	5435	1786	20.727	0.000330	0.9772	90.79
		17	4132	1358
		61	1.8728	6156
[56]	MOHHO	23	519	171	30.2	0.001	0.98	86.5742
		60	1176	387
		62	1179	387
[56]	MOIHHO	13	1083	341	28.9	0.0003	0.98	87.1521
		61	799	263
		63	1229	404

. 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).

5.2.3. Three DG Units at 0.82 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. 69-bus, three DG units, 0.82 pf.

Reference	Method	Location	Size		APL (kW)	VD	VSI	LR (%)
			(kW)	(KVAR)
proposed	MRFO	11	505.5087	352.8473	4.2952	0.00010	0.9773	98.09
		61	1692.7	1181.5
		18	382.3651	266.8925
[51]	QODELFA	11	505.8	353.1	4.297	0.00010	0.9771	98.09
		18	385.9	258.9
		61	1693.9	1182.3

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.

5.2.4. Three DG Units at Optimum Power Factor

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. 69-bus, three DG units at optimum pf.

Reference	Method	Location	DG Capacity and PF			APL (kW)	VD	VSI	LR (%)
			(kW)	(KVAR)	(pf)
Proposed	MRFO	17	388.009	254.684	0.8360	4.2775	0.000106	0.9773	98.0984
		11	494.701	365.035	0.8047
		61	1680.9	1203.00	0.8132
[52]	SFSA	11	566.9	397.0	0.819	4.298	0.000116	0.9773	98.0892
		21	336.0	222.7	0.833
		61	1675.2	1178.8	0.818
[56]	MOHHO	15	332	846	0.37	21.8	0.0008	0.98	90.3085
		60	314	838	0.35
		61	1784	335	0.98
[56]	MOIHHO	13	1064	779	0.81	13.9	0.0005	0.991	93.8206
		49	1235	403	0.95
		62	1610	1181	0.81

. 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. Comparing results for different power factors.

Method	Location	DG Capacity and PF			APL (kW)	VD	VSI	LR (%)
		(kW)	(KVAR)	(Pf)
Optimum pf	17	388.009	254.684	0.8360	4.2775	0.000106	0.9773	98.0984
	11	494.701	365.035	0.8047
	61	1680.9	1203.00	0.8132
0.82 pf	11	505.5087	352.8473	0.82	4.2952	0.00010	0.9773	98.09
	61	1692.7	1181.5	0.82
	18	382.3651	266.8925	0.82
0.95 pf	11	598.0106	196.5566	0.95	20.7702	0.00016	0.9772	90.77
	18	425.9067	139.9888	0.95
	61	1895.7	623.1016	0.95
Unity pf	18	530.0261	0	1	72.7496	0.0016	0.9522	67.6582
	11	619.4599	0	1
	61	1920.7	0	1

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.

6. 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.