A novel nested dynamic programming (nDP) algorithm for multipurpose reservoir optimization

Historically, the two most widely used methods for optimal reservoir operation have been dynamic programming (DP) and stochastic dynamic programming (SDP). These two methods suffer from the so-called ‘dual curse’ which prevents them to be used in reasonably complex water systems. The first one is the ‘curse of dimensionality’ that denotes an exponential growth of the computational complexity with the state–decision space dimension (Bellman 1957). The second one is the ‘curse of modelling’ that requires an explicit model of each component of the water system (Bertsekas & Tsitsiklis 1995) to anticipate the effect of each system's transition. The literature offers various strategies to overcome the curse of dimensionality such as successive approximations, incremental DP, and differential DP (Jacobson & Mayne 1970; Bellman 1986; Anvari et al. 2013; Li et al. 2013). The application of various DP and SDP methods in optimal reservoir operation are reviewed in Yeh (1985) and for multi-reservoir systems in Labadie (2004). One of the additional issues arises is when assessment of the constraints and objective functions requires running complex models (for example, river, a model to calculate water levels or flows downstream). One of the possible solutions is the use of fast approximate surrogate models (e.g., artificial neural networks replicating the behaviour of a hydrodynamic river model (Solomatine & Torres 1996)) – this reduces the total running time but still for real-life problems the use of DP can be prohibitively long.

This paper addresses the problem of optimal reservoir operation concerning multiple objectives that are related to: (1) reservoir releases to satisfy several downstream users competing for water with dynamically varying demands; (2) deviations from the target minimum and maximum reservoir water levels; and (3) hydropower production that is a combination of the reservoir water level and the reservoir releases.

Addressing such a problem with classical DP requires a reasonably high level of discretization of the reservoir storage volume, which in combination with the required releases discretization for meeting the demands of downstream users leads to computationally expensive formulations and causes the curse of dimensionality.

We present an alternative approach, in which at each transition of the classical DP an additional optimization algorithm is run to identify the optimal releases (allocations) to individual users. Because this second optimization algorithm is ‘nested’ inside the classical DP algorithm we name this method ‘nested dynamic programming (nDP)’. Depending on the way we formulate the objective function related to deficit in the allocation problem in the nested optimization, two methods are implemented: (1) Simplex for linear allocation problems and (2) quadratic Knapsack method in the case of non-linear problems.

The presented nDP algorithm was coded as a prototype application in Java using Eclipse. The DP part of the nDP was coded by using the example from the book by Loucks & Van Beek (2005). The nested optimization algorithms, Simplex and quadratic Knapsack methods are developed as separate modules. The nDP algorithm was tested at the multipurpose reservoir Knezevo of the Zletovica hydro-system located in the Republic of Macedonia, with the purpose of urban water supply, agriculture, ensuring ecological flow and generation of hydropower.

The article is organized in six sections. After the Introduction, the second section presents the problem, formulation of variables, objective functions and the nDP algorithm. The third section introduces the case study, which is followed by a section presenting the implementation of the nDP for this case. Results and discussion follow, then a section including the comparisons of the nDP with others DP methods, which is followed by conclusions.

where s_t is the state vector representing discrete reservoir storage volume at the beginning of the period t; T is the number of stages in the sequential decision process; is the state value function; is the actions or decision variables vector during the period t; is the reward from period t when the current state is s_t, the action a_t is executed and the resulting state is s_t+1.

The nDP contains a nested optimization algorithm inside the DP algorithm that optimally allocates the total reservoir release r_t to different users corresponding to their demands d_it, as shown in Figure 1.

Assuming a backward-moving DP, at each nDP transition, the reservoir beginning and end states are known. From the mass balance Equation (2), the release r_t can be calculated and then nested optimization algorithm is run to identify allocation of r_t between n users. The inputs of the nested optimal allocation algorithm are the reservoir release, the users' demands and their relative importance, while the decision variables are the water volumes given to users (called subsequently the ‘users' releases') for satisfying their demands. The overall objective function may include terms related to users' releases, deviations from target reservoir levels, or both releases and levels (e.g., hydropower). At the end of the transition, all decision variables are estimated and the reward is calculated.

The nDP pseudo-code is presented below:

1. Discretize storage s_t and s_t+1 in m intervals, i.e., s_i,t (i = 1, 2, …, m), s_j,t+1(j = 1, 2, …, m).

2. Set time at t = T − 1.

3. Set reservoir level i = 1(for time step t).

4. Set reservoir level j = 1(for time step t + 1).

5. Calculate the first group of the objective functions (related to deviations from target reservoir levels).

6. Calculate the total release r_t using Equation (2).

7. If total release r_t is bigger or equal to total demand then return 0, and go to step 9.

8. Execute the nested optimization algorithm to allocate the total release to all users {r_1t, r_2t, … r_nt} in order to meet their individual demands.

9. Calculate the second group of the objective functions (related to users' releases).

10. Using the reservoir levels and the user releases, calculate the third group of the objective functions (related to hydropower production).

11. Combine the objective functions from steps 5, 8 and 9 into the aggregate objective function V(s_t).

12. j = j + 1.

13. If j ≤ m, go to step 5.

14. Select the optimal actions (decision variables) {a_1t, a_2t, … a_nt}_opt, which consist of the optimal transition {s_t+1}_opt and the users' releases {r_1t, r_2t … r_nt}_opt that give minimal value of V(s_t).

15. i = i + 1.

16. If i ≤ m, go to step 4.

17. If t = 0, stop.

18. t = t − 1.

19. Go to step 3.

The optimal allocation algorithm is incorporated (nested) in the DP method (step 7) and directly updates the state value function V(s_t) at each time step consequently changing the optimal reservoir policy and solving the multi-objective optimization problem.

The action vector {a_1t, a_2t, … a_nt} consists of the transition state s_t+1, and the users' releases {r_1t, r_2t, … r_nt}. The corresponding total release r_t, at each transition is calculated from Equation (3). The e_t evaporation losses can be calculated at each transition using the reservoir area that is a function of the reservoir storage volume s_t, and other factors.

Depending on the problem formulation various methods can be used to optimally allocate the total reservoir release r_t between n water users. Optimality could be understood in different ways and this is reflected in the objective function. We have chosen to use two methods: Simplex method in the case of linear problem and quadratic Knapsack algorithm for quadratic problem. Each water user is described with its demand d_it and corresponding weight W_it at time step t. For the nested optimal allocation the following variables are relevant: d_1t, d_2t, … d_nt are users' demands; W_1t, W_2t, … W_nt are the corresponding demands' weights; r_t is the reservoir release; r_1t, r_2t, … r_nt are the users' releases, v is the release discretization value.

The hydro system Zletovica is located in the eastern part of the Republic of Macedonia and it uses the water resources of the river Zletovica and its tributaries. It is composed of the reservoir Knezevo, several water distribution canals used for delivering water to downstream users and associated infrastructure structures, as shown in the schematic representation in Figure 2. The reservoir Knezevo is a multipurpose reservoir, and its main objective is to provide drinking water to several towns and populated areas in the region, as well as to provide environmental flows in Zletovica, water for agriculture and hydropower (in this order of decreasing priority).

There are five towns in this region (Kratovo, Probishtip, Zletovo, Shtip and Sveti Nikole) and two large agricultural irrigation regions named upper and lower zone. The town Kratovo is not modelled because the intake for this town is upstream from the Knezevo reservoir and only a small water quantity is required. The region is characterized with mountainous topography and the system is designed to include several small hydropower plants, most of which are located downstream of the reservoir Knezevo as derivational (run of river) power plants that utilize the natural head differences created by the topography. The hydropower system is still in development and according to the feasibility study report (GIM 2010) the plan is to build eight hydropower plants.

The GIM (2008) report contains a detailed hydrological model of the Zletovo river basin with four river flow measurement points on the river Zletovica. The monthly time series data of the flow measurement points from the year 1951–1990 are used in this research. There is a significant tributary inflow between the Knezevo reservoir and the last branching point. First, the tributary inflow is used to satisfy the water demand objectives, and if additional water quantities are needed then they are released from the reservoir.

The action vector a_t consist of six actions or decision variables: s_t+1, r_3t, r_4t, r_5t, r_6t, r_7t which are the next optimal reservoir state and water user releases at each time step. Using these decision variables, it is possible to calculate all other variables and objective functions.

Since in Equation (9) various objective function types are aggregated into one, their values need to be adjusted so they are comparable in order of magnitude. The approach used was to assign the appropriate coefficients (c_i in Equation (9)). A similar approach to balancing was taken in other previous research studies, e.g., Pianosi & Soncini-Sessa (2009) and Rieker (2010).

The nested algorithm is selected in the beginning of optimization. The nDP is constructed in a way that first the tributary flow (downstream of the Knezevo reservoir) is optimally divided between users' demands using the previously selected nested optimization algorithm. Afterwards, the remaining, still unsatisfied user demands are met from the reservoir releases by applying the optimization as described above.

The nDP algorithm was tested on the Zletovica hydro system using 40-year monthly data from 1951 to 1990, which resulted in 480 time steps. The reservoir operation volume is 23 × 10⁶ m³ which was discretized in 115 equal levels (200 × 10³ m³ each). The target for the recreation reservoir level was set at 1,035 [amsl], for the flood control reservoir level at 1,058 [amsl] and for the hydropower production at 1.5 GW-h per month. The targets for the consumptive users were set to their predetermined monthly demands.

To test and evaluate the proposed nDP algorithm and to compare results using the two nested optimization methods, Simplex and quadratic Knapsack, four different scenarios were considered, as shown in Table 1. They are based on two sets of weights for objective function and two formulations for the deficits in the nested optimization problem: (a) linear (hence the Simplex method, Equation (4)) named L₁ and L₂; and (b) quadratic (the quadratic Knapsack method, Equation (8)) named Q₁ and Q₂. In the quadratic Knapsack the discretization of allocation was set to 50. The higher weights are set for urban water supply objectives (w₃ and w₅) followed by ecological flow, irrigation, reservoir water level objectives and hydropower. The second set of weights gives higher priority to one of the urban water supply objectives (O₅) at the expense of the two irrigation objectives (O₄ and O₆).

Depending on the formulation that is used in the nested optimization, the obtained differences in results can be observed in Figures 3,⁴–5. The quadratic formulation in Q₁ and Q₂ scenarios influence the D₁ and D₂ deviations, producing better results compared with the L₁ and L₂ scenarios, respectively, as can be seen in Figure 3. The D₃ and D₅ deficits in the L₁ and L₂ scenarios are at the expense of D₄ and D₆ (agriculture) that have higher deficits, as shown in Figure 4. When comparing L₁ and L₂ with Q₁ and Q₂ scenarios, it can be seen that D₃ and D₅ have less deficits when a linear formulation is used. The quadratic formulation (quadratic Knapsack) makes a more balanced distribution between the water demand users, which is due to its objective function, shown in Equation (8). The quadratic formulation objective function alleviates the D₆ peak deficits that occur in L₁ and L₂ scenarios, at the expense of D₃, D₄ and D₅. The quadratic formulation has better results concerning D₈ hydropower deficits, as shown in Figure 5.

These differences in results when using linear or quadratic formulations in the nested algorithm come from the properties of the objective functions as described in Equations (4) and (8). Linear objective functions have the known characteristic that optimal solutions are always corner-point feasible solutions, which means that when the objective function is formulated as a weighted sum of separate linear objectives, the optimal solution will favour objectives (‘users' in our case) with higher weights. Therefore, D₅ and D₃ (which are with the highest weights) show smaller deficits for L₁ and L₂ than with Q₁ and Q₂. The quadratic formulation allows for a more balanced solution overall (with better trade-offs) and in fact alleviates huge deficits obtained with linear formulation for users with lower weights. However, linear formulation performs faster, and depending on the problem, targets and the ‘desired’ optimal reservoir operation, the modeller can choose between the linear and quadratic formulation. It is also important to note that the proposed solution is general in the sense that it allows for inclusion of other nested algorithms, and it is not restricted to the two proposed and tested here.

The w₅ (corresponding to the urban water supply) is the highest weight in the L₂ and the Q₂ scenarios shown in Table 1, and consequently, it is the most satisfied water demand user, as shown in Figure 4. The D₃ deficit in the L₁ scenario is zero, but it has a significant increase due to the w₅ increase in the L₂ scenario, although w₃ weight is unchanged, as shown in Figure 4. The D₄ and D₆ deficits are slightly increased in L₂, compared to the L₁ scenario. The D₃ and D₅ deficits are lower in Q₂ compared to the Q₁ scenario. These deficit reductions are achieved at the expense of D₄ and D₆ deficits. When comparing L₁, L₂, Q₁ and Q₂ scenarios considering D₁, D₂, D₃, D₇ and D₈ deficits and their weights, it is obvious that the results are different even without changes of the corresponding weights, as shown in Figures 4,⁵–6. This happens because the change in one or more weights, in this case w₄, w₅ and w₆, modifies the main objective function and, consequently, produces different results. The presented scenarios demonstrate that by changing the weights in accordance with the user preferences it is possible to create different optimal reservoir operations.

As it is difficult to present the 40 years nDP optimization results' graphs, a 3-year (1985–1987) sample was selected from the Q₂ scenario to show some characteristic results. Figure 6 shows the variations of reservoir inflow, tributaries inflow, reservoir volume and reservoir outflow. Figure 7 shows the recreation and flood target’ deviations. Figure 8 shows water users' demands. Figure 9 shows water users' deficits and Figure 10 shows hydropower deficits. This 3-year period is a combination of two dry and one wet year in between them.

This 3-year period is characterized by the significant reservoir and tributary inflows in the spring months (from February until May) due to high precipitation and snowmelt, which is relatively small in the other periods, as shown in Figure 6. The reservoir releases are rising between June and October because of the increased demand of the agricultural users, as shown in Figures 6 and 8. The town users (d_3t and d_5t) and the ecology (d_7t) demands are relatively constant, while the agriculture users' (d_4t and d_6t) demands are variable reaching a maximum in the summer months, as shown in Figure 8. The reservoir inflow and the tributary inflow in 1986 are larger and more widely distributed between January and June compared to 1985 and 1987 when there are high reservoir inflows in April and May, as shown in Figure 6. The difference between 1986 and the other 2 years, 1985 and 1987, is shown in the reservoir release shown in Figure 6, the water users' deficits in Figure 9 and hydropower deficits shown in Figure 10. The majority of users' deficits occur from April until October in 1985 and 1987, as shown in Figure 9. The deficits D₄ and D₆ (agriculture) are larger than D₃ and D₅ (towns) deficits in 1985 and 1987, because of the smaller weights w₄ and w₆ and their larger demand, as shown in Figures 8 and 9. If the linear formulation results were presented using the L₁ or the L₂ scenarios, the difference between the agriculture and the town's deficits would be even larger.

The reservoir satisfies (almost) all objectives in 1986, as shown in Figures 9 and 10. There are reservoir flood deviations (D₂) in May 1985 and 1987 because of peak reservoir inflows, as shown in Figures 6 and 7. On the other side, there are recreational level deviations (D₁) in October and December 1985. This coincides with the reservoir volume shown in Figure 6, which is lowest in periods between November and January. The last months of 1987 suggest that probably the next year is very dry. The reservoir stores additional volume in winter months (October–December 1987) at the expense of water demand users and the hydropower, as shown in Figures 6, 9 and 10.

The nDP algorithm was tested with a variable reservoir volume discretization using the Q₂ scenario settings. The reservoir storage volume was discretized: (1) from 0 to 5 × 10⁶ m³ and from 22 × 10⁶ to 23 × 10⁶ m³ in intervals of 200 × 10³ m³; and (2) the rest was discretized in intervals of 150 × 10³ m³ (resulting in the 144 discretization levels). The higher discretization has directly increased the computational resources and the running time, but it has improved the results (1–2%). This test shows that the nDP algorithm allows for higher levels of reservoir level discretization without prohibitive computational time, and for each problem, a balance between the desired level of accuracy and the computational time and resources can be determined. In addition, the nDP algorithm was tested with variable objective weights at different time steps.

In designing the optimal reservoir operation, one needs to be aware of the objectives' interdependence, in our case, users' releases, reservoir level and hydropower. The hydropower production is a function of the total users' releases (reservoir release) and the reservoir level. The hydropower weight influences both the total users' releases and the reservoir level. An additional scenario was executed based on the Q₂ scenario settings, where the highest priority (weight) was set to the hydropower. The scenario result showed that the reservoir was filling to the highest possible level and afterwards made significant releases, increasing the hydropower production. The increased hydropower production directly influenced the users' releases and the reservoir levels. This shows that although the main objective function is composed of a sum of weighted objectives, there is additional complexity because of the objectives' interdependencies, which means that changing the objective weight influences other dependent objectives and the overall results. The desired optimal reservoir operation can be designed by tuning the weights.

The nDP makes an exhaustive search over all possible states and actions, executing the nested optimization at each transition, and making a significant number of calculations. For the presented set-up, the scenario execution time was under 5 minutes on a standard desktop computer. As expected, the execution time of the Q₁ and the Q₂ scenario was longer because of the more complex calculation in the quadratic Knapsack than in the Simplex optimization algorithm. The minor obstacle for the nDP implementation was the memory requirements that in our case were solved by increasing the maximum memory allocation pool for the Java Virtual Machine.

The presented case study has five water demand objectives. The classical DP approach would model these objectives as five different releases. Since the DP algorithm is an exhaustive search, all possible states and actions need to be evaluated. Let us assume that the discretization is the same and that the five releases are discretized on 10 × 10³ m³.

The DP algorithm makes a transition between all states. Let us consider only one transition from full to empty reservoir. In this simple calculation, the possible combinations of five releases are below 2,300⁵. There are 480 time steps and over 3,000 transition at each time step. A simple assessment shows that the total number of DP transitions is around 10²³, which makes this problem computationally unsolvable. For the applications considered, the usage of classical DP becomes impossible when there are several decision variables.

To decrease the optimization problem dimension and to allow for DP implementation a simplified sub-optimal approach can be taken: the water demand objectives can be combined into a single aggregated objective so that at the first stage the standard DP optimization can be employed. The resulting releases are quasi-optimal since only the aggregated demand is optimized rather than the individual ones. At the second stage, these quasi-optimal reservoir releases are distributed between the different water demand users for each time step according to their importance. Such an approach can be called the aggregated water demand (AWD) DP algorithm. It would be good to compare AWD DP with nDP but using the case considered in this paper it appeared to be difficult.

The new optimization problem now has four objectives, (1) recreation level target, (2) flood level target, (3) (aggregated) water demand user and (4) hydropower, and it can be solved by AWD DP. In the first stage, the DP algorithm is executed with the four objectives and the total optimal reservoir releases identified. In the second stage, these releases are be distributed at each time step between the five users, forming allocations ⁠, ⁠, ⁠, and ⁠.

Let us assume that hydropower is a priority, and only the agriculture-related allocations and are used for producing hydropower. In this case, the AWD approach cannot be applied because it is difficult to create an objective function in DP that will represent the hydropower weight when all water demand users are aggregated (since hydropower is generated only by two out of four allocations). The solution for this problem would be to divide the water demand users and their allocations into two groups, (a) ⁠, and (b) ⁠, ⁠, ⁠, and in this way to represent the hydropower in the objective function. This will introduce the second decision variable in DP that considerably increases the dimensionality of the problem. If we make this case more complicated and assume that the water demand user is very important compared to the other users, and needs to be represented as a separate decision variable in DP, then the application of the classical AWD would become even more computationally demanding.

An alternative is to use the nested approach, as introduced in this paper, which makes it easy to introduce a large number of separate objective functions (e.g., hydropower or d_3t prioritization as discussed above).

This article presented a novel nested approach to solving a DP problem that allows for several objective functions at each time step. The use of the nDP algorithm was implemented in the Zletovica river basin case study with eight objectives and six decision variables. The nDP algorithm has the following advantages:

• It effectively alleviates the curse of dimensionality of a standard DP algorithm.

• It allows for considering all objective functions at each time step and hence allows for deeper optimization if compared to the sub-optimal approach where all water demands are aggregated (AWD DP).

• Computationally, nDP runs fast on standard personal computers. In the presented case study, optimization for one scenario required no more than 5 minutes on a standard desktop PC.

• The algorithm allows for employing dense and variable discretization of the reservoir volume and release.

• It supports using a variable (time-dependent) weight at each time step for every objective function.

• It provides a framework for including many objectives in the nested optimization algorithm without a significant change in the source code or an increase in the computational expenses.

• Different optimization algorithms can be used in the nesting for water allocation; however, since nested optimization has to be repeated multiple times (for each transition of DP) the algorithm used for this purpose needs to be fast. Thus, in the presented nDP framework we used Simplex and quadratic Knapsack algorithms which showed excellent performance.

• The nDP framework can support complex non-linear objective functions.

• The presented nesting approach can be applied not only in the deterministic DP context, but also in SDP, reinforcement learning and other optimization algorithms.

One of the limitations of the presented study is the fully deterministic formulation, without explicit consideration of uncertainties, so characteristic of the reservoir operation. Therefore, further research will be focussed on implementing the nDP approach in a SDP setting. In addition, it is planned to test the nested optimization approach in the problems solved by reinforcement learning.