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

Abstract

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 modification 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 different cases. Hence, an effective optimization method was provided for the multi-objective operation of hydropower system.

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,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.

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,9,10,11,12,13], like linear programming [14], quadratic programming [15], dynamic programming [16,17,18], Lagrange relaxation, and network optimization [19,20,21]. However, the hydropower operation problem is usually modeled with nonlinear characteristic curves, physical constraints, or objective functions [22,23,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,29,30]. In recent years, with the booming development of computer technology, many evolutionary algorithms have been proposed to resolve these kind of problems [31,32,33], like genetic algorithm (GA) [34], differential evolution (DE) [35,36], particle swarm algorithm (PSO) [37,38,39,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,45,46]. However, due to the premature convergence, evolutionary algorithms often fall into local optima, which have limited their widespread applications in practical engineering [47,48,49,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,52,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.

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:

$$X_i(k)=(x_i^1(k),\cdots ,x_i^d(k),\cdots ,x_i^D(k))$$

(1)

where $x_i^d(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:

$$F_{in}^d(k)=G(k)\frac{M_i(k)\times M_n(k)}{R_{in}(k)+\phi }(x_i^d(k)-x_n^d(k))$$

(2)

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)=\left|\right|X_i\left(k\right),X_n(k{\left)\right||}_2$ is the Euclidian distance between the ith agent and the nth agent at the kth iteration. $\phi$ 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:

$$G(k)=G_0\times \exp (-\alpha \frac{k}{\overline{k}})$$

(3)

where $G_0$ is the initial gravitational constant. $\alpha$ is the attenuation coefficient. $\overline{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:

$$\{\begin{array}{l}{m}_i(k)=\frac{fi{t}_i(k)-worst(k)}{best(k)-worst(k)}\\ M_i(k)=m_i(k)/{\displaystyle \sum _{i=1}^N{m}_i(k)}\end{array}$$

(4)

where $fi{t}_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.

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_i^d(k)$ for short) is obtained by:

$$F_i^d(k)={\displaystyle \sum _{j=1,j\ne i}^{Kbest}ran{d}_j{F}_{ij}^d(k)}$$

(5)

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 ($ac{c}_i^d(k)$ for short) can be expressed as follows:

$$ac{c}_i^d(k)=\frac{F_i^d(k)}{M_i(k)}$$

(6)

The velocity and position values of the ith agent in the dth dimension at the k+1th iteration ($x_i^d(k+1)$ and $v_i^d(k+1)$ for short) can be updated as below:

$$\{\begin{array}{l}{v}_i^d(k+1)=ran{d}_i\times v_i^d(k)+ac{c}_i^d(k)\\ x_i^d(k+1)=x_i^d(k)+v_i^d(k+1)\end{array}$$

(7)

where rand_i 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 $x\in \left[a,b\right]$, the corresponding opposite point (x¹ 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:

$$\{\begin{array}{l}r{e}_i^d(k)=U{b}^d+L{b}^d-a_i^d(k)\\ a_i^d(k)=c_1\times x_i^d(k)-c_2\times rand\times [gBes{t}_{}^d(k)-x_i^d(k)]\end{array}$$

(8)

where $r{e}_i^d(k)$ is the ith opposite agent in the dth dimension at the kth iteration. $U{b}^d$ and $L{b}^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]. $gBes{t}_{}^d(k)$ is the dth dimension of the global optimal agent at the kth 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 $\upsilon$; second, the first cbest agents in $\upsilon$ directly enter the next cycle, while the left agents in $\upsilon$ 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)=pBes{t}_{\lambda }^d(k)+r_1\times (pBes{t}_i^d(k)-gBes{t}_{}^d(k))$$

(9)

where $B_i^d(k)$ is the position of the ith mutated agent in the dth dimension at the kth iteration. $pBes{t}_i^d(k)$ is the best-known position of the ith agent in the dth dimension at the kth iteration. $\lambda$ is the number randomly selected from the index set $\left\{1,2,\cdots ,N\right\}$. r₁ 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.

$$x_i^d(k)=\{\begin{array}{l}U{b}^d-r_2^{}\times {\Delta }_i^d(k)\hspace{1em}\mathrm{if}\left(x_i^d(k)>U{b}^d\right)\hfill \\ L{b}^d+r_2^{}\times {\Delta }_i^d(k)\hspace{1em}\mathrm{if}\left(x_i^d(k)<L{b}^d\right)\hspace{1em}\hfill \end{array}$$

(10)

$${\Delta }_i^d(k)=\{\begin{array}{l}{x}_i^d(k)-U{b}^d\hspace{1em}\mathrm{if}\left(x_i^d(k)>U{b}^d\right)\hfill \\ L{b}^d-x_i^d(k)\hspace{1em}\mathrm{if}\left(x_i^d(k)<L{b}^d\right)\hspace{1em}\end{array}$$

(11)

where $r_2^{}$ 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. Detailed information of the twelve benchmark functions.

Function	D	Range	fmin
$F_1(x)={\displaystyle \sum _{i=1}^n{x}_i{}^2}$	30	[−100,100]	0
$F_2\left(x\right)={\displaystyle \sum _{i=1}^n\left|x_i\right|}+{\displaystyle \prod _{i=1}^n\left|x_i\right|}$	30	[−10,10]	0
$F_3(x)={\displaystyle \sum _{i=1}^n({\displaystyle \sum _{j=1}^i{x}_j}}{)}^2$	30	[−100,100]	0
$F_4\left(x\right)=\max \left\{\left|x_i\right|,1\le i\le n\right\}$	30	[−100,100]	0
$F_5(x)={\displaystyle \sum _{i=1}^{n-1}[100{(x_{i+1}^{}-x_i^2)}^2+{(x_i^{}-1)}^2]}$	30	[−30,30]	0
$F_6(x)={\displaystyle \sum _{i=1}^n{(x_i+0.5)}^2}$	30	[−100,100]	0
$F_7(x)={\displaystyle \sum _{i=1}^{n}i{x}_i^4}+random[0,1)$	30	[−1.28,1.28]	0
$F_8(x)={\displaystyle \sum _{i=1}^n-x_i}\sin (\sqrt{\left|x_i\right|})$	30	[−500,500]	−418.9 × D
$F_9(x)={\displaystyle \sum _{i=1}^n[x_i{}^2-10\cos (2\pi x_i)+10]}$	30	[−5.12,5.12]	0
$\begin{array}{l}{F}_{10}(x)=-20\exp (-0.2\sqrt[]{\frac{1}{n}{\displaystyle \sum _{i=1}^n{x}_i^2}})-\exp (\frac{1}{n}{\displaystyle \sum _{i=1}^n\cos (2\pi x_i))}\\ \hspace{1em}\hspace{1em}\hspace{1em}+20+e\end{array}$	30	[−32,32]	0
$F_{11}(x)=\frac{1}{4000}{\displaystyle \sum _{i=1}^n{x}_i^2-{\displaystyle \prod _{i=1}^n\cos (\frac{x_i}{\sqrt{i}}}})+1$	30	[−600,600]	0
$\begin{array}{l}{F}_{12}\left(x\right)=\frac{\pi }{n}\{10{\sin }^2(\pi y_1)+{\displaystyle \sum _{i=1}^{n-1}}{(y_i-1)}^2[1+10{\sin }^2(\pi y_{i+1})]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{(y_n-1)}^2\}+{\displaystyle \sum _{i=1}^{n}u(x_i},10,100,4)\\ \hspace{1em}\hspace{1em}\hspace{1em}{y}_i=1+\frac{x_i+1}{4},u(x_i,a,k,m)=\{\begin{array}{l}k{(x_i-a)}^m{x}_i>a\\ 0a\le x_i<a\\ \hspace{1em}k{(-x_i-a)}^m{x}_i<-a\end{array}\end{array}$	30	[−50,50]	0

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₁–F₈) with one global optimum, and multimodal functions (F₉–F₁₂) 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.

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₁ and c₂) 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. Statistical indexes of 12 algorithms for all the test functions.

Function	Item	CS	MCS	LSA	GWO	FA	WOA	ALO	DE	PSO	SCA	GSA	EGSA
F1	Ave.	9.06 × 101	1.01 × 100	4.81 × 10−8	6.59 × 10−28	1.20 × 10−2	1.41 × 10−30	2.59 × 10−10	7.80 × 10−6	4.58 × 10−7	2.87 × 10−35	4.00 × 10−9	6.96 × 10−134
	Std.	2.62 × 101	2.72 × 10−1	3.40 × 10−7	6.34 × 10−5	4.30 × 10−3	4.91 × 10−30	1.65 × 10−10	3.14 × 10−6	3.51 × 10−7	1.26 × 10−34	2.92 × 10−9	3.70 × 10−133
F2	Ave.	9.70 × 100	1.81 × 10−1	3.68 × 10−2	7.18 × 10−17	3.73 × 10−1	1.06 × 10−21	1.84 × 10−6	3.73 × 10−4	1.19 × 10−4	2.21 × 10−27	5.48 × 10−5	5.21 × 10−69
	Std.	1.98 × 100	3.31 × 10−2	1.56 × 10−1	2.90 × 10−2	1.01 × 10−1	2.39 × 10−21	6.58 × 10−7	9.05 × 10−5	8.03 × 10−5	6.99 × 10−27	2.24 × 10−5	2.09 × 10−68
F3	Ave.	3.84 × 103	4.62 × 102	4.32 × 101	3.29 × 10−6	1.81 × 103	5.39 × 10−7	6.06 × 10−10	2.94 × 104	2.33 × 100	9.74 × 101	3.61 × 102	4.30 × 10−119
	Std.	7.24 × 102	1.23 × 102	2.99 × 101	7.91 × 101	6.60 × 102	2.93 × 10−6	6.34 × 10−10	4.13 × 103	1.14 × 100	2.24 × 102	1.45 × 102	2.29 × 10−118
F4	Ave.	7.23 × 100	1.73 × 100	1.49 × 100	5.61 × 10−7	7.67 × 10−2	7.26 × 10−2	1.36 × 10−8	1.43 × 100	5.35 × 10−1	2.66 × 100	9.48 × 10−2	5.58 × 10−70
	Std.	6.76 × 10−1	5.12 × 10−1	1.30 × 100	1.32 × 100	1.46 × 10−2	3.97 × 10−1	1.81 × 10−9	3.21 × 10−1	1.33 × 10−1	4.16 × 100	3.92 × 10−1	2.12 × 10−69
F5	Ave.	4.98 × 100	5.81 × 101	6.43 × 101	2.68 × 101	1.28 × 102	2.79 × 101	3.47 × 10−1	4.45 × 102	1.64 × 102	2.74 × 101	3.45 × 101	2.69 × 101
	Std.	1.75 × 100	3.31 × 101	4.38 × 101	6.99 × 101	2.79 × 102	7.64 × 10−1	1.10 × 10−1	1.55 × 102	1.55 × 102	3.85 × 10−1	3.65 × 101	9.72 × 10−2
F6	Ave.	4.31 × 104	4.44 × 103	3.34 × 100	8.17 × 10−1	0	3.12 × 100	2.56 × 10−10	7.94 × 10−6	9.29 × 10−7	5.30 × 10−2	4.99 × 10−9	8.23 × 10−15
	Std.	7.21 × 103	7.52 × 102	2.09 × 100	1.26 × 10−4	0	5.32 × 10−1	1.09 × 10−10	3.39 × 10−6	1.93 × 10−6	4.98 × 10−2	3.74 × 10−9	4.61 × 10−15
F7	Ave.	2.46 × 10−2	9.10 × 10−3	2.41 × 10−2	2.21 × 10−3	3.52 × 10−2	1.43 × 10−3	4.29 × 10−3	1.35 × 10−1	3.55 × 100	7.48 × 10−3	3.05 × 10−2	4.76 × 10−4
	Std.	7.90 × 10−3	2.20 × 10−3	5.72 × 10−3	1.00 × 10−1	2.40 × 10−2	1.15 × 10−3	5.09 × 10−3	2.93 × 10−2	4.79 × 100	6.96 × 10−3	1.59 × 10−2	2.32 × 10−4
F8	Ave.	−8.98 × 103	−9.80 × 103	−8.00 × 103	−6.12 × 103	−5.90 × 103	−5.08 × 103	−1.61 × 103	−7.95 × 103	−6.39 × 103	−5.66 × 103	−2.75 × 103	−1.19 × 104
	Std.	1.98 × 102	5.31 × 102	6.69 × 102	4.09 × 103	6.56 × 102	6.96 × 102	3.14 × 102	3.33 × 102	1.29 × 103	3.52 × 102	4.00 × 102	3.13 × 102
F9	Ave.	2.94 × 102	1.35 × 102	6.28 × 101	3.11 × 10−1	2.63 × 101	0	7.71 × 10−6	1.32 × 102	7.87 × 101	1.87 × 10−9	1.67 × 101	0
	Std.	1.43 × 101	2.16 × 101	1.49 × 101	4.74 × 101	9.15 × 100	0	8.45 × 10−6	8.78 × 100	2.91 × 101	1.02 × 10−8	3.82 × 100	0
F10	Ave.	1.93 × 101	1.21 × 101	2.69 × 100	1.06 × 10−13	5.12 × 10−2	7.40 × 100	3.73 × 10−15	1.37 × 10−3	5.99 × 10−4	1.34 × 101	3.25 × 10−5	3.64 × 10−15
	Std.	3.50 × 10−1	7.52 × 10−1	9.11 × 10−1	7.78 × 10−2	1.37 × 10−2	9.90 × 100	1.50 × 10−15	8.85 × 10−4	5.43 × 10−4	9.64 × 100	1.10 × 10−5	1.08 × 10−15
F11	Ave.	2.12 × 102	8.32 × 100	7.24 × 10−3	4.49 × 10−3	5.84 × 10−3	2.89 × 10−4	1.86 × 10−2	3.36 × 10−3	6.98 × 10−3	1.86 × 10−5	4.34 × 100	0
	Std.	3.97 × 101	1.54 × 100	6.70 × 10−3	6.66 × 10−3	1.43 × 10−3	1.59 × 10−3	9.55 × 10−3	1.06 × 10−2	8.16 × 10−3	1.02 × 10−4	1.66 × 100	0
F12	Ave.	1.47 × 100	1.38 × 10−1	3.58 × 10−1	5.34 × 10−2	2.40 × 10−4	3.40 × 10−1	9.74 × 10−12	6.87 × 10−1	2.01 × 10−2	8.52 × 10−3	9.59 × 10−2	5.30 × 10−17
	Std.	3.61 × 10−1	2.86 × 10−1	7.44 × 10−1	2.07 × 10−2	1.00 × 10−4	2.15 × 10−1	9.33 × 10−12	3.05 × 10−1	2.26 × 10−2	4.69 × 10−3	1.75 × 10−1	1.92 × 10−17

Note: Ave. = average; Std. = standard deviation.

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₅, the WOA performance was tied with EGSA in F₉ 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₄–F₅, 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₃, 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.

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. The results of Wilcoxon test for a single problem for statistical significance level at alpha = 0.05.

Function	EGSA–CS	EGSA–MCS	EGSA–LSA	EGSA–GWO	EGSA–FA	EGSA–WOA	EGSA–ALO	EGSA–DE	EGSA–PSO	EGSA–SCA	EGSA–GSA
F1	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F2	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F3	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F4	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F5	CS	EGSA	EGSA	GWO	EGSA	EGSA	ALO	EGSA	EGSA	EGSA	EGSA
F6	EGSA	EGSA	EGSA	EGSA	FA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F7	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F8	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F9	EGSA	EGSA	EGSA	EGSA	EGSA	Tie	EGSA	EGSA	EGSA	EGSA	EGSA
F10	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F11	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
F12	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA	EGSA
W/T/L	11/0/1	12/0/0	12/0/0	11/0/1	11/0/1	11/1/0	11/0/1	12/0/0	12/0/0	12/0/0	12/0/0

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. The results of Wilcoxon signed-rank test for statistically significance level at alpha = 0.05.

Item	EGSA–CS	EGSA–MCS	EGSA–LSA	EGSA–GWO	EGSA–FA	EGSA–WOA	EGSA–ALO	EGSA–DE	EGSA–PSO	EGSA–SCA	EGSA–GSA
R+	72	78	78	70	77	77.5	67	78	78	78	78
R−	6	0	0	8	1	0.5	11	0	0	0	0
p-value	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
Significant	√	√	√	√	√	√	√	√	√	√	√

, 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 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.

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₁–F₄ 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₉ and F₁₁, and better results than other methods in F₈, F₁₀, and F₁₂. Thus, the above analysis fully proves that the presented method has a superior convergence speed and global search ability.

4. EGSA for the Multi-Objective Operation of Cascade Hydropower Reservoirs

4.1. Mathematical Model

4.1.1. 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:

$$E=\max {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{n=1}^N{P}_{n,t}^{}\times {\Delta }_t}}$$

(12)

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. ${\Delta }_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:

$$F=\min \sqrt{\frac{1}{2}{\displaystyle \sum _{t=1}^T{\left[Loa{d}_t-{\displaystyle \sum _{n=1}^N{P}_{n,t}}\right]}^2}}$$

(13)

where $Loa{d}_t$ is the load demand of the power system at the tth period.

4.1.2. Physical Constraints

(1) Water Balance Equation:

$$V_{n,t+1}=V_{n,t}+3600{\Delta }_t\times \left[I_{n,t}+{\displaystyle \sum _{i=1}^{Num}(Q_{i,t}+S_{i,t})}-Q_{n,t}-S_{n,t}\right]$$

(14)

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:

$${\underset{\_}{V}}_{n,t}\le V_{n,t}\le {\overline{V}}_{n,t}$$

(15)

where ${\overline{V}}_{n,t}$ and ${\underset{\_}{V}}_{n,t}$ denote the maximum and minimum storage volumes of the nth hydroplant at the tth period, respectively.

(3) Turbine Discharge Constraint:

$${\underset{\_}{Q}}_{n,t}\le Q_{n,t}\le {\overline{Q}}_{n,t}$$

(16)

where ${\overline{Q}}_{n,t}$ and ${\underset{\_}{Q}}_{n,t}$ denote the maximum and minimum turbine discharges of the nth hydroplant at the tth period, respectively.

(4) Total Discharge Constraint:

$${\underset{\_}{q}}_{n,t}\le q_{n,t}\le {\overline{q}}_{n,t}$$

(17)

where ${\overline{q}}_{n,t}$ and ${\underset{\_}{q}}_{n,t}$ denote the maximum and minimum total discharges of the nth hydroplant at the tth period, respectively.

(5) Power Output Constraint:

$${\underset{\_}{P}}_{n,t}\le P_{n,t}\le {\overline{P}}_{n,t}$$

(18)

where ${\overline{P}}_{n,t}$ and ${\underset{\_}{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:

$$\{\begin{array}{l}{Z}_{n,0}=Z_n^{beg}\\ Z_{n,T}=Z_n^{end}\end{array}$$

(19)

where $Z_n^{beg}$ and $Z_n^{end}$ represent the initial and final water level of the nth hydroplant, respectively.

4.2. 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:

$$A={\left(a_{i,j}\right)}_{m\times J}=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1J}\\ a_{21}& a_{22}& \cdots & a_{2J}\\ ⋮& ⋮& a_{i,j}& ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mJ}\end{array}\right]$$

(20)

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:

$$B={\left(b_{i,j}\right)}_{m\times J}=\left[\begin{array}{cccc}{b}_{11}& b_{12}& \cdots & b_{1J}\\ b_{21}& b_{22}& \cdots & b_{2J}\\ ⋮& ⋮& b_{i,j}& ⋮\\ b_{m1}& b_{m2}& \cdots & b_{mJ}\end{array}\right]$$

(21)

$$b_{ij}=\{\begin{array}{l}\frac{a_j^{\max }-a_{ij}}{a_j^{\max }-a_j^{\min }}\hspace{1em}\mathrm{minimization} \mathrm{attribute}\\ \frac{a_{ij}-a_j^{\min }}{a_j^{\max }-a_j^{\min }}\hspace{1em}\mathrm{maximization} \mathrm{attribute}\end{array},i=1,2,\cdots ,m;j=1,2,\cdots ,J$$

(22)

where $a_j^{\max }$ and $a_j^{\min }$ are the maximum and minimum values of the jth attribute.

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:

$$R={\left(r_{i,j}\right)}_{m\times J}=\left[\begin{array}{cccc}{r}_{11}& r_{12}& \cdots & r_{1J}\\ r_{21}& r_{22}& \cdots & r_{2J}\\ ⋮& ⋮& r_{i,j}& ⋮\\ r_{m1}& r_{m2}& \cdots & r_{mJ}\end{array}\right]$$

(23)

$$r_{ij}=\frac{b_{ij}}{\sqrt{{\displaystyle \sum _{i=1}^m{b}_{ij}^2}}},i=1,2,\cdots ,m; j=1,2,\cdots ,J$$

(24)

where $r_{ij}$ is the normalized value of the jth attribute in the ith scheme.

Step 4: Set the weight vector $w=\left[w_1,\cdots ,w_j,\cdots ,w_J\right]$, where $w_j$ is the weight of the jth attribute and there is ${\displaystyle \sum _{j=1}^J{w}_j}=1$. Then, calculate the weighted normalized decision-making matrix V by:

$$V={\left(v_{i,j}\right)}_{m\times J}={\left(w_j\times r_{i,j}\right)}_{m\times J}$$

(25)

Step 5: Determine the positive ideal scheme $C^{+}$ and negative ideal scheme $C^{-}$ by:

$$\{\begin{array}{c}{C}^{+}=\left[c_1^{+},c_2^{+},\dots ,c_J^{+}\right]=\max \left\{v_{i,j}|i=1,2,\cdots ,m\right\}\\ C^{-}=\left[c_1^{-},c_2^{-},\dots ,c_J^{-}\right]=\min \left\{v_{i,j}|i=1,2,\cdots ,m\right\}\end{array}j=1,2,\dots ,J$$

(26)

Step 6: Calculate the Euclidean distances between the decision scheme and the positive (negative) ideal scheme, which could be expressed as below:

$$\{\begin{array}{c}{D}_i^{+}=\sqrt{{\displaystyle \sum _{j=1}^J{(v_{ij}-c_j^{+})}^2}}\\ D_i^{-}=\sqrt{{\displaystyle \sum _{j=1}^J{(v_{ij}-c_j^{-})}^2}}\end{array},i=1,2,\dots m$$

(27)

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

$$f_i=\frac{D_i^{-}}{D_i^{-}+D_i^{+}},i=1,2,\dots m$$

(28)

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.

$$f_{best}=\max \{f_1,f_2,\dots ,f_m\}$$

(29)

4.3. Details of EGSA for Multi-Objective Operation of Cascade Hydropower Reservoirs

4.3.1. 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:

$$X_i(k)=\left[\begin{array}{cccc}{q}_{1,1}& q_{1,2}& \cdots & q_{1,T}\\ q_{2,1}& q_{2,2}& \cdots & q_{2,T}\\ ⋮& ⋮& q_{n,t}& ⋮\\ q_{N,1}& q_{N,2}& \cdots & q_{N,T}\end{array}\right]$$

(30)

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:

$$q_{n,t}={\underset{\_}{q}}_{n,t}+random\times \left({\overline{q}}_{n,t}-{\underset{\_}{q}}_{n,t}\right)$$

(31)

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

4.3.2. 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

$$V_{n,t}=\max \left\{{\underset{\_}{V}}_{n,t},\min \left\{V_{n,t-1}+3600{\Delta }_{t-1}\times \left[I_{n,t-1}+{\displaystyle \sum _{i=1}^{Num}(Q_{i,t-1}+S_{i,t-1})}-Q_{n,t-1}-S_{n,t-1}\right],{\overline{V}}_{n,t}\right\}\right\}$$

(32)

Step 3: Set L = L + 1, and then calculate the violation value ${\eta }_n$ of the final storage by Equation (33). Then, if $\left|{\eta }_n\right|$ is smaller than the accuracy $\epsilon$ or L is larger than the maximum iteration, go to Step 5; otherwise, turn to Step 4.

$${\eta }_n=f_n\left(Z_{n,T}\right)-f_n\left(Z_n^{end}\right)$$

(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.

$$q_{n,t}=\max \left\{{\underset{\_}{q}}_{n,t},\min \left\{q_{n,t}+\frac{{\eta }_n}{T\times 3600{\Delta }_t},{\overline{q}}_{n,t}\right\}\right\}$$

(34)

Step 5: Set n = n+1. If n ≤ N, turn to Step 2 for the new hydroplant; otherwise, stop the iteration.

4.3.3. 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\left[X_i(k)\right]$ 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^’).

$$Viol\left[X_i(k)\right]={\displaystyle \sum _{l_1=1}^{L_g}{c}_{l_1}·\max \left\{g_{l_1}\left(X_i(k)\right),0\right\}}+{\displaystyle \sum _{l_2=1}^{L_e}{c}_{l_2}·{\left[e_{l_2}\left(X_i(k)\right)\right]}^2}$$

(35)

$$\{\begin{array}{l}{E}^{\prime }\left[X_i(k)\right]={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{n=1}^N{P}_{n,t}^{}\times \Delta t}}-Viol\left[X_i(k)\right]\\ F^{\prime }\left[X_i(k)\right]=\sqrt{\frac{1}{2}{\displaystyle \sum _{t=1}^T{\left[Loa{d}_t-{\displaystyle \sum _{n=1}^N{P}_{n,t}}\right]}^2}}+Viol\left[X_i(k)\right]\end{array}$$

(36)

where $g_{l_1}(·)\le 0$ is the l₁th inequality constraint; $e_{l_2}(·)\le 0$ is the l₂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.

4.3.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.

5. Case Studies

5.1. 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 Figure 5 and Figure 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.

5.2. Case Study 1: Power Generation of Cascade Hydropower Reservoirs

5.2.1. 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. Statistical results of the 5 methods in 20 independent runs for 4 cases of power generation (10⁴ kW·h).

Case	Method	Best	Worst	Average	Standard Deviation	Range
Case 1	DE	5388.10	5386.25	5387.19	0.56	1.85
	PSO	5393.67	5391.27	5392.49	0.57	2.40
	SCA	5388.26	5384.34	5385.65	0.89	3.92
	GSA	5388.37	5385.34	5386.76	0.88	3.03
	EGSA	5396.87	5396.81	5396.86	0.02	0.06
Case 2	DE	6439.12	6434.22	6436.14	1.43	4.90
	PSO	6445.48	6442.35	6444.12	0.79	3.13
	SCA	6439.06	6435.25	6437.07	1.03	3.81
	GSA	6440.17	6437.00	6438.55	0.79	3.17
	EGSA	6450.34	6450.33	6450.33	0.01	0.01
Case3	DE	4977.94	4976.38	4977.17	0.40	1.56
	PSO	4981.62	4979.71	4980.58	0.55	1.91
	SCA	4975.51	4973.58	4974.54	0.58	1.93
	GSA	4976.74	4973.83	4975.25	0.71	2.91
	EGSA	4984.56	4984.54	4984.55	0.01	0.02
Case4	DE	4433.08	4431.68	4432.32	0.44	1.40
	PSO	4435.89	4434.23	4435.11	0.49	1.66
	SCA	4431.51	4428.62	4429.58	0.70	2.89
	GSA	4431.61	4429.07	4430.37	0.64	2.54
	EGSA	4438.33	4438.32	4438.33	0.01	0.01

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.

5.2.2. Comparison of the Optimal Results Obtained by Different Evolutionary Algorithms

Table 6. Detailed results in the best solutions by different methods (10⁴ kW·h).

Runoff	Method	HJD	DF	SFY	WJD	GPT	Sum
Case 1	DE	726.73	646.19	527.41	1106.65	2381.12	5388.10
	PSO	728.66	647.76	526.51	1107.99	2382.75	5393.67
	SCA	727.92	644.60	527.68	1106.77	2381.29	5388.26
	GSA	728.35	647.35	526.77	1107.15	2378.75	5388.37
	EGSA	728.84	646.94	530.33	1108.10	2382.66	5396.87
Case 2	DE	869.46	770.94	632.25	1322.47	2844.00	6439.12
	PSO	871.70	774.09	627.59	1325.38	2846.72	6445.48
	SCA	871.15	771.70	628.34	1322.76	2845.11	6439.06
	GSA	871.48	774.01	627.42	1324.03	2843.23	6440.17
	EGSA	871.58	773.17	633.84	1324.91	2846.84	6450.34
Case 3	DE	234.57	400.46	461.27	1052.33	2829.31	4977.94
	PSO	235.21	402.56	459.19	1054.13	2830.53	4981.62
	SCA	233.99	401.01	459.94	1052.01	2828.56	4975.51
	GSA	235.02	402.45	457.71	1053.20	2828.36	4976.74
	EGSA	235.50	401.67	462.62	1054.23	2830.54	4984.56
Case 4	DE	208.46	355.81	411.02	936.99	2520.8	4433.08
	PSO	208.58	358.21	409.24	938.46	2521.4	4435.89
	SCA	208.26	356.74	409.21	937.08	2520.22	4431.51
	GSA	208.82	357.66	408.43	937.47	2519.23	4431.61
	EGSA	209.13	357.35	411.88	938.51	2521.46	4438.33

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⁴ kW·h, 3.2 × 10⁴ kW·h, 8.61 × 10⁴ kW·h, and 8.5 × 10⁴ 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.

5.2.3. Convergence Analysis of Different Evolutionary Algorithms

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.

5.3. Case Study 2: Peak Operation of Cascade Hydropower Reservoirs

5.3.1. 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. Statistical results of the 5 methods in 20 independent runs for 4 cases of peak operation.

Season	Item	Best	Worst	Average	Standard Deviation	Range
Spring	DE	32,130.69	32,183.80	32,153.90	16.22	53.11
	PSO	32,052.91	32,127.22	32,084.17	20.24	74.31
	SCA	32,160.67	32,230.89	32,194.00	22.32	70.22
	GSA	32,173.70	32,224.22	32,199.76	14.41	50.52
	EGSA	31,986.47	31,986.62	31,986.53	0.04	0.15
Summer	DE	33,006.63	33,093.05	33,051.56	23.51	86.42
	PSO	32,962.41	33,057.03	32,993.63	22.72	94.62
	SCA	33,032.40	33,129.77	33,100.08	24.05	97.37
	GSA	33,069.24	33,134.56	33,097.29	19.83	65.32
	EGSA	32,888.05	32,888.16	32,888.10	0.03	0.11
Autumn	DE	36,681.47	36,752.01	36,719.49	17.51	70.54
	PSO	36,595.43	36,714.71	36,667.14	30.32	119.28
	SCA	36,745.91	36,840.43	36,781.14	24.73	94.52
	GSA	36,738.72	36,835.34	36,802.57	27.25	96.62
	EGSA	36,536.95	36,537.06	36,537.00	0.03	0.11
Winter	DE	35,825.89	35,890.14	35,857.52	14.79	64.25
	PSO	35,762.63	35,850.26	35,811.97	23.97	87.63
	SCA	35,822.72	35,950.51	35,896.25	32.32	127.79
	GSA	35,879.31	35,965.39	35,931.63	22.63	86.08
	EGSA	35,671.70	35,671.80	35,671.74	0.03	0.10

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.

5.3.2. 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. Detailed results in the best solutions as obtained by the different methods (MW).

Season	Method	Item	Peak	Valley	Peak-valley	Average	Standard Deviation
Spring		Original	13,477.93	10,101.60	3376.33	11,910.93	1281.95
	DE	Optimization	10,952.90	7716.77	3236.13	9242.92	789.32
		Reduction	2525.03	2384.83	140.20	2668.01	492.63
	PSO	Optimization	10,050.98	8169.90	1881.08	9231.64	640.11
		Reduction	3426.95	1931.7	1495.25	2679.29	641.84
	SCA	Optimization	10,752.18	7809.30	2942.88	9239.61	926.13
		Reduction	2725.75	2292.3	433.45	2671.32	355.82
	GSA	Optimization	10,737.30	7592.29	3145.01	9229.75	1058.60
		Reduction	2740.63	2509.31	231.32	2681.18	223.35
	EGSA	Optimization	9428.01	8987.22	440.79	9232.28	165.13
		Reduction	4049.92	1114.38	2935.54	2678.65	1116.82
Summer		Original	14,119.78	10,342.98	3776.80	12,170.94	1302.21
	DE	Optimization	11,151.13	8535.00	2616.13	9505.54	670.54
		Reduction	2968.65	1807.98	1160.67	2665.40	631.67
	PSO	Optimization	10,610.47	8392.24	2218.23	9492.57	673.27
		Reduction	3509.31	1950.74	1558.57	2678.37	628.94
	SCA	Optimization	10,749.92	7591.09	3158.83	9500.12	839.87
		Reduction	3369.86	2751.89	617.97	2670.82	462.34
	GSA	Optimization	11,379.10	7774.71	3604.39	9493.04	1028.32
		Reduction	2740.68	2568.27	172.41	2677.90	273.89
	EGSA	Optimization	9759.71	9234.99	524.72	9492.46	172.66
		Reduction	4360.07	1107.99	3252.08	2678.48	1129.55
Autumn		Original	14,773.84	10,786.00	3987.84	13,220.16	1528.83
	DE	Optimization	12,615.24	8928.53	3686.71	10,553.86	880.28
		Reduction	2158.6	1857.5	301.1	2666.3	648.6
	PSO	Optimization	11,234.24	9366.35	1867.89	10,547.10	613.62
		Reduction	3539.6	1419.65	2119.95	2673.06	915.21
	SCA	Optimization	11,891.57	8142.23	3749.34	10,552.90	1099.32
		Reduction	2882.27	2643.77	238.50	2667.26	429.51
	GSA	Optimization	12,649.33	8347.14	4302.19	10,541.04	1193.12
		Reduction	2124.51	2438.86	−314.35	2679.12	335.71
	EGSA	Optimization	10,760.16	10,182.25	577.91	10,545.07	221.99
		Reduction	4013.68	603.75	3409.93	2675.09	1306.84
Winter		Original	14,913.26	11,028.24	3885.02	12,971.28	1489.76
	DE	Optimization	11,709.53	8772.72	2936.81	10,303.38	955.49
		Reduction	3203.73	2255.52	948.21	2667.90	534.27
	PSO	Optimization	11,884.52	8637.32	3247.20	10,293.13	970.91
		Reduction	3028.74	2390.92	637.82	2678.15	518.85
	SCA	Optimization	12,391.09	8667.77	3723.32	10,301.54	1052.86
		Reduction	2522.17	2360.47	161.70	2669.74	436.90
	GSA	Optimization	11,946.42	8414.17	3532.25	10,293.50	1174.01
		Reduction	2966.84	2614.07	352.77	2677.78	315.75
	EGSA	Optimization	10,568.93	10,012.57	556.36	10,295.50	208.98
		Reduction	4344.33	1015.67	3328.66	2675.78	1280.78

Note: Reduction = original–optimization of the target method.

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

5.3.3. Rationality Analysis of the Best Results Obtained by the Different Evolutionary Algorithms

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.

5.4. Case Study 3: Mutli-Objective Operation of Cascade Hydropower Reservoirs

5.4.1. Comparative Analysis of the Optimal Results Obtained by the Different Methods with 100 Weights

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 executed in 100 different weight combinations. Specifically, the weight w₁ for power generation was increased from 0 to 1.0 at the same interval of 0.01, while the weight w₂ for peak operation was set as 1.0 – w₁. From Figure 11, it was observed that the total generation was gradually decreasing with the 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.

5.4.2. Rationality Analysis of the Best Results Obtained by the Different Evolutionary Algorithms

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.

6. 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.