Stochastic optimal operation of reservoirs based on copula functions

Journal of Hydrology 557 (2018) 265–275

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Research papers

Stochastic optimal operation of reservoirs based on copula functions Xiao-hui Lei a, Qiao-feng Tan b,⇑, Xu Wang a, Hao Wang a, Xin Wen b, Chao Wang a, Jing-wen Zhang c a

State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China c State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China b

a r t i c l e

i n f o

Article history: Received 21 August 2017 Received in revised form 27 October 2017 Accepted 13 December 2017 Available online 14 December 2017 This manuscript was handled by G. Syme, Editor-in-Chief, with the assistance of Bellie Sivakumar, Associate Editor Keywords: Stochastic dynamic programming (SDP) Copula function Reservoir operation Uncertainty Value function

a b s t r a c t Stochastic dynamic programming (SDP) has been widely used to derive operating policies for reservoirs considering streamflow uncertainties. In SDP, there is a need to calculate the transition probability matrix more accurately and efficiently in order to improve the economic benefit of reservoir operation. In this study, we proposed a stochastic optimization model for hydropower generation reservoirs, in which 1) the transition probability matrix was calculated based on copula functions; and 2) the value function of the last period was calculated by stepwise iteration. Firstly, the marginal distribution of stochastic inflow in each period was built and the joint distributions of adjacent periods were obtained using the three members of the Archimedean copulas, based on which the conditional probability formula was derived. Then, the value in the last period was calculated by a simple recursive equation with the proposed stepwise iteration method and the value function was fitted with a linear regression model. These improvements were incorporated into the classic SDP and applied to the case study in Ertan reservoir, China. The results show that the transition probability matrix can be more easily and accurately obtained by the proposed copula function based method than conventional methods based on the observed or synthetic streamflow series, and the reservoir operation benefit can also be increased. Ó 2017 Published by Elsevier B.V.

1. Introduction The reservoir hydropower operation is generally formulated as a nonlinear, stochastic optimization problem. To guide the reservoir operation efficiently, the reservoir operators usually predefine an operating policy in the design stage and use it in the real-time operation. According to the methods used to deal with stochastic inflow, the operating policy can be obtained by 1) implicit stochastic optimization (Chandramouli and Deka, 2005; Chandramouli and Raman, 2001; Mehta and Jain, 2009; Mousavi et al., 2005; Panigrahi and Mujumdar, 2000; Willis et al., 1984; Young, 1967), which can convert a stochastic optimization problem into a deterministic one to operate the reservoir under several equally likely inflow scenarios and then extract the operating rules from the resulting set of optimal operating data; 2) parameterization-simula tion-optimization (Koutsoyiannis and Economou, 2003; Koutsoyiannis et al., 2010; Momtahen and Dariane, 2007; Nalbantis and Koutsoyiannis, 1997; Oliveira and Loucks, 1997; Tan et al., 2017a), which first predefines a shape for the rule curve based on some parameters and then use heuristic strategies to look ⇑ Corresponding author. E-mail address: [email protected] (Q.-f. Tan). https://doi.org/10.1016/j.jhydrol.2017.12.038 0022-1694/Ó 2017 Published by Elsevier B.V.

for the combination of parameters that provide the best reservoir operating performance under possible inflow scenarios; and 3) explicit stochastic optimization (Chou et al., 2016; Harboe, 1993), which converts the inflow with different probabilities directly into the optimization problem. Stochastic dynamic programming (SDP) is one of the most popular explicit stochastic optimization methods (Celeste and Billib, 2009). Previous studies on SDP have focused mainly on the inflow uncertainties and the curse of dimensionality. The classical SDP model (Loucks et al., 1981) does not take into account the forecast information in addressing the inflow uncertainties in reservoir operation. Given the significant effect of the forecast uncertainty on the reservoir operation, it is explicitly incorporated into some reservoir operation models (Karamouz and Vasiliadis, 1992; Mujumdar and Nirmala, 2007; Stedinger et al., 1984; Xu et al., 2014; Kim and Palmer, 1997). Karamouz and Vasiliadis (1992) proposed a Bayesian SDP (BSDP) model, in which a Bayesian approach was incorporated within the SDP formulation. BSDP is different from the classical SDP in choosing state variables. Bayesian Decision Theory can easily incorporate new information by updating prior probabilities to posterior probabilities, thus reducing the effect of natural and forecast uncertainties on the model. Kim and Palmer (1997) proposed a BSDP model to investigate the value

266

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

of seasonal inflow forecasts in hydropower generation, based on which the monthly operating policies were derived for Skagit Hydropower System. Mujumdar and Nirmala (2007) used the BSDP model to derive an operating policy for a multi-reservoir hydropower generation system. Xu et al. (2014) proposed a new twostage BSDP (TS-BSDP) model for real-time operation of cascaded hydropower systems to handle varying uncertainty of flow forecasts from Quantitative Precipitation Forecasts. In the last decades, a number of aggregation-disaggregation models have been proposed to solve the curse of dimensionality. Mujumdar and Nirmala (2007) aggregated individual inflows of reservoirs in a cascaded hydropower system and used a BSDP model to represent the forecast uncertainty of the aggregate inflow. Tejada-Guibert et al. (2010) and Mujumdar and Nirmala (2007) proposed SDP models in which only the inflows of individual reservoirs were aggregated, and the storage of individual reservoirs and the aggregate flow were taken as state variables. Xu et al. (2014) proposed a TS-BSDP model with aggregate flow and storage of individual reservoirs as state variables to reduce the complexity of the optimal problem. The curse of dimensionality can also be avoided by approximating the value function by the iteration method (Cervellera et al., 2006; Drouin et al., 1996; Gal, 1979; Lee and Lee, 2006). Although the previous researches about SDP have made great improvements. There are two aspects the previous study does not pay enough attention. On the one hand, the transition probability matrix is usually based on observed or synthetic streamflow series in which the probability is replaced with the frequency. Consequently, the precision of the transition probability matrix depends critically on the representativeness and size of the sample or the reliability of the synthetic streamflow series. The transition probability matrix can be distorted in the case that the observation occurs by chance or the synthetic streamflow series cannot well characterize the original data. On the other hand, the value function in the last stage is ignored in most previous studies. That is to say, there is no difference in future value no matter where the ending water storage of the last period is. This is obviously conflicting with the actual situation. As we all know, different carryover storages will bring different future benefits. The copula function has been widely used in hydrology and water resources fields due to its ability to set up the joint distribution of multiple variables (Favre et al., 2004; Jenq-Tzong et al., 2010; Karmakar and Simonovic, 2009; Singh and Zhang, 2006; Reddy and Ganguli, 2012; Tan et al., 2017b). In this study, a new method based on copula function to calculate the transition probability matrix was proposed and the results were compared with the conventional statistical method. Meanwhile, the value in the last stage was calculated from a simple recursive equation by stepwise iteration, and the value function was fitted with a linear regression model. These improvements were incorporated into the classic SDP and then applied to guide the hydropower operation of Ertan reservoir, China. The next section introduces the proposed method. A case study is presented in Section 3, and the conclusions are drawn in Section 4.

where f ðst1 ; qt Þ is the maximum expected hydropower generation from period t to the entire planning horizon (T), in which t = 1, 2, . . . T; Bt ðst1 ; qt ; st Þ is the hydropower output in period t, in which st1 and st are the beginning and ending storages, and qt is the inflow; Eqt is the expectation operator; and Dt is the decision interval. The following constraints should be included in the model:

2. Methodology

f t ðHt Þ ¼ maxfEqt jHt Bt ðst1 ; qt ; st Þ þ Eqt jHt  EHtþ1 jqt ;Ht  f tþ1 ðHtþ1 Þ

2.1. Development of recursive equation

where Eqt jHt is the conditional expectation operator for the inflow qt in period t for a specific state Ht ; and EHtþ1 jqt ;Ht is the conditional expectation operator for the state in the next stage, which is depended on the state and decision in period t. In this study, the inflow and the beginning storage of period t (qt and st1 ) are taken as the state variables, and the ending storage of period t (st ) is taken as the decision variable. Therefore, the recursive equation is defined as:

(1) Water balance constraint.

st ¼ st1 þ ðqt  rt ÞDt (2) Water storage constraint.

st;min 6 st 6 st;max

( ) T X f ðst1 ; qt Þ ¼ max Eqt ðBt ðst1 ; qt ; st ÞÞDt t¼1

ð1Þ

ð3Þ

(3) Flow constraint.

r t;min 6 r t 6 r t;max ; Q t 6 Q max

ð4Þ

(4) Hydropower output constraint.

Nt 6 Nmax

ð5Þ

(5) Relationship curves between reservoir water storage and its water level.

Z t ¼ f ðV t Þ and V t ¼ f

1

ðZ t Þ

ð6Þ

(6) Relationship curves between tail water level and its generation discharge.

Q t ¼ gðZ dr;t Þ and Z dr;t ¼ g 1 ðQ t Þ

ð7Þ

where rt rt;max ; rt;min st;max ; st;min Q t , Q max N t , N max Z t , Z dr;t f ðÞ gðÞ

The water release of the reservoir in period t, m3/ s The upper and lower limits of the reservoir release in period t, m3/s The upper and lower limits of the reservoir storage in period t, m3 The generation discharge and its upper limit in period t, m3/s The electric output in period t and the installed capacity, kW.h The ending and tail water levels of the reservoir in period t, m The function between reservoir water level and its water storage The function between tail water level and its generation discharge

2.1.2. Recursive equation considering the expected future value In the classical SDP, the optimal reservoir decision can be derived by solving the following recursive equation (TejadaGuibert et al., 2010): st

2.1.1. Objective function and constraints This study focuses on the reservoir whose main purpose is to maximize the total hydropower generation. Thus, the objective function is as follows:

ð2Þ

ð8Þ

267

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

( f t ðst1 ; qt Þ ¼ max Bt ðst1 ;j qt ; st Þ þ st

X ptþ1 ðj qtþ1 jqt Þf tþ1 ðst ; qtþ1 Þ

) gT +1 ( sT ,i )=0,i =1,2,...,N

Storage dispersion into N points

j

ð9Þ

Calculate recursively using Eq.(13)

where j is the inflow class index; ptþ1 ðj qtþ1 jqt Þ is the conditional probability of qtþ1 2 j based on qt . In the last period, the recursive equation is:

(

f T ðsT1 ; qT Þ ¼ max BT ðsT1 ;j qT ; sT Þ þ sT

) X pTþ1 ðj qTþ1 jqT Þf Tþ1 ðsT ;qTþ1 Þ

f1 ( s0, i ),i =1,2,...,N

fT +1 ( sT ,i )=f1 ( s0,i )

No

Convergence criterion met?

j

ð10Þ In previous studies, the value beyond period T (hereafter referred to as future value function) is often ignored, and thus defaulting f Tþ1 ðsT ; qTþ1 Þ ¼ 0. It means that there is no difference in future value no matter where the ending water storage of period T is. This is obviously conflicting with the actual situation. In fact, different carryover storages will bring different future benefits. In this study, the expected future value is simplified as a function of the storage:

X pTþ1 ðj qTþ1 jqT Þf Tþ1 ðsT ; qTþ1 Þg ¼ EqTþ1 ðf Tþ1 ðsT ; qTþ1 ÞÞ ¼ g Tþ1 ðsT Þ

Yes Calculate the minimal expected future value gT +1 min

gT +1 ( sT ,i )=gT +1 ( sT ,i )-gT +1,min i 1, 2,..., N

Output the function gT +1 ( sT )

Fit the relationship of sT ,i with gT +1 ( sT ,i ) i 1, 2,..., N

Fig. 1. The flowchart to calculate the expected future value function.

j

ð11Þ The long-term inflow is a typical periodic Markov process with a known period T. Accordingly, the reservoir operation is a controlled Markov process. Su and Deininger (1972) have proved that the recursive calculation results in a constant difference between the values of adjacent cycles at the same period and state, which is the optimal expected benefit of the cycle. Since the discrete state of the beginning storage in the first period is the same as that of the ending storage in the last period, the following equation can be obtained:

g 1 ðs0 Þ ¼ g Tþ1 ðsT Þ þ C; 8 i ¼ 1; 2; . . . ; N

ð12Þ

where i is the storage discrete index; N is the dispersion number, and C is a constant which is the optimal expected benefit of the cycle. Based on the characteristics of the periodic Markov process, the expected future value function can be calculated by the following steps (see Fig. 1): Step 1: Assuming that the decision (ending storage) in period t only depends on the beginning storage in period t, and develop the following recursive equation:

( X

g t ðst1 Þ ¼ max st

)

pt ðj qt ÞBt ðst1 ;j qt ; st Þ þ g tþ1 ðst Þ

ð13Þ

j

Step 2: Disperse the beginning and ending storages in each period into N points between the upper and the lower limits of the water storage; Step 3: Assuming that the expected future value in each storage state is 0, that is g Tþ1 ðsT;i Þ ¼ 0; i ¼ 1; 2; . . . ; N; Step 4: Calculate recursively using Eq. (13) to obtain the value in each beginning state of period 1, thus obtaining the g 1 ðs0;i Þ; i ¼ 1; 2; . . . ; N; Step 5: Given the convergence criterion:

g 1 ðs0;i Þ  g Tþ1 ðsT;i Þ ¼ C; 8 i ¼ 1; 2; . . . ; N

ð14Þ

If the convergence criterion is met, go to Step 6; otherwise, let g Tþ1 ðsT;i Þ ¼ g 1 ðs0;i Þ and return to Step 4; Step 6: Calculate the minimal expected future value g Tþ1;min by Eq. (15):

g Tþ1;min ¼ ming Tþ1 ðsT;i Þ; i ¼ 1; 2 . . . N i

ð15Þ

Step 7: Update g Tþ1 ðsT;i Þ ði ¼ 1; 2; . . . ; NÞ by minus the g Tþ1;min ; Step 8: Fit the relationship between sT;i and g Tþ1 ðsT;i Þ, and obtain the value function g Tþ1 ðsT Þ. Finally, the recursive equation of the last period considering the expected future value is:

  f T ðsT1 ; qT Þ ¼ max BT ðsT1 ;j qT ; sT Þ þ g Tþ1 ðsT Þ sT

ð16Þ

As the iteration number increases, the expected future value will accumulate accordingly. For SDP, it is the difference between expected future values in different ending storage states, rather than the values themselves, that can have an effect on the operating results. Therefore, the Step 6 and Step 7 are used to avoid too large expected future values without changing its role in stochastic optimal operation. 2.2. Transition probability matrix based on copula functions 2.2.1. Conditional probability expression Assuming that the marginal distribution of inflows qt and qtþ1 are Fðqt Þ and Fðqtþ1 Þ, respectively. According to Sklar’s theorem (Sklar, 1959), their joint distribution F qt ;qtþ1 ðqt ; qtþ1 Þ is:

F qt ;qtþ1 ðqt ; qtþ1 Þ ¼ Cðu; v Þ ¼ CðFðqt Þ; Fðqtþ1 ÞÞ

ð17Þ

u ¼ Fðqt Þ

ð18Þ

v ¼ Fðqtþ1 Þ

ð19Þ

Given x1 < qt 6 x2 , the conditional probability of qtþ1 6 y2 is:

Pðqtþ1 6 y2 jx1 < qt 6 x2 Þ ¼ ¼

F qt ;qtþ1 ðx2 ; y2 Þ  F qt ;qtþ1 ðx1 ; y2 Þ F qt ðx2 Þ  F qt ðx1 Þ

Cðu2 ; v 2 Þ  Cðu1 ; v 2 Þ u2  u1

ð20Þ

where CðÞ is a kind of copula function; u ¼ Fðqt Þ and v ¼ Fðqtþ1 Þ are the marginal distributions of qt and qtþ1 ; ui ¼ Fðxi Þ and v i ¼ Fðyi Þ, i = 1, 2. Given the inflow interval qt 2 ðx1 ; x2  in the period t, the condition probability of qtþ1 2 ðy1 ; y2  can be calculated by using Eq. (21).

Pðy1 < qtþ1 6 y2 jx1 < qt 6 x2 Þ ¼ Pðqtþ1 6 y2 jx1 < qt 6 x2 Þ  Pðqtþ1 6 y1 jx1 < qt 6 x2 Þ

ð21Þ

268

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

2.2.2. Joint distribution based on Archimedean copulas The copula function C is a bivariate Archimedean copula, if it satisfies the following condition:

Cðu; v Þ ¼ u1 ðuðuÞ þ uðv ÞÞ u; v 2 ½0; 1

ð22Þ

where uðÞ is the generator function of the copula, and u1 is its inverse function. The uðÞ is a continuously decreasing convex function, whose uð1Þ ¼ 0 and uð0Þ ¼ 1. In this study, three widely used Archimedean copula functions are applied, including, GumbelHougaard (hereafter referred to as Gumbel), Clayton and Frank copulas. The explicit expression of these copula functions and their generator functions are listed in Table 1. There is one dependence structure parameter in the Archimedean copulas. Nelsen (2006) has shown that, if X and Y are the random variables with a Archimedean copula C generated by uðtÞ, the population version of Kendall’s tau (s) for X and Y can be given by

s¼1þ4

Z

1 0

uðtÞ dt u0ðtÞ

ð23Þ

The relationship between the dependence structure parameter h and Kendall’s s of the three copulas can be derived from Eq. (23) (see Table 1). Let (X1, Y1) and (X2, Y2) be independent and identically distributed random vectors, each with joint distribution function C. Then, the Kendall’s s is defined as the probability of concordance minus that of discordance:

s ¼ sX;Y ¼ P½ðX 1  X 2 ÞðY 1  Y 2 Þ > 0  P½ðX 1  X 2 ÞðY 1  Y 2 Þ < 0 ð24Þ The Kendall’s s can be calculated using observation samples, and then the copula parameter can be estimated by using its functional relationship with Kendall’s s. 2.2.3. Selection of an appropriate copula Three indexes are used to evaluate the appropriateness of the hypothesized copulas, including Cramer-von Mises statistic (Sn), P-value (P) and Akaike information criterion (AIC). Cramer-von Mises statistic can be calculated by the following equation:

Sn ¼

n X

fC n ðU i;n ; V i;n Þ  C hn ðU i;n ; V i;n Þg2

ð25Þ

i¼1

where n is the number of observation sample; ðU i;n ; V i;n Þ is the marginal distribution function value of the observation sample ðX i;n ; Y i;n Þ; C n ðU i;n ; V i;n Þ is the empirical copula value of ðU i;n ; V i;n Þ; hn is the dependence structure parameter obtained from the n observation samples; and C hn ðU i;n ; V i;n Þ is the theoretical copula value of ðU i;n ; V i;n Þ. P-value is an attached index of Sn introduced by Genest et al. (2009), and it can be calculated by parametric bootstrapping method as follows:

Table 1 Copula function Chðu; v Þ, generating function uðtÞ and functional relationship of Kendall’s sh with copula parameter for various Archimedean copulas. Copula family Goumbel Clayton Frank

Chðu; v Þ

uðtÞ

 h i1  h h h exp  ð ln uÞ þ ð ln v Þ  1 max uh þ v h  1 h ; 0  hu hv 1Þ  1h ln 1 þ ðe 1Þðe eh 1

Dk ðxÞ is the Debye function, Dk ðxÞ ¼ xkk

Rx

tk 0 et 1 dt

Kendall’s

ð ln tÞ

ðh  1Þ=h

ðt h  1Þ=h

h=ðh þ 2Þ

h

eht 1

 ln eh 1

sh

1 þ 4h ½D1 ðhÞ  1

for any positive integer k.

(6) The P-value can be obtained byp ¼ N1 lðÞ is a logical indicator.

PN

k k¼1 lðSn

P Sn Þ; where

More information about Sn and P-value can be found in Genest et al. (2009) and the references therein. AIC can be calculated by the following expression:

AIC ¼ n ln

  Sn þ2k n

ð27Þ

The smaller the AIC and Sn values are, and the higher the P-value is, the more strongly the assumed copula is accepted. 3. Case study 3.1. Ertan reservoir Ertan reservoir is a key project for Sichuan power grid in China, which is located in the lower reaches of the Yalong River basin in Sichuan province of southern China (Fig. 2). The main purpose of Ertan reservoir is to maximize the total power generation by rational utilization of the water resources. In order to ensure the safety of the middle and lower reaches of Yangtze River in flood season, it is required to reserve 16 billion m3 flood control capacity in June and July. The water level in these two months should be maintained lower than 1192 m. The main characteristic parameters of Ertan reservoir are listed in Table 2. The 55-year observed inflows from 1958 to 2012, collected on ten-days1 and monthly basis is used to calculate state-transition matrix and simulate the operation of Ertan reservoir. During the low-flow period from November to March, the reservoir is mainly recharged by groundwater with low uncertainty; while during the high-flow period from April to October, the reservoir is mainly recharged by rainfall with relatively higher uncertainty. Thus, in this study, each year is divided into 26 periods, with each period corresponding to a month in the low-flow period from November to March, whereas to approximately 10 days in the high-flow period from April to October. 3.2. Results and discussions

(1) Give the value of N and initialize k = 1; (2) Generate random samples ðU k1;n ; V k1;n Þ; . . . ; ðU kn;n ; V kn;n Þ from the copula function C hn ; (3) Derive the dependence structure parameter hkn and empirical value C kn ðU ki;n ; V ki;n Þ k k ðU 1;n ; V 1;n Þ; . . . ; ðU kn;n ; V kn;n Þ; copula

from

pseudo-observations

(4) Calculate the Cramer-von Mises statistic Skn :

Skn ¼

n X 2 fC kn ðU ki;n ; V ki;n Þ  C hkn ðU ki;n ; V ki;n Þg

3.2.1. Marginal distribution In China, Pearson Type III (PTIII) distribution is commonly used to characterize the distribution of hydrologic variables (Yun and Singh, 2008). In this study, the PTIII distribution is assumed to fit the marginal distribution, and then the Kolmogorov-Smirnov (KS) test (Massey, 1951) is used to evaluate the goodness-of-fit. The K-S statistic is defined as:

ð26Þ

i¼1

(5) If k < N, k ¼ k þ 1; return to Step (2); otherwise go to Step (6);

1 In China, a month is divided into three periods: the first 10-day (early), the middle 10-day (mid) and the rest day (late) of the month. Taking May as an example, the first 10-day of May is denoted by early May, the middle 10-day of May is mid May, and the last 10-day of May is late May.

269

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

Yalong River Basin N

Guandi

S

Ertan Tongzilin

Jinping-II Jinping-I

Legend Reservoir River

Fig. 2. The geographic location of Ertan reservoir.

Table 2 Basic parameters of Ertan reservoir. Parameters

Value

Parameters

Value

Installed capacity (MW) Firm output (MW) Maximum power flow (m3/s) Minimum release (m3/s)

3300 1028 2202

Normal pool level (m) Dead water level (m) Dead storage capacity (106 m3) Total storage capacity (106 m3)

1200 1155 2400

Dmax ¼ maxðjFðxÞ  GðxÞjÞ

20

5793

ð28Þ

where FðxÞ and GðxÞ are the theoretical and empirical cumulative distribution functions (CDFs), respectively. Given the significance level a, the critical value of KS test (Da ) can be obtained. The assumption is valid only when Dmax 6 Da . The smaller the Dmax is, the better the fit will be. Table 3 shows the parameters and goodness-of-fit statistics of the marginal distributions. The parameters of PTIII distribution, including  x, C v and C s , are optimized by genetic algorithm (GA), whose objective is to minimize the sum of squares of deviations between the theoretical and empirical frequencies. Given

a ¼ 0:05, the critical value D0:05 ¼ 0:1831. Table 3 shows that Dmax < D0:05 for all periods, indicating that the PTIII distributions are satisfactory. 3.2.2. Joint distribution The joint distribution of inflows in adjacent periods is established by the above three copula functions (Gumbel, Clayton and Frank copulas). The dependence structure parameter h is estimated based on its functional relationship with Kendall’s tau (see Table 1). Based on the three selection criteria, the final copula function for each period is selected (see Table 4). Goumbel copula can be selected for 11 periods, Clayton copula for 3 periods, and Frank copula for 12 periods, respectively. However, it is noted that Frank copula can also achieve satisfactory performance for the 3 periods selected Clayton copula. Thus, the Clayton copula seems to be not as effective as the other two copulas in this study; while the performance of Goumbel copula is as good as that of Clayton copula, both of which have relatively low Sn and AIC values and

high P-value. Thus, we can use either Goumbel or Frank copula for each period. In this study, we use the selected copula for each period. 3.2.3. Transition probability matrix The inflow in each period is divided into 11 classes according to the CDF thresholds: 5%, 10%, 20%, 30%, 40%, 50%, 60%, 70%, 80% and 90%. And the following three schemes are used to calculate the transition probability matrix: Scheme I: Conventional method based on measured streamflow series (hereafter referred to as Con-MI). Scheme II: Conventional method based on synthetic streamflow series (hereafter referred to as Con-MII), in which the synthetic streamflow series are generated, and then the conventional method is used. Scheme III: Copula function method based on measured streamflow series (hereafter referred to as Copula-M), in which the joint distribution of adjacent periods is established, and then the conditional probability formula is derived to calculate the transition probability matrix. In Con-MII, to characterize the nonstationary properties and the periodic dependence structure of the seasonal hydrologic time series, the Periodic Auto Regression Moving Average (PARMA) models are used to generate the synthetic streamflow series. The PARMA (p, q) model can be expressed as (Salas, 1993)

Y v ;s ¼

p q X X /i;s Y v ;si þ ev ;s  aj;s ev ;sj i¼1

ð29Þ

j¼1

where Y v ;s represents the streamflow process for year v and season

s, and it is normally distributed with mean zero and variance r2s ðYÞ for each season s; the em;s is the uncorrelated noise term which is also normally distributed with mean zero and variance r2s ðeÞ for each season s; f/1;s ; . . . ; /p;s g are the periodic autoregressive parameters; and fa1;s ; . . . ; aq;s g are the periodic moving average parameters. If the number of seasons or periods is x, then a PARMA(p,q) model consists of x numbers of individual ARMA(p,

q) models, where the dependence is across seasons instead of years. In this study, the PARMA (1, 1) model is used to simulate the inflow of Ertan reservoir for 5000 years and the Least Squares

270

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

Table 3 The parameters and goodness-of-fit statistics of marginal distributions. Parameter*

Period

 x Cv Cs Dmax

1

2

3

4

5

6

7

8

9

10

11

12

13

511 0.11 0.27 0.081

448 0.1 0.24 0.049

437 0.09 0.19 0.055

475 0.12 0.24 0.060

518 0.17 0.34 0.105

608 0.19 0.39 0.066

699 0.21 0.41 0.091

807 0.23 0.45 0.069

1001 0.28 0.56 0.066

1365 0.35 0.7 0.058

1946 0.34 0.68 0.080

2812 0.34 0.67 0.057

3517 0.33 0.67 0.082

14

15

16

17

18

19

20

21

22

23

24

25

26

3763 0.34 0.68 0.075

3392 0.3 0.6 0.050

3397 0.38 0.76 0.078

3404 0.45 0.9 0.095

3578 0.38 0.76 0.065

3713 0.41 0.82 0.090

3585 0.34 0.69 0.068

3237 0.29 0.57 0.066

2772 0.3 0.59 0.079

2168 0.26 0.53 0.074

1677 0.22 0.43 0.064

1105 0.16 0.32 0.050

690 0.14 0.27 0.063

Parameter

Period

 x Cv Cs Dmax *

 x, C v and C s are the mean, variation coefficient and skewness coefficient.

Table 4 The parameters and goodness-of-fit statistics of joint distributions. Goumbel

Selected copula

P-value

AIC

h

Sn

P-value

AIC

h

Sn

P-value

AIC

6.2 3.5 1.6 1.8 1.8 1.6 1.5 2.0 1.7 2.2 1.8 1.4 1.7 1.6 1.7 2.0 1.6 2.1 1.6 2.1 2.4 2.7 3.4 4.8 5.2 6.4

0.016 0.022 0.046 0.025 0.024 0.030 0.037 0.024 0.026 0.026 0.030 0.022 0.045 0.020 0.021 0.024 0.021 0.053 0.079 0.018 0.043 0.024 0.020 0.024 0.019 0.029

1.000 0.998 0.924 1.000 0.996 0.994 0.980 0.998 0.994 0.996 0.996 0.998 0.958 0.998 0.996 0.998 0.998 0.896 0.710 1.000 0.928 0.998 1.000 0.998 0.998 0.982

446 428 388 422 425 410 399 424 419 419 411 427 389 434 430 424 432 380 358 439 391 424 434 424 435 413

5.6 3.3 1.1 1.2 1.3 1.4 1.3 1.5 1.4 1.7 1.8 0.6 0.9 0.7 1.2 1.2 1.0 2.1 1.5 1.5 1.8 2.0 2.6 6.3 6.8 5.7

0.093 0.077 0.056 0.049 0.042 0.020 0.016 0.064 0.024 0.078 0.037 0.041 0.101 0.051 0.037 0.075 0.050 0.079 0.035 0.074 0.126 0.126 0.119 0.042 0.044 0.108

0.682 0.738 0.908 0.940 0.966 1.000 1.000 0.850 0.998 0.794 0.984 0.964 0.602 0.926 0.980 0.784 0.908 0.720 0.984 0.752 0.518 0.476 0.504 0.960 0.930 0.564

349 359 377 384 392 433 445 370 423 359 399 394 345 382 400 361 383 358 403 362 332 332 336 393 390 341

22.4 12.1 4.5 4.6 5.0 4.1 4.2 6.0 4.5 7.4 5.2 2.3 4.8 3.4 4.3 4.6 4.3 7.7 4.9 6.3 9.0 9.0 12.5 20.1 22.3 24.9

0.024 0.026 0.028 0.031 0.020 0.022 0.023 0.031 0.019 0.020 0.031 0.030 0.028 0.035 0.016 0.047 0.023 0.025 0.036 0.018 0.023 0.032 0.017 0.016 0.013 0.029

0.992 1.000 0.998 0.984 1.000 1.000 1.000 0.984 1.000 0.998 0.995 0.998 0.994 0.982 0.998 0.928 0.998 0.998 0.980 1.000 0.996 0.988 1.000 1.000 0.998 0.988

424 420 416 409 435 427 426 410 435 435 410 410 414 403 445 387 426 420 401 438 427 408 442 445 457 414

20

Average difference (m 3/s)

Frank

Sn

10

0

-10

-20

5

10

15

20

25

80

60

40

20

0

Goumbel Goumbel Frank Goumbel Frank Clayton Clayton Goumbel Frank Frank Goumbel Goumbel Frank Goumbel Frank Goumbel Goumbel Frank Clayton Goumbel Frank Goumbel Frank Frank Frank Frank

0.03

Variation coefficient difference

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

Clayton

h

Standard deviation difference (m 3/s)

Period

5

Period

10 15 Period

(a)

(b)

20

25

0.025 0.02 0.015 0.01 0.005 0

5

10

15 Period

20

25

(c)

Fig. 3. The main statistic parameters of the observed and synthetic streamflow series: (a) average difference; (b) standard deviation difference; and (c) variation coefficient difference.

(LS) method is used to estimate its parameters. Fig. 3 shows the main parameters of the observed and synthetic streamflow series. The inflow in dry season is relatively stable, and thus the synthetic

streamflow series can hold the statistical properties of the original series; while the inflow in flood season has higher uncertainty, resulting in a larger difference in the average, standard deviation

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

and variation coefficient between the two series, especially from periods 6 to 23. On the whole, the synthetic streamflow series can well characterize the statistical properties of the original series. Fig. 4 shows the transition probability matrices obtained in dry (January and November) and wet (late June and mid September) periods by the three schemes, respectively. Because the inflows in dry periods are mainly recharged by groundwater with lower uncertainty, there is an significant autocorrelation of inflows between adjacent periods, and thus the transition probabilities are centered on the one of the diagonal; while the inflows in flood periods are mainly recharged by rainfall with relatively higher uncertainty, the transition probabilities are distributed in a wider region. For Con-MII and Copula-M schemes, the transition probabilities are centered on the one of the diagonal, and decrease along with the other diagonal. Due to the autocorrelation of inflows, the transition probabilities for inflows of adjacent periods to fall into the same class are higher, and the larger the inflow class differences between two adjacent periods, the smaller the transition probabilities will be. The advantage of Con-MI is that it does not relay on a specific parametric form of the joint distribution or a simulation model to generate the synthetic streamflow series. However, as its transition probability is calculated based on the assumption that the frequency is equal to the probability, the representation errors cannot be avoided due to the limited number of observed samples. On the one hand, many values in its transition probability matrices are zeros, some of which appear to be unreasonable. For example, when the inflow of middle September is in class 8, the inflow of late September just has two situations, with a probability of 0.6 falling into class 7 and a probability of 0.4 falling into class 11. However, due to the stochastic and fuzzy nature of hydrological phenomena, the inflow of late September is also likely to fall into other classes, especially class 8–10. Maybe the inflows in these classes happened in the past, but due to the limited age of the hydrological stations, they were not recorded. On the other hand, the transition probability cannot be calculated when there is no observed inflow data in some classes (Fig. 4(d), whose transition probability is represented by 0.5). Let Pðqtþ1 jqt 2 iÞ represent the conditional probability of qtþ1 when qt is in class i. If there is no data meeting qt 2 i, Pðqtþ1 jqt 2 iÞ cannot be calculated by the Con-MI scheme. In this study, Pðq7 jq6 2 1Þ, Pðq15 jq14 2 1Þ, Pðq19 jq18 2 1Þ, Pðq21 jq20 2 6Þ and Pðq23 jq22 2 8Þ cannot be calculated. By generating the synthetic streamflow series or building the marginal distributions to interpolate and extrapolate the inflows, Con-MII and Copula-M can avoid the representation errors of Con-MI. The behavior of Copula-M and Con-MII looks similar. Both of them can avoid the representation errors of Con-MI. On the one hand, their transition probabilities are distributed more uniformly with less zero values. Meanwhile, the situation that the transition probabilities are unable to be calculated does not happen. However, it should be noted that Copula-M has three advantages over Con-MII. (1) Calculation is much simpler than the conventional methods. As it only needs to calculate the structure parameters of the joint distribution between adjacent periods, then the explicit expression for calculating the transition probability can be obtained. The structure parameter h of the Archimedean copulas can be estimated based on Kendall’s s using some simple equations (see Table 1). Then, the transition probability matrices can be obtained easily using Eqs. (20) and (21). However, for Con-MII, after generating the synthetic streamflow series, the transition probability is still calculated using the conventional method.

271

(2) Generation of synthetic streamflow series can be avoided. The use of PTIII distribution as the marginal distribution allows for the implicit interpolation and extrapolation of the inflows without depending on other synthetic streamflow generating models. This is consistent with the inflow class division, which also relies on the PTIII distribution. Thus, the uncertainty resulting from the simulation model can be reduced. In addition, in the BSDP model, we should consider uncertainties resulting from both natural streamflow and forecast. It is very difficult to generate the synthetic natural and forecasted inflow series at the same time without changing the characteristics of the original ones due to the complicated correlation relationship between them. The Copula-M can avoid the generation of synthetic streamflow series, and thus it can be directly used in BSDP. (3) Data storage can be reduced obviously. For Con-MII, the transition probability matrices in each period need to be stored, and thus the data storage grows exponentially with the number of the inflow class. However, for Copula-M, only the parameters of marginal and joint distributions ( x, C v ; C s and h) in each period need to be stored. 3.2.4. Future value function The beginning and ending storages in each period are dispersed into 20 points, and then a sample set of ðsTþ1;i ; g Tþ1 ðsT;i ÞÞ; i ¼ 1; 2 . . . ; 20 is obtained by the stepwise iteration method introduced in Section 2.1.2. Through a simple linear fitting, the final future value function is:

gðsTþ1 Þ ¼ 0:0041  sTþ1  9:7383

ð30Þ

With the increase of the carryover storage, the future value increases gradually (see Fig. 5). It means that the more water is stored in the reservoir, the more hydropower generation may win in the future due to the higher water head or more release. As is shown in Fig. 5, when the reservoir is at dead water level with a carryover storage of 2400 million m3, the expected future value is 0. While when the reservoir is at normal water level with a carryover storage of 5793 million m3, it can generate more 14:1734  108 kW.h hydropower. However, the more water stores in the reservoir, the less water is released that may result in a decrease of hydropower generation in the current year. The problem of how much water to withhold from immediately beneficial deliveries, retaining that water in storage, is known as ‘‘ hedging’’. For water-supply operating problem, the operating benefit has a relatively simple relationship with water volume. Some analytical hedging rules considering the balance between beneficial release and carried-over storage value have been derived (Draper and Lund, 2004; Shiau, 2011; You and Cai, 2008a,b; Zhao et al., 2010). However, due to the coupling relation of water head and release, the future value function is much complicated for hydropower generation operation. In this study, the expected future value function is obtained by a simple stepwise iteration method, and used in SDP to balance the benefit between the current year and the future. This function not only can be used in SDP, but also in Deterministic Dynamic Programming algorithm (DDP). For example, in making an annual hydropower generation schedule, the year-end water level can be controlled dynamically using Eq. (30) according to the inflow of the year, which is especially important for reservoirs with multi-year regulation performance. 3.2.5. Operating policy Based on the transition probability matrices obtained by the three schemes, the operating policy can be obtained by backward recursion. In this study, the SDP based on the three schemes are

272

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

Con-MI

Con-MII

1

Probability

Probability

0.8 0.6 0.4 0.2

Copula-M

0.8

0.8

0.6

0.6

Probability

(a)

0.4 0.2

0.4 0.2

0

0

0 1

1

2

3

Inf l ow

4

5

6 cl a 7 ss 8 of 9 pe ri o 10 d2 11 (F eb )

1

4

6

5

1 riod f pe ss o w cla o fl In

2

3

9

8

7

1 3

4 Inf 5 l ow 6 cl a 7 ss 8 of pe 9 ri o d 2 1011 (F eb )

11

10

2

) (Jan

1

3

4

5

f ss o w cla Infl o

2

6

7

d perio

10

9

8

2

Inf 3 l ow 4 cl a 5 ss 6 7 of pe ri o 8 9 d2 (F 1011 eb )

11

an ) 1 (J

4

6

5

riod f pe ss o w cla Infl o

2

1

3

10

9

8

7

11

an ) 1 (J

(b)

Probability

Probability

0.6 0.4

0.8

0.8

0.6

0.6 Probability

1 0.8

0.4 0.2

0.2

0.4 0.2

0

0 1

0 1

2

3 Inf l ow 4 cl a 5 6 ss 7 of pe ri o 8 9 d2 6 ( 10 De 11 c)

2

3

4

f ss o w cla Infl o

1

6

5

d perio

9

8

7 25 (

10

1

2

Inf 3 l ow 4 cl a 5 ss 6 of pe 7 8 ri o 9 d2 6 ( 10 De 11 c)

11

) Nov

6

1

9

8

7

10

11

2

3

4 Inf l ow 5 6 cl a 7 ss of 8 pe 9 ri o d 2 10 11 6( De c)

v) 5 ( No 4 d 25 3 perio 2 s of s la wc Infl o

5

1

3

2

w Infl o

4 d perio s of clas

25

10

9

8

7

6

11

v) ( No

(c)

Probability

Probability

0.6 0.4

0.4

0.5

0.3

0.4

Probability

1 0.8

0.2 0.1

0.3 0.2 0.1

0.2 0

0 1 Inf 2 l ow 3 cl a 4 ss 5 of 6 pe 7 ri o d1 8 9 3( ea 10 rly 11 Ju l)

1

11

9

Con-M(I)

Inf 1 2 l ow 3 cl a ss 4 5 of pe 6 7 ri o d1 8 3( 9 ea rly 10 11 Ju l)

11 10 un) J e t 6 (la 5 d 12 4 perio 3 2 s of 1 clas w o Infl 7

8

9

Con-M(II) 0.6

0.5

0.4

0.4

0.2

1

ow

3

4

c la ss 5 6 of eri 7 od 8 21 (la 9 10 te S 11 ep )

6

7

d 4 er io 3 ss of p 2 la 1 flow c In 5

8

20 (

9

m

10

p id Se

)

11

Infl

11

) Jun

0 1

2

lat e 5 12 ( 4 riod 3 f pe 2 ss o la 1 c w Infl o

0.2

0

-0.5

Infl

Probability

0.6

0

10

9

8

7

6

Copula-M

1

Probability

Probability

(d)

10 ) Jun 7 lat e 6 5 12 ( d 4 erio p f 3 ss o 2 w cla Infl o 8

0

Inf 1 l ow 2 cl a 3 ss 4 of 5 pe 6 ri o 7 d1 8 3( 9 ea rly 10 Ju 11 l)

ow

2

3

c la ss o

Infl 4

5

f er

6

io d

7

21

8

( la

9

10

te S 11 ep )

1

2

Inflo

3

4

5

ss o w c la

6

7

8

20 r iod f pe

9

(m

10

11

p id Se

)

1

ow

2

3

c la ss 4 5 of er i 6 od 7 21 8 ( la 9 te S 10 e p 11 )

1

2

3

Inflo

4

5

ss w c la

6

7

8

9

d 20 er io of p

10

11

Se (m id

p)

Fig. 4. The transition probability matrices obtained by the three schemes in (a) January; (b) November; (c) late June; and (d) mid September. The three columns show the results obtained by Con-MI, Con-MII and Copula-M respectively.

referred to as SDP (Con-MI), SDP (Con-MII) and SDP (Copula-M), respectively. For the transition probabilities that are unable to be calculated using Con-MI, we assume that they occur in each class with the same probability just as the traditional treatment. Just for example, Fig. 6 shows the operating policy in January (dry season) and late June (flood season) obtained by SDP (Copula-M).

When the beginning storage of the reservoir and the class of the inflow are known, the ending storage can be obtained using the operating policy. It is found that the operating policy is more sensitive to inflow class division in flood season than that in dry season. In this study, we set the same number of inflow class division according to the frequency thresholds. In fact, the inflow can be

273

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

Future value (108 kW.h)

15

(5793,14.1734) 10

y = 0.0041x - 9.7383

5 (2400,0) 0 2000

2500

3000

3500

4000

4500

5000

5500

6000

Carryover storage (106 m3) Fig. 5. The expected future value function of the Ertan reservoir.

dispersed with a larger step in dry season to avoid the curse of dimensionality, but with a smaller step in flood season to provide a more detailed operating policy and reduce the effect of inflow uncertainty. 3.2.6. Simulation operation results The operating policies are used to guide the operation of inflows from 1959 to 2012, and the performance of the hydropower system is evaluated using two indicators, including average annual hydropower generation (AAHG) and average annual spillage (AAS) (Table 5). Con-MII can overcome the representation errors of Con-MI by generating the synthetic streamflow series. As a results, the AAHG of the whole year can be improved from 168:67  108 kW:h to 168:94  108 kW:h, with an increase of 0:27  108 kW:h. Copula-M can solve the representation errors of Con-MI without depending on some other simulation models by building a continuous marginal distribution of inflows. The hydropower generation of SDP (Copula-M) is increased by 0:34  108 and 0:07  108 kW:h compared with that of SDP (Con-MI) and SDP (Con-MII), respectively, and this increase is observed mainly in flood season. Fig. 3 shows that the synthetic streamflow series can hardly hold the characteristics of the original series in flood season, and thus the transition probability matrices obtained by Con-MII in flood season have relatively larger errors and the corresponding power generation of SDP is lower than that of Copula-M. Fig. 7 shows the average storage process of the Ertan reservoir guided by the three operating policies. As the transition probability matrices of Copula-M and Con-MII are quite similar, the storage

5000

3

5500

i=1 i=2

4500

6

6

3

Reservoir ending storage (10 m )

6000

Reservoir ending storage (10 m )

processes of SDP (Copula-M) and SDP (Con-MII) are quite close. In this study, Pðq7 jq6 ¼ 1Þ, Pðq15 jq14 ¼ 1Þ, Pðq19 jq18 ¼ 1Þ, Pðq21 jq20 2 6Þ and Pðq23 jq22 ¼ 8Þ are assumed to occur in any inflow class with the same probability. However, this assumption may have an effect on the final operating decision. The storage of SDP (Con-MI) is lower than that of SDP (Copula-M) and SDP (ConMII), especially from periods 15 to 23, which can explain why the AAS of the SDP (Con-MI) is lower than that of the other schemes (see Table 5). To further analyze the difference of operating results obtained by SDP (Copula-M) and SDP (Con-MI), Fig. 8 shows the ratio of the storage and output obtained by SDP (Copula-M) to that obtained by SDP (Con-MI). As we can see, the storage of the SDP (Copula-M) is always higher than that of SDP (Con-MI) in flood season. A higher storage results in a higher hydropower water head, and thus the output of SDP (Copula-M) is higher in the most periods of the flood season (except in periods 7 and 13) and the AAS is relatively higher (see Table 5). In period 7, the reservoir storage declines (see Fig. 7). The storage of SDP (Con-MI) falls sharply, resulting in a higher release and output in this period. However, in the following periods 8 to 12, as the storage of SDP (CopulaM) is higher, its output is higher than that of SDP (Con-MI). In period 13, the reservoir is impounding (see Fig. 7). The SDP (CopulaM) can guide the reservoir to store water as soon as possible, and thus the release and output are lower than that of SDP (Con-MI) in this period. From periods 15 to 23, the reservoir storage of SDP (Copula-M) is obviously higher than that of SDP (Con-MI). In these periods, the reservoir storage is high, and the output is mainly controlled by the reservoir hydropower water head. Therefore, the output of SDP (Copula-M) is always higher than that of SDP (Con-MI). In a short, the SDP (Copula-M) can reduce the release speed in the decline stage while speed up to store water in the impounding stage, making it possible for the reservoir to use the high hydropower water head to generate more hydropower. Similarly, Fig. 9 shows the ratio of the ending storage and output obtained by SDP (Copula-M) to that obtained by SDP (Con-MII). The storage of SDP (Copula-M) is always higher than that of SDP (Con-MII) in flood season, resulting in a relative higher AAS (see Table 5). Meanwhile, the output of SDP (Copula-M) is higher than that of SDP (Con-MII) due to the higher hydropower water head except in periods 7, 10, 17 and 19. In these four periods, the release is decreased to store more water for future use, thus the ending storage ratio increases while the hydropower generation decreases.

5000 4500 i=1,2,...7

4000

i=8,9 i=10

3500 3000 2500

i=11

3000

3500

4000

4500

5000 6

5500 3

Reservoie begining storage (10 m )

(a)

6000

i=3 4000

i=4 i=5 i=6 i=7

3500

i=8

3000

i=9 i=10 i=11

2500

2000

2500

3000

3500

4000

4500 6

3

Reservoie begining storage (10 m )

(b)

Fig. 6. The operating policy obtained by SDP (Copula-M) in (a) January; and (b) late June.

5000

274

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

Table 5 The application results of the three operating policies. Average annual hydropower generation (108 kW.h)

Operating policy

SDP (Con-MI) SDP (Con-MII)) SDP (Copula-M)

The whole year

Dry season

Flood season

168.67 168.94 169.01

36.55 36.55 36.55

132.13 132.40 132.46

SDP (Con-MI)

SDP (Copula-M)

94.70 96.27 96.40

SDP (Con-MII)

6000

6

3

Ending storage (10 m )

Average annual spillage (108 m3)

5000

4000

3000 5

10

15

20

25

20

25

Period Fig. 7. Average storage process of the Ertan reservoir.

1.1

The storage ratio of SDP (Copula-M) to that of SDP (Con- MI)

Ratio

The Output ratio of SDP (Copula-M) to that of SDP (Con- MI) Ratio=1

1.05

1

5

10

15 Period

Fig. 8. The ratio of the storage and output of SDP (Copula-M) to that of SDP (Con-MI).

1.02

Ratio

1

0.98 The Output ratio of SDP (Copula-M) to that of SDP (Con- MII) The storage ratio of SDP (Copula-M) to that of SDP (Con- MII)

0.96

Ratio=1

5

10

15

20

25

Period Fig. 9. The ratio of the storage and output of SDP (Copula-M) to that of SDP (Con-MII).

4. Conclusions In this study, we proposed a stochastic optimization operation model for hydropower generation reservoirs, in which the Copula function was used to calculate the transition probability matrices and a stepwise iterative method was used to calculate the expected

future value function. The transition probability matrices obtained by Copula function method were compared with that obtained by the conventional method based on observed and synthetic streamflow series, respectively. Then, the transition probability matrices were incorporated into SDP to obtain the operating polices. The following are the conclusions drawn from the present study.

X.-h. Lei et al. / Journal of Hydrology 557 (2018) 265–275

(1) Both the Copula-M and Con-MII can avoid the representative error of the Con-MI in calculating the transition probability. Meanwhile, in some special case, where the transition probability is unable to be calculated due to the lack of measured samples, the assumption that the inflow occurs in any class with the same probability will reduce the operating benefit; (2) Compared with Con-MII, there are three advantages of using Copula-M to calculate the transition probability. The calculation process is much simpler, the generation of synthetic streamflow series can be avoided and the data storage can be reduced obviously. (3) SDP (Copula-M) can improve the hydropower generation of SDP (Con-MII) and SDP (Con-MI), especially in flood season. This study is based on the assumption of perfect prediction of the runoff in the current period, and thus the actual operating performance is closely related to the prediction ability. Considering the forecast error and the non-stationarity of the future inflow, the obtained operating policy may not be as effective as expected. Much work has been done about considering the forecast uncertainty into SDP (Karamouz and Vasiliadis, 1992; Mujumdar and Nirmala, 2007; Stedinger et al., 1984; Xu et al., 2014). In this study, we choose the simplest recursive equation and focus on a single reservoir to avoid too much complicated situations. But it has proven that the Copula-M is effective to calculate the transition probability matric and can lead to a hydropower increase, so it can be directly used in other more complicated situations. Acknowledgement This paper was jointly supported by National Key Research and Development Project (2016YFC0402208, 2016YFC0401903, 2016YFC0400903, 2016YFC0402204 and 2016YFC0402201), National Key Technology R&D Program (2015BAB07B03), and State Key Laboratory Simulation and Regulation Water Cycle in River Basin (2016CG05). References Celeste, A.B., Billib, M., 2009. Evaluation of stochastic reservoir operation optimization models. Adv. Water Resour. 32 (9), 1429–1443. Cervellera, C., Chen, V.C.P., Wen, A., 2006. Optimization of a large-scale water reservoir network by stochastic dynamic programming with efficient state space discretization. Eur. J. Oper. Res. 171 (3), 1139–1151. Chandramouli, V., Deka, P., 2005. Neural network based decision support model for optimal reservoir operation. Water Resour. Manage. 19 (4), 447–464. Chandramouli, V., Raman, H., 2001. Multi-reservoir modeling with dynamic programming and neural networks. J. Water Resour. Plann. Manage. 127 (2), 89–98. Chou, K., Jiingyun You, G., Jang, J.H., 2016. Using ensemble streamflow prediction in the reservoir operation during drought by implicit and explicit stochastic optimization: case study in Shihmen Reservoir, EGU General Assembly Conference. Draper, A.J., Lund, J.R., 2004. Optimal Hedging and Carryover Storage Value. J. Water Resour. Plann. Manage. 130 (1), 83–87. Drouin, N., Gautier, A., Lamond, B.F., Lang, P., 1996. Piecewise affine approximations for the control of a one-reservoir hydroelectric system. Eur. J. Oper. Res. 89 (1), 53–69. Favre, A.C., Adlouni, S.E., Perreault, L., Thiémonge, N., Bobée, B., 2004. Multivariate hydrological frequency analysis using copulas. Water Resour. Res. 40 (1), 290– 294. Gal, S., 1979. Optimal management of a multireservoir water supply system. Water Resour. Res. 15 (15), 737–749. Genest, C., Rémillard, B., Beaudoin, D., 2009. Goodness-of-fit tests for copulas: a review and a power study. Insurance Mathematics Econ. 44 (2), 199–213. Harboe, R., 1993. Explicit Stochastic Optimization. Springer, Netherlands, pp. 295– 306. Jenq-Tzong, S., Wang, H.Y., Chang-Tai, T., 2010. Bivariate frequency analysis of floods using copulas. Jawra J. Am. Water Resour. Assoc. 42 (6), 1549–1564.

275

Karamouz, M., Vasiliadis, H.V., 1992. Bayesian stochastic optimization of reservoir operation using uncertain forecasts. Water Resour. Res. 28 (5), 1221–1232. Karmakar, S., Simonovic, S.P., 2009. Bivariate flood frequency analysis. Part 2: a copula-based approach with mixed marginal distributions. J. Flood Risk Manage. 2 (1), 32–44. Kim, Y.O., Palmer, R.N., 1997. Value of seasonal flow forecasts in Bayesian stochastic programming. J. Water Resour. Plann. Manage. 123 (6), 327–335. Koutsoyiannis, D., Economou, A., 2003. Evaluation of the parameterizationsimulation-optimization approach for the control of reservoir systems. Water Resour. Res. 39 (6), 250–252. Koutsoyiannis, D., Efstratiadis, A., Karavokiros, G., 2010. A decision support tool for the management of multi-reservoir systems. Jawra J. Am. Water Resour. Assoc. 38 (4), 945–958. Lee, J.H., Lee, J.M., 2006. Approximate dynamic programming based approach to process control and scheduling. Comput. Chem. Eng. 30 (10–12), 1603–1618. Loucks, D.P., Stedinger, J.R., Haith, D.A., 1981. Water Resource System Planning and Analysis. Massey, F.J., 1951. The Kolmogorov-Smirnov test for goodness of fit. J. Am. Stat. Assoc. 46 (253), 68–78. Mehta, R., Jain, S.K., 2009. Optimal operation of a multi-purpose reservoir using neuro-fuzzy technique. Water Resour. Manage. 23 (3), 509–529. Momtahen, S., Dariane, A.B., 2007. Direct search approaches using genetic algorithms for optimization of water reservoir operating policies. J. Water Resour. Plann. Manage. 133 (3), 202–209. Mousavi, S.J., Ponnambalam, K., Karray, F., 2005. Reservoir operation using a dynamic programming fuzzy rule-based approach. Water Resour. Manage. 19 (5), 655–672. Mujumdar, P.P., Nirmala, B., 2007. A Bayesian stochastic optimization model for a multi-reservoir hydropower system. Water Resour. Manage. 21 (9), 1465–1485. Nalbantis, I., Koutsoyiannis, D., 1997. A parametric rule for planning and management of multiple-reservoir systems. Water Resour. Res. 33 (9), 2165– 2177. Nelsen, R.B., 2006. An Introduction to Copulas. Springer, xx, 315 pp. Oliveira, R., Loucks, D.P., 1997. Operating rules for multireservoir systems. Water Resour. Res. 33 (4), 839–852. Panigrahi, D.P., Mujumdar, P.P., 2000. Reservoir operation modelling with fuzzy logic. Water Resour. Manage. 14 (2), 89–109. Reddy, M.J., Ganguli, P., 2012. Bivariate flood frequency analysis of upper Godavari river flows using archimedean copulas. Water Resour. Manage. 26 (14), 3995– 4018. Salas, J.D., 1993. Analysis and modeling of hydrologic time series. Handbook of Hydrology. Shiau, J.T., 2011. Analytical optimal hedging with explicit incorporation of reservoir release and carryover storage targets. Water Resour. Res. 47 (1), 238–247. Singh, V.P., Zhang, L., 2006. Bivariate flood frequency analysis using the copula method. J. Hydrol. Eng. 11 (2), 150–164. Sklar, A., 1959. Fonctions de repartition a n dimensions et Leurs Marges. Publ. Inst. Stat. Univ. Paris 8, 229–231. Stedinger, J.R., Sule, B.F., Loucks, D.P., 1984. Stochastic dynamic programming models for reservoir operation optimization. Water Resour. Res. 20 (11), 1499– 1505. Su, S.Y., Deininger, R.A., 1972. Generalization of White’s method of successive approximations to periodic markovian decision processes. Oper. Res. 20 (2), 318–326. Tan, Q.F. et al., 2017a. Derivation of optimal joint operating rules for multi-purpose multi-reservoir water-supply system. J. Hydrol. 551, 253–264. Tan, Q.F. et al., 2017b. The dynamic control bound of flood limited water level considering capacity compensation regulation and flood spatial pattern uncertainty. Water Resour. Manage. 31, 1–16. Tejada-Guibert, J.A., Johnson, S.A., Stedinger, J.R., 2010. The value of hydrologic information in stochastic dynamic programming models of a multireservoir system. Water Resour. Res. 31 (10), 2571–2579. Willis, R., Finney, B.A., Chu, W.S., 1984. Monte Carlo optimization for reservoir operation. Water Resour. Res. 20 (9), 1177–1182. Xu, W., Zhang, C., Peng, Y., Fu, G., Zhou, H., 2014. A two stage Bayesian stochastic optimization model for cascaded hydropower systems considering varying uncertainty of flow forecasts. Water Resour. Res. 50 (12), 9267–9286. Young, G.K., 1967. Finding reservoir operating rules. J. Hydraulics Div. 93 (6), 297– 322. You, J.Y., Cai, X., 2008a. Hedging rule for reservoir operations: 1. A theoretical analysis. Water Resour. Res. 44 (1), 186–192. You, J.Y., Cai, X., 2008b. Hedging rule for reservoir operations: 2. A numerical model. Water Resour. Res. 44 (1), 568–569. Yun, R., Singh, V.P., 2008. Multiple duration limited water level and dynamic limited water level for flood control, with implications on water supply. J. Hydrol. 354 (1–4), 160–170. Zhao, J., Cai, X., Wang, Z., 2010. Optimality conditions for a two-stage reservoir operation problem. AGU Fall Meeting, 532–560.