Characteristics Analysis and Fuzzy Fractional-Order PID Parameter Optimization for Primary Frequency Modulation of a Pumped Storage Unit Based on a Multi-Objective Gravitational Search Algorithm

: Compared with conventional hydropower units, the pumped storage unit has the characteristics of diverse working conditions and frequent switching. Therefore, the stability and regulation quality of the primary frequency modulation transition process of the regulating system is very important. Due to the “S” characteristics of the pumped storage unit (PSU), the pumped storage unit regulating system has a strong nonlinearity, and the conventional proportional-integral-derivative (PID) controller cannot provide high-quality control under low water head conditions. In this paper, the nonlinear PSU model with an elastic water hammer e ﬀ ect is studied, and the fuzzy fractional-order PID (FFOPID) controller is designed to improve the stability of the system. The membership function and the control parameters of the fractional-order PID are optimized based on the multi-objective gravitational search algorithm (MOGSA). The experimental results show that the optimized design of the FFOPID controller has better control quality than the traditional PID controller, the fractional-order PID (FOPID) controller, and the fuzzy PID controller (FPID) when the system is disturbed by the rotating speed under low water head.


Introduction
With the constant cost of traditional energy sources and the increasing pollution of the environment, the use of renewable energy has received extensive attention. Due to the simple adjustment process, fast adjustment speed, and economy of the pumped storage unit (PSU), it mainly undertakes the task of peaking and frequency modulation in the power system [1,2]. In recent years, with the formation of regional power grids, equipment fault diagnosis and the primary frequency modulation of the pumped storage unit regulating system has played a vital role in the safe and stable operation of the power system [3][4][5][6]. In previous primary frequency modulation studies, the modeling of PSU usually adopts the simplified linear model and the water diversion system model is based on the rigid water hammer theory [7,8]. However, the model of PSU is a complex nonlinear system [9]. In addition, the impact of elastic water hammer on pipeline must be considered for the pumped storage unit regulating system with a long water diversion system [10]. Therefore, it is necessary to establish a

The Model of the Water Diversion System
The elastic water hammer in the pressure pipeline considers the compressibility of water and the elasticity of the pipe. The basic equation of the water hammer consists of two equations: The motion equation and the continuous equation in the water hammer process [32,33]. The basic equation of the water hammer is the mathematical expression of the law of constant flow with pressure, which reflects the change law of flow velocity and water head in the hydraulic transition process. The differential segment with length dx taken out of the pipeline system is shown in Figure   2. The differential segment dx is subject to the friction force AD F against the fluid from the pipe wall on surface AD , the pressure AB F on surface AB , the pressure BC F on the fluid from the pipe wall on surface BC , the pressure CD F on surface CD , and the gravity G . Before deriving Equation (1), the following assumptions must be made: The flow in the pipeline system is onedimensional flow; the flow in each section of the pipeline system is assumed to be gradual flow; the resistance loss formula for calculating the constant flow in the pipeline can also be used to calculate the unsteady flow. According to Newton's second law, the force equation of the differential segment dx is as follows:

The Model of the Water Diversion System
The elastic water hammer in the pressure pipeline considers the compressibility of water and the elasticity of the pipe. The basic equation of the water hammer consists of two equations: The motion equation and the continuous equation in the water hammer process [32,33]. The basic equation of the water hammer is the mathematical expression of the law of constant flow with pressure, which reflects the change law of flow velocity and water head in the hydraulic transition process. The differential segment with length dx taken out of the pipeline system is shown in Figure 2.

The Model of the Water Diversion System
The elastic water hammer in the pressure pipeline considers the compressibility of water and the elasticity of the pipe. The basic equation of the water hammer consists of two equations: The motion equation and the continuous equation in the water hammer process [32,33]. The basic equation of the water hammer is the mathematical expression of the law of constant flow with pressure, which reflects the change law of flow velocity and water head in the hydraulic transition process. The differential segment with length dx taken out of the pipeline system is shown in Figure   2. The differential segment dx is subject to the friction force AD F against the fluid from the pipe wall on surface AD , the pressure AB F on surface AB , the pressure BC F on the fluid from the pipe wall on surface BC , the pressure CD F on surface CD , and the gravity G . Before deriving Equation (1), the following assumptions must be made: The flow in the pipeline system is onedimensional flow; the flow in each section of the pipeline system is assumed to be gradual flow; the resistance loss formula for calculating the constant flow in the pipeline can also be used to calculate the unsteady flow. According to Newton's second law, the force equation of the differential segment dx is as follows: where P is the pressure at the center of the AB section; A is the flow area of the pipeline; 0  is the shear stress of the pipe wall acting on the fluid periphery; D is pipe diameter;  is the fluid bulk density;  is the angle between the pipe centerline and the horizontal line; and V is the average velocity. The differential segment dx is subject to the friction force F AD against the fluid from the pipe wall on surface AD, the pressure F AB on surface AB, the pressure F BC on the fluid from the pipe wall on surface BC, the pressure F CD on surface CD, and the gravity G. Before deriving Equation (1), the following assumptions must be made: The flow in the pipeline system is one-dimensional flow; the flow in each section of the pipeline system is assumed to be gradual flow; the resistance loss formula for calculating the constant flow in the pipeline can also be used to calculate the unsteady flow. According to Newton's second law, the force equation of the differential segment dx is as follows: where P is the pressure at the center of the AB section; A is the flow area of the pipeline; τ 0 is the shear stress of the pipe wall acting on the fluid periphery; D is pipe diameter; γ is the fluid bulk density; α is the angle between the pipe centerline and the horizontal line; and V is the average velocity.
Since the shear stress τ 0 of a non-constant flow is a very complicated and unresolved problem, we have to borrow the research results of the constant flow to approximate the shear stress in the unsteady flow. The shear stress direction is always opposite to the flow velocity. In the actual work, we use the water head H as the main parameter, so we change P to H. In the water diversion system Energies 2020, 13, 137 4 of 20 of a hydropower station, the rate of change of the flow rate Q with space is small, so it is neglected. Moreover, the rate of change of the density p and the flow area of the pipeline A along x is much smaller than the rate of change of the flow rate V, so it is ignored. According to Equation (1), the basic equations describing water hammer can be obtained as the motion equation and the continuous equation: where Q is the flow rate of a certain section at time t; |Q| is the absolute value of the flow, which is positive or negative depending on the direction of the flow; H is the water head at time t of the section at a corresponding datum; g is the acceleration of gravity; a is the propagation velocity of water hammer wave; and f is the pipeline friction coefficient. By integrating Equations (2) and (3) along the characteristic line C + and C − , as shown in Figure 3, it can be converted into the following characteristic equation [33]: where The C a and R are constants; ∆t for time step, it is necessary to meet the Courant stability condition ∆t < ∆x (a+V) . For the internal nodes of the pipeline, the value of flow Q P can be solved: The value of H P can be obtained by Equation (4) or Equation (5). In this way, the pressure and flow rate of each node in the pipeline can be obtained at all times.
we have to borrow the research results of the constant flow to approximate the shear stress in the unsteady flow. The shear stress direction is always opposite to the flow velocity. In the actual work, we use the water head H as the main parameter, so we change P to H . In the water diversion system of a hydropower station, the rate of change of the flow rate Q with space is small, so it is neglected. Moreover, the rate of change of the density p and the flow area of the pipeline A along x is much smaller than the rate of change of the flow rate V , so it is ignored. According to Equation (1), the basic equations describing water hammer can be obtained as the motion equation and the continuous equation: where Q is the flow rate of a certain section at time t ; Q is the absolute value of the flow, which is positive or negative depending on the direction of the flow; H is the water head at time t of the section at a corresponding datum; g is the acceleration of gravity; a is the propagation velocity of water hammer wave; and f is the pipeline friction coefficient.
By integrating Equations (2) and (3) along the characteristic line + C and -C , as shown in Figure 3, it can be converted into the following characteristic equation [33]: . The a C and R are constants; t  for time step, it is necessary to meet the Courant stability condition For the internal nodes of the pipeline, the value of flow P Q can be solved: The value of P H can be obtained by Equation (4) or Equation (5). In this way, the pressure and flow rate of each node in the pipeline can be obtained at all times.

FOPID Controller
The transfer function of conventional integer solution PID controller can be described as: where U(s) is the output of the controller, and E(s) is the error input of the controller. The traditional integer-order PID controller adjusts the system performance by K p , K i , and K d . The fractional-order PI λ D µ controller has the same three parameters as the PID controller: K p , K i , and K d . Two new parameters are added: Integral order λ and differential order µ. The conventional PID controller is just fixed value on the coordinate axis, while the fractional order PI λ D µ controller is fixed value on the entire coordinate plane. The transfer function of FOPID controller can be described as [34]: Depending on the order of the control object, the λ and µ ranges are adjusted. λ mainly affects the steady-state accuracy of the system, while µ mainly affects the overshoot of the system. With the addition of two degrees of freedom, the FOPID controller is more flexible and has a wider range of control. The most common method to realize fractional-order transfer functions in simulation or practical research is to approach them by integer-order transfer functions [31]. In [35], Oustaloup proposed an approximate method using the recursive distribution of N poles and N zeros. To achieve a balance between the complexity and accuracy of the FFOPID controller, a 5th order Oustaloup's recursive approximation is performed for the integral-differential operator in the selected frequency band of ω ∈ 10 −2 , 10 2 rad/s. The high-order filter is as follows: where the poles ω k , zeros ω k , and gain K of the filter are obtained from the following formula: where α is the order of differ-integration and α ∈ (0, 1), (2N + 1) is the order of filter. The Laplace transform α-th derivative with α ∈ R + of a signal x(t) relaxed at t = 0 is as follows:

FFOPID Controller
Based on fuzzy set theory and fuzzy reasoning, fuzzy control transforms the knowledge and control experience expressed by expert language into fuzzy rules, which are controlled by computer. When the rotational speed of PSU deviates from the rated value, the FFOPID controller in the system operation process continuously calculates e (deviation) and e c (rate of deviation change), according to the fuzzy rule principle online adjustment K p , K i , K d , so that the system can reach a stable state. The basic structure of FFOPID controller is shown in Figure 4.

Fuzzy inference
The knowledge base  The modules MF, R, and fd in the "knowledge base" in Figure 4 are, respectively, the membership function, the fuzzy rule, and the sharpening algorithm, which are obtained offline. The core modules of the fuzzy controller: Fuzzy module (D/F), approximate reasoning module (A * o R), and sharpening module (F/D). The D/F module completes the operation of converting the clear quantity into the fuzzy quantity. In this experiment, the fuzzy set of the input quantity and the output quantity is [NB, NM, NS, ZO, PS, PM, PB], respectively, represent negative big, negative medium, negative small, zero, positive small, positive medium, and positive big. The membership function of the input quantity is shown in Figure 5a,b. The membership function of the output quantity is represented by a Gaussian function, because the membership function of this kind is smooth, non-zero at all points, which may be more suitable to describe a complicated fuzzy relationship. The position parameters {cen p , cen i , cen d } and scale parameters { → wid p , → wid i , → wid d } of Gaussian function are set as optimization variables. In this paper, the location parameters and the scale parameters of the Gaussian function are optimized. The A * o R module completes the approximate inference operation of the input fuzzy quantity through fuzzy rules and obtains the output fuzzy quantity U. The settings of the fuzzy rules are shown in Tables 1-3. The F/U module completes the conversion of fuzzy quantity U to clear quantity. The relationship is as follows: In the formula, K p0 , K i0 , and K d0 are the initial parameters. e is the deviation between the actual speed and the rated speed of the PSU. Because the FOPID controller has {K p0 , K i0 , K d0 , λ, µ} five parameters that need to be designed, coupled with the introduction of the fractional-order, makes the design more difficult. So, this paper uses the MOGSA algorithm to set the initial parameters of the FOPID controller.

The Model of Servo Mechanism
According to the output u controlled by the governor, the servo mechanism of the governor operates the guide vanes of the PSU to control the flow of the unit. The servo mechanism is composed of assistant servomotor and main servomotor. Due to the existence of a dead zone module and

The Model of Servo Mechanism
According to the output u controlled by the governor, the servo mechanism of the governor operates the guide vanes of the PSU to control the flow of the unit. The servo mechanism is composed of assistant servomotor and main servomotor. Due to the existence of a dead zone module and saturated zone module, it is a typical nonlinear link. The structure is shown in Figure 6.
where u is the output signal of the controller; k 0 is the amplification factor; T yb is the time constant of the assistant servomotor; T y is the main servomotor time constant; y 1 is the main pressure valve opening; and y is the guide vane of the PSU.
where u is the output signal of the controller;

Model of the Pumped Storage Unit
There is a nonlinear relationship between the relative value of the change of the guide vane opening degree of the PSU and the relative value of the displacement of the servomotor. At present, no specific mathematical model can be established. Therefore, the simulation calculation of the PSU should be based on the full characteristic curve [36]. The full characteristic curve of the PSU is based on the vane opening  as an intermediate variable to establish the flow characteristic curve 11 11 QN  and torque characteristic curve 11 11 M N  , as shown in Figure 7. The different combinations of these operating parameters divide the full characteristic curve into five operating conditions, including: The pump zone, the pump brake zone, the turbine zone, the turbine brake zone, and the reverse pump zone.
When the frequency is adjusted in the working condition of the turbine, especially in the case of external load reduction under low water head, the unit speed of the turbine is shown as follows: where 11 Q is the unit flow, 11 M is the unit torque, 11 n is the unit rotation speed, n is the rotation speed, and D is the nominal diameter of the PSU.
It can be seen from Figure 7 that the flow characteristic curve and torque characteristic curve of the PSU are overlapped and twisted in the area, with a relatively large rotational speed, which is the "S" characteristic area. The nominal diameter D of the PSU is constant, and the reduction of external

Model of the Pumped Storage Unit
There is a nonlinear relationship between the relative value of the change of the guide vane opening degree of the PSU and the relative value of the displacement of the servomotor. At present, no specific mathematical model can be established. Therefore, the simulation calculation of the PSU should be based on the full characteristic curve [36]. The full characteristic curve of the PSU is based on the vane opening θ as an intermediate variable to establish the flow characteristic curve Q 11 − N 11 and torque characteristic curve M 11 − N 11 , as shown in Figure 7. The different combinations of these operating parameters divide the full characteristic curve into five operating conditions, including: The pump zone, the pump brake zone, the turbine zone, the turbine brake zone, and the reverse pump zone.
When the frequency is adjusted in the working condition of the turbine, especially in the case of external load reduction under low water head, the unit speed of the turbine is shown as follows: where Q 11 is the unit flow, M 11 is the unit torque, n 11 is the unit rotation speed, n is the rotation speed, and D is the nominal diameter of the PSU. It can be seen from Figure 7 that the flow characteristic curve and torque characteristic curve of the PSU are overlapped and twisted in the area, with a relatively large rotational speed, which is the "S" characteristic area. The nominal diameter D of the PSU is constant, and the reduction of external load leads to the increase of the unit rotation speed n. Then, the unit rotation speed n 11 with lower water head H is relatively large, which makes the PSU easily run into the reverse pump zone under the influence of "S" characteristics, thus causing the rotation speed n to fluctuate around the network frequency value.
In order to solve the single input-multiple out problem in "S" area, the Suter transformation method was adopted to convert the full characteristic curve of the PSU in the four quadrants into two inquire transformation curves. The transformed whole characteristic curve was transformed into WM(x, y) and WH(x, y) curves on the x axis, but there were still some problems such as uneven Energies 2020, 13, 137 9 of 20 opening distribution [37]. In view of the above problems, the improved Suter transformation was adopted [38]. The transformation relations are as follows.
where a = n n r ; ; C y 1 is the coefficient, 0.1-0.3; C y 2 is the coefficient, 0.4-0.6; and the subscript r represents the rating. According to the improved Suter transformation, the calculation model of pump turbine is expressed as follows: In Equation (14), h 11 and m 11 represent the interpolation of the characteristic curve of the PSU. The improved Suter transformation can effectively solve the obstacle of single input-multiple output in the "S" area. In order to solve the single input-multiple out problem in "S" area, the Suter transformation method was adopted to convert the full characteristic curve of the PSU in the four quadrants into two inquire transformation curves. The transformed whole characteristic curve was transformed into WM( x, y ) and WH( x, y ) curves on the x axis, but there were still some problems such as uneven opening distribution [37]. In view of the above problems, the improved Suter transformation was adopted [38]. The transformation relations are as follows.

The Model of Generator and Loads Unit
The dynamic equilibrium relationship between the PSU's active torque m t and the generator's load torque m g is described by the motion equation [39]. According to the dynamic principle, the motion equation can be established as follows: where T a is inertia time of the generator, and e g is the self-regulating coefficient of the generator. The gravitational search algorithm (GSA) optimization algorithm mainly uses the gravity law between two objects to guide the motion of each particle to optimize the search for the optimal solution, the position of each particle corresponds to a solution to the problem. In the whole search space, the particle with the largest inertia mass represents the optimal solution to the problem. The optimization process of GSA is as follows.
Initialize each parameter and population position, suppose there is a particle in a separate system. Define the position of the particle as follows: where x d i is the position of the ith mass in dth dimension, and N is the size of the population. The inertial mass is calculated by the fitness value of the particle: where f it i (t) represent the fitness value of the agent i at the time t, the best and worst fitness values are defined below (for a minimization problem): Calculate the joint force of gravitation on the dth dimension space of ith agent at the time t: where F d ij (t) is the magnitude of the gravitational force on particle i from particle j at t moment.
where the G(t) is the gravitational constant at t time, R ij is the distance between the ith agent and dth agent, ε is a constant which is set for avoiding the divisor equal to zero.
According to the law of motion, the acceleration a d i (t) of agent i at t time, and in direction dth is given as follows: Then, calculate the velocity and position value of the agent i at the next moment: Energies 2020, 13, 137 11 of 20

Multi-Objective Optimization
Different from the single objective particle swarm optimization, the multi-objective particle swarm optimization has no single optimal solution, but a group of non-dominant solutions. The flowchart of the MOGSA algorithm is shown in Figure 8 and the detailed process of MOGSAs algorithm is as follows: e. tep 6: For each individual population, update velocity and next position, respectively. tep 7: Evaluate the fitness value of each agent's new position. After, all of the new individ posed to uniform mutation operator and those individuals among the resulting populat is considered as non-dominated ones, are added to the archive. Finally, the keep nality removes those solutions that have recently become dominated in the archive. tep 8: Update the archive and check the number of non-dominant solutions. If the numbe ts in the archive gets to its limit, one of elements should be omitted from the archive.  Step 1: Set the maximum number of iterations N, population size nPop, Pareto Front size nRep, dimension of search space D, and search space upper S max and lower limits S min .
Step 2: The population positions randomly initialized in the search space and then the velocity vector corresponding to each particle is initialized to zero.
Step 3: Evaluate the fitness value of each agent on each object fitness function; next, update the archive which has a grid structure. The grid structure is created as follows: Each dimension in the objective space is divided into equal divisions, where the agent i denotes dimension index, and hence, for a k-objective optimization problem, there will be different segments.
Step 4: Update the gravitational constant the acceleration for each agent, calculate the mass (active mass) corresponding to each solution in the archive, and assign a uniform mass (passive mass) to the search particles.
Step 5: By calculating the agent mass, the optimal non-dominant solution is selected from the archive.
Step 6: For each individual population, update velocity and next position, respectively.
Step 7: Evaluate the fitness value of each agent's new position. After, all of the new individuals are exposed to uniform mutation operator and those individuals among the resulting population, which is considered as non-dominated ones, are added to the archive. Finally, the keeping functionality removes those solutions that have recently become dominated in the archive.
Step 8: Update the archive and check the number of non-dominant solutions. If the number of elements in the archive gets to its limit, one of elements should be omitted from the archive.
Step 9: Determine whether the number of iterations N is reached; if the conditions are not met, loop operation Step 3-Step 9. Otherwise, output the non-dominant solution stored in the archive.

Objective Functions
In order to obtain controller optimization parameters that can satisfy the comprehensive index of the unit speed and the vane opening degree, two objective functions are set. The first objective function considers the integral time absolute error (ITAE) of rotational speed, oscillation level of the unit speed (OUS), and steady-state error (SSE) of the unit. The second objective function considers the ITAE of guide vane opening and the degree of oscillation of the vane opening (OVO). The integral of time multiplied absolute error ITAE criterion is generally used to minimize time multiplied absolute error of the control system [40].
where OUS = Max(x) − Min(x); SSE = x(N) − c; OVO = Max(y) − Min(y). x(N) is the rotation speed of the unit after stabilization, c is the rated speed. Max(x) and Min(x), respectively, represent the maximum and minimum values of the rotation speed in the adjustment process. Max(y) and Min(y), respectively, represent the maximum and minimum values of guide blade opening in the regulation process. w 1 , w 2 , and w 3 are weight coefficients to set the value of the different indexes within the same preset range. The importance of the three control indexes ITAE x , OUS, and SSE in J 1 and the two control indexes ITAE y and OVO in J 2 have the same importance for the objective function value. Therefore, when selecting the weight coefficients of w 1 ,w 2 , and w 3 , matching with the magnitudes of ITAE x and ITAE y should be considered, so that the actual values of the five control indicators are of the same magnitude.

Optimization Variables
In order to obtain good stability in the regulation system, the controller parameters of the governor are usually optimized by the optimization algorithm. Therefore, this paper adopts the MOGSA algorithm to, respectively, optimize the parameters of PID control, FOPID controller, and FFOPID controller. The optimization variables of the different controllers are listed in Table 4. Table 4. Optimization variables of different controllers.

Controller Optimization Variables
PID controller Definition: PID is proportional-integral-derivative Definition: FOPID is fractional-order proportional-integral-derivative Definition: FPID is fuzzy proportional-integral-derivative Definition: FFOPID is fuzzy fractional-order proportional-integral-derivative

Experiences and Results Analysis
In this section, the nonlinear model of the pumped storage unit regulating system with elastic water hammer effect was simulated in MATLAB. In order to show that the FFOPID controller has a better control effect than the traditional controller, they were successively put into the multiple simulation experiments for comparison.

Model Parameters
In the following simulation experiments, parameters of MOGSA were set as: N = 100, nPop = 50, nRep = 30, G0 = 5, nGrid = 10. Optimization variables were searched within specified ranges, the adopted search range for the PID controller parameters were restricted to K p ,  Table 5.

Comparison of the Performance of Different Controllers
In this part, the simulation experiments show that the pumped storage unit regulating system is assumed to be disturbed by frequency under the 25% guide vane opening condition at low water head. Furthermore, experiments were conducted under 526 m water head condition. The simulation time was 100 s, and the relative values of 0.01 and −0.01 speed disturbance were added at 5 s. Four sets of comparative experiments were carried out in the above conditions to compare the performance of the different controllers, and optimization parameters of the group that can obtain the best speed response curve were selected from the Pareto optimal solution set, respectively, as illustrated in Figure 9. The specific results of the optimization parameters are shown in Table 6, and Figure 10 shows the membership function of the output variable optimized by the MOGSA algorithm.
The comparison of speed curves is shown in Figure 11. It can be seen from Figure 11 that when the speed disturbance occurs, the response curve with the FFOPID controller fluctuates less and the adjustment time is shorter. In order to show the optimization effect of the FFOPID controller, the experimental results of control indices are shown in Table 7. By comparing the corresponding control indices of the PID controller, FOPID controller, FPID controller, and FFOPID controller, the performance of these controllers can be fully displayed. It can be seen from Table 7 that when the PSU is disturbed by the relative speeds of 0.01 and −0.01 under the FFOPID controller, the control indexes ITAE x are 42.54 and 28.45, respectively, which are smaller than other controllers. The smaller ITAE x indicates that the speed curve of the PSU has better dynamic characteristics. From Figure 11, it can be seen that when the PSU is disturbed in 5 s, the speed curve enters into the oscillation area. Under the action of the FFOPID controller, the speed fluctuation is small, and the control index OUS value is lower than that of the traditional controller. According to the regulation time, the stability time values of the two disturbances are 8.66 s and 8.32 s, respectively, which are far less than the PID, FOPID, and FPID controllers. These indices show that the FFOPID controller is more suitable for improving the control performance of the nonlinear pumped storage unit regulating system with elastic water hammer effect, compared with other controllers. The transient processes of the pumped storage unit regulating with different controllers corresponding to rotational speed disturbance of 0.01 are shown in Figure 12, in which mechanical torque, guide vane opening, and characteristic curves of the water hammer of the surge chamber can be seen.     (7.7436,9.7880,9.  The transient processes of the pumped storage unit regulating with different controllers corresponding to rotational speed disturbance of 0.01 are shown in Figure 12, in which mechanical torque, guide vane opening, and characteristic curves of the water hammer of the surge chamber can be seen.
As shown in Figure 12a, the FFOPID controller can control the mechanical torque to a stable state faster than the PID controller and FOPID controller. It can also be seen from Figure 12b that under the action of the FFOPID controller, the guide blade opening can be balanced quickly and rapidly, so the corresponding rotational speed curve first reaches a stable state. Therefore, it can be considered that the FFOPID controller is better. Meanwhile, the results show that the oscillation amplitude of the water hammer in the surge chamber obtained by the FFOPID controller is smaller than that obtained by other controllers, as shown in Figure 12c,d. The above conclusion indicates that the FFOPID controller can provide more stable operating conditions.

Robustness Analysis of the Proposed Method
The low-frequency disturbance is the key factor to the stability of the pumped storage unit regulating system under the condition of low water head. Furthermore, the water hammer effect caused by the controller's poor robustness is an important cause of PSU vibration. Therefore, in this paper, the robustness test of the FFOPID controller was carried out by adjusting the value of a rotational speed disturbance of 0.004, −0.004. Figure 13, respectively, gives the rotational speed As shown in Figure 12a, the FFOPID controller can control the mechanical torque to a stable state faster than the PID controller and FOPID controller. It can also be seen from Figure 12b that under the action of the FFOPID controller, the guide blade opening can be balanced quickly and rapidly, so the corresponding rotational speed curve first reaches a stable state. Therefore, it can be considered that the FFOPID controller is better. Meanwhile, the results show that the oscillation amplitude of the water hammer in the surge chamber obtained by the FFOPID controller is smaller than that obtained by other controllers, as shown in Figure 12c,d. The above conclusion indicates that the FFOPID controller can provide more stable operating conditions.

Robustness Analysis of the Proposed Method
The low-frequency disturbance is the key factor to the stability of the pumped storage unit regulating system under the condition of low water head. Furthermore, the water hammer effect caused by the controller's poor robustness is an important cause of PSU vibration. Therefore, in this paper, the robustness test of the FFOPID controller was carried out by adjusting the value of a rotational speed disturbance of 0.004, −0.004. Figure 13, respectively, gives the rotational speed response of the pumped storage unit regulating system when subjected to low-frequency disturbance at 526 m and 530 m water head during operation. It can be seen that the control performance of the FFOPID controller under low-speed disturbance is obviously better than that of traditional controllers. The experimental results show that the proposed FFOPID controller of the pumped storage unit regulating system shows sufficient robustness to hydraulic parameter variation of the system under primary frequency modulation.

Robustness Analysis of the Proposed Method
The low-frequency disturbance is the key factor to the stability of the pumped storage unit regulating system under the condition of low water head. Furthermore, the water hammer effect caused by the controller's poor robustness is an important cause of PSU vibration. Therefore, in this paper, the robustness test of the FFOPID controller was carried out by adjusting the value of a rotational speed disturbance of 0.004, −0.004. Figure 13, respectively, gives the rotational speed response of the pumped storage unit regulating system when subjected to low-frequency disturbance at 526 m and 530 m water head during operation. It can be seen that the control performance of the FFOPID controller under low-speed disturbance is obviously better than that of traditional controllers. The experimental results show that the proposed FFOPID controller of the pumped storage unit regulating system shows sufficient robustness to hydraulic parameter variation of the system under primary frequency modulation.

Conclusions
In order to improve the stability of the pumped storage unit regulating system when it encounters rotation speed disturbance in low water head operation, this paper presented an FFOPID controller based on the fuzzy theory and FOPID controller. Furthermore, a multi-objective optimization strategy combining fuzzy rules and control parameters was applied to the primary

Conclusions
In order to improve the stability of the pumped storage unit regulating system when it encounters rotation speed disturbance in low water head operation, this paper presented an FFOPID controller based on the fuzzy theory and FOPID controller. Furthermore, a multi-objective optimization strategy combining fuzzy rules and control parameters was applied to the primary frequency regulation of the pumped storage unit regulating system. In the simulation experiment, the FFOPID controller was compared with the FOPID controller, PID controller, and FPID controller. In order to verify its robustness, different speed disturbance tests were carried out under the condition of low water head. The experimental analysis results indicated that the FFOPID controller could effectively restrain the oscillation of rotational speed at low water head. In future work, we will apply the FFOPID controller to more operating conditions of the pumped storage unit regulating system to verify its control performance.