Modified Cuckoo Search Algorithm for Solving Nonconvex Economic Load Dispatch Problems

This paper presents the application of Modified Cuckoo Search Algorithm (MCSA) for solving Economic Load Dispatch (ELD) problems. The MCSA method is developed to improve the search ability and solution quality of the conventional CSA method. In the MCSA, the evaluation of eggs has divided the initial eggs into two groups, the top egg group with good quality and the abandoned group with worse quality. Moreover, the value of the updated step size in MCSA is adapted as generating a new solution for the abandoned group and the top group via the Levy flights so that a large zone is searched at the beginning and a local zone is foraged as the maximum number of iterations is nearly reached. The MCSA method has been tested on different systems with different characteristics of thermal units and constraints. The result comparison with other methods in the literature has indicated that the MCSA method can be a powerful method for solving the ELD.


Introduction
Economic Load Dispatch (ELD) problem is one of the major optimization issues in power system operation.The objective of the ELD problem is to allocate the power demand among committed generators in the most economical manner while all physical and operational constraints are satisfied.The cost of power generation, particularly in fossil fuel plants, is very high, and the economic dispatch can help save a significant amount of cost [1].In the past, the problem was very simple due to considering a quadratic fuel cost function for the objective with linear constraints [2].However, it is more realistic to represent the generation cost function for fossil fuel fired generators as a segmented piece-wise quadratic function or with Valve Point Effects (VPE) [3].In addition, generating units may be supplied with Multiple Fuel (MF) sources such as gas and oil to produce electricity [4].Therefore, the fuel cost function of thermal units with multiple fuels can be represented as a piecewise nonlinear function.Moreover, the generators should also be satisfied with operating conditions such as upper and lower generation limits and prohibited operating zones [5].Over the past decades, a large number of methods using mathematical programming were widely applied to the ELD problem such as Dynamic Programming (DP) [2] and [6], Newton's method [2], lambda iteration method [7], gradient method [8], Linear Programming (LP) [9], Lagrangian relaxation algorithm [10], and Quadratic Programming (QP) [11].However, the applicability of the conventional methods to the problem is limited to systems with a convex objective function.Consequently, these methods can only be applied to the systems where the cost function of each generator is approximately represented by a simple quadratic function and the effects of valve-points are ignored [12].Recently, several artificial intelligence-based methods have been developed for solving the ELD problems such as Evolutionary Programming (EP) [1], [13] and [14], Hopfield Neural Network (HNN) [15] and [16], Genetic Algorithm (GA) [17], [18] and [19], Biogeography-Based Optimization (BBO) [20], Particles Swam Optimization (PSO) [21] and [22], Differential Evolution algorithm (DE) [12] and [23], and Harmony Search (HS) [24].The GA and EP methods have many advantages such as global search capability, robust and effective constraints handling capacity, reliable performance and minimum information requirement, making them a potential choice for solving the ELD problems [12].However, the disadvantage of GA is that it is sensitive to the selection of control parameters [25].HLN method can be implemented on large-scale problems but it still suffers many drawbacks such as local optimum solution, long computation time, and difficult implementation [26].PSO is efficient for solving optimization problems; however, the performance of the traditional PSO greatly depends on its parameters and it often suffers from the problem of being trapped in local optima [27].
The potential of DE is the fast convergence characteristic [28] and the advantage can yield higher probability of searching toward a local optimum or getting premature convergence.This drawback could be overcome by employing a larger population.However, larger population will lead to longer time to estimate he fitness function [29].For HS [24], the jazz improvisation seeks to find a musically pleasing harmony similar to the optimum design process which seeks to find an optimum solution.However, HS uses a stochastic random search instead of a gradient search and the derivative information is unnecessary.In addition, the HS may get the premature convergence in the performance [30].In fact, each method has advantages and disadvantages for implementation to solve the ELD problems.Therefore, combined methods have been suggested to take advantages of each method and enhance their search ability.Several hybrid methods have been proposed for solving the ELD problems such as hybrid GA-PSO method [31], Hybrid Stochastic Search (HSS) [32], hybrid PSO-SQP [33], and hybrid genetic algorithm [34].
Generally, these hybrid methods can obtain better solution quality than each member method.However, these methods suffer a difficulty of the proper selection of many controllable parameters for dealing with different problems.
Cuckoo Search Algorithm (CSA) is a recently developed meta-heuristic algorithm inspired from the reproduction strategy of cuckoo species in the nature for solving optimization problems [35].CSA has been widely and successfully applied to many optimization problems in engineering such as economic load dispatch [35] and hydrothermal scheduling [36] since it uses the Levy flights as a high class random walk to generate a new solution, and new solution generations are performed in each iteration.The effectiveness of the CSA method is better than GA and PSO on many benchmarked functions given in [35] and better than PSO, DE, GA, BBO shown in [26].However, the conventional CSA still suffers a slow convergence for complex and large-scale systems due to the disadvantage of fixed value of updated step size parameter.Therefore, based on the drawback of CSA, a new modified CSA (MCSA) method has been developed to enhance its search ability and improve solution quality [37].The MCSA method is developed fulfilling the changes in the first new solution generation via the Levy flights to improve the quality of the solutions.The evaluation of eggs has divided the initial eggs into two groups, the top egg group with good quality and the abandoned group with worse quality.The value of the updated step size in MCSA is adapted as generating a new solution for the abandoned group and the top group via the Levy flights so that a large zone is searched at the beginning and a local zone is foraged as the maximum number of iterations is nearly reached.In addition, each new solution in the top group is also generated by using the newly obtained one from each two eggs, and the new one in the abandoned group is produced by the newly gained one from the best egg and each individual egg.As a result, the MCSA has been demonstrated to be more efficient than the conventional CSA as they have been tested on many benchmarked functions [37].
In this paper, MCSA is applied for solving the ELD problems considering generator and system constraint characteristics such as valve point loading effects, multiple fuel options, prohibited operating zones, ramp rate constraints, spinning reserve, and power losses.The MCSA method has been tested on different systems with different characteristics of thermal units and constraints.The results obtained by MCSA have been compared to those from other methods available in literature.

ELD with Quadratic Fuel Cost and Valve Point Effects of Thermal Units
In the classical ED problem, the objective of the problem is written as: when representing the characteristic of power outputfuel cost of thermal units, there are two subcases forming its mathematical function: the quadratic fuel cost function where valve point loading effects are neglected, and the nonconvex fuel cost function where the effects are taken into account.The former is approximately represented as in the following Eq.( 2), whereas the latter is built as in the Eq. ( 3) below [2] and [3].
where N is the number of generators; P i is the real power output of generator i; P i,min is the minimum power output of unit i; and a i , b i , c i , e i , f i are fuel cost coefficients of unit i.
Subject to:

1) Real Power Balance
The total real power output of generating units satisfies total load demand plus system power losses: and the total power loss P L is calculated using Kron's formula [36].

2) Generator Capacity Limits
The power output of each unit must be satisfied by: where P i,min and P i,max are maximum and minimum power outputs of unit i.

ELD with Multiple Fuel Options and Valve Point Effects of Thermal Unit
Since the thermal units can be driven by using Multiple Fuels (MF) and the Valve Point loading Effects (VPF) are also considered, the fuel cost function can be represented as follows [19]: where a ij , b ij , c ij are fuel cost coefficients for fuel type j of unit i and P ij,min and P ij,max are lower and upper limits for fuel j of unit i, respectively; and e ij , f ij are fuel cost coefficients for fuel type j of unit i reflecting valve-point effects.

1) Prohibited Operated Zones
Prohibited operating zones are infeasible domains that thermal units are not allowed to work with.In fact, each thermal unit can have several prohibited operating zones considered during its operation.Along the power output-fuel cost curve, several infeasible zones are included and their power output must be out of the zones.Consequently, this difficulty is also a challenge to validate the performance of an optimization algorithm, and it is considered in the paper.
For generating units with POZ, their entire feasible operating zones are decomposed in feasible sub-regions and their feasible operating points should be in one of the sub-regions as follows: where n i is the number of prohibited zones of unit i; Ω is the set of units with POZ; and P u ik and P l ik are upper and lower bounds for prohibited zone k of unit i.

2) Spinning Reserve Constraint
For safety operation, a spinning reserve for the system is required.The spinning reserve constraint for all units is defined as: where the operating margin of each generating unit S i is determined by: where S i,max is maximum spinning reserve contribution of unit i, and S R is the total system spinning reserve requirement.

Modified Cuckoo Search Algorithm
The conventional CSA is comprised of two new solution generations via Levy flights and via the action of alien egg discovery.The new solutions generated via Levy flights are obtained as below [26]: where X best and X i are the best eggs, and the i th egg among the number of eggs; α > 0 is an updated step size.The value of α has a significant influence on the final solution because it will lead to different new solutions as it is set to different values.If this parameter is set to a high value, there is a huge difference between the old and new solutions, and the optimal solution is either obtained fast or omitted.As the current iteration is high, the new obtained solution should be searched nearby the previous solution.However, the method has to find the optimal solution in a large search zone for this set value which may not reach the best optimal solution.Based on the analyzed drawback of the conventional CSA, it is clearly better to search the optimal solution in a small zone as the iteration counter is increased to the maximum number of iterations, which is predetermined for the iterative process [37].The modified version of CSA was developed by Walton et al in 2013 [37] by focusing on the updated step size and improving the quality of all solutions which are generated via Levy flights.In the modified CSA, before applying the first new solution generation via Levy flights, all the eggs are ranked and classified into a top group with better quality eggs, and an abandoned group with worse quality eggs.The abandoned group only focuses on the step size α, which is a variable with the value decreased as the current iteration is increased.It is more complicated when a new solution in the top group is newly generated because it needs information exchange between each two eggs, one is randomly picked and one is picked in order.There are three cases of the two picked eggs as follows: • the chosen egg and random egg are the same, • both eggs possess the same fitness, • the random one has lower or higher fitness than the predetermined egg.
In addition, each new egg in the abandoned group is also produced by the newly gained one from the best egg and each individual egg, and the updated step size during the search process for the abandoned group is adapted to improve the quality of each solution.

New Solution Generation for MCSA
As described above, in the Modified CSA all nests are first sorted in the descending order based on their fitness function value and then classified into two groups.1) The First New Solution Generation via Levy Flights • Generation of new solution for the abandoned group: Based on the modification applied to the abandoned eggs (d = N otop+1 , . . ., N p and N otop and N p are respectively the number of eggs in the top group and in the initial population), the optimal path for the Levy flights is calculated using Mantegna's algorithm as follows: where rand 1 is the distributed random number in [0, 1], the step size α is determined by where 1/ √ G is the number current iteration, and ∆X is obtained by: • Generation of new solution for the top egg group: The modification applied to the eggs in the top group (d = 1, • • • , N otop ) is described in Subsection 3.1.The optimal path for the Levy flights is calculated using Mantegna's algorithm as follows: where rand 2 is the distributed random numbers in [0, 1].The value of α and ∆X will be determined depending on considered cases as follows: -Case 1 : The same egg is picked twice: where α = 1/G 2 .-Case 2 : Both eggs have the same fitness value function: -Case 3 : The random egg has lower fitness than egg d, or the random egg has higher fitness than egg d, For Cases 2 and 3, α is set to 1 and ϕ = (1+ √ 5)/2.

2)
The Second New Solution Generation via Discovery of Alien Eggs Similar to the conventional CSA, the second new solution generation via discovery of alien egg is also employed in the Modified CSA but all eggs of the top group and abandoned group are integrated into one group first.The new solution due to this action can be found as follows: where K is the updated coefficient determined based on the probability of a host bird to discover an alien egg in its nest [26] and ∆X dis d is the increased value [26].

Calculation of Power Output for Slack Thermal
To guarantee that the equality constraint shown in Eq. ( 4) is always satisfied, a slack generating unit is arbitrarily selected and therefore its power output will be dependent on the power output of remaining N − 1 generating units in the system.The detail calculation of slack variables can be found in [37].

Implementation of MCSA for ELD
The MCSA method is implemented for solving ELD problem as follows.
Initialization: A population of N p host nests is represented by X d (d = 1, . . ., N p ) where each solution vector of variables given by X d = [P 2,d , P 3,d , . . ., P N,d ].The power output of the thermal units in N p nests are randomly initialized satisfying P i,min ≤ P i,d ≤ P i,max .
Based on the initialized nests, the fitness function to be minimized corresponding to each nest for the considered problem is calculated as: where K s and K r are penalty factors; P 1,d is the power output of slack thermal unit 1; and S di is calculated from Eq. (9) and Eq. ( 10).
The limit for slack thermal unit 1 in Eq. ( 20) is determined as follows: where P 1,max and P 1,min are the maximum and minimum power outputs of slack thermal unit 1, respectively.

1) Evaluation of the New Solution
All the initial nests are evaluated, ranked in descending order and newly generated via Levy flights in Subsection 1) are redefined as they violate the bounds as follows: otherwise. ( The slack thermal unit P 1,d is calculated from Subsection 3.3.The fitness function value of the new egg is calculated using Eq. ( 20) and then compared to that from the old egg.The egg with better fitness function value is considered as the new solution The second new solution generation is performed by the action of alien egg discovery in Subsection 2) and the new solution is then redefined by using Eq. ( 22).The slack thermal unit 1 is also obtained by Subsection 3.3.The value of the fitness function is recalculated using Eq. ( 20) and the nest corresponding to the best fitness function is set to the best nest X best of the population.

2) Stopping Criterion
In the MCSA, the stopping criterion is based on the maximum number of iterations.That means the algorithm is terminated as the current iteration is equal to the maximum number of iterations.
Generally, the stopping criteria for methods solving optimization problems are usually based on iterative error of two consecutive iterations, constraint mismatch, and maximum number of iterations.In fact, depending on the application of different solution methods for solving optimization problems, the stopping criteria may be used in different ways as long as the obtained final solution is a feasible one.As mentioned in the paper, the stopping criterion of the proposed MCSA method is based on the maximum number of iterations.The proposed MCSA method handle equality constraint in Eq. ( 4) by using slack variables so that the equality constraint is always satisfied by this c 2016 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING way.For inequality constraint in Eq. ( 5), the proposed methods can handle them in fitness function described in Eq. ( 20) and limit check after each iteration to guarantee that the obtained solution is a feasible one.The iterative error of two consecutive iterations is not consider for the proposed method since it is the population based method using random search and it may happen that the obtained solution for many iterations is not improved.That means, the solution obtained from many iterations is the same and it cannot be used as a stopping criterion since it may lead to early termination of algorithm with non-optimal solution.In fact, the stopping criteria of population based methods are always based on the maximum number of iterations.However, a small number of iterations may lead to non-optimal solution.On the contrary, a large number of iterations will lead to long computational time.Therefore, for a proper selection of maximum number of iterations for each system, there should be a need of experiments.Above all, optimal solutions from different methods can be used for comparison to each other as long as the solutions are feasible; that means all problem constraints are satisfied and different stopping criteria may be used.Therefore, the optimal solutions from population based methods can be compared to those from conventional methods.This comparison is very popular in literature.

3) Overall Iterative Algorithm
The overall procedure of the proposed MCSA for solving ELD problem is described as follows. • Step 1: Select parameters for MCSA.Initialize population of host nests.
• Step 2: Calculate the slack thermal unit and set iteration counter G = 1.
• Step 3: Evaluate the fitness function to choose X abandoned d , X top d , X topr and X best .
• Step 4: Generate new solutions for abandoned eggs via Levy flights.
-Calculate thermal power outputs -Calculate the fitness function in Eq. ( 20).
• Step 6: Generate new solutions for top eggs via Levy flights.
-Calculate thermal power outputs.
• Step 8: Put new eggs generated in Steps 4 and 6 in a group of eggs.
-Calculate all thermal power outputs.
-Evaluate fitness function to choose new X best .
• Step 11: If G is less than the maximum number of iterations, G = G + 1 and return to Step 3. Otherwise, stop.

Numerical Results
The MCSA has been tested on different systems corresponding to the formulated problems including 13unit system considering valve point loading effects with two load demands of 1800 MW and 2520 MW, 20-unit system with quadratic cost function and transmission losses, systems up to 160 units considering valve point loading effects and multiple fuel options, and systems up to 90 units considering prohibited operating zones and spinning reserve.The algorithm is coded in Matlab platform and run 100 independent trials for each case on a 1.8 GHz PC with 4 GB of RAM.

13-Unit System with Valve Point Loading Effects
The system consists of 13 thermal units with the valvepoint effects from [1,33] supplying to two load demands of 1800 MW and 2520 MW.The power losses for this system are neglected.For implementation of MCSA, the number of nests and the maximum number of iterations are respectively set to 24 and 10000 for the two load demands whereas the probability pa is tuned in the range from 0.1 to 0.9 with a step of 0.1.
The obtained results including minimum cost, average cost, maximum cost, standard deviation cost and average computational time from the MCSA method with different values of pa are given in Tab. 1 and Tab. 2. As observed from the tables, the MCSA method can obtain the best solution at p a = 0.8 for the load demand of 1800 MW and p a = 0.8-0.9 for the load demand of 2520 MW.The minimum cost and computational time from the MCSA method are compared to those from other methods as in Tab. 3. Obviously, the MCSA method can obtain better best total cost c 2016 ADVANCES IN ELECTRICAL AND ELECTRONIC ENGINEERING than other methods for the both cases.Moreover, the computational time from MCSA is also faster than that from other methods.[19] 17963.9858.28 --HGA [34] 17964.8115.33 24172.25 15.33 GA-PSO [31] 17968.50---PSO-SQP [33] 17969.9333.97 24261.05-DE [12] 17963.83-24169.92-QPSO [21] 17969.01---HSS [

20-Unit System with Quadratic Cost Function and Transmission Losses
The test system having 20 thermal units with quadratic fuel cost and transmission losses from [39] are considered in this Section.The parameters of MCSA including the number of nest, the maximum number of iteration and the probability pa for this system are set to 12, 500 and 0.4, respectively.The obtained result including minimum cost, average cost, maximum cost, standard deviation cost and computational time are given in Tab. 4 and the result comparison from MCSA to other methods is shown in Tab. 5.As observed, the MCSA method can obtain better solution quality but take longer computational time than Lambda-iteration method [39], Hopfield model [39] and BBO [20].However, the computational time by all methods for the system is less than one second.

Systems with Multiple Fuel Options and Valve Point Loading Effects
In this Section, the MCSA method is tested on the systems with 10, 20, 40, 80 and 160 units in [19].The number of nest and probability pa are respectively set to 12 and 0.1 for all systems and the maximum number of iterations is respectively set to 1000, 2000, 3000, 5000 and 8000 for the 10-unit, 20-unit, 40-unit, 80unit and 160-unit systems.The minimum cost and computational time obtained by the MCSA method and the CSA [26] are given in Tab. 6. Obviously, the MCSA method can obtain better solutions than the CSA method [26] for the systems.However, the computational time from the CSA method is faster than the MCSA method.Note that the CSA [26] is performed on a 2.1 GHz PC with 2 GB of RAM.Moreover, the average costs obtained by the MCSA method for the systems are also compared to those from CGA_MU [19] and IGA_MU [19] as given in Tab. 7. As observed from the table, the MCSA method is better than other methods in terms of fuel cost and computational time.Note the computational times for both CGA and IGA AMUM methods were run on a PIII-700 PC.The test systems from [17] comprise 15, 30, 60 and 90 units.The maximum number of iterations for the MCSA for the systems is set to 500, 2000, 2500 and 3000 for the systems, respectively.The probability pa is fixed at 0.25 for all systems.The results obtained by the MCSA method and CSA [26] are given in Table 8 to compare the effectiveness between the two methods.
Tab. 9: Result comparision for the systems with prohibited operating zones and spinning reserve.Obviously, the MCSA method can obtain approximate or better cost but longer computational time than CSA [26] for all systems.In addition, Tab. 9 shows a comparison of the average total cost and computational time obtained by the MCSA for the systems to those from CGA and IGAMUM in [17].As observed, the MCSA can obtain better average cost and computational time that CGA [17] and IGAMUM [17] where the computational times for CGA and IGAMUM were from a PIII-700 PC.

Conclusion
In this paper, the MCSA method has been successfully applied for solving ELD problems with different objective functions such as quadratic fuel cost function, nonconvex fuel cost function, and multiple fuel cost function of thermal units considering different constraints such as transmission losses, prohibited operating zones and generation limits.The main modifications of the MCSA method based on the conventional CSA method are classification of nests in two groups based on their fitness function value and the information exchange in each group to enhance its search ability and speed up convergence process.The advantages of the MCSA over other considered methods are few control parameters, fast convergence and high quality solution.The strong points have pointed out by testing on several systems and the obtained results have indicated that the MCSA method is more effective than the compared methods.Therefore, the MCSA can be a very favorable method for solving ELD problems The nests with high fitness function value X abandoned d are put in abandoned group, and the other ones X top d are put in a top group.A nest which is randomly picked among the X top d nests is called X topr .Besides, another one with the best quality is named X best d .Two new solution generations are respectively obtained as below.
Tab. 1: Result obtained by MCSA for the 13-unit system with load demand of 1800 MW by different values of pa.Tab.2: Result obtained by MCSA for the 13-unit system with load demand of 2520 MW by different values of pa.
Result comparision for the 20 untis system.
Tab. 7: Result comparision for average cost of the systems with multiple fuel option and value point loading effects.