Multi-Objective Operation of Cascade Hydropower Reservoirs Using TOPSIS and Gravitational Search Algorithm with Opposition Learning and Mutation

: In this research, a novel enhanced gravitational search algorithm (EGSA) is proposed to resolve the multi-objective optimization model, considering the power generation of a hydropower enterprise and the peak operation requirement of a power system. In the proposed method, the standard gravity search algorithm (GSA) was chosen as the fundamental execution framework; the opposition learning strategy was adopted to increase the convergence speed of the swarm; the mutation search strategy was chosen to enhance the individual diversity; the elastic-ball modiﬁcation strategy was used to promote the solution feasibility. Additionally, a practical constraint handling technique was introduced to improve the quality of the obtained agents, while the technique for order preference by similarity to an ideal solution method (TOPSIS) was used for the multi-objective decision. The numerical tests of twelve benchmark functions showed that the EGSA method could produce better results than several existing evolutionary algorithms. Then, the hydropower system located on the Wu River of China was chosen to test the engineering practicality of the proposed method. The results showed that the EGSA method could obtain satisfying scheduling schemes in di ﬀ erent cases. Hence, an e ﬀ ective optimization method was provided for the multi-objective operation of hydropower system.


Introduction
With the merits of low pollution emission and high work efficiency, hydropower is seen as one of the most important renewable energy sources and usually shoulders many different kinds of management requirements, in practice [1][2][3]. For hydropower enterprises, it is preferred to make full use of the water resources to obtain the maximum economic benefit [4]; for power systems, the flexible hydropower generators are often asked to respond to the peak loads and smooth the energy demand curves [5]. Under such a background, the operation optimization of cascade hydropower reservoirs balancing a power generation and peak operation has become one of the most important tasks in both water resource systems and modern power systems [6]. However, as far as we know, there are few reports that have addressed this engineering problem, and therefore, the goal of this research was to refill this huge gap between theoretical research and engineering requirement. Without a loss of generality, it is assumed that there are N agents in the evolutionary swarm to determine the optimal solution for the minimization-optimization problem with D variables, and then the position of the ith agent at the kth iteration ( ( ) i X k for short) can be expressed as follows: x k is the position of the ith agent in the dth dimension at the kth iteration.
In the evolutionary process, the gravitational force of the nth agent acting on the ith agent in the dth dimension at the kth iteration can be expressed as follows:  R k X k X k is the Euclidian distance between the ith agent and the nth agent at the kth iteration. ϕ is a small constant used to prevent the denominator being 0. ( ) G k is the gravitational constant of the swarm at the kth iteration, which is dynamically updated by: where 0 G is the initial gravitational constant. α is the attenuation coefficient. k is the maximum number of iterations. After obtaining the latest positions of all agents, the mass of the ith agent at the kth iteration ( ( ) i m k for short) can be expressed as follows: (4) where ( ) i fit k is the fitness of the ith agent at the kth iteration. ( ) best k and ( ) worst k are the best and worst fitness values of the swarm at the kth iteration, respectively.
In order to reduce the negative influence of agents that are not conducive to the evolution of the swarm, the resultant force of the ith agent in the dth dimension at the kth iteration ( ( ) Without a loss of generality, it is assumed that there are N agents in the evolutionary swarm to determine the optimal solution for the minimization-optimization problem with D variables, and then the position of the ith agent at the kth iteration (X i (k) for short) can be expressed as follows: where x d i (k) is the position of the ith agent in the dth dimension at the kth iteration. In the evolutionary process, the gravitational force of the nth agent acting on the ith agent in the dth dimension at the kth iteration can be expressed as follows: where M i (k) and M n (k) are the gravitational masses of the ith agent and the nth agent at the kth iteration, respectively. R in (k) =||X i (k), X n (k) || 2 is the Euclidian distance between the ith agent and the nth agent at the kth iteration. ϕ is a small constant used to prevent the denominator being 0. G(k) is the gravitational constant of the swarm at the kth iteration, which is dynamically updated by: where G 0 is the initial gravitational constant. α is the attenuation coefficient. k is the maximum number of iterations. After obtaining the latest positions of all agents, the mass of the ith agent at the kth iteration (m i (k) for short) can be expressed as follows: where f it i (k) is the fitness of the ith agent at the kth iteration. best(k) and worst(k) are the best and worst fitness values of the swarm at the kth iteration, respectively.
Water 2019, 11, 2040 4 of 28 In order to reduce the negative influence of agents that are not conducive to the evolution of the swarm, the resultant force of the ith agent in the dth dimension at the kth iteration (F d i (k) for short) is obtained by: where rand j is the random number uniformly distributed in the range of [0,1]. Kbest is the number of agents with better fitness values. Based on the motion law in the classical Newtonian mechanics, the acceleration of the ith agent in the dth dimension at the kth iteration (acc d i (k) for short) can be expressed as follows: The velocity and position values of the ith agent in the dth dimension at the k+1th iteration (x d i (k + 1) and v d i (k + 1) for short) can be updated as below: where rand i is the random number uniformly distributed in the range of [0,1].

Opposition Learning Strategy to Improve the Convergence Speed of the Swarm
In practice, a large number of situations can be expressed by the opposition concept, like west-east, south-north, up-down, and left-right. Inspired by this case, opposition learning was proposed for intelligent computing. For the point x ∈ [a, b], the corresponding opposite point (x 1 for short) of x could be easily obtained by x 1 = a + b − x. Based on previous references, when the optimal objective function value was unknown, the opposite directions of the current solution could increase the probability of finding better solutions [56][57][58]. The opposition learning strategy could enlarge the search region in the three-dimensional space. Thus, the modified opposition learning strategy based on the social learning mechanism in PSO is proposed to improve the convergence speed of the swarm, which could be expressed as below: where re d i (k) is the ith opposite agent in the dth dimension at the kth iteration. Ub d and Lb d are the upper and lower limits of the dth dimension. c 1 and c 2 are two learning factors. rand is the random number uniformly distributed in the range of [0,1]. gBest d (k) is the dth dimension of the global optimal agent at the kth iteration.

Partial Mutation Strategy to Enhance the Individual Diversity
Generally, most agents at a later evolutionary stage are often attracted by the heavier individuals, leading to the loss of swarm diversity [47]. As a result, the GSA method easily falls into the local optima, and it is necessary to find some ways to help the agents search in different regions of the problem space. In some evolutionary algorithms represented by DE and GA, the elitism strategy is often employed to conserve better solutions, while the mutation strategy is used to promote the swarm diversity [59]. Inspired by this, a new partial mutation strategy is introduced to achieve the above goal. Specifically, the agents produced by the opposition learning strategy are first merged with the original parent swarm to form the offspring swarm υ; second, the first cbest agents in υ directly enter the next cycle, while the left agents in υ are used to generate the mutated individuals using Equation (9). In this way, some elite agents are maintained, while the search directions of other agents are promoted, which effectively enhance the individual diversity.
where B d i (k) is the position of the ith mutated agent in the dth dimension at the kth iteration. pBest d i (k) is the best-known position of the ith agent in the dth dimension at the kth iteration. λ is the number randomly selected from the index set {1, 2, · · · , N}. r 1 is the random number uniformly distributed in the range of [−0.5, 0.5].

Elastic-Ball Modification Strategy to Promote Solution Feasibility
During the random search process of the algorithm, the newly-obtained agent falls out of the feasible space with a certain probability. In the conventional method, the infeasible agents are often forced to be equal to the boundary value. As a result, the number of the agents with boundary values increase gradually with the increasing number of iterations, which might produce a certain negative influence on the evolution of the swarm. To alleviate this problem, an improved elastic-ball strategy was proposed to modify the infeasible agents [60]. Specifically, the infeasible agents were first modified to the feasible position using Equations (10) and (11), and if the newly-obtained agent was still infeasible, the corresponding position would be randomly initialized in the feasible space.
where r 2 is the random number uniformly distributed in the range of [0,1].

Execution Procedure of the Proposed EGSA Method
In the EGSA method, three modified strategies (the opposition learning strategy, mutation search strategy, and the elastic-ball modification strategy) were introduced to enhance the global search ability of the standard GSA method in the complex nonlinear constrained optimization problems. Specially, the GSA module could well guarantee the search performance of the swarm; the opposition learning strategy could increase social interaction ability of the agents by exploring the reverse areas; the mutation search strategy was used to balance the convergence speed and the swarm diversity; while the elastic-ball modification strategy could effectively guarantee the feasibility of the newly-generated agents. By combining the advantages of the above strategies, the EGSA method had a stronger ability to enhance the local exploration and global search, in comparison with the conventional GSA method, which would help alleviate the premature convergence problem. Then, the brief execution procedure of the EGSA method was summarized as below: Step 1: Set the values of the computational parameters and then randomly generate the initial swarm in the problem space.
Step 2: Calculate the fitness values of all the agents in the current population, and then update the personal best-known of each agent and the global best-known agent of the swarm.
Step 3: Calculate the correlated variables (like the gravitational coefficient, mass, and acceleration) to update the velocity and position values of all the agents.
Step 4: Execute the opposition learning strategy to increase the convergence speed of the swarm.
Step 5: Execute the partial mutation search strategy to enhance the individual diversity.
Step 6: Execute the elastic-ball modification strategy to promote the solution feasibility.
Step 7: Repeat Step 2-6 until the stopping criterion is met, and then the global optimal position is regarded as the final solution of the optimization problem.

Benchmark Functions
To verify the performance of the EGSA algorithm, 12 classic benchmark functions are chosen for numerical experiments. Table 1 shows the details of the selected functions, where D is the number of variables; range is the search scope for each variable; and f min is the optimal objective value in theory. For all test functions, the optimal value of F 8, −418.9 × D, varies with the dimension; while the optimal values for other functions are 0. Meanwhile, the benchmark functions are divided into two categories-unimodal functions (F 1 -F 8 ) with one global optimum, and multimodal functions (F 9 -F 12 ) with multiple local optimums. The unimodal functions are used to test the convergence speed of the algorithm, while the multimodal function can test the ability of jumping out of the local optimum, which can fully verify the performance of the evolutionary methods. Table 1. Detailed information of the twelve benchmark functions.

Parameters Settings
For the sake of fair comparison, the results of 11 famous evolutionary algorithms are introduced in this section, including differential evolution (DE), particle swarm optimization (PSO), sine cosine algorithm (SCA), gravitational search algorithm (GSA), ant lion optimizer (ALO) [61], cuckoo search algorithm (CS) [62], modified cuckoo search algorithm (MCS) [63], lightning search algorithm (LSA) [64], grey wolf optimizer (GWO) [65], firefly algorithm (FA) [66], and whale optimization algorithm (WOA) [67]. The results of these five methods (DE, PSO, SCA, GSA, and EGSA) were developed in the JAVA procedures, while the results of the other methods were taken from the corresponding literature. For the five developed methods, the swarm size and maximum iteration were set as 50 and 1000; while the other parameters were set as: DE: The scaling factor was set as 0.5 and the crossover probability was set as 0.6. PSO: The inertia weight w was linearly decremented from 0.9 to 0.3; while two learning factors (c 1 and c 2 ) were set as 2.0, respectively. SCA: The computational constant a was set as 2.0. GSA: The attenuation coefficient was set as 20; the initial gravitational constant was set as 100. EGSA: The attenuation coefficient was set to 20; the initial gravitational constant was set to 100; while the value of cbest was set to 0.7.

Result Comparison in Multiple Runs
For the sake of alleviating the adverse influence of random initial seeds, all approaches were independently executed 30 times. Table 2 shows the average (Ave.) and standard deviation (Std.) of 12 algorithms for all the test functions, where the number of variables was set as 30. It could be seen that for all the functions, EGSA was always superior to the GSA method in terms of both the average and standard deviation; as compared with other methods, the EGSA method could obtain satisfying results in most functions. For instance, EGSA outperformed ALO in the average of all test functions except for F 5 , the WOA performance was tied with EGSA in F 9 and defeated by EGSA in the other functions. To sum up, the proposed method could obtain better results than the other evolutionary algorithms in the 12 test functions.

Box and Whisker Test
In this section, the box and whisker test was employed to demonstrate the objective distribution of the solutions since it can show the abundant information (including minimum, second and third quartile, median, and maximum) on the studied data. Figure 2 shows the detailed data obtained by the three methods in 30 independent runs. It could be clearly observed that compared to the other methods, the proposed method had a smaller distribution and better values in all indices. For instance, the standard GSA algorithm had different degrees of outliers on F 4 -F 5 , which indicated that this method easily fell into the state of premature convergence in the search process; the SCA algorithm exhibited outliers in most functions but had a distribution of relatively discrete solutions in F 3 , which indicated that SCA was not ideal in terms of robustness. Additionally, the EGSA method had a more concentrated solution distribution and fewer outliers in all functions, demonstrating its satisfying robustness and search ability. Thus, the performances of the EGSA method were superior to the comparative methods in the 12 functions. The scaling factor was set as 0.5 and the crossover probability was set as 0.6. PSO: The inertia weight w was linearly decremented from 0.9 to 0.3; while two learning factors (c1 and c2) were set as 2.0, respectively. SCA: The computational constant a was set as 2.0. GSA: The attenuation coefficient was set as 20; the initial gravitational constant was set as 100. EGSA: The attenuation coefficient was set to 20; the initial gravitational constant was set to 100; while the value of cbest was set to 0.7.

Result Comparison in Multiple Runs
For the sake of alleviating the adverse influence of random initial seeds, all approaches were independently executed 30 times. Table 2 shows the average (Ave.) and standard deviation (Std.) of 12 algorithms for all the test functions, where the number of variables was set as 30. It could be seen that for all the functions, EGSA was always superior to the GSA method in terms of both the average and standard deviation; as compared with other methods, the EGSA method could obtain satisfying results in most functions. For instance, EGSA outperformed ALO in the average of all test functions except for F5, the WOA performance was tied with EGSA in F9 and defeated by EGSA in the other functions. To sum up, the proposed method could obtain better results than the other evolutionary algorithms in the 12 test functions.

Box and Whisker Test
In this section, the box and whisker test was employed to demonstrate the objective distribution of the solutions since it can show the abundant information (including minimum, second and third quartile, median, and maximum) on the studied data. Figure 2 shows the detailed data obtained by the three methods in 30 independent runs. It could be clearly observed that compared to the other methods, the proposed method had a smaller distribution and better values in all indices. For instance, the standard GSA algorithm had different degrees of outliers on F4-F5, which indicated that this method easily fell into the state of premature convergence in the search process; the SCA algorithm exhibited outliers in most functions but had a distribution of relatively discrete solutions in F3, which indicated that SCA was not ideal in terms of robustness. Additionally, the EGSA method had a more concentrated solution distribution and fewer outliers in all functions, demonstrating its satisfying robustness and search ability. Thus, the performances of the EGSA method were superior to the comparative methods in the 12 functions.

Wilcoxon Nonparametric Test
To achieve statistically meaningful conclusions, the Wilcoxon nonparametric test was employed to test the significance level of various methods. Table 3 shows the statistical results of various methods. If the results of the EGSA was better than the control method, it was recorded as a win; if two methods were tied, it was recorded as a tie; otherwise, it was recorded as a loss. From Table 3, it can be seen that the proposed method could produce better solutions as compared to the other methods, due to the larger number of the 'Win' symbol. Then, the multiple-problem statistical results obtained by the Wilcoxon nonparametric test with a significance level of 0.05 are given in Table 4, where the mean objective value obtained by each method is chosen as the data sample. From Table  4, it can be found that the proposed method could achieve higher R + than Rin all comparisons, while all values of p were smaller than 0.05. This case proved that the EGSA method outperformed the other methods in a statistical manner.

Wilcoxon Nonparametric Test
To achieve statistically meaningful conclusions, the Wilcoxon nonparametric test was employed to test the significance level of various methods. Table 3 shows the statistical results of various methods. If the results of the EGSA was better than the control method, it was recorded as a win; if two methods were tied, it was recorded as a tie; otherwise, it was recorded as a loss. From Table 3, it can be seen that the proposed method could produce better solutions as compared to the other methods, due to the larger number of the 'Win' symbol. Then, the multiple-problem statistical results obtained by the Wilcoxon nonparametric test with a significance level of 0.05 are given in Table 4, where the mean objective value obtained by each method is chosen as the data sample. From Table 4, it can be found that the proposed method could achieve higher R + than Rin all comparisons, while all values of p were smaller than 0.05. This case proved that the EGSA method outperformed the other methods in a statistical manner.     Figure 3 shows the convergence trajectories of five algorithms for 12 functions with 30 variables. It can be seen that the proposed method could constantly improve the solution's quality from the beginning to the end, and could find the best solutions in the end. Taking the functions F 1 -F 4 as examples, three methods (GSA, PSO, and DE) quickly fail into local optima because their best-so-far solutions change slightly in the evolutionary process; both EGSA and SCA converge to a smaller objective at the initial stage, but EGSA can quickly find better solutions at iteration 800, while SCA fails to make it. Additionally, for multimodal functions, EGSAs achieve the global optimal solution in both F 9 and F 11 , and better results than other methods in F 8 , F 10 , and F 12 . Thus, the above analysis fully proves that the presented method has a superior convergence speed and global search ability.  Figure 3 shows the convergence trajectories of five algorithms for 12 functions with 30 variables. It can be seen that the proposed method could constantly improve the solution's quality from the beginning to the end, and could find the best solutions in the end. Taking the functions F1-F4 as examples, three methods (GSA, PSO, and DE) quickly fail into local optima because their best-so-far solutions change slightly in the evolutionary process; both EGSA and SCA converge to a smaller objective at the initial stage, but EGSA can quickly find better solutions at iteration 800, while SCA fails to make it. Additionally, for multimodal functions, EGSAs achieve the global optimal solution in both F9 and F11, and better results than other methods in F8, F10, and F12. Thus, the above analysis fully proves that the presented method has a superior convergence speed and global search ability.

Objective Functions
With the advantages of low pollutant emission and high operational efficiency, hydropower now plays an increasingly important role in the energy system throughout the world [68][69][70]. As a result, the cascade hydropower reservoirs are often responsible for a variety of operational requirements. Considering the practical requirements of generation enterprises and electrical power systems in China, the generation benefit and peak operations were chosen as the focus of this paper. (

1) Power Generation
For almost all hydropower enterprises, the economic benefit is an important indicator of performance check under the market environment [71][72][73]. Generally, the same amount of available water resource can produce numerous scheduling schemes with different generation benefits. Hence, the first objective was chosen to find the best scheme that could maximize the total power generation of all hydropower reservoirs, which could be expressed as below: where E is the power generation of all hydropower reservoirs. T is the number of periods. N is the number of hydropower reservoirs. P n,t is the output of the nth hydroplant at the tth period. ∆ t is the total hours of the tth period.
(2) Peak Operation In modern power systems, the hydropower system is often asked to respond to the peak loads of power grid to produce a relatively smooth load series for other kinds of energy systems (like thermal, solar, and wind) [74][75][76]. In this way, the total operational cost of the electrical power system in the long run can be sharply reduced [77][78][79]. Hence, the second objective function was chosen to minimize the mean square deviation of the residual load curve obtained by subtracting all power outputs of hydropower reservoirs from the original energy load curve, which could be expressed as: where Load t is the load demand of the power system at the tth period.

Physical Constraints
(1) Water Balance Equation: where V n,t , I n,t , Q n,t , S n,t are the storage volume, local inflow, turbine discharge and spillage of the nth hydroplant at the tth period, respectively. Num is the number of upstream hydroplants directly connected to the nth hydroplant.
(2) Reservoir Storage Constraint: where V n,t and V n,t denote the maximum and minimum storage volumes of the nth hydroplant at the tth period, respectively.
(3) Turbine Discharge Constraint: where Q n,t and Q n,t denote the maximum and minimum turbine discharges of the nth hydroplant at the tth period, respectively.
(4) Total Discharge Constraint: q n,t ≤ q n,t ≤ q n,t where q n,t and q n,t denote the maximum and minimum total discharges of the nth hydroplant at the tth period, respectively. (5) Power Output Constraint: P n,t ≤ P n,t ≤ P n,t where P n,t and P n,t denote the maximum and minimum power outputs of the nth hydroplant at the tth period, respectively. (6) Initial and Final Water Levels Constraint: where Z beg n and Z end n represent the initial and final water level of the nth hydroplant, respectively.

Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS)
As mentioned above, the multi-objective operation of cascade hydropower reservoirs is a typical multiple-criteria decision-making problem. Thus, it was necessary to find some tools to transform the original multi-objective problem into the single-objective optimization problem that can be addressed by EGSA. After a comprehensive comparison, TOPSIS was introduced to achieve this goal. In TOPSIS, the best scheme had the smallest distance from the positive ideal scheme and the largest distance from the negative ideal scheme, respectively. Then, the procedure of TOPSIS was give as below: Step 1: Create an initial decision evaluation matrix with m schemes and J attributes. For the target problem, each scheme is composed of two attributes-the power generation E and the peak operation F. Then, the decision-making matrix A can be expressed as follows: where a i,j is the original value of the jth attribute in the ith scheme.
Step 2: Obtain the modified decision evaluation matrix to reduce the potential impact of various attribute properties (like measurement unit and tendency); given below: where a max Step 3: Obtain the normalized decision-making matrix to reduce the possible numerical problem caused by the feature differences (like units or magnitude), which could be expressed as below: where r ij is the normalized value of the jth attribute in the ith scheme.
Step 4: Set the weight vector w = w 1 , · · · , w j , · · · , w J , where w j is the weight of the jth attribute and there is J j=1 w j = 1. Then, calculate the weighted normalized decision-making matrix V by: Step 5: Determine the positive ideal scheme C + and negative ideal scheme C − by: Step 6: Calculate the Euclidean distances between the decision scheme and the positive (negative) ideal scheme, which could be expressed as below: where D + i or D − i is the distance between the ith scheme and the position (or negative) ideal scheme.
Step 7: Calculate the closeness of each decision plan to the negative ideal scheme by Step 8: Use the obtained closeness values to sort all the schemes, and the scheme possessing the maximum closeness will be treated as the optimal decision scheme.

Individual Structure and Swarm Initialization
To enhance the execution efficiency of reservoir operation, the total discharge of the hydropower reservoir is selected as the decision variable for encoding, and then the ith agent at the kth iteration X i (k) can be expressed as below: During the initialization phase, the elements of the solution X i (k) are randomly generated in the preset of the total discharge scopes of all hydroplants, which could be expressed as follows: where random is the random number uniformly distributed in the range of [0,1].

Heuristic Constraint Handling Method
During the evolutionary process, it is possible that some agents violate the equality or inequality constraints imposed on the hydropower system. In order to effectively modify the infeasible agents, this section presents a heuristic constraint handling method that tries to equally allocate the violated final storage volume to the turbine discharges of all adjusting periods. The detailed procedures for the solution X i (k) are given as below: Step 1: Set n = 1 and determine the necessary parameters for the constraint processing.
Step 2: Set the counter L = 1, and then use the water balance equation to calculate the storage of the nth hydroplant at the tth period by Step 3: Set L = L + 1, and then calculate the violation value η n of the final storage by Equation (33). Then, if η n is smaller than the accuracy ε or L is larger than the maximum iteration, go to Step 5; otherwise, turn to Step 4.
η n = f n (Z n,T ) − f n Z end n (33) where f n (·) is the nonlinear stage-storage curve of the nth hydropower reservoir.
Step 4: Use Equation (34) to obtain the modified total discharges of the target reservoir, and then calculate the corresponding reservoir storage volume by Equation (32) before turning to Step 3.
Step 5: Set n = n+1. If n ≤ N, turn to Step 2 for the new hydroplant; otherwise, stop the iteration.

Calculation of the Modified Objective Values
Generally, the agents modified by the above heuristic constraint handling procedures would satisfy most of the physical constraints. Nevertheless, the agents might still violate some constraints due to reasons like unreasonable operation trajectory and limited adjusted times. To eliminate the negative effect of the infeasible agents, the value of the constraint violation involved in the solution X i (k) (Viol[X i (k)] for short) was obtained by Equation (35), and then merged into two original objective values (E and F) by Equation (36) to obtain the modified objective values (E ' and F ' ).
where g l 1 (·) ≤ 0 is the l 1 th inequality constraint; e l 2 (·) ≤ 0 is the l 2 th equality constraint; c l 1 and c l 2 are the penalty coefficients for the relevant inequality or equality constraint. L g and L e denote the number of inequality and equality constraints, respectively.

Execution Procedures of the EGSA Method for the Target Problem
The flowchart of the developed EGSA method for solving multi-objective operation of cascade hydropower reservoirs is given in Figure 4.

Execution Procedures of the EGSA Method for the Target Problem
The flowchart of the developed EGSA method for solving multi-objective operation of cascade hydropower reservoirs is given in Figure 4.

Engineering Background
In this section, five reservoirs of the Wu hydropower system in Southwest China were chosen to test the performance of the proposed method, including Hongjiadu (HJD), Dongfeng (DF), Suofengying (SFY), Wujiangdu (WJD), and Goupitan (GPT). Since being put into operation, the Wu hydropower system has played a vitally important role in the environment protection and economic development of the Guizhou province. Figure 5 shows the topological structure of the selected hydropower reservoirs, while the typical load demand curves of the four seasons are shown in Figure  6. From Figures 5 and 6, it could be observed that in the Wu hydropower system, there were tight hydraulic and electrical connections from upstream to downstream reservoirs, while four load curves with different features further increased the optimization difficulty.

Engineering Background
In this section, five reservoirs of the Wu hydropower system in Southwest China were chosen to test the performance of the proposed method, including Hongjiadu (HJD), Dongfeng (DF), Suofengying (SFY), Wujiangdu (WJD), and Goupitan (GPT). Since being put into operation, the Wu hydropower system has played a vitally important role in the environment protection and economic development of the Guizhou province. Figure 5 shows the topological structure of the selected hydropower reservoirs, while the typical load demand curves of the four seasons are shown in Figure 6. From Figures 5 and 6, it could be observed that in the Wu hydropower system, there were tight hydraulic and electrical connections from upstream to downstream reservoirs, while four load curves with different features further increased the optimization difficulty.

Robustness Testing of Different Evolutionary Algorithms
In this section, to verify the robustness of the developed method, four famous methods (DE, PSO, SCA, and GSA) were introduced for comparative analysis. Table 5 shows the detailed statistical results of the 5 methods in 20 independent runs for 4 cases, where the swarm size and maximum iteration were set as 50 and 500, respectively. It can be clearly seen that the solutions obtained by the EGSA method were more stable than the other methods. For instance, in Case 1, the EGSA method could make about 96.4%, 96.5%, 97.8%, and 97.7% reductions in the standard deviations of the DE, PSO, SCA, and GSA. Hence, the feasibility of the EGSA method was fully demonstrated in this case.

Robustness Testing of Different Evolutionary Algorithms
In this section, to verify the robustness of the developed method, four famous methods (DE, PSO, SCA, and GSA) were introduced for comparative analysis. Table 5 shows the detailed statistical results of the 5 methods in 20 independent runs for 4 cases, where the swarm size and maximum iteration were set as 50 and 500, respectively. It can be clearly seen that the solutions obtained by the EGSA method were more stable than the other methods. For instance, in Case 1, the EGSA method could make about 96.4%, 96.5%, 97.8%, and 97.7% reductions in the standard deviations of the DE, PSO, SCA, and GSA. Hence, the feasibility of the EGSA method was fully demonstrated in this case.

Robustness Testing of Different Evolutionary Algorithms
In this section, to verify the robustness of the developed method, four famous methods (DE, PSO, SCA, and GSA) were introduced for comparative analysis. Table 5 shows the detailed statistical results of the 5 methods in 20 independent runs for 4 cases, where the swarm size and maximum iteration were set as 50 and 500, respectively. It can be clearly seen that the solutions obtained by the EGSA method were more stable than the other methods. For instance, in Case 1, the EGSA method could make about 96.4%, 96.5%, 97.8%, and 97.7% reductions in the standard deviations of the DE, PSO, SCA, and GSA. Hence, the feasibility of the EGSA method was fully demonstrated in this case.  Figure 7 draws the box and whisker test of the best-so-far solutions obtained by 5 algorithms in 4 cases. It can be observed that the performances of four control methods in the power generation of hydropower system are relatively stable but still inferior to the EGSA method possessing the smallest variation ranges in the final solutions. Therefore, the proposed method proves to be an effective tool to deal with the power generation of a hydropower system.  Figure 7 draws the box and whisker test of the best-so-far solutions obtained by 5 algorithms in 4 cases. It can be observed that the performances of four control methods in the power generation of hydropower system are relatively stable but still inferior to the EGSA method possessing the smallest variation ranges in the final solutions. Therefore, the proposed method proves to be an effective tool to deal with the power generation of a hydropower system.  Table 6 shows the detailed results in the best solutions obtained by different methods. It can be seen that in all cases, EGSA can always find the best results among all algorithms. For instance, compared with DE, PSO, SCA, and GSA, the EGSA method could make about 8.77 × 10 4 kW·h, 3.2 × 10 4 kW·h, 8.61 × 10 4 kW·h, and 8.5 × 10 4 kW·h improvements in power generation in Case 1, respectively. Additionally, among all hydroplants, the GPT hydroplant with a huge installed capacity produced the largest proportion of power generation in all solutions, which was in line with the actual operation situation of the hydropower system [80][81][82][83]. Thus, the above analysis demonstrated the effectiveness and rationality of the scheduling process obtained by the proposed method.   Table 6 shows the detailed results in the best solutions obtained by different methods. It can be seen that in all cases, EGSA can always find the best results among all algorithms. For instance, compared with DE, PSO, SCA, and GSA, the EGSA method could make about 8.77 × 10 4 kW·h, 3.2 × 10 4 kW·h, 8.61 × 10 4 kW·h, and 8.5 × 10 4 kW·h improvements in power generation in Case 1, respectively. Additionally, among all hydroplants, the GPT hydroplant with a huge installed capacity produced the largest proportion of power generation in all solutions, which was in line with the actual operation situation of the hydropower system [80][81][82][83]. Thus, the above analysis demonstrated the effectiveness and rationality of the scheduling process obtained by the proposed method.   Figure 8 shows the convergence trajectories of 5 algorithms in different cases. It could be observed that in four cases, SCA failed into a local optima at the early stage, and then started to improve the quality of solutions as the iteration number reached 300, but still failed to find out the satisfying solutions at the end; due to the loss of individual diversity [84], three other methods (PSO, GSA, and DE) outperformed the SCA method but were still inferior to the proposed method from the beginning to end. Additionally, the proposed method could quickly converge to the scheduling schemes that were close to the optimal solution at the early stage, and then change slightly with an increasing number of iterations. Thus, it could be concluded that three modified strategies could effectively enhance the performance of the standard GSA method.  Figure 8 shows the convergence trajectories of 5 algorithms in different cases. It could be observed that in four cases, SCA failed into a local optima at the early stage, and then started to improve the quality of solutions as the iteration number reached 300, but still failed to find out the satisfying solutions at the end; due to the loss of individual diversity [84], three other methods (PSO, GSA, and DE) outperformed the SCA method but were still inferior to the proposed method from the beginning to end. Additionally, the proposed method could quickly converge to the scheduling schemes that were close to the optimal solution at the early stage, and then change slightly with an increasing number of iterations. Thus, it could be concluded that three modified strategies could effectively enhance the performance of the standard GSA method.

Robustness Testing of Different Evolutionary Algorithms
In this section, the proposed method was applied to deal with the peak operation of a hydropower system. Table 7 shows the statistical results of the 5 methods in 20 independent runs for

Robustness Testing of Different Evolutionary Algorithms
In this section, the proposed method was applied to deal with the peak operation of a hydropower system. Table 7 shows the statistical results of the 5 methods in 20 independent runs for 4 cases. It was observed that all of the methods could achieve the goal of reducing the peak loads of a power system due to the reductions in the statistical indices; while the proposed method had the performances among all methods. For instance, the average improvements in the range of the objective values of EGSA was about 99.7%, compared to the three other methods in the 4 cases. Hence, the feasibility of the EGSA method in addressing the peak operation of the hydropower system was proved in this case.  Figure 9 shows the box and whisker test of 5 algorithms in different cases. It can be seen that the distribution ranges of the solutions obtained by the proposed method were much smaller than the other methods, demonstrating the effectiveness of the constraint handling method and the modified strategies. Thus, it can be concluded that the EGSA method could produce stable scheduling schemes for the complex hydropower system peak operation problem.
Water 2019, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/water 4 cases. It was observed that all of the methods could achieve the goal of reducing the peak loads of a power system due to the reductions in the statistical indices; while the proposed method had the performances among all methods. For instance, the average improvements in the range of the objective values of EGSA was about 99.7%, compared to the three other methods in the 4 cases. Hence, the feasibility of the EGSA method in addressing the peak operation of the hydropower system was proved in this case.  Figure 9 shows the box and whisker test of 5 algorithms in different cases. It can be seen that the distribution ranges of the solutions obtained by the proposed method were much smaller than the other methods, demonstrating the effectiveness of the constraint handling method and the modified strategies. Thus, it can be concluded that the EGSA method could produce stable scheduling schemes for the complex hydropower system peak operation problem.  Figure 9. Box and whisker test of 5 algorithms in different cases for peak operation.

Comparison of the Optimal Results Obtained by Different Evolutionary Algorithms
To clearly show the feasibility of the EGSA method, the detailed results in the best solutions as obtained by different methods are given in Table 8. It was found that the EGSA method could achieve satisfactory results. For instance, the proposed method could bring about a 30.9% reduction in the peak values in Summer, which was obviously better than 21.0% of DE and 24.9% of PSO, 23.9% of SCA and 19.4% of GSA. Thus, this simulation clearly proved the superiority of the proposed method in responding to the peak operation requirement of the power system.  Figure 10 shows the detailed output process obtained by the different methods in four cases. It was found that there were often two peak periods (morning and evening) in the original load curves, and obvious differences in the load features (like peak times or numbers) of the four load curves, which would increase the optimization difficulty of the peak operation. The two methods (GSA and EGSA) exhibited excellent responses to the load changes by collecting the hydropower generation at the peak periods and storing the energy at valley periods. The EGSA method was superior to the GSA method due to its smoother residual load curves; besides, the outputs of all the hydroplants varied in the feasible range between the minimum output limit and the installed capacity. Thus, the EGSA method could produce feasible solutions for the peak operation of a hydropower system.  Figure 10 shows the detailed output process obtained by the different methods in four cases. It was found that there were often two peak periods (morning and evening) in the original load curves, and obvious differences in the load features (like peak times or numbers) of the four load curves, which would increase the optimization difficulty of the peak operation. The two methods (GSA and EGSA) exhibited excellent responses to the load changes by collecting the hydropower generation at the peak periods and storing the energy at valley periods. The EGSA method was superior to the GSA method due to its smoother residual load curves; besides, the outputs of all the hydroplants varied in the feasible range between the minimum output limit and the installed capacity. Thus, the EGSA method could produce feasible solutions for the peak operation of a hydropower system.  The focus of this section is on the demonstration of the feasibility of the proposed method in addressing the multi-objective optimization problems of a hydropower system. Figure 11 draws the distributions of the objective functions obtained by different methods, where each method was Water 2019, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/water increasing objective value of peak operation, which implied that there was an obvious conflict between power generation and peak operation. Additionally, the solutions of PSO and GSA were dominated by that of the EGSA method, which meant that the EGSA method had a higher probability to obtain the Pareto optimal solutions than the PSO and GSA method. Thus, the proposed method can generate the near-optimal Pareto solutions to approximate the relationship between power generation and peak shaving in practice.  Figure 12 shows the detailed results of three typical schemes obtained by the EGSA method in Spring, where Scheme 1 showed the maximum power generation, Scheme 3 had the best peak operation performance, while Scheme 2 could achieve a balance between power generation and peak operation. It could be seen that for all the hydropower reservoirs, the water levels varied in the preset range between a dead water level and a normal water level while the power outputs were smaller than the installed capacity, demonstrating the effectiveness of the heuristic constraint handling method in guaranteeing the feasibility of the solution. Meanwhile, there were obvious differences in the peak operation and power generations of the three typical schemes, which indicated that it was necessary to choose a suitable scheme based on the actual requirement. Thus, this case again proved the practicability of the EGSA method in solving the multi-objective operation problems.  Figure 12 shows the detailed results of three typical schemes obtained by the EGSA method in Spring, where Scheme 1 showed the maximum power generation, Scheme 3 had the best peak operation performance, while Scheme 2 could achieve a balance between power generation and peak operation. It could be seen that for all the hydropower reservoirs, the water levels varied in the preset range between a dead water level and a normal water level while the power outputs were smaller than the installed capacity, demonstrating the effectiveness of the heuristic constraint handling method in guaranteeing the feasibility of the solution. Meanwhile, there were obvious differences in the peak operation and power generations of the three typical schemes, which indicated that it was necessary to choose a suitable scheme based on the actual requirement. Thus, this case again proved the practicability of the EGSA method in solving the multi-objective operation problems.

Rationality Analysis of the Best Results Obtained by the Different Evolutionary Algorithms
Water 2019, 11, x; doi: FOR PEER REVIEW www.mdpi.com/journal/water increasing objective value of peak operation, which implied that there was an obvious conflict between power generation and peak operation. Additionally, the solutions of PSO and GSA were dominated by that of the EGSA method, which meant that the EGSA method had a higher probability to obtain the Pareto optimal solutions than the PSO and GSA method. Thus, the proposed method can generate the near-optimal Pareto solutions to approximate the relationship between power generation and peak shaving in practice.  Figure 12 shows the detailed results of three typical schemes obtained by the EGSA method in Spring, where Scheme 1 showed the maximum power generation, Scheme 3 had the best peak operation performance, while Scheme 2 could achieve a balance between power generation and peak operation. It could be seen that for all the hydropower reservoirs, the water levels varied in the preset range between a dead water level and a normal water level while the power outputs were smaller than the installed capacity, demonstrating the effectiveness of the heuristic constraint handling method in guaranteeing the feasibility of the solution. Meanwhile, there were obvious differences in the peak operation and power generations of the three typical schemes, which indicated that it was necessary to choose a suitable scheme based on the actual requirement. Thus, this case again proved the practicability of the EGSA method in solving the multi-objective operation problems.

Conclusions
In this research, an enhanced gravitational search algorithm coupled with the Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) was developed to deal with the complex multi-objective operation problem of cascade hydropower reservoirs balancing the power generation and peak operation benefits. The proposed method and several famous evolutionary methods were used to solve the 12 famous benchmark functions and the optimal operation of the Wu hydropower system in China. The critical findings obtained from the experimental results are given as below: (1) Due to the loss of the population diversity, the conventional GSA method suffered from severe premature convergence shortcomings. The proposed method based on the three modified strategies (opposite learning strategy, partial mutation strategy and elastic-ball modification strategy) could effectively improve the convergence speed, swarm diversity, and solution feasibility of the standard GSA method, respectively. (2) For the original complex multi-objective optimization problem, balancing power generation and peak operation requirements, the famous TOPSIS method was used to transform it into the relatively simple single-objective problem, which could help make an obvious reduction in the modeling difficulty of the multi-objective decision. (3) There was a competitive relationship between the generation benefit of the hydropower enterprise and the peak operation of the power system. In other words, an increasing value of one objective would obviously reduce another objective value. Thus, it was necessary for operators to carefully determine the scheduling schemes so as to effectively balance the practical requirements of power generation enterprises and power grid companies.

Conclusions
In this research, an enhanced gravitational search algorithm coupled with the Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) was developed to deal with the complex multi-objective operation problem of cascade hydropower reservoirs balancing the power generation and peak operation benefits. The proposed method and several famous evolutionary methods were used to solve the 12 famous benchmark functions and the optimal operation of the Wu hydropower system in China. The critical findings obtained from the experimental results are given as below: (1) Due to the loss of the population diversity, the conventional GSA method suffered from severe premature convergence shortcomings. The proposed method based on the three modified strategies (opposite learning strategy, partial mutation strategy and elastic-ball modification strategy) could effectively improve the convergence speed, swarm diversity, and solution feasibility of the standard GSA method, respectively. (2) For the original complex multi-objective optimization problem, balancing power generation and peak operation requirements, the famous TOPSIS method was used to transform it into the relatively simple single-objective problem, which could help make an obvious reduction in the modeling difficulty of the multi-objective decision. (3) There was a competitive relationship between the generation benefit of the hydropower enterprise and the peak operation of the power system. In other words, an increasing value of one objective would obviously reduce another objective value. Thus, it was necessary for operators to carefully determine the scheduling schemes so as to effectively balance the practical requirements of power generation enterprises and power grid companies.