Optimal power flow solution of power system incorporating stochastic wind power using Gbest guided artificial bee colony algorithm

Electrical Power and Energy Systems 64 (2015) 562–578

Contents lists available at ScienceDirect

Electrical Power and Energy Systems journal homepage: www.elsevier.com/locate/ijepes

Optimal power flow solution of power system incorporating stochastic wind power using Gbest guided artificial bee colony algorithm Ranjit Roy ⇑, H.T. Jadhav Dept. of Electrical Engineering, S. V. National Institute of Technology, Surat, India

a r t i c l e

i n f o

Article history: Received 4 August 2013 Received in revised form 24 June 2014 Accepted 6 July 2014

Keywords: Gbest guided artificial bee colony Optimal power flow Weibull probability distribution function Wind power

a b s t r a c t This paper focuses primarily on implementation of optimal power flow (OPF) problem considering wind power. The stochastic nature of wind speed is modeled using two parameter Weibull probability density function. The economic aspect is examined in view of the system wide social cost, which includes additional costs like expected penalty cost and expected reserves cost to account for wind power generation imbalance. The optimization problem is solved using Gbest guided artificial bee colony optimization algorithm (GABC) and tested on IEEE 30 bus system. The simulation results obtained using proposed method are compared with other methods available in the literature for a case of conventional OPF as well as OPF incorporating stochastic wind. Subsequently an extensive simulation study is conducted to investigate the effect of wind power and different cost components on optimal dispatch and emission. Numerical simulations indicate that the operation cost of system and emission depends upon the transmission system bottlenecks and the intermittency of wind power generation. Ó 2014 Elsevier Ltd. All rights reserved.

Introduction The majority of the world’s power plants produce electricity with natural gas, coal or oil as primary resource to drive electrical generators. The human induced greenhouse gas emissions particularly carbon dioxide is the main cause of global warming. With availability of economic incentives and other policy mechanisms, there is a growing trend towards wind energy installations all around the world. The electrical power systems throughout the world are undergoing significant changes due to increased penetration level of wind power creating new challenges to system planning and operation. As a result of this an optimal power flow (OPF) is becoming an important tool for power system operators both in planning and operating stage. The main aim of an OPF is to determine the optimal values of control variables of power systems for economic operation under steady state, while satisfying set of equality and inequality constraints. The conventional OPF problem, considering fossil fuel based power plants in system has been extensively studied in the past [1–21]. Moreover, the previous studies associated with the wind–thermal coordination and economics are mostly based on deterministic approach assuming perfect forecast [22–27].

⇑ Corresponding author. Tel.: +91 261 2201664; fax: +91 261 2227334. E-mail address: [email protected] (R. Roy). http://dx.doi.org/10.1016/j.ijepes.2014.07.010 0142-0615/Ó 2014 Elsevier Ltd. All rights reserved.

However, since, in contrast to conventional fossil fuel based power generation, the wind power has random nature, it is necessary to suitably include the wind power uncertainty in optimization problem to determine optimal values of decision variables. Miranda and Hang in [28] proposes a cost model to include wind-powered generators as independent sources, using concepts from the fuzzy set theory. The authors of [28] added a penalty cost in classical economic dispatch problem representing compensation payment to private owners of wind farms for not using the available wind power capacity. In another similar work, an additional cost term for overestimation of the available wind power has been included in dispatch problem along with use probability functions to characterize the wind speed profiles in [29]. To account for wind power uncertainty, a model similar to [28] is analyzed using particle swarm optimization algorithm in [30]. The pioneering work by Hetzer [29], which includes over-estimation and under-estimation of available wind power in classical economic dispatch model, has been studied extensively during recent years in [31–41]. In all these studies the uncertainty of wind power is modeled by Weibull probability distribution function. Although all these studies provides valuable insight into the economic dispatch strategies for wind integrated power systems, none of them except [29] addresses the effect of change in wind speed profile and wind power cost coefficients on optimal dispatch schedule of power plants. Moreover, though the impact of different wind speed profile and cost coefficients corresponding to overestimation and

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

underestimation of available wind power is presented in [29], the reactive power capability of the wind farms, transmission line thermal loading limits and bus voltage constraints are not included in the study. This paper solves an optimal power flow (OPF) problem of wind–thermal of power system for network operating conditions similar to those reported in [42–45] with the wind speed profile characterized by Weibull probability function [46]. In [42] the OPF model, which accounts for the additional costs of managing wind intermittency, is based on probability/relative frequency histograms of forecasting error. However, the study in [42] considers induction generator based wind farm and hence the reactive power capability of variable speed wind turbine is not included in OPF formulation. The reactive power supply capability of variable speed wind turbines e.g. doubly fed induction generator (DFIG) [47] has been utilized for bus voltage improvement in [43–45] with other operating conditions assumed same as given in [42]. In [43,45] wind–thermal optimal power flow (WTOPF) problem is implemented on IEEE 30 bus test system by applying Bacterial Foraging and Modified Bacteria Foraging algorithms respectively. While authors of [44] have implemented wind–thermal OPF on IEEE 39 bus test system by applying self adaptive evolutionary programming technique. Although the work in [42–45] which is based on cost model suggested in [29] appears to be more practical, none of these papers compares the results with published literature to validate the performance of different techniques employed. Moreover, the effect of change in wind speed profile and wind power cost coefficients on optimal dispatch schedule of power plants is not examined in [42–45]. The proposed work assumes that the power network is integrated with wind a farm which consists of variable speed wind turbines capable of supplying reactive power to grid. The unique contributions of this work are given:  The conventional OPF problem is implemented for standard IEEE 30 bus test system by GABC algorithm.  The OPF formulation considering stochastic wind power, reactive power capability of wind turbines and emission constraint is developed and implemented by GABC algorithm.  The OPF framework developed for wind–thermal system is analyzed considering N  1 contingency criteria.  The effect of different wind speed profiles, wind power cost coefficients and emission constraint on optimal dispatch is examined for OPF framework developed for wind–thermal system. The simulation results obtained for conventional OPF problem as well as OPF incorporating stochastic wind power are compared with other methods available in literature. The rest of the paper is organized as follows. Section ‘OPF problem formulation’ gives the mathematical formulation of the OPF problem considering wind power; Section ‘Optimization algorithms’ provides overview of Gbest guided artificial bee colony algorithm and its implementation for wind thermal optimal power flow; Section ‘Numerical results and discussion’ contains simulation results and discussion of different test cases and scenarios. Conclusions are summarized in Section ‘Conclusion’. OPF problem formulation The primary objective of OPF is to optimize the settings of control variables to meet certain objectives while satisfying set of equality and inequality constraints. In general, the OPF problem can be mathematically expressed as follows:

Minimize f ðx; uÞ

ð1Þ

563

subject to,

gðx; uÞ ¼ 0

ð2Þ

hðx; uÞ 6 0

ð3Þ

where f is the objective function to be minimized, g is set of equality constraints representing nodal power injections, and h is set of inequality constraints. The vector u consists of independent variables or control variables and vector x consists of dependent variables or state variables. Control variables The control variables, u of wind–thermal OPF problem are set of variables which are determined by optimization algorithm. These variables are listed below and defined in Eq. (4). (i) pG, real power generation of thermal units at PV buses except slack bus (ii) w, real power generation of wind farms (iii) T, tap settings of transformer (iv) QC, shunt VAR compensation.

uT ¼ ½pG2 . . . pGNG ; w1 . . . wNW ; V G1 . . . V GNG ; Q C 1 . . . Q C NC ; T 1 . . . T NT  ð4Þ where NG, NW, NC and NT are the number of thermal units, the number of wind farms, number of VAR compensators and the number of regulating transformers, respectively. State variables These are the set of variables describing the mathematical state of system. The set of state variables for the OPF problem formulation are given below and defined in Eq. (5). (i) pG1 ; slack bus real power (ii) VL, load (PQ) bus voltage magnitude (iii) QG, generator reactive (iv) Qw, wind farm reactive power (v) Sl, transmission line loading (in MVA)

xT ¼ ½pG1 ; V L1 . . . V LNL ; Q G1 . . . Q GNG ; Q w1 . . . Q wNw ; Sl1 . . . Slnl 

ð5Þ

where NL, and nl indicates number of load buses, and number of transmission lines, respectively. Technical constraints The wind thermal optimal power flow problem which minimizes the objective function given by Eq. (1) should satisfy set of equality and inequality constraints described in the following sections: Equality constraints The equality constraints are typically defined by active and reactive power balance equations at each load as given in Eq. (6).

9 X Pk ¼ 0 = at each node; k ðnet power injectionsÞ X ; Qk ¼ 0

ð6Þ

where Pk and Qk are net active and reactive power injections at kth node. Inequality constraints These constraints consists of active and reactive power outputs of thermal units, active power outputs of wind farms, voltage at generator and load buses, transformer taps, shunt injections, reactive power output limits of wind farms and transmission line loading limits as defined in Eqs. (7)–(14),

564

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

(a) Generator constraints: The real and reactive power of thermal units and wind farms, generator voltage magnitudes should be restricted within their upper and lower bounds:

pGi;min  pGi  pGi;max ;

i ¼ 1; . . . ; NG

ð7Þ

i ¼ 1; . . . NG

ð8Þ

Q Gi;min  Q Gi  Q Gi;max ;

0  wi  wr;i ;



i ¼ 1; . . . NW

ð9Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðSw;i Þ2  ðwr;i Þ2 6 Q w;i 6 ðSw;i Þ2  ðwr;i Þ2 i ¼ 1; . . . NW

V Gi;min  V Gi  V Gi;max ; i ¼ 1; . . . NG ð10Þ where Sw is apparent electric power of wind farm. (b) Transformer constraints: The tap settings of transformer should to be restricted within their specified lower and upper bounds as follows:

T i;min 6 T i 6 T i;max ;

i ¼ 1; . . . NT

i ¼ 1; . . . NC

ð12Þ

(d) Security constraints: The optimal settings obtained while minimizing objective function Eq. (1) should not violet the upper and lower limits of voltage magnitude at load buses. Moreover, the complex power flow through each transmission line should not exceed its capacity limits. These constraints can be mathematically formulated as follows:

V min 6 V Li 6 V max Li Li ; Sli  Smax li ;

i ¼ 1; . . . ; NL

i ¼ 1; . . . ; nl

ð13Þ ð14Þ

where Sl and Smax indicates apparent power flow and maximum apparent power flow limit of transmission line. It can be observed that the control variables are self-constrained while the dependent variables namely load bus voltage magnitude; real power generation output at slack bus, transmission line loading, reactive power generation output of thermal units and wind farms can be included into an objective function as quadratic penalty terms. In these terms, a penalty factor multiplied with the square of the difference between actual value and limiting value of dependent variable is added to the objective function and any unfeasible solution obtained during optimization process is ignored [7]. Mathematically, the augmented objective function can be expressed as follows:

faug



xmax xmin

ð16Þ

It is assumed that the wind farms are owned by independent power producers. The objective of the wind–thermal OPF formulation is to minimize the expected cost of operation (TC) with all system constraints satisfied. The different cost components includes; (i) the fuel cost of coal fired generators, (ii) expected cost of wind power to be purchased, (iii) expected penalty cost for not utilizing available wind power due to network congestion (iv) expected reserve cost due the shortage of wind power (v) emission cost in $/h [34]. The overall cost function of wind thermal power system is given in Eq. (17);

f ðp; w; eÞ ¼

NG NW NW X X X fi ðpGi Þ þ C w;i ðwi Þ þ C Pw;i ðwi;av l  wi Þ i

i

i

NW X þ C rw;i ðwi  wi;av l Þ þ C Tax E

ð17Þ

i

where fi(pGi) represents the fuel cost of ith thermal unit in $/MW, pGi represents scheduled generation of ith thermal unit in MW, wi represents scheduled generation of ith wind farm in MW, wi,av represents actual power generation of ith wind farm in MW, CTax represents carbon tax value in $/ton, E is total emission in ton, P 2 ð NG i ðai þ bi pGi þ ci pGi ÞÞ. The terms a; b and c indicates emission coefficients of thermal units. The first term of objective function (17) is fuel cost function of a thermal generators and it can be expressed in either quadratic form or considering valve loading effects. The second term stands for direct cost of wind power purchased from wind farm operator. The wind farms are assumed to be owned by non-utility operator. Though the wind energy conversion systems have low variable cost, the operation and maintenance cost is substantial. The third term corresponds to penalty cost which is due to underestimation of available wind power that results in wind spillage and hence loss of low-priced available wind power. The fourth term represents reserve cost which is due to overestimation of available wind power. This term becomes significant when the wind power is comparable part of overall power capacity of power network in which the system incur additional costs to provide reliable backup for the wind turbines [48]. The last term corresponds to emission cost or carbon tax for penalizing thermal power producers for minimizing air pollution. Fuel cost curves There are three different types of fuel cost curves namely, quadratic, quadratic with valve point loading effects and piecewise quadratic and defined respectively in Eqs. (18)–(20),

fi ðpGi Þ ¼

NL   2 2 X ¼ f þ kPG PG1  Plim þ kv V Li  V lim G1 Li

x > xmax x < xmin

Objective function

ð11Þ

(c) Shunt VAR compensator constraints: The shunt capacitor VAR compensators should be restricted within their lower and upper bounds as follows:

Q C;i;min 6 Q C;i 6 Q C;i;max ;

xlim ¼

N X ðai þ bi pGi þ ci p2Gi Þ

ð18Þ

i¼1

i¼1 NG  NW  2 2 X X þ kQG Q Gi  Q lim þ kQW Q wi  Q lim Gi wi i¼1

þ ks

Nl  2 X Sli  Smax li

fi ðpGi Þ ¼

N X    ðai þ bi pGi þ ci p2Gi Þ þ jdi sin fi pGi;min  pGi j

ð19Þ

i¼1

i¼1

ð15Þ

i¼1

where kPG ; kv ; kQG ; kQW and ks are penalty factor terms and xlim is limit value of dependant variable x. If x is higher than upper limit, xlim is set to upper limit, likewise if x is less than lower limit, xlim is set to lower limit value. These conditionals are defined in Eq. (16);

fi ðpGi Þ ¼

8 a þ bi1 pGi þ ci1 p2Gi > > > i1 > > > > < ai2 þ bi2 pGi þ ci2 p2Gi > >  > > > > > : a þ b p þ c p2 ih ih Gi ih Gi

9 pGi;min 6 pGi 6 pGi1 > > > > > > = pGi1 6 pGi 6 pGi2 >  pGih1 6 pGi 6 phi;max

> > > > > > > ;

ð20Þ

565

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

where ai, bi, ci, di and fi are the cost coefficients of ith thermal unit, and aik, bik, and cik are the fuel cost coefficients of ith thermal unit for fuel type h [8]. Cost function for wind farm

The terms fv(v) is the fraction of time for which the wind has velocity, m/s. The cumulative distribution function, Fv(v) is the probability that the wind velocity is the same or less than v m/s and is given by the equation, k

F v ðv Þ ¼ 1  eðv =cÞ

ð25Þ

Usually wind farms are owned by private parties and grid operators pay for wind power purchased from wind farm operator based on the power purchase agreement. A cost function for power derived from the wind-farm operator is defined in (21) [29],

The cumulative distribution function is the integral of Weibull pdf.

C w;i ðwi Þ ¼ di wi

The power extracted, w by a wind turbine is presented as a fraction of the upstream wind power in watts as follows [52],

ð21Þ

where di is the direct cost coefficient for ith wind farm. The analysis in literature shows that the total costs of wind generation is around 57% of the total thermal cost [49]. Cost due to underestimation of available wind power An excess of wind generation due to underestimation of wind power may cause system operating problems such as transmission line congestion that leads to wind generator power curtailment during normal operation. Wind farm operators are rewarded depending upon electricity market structure. This amount is generally 15–100% of the amount that the wind farm operator could have earned otherwise [50]. The penalty cost for not using available wind power for ith wind farm is formulated below as:

C pw;i ðwi;av  wi Þ ¼ kp;i ðwi;av l  wi Þ ¼ kp;i

Z

wr;i wi

ðw  wi ÞfW ðwÞdw ð22Þ

where kp,i is penalty cost coefficient due to underestimation of wind power for ith wind farm in $/h, wi,avl is generated or available wind power of ith wind farm in MW and fW(w) is a probability density function of wind power output which is discussed in detail in Sections ‘Wind power characterization based on Weibull probability distribution function’, ‘Mapping wind speed to wind power’, ‘Wind speed probabilities for different wind speeds’ and ‘Analysis of overestimation and underestimation costs’.

The reserve cost is due to shortage of wind power so the system operator has to buy power from other costly reserves to satisfy the load demand. The reserve cost as a result of overestimation of available wind power from wind farm is formulated as:

Z

w¼

1 qAv 3  C p 2

ð26Þ

where q is the air density and A is swept area of wind turbine blade. In this paper, the value power coefficient Cp is assumed to be 0.3. The output of the wind powered generator for a given wind speed input is given by Eqs. (27)–(29):

v < v in

w ¼ 0; w¼

1 qAV 3  C p ; 2

v > vo

ð27Þ

v in  v  v r

ð28Þ

and

vr  v  vo

w ¼ wr

ð29Þ

where vin is cut-in speed, vo is cut out speed and vr is rated speed of wind turbine system in m/s. Wind speed probabilities for different wind speeds The probability of the wind power being zero, rated and in between zero and rated is given by Eqs. (30)–(32) respectively [29].

Prfw ¼ 0g ¼ F v ðv in Þ þ ð1  Fðv o ÞÞ   v in k v o k ¼ 1  exp  þ exp  c c

ð30Þ

  k   k v in vo  exp  Prfw ¼ wr g ¼ F v ðv o Þ  Fðv r Þ ¼ exp  c c

Cost due to overestimation of wind power

    C rw;i wi  wi;av ¼ kr;i wi  wi;av l ¼ kr;i

Mapping wind speed to wind power

ð31Þ "

fw ðwÞ ¼

 k=3  k=3 k 2 1 2w ððk=3Þ1 w Þexp  3ck A ck A

# ð32Þ

wi

ðwi  wÞfW ðwÞdw

Analysis of overestimation and underestimation costs

0

ð23Þ where kr,i is reserve cost coefficient due to overestimation of wind power for ith wind farm in $/h. The value of kr can be as low as the cost coefficient b of thermal units [51] or it can be more than it as it is power taken from generating units having small start up time but may be costly e.g. gas power plant.

C rw ðw1  W 1 Þ ¼ kr ðw1  W 1 Þ ¼ ½s1 þ s2 

Wind power characterization based on Weibull probability distribution function As shown in [46], the characteristic of wind speed variations can be better described by the Weibull probability distribution function, f(v)(v) with shape factor k and scale factor c. During any given time interval, the probability of wind speed being v m/s can be well described by the following equation.

fv ðv Þ ¼

  ðk1Þ k k v eðv=cÞ ; c c

0
The analysis of overestimation cost is carried out by assuming equivalent wind farm of equivalent rated power wr. Let the forecasted wind power of a wind farm is assumed be w1 (0 6 w1 < wr) and power actually produced is assumed to be W1. The overestimation cost is given by Eq. (33),

ð24Þ

where; s1 ¼ w1 PðW 1 ¼ 0Þ    k   k

v in vo þ exp  ¼ w1 1  exp  c c s2 ¼

Z 0

ð33Þ

ð34Þ

w1

ðw1  wÞfW ðwÞdw

ð35Þ

Also, the equations for underestimation can be written as

C pw ðW 1  w1 Þ ¼ kp ðW 1  w1 Þ ¼ s3 þ s4

ð36Þ

566

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

where 





s3 ¼ ðwr  w1 ÞP W 1 ¼ wr ¼ w1

s4 ¼

Z

  k   k

vr vo exp   exp  c c

ð37Þ

wr

ðw  w1 ÞfW ðwÞdw

ð38Þ

w1

In this paper, the OPF problem is solved for two scenarios by applying GABC algorithm. In first scenario, the conventional OPF problem is solved for standard IEEE 30 bus test system by ignoring all terms of objective function given in Eq. (17), except first term. In this scenario, five different cases have been examined with different cost curves given by Eqs. (18)–(20). In second scenario, an OPF problem is solved for modified IEEE 30 bus wind–thermal power system with objective function represented by Eq. (17). For this scenario, the first term of objective function given by Eq. (17) is replaced by Eq. (18) and the term corresponding to emission cost of thermal plants is ignored. Optimization algorithms Artificial bee colony algorithm The colony of artificial bees consists of three types of bees namely employed bees (Ne), onlooker bees (No) and scout bees [53]. All bees that are associated with a particular food source (probable solution) are called employed bees while the bees waiting in the hive for employed bees are called onlooker bees. The scout bees search for new food sources in a random manner. The employed bees carry the loads of nectar from the food source to the hive and share the information about the position of food source with onlooker bees by performing specific type of dance called waggle dance. All onlookers observe the numerous dances performed by employed bees and based on the information received from them the onlookers choose a best food source to follow depending on the probability (Pi) value associated with that food source which is calculated using Eq. (39),

fit Pi ¼ PNe i j fit j

ð39Þ

where i e {1, 2, . . .., SN}, SN represents number of solutions (food source positions) generated initially in random manner, fiti is the fitness value of the ith solution, which depends on quality (fitness) of food source position i, Ne is the number of food sources that are equal to the number of employed bees and j e {1, 2, . . ., D} and D indicates number of decision parameters in optimization problem. Eq. (39) indicates that the food sources with better fitness attract more onlooker bees compared to food sources with inferior quality. This means more onlookers are send to location where the probability of getting good quality food source is more. The onlookers or scout bees become employed bees once they find a food source and starts exploiting it. An employed bee becomes scout bee, once the associated food source is completely exhausted. The scout searches for new food source randomly according to Eq. (40).

xi;j ¼ xmin;j þ randðxmax;j  xmin;j Þ

ð40Þ

where xmin,j and xmax,j represents upper and lower limits of jth decision variable of ith solution and the term rand represents random number in the interval 0 and 1. The complete solution vector for ith food source in the population can be represented by xi = {xi1, xi2, . . ., xiD}. In ABC algorithm, the entire population represents a set of possible solutions, xi which also represents position of food source. The solution search process in ABC algorithm is carried out in four steps namely initialization, employed bee phase, onlooker bee phase and

scout bee phase. During initialization, a set of candidate solutions are randomly generated using (40). Each employed and onlooker bee searches for a new food source vij by modifying the old food source xij using Eq. (41):

v ij ¼ xij þ ;ij ðxij  xkj Þ

ð41Þ

where ;ij is a random number in the interval [1, 1], vij is a new candidate solution (food source), xij is current food source (old food source), xkj is food source in the vicinity of xij, j is randomly chosen index such that (i – k), k e {1, 2, . . .., SN}. The new food source vij is compared with xij i.e., the nectar amount of new food source is compared with old one. The bee remembers new food source if it is better than old one otherwise old food source is retained in the memory i.e. greedy selection process is applied. Each food source in the ABC algorithm has a trail counter (TL) which is set or reset depending upon whether the newly found food source is better than old one or not. If a food source associated with employed and onlooker bees does not improve after a certain number of trail limits, the food source is discarded and subsequently the associated bee becomes a scout bee. The scout bee so generated randomly searches a new food source according to Eq. (40) and again becomes an employed bee. The trail limit is a control parameter of ABC algorithm and it controls generation of scout bees. At the end of each search cycle the trail limit value set by user is compared with trail limit counter. If the value of trail limit counter exceeds the predefined trail limit value, the employed bee becomes a scout bee. The process is repeated for maximum number of cycles (MNC) defined by user. Gbest guided artificial bee colony algorithm (GABC) The basic ABC algorithm does not use the global best solution of population in each generation as in case of particle swarm optimization (PSO) algorithm. The ABC algorithm has good global search capability but poor local search capability, while the PSO has good local search ability but poor global search ability. In order to enhance the exploitation capability of ABC algorithm and to improve the ability of ABC to apply the knowledge of the previous good solutions to find even better candidate solutions, an additional term is introduced in Eq. (41) [54]. The modified Eq. (42) directs the search trajectory towards global optima.

v ij ¼ xij þ ;ij ðxij  xkj Þ þ Cðxj  xij Þ;

k–i

ð42Þ

where C is a random number preferably chosen in the interval [0, 2] and the term, xj represents jth element of the global best solution in current generation. Improved GABC algorithm To further strengthen the exploitation capability of ABC algorithm, equation (39) is modified as in (43) [41],

Pi ¼

0:9  fiti þ 0:1 fitbest

ð43Þ

where fitbest is the fitness value of global best solution i.e. the best solution in current solution vector. Thus the new search process takes advantage of knowledge of the global best solution to modify the search path as in case of PSO. The tuning parameter C needs to be properly chosen so that the new solution is not driven away from desired candidate solution. Implementation of GABC to the wind–thermal OPF problem This part describes the application of the proposed algorithm for solving the wind–thermal OPF problem. To implement the

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

proposed algorithm for wind thermal OPF problem, it is necessary to execute the steps (1)–(15). It is to be noted that the procedure outlined below is also applicable to solve conventional OPF problem in that the terms corresponding to wind power and emission are neglected. Step 1: Input the power system data, limiting values of control variables and parameters of optimization algorithm namely SN, Ne, No, TL, C, MNC. Set cycle number, cycle = 1. Step 2: (Initialization) an initial solution of all control variables is randomly generated in the feasible search space using Eq. (40) and presented below in matrix form:

2 6 6 6 Initial solutionsðPopulationÞ ¼ 6 6 6 4

U1

3

7 U2 7 7 7 .. 7 . 7 5

ð44Þ

U Ne  U i ¼ ui;1 ; ui;2 ; . . . ; ui;D

567

Step 7: Calculate probability value for each solution obtained in Step 4 using Eq. (43). Step 8: Allocate onlookers to further improve solutions obtained in Step 4 depending upon the probability value of each solution (obtained in Step 4). Step 9: (Onlooker bee phase) modified solution vectors obtained in Step 4 are further altered in this phase by following procedure as discussed in Step 4. Step 10: Using these modified values of control variables obtained in Step 9, execute Step 3. Step 11: (Greedy selection) retain the best solution obtained in Step 6 & 10 and discard the inferior one. Step 12: (Scout phase) Check for solution which has not improved for predefined trails. If exists, replace it with a new randomly generated solution using Eq. (40). Step 13: Retain the best solution found so far and increment counter of cycle, cycle = cycle + 1. Step 14: If cycle < maximum number of cycles (MNC), go to Step 4. Step 15: STOP, otherwise and print best results.

ð45Þ

where Ne indicates number of solutions (number of employed bees) and D represents number of control variables.

 U i ¼ pG;i ; wi ; V G;i ; Q C;i ; T i

ð46Þ

Different steps involved to solve the optimal power flow problem by GABC algorithm are presented in the flowchart shown in Fig. 1. The flow chart (2), which is N–R power flow program, is repeated within flowchart (1) for given number cycles.

where

 pGi ¼ pG;i;1 ; pG;i;2 ; . . . pG;i;j . . . ; pG;i;NG ;

j–slack bus

 wi ¼ wi;1 ; wi;2 ; . . . ; wi;NW

ð48Þ

 V Gi ¼ V G;i;1 ; V G;i;2 ; . . . ; V G;i;NG 

Q C;i ¼ Q C;i;1 ; Q C;i;2 ; . . . ; Q C;i;NC

ð47Þ

ð49Þ

 T i ¼ T i;1 ; T i;2 ; . . . ; T i;NT

ð50Þ ð51Þ

Step 3: with these values of control variables run Newton– Raphson power flow program and evaluate augmented objective function given by Eq. (15). The term f of Eq. (15) is replaced by the term f(p, w, e) given by Eq. (17). Step 4: (Employed bee phase) randomly generated solutions in previous step are modified in this step by using Eqs. (52)–(56) which are given below for different control variables.

Numerical results and discussion In this section the optimal power flow problem is implemented using GABC algorithm on standard as well as modified IEEE 30-bus test system for three types of studies. In first case, the simulation studies are carried out on a standard IEEE 30-bus test system by solving conventional OPF problem while in second case, the OPF problem is implemented on wind–thermal power system for given wind speed profile and wind power cost coefficients. The purpose of both these studies is to validate the results obtained by GABC algorithm by comparing them with the results available in the literature. In the third case, the OPF studies are presented to analyze the effect of different wind speed profiles and wind power cost coefficients on optimal cost of operation of wind–thermal power system. In all studies carried out, 100 test runs were performed while solving the OPF problem. Optimal power flow without incorporating wind power

pij ðmÞ ¼ pij þ ;ij ðpij  pkj Þ þ Cðpbest  pij Þ

ð52Þ

wij ðmÞ ¼ wij þ ;ij ðwij  wkj Þ þ Cðwbest  wij Þ

ð53Þ

  V G;ij ¼ V G;ij þ ;ij V G;ij  V G;kj þ CðV G;best  V Gij Þ

ð54Þ

  Q C;ij ðmÞ ¼ Q C;ij þ ;ij Q C;ij  Q C;kj þ CðQ C;best  Q C;ij Þ

ð55Þ

  T ij ðmÞ ¼ T ij þ ;ij T ij  T kj þ CðT best  T ij Þ

ð56Þ

where (k – i), pij(m), wij(m), VG,ij(m), QC,ij(m) and Tij(m) represents respectively the modified value of real power of thermal unit, real power of wind farm, generator bus voltage, shunt compensation VAR and tapping of transformer. Similarly the terms pbest, wbest, VG,best, QC,best and Tbest represents values of control variables corresponding to best solution obtained in previous step. The term C is tuning parameter which is set to 1.5 in this paper [54]. Step 5: Using these modified values of control variables execute Step 3. Step 6: (Greedy selection) retain best solution obtained in Step 3 & 5 and discard inferior one.

In this part of study, the performance of the GABC algorithm is tested on the standard IEEE 30-bus test system without incorporating wind power where the objectives are: minimization of generation fuel cost for different fuel cost curves, voltage profile improvement, voltage stability enhancement and voltage stability enhancement during contingency condition. The test system consists of six thermal units at 1, 2, 5, 8, 11 and 13 nodes inter connected with 41 transmission lines [40]. The network system has four transformers with off-nominal tap ratio at lines 6–9, 6–10, 4–12 and 28–27. In addition, buses 10, 12, 15, 17, 20, 21, 23, 24 and 29 are selected as shunt VAR compensation buses. The total system demand is (2.834 + j1.262) p.u. on base of 100 MVA. The lower limit of voltage magnitude for all busses is set to 0.95 p.u. The upper limit of voltage magnitude is set to 1.05 p.u. for all load busses, while the upper limit of all other generator busses is set to 1.1 p.u. The lower and upper limits of all transformer taps are respectively set to 0.9 and 1.1 p.u. In the past, several mathematical techniques have been applied to solve conventional OPF problem such as gradient based methods, Newton-based method, simplex method, sequential linear

568

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

Start Enter power system data, limits of control variables and algorithm parameters Generate initial solutions randomly using (40) and evaluate each solution using flowchart (2)

Start with solutions from flowchart (1)

Run Newton-Raphson power flow program

cycle=cycle+1

Send employed bees to modify each solution using (42) and evaluate each modified solution using flowchart (2)

Evaluate the objective function using (15)

Employ the greedy selection mechanism between new (modified) and old solution

Objective function value

Evaluate the probability values for the solution using (43)

(2)

Allocate onlookers depending upon probability value and modify each solution using (42). Evaluate each modified solution using flowchart (2) Employ the greedy selection mechanism between solution found by onlooker bee and old solution Memorize best solution (position of food source) Find the abandoned solution and replace it with new solution using (40) NO

cycle >MCN Yes

Display the final solution Stop

(1) Fig. 1. Flowchart of GABC algorithm.

programming, sequential quadratic programming, and interior point methods. However most of these methods suffer due to non-linear and discrete objective functions. In order to overcome the limitations of classical mathematical techniques, a wide variety of non-deterministic and heuristic optimization methods have been reported in literature during last two decades. Some of these techniques are: differential evolution, (DE) [4], biogeography-based optimization, (BBO) [5], multi-agent based differential evolution, (MADE) [6], gravitational-search-algorithm, (GSA) [7], artificial bee colony algorithm (ABC) [8], evolutionary programming (EP) [9], modified differential evolution (MDE) [10], modified shuffled frog leaping algorithm (MSFLA) [11], enhanced genetic algorithm (EGA) [12], improved genetic algorithm (IGA) [13], evolving ant direction differential evolution (EADDE) [14], particle swarm optimization with differentially perturbed velocity hybrid algorithm with adaptive acceleration coefficient (APSODV) [15], genetic evolving ant direction PSODV hybrid algorithm (GEADPSODV) [16], particle swarm optimization (PSO) [17], improved evolutionary programming (IEP) [18], benders decomposition and special ordered set method (BD-SOS) [19], improved bacterial foraging

(IBF) [20], harmony search (HS) [55], Parallel PSO [56], krill herd algorithm (KHA) [57], quasi-oppositional teaching learning based optimization (QOTLBO) [58], hybrid fuzzy particle swarm optimization–Nelder Mead algorithm (FPSO–NM) [59] and modified bacteria foraging algorithm (MBFA) [45]. In the following section, the optimal results obtained by GABC algorithm for conventional OPF problem are compared with above methods reported in literature. Quadratic fuel cost curve (Case 1) To study this case, the objective function Eq. (17) is simplified by ignoring all terms except fuel cost function. The fuel cost function for this case is given by Eq. (18). The minimum and maximum limits of real power generation and their cost coefficients are taken from [10]. The optimum control parameter settings obtained by GABC and other algorithms are presented in Table 1. The fuel cost obtained by proposed method is compared with other methods in the literature. The results of this comparison are given in Table 2. The best cost calculated by GABC for this case is 800.0963 $/h which is better than that of original ABC algorithm [8], and other

569

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578 Table 1 Optimal power flow without considering wind power (Case 1). Control variables

Case 1

Case 2

Case 3

Case 4

Case 5

P1 P2 P5 P8 P11 P13 V1 V2 V5 V8 V11 V13 T11 T12 T13 T14 QC10 QC12 QC15 QC17 QC20 QC21 QC23 QC24 QC29 Fuel cost ($/h) Power loss (MW) Voltage deviations L-index

177.2360 48.7010 21.3530 21.1340 11.8980 12.0000 1.0910 1.0710 1.0400 1.0450 1.0600 1.0490 1.0260 0.9300 0.9690 0.9800 3.0520 2.5110 4.3990 5.0000 4.3800 5.0000 3.1310 5.0000 2.2990 800.0963 8.9340 1.009 –

195.2998 52.0571 15.0000 10.0000 10.0000 12.0000 1.0930 1.0702 1.0260 1.0504 1.0693 1.0623 1.0450 0.9820 1.0170 1.0050 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 931.7450 10.9570 0.4575 –

139.9540 54.9960 24.1000 35.0000 19.9130 16.2520 1.0850 1.0670 1.0400 1.0520 1.0940 1.0620 1.0220 0.9270 0.9770 0.9690 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 647.0300 6.8160 0.8010 –

175.8874 48.8577 21.6348 22.3929 12.4071 12.0000 1.0372 1.0224 1.0056 1.0012 1.0223 1.0038 1.0370 0.9000 0.9560 0.9750 5.0000 0.0000 3.0054 0.0000 5.0000 5.0000 5.0000 5.0000 3.7582 803.5785 9.7800 0.1007 –

176.455178 48.859400 21.762600 21.431600 12.161700 12.000000 1.072200 1.062100 1.056700 1.043800 1.028600 1.049000 1.0125 0.9000 1.0000 0.9750 5.0000 5.0000 3.0001 4.9999 4.9999 3.9999 3.9999 4.9999 2.9999 801.5821 9.2700 0.960916 0.1370

Table 2 Comparison of fuel costs (Case 1).

810 Fuel cost ($/h)

EP [9] MDE [10] IBF [20] MSFLA [11] EGA [12] IGA [13] ABC [8] Parallel PSO [49] BD-SOS [19] EADDE [14] APSODV [15] GEADPSODV [16] GABC

802.620 802.376 802.325 802.287 802.060 800.805 800.660 800.640 800.321 800.2041 800.1159 800.1010 800.0963

similar approaches reported in recent literatures. The convergence graph is shown in Fig. 2. The fuel cost obtained by other methods namely DE [4], BBO [5], MADE [6], GSA [7], (HS) [55], KHA [57], QOTLBO [58] and HFPSO–NM [59] are respectively 799.2891 $/h, 799.1116 $/h, 798.7300 $/h, 798.6751 $/h, 798.800 $/h, 799.0311 $/h, 798.9152 $/h and 794.9545 $/h. Though the fuel cost value for these cases is less than GABC approach, none of these methods satisfies all technical constraints mentioned in Section ‘OPF problem formulation’. In particular, in case of [4], the voltage magnitude at all the load buses is not within specified limit except at bus number 26, and 30. While in [7] the best solution is not true as it results in violation of reactive power limits of generators at bus number 2 and 8. These values are found to be, respectively, 50.8539 MVAr and 113.6966 MVAr (the corresponding limits are 20 and 50 MVAr). Furthermore, there is severe violation of voltage limits at all the load buses except at bus number 4, 7 and 28; and the actual value of slack bus power is found to be 177.8070 MW. The best solution given in [55] suffers from voltage

808

Fuel Cost ($/h)

Algorithm

806 804 802 800 0

5

10

15

20

25

30

35

40

45

50

Number of Cycles Fig. 2. Convergence characteristic of IEEE30-bus system for case 1.

magnitude violations at all the load buses except at bus 7, and also the calculated value of slack bus power is 177.3442 MW (value reported in paper is 176.5883 MW). In case of [57,58] all load bus voltage magnitudes are greater than 1.05 p.u. and hence these optimal results do not satisfy inequality constraints and hence cannot be accepted. The best fuel cost published in [59] is 794.9545 ($/h) but according to generation schedule provided in this paper the fuel cost is found to be 800.9195 ($/h). Moreover, the shunt VAR compensation is more than its limiting value of 5 MVAr. Also, the generator and load bus voltages are not within their upper and lower limits rather all generator bus voltages are greater than maximum limit of 1.1 p.u. while all load bus voltages are greater than their maximum limit of 1.05 p.u. The comparison of GABC with other algorithms in terms of minimum, maximum and average fuel cost is given in Table 3. It is evident from Fig. 3, that the frequency of getting best value of fuel cost less than average fuel cost of 800.152605 $/h by GABC algorithm for 100 trial run is better than other methods.

570

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

Table 3 Comparison table for Case 1.

Table 4 Comparison of fuel costs (Case 2).

Algorithm

Max

Min

Average

Simulation time

Algorithm

Fuel cost ($/h)

GABC HFPSO-NM [59] KHA [57] QOTLBO [58] GSA [7] BBO [5] EP [9] SA IEP [18] DE [4] PSO [17] MDE [10] EADDE [14] EADHDE [14] EGA [12]

800.3420 NA 799.1486 801.1229 799.0284 799.2042 805.6100 NA 802.5810 NA NA 802.4040 800.2784 NA 802.14

800.0963 794.9545 799.0311 798.9152 798.6751 799.1116 802.62 799.45 802.4650 799.2891 800.41 802.376 800.2041 800.1579 802.06

800.1526 NA 799.0756 799.0037 798.9131 799.1985 803.5100 NA 802.5210 NA NA 802.3820 800.2412 NA NA

2.76 NA 2.97 NA 10.7582 11.02 51.4 760 5940.78 NA NA 23.25 3.32 NA 76

ABC [8] IEP [18] PSO [17] MDE [10] GSA [7] MADE [6] BBO [5] GABC

945.4495 932.7642 932.3942 930.7930 929.7240 929.6832 919.7647 931.7450

980 975 970 965

Fuel Cost ($/h)

Quadratic fuel cost curve with valve point loadings (Case 2) This case is similar to case 1 however a sine component is added to the cost curves of the generating units at buses 1 and 2 to simulate the valve point loading effects [15]. Moreover the shunt VAR compensation at buses as in case 1 is not considered in this study. The minimum fuel cost obtained from the proposed approach is 931.7450 $/h is shown in Table 1, which is less than PSO [17] and IEP [18]. Table 4 gives fuel cost comparison with other approaches which shows that the fuel cost by GABC is more than BBO [5], MDE [10], GSA [7], MADE [6], KHA [57], and HFPSO-NM [59]. However in case of [5] there is violation of real power balance constraint as the sum of real power of generating units and real power loss is reported as 294.464 MW and 12.18 MW respectively while the total load is 283.4 MW. Also, the true value of slack bus power obtained from power flow program is 200.1576 MW, and the apparent power in the transmission line between bus numbers 1 and 2 is 136.5234 MVA (loading limit is 130 MVA). In case of [10] the reactive power output of generators at bus numbers 8, 11 and 13 are found to be respectively 68.5531 MVAr, 56.8561 MVAr, and 65.8913 MVAr violating their specified limits. Moreover, the voltage magnitudes at bus numbers 24, 26 and 30 are below 0.95 p.u. and the apparent power in the transmission line between bus number 1 and 2 is found to be 130.2920 MVA. For the best solution reported in [7], the calculated value of slack bus power is 202.3699 MW (upper limit is 200 MW) and the reactive power outputs of the generators at bus number 2, 5, and 11 is found to be respectively 82.7211 MVAr, 119.6009 MVAr, and 11.8231 MVAr, which violates their specified limits. Also, the voltage magnitude at bus number 18, 19, 20, 23, 24, 26, 29, and 30 is below their lower limit of 0.95 p.u., and the apparent power flow through the

960 955 950 945 940 935 930 0

10

20

30

40

50

60

70

80

90

100

Number of Cycles Fig. 4. Convergence characteristic of IEEE30-bus system for case 2.

transmission line between bus number 1 and 2 is 176.0764 MVA. While in case of [6] the voltage at all load buses except bus number 12, 14, 15, 16, 17 and 18 is much less than minimum value of 0.95 p.u. The fuel cost obtained by (APSODV) [15], (EADDE) [14] and (KHA) [57] is 930.6872 $/h, 930.7459 $/h and 925.994 $/h respectively, but in these papers shunt VAR compensation at buses as in Case 1 is considered and therefore these results cannot be considered for comparison. The fuel cost obtained by (HFPSO-NM) [59] is better than GABC approach. But, generator active power output of generator at bus No. 8 and 13 is 4.652 MW and 7.087 MW which is less than minimum limit of power output at these buses. The variation of the total fuel cost with respect to number of cycles is shown in Fig. 4. The minimum, maximum and average values of fuel cost obtained by GABC algorithm are given in Table 5 and compared with other methods.

Fig. 3. Fuel cost obtained by GABC algorithm for 100 different trials.

571

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578 Table 5 Comparison table for Case 2. Algorithm

Max

Min

Average

Simulation time

GABC KHA [57] HFPSO-NM [59] GSA [7] BBO [5] EP [9] IEP [18] MDE [10]

933.3246 929.7234 NA 932.0487291 919.8876 926.68 958.263 954.073

931.7450 925.994 912.0229 929.7240472 919.7647 919.89 953.573 930.793

932.5348 927.0338 NA 930.9246338 919.8389 921.45 956.46 942.501

2.94 3.08 NA 9.8374 11.15 NA 5614.98 41.85

Piecewise quadratic fuel cost curve (Case 3) In this case, the generators of buses 1 and 2 have multiple fuel option. They operate in two different operating regions having two different types of fuel cost characteristics. The fuel cost coefficients of these generators are available in [9]. The fuel cost coefficients of other four generators connected to bus number 5, 8, 11 and 13 have the same values as in case 1. The proposed approach is applied to this case considering no shunt VAR compensation. The results obtained for optimal settings of control variables for this case study are given in Table 1 shows that, the GABC gives best solution while minimizing fuel cost for this OPF problem. The best fuel cost result obtained by GABC technique is compared with other techniques in Table 6. The convergence graph of GABC algorithm for this OPF problem is shown in Fig. 5. From the results in Table 1, it can be seen that the minimum fuel cost is 647.0300$/ h, which is less in comparison with values reported in the literature. The minimum fuel cost obtained by APSODV [15], EADDE [14] and GEADPSODV [16] is 629.4619 $/h, 629.7223 $/h and 629.4801 $/h respectively but in these papers shunt VAR compensation at load buses as in case 1 is considered and therefore these results cannot be considered for comparison. While the minimum cost obtained by MADE [6] and GSA [7] is respectively 646.848066 $/h and 644.9962 $/h but there is severe voltage violation of all load bus voltage limits in case of [6]. On other hand, in case of [7], the voltage magnitude at bus number 12 obtained by NR load flow program is found to be 1.0555 p.u., which is more than the upper limit of 1.05 p.u. The minimum, maximum and average values of fuel cost obtained by GABC algorithm for this case are given in Table 7 and compared with other well established method. Voltage profile improvement (Case 4) The objective function J, chosen in this case is a suitable linear combination of fuel cost and voltage deviation. This objective function balances the two objectives in such a way that one objective should not dominate the other. The new objective function may be expressed as in Eq. (57).

J¼

NPQ N X X ai þ bi pGi þ ci p2Gi þ g jV Lk  1j i¼1

ð57Þ

k¼1

where NPQ is the number of load buses and g is weighting factor which is set 100 in this study. The overall objective function is minimized after satisfying all the constraints given by Eqs. (6)–(14) mentioned in Section ‘Optimization algorithms’. The result obtained from the GABC method has been compared with BBO [5], PSO [17], DE [4] and GSA [7], KHA [57] and HFPSO–NM [59]. It is clear from the result that the total fuel cost has reduced to 803.5785 $/h and voltage deviation has reduced to 0.1007 p.u. as shown in Table 8. Its convergence characteristic is shown in Fig. 6. The voltage deviation 0.093269 obtained by GSA [7] is reported to be less than proposed method but by definition the actual value calculated is found to be 0.105287. Although HFPSO–NM [59] gives better results but, reactive power compensation at bus number 15, 20 and 21

Table 6 Comparison of fuel costs (Case 3). Algorithm

Fuel cost ($/h)

DE [4] IEP [18] ABC [8] MDE [10] EP [9] BBO [5] PSO [17] GSA [7] GABC

650.822 649.312 649.0855 647.846 647.79 647.7437 647.69 646.8481 647.03

exceeds 5 MVAr limits. The shunt power compensator maximum limit is 0–5 MVAr. The voltage performance index obtained by KHA [57] is 0.0996 which is better than GABC but on the other hand fuel cost obtained is 804.6337 $/h which is greater than GABC algorithm. The minimum, maximum and average values of fuel cost obtained by GABC algorithm for this case are given in Table 9 and compared with other methods reported in literature. Voltage stability enhancement with fuel cost minimization (Case 5) The increased loading of power transmission system can cause voltage instability and voltage collapse. The voltage collapse proximity indicator is measure of closeness of the system operating point to voltage collapse. The determination of L-index of the load buses in the system is most popular method of analyzing the voltage collapse proximity indicator. The bus with the highest L-index value is considered to be the most vulnerable bus in the system. The L-index calculation procedure can be found in [60]. The simultaneous optimization of fuel cost and the voltage stability enhancement leads to the following single objective function,

J¼

N X ai þ bi pGi þ ci p2Gi þ gðmaxðLk ÞÞ

ð58Þ

i¼1

where k indicates number of load buses, g is user defined weighting factor which is set to 6000 in this study and Lk is L-index of the Kth bus. The optimal setting of control variables for this case study is given in Table 1 while the comparison of optimal solution obtained by GABC with other methods is given in Table 10. The best result given in [4] is an infeasible solution because the reactive powers of the generators at buses 8, 11, and 13 are 122.6995 MVAr, 19.5035 MVAr, and 33.1927 MVAr, respectively; which violate their limits given in [61]. Also, the voltage magnitudes at all the load buses except at bus 4 violate their upper limits, and the line flow in transmission line between bus 6 and bus 8 is found to be 78.0045 MVA. For the optimum control variables given in [5], the voltage magnitudes at buses 27 and 29 were recalculated and it was found that these values are 1.0664 p.u. and 1.0634 p.u., respectively. These values are higher than their upper limits. In addition, the line flow in transmission line connecting bus 6 and

572

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578 Table 9 Comparison table for Case 4.

675

Fuel Cost ($/h)

670

665

660

655

Algorithm

Max

Min

Average

Simulation time

GABC KHA [57] HFPSO-NM [59] GSA [7] BBO [5] PSO [17] DE [4]

0.1097 NA NA 0.094171 0.1207 NA NA

0.1007 0.0996 0.09764 0.093269 0.102 0.0891 0.1357

0.1052 NA NA 0.093952 0.1105 NA NA

2.98 3.22 NA 11.5873 13.23 NA NA

650

Table 10 Comparisons of the results obtained for Case 5. 645 0

10

20

30

40

50

60

70

80

90

100

Number of Cycles Fig. 5. Convergence characteristic of IEEE30-bus system for case 3.

Table 7 Comparison table for Case 3. *

Algorithm

Max

Min

Average

Simulation time

GABC HFPSO-NM [59] GSA [7] BBO [5] EP [9] IEP [18] DE [4] PSO [17] MDE [10]

647.1234 NA 646.938163 647.7928 652.67 651.125 NA 647.87 650.664

647.03 647.618 646.848066 647.7437 647.79 649.312 650.8224 647.69 647.846

647.0767 NA 646.896273 647.7645 649.67 650.217 NA 647.73 648.356

3.48 NA 10.2716 11.94 51.6 6025.6 NA NA 37.05

Table 8 Comparison of different methods for voltage profile improvement cost (Case 4). Algorithm

Fuel cost ($/h)

Vol. deviation

BBO [5] PSO [17] DE [4] GSA [7] GABC

804.9980 806.3800 805.2620 804.3150 803.5790

0.102 0.0891 0.1357 0.09327 0.1007

880 870

Algorithm

L-index

Fuel cost

ABC DE [4] BBO [5] GSA [7] HS [48] PSO [17] GABC

0.1379 0.1219* 0.1104* 0.1162* 0.1075* 0.1246* 0.137

801.6650 807.5272 805.7252 806.6013 799.9401 801.7600 801.5821

These solutions violate voltage limit constraints.

11 and 13 come 10.9947 MVAr, and 16.0325 MVAr, respectively; which violate their limits given in [60]. In addition, the voltage magnitudes at buses 9, 10, 12 and 27 violate their upper limits. The best solution given in [55] is an infeasible solution because reactive powers of the generators at buses 8, 11 and 13 were given as 70.1633 MVAr, 21.8574 MVAr, and 24.8971 MVAr, respectively; which violate their corresponding limits given in [60]. In addition, the bus voltage magnitudes at all the PQ buses are higher than their upper limits, and the true value of the line flow corresponding to the transmission line connecting bus 6 and bus 8 is 37.3929 MVA (upper limit is 36 MVA). The best result obtained by KHA [57], QOTLBO [58] and HFPSO– NM [59] overwhelms the results obtained by GABC algorithms but in all these methods, the load bus voltage constraints are violated. In case of KHA [57], all load bus voltages are greater than 1.05 p.u., while in case of QLTB [58], except at load bus 7 and 28, all load bus voltages are greater than 1.05 p.u. Moreover, in case of HFPSO–NM [59], except load buses 4, 6, 7 and 28, all load buses voltages are greater than 1.05 p.u. The minimum, maximum and average values of fuel cost obtained by GABC algorithm for this case are given in Table 11 and it is compared with other methods reported in literature. The convergence characteristic is shown in Fig. 7.

Fuel cost ($/h)

860

Optimal power flow considering wind power 850 840 830 820 810 800 0

10

20

30

40

50

60

70

80

90

100

Number of Cycles Fig. 6. Convergence characteristic of IEEE30-bus system for case 4.

bus 8 is 38.1276 MVA. It is found that for the optimal control settings obtained in [7], the voltage magnitude at all load buses violates the upper limit of all load buses by large amount. The optimum control variables given in [17] represent an infeasible solution because the reactive powers of the generators at buses

In this section, GABC algorithm is used to solve OPF problem for power system incorporating stochastic wind power in addition to thermal generators. This study is divided into two parts. In first part, the wind power modeled using Weibull probability distribution function is included in objective function in the form of terms like overestimation cost, underestimation cost and direct cost. The objective function is given by Eq. (17) is minimized subject to all technical constraints given in Section ‘OPF problem formulation’. It is to be noted that the term corresponding to emission in Eq. (17) is not included for this case. The fuel cost curve given by Eq. (18) is assumed to model output of thermal units. In second part, the OPF problem is solved to minimize the objective function defined in first part for different line outage scenarios (N  1 contingency conditions). This study intends to investigate the performance of GABC algorithm in terms of its ability set optimal values of all control variables responsible to maintain the bus voltages close to nominal value.

573

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578 Table 11 Comparison table for Case 5.

910 900

Max

Min

Average

Simulation time

890

GABC KHA [59] HFPSO-NM [57] QOTLBO [58] GSA [7] BBO [5] PSO [17] DE [4]

0.1410 NA NA NA 0.12284 0.1214 NA NA

0.1370 0.1092 0.1023 0.1256 0.116247 0.1104 0.1246 0.1219

0.139 NA NA NA 0.120538 0.1186 NA NA

2.89 3.16 NA NA 13.6378 16.29 NA NA

880

Total cost ($/h)

Algorithm

870 860 850 840 830 820 810 0

5

10

15

870

20

25

30

35

40

45

50

Number of cycles Fig. 8. Convergence characteristics of GABC for Comparative study with [43].

860

Fuel cost ($/h)

850

generation schedule and convergence characteristic for this simulation case is in Table 13 and Fig. 9. From Fig. 9, it can be noted that the number of iterations taken by proposed method is less than 30 while for same network operating condition the number of iterations taken by MBFA [45] and GA [45] to track optimal solution is more than 10,000. This shows that the simulation results obtained by GABC algorithm for wind–thermal power system are better than BFA [43], MBFA [45] and GA [45].

840 830 820 810 800 0

10

20

30

40

50

60

70

80

90

100

Number of Cycles Fig. 7. Convergence characteristic of IEEE30-bus system for case 5.

Optimal cost evaluation under normal condition In this section an OPF problem is solved for wind–thermal system and the optimal results obtained for modified IEEE 30 bus test system [62] by proposed method are compared with BFA [43] and MBFA [45] algorithms. The standard test system consists of six thermal units at bus number 1, 2, 5, 8, 11 and 13. The system is modified by replacing three thermal units at bus number 5, 11 and 13 by DFIG based wind farms with reactive power capability. The optimal value of cost of system operation obtained by GABC was found to be 819.2931 $/h when the test system data given in [43] was considered. The generation schedule and convergence graph for this simulation case is given in Table 12 and Fig. 8 respectively. It was observed that the optimal cost of operation of wind– thermal system obtained by BFA [43] is higher than GABC by 128.20 $/h. The optimal value of cost of system operation obtained by GABC was found to be 734.5377 $/h for scenario 1 mentioned in [45] when the test system data given in [45] was considered. On other hand, the cost of system operation given for BFA [45] and GA [45] algorithms are 2062.52 $/h and 2063.34 $/h respectively. It means that the optimal results obtained by GABC algorithm are much better than MBFA [45] and GA [45]. The optimum

System analysis under (N  1) contingency states To investigate the effectiveness of GABC algorithm in terms of its ability to deal with OPF problem under contingency conditions when system contains stochastic wind power, it was applied to system by creating line outage conditions. The test system and data given in [45] is used to study this simulation case. There are five contingency conditions considered in the analysis presented in this section. For these five conditions the voltage violation index (VVI) is calculated. Table 14 shows different contingency conditions and corresponding (the voltage violation index) VVI for different algorithms. Figs. 10–14 show the variation in voltages under 10–20, 15–18, 15–23, 12–15 and 6–28 line outage conditions. It is evident from Figs. 10–14 that the variation in voltage profile under contingency analysis for GABC algorithm is less than BFA [45] and GA [45]. Moreover, it can be noted that voltage magnitude of all buses obtained under different contingency conditions by applying GABC algorithm are close to the nominal voltage before contingency. This indicates that GABC performs better than BFA [45] and GA [45] under contingency conditions. Effect of different wind speed profiles, wind power cost coefficients and emission constraint on OPF solutions In this section the OPF solutions are obtained for wind thermal power system by applying GABC algorithm however the main focus of this study is to analyze the impact of different variables

Table 12 Comparative study with [43]. Bus no.

WT1 WT2 WT3 TG1 TG2 TG3 Fuel cost

5 11 13 1 2 8 ($/h)

NR – Not reported.

GABC

BFA [43]

Voltage

Power

Reserved real power

Excess power

Voltage

Power

Reserved real power

Excess power

1.0334323 1.0801729 1.0776187 1.0543143 1.0476410 1.0437192 819.2931

60.000000 60.000000 59.999979 50.219319 20.581065 35.000000

25.140793 3.4623820 15.077363 – – –

0.4386995 0.0082242 0.0029253 – – –

1.0092 1.082 1.071 1.043 1.0089 1.06 947.5

50.729 40.405 39.162 56.530 34.285 65.956

25.0195 7.0844 11.5768 – – –

3.1101 11.6107 12.5528 – – –

574

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

Table 13 Comparative study with [45]. GABC Bus no. WT1 5 WT2 11 WT3 13 TG1 1 TG2 2 TG3 8 Fuel cost ($/h)

BFA [45]

GA [45]

Voltage

Power

Reserved real power

Excess power

Voltage

Power

Voltage

Power

1.033642 1.071041 1.080754 1.059073 1.052726 1.035523 734.5377

50 40 40 69.351162 62.843432 24.641747

20.156731 2.3064850 10.050941 – – –

0.355943 0.004178 0.001490 – – –

– – – – – – 2062.52

48.80 40.00 37.80 59.26 61.70 38.00

– – – – – – 2063.34

50.00 40.00 38.00 59.46 62.00 40.00

860 Nominal voltage (BFA)[45] Voltage during contingency (BFA)[45]

840

Nominal voltage (GABC)

1.1

Voltage (p.u.)

Total cost ($/h)

Voltage during contingency (GABC)

820 800 780 760 740

1.05 1 0.95 0.9

720 0

5

10

15

20

25

30

35

40

45

50

0.85 0

Number of cycles

5

10

15

Fig. 9. Convergence characteristics of GABC for Comparative study with [45].

Voltage violation index (p.u.)

10–20 15–18 15–23 12–15 6–28

GABC

BFA [45]

GA [45]

0.277672 0.191803 0.211411 0.134684 0.245614

0.4994 0.5912 0.5683 0.3319 0.4007

0.5287 0.6372 0.6155 0.3655 0.4216

Voltage (p.u.)

1 2 3 4 5

30

Nominal voltage (BFA)[45] Voltage during contingency (BFA)[45] Nominal voltage (GABC) Voltage during contingency (GABC)

1.1

Line outage no.

25

Fig. 11. Voltage profile comparison during 15–18 line outage.

Table 14 Contingency analysis for GABC algorithm. Sr. no.

20

Bus number

1.05 1 0.95 0.9 0.85 0

10

15

20

25

30

Fig. 12. Voltage profile comparison during 15–23 line outage.

1.05 Nominal voltage (BFA)[45] Voltage during contingency (BFA)[45] Nominal voltage (GABC) Voltage during contingency (GABC)

1 1.1 0.95 0.9 0.85 0

5

10

15

20

25

30

Bus number

Voltage (p.u.)

Voltage (p.u.)

1.1

5

Bus number

Nominal voltage (BFA)[45] Voltage during contingency (BFA)[45] Nominal voltage (GABC) Voltage during contingency (GABC)

1.05 1 0.95

Fig. 10. Voltage profile comparison during 10–20 line outage.

0.9

associated with OPF formulation on optimal solutions. An IEEE 30bus system discussed in Section ‘Optimal power flow without incorporating wind power’ is used as a test system to investigate the effect of different variables such as wind speed scale parameter ‘‘c’’, wind power cost coefficients (kp and kr) and emission constraint on OPF solutions. Two extra wind farms are placed on two separate buses namely bus-26 and 30. These buses are selected based on their voltage stability limits [40]. The voltage

0.85 0

5

10

15

20

25

30

Bus number Fig. 13. Voltage profile comparison during 12–15 line outage.

stability of these buses is highest; hence the wind turbines are placed at these buses. It is assumed that each wind farm represents the aggregate generation model and consists of several pitch

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578 1.1 Nominal voltage (BFA)[45] Voltage during contingency (BFA)[45] Nominal voltage (GABC) Voltage during contingency (GABC)

Voltage (p.u.)

1.05

1

0.95

0.9

0.85 0

5

10

15

20

25

30

Bus number Fig. 14. Voltage profile comparison during 6–28 line outage.

controlled DFIG based wind turbines capable of supplying active as well as reactive power to the grid. Since there are many variables in the model, and to arrive at some meaningful conclusion, the analysis below considers the effect of the variation of one variable on different cost components of objective function and the penetration level of wind power. In both, the text and figures that follow, the abbreviations UE, OE, FC and TC indicates respectively underestimation cost, overestimation cost, fuel cost and total cost of operation all in $/h. The notations W26 and W30 represent active power outputs from wind farms at node 26 and 30 respectively. In this study, the following parameters of wind turbine are considered: electrical power wr = 2 MW; apparent electric power S = 2.2 MVA; air density q = 1.225 kg/m3; cut-in, rated and cut-out wind speed: 4 m/s, 12 m/s and 25 m/s, respectively. Optimal power flow as a function of the Weibull c factor The cost coefficients associated with wind power due to underestimation and overestimation of wind power are as follows, d1 = d2 = 1.3, kp1 = kp2 = 1, kr1 = kr2 = 4. Notice that the direct cost of wind power is less than average cost of thermal power while the penalty cost for not using available wind power is assumed to be less than direct power [51]. The wind farm capacity at each node is assumed to be 30 MW. The wind speed at given site follows Weibull distribution with shape factor k = 2. The scale factor c is varied from 3 to 60 to assess the effect of increased penetration of wind power on TC, FC, UE, OE and MVA loading of T25–24. From Fig. 15, it can be observed that, for larger value of scale factor c, a larger amount of higher wind speeds will be likely and hence power output from both wind farms increases rapidly. However

575

the amount of wind power injected into the network by W26 is less than W30 due to thermal loading limit of transmission line T25–24. The wind power curtailment due to transmission line bottleneck results in increase in overestimation cost and hence overall operating cost (TC) and fuel cost (FC), even if the wind speed is favorable. The power injected by W30 is almost its rated value in range of value of c from 6 to 29 which is due to availability of two transmission lines, T30–27 and T30–29 to carry power into network. Notice that the total operational cost reduces as scale factor c increases due to low cost of wind power for certain values of scale factor, after certain point however the imbalance cost increases and hence total operational cost also increases. Moreover, for very high value of Weibull parameter the probability of favorable wind speed drops and hence the wind power injection tends to very small value. This in turn leads to small wind power imbalance cost. The right hand side vertical axis of Fig. 15 represents line loading in MVA, wind power outputs from W26, W30 and total wind power W in MW. The convergence characteristics for different values of ‘c’ is given in Fig. 16. While the optimal cost of operation for 100 trail runs is shown in Fig. 17. From Fig. 16 it can be seen that the number of iterations required to track optimal solution is more in case of wind–thermal system as compared to conventional OPF. This is due to the fact that the objective function for wind–thermal system includes additional cost functions to account for direct cost of wind power and imbalance costs, Moreover, the average cost of operation obtained for this case after 100 trial run is found to be 752.7088 $/h. According to Fig. 17, the proposed method has ability to consistently track the optimal solution since the optimal results obtained for many trail runs are less than average cost of operation. Optimal power flow as a function of the reserve cost coefficient kr The reserve cost, that is, the cost of providing spinning reserve capacity by the coal-fired power plants, to account for the forecast error in wind speed, depends on reserve cost coefficient. The following cost coefficients are assumed for both wind farms, d1 = d2 = 1.3, kp1 = kp2 = 1, k = 2 and c = 10. The wind farm capacity at each node is assumed to be 20 MW. As the reserve cost coefficient increases from its base value of 4, the output of both wind farms W26 and W30 begin to fall as shown in Fig. 18. This is as expected from the model, in which, increasing the reserve cost coefficient compels the network operator to be more conservative in planning wind power as a greater cost will be paid for overestimating the amount of wind power to be scheduled in the given time interval. The reduction in low cost wind power results in increase in fuel cost (FC) as well as total operational cost (TC)

Fig. 15. Overestimation cost (OE), underestimation cost (UE) and wind power penetration vs. Weibull ‘c’ parameter.

576

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578 900 Without wind source

880

With wind source C=4 C=8 C=12 C=24 C=32

Total cost ($/h)

860 840 820 800 780 760 740 720

0

10

20

30

40

50

60

70

80

90

100

Number of cycles Fig. 16. Convergence characteristics of IEEE 30 bus test system with additional two wind farms for different values of ‘c’.

753.6

average value

Total cost ($/h)

753.4 753.2 753 752.8 752.6

respectively. The wind farm capacities are same as above. According the cost model given in Eq. (17), the increase in the penalty cost coefficient implies that the system operator should schedule more amount of wind power for minimizing wind spillage and hence to reduce operating cost. The increase in penalty cost can be viewed as policy to encourage clean energy to minimize emission. Fig. 19 shows that the wind power from W26 and W30 increases towards their rated value of 20 MW as the penalty cost coefficient increases. However the wind power curtailment at node 26 becomes indispensable due to line loading limit of T25–24 line. The vertical axis on right hand side of Fig. 19 indicates the line loading limit of transmission line between bus number 24 and 25 (T25–24) in MVA while the left hand side of vertical axis of Fig. 19 indicates the wind power output of wind farms, W26 and W30, and total wind power, W in MW. As shown in Fig. 20, the wind spillage at node 26 causes increase in overestimation cost by large amount compared to underestimation cost as the wind speed forecast is favorable due to larger value of scale factor (c = 10). Fig. 21 shows that the fuel cost reduces with increases in kp as more wind power is injected at lower cost to meet the load demand, however overall operation cost (TC) goes high due to increase in overestimation cost. The vertical axes on left hand and right hand side of Fig. 21, indicates, respectively the total cost of operation, TC ($/h) and fuel cost, FC ($/h).

752.4 752.2 752

Number of trials Fig. 17. Overall cost of operation (TC) obtained by GABC algorithm for wind– thermal system for 100 different trial.

but the power imbalance costs (OE + UC) falls rapidly as wind power penetration tends to decrease for large value of reserve cost. For very high values of reserve cost coefficient the wind power imbalance cost becomes constant and depends on penalty cost coefficient kp. The vertical axis on right hand side of Fig. 18 indicates wind power output from W26, W30 and total wind power injection W in MW. Optimal power flow as a function of the penalty cost coefficient kp The cost coefficients for wind power, overestimation of wind power for both wind farms are as follows, d1 = d2 = 1.3, kr1 = kr2 = 4. The shape factor and scale factor are set to k = 2 and c = 10

Optimal power flow with emission constraint To address environment issues carbon tax has been implemented in many countries as a policy tool to reduce carbon emissions. In this simulation study an additional emission constraint is included in cost function given by Eq. (17) and the effect of wind power on emission pollutants is examined [34]. The value of carbon tax is assumed to be $ 20. The other parameters are the same as in Section ‘Optimal power flow as a function of the Weibull c factor’. From Fig. 22 it is clear that as the contribution of wind power increases in the system, the overall cost (TC) and emission (E) drops however for larger value of scale parameter (c) the value of TC and E increases. This is mainly due to increase in overestimation cost for larger values of scale parameter (c) as shown in Fig. 23. The main reason for rise in overestimation cost is the due to line loading limit of T25–24 line which results in wind power curtailment at bus number 26. The values of wind power at bus number 26 for different values of scale factor are indicated on the right-hand side vertical axis of Fig. 23.

Fig. 18. Wind power penetration and different costs vs. reserve cost coefficient, kr.

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

577

Fig. 22. Fuel cost (FC) and emission (E) vs. Weibull ‘c’ parameter.

Fig. 19. Wind power penetration and MVA line loading of T24–25 vs. penalty cost coefficient, kp.

Fig. 23. Overestimation cost (OE), underestimation cost (UE) and wind power W26 vs. Weibull ‘c’ parameter.

Conclusion

Fig. 20. Overestimation cost (OE) and underestimation cost (UE) vs. penalty cost coefficient, kp.

This paper presents optimal power flow study considering probabilistic nature of wind power. The wind power intermittency is modeled by two parameter weibull probability function. The optimization problems are solved by recently developed Gbest artificial bee colony algorithm. The simulations results obtained for standard IEEE-30 bus system without considering wind power are compared with other methods available in the literature. Moreover, the optimal scheduling of wind turbine systems and thermal units is compared with other methods for similar OPF framework. From simulation results, it is found that the proposed method based on GABC algorithm gives better accuracy of results compared to other well established methods tried in the past with faster convergence and better solution quality. Finally the optimal power flow study considering stochastic wind power shows that the optimal value of wind power from a particular wind farm not only depends on the values of the reserve and penalty cost coefficients associated with the wind farms but also on its location in power network as the transmission system capacity can limit wind power injection. References

Fig. 21. Overall cost of operation (TC) and fuel cost (FC) vs. penalty cost coefficient, kp.

[1] Carpentier J. Optimal power flows. Int J Electr Power Energy Syst 1979;1:3–15. [2] Momoh JA, El-Hawary ME, Adapa R. A review of selected optimal power flow literature to 1993 part i: nonlinear and quadratic Programming Approaches. IEEE Trans Power Syst 1999;14:96–103.

578

R. Roy, H.T. Jadhav / Electrical Power and Energy Systems 64 (2015) 562–578

[3] Monoh JA, Ei-Hawary ME, Adapa R. A review of selected optimal power flow literature to 1993 part ii: newton, linear programming and Interior Point Methods. IEEE Trans Power Syst 1999;14:105–11. [4] Abou El Ela AA, Abido MA, Spea SR. Optimal power flow using differential evolution algorithm. Electr Power Syst Res 2010;80:878–85. [5] Bhattacharya A, Chattopadhyay PK. Application of biogeography-based optimisation to solve different optimal power flow problems. IET Gener Transm Distrib 2011;5:70–80. [6] Sivasubramani S, Swarup KS. Multiagent based differential evolution approach to optimal power flow. Appl Soft Comput 2012;12:735–40. [7] Duman S, Güvenç U, Sönmez Y, Yörükeren N. Optimal power flow using gravitational search algorithm. Energy Convers Manag 2012;59:86–95. [8] Rezaei Adaryani M, Kkarami A. Artificial bee colony algorithm for solving multi-objective optimal power flow problem. Int J Electr Power Energy Syst 2013;53:219–30. [9] Yuryevich J, Wong KP. Evolutionary programming based optimal power flow algorithm. IEEE Trans Power Syst 1999;14:1245–50. [10] Sayah S, Zehar K. Modified differential evolution algorithm for optimal power flow with non-smooth cost functions. Energy Convers Manag 2008;49:3036–42. [11] Niknam T, Narimani M Rasoul, Jabbari M, Malekpour AR. A modified shuffle frog leaping algorithm for multi-objective optimal power flow. Energy 2011;36:6420–32. [12] Bakirtzis AG, Biskas PN, Zoumas CE, Petridis V. Optimal power flow by enhanced genetic algorithm. IEEE Trans Power Syst 2002;17:229–36. [13] Lai LL, Ma JT, Yokoyama R, Zhao M. Improved genetic algorithms for optimal power flow under both normal and contingent operation states. Int J Electr Power Energy Syst 1997;19:287–92. [14] Vaisakh K, Srinivas LR. Evolving ant direction differential evolution for OPF with non-smooth cost functions. Eng Appl Artif Intell 2011;24:426–36. [15] Vaisakh K, Srinivas LR. Adaptive PSODV algorithm for OPF with non-smooth cost functions and statistical analysis. Simul Model Pract Theory 2011;19:1824–46. [16] Vaisakh K, Srinivas LR, Meah K. Genetic evolving ant direction PSODV hybrid algorithm for OPF with non-smooth cost functions. Electr Eng 2013;95:185–99. [17] Abido MA. Optimal power flow using particle swarm optimization. Int J Electr Power Energy Syst 2002;24:563–71. [18] Ongsakul W, Tantimaporn T. Optimal power flow by improved evolutionary programming. Electr Power Components Syst 2006;34:79–95. [19] Amjady N, Ansari MR. Non-convex security constrained optimal power flow by a new solution method composed of Benders decomposition and special ordered sets. Int Trans Electr Energy Syst 2013. [20] Amjady N, Fatemi H, Zareipour H. Solution of optimal power flow subject to security constraints by a new improved bacterial foraging method. IEEE Trans Power Syst 2012;27:1311–23. [21] Kumar S, Chaturvedi DK. Optimal power flow solution using fuzzy evolutionary and swarm optimization. Int J Electr Power Energy Syst 2013;47:416–23. [22] Denny E, O’Malley M. Wind generation, power system operation, and emissions reduction. IEEE Trans Power Syst 2006;21:341–7. [23] Arabali A, Ghofrani M, Etezadi-Amoli M. Cost analysis of a power system using probabilistic optimal power flow with energy storage integration and wind generation. Int J Electr Power Energy Syst 2013;53:832–41. [24] Lee JC, Lin WM, Liao GC, Tsao TP. Quantum genetic algorithm for dynamic economic dispatch with valve-point effects and including wind power system. Int J Electr Power Energy Syst 2011;33:189–97. [25] Chang YC, Lee TY, Chen CL, Jan RM. Optimal power flow of a wind-thermal generation system. Int J Electr Power Energy Syst 2014;55:312–20. [26] Yuan Y, Zhang X, Ju P, Qian K, Fu Z. Applications of battery energy storage system for wind power dispatchability purpose. Electr Power Syst Res 2012;93:54–60. [27] Guo CX, Bai YH, Zheng X, Zhan JP, Wu QH. Optimal generation dispatch with renewable energy embedded using multiple objectives. Int J Electr Power Energy Syst 2012;42:440–7. [28] Miranda V, Hang PS. Economic dispatch model with fuzzy wind constraints and attitudes of dispatchers. IEEE Trans Power Syst 2005;20:2143–5. [29] Hetzer J, Yu DC, Bhattarai K. An economic dispatch model incorporating wind power. Energy Convers, IEEE Trans 2008;23:603–11. [30] Wang L, Singh C. Balancing risk and cost in fuzzy economic dispatch including wind power penetration based on particle swarm optimization. Electr Power Syst Res 2008;78:1361–8. [31] Liu X, Xu W. Minimum emission dispatch constrained by stochastic wind power availability and cost. IEEE Trans Power Syst 2010;25:1705–13. [32] Liu X, Member S, Xu W. Economic load dispatch constrained by wind power availability: a here-and-now approach. IEEE Trans Sustain Energy 2010;1:2–9. [33] Liu X. Economic load dispatch constrained by wind power availability: a waitand-see approach. IEEE Trans Smart Grid 2010;1:347–55.

[34] Yao F, Dong ZY, Meng K, Xu Z, Iu HH, Wong KP. Quantum-inspired particle swarm optimization for power system operations considering wind power uncertainty and carbon tax in Australia. IEEE Trans Ind Informatics 2012;8:880–8. [35] Zhao J, Wen F, Dong ZY, Xue Y, Wong KP. Optimal dispatch of electric vehicles and wind power using enhanced particle swarm optimization. IEEE Trans Ind Informatics 2012;8:889–99. [36] Reddy SS, Panigrahi BK, Kundu R, Mukherjee R, Debchoudhury S. Energy and spinning reserve scheduling for a wind-thermal power system using CMA-ES with mean learning technique. Int J Electr Power Energy Syst 2013;53:113–22. [37] Morshed MJ, Asgharpour A. Hybrid imperialist competitive-sequential quadratic programming (HIC-SQP) algorithm for solving economic load dispatch with incorporating stochastic wind power: a comparative study on heuristic optimization techniques. Energy Convers Manage 2014;84:30–40. [38] Zhang Y, Yao F, Iu HH-C, Fernando T, Wong KP. Sequential quadratic programming particle swarm optimization for wind power system operations considering emissions. J Mod Power Syst Clean Energy 2013;1:231–40. [39] Azizipanah-Abarghooee R, Niknam T, Roosta A, Malekpour AR, Zare M. Probabilistic multiobjective wind-thermal economic emission dispatch based on point estimated method. Energy 2012;37:322–35. [40] Mondal S, Bhattacharya A, Nee Dey SH. Multi-objective economic emission load dispatch solution using gravitational search algorithm and considering wind power penetration. Int J Electr Power Energy Syst 2013;44:282–92. [41] Jadhav HT, Roy R. Gbest guided artificial bee colony algorithm for environmental/economic dispatch considering wind power. Expert Syst Appl 2013;40:6385–99. [42] Jabr RA, Pal BC. Intermittent wind generation in optimal power flow dispatching. IET Gener Transm Distrib 2009;3:66. [43] Mishra S, Mishra Y, Vignesh S. Security constrained economic dispatch considering wind energy conversion systems. IEEE Power Energy Soc Gen Meet; 2011. [44] Shi L, Wang C, Yao L, Ni Y, Bazargan M. Optimal power flow solution incorporating wind power. IEEE Syst J 2012;6:233–41. [45] Panda A, Tripathy M. Optimal power flow solution of wind integrated power system using modified bacteria foraging algorithm. Int J Electr Power Energy Syst 2014;54:306–14. [46] Patel, Mukund R. Wind and solar power systems: design, analysis, and operation. CRC Press; 2012. [47] Jadhav HT, Roy R. A comprehensive review on the grid integration of doubly fed induction generator. Int J Electr Power Energy Syst 2013;49:8–18. [48] Alhasawi FB, Milanovic JV. Techno-economic contribution of FACTS devices to the operation of power systems with high level of wind power integration. IEEE Trans Power Syst 2012;27:1414–21. [49] Chade Ricosti JF, Sauer IL. An assessment of wind power prospects in the Brazilian hydrothermal system. Renew Sustain Energy Rev 2013;19: 742–53. [50] Aparicio N, MacGill I, Rivier Abbad J, Beltran H. Comparison of wind energy support policy and electricity market design in Europe, the United States, and Australia. IEEE Trans Sustain Energy 2012;3:809–18. [51] Morales JM, Conejo AJ, Liu K, Zhong J. Pricing electricity in pools with wind producers. IEEE Trans Power Syst 2012;27:1366–76. [52] Slootweg JG, De Haan SWH, Polinder H, Kling WL. General model for representing variable speed wind turbines in power system dynamics simulations. IEEE Trans Power Syst 2003;18:144–51. [53] Karaboga D, Basturk B. A powerful and efficient algorithm for numerical function optimization: Artificial bee colony (ABC) algorithm. J Glob Optim 2007;39:459–71. [54] Zhu G, Kwong S. Gbest-guided artificial bee colony algorithm for numerical function optimization. Appl Math Comput 2010;217:3166–73. [55] Sivasubramani S, Swarup KS. Multi-objective harmony search algorithm for optimal power flow problem. Int J Electr Power Energy Syst 2011;33: 745–52. [56] Kim J-Y, Mun K-J, Kim H-S, Park JH. Optimal power system operation using parallel processing system and PSO algorithm. Int J Electr Power Energy Syst 2011;33:1457–61. [57] Roy PK, Paul C. Optimal power flow using krill herd algorithm. Int Trans Electr Energy Syst 2014. [58] Mandal B, Kumar Roy P. Multi-objective optimal power flow using quasioppositional teaching learning based optimization. Appl Soft Comput 2014;21:590–606. [59] Joorabian M, Afzalan E. Optimal power flow under both normal and contingent operation conditions using the hybrid fuzzy particle swarm optimisation and Nelder-Mead algorithm (HFPSO-NM). Appl Soft Comput J 2014;14:623–33. [60] Kessel P, Glavitsch H. Estimating the voltage stability of a power system. IEEE Trans Power Deliv 1986;1:346–54. [61] http://shodhganga.inflibnet.ac.in/bitstream/10603/1221/18/18_appendix.pdf. [62] Pai MA. Computer methods in power system analysis. TMH Publishers; 1996.