Long-Term Hydropower Generation Scheduling of Large-Scale Cascade Reservoirs Using Chaotic Adaptive Multi-Objective Bat Algorithm

: With growing concerns over renewable energy, the cascade hydropower reservoirs operation (CHRO), which balances the development of economic beneﬁts and power supply security, plays an increasingly important role in hydropower systems. Due to conﬂicting objectives and complicated operation constraints, the CHRO problem considering the requirements of maximizing power generation beneﬁt and ﬁrm power output is determined as a multi-objective optimization problem (MOP). In this paper, a chaotic adaptive multi-objective bat algorithm (CAMOBA) is proposed to solve the CHRO problem, and the external archive set is added to preserve non-dominant solutions. Meanwhile, population initialization based on the improved logical mapping function is adopted to improve population diversity. Furthermore, the self-adaptive local search strategy and mutation operation are designed to escape local minima. The CAMOBA is applied to the CHRO problem of the Qingjiang cascade hydropower stations in southern China. The results show that CAMOBA outperforms the multi-objective bat algorithm (MOBA) and non-dominated sorting genetic algorithms-II (NSGA-II) in di ﬀ erent hydrological years. The spacing (SP) and hypervolume (HV) metrics verify the excellent performance of CAMOBA in diversity and convergence. In summary, the CAMOBA is demonstrated to get better scheduling solutions, providing an e ﬀ ective approach for solving the cascade hydropower reservoirs operation (CHRO). including MOBA and NSGA-II. The experimental results indicate that there is a competitive relationship between the two targets in wet, normal and dry years. Compared to the MOBA and NSGA-II, the CAMOBA


Introduction
In recent years, with the increasing awareness of environmental protection and sustainable development, renewable energy has attracted more and more people's attention [1]. Hydropower generated by fast running or falling water has become an important part of renewable energy [2,3]. The total generation of hydropower in China exceeded 1 trillion kWh by the end of 2015, while the total installed capacity of hydropower reached 320 million kilowatts [4]. Therefore, the operation of cascade hydropower stations is becoming an important issue in the optimization of power systems in China. Generally, properly increasing the output of hydropower stations in dry seasons is of great significance to meet the load balance of power systems and improve the safety of power supplies. However, increasing the output in dry seasons will cause the water level to drop faster. The loss of water head will reduce the total annual generating capacity of the power station. Therefore, in order to meet actual requirements of the cascade hydropower reservoirs operation (CHRO) problem, this paper establishes the reservoir operation model of cascade hydropower stations with consideration of annual power generation benefit and firm power output. the comparison between results of the proposed method and other algorithms. In the end, Section 6 outlines the conclusions of this work.

Objective Function
A multi-objective model is developed to optimize the CHRO problem that needs to meet two objective functions during scheduling.
(1) Maximizing the annual power generation of the hydropower system: (2) Maximizing the firm power output: where N k,t is the power output of the kth reservoir at tth period; L k is the synthetic output coefficient of the kth reservoir; H k,t is the hydraulic head of the kth reservoir at the tth period (m); Q LEC k,t is the generation flow of the kth reservoir at the tth period (m 3 /s); ∆t is the time duration of a single period.

Constraints
Objective functions are subject to the following constraints when solving the CHRO problem.
(1) Reservoir water balance equation: SV k,t = SV k,t−1 + (I k,t − Q OUT k,t ) × ∆t k = 1, 2, · · · , K, t = 1, 2, · · · , T (2) Reservoir storage conversion: SV k,t = sv k (Z k,t ) k = 1, 2, · · · , K, t = 1, 2, · · · , T + 1 (4) (3) Reservoir water head: Z k,t = sv −1 k SV k,t + SV k,t−1 /2 k = 1, 2, · · · , K, t = 1, 2, · · · , T Z down k,t = z k Q OUT k,t k = 1, 2, · · · , K, t = 1, 2, · · · , T H k,t = Z k,t − Z down k,t k = 1, 2, · · · , K, t = 1, 2, · · · , T (4) Reservoir water level constraint: Z min k,t ≤ Z k,t ≤ Z max k,t k = 1, 2, · · · , K, t = 1, 2, · · · , T+1 (5) Reservoir outflow constraint: (6) Power output constraint: where SV k,t and SV k,t−1 are storage capacity of the kth reservoir at the tth and (t−1)th period, respectively (m 3 ); I k,t and Q OUT k,t are inflow and outflow of the kth reservoir at the tth period, respectively (m 3 /s); Z k,t is the water level of the kth reservoir at the tth period (m); sv k (Z k,t ) is the storage-capacity curve of the kth reservoir; z k Q OUT k,t is the function between the outflow and tail water level of the kth reservoir; Z down k,t is the tail water level of the kth reservoir (m); Z min k,t and Z max k,t are the lowest and highest water level of the kth reservoir at the tth period, respectively (m); Q OUT k,t is the outflow of the kth reservoir at the tth period (m 3 /s); Q S k,t is the spillage of the kth reservoir at the tth period (m 3 /s); Q min k,t and Q max k,t are minimum and maximum outflow of the kth reservoir at the tth period, respectively (m 3 /s); q k (Z k,t ) is the function between the water level and maximum outflow of the kth reservoir; N max k,t is the maximum power output of the kth reservoir at the tth period.

Overview of the Multi-Objective Bat Algorithm
The bat algorithm (BA) is a new heuristic algorithm based on the microbat echolocation behavior and proposed by Yang in 2010 [12]. In order to simplify and facilitate application, the algorithm adopts three idealized conditions: (1) The bat uses the echolocation principle to confirm the distance and accurately distinguish between obstacles and prey; (2) bats fly with speed V j and a fixed frequency F min to search for prey at coordinate position (solution) X j . They can flexibly amend the frequency F and adjust the pulse emission rate r j based on their distance from the target; (3) assume the loudness A j changing from a maximum value to a small value.
In the BA, new solutions and velocity updates can be obtained from Equations (13) and (14).
where β ∈ [0, 1] is a random vector; j is the number of bats; g is the iteration number; X * is the current global best location; F ∈ [F min , F max ] is the frequency; V g j is the velocity of the jth bat in the gth generation; X g j is the position (solution) of the jth bat in the gth generation. For each bat, if rand > r j , the new solution X g j,new is generated around the current best solution to complete the local search according to the Equation (15).
where X g j,old is the best solution of the jth bat in the first to gth generation; X g j,new is the new solution of the jth bat in the gth generation; ε ∈ [−1, 1] is a random vector; A g is the mean loudness of all the bats in the gth generation.
If rand > A j and the fitness of X g j,new is better than that of X * , X * is replaced by X g j,new . Then reduce A j and increase r j according to Equations (16) and (17). where α and γ are constants; A g+1 j is the loudness of the jth bat in g+1th generation; r g+1 j is the pulse emission rate of the jth bat in the g+1th generation; g is the iteration number.
In 2011, Yang proposed the MOBA algorithm based on BA to solve multi-objective optimization problems. The weighted sum is used to combine all objectives f m into a single objective in MOBA according to Equation (18) [17].
where ω m is the weight of the mth objective. All weights are generated randomly from the uniform distribution (0, 1).

Implementation of CAMOBA for solving the CHRO problem
Assuming that the D dimensional real space is the search space of the optimization problem, the algorithm is a generation population Where j is the number of individuals in the population, j = 1, 2, · · · , NG; g is the iteration number, g = 1, 2, · · · , g max . In the CHRO problem, D is the number of operation periods. The water levels of all reservoirs are chosen as the decision variable and encoded. The solution X g j is shown as follows.
where Z g k,t is the water level of the kth reservoir at the tth period in the gth generation (m).

External Archive Set Maintenance and Updating Operation
Zitzler [22] proposed an external population (called an archive set) in 1999 to preserve the non-dominant solutions found during evolution. The external archive set is added and named EAS( ) in this paper. The size of the external archive set is generally a constant N E due to the limitation of the computation source. The following strategies are used to incorporate all non-dominant solutions found in the current population g (called CSet (g)) into EAS( ): For each solution n g,j in CSet (g): (1) If EAS ( ) is vacant, put n g,j into EAS( ) directly.
(2) Compare n g,j with other solutions in EAS( ). If n g,j is not dominated by any solution in EAS ( ), add n g,j into EAS ( ) and delete the solution dominated by n g,j in EAS( ). are given a larger distance value to maintain the diversity of the external archive set [23]. The crowding distance can be expressed as: where Crowd j is the crowding distance of the jth solution; Fit m j is the fitness of the jth solution of the mth objective; min{Fit m } and max{Fit m } are the minimum and maximum fitness of the mth objective, respectively.

Initial Population Generation
The diversity and distribution of the initial population can influence the final optimal solutions of the algorithm. In MOBA, the initial population is usually generated randomly, which may lead to repeated solutions occupying memory space, and initial solutions may be concentrated in a limited space ( Figure 1a). In recent years, chaotic sequences have gradually replaced random sequences and achieved good results in many applications due to nonlinear phenomena such as ergodicity, randomness and certainty of chaos [24]. Thus, the initial population is generated based on chaos theory to enhance distribution diversity and uniformity of the population (Figure 1b). There are many models to generate chaotic sequences, and the improved logical mapping function [25] is adopted in this paper to generate the chaotic sequences. Its equation is as follows: where y q d is a chaotic variable and y q d ∈ (−1, 1). Equation (21) is chaotic state when y 1 d 0.
Water 2019, 12, x FOR PEER REVIEW 6 of 16 The diversity and distribution of the initial population can influence the final optimal solutions of the algorithm. In MOBA, the initial population is usually generated randomly, which may lead to repeated solutions occupying memory space, and initial solutions may be concentrated in a limited space (Figure 1(a)). In recent years, chaotic sequences have gradually replaced random sequences and achieved good results in many applications due to nonlinear phenomena such as ergodicity, randomness and certainty of chaos [24]. Thus, the initial population is generated based on chaos theory to enhance distribution diversity and uniformity of the population (Figure 1(b)). There are many models to generate chaotic sequences, and the improved logical mapping function [25] is adopted in this paper to generate the chaotic sequences. Its equation is as follows: where max

Self-Adaptive Local Search Strategy
Due to its properties of stable and randomness tendency, the normal cloud model [26] is integrated into CAMOBA to maintain the diversity of solutions. Each cloud is determined by three main parameters: Entropy parameter n E , expectation parameter x E and hyper entropy parameter Generate D different initial values, and then generate D chaotic sequences with different trajectories by iteration using Equation (21). Convert chaotic variables y q d to the value interval of the decision variables [x min d , x max d ], and initial positions x 1 j,d of the bat algorithm are generated by Equation (22).
where x max d and x min d are the maximum and minimum value of the dth decision variable, respectively; x 1 j,d is the initial position of the jth bat in the dth dimension.

Self-Adaptive Local Search Strategy
Due to its properties of stable and randomness tendency, the normal cloud model [26] is integrated into CAMOBA to maintain the diversity of solutions. Each cloud is determined by three main parameters: Entropy parameter E n , expectation parameter E x and hyper entropy parameter H e . The cloud operator L(E x , E n , He) of the cloud model for the current best solution is formulated as follows: (23) where, E x , E n and H e are the expectation value, entropy and hyper-entropy, respectively.
Steps of cloud drops generation are given as follows: Step 1: Generate a random normal distribution value E n with the expected value of E n and hyper entropy H e .
Step 2: Generate a random normal distribution value x p with the expected value of E x and hyper entropy E n generated in step 1.
Step 4: One drop x p , u p is attained. Steps 1 to 3 are repeated until enough cloud drops are generated.
Step 5: Each cloud drop represents a potential solution. If the newly generated cloud drop (bat) dominates the current best solution X g j,old , the X g j,old will be replaced. The self-adaptive local search strategy is described as in Equation (24). In order to judge whether the bat individual is in a stable state, two variables stable (g) and move (g) are introduced. The stable state is defined as: Bat individual j satisfies stable(g) < 10 −6 & move(g) < 10 −6 in successive 0.05g max generations.
where, f X g j,old is the fitness of bat individual j; f X * j,old is the optimal fitness of bat individual j in iteration; stable(g) is the gap between the current fitness of the jth bat and the optimal fitness in iteration; move(g) is the fitness increment of the jth bat in the gth generation; g max is the maximum iteration.

Mutation Operation
In the multi-objective bat algorithm (MOBA), the decreased diversity may cause precocity convergence (Figure 1a). Therefore, the mutation operation, as shown in Figure 1b, is proposed to diversify bats and enhance the ability of the global search, which can be expressed as: where, X g R 1 , EAS , X g R 2 , EAS and X g R 3 , EAS are three archive individuals randomly chosen from the external archival set EAS( ). δ is the mutation constant that is used to regulate the size of the disturbance and enhance the performance of the bat. Then, the bat X g j, EAS is compared with other solutions in EAS( ) to update the EAS( ).

Procedures of CAMOBA for Solving the CHRO Problem
The flowchart of CAMOBA for solving the CHRO problem is presented in Figure 2.

Case Study Description
The three hydropower stations on the Qingjiang River in China are used to verify the practical

Case Study Description
The three hydropower stations on the Qingjiang River in China are used to verify the practical feasibility of CAMOBA for solving the CHRO problem. The Qingjiang River is the first main tributary of the Yangtze River below the Three Gorges. The 423 km long Qingjiang River has a catchment area of 1.7 × 10 3 km 2 . The Qingjiang cascade hydropower project plays an important role in promoting social and economic development of the Yangtze River with multiple functions, including hydropower generation, flood control, ecological protection, etc. Figure 3    The hydrological year can be divided into three types, wet, normal and dry years, according to the annual runoff [27]. According to the runoff data of the Qingjiang River basin from 1971 to 2006, the typical representatives of June 1995 to May 1996, June 2004 to May 2005 and June 1985 to May 1986 are selected as the wet, normal and dry years, respectively. The whole operation time is divided into 12 periods, each of which is one month. The monthly inflow data of the Qingjiang cascade reservoirs in these three typical years are shown in Figure 4.

Parameter Settings
To testify the feasibility and effectiveness of CAMOBA in the application of the CHRO problem, CAMOBA is compared with MOBA [17] and NSGA-II [24]. NSGA-II is a

Parameter Settings
To testify the feasibility and effectiveness of CAMOBA in the application of the CHRO problem, CAMOBA is compared with MOBA [17] and NSGA-II [24]. NSGA-II is a high-performance multi-objective evolutionary algorithm that has been widely used in various disciplines [28]. The three optimization methods perform 10 independent runs. The maximum iteration g max is set at 1000 for the three algorithms. For CAMOBA and NSGA-II, the population size NG and the external archive set size N E are set at 200 and 30. For MOBA, the population size NG is 30. The maximum iteration of the chaos operator q max is 200 in CAMOBA. The mutation constant δ is 0.1 in CAMOBA. The other parameters of CAMOBA and MOBA are taken from reference [17]. The frequency F is varied from 0 to 1; α and γ are set at 0.9. The parameters of NSGA-II take the recommended values specified in reference [29]; the probabilities of crossover and mutation are set at 0.8 and 0.33, respectively.

Results and Discussion
The Pareto optimal fronts of different hydrological years calculated by the three methods are drawn in Figure 5. It can be seen from Figure 5 that the CHRO problem can be solved by CAMOBA, MOBA and NSGA-II. Meanwhile, it can be clearly observed that there is a competitive relationship between two targets in wet, normal and dry years, and the greater the annual power generation, the smaller the firm power output. The annual power generation in a wet year is larger than that in a normal year, and the annual power generation in a dry year is the smallest. Moreover, it can be easily concluded that when solving the complex CHRO problem by NSGA-II and MOBA, the resulting Pareto solution is dominated by CAMOBA, which shows that CAMOBA can generate more annual power generation with the same firm power output. For example, when the maximal firm power output is considered, the proposed CAMOBA method can increase the firm power output by about 21.1423 MW and 32.5989 MW, respectively, compared to the MOBA and NSGA-II in a normal year while increasing the annual power generation 0.0668 × 10 8 kWh and 0.1551 × 10 8 kWh, respectively. Meanwhile, it can be seen that CAMOBA has a more widely distributed optimal solution than MOBA and NSGA-II, which means that CAMOBA performs better in terms of solution diversity. The maximum (Max), mean (Mean) and standard deviation (Std) of the maximum annual power generation and firm power output are listed in Table 1. As shown in Table 1, it can be clearly seen that solutions of CAMOBA are close to the best solution, while all the indexes of CAMOBA are better than those of MOBA and NSGA-II. For example, the standard deviations of CAMOBA are 0.0173 and 29.5803 with respect to two objectives and less than that of MOBA and NSGA-II in a normal year. Thus, the results show that CAMOBA can provide better scheduling solutions compared to MOBA and NSGA-II when dealing with the complex CHRO issue.
in a normal year while increasing the annual power generation 0.0668 × 10 8 kWh and 0.1551 × 10 8 kWh, respectively. Meanwhile, it can be seen that CAMOBA has a more widely distributed optimal solution than MOBA and NSGA-II, which means that CAMOBA performs better in terms of solution diversity. The maximum (Max), mean (Mean) and standard deviation (Std) of the maximum annual power generation and firm power output are listed in Table 1. As shown in Table  1, it can be clearly seen that solutions of CAMOBA are close to the best solution, while all the indexes of CAMOBA are better than those of MOBA and NSGA-II. For example, the standard deviations of CAMOBA are 0.0173 and 29.5803 with respect to two objectives and less than that of MOBA and NSGA-II in a normal year. Thus, the results show that CAMOBA can provide better scheduling solutions compared to MOBA and NSGA-II when dealing with the complex CHRO issue.   The spacing (SP) [30] and hypervolume (HV) [31] metrics are used to further compare the performance in diversity and convergence of the proposed CAMOBA method with MOBA and NSGA-II. The three optimization methods perform 10 independent runs. The maximum (Max), minimum (Min), mean (Mean), standard deviation (Std) and average execution times of SP and HV are listed in Table 2, and best results are indicated in bold type. As the average computational time of the three algorithms listed in Table 2 show, CAMOBA is faster than NSGA-II but slightly slower than MOBA. CAMOBA consumes more computational time than MOBA because of the addition of the external archive set and the mutation operation to the proposed CAMOBA method. However, the primary purpose of algorithms is to obtain better Pareto optimal fronts. The running time of CAMOBA in Table 2 is reasonable and acceptable. Besides, the maximum and minimum SP results of NSGA-II are slightly better than CAMOBA in a wet year. However, the mean and standard deviation SP results of CAMOBA are significantly better than MOBA and NSGA-II in the three hydrological years. The observations indicate that CAMOBAs are superior to MOBA and NSGA-II in most of the runs in terms of diversity. Meanwhile, the HV results of CAMOBA are significantly better than those of MOBA and NSGA-II in all typical years, indicating that CAMOBA is superior to the other two algorithms in terms of convergence. The SP and HV convergence curves of the three algorithms are provided in Figure 6. The whole convergence process and the detailed convergence trajectory in iteration 1 to 200 are shown in the first and second rows, respectively. As can be seen from Figure 6, CAMOBA is the first one that converges to a good level and remains stable in all hydrological years. Therefore, the feasibility and superiority of the proposed CAMOBA method for solving the multi-objective CHRO problem in wet, normal and dry years within a reasonable execution time have been verified. For simplicity, only the results of the normal year are presented in detail. The reservoir operation processes of the 1st, 15th and 30th schemes (Table 3) obtained by CAMOBA are displays in Figure 7. From Figure 7, it can be seen that there is significant difference between the three scheduling schemes. The main focus of the 1st scheme is maximizing annual power generation. The monthly power generation of the 1st scheme varies greatly in a year. In the dry season (from December to February), reservoirs in the 1st scheme reduce outflows to reserve water, and the power output is lower than other schemes at this time. When inflow is large in March and April, a large amount of water is discharged from reservoirs in the 1st scheme to generate more power generation than other schemes. The 30th scheme focuses on maximizing the firm power output. Reservoirs in the 30th scheme discharge more water than other schemes in the dry season to generate more power output and improve the power supply safety of the power system. It can be seen that the monthly power generation of the 30th scheme has the smallest variation in a year. The 15th scheme is a compromise between annual power generation and firm power output. It is conducive to the balanced development of economic benefits and power supply security.  Note:↓denotes that the smaller value is better.↑means that the bigger value is better.

Conclusions
The CHRO is a complicated multi-objective optimization problem that simultaneously considers annual hydropower generation and firm power output. In this paper, CAMOBA has been proposed to handle the CHRO problem. In CAMOBA, the external archive set is added to preserve the non-dominant solutions. The population initialization based on chaos theory is adopted to improve population diversity. The self-adaptive local search strategy based on the normal cloud model is proposed to update solutions. The mutation operation is designed to mitigate premature convergence. Finally, CAMOBA is applied to the CHRO problem of the Qingjiang cascade hydropower stations in southern China. The case verifies the validity and feasibility of the proposed CAMOBA method. Its superiority in convergence and diversity is verified by comparing the results with other algorithms, including MOBA and NSGA-II. The experimental results indicate that there is a competitive relationship between the two targets in wet, normal and dry years. Compared to

Conclusions
The CHRO is a complicated multi-objective optimization problem that simultaneously considers annual hydropower generation and firm power output. In this paper, CAMOBA has been proposed to handle the CHRO problem. In CAMOBA, the external archive set is added to preserve the non-dominant solutions. The population initialization based on chaos theory is adopted to improve population diversity. The self-adaptive local search strategy based on the normal cloud model is proposed to update solutions. The mutation operation is designed to mitigate premature convergence. Finally, CAMOBA is applied to the CHRO problem of the Qingjiang cascade hydropower stations in southern China. The case verifies the validity and feasibility of the proposed CAMOBA method. Its superiority in convergence and diversity is verified by comparing the results with other algorithms, including MOBA and NSGA-II. The experimental results indicate that there is a competitive relationship between the two targets in wet, normal and dry years. Compared to the MOBA and NSGA-II, the CAMOBA method proposed in this paper can generate more annual power generation with the same firm power output within a reasonable computational time, so it provides a new approach for solving the CHRO problem.