Fractional-PID and Its Parameter Optimization for Pumped Storage Units Considering the Complicated Conduit System

: Speed governing control is signiﬁcant in ensuring the stable operation of pumped storage units. In this study, a state-space equation mathematical model of the pumped storage governing system considering the complex hydraulic pipeline structure of the pumped storage plant is proposed to describe the system’s dynamic behaviors under small disturbance conditions. Considering the frequent operating condition transitions and the complicated nonlinear dynamic characteristics of the pumped storage units, the fractional-order PID (FOPID) scheme that possesses a higher degree of control freedom than the traditional PID scheme is discussed in detail. To optimize the control parameters of the unit governor, an improved gravitational search algorithm (IGSA) that combines the basic searching mechanisms of the gravitational search algorithm and chaotic search, elastic sphere boundary treatment, and elite guidance strategy is developed. Comparative studies have been carried out under frequency and load disturbance conditions. Simulation results indicate that the control performance of FOPID is better than that of PID under diverse operating conditions and the proposed IGSA has satisfactory parameter optimization capability.


Introduction
With the gradual increase in the proportion of new energy sources such as wind power and solar energy in the electrical network [1,2], pumped storage is becoming more and more important as an effective regulation power source [3].As the core control system of pumped storage units, the pumped storage governing system (PSGS) undertakes the important tasks of manipulating the start-up and shutdown of the unit, the working condition conversions, and the peak regulation and frequency regulation [4].PSGS is a nonlinear complex system, which may present complicated dynamic behavior during operation.Therefore, the establishment of an accurate and effective simulation model [5] as well as the optimization of governor control parameters [6] can achieve a better control effect on the unit and ensure the stable operation of the unit.
With the development of pumped storage power plants, extensive studies have been reported on the modeling and control of PSGS.A typical PSGS consists of a hydraulic system, pump-turbine, generator, governor and electro-hydraulic servomechanism [7].In system modeling, the PSGS is similar to the conventional hydraulic turbine governing system (HTGS) in turbine mode.Usually, time-domain models are utilized to simulate the evolution of state variables during the transition process.The state-of-the-art methods include the method of characteristics (MOC) [8,9], the finite difference method [8,10], and the differential equation model [10], etc.Under small disturbance conditions, the system model ignores the system nonlinearities and selects a rigid water hammer model for the pipes and six-parameter model for the pump-turbine [11]; hence, the PSGS model can be described by linearized state space equations [12].For the speed governor, PID controller is the most commonly used for PSGS due to its structural simplicity and applicability.With the in-depth study of PID control laws, PID-like control laws such as nonlinear PID [13] and FOPID [14] are increasingly used.FOPID is widely used in various engineering fields [15], e.g., the chemical industry [16], nuclear power [17], and aerospace [18].Compared with traditional PID control, FOPID control possesses two more adjustable parameters, i.e., the differential order and the integral order.It allows the FOPID controller to have better robustness and control performance than the PID controller, but at the same time increases the difficulty of parameter optimization [19].
In pumped storage regulation systems, the governor is a very important component.The system's control performances are mainly decided by the tunable control parameters, and the optimization of these parameters can prominently enhance the system stability level.For the controller parameter optimization problem, the traditional rectification methods include the orthogonal test approach [20] and the simplex approach [21].Although these methods are simple to operate, they are not accurate enough to obtain the optimal parameters under complex operating circumstances.The meta-heuristic algorithms [22] (MA) have achieved better results in optimizing the controller parameters of hydropower units.MA include the Genetic Algorithm [23] (GA) based on Darwinian evolutionary theory, the Particle Swarm Optimization algorithm [24] (PSO) based on bird predation, the Differential Evolution Algorithm [25] (DEA) that uses the differences between random vectors to generate new vectors, the Gradient Descent Algorithm [26] (GDA) that uses the gradient search algorithm to keep approaching the optimal result, the Pattern Search Algorithm [27] (PSA) with axial exploration in the direction of feasible descent, and the Simulated Annealing Algorithm [28] (SA) based on the annealing process of solids.The Sine Cosine Algorithm (SCA) achieves close to the optimal result by the exact solution of the sine-cosine function [29], and the Gravity Search Algorithm [30] (GSA) based on gravity, etc.Unlike GA and PSO, which are based on biological phenomena, GSA is an optimization method based on the laws of gravity and mass interaction in physics.GSA has been verified to be more efficient than other optimization methods in optimizing the controller parameters [31], but still suffers from localized and premature convergence problems.However, almost all of the above research works have used the conventional PID controller which ultimately leads to a sub-optimal response.The superiority of the FOPID controllers are demonstrated in the results of the research work in this paper.
The accuracy and efficiency of optimization algorithms are crucial in solving engineering problems.In recent years, scholars have performed various research works in related fields and proposed many effective improvement ideas, including boundary processing, strategy improvement, algorithm fusion and parameter improvement.To improve the capability of GSA, Li et al. [32] proposed an improved gravitational search algorithm combining PSO and GSA.Sarafrazi et al. [33] proposed fusing GSA with a physics-inspired Kepler algorithm to speed up the algorithm's search and improve accuracy.Zhang et al. [34] designed and proposed a hybrid strategy based on Cauchy and Gaussian mutation to improve the exploration of GSA.Yin et al. [35] incorporate a cross-search in the GSA to improve the algorithm's development capability.He et al. [36] introduced the concept of repulsive force and proposed an improved gravitational search algorithm under the joint effect of repulsive and gravitational forces.Tian et al. [37] combined the Water Cycle Algorithm (WCA) with the GSA and use the concepts of watersheds and evapotranspiration to enhance the search capability.Based on previous research, this paper proposes an improved GSA (IGSA) in which the following four improvement strategies are incorporated into the GSA.Firstly, the chaotic operator is added to increase the diversity of the population and the randomness of the search; secondly, the adaptive gravitational decay factor is introduced to improve the change rule of the gravitational constant; subsequently, the elastic sphere strategy is carried out in the searching boundary treatment process; finally, the population of elite particles are applied to accelerate the convergence rate of the algorithm.
Water 2023, 15, 3851 3 of 23 When using MA to optimize controller parameters, defining a proper objective function is crucial for the optimization of results.Usually, the input of the objective function is the parameter to be optimized, and the output value is the fitness value.In this paper, we select an objective function that comprehensively takes account of the rotational speed error and the water pressure fluctuation during the transition process.
Most of the traditional models do not take into account the hydraulic structures of the pumped storage power plant, such as the surge tanks, and cannot reflect the transient processes of the crucial components.In order to improve the dynamic description ability of the model, this paper establishes an accurate state-space model of the PSGS which fully considers the hydraulic characteristics of the pipes, surge tanks and pump-turbine.In order to obtain better control performance, the FOPID controller with higher control degrees of freedom is chosen and compared to the conventional PID controller.In addition, the standard GSA algorithm was improved in this study to better optimize the control parameters.
The rest of this paper is organized as follows: In Section 2, the mathematical equations of the PSGS derived from the state-space equations are established.In Section 3, the fractional order calculus and the structure of the FOPID controller are introduced.Subsequently, Section 4 introduces the GSA algorithm and the proposed improvements.Then, the simulation results are analyzed in Section 5. Finally, the conclusions of this study are condensed in Section 6.

Mathematical Model
The hydroelectric power plant is an important component in the modern power system that converts water energy into electricity.Pumped storage power plant is a special kind of hydroelectric power plant, which uses water as an energy storage medium to store and manage electrical energy through the mutual conversion of electrical energy and potential energy.
The PSGS is composed of upstream and downstream reservoirs, penstock, surge tanks, pump-turbine, generator and speed governor.Water flows from the upstream reservoir through the diversion pipeline, upstream surge tank and penstock to the pump-turbine inlet, the kinetic energy of the running water drives the runner of the pump-turbine to rotate and generates mechanical energy, the shaft of the pump-turbine is connected to the generator and drives the generator rotor to rotate synchronously.The speed governor guarantees the system's frequency stability by regulating the guide vane opening.
This paper divides the PSGS into five parts: hydraulic system, pump-turbine, generator and load, speed governor, and hydraulic servomechanism, with combined modeling of the PSGS based on mathematical models of each component.

Hydraulic Systems' Modeling
A typical layout of a pumped storage plant with two surge tanks is shown in Figure 1.In order to study the control strategy for the speed governor, this work divides the PSGS into five parts: the hydraulic system, pump-turbine, generator and load, speed governor, and hydraulic servomechanism.
As shown in Figure 1, the water level of the upstream reservoir in the pumped storage power plant is H u (in meters), the cross-sectional areas of the diversion tunnel, penstock and downstream tailrace tunnel are A 1 ,A 2 ,A 3 ,A 4 (in square meters), the lengths are L 1 ,L 2 ,L 3 ,L 4 (in meters), the cross-sectional areas of the upstream and downstream surge tanks are A s1 ,A s2 , and the water levels are H s1 ,H s2 .Finally, a downstream reservoir with a water level of H d can be seen on the far right of the Figure 1, where the reservoir levels H u and H d are considered as constants.

Overall Pipeline Modeling
In applied hydraulic transients [38], if the pipe length is less than 600-800 m, or the governing system operates in small fluctuations conditions, the error of the rigid water hammer model can meet the engineering requirements.
If head loss is considered, the rigid water strike transfer function can be obtained as follows: where  is water inertia time constant; ℎ is the head loss of the tunnel.As depicted in Figure 1, the overall pipeline is divided into three parts: the upstream diversion tunnel, the penstock and the downstream tailwater tunnel.Using the rigid water strike transfer function, the mathematical model of the dynamic relationship between these three sections of the tunnel is given below.

Upstream diversion tunnel:
As displayed in Figure 1, the upstream diversion tunnel connects the upstream reservoir and the upstream surge tank.Since the upstream reservoir  is constant, the water pressure deviation remains constant.The end of the diversion tunnel is connected to the upstream surge tank, so the pressure deviation of the surge tank can be considered the same as the pressure deviation of the diversion tunnel.The mathematical model of the upstream diversion tunnel is derived from the rigid water strike Equation (1): where:  -water inertia time constant of the diversion tunnel; ℎ -relative head loss of the diversion tunnel; ℎ ,ℎ -relative water pressure deviation of the diversion tunnel and the upstream surge tank;

Overall Pipeline Modeling
In applied hydraulic transients [38], if the pipe length is less than 600-800 m, or the governing system operates in small fluctuations conditions, the error of the rigid water hammer model can meet the engineering requirements.
If head loss is considered, the rigid water strike transfer function can be obtained as follows: where T w is water inertia time constant; h f is the head loss of the tunnel.As depicted in Figure 1, the overall pipeline is divided into three parts: the upstream diversion tunnel, the penstock and the downstream tailwater tunnel.Using the rigid water strike transfer function, the mathematical model of the dynamic relationship between these three sections of the tunnel is given below.

1.
Upstream diversion tunnel: As displayed in Figure 1, the upstream diversion tunnel connects the upstream reservoir and the upstream surge tank.Since the upstream reservoir H u is constant, the water pressure deviation remains constant.The end of the diversion tunnel is connected to the upstream surge tank, so the pressure deviation of the surge tank can be considered the same as the pressure deviation of the diversion tunnel.The mathematical model of the upstream diversion tunnel is derived from the rigid water strike Equation (1): where: T w1 -water inertia time constant of the diversion tunnel; h f 1 -relative head loss of the diversion tunnel; h 1 ,h s1 -relative water pressure deviation of the diversion tunnel and the upstream surge tank; q 1 -the relative flow deviation of the diversion tunnel.

Penstock:
As presented in Figure 1, the penstock connects the upstream surge tank and pumpturbine inlet, and the tailrace tunnel connects the pump-turbine outlet and downstream Water 2023, 15, 3851 5 of 23 surge tank.The inertia of water flow is also considered in the penstock, and the mathematical model of the penstock is deduced as follows: where: T w2 ,T w3 -water inertia time constant of the penstock and the tailrace tunnel; h f 2 ,h f 3 -relative head loss of the penstock and the tailrace tunnel; h s2 ,h 2 ,h 3 -relative water pressure deviation of the upstream surge tank, the penstock and the tailrace tunnel; q t -relative flow deviation of the pump-turbine.

Downstream tailwater tunnel
As displayed in Figure 1, the tailrace tunnel connects the downstream surge tank and the downstream reservoir, the downstream reservoir H d is constant and the water pressure deviation is kept constant, the mathematical model of the downstream tailrace tunnel is deduced as follows: where T w4 -water inertia time constant of the downstream tailrace tunnel; h f 4 -relative head loss of the downstream tailrace tunnel; h s2 -relative water pressure deviation of the downstream tailrace tunnel; q 4 -relative flow deviation of the downstream tailrace tunnel.

Surge Tank Modeling
As shown in Figure 1, the upstream and downstream surge tanks are straight cylinder surge tanks.The surge tank plays the role of reflecting water hammer waves and reducing water hammer pressure.In omitting the surge tank inlet damping, the straight cylinder surge tank can be described as: where: ∆H-water level in the surge tank changes, and water level rises as positive; ∆Q-flow in and out of the surge tank, inflow is positive; A s -cross-sectional area of the surge tank.
Taking the relative value of Equation (5): where: T j -time constant of the surge tank; h s -relative water pressure deviation; q s -relative flow deviation.
Considering the flow continuity, the relationship between the flows in the upstream and downstream surge tanks can be deduced as follows [38]:    q 1 = q s1 + q 2 q 3 = q s2 + q 4 q 2 = q t = q 3 (7) where q 2 ,q 3 ,q s1 ,q s2 , respectively, are the relative flow deviations of the penstock, the tailrace tunnel, the upstream surge tank, and the downstream surge tank.

Pump-Turbine Modeling
For the study of the mathematical model of the pump-turbine, a six-parameter linear model can be used to describe the dynamic characteristics of the hydraulic turbine for small fluctuation problems, and the parameters of the pump-turbine are shown in Appendix A, and near a stable operating point, the linear model of the hydraulic turbine can be expressed as [11]: m t = e x x + e y y + e h h t q t = e qx x + e qy y + e qh h t (8) where: m t -relative value of hydraulic turbine torque deviation; h t -relative value of hydraulic turbine head deviation; x-relative value of rotation speed deviation; y-relative value of guide vane opening deviation; The six transfer functions are e x = ∂M t ∂x , e y = ∂M t ∂y , e h = ∂M t ∂h t , e qx = ∂Q ∂x , e qy = ∂Q ∂y , e qh = ∂Q ∂h t .

Generator Modeling
In the modeling of the governing system of a pumped storage power plant, the most frequently used model for the synchronous generator is the first-order model.That is, considering only the dynamic response process of rotor motion, and treating the generating motor as a rotating rigid shaft, which is coaxial with the pump-turbine and connected by a coupling, and the parameters of the generator are shown in Appendix B. The dynamic model considering the load is expressed as [11]: where: m g0 -load torque; ∂x -load self-regulation factor; T a -time constant of mechanical inertia of the generator, expressed as: where GD 2 is the flywheel torque of the rotating part of the unit; n r is the rated speed of the unit; P r is the rated output of the unit.The transition process of the pumped storage power plant governing system, the change in m g0 is usually regarded as a load disturbance to the PSGS.

Modeling of Electro-Hydraulic Servomechanism
The electro-hydraulic servomechanism is the actuator of the governor.It receives output signals from the controller and performs electro-hydraulic conversion and hydraulic amplification of the electrical control signals.Eventually, the electrical signal is converted into the displacement signal of the receiver.The electro-hydraulic servo mechanism is usually simplified to a first-order inertia link, as stated in Equation (10) [39], where: u-control signal output by the controller; y-stroke of the servo motor; T y -servomotor response time.

PSGS Model
Based on the separate mathematical models of each link of the governing system of the pumped storage power plant established earlier, a block diagram of the transfer function of the pumped storage power plant governing system can be obtained, as presented in Figure 2: As shown in Figure 2, the physical state variables , ,  , ℎ ,  , ℎ ,  , ℎ that can be measured, are chosen to describe the overall PSGS model in the form of the state space equations, as displayed in Equation (11), Based on the mathematical model of the penstock, the equation is discretized, and the flow values obtained using adjacent discrete sampling times are used to calculate Δ, then, ℎ ,ℎ can be expressed as: As shown in Figure 2, the physical state variables x, y, q 1 , h s1 , q t , h s2 , q 4 , h t that can be measured, are chosen to describe the overall PSGS model in the form of the state space equations, as displayed in Equation (11), e qx e n T a e qh x + e qy T y e qh − e qx e y T a e qh y +  (11) Based on the mathematical model of the penstock, the equation is discretized, and the flow values obtained using adjacent discrete sampling times are used to calculate ∆q, then, h 2 ,h 3 can be expressed as: Water 2023, 15, 3851 8 of 23

Definition of Fractional Calculus
Fractional calculus [40] is a mathematical theory that studies the properties of calculus and integral operators of arbitrary order, and it extends the order of calculus to the field of fractions and even plurals.Fractional calculus is a generalization of differentiation and integration of non-integer order, and its operation basic operation operator is Water 2023, 15, x FOR PEER REVIEW 8 of 23

Definition of Fractional Calculus
Fractional calculus [40] is a mathematical theory that studies the properties of calculus and integral operators of arbitrary order, and it extends the order of calculus to the field of fractions and even plurals.Fractional calculus is a generalization of differentiation and integration of non-integer order, and its operation basic operation operator is a t D α (, is the upper and lower bounds of the operation operator,  is the order of the calculus).
The calculus operation operator unifies the differentiation and integration, and a variety of definitions of fractional calculus emerge from the theory of fractional calculus, such as the ℎ definition, the   −  (  −  ) definition, the  − ( − ) and the  definition etc.
In this paper, using the  −  definition, the calculus of order  for a certain continuously derivable function () is: where .[ ] is the largest integer smaller than the real number , and ℎ is the calculation step.

FOPID Controller Structure
Compared with the traditional integer-order PID, FOPID controller introduces differential order  and integral order , which has two more degrees of freedom and can achieve better control effects.The PID and FOPID structure diagrams are given in Figure 3.The transfer function of the FOPID controller can be expressed as: Its time domain expression is: where () is the error signal, () is the controller output signal,  , , are the controller gain parameters, ,  are the fractional calculus orders, and ,  ∈ [0,2].When  =  = 1, the FOPID controller is transformed into the integer-order PID controller.(a, t is the upper and lower bounds of the operation operator, α is the order of the calculus).
The calculus operation operator unifies the differentiation and integration, and a variety of definitions of fractional calculus emerge from the theory of fractional calculus, such as the Cauchy definition, the Gr In this paper, using the G − L definition, the calculus of order α for a certain continuously derivable function f (t) is: lgorithm of the Speed Governor tional Calculus lus [40] is a mathematical theory that studies the properties of calcurators of arbitrary order, and it extends the order of calculus to the even plurals.Fractional calculus is a generalization of differentiation on-integer order, and its operation basic operation operator is a t D α d lower bounds of the operation operator,  is the order of the calcueration operator unifies the differentiation and integration, and a vaof fractional calculus emerge from the theory of fractional calculus, ℎ definition, the   −  (  −  ) definition, the ( − ) and the  definition etc. sing the  −  definition, the calculus of order  for a certain confunction () is: where is the largest integer smaller than the real number t−α h , and h is the calculation step.

FOPID Controller Structure
Compared with the traditional integer-order PID, FOPID controller introduces differential order µ and integral order λ, which has two more degrees of freedom and can achieve better control effects.The PID and FOPID structure diagrams are given in Figure 3.The transfer function of the FOPID controller can be expressed as: such as the ℎ definition, the   −  (  −  ) definition,  − ( − ) and the  definition etc.
In this paper, using the  −  definition, the calculus of order  for a certain c tinuously derivable function () is: .[ ] is the largest integer smaller than the number , and ℎ is the calculation step.

FOPID Controller Structure
Compared with the traditional integer-order PID, FOPID controller introduces ferential order  and integral order , which has two more degrees of freedom and achieve better control effects.The PID and FOPID structure diagrams are given in   Its time domain expression is: where e(t) is the error signal, u(t) is the controller output signal, K p ,K i ,K d are the controller gain parameters, λ, µ are the fractional calculus orders, and λ, µ ∈ [0, 2].When λ = µ = 1, the FOPID controller is transformed into the integer-order PID controller.[30].The law of gravity is a law that explains the relationship between objects interacting with each other.In the gravitational search algorithm, gravity is equivalent to an information transfer tool that enables information sharing among individuals and optimal search by the group under the effect of gravity.
Assuming that there are N particles, the position and velocity of the ith individual in the Ddimensional space are denoted as where x d i and v d i , respectively, denote the position and velocity components of individual i in the d-dimensions, and the position denotes the solution of the problem.The position and velocity are first initialized in the solution space and velocity space, the objective function values of each individual are calculated and evaluated, the mass, gravitational force and acceleration of each individual are calculated, and finally the velocity and position of the individual are updated [41].
According to Newton's gravitation theory, in the d dimensions, the gravitational force of individual j on i is expressed as: where G(t), R ij (t), respectively, are the universal gravitational constants at the tth iteration and the Euclidean distances of individuals i and j, with the following expressions: where G 0 , α and ε are constants, t is the current number of iterations and T is the maximum number of iterations, i,j {1, 2, • • • , N} and i = j.
In the d dimensions, the combined force on individual i is: where rand j is a random variable obeying a uniform distribution between [0, 1], and Kbest is the collection of the top k individuals with the optimal fitness value and maximum quality.
The updated velocity and position can be obtained as:

Improved Gravitational Search Algorithm (IGSA)
In order to improve the global search capability of GSA and avoid falling into a local optimum, the following four improvement strategies are incorporated in this paper: Chaos operator.Chaos [42] is a seemingly random movement that cannot be repeated precisely that occurs in nature.Chaotic motion is a complex state of motion unique to deterministic nonlinear dynamical systems, it has characteristics such as class randomness, initial value sensitivity, and ergodicity.
Chaos mapping is a mapping (evolutionary function) that exhibits some chaotic behavior, some common chaos mappings include Logistic map, Sinusoidal map, Tent map, Sine map and Bernoulli map.In this paper, we introduce the chaos operator to improve the GSA by choosing the ergodic logistic mapping, namely: where control parameters u ∈ [0, 4], r ∈ (0, 1) and r 0 / ∈ [0.25, 0.5, 0.75, 1], r t ∈ (0, 1) is the number of chaos generated in the tth iteration.
In the position update phase, chaos sequences are introduced to improve the convergence of the algorithm, while chaos perturbations can help individuals escape from their current position when they fall into a local optimum.The steps are as follows: (1) A d-dimensional random vector c d 1 ∈ [0, 1] is generated, and the control parameter u = 4, that is a fully chaos state is reached by Equation ( 20): , generating a chaos vector with k denoting the number of individuals; (2) The generated chaos vector is added to the search space, at which point the speed update equation in the improved algorithm is: where ξ is the factor that controls the range of chaos; (3) Using the improved velocity update equation to calculate the position at the next moment: Adaptive universal gravitational constant decay factor.
From the gravitational force calculation Equation ( 16), we can see that the universal gravitational constant G(t) is positively related to the gravitational force value, and its value plays an important role in the calculation of the gravitational force.When G(t) takes a larger value, the search range of the algorithm is wider to avoid falling into the local optimum, while when G(t) takes a smaller value, the algorithm can better converge to the global optimum.From Equation (17), it can be seen that the universal gravitational constant decay factor α in the standard GSA is taken as a constant, which limits the performance of the algorithm.Therefore, an adaptive universal gravitational constant decay factor is proposed: where α 0 is the initial value of the gravitational decay factor, w and δ are scaling factors, w ∈ [0, 1], δ > 1, θ is the shift factor, which is taken in this paper as α 0 = 20, w = 1, δ = 100, θ = 0.35.
As the number of iterations increases, the universal gravitational constant shows a nonlinear decrease, then the improved algorithm has a large gravitational constant at the beginning to enhance the search ability, and the algorithm reduces the gravitational constant at the later stage to accelerate the convergence of the algorithm.
Elastic ball boundary treatment.Boundary constraints are added when initializing individuals, but as the iteration proceeds, there may be individuals that cross the boundary after the position update, and for the individuals that cross the boundary, the elastic ball boundary is handled by: It can be seen from Equation ( 24) that the elastic ball boundary treatment increases the individual positions, enhances the search capability of the algorithm, and reduces the risk of the algorithm falling into a local optimum.
The velocity update formula of GSA only considers the effect of acceleration, while the velocity update in the particle swarm seeking algorithm takes into account the individual memory and group information exchange, it is more convergent for particle seeking.Therefore, the velocity update of IGSA combined with the chaos arithmetic and PSO velocity update mechanism is obtained: where c 1 , c 2 are the adaptive learning factors that vary with the number of iterations, F d best and X d best denote the individual optimal position and the global optimal position in the d-dimensions, respectively.
Because the addition of the optimal position in the velocity update equation will reduce the search capability of the algorithm, the adaptive learning factor is added to balance the global search capability and local optimization capability of the algorithm.c 1 , c 2 update equations [43] are: Based on the above improvement strategies, a new improved gravitational search algorithm (IGSA) is proposed in this paper with the following steps: Step 1: Random initialization of the population and setting parameters.The setup parameters include the population size N, the maximum number of iterations T, the boundary of the position [Lb, Ub], the initial value of the universal gravitational force G 0 and the initial value of the gravitational decay factor α 0 .The population individual positions are randomly initialized.
Step 2: It is determined whether the updated individual position is out of bounds.The transgression uses the elastic ball boundary strategy of Equation (24) to assign a new position.
Step 3: Calculate the fitness value of all individuals in the population, and determine whether the fitness value is NaN (Not A Number) or not, if it is, then randomly initialize the position of the individual and calculate the fitness value again, and repeat this step until the fitness value is an output table value.
Step 4: Update global optimal position X d best and individual optimal position F d best , if f zbest(t) < f zbest, then X d best = zbest(t), f zbest= f zbest(t).
Step 5: Calculating the individual mass M i , the gravitational parameter G is calculated according to Equations ( 17) and ( 23), the gravitational force F d i is calculated according to Equations ( 16) and (18), finally calculating the acceleration a d i .
Step 6: The adaptive learning factor is calculated according to Equations (26), and the speed and position are updated according to Equations ( 25) and (19).
Step 7: The number of iterations t = t + 1, if t < T, go to step 2, otherwise end the loop and obtain the optimal objective function value f zbest and the optimal position vector X d best .The algorithm flow chart of IGSA is indicated in Figure 4:

Objective Function
Time multiplied by the integral of the absolute value of the error (ITAE) is as an evaluation index for the performance of the speed regulation system of h units.The objective function proposed in this paper takes into account the rotat error value and the water pressure deviation.The calculated ratio of ITAE to of samples is taken as the first output, and the second output is the absolute va pressure deviation minimum value during the transition process.By adjusting of both during the parameter optimization process to ensure that the two ou are approximately equal.The objective function can be presented as:

Objective Function
Time multiplied by the integral of the absolute value of the error (ITAE) is often used as an evaluation index for the performance of the speed regulation system of hydropower units.The objective function proposed in this paper takes into account the rotational speed error value and the water pressure deviation.The calculated ratio of ITAE to the number of samples is taken as the first output, and the second output is the absolute value of water pressure deviation minimum value during the transition process.By adjusting the weights of both during the parameter optimization process to ensure that the two output values are approximately equal.The objective function can be presented as: Its discrete expression is given by: where w 1 , w 2 are weighting factors, N is the number of samples, T indicates the time sequence and h tmin indicates the minimum value of water pressure deviation.

Simulation Parameter Setting
In this section, the HTGS model proposed in Section 2.5 is simulated in the MATLAB R2020b software environment, the simulation environment is under a 2.4 GHz inteli7 CPU and 16G RAM.Tables 1 and 2 show the simulation parameters of HTGS and the transfer coefficients of the pump-turbine, respectively.The values of the boundary on the FOPID and PID tuning parameters are given in Table 3.
Based on the above unit parameters, the frequency disturbance and load disturbance simulation tests are performed in this study.

Load Disturbance Scenario
According to the model established in this paper, under a 10% load disturbance, the optimization performances of four parameter optimization algorithms (i.e., SA, GA, PSO and IGSA) and the control effects of the PID and FOPID are compared, respectively.

Comparison of SA, GA, PSO and IGSA
The parameters of SA are set as follows: initial temperature T s = 100, termination temperature T e = 0.01, temperature decay factor α = 0.95, maximum number of iterations T = 100, maximum number of iterations at temperature T k Lk = 30.
The parameters of GA are set as follows: population size N = 30, maximum number of iterations T = 100, crossover probability P c = 0.7, variance probability P m = 0.01.
The parameters of PSO are set as follows: population size N = 30, maximum number of iterations T = 100, learning factor c 1 = c 2 = 2, velocity weight w = 0.6.
The parameters of IGSA are set as follows: population size N = 30, maximum number of iterations T = 100, G 0 = 20, the gravitational constant decay factor parameter is set the same as revealed in Section 4.2.
When the weight factor w 1 = 0.05, w 2 = 0.95, after 100 iterations, the convergence curves of the fitness functions of SA, GA, PSO and IGSA are shown in Figure 5: According to the fitness convergence trend displayed in Figure 5, it can be seen th among the four algorithms, the final fitness value of GA is the largest and that of IGSA the smallest.From the convergence curve, SA starts with the largest fitness value and d creases the fastest, PSO decreases the fitness more slowly during the descent process, an GA converges the slowest.Comparing the final fitness values, it can be seen that SA, G and PSO easily fall into the local optimum and have difficulty in jumping out of it, whi the final fitness value of IGSA is smaller.This shows that as the number of iterations in creases, the fitness value of IGSA has been decreasing.The chaos perturbations added the position update phase enables IGSA to jump out of the local optimum continuousl which indicates that IGSA has better global exploration ability.
Table 4 shows the dynamic indicators for SA, GA, PSO and IGSA with optimal con trol parameters:  According to the fitness convergence trend displayed in Figure 5, it can be seen that among the four algorithms, the final fitness value of GA is the largest and that of IGSA is the smallest.From the convergence curve, SA starts with the largest fitness value and decreases the fastest, PSO decreases the fitness more slowly during the descent process, and GA converges the slowest.Comparing the final fitness values, it can be seen that SA, GA and PSO easily fall into the local optimum and have difficulty in jumping out of it, while the final fitness value of IGSA is smaller.This shows that as the number of iterations increases, the fitness value of IGSA has been decreasing.The chaos perturbations added in the position update phase enables IGSA to jump out of the local optimum continuously, which indicates that IGSA has better global exploration ability.
Table 4 shows the dynamic indicators for SA, GA, PSO and IGSA with optimal control parameters: Figure 6 compares the relative values of rotation speed deviation x, guide vane opening deviation y and hydraulic turbine head deviation h t variation curves under the optimal control parameters of SA, GA, PSO and IGSA during load disturbance: Combining Table 4 and Figure 6, the dynamic indexes under the optimal contro rameters of IGSA are significantly better than the other three algorithms.In Table 4, has the smallest objective function, stabilization time and absolute value of the mini value of water pressure deviation.It can also be seen in Figure 6 that the unit speed ation transition process optimized by the IGSA algorithm is significantly improved pared to the other three algorithms.The speed increase process is smoother, the stab tion time is shorter, and the absolute value of the minimum value of the water pre deviation is also the smallest, which achieves better control parameter optimization.O all, the transition process of  ,  , ℎ under the optimal control parameters of IG obviously improved.

Comparison of PID and FOPID
In order to prove that FOPID control has a better control effect than PID cont the parameter optimization of pumped storage speed regulation systems, tests sele optimal control parameters for both controls.Combining Table 4 and Figure 6, the dynamic indexes under the optimal control parameters of IGSA are significantly better than the other three algorithms.In Table 4, IGSA has the smallest objective function, stabilization time and absolute value of the minimum value of water pressure deviation.It can also be seen in Figure 6 that the unit speed deviation transition process optimized by the IGSA algorithm is significantly improved compared to the other three algorithms.The speed increase process is smoother, the stabilization time is shorter, and the absolute value of the minimum value of the water pressure deviation is also the smallest, which achieves better control parameter optimization.Overall, the transition process of x, y, h t under the optimal control parameters of IGSA is obviously improved.

Comparison of PID and FOPID
In order to prove that FOPID control has a better control effect than PID control in the parameter optimization of pumped storage speed regulation systems, tests select the optimal control parameters for both controls.Combining Table 4 and Figure 6, the dynamic indexes under the optimal cont rameters of IGSA are significantly better than the other three algorithms.In Table 4 has the smallest objective function, stabilization time and absolute value of the min value of water pressure deviation.It can also be seen in Figure 6 that the unit spee ation transition process optimized by the IGSA algorithm is significantly improve pared to the other three algorithms.The speed increase process is smoother, the sta tion time is shorter, and the absolute value of the minimum value of the water pr deviation is also the smallest, which achieves better control parameter optimization all, the transition process of  ,  , ℎ under the optimal control parameters of I obviously improved.

Comparison of PID and FOPID
In order to prove that FOPID control has a better control effect than PID con the parameter optimization of pumped storage speed regulation systems, tests sel optimal control parameters for both controls.Figure 7 compares ,  and ℎ va curves of the model under the load disturbance of the two controllers: Table 5 shows the dynamic indicators for PID and FOPID with optimal cont rameters: Table 5 shows the dynamic indicators for PID and FOPID with optimal control parameters: In Figure 7, it is obvious that the dynamic performance of the FOPID controller is significantly better than that of the PID controller during load disturbance.From the data presented in Table 5, it can be concluded that the values of the objective function, stabilization time and absolute value of the minimum value of water pressure deviation under the FOPID control are reduced by 0.0053, 19 and 0.0042, respectively, which verifies the superiority of the FOPID control.

Frequency Disturbance Scenario
Under a 5% frequency disturbance, this section compares four different parameter optimization algorithms of SA, GA, PSO and IGSA and the control effects of the two controllers PID and FOPID.In Figure 7, it is obvious that the dynamic performance of the FOPID controller is significantly better than that of the PID controller during load disturbance.From the data presented in Table 5, it can be concluded that the values of the objective function, stabilization time and absolute value of the minimum value of water pressure deviation under the FOPID control are reduced by 0.0053, 19 and 0.0042, respectively, which verifies the superiority of the FOPID control.

Frequency Disturbance Scenario
Under a 5% frequency disturbance, this section compares four different parameter optimization algorithms of SA, GA, PSO and IGSA and the control effects of the two controllers PID and FOPID.According to the fitness convergence trend displayed in Figure 8, in the 5% frequency disturbance, SA starts with the largest fitness value and the fitness value falls into the local optimum after two decreases, PSO falls into the local optimum after about seven iterations, and GA has the largest fitness value after the continuous optimization search.Comparing the final fitness values of the four algorithms, it is observed that SA, GA and PSO easily fall into the local optimum and have difficulty jumping out of it.The proposed IGSA has a wider search range and better search capability, and the optimal fitness value decreases continuously during the iterative process.It proves that the IGSA has better applicability and explorability.
Table 6 shows the dynamic indicators for SA, GA, PSO and IGSA with optimal control parameters:   9 compares ,  and ℎ variation curves under the optimal control parameters of SA, GA, PSO and IGSA during frequency disturbance: Combining Table 6 and Figure 9, the dynamic indexes under the optimal control parameters of IGSA are significantly better than the other three algorithms.In Table 6, the values of the objective function, overshoot, stabilization time and absolute value of the minimum value of water pressure deviation of IGSA are lower than the other three algorithms.From Figure 9, the speed deviation transition process of the unit optimized by the IGSA algorithm is significantly improved compared to the other three algorithms.The speed increase process is smoother and reaches stability faster, the speed overshoot is almost nothing, and the absolute value of the minimum value of the hydraulic pressure deviation is also the smallest, which achieves a better optimization of the control parameters.In conclusion, the transition process of , , ℎ under the optimal control parameters of IGSA is obviously improved.Combining Table 6 and Figure 9, the dynamic indexes under the optimal control parameters of IGSA are significantly better than the other three algorithms.In Table 6, the values of the objective function, overshoot, stabilization time and absolute value of the minimum value of water pressure deviation of IGSA are lower than the other three algorithms.From Figure 9, the speed deviation transition process of the unit optimized by the IGSA algorithm is significantly improved compared to the other three algorithms.The speed increase process is smoother and reaches stability faster, the speed overshoot is almost nothing, and the absolute value of the minimum value of the hydraulic pressure deviation is also the smallest, which achieves a better optimization of the control parameters.In conclusion, the transition process of x, y, h t under the optimal control parameters of IGSA is obviously improved.Table 7 shows the dynamic indicators for PID and FOPID with optimal c rameters: It is evident from Figure 10 that the FOPID controller achieves better resul control at frequency disturbance.By analyzing the results of each parameter in is found that the values of the objective function, overshoot, stabilization time lute value of the minimum value of water pressure deviation are reduced by 0.01 13.7 and 0.0255, respectively.Combined with Table 7 and Figure 10, it is found t frequency disturbance, the adjustment quality improvement under FOPID contr obvious.

Results Analysis
Through the frequency disturbance and load disturbance tests on the mo paring the four algorithms of SA, GA, PSO and IGSA, it is proved that the op results of IGSA are better than those of SA, GA and PSO, and IGSA has better se ity, while the FOPID parameters optimized by IGSA can make the pumped sto lation system have better dynamic performance; comparing the control effec FOPID and PID controllers, it can be seen that the dynamic performance obtaine the action of the FOPID controller is better than that of the PID controller.Table 7 shows the dynamic indicators for PID and FOPID with optimal control parameters: It is evident from Figure 10 that the FOPID controller achieves better results in PSGS control at frequency disturbance.By analyzing the results of each parameter in Table 6, it is found that the values of the objective function, overshoot, stabilization time and absolute value of the minimum value of water pressure deviation are reduced by 0.0198, 43.7%, 13.7 and 0.0255, respectively.Combined with Table 7 and Figure 10, it is found that under frequency disturbance, the adjustment quality improvement under FOPID control is more obvious.

Results Analysis
Through the frequency disturbance and load disturbance tests on the model, comparing the four algorithms of SA, GA, PSO and IGSA, it is proved that the optimization results of IGSA are better than those of SA, GA and PSO, and IGSA has better search ability, while the FOPID parameters optimized by IGSA can make the pumped storage regulation system have better dynamic performance; comparing the control effects of both FOPID and PID controllers, it can be seen that the dynamic performance obtained through the action of the FOPID controller is better than that of the PID controller.

Conclusions
This paper establishes a simulation model of a pumped storage speed regulation system with double surge tanks and compares the control effects of two controllers PID and FOPID.The IGSA is proposed to optimize the governor parameters.By comparing the convergence curves of the fitness function under load disturbance and frequency disturbance, it is verified that the IGSA proposed in this paper has better search capability than SA, GA and PSO, and it can address the algorithm precociousness and the local optimum trapping problem.In addition, under load disturbance and frequency disturbance, comparing the control effects under two different controllers PID and FOPID, it is found that FOPID has better comprehensive control performances.
In summary, the main contribution of this paper lies in the following two aspects: (1).A new PSGS state-space equation model which fully considers the hydraulic characteristics of the pipes, surge tanks and pump-turbine and the eletromechanical behaviors of the generator and hydraulic servomechanism is proposed.
(2).An improved GSA optimizer combining the basic searching mechanisms of the gravitational search algorithm and chaotic search, elastic sphere boundary treatment, and elite guidance strategy is developed for the control parameter optimization of the FOPID scheme.Through comparative case studies under load disturbance and frequency disturbance, it is proved that the proposed IGSA shows superiority over parameter optimizers on governor parameter optimization.
Future work will also focus on comparing simulation results with experimental results.relative flow deviation q 1 , q 2 , q 3 , q 4 , q s1 , q s2 relative flow deviation of diversion tunnel, penstock, tailrace tunnel, downstream tailrace tunnel, upstream surge tank and downstream surge tank ∆H water level in the surge tank changes ∆Q flow in and out of surge tank A s cross-sectional area of surge tank T j time constant of surge tank T j1 , T j2 time constant of upstream surge tank and downstream surge tank m t relative value of hydraulic turbine torque deviation h t relative value of hydraulic turbine head deviation q t , q t−1 relative flow deviation of the pump-turbine at moments t and t − 1 x, y relative value of rotate speed deviation and guide vane opening deviation y (servo mechanism) stroke of the servo motor e x , e y , e h partial derivatives of the torque with respect to head, guide vane and turbine speed e qx , e qy , e qh partial derivatives of the flow with respect to head, guide vane and turbine speed T a time constant of mechanical inertia of the generator m g0 load torque e g load self-regulation factor e n synthetic self-regulation coefficient GD 2  flywheel torque of rotating part of the unit n r rated speed of the unit P r rated output of the unit u control signal output by the controller T y servomotor response time FOPID:

EVIEW 4 of 23 Figure 1 .
Figure 1.Typical layout of a pumped storage plant with two surge tanks.

Figure 1 .
Figure 1.Typical layout of a pumped storage plant with two surge tanks.

Water 2023 , 23 Figure 2 .
Figure 2. Block diagram of the transfer function of the pumped storage power plant governing system.

Figure 2 .
Figure 2. Block diagram of the transfer function of the pumped storage power plant governing system.

Figure 3 .
Figure 3.Comparison of the PID and FOPID structure diagram.
is the largest integer smaller than the real is the calculation step.r Structure h the traditional integer-order PID, FOPID controller introduces difd integral order , which has two more degrees of freedom and can ol effects.The PID and FOPID structure diagrams are given in Figure ion of the FOPID controller can be expressed as: () = () () =  +   +   (14) expression is: () =  () +   () +   () (15) ror signal, () is the controller output signal,  , , are the conters, ,  are the fractional calculus orders, and ,  ∈ [0,2].When D controller is transformed into the integer-order PID controller. of the PID and FOPID structure diagram.

Figure 3 .
Figure 3.Comparison of the PID and FOPID structure diagram.

Figure 3 .
Figure 3.Comparison of the PID and FOPID structure diagram.

Water 2023 ,Figure 5 .
Figure 5. Convergence curves of the fitness functions under load disturbance.

Figure 5 .
Figure 5. Convergence curves of the fitness functions under load disturbance.

Figure 6 .
Figure 6.Variation curves of , , ℎ of four algorithms under load disturbance.

Figure 7 comparesFigure 6 .
Figure 6.Variation curves of n, y, h t of four algorithms under load disturbance.
Figure 7 compares x, y and h t variation curves of the model under the load disturbance of the two controllers:

Figure 6 .
Figure 6.Variation curves of , , ℎ of four algorithms under load disturbance.

Figure 7 .
Figure 7. Variation curves of , , ℎ of two controllers under load disturbance.

Figure 7 .
Figure 7. Variation curves of n, y, h t of two controllers under load disturbance.
5.3.1.Comparison of SA, GA, PSO and IGSAThe SA, GA, PSO and IGSA parameters are set the same as in Section 5.2.1.When the weight factor w 1 = 0.115, w 2 = 0.885, after 100 iterations, the convergence curves of the fitness functions of SA, GA, PSO and IGSA are displayed in Figure8: 5.3.1.Comparison of SA, GA, PSO and IGSAThe SA, GA, PSO and IGSA parameters are set the same as in Section 5.2.1.When the weight factor  = 0.115,  = 0.885, after 100 iterations, the convergence curves of the fitness functions of SA, GA, PSO and IGSA are displayed in Figure8:

Figure 8 .Figure 8 .
Figure 8. Convergence curves of the fitness functions under frequency disturbance.

Figure 9 .
Figure 9. Variation curves of , , ℎ of four algorithms under frequency disturbance.

Figure 9 .
Figure 9. Variation curves of n, y, h t of four algorithms under frequency disturbance.

Figure 10 Figure 10 .
Figure 10 compares x, y and h t variation curves of the model under the frequency disturbance of the two controllers:

Figure 10 .
Figure 10.Variation curves of n, y, h t of two controllers under frequency disturbance.

]
studies the properties of calcu-, and it extends the order of calculus to the calculus is a generalization of differentiation operation basic operation operator is a t D α ration operator,  is the order of the calcuhe differentiation and integration, and a vaerge from the theory of fractional calculus,  −  (  −  ) definition, the definition etc. , the calculus of order  for a certain conis the largest integer smaller than the real rder PID, FOPID controller introduces difhas two more degrees of freedom and can PID structure diagrams are given in Figure er can be expressed as:+   +  (14)  () +   () (15) ntroller output signal,  , , are the conal calculus orders, and ,  ∈ [0,2].When d into the integer-order PID controller.ture diagram.

Kbest
the collection of the top k individuals with the optimal fitness value and maximum quality a d i (t) the acceleration of individual i in the d-dimensions r t the number of chaos generated in the tth iteration u chaos control parameters c d 1 a d-dimensional random vector ξ factor that controls the range of chaos α 0 initial value of the gravitational decay factor w, δ scaling factors θ shift factor Ub(d), Lb(d) the upper and lower boundaries of the positions of the individuals in d-dimension c 1 , c 2 adaptive learning factors F d best , X d best individual optimal position and the global optimal position in the d-dimensions f fitness value w 1 ,w 2 weighting factors N, T the number of samples and time sequence h tmin the minimum value of water pressure deviation Inspired by the law of universal gravity, Rashedi et al. proposed a new swarm intelligence optimization algorithm-the gravitational search algorithm (GSA) in 2009

Table 2 .
5, 0.6 Head loss of the tunnel h f 1 , h f 2 , h f 3 , h f 4 0.0026, 0.003, 0.003, 0.0015 Time constant of the surge tan k T j1 , T j2 Transfer coefficients of water pump-turbine.

Table 3 .
The values of the boundary on the FOPID and PID tuning parameters.

Table 4 .
Dynamic indicators of SA, GA, PSO and IGSA optimal control parameters under load di turbance.
Figure 6 compares the relative values of rotation speed deviation , guide vane open ing deviation  and hydraulic turbine head deviation ℎ variation curves under the op timal control parameters of SA, GA, PSO and IGSA during load disturbance:

Table 4 .
Dynamic indicators of SA, GA, PSO and IGSA optimal control parameters under load disturbance.

Table 5 .
Dynamic indicators of PID and FOPID optimal control parameters under load disturbance.

Table 5 .
Dynamic indicators of PID and FOPID optimal control parameters under load disturbance.

Table 6 .
Dynamic indicators of SA, GA, PSO and IGSA optimal control parameters under frequency disturbance.Figure9 comparesx, y and h t variation curves under the optimal control parameters of SA, GA, PSO and IGSA during frequency disturbance:

Table 6 .
Dynamic indicators of SA, GA, PSO and IGSA optimal control parameters under frequency disturbance.

Table 7 .
Dynamic indicators of PID and FOPID optimal control parameters under freq turbance.

Table 7 .
Dynamic indicators of PID and FOPID optimal control parameters under frequency disturbance.
Author Contributions: Conceptualization, Y.Z.(Yang Zheng)and X.Z.; writing-original draft preparation, X.Z.; methodology, B.X.; writing-review and editing, W.L. and Y.Z.(Yidong Zou); visualization, J.C.All authors have read and agreed to the published version of the manuscript.This work was supported by Hubei Provincial Nature Science Foundation of China [Grant No. 2022CFD165] and the Fundamental Research Funds for the Central Universities (Grant No. 2042022kf1022).For the sake of information security, the original data used in this paper will not be disclosed.Conflicts ofInterest: B.X. was employed by the company Technology and Research Center, China Yangtze Power Co., Ltd.The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.T w2 , T w3 , T w4 water inertia time constant of diversion tunnel, penstock, tailrace tunnel and downstream tailrace tunnel h s relative water pressure deviation h s1 , h s2 relative water pressure deviation of upstream surge tank and downstream tailrace tunnel h 1 , h 2 , h 3 relative water pressure deviation of the diversion tunnel, penstock, tailrace tunnel h f 1 , h f 2 , h f 3 , h f 4 relative head loss of diversion tunnel, penstock, tailrace tunnel and downstream tailrace tunnel q s Funding: