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.

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.

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

The objective of maximizing HE is described as:
(1)
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.
(2)
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.
The decision vector consists of forebay elevations of M reservoirs at the beginning of T periods: . The outflow discharges are obtained based on the water balance equation:
(3)
(4)
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.
In a multi-reservoir system, the inflow of a reservoir is usually composed of the outflow of its upstream reservoir and the lateral inflow from the local catchment. In the case of a large time step (e.g., 10 days or a month), the water transport time can often be ignored. If this is the case, the reservoir inflow takes the following form:
(5)
where is the lateral inflow discharge into reservoir i in period t.

Constraints

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

  • (1)
    Forebay elevation constraint:
    (6)
    where and are the lower and upper bounds, respectively, on the forebay elevation of reservoir i at the beginning of period t.
  • (2)
    Outflow discharge constraint:
    (7)
    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)
    Power output constraint:
    (8)
    where and are the firm output and installed capacity of reservoir i, respectively.

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.

Figure 1

Function evaluation procedure of a hydropower scheduling problem.

Figure 1

Function evaluation procedure of a hydropower scheduling problem.

Close modal

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.

Figure 2

Flowchart of population evolution with the ECHT. PAR, parent population, OFF, offspring population.

Figure 2

Flowchart of population evolution with the ECHT. PAR, parent population, OFF, offspring population.

Close modal

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

To conveniently describe the CHTs concerned, the general concepts of constrained optimization are presented first. Without loss of generality, the minimization problem is considered here:
(9)

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

For the infeasible individual, the violation amount of the jth constraint is
(10)
where δ is the tolerance value (e.g., 0.001 or 0.0001).
The overall constraint violation is defined as
(11)
(12)

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.

The NRAM evaluates feasible individuals in the population according to the normalized objective function value:
(13)
(14)
where is the fitness, is the normalized objective function value, and are the minimum and maximum objective function values, respectively, of feasible individuals in the current population.

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 .

For infeasible individuals in the same nondomination front, those with low constraint violations and favorable objective function values are preferred in the case of small and large feasibility ratios, respectively. The modified nondomination rank is defined as:
(15)
where is the feasibility ratio of the current population, is the maximum nondomination rank of infeasible individuals in the current population.
At this point, a feasible (infeasible) individual in the population can be compared with other feasible (infeasible) individuals based on . When the two types of individuals are mixed together, the feasibility ratio is used to avoid abandoning infeasible individuals completely and to protect feasible individuals when they are rare. In addition, the iteration rate of the search (denoted η) is introduced to facilitate exploration or exploitation by adjusting the population composition. The fitness of the infeasible individual is finally formulated as:
(16)

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.

In the BRM, feasible individuals are ranked according to their objective function values, while infeasible individuals are ranked based on their penalty function values. For the infeasible individual, penalty function is defined as:
(17)
(18)
(19)
where is the number of constraints not violated by the current population, Npop is the population size, IP is the infeasible part of the current population, is the rank of 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.

After feasible and infeasible individuals are ranked separately, the infeasible queue is merged into the feasible queue. The fitness in the BRM is defined as:
(20)
(21)
(22)
(23)
where FP is the feasible part of the current population. The coefficient is introduced to push the infeasible queue toward the end of the merged queue, which gives priority to the best feasible individuals. The coefficient is responsible for expanding the infeasible queue within the merged queue.

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

Adaptive penalty function

The APF is a two-phase CHT that adopts different strategies to assign fitness at different stages of the search process. Generally, in the first search stage, all individuals in the population are infeasible. The first priority of the search is to locate a feasible space. Fitness at this stage is thus defined as:
(24)
Once feasible space is located (i.e., the second search stage), the fitness of the feasible individual equals the normalized objective function value :
(25)
(26)
where and are the minimum and maximum objective function values, respectively, of the current population.
For the infeasible individual, the fitness is formulated as:
(27)
(28)
(29)

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.

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 × 106 km2 (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 × 104 km2 (Zhou et al. 2016).

Figure 3

Map of the study area.

Figure 3

Map of the study area.

Close modal

TGR is one of the world's largest water resources projects (Hu et al. 2018a). With a drainage area of 1.0 × 106 km2, the TGR has a mean annual runoff of 4.51 × 1011 m3 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 × 1010 m3. 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 × 1010 m3 to store the floodwater temporarily. The installed capacity and firm output of the TGR hydropower system are 2.25 × 107 kW and 4.99 × 106 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 × 108 m3 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).

Table 1

Summary of main characteristics of the TGR, GR, SBYR and GHYR

ParameterTGRGRSBYRGHYR
Location Yangtze River Yangtze River Qing River Qing River 
Gross storage (m33.93 × 1010 7.11 × 108 4.31 × 109 3.12 × 109 
Flood control storage (m32.22 × 1010 5.00 × 108 5.00 × 108 
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 × 107 2.74 × 106 1.84 × 106 1.21 × 106 
Firm output (kW) 4.99 × 106 1.04 × 106 3.11 × 105 1.87 × 105 
ParameterTGRGRSBYRGHYR
Location Yangtze River Yangtze River Qing River Qing River 
Gross storage (m33.93 × 1010 7.11 × 108 4.31 × 109 3.12 × 109 
Flood control storage (m32.22 × 1010 5.00 × 108 5.00 × 108 
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 × 107 2.74 × 106 1.84 × 106 1.21 × 106 
Firm output (kW) 4.99 × 106 1.04 × 106 3.11 × 105 1.87 × 105 

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 × 104 km2 and a mean annual inflow of 296 m3/s. The reservoir has a storage capacity of 4.31 × 109 m3 at the normal pool level (400 m), capable of multi-year regulation. The firm output of the SBYR hydropower plant is 3.11 × 105 kW. For the GHYR, the drainage area and mean annual inflow are 1.44 × 104 km2 and 403 m3/s, respectively. Its normal pool level is 200 m, corresponding to a gross storage of 3.12 × 109 m3. GHYR's hydropower system has a firm output of 1.87 × 105 kW. The flood control storage is 5.00 × 108 m3 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.

Table 2

Test cases designed for ECHT evaluation

Test caseReservoir systemInflow 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 caseReservoir systemInflow 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 
Figure 4

Different scenarios of (a) the TGR inflow, (b) the SBYR inflow and (c) the corresponding lateral inflow to the GHYR.

Figure 4

Different scenarios of (a) the TGR inflow, (b) the SBYR inflow and (c) the corresponding lateral inflow to the GHYR.

Close modal

For each test case, the feasible space accounts for a very small part of the search space; none of the 1 × 108 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 Pf 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, nine and eight 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.

Table 3

Summary of the ECHT results

Test caseMaximum (1011 kW·h)Mean (1011 kW·h)Minimum (1011 kW·h)Test caseMaximum (109 kW·h)Mean (109 kW·h)Minimum (109 kW·h)
T1 1.016720a 1.016646a 1.016371a S1 9.400428a 9.392829a 9.383760a 
T2 1.043562a 1.043315a 1.042801a S2 11.057987a 11.043160a 10.989475 
T3 1.019856a 1.019674a 1.019500a S3 11.772772a 11.747224a 11.672634a 
T4 1.068335a 1.067998a 1.067205a S4 7.876145a 7.858966a 7.828219a 
T5 1.032983a 1.032760a 1.031664a S5 7.894330a 7.876874 7.841687 
Test caseMaximum (1011 kW·h)Mean (1011 kW·h)Minimum (1011 kW·h)Test caseMaximum (109 kW·h)Mean (109 kW·h)Minimum (109 kW·h)
T1 1.016720a 1.016646a 1.016371a S1 9.400428a 9.392829a 9.383760a 
T2 1.043562a 1.043315a 1.042801a S2 11.057987a 11.043160a 10.989475 
T3 1.019856a 1.019674a 1.019500a S3 11.772772a 11.747224a 11.672634a 
T4 1.068335a 1.067998a 1.067205a S4 7.876145a 7.858966a 7.828219a 
T5 1.032983a 1.032760a 1.031664a S5 7.894330a 7.876874 7.841687 

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 × 106 kW h and 5.0 × 108 kW h.

Figure 5

Optimized hydroelectricity of ECHT minus that of an individual CHT in terms of (a) maximum result, (b) mean result and (c) minimum result in 25 runs. The better results of NRAM on S5 in (b) and S2 and S5 in (c) are not shown. If a CHT does not return any feasible solution in all runs, there is nothing to be shown.

Figure 5

Optimized hydroelectricity of ECHT minus that of an individual CHT in terms of (a) maximum result, (b) mean result and (c) minimum result in 25 runs. The better results of NRAM on S5 in (b) and S2 and S5 in (c) are not shown. If a CHT does not return any feasible solution in all runs, there is nothing to be shown.

Close modal

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.

Figure 6

Bar plots of feasible rates of the different methods.

Figure 6

Bar plots of feasible rates of the different methods.

Close modal
Figure 7

Box plots of numbers of generations capturing the first feasible individual across 25 runs on (a) TGR-GR cases and (b) SBYR-GHYR cases.

Figure 7

Box plots of numbers of generations capturing the first feasible individual across 25 runs on (a) TGR-GR cases and (b) SBYR-GHYR cases.

Close modal

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.

Figure 8

Scheduled (a) TGR levels and (b) GR levels obtained with the ECHT.

Figure 8

Scheduled (a) TGR levels and (b) GR levels obtained with the ECHT.

Close modal

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

Scheduled (a) SBYR levels and (b) GHYR levels obtained with the ECHT.

Figure 9

Scheduled (a) SBYR levels and (b) GHYR levels obtained with the ECHT.

Close modal

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.

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 × 105 and 6.0 × 108 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.

Figure 10

Bar plots of feasible rates of the different methods (under a smaller NFE).

Figure 10

Bar plots of feasible rates of the different methods (under a smaller NFE).

Close modal
Figure 11

Optimized hydroelectricity of ECHT minus that of an individual CHT in terms of (a) maximum result, (b) mean result and (c) minimum result in 25 runs (under a smaller NFE). The better results of NRAM and SR on S2 in (c) are not shown. If a CHT does not return any feasible solution in all runs, there is nothing to be shown.

Figure 11

Optimized hydroelectricity of ECHT minus that of an individual CHT in terms of (a) maximum result, (b) mean result and (c) minimum result in 25 runs (under a smaller NFE). The better results of NRAM and SR on S2 in (c) are not shown. If a CHT does not return any feasible solution in all runs, there is nothing to be shown.

Close modal

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.

Figure 12

Scatter plots of objective function values and constraint violations for the designed test cases.

Figure 12

Scatter plots of objective function values and constraint violations for the designed test cases.

Close modal

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.

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.

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

Banos
R.
Manzano-Agugliaro
F.
Montoya
F. G.
Gil
C.
Alcayde
A.
Gomez
J.
2011
Optimization methods applied to renewable and sustainable energy: a review
.
Renewable & Sustainable Energy Reviews
15
(
4
),
1753
1766
.
Clarkin
T.
Raseman
W.
Kasprzyk
J.
Herman
J. D.
2018
Diagnostic assessment of preference constraints for simulation optimization in water resources
.
Journal of Water Resources Planning and Management-ASCE
144
(
8
),
04018036
.
Coello
C. A.
2002
Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: a survey of the state of the art
.
Computer Methods in Applied Mechanics and Engineering
191
(
11–12
),
1245
1287
.
de Paula Garcia
R.
de Lima
B. S. L. P.
de Castro Lemonge
A. C.
Jacob
B. P.
2017
A rank-based constraint handling technique for engineering design optimization problems solved by genetic algorithms
.
Computers & Structures
187
,
77
87
.
Haddad
O. B.
Moravej
M.
Loáiciga
H. A.
2015
Application of the water cycle algorithm to the optimal operation of reservoir systems
.
Journal of Irrigation and Drainage Engineering
141
(
5
),
04014064
.
Houck
C. R.
Joines
J.
Kay
M. G.
1998
A Genetic Algorithm for Function Optimization: A MATLAB Implementation
.
North Carolina State University
,
Raleigh, NC
,
USA
.
Hu
T. F.
Mao
J. Q.
Pan
S. Q.
Dai
L. Q.
Zhang
P. P.
Xu
D. D.
Dai
H. C.
2018a
Water level management of lakes connected to regulated rivers: an integrated modeling and analytical methodology
.
Journal of Hydrology
562
,
796
808
.
Hu
T. F.
Mao
J. Q.
Tian
M. M.
Dai
H. C.
Rong
G. W.
2018b
New constraint-handling technique for evolutionary optimization of reservoir operation
.
Journal of Water Resources Planning and Management-ASCE
144
(
3
),
04017097
.
Labadie
J. W.
2004
Optimal operation of multireservoir systems: state-of-the-art review
.
Journal of Water Resources Planning and Management-ASCE
130
(
2
),
93
111
.
Li
F. F.
Shoemaker
C. A.
Wei
J. H.
Fu
X. D.
2013
Estimating maximal annual energy given heterogeneous hydropower generating units with application to the Three Gorges system
.
Journal of Water Resources Planning and Management-ASCE
139
(
3
),
265
276
.
Mallipeddi
R.
Suganthan
P. N.
2010
Ensemble of constraint handling techniques
.
IEEE Transactions on Evolutionary Computation
14
(
4
),
561
579
.
Mallipeddi
R.
Suganthan
P. N.
Pan
Q. K.
Tasgetiren
M. F.
2011
Differential evolution algorithm with ensemble of parameters and mutation strategies
.
Applied Soft Computing
11
(
2
),
1679
1696
.
Mao
J. Q.
Zhang
P. P.
Dai
L. Q.
Dai
H. C.
Hu
T. F.
2016
Optimal operation of a multi-reservoir system for environmental water demand of a river-connected lake
.
Hydrology Research
47
(
S1
),
206
224
.
Mezura-Montes
E.
Coello
C. A.
2011
Constraint-handling in nature-inspired numerical optimization: past, present and future
.
Swarm and Evolutionary Computation
1
(
4
),
173
194
.
Naeini
M. R.
Yang
T.
Sadegh
M.
AghaKouchak
A.
Hsu
K.
Sorooshian
S.
Duan
Q.
Lei
X.
2018
Shuffled Complex-Self Adaptive Hybrid EvoLution (SC-SAHEL) optimization framework
.
Environmental Modelling & Software
104
,
215
235
.
Nicklow
J.
Reed
P.
Savic
D.
Dessalegne
T.
Harrell
L.
Chan-Hilton
A.
Karamouz
M.
Minsker
B.
Ostfeld
A.
Singh
A.
Zechman
E.
2010
State of the art for genetic algorithms and beyond in water resources planning and management
.
Journal of Water Resources Planning and Management
136
(
4
),
412
432
.
Poláková
R.
2017
L-SHADE with competing strategies applied to constrained optimization
. In:
2017 IEEE Congress on Evolutionary Computation (CEC)
,
San Sebastian
, pp.
1683
1689
.
Rodrigues
M.
de Lima
B. S. L. P.
Guimarães
S.
2016
Balanced ranking method for constrained optimization problems using evolutionary algorithms
.
Information Sciences
327
,
71
90
.
Runarsson
T. P.
Yao
X.
2000
Stochastic ranking for constrained evolutionary optimization
.
IEEE Transactions on Evolutionary Computation
4
(
3
),
284
294
.
Runarsson
T. P.
Yao
X.
2005
Search biases in constrained evolutionary optimization
.
IEEE Transactions on Systems, Man, and Cybernetics – Part C: Applications and Reviews
35
(
2
),
233
243
.
Schardong
A.
Simonovic
S. P.
2015
Coupled self-adaptive multiobjective differential evolution and network flow algorithm approach for optimal reservoir operation
.
Journal of Water Resources Planning and Management-ASCE
141
(
10
),
04015015
.
Tessema
B.
Yen
G. G.
2009
An adaptive penalty formulation for constrained evolutionary optimization
.
IEEE Transactions on Systems, Man, and Cybernetics – Part A: Systems and Humans
39
(
3
),
565
578
.
Trivedi
A.
Sanyal
K.
Verma
P.
Srinivasan
D.
2017
A unified differential evolution algorithm for constrained optimization problems
. In:
2017 IEEE Congress on Evolutionary Computation (CEC)
,
San Sebastian
, pp.
1231
1238
.
Yeh
W. W. G.
1985
Reservoir management and operations models: a state-of-the-art review
.
Water Resources Research
21
(
12
),
1797
1818
.
Zamuda
A.
2017
Adaptive constraint handling and success history differential evolution for CEC 2017 constrained real-parameter optimization
. In:
2017 IEEE Congress on Evolutionary Computation (CEC)
,
San Sebastian
, pp.
2443
2450
.
Zhou
Y. L.
Guo
S. L.
Liu
P.
Xu
C. Y.
Zhao
X. F.
2016
Derivation of water and power operating rules for multi-reservoirs
.
Hydrological Sciences Journal
61
(
2
),
359
370
.