Research Article  Open Access
Thang Trung Nguyen, Dieu Ngoc Vo, Hai Van Tran, Le Van Dai, "Optimal Dispatch of Reactive Power Using Modified Stochastic Fractal Search Algorithm", Complexity, vol. 2019, Article ID 4670820, 28 pages, 2019. https://doi.org/10.1155/2019/4670820
Optimal Dispatch of Reactive Power Using Modified Stochastic Fractal Search Algorithm
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 Lindex 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 30bus system and IEEE 118bus 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 PSObased method group, GAbased method group, DEbased 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, onload 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 loopgenetic 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], quasioppositional differential evolution (QODE) [22], multiobjective differential evolution (MODE) [23], a novel hybrid chaotic artificial bee colony, and differential evolution algorithm (CABCDE) [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 timevarying 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 30bus system and compared to three other DE methods such as conventional DE, selfadaptive 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. CABCDE 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). CABCDE has been compared to the standard DE via the results obtained from IEEE 14bus system and IEEE 30bus 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], comprehensivelearning particle swarm optimization (CLPSO) [27], improved pseudogradient searchparticle swarm optimization (IPGSPSO) [28], hybrid particle swarm optimization and imperialist competitive algorithms (HPSOICA) [29], fuzzy adaptive heterogeneous comprehensivelearning particle swarm optimization (FAHCLPSO) [30], and hybrid cuckoo search and particle swarm optimization (CSPSO) [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 searchparticle swarm optimization (IPGSPSO), 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, IPGSPSO 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 timevarying inertia weight (PSOTVIW), PSO with timevarying acceleration coefficients (PSOTVAC), selforganizing PSO with TVAC (SPSOTVAC), PSO with constriction factor PSOCF, pseudogradient search PSO (PGSPSO), stochastic weight tradeoff PSO (SWTPSO), and pseudogradient search and SWTPSO (PGSWTPSO). HPSOICA was a combination of the imperialist competitive algorithm and improved PSO. In HPSOICA, 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 ComprehensiveLearning Particle Swarm Optimization (HCLPSO) aiming to avoid trapping in local optimal zones. Via the comparison of results obtained from IEEE 30bus system, IEEE 118bus system, and IEEE 354bus 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 (IGSAFCSS) [35], quasioppositional 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 (HNMSFA) [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], mothflame 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 twentythree 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 (CMAES), 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 walkbased diffusion (SFSDE) [57], hybrid stochastic fractal search and pattern search technique (HSFSPS) [58], hybrid SFS based on different mutation models of DE (HSFSDE) [59], improved SFS (ISFS) [60], and chaosbased SFS (CSFS) [61]. In general, almost all improved methods have used different mutation models of DE available in [62]. SFSDE has applied currenttobest mutation model for solutions with a probability of 50%. If random number was lower than 0.5, the currenttobest 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. HSFSPS 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. HSFSPS has been employed for designing different controllers based on PI, PID, and cascaded PIPD. Final evaluation in the study was that cascaded PIPD controller, which was yielded by applying HSFSPS, was superior to remaining controllers applying HSFSPS. So, the performance of HSFSPS is now still a question for readers because it was not proven to be better than any optimization algorithms, even SSFS and PS. HSFSDE in [59] has abandoned original diffusion process of SSFS, and added mutation technique based on targettobest/1 in [62] for new diffusion process. In addition, selfadaptive 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. HSFSDE 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 HSFSDE 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 33bus 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 69bus system and 118bus 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 (SSFSGauu) and SSFS with Levy random walk (SSFSLevy) 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 Lindex. 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 Lindex: 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 Lindex at each bus in power system. The enhancement of voltage stability is carried out by reducing the highest value of Lindex at one bus in the power system. Namely, Lindex at each bus i and the voltage stability enhancement are mathematically formulated by [11]:where L_{i} is the Lindex 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 Lindex. 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 30bus system and IEEE 118bus system with three different single objectives, power loss, voltage deviation, and Lindex. 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 30bus system and IEEE 118bus 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 30Bus 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 (SSFSGauss) 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 Lindex 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 SSFSGauss and SSFSLevy 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 Lindex 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 SSFSGauss. The worst minimums for power loss, voltage deviation, and Lindex obtained by the proposed method are, respectively, 4.5408, 0.0897, and 0.1253 while those of SSFSGauss are, respectively, 4.584, 0.093, and 0.12484, and those of SSFSLevy 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 SSFSLevy for all cases and better than the best performance of SSFSGauss for power loss and voltage deviation although the worst minimum of Lindex from the proposed method is higher than the minimum of SSFSGauss, but the second worst minimum of the proposed method is equal to such value of SSFSGauss. 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 SSFSGauss and SSFSLevy 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 Lindex 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 Lindex 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. SSFSGauss 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 30bus system.
5.1.2. Result Survey Obtained from IEEE 118Bus 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 Lindex optimization case. Clearly, only for Lindex optimization case, the proposed method gets the second best maximum and the worst standard deviation. Compared to SSFSGauss and SSFSLevy, the minimum reduction of Lindex 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 SSFSGauss and SSFSLevy 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 Lindex 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 SSFSGauss and SSFSLevy. 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 Lindex 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 118bus system.

Optimal solutions found by the proposed method are reported in Appendix.
5.2. Comparison of Results Obtained from IEEE 30Bus 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 Lindex 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 Lindex 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 IPGSPSO 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 Lindex of these methods. As a result, Lindex of optimal solutions reported by IPGSPSO and QOTLBO in [36] remains unchanged and all dependent variables from IPGSPSO 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 Lindex 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 IPGSPSO has better solution than the proposed method while ALO in [45] has the same solution quality as the proposed method. However, IPGSPSO and ALO have used very high values of Max_{Iter} and N_{spi}, namely, 200 and 20 for IPGSPSO, 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 IPGSPSO 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 IPGSPSO but mean and standard deviation from the proposed method are lower than those from IPGSPSO.
In summary, the proposed method can find better results than most methods for the case excluding IPGSPSO with better solution and ALO with the same solution quality, but the proposed method has faster convergence than all methods excluding IPGSPSO. Thus, we can conclude that the proposed method is very effective for the case of Lindex optimization.
5.3. Comparison of Result Obtained from IEEE 118Bus Power System
In comparison with other methods for IEEE 118bus system, obtained results for finding solutions of power loss optimization, voltage deviation and Lindex 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 Lindex 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 IGPPSO 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 Lindex optimization. Thus, the proposed method is slower to be convergent to high quality solution than three PSO methods in [28] for Lindex 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 Lindex 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 Lindex 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.
References
 W. VillaAcevedo, J. LópezLezama, and J. ValenciaVelásquez, “A novel constraint handling approach for the optimal reactive power dispatch problem,” Energies, vol. 11, no. 9, p. 2352, 2018. View at: Publisher Site  Google Scholar
 T. Nguyen, D. Vo, N. Vu Quynh, and L. Van Dai, “Modified cuckoo search algorithm: a novel method to minimize the fuel cost,” Energies, vol. 11, no. 6, p. 1328, 2018. View at: Publisher Site  Google Scholar
 D. C. Yu, J. E. Fagan, B. Foote, and A. A. Aly, “An optimal load flow study by the generalized reduced gradient approach,” Electric Power Systems Research, vol. 10, no. 1, pp. 47–53, 1986. View at: Publisher Site  Google Scholar
 D. I. Sun, B. Ashley, B. Brewer, A. Hughes, and W. F. Tinney, “Optimal power flow by Newton approach,” IEEE Transactions on Power Apparatus and Systems, vol. 103, no. 10, pp. 2864–2880, 1984. View at: Publisher Site  Google Scholar
 S. Granville, “Optimal reactive dispatch through interior point methods,” IEEE Transactions on Power Systems, vol. 9, no. 1, pp. 136–146, 1994. View at: Publisher Site  Google Scholar
 E. Rezania and S. M. Shahidehpour, “Real power loss minimization using interior point method,” International Journal of Electrical Power & Energy Systems, vol. 23, no. 1, pp. 45–56, 2001. View at: Publisher Site  Google Scholar
 J. Zhu and X. Xiong, “Optimal reactive power control using modified interior point method,” Electric Power Systems Research, vol. 66, no. 2, pp. 187–192, 2003. View at: Publisher Site  Google Scholar
 O. Alsaç, J. Bright, M. Prais, and B. Stott, “Further developments in LPbased optimal power flow,” IEEE Transactions on Power Systems, vol. 5, no. 3, pp. 697–711, 1990. View at: Publisher Site  Google Scholar
 B. Stott and J. L. Marinho, “Linear programming for powersystem network security applications,” IEEE Transactions on Power Apparatus and Systems, vol. 98, no. 3, pp. 837–848, 1979. View at: Publisher Site  Google Scholar
 A. M. Chebbo and M. R. Irving, “Combined active and reactive dispatch. i. problem formulation and solution algorithm,” IEE ProceedingsGeneration, Transmission and Distribution, vol. 142, no. 4, pp. 393–400, 1995. View at: Publisher Site  Google Scholar
 F.C. Lu and Y.Y. Hsu, “Reactive power/voltage control in a distribution substation using dynamic programming,” IEE Proceedings Generation, Transmission and Distribution, vol. 142, no. 6, pp. 639–645, 1995. View at: Publisher Site  Google Scholar
 V. H. Quintana and M. SantosNieto, “Reactivepower dispatch by successive quadratic programming,” IEEE Transactions on Energy Conversion, vol. 4, no. 3, pp. 425–435, 1989. View at: Publisher Site  Google Scholar
 N. Grudinin, “Reactive power optimization using successive quadratic programming method,” IEEE Transactions on Power Systems, vol. 13, no. 4, pp. 1219–1225, 1998. View at: Publisher Site  Google Scholar
 V. de Sousa, E. Baptista, and G. da Costa, “Optimal reactive power flow via the modified barrier Lagrangian function approach,” Electric Power Systems Research, vol. 84, no. 1, pp. 159–164, 2012. View at: Publisher Site  Google Scholar
 Y. J. Cao and Q. H. Wu, “Optimal reactive power dispatch using an adaptive genetic algorithm,” in Proceedings of the Second International Conference on Genetic Algorithms in Engineering Systems, pp. 117–122, Glasgow, UK, 1997. View at: Publisher Site  Google Scholar
 D. T. Khanmiri, N. Nasiri, and T. Abedinzadeh, “Optimal reactive power dispatch using an improved genetic algorithm,” International Journal of Computer and Electrical Engineering, vol. 4, no. 4, pp. 463–466, 2012. View at: Publisher Site  Google Scholar
 B. B. Pal, P. Biswas, and A. Mukhopadhyay, “GA based FGP approach for optimal reactive power dispatch,” Procedia Technology, vol. 10, pp. 464–473, 2013. View at: Publisher Site  Google Scholar
 M. S. Alam and M. De, “Optimal reactive power dispatch using hybrid loopgenetic based algorithm,” in Proceedings of the 2016 19th National Power Systems Conference (NPSC), pp. 1–6, Bhubaneswar, India, December 2016. View at: Publisher Site  Google Scholar
 K. Rayudu, G. Yesuratnam, and A. Jayalaxmi, “Improving voltage stability by optimal reactive power dispatch based on genetic algorithm and linear programming technique,” in Proceedings of the 2016 International Conference on Electrical, Electronics, and Optimization Techniques, ICEEOT 2016, pp. 1357–1362, March 2016. View at: Google Scholar
 A. A. A. E. Ela, M. A. Abido, and S. R. Spea, “Differential evolution algorithm for optimal reactive power dispatch,” Electric Power Systems Research, vol. 81, no. 2, pp. 458–464, 2011. View at: Publisher Site  Google Scholar
 H. Singh and L. Srivastava, “Modified differential evolution algorithm for multiobjective VAR management,” International Journal of Electrical Power & Energy Systems, vol. 55, pp. 731–740, 2014. View at: Publisher Site  Google Scholar
 M. Basu, “Quasioppositional differential evolution for optimal reactive power dispatch,” International Journal of Electrical Power & Energy Systems, vol. 78, pp. 29–40, 2016. View at: Publisher Site  Google Scholar
 M. Basu, “Multiobjective optimal reactive power dispatch using multiobjective differential evolution,” International Journal of Electrical Power & Energy Systems, vol. 82, pp. 213–224, 2016. View at: Publisher Site  Google Scholar
 Y. Li, X. Li, and Z. Li, “Reactive power optimization using hybrid CABCDE algorithm,” Electric Power Components and Systems, vol. 45, no. 9, pp. 980–989, 2017. View at: Publisher Site  Google Scholar
 M. Kumari and M. Sydulu, “Improved Particle Swarm Algorithm Applied to Optimal Reactive Power Control,” in Proceedings of the 2006 IEEE International Conference on Industrial Technology, pp. 1873–1878, Mumbai, India, December 2006. View at: Publisher Site  Google Scholar
 G. Cai, Z. Ren, and T. Yu, “Optimal reactive power dispatch based on modified particle swarm optimization considering voltage stability,” in Proceedings of the 2007 IEEE Power Engineering Society General Meeting, pp. 1–5, Tampa, FL, USA, June 2007. View at: Publisher Site  Google Scholar
 K. Mahadevan and P. S. Kannan, “Comprehensive learning particle swarm optimization for reactive power dispatch,” Applied Soft Computing, vol. 10, no. 2, pp. 641–652, 2010. View at: Publisher Site  Google Scholar
 J. Polprasert, W. Ongsakul, and V. N. Dieu, “Optimal reactive power dispatch using improved pseudogradient search particle swarm optimization,” Electric Power Components and Systems, vol. 44, no. 5, pp. 518–532, 2016. View at: Google Scholar
 M. Mehdinejad, B. MohammadiIvatloo, R. DadashzadehBonab, and K. Zare, “Solution of optimal reactive power dispatch of power systems using hybrid particle swarm optimization and imperialist competitive algorithms,” International Journal of Electrical Power & Energy Systems, vol. 83, pp. 104–116, 2016. View at: Publisher Site  Google Scholar
 E. Naderi, H. Narimani, M. Fathi, and M. R. Narimani, “A novel fuzzy adaptive configuration of particle swarm optimization to solve largescale optimal reactive power dispatch,” Applied Soft Computing, vol. 53, pp. 441–456, 2017. View at: Publisher Site  Google Scholar
 J. Ding, Q. Zhang, and Y. Ma, “Optimal reactive power dispatch based on the CSPSO algorithm,” in Proceedings of the 13th Conference on Industrial Electronics and Applications (ICIEA), IEEE, Wuhan, China, 2018. View at: Publisher Site  Google Scholar
 S. Duman, Y. Sönmez, U. Güvenç, and N. Yörükeren, “Optimal reactive power dispatch using a gravitational search algorithm,” IET Generation, Transmission & Distribution, vol. 6, no. 6, pp. 563–576, 2012. View at: Publisher Site  Google Scholar
 P. K. Roy, B. Mandal, and K. Bhattacharya, “Gravitational search algorithm based optimal reactive power dispatch for voltage stability enhancement,” Electric Power Components and Systems, vol. 40, no. 9, pp. 956–976, 2012. View at: Publisher Site  Google Scholar
 S. Duman, Y. Sonmez, U. Guvenc, and N. Yorukeren, “Application of gravitational search algorithm for optimal reactive power dispatch problem,” in Proceedings of the International Symposium on INnovations in Intelligent SysTems and Applications (INISTA '11), pp. 519–523, June 2011. View at: Publisher Site  Google Scholar
 G. Chen, L. Liu, Z. Zhang, and S. Huang, “Optimal reactive power dispatch by improved GSAbased algorithm with the novel strategies to handle constraints,” Applied Soft Computing, vol. 50, pp. 58–70, 2017. View at: Publisher Site  Google Scholar
 B. Mandal and P. K. Roy, “Optimal reactive power dispatch using quasioppositional teaching learning based optimization,” International Journal of Electrical Power & Energy Systems, vol. 53, pp. 123–134, 2013. View at: Publisher Site  Google Scholar
 M. Ghasemi, M. Taghizadeh, S. Ghavidel, J. Aghaei, and A. Abbasian, “Solving optimal reactive power dispatch problem using a novel teachinglearningbased optimization algorithm,” Engineering Applications of Artificial Intelligence, vol. 39, pp. 100–108, 2015. View at: Publisher Site  Google Scholar
 A. Rajan and T. Malakar, “Optimal reactive power dispatch using hybrid Nelder–Mead simplex based firefly algorithm,” International Journal of Electrical Power & Energy Systems, vol. 66, pp. 9–24, 2015. View at: Publisher Site  Google Scholar
 K. Rayudu, G. Yesuratnam, and A. Jayalaxmi, “Artificial Bee Colony algorithm for optimal reactive power dispatch to improve voltage stability,” in Proceedings of the 2016 IEEE International Conference on Circuit, Power and Computing Technologies, ICCPCT 2016, pp. 1–7, March 2016. View at: Google Scholar
 K. Abaci and V. Yamaçli, “Optimal reactivepower dispatch using differential search algorithm,” Electrical Engineering, vol. 99, no. 1, pp. 213–225, 2017. View at: Publisher Site  Google Scholar
 A. Rajan and T. Malakar, “Exchange market algorithm based optimum reactive power dispatch,” Applied Soft Computing, vol. 43, pp. 320–336, 2016. View at: Publisher Site  Google Scholar
 A. Mukherjee and V. Mukherjee, “Chaos embedded krill herd algorithm for optimal VAR dispatch problem of power system,” International Journal of Electrical Power & Energy Systems, vol. 82, pp. 37–48, 2016. View at: Google Scholar
 M. H. Sulaiman, Z. Mustaffa, M. R. Mohamed, and O. Aliman, “Using the gray wolf optimizer for solving optimal reactive power dispatch problem,” Applied Soft Computing, vol. 32, pp. 286–292, 2015. View at: Publisher Site  Google Scholar
 A. A. Heidari, R. Ali Abbaspour, and A. Rezaee Jordehi, “Gaussian barebones water cycle algorithm for optimal reactive power dispatch in electrical power systems,” Applied Soft Computing, vol. 57, pp. 657–671, 2017. View at: Publisher Site  Google Scholar
 S. Mouassa, T. Bouktir, and A. Salhi, “Ant lion optimizer for solving optimal reactive p