## Abstract

Reservoir scheduling based on evolutionary algorithms needs to handle potentially stringent physical and operational constraints. Both generic and reservoir scheduling problem-specific constraint-handling techniques (CHTs) have their own merits and limitations. No CHT currently available can yield better solutions than the others consistently. To ensure good reservoir operation schedules, we develop an ensemble of CHTs (ECHT) that can utilize the advantages of different individual CHTs. In the ensemble, each CHT has its own population. In every generation, the different offspring populations are mixed together and evaluated. Each CHT then assigns fitness to all individuals and selects some of them to form its new parent population. The ECHT has been tested against long-term hydropower scheduling of two large-scale reservoir systems in China. Results show that the ECHT outperforms the state-of-the-art CHTs, and its probability of returning feasible solutions is much higher. The reservoir levels optimized with the ECHT are well suited for hydropower generation, which also reduce the chance of reservoir spilling.

## INTRODUCTION

Efficient reservoir management is necessary to obtain more economic, social and environmental benefits (Banos *et al.* 2011). Reservoir scheduling involves the optimization of a problem with one or more objectives, which is often dynamic, high-dimensional, multi-stage, nonlinear and nonconvex (Labadie 2004). The optimization is subject to a range of physical and operational constraints, such as limitations on reservoir storage, water release and power output, as well as the conservation of mass. For a multi-reservoir system, the hydraulic connection between adjacent reservoirs and the water transport delay further complicate the problem. A variety of optimization techniques has been applied to reservoir scheduling (Yeh 1985). Among the commonly used techniques, evolutionary algorithms (EAs) have gained increasing attention mainly due to their general applicability (Schardong & Simonovic 2015). For instance, the problems to be optimized do not have to be temporally separable, linear or differentiable (Labadie 2004). Even though EAs cannot guarantee to find the global optima, they, in most cases, can yield encouraging optimization results. In order to tackle constraints of reservoir scheduling problems, EAs have to work with constraint-handling techniques (CHTs), and some of the EAs have built-in CHTs (e.g., Poláková 2017; Trivedi *et al.* 2017; Zamuda 2017).

CHTs used with EAs for single-objective constrained optimization problems can be broadly classified into four categories: (1) penalty function techniques, (2) techniques preferring feasible individuals, (3) multi-objective optimization techniques and (4) techniques handling feasible and infeasible individuals separately (Coello 2002; Mezura-Montes & Coello 2011). In general, the penalty function techniques add an additional term to the original objective function to penalize infeasible individuals for violating any constraint. Although canonical penalty functions are easy to implement, tedious parameter tuning is necessary to obtain reasonable penalty coefficients (Runarsson & Yao 2000). By contrast, parameter-less adaptive penalties adjust themselves based on the information exploited from the population. For the second category, feasible individuals are simply preferred over infeasible ones. In some cases, constraint violation below a certain threshold is ignored at the beginning of the search and the threshold gradually decreases to zero in the search process. These CHTs have essentially no free parameter to adjust. Nevertheless, an additional mechanism for maintaining population diversity is required to avoid local optima (Coello 2002). The third category treats constraint satisfaction as an independent objective and uses multi-objective optimization concepts to tackle the problem. These techniques avoid parameter tuning and achieve a balance between objective optimization and constraint satisfaction. They can, however, be computationally inefficient (Runarsson & Yao 2005). CHTs in the fourth category handle feasible and infeasible individuals differently, mainly because comparing individuals of the same kind is relatively straightforward. The separately evaluated individuals are then merged by considering the search stage and population composition. Such techniques can incorporate various strategies for fitness assignment conveniently.

In the field of reservoir scheduling, different types of CHTs have been successfully applied (Nicklow *et al.* 2010). Penalty functions are most frequently used due to their simplicity (e.g., Haddad *et al.* 2015; Zhang *et al.* 2016; Afshar & Hajiabadi 2018), followed by repair methods that modify infeasible individuals into feasible ones (e.g., Niu *et al.* 2018; Wang *et al.* 2018) and methods that avoid exploring the infeasible space completely (e.g., Afshar & Moeini 2008; Moeini & Babaei 2017). However, penalty functions involve multiple free parameters to adjust; the other two types are useful only when the constraints are simple so that feasible individuals can be easily obtained or infeasible space can be easily separated from feasible space. More recently, a CHT adopting new adaptive strategies for fitness assignment was developed by the authors for solving reservoir scheduling problems (Hu *et al.* 2018b). This CHT (described later) showed encouraging performance compared with the existing ones. However, it cannot guarantee a feasible solution under extremely dry conditions, meaning that there is still room for improvement in the strategy of locating feasible space.

It is obvious that different CHTs have their respective merits and limitations. No CHT can outperform all others on problems of different types, and some can be more effective at certain stages of the search. In this context, an ensemble of CHTs (ECHT) was proposed to utilize the advantages of different techniques (Mallipeddi & Suganthan 2010). In this approach, each CHT in the ensemble maintains its own population and the different populations are merged and rebuilt in every generation (Mallipeddi & Suganthan 2010). The different CHTs assign fitness to the merged population separately, meaning promising individuals carrying key information to guide the search direction can be better recognized. In addition, the processes of locating feasible space and approaching global optima would accelerate. The ECHT has been reported to have better performance than the individual CHTs from the ensemble (Mallipeddi & Suganthan 2010). Apart from CHTs, the ensemble strategy has also been used to improve the performance of evolutionary operators and EAs (Mallipeddi *et al.* 2011; Naeini *et al.* 2018).

To the best of our knowledge, the method of CHT ensemble has not been reported in the field of water resources management. This study, inspired by the work of Mallipeddi & Suganthan (2010), develops a new ECHT to tackle the potentially stringent constraints in reservoir scheduling problems. The ECHT takes advantage of different individual CHTs' merits so that it can better guide the search direction and help make better reservoir operation schedules, especially under dry conditions. Multiple test cases for ECHT evaluation are designed based on different reservoir systems and inflow conditions. Four individual CHTs are also applied for performance comparison.

## PROBLEM FORMULATION

This study adopts the long-term hydropower scheduling as an example to demonstrate the applicability of the proposed ECHT in solving reservoir scheduling problems.

The scheduling of a reservoir system can be considered as a two-step process. The long-term scheduling on a monthly basis is conducted in the first place. Its results are then used as the boundary for more detailed scheduling (e.g., daily or weekly) (Labadie 2004). The start time of long-term scheduling can change every month with reservoir system status and inflow data updated. Future reservoir inflow over a long time span can be obtained based on long-term runoff forecast.

Long-term hydropower scheduling of a multi-reservoir system is a constrained optimization problem. It attempts to optimize water releases over time to boost hydroelectricity (HE) generation. Meanwhile, various constraints on reservoir operation must be satisfied. For a general hydropower scheduling problem, box constraints on reservoir forebay elevation, outflow discharge and power output are inevitable. However, other types of constraints can be eliminated by choosing proper decision variables. Long-term hydropower scheduling in this study is formulated as follows.

### Objective function

*HE*is described as: where

*M*is the number of reservoirs being considered,

*T*is the number of scheduling periods, is the power output of reservoir

*i*in period

*t*, and Δ

*t*is the time step. Reservoir

*i*is located upstream of reservoir

*i*+ 1. where is the conversion factor of reservoir

*i*in period

*t*, is the discharge through the turbines of reservoir

*i*for electricity generation in period

*t*and , and are the average forebay elevation, tailwater elevation and average head loss, respectively, of reservoir

*i*in period

*t*. , where is the forebay elevation of reservoir

*i*at the beginning of period

*t*.

*M*reservoirs at the beginning of

*T*periods: . The outflow discharges are obtained based on the water balance equation: where is the storage of reservoir

*i*at the beginning of period

*t*and , , , and are the inflow discharge, outflow discharge, flow discharge taken out for agricultural and industrial uses, discharge of reservoir storage loss (e.g., evaporation and seepage) and spilled discharge, respectively, of reservoir

*i*in period

*t*. The spilled discharge is part of the outflow that exceeds the maximum discharge through reservoir turbines. As can be seen, the water balance equation no longer needs to be tackled as an independent equality constraint owing to decision variable selection.

*i*in period

*t*.

### Constraints

Reservoir operation is subject to multiple constraints. The most common ones are listed as follows:

- (1)
- (2)Outflow discharge constraint: where is the lower bound on the outflow discharge of reservoir
*i*in period*t*that meets the demands of downstream navigation, water supply, irrigation, water quality control, etc. is the corresponding upper bound, which is limited by both reservoir discharge capacity and downstream flood control requirements. - (3)

Figure 1 shows the calculation procedure of the key variables mentioned above. With forebay elevation , reservoir storage is estimated based on the forebay elevation–reservoir storage relationship. Then outflow discharge is calculated based on water balance equation (Equations (3) and (4)). By information of the maximum turbine flow rate, discharge for electricity generation and spilled discharge are determined. With forebay elevation , average forebay elevation is clear. With outflow discharge , tailwater elevation is calculated based on outflow discharge–tailwater elevation relationship. With , and , power output and *HE* are obtained (Equations (1) and (2)). According to Equations (6)–(8), violation amount of each constraint is then calculated.

## METHODS FOR CONSTRAINT HANDLING

In this section, we first introduce the realization of the ECHT and then describe four state-of-the-art CHTs briefly, namely, the nondomination rank-based adaptive method (NRAM) devised by Hu *et al.* (2018b), the balanced ranking method (BRM) by Rodrigues *et al.* (2016), the stochastic ranking (SR) by Runarsson & Yao (2000) and the adaptive penalty function (APF) by Tessema & Yen (2009).

The four individual CHTs involve different ideas of how to tackle the constraints of an optimization problem. The NRAM evaluates feasible and infeasible individuals differently and involves non-dominated sorting in infeasible individual evaluation. It thus belongs to the third and fourth categories of CHTs at the same time. The BRM ranks different types of individuals separately and it falls into the fourth category. The idea of SR is so unique that it does not resemble any pre-defined CHT category. As for the APF, it clearly falls into the first CHT category.

The NRAM is a newly devised CHT for solving reservoir scheduling problems. The technique is free of parameter tuning and independent of the problem being optimized and the EA being used. The NRAM shows encouraging performance compared with existing CHTs, yet its computational efficiency is relatively low (Hu *et al.* 2018b). The BRM is a relatively new CHT, which is also parameter-less and independent of any particular EA. This technique is easy to use and able to provide satisfactory optimization results (de Paula Garcia *et al.* 2017). The SR is straightforward to implement and highly competitive with more sophisticated CHTs. This technique requires parameter tuning and its performance is sensitive to parameter value selection (Rodrigues *et al.* 2016). The APF is parameter-free, simple to implement and applicable to problems with very small feasible space. In some cases, however, the APF can penalize infeasible individuals insufficiently so that feasible individuals may not survive into the next generation (Hu *et al.* 2018b).

### Ensemble of CHTs (ECHT)

Similar to an individual CHT, the ECHT initially proposed by Mallipeddi & Suganthan (2010) assigns fitness to the population based on objective function values and constraint violations. This approach works at both stages of the search process, namely (1) the stage of locating feasible space and (2) the stage of converging to global optima. If all CHTs composing the ensemble assign fitness to the population as a whole, the ECHT can be coupled with EAs evaluating the population together, rather than those focusing on pair-wise comparisons.

Figure 2 presents the procedure of population evolution when an EA is coupled with the ECHT. An individual CHT in the ensemble (i.e., NRAM, BRM or SR as shown in Figure 2) has its own population. In every generation, offspring populations corresponding to the different CHTs are generated by applying evolutionary operators (e.g., selection, mutation and crossover) to the parent populations. Then, the offspring populations are mixed together and evaluated in terms of objective function value and constraint violations (see Equations (1)–(8)). Finally, each CHT assigns fitness values to all the individuals in the mixed population and selects the best third of them to make up a new population. In the next generation, the new populations of the different CHTs turn parent populations and the procedure described above runs again. The whole evolution process stops when the maximum number of generations is reached.

In order to fulfill the potential of the ECHT, CHTs composing the ensemble should have different theoretical foundations to achieve methodological diversity (Mallipeddi & Suganthan 2010). In addition, the CHTs need to have good performance in terms of the optimization results and the chance of returning feasible solutions. Three CHTs, i.e., NRAM, BRM and SR, are selected to compose the ensemble in this study. The APF, however, is ignored due to its relatively poor performance in our previous study (Hu *et al.* 2018b).

The goal of Equation (9) is to find the optimal **x** that minimizes the objective function and satisfies *q* inequality constraints and *m–q* equality constraints

### Nondomination rank-based adaptive method

Given that the constraints of reservoir scheduling problems are potentially stringent and the global optima often lie on the boundary between feasible and infeasible spaces, the NRAM emphasizes exploiting information from infeasible individuals and preserving them throughout the search. In addition, the population composition is dynamically adjusted to facilitate exploration and exploitation at the beginning and end of the search, respectively.

When infeasible space is explored, objective optimization and constraint satisfaction are deemed as bi-objective optimization with the objective vector , where is the normalized objective function value of the infeasible individual defined in a similar way as in Equation (14). In the NRAM, infeasible individuals are compared mainly according to the nondomination rank obtained based on .

*η*) is introduced to facilitate exploration or exploitation by adjusting the population composition. The fitness of the infeasible individual is finally formulated as:

### Balanced ranking method

The BRM assigns fitness by ranking the population. This technique involves a separate ranking of feasible and infeasible individuals and merging of all ranked individuals afterward.

**x**among IP according to , is the number of infeasible individuals in the current population.

According to Equations (17)–(19), if the current population violates most of the constraints ( has a small value), objective function value does not have much influence over the penalty value of an infeasible individual. However, when a majority of the individuals are feasible ( has a large value), infeasible individuals with favorable objective function values will gain an edge in penalty function evaluation.

### Stochastic ranking

The SR assigns fitness by sorting the population as in the BRM. This technique was proposed to balance the dominance of objective function value and degree of constraint violation in the comparison of individuals. The SR sorts the current population using a bubble-sort-like procedure. The probability of comparing two adjacent individuals according to objective function value is 1 if both individuals are feasible; otherwise, a user-defined probability *P*_{f}.

### Adaptive penalty function

and are two terms used to penalize infeasible individuals for constraint violation. gives priority to infeasible individuals having low constraint violations (favorable objective function values) when the current population has a small (large) feasibility ratio.

## CASE STUDY

We applied the ECHT to the hydropower scheduling of two large-scale reservoir systems in China: the Three Gorges Reservoir-Gezhouba Reservoir (TGR-GR) system and the Shuibuya Reservoir-Geheyan Reservoir (SBYR-GHYR) system. Multiple test cases were designed based on the two reservoir systems under different inflow conditions.

### Study area

The aforementioned four reservoirs are located on the Yangtze River or its tributary, the Qing River (Figure 3). The Yangtze River is the largest river in China and the third largest in the world. It has a main channel of 6,300 km and a drainage area of 1.8 × 10^{6} km^{2} (Mao *et al.* 2016). Among the Yangtze River's tributaries located downstream of the TGR, the Qing River is one of the largest, which reaches a length of 423 km and drains an area of 1.76 × 10^{4} km^{2} (Zhou *et al.* 2016).

TGR is one of the world's largest water resources projects (Hu *et al.* 2018a). With a drainage area of 1.0 × 10^{6} km^{2}, the TGR has a mean annual runoff of 4.51 × 10^{11} m^{3} at the dam site (Li & Qiu 2016). Table 1 summarizes the main characteristics of the TGR and the other three reservoirs. The normal pool level of the TGR is 175 m, corresponding to a storage capacity of 3.93 × 10^{10} m^{3}. The TGR's forebay elevation is kept at the flood control water level (145 m) throughout the flood season, thereby vacating a room of 2.22 × 10^{10} m^{3} to store the floodwater temporarily. The installed capacity and firm output of the TGR hydropower system are 2.25 × 10^{7} kW and 4.99 × 10^{6} kW, respectively. The Gezhouba Dam, located 38 km downstream from the Three Gorges Dam, is a run-of-the-river dam. The GR has a storage capacity of 7.11 × 10^{8} m^{3} and a normal pool level of 66 m. However, in order to smooth the TGR outflow, the GR forebay elevation is allowed to fluctuate between 63.0 m and 66.5 m within a day (Li *et al.* 2013).

Parameter | TGR | GR | SBYR | GHYR |
---|---|---|---|---|

Location | Yangtze River | Yangtze River | Qing River | Qing River |

Gross storage (m^{3}) | 3.93 × 10^{10} | 7.11 × 10^{8} | 4.31 × 10^{9} | 3.12 × 10^{9} |

Flood control storage (m^{3}) | 2.22 × 10^{10} | - | 5.00 × 10^{8} | 5.00 × 10^{8} |

Crest elevation of the dam (m) | 185 | 70 | 409 | 206 |

Normal pool level (m) | 175 | 66 | 400 | 200 |

Flood control water level (m) | 145 | - | 391.8 | 192.2 |

Installed capacity (kW) | 2.25 × 10^{7} | 2.74 × 10^{6} | 1.84 × 10^{6} | 1.21 × 10^{6} |

Firm output (kW) | 4.99 × 10^{6} | 1.04 × 10^{6} | 3.11 × 10^{5} | 1.87 × 10^{5} |

Parameter | TGR | GR | SBYR | GHYR |
---|---|---|---|---|

Location | Yangtze River | Yangtze River | Qing River | Qing River |

Gross storage (m^{3}) | 3.93 × 10^{10} | 7.11 × 10^{8} | 4.31 × 10^{9} | 3.12 × 10^{9} |

Flood control storage (m^{3}) | 2.22 × 10^{10} | - | 5.00 × 10^{8} | 5.00 × 10^{8} |

Crest elevation of the dam (m) | 185 | 70 | 409 | 206 |

Normal pool level (m) | 175 | 66 | 400 | 200 |

Flood control water level (m) | 145 | - | 391.8 | 192.2 |

Installed capacity (kW) | 2.25 × 10^{7} | 2.74 × 10^{6} | 1.84 × 10^{6} | 1.21 × 10^{6} |

Firm output (kW) | 4.99 × 10^{6} | 1.04 × 10^{6} | 3.11 × 10^{5} | 1.87 × 10^{5} |

The Shuibuya Dam and Geheyan Dam on the Qing River have a distance of 92 km between each other. The SBYR has a contributing area of 1.09 × 10^{4} km^{2} and a mean annual inflow of 296 m^{3}/s. The reservoir has a storage capacity of 4.31 × 10^{9} m^{3} at the normal pool level (400 m), capable of multi-year regulation. The firm output of the SBYR hydropower plant is 3.11 × 10^{5} kW. For the GHYR, the drainage area and mean annual inflow are 1.44 × 10^{4} km^{2} and 403 m^{3}/s, respectively. Its normal pool level is 200 m, corresponding to a gross storage of 3.12 × 10^{9} m^{3}. GHYR's hydropower system has a firm output of 1.87 × 10^{5} kW. The flood control storage is 5.00 × 10^{8} m^{3} for both the SBYR and GHYR.

For both reservoir systems, the scheduling horizon was a year and the time step was one-third of a month. The lateral inflow to the GR was ignored due to the short distance between the Three Gorges Dam and Gezhouba Dam. For each reservoir, the final forebay elevation at the end of the scheduling was assumed to equal the initial forebay elevation. Such an assumption is somewhat unrealistic. However, it should not influence the evaluation of the ECHT in terms of the ability to locate feasible space and approach global optima.

### Test case design

Only reservoir inflows in dry conditions were used in the design of test cases. For a multi-reservoir system, it is difficult to allocate over time the limited water resources (current reservoir storage and future inflows) in dry years with all the constraints satisfied (Labadie 2004). Such a problem can challenge the effectiveness of a CHT in terms of whether the search can locate feasible space and then approach global optima as rapidly as possible.

Inflow data from 1901 to 2000 for TGR-GR system and data from 1975 to 1987 for SBYR-GHYR system were collected. For each reservoir system, historical inflow data in the five driest years were identified as five inflow scenarios. Note that each inflow scenario allowed at least one of the constraint-handling methods to find a feasible solution. Figure 4 presents the different TGR inflow scenarios (QT1–QT5), SBYR inflow scenarios (QS1–QS5) and the corresponding lateral inflows to the GHYR. Based on these inflow scenarios, 10 test cases (T1 to T5 and S1 to S5) for ECHT evaluation were designed and listed in Table 2.

Test case | Reservoir system | Inflow scenario |
---|---|---|

T1 | TGR-GR | QT1 |

T2 | TGR-GR | QT2 |

T3 | TGR-GR | QT3 |

T4 | TGR-GR | QT4 |

T5 | TGR-GR | QT5 |

S1 | SBYR-GHYR | QS1 |

S2 | SBYR-GHYR | QS2 |

S3 | SBYR-GHYR | QS3 |

S4 | SBYR-GHYR | QS4 |

S5 | SBYR-GHYR | QS5 |

Test case | Reservoir system | Inflow scenario |
---|---|---|

T1 | TGR-GR | QT1 |

T2 | TGR-GR | QT2 |

T3 | TGR-GR | QT3 |

T4 | TGR-GR | QT4 |

T5 | TGR-GR | QT5 |

S1 | SBYR-GHYR | QS1 |

S2 | SBYR-GHYR | QS2 |

S3 | SBYR-GHYR | QS3 |

S4 | SBYR-GHYR | QS4 |

S5 | SBYR-GHYR | QS5 |

For each test case, the feasible space accounts for a very small part of the search space; none of the 1 × 10^{8} randomly generated solutions satisfies all the constraints. The designed test cases can thus provide a testbed for the different methods for constraint handling.

### Experiment setup

Evolutionary operators of EA and CHTs both contribute to the optimization process. EA operators (e.g., mutation and crossover) generate new individuals in each generation, while the CHT assigns fitness to these individuals. Since this study focused on CHTs, a genetic algorithm (GA) with elementary evolutionary operators (Houck *et al.* 1998) was adopted as the search engine. This should, ideally, accentuate the differences between the ECHT and the individual CHTs.

In this study, GA population size and maximum number of generations both had large values so that the search could converge. The population size of GA was 210 when coupled with an individual CHT. As for the ECHT, GA population size was 210 as well (70 for each of the three CHTs in the ensemble). For TGR-GR test cases, the maximum number of generations was 1,500, while the corresponding value for SBYR-GHYR test cases was 10,000 as feasible space became much more difficult to find. The other GA parameter values were from our previous research (Hu *et al.* 2018b). In addition, the elitist strategy was used to preserve the best feasible individual found in the search process. For each constraint-handling method, 25 runs with random seeds were performed on every test case, which was useful for alleviating the uncertainty in optimization results.

The SR is the only CHT considered in this study that has a free parameter. For the present application, SR parameter *P*_{f} was set to 0.45, which lay within the range recommended by Runarsson & Yao (2000).

### Methods' performance comparison

To characterize the performance of a constraint-handling method, the maximum, mean and minimum values of the 25 optimization results on each test case were calculated. In addition, the ratio of runs that returned feasible solutions to all runs (i.e., feasible rate) and the number of generations capturing the first feasible individual were also considered. The last two indicators were used to reflect the probability and efficiency of a method to locate feasible space. They are of great importance especially in real-world applications, which require the search algorithm to yield satisfactory feasible solutions within a limited number of runs.

Table 3 summarizes the ECHT results on all designed test cases. According to the maximum, mean and minimum results, the ECHT outperforms all the individual CHTs on 10, 9 and 8 test cases, respectively. In addition, the optimization results of the ECHT are highly stable. The relative difference between the best and worst results is less than 0.15% among TGR-GR cases and 0.90% among SBYR-GHYR cases. On nine test cases out of 10, the mean result is greater than the average of the maximum and minimum results, indicating that the results of the ECHT tend to be distributed near the best result rather than the worst.

Test case | Maximum (10^{11} kW·h) | Mean (10^{11} kW·h) | Minimum (10^{11} kW·h) | Test case | Maximum (10^{9} kW·h) | Mean (10^{9} kW·h) | Minimum (10^{9} kW·h) |
---|---|---|---|---|---|---|---|

T1 | 1.016720^{a} | 1.016646^{a} | 1.016371^{a} | S1 | 9.400428^{a} | 9.392829^{a} | 9.383760^{a} |

T2 | 1.043562^{a} | 1.043315^{a} | 1.042801^{a} | S2 | 11.057987^{a} | 11.043160^{a} | 10.989475 |

T3 | 1.019856^{a} | 1.019674^{a} | 1.019500^{a} | S3 | 11.772772^{a} | 11.747224^{a} | 11.672634^{a} |

T4 | 1.068335^{a} | 1.067998^{a} | 1.067205^{a} | S4 | 7.876145^{a} | 7.858966^{a} | 7.828219^{a} |

T5 | 1.032983^{a} | 1.032760^{a} | 1.031664^{a} | S5 | 7.894330^{a} | 7.876874 | 7.841687 |

Test case | Maximum (10^{11} kW·h) | Mean (10^{11} kW·h) | Minimum (10^{11} kW·h) | Test case | Maximum (10^{9} kW·h) | Mean (10^{9} kW·h) | Minimum (10^{9} kW·h) |
---|---|---|---|---|---|---|---|

T1 | 1.016720^{a} | 1.016646^{a} | 1.016371^{a} | S1 | 9.400428^{a} | 9.392829^{a} | 9.383760^{a} |

T2 | 1.043562^{a} | 1.043315^{a} | 1.042801^{a} | S2 | 11.057987^{a} | 11.043160^{a} | 10.989475 |

T3 | 1.019856^{a} | 1.019674^{a} | 1.019500^{a} | S3 | 11.772772^{a} | 11.747224^{a} | 11.672634^{a} |

T4 | 1.068335^{a} | 1.067998^{a} | 1.067205^{a} | S4 | 7.876145^{a} | 7.858966^{a} | 7.828219^{a} |

T5 | 1.032983^{a} | 1.032760^{a} | 1.031664^{a} | S5 | 7.894330^{a} | 7.876874 | 7.841687 |

*Note*: ^{a} indicates that the ECHT outperforms all the individual CHTs.

The differences in optimized HE between the ECHT and the individual CHTs are shown in Figure 5. The ECHT appears to have a noticeable advantage over the other methods. The NRAM produces relatively good results among the four individual CHTs. For the BRM and SR, the former performs better in terms of the maximum result; the two CHTs become close when compared against the mean and minimum results. The APF yields feasible solutions on only two test cases and its results are shown to be far less competitive than the others. In terms of the best optimization result, the differences between the ECHT and the individual CHTs generally lie between 5.0 × 10^{6} kW h and 5.0 × 10^{8} kW h.

Figure 6 compares the feasible rates of the different methods for constraint handling. The average feasible rates of the ECHT, NRAM, BRM, SR and APF are 0.83, 0.72, 0.21, 0.27 and 0.05, respectively. On average, the ECHT has an advantage of 0.10, 0.62, 0.56 and 0.78 over the NRAM, BRM, SR and APF, respectively. The ECHT has a feasible rate of 1.00 on five test cases; the NRAM, second only to the ECHT, has a feasible rate of 1.00 on three. It is worth noting that the ECHT and NRAM are the only two methods that manage to return a feasible solution on every test case; the BRM, SR and APF fail on five, three and eight test cases, respectively. Figure 7 depicts generation numbers that return the first feasible individual across 25 runs of each method on every test case. The generation numbers of the ECHT are relatively small on most of the test cases, meaning this method can locate feasible space more efficiently. The SR and APF generally have the largest generation numbers among the five methods.

### Scheduled reservoir levels

The TGR and GR water levels corresponding to the greatest HE obtained with the ECHT are shown in Figure 8. The TGR levels rise rapidly from the flood control water level (i.e., 145 m) since mid-September and reach the normal pool level (i.e., 175 m) by the end of October. Generally, the reservoir levels remain unchanged until the next January. In the following periods, the TGR levels gradually decline in response to the decreased inflows. The reservoir levels would rise in late March and April due to the increased inflows. Afterward, they gradually fall back to 145 m by the end of early June. The TGR levels rise to or stay at the maximum water levels whenever possible and fall only when releasing the inflow cannot satisfy all the constraints. Such an operation strategy leads to large water heads that facilitate the TGR hydropower generation.

As shown in Figure 8(b), the GR levels lie on the upper bound throughout the scheduling horizon, in order to maximize the water heads. Releasing the GR inflows (identical to the TGR outflows in this study) is sufficient to satisfy all the underlying constraints. The GR has such a small gross storage that lowering the reservoir level goes against increasing hydropower generation, even though the corresponding outflow would increase slightly.

The SBYR and GHYR water levels corresponding to the greatest HE obtained with the ECHT are presented in Figure 9. According to Figure 9(a), the SBYR levels gradually decline to guarantee the firm output in the dry season (November to the next March). Then the water levels show rising trends in April due to the increased inflows. In the following periods, peak discharges will occur under several inflow scenarios. To avoid reservoir spilling, the SBYR levels decline first and then rise during peak discharges. By the end of July, the SBYR levels reach the upper bound (i.e., 391.8 m). Afterward, the water levels increase slowly or even decrease sometimes to release more water and raise the GHYR levels. When the SBYR levels reach the normal pool level (i.e., 400 m), they generally remain unchanged to maintain the water heads for a period of time and then decline gradually due to the decreased inflows.

Figure 9(b) shows that the GHYR water levels lie predominantly on the upper bound for large water heads, which is similar to the GR. In some cases, the GHYR lowers its forebay elevation before peak discharges to avoid reservoir spilling.

## DISCUSSION

The results from this study show that the developed ECHT outperforms four representative individual CHTs: the NRAM, BRM, SR and APF. The ensemble method is found to be capable of producing highly competitive solutions to hydropower scheduling problems with stringent constraints. More importantly, the ECHT has a much higher probability of returning feasible solutions than the simple CHTs, which is of great value in real-world applications.

The better performance of the ECHT can be explained in two aspects. First, promising individuals in the population play a more important role in the search process with the ECHT than with an individual CHT. In the ensemble method, after all offspring individuals are merged and evaluated, each CHT in the ensemble selects some individuals it prefers for the next generation. Promising individuals that carry key information for guiding the search direction can be better recognized, while the worst individuals are unlikely to survive into the next generation. Second, the ECHT exploits the advantages of different CHTs effectively. Under certain circumstances, each CHT can be more effective in terms of the probability and efficiency of locating feasible space, the ability to converge to global optima, etc. The CHTs in the ensemble assign fitness to the merged population separately, meaning they jointly contribute to the adjustments of the search direction.

It is also important to note that the advantages of the ECHT come at the cost of longer computational time since it spends more time assigning fitness to the population in each generation. However, GA coupled with the ECHT requires less computational time than the total time of GA coupled with each CHT in the ensemble since every function evaluation is effectively used in the ECHT.

This study used a large GA population size and a large maximum number of generations to fulfill the potential of the ECHT and individual CHTs. To further verify the better performance of the ECHT, we also adopted another set of parameter values, i.e., multiplying the original population size by 0.5 and the original maximum number of generations by 0.7. As a result, the number of function evaluations (NFE) became 65% smaller.

Figure 10 shows the feasible rates of different methods under a smaller NFE. The ECHT fails in returning any feasible solution on two test cases, while the NRAM, BRM, SR and APF fail on four, four, six and eight test cases, respectively. Compared with Figure 6, the decline in feasible rates of the ECHT, NRAM and SR most likely results from the smaller NFE. The average feasible rates of BRM and APF, however, do not change much. All the methods fail on S4 and S5 and the ECHT is the only method that returns feasible solutions on S1 and S3. On the rest of the test cases, the ECHT has a larger or equal feasible rate than the four individual CHTs. Figure 11 presents the differences in optimized HE between the ECHT and the individual CHTs under a smaller NFE. Without considering S4 and S5, the ECHT outperforms all the individual CHTs on seven out of eight test cases. The differences in maximum HE between the ECHT and the individual CHTs lie within the range of 6.5 × 10^{5} and 6.0 × 10^{8} kW h. Compared with Figure 5, the advantage of the ECHT in terms of the maximum result becomes larger especially over NRAM, SR and APF. The above results demonstrate the better performance of the ECHT under a smaller NFE.

In this study, the ECHT shows the best performance in terms of returning feasible solutions. The NRAM comes second, followed by the BRM and SR. The APF has the smallest feasible rate. The different individual CHTs adopt different strategies for infeasible individual evaluation before locating feasible space. The NRAM, BRM and SR use both objective function value and constraint violation, while the APF merely uses the latter information for fitness assignment. Even though the APF's strategy appears to be the most efficient, this study demonstrates that using objective function value in infeasible individual evaluation can facilitate the finding of feasible space. It can be attributed to the fact that constraint satisfaction does not conflict with objective function optimization, which results from the problem formulation in this study.

Figure 12 presents the relationships between objective function value and constraint violation from 2,000 randomly generated infeasible individuals for each test case. Objective function value shows an overall increasing trend when constraint violation decreases. It suggests that information of objective function value can help locate feasible space. NRAM, BRM and SR use such information in different ways. Objective function value and constraint violation are almost equally important in the NRAM. The only difference between them is that, for infeasible individuals with identical nondomination rank, those with lower constraint violations are preferred. The BRM and SR rank infeasible individuals with both types of information considered but constraint violation is more important. In summary, objective function value plays a more important role in the NRAM than in the BRM and SR, which may explain the larger feasible rate of the NRAM.

The development of the ECHT in this study considers both methodological diversity and including a CHT well suited for reservoir scheduling problems (i.e., the NRAM). Given the extremely small feasible space of the designed test cases, the NRAM should be included to improve the chance of returning feasible solutions, even though its computational time is longer than some other methods. On the other hand, when applied to simpler problems with shorter scheduling horizon or wetter hydrological scenario, the ECHT should highlight methodological diversity and avoid time-consuming CHTs so that reservoir operation schedules can be made in a shorter period of time.

This study clearly demonstrates that the ECHT can be used in place of an individual CHT to improve the optimization results of an EA. However, it should be noted that some other factors (e.g., optimization problem formulation, search algorithm selection and parametrization) may affect EA results as well (Clarkin *et al.* 2018). More research efforts can be put on the evaluation and comparison of contributions of the different factors.

## CONCLUSIONS

This study developed an ECHT to deal with the potentially stringent constraints of reservoir scheduling problems. The ECHT was tested against the long-term hydropower scheduling of two large-scale reservoir systems in China. Four state-of-the-art individual CHTs, i.e., the NRAM, BRM, SR and APF, were also applied for performance comparison. Based on the results and discussion, the following conclusions can be drawn: The ECHT outperforms the individual CHTs when applied to highly constrained reservoir scheduling problems. The ensemble method has a better chance of returning highly competitive feasible solutions, which is of great importance in real-world applications. In addition, GA coupled with the ECHT requires less computational time than the total time of GA coupled with each CHT in the ensemble. The reservoir levels optimized with the ECHT are well suited for hydropower generation and the probability of reservoir spilling is reduced.

## ACKNOWLEDGEMENT

The authors thank the anonymous reviewers for their constructive comments. This research is funded by the National Key Research and Development Program of China (2017YFC0405301, 2018YFC0407606), China Postdoctoral Science Foundation (2019M651890), the Central Public-interest Scientific Institution Basal Research Fund (Y519008) and the National Natural Science Foundation of China (51379059, 41601376).