Stochastic mid-term generation scheduling incorporated with wind power

Electrical Power and Energy Systems 64 (2015) 937–946

Contents lists available at ScienceDirect

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

Stochastic mid-term generation scheduling incorporated with wind power Taher Niknam a,⇑, Hamid Reza Massrur b a b

Department of Electrical and Electronics Engineering, Shiraz University of Technology, Shiraz, Iran Department of Electrical Engineering, Shiraz Branch, Islamic Azad University, Shiraz, Iran

a r t i c l e

i n f o

Article history: Received 6 September 2013 Received in revised form 7 July 2014 Accepted 28 July 2014

Keywords: Mid-term generation scheduling Point estimate methods Wind power generation Modified gravitational search algorithm Self-adaptive strategy

a b s t r a c t A challenge now facing system operator is how to schedule optimally the generation units in a wind integrated power system over a one year time horizon considering the effects of wind forecasting and variability; also, regarding the effects of load uncertainty. By the same token, this paper first develops a new formulation for Stochastic Mid-term Generation Scheduling (SMGS). In the formulation, 2m + 1 point estimate method is developed to accurately estimate the output variables of Mid-term Generation Scheduling (MGS) problem. Then, the formulation is combined with adaptive modified gravitational search algorithm and a novel self-adaptive wavelet mutation strategy for the establishment of new robust algorithm for the present problem. It is noteworthy to say that the classical methods considered certain wind information in the deterministic solution of the MGS problem which is not the realistic approach. However, this study improves modeling of wind–thermal system in the MGS problem by considering possible uncertainties when scheduling the generators of power system. The proposed model is capable of taking uncertainty of load and wind into account. The proposed method is applied on two test cases and the numerical results confirmed the efficiency and stability of the proposed algorithm. Ó 2014 Elsevier Ltd. All rights reserved.

Introduction In order to supply a high quality electrical power for the consumer in an economic and secure manner, the System Operator (SO) in the control of the power systems, encounters numerous technical and economic difficulties. Economic operation deals with the Generation Scheduling (GS) in the electrical power system to produce coordination between generation units. The GS problem, from the time perspective, is divided into three short-term, midterm and long-term stages [1]. The MGS problem for a power system is utilized to specify the schedules of generation units which minimize the total generation and operation costs over a scheduling horizon while satisfying the load demand, spinning reserve and operating constraints of generation units and transmission lines constraint. The MGS problem represents many activities such as maintenance scheduling, fuel budgeting, rate forecasting, emission and production costs. Typically, the scheduling horizon of MGS is one year and the decision stage is a week or a month [2]. In each rudimentary time interval of MGS problem all the variables of the problem (e.g. load, thermal generation, etc.) are considered constant with respect to time. ⇑ Corresponding author. Tel.: +98 711 7264121; fax: +98 711 7353502. E-mail address: [email protected] (T. Niknam). http://dx.doi.org/10.1016/j.ijepes.2014.07.076 0142-0615/Ó 2014 Elsevier Ltd. All rights reserved.

Within the last two decades, various optimization models for short-term generation scheduling have been competently advanced [3–8]. But, less consideration has been specifically allocated to mid-term generation scheduling. The MGS is at least as significant as short-term GS because the SO caries out some plans such as maintenance scheduling, budgeting, fuel planning and unit commitment planning through the MGS. Furthermore, the shortterm generation scheduling is based upon MGS outcomes. Ref. [9] suggested a stochastic generation scheduling model to treat the mid-term fuel planning issue. A model for the solution of stochastic mid-term scheduling of a price-maker hydro producer is offered in [10]. It should be emphasized that despite the works carried out in area of medium/long-term generation scheduling to date, some issues still remain to be solved pleasantly. One of these problems is the produced challenge from the combination of wind power systems and power system. A large-scale integrated system offers many issues to the midterm operation and scheduling of the electrical power system because the wind power is difficult to predict and highly intermittent. Specially, in a time horizon of one year, uncertainties of various sorts emerge in a wind–thermal system which cause changing in the power generation of units. Unforeseen increase or decline in power production of a wind farm system may need start units to cover the overflow or shortage unscheduled power produc-

938

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

Nomenclature Indices b c d FE g i, j k p t u w

dimension index uncertain output variable uncertain input variable function evaluation index thermal unit index mass index generation index concentration location index time interval index (month) moment index wind farm index

Constants FEmax the maximum number of function evaluations Klg sensitivity coefficient for flow of line l with respect to power output of thermal unit g (see Appendix A for further information) Klw sensitivity coefficient for flow of line l with respect to power output of wind farm w (see Appendix A for further information) NB number of system buses NBg number of generation buses Ng number of thermal units Nh number of hours at time study Nl number of transmission lines Nt number of hours at time t Nw number of wind farms OMFCT(g) operation and maintenance fixed cost of thermal unit g ($/MW yr) OMFCW(w) operation and maintenance fixed cost of wind farm w ($/MW yr) OMVCT(g) operation and maintenance variable cost of thermal unit g ($/MW h) OMVCW(w) operation and maintenance variable cost of wind farm w ($/MW h) randi, randj random function generator in the range [0, 1] Prat the rated output wind power PGmax,g upper limit of thermal unit g (MW) PGmin,g lower limit of thermal unit g (MW)

tion [11]. Load and wind power uncertainty can cause the weak performance of MGS in terms of level of reliability in supplying loads and fuel consumption. If uncertain data about the wind power and load demand are employed in deterministic conventional simulation the power generation performance cannot be precisely estimated. The generation schedule based upon deterministic cost function yields the lowest expected cost. However, this cost is linked with a relatively large contradiction with reality that can establish a risk. Therefore, for mid-term operation of power system, it becomes essential to consider the stochastic nature of the wind power and the load. Then, a stochastic representation of these must be used in the optimization model to prevent the serious risk for the power generation of units in the MGS problem. In this regard, the operators and planners who have to carry out such a critical task need an appropriate optimization tool. This tool is the uncertainty analysis with stochastic models. The aim of uncertainty analysis is to recognize the various fountainheads of uncertainty and to quantify their impacts on the final model outcomes [12]. In the technical literature, there are a number of techniques to cope with problems under uncertainty. These techniques may be

PLmax,l PLmin,l PWmax,w T

vin vout vrat lzd rzd

maximum power flow limit of transmission line l (MW) minimum power flow limit of transmission line l (MW) maximum generation of wind farm w (MW) number of periods under study the cut-in wind speed the cut-out wind speed the rated wind speed mean value of uncertain input variable zd variance value of uncertain input variable zd

Variables ad,p vector of uncertain input variables for pth concentration location of dth random variable F bij ðkÞ the force applying between the agents i and j in the bth dimension at generation k G(k) the gravitational constant at generation k J() total cost function of MGS problem mi(k) the inertial mass at generation k Mi(k) the active gravitational mass at generation k NFE number of function evaluation of a success run to reach the pre-determined value P expected value of system load demand at time t (MW) load ðtÞ PGD g ðtÞ expected value of load portion of thermal unit g at time t (MW) PGR g ðtÞ expected value of reserve portion of thermal unit g at time t (MW) PW  expected value of generation of wind farm w at time t w ðtÞ (MW) u(g, t) decision variables of unit g at time t Ug(t) commitment state of unit g at time t (on = 1, off = 0) Vw(t) commitment state of wind farm w at time t (on = 1, off = 0) v elbi ðkÞ the velocity of the agent i, in the bth dimension and at generation k W av ;w ðtÞ expected value of maximum available wind power of wind farm w at time t (MW) x(g, t) state variables of unit g at time t Superscript  expected value of random variables

classified into three principal groups [13]: Monte Carlo simulation, analytical methods and approximate methods. The main drawback of the Monte Carlo method is the huge number of simulations required to achieve convergence. The Monte Carlo method uses deterministic routines to solve the stochastic problem in each simulation. Analytical methods are computationally more efficient, but they need some mathematical assumptions in order to simplify the problem [14]. Therefore, Point Estimate Method (PEM) is used in this paper to estimate the information about the stochastic nature of wind power and load. The PEM is one of the most powerful approximate methods. The PEM utilizes the approximated lowerorder probability moments to assess the uncertain parameters [14]. It was initially presented by Rosenblueth [15]. In this paper 2m + 1 point estimate method is used to calculate the moments of the total operation costs which is a function of several random variables such as power of wind farms and load demand. The proposed model regards the random nature of the wind power and the load demand. Generally speaking, the MGS problem for a wind–thermal system is commonly a large-scale stochastic non-convex and

939

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

non-linear optimization problem due to the large number of thermal and wind units and decision interval. Therefore, Gravitational Search Algorithm (GSA) [16] is proposed to solve the SMGS problem. The GSA is a new population-based optimization algorithm based on the law of gravity and mass interactions. In the GSA, the searcher agents are a set of masses which interact with each other based on the Newtonian gravity and the laws of motion. The performance of the original GSA depends on its parameters such as gravitational constant. This dependency causes that the algorithm being trapped in local optima. Hence, a new Self Adaptive Mutation Strategy (SAWM) is combined with the original GSA which causes to present the powerful algorithm named Adaptive Modified Gravitational Search Algorithm (AMGSA). The SAWM effectively improves the performance of the original GSA algorithm and increases the robustness of algorithm to solve the SMGS problem. The SAWM strategy dynamically varies the mutation for algorithm along the search procedure which cause the intelligence search and impede the algorithm to fall in the local optima. On this basis, the AMGSA is able to achieve faster convergence without getting trapped in the local optima.

constant percentage of total wind power production (SRW) which regarded as an excess reserve needed to compensate for the inaccuracy made by the mismatch between the actual wind power generation and the forecasted wind generation. The reserve constraint relevant to power system is: Ng Nw X X PGRg ðtÞ  U g ðtÞ P SRW  PW w ðtÞ  V w ðtÞ þ PR ðtÞ W¼1

g¼1

ð5Þ

t ¼ 1; 2; . . . ; T The wind power is inherently non-dispatchable because its primary source of energy is uncontrollable. Due to this reason, the wind power generation is different from the power generated by thermal generation in terms of its restricted controllability. This difference is in contradiction with the maximum generation of wind farm and the available wind power of wind farm. Hence, the unequal constraint between available wind power and the actual wind power generation must be fulfilled:

PW w ðtÞ 6 W av ;w ðtÞ

ð6Þ

t ¼ 1; 2; . . . ; T Problem formulation The objective function of the MGS problem is the minimization of the total operation costs of the generating units within the scheduling horizon, while the minimization is limited to several constraints. The objective function is elucidated as the total generation costs comprising the operation and maintenance expenditure and the fuel expenses and its minimization equation is specified as follows: Ng n  T X  o   X Min J xðg; tÞ;uðg; tÞ; ad;p ¼ H PGDg ðtÞ  Nt U g ðtÞ

PLmin;l 6 PLl ðtÞ 6 PLmax;l

t¼1 g¼1

þ

PLl ðtÞ ¼

Ng n T X X PGDg ðtÞ t¼1 g¼1



Uncertainty characterization based on 2m + 1 PEM

t¼1 w¼1

A quadratic cost function will be imagined for the thermal generator that is obtained by:

ð2Þ

where Ag, Bg, Cg are the cost coefficients for the gth thermal units. The essential trouble of the MGS problem is the equilibrium between the generation and demand:

g¼1

w¼1

ð3Þ

t ¼ 1; 2; . . . ; T The minimum and maximum generation unit limits of the thermal generators also must be fulfilled:

PGmin;g 6 PGDg ðtÞ þ PGRg ðtÞ 6 PGmax;g

ð7Þ

l ¼ 1; 2; . . . ; Nl

 T X Nw  X Nt þ ð1Þ PW max;w  OMFCWðwÞ  Nh t¼1 w¼1

Ng Nw X X PGDg ðtÞ  U g ðtÞ þ PW w ðtÞ  V w ðtÞ ¼ Pload ðtÞ

w¼1

t ¼ 1; 2; . . . ; T o

T X Nw X  PW w ðtÞ  OMVCWðwÞ  Nt V w ðtÞ

HðPGDg ðtÞÞ ¼ Ag þ Bg  PGDg ðtÞ þ C g  ðPGDg ðtÞÞ2

Ng Nw X X K lg  PGDg ðtÞ  U g ðtÞ þ K lw  PW w ðtÞ  V w ðtÞ g¼1

þ PGRg ðtÞ  OMVCTðgÞ  Nt U g ðtÞ  Ng  T X X Nt PGmax;g  OMFCTðgÞ  þ Nh t¼1 g¼1 þ

In power system, transmission lines interconnect generation units located in multiple areas. Accordingly, the obtained schedule may cause some transmission lines to be overloaded if the network topology and transmission lines constraints are not taking into account in the main SMGS problem [17]. Hence, Based on DC load flow model, the transmission lines constraints will be modeled in the SMGS problem in order to decrease computation time. So, the following transmission constraint must be satisfied:

ð4Þ

In this study, there are two elements that institute the operation reserve power system. The first one is taken into account as a fixed percentage of the total load demands (PR(t)). The second one is the

The PEM is a numerical method applied to solve the probabilistic non-linear problems. The PEM calculates the statistical information of the output variables of MGS problem using the solution set of the replacement deterministic MGS problems. In the PEM method each random PDF of wind power and load is replaced by restricted locations. In this paper, PEM is used to estimate the information about the stochastic nature of wind power and load. For the first time, the PEM was developed by Rosenblueth [15]. Rosenblueth’s method required the 2m number of simulations of which the m is the number of random variables. This huge number of simulations of 2m method was computationally ineffective. Since then, several methods have been presented to improve the number of simulations to be performed [18,19]. Finally, an effective PEM were introduced by Hong which is consisted Km and Km + 1 schemes (K is a parameter depending on the type of Hong’s PEM schemes). Hong’s PEM is able to handle the symmetric and asymmetric variables [20]. The 2m + 1 scheme in comparison to 2m scheme is more accurate due to its use the Kurtosis of the input random variables. Moreover, the standard locations of PEM depend on the number of input random variables. This causes that the 2m scheme accuracy decreases when the total number of input random variables increases. Concisely, the idea behind using PEM in this paper is to efficiently and accurately estimate the uncertainty concerned with the output variables of MGS problem.

940

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

The PEM plays as an estimator of the moments of a random output variable S which is a function (F) of m uncertain input variables (zd) [14] as similar to the following equation:

S ¼ F ðz1 ; z2 ; . . . ; zd ; . . . ; zm Þ

p ¼ 1; 2; 3

ð9Þ

where nzd;p is the standard location. The values of the standard locations are computed as:

nzd;p nzd;3

kz ¼ d;3 þ ð1Þ3p 2 ¼0

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 kzd;4  k2zd;3 p ¼ 1; 2 4 ð10Þ

The variables kzd;3 and kzd;4 are the third and fourth central moments of zd which are defined by following equations:

E



kzd;3 ¼

zd  lzd  3

3

rzd

; kzd;4

h i E ðzd  lzd Þ4 ¼  4

rzd

ð11Þ

where E is the expectation operator. The weighting factor corresponding to the location zd,p can be computed as the following equation:

xzd;p ¼

ð1Þ3p   nzd;p nzd;1  nzd;2

xzd;3 ¼

1 1  m kzd;4  k2z d;3



ð12Þ



lz1 ; lz2 ; . . . ; zd;p ; . . . ; lzm p ¼ 1; 2; 3

ð13Þ

It can be concluded from Eqs. (9) and (10) that the third location of every input random variable zd is equal to its mean zd;3 ¼ lzd . Consequently, the m of the 3m locations are to be performed at the same point ðlz1 ; lz2 ; . . . ; lzd ; . . . ; lzm Þ. Therefore, it is adequate to apply only one evaluation of F function at this location. The expected value of the uth moment of the cth output random variable is computed as follows: m X 3 m X 2 X   X     E Suc ¼ xzd;p  Scðd;pÞ u ¼ xzd;p  Scðd;pÞ u d¼1 p¼1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EðS2c Þ  ðEðSc ÞÞ2 ¼ EðS2c Þ  l2Sc

ð15Þ

Adaptive modified gravitational search algorithm Standard gravitational search algorithm The Gravitational Search Algorithm (GSA) is based on the Newton’s law of gravity and law of motion [21] which developed by Rashedi in 2009 [16]. In GSA, agents are regarded as objects and their performances are evaluated by their masses. In this algorithm, each mass by a gravitational force attracts every other mass. Accordingly, position of each agent is updated after calculating acceleration and velocity of each agent by using law of motion. To explain the GSA, contemplate a space with N agents (masses). In this context, the position of the ith agents is located as:

  X i ¼ x1i ; . . . ; xbi ; . . . xni for i ¼ 1; 2; . . . ; N

ð16Þ

where xbi introduces the location of ith agent in the bth dimension into search space. The mass of each agent is computed as the following equations:

Fit i ðkÞ  worstðkÞ bestðkÞ  worstðkÞ mi ðkÞ Mi ðkÞ ¼ PN j¼1 mi ðkÞ mi ðkÞ ¼

ð17Þ ð18Þ

where Fiti(k) is the objective function of SMGS related to Eq. (1). Also, best(k) and worst(k) in minimization problems are the minimum and maximum the fitness values at generation k. The gravitational force on mass ith from mass jth is expressed as follows:

F bij ðkÞ ¼ GðkÞ

 M i ðkÞ  Mj ðkÞ  b xj ðkÞ  xbi ðkÞ Rij ðkÞ þ e

ð19Þ

where G(k) is a function of the initial value G0 and function evaluation FE. Therefore, G(k) is defined as:

p ¼ 1; 2 and

The F function is calculated for each achieved point. In each calculation of F function, one of the input random variable is fixed to one of its concentration locations and the other input random variables are assumed to be fixed to their mean value.

Scðd;pÞ ¼ F c

rSc ¼

ð8Þ

It is worthy of mention that in the SMGS problem the variable S is equal to total cost of operation of power system related to Eq. (1) and the variable m is the value of wind power of wind farms and load for every month. The 2m + 1 PEM utilizes three probability concentration locations to replace Probability Distribution Function (PDF) of input random variables by matching the first four moments of input random variables containing the mean, variance, coefficients of Skewness, and Kurtosis. Each concentration location contains two pairs (zd,p, wd,p), k = 1, 2; where zd,p and wd,p are the location and the weighting factor, respectively. The weighting factor ascertains the influence of the corresponding location in estimating the statistical moments of the output random variables. The three concentration locations for each random variable are calculated as follows:

zd;p ¼ lzd þ nzd;p  rzd

lSc ¼ EðSc Þ

 FE GðkÞ ¼ G0  exp D  FEmax

where D is a constant term. Rij(k) is the Euclidian distance between two masses i and j which is as follows:

  Rij ðkÞ ¼ X i ðkÞ; X j ðkÞ2

F bi ðkÞ ¼

X

d¼1

The mean and the Standard Deviation (SD) values for cth random output variable are as:

randj F bij ðkÞ

ð22Þ

j2qbest;j–i

qbest is the set of first q agents in the population that have the biggest mass and the best fitness value. qbest is a function of time. qbest initialized to total number of agents at the beginning and is decreased linearly to 1 with time. The acceleration  of the agent i, b in the bth dimension and at generation k acceli ðkÞ is: b

ð14Þ

ð21Þ

by a randomly weighted sum of the bth components of the forces from other qbest agents, the total force applying on the agent i in   the bth dimension F bi ðkÞ is computed below:

d¼1 p¼1

m h  u iX þ F lz1 ; lz2 ; . . . ; lzd ; . . . ; lzm xzd;3

ð20Þ

acceli ðkÞ ¼

F bi ðkÞ M i ðkÞ

ð23Þ

Consequently, the updated velocity of an agent is regarded as its acceleration added to a fraction of its current velocity. Hence, its updated velocity and position are calculated as:

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

v elbi ðk þ 1Þ ¼ randi  v elbi ðkÞ þ accelbi ðkÞ b xbnew;i ðk þ 1Þ ¼ xbi ðkÞ þ v eli ðk þ 1Þ

ð24Þ ð25Þ

After the new mutant agent was generated, the fitness value of the  mutant agent xmut;i ðk þ 1Þ should be compared with the fitness   value of the existing vector xnew;i ðk þ 1Þ for the replacement operation:

(

Self Adaptive Wavelet Mutation Strategy (SAWM)

X new;i ðk þ 1Þ ¼ In this paper a novel powerful SAWM strategy is proposed for improving the performance of the original GSA and to achieve faster convergence without getting trapped in the local optima. The SAWM operation dynamically varies the mutating space along the search procedure based on the Wavelet theory [22]. In the SAWM strategy, all elements of each agent will have a chance to be mutated, managed by a probability of mutation (Pmut). The value of Pmut is between 0 and 1. For each element of particle, a random number between 0 and 1 will be produced. If the random number is less than or equal to Pmut , a mutation will occur on corresponding element. This method is based on the fact that the total number of agents for their performance promotion tend to move in a similar way with the best(k) and avoid to move in the path of   worst(k). A mutant agent xbmut ðk þ 1Þ is generated as: ( b xnew;i ðk þ 1Þ þ -ðbest ðkÞ  xbnew;i ðk þ 1ÞÞ if - > 0 xbmut;i ðk þ 1Þ ¼ ð26Þ xbnew;i ðk þ 1Þ þ -ðworst ðkÞ  xbnew;i ðk þ 1ÞÞ if - 6 0 1 u - ¼ wh;0 ðuÞ ¼ pffiffiffi w ð27Þ h h where w(u) is Morlet wavelet function which is defined as follows: ð/Þ2 2

wðuÞ ¼ exp

cosð5uÞ

ð28Þ

Thus - is:

1

h

2

ðuh Þ



- ¼ pffiffiffi exp

2

 u cos 5 h

ð29Þ

A larger value of |-| makes large changes in the mutation and a large searching space. Conversely, a small value of |-| makes a large searching space. If r is positive, the mutant agent moves toward the location of global best solution and if r is negative, the mutant agent recedes from the location of worst solution. u can be randomly generated from the interval [2.5 h, 2.5 h] due to over 99% of the total energy of the mother wavelet function is located in the interval [2.5, 2.5]. The value of the dilation parameter (h) is set to vary with time in order to meet the fine-tuning purpose. In this context, at the early function evaluations of the evolution, the mutation operator executes a broader search in the search space and at the later steps, the search is restricted around the global location. Specifically, the value of h at the early function evaluations is small which causes the maximum value of - will become large and with the passage of time the value of h is larger that the maximum value of - will become very small. As a result, h is given as:

h ¼ exp lnðgÞð1FE=FEmax Þ

b

þlnðgÞ

ð30Þ

where g and b are the upper limit and shape parameter of increasing function of h, respectively. The value of b critically affects the performance of the AMGSA. Hence, In order to reach a desirable balance between the exploration and exploitation capability of the AMGSA, the value of b should be chosen correctly. Therefore, in the initial stages of the search procedure, a small value of b should be used to increase the step size of the mutation in order to reach the optimal point and a large value of b is employed to fine-tune faster at the end of search procedure. Thus, b is proposed as follows:

b ¼ bmin þ

bmax  bmin  FE FEmax

ð31Þ

941

    X mut;i ðk þ 1Þ if Fit X mut;i ðk þ 1Þ 6 Fit X new;i ðk þ 1Þ X new;i ðk þ 1Þ otherwise ð32Þ

Implementation of AMGSA to SMGS problem In this section, step-by-step methodology for solving the SMGS problem is provided. The flowchart of the whole process of proposed method is given in Fig. 1. Step 1: Input the preliminary data of the SMGS problem. Step 2: Execute 2m + 1 PEM scheme to determine the three concentration locations and weighting factors of each random variable. Step 3: Characterize the input variables and zd,p of F function for dth random variable. Step 4: Randomly initialize the locations of the masses for the population. The whole agent must satisfy the quality and inequality constraints of MGS problem. Step 5: Compute objective functions for each agent. Step 6: Update G(k), best(k), worst(k) and Mi(k) for i = 1, 2, . . . , N. Step 7: Calculate the total force in different directions. Step 8: Calculate the acceleration and the velocity of each agent. Step 9: Update the location of masses. Step 10: Perform SAWM strategy on all the existing agents. The new agent set must satisfy the quality and inequality constraints of MGS problem. Step 11: Repeat step 5 to step 10 until the current function evaluation number reaches the pre-specified maximum number of function evaluations. Step 12: Go to Step 3 until the whole random variable is considered. Step 13: Appraise the first and the second moments of the objective functions based on 2m + 1 PEM scheme. Case studies In this section, two case studies are investigated for the SMGS problem to display the usefulness of the proposed method in the practical applications. The scheduling horizon of SMGS for all cases is considered for one year and the decision stage is on the monthly basis. The studies have been carried out on MATLAB 7.10 using a Pentium P4, Dual core 2.21-GHz personal computer with 1 GB of RAM. Description of the case studies Case 1: the first case consists of 10 thermal units and two wind farms (units 11 and 12). The monthly total load of this test system is 1500. The input data for thermal units are presented in [23,24]. Case 2: in this case, in order to illustrate the scalability of the proposed approach, the large-scale 42-unit is selected. This case comprises 40 thermal units and two wind farms (units 41 and 42) and the monthly peak load is assumed to be 9500 MW. The system data is taken from [24,25]. The PDF load demand is considered as normal distribution function. The mean and standard deviation values for monthly peak load in each interval for these cases are featured in [24]. It is wor-

942

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

Fig. 1. Flowchart of the proposed stochastic method using AMGSA algorithm.

thy of mention that due to none existence of any test case to be available the pertinent data for its generators characteristics and network topology, the transmission lines constraints for existing test cases have not been considered in simulation. However, the proposed model has the capability for the application of these constraints. General experimental setting In order to ensure a good convergence of the AMGSA, the parameters G0 ; D and e are 100, 20 and 0.00001, respectively. The mutation probability and parameters g; bmax and bmin of SAWM strategy are set to 0.01, 1000, 5 and 0.5, respectively which have chosen by trial and error through experiments for good performance. It should be noted that the first component of reserve requirement PR(t) is presumed to be 5% of total load demand in

each time interval and the second one SRW is fixed to 10% of total wind power generation for two case studies [24]. Wind turbine generation pattern modeling In this study, for a realistic assessment of wind power, the wind speeds are parameterized using Weibull PDF [26]. Once the uncertain nature of the wind is described as a random variable, the output power of the wind farm through a conversion from wind speed to wind power may also be described as a random variable. A simplified linear piecewise function is employed to expand the relation between wind speed and output of wind generator [27]. For the sake of simplicity, the employed data for the wind farms of all test cases are identical. The Weibull parameters of the wind speed in each wind farm for all test cases are presumed to be c = 15, k = 2.2. The wind power generation unit’s model will be

943

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

employed to model a wind farm rated at 80 MW, with cut-in, rated and cut- out wind speeds of 2.5, 14, and 25 m/s, respectively. Comparison AMGSA with other algorithms Statistical comparative The performance of the original GSA, the Differential Evolution (DE), Differential Evolution with SAWM strategy (DE-SAWM), the Particle Swarm Optimization (PSO) and the proposed AMGSA on solving the SMGS problem for two cases is evaluated. The values of parameters CR and F are taken into account for DE algorithm are 0.88 and 0.47, respectively. The value of parameters C1 and C2 of PSO algorithm are 1.2 and 2.3, respectively. The employed setting for GSA and DE-SAWM are the setting used for AMGSA and DE, respectively. For all of the aforementioned test cases, the algorithms perform 25 independent runs. The numbers of the population of all algorithms are equal to 30 and 50 for the test Cases 1 and 2, respectively. The maximum number of Function Evaluations (FEmax) is set to 30,000 and 50,000 for the test Cases 1 and 2, respectively. Table 1 gives the statistical comparative results of the proposed algorithm versus the other methods such as mean, standard deviation, best, worst solutions and the Deviation of the Best (DB) and Deviation of the Worst solutions (DW) [28]. The best results are typed in bold. Comparing the mean solutions of the proposed method with other methods corroborates the superiority of AMGSA over the other methods. This Table shows that the proposed method reaches a lower total cost than GSA (M$255.112 opposed to M$255.431) in lower SD, DB and DW. In Table 1, the best solution of the AMGSA without any exception is better than the best solution of all other algorithms in two cases. These comparisons unveil the capability of the AMGSA to solve the applicable MGS problem in different test cases. The corresponding SD of AMGSA, GSA, DE-SAWM, DE and PSO are T$0.095, T$0.132, T$0.155, T$0.173 and T$0.235 for test Case 2. Referring to property of the wavelet theory, the summation of the positive value of - is equal to the summation of the negative value of - when u is randomly generated and the number of samples is large. Therefore, the total number of positive mutation of SAWM and the total number of negative mutation of SAWM throughout the evolution are approximately the same. By considering the smaller standard deviations of AMGSA than other algorithm, perception that this property of the wavelet provides a quality and stable solution. Moreover, Table 1 displays the CPU execution time of different optimization algorithms for test cases. As shown in this table, in the comparison with other optimization algorithm, the proposed AMGSA after the GSA has the lower CPU time in the same condition. But, by spending a few times than the GSA, the AMGSA can reach the significant and salient lower

operation cost than GSA and other optimization algorithms. This significant lower operation cost is evident in both test cases. On the other hand, mid-term generation scheduling is carried on the one year with monthly intervals and also CPU time of AMGSA and entire algorithms are very lower than one month. Hence, the operation cost is in high level of consideration, in trade off CPU time and operation cost witch yielded by algorithms. Accordingly, it is rational and is suggested that AMGSA is befitting algorithm in solving SMGS problem. Figs. 2 and 3 represent the mean convergence characteristics of the proposed algorithm and other methods for the two test cases. As can be inferred from these figures, it is evident that AMGSA algorithm has the capability to relax the stagnation and skip the local optima point with smaller standard deviations, demonstrating better efficiency and stability than the other three algorithms. t-Test To appraise the significant difference between AMGSA and the other optimization algorithms the t-test method is implemented [29]. If the first algorithm is better than the second one, the t-value will be positive and it is negative if the first algorithm is trivial. In Table 2, the t-values with 24° of freedom at a 0.05 level of significance are supplied to statistically display the comparison between AMGSA and the other algorithms. From this Table, we can perceive that the AMGSA is the best in terms of t-values. It is obvious that the t-values between the AMGSA and other optimization algorithms are enormously higher than 2.06, implying that the AMGSA is significantly better with a 98% confidence level, than other algorithms. Empirical distribution of normalized success performance In order to conclusively compare the overall performances of all algorithms on the whole test case; we plot the empirical distribution of normalized success performance [30]. We first computed the successful runs of all algorithms on each test case. The successful run of an algorithm means that this algorithm can discover at least one solution whose fitness value is not worse than pre-determined value. The success rate of algorithm is computed as the number of successful runs divided by the total number of independent runs. Moreover, the Average Number of Function Evaluations (ANFE) that needed finding the pre-designate value, is docketed when an algorithm solves the SMGS problem with 100% success rate [29]. The ranking of the optimization algorithms is based on the Successful Performance (SP) of an algorithm, which is specified as:

SP ¼ meanðNFE of successful runsÞ 

Total number of all runs Number of successful runs

ð33Þ

Table 1 Results obtained by different algorithms for Cases 1 and 2. Text Case

Solution Technique

Best

Mean

Worst

SD

DB (%)

DW (%)

CPU time (h)

12-unit

Total operation cost (M$) PSO DE DE-SAWM GSA AMGSA

269.458 263.815 255.993 255.431 255.112

271.454 266.312 260.549 256.175 255.769

273.475 268.568 262.415 257.148 257.148

0.531 0.364 0.319 0.167 0.038

0.735 0.938 1.749 0.290 0.257

0.745 0.847 0.716 0.623 0.539

4.582 4.568 4.597 4.515 4.576

Total operation cost (T$) PSO DE DE-SAWM GSA AMGSA

17.964 16.958 16.745 16.323 15.714

19.348 17.258 16.869 16.451 15.825

20.962 17.839 17.498 17.004 16.347

0.235 0.173 0.155 0.132 0.095

7.153 1.738 0.735 0.778 0.701

8.342 3.367 3.729 3.361 3.299

28.642 28.374 29.151 27.361 27.817

42-unit

The best results of each indices are typed in bold.

944

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

Fig. 2. Mean convergence characteristics of AMGSA, GSA, DE-SAWM, DE and PSO for Case 1.

Fig. 3. Mean convergence characteristics of AMGSA, GSA, DE-SAWM, DE and PSO for Case 2.

Table 2 t-Value between AMGSA and other optimization methods. Test case

t-Value between AMGSA and PSO

t-Value between AMGSA and DE

t-Value between AMGSA and DE-SAWM

t-Value between AMGSA and GSA

Case I Case II

147.32 74.09

144.04 40.53

74.39 32.67

11.85 22.76

Next, all SPs are divided by the best SP on respective test case to achieve the normalized SP. The outcomes of all test cases are used where at least one optimization algorithm was successful in at least one run. Tables 3 and 4 present detailed information on the number of function evaluations needed to reach the tolerance limit and the success rate. As can be observed from these tables, the AMGSA performs much better with success rates of 100% on all test cases. It is worthy of mention that for the algorithms that their success rate equal to zero, the mean and SD of FEs will not be defined for them. The graph of empirical cumulative distribution for all test cases is plotted in Fig. 4. In this figure, small values of normalized successful performance and large values of the empirical distribution are

Table 3 Performance of various optimization algorithms for Case 1. Solution technique

Mean of FEs

SD of FEs

ANFE

Success rate (%)

PSO DE DE-SAWM GSA AMGSA

– – 27,810 23,301 19,721

– – 3543 2523 2307

– – – 23,301 19,721

0 0 12 100 100

The best results of each indices are typed in bold.

Table 4 Performance of various optimization algorithms for Case 2. Solution technique

Mean of FEs

SD of FEs

ANFE

Success rate (%)

PSO DE DE-SAWM GSA AMGSA

– – – 45,023 37,382

– – – 2691 2480

– – – – 37,382

0 0 0 32 100

The best results of each indices are typed in bold.

superior. The first algorithm that attains the top of the graph will be regarded as the best optimization algorithm. It is evident that the empirical distribution of AMGSA attains its maximum of one at the smaller value of normalized successful performance than other algorithms, demonstrating a higher efficiency in solving SMGS problems.

Analysis of stochastic results This part investigates the stochastic results which founded by AMGSA for Cases 1 and 2. The elaborate illustration of the best dispatch gained by AMGSA for SMGS problem for test Case 1 is also given in Table 5. This Table indicated that the units with smaller

945

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

Fig. 4. Empirical distribution of normalized success performance on Cases 1 and 2.

Table 5 Best dispatch of SMGS problem found by AMGSA for Case 1. Time (Month)

1 2 3 4 5 6 7 8 9 10 11 12

Power of generation units (MW) Unit 1

Unit 2

Unit 3

Unit 4

Unit 5

Unit 6

Unit 7

Unit 8

Unit 9

Unit 10

455 455 455 455 455 455 455 455 455 455 455 455

455 455 383.180 455 455 455 455 431.618 410.215 455 455 455

130 130 130 130 130 130 130 124.190 123.016 130 130 130

130 130 127.384 130 130 130 130 125.867 130 130 130 130

113.557 111.562 0 52.607 128.789 130 106.481 25 25 79.877 162 162

0 20.726 0 0 20.125 22.258 25.476 0 0 50.566 38.290 80

0 0 0 0 0 0 0 0 0 0 0 25.981

13.571 0 0 0 0 0 0 0 0 0 11.988 20.370

0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0

Table 6 Stochastic results founded by AMGSA for Cases 1 and 2. Test case Case I Case II

Expected

SD

0.95% confidence intervals

255.769

1.030

0.403

15.825

0.145

0.056

VSS (%)

Total operation cost (M$) 1.051 Total operation cost (T$) 1.780

Relative error (%)

Coefficient of variation (%)

0.157

0.402

0.353

0.916

Fig. 5. Cumulative probability distribution of the total yearly operation for Cases 1 and 2.

average power generation costs are loaded first and leaving the high-cost units are loaded for the power generation of load peak. Table 6 is exhibited the indexes such as, first moment, standard deviation, the value of expected, relative error [31], 95% confidence intervals, coefficient of variation [32] and the Value of the Stochastic Solution (VSS) [33] for the outcome of each test case. The first order and the standard deviation of the operation costs of the best solutions for Case 1 are M$255.769 and M$1.030, respectively.

The VSS index specifies how much the SO, in order to acquire a realization of the future of the uncertainty variables by the stochastic PEM, should spend. For example, the SO would be inclined to pay 1.051% of deterministic cost more than deterministic cost to reveal uncertainties of cost load demand and wind farm’s power generation in the SMGS problem of test Case 1. Although, considering the uncertainties imposes extra cost to the operation and management system, it provides a more realistic condition for the parameters of power system. It should be maintained that lowering the value of

946

T. Niknam, H.R. Massrur / Electrical Power and Energy Systems 64 (2015) 937–946

confidence interval of operation costs expresses the accurate level of the expected value of operation cost. The values of confidence interval for test Cases 1 and 2 are 0.403 and 0.056, respectively. Also, the more effectiveness of the proposed method, is causes the lower rate of the relative error. Finally, based on the Grame-Charlier expansion [34] and by utilizing the first and the second moments of the yearly operation costs which were founded by the 2m + 1 PEM, the Cumulative Distribution Function (CDF) of yearly operation cost is presented in Fig. 5. The CDF of operation costs helps the SO to understand how the uncertain variables affect the system. Conclusions This paper presented an effective stochastic 2m + 1 PEM model to address the influence of load demand and wind power generation uncertainties on the yearly generation scheduling of the power system. The power system operators, by using this robust stochastic generation scheduling can accommodate the intermittent nature of the wind generation. In order to solve the proposed stochastic model of MGS problem, the proposed evolutionary algorithm named AMGSA was employed. Moreover, In order to improve the robustness and the efficiency of the proposed approach, a novel Self-Adaptive Wavelet Mutation based on Wavelet theory have been implemented with the AMGSA. We have compared the performance of AMGSA with the other optimization algorithms and concluded that AMGSA was more efficient in acquiring better quality solutions which are more stable with the relatively smaller standard deviation. It had higher success rates and superior empirical distribution of success performance. The Application of the stochastic method demonstrated that although the consideration of uncertainties imposes extra cost on the scheduling system than the earlier deterministic approaches, it provides a more realistic understanding of the power system’s parameters. Future works will concentrate on modeling issues in connection to uncertainty derived from electricity and fuel price forecasting of power system. Appendix A. Calculation of kln The DC load flow elucidates the net real power injection P to bus voltage angles d by disregarding transmission line resistances and supposing voltage magnitudes to be 1 per unit. The matrix equation is given as follows:

P ¼ B0  d ! d ¼ A0  P

ð34Þ 0 matrix, bij ¼  x1ij , for i – ref Pn 0 0 bii ¼  j¼1;j–1 for i – ref; bij

0

where B is the bus susceptance and 0 j – ref; bij ¼ 0 for i = ref or j = ref; ¼0 for i = ref. B0 is a singular matrix. Hence, if there are NB buses, we only have NB  1 linearly independent equations. Accordingly, we have a (NB  1)  (NB  1) sub-matrix of B0 , which can be inverted. A0 represents the inverse of the sub-matrix of B0 plus a zero row and column corresponding to the reference bus. The power flow Pi,j is given as:

Pi;j ¼

di  dj xij

ð35Þ

If we use ai and aj symbolizing the ith and jth rows of A0 , then;

 Pi;j ðtÞ ¼

NB X ai  aj  PðtÞ ¼ kln  Pn ðtÞ xij n¼1

ð36Þ

therefore

K ln ¼ kij ¼

ain  ajn xij

ð37Þ

References [1] Fosso B, Gjelsvik A, Haugstad A, Mo B, Wangensteen I. Generation scheduling in a deregulated system. The Norwegian case. IEEE Trans Power Syst 1999;14:75–81. [2] Quintana VH, Chikhani AY. A stochastic model for mid-term operation planning of hydro-thermal systems with random reservoir inflows. IEEE Trans Power App Syst 1981;100:1119–27. [3] Basu M. Improved differential evolution for short-term hydrothermal scheduling. Int J Elect Power Energy Syst 2014;58:91–100. [4] Lakshmi K, Vasantharathna S. Gencos wind–thermal scheduling problem using artificial immune system algorithm. Int J Elect Power Energy Syst 2014;54:112–22. [5] Gjorgiev B, Kanceva D, Cepin M. A new model for optimal generation scheduling of power system considering generation units availability. Int J Elect Power Energy Syst 2013;47:129–39. [6] Mahari A, Zare K. A solution to the generation scheduling problem in power systems with large-scale wind farms using MICA. Int J Elect Power Energy Syst 2014;54:1–9. [7] Zhou J, Liao X, Ouyang S, Zhang R, Zhang Y. Multi-objective artificial bee colony algorithm for short-term scheduling of hydrothermal system. Int J Elect Power Energy Syst 2014;55:542–53. [8] Jiekang W, Zhuangzhi G, Fan W. Short-term multi-objective optimization scheduling for cascaded hydroelectric plants with dynamic generation flow limit based on EMA and DEA. Int J Elect Power Energy Syst 2014;57:189–97. [9] Kumar ABR, Vemuri S, Ebrahimzadeh P, Farahbakhshian N. Fuel resource scheduling – the long-term problem. IEEE Trans Power Syst 1986;1:145–51. [10] Baslis CG, Bakirtzis AG. Mid-term stochastic scheduling of a price-maker hydro producer with pumped storage. IEEE Trans Power Syst 2011;26:1856–65. [11] Wei Q, Jinfeng L, Christofides PD. Supervisory predictive control for long-term scheduling of an integrated wind/solar energy generation and water desalination system. IEEE Trans Control Syst Technology 2012;20:504–12. [12] Franceschini S, Marani M, Tsai C, Zambon F. A perturbance moment point estimate method for uncertainty analysis of the hydrologic response. Adv Water Resour 2012;40:46–53. [13] Tung YK, Yen BC. Hydro-systems engineering uncertainty analysis. New York: McGraw-Hill; 2005. [14] Caro E, Morales JM, Conejo AJ, Minguez R. Calculation of measurement correlations using point estimate. IEEE Trans Power Del 2010;25:2095–103. [15] Rosenblueth E. Point estimates for probability moments. Proc Nat Acad Sci 1975;72:3812–4. [16] Rashedi E, Nezamabadi-pour H, Saryazdi S. GSA: a gravitational search algorithm. Inform Sci 2009;179:2232–48. [17] Marwali MKC, Shahidehpour SM. A deterministic approach to generation and transmission maintenance scheduling with network constraints. Electr Power Syst Res 1998;47:101–13. [18] Li KS. Point estimate method for calculating statistical moments. J Eng Mech 1992;118:1506–11. [19] Harr ME. Probabilistic estimates for multivariate analysis. Appl Math Model 1989;13:313–8. [20] Hong HP. An efficient point estimate method for probabilistic analysis. Reliab Eng Syst Saf 1998;59:261–7. [21] Holliday D, Resnick R, Walker J. Fundamentals of physics. New York: John Wiley & Sons; 1993. [22] Rioul O, Vetterli M. Wavelets and signal processing. IEEE Signal Process Mag 1991;8:14–38. [23] Cheng CP, Liu CW, Liu CC. Unit commitment by Lagrangian relaxation and genetic algorithms. IEEE Trans Power Syst 2000;15:707–14. [24] Siahkali H, Vakilian M. Electricity generation scheduling with large-scale wind farms using particle swarm optimization. Elect Power Syst Res 2009;79:826–36. [25] Chen PH, Chang HC. Large-scale economic dispatch by genetic algorithm. IEEE Trans Power Syst 1995;10:1919–26. [26] Xu M, Zhuan X. Identifying the optimum wind capacity for a power system with interconnection lines. Int J Elect Power Energy Syst 2013;51:82–8. [27] Sobolewski RA, Feijóo AE. Estimation of wind farms aggregated power output distributions. Int J Elect Power Energy Syst 2013;46:241–9. [28] Amjady N, Nasiri-Rad H. Nonconvex economic dispatch with AC constraints by a new real coded genetic algorithm. IEEE Trans Power Syst 2009;24:1489–502. [29] Wang Y, Li B, Weise T, Wang J, Yuan B, Tian Q. Self-adaptive learning based particle swarm optimization. Inf Sci 2011;181:4515–38. [30] Qin AK, Huang VL, Suganthan PN. Differential evolution algorithm with strategy adaptation for global numerical optimization. IEEE Trans Evolut Comput 2009;13:398–417. [31] Wu L, Shahidehpour M, Li T. Stochastic security constrained unit commitment. IEEE Trans Power Syst 2007;22:800–11. [32] Niknam T, Azizipanah-Abarghooee R, Narimani MR. An efficient scenariobased stochastic programming framework for multiobjective optimal microgrid operation. Appl Energy 2012;99:455–70. [33] Niknam T, Massrur HR, Bahmani-Firouzi B. Stochastic generation scheduling considering wind power generators. J Renew Sustain Energy 2012;4:063119. http://dx.doi.org/10.1063/1.4767930. [34] Zhang P, Lee S. Probabilistic load flow computation using the method of combined cumulants and Grame-Charlier expansion. IEEE Trans Power Syst 2004;19:676–82.