Single and Multiobjective Optimal Reactive Power Dispatch Based on Hybrid Artiﬁcial Physics–Particle Swarm Optimization

: The optimal reactive power dispatch (ORPD) problem represents a noncontinuous, nonlinear, highly constrained optimization problem that has recently attracted wide research investigation. This paper presents a new hybridization technique for solving the ORPD problem based on the integration of particle swarm optimization (PSO) with artiﬁcial physics optimization (APO). This hybridized algorithm is tested and veriﬁed on the IEEE 30, IEEE 57, and IEEE 118 bus test systems to solve both single and multiobjective ORPD problems, considering three main aspects. These aspects include active power loss minimization, voltage deviation minimization


Introduction
Optimal reactive power dispatch (ORPD) is an important problem in power system operation. The economy of grid operation has two main aspects to consider: active (Watt) and reactive (Var) power control problems. The Watt problem concerns regulating and controlling the output of the generation units to reduce the overall costs of production. On the other hand, the Var is considered a more complex problem due to the nature of control variables involved in its operation, where it focuses on different voltage control aspects of the grid components (i.e., tap-changing transformers, reactive compensators, etc.) to reduce the overall grid losses, and improve voltage balance. ORPD is considered a pivotal problem in this manner, which aims to solve highly constrained, nonconvex, and nonlinear optimization problems that possess both discrete and continuous control variables to achieve important goals, such as minimizing active power losses and voltage deviations, while improving the voltage stability index of the grid. These types of operational issues emerge due to the complexity arising in grid modernization. Specifically, an optimal reactive power dispatch is essential to help maintain the voltage level in loading conditions by reducing the voltage deviation and power quality issues that emerge from the stochastic fluctuations of the power output. The latter due to the unpredictability of sources such as renewable energy and electric vehicle integration.
Recent years have witnessed growing attention on the metaheuristics and population-based techniques to solve various problems in power system operation and control. These modern approaches

Mathematical Formulation of the ORPD
The hybrid APOPSO algorithm developed in this work aims to individually and simultaneously minimize three main objectives in the ORPD problem subjected to equality and inequality constraints, namely, minimizing total active power loss, reducing voltage deviations, and improving voltage stability index (the L-index).

General Optimization Problem Formulation
The general mathematical formulation of an objective function can be expressed as subject to E (x, y) = 0 (2) where F is the objective function to be minimized, E and I are nonlinear constraints represented as vectors that resemble both the independent (control) and dependent variables of the problem.
Specifically, x is a vector containing the control variables of the ORPD which namely include the reactive power compensators G C , dynamic tap-setting of transformers T, and voltage levels at generation units Vg.
x can be expressed as y is a vector representing the dependent variables, which include slack generator P G slack , voltage levels at transmission lines V L , reactive power from generation units Q G , and the apparent power S L . y can be expressed as y = P G slack ; V L 1 . . . ., V L N Load ; Q G 1 . . . . , Q G Ng ; S L 1 . . . . , S L N Tr ]

Minimization of MW Losses Function
The fitness function formulated to reduce the overall MW losses in the system can be expressed as (6) where N TL is the number of transmission lines in the system, G M is the conductance of the transmission lines between the i th and the j th buses, and ∅ ij is the phase angle between buses i and j.

Minimization of Voltage Deviations
It is highly significant to record and notice all the voltage deviations across all the system's buses from reference points. This is done to ensure proper consideration of the voltage limits in optimal reactive power planning and operations, and not only as constraints, as bus voltages may operate at their maximum limits without invoking any violation, yet leading to improper reactive power reserves that could cause outages and faults. The fitness function of voltage deviation minimization can be written as where V LK is the load voltage at the k th load bus, V Desired LK is the desired load voltage at the k th load bus that is usually set to be 1.0 pu, Q GK is the reactive power from generators at the k th load bus, Q Lim KG is the generator reactive power limit, while N LB and N g are the number of the load buses and generation units in the system, respectively.

Minimization of Voltage Stability Index
The L-index is utilized in this work as a metric indicator of voltage stability performance. This index is newly introduced and presented in [24]. It assesses the steady-state voltage levels for any node in the test system and provides highly consistent results compared to other voltage stability indices, e.g., voltage collapse prediction index (VCPI), equivalent node voltage collapse index (ENVCI), fast voltage stability index (FASI), etc. [25]. Such indices are crucial factors to be considered in the planning and operation of an electrical network. L-index can be expressed as and such that where N PV is the number of the PV (generator type) buses in the system.

Equality Constraints
Equality constraints of the ORPD are the usual nonlinear power flow equations which provide voltage levels and angles at each node in the test system, expressed as where P Gi and P Di are the real power generation and demand at the system respectively while G ij and B ij are the real and imaginary entries of the bus admittance matrix corresponding with the i th row and j th column, for i = 1, . . . , N Total buses .
where Q Gi and Q Di are the reactive power generation and demand at the system respectively, for i = 1, N Total buses .

Inequality Constraints
It should be noted that both dependent and independent variables must operate within specific limits imposed on them to ensure proper operation. Therefore, we carefully consider both equality and inequality constraints in our modeling. The boundary of operation for the independent (control) variables applies to (i) the slack generator output limit, (ii) reactive power limit of the generation units, (iii) the load buses voltage limits, and (iv) the apparent power flow limit. These operational inequalities can be represented as P min where N g , N LB , and N TL are the numbers of generators, load buses, and transmission lines in the system. The dependent variables, such as (i) the generation unit voltage levels, (ii) reactive power compensators, (iii) the position of the tap transformers, are also operationally restricted and limited as follows where N C and N RT are the numbers of the compensator devices and regulating transformers, respectively.

Multiobjective Fitness Function
The previous fitness functions are incorporated into a multiobjective optimization function with penalty factors established to consider the dependent variables into the objective function minimization: where x v , x g , and x f are penalty factors incorporated to enforce the limits on the control variables to avoid any violation to the voltage deviation and stability index levels, assumed in this work to be 100.
The limit values are defined as

Mathematical Framework of the Metaheuristic Algorithms
A thorough discussion on the APO and PSO algorithms and their hybridization is presented in this section.

Artificial Physics Optimization (APO)
APO, as a naturally inspired metaheuristic methodology, is well presented in [20,26]. APO is based on the idea that an exerted force may result in either attractive or repulsive aggregation of physical entities (namely the particles or solutions) leading to a movement that represents the search to find local and global optima. Specifically, the process is based on three main observations: initializations, calculation of force, and motion of particles. At the initialization step, particles are Energies 2019, 12, 2333 6 of 24 sampled stochastically within a multidimensional decision space. The central presumption of APO is based on treating the particles (possible solutions) as physical entities that exhibit mass, position, and velocity, with the mathematical representation of the mass mapped as the fitness function. The mathematical representation of the mass (fitness) function is expressed as follow: Equations (24) and (25) can be mapped into the interval (0, 1) through an elementary transformation function. The mass functions can be rewritten as where the function f (x best ) is the objective function value at the position of the best received value for the individual (swarm particle), while f (x worst ) refers to the function value of the worst individual swarm reported: where S = {1; population of N agents} Once each particle's mass is identified, a velocity vector will be produced. The inevitable changes in velocity in the iteration process are controlled by the level and amount of force exerted on the particle, which is the second stage of the algorithm; calculation of the force, which is based on the mass of the particle and its distance from its neighbors. The force exerted on a particle i via another particle j can be found via: where F ij,k is the k th force quantity enforced on particle i via particle j in their dimensions; x i,k and x j,k are the k th dimension coordinates for the swarm particles i and j; r ij, k is the distance between these coordinates. Sgn(r) represents the signum function, while G(r) denotes the gravitational factor that follow the changes iteratively with r ij, k . Both of them can be expressed as: The g can be assumed as any value to provide simplicity and flexibility when experimenting. In our studies, we assumed these values based on studies presented in [25]. The total force exerted on all particles can be rewritten mathematically as: The third stage is understanding the motion principles of the particles in the decision space, where the computed force is utilized to determine the velocity of the particles that are used to find (and then update iteratively) the respective positions of the particles. Such motions are set in either two-or three-dimensional space, in which particles can be locally spotted, and can be mathematically represented as where V i, k and x i,k are the k th components of particle i's velocity and distance at iteration t. Beta is a uniformly distributed random number distributed on the interval (0, 1), while w is the user-specified inertia weight that can be iteratively updated, usually between 0.1 to 0.99. The inertia influences how two velocity values iteratively change. Larger values of w is a good indication of greater velocity changes, while small values is only used when we only want to facilitate a local search. Each particle identifies the information of its neighbors, while the physical attractiveness/repulsiveness rule serves as the search strategy in this algorithm to guide the population to search within the region of a possible solution in accordance with their fitness function. The high accuracy and ability to map the particle's mass as a fitness function influence the whole optimization process, to which the relationship is proportionally related; the more accurately the objective function is designed, the bigger mass will be produced, which leads to a higher level of attractiveness, or in other words, more optimized searching strength, as particles will be naturally attracted to higher masses. The iteration process in the APO leads to the updating of all particles' positions, and accordingly, the objective fitness function is adjusted to those new positions. Then, the fitness function identifies a new best individual and marks its position vector as the best solution. In this way, the second and third steps of the algorithm, force calculation and motion, are iteratively performed until a stopping criterion is achieved. Such criteria may be a predetermined number of executed iterations or reaching several successive iterations with no difference in the value of the best obtained particle position.

Particle Swarm Optimization (PSO)
PSO is a population-based, bio-inspired metaheuristic algorithm that was established by Kennedy and Eberhart in 1995 [19]. It is based on the concept of evolutionary computational method, where a system of study is started with an initial population of randomized solutions, updated iteratively in the process of searching for the local and global optima. The candidate solutions, known as particles, fly in the decision space with the velocity obtained in its previous best solutions, as well as its group's best results. Both the velocity and position of each particle are updated accordingly using the following mathematical formulas: where X ij (t) and V ij (t) are vector representations in the solution space for both the velocity and position of particle i, while Pbest and gbest are the best individual and global optimal obtained solutions. The performance of PSO as a validated and well-proven metaheuristic technique is widely spread in the literature in different fields of study. This is due to its powerful searching capacity and premature convergence without the need to find local optimal. Figure 1 shows the basic concept of the searching methodology and motion principle for particle i in PSO, where V(t), Xm, and X are three vectors describing the coordinates of the best solution in the decision space. the searching methodology and motion principle for particle i in PSO, where V(t), Xm, and X are three vectors describing the coordinates of the best solution in the decision space.

Hybridization of APO and PSO to Solve the ORPD Problem
The primary goal of establishing a hybridization of APO and PSO is to combine their individual strengths to form an optimized algorithm that utilizes the global search capabilities of APO with the strong local exploratory search performance of PSO, while improving its convergence performance. In other words, such hybridization aims to form a successful partnership among the local and global searching capabilities of the two algorithms to overcome any shortages each one may face if performed alone. Talbi et al. (2009) provide extensive analysis of the concept of integrating two metaheuristic techniques, which can be achieved on either lowly or highly heterogeneous integration [27,28]. This work combines the two algorithms as a low-heterogeneity routine. Successful implementation of the hybridization requires modifying the particles' velocity and position equations, as follows: , ( + 1) = , ( ) + , ( + 1) In our proposed APOPSO to solve the ORPD, we first define the dependent and control variables with their respective limits over a defined fitness function. Then, we randomly initialize the input values of the population (particles or swarms). Each particle represents a candidate solution. After the initialization step, we establish the best and worst values of the load flow and rank the obtained results. After that, the mass function given in Equation (26) will be assessed according to those results, and a force calculated using Equations (29) and (33) will be exerted on the particles, then their velocity and positions are updated based on Equations (36) and (37) and the distance between them according to (30), then ranking the newly produced results according to their fitness values. The process is repeated iteratively, and in each iteration we check whether there is a violation that occurred at any level to ensure proper operation within the limits. The iteration process stops updating the velocities and positions once an ending criterion is met. Figure 2 shows the flowchart of the proposed APOPSO algorithm. The pseudocode of the combined algorithm is as follows: Figure 1. Criteria of the searching in the particle swarm optimization (PSO).

Hybridization of APO and PSO to Solve the ORPD Problem
The primary goal of establishing a hybridization of APO and PSO is to combine their individual strengths to form an optimized algorithm that utilizes the global search capabilities of APO with the strong local exploratory search performance of PSO, while improving its convergence performance. In other words, such hybridization aims to form a successful partnership among the local and global searching capabilities of the two algorithms to overcome any shortages each one may face if performed alone. Talbi et al. (2009) provide extensive analysis of the concept of integrating two metaheuristic techniques, which can be achieved on either lowly or highly heterogeneous integration [27,28]. This work combines the two algorithms as a low-heterogeneity routine. Successful implementation of the hybridization requires modifying the particles' velocity and position equations, as follows: In our proposed APOPSO to solve the ORPD, we first define the dependent and control variables with their respective limits over a defined fitness function. Then, we randomly initialize the input values of the population (particles or swarms). Each particle represents a candidate solution. After the initialization step, we establish the best and worst values of the load flow and rank the obtained results. After that, the mass function given in Equation (26) will be assessed according to those results, and a force calculated using Equations (29) and (33) will be exerted on the particles, then their velocity and positions are updated based on Equations (36) and (37) and the distance between them according to (30), then ranking the newly produced results according to their fitness values. The process is repeated iteratively, and in each iteration we check whether there is a violation that occurred at any level to ensure proper operation within the limits. The iteration process stops updating the velocities and positions once an ending criterion is met. Figure 2 shows the flowchart of the proposed APOPSO algorithm. The pseudocode of the combined algorithm is as follows: Step 1: Read and evaluate the input data [Tr, Qc, . . . .]. All values must be normalized in per unit system.
Step 2: 2.1: Define the independent (control) variables X within their specific boundary levels; 2.2: Define the dependent variables Y within their specific boundary levels; 2.3: Define the fitness function with its associated penalty factors. Step 3: Generate an initial randomized population with N agents in the decision space. Specify the desired number of iterations to be performed. It should be noted that the initial positions of the population must be strictly within their boundary levels.
The initialized value of the k th control parameter in an i th particle (candidate solution) can be found using the following mathematical expression: Energies 2019, 12, x FOR PEER REVIEW 9 of 24 Figure 2. The proposed algorithm to solve the optimal reactive power dispatch (ORPD) problem.
Step 1: Read and evaluate the input data [Tr, Qc, ….]. All values must be normalized in per unit system.
Step 2: 2.1: Define the independent (control) variables X within their specific boundary levels; 2.2: Define the dependent variables Y within their specific boundary levels; 2.3: Define the fitness function with its associated penalty factors.
Step 3: Generate an initial randomized population with N agents in the decision space. Specify Rand is a number randomly allocated in the interval [0-1], while x d i,max and x d i,min are the boundary limits of the control variable d. The i th particles corresponding to the optimal dispatch problem can be rearranged in a vector form as follows: Step 4: Run the system's load flow to calculate the transmission line losses. Calculate the fitness values of all candidate solutions using the mass function. Select the minimally obtained result as best.
Step 5: Check if the control variables are within their boundary limits. If yes, proceed to step 6. If no, then penalize using the penalty function in Equation (21). The penalization is considered only for the multiobjective case studies.
Step 6: Evaluate the mass function's best and worst values using Equation (26). Calculate the force based on F i,k as in Equation (33).
Step 7: Update the velocity and position of the particles according to the modified PSO Equations (38) and (39).
Step 8: Evaluate the fitness function by newly obtained population information. Check if the control variables are within their boundary limits. If yes, proceed to Step 9. If no, then penalize using the quadratic penalty function in Equation (21).
Step 9: Record the best fitness values and rank the best obtained solutions.
Step 10: Repeat Steps 4-9 until a stopping criterion is achieved.
Step 11: Print best results and end.

Simulation and Results
This section provides an analysis of the obtained results on the three (IEEE 30, IEEE 57, and IEEE 118) test systems to verify the capacity of the proposed algorithm to solve the ORPD problem. Table 1 shows the parameters of the test systems used in our study. The test system's detailed line data and parameters can be found in [29][30][31]. The results were obtained utilizing a developed Matlab code that was run in an i7 Core, 3 GHz, 16 GB RAM computer with Matlab 2018-b. For each test system, the main outcomes to be recorded is the influence of the algorithm on the minimization of the MW losses, the minimization of the voltage deviation, and the voltage stability index (VSI) improvement (L-index), first individually and then concurrently. We implemented the APOPSO algorithm on a total population of 200 particles, with a maximum iteration run of 100.

IEEE 30 Bus Test System
The IEEE 30-bus test system has a total of six synchronous generators located at buses 1 (the slack bus), 2, 5, 8, 11, and 13, 41 transmission lines, four power transformers with off-nominal tap-settings at lines 6-9, 6-10, 4-12, and 28-27, and nine reactive power compensators at buses 10,12,15,17,20,21,23,24, and 29, respectively. This test system has a total active and reactive consumption of 2.834 and 1.262 per unit on a 100-MVA base. It has a total of 19 control variables that we considered in our study: six generator voltage inputs, four transformers tap-settings, and nine reactive power compensation devices. For this study to be precise, we constrained the level of the system's voltages to be strictly within limits of 0.95 and 1.10 p.u. The lower and upper limits of the transformer's tap settings are set to be between 0.9 and 1.1 p.u., whereas the reactive compensation devices should be ranged from 0 to 5 MVAR. Failure for a result to be within these limits will result in the use of the penalty factors we introduced in the multiobjective fitness function in Equation (21). The test system data was obtained from [29].
We measure the strength of our proposed algorithm by initially investigating each targeted objective (loss minimization, Vd minimization, and VSI improvement) individually. Table 2 shows the results obtained for the hybrid APOPSO when applied on the IEEE 30 bus system. The table illustrates the results when the simulation performed via APO and PSO separately at first before their hybridization. The results show great robustness over their integration throughout 100 consecutive trials. The initial values were proposed from previous literature conducted on the same test system, yet for different studies [32][33][34]. The optimal results obtained are then compared with previous findings reported in the literature in an attempt to illustrate the high consistency and power of the APOPSO algorithm. Table 3 presents the values for the three obtained objectives in comparison with other results reported by differential evolution algorithm (DE) [10], gravitational search algorithm (GSA) [15], particle swarm optimization with agent leader algorithm (ALCPSO) [8], and comprehensive learning particle swarm optimization (CLPSO) [35] in regard to the control variables of this study. The results show the superiority of the APOPSO algorithm over the reported results from these studies. For example, the power flow results in [8,15] specify values that are not considered optimal at some buses in the system. For instance, reactive power outputs from compensators at buses 15, 17, 20, 21, and 23 are either at or near the violation of their minimum limits, where they are too small to be operationally feasible, compared with the values obtained in DE, CLPSO, and APOPSO. The overall findings for each case study considered demonstrate APOPSO to be superior to the reported algorithms.
The final values at 100, 80, and 50 trials, respectively, for a population of 200 objects are as follow: Ploss = 4.3982 MW; Vd = 1.0477; L-index = 0.1267. Table 4 shows the statistical analysis of our simulation, with the best, worst, mean, and standard deviations over 100 trials. Figure 3 illustrates the Vd throughout these trials for the three test systems. The convergence performance of the PSO significantly improved with the hybridization, and Figure 4 shows the algorithm's fast convergence characteristic towards the optimal results considering the loss minimization case on the IEEE 30 bus system, while Figure 5 provides an insight into the statistical accuracy of the algorithm, presenting the Weibull distribution of the VSI values in 100 trials around the mean value.  convergence performance of the PSO significantly improved with the hybridization, and Figure 4 shows the algorithm's fast convergence characteristic towards the optimal results considering the loss minimization case on the IEEE 30 bus system, while Figure 5 provides an insight into the statistical accuracy of the algorithm, presenting the Weibull distribution of the VSI values in 100 trials around the mean value.   convergence performance of the PSO significantly improved with the hybridization, and Figure 4 shows the algorithm's fast convergence characteristic towards the optimal results considering the loss minimization case on the IEEE 30 bus system, while Figure 5 provides an insight into the statistical accuracy of the algorithm, presenting the Weibull distribution of the VSI values in 100 trials around the mean value.    The multiobjective values obtained for MW loss, voltage deviation and L-index are simultaneously shown in the rightmost column in Table 3, whereas Figure 6 depicts the produced Pareto optimal values for the multiobjective fitness function on the IEEE 30 test system. In addition to its superiority compared in the reported literature, the results obtained from APOPSO show no violations at any dependent variables in the system. It should be mentioned that we considered the nondomination criteria in the sorting and crowding of the distances when it comes to solving the multiobjective equation. Specifically, fuzzification of each fitness function incorporated in Equation (21) and applied on particle Z can be carried according to: where the maximum and minimum limits correspond to the objective function of the i th objective function, respectively. The normalization of contribution from each fitness function on particle Z can be calculated as: where R is the number of nondominated obtained results, and N is the total number of fitness (objective) functions.  The multiobjective values obtained for MW loss, voltage deviation and L-index are simultaneously shown in the rightmost column in Table 3, whereas Figure 6 depicts the produced Pareto optimal values for the multiobjective fitness function on the IEEE 30 test system. In addition to its superiority compared in the reported literature, the results obtained from APOPSO show no violations at any dependent variables in the system. It should be mentioned that we considered the nondomination criteria in the sorting and crowding of the distances when it comes to solving the multiobjective equation. Specifically, fuzzification of each fitness function incorporated in Equation (21) and applied on particle Z can be carried according to: Where the maximum and minimum limits correspond to the objective function of the i th objective function, respectively. The normalization of contribution from each fitness function on particle Z can be calculated as: where R is the number of nondominated obtained results, and N is the total number of fitness (objective) functions. Figure 6. Pareto optimal front for the IEEE 30 test system. Figure 6. Pareto optimal front for the IEEE 30 test system.

IEEE 57 Bus Test System
The IEEE 57-bus system consists of seven synchronous generators located at buses 1, 2, 3, 6, 8, 9, and 12, 80 transmission lines, three reactive power compensators located at buses 18, 25, and 53, and 15 transformers with off-nominal tap ratio. Line data, bus data, variable limits, and the initial values of the control variables are given in [30,31]. Twenty-five variables are set within the decision space to be investigated using the APOPSO algorithm, which are seven generator voltages, three reactive compensators, and 15 power transformers' tap-settings. Details about the tests system parameters are shown in Table 1. Table 5 presents the values obtained when applying the APOPSO algorithm on the IEEE 57 test system. It also shows the multiobjective results in the rightmost column when implementing the algorithm over the defined multi-objective function MOF function in Equation (21) to solve for the three cases with defined penalties for violation of the limits. It shows the capacity of APOPSO to efficiently optimize nonlinear, constrained problems in complex systems.
As the results show, APOPSO demonstrates greater search capabilities in finer results obtained compared with GSA [15], multiobjective differential evolution algorithm (MODE) [36], and chaotic krill herd optimization (CKHA) [37]. Further comparisons show the dominance of APOPSO over other reported findings that are not presented for simplicity of presentation. Table 4 shows the statistical analysis of the obtained results for the IEEE 57 test system over 100, 75, and 50 trials for MW loss minimization, Vd minimization, and voltage stability improvement respectively. While Figure 3 plots the Vd minimization values obtained over the trials, Figure 7 presents the Pareto optimal obtained for the IEEE 57 test system. The eminence search capabilities of the APOPSO is perfectly illustrated in the Weibull distribution presented in Figure 5, which shows the precision of the algorithm to determine the optimum values around the mean obtained over 50 trials. As the results show, APOPSO demonstrates greater search capabilities in finer results obtained compared with GSA [15], multiobjective differential evolution algorithm (MODE) [36], and chaotic krill herd optimization (CKHA) [37]. Further comparisons show the dominance of APOPSO over other reported findings that are not presented for simplicity of presentation. Table 4 shows the statistical analysis of the obtained results for the IEEE 57 test system over 100, 75, and 50 trials for MW loss minimization, Vd minimization, and voltage stability improvement respectively. While Figure 3 plots the Vd minimization values obtained over the trials, Figure 7 presents the Pareto optimal obtained for the IEEE 57 test system. The eminence search capabilities of the APOPSO is perfectly illustrated in the Weibull distribution presented in Figure 5, which shows the precision of the algorithm to determine the optimum values around the mean obtained over 50 trials.

IEEE 118 Bus Test System
For us to confirm the capacity performance and robustness of the proposed APOPSO algorithm for solving the ORPD problem on a large scale, we applied it on the IEEE 118 test system. This system's line data and parameters were adopted from [31]. It has 54 synchronous generator-buses, 64 load-buses, 12  The results for the multiobjective function in Equation (21) after its fuzzification are presented in the rightmost column of the table, and its Pareto optimal is shown in Figure 8. The results present evidence of the strength and robustness of APOPSO in solving real-life, highly complex scenarios, as the produced outcomes show good optimization for the highly-constrained, nonlinear, multiobjective fitness function in a large complex system like the IEEE 118. Table 4 presents the statistical analysis of the simulation process for each objective, while Table 6 shows the results obtained for each fitness function when simulated individually and simultaneously. No violations have been recorded in the dependent variables throughout the study. Figure 9 presents an illustration of the range of values obtained over the 100 trials for the MW loss reduction objective for this test

IEEE 118 Bus Test System
For us to confirm the capacity performance and robustness of the proposed APOPSO algorithm for solving the ORPD problem on a large scale, we applied it on the IEEE 118 test system. This system's line data and parameters were adopted from [31]. It has 54 synchronous generator-buses, 64 load-buses, 12 reactive power compensators, nine tap-setting power transformers, and 186 transmission lines. The bus voltages levels are constrained to remain within the strict boundary of the per unit range [0.95, 1.1]. Details of the test system's parameters are shown in Table 1. The table presents the values of the simulation results on the IEEE 118 test system.
The results for the multiobjective function in Equation (21) after its fuzzification are presented in the rightmost column of the table, and its Pareto optimal is shown in Figure 8. The results present evidence of the strength and robustness of APOPSO in solving real-life, highly complex scenarios, as the produced outcomes show good optimization for the highly-constrained, nonlinear, multiobjective fitness function in a large complex system like the IEEE 118. Table 4 presents the statistical analysis of the simulation process for each objective, while Table 6 shows the results obtained for each fitness function when simulated individually and simultaneously. No violations have been recorded in the dependent variables throughout the study. Figure 9 presents an illustration of the range of values obtained over the 100 trials for the MW loss reduction objective for this test system. The results show superior consistency in getting the results within a close range over a large number of iterations.

Conclusions
A new hybrid technique for solving the ORPD problem was established in this paper. The proposed APOPSO algorithm was tested and verified on three IEEE test systems: the 30, 57, and 118 bus systems. The results demonstrated the excellent search capacity of the PSO and APO algorithms when they are integrated. This process creates a robust hybrid algorithm to solve the ORPD optimization problem. Furthermore, the consistency of the optimized results in large complex systems such as the IEEE 57 and 118 bus systems prove that the combined algorithm overcomes some of the difficulties the PSO traditionally faces. For instance, being trapped in local optima, which usually lead to slower convergence. Comparison with previously reported ORPD results based on different metaheuristic methodologies verified the obtained results of the APOPSO algorithm, which shows superiority in most of the obtained results for MW loss minimization, voltage deviations minimization, and voltage stability improvement. The results obtained were based on the solution of both single and multiobjective fitness functions. The results of the hybrid algorithm produced no violations of any of the constraints placed on the dependent variables. Overall, the APOPSO algorithm shows great potential for various types of studies. Future studies may include the incorporation of electric vehicles to improve the voltage profiles' test system during specific operation criteria of the ORPD problem. The results of [38][39][40] could be used for such a study. Another aspect to investigate is the integration of flywheel energy storage systems (FESS) in the ORPD problem when integrated into the grid [41].

Conflicts of Interest:
The authors declare no conflict of interest. Number of the generation units in the system N LB Number of the load buses in the system N TL Number of the transmission lines in the system N PV

ORPD
The number of the PV (generator type) buses F ij,k The k th force quantity enforced on particle i via particle j in their dimensions F i,k The total force exerted on all particles V i, k The k th component of particle i 's velocity at iteration t x i, k The k th component of particle i 's distance at iteration t µ Z The normalization (fuzzification) factor from each fitness function on particle Z Pbest Individual swarm's best local solution gbest Best global obtained solution V LK The load voltage at the k th load bus V Desired

LK
The desired load voltage at the k th load bus G M The conductance of the transmission lines between the i th and the j th buses ∅ ij The phase angle between buses i and j