Self-optimization system dynamics simulation of reservoir operating rules

Operating rules have been used widely in the reservoir long-term operation duo to its characteristics of coping with inflow uncertainty and easy implementation. And implicit stochastic optimization (ISO) has been widely applied to derive reservoir operation rules, based on linear regression or nonlinear fitting method. However, the maximum goodness-of-fit criterion of fitting method may be unreliable to determine the effective rules. Therefore, this paper develops a self-optimization system dynamics (SD) simulation of reservoir operation for optimizing the operating rules, by taking advantages of feedback loops in SD simulation. A deterministic optimization operation model is firstly established, and then resolved using dynamic programming (DP). Simultaneously, the initial operating rules (IOR) are derived using the linear fitting method. Finally, the refined optimal operating rules (OOR) are obtained by improving the IOR based on the self-optimization SD simulation. China’s Three Gorges Reservoir is used as a case study. The results show that the SD simulation is competent in simulating a complicated hydropower system with feedback and causal loops. Moreover, it makes a contribution to improve the IOR derived by fitting method within an ISO frame. And the OOR improve effectively the guarantee rate of power generation on the premise of ensuring power generation.


Introduction
Reservoir operation is an extremely complicated partly due to the uncertainty of reservoir inflow and the stochastic fluctuation of load demand [1]. Operating rules have been used widely in the reservoir optimization operation. And the operating rules generally specify decision variables (water discharge or power output) as a funtion of available information, such as current water level and reservoir inflow [1][2][3][4]. Taking random factors into account, especially stochastic nature of inflow, the explicit stochastic optimization (ESO) has been proposed, and often used in studies to derive the operating rules [5]. However, the ESO suffers from curse of dimensionality. Therefore, the implicit stochastic optimization (ISO) is an efficient alternative to ESO. ISO is in fact a deterministic approach, but based on the long representative hydrologiclal data, the past records or the synthetic streaflow, most stochastic aspects of reservoir operation problem are implicitly considered [2,5]. As a simple and efficient method, ISO has been widely used to deduce opermal reservoir operating rules.
Within an ISO framework, fitting method or simulation-based optimization can be employed to derive reservoir operating rules. Nevertheless, the fitting method based on maximum goodness-of-fit criterion is not always reliable or appropriate to determine the operating rules parameters [4,6] because of partly the existence of outliers. Therefore, a more feasible approach which combines simulation and optimization has been suggested and used to derive operating rules [7]. But the conventional simulation-based optimization model developed using high level programming languages do not provide insights into the dynamics of reservoir system [8]. Beside, for the conventional simulation-based optimization, model user or operator can't easily participate in model construction. Therefore, there is still a gap between theoretical research and practice.
System Dynamics (SD) was first proposed by Forrester in the 1950s, which is good at simulating complex system by establishing appropriate feedback loops that indicate causal relationship among system elements [9]. SD simulation is more flexible and more transparent inherently than conventional simulation method, and is helpful for increasing participation of model users in the model development [8]. Therefore, based on SD simulation, modeler or operator can better understand and analyze dynamic behavior of complex system, such as reservoir operation system. Although SD has been widely used in complex systems, such as society, economy [10] and water resources [8,9,11,12], there are few studies on operating rules optimization in SD simulation. Hence, in this paper, for improving the operating rules derived by fitting method, a selfoptimization SD simulation of reservoir operating rules is developed.
The remainder of this paper is organized as follows: In Section 2, the deterministic reservoir operation model, reservoir operating rules derivation and self-optimization SD dynamics simulation model are described. Section 3 describes a case study application to China's Three Gorges Reservoir and results are discussed. In Section 4, conclusions are given finally.

Methodology
As shown in Fig. 1, the methodology consists of three components, the deterministic reservoir optimization operation, the linear operating rules derivation, and the self-optimization SD simulation of reservoir operating rules respectively. A two-objective reservoir operation model is constructed with the maximum average annual power generating and the maximum guarantee rate, and dynamic programming (DP) algorithm is used to solve the model with the historical reservoir inflow series as the input. Then, the linear operating rules (as the initial operating rules, short as IOR) are deduced by the linear fitting method. Finally, the IOR are improved through self-optimization SD simulation, and the optimal operating rules (OOR) are identified.

Deterministic reservoir operation
The optimal solution of deterministic reservoir operation is essential for the deducing of linear operating rules using the linear regression technique. Therefore, the deterministic reservoir optimization operation model is established in this section, and the model objectives, constraints and the solution algorithm are introduced.

Objective function
(1) Maximize the average annual power generation , , where E is the average annual power generation; T and M are the number of years and the number of time steps per years, respectively; , i j t Δ is the interval of time period; N i,j is the power output (kW) during time period j in the ith year, and is calculated as follow: where K is the comprehensive efficiency coefficient; H i,j and Q i,j are the effective water head and the water discharge, respectively; 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 constraints are described as follows: (1) Reservoir water storage, water discharge and power output constraints [13] min max , min max , min max , where V i,j is the reservoir storage during time period j in the ith year; Q i,j is the water discharge during time period j in the ith year; N i,j is the power output during time period j in the ith year; min (2) Initial and terminal water level constraints [13] begin end ,0 , where begin i V and end i V are the initial and terminal reservoir storage, respectively.
(3) Water balance equation [13] , 1 , where I i,j is the reservoir inflow during time period j in the ith year. Because the water loss caused by reservoir seepage and evaporation is negligible and therefore not considered here.

Solution
In order to solve the two-objective problem, we can simplify the two-objective model into a single-objective model by adding a penalty function, which can be expressed as follows: where α , λ and γ are all the penalty coefficient; β is the penalty quantity. The above model is a deterministic reservoir optimal operation problem when both the reservoir inflow and the reservoir initial and final water level are known, so it can be solved directly by DP algorithm.

Reservoir operating rules derivation
As the linear operating rules are the simplest and can effectively promote the reservoir optimization operation, therefor, the two-parameter linear operating rules are adopted here and as follow: where Q i,j is the reservoir discharge to be determined; is the available water; a j and b j are the operating rules parameters to be determined firstly using linear regression and to be improved then using selfoptimization SD simulation.

Self-optimization system dynamics simulation
In order to overcome the defects of linear fitting method and traditional simulation-based optimization in derivation of the operating rules, such as low efficiency, inaccuracy, unreliability and opacity, a self-optimization SD simulation model is proposed and established and described as follows.

System dynamics simulation
Vensim is an object-oriented simulation environment based on the principle of SD, which allows you to efficiently and simply simulate, analyse, and optimize dynamic systems [9]. The Vensim Software provides several general modelling blocks, which are stock, flow and connecter, respectively. Flow is rate variable like as inflow and outflow of reservoir system, and stock is accumulated variables or state variable changed by the flow and represents the state of dynamics system, such as the storage of reservoir system. Connecter represents the relationship among variables, and is denoted as line with an arrow. Obviously, it is appropriate to carry out reservoir SD simulation using the Vensim simulation environment. A causal loop diagram of reservoir system is shown in Fig. 3 (a) that represents the interrelationship between various variables and helps user understand the dynamic behaviour of reservoir system. In Fig. 3 (a), the lines with an arrow are the connecter, and indicate clearly the cause relationships among inflow, storage, head and so on. The sign (+ or -) at the end of the lines suggests 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, spill. The sign * indicates that the correlativity is not clear and the feedback interlinkage may present different polarities at different time periods.

Self-optimization system dynamics modelling
Compared with the fitting method and the conventional simulation-based optimization method, the SD simulation is more flexible and more transparent inherently, and is helpful for increasing participation of model users in the model development, and is able to represent clearly complex causal relationship among system elements using the causal loop diagram. In addition, causal tracing can indicate which type of behaviour is caused by which variable in the reservoir system, and help user figure out which elements of reservoir system influence the power generation or guarantee rate that are the performance index to be refined by SD simulation, and effectively indicate the direction of improvement of IOR. In other words, by adding feedback factors, the power generation or the guarantee rate, the IOR can be improved easily based on feedback loops in SD simulation model. Therefore, a selfoptimization SD simulation model of reservoir operating rules for improving the IOR derived by fitting method is proposed and constructed, and it can directly optimize maximum power generation and maximum guarantee rate by adjusting operating rules parameters iteratively based on feedback information in SD simulation.
Initial operating rules (IOR) Self-optimization system dynamics simulation of operating rules Analyze the feedback information from feedback factor Identify the improvement direction of operating rules Adjust the operating rules parameters   The stock and flow diagram of self-optimization SD model of a reservoir system is shown in Fig. 3 (b). In Fig.  3 (b), the red lines represent the added feedback loops, and begin with the feedback factors, power generation, guarantee rate, reservoir level and discharge. The iterative framework of self-optimizing SD simulation of operating rules is shown as Fig. 2. 3 Case study

The Three Gorges Reservoir
The three gorges reservoir (TGR) is the largest and the most crucial comprehensive reservoir in the Changjiang River, the longest river of China. The average annual runoff at dam site is 4.51×10 11 m 3 /year. The reservoir receives inflow from a catchment area of 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, power 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 streamflow records (from 1882 to 2011) from the Yichang hydrologic station, which is located approximately 40 km downstream of the TGR, were used for the case study. All streamflow records have been restored, therefor, the streamflow statistical feature is consistent. For this case study, the time period of deterministic optimization operation model could be selected as a daily time step, 10 days' time step, a weekly 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 other, considering the daily time step increases the computational complexity and the monthly time step ignores the runoff variation within a month.
Other parameters of the operation model are set as follows. The initial and final 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 m3/s to guarantee downstream navigation and ecological needs. The comprehensive efficiency coefficient K in Eq.
The other data sets of the TGR including reservoir forebay elevation-storage curves, reservoir outflowtailwater elevation curves and head-maximum output curves are used for this case study. Based on the above parameters and data, the self-optimizing SD simulation was applied to the TGR. And some formulations of variables in Fig. 3 are formulated as Table 1.

Optimal operation solution and IOR
The DP algorithm is used to solve the optimal trajectory of the deterministic optimization scheduling model. As shown in Table 2, average annual power generation of 921.6212×108 kWh and average guarantee rate of 99.08% are produced. According to the Table 1. Some variables formulation used in the self-optimization system dynamics simulaiton.  optimal scheduling solution, the scatter plot of the optimal water discharge vs. available water is developed (Fig. 4) and the functional form of operating rules is discussed and analysed. In Fig. 4, an almost linear relationship between water discharge and available water is observed within the feasible spaces that are denoted by the four areas enclosed by lines with different colour. Because the time period of the optimal operation model may be 8, 9, 10 or 11 days, there are four different feasible region corresponding with the four time period. Obviously, the optimal discharge can be fitted by linear operating rules for the flood season, and by piecewise linear operating rules for the non-flood season. The initial linear operating rules (IOR) are derived using linear fitting method based on the above optimal operation solutions, and reused for the reservoir operation. The operation results for the IOR are illustrated in Table  2. Fig. 5 shows the detailed linear operating rules of the last 11 days of August and the last 10 days of April. As shown in Fig. 5(a), during the last 11 days of August, the drived linear operating rules indicates that the decision variable, water discharge, is equal to reservoir inflow. In fact, during flood season, the operating rules functions are the lower reservoir bound in Fig. 4. During the last 10 days of April, a threshold value is defined in Fig. 5 (b). When the available water is less than threshold value, the water discharge is determined by rule 1 in Fig. 5 (b) to ensure to generate guaranteed output, but often failed (look at the blue bar in Fig. 7). When the available water exceeds the threshold, the water discharge is determined according to the rule 2 in Fig. 5 (b) so as to more power generation is generated. However, within an ISO framework, the operating rules derived by fitting method with the maximum goodness-of-fit criterion may be unreliable to guide reservoir operation because of the existence of outliers partly. Therefore, a self-optimization SD simulation (Fig. 3) that is established above is applied to refine the IOR, and the results are illuminated in next section. The last 10 days of April, in no n-flood season.

Optimal operating rules (OOR)
The IOR deduced by linear fitting method in the beginning is used as inputs for the self-optimization SD simulation. And the OOR are obtained by refining the IOR using SD simulation, according to the iterative step in Fig. 2. Then, the average annual power generation and guarantee rate for the IOR, OOR and optimal solution are listed and compared in Table 2. Fig. 6 shows the OOR and IOR for the middle 10 days of April. Further, the three solution all are reused for the reservoir operation with the inflow of 1998 year as the input and Fig. 7 presents the results, power output and water discharge, from the first 10 days of January to the middle 10 days of April. Compared with the IOR, during the middle 10 days of April, according to the rule1, the OOR increase the water discharge to ensure guaranteed output is generated, as shown in Fig. 6. And the observation can be illustrated by the area enclose by ellipse in Fig. 7. In fact, based on Fig.7, the inference that OOR will increase water discharge in order to improve guarantee rate in  non-flood season can be made. Therefore, the OOR is more reliable than the IOR. According to Table 2, the optimal solution produce more power generation and better guarantee rate than the IOR and OOR; the guarantee rate of OOR is better than IOR, while there is a little difference between power generation of OOR and IOR. Obviously, the self-optimization SD simulation is appropriate to refine the IOR, and the OOR is more effective and more reliable than IOR. In addition, the SD simulation is helpful for increasing participation of model users in the simulation of operating rules and increase users confident in operating rules.

Discussion of results
The self-optimization SD simulation of reservoir operating rules for the TGR expands the application of SD approach and broadens the simulation-based optimization method that has been widely used to derive the operating rules.
From the research above, the fitting method with the maximum goodness-of-fit criterion may be unreliable to derive operating rules, which may lead to lower guarantee rate, but a new simulation-based method, the selfoptimation SD simulation of operating rules, is verified to be an effective and reliable way to refine the initial operating rules (IOR) that are generally obtained by fitting method in the beginning. Moreover, the self-optimization SD simulation model is able to increase participation of operators in the model development and help operators understand causal realtionship among variables of reservoir system, and promote the OOR that ensure both power generation and guarantee rate to be produced effectively and reliably. Finally, IOR are improved using the self-optimization SD simulation, and the OOR are obtained. Compared with the IOR, the OOR will increase water discharge duing the non-flood season to increase guarantee rate. Compared with the increased guarantee rate, the reduction in power generation is small.

Conclusions
This study focuses on realizing the self-optimization SD simulation of operating rules so as to improve initial operating rules derive by fitting method. The deterministic reservoir optimization model with two-objective function is constructed. And the DP algorithm is used to solve the optimal model, based on 130-year inflow records (from 1882 to 2011). The linear fitting method is employed to fit the IOR from the optimal solution obtained by DP. Finally, the self-optimization SD simulation model is established to refine the IOR and obtain the OOR by adding feedback loops that begin with the feedback factors, guarantee rate, discharge, power generation and water level and end with the operating rules. Overall, from this study, conclusions can be drawn a follows: (1) the operating rules derived by fitting method with the maximum goodness-of-fit criterion may be unreliable to guide reservoir operation because of the existence of outliers partly. (2) with guarantee rate and power generation as feedback factors, the selfoptimization SD simulaion can effectively imporve the IOR and obtain the OOR. (3) compared with the IOR, the guarantee rate of the OOR is increased by 3.2% , while the average annual power generation is reduced by only 0.03%. The self-optimization SD simulation model has been successfully applied to the TGR, The OOR are obtained and can be used to guide reservoir operation. Nevertheless, the application of self-optimization SD simulation in cascade reservoir operating rules requires further researth.