**1. 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–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.

From a mathematical point of view, the target problem was a typical multi-objective constrained optimization problem with a set of complex equality or inequality constraints [7], and many classical methods have been successfully developed by scholars, during the past decades [8–13], like linear programming [14], quadratic programming [15], dynamic programming [16–18], Lagrange relaxation, and network optimization [19–21]. However, the hydropower operation problem is usually modeled with nonlinear characteristic curves, physical constraints, or objective functions [22–24]. The above-mentioned traditional methods might fail to address the complexity due to various defects, like dimensionality problem [25], high computational burden [26], duality gap [27], or parameter tuning [28–30]. In recent years, with the booming development of computer technology, many evolutionary algorithms have been proposed to resolve these kind of problems [31–33], like genetic algorithm (GA) [34], differential evolution (DE) [35,36], particle swarm algorithm (PSO) [37–40], cuckoo search (CS) [41], Covariance Matrix Adaptation Evolution Strategy with a Directed Target to Best Perturbation (CMA-ES-DTBP) [42], and a clustered adaptive teaching–learning-based optimization (CATLBO) [43]. Compared with traditional methods, the evolutionary algorithms can produce satisfying solutions in most cases, regardless of the problem features (like continuity or nonconvexity) [44–46]. However, due to the premature convergence, evolutionary algorithms often fall into local optima, which have limited their widespread applications in practical engineering [47–50]. Hence, it is necessary to find some effective modified strategies to enhance the performances of the evolutionary algorithms in the hydropower operation problems.

Gravitational search algorithm (GSA) is a famous evolutionary algorithm based on the laws of gravity and mass interactions in Newtonian mechanics [51–53]. Although GSA outperforms the standard PSO and GA methods in many problems, certain defects (like exploration–exploitation imbalance and local convergence) still exist. For improving the performance of the standard GSA method, this study proposes an enhanced GSA algorithm (EGSA) where three modified strategies (opposition learning strategy, mutation search strategy, and elastic-ball modification strategy) are integrated to alleviate the premature convergence problem. Several famous test functions are used to verify the feasibility of the developed method in solving the global optimization problems, and the results demonstrate that the modified strategies could effectively enhance the convergence speed and the global search ability of the population at the same time. Then, two specific strategies were embedded into the EGSA method to help address the multi-objective operation problem of a hydropower system—(i) the practical constraint handling strategy for solution feasibility and (ii) the Technique for Order Preference by Similarity to an Ideal Solution (TOPSIS) for multi-objective decision. Finally, the EGSA method was applied to address the operation optimization of cascade hydropower system in the Wu River of China. The simulations showed that the proposed method could produce better scheduling schemes in different application scenarios.

The rest of this paper is organized as below—the EGSA method is presented in Section 2; several benchmark functions are used to verify the performance of the EGSA method in Section 3; Section 4 gives the mathematical model for the multi-objective operation of cascade hydropower reservoirs; Section 5 compares the results of different methods; and the conclusions are drawn in the end.

#### **2. Enhanced Gravitational Search Algorithm (EGSA)**

#### *2.1. Gravitational Search Algorithm (GSA)*

GSA is a novel swarm-based evolutionary algorithm inspired by the universal gravitation and mass interactions in classical Newton's Mechanics [54]. In GSA, the potential solution for the target optimization problems is considered to be an agent with a certain quality, and all the agents are attracted by other agents via the so-called gravity force, which can effectively help the swarm share the acquired beneficial information about the surrounding [55]. Figure 1 draws the sketch map of the GSA method. It can be seen that during the evolutionary process, the gravity force of the two agents is directly and inversely proportional to their qualities and distances, respectively; while the agents with better performances usually have larger attraction forces and slower velocity than those agents with smaller quality. By these means, all agents gradually finish the transition from the global exploration to local exploitation in the decision space, which help the swarm obtain the global optima.

**Figure 1.** Sketch map of the gravitational search algorithm (GSA) method for a two-variable problem.

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 *i*th agent at the *k*th iteration (*Xi*(*k*) for short) can be expressed as follows:

$$X\_i(k) = (x\_i^1(k), \dots, x\_i^d(k), \dots, x\_i^D(k))\tag{1}$$

where *x<sup>d</sup> <sup>i</sup>* (*k*) is the position of the *i*th agent in the *d*th dimension at the *k*th iteration.

In the evolutionary process, the gravitational force of the *n*th agent acting on the *i*th agent in the *d*th dimension at the *k*th iteration can be expressed as follows:

$$F\_{\rm in}^d(k) = G(k)\frac{M\_i(k) \times M\_n(k)}{R\_{\rm in}(k) + q} (\mathbf{x}\_i^d(k) - \mathbf{x}\_n^d(k))\tag{2}$$

where *Mi*(*k*) and *Mn*(*k*) are the gravitational masses of the *i*th agent and the *n*th agent at the *k*th iteration, respectively. *Rin*(*k*) =||*Xi*(*k*), *Xn*(*k*)|| <sup>2</sup> is the Euclidian distance between the *i*th agent and the *n*th agent at the *k*th iteration. ϕ is a small constant used to prevent the denominator being 0. *G*(*k*) is the gravitational constant of the swarm at the *k*th iteration, which is dynamically updated by:

$$\mathbf{G}(k) = \mathbf{G}\_0 \times \exp(-\alpha \frac{k}{\overline{k}}) \tag{3}$$

where *G*<sup>0</sup> 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 *i*th agent at the *k*th iteration (*mi*(*k*) for short) can be expressed as follows:

$$\begin{cases} m\_i(k) = \frac{fit\_i(k) - \text{nworst}(k)}{best(k) - \text{nørst}(k)}\\\ M\_i(k) = m\_i(k) \Big| \sum\_{i=1}^{N} m\_i(k) \end{cases} \tag{4}$$

where *fiti*(*k*) is the fitness of the *i*th agent at the *k*th iteration. *best*(*k*) and *worst*(*k*) are the best and worst fitness values of the swarm at the *k*th 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 *i*th agent in the *d*th dimension at the *k*th iteration (*F<sup>d</sup> <sup>i</sup>* (*k*) for short) is obtained by:

$$F\_i^d(k) = \sum\_{j=1, j \neq i}^{Kbest} rnnd\_j F\_{ij}^d(k) \tag{5}$$

where *randj* 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 *i*th agent in the *d*th dimension at the *k*th iteration (*accd <sup>i</sup>* (*k*) for short) can be expressed as follows:

$$acc\_i^d(k) = \frac{F\_i^d(k)}{M\_i(k)}\tag{6}$$

The velocity and position values of the *i*th agent in the *d*th dimension at the *k*+1th iteration (*x<sup>d</sup> <sup>i</sup>* (*<sup>k</sup>* <sup>+</sup> <sup>1</sup>) and *vd <sup>i</sup>* (*k* + 1) for short) can be updated as below:

$$\begin{cases} \upsilon\_i^d(k+1) = rand\_i \rtimes \upsilon\_i^d(k) + acc\_i^d(k) \\ \ x\_i^d(k+1) = \mathfrak{x}\_i^d(k) + \upsilon\_i^d(k+1) \end{cases} \tag{7}$$

where *randi* is the random number uniformly distributed in the range of [0,1].

#### *2.2. 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 *<sup>x</sup>* <sup>∈</sup> [*a*, *<sup>b</sup>*], the corresponding opposite point (*x*<sup>1</sup> for short) of *<sup>x</sup>* could be easily obtained by *<sup>x</sup>*<sup>1</sup> = *<sup>a</sup>* + *<sup>b</sup>* <sup>−</sup> *<sup>x</sup>*. 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–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:

$$\begin{cases} \begin{array}{c} re\_i^d(k) = \mathsf{U}b^d + \mathsf{L}b^d - a\_i^d(k) \\\ a\_i^d(k) = c\_1 \times \mathbf{x}\_i^d(k) - c\_2 \times rand \times \left[ \mathsf{gBest}^d(k) - \mathsf{x}\_i^d(k) \right] \end{array} \tag{8}$$

where *re<sup>d</sup> <sup>i</sup>* (*k*) is the *<sup>i</sup>*th opposite agent in the *<sup>d</sup>*th dimension at the *<sup>k</sup>*th iteration. *Ub<sup>d</sup>* and *Lbd* are the upper and lower limits of the *d*th dimension. *c*<sup>1</sup> and *c*<sup>2</sup> are two learning factors. *rand* is the random number uniformly distributed in the range of [0,1]. *gBestd*(*k*) is the *d*th dimension of the global optimal agent at the *k*th iteration.

#### *2.3. 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.

$$B\_i^d(k) = pBest\_{\lambda}^d(k) + r\_1 \times (pBest\_i^d(k) - \mathcal{gBest}^d(k))\tag{9}$$

where *Bd <sup>i</sup>* (*k*) is the position of the *<sup>i</sup>*th mutated agent in the *<sup>d</sup>*th dimension at the *<sup>k</sup>*th iteration. *pBestd <sup>i</sup>* (*k*) is the best-known position of the *i*th agent in the *d*th dimension at the *k*th iteration. λ is the number randomly selected from the index set {1, 2, ··· , *N*}. *r*<sup>1</sup> is the random number uniformly distributed in the range of [−0.5, 0.5].

#### *2.4. 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.

$$\mathbf{x}\_{i}^{d}(k) = \begin{cases} \, \text{l}lb^{d} - r\_{2} \times \Delta\_{i}^{d}(k) \, \text{if} \begin{pmatrix} \mathbf{x}\_{i}^{d}(k) > \text{l}lb^{d} \\ \mathbf{I}b^{d} + r\_{2} \times \Delta\_{i}^{d}(k) \, \text{if} \begin{pmatrix} \mathbf{x}\_{i}^{d}(k) < \text{l}lb^{d} \end{pmatrix} \end{cases} \tag{10}$$

$$\Delta\_i^d(k) = \begin{cases} \mathbf{x}\_i^d(k) - \mathsf{U}lb^d \text{ if} \{ \mathbf{x}\_i^d(k) > \mathsf{U}lb^d \} \\\ \mathsf{L}b^d - \mathsf{x}\_i^d(k) \text{ if} \{ \mathbf{x}\_i^d(k) < \mathsf{L}lb^d \} \end{cases} \tag{11}$$

where *r* <sup>2</sup> is the random number uniformly distributed in the range of [0,1].

#### *2.5. 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.

#### **3. Numerical Experiments to Verify the Performance of the EGSA Method**

#### *3.1. 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.

#### *3.2. 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*<sup>1</sup> 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.

#### *3.3. Comparison with Other Evolutionary Algorithms in Small-Scale Problems*

#### 3.3.1. 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*<sup>9</sup> 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.

#### 3.3.2. 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.

**Figure 2.** *Cont*.

**Figure 2.** Box and whisker test of three algorithms for 12 functions with 30 variables.

#### 3.3.3. 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<sup>+</sup> than R– in 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.



*Water* **2019**, *11*, 2040



R+ R− *p*-value Significant

√√

72

6

6.84 × 10−3

4.88 × 10−4

4.88 × 10−4

 √

 √

 √

 √

 √

 √

 √

 √

 √

1.22 × 10−2

9.77 × 10−4

9.77 × 10−4

2.69 × 10−2

4.88 × 10−4

4.88 × 10−4

4.88 × 10−4

4.88 × 10−4

 78

 0

 78

 0

 70

 8

 77

 1

 77.5

 0.5

 67

 11

 78

 0

 78

 0

 78

 0

 78

 0

**Item**

 **EGSA–CS**

 **EGSA–MCS**

 **EGSA–LSA**

 **EGSA–GWO**

 **EGSA–FA**

 **EGSA–WOA**

 **EGSA–ALO**

 **EGSA–DE**

 **EGSA–PSO**

 **EGSA–SCA**

 **EGSA–GSA**

#### 3.3.4. Convergence Analysis

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*<sup>4</sup> 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*<sup>9</sup> 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.** Convergence trajectories of the different methods for 12 benchmark functions with 30 variables.
