Enhanced harmony search algorithm for sustainable ecological operation of cascade hydropower reservoirs in river ecosystem

With the merits of superior performance and easy implementation, the harmony search (HS), a famous population-based evolutionary method, has been widely adopted to resolve global optimization problems in practice. However, the standard HS method still suffers from the defects of premature convergence and local stagnation in the complex multireservoir operation problem. Thus, this study develops an enhanced harmony search (EHS) method to improve the HS’s search ability and convergence rate, where adaptive parameter adjustment strategy is used to enhance the global search performance of the swarm, while the elite-learning evolutionary mode is used to improve the converge trajectory of the population. To verify its practicability, EHS is applied to solve numerical optimization and multireservoir operation problems. The results show that EHS can produce better results than several existing methods in different cases. For instance, the mean objective of EHS is improved by about 23.9%, 28.7% and 26.8% compared with particle swarm optimization, differential evolution and gravitational search algorithm in 1998–1999 typical runoff case. Hence, an effective optimizer is developed for sustainable ecological operation of cascade hydropower reservoirs in river ecosystem.


Introduction
In recent years, to guarantee the energy and water demand of modern society, a large number of hydropower reservoirs are being put into operation throughout the world [1][2][3][4][5]. In practice, reservoirs not only effectively increase the economic and social benefits [6][7][8][9], but also bring some negative impacts to river ecosystem health, like runoff characteristics variation and topographical change [10][11][12][13]. As a result, many aquatic organisms may die or even become extinct when their basic water demands are not well satisfied for a long time. In order to alleviate this problem, the environmental and ecological protection is gaining more and more attention all over the world [14][15][16]. Recently, many engineering measures (like dredging waterway and water conservancy projects) and non-engineering measures (like laws and regulations) have been attempted [17][18][19][20][21][22]. Based on the previous researches, ecological operation of cascade hydropower reservoirs is one of the most effective and economic measures in promoting the eco-civilization construction [23][24][25][26]. Thus, the goal of this study is to determine the suitable operation policies of multireservoir system to minimize the total ecological shortage during the scheduling horizon.
From the mathematical viewpoint, the ecological operation of cascade hydropower reservoirs belongs to a high-dimensional and multi-stage constrained optimization problem with a variety of physical constraints [27][28][29][30]. Over the past decades, many methodologies have been successfully proposed to address engineering optimization problems [31][32][33][34], like linear programming [35][36][37], nonlinear programming [38][39][40], dynamic programming [41][42][43], network-based optimizer [44][45][46] and decomposition coordination tool [47]. These traditional methods have enjoyed a great deal of success in water resources system [48], but still suffer from some drawbacks in large and complex engineering problems, like long execution time, huge memory usage, parameter tuning and duality gap [49][50][51][52]. To relieve the defects rooted in the traditional methods, many evolutionary algorithms (EAs) inspired by biological learning, natural law or social behaviors are developed, like genetic algorithm (GA) [53,54], particle swarm optimization [55][56][57], gray wolf optimizer [58], differential evolution [59], ant colony optimization algorithm [60], gravitational search algorithm [61] and cooperation search algorithm [62,63]. Without requiring gradient-dependent information, EAs can provide satisfying solutions for nonconvex or even discontinuous problems, promoting their extensive applications in numerous engineering fields [64]. However, due to the loss of swarm diversity, the original EAs may suffer from the local stagnation shortcomings or premature convergence in varying degrees [65]. Hence, scientists are devoted to developing effective tools to enhance the performances of EAs in engineering problems.
As a typical metaheuristic EAs, the harmony search (HS) tries to mimic the musician's behavior in seeking for appropriate harmony [66]. With the advantages of low computation burden, easy implementation, and satisfying adaptability, HS has been widely employed in many engineering problems, like parameter identification and feature selection [67][68][69]. Nevertheless, the standard HS method usually fails to improve the obtained solutions at the late search process, resulting in the premature convergence problem [70][71][72]. To overcome the HS disadvantage, this paper aims at introducing an enhanced harmony search (EHS) based on two modified strategies: the adaptive parameter adjustment strategy for improving global search capability of the population, and the elite-learning evolutionary mode for improving the convergence rate of the swarm. To test its performance, EHS is used to solve numerical functions and ecological operation of multireservoir system. The results demonstrate the superiority of the proposed method in providing satisfying solutions. Finally, the main novelty of this paper lies in: an effective EHS method using multiple improvement strategies is provided to obtain satisfying scheduling schemes for reservoir operation.
This study is organized as below: section 2 presents the technical details of EHS. Section 3 tests the EHS performance in benchmark functions. Section 4 gives the details of EHS in solving multireservoir operation problem. Section 5 tests the EHS feasibility in a real-world hydropower system. Section 6 presents the conclusions.

HS algorithm
HS is a novel population-based derivative-free metaheuristic algorithm developed to resolve the global optimization problem [73][74][75]. In HS, each decision variable is regarded as the pitch of a musical instrument; the pitch range of musical instrument is seen as the changing range of each variable; the solution with numerous variables is regarded as a harmony vector; a set of solutions are generated to search for the global optimum. In the evolutionary process, the pitches with better objectives are conserved in harmony memory (HM), and then the swarm can iteratively approach the optimal solution under the guidance of elite solutions [76][77][78]. For the sake of simplicity, the goal is set to find the best solution for the d-variable minimization problem, and then the search procedures of the HS method can be summarized as below: Step 1: Define the values of the computational parameters, like the harmony memory size (HMS), harmony search considering rate (HMCR), pitch adjusting rate (PAR), bandwidth (bw) and maximum iteration (MI).
Step 2: Randomly generate the initial HM in the feasible search space by equation (1). Then, each harmony x i denotes a d-dimensional solution for the target optimization problem where x i,j is the jth element in the ith solution. φ (a, b) is used to generate the random number uniformly distributed in [a, b].x j and x j are the upper and lower limit of the jth element.
Step 3: Improve the elements of all the current solutions vectors. After the initialization process, all harmony vectors are dynamically updated by three operators, including memory consideration, pitch adjustment and random selection. The detailed search process for each solution is presented as below: Step 4: Update the HM. Once the performance of the newly-obtained solutions x new is better, x new replaces the worst solution in HM.
Step 5: Repeat steps 3 and 4 to gradually improve the quality of the entire swarm until the terminal condition is satisfied. Then, the best harmony is seen as the final solution for the target problem.

EHS algorithm
For the standard HS method, the swarm may converge to local optima with a high probability due to the loss of individual diversity. To overcome this defect, EHS based on adaptive parameter adjustment For each dimension j = 1 to d then If φ (0, 1) < HMCR then //memory consideration If φ (0, 1) < PAR then //pitch adjustment xnew ,j = x k,j , k∈ [1,HMS] //randomly choose from HM Else and elite-learning evolutionary mode is proposed in this section.

Adaptive parameter adjustment strategy
To enhance the HS performance, the adaptive parameter adjustment strategy is used to update dynamically two key parameters used to exploit the local optima, PAR and bw. Specially, PAR is gradually changed from the minimum value PAR min to the maximum value PAR max by equation (2), which can help enlarge the search space and avoid falling into local minimum where gen denotes the current iteration number. Generally, a large value of bw at the early evolutionary stage will improve the solution diversity while a small value of bw at the end will enhance the refined search around the optimal solution [79]. In order to achieve this goal, bw is decreased exponentially from the maximum bandwidth bw max to the minimum bandwidth bw min , which can be expressed as below:

Elite-learning evolutionary mode
Generally, the global-best agents have large probability to approximate the theoretically optimal solution. Based on this point, the famous particle swarm optimization makes full use of the global best-known solution of the swarm as well as the personal bestknown solution of each particle to enhance the global exploitation and local exploration. Obviously, the standard HS method fails to do the operation in the search process. Thus, the elite-learning evolutionary mode is introduced to promote information interaction and experience interchange [80][81][82], which can be expressed as: where pbest i is the ith personal best-known solution. gbest ind is the randomly-chosen solution chosen from m global best-known solutions found by far, and there is ind ∈ {1, 2, · · · , m}.

Overall execution procedures
The detailed process of EHS for each solution is presented as below: Generating a new solution in EHS.

Statistical results analysis
To verify the EHS performances, the results of several methods in [81] are introduced for comparison, including GA, particle swarm optimization (PSO), biogeography-based optimization (BBO), flower pollination algorithm (FPA), moth-flame optimization algorithm (MFO), bat algorithm (BAT) and differential evolution (DE). The results of HS and EHS are obtained in the same parameters setting: HMS = 50, MI = 1000, HMCR = 0.9, PAR = 0.3, bw = 0.001. To avoid the potential negative effect of random seeds, each method is run 20 independent times. The data in tables 1 and 2S shows that for some functions (F 1 -F 6 , F 10 -F 11 and F 13 -F 19 ), EHS can find the best average objective value among all the methods; meanwhile, for other test functions, the EHS results are better than or close to other methods. Thus, EHS can yield satisfying near-optimal solutions for 24 benchmark functions.

Box plot analysis
The box plot in figures 1 and 1(S) shows that the objective values of EHS exhibits fewer outliers and smaller distributions than the traditional HS method in most test functions. Hence, the feasibility of two modified strategies in improving the HS performance is proved in this section.

Convergence trajectory comparison
The convergence trajectory of the two methods for several classic benchmark functions shows that the standard HS method cannot find satisfying results which are close to the best objective value from beginning to end (figure 2, and see supporting information figure 2(S) for details). Besides, EHS can quickly find high-quality solutions at the early search stage and then continually improve the objective values with the increasing iterations. The reasons for the superior performance of EHS lies in: the adaptive parameter adjustment strategy can dynamically change the search range to improve the population diversity, while the elite-learning evolutionary mode helps the swarm gain abundant information from excellent agents to enhance the convergence rate. To be mentioned, for problems F 1 -F 13 , the differences of two methods can be easily distinguished since the objective values vary in large range. For instance, the objective value of EHS is reduced from the initial value at 10 4 magnitude to the final value at 10 −8 magnitude. Meanwhile, the other problems (F 14 -F 24 ) have fixed variables and the variation range of the objective values are small. Thus, EHS has larger probabilities to produce better results for numerical functions as compared to the HS method Figures 3 and 3(S) shows that all the solutions are randomly distributed in the problem space at the initial stage, and then all the solutions gradually converge to the promising area around global optima; besides, the average fitness values are quickly reduced as the iterations increase. Hence, EHS achieves a compromise between exploration and exploitation regardless of problem features.

Wilcoxon nonparametric statistical test
Here, the Wilcoxon nonparametric statistical test is employed to examine the EHS performance because it does not require the distribution information of the studied samples [83]. From table 2, it can be found that the R + value is larger than the R − value in all the comparisons; besides, all the p values in all the comparisons are also smaller than 0.05. Thus, EHS is statistically significantly better than other methods.

Benchmark functions Set II: 28 CEC2013 problems
In this section, 28 CEC2013 problems with ten variables are used to test the EHS performance, which can be roughly divided into three different kinds of groups: unimodal functions (F 1 -F 5 ), multimodal functions (F 6 -F 20 ) and composition functions (F 21 -F 28 ). It becomes much more difficult for EAs to find feasible solutions due to the rough and uneven response surface of the CEC2013 functions (figure 4(S)). The data in tables 3 and 3S show that for unimodal functions (F 1 -F 5 ), EHS is always superior to the standard HS method in all the indexes; for other functions, EHS has stronger abilities than HS to jump out of local optima. Meanwhile, the average execution time of HS and EHS are neck and neck, demonstrating satisfying execution efficiency. Thus, EHS can yield feasible solutions for the global optimization problem.

Mathematical model 4.1.1. Objective function
Under the market environment, the goal of hydropower enterprise is often set to maximize the generating benefit and then reservoirs tend to raise the water levels by reducing the outflows at some periods important for aquatic organism, which may break the ecological equilibrium. With the growing sustainable development awareness, the ecological operation of multiple hydropower reservoirs is gaining increasing attention throughout the world. To effectively protect the river ecosystem health, the operational goal is often set to minimize the total water shortage volume of all the reservoirs involved in the considered hydropower system. Then, the objective function can be mathematically expressed as below:  where F is the total water shortage volume of all the involved reservoirs. N and T are the number of reservoirs and periods. E max i,t and E min i,t are the maximum and minimum ecological flow of reservoir i at period t. O i,t and y i,t are the total outflow and improper flow of reservoir i at period t.

Physical constraints (a) Water balance equations
where V i,t , q i,t , I i,t , Q i,t and S i,t denote the storage capacity, local inflow, total inflow, power discharge and surplus water of reservoir i at period t. U i denotes the set of upstream reservoirs directly connected to reservoir i.
(b) Initial and final forebay water level equations where Z where H i,t is the hydraulic head of reservoir i at period t. Z i,t and d i,t denote the forebay and downstream water levels of reservoir i at period t. (d) Forebay water level limits where Q min i,t and Q max i,t are the minimum and maximum power discharges of reservoir i at period t.
where O min i,t and O max i,t are the minimum and maximum total outflows of reservoir i at period t. (g) Power output limits where P min where f 1 i (·), f 2 i (·) and f 3 i (·) are nonlinear stagestorage, stage-discharge and generation curve of reservoir i.

Detailed technical procedures 4.2.1. Solution structure and swarm initialization
For each hydropower reservoir, the detailed results in the entire scheduling horizon can be directly deduced when the water level per period is known. As a result, the decision variables are chosen as the water levels of all the reservoirs. Then, each solution in EHS becomes a N × T matrix presented as below: For the initial swarm, the variables of each individual solution will be randomly determined in the viable search space between minimum and maximum water levels, which can be described as:

Constraint handling method
Generally, the operational scheme of hydropower reservoir is affected by many factors, like the outflows of upstream reservoirs, transmission capacity of power system and water demands of production departments. In other words, a slight modification of reservoir operation such as changing water levels may affect the final objective value of cascade hydropower reservoirs. Hence, an effective constraint handling method is of great importance to guarantee the solution feasibility [84,85]. For the reservoir operation problem, the considered constraints can be roughly divided into two categories: equality constraints (like water balance equation) or inequality constraints (like boundary limits of water level or total outflow). During the evolutionary process, the infeasible variables of each solution are forced to locate at the nearest boundary value before using equation (17) to obtain its fitness value, where the constraint violation value θ is merged into the objective value. Then, infeasible solutions are dominated by feasible solution with a large probability, and the infeasible solution with larger θ value is dominated by infeasible solution with smaller θ values

Overall execution framework
The execution procedures of EHS for sustainable ecological operation of cascade hydropower reservoirs in river ecosystem are summarized as below: Step 1: Define the basic computational parameters and then use the method in section 4.2.1 to randomly generate the initial swarm in the problem space.
Step 2: Use the constraint handling method section 4.2.2 to evaluate the fitness values of all the solutions in the current swarm.
Step 3: Update the personal best-known position for each solution and the global best-known positions for the swarm. Then, the EHS method in section 2.2 is used to update the positions of all the solutions.
Step 4: If the terminal condition is not met, go to step 2 for the next cycle; otherwise, stop the iterations and then the best solution is treated as the final solution for the ecological operation of cascade hydropower reservoirs.

Engineering background
In this section, five hydroplants located on Wu River in southwest China are used to verify the EHS feasibility in reservoir operation, including Hongjiadu (HJD), Dongfeng (DF), Suofengying (SFY), Wujiangdu (WJD) and Goupitan (GPT). As one of China's largest hydropower bases, Wu River with abundant resources and advantageous development condition plays an irreplaceable role in guaranteeing the sustainable economic and social development of Guizhou province. To reduce the potential negative effect of future uncertain climate change, operators in Wu River are paying increasing attention to the ecological operation of cascade reservoirs in recent years. Then, it is important to develop scientific optimization models and solution tools to find the scheduling schemes of cascade reservoirs under the changing environments.
To satisfy practical necessity of the Wu River, three different kinds ecological flows (minimum, suitable and ideal) for each reservoir are obtained by interpolating in the hydrological frequency curve. Three ecological flows rise at first and then descend in the scheduling horizon while the maximum or minimum values are achieved at the flood seasons ( figure 4). Several evolutionary methods (DE, PSO, GSA and EHS) are introduced for comparison, while the entire scheduling horizon is set as 1 year with 12 identical periods (months).

Case study 1 5.2.1. Statistical results analysis
The data in table 4 show that the developed method outperforms several control methods with respect to all the statistical indexes under all the cases. For instance, in the minimum ecological flow case, EHS betters PSO, DE and GSA with about 5.7%, 31.2% and 21.7% improvements in the best objective value when the frequencies of local runoff are set 95%; in the suitable ecological flow case, the mean objective value obtained by EHS with 75% local runoff is improved by about 6.4%, 17.8% and 7.1% compared with PSO, DE and GSA; in the ideal ecological flow case, the standard deviation of the objective value obtained by EHS with 85% runoff is reduced by about 91.2%, 71.7% and 90.9% as compared with PSO, DE and GSA. Besides, the execution time of EHS in different cases is close to 4 s and smaller than other methods, demonstrating superior search efficiency and convergence rate. As a result, EHS proves to be an effective tool for reservoir operation with the goal of maintaining river ecosystem health.

Box plot analysis
From figure 5, it can be found that in the minimum ecological flow cases, the DE method fails to produce satisfying solutions while three other methods (EHS, PSO and GSA) can find the best scheduling schemes without ecological flow shortages, demonstrating the importance of modified strategies; as the total runoff volume decreases, three control methods (DE, PSO and GSA) fall into local optima while EHS can provide solutions with less ecological shortage, demonstrating the superiority of EHS. Besides, the objective values obtained by EHS are much concentrated than three other methods. Thus, the robustness of EHS in reservoir operation is successfully proved.

Scheduling process analysis
For each reservoir, the total ecological deficient grows rapidly with the increase of the required ecological flows; while the maximum ecological shortage tends to occur at the flood season where the required ecological flows are relatively large (figure 6). From the above analysis, it can be concluded that EHS can yield reasonable solutions for reservoir operation problem.

Case study 2 5.3.1. Statistical results analysis
In this section, EHS is employed to resolve the reservoir ecological operation model where the typical runoff is chosen as the input information. Table 5 gives the statistical results of several methods in different runoff cases, where the goal is to satisfy the ideal ecological water requirements as far as possible. It can be seen that the objective values obtained by EHS are better than other methods in all the cases, demonstrating its feasibility in engineering application. For instance, the mean objective value of EHS is improved by about 23.9%, 28.7% and 26.8% compared with PSO, DE and GSA in the 1998-1999 typical runoff. Thus, the superiority of EHS in guaranteeing the scheme's quality under different operational environmental is proved again.

Box plot analysis
The solutions of three control methods vary in a relatively wide range in different scenarios, exhibiting unsatisfying stability in the reservoir ecological operation problem ( figure 7). Besides, the solutions by EHS nearly become a straight line, much more concentrated than three control methods. Thus, EHS can produce stable and reliable scheduling schemes for the multi-constraint reservoir operation optimization problem.

Scheduling process analysis
From figure 8, the flowing interesting phenomenon can be deduced: (a) generally, the ecological shortages of the Wu hydropower system vary with the total inflow of all the reservoirs. In other words, the larger the runoff volume, the smaller the ecological shortages. Then, river ecosystem will be destroyed gradually when the ecological water requirement is not well met for a relatively long time. (b) Due to the limited water inflow and high ecological requirement, the middle reservoirs (like SFY and DF) are more likely to face the serious water scarcity, especially during flood season. This case implies that operators should pay attention to the total discharger of cascade reservoirs in the scheduling horizon. Therefore, the feasibility of EHS in the ecological operation of cascade reservoirs problem is fully proved.

Discussions
For the above simulations, it can be clearly seen that the EHS approach can produce satisfying results in numerical experiments and engineering applications. The reasons for the superiority of the proposed method are contributed to the dynamic combination of the classical HS optimizer and two effective modified strategies. In the standard HS method, the relatively simple search module is used to search in the area around the elite solutions and incrementally update the positions of all the harmony vectors.
Although this mode has the merits of easy implementation and flexibility, all the solution's elements will gradually become similar as the iteration proceeds. As a result, the swarm easily suffers from local convergence and then falls into local optima in the last round with a high probability. On the other hand, EHS makes use of two effective strategies to balance global exploration and local exploitation. Specially, the adaptive parameter adjustment strategy can effectively improve the global search ability of the swarm, while the elite-learning evolutionary mode can sharply enhance the convergence rate of the swarm. Besides, the constraint handling tool can effectively increase the probability of yielding feasible solution by dynamically adjusting the scheduling process of the target solution. With the aid of several hybrid strategies, the proposed method possesses the merits of low computational burden, high execution efficiency and gradient information avoidance. Thus, the proposed method is able to obtain high-quality solutions for the complex global optimization problems. Meanwhile, it should be pointed out that as the population diversity decreases gradually, EHS may be trapped in local optima when used to solve other engineering optimization problems. Thus, to improve the performance of EHS in practice, the future research works can be roughly directed into three aspects: the first is to find more effective search modes helping the population explore in the promising areas; the second is to develop the parallelization optimization modules into the iterative calculation process to reduce the execution time; the last is to develop the problem-based constraint handling tool for guaranteeing the solution quality.

Conclusions
As a classical population-based evolutionary method, the HS method still suffers from the drawbacks of premature convergence as used to solve the complex multireservoir operation problem. To enhance the global search capability of the HS method, this paper develops an EHS method coupled with adaptive parameter adjustment strategy and elite-learning evolutionary mode. To fully verify the performance of EHS, a comprehensive simulation test is firstly conducted on two kinds of test functions, and the results demonstrate that compared with the standard HS method, the developed method has superior performances in both solution quality and convergence speed. Then, EHS is used to solve the nonlinear and multistage multireservoir ecological operation optimization problem. The results indicate that the developed method yields competitive results compared with several existing methods in different scenarios. Thus, the proposed method is an alternative tool for solving engineering optimization problems.

Data availability statement
The data generated and/or analyzed during the current study are not publicly available for legal/ethical reasons but are available from the corresponding author on reasonable request.