Dynamic environmental economic dispatch of hybrid renewable energy systems based on tradable green certificates

Journal Pre-proof Dynamic environmental economic dispatch of hybrid renewable energy systems based on tradable green certificates Xiaozhu Li, Weiqing Wang, Haiyun Wang, Jiahui Wu, Xiaochao Fan, Qidan Xu PII:

S0360-5442(19)32394-1

DOI:

https://doi.org/10.1016/j.energy.2019.116699

Reference:

EGY 116699

To appear in:

Energy

Received Date: 16 April 2019 Revised Date:

10 October 2019

Accepted Date: 4 December 2019

Please cite this article as: Li X, Wang W, Wang H, Wu J, Fan X, Xu Q, Dynamic environmental economic dispatch of hybrid renewable energy systems based on tradable green certificates, Energy (2020), doi: https://doi.org/10.1016/j.energy.2019.116699. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Dynamic Environmental Economic Dispatch of Hybrid Renewable Energy Systems based on Tradable Green Certificates Xiaozhu Li1, Weiqing Wang1, Haiyun Wang1, Jiahui Wu1, Xiaochao Fan1, Qidan Xu2 (1. Engineering Research Center of Ministry of Education for Renewable Energy Generation and Grid Connection Technology, Xinjiang University, Urumqi 830047, Xinjiang, China. 2. State Grid Xinjiang Electric Power Co., Ltd. Economic and Technological Research Institute, Urumqi 830000, Xinjiang, China） Abstract: China is accelerating the formulation of a renewable portfolio standard to promote the development and utilization of renewable energy. The implementation of the standard and corresponding tradable green certificates will change the scheduling procedures of thermal power and other renewable energy generators in the electric power system. This paper formulates a multi-objective dynamic economic emission dispatch model for wind-solar-hydro power under tradable green certificates. Weibull, Beta and Gumbel distributions are used to simulate the random behavior of wind, solar and hydro power, respectively, thus taking into account the overestimation and underestimation costs caused by the randomness of renewable energy. Moreover, a new multi-objective moth-flame optimization (MOMFO) based solution is proposed according to the characteristics of the dispatch model. The MOMFO is based on the external elite reservation and Pareto dominance. To improve the search capability, the strategy of average entropy initialization, Lévy flight and variable scale chaotic strategy are introduced. A heuristic dynamic relaxation method is proposed, which is used to handle the complex constraints of the model. Finally, the feasibility of the proposed scheduling model and the effectiveness of the algorithm are illustrated using examples of different scheduling situations based on IEEE 39-bus standard system. Keywords: Dynamic economic environment dispatch, Renewable portfolio standard, Tradable green certificate, Hybrid renewable energy system, Uncertainty, multiobjective moth-flame optimization NW

Nomenclature Ns

Number of wind power plants

Numbers of wind-powered generator.

Scheduled / available (actual) wind power from the jth wind PWi,t, PWi,tur,t

NR

Number of small-hydro power plants

power generator at the time-instant t.

Direct / penalty / reserve cost coefficient for jth wind power FWi, kpWi, krWi

Direct / penalty / reserve cost coefficient for lth small-hydro FRi, kpRi, krRi

generator.

power generator.

Direct / penalty / reserve cost coefficient for kth solar power FSi, kpSi, krSi

Scheduled / available(actual) small-hydro power from the lth PRi,t, PRi,ture,t

generator.

small-hyro power generator at the time-instant t.

Scheduled / available(actual) solar power from the kth

QW-TGC,t, QH-TGC,t,

Transaction quantity of TGC

for

wind power,

Solar-power generator at the time-instant t.

QS-TGC,t

hydropower and solar power at the time-instant t.

Load / loss at the time-instant t.

WSin, WSoff, WSR

Cut-in, cut-off and rated wind speed.

PSi,t, PSi,tur,t Pload,t, Ploss,t WSt

Wind speed at time interval t. 2

kt , ct

Shape and scale parameters of Weibull distribution at time interval t

SIt

Solar irradiance, W/m .

λ γ

Shape parameters of the gummble distribution function

at, bt

Shape parameters of the beta distribution function at time interval t.

RIt

Stream velocity at time interval t

PWR

Rated power output of wind generation unit.

PSR

Rated power output of PV

1. Introduction

time, renewable energy power will be sold as a commodity in

To promote the development and utilization of renewable

the trading market to ensure the subject of obligations can

energy and ensure the upgrading of energy structure according

complete their quota of renewable energy production and use.

to the “13th Five-Year Plan” for China's electric power

The three editions of RPS draft are detailed and clear in the

development, the State Energy Administration issued the first

formulation, implementation, supervision, assessment, the

three editions of “renewable portfolio standard (RPS) and

development of RPS policy, which has been continuing for

assessment measures drafts” on March 23, September 13 and

many years, is about to be finalized. In future, the share of

November 13, 2018. A tradable green certificate (TGC) is a

China's renewable energy consumption will be seeing a steep

complementary transaction system of the RPS. At the appointed

increase.

Dynamic environmental economic dispatch (DEED) is a

consists of addition of individual objective functions, not a real

common multi-objective optimization problem with high

multi-objective

latitude, non-linearity and strong coupling in power systems.

economic dispatch model including wind-solar-small hydro

[1][2]

Traditional methods such as price penalty factors

algorithm.

A

multi-objective

stochastic

, weighted

power is proposed in [22], but the dynamic changes of load

and semi-positive definite programming strategy , weight sum

demand and renewable are not taken into account in the model.

[3]

[4]

mode

etc. But the non-differentiable or non-convex and can

Ref. [23] proposed a dynamic day-ahead stochastic scheduling

not be effectively solved in these methods, meanwhile the

for thermal-hydro-wind-PV systems, but the model is a single

solution is sensitive to the selection of initial value and easy to

objective only considered the economic benefits. In [23],

fall into local optimum. As a result more and more scholars

dynamic characteristics are considered for establishment of a

regard DEED problem as a real multi-objective optimization

scheduling model for wind-solar system and an improved

problem

[5-16]

to

solve

.

For

example,

fast

multi-objective genetic algorithm 2 (NSGA-II) [10]

annealing

and

[5][6]

elitist

fireworks algorithm is used to solve the model. Ref. [24]

, simulated

presented a multi-energy scheduling scheme based on

, multi-objective particle swarm optimization

[11]

decomposition

based

which

increases

renewable

energy

Multi-objective

penetration while reducing operating costs. At present, most of

, these algorithms can

the research studying the impact of RPS on power systems is

effectively solve the traditional DEED problem and have been

concentrated on the macro level. Refs. [25][26][27] studied the

successfully applied. With the large-scale renewable energy

impact of RPS on the development of China's wind power

(MOPSO)

,

biogas-solar-wind,

)[13]

evolutionary algorithms (MOEA/D

interconnection, Energy Internet

[14]

and various multi-energy

industry based on a system dynamics model. Ref. [28]

systems have been proposed one after another. Its randomness

established the maximum benefit model under the RPS. Ref.

and volatility will bring challenges to the traditional

[29] constructed the economic dispatch model under RPS, but

dispatching system, the dispatching model becomes more and

only wind power grid connection is considered.

more complex. So some algorithms for solving multi- objective non-linear, non-convex,

To meet the requirements of RPS while taking into account

strongly coupled and complex

the economic environment and other factors can reduce the cost

constrained models are widely used, like hybrid biogeography

and emissions, which is the key to effectively promote the

based optimization with brain storm optimization

[15]

, shuffle

frog leaping algorithm and particle swarm optimization efficient fitness-based differential evolution [18]

algorithm

[17]

[16]

, an

, Squirrel search

.

reform of the national electric power system and develop new energy science. This is the motivation for the research presented in this paper. The novel contributions are as follows: 1) Considering the requirements of the third edition of

Simultaneously, the mathematical model about power

RPS and Assessment Measures (Draft for Comments)[29], we

system dispatching involving renewable energy has become the

propose a short-term DEED model for hybrid renewable

subject of interest for scholars in China and abroad. Ref. [19]

energy system including wind, solar, and hydro (cascade and

proposed an economic dispatch model by thermal - wind which

runoff) power based on the TGC (HDEED-TGC).

included the impact of the estimated and underestimated wind

2) A

new

multi-objective

moth-flame

optimization

power forecasts on the cost, proposed a new bat algorithm with

(MOMFO) is proposed based on the characteristics of the

black hole model and chaotic search to solve the economic

HDEED-TGC.

dispatch problem. But the effect of non-convex characteristics due to valve point effect is not considered. Ref. [20] integrates

3) A heuristic dynamic relaxation method is proposed, which is used to handle the complex constraints of the model.

wind and pumped storage with thermal units to established the

4) Using MOMFO to solve the model for HDEED-TGC,

optimal economic dispatch model, and the comparison results

the optimal generator output and the number of TGC

show that the coordinated systems effectively improves the

transactions under the green certificate trading mode are

profits and decreases emission. Ref. [21] established a model

obtained. It also explores the extent of the penetration of

for

probability

renewable energy into the power grid without using the

distributions, and a two-stage compensation algorithm is used

transaction price of TGC, which has a certain reference value

to solve the model. However, the overall objective function

for refining the renewable energy quota system.

dynamic

multi-wind

plants

based

on

This paper is organized as follows: Section 2 presents the

solar and hydro (cascade and runoff) power units. The single

proposed HDEED-TGC model. Section 3 describes the

line diagram of the adapted network is depicted in Fig. 1. The

formulation of stochastic characteristics of the renewable

network has four thermal generators at buses 30, 31, 32, 33.

sources. Section 4 introduces the moth-flame algorithm and the

The hydropower output usually includes runoff and cascade

proposed MOMFO used to solve HDEED-TGC model. Section

power stations. The power station of runoff is related to the

5 discusses the case studies and simulation results, and

flow called small hydropower, as it only provides several

conclusions are given in Section 6.

megawatts of power. The output of cascade power station is

2. Modelling of HDEED-TGC

usually related to storage capacity and drainage. A cascade and

In the model, the uncertainties of renewable energy

a small hydropower unit as connected at bus 37, the remaining

output are taken into account, and the output of thermal units

cascade units are connected to buses 34, 35 and 36. The PV

and renewable energy are reasonably arranged to optimally

unit supplies power to bus 38 while the output of the wind farm

coordination. Complete the RPS issued by the State Energy

is connected to bus 39. The number and rated power of each

Administration while giving consideration to maximizing

unit are shown in Table 1. Next, we formulate the

economic profits and environment. For the purpose of our study,

HDEED-TGC optimization problem, including the objective

a standard IEEE 39-bus system is adapted to include wind,

functions and the constraints involved.

G 30

37 26 25

28

29

2 27

38

1 18

3

17 21

16

39 15

G 4

24

36

13

5 6

9

G

Thermal unit

14

Cascade hydropower station

23

12 19 11

7 31

20

22

Solar PV arry

10 34 8

33

35

32 G

G

G

G

Small-hydro generaator

G

Wind generator

Fig.1 Adapted IEEE 39-bus system Table 1. Rated Power of Each Unite for Modifed IEEE 39-bus System Thermal plant (bus 30、31、32、33) No. of unit

Rated power of each unit

Hydro power (bus 34、35、36、37） No. of unit

PGR1=PGR2=470MW 4

Rated power of each unit

Solar PV plant(bus 38)

No. of turbines

Total rated power

300MW

100

300MW

PHR=500MW; 5

PGR3=340MW, PGR4=300MW

Wind power plant (bus 39)

Rated power PSR

PRR=5MW

2.1 Objective functions There are two main objective functions, which are described as follows: 2.1.1 Minimization of cost

power stations, and the TGC cost are shown in (1). T

(

)

min FC = ∑ F ( PGi ,t ) + F ( PWi ,t )t + F ( PSi ,t )t + F ( PHi ,t )t + F ( PRi ,t )t + FTGC,t (1) t =1

where T is the total scheduling period; t is time interval; F(PG)t

The first objective is the minimization of cost in the

is the fuel cost of thermal generators at the time-instant t;

dispatching cycle. It includes the fuel cost of running the

F(PWi,t), F(PSi,t), F(PHi,t), F(PRi,t) are the operating cost of wind,

thermal generators, operating costs of wind, solar and hydro

solar, cascade and small-hydro power at t, respectively; FTGC,t

is the cost of TGC at t.

ith reservoir at time t. The output of the cascade generators is a

The fuel cost of a thermal generator follows a quadratic relationship with the generated output power. Considering the valve point loading effect, a pulse will be superimposed on the energy consumption curve of the generator[30] It can be expressed as follows: NG

2 2 N H C V  1i Hi ,t + C2i QHi ,t + C3iVHi ,t QHi ,t  F ( PHi ,t ) = FHi ⋅ ∑   i =1  + C4 iVH i ,t + C5i QH i ,t + C6 i  

(6)

where C1i - C6i are power generation coefficients of ith unit; FHi

{

F ( PG i ,t ) = ∑ ai + bi PGi ,t + ci PGi ,t 2 + d i sin ei ( PGmin − PGi ,t  i i =1

function of drainage and storage capacity given by (6)[31].

}

(2)

represent the direct cost coefficient for ith unit. The TGC includes a mechanism for purchase and sale.

where, NG is the numbers of thermal generators; ai bi ci di

According to the requirements of the third edition of RPS and

ei are fuel cost coefficients of the ith thermal generator. PGi,t is

Assessment

active power of the ith generator at t; PGi,t

min

is minimum of

active output.

Measures

(Draft

for

Comments)[29],

the

transactions involving the TGC are divided into non-water renewable energy and hydropower energy. To promote the

The operation cost of renewable energy power generation

diversification of renewable energy, a certificate issuance

consists of three components: a direct cost associated with

principle according to different types of renewable energy

generating scheduled power, a penalty cost for underestimation

technology is proposed in this paper, a unit of TGC can be

and a reserve cost for overestimation. The direct cost is borne

obtained for 1 MWh of wind or hydropower and 2 unit of TGC

by the entity generating the renewable energy, e.g. owner of a

for 1 MWh solar power. So the total transaction cost of the

wind farm, etc. Generally, the renewable power energy

TGC is as follows:

generating entity is owned by a private operator, and the system

FTGC,t = M TGC (QW-TGC,t + QH-TGC,t ) + 2M TGC (QS-TGC,t )

operator has to pay the owner for the power. The penalty and

where MTGC represent the transaction cost renewable energies;

reserve costs are related by a linear difference between the

Q represent the transaction quantity corresponding to different

available and scheduled power. So the operation cost of wind,

types of renewable energies are defined in the beginning.

solar and small-hydro power can be expressed as follows:

2.1.2 Minimum emission

Eq.(3)-(5), the corresponding symbols are defined in the beginning of this article.

(7)

The second objective is the minimization of harmful gases emitted from only the fossil energy-based thermal unit.

NW

NW

NW

For the sake of generality, the emission characteristics of all

i =1

i =1

i =1

harmful gases are considered numerically. The functional

F ( PWi ,t ) = ∑ FWi ( PWi ,t ) + ∑ k rWi ( PWi ,tur ,t −PWi ,t ) + ∑ kpWi ( PWi ,t − PWi ,tur ,t ) (3) NS

NS

NS

relationship between the emission of harmful gases and the

i =1

i =1

i =1

active output of thermal units can be approximated as a sum of

F ( PSi ,t ) = ∑ FSi ( PSi .t ) + ∑ kpSi ( PSi ,tur ,t − PSi ,t ) + ∑ krSi ( PSi ,t − PSi ,tur ,t ) (4) NR

NR

NR

i =1

i =1

i =1

F ( PR i ,t ) = ∑ FRi ( PRi .t ) + ∑ kpR i ( PR i ,tur ,t − PRi ,t ) + ∑ krR i ( PR i ,t − PRi ,tur ,t ) (5) IH1,t

T

t =1 i =1

PH2,t Reservoir 2 VH2,t

Reservoir 1 VH1,t QH1,t

QH2,t

IH3,t Reservoir 3 VH3,t

NG

{

}

FE = ∑∑ (α i + β i PGi ,t + γ i PG2i ,t ) + ξ i exp(λi PGi ,t )

IH2,t

PH1,t

quadratic and exponential functions, as shown in (8).

PH3,t

(8)

where αi, βi, γi, ζi, λi represent the emission characteristic coefficient of the ith generator. 2.2 Constraints The optimization model must satisfy different constraints,

QH3,t IH4,t

PH4,t

which are described in the following.

Reservoir 4 VH4,t QH4,t

Fig. 2 The analysis of energy conversion at time interval t

The cascade hydropower has three-stage stations with four cascade generators. The analysis of energy conversion at a time interval t is shown in Fig. 2, where IHi,t, QHi,t, VHi,t are inflow volume, water discharge and storage volume from the

1)The power balance constraint in (9), the corresponding symbols are defined in the beginning of this article. NW NS NH  NG  ∑ PGi ,t + ∑ PWi ,t + ∑ PSi ,t + ∑ PHi ,t +   i=1  i =1 i =1 i =1  NR  = Pload,t + Ploss,t  P +Q  + Q + Q ∑ R W-TGC, t H-TGC, t S-TGC, t  i=1 i ,t 

(9)

The ramp rate limits for each thermal generator is shown in

(10), where DRi,t, URi,t are up / down ramp limits of the ith

made up of inflow volume, water discharge, storage volume at

thermal generator at t.

time t and water discharge of upstream reservoirs. Rui is the

DRi ,t ≤ PGi ,t − PG i ,t −1 ≤ U R i .t

(10)

number of upstream reservoirs of the ith reservoir; τni represent

2)To maintain the stability of operation, the active output

the delay of water delivery from reservoir nth to ith.

of each generator must be limited. The output / planned power

4)The TGC transactions must satisfy a quota ratio within 24 hours of the total dispatching cycle, including the total

of each generators must be within their limits. PGmin ≤ PG i ,t ≤ PGmax i i

(11)

PWmin ≤ PWi ,t ≤ PWmax i i

(12)

renewable energy and the non-water energy quota ratio. These quotas are set by the government ,so the constraints for them are given in (19) and (20). NS NH NR  NW  T ∑ PWi ,t + ∑ PSi ,t + ∑ PHi ,t + ∑ PR i ,t   i=1  ≥ Index ⋅ ∑ Ploadt i =1 i =1 i =1 ∑ t =1  t =1  +QW-TGC,t + QH-TGC,t + QS-TGC,t  T

PSmin ≤ PSi ,t ≤ PSmax i i

(13)

PHmin ≤ PHi ,t ≤ PHmax i ,t i ,t

(14)

PRmin ≤ PR i ,t ≤ PRmax i i

NS  NW  T ∑ PWi ,t + ∑ PSi ,t +   i=1  ≥ Index nh ⋅ ∑ Ploadt i =1 ∑ t =1  t =1  QW-TGC,t + QS-TGC,t 

(19)

T

(20)

(15) where Index, Indexnh respect the total RPS and non-water ratio.

3)Hydropower constraints consist of reservoir capacity

2.3 Dynamic environmental economic dispatch model based

constraint(16), drainage constraints(17) and water dynamic

on the TGC

balance limits. The water dynamic balance limits are simplified

Combining the objective functions and constraints given

in the form of mathematical constraints without considering

earlier, the structure of HDEED-TGC is shown in Fig. 3. Eq.

water spillage of hydro plant in (18). The initial and final

(21) shows the decision variable xi for the optimization model.

storage capacity are set to Vi,ini and Vi,end.

This variable is a matrix, which made up of five individual

VHmin ≤ VHi ,t ≤ VHmax i i

(16)

QHmin ≤ QHi ,t ≤ QHmax i i

(17)

matrices and three column vectors. Each matrix is indicated by a set of columns connected by dotted lines. These matrices and columns represent the following decision variables: PGi,t, PWi,t, PSi,t, QHi,t, PRi,t, QW-TGC,t, QH-TGC,t, QS-TGC,t.

Rui

VHi ,( t +1) = VHi ,t + I Hi ,t − QHi ,t + ∑ (QH n ,( t −τ ) )

(18)

ni

i =1

In (18) the water dynamic balance limits for each reservoir is  PG1,1   PG xi =  1,2  M P  G1,T

L PGN L PGN

G ,1

PR1,1

L PR N

G ,2

PR1,2

L PR N

M PR1,T

L M L PR N

L M L PGN

G ,T

Thermal sub-system

R ,1

QH1,1

L QH N

R ,2

QH1,2

L QH N

M

L M L QH N

R ,T

QH1,T

Hydro sub-system

H ,2

L PWN

M

L M L PWN

PW1,T

H ,T

Wind sub-system

L PSN ,1

QW-TGC,1

QH-TGC,1

W ,2

L PSN

QW-TGC,2

QH-TGC,2

L M L PSN ,T

M

M

W ,T

M PS1,T

QW-TGC,T

QH-TGC,T

S

S ,2

S

PV sub-system PS1

QS-TGC,1   QS-TGC,2   M  QS-TGC,T  

(21)

TGC QW-TGC

PV 1

Wind plant 1 PV 2

QH-TGC

PS2

QH-TGC

…

PH1

PG(NG)

Thermal plant NG

PS1,1 PS1,2

PW1

PW(NW)

Reservoir1

Fuel

W ,1

QW-TGC

Small-hyro PR(NR)

…

…

Thermal plant 2

L PWN

…

Thermal plant 1 PG2 Fuel

PW1,1 PW1,2

PR1

PG1

Fuel

H ,1

PH(NH)

QS-TGC

PS(NS)

Wind plant NW

Reservoir2

QS-TGC

Sell

PV 3

Buy

DEED-RQ system MOMFO algorithm

Dynamic scheduling scheme

constraint

Fig.3 Schematic outline for the scheduling system

3. Random model of renewable sources The output of wind, solar and small hydropower station

have a random behavior, each are dependent upon wind speed, solar radiation and river discharge. In this section, the

probability density functions (PDF) of Weibull distribution[22], Beta distributio

[23]

and Gumbel distribution

[23]

are used to

describe the wind speed, solar radiation and river flow velocity. Consequently, the PDF of each renewable energy output are obtained by the relationship between wind speed, solar

Beta distribution can be used to describe the distribution of random solar radiation[23]. Its PDF is expressed in (27) fSIt (SIt ) = ( Γ(at + bt ) / ( SImax ⋅Γ(at )Γ(bt )) ) (SIt / SImax )at −1(1 − SIt / SImax )bt −1 at = SIt

2

((1− SI ) / σ t

2 t

)

−1/ SIt , bt = at (1/ SIt −1)

(27)

radiation, river flow speed and output power. The linear

The energy conversion relationship between the PV output

difference between available / scheduled power of renewable

power and solar radiation is given by (28)[22].

cumulative probability function of output power. Symbols in

2  P ⋅ ( SI t / SI max SI C ), 0 < SI t < SI C PSt ( SI t ) =  SR  PSR ⋅ ( SI t / SI max ), SI t ≥ SI C

formulas are dealt with in this section are defined in the

where, based on the analysis of sample data, SImax is the

beginning of this article.

maximum solar radiation set as 1000 W/m2 and SIC is the

3.1 Hourly characteristics of wind source

fixed-point set as 120 W/m2. The PDF of PV output is shown in

energy as outlined in Section 2 are expressed in the form of a

Weibull distribution is the most suitable model for describing wind speed distribution, the PDF is as follows: fWS ,t (WSt ) = ( kt / ct ) (WSt / ct )

kt −1

exp( − (WSt / ct ) )

expressed in (30) (31). (22)

Assuming the relationship between output power and wind speed of wind turbines to be a piecewise linear function, the output can be expressed as shown in (23)[22]. In this paper, we consider Enercon E82-E4 turbines, whose cut-in, cut-off and rated wind speeds are 3 m/s, 25 m/s and 16 m/s.  0, WS t < WS in ,WS t > WS off  WS  t − WS in , WS in ≤ WS t ≤ WS R PWt (WS t ) =  PWR WS R − WS in   PWR , WS R ≤ WS t ≤ WS off 

(29). Hence, the linear difference between the available and scheduled power of the PV as outlined in Section 2 can be

kt

kt = (σ t / WS t ) −1.086 , ct = WS t / (Γ (1 + 1 / kt ))

PSt SI C at2-1  1 Γ(at + bt ) ⋅  SI 2 ⋅ Γ ( a ) Γ ( b ) P t t SR SI max  max  PSt SI C 12 bt −1 SI maxSIC 12 P SI × 1 − , 0 < PSt < SR C PSR SI max PSt PSR SI max  f PS ( PSt ) =  t Γ(at + bt ) PSt at −1 (29)  1  SI ⋅ Γ(a )Γ(b ) ( P ) max SR t t   PSt bt −1 SI max P SI ) , PSt ≥ SR C ×(1 − P P SImax SR SR 

(23)

speed, the PDF of unit output power is also discontinuous as can be observed in (24). Hence, the linear difference between the available and scheduled power of the wind turbine as outlined in Section 2 can be expressed in (25)(26), where Pr(.) denotes the probability of events in parentheses. WSoff kt WSin kt  Pr( PWt = 0) = 1 − exp(−( c ) ) + exp(−( c ) ) t t   WSoff kt WSR kt ) ) − exp(−( ) ) Pr( PWt = P WR ) = exp(−( ct ct  f PW,t ( PWt ) =  kt −1  kt (WSR − WSin ) WS + PWt (WS − WS )   in R in  kt  ct * P WR  (24) P WR   WS + ( P / P ) * ( WS − WS )    in Wt WR R in )kt −1  else *exp  −( ct    PWRi

(PWt − PWi,t ) fPW (PWt )d(PWt ) + Pr(PWt = P WR )*(P WR −PWi,t ) (25)

PWi ,t

t

( PWi ,t − PWi ,tur ,t ) = ∫

PWi ,t

0

( PWi ,t − PW t ) f PW ( PWt )d ( PWt ) + Pr( PWt = 0) * PWi ,t (26) t

( PSi ,tur ,t − PSi ,t ) = ∫

PSR i

( PSi ,t − PSi ,tur ,t ) = ∫

PSi ,t

PSi ,t

0

As the output wind power is a discontinuous function of wind

(PWi,tur,t − PWi ,t ) = ∫

( PSt − PSi ,t ) f PS ( PSt )d ( PSt )

(30)

( PSi ,t − PSt ) f PS ( PSt )d ( PSt )

(31)

t

t

3.3 Hourly characteristics of small-hydro source The distribution of river velocity can be described by Gumbel distribution[22], the PDF is shown in (32).  1  RI − λ   RI t − λ   f RIt ( RI t ) =   exp  t  exp  − exp   γ   γ   γ  

(32)

The output power of a small hydropower generator is highly correlated with river flow velocity and head, Eq. (33) is a simplification of the relationship between them, where η is the efficiency of turbine-generator assembly set as 0.85, ρ is the water density set as 1000 kg/m3, H is effective pressure head set as 25 m and g is the acceleration due to gravity equal to 9.81 m/s2. PR t ( RI t ) = ηρ gRI t H

(33)

The PDF of the small hydro output PRt can be shown in (34), where h=ηρgH. f PR ( PR t ) = t

3.2 Hourly characteristics of solar source

(28)

PR − hλ PR − hλ 1 exp( t )exp(− exp( t )) hγ hγ hγ

(34)

Therefore, the linear difference between the available /

scheduled power can be expressed in (35)(36). ( PR i ,tur ,t − PR i ,t ) = ∫

PRR i

PR i ,t

( PR i ,t − PRi ,tur ,t ) = ∫

PR i ,t

0

Gumbel distribution describing are approximately regarded as constant in each dispatching period. These parameter are λ = 15

( PR t − PR i ,t ) f PR ( PR t )d ( PRt )

(35)

( PRi ,t − PR t f PR ( PR t )d ( PR t )

(36)

t

and γ = 1.2. Moreover 8760 samples of solar radiation and wind speed data in a region were selected for the whole year[32].

t

In this area the average daily temperature of the whole year is about 25

3.4 Selection of parameters of each PDF

. The control parameters kt and ct of the Weibull

The river velocity is limited by regional conditions. It is

distribution, and at, and bt of the Beta distribution obtained

closely related to the length of the river, the number of

after running Monte-Carlo scenarios and each parameter are

tributaries, also affected by climatic conditions. However, these

given in Table 2. The PDFs of Weibull and Beta distributions

conditions within a total dispatching period of 24 h is not

fitted to each dispatching period are shown in Figs. 4 and 5.

change quite obviously. Therefore, the control parameters in 0.08 0:00:00 1:00:00 2:00:00 3:00:00 4:00:00 5:00:00 6:00:00 7:00:00 8:00:00 9:00:00 10:00:00 11:00:00 12:00:00 13:00:00 14:00:00 15:00:00 16:00:00 17:00:00 18:00:00 19:00:00 20:00:00 21:00:00 22:00:00 23:00:00

0.07

0.06

0.05

0.04

0.03

0.02

0.01

0 0

5

10

15

20

25

30

Wind seed (WS,m/s) distribution for the windfram in different time periods

Fig.5 The PDF of solar irradiance in each period

Fig.4 The PDF of wind speed in each period

Table 2 The control parameters of Weibull and Beta in each period Time(h)

0

1

2

3

4

5

6

7

8

9

10

11

Wbl-k

1.474

1.476

1.501

1.539

1.554

1.586

1.580

1.600

1.654

1.745

1.819

1.854

Wbl-c

10.018

10.050

10.202

10.366

10.448

10.669

10.676

10.933

11.394

11.789

11.970

12.045

Bata-a

0.000

0.000

0.000

0.000

0.000

0.000

0.000

4.985

7.485

6.542

5.882

4.061 1.238

Beta-b

0.000

0.000

0.000

0.000

0.000

0.000

0.000

33.808

16.318

6.426

3.054

Time(h)

12

13

14

15

16

17

18

19

20

21

22

23

Wbl-k

1.874

1.911

1.956

1.954

1.914

1.784

1.559

1.403

1.327

1.329

1.373

1.445

Wbl-c

11.984

12.013

12.085

12.004

11.707

11.156

10.382

9.838

9.510

9.461

9.654

9.954

Bata-a

3.382

3.268

3.804

4.066

3.804

2.472

0.000

0.000

0.000

0.000

0.000

0.000

Beta-b

1.085

0.931

1.668

3.388

6.853

12.188

0.000

0.000

0.000

0.000

0.000

0.000

4. Multi-objective Moth-flame Optimization

positioning around a flame (37), which simulates the helical

For multi-objective optimization problems, each objective

flight path of moths around the flame, where Mi is the position

influences and restricts the others. Therefore, our aim is to find

of the ith moth represent the decision variable to be optimized;

a set of solutions (Pareto optimal solutions) that are optimal for

Fj is the jth flame represent the best current position, N is the

each objective. So it can be described as obtaining a set of

current number of iterations, Di represents the distance between

Pareto optimal solutions with sufficient diversity and uniform

the ith moth and its corresponding jth flame as shown in (38), b

distribution[33].

represents the logarithmic helix shape constant and t is a

4.1 Standard Moth-flame Optimization

random number between [r, 1], indicating the distance from the

The Moth-flame Optimization (MFO) was proposed in

next generation of Mi to its corresponding flame (t = -1, the

[34], which is based on the lateral orientation of a moth at a

nearest, while t = 1, the farthest), where r is a convergence

fixed angle to the moonlight at night. When the moth comes

constant and decreases linearly with

across artificial lights, it will fly in a spiral direction toward the flame, inspired by this phenomenon the MFO is appeared. In MFO, each moth searches for a better position by transverse

N between [-2, -1].

M iN +1 = Di ⋅ ebt ⋅ cos(2π t ) + Fj

(37)

Di = |Fj − M i |

(38)

The number of flames is reduced by the N, shown in (39). flame.no = round (flame.max − N *

flame.max − 1 ) Gmax

population[36]. According to this strategy, the position of ith moth first updated by (37) and updated using the Lévy flight

(39)

strategy once again, which is shown in (43).

where Gmax is the maximum number of iterations, flame.max is

M iN +1 = M iN + α ∗ L(λ ) ∗ M iN -best , i = 1,2,… dim

the flames’. When the number of moths is less than the flames,

where, * is point-to-point multiplication; а is step size control

the ith moth follows the ith flame, otherwise, the additional

factor set as 0.01; best is the optimal position of moth; L(λ) is

moths follow the last flame.

the PDF of Lévy distribution is shown in (44), where λ is a

4.2 Proposed multi-objective moth-flame optimization

constant between [1,3].

We propose a new algorithm to further improve the effectiveness and competitiveness. The novel features of it are

L(λ ) =

as follows : •An average entropy strategy to initialize. •A Lévy flight strategy to update the moth position.

Γ(1 + λ *sin(

πλ

(43)

1

λ

)

2 λ −1 ( ) 1+ λ Γ( )*λ *2 2 2

(44)

4.2.3 Variable scale chaotic

•A variable-scale chaotic strategy .

Using the randomness and ergodicity of chaotic variables,

•Based on external elite and Pareto dominance

a chaotic mutation strategy is applied to the variable best in

•Crowding entropy.

(43). The logistic chaotic mapping is shown in (45)[37], where

•Modifying the selection dominant body.

µ=4, logistic.max is the maximum number of iterations of the

4.2.1 Initialization based on average entropy

chaotic search. The steps are as follows:

A uniformly distributed initial population in the search

Step 1: Set xi1 as a random number between [0, 1].

space can be obtained using the average entropy population

Step 2: xi2 is derived using (45).

initialization strategy given in [30]. Considering the initial

Step 3: Variation scale of the chaotic search is derived by

population to have NP individuals and the dimension of

(46), where Mimax / Mimin are upper / lower limits of the ith .

population is dim, according to the information theory, the

Step 4: Calculate the latest position of the ith moth using

group entropy H is the sum of individual entropies, as shown in

(47). If the fitness value of the latest position is better than that

(40). In (41), m is the number of known individuals, k is the

before, the update on dimension i is retained. The process is

new individual and Pik is the degree of similarity between any

repeated until n reaches logistic.max.

two-dimensional variables. In (42), Aj and Bj are upper and lower bounds of the jth variable. H=

1 dim ∑ H j , i = 1 2 ……dim dim j =1

xin+1 = µ xin (1 − x in ) n = 1,2,…,logistic.max

(45)

L = M imin + xin+1 ( M imax − M imin )

(46)

(40) M in+1 = (1 − λg ) M in + λg ⋅ L

m

1 Hj = ∑ (− Pik log Pik ) m + 1 i =1

(

Pik = 1 − (|x j (i ) − xk (i )|) / ( B j − A j )

(41)

)

(42)

λg = 1 − ( ( N − 1) / N )

β

(47)

4.2.4 External elite archives The proposed MOMFO uses the elite external archive strategy to preserve the non-dominant solutions found during

First, m known individuals are randomly generated, then

the evolution process. In the evolutionary process, there are

random generation of a new individual,

three situations described as follows:

calculated the

average entropy by (40)-(42), If it is greater than the set threshold (0.096), the individual is added to the initial population, otherwise abandoned. This process in repeated until

1)The test individual is dominated by any individual in the external archive and refuses to enter the external archive. 2)When an individual dominates an individual in an

NP individuals are obtained.

external archive, the individual is removed from the external

4.2.2 Lévy flight strategy

archive.

Lévy flight strategy was proposed by Paul Lévy, a French mathematician

[35]

. Shlesinger MF applied it to the optimization

of intelligent algorithms, thus increasing the diversity of

3)If the test individual and the individual in the external archive do not dominate each other, the test individual enters the external archive.

When the maximum capacity of external archives is reached,

the

diversity preservation

strategy based

its two adjacent for the jth objective function. The calculation

on

of congestion entropy is shown in (49).

congestion entropy is adopted to tailor Pareto solution. This strategy is described in the next sub-section.

Deb et al. used congestion distance to tailor Pareto solution set in NSGA-II

cij Eij f jmax − f i min

k

dlij log 2 ( plij ) + duij log 2 ( puij )

j =1

f jmax − f i min

= −∑

(49)

where fjmax / fjmin are the upper / lower bounds of the jth

4.2.5 Diversity preservation based on congestion entropy [31]

k

j =1

CEi = ∑

objective function, and k is the number of objective functions.

. However, this method does not

4.2.6 Selection of dominant body

consider the distribution between adjacent and cannot reflect

Based on Pareto dominance theory, the relationship

the congestion between different solutions. The congestion is

between two bodies can be classified as follow, when not

evaluated by combining congestion distance and distribution

dominate each other, the highest crowding degree is chosen.

entropy into a congestion entropy, the definition of it is as

1) xi is dominant over ui expressed as xi f ui .

follows:

2) ui is dominant over xi expressed as ui f xi . 3) None of them is dominant expressed as ui pxi ∧ xi pui .

Eij = −[ plij log 2 ( pl ij ) + puij log 2 ( puij )] plij =

dlij cij

, puij =

duij cij

(48)

, cij = dlij + duij

4.2.7 Algorithm implementation The implementation of the MOMFO is shown in the

where dlij and duij are the distance between the ith solution and

following pseudo-code.

MOMFO-Algorithm flow Input

N Gmax dim ub

lb(Variable range)

Output M(final population) PF(Pareto front) 1.

M←Population initialization uing Mean Entropy(Eq.(40)-(42))

2.

Update external elite archive A

3.

While termination criterion not fulfilled do

4.

{Boundary treatment

5.

Update flame no. sing Eq. (39), OM=fitnessfuncion(M);

6.

If iteration==1 {F=sort(M); OF=sort(OM);}

7.

Else {F=sort(Miteration-1 M);OF=sort(OMiteration-1

8.

Update M using Eq.(37)

9.

Mbest←Logistoc variable scale chaotic variation using Eq.(45)-(47)

10. 11. 12.

OM);}end if

for each search agent {Update the position of the current search agent using Lévy flight (Eq.(43))} end for Update external elite archive A} end while

13.M=A; FP=fitnessfuncion(M); return M and FP; Table 3 Algorithm settings MOMFO

MOSADDE

NSGA-

SPEA2

MOPSO

NP

100

100

100

100

100

Amax

100

100

100

100

100

flame.max

100

variable scale shrinkage factor β

1

cross distribution index ηc

20

20

variation distribution index ηm

20

20

crossover probability Pc

[0.1 0.9]

0.9

0.9

mutation probability Pm

0 0.5

1/dim

1/dim

mutation rate

0.5

adaptive mesh dimension

250

4.3 Validation of MOMFO algorithm Ten multi-objective classical test functions are selected to

verify the performance of MOMFO, the test group include a two-objective optimization and scalable DTLZ series problems

[38]

. For DTLZ, considering the characteristics of HDEED-TGC,

we

just

consider

in

between the Pareto frontier and the real frontier, the smaller the

three-objective optimization. Three classical multi-objective

GD is, the closer to the real frontier. The △ is used to evaluate

algorithms, NSGA-II [38]

2 (SPEA2)

[38]

performance

[38]

differential

[38]

of

MOMFO

, strength Pareto evolutionary algorithm

and MOPSO

multi-objective (MOSADDE)

the

performance of the algorithm. The GD evaluates the proximity

the breadth and diversity, a smaller value is better. The HVR

, and the improved

evolution

evaluates the volume covered in the target space over a range

algorithm

of [0, 1], and is close to 1 as possible.

were compared with MOMFO. The algorithm

4.3.2 Test result

settings are shown in Table 3. The test functions and their

To avoid randomness, each algorithm is executed 30

parameter settings are shown in Appendix A.

times, using median X-z and quartile difference IQR of the 30

4.3.1 Performance index

results to compare the statistical results as shown in Table 4.

Generational

distance

(GD),

Spread

(

)

and

For a more intuitive comparison of the frontier obtained by

Hypervolume Ratio (HVR) are selected to measure the

MOMFO with the real Pareto frontier, see Appendix B.

Table 4 Contrast results of performance indicators function Fonseca

MOMFO(GD)

NSGA-Ⅱ(GD)

MOSADDE(GD)

SPEA2(GD)

MOPSO(GD)

X-z

IQR

X-z

IQR

X-z

IQR

X-z

IQR

X-z

IQR

2.799E-04

3.874E-05

1.479E-04

1.1E-05

3.116E-04

4.1E-05

2.281E-04

2.6E-05

3.274E-04

5.6E-04

Schaffer

1.829E-05

1.266E-06

2.356E-04

1.8E-05

2.465E-04

4.9E-05

2.752E-04

9.0E-05

2.443E+00

5.4E-04

ZDT1

1.017E-04

8.620E-05

1.645E-04

6.9E-05

2.196E-04

6.6E-05

1.575E-03

2.1E-02

1.692E-02

1.1E-02

ZDT2

9.779E-05

1.282E-05

5.396E-05

9.6E-06

1.687E-04

4.6E-05

4.995E-03

5.4E-03

5.298E-02

4.1E-02

ZDT3

1.747E-04

1.754E-05

3.110E-04

4.1E-05

3.532E-04

1.3E-05

1.701E-03

4.2E-03

4.869E-02

9.0E-03

ZDT4

1.041E-05

3.340E-06

1.027E-04

6.7E-05

5.185E-04

1.4E-04

1.128E-01

6.5E-03

2.097E-01

2.5E-01

ZDT6

2.201E-02

3.863E-02

4.176E-04

1.5E-05

7.943E-03

1.1E-03

1.915E-03

1.6E-04

1.992E-02

5.3E-02

DTLZ2

2.081E-03

5.721E-04

5.901E-04

3.3E-05

1.243E-03

2.2E-04

5.183E-03

3.2E-03

1.319E-02

6.3E-03

DTLZ3

1.098E-03

9.370E-05

7.368E-04

3.7E-05

4.692E-03

1.6E-02

6.092E-02

9.8E-02

3.652E+01

5.4E+01

DTLZ4

6.923E-03

1.019E-03

5.364E-04

4.2E-05

1.200E-03

4.1E-04

3.489E-03

3.9E-03

1.063E-02

3.6E-03

function

MOMFO(spr)

NSGA-Ⅱ(spr)

MOSADDE(spr)

X-z

IQR

X-z

IQR

X-z

IQR

SPEA2(spr) X-z

MOPSO(spr)

IQR

X-z

IQR

Fonseca

1.344E-03

2.296E-04

1.156E-01

1.1E-02

3.769E-01

3.1E-02

1.788E-01

1.3E-02

7.781E-01

3.4E-01

Schaffer

2.480E-03

3.496E-04

1.408E-01

1.6E-02

2.862E-01

2.5E-02

2.741E-01

2.4E-02

7.373E-04

5.3E-02

ZDT1

5.346E-04

5.346E-04

1.325E-01

6.7E-03

5.077E-01

3.0E-02

2.676E-01

1.8E-01

2.949E-01

2.7E-02

ZDT2

5.403E-04

1.193E-04

1.256E-01

1.5E-02

4.831E-01

2.3E-02

4.322E-01

3.2E-01

2.851E-01

1.6E-02

ZDT3

6.968E-04

6.898E-05

4.627E-01

8.8E-03

5.881E-01

4.8E-02

4.754E-01

6.0E-02

6.150E-01

4.1E-03

ZDT4

5.407E-04

9.785E-05

1.169E-01

1.0E-02

3.373E-01

2.3E-02

6.678E-01

6.8E-01

3.145E-01

5.2E-02

ZDT6

5.095E-04

1.734E-02

1.332E-01

1.4E-02

4.946E-01

4.4E-02

2.363E-01

2.8E-01

1.156E+00

1.3E-01

DTLZ2

2.029E-03

3.026E-04

6.034E-01

5.0E-02

8.246E-01

8.7E-02

2.413E-01

5.4E-02

7.868E-01

4.0E-01

DTLZ3

7.656E-03

1.467E-03

5.870E-01

4.3E-02

8.879E-01

7.3E-02

2.188E-01

5.0E-02

DTLZ4

4.183E-03

6.659E-04

1.594E+00

3.6E-01

7.749E-01

8.3E-02

2.413E-01

5.6E-02

function Fonseca

[38]

—————— 1.594E+00

3.7E-01

MOSADDE(HVR)

NSGA-Ⅱ(HVR)

X-z

IQR

X-z

IQR

X-z

IQR

X-z

IQR

X-z

IQR

9.999E-01

2.745E-06

9.980E-01

3.8E-04

9.798E+02

1.0E-03

9.853E-01

9.4E-04

9.641E-01

3.9E-02

MOMFO(HVR)

SPEA2(HVR)

MOPSO(HVR)

Schaffer

9.989E-01

5.645E-05

9.973E-01

5.0E-05

9.969E-01

1.3E-05

9.968E-01

1.9E-04

9.967E-01

1.1E-03

ZDT1

9.999E-01

2.924E-06

9.971E-01

4.6E-04

9.947E-01

3.2E-04

9.914E-01

9.3E-04

9.423E-01

9.4E-02

ZDT2

9.999E-01

5.983E-06

9.951E-01

7.0E-04

9.955E-01

2.4E-04

9.877E-01

4.1E-03

5.616E-01

3.1E-01

ZDT3

9.985E-01

2.557E-05

9.979E-01

3.1E-04

9.971E-01

3.0E-04

9.863E-01

4.6E-03

6.278E-01

6.4E-02

ZDT4

9.999E-01

6.638E-06

9.995E-01

1.0E-05

9.895E-01

1.4E-03

6.212E-01

2.3E-01

6.879E-01

3.3E-01

ZDT6

8.462E-01

1.414E-04

9.994E-01

4.0E-05

8.667E-01

1.8E-02

9.389E-01

7.2E-03

9.496E-01

3.0E-03

DTLZ2

9.798E-01

1.277E-05

9.790E-01

2.4E-03

9.591E-01

1.5E-03

9.961E-01

6.4E-04

9.858E-01

8.0E-03

DTLZ3

9.395E-01

5.940E-05

9.790E-01

1.5E-03

9.277E-01

1.7E-01

8.610E-01

1.6E-01

0.000E+00

0.0E+00

DTLZ4

9.695E-01

4.628E-05

9.794E-01

3.6E-03

9.692E-01

6.5E-01

9.598E-01

6.6E-02

9.319E-01

1.5E-02

Note:—— representation can't converge to Pareto optimal frontier.

From Table 3, according to the GD, compared with the MOSADDE, the MOMFO achieves better results for

two-objective

optimization

except

ZDT6.

For

the

three-objective optimization, the results obtained by the

MOMFO are not as good as those obtained by MOSADDE.

constraint is calculated as ε_nh. Then, at each time, the

Compared with the SPEA2 and MOPSO, MOMFO achieves

relaxation values for PWi,t, PSi,t and QW-TGC,t, QSI-TGC,t. The

better results by 1-2 orders of magnitude. If the values of △ are

distribution of ε_nh according to the relaxation values is shown

compared, the MOMFO achieves better results than the

in (52), where P r-ti is the relaxation of the ith decision

MOSADDE, NSGA-II, SPEA2 and MOPSO by 1-3 orders of

variables at t, P r-ti =Pt ,max − Pt ,i . Each decision variable is adjusted i

magnitude. Thus, the diversity preservation strategy based on

according to the value of

congestion entropy can reflect the congestion more accurately

constraints of RPS are treated as the same, just adjusted the

and improve the diversity of solutions. The hypervolume is a

volume of hydropower TGC transactions.

simultaneous measure of convergence and diversity, like the

N   ∆Pti =  P r-ti / ∑ P r-ti  ε _ nh   i =1

MOPSO cannot converge to DTLZ3, so the corresponding HVR value is 0. The only exception is ZDT6, the convergence

Pt,i given by (52). The total quota

(52)

5.1.3 Handling power balance constraints

of MOMFO in ZDT6 is insufficient, so the HVR value is small;

There are numerous methods to deal with the power

However, the other algorithms also have very small HVR

balance constraints in DEED problem[40][44], but the diversity of

values which can be further proved in Appendix B. In summary,

solutions are lost. We use the following three ways to handle:

MOMFO

has

optimization

advantages problems,

in

solving

especially

multi-objective

for

two-objective

optimization problems.

5. Application Problem

of

1)The dynamic relax-able constraint processing method is adopted to revise the output of thermal units. 2)The minimum violation for constraint is regarded as the

MOMFO

in

HDEED-TGC

best in (33). 3)The degree of violation for constraint in each t is added to the objective function as the third.

5.1 Handling the complex constraints Based on HDEED-TGC contains complex equality and

The steps for the first way described as follows:

inequality constraints, the MOMFO algorithm is improved and

Step 1:Set the range of infeasible solution PGi,t at the

a processing method suitable for handling complex constraints

initial time (only considering generator output constraints),

in the model is added. 5.1.1 Handling the constraints related to cascade hydropower station

PGmin = PGmin , PGmax = PGmax i ,1 i i ,1 i

Step 2: Calculate the violation for balance constraints Dt.,

A two-stage treatment method for water balance

where ε(0) = 1e-6 and gen is the current number of MOMFO

constraints is given in the following: Stage 1: Calculate the violation of reservoir capacity constraints for each time by (50) and revise the drainage using

Stage 2: Reuse (50), select any one time t and using (51) to revise the drainage of each reservoir at t, revise again by (17).

T

((

n =1 t =1

)

QHi ,t = QHi ,t + ∆VHi / T / ∆t

)

ni

t =1

t = 1,2,L,T

by (54), where PG = PGmax − PG i ,t

(50) ∆PGi ,t

(51)

After completing the treatment of water balance constraints, the storage capacity VHi,t may still violate the (16). A heuristic strategy is used to deal with storage constraints

[40]

.

The specific processing method is given in Fig.6. 5.1.2 Renewable Portfolio Standard constraints First, the violation degree of the non-water RPS

(53)

Step 3: Calculate the relaxation values and revised values

T

∆VHi = Vi ,ini − Vi ,end + ∑ I Hi ,t + ∑∑ (QH n ,( t −τ ) ) − ∑ QHi ,t t =1

ε (0) ∗ (1 − ( N / G max )),0 < N < G max ε N =

i ,t

Repeat stage 2 until the constraint is satisfied. Rui

iterations, continue to Step 3, otherwise jump to Step 5. 0, N ≥ G max

(51) at each t, then revise again using (17).

T

at t. If the violation is greater than a minimum value ε (53),

i ,t

and PG = PG − PGmin . i ,t

i ,t

NG    PGi ,t / ∑ PGi ,t  Dt , Dt < 0 i =1   = N  P / G P  D , D > 0 Gi ,t  t t  Gi ,t ∑ i =1  

i ,t

(54)

Step 4: Calculate P’Gi,t as the solution of this t by (55).  PGi ,t + ∆PGi ,t , Dt < 0 PG' i ,t =   PGi ,t − ∆PGi ,t , Dt > 0

(55)

Step 5: Setting t=t+1, the upper/ lower limits of PGi,t are expressed in (56) for the next period t.

PGmax = min[ PGmax ,( PGi ,t −1 + U R i ,t -1 )] i ,t i min G i ,t

P

= max[ P

min Gi

Step 6: Repeat Steps 2-5, until t = T.

,( PGi ,t −1 − DR i ,t -1 )]

(56) Start Calculating V by(18) for each reservoir at each time t=1 VHimin≤VHi,t ≤VHimax

QHi,t-1=QHi,t-1t=t-1

Qi,t2

N VHi,t>VHimax

Y max

Qi,t2=△Qi,t-(QHimin-QHi,t) QHi,t=QHimin; Y

△Qi,t=(VHi -VHi,t) QHi,t=QHi,t-△Qi,t QHi,t+1=QHi,t+1+ Qi,t

QHi,t, QHi,t+1 satisfied (17)

N

QHi,t-1=QHi,t-1t=t-1

N min

△Qi,t=(VHi,t -VHi ) QHi,t=QHi,t+△Qi,t QHi,t+1=QHi,t+1-△Qi,t

Y t=1

N

Y

Y

QHi,t, QHi,t+1 satisfied (17)

Y

Y t=1

QHi,t>QHimax

N t=t+1

N t=t+1 N

N

Y Qi,t1=△Qi,t-(QHi,t-QHimax) QHi,t=QHimax

t>T

t=T N QHi,t+1=QHi,t+1+ Qi,t1 t=t+1

Y

QHi,t
Y Qi,t1=△Qi,t-(QHi,t-QHimax) QHi,t=QHimax

t=t+1

Y

N

Qi,t2=△Qi,t-(QHimin-QHi,t) QHi,t=QHimin;

N

QHi,t
QHi,t>QHimax

Qi,t2

Y

t=T N QHi,t+1=QHi,t+1+ i=i+1

End

Qi,t2

Fig.6 Flow chart of heuristic strategy to deal with storage constraints

5.2 Optimal compromise solution extraction based on fuzzy

Kron loss coefficient method is used to calculate the network

theory

loss in (58).

Starting from reality, we need to find just one solution

NG NG

NG

j =1 i =1

i =1

Ploss = ∑∑ PG i ,t Bij PG j ,t + ∑ Bi 0 PGi ,t + B00

that satisfies the objective functions. The optimal compromise

(58)

solution is extracted using fuzzy mathematics, linear function is

where Bij, Bi0 and B00 are the transmission loss coefficients of

chosen as the membership. Eq. (57) is used to normalize the

the system. The values for consumption, emission characteristic

µxi

coefficient, output power limit and system load of each thermal

values of each objective functions after fuzzification, where

is the xth solution after the fuzzy representation at ith objective

generator in the system can be found in [44].

function, M and Amax are the number of objective functions and

The parameters of the MOMFO algorithm are: NP=55,

non-dominant solutions. The solution that corresponds to the

Gmax=20000, Amax =55, flame.max=55 and chaotic variable

x

maximum µ is the optimal compromise. 

M

 Amax

M

scale shrinkage factor β = 1. Figure 7 depicts the relevant 

situation of using MOMFO to solve, Figure 7(a) shows the

µ x =  ∑ µix /  ∑∑ µi j   i =1 j =1 i =1 



(57)



non-dominant solutions and it can be seen that the Pareto optimal solutions are more uniformly distributed in the Pareto

5.3 Analysis and conclusion Three different cases of the 10-machine system are

front. Figure 7(b) describes the convergence of each objective

considered. This example is simulated on Windows 7PC (2.90

functions with the iterations. The optimal compromise solution

GHz, 4 GB RAM) using MATLAB R2015b.

given in Table 5, the power balance can be checked in each

5.3.1 Case 1:Traditional 10-machine system model

section. To illustrate the competitiveness of MOMFO, the

The example system consists of a dynamic system

results of cost and emission corresponding to the compromise

composed of 10 thermal generators. Considering the valve

solution are compared with other algorithms for the same

point effect, power balance constraints and ramp rate limits, the

model. The results are shown in Table 6.

Table 5 Optimal compromise solutions of Case1 for MOMFO T

TG1

TG2

TG3

TG4

TG5

TG6

TG7

TG8

TG9

TG10

Total

Pload

Ploss

1

208.2868

153.7434

87.9388

134.5217

132.4933

96.1325

81.9522

81.2323

67.5185

12.5178

1056.337

1036

20.3371

2

158.3082

141.3521

121.8625

139.9413

154.0515

82.7073

111.8472

100.1208

80.0000

42.4250

1132.616

1110

22.6160

3 4

150.0000 202.1476

210.5687 203.0931

174.6413 185.0831

136.4403 129.6849

201.8799 205.7610

130.9575 180.9575

81.8472 104.2665

90.2619 115.4170

69.1691 76.6655

41.4902 39.5470

1287.256 1442.623

1258 1406

29.2560 36.6233

5 6

204.1202 242.2816

211.2923 225.9521

177.7308 193.9570

168.1833 191.8853

212.8707 230.3723

218.2201 251.6767

107.3291 122.8028

96.4607 109.0577

71.5773 72.3933

52.6675 36.8802

1520.452 1677.259

1480 1628

40.4520 49.2589

7

249.3148

261.3326

240.7596

203.9144

206.6505

221.0667

130.0000

120.0000

79.7945

43.4654

1756.298

1702

54.2984

8 9

205.2519 273.3921

257.8501 289.9026

292.1582 291.1204

216.8838 265.8409

243.0000 242.2090

245.9158 252.4261

127.5278 128.1525

114.4450 119.3429

76.4081 79.6928

55.0000 51.7164

1834.441 1993.796

1776 1924

58.4406 69.7956

10 11

297.3517 324.4469

326.6217 362.6937

326.8739 336.6223

280.4900 292.7147

240.5834 239.7102

256.0497 259.8485

127.9977 128.3408

117.7472 120.0000

77.6903 76.8974

48.6077 50.3484

2100.013 2191.623

2022 2106

78.0133 85.6230

12 13

337.0188 314.8044

388.6467 356.5369

335.5291 315.3104

297.7982 300.0000

241.6065 237.9913

258.6495 258.0681

129.1151 125.5264

118.8062 115.6256

77.8054 78.9891

54.9203 51.8029

2239.896 2154.655

2150 2072

89.8958 82.6551

14

290.2647

292.6143

312.7499

263.3734

230.9103

232.3070

124.3921

120.0000

72.6646

55.0000

1994.276

1924

70.2763

15 16

236.0304 166.5872

226.9841 186.4861

232.7499 254.1606

269.3432 221.8113

234.0584 202.7807

259.0513 228.3221

128.1340 128.9860

115.1094 86.1781

79.3733 67.4901

53.4730 55.0000

1834.307 1597.802

1776 1554

58.3069 43.8022

17 18

160.7940 185.5098

155.2640 220.0972

189.3953 242.6250

219.1333 181.0250

225.2967 239.5795

236.9516 238.5868

117.2803 128.0609

107.9683 119.5043

67.6519 67.9733

39.5841 53.5850

1519.319 1676.547

1480 1628

39.3193 48.5468

19

226.8256

219.6349

269.3953

240.6524

243.0000

258.4040

121.7053

119.5396

80.0000

55.0000

1834.157

1776

58.1571

20 21

268.6530 274.0615

290.8018 295.3137

322.3916 289.3344

280.8892 276.9121

243.0000 240.4966

259.8661 238.7229

129.4812 130.0000

120.0000 120.0000

79.6562 78.3221

50.4993 50.7471

2045.238 1993.910

1972 1924

73.2384 69.9105

22 23

206.3200 158.8492

218.6988 140.6377

249.0805 177.2168

228.1338 197.5307

220.9515 192.8129

208.4040 202.6989

125.7078 96.5571

109.2698 95.9891

64.5730 73.1345

45.6143 28.4191

1676.753 1363.846

1628 1332

48.7534 31.8460

24

153.4721

141.0788

169.9949

178.7087

173.5693

162.3713

96.2290

79.5722

38.2859

15.9301

1209.212

1184

25.2124

Table 6 Comparisons of fuel cost and emission for optimal compromise solutions of various algorithms

Cost Emssion (1b) 1.95

MOMFO

Pro_NSG A-Ⅱ[30]

EMODEDCH[47]

Pro _PSO[40]

Pro_MOE A/D[46]

IMOEA/D _CH[46]

MAMOD E[45]

IBFA[46]

RCGA/NS GA-Ⅱ[46]

DSPPMOGA-DE

2488134.68 286009.59

2555180.8 299140.86

2530558.9 299124.96

2508637.5 296807.30

2516800.7 297015.14

2514600 298360

2514113 302742

2517117 299037

2522600 309940

2524513.06 288724.82

105

Pareto Optimal Front

105 3.05

106

2.88

Pareto optimized by MOMFO

Extreme Point: Economic Optimum

Convergence process of economy Convergence process of emssion

1.85

2.95

2.44

1.8 Extreme Point: Emission Optimum

1.75 2.43

Emission function

Cost function

Emission (1b)

1.9

2.44

2.45

2.46

2.47

2.48

2.49

2.50

2.51

Cost (＄)

2.52 106

2.4

2.8 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

Iterations

Fig.7.(a) The Pareto front

2 104

Fig.7.(b) The convergence curves Fig.7 The solution process of MOMFO for Case 1

From Table 6 it can be noted that MOMFO reduces the

coefficients of thermal units are shown in Table 7. Table 8

total generation cost and emission. It is obvious that the

shows the direct/ penalty/ reserve cost coefficients of hydro/

MOMFO provides the most satisfactory solution for solving the

solar/ wind power[48]. The initial and final reservoir capacities,

DEED problem in a traditional power system. It has a

reservoir capacity limit and inflow of cascade reservoirs in each

significant advantage in terms of fuel cost or pollution

operation period are the same as those used in [31]. The system

emissions and can provide better scheduling scheme.

load and MOMFO settings are the same as in Case 1.

5.4.2 Case2

Modified 10-machine system without the RPS

In this case, the hybrid renewable energy system under RPS is not implemented. The cost/ emission characteristic Table 7 Cost and emission coefficient of thermal generators Generator

Bus

a

b

c

d

e

α

β

γ

ζ

λ

G1 G2

30 31

786.7988 451.1591

38.5379 46.1591

0.1524 0.1058

450 600

0.041 0.036

103.3908 103.3908

-2.4444 -2.4444

0.0312 0.0312

0.5035 0.5035

0.0207 0.0207

G3 G4

32 33

1049.9977 1243.5311

40.3965 38.3055

0.028 0.0354

320 260

0.028 0.052

300.391 300.391

-4.0695 -4.0695

0.0509 0.0509

0.4968 0.4968

0.0202 0.0202

Table 8 Direct, penalty and reserve cost coefficient for stochastic renewable energies Solar PV plant(bus 38)( /MW)

Small hydro power(bus 37)( /MW) Direct

Penalty

Reserve

Direct

Penalty

Fr=8

Kpr=1.5

Krr=70

Fsi=9

Kpsi=1.5

Wind power plant (bus 39)( /MW)

Reserve

Direct

Penalty

Reserve

Krsi=78

Fw=8

Kpw=1.5

Krw=75

The optimal compromise solution to the generators’

depicts the relevant situation of using MOMFO to solve,

output are shown in Fig. 8, the power balance constraint can be

Figure11(a) shows the non-dominant solutions and it can be

checked in each section. The corresponding drainage and

seen that the Pareto optimal solutions are more uniformly

reservoir capacities of each cascade hydropower station are

distributed. Figure11(b) describes the convergence of each

shown in Figs. 9-10. It can be observed that the output of each

objective functions with the iterations and it can be seen that

units, the discharge and capacity of reservoirs are within the

the solution can converge effectively.To verify the accuracy of

feasible range.

the solution, the optimal compromise solution to the generators’

2500

output by MOMFO and the number of TGC transactions in

2000

each period are shown in Fig.12.

1500

) W M ( 1000 ut tpu O 500

105

1.55

Pareto Optimal Front Pareto optimized by MOMFO

Extreme Point: Economic Optimum

0 1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Time (h)

TG1 HG1 PV

TG2 HG2 WG

TG3 HG1 Pload+Ploss

TG4 HG4+S-h

Fig. 8 Each generators output and verification of power balance

Reservoir1

Reservoir2

1.45

1.4 Extreme Point: Emission Optimum

1.35 1.89

constraint 30

Emission (1b)

1.5

1.9

1.91

1.92

1.93

1.94

1.95

1.96

1.97

Cost ( ＄)

Reservoir3

1.98 106

Fig.11.(a) The Pareto front of MOMFO

Reservoir4

105 3.4

6

10

2.3

25

convergence process of economic convergence process of emission

)h / 20 3 m 4 01 (e 1 5 gr ah sci 1 0 D

2.6

Cost function

Emission fuction

2.1

5 1.9

1.8

0 1

2

3

4

5

6

7

8

9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

T im e (h)

Fig.9 Hourly water discharge rates of the best compromise solution for Case2

1.7

1 0

0.4

Iterations

0.8

1.6

1.2

2

Fig.11.(b) The convergence curves for MOMFO

2 00 1 80

Fig.11 The solution process of MOMFO for Case 3

1 60

2500

)h 1 40 3/ 4m1 20 01 ( 1 00 e m ul 80 o V 60

2000 ) W1500 (M t u pt 1000 u O 500

40 20

R eser vo ir1

R eser vo ir2

R eser vo ir3

R eser v oir4

0

0 1

2

3

4

5

6

7

8

9

10 11 1 2 13 14 15 1 6 17 18 19 2 0 21 22 23 2 4

Time ( h)

Fig.10 Hourly reservoir volumes of the best compromise solution for Case2

5.4.2 Case3: Modified 10-machine system under the RPS This case consists of a short-term DEED model of hybrid

1

3

-500

5

7

9

11

13

15

17

19

21

23

Time (h) TG1 HG1 PV TGC-H

TG2 HG2 WG Pload-Ploss

TG3 HG3 TGC-W

TG4 HG4+S-h TGC-PV

Fig.12 Optimal compromise solutions of Case3

renewable energy system based on the TGC. The transaction

Figure 13 compares the costs and emissions of three Cases.

cost of TGC is $45/MW, and the quotas of total renewable

Figure 14 is the operation cost of thermal/ wind/solar/

energy and non-hydropower’s are set at 25% and 15%, and

hydropower under the optimal compromise solution of Case 3.

others parameters as the same as those in Case 2. Figure 11

Comparing the data in Figure 13, it can be concluded that the

addition of RPS is beneficial to the economic and environmental benefits. The reason is that, for the system without the RPS, the renewable energy will be purchased preferentially under the dual-objective scheduling with minimal pollution emission, due to its non-polluting characteristics. When the system under the RPS, the dispatching system not only meets the requirements of power consumption and daily

Fig.16 PV power planned under different transaction prices of TGC ) 100 W (M 80 s 60 n iot ca 40 sn ar 20 t C 0 G T f -20 o e -40 um l -60 ov la -80 t o -100 T

1

operation, but also meets the quotas of RPS. The TGC is an equivalent converter of renewable energy and the trade of both

5

7

9

11

13

15

17

19

21

23

Time (h) TGC=30

TGC=35

TGC=40

TGC=50

TGC=55

TGC=60

TGC=45

Fig.17 Total volume of TGC transactions under different transaction

renewable energy and TGC promote the consumption of renewable energy.

3

prices of TGC

To study the influence of different transaction prices of the TGC on HDEED-TGC, the prices are varied between 30

300.00

0 0 0 0 ?1 250.00 ?

T otol c o st

T o to l em s sion

and 60 in a step-size of 5. Comparisons of wind/ PV

200.00

output ,and total power of TGC transactions under different

150.00

prices are shown in Figure 15-17. It is evident that when the

100.00

price of TGC falls, wind power and PV plan output line moving 50.00

downward at peak and valley, the proportion of renewable 0.00 C ase1

C ase2

Case3

energy consumption in the total electricity decreases. Conversely, when the price of TGC rises, the planned curve of

Fig.13 Comparison of the results of Case1, Case2 and Case3 90000

Themal

PV

Wind

Cascade-hydro

Small-hydro

200

80000

180

70000

160 140

60000

wind power and PV output rises, leading to an increase in the proportion of renewable energy consumption. This can be explained as follows: when the price of the TGC is low, the

120

ts 50000 o C 40000

100 80

30000

transactions are mainly purchased and the quota limit can be met for a lower consumption of renewable energy source. On

60

20000

40

10000

20

0

0 1

2

3

4

5

6

7

8

9

10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 Time(h)

the other hand, when the price is high, the system will consume a large proportion of renewable energy and try to exceed the quota requirements of electricity sales to obtain economic benefits. Therefore, the introduction of the TGC trading system

Fig.14 Break-up of total cost for Case 3

increases the flexibility of the power grid performance.

350

6. Conclusion

) 300 W (M250 de n na 200 l p re 150 w o p dn 100 i W

Since China is at a critical stage of promoting the reform of the national electric power system and the effective development of new energy science, the policy for renewable

50

TGC=30

TGC=35

TGC=40

TGC=50

TGC=55

TGC=60

TGC=45

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

Time (h)

energy quota system is constantly being updated and elucidated. In this paper, the requirements of the latest edition of RPS and Assessment Measures are considered, a short-term DEED

Fig.15 Wind power planned under different transaction prices of TGC

model of hybrid power supply system is established. Aiming at the high latitude, non-linear, non-convex and complex

350

TGC=30

TGC=35

TGC=40

TGC=50

TGC=55

TGC=60

TGC=45

) 300 W (M250 de n 200 anl p re 150 w o p 100 V P

constraints involved in the model, MOMFO is designed to solve the problem. The performance of MOMFO is validated over ten multi-objective test functions and the test system for a traditional 10-machine model. The result shown that MOMFO has obvious advantages over other algorithms for solving the

50 0 1 2

3 4 5 6

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

Time (h)

DEED. Then MOMFO is further applied in scenarios with and

without the RPS, explored the extent of the penetration of

are also small. Furthermore, when determining the penalty and

renewable energy into the power grid under different scenarios

storage costs of wind, solar and small hydropower, the

and prices of TGC. It an be used as a reference and justification

influence of spinning reserve source on the price is not

for refining the renewable energy quota system.

considered. Based on the third edition of RPS and Assessment

In this paper, it is assumed that the annual climate change of the selected data sample is small, the seasonal characteristics

Measures, these factors mentioned above will be studied in future work including demand-side response constraints.

are not considered, the temperature and light intensity changes

Appendix Appendix A: A test function No.

Function name

Formula

Variable range

Dim

x ∈ [− 4, 4]

n=3

x ∈ [ −10 3 ,10 3 ]

n=1

ZDT1

 f ( x) = x 1  1 min =  f 2 ( x ) = g ( x )[ 1 − x1 / g ( x ) ]  n  g ( x ) = 1 + 9 ( ∑ i = 2 x i ) /( n − 1)

x i ∈ [ 0 ,1]

n=30

ZDT2

 f ( x) = x 1 1  min =  f 2 ( x ) = g ( x )[ 1 − ( x1 / g ( x )) 2 ]  n  g ( x ) = 1 + 9 ( ∑ i = 2 x i ) /( n − 1)

x i ∈ [ 0 ,1]

n=30

ZDT3

 f ( x) = x 1  1 min =  f 2 ( x ) = g ( x )[ 1 − x 1 / g ( x ) − x1 sin( 10 π x1 ) / g ( x )]  n 0 . 25  g ( x ) = 1 + 9 [( ∑ i = 2 x i ) /( n − 1)]

x i ∈ [ 0 ,1]

n=30

ZDT4

 f ( x) = x 1  1 min =  f 2 ( x ) = g ( x )[ 1 − x 1 / g ( x ) ]  n 2  g ( x ) = 1 + 10 ( n − 1) + ∑ i = 2 [ x i − 10 cos( 4 π x i )

ZDT6

 f ( x ) = 1 − exp( − 4 x ) sin 6 ( 6 π x ) 1 1  1 min =  f 2 ( x ) = g ( x )[ 1 − ( f 1 ( x ) / g ( x )) 2 ]  n 0 . 25  g ( x ) = 1 + 9 [ ∑ i = 2 x i /( n − 1)]

1

Fonseca

  n 2   f 1 ( x ) = 1 - exp  − ∑ ( x i − 1 / n )    i =1  min =  n 2   f ( x ) = 1 - exp  − ∑ ( xi + 1 / n )   2  i =1 

2

Schaffer

 f ( x ) = x 2 min =  1  f 2 ( x ) = ( x − 2 ) 2

3

4

5

6

7

8

9

10

DTLZ2

DTLZ3

DTLZ4

 f 1 ( x ) = (1 + g ( X M )) cos( x 1π / 2 ) cos( x 2 π / 2 )  f ( x ) = (1 + g ( X )) cos( x π / 2 ) sin( x π / 2 ) M 1 2  2 min =  f ( x ) = (1 + g ( X M )) sin( x 1π / 2 )  3 g(X M ) = ∑ ( x i − 0 .5 ) 2 xi ∈ X M   f 1 ( x ) = (1 + g ( X M )) cos( x 1π / 2 ) cos( x 2 π / 2 )  f ( x ) = (1 + g ( X )) cos( x π / 2 ) sin( x π / 2 ) M 1 2  2 min =  f 3 ( x ) = (1 + g ( X M )) sin( x 1π / 2 )   g ( X M ) = 100 [ X M + ∑ ( x i − 0 . 5 ) 2 − cos( 20 π ( x 1 − 0 . 5 ))] xi ∈ X M 

 f 1 ( x ) = (1 + g ( X M )) cos( x 1 a π / 2 ) cos( x 2 a π / 2 )  a a  f 2 ( x ) = (1 + g ( X M )) cos( x 1 π / 2 ) sin( x 2 π / 2 ) min =  a  f 3 ( x ) = (1 + g ( X M )) sin( x 1 π / 2 ) g(X ) = ∑ xi∈ X M ( x i − 0 .5 ) 2 M 

x1 ∈ [ 0 ,1] x i ∈ [ -5 ,5 ] i ≠ 1

n=10

x i ∈ [ 0 ,1]

n=10

x i ∈ [ 0 ,1]

n=12 |XM|=10

x i ∈ [ 0 ,1]

n=12 |XM|=10

x i ∈ [ 0 ,1]

|XM|=10 a=100

n=12

Appendix B: Comparison of real frontier and MOMFO frontier

Ture frontie

Ture frontie

MOMFO

MOMFO

Ture frontie Ture frontie MOMFO MOMFO

a. Fonseca

b. Schaffer

c. ZDT1

d. ZDT2 Pareto Optimal Front 8 7

Ture frontie Ture frontie 6

MOMFO

5

MOMFO

4 3 2 1 0 0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

f1

f. ZDT4

e. ZDT3

g. ZDT6 Pareto Optimal Surface

Ture frontie

1 0.8 0.6 0.4 0.2 0 1

1 0.8

0.5

0.6 0.4

f2

h. DTLZ2

0

0.2 0

f1

i. DTLZ3

g. DTLZ4

Acknowledgement

considering wind power penetration[J]. International Journal of Electrical

This work was supported by National Natural Science

Power & Energy Systems, 2013, 44(1):282-292.

Foundation of China (51667020,51666017). Tianshan Cedar

[2]. Elattar E E . Modified harmony search algorithm for combined

Project in Xinjiang Autonomous Region (2017XS02). Natural

economic emission dispatch of microgrid incorporating renewable

Science Projects of Scientific Research Program of Universities

sources[J]. Energy, 2018:S0360544218312027-.

in Xinjiang Autonomous Region (XJEDU2017I002), Key R&D

[3]. Jubril A M , Komolafe O A , Alawode K O . Solving Multi-Obj

Program

ective Economic Dispatch Problem Via Semidefinite Programming[J]. I

Projects

in

Xinjiang

Autonomous

Region

(2016B0201), Opening of Laboratory in Xinjiang Autonomous

EEE Transactions on Power Systems, 2013, 28(3):2056-2064.

Region (2018D04005).

[4]. Elsakaan A A , El-Sehiemy R A , Kaddah S S , et al. An enhan

References

ced moth-flame optimizer for solving non-smooth economic dispatch p

[1]. Mondal S , Bhattacharya A , Nee Dey S H . Multi-objective economic

roblems with emissions[J]. Energy, 2018:S0360544218311538-.

emission load dispatch solution using gravitational search algorithm and

[5]. Mustafa M N , Holanda B U , Moya Rodríguez Jorge Laureano,

et al. A new approach to economic-emission load dispatch using NSG

[18]. M. Basu . Squirrel search algorithm for multi-region combined h

A II. Case study[J]. International Transactions on Electrical Energy Sy

eat and power economic dispatch incorporating renewable energy sourc

stems, 2018:e2626-.

es[J]. Energy, 2019, 182:296-305.

[6]. Nascimento M H R , Marcus Vinícius Alves Nunes, Jorge Laurea

[19]. Liang H , Liu Y , Shen Y , et al. A Hybrid Bat Algorithm for

no Moya Rodríguez, et al. A new solution to the economical load dis

Economic Dispatch with Random Wind Power[J]. IEEE Transactions o

patch of power plants and optimization using differential evolution[J].

n Power Systems, 2018:1-1.

Electrical Engineering, 2016, 99(2):1-11.

[20]. Nazari M E , Ardehali M M . Optimal coordination of renewabl

[7]. Jiang, Xingwen. Dynamic environmental economic dispatch using

e wind and pumped storage with thermal power generation for maxim

multiobjective;differential evolution algorithm with expanded double sel

izing economic profit with considerations for environmental emission b

ection and;adaptive random restart[J]. International Journal of Electrical

ased on newly developed heuristic optimization algorithm[J]. Journal o

Power & Energy Systems, 2013, 49(1):399-407.

f Renewable and Sustainable Energy, 2016, 8(6):065905.

[8]. Nascimento M H R , Marcus Vinícius Alves Nunes, Jorge Laureano

[21]. Xie M , Xiong J , Ke S , et al. Two-Stage Compensation Algor

Moya Rodríguez, et al. A new solution to the economical load dispatch of

ithm for Dynamic Economic Dispatching Considering Copula Correlati

power plants and optimization using differential evolution[J]. Electrical

on of Multi-wind Farms Generation[J]. IEEE Transactions on Sustaina

Engineering, 2016, 99(2):1-11.

ble Energy, 2016, PP(99):1-1.

[9]. Zhu Z , Wang J , Baloch M H . Dynamic economic emission dis

[22]. Partha P. Biswas a ,P.N. Suganthan a , B.Y. Qu b ,etc.Multiobjec

patch using modified NSGA-II: Dynamic Economic Emission Dispatch

tive economic-environmental power dispatch with stochastic wind-solar-

[J]. International Transactions on Electrical Energy Systems, 2016, 26

small hydro power[J].Energy, 2018, 150(2):1039-1057.

(12):2684-2698.

[23]. Yue Yin. Tianqi Liu, Chuan He. Day-ahead stochastic coordinated

[10]. Jorge Laureano Moya Rodríguez. Solution to economic emission

scheduling for thermal-hydro-wind-photovoltaic systems[J]. Energy, 2019,

load dispatch by simulated annealing: case study[J]. Electrical Enginee

187:115944.

ring, 2017:1-13.

[23]. Kumar J V , Chandra P V , Nikhil G , et al. Integration of ren

[11]. Coello C A C, Lechuga M S. MOPSO: a proposal for multiple

ewable energy sources in dynamic economic load dispatch problem usi

objective particle swarm optimization[C]// wcci. IEEE Computer Society,

ng an improved fireworks algorithm[J]. IET Renewable Power Generat

2002:1051-1056.

ion, 2018, 12(9):1004-1011.

[12]. Nebro A J, Durillo J J, Garcia-Nieto J, et al. SMPSO: A new P

[24]. Bin Zhou ,Da Xu, Canbing Li ,etc.Optimal Scheduling of Biogas

SO-based metaheuristic for multi-objective optimization[C]// IEEE, 200

–Solar–Wind Renewable Portfolio for Multicarrier Energy Supplies[J].I

9:66-73.

EEE TRANSACTIONS ON POWER SYSTEMS, 2018, 33(6):6229-623

[13]. Guo C X , Zhan J P , Wu Q H . Dynamic economic emission

9.

dispatch based on group search optimizer with multiple producers[J].

[25]. Zhang Y Z, Zhao X G, Ren L Z, et al. The development of the

Electric Power Systems Research, 2012, 86(4):8-16.

renewable energy power industry under feed-in tariff and renewable

[14]. Yang F, Bai C, Zhang Y. Research on the Value and Implementa

portfolio standard: A case study of China’s wind power industry[J]. Jo

tion Framework of Energy Internet[J]. Proceedings of the Csee, 2015,

urnal of Cleaner Production 2017,

35(14):3495-3502.

[26]. Zhang Y Z, Zhao X G, Ren L Z, et al. The development of Chi

[15]. Guojiang X , Dongyuan S . Hybrid biogeography-based optimizat

na’s biomass power industry under feed-in tariff and renewable portfol

ion with brain storm optimization for non-convex dynamic economic d

io standard: A system dynamics analysis[J]. Energy,2017(139):947-961

ispatch with valve-point effects[J]. Energy, 2018, 157:424-435.

[27]. Xin-Gang Z , Yu-Zhuo Z , Ling-Zhi R , et al. The policy effect

[16]. Narimani H , Razavi S E , Azizivahed A , et al. A multi-objecti

s of feed-in tariff and renewable portfolio standard: A case study of

ve framework for multi-area economic emission dispatch[J]. Energy, 2

China's waste incineration power industry[J]. Waste Management,2017.

018.

[28]. Xin-Gang Z, Xia-Meng L, Bei L, et al Analysis of the influence

168 :1262-1276.

[17]. Xin Shen , Dexuan Zou , Na Duan, Qiang Zhang . An efficient

of renewable energy quota system on the optimization of power supp

fitness-based differential evolution algorithm and a constraint handling

ly structure in China[J]. Renewable Energy Resources, 2013,31(04):1-5.

technique for dynamic economic emission dispatch[J]. Energy, 2019,

[29]. Ji L, Yi Z, Yu-Zhuo Z, et al Energy-saving Economic Dispatch

186:115801.

of Wind Power Grid Connection Based on Renewable Energy Quota

System[J/OL].Power System Technology:1-7[2019-03-10].https://doi.org/1 0.13335/j.1000-3673.pst.2018.2187.http://www.nea.gov.cn/2018-11/15/c_13

heduling Based on Improved PSO Algorithms[J]. Heilongjiang Electric Power, 2017, 39(2):120-124.

7607356.htm

[40]. Hao T. Short-term optimal dispatching method for hydrothermal p

[30]. Zhi-Jian Z, Jie W. Dynamic Environmental Economic Dispatch of

ower plants based on gravity search algorithm[D].Huazhong University

Power System Based on Improved NSGA-II[J]. Electric Power Auto

of Science and Technology. 2014.

mation Equipment, 2017, 37(2):176-183.

[41]. Zhang Y, Gong D, Geng N, et al. Hybrid bare-bones PSO for d

[31]. Basu M . An interactive fuzzy satisfying method based on evolut

ynamic economic dispatch with valve-point effects[J]. Applied Soft Co

ionary programming technique for multiobjective short-term hydrotherm

mputing, 2014, 18(C):248-260.

al scheduling[J]. Electric Power Systems Research, 2004, 69(2-3):277-2

[42]. Zaman M F , Elsayed S M , Ray T , et al. Evolutionary Algori

85.

thms for Dynamic Economic Dispatch Problems[J]. IEEE Transactions

[32]. Xuena A N, Zhang S, Xue L I. Equilibrium Analysis of Oligopo

on Power Systems, 2016, 31(2):1486-1495.

listic Electricity Markets Considering Tradable Green Certificates[J]. A

[43]. Elaiw A M , Xia X , Shehata A M . Hybrid DE-SQP and hybri

utomation of Electric Power Systems, 2017, 41(9):84-89.http://dkasolarc

d PSO-SQP methods for solving dynamic economic emission dispatch

entre.com.au/source/yulara/yulara-1-fixed

problem with valve-point effects[J]. Electric Power Systems Research,

[33]. Wang Y N , Wu L H , Yuan X F . Multi-objective self-adaptive

2013, 103(8):192-200.

differential evolution with elitist archive and crowding entropy-based diversity measure[J]. Soft Computing, 2010, 14(3):193-209.

[44]. Wang Y N , Wu L H , Yuan X F . Multi-objective self-adaptive differential evolution with elitist archive and crowding entropy-based

[34]. Mirjalili S . Moth-flame optimization algorithm: A novel nature-i

diversity measure[J]. Soft Computing, 2010, 14(3):193-209.

nspired heuristic paradigm.[J]. Knowledge-Based Systems, 2015, 89:22

[45]. Yong-Shen Z, Jie W, Bo-Yang Z , et al. Multi-objective dynamic

8-249. [35]. Schreier A L, Grove M. Ranging patterns of hamadryas baboons: random walk analyses[J]. Animal Behaviour, 2010, 80(1):75-87. [36]. Shlesinger M F. Mathematical physics: search research.[J]. Nature, 2006, 443(7109):281-282.

environmental economic dispatch with wind farms[J]. Power System Technology, 2015, 39(5):1315-1322. [46]. Xin-Wen J, Jiann-Zhong Z, Wang H, et al. Modeling and Solutio n of Power System Dynamic Environmental Economic Dispatch[J]. Po wer System Technology, 2013, 37(2):385-391.

[37]. Shi Y, Yu P. On chaos of the logistic maps[J]. Dynamics of Con

[47]. Bai-Hao Q. Research on Power System Dispatching with Wind P

tinuous Discrete & Impulsive Systems, 2007, 14(2):175-195.

ower and Electric Vehicle Based on Evolutionary Computing.[D].zhong

[38]. Liang-Hong W, Yao-Nan W. Dynamic Differential Evolution Algo

yuan university of technology. 2018

rithms and Their Applications[M]. Science Press, 2014.

[48]. Nang-Nang Z, Study of short-term generation of renewable energ

[39]. Guo-Dong W. Research on Dynamic Environmental Economic Sc

y quota under the[D]. Shanghai Jiao Tong University. 2014.








Highlights A new dynamic environmental economic dispatching model is proposed for hybrid energy A multi-objective moth-flame optimization(MOMFO) is proposed for dispatching model MOMFO is effective and competitive compared with others multi-objective algorithms An efficient heuristic dynamic relaxation method is developed for complex constraints The influence of different trading systems on scheduling results in model is analyzed

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: