Performance of Various Metaheuristic Techniques for Economic Dispatch Problem with Valve Point Loading Effects and Multiple Fueling Options

The paper presents the metaheuristic approaches based on pattern search and simulated annealing hybridized with sequential quadratic programming, a powerful nonlinear adaptive technique for estimation of the finest combination of generated power in a given system at lowest operating cost while sustaining the operating condition of system efficiently. The fuel cost is minimized by satisfying the nonlinear operating conditions of thermal units mainly based on generation capacity constraints, generator ramp limit, power balance constraints, and valve point loading effect and by keeping in view the prohibited operating zones, respectively. About the optimization, a comparative study is made for pattern search (PS) and simulated annealing (SA), as a viable global search technique and sequential quadratic programming, an efficient local optimizer and their hybrid versions. The applicability, stability, and reliability of the designed approaches are validated through comprehensive statistical analysis based on Monte Carlo simulations.


Introduction
The performance and designing of electric power system serve two purposes; one is to maintain the economy and the other is to include reliability of the system [1].Economic dispatch problem (EDP) is one of the significant optimization problems in industrial independent power system and acquires increasing importance as the operating cost and electric energy demand is increasing globally around the planet [2,3].The procedure involves dividing the power demand of generating system for online thermal units in such manner that their operating cost is optimal while satisfying the power demand to the user and constraint adequately [3,4].Since the operating cost of generating unit is inflated an optimal EDP can lead to prominent cost saving with reducing the atmospheric emission of thermal power plant [4].In the energy crises age, a number of engineers and scientists are taking interest in the optimized solution of this industrial problem [5] by keeping in view the constraints like generator ramp limit, power balance, and valve point loading effect (VPLE) and prohibited operating zones [6,7].However the facts of EDP with multiple fueling option (MFO) and MFO combined with VPLE conditions is not well addressed simultaneously, which are the motivations for taking this study with metaheuristic computing approach.
Traditional numerical models and programming strategies have been developed, like lambda iteration method, quadratic programming, Lagrangian relaxation algorithm, base point participation factor, gradient method linear programming, reduced gradient method, piecewise linear cost function, Newton method, and dynamic programming [8,9], but all these EDP tools are flexible and applicable if the function between operating cost and fuel or heat are piecewise linear and increasing categorically [10].Modern approaches use artificial neural network models which show dominance on classical techniques [11], while its parameters can overwrite during the learning in complex constraints based EDP problems.In order to overcome this issue the

Economic Dispatch Problem Formulation
The objective of EDP is to find the finest combination of generated power in a given system at lowest operating cost while sustaining the operating condition of system efficiently by satisfying the nonlinear operating conditions formulated as follow.

Fuel Cost Minimization.
Assume we have a system where n thermal unit committed generating real power   and total fuel cost is   (  ); usually the fuel cost curve of committed units is assumed to be quadratic function of the real power output of generators.Consider where NG is the number of committed generators and   is real output power of generator in megawatt (MW) and can be

Generation Capacity Constraint.
For the stable operation of each thermal unit the power generated, that is,   , must be restricted by minimum power that must be generated and maximum power that a unit can deliver as follows: where  0  , UR  , and DR  are the previous operating point, upper ramp limit, and down ramp limit of the th generator, respectively.

Power Balance
Constraint.Usually all the power plants are not working on same place they are spread geographically and also they must have some transmission power losses.The total power generated by committed thermal unit must be equal to the power demand of electric network plus the losses associated in power transmission [22]: where   is the power demand, NG is number of generators, and   is the transmission power loss collective equation that is also known as coordination equation.The most approximate and important method for transmission network losses as a function of generator power is  coefficient method based on fact that, under normal operating conditions, the transmission loss is quadratic in the injected bus real power.Equation (4b) represents the general form of the loss formula using  coefficient method.Here NB,     ,     , and   are the number of buses, the real power, reactive power of th and th buses, and elements of impedance matrix, respectively.However in this study  coefficient method will be used so   are the loss constant of coefficient .

Valve-Point Loading Effect (VPLE).
For a power plant having NG committed units have multivalve opening sequence to meet the fluctuating loads demand.VPLE introduce ripples in the input-output characteristics of generating units resultantly making the heat curve rate nonconvex with increase in unit loading.Heat rate curve decreases between the opening points of two valves as shown in Figure 1 [23].VPLE introduces ripples in the heat rate curve in multi valve opening turbine as shown in Figure 1, and for a perfect modeling of EDP a filter in the form of sinusoidal function must be added in objective function in this paper as given below: where   and   are the coefficient of th unit that represents valve-point effect.

Prohibited Operating Zones.
According to [3] practical generating units contain forbidden operating regions also known as "prohibited operating zones" due to machinery errors such as boiler pumps motors faults which ultimately leads input-output curve nonlinear and this can be modeled as follows: where   is the number of prohibited zones in the th generator curve,  is the index of prohibited zone, and    and    are lower and upper limit of th prohibited zone of the th generator.

EDP with Multiple-Fueling Option (MFO).
Another significant consideration involved in designing of EDP optimal cost is unit multiple fueling.A thermal unit has multiple cost rate curves because unit is fueled by different sources and needs to operate the unit on lower curve which result in producing hybrid cost function (HCF) and HCF contains information about unit operation or type of fuel used in generating operation.Practically multiple fuel sources are supplied to the thermal units and a single thermal unit can be represented by multiple quadratic cost functions; each cost function represents the type of fuel used and each unit must pin point the economic fuel to burn and fuel cost function formulation is as follows [23]: where   ,   , and   are the cost coefficient of generator th using  fuels.

EDP Considering Both VPLE and MFO.
To obtain the realistic and transparent EDP solution in this paper both MFO and VPLE model are considered.Mathematical modeling for this consideration is combination of both cost functions given in ( 6) and (8).Consider

Proposed Methodology
The brief introduction of optimizers used for optimization of EDP is narrated along with the necessary detail of hybrid algorithm in the form of procedural steps.
3.1.Simulated Annealing.Simulated annealing is heuristic as well as probabilistic method proposed by Scott Kirkpatrick, C. Daniel Gelatt, and Mario P. Vecchi in 1983 for localizing better approximation to find global minimum of given cost function that may be surrounded with several local minima in a large search space [24].SA is more efficient as compared with hill climbing and is a good choice for engineers and researchers to handle complex nonlinear systems especially problems like EDP due to its coding leniency to get guarantee finest optimal solutions.The method works on the principle of randomly initializing of finite solution space of specific size  on the basis of scattered neighborhood points that mimic from the principle of thermal annealing of solids [25].The decreasing of temperature at a specific rate in each iteration resembles the optimization of the merit function of EDP subject to varying the heat values as adaptive parameters of a specific problem.The candidate solution is updated based on the values of the fitness in terms of mean square error of the previous possible solution along with the index term; on this basis the iterative process matures towards the region of convergence of the fitness function.SA belongs to the family of global search techniques so it has various start points in the start of the algorithm to avoid local minimum.Usually the number of iterations of SA is very large while the computational budget is quite reasonable for most of the test bench problem.Although, it has been applied in various areas of engineering, some key applications are found in crystallographic refinement, neural computing [26], optimal reconfiguration of distribution networks, wireless communication, phase balancing, van der pol oscillators involved in mechanical engineering [27], and combinatorial optimization problem [28].The generic flow diagram used for the simulated annealing algorithm is provided in Figure 2.

Pattern Search.
Pattern search (PS) is numerical optimization method proposed by Hooke and Jeeves [29] for statistical problems and main advantages of PS are that it is efficient and computational structure is simple to understand and also it does not required gradient of cost function so it can be used in functions that are not constant or not differentiable and method is also known as direct-search, derivative-free approaches.Pattern search algorithm searches for a sequence of points to reach an optimal solution by initializing a set of points called mesh.The mesh formation procedure is involved by adding a pattern vector to the initial point and if PS acquires a point in mesh that refines the objective function then the new point becomes the current point at the next step of the algorithm [30].
Usually the practical nonlinear optimization problem of power system lies outside the reach of the standard optimization techniques.The practical EDP problems can be optimized by PS which is computationally efficient and easy to implement as compared to other heuristic techniques [31].PS initializes by computing sequence of point that may or may not be reached optimally and established a set point called mesh around these given points.This current point could be a starting point provided by user or extracted from previous step of algorithm.Current points are then added to a pattern vector which is scalar multiple of a set of vectors for mesh formation.If an objective function is improved by point in mesh then the new point becomes the current point at the next iteration [32].
Pattern search is efficient optimization techniques applied in areas of parallel computing, pattern character recognition, graphic user interference thermal control systems, nonlinear and optimal control theory [33], and so forth.

Sequential Quadratic Programming.
The SQP algorithm is one of the nonlinear programming techniques used extensively by the research community for constraints optimization problems.The strength of the scheme is well-established in terms of efficiency, accuracy, and percentage of successful solutions over a large number of test problems.Thum and Schilling [34] have provided excellently a brief introduction and description of SQP algorithms.As the physics of SQP is based on local search techniques so it is more probable to get stuck in the local minimum so some of the iterations during Monte Carlo simulation can show divergence or least accuracy.However the rate of convergence and level of accuracy for successful runs is dramatically high for SQP algorithms.The detailed history, mathematical formulation, importance, and application related with SQP methods can be seen in [35,36].The SQP methods have been applied on broad disciplines of engineering and applied science from the day of their birth to today; for instance, few recently published articles addressing SQP technique are referred to for interested readers in [37,38].Now the intention is to formulate the overall procedure and the optimization of the unknown designed parameters of the units generation along with the fuel cost; the process of the proposed scheme is provided in Figure 3.The modeling is performed based on relations provided in (1) and relations of constraints, respectively.The optimization is performed using global versions of simulated annealing and pattern search hybrid with sequential quadratic programming, an efficient local search technique.In Figure 3(b) the designed parameters obtained from global optimizer are provided as a start point to SQP for further tuning of the design.The algorithm used for the optimization of EDP in the form of the logical steps based on hybrid computational technique is as follows.
Step 1 (population generation).A random and well scattered population of size  is generated with agents of length exactly equal to the number of power generating units .
Step 2 (initialization).Initialize and assign the general and specific values for parameter setting involved in pattern search and simulated annealing algorithm as provided in Table 1.Step 3 (fitness evaluation).Find the fitness of each agent using the fitness evaluation function  as given below in the form of mean square error (MSE): where Ĉ op is the approximated operating cost for th generator.
Step 4 (termination criteria).The algorithm is terminated on the basis of the following criteria.Step 5 (reproduction).New population is generated using the operations evolved in PS and SA as per the setting given in Table 1 and then go to Step 3.
Step 6 (refinement via SQP method).One of the best results obtained by PS and SA is provided to SP as a start point for further refinement, and the parameter setting and values used for SQP algorithm are provided in Table 2.
Step 7 (storage).Store the global best of this cycle and repeat the procedure from Steps 3 to 7, to get sufficient numbers of independent runs to have sufficient data for performing analysis.
Step 8 (statistical analysis).The statistical analysis is made on the basis of four fundamental performance modes like best, worst, mean, and standard deviation and skewness.

Simulation and Results
In this section, a scenario based simulation is performed based on three case studies composed of thirteen, fifteen, and twenty generating units, respectively, for EDP problem in which the objective functions were ramp rate limit, prohibited operating zones in the power system operation; MFO and VPLE and transmission losses are also employed to demonstrate efficiency of proposed schemes using MATLAB built in subroutines of the optimizers like PS, SA, SQP, PS-SQP, and SA-SQP, respectively, whose parameters setting and values are provided in Tables 1 and 2. The reference data sheet for thirteen thermal units [37] with load demand of 3080 MW is provided in Appendix A. The data set used for fifteen thermal units [19][20][21] along with generating unit capacity, coefficients, ramp rate limits, and prohibited zones is provided in Appendix B while the load demand expected to be determined is 2630 MW.However, the data sheet for twenty thermal units of generating unit capacity and coefficients [38] are provided in Appendix C, where system supplies a total load demand of 2500 MW.One of the best designed parameters is in fact the power generating values of generators and is presented in Figure 4 at semilog scale in order to widen the small effects in the power values.The lowest cost of the power generation can be obtained by taking the generating values from the hybrid schemes based on PS-SQP and SA-SQP, respectively, from Figures 4(a), 4(b), and 4(c) in fitness evaluation function provided in (10) along with the subsequent relations given in (1) and expressions of EDP constraints.The -axis in Figure 4 defines the number of thermal units taken in the study and axis demonstrates the solvers PS, SA, SQP, PS-SQP, and SA-SQP, respectively, while -axis provides the optimized values of the power for thermal units.
The optimized final values along with computational time for PS, SA, SQP, PS-SQP, and SA-SQP are provided in Table 3.The time analysis carried out for this paper is on Intel(R) Core (TM) i5-3210M CPU @ 2.50 GHz, 2.5 GHz, 4.0 GB of RAM, and running with MATLAB version 2012a.It is quite evident from the table that level of accuracy is dominated in case of hybrid schemes as compared to global search and local search techniques separately.However the computational budget of hybrid approach is a little bit more as compared to other presented schemes but this effect can be over shed due to the high accuracy attained by hybrid schemes.The local search scheme seems to be outperformed in terms of both accuracy and computational cost; however one cannot rely on one of these best values of SQP as the local search techniques are more probable to get stuck in the local optima due to bad initial guess.Another observation is that as the number of  generating units increases the level of accuracy reduces while computational time increases due to increase in the EDP complexity.
The reliability of the stochastic algorithms is being tested by a comprehensive statistical analysis in terms of minimum of cost (Min cost ), maximum value of cost (Max cost ), mean cost value (Mean cost ), and standard deviation of the cost (STD cost ).The analysis is not only performed on cost value but also performed on computational time needed for the optimization in term of mean CPU (Mean CPU ) time.The results of the Monte Carlo simulations for stochastic solvers are narrated in Tables 4 and 5 for fifteen and twenty units, respectively, while for thirteen units the same analysis is provided in Appendix D.
It is quite evident from Tables 4 and 5 that the results of the proposed schemes are in good agreement in terms of reliability, effectiveness, economical cost, and computational budget required in terms of time complexity; however the hybrid approaches PS-SQP and SA-SQP show more supremacy than that of PS, SA, and SQP.The computational budget of hybrid schemes seems to computationally heavy of order of just few seconds; however this aspect seems to be minimal at the cost of level of accuracy achieved.
As the proposed scheme is based on heuristic computation one of the best results cannot define the effectiveness and reliability of the schemes; therefore in order to see the effectiveness of given stochastic methods the results are obtained for hundred independent runs rather on a single successful run.The fitness given in (10) obtained on all independent iterations is plotted in Figure 5 for three variations of generating units.
The percentage of convergence is also calculated by fixing a certain level of the threshold ( ≤ 10 −10 ) in the level of fitness obtained by all proposed algorithms and the results  9.

Figure 3 :
Figure 3: Process of proposed solution (a) overall procedure and (b) optimization of fitness function using hybrid computation.
(a) The predefined fitness value is achieved.(b) The total number of generations is executed.(c) Predefined value of the function tolerance is achieved.(d) The desired function evaluations are performed.If any of the above criteria are met then go to Step 6 otherwise simulate Step 5.

Figure 4 :
Figure 4: Analysis of the optimizers for 100 independent runs for (a) thirteen generating units, (b) fifteen generating units, and (c) twenty generating units.

Table 10 :
Optimal solution of the thirteen, fifteen, and twenty thermal units, respectively, using solvers involving PS, SA, SQP, PS-SQP, and SA-SQP.

Table 1 :
Parameter setting/values using PS and SA.

Table 2 :
Parameter setting/values using SQP method.

Table 3 :
Fitness evaluation values and their respective computational time for proposed schemes.

Table 4 :
Reliability based on 100 independent runs for fifteen units system with   = 2630 MW.

Table 5 :
Reliability based on 100 independent runs for twenty units system with   = 2500 MW.