Daily hydrothermal scheduling with economic emission using simulated annealing technique based multi-objective cultural differential evolution approach

Energy 50 (2013) 24e37

Contents lists available at SciVerse ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Daily hydrothermal scheduling with economic emission using simulated annealing technique based multi-objective cultural differential evolution approach Huifeng Zhang, Jianzhong Zhou*, Na Fang, Rui Zhang, Yongchuan Zhang School of Hydropower and Information Engineering, Huazhong University of Science and Technology, Wuhan 430074, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 4 March 2012 Received in revised form 13 November 2012 Accepted 1 December 2012 Available online 5 January 2013

This paper develops a simulated annealing based multi-objective cultural differential evolution (SAMOCDE) for solving daily hydrothermal operation scheduling with economic emission (SHOSEE) while considering water transport delay and transmission loss of system load, and it is formulated for optimizing thermal cost and emission issue with complex non-convex constraints. The simulated annealing technique is integrated into the multi-objective differential computation model motivated by redesigned cultural knowledge structures, and it can be used to avoid the premature convergence by properly controlling the population space evolution. For testifying the convergence ability and diversity distribution of the proposed algorithm, various benchmark functions are used to verify the stable efficiency of complex constraint-handling and high-dimensional search ability. The implementation of the proposed SA-MOCDE on SHOSEE problem proves that more desirable fuel cost and emission cost are obtained by SA-MOCDE in comparison to other established method recently, and the analysis on the running performance of convergence and diversity metric also reveals that SA-MOCDE can be a promising alternative for solving SHOSEE problem.
2012 Elsevier Ltd. All rights reserved.

Keywords: Hydrothermal system Transmission loss Simulated annealing technique Cultural knowledge structure Premature convergence Convergence and diversity metric

1. Introduction The daily hydrothermal optimal scheduling plays an important role in inter-connected power system operation. The main goal of hydrothermal scheduling problem is to minimize the total fuel cost of thermal plants while properly dispatching the power generation of hydro plants and thermal units with satisfying various types of highly nonlinear coupled constraints, which increases the difficulty of optimization search process. In the past few years, various mathematical methods have been used to solve this problem, such as dynamic programming (DP) [1,2], non-linear programming (NP) [3], lagrangian relaxation [4] and network flow method [5]. They can tackle with the hydrothermal problem to some extent, but their optimization process mainly depends on the initial solution and easily falls into local optima. Then different heuristic problem-solving techniques, such as genetic algorithm [6,7], evolutionary programming [8], particle swarm optimization [9] and artificial immune system [10], have

* Corresponding author. Tel.: þ86 (0)2 78 7543 127. E-mail addresses: [email protected] mail.hust.edu.cn (J. Zhou).

(H.

Zhang),

0360-5442/$ e see front matter
2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.energy.2012.12.001

jz.zhou@

been applied on the hydrothermal operation system. Those methods can overcome the drawbacks that mathematical methods depend much on the convexity and continuity requirements of operation functions, and they can also have the global convergence ability with relative low computation effort. Since more concerns on environmental pollution are taken into consideration, the cheapest fuel cost and stable power supply cannot satisfy the society demand. Generally, several measurements can be taken to reduce the emission issue of gaseous pollutants [11e13], such as installation of cleaning equipment, replacement of the aged fuelburners, and emission dispatching. However, the first two choices need replacement of new equipments or modification of those existing ones, which can be only considered as long-term choice. Since conventional pure scheduling no longer satisfies this requirement, then economic scheduling can be extended to a multi-objective optimization problem (MOP) with considering both fuel cost and emission issue. Then, several weighted and constrained approaches were proposed to solve this multi-objective optimization problem by converting it into single-objective optimization problem [14], such as the interactive fuzzy satisfying method [15], goal-attainment method [16] and quantum-behaved particle swarm optimization method [17]. However, it lacks flexibility and engineering practicability at different degrees, generally there are two drawbacks: (1) the weights of those

H. Zhang et al. / Energy 50 (2013) 24e37

objectives are difficult to obtain. (2) Only one solution is obtained by each optimization, once the trade-off weights are changed, the whole optimization will be set up again. Afterwards, some evolutionary algorithms (EAs) have been extended to multi-objective evolutionary algorithms (MOEAs) for solving this problem. Among those so-called MOEAs, the most popular ones are Deb’s non-dominated sorting genetic algorithm (NSGA-II) [18], Zitzler’s enhanced strength pareto evolutionary algorithm (SPEA-II) [19], multi-objective particle swarm optimization (MOPSO) [20,21] and multi-objective differential evolutionary algorithm (MODE) [22e24]. All these approaches optimize the hydrothermal problem in single simulation run, and produce a set of nondominated solutions for decision making. Furthermore, those MOEAs provide the global search for various objectives, no matter whether they are non-differentiable or discontinuous. In this paper, a developed simulated annealing based multi-objective cultural differential evolution (SA-MOCDE) is proposed to solve the SHOSEE problem. Differential evolution (DE) is a simple yet powerful population based evolutionary algorithm [25], it provides an effective heuristic method for optimizing those particularly non-linear and nondifferentiable objectives. Recently it has been used to deal with those MOPs due to its population-based methodology and independence of problem representation, such as Pareto differential evolution (PDE) algorithm [26], Pareto differential evolution approach (PDEA) [27], differential evolution for multi-objective optimization (DEMO) [23] and adaptive differential evolution algorithm (ADEA) [24] are proposed. Though they have great ability of obtaining those optimal solutions, premature convergence is usually suffered at different degrees especially when it converges to certain optima. To solve the SHOSEE problem efficiently, an improved multiobjective differential evolutionary algorithm based on culture belief model is proposed to tackle with this problem. Culture algorithm provides the computational model for multi-objective differential evolutionary algorithm, all the evolutionary operators are integrated into the population space, and the obtained population experience is stored in the belief space during the population space evolution. Furthermore, in order to avoid the premature problem, simulated annealing technique is embedded into the DE operator to increase the population diversity especially when it converges to local optima. The efficiency of SA-MOCDE is proved with the application on the SHOSEE problem, and the obtained desirable results reveal that the proposed SA-MOCDE can be a promising alternative for solving SHOSEE problem. This paper is organized as follows: The problem formulation of SHOSEE is introduced in Section 2, some basic definitions about multiobjective optimization is presented in Section 3, the proposed SAMOCDE and its implementation are presented in Section 4 and Section 5, the application of SA-MOCDE on the hydrothermal system is presented in Section 6, and the research results are concluded in Section 7.

segments of quadratic functions of the active power output of the generator, and the fuel cost can be described as: 2 fit ðPsit Þ ¼ ai þ bi *Psit þ ci *Psit

2.1. Economic scheduling In the economic scheduling model, the major problem of hydrothermal scheduling is to minimize the total fuel cost with considering various hydrothermal constraints. Generally, the fuel cost curve for any thermal generating unit can be represented by

(1)

where Psit is the thermal output generated by the i-th thermal unit at t-th time period, and ai, bi, ci are the fuel cost coefficients of i-th thermal unit. However, since thermal power plants consist of generating units with multi-valve steam turbines in order to incorporate flexible operational facilities in reality, the above formulation (1) of fuel cost needs to be modified for more practicable and more accuracy. The effect of valve-point effect loading may be considered by adding a sinusoidal function to the quadratic cost function described above, and the modified fuel cost function can be presented as follows:


o
n 

2 min fit0 ðPsit Þ ¼ ai þ bi *Psit þ ci *Psit þ
di *sin ei * Psi  Psit


(2)

where Psit is the thermal output generated by the i-th thermal unit min is the minimum output of i-th thermal unit, at t-th time period, Psi and ai, bi, ci, di, ei are the fuel cost coefficients of i-th thermal unit. For a given hydrothermal system, the total fuel cost is presented as the summation of all thermal units during the entire scheduling period, which is expressed as:

Min F1 ¼

XNs XT i¼1

f 0 ðP Þ t¼1 it sit

(3)

where Ns denotes the number of thermal units in the hydrothermal system, T is the length of entire operation period. 2.2. Emission scheduling Due to the increase concern over the environmental consideration, society demands adequate and secure electricity not only at cheapest price but also at the minimum pollution level. The main task of emission scheduling problem is to minimize the thermal emission at minimum level, the amount of emission pollutants mainly depends on the generated thermal output. Here, the emission of nitrogen oxides (NOx) is considered as the major effect from the viewpoint of environment conservation, and the emission amount of one thermal unit during each time period can be formulated as follows: 2 eit ðPsit Þ ¼ ai þ bi *Pit þ gi *Psit þ hi *expðdi *Psit Þ

(4)

where ai, bi, gi, hi, di are the emission coefficients of the i-th thermal unit, and eit(Psit)denotes the emission rate of the i-th thermal unit at t-th time period. The emission scheduling problem can be described as the attempt to minimize the total emission amount of all the thermal units in the hydrothermal system during the entire operation period, which can be formulated as:

2. Problem formulation of SHOSEE The SHOSEE is formulated as a multi-objective optimization problem for minimizing the thermal cost and economic emission synchronously by utilizing water resource at each time period while satisfying the hydraulic relationships and various equality and inequality constraints [28].

25

Min F2 ¼

XNs

XT t¼1

i¼1

½eit ðPsit Þ

(5)

2.3. Constraints (1) System load balance limits The total active power generation must keep consistence with the summation of system load and transmission loss at each time interval over the scheduling horizon, which is generally formulated as:

XNs

P i¼1 sit

þ

XNh

P j¼1 hjt

¼ PDt þ PLt

(6)

26

H. Zhang et al. / Energy 50 (2013) 24e37

where Phjt is the output generated by the j-th hydro-plant at t-th time interval, PDt is the system load at t-th time interval, PLt is the transmission loss at t-th time interval, and Nh is the number of hydro plants in the hydrothermal system. The transmission loss is calculated in the form of quadratic function with B loss coefficients as follows:

PLt ¼

Ns X Ns X

Psit Bij Psjt þ

i¼1 j¼1

Ns X

B0i Psit þ B00

(7)

i¼1

where Bij, B0i, B00 are the B loss coefficients, Bij is the element of matrix B, and all these coefficients are given in reference [29]. Conventionally, the hydro power is closely related to reservoir discharge as well as net water head, and the net water head can also be expressed as the function of reservoir storage. Then the generated hydro power of any hydro plant in each time period can be formulated in terms of water discharge and reservoir storage, which is described as follows: 2 2 Phit ¼ C1i *Vhit þ C2i *Qhit þ C3i *Vhit *Qhit þ C4i *Vhit þ C5i *Qhit þ C6i

" Vhit ¼ Vhi;t1 þ Ii;t  Qhit  Si;t þ i ¼ 1; 2; :::; Nh ; t ¼

(2) Generation limit constraints

k¼1

Qhk;tsk þ Sk;tsk

1; 2; :::; T; 0  ski

i

i



#

Dt

ð12Þ

T

where Ii,t is water inflow of the i-th hydro plant at t-th time period, Si,t is spillage volume of the i-th hydro plant at t-th time period, Qhk;tsi and Sk;tsi represent the delayed discharge rate and spillage water of k-th hydro plant at t-th time period with the transmission delay ski between the i-th hydro plant and the k-th upstream hydro plant, Dt represents the time period length. (6) Initial and final reservoir storage limits

Vhi0 ¼ Vhi;begin ; VhiT ¼ Vhi;end

ði ¼ 1; 2; :::; Nh Þ

(13)

where Vhi,begin is the initial storage volume of the i-th hydro plant, Vhi,end is the final storage volume of the i-th hydro plant after T scheduling periods.

(9)

3. Multi-objective optimization problem In the tives are generally described

real world, non-commeasurable and conflictive objecoften involved in the optimization problems, and the multi-objective optimization problem can be as follows:

where f1,f2,.,fm are the objectives of MOP, and there are at least two conflictive objectives among them, m is the number of objectives, gj  0, hk ¼ 0 represent inequality and equality

MinFðx1 ; x2 ; .; xn Þ ¼ fmin f1 ; min f2 ; .min fi .; min fm g gj ðx1 ; x2 ; .; xn Þ  0 ði ¼ 1; 2; .; m; j ¼ 1; 2; .; S; k ¼ 1; 2; .; KÞ hk ðx1 ; x2 ; .; xn Þ ¼ 0

subject to

min ; P max represent the minimum and maximum power where Phi hi

generation limits of the i-th hydro plant, Psjmin ; Psjmax represent the minimum and maximum power generation limits of the j-th thermal unit. (3) Reservoir storage limits

min max Vhi  Vhit  Vhi

Ni  P

ð8Þ

where C1i, C2i, C3i, C4i, C5i, C6i present the power generation coefficients of the i-th hydro plant, Vhit and Qhit denote the reservoir storage and water discharge of the i-th hydro plant at t-th time period.

min  P max Phi hit  Phi min max Psj  Psjt  Psj

(5) Water dynamic balance constraint

(10)

(14)

constraints in MOP, S is the number of inequality constraints, K is the number of equality constraints. Since the objective functions are non-commensurable and conflicted with each other, it is impossible to find single optimal solution along with the population evolution. Instead of single solution, a set of non-dominated solutions are generated at each generation, and it doesn’t exist that a solution is completely superior or inferior to the other ones, which can be properly described by the following dominance relationship [30]. If there are two feasible solutions x1, x2, and x1 dominates x2 (x1 _x2 ) when the following conditions are satisfied:

min ; V max represent the minimum and maximum storage where Vhi hi limits of the i-th hydro-plant.

ci ¼ 1; 2; .; m; fi ðx1 Þ  fi ðx2 Þ

(4) Reservoir discharge limits

min max Qhi  Qhit  Qhi

(11)

min ; Q max represent the minimum and maximum water where Qhi hi discharge of the i-th hydro-plant.

(15)

where x1 dominates x2 (or x1 has higher dominance order than x2). If the formulation (15) is not always satisfied, both x1, x2 have same dominance order. If no solution in the current population dominates x1, then it is defined as pareto solution or non-dominated solution, and the objective value of pareto solution forms the pareto front of this MOP.

H. Zhang et al. / Energy 50 (2013) 24e37

4. The simulated annealing based multi-objective cultural differential evolutionary algorithm for solving SHOSEE problem The culture algorithm (CA) was proposed by Reynolds [31], since it was efficient for dealing with single objective optimization problem, it has been extended to solve multi-objective optimization problem. In this paper, multi-objective differential evolutionary algorithm is integrated into the computation model provided by CA. Furthermore, in order to avoid the premature problem, the simulated annealing technique is also embedded into the population space evolution, and the obtained population experience is stored in the belief space to guide population space evolution. And the efficiency of the proposed SA-MOCDE is proved by properly dealing with various benchmark test problems.

4.1. Overview of culture algorithm (CA) Some social scientists have suggested that culture could be transmitted within population as an inheritance mechanism. Then Reynolds proposed a cultural algorithm model [31], in which culture evolution is taken as an inheritance mechanism that operates at two levels: the micro-evolutionary and the macroevolutionary levels. At the micro-evolutionary level, the individuals in the population are presented as the behavior traits, which can be inherited from generation to generation by some motivated operators. While at the macro-evolutionary level, the individual experience of the micro-evolutionary level is presented, it can be merged to form the population “mappa” with motivated operators, and the individual experience and behavior trait can be properly converted to each other by some communication links. Corresponding with above two evolutionary levels, there are two main evolutionary spaces in the culture evolution: population space and belief space, which can be seen in Fig. 1. The population space takes charge of storing the evolutionary population and evolutionary operators, while the belief space is mainly in charge of storing the population experience or knowledge information obtained by population space evolution. Moreover, there are some communication protocols within/ between belief space and population space, such as generate (), objective (), select (), accept (), update () and influence (). Evolutionary population is produced by generate () function, the objective values are calculated by objective () function, and the optimal individual is selected by select () function. If the individuals are elite

Update ()

Belief space

Accept ()

Select ()

Influence ()

Population space

Generate ()

Fig. 1. The framework of culture algorithm.

Objective ()

27

ones, accept the individual experience by accept () function, then the belief space can be updated by update () function, and in turn, the newly obtained population experience can guide the population space evolution by the influence () function, and generates the individuals for next generation. 4.2. Modified population space evolution This paper modifies population space evolution mainly in following two aspects: elite retention mechanism and improved crossover rate. The elite retention mechanism retains the elite individual as the evolution proceeds, which can promote the population convergence further. And the crossover rate is developed to an adaptive parameter, which can properly control the population convergence as the population diversity changes. 4.2.1. Elite retention mechanism In the MOEAs, all the non-dominated individuals are arranged to be stored into the archive set. While the size of archive set is limited, and when the amount of non-dominated individuals exceeds this given size, a retention mechanism is needed to keep the elitism individuals [32]. 4.2.2. Modified crossover rate In differential evolution, the crossover rate is a constant, which never satisfies the population evolution requirement due to the changeable population diversity and convergence speed [33]. Here, an adaptive crossover rate is developed to control the population evolution. To accelerate the convergence speed [24], the modified crossover rate can be expressed as follows:

Pc ¼ Pc0 $2e

ð1GenNum=ðGenNowþ1ÞÞ

(16)

where GenNow represents the current generation number and Pc0 denotes the initial crossover rate when the generation number GenNow ¼ 0. As the generation number GenNow increases, the crossover rate also increases, which means that the diversity of the evolutionary population is decreased, when GenNow reaches the GenNum, the diversity of evolutionary population decreases to the nadir while the convergence speed increases to the climax. 4.3. Redefined knowledge structure of belief space Culture belief is mainly obtained from population experience along with population space evolution, and it contains a set of schemata, which are also known as cultural knowledge structures. Saleem has proposed some knowledge structures [34], the most famous structures of which are situational knowledge, normative knowledge, historic knowledge and domain knowledge structures, etc. In this paper, the situational knowledge, normative knowledge and historic knowledge structures are presented as computational models for the proposed algorithm and these computation models can be redesigned properly to solve the SHOSEE problem. 4.3.1. Situation knowledge structure The situational knowledge structure mainly takes charge of storing elite population experience as population evolution proceeds, and this population experience also guides the population space evolution, which can promote the population convergence by influence () function. In the multi-objective differential evolution, the obtained non-dominated individuals are properly stored in the situational knowledge structure. In the mutation operator, two individuals are randomly selected from situational knowledge structure, can lead the current individual to generate elite individual with population evolution.

28

H. Zhang et al. / Energy 50 (2013) 24e37

Situational knowledge is the basis of other knowledge, those non-dominated population experience in current population space is added into situational knowledge structure by accept () function. After population space evolution, the feasible individual experience that has higher dominance order is stored in the situational knowledge structure with dominance relationship introduced in Section 3. However, the situational knowledge structure has certain size, if the amount of those non-dominated individuals exceeds this size, and a truncation strategy is needed to retain these elite individuals, which has been introduced in Section 4.2.1.

  P DTr ¼ eDTr =Tr

According to formulation (20), when population evolution is close to convergence, the accept probability will become large, and the population diversity can also increase. In this paper, SA is embedded into situational knowledge for preventing local convergence, and the situational knowledge structure can be seen in Table 1. Since the accept probability of simulated annealing technique can affect the selection mechanism of DE, the selection operator in Section 4.1 can be modified as:

( 4.3.2. Normative knowledge structure The interval information of decision variables is mainly stored in the normative knowledge structure, which can be properly interpreted as legal restraint in human society. Since hydrothermal operation is often presented as a complex-constrained problem, normative knowledge structure can be properly designed to deal with those non-linear constraints here. In the normative knowledge structure, the upper and lower bounds of decision variables are properly stored as constraint intervals, and we adjust the individual experience mainly by moving those variables into these intervals. For properly controlling the violation information of function value, if the function value Fi violates its constraint limits, the returned information can be obtained by formulation (17):

fi ¼ ðFi  Fimin Þ=ðFimax  Fimin Þ

(17)

where Fimin, Fimax are the constraint limits of Fi. When it returns fi < 0, it means Fi < Fimin, when it returns fi > 1, then Fi > Fimax, and if 0  fi  1 is obtained, it means that Fi is under control. However, those decision variables cannot avoid violating constraint limits along with the constraint-handling process. Then in the situational knowledge structure we force these variables to move into the feasible domain by formulation (18) as follows:

x’i

8 < li ; if xi < li ¼ x ; if li  xi  ui : i ui ; if xi > ui

(18)

4.3.3. Historic knowledge structure The historic knowledge was originally proposed to find patterns in the environmental changes, and then used to solve dynamic objective functions [35]. Since all DEs suffer premature problem at different degrees, the extended MODE is also not an exception. Here, we redesign the historic knowledge structure with a simulated annealing technique to overcome this premature problem by properly controlling the population convergence state. Simulated annealing (SA) was originated from statistic mechanic, and it was used by Metropolis to simulate a collection of atoms in equilibrium at a given temperature [36], and then was first proposed for solving combinatorial optimization problem by Kirkpatrick [37]. The principle of SA is simple: Given a temperature parameter Tr, compare the objective value of current individual xg and the generated individual ugþ1, the total deviation of these objective values is defined as:

DTr ¼

 
Xm

ðjÞ  gþ1  ðjÞ g
u  f xi

f T T i r j¼1 r

(19)

where fTr($) is the j-th objective value under temperature Tr (Tr > 0), r (r ¼ 1,2,.,R) represents the iteration number as temperature changes, and Trþ1 ¼ aTr (0
(20)

xgþ1 i

¼

gþ1

ui

;

g

xi ;

if



gþ1

ui

g

_xi



    or P DTr > rnd

else

(21)

where rnd() is the random parameter in (0,1). Before utilizing this selection mechanism, check the convergence degree every 20  xgi khε (ε is the permitted accuracy) is found, generations, if kxgþ1 i let the modified selection operator replace the previous selection operator, set the temperature T0 and take the simulated annealing technique introduced above until the population evolution is ended. 4.4. Simulation of the proposed SA-MOCDE on benchmark problems For better testing the performance of SA-MOCDE, some benchmark problems are used to verify the optimization efficiency both on convergence and diversity property by Deb’s convergence and diversity metrics [18]. In addition, the proposed SA-MOCDE is also taken into comparison with other alternatives established recently on these benchmark problems. Firstly, ZDTs are taken for verifying the diversity distribution and convergence ability of SA-MOCDE, and these test functions are presented in Table 2. Here, ZDT1, ZDT3, ZDT4 and ZDT6 are used to testify the search ability mainly by whether they are continuous, uniform or convex. Main parameters of SA-MOCDE are set as follows: Initial crossover rate is set to 0.8, and the control parameter Pc0 is set to 0.5. To match the settings of the algorithms used for comparison, the population size N is set to 100, the algorithm runs for 250 generations, and the size of archive set NA is set to 100. The obtained pareto fronts by SA-MOCDE on these ZDTs are depicted in Fig. 2, it is seen that the obtained pareto fronts are very close to the true pareto fronts presented by red lines. In addition, the obtained pareto fronts are distributed uniformly, which reveals that the proposed SA-MOCDE has good diversity distribution in dealing with ZDTs’ problems. Furthermore, the convergence metric and the diversity metric are used to verify the performance of these ZDTs with comparisons to other established methods, and the obtained results are presented in Tables 3 and 4. Here, NSGA-II, SPEA2, PMODE, DEMO and ADEA are taken into comparison, the upper cell points out the average value of convergence or diversity after 10-time calculation, and the lower cell denotes the variance of these calculations. From Tables 3 and 4, both the obtained convergence value and diversity value are smaller and the variances of them are also lower than that of other comparative methods, which reflects that the proposed SAMOCDE have both good convergence ability and diversity distribution as well as relative stable search ability for solving ZDTs. Secondly, some constrained test problems are used to verify the search ability of the proposed SA-MOCDE. Binh (2) and Kita are two

Table 1 Historic knowledge structure.

DT1

DT2

.

DTR

P(DTr)

T1

T2

.

TR

Accept probability

H. Zhang et al. / Energy 50 (2013) 24e37

29

Table 2 Benchmark test problems. Problems

Variable bounds

Objective function

Variable number

Optimal solution

ZDT1

0  xi  1 i ¼ 1,2,.,m

m ¼ 30

ZDT3

0  xi  1 i ¼ 1,2,.,m

ZDT4

0  x1  1 5  x1  5; i ¼ 2,3,.,m 0  xi  1 i ¼ 1,2,.,m

f1(x) ¼ x1 pffiffiffiffiffiffiffiffiffi f2 ðxÞ ¼ gð1  f1 =g Þ P g ¼ 1 þ 9* m i ¼ 2 xi =ðm  1Þ f1(x) ¼ x1 pffiffiffiffiffiffiffiffiffi f2 ðxÞ ¼ gð1  f1 =g  f1 =g*sinð10pf1 ÞÞ Pm g ¼ 1 þ 9* i ¼ 2 xi =ðm  1Þ f1(x) ¼ x1 pffiffiffiffiffiffiffiffiffi f2 ðxÞ ¼ gð1  f1 =g Þ P 2 g ¼ 1 þ 10ðm  1Þ þ m i ¼ 2 ðxi  10*cosð4pxi ÞÞ f1 ðxÞ ¼ 1  expð4xi Þsin6 ð6px1 Þ f2 ðxÞ ¼ gð1  ðf1 =gÞ2 Þ P g ¼ 1 þ 9* m i ¼ 2 xi =ðm  1Þ

x1 ˛ [0,1] xi ¼ 0 i ¼ 2,3,.,m x1 ˛ [0,1] xi ¼ 0 i ¼ 2,3,.,m x1 ˛ [0,1] xi ¼ 0 i ¼ 2,3,2,m x1 ˛ [0,1] xi ¼ 0 i ¼ 2,3,.,m

ZDT6

ZDT1

1

Pareto front by SA-MOCDE True pareto front

m ¼ 10

ZDT3 Pareto front by SA-MOCDE True pareto front

0.5

0.6

0

0.4

-0.5

0.2 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

ZDT4

1

-1

0.6

0.4

0.4

0.2

0.2 0.1

0.2

0.3

0.4

0.5

0.1

0.2

0.3

0.4

0.6

0.7

0.8

0.9

0.5

1

0 0.2

0.6

0.7

0.8

0.9

ZDT6 Pareto front by SA-MOCDE True pareto front

0.8

0.6

0

0

1

Pareto front by SA-MOCDE True pareto front

0.8

0

m ¼ 10

1

0.8

0

m ¼ 10

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 2. Pareto front of ZDT1, ZDT3, ZDT4, ZDT6.

Table 3 The convergence metric g. Problems NSGA-II SPEA2 (real-coded)

PMODE

DEMO/ (parent)

ADEA

SA-MOCDE

ZDT1

0.005800 0 0.021560 0 0.638950 0.500200 0.026230 0.000861

0.001083 0.000113 0.001178 0.000059 0.001037 0.000134 0.000629 0.000044

0.002741 0.000385 0.002741 0.000120 0.100100 0.446200 0.000624 0.000060

0.000063 0.000024 0.000034 0.000079 0.002328 0.000039 0.000032 0.000032

ZDT3 ZDT4 ZDT6

0.033482 0.004750 0.114500 0.007940 0.513053 0.118460 0.296564 0.013135

0.023285 0 0.018409 0 4.9271 2.703 0.232551 0.004945

number is set as 300, and all the other parameters are set as the same as that on ZDTs. In Figs. 3 and 4, the red points denote the obtained pareto front and blue points make up the feasible domain of these test problems, and the pareto front is on the edge of this feasible domain. Thirdly, high-dimensional optimization problems are used to verify the search ability of the proposed SA-MOCDE. Here, the proposed SA-MOCDE is applied on Tamaki and two DTLZs problems [38], and we take 12 decision variables in DTLZ5, and take M ¼ 3

constrained test problems, the obtained feasible domain and pareto front are clearly presented in Figs. 3 and 4. Here, the main parameters are set as follows: The size of evolutionary population is set as 200, the size of archive set is set as 100, maximum generation Table 4 The diversity metric D. Problems NSGA-II SPEA2 (real-coded)

PDEA

DEMO/ (parent)

ADEA

SA-MOCDE

ZDT1

0.298576 0.000742 0.623812 0.000225 0.840852 0.035741 0.473074 0.021721

0.325237 0.030249 0.309436 0.018603 0.359905 0.037672 0.442308 0.021721

0.382890 0.001435 0.525770 0.043030 0.436300 0.110000 0.361100 0.036100

0.074235 0.000068 0.297641 0.000429 0.096258 0.000276 0.078934 0.002312

ZDT3 ZDT4 ZDT6

0.390307 0.001876 0.738540 0.019706 0.702612 0.064648 0.668025 0.009923

0.154723 0.000874 0.469100 0.005265 0.823900 0.002883 1.04422 0.158106

Fig. 3. Pareto front obtained by SA-MOCDE for Binh(2).

30

H. Zhang et al. / Energy 50 (2013) 24e37 9 1

8

0.8

7

0.6 f3

6

0.4

5

0.2

4

0 0 0.2

3

0.4 0.6

2

0.8 1

1 -40

-35

-30

-25

-20

-15

-10

-5

0

5

0

0.2

0.1

0.4

0.3

10

0.5

0.7

0.6

0.8

0.9

1

f2

f1

Fig. 4. Pareto front obtained by SA-MOCDE for Kita. Fig. 6. Pareto front obtained by SA-MOCDE for DTLZ5.

and 20 decision variables. Since DTLZ7 has 2M1(M is the number of objective functions) disconnected Pareto-optimal areas in the search space, population size is set to 500 and the size of archive set is also set to 500, and the pareto fronts of DTLZs and Tamaki are obtained by SA-MOCDE after 500 generations, and the obtained results are presented in Figs. 5, 6 and 7. It is seen that those obtained non-dominated points are uniformly distributed and converge well to the true pareto fronts, which are presented as curve surface or curve line in three-dimensional space. From above result comparisons, the proposed SA-MOCDE can have great search ability in solving ZDTs with good performance both on convergence and diversity distribution, proper constrainthandling ability in dealing with constrained optimization problems and good multiple-objective optimization ability on threedimensional space.

5.1. The solution encoding strategy Due to the high-dimensional non-linear characteristics of SHOSEE problem, a set of water discharge and a set of thermal output are taken as decision variables. For T scheduling periods on the schedule horizon, and there are T schedule periods of water discharge in Nh hydro plants and T schedule periods of thermal output in Ns thermal units, which can be encoded in formulation (25).

2

Qh1;0 Qh2;0 6 Qh1;1 Qh2;1 X¼6 4 « « Qh1;T1 Qh2;T1

/ QhNh ;0 Ps1;0 Ps2;0 / QhNh ;1 Ps1;1 Ps2;1 « « « « / QhNh ;T1 Ps1;T1 Ps2;T1

3 / PsNs ;0 / PsNs ;1 7 7 5 « « / PsNs ;T1 (25)

5. The implementation of SA-MOCDE on hydrothermal operation system scheduling Since the SHOSEE is a complex-constrained problem in the application, proper implementation strategy mainly depends on the constraint-handling technique during the optimization process. The proposed SA-MOCDE is applied on the hydrothermal operation system by pruning the infeasible regions and promoting the promising regions. Here, we can focus on the implementation strategy of SHOSEE problem in detail.

5.2. Population initialization strategy Before the constraint-handling strategy of SHOSEE problem, the initial population can be randomly generated in the feasible interval, which is presented as follows:

xi ¼ li þ Uð0; 1Þ$ðui  li Þ

1

6 5

0.6

4

f3

0.8

3

0.4

2 0

0.2 0 1

(26)

0 0.2 0.4 0.9

0.8

0.7

0.6

0.6 0.5

0.4

0.3

0.2

0.1

0.8 0

1

Fig. 5. Pareto front by SA-MOCDE for Tamaki.

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.9

0.8

0.7

0.6

0.5 f1

0.4

Fig. 7. Pareto front obtained by SA-MOADE for DTLZ7.

0.3

0.2

0.1

0

H. Zhang et al. / Energy 50 (2013) 24e37

31

Initialization: Initialize the population space, belief space and set the generation number g 0

DE Operator

Population space

Evolution

Crossover

Constraint handling

Mutation

Selection

Accept()

Influence()

Accept the nondominated individuals

g mod 20 xi g

1

0 and

xi g

Check the constraint boundary

Simulated Annealing

Situational Knowledge

Historic Knowledge

Normative knowledge

Belief Space g

g 1

g > gmax

NO

Yes Print the result and termination Fig. 8. The flowchart of culture belief based multi-objective hybrid differential evolutionary algorithm.

where ui, li are the upper and lower bounds of decision variable, U(0,1) is a random generated parameter of uniform distribution in (0,1), all the individuals can be properly initialized under control. 5.3. Constraint handling techniques The constraint handling procedure is a key step for solving SHOSEE problem due to its various nonlinear coupled inequality and equality constraints, and the SHOSEE has two sets of equality constraints and four sets of inequality constraints. Compared with the inequality constraints, the equality constraints are more difficult to be properly handled, the system load balance constraint is time independent while terminal reservoir storage is time-related. Since the modification of water discharge and reservoir storage will

affect system load allocation, water dynamic balance constraint must be properly handled first. Since reservoir storage is closely related to the water discharge and water inflow of hydro plant, the terminal reservoir storage can be properly controlled by modifying the water discharge process. Here, the constraint of terminal reservoir storage is handled with coarse and fine tuning iterations max, Jmax, coarse adjustment allocates the total violation equally to the hydro plants, and fine tuning adjustment modifies the water discharge and reservoir storage violation. After the modification of water discharge, the output of hydro plant also changes. To maintain the system load balance, the constraint adjustment is mainly modified by thermal power allocation. The constraint handling technique of system load balance is similar to the method used to handle terminal reservoir

32

H. Zhang et al. / Energy 50 (2013) 24e37

storage constraint with coarse and fine tuning iterations max, Jmax, and the adjusting technique of terminal reservoir storage and system load balance has been introduced in literature [39]. After the constraint-handling modification on infeasible individual, the modified individual can still be infeasible due to its finite iterations. Generally, the permitted accuracy can be used to distinguish the feasibility of individual Xi with its total violation characteristic, which is presented in the formula (27).

Total violationðXi Þ ¼

XNh k¼1

jDVk ðXi Þj þ

XT t¼1

jDPt ðXi Þj

(27)

To obtain the dominance relationship between two selected individuals, the feasible selection strategy proposed by Deb can be followed in literature [40]. 5.4. The flowchart of SA-MOCDE on the application of hydrothermal system The implementation of the proposed SA-MOCDE on hydrothermal system can be described as the following flowchart in Fig. 8. 6. Case study 6.1. Description of hydrothermal system operation Here, the hydrothermal system consisting of a multi-chain cascade of four hydro plants and three thermal units is studied to demonstrate the feasibility and effectiveness of the proposed SAMOCDE for solving daily hydrothermal scheduling with economic emission problem. The total scheduling period is 24 h with an hour interval for each scheduling period, and the hydraulic relationships among those four hydro plants are shown in Fig. 9, where four hydro-plants are connected in series and parallel. While the generated hydro power is closely related to the thermal power by system load requirements as it is shown in formulation (6), and the hydro power depends on the reservoir storage and water discharge. Thus, the deviation of reservoir storage or water discharge can affect the power dispatch of thermal units, which influences the total fuel cost and emission issue. The case study consists of two test systems, for simplicity system transmission loss isn’t taken into consideration in test system 1, while in test system 2, it is completely considered. Moreover, the obtained scheduling results are taken into

comparison with that obtained by the comparative method, and all the details of hydrothermal system are given in literature [15,41]. 6.2. Test system 1 6.2.1. The parameter setting The main parameters are set as follows: The generation number is set to 1000, and the size of evolutionary population and archive set are set to 50 and 30. The initial crossover rate Pc0 is set to 0.4, and both the tuning iteration number max, Jmax are set to 10. The permitted accuracies in handling terminal reservoir storage εv and system load balance εp are set to 0.01, and the permitted total accuracy εtotal is set to 0.1. 6.2.2. Results and comparison In this test system, the transmission loss and spillage water are neglected for simplicity, and the implementation of the proposed SA-MOCDE is applied on this test system while comparing the obtained results with MODE [22] and HMOCA [39] under same parameter settings. 30 non-dominated schemes obtained by SA-MOCDE and MODE are presented in Table 5, the total computation time is 743 s in comparison to 752 s by MODE, which means that it takes about 24.7 s while 25.0 s with MODE to generate one scheme. It is seen that both the thermal cost and emission cost obtained by SAMOCDE are smaller than that of MODE. In addition, the nondominated schemes obtained by SA-MOCDE, MODE and HMOCA [39] are presented as pareto fronts depicted in Fig. 10. It is found that the obtained pareto front by SA-MOCDE has more desirable distribution than that of MODE, while HMOCA shows great ability of optimizing emission cost in spite of low efficiency on thermal cost. For further analysis on these schemes, scheme (15) is taken as the compromise scheme, and more details about compromise scheme are shown in Table 6, where water discharge, hydro output and thermal output of hydrothermal system are presented. According to the details of hydrothermal system in literature [41], water discharge, hydro output and thermal output at each scheduling period are properly under control, and the summation of hydro output and thermal output keeps consistence with the system load at each time period. Moreover, reservoir storages of four hydro plants are depicted in Fig. 11, and it is found that the reservoir storage at each time period is also under control. In Table 7, literature [17] and literature [42] provide pure scenarios on this hydrothermal system, it is seen that economic

Fig. 9. Hydraulic relationships among the four hydro-plants.

H. Zhang et al. / Energy 50 (2013) 24e37

33

Table 5 The comparison of non-dominated schemes between SA-MOCDE and MODE for test system 1. Scheme SA-MOCDE

MODE

Scheme SA-MOCDE

Thermal cost ($) Emission cost (lb) Thermal cost ($) Emission cost (lb) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

42,038 42,052 42,122 42,165 42,180 42,200 42,287 42,312 42,367 42,425 42,508 42,624 42,690 42,924 42,988

17,450 17,321 17,216 17,134 17,073 16,938 16,825 16,737 16,651 16,623 16,601 16,527 16,494 16,455 16,440

42,252 42,273 42,316 42,434 42,483 42,621 42,683 42,715 42,849 42,997 43,010 43,107 43,217 43,307 43,489

17,576 17,418 17,381 17,323 17,207 17,138 17,089 16,957 16,924 16,863 16,792 16,749 16,669 16,592 16,519

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

SA-MOCDE MODE

Emission cost (lb)

HMOCA

17000

16500

16000

15500 41000

43000

45000 Thermal cost ($)

47000

43,032 43,252 43,319 43,481 43,895 44,099 44,204 44,355 44,443 44,561 44,754 44,802 45,041 45,151 45,503

16,371 16,335 16,311 16,263 16,203 16,170 16,159 16,139 16,126 16,112 16,095 16,078 16,064 16,040 16,004

43,586 43,627 43,729 43,816 43,904 44,092 44,113 44,217 44,393 44,593 44,729 44,914 45,017 45,032 45,019

16,474 16,447 16,392 16,363 16,338 16,307 16,282 16,230 16,184 16,149 16,104 16,079 16,077 16,068 16,072

load scheduling (ELS) is pure thermal cost scheduling with fuel cost 43,500($) in 72.96 s, economic emission scheduling (EES) is pure emission cost scheduling with emission rate 18,257(lb) in 72.74 s, and combined economic emission scheduling (CEES) is economic cost scheduling on considering emission cost with thermal cost 44,914($) and emission rate 19,615(lb) in 74.97 s. While in comparison to those obtained results in literature [17], the thermal cost obtained by SA-MOCDE is 1005($) more than that of ELS, the obtained emission cost is 363(lb) lower than that of EES, and the obtained thermal cost is 57($) less than that of CEES while emission cost is 563(lb) lower than that of CEES, which means that the proposed SA-MOCDE is more efficient than that in literature [17] and literature [42].

18000

17500

MODE

Thermal cost ($) Emission cost (lb) Thermal cost ($) Emission cost (lb)

49000

Fig. 10. The pareto front obtained by SA-MOCDE, MODE and HMOCA in test system 1.

6.2.3. Analysis on running performances of convergence and diversity Here, we can focus on the convergence and diversity performance of SA-MOCDE in solving this hydrothermal problem [43]. The convergence metric CðgÞ mainly provides the relative clustering degree between archive set and the true pareto front in the

Table 6 The details of compromise scheduling result by SA-MOCDE for test system 1. Hours

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Water discharge (104*m3/s)

Hydro power (MW)

Thermal power (MW)

Load (MW)

Qh1

Qh2

Qh3

Qh4

Ph1

Ph2

Ph3

Ph4

Ps1

Ps2

Ps3

8.832 8.136 6.998 8.117 7.833 6.055 9.513 9.829 9.675 8.191 9.623 8.901 9.75 9.646 9.248 6.805 9.567 9.186 9.247 8.951 5.5 5.398 5 5

6 6.03 6 6.002 6.059 6 6.33 8.462 8.719 8.163 8.513 8.048 9.207 9.694 8.612 6.627 10.181 10.817 11.111 10.874 8.666 12.154 11.252 8.479

22.652 26.428 23.745 20.536 18.456 19.041 19.109 17.174 17.602 19.101 17.929 17.573 18.454 17.684 17.553 16.406 15.387 14.968 14.652 15.711 11.334 11.857 12.563 11.499

6 6 6 6 6 6 9.619 16.573 17.281 18.185 19.109 17.174 17.602 19.101 17.929 17.573 18.454 17.954 19.939 20 19.337 18.998 19.826 16.308

80.057 76.315 68.937 76.73 74.606 61.336 84.009 85.029 83.994 75.996 84.829 81.736 86.571 86.426 84.974 69.357 87.928 85.696 85.767 83.54 58.144 57.501 54.302 54.705

49 50.377 51.279 52.935 54.908 55.459 58.326 71.917 72.605 68.992 71.495 69.006 75.59 77.45 71.265 59.011 80.435 81.491 79.758 76.098 63.576 79.167 73.664 58.426

31.824 0.171 12.412 26.686 34.09 30.681 29.985 36.252 33.525 27.141 31.01 32.827 30.125 34.463 36.06 40.809 44.817 47.509 48.799 47.824 52.09 54.998 57.479 57.765

131.88 129.027 125.744 121.625 115.822 131.396 195.038 276.407 285.677 293.4 300.482 286.793 290.019 300.428 292.412 289.809 296.109 292.587 305.361 303.27 295.558 289.318 289.906 260.93

109.792 175 174.723 105.706 121.322 174.557 174.9 175 175 174.966 174.966 175 175 175 174.953 174.973 175 174.926 175 174.921 174.921 115.11 108.333 102.713

209.114 209.603 124.751 134.625 127.484 208.079 268.764 213.499 217.56 209.906 209.881 275.121 228.463 215.372 210.133 209.78 218.744 214.188 230.423 210.188 124.925 124.94 125.412 122.403

138.334 139.507 142.154 131.693 141.767 138.492 138.978 151.896 221.64 229.599 227.337 229.517 224.232 140.861 140.202 216.261 146.967 223.603 144.892 154.159 140.787 138.966 140.903 143.057

750 780 700 650 670 800 950 1010 1090 1080 1100 1150 1110 1030 1010 1060 1050 1120 1070 1050 910 860 850 800

34

H. Zhang et al. / Energy 50 (2013) 24e37 Convergence of SA-MOCDE

Convergence of MODE

Convergence of HMOCA

Diversity of SA-MOCDE

Diversity of MODE

Diversity of HMOCA

1 0.9

Evolutionary performance

0.8

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

200

400

600

800

1000

Generations

Fig. 12. The running performance of convergence and diversity for test system 1. Fig. 11. The trajectory of four reservoir storages in test system 1.

Table 7 Scheduling results in Refs. [17,42]. Objective

Ref. [42]

Ref. [17]

ELS

EES

CEES

ELS

EES

CEES

Fuel cost ($) Emission (lb) Comp. time (s)

43,500 21,092 72.96

51,449 18,257 72.74

44,914 19,615 74.97

41,983 24,482 112.00

44,432 16,803 112.56

43,045 17,003 120.00

g-th generation, and it is also within the range of [0,1]. While the diversity metric DðgÞ describes the diversity distribution of obtained pareto-optimal individuals on the (M1)-dimensional plane in the g-th generation, and it is also limited in [0,1]. Since the true pareto front is impossible to be known in solving the SHOSEE problem, it can be replaced by P * ¼ non dominatedðW1000 g¼1 Pg Þ, where Pg represents the obtained pareto front in the g-th generation. Here, 1000 generations are divided into 20 intervals with 50 generations at each interval, and the true pareto front can be presented as P * ¼ non dominatedðW20 g¼1 P50*g Þ for simplicity. Then, the running performance of convergence and diversity metric is shown in Fig. 12. In Fig. 12, it is seen that the population convergence decreases and diversity distribution increases as differential evolution proceeds. Since the simulated annealing technique adjusts the population diversity every 20 generations, convergence process of MODE may have rapid convergence speed at the beginning and reach a stable state early, while the trajectory of proposed SAMOCDE converges more slowly than that of MODE, which can avoid premature convergence problem efficiently. Moreover, HMOCA has better diversity distribution than SA-MOCDE in the whole evolution process, but it has lower convergence rate due to its low convergence ability on thermal cost. And simultaneously the diversity trajectory of SA-MOCDE has better distribution than that of MODE throughout the whole evolutionary process.

6.3. Test system 2 6.3.1. The parameter setting Since the complexity of test system increases, the main parameters are set as follows: The generation number is set to

2000, and the size of evolutionary population and archive set are set to 100 and 30. The initial crossover rate Pc0 is set to 0.4, and both ’ ’ ; Jmax are set to 20. The permitted the tuning iteration number Imax accuracies in handling terminal reservoir storage εv and system load balance εp are set to 0.01, and the permitted total accuracy is set to 0.1. 6.3.2. Results and comparison In this test system, the transmission loss of system load is taken into consideration, and it converts the system load balance into a nonlinear balance according to formulation (7) and formulation (8), which increases the difficulty for properly handling the system load balance. Here, we still utilize the constraint-handling technique introduced in Section 5.3, the thermal output can be properly controlled in an acceptable accuracy. Here, the MODE is taken into comparison with the same parameter setting while utilizing the constraint-handling technique. In Table 8, the total computation time of generating 30 nondominated schemes is 1092 s in comparison to 1098 s by MODE, which means that it takes about 36.4 s while 36.6 s with MODE to generate one scheme. It is seen that the thermal cost is sorted in ascending order while the emission cost is in descending order, and it is also found that the obtained thermal cost and emission cost are smaller than that of MODE. From Fig. 13, the pareto front obtained by SA-MOCDE is under that obtained by MODE, and the nondominated points obtained by SA-MOCDE are well distributed while those points obtained by MODE are disorderly distributed. In comparison to those obtained pareto front in literature [39], SAMOCDE and MODE have higher efficiency on thermal cost optimization, while HMOCA has lower emission cost in relative narrow scope, which is merely range in between 17,000 and 18,000. In addition, scheme (15) is also taken as compromise scheme for analysis on the scheduling results, and more details about the compromise scheme are shown in Table 9. It can be found that water discharge, hydro output and thermal output are properly under control, and the summation of hydro output and thermal output keeps the balance with system load and transmission loss. Moreover, the trajectory of reservoir storage is depicted in Fig. 14, where all the reservoir storage and terminal storage satisfy the constraint limits. And it is seen that the water storages of reservoir 1 and reservoir 2 are lower than that of reservoir 3 and reservoir 4. Since few literature have taken research on hydrothermal system with considering the transmission loss, we can focus on

H. Zhang et al. / Energy 50 (2013) 24e37

35

Table 8 The comparison of non-dominated schemes between SA-MOCDE and MODE for test system 2. Scheme SA-MOCDE

MODE

Scheme SA-MOCDE

Thermal cost ($) Emission cost (lb) Thermal cost ($) Emission cost (lb) 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

42,167 42,176 42,189 42,206 42,286 42,354 42,422 42,505 42,594 42,678 42,812 42,890 43,013 43,097 43,165

19,981 19,653 19,376 19,206 19,002 18,731 18,501 18,323 18,105 17,796 17,702 17,625 17,521 17,498 17,464

42,248 42,289 42,326 42,354 42,471 42,519 42,667 42,756 42,812 42,855 42,878 42,915 43,022 43,127 43,245

20,014 19,376 19,206 18,971 18,900 18,694 18,523 18,235 18,132 17,964 17,796 17,753 17,724 17,581 17,523

16 17 18 19 20 21 22 23 24 25 26 27 28 29 30

Emission cost (lb)

SA-MOCDE MODE

19500

HMOCA

43,214 43,350 43,467 43,556 43,720 43,811 43,932 44,120 44,155 44,165 44,321 44,542 44,641 44,956 45,264

17,354 17,235 17,132 17,056 17,039 17,009 16,931 16,929 16,897 16,863 16,832 16813 16,841 16,747 16,713

43,324 43,527 43,650 43,756 43,820 43,865 43,925 43,974 44,024 44,178 44,355 44,365 44,542 45,156 45,368

17,493 17,426 17,235 17,156 17,139 17,065 17,042 17,021 17,001 16,999 16,997 16,963 16,913 16,857 16,819

deviation in comparison to test system 1. Due to non-linear characteristics of system load constraint, it increases the difficulty for proper constraint-handling, which can also lead an increase on the total fuel cost as well as emission cost in comparison to test system 1. From Table 9, it is found that the transmission loss is under 10 MW, and the successful transmission rate of each time period is higher than 95%. In addition, the obtained trajectory of water discharge and reservoir storage totally keeps consistent with that of test system 1.

20500 20000

MODE

Thermal cost ($) Emission cost (lb) Thermal cost ($) Emission cost (lb)

19000 18500 18000 17500 17000 16500 41500

43500

45500 Thermal cost ($)

47500

Fig. 13. The pareto front obtained by SA-MOCDE, MODE and HMOCA in test system 2.

6.3.3. Analysis on running performances of convergence and diversity Similarly, 2000 generations are divided into 20 intervals with 100 generations each interval, and the true pareto front can be replaced by P * ¼ nondominatedðW20 g¼1 P100*g Þ for simplicity, and the running performance of convergence and diversity of SA_MOCDE and MODE is shown in Fig. 15. It is seen from Fig. 15 that the population convergence decreases as the evolution proceeds, while the diversity indictor increases in

Table 9 The details of compromise schedule result by SA-MOCDE for test system 2. Hours

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Water discharge (104*m3/s)

Hydro power (MW)

Thermal power (MW)

Q1

Q2

Q3

Q4

Ph1

Ph2

Ph3

Ph4

Ps1

Ps2

Ps3

8.828 7.347 6.444 8.928 8.56 7.9 10.095 8.337 9.742 9.468 9.625 9.367 10.323 9.203 10.261 7.225 6.346 8.432 6.868 7.558 6.669 5.694 5.666 6.115

6.136 6 6 6 6.396 6.471 7.497 6 9.389 7.509 8.876 7.385 9.958 11.717 9.735 7.887 7.186 11.365 8.886 10.887 10.332 10.728 11.3 8.361

23.448 23.595 21.889 20.27 17.366 18.321 18.387 19.419 17.977 18.079 15.681 17.918 17.396 17.001 17.344 17.76 17.923 16.332 14.772 14.652 11.735 12.962 13.566 13.334

6 6 6 6 6 6 6.618 13.384 17.366 18.623 18.262 19.242 18.755 18.494 16.204 16.202 17.396 19.185 18.03 19.959 19.826 19.646 19.943 19.196

80.03 71.162 65.011 81.944 79.258 74.453 86.063 76.186 83.902 82.665 83.986 83.498 88.289 83.05 89.174 71.636 65.558 80.903 69.992 75.041 68.226 60.365 60.428 64.623

49.963 50.086 51.22 52.861 57.255 58.596 65.727 54.754 76.457 64.814 73.747 64.867 79.806 86.372 76.053 65.42 61.078 82.866 67.931 76.272 71.935 72.729 73.798 57.692

27.087 20.631 24.892 29.656 39.34 35.261 35.26 30.671 35.141 34.601 41.193 35.104 38.085 40.769 41.349 41.279 41.984 48.524 52.146 52.432 53.979 56.251 57.349 58.463

131.88 129.027 125.744 121.625 115.822 132.099 155.2 247.326 288.258 297.259 294.482 301.16 298.15 295.621 277.963 277.463 288.483 300.968 291.009 302.808 299.738 296.675 294.871 284.891

112.014 164.961 172.408 174.547 116.251 174.907 175 175 174.753 174.927 175 174.954 175 174.93 174.988 174.479 166.857 175 175 174.955 154.323 102.746 102.658 102.722

209.404 208.476 125.084 125.033 125.063 124.903 209.725 209.824 209.805 209.799 209.849 209.68 209.77 127.148 209.667 209.771 206.671 209.584 208.632 209.611 124.828 125.141 124.926 124.74

143.348 140.245 139.365 67.225 139.866 204.609 229.499 222.572 228.122 222.259 228.195 288.598 227.333 227.396 145.696 226.355 225.547 228.613 211.372 164.081 140.401 148.891 138.629 109.17

PL (MW)

Load (MW)

3.726 4.588 3.724 2.891 2.855 4.828 6.474 6.333 6.438 6.324 6.452 7.861 6.433 5.286 4.89 6.403 6.178 6.458 6.082 5.2 3.43 2.798 2.659 2.301

750 780 700 650 670 800 950 1010 1090 1080 1100 1150 1110 1030 1010 1060 1050 1120 1070 1050 910 860 850 800

36

H. Zhang et al. / Energy 50 (2013) 24e37

constraints. In this paper, a simulated annealing technique based multi-objective cultural differential evolutionary algorithm is proposed to solve this SHOSEE problem. The simulated annealing technique is incorporated into the cultural computational model, and it prevents the premature problem by proper controlling evolutionary population diversity. Furthermore, the crossover rate is also modified to an adaptive parameter, which can properly control the population convergence. Some benchmark problems are used to testify the convergence ability and diversity distribution of pareto front by SA-MOCDE, and the efficiency of SA-MOCDE is also verified by three-dimensional, complex constrained test functions further. Ultimately, the proposed SA-MOCDE is applied on the SHOSEE problem with comparison to other alternatives, the obtained results have more desirable thermal cost and emission cost, and the analysis on the running performance of convergence and diversity also reveals that SA-MOCDE can be a robust and effective alternative for solving SHOSEE problem. Acknowledgment Fig. 14. The trajectory of four reservoir storages in test system 2.

Convergence of SA-MOCDE

Convergence of MODE

Convergence of HMOCA

Diversity of SA-MOCDE

Diversity of MODE

Diversity of HMOCA

1

This work is granted by the Public Welfare Industry of the Ministry of Water Resources (No. 201001080), the special research fund for high school doctoral program (No. 20100142110012), and National Natural Science Foundation of China (No. 51107047). References

0.9

Evolutionary performance

0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0

200

400

600

800

1000 1200 Generations

1400

1600

1800

2000

Fig. 15. The running performance of convergence and diversity for test system 2.

the evolution process. In addition, since simulated annealing technique has been utilized to increase the population diversity, the trajectory of population convergence obtained by SA-MOCDE converges more slowly and reaches lower convergence degree than that of MODE and HMOCA, and generally the trajectory of evolutionary diversity obtained by SA-MOCDE has more desirable distribution than that obtained by MODE and HMOCA especially when population comes to converge. From the above comparison results, it is found that the proposed SA-MOCDE can obtain better scheduling results than other comparative alternatives, and the analysis on the running performance proves that the proposed SA-MOCDE also has desirable convergence and diversity performance in solving hydrothermal scheduling problem. 7. Conclusion Daily hydrothermal optimization scheduling with economic emission is a non-linear, multi-phase and complex multi-objective operational problem with considering coupled dynamic

[1] Ruzic S, Rajakovic N, Vuckovic A, Utility E. A flexible approach to short-term hydro-thermal coordination. I. Problem formulation and general solution procedure. IEEE Trans Power Syst 1996;11(3):1564e71. [2] Wollenberg B, wood A. Power generation, operation and control. John Wiley &Sons; 1996. p. 264e327. [3] Saha T, Khaparde S. An application of a direct method to the optimal scheduling of hydrothermal system. IEEE Trans Power Appl Syst 1978:977e83. [4] Guan X, Peter B. Nonlinear approximation method in Lagrangian relaxation based algorithms for hydrothermal scheduling. IEEE Trans Power Syst 1998; 13(1):226e35. [5] Brannud H, Bubenko JA, Sjelvgren D. Optimal short term operation planning of a large hydrothermal power system based on a nonlinear network flow concept. IEEE Trans Power Syst 1986;1(4):75e82. [6] Orero S, Irving M. A genetic algorithm modelling framework and solution technique for short term optimal hydrothermal scheduling. IEEE Trans Power Syst 1998;13(2):501e18. [7] Kumar S, Naresh R. Efficient real coded genetic algorithm to solve the nonconvex hydrothermal scheduling problem. Int J Elect Power Energy Syst 2007;29(10):738e47. [8] Hota P, Chakrabarti R, Chattopadhyay P. Short-term hydrothermal scheduling through evolutionary programming technique. Electr Power Syst Res 1999; 52(2):189e96. [9] Yu B, Yuan X, Wang J. Short-term hydro-thermal scheduling using particle swarm optimization method. Energy Convers Manage 2007;48(7):1902e8. [10] Basu M. Artificial immune system for fixed head hydrothermal power system. Energy 2011;36:606e12. [11] El-Keib AA, Ma H, Hart JL. Economic dispatch in view of the clean air act of 1990. IEEE Trans Power Syst 1994;9(2):972e8. [12] Talaq JH, El-Hawary F, El-Hawary ME. A summary of environmental/economic dispatch algorithms. IEEE Trans Power Syst 1994;9(3):1508e16. [13] Helsin JS, Hobbs BF. A multiobjective production costing model for analyzing emission dispatching and fuel switching. IEEE Trans Power Syst 1989;4(3): 836e42. [14] Vaisakh K, Srinivas LR. A genetic evolving ant direction DE for OPF with nonsmooth cost functions and statistical analysis. Energy 2010;35:3155e71. [15] Basu M. An interactive fuzzy satisfying method based on evolutionary programming technique for multiobjective short-term hydrothermal scheduling. Electr Power Syst Res 2004;69:277e85. [16] Basu M. A simulated annealing-based goal-attainment method for economic emission load dispatch of fixed head hydrothermal power systems. Electr Power Energy Syst 2005;27(2):147e53. [17] Mandal KK, Chakraborty N. Daily combined economic emission scheduling of hydrothermal systems with cascaded reservoirs using self organizing hierarchical particle swarm optimization technique. Expert Syst Appl 2012;39: 3438e45. [18] Deb K, Pratap A. A fast and elitist multi-objective genetic algorithm: NSGA-II. IEEE Trans Evol Comput 2002;6(2):182e97. [19] Zitzler E, Laumanns M, Thiele L. SPEA2: improving the strength Pareto evolutionary algorithm. Technical report 103, computer engineering and

H. Zhang et al. / Energy 50 (2013) 24e37

[20] [21]

[22]

[23] [24] [25]

[26]

[27]

[28]

[29]

[30] [31]

networks laboratory (TIK), Swiss Federal Institute of Technology (ETH) Zurich, Gloriastrasse 35, CH-8092 Zurich, Switzerland; 2001. http://www. kddresearch.org/Courses/Spring-2007/CIS830/Handouts/P8.pdf. Wang Y, Yang Y. Particle swarm with equilibrium strategy of selection for multi-objective optimization. Eur J Oper Res 2010;200(1):187e97. Tripathi PK, Bandyopadhyay S, Pal SK. Multi-objective particle swarm optimization with time variant inertia and acceleration coefficients. Inform Sci 2007;177(22):5033e49. Xue F, Sanderson AC, Graves RJ. Pareto-based multi-objective differential evolution. In: Proceedings of the 2003 congress on evolutionary computation, vol. 2. Canberra, Australia: IEEE Press; 2003. p. 862e9. Rolic T, Filipic B. DEMO: differential evolution for multi-objective optimization, lecture notes in computer science. Berlin: Springer; 2005. p. 520e33. Qian WY, Li AJ. Adaptive differential evolution algorithm for multi-objective optimization problems. Appl Math Comput 2008;201(1):431e40. Storn R, Price K. Differential evolution: a simple and efficient adaptive scheme for global optimization over continuous spaces. Technical Report TR-95e012, Berkeley, USA: International Computer Science Institute; 1995. http://www. icsi.berkeley.edu/ftp/pub/techreports/1995/tr-95-012.pdf. Abbass HA, Sarker R, Newton C. PDE: a pareto-frontier differential evolution approach for multi-objective optimization problems. In: Proceedings of the congress on evolutionary computation 2001 (CEC2001), vol. 2. Piscataway (NJ): IEEE Service Center; 2001. p. 971e8. Madavan NK. Multi-objective optimization using a Pareto differential evolution approach. In: Proceeding of the congress on evolutionary computation (CEC2002), vol. 2. Piscataway (NJ): IEEE Service Center; 2002. p. 1145e50. Sivasubramani S, Shanti Swarup K. Hybrid DE-SQP algorithm for non-convex short term hydrothermal scheduling problem. Energy Convers Manage 2011; 52:757e61. Lakshminarasimman L, Subramanian S. Short-term scheduling of hydrothermal power system with cascaded reservoirs by using modified differential evolution; 2006. Coello CAC, Pulido GT, Lechuga MS. Handling multiple objectives with particle swarm optimization. IEEE Trans Evol Comput 2004;8(3):256e79. Reynolds R. An introduction to cultural algorithms. In: Proceedings of the 3rd annual conference on evolutionary programming; 1994. p. 131e9.

37

[32] Zitzler E, Laumanns M, Thiele L. SPEA2: improving the strength Pareto evolutionary algorithm. Technical report 103, computer engineering and networks laboratory (TIK), Swiss Federal Institute of Technology (ETH) Zurich, Gloriastrasse 35, CH-8092 Zurich, Switzerland; 2001. http://www. kddresearch.org/Courses/Spring-2007/CIS830/Handouts/P8.pdf. [33] Youlin L, et al. Short-term scheduling optimization for hydro-thermal power systems based on adaptive hybrid differential evolution algorithm. Power Syst Technol 2009;33(13):32e6. [34] Saleem SM. Knowledge-based solution to dynamic optimization problems using cultural algorithms, PhD thesis, Wayne State University, Detroit, Michigan; 2001. [35] Gelbard Roy, Carmeli Abraham, Bittmann Ran M, Ronen Simcha. Cluster analysis using multi-algorithm voting in cross-cultural studies. Expert Syst Appl 2009;36(7):10438e46. [36] Metropolis N, Rosenbluth A, Rosenbluth M, Teller A, Teller E. Equation of state calculation by fast computing machines. J Chem Phys 1953;21:1087e92. [37] Kirkpatrick S, Gelatt Jr CD, Vecchi MP. Optimization by simulated annealing. Science May 1983;220(4598):671e80. [38] Deb K, Lothar T. Scalable test problems for evolutionary multi-objective optimization. TIK-Report NO.112. Computer engineering and networks laboratory (TIK); 2001. http://www.tik.ee.ethz.ch/sop/publicationListFiles/ dtlz2001a.pdf. [39] Youlin Lu, Jianzhong Zhou, Hui Qin, et al. A hybrid multi-objective cultural algorithm for short-term environmental/economic hydrothermal scheduling. Energy Convers Manage 2011;52(5):2121e34. [40] Deb K. Multi-objective optimization using evolutionary algorithms. Chichester, UK: Wiley; 2001. p. 291e5. [41] Hota PK, Barisal AK, Chakrabarti R. An improved PSO technique for short-term optimal hydrothermal scheduling. Electr Power Syst Res 2009;79(7):1047e 53. [42] Mandal KK, Chakraborty N. Short-term combined economic emission scheduling of hydrothermal power systems with cascaded reservoirs using differential evolution. Energy Convers Manage 2009;50(1):97e104. [43] Deb K, Jain S. Running performance metrics for evolutionary multi-objective optimization. In: Proceedings of the fourth Asia-Pacific conference on simulated evolution and learning; 2002. p. 13e20.