Deterministic‐like solution to the non‐convex economic dispatch problem

Funding information Advanced Power and Energy Centre (APEC) at Khalifa University of Science and Technology, Abu Dhabi, United Arab Emirates, Grant/Award Number: RCII-006-2018 Abstract This paper proposes a novel stochastic solution method for the static nonconvex power economic dispatch problem. All practical features such as valve point effects, ramp rate limits, prohibited operating zones, multiple fuel options, spinning reserve constraints, and transmission losses are considered. To develop the proposed algorithm, two modifications that significantly enhance the exploitation capability of the artificial bee colony algorithm are introduced. The first modification concerns onlooker and employed bees to focus their search around the best solution found so far. The second modification optimizes the escaping behavior from local minimums using scout bees. The proposed modifications form a multi-level concentric search around the best solution found. Finally, combining the modified artificial bee colony algorithm with Levy flight cycles, which improve the escaping capability from local minimums, leads to the proposed algorithm. The algorithm parameters have been tuned to provide a deterministic-like solution with ten benchmark problems. The improvement added by each proposed modification to the artificial bee colony algorithm has been confirmed using the Wilcoxon rank-sum test. The obtained results by the proposed algorithm are superior compared to those reported in the literature. Moreover, estimated probabilities of more than 99.9% to obtain the global optimal solution are achieved.


INTRODUCTION
The economic dispatch problem (EDP) is a vital operational problem in power systems. High non-convexities and nonlinearities levels due to considering practical aspects such as valve point effects, prohibited operating zones, ramp rate limits, multiple fuel options, and transmission losses render the problem one of the most challenging optimization problems in the power system literature. Due to the high penetration levels of renewable energy resources, it is expected that the problem will be solved within smaller dispatching intervals in order to update the operating set points while following the most updated forecasts. obtaining the global optimal solution of the non-convex EDP can lead to significant monetary losses, especially for large scale and highly non-convex systems [1,2]. Several techniques have been proposed to solve the non-convex EDP. These techniques can be aggregated into two major groups: deterministic techniques and stochastic ones. Considering the deterministic approaches, a fast˘-Iteration method is introduced to solve the EDP while considering the prohibited operating zones in [3]; however, this method cannot deal with generation units incorporating valve point effects and multiple fuel options. In [4], the non-convex EDP has been formulated as a rank-relaxed semi-definite programming problem, which is solved using the branch-and-bound technique combined with the method of convex iteration, yet this approach does not consider prohibited operating zones or transmission losses. In addition, it can only approach the global optimal solution with limited accuracy. A steepest decline method is proposed in [5] mainly for solving the EDP with valve point effects while considering transmission losses. Unfortunately, this algorithm may fail to find the global optimal solution for some instances of the problem if the algorithm is not properly tuned for the considered instance.
To summarize, there is no deterministic technique known up to date that can handle all the practical features of the non-convex EDP simultaneously while providing the global optimal solution with high accuracy. On the other hand, stochastic heuristic algorithms can be applied while considering all the practical features of the non-convex EDP. Furthermore, it can be noted that approximately all the known global optimal solutions of the non-convex EDP up to date have been discovered using stochastic heuristic algorithms. One of the powerful stochastic heuristic algorithms that have been proposed as a generic optimization solver is the artificial bee colony (ABC) algorithm. It has been shown in [6] that the ABC algorithm is capable of providing better statistical results than those of the particle swarm optimization (PSO) considering several benchmark test functions. On the other hand, the Levy flights which mimic the flight behaviour of many insects and animals have been proposed to improve the performance of several stochastic algorithms by enhancing their capability for escaping from local minimums [7][8][9]. To improve the performance for variants of the ABC algorithm, these variants have been augmented with the Levy flights in [10][11][12]. On the other hand, it is known that the ABC algorithm has good exploration and poor exploitation capabilities [13,14]. To improve the exploitation capability of the ABC algorithm, it is proposed in [13] modifying the update equations of the onlooker and employed bees so that these bees are guided by the best solution found so far (Gbest). The produced algorithm is known as Gbest-guided ABC (GABC). Due to the need for improving the exploitation capability, it is proposed in this manuscript to combine a new Gbest-guided version of the ABC algorithm with the Levy flights (GABC-LF) to produce better results than those obtained by combining the basic ABC algorithm with the Levy flights (ABC-LF). However, up to the authors' knowledge, neither the ABC-LF nor any other algorithm produced by combining a variant of the ABC algorithm with the Levy fights has been applied for solving the non-convex EDP before. Furthermore, it is observed that enhancing further the exploitation capability will improve the performance of the GABC-LF. The exploitation capability has been further enhanced by adding scaling factors to the update equations of the onlooker and employed bees. In addition, the scout bee has been modified to generate random solutions within a limited volume around the Gbest. The resultant algorithm is designated as Gbest-focused ABC-LF (GFABC-LF). This algorithm is powerful enough to provide a deterministiclike solution to several benchmark problems of the static non-convex EDP.
On the other hand, the operational mechanism of the stochastic heuristic algorithms cannot be proved mathematically. Consequently, these algorithms are compared based on a statistical analysis of their results delivered using several trials of the algorithm. The number of trials adopted in the literature to get the best possible solution differs from a publication to another. For example, 100 trials are used in [15][16][17][18][19], while 50 trials are utilized in [20,21]. Therefore, the best possible solution is obtained within a total computational time equal to the average computational time per a single trial multiplied by the total number of trials. The total computational time can be reduced by decreasing the number of trials. This can be accomplished if the proposed algorithm is powerful enough to achieve a high probability of obtaining the global optimal solution. Yet, it is observed in the literature that the published statistical results typically do not discuss the probability of obtaining the optimal solution. In addition, the effect of increasing the total number of objective function evaluations (TOFEs) on this probability is not investigated, although this investigation can reveal important information that can help in reducing the overall computational time. The main contribution of this paper is to present a new solution method (GFABC-LF) that is capable to deliver a high probability of achieving the global optimal solution while aggregating the following advantages: -Unlike deterministic methods, the proposed solution method provides the global optimal solution with high accuracy while considering all the practical features of the non-convex EDP. -Compared to previous stochastic algorithms, the proposed solution method is capable of achieving the global optimal solution with a high frequency at acceptable values of TOFEs and total computational time. Furthermore, the proposed method shows a high degree of robustness against any variability, which results from considering different problems with different specifications.
The capability of obtaining the global optimal solution with a high frequency is translated into an estimated probability of hitting the global optimal solution at a certain TOFE value, and it is shown that this probability can be controlled and exploited for a wide range of benchmark problems including the most challenging problems to achieve a deterministic-like behaviour within a reasonable computational time. An important benefit of introducing a deterministic-like behaviour is to reduce the total computational time during the real-time application. This time reduction is beneficial to confirm obtaining the global optimal solution.
The paper is organized as follows. The formulation of the non-convex EDP is presented in Section 2. Section 3 reviews the ABC algorithm and the Levy flights mechanism. The proposed solution method is introduced in Section 4, and the simulation results are presented in Section 5. Section 6 presents a sensitivity analysis and discusses the parameters tuning for the proposed algorithm. Section 7 accommodates a discussion of related practical aspects. The future work and the dynamic nonconvex EDP are discussed in Section 8. Finally, the conclusion is delivered in Section 9.

MATHEMATICAL FORMULATION OF THE PROBLEM
The standard formulation of the problem is adapted to allow comparing the proposed method with previously proposed solution techniques using well-known benchmark problems. This formulation is as follows [15,16]: subject to: where D is the total number of generating units, j is an index for these units, P j is the output power from unit j, C j expresses the cost of power generated by unit j, C T denotes the total operational cost, P Load stands for the total system load, P Loss signifies the total transmission losses, k is an index for prohibited operating zones, K j denotes the total number of prohibited operating zones for unit j, P j 0 stands for the power output of unit j at the pervious time interval, Ω is the set of units having prohibited operating zones, P l j,k represents the lower limit of the kth prohibited zone of unit j, P u j,k stands for the upper limit of the kth prohibited zone of unit j, DR j and U R j are the downramp and upper-ramp rate limits of unit j, respectively, P min j and P max j are the lower and upper generation limits of unit j, correspondingly, S T expresses the total required spinning reserve of the system, and S max j is the maximum spinning reserve of unit j.
If valve point effects and multiple fuel option are not considered, the generation cost is usually modelled using a convex quadratic cost function as follows: In (6), a j , b j , and c j are cost coefficients specific to unit j. Considering valve point effects leads to a non-convex cost function, which can be expressed using: where e j and f j are cost coefficients for modelling valve point effects. When multiple fuel options are considered in addition to valve point effects, the cost function is modelled as follows [22]: where a j,n , b j,n , c j,n , e j,n , and f j,n are cost coefficients of unit j using fuel type n. P min j,n and P max j,n stands for the lower and upper limits of unit j using fuel n. In (2), total transmission losses, P Loss , can be approximated with Kron's loss formula which is stated as follows: B m, j , B 0 j , and B 00 are loss coefficients. Equation (3) forms the upper and lower operational bounds for unit j while considering ramp rate limits. On the other hand, if ramp rate limits are not considered, Equation (3) is reduced to (10).
Finally, (4) represents additional bounds due to considering prohibited operating zones, and (5) expresses the spinning reserve constraint.

ABC ALGORITHM AND LEVY FLIGHTS
The proposed solution method is developed based on the wellknown ABC algorithm while adapting Levy flights. It is known that the basic ABC algorithm has good exploration capability and is poor at exploitation [13,14], so enhancing the exploitation capability of the ABC will create a stronger algorithm. On the other hand, the Levy flights have been combined recently with several algorithms in order to improve their performance by enhancing their escaping capability from local optimums [7,8]. Therefore, before presenting the proposed algorithm, the basic ABC algorithm and the Levy flights mechanism are reviewed in the following subsections.

ABC algorithm
The ABC algorithm has been proposed by Karaboga in [6]. In this algorithm, candidate solutions for the problem are represented as food source positions. Fitness function values indicate the amount of food (nectar) in corresponding food source positions. The algorithm has three moving patterns used to search for better new solutions. These patterns are referred to be done by three types of bees: employed, onlooker, and scout bees. Before applying the moving patterns, an initial population of food sources is generated using the following formula [23]: where P i, j denotes element j of the solution vector P i (i ϵ {1, 2, ….., N}, where N is the population size, and j ϵ {1, 2, ….., D}, D indicates the problem dimension). P min j and P max j are the lower and upper limits of dimension j. rn i, j is a random number uniformly distributed between 0 and 1. After the initialization, the three types of bees are repetitively applied in order to search for new better solutions. The number of employed bees is equal to the population size. Each employed bee works to update one solution within the population using the following expression: where i, k ϵ {1, 2, ….., N}, and i ≠ k, j ϵ {1, 2, ….., D} is randomly chosen for each i.̇rn i, j is a random number uniformly distributed between −1 and 1. The fitness values of the new solutions are calculated using: where f i is the fitness value corresponding to the solution i. C i T stands for the total cost of solution i. Since C i T is always positive for the EDP, Equation (13) is sufficient in a minimization problem. If a new solution has a fitness value higher than the old solution, then this new solution will replace the old one; otherwise, the old solution is retained. This selection process is denoted as the greedy selection [24]. After the application of the employed bees, the onlooker's phase starts by assigning a probability, p i , for each solution as follows: Thereafter, each onlooker bee chooses a solution based on the probability value of this solution. Once an onlooker bee chooses a solution, it generates a new one using (12). The fitness value is evaluated for the new solution, and a greedy selection process is applied as explained with the employed bees. There is a counter associated with each solution. This counter is reset to zero if the newly proposed solution by the employed or onlooker bee replaced the old solution. Otherwise, the counter is incremented by one. In each iteration, if the maximum counter value in the population exceeded a certain predefined value, known as the limit, the scout bee is applied to replace the associated solution with a randomly generated one using (11).

Levy flights mechanism
It is observed that the flight behaviour of many insects and animals follows certain characteristics of what is known as Levy flights [25]. The Levy flights perform random walk with step lengths drawn from the Levy probability distribution [26]. Usually, the Mantegna algorithm is adapted to generate steps from a symmetric Levy stable distribution. Using this algorithm, the step length is calculated using: where r x and r y are two normally distributed random variables, The standard deviation values are given as follows: where is a constant and ( * ) expresses the gamma distribution function. After computing the step length, the element j of solution i is generated using the following formula: In (18), P i, j and P new i, j are the jth components of old and new solutions, respectively, is a scaling factor, rn i, j is a normally distributed random variable with zero mean and standard deviation of 1 evaluated in correspondence to element j of solution i. P Gbest i, j is the jth dimension of the best solution obtained so far. The best solution found so far is also denoted as Gbest.

PROPOSED SOLUTION METHOD
Since the ABC algorithm has good exploration and poor exploitation capabilities, enhancing its exploitation will improve the performance while benefiting from the good exploration.
To improve exploitation, two enhancements to the basic ABC algorithm are recommended. The first enhancement concerns the update equation used by the onlooker and employed bees. The second improvement modifies the scout bee performance. After improving the exploitation, an optimized number of Levy flight cycles is combined with the modified ABC algorithm to further improve the overall performance of the algorithm.

Fine search equation (modified update equation)
In order to improve the exploitation capability of the ABC algorithm, the following equation is proposed for the onlooker and employed bees instead of (12): where rn i, j is a random number that obeys the uniform distribution between 0 and 1. r 1 , r 2 ϵ {1, 2, ….., N} are two integer numbers selected randomly while achieving r 1 ≠ r 2 , 1,2 denote two focusing factors, 1 and 2 , with the purpose of adjusting the focusing level for the search process around the Gbest. 1 is used with employed bees, while 2 for onlookers. These factors can have the same value. However, different values are suggested and used in order to generate solutions around the Gbest over two different levels. Small values of 1,2 perform a fine search around the Gbest. Since (19) is executed for each dimension separately with a small value for the focusing factor, a small change that causes a decrement in the cost function over a certain dimension will have a high chance of being caught by the algorithm. In the ABC algorithm [6] and the GABC algorithm [13], the scaling parameters 1,2 do not exist in the update equations. The non-existence of 1,2 is equivalent to having these parameters with a fixed value of 1 in the update equations of these algorithms. Equation (19) provides better statistical results for the non-convex EDP if relatively small values of 1,2 such as 0.01 or 0.1 are adopted in the proposed algorithm. Since small values of 1,2 introduce small changes to the Gbest and enhance the exploitation capability, Equation (19) is designated as a fine search equation.

Improved equation for the scout bee
In the basic ABC algorithm, the scout bee uses (11) to replace the exhausted solution with a new one. Therefore, exhausted solutions are replaced by entirely random ones that neglect the information gained by the algorithm as the search progresses. To make the scout bee performing better, the following equation is proposed to exploit the information gained by the algorithm so far through generating random solutions around Gbest as follows: wherern i, j is a random variable uniformly distributed over the interval [−0.5, 0.5], expresses a scaling factor which adds an upper bound for the distance between the newly generated solution and Gbest. A small value for the scaling factor will reduce the ability of the algorithm to escape from local minimums. While a large value results in the opposite, it will incur information loss. Therefore, a compromise value has to be selected to optimize the performance. The term Gbest-focused ABC (GFABC) is used to describe the ABC algorithm when Equations (19) and (20) are utilized.

Constraints handling mechanism
Constraints handling mechanisms can affect the statistical results of any heuristic algorithm. The adapted mechanism is a generic one and is used to handle the equality constraint while satisfying the upper and lower bounds of decision variables. In order to satisfy the generation limits, the unit that violates its generation limit is forced to work at that violated limit as follows: For the case of considering ramp rate limits, a new set of upper and lower bounds for the units are formed using Equation (3) before applying (21). The power balance constraint is satisfied through selecting a unit randomly to carry the mismatch between the total power generation and total load. If this unit cannot carry the mismatch solely, the unit is fixed at its generation limit that minimizes the mismatch, and then another different unit is picked randomly to carry the remaining mismatch. The process is repeated as shown in Algorithm 1 till the mismatch is equal to zero. In Algorithm 1, the power mismatch is calculated using the following formula: To incorporate transmission losses, few modifications are added to Algorithm 1. The first modification is to calculate the mismatch as follows: The second modification incorporates utilizing additional inner loops, which iterate between updating the output power of the selected generation units, calculating the transmission losses using (9), and calculating the mismatch using (23) until a certain tolerance is met. In the next algorithm, N x stands for the number of solutions passed to the algorithm. 3: Calculate mismatch i using (22). 4: If mismatch i > 0 5: Select a unit randomly, unit m.

6:
While Select a different unit to be m.

GFABC with Levy flights
To further improve the performance, Levy flight cycles are combined with the GFABC algorithm to produce the proposed algorithm, GFABC-LF. Before starting the algorithm, two vectors are created: the cost vector and the counter vector, each has a dimension of N × 1, where N is the population size, and each element of these vectors corresponds to one solution in the population. The cost vector stores the cost of each solution in the population, and the counter vector contains the number of successive iterations for which the corresponding solution has not been improved. If any solution is improved by the onlooker or employed bee, the corresponding value in the counter vector is reset to zero; otherwise, it is incremented by one. For each iteration, the maximum value in the counter vector is checked against a fixed value known as the limit. If this maximum value exceeded the limit, the corresponding solution is replaced with the scout bee. Algorithm 2 lists the steps of updating the cost vector using the greedy selection process, in addition to updating the counter vector. In Algorithm 2, cost i and counter i express element i in the cost and counter vectors, correspondingly, and N x signifies the number of solutions passed to the algorithm.  (13) and (14).
14: Apply the scout bee for the solution that has the maximum counter value and exceeding the limit, replace this solution using (20), execute Algorithm 1 for the new solution, and update relative cost i value.

19: End While
The final proposed technique is detailed in Algorithm 3, in which the total number of Levy flight cycles per each iteration is designated as N L , and an index t is used to point at each cycle of the Levy flights. The update of Gbest is done as follows: The pattern of jumps added to a starting point of (0, 0) using 50 Levy cycles

FIGURE 2
The pattern of jumps added to a starting point of (0, 0) using 500 Levy cycles It should be pointed out that there is more than one possible way to combine the Levy flight cycles with the GFABC algorithm. One of these ways is to add the Levy cycles before the scout bee phase while updating the counter vector within each iteration of the Levy cycles. This approach leads to relatively competitive results under the condition of using a much higher limit value compared to that used with Algorithm 3. However, after investigating several options for the merge process, the one described in Algorithm 3 is selected.

Illustration of the proposed modifications
This subsection elaborates on the operation and motivation of the proposed modifications. There are four moving patterns used by the GFABC-LF: three are made by the three types of bees, and the fourth one is due to the Levy flights. These moving patterns are centralized by the Gbest (P Gbest ), the best solution found so far. Once a change made by any moving pattern leads to a better cost value than that of the Gbest, the Gbest is updated with the corresponding better solution. The Levy flights is a strong mechanism to generate new random solutions. This mechanism has achieved a lot of success and attracted wide attention since it has been proposed for the first time [8]. Figure 1 presents a two-dimensional plot for the changes added to a starting point of (0, 0) using 50 Levy cycles, while Figure 2 utilizes 500 Levy cycles. It is clear from Figures 1 and 2 that  Figures 1 and 2, it can be observed that larger jumps appear in Figure 2 than those in Figure 1. This demonstrates that utilizing more Levy cycles increases the range of the possible changes introduced by Levy flights. This can help in increasing the diversity in the population and the chances of escaping from local minimums that require relatively large jumps. Therefore, more than one Levy cycle are utilized in each iteration of the algorithm. In the proposed algorithm, the best solution found so far (Gbest) is used instead of the starting point (0, 0). On the other hand, the modifications proposed for the ABC equations aims at improving the exploitation capability of the algorithm through performing a multi-level concentric search around the P Gbest j as clarified in Figure 3. This figure shows three zones for modifying the output power of a unit with valve point effects using the three types of bees. In Figure 3, the solid black line represents the cost function of a typical generation unit with valve point effects. The centre of the three circles (represented by X) indicates the cost of the P Gbest j , which is the jth output power of the best solution obtained so far. The radius of the smallest dotted blue circle reflects the limited random change added to the P Gbest j by the onlooker bee in Equation (19) and this change has a maximum value of 2 (P r1, j − P r2, j ), the radius of the solid green circle reflects the limited random change added to the P Gbest j by the employed bee in Equation (19), and this limited change is controlled by 1 . 1 and 2 are meant to take small values in order to improve the exploitation capability, and it is thought that selecting different small values for 1 and 2 achieves better exploitation over two levels and balances the exploitation accuracy and speed. In Figure 3, the radius of the largest dashed red circle images the controlled change introduced by the scout bee to the P Gbest j . This change has a maximum value of × (P max j − P min j ) as conferred by Equation (20). The value of is chosen large enough to allow for escaping from local minimums by the scout bees without jumping away from the current achieved search progress. For a certain unit, when, for example, a change is made towards increasing the generation of this unit by the GFABC-LF algorithm, this increase is equally balanced and substituted by a decrease in the output power from another randomly selected unit using the constraint handling mechanism. For any change made by the four moving patterns of the algorithm, if a solution is found with a total cost lower than the cost of P Gbest j , it will replace P Gbest j immediately, so as the search progresses, the centre of the three circles (P Gbest j ) moves on the black line in Figure 3 for all the units in an iterative way until the global optimal solution is obtained. Large jumps, as indicated in Figure 3, can occur in the proposed algorithm thanks to the Levy flights and the constraints handling mechanism. On the other hand, regarding the ability to escape from local minimums, the algorithm has a strong ability to escape from local minimums since in addition to the scout bee, the Levy flight cycles also provide a strong mean for escaping from local minimums, and an optimized number of five Levy flight cycles in each iteration are used instead of just one cycle as with the original Cuckoo search algorithm [27].

5.1
Case studies setup

Parameters setting
For all case studies, the parameters of the proposed algorithm have been fixed to the values in Table 1. These values have been obtained after several experiments to get a satisfactory performance over a large number of benchmark problems. The adopted methodology to obtain the values in Table 1 is discussed in Section 6. To compare the proposed algorithm with other algorithms, the TOFE value per each trial has been considered. A personal computer with core i7 (3.4 GHz) processor and 8 GB RAM is utilized for performing the case studies.

Accuracy measurement
When executing n trials of the algorithm, the accuracy of approaching the global optimal solution is measured using the following expression: where X L indicates the number of times, out of n trials, the global optimal solution is obtained within an accuracy of 1e −L , and L ∈ Z + determines the accuracy level, r is an index for trials, C r BT denotes the best total cost obtained in trial r, C * T stands for the best-known cost or the total cost associated with the global optimal solution. The comparator operator in (25) returns a value of 1 if its left-hand side (LHS) is lower than or equal to the right-hand side (RHS), otherwise returns zero. For example, X 4 expresses the number of trials at which the total cost is lower than or equal to the cost of the best-known solution plus 1e −4 . The values of C * T for the considered benchmark problems have been evaluated by the proposed algorithm using a large number of trials with adequate simulation time, as discussed in the case studies.

Benchmark problems
In order to demonstrate the powerful capability of getting the global optimal solution and the robustness of the proposed algorithm while solving the EDP, 10 benchmark problems are considered. The specifications of these problems are listed in Table 2. The first column is used to assign a number for each considered problem, and the second column indicates the problem dimension. The third column gives the total system load of each problem. In the fourth column, the practical features associated with each problem are listed using the following abbreviations: VP for valve point effects, POZ for prohibited operating zones, MF for multiple fuel options, RRL for ramp rate limits, SRC for spinning reserve constraint, and TLs for transmission losses. The fifth column in Table 2 presents references for the data accompanying each benchmark problem. Finally, the last column displays the C * T values obtained by the proposed algorithm.

Simulation scenarios
Based on selecting different TOFE values, there is an enormous number of scenarios for which the results of the proposed algorithm can be displayed. However, in most of the case studies, two scenarios have been utilized. The first one is to simulate the algorithm with a low TOFE value, which is low enough to achieve hitting the global optimal solution at least once out of the n trials considered. The second scenario is to simulate the algorithm with a TOFE value that is high enough to hit the global optimal solution n times out of n trials with a certain predefined accuracy level. GFABC-LF 1 stands for the proposed algorithm simulated with scenario 1 (sc. 1), while GFABC-LF 2 indicates the use of scenario 2 (sc. 2).

Case study 1 (Validation of the proposed improvements)
The purpose of this case study is to demonstrate the effect of each proposed improvement to the ABC algorithm. The proposed algorithm is compared with the basic ABC algorithm, GFABC −1 , which designates the proposed algorithm without Levy flights and with (11) instead of (20), and GFABC which expresses the proposed algorithm without Levy flights. Approximately, the same TOFE value is used by all the algorithms. Benchmark problem no. 5 is utilized in this case study. Figure 4 displays the summary statistics of the four algorithms using box plots. By comparing the summary statistics of GFABC −1 and ABC, it is evident that the former has lower minimum, average, and maximum cost values, which implies an improvement in the performance of the algorithm using only Equation (19). Adding (20), to produce GFABC −1 , has the effect of improving the average, standard deviation (std.) and slightly improves the minimum cost value. Finally, combining the Levy flights with the GFABC brings the minimum, average, and maximum cost values very close to each other while the minimum cost is equal to the cost of the known global optimal solution to this problem. Each algorithm is executed 100 times to produce the statistical results in this case study. To demonstrate that the improvements in the descriptive statistics are not a result of sampling errors, equivalently, these improvements did not occur by chance, the Wilcoxon rank-sum test is utilized to compute the p-value for each couple of algorithms. The p-value is used to decide whether the null hypothesis should be accepted or rejected. The null hypothesis is accepted if the p-value > .05. Table 3 presents the p-value of the Wilcoxon rank-sum test for each couple of algorithms. From Table 3, it is clear that the p-value for any two algorithms is much lower than .05. Therefore, each proposed improvement leads to a different algorithm. This conclusion and the summary statistics displayed in Figure 4 confirm the improvements added to the ABC algorithm using the proposed modifications.

Case study 2 (Effect of increasing the TOFE value and Comparison with PSO variants)
PSO algorithm and its improved variants are among the most common and successful stochastic heuristic algorithms in the literature. In general, PSO algorithms are known for their high convergence rate. In this case study, two improved variants of the PSO algorithm, which have been applied for solving the non-convex EDP, are compared with the GFABC-LF considering the effect of increasing the TOFE value on the probability of hitting the global optimal solution. The algorithms considered are PSO with both chaotic sequences and crossover (CCPSO) [15] and self-tuning improved random drift particle swarm optimization (ST-IRDPSO) [24]. The effect of increasing the TOFE value is compared for the three algorithms by plotting an improvement curve. The improvement curve plots the number of times at which the global optimal solution within a certain accuracy is obtained against the TOFE value. The TOFE value is raised by increasing both the number of iterations and population size, simultaneously, for the three algorithms. Benchmark problem no. 5 is adapted, and 100 trials are executed for each considered TOFE value within the range of 1e +3 to 10e +6 . Figure 5 shows the improvement curves for the three algorithms. The three algorithms in Figure 5 showed an increasing trend, which is hardly observable for the CCPSO and the ST-IRDPSO due to the stochastic nature of these algorithms and the very slow improvement rate associated with them. On the other hand, the value of X 4 increases rapidly for the GFABC-LF after 0.5e +6 . It continues with the increase however with smaller rates as the TOFE increases more than 1.5e +6 . The curves show that much powerful performance can be obtained by the GFABC-LF than that associated with simulating it using a small value of TOFE, and the algorithm can tackle problems with a higher difficulty level. Furthermore, compared to CCPSO and ST-IRDPSO, the improvement curve of the GFABC-LF shows more robustness against changes in the TOFE value.

Case study 3 (Demonstration of the deterministic-like solution)
The purpose of this case study is to demonstrate the deterministic-like characteristics as well as the robustness of the proposed solution method by applying it for solving seven benchmark problems while achieving an estimation of 100% probability of hitting the global optimal solution within a reasonable time. The minimum accuracy level is defined to be 1e −3 . Five hundred trials are executed to achieve a more reliable estimation of the probability than that with 100 trials. Two scenarios are used as discussed previously, and the event of hitting the global optimal solution 500 times of 500 trials is recorded with the second scenario, sc. 2. Table 4 presents the statistical results of this case study. This table demonstrates that an esti- mated probability of 100% of hitting the global optimal solution can be achieved with scenario 2 in less than 50 s for the 140-unit benchmark system.

Case study 4 (Large-scale system)
A benchmark system is a system with a complete and wellpresented data, and this system is usually adapted for comparing the results of several solution techniques presented by different publications. The largest scale benchmark system observed in the literature that satisfies the above conditions is the 160-unit benchmark system [28]. The data of this system is obtained by replicating the data of the 10-unit benchmark system in [28]. Table 5 presents the results of comparing the GFABC-LF with the conventional genetic algorithm with multiplier updating (CGA-MU), improved genetic algorithm with multiplier updating (IGA-MU), fuzzy-based hybrid particle swarm optimization-differential evolution (FBHPSO-DE),

FIGURE 6
Convergence of the total cost for the 160-unit system pseudo-inspired chaotic bat algorithm (PI-CBA), cuckoo search algorithm (CSA), one rank cuckoo search algorithm (ORCSA), and chaos mutation firefly algorithm (CMFA). In this case study, the TOFE value of the GFABC-LF is selected such that the total computational time of the GFABC-LF is nearly equal to that of the CMFA. The results in Table 5 demonstrates that the best statistical results are obtained by the GFABC-LF. Furthermore, the average cost value is approximately equal to the minimum cost value for the GFABC-LF. In Table 5, NR stands for not reported in the corresponding literature. One hundred independent trials have been used to record the statistical results of the GFABC-LF in Table 5. Figure 6 displays two convergence curves for the total cost of the 160-unit benchmark system as developed by the GFABC-LF. In Figure 6, the convergence to the minimum cost is achieved within 500 iterations. On the other hand, 320 units system can be formed by doubling the size of the 160 units system. The results of applying the GFABC-LF for this system is presented in Table 6. Increasing the TOFE value will lead to better statistical results than those reported in Table 6, but how much the TOFE value should be increased can be determined based on the required accuracy, the length of the considered dispatch interval in the real system, and the available computational power in the real dispatching centre.

Case study 5 (Highly non-convex problems)
Highly non-convex problems due to the existence of valve point effects with relatively large peaks in all the units are considered in this case study. Hitting the global optimal solution, within an accuracy of 1e −3 , 500 times of 500 trials is not as easy as that with the case study 3, in which the population size did not exceed 50 for all the problems. With highly non-convex problems, it is observed that a higher population size provides better results since it enhances the exploration capability further. Con-  Figure 5, increasing the TOFE value to 17.5 e +6 leads to hitting the global optimal solution with a probability of 98.033% (2941 times out of 3000 trials) with an average time per single run of 266.26 s. A better solution is to form a group of six parallel runs, each single run has a TOFE value of 4.2 e +6 . This achieves a probability of 99.995% with a total time of 92.984 s using the parallel computing toolbox in MATLAB on one personal computer. The previous probability is estimated as follows. With a total number of single runs equal to 6 × 500 (trials) = 3000, the algorithm succeed to hit the global optimal solution, with an accuracy of 1e −3 , 2428 times out of 3000 trials, which gives a probability of 80.933%, and therefore, a group of six single runs will have a success probability of 99.995195% as listed in Table 7, last column. The same is applicable for estimating the probability associated with a group of six parallel runs for problem 6. This approach is more accurate for estimating the probability than that based on considering 500 trials only. Table 7 displays the results, for a group of six parallel runs on a single PC for benchmark problems 5 and 6.

SENSITIVITY ANALYSIS AND PARAMETERS TUNING
To demonstrate the effect of changing each parameter value on the performance of the algorithm, problem no. 5 in Table 2 is considered. The maximum number of iterations is fixed to 2000 and the population size to 50. Twenty trials are utilized to record the statistical results. In the following paragraphs, any unmentioned value for a parameter implies that this parameter adopts the corresponding value in Table 1. To show the effect of varying the values of 1 and 2 , the algorithm performance is assessed while varying the values of both the parameters simultaneously. The considered values for 1 and 2 are 0.001, 0.005, 0.01, 0.05, 0.1, 0.5, 1, 2.5, 5, 7.5, and 10. Figure 7 plots the average cost and standard deviation values for each considered value of the parameters. Figure 7 shows that the values of 1 and and onlooker bees), the exploitation capability of the algorithm is improved. The best average cost values are obtained with values of 1 and 2 between and equal to 0.01 and 0.1. It is also observed from Figure 7 that as the values of 1,2 decrease below 0.01, the average cost starts to increase. A possible interpretation is that minimal values of 1,2 make progress towards the global optimal solution slower than required within the defined number of iterations. In Figure 7, both 1 and 2 assume the same value. The statistical results can further be improved if 1 is different than 2 . To illustrate this, Table 8 lists the average cost and standard deviation values for some cases of equal values for 1 and 2 and other cases of unequal values. The cases of unequal values are selected such that 1 is larger than 2 by 10 times the value of 2 . Comparing the results of unequal values with those of equal values indicates that a slight enhancement in the average cost and standard deviation values occurs when different values for 1 and 2 are considered. This slight enhancement may be neglected for the above example, but when 500 trials are considered, this enhancement becomes more apparent as a higher number of hitting the global optimal solution. A higher number of hitting the global optimal solution is beneficial for reducing the total number of runs, and therefore, the solution time during the real-time application.
To demonstrate the effect of changing the value of the scout bee scaling factor on the performance of the algorithm, Figure 8 plots the average and standard deviation values for different values of .
As shown in Figure 8, the best average cost and standard deviation values occur when = 0.25. Yet, this is not the best value for . After testing more values close to = 0.25, it has been noted that = 0.27 provides better statistical results for several benchmark problems compared to the value of 0.25. The results in Figure 8 confirm the need for a compromise value for , as explained in Section 4.2 of the paper.  For the limit parameter, the range of the tested values is determined based on the range of the typical limit values used with the ABC algorithm. Figure 9 displays the effect of varying the value of the limit parameter on the performance of the GFABC-LF. Unlike with 1,2 and , there is no clear trend in Figure 9, and therefore a possible simple interpretation to explain the change in the performance of the algorithm as a function of the limit value is challenging to attain. But from Figure 9, the limit value of 50 achieves the best average cost and standard deviation values followed by the value of 200. Consequently, the value of 50 is utilized for the GFABC-LF in the previous case studies. Values close to the limit value of 50, such as 40 and 60, did not improve the results obtained with the limit value of 50.
The scaling parameter in Equation (18) has been fixed to 0.001 in [10] and 0.01 in [11], but fixing the value of to 0.1 has led to better statistical results with the GFABC-LF. Figure 10 depicts the average cost and standard deviation values corresponding to the values of 0.001, 0.01, 0.1, and 1 for . Large values for such as 1 make the Levy flights very aggressive generating solutions outside the limits of the decision variables [37]. On the other hand, the parameter has been assigned the value  [16] since deviating from this value does not lead to an improvement in the statistical results of the GFABC-LF. Fifteen Levy flight cycles per each iteration are recommended in [10] for the ABC-LF algorithm; however, with the GFABC-LF, five cycles were found to provide competitive statistical results with much lower computational time than that with the case of considering 15 cycles. Figure 11 displays the average cost and the average computational time per each run for the GFABC-LF as a function of the number of Levy flight cycles. Figure 11 shows that the average cost value attained with N L = 5 is comparable with the average cost value obtained with a larger number of Levy cycles, but five Levy cycles require less computational time than a higher number of Levy cycles. On the other hand, if only one Levy cycle is used as in the Cuckoo search algorithm, a remarkable increase in the average cost is obtained.
The above sensitivity analysis justifies the parameters setting mentioned in Table 1 and supports the interpretation of the algorithm operation provided in the context of Section 4. To attain the values of the parameters reported in Table 1, extensive sensitivity analysis using several benchmark problems has been executed. The parameters setting is adjusted to provide satisfactory results considering a wide range of problems. Satisfactory results mean obtaining a deterministic-like solution within an acceptable computational time for several benchmark problems. To clarify further the tuning procedure, a set of promising parameter values has been defined for each problem based on several experiments. In each experiment, one of the parameters is improved while fixing the other parameters. Once a parameter value is improved, it is fixed, and another parameter is considered for the improvement while fixing the other parameters. The experiments have been repeated until satisfactory results are obtained with the most challenging benchmark problems. By comparing the parameters' values specified for the different problems and the range of values for each parameter, insight has been developed about the most suitable parameters setting for several problems. This creates a nonoptimal parameter setting for each problem, but an acceptable one for all the considered problems, and since the algorithm is powerful thanks to the strong balance between the exploitation and exploration and the good capability of escaping from local minimums, this parameter setting was sufficient to provide a deterministic-like solution (hitting the global optimal solution 500 times out of 500 trials) for all the considered benchmark problems.

PRACTICAL ASPECTS AND DISCUSSION
The longest simulation time is achieved with the 80-unit system, problem no. 6, since this system has valve point effects with relatively large peaks in all the units and therefore has the highest difficulty level. Compared to this system, the Korean power system with 140 thermal units, problem no. 8, has only 12 generators with valve point effects, and hence it has a much lower difficulty level. Therefore, each problem has a certain level of difficulty that can be estimated offline. In real systems with a difficulty level equivalent to that of the 80-unit system, the following strategy is recommended. The improvement curves for different load values of the system have to be plotted offline, and the worst improvement curve is selected. From this curve, a suitable value of TOFE and a number of parallel runs can be chosen to maximize the probability of hitting the global optimal solution based on the available offline computational power. After achieving the highest possible probability offline, any further increase in the number of parallel runs and the TOFE value for each single run will render the algorithm more and more reliable during the online operation. The only limitation to getting an extremely reliable algorithm is the available computational power during the online operation. Only one personal computer with traditional specifications has been used to estimate a probability of more than 99.9% including the most challenging benchmark systems; however, this will not be the case in a real dispatching centre, and much more computational power is expected to be available for systems such as those containing 80 thermal units with valve point effects.

DYNAMIC NON-CONVEX EDP AND FUTURE WORK
The dynamic non-convex EDP is much more challenging than the static non-convex EDP. The static non-convex EDP, also known as the practical EDP [28,38], is a non-linear, non-smooth, non-convex, and NP-hard optimization problem [39,40]. For the same system, the dynamic EDP implies multiplying the size of the static EDP by 24 to model the correlated generation dispatch for 24 h. Typically, in the previous literature, special strategies for generating an initial population and treating the dynamic nature of the problem are proposed besides applying any suggested meta-heuristic algorithm [41][42][43]. These strategies can have much more impact on the solution quality of the problem than the utilized meta-heuristic algorithm. For example, several meta-heuristic algorithms have been introduced before for solving the dynamic non-convex EDP such as the enhanced bee swarm optimization algorithm (EBSO) [44], modified hybrid evolutionary programming-sequential quadratic programming (MHEP-SQP) [45], and time-varying acceleration coefficients iteration particle swarm optimization (TVAC-IPSO) [46]. However, with simple meta-heuristic techniques such as the genetic algorithm and the differential evolution, the authors of [41] presented much better solutions to the dynamic non-convex EDP thanks to the proposed heuristics for treating the constraints. The proposed algorithm has been tested for solving the 10-unit benchmark problem without transmission losses. The data of this problem is available in [47]. Since it is not the main purpose of this paper to introduce a solution method to the dynamic non-convex EDP, a simple strategy is used to handle the constraints, and a random initial population is utilized. A solution for the first hour is developed by solving the static non-convex EDP for this hour. This solution is used to provide initial output power for the units. Then, a random initial population is generated for the remaining 23 h. Each individual of the population has 230 decision variables (23 × 10). The GFABC-LF has been applied for optimizing the power dispatch for the whole operational cycle while forcing the ramp rate limits by penalizing the solution that violates these limits. The obtained statistical results are presented in Table 9.
The same parameters setting used with the static variant of the problem is considered to obtain the results in Table 9, i.e. the parameters setting is not optimized for the dynamic case. However, introducing an efficient strategy for treating the dynamic nature of the problem and a procedure for creating an initial population of feasible solutions can be even more important than optimizing the parameters setting for the dynamic non-convex EDP. Several strategies for treating the ramp rate limits considering the entire operational cycle and several methods for generating a random initial population of feasible solutions for the dynamic non-convex EDP are adopted in the previous literature [41][42][43]. Finding a new strategy for treating the dynamic nature and a strategy to generate an initial feasible population for the dynamic non-convex EDP is an open field of research and may require more than one complete paper, so proposing a method for solving the dynamic non-convex EDP is out of this paper scope and is considered for future research.

CONCLUSION
Advanced optimization algorithm is proposed for solving the static non-convex EDP. In this context, the problem formulation included the valve point effects, ramp rate limits, prohibited operating zones, multiple fuel options, spinning reserve constraints, and transmission losses to obtain precise ED solution.
To develop the proposed algorithm, improving the exploitation capability of the ABC algorithm and combining the improved ABC with the Levy flights are proposed. To improve the exploitation capability, scaling factors are added to the update equations of the employed and onlooker bees to focus their search around the best solution found so far (Gbest). Furthermore, the scout bee is modified to generate random solutions within a limited volume around the Gbest. These proposed modifications establish a multi-level concentric search around the Gbest, producing a much stronger exploitation capability for the algorithm. On the other hand, five Levy flight cycles are added at each iteration to strengthen the capability of escaping from local minimums. In addition, the increased diversity in the population introduced by Levy flights enhances further the exploration capability of the algorithm. Consequently, the resultant good balance between the strong exploitation and exploration capabilities justifies the strength of the proposed algorithm. The effect of each proposed improvement added to the basic ABC algorithm has been separately investigated and validated. The proposed algorithm is capable of hitting the global optimal solution with a high frequency and is robust since it achieves this with 10 benchmark problems using one parameter setting. The algorithm shows a clear improvement rate in the performance as a function of the simulation time. Consequently, improvement curves can be plotted and exploited in maximizing the reliability of the proposed algorithm for the most challenging benchmark problems. Large-scale problems are considered. Although these problems have valve point effects and multiple fuel options in all the units, these problems did not introduce the highest difficulty level to the proposed algorithm, since the modelled valve point effects have smaller peaks than those with the 80-unit benchmark system. Finally, an estimation of 100% of hitting the global optimal solution is achieved with most of the considered benchmark problems, and more than 99.9% is obtained with the most difficult ones. Developing an appropriate strategy to handle the ramp rate limits for an entire operational cycle of 24 h and developing a suitable initialization strategy are considered for future research while extending the application of the proposed algorithm to the dynamic non-convex EDP.