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

With the implementation of the most stringent water resources management system and the advancement of the construction process of reservoir terrace basins, the research and application of the theory and method of joint operation of reservoir groups are becoming more and more important. The Differential Evolution Adaptive Metropolis (DREAM) algorithm is a sampling algorithm based on the Markov Chain Monte Carlo method proposed in recent years. The algorithm satis ﬁ es ergodicity and is good at handling problems with multivariate nonlinearity, high dimensionality, and multi-peaks, and as such the algorithm is a new global optimization solution. This paper elaborated the solution mechanism of the standard DREAM algorithm, and the algorithm was applied to the optimal operation model of the reservoir group in Jialing River. First, we optimized and analyzed the multi-objective supply operation model of the reservoir group water in Jialing River. Then the multi-attribute decision-making and evaluation index system of water supply operation rules for the reservoir group to assess the optimization of the operation model was adopted. Finally based on the results of the evaluation, the best water supply operation scheme for the reservoir group of Jialing River was selected. The results show that the Baozhusi Reservoir can fully meet the planned water supply requirements in dry years, while the two reservoirs, Shengzhong and Tingzikou, need to be weighed against the evaluation indexes of water supply operation. The research provides a theoretical basis for the DREAM algorithm in the optimal operation of the reservoir group and the actual operation scheme for the reservoir group of Jialing River.


INTRODUCTION
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 8 kW, of which 4.93 × 10 8 kW are technically developable installed capacity and 3.78 × 10 8 kW are economically developable capacity (Zhao et al. ).China ranks first in the world in terms of both the storage capacity of hydropower resources and the potential exploitation of hydropower resources (Zheng ).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. ).Although there are many reservoirs, the utilization efficiency of water resources is rather low (Shen et al. ).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. ; Jiang et al. ; Tan et al. ).Reservoir operation is a strongly constrained, non-linear, multi-stage combined optimization problem (Labadie ).The main objective is to establish a reasonable optimal operation model and to choose the most appropriate approach (Chen et al. ).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  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.

MATERIALS Study area
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. ).Its basin area is 15.98 × 10 4 km 2 , 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. ).The locations of these reservoirs are shown in Figure 1, and the main characteristic parameters of the three reservoirs are shown in Table 1.

Data
The data used in this paper included annual runoff data, reservoir characteristics, and regional water demand plan-

Bayesian principle and Markov Chain Monte Carlo method
Bayesian statistics play an important role in mathematical statistics.The main point of its view is that for the observed value Y ¼ [y 1 , y 2 , …, y n ], any unknown parameter X ¼ [x 1 , x 2 , …, x n ] can be regarded as a random variable, which can be described by a probability distribution (Kadane ).The prior distribution h(X ) can be determined according to the prior information (Marin et al. ), and then the posterior distribution can be obtained by inference estimation, namely: For the i-th component x i in X, the probability value where X [Ài] is the remaining variables after removing the i-th . Thus, the conditional probability distribution can be obtained.
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 ).The core idea is to construct a Markov chain (q 0 , q 1 , …, 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 ).
Assume that X t is a random variable at time t, and its state space is composed of all possible X values.The value of the X state at the next time is related only to the X t at the current time, and not related to other periods (Koch ).This can be expressed as: The above is a Markov process, and the random sequence (X 0 , …, X n ) generated by the Markov process is a Markov chain.As such, the probability of a Markov chain The Metropolis-Hasting algorithm is the main method in the Markov Chain Monte Carlo algorithm.This method was proposed by Metropolis in 1953 and then further improved by Hasting in 1970 (Betancourt ).This method obtains the series of samples by establishing a Markov chain, and the sampling calculation steps are as follows: (1) Suppose the samples are from the distribution π(x), where π(x) ¼ h(x)/k, k is a constant, and it is difficult to be solved.Choose any transition probability q(x, y), and initialize the value of x.This can be expressed as: When the process is x t ¼ x, a potential transition process x → y is generated by q(x, y), where q(x, y) is called the proposal distribution, and its purpose is to ensure the stationary distribution of π(x) in p(x, y).Among them, α(x, y) ¼ h(y)q(y, x) h(x)q(x, y) .
(3) Randomly extract a number u from [0, 1], then: (4) Determine whether the termination condition is met; if not, return to step (1) and continue.
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 1 , x 2 , …, x n ) can be obtained, and the series of samples follows the distribution of the function π(x).

DREAM algorithm and its solving steps
The DREAM algorithm is based on the modified Differential Evolution Markov Chain (DE-MC) algorithm (Braak ).
The algorithm first uses sampling techniques to obtain the (2) Run the optimization model, and calculate the posterior probability distribution π(x i ) that the parameters obey.
(3) For each Markov chain, a potential transition is generated according to Equation (7): where z i 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 distri- (4) Use cross-probability P c to cross-process individuals among populations.
(5) Calculate the probability distribution π(z i ) and acceptance probability α(x i , z i ) that the parameters follow, where y refers to the observed result of the model: (6) According to the acceptance probability α(x i , z i ), judge whether to accept z i or reject z i .If α(x i , z i ) ! u, the new candidate sample is accepted; otherwise, it is not accepted.End the t-generation evolution of N parallel Markov chains.
(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: where g is the evolution algebra of each Markov chain, q is the number of Markov chains used for evaluation, B/g is the variance of the average of q Markov chains, and W is the average of the variance of q Markov chains.
If, for every dimension of the Markov chain, there is R j < 1.2, it means that the parameters converge to a stable posterior distribution, and the calculation is terminated; otherwise, the calculation continues.

Mathematical model for optimal operation of cascade reservoirs
The cascade reservoir group operates jointly as a whole, which not only has the compensation benefit of reservoir During the scheduling of water supply operation of reservoir groups, there are two goals: the total water shortage rate of the system and the damage vulnerability of the water supply task in a single reservoir period (Zhou ).The objective functions are as follows: Reservoir water balance constraint: Restrictions on upper and lower limits of reservoir storage capacity: Restrictions on upper and lower limits of reservoir supply capacity: R m,t !0 (14) Rules of water supply operation: where S m,t and S m,tþ1 represent the beginning and ending water storage capacity of the m-th reservoir in the t-th period; I m,t represents the incoming water volume of the m-th reservoir in the t-th period; E m,t represents the losses of the m-th reservoir in the t-th period, that is, evapotranspiration and seepage; SP m,t represents the surplus water in the m-th reservoir in the t-th period; S min m,t and S max m,t respectively represent the upper limit and lower limit of the water storage capacity of the m-th reservoir in the t-th period; D min m,t represents the minimum water demand of the m-th reservoir in the t-th period; A j represents a plan for the water supply operation rules of the reservoir; SWA m,t , EWA m,t , and HF m,t are the three parameters of the m-th reservoir water supply restriction rules; PE m,t is the distribution coefficient of common water supply in parallel reservoirs, namely Shengzhong Reservoir and Tingzikou Reservoir; the limiting conditions of each parameter are 0 SWA m,t D m,t , D m,t EWA m,t D m,t þ K m , 0 HF m,t 1 and 0 PE m,t 1; and K m is the regulating ability of the m-th reservoir, as shown in Table 1.

Evaluation index system of water supply operation rules for reservoir group
The reliability of water supply is generally defined as the probability or frequency that the water supply status of the reservoir is normal during the analysis period of reservoir operation (that is, no damage to water supply, and the water supply of the reservoir is equal to the capacity of water storage) (Ji et al. ).It is usually obtained by dividing the number of hours of normal water supply by the total number of hours of reservoir operation analysis periods, namely where α is the water supply reliability of the reservoir (Lu where γ is the reservoir water supply damage recovery resilience; where v is the damage vulnerability of reservoir water supply; and DR t ¼ (D t À R t )/D t is the relative water shortage in the t-th period.
The deficit percent day index (DPD) is used to indicate the two characteristics of the duration and intensity of water shortage (Chen & Chen ): where NDC represents the number of days of water shortage event; and DDR represents the daily averaged water shortage rate (%) of water shortage time.
The generalized shortage index (GSI) can not only consider the duration and intensity information of specific water shortage events during the year, but also the annual water shortage rate (Hsu ).The calculation formula is shown below: where β is a constant index, and generally is 2; DY i represents the number of days in year i; and DPDa i represents the sum of DPD values for all water shortage events in year i.

RESULTS AND DISCUSSION
Solution for optimal operation of cascade reservoir 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 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  and Tingzikou Reservoir, but also reduce the size of the generalized water shortage index GSI.

CONCLUSIONS
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 al. ).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. ), particle swarm optimization (Zhang et al. ; Stoppato et al. ), ant colony algorithm (Liao et al. ), differential evolution algorithm (Mohanty et al. ; Wen et al. ) and cuckoo algorithm (Kanagaraj et al. ), 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. ).The DREAM algorithm is a sampling algorithm based on the Markov Chain Monte Carlo (MCMC) method developed by Vrugt ().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 multipeak cases (Keating et al. ; Wöhling & Vrugt ; Laloy et al. ; Linde & Vrugt ; Lochbühler et al. ).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. ).However, this algorithm is rarely used in optimal reservoir operation.
ning 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.

Figure 1 |
Figure 1 | The location of reservoirs for Jialing River.
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. ; Vrugt & Beven ): (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 ): (1) Determine the number n of parameters to be identified by the model, and randomly generate initial populations {x i , 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.
bution U(Àb, b); ε is variance and obeys the normal distribution N(0, b 0 ), b and b 0 are both self-defined minimal values; δ represents the number of parallel chains used to generate the candidate sample; γ(δ, d 0 ) is the scale factor; r i ( j) and r 2 (n) are the serial numbers of randomly selected parallel chains.
capacity and hydrology but is also needed to meet the demand of industrial and agricultural water (Ni et al. ).Shengzhong and Tingzikou share the water demand for industrial and agricultural water below Shengzhong and Tingzikou Reservoir, while Baozhusi Reservoir independently covers the water demand for industrial and agricultural areas between Baozhusi and Tingzikou.
11) where m represents the number of a specific reservoir; t represents the time; D m,t represents the planned water supply of the m-th reservoir in the t-th period to meet the downstream industrial and agricultural water demand; and R m,t represents the actual water supply of the m-th reservoir in the t-th period.Constraints on the optimal operation of cascade reservoirs are mainly classified into equality constraints and inequality constraints (Bai et al. ).Equality constraints include water balance and flow balance constraints; inequality constraints include water level, discharge flow, and output constraints (Teegavarapu et al. ).
represents the number of times the reservoir water supply is recovered from the damaged state to the normal state; Q T t¼1 (R t < D t ) represents the number of hours of periods during which the reservoir water supply is in a damaged state; and T is the total periods during the reservoir operation analysis period.The damage vulnerability of water supply can generally be expressed by the damage vulnerability of reservoir water supply in a single period (Liu et al. ).The damage vulnerability of reservoir water supply in a single period represents the maximum relative water shortage in a single period during the reservoir operation analysis period (Christodoulou & Agathokleous ).The damage to the water supply of the reservoir under continuous drought conditions is inevitable.Water shortage may exist in several periods at the same time.The damage vulnerability of water supply (Guo et al. ) in a single period can be used to measure the severity of the water shortage, namely

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 8 m 3 .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.

Figure 2 |
Figure 2 | Multi-objective optimization process of water supply operation rules of reservoir group of Jialing River.

Figure 3 |
Figure 3 | Data of average annual incoming water volume and planned annual water consumption of dry year of each reservoir.

figure
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.

Figure 4 |
Figure 4 | The convergence of the single-objective function.
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.

Table 1 |
Main characteristic parameters of the reservoirs in Jialing River al. ); D t is the water storage capacity of the user; R t ) refers to the conditional probability that the water supply of the reservoir will recover from the damaged state to the normal state when the water supply is damaged during a certain period.It is generally expressed by the average frequency of the water supply of the reservoir recovered from the damaged state to the normal water supply state during a certain operation analysis period(Moy et al.