Multi-objective optimal operation of reservoir group in Jialing River based on DREAM algorithm Open Access

China is rich in hydropower resources with great potential for further development. According to the results of the census of hydropower resources, the theoretical reserves of river hydropower in China are 6.76 × 10⁸kW, of which 4.93 × 10⁸kW are technically developable installed capacity and 3.78 × 10⁸kW are economically developable capacity (Zhao et al. 2020). China ranks first in the world in terms of both the storage capacity of hydropower resources and the potential exploitation of hydropower resources (Zheng 2007). With the rapid development of society and the economy, almost all large and medium-sized rivers in China have adopted or planned to adopt cascade reservoir development. But the lack of unified river basin planning leads to disorder in reservoir development and utilization (Sun et al. 2019). Although there are many reservoirs, the utilization efficiency of water resources is rather low (Shen et al. 2019). Therefore, how to optimize the operation of reservoir groups and increase the utilization rate of hydropower resources is of great significance for realizing the maximization of China's economic and social benefits and maintaining the sustainable development of the ecological environment.

Reservoir operation is the optimal scheduling of reservoirs or groups of reservoirs based on the process of water coming from each reservoir and the existing storage capacity, using computer technology and optimization methods to solve the optimal operation scheme (Lu et al. 2015; Jiang et al. 2017; Tan et al. 2020). Reservoir operation is a strongly constrained, non-linear, multi-stage combined optimization problem (Labadie 2004). The main objective is to establish a reasonable optimal operation model and to choose the most appropriate approach (Chen et al. 2008). Two types of methods, including traditional methods and heuristic methods, have been widely applied to the optimal operation of reservoirs, while there are some problems in traditional methods, such as unstable convergence results, complex algorithms, and the ‘curse of dimensionality’ (Xiong et al. 2019). With the development of computer technology, a class of modern intelligent optimization algorithms based on biological evolution and artificial intelligence with strong robustness and versatility has been developed. These algorithms, such as the genetic algorithm (Nicklow et al. 2010), particle swarm optimization (Zhang et al. 2013; Stoppato et al. 2014), ant colony algorithm (Liao et al. 2012), differential evolution algorithm (Mohanty et al. 2014; Wen et al. 2018) and cuckoo algorithm (Kanagaraj et al. 2013), are simple in principle, easy to implement, and have a strong capability of parallel search and global optimization, which have been widely used in the optimization of reservoir operations. The fast solution method of a multi-objective and large system and a parallel algorithm using advanced computer technology for improving solution efficiency and global optimal solution ability are the critical issues of the technical research on the optimal operation of reservoir groups.

The DREAM algorithm (Differential Evolution Adaptive Metropolis algorithm) is based on Bayes' theorem. It can adaptively adjust the sampling direction and step size of the Markov chain during the sampling process to maximize the utilization of all data and information, and obtain clear posteriori distribution (Vrugt et al. 2008). The DREAM algorithm is a sampling algorithm based on the Markov Chain Monte Carlo (MCMC) method developed by Vrugt (2016). It combines the adaptive Metropolis algorithm and the differential evolution algorithm and satisfies ergodicity. The DREAM algorithm will outperform other similar sampling methods, and is better than the general optimization methods in terms of optimization in high-dimensional space when solving multivariable nonlinear cases, high-dimensional cases and multi-peak cases (Keating et al. 2010; Wöhling & Vrugt 2011; Laloy et al. 2013; Linde & Vrugt 2013; Lochbühler et al. 2014). Due to its advantages of high sampling efficiency and stable performance, it has been widely used in parameter optimization and uncertainty research of hydrological basin models (Laloy et al. 2015). However, this algorithm is rarely used in optimal reservoir operation.

In this paper, a mathematical model of reservoir group optimization and scheduling in the river basin is solved by the DREAM algorithm, and the best water supply scheduling rule scheme for the Jialing River reservoir group is selected based on the reservoir group water supply scheduling rule evaluation index system and multi-attribute decision. Furthermore, the applicability of the DREAM algorithm in optimal reservoir operation is discussed.

The Jialing River is a primary tributary on the left bank of the upper reaches of the Yangtze River. It flows through Shaanxi Province, Gansu Province, Sichuan Province, and Chongqing Municipality. Its mainstream is 1,120 km long, with a drop of 2,300 m, and an average gradient of 2.05‰ (Lei et al. 2006). Its basin area is 15.98 × 10⁴km², accounting for 9% of the basin of the Yangtze River. The Jialing River has obvious flood seasons and non-flood seasons. The flood season of the mainstream is from June to August, and the non-flood season is from September to May in the following year. There are three large reservoirs, including Baozhusi Reservoir, Tingzikou Reservoir and Shengzhong Reservoir, in the mainstreams and tributaries of the middle reaches of the Jialing River (Gao et al. 2017). The locations of these reservoirs are shown in Figure 1, and the main characteristic parameters of the three reservoirs are shown in Table 1.

The data used in this paper included annual runoff data, reservoir characteristics, and regional water demand planning data for a total of 50 years from 1961 to 2010. They are collected from The Jialing River Basin Comprehensive Planning (2013–2025), which was published by Changjiang Water Resources Commission of the Ministry of Water Resources in 2013. Considering the optimization goal is to minimize the water shortage of the water supply, the data in dry years are used for the annual runoff data.

The Markov Chain Monte Carlo method can sample any probability density function p(x) that obeys any form, which solves the sampling problem of any form of the probability density function p(x) (Lavine 1996). The core idea is to construct a Markov chain (q₀, q₁, …, q_k) to start sampling from any state, transfer the state according to the constructed Markov chain, and after iterative calculation, the sample to be drawn can be finally obtained. It should be noted that the MCMC algorithm cannot directly sample from p(x), but p(x) can be estimated by a normalized constant (Zhao 2016).

Since the transition probability from state x_t to state x_t+1 in the above equation only depends on x_t, a Markov chain is generated. After enough iteration calculations, the sample series will tend to be stable, so that the sample series (x₁, x₂, …, x_n) can be obtained, and the series of samples follows the distribution of the function π(x).

The DREAM algorithm is based on the modified Differential Evolution Markov Chain (DE-MC) algorithm (Braak 2006). The algorithm first uses sampling techniques to obtain the information of parameter samples according to the recommended distribution of the samples, and then the DE algorithm is added to the calculation process of the Markov chain. In the DREAM algorithm, the DE-MC algorithm is improved in the following aspects (Yang et al. 2016; Vrugt & Beven 2018): (1) In the DREAM algorithm, adaptive random subspace sampling technology is used to replace the random walk Metropolis sampling in the DE-MC algorithm. (2) In the DREAM algorithm, in order to ensure the diversity of parallel sequence samples, the crossover probability can be adjusted adaptively in the initialization stage. (3) The operation of useless chains is added to the DREAM algorithm. That is, if the search interval of the parallel sequence of the algorithm is meaningless to the obtained optimal solution, the useless chain can be removed, and the convergence rate of the algorithm can be improved.

The solution steps of the DREAM algorithm are as follows (Vrugt & Ter Braak 2011):

• (1)

  Determine the number n of parameters to be identified by the model, and randomly generate initial populations {xⁱ, i = 1, …, N} according to the prior distribution of the parameters; x is an N × d dimensional matrix, where i refers to the Markov chain to which the parameter belongs, N is the number of Markov chains, and d is the dimension of the parameter vector. Start the t-generation evolution of N parallel Markov chains.

• (2)

  Run the optimization model, and calculate the posterior probability distribution π(xⁱ) that the parameters obey.

• (3)
  For each Markov chain, a potential transition is generated according to Equation (7):

  

  (7)

  where zⁱ is the individual variation; γ is the scaling factor, and under normal circumstances, γ = 2.38/(2δd)^0.5; x is the model parameter; e and ɛ are both randomly generated small numbers, and e obeys a uniform distribution U(−b, b); ɛ is variance and obeys the normal distribution N(0, b′), b and b′ are both self-defined minimal values; δ represents the number of parallel chains used to generate the candidate sample; γ(δ, d′) is the scale factor; r_i(j) and r₂(n) are the serial numbers of randomly selected parallel chains.
• (4)

  Use cross-probability P_c to cross-process individuals among populations.

• (5)
  Calculate the probability distribution π(zⁱ) and acceptance probability α(xⁱ, zⁱ) that the parameters follow, where y refers to the observed result of the model:

  

  (8)
• (6)

  According to the acceptance probability α(xⁱ, zⁱ), judge whether to accept zⁱ or reject zⁱ. If α(xⁱ, zⁱ) ≥ u, the new candidate sample is accepted; otherwise, it is not accepted. End the t-generation evolution of N parallel Markov chains.

• (7)

  Remove useless chains according to the IQR (Inter-Quartile Range) method.

• (8)
  Use the sampling points of the last 50% of each Markov chain to calculate S_R (proportional reduction factor) to determine whether the Markov chain has converged:

  

  (9)

If, for every dimension of the Markov chain, there is < 1.2, it means that the parameters converge to a stable posterior distribution, and the calculation is terminated; otherwise, the calculation continues.

Constraints on the optimal operation of cascade reservoirs are mainly classified into equality constraints and inequality constraints (Bai et al. 2015). Equality constraints include water balance and flow balance constraints; inequality constraints include water level, discharge flow, and output constraints (Teegavarapu et al. 2013).

The DREAM algorithm was applied to optimize the water supply operation model of the reservoir group of Jialing River. The basic processes are shown in Figure 2. The hydrologic frequency analysis is based on the data from Beibei Hydrologic Station covering the period from 1961 to 2010, and the P-III frequency curve is used for the adaptation. The annual runoff was 646 × 10⁸ m³. The years with a design guarantee rate of 62.5–87.5% are classified as partial dry years, and the years with a rate greater than 87.5% are classified as extremely dry years, resulting in 16 partial dry years and four extremely dry years. The data of these 20 years were used as the input for the model. The annual incoming water volume of the dry years and the planned annual design water consumption are shown in Figure 3. The annual incoming water volume reaches the maximum in July, August, and September, while the water demand is greater in May, June, and July. The user's water demand has a certain elasticity. Within the elastic range, water-saving and other measures can be used to alleviate the lack of water supply in the reservoir, and a small amount of water shortage will not cause greater losses. However, deep water shortage beyond the elastic range of user demand will cause significant economic losses, so the value of the variable HF_m,t should not be too large. In this paper, the range of the variable HF_m,t is set to [0.1, 0.3]. Under the condition of monthly calculation, for a single reservoir, the water supply restriction rules have 36 parameters, and as such the number of parameters of the water supply operation multi-objective optimization model of the reservoir group of Jialing River is 108.

In the DREAM model, there are ten Markov chains in total, and the number of parallel chain evolutions is 4,000. It is assumed that the prior distribution of each parameter is a uniform distribution. The initial population is randomly generated according to the prior distribution of parameters, which is a 108 × 10 matrix. After that, the Markov chain iterative evolution calculation process starts. The joint operation optimization model of the reservoir group system is run to calculate the posterior probability distribution of the objective function and parameters. For each Markov chain, a potential transition is generated according to Equation (7). Then, according to the crossover and variation between different Markov chain data, the cross-probability is set to 0.8, and the variation probability is set to 1%. The transition probability of the obtained population is calculated and it is determined whether the Markov chain transfers, and if the transition probability is greater than the judgment standard, then the new sample is accepted. Thus the iterative evolution calculation of this Markov chain is completed. According to the IQR method, the useless chains in the model are statistically removed, and the last 50% of the sampling points of each Markov chain are taken to calculate its value to determine whether the Markov chain has converged. If it converges, the calculation is stopped, otherwise the next iteration is entered.

First, the single objective optimization calculation of the water supply operation rules of the reservoir group in Jialing River was carried out to obtain the extreme value of the single objective function. After analyzing the convergence of the parameter sampling results, it is found that when the parallel Markov chain evolution algebra reaches 100 generations, the objective function results are basically stable, and the convergence judgment index S_R of each parameter is less than 1.2, which indicates that all parameters converge to a stable posterior distribution. According to the results of the optimization calculation, the minimum value of the system's total water shortage rate MDR is 8.3%, and the minimum value of the damage vulnerability TDR of the water supply task during a single reservoir period is 34.4%. The convergence process of the objective function is shown in Figure 4. As can be seen from the data in the figure, the convergence efficiency of the DREAM algorithm is high.

The multi-objective optimization calculation of water supply operation rules of the reservoir group of Jialing River was carried out. The convergence of the parameter sampling results was analyzed, and it was found that when the parallel Markov chain evolution reaches 1,000 generations, the convergence index S_R of each parameter is less than 1.2, indicating that all parameters converge to a stable posterior distribution. The ten bits of Markov chain convergence were selected for plotting, with a total of 100 sets of data. According to the results of the optimization calculation, the water shortage rate of the system's total water supply ranges from 8.38% to 10.05%, and the damage vulnerability of the water supply task during a single reservoir period ranges from 34.73% to 62.50% (see Figure 5). The figure shows a set of alternative schemes for the optimal water supply operation rules of the three reservoirs under the runoff conditions of dry years. In dry years, reservoir operation management decision-makers can choose specific preference schemes from the set of alternative schemes to guide the operation of each reservoir. In addition, the results of the last 1,000 generations of statistical Markov chain sampling results were statistically analyzed. The result of the histogram of the objective function is shown in Figure 6. If the shape of the histogram is high, the optimal value of the parameter is easy to obtain; if the histogram of the posterior distribution is flat, the optimal value of the parameter is difficult to determine, and the uncertainty of the parameter is large.

The index system composed of the evaluation index system of water supply operation rules was used to perform the optimal calculation on the multi-objective optimization results. The results are shown in Table 2. If the water supply operation of Baozhusi Reservoir in dry years is guided by any operation rule scheme, the reliability α of water supply and resilience γ of the reservoir are both 1, and the vulnerability v of the water supply of the reservoir and the user's generalized shortage index (GSI) are equal to 0. It shows that Baozhusi Reservoir can fully meet the planned water supply requirements in dry years, and there is no preference for any scheme. Correspondingly, the four evaluation indexes of the water supply performance of Shengzhong Reservoir and Tingzikou Reservoir fluctuate greatly, which is influenced by different schemes, and the number of alternative operation schemes needs to be further reduced.

The water operation schemes of Shengzhong Reservoir and Tingzikou Reservoir are ranked and selected according to the reservoir evaluation indexes. The final preference scheme for Shengzhong Reservoir is shown as Scheme 1 of Table 3, and the final preference scheme for Tingzikou Reservoir is shown as Scheme 2 of Table 3. It can be found that the finally selected preference scheme is the result after weighing among the evaluation indexes. When the water supply reliability α is high, the water supply damage recovery ability γ is strong, the water supply damage vulnerability v is low, and the generalized water shortage index GSI is small. Comparing with the single-objective optimization results, both Scheme 1 and Scheme 2 can not only increase the reliability of water supply α and resilience γ of Shengzhong Reservoir and Tingzikou Reservoir, but also reduce the size of the generalized water shortage index GSI.

In this paper, the DREAM algorithm was applied to research on the multi-objective optimal operation of cascade reservoirs. Comparing with the traditional optimal operation of cascade reservoirs, the DREAM algorithm satisfies ergodicity based on the combination of the adaptive Metropolis algorithm and the differential evolution algorithm, which makes it more suitable for the problems of joint operation of reservoir groups troubled with multi-variable nonlinearity and high dimensionality. The DREAM algorithm provides a compelling new way for research on the multi-objective optimal operation of cascade reservoirs, and was applied to optimize the practical multi-objective model of water supply operation for the three reservoirs of Jialing River: Baozhusi Reservoir, Shengzhong Reservoir, and Tingzikou Reservoir. The evaluation index system and multi-attribute decision-making of the water supply operation rules of the reservoir group are further applied to select the best scheme of water supply operation rules for the reservoir group of Jialing River. It is found that the Baozhusi Reservoir can fully meet the planned water supply needs in dry years, and the two reservoirs, Shengzhong and Tingzikou, need to weigh the water supply operation evaluation indexes. Therefore, the elaborated theoretical basis of the DREAM algorithm and the available operation scheme for the reservoir group of Jialing River presented in this research introduce a new feasible solution for problems of optimal operation of cascade reservoirs.

The authors are grateful to the financial support from Natural Science Foundation Project of CQ CSTC.

Methodology, W.D.; Investigation, W.D.; Resources P.P.; Funding acquisition, C.Z.; Software, W.D.; Supervision, X.Z.; Validation, C.Z.; Writing – original draft, P.P.; Writing – review & editing, S.Y.

This research was funded by the National Key R&D Program of China grant number 2016YFC0402004 and Natural Science Foundation Project of CQ CSTC grant number (KJQN201800712).

The authors declare no conflict of interest.

All relevant data are included in the paper or its Supplementary Information.