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

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.


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 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.

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.

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.

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].

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.

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]. 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; I(t) is the average inflow of the inflow sequences I(t) = (I 1 (t), I 2 (t), . . . , I T (t)); S I (t) is the standard deviation of the I(t); X(t) = (X 1 (t), X 2 (t), . . . , X T (t)) is the standardized skewness sequences; C SX (t) is the skewness coefficient of X(t); and Y(t) = (Y 1 (t), Y 2 (t), . . . , Y T (t)) is the standardized normal sequences.

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: where, Y i,t and Y i (t) are equivalent; ϕ 0,t and ϕ 1,t are the regression coefficients of the t-th time period; ε 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 σ 2 t . Obviously, Model (3) has 108 parameters. After parameters ϕ 0,t , ϕ 1,t and σ 2 t 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).

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.

Objective Function
(1) Maximize the average annual power generation [15] 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; ∆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: 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] where P γ is the guarantee rate of hydropower generation; N g is the guaranteed output and other symbols are the same as the above.

Constraints
The model constraints are as follows: (1) Reservoir storage, water release and power output constraints [54] 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 min where V 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.

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: where α, λ and γ are the penalty coefficient; β 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].

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.
(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.
where R i,j is the reservoir water release to be determined by operating rules; 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.

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.

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 Figures 4 and 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 stockflow diagram is shown in the case study. In addition, in Figures 4 and 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.    In Figures 4 and 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 stockflow diagram is shown in the case study. In addition, in Figures 4 and 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.  In Figures 4 and 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 Figures 4 and 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.

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 11 m 3 /year. The catchment area of the reservoir is about 10 6 km 2 . With a flood control capacity of 2.215 × 10 10 m 3 and a usable storage capacity of 1.65 × 10 10 m 3 , 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 5 kW and 2 units of 5 × 10 4 kW, with a total installed capacity of 2.25 × 10 7 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 3 /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 levelwater release ability curves are used for this case study.
Water 2019, 11, x FOR PEER REVIEW 10 of 20

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 11 m 3 /year. The catchment area of the reservoir is about 10 6 km 2 . With a flood control capacity of 2.215 × 10 10 m 3 and a usable storage capacity of 1.65 × 10 10 m 3 , 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 5 kW and 2 units of 5 × 10 4 kW, with a total installed capacity of 2.25 × 10 7 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 3 /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 levelwater release ability curves are used for this case study.

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

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. 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.  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.
Water 2019, 11, x FOR PEER REVIEW 11 of 20 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.

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 Figures 8 and 9. Note, the available water is approximately represented in two ways, I i,j + V i,j /∆t i,j and I i,j−1 + V i,j /∆t i,j , which has been discussed in Section 2.3.
Water 2019, 11, x FOR PEER REVIEW 12 of 20

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 Figures 8 and 9. Note, the available water is approximately represented in two ways, , Δ , which has been discussed in Section 2.3.    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 Figures 8 and 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 /∆t i,j or I i,j−1 + V i,j /∆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.
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. (b) 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. Figures 8a and 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. Figures 8b and 9b show this correlation in each of several time periods, which is able to show this correlation more clearly than Figures 8a and 9a.
In Figures 8a and 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 Figures 8 and 9, there is a strong linear correlation between optimal water release and available water no matter the available water is approximated by 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 Figures 8 and 9, operating rule A and operating rule B are derived respectively (Table 1). 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 Figures 8 and 9. Figure 9. Scatter plots of optimal water release versus available water (I i, j−1 + V i, j /∆t i, j ). (a) 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. (b) 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.
According to the optimal operation results shown in Figures 8 and 9, operating rule A and operating rule B are derived respectively (Table 1). 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 Figures 8 and 9.

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 Figure10, 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 nonflood 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 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. 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. Simultaneously, water release, reservoir storage and power output process of TGR are drawn using Vensim software (Figures 11-13). 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. Simultaneously, water release, reservoir storage and power output process of TGR are drawn using Vensim software (Figures 11-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 Figures 11-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 decisionmaking 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.  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 Figures 11-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.

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 Figures 8 and 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 (Figures 11-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.

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.