Optimal Reactive Power Generation for Transmission Power Systems Considering Discrete Values of Capacitors and Tap Changers

: In this paper, an improved coyote optimization algorithm (ICOA) is developed for determining control parameters of transmission power networks to deal with an optimal reactive power dispatch (ORPD) problem. The performance of ICOA method is superior to its conventional coyote optimization algorithm (COA) thanks to modiﬁcations of two new solution generations of COA. COA uses a center solution to generate an update step size in the ﬁrst solution generation and produced one new solution by using random factors to diversify the search space in the second solution generation. By tackling the drawbacks of COA, ICOA can reduce control parameters and computation steps, shorten execution time, and provide better results. ICOA is compared to its conventional COA for three standard IEEE systems of 30-, 57, and 118-buses with continuous and discrete control variables. Moreover, three other algorithms such as water cycle algorithm (WCA), salp swarm algorithm (SSA), and sunﬂower optimization algorithm (SFOA) have been also implemented for further investigation of the real performance of the proposed method. All the applied methods are metaheuristic algorithms based on population and randomization. The result comparison from the test systems has indicated that ICOA can provide higher solution quality than other methods with reasonable execution time. Therefore, ICOA is a reliable tool for ﬁnding optimal solutions of the ORPD problem.


Motivation and Incitement
Nowadays, ORPD plays an extremely important role and is acknowledged by the interest of engineers and managers as well as a number of researchers in power system planning and operation areas [1]. The ORPD problem has many complicated constraints and different single objectives [2]. The main target of the ORPD problem is to cut the active power loss on conductors and improve the profile of load voltage, while keeping operating parameters of electric components within allowable ranges. In the ORPD, active power generation of PPs is predetermined, while the voltage of PPs is tuned reasonably. In addition, tap of transformers and generation of capacitors are also tuned for keeping values of reactive power generation of PPs, currents of conductor, and load voltages within allowable operating ranges [3]. Therefore, finding solutions of ORPD is one of the greatest challenges for scientists, and optimization algorithms become believable mathematical tools.
Methods miming the rule of natural selection or heredities, relating to the terms of genetics, crossover, and mutation selection, are called GA. AGA [9] is a reform of the genetic algorithm. SARCGA has been recommended in [10] by using simulated binary crossover, binary tournament selection, and polynomial mutation distribution of RCGA. In [11], HLGA is an incorporation of loop method and GA in which striking global search characteristic of GA has been used to exploit large solution spaces, whilst the loop method has played a task of search within the predetermined zones found by GA. In [13], authors have solved the ORPD problem by applying standard GA with Roulette wheel-based selection mechanism and improved GA with classification selection mechanism. A new different version, called SGA [14], was introduced in 2018 for handling the ORPD problem considering power-loss minimization only. A different optimization technique, which is called particle swarm optimization, is also not less famous than GA. It includes original PSO [15] and lots of variants such as CLPSO [16], ALC-PSO [17], PSO [18], PSO-IPG [19], HPSO [20], HPSO-ICA [21], PSO-GT [22], and HPSO-TS [23]. HMAPSO is a combination of GA, multiagent system and the excellent characteristic of PSO. Therefore, the method is capable of finding solutions better than other implemented methods such as GA, PSO, HPSO, and MPSO. PSO-IPG [19] is more effective than PSO-AIW, PSO-AAC, SPSO-AAC, PSO-CF, PSO-PG, PSO-SWT, and PSO-PGSWT.
A large DE group consisting of standard DE [24], DE-AS [25], MTLA-DDE [26], QODE [27], and IABC-DE [28] has been reapplied or newly recommended for the ORPD problem. The purpose of such methods is to overcome some shortcomings of DE, such as premature convergence to not very high-quality solutions and easily getting stuck into less effective solutions by improving main stages of DE or combination of DE and other methods with more potential characteristics. In DE-AS [25], authors have modified mutation and selection mechanisms of DE by applying VSMF and PST rule used in the ant system. Using both VSMF and PST has supported the method to achieve not only good solution quality but also avoid consuming time manner better than PSO and DE through IEEE 30-bus system with minimizing power-loss objective function.
In general, authors have tried to develop better optimization algorithms for finding better optimal solutions than those of previously publicized algorithms. Approximately all demonstrations have been based on quality of the most optimal solutions and convergence characteristic of the best run among the number of runs. As compared to previous methods shown in earlier studies, algorithms have not been proved persuasively. The settings of control parameters or the comparisons of control parameters have not been clarified or just have just been analyzed sketchily. Some studies were based on execution time to conclude potential search speed of applied methods; however, they have forgotten that the same methods run on different computers would result in much different computational time.

Contributions and Paper Organization
In this paper, the concerned ORPD problem with three different standard IEEE systems consisting of 30-, 57-, and 118-buses was solved by ICOA, which was first developed by replacing two lowly effective techniques of standard COA. COA, a novel populationbased metaheuristic algorithm, was developed by Juliano and Leandro in 2018 [56]. COA achieved a smaller average, especially for hybrid benchmark functions and cases with high number of variables. However, the real performance of COA for engineering problems such ass ORPD has not been demonstrated persuasively, since there has been no application so far. However, COA has not achieved good results for other problems in engineering, and it has been suggested to be modified [57][58][59]. For the cases of solving benchmark functions, many functions have real optimal solutions with zero values, and construction of COA has produced solutions around the real ones. In fact, in the first generation of each iteration, COA produced center solutions, and then the new solution was newly produced based on the center ones. The center ones are around the middle point of upper and lower boundaries. The feature unintentionally supports COA to exploit the search space, wherein real global optimum solutions are existing for cases of benchmark functions. Unfortunately, the feature is not effective for the ORPD problem because the real optimum solutions of the ORPD problem are not around the center ones. Therefore, the proposed ICOA method applied the first modification on the first generation by canceling center solutions and using the most effective solution so far. The second modification has been proposed due to low efficiency of the second generation. In the second generation, COA produced one new solution for each pack based on randomization, but there must be selection condition for the randomization. Thus, the second generation is time consuming, but its effectiveness is not high. Basically, the feature is appropriate for the case that approximately all current solutions are lumped in local optimum zones. It can be supported to refresh almost all solutions except the best one, resulting in a high possibility of explore other zones effectively and quickly rather than thank to nature ability. Therefore, in the second modification, we have used another strategy for producing such new solutions of each pack. Thanks to the application of two modifications, the proposed ICOA approach can solve ORPD more effectively and robustly than COA in terms of better stabilization of found solutions as well as finding more effective solutions faster. In addition, the proposed ICOA approach is also competed with other applied metaheuristic algorithms including original ones, improved versions of original ones together with other hybrid methods mentioned above. The IEEE 30-, 57-, and 118-bus power systems with three single objectives, being the TPL, TVD, and L index, are used for comparison and conclusion. The main contributions of the study can be summarized as follows: • Show the structure of conventional COA and indicate its strong points and weak points. Based on found strong and weak points, good features of COA are retained, but bad features are changed by applying modifications, • Simulate the result for the modified method and its conventional method for all study cases to indicate the impact of modifications on the performance of the proposed method, • Provide a strong optimization approach to get good solutions for the ORPD problem. Power loss, voltage deviation, and voltage stability of the applied systems become better by the solutions provided by the proposed approach, • The proposed method can solve the ORPD problem faster than other compared methods.
In addition to the introduction section, the other parts of the paper are as follows: Description of concerned the ORPD problem is shown in Section 2. The original COA and the proposed ICOA method are explained in detail in Section 3. The whole computation procedure of solving the concerned the ORPD problem by using the proposed ICOA approach is given in Section 4. Section 5 shows the comparison of results obtained by the proposed methods and other methods. Finally, the conclusion is given in Section 6. In addition to references, optimal solutions obtained by the proposed ICOA method are also included at the end of the paper.

Optimal Reactive Power Dispatch Problem
In the ORPD problem, systems with electric components shown in Figure 1 are considered. Due to the given active power generation of power plants, three single objective functions are considered to be minimized, while constraints are strictly imposed on electric components and other factors. Basically, constraints of electric components are represented by inequalities while those of other factors are expressed by equalities. The expression can be observed by using the following mathematical functions:

Objective Functions
Minimization of total active power losses (TPL) in all branches: Total energy losses in transmission power network are significant because current magnitude flowing through conductors are high and the length of transmission lines can be up to hundreds of kilome-ters. In order to reduce the cost of energy losses, the minimization of total power losses is a unique solution so far. The total power losses can be found by the following formula: Minimization of total voltage deviation (TVD): In transmission power network, voltage magnitude of load buses and generator buses are the primary concerns for stable operation of the network. The voltage of generators can be controlled by the adjusting exciter, while voltage of loads is dependent on parameters of other components in transmission power network. The good voltage profile is a target of operation of power systems, and it can be achieved by minimizing total voltage deviation. The total voltage deviation is calculated by the model below: Minimization of voltage stability-L index: In transmission power networks, voltage of load is not stable and fixed at a value due to the change of consumed reactive power and active power. This can lead to different cases of power system status, high voltage, and voltage collapse. Basically, L index is in a range from 0 to 1 in which 0 is the best case and 1 is the worst status of power system. Therefore, the voltage stability enhancement can be achieved if L index is very close to 0 [11]. In order to get the target, L index at each load bus is determined first, and then the highest L index of a bus is selected for minimization. L i index of the ith bus is first calculated by [11]: Then, L index is determined by choosing L i with the highest value (i.e., L index = max(L i )).

A Set of Constraints
The set of constraints imposed on the ORPD is described as follows: Lower bound and upper bound constraints of generators: As explained in Section 2.2, active power output of generators are known and always within lower bound and upper bound. Thus, the control of each generator is only to monitor and adjust reactive power and voltage magnitude within an allowed range. Generators can be working stably if the following rules are always obeyed: V Gen,min ≤ V Gen,i ≤ V Gen,max ; i = 1, . . . , NG Lower bound and upper bound constraints of shunt capacitor banks: At some buses in transmission power network, shunt capacitor banks are additionally installed with intent to supply reactive power to loads and reduce voltage drop in transmission lines. The reactive power output of shunt capacitor banks must follow the rule below: In the study, the value of reactive power output of shunt capacitors is classified into two different cases. It is a continuous variable in the first case as long as it is within lower and upper bounds as shown in Equation (6), but it is discrete variable with a step of 100 KVAR in the second case.
Lower bound and upper bound constraints of tap changer of transformers: During the operation of transformers, voltage of secondary side is always controlled by setting tap changer, while voltage of primary side is dependent on current flowing through transmission lines, voltage of generators, and reactive power of shunt capacitor banks. The setting of tap changer of transformer at the bus must obey the following inequality: Similar to shunt capacitors, tap changer of transformers is also considered to be a continuous variable and discrete variable for two study cases. In the second case, it is constrained by each change step of 0.01 pu.
Capacity constraint of transmission lines: On the contrary to other devices, transmission lines are conductors, which are not imposed on lower bound, but upper bound constraint is strictly supervised. Apparent power flowing through transmission lines should not be higher than its capacity. If the following inequality can be accurately responded, the transmission line is working stably.
Lower bound and upper bound constraints of voltage at load buses: In addition to two objectives regarding voltage such as total voltage deviation and L index, voltage at each load bus is also imposed on lower and upper bounds for the purpose of stable operation of load. The requirement of load bus voltage is as follows: Power balance constraint at each bus: In addition to constraints from electric components in transmission lines, the ORPD problem also considers active and reactive power balance at each bus as a serious concern due to issues regarding frequency and voltage. Power generated by generator and capacitor, power consumed by load, and power of connected transmission line have the following relationship:

Original Coyote Optimization Algorithm
In the COA, Coyote community is equally disunited into N pack packs, and each pack contains N coyote coyotes. Therefore, the population size is the result of multiplication of N pack and N coyote . The initial generation for each pack among the population can be performed by using the following popular model: where: Sol min = x j,min ; j = 1, . . . , NCV Sol max = x j,max ; j = 1, . . . , NCV where Sol d,p is the social condition of the dth coyote in the pth pack corresponding to the ith solution in the pack p. Sol min and Sol max represent upper and lower bounds of the social conditions; x j,min and x j,max are the lowest value and highest value of the jth control variable; and NCV is number of control variables. The initialization step of COA is ended by the two following tasks: calculate quality of each solution Fitness d,p , and determine the best solution in each pack, Sol best,p .
In order to produce new social condition for each coyote, COA produces a center condition for each pack that is similar to a center solution, Sol center,p . The center solution Sol center,p collects center variables in each pack. Thus, there are two cases for each center variable in the center solution due to the number of coyotes in each pack. For the case of odd number of coyotes, the center variable of the coyote (N coyote + 1)/2 is picked up, but for the case of even number of coyotes, the variable of coyote N coyote /2 is taken. After determining the center solution for each pack, solutions are newly updated by the following rule: It is not sure that each variable x new j,d,p of new solution Sol new d,p always exactly satisfies lower limit x j,min and upper limit x j,max . Thus, each variable must be verified and corrected by using the following equation: For replacing the worst social condition in each pack, COA produces a new social condition, and then its quality is evaluated. In case that its quality is more effective than the worst condition, which is being existed in the pack, the worst one is replaced. Otherwise, the worst condition cannot be eliminated unintentionally. The new condition Sol new p is also corresponding to a new solution in which the jth variable x j of Sol new p can be found by the following way: where x new j,r1,p and x new j,r2,p are two variables in two randomly picked solutions; x j,r is the jth variable obtained by randomization within lower and upper bounds. Pro 1 and Pro 2 are, respectively, the scatter and association probabilities and obtained by: Before terminating an iteration of search process, COA performs the action of leaving current pack and joining another pack for coyotes. The phenomenon results in the exchange of solutions among available packs. Each two randomly selected solutions in two random packs are exchanged if the following condition is satisfied: The exchange of solutions among different packs can avoid lumping all solutions into one zone, restricting search space and leading to local optimum easily. However, the effectiveness of the operator is very dependent on the random number and the value of N coyote . If N coyote is less than 10, the possibility that exchange phenomenon occurs is less than 50%. If N coyote is not smaller than 14, the solution exchange certainly takes place.
The whole computation process of COA is expressed in Figure 2.

Improved Coyote Optimization Algorithm
In each iteration, the search process of COA is comprised of two new solution generations, in which the first generation produces N pop solutions by using Equation (15) but the second generation produces only one solution by using Equation (19). Thus, there are (N pop + 1) new solutions for each iteration, and there are [(N pop + 1)*N Iter ] new solutions for one run. The information indicates that if the function of the two equations is ineffective or insignificantly effective, the best optimal solution of COA is of very low quality, and the considered optimization problem has been solved unsuccessfully. Consequently, in the paper, we have investigated the performance of the two methods of producing new solutions and proposed two new methods for the two major new solution generations. The new methods are presented and explained in the following sections.

The First Modification
In the first modification, we focus on the function of Equation (15) with intent to produce new solutions with higher quality. As seen from Equation (15), COA has employed a center solution Sol center,p for each pack p to update new solutions because in natural behavior of coyote, a center tendency or a center social condition of all coyotes can lead to a better new social condition for each coyote in the flock. Each coyote can find new social conditions more effectively by using the center social condition from the whole population. However, center variables are impossible to form a promising solution. Furthermore, the use of the center variables, which are in charge of producing new solutions in Equation (15), is completely ineffective. Center variables included in one solution, and a center social condition are totally different in the nature phenomenon and in math. However, the method was successful in [56] in dealing with benchmark optimization functions, because many functions have optimal solutions with "0 value" variables. It should be noted that the "0 value" solutions are in the middle of the lower bound and upper bound of variables. Unintentionally, center solution is very useful for finding "0 value" solutions. Unlikely, the manner is not useful for the ORPD problem, because the optimal solutions of the ORPD problem are not around the center ones. Consequently, we suggest the model in Equation (15) should be changed by replacing center solution with a global best solution, which has been found so far. The search is performed as follows: The equation can take advantage of the information exchange between current solutions and the so-far best solution. By using the modification, COA can take several advantages as follows: 1.
Reduce determination of N pack center solutions for N pack : The task of finding center solutions uses high simulation time, because it is comprised of many steps such as ranking each variable for the whole population and determining center variables from result of rank. The advantage is highly significant for large scale problem such as the ORPD problem due to high number of control variables and high population size.

2.
Many new solutions with high quality are found: The superiority can reduce the number of iterations, reducing the number of iterations, but the best optimal solution is still much more effective than that found by origin COA method.
The performance of the first modification is investigated by comparing the original COA and COA with the first modification, called COA1.

The Second Modification
In the second modification, we concentrate on the performance of the second new solution generation, which is accomplished by using Equation (19) above. The review on Equation (19) can indicate that the use of x new j,r1,p and x new j,r2,p has the same function, because the former is taken from a random solution and the latter is also taken from anther random one among the population. Thus, the conditions of (µ 4 < Pro 1 ) and (µ 4 < Pro 1 + Pro 2 ) are meaningless, but they result in more computation steps and time-consuming aspect. In addition, Equation (19) uses random variables within lower and upper bounds in the case that the two conditions above are denied completely. Some variables in the new solution are taken from two random solutions, while rest of variables in the new solution are randomly produced. Clearly, the combination is not from certain and believable arguments. Thus, in the second modification, we propose using the following model for finding one new solution for each pack.
By using the modification, COA can take advantages as follows: 1.
Reduce the calculation of Pro 1 and Pro 2 .

3.
Reduce random generation of variable x j,r in Equation (19).

4.
Shorten simulation time due to the reduction of three tasks above.

5.
N pack solutions with high quality are produced: The advantages can enable to reduce the number of iterations, leading to shorter simulation time, but the best optimal solution is still much more effective than that of original COA method.
The whole computation process of the proposed ICOA approach is given in Figure 3.

Population Initialization
In order to carry out the optimal solution search program of the proposed ICOA method, the first step is to produce a set of initial solutions, called population. On the contrary to other metaheuristic algorithms, ICOA as well as COA have N pack packs, and there are N coyote particles in each. Thus, the initial population generation is executed for each pack as the following model.
In the equation above, Sol min and Sol max are the minimum and maximum values of control variables of the ORPD problem and can be represented as below: Finding solutions of the ORPD problem is also supported by using Matpower program for obtaining remaining variables, which are different from control variables available in solution Sol d,p . Each solution Sol d,p is transferred to Matpower program for running, and then a set of dependent variables is obtained, consisting of apparent flow through each transmission line S branch,ij , reactive power output of generators Q Gen,has , and voltage of load buses V load,i . The three variables are imposed on inequality constraints (4), (8), and (9), and the three constraints are of course included in fitness function. However, other inequality constraints (5)-(7) together with two equality constraints (10) and (11) are not taken into account in fitness function. In fact, V Gen,has and TT i accompanied with Q cap,has are control variables, which are considered in each social condition of coyotes. Thus, such variables' boundaries are always monitored and corrected within their limits, as shown in Section 4.3 below. In another side, equality constraints (10) and (11) have been taken into account in the Newton-Raphson method, which is run in Matpower program. Consequently, before calculating fitness function, such three dependent variables are checked and penalized as follows: As a result, fitness function for evaluating quality of solutions can be constructed by the following equation: After having fitness function, the solution with the lowest fitness function in each pack is set to local optimal solution Sol best,p , and the best solution among all packs is set to the global best solution Sol Gbest .

Processes of Newly Updated Solutions
The proposed ICOA method updates new solutions two times. In the first time, N coyote solutions in N pack packs are newly produced by using Equation (23), but only one solution in each pack is updated for the second time by using Equation (24). As mentioned in Section 3, Sol new d,p and Sol new p are, respectively, new solutions in the first and the second updated solution processes.

Correction of Violation of New Solutions
After each updated solution process, correction is executed for handling violation of all control variables. The conventional method that has been applied so far is called boundary constraint handling method. By using the method, control variables are assigned to lower bound/upper bound if they are less/higher than lower/upper bound. Namely, the method can be implemented by the following models:

Termination of Iterative Search Algorithm
At the end of each iteration, one solution with the lowest fitness function is identified and set to the so-far best solution. In case that the current iteration is the last iteration, the so-far best solution is the best optimal solution for the run, and it is stored as a candidate. After finishing a predetermined number of runs, we have a set of candidates, and the best candidate is found by determining the lowest fitness function. On the other hand, the performance of the proposed method is also judged via minimum, average, and maximum fitness function values together with the standard deviation of the number of runs. Consequently, the termination condition of each run is dependent on the maximum number of iterations while the termination condition of each study case is dependent on the maximum number of runs. The maximum number of iterations N Iter is determined based on the number of control variables, which is dependent on the sizing of considered systems. The maximum number of runs N run is selected depending on both the performance of implemented methods and the sizing of considered systems. Normally, finding highquality solutions for small-scale systems is much easier and faster than for larger-scale systems. In dealing with three standard IEEE systems with 30, 57, and 118 buses, the proposed ICOA can find high-quality solutions by using N run = 50. Therefore, the selection is adopted for the purpose of reducing simulation time.

The Entire Computation Procedure
Finding ORPD problem solutions can be successfully accomplished by applying the flowchart in Figure 4.

Numerical Results
The efficiency of the proposed ICOA method has been simulated by using three standard tests with the IEEE 30-bus, IEEE 57-bus, and IEEE 118-bus systems considering three single-objective functions consisting of power-loss minimization, total-voltage deviation, and L index. The optimal solution search program of the proposed method has been coded on MATLAB program language and executed on a personal computer with a 2.4 GHz processor and a 4 GB RAM. In addition to comparisons with other methods available in previously published studies, we have also employed other state-of-the-art and popular methods such as salp swarm algorithm (SSA) [60], sunflower optimization algorithm (SFOA) [61], water cycle algorithm (WCA) [62], standard coyote optimization algorithm (COA) [56], improved coyote optimization algorithm with the first modification (COA1), and improved coyote optimization algorithm with the second modification (COA2). For finding optimal solutions, each method has been run 50 times in independent trials for each objective of each system. The setting of the adjustment parameters of the proposed method and other ones as well as test cases are summarized in Table 1.
It is noted that for each study system corresponding to each objective, the proposed ICOA has been run two different cases as mentioned in Section 2.
Case 1: All control variables are set to continuous variables (CV). Case 2: Reactive power output of shunt capacitors and tap changer of transformers are set to discrete variables (DV) in which a change step of reactive power output is 0.1 MVAR and that of tap changer is 0.01 pu.
The two cases are indicated in each comparison table because it is very important in evaluating the performance of methods. Basically, methods using continuous variables tend to obtain smaller fitness function than using discrete variables. For comparison with other implemented methods, we run ICOA and these methods with continuous variables only. V Gen,has (i = 1, . . . , 6) Q cap,has (i = 1, . . . , 9) TT i (i = 1, . . . , 4) The IEEE 57-bus system [63] 7 generator buses, 50 load buses, 80 branches, 3 VAR compensators, 15 transformers.

Performance Comparison Criteria for Different Methods
In the evaluation of the performance of optimization tools, while solving ORPD problem solutions, the sole key factor that has been used to point out the performance of compared methods is the best optimal solution reflected via fitness function [23,24]. Some of the previous studies [19,61] have used more comparison factors such as execution time in seconds, but computers, which have been used for simulation, have not been mentioned. It is obvious that computers with a higher processor spend less time for the same condition of control parameters or Matlab code. There is also another view point about the comparison of performance depending on the setting of adjustment parameters [64]. The higher population size and number of iterations can produce the higher number of new solutions, but it takes a long time to find optimal solutions. However, the statement is not completely accurate because there are some methods with complicated and time-consuming computation steps, although these methods and other ones produce the same number of new solutions for per iteration and per run. The simplest example for demonstrating the point is the comparison between PSO and DE. The two methods produce a set of new solutions equaling population size per iteration, and the same number of new solutions per run, but PSO has a smaller number of computation steps and simpler computation steps than those of DE. PSO updates velocity and position for all particles at the same time, leading to short simulation time, but DE must perform mutation for one solution by one solution, and then the crossover operation must be conducted. Thus, DE tends to spend longer simulation time.
Consequently, we suggest applying acceptably clear and fair comparison criteria for the evaluation of the performance of the proposed method and other ones. We focus on the following main issues: The best optimal solution: The best optimal solution comparison can point out the robustness of the method. The method with the better solution is considered to be stronger in the search process.
The solution search speed: It is also a very important factor for providing an accurate assessment of methods performance. A method with better solutions but with much longer simulation time is not a method with more potential. A method that finds better solutions and uses shorter simulation time is considered to be a better method, if the time comparison is carried out in the same conditions such as the same program language, the same program way, and the same computer. However, the expected conditions are hardly ever reached. For dealing with the issue, studies in [64,65] have proposed different comparison criteria for evaluating the optimal solution search speed of compared methods. Authors [64] have suggested a persuasive argument in regard to N pop and N Iter that the two factors directly influence the performance of implemented methods. The setting with higher values for the two factors can improve optimal solutions considerably, but convergence process to these solutions is time consuming and vice versa. Thus, they have suggested that the total number of newly produced solutions N spr or the number of quality solution evaluations over one run N sqe should be employed as the main comparison factor. Methods with the same or smaller number of quality solution evaluations but finding better optimal solutions become more powerful methods. Due to the suggestion, we construct different N sqe for different compared methods as follows: The three abovementioned formulas are applied for different methods with a different number of generations and different number of new solutions per generation. Formula (35) can calculate the number of solution quality evaluations for methods producing the number of new solutions equaling population size such as PSO, DE, SSA, GWO, BBO, etc. Formula (36) provides the number of solution quality evaluations for methods having two generations per iteration, and there is a population size of new solutions per generation such as ABC, TLBO, QOTLBO, CSA, etc. The model of Formula (37), which is applied for COA and the proposed method, is different than that of (35) and (36), because COA and ICOA have different search procedures than other methods. (N coyote .N pack ) is the number of quality solution evaluations for the first generation, while N pack is the number of solution quality evaluations for the second generation. On the other hand, (N coyote .N pack ) is equal to population size, and Formula (37) can be rewritten as follows: With the form, if we set (N coyote .N pack ) equally to the population size of other methods, COA and the proposed ICOA method produce a higher number of new solutions and spend more time than other methods with Formula (35) but the two methods produce a smaller number of new solution and use a shorter time for evaluating solution than other methods with Formula (36).
Although unlikely, the authors in [65] have suggested using computation time for fair comparison of performance. They did not consider N sqe , but they did consider the processor of computers as a strict factor. They have defined a time index of the compared method has, Ti i by using the following model: In Equation (39), GH i and GH (in GHz) are the processor of computers, which have been used to run the compared method has and the proposed method. Simtime i and Simtime (in second) are the simulated time of the compared method has and the proposed method. According to the viewpoint of the authors, the compared method has is slower than the proposed method if Ti i is higher than 1 and vice versa.
The stable search ability: The stable search ability is not the most important factor for evaluating the performance of method, but it can be used for further investigation in the case that two compared methods have the same optimal solution and the same search speed. Normally, the stability of the solution search is reflected via the average fitness of all runs. On the other hand, a better comparison can be made if the best fitness of all runs is plotted in the same figure. For comparison with other methods available in previous studies, we compare the average, but for implemented methods in the study, we plot all runs in the same figure, and a clear view can point out which method is superior to others. The applied comparison criteria cannot provide a fair comparison with the accuracy of 100%, but a general comparison with acceptable accuracy can be useful for concluding the performance of the proposed method as compared to other ones.

Testing Performance of COA and ICOA on Different Adjustment Parameters
In the section, we have presented the result obtained by COA and the proposed ICOA method by setting different values to the number of coyotes in each pack N coyote and the number of packs N pack . The task aims to find the most appropriate values for the two adjustment control parameters, N coyote and N pack , so that the proposed ICOA can get the highest performance with the best objective and the fastest search ability. In addition, a comparison of performance between COA and the proposed method as setting the same parameter was also investigated.
The considered ORPD problem is a complicated problem, and it takes optimization tools very much time to find optimal solutions due to the function of Matpower program. For judging the quality of one solution, the program must be run for dependent variables, and then fitness function can be calculated. Thus, the number of new solutions produced in each iteration N spi and the number of new solutions produced in each run N spr are criteria for comparing the search speed of different applied methods. We have set different values (2,3,4,5) to N pack and the values of N coyote have been 9, 6, 4, and 3, accordingly. Although the different setting results in different population size, N pop = 18, 18, 16, 15, N spi = 20 is the same for all settings excluding the case of N pack = 3 and N coyote = 6 with N spi = 21. Over 50 runs, the best, average, and worst fitness values together with the standard deviation and average execution time are reported in Tables 2-4 for TPL, TVD, and L index objectives, respectively. The best fitness values from the three tables can reveal the small impact of the two control parameters on the performance of the proposed ICOA, once the difference between the smallest best value and the highest best value is very tiny. For instance, the smallest power loss and the highest power loss are 4.5128 and 4.5238 MW. Similarly, the two values for TVD objective and L index objective are 0.0888 and 0.9833 pu, and 0.1242 and 0.1250 pu, respectively. However, there is a coincidence among the three cases since the smallest best fitness are found by the setting of N coyote = 4 and N pack = 4, and the highest best fitness are found by the selection of N coyote = 9 and N pack = 2. Clearly, a small number of packs could not produce promising quality solutions. However, it does not mean that a high number of packs is a good selection for a better comparison if the number of coyotes is not high enough. The setting of N coyote = 3 and N pack = 5 could not achieve the best result among all settings. In addition, the survey of standard deviation can show the best stability of search by setting N coyote = 9 and N pack = 2, while the second-best search stability is recorded by setting N coyote = 4 and N pack = 4. The best performance and the second-best stability can lead to a conclusion that the setting of N coyote = 4 and N pack = 4 is the most appropriate parameters for the IEEE 30-bus, and they are also applied for the IEEE 57-bus and 118-bus systems.  In regard to comparisons with COA method, the results of the proposed method are more favorable since almost all minimum values of the proposed methods are less than those of COA method, and the success rate (SR) of the proposed method is always higher than SR of COA. In fact, with different setting cases for N coyote and N pack , COA found different minimum TPL values to be 4.789, 4.682, 4.5777, and 4.595 MW, while those from ICOA are, respectively, 4.5178, 4.5154, 4.5128, and 4.5238 MW. The numbers indicate that the proposed method could reach less TPL by 0.2712, 0.1666, 0.0649, and 0.0712 MW, which are corresponding to the improvement level of 5.66%, 3.558%, 1.418%, and 1.55%. The improvement level of ICOA over COA is, respectively, 25.11%, 22.29%, 27.86%, and 21.69% for TVD case and is, respectively, 5.02%, 4.67%, and 0.72% for L index case. Clearly, the proposed ICOA approach could find a better optimal solution than COA for the same settings of control parameters. In relation to capability of handling constraints, the proposed method reported SR with 100% for TPL case and TVD case, and SR with 98% for L index case while COA successfully handled the system with SR less than 100%. Namely, SR is 92% for TPL case, 90% for TVD case, and 76% for L index case.

Investigation of the Proposed Method's Performance on the IEEE 30-Bus System
In the section, the results obtained by the proposed ICOA for the standard IEEE 30-bus system are compared to those from other implemented methods consisting of WCA, SSA, SFOA, COA, COA1, and COA2, and other methods available in previous studies. As mentioned in Section 5.1 above, the best optimal solution with the lowest fitness function is the first comparison factor for concluding the effectiveness of methods. Secondly, N sqe and Ti are the second factor for concluding the search speed of the proposed method. Finally, the average fitness of all runs is compared to other existing methods, and 50 runs of all implemented methods are plotted in the same figure for comparison of search stability. All the values regarding the three factors are reported in Tables 5-7 for TPL, TVD, and L index objectives. The minimum TPL in Table 5 shows that the proposed ICOA method can find the best solution for the case of continuous variables and discrete variables. The power loss of ICOA is the lowest value with 4.5128 MW for the case of continuous variables and 4.5138 MW for the case of discrete variables, while the best power loss in other studies is much higher.  For the case of continuous variables, the study [23] employed three methods in which the best one, HPSO-TS, found the most promising solution with 4.5213 MW. The study [19] implemented eight improved PSO approaches, and the best one determined the best solution with 4.52560 MW. The study [50] employed three methods in which the best one reported 4.5943 MW for power loss. Clearly, the proposed method found the best solution with 4.5128 MW, while the second-best method, GSA [29], reported 4.51431 MW, and the worst method, TS [23], reported 4.9203 MW. The further comparison indicates that the proposed ICOA method can improve the result from 0.01% to 8.3%. For the case of discrete variables, the study [37] employed three methods in which the best HAS method reported 4.9059 MW. The best ALO method among four applied methods in [44] reported 4.59 MW, while that of JA [47] was 4.5495 MW. Clearly, these methods found less effective solutions than the proposed method, and the improvement of result from the proposed method over other ones was from 0.79% to 8.6%.
The comparison with six other implemented methods also indicates the outstanding performance of the proposed method. The best power loss of WCA, SSA, SFOA, COA, COA1, and COA2 are 4.58218, 4.7175, 4.7364, 4.5777, 4.5340, and 4.5328 MW, which are higher than that of the proposed method by 0.068, 0.204, 0.222, 0.064, 0.020, and 0.019 MW corresponding to the improvement of 1.488%, 4.314%, 4.696%, 1.392%, 0.441%, and 0.415%. Clearly, COA1 and COA2 can improve the result relatively well as compared with COA, while the improvement of the proposed method over COA is highly significant. Thus, it can lead to a conclusion that the proposed ICOA method is really highly effective for minimizing TPL objective of the IEEE 30-bus system.
For comparison of the TVD objective shown in Table 6, all compared methods have considered only continuous values for all control parameters. Therefore, the comparison is only carried out for continuous variables. However, we have also reported obtained results from the proposed method for the case of discrete variables. As seen from Table 6, the best TVD of the proposed method, which is 0.0888 pu, is better than that of six other implemented methods and almost all methods in other studies excluding GSA [29] and QOTLBO [33], which reported 0.0676 and 0.0856 pu. Because taking control variables included generator voltages, tap setting of transformers, and power output of shunt capacitor banks of GSA and QOTLBO as the input parameters to run the Matpower program, the recalculated value of TVD is equal to 0.18622 and 0.10313 pu. Clearly, the two recalculated values are much higher than 0.0888 pu from the proposed method. The comparisons with other remaining methods imply that the proposed method can improve the best optimal solution quality from 0.45% to 57%, and the improvement is 27.2% as compared to COA method. Consequently, it can lead to a conclusion that the proposed method is very robust for finding solutions of the IEEE 30-bus system with TVD objective.
According to Table 7 showing result comparison of L index objective, GSA in [29] for the case of continuous variables and four methods in [44] for the case of discrete variables reported the best results with the lowest L index, approximately 0.11 pu, while that of the proposed method is 0.1242 and 0.12437 pu, corresponding to continuous and discrete variables, respectively. However, the recalculated value of L index is 0.1247 and 0.1241 pu for GSA and ALO, which are higher and approximately equal to that of the proposed method. The verification of reported solutions from PSO-IPG [19] and QOTLBO [33] confirms that the two methods have reported accurate solutions; however, the two methods have the same solution quality with the proposed ICOA method. As compared to other remaining methods, the L index values indicate the proposed method can find the optimal solution with the improvement up to 17.1% and that value is 0.8% as compared to the COA method. Thus, the proposed method can find the solution with the same quality as other leading methods for the case of L index optimization for the IEEE 30-bus system. For the comparison of solution search speed, the values of N sqe and Ti are also summarized in Tables 5-7. Table 5 indicates that the proposed method is ranked the fastest method together with six other implemented methods by spending 2000 solution quality evaluations, while the values from other methods in other studies are 4000, 8000, 10,000, 15,000, 18,000, 20,000, and 75,000. For comparison of converted simulation time Ti, the proposed is faster than other methods approximately from 1.3 times to 15 times. Similarly, Tables 6 and 7 also display the strong search speed of the proposed method once other methods have used N sqe = 4000, 10,000, 20,000, and 75,000 for TVD objective and N sqe = 4000, 8000, 10,000, 20,000, and 25,000 for L index objective. Accordingly, the proposed method is faster than these compared methods approximately from 1.6 to 3.5 times for TVD objective and approximately from 2 to 11 times for L index objective. All in all, the proposed method is more powerful than other compared methods available in previous studies in terms of solution search speed.
Compared to other implemented methods, the proposed method and these methods have approximately the same search speed due to the same setting of adjustment parameters; however, the proposed method has converged to optimal solutions with high quality, while others have not yielded the achievement. As shown in Section 5.2, the proposed method has reached a success rate of up to 100% for TPL and TVD cases and 98% for L index case, while the implemented methods have not reached the same success rate. For TPL case, only WCA has achieved SR with 100%, while SSA and SFOA have reached SR with 89% and 87%. SR of the three methods was, respectively, 100%, 92%, and 93% for TVD case and was, respectively, 90%, 78%, and 81% for L index case. Clearly, the proposed method is more capable of handling all constraints of the IEEE 30-bus system. One again, the conclusion can be confirmed via the best fitness of 50 runs plotted in Figures 5-7 for TPL, TVD, and L index objective, respectively. In Figure 5, all applied methods have high fluctuations, but only the proposed ICOA reaches many better solutions with less power loss than other methods. In Figures 6 and 7, the proposed ICOA has smaller fluctuation than other methods, and it also find solutions with much lower values of TVD and L index than others. Clearly, the figures show that the proposed method is more stable than other ones for finding solutions, since almost all solutions from the proposed method have lower fitness than other methods. Figure 8 shows the voltage of all loads for TVD objective. It is noted in the figure that voltage of loads closest to 1.0 pu is the best, with the smallest deviation, because 1.0 is an expected value. Voltage from the first load to the 24 loads indicates that the proposed method always keeps a small deviation with the expected voltage 1.0 pu. Thus, the proposed ICOA has the best value of TVD, as shown in Table 6 and Figure 6.
Optimal solutions obtained by the proposed method for IEEE 30-bus system are given in Table A1 in Appendix A.

Investigation of the Proposed Method's Performance on IEEE 57-Bus System
In this section, three different objectives of the IEEE 57-bus are considered to be optimized. Tables 8-10 show the information comparisons of all compared approaches with the proposed ICOA approach for TPL, TVD, and L index optimization. For the three objective comparisons, all the compared methods in other studies have been run only by setting control variables to continuous values. Thus, the comparisons with these methods and four other implemented methods are only performed for the case of continuous variables. However, all results obtained by the proposed methods for the case of discrete variables are also reported.
For the comparison with four implemented WCA, SSA, SFOA, and COA approaches, the high performance of the proposed ICOA approach can be identified clearly and exactly, since approximately all results from the proposed ICOA approach are superior to those from these methods. The proposed approach can reduce TPL to 22.376 MW, while the four methods found the best solutions with 26.0402, 25.3854, 26.6541, and 24.5358 MW. Clearly, the proposed approach can obtain less TPL than these methods by 3.6642, 3.0094, 4.2781, and 2.1598 MW. Similarly, the proposed ICOA approach continues to be more effective than these methods for TVD optimization and L index optimization, since the best TVD and the best L index are, respectively, less than these methods by 0.126, 0.335, 0.186, and 0.066 pu,  and 0.02721, 0.03831, 0.03141, and 0.02401 pu. The values can result in the improvement level of the proposed ICOA approach over these methods. Namely, the improvement for TVD is 20.790%, 55.346%, 30.772%, and 10.907% and for L index is 10.81%, 15.22%, 12.47%, and 9.5%, respectively. The analysis of the best TPL, TVD, and L index can confirm the high efficiency of the proposed approach as compared to the four implemented methods. For judging the stabilization of search ability of the proposed approach, the average for 50 runs of the three optimization cases is compared among all methods. 24.4922 ----X -PSO [45] 24.3826 ----X -CSA [45] 24.2619 ----X -WCA [46] 24.82 ----X -GBBWCA [46] 23.27 ----X -MFO [48] 24.25293 --9000 -X -MOGWA [49] 21.171 ------GSA [29] 23 It can be realized that the search process of the proposed approach is more stable than the approaches.  [43] 0.6589 --100,000 -X -WCA [46] 0.6631 --50,000 -X -GBBWCA [46] 0.6501 --50,000 -X - The performance of the proposed ICOA approach is also evaluated effectively and strongly as comparisons with other approaches are carried out and analyzed. In fact, the proposed ICOA approach can find out more effective optimal solutions than almost all methods for TPL optimization and all methods for TVD optimization. Only MOGWA [49] has shown TPL of 21.171 MW, which is less than that of the proposed ICOA approach; however, the verification implies that the optimal solution has not been shown in the study [49]. Other comparisons with the second-best solution and the worst solution indicate that the proposed ICOA approach can reduce TPL from 0.894 to 3.264 MW, corresponding to the improvement level from 3.99% to 14.59%. The proposed ICOA approach also finds less TVD than others from 0.0433 to 0.1258 pu, corresponding to the improvement level from 7.16% to 20.79%. The average of 50 runs found by the proposed ICOA approach is less than almost all methods excluding ALC-PSO [17]. However, it should be pointed out that the method has used a much higher value for N sqe than the proposed ICOA approach. Namely, N sqe is 30,000 for ALC-PSO [17] but 6000 for the proposed ICOA approach. The reported values of Ti also lead to the same manner, since the proposed method can be faster than other methods by approximately 15-62 times for TPL optimization and approximately 21 times for TVD optimization. Consequently, there are enough evidences consisting of the best optimal solution, the best stabilization of search, and the fastest search ability to conclude that the proposed ICOA approach is the most effective solution for the ORPD problem with IEEE 57-bus system.
Optimal solutions obtained by the proposed method for IEEE 57-bus system are given in Table A2 in Appendix A.

Investigation of the Proposed Method's Performance on the IEEE 118-Bus System
In this section, the proposed ICOA is applied for the purpose of reducing TPL, improving voltage profile, and enhancing voltage stabilization of the IEEE 118-bus system. In addition, the application of SSA, WCA, SFOA, COA, COA1, and COA2 methods for the system is also accomplished for further investigation of performance. The main comparison data consisting of minimum fitness, average fitness, N sqe , and Ti together with others are reported in Table 11 for TPL objective, in Table 12 for TVD objective, and in Table 13 for L index objective. Other compared methods have set control variables to continuous values only and ignored the case of discrete values. Thus, analysis as well as discussion on results is only performed for the case of continuous variables. However, results from the proposed method for the case of discrete values are also reported in each comparison table.  [33] 112.2789 113.7693 -10,000 -X -TLBO [33] 116.4003 121.3902 -10,000 -X -Comparisons of minimum and average fitness values among the proposed and implemented SSA, WCA, SFOA, COA, COA1, and COA2 methods shows that the proposed method is more effective in finding the best optimal solutions and more stable in the whole search process because the two values of the proposed method are less than those of these implemented ones for three objectives. The same values of N sqe and Ti confirm that all methods have been run on the same conditions and fair comparisons have been executed. Nevertheless, the proposed method is always superior to the methods, as observed in Figures 9-11. Approximately all runs of the proposed method could find better quality solutions for different objectives, and the fluctuations of the search process are very small, but very high fluctuations can be seen from the curves of other ones. Moreover, Figure 12 shows that nearly all loads of ICOA have voltage approximate with the reference value, 1 pu, while those from other ones are very different than 1 pu. The reported SR values in the three tables are useful for demonstrating the high potential of dealing with complicated constraints of the proposed method. For the three study cases, the proposed method handled all constraints as successfully as/more successfully than WCA, but its capability of dealing with constraints was better than other ones, especially much better than COA.    The comparisons with other existing methods available in other studies can point out that the proposed method has a worse result than PSO-IPG [19], SARCGA [10], and QOTLBO [33] for TPL objective, and three PSO methods in [19] consisting of PSO-SWT, PSO-PGSWT, and PSO-IPG for L index objective. Nevertheless, it needs to confirm that not all the methods are really superior to the proposed method for the two objective functions. The verification of reported solutions can lead to the conclusion that SARCGA is not a real effective method, because there was no solution reported in [10]. QOTLBO [33] is only more effective than the proposed method for TPL objective but it is less effective for two cases of TVD and L index objectives. However, it should be noted that QOTLBO used a higher number of solution quality evaluations, 10,000, while that of the proposed method was only 8000. Thus, it can be said that the proposed method is more effective than QOTLBO. Among the rest of methods in [19], PSO-IPG is more effective than the proposed method for two cases of TPL and L index objective, but it is less effective than the proposed method for TVD objective. PSO-SWT and PSO-PGSWT used the same number of solution quality evaluations as the proposed method, but the two methods found less effective solutions for two cases of TPL and TVD objectives. However, it should be noted that these methods in [19] were run by setting upper bound of the load bus voltage to 1.1 pu, while the proposed method and other ones were run on the bound of 1.05 pu. Thus, we have continued to implement the proposed method by assigning the upper bound of load bus voltage to 1.1 pu. All results obtained for three single objectives are shown in Table 14. It can be seen that the best TPL, TVD, and L index from the method for the same comparison condition with methods in [19] (i.e., continuous variables, the same upper bound of load bus voltage of 1.1 pu) are 111.2344 MW, 0.1561 pu, and 0.5678 pu; meanwhile, these values of the best PSO-IPG method in [19] are 115.06 MW, 0.162 pu, and 0.0568 pu. Clearly, the proposed ICOA method can improve performance of up to 3.32% for TPL objective and 3.64% for TVD objective over PSO-IPG, while the proposed method has the same optimal solution quality with PSO-IPG for L index case. For comparison with other methods shown in Tables 11-13, the improvement level of the proposed method can be from 0.37% to 12.93% for TPL objective, from 0.93% to 92.82% for TVD, and from 0.15% to 56.33% for L index.  [33] 0.0608 0.0631 -10,000 -X -TLBO [33] 0.0613 0.0626 -10,000 -X - All methods in [19] were run by setting N pop = 40 and N Iter = 200 corresponding to N sqe = 8000, which is similar to that of the proposed method but their execution time was reported completely differently. For instance, the simulation time for TPL objective was approximately about 95 s. However, for TVD objectives, some methods reported about 75 s, and some other methods reported about 50 s, whilst one method, PSO-PGSWT, reported over 110 s. For L index, simulation time was over 110 s for five methods, but it was under 60 s for three other ones. The unreasonable issue is that methods using pseudogradient search reported the shortest time, even though pseudogradient search is complicated and it takes more time to determine the better direction for velocity. The comparisons of search speed easily led to the conclusion that the search speed of the proposed method is ranked in the leading method group, with once only methods in [19] having been executed for solving the IEEE 118-bus system by setting the same number of solution quality evaluations, which is less than that of all other methods. By considering Ti, it can be indicated that Ti of other methods is higher than 1, i.e., they were slower to converge to their optimal solutions. Consequently, it can conclude that the proposed method is effective for solving the IEEE 118-bus system.
Optimal solutions obtained by the proposed method for IEEE 118-bus system are given in Tables A3-A8 in Appendix A.

Conclusions
This paper has presented an ICOA with intent to enhance convergence speed to the global optimum solution with the lowest value of fitness function as solving different systems of reactive power dispatch problem. The proposed method has been developed successfully in dealing with all constraints of the problem and in finding promising solutions with higher quality than those of the standard coyote optimization algorithm. In terms of the optimal solution search function for a general optimization problem, the implementation of the proposed method is simpler than that of COA, because the determination of center solutions in each pack is no longer essential in the proposed method and the offspring generation process is suggested to be changed in the proposed method. In COA, finding center solutions is time consuming, and the quality of center solutions is low, while offspring determination also needs high simulation time due to the process of finding control variables one by one. About the result comparison for the three IEEE power systems for the ORPD problem, each of two modifications was demonstrated to be more effective than the original procedure of COA via the comparisons of COA1 (COA with the first modification) and COA2 (COA with the second modification) with the standard COA. All the best optimal solutions from the two modified versions have better quality than those of COA for all objectives of three power systems. Furthermore, the success rate of the proposed method is always higher than that of COA. As compared to other methods consisting of implemented methods such as WCA, SSA, and SFOA, and other existing methods available in previous studies, the performance of the proposed method was thoroughly evaluated by focusing on three comparison criteria: the best optimal solution quality, the search speed, and the stabilization. The proposed method found more promising results than all methods for the IEEE 30-bus system in terms of better quality of the best optimal solution, faster search speed with the lower number of solution quality evaluations, and good stabilization. For the IEEE 57-bus and IEEE 118-bus systems, the proposed ICOA approach was more effective than all other methods. Consequently, it indicated that the proposed method is a powerful optimization tool for dealing with the ORPD problem, and it can be one of the most effectively alternative methods for other optimization problems in power systems. Data Availability Statement: Data of the study are taken from [63]. The best solution and the worst solution in the pack p Sol best,r1 , Sol best,r2 , Sol best,r3 , Sol best,r4

Conflicts of
The best solutions in all packs picked up randomly Shunt admittance of branch ij µ 1 , µ 2 , µ 3 , µ 4 , µ 5 , µ 6 , µ 7 Random numbers in range from 0 to 1 α i , α j Phase angles of voltage at bus i and bus j φ 1 , φ 2 , φ 3 Penalty factors in fitness function Appendix A