#### Abstract

This paper applies a proposed modified stochastic fractal search algorithm (MSFS) for dealing with all constraints of optimal reactive power dispatch (ORPD) and finding optimal solutions for three different cases including power loss optimization, voltage deviation optimization, and L-index optimization. The proposed MSFS method is newly constructed in the paper by modifying three new solution update mechanisms on standard stochastic fractal search algorithm (SSFS). The first modification is to keep only one formula and abandon one formula in the diffusion process while the second modification and the third modification are used in the first update and the second update. In two updates of SSFS, solutions with low quality are updated with high probability while other solutions with high quality do not get chances to be updated. This manner results in the fact that some promising solutions around the high quality solutions can be missed. In order to tackle this restriction, the second modification of MSFS is to newly update the worst solutions in the first update and the best solutions in the second update. In the third modification, all existing formulas of SSFS in the two updates are abandoned and the same new proposed technique is used for updating such solutions in two updates. Compared to SSFS, the three modifications can bring advantages to MSFS such as using smaller number of produced solutions per iteration, spending shorter execution time, finding better optimal solutions, and owning more stable search ability. Furthermore, the proposed method also sees its effectiveness and robustness over SSFS by testing on IEEE 30-bus system and IEEE 118-bus system with three different single objectives for each system. The proposed method can find less minimum, average, and maximum for all the cases in addition to faster search speed. Besides, the proposed method is also compared to other methods such as PSO-based method group, GA-based method group, DE-based method group, and other recent methods. Result comparisons also indicate that the proposed method can be more efficient than almost all these methods with respect to less minimum and smaller values of control parameters. As a result, evaluation of the performance of the proposed method is that it should be used for seeking solutions of ORPD problem.

#### 1. Introduction

Optimal reactive power dispatch (ORPD), a separate problem of optimal power flow (OPF) computation, is one of the material optimization problems for secure and stable operation in modern power systems. The main task of such ORPD problem is to find optimum settings of all control parameters such as reactive power output of generators, on-load tap changer of transformers, and shunt capacitor power outputs. In addition, dependent parameters as voltage of load buses, apparent power flow through transmission lines, and reactive power of generators must be within allowable specified range [1, 2]. In addition, different objective functions of ORPD problem comprising total transmission line power losses (*P*_{Loss}), total voltage deviation (TVD), and amelioration of voltage stability index (VSI) are always optimized. When reaching all considered objectives, the considered power systems are working economically and stably.

In the first stage, many traditional numerical methods called deterministic methods such as gradient search (GS) [3], Newton method (NM) [4], interior point method (IPM) [5–7], linear program (LP) [8–10], dynamic programming method (DPM)[11], quadratic programming method (QPM) [12, 13], and Lagrangian method (LM) [14] have found optimal solutions with acceptable quality, but these methods have met many disadvantages such as highly time consuming manner, high number of iterations, huge number of computation processes, incapability of handling nondifferentiable constraints and objective functions, and easily falling into local optimum solution zone. In this regard, the ones should be retired to make room for new algorithms that have better search ability. Recently, many metaheuristic methods inspired from nature phenomenon or behavior of animals have been more widely and successfully applied for solving such ORPD problem. Many methods have been continually grown and become a big family of methods like the variants of genetic algorithm (GA) [15–19], variants of differential evolution (DE) [20–24], variants of particle swarm optimization (PSO) [25–31], variants of gravitational search algorithm (GSA) [32–35], and many new standard methods [36–49]. In adaptive genetic algorithm (AGA) [15], the method changed both mutation probability and crossover probability based on comparison of the maximum fitness value and average fitness value of the population to enhance global search quality and fast convergence speed. Authors in [16] proposed an improved genetic algorithm (IGA) by using real code for continuous variables and binary code for discrete variables to decrease computation time. However, results of IGA method compared to those of real GA and Binary GA showed that its computation time was insignificantly improved. A version of GA called an efficient genetic algorithm (EGA) in [17] has applied Misum fuzzy goal programming model to find important weights of the goal function. This function was considered as the fitness function of GA for solving ORPD problem. As a result, EGA has decreased computation process and avoided approximation error of classical linearization methods. Alam and De proposed hybrid loop-genetic based algorithm (HLGA) in [18]. HLGA has used the global search feature of GA to find the good search space and then implemented the local search process utilizing loop optimization method within the identified spaces. When comparing to GA, the method has been faster convergent to optimal solutions, but there was no improvement of such obtained results. Thus, HLGA has still coped with a disadvantage so that solution quality was low. DE group has been introduced to ORPD problem such as standard differential evolution (DE) [20], modified differential evolution algorithm (MDE) [21], quasi-oppositional differential evolution (QODE) [22], multiobjective differential evolution (MODE) [23], a novel hybrid chaotic artificial bee colony, and differential evolution algorithm (CABC-DE) [24]. These improved DE methods were formed mainly based on the improved mutation, crossover, and selection operation or the combination of DE and other algorithms for surmounting several drawbacks of DE such as easily getting in stuck in local optimal zones and premature convergence. MDE has applied time-varying chaotic mutation rate* F* and crossover rate* CR* to replace standard operations of DE for time saving and reduction of setting the parameters. This method has been also proven on IEEE 30-bus system and compared to three other DE methods such as conventional DE, self-adaptive DE (SaDE), and E with integration of new mutation and crossover strategies (NMCDE). The comparison results of QODE and other methods showed that the standard deviation of active power loss, voltage deviation, and voltage stability index for all test systems was very small. So, readers cannot see the effectiveness and robustness of QODE over other ones. CABC-DE was an improved version of DE based on the combination of basic characteristic of DE and chaotic artificial bee colony algorithm (CABC). In such method, search capacity of DE was significantly improved by applying the foraging behaviors of observation bees and the search operations of detective bees from artificial bee colony algorithm (ABC). CABC-DE has been compared to the standard DE via the results obtained from IEEE 14-bus system and IEEE 30-bus systems. Simple comparison via power loss could see that this method could get good solution with small populations than classical DE. PSO and its improved versions have been applied for ORPD problem such as the improved particle swarm optimization (IPSO) [25], modified particle swarm optimization (MPSO) [26], comprehensive-learning particle swarm optimization (CLPSO) [27], improved pseudogradient search-particle swarm optimization (IPGS-PSO) [28], hybrid particle swarm optimization and imperialist competitive algorithms (HPSO-ICA) [29], fuzzy adaptive heterogeneous comprehensive-learning particle swarm optimization (FAHCLPSO) [30], and hybrid cuckoo search and particle swarm optimization (CS-PSO) [31]. In conventional PSO, the update operation of each particle was dependent on the best position of each particle (*P*_{best}) and the best position of the whole particles (*G*_{best}). With searching process, most of particles always updated around in specified position. Thus, PSO only explored zones near the best particle, which could be in local optimal zones. In order to tackle the drawback, IPSO used a new update operation that each particle would be updated by the best position of it, the best position of the groups, and the best position of the neighboring particles. By using the mechanism, capacity of finding the global solution of IPSO was enhanced better and its convergence speed was faster than PSO. MPSO applied the Michalewicz’s nonuniform mutation theory into PSO to diversify particle position search processes. Adding new strategy, MPSO decreased the possibility of getting stuck in local zones. However, the effectiveness of MPSO has been only investigated for transmission line power loss objective while other cases like total voltage deviation (*TVD*) and voltage stability index have not been demonstrated. CLPSO presented in [27] has proposed three modifications relying on the structure of PSO. Firstly, the velocity of each particle could be updated by using its* P*_{best} or potentially other particle* G*_{best} in the group. Secondly, a new formula has been introduced to adaptively obtain the inertia weight values. Lastly, a new range has been suggested to limit the velocity value. With above modifications, CLSO has overcome the premature convergence of PSO. Another improved version of PSO, called the improved pseudogradient search-particle swarm optimization (IPGS-PSO), has been applied to solve ORPD problem with three independent objectives comprising power loss, voltage deviation, and voltage stabilization index. Pseudogradient theory has been added to PSO to find the best velocity direction for particles. So, IPGS-PSO could determine better direction for velocity update. As a result, the method was a more effective method than other PSO methods such as PSO with time-varying inertia weight (PSO-TVIW), PSO with time-varying acceleration coefficients (PSO-TVAC), self-organizing PSO with TVAC (SPSO-TVAC), PSO with constriction factor PSO-CF, pseudogradient search PSO (PGS-PSO), stochastic weight trade-off PSO (SWT-PSO), and pseudogradient search and SWT-PSO (PGSWT-PSO). HPSO-ICA was a combination of the imperialist competitive algorithm and improved PSO. In HPSO-ICA, authors supplemented an independent country into imperialist competitive algorithm to enhance the searching capability of the method. By this work, this method has spent more time than PSO and ICA searching solutions. FAHCLPSO has applied the Mamdani fuzzy rules to adjust the inertia weight of Heterogeneous Comprehensive-Learning Particle Swarm Optimization (HCLPSO) aiming to avoid trapping in local optimal zones. Via the comparison of results obtained from IEEE 30-bus system, IEEE 118-bus system, and IEEE 354-bus system with the objective of power loss minimization and voltage deviation optimization, FAHCLPSO has shown its easy application and its highly computational efficiency when comparing other PSO methods and metaheuristic algorithms. In addition to the presence of the three largest groups, some small groups have been also introduced to tackle ORPD problem such as the gravitational search algorithm (GSA) [32–34], improved GSA with feasible conditional selection strategies (IGSA-FCSS) [35], quasi-oppositional teaching learning based optimization (QOTLBO) [36], teaching learning based optimization (TLBO) [36], modified Gaussian barebones based TLBO (MGBTLBO) [37], and Gaussian barebones based TLBO (GBTLBO) [37]. From 2015 to 2017, a high number of methods were employed for ORPD problem such as the hybrid Nelder–Mead simplex based firefly algorithm (HNMS-FA) [38], Artificial Bee Colony Algorithm (ABC) [39], differential search algorithm (DSA) [40], exchange market algorithm (EMA) [41], chaotic krill herd algorithm (CKHA) [42], gray wolf optimizer (GWO) [43], Gaussian barebones water cycle algorithm (GBBWCA) [44], ant lion optimizer (ALO) [45], moth-flame optimization technique (MFOT) [46], whale optimization algorithm (WOA) [47], Ant Colony Optimization Algorithm (ACOA) [48], and backtracking search algorithm (BTSA) [49]. All in all, most of these methods had a strong search ability and outperformed deterministic algorithms, original metaheuristic algorithms in terms of solution quality, computing time, and convergence speed.

In the paper, we have proposed a modified stochastic fractal search algorithm (MSFS) for handling constraints and finding three different single objectives of ORPD problem. Standard stochastic fractal search algorithm (SSFS) was first introduced in 2015 by Salimi [50]. The author applied SSFS for twenty-three benchmark functions and compared results with other popular metaheuristic methods such as cuckoo search algorithm, modified cuckoo search algorithm, backtracking search optimization algorithm, derandomized evolution strategy with covariance matrix adaptation (CMA-ES), animal migration optimization, PSO, DE, ABC, and GSA. The minimum of all benchmark functions found by SSFS was less than those from all methods, but it should be pointed out that the author has used the same population size (*N*_{pop}) and maximum number of iteration (*Max*_{Iter}) for all compared methods to SSFS. It is clear that the selection is not fair for these implemented methods because total number of produced solutions of these methods for each run is not equal. Especially, SSFS has used very high number of diffusions (*N*_{d}) equaling to 10 while the total produced solutions per iteration (*N*_{spi}) by SSFS is (). For other methods, is much smaller than that of SSFS; for example, it is approximately equal to () for CSA, () for PSO, and DE. Consequently, the real performance of SSFS could not be demonstrated persuasively. In spite of the limitation, SSFS has been also applied for different optimization problems in electrical engineering fields such as status estimation of power system [51], DC motor control [52], environment and economic dispatch [53], PID voltage regulator [54], wind integrated multiobjective power dispatch [55], and position determination of Distributed Generators in distribution power network [56]. These studies have tried to prove the outstanding performance of SSFS over other ones via comparison of output parameters such as the fuel cost, emission, speed response, and obtained voltage; however, they have forgotten that convergence speed reflected via comparison of* N*_{spi} and* Max*_{Iter} is also a very important comparison criterion in addition to such objectives. Thus, other authors have investigated the real performance of SSFS and they have identified the algorithm has needed more improvements for avoiding drawbacks such as low effectiveness for large scale problem, complicated and nondifferentiable objective functions, and high probability of finding local optimums. As a result, improved versions of SSFS have been introduced to readers with difference applications in mathematical and electrical engineering fields. These methods are SFS with DE and Gaussian walk-based diffusion (SFS-DE) [57], hybrid stochastic fractal search and pattern search technique (HSFS-PS) [58], hybrid SFS based on different mutation models of DE (HSFS-DE) [59], improved SFS (ISFS) [60], and chaos-based SFS (CSFS) [61]. In general, almost all improved methods have used different mutation models of DE available in [62]. SFS-DE has applied current-to-best mutation model for solutions with a probability of 50%. If random number was lower than 0.5, the current-to-best mutation model has been employed. Otherwise, the standard diffusion process of SSFS has still been applied. This method has been tested on thirty benchmark functions and compared to other six methods including different improved versions of DE and SSFS. The authors have set control parameters for all such implemented methods based on fair comparison with the same number of fitness evaluations. The method has found better results than all the compared methods for half of benchmark functions and it has gotten worse results for another half. Clearly, the performance of the method is still a question for readers to decide if it is really effective for their studied problems. HSFS-PS was a combination of SSFS and pattern search (PS) algorithm in which SSFS was completely and successfully carried out first and then optimal solution was used to run PS algorithm. Thus, compared to SSFS and PS, the method was more complicated and highly time consuming. HSFS-PS has been employed for designing different controllers based on PI, PID, and cascaded PI-PD. Final evaluation in the study was that cascaded PI-PD controller, which was yielded by applying HSFS-PS, was superior to remaining controllers applying HSFS-PS. So, the performance of HSFS-PS is now still a question for readers because it was not proven to be better than any optimization algorithms, even SSFS and PS. HSFS-DE in [59] has abandoned original diffusion process of SSFS, and added mutation technique based on target-to-best/1 in [62] for new diffusion process. In addition, self-adaptive crossover technique in [63] has been applied to retain new solutions from such new diffusion process. This method together with PSO, DE, and SSFS has been implemented for tuning parameter of chaotic system based on fractional order. HSFS-DE together with such three compared ones has been fairly implemented by thoroughly selecting population size as well as the number of iterations. Simulation results have indicated that HSFS-DE could reach the highest accuracy and fastest response with very low oscillation for Chen, Lorenz, and Chua systems. The proposed modification on ISFS in [60] and CSFS in [61] is totally different from those in [57–60]. In fact, ISFS has suggested using sigmoid, inverse cosine, exponential, hyperbolic cosine, and hyperbolic cosecant function for calculating standard deviation *δ* in two formulas of diffusion process while CSFS has tried to apply different chaotic maps for replacing uniform random number within zero and one, which was used in the first formula of diffusion process. ISFS has been applied to different IEEE systems for comparison with SSFS, but the exact name of the considered problem has not been mentioned in [60]. The study also pointed out that ISFS was superior to SSFS in terms of better results and much faster execution time; however, control parameters consisting of population size* N*_{pop}, number of diffusion* N*_{d}, and maximum number of iteration* Max*_{Iter} have not been reported in the study. Clearly, the outstanding performance of ISFS was not demonstrated persuasively and it can be wondered if readers should use it as the most reasonable choice. Ten versions of chaotic maps have been employed, forming ten CSFS methods in [61]. Three main control parameters have been set equally for both SSFS and ten CSFS methods, and result comparisons have shown that the performance of such ten methods has not been stable and superior to SSFS because for some study cases SSFS found worse results than them but for other study cases SSFS found better results than them. The unstable performance between SSFS and such ten CSFS is due to the impact of randomization in formulas in diffusion, the first and the second update processes. Furthermore, the study has stated that such ten methods have been more effective and robust than PSO; however, the authors have set smaller values for control parameters for PSO. For example, for a 33-bus system* N*_{pop}=100 and* Max*_{Iter}=100 were set for PSO but* N*_{pop}=50,* Max*_{Iter}=100 and* N*_{d}=5 were set for such ten CSFS methods. Thus, through the whole search CSFS has produced approximately (505+250)100=35,000 solutions while PSO has produced (100100)=10,000 solutions. Clearly, CSFS could take more advantages for finding better results than PSO. The selection manner continued to be used for 69-bus system and 118-bus system, and CSFS methods continued to show better results than PSO. In summary, authors have proposed the different improved versions of SSFS to get the best performance of their suggested methods, but they have forgotten to select fair comparison criteria with intent to give an accurate evaluation of performance between implemented. In the paper, we have proposed a modified stochastic search algorithm (MSFS) by applying three modifications on three new solution generations of SSFS. In the first generation, we have used only one formula and canceled one formula while in the second and third new solution generations, we have proposed a new technique for generating new solutions. In SSFS, there are two different formulas used in each new solution generation by a condition of comparison between a random number and a predetermined number, which are normally selected to be 0.5. The fact that the second formula in each generation is less effective than the first formula along with use of both formulas leads to time consuming manner as well as high number of computations. Thus, in the new proposed technique we have given only one new formula with greatly effective performance. In addition, we have changed old solutions, which are selected to be newly updated. Namely, SSFS updates approximately all solutions for the first and the second update processes, but we have suggested updating some worst solutions for the first update process but updating other best solutions for the second update. The sum of updated solutions in the first and the second updates is equal to* N*_{pop} while that number of SSFS is . By applying the three modified techniques, our contributions are as follows:(i)Indicate disadvantages of formulas used in SSFS.(ii)Propose to keep high efficient formulas and abandon ineffective ones.(iii)Cancel the update of* N*_{pop} solutions and calculation of* N*_{pop} evaluation functions, leading to reduction of high number of computations as well as execution time.(iv)Propose new techniques for producing new solutions in the first and the second update processes.(v)Find better results for different systems of ORPD problem.

In the paper, the proposed method together with SSFS using Gaussian random walk (SSFS-Gauu) and SSFS with Levy random walk (SSFS-Levy) have been applied for seeking three different objectives of two standard IEEE systems with 30 and 118 buses of ORPD problem. The results obtained by the three methods consisting of minimum, mean, maximum, and standard deviation in addition to control parameters are compared with each other. Furthermore, the results of the proposed method are also compared to those from other existing methods mentioning introduction for further investigation of the proposed method performance.

#### 2. Optimal Reactive Power Dispatch Problem Formulation

ORPD problem can be mathematically expressed in terms of objective functions and a set of constraints. Normally, ORPD considers three different targets consisting of minimization of total active power loss, minimization of total voltage deviation of all load buses, and minimization of voltage stability L-index. Besides, all constraints in transmission power network such as voltage bounds of all load buses, upper bound of apparent power flow through transmission branches, reactive power bounds and voltage bounds of generators, and real and reactive power flow balance at each bus. The detail of objectives and constraints is as follows.

##### 2.1. Objective Functions

(1)* Reduction of total active power loss: *total active power loss* P*_{Loss} in all branches of high voltage transmission power networks is significant, causing ineffectiveness for power system. Thus, reducing total reactive power loss is one of the most important targets when operating power system. The objective is described by [20]:(2)* Reduction of total voltage deviation of all load buses: *power quality consists of two most important factors, voltage and frequency, in which frequency can be controlled by tuning reactive power balance while voltage control is dependent on many factors such as transformer, reactive power of generator, and total reactive power loss in branches and capacitors. Voltage control is effective when total voltage deviation is very small, nearby zero. The objective can be described by the following model [28]: where* V*_{loadi} is voltage of load bus and 1 Pu is the expected value of each load bus.(3)* Reduction of L-index: *voltage instability is one of the most harmful phenomenon for power system, which can cause voltage collapse gradually even immediately. Thus, the enhancement of voltage stability is one of the most important targets as two objective mentioned above. The enhancement of voltage stability is equivalent to minimization of voltage stability indicator, normally called L-index at each bus in power system. The enhancement of voltage stability is carried out by reducing the highest value of L-index at one bus in the power system. Namely, L-index at each bus* i *and the voltage stability enhancement are mathematically formulated by [11]:where* L*_{i} is the L-index value of bus* i* and* Y*_{ij} is the mutual admittance between buses* i* and* j*.

##### 2.2. Constraints

*Active and reactive power flow balance at each bus: * at each bus, active power flow is the sum of three terms, generator, load, and lines while reactive power flow at each bus is the sum of four terms, generator, load, lines, and capacitors. The balance of all the terms must be exactly met by the following models [27]: In addition, all electric components in power systems such as generators, transformers, conductors, capacitor banks, and loads are also constrained by security work capability. Thus, electrical bounds of these components are also considered seriously during operation process. In fact, they are subject to the following rules [27]:

#### 3. The Proposed Modified Stochastic Fractal Search Algorithm (MSFS)

##### 3.1. Overview on Standard Stochastic Fractal Search Algorithm (SSFS)

SSFS comprised three new solution generations corresponding to three main stages, diffusion process, the first update process, and the second update process in which diffusion process is the first generation, the first update is the second new solution generation and the second update is the third new solution generation. In the diffusion process, each point will be diffused into* N*_{d} different points. Therefore, there are* N*_{d} times number of points, which are newly updated in the diffusion. SSFS uses Gaussian walk to produce new solutions for diffusion process in which there are two options for each solution. If random number of the* dth* solution* rand*_{d} are less than* walk*, which is selected within 0 and 1, the first option is selected. Otherwise, the second option is selected. The diffusion process is described by the following model: where point_{best} is the best point among the current population;* rand*_{d} is a random number for* point*_{d}; and *α* is obtained bywhere* CI* is current iteration.

The first update process is carried out by using only formula (15) and the second update process is based on (16). However, not every solution in the whole population is updated but depends on another condition. Before applying the first and second update processes, the whole population is evaluated and ranked in descending order of fitness function value. Then, the order of solutions is set to their rank, solution* d* corresponding to* rank*_{d}. For example, the worst solution is with* rank*_{1}= 1 and the best solution is with rank_{Npop}=*N*_{pop}. If a random number* λ*_{d} within 0 and 1 of solution* d* is higher than* rank*_{d}/*N*_{pop}, solution* d* is newly updated by using (15) or (16). Otherwise, it remains unchanged. The whole search process of SSFS can be described in flowchart of Figure 1.

##### 3.2. The Proposed Modified Stochastic Fractal Search Algorithm

In the proposed MSFS, we carry out three modifications on diffusion process, the first update and the second update. The modifications are described in detail as follows:

###### 3.2.1. The First Modification

In the first modification, we suggest using the first model for all points by setting walk to 1. It is clearly seen that the first model considers the impact of each point and the best point for producing updated step size while the second model only utilizes Gaussian random walk. Control variables in each point are dependent on different considered optimization problem, so Gaussian random walk cannot result in promising solutions with high quality. As a result, diffusion is carried out by using the following formula:

Furthermore, using only the proposed formula can reduce the selection of values for parameter* walk* and shorten simulation time that SSFS has coped with so far.

###### 3.2.2. The Second Modification

In the second modification, we change update strategy in the first and second update processes. SSFS updates almost all solutions depending on a comparison of random number for the two update processes shown in step 10 and step 16 in Figure 1. Thanks to the steps, it can be seen that almost all the worst solutions are newly updated with very high probability while the best solutions are hardly ever newly updated for the two processes. Clearly, the manner can miss promising solutions, which may be around the best solutions. Furthermore, the search ability of (15) in the first update, and (16) in the second update is less effective. Thus, the generation of* N*_{pop} solutions in the first update and the generation of* N*_{pop} solutions in the second update are time consuming and ineffective. In order to avoid overlapping the disadvantages of SSFS, the second proposed modification is applied. We choose a fraction* P*_{a} of the population, (*P*_{a}*N*_{pop}) solutions with the low quality to be updated in the first update process while [(*1* -*P*_{a})*N*_{pop}] solutions with high quality will be updated for the second update process. Besides, we also propose other formulas, which are different from formulas (15) and (16), for updating new solutions in the first and the second update processes. The change will be presented in the third modification below.

###### 3.2.3. The Third Modification

In the third modification, we propose a new technique for updating new solutions of the first and the second update processes. SSFS uses formulas (15) and (16), which are similar to mutation operation of classical DE in [20] and coping with different limitations such as easily trapping in local optimal zones and hardly ever converging to global optimum. Furthermore, both formula (15) and the first model in formula (16) are totally different from other metaheuristic since new solutions are updated by subtracting an updated step size from old solutions while approximately all methods add an updated step size to old solutions for finding new solutions, namely, PSO [25–31] with velocity and position calculation formulas, Cuckoo search algorithm [66] with Levy flight formula and strange egg discovery formula, Bat algorithm [67] with formulas calculating frequency, velocity and position, DE [20] with mutation formula, and GA [65] with polynomial mutation formula, etc. Improved social spider optimization (ISSO) was proposed in [64] by performing several modifications on social spider optimization (SSO) in which one main modification was to abandon one formula with subtracting updated step sizes and keep another with adding updated step size. ISSO was more effective than SSO in terms of capability of finding the best solutions and shortening execution time. The review on many metaheuristic algorithms indicates that SSFS can miss the promising solutions with high quality. As a result, in the proposed technique, we suggest formulas like (15) and the first model in formula (16) should be canceled and replaced with more effective formulas.

At the beginning, we introduce two new formulas (18) and (19), which are used to measure quality of each solution among the current population. In the two formulas,* ∆FF*_{d} is fitness function deviation between fitness function of the* dth* point* FF(Point*_{d}*)* and that of best point* FF(Point*_{best}*)* while* ∆FF*_{mean} is the fitness function deviation between mean fitness function of population* FF*_{mean} and fitness function of the best point* FF(Point*_{best}*)*. We use the two objects to measure the quality of each considered solution wherein solution* d* with* ∆FF*_{d}>* ∆FF*_{mean} is considered to be bad solution and solution* d* with another case is evaluated to be good solution. A strategy of finding new optimal solutions is plotted in Figure 2 by the presence of five current points (consisting of two good points in blue, two bad points in black, and the best point in red) and four newly updated points (two of bad points in yellow and two of good point in green). The figure shows that good points are newly updated around themselves while newly updated bad points can be found nearby the best point. Thus, search space around many good points is expanded while narrow space around the current best solution is also exploited strongly. The search strategy can explore global search zones and also exploit local search zones effectively and quickly. For simplicity of implementing the search technique, a pseudocode is shown in Algorithm 1 and the whole search of the proposed MSFS method is described in Figure 3.

*Algorithm 1. *New solution search strategy for the first and the second update processes of MSFS For* d*=1:*N*_{pop} If ∆*FF*_{d}>∆*FF*_{mean} else end end

#### 4. The Application of the Proposed Method for Solving ORPD Problem

##### 4.1. Selection of Control Variables and Producing a Set of Points

Similarly to optimal power flow problem, ORPD problem is also supported by using Matpower program for obtaining the remaining variables in a considered transmission power network. Matpower needs input data consisting of active power of generators, which is given data, and a set of control variables such as generator voltage* V*_{Gi}, transformers tap setting* T*_{i}, and reactive power of capacitors . Thus, control variables included in each point* d* (solution* d*) are as follows:

The first step is to randomly produce a set of points so that control variables must be within their lower limit and upper limit. The following model expresses initialization step.where* CV*_{min} and* CV*_{max} are respectively lower limit and upper limit of control variables shown in (22).

##### 4.2. Construction of Solution Quality Evaluation Function

Solution quality evaluation function is normally called fitness function when applying metaheuristic algorithms. Basically, the function covers two main terms that are objective function and penalty function wherein objective function is in charge of evaluating effectiveness of control variable regarding considered objective and penalty function undertakes evaluating violation of dependent variables. In problem formulation section, there are two equality constraints and six other inequality constraints; however, not all of these constraints are included in fitness function for evaluating quality. In fact,* V*_{Gi},* T*_{i} together with* Q*_{c} are control variables, so they are always verified and corrected within their limits. On the other hand, formulas (5) and (6) are equality constraints, which have already been taken into account and solved by Matpower program. Thus, inequality constraints (8), (9), and (11) and equality constraints (5) and (6) are not terms in fitness function as a result. Predetermined active power of all generators and control variables are added to Matpower program for running, and then dependent variables consisting of reactive power of generator* Q*_{Gi}, voltage of load buses* V*_{loadi}, and apparent power flow of transmission lines* S*_{l} are determined in addition to three single objectives including power loss, total voltage deviation, and L-index. Finally, quality evaluation function of solution* d* or point* d* is expressed by where* K*_{Q},* K*_{V}, and* K*_{S} are penalty factors corresponding to violation of generator reactive power, voltage of load buses, and apparent power flow of transmission lines;* ∆Q*_{Gi,d},* ∆V*_{loadi,d}, and* ∆S*_{l,d} are violation of reactive power of generator* i*, voltage of load bus* i*, and apparent power flow of transmission line* l *corresponding to solution* d*. Normally, the selection of penalty factors is based on experiment while violation objects are exactly determined by the set of following models:

##### 4.3. Stopping Condition for Iterative Algorithm

The whole search process of the proposed method is stopped since the predetermined maximum iteration number* Max*_{Iter} is reached. At the time, optimal solution for the current run is stored in addition to its evaluation function value. Similarly, after each 100 runs for optimizing one objective function, the best optimal solution is determined and then results are analyzed based on the best solution, the mean solution, the worst solution as well as standard deviation of all solutions.

##### 4.4. The Entire Implementation of MSFS for ORPD Problem

The implementation of the proposed method for finding solutions of ORPD problem can be presented in the following.

*Algorithm 2. *The application of MSFS for solving ORPD problem:*Step 1.* Determine control parameters consisting of population size* N*_{pop}, the number of diffusion* N*_{d}, the maximum number of iterations* Max*_{Iter} and* P*_{a}.*Step 2.* Produce a set of initial points as shown in (23).

Step *3*. Start running Matpower 4.1 for finding dependent variables. *Step 4.* Calculate solution quality evaluation function value as (24)-(27).(i)Determine solution with the lowest value and name it the best point* Point*_{best}.(ii)Start loop search algorithm by setting current iteration number* CI* to 1. *Step 5.* Carry out diffusion process by running (17). *Step 6.* Start running Matpower 4.1 for finding dependent variables. *Step 7.* Calculate solution quality evaluation function value as (24)-(27).(i)Compare such value of each old solution and that of new solutions to keep only one best solution.*Step 8.* Sort all kept solutions in descending order of solution quality evaluation function value.*Step 9.* Choose the first (*P*_{a}*N*_{pop}) solutions.(i)Carry out the first update process by using Algorithm 1 to produce new solutions of the selected solutions. *Step 10.* Start running Matpower 4.1 for finding dependent variables. *Step 11.* Calculate solution quality evaluation function value as (24)-(27).(i)Compare such value of each old solution and that of its new sole solution to keep better one.*Step 12.* Sort all kept solutions in ascending order of solution quality evaluation function value.*Step 13.* Choose the first (*P*_{a}*N*_{pop}) solutions.(i)Carry out the second update process by using Algorithm 1 to produce new solutions of the selected solutions. *Step 14.* Start running Matpower 4.1 for finding dependent variables. *Step 15.* Calculate solution quality evaluation function value as (24)-(27).(i)Compare such value of each old solution and that of new solution to keep better one.*Step 16.* In all current solutions, determine one solution with the lowest solution quality evaluation function value and name it the best solution* Point*_{best.}*Step 17.* Check stopping condition: if current iteration equals the maximum iteration number* Max*_{Iter}, stop iterative algorithm. Otherwise, back to step 5 and set* CI*=*CI*+1.

#### 5. Numerical Results

In the section, we apply both SSFS and the proposed MSFS methods for finding solutions of two sample power systems, IEEE 30-bus system and IEEE 118-bus system with three different single objectives, power loss, voltage deviation, and L-index. The two methods are coded on Matlab 2016a and a personal laptop with 2.4 Ghz processor and 4 Gb of ram. For finding the best solution for each objective functions among three considered objectives, each method is run 100 and then results in terms of the best objective, the mean objective, the worst objective, and standard deviation and simulation time on average are also reported. Three main parameters consisting of the maximum number of iterations, the maximum number of diffusion, and the population size for each case are selected as shown in Table 1. Besides, the number of produced solutions per run (*N*_{nspr}) of SSFS and the proposed method is also reported in the table by using the following formulas:

##### 5.1. Testing Performance of SSFS and MSFS by Using Different Random Walks and Control Parameters

In the first subsection, we test performance of SSFS method when applying two different random walks consisting of Gaussian random walk and Levy random walk in which Gaussian random walk model was presented in (13) and Levy random walk was presented in [50]. Besides, we also perform performance investigation of the proposed MSFS method by using different* P*_{a} values. We set four* P*_{a} values consisting of 0.8, 0.6, 0.4, and 0.2.

For implementation of the survey, IEEE 30-bus system and IEEE 118-bus system, which can be reached by referring to [68], are employed. The first system consists of 6 generator buses, 24 load buses, and 41 transmission lines, 9 capacitor banks, and 4 transformers. The second system is constructed by more number of generators, load buses, transmission lines, capacitor banks, and transformers. The number of control variables for the two different systems is also different wherein 19 control variables are included for the first system including 6 voltage values of 6 generators, 9 reactive power output values of 9 shunt capacitor banks, and 4 tap setting values of 4 transformers, and 73 control variables are included for the second system including 54 voltage values of 54 generators, 14 reactive power output values of 14 shunt capacitor banks, and 9 tap setting values of 9 transformers.

###### 5.1.1. Result Survey Obtained from IEEE 30-Bus System

The results obtained by the proposed method with the setting of different* P*_{a} are reported in Table 2. In addition, the results from the proposed method, SSFS with Gaussian random walk (SSFS-Gauss) and SSFS with Levy random walk, are also summarized in Table 3. Bold numbers, which are the best minimum, mean, maximum, and standard deviation values in Table 2, reveal that the proposed method can find the best optimal solutions for all three objectives by setting* P*_{a} =0.6; however, the search stability of* P*_{a} =0.6 is not the highest stable among four values since it cannot reach the best average and standard deviation for power loss, and the best average, mean, maximum, and standard deviation for L-index objective. In spite of the limitation, the best* P*_{a} is still 0.6.

In Table 3, we also make the best minimum, average, maximum, and standard deviation numbers in bold for better observation and performance comparison. It is clearly visible that most bold numbers, especially all the best minimum values for three objectives, are belonging to the proposed MSFS method. Between two SSFS methods, SSFS with Gaussian random walk can get better results for approximately all three cases excluding standard deviation for power loss objective. The improvement level of best solutions obtained by the proposed method over SSFS-Gauss and SSFS-Levy are calculated by using (30) to be 1.52% and 3.17% for power loss optimization, 6.02% and 10.27% for voltage deviation optimization, and 0.35% and 1.11% for L-index optimization. Obviously, the improvement is significant and necessary for operating power system economically and stably because cost saving in $ yearly is a very high amount, and better voltage stabilization can enable loads such as industrial zones and manufactures work effectively.

In addition, we also compare the worst and the second worst minimum numbers of the proposed method in Table 2 to the best minimum value of SSFS-Gauss. The worst minimums for power loss, voltage deviation, and L-index obtained by the proposed method are, respectively, 4.5408, 0.0897, and 0.1253 while those of SSFS-Gauss are, respectively, 4.584, 0.093, and 0.12484, and those of SSFS-Levy are, respectively, 4.6619, 0.0974, and 0.1258. Clearly, the worst performance of the proposed method is also better than the best performance of SSFS-Levy for all cases and better than the best performance of SSFS-Gauss for power loss and voltage deviation although the worst minimum of L-index from the proposed method is higher than the minimum of SSFS-Gauss, but the second worst minimum of the proposed method is equal to such value of SSFS-Gauss. Furthermore, we also plot 100 values of quality solution evaluation function and optimal solution search characteristic of the best run among 100 runs obtained by the three methods for the three single objectives in Figures 4–9. Figures 4, 5, and 6 show that the proposed method can find many solutions with lower evaluation function value than SSFS-Gauss and SSFS-Levy for all three objectives. For fluctuation comparison, the proposed method can get the best stabilization over the two SSFS methods only for voltage deviation objective, but it has high fluctuation for power loss and L-index objectives. The limits of MSFS can be easily understood because it has used smaller* N*_{nspr}, which was shown in Table 1. In fact, SSFS methods have employed approximately* N*_{nspr}=2,000, but the proposed method has used* N*_{nspr} =1,500 for power loss and L-index optimization cases. It took the proposed method around 3 seconds but SSFS method around 4 seconds. On the other hand, the performance of applied methods is always evaluated by quality of the best optimal solutions and capability of finding global optimum or near global optimum. SSFS-Gauss can reach lower fluctuation than the proposed method, but it cannot find near global optimum; it falls into local optimal zone and it hardly very jumps out the local zones for moving to global optimum or near global optimum. One more time, faster convergence ability to optimum of the proposed method compared to SSFS methods is also observed in Figures 7, 8, and 9.

In summary, the proposed method is superior to SSFS methods in terms of capability of finding higher quality solutions and capability of jumping out search zones with local optimum for moving to promising search space with more effective solutions in addition to faster search ability. Thus, it can be indicated that the proposed modifications in the paper are reasonable for finding solutions of IEEE 30-bus system.

###### 5.1.2. Result Survey Obtained from IEEE 118-Bus System

In order to execute the proposed method for the system,* P*_{a}=0.6, which was the best performance value for the first system above, is used. The results obtained by the three methods are reported in Table 4 for comparison. In the table, the best minimum, average, maximum, and standard deviation numbers are also remarked in bold for better observation. Bold numbers indicate that the proposed MSFS method can find the best optimal solutions with the lowest values of objective function for three optimization cases. Furthermore, the proposed method also reaches the best average, maximum, and standard deviation for power loss optimization case and voltage deviation optimization case in addition to the best mean for L-index optimization case. Clearly, only for L-index optimization case, the proposed method gets the second best maximum and the worst standard deviation. Compared to SSFS-Gauss and SSFS-Levy, the minimum reduction of L-index from the proposed method is not much but that for power loss and voltage deviation is significant. Namely, the improvement level of the proposed method over SSFS-Gauss and SSFS-Levy are, respectively, 3.31% and 5.7% for power loss optimization, 13.08% and 24.8% for voltage deviation optimization, and 0.05% and 0.082% for L-index optimization. In addition to tables and numbers, Figures 10–15 are also used as evidence for proving the superiority of the proposed MSFS over SSFS methods. Unlike Figures 4–6, Figures 10–12 confirm outstanding performance of the proposed MSFS compared to SSFS methods. The figures show that all runs of the proposed method for power loss and voltage deviation optimization cases can find better solutions than both SSFS-Gauss and SSFS-Levy. Furthermore, the proposed method also copes with very low fluctuations for the two cases while that of SSFS methods is much worse. Although results of the proposed method for L-index optimization shown in Figure 12 are not as good as those for two other cases, the proposed method can find many solutions with lower evaluation function value and its fluctuation is also smaller. Figures 13–15 can indicate that solutions of the proposed method at the iteration are much more effective than those of SSFS methods at the last iteration, the iteration. On the other hand, it is also noted that the proposed method used* N*_{nspr}=9,000, but SSFS methods used approximately* N*_{nspr} =12,000. The simulation execution time of the proposed method was around 60 seconds while that of SSFS method was around 80 seconds. It is clearly visible that the proposed method is much more effective and robust than SSFS methods for the IEEE 118-bus system.

Optimal solutions found by the proposed method are reported in Appendix.

##### 5.2. Comparison of Results Obtained from IEEE 30-Bus Power System

In this subsection, the best results of the proposed MSFS method in terms of the best minimum, average, maximum and standard deviation in addition to other important parameters such as maximum number of iterations* Max*_{Iter} and number of newly produced solutions per iteration* N*_{spi} are compared to other methods. Tables 5, 6, and 7 summarize comparisons of power loss, voltage deviation, and L-index obtained by different methods. For better comparison of the minimum, we also use (30) to calculate improvement level (ImLe) of the proposed method over other compared methods. For observation from ImLe, it can be seen that some values are represented by symbol “+”and some values are represented by symbol “-”. The two symbols have different meaning and different evaluation of the proposed method performance. Symbol “+” means that the proposed method has found better solution than such compared methods while symbol “-” points out that such compared methods could find better solutions than the proposed method. However, there have been some cases reporting solutions with constraint violation or incorrect objective. Thus, for the case we substitute reported solutions of compared methods for running Matpower program and check solution as well as real objectives. In Table 5, the symbol “-” can be seen for comparison with FAHCLPSO [30] and MFOT [46] wherein the minimum of the two methods is respectively 4.4877 MW and 4.5128 MW corresponding to ImLe of -0.59% and -0.03%. We have substituted optimal solutions reported in [30, 46] to Matpower for running and obtaining all remaining variables. The recalculated power loss of FAHCLPSO is 6.8372MW but that of MFOT is still 4.5128 MW. Thus, it indicates that the proposed method found better solution than FAHCLPSO but less effective optimal solution than MFOT for the case of minimizing power loss. However, MFOT is not the strongest method because it has used* Max*_{Iter}*=*300 and* N*_{psi}*=*30 while the proposed method has employed* Max*_{Iter}*=*50 and* N*_{psi}*=*30. Clearly, MFOT had to take approximately 6 times execution time of the proposed method. For further investigation of performance in comparison with MFOT, we also set* Max*_{Iter}*=*300 and* N*_{psi}*=*30 for running the proposed method and results obtained are also reported in the last row in the table. The minimum obtained by the setting is also 4.5128 MW equaling that of MFOT method. Thus, it can conclude that MFOT is not more effective than the proposed method for the case. As compared to other methods with symbol “+” of ImLe, it is clear that the number of such methods is high. They are GA variants, PSO variants, DE variants, TLBO variants, standard GSA, and standard ALO. The improvement level of the proposed method over these methods can be up to 7.37%. Furthermore, the proposed method is also the fastest tool because it uses the smallest values of* Max*_{Iter}= 50 and* N*_{spi} = 30 while those from others are much higher, especially the comparison with MDE in [21] using* Max*_{Iter}*=*1500 and* N*_{psi}*=*15, PSO variants in [28] using* Max*_{Iter}*=*200 and* N*_{psi}*=*30, and GSA in [32] using* Max*_{Iter}*=*200 and* N*_{psi}*=*100. In summary, for the case of power loss optimization the proposed method can find better optimal solution than all methods and it is also the fastest optimization tool using the smallest number of iterations and the smallest population size.

For the voltage deviation minimization reported in Table 6, it should be noted that FAHCLPSO [30], GSA [32], and QOTLBO [36] have reported very good voltage deviation corresponding to ImLe of -102.55%, -29.29%, and -2.1%. However, as substituting their optimal solution for running Matpower, we have obtained 1.3442, 0.1862, and 0.1031 for FAHCLPSO, GSA, and QOTLBO, respectively. Clearly, the true values are much higher than real value of the proposed method, which is 0.0874. For comparison with other methods, the proposed method could find better results with faster manner. Namely, the improvement level of the proposed method over other ones can be from approximately 1% to 57% and its search process is very fast with* Max*_{Iter}*=*100 and* N*_{psi}*=*30, but those from others are 100 and 50, 1500 and 15, 200, and 20, 100 and 100. In summary, for voltage deviation objective, the proposed method is the best one for finding optimal solution with the fastest speed.

For comparison with L-index obtained by other methods shown in Table 7, the proposed method can find better optimal solution with improvement level up to +17.01% excluding the comparison with IPGS-PSO in [28], GSA in [32], QOTLBO in [36], and BA, GWO, ABC and ALO in [45]. For more exact comparison, we have also checked solutions and recalculated real L-index of these methods. As a result, L-index of optimal solutions reported by IPGS-PSO and QOTLBO in [36] remains unchanged and all dependent variables from IPGS-PSO are within lower and upper limits but reactive power of generator 4 of QOTLBO is 79.1274 MVAR, which is higher than upper limit of 60 MVAR. On the contrary, recalculated L-index obtained by other methods changes and is much higher than reported values. Namely, those of GSA, ALO, BA, GWO, and ABC are, respectively, 0.1247, 0.1244, 0.1264, 0.1257, and 0.1264. Clearly, only IPGS-PSO has better solution than the proposed method while ALO in [45] has the same solution quality as the proposed method. However, IPGS-PSO and ALO have used very high values of* Max*_{Iter} and* N*_{spi}, namely, 200 and 20 for IPGS-PSO, and 100 and 40 for ALO. Clearly, ALO is slower than the proposed method for converging to the solution with the same quality. For further comparison with IPGS-PSO in [28], we set* Max*_{Iter} and* N*_{spi} to 130 and 30 for running the proposed method and obtained results reported in the last row of Table 7 indicate that the proposed method can find approximate quality solution with IPGS-PSO but mean and standard deviation from the proposed method are lower than those from IPGS-PSO.

In summary, the proposed method can find better results than most methods for the case excluding IPGS-PSO with better solution and ALO with the same solution quality, but the proposed method has faster convergence than all methods excluding IPGS-PSO. Thus, we can conclude that the proposed method is very effective for the case of L-index optimization.

##### 5.3. Comparison of Result Obtained from IEEE 118-Bus Power System

In comparison with other methods for IEEE 118-bus system, obtained results for finding solutions of power loss optimization, voltage deviation and L-index are included in Tables 8, 9, and 10, respectively. The obtained result indicates that the proposed method can find better results than most methods but its results are also worse than several ones. Namely, for power loss optimization reported in Table 8, the proposed method has lower minimum, mean, and maximum power loss than two PSO methods in [27] and eight PSO methods in [28], TLBO in [36], HEP in [65] and MFOT in [46] corresponding to the improvement level from 0.04 to 0.13%. But the proposed method has higher power loss than QOTLBO in [36] and SARCGA in [65] corresponding to the ImLe of -0.02% and -0.01%. The verification of available solutions shows that QOTLBO in [36] reported correct power loss and valid solution, but there was no solution given in [65] for SARCGA. For the case of optimizing voltage deviation in Table 9, the proposed method is greatly superior to all methods including ten PSO methods in both [27, 28], and TLBO and QOTLBO in [36] corresponding to the ImLe from 0.1% to 0.93%.

For L-index comparison in Table 10, the proposed method has better optimal solution quality than two PSO methods in [27], four PSO methods in [28], and two TLBO methods in [36], but it gets worse solutions than three other PSO methods in [28] and approximate solution with another PSO method in [28]. The best performance of the proposed method is more effective than the worst method with the improvement level of 0.57% while the best method IGP-PSO can improve result over the proposed method only with 0.07%. For comparison with other methods in terms of convergence speed, it can be noted that the proposed method is always faster than two PSO methods in [27] for all cases, MFOT in [46] for power loss optimization case and two TLBO methods in [36] for the last two cases because the proposed method has used less values of control parameters than these methods, but it could get better results than these methods. Namely, the numbers of iterations and produced solutions per iteration of the proposed method are 200 and 45 but those are 200 and 120 for methods in [27], 1,000 and 30 for MFOT in [46], and 100 and 100 for TLBO methods in [36]. For comparison with methods in [28], the proposed method has used equal number of iterations but slightly higher number of produced solutions, namely, 45 compared to 40 but it can improve results significantly for power loss and voltage deviation optimization cases and worse result for L-index optimization. Thus, the proposed method is slower to be convergent to high quality solution than three PSO methods in [28] for L-index optimization case.

In summary, the proposed method can find better results than approximately all methods excluding a few ones and it also converges to better results faster by using smaller number of (*Max*_{Itr}*N*_{spi}). Thus, it can conclude that the proposed method is very efficient for searching optimal solutions of ORPD problem. Optimal solutions corresponding to their objective are given in Appendix.

#### 6. Conclusions

This paper has proposed an effective improved version of standard stochastic fractal search algorithm for seeking optimal solutions of ORPD problem considering two standard IEEE systems as study cases. The proposed MSFS method is a metaheuristics based powerful optimization search tool with applying modifications on three generations of producing new solutions in the iterative technique. In the first generation, the proposed method has chosen only one most effective formula among two available ones for producing high quality solutions. Similarly, the second and the third generations have also been carried out by keeping a more effective formula and abandoning a less effective formula. In addition, the kept formula in the last generations has been modified with intent to exploit performance of highly efficient space and avoid spending time searching ineffective zones. By using the proposed modification, the proposed method could reduce a number of computations such as canceling three formulas in three generations, canceling selection of parameter “walk” and comparison steps between random number and 0.5 in the second and the third generations. On the other hand, we have also tried the performance of Gaussian random walk and Levy random walk for SSFS, and we have concluded that the proposed method needed Gaussian random walk rather Levy random walk. The obtained result comparisons with SSFS have given persuasive evidences for conclusion that the proposed method could find much better solutions, owned more stable search capability and converged to better results with faster speed. Further evaluation of the proposed method has been done by comparing with a high number of existing methods including PSO family, DE family, GA family, TLBO family, and other recent methods. Comparisons have indicated that the proposed method could find solutions with less power loss, voltage deviation, and L-index than those from nearly all methods. Furthermore, search speed of the proposed method was also faster than all methods excluding a few improved versions of PSO for L-index optimization case. As a result, it can be pointed out that the proposed method deserves to be applied for ORPD problem as one of the most effective optimization tools.

#### Appendix

#### Nomenclature

: | The fraction of the best solutions selected in the second update |

: | Maximum number of iterations |

: | Number of buses in the considered power system |

: | Number of buses with the presence of capacitor banks |

: | Maximum number of diffusion |

: | Number of generators with the presence of generators |

: | Number of transmission branches |

: | Number of load buses |

: | Population size or a set of points |

: | Number of new solutions produced per iteration |

: | Number of buses with transformer |

: | The fraction of the worst solutions selected in the first update |

: | Active and reactive power of load installed at bus , respectively |

: | Total power loss in all transmission branches |

: | Randomly picked points in the current population |

: | Reactive power output of capacitor banks installed at bus |

: | Lower and upper reactive power output limitations of capacitor bank installed at bus , respectively |

: | Lower and upper reactive power output limitations of generator at bus , respectively |

: | The rank of solution |

: | Maximum apparent power flow of branch |

: | Lower and upper tap changer limitations of all transformers, respectively |

: | Total voltage deviation of all load buses |

: | Lower and upper voltage magnitude limitations of all generators, respectively |

: | Working voltage magnitude of buses and , respectively |

: | Lower and upper voltage magnitude limitations of all loads, respectively |

: | Two random number within 0 and 1 of , respectively |

: | Real and imaginary terms of admittance of branch , which is connecting buses and , respectively |

: | Angles of voltage at buses and , respectively. |

#### Data Availability

No data were used to support this study. All data are provided in full in the results section of this paper.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.