Identifying Efficient Operating Rules for Hydropower Reservoirs Using System Dynamics Approach—A Case Study of Three Gorges Reservoir, China

Abstract

In the long-term operation of hydropower reservoirs, operating rules have been used widely to decide reservoir operation because they can help operators make an approximately optimal decision with limited runoff information. However, the problems faced by reservoir managers is how to make and select an efficient operating rule properly. This study focuses on identifying efficient and reliable operating rules for the long-term operation of hydropower reservoirs using system dynamics (SD) approach. A stochastic hydrological model of reservoir inflow time series was established and used to generate a large number of inflow scenarios. A deterministic optimization operation model of hydropower reservoirs was constructed and then resolved using dynamic programming (DP) algorithm. Simultaneously, within implicit stochastic optimization (ISO) framework, different operating rules were derived using linear fitting methods. Finally, the most efficient one of the existing operating rules was identified based on SD simulation coupled with the operating rules. The Three Gorges Reservoir (TGR) in central China was used as a case study. The results show that the SD simulation is an efficient way to simulate a complicated reservoir system using feedback and causal loops. Moreover, it can directly and efficiently guide reservoir managers to make and identify efficient operating rules for the long-term operation of hydropower reservoirs.

1. Introduction

Reservoir operation is an extremely complicated and comprehensive problem partly due to the uncertainty of reservoir inflow and the stochastic fluctuation of water use demand [1,2,3,4]. In fact, the topic of reservoir operation has been extensively studied and discussed [4,5,6,7,8,9,10,11,12,13,14]. For tackling the uncertainty of reservoir operation, operating rules have been used widely to provide guidelines for reservoir operation decisions because of their capacity of coping with inflow uncertainty and characteristics of easy implementation [1,2,15,16,17,18,19,20]. Here, the operating rules are a set of rules for determining the quantities of water to be stored and to be released from a reservoir under various conditions and demands [9].

Generally, the operating rules can be presented in the form of graphs, tables or functions [21,22] and specify operational decisions (e.g., water release or power output) during each period as a function of appropriate available information (e.g., current or previous reservoir water level, current or previous reservoir inflow and the time of year) [1,21,23]. There are two ways to generate operating rules considering the uncertainty of decision input information—explicit stochastic optimization (ESO) [24] and implicit stochastic optimization (ISO) [17]. These two frameworks deal explicitly and implicitly with the stochasticity of inflow respectively [2,7] and are shown in Figures 2 and 5 in the literature [7]. Bellman stochastic dynamic programming (SDP) is the most famous ESO method [25] and has been the most widely used [26,27]. Unfortunately, the SDP suffers from the curse of dimensionality, namely, the computational complexity increases exponentially with the increase of dimension (i.e., number of reservoirs) [2]. According to the literatures [23,27,28,29], the ISO method is an efficient alternative to ESO for considering uncertainty of reservoir inflow. Compared with ESO, ISO has more computational advantages and has been widely used to derive operating rules for the long-term operation of reservoirs [15,29,30,31,32]. ISO is actually a deterministic optimization approach. But based on the long representative hydrological data (observed or synthetic inflow scenarios) [2], most stochastic aspects of reservoir operation problems are implicitly considered by ISO method [8,33]. As a simple and efficient method, the ISO method was employed to generate operating rules in different functional forms using synthetic inflow generated through hydrologic simulation in this paper.

System dynamics (SD) was firstly proposed by Forrester in the 1958 [34], which is good at simulating complex system by establishing appropriate feedback loops that indicate causal relationship among system elements [35,36]. SD simulation is more flexible and more transparent inherently than conventional simulation method [37] and is helpful for increasing participation and confidence of model users in the model development [37,38]. Therefore, based on SD simulation, modelers and operators can better understand and analyze dynamic behavior of complex systems, such as reservoir systems [35,39]. As an efficient and simple approach, SD approach has been applied in sociology and economics [40,41,42,43] and water resources systems [38,44,45,46] and also become popular in reservoir operation modeling [35,36,37,47].

Using the SD approach, Xu et al. [38] constructed a water resources SD model to analyze the water resource sustainability of the Yellow River in China. King et al. [48] used SD simulations to assess hydropower system safety. Yuan et al. [35] developed a self-optimization system dynamics simulation model of real-time short term cascade hydropower system and extracted operation rules for short-term hydropower reservoir operation by SD simulation. Teegavarapu et al. [47] proposed a simulation model to understand the dynamic behavior of a hydraulically coupled multiple reservoirs system and help make decisions in real-time operation. Felfelani et al. [36] established an SD model to simulate different hedging rules for exploring the most effective operation policy of long-term reservoir water supply dispatching. Ahmad et al. [37] used the SD approach to a reservoir operation problem for flood control and investigated the impacts on flood control capacity of reservoirs by simulation a gated spillway and an un-regulated spillway. These studies indicate that the SD approach is an efficient tool for the simulation of reservoir systems.

Within an ISO framework, given the different function forms of operating rules and the optimal solutions of deterministic reservoir optimization dispatching, the different operating rules of hydropower reservoirs can be obtained using linear or nonlinear fitting methods [15,21,22,31,32,39]. Then, the question is how do we know which one of these operating rules is the most efficient? To solve this problem, a simulation model can be constructed to identify the most efficient one from the known operating rules [29]. There are two ways to constructed the simulation models, traditional method based on high level programming language and the SD method. However, the traditional simulation models developed based on high level programming language are not able to provide insights into the dynamics of the reservoir system and the model structure is rigid and the interaction among different system elements within the models is not transparent [47]. As mentioned above, the SD approach is a promising alternative for simulation of reservoir systems [36,37] and it can model the reservoir operation systems in a realistic way [37]. Moreover, in the process of SD modelling, reservoir operators can effectively participate in the model development, which can enhance their confidence in the SD model [37]. Therefore, inspired by the literature [29], in order to identify the most efficient operating rules for the long-term operation of hydropower reservoirs from different operating rules, an SD simulation model is proposed in this paper. Major contributions are outlined as follows:

(1) An SD simulation model for the long-term operation of hydropower reservoir is established, including a causal loop diagram and a stock flow diagram, which can efficiently describe the dynamic characteristics of the hydropower reservoirs.

(2) The SD simulation model is applied to a real-world case study at Three Gorges Reservoir (TGR), which is the largest and most important reservoir along China’s longest river. Two simple operating rules for long-term operation of TGR are derived by ISO method and simulated using SD model and the most efficient one is identified. In this application, the SD model is coupled with the operating rules, which enables SD simulation to efficiently identify the efficient operating rules for hydropower reservoirs. And the purpose of this application is to demonstrate the role of SD approach in simulating the dynamics of reservoir and identifying efficient operating rules of hydropower reservoirs.

The SD model developed in this study is different from the previous research work [35,36,37,47]. The literature [35] aims to realize self-optimization of real-time short-term operation of cascade hydropower system by the SD simulation. The literature [37] focuses on the simulation of a single multipurpose reservoir for flood management purposes using the SD model. The literature [47] focus on modelling the dynamics of a hydraulically coupled hydropower reservoir system operation, which is able to provide real-time decision support for short-term operation of reservoirs. The literature [36] aims to obtain the effective hedging rules for long-term water supply operation of reservoirs by SD simulation. While the purpose of this study is to identifying the efficient operating rules for long-term operation of hydropower reservoirs by SD simulation coupled with operating rules.

The remainder of this paper is organized as follows—in Section 2, the hydrological simulation model is established to generate synthetic inflow series. The deterministic optimization operation model of hydropower reservoirs is then constructed to obtain optimal operation trajectories with the synthetic inflow series as input. Next, different operating rules are defined. Finally, an SD simulation model is developed to identify the most efficient one from given operating rules. In Section 3, a case study of the TGR in central China is described. The results are then presented and discussed. In Section 4, conclusions are given finally.

2. Methods

As shown in Figure 1, the methodology consists of four modules, the hydrologic time series modelling, the deterministic reservoir optimization operation, the reservoir operating rules derivation and the SD simulation of hydropower reservoir operation respectively. A first-order seasonal autoregressive model (SARM) is established based on the historical reservoir inflow series and used to generate long series reservoir inflow. Next, a bi-objective reservoir optimization operation model is constructed considering maximizing average annual power generation and maximizing generation guarantee rate and DP algorithm is employed to solve the model with the simulated synthetic inflow series as input. Then, within the ISO framework, different linear operating rules are derived by the linear fitting method. Finally, an SD simulation model of long-term operation of hydropower reservoir is developed and the most efficient of the known operating rules is identified ultimately based on simulation results.

2.1. Hydrological Time Series Modelling

Generally, the historical inflow sequence is too short to extract the operating rules reliably. In order to generate a long-term reservoir inflow sequence that reflects both features of historical inflow and characteristic of possible future inflow, a hydrological simulation model needs to be built [49]. Obviously, due to the periodic revolution and rotation of the earth, the climate also shows periodic characteristics. That is to say, the ten-day (in this study, a year is divided into 36 periods) inflow of reservoir have periodic characteristics. Therefore, a first-order SARM is used to characterize the inflow time series of reservoir and to generate a large number of synthetic reservoir inflow scenarios in this study. And the pre-processing of the historical inflow series and the structure of first-order SARM are introduced in this section.

2.1.1. Pre-Processing of Historical Inflow Series

The Pearson type III (P-III) distribution is recommended by the Ministry of Water Resources of China as the distribution of hydrological elements [50]. That is to say, hydrological elements including reservoir inflow obey skewed distribution. However, the processing of skewness distribution is inconvenient. Thus, in order to establish the first-order SARM more conveniently and simulate reservoir inflow effectively, it is necessary to convert the original skewness sequences into the standard normal sequences. And the task consists of two steps—(1) Converting the original skewness inflow sequences into standardized skewness sequences; (2) Converting the standardized skewness sequences into standardized normal sequences. As shown below, step (1) can be accomplished by applying Equation (1) and step 2 can be accomplished by applying Equation (2) [50].

$$X_i(t)=\frac{I_i(t)-\overline{I}(t)}{S_I(t)}$$

(1)

$$Y_i(t)=\frac{6}{C_{SX}(t)}\left\{{\left[\frac{C_{SX}(t)}{2}{X}_i(t)+1\right]}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right.}-1\right\}+\frac{C_{SX}(t)}{6},$$

(2)

where I_i(t) is the inflow of the t-th time period in the i-th year and the time interval of each period is 10 days; $\overline{I}(t)$ is the average inflow of the inflow sequences I(t) = (I₁(t), I₂(t), …, I_T(t)); S_I(t) is the standard deviation of the I(t); X(t) = (X₁(t), X₂(t), …, X_T(t)) is the standardized skewness sequences; C_SX(t) is the skewness coefficient of X(t); and Y(t) = (Y₁(t), Y₂(t), …, Y_T(t)) is the standardized normal sequences.

2.1.2. First-Order Seasonal Autoregressive Model

In order to generate long-term reservoir inflow sequences, there are two simple and effective ways to do this. One is the first-order SARM [51], the other is disaggregation model [52]. This paper focuses on the first way. According to the Y_i(t) in Section 2.1.1, the first-order SARM is as follow:

$$\{\begin{array}{l}{Y}_{i,1}={\phi }_{0,1}+{\phi }_{1,1}{Y}_{i-1,36}+{\epsilon }_{i,1},{\epsilon }_{i,1}~N(0,{\sigma }_1^2)\\ Y_{i,2}={\phi }_{0,2}+{\phi }_{1,2}{Y}_{i,1}+{\epsilon }_{i,2},{\epsilon }_{i,2}~N(0,{\sigma }_2^2)\\ ⋮⋮⋮⋮⋮\\ Y_{i,36}={\phi }_{0,36}+{\phi }_{1,36}{Y}_{i,35}+{\epsilon }_{i,36},{\epsilon }_{i,36}~N(0,{\sigma }_{36}^2)\end{array},$$

(3)

where, Y_i,t and Y_i(t) are equivalent; ${\phi }_{0,t}$ and ${\phi }_{1,t}$ are the regression coefficients of the t-th time period; ${\epsilon }_{i,t}$ is the random error term of the t-th time period in the i-th year and obeys a normal distribution with a mean of 0 and a variance of ${\sigma }_t^2$. Obviously, Model (3) has 108 parameters.

After parameters ${\phi }_{0,t}$, ${\phi }_{1,t}$ and ${\sigma }_t^2$ of the Model (3) are estimated according to the observation value of reservoir inflow, we can use the Model (3) to generate synthetic inflow. Be careful, the simulation value from the Model (3) is Y and the real simulation value of reservoir inflow is Y’s inverse transformation, that is, the inverse transformation of Equations (1) and (2).

2.2. Deterministic Reservoir Optimization Operation Model

Within an ISO framework, the optimal solutions of deterministic reservoir optimization operation is requisite for the deducing of operating rules [17,53]. Hence, a deterministic optimization operation model of hydropower reservoirs is constructed in this section and the model objectives, constraints and the solution algorithm are introduced.

2.2.1. Objective Function

(1) Maximize the average annual power generation [15]

$$\max E=\frac{1}{T}{\displaystyle \sum _{i=1}^T{\displaystyle \sum _{j=1}^M{N}_{i,j}\Delta t_{i,j}}},$$

(4)

where E is the average annual power generation of reservoirs; T and M are the number of years and the number of time periods per year, respectively; $\Delta t_{i,j}$ is the time interval of j-th time period in the i-th year; N_i,j is the power output of j-th time period in the i-th year and is determined as follow:

$$N_{i,j}=K{Q}_{i,j}{H}_{i,j},$$

(5)

where K is the comprehensive generation coefficient; H_i,j is the net water head of the j-th time period in the i-th year; Q_i,j is the water release passing turbines of the j-th time period in the i-th year, constrained by the turbine's water passing capacity and generator’s power generating capacity.

(2) Maximize the guarantee rate of hydropower generation [15]

$$\max P_{\gamma }=\frac{1}{T\times M}{\displaystyle \sum _{i=1}^T{\displaystyle \sum _{j=1}^M\Delta m_{i,j}(N_{i,j})}},$$

(6)

$$\Delta m_{i,j}(N_{i,j})=\{\begin{array}{l}1N_{i,j}\ge N_g\\ 0N_{i,j}<N_g\end{array}$$

(7)

where $P_{\gamma }$ is the guarantee rate of hydropower generation; N_g is the guaranteed output and other symbols are the same as the above.

2.2.2. Constraints

The model constraints are as follows:

(1) Reservoir storage, water release and power output constraints [54]

$$\{\begin{array}{l}{V}_j^{\min }\le V_{i,j}\le V_j^{\max }\\ R_j^{\min }\le R_{i,j}\le R_j^{\max }\\ N_j^{\min }\le N_{i,j}\le N_j^{\max }\end{array}$$

(8)

$$R_{i,j}=Q_{i,j}+S_{i,j}$$

(9)

where V_i,j is the reservoir storage of the j-th time period in the i-th year; R_i,j is the water release of the j-th time period in the i-th year; S_i,j is the spill water not passing turbines during the j-th time period in the i-th year; $V_j^{\min }$ and $V_j^{\max }$ are the minimum and maximum storage in the j-th time period, respectively; $R_j^{\min }$ and $R_j^{\max }$ are the minimum and maximum water release in the j-th time period, respectively; $N_j^{\min }$ and $N_j^{\max }$ are the minimum and maximum power output in the j-th time period, respectively.

(2) Initial and terminal water level constraints [54]

$$V_{i,0}=V_i^{\mathrm{begin}}，V_{i,M}=V_i^{\mathrm{end}},$$

(10)

where $V_i^{\mathrm{begin}}$ and $V_i^{\mathrm{end}}$ are the initial and terminal reservoir storages, respectively.

(3) Water balance equation [1,54]

$$V_{i,j}=\{\begin{array}{l}{V}_{i,j-1}+(I_{i,j-1}-R_{i,j-1})\Delta t_{i,j-1},j=2,3,...,M\\ V_{i-1,M}+(I_{i-1,M}-R_{i-1,M})\Delta t_{i-1,M},j=1\end{array},$$

(11)

where I_i,j is the reservoir inflow of the j-th time period in the i-th year. Because the water loss caused by reservoir seepage and evaporation is negligible and therefore not considered here.

2.2.3. Solution

To solve the bi-objective problem, we can simplify the bi-objective model into a single-objective model by adding a penalty function [15], which can be described as follows:

$$\max E=\frac{1}{T}{\displaystyle \sum _{i=1}^T{\displaystyle \sum _{j=1}^M(N_{i,j}-\alpha \beta )\Delta t_{i,j}}}$$

(12)

$$\{\begin{array}{l}\alpha =\{\begin{array}{l}0，N_{i,j}\ge N_g\\ 1，N_{i,j}<N_g\end{array}\\ \beta =\lambda {(N_{i,j}-N_g)}^{\gamma }\end{array}$$

(13)

where $\alpha$, $\lambda$ and $\gamma$ are the penalty coefficient; $\beta$ is the penalty quantity. Obviously, the above model is a deterministic reservoir optimization operation problem when the reservoir inflow, the reservoir initial level and the terminal water level all are known. Therefore, DP algorithm is able to solve the optimization model conveniently. For more information about DP, please refer to the Bellman’s research work [55].

2.3. Reservoir Operating Rules Derivation

Since the implicit stochastic optimization (ISO) method was proposed [17], which is an efficient alternative to explicit stochastic optimization [27,28], various forms of operating rules have been developed, including linear operating rules [17,53] and nonlinear operating rules, such as artificial neural networks [2]. According to literature [1,16,23], among these operating rules, linear operating rules are the simplest, simultaneously can effectively improve reservoir operation performance. Hence, the simple linear operating rules are used in this study. Actually, within an ISO framework, we can obtain different operating rules including linear or nonlinear by using different data mining method, but, without loss of generality, the linear operating rules are used as the example to verify the performance of SD simulation model in this paper. In this study, only two forms of linear operating rules are introduced and used.

(1) Operating rules considering the decision information including predicted inflow of reservoir in the current time period and reservoir storage in the current time period, abbreviated as operating rule A.

$$R_{i,j}=a_j(I_{i,j}+{V_{i,j}/\Delta t}_{i,j})+b_j,\hspace{1em}j=1,2,...,M.$$

(14)

(2) Operating rules considering the decision information including actual inflows of reservoir in the last time period and reservoir storage in the current time period, abbreviated as operating rule B.

$$R_{i,j}=a_j(I_{i,j-1}+{V_{i,j}/\Delta t}_{i,j})+b_j,\hspace{1em}j=1,2,...,M,$$

(15)

where R_i,j is the reservoir water release to be determined by operating rules; $(I_{i,j}+{V_{i,j}/\Delta t}_{i,j})$ and $(I_{i,j-1}+{V_{i,j}/\Delta t}_{i,j})$ represent approximately the available water of the j-th time period in the i-th year; a_j and b_j are operating rule parameters to be determined by linear fitting method based on the deterministic optimization operation trajectories (reservoir inflow process, reservoir storage process and water release process).

An intuitive description of linear operating rules is shown in Figure 2 [16,23]. The grey area is the feasible region for decision-making of reservoir water release, which is enclosed by the reservoir storage bounds and water release bounds. The three lines with different colors represent operating rules of different time periods, such as the first, second and third time period of 36 time periods in each year. In addition, when available water is less than the critical value, water release is determined by using the rule 2. While available water is greater than critical value, water is released according to the operating rule 1. Here, the form of rules 1 and 2 is determined by Equation (14) or Equation (15). That is to say, rules 1 and 2 are concrete implementations of rule A or rule B.

In order to identify the most efficient operating rules for the long-term operation of hydropower reservoirs from different operating rules, an SD simulation model is introduced in the Section 2.4.

2.4. System Dynamics Simulation Model of Hydropower Reservoir Operation

SD simulation, a simulation technique based on feedback of system elements, has attracted much attention from water resources researchers and is considered as one of the efficient tools in reservoir operation modeling [36]. Therefore, in order to evaluate different operating rules effectively, an SD simulation model is developed as follows.

The process of SD modeling includes the following steps—define the dynamic hypothesis; develop the causal loop diagrams; construct the stock-flow diagrams [43]. In addition, the SD model mainly includes three kinds of variables, stocks, flow and auxiliary variables, respectively. Stocks are accumulated variables or state variables changed with the flow and represent the state of dynamic systems, such as the storage of reservoir systems. Flow is rate variable like as inflow and outflow of reservoir systems. Auxiliary variables are equivalent to intermediate variables. Just as Professor Forrester used water in bathtubs to illustrate stocks and flow when he created system dynamics, reservoir systems have a typical structure of stocks and flow like bathtubs.

The dynamic characteristics of reservoirs are obvious. Reservoir storage changes with inflow and outflow. When net inflow is positive, reservoir storage increases, while when net inflow is negative, reservoir storage decreases. Based on the above dynamic characteristics, we develop the causal loop diagram (Figure 3) and the stock-flow diagram (Figure 4) for hydropower reservoirs using the Vensim Soft [56], which is an object-oriented simulation environment based on the principle of SD and allows you to efficiently and simply simulate, analyze and optimize dynamic systems.

In Figure 4 and Figure 5, blue or purple lines with an arrow are causal chains or connecter and represent the cause relationship among variables. In Figure 4, the sign (+ or −) at the end of the lines indicate the type of effect one variable has on the other, positive feedback or negative feedback. For example, an increase in the independent variable, reservoir level, leads to an increase in the dependent variable, head. The sign * represents that water release is determined by reservoir storage, inflow and operating rules and the polarity of causal chain among them is uncertain due to the complexity of operating rules. Actually, the stock-flow diagram here is simplified. Detailed stock-flow diagram is shown in the case study. In addition, in Figure 4 and Figure 5, all variables have physical explanations and the model structure also satisfies the physical mechanism (water balance) strictly. Therefore, the verification of model parameters and structure can be simplified.

Finally, it is worth emphasizing that the reservoir inflow and operating rules in stock-flow diagram and causal loop diagram are exogenous variables and others are endogenous variables. Moreover, these two variables have the characteristics of uncertainty. That is to say, the input of the SD model for hydropower reservoirs is uncertain. While its structure and parameters are certain, because it satisfies the physical mechanism (water balance) strictly.

3. Case Study

3.1. The Three Gorges Reservoir

The TGR is the largest and the most crucial comprehensive reservoir, which lies in the upper reaches of the Yangtze River, the longest river of China (Figure 5). The average annual inflow at dam site is approximately equal to 4.51 × 10¹¹ m³/year. The catchment area of the reservoir is about 10⁶ km². With a flood control capacity of 2.215 × 10¹⁰ m³ and a usable storage capacity of 1.65 × 10¹⁰ m³, the TGR plays an important role in flood control, hydropower generation, water supply and navigation. The hydropower station is equipped with 32 units of 7 × 10⁵ kW and 2 units of 5 × 10⁴ kW, with a total installed capacity of 2.25 × 10⁷ kW.

The 130-year reservoir inflow records (from 1882 to 2011) from the Yichang hydrologic station, which is located approximately 40 km downstream of the TGR, are used to estimate the parameters of the first-order SARM. And the inflow sequence from 1996 to 2011 is used as the simulation input of the SD simulation model. All streamflow records have been restored, therefor, the streamflow statistical feature is consistent. For this case study, the time interval of deterministic optimization operation model could be selected as a daily time step, 10 days’ time step or a monthly time step. But a 10 days’ time step (actually, 8 or 9 days (February), 10 days or 11 days, a traditional Chinese measure of time) is more appropriate than others, considering the daily time step increases the computational complexity and the monthly time step ignores the runoff variation within a month [1,16].

In addition, other parameters of the deterministic reservoir optimization operation model are set as follows. The initial and terminal water levels of the TGR are set to the flood limit level of 145 m. For the flood season (11 June to 10 September), the water level is controlled as the flood limit level. For the non-flood season, the minimum and maximum water level are set as 145 m and 175 m, respectively. The minimum discharge is set as 5000 m³/s to guarantee downstream navigation and ecological needs. The comprehensive generation coefficient K in Equation (5) is set as 8.8.

The other data sets of the TGR including reservoir level-storage curves, reservoir outflow- tail water elevation curves, water head-power generation water release ability curves and reservoir level- water release ability curves are used for this case study.

3.2. Results

3.2.1. Simulated Reservoir Inflow Series Based on the First-Order SARM

The parameters of the first-order SARM (Equation (3)) is calibrated based on the 130-year reservoir inflow records (from 1882 to 2011) from the Yichang hydrologic station (restricted to article space, the 108 model parameters of SARM are not listed in detail). Then, 10,000 years of synthetic inflow series is generated with a ten-day time step using the calibrated model (the 10,000 years of synthetic inflow sequence generated here is only used to verify the validity of the SARM model). The basic statistical parameters of the observed and simulated inflow data, including the mean value, standard deviation, skewness coefficient and lag-1 correlation coefficient, are calculated, as shown in Figure 6. It can be seen that the mean value, standard deviation, skewness and lag-1 correlation coefficient of the simulated data fit those of observed data well. Obviously, there is no significant difference between the main statistical parameters of the simulated sequence and the corresponding statistical parameters of the observed inflow sequence. Therefore, the first-order SARM passes the practical test and is applicable.

The 1000 years of synthetic inflow sequence (Figure 7) is regenerated after the verification of SARM is completed, which is used as input to deterministic reservoir operation model. Of course, ISO approach is more effective considering a longer inflow series but it will take a long time to solve the reservoir optimization operation model. Actually, 1000 years of synthetic inflow series is enough to reflect uncertainty of reservoir inflow.

3.2.2. Optimal Solution of Deterministic Optimization Operation and Linear Operating Rules

With the 1000 years of synthetic inflow series (Figure 7) as input, the deterministic reservoir optimization scheduling model of TGR is solved by using DP algorithm. Then the optimal trajectories (water release process, reservoir storage process and available water process) are obtained, which are used to derive operating rules A and B. In order to verify the rationality of linear operating rules A and B, the scatter plots of optimal water release versus corresponding available water is drawn in Figure 8 and Figure 9. Note, the available water is approximately represented in two ways, $I_{i,j}+{V_{i,j}/\Delta t}_{i,j}$ and $I_{i,j-1}+{V_{i,j}/\Delta t}_{i,j}$, which has been discussed in Section 2.3.

Figure 8a and Figure 9a show correlation between optimal water release and available water in each of 36 time periods (a year is divided into 36 time periods), which show correlation between optimal water release and available water in a comprehensive way. Figure 8b and Figure 9b show this correlation in each of several time periods, which is able to show this correlation more clearly than Figure 8a and Figure 9a.

In Figure 8a and Figure 9a, the purple plus signs and green plus signs denote the correlation between optimal water release and available water during flood season and non-flood season, respectively. In Figure 8b, the blue plus signs denote this correlation during the last 11 days of October and the red plus signs denote this correlation during the first 10 days of November. In Figure 9b, the black plus signs denote this correlation during the last 10 days of June, the blue plus signs denote this correlation during the last 10 days of September and the red plus signs denote this correlation during the last 10 days of April. Obviously, according to Figure 8 and Figure 9, there is a strong linear correlation between optimal water release and available water no matter the available water is approximated by $I_{i,j}+{V_{i,j}/\Delta t}_{i,j}$ or $I_{i,j-1}+{V_{i,j}/\Delta t}_{i,j}$. That is to say, linear operating rule A and linear operating rule B both are rational and applicable for long-term operation of TGR.

According to the optimal operation results shown in Figure 8 and Figure 9, operating rule A and operating rule B are derived respectively (

Table 1. Derived operating rules A and B for TGR in January and August.

Time	Operating Rules	Rule 1		Critical Value (m3/s)	Rule 2
		aj	bj		aj	bj
The first 10 days of January	Operating rule A	1.00000	−45,486.1	50,640	−0.03103	6742.2
	Operating rule B	0.01589	4428.9	50,160	−0.02303	6365.6
The middle 10 days of January	Operating rule A	1.00000	−45,486.1	50,700	−0.03231	6810.0
	Operating rule B	0.68758	−29,870.7	50,910	−0.02608	6518.2
The last 11 days of January	Operating rule A	1.00000	−41,351.0	46,520	−0.03698	6893.4
	Operating rule B	0.17476	−2948.9	46,270	−0.03589	6853.8
The first 10 days of August	Operating rule A	0.99439	−19,586.6	/	/	/
	Operating rule B	0.48547	4707.5	/	/	/
The middle 10 days of August	Operating rule A	1.00000	−19,849.5	/	/	/
	Operating rule B	0.56703	1026.1	/	/	/
The last 11 days of August	Operating rule A	1.00000	−18,045.0	/	/	/
	Operating rule B	0.49145	4454.7	/	/	/

). They will be coupled with an SD model to simulate the long-term power generation operation of TGR in next section. Due to limited space, Table 1 shows only the operating rules A and B for TGR in the January and August. For operating rules in other time periods, the contents are similar to Table 1. Finally, it should be emphasized that both operating rules A and B are linear in flood season and piecewise linear in non-flood season, which can be proved by Figure 8 and Figure 9.

3.2.3. Identifying the More Efficient Operating Rules Based on System Dynamics Simulation

According to the causal loop diagram (Figure 3) and rough stock-flow diagram (Figure 4), a practical SD simulation model is constructed (Figure 10). In Figure 10, the red part represents any operating rules of hydropower reservoirs. But here, no loss of generality, we only construct this part in the form of linear operating rules that are described in Equations (14) and (15). Coefficients a-f and b-f illustrate the operating rules of flood season, that is, the operating rules of flood season are linear. And coefficients a1-nonf, b1-nonf, a2-nonf, b2-nonf and critical point indicate the operating rules of non-flood season, that is, the operating rules of non-flood season are piecewise linear. In addition, the meaning of other variables in the stock-flow diagram is just like their names, for example, get reservoir level by storage is a table function, power generation and reservoir storage are stocks or state variables, reservoir inflow and stage power generation both are flows or rate variables and others are auxiliary variables.

In order to identify the most efficient one from operating rules A and B, the reservoir inflow that is exogenous variable of the SD model is set as the observed inflow of TGR from 1996 to 2011. Obviously, the time boundary of the SD model is 1996–2011 and the spatial boundary of the model is TGR. Coupling operating rule A and operating rule B with the SD model respectively, we can get the simulation results shown in

Table 2. Comparison of simulation results of three schemes (the operation horizon is 1996–2011).

Scheme	Average Annual Power Generation (108 kWh)	Guarantee Rate (%)	Spill Water (106 m3)
Operating rule A	849.67	95.31	34.72
Operating rule B	846.22	93.06	34.72
DP algorithm	853.84	98.09	0

. In addition, Table 2 also contains deterministic optimization operation results of TGR with operation horizon of 1996–2011.

As shown in Table 2, according to the two indices of average annual power generation and guarantee rate, operating rule A is identified as more efficient than operating rule B by SD simulation. In addition, the two operating rules both are worse than deterministic optimization operation. In fact, this is inevitable, because deterministic reservoir operation utilizes all inflow information (that is, all inflow information is known), however, operating rules only can utilize available information, such as previous inflow and reservoir level information or forecast information (that is, the available information is limited). Through SD simulation, we can draw the conclusion that operating rule A is more efficient than operating rule B for long-term operation of TGR. More importantly, due to the transparency of SD model, not only the modeler is convinced of such results but also the reservoir operators are confident of this conclusion.

To illustrate the difference between operating rules A and B in more detail, operation results of two operating rules for one hydrological year are obtained based on SD simulation (here, the hydrological year refers to the period from 11 June 2001 to 11 June 2002) and shown in

Table 3. Operation results of operating rules A and B (the operating horizon is 11 June 2001 to 11 June 2002).

Scheme	Annual Power Generation (108 kWh)	Guarantee Rate (%)	Spill Water (106 m3)
Operating rule A	824.67	86.5	0
Operating rule B	813.99	86.5	0

. Simultaneously, water release, reservoir storage and power output process of TGR are drawn using Vensim software (Figure 11, Figure 12 and Figure 13).

Table 3 shows that the annual power generation of operating rule A is larger than that of operating rule B, while the other indicator is not different. This result is explained as follows. As shown in Figure 11, Figure 12 and Figure 13, in the flood season (from 1-th to 9-th time period), operation results of the two operating rules are the same. However, at the beginning of non-flood season (10-th to 13-th time period), the water release decision of operating rule A is less than that of operating rule B and thus the power output of operating rule A is smaller than that of operating rule B. But rule B will consume the water in the TGR prematurely (in fact, during beginning of non-flood season, the reservoir level is low and such a water release decision is not conducive to power generation) and affect the annual power generation. Finally, at the end of non-flood season (30-th to 36-th time period), the water release decision of rule A is greater than that of rule B, while maintaining a higher water level than that of rule B and so the annual power generation is greater than that of rule B. In short, the decision-making mode of rule A and rule B in flood season is the same but there are differences in non-flood season. At the beginning of the non-flood season, the reservoir water level of TGR is low, rule A keeps a smaller water release decision while rule B keeps a larger water discharge decision. At the end of the non-flood season, the reservoir water level of TGR is high, rule A keeps a greater water release decision while rule B keeps a smaller water discharge decision. Obviously, according to the above discussion, operating rule A is more beneficial to hydroelectric power generation of TGR than operating rule B. That is to say, the most efficient operating rules are identified successfully from given operating rules A and B by SD simulation.

3.3. Discussion of Results

Coupling operating rules with SD model to identify the most efficient one from given operating rules A and B for long-term operation of TGR is a typical application. Through SD simulation, operating rule A is identified as the most efficient one from operating rule A and rule B derived by ISO method. This application research demonstrates the role of SD approach in simulating the dynamics of reservoir and identifying efficient operating rules of hydropower reservoirs. In addition, compared with the traditional simulation model based on high level programming language, the transparency of SD model makes the identified efficient operating rules more convincing to model users or reservoir operators.

It is worth emphasizing that SD simulation is mainly used to help reservoir operators to make and identify more efficient operating rules, while the operating rules are decision-making tools to guide reservoir operators to make water release decisions on operation of the hydropower reservoir. In this study, the function of SD model is to identify this efficient decision-making tool. Obviously, according to the research results in this paper, SD approach can be competent for this task.

In addition, from the study above, there may be following several problems. Firstly, can the first-order SARM effectively simulate reservoir inflow? Obviously, from the main statistical parameters (mean, standard deviation, skewness coefficient, lag-1 correlation coefficient) (Figure 6), this method is efficient but it is not qualified for the higher order statistical parameters. Fortunately, the first-order model can satisfy this research. Secondly, the optimization operation of hydropower reservoir system is a complex non-linear problem, as shown in Section 2.2. Therefore, the operating rules function of reservoir operation are essentially non-linear. But according to Figure 8 and Figure 9, the linear operating rules are appropriate for the TGR. Despite the above problems, it does not affect the focus of this paper, because operating rules are just an exogenous input of SD model. No matter what the exogenous input is, the SD model can identify its advantages and disadvantages through simulation of reservoir operation and then help reservoir operators to make and identify more efficient operating rules for hydropower reservoirs.

There is a problem that is also worth explaining, that is, model uncertainty. Generally speaking, models involve three parts—model inputs, model itself and model outputs. Uncertainty may lie in any of the three parts. In this paper, For the SD model proposed in this paper, the model inputs include reservoir inflow and operating rules and the model outputs can be any variable you want to observe. Because the SD model itself (Figure 10) strictly satisfies the physical mechanism of reservoir operation (i.e., water balance principle), there is no uncertainty in the model. Therefore, the uncertainty mainly exists in model inputs, which makes the model outputs uncertain. The uncertainty of model inputs including reservoir inflow and operating rules is discussed below.

The uncertainty of operating rules includes two aspects—structural uncertainty and parameter uncertainty. The structure of operating rules can be linear or nonlinear. For example, in this paper, the structure of operating rules A and B is assumed to be linear. Because the operating rules are generated by mining the deterministic optimization operation trajectories of reservoirs, the problems of parameter uncertainty is also inevitable. However, the uncertainty in operating rules is the problem that should be dealt with when making the operating rules, not the problems that this paper should pay attention to. That is to say, as long as there is an operating rule that can effectively consider the uncertainty, the SD model proposed in this paper can be coupled with it to simulate the operation of reservoir, so as to evaluate its performance of guiding the operation of reservoir.

In this paper, as an exogenous input for the SD model, the future value of reservoir inflow is indeed uncertain. In order to avoid this problem, this paper uses historical reservoir inflow as simulation input of the SD simulation. Of course, if we can build a probability prediction model to describe the uncertainty of reservoir inflow, the SD model will be better but it will be a tough job.

Finally, combined with the above research results, we can get the following conclusions—(1) for the TGR, the operating rule A is identified as more efficient than operating rule B using SD simulation; (2) through the graphic analysis (Figure 11, Figure 12 and Figure 13), it is suggested that TGR should be scheduled according to the following basic rules, that is, at the early non-flood season, the water release should be reduced appropriately and at the end of non-flood season, the water release should be increased appropriately, so as to maintain high head operation of the reservoir throughout year and increase the power generation.

4. Conclusions

This study focuses on identifying an efficient operating rule for long-term operation of hydropower reservoirs by SD simulation coupled with operating rules. The first-order SARM (hydrological simulation model) is established and its parameters are estimated based on historical inflow records. Simultaneously, the synthetic inflow series is generated with a ten-day time step using the calibrated hydrological simulation model. The deterministic reservoir optimization model with two-objective function is constructed. And the DP algorithm is used to solve the optimal model based on synthetic inflow series. Within an ISO framework, the linear fitting method is employed to fit the operating rules A and B according to the optimal solutions obtained by DP. Then, the SD simulation model for long-term operation of hydropower reservoirs is developed to identify the most efficient operating rules for reservoirs from given operating rules. Finally, the SD simulation model is applied to a real-world case study at Three Gorges Reservoir (TGR). From this study, conclusions can be drawn as follows—(1) the concept of stocks and flow in SD is very consistent with variables of reservoirs system and the causal relationship among variables about reservoirs can be clearly described by SD model; (2) by coupling with the operating rules, the SD simulation can effectively simulate the dynamics of hydropower reservoirs and identify efficient operating rules for long-term operation of hydropower reservoirs; (3) compared with the traditional simulation models based on high level programming language, the transparency of SD model makes the identified efficient operating rules more convincing to model users or reservoir operators.