Study on reservoir optimal operation based on coupled adaptive ε constraint and multi strategy improved Pelican algorithm

The optimal operation of reservoir groups is a strongly constrained, multi-stage, and high-dimensional optimization problem. In response to this issue, this article couples the standard Pelican optimization algorithm with adaptive ε constraint methods, and further improves the optimization performance of the algorithm by initializing the population with a good point set, reverse differential evolution, and optimal individual t-distribution perturbation strategy. Based on this, an improved Pelican algorithm coupled with adaptive ε constraint method is proposed (ε-IPOA). The performance of the algorithm was tested through 24 constraint testing functions to find the optimal ability and solve constraint optimization problems. The results showed that the algorithm has strong optimization ability and stable performance. In this paper, we select Sanmenxia and Xiaolangdi reservoirs as the research objects, establish the maximum peak-cutting model of terrace reservoirs, apply the ε-IPOA algorithm to solve the model, and compare it with the ε-POA (Pelican algorithm coupled with adaptive ε constraint method) and ε-DE (Differential Evolution Algorithm) algorithms, the results indicate that ε. The peak flow rate of the Huayuankou control point solved by the IPOA algorithm is 12,319 m3/s, which is much lower than the safe overflow flow rate of 22,000 m3/s at the Huayuankou control point, with a peak shaving rate of 44%, and other algorithms do not find effective solutions meeting the constraint conditions. This paper provides a new idea for solving the problem of flood control optimal operation of cascade reservoirs.


Flood control operation model of cascade reservoirs
Objective function.In order to reduce the flood control burden of the downstream reservoirs and the safety of the downstream flood control points, the flood control operation model of cascade reservoirs is established based on the maximum peak clipping criterion.The maximum peak clipping criterion not only ensures that the maximum peak flow at the downstream control point is reduced, but also makes the discharge flow of the reservoir more uniform, reducing the flood risk of the basin.
As shown in Fig. 1, suppose that the basin flood control system has N cascade reservoirs, 1 flood control point, the inflow of the first reservoir and the interval flood flow between reservoirs are known, the number of dispatching periods is T, and the maximum peak clipping objective function of the downstream control point is as follows: where q n,t is the discharge flow of the n-th reservoir in the t period; Q n,t is the interval flood flow between the n-th reservoir and the (n + 1)-th reservoir; t is the operation interval.
(1) Where V n,t−1 , V n,t is the initial storage capacity and final storage capacity of the nth reservoir in the t period;I n,t is the inflow of the n-th reservoir in the t period;q n,t is the outflow of the n-th reservoir in the t period;Q max is the maximum downstream flow allowed to ensure the safety of the downstream river;q Z n,t is the maximum discharge capacity of the nth reservoir when the initial water level is Z n,t ;Z n,0 is the initial water level at the initial time of the n-th reservoir operation;Z n,begin is the starting water level of the n-th reservoir;Z n,T is the water level at the end of the n-th reservoir operation period.Z n,end is the expected water level at the end of the n-th reservoir operation period.

IPOA algorithm based on adaptive ε-constraint
POA algorithm.In the Pelican optimization algorithm, the behavior and strategies of the pelicans during attack and hunting were simulated to update candidate solutions.The hunting process is divided into two stages: approaching the prey (exploration stage) and flying on the water surface (development stage).

Initialization.
In the equation, X ij is the position of the j-th dimension of the i-th pelican; N is the number of populations; D is the dimension of the decision variable; lb and ub are the lower and upper bounds of the decision variables respectively.
Phase 1: Exploration phase.In the first stage, the pelican determines the location of its prey and then moves towards this designated area.The mathematical model is as follows: In the formula, X new1 i,j is the position of the j-th dimension of the i-th pelican based on the first stage update; I is a random integer of 1 or 2; P j is the position of the prey in the jth dimension; F p is the objective function value of the prey; F i is the objective function value of the i-th pelican.

Phase 2: Development phase.
In the formula, X new2 i,j is the position of the j-th dimension of the i-th pelican based on the second stage update; R is a random integer of 0 or 2; iter is the number of current iterations; Maxiter is the maximum number of iterations.

IPOA algorithm.
In order to improve the performance of the POA algorithm, the following improvements are made in this paper on the basis of the POA algorithm.
Goodpoint set principle.The standard POA algorithm uses a random method to initialize the pelican population, which is highly randomized and some better pelican individuals are easily ignored.In this paper, we adopt the good point set theory proposed by Chinese mathematician Hua Luogeng to initialize the pelican population.Good point set is an effective method that can select points uniformly, and compared with the random method, the points taken using the good point set method can be more uniformly distributed in the search space.The literature 17,18 demonstrated that the population initialized with good point set theory is more uniform and can increase the diversity of the population during the initialization process.The formula is as follows: where p is the smallest prime number satisfying Fusion of reverse differential evolution strategy.The main idea of the reverse learning strategy 19 is that when searching for the optimal solution, the current solution and the reverse solution are searched simultaneously, and the optimal solution is selected by comparing the fitness values of the current solution and its reverse solution.The initial population can increase the diversity of the population 20 by adding a reverse population through the reverse learning strategy, and the reverse population solving formula is as follows: where X is the reverse solution.X is the current solution.Differential evolution algorithm 21 comes from the genetic algorithm proposed earlier, and also has the evolution process of crossover, mutation and selection.Differential evolution of pelican population after reverse learning is carried out as follows: First, each pelican individual of the current population and the reverse population was subjected to a mutation operation by Eq. ( 13) to obtain mutant individuals.where x i is the current pelican individual; u i is the mutant individual corresponding to the current pelican individual; K is the scaling factor; x r1 , x r2 are two pelican individuals randomly selected.
Then, a new pelican individual is generated by the crossover operation of Eq. ( 14).
where CR is the crossover probability.Finally, more suitable individual pelicans were selected by comparing the magnitude of the fitness values, as shown in Eq. (15).
Optimal individual adaptive t-distribution variation.The convergence of the algorithm to the local extremum depends on the optimal position of the individual 22 .Therefore, in this paper, the adaptive t-distribution variation strategy in the literature 23 is applied to the variation of the optimal pelican individual, and the current number of iterations is used as the degree of freedom of the t-distribution.At the beginning of the iteration, the t-distribution variation is close to the Coasean distribution variation, which is conducive to enhancing the search ability of the pelican individual at the global level and increasing the diversity of the population; at the end of the iteration, the t-distribution variation is close to the Gaussian distribution variation, which can enhance the search ability of pelican individuals near the optimal point and accelerate the convergence speed of the algorithm.The optimal individual adaptive t-distribution variance is formulated as follows: where X t best are mutated pelican individuals; X best are currently the best pelican individuals; iter is the number of current iterations; t represents t-distribution.
Adaptive ε constraint method.The ε-constraint method 24 is a method proposed by Takahama for solving constrained optimization problems, which retains infeasible individuals with low constraint violation and low objective function values by relaxing the constraints, and gives these excellent solutions a chance to enter the next generation population, which in turn explores to more regions and finds better objective function values.To overcome the problem that the traditional ε-constraint method tends to fall into local optimal solutions, an adaptive ε-constraint method is proposed in the literature 14 with the following improvements: 1. Improving the individual comparison criterion to make full use of good infeasible individuals to explore to more solution space, thus increasing the population diversity.The specific criterion is shown in Eq. ( 17): (11) where f(x) is the function fitness value;G(x) is the constraint violation value; Ps is a random number on the interval [0.9,1].2. Adaptive ε adjustment strategy, the ε value of the traditional ε constraint method depends only on the number of iterations, the adaptive ε value fully considers the relationship between the objective function value and the size of constraint violation and the proportion of feasible individuals in the population, and makes adaptive adjustments at each generation.ε equation is as follows: where Te is the number of truncated evolutionary iterations; is the proportion of viable individuals in the population.The value of Te should be appropriate, too large will make a large number of infeasible individuals appear in the late iteration and affect the population convergence to a feasible solution; too small will eliminate a large number of infeasible individuals in the early iteration and easily fall into a local optimum solution.
Computational flow of the ε-IPOA algorithm.The steps of the ε-IPOA algorithm are as follows, and the flowchart is shown in Fig. 2.
Step 1: Initialize the parameters, including the number of populations N, the maximum number of iterations T, the number of truncated evolutionary iterations Te, α min , α max , the boundary conditions and the dimensionality of the decision variables.
Step 2: Generate the pelican population by initializing the set of good points according to Eq. ( 11) and calculate the fitness value.
Step 3: Perform the reverse differential evolution operation on the pelican population according to Eqs. ( 12)-( 15).Step 4: Performing t-distribution variation on the optimal pelican individuals according to Eq. ( 16).
Step 5: Updating the position of the next generation of pelican individuals and calculating the fitness value.
Step 6: Comparing individuals and performing selection according to Eq. ( 17).
Step 7: Determine whether the condition to terminate the iteration is met, if so, go to step 8, otherwise go to step 3. Step 8: Output the optimal pelican individual and the optimal fitness value.

Simulation test of ε-IPOA algorithm.
To verify the effectiveness and feasibility of the POA algorithm, this article selects 24 test functions from the internationally recognized cec2006 test function set for simulation experiments, and compares them with DE and POA algorithms.Among them, each test function runs independently 50 times, with a population of 200, a maximum iteration count of 10,000, a maximum function evaluation count of 500,000, a truncated evolution iteration count of 1000, and a tolerance value of 0.0001 for equation constraint violations.The experimental results are shown in Table 3.
The experimental results are shown in the table above, and the bolded font indicates the optimal effect.From Table 1, it can be seen that: 1. Adaptive ε-constraint method can effectively assist the IPOA algorithm in handling constrained optimization problems.2. The three algorithms run independently on each function for 50 times.When the evaluation times are consistent, the average value of the results obtained by ε-IPOA algorithm is lower than that of the other two algorithms.This shows that when the Time complexity is consistent, the performance of ε-IPOA algorithm is better than that of ε-POA algorithm and ε-DE algorithm, and it has stronger global search ability.Except for the g20 and g22 functions, ε-IPOA algorithm has found the optimal solution that satisfies the constraint conditions. (17 3. The three algorithms were independently run 50 times on each function, and the standard deviation of the results obtained by ε-IPOA algorithm was smaller than that of the other two algorithms, demonstrating the good robustness of ε-IPOA algorithm.4. ε-IPOA algorithm only achieves a standard deviation of 0 when solving functions g01 and g12, indicating that this algorithm still has potential for development and needs to further improve its exploration and development capabilities.

Case analysis
Study area.The Yellow River is the second largest river in China, with a total field of 5464 km and a basin area of 795,000 km 2 , originating from the Bayankara Mountains on the Qinghai-Tibet Plateau, flowing from west to east that through nine provinces and regions, including Qinghai, Sichuan, Gansu, Ningxia, Inner Mongolia, Shanxi, Shaanxi, Henan and Shandong, and injecting into the Bohai Sea in Kenli County, Shandong Province.In this paper, the area from Sanmenxia to Huankou in the middle and lower reaches of the Yellow River is selected as the study area, and flood control and scheduling research is carried out for Sanmenxia and Xiaolangdi tandem reservoirs.The geographical location is shown in Fig. 3.The tandem reservoir flood control and scheduling model takes the storage capacity of each reservoir at each time period as the decision variable, contains constraint constraints such as water balance, safe river discharge and gate discharge flow, and is solved by the ε-IPOA algorithm.Sanmenxia Reservoir and Xiaolangdi Reservoir are both backbone reservoirs on the main stream of the middle and lower reaches of the Yellow River, mainly for flood control, taking into account the role of irrigation and power generation.Sanmenxia Reservoir controls a basin area of 688,000 km 2 , accounting for 91.5% of the total basin area of the Yellow River, and controls the flooding in the area above Sanmenxia.Xiaolangdi Reservoir controls a watershed area of 694,000 km 2 , accounting for 92.3% of the total area of the Yellow River basin, and is the only large comprehensive water conservancy project with large reservoir capacity in the middle and lower reaches of the Yellow River, except for Sanmenxia.The characteristic parameters of the two reservoirs are shown in Table 2.    www.nature.com/scientificreports/Flood process analysis.The floods in the area from Sanmenxia to Huayuankou mainly come from the upstream water of Sanmenxia reservoir, the floods between Sanmenxia and Xiaolangdi reservoir area and the floods between Xiaolangdi reservoir and Huayuankou control point.These three kinds of floods rise fiercely with high peaks and large amounts, causing a great threat to the downstream.In this paper, we select the 1000year flood of 1958, which is the largest flood in the Yellow River since the availability of measured hydrological data, mainly caused by persistent heavy rainfall.The flood process is shown in Table 3.
Figure 4 shows the flood evolution process from Sanmenxia reservoir to Huayuankou interval.It is assumed that the inlet flow of Sanmenxia reservoir is Q 1 , the interval flood from Sanmenxia to Xiaolangdi reservoir is Q 2 , and the interval flood from Xiaolangdi reservoir to the control point of Huayuankou is Q 3 .The flood process of Huayuankou section is composed of floods in each area through the action of river evolution and reservoir regulation.
The process of flood water moving from upstream to downstream in a river is called flood evolution.The study of flood evolution allows staff at downstream sites to forecast the flood process at downstream sites based on the flow process measured at upstream sites, providing a basis for downstream flood forecasting and flood control.The commonly used calculation methods are hydrological method and hydrodynamic method, the hydrological method is simple to calculate and requires less information, the hydrodynamics is limited by the measured topographic information 25 , this paper selects the hydrological method based on the tank storage equation and water balance principle-Muskingum method 26 calculation, the parameter values used are shown in Table 4, Where K is the tank storage coefficient, x is the flow specific gravity coefficient, and △t is the time interval (hour).

Results and discussion
This paper aims to minimize the peak flow of Huayuankou section, and adopts ε-IPOA algorithm solvethe model.In the solution, the population size and the maximum number of iterations are 200 and 30w, respectively, and the number of truncation iterations Te is 1000.based on the measured hydrological data, the flood calendar time and reservoir operation period in this paper are 13d, and the calculation period is 4 h.

Analysis of operation results.
According to the maximum peak-clipping objective function of the garden mouth section established in "Objective function", the ε-IPOA algorithm proposed in this paper is used to solve the flood control operation model of Sanmenxia and Xiaolangdi cascade reservoirs, and the flood process of the Huayuankou section is shown in Fig. 5 and Table 5, and the operation results of Sanmenxia and Xiaolangdi reservoirs are shown in Table 6, Figs. 6 and 7.
It can be seen from Tables 6 and 7 that after the joint dispatching, the peak flow of Huayuankou section is 12,319 m 3 /s, which does not exceed the controlled discharge of 22,000 m 3 /s of Huayuankou, and the peak clipping rate has reached 44%.Through the joint operation of the two reservoirs, the peak clipping rate have reached 27.1% and 64.1% respectively, which is conducive to the safety of the reservoir itself and reduces the flood control pressure in the downstream.Figures 6 and 7 are the operation process charts of Sanmenxia Reservoir and Xiaolangdi Reservoir respectively.It can be seen from the figure that the water level processes of the two reservoirs are between the flood limit water level and the flood control high water level, indicating that the solution results meet the constraints and reach the feasible solution.The inflow flood of Xiaolangdi Reservoir has reached a "double peak", with the maximum peak flow exceeding 25,000 m 3 /s, which greatly exceeds the safe discharge of the river channel and increases the risk of flood disaster.After regulation, the discharge flow of Xiaolangdi Reservoir is stable within 10,000 m 3 /s, and the flood process is generally stable, ensuring the safety of flood discharge of the downstream river channel.

Comparison of operation results. Algorithm comparison.
To verify the feasibility and applicability of the ε-IPOA algorithm, the ε-POA and ε-DE algorithms are chosen to solve the above model in this paper.To ensure the fairness of the algorithms, the parameters and initial conditions of the ε-POA and ε-DE algorithms are kept the same as those of the ε-IPOA algorithm.Unfortunately, neither algorithm found a feasible solution for the model, further illustrating the superiority of the ε-IPOA algorithm in solving the reservoir scheduling problem.Figures 8 and 9 show the solution results of the ε-POA and ε-DE algorithms for the Sanmenxia reservoir, respectively.
Comparison with single reservoir scheduling results.In order to illustrate more intuitively the effect of joint operation than single reservoir operation, this paper also calculates the flooding process of the Huayuankou section under single reservoir operation, and the results are shown in Tables 7 and 8.
It can be seen from the above table that under the two modes, the peak flow of Huayuankou section does not exceed the control flow, but the peak clipping effect of Huayuankou section under the joint dispatching mode is better than that of single reservoir dispatching.By comparing Tables 6 and 8, under the single reservoir operation mode, the water levels of Sanmenxia and Xiaolangdi reservoirs have not recovered to the initial water level at the end of the operation period, which is not conducive to coping with the arrival of the next flood in the flood season.In addition, under the single reservoir operation mode, the peak clipping rate of Sanmenxia and Xiaolangdi reservoirs is lower than that of joint operation.Therefore, the peak clipping effect of the joint operation mode is the best, which can play a more important role in flood control of the reservoir.angdi reservoir and Huayuankou control point are 27.1%, 64.1% and 44%, respectively, which effectively ensure the flood safety of the Huayuankou section.And the ε-POA algorithm and ε-DE algorithm with which for did not find a feasible solution.In addition, the results of the optimal operation of single reservoir are compared, and the joint operation is better than the single reservoir operation.3. The results indicate that the ε-IPOA algorithm is feasible for solving the optimization operation problem of reservoir flood control, and can effectively address the strong constraints, multi-stage, and high-dimensional problems of reservoir operation models.This algorithm provides a new approach to solve the optimization scheduling problem of reservoir groups.4. In the future, the performance of this algorithm will be further optimized through other strategies and applied to more complex series parallel hybrid reservoir groups.In addition, this algorithm is also applied to other engineering constrained optimization problems.

Figure 3 .
Figure 3. Geographical location map of the research object.

Table 1 .
Comparison results of constraint algorithms.

Table 3 .
Flooding process.Flood evolution of Sanmenxia and Xiaolangdi reservoirs.

Table 6 .
Calculation results of each reservoir.

Table 7 .
Comparison of operation results.

Table 8 .
Single reservoir operation results.