## Abstract

Conventional reservoir operation emphasizes power generation (PG) with ignoring downstream river ecosystem and water supply benefits for sustainable development. Compared with the model defining water supply and ecological flow requirements as constraints, a novel long-term multi-objective scheduling model in complex parallel reservoir system (LTMOSCPRS) is developed to assess and achieve win–win and sustainable development for energy, water supply, and ecological benefits. The suitable and ideal ecological water requirements are calculated based on the requirement level index. Afterward, the novel multi-objective shuffled frog leaping algorithm (NMOSFLA) including renewed frog grouping, local search, and external elite frog set mutation strategies is proposed. Results indicate that three optimization objectives expose a mutual competing relationship. The benefit of the river ecosystem will increase at a loss of PG and water supply guarantee rate (WSGR). Therefore, the parallel reservoir system should be adjusted to improve the benefits of WSGR and ecological water spill and shortage (EWSS) with minimizing the loss of PG, simultaneously. Finally, the NMOSFLA is verified to outperform other compared methods at the solution diversity and convergence which is evaluated by multiple indicators. Overall, the NMOSFLA provides efficient reservoir operation schemes for decision-makers to select optimal trade-off schemes and feasible ways to solve the LTMOSCPRS.

## HIGHLIGHTS

Propose a novel multi-objective shuffled frog leaping algorithm (NMOSFLA) based on multiple improvement strategies.

Adopt the NMOSFLA to efficiently solve long-term multi-objective scheduling of complex reservoir systems.

Calculate the suitable and ideal ecological water requirements based on the requirement level index (RLI).

Verify that the NMOSFLA outperforms other methods on the solution diversity and convergence.

## INTRODUCTION

Reservoirs are one of the most efficient infrastructures for water resources management and planning (Zhang *et al.* 2016; Ming *et al.* 2017). Due to the increasing trend of water shortage, regulating reservoir-based water resources systems has become very important (Mansouri *et al.* 2017). Efficient reservoir operation strategies played a considerable role in flood control, electricity generation, water supply, and shipping benefits (Wu *et al.* 2012). The complexity of the reservoir system increases along with the increasing demand for water resources and energy. Conventional reservoir operation focuses on maximizing economic benefits without taking into serious attention the ecosystem health simultaneously (Labadie 2014). Because of the fast-growing water supply, energy demand, and ecological protection awareness, conventional reservoir operation regulations should be adjusted adaptively to coordinate economic and ecological benefits effectively.

Historically, multi-objective problems were decomposed into single-objective optimization problems and solved by methods including linear programming (LP) and nonlinear programming (NLP), dynamic programming (DP), and progress optimality algorithm (POA) (Feng *et al.* 2017; Ding *et al.* 2018; Zhang *et al.* 2019). However, DP may suffer from an exponential increase in computation and requires more computation spaces when solving complex high-dimension problems (Ming *et al.* 2015). Traditional multi-objective optimization techniques (weighting and E-constraint) are used to decompose the multi-objective problem and solved it by classical optimization algorithms (Guo *et al.* 2018). Li & Qiu (2016) gave corresponding weights to targets including power generation (PG), flood control, and ecology, and adopted a genetic algorithm to optimize the multi-objective operation of the Longyangxia reservoir. Long & Mei (2017) analyzed the scheduling scheme of the reservoir located in Jinsha River under different ecological flow constraints. These methods can simplify processes, especially when complex constraints are attached to the model. Nevertheless, several computations are required for obtaining a group of non-inferior solutions. Besides, the weight distribution between the optimization goals we strived for is inclined to be subjective to some extent, resulting in a large gap between the ideal solutions and the solutions we optimized. Due to multiple reservoir operation objectives, it is difficult to find a single solution that optimizes all goals simultaneously. Instead, multiple solutions exist in the form of trade-offs, also referred to as Pareto-optimal solutions.

To this end, multi-objective optimization measures such as multi-objective artificial bee colony (ABC), multi-objective differential evolution algorithm (MODE), and NSGA-II (Ashtiani *et al.* 2015; Aboutalebi *et al.* 2017; Pérez *et al.* 2017) have been increasingly applied to an optimization problem in water resources systems. Reddy & Kumar (2006, 2007) adopted a multi-objective genetic algorithm (MOGA) and multi-objective particle swarm optimization (MOPSO) to find optimal schemes for multi-objective reservoir systems including PG, agricultural irrigation, and downstream water quality requirements. Chang & Chang (2009) developed an integrated simulation model of two parallel reservoirs to formulate joint operating strategies based on independent operating policies for each reservoir of the system and used NSGA-II to evaluate the optimal joint operating strategies. Fang *et al.* (2018) proposed a multi-objective differential evolution-chaos shuffled frog leaping algorithm (MODE-CSFLA) for water resources system optimization to maximize water supply benefit. Liu *et al.* (2018) disclosed flood operation rules for reservoirs using a multi-objective cultured evolutionary algorithm based on decomposition (MOCEAD). Among previous studies, few studies involved reservoir optimal operation based on target objectives including PG, water supply, and ecological benefits simultaneously, especially parallel reservoir systems. Furthermore, the computation complexity increases when the number of reservoirs increases such as parallel and series connections in the multi-reservoir system (Azizipour *et al.* 2016). Therefore, it is necessary to develop a new effective method for solving a long-term multi-objective scheduling models in a complex parallel reservoir system (LTMOSCPRS).

To this end, the shuffled frog leaping algorithm (SFLA), as one of the biological intelligent population optimization algorithms, is applied to this topic which was primarily developed in 2003 (Eusuff & Lansey 2003). In comparison with other algorithms, the SFLA is a method of meme evolution based on a cooperative search strategy. The population in SFLA is divided into a few independent memeplexes or sub-populations. The local elite frog in each sub-population leads others to evolve from different directions independently. Afterward, frogs in respective sub-populations will shuffle and exchange information with each other. In recent years, the SFLA has been increasingly used for power system resource allocation, water resources system and optimization problems in the transportation system (Hidalgo *et al.* 2015; Orouji *et al.* 2016). Beykverdi & Ashouri (2014) proposed ICA-SFLA which incorporates an imperialistic competitive algorithm and the Nelder-Mead simplex algorithm. Sun *et al.* (2016) developed a new parallel normal cloud mutation SFLA and applied it to cascade reservoirs’ optimal operation. Meanwhile, researchers extended SFLA to multi-objective problems and proposed the multi-objective SFLA (MOSFLA) based on Pareto-optimality. He *et al.* (2014) proposed a multi-objective SFLA based on matrix binary encoding to solve cooperative task assignments. Zhang *et al.* (2014) developed the MOSFLA which incorporates the search mechanism of a quantum algorithm and applied it to the oilfield development plan.

The MOSFLA above has been verified to outperform another algorithm such as NSGA-II by researchers (Tian 2015; Fang *et al.* 2018). In this paper, to further enhance the performance of MOSFLA, a novel multi-objective SFLA (NMOSFLA) with an elite frog evolution mechanism, updated frog grouping and local search strategies are proposed. The NMOSFLA is used to solve the muli-objective operation model for the parallel reservoir system located in Yi and Luo rivers. Minimizing ecological water shortage and spill (EWSS), maximizing PG and water supply guarantee rate (WSGR) are defined as three target objectives. Suitable and ideal ecological water requirements are calculated based on the requirement level index (RLI). With help of the NMOSFLA, the trade-off relationship between target objectives is analyzed and performance indicators including the inverted generational distance (IGD), generational distance (GD), spacing (SP), and hypervolume (HV) are selected to demonstrate the effectiveness and feasibility of the proposed methods.

## THE OVERVIEW OF THE PARALLEL RESERVOIRS SYSTEM

Guxian reservoir is located on the Luo River and the catchment controlled by the reservoir is 5,370 km^{2}, accounting for 44.6% of the catchment controlled by Luo River. With a total designed storage capacity of 1.175 billion m^{3} and a beneficial reservoir capacity of 4.93 × 10^{8} m^{3}, the Guxian reservoir plays an important role in flood mitigation, power generation, irrigation, and water supply. Luhun reservoir, located in Songxian county, is also a large water conservancy project integrating the function of flood mitigation, power generation, irrigation, water supply, and fish breeding. The total storage capacity of the Luhun reservoir is 1.32 billion m^{3} and the flood mitigation capacity is 6.77 × 10 m^{3}.

*system 1*containing C1–C5 control sections.

*System 2*includes C6–C9 control sections which are in Yi River. The downstream section C10 is

*system 3*which contains the virtual reservoir and C10 control section.

Characteristic water level and parameters . | Reservoirs . | |
---|---|---|

Guxian . | Luhun . | |

Normal storage water level (m) | 534.8 | 319.5 |

Flood control level (m) | 529.3 | 317 |

Dead water level (m) | 495 | 298 |

Design flood water level (m) | 548.55 | 327.5 |

Total storage capacity (10^{8} m^{3}) | 11.75 | 13.2 |

Beneficial reservoir capacity (10^{8} m^{3}) | 4.93 | 5.83 |

Flood control storage | 6.01 | 6.77 |

Installed capacity (MW) | 60 | 10.65 |

Annual utilization hours of installed capacity (h) | 2,458.8 | 2,533 |

Comprehensive efficiency coefficient | 8.5 | 8.5 |

Characteristic water level and parameters . | Reservoirs . | |
---|---|---|

Guxian . | Luhun . | |

Normal storage water level (m) | 534.8 | 319.5 |

Flood control level (m) | 529.3 | 317 |

Dead water level (m) | 495 | 298 |

Design flood water level (m) | 548.55 | 327.5 |

Total storage capacity (10^{8} m^{3}) | 11.75 | 13.2 |

Beneficial reservoir capacity (10^{8} m^{3}) | 4.93 | 5.83 |

Flood control storage | 6.01 | 6.77 |

Installed capacity (MW) | 60 | 10.65 |

Annual utilization hours of installed capacity (h) | 2,458.8 | 2,533 |

Comprehensive efficiency coefficient | 8.5 | 8.5 |

## THE MODEL ESTABLISHMENT FOR THE LTMOSCPRS

To reflect the water supply benefit, the WSGR is proposed which can be elaborated into the ratio between water supply and corresponding demand. Moreover, multiple objectives also contain the maximum power generation for the hydropower station located in the reservoir and the minimum EWSS. Both ecological water spills and shortages can lead to adverse effects on ecological health and stability. The optimization objectives are described by Equations (1)–(3), respectively.

### Objective function

*C*is the quantity of control cross-sections; and

*T*denote total of reservoirs and scheduling time periods, respectively; is the hydropower output of the th hydropower station; is the efficiency coefficient of the unit in the

*i*th hydropower station; is the power generation flow for the th reservoir at th time period; is the reservoir water release corresponding to the th control cross-section; is the ecological water requirements downstream corresponding to the th control cross-section.

### Constraints

## THE NMOSFLA

### Review of the traditional single-objective SFLA

Inspired by the foraging behavior about frogs in the swamp, the SFLA is a metaheuristic optimization method based on particle swarm intelligence (PSO) as a local search tool and the idea of competitiveness and mixing information from parallel local from shuffled complex evolution (SCE) (Duan *et al.* 1992; Eusuff & Lansey 2003). The frog population is partitioned into several parallel sub-populations referred to as memeplexes, of which frogs update prey location and exchange information within each memeplex independently to some extent. This can be regarded as the local search strategy. Then, information communication between frogs in memeplexes is conducted through the shuffling process and global prey information is exchanged and preserved. The details of SFLA are elaborated as follows:

*Step 1.*Initializing frog population and calculating the fitness of each frog. Frogs are ranked in descending order according to their respective fitness.*Step 2.*The frog population (size ) is grouped into*m*sub-populations that contain*h*frogs. The frogs are divided into corresponding sub-populations in sequence until the*m*th frog is grouped in the*m*th sub-population. Then, the (*m*+ 1)th frog and the (*m*+ 2)th frog are arranged into the first and second sub-populations, respectively, once again,…, until all frogs are grouped into corresponding sub-populations (Sun*et al.*2016).*Step 3.*The worst frogs in each sub-population are updated according to the following equations which are used to update frog leaping step and prey locations, respectively.where is a random number in the interval [ 0,1]; and are the best and worst frog positions in each sub-population, respectively; is the renewed worst frog position. is leaping step of the th frog; and are lower and upper bounds of leaping step, respectively.*Step 4.*If the position of the new frog generated is inferior, the will be replaced by the global optimal frog to reproduce the new frog. A frog position is produced randomly if the new frog is still inferior than before. Then, the sub-populations are shuffled to exchange information carried by frogs from each sub-population and preserve the current global optimum (Luo*et al.*2015). Finally, executing local evolution and global shuffling processes repeatedly until iterations or convergence precision requirements are satisfied.

### External Elite Frog Set mechanism

#### The External Elite Frog Set maintenance

The External Elite Frog Set (EEFS) preserved the non-dominated solutions optimized in each iteration. These frogs can guide and modify the evolution direction of the frog population. If the new elite frog generated and the original frog dominates each other, both frogs are included in the EEFS for an alternative optimal solution. The details of the update mechanism are as follows:

- •
The alternative frog will be added directly to the elite frog set if the frog set is empty.

- •
The new alternative frog will be moved into the elite set if it dominates other frogs in the set. The frogs dominated by the new frog will be deleted.

- •It is averse to the search efficiency if the scale of EEFS gradually grows as iteration proceeds. To this end, we developed a crowding distance strategy (CDS) to control the elite frog quantity. The frog set scale is defined as . First, we calculated the crowding distance of frogs and rank them in descending order. Then, the redundant low-ranking frogs are eliminated from the set. Finally, the frogs with large crowding distance are retained to maintain diversity and distribution uniformity in the elite frog set. The crowding distance is calculated by the following equation.where denotes the crowding distance of the
*m*th elite frog. , and represent fitness of*m*− 1,*m*, + 1th frog corresponding to the th optimization objective, respectively. The and are the best and worst fitness in the elite frog set, respectively. Usually, all the elite frogs are selected into the elite frog set and ranked by their crowding distances to eliminate the frogs with lower crowding distance. That may make optimal solutions mainly distributed in the narrow space even though the number of judgments of dominance relation can be reduced. Finally, all frogs in this narrow space may be excluded from the EEFS. To this end, an improved maintenance strategy is proposed to add into elite frogs successively and eliminated frogs with intensive distribution density.

#### The EEFS with mutation strategy

At the initial iteration stage, the frog population can search the local area of feasible domain space. Thus, the optimal solutions acquired by the frog population are limited and the EEFS can retain the elite solution information efficiently. However, as the iteration progresses, the amount of non-dominated solutions increases sharply. The partial optimal solution information is lost simultaneously even though the CDS helps to control the quantity of optimal solutions. Moreover, the frog evolution is easy to trap into local optimum with a strategy focusing on locally (frog sub-population) and globally optimal frogs.

*et al.*2002). These selected frogs with the mutation are temporarily stored in the elite frog redundant set (EFRS) and the scale of is regulated by following Equation (15). The frogs with a mutation in EFRS can increase the population diversity of EEFS and avoid local optimum to a certain extent. The population disturbance decision condition is described in Equation (16).where denotes the selected frogs from sorted EEFS; denotes the scale of EFRS;

*H*denotes the dimension of optimization problem; and represent upper and lower bounds of frog in the th dimension, respectively; denotes the mutation coefficient which is correlated with mutation distribution index ; denotes the random number in the interval [0,1]. denotes the controlling coefficient and denotes the scale of EEFS.

*M*and are population scale and leaping step of the th frog, respectively.

The details of the mutation strategy are:

*Step 1*If the frog scale of sorted EEFS is less than , the total frogs will be copied to the EFRS. Otherwise, frogs are randomly selected from elite frogs to join EFRS.*Step 2*Mutation operation for frogs in the EFRS is conducted according to Equations (13) and (14). Confirm that the frog population disturbance decision condition is satisfied or not according to Equation (16). If the frog leap step*d*is lower than the leaping step limitation, the random frogs are selected to replace the frogs with the worst fitness and quality. Otherwise, turn to the next step.*Step 3*Frogs with mutation operation in EFRS are moved into the EEFS according to the maintenance strategy elaborated in*section 4.2.1*.*Step 4*The EFRS will be eliminated and the mutation operation is terminated.

### The renewed grouping strategy based on the clustering method

For the multi-objective optimization problem, the non-dominated frogs in the population represent the potential optimal solution and guide other individuals to approach the optimal region. The frog sub-population is usually small leading to a decline in local search efficiency if the non-dominated frogs overspread in the sub-population. Moreover, several sub-populations will search in a similar feasible region if the space of non-dominated frogs from different sub-populations exists in overlapping areas. Thus, computation resources are wasted, and the algorithm is inclined to be prematurely convergent to the local Pareto-optimal front.

To this end, we proposed a modification strategy that the non-dominated and dominated frogs in the population are grouped separately. Non-dominated frogs with similar geometric positions in the population are classified into the same sub-population by clustering method so that sub-populations contain the best individuals from different regions. Furthermore, in terms of the dominated frogs, proximity between non-dominated and dominated frogs is introduced to accomplish the grouping task.

#### The non-dominated frog grouping

For the frog population , *m* sub-populations are divided which contain *h* frogs. The non-dominated frogs are gathered in a set which are grouped with a clustering measure. In other words, non-dominated frogs are divided into non-overlapping regions according to respective geometric positions. For the grouping process, two situations should be considered:

- (1)
The number of non-dominated frog is represented by

*a*and*m*denotes the total sub-populations. If , each non-dominated frog represents a sub-population. Then, the size of the sub-population is adjusted to . - (2)
If , the clustering (Kanungo

*et al.*2002) is used to find an optimal grouping of . The quantity of categories in the clustering method is defined as the quantity of sub-populations*m*. First, select the total of*m*frogs from the set , and the corresponding vector is reckoned as the center of clustering . Euclidean distance is adopted to quantify the distance between the target vector of the frog individual and the center of clustering. Then, the following steps are necessary to be conducted.

*Step 1:* Select the nearest cluster center for each frog individual as its category information *.* Then, frogs with the common information are assigned to the same sub-population.

*Step 3:*The mean square deviation which can be obtained by Equation (18) is defined as the standard measure function.

*Step 4:* Repeat *Steps 1* − *3* until the terminal condition is satisfied. The set contains the optimal grouping strategy for the non-dominated frogs.

#### The dominated frog grouping

Dominated frogs are assigned to corresponding which are ranked in descending order. The larger indicates that is more inclined to be the better individual (closer to ). Then, frogs will be assigned to sub-populations within their size limitation based on principles of conventional grouping strategy (Step 3 in *Section 4.1*).

In terms of circumstance in *Section 4.3.1* that individuals without dividing exist when all sub-populations have reached maximum size, these frogs will be assigned into different sub-populations randomly. Moreover, to increase population diversity and avoid unequal distribution, one individual is assigned to join a sub-population randomly selected from other sub-populations.

#### Grouping strategy at the end of the iteration

The non-dominated frogs gradually increase at the end of the iteration. The population diversity decreases even though the non-dominated frogs occupy the whole population. If the frog population is still grouped by the clustering method, individuals in each sub-population are gathered in a narrow area of the target space. In other words, the algorithm may be trapped in local optima under this circumstance. To this end, the frog grouping strategy should be adjusted to reactive the evolution ability of the frog population. Based on the clustering method proposed in *Section 4.3.1*, the best clustering center is located and we calculated the degree of approximation between each frog and the clustering center . The frogs are rearranged in accordance with the descending order of and assigned to corresponding sub-populations.

### The modified local search strategy

The crowding distance and distance from the optimal Pareto front are introduced concerning the elite frog selection strategy (Yang *et al.* 2019b). The elite frog is obtained from the two-selection decision to conduct the evolution and improve the diversity of the frog population. Details are as follows.

*Step 1.*Define the fitness of as . The roulette wheel is adopted to choose the smaller frog with great probability for the elite individual. The fitness of is calculated by:

*Step 2.*To avoid the excessively dense distribution of frogs selected, we screened out the individuals with large crowding distance . The roulette wheel is used to find the targeted frog with a large probability. The is assessed to determine the dominance relationship with the original worst frog . The worst frogs will be eliminated based on three cases of elimination decision.

*Case 1:* The will replace the if the former dominates the latter.

*Case 2:*Due to a markedly decrease in population diversity, may be dominated by . Therefore, a more robust frog should be found to replace the original frog with reactive search performance. The new frog is obtained based on the following equation.where is a random number in the range of [0,1]; denotes adjustment coefficient; and are corresponding initial and terminal values; represents controlling factor ranging from 0.4 to 0.7.

*t*and are current iteration numbers for local search in sub-population and global shuffling.

*T*and

*L*represent the maximum iterations for local search and global shuffling, respectively. The larger helps to escape from the local area by searching in a larger space, while in the later iteration stage a smaller maintains rapid convergence of the NMOSFLA.

*Case 3:* If and are in the same dominant position, let replace or keep unchanged with equal probability.

### The flow chart of the NMOSFLA

## CASE STUDY

### Ecological flow acquisition

In this work, based on the RLI, the suitable and ideal ecological flows for the control sections located in Yi and Luo rivers are determined considering the annual runoff and ecological health extent. The RLI can reflect the effects of factors including vegetation coverage, biological diversity, and water quality on the ecological flow (Lyu *et al.* 2016).

*et al.*2019a). In terms of the flood season, it will raise to 15% of the average annual runoff to guarantee fish breeding and aquatic ecosystem. Based on the BEF, the requirement level indices of suitable ecological flow (SEF) and ideal ecological flow (IEF) and corresponding ecological flow requirements are calculated by the following equations. The results can be reached in Supplementary Tables S1 and S2 (Appendix A).where , denotes suitable, ideal ecological flow, respectively; is the average monthly runoff; , and denote requirement level indices of basic, suitable and ideal ecological flow, respectively; is correction coefficient of SEF; represents corresponding characteristic index; and is determined by entropy evaluation method, representing membership degrees of variation coefficient and characteristic index of SEF to RLI of SEF, respectively. and are the annual variation coefficient and characteristic index of IEF, respectively. and denote membership degrees of and to the RLI of IEF, respectively.

### The encoding method and constraint handling

#### The encoding method

*T*is the scheduling time period (month). The and denote reservoir stage at the th time period of the Luhun and Guxian reservoirs, respectively. The and represent water supply flows at the th period of the Luhun and Guxian reservoirs, respectively.

#### The constraint handling strategy

In this paper, the water level corridor and penalty function are adopted simultaneously to handle constraints. The reservoir water level and water supply flow generated are restricted within the feasible zone to satisfy the constraints mentioned. The reservoir discharge and hydropower output constraints can be converted to water level constraints using the water balance formula and water level corridor. The conversion steps are elaborated in Yang *et al.* (2019b).

The convergent speed can be accelerated when the search zone is narrowed into the water level corridor. Nevertheless, the frog locations generated in the corridor could still be unfeasible (Zhou *et al.* 2010). To this end, the penalty function is used to eliminate infeasible frogs and further improve frog overall quality.

### The main parameter settings for the NMOSFLA

The control parameters for the NMOSFLA are selected as follows with contrastive analysis. The number of local iterations in sub-population and global shuffling iterations are set to 20 and 200, respectively. The total iterations of chaotic mapping initialization are set to 100, the frog population size *M* is 200 and the corresponding sub-population or memeplex size *P* is predefined as 20. The scale of EEFS = 100. The interval of the frog leap step is set to [–4,4]. The controlling coefficient of the size of EFRS is defined as 0.50. Mutation distribution index = 5. The initial and terminal adjustment coefficient and are 0.90 and 0.40, respectively. The controlling factor = 0.60.

## RESULTS AND DISCUSSION

*P*= 50%). The results of the wet year (

*P*= 20%) and dry year (

*P*= 80%) are listed in Supplementary material, Appendix B (Figures S1 and S2).

It is obvious that the remarkable competitive relationship between the power generation, ecological benefit, and water supply is presented. It implies that the violation of ecological and water supply benefits will increase as the power generation is improved. Moreover, the less evenly distribute non-dominated solutions indicate that power generation, ecological benefit, and water supply are difficult to balance simultaneously irrespective of the reservoir inflow process. In terms of the EWSS and power generation, the solutions distribute more evenly than others which verifies the more remarkable competitive relationship between the two objectives. For the dry year, solutions are less uniformly distributed than that in the normal and wet years, implying that relative finite optimization space remains to take power generation and ecological benefit into full consideration due to the insufficient inflow. Furthermore, take the Guxian reservoir system, for instance, minimum PG and WSGR increase as the inflow becomes larger. The maximum PG is 1.786 × 10^{8}, 1.363 × 10^{8}, and 0.992 × 10^{8} kW h while the WSGR is 1.091, 1.081, and 1.014 for the wet, normal, and dry years, respectively. However, EWSS for downstream requirements has different trends as reservoir inflow increases. In detail, the ranges of EWSS corresponding to the normal and dry years are less than that in the wet year. It is smaller for the typical dry year with a minimum EWSS value of 5.820 × 10^{8} m^{3}, which indicates that the effect of power generation and water supply on the ecological benefit is relatively less in the dry year.

Furthermore, Figures 3 and 4 also demonstrate that mutual impact is remarkable among power generation, ecological and water supply benefits. For the Guxian reservoir system (the wet year), the power generation increases from 1.786 to 1.819 × 10^{8} kW · h when EWSS changes from 10.872 to 12.742 × 10^{8} m^{3} at a relatively steady rate. When the WSGR improves from 0.936 to 1.091, power generation decreases slowly at the range of 0.936 to 1 and then sharply declines from 1.813 to 1.786 × 10^{8} kW · h. At this stage, the WSGR improves at more cost of power generation benefit. A similar trend is presented between the EWSS and WSGR which reflects the mutual competitive relation. For the normal year, when WSGR improves from 0.907 to 1.081, corresponding power generation only decreases from 1.388 to 1.363 × 10^{8} kW · h with a relatively steady rate. Compared with wet years, the water supply benefit can be achieved with a higher guarantee rate with less power generation loss. A similar conclusion can be reached between EWSS and WSGR. Reservoir operation for a normal year has less impact on the power generation and ecological benefits as satisfying water supply benefit. For the dry year, due to the insufficient reservoir inflow, power generation is impaired considering water supply and ecological protection simultaneously. When the WSGR objective raises from 0.936 to 0.953, the influence on power generation is negligible, while a sharp decrease in power generation is incurred with the continuous increase of WSGR. Therefore, the operation regulation of the parallel reservoir system should be adjusted to improve the water supply and ecological benefits such as decreasing the loss of power generation.

It can be seen from Figures 5 and 6(b), the regulative space in the major flood season (July and August) is smaller than that in other months due to flood control task conducted by reservoirs. The differences between schemes (schemes *WSGR*, *PG,* and *EWSS*) focus on the dry season and water supply period after the flood season. As can be seen from Figures 5, 6(a) and 6(d), the reservoir discharge process of parallel reservoirs indicates that relatively abundant natural inflow in a typical wet year results in large water spills in May before major flood season. Still, the ecological water shortage in the dry season and water supply period at end of the flood season can be remarkably alleviated. For the normal year, the impact on ecology mainly reflects in form of water shortage at a later stage of flood season and dry season, especially January, February, and April. This can be attributed to the water supply task in these months at a loss of ecological benefit and partial power generation. In the future, with the increasing trend of water demand, the power generation and ecological benefit will be further sacrificed to guarantee the water supply benefit. Therefore, the operational requirements for reservoirs are inclined to be harsher. From May to September, the water shortage situation turns into excessive water spills due to an increase of water inflow and reservoir stage demand. For the dry year, the ecological water spill presents a downtrend, and water shortage in the dry season intensifies due to inadequate water inflow and an increase of water supply.

During the flood season of parallel reservoirs starting in July, flood mitigation is the premier task. According to the operation regulations, reservoirs begin to fill during the flood recession of September. The elevated water stage of the Guxian reservoir cannot exceed 534.8 m at the end of this month. To this end, the range of the water stage is from 524 to 534.8 m. For the wet and normal years, the water stage can reach to 534 around at the end of September and November without shortage of EWSS. While for the dry year, the water stage rises to 532 m around to satisfying water supply benefit and EWSS at the end of September at the loss of power generation.

More findings can be reached if we inspect Figures 5 and 6. Due to emphasis on EWSS benefit, the scheme *EWSS* for parallel reservoirs increase reservoir discharge and water stage drawdown occurs to alleviate ecological water shortage at the stage of the dry season. During the reservoir drawdown period, the water stage is not necessary to intensely decline, and corresponding water discharge decreases to avoid excessive ecological water spill. During the flood season, the risk of downstream ecology destruction is alleviated due to mitigating floods by the reservoir. These measures contribute to protecting the health of the downstream river ecological system. Nevertheless, the quick drawdown of the water stage results in the operation at a lower water stage for hydro-units, which obviously decreases the power generation benefit. The scheme *PG* focus on the power generation benefit. For these schemes, they cut down reservoir water discharge to maintain a higher water stage. Therefore, ecological water shortage inevitably amplifies because of power generation priority especially in the dry season (from January to April) and October. For the schemes *WSGR,* from January to April, the water supply task is relatively urgent due to less reservoir inflow compared with flood season. Thus, the ecological benefit is sacrificed to a certain extent. In terms of the dry year with less inflow, the regulated space is limited even though the *WSGR* is impaired in comparison with other years.

As discussed above, current operation regulations have difficulty balancing economic and ecological benefits to a certain extent. Hence, conducting the optimized scheduling of parallel reservoirs system is beneficial to protect ecological health and water supply benefit simultaneously. With the NMOSFLA, a set of non-dominated schemes can be generated, and decision-makers could select a compromise and acceptable scheduling scheme according to the requirements of the reservoir regulation to be implemented.

To further verify the advantages and efficiency of the NMOSFLA on multi-objective optimal operation for parallel reservoirs and test functions from ZDT and DTLZ packages, the IGD and GD indicators are used to evaluate scheduling solution PF (Li *et al.* 2017). The multi-objective genetic and evolution algorithms such as the fast-elitist non-dominated sorting genetic algorithm (NSGA-II), multi-objective evolutionary algorithm based on decomposition (MOEA/D), and multi-objective particle swarm optimization/based on decomposition (MOPSO/dMOPSO) are selected as a contrast. The indicator of IGD is to comprehensively evaluate the convergence and diversity of the algorithm. For the indicator of GD (Generational Distance), it is used to evaluate the distance between the real and searched Pareto fronts. The parameters including population scale and maximum iterations are equal. A total of the 20 simulations are conducted independently to obtain the final PFs and compute corresponding IGD, GD, HV, and SP. The results comparison (long-term parallel reservoirs scheduling) are listed in Tables 2 and 3.

. | . | The mean of IGD . | |||||
---|---|---|---|---|---|---|---|

Reservoir system . | Typical year . | NSGA-II . | MOEA/D . | dMOPSO . | MOEADDE . | MOPSO . | NMOSFLA . |

Guxian | The wet year | 2.35 × 10^{−02} | 8.65 × 10^{−02} | 4.60 × 10^{−03} | 4.84 × 10^{−02} | 2.24 × 10^{0} | 1.59 × 10 ^{−03} |

The normal year | 7.04 × 10^{−02} | 4.33 × 10^{−01} | 9.19 × 10 ^{−03} | 8.18 × 10^{−02} | 7.33 × 10^{0} | 1.96 × 10^{−02} | |

The dry year | 1.86 × 10^{−01} | 4.12 × 10^{−01} | 1.06 × 10^{−02} | 2.33 × 10^{−01} | 1.66 × 10^{1} | 1.00 × 10 ^{−02} | |

Luhun | The wet year | 6.90 × 10^{−03} | 8.42 × 10^{−02} | 6.57 × 10^{−03} | 2.71 × 10^{−02} | 4.85 × 10^{0} | 3.57 × 10 ^{−03} |

The normal year | 1.46 × 10^{−02} | 1.90 × 10 ^{−03} | 7.37 × 10^{−03} | 5.70 × 10^{−02} | 1.15 × 10^{1} | 3.04 × 10^{−03} | |

The dry year | 3.10 × 10^{−03} | 4.27 × 10^{−02} | 1.52 × 10 ^{−03} | 1.20 × 10^{−02} | 7.09 × 10^{0} | 2.59 × 10^{−03} | |

. | . | The STD of IGD . | |||||

Guxian | The wet year | 3.05 × 10^{−02} | 5.43 × 10^{−02} | 1.70 × 10^{−03} | 8.23 × 10^{−02} | 2.74 × 10^{0} | 9.03 × 10 ^{−04} |

The normal year | 6.10 × 10^{−02} | 1.09 × 10^{−01} | 3.40 × 10^{−03} | 2.16 × 10^{−01} | 6.25 × 10^{0} | 1.26 × 10 ^{−03} | |

The dry year | 1.47 × 10^{−01} | 2.00 × 10^{−01} | 4.12 × 10 ^{−03} | 3.53 × 10^{−01} | 1.61 × 10^{1} | 1.57 × 10^{−02} | |

Luhun | The wet year | 8.34 × 10^{−02} | 7.69 × 10^{−02} | 5.37 × 10^{−03} | 2.21 × 10^{−02} | 6.63 × 10^{0} | 3.59 × 10 ^{−03} |

The normal year | 2.22 × 10^{−01} | 1.73 × 10^{−01} | 6.17 × 10^{−03} | 5.94 × 10^{−02} | 1.97 × 10^{1} | 5.31 × 10 ^{−03} | |

The dry year | 5.93 × 10^{−01} | 3.89 × 10^{−01} | 1.60 × 10^{−02} | 1.60 × 10 ^{−02} | 8.22 × 10^{0} | 7.86 × 10^{−02} |

. | . | The mean of IGD . | |||||
---|---|---|---|---|---|---|---|

Reservoir system . | Typical year . | NSGA-II . | MOEA/D . | dMOPSO . | MOEADDE . | MOPSO . | NMOSFLA . |

Guxian | The wet year | 2.35 × 10^{−02} | 8.65 × 10^{−02} | 4.60 × 10^{−03} | 4.84 × 10^{−02} | 2.24 × 10^{0} | 1.59 × 10 ^{−03} |

The normal year | 7.04 × 10^{−02} | 4.33 × 10^{−01} | 9.19 × 10 ^{−03} | 8.18 × 10^{−02} | 7.33 × 10^{0} | 1.96 × 10^{−02} | |

The dry year | 1.86 × 10^{−01} | 4.12 × 10^{−01} | 1.06 × 10^{−02} | 2.33 × 10^{−01} | 1.66 × 10^{1} | 1.00 × 10 ^{−02} | |

Luhun | The wet year | 6.90 × 10^{−03} | 8.42 × 10^{−02} | 6.57 × 10^{−03} | 2.71 × 10^{−02} | 4.85 × 10^{0} | 3.57 × 10 ^{−03} |

The normal year | 1.46 × 10^{−02} | 1.90 × 10 ^{−03} | 7.37 × 10^{−03} | 5.70 × 10^{−02} | 1.15 × 10^{1} | 3.04 × 10^{−03} | |

The dry year | 3.10 × 10^{−03} | 4.27 × 10^{−02} | 1.52 × 10 ^{−03} | 1.20 × 10^{−02} | 7.09 × 10^{0} | 2.59 × 10^{−03} | |

. | . | The STD of IGD . | |||||

Guxian | The wet year | 3.05 × 10^{−02} | 5.43 × 10^{−02} | 1.70 × 10^{−03} | 8.23 × 10^{−02} | 2.74 × 10^{0} | 9.03 × 10 ^{−04} |

The normal year | 6.10 × 10^{−02} | 1.09 × 10^{−01} | 3.40 × 10^{−03} | 2.16 × 10^{−01} | 6.25 × 10^{0} | 1.26 × 10 ^{−03} | |

The dry year | 1.47 × 10^{−01} | 2.00 × 10^{−01} | 4.12 × 10 ^{−03} | 3.53 × 10^{−01} | 1.61 × 10^{1} | 1.57 × 10^{−02} | |

Luhun | The wet year | 8.34 × 10^{−02} | 7.69 × 10^{−02} | 5.37 × 10^{−03} | 2.21 × 10^{−02} | 6.63 × 10^{0} | 3.59 × 10 ^{−03} |

The normal year | 2.22 × 10^{−01} | 1.73 × 10^{−01} | 6.17 × 10^{−03} | 5.94 × 10^{−02} | 1.97 × 10^{1} | 5.31 × 10 ^{−03} | |

The dry year | 5.93 × 10^{−01} | 3.89 × 10^{−01} | 1.60 × 10^{−02} | 1.60 × 10 ^{−02} | 8.22 × 10^{0} | 7.86 × 10^{−02} |

. | . | The mean of GD . | |||||
---|---|---|---|---|---|---|---|

Reservoir system . | Typical year . | NSGA-II . | MOEA/D . | dMOPSO . | MOEADDE . | MOPSO . | NMOSFLA . |

Guxian | The wet year | 1.22 × 10^{−03} | 2.15 × 10^{−02} | 9.24 × 10^{−04} | 4.42 × 10^{−02} | 3.52 × 10^{−01} | 5.60 × 10 ^{−04} |

The normal year | 3.66 × 10^{−03} | 4.73 × 10^{−02} | 1.80 × 10 ^{−04} | 3.18 × 10^{−02} | 4.48 × 10^{−01} | 6.89 × 10^{−04} | |

The dry year | 6.00 × 10^{−03} | 8.13 × 10^{−03} | 2.07 × 10 ^{−03} | 1.68 × 10^{−01} | 8.28 × 10^{−01} | 3.53 × 10^{−03} | |

Luhun | The wet year | 5.03 × 10^{−04} | 1.50 × 10^{−03} | 1.41 × 10^{−05} | 4.79 × 10^{−03} | 4.22 × 10^{−01} | 4.79 × 10 ^{−06} |

The normal year | 1.11 × 10 ^{−05} | 3.26 × 10^{−03} | 3.63 × 10^{−05} | 9.94 × 10^{−03} | 1.11 × 10^{0} | 5.01 × 10^{−05} | |

The dry year | 2.46 × 10^{−04} | 7.09 × 10^{−03} | 9.31 × 10^{−05} | 2.07 × 10^{−02} | 2.90 × 10^{0} | 5.25 × 10 ^{−05} | |

. | . | The STD of GD . | |||||

Guxian | The wet year | 6.15 × 10^{−03} | 4.87 × 10^{−03} | 7.32 × 10^{−04} | 5.51 × 10^{−02} | 4.69 × 10^{−01} | 4.52 × 10 ^{−04} |

The normal year | 1.23 × 10^{−02} | 9.73 × 10^{−03} | 1.46 × 10 ^{−04} | 6.19 × 10^{−02} | 6.00 × 10^{−01} | 6.28 × 10^{−04} | |

The dry year | 1.74 × 10^{−02} | 1.79 × 10^{−03} | 1.77 × 10 ^{−04} | 1.01 × 10^{−01} | 1.22 × 10^{0} | 7.86 × 10^{−04} | |

Luhun | The wet year | 2.83 × 10^{−04} | 3.49 × 10^{−02} | 7.79 × 10^{−05} | 2.53 × 10^{−03} | 7.51 × 10^{−02} | 2.34 × 10 ^{−05} |

The normal year | 6.44 × 10^{−04} | 1.03 × 10^{−01} | 2.14 × 10^{−04} | 7.13 × 10^{−03} | 2.07 × 10^{−01} | 2.83 × 10 ^{−05} | |

The dry year | 1.46 × 10 ^{−04} | 3.06 × 10^{−01} | 5.88 × 10^{−04} | 2.01 × 10^{−02} | 5.73 × 10^{−01} | 3.43 × 10^{−04} |

. | . | The mean of GD . | |||||
---|---|---|---|---|---|---|---|

Reservoir system . | Typical year . | NSGA-II . | MOEA/D . | dMOPSO . | MOEADDE . | MOPSO . | NMOSFLA . |

Guxian | The wet year | 1.22 × 10^{−03} | 2.15 × 10^{−02} | 9.24 × 10^{−04} | 4.42 × 10^{−02} | 3.52 × 10^{−01} | 5.60 × 10 ^{−04} |

The normal year | 3.66 × 10^{−03} | 4.73 × 10^{−02} | 1.80 × 10 ^{−04} | 3.18 × 10^{−02} | 4.48 × 10^{−01} | 6.89 × 10^{−04} | |

The dry year | 6.00 × 10^{−03} | 8.13 × 10^{−03} | 2.07 × 10 ^{−03} | 1.68 × 10^{−01} | 8.28 × 10^{−01} | 3.53 × 10^{−03} | |

Luhun | The wet year | 5.03 × 10^{−04} | 1.50 × 10^{−03} | 1.41 × 10^{−05} | 4.79 × 10^{−03} | 4.22 × 10^{−01} | 4.79 × 10 ^{−06} |

The normal year | 1.11 × 10 ^{−05} | 3.26 × 10^{−03} | 3.63 × 10^{−05} | 9.94 × 10^{−03} | 1.11 × 10^{0} | 5.01 × 10^{−05} | |

The dry year | 2.46 × 10^{−04} | 7.09 × 10^{−03} | 9.31 × 10^{−05} | 2.07 × 10^{−02} | 2.90 × 10^{0} | 5.25 × 10 ^{−05} | |

. | . | The STD of GD . | |||||

Guxian | The wet year | 6.15 × 10^{−03} | 4.87 × 10^{−03} | 7.32 × 10^{−04} | 5.51 × 10^{−02} | 4.69 × 10^{−01} | 4.52 × 10 ^{−04} |

The normal year | 1.23 × 10^{−02} | 9.73 × 10^{−03} | 1.46 × 10 ^{−04} | 6.19 × 10^{−02} | 6.00 × 10^{−01} | 6.28 × 10^{−04} | |

The dry year | 1.74 × 10^{−02} | 1.79 × 10^{−03} | 1.77 × 10 ^{−04} | 1.01 × 10^{−01} | 1.22 × 10^{0} | 7.86 × 10^{−04} | |

Luhun | The wet year | 2.83 × 10^{−04} | 3.49 × 10^{−02} | 7.79 × 10^{−05} | 2.53 × 10^{−03} | 7.51 × 10^{−02} | 2.34 × 10 ^{−05} |

The normal year | 6.44 × 10^{−04} | 1.03 × 10^{−01} | 2.14 × 10^{−04} | 7.13 × 10^{−03} | 2.07 × 10^{−01} | 2.83 × 10 ^{−05} | |

The dry year | 1.46 × 10 ^{−04} | 3.06 × 10^{−01} | 5.88 × 10^{−04} | 2.01 × 10^{−02} | 5.73 × 10^{−01} | 3.43 × 10^{−04} |

Upon inspection of Tables 2 and 3, the IGD and GD indicators for the NMOSFLA is prominent than those obtained by compared methods. For the wet year of Guxian reservoir, the mean of IGD obtained by NMOSFLA and compared methods are 1.59 × 10^{−03}, 2.35 × 10^{−02}, 8.65 × 10^{−02}, 4.60 × 10^{−03}, 4.84 × 10^{−03} and 2.24 × 10^{0}. For the wet year of Luhun reservoir, mean of IGD are 3.57 × 10^{−03}, 6.90 × 10^{−03}, 8.42 × 10^{−02}, 6.57 × 10^{−03}, 2.71 × 10^{−02} and 4.85 × 10^{0}. Similarly, the mean of GD for Guxian reservoir (wet year) are 5.60 × 10^{−04}, 1.22 × 10^{−03}, 2.15 × 10^{−02}, 9.24 × 10^{−04}, 4.42 × 10^{−02} and 3.52 × 10^{−01}. Corresponding results for Luhun reservoir are 4.79E^{−06}, 5.03E^{−04}, 1.50E^{−03}, 1.41E^{−05}, 4.79E^{−03} and 4.22E^{−01}. In terms of STD of GD, for the Luhun reservoir, NMOSFLA outperforms others except for the dry year. Results discussed imply that NMOSFLA gets closer to the most preferred scheduling solutions for the parallel reservoirs system in comparison with other methods.

Comparison results (test functions from ZDT and DTLZ packages) are shown in Supplementary material, Appendix A: Tables S3 and S4. Similarly, results indicate that the NMOSFLA outperforms on most of the test functions in terms of IGD (ZDT1, 2; DTLZ1, 2, 3) and GD (ZDT1, 2, 3, 6; DTLZ1, 4) than that of other rivals. The results justify that the convergence and diversity of an approximation set obtained by the NMOSFLA are better than other algorithms to some extent.

Furthermore, we applied the modified algorithm (NMOSFLA) to solve the long-term multi-objective scheduling in series connection reservoirs (LTMOSSCR). The series connection reservoirs – Three Gorges-Gezhouba projects (TGGP) are used as a typical case to verify the performance of the NMOSFLA. The TGGP, located at the end of the upper reaches of the Yangtze River in China, consists of the Three Gorges hydropower project (TGP) and Gezhouba hydropower project (GP) and large-scale hydro-units. Here, we introduced the SP and hypervolume (HV) metrics to further demonstrate the advantage of the proposed NMOSFLA. The smaller SP and greater HV values usually represent that the algorithm has a better performance in diversity and convergence, respectively. The maximum, mean, minimum, and standard deviation (STD) values of SP and HV are shown in Table 4. The bolded values we marked represent the most excellent performance among all algorithms. It is known that both SP and HV require a set of Pareto-optimal front as a reference set. In this paper, the reference set was obtained by DP with the constraint handling method. We set the water supply flow and ecological flow as constraints to convert the multi-objective optimization problem to a single optimization problem. Then, the DP is repeatedly run to obtain a group of non-inferior solutions.

Typical year . | Algorithm . | SP . | HV . | ||||||
---|---|---|---|---|---|---|---|---|---|

Max . | Mean . | Min . | Std . | Max . | Mean . | Min . | Std . | ||

Wet year | NMOSFLA | 0.0194 | 0.0104 | 0.0081 | 0.0011 | 0.9061 | 0.9022 | 0.8975 | 0.0020 |

dMOPSO | 0.0721 | 0.0498 | 0.0191 | 0.0279 | 0.8608 | 0.7842 | 0.6030 | 0.1051 | |

MOEA/D | 0.0847 | 0.0594 | 0.0303 | 0.0173 | 0.9841 | 0.8234 | 0.6709 | 0.0812 | |

MOPSO | 0.0854 | 0.0214 | 0.0085 | 0.0119 | 0.9714 | 0.7959 | 0.5168 | 0.1335 | |

NSGA-II | 0.0167 | 0.0095 | 0.0076 | 0.0007 | 0.8976 | 0.8866 | 0.8790 | 0.0048 | |

Normal year | NMOSFLA | 0.0165 | 0.0132 | 0.0083 | 0.0014 | 0.8496 | 0.8343 | 0.8238 | 0.0089 |

dMOPSO | 0.1010 | 0.0395 | 0.0163 | 0.0205 | 0.8652 | 0.7068 | 0.5341 | 0.0720 | |

MOEA/D | 0.0655 | 0.0573 | 0.0305 | 0.0068 | 0.8607 | 0.8044 | 0.7414 | 0.0226 | |

MOPSO | 0.1274 | 0.0251 | 0.0096 | 0.0344 | 0.9320 | 0.7250 | 0.5251 | 0.1155 | |

NSGA-II | 0.0158 | 0.0142 | 0.0085 | 0.0019 | 0.8332 | 0.8217 | 0.8120 | 0.0074 | |

Dry year | NMOSFLA | 0.0146 | 0.0111 | 0.0098 | 0.0012 | 0.8408 | 0.8320 | 0.7443 | 0.0206 |

dMOPSO | 0.1289 | 0.0604 | 0.0139 | 0.0396 | 0.7021 | 0.6189 | 0.4078 | 0.0984 | |

MOEA/D | 0.2126 | 0.0902 | 0.0072 | 0.0641 | 0.8287 | 0.7255 | 0.7504 | 0.3030 | |

MOPSO | 0.0172 | 0.0119 | 0.0092 | 0.0021 | 0.8172 | 0.7332 | 0.4454 | 0.1233 | |

NSGA-II | 0.0877 | 0.0123 | 0.0277 | 0.0177 | 0.8733 | 0.7899 | 0.3018 | 0.1412 |

Typical year . | Algorithm . | SP . | HV . | ||||||
---|---|---|---|---|---|---|---|---|---|

Max . | Mean . | Min . | Std . | Max . | Mean . | Min . | Std . | ||

Wet year | NMOSFLA | 0.0194 | 0.0104 | 0.0081 | 0.0011 | 0.9061 | 0.9022 | 0.8975 | 0.0020 |

dMOPSO | 0.0721 | 0.0498 | 0.0191 | 0.0279 | 0.8608 | 0.7842 | 0.6030 | 0.1051 | |

MOEA/D | 0.0847 | 0.0594 | 0.0303 | 0.0173 | 0.9841 | 0.8234 | 0.6709 | 0.0812 | |

MOPSO | 0.0854 | 0.0214 | 0.0085 | 0.0119 | 0.9714 | 0.7959 | 0.5168 | 0.1335 | |

NSGA-II | 0.0167 | 0.0095 | 0.0076 | 0.0007 | 0.8976 | 0.8866 | 0.8790 | 0.0048 | |

Normal year | NMOSFLA | 0.0165 | 0.0132 | 0.0083 | 0.0014 | 0.8496 | 0.8343 | 0.8238 | 0.0089 |

dMOPSO | 0.1010 | 0.0395 | 0.0163 | 0.0205 | 0.8652 | 0.7068 | 0.5341 | 0.0720 | |

MOEA/D | 0.0655 | 0.0573 | 0.0305 | 0.0068 | 0.8607 | 0.8044 | 0.7414 | 0.0226 | |

MOPSO | 0.1274 | 0.0251 | 0.0096 | 0.0344 | 0.9320 | 0.7250 | 0.5251 | 0.1155 | |

NSGA-II | 0.0158 | 0.0142 | 0.0085 | 0.0019 | 0.8332 | 0.8217 | 0.8120 | 0.0074 | |

Dry year | NMOSFLA | 0.0146 | 0.0111 | 0.0098 | 0.0012 | 0.8408 | 0.8320 | 0.7443 | 0.0206 |

dMOPSO | 0.1289 | 0.0604 | 0.0139 | 0.0396 | 0.7021 | 0.6189 | 0.4078 | 0.0984 | |

MOEA/D | 0.2126 | 0.0902 | 0.0072 | 0.0641 | 0.8287 | 0.7255 | 0.7504 | 0.3030 | |

MOPSO | 0.0172 | 0.0119 | 0.0092 | 0.0021 | 0.8172 | 0.7332 | 0.4454 | 0.1233 | |

NSGA-II | 0.0877 | 0.0123 | 0.0277 | 0.0177 | 0.8733 | 0.7899 | 0.3018 | 0.1412 |

The results imply that the NMOSFLA outperforms other algorithms on the SP criterion in terms of the normal and dry years while dominating other rivals in terms of HV for 3 typical years. As discussed above, the NMOSFLA outperforms other methods on convergence and diversity of schemes, implying that it is an effective method for solving multi-objective scheduling problems both in series connection and parallel reservoirs.

## CONCLUSION

The conservation of river ecosystems and water supply benefits affected by reservoir operation has attracted increasing attention recently. To realize the harmonious development of power generation, water supply and ecological protection tasks, the model for the LTMOSCPRS is established. The model takes maximum power generation, water supply guarantee rate, and minimum volume of ecological water shortage and spill as three optimization objectives. Moreover, the NMOSFLA is applied to solve the problem. Results in the form of set of effective non-dominated schemes between each objective can be obtained. It was found that target objectives expose a competing relationship which means decision-makers need to select a compromised scheme according to operation regulations and requirements. Furthermore, a loss of power generation benefit would emerge if the water supply guarantee rate is maintained at a higher level and an adequate volume of downstream ecological water is provided by parallel reservoirs. Finally, in comparison with other algorithms, the NMOSFLA can find more feasible operation schemes based on modification strategies. The economic and ecological benefits are coordinated to reach their balance. The performance of the NMOSFLA is verified through test functions, two cases of series connection reservoirs, and parallel reservoirs. The NMOSFLA outperforms solution diversity and convergence compared with other selected algorithms.

## ACKNOWLEDGEMENTS

This work is supported by the National Natural Science Foundation of China (Grant No. 52109034); the Chinese Universities Scientific Fund (Grant Nos 2452021085 and 2452021086); Sichuan University, State Key Laboratory of Hydraulics and Mountain River Engineering (Grant No. SKHL2112); Key Laboratory of Resource Environment and Sustainable Development of Oasis, Gansu Province (Grant No. GORS202101), and the Research and Extension Project of Hydraulic Science and Technology in Shanxi Province (Grant No. 2017DSW02). The authors appreciated the insightful comments and suggestions from editors and anonymous reviewers.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.