A New Reservoir Operation Chart Drawing Method Based on Dynamic Programming

: A reservoir operation chart is an important tool in guiding actual reservoir operation at present. There are mainly two kinds of methods in drawing the operation chart, i.e., conventional methods and optimization methods, but each of them has some shortcomings, such as the repeated empirical inspection and correction of conventional methods, and the sensitivity to the initial trajectories of some optimization algorithms. In view of this, based on the principle of dynamic programming (DP), this paper coupled the reservoir operation chart drawing model and the DP model, and proposed a new reservoir operation chart drawing method which has no empirical inspection and correction, no requirement for initial solution, no problem of premature convergence and local convergence. In addition, this method can guarantee the global convergence of the results to a certain extent because of the global convergence of DP. Ya Yangshan reservoir in the Li Xianjiang River of China was selected as the research object to derive the operation chart by the drawing method. The simulation results show that the proposed method in this paper presents better performance compared with the conventional method on power generation, guaranteed output, and assurance rate, especially on the latter, which has a 2.68% increase. In addition, compared with the deterministic optimization results of DP, it is found that the results of the proposed method are very close to that of deterministic DP, the differences are only 1.8 GWh (0.36% decline) and 1.6 GWh (0.32% decline). So, the validity and rationality of the proposed method are further veriﬁed by the simulation results.


Introduction
As one of the clean energy sources, hydropower energy is the most important renewable source for generation of electricity worldwide and can be developed on a large-scale [1,2]. It can reduce the use of fossil fuels for electricity generation, and reduce worldwide CO 2 -emmissions. Hydropower contributed 16.5% to the world electricity generation in 2012, while the other renewables contributed only 5.2% [3]. Compared with other energy sources, hydropower enjoys exceptional advantages [4,5] which are clean, pollution-free, quick in output, and can quickly adapt to load changes of power system [6]. In addition, water is the main input for hydropower station to produce the electricity [7], and it is present and usable all over the year in contrast to wind and solar that are intermittent technologies and only usable when these resources are available, so hydropower energy is a superior sustainable energy source to help maintain sustainable growth and quality of life [8]. These advantages of hydropower mean that reservoir operation optimization research has the attention of many scholars all over the world [9][10][11]. However, in the actual long-and mid-term reservoir operation, operation Dynamic programming (DP) is a powerful optimization technique [32,33], and the most striking feature of DP is that it can ensure the global optimum and no initial solution is needed. Besides, DP has no requirements on non-convex and unsmooth nature of the optimization problems. This makes DP obtain a high popularity in the classical optimization algorithms [34]. However, although many methods were used to optimize the reservoir operation chart as mentioned above, DP was rarely used to derive the reservoir operation chart directly. Therefore, the optimality of DP has not been well applied to the reservoir operation chart optimization.
In view of this, the reservoir operation chart drawing model and the DP model are coupled in this research, and this paper proposed a new reservoir operation chart drawing method based on DP. Compared with the conventional drawing method, it is more direct and effective to directly calculate the optimal operation chart through the DP model, and there is no repeated empirical inspection and correction. Compared with other optimization algorithms, there is no requirement for an initial solution, and there is no problem with premature convergence and local convergence too. In addition, the proposed method can guarantee the global convergence of the results to a certain extent because of the global convergence of DP. The remaining parts are organized as follows in this paper. Section 2 will present the drawing model of reservoir operation chart based on DP, including the specific calculation process and the detailed steps in the actual application. Section 3 will show the case study of Ya Yangshan reservoir in Li Xiangjiang River, and the results will be presented and analyzed in Section 3.2. Finally, the conclusions of this paper will be summarized and provided in Section 4.

DP Based Reservoir Operation Chart Drawing Model
There are three types of output in the actual hydropower station dispatching, i.e., guaranteed output, increased output, and reduced output. Correspondingly, there are three kinds of operation zones in the operation chart, which are the guaranteed output zone, the increased output zone and the reduced output zone. So, there are three types of operation curves in a reservoir operation chart, i.e., the basic operation curves (including the upper and lower basic operation curve), the increased output curves, and the reduced output curves. The upper basic operation curve and lower basic operation curve corresponds to the upper and lower boundary of guaranteed output zone.
The conventional reservoir operation chart is shown in Figure 1, in which the number of the increased and the reduced output curves can be more than one according to the needs. Here, for the simplified representation, only one increased output curve and one reduced output curve are provided. In Figure 1, the abscissa axis represents the time (month), and the ordinate axis represents the reservoir water level whose upper limit is the flood control level in the flood season the normal storage level in the non-flood season [35]. The curves in Figure 1 from top to bottom are in turn the increased output curve, the upper basic curve, the lower basic curve and the reduced output curve. These curves divide the whole dispatching map area into several corresponding output zones. The detailed process of how to use the DP model to directly calculate the optimal reservoir operation chart will be described in the following.

3
In view of this, the reservoir operation chart drawing model and the DP model are coupled in this research, and this paper proposed a new reservoir operation chart drawing method based on DP. Compared with the conventional drawing method, it is more direct and effective to directly calculate the optimal operation chart through the DP model, and there is no repeated empirical inspection and correction. Compared with other optimization algorithms, there is no requirement for an initial solution, and there is no problem with premature convergence and local convergence too. In addition, the proposed method can guarantee the global convergence of the results to a certain extent because of the global convergence of DP. The remaining parts are organized as follows in this paper. Section 2 will present the drawing model of reservoir operation chart based on DP, including the specific calculation process and the detailed steps in the actual application. Section 3 will show the case study of Ya Yangshan reservoir in Li Xiangjiang River, and the results will be presented and analyzed in Section 3.2. Finally, the conclusions of this paper will be summarized and provided in Section 4.

DP Based Reservoir Operation Chart Drawing Model
There are three types of output in the actual hydropower station dispatching, i.e., guaranteed output, increased output, and reduced output. Correspondingly, there are three kinds of operation zones in the operation chart, which are the guaranteed output zone, the increased output zone and the reduced output zone. So, there are three types of operation curves in a reservoir operation chart, i.e., the basic operation curves (including the upper and lower basic operation curve), the increased output curves, and the reduced output curves. The upper basic operation curve and lower basic operation curve corresponds to the upper and lower boundary of guaranteed output zone.
The conventional reservoir operation chart is shown in Figure 1, in which the number of the increased and the reduced output curves can be more than one according to the needs. Here, for the simplified representation, only one increased output curve and one reduced output curve are provided. In Figure 1, the abscissa axis represents the time (month), and the ordinate axis represents the reservoir water level whose upper limit is the flood control level in the flood season the normal storage level in the non-flood season [35]. The curves in Figure 1 from top to bottom are in turn the increased output curve, the upper basic curve, the lower basic curve and the reduced output curve. These curves divide the whole dispatching map area into several corresponding output zones. The detailed process of how to use the DP model to directly calculate the optimal reservoir operation chart will be described in the following.  Firstly, separate the four curves in Figure 1 and place them in four separated coordinates, and get the feasible water level range of each operation curve, and discretize it, as shown in Figure 2.  Figure 1 is taken out separately and placed in a separate coordinate system. Figure 2b represents that the upper basic operation curve in Figure 1 is taken out separately and placed in a separate coordinate system. Figure  Firstly, separate the four curves in Figure 1 and place them in four separated coordinates, and get the feasible water level range of each operation curve, and discretize it, as shown in Figure 2. Figure 2a represents that the increased output curve in Figure 1 is taken out separately and placed in a separate coordinate system. Figure 2b represents that the upper basic operation curve in Figure 1 is taken out separately and placed in a separate coordinate system. Figure 2c,d have similar meanings to Figure 2a   Secondly, from Figure 2, it can be seen that, corresponding to each discretized water level point of the basic operation curves, the increased output curve and the reduced output curve at the beginning of a stage, a water level combination can be obtained, and another similar combination by the discretized water level points of each curve at the end of this stage can be obtained too. The two combinations can constitute a stage-operation-chart if their position relationship is satisfied (no cross), as shown in Figure 3. Secondly, from Figure 2, it can be seen that, corresponding to each discretized water level point of the basic operation curves, the increased output curve and the reduced output curve at the beginning of a stage, a water level combination can be obtained, and another similar combination by the discretized water level points of each curve at the end of this stage can be obtained too. The two combinations can constitute a stage-operation-chart if their position relationship is satisfied (no cross), as shown in Figure 3. All the possible discretized water level combinations of the four operation curves over the whole operation period are shown in Figure 4, where C(k, t) represents the kth combination in the tth stage, k = 1, 2, …, M n − 1, M n ; n is the count of the operation curves in the operation chart; M is the number of discrete points of the feasible water level range; t is the index of stage, and t = 1, 2, …, T; T is the  All the possible discretized water level combinations of the four operation curves over the whole operation period are shown in Figure 4, where C(k, t) represents the kth combination in the tth stage, k = 1, 2, . . . , M n − 1, M n ; n is the count of the operation curves in the operation chart; M is the number of discrete points of the feasible water level range; t is the index of stage, and t = 1, 2, . . . , T; T is the total number of stages over the entire planning horizon.  All the possible discretized water level combinations of the four operation curves over the whole operation period are shown in Figure 4, where C(k, t) represents the kth combination in the tth stage, k = 1, 2, …, M n − 1, M n ; n is the count of the operation curves in the operation chart; M is the number of discrete points of the feasible water level range; t is the index of stage, and t = 1, 2, …, T; T is the total number of stages over the entire planning horizon.

M 4 combination
Stage  On this basis, the reverse recursion calculation for obtaining the reservoir operation chart can be carried out directly according to the DP model. Two recursion procedures are used to calculate the reservoir operation chart in DP, which are the reverse recursion procedure and the chronological order recursion procedure. Starting from the Tth stage, the output or power generation can be calculated up to the first stage in the reverse recursion procedure, and then, through the chronological order recursion procedure, the optimal discretized water level combination variation of each stage On this basis, the reverse recursion calculation for obtaining the reservoir operation chart can be carried out directly according to the DP model. Two recursion procedures are used to calculate the reservoir operation chart in DP, which are the reverse recursion procedure and the chronological order recursion procedure. Starting from the Tth stage, the output or power generation can be calculated up to the first stage in the reverse recursion procedure, and then, through the chronological order recursion procedure, the optimal discretized water level combination variation of each stage can be gained. For the operation chart calculation with one year runoff data, the recursive equation for the tth stage can be expressed as [36]

C(1,T)
where OCB t (C k1 t−1 ) is the optimal cumulative benefit of the beginning combination C k1 t−1 at the tth stage; OCB t+1 (C k2 t ) is the optimal cumulative benefit of the beginning combination C k2 t at the (t + 1)th stage; N t () is the output function of an operation stage and it is related to C k1 t−1 and C k2 t ; k1 represents the index of a discretized water level combination at the beginning of a stage; k2 represents the index of a discretized water level combination at the end of a stage; k1 = 1, 2, . . . , M n − 1, M n ; k2 = 1, 2, . . . , M n − 1, M n ; C represents a discretized water level combination; C t−1 k1 is equal to the C(k1, t − 1) in Figure 4, and C t k2 is equal to the C(k2, t) in Figure 4. The optimal cumulative benefit here means the sum of output or power generation in the optimal operation process from stage t to stage T.
In calculating the operation chart by DP with the Y years runoff data, for each discretized water level combination in the reverse recursion process, the output calculation of runoff data of all the Y years in the current stage needs to be repeatedly implemented to obtain the optimal cumulative benefit (average annual power generation) of current stage. Then, take the average annual beginning water level of Y years as the final beginning water level of current stage. So, the recursive equation for the tth stage calculation can be represented as In addition, in the reverse recursion process, there is a formula that reflects the beginning water level relationship of each stage, which can be shown by the Formula (3) This formula means that the beginning water level Z t−1 that corresponding to the discrete water level combination C t−1 k1 in the tth stage is the function of the end water level Z t and the average annual power generation ∑Nt −1 (C t−1 k1 , C t−1 optk2 )/Y which is corresponding to the optimal candidate path, where the optk2 is determined by formula (2), and F represents a function. Before the reverse recursion calculation, the data known are the long series of runoff data Q j t (j = 1, 2, . . . , Y; t = 1, 2, . . . , T) and the fixed water level at the end of the whole operation period. The fixed water level is usually the dead level. Then, the specific steps of calculating reservoir operation chart based on DP can be summarized as follows.
Step 1: Divide the whole operation period into T stages, and get the discrete water level combinations based on the feasible range (Z lower~Zupper ) of each operation curve in each stage. There will be M discrete points for each operation curve in each stage, and M n discrete water level combinations if the count of the operation curves in the operation chart is n.
Step 2: Start reverse recursion procedure from stage T. For any stage-operation-chart in the current stage, for example, the stage-operation-chart constituted by the combination C(3, T − 1) and C(1, T) in Figure 4, the optimal cumulative benefit OCB(3, T − 1), the optimal candidate path OCP(3, T − 1) and the corresponding beginning water level Z(3, T − 1) for combination C(3, T − 1) can be figured out by an assumption calculation. The known inflow Q T j (j = 1, 2, . . . , Y) and the fixed end water level (dead level) of combination C(1, T) are used in this calculation. Save the obtained OCB(3, T − 1), OCP(3, T − 1), and the corresponding Z(3, T − 1) for combination C(3, T − 1), and go to Step 3.
Step 4: Set t = t − 1, go to the next stage's calculation. Unlike the optimal cumulative benefit OCB(k, T − 1) in Step 2, which only represents the benefit of current stage, the OCB(k, t) is the maximal sum of the benefits from present stage t to last stage T, and it can be calculated by Formula (2). In addition, the end water level of each combination in the tth stage is no longer the dead level, but the beginning water level Z(k, t + 1) that determined in the (t + 1)th stage's calculation.
Step 5: The reverse recursion procedure is over when the first stage's calculation has been done. Then, based on the saved optimal candidate paths, trace back the optimal path from the first stage to the Tth stage to obtain the optimal combination trajectories {C t }, and obtain the reservoir operation chart further by {C t }.
The whole flowchart of calculating the reservoir operation chart based on DP is provided in Figure 5. Step 5: The reverse recursion procedure is over when the first stage's calculation has been done. Then, based on the saved optimal candidate paths, trace back the optimal path from the first stage to the Tth stage to obtain the optimal combination trajectories {Ct}, and obtain the reservoir operation chart further by {Ct}.
The whole flowchart of calculating the reservoir operation chart based on DP is provided in Figure 5. No t=t-1

Yes
Calculate annual average N t and water level Z t by inflow Q t (j=1,2,…,Y) and the stageoperation-chart built by C t-1 k1 and C t k2 Calculate ,and the corresponding optimal k2 from the range of 1, 2,…, M n .
Carry out the chronological order recursion according to the optimal cumulative benefit array OCB , and get the optimal combination of each stage, i.e., {C t , t=1,2,…,T} Output the whole operation chart Input data, including the number of stages and states, and the boundary conditions Build the discrete combination space according to Fig.4.
Save the optimal k2 and the corresponding and Z t    In the output calculation of each stage-operation-chart, there are two aspects that require attention. Firstly, there is an assumption calculation in the derivation of the beginning water level of each combination, that is, assume the unknown beginning water level is located in one of the five operation zones (i.e., zone A, B, C, D, and E) in Figure 3. If it is zone A, the output can be determined as N A , then, the beginning water level Z(k, t) can be figure out by an iterative calculation using the known inflow Q T j (j = 1, 2, . . . , Y), N A and the end water level Z(k, t + 1). If the obtained beginning water level Z(k, t) is located in zone A, the hypothesis is right. Otherwise, suppose another output situation and continue, until all the possible output situations are traversed (as shown in Figure 6). At this moment, there may be more than one output situations that meet the hypothesis. So, there is a comparison among these output conditions, and finally keep one output situation and the corresponding beginning water level according to the principle of maximizing the cumulative benefit. While, if there is no output situation that meets the hypothesis, the assumed output that has the minimum error can be chose as the final output, but this situation should be punished at a certain extent, so that to avoid the combinations which have this situation being selected at last. The calculation process mentioned above is shown in the left part of Figure 5. a comparison among these output conditions, and finally keep one output situation and the corresponding beginning water level according to the principle of maximizing the cumulative benefit. While, if there is no output situation that meets the hypothesis, the assumed output that has the minimum error can be chose as the final output, but this situation should be punished at a certain extent, so that to avoid the combinations which have this situation being selected at last. The calculation process mentioned above is shown in the left part of Figure 5. Secondly, In the process of calculating operation chart by DP, with the goal of maximizing the power generation, the assurance rate constraint also need to be taken into consideration. Therefore, in the output calculation of each stage, the situation in which the output is lower than the guaranteed output should be punished, so that to avoid the combinations which have this kind of situation being selected at last, and try to control the output damage depth when it is inevitable.
In the output calculation, the following constraints need to be satisfied [37].
(1) Water volume balance: (2) Volume limits of reservoir: (3) Discharge flow requirements of downstream reservoirs: Secondly, In the process of calculating operation chart by DP, with the goal of maximizing the power generation, the assurance rate constraint also need to be taken into consideration. Therefore, in the output calculation of each stage, the situation in which the output is lower than the guaranteed output should be punished, so that to avoid the combinations which have this kind of situation being selected at last, and try to control the output damage depth when it is inevitable.
In the output calculation, the following constraints need to be satisfied [37].
(1) Water volume balance: (2) Volume limits of reservoir: (3) Discharge flow requirements of downstream reservoirs: (4) Power generation limits: Constraint of the expected output is considered in the power generation limits, and it can be expressed by Formula (8) where H t is the average water head in the tth stage, unit: m; I t is the inflow rate in the tth stage, unit: m 3 /s; k is the efficiency coefficient; L t is the discarded water outflow rate in the tth stage, unit: m 3 /s; N t is the output in the tth stage, unit: kW; N t,min is the lower output limit in the tth stage, unit: kW; N t,max is the upper output limit in the tth stage, unit: kW; N t,exp is the expected output in the tth stage, unit: kW; q t is the outflow rate through the turbines of hydropower station in the tth stage, unit: m 3 /s; Q t is the whole outflow rate of reservoir in the tth stage, and Q t = q t + L t , unit: m 3 /s; Q t,min is the lower limit of Q t in the tth stage, unit: m 3 /s [38]; Q t,max is the upper limit of Q t in the tth stage, unit: m 3 /s; V t−1 is the volume of the reservoir at the beginning of the tth stage, unit: m 3 ; V t,min is the lower Energies 2018, 11, 3355 9 of 17 limit of V t in the tth stage, unit: m 3 ; V t,max is the upper limit of V t in the tth stage, unit: m 3 ; ∆t is the duration of an operation stage, which is 2,628,000 s in this paper. After the reservoir operation chart is obtained, it can be used to guide the actual reservoir operation, that is, the operation decision can be made according to the reservoir's current water level and the operation zone where the water level located in the chart. There are several kinds of output situations. If the water level located in the increased output zone, the decision is increased output. If the water level located in the guaranteed output zone, the decision is guaranteed output. If the water level located in the reduced output zone, the decision is reduced output. The purpose of drawing reservoir operation chart is to guide the actual reservoir operation, but its rationality needs to be tested by simulation operation using a long series of historical runoff data. The steps of the simulation can be summarized as follows.
Step 1: According to the reservoir's current water level, determine the output N t of current stage.
Step 2: According to the reservoir's initial state of current stage and the inflow rate, carry out the output calculation by the determined output N t , and obtain the end state of the reservoir, i.e., the water level at the end of this stage.
Step 3: Determine whether the obtained end state of current stage satisfy the constraint or not, that is to check whether the water level at the end of this stage is between dead level and normal level (or flood control level in flood season) or not. If the constraint is satisfied, operate the hydropower station by N t in this stage, otherwise, go to Step 4.
Step 4: If the water level is higher than normal level (or flood control level in flood season) at the end of this stage, then set the water level equal to normal level (or flood control level), and calculate the actual output of this stage [39]. If water level is lower than the dead level at the end of this stage, then set the water level equal to the dead level and calculate the actual output of this stage.
Step 5: Finally, get the average annual power generation and assurance rate by statistics, and analyze the variations of output process and water level process, so that to verify the rationality of the obtained reservoir operation chart.

Data and Brief Introduction to Research Object
Li Xianjiang River is a tributary of Red River, and it is located in Yunnan province of China. There is a 473 km main channel, 1790 m natural drop, and 19,309 km 2 control area in China. At the border of China, the annual average runoff is 460 m 3 /s. There are eight hydropower stations on the main channel, i.e., Tu Kahe, Ge Lantan, Xin Anzhai, Ju Pudu, Long Ma, Xin Pingzhai, Shi Menkan, and Ya Yangshan from downstream to upstream. Three hydropower stations in these hydropower stations have the seasonal regulation performance, which are Long Ma, Shi Menkan, and Ya Yangshan. The Ya Yangshan is the topmost and the dominant hydropower station of the main channel, whose operation level can determine the whole power generation benefit of the cascade reservoirs. Therefore, Ya Yangshan hydropower station is selected as the research object in this paper. Location of the selected hydropower station in Li Xianjiang River basin is shown in Figure 7 [40], and its basic parameters are shown in Table 2. The available long series of runoff data for this hydropower station is from 1957 to 2000, a total of 43 years.

10
Yangshan. The Ya Yangshan is the topmost and the dominant hydropower station of the main channel, whose operation level can determine the whole power generation benefit of the cascade reservoirs. Therefore, Ya Yangshan hydropower station is selected as the research object in this paper. Location of the selected hydropower station in Li Xianjiang River basin is shown in Figure 7 [40], and its basic parameters are shown in Table 2. The available long series of runoff data for this hydropower station is from 1957 to 2000, a total of 43 years.

Results of the Two Methods
As can be seen from Table 2, the guaranteed output of Ya Yangsha reservoir is 23.2 MW, and the assurance rate is 95%. However, in the actual calculation by conventional method [41] or DP based optimization method (as shown in Figure 5), it was found that the requirements of the guaranteed output and assurance rate cannot be satisfied at the same time. As we know, the assurance rate and power generation of a hydropower station are usually gradually reduced with the increase of the guaranteed output, so it is an effective way to reduce the guaranteed output to improve the assurance rate and power generation. There is an inverse relationship between the assurance rate and guaranteed output in the hydropower station operation. Therefore, for facilitating the comparative analysis, this paper firstly provides the results of the conventional method, in which only one indicator (the guaranteed output or assurance rate) is satisfied. Then, the results of DP based method will be provided in the same form.
For the conventional method, when the guaranteed output of the reservoir is 23.2 MW, the simulation results of the conventional reservoir operation chart through a long series of historical runoff can be obtained. The assurance rate is 93% and the power generation is 496.3 GWh. The corresponding reservoir operation chart is shown in Figure 8. Where, the 1.0 × 23.2 MW corresponds to the basic operation curve, the 1.1 × 23.2 MW corresponds to the increased output curve whose output is 1.1 times the guaranteed output, the 1.2 × 23.2 MW corresponds to the increased output curve whose output is 1.2 times the guaranteed output, and the 0.9 × 23.2 MW corresponds to the reduced output curve whose output is 0.9 times the guaranteed output.

Results of the Two Methods
As can be seen from Table 2, the guaranteed output of Ya Yangsha reservoir is 23.2 MW, and the assurance rate is 95%. However, in the actual calculation by conventional method [41] or DP based optimization method (as shown in Figure 5), it was found that the requirements of the guaranteed output and assurance rate cannot be satisfied at the same time. As we know, the assurance rate and power generation of a hydropower station are usually gradually reduced with the increase of the guaranteed output, so it is an effective way to reduce the guaranteed output to improve the assurance rate and power generation. There is an inverse relationship between the assurance rate and guaranteed output in the hydropower station operation. Therefore, for facilitating the comparative analysis, this paper firstly provides the results of the conventional method, in which only one indicator (the guaranteed output or assurance rate) is satisfied. Then, the results of DP based method will be provided in the same form.
For the conventional method, when the guaranteed output of the reservoir is 23.2 MW, the simulation results of the conventional reservoir operation chart through a long series of historical runoff can be obtained. The assurance rate is 93% and the power generation is 496.3 GWh. The corresponding reservoir operation chart is shown in Figure 8. Where, the 1.0 × 23.2 MW corresponds to the basic operation curve, the 1.1 × 23.2 MW corresponds to the increased output curve whose output is 1.1 times the guaranteed output, the 1.2 × 23.2 MW corresponds to the increased output curve whose output is 1.2 times the guaranteed output, and the 0.9 × 23.2 MW corresponds to the reduced output curve whose output is 0.9 times the guaranteed output. When the assurance rate of the reservoir is asked to not less than 95% in the conventional method, the simulation results of the conventional method show that, the guaranteed output is 22.4 MW, and power generation is 496.7 GWh. The corresponding reservoir operation chart is shown in Figure 9.  When the assurance rate of the reservoir is asked to not less than 95% in the conventional method, the simulation results of the conventional method show that, the guaranteed output is 22.4 MW, and power generation is 496.7 GWh. The corresponding reservoir operation chart is shown in Figure 9.
For the DP based method, when the guaranteed output of the reservoir is 23.2 MW, the reservoir operation chart obtained by the proposed method in this paper is shown in Figure 10. Its simulation results through a long series of historical runoff show that, the corresponding assurance rate is 94%, and the power generation is 498.9 GWh.
12 Figure 9. Conventional reservoir operation chart considering the assurance rate constraint.
For the DP based method, when the guaranteed output of the reservoir is 23.2 MW, the reservoir operation chart obtained by the proposed method in this paper is shown in Figure 10. Its simulation results through a long series of historical runoff show that, the corresponding assurance rate is 94%, and the power generation is 498.9 GWh. When the assurance rate of the reservoir is asked to not less than 95% in the DP based method, the reservoir operation chart obtained by the proposed method is shown in Figure 11. Its simulation results show that the corresponding guaranteed output is 23.0 MW, and the power generation is 501.2 GWh.   For the DP based method, when the guaranteed output of the reservoir is 23.2 MW, the reservoir operation chart obtained by the proposed method in this paper is shown in Figure 10. Its simulation results through a long series of historical runoff show that, the corresponding assurance rate is 94%, and the power generation is 498.9 GWh. When the assurance rate of the reservoir is asked to not less than 95% in the DP based method, the reservoir operation chart obtained by the proposed method is shown in Figure 11. Its simulation results show that the corresponding guaranteed output is 23.0 MW, and the power generation is 501.2 GWh.  When the assurance rate of the reservoir is asked to not less than 95% in the DP based method, the reservoir operation chart obtained by the proposed method is shown in Figure 11. Its simulation results show that the corresponding guaranteed output is 23.0 MW, and the power generation is 501.2 GWh. Figure 11. DP based reservoir operation chart considering the assurance rate constraint.
Whether the actual dispatcher takes the assurance rate or the guaranteed output as the main control factor to decide the final reservoir operation chart, the actual operation requirements and situations of hydropower station are needed. For the hydropower station operation in Li Xianjiang River, considering the demands of power grid security and stability, there is more focus on the reliability of hydropower generation. Thus, this paper takes the assurance rate as the main control factor to determine the final reservoir operation chart. In view of this, the second case mentioned above conforms more to the actual situation. Whether the actual dispatcher takes the assurance rate or the guaranteed output as the main control factor to decide the final reservoir operation chart, the actual operation requirements and situations of hydropower station are needed. For the hydropower station operation in Li Xianjiang River, considering the demands of power grid security and stability, there is more focus on the reliability of hydropower generation. Thus, this paper takes the assurance rate as the main control factor to determine the final reservoir operation chart. In view of this, the second case mentioned above conforms more to the actual situation.

Contrastive Analysis
In order to facilitate the comparative analysis, the results of the two methods under the same boundary conditions are put together. For example, when the guaranteed output of the conventional method and the proposed method are both the same 23.2 MW, the simulation results of each method are shown in Table 3. When the assurance rate of the two methods are both asked to not less than 95%, the simulation results of each method are shown in Table 4.  Table 3, it can be seen that, compared with the conventional method, the assurance rate and power generation of the proposed method in this paper both have a certain degree of improvement when the guaranteed output of the two methods are the same 23.2 MW, and the growth is 0.52% and 1.08%, respectively. From Table 4, it can be seen that, when the assurance rate of the two methods are both asked to not less than 95%, the power generation and guaranteed output of the proposed method also have a certain degree of improvement compared with the conventional method, and the growth is 0.91% and 2.68%, respectively. So, from the aspects of assurance rate, guaranteed output, and power generation, all the simulation results indicate that the proposed method based on DP in this paper is better than the conventional method.
Besides, in comparing the two situations-i.e., the situation of the guaranteed output of the two methods are the same 23.2 MW and the situation of the assurance rate of the two methods are both asked to not less than 95%-it is found that the latter situation has a better result on the annual power generation, that is, 0.91% growth is greater than 0.52% growth. This result is consistent with that this paper takes the assurance rate as the main control factor to determine the final reservoir operation chart.
Compared the Figure 10 with Figure 8, it can be found that the guaranteed output zone of the Figure 10 has relatively expanded a little, and the reduced output zone has narrowed down. This change brings in a corresponding improvement for the assurance rate, which is consistent with the simulation results in Table 3 (increased from 93% to 94%).
Compared the Figure 11 with Figure 9, it can be seen that the guaranteed output zone of the Figure 11 has relatively expanded, and the reduced output zone has narrowed down. Generally, this change will make the assurance rate increase if the guaranteed output unchanged. However, the guaranteed output increased, which is changed to 23.0 MW from 22.4 MW at this moment, so the assurance rate remains unchanged.
Similarly, compared the Figure 11 with Figure 10, it can be found that the guaranteed output zone of the Figure 11 has relatively expanded, and the reduced output zone has narrowed down, which will make the assurance rate increase if the guaranteed output fixed. However, at this moment, the change of the operation zones cannot make the assurance rate reach 95%, therefore, the guaranteed output needs to be reduced (changed from 23.2 to 23.0), so that to make the assurance rate reach 95%.
In addition, in order to compare the results with that of deterministic DP, this paper implements the direct optimization by DP when the assurance rate is asked to not less than 94% and 95% respectively, and the guaranteed output is fixed as 23.2 MW and 23.0 MW correspondingly. The optimization results are as follows. When the assurance rate is 94% and the guaranteed output is 23.2 MW, the annual average power generation of DP is 500.7 GWh, and when the assurance rate is 95% and the guaranteed output is 23.0 MW, it is 502.8 GWh.
For the proposed method in this paper, the power generation is 498.9 GWh when the assurance rate is 94% (as shown in Table 3), and the power generation is 501.2 GWh when the assurance rate is 95% (as shown in Table 4), which are both slightly less than the results of DP. However, the difference is small, which is only 1.8 GWh (0.36%) and 1.6 GWh (0.32%), respectively. This result indicates that the proposed method in this paper has maintained the global convergence of DP well to a certain extent.
Taking the proposed method as an example, the average annual water level variation of Ya Yangshan reservoir is shown in Figure 12 when the assurance rate is 95%, and the average annual output variation is shown in Figure 13.
14 change brings in a corresponding improvement for the assurance rate, which is consistent with the simulation results in Table 3 (increased from 93% to 94%).
Compared the Figure 11 with Figure 9, it can be seen that the guaranteed output zone of the Figure 11 has relatively expanded, and the reduced output zone has narrowed down. Generally, this change will make the assurance rate increase if the guaranteed output unchanged. However, the guaranteed output increased, which is changed to 23.0 MW from 22.4 MW at this moment, so the assurance rate remains unchanged.
Similarly, compared the Figure 11 with Figure 10, it can be found that the guaranteed output zone of the Figure 11 has relatively expanded, and the reduced output zone has narrowed down, which will make the assurance rate increase if the guaranteed output fixed. However, at this moment, the change of the operation zones cannot make the assurance rate reach 95%, therefore, the guaranteed output needs to be reduced (changed from 23.2 to 23.0), so that to make the assurance rate reach 95%.
In addition, in order to compare the results with that of deterministic DP, this paper implements the direct optimization by DP when the assurance rate is asked to not less than 94% and 95% respectively, and the guaranteed output is fixed as 23.2 MW and 23.0 MW correspondingly. The optimization results are as follows. When the assurance rate is 94% and the guaranteed output is 23.2 MW, the annual average power generation of DP is 500.7 GWh, and when the assurance rate is 95% and the guaranteed output is 23.0 MW, it is 502.8 GWh.
For the proposed method in this paper, the power generation is 498.9 GWh when the assurance rate is 94% (as shown in Table 3), and the power generation is 501.2 GWh when the assurance rate is 95% (as shown in Table 4), which are both slightly less than the results of DP. However, the difference is small, which is only 1.8 GWh (0.36%) and 1.6 GWh (0.32%), respectively. This result indicates that the proposed method in this paper has maintained the global convergence of DP well to a certain extent.
Taking the proposed method as an example, the average annual water level variation of Ya Yangshan reservoir is shown in Figure 12 when the assurance rate is 95%, and the average annual output variation is shown in Figure 13.  . Figure 13. Average annual output variation of Ya Yangshan reservoir.
From Figure 12, it can be seen that the reservoir filling up quickly at the beginning stages of the storage period, and then maintains the operation with a relatively high water level, so that to increase the water head efficiency. In the supply period, the reservoir's water level drops homogeneous, which makes the output uniformed, so that to avoid the concentrated destruction of output or water abandonment at the last few stages of the entire planning horizon. Correspondingly, from Figure 13, From Figure 12, it can be seen that the reservoir filling up quickly at the beginning stages of the storage period, and then maintains the operation with a relatively high water level, so that to increase the water head efficiency. In the supply period, the reservoir's water level drops homogeneous, which makes the output uniformed, so that to avoid the concentrated destruction of output or water abandonment at the last few stages of the entire planning horizon. Correspondingly, from Figure 13, it can be seen that the reservoir basically maintains the guaranteed output in the supply period, except for the increased output of the last stage because of emptying the reservoir. Thus, it can be concluded that the average annual operation process is conformed to the actual operation principle, which reflects the consistency of actual facts with the simulation results of the obtained reservoir operation chart. Hence, the rationality of the proposed method in this paper is further proved.

Summary and Conclusions
Reservoir operation chart plays a very important role in the actual operation of reservoir, and this paper presents a new method to calculate the reservoir operation chart based on DP. Compared with the conventional drawing method, it is more direct and effective to calculate the optimal operation chart through the DP model directly, and to a certain extent, it can guarantee the global convergence of the results because of the global convergence of DP. This paper took the Ya Yangshan reservoir of Li Xianjiang River in southwest China as an example to derive the operation chart by the proposed method, and the following conclusions can be summarized: (1) Through the simulation of case study, the results show that the proposed method presents better performance compared with the conventional method, 2.68% growth on assurance rate and 0.91% growth on power generation can be obtained by the proposed method. Besides, the simulation operation processes (water level and output) of the proposed method are conformed to the actual operation principle, this reflects the consistency of actual facts with the simulation results of the obtained reservoir operation chart. So, the validity and rationality of the proposed method are verified by the simulation results.
(2) Compared with the direct optimization results of DP, although the power generations of the proposed method under two assurance rates are both slightly less than that of DP, the differences are small, which is only 0.36% and 0.32% decline, respectively. This result shows that the proposed method in this paper has maintained the global convergence of DP to a certain extent.
(3) However, the advantage of the proposed method is not very significant compare to the conventional method, and the increase on power generation is only 0.91%. The possible reason may be the simplified treatments in the reverse recursion procedure, which cannot give full play to the global convergence of DP.
In summary, the study in this paper has achieved some good results. The proposed method can avoid the shortage of other methods in drawing and optimizing the reservoir operation chart, and make up for the deficiency of the research on DP based reservoir operation chart optimization. However, some work still needs to be done in the future. For improving the significance of this method, fine processing is needed in the reverse recursion procedure, and for the validity of the proposed method for other hydropower stations, more case studies need to be done, especially the joint operation chart derivation of cascade reservoirs. In addition, in the derivation of joint operation chart of cascade reservoirs by DP, there may be a long computing time, so the parallel computing techniques or other measures need to be adopted.