Optimal allocation model and method for parallel ‘reservoir and pumping station’ irrigation system under insufficient irrigation conditions

The effective regulation and storage of water resources by reservoirs in arid and semi-arid areas is important for alleviating water resource shortages. In this paper, multiple irrigation reservoirs and pumping stations are evaluated, according to their special hydraulic connections and used to establish the water resources optimal allocation model of the parallel ‘reservoir and pumping station’ irrigation system. The mathematical model takes the maximum total income of the whole irrigation area as the objective function, the water supply and spill of the reservoir and the water replenishment of the pumping station as the decision variables, and the system water supply, agricultural water rights, reservoir operation criteria as the constraints. A new multilevel decomposition aggregation dynamic programming (MDADP) algorithm is proposed to solve the complex nonlinear model and is compared with the real-coded genetic algorithm, particle swarm optimization, cat swarm optimization and whale optimization algorithm. From the analysis of optimality and the applicability of the algorithm, MDADP was found to be more suitable than the above heuristic models for solving problems in the field of the optimal allocation of water resources.


Introduction
As global warming and the increase of human activities, the amount of available water resources in the world is decreasing sharply, while the distribution of water resources was becoming increasingly uneven (Ashofteh et al. 2012;Haddeland et al. 2014). Agriculture accounts for approximately 70% of the world's fresh water resources uses (Li et al. 2018), and the serious shortage of water for irrigation was of great concern.
Common water resource projects include reservoirs, ponds, pumping stations and wells. The joint operation of these projects can effectively ensure irrigation in dry years. Therefore, a lot of research focused on the joint regulation of multi-reservoirs (Ashrafi and Dariane 2017;Wan et al. 2017;Louati et al. 2011) and surface water and groundwater (Wan et al. 2018;Alizadeh et al. 2017;Ketabchi and Ataie-Ashtiani 2015). In addition, considering river agricultural water rights, pumping stations can also replenish water resources outside of the local region to supplement irrigation reservoirs or irrigated areas, but the joint investigation of multiple reservoirs and multiple pumping stations was rare.
Many kinds of constraints, such as the water balance equation, hydraulic linkage and dispatching criterion, should be considered in the joint dispatching model of multi-water source engineering. Therefore, due to these coupling constraints and the size of the system, the galway coupling model was usually applied, which presents a nonlinear coupling problem (Pérez-Díaz et al. 2010;Mo et al. 2013;Afshar 2013). To solve this problem, many methods, including linear programming (LP) (Reis et al. 2006;Azamathulla et al. 2008), nonlinear programming (NLP) (Martin 1983;Lund and Ferreira 1996;Barros et al. 2003), dynamic programming (DP) (Kumar and Baliarsingh 2003;Goor et al. 2011;Li et al. 2011) and progressive optimization methods (POA) (Cheng et al. 2012), have been proposed and discussed in past decades. In addition to the above classical methods of function optimization, various heuristic models were also proposed, such as the genetic algorithm (Cinar et al. 2010;Hincal et al. 2011), ant colony algorithm (Jalali et al. 2007;Deep et al. 2009), simulated annealing algorithm (Georgiou et al. 2006), and particle swarm optimization algorithm Wang et al. 2010). These function optimization methods and intelligent optimization algorithms were widely used in the water resources optimal scheduling field. However, due to the existence of complex constraints in the system, the complexity of the mathematical model required has increased. There are drawbacks when applying the above methods/algorithms to solve the problem. For example, although LP can solve high-dimensional problems in theory, it needs long operation time; NLP cannot deal with non-convexity problems and the convergence is poor (Afshar 2013); DP was easily affected by 'Curse of dimensionality' (Cheng et al. 2012). It was difficult to obtain an initial feasible solution for a complex system, which leads to the local optimal solution of the POA, and the common heuristic methods have difficulty dealing with the equality constraint and the judgment constraint (Cheng et al. 2013), there was no strict mathematical proof of the optimal solution (Samanipour and Jelovica 2020).
The decomposition-aggregation method, also known as the two-level algorithm, has been widely used to solve complex large-scale system problems. Through decomposition, the large-scale system was divided into a series of subsystems, which can be solved by the traditional function optimization method, and then a series of subsystems are aggregated by constructing regression equations, from which the global optimal solution of the large-scale system is obtained. On this basis, Gong and Cheng (2018) and Gong et al. (2019) took the coordination variable in the process of decomposition-aggregation as the decision variable, and aggregated a series of subsystems into a dynamic programming model of large-scale systems according to the recursive principle of DP, thus proposing a large-scale system decomposition aggregation dynamic programming (DADP) method.
In this research, an irrigation system composed of many irrigation reservoirs and pumping stations was studied, from which an optimal allocation model of water resources for a 'reservoir and pumping station' irrigation system was proposed. The model established in this paper has the following characteristics: firstly, when the total amount of agricultural water rights was determined, the allocation of water rights in sub-irrigation areas was taken as a decision-making variable. Then, the optimal allocation model of parallel 'reservoir and pumping station' irrigation systems can be decomposed into N 'single reservoir and single pumping station' irrigation system. Secondly, the irrigation water quantity of each subirrigation area for two kinds of dry crops that are rotated was taken as the decision variable in the water-deficient irrigation area, then the optimal water allocation model of a 'single reservoir and single pumping station' irrigation system can be decomposed into two optimal water allocation models of a single crop irrigation system. Therefore, based on the DADP algorithm, a multilevel decomposition aggregation dynamic programming (MDADP) algorithm was proposed to solve the model of a parallel 'reservoir and pumping station' irrigation system. In addition, two classic heuristic [real-coded genetic algorithm (RGA) and particle swarm optimization (PSO)] and two new heuristic [cat swarm algorithm (CSO) and whale optimization algorithm (WOA)] models were used and compared with the proposed MDADP algorithm, and the optimality and the applicability of these algorithms were verified with an example.

Model construction
The water resource shortage area north of the Huai River in Jiangsu Province, China, were used as the study area; this irrigation area is usually planted with wheat and corn. The two dry crops undergo rotation cultivation (the growth stages do not overlap). The entire growth period of wheat is usually from October to May, while the entire growth period of corn is usually from June to September. Therefore, the combined total duration of the growth period for the two dry crops is consistent with the whole hydrological year of the irrigation reservoir (i.e., from October to September). Therefore, in the process of mathematical model optimization, the growth stage of the dry crops was taken as the optimal operation stage of the irrigation reservoir and pumping station. The parallel 'reservoir and pumping station' system is shown in Fig. 1.
where: i is the number of the reservoir or pumping station; N is the total number of reservoirs or pumping stations; j is the crop type (j = 1,2); k is the crop growth stage; X i,j,k is the actual water supply of i reservoir to the k growth stage of j crops, MCM; Y i,j,k is the water replenishment of i pumping station at the k growth stage of j crops, MCM; LS i,j,k is the inflow of i reservoir at the k growth stage of crop j, MCM; PS i,j,k is the water spill of i reservoir at the k growth stage of crop j, MCM; EF i,j,k is the evaporation capacity of i reservoir at the k growth stage of crop j, MCM. (Note: MCM is million cubic meters).

Objective function
The water supply strategy of irrigation reservoirs is usually optimized as an effective method to ensure the efficient utilization of limited water resources. There is a complex relationship between the reservoir water supply and crop yield. On one hand, it shows the impact of full and insufficient irrigation on the maximum and actual yield of final crops; on the other hand, it shows the effect of water stress caused by water deficit at a certain growth stage on crop transpiration and actual yield. The relationship between water and yield proposed by Jensen (1968) and Doorenbos and Kassam (1979) is whidely used to represent this relationship. Based on the water production function given by Jensen, the objective function of maximizing the total agricultural income in the whole irrigation area was established in the current study. Simultaneously, to facilitate the joint operation calculation of the parallel 'reservoir and pumping station' irrigation system, the ratio of the actual evapotranspiration generation and the potential evapotranspiration generation in the model was converted into the ratio of the reservoir water supply and the crop water demand (Wardlaw and Barnes 1999).
For irrigated areas where two major crops are rotated, the objective function can be expressed as: where: G is the total income of the whole irrigation area, 10 4 rmb; M j is the total number of j crop growth stages; (Y m ) i,j is the maximum yield of j crops when the reservoir is fully supplied with water, kg/hm 2 ; h i,j,k is the sensitivity index of the k growth stage yield of the j crop irrigated by i reservoir under water shortage; YS i,j,k is the maximum water demand of the k growth stage of the j crop irrigated by reservoir i, MCM; A i,j is the planting area of the j crop irrigated by reservoir i, hm 2 ; Pi,j is the unit price of the j crop irrigated by reservoir i, kg/rmb. (Note: rmb is the legal currency of China).

Constraints
(1) System water supply constraints: 1) Total water supply constraints of the system: 2) "Single reservoir-single station" irrigation system annual total water supply constraints: 3) Water supply constraints during the whole growth period of a single crop: where: SK i is the total annual water supply of i reservoir, MCM; BZ i is the total water replenishment of i pumping station, that is, the maximum agricultural water right that the pumping station can draw in a year, MCM; SW is the total 2, ..., N;j = 1, 2) 1 3 199 Page 4 of 16 water supply of the system, MCM; IW i,j is the total irrigation water of the whole growth period of j crop irrigated by i reservoir, MCM.
(2) Pumping station constraints: 1) Total water replenishment of pumping stations: 2) Allowable water replenishment of each pump station: 3) Capacity constraint of pumping station: where, SQ is the total amount of water replenishment by the pumping station (water right); Q i is the design flow of the i pump station, m 3 /s; N i,j,k is the maximum operating time of the i pump station in the k growth stage of the j crop, h.
(3) Constraints of reservoir operation criteria: Among these, V i,j,k was determined by the water balance equation: where, V i,j,k is i reservoir water storage at the end of the k growth stage, which should be the corresponding time between the maximum and minimum water storage. 1) If V i,j,k ≤ V min i,j,k , then the reservoir needed to replenish water through the pumping station, and replenish water to the lower limit of the reservoir water storage, The water replenishment of the pumping station is shown in formula (10).
The water spill of the reservoir is shown in formula (11).
2) If V i,j,k ≥ V max i,j,k , then the pumping station did not need to replenish water, and the water spill of the reservoir is shown in formula (12).
k , the reservoir did not need to replenish water or spill water.
(4) Crop water demand constraints: where: YS min i,j,k is the minimum water demand of the j crop in the k growth stage irrigated by the i reservoir, MCM.
(5) Initial and boundary conditions: where: V 0 is the initial storage capacity of the reservoir; V end is the ending storage capacity after reservoir regulation and storage.

Multilevel decomposition aggregation dynamic programming (MDADP)
The optimal allocation model of water resources for a parallel 'reservoir and pumping station' irrigation system established in this study has the following characteristics. Firstly, to distribute the total amount of water raised by the pumping station (the total water right) to each sub-irrigation area, the water resources optimal allocation model of parallel 'reservoir and pumping station' irrigation system can be decomposed into N 'single reservoir and single pumping station' irrigation system; Secondly, the optimal water allocation model of a 'single reservoir and single pumping station' irrigation system can be decomposed into two optimal water allocation models for dry crop irrigation.
Considering that the optimal allocation model of water resources for the parallel 'reservoir and pumping station' irrigation system is a stage divisible high dimensional dynamic programming model, the conventional function optimization method is subject to problems such as 'Curse of dimensionality' (Cheng et al. 2012), which is difficult to solve. Therefore, the use of DADP to solve the high dynamic programming model is proposed. As shown in Fig. 2, DADP combines the decomposition-aggregation (DA) and DP (Bellman and Dreyfus 1964) approaches: first, the largescale system model is decomposed into a series of subsystem optimization models, which are solved by DP, then, the coordination variables of each subsystem are taken as decision variables. According to the recursive principle of DP, a series of subsystem models can be aggregated into a large-scale system DP model, and the global optimal solution can be obtained by examining the optimal solution of each subsystem from the optimal solution of the DP model of a large-scale system.
Obviously, the DADP algorithm only needs one decomposition and aggregation operation to solve the conventional V 0 =V end large-scale system mathematical model. However, the parallel 'reservoir and pumping station' irrigation system built in the current study has multi-level decomposition, so a MDADP was proposed based on the DADP algorithm, the algorithm flow of which is shown in Fig. 3. The MDADP solution steps are as follows: Step1: first-level decomposition of large systems. By determining the total amount of agricultural water rights in the transit river, and discretizing the allowable allocation of water rights in each sub-irrigation area, the first level decomposition can decompose the water resources optimal allocation model of the parallel 'reservoir and pumping station' irrigation system into N "single reservoir and single pumping station" irrigation systems (level I subsystem model); the model is as follows: Objective function: where: F i is the annual output value of the irrigation area where i reservoir is located.

Constraints:
Step2: second-level decomposition of large systems. In the current study, the effects of water stress on the final yield of dry crops was studied, and the water competition of two types of dry crops under the condition of limited water resources was also considered. The optimal water allocation model of any 'single reservoir and single pumping station' irrigation system was decomposed into two optimal water allocation models of a single dry crop (level II subsystem model), which were as follows: Objective function: where: F i,j is the output value of j crop irrigated by i reservoir.

Constraints:
In addition, the level II subsystem model must satisfy the constraints of the dispatching criterion, pumping station capacity, crop water demand and reservoir boundary conditions.
Step3: Level II subsystem solution. The decomposed subsystems belonged to a one-dimensional DP model, which could be solved by the classical DP method. The state variables λ i,j,k are the total water supply of the reservoir to a single crop during the whole growth period, and were discrete in a fixed step d in [0, IW i,j ]. For each state variable λ i,j,k , the decision variable X i,j,k was discretized in its feasible region [0, λ i,j,k ] by the step length d. The optimal water supply process of a reservoir to each growth stage of dry crops in the subsystem can was obtained by using the sequential recursion and the state transfer equationλ i,j,k-1 = λ i,j,k − X i,j,k (i = 1, 2, …, N; j = 1, 2; k = 1, 2, …, M j ) and the value of the objective function F i,j , according to the optimal water supply process of the reservoir, the optimal spill water process of the reservoir PS i,j,k (i = 1, 2, …, N; j = 1, 2; k = 1, 2, …, M j ) and the optimal water replenishment process of the pumping station Y i,j,k (i = 1, 2, …, N; j = 1, 2; k = 1, 2, …, M j ). Step4: Level II subsystem aggregation. A series of corresponding relationships of IW i,j ~ f i, j (IW i,j ) (i = 1, 2, …, N; j = 1, 2) could be obtained by solving the subsystems of level II, according to the recursive principle of DP, then N level II aggregation models were obtained by aggregating the two level II subsystems which that belonged to the same subarea.
Objective function:

Constraints:
Step5: Level I subsystem aggregation. The water rights BZ i (i = 1, 2, …, N) of each sub-irrigation area were discretized, and the discrete values were substituted into the level II aggregation model, then a series of corresponding relationships of BZ i ~ F i (BZ i ) (i = 1, 2, …, N) could be obtained. Similarly, according to the recursive principle of DP, N level II aggregation models could be re-aggregated into level I aggregation models, resulting in the level I aggregation model becoming a large-scale system model.

Constraints:
By using the classical DP method to solve the level I aggregation model, the optimal annual output value G * of the whole irrigation area and the optimal allocation water right BZ i * (i = 1, 2, …, N) of each sub-irrigation area could be obtained according to the results of BZ i * (i = 1, 2, …, N). From the reexamination of the level II aggregation model, the optimal total irrigation water for different dry crops could be obtained by IW i,j * (i = 1, 2, …, N; j = 1, 2); Finally, according to the optimal total irrigation water quantity of dryland crops, the optimal water supply process X i,j,k * (i = 1, 2, …, N; j = 1, 2; k = 1, 2, …, M j ) and water spill process PS i,j,k * (i = 1, 2, …, N; j = 1, 2; k = 1, 2, …, M j ) of the reservoirs in the whole growth period of the different dryland crops could be obtained, along with the optimal replenishment water process of the pumping station Y i,j,k * ( i = 1, 2, …, N; j = 1, 2; k = 1, 2, …, M j ).

Real-coded genetic algorithm (RGA)
Wright (1991) proposed a RGA, which is a real-to-individual coding method that can achieve better optimization results than binary-coded genetic algorithms.
For a given optimization problem with m variables, the n initial solution sets x 1 i , x 2 i , ..., x m i are usually generated in a random manner when the RGA is used to solve the problem, where the process of population iteration, selection, crossover and mutation are adopted to ensure the evolution of the population. In this study, a roulette wheel was used for the selection operation, and random single-point method was used for cross operation. First, a sequence number t smaller than m was randomly generated, and then the sequence numbers corresponding to a pair of genes t to m were exchanged with each other to complete the crossover operation. For example: a pair of genomes X 1 = x 1 1 , ..., x t 1 , ..., x m 1 and X 2 = x 1 2 , ..., x t 2 , ..., x m 2 will become X 1 � = x 1 1 , ..., x t 2 , ..., x m 2 and X 2 � = x 1 2 , ..., x t 1 , ..., x m 1 through the cross operation; The mutation operation was completed by random single point mutation.
Considering that the mathematical model established in this study contains many constraints, in the process of solving RGA, the penalty term was constructed to deal with many these constraints, as follows: where, P 1 , P 2 , P 3 , and P 4 are the penalty functions of the upper and lower limits of reservoir capacity, the water right constraints of sub-regions and the agricultural total water right constraints of the whole irrigation area respectively, and μ 1 , μ 2 , μ 3 , and μ 4 are the penalty factors respectively, in the actual solution process, the penalty factors are all taken as − 10,000,000.
By integrating the objective function of the model with the penalty term, the fitness function of the RGA was constructed as follows: where: W is the fitness function value.

Particle swarm optimization (PSO)
J. Kennedy first proposed the PSO algorithm, which originated from bird predation behavior research. The algorithm abstracts each bird as a particle, and the location of the food is the optimal solution of the problem. The distance between the bird and the food is the objective function of the current particle. Each particle s has its own position vector θ s (containing components X, Y, PS, etc.) and velocity vector υ s , which are updated according to Eqs. (29) and (30) during each iteration, until the iteration termination condition is met. Similarly, in the process of solving PSO, the penalty term shown in the constructive formulae (24)-(27) is needed to deal with the constraint.
where: ω is the inertia factor, is a non-negative number, usually the size of ω, and the global and local search ability of the algorithm; c 1 and c 2 as acceleration factors.

Cat swarm algorithm (CSO)
Chu and Tsai (2007) proposed a CSO by observing the daily behavior of cats. Each cat of the algorithm corresponds to a set of solutions to the optimization problem, and the cat attributes consist of speed, fitness, and tracking or searching mode flag values (usually 0 or 1). Each cat is in its initial position, and the cat's signature values are used to determine whether the cat is in tracking or searching mode. If, in tracking mode, the position of the cat is changed by changing the speed of each dimension of the cat according to Eqs. (31) and (32); if in search mode, n cats first copy their position and then perform mutation operations, making each cat occupy a new position; the fitness is then calculated again from the 2n cat position to select the cat with the highest fitness to replace the current cat, completing the position update. Finally, the cat is placed in a new position and the optimal cat is retained until the termination condition is satisfied. Similarly, in the process of solving CSO, the penalty term shown in the constructive formula (24) ~ (27) are needed to deal with the constraint.
where: v d k (t + 1) is the d-th dimensional velocity vector of the k-th cat after updating; l is the dimension of the solution; X d best (t) is the d-th dimension component of the current optimal solution; X d k (t) is the d-th dimension component of the current position of the k-th cat; c is a constant and a random number of [0,1].

Whale optimization algorithm (WOA)
The WOA is a kind of heuristic optimization algorithm that simulates the hunting behavior of the humpback whale (Mirjalili and Lewis 2016). The algorithm mainly includes three important stages: encircling prey, catching prey with a bubble net, and searching for prey. In the algorithm, each humpback whale represents a set of solution sets, and the global optimal solution is obtained by continuously updating the whale's position in the solution space.
(1) Stalking prey. The WOA algorithm first assumes the target prey or an approximate optimal solution, i. e. the optimal whale position is defined, and then the other whales update their positions to the optimal whale according to the formulae (33)-(36).
where: D is the position of the current individual and the optimal solution; t is the number of iterations; A and C are the vector coefficients; X * (t) is the position vector of the current optimal solution. In the iterative process, a drops linearly from 2 to 0, and r 1 and r 2 are [0,1] random numbers.
(2) Bubble net predation. As the whales feed in the bubble net, they move along a spiral path. Therefore, the position update between the whale and prey is expressed by a logarithmic spiral equation, as shown in Eqs. (37) and (38). where: b is the helix shape parameter and l is the random number of the uniform distribution in the range [-1, 1]. Usually in WOA, the probability p is used to choose between bubble net predator and stalking predator.
(3) Search for prey. To ensure that the whale is fully searched in solution space, when |A|≥ 1, the searching individual moves toward the whale in a random position, as shown in Eqs. (39) and (40).
Similarly, the WOA needs to construct the penalty term shown in Eqs. (24) to (27) to deal with the constraint.

Study area
Xinyi City is located at the northern edge of Jiangsu Province, China. Its geographical position is 117° 59ʹ-119° 39ʹE and 34° 06ʹ-34°26ʹN. It belongs to an arid and water deficient area in Huaibei, Jiangsu Province. The location of the research area is shown in Fig. 4. There are two irrigation areas in the Yibei district and the Gao'a district. The irrigation area is limited by the natural conditions in the hilly and mountainous areas. The crops are mainly dry crops, usually wheat and corn varieties. The wheat planting area in the Gao'a irrigation area is 1.54 × 10 3 hm 2 , while the corn planting area is 1.63 × 10 3 hm 2 . The wheat planting area in the Yibei irrigation area is 1.26 × 10 3 hm 2 , while the corn planting area is 1.47 × 10 3 hm 2 .
The distribution of rainfall in the study area is uneven, the rainfall in the flood season accounts for 68% of the annual rainfall. The average annual evaporation is 968.0 mm, and the irrigation water use coefficient is 0.65. The main irrigation water source of the Gaotang (GT) Reservoir is 15 MCM and the lower limit of reservoir capacity is 3 MCM in the Gao'a irrigation area. The main irrigation water source of the Shadun (SD) Reservoir is 7 MCM and the lower limit of the reservoir capacity is 3 MCM in the Yibei irrigation area.
The two main dry crops (wheat and corn) planted in the study area, which was affected by topographic factors, underwent crop rotation. Considering that the growing cycle of wheat is usually from October to May while that of corn is usually from June to September, the total duration of the growth period of the two crops is basically consistent with that of the whole hydrological year (from October to September of the following year). Therefore, in this study the regulation and storage periods of the GT reservoir and SD reservoir were divided according to the growth stages of the crops. The inflow, rainfall, measured evaporation of E 601 and crop water demand of the reservoir are shown in Fig. 5. The crop parameters for the irrigation areas are shown in Tables 1 and 2.
Because of the high topography of the study area, the water demand of the irrigation area cannot be met only by gravity diversion and reservoir impoundment. Therefore, through the transformation, upgrading, and modernization of the irrigation district, new water source projects (such as pump stations) are used to ensure the irrigation needs of the district are met in drought years. The existing Shengli (SL) Fig. 4 Location of study area and Haohu (HH) pumping stations, which are designed to operate at 1m 3 /s for 20 h/day, provide supplementary irrigation to the GT and SD reservoirs by drawing water from Shu River. According to Xinyi's water resources planning data, the total amount of river agricultural water rights in a typical year is 12 MCM.
The evaporation loss water volume of the two irrigation reservoirs was determined according to the evaporation depth and the average water area of the reservoir in this period. The evaporation depth was determined by the measured evaporation data of the E 601 evaporator and corrected by the conversion coefficient w i,j . while the water area was determined according to the reservoir area storage capacity relationship function provided by the reservoir management personnel.
where: EF i,j,k is the evaporation capacity of reservoir i at the k growth stage of crop j, MCM; E i,j,k is the evaporation capacity of the evaporator, mm; w i,j,k is the conversion coefficient of the water surface evaporation of i reservoir at the k growth stage of j crops; α and β are the coefficients (GT Reservoir: α = 1.963 × 10 -3 , β = 2.013; SD Reservoir: α = 2.117 × 10 -3 , β = 1.863).

Evaluation indicators
The reliability of water supply guarantee rate of reservoir and vulnerability index of water shortage index (Chanda et al. 2014) are used to further evaluate the water supply guarantee ability of irrigation system after optimized solution of five algorithms. The water supply guarantee rate represents the satisfaction degree of irrigation system water demand, and is expressed as the ratio of actual water supply to water demand, and calculated according to Formula (42); the water shortage index represents the severity of water shortage in irrigation system, and is calculated according to Formula (43).

Algorithm optimality
In order to compare the performance of the MDADP, RGA, PSO, CSO and WOA algorithms, their optimality was compared from four aspects: the optimality of the objective function, the stability of the algorithm, the convergence of the algorithm, and the solving speed of the algorithm.
To compare the optimality of the objective function, the sensitivity analysis of the parameters of the five algorithms was needed first. As shown in Table 3, the optimal value of the objective function for the MDADP was 6.76E + 07 rmb, with no parameters to be calibrated. For the RGA, the optimal value of the objective function was 6.68E + 07 rmb, with a population size of 100 and a crossover rate of 0.6 for pm, while the variation rate pe was 0.2. For PSO, the optimal value of the objective function was (6.67E + 07) rmb, at which time the population size was 60, the inertia weight ω was 0.7, and the acceleration factor (c 1 = c 2 ) was 1.2. For the CSO, the optimal value of the objective function is 6.67E + 07 rmb, when the population size was 80, MR is 0.4, CDC and SRD were both Objective function 20 6.43E+07 1 6.49E+07 1 6.70E+07 40 6.56E+07 2 6.70E+07 2 6.63E+07 60 6.33E+07 3 6.52E+07 3 6.65E+07 80 6.62E+07 4 6.37E+07 4 6.48E+07 100 6.70E+07 5 6.25E+07 5 6.37E+07 0.3. For the WOA, when the population size was 100, the optimal value of the objective function is 6.70E + 07 rmb, the initial value of convergence factor a was 2, and the shape parameter b was 1. Although the objective function values of the above algorithms are close, the MDADP was still slightly better than the other heuristic algorithm. At the same time, the MDADP did not need calibration parameters in the solution process, while the other four intelligent algorithms all had 2-3 parameters needing sensitivity analysis. The essence of the MDADP algorithm is that many dynamic programs are nested in each other and the result is not affected by the additional parameters, so it has better operability.
According to the parameters measured in Table 3, the five algorithms were run 10 times to verify the stability of each algorithm. As shown in Table 4 and over 10 runs, the MDADP achieved an optimal value for the objective function, and still outperformed the RGA, PSO, CSO, and WOA for the optimal solution, the worst solution, the mean, and the standard deviation. Therefore, the results show that the proposed MDADP has better stability than the other algorithm. Figure 6a shows the convergence trends of the five algorithms. It is known that the running of the MDADP algorithm is independent of the number of iterations of the algorithm, and the optimal objective function value can be obtained in one run. Therefore, the MDADP converged much faster than the RGA, PSO, CSO, and WOA. Figure 6b shows the run time of the five algorithms. Although running time of the RGA, PSO, CSO and WOA algorithms increased as the number of iterations increased, they still saved compared with the MDADP. Therefore, the MDADP algorithm has the disadvantages of requiring large storage and a long computation time.

Algorithm applicability
The applicability of the MDADP, RGA, PSO, CSO and WOA algorithms was compared and analyzed. The operation process and results of the parallel 'reservoir and pumping station' irrigation system are considered here. Figure 7a and b show the water supply process of the SD and GT reservoirs after the optimization solution. It was obvious that the water supply trends of the two reservoirs obtained by the five algorithms were basically similar; there were only subtle differences in how much water was supplied at each stage. Table 5 shows that the MDADP algorithm obtained the largest total water supply and the smallest water shortage. Compared with the RGA, PSO, CSO, and WOA algorithms, the MDADP algorithm increased the water supply by 0.3%, 0.9%, 2.5% and 1.0%, respectively, and reduced the water shortage by 2.8%, 7.3%, 18.4% and 8.4%, respectively, while the total income of irrigation area increased by 1.2%, 1.3%, 1.3% and 0.9%, respectively. Figure 7c and d show the changes in the storage capacity of the two reservoirs. The capacity curves obtained by the MDADP were lower than those by the RGA, PSO, CSO and WOA, especially for the SD Reservoir. It was also evident that at the end of the storage period of the SD reservoir, the MDADP guaranteed V end = V 0 , while the other four heuristic algorithms were all V end > V 0 . The re-analysis of Fig. 7e and f shows that the main reasons leading to the termination of the reservoir was that the agricultural water rights allocated to the SD Reservoir by the RGA, PSO, CSO, and WOA are too much. However, the total water supply of the final irrigation system was less than that of the MDADP, indicating that the excess water rights allocated by the RGA, PSO, CSO, and WOA algorithms were stored in the reservoir, resulting in V end > V 0 . It also meant that heuristic irrigation systems did not work as well as reservoirs and pumping stations, and scarce water was not used effectively.
In contrast, usually the initial set of solutions of a heuristic algorithms was randomly generated in the feasible region, whereas for a linear constraint containing d decision variables, such as x 1 + x 2 + ... + x d ≤ B(B ∈ R) , if the intelligent algorithm wants to deal with this constraint, it needs to construct a penalty function to ensure that the sum of the d decision variables was less than or equal to B. Obviously, satisfying the condition that the sum of d random numbers is equal to B is the strict, therefore, most of the solutions in the initial solution set ultimately satisfied the condition that the sum of d random numbers was less than B. The above was also the key reason for V end > V 0 after solving by the RGA, PSO, CSO and WOA algorithms. The MDADP algorithm guaranteed V end = V 0 , mainly because the essence of this algorithm was the multiple nested solution process of the DP algorithm, and according to the state transition equation of DP, the sum of d random numbers was guaranteed to be equal to B. Table 6 shows the water supply guarantee rate and water shortage index of two reservoirs obtained by the five algorithms. It was clear that the irrigation system should be as large as possible in terms of the guaranteed water supply and the water shortage index should be as small as possible. Compared with the RGA, PSO, CSO, and WOA, the MDADP improved the water supply guarantee rate by 0.2%, 7.3%, 11.2% and 1.7%, respectively, and reduced the water shortage index by 3.9%, 19.4%, 2.0% and 26.0%, respectively. For GT Reservoir, the water supply guarantee rate of the MDADP was increased by 2.0%, 14.1%, and 0.5% compared with RGA, PSO and WOA, respectively, while it decreased by 4.0% compared with the CSO. Compared with RGA and WOA, the water deficit index decreased by 11.2% and 26.7%, and increased by 30.3% compared with PSO and CSO. Therefore, the MDADP algorithm can be used to solve the model to obtain a higher water supply guarantee rate and a smaller water shortage index.
Overall, it was found that the proposed MDADP algorithm was more reasonable than the RGA, PSO, CSO and WOA algorithms by analyzing the operation process and operation results of the parallel 'reservoir and pumping station' irrigation system. It was more suitable for solving the optimal allocation model of water resources in the joint operation of reservoirs and pumping stations.

Conclusion
Research on the optimal operation of reservoir water resources is of great significance importance for the comprehensive development and management of water resources. By focusing on several irrigation reservoirs and pumping stations, defining the hydraulic connection between reservoirs and pumping stations, and accounting for the restriction of river agricultural water rights, a parallel 'reservoir and pumping station' irrigation system model for the optimal allocation of water resources was established, and a new MDADP algorithm was proposed and compared with the RGA, PSO, CSO and WOA heuristic algorithms.
The optimality analysis showed that the MDADP performed better than the RGA, PSO, CSO and WOA in the optimality, stability and convergence of the objective function, moreover, the MDADP algorithm did not need the calibration parameter, and the operability of the algorithm was also good. However, the MDADP algorithm also had the disadvantages of requiring large storage and a long computation time, and needs additional methods to improve its performance. The applicability analysis showed that the MDADP algorithm was more effective than the above heuristic algorithms in improving the efficiency of water resources in terms of the scheduling process and the scheduling results of irrigation systems, increasing the total revenue of the whole irrigation area. Therefore, the MDADP algorithm is a practical optimization algorithm to solve the complex water source engineering joint dispatching problem, and popularization and application values.   Fig. 7 a Is the water supply process of SD reservoir; b is the water supply process of GT reservoir; c is the storage capacity curve of SD reservoir; d is the storage capacity curve of GT reservoir; e is the water replenishment process of HH pumping station; f is the water replenishment process of SL pumping station