A Decomposition and Dynamic Programming Aggregation Method for the Optimal Water Allocation of Reservoirs in Series

Research on water allocation of multiple reservoirs with the purpose of reducing water spills and improving the local runoff utilization is a matter of great concern in humid areas with uneven temporal and spatial distributions of water resources. An optimization model for a system of reservoirs in series is developed to minimize water shortages. Several constraints restrict the objective function, including available water, operation rules and water rights for replenishment of the reservoirs with water. The model features multiple dimensions with a single coupling constraint of the large-scale system. A decomposition and dynamic programming aggregation method (DDPA) is proposed; the subsystem models and the aggregation model are both solved with the classical one-dimensional dynamic programming. Compared with the conventional decomposition-coordination method, the proposed method is concise but reliable because it can directly use the results of subsystems to form the one-dimensional dynamic programming aggregation model, avoiding the iterative calculations according to the coordinating function. Compared with the meta-heuristic algorithms, the proposed method is more efficient because it is independent of any algorithm parameter. The proposed method may provide a new reference for solving similar multi-reservoir optimization models.


Introduction 1
Reservoirs are the main water sources in the hilly regions of southern China and 2 southeast Asia. The mean annual rainfall is over 1,000 mm in these humid regions, but 3 70-80% of that occurs during the flood season due to the monsoon climate (Quinn et al. 4 2018). Although the annual total inflows into the local reservoirs are usually more than 5 what is needed to meet demands, a large amount of water is drained due to uneven 6 temporal and spatial distribution of inflows, which results in seasonal water shortages. 7 Because of this shortage, it is necessary to replenish a reservoir from other reservoirs 8 during the dry season (Ming et al. 2017; Reca et al. 2015), thus gradually forming a system 9 of reservoirs in series. 10 The basic principle of joint operation for a multi reservoir system is to redistribute 11 water resources through the hydraulic connections among connected reservoirs (Ahmad 12 et al. 2014). This setting can maximize the storage capacity of each reservoir and reduce 13 loss of water. Thus, the operation purpose of reservoirs in series is to reduce water spills 14 and improve runoff utilization; the decision variable is the amount of water supply in each 15 period, which is subject to the reservoir's capacities including water storage and water 16 supply (Rani et  be added into constraints, which inevitably complicate the model (Ehteram et al. 2017;20 Ashrafi and Dariane 2017). 21 water allocation models of reservoirs. Dynamic programming based on Bellman's 23 principle (Bellman and Dreyfus 1964) is highly applicable to this type of multi-stage 24 decision-making process. However, the 'curse of dimensionality' is induced when dealing 25 with a large number of reservoirs (Cheng et al. 2017; Chen et al. 2016). Decomposition 26 (Turgeon 1981) is the mainstream idea to dealing with multidimensional optimization 27 problems when using dynamic programming. The decomposition-coordination method 28 the optimal state of the whole system through the successive iterative calculations between 30 the large-scale system and subsystems according to the coordinating variables. Therefore, 31 the convergence performance of the coordination variable has a great impact on the final 32 results. In addition, the Lagrange multiplier method (Jafari and Alipoor 2011) is usually 33 used to deal with the constraints when adopting the conventional decomposition-34 coordination method. However, it is restricted by the differentiability and convexity of the 35 objective function and it may have difficulty in handling some complex constraints, for 36 instance, the constraint which contains if statements. 37 Recently, meta-heuristic algorithms such as the genetic algorithm (Allawi et al. 38 2018a), the particle swarm optimization algorithm (Chen et al. 2018) and the ant colony 39 algorithm (Moeini et al. 2013) have become the most popular methods for the 40 optimization of reservoir operation models. The greatest advantage of meta-heuristic 41 algorithms is that they can solve the multidimensional model without the decomposition 42 operations. However, the globally optimal final results cannot be guaranteed because 43 decision variables are randomly sampled and updated within the feasible regions according to the specific iteration rules, which is a common limitation among these 45 algorithms (Hossain and El-shafie 2013; Allawi et al. 2018b). Furthermore, when the 46 meta-heuristic algorithms are adopted, a constrained optimization problem usually 47 their own supply objects. During the operation period, each reservoir in the system can be 67 replenished with water from another reservoir or from an outside river by pumping or 68 gravity. 69 70 Note: (a) a system of reservoirs in series; (b) a subsystem of a reservoir; (c) a subsystem 71 of a reservoir and a pumping station 72

Figure 1 A system of reservoirs in series 73
In Figure 1: and EFi,t (L 3 ) are water supply, 74 water demand, local inflow, water replenishment, water spill and evaporation of reservoir 75 i in period t, respectively, N is the total number of reservoirs in the system, i is the sequence 76 number of reservoir, i = 1, 2,…, N, and t is the sequence number of each period. 77 1.2 Optimization model 78

Objective function 79
The operation cycle of the system in this study is one year, so the objective function is to minimize the annual sum of the squared water shortage of each reservoir ( where F is the annual sum of squared water shortage of each reservoir in each period and 84 T is the total number of periods. 85

Constraints 86
① Annual available water of the whole system 87 The annual available water of the whole system includes the total available water of 88 the reservoirs and the available volume from outside water sources limited by the water 89 rights, which can be expressed as Equation (2): 90 where SK is the total annual available water of the reservoirs and BZ is the maximum 92 volume that can be acquired from outside water sources. 93 ② Annual available water of reservoir i 94 The annual available water of each reservoir is the sum of local inflows and water 95 replenishments from other water sources, which can be expressed as Equation (3): 96 The lower and upper bounds of each reservoir's water storage in each period can be 102 expressed as Equation (5): 103 where , According to the water balance equation, the water storage of reservoir i in period t 107 can be determined by Equation (6): 108 The water replenishment and water spill of each period can be determined according to , then excess water should be released. Yi,t and PSi,t can be determined by 116 Equation (9) and Equation (10), respectively. 117 , then both of Yi,t and PSi,t should be zero, shown as Equation (11). 120 , , ④ Maximum water supply can be expressed as Equation (13)  124 , ⑤ Maximum water replenishment is restricted by the pumping capacity or water 126 conveyance capability, which can be expressed as Equation (14): 127 It is also necessary to introduce non-negativity constraints to avoid negative values, 135 which would be physically impossible. 136

Solving methods 137
In the above model, the constraints from to can be transformed into the feasible 138 regions of the decision variable Xi,t. Therefore, it is a dynamic programming problem of 139 N+1 dimensions, which is restricted by the constraints and . A decomposition and 140 dynamic programming aggregation method for this optimal water allocation of reservoirs 141 in series is proposed in this paper, which solves both of the subsystem models and the 142 aggregation model with dynamic programming. In order to evaluate the solving 143 performance, two commonly used meta-heuristic algorithms, the genetic algorithm and particle swarm optimization algorithm, are adopted simultaneously. 145

Decomposition and dynamic programming aggregation 146
① Decomposition of the system 147 The 'curse of dimensionality' will inevitably arise if using dynamic programming 148 (DP) to solve the model directly, hence it is necessary to first reduce the dimensions. The 149 large-scale system will then be decomposed into several subsystems and each subsystem 150 only consists of a reservoir or a reservoir and a pumping station as shown in Figure 1 The optimization model of the subsystem is expressed as Equation (15) and (16),  154 where the objective function is expressed as 155 156 and the annual available water of subsystem is expressed as 157 where fi is the annual sum of squared water shortage in each period of subsystem i. 159 Moreover, the lower and upper bounds of water storage, water balance equation, 160 initial and boundary conditions should also be imposed on each subsystem, which can be 161 considered in the operation rule. The operation rule is integrated into recursive procedures 162 of DP to correct water storage Vi,t and obtain water spill PSi,t and water replenishment Yi,t 163 of each period simultaneously. In this manner, one-dimensional dynamic programming 164 can be adopted to solve the subsystem models as shown in Figure 2(a).

③ Aggregation of the system 166
For subsystem i, the annual available water Wi can be discretized in its feasible region, 167 and the objective function value fi can be obtained under each Wi. After that, the 168 aggregation model can be formed directly using the optimization results of subsystems 169 fi(Wi), which can be expressed as Equation (18) and (19). 170 Objective function: Coupling constraint: The aggregation model is also a typical one-dimensional dynamic programming 175 model, where the decision variable is the annual available water of subsystem Wi . The 176 inverse recursion of DP should be applied as shown in Figure 2

(b) Recursion process of DP for the aggregation model
Test the operation rule 184 Note: d is a certain increment; the aggregation model is optimized in reverse order so that 185 Yi+1 is a known quantity when Vi is calculated according to Equation (6); M is a positive 186 number, which is large enough. 187  The optimizing procedure of PSO starts with a number of initial particles provided 216 stochastically in the feasible region. Each particle's position vector and velocity vector 217 are updated according to Equations (23) and (24) in each iteration until the termination 218 criterion is met: 219 where υs is the velocity vector of particle s, θs is the position vector of particle s, k is the 222 number of iteration times, ω is the inertia weight, ps is the best position of particle s, g is 223 the global best solution among all the particles, c1 and c2 are acceleration coefficients and 224 r1 and r2 are random values that vary in [0, 1]. 225

Case study 226
In order to compare the performances of the three algorithms, the model is applied 227 to a dual-reservoir-and-dual-pumping-station system in Nanjing, China, which consists of 228 the SH reservoir, the HWB reservoir and two pumping stations, HZ and XZ, as shown in 229 The system provides irrigation water and each reservoir provides water 233 independently to its own irrigation area. The average topography of the HWB reservoir is 234 much higher than that of the SH reservoir. During water shortages, HZ pumping station can replenish the HWB reservoir with water from the SH reservoir, while the SH reservoir 236 can be replenished with water from the BL river through the XZ pumping station.  Tables 1 and 2,  241 respectively. 242  Notes: The annual water rights of the XZ pumping station, which means the maximum 246 volume that can be pumped from the BL river, are allocated by the local water authority; 247 the HZ pumping station is an internal pumping station in the system and is not defined or 248 limited by any water rights. 249

② Inflows and water demands 250
The inflow and water demand in each month at 75% probability of exceedance are 251 shown in Table 3. 252

③ Evaporation 254
Evaporation loss is a function of evaporation depth and average free water surface of 255 a reservoir in each period, which can be described as Equation (25). The average free 256 water surface of each reservoir in a specific period is determined by its surface-volume 257 relationship, which was provided by Liuhe Water Authority. 258 where EFt (10 4 m 3 ) is the evaporation loss of a reservoir in period t, Et (mm) is the evaporation depth of E601 evaporator in period t as shown in Table 4, kt is the correction 261 coefficient for period t as shown in Table 4, Vt (10 4 m 3 ) is the average water storage in 262 period t and α and β are the reservoir coefficients. For the SH reservoir, α = 1.194×10 -3 , β 263 = 2.575, and for the HWB reservoir, α = 1.657×10 -3 , β = 1.862. 264

Optimization Results 267
The results obtained with different algorithms are shown in Table 5. The conventional 268 operation of the system was conducted in Excel-2010 according to the standard operation 269 policy (SOP), while the other three optimization methods were programmed in Visual 270 Basic language based on a PC with i3 CPU/3.70 GHz/4 GB ram. 271 Compared with the results using SOP, the water allocation of the system can be 272 optimized by using the optimization methods. By increasing the water diverted from the 273 SH reservoir to the HWB reservoir, the total water spill of the SH reservoir and the total 274 water shortage of the HWB reservoir were reduced. Among the three optimization 275 methods, DDPA minimized the objective function better than the other algorithms. The 276 value of objective function obtained with DDPA was decreased by 7.5% and 5.6% when 277 compared to GA and PSO, respectively. 278 Although the amount of water replenishment from the BL river was same, the total 279 water spill of the system obtained with DDPA was reduced by 2.0% while the total water 280 shortage was reduced by 5.3% and 10.0% when compared to GA and PSO, respectively.
It is indicated the utilization rate of local runoff in the DDPA results is relatively higher 282 on premise of the same external water rights. 283 On the other hand, the operating time of DDPA was 3.85 times and 3.10 times than 284 GA and PSO respectively, though the results obtained with DDPA were better. As shown 285 in Table 5, GA and PSO are more efficient in each operation when the algorithm 286 parameters are determined. However, the operating time is quite longer when considering 287 the tests of the appropriate algorithm parameters of GA and PSO, which is discussed in 288 the next section. 289 where Reli (%) is the volumetric reliability index of reservoir i and Vuli (%) is the 308 vulnerability index of reservoir i. 309 The volumetric reliability and vulnerability indexes of different methods are shown 310 in Table 6. For the SH reservoir, the water demand can be met entirely using the SOP or 311 optimization methods because the available water in the reservoir is much more than the 312 water demand. For the HWB reservoir, the volumetric reliability index of DDPA is 99.51%, 313 which is 1.5%, 0.02% and 0.05% more than that of SOP, GA and PSO, respectively. The 314 vulnerability indexes of DDPA and GA are both 8.20%, which are 20.93% and 0.54% less 315 than that of SOP and PSO, respectively. Thus, the performance of DDPA is relatively better because it can simultaneously achieve a high percentage of volumetric reliability 317 index and a low percentage of vulnerability index. 318  Bellman's principle can guarantee the global optimal results. Therefore, the final optimal 344 results achieved by DDPA are always the same while the actual operation time just varies 345 within a minor range. Correspondingly, if considering the time consumed for testing the algorithm parameters, the operation time of GA or PSO would be much longer than that 347 of DDPA as shown in Table 5. 348

Conclusion 349
This study developed an optimization model for the system of reservoirs in series to 350 minimize water shortage, which is restricted by available water, operation rules and water 351 rights for replenishment with water. A decomposition and dynamic programming 352 aggregation method was proposed to resolve this model. The primary characteristic of the 353 method is that the subsystem models and the aggregation model are both solved with the 354 classical one-dimensional dynamic programming. Therefore, it can reduce the dimension 355 by transforming the N+1 dimensional dynamic programming into the N+1 iterative 356 calculations of one-dimensional dynamic programming. 357 Compared with the conventional decomposition-coordination method, the proposed 358 method is concise but reliable because it can directly use the results of subsystems to form 359 the one-dimensional dynamic programming aggregation model, avoiding the iterative 360 calculations according to the coordinating function. Compared with the meta-heuristic 361 algorithms, the proposed method is more efficient because it is independent of any 362 algorithm parameter. 363

Declarations 364
Acknowledgements: the authors wish to acknowledge the managers of Liuhe Water 365 Authority, Nanjing, China, for their support of providing reservoir, irrigation and 366 hydrology data. 367 Funding: this work was supported by the National Natural Science Foundation of China   Flow chart of DDPA Note: d is a certain increment; the aggregation model is optimized in reverse order so that Yi+1 is a known quantity when Vi is calculated according to Equation (6); M is a positive number, which is large enough.

Figure 3
Location of the system Note: The designations employed and the presentation of the material on this map do not imply the expression of any opinion whatsoever on the part of Research Square concerning the legal status of any country, territory, city or area or of its authorities, or concerning the delimitation of its frontiers or boundaries. This map has been provided by the authors.