Optimization of Guide Vane Closing Schemes of Pumped Storage Hydro Unit Using an Enhanced Multi-Objective Gravitational Search Algorithm

The optimization of guide vane closing schemes (OGVCS) of pumped storage hydro units (PSHUs) is a cooperative control and optimal operation research field in renewable energy power generation technology. This paper presents an OGVCS model of PSHUs considering the rise rate of the unit rotational speed, the specific node pressure of each hydraulic unit, as well as various complicated hydraulic and mechanical constraints. The OGVCS model is formulated as a multi-objective optimization problem to optimize conflicting objectives, i.e., unit rotational speed and water hammer pressure criteria. In order to realize an efficient solution of the OGVCS model, an enhanced multi-objective bacterial-foraging chemotaxis gravitational search algorithm (EMOBCGSA) is proposed to solve this problem, which adopts population reconstruction, adaptive selection chemotaxis operator of local searching strategy and elite archive set to efficiently solve the multi-objective problem. In particular a novel constraints-handling strategy with elimination and local search based on violation ranking is used to balance the various hydraulic and mechanical constraints. Finally, simulation cases of complex extreme operating conditions (i.e., load rejection and pump outage) of a ‘single tube-double units’ type PSHU system are conducted to verify the feasibility and effectiveness of the proposed EMOBCGSA in solving OGVCS problems. The simulation results indicate that the proposed EMOBCGSA can provide a lower rise rate of the unit rotational speed and smaller water hammer pressure than other methods established recently while considering various complex constraints in OGVCS problems.


Introduction
In recent years, in order to improve the global greenhouse effect and carbon dioxide emissions, the use of wind power, solar power generation, biomass energy and other renewable energy sources has been expanding to maintain China's rapid development momentum [1,2].Grid installed capacity and new energy power generation has continued to grow, which effectively alleviates the pressure on economic development of coal, oil and other fossil energy dependence, contributes to the energy structure reform and the development of green energy.It is estimated that by the year 2020, the installed wind power, photovotaic and pumped storage capacity will reach 210 million kW, 110 million kW and 100 million kW, respectively.However, wind energy and solar energy are intermittent resources, which exacerbates the contradiction between the development and utilization of new energy sources.
As special energy storage power supplies, wind power-pumped storage plants (PSPs) and solar power-PSPs are the most commonly used centralized and large-scale renewable energy complementary operation means [3][4][5][6][7].The flexibility of PSP makes up for the randomness and heterogeneity of wind power and solar power generation, which is helpful to improve the reliability of the power grid and promote the integration of renewable energy forms.
The OGVCS is an important means to solve the regulation and guarantee calculation of the pumped storage power plant, and is also the preferred method to optimize the hydraulic transient process.During the dynamic process of of PSHUs with the non-optimal closing law under extreme conditions (i.e., load rejection or pump outage), the rise of the unit rotational speed and the water hammer pressure will exceed the maximum design value, which can cause the runaway of the PSHU, abnormal vibrations or active guide vane asynchronous and other negative phenomena [8,9].OGVCS plays an important role in ensuring the security and stability of the power grid, and several methods have been proposed to deal with the OGVCS problem.Based on the characteristics of the rigid water-column pressure during the transient process, [10] proposed a two-phase guide-vane closing scheme and three-phase valve-closing schemes to control the pulsating pressures and the runaway speed problems.Reference [11] analysed the effect of valve closure on the water-hammer pressure and laid the theoretical foundation for improving the turbine guide-vane closure.In order to avoid the operating point from being in the "S" characteristic area, Kuwabara et al. [12] proposed a curved closing scheme that could effectively decrease the water-hammer pressure.Additionally, appropriate two-phase closing schemes [13,14] and misaligned guide-vane methods [15][16][17] also have been applied to control the fluctuation of PSHUs.Three-phase valve-closing and curved closing schemes exert high demands on the governor servomotor, which is difficult to realize in a PSP.Meanwhile, the misaligned guide-vane method can significantly increase the pulsating pressure and the runner radial forces during the PSHU start-up process [17].Therefore, optimal closing schemes are required to realize optimal coordinated operation of PSP and new energy generation technologies based on a deep analysis of the operational mechanism and flow characteristics of PSHUs.
Taking into full consideration the complex hydraulic, mechanical and electrical coupling characteristics and nonlinear dynamic response process of PSHUs, a multi-objective optimization model of guide vane closing schemes is established in this paper, which includes hydraulic and mechanical multiple constraint factors.The OGVCS problem is a complex, multi-objective and multi-constraint optimization problem, which aims to reasonably and simultaneously control the rotational speed and water hammer pressure of volutes, draft tubes, surge tanks and so on.Multi-objective intelligent optimization algorithms are an effective way to solve complex multi-objective problems, and many intelligent optimization algorithms, such as particle swarm optimization (PSO), gravitational search algorithm (GSA) and so on, are widely used in solving such multi-objective optimization problems.Multi-objective particle swarm optimization (MOPSO) [18,19], non-dominated sorting genetic algorithm-II (NSGA-II) [20], multi-objective differential evolution (MODE) [21,22], multi-objective gravitational search algorithm (MOGSA) [23,24] and multi-objective bee colony optimization algorithm (MOBCO) [25,26], have all been proposed to solve complex multi-objective optimization problems with practical and efficient modeling of coupling constraints.However, premature minimization phenomena and local convergence are still common obstacles to the performance of these stochastic searching algorithms.The so-called bacterial-foraging chemotaxis gravitational search algorithm (BCGSA) is enhanced by Pbest-Gbest-guided movement, adaptive elastic-ball method and the chemotaxis operator strategy of the bacterial-foraging algorithm [27].The global exploration and local exploitation performance of BCGSA have been proved in [27].Inspired by the NSGA-II, EMOBCGSA has been introduced to deal with the OGVCS by reconstructing the optimal structure and method.Finally, cases of complex extreme operating conditions of 'single tube-double unit'-type PSHU systems are studied to verify the feasibility and effectiveness of the proposed EMOBCGSA in solving OGVCS problems.
The remainder of this paper is organized as follows.Section 2 introduces numerical calculation modeling of a "single tube-double unit"-type PSHU system.Section 3 establishes a multi-objective Energies 2017, 10, 911 3 of 23 optimization model of guide vane closing schemes and multiple constraints are introduced simultaneously.Section 4 delineates the general procedure and multi-objective improvements in BCGSA.Section 5 illustrates the practical solving procedure for OGVCS problems and experimental results, along with a few discussions, respectively.The conclusions are summarized in Section 6.

Numerical Calculation Model of the Pumped Storage Hydro Unit System
The "single tube-double unit"-type PSHU system is a typical kind of power generation system that couples hydraulic, mechanical and electrical factors, simplified mainly into four parts, namely reservoir, pressure water pipeline, surge tank and pumped storage hydro unit.Figure 1 shows a schematic diagram of the "single tube-double units"-type PSHU system used in this paper, where the PSHU system is divided into upstream and downstream reservoirs, eight pressure water pipelines, upstream and downstream surge tanks, and two PSHUs.In this section, mathematical models of each link are described.
Energies 2017, 10, 911 3 of 23 BCGSA.Section 5 illustrates the practical solving procedure for OGVCS problems and experimental results, along with a few discussions, respectively.The conclusions are summarized in Section 6.

Numerical Calculation Model of the Pumped Storage Hydro Unit System
The "single tube-double unit"-type PSHU system is a typical kind of power generation system that couples hydraulic, mechanical and electrical factors, simplified mainly into four parts, namely reservoir, pressure water pipeline, surge tank and pumped storage hydro unit.Figure 1 shows a schematic diagram of the "single tube-double units"-type PSHU system used in this paper, where the PSHU system is divided into upstream and downstream reservoirs, eight pressure water pipelines, upstream and downstream surge tanks, and two PSHUs.In this section, mathematical models of each link are described.

Model of Pressure Water Pipeline
The well-known characteristics method is employed in the modeling of pressure water pipelines.The basic motion equation and continuous unsteady flow equation in a pressure pipeline can be expressed by Equations ( 1) and (2) [28].Equations ( 1) and ( 2) can be solved using the method of characteristics, which transforms the formulas into a simplified equation set in the range of the characteristic line dx a dt   .The simplified equation set is described as by Equations ( 3) and (4): Motion equation: Continuous equation: where V is the flow velocity, H is the piezometric head, a is the water hammer velocity, f is the friction coefficient, D is the pipeline diameter,  is the angle between the line and a horizontal surface of each section of the pipeline center, V and H are functions of the pipe length x and time t , respectively, A is the pipe section area, and Q is the water flow in the pipe section.
The difference network is constructed using the above simplified equation set, then the finite difference method is used to solve the problem.Figure 2 shows the characteristic line difference mesh.

Model of Pressure Water Pipeline
The well-known characteristics method is employed in the modeling of pressure water pipelines.The basic motion equation and continuous unsteady flow equation in a pressure pipeline can be expressed by Equations ( 1) and (2) [28].Equations ( 1) and ( 2) can be solved using the method of characteristics, which transforms the formulas into a simplified equation set in the range of the characteristic line dx dt = ±a.The simplified equation set is described as by Equations ( 3) and (4): Motion equation : Continuous equation : where V is the flow velocity, H is the piezometric head, a is the water hammer velocity, f is the friction coefficient, D is the pipeline diameter, α is the angle between the line and a horizontal surface of each section of the pipeline center, V and H are functions of the pipe length x and time t, respectively, A is the pipe section area, and Q is the water flow in the pipe section.The difference network is constructed using the above simplified equation set, then the finite difference method is used to solve the problem.Figure 2 shows the characteristic line difference mesh.For Figure 2, the length of L pipeline is divided into N segments, and the length of each segment is ∆x = L/N, while the time step of the differential network is ∆t = ∆x/a.
For Figure 2, the length of L pipeline is divided into N segments, and the length of each segment is , while the time step of the differential network is / t x a    .Diagonal AP satisfies the condition x a t    and diagonal BP satisfies the condition x a t     .Equations ( 5) and ( 6) are the integral of motion equation between A and P and the integral of continuous equation between B and P, respectively: Finally, the correct processing and analog computation of boundary conditions for the PSHU system in the transient process of the research and control are very important.The boundary conditions are a conduit or a connected terminal, including upstream reservoir, upstream surge tank, series nodes, bifurcated pipe nodes, ball valve, PSHU, downstream surge tank and downstream reservoir.Besides the surge tank and PSHU, the other boundary conditions as well as the spiral case and tail pipe usually can be considered as pressure water pipelines.

Surge Tank Model
There are many types of surge tanks in the built pumped storage power stations, including impedance type, differential type and air cushion type.Based on the actual situation of pumped storage power stations in China, the impedance surge tank is adopted by the PSHU system in this paper.The schematic diagram of the impedance surge tank is shown in Figure 3, from which the impedance surge tank is connected with the water pipeline system through a small impedance orifice, which has the advantages of small volume and simple structure.The corresponding basic equations can be described by Equation (7) [29][30][31]: where S H , S Q are the surge tank bottom pressure and flow, 1 A is the area of the impedance orifice, C H , 2 A are the elevation and area of the surge tank, R K is the bottom orifice flow loss coefficient.Diagonal AP satisfies the condition ∆x = a∆t and diagonal BP satisfies the condition ∆x = −a∆t.Equations ( 5) and ( 6) are the integral of motion equation between A and P and the integral of continuous equation between B and P, respectively: Finally, the correct processing and analog computation of boundary conditions for the PSHU system in the transient process of the research and control are very important.The boundary conditions are a conduit or a connected terminal, including upstream reservoir, upstream surge tank, series nodes, bifurcated pipe nodes, ball valve, PSHU, downstream surge tank and downstream reservoir.Besides the surge tank and PSHU, the other boundary conditions as well as the spiral case and tail pipe usually can be considered as pressure water pipelines.

Surge Tank Model
There are many types of surge tanks in the built pumped storage power stations, including impedance type, differential type and air cushion type.Based on the actual situation of pumped storage power stations in China, the impedance surge tank is adopted by the PSHU system in this paper.The schematic diagram of the impedance surge tank is shown in Figure 3, from which the impedance surge tank is connected with the water pipeline system through a small impedance orifice, which has the advantages of small volume and simple structure.The corresponding basic equations can be described by Equation (7) [29][30][31]: where H S , Q S are the surge tank bottom pressure and flow, A 1 is the area of the impedance orifice, H C , A 2 are the elevation and area of the surge tank, K R is the bottom orifice flow loss coefficient.

Pumped Storage Hydro Unit Model
The PSHU is composed of a pump-turbine and generator, which are the key components of a PSHU system.At present, pump-turbine modeling based on the characteristic curves has been widely used.In this paper, the pump-turbine is expressed with the torque function and flow function for state variables, including guide vane opening, generator rotational speed and water head, shown as Equation ( 8) [32].The analytic expressions of the nonlinear functions of ( ) to obtain, however, based on the characteristic curves (Figure 4), torque and flow of the pump-turbine at a certain time can be calculated by means of interpolation or nonlinear function fitting.In order to overcome the obstacle of single input-multiple output during the interpolation calculation process, the logarithmic-curve-projection method is adopted for the mathematical transformation of the characteristic curves [33]: where t m and q indicate the torque and flow of the pump-turbine, respectively.In this paper, only the PSHU rotational speed change of standalone running in the isolated power grid is considered [34].The dynamic equation of the synchronous generator considering the load characteristics is simplified to give Equation ( 9):

Pumped Storage Hydro Unit Model
The PSHU is composed of a pump-turbine and generator, which are the key components of a PSHU system.At present, pump-turbine modeling based on the characteristic curves has been widely used.In this paper, the pump-turbine is expressed with the torque function and flow function for state variables, including guide vane opening, generator rotational speed and water head, shown as Equation ( 8) [32].The analytic expressions of the nonlinear functions of f m (•) and f q (•) are difficult to obtain, however, based on the characteristic curves (Figure 4), torque and flow of the pump-turbine at a certain time can be calculated by means of interpolation or nonlinear function fitting.In order to overcome the obstacle of single input-multiple output during the interpolation calculation process, the logarithmic-curve-projection method is adopted for the mathematical transformation of the characteristic curves [33]: where m t and q indicate the torque and flow of the pump-turbine, respectively.

Pumped Storage Hydro Unit Model
The PSHU is composed of a pump-turbine and generator, which are the key components of a PSHU system.At present, pump-turbine modeling based on the characteristic curves has been widely used.In this paper, the pump-turbine is expressed with the torque function and flow function for state variables, including guide vane opening, generator rotational speed and water head, shown as Equation ( 8) [32].The analytic expressions of the nonlinear functions of ( ) to obtain, however, based on the characteristic curves (Figure 4), torque and flow of the pump-turbine at a certain time can be calculated by means of interpolation or nonlinear function fitting.In order to overcome the obstacle of single input-multiple output during the interpolation calculation process, the logarithmic-curve-projection method is adopted for the mathematical transformation of the characteristic curves [33]: where t m and q indicate the torque and flow of the pump-turbine, respectively.In this paper, only the PSHU rotational speed change of standalone running in the isolated power grid is considered [34].The dynamic equation of the synchronous generator considering the load characteristics is simplified to give Equation ( 9): In this paper, only the PSHU rotational speed change of standalone running in the isolated power grid is considered [34].The dynamic equation of the synchronous generator considering the load characteristics is simplified to give Equation ( 9): Energies 2017, 10, 911 6 of 23 where x is the generator rotational speed, T a is the inertial time constant of the generator, e g is the adjusting coefficient of the generator, and mg 0 is the load disturbance which represents the load changes.

Optimization of Guide Vane Closing Schemes Problem Formulation
The dynamic response characteristics of the unit and the hydraulic influence of the large fluctuation condition of the PSHU are analyzed under the conditions of different guide vane closing laws.Considering the increase of unit rotational speed and pressure of each node in the hydraulic unit as the two objectives, the optimal model of the guide vane closing law is established.

Rotational Speed Objective
The rise of the rotational speed is one of the main characteristics of the pumped storage units under load rejection or pump power off conditions.Under load rejection conditions, a non-optimal guide vane closure law will lead to the runaway phenomenon, causing mechanical damage, vibrations and noise in the PSHU.The minimize rotational speed objective Obj x of the OGVCS problem can be expressed as follows: where N pu is the number of units in a hydraulic unit, x i is the rotational speed of i unit during the transition process, x r,i is the rated rotational speed of the unit under stable conditions.

Water Hammer Pressure Objective
In order to minimize the pressure of each hydraulic unit, comprehensively considering the volute end pressure, tail pipe inlet pressure and surge tank water level value, the water hammer pressure objective function Obj pre is established as follows: where P vol_e,i , P dra_s,i are the maximum value of the volute end pressure and tail pipe inlet pressure for i unit, L sur_up , L sur_down are the maximum water level of the upstream surge tank and downstream surge tank, w P , w L are the weight coefficients of the water hammer pressure and water level in the surge tank.

Multiple Constraints
The OGVCS problem should satisfy the following equality and inequality constraints: (1) Rising rotational speed constraint To guarantee regulation calculation and transient analysis of a PSHU system, there is a clear maximum constraint value constant x for the rise value of rotational speed in all kinds of extreme cases: (2) Limit of rotational speed fluctuation In this paper, the limit of the times of rotational speed fluctuation is introduced in order to achieve a good dynamic quality of the large fluctuation condition process: where constant obj xr is the rising rate constant for dynamic quality requirements, N x f is the number of fluctuations, and constant x f is constraint constant of the number of fluctuations.
(3) Volute pressure constraint The maximum corrected value and the pressure constraint of the volute inlet pressure considering the pressure fluctuation and the calculation error are described as follows: where P vol_s,i , P vol,i and Pm vol_s,i are the maximum calculated value, initial value and maximum corrected value of the volute inlet pressure for i unit, H n is the net head, and constant Pm vol_s is the constraint constant of P vol_s,i .
(4) Draft tube inlet pressure constraint where P dra_s,i , P dra,i and Pm dra_s,i are minimum calculated value, initial value and minimum corrected value of the draft tube inlet pressure for i unit.
(5) The surge water level limits where L sur_up and l sur_up are the maximum and minimum values of the surge water level of the upstream surge tank, L sur_down and l sur_down are maximum and minimum values of the surge water level of the downstream surge tank, and the other terms are the corresponding constraint constants.
(6) Speed governor oil velocity limit If the guide vanes are closed within a short time, the control accuracy of the governor curve slope has higher requirements.In order to meet such a short closing time demand, the servomotor oil speed must be large enough, which can cause major safety issues.Therefore, this paper takes into account the limiting factor of the slope control of the governor and transforms it into the guide vane closing rate: where ∆Y and Y _max are change value and rated maximum value guide vane opening, t is the guide vane closing time, and Tr is the minimum closing time limit.

Overview of Bacterial-Foraging Chemotaxis Gravitational Search Algorithm
BCGSA is an improvement of the standard gravitational search algorithm (GSA) by incorporating the P best -G best -guided movement, adaptive elastic-ball method and chemotaxis operator strategies.Then the specific steps of BCGSA are as follows: Step 1: Initialization.Randomly initialize the agent of population position and velocity.
Step 2: Fitness calculation and information save.Calculate the fitness of agents according to their initial position, storing the current position of each agent P best (t) as the best record position of the agent and position of best agent G best (t) as the best position in global.
Step 4: Calculate the gravitational constant in the current iteration and acceleration for each agent.
Step 5: Update agents' velocity v d i and position x d i with equations ( 18) and ( 19) [35].
where r 1 , r 2 and r 3 are random variables in the range [0, 1], c 1 and c 2 are learning genes in the range [0, 2], P ibest (t) is the best position that i-th agent has ever suffered until time t, G best (t) is the global best position in the agents until time t.
Step 6: Judge whether the new position of the agent is beyond the boundary.Invoke the adaptive elastic-ball program (Equation ( 20)) to handle the new position which is against the boundary.Besides, there are a few rest agents against the boundary whose position will be reset by Equation (21) [36]: where ζ is the adaptive attenuation coefficient of reflective power, Ub(d) and Lb(d) are upper and lower boundary limit in the d-th dimension, respectively.Step 7: Evaluate the fitness in accordance with the new position of each agent, storing position of the best agent X best and the worst agent X worst in the current iteration.
Step 8: The chemotaxis operator [37] is applied for X best and X worst .] according to Equations ( 22) and (23).
x(i, j + 1) = x(i, j) + C(i)φ(j) Energies 2017, 10, 911 where X best is the agent of the best fitness, X worst is the agent of the worst fitness, x(i, j) is the position of the i-th agent in the j-th chemotaxis.(d) Swim.Compute fitness of obtained new agent x(i, j + 1), then compare f itness(i, j + 1) and f itness(i, j).If the f itness(i, j + 1) is better, save the new agent as X best new or X worst new and start agent swimming, then let m = 0 and repeat Equation (23)

The Basic Definition of Multi-Objective Optimization Problem
Generally, the multi-objective minimize optimization problem of n dimensions decision variable and m dimensions subunit objectives can be described as follows: where [ f 1 , f 2 , . . ., f m ] are m dimensions of objective vectors, F(x) is a mapping function from decision space to object space, g j (x) ≥ 0 (j = 1, 2, . . ., J) defines the J inequality constraints, h k (x) = 0 (k = 1, 2, . . ., K) defines the K inequality constraints, x = (x 1 , x 2 , . . ., x n ) is the m dimensions decision variable space.
Due to the conflictive relationship among multiple subobjectives, a set of non-dominated solutions exist rather than only one optimal solution.Thus, the key concept of Pareto optimality [38,39] associated with the trade-off between conflictive objectives is adopted: (1) Pareto dominance relation X f is a feasible solution set, x 1 and x 2 are feasible solutions, [x 1 , x 2 ] ∈ X f , if x 1 dominates x 2 , note as x 1 x 2 , if and only if: (2) Pareto optimal solution For a given multi-objective optimal problem, if x * ∈ X f is the Pareto optimal solution, if and only if: (3) Pareto optimal solution set The set of all Pareto optimal solutions is called the Pareto optimal solution set P * : (4) Pareto optimal front Energies 2017, 10, 911 10 of 23 The Pareto optimal solution set corresponding to the objective vector value in the objective domain space is called the Pareto optimal front PF * : The BCGSA is based on the single objective optimization design, which does not have the ability to deal with multi-objective optimization problems.Therefore, based on the BCGSA, this section reconstructs the optimal structure and mechanism of solving multi-objective optimization problems.

Population Reconstruction
In dealing with multi-objective optimization problems, it is more complex to determine the non-inferiority relation among population individuals and it is the key to improve the efficiency and ergodicity of the algorithm to establish a reasonable non-inferiority relation among individuals.Therefore, this paper further introduces a fast non-dominated sorting method, the crowding-distance concept, and then completes the whole population structure reconstruction.
(1) Non-dominated sorting In order to realize the distribution and diversity of population of the BCGSA algorithm in dealing with multi-objective optimization problems in this group, inspired by NSGA-II [20], the fast non-dominated sorting method is applied to divide the population into a number of levels, the specific steps are as follows: Assuming that the population size is 7, the objective number is 2, and the layering schematic is shown in Figure 5.

Enhanced Multi-Objective Bacterial-Foraging Chemotaxis Gravitational Search Algorithm
The BCGSA is based on the single objective optimization design, which does not have the ability to deal with multi-objective optimization problems.Therefore, based on the BCGSA, this section reconstructs the optimal structure and mechanism of solving multi-objective optimization problems.

Population Reconstruction
In dealing with multi-objective optimization problems, it is more complex to determine the noninferiority relation among population individuals and it is the key to improve the efficiency and ergodicity of the algorithm to establish a reasonable non-inferiority relation among individuals.Therefore, this paper further introduces a fast non-dominated sorting method, the crowding-distance concept, and then completes the whole population structure reconstruction.
(1) Non-dominated sorting In order to realize the distribution and diversity of population of the BCGSA algorithm in dealing with multi-objective optimization problems in this group, inspired by NSGA-II [20], the fast non-dominated sorting method is applied to divide the population into a number of levels, the specific steps are as follows: Step  (2) Crowding-distance calculation Crowding-distance is an important index to evaluate the distribution density of solution space, which can be described as follows: Crowding-distance is an important index to evaluate the distribution density of solution space, which can be described as follows: where Γ[p] distance is the crowding-distance of individual p, f p,l is the fitness of objective l of individual p, L is the number of objective, f max l and f min l are the maximum and minimum fitness for individuals in each layer of the L dimensional objective.
Thus, each individual consists of two attributes: the layer number and the crowding-distance.Depending on these two attributes can be defined between individual preference relations.i is superior to j, if and only if: (3) Individual quality calculation For EMOBCGSA, the quality of the individual cannot be calculated on the value of the objective function, but it is obtained by the rank of the non-dominated layer number.For l objective functions of multi-objective optimization problem, the calculation process of individual quality given as follows: Step 1: In the rank = 1, the fitness values of the individuals with the smallest value of l objective functions and the maximum crowding-distances are set to 1, and the fitness value of other individuals is 2.
Step 2: In the rank = 2, set all the individual fitness values to 3.
Step 3: Repeat Step 2, individual fitness value of rank ≥ 3 is set to (rank + 1), thus completing all layers of individual fitness assignment.

Multi-Objective Adaptive Chemotaxis Operation
For multi-objective adaptive chemotaxis operations, the selections of G best and G worst are the key to speeding up the convergence of EMOBCGSA.In order to avoid the premature convergence and maintain the diversity of the population, the concept of constraint violation, as shown in Equation (32), is introduced in the selection process.Combined with individual crowding-distance, the probability of individual selection in rank = 1 and rank = max is calculated by Equation (33).By means of Equations ( 32) and (33) to adaptively calculate the probability of individual selection, then the tracking targets (G best and G worst ) of multi objective chemotaxis operation are determined by the roulette wheel method: where voilation s is the constraint violation of individual s, value ic is the value of s in ic constraint, Lu ic and Ld ic are the upper and lower limits of constraint ic, nc is the number of constraints; p s is the selection probability of s, ε is the feasible margin of constraint failure depth evaluation (if voilation s > ε, and individual s is infeasible solution), N v and N d are the number of infeasible solutions and feasible solutions, Sp1 and Sp2 are the initial selected probability constants of the infeasible and feasible solutions.

Multi-Objective Velocity Update Strategy
For the velocity update Equation ( 18) of EMOBCGSA, G best and P best need to be redefined when dealing with multi-objective problems.The selection mechanism of G best has been introduced in detail in Section 4.3.2, and P d ibest (t) is the location of the optimal memory of the individual i on the d dimension of the t th generation.The selection method of P best can be carried out according to Equation (30).

The Update and Maintenance of the Elite Archive Set
The elite archive set (EAS) is used to store the non-dominated solutions obtained during EMOBCGSA evolution, and provide guidance for the evolution of a population.In the actual operation of the algorithm, the newly generated individual np (1) If the current capacity of EAS is zero, the elite candidates are added directly to the EAS.
(2) If the EAS is not empty, and the elite candidate is not dominated by any elite individual in the original EAS, the elite candidate is added to the EAS.At the same time, the initial elite individual who is dominated by the candidate elite is removed from the EAS.(3) If the number of elite individuals stored in EAS has exceeded NQ, the EliteSet truncation method is used to maintain it.
All the individuals will be selected by the non-dominated sorting, and the first NQ individuals with larger crowding-distance will be stored in EAS, where the redundant individuals will be removed.

Performance Test of Enhanced Multi-Objective Bacterial-Foraging Chemotaxis Gravitational Search Algorithm
In order to verify the performance of the EMOBCGSA in dealing with multi-objective problems, the test function set of ZDT is selected to test and evaluate the performance of the proposed algorithm [40].In addition, the convergence and diversity indexes are introduced to evaluate the performance of the algorithm, respectively, such as Equations ( 34) and (35).
where N is the size of the Pareto solution set, x is the arbitrary nondominated solution of the Pareto solution set, x is the nearest pairing solution from the x on the real Pareto front, d i is the Euclidean distance between two consecutive non-dominated solutions, d is the average value of all d i , d b and d e respectively show that the algorithm obtains the Euclidean distance between the ends of the Pareto solution and the real Pareto front.The convergence and diversity index of the test function set are calculated, and the results are compared with those of NSGA-II, MOPSO and multi objective evolutionary algorithm (MOEA).The mean and variance of the two indexes after the algorithms were run 20 times are shown in Tables 1  and 2. From Table 1, we can see that the EMOBCGSA test results are better than those of the other methods in convergence index, mean and variance, and the algorithm has strong convergence and robustness.In the distribution of diversity of non-dominated solutions, the results of EMOBCGSA in Table 2 are superior to the other optimization algorithms in terms of diversity index, mean and variance.The variance is small, and the calculation results are also stable.

Proposed Enhanced Multi-Objective Bacterial-Foraging Chemotaxis Gravitational Search Algorithm Approach for Optimization of Guide Vane Closing Schemes Problem
Depending on the mathematical description of the optimization objectives and constraints, OGVCS is a multi-objective optimization problem with multiple variables and multiple nonlinear constraints.The objectives and constraints involve a lot of hydraulic and mechanical factors, so it is important to have an efficient algorithm to solve the problem.In this section, the procedure of the proposed EMOBCGSA for solving OGVCS problems with complex constraints is described, and the corresponding constraints-handling strategy will be discussed later.The generalized procedure can be given as follows.From Figure 6, T and [ , , ] Tf Y Tt are the decision variables of the two kinds of closing law, respectively.Taking the broken line closure as an example and assuming that the optimized 'single tube-PT N units' type PSHU system, the individual a can be constructed as matrix a U with , a U is described as Equation ( 36), Np is the size of population.Individuals can be randomly initialized by Equation (37): where

Constraint-Handling Strategy
In this section, based on the concept of constraint violation, a novel constraints-handling strategy with elimination and local search based on violation ranking is proposed.The constraint-handling strategy is used to balance the various hydraulic and mechanical constraints of the OGVCS problem.The speed governor oil velocity limit of the governor can be dealt with after the initialization and updating of the individual, and the additional constraint is to judge whether the individual has violated the constraint condition after the OGVCS problem is solved.Therefore, the constrainthandling strategy can be divided into two parts: pre-correction and post-correction.
(1) Pre correction Pre-repair is mainly aimed at tackling the speed governor oil velocity constraints.The processing steps are as follows: From Figure 6, T and [T f , Y, Tt] are the decision variables of the two kinds of closing law, respectively.Taking the broken line closure as an example and assuming that the optimized 'single tube-N PT units' type PSHU system, the individual a can be constructed as matrix U a with 3 × N PT , U a is described as Equation ( 36), N p is the size of population.Individuals can be randomly initialized by Equation (37): where T f max 1,n and T f min 1,n are the upper and lower limits of T f , Y max 1,n and Y min 1,n are the upper and lower limits of Y, Tt max 1,n and Tt min 1,n are the upper and lower limits of Tt.

Constraint-Handling Strategy
In this section, based on the concept of constraint violation, a novel constraints-handling strategy with elimination and local search based on violation ranking is proposed.The constraint-handling strategy is used to balance the various hydraulic and mechanical constraints of the OGVCS problem.The speed governor oil velocity limit of the governor can be dealt with after the initialization and updating of the individual, and the additional constraint is to judge whether the individual has violated the constraint condition after the OGVCS problem is solved.Therefore, the constraint-handling strategy can be divided into two parts: pre-correction and post-correction.
(1) Pre correction Pre-repair is mainly aimed at tackling the speed governor oil velocity constraints.The processing steps are as follows: Step 1: For individual i, calculate the closing rate for each segment ∆Y/t, judge the relationship between ∆Y/t with the rated closing rate Y_max/Tr.Step 2: If ∆Y/t > Y_max/Tr, the closing rate of the guide vanes assigned to this segment is Y_max/Tr, reverse calculation closing time t_new, then replace t with t_new.
Step 3: Repeat Steps 1 and 2, until all of the individuals constraints have been repaired.
(2) Post-correction The OGVCS problem is solved by the individuals who have been pre-repaired.According to the calculation value of the transient process of each hydraulic unit of PSHU system and Equation (32), we can calculate the total amount of constraint violation voilation i .If voilation i = 0, according to eliminate and local search strategy based on violation ranking for post repair, detailed steps are as follows: Step 1: Conduct a descending sort for individuals of voilation i = 0.
Step The process of OGVCS problem solving with EMOBCGSA is described as follows, and the flowchart is illustrated in Figure 7: Step 1: Set guide vane closing mode and corresponding parameters; initialize the control parameters [NQ, N p, G max , G 0 , β, c 1 , c 2 , Nc, Ns] of the EMOBCGSA, set the evolutionary current algebra g = 1.
Step 2: According to Equations (34) and (35) construct and initialize the individuals of a population.
Step 3: Complete the constraints pre-repair of individuals; simulated solutions of the OGVCS problem, Obj x and Obj pre are calculated by Equations ( 10) and (11), respectively; complete the constraints post-repair; the non-dominated sorting of the population is conducted to determine the number of layer and crowding-distance of the population, and all the non-dominated individuals of the rank = 1 layer to meet the constraints that were added to EAS.Step 5: If g < G max , set g = g + 1, repeat Steps 3 and 4, and update and maintenance of EAS; otherwise output the current EAS as the Pareto optimal front, OGVCS problem calculation completed.

Experiment Scenarios and Parameters Setting
In this section, the "single tube-double unit"-type PSHU system model is constructed to verify the feasibility and effectiveness of EMOBCGSA for solving the OGVCS problem.The model parameters are calculated on the basis of the actual parameters of a Chinese pump storage plant (PSP).The actual physical parameters of the PSP and PSHU are given in Table 3.The flow and torque characteristic curves of the TPV32-LJ-385 pump-turbine are shown in Figure 4. Additionally, the simulation conditions of the experiment are set as shown in Table 4. Four scenarios of guide vane closing modes are used for model validation experiments.

Experiment Scenarios and Parameters Setting
In this section, the "single tube-double unit"-type PSHU system model is constructed to verify the feasibility and effectiveness of EMOBCGSA for solving the OGVCS problem.The model parameters are calculated on the basis of the actual parameters of a Chinese pump storage plant (PSP).The actual physical parameters of the PSP and PSHU are given in Table 3.The flow and torque characteristic curves of the TPV32-LJ-385 pump-turbine are shown in Figure 4. Additionally, the simulation conditions of the experiment are set as shown in Table 4. Four scenarios of guide vane closing modes are used for model validation experiments.Based on the above results analysis, the proposed method can get the optimal solution set, which makes the PSHU run safely and stably, and the water diversion system has the smallest fluctuation, which provides the best decision space for the decision maker.Additionally, in order to verify the availability of the proposed constraints-handling strategy, the best Objpre and Objx schemes of scenario 3 is taken as the compromise schemes in our study.The Based on the above results analysis, the proposed method can get the optimal solution set, which makes the PSHU run safely and stably, and the water diversion system has the smallest fluctuation, which provides the best decision space for the decision maker.Additionally, in order to verify the availability of the proposed constraints-handling strategy, the best Objpre and Objx schemes of scenario 3 is taken as the compromise schemes in our study.The Additionally, in order to verify the availability of the proposed constraints-handling strategy, the best Obj pre and Obj x schemes of scenario 3 is taken as the compromise schemes in our study.The simulation results of the PSHU transition process are shown in Figure 10.It can be seen clearly from the chart that the transition process of the corresponding scheme satisfies the above definition of multiple constraints and has a good dynamic response process of the PSHU.The corresponding index values of scenario 3 are shown in Tables 6 and 7. From these tables, it is evident that the maximum water hammer pressure of each hydraulic unit and the rotational speed rise value can meet the requirements of the adjustment and guarantee calculation.In particular, the minimum Obj x is only 27.5%.6 and 7. From these tables, it is evident that the maximum water hammer pressure of each hydraulic unit and the rotational speed rise value can meet the requirements of the adjustment and guarantee calculation.In particular, the minimum Objx is only 27.5%.

Pump Outage Condition
Under the pump outage condition, comparative analysis of EMOBCGSA Pareto non-dominated set by the four scenarios, is shown in Figure 11.From Figure 11, it is evident that the non-dominated solution set obtained in scenario 4 is closer to the true Pareto optimal front, and has the sets of best Objpre, Objx and the most balanced of the two objectives.In Table 8, the schemes at the inflection point of the Pareto optimal front by the four scenarios for the comparative analysis are presented (15  Under the pump outage condition, comparative analysis of EMOBCGSA Pareto non-dominated set by the four scenarios, is shown in Figure 11.From Figure 11, it is evident that the non-dominated solution set obtained in scenario 4 is closer to the true Pareto optimal front, and has the sets of best Obj pre , Obj x and the most balanced of the two objectives.In Table 8, the schemes at the inflection point of the Pareto optimal front by the four scenarios for the comparative analysis are presented (15 schemes in scenario 1, 14 schemes in scenario 2, 17 schemes in scenario 3 and 22 schemes in scenario 4).The simulation results of the PSHU transition process under pump outage conditions and the corresponding index values are presented in Figure 12 and Tables 9 and 10.From the guide vane closing laws of Figure 12a and index values, the draft inlet pressure increases at first and then decreased and becomes stabilized; The rotational speed extreme value depends on the guide vane closing time; if the closing time is short, the unit will only run between the pump operating area and pump brake operating area, and encounter no rotating speed reversal phenomenon.Therefore, the closing time of the guide vane is extremely important to the transient process of the PSHU under pump outage conditions.

Figure 1 .
Figure 1.Structure of a 'single tube-double unit'-type pumped storage hydro unit system.

Figure 1 .
Figure 1.Structure of a 'single tube-double unit'-type pumped storage hydro unit system.

Figure 3 .
Figure 3. Schematic diagram of the impedance surge tank.

Figure 3 .
Figure 3. Schematic diagram of the impedance surge tank.
parameters Nc, Ns and C(i).Where Nc is the number of chemotaxis steps, Ns is the number of swim steps.(b)Tumble.Generate a random vector∆(i) ∈ R D with each element ∆ d (i), d = 1,2, 3, . . ., D, a random number on [−1, 1].(c) Move.Define unit length random direction φ(j) in j-th chemotaxis, then updates the best agent X best = [x best 1 , x best 2 , . . ., x best d , . . ., x best D ] and worst agent X worst =

Step 1 :
If any individual x d i is not dominant in the current population, the individual is Pareto non-dominated, then its layer rank is 1, and all the individuals in the rank = 1 constitute the Pareto optimal front.Step 2: In the whole population, the individuals of rank = 1 are excluded; then repeat Step 1 until obtain the individual solution set of the second layer, and the individual attributes of the rank = 2 are given.Step 3: By analogy, repeat Step 2 until the completion of the entire population of non-dominated type.

1 :
If any individual d i x is not dominant in the current population, the individual is Pareto non-dominated, then its layer rank is 1, and all the individuals in the rank = 1 constitute the Pareto optimal front.Step 2: In the whole population, the individuals of rank = 1 are excluded; then repeat Step 1 until obtain the individual solution set of the second layer, and the individual attributes of the rank = 2 are given.Step 3: By analogy, repeat Step 2 until the completion of the entire population of nondominated type.Assuming that the population size is 7, the objective number is 2, and the layering schematic is shown in Figure 5.

Figure 5 .
Figure 5.The layering schematic of population.

Figure 5 .
Figure 5.The layering schematic of population.
reciprocal domination.EAS update and maintenance strategies are performed as follows (NQ is the capacity of EAS): At present, the common guide vane closing mode is straight line closure, two segment type broken line closure and three segment delay broken line closure.Because the governor of PSHU does not have the function of delay in the pump condition, this paper does not consider the three segment closing law.Figure6shows the schematic diagram of the straight line closure and the broken line closure.Energies 2017, 10, 911 14 of 235.1.1.Population InitializationAt present, the common guide vane closing mode is straight line closure, two segment type broken line closure and three segment delay broken line closure.Because the governor of PSHU does not have the function of delay in the pump condition, this paper does not consider the three segment closing law.Figure6shows the schematic diagram of the straight line closure and the broken line closure.

Figure 6 .
Figure 6.Schematic diagram of the straight line closure and the broken line closure.
Ttare the upper and lower limits of Tt .

Figure 6 .
Figure 6.Schematic diagram of the straight line closure and the broken line closure.

Figure 7 .
Figure 7.The flowchart of EMOBCGSA for solving OGVCS problem.

Figure 7 .
Figure 7.The flowchart of EMOBCGSA for solving OGVCS problem.

2 /Figure 9 .
Figure 9.The Pareto front obtained by different scenarios under load rejection.

2 /Figure 9 .
Figure 9.The Pareto front obtained by different scenarios under load rejection.

Figure 9 .
Figure 9.The Pareto front obtained by different scenarios under load rejection.

Energies 2017 ,
10, 911 19 of 23 simulation results of the PSHU transition process are shown in Figure 10.It can be seen clearly from the chart that the transition process of the corresponding scheme satisfies the above definition of multiple constraints and has a good dynamic response process of the PSHU.The corresponding index values of scenario 3 are shown in Tables

Figure 10 .
Figure 10.The transient process of minimum Objx and Objpre schemes of scenario 3: (a) guide vane closing laws; (b) rotational speed; (c) volute end pressure; (d) draft tube inlet pressure; (e) upper surge water level and (f) down surge water level.

Figure 10 .
Figure 10.The transient process of minimum Obj x and Obj pre schemes of scenario 3: (a) guide vane closing laws; (b) rotational speed; (c) volute end pressure; (d) draft tube inlet pressure; (e) upper surge water level and (f) down surge water level.

Energies 2017 ,
schemes in scenario 1, 14 schemes in scenario 2, 17 schemes in scenario 3 and 22 schemes in scenario 4).The simulation results of the PSHU transition process under pump outage conditions and the corresponding index values are presented in Figure12and Tables9 and 10.From the guide vane closing laws of Figure12aand index values, the draft inlet pressure increases at first and then decreased and becomes stabilized; The rotational speed extreme value depends on the guide vane closing time; if the closing time is short, the unit will only run between the pump operating area and pump brake operating area, and encounter no rotating speed reversal phenomenon.Therefore, the closing time of the guide vane is extremely important to the transient process of the PSHU under pump outage conditions.

Figure 11 .
Figure 11.The Pareto front obtained of EMOBCGSA by different scenarios under pump outage conditions.

Figure 11 .
Figure 11.The Pareto front obtained of EMOBCGSA by different scenarios under pump outage conditions.
following this tumble direction until m = Ns; else let n=Ns directly, this is the end of the while statement.(e) Repeat (a)-(d) until Nc reaches the stop criteria.Replace G best and G worst by X new , only when X best new and X worst new are better.Step 9: Compare the obtained fitness of a new position x i (t + 1) with fitness of P ibest (t) and G best (t) while i changes from 1 to N. If x i (t + 1) has a better fitness value, replace position of P ibest (t) and G best (t) by x i (t + 1).Step 10: Repeat Steps 3-9 until the stop criterion is reached.
2: Introduce coefficient Cv of constraint violation, if voilation i < Cv, start the local iterative search procedure.Then, randomly generated ∆y i , and the position of the inflection point of two segment type broken line closures is updated to Y i + ∆y i , and we calculate the OGVCS problem.Repeat the above local iterative search procedure until voilation i = 0 or the required local search steps are met.Step 3: Repeat Step 1, complete local iterative search for individuals who are voilation i < Cv, eliminate individuals who is voilation i = 0. 5.1.3.Outline of Enhanced Multi-Objective Bacterial-Foraging Chemotaxis Gravitational Search Algorithm for Solving Optimization of Guide Vane Closing Schemes Problems Calculating the G(t), P ibest (t) and G best (t).The M i (t) of the individual is calculated by the number of ranks and crowding-distance.(2) Multi-objective chemotaxis operation.Determine the tracking target G best (t) and G worst (t) of the multi-objective chemotaxis operation, update G best (t) and G worst (t)

Table 3 .
Actual physical parameters of the pump storage plant and pumped storage hydro unit.

Table 3 .
Actual physical parameters of the pump storage plant and pumped storage hydro unit.O and 13.6 mH 2 O at best Obj pre , 19.17%, 17.49% and 1.01% at best Obj x Based on the above results analysis, the proposed method can get the optimal solution set, which makes the PSHU run safely and stably, and the water diversion system has the smallest fluctuation, which provides the best decision space for the decision maker.

Table 6 .
The index values of the minimum Objx and Objpre schemes of scenario 3 (I).

Table 7 .
The index values of the minimum Objx and Objpre schemes of scenario 3 (II).

Table 6 .
The index values of the minimum Obj x and Obj pre schemes of scenario 3 (I).

Table 7 .
The index values of the minimum Obj x and Obj pre schemes of scenario 3 (II).

Table 8 .
The detailed Pareto optimal schemes obtained of EMOBCGSA by different scenarios under pump outage conditions.

Table 8 .
The detailed Pareto optimal schemes obtained of EMOBCGSA by different scenarios under pump outage conditions.