Multi-objective optimal strategy for generating and bidding in the power market

Energy Conversion and Management 57 (2012) 13–22

Contents lists available at SciVerse ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Multi-objective optimal strategy for generating and bidding in the power market Chunhua Peng ⇑, Huijuan Sun, Jianfeng Guo, Gang Liu Department of Electrical & Electronics Engineering, East China Jiaotong University, Nanchang, Jiangxi Province 330013, China

a r t i c l e

i n f o

Article history: Received 11 August 2011 Received in revised form 6 December 2011 Accepted 6 December 2011 Available online 2 January 2012 Keywords: Multi-objective optimization Differential evolution Non-dominated sorting Electricity market Generation bidding

a b s t r a c t Based on the coordinated interaction between units output and electricity market prices, the benefit/risk/ emission comprehensive generation optimization model with objectives of maximal profit and minimal bidding risk and emissions is established. A hybrid multi-objective differential evolution optimization algorithm, which successfully integrates Pareto non-dominated sorting with differential evolution algorithm and improves individual crowding distance mechanism and mutation strategy to avoid premature and unevenly search, is designed to achieve Pareto optimal set of this model. Moreover, fuzzy set theory and entropy weighting method are employed to extract one of the Pareto optimal solutions as the general best solution. Several optimization runs have been carried out on different cases of generation bidding and scheduling. The results confirm the potential and effectiveness of the proposed approach in solving the multi-objective optimization problem of generation bidding and scheduling. In addition, the comparison with the classical optimization algorithms demonstrates the superiorities of the proposed algorithm such as integrality of Pareto front, well-distributed Pareto-optimal solutions, high search speed. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction In order to improve efficiency of electric power production in a coordinated and sustainable way for the energy production and environment protection, it is very important for power industry, especially for generation side, to focus on marketization, energy saving and emissions reduction. By the working of the competitive mechanism with the purposes of energy saving, environmental protection and economy, the generation order and time of each kind of generation unit can be determined and then the units with low energy consumptions, emissions and costs would have priority to be chosen. It can promote healthy development of electric power production. In electric power market environment, in order to maximize selling profits, the general method taken by generation companies can be described as follows. First, determine the electricity price trend by electricity price short-term forecasting model [1]. Second, determine reference of bidding price of each period. And then globally optimize unit output based on the reference to achieve optimal generation bidding schedule by period [2]. However, the uncertain behavior of supply and demand market causes uncertain fluctuation of electricity price. That results in the contradiction of high bidding price with high risk and low bidding price with low benefits. To guarantee the success of bidding for generation companies, the bidding risk should be reasonably estimated to ensure the balance between risk and profit [3,4]. Moreover, as a main air ⇑ Corresponding author. Tel.: +86 07917046236; fax: +86 07917046184. E-mail address: [email protected] (C. Peng). 0196-8904/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.enconman.2011.12.006

pollution source, thermal power plants should minimize emissions to satisfy the requirement of environmental protection [5,6]. Therefore, for coordinated and sustainable development of energy production and environment, and also for the balance between profit and social responsibility, it is important to design optimal bidding strategy with minimal bidding risk, minimal emissions and maximal profit in electric power market environment. Many researches on bidding strategy for generation side are based on game theory [7,8]. Oligopoly game is mostly used to find market balanced point to achieve optimal bidding strategy. For using game theory, models are usually simplified and problems are transferred to game with complete information. That is, in the process, cost functions and bidding strategy of each unit were regarded as common knowledge. However, such information is confidential to most of generation companies and can hardly be achieved. That is to say, those models are more suitable for analyzing market potential than for determining actual generation bidding schedule in electric power market. For generation dispatch, current researches mostly focused on the optimization of load dispatch considering power flow, purchase cost and environment protection from the view of whole power system [9]. In some references, which analyze the relation model between emissions and output, the emission models are established either separately according to the main harmful gases such as SO2, CO2, NOx, or comprehensively [10]. However, they all ignored both bidding risk and global optimization of unit output in the dealing day. Studies on multi-objective optimization of unit operation and generation bidding that aim at realizing energy saving, environmental protection, low bidding risk and high profit while comprehensively considering

14

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

the dynamic relationship among forecasting price in the dealing day, bidding risk of generation companies, emission control and unit output, are still not enough. So, in order to balance the relationship among generation, electricity market price and environmental protection and to guarantee flexibility and effectivity of bidding strategy of generation side, a benefit/risk/emission generation optimization model based on market price forecasting curve in the dealing day for maximizing profit and minimizing both risk and emissions, with consideration of atmospheric emission and bidding risk coefficient, is established. The multi-objective optimal unit operation plan can be achieved from the model and thus to form the optimal generation bidding schedule of generation company. Since the benefit/risk/emission generation optimization model has characteristics of unsmooth, non-convex and objectives coupling [2,3,9], and the scales of each objective are different, there is usually no such absolute optimal solution which can optimize all objectives at the same time. Only compromised solutions could be made to optimize each objective as much as possible through coordination and compromise. It is difficult to solve such multiobjective problems by traditional weighted methods which convert a multi-objective optimization problem (MOP) into a single objective optimization problem. Following the development of multi-objective evolution algorithm, a series of multi-objective optimization algorithms, which are based on genetic algorithm, particle swarm optimization, ant colony optimization, etc., have been proposed and have made considerable progress [11,12]. However, while applying them in the problem of benefit/risk/emission comprehensive generation optimization, we found there are some drawbacks such as unstable convergence, unevenly distribution of non-dominated solutions and low speed. So, a more effective and simpler algorithm, the improved multi-objective differential evolution algorithm, is proposed in this paper for this problem. The validity of the proposed algorithm is verified by simulations and applications. 2. Multi-objective problem of generation bidding 2.1. Objective functions 2.1.1. Objective function of bidding risk Electricity market price is the organizational core of electricity exchange in market environment. Market participants should choose their bids in order to maximize their profits, and the bidding reference can be determined by short-term forecasting of electricity market price. The characteristic of competitive bidding in power market will lead to the contradiction of high bidding price with high risk and low bidding price with low profit. Here, as a schedule cycle, one dealing day is divided into 24 periods. After statistically analyzing the historical data of electricity price by period, it can be considered that the actual market marginal price in   period n obeys a normal distribution of N ln ; r2n [13,14], where ln is the forecasted market marginal electricity price in period n, r2n is the variance of price in period n and can be achieved by analyzing historical data. Assuming the actual bidding price is pn, the objective function of minimal bidding risk in period n can be calculated by the probability of bidding failure (when pn is higher than the actual market marginal price):

1 min:v n ¼ Fðpn jln ; rÞ ¼ pffiffiffiffiffiffiffi r 2p

Z

pn



e

ðxln Þ2 2r 2

dx

price should be normally lower than forecasted price. That means

vn should be less than 0.5, otherwise the probability of bidding failure would increase rapidly. 2.1.2. Objective function of generation profit In electricity market, generation companies can determine optimal output of each unit according to unit operation status and electricity market price. Therefore, the relationship between each unit is relatively independent. Each unit can be modeled separately to realize generation optimization. In this paper, one dealing day is considered as a dispatch cycle and divided into 24 periods. After comprehensively considering generation fuel cost, unit startup/shutdown cost, fixed cost and bidding risk, the objective function of maximizing generation profit for one unit in one dealing day is [2,15]:

min:  Gðu; qÞ ¼ 

n¼1

 Cfa  un1 ð1  un ÞD  un ð1  un1 ÞSg

ð2Þ

where G is generation profit of one unit in one dealing day; xor() is exclusive-or operation; qn is actual unit average output in period n (equals to generated energy in period n in numerical); qmax is maximum power output of the unit; rn is spinning reserve capacity price; un is 0/1 variable which represents the operation status of the unit where 0 means shutdown and 1 startup. If un1 is 1 while un is 0, then there will be shutdown cost. And if un1 is 0 while un is 1 then there will be startup cost. Cfa is average fixed cost and can be achieved by dividing the total fixed cost into parts corresponding to periods; D is shutdown cost and can be regarded as a constant; fuel cost Cn and unit startup cost S can be calculated by (3) and (4), respectively [16].

C n ¼ aq2n þ bqn þ c þ jg sinðhðqn  qmin ÞÞj

ð3Þ

S ¼ K 0 þ K 1 ð1  eT=s Þ

ð4Þ

where a, b and c represent characteristics parameters of fuel cost curve; qmin is minimum power output of the unit; g and h are valve point effect parameters; K0 is the turbine startup cost and K1 is the startup cost while the boiler has been fully cooled; T is shutdown duration and s is boiler cooling time constant. The pn in (2) can be represented as bidding risk coefficient by the inverse function of (1):

pn ¼ F 1 ðv n jln ; rÞ ¼ fpn : Fðpn jln ; rÞ ¼ v n g

ð5Þ

For the market with marginal electricity price as its unified purchase price, the unit outputs can be optimized only by the forecasted market marginal prices of 24 periods. 2.1.3. Objective function of atmospheric pollutants The atmospheric pollutants caused by power plant mainly include CO2, SO2, NOx, and so on. The relationship between each atmospheric pollutant and generation can be modeled separately. Here, we adopt the atmospheric pollutants comprehensive emission model [9]. The objective function that minimizes atmospheric pollutants E in one dealing day is:

ð1Þ

1

where vn is defined as coefficient of bidding risk in period n, which represents the probability of bidding failure and is valued between 0 and 1. From (1), it is clear that there are two ways to decrease bidding risk, one is to decrease the bidding price and the other is to increase the forecasting accuracy of electricity price. Actual bidding

24 X fxorðqn ; 0Þ½qn pn þ ðqmax  qn Þr n  C n 

min: E ¼

24 n o X   102 aq2n þ bqn þ c þ n expðkqn Þ

ð6Þ

n¼1

where ai, b, c, n and k are coefficients of emission characteristics of the ith generator and can be achieved by applying least squares method on monitoring data of unit emission.

15

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

2.2. Constraints

3. Design of multi-objective differential evolution algorithm

(1) Generation limits constraint:

qmin 6 q 6 qmax

3.1. Description of multi-objective optimization problem

ð7Þ

where qmin and qmax are the minimum power output and capacity of the unit, respectively.

8 (1) h Minimum i online time and offline time constraints: > < T Un1  T Umin ½un1  un  P 0 h i ; ðn ¼ 1; 2; . . . ; 24Þ > : T Dn1  T Dmin ½un  un1  P 0

minfi ðXÞ; ð8Þ

where T Un1 and T Dn1 are continuous online time and offline time of the unit until period n1 respectively, T Umin and T Dmin are minimum online time and minimum offline time of the generating unit respectively, and all in units of 1 h normally. To avoid that the unit output changes too rapidly or the unit is always in the process of start-up or shut-down, the output variation rate R has to be constrained:



Rd 6 R 6 Ru ; if q P qmin R0 6 jRj 6 R1 ; if 0 6 q < qmin

2.3. Analysis of average output in one period Suppose the power output at the beginning of period n is qn,b and current status (online or offline) has been maintained for Tc, after comprehensively considering all the constraints, on/off states of generating units at the beginning of period n and the range of output at the end of period n (that is qn,e) can be determined as follows. Also, qn,e can be used as the initial power output at the beginning of the next period (qn+1,b).  When qn,b P qmin, the unit is online, and is allowed to shutdown if TcPT Umin and then the range of qn,e is [max(qn,b  Rd, 0), min(qmax, qn,b + Ru)]; Otherwise, the unit is not allowed to shutdown and the range of qn,e is [max (qn,b  Rd, qmin), min (qmax, qn,b + Ru)].  When 0 < qn,b < qmin, the unit should be in the progress of shutdown if qn,b < qn1,b and the range of qn.e is [max(qn,b  R1, 0), max(qn,b  R0, 0)]; if qn,b > qn1,b, the unit should be in the progress of startup and the range of qn.e is [min(qn,b + R0, max(qmin, qn,b + Ru)), min(qn,b + R1, max(qmin, qn,b + Ru))].  When qn.b = 0, the unit is offline. If Tc P T Dmin , the unit is allowed to startup and the range of qn,e is [min(R0, qmin), min(R1, qmin)] or 0, otherwise the unit is not allowed to startup and qn.e is 0. Usually, electricity price and unit output are reported as schedules by period. So unit output should be discretized. Output in one period, fuel cost, and environmental cost can be calculated by average output. In reality, the unit output should be as smooth as possible for load-follow and automatic generation control (AGC). Assume the output changes linearly from the beginning to the end in one period, the average output qn of period n is:

qn;e þ qn;b qn;e þ qn1;e ¼ 2 2

8 d > < X ¼ ðx1 ;x2 ;  ;xd Þ 2 R subject to g j ðXÞ 6 0; j ¼ 1; 2;. .. ;J > : hk ðXÞ ¼ 0; k ¼ 1; 2; .. .; K

ð10Þ

The output at the end of each period can be computed by multiobjective optimization algorithm.

i ¼ 1; 2;. .. ;Nobj

ð11Þ where fi(X) is objective functions; X is d-dimension decision-making vector; Nobj is the number of objective functions; gj(X) is inequality constraint function and hk(X) is equality constraint function. To evaluate the superiority or inferiority of solutions of MOP objectively, the following definitions are frequently used. For decision-making vector A and B,

ð9Þ

where Ru and Rd are the highest increasing and decreasing variation rate in normal operation, respectively; R1 and R0 are the upper limit and lower limit of variation rate while unit is in the process of start-up or shut-down. The multi-objective optimization model of benefit/risk/emission generation bidding in one dealing day is formed by integrating all objective functions and constraints mentioned above.

qn ¼

Taking a multi-objective minimization problem with a set of constraints as an example, the mathematical description of MOP can be described as:

(1) Pareto Dominance: A  B if and only if



fi ðAÞ 6 fi ðBÞ; 8i 2 f1; 2; . . . ; Nobj g fj ðAÞ < fj ðBÞ; 9j 2 f1; 2; . . . ; Nobj g

ð12Þ

(2) Pareto optimal or Pareto non-dominated: Solution A is Pareto optimal (Pareto non-dominated) solution if and only if

:9X 2 Rd : X  A

ð13Þ

Eqs. (12) and (13) should meet the constraints given in (11). The set of all Pareto optimal solutions is called Pareto optimal set. The values of objective functions corresponding to the Pareto optimal solution are called non-dominated objective vector. The area formed by all non-dominated objective vectors is called Pareto front. To solve MOP actually is to find Pareto optimal solutions as much as possible and make their corresponding objective vectors uniformly distribute on Pareto front. 3.2. Pareto non-dominated sorting and selection In recent years, multi-objective evolutionary algorithms (MOEAs) have become a research focus because they need not measure weighting relationship between each objective accurately and can find multiple non-dominated solutions with good distribution in a single run [17]. According to current existing research articles, from the perspective of algorithm frameworks, MOEAs can be divided into MOEAs based on non-dominated sorting(NSGA-II) [11], decomposition-based MOEAs(MOEA/D) [18], preference-based MOEAs [19], indicator-based MOEAs [20], memetic MOEAs [21], coevolutionary MOEAs [22], and so on. One of the most effective and commonly used MOEA is NSGA-II. It improves conventional non-dominated sorting genetic algorithm (NSGA) by adopting fast sorting of Pareto non-dominated solutions, elitism preservation and selection operator based on sorting grade of each solution and crowding distance. And its performance has been improved greatly. More details about NSGA-II procedure and Pareto non-dominated sorting are in [11]. At the same Pareto non-dominated sorting grade, the greater the crowding distance value of individual is, the looser the area where it situated is, and that means that the more valuable the individual in the area is and the more it should be kept in evolution process. The selection operation in NSGA-II works according to this principle. Suppose individual A and C are individuals before and after individual B, respectively, then the crowding distance (it

16

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

means the dispersion degree) Dc(B) of individual B in NSGA-II can be calculated by the equation below.

Dc ðBÞ ¼

Nobj X

jfi ðAÞ  fi ðCÞj

ð14Þ

i¼1

where fi(A) and fi(C) are the ith objective function values of individual A and C, respectively. The crowding distances of individuals on border are defined as infinite to assure them to be selected into next generation unconditionally. Though the equation is concise and easy to use, it still has limitations. In fact, the dispersion degree of individual B is not only related to the size of its neighborhood but also to its position in its neighborhood. So Eq. (14) would perhaps lead to elimination of some well-distributed individuals and preservation of some bad-distributed individuals in evolution [23]. The diversity and distribution of solutions would deteriorate along with the development of evolution. Thus the solutions cannot be convergent to Pareto front evenly and accurately. Therefore, as shown in Fig. 1, the equation for calculating the crowding distance of individual B between A and C can be improved as follows: Nobj X Dc ðBÞ ¼ ðjfi ðAÞ  fi ðCÞj  jfi ðBÞ  fi ðOÞjÞ i¼1

¼

Nobj X ðjfi ðAÞ  fi ðCÞj  0:5 þ min½jfi ðAÞ  fi ðBÞj; jfi ðBÞ

only a few lines of code and takes very few control parameters in comparison with most of evolutionary algorithms. For most of non-linear optimization problems, it is proved to be better than those algorithms such as genetic algorithm, particle swarm optimization algorithm, evolution strategy, adaptive simulated annealing [24–27]. Therefore, this paper proposes the non-dominated sorting differential evolution algorithm by substituting genetic operation in NSGA-II algorithm with differential evolution operation. 3.3. Differential evolution algorithm Just like GA, DEA also involves computational operations such as mutation, crossover and selection. The detailed description of DEA procedure can be found in [24]. In differential evolution, suppose the population size is Np, and each individual vector of the Gth generation is X i;G ¼ ðxi;1 ; xi;2 ; . . . ; xi;d Þ, the middle individual vector Yi,G+1 = (yi,1, yi,2, . . . , yi,d) can be produced by mutation through (16)–(18), and so on. Some typical DE mutation strategies are summarized in [25].

Y i;Gþ1 ¼ X r3 ;G þ F  ðX r1 ;G  X r2 ;G Þ

ð16Þ

Y i;Gþ1 ¼ X best;G þ F  ðX r1 ;G  X r2 ;G Þ

ð17Þ

Y i;Gþ1 ¼ X i;G þ F  ðX best;G  X i;G Þ þ F  ðX r1 ;G  X r2 ;G Þ

ð18Þ

i¼1

 fi ðCÞjÞ

ð15Þ

where O is the center point between individual A and C, fi(B) and fi(O) are values of the ith objective function of individual B and O, respectively. Eq. (15) can comprehensively reflect that the crowding distance of individual B is related not only to neighborhood size of the objective function (represented by |fi(A)  fi(C)|) but also to the distance between the individual and the center of neighborhood (represented by |fi (B)  fi(O)|). That is, |fi (A)  fi(C)| "or |fi (B)  fi(O))|;, Dc(B)". Moreover, since NSGA-II adopts the crossover and mutation mechanism of genetic algorithm, when we used it to solve the complicated MOP of this paper, some disadvantages of genetic algorithm such as instable convergence, slow and easy to premature exist in NGSA-II as well. Differential evolution (DE) is now regarded as one of the most powerful stochastic real-parameter optimization algorithms. Unlike traditional evolutionary algorithms, the mutation operator in DE is realized by the scaled differences of randomly selected and distinct population members. Therefore, no separate probability distribution has to be used, which makes the evolution completely self-organizing. Additionally, DE is very simple and easy to use as it requires

f2

A B

O C

A’

B’ O’

C’

Fig. 1. Crowding distance calculation.

f1

where i; r1 ; r 2 ; r 3 2 f1; 2; . . . ; N p g; r1, r2 and r3 are randomly chosen and r1 – r2 – r3; Xbest,G is the best individual vector of the Gth generation; mutation scaling factor F is a real constant from the interval [0, 2]. After analysis and comparison, we improved mutation strategy by introducing jitter variation into (17), as shown in (19).

Y i;Gþ1 ¼ X best;G þ ½F þ ð1  FÞ  rand  ðX r1 ;G  X r2 ;G Þ

ð19Þ

where rand is a uniformly random real number from the interval [0, 1]. However, the difference vector ðX r1 ;G  X r2 ;G Þ in (19) usually tends to be zero and then lose effectiveness because of small difference of each individual in later stage of the evolution. To avoid such situation to occur, we improved the mutation strategy by borrowing the polynomial mutation operator of MOEA/D-DRA described in [28], as shown in (20).  if Ra < 0:2  0:1  jxr1 ;j  xr2 ;j j xbest;j þ kj  ðaj  bj Þ; yi;j ¼ xbest;j þ ½F þ ð1  FÞ  rand  ðxr1 ;j  xr2 ;j Þ; otherwise ( 1 g þ1 if Rb < 0:5 ð2  Rb Þ  1; with kj ¼ 1 1  ð2  2  Rb Þgþ1 ; otherwise ð20Þ where j e {1, 2, . . . , d}; Ra and Rb are uniformly random number from the interval [0, 1]; the distribution index g is a control parameter; aj and bj are the lower and upper bounds of the jth decision variable, respectively. In this mutation strategy, a mutation probability 0:2  0:1  jxr1 ;j  xr2 ;j j is designed to make dynamic selection of DE mutation operator and polynomial mutation operator. If the element difference of an individual is less than 2, then the smaller this difference is, the greater the occurrence probability of polynomial mutation is. Such improvement can effectively increase the diversity of individuals and avoid evolution stagnation in later stage of evolution. Then, through crossover between objective individual vector Xi,G and middle individual vector Yi,G+1 by (21), a candidate individual of next generation Zi,G+1 is produced to maintain diversity of population.

17

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

( Z i;Gþ1 ¼ ðzi;1 ; zi;2 ;    ; zi;d Þ; zi;j ¼

xi;j ;

Rj > C R

ð21Þ

yi;j ; otherwise

In the above equation, if hi,j is zero, then hi,j ln hi,j = 0. The weight for each objective can be defined as:

1  Hj

where Rj is a uniformly random real number from the interval [0, 1]; crossover probability factor CR is a real constant from the interval [0, 1].

xj ¼

3.4. General best solution

And then, use Eq. (24) to calculate the general satisfaction degree of each solution in the Pareto optimal set.

Usually, only one solution will be put into practice, the decision makers have to choose a general best solution from the Pareto optimal set, therefore a multi-objective optimization decision-making is necessary in the end [29]. For better comparison, the satisfaction of each solution can be evaluated by fuzzy set theory. Actually, for a Pareto optimal solution, the closer the value of an objective function against the minimum, the higher its satisfaction degree corresponding to the objective function is, and vice versa. To simplify the process of assessing satisfaction degree for each solution, it can be regarded that satisfaction degree and value of objective function is approximately in inverse proportion. So we can use a linear fuzzy membership function, (21), to evaluate the satisfaction degree for each Pareto optimal solution corresponding to each objective [9].

si;j ¼

8 1; > > > max < fj

fi;j 6 fjmin fi;j

fjmax fjmin > > > : 0;

; fjmax > fi;j > fjmin

ð21Þ

fi;j P fjmax

ði ¼ 1; 2;    ; m; j ¼ 1; 2;    ; Nobj Þ where m and Nobj are the number of Pareto optimal solutions and objective functions respectively; fi,,j is the jth objective function value for the ith Pareto optimal solution; fjmax and fjmin are the maximum and minimum value of the jth objective function respectively; si,j is 1 means complete satisfaction to the jth objective function value while 0 means complete dissatisfaction. The linear membership shape is shown in Fig. 2. In information theory, entropy can be used to measure the amount of useful information provided by data [30]. The basic idea of entropy weighting method is that the smaller the difference degree of the data associated with an objective is, the larger its entropy is, in other words, the less the amount of information provided by the objective is, the smaller the weight of this objective is in general evaluation, and vice versa. This is consistent with the actual decision-making ideas. So we can use all si,j to constitute a decision-making information matrix (m  Nobj), and then use the decision-making information matrix to determine the weight of each objective in general evaluation by entropy weighting method. The entropy weighting method can effectively avoid subjectivity of weights assignment, so the evaluation results will be more realistic [29]. Based on the principle of entropy weighting method, the weight of each objective H can be defined as (22). m P

Hj ¼  i¼1

ðhi;j ln hi;j Þ

;

ln m

hi;j ¼ si;j =

m P

si;j

ð22Þ

i¼1

ð23Þ

N obj

Nobj 

P

Hj

j

Vi ¼

N obj X

xj si;j ; ði ¼ 1; 2;    ; mÞ

Finally, the Pareto optimal solution that has the maximum V is chosen to be the general best solution. 3.5. Multi-objective differential evolution algorithm procedure There are some researchers who have proposed to solve MOP using differential evolution algorithm [31]. However, in the nondominated sorting process of these algorithms, crowding distance of each individual is calculated by (14) just like NSGA-II. And when these algorithms are used in solving a more complex multi-objective problem just like the benefit/risk/emission comprehensive generation optimization investigated in this paper, it is difficult to get a complete, well distributed and accurate Pareto front. In addition, these algorithms did not give an extraction method of the general best solution, but just a Pareto optimal set, which restrict their practicability. Therefore, based on the improved Pareto non-dominated sorting and differential evolution algorithm mentioned above, an improved non-dominated sorting differential evolution algorithm (INSDEA) is proposed. The procedure of INSDEA is shown as Fig. 3, where Gmax is maximum generation, F is mutation scaling factor, CR is crossover probability factor. The initial parent population U0 is formed by m real-num-

Start Input unit parameters and data of forecasted electricity prices, and set Gmax, F and CR appropriately Improved Pareto non-dominated sorting

Initialization: evolutional generation g = 0. Make initial parent population matrix Ug (m 25) and initial child population Sg (empty)

Parent population update: Ug U g 1 Child population update: a) Tournament selection b) Differential evolution

Population mixture: Mg Ug Sg, then calculate objective function values of each individual in Mg

ðj ¼ 1; 2;       ; Nobj Þ g=g+1

Sg

No

S

g 1

g =Gmax ?

Yes General best solution selection: Integrated fuzzy set theory and entropy weighting method

s 1

0

ð24Þ

j¼1

f min

f max

Fig. 2. The linear membership shape.

f

End Fig. 3. Flowchart of the proposed INSDEA.

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

ber encoding individuals with 25 digits. The 1st–24th digits of each individual are output values at the end of each period qn.e (n = 1, 2, . . . , 24) which are randomly produced in sequence according to initial unit status of the first period and all constraints, and the 25th code is the bidding risk coefficient which is randomly produced between 0 and 1. The cyclic iterative process consists of four stages: Stage 1 (Population mixture): combine child population Sg and parent population Ug into a temporary population Mg. All repeated individuals are forced to be locally mutated to assure all the individuals in the population are different from each other. Calculate values of each objective function of all newborn individuals. Stage 2 (Pareto non-dominated sorting): According to Pareto non-dominated sorting strategy, find the Pareto non-dominated individual set in current population to compose Ps(1) by comparing objective function values of each individual. Remove all the individuals in Ps(1) from the current population. And then find new Pareto non-dominated individual set in remained population to compose Ps(2). Repeatedly to do so until all individuals are finished sorting. Then calculate crowding distances of each individual in each grade by (15). Stage 3 (Parent population update): fill Ps(1), Ps(2)... into an empty population in turn until the population scale would exceed m if further filling. And then fill individuals into Ps(i) in a descending order according to crowding distance till the population scale reach m. Thus form new parent population Ug+1. Stage 4 (Child population update): a. (Tournament selection): applying random tournament selection to generate optimal population from Ug+1 based on the principle that the smaller of the sorting grade the prior at different grades and the greater of crowding distance the prior at the same grade. The size of optimal population is generally half of the size of parent population. b. (Differential evolution): applying mechanism of mutation and crossover of DE analyzed above as (20) and (21) to produce new child population Sg+1. Stage 5 (General best solution selection): Finally, solutions in the parent population at the last iteration are regarded as the Pareto optimal set, and then the general best solution for this MOP can be extracted from the optimal set according to the method described in Section 3.4.

Table 1 Units parameters. a ($/MW2 h) 0.012

b ($/MW h) 10.91

c ($/h) 720.63

g ($/h) 56

h (rad/MW) 0.082

K0 ($) 0

K1 ($) 2910

s (h)

Cfa ($/h) 1295

D ($) 120

a (t/MW2 h)

b (t/MW h) 5.849e2

c (t/h)

6.213e4

8.510

n (t/h) 5e4

k (MW h) 0.019

4

qmin (MW)

qmax (MW)

160

350

TD min (h) 3

TU min (h) 4

Rd (MW/h) 120

Ru (MW/h) 70

R0 (MW/h) 60

R1 (MW/h) 150

80 Market margin price

Price ($/MWh)

18

60

Spinning reserve price

40 20 0 0

6

12

18

24

Hours Fig. 4. Forecasted electricity prices curves.

Table 2 r2n (variances of prices). Hour r2 ($/MW h)

1 0.20

2 0.44

3 0.09

4 0.26

5 0.29

6 0.55

Hour r2 ($/MW h)

7 1.99

8 1.98

9 1.61

10 1.62

11 0.08

12 1.01

Hour r2 ($/MW h)

13 1.01

14 0.93

15 0.68

16 0.41

17 0.60

18 1.49

Hour r2 ($/MW h)

19 0.25

20 1.43

21 0.47

22 0.89

23 0.23

24 0.23

4. Example and analysis In order to give a comprehensive illustration of the effectiveness of the proposed INSDEA in solving the multi-objective optimization problem of generation bidding and scheduling, a set of experiments on one unit in a generation plant are developed and performed in this section. The parameters of the unit are shown in Table 1. The initial unit status is: having been running for 3 h; the power output is 200 MW. The forecasted values of market margin electricity prices and spinning reserve prices in one dealing day are shown in Fig. 4 and the variances of electricity price in each period are shown in Table 2. All the experiments are done on a personal computer with Intel Core2 T4500 2.3-GHz processor, 2-GB 667-MHz DDR2 memory, and the programming environment for all operations and plots is Matlab 7.4. Aiming at solving two-objective optimization problems, one with profit and risk as its objectives and the other with profit and emissions as its objectives, the INSDEA proposed in this paper, NSDEA (without any improvement on the crowding mechanism) and NSGA-II are adopted. The parameters for the first two algorithms are specified as: Gmax is 300, m is 200, F is 0.3, g is 20 and

CR is 0.5, as for NSGA-II, Gmax is 600, m is 200, crossover probability is 0.95 and mutation probability is 0.05. Furthermore, some recent MOEAs such as SNOV-DS-MODE [32] and 2LB-MOPSO [33] are also applied in the same problem for comparison. After several tests, the parameters of these two algorithms were determined respectively as followed. For SNOV-DS- MODE, Gmax is 600, m is 300, P (scanning percentage) is 0.9, n_bins (the number of bins) is 200, and the other parameters are same as those of INSDEA; for 2LBMOPSO, Gmax is 600, m is 300, count (the number of iterations the particle fails to contribute a solution to the external archive) is 5 and n_bins is 10. Fig. 5 shows the Pareto front of optimization problem with profit and risk as its two objectives. Fig. 6 shows the Pareto front of optimization problem with profit and emissions as its objectives while bidding risk coefficient is 0.5. Shapes of the two kinds of Pareto fronts are different due to their different functional relationship. Fig. 5 indicates that bidding risk would increase with maximum expected profit without considering emissions. However, after optimizing, the output schedule can be adjusted flexibly

19

INSDEA

0.8

NSDEA

0.6

NSGA-II

0.4

INSDEA

0

20

105

110

115

120

125

130

135

Profit (1000$) 1 SNOV-DS-MODE

0.8

0.4 0.2 105

110

115

25 20 15 10 5 0

120

125

130

135

80

100

120

0

20

2LB-MOPSO

0.8

40

60

80

100

120

25 20 15 10 5 0

100

120

NSGA-II

0

20

40

60

80

Profit (1000$) Emission (ton)

1 0.6 0.4

25 20 15 10 5 0

0.2

SNOV-DS-MODE

0

20

40

60

80

100

120

100

120

Profit (1000$) 105

110

115

120

125

130

135

Profit (1000$) Fig. 5. The Pareto front of risk-profit.

Emission (ton)

Risk coefficient

60

NSDEA

Profit (1000$)

0

40

Profit (1000$)

0.6

0

Emission (ton)

0.2 0

Risk coefficient

25 20 15 10 5 0

Profit (1000$)

Emission (ton)

Risk coefficient

1

Emission (ton)

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

25 20 15 10 5 0

2LB-MOPSO

0

20

40

60

80

Profit (1000$) followed by the change of bidding tactics to realize maximum of profit and the sensitivity about bidding risk can be decreased dramatically. For those conservative generation companies, they can get high profit applying bidding tactic with low bidding risk. From Fig. 6, it is clear that emissions would increase with maximum expected profit. To evaluate the performance and efficiency of the algorithms employed in our cases, performance indicators (PIs) is used. At present, a series of nonparametric statistical tests have been developed for comparing the performance of evolutionary and swarm intelligence algorithms [34]. In this paper, the distribution and spread PIs are defined as the diversity metric (D) and maximum spread (Smax), which have been described in detail in [11,35], respectively. D is used to measure the nonuniformity in the distribution, and Smax is used to measure the distance between the boundary solutions in the Pareto-approximation set. However, since the optimal solution set for our multi-objective optimization problem is unknown, it is hard to use convergence metric [11] to evaluate accuracy PI. Therefore, a volume-based accuracy PI, also called Hypervolum (H) [36], is introduced. The basic idea of Hypervolum is that the larger the space that the solutions can dominate in the fitness space, the better the corresponding algorithm is. Table 3 shows the mean and variance of the accuracy, distribution and spread PIs obtained by INSDEA, NSDEA, NSGA-II, SNOV-DS-MODE and 2LB-MOPSO, respectively. All results are average of 10 runs, and calculations of all PIs are on normalization basis. From Table 3, we can find that INSDEA has the largest mean Smax and H and lowest mean D. Based on comprehensive analysis of Figs. 5 and 6 and Table 3, we can also find that INSDEA are superior to the other four algorithms in solutions distribution, searching speed, convergence stability, accuracy and integrality of Pareto front. Moreover, both SNOV-DS-MODE and 2LB-MOPSO

Fig. 6. The Pareto front of profit-emission (v = 0.5).

are very sensitive to parameters, especially P, n_bins and count. Those parameters have great impact on simulation results and the results are disappointment unless higher Gmax and m are set. So, even though the implementation procedure of SNOV-DS(Summation of Normalized Objective Values and Diversified Selection) is simpler than that of NS(Non-dominated Sorting and Selection), with higher Gmax and m, SNOV-DS-MODE has no obvious advantage in terms of convergence speed for solving the optimization problem of generation bidding and scheduling, and because without NS, it is no guarantee that all solutions are Pareto optimal solutions. Considering three optimization objectives of profit, bidding risk and emissions, Fig. 7 shows the Pareto fronts achieved by INSDEA and NSGA-II. Gmax = 2000 and m = 600 in both algorithms. Furthermore, the equatorial projections of emission and profit achieved from Fig. 7 are showed in Fig. 8 in order to more clearly show the difference. From these figures, Pareto front achieved by INSDEA is more complete and more smooth and better distributed than that achieved by NSGA-II. The equatorial projection of bidding risk and profit achieved by INSDEA is shown in Fig. 9. It indicates that high profit could be achieved with low bidding risk. From the up front of the equatorial projection, it shows that the maximum expected profit would increase with bidding risk without considering emissions as well. Generally, to avoid bidding failure, bidding risk coefficients should be controlled to less than 0.6 according to actual situation. Fig. 10 shows general best curves of the unit average output schedule achieved by INSDEA and NSGA-II, respectively. Corresponding

20

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

Table 3 Performance measures. Risk-profit

INSDEA ± NSDEA ± NSGA-II ± SNOV-DS-MODE ± 2LB-MOPSO ±

Profit-emission

H

D

Smax

H

D

Smax

0.12188 0.00216 0.11495 0.00539 0.09187 0.00641 0.11954 0.00547 0.12023 0.00596

0.43832 0.00566 0.74903 0.01608 0.44961 0.01964 0.53846 0.01839 0.48454 0.01742

1.01486 0.03645 1.01471 0.04296 0.99918 0.04905 1.00849 0.05182 1.01376 0.04943

0.21682 0.00320 0.20749 0.00615 0.18920 0.00483 0.21054 0.00692 0.21354 0.00659

0.35242 0.00514 0.43472 0.01258 0.56416 0.01621 0.81087 0.01446 0.66215 0.01096

1.45711 0.04281 1.29448 0.06649 1.18602 0.06837 1.23188 0.06813 1.26979 0.06686

150

Profit (1000$)

Profit (1000$)

150 100 50 0 -50 1

100 50 0

INSDEA 0.8

0.6

0.4

0.2

Risk coefficient

0

5

0

10

20

15

-50 0

25

0.2

0.4

0.6

0.8

1

Risk coefficient

Emission (ton)

Fig. 9. The equatorial projection of risk-profit by INSDEA.

400

100 50 0 1 0.8 0.6 0.4 0.2

Risk coefficient 0

NSGA-II

5

0

10

20

15

25

Emission (ton)

Output (MW)

Profit (1000$)

150

300 200 100 0

Output by INSDEA Output by NSGA-II 0

6

12

18

24

Hours

Fig. 7. The Pareto front of three objectives.

Fig. 10. General best curves of the unit average output schedule.

Profit (1000$)

150 Table 4 Generating and bidding schedule.

100 50 0

INSDEA

-50 0

5

10

15

20

25

Emission (ton)

Profit (1000$)

150 100

Hour Output (MW) Price ($/MW h) Hour Output (MW) Price ($/MW h) Hour Output (MW) Price ($/MW h) Hour Output (MW) Price ($/MW h)

1 247.89 25.59 7 167.85 37.60 13 190.62 33.51 19 350 55.07

2 217.36 20.56 8 237.85 43.79 14 167.85 28.81 20 289.20 47.46

3 205.32 23.55 9 307.85 44.50 15 237.85 26.46 21 350 42.70

4 174.17 22.14 10 333.09 47.07 16 286.93 25.42 22 230 35.29

5 169.25 22.35 11 350 48.29 17 240.47 35.20 23 199.63 28.89

6 163.20 24.70 12 230 41.49 18 310.47 53.26 24 138.36 23.98

50 0 -50

NSGA-II

0

5

10

15

20

Emission (ton) Fig. 8. The equatorial projection of emission-profit.

25

with the former schedule, the generation bidding risk coefficient is 0.337 while generation expected profit of the day is 91,458 $ and total emissions 10.024 ton. Corresponding with the latter schedule, the generation bidding risk coefficient is 0.487 while generation expected profit of the day is 87,881 $ and total emissions 10.541 ton. It is clear that the generation company can achieve higher

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

profit with lower bidding risk and lower emissions after optimization by INSDEA. The bidding price of each period can be determined according to each coefficient risk by (5). The generating and bidding schedule achieved by INSDEA is shown in Table 4. For those generation companies which have several units, the first step is to calculate optimal output scheme of each unit in each period by INSDEA. Then add them all together by period to get the final generating and bidding plans of the whole company. Moreover, in electric power market, it is permitted that the generation companies provide more than one output-bidding schedule. That means the generation companies can choose several output-bidding schedules with different bidding risk coefficients from Pareto optimal solutions to report to the dispatch center, thus to increase space of decision-making and reduce bidding risk. Since all unit outputs are determined considering bidding risk and based on electricity market price fluctuation, the generating and bidding schedules cannot only ensure units to be chosen to generate electricity in the dealing day but also achieve ideal profit.

5. Conclusions This paper has proposed a hybrid multi-objective evolutionary algorithm (INSDEA) to handle the benefit/risk/emission generation bidding and scheduling optimization problem. In INSDEA, a new calculation method of crowding distance which takes into account the distance between the individual and the center of its neighborhood is designed. It can effectively avoid the elimination of welldistributed individuals and the preservation of bad-distributed individuals during the stage of Pareto non-dominated sorting. The diversity and distribution of solutions with new crowding distance calculation is much better. On the other hand, an improved differential evolution operation with the introduction of jitter variation in mutation strategy is substituted for the genetic operation adopted by NSGA-II. It can effectively overcome premature convergence and search bias problems and improve the stability and convergence speed in the process of the optimization. Besides focusing on an organic integration on Pareto non-dominated sorting and differential evolution algorithm, fuzzy set theory is used to evaluate the satisfaction degree for each Pareto optimal solution corresponding to each objective, entropy weighting method is employed to determine the weight of each objective to extract the general best solution objectively, and thus to form a reasonable method for extracting the general best solution. Several optimization runs of the proposed approach have been carried out on a thermal unit. The results show that the proposed approach is efficient for solving multi-objective optimization problem of generation bidding for maximizing profit and minimizing both risk and emissions. It can obtain Pareto-optimal set easily and quickly, and the non-dominated solutions in the obtained Pareto-optimal set are well distributed and have satisfactory diversity characteristics. In addition, the general best scheme integrated generation bidding and scheduling can be selected automatically. The proposed method has been applied in the Decision-making Support System for Energy-saving Generating and Power Market Bidding developed by Electric Power Research Institute of Jiangxi Province. In this system, INSDEA is mainly employed in guiding the regional generation units in Jiangxi to achieve environmental/economic load dispatch at the early development stage of China electric power market. The performance is quite great in energy saving and emission reduction. Furthermore, in the future, the system will allow all thermal generating units in the region using this algorithm to achieve benefit/risk/emission generation bidding and scheduling optimization at the maturity stage of electric power market. We are going to broaden its application followed by the development of power market reform in China.

21

Acknowledgments This work was supported by the National Natural Science Foundation of China (51167005, 60964004), the Natural Science Foundation of Jiangxi Province (2009GZS0016) and the Science Research Founds of Jiangxi education department (GJJ11114). References [1] Catalão JPS, Pousinho HMI, Mendes VMF. Short-term electricity prices forecasting in a competitive market by a hybrid intelligent approach. Energy Convers Manage 2011;52:1061–5. [2] Arroyo JM, Conejo AJ. Optimal response of a thermal unit to an electricity spot market. IEEE Trans Power Syst. 2000;15:1098–104. [3] Ni E, Luh PB, Rourke S. Optimal integrated generation bidding and scheduling with risk management under a deregulated power market. IEEE Trans Power Syst 2004;19:600–9. [4] Boonchuay C, Ongsakul W. Optimal risky bidding strategy for a generating company by self-organising hierarchical particle swarm optimization. Energy Convers Manage 2011;52:1047–53. [5] Vahidinasab V, Jadid S. Normal boundary intersection method for suppliers’ strategic bidding in electricity markets: an environmental/economic approach. Energy Convers Manage 2010;51:1111–9. [6] Cai J, Ma X, Li Q, Li L, Peng H. A multi-objective chaotic particle swarm optimization for environmental/economic dispatch. Energy Convers Manage 2009;50:1318–25. [7] Spear SE. The electricity market game. J Econ Theory 2003;109:300–23. [8] Shayanfar H A, Saliminia Lahiji A, Aghaei J, Rabiee A. Generation expansion planning in pool market: a hybrid modified game theory and improved genetic algorithm. Energy Convers Manage 2009;50:1149–56. [9] Abido MA. Multiobjective evolutionary algorithms for electric power dispatch problem. IEEE Trans Evol Comput 2006;10:315–29. [10] Talaq JH, Ferial, EI-Hawary ME. A summary of environmental/economic dispatch algorithms. IEEE Trans Power Syst 1994;9:1508–16. [11] Deb K, Pratap A, Agarwal S. A fast and elitist multi-objective genetic algorithm: NSGA-II. IEEE Trans Evol Comput 2002;6:182–97. [12] Leticia C, Susana E, Carlos A. A particle swarm optimizer for multi-objective optimization. J Comput Sci Technol 2005;5:204–10. [13] Soleymani S, Ranjbar AM, Shirani AR. New approach to bidding strategies of generating companies in day ahead energy market. Energy Convers Manage 2008;49:1493–9. [14] Rahimiyan M, Mashhadi HR. Risk analysis of bidding strategies in an electricity pay as bid auction: a new theorem. Energy Convers Manage 2007;48:131–7. [15] Wu QH, Guo J. Optimal bidding strategies in electricity markets using reinforcement learning. Elect Power Compon Syst 2004;32:175–92. [16] Bajpai P, Singh SN. Fuzzy adaptive particle swarm optimization for bidding strategy in uniform price spot market. IEEE Trans Power Syst 2007;22:215–60. [17] Zhou A, Qu B-Y, Li H, Zhao S-Z, Suganthan PN, Zhang Q. Multiobjective evolutionary algorithms: a survey of the State-of-the-art. Swarm Evol Comput 2011;1:32–49. [18] Zhang Q, Li H. MOEA/D: a multiobjective evolutionary algorithm based on decomposition. IEEE Trans Evol Comput 2007;11:712–31. [19] Thiele L, Miettinen K, Korhonen PJ, Luque JM. A preference-based evolutionary algorithm for multi-objective optimization. Evol Comput 2009;17:411–36. [20] Basseur M, Zitzler E. Handling uncertainty in indicator-based multiobjective optimization. Int J Comput Intel Res 2006;2:255–72. [21] Lara A, Sanchez G, Coello Coello CA, Schutze O. HCS: a new localsearch strategy for memetic multiobjective evolutionary algorithms. IEEE Trans Evol Comput 2010;14:112–32. [22] Goh CK, Tan KC. A competitive-cooperative coevolutionary paradigm for dynamic multiobjective optimization. IEEE Trans Evol Comput 2009;13:103–27. [23] Peng C, Sun H, Guo J. Multi-objective optimal PMU placement using a nondominated sorting differential evolution algorithm. Int J Elect Power Energy Syst 2010;32:886–92. [24] Storn R, Price K. Differential evolution-a simple and efficient heuristic for global optimization over continuous space. Global Optim 1997;11:341–59. [25] Tasoulis DK, Pavlidis NG, Plagianakos VP, Vrahatis MN. Parallel differential evolution. In: Proc. 2004 IEEE Cong Evolutionary Computation, CEC’04, Portland, OR; 2004. [26] Coelho LS, Mariani VC. Improved differential evolution algorithms for handling economic dispatch optimization with generator constraints. Energy Convers Manage 2007;48:1631–9. [27] Paterlini S, Krink T. Differential evolution and particle swarm optimisation in partitional clustering. Comput Stat Data Anal 2006;50:1220–47. [28] Zhang Q, Liu W, Li H. The performance of a new version of MOEA/D on CEC09 unconstrained MOP test instances. In: Proc Congress on Evolutionary Computation, 2009, Norway; 2009. [29] Li X. Study of multi-objective optimization and multi-attribute decisionmaking for economic and environmental power dispatch. Elect Power Syst Res 2009;79:789–95. [30] Shannon CE, Weaver W. The mathematical theory of communication. Urbana: The University of Illinois Press; 1949.

22

C. Peng et al. / Energy Conversion and Management 57 (2012) 13–22

[31] Mezura-Montes E, Reyes-Sierra M, Coello Coello CA. Multi-objective optimization using differential evolution: a survey of the State-of-the-Art. In: Chakraborty UK, editor. Advances in differential evolution. Berlin: Springer; 2008. p. 173–96. [32] Qu BY, Suganthan PN. Multi-objective evolutionary algorithms based on the summation of normalized objectives and diversified selection. Inform Sci 2010;180:3170–81. [33] Zhao SZ, Suganthan PN. Two-lbests based multi-objective particle swarm optimizer. Eng Optim 2011;43:1–17.

[34] Derrac J, García S, Molina D, Herrera F. A practical tutorial on the use of nonparametric statistical tests as a methodology for comparing evolutionary and swarm intelligence algorithms. Swarm Evol Comput 2011;1:3–18. [35] Zitzler E, Deb K, Thiele L. Comparison of multiobjective evolutionary algorithms: empirical results. Evol Comput 2000;8:173–95. [36] Zitzler E, Thiele L. Multiobjective evolutionary algorithms: a comparative case study and the strength Pareto approach. IEEE Trans Evol Comput 1999;3:257–71.