Multiobjective Long-Term Generation Scheduling of Cascade Hydroelectricity System Using a Quantum-Behaved Particle Swarm Optimization Based on Decomposition

Multiobjective long-term generation scheduling (MOLTGS) plays a vital role in coordinating the contradiction between the generation and reliability of cascade hydroelectricity system (CHS). In this paper, a multiobjective quantum-behaved particle swarm optimization based on decomposition (MOQPSO/D) is presented to solve the MOLTGS problem of maximizing total hydroelectricity generation and firm hydroelectricity output. In MOQPSO/D, the improved logistic map is adopted to initialize the population in the feasible space, where various equality and inequality constraints are handled by a constraint handling method based on conversion and repair strategies. An improved Tchebycheff decomposition with a modified generator of direction vectors, a modified QPSO operator and normal cloud mutation are the key components of MOQPSO/D in addressing the nonlinearity and complexity of MOLTGS. The feasibility and effectiveness of MOQPSO/D are verified first on numerical experiments of benchmark instances and then on real engineering examples of the Three Gorges (TG) and Gezhouba (GZB) plants within three representative years. The results show that MOQPSO/D can not only perform better robustness, convergence and diversity performance than other seven competitors but also output a group of Pareto optimal schemes within a reasonable amount of time. So, this paper provides a new effective alternative to solve the MOLTGS problem of CHS when considering both generation and reliability objectives.


I. INTRODUCTION
With the ongoing cascade hydroelectricity development in China, hydroelectricity has claimed an ever-growing share of the nation's energy market over the past generation. By the end of 2015, approximately 80% of the technical exploitable hydroelectricity, accounting for 20.9% of the total installed capacity in China, had been used, with more than 1.0 trillion kW·h of hydroelectricity generation per year [1]. The development of appropriate long-term generation scheduling (LTGS) schemes for cascade hydroelectricity The associate editor coordinating the review of this manuscript and approving it for publication was Sotirios Goudos . systems (CHSs) is essential to improve efficiency [2]. Arguably, maximizing the generation benefit is the primary objective of LTGS for CHSs [3]. However, this objective tends to result in CHSs storing more water to obtain a high forebay level and high generation efficiency for as long as possible within the scheduling horizon. Because of the seasonal variations in inflow caused by the monsoonal climate in China, during some periods, water will be insufficient for hydroelectricity generation, leading to poor system reliability [4]- [7]. To balance the conflicting objectives of generation benefit and system reliability within the constrained search space, LTGS has become a multiobjective optimization problem (MOP), i.e., multiobjective LTGS (MOLTGS).
The nonlinearity of both constraints and objectives and the complexity of the CHS make the MOLTGS problem difficult to solve.
Since more than one objective should be considered, MOLTGS problems cannot be solved directly by single-objective optimizers, such as mathematical programming (linear, nonlinear and dynamic) and heuristic algorithms (HAs) [8], [9]. Thus, some conversion processes can be adopted to transform the MOP into a mono-objective problem with constraints, weights or penalties [10]. Wang et al. [11] proposed a constraint technique for transforming the vector optimization problem involving multiple stakeholders into a scalar problem with hydropower as the objective and the other two factors as constraints. In the study conducted in [12], a four-objective optimization problem was processed, where two objectives were treated as constraints and the other two were merged using weight coefficients. Yang et al. [13] integrated the generation and reliability objectives using a penalty function for a hydro/photovoltaic hybrid system. However, none of these methods are perfect, even though there have been many successful applications; for example, in some cases, the solution space shrinks, and as a result, some preferred solutions are ignored, the weights of objectives cannot be precisely obtained, and optimal solutions cannot approximate the nonconvex portions of the true Pareto front. Furthermore, single-objective optimizers also face difficulties because they yield only one solution within one simulation run, but one optimal solution cannot normally satisfy all conflicting objectives [14], [15]. Therefore, various multiobjective HAs (MOHAs), which endow HAs (like immune optimization, artificial fish swarm algorithm and kidney algorithm) with multiobjective searching capability, have been proposed by researchers in recent decades [16]- [21]. Qin et al. [22] integrated differential evolution into a cultural framework to develop a hybrid multiobjective optimizer and then adopted this method to address the safety issues of a dam and its upstream and downstream floods. MOHAs treat the multiple objectives of MOLTGS as a whole and adopt fitness assignment mechanisms (like alternating objective-based fitness assignment and domination-based fitness assignment) to associate each solution with a relative fitness value to measure its quality. However, these algorithms require considerable efforts to iterate solutions due to high computational cost of the fitness assignment, and diversity maintenance techniques are indispensable to overcome premature convergence and prevent the rapid collapse of population diversity during the search [23], [24].
Quantum-behaved particle swarm optimization (QPSO) is an enhanced swarm intelligence optimizer, in which particle swarm optimization (PSO) comes equipped with quantum mechanics; this technique was first proposed by Sun et al. [25]. The bound state of the particles in QPSO displays quantum behavior, making it impossible to capture the position and velocity simultaneously [25]. QPSO has better global search ability than PSO because the particles can reach any arbitrary position within the search scope via the guidance of a local attractor and a global point denoted as the mean best [26]. In addition, QPSO has few control parameters in its position update equation and can guarantee global convergence with a probability of 1, which makes it much simpler and more accurate than pure PSO. QPSO was designed for MOLTGS and has achieved great success [27], [28]. Moreover, decomposition-based MOHAs are widely developed since decomposition-based fitness assignment was proposed by Zhang and Li [24]. In such MOHAs, a uniform distribution of weight vectors aggregates the objectives of an MOP into a set of interoperable mono-objectives, and then these mono-objectives are optimized simultaneously by single-objective HA with neighborhood information [24]. Many MOPs have been solved by decomposition-based MOHAs because of their simplicity and effectiveness [23], [29], [30]. Conventional Tchebycheff decomposition and simulated binary crossover (SBX) are widely used decomposition and evolution approaches [31], which bring decomposition-based MOHAs two shortcomings: uniformly distributed solutions cannot always be obtained by conventional Tchebycheff decomposition [32], and the new solutions generated by SBX are often poor in quality [33], [34]. Under such circumstances, this paper proposes a new variation of multiobjective QPSO (MOQPSO), called MOQPSO based on decomposition (MOQPSO/D), in which uniformization of subproblem improvement regions and nondimensionalization of objective functions are addressed by an improved Tchebycheff decomposition with a modified generator of direction vectors [30]; the possibility of generating an excellent solution is greatly enhanced by a modified QPSO operator, and the ability to avoid local optima is improved with the use of normal cloud mutation [35]. This paper is the first to integrate the above techniques to address MOPs. Finally, with the initial population generation and constraint handling method, the applicability of MOQPSO/D to 21 benchmark instances and 3 MOLTGS problems of the Three Gorges (TG) and Gezhouba (GZB) plants is presented to demonstrate the feasibility and effectiveness.
The remainder of the paper is structured as follows. The mathematical model of the MOLTGS problem and the implementation of MOQPSO/D are discussed in section II and III, respectively. Then, the applications of MOQPSO/D to numerical experiments and engineering practices are presented in section IV and V, respectively, followed by the conclusion of this paper.

II. MATHEMATICAL MODEL A. OBJECTIVE FUNCTIONS
Total hydroelectricity generation and firm hydroelectricity output of a CHS are defined, respectively, as the total amount of system hydroelectricity generation and the minimum system hydroelectricity output in all scheduling intervals, which are selected as objective functions. Maximizing both objectives simultaneously would strike a balance between generation benefit and system reliability, being of great significance for generating as much and stable hydroelectricity as possible with the available water resources [5].
where f 1 and f 2 refer to the total hydroelectricity generation and firm hydroelectricity output of the CHS for the entire scheduling horizon, respectively; t and m are the interval and plant indices, respectively; T and M are the total numbers of intervals and plants, respectively; C m is the conversion coefficient of water and electricity in plant m; H m,t denotes the difference between the average forebay and tail levels of plant m in interval t; P m,t is the hydroelectricity output of plant m in interval t; Q P m,t is the water discharge rate through the powerhouse of plant m in interval t, which is the difference between the total water discharge rate Q O m,t and spillage water discharge rate Q S m,t ; and t denotes the length of interval t.

B. CONSTRAINTS
The feasible ranges of variables are delimited by the following rigid constraints, and between these constraints, tight coupling makes it possible to cast each feasible range into the decision space to find a feasible decision space without violating any of the constraints [30]. The variables mentioned here are average values in all scheduling intervals except forebay level and storage.

1) WATER BALANCE EQUATION
where V m,t denotes the initial storage volume of plant m in interval t, while V m,t+1 denotes that in the next interval; and Q I m,t is the water inflow rate of plant m in interval t.

2) HYDRAULIC CONNECTION
where Q L m+1,t is the local water inflow rate of downstream plant m + 1 in interval t.   where Z ini m and Z ter m are the initial and terminal boundary conditions of forebay level Z m,t (t = 1, 2, . . . , T + 1), respectively.

III. IMPLEMENTATION OF MOQPSO/D FOR MOLTGS A. QUANTUM-BEHAVED PARTICLE SWARM OPTIMIZATION
The position and velocity of the particles in QPSO cannot be determined simultaneously according to the uncertainty principle [25]. Hence, without using position and velocity, QPSO adopts a wave function ψ (Schrödinger equation) to depict the particle state and the Monte Carlo method to determine the particle positions. The position update mechanism is as follows: where x k n represents the nth particle at iteration k, while x k+1 n represents that at the next iteration; α k is the contractionexpansion (CE) coefficient for tuning the velocity of convergence [26]; r 1 and r 2 are two numbers randomly generated in range [0, 1]; K is the maximum number of iterations; and mb k and la k n represent the mean best of the population and the local attractor of particle n at iteration k, respectively, which can be calculated by taking the global best gb k n = gb k n,1 , gb k n,2 , . . . , gb k n,D and personal best pb k n = pb k n,1 , pb k n,2 , . . . , pb k n,D as follows: la k n = la k n,1 , la k n,2 , . . . , la k where D is the dimension of the particles, and ϕ 1 and ϕ 2 are two numbers randomly generated in range [0, 1]. Since the acceleration coefficients c 1 and c 2 are always set to the same value in QPSO [36], ϕ becomes a random number in (0, 1).
Pure QPSO is designed to handle mono-objective optimization problems, and its performance has been investigated in many studies [36]- [38]. This operator was developed to solve MOPs in MOQPSO/D, in which the definitions of the personal best and global best were modified for use with an improved Tchebycheff decomposition. Moreover, the nondominated solutions found by MOQPSO/D are archived into an external set, filtered by the crowding distance [39] and updated by a jumping-out mechanism [40]. Additional details are given in section III-E.

B. INITIAL POPULATION GENERATION
The storages at the end of each scheduling interval are used as decision variables when solving the MOLTGS problem of a CHS. For T intervals over the entire scheduling horizon with M hydroelectricity stations, there are D = M (T +1) schedules for storages. Thus, a population of N particles, with each particle representing a solution to the MOLTGS problem, is expressed as follows, Eq. (12), as shown at the bottom of this page, where k is the current iteration number, which is set to 0 during initialization.
For the merits of ergodicity, symmetry and randomicity, chaotic sequences produced by improved logistic map can endow the initial population with desirable quality and diversity [41], [42]. The formulae for rescaling the chaotic sequences into the decision space are stated as: where u n = [u n,1 , u n,2 , . . . , u n,D ] denotes the chaotic sequences with u n,d ∈ (−1, 1) and u n,d / ∈ {-0.5, 0, 0.5}, and is the search range of the dth decision variable.

C. CONSTRAINTS HANDLING METHOD
In the MOLTGS problem, two kinds of constraints should be addressed simultaneously: inequality constraints and equality constraints. Conventional constraint handling method often uses a repair strategy to handle the storage constraint when initializing the population or treating newly generated solutions; this approach forces the equality constraints to be satisfied when calculating total water releases and hydroelectricity outputs and adopts a penalty function to address the total water release and hydroelectricity output constraints by reducing the fitness of infeasible solutions [43]. However, although infeasible solutions will have a greater chance of being eliminated than feasible solutions, some search ability has been wasted in the infeasible portions of the search space.
Since the inequality and equality constraints of the MOLTGS problem are highly coupled, the following constraint handling method can convert the total water outflow and hydroelectricity output constraints into storage constraints [44]. Take the mth hydroelectricity plant as an example and denote v m () as the storage volume corresponding to forebay level The analysis occurs in a positive sequence under the assumption that plant m runs on the lower and upper total water release boundaries from the first scheduling interval to the t-1th scheduling interval; then, the positive storage range [V Then, the storage limits are updated with the following formulae: The repair strategy will consider the updated storage scope for each two adjacent scheduling intervals. Storage generated within the updated limits may also constitute a violation of the total water release constraint. Therefore, taking the lower and upper total water release boundaries for Eq. (3), Eq. (17) is adopted to determine the repaired storage by leaving the storages before and after unchanged. A similar procedure can also be adopted to convert the hydroelectricity output H. Hu, K. Yang: MOLTGS of CHS Using a QPSO Based on Decomposition constraint. (17) where V k m,t is plant m's storage volume in interval t at iteration k.

D. IMPROVED TCHEBYCHEFF DECOMPOSITION
While conventional Tchebycheff decomposition has been widely applied to design decomposition-based MOHAs with a uniform distribution of weight vectors, there is no guarantee that the obtained solutions are uniformly distributed, especially for some complex MOPs such as MOLTGS [32]. Moreover, the objective functions of the MOLTGS problem are not nondimensionalized in conven-  tional Tchebycheff decomposition. Therefore, an improved Tchebycheff decomposition is presented here to aggregate the objective functions in Eqs. (1) and (2) into a group of interoperable mono-objective functions, which can be stated as follows: where λ n = λ n,1 , λ n,2 , . . . , λ n,L with λ n 2 = 1 and λ n,l > 0(l = 1, 2, . . . , L) refers to a direction vector, ε = 0.00001 is added to prevent dividing by zero, and L = 2 is the number of objective functions. For the MOLTGS problem formulated in section II, z W l and z B l , respectively, denote the minimum and maximum function values of the corresponding objectives.
The improved Tchebycheff decomposition in Eq. (18) assigns uniform improvement regions for subproblems, leading to desirable diversity in the search results [30], [32]. However, the original generator of direction vectors is not suitable for the MOLTGS problem because the true Pareto front cannot be obtained in advance. A modified generator of direction vectors is proposed based on a quarter of a unit circle. As shown in Eq. (19), a uniform distribution of the direction vector can be generated without knowing the true

E. MODIFIED QUANTUM-BEHAVED PARTICLE SWARM OPTIMIZATION OPERATOR
According to Eqs. (9), (10) and (11), the flight of particles is steered by the personal best and global best cooperatively. MOQPSO/D defines the personal best as the previous best position to the corresponding subproblem for each individual in the swarm with regard to the aggregation function of the improved Tchebycheff decomposition. The previous best is initialized at the same position as the particle (i.e., pb 0 n = x 0 n ) since there is no previous movement in the initial population, and updating occurs by comparing the current previous best with the new offspring generated by the particle [30]. The global best of particles should be defined carefully in MOQPSO/D since different particles are associated with different direction vectors and different subproblems. One popular way, which is also used here, is to define the global best of a given particle as the position with the best aggregation value to the corresponding subproblem in the swarm [30]. This means is coupled with an external archive in which the nondominated solutions visited by all particles in previous iterations are collected; that is, the information  from the external archive serves to define the global best. This external archive first collects all the nondominated particles during initialization, and it is updated by new nondominated offspring in the following generation. During evolution, even though dominated particles will be removed from the archive, the number of archived solutions will rapidly increase with an increasing number of iterations. Therefore, a finite maximum size A is applied to filter redundant nondominated solutions with small crowding distances when the archive has no space for new nondominated offspring [30]. In addition, the elitist learning strategy, a jumping-out mechanism, is also applied to update the external archive; the pseudocode is provided in [40].

F. NORMAL CLOUD MUTATION
Under the guidance of the mean best and local attractor, the swarm will quickly converge to the current best position and rapidly lose diversity [35]. If the current best position is a local Pareto-optimal front (POF), then all particles will easily get stuck in it if the algorithm performs only a QPSO-based search [28]. Therefore, a perturbation search operator, i.e., normal cloud mutation, was designed for MOQPSO/D to address population diversity issues. This mutation causes a substitute for each selected variable to escape from local fronts. The number of particles affected by normal cloud mutation is controlled by a mutation rate mr [45], which decreases nonlinearly in the course of the process. For each solution x k+1 n , if flip((1 − (k + 1) K ) 5 / mr ) is valid, then the randomly selected variable x k+1 n,d mutates with the following: where r represents a random factor from 0 to 1. Cloud drop (x C , y C ) is calculated by the eigenvector (Ex, En, He).

IV. NUMERICAL EXPERIMENTS A. EXPERIMENTAL SETTINGS
In this section, seven popular MOHAs, i.e., AGE-MOEA [46], A-NSGA-III [47], CMOPSO [48], ENS-MOEA/D [49], KnEA [50], SMPSO [51] and MOQPSO [27], along with MOQPSO/D are applied to 21 dual-objective instances. These benchmarking instances are selected from commonly employed ZDT [52], DTLZ [53] and WFG [54] test suites, offering a comprehensive set of characteristics that are likely to make the performance of the compared algorithms fully tested. The parameter setting of each instance is identical VOLUME 8, 2020  to that in [48], including the population size, the maximum number of iterations and the number of decision variables. A total of 30 independent runs are executed for each algorithm on each test problem, and in each independent run Hypervolume [52] and Spacing [55] are used to measure the quality of the visited POFs from both convergence and diversity perspectives. Generally, one Hypervolume indicates perfect performance in both convergence and diversity, while zero Spacing indicates ideal diversity performance. Note that the reference point for Hypervolume is the worst function values of the corresponding objectives in the visited POFs.
All eight algorithms are coded in MATLAB to obtain fair results. It should be noted that the source codes of AGE-MOEA, A-NSGA-III, CMOPSO, ENS-MOEA/D, KnEA and SMPSO are downloaded from the website of PlatEMO [56], and the source code of MOQPSO is developed by the authors according to [27]. Unless otherwise specified, the parameters for competitors are set to values given in the corresponding references. For ENS-MOEA/D cases, the neighborhood sizes are {5, 10, 15, 20} for all instances.
For MOQPSO/D cases, a time-varying CE coefficient with a value that decreases linearly from 1.0 to 0.5 with the number of iterations is adopted for the modified QPSO operator [26]; a rate of 1/D, where D refers to the number of decision variables, is used to limit normal cloud mutation; and the maximum external archive size is equal to the population size. Since there is no constraint conversion involved in benchmark instances, the repair strategy is only required when the decision variable exceeds the predefined scope. The above algorithms were executed on a 3.60-GHz personal PC.

B. PERFORMANCE COMPARISONS
Tables 1 and 2 provide the statistical characteristics (mean and standard deviation) of Hypervolume and Spacing values, respectively, for each method. Furthermore, the Friedman test [57] is carried out to rank the eight contenders respectively according to their Hypervolume and Spacing performance. The ranks are also enlisted in Tables 1 and 2. Table 1 shows that MOQPSO/D statistically performs better in improving Hypervolume than all the rest on 8 out VOLUME 8, 2020 while MOQPSO/D is inferior to more than one peers on the four remaining problems (ZDT3, WFG4, WFG6 and WFG8) but still highly competitive. In Table 2 The above-mentioned Friedman ranks are acquired based solely on the mean values. Therefore, the stability of all candidate algorithms should be checked according to their standard deviation values in Tables 1 and 2. Only in rare cases may the stability of MOQPSO/D be worrying in Table 1, including ZDT 4, WFG 1 and WFG2. In Table 2, the alarming cases are ZDT4, ZDT6, DTLZ1, DTLZ3 and WFG1. Such stability performance is not the best among the candidate algorithms but still endows MOQPSO/D with core competitiveness. The superiority of MOQPSO/D drawn from the accuracy analysis is not swayed by the stability comparison.

V. CASE STUDY A. ENGINEERING BACKGROUND
The MOLTGS problem of the TG and GZB, two adjacent hydroelectricity infrastructures at the segmentation point of the Yangtze's upper and middle reaches (Fig. 4), China, serves to ascertain whether MOQPSO/D is feasible and effective. The former is the world's largest plant, with 34 hydroelectric turbines installed, and the latter is the world's largest run-of-river plant, with 22 hydroelectric turbines installed, both of which are operated for multiple water use sectors. Maximum total hydroelectricity generation and firm hydroelectricity output of the TG and GZB are selected as objectives, and the requirements of other sectors are transformed into constraints. In general, the TG's forebay level must not exceed the flood-limited forebay level of 145 m (from 144.9 m to 146.5 m, in practice) from June to September and 175 m in the other months for controlling the latent flood peak, and it must drop from no lower than 155 m to the flood-limited forebay level in May and rise from the flood-limited forebay level to 175 m in October. The initial and terminal boundary conditions of the TG's forebay level are both 175 m. The GZB is a reverse regulation project of the TG and is located 38 kilometers downstream from it. The main functions of the GZB are hydropower generation and navigation improvement. The forebay level is assumed to be 66 m in all months because of its small regulation ability. Currently, balancing generation benefit and system reliability is an urgent problem in the TG and GZB since water inflows are unevenly distributed.
The daily inflows of the TG are aggregated and averaged to form monthly time series within three representative years, whose annual inflows correspond to the 25, 50 and 75 percentiles of those from 1950 to 2018, representing high, moderate and low flow scenarios, respectively. The search scopes are [4990, 22500] MW and [1040, 2735] MW for the hydroelectricity outputs of the TG and GZB, respectively. The lower total water release bound is set to 5000 m 3 /s for both plants to satisfy the requirements of navigation and ecological protection, and the upper total water release bound is calculated by considering both the flood control requirement and the spillway and powerhouse outlet capacities. The other characteristics of the TG and GZB were given in [58].

B. PARAMETER SETTING
An investigation of optimal parameter setting is of great significance to unleashing the algorithm's full potential on a given problem, especially when solving a real engineering problem such as MOLTGS. A two-dimensional grid search [59] is employed here to screen out the optimal population size and maximum iteration number from sets {100, 300, 500, 1000} and {300, 500, 1000, 2000}, respectively, under the high flow scenario. For fair comparisons, A = 100 is fixed for the external archive during the search, and the other parameters are determined in the same way as the numerical experiments adopt. Table 3 lists the search results where the performance of each parameter combination is measured in terms of the mean total hydroelectricity generation, firm hydroelectricity output, Hypervolume, Spacing and computational time values. Each value in Table 3 is calculated by running MOQPSO/D with the corresponding parameters 30 times independently. It can be concluded from Table 3 that both a greater maximum iteration and a larger population are likely to make MOQPSO/D produce a better POF but introduce a more expensive computational cost, and the performance gain of large population is smaller than that of great iteration. The potential of large population may be restricted by the maximum external archive size. However, the selection of A = 100 is not arbitrary: MOQPSO/D can offer a suitable number of Pareto-optimal schemes (POSs) VOLUME 8, 2020  Fig. 5. This figure demonstrates that two optimization objectives collide with each other since an increase in total hydroelectricity generation will inevitably decrease firm hydroelectricity output in these fronts, and the ideal solution is shown in the top right corner of each subfigure. Visual comparisons suggest that the POFs of MOQPSO/D are widely spread and evenly distributed in the objective space and close to the top right corner, indicating that the present algorithm outperforms the other seven peers.
Moreover, to assess the performance in detail, the mean and standard deviation of the total hydroelectricity generation, firm hydroelectricity output, Hypervolume and Spacing are employed to measure the robustness, convergence and diversity of the final POFs achieved by the corresponding algorithms for the three flow scenarios. These statistical indicators are reported in Tables 4-7, and in the same time, the corresponding Friedman ranks are also reported in each table using the mean values. The data given in bold highlight the best results for each representative regime, i.e., the greatest total hydroelectricity generation, firm hydroelectricity output and Hypervolume and smallest Spacing values. ENS-MOEA/D is first excluded from being applied to the engineering problems because of its clear disadvantages in convergence and diversity. The mean objective values in Tables 6 and 7 also indicate that ENS-MOEA/D tends to offer one or a few schemes in each year. In the same manner, AGE-MOEA and KnEA are then knocked out. Even though they rank batter than MOQPSO/D in terms of Hypervolume, their poor performance in two objectives means that AGE-MOEA and KnEA always cannot guarantee to provide a POF as wide spread as MOQPSO/D achieves, which is of great meaning to broadening the decision makers' options. Finally, when competing with A-NSGA-III, CMOPSO, SMPSO and MOQPSO, the performance of MOQPSO/D is preferable, without question, because MOQPSO/D owns 2.8, 1.3, 2.3 and 1.0 average ranks, respectively in Tables 4, 5, 6 and 7, claiming the first places among the remaining algorithms. Meanwhile, the small standard deviations of MOQPSO/D in these tables, especially in Table 7, reconfirm its superiority from the view of stability. The evolution curves of the average Hypervolume and Spacing during the entire iteration process are presented in Fig. 6. It can be seen that the MOQPSO/D curves (red) converge to stable values with promising speed at about 400 iterations and plateaus at the third, third and second best levels at later iterations in Fig. 6(a), (c) and (e), respectively. As for Spacing, the advantages of the red lines are obvious. Additionally, the fluctuations in the red curves demonstrate that MOQPSO/D can avoid local optima by using normal cloud mutation. Therefore, the above analysis suggests that the developed alternative surpasses the other seven techniques for the MOLTGS problems of the TG and GZB within three representative years with respect to robustness, convergence and diversity metrics.
Finally, to assess the computational efficiency of MOQPSO/D, the average execution times of the corresponding approaches over 30 independent runs are recorded in VOLUME 8, 2020       method. However, ascertaining their contributions in this paper is unnecessary because they have been proved in many previous articles [4], [7].

D. RESULTS OF GENERATION SCHEDULING
Three typical schemes (corresponding to the objective function maximum for total hydroelectricity generation, the best compromise considering both objective functions and the objective function maximum for firm hydroelectricity output) reached by MOQPSO/D are also shown in Fig. 5 (with cyan markers of different shapes) for each representative flow scenario, and they reflect the generation-optimal scheme (GOS), compromise-optimal scheme (COS) and reliability-optimal scheme (ROS). The forebay level, hydroelectricity output and total water release trajectories of the selected typical schemes are shown in Fig. 8 for the representative flow scenarios. Note that the forebay level of the GZB and total water release of the TG are not presented here because of the negligible regulation ability of the GZB. Obviously, these forebay levels, hydroelectricity outputs and total water releases vary within the boundary constraints during the scheduling horizon, indicating the effectiveness of the initial population generation and constraint handling method. Their changes are similar within each representative year, with obvious differences from January to May. Thus, there is a competitive relation between total hydroelectricity generation and firm hydroelectricity output in this period, and it is necessary to choose a suitable scheme for balancing these two objectives. Clearly, the GOS releases less water from January to April to maintain a high forebay level, leading to low hydroelectricity output in these months and high hydroelectricity output in May. This operation process is generally unstable but can yield maximum total hydroelectricity generation. In contrast, the ROS discharges more water to increase hydroelectricity output in the first four months, thereby maintaining system stability but reducing total hydroelectricity generation by lowering the forebay level. As a compromise, the operation curves of the COS lie at nearly an equal distance from those of the GOS and ROS. Fig. 8 also shows that the range of options for the forebay level, hydroelectricity output and total water release decreases with decreasing inflow. The quantitative VOLUME 8, 2020 relationships between the two objectives are considered in Table 8, where the bracketed values are relative increases (↑) and reductions (↓) compared with the GOS for the corresponding objectives. A 1.00% decrease in total hydroelectricity generation can result in increases of at least 16.30%, 16.01% and 11.19% in firm hydroelectricity output in high, moderate and low flow years, respectively.

VI. CONCLUSION
MOLTGS plays a prominent role in guiding the operation of CHSs. Since more than one objective should be simultaneously optimized without violating constraints, MOLTGS has become a very complex MOP. Maximum total hydroelectricity generation and firm hydroelectricity output are selected as objective functions to develop a mathematical model for solving the MOLTGS problem, and a novel MOQPSO, namely, MOQPSO/D, is proposed to solve the mathematical model for the first time. During evolution, an improved Tchebycheff decomposition with a modified generator of direction vectors is introduced to divide the MOLTGS problem into a group of scalar subproblems, assigning uniform improvement regions to each subproblem and nondimensionalizing the objective functions. All subproblems are optimized simultaneously with a modified QPSO operator and normal cloud mutation, where the personal best and global best of each particle are defined with regard to the corresponding aggregation function. After initializing the population with improved logistic map and handling constraints with conversion and repair strategies, 21 dual-objective benchmarks and 3 MOLTGS problems of the TG and GZB are used to ascertain the MOQPSO/D's feasibility and effectiveness from the perspectives of both numerical experiment and engineering practice. The obtained results show that the developed MOQPSO/D method can obtain POSs with better performance than the other seven approaches in terms of robustness, convergence, diversity and computational efficiency. Therefore, a new superior alternative is provided for solving the MOLTGS problem by considering both generation benefit and system reliability. Future analyses should focus on the applicability of MOQPSO/D to more complex hydroelectricity systems. Meanwhile, the MOLTGS problems with more objectives or on finer time scales are also of interest to MOQPSO/D.