Peak Operation Problem Solving for Hydropower Reservoirs by Elite-Guide Sine Cosine Algorithm with Gaussian Local Search and Random Mutation

In recent years, growing peak pressure is posing a huge challenge for the operators of electrical power systems. As the most important clean renewable energy, hydropower is often advised as a response to the peak loads in China. Thus, a novel hybrid sine cosine algorithm (HSCA) is proposed to deal with the complex peak operation problem of cascade hydropower reservoirs. In HSCA, the elite-guide evolution strategy is embedded into the standard sine cosine algorithm to improve the convergence rate of the swarm. The Gaussian local search strategy is used to increase the diversity of the population. The random mutation operator is adopted to enhance the search capability of the individuals in the evolutionary process. The proposed method is applied to solve the complex peak operation problem of two hydropower systems. The simulations indicate that in different cases, HSCA can generate the scheduling results with higher quality than several benchmark methods. Hence, this paper provides a feasible method for the complex peak operation problem of cascade hydropower reservoirs.


Introduction
Economic development and social progress are closely linked with the energy industry. With the rapid promotion of the people's living standard, the energy consumption grows rapidly all over the world [1]. As the world's largest developing country [2], China's energy demand has also been booming over the past decades [3]. The peak operation of power systems is attracting more and more attention from scholars and engineers [4]. Among all the energy resources, the coal-fired units often produce a large amount of pollution in the form of gas and waste. Nuclear plants cannot respond to the rapid load change at peak periods due to their relatively low flexibility. The power generation of renewable energy (like wind and solar) is usually affected by the weather situation [5][6][7]. Compared with the above energies, hydropower is seen as one of the most important environment-friendly and pollution-free energy sources since it has the advantages of low operating cost and multiple comprehensive benefits (like flood control, power generation, water supply, shipping and irrigation) [8][9][10]. Meanwhile, the working conditions of hydropower plants can rapidly change in a very short period of time [11][12][13].
Thus, hydropower is playing an irreplaceable role in the peak operation of electric power systems, especially in China.
Generally, the peak operation of hydropower reservoirs aims at finding the optimal scheduling solution that can effectively reduce the maximum loads of power systems. In such a way, a smooth residual load is left for other kinds of power plants (like thermal generators) so that the frequent change of the low-efficiency generating units can be sharply reduced, which will cut down the total operational cost of power systems in the long run [14]. In practice, a variety of constraints should be considered in the peak operation of cascade hydropower reservoirs [15]. Then, a slight change of the decision variables at any period may lead to a huge change on the scheduling results of the hydropower system [16], making it become much more difficult to obtain a satisfying solution [17][18][19]. Over the past decades, many classical methods have been developed by scholars to determine a feasible decision-making solution for the complicated hydropower peak operation problem [20], such as linear programming [21], dynamic programming [22] and nonlinear programming [23]. Even though a great deal of success has been achieved in practice, these methods may fail to solve the non-differentiable problems due to the strong nonlinearity and tight coupling relationship in hydropower reservoirs [24][25][26]. Hence, with the merits of strong flexibility and adaptation, and the evolutionary methods inspired by biological evolution in nature are attracting more and more attention in a lot of engineering fields. The representative methods are the genetic algorithm (GA) [27], differential evolution (DE) [28] and the particle swarm optimization algorithm (PSO) [29]. Overall, evolutionary methods can generate multiple random initial solutions for optimization and find feasible solutions for the complicated optimization problems, but are often prone to defects such as premature convergence and optimization stagnation [30][31][32][33]. Thus, it is of great necessity to develop more effective methods for the peak operation of cascade hydropower reservoirs.
The sine cosine algorithm (SCA) is a brand new evolutionary algorithm inspired by the sine and cosine functions [34,35]. The key point of SCA is to iteratively explore unknown problem space through a simple but powerful mathematical structure and model [36,37]. Since its origin, the strong search ability of the SCA method has been verified in many experiments [38]. However, there are few reports about using SCA to deal with the hydropower peak operation problem. Thus, in order to fill the research gap, SCA is introduced to solve the reservoir operation problem for the first time. Unfortunately, the SCA algorithm is not a panacea that can be applied to all the problems due to the inherent premature convergence in intelligence algorithms [39,40]. In order to enhance the SCA performance, a novel elite-guide sine cosine algorithm based on Gaussian local search and a random mutation operator is proposed in this research, which is named hybrid sine cosine algorithm (HSCA). Then, several benchmark functions are used to verify the feasibility of the HSCA method in the global optimization problems, while two hydropower systems are employed to further test the performances of the HSCA method in the complex hydropower peak operation problem. The results show that mixing strategies can effectively improve the local exploration and global search abilities of the standard SCA method.
Finally, to better understand this study, the key contributions are summarized as below. The SCA method is employed to deal with the hydropower peak operation problem, which may be the first research report so far. A novel HSCA method based on the elite-guide evolution strategy, Gaussian local search strategy and random mutation strategy is proposed to enhance the SCA performance, which can produce satisfying optimization results for the complex problems. A practical constraint handling strategy is proposed to guarantee the solution quality in the evolutionary process. To sum up, this paper has the potential to provide certain technical reference for engineering problems in other places in the world.

Hybrid Sine Cosine Algorithm (HSCA)
In this section, the optimization problem and the detailed information of the developed HSCA method are given, including the elite-guide evolution strategy, Gaussian local search strategy and random mutation strategy.

Optimization Problem
In practice, global optimization widely exists in a variety of engineering fields. Without loss of generality, it is assumed that the goal of the optimization problem is to seek out the best decision variable vector X to minimize the objective function while satisfying a set of complicated equality and inequality constraints, which can be described as: where f (X) is the objective function, g y (X) ≤ 0 is the yth inequality constraint. h m (X) = 0 is the mth equality constraint. Y and M are the number of inequality and equality constraints. X and X are the upper and lower limits of the decision variable vector X.
When the evolutionary algorithms are employed to solve the problem presented in Equation (1), an improved solution will be generated on the basis of the newly obtained solutions. In the iterative process, the fitness value of solution X k i can be obtained by Equation (2) where the total violation value is added into the objective function: where X k i is the ith solution at the kth iteration, F(X k i ) and f (X k i ) are the fitness value and objective function of X k i . g y (X k i ) and φ y are the yth inequality constraint value and its penalty coefficient. h m (X k i ) and λ m denote the mth equality constraint value and its penalty coefficient. It is to be mentioned that if X k i is out of the boundary limit, it will be forced to be in the feasible space:

Sine Cosine Algorithm (SCA)
The standard SCA method starts from the random determination of a set of solutions in the search space. Then, the oscillating mathematical model based on disturbance rate is used to dynamically update the positions of all the individuals, and then the search regions can be explored with different sine and cosine return values [41][42][43]. During the evolutionary process, the greedy selection strategy is used to guarantee that the population can move in the better evolutionary direction, the evolution formula for the individuals in SCA is given as below: where r 1 is the variable used to balance the local exploration and global exploitation of the swarm. r 2 is a random number uniformly in the range of [0, 2π], which can control the search range of the sine and cosine functions in SCA. r 3 is a random number in the range of [0, 2], which can affect the search direction by adding the stochastic perturbation to the global optimal position. r 4 is a random number in the range of [0, 1]. k is the maximum iteration. I is the swarm size. gBest k is the position of the global optimal solution at the kth iteration, which is defined below: 2.3. Elite-Guide Evolution Strategy to Improve the Convergence Speed Generally, the global optima are the best position found by the swarm from the beginning to the current cycle. When the individuals search in the neighborhoods of the global optima, it is possible to find a solution with better performance [44][45][46]. Thus, the elite-guide evolution strategy is proposed to improve the convergence rate of the population. The combined individual balancing the current individual and the global optima is firstly generated, and then the positions of all the individuals are dynamically updated to improve the converged rate. Then, the above analysis can be expressed as: where M k i is the position of the ith combined individual at the kth iteration. a is the combination coefficient in the range of [0,1].

Gaussian Local Search Strategy to Improve the Individuals' Search Capability
As mentioned above, it is possible to find a better position around the global optimal individual in the evolution process [47][48][49][50]. Therefore, the Gaussian local search strategy in Figure 1 is used to improve the individuals' performance. Based on the Gaussian distribution, a random perturbation is employed to generate a Gaussian individual around the global optima, and then the greedy selection is used to ensure that the better solution is chosen, which can help the swarm jump out of local optima and improve the convergence rate. Then, the Gaussian local search strategy can be expressed as: where Gauss(0, 1) is the random number obeying the Gaussian distribution. B k is the position of the Gaussian individual at the kth iteration.

Elite-Guide Evolution Strategy to Improve the Convergence Speed
Generally, the global optima are the best position found by the swarm from the beginning to the current cycle. When the individuals search in the neighborhoods of the global optima, it is possible to find a solution with better performance [44][45][46]. Thus, the elite-guide evolution strategy is proposed to improve the convergence rate of the population. The combined individual balancing the current individual and the global optima is firstly generated, and then the positions of all the individuals are dynamically updated to improve the converged rate. Then, the above analysis can be expressed as: (1 )

Gaussian Local Search Strategy to Improve the Individuals' Search Capability
As mentioned above, it is possible to find a better position around the global optimal individual in the evolution process [47][48][49][50]. Therefore, the Gaussian local search strategy in Figure 1 is used to improve the individuals' performance. Based on the Gaussian distribution, a random perturbation is employed to generate a Gaussian individual around the global optima, and then the greedy selection is used to ensure that the better solution is chosen, which can help the swarm jump out of local optima and improve the convergence rate. Then, the Gaussian local search strategy can be expressed as: (1 Gauss(0,1)) g B e s t gBest gBest (10) Where Gauss(0,1) is the random number obeying the Gaussian distribution. k B is the position of the Gaussian individual at the kth iteration.

Random Mutation Strategy to Increase the Population Diversity
Similar to the other evolutionary methods, the SCA method is often trapped into local optima because the population diversity gradually decreases during the evolution process [51]. Hence, it is important to find some strategies that can help individuals search in different directions [52,53]. To achieve this goal, a novel random mutation strategy is proposed to effectively increase the diversity of the swarm. Firstly, the global optimal solution is used to generate a mutant individual for each individual so as to provide multiple search directions. Then the parent solution will be replaced by the mutant individual when certain criteria is satisfied. Figure 2 shows the sketch map of the random mutation strategy, which can be expressed as follows: where MU k i is the position of the ith mutant individual at the kth iteration. P x is the preset mutation probability, p is the index randomly selected from the set of {1, 2, · · · , I}. r 5 and r 6 are two random numbers in the range of [0, 1]. F is the average fitness value of all the individuals in the current population. Similar to the other evolutionary methods, the SCA method is often trapped into local optima because the population diversity gradually decreases during the evolution process [51]. Hence, it is important to find some strategies that can help individuals search in different directions [52,53]. To achieve this goal, a novel random mutation strategy is proposed to effectively increase the diversity of the swarm. Firstly, the global optimal solution is used to generate a mutant individual for each individual so as to provide multiple search directions. Then the parent solution will be replaced by the mutant individual when certain criteria is satisfied. Figure 2 shows the sketch map of the random mutation strategy, which can be expressed as follows: Where k i MU is the position of the ith mutant individual at the kth iteration.
x P is the preset mutation probability, p is the index randomly selected from the set of  

Execution Procedure of the Proposed HSCA Method
In the HSCA algorithm, three modified strategies are used to improve the performance of the standard SCA method in high-dimensional and multi-constraint problems, including the convergence speed, swarm diversity and solution stability. Figure 3 shows the sketch map of the HSCA algorithm while Table 1 shows the execution procedures of the HSCA method. It can be found that firstly, the SCA method can help individuals search in the problem space. Secondly, the random mutation strategy is used to generate the mutant solutions and replace the poor quality individuals with a certain probability. Finally, the Gaussian local search strategy is used to enhance the exploration capabilities as well as the quality of the global optimal position in the current swarm. Step 1: Set the computational parameters and generate the initial swarm in the problem space.
Step 2: Use Equation (2)~(3) to evaluate the fitness values of all individuals in the current swarm.
Step 3: Use Equation (6) to update the global optimal position of the swarm.
Step 4: Use Gaussian local search strategy in Equation (9)~(10) to enhance the swarm's search ability.

Execution Procedure of the Proposed HSCA Method
In the HSCA algorithm, three modified strategies are used to improve the performance of the standard SCA method in high-dimensional and multi-constraint problems, including the convergence speed, swarm diversity and solution stability. Figure 3 shows the sketch map of the HSCA algorithm while Table 1 shows the execution procedures of the HSCA method. It can be found that firstly, the SCA method can help individuals search in the problem space. Secondly, the random mutation strategy is used to generate the mutant solutions and replace the poor quality individuals with a certain probability. Finally, the Gaussian local search strategy is used to enhance the exploration capabilities as well as the quality of the global optimal position in the current swarm.
Step 5: Use the random mutation strategy in Equation (11)~(12) to generate the variant individuals.
Step 6: For each individual, use Equation (7) to update the positions after using Equations (5) and (8) to calculate the values of two intermediate variables (r1 and k i M ).
Step 8: If the maximal iteration is met, stop the iterative process and then export the global optimal position as the best solution for the target problem. Otherwise go to Step 2 for another cycle.

Numerical Experiments to Verify the Feasibility of the HSCA Method
In this section, six classical benchmark functions were used to testify the performance of HSCA, while four classic algorithms were introduced for comparative analysis.

Benchmark Functions
In this section, six famous benchmark functions in Ref [54] are employed to test the feasibility of the HSCA method. Table 2 shows the detailed information of test functions. The selected benchmark functions can be roughly classified into two different groups, unimodal functions (F1-F4) used to test the convergence ability and multimodal functions (F5-F6) used to testify the ability of jumping out of local optima.

Parameter Settings
Here, four methods are introduced for comparison and analysis, including genetic algorithm, differential evolution algorithm, particle swarm optimization algorithm and sine cosine algorithm.  Step 1: Set the computational parameters and generate the initial swarm in the problem space.
Step 2: Use Equations (2) and (3) to evaluate the fitness values of all individuals in the current swarm.
Step 3: Use Equation (6) to update the global optimal position of the swarm.
Step 4: Use Gaussian local search strategy in Equations (9) and (10) to enhance the swarm's search ability.
Step 5: Use the random mutation strategy in Equations (11) and (12) to generate the variant individuals.
Step 6: For each individual, use Equation (7) to update the positions after using Equations (5) and (8) to calculate the values of two intermediate variables (r 1 and M k i ).
Step 7: If the maximal iteration is met, stop the iterative process and then export the global optimal position as the best solution for the target problem. Otherwise go to Step 2 for another cycle.

Numerical Experiments to Verify the Feasibility of the HSCA Method
In this section, six classical benchmark functions were used to testify the performance of HSCA, while four classic algorithms were introduced for comparative analysis.

Benchmark Functions
In this section, six famous benchmark functions in Ref [54] are employed to test the feasibility of the HSCA method. Table 2 shows the detailed information of test functions. The selected benchmark functions can be roughly classified into two different groups, unimodal functions (F 1 -F 4 ) used to test the convergence ability and multimodal functions (F 5 -F 6 ) used to testify the ability of jumping out of local optima. Table 2. Details of the benchmark functions.

Function
Range Best Function

Parameter Settings
Here, four methods are introduced for comparison and analysis, including genetic algorithm, differential evolution algorithm, particle swarm optimization algorithm and sine cosine algorithm. For the sake of comparison, the swarm size and maximum iterations of four algorithms are set to 50 and 500 based on numerical experiments. Besides, the other parameters are defined as: HSCA: The combination coefficient a and mutation probability P x are set to 0.5 and 0.01. GA: The crossover probability and mutation probability are set to 0.6 and 0.2, respectively. DE: The crossover probability and scaling factor are set to be 0.8 and 0.5, respectively. PSO: The inertia weight (w) is decreased linearly from the initial value 0.9 to the final value 0.4, while the values of two learning factors (c 1 and c 2 ) are set as 2.0, respectively. The data listed in Tables 3 and 4 show that HSCA achieves better results than four other methods. Especially in the experiments of three benchmark functions (F 1 , F 2 and F 5 ), HSCA is able to find the optimal objective value. For three other functions (F 3 , F 4 and F 6 ), HSCA can still seek out the best solutions among all the methods, demonstrating the feasibility of the developed modified strategies in enhancing the SCA performance. From Figures 4 and 5, it can be seen that for three functions (F 1 , F 5 and F 6 ), the search abilities of GA and DE gradually become poor. PSO using historical memory of individuals and swarm can obtain better results than GA and DE but is still inferior to HSCA. The convergence rate of SCA outperforms three other methods. For F 3 and F 4 , DE and SCA produce satisfying performances at the initial phase, while the advantage of PSO gradually emerges as the search process increases. Meanwhile, the convergence rate of HSCA always stands out among five algorithms regardless of the function form. Thus, it can be concluded that the modified strategies can effectively accelerate the convergence rate and search capability, and HSCA is a feasible tool for solving global optimization problems.  The data listed in Tables 3-4 show that HSCA achieves better results than four other methods. Especially in the experiments of three benchmark functions (F1, F2 and F5), HSCA is able to find the optimal objective value. For three other functions (F3, F4 and F6), HSCA can still seek out the best solutions among all the methods, demonstrating the feasibility of the developed modified strategies in enhancing the SCA performance. From Figures 4-5, it can be seen that for three functions (F1, F5 and F6), the search abilities of GA and DE gradually become poor. PSO using historical memory of individuals and swarm can obtain better results than GA and DE but is still inferior to HSCA. The convergence rate of SCA outperforms three other methods. For F3 and F4, DE and SCA produce satisfying performances at the initial phase, while the advantage of PSO gradually emerges as the search process increases. Meanwhile, the convergence rate of HSCA always stands out among five algorithms regardless of the function form. Thus, it can be concluded that the modified strategies can effectively accelerate the convergence rate and search capability, and HSCA is a feasible tool for solving global optimization problems.

Function
Item

HSCA for the Peak Operation of Cascade Hydropower Reservoirs
In this section, the details of the HSCA method for the peak operation of cascade hydropower reservoirs are given, like the optimization model, swarm initialization and the constraint handling method.

Objective Function
As mentioned above, the dynamic balance between load demand and energy consumption at peak periods is posing an enormous challenge for operators in electric power systems. Thus, the goal of this paper is chosen to determine the scheduling process of all the hydropower plants to obtain a smooth residual load curve, which can be expressed as below: where T is the number of periods. L t is the load demand of the power system at tth period. N is the number of hydropower reservoirs. P n,t is the output of nth hydropower plant at tth period.

Operation Constraints
In the modeling process, the following equality and inequality constraints should be strictly satisfied to guarantee the safe operation of power system in the entire scheduling horizon [55][56][57][58].
(1) Water balance equation: where NU n is the number of upstream reservoirs directly connected to the nth reservoir. V n,t = f n (Z n,t ), I n,t , Q n,t and S n,t are the storage volume, local flow, turbine discharge and spillage of the nth reservoir at the tth period. f n (·) denotes the nonlinear relationship between forebay water level and storage volume of the nth reservoir. Z n,t is the forebay water level of the nth reservoir at the tth period.
(2) Water level constraint: where Z down n,t and Z up n,t are the minimum and maximum forebay water levels of the nth reservoir at the tth period, respectively.
where q down n,t and q up n,t are the minimum and maximum turbine discharges of the nth reservoir at the tth period, respectively.
where P down n,t and P up n,t are the minimum and maximum power outputs of the nth reservoir at the tth period, respectively. (5) Total output constraints: where L down n,t and L up n,t are the minimum and maximum total outputs of hydropower system at the tth period, respectively. (6) Initial and final water levels constraints: where Z begin n and Z end n are the initial and final water levels of the nth reservoir, respectively.

Individual Encoding and Swarm Initiation Strategies
When solving the peak operation problem of cascade hydropower reservoirs, the appropriate decision variable will be useful to accelerate the execution efficiency of the algorithm. Here, the decision variables are chosen as the total discharges of each hydropower reservoir that can be used to obtain the detailed scheduling processes (like output and water level). Then, the elements for the ith individual at the kth iteration can be described as below: X k i = [q 1,1 , q 1,2 , · · · , q 1,T · · · , q n,t , · · · , q N,1 , q N,2 , · · · , q N,T ] In the initialization stage, the decision variable q n,t of each hydropower reservoir is randomly generated in the feasible range between q down n,t and q up n,t ,which can be described as follows: where Random is the random number in the range of [0,1].

Constraint Handling Method
Due to plenty of physical constraints in practice, some individuals may fall into the infeasible zones during the evolutionary process [53,59]. To address this case, an effective constraint handling method based on penalty coefficient and heuristic modification strategies is proposed, where the variables violating the feasible zone are forced to be the boundary range, and then the constraint violation is merged into the objective function. The detailed procedure for solution X k i is as below: Step 1: Set n = 0 and start the constraint handling process.
Step 2: Set n = n + 1. If n > N, turn to Step 7, otherwise go to Step 3 to deal with the constraints of the nth hydropower reservoir.
Step 3: Define the counter c = 1. The total discharges q n,t are substituted into the water balance equation to calculate the feasible storages, which can be expressed as below: Step 4: Set c = c + 1. Calculate the difference ∆V between the obtained final storage V n,T and the predetermined terminal storage V end n . If ∆V is smaller than the minimum accuracy µ or c exceeds the constraint adjustment time c max , go to Step 6, otherwise go to Step 5.
Step 5: Use Equation (24) to calculate the modified total discharges in the feasible range. Then, Equation (22) is used to recalculate the corresponding reservoir storage volume before going back to Step 4.
Step 6: Record the outputs and constraint violations of the nth reservoir and turn to Step 2.
Step 7: Calculate the objective function and constraint violation values for the target solution, which can be described as below: where F X k i is the fitness value of solution X k i . J is the number of constraints. λ j and A j are the penalty factor and constraint violation value of the jth constraint, respectively.

Overall Execution Procedures
The execution procedure of HSCA for the peak operation of cascade hydropower reservoirs is given in Figure 6.
Step 6: Record the outputs and constraint violations of the nth reservoir and turn to Step 2.
Step 7: Calculate the objective function and constraint violation values for the target solution, which can be described as below: Where ( )

Overall Execution Procedures
The execution procedure of HSCA for the peak operation of cascade hydropower reservoirs is given in Figure 6.

Case Studies
In this section, the HSCA method was applied to determine the optimal scheduling process of two multi-reservoir hydropower systems, and the detailed procedures are given in Section 4.2. In the first hydropower system, it is assumed that there are linear relationships and no units in the characteristic curves to reduce the modeling complexity. In the second hydropower system, all the

Case Studies
In this section, the HSCA method was applied to determine the optimal scheduling process of two multi-reservoir hydropower systems, and the detailed procedures are given in Section 4.2. In the first hydropower system, it is assumed that there are linear relationships and no units in the characteristic curves to reduce the modeling complexity. In the second hydropower system, all the prototype physical constraints are used to verify the practicality of the developed methods.

Mature Multi-Reservoir System
In this section, the mature multi-reservoir system problem is used to verify the feasibility of the proposed method. The topology of the multi-reservoir system is shown in Figure 7, while basic information can be found in Tables 5 and 6. It is to be mentioned that the storage capacity changes of the 7th and 10th reservoirs are larger than those of other reservoirs. Then, five swarm-based methods mentioned in Section 3.2. are introduced for comparison. To guarantee the computational efficiency, the swarm size and maximal iterations are set as 50 and 200 for the mature multi-reservoir system. To avoid the possible solution randomness problem, four cases with different final storages in Table 7 are used for testing, and each method is independently executed 20 times.
Energies 2019, 9, x 12 of 24 In this section, the mature multi-reservoir system problem is used to verify the feasibility of the proposed method. The topology of the multi-reservoir system is shown in Figure 7, while basic information can be found in Tables 5-6. It is to be mentioned that the storage capacity changes of the 7th and 10th reservoirs are larger than those of other reservoirs. Then, five swarm-based methods mentioned in Section 3.2. are introduced for comparison. To guarantee the computational efficiency, the swarm size and maximal iterations are set as 50 and 200 for the mature multi-reservoir system. To avoid the possible solution randomness problem, four cases with different final storages in Table  7 are used for testing, and each method is independently executed 20 times.

Convergence Process Analysis
The convergence trajectories of five methods in four cases with different terminal storages are drawn in Figure 8. It can be found that as the iteration processes, the performances of DE and GA are gradually improved and the objective values of PSO and SCA change slightly, which means that the four methods are lacking search ability for the multi-constrained reservoir operation problems. On the other hand, the objective values of the HSCA method always outperform the other algorithms from beginning to end, while the entire search process presents a smooth decline trend at the early stage. Thus, it can be concluded that the proposed method has strong search efficiency in solving the high-dimensional cascade reservoir operation problem.

Convergence Process Analysis
The convergence trajectories of five methods in four cases with different terminal storages are drawn in Figure 8. It can be found that as the iteration processes, the performances of DE and GA are gradually improved and the objective values of PSO and SCA change slightly, which means that the four methods are lacking search ability for the multi-constrained reservoir operation problems. On the other hand, the objective values of the HSCA method always outperform the other algorithms from beginning to end, while the entire search process presents a smooth decline trend at the early stage. Thus, it can be concluded that the proposed method has strong search efficiency in solving the high-dimensional cascade reservoir operation problem.  Figure 9 shows the optimal objective values of five algorithms for four cases. It can be seen that the objective value of the proposed HSCA algorithm is the best in different cases. In the comparison of SCA, the objective values are reduced by about 240.64, 225.65, 195.98, and 191.66, respectively. Thus, HSCA has strong search ability in the peak operation of cascade hydropower reservoirs.  Figure 9 shows the optimal objective values of five algorithms for four cases. It can be seen that the objective value of the proposed HSCA algorithm is the best in different cases. In the comparison of SCA, the objective values are reduced by about 240. 64, 225.65, 195.98, and 191.66, respectively. Thus, HSCA has strong search ability in the peak operation of cascade hydropower reservoirs.  Figure 9. Objective value of five algorithms for the mature hydropower system in four cases. Figure 10 shows the box plots of five methods obtained in 100 independent runs for the mature hydropower system in four cases. The box plots are composed of the lower quartile, the median, the upper quartile, the upper and lower limits of the objective values, where the median is the number in the middle position in ascending order. In general, the stability of the dispersive box plot is poor. From Figure 10, it can be found that in four cases, GA and SCA show a relatively large distribution. The distributions of PSO and DE are relatively small, while HSCA shows a highly concentrated distribution. Besides, similar results between PSO and HSCA can be observed, and the reasons are analyzed as follows. Several modified strategies contribute to the best search ability of the HSCA method. In the evolutionary model of PSO, the tendency following preceding velocity and the impacts of the personal best and global best solution are simultaneously considered, making it have better performance than GA and SCA. Thus, this case fully demonstrates that although the initial swarm is generated from different random seeds, the developed HSCA method always achieves stable results in reducing the peak loads of power systems.   Figure 10 shows the box plots of five methods obtained in 100 independent runs for the mature hydropower system in four cases. The box plots are composed of the lower quartile, the median, the upper quartile, the upper and lower limits of the objective values, where the median is the number in the middle position in ascending order. In general, the stability of the dispersive box plot is poor. From Figure 10, it can be found that in four cases, GA and SCA show a relatively large distribution. The distributions of PSO and DE are relatively small, while HSCA shows a highly concentrated distribution. Besides, similar results between PSO and HSCA can be observed, and the reasons are analyzed as follows. Several modified strategies contribute to the best search ability of the HSCA method. In the evolutionary model of PSO, the tendency following preceding velocity and the impacts of the personal best and global best solution are simultaneously considered, making it have better performance than GA and SCA. Thus, this case fully demonstrates that although the initial swarm is generated from different random seeds, the developed HSCA method always achieves stable results in reducing the peak loads of power systems. analyzed as follows. Several modified strategies contribute to the best search ability of the HSCA method. In the evolutionary model of PSO, the tendency following preceding velocity and the impacts of the personal best and global best solution are simultaneously considered, making it have better performance than GA and SCA. Thus, this case fully demonstrates that although the initial swarm is generated from different random seeds, the developed HSCA method always achieves stable results in reducing the peak loads of power systems.  Figure 10. Results of different methods for the mature hydropower system in four cases. Figure 10. Results of different methods for the mature hydropower system in four cases.

Engineering Background
In this section, five reservoirs located in the mainstream of the Wu River in China are selected to further test the performance of the HSCA method, including Hongjiadu (HJD), Dongfeng (DF) and Suofengying (SFY), Wujiangdu (WJD) and Goupitan (GPT). Wu River is regarded as one of the most important hydropower bases and plays an important role in boosting the economic development of Guizhou province in southwest China. Figure 11 and Table 8 show the topological structure and basic information of the studied hydropower reservoirs, respectively. Next, the five methods in Section 3.2. are employed for comparison. To balance the computational accuracy and execution time, the swarm size and maximal iterations are set as 150 and 500. The scheduling horizon is set as 1 day with 24 identical intervals.

Engineering Background
In this section, five reservoirs located in the mainstream of the Wu River in China are selected to further test the performance of the HSCA method, including Hongjiadu (HJD), Dongfeng (DF) and Suofengying (SFY), Wujiangdu (WJD) and Goupitan (GPT). Wu River is regarded as one of the most important hydropower bases and plays an important role in boosting the economic development of Guizhou province in southwest China. Figure 11 and Table 8 show the topological structure and basic information of the studied hydropower reservoirs, respectively. Next, the five methods in Section 3.2. are employed for comparison. To balance the computational accuracy and execution time, the swarm size and maximal iterations are set as 150 and 500. The scheduling horizon is set as 1 day with 24 identical intervals.  Figure 11. Topological structure of the multi-reservoir system in the Wu River.  Table 9 shows the real-world simulation results of five algorithms in spring. From Table 9, it can be found that HSCA can get better results than four other methods in terms of all the statistical indexes. For instance, as compared with GA, DE, PSO and SCA, the peak loads in the optimization results of the HSCA method are reduced by about 1300.82, 628.00, 190.85, 1148.04 MW, respectively. Obviously, the proposed method can obtain smoother residual load curves and stand out among the compared algorithms.  Figure 11. Topological structure of the multi-reservoir system in the Wu River.  Table 9 shows the real-world simulation results of five algorithms in spring. From Table 9, it can be found that HSCA can get better results than four other methods in terms of all the statistical indexes. For instance, as compared with GA, DE, PSO and SCA, the peak loads in the optimization results of the HSCA method are reduced by about 1300.82, 628.00, 190.85, 1148.04 MW, respectively. Obviously, the proposed method can obtain smoother residual load curves and stand out among the compared algorithms.  Table 10 shows the optimization results of the HSCA algorithm in different load requirements. The data in Table 10 shows that HSCA can achieve satisfying scheduling results in three seasons. For instance, compared with the original loads, the peak values of the residual loads obtained by HSCA was improved by 26.79%, 21.68% and 25.71%, while the average standard deviations are reduced by more than 80%. Thus, this case fully demonstrates that HSCA can undertake the peak operation task of the real-world hydropower system in practice. To further show the solution rationality, Figure 12 shows the best scheduling results obtained by the HSCA method for the Wu hydropower system in four seasons. For the original load curve in four seasons, the loads in the valley periods are relatively stable at the beginning and then increase at both morning and evening peak periods. Obviously, for the optimization method, it becomes much more difficult to find the feasible scheduling solution due to the fast changing characteristic of the energy load curve. The HSCA method can collect hydropower generation during peak periods to respond to the load change and decrease the outputs during valley periods to store energy in reservoirs. In the hydropower system, HJD mainly undertakes the peak operation tasks at the early stages. The outputs of three downstream reservoirs (DF, SFY and WJD) show a relatively stable trend, while GPT plays an important role in providing generation at the later stages. In this way, the satisfactory residual loads are produced in four cases with varying load requirements. Thus, it can be concluded that HSCA has strong capability in solving the hydropower peak operation problem with multiple spatial-temporal coupling constraints.
Energies 2019, 9, x 17 of 24 be concluded that HSCA has strong capability in solving the hydropower peak operation problem with multiple spatial-temporal coupling constraints. In the HSCA method, six computational parameters may have certain effects on the optimized results. The orthogonal experiment design is used to test the role of six parameters [58]. Table 11 shows the results of the HSCA method with different parameter combinations in summer. In Table  11, the level denotes the possible value of the target parameter (the swarm sizes in level 1~3 are 100, 150 and 200), while the objective value denotes the average value of the standard deviation of the HSCA method in five independent runs. It can be found that the maximum iterations and the constraint adjustment time have the greatest influence on the final result. The reasons lie in that two parameters with larger values can effectively help the swarm produce feasible solutions with better objective values during the evolutionary process. Thus, the computation parameters have different effects on the objective value obtained by the HSCA method, and it is of great significance to carefully choose the parameter combination based on the characteristics of engineering problems in practice.  In the HSCA method, six computational parameters may have certain effects on the optimized results. The orthogonal experiment design is used to test the role of six parameters [58]. Table 11 shows the results of the HSCA method with different parameter combinations in summer. In Table 11, the level denotes the possible value of the target parameter (the swarm sizes in level 1~3 are 100, 150 and 200), while the objective value denotes the average value of the standard deviation of the HSCA method in five independent runs. It can be found that the maximum iterations and the constraint adjustment time have the greatest influence on the final result. The reasons lie in that two parameters with larger values can effectively help the swarm produce feasible solutions with better objective values during the evolutionary process. Thus, the computation parameters have different effects on the objective value obtained by the HSCA method, and it is of great significance to carefully choose the parameter combination based on the characteristics of engineering problems in practice.  Figure 13 shows the box plot of five algorithms in different seasons, where the results of each method are obtained in 20 independent runs where the initial swarm per run is generated from a newly obtained random seed. As showed in Figure 13, although three methods (GA, DE and SCA) exhibit satisfying stability in the peak operation of hydropower systems, the results of the HSCA method have the best performance due to the smallest fluctuation in the objective function. Thus, this case proves that from different initial swarms, HSCA can generate a stable scheduling process for high-dimensional and multi-constrained engineering problems.
Energies 2019, 9, x 18 of 24 newly obtained random seed. As showed in Figure 13, although three methods (GA, DE and SCA) exhibit satisfying stability in the peak operation of hydropower systems, the results of the HSCA method have the best performance due to the smallest fluctuation in the objective function. Thus, this case proves that from different initial swarms, HSCA can generate a stable scheduling process for high-dimensional and multi-constrained engineering problems.  Figure 13. Box plot test for the real-world cascade hydropower system in four seasons. Figure 14 shows the results of the remaining electrical loads obtained by three methods in four seasons, where the difference represents the maximum peak load minus the minimum valley load. As shown in Figure 14, there are no obvious differences in the average loads of SCA and HSCA, but the peak-valley differences and standard deviations of HSCA are much smaller than SCA. Generally, the smaller values of the peak-valley difference and standard deviation denote a smoother load series [60]. Therefore, it can be concluded that the HSCA approach has a stronger optimization ability in using the hydropower generation to reduce the peak loads of the power system.  Figure 14 shows the results of the remaining electrical loads obtained by three methods in four seasons, where the difference represents the maximum peak load minus the minimum valley load. As shown in Figure 14, there are no obvious differences in the average loads of SCA and HSCA, but the peak-valley differences and standard deviations of HSCA are much smaller than SCA. Generally, the smaller values of the peak-valley difference and standard deviation denote a smoother load series [60]. Therefore, it can be concluded that the HSCA approach has a stronger optimization ability in using the hydropower generation to reduce the peak loads of the power system. Energies 2019, 9, Table 12 shows the statistical results of different algorithms in 12 cases where the operation conditions (like local inflow and load demand) per case are different, including the peak load, peak-valley difference, standard deviation and load rate. Figure 15 shows the hydropower output and residual loads obtained by HSCA and SCA. From the data listed in Table 12, it can be found that the results obtained by the SCA algorithm in 12 cases are inferior to HSCA since the peak-valley difference and peak load are improved by about 300 MW and 700 MW, while the average load rates are increased by about 96%. This fully demonstrates that HSCA can solve the complex hydropower peak operation problem. Meanwhile, from Figure 15, it can be found that the residual loads obtained by SCA are unstable in the scheduling horizon, while the hydropower outputs of HSCA can respond to the load change in a better way and a smoother residual load curve is left for other energy sources represented by thermal plants. Thus, this case proves that the HSCA method has a strong ability to satisfy the peak loads of the power system under uncertainty.   Table 12 shows the statistical results of different algorithms in 12 cases where the operation conditions (like local inflow and load demand) per case are different, including the peak load, peak-valley difference, standard deviation and load rate. Figure 15 shows the hydropower output and residual loads obtained by HSCA and SCA. From the data listed in Table 12, it can be found that the results obtained by the SCA algorithm in 12 cases are inferior to HSCA since the peak-valley difference and peak load are improved by about 300 MW and 700 MW, while the average load rates are increased by about 96%. This fully demonstrates that HSCA can solve the complex hydropower peak operation problem. Meanwhile, from Figure 15, it can be found that the residual loads obtained by SCA are unstable in the scheduling horizon, while the hydropower outputs of HSCA can respond to the load change in a better way and a smoother residual load curve is left for other energy sources represented by thermal plants. Thus, this case proves that the HSCA method has a strong ability to satisfy the peak loads of the power system under uncertainty.

Result Discussion
From the above simulated results, it can be clearly found that the developed method has better performances than several benchmark methods (like GA, PSO and SCA). The standard GA method uses the crossover, selection and mutation operators to simulate the survival of the fittest in nature, and the parameter choices of three operators and the quality of the initial swarm have an important effect on the final solution. In the standard PSO method, each potential solution is seen as a particle with velocity and position vectors, and all the particles tend to fall into the neighborhood of the best-known solution of the swarm due to the loss of population diversity. In the SCA method, the sine and cosine functions are employed to guide the evolutionary process while the role of elitist individuals is ignored at a certain degree, resulting in the premature convergence problem. On the other hand, several modified strategies are helpful to guarantee the superior performance of the developed HSCA method. Firstly, the elite-guide strategy embedded in the updating equation of the standard method takes full use of the elitist individuals to provide multiple search directions, which can effectively enhance the convergence rate of the swarm in the evolutionary process. Secondly, the Gaussian local search strategy is used to generate a temporary solution for the global best-known position, which can dynamically update the shared knowledge of the swarm and then improve the search ability of all the individuals. Thirdly, the random mutation strategy is used to provide the chance to improve the solution quality for each individual, which can increase the diversity of the swarm and avoid falling into local optima. Finally, the heuristic constraint handling method can iteratively adjust the scheduling processes of all the hydropower reservoirs from upstream to downstream, which can alleviate the negative effects of the infeasible solutions. In this way, the HSCA algorithm is able to produce solutions with higher quality in different cases as used to handle the peak operation of cascade hydropower reservoirs. Besides, it can be implied that the solution quality of the evolutionary algorithms can be improved by incorporating multiple beneficial strategies in the iterative search process [56,61,62], which can provide technical reference for similar engineering problems.

Conclusions
In recent years, the peak operation of cascade hydropower reservoirs is posing a great challenge for operators in China's electrical power system. As a brand new evolutionary method, the sine cosine algorithm (SCA) has not been applied in the hydropower operation field. In order to fill this gap, the performance of the SCA method is tested in the studied engineering problem for the first time, but it was found that the standard SCA method still suffers from the premature convergence problem. Thus, for the sake of improving the SCA performance, this research develops a hybrid sine cosine algorithm (HSCA) linking the merits of several modified strategies into SCA. The HSCA method is applied to solve several test functions and the hydropower peak operation problems. The following conclusions can be deduced from the experimental results: (1) due to the premature convergence problem, the standard SCA, GA and PSO methods have poor performances in the peak operation of cascade hydropower reservoirs, (2) the HSCA method can obtain better solutions by reasonably allocating hydropower generation in the scheduling horizon, where the elite-guide evolution strategy, Gaussian local search strategy and random mutation operator can effectively improve the convergence speed and search ability of the swarm, (3) a practical constraint handling strategy is proposed to guarantee the solutions' quality in the evolutionary process. Thus, a feasible tool is provided to solve the peak operation of cascade hydropower reservoirs.
It should be mentioned that even though the feasibility of the developed method is verified in this study, there are still some shortcomings deserving attention. One is the inevitable premature convergence problem rooted in evolutionary algorithms. The other is the practicality of the constraint handling technique in the elaborate management of cascade hydropower reservoirs where more complex spatial-temporal coupling constraints need to be considered. Hence, in the future, the research can be carried out in two aspects. The first is to enhance the global search ability of the HSCA method by developing more effective strategies or introducing new ideas from other methods. The other is to develop more practical constraint handling tools based on the actual characteristics of the target engineering problems.