Hopf Bifurcation and Parameter Sensitivity Analysis of a Doubly-Fed Variable-Speed Pumped Storage Unit

: The doubly-fed variable speed pumped storage unit is a storage system suitable for joint operation with renewable energy sources to smooth the imbalance between renewable energy supply and electricity demand. However, its working principle and operation control are more complex than those of constant speed pumped storage. In this study, a nonlinear model of doubly-fed variable speed pumped storage units (VSPSUs) considering nonlinear characteristics of the head loss is established. The study ﬁnds that a supercritical Hopf bifurcation occurs in the system, and the area enclosed by the lower side of the bifurcation line and the coordinate axis is the stability domain of the system. The active power step perturbation from − 0.3 to 0.3 will gradually reduce the area of the stability domain and narrow the adjustable range of the control parameters. In addition, the sensitivity of the model full state variables and the primary and secondary relationships to the changes of subsystem parameters is analyzed systematically using the trajectory sensitivity. It is found that there is a large difference in the sensitivity of different state variables to the parameters. The state variables are much more sensitive to the transfer coefﬁcient of hydraulic turbine torque to guide vane opening, the unit inertia time constant, and the controller proportional gain change than other parameters, which are deﬁned as highly sensitive parameters. The receiver response time constant and the turbine ﬂow-to-head transfer coefﬁcient are the corresponding low-sensitivity parameters.


Introduction
As power system operators pursue low-carbon investments and long-term energy sustainability, large-scale development and a high proportion of variable renewable energy (VRE) access in the power system will become inevitable in the coming decades. In fact, renewables remained the fastest growing source of energy in buildings, increasing 4.1% annually on average between 2009 and 2019, and reached their highest recorded share in the global electricity mix in 2020-an estimated 29% [1]. The higher the share of generation from variable renewable energy sources, the more flexible the power system must be to maintain a daily or annual balance between supply and demand. To manage the long-term imbalance between VRE supply and electricity demand, scholars have investigated solutions to smooth the imbalance between VRE supply and electricity demand and evaluated their role in deep decarbonization of power generation using advanced power system investment and operation models. These systems include nuclear power plants and natural gas plants equipped with carbon capture and storage, flexible demand, battery energy storage, and long-term storage technologies [2][3][4]. Energy storage technologies can provide long-term and seasonal energy conversion, allow for the temporal separation of energy generation and consumption, and address the intermittency of variable renewables. Dimanchev et al. used a detailed capacity expansion and dispatch model to verify the optimal role of reservoir hydro as an energy storage resource. Therefore, pumped storage units (PSUs), which use both upstream and downstream reservoirs for energy conversion and storage, are currently large-scale energy storage systems that are widely used around the world [5]. PSUs do not involve chemical conversion processes, which may easily generate pollution. They are highly adjustable to flexible and variable operating conditions, and have a combined operating efficiency of 70-85% [6]. Therefore, PSUs as energy storage devices are well suitable for combined operation with renewable energy in hybrid energy power systems [7].
There are two common types of PSU, namely constant speed pumped storage units (CSPSUs) and variable speed pumped storage units (VSPSUs) [8]. Currently, most of the PSUs are of constant speed. However, for CSPSUs, the speed regulation of the units is slow in the generation condition, and the input power of the unit is not adjustable in the pumping condition. Compared with the CSPSU, recent studies have demonstrated that the VSPSU not only inherits the advantages but also has the following additional advantages: (a) the active power of VSPSUs can be controlled quickly; (b) the frequency stability of the grid can be greatly improved; (c) and the efficiency of the units can be improved in both power generation and pumping conditions [9][10][11]. However, the working principle and operation control of VSPSUs are more complicated than those of CSPSUs. Therefore, the understanding of the dynamic behavior of VSPSUs is now in urgent need of enhancement, and the stability of VSPSUs has significant research value.
At present, the research on doubly-fed variable-speed pumped storage units (VSPSUs) is focused on the following aspects: mathematical modeling and simulation, dynamic characteristics and operational stability and reliability, and optimal control. Their contributions are mainly as follows. Gao et al. developed a fast and high-precision model of a VSPSU, which can better characterize VSPSUs [12]. Mohanpurkar et al. achieved real-time co-simulation of hydrodynamics and electrical events for adjustable-speed pumped storage hydro [13]. Zhao et al. established a model for the power generation and pumping conditions of VSPSUs and verified the model and the proposed control strategy using the MATLAB/Simulink simulation platform [14]. Xie et al. studied an electromechanical transient model of the VSPSU, preserving the dynamic process of the motor flux linkage and the dynamic characteristics of the controller inner loop [15]. Zhao et al. established the output power model of VSPSUs, and found that the VSPSU has good compensation ability for the power fluctuation of a power system [16]. Yang et al. evaluated the performance of the VSPSU in mitigating the power regulation of wind power variations and found that it excels in promoting power system stability [17]. Joseph et al. investigated the dynamic performance of VSPSUs in the event of power and control circuit failures as well as converter and sensor failures, providing ideas for the stable operation of VSPSUs [18,19]. Damdoum et al. proposed new simple fault ride-through strategies to reduce the negative impacts of grid fault occurrence on the doubly-fed induction machine pumped storage system [20]. Madeira et al. investigated the effect of sudden changes in excitation capacitances, resistive loads, or recovered head on the operation of a pump working as a turbine-self-excited induction generator model (PAT-SEIG) system using a constructed PAT-SEIG model [21]. Pagaimo et al. studied the transient characteristics of series-connected pumps working as turbines in an off-grid system and found that one change in the first pump of a turbine group will significantly affect the other group dynamics [22]. Gao et al. studied the stability of a VSPSU regulation system and found that the stability domain of the pumping condition is larger than that of the generation mode when the parameters of the power controller are fixed [23]. The above literature provides a basis for studying the stability of VSPSUs but does not investigate the effect of the intrinsic nonlinear characteristics of the system components on the stability of the whole unit.
Meanwhile, Cruz et al. showed that stability robustness has a close relationship with sensitivity robustness in nonlinear feedback systems [24]. Therefore, parameter sensitivity analysis is also an important part of system stability analysis. Some of the current studies focus on the influence of individual subsystem parameters (e.g., turbine system parameters [25]) on individual system characteristics (e.g., unit vibration [26]), and some scholars have also studied the influence of unit parameters on control characteristics using the zero-pole variation of the closed-loop transfer function and have used it to design the unit parameters [27]. Liu et al. analyzed the sensitivity of the hydropower unit damping to the system parameters and revealed the influence mechanism of the parameters on ultralow-frequency oscillations [28]. These studies lack a systematic study of the primary and secondary relationships between the effects of parameters on the system, and there are few examples of their application to VSPSUs.
Therefore, this paper focuses on the nonlinear bifurcation and parameter sensitivity characteristics of the VSPSU. In Section 2, the process of deriving the models of each subsystem is described, and a fifth-order nonlinear model of the VSPSU considering pressure head loss is established. A parameter sensitivity analysis method is also introduced in this section. In Section 3, nonlinear bifurcation analysis of the VSPSU under PI control strategy is performed to find its stability domain and analyze the bifurcation types of the system. The dynamic response of the system under active power perturbation is simulated numerically, and the conclusions obtained from the theoretical analysis are verified by the simulation. In Section 4, the effects of different active power steps on the stability domain of the system and the sensitivity of the system to the changes of key parameters of each subsystem are analyzed. Conclusions are drawn in the final section.

Mathematical Modeling and Parameter Sensitivity Analysis Method of Doubly-Fed Variable Speed Pumped Storage Units
In the field of automatic control, we often use the state space method to analyze systems, where the dynamic characteristics of the system are described by a number of first-order differential equations composed of state variables. They can reflect the changes of all independent variables of the system, so that all the internal motion states of the system can be determined at the same time, and they can also handle the initial conditions easily. In addition, they can be used in nonlinear systems, time-varying systems, multipleinput-multiple-output systems, and stochastic processes because they can be analyzed and designed with computers and controlled in real time. A VSPSU is a nonlinear system with hydraulic-mechanical-electrical coupling, for which a reliable and realistic mathematical model can be developed using the state space method. The VSPSU mainly consists of a governor, electro-hydraulic servo system, pump turbine, water diversion system and doubly-fed induction motor (DFIM), and the structure diagram of a VSPSU is shown in Figure 1. ity analysis is also an important part of system stability analysis. Some of the curren studies focus on the influence of individual subsystem parameters (e.g., turbine system parameters [25]) on individual system characteristics (e.g., unit vibration [26]), and som scholars have also studied the influence of unit parameters on control characteristics us ing the zero-pole variation of the closed-loop transfer function and have used it to desig the unit parameters [27]. Liu et al. analyzed the sensitivity of the hydropower un damping to the system parameters and revealed the influence mechanism of the param eters on ultra-low-frequency oscillations [28]. These studies lack a systematic study of th primary and secondary relationships between the effects of parameters on the system and there are few examples of their application to VSPSUs. Therefore, this paper focuses on the nonlinear bifurcation and parameter sensitivit characteristics of the VSPSU. In Section 2, the process of deriving the models of eac subsystem is described, and a fifth-order nonlinear model of the VSPSU considerin pressure head loss is established. A parameter sensitivity analysis method is also intro duced in this section. In Section 3, nonlinear bifurcation analysis of the VSPSU under P control strategy is performed to find its stability domain and analyze the bifurcatio types of the system. The dynamic response of the system under active power perturba tion is simulated numerically, and the conclusions obtained from the theoretical analysi are verified by the simulation. In Section 4, the effects of different active power steps o the stability domain of the system and the sensitivity of the system to the changes of ke parameters of each subsystem are analyzed. Conclusions are drawn in the final section.

Mathematical Modeling and Parameter Sensitivity Analysis Method of Doubly-Fed Variable Speed Pumped Storage Units
In the field of automatic control, we often use the state space method to analyz systems, where the dynamic characteristics of the system are described by a number o first-order differential equations composed of state variables. They can reflect th changes of all independent variables of the system, so that all the internal motion states o the system can be determined at the same time, and they can also handle the initial con ditions easily. In addition, they can be used in nonlinear systems, time-varying systems multiple-input-multiple-output systems, and stochastic processes because they can b analyzed and designed with computers and controlled in real time. A VSPSU is a non linear system with hydraulic-mechanical-electrical coupling, for which a reliable and re alistic mathematical model can be developed using the state space method. The VSPSU mainly consists of a governor, electro-hydraulic servo system, pump turbine, water d version system and doubly-fed induction motor (DFIM), and the structure diagram of VSPSU is shown in Figure 1.

Modeling of Hydraulic Subsystem
In building the mathematical model of a VSPSU, we need to consider a combinatio of hydraulic-mechanical-electrical factors. The hydraulic part of the factors is mainl represented by the pressure pipe used to transmit water, whose dynamic characteristic

Modeling of Hydraulic Subsystem
In building the mathematical model of a VSPSU, we need to consider a combination of hydraulic-mechanical-electrical factors. The hydraulic part of the factors is mainly represented by the pressure pipe used to transmit water, whose dynamic characteristics can be described using an ordinary differential equation. Assuming that the walls of the pressure pipe and the water in the pressure pipe are rigid, the water hammer pressure due to the variation of the guide vane opening of the pump turbine can be calculated by the rigid water column theory. Based on Newton's second law of motion, the fluid motion in the pressure pipe can be described by the following equation: where L t is the length of penstock, V t0 is the initial value of flow velocity in penstock, H is the turbine working head, H 0 is the initial value of the turbine working head, g is acceleration of gravity, Q t is the reference flow rate of the turbine, Q t0 is the initial value of the reference flow rate of the turbine. h t is the head loss of penstock, h t0 is the initial value of the head loss of penstock.
Since the state variables in this paper are expressed in the form of relative values of deviation, the deviation of the above equation regarding flow rate and water pressure can be simplified by h = H−H 0 H 0 and q = Q t −Q t0 Q t0 . L t V t0 gH 0 represents the characteristics of the inertia of water flow in the pressure pipe. We can use a time constant of water inertia T W to simplify the substitution. Then the simplified equation is obtained as follows: By associating with the flow rate at this moment, we can obtain the differential equation containing only h and q. For the head loss in penstock, the expression for the head loss is h t = αQ 2 t . At the initial moment, that is t = 0, when there is h t0 = αQ 2 t0 , then we , using the relationship between Q t and q, and h t = h t0 (1 + 2q + q 2 ) can be obtained by substituting Q t = Q t0 + Q t0 q. Therefore, we can obtain the equation at the initial moment: h t − h t0 = h t0 (2q + q 2 ).
The above equations can be combined to obtain the dynamic equation considering the pressure pipe head loss as follows [29]: This shows that the head loss of the penstock is nonlinear.

Modeling of Mechanical Subsystem
The mechanical part of the factors is reflected in the prime mover of hydroelectric power generation, i.e., the pump turbine, and the governor used to control the turbine speed and thus ensure the frequency stability of the power generation. In the modeling of the pump turbine, to represent the dynamic characteristics of the turbine in a relatively simple form, and considering that the variation of the turbine speed, head and guide vane opening is small in the study of this paper, we assume that the nonlinear relationship of the pump turbine is linear. The moment equation and flow equation of the pump turbine are used for further analysis and research, and the specific expressions of the equations are as follows: m t = e h h + e x x + e y y q = e qh h + e qx x + e qy y where x = n−n 0 n 0 and y = a−a 0 a 0 ; e h , e x , e y are the transfer coefficients of pump turbine torque to head, speed, guide vane opening respectively, e qh , e qx , e qy are the transfer coefficients of pump turbine flow to head, speed, guide vane opening respectively, n is the pump turbine speed, n 0 is the initial value of speed, a is the guide vane opening of pump turbine, a 0 is the initial value of guide vane opening.
The governor is the controller of the unit, and its control law commonly used in industrial applications is mainly PI control or PID control. PI control is used in this paper. In order to achieve the control goal, an electrohydraulic servo system is required to convert the output of the PI control strategy into the action of the guide vane, i.e., the guide vane opening. This conversion is generally in the form of integration, supplemented by the conversion coefficient of the response. The equations of this control strategy and the electrohydraulic servo system are: where u is the governor regulation output, K p is the proportional gain, K i is the integral gain, and T y is the receiver response time constant.

Modeling of Electrical Subsystem
The electrical factors are concentrated in the DFIM of a VSPSU. The mathematical model of the DFIM consists of voltage equations, magnetic chain equations, torque equations and equations of motion. The mathematical model of the DFIM in the conventional three-phase stationary coordinate system is a very complex and strongly coupled timevarying system. To facilitate the analysis and solution of the operating characteristics, the mathematical model of the DFIM is transformed by using the coordinate transformation method, and the mathematical model in the three-phase stationary abc coordinate system is equated to that in two identical rotating dq coordinate systems.
The transformed voltage and magnetic chain equations are as follows: where u ds , u qs , u dr , u qr are the d, q axis voltage of stator and rotor, respectively; R s , R r is the stator and rotor resistance, respectively; i ds , i qs , i dr , i qr are the d, q axis current of stator and rotor, respectively; Ψ ds , Ψ qs , Ψ dr , Ψ qr are the d, q axis flux linkage of stator and rotor, respectively; ω s denotes the synchronous speed, s r = ω s −ω r ω s denotes the relative speed deviation. L m is the mutual inductance between stator and rotor; L ss , L rr are the inductance of stator and rotor respectively; subscripts s and r represent the stator and rotor components respectively; subscripts d and q represent the d-axis and q-axis components respectively.
The equation of motion and the equation of torque can be described as follows: where T a is the unit inertia time constant; T m is the active torque of the turbine; and T e is the load torque of the generator; p n is the number of poles. When the model of the DFIM is converted to two identical steps in a rotation coordinate system, vector control is generally used as the control method. We can let the d-axis coincide with the stator voltage vector, while ignoring the stator resistance which is much smaller than the reactance, and finally obtain the stator magnetic chain vector in steady state with 90 • phase difference from the stator voltage vector. After orientation, a series of updated system equations can be obtained. The most important power equation and torque equation will be updated in Equation (9): where U s is effective value of grid voltage; P s is active power; Q s is reactive power.
x − e qy e qh By combining the above equations and selecting the independent state variables q, x, y, u, i dr that represent the characteristics of the VSPSU, with appropriate variable substitutions and simplifications, we can obtain a fifth-order nonlinear model for a VSPSU as Equation (10).

Parameter Sensitivity Analysis Method
In this paper, the definition of trajectory sensitivity and average trajectory sensitivity in the literature [30] are used to analyze the parameters in the system accordingly, and the nonlinear model of a VSPSU is expressed as Equation (10). Under certain operating conditions, taking the rotational speed of VSPSU as an example, the trajectory sensitivity of the system state variables with respect to a parameter can be defined as Equation (11): where, x 0 is the rotational speed of the VSPSU when the value of the parameter θ i is θ i0 , ∆θ i is the variation of the θ i , and θ r is the other parameters in the model except for the θ i . In order to facilitate the comparison and analysis of the parameter trajectory sensitivity, the relative trajectory sensitivity is defined, and its absolute average value is taken as the index for quantitative analysis. The specific expression is shown in Equation (12): where t 1 , t 2 are the start and end moments of a dynamic process, respectively, which are reasonably chosen according to the response dynamic curve of the system.

Hopf Bifurcation Theory
Power systems consist of various nonlinear parts, including synchronous generators (providing active and reactive power for the networks), load buses (representing power consumers), and distribution transmission lines. The dynamic characteristics of those systems are represented by the swing equation of synchronous generators, and derivative and algebraic equations of dynamic load (such as induction motors), thus proving the nonlinear load characteristics of the power system [31]. In actual engineering problems, the parameters on which the dynamic system depends often undergo arbitrarily small changes (also known as perturbations). When these parameters change, the phase trajectory topology of the unstable system will undergo essential changes. This situation is called bifurcation. As the load conditions change, the system with nonlinear derivative equations undergoes a qualitative change at the bifurcation point. Hopf bifurcation and saddle node bifurcation are considered as typical bifurcations in power systems. Hopf bifurcation theory is applied in this paper to study the nonlinear bifurcation characteristics of VSPSUs.
Consider a general nonlinear system .
x as follows: .
where x is the n dimensional state variable, and µ is the m dimensional parameter vector. When the parameter µ passes a certain critical value µ 0 , the phase trajectory topology of the system suddenly changes at µ = µ 0 , and the system is called bifurcated at µ = µ 0 , and µ 0 is called the bifurcation value. The totality of bifurcation values is called the bifurcation set. When the parameters change, the dynamic behavior of the nonlinear system switches between a stable equilibrium point and a stable limit cycle. This dynamic evolution process is called Hopf bifurcation.
Assume that x = x 0 is the equilibrium point of the system: that is, f (x 0 , µ) = 0 is satisfied. Suppose the Jacobian matrix of the system at the equilibrium point x 0 is J(µ) = DF x (x 0 , µ). Expand the characteristic equation det(sI − J(µ)) = 0 of the Jacobian matrix and arrange them in descending powers to obtain Equation (14): where c i (µ)(i = 1, 2, · · · , n) is the coefficient of the characteristic equation, and for a certain bifurcation parameter value µ = µ c , if the following conditions are met, the n-dimensional nonlinear dynamic system has Hopf bifurcation: there is c i = 0. At this moment, the characteristic equation has a pair of pure imaginary roots at µ = µ c ; Then the system will have Hopf bifurcation at µ = µ c , and the period of periodic motion at this time is T = 2π ω . In addition, the transversal coefficient σ (µ c ) can also be used to distinguish the type of bifurcation. When σ (µ c ) > 0, the Hopf bifurcation that occurs is supercritical, and at µ < µ c , the equilibrium point of the system is the stable focus, and for a sufficiently small µ which satisfies µ > µ c , the system will bifurcate from the equilibrium position and perform periodic motion. In the phase trajectory, it will produce a stable limit cycle, and then the system will enter a continuous oscillation. When σ (µ c ) < 0, the Hopf bifurcation that occurs is subcritical. At this time, the system shows stable limit cycles and stable focal motion characteristics at the left and right sides of µ c respectively.

Hopf Bifurcation Analysis of VSPSU
For the fifth-order VSPSU nonlinear model established in this paper, firstly, the equilibrium point of the nonlinear dynamic system can be described as and it can be calculated based on the nonlinear state equation f (X B , µ) = 0. The result is shown in Equation (15).
Secondly, the Jacobian matrix of the nonlinear dynamic system is formulated at the equilibrium point to determine the existence of Hopf bifurcation. The Jacobian matrix of the system at its equilibrium point X B is shown in Equation (16). The detailed partial derivatives in the Jacobi matrix are given in Appendix A.
According to the Hopf bifurcation theory introduced earlier, the characteristic equation of J(µ) can be formulated from det(sI − J(µ)) = 0 and the result is shown in Equation (17): where s is the eigenvalue of det(sI − J(µ)) = 0, and the expressions of c i (i = 1, 2, . . . , 5) are presented in Appendix B.
For the VSPSU nonlinear system in this paper, the bifurcation point where the Hopf bifurcation occurs is the edge point between the stability and the instability of the system. In the coordinate plane composed of several parameters (K p , K i ), we consider the set of bifurcation points in different states will form a curve in the coordinate plane. The bifurcation line can divide the stable region and unstable region of the system. Therefore, the position of the bifurcation line can largely reflect the dynamic characteristics of the system. Assuming that the system is affected by perturbation (the aforementioned perturbation), in the stable region, the dynamic response of the system will converge to a new steady state; on the bifurcation line, the response of the system will oscillate with the same amplitude; and in the unstable region, the response of the system will diverge, and a new steady state cannot be obtained. According to those principles, the bifurcation line of the system can be drawn.
Select K i as the bifurcation parameter of the system, and the corresponding bifurcation point can be expressed as µ c = K * i ; select the proportional-integral parameter plane, and K p , K i are the abscissa and ordinate, respectively, to determine the bifurcation line of the system. The disturbance is set to a 10% step drop of the rated value of the active power at the time t = 0 s; that is, P re f = −0.1. The rest of the parameters involved in the system are shown in Table 1. Use the inequality in the aforementioned criterion to determine the feasible range of the parameter, and the equation to obtain the relationship between the parameters. The determined bifurcation line of the system is shown in Figure 2.  At the same time, applying the aforementioned definition of transversal coefficient to the system described in this paper, the available system's transversal coefficient expression is shown in Equation (18):  At the same time, applying the aforementioned definition of transversal coefficient to the system described in this paper, the available system's transversal coefficient expression is shown in Equation (18): where c i = dc i dµ , (i = 1, · · · , 5), and the pure virtual characteristic root at µ = µ c is s 1,2 = ±i c 5 −c 1 c 4 c 3 −c 1 c 2 , and the values of the transversal coefficient corresponding to all bifurcation points are shown in Figure 3.  It can be clearly found that all the transversal coefficients are bigger than 0, which means that the bifurcation in the fifth-order VSPSU nonlinear system is supercritical. It can be clearly found that all the transversal coefficients are bigger than 0, which means that the bifurcation in the fifth-order VSPSU nonlinear system is supercritical. That is, the area enclosed by the coordinate axis under the bifurcation line is the stable region of the system.

Numerical Stability Simulation
In order to verify the accuracy of the stable region drawn above, 5 points S1-S5 are selected in Figure 2 for the numerical simulation of the dynamic response, and the values of K p , K i of the selected five points are shown in Table 2. Using the ode45 solving function in MATLAB software to solve the nonlinear derivative equations for numerical simulation, we can calculate the dynamic response process of the characteristic state variables q, x, y corresponding to the five different points in Table 2 under the active power disturbance, and the phase trajectory of the dynamic response of these variables as well. The results are shown in Figure 4.
From the results in Figure 4, we can conclude that the dynamic responses and phase trajectories of the selected four state point variables are consistent with the results obtained from the Hopf bifurcation theory analysis. S1 and S2 are located in the stable domain of the coordinate plane formed by the parameters. After the system is disturbed by the change of active power, the dynamic response of characteristic variables (i.e., flow, frequency, opening) needs to go through several cycles of response attenuation and oscillation before reaching the final equilibrium state. Correspondingly, the three-dimensional phase composed of characteristic variables will have several trajectories before reaching the equilibrium point; S3 is located on the bifurcation line of the coordinate plane formed by the parameters, the dynamic response of the characteristic variables will enter a continuous oscillation state with constant amplitude, and the corresponding phase trajectory is a stable limit cycle. For S4 and S5, which are located in the unstable domain, the characteristic variables will first go through several cycles of divergent response state before entering a stable constant amplitude oscillation, and the phase space trajectory is divergent first and then enters a stable limit cycle.
Comparing S1 and S2, the smaller the value of K i and the larger the distance ( K i − K * i ) from the bifurcation line, the shorter the time for the system to return to the equilibrium state after the disturbance. This means that the further the points in the stability domain are from the bifurcation line, the better the stability of the system. Similarly, comparing S4 and S5, the larger the value of K i and the greater the distance ( K i − K * i ) from the bifurcation line, the shorter the time from divergence to the equal amplitude oscillation, and the larger the amplitude of the oscillation in the constant amplitude oscillation. In other words, the system becomes less stable when the points in the unstable domain are far from the bifurcation line. It can be seen that the improper selection of the control parameters of the system will cause the CSPSU to enter a state of constant amplitude oscillation. In order to protect the safe and stable operation of the system, it is necessary to guide the selection of the PI controller parameters according to the reference stability domain.
At the same time, by substituting the pure imaginary characteristic root of the system into the calculation formula of the oscillation period and then taking the reciprocal, the oscillation frequency of the point on the bifurcation line can be obtained. The specific formula is shown in Equation (19).
Calculate the oscillation frequency corresponding to all the bifurcation points, and carry out the dynamic response simulation at the same time to obtain the amplitude value of the state variable x in the system during constant amplitude oscillation. The result is shown in Figure 5. The three equal-amplitude oscillation points P1, P2, P3 are selected to verify the results in Figure 5, and the control parameters are shown in Table 3, while the simulated curves are obtained as shown in Figure 6. That is, the area enclosed by the coordinate axis under the bifurcation line is the stable region of the system.

Numerical Stability Simulation
In order to verify the accuracy of the stable region drawn above, 5 points s1-s5 are selected in Figure 2 for the numerical simulation of the dynamic response, and the values of 、 p i K K of the selected five points are shown in Table 2. Using the ode45 solving function in MATLAB software to solve the nonlinear derivative equations for numerical simulation, we can calculate the dynamic response process of the characteristic state variables 、 、 q x y corresponding to the five different points in Table 2 under the active power disturbance, and the phase trajectory of the dynamic response of these variables as well. The results are shown in Figure 4.    From the results in Figure 4, we can conclude that the dynamic responses and phase trajectories of the selected four state point variables are consistent with the results obtained from the Hopf bifurcation theory analysis. S1 and S2 are located in the stable domain of the coordinate plane formed by the parameters. After the system is disturbed by the change of active power, the dynamic response of characteristic variables (i.e., flow, frequency, opening) needs to go through several cycles of response attenuation and oscillation before reaching the final equilibrium state. Correspondingly, the three-dimensional phase composed of characteristic variables will have several trajectories before reaching the equilibrium point; S3 is located on the bifurcation line of the coordinate plane formed by the parameters, the dynamic response of the characteristic variables will enter a continuous oscillation state with constant amplitude, and the corresponding phase trajectory is a stable limit cycle. For S4 and S5, which are located in the unstable domain, the characteristic variables will first go through several cycles of divergent response state before entering a stable constant amplitude oscillation, and the phase space trajectory is divergent first and then enters a stable limit cycle. Comparing S1 and S2, the smaller the value of i K and the larger the distance ( * ii KK − ) from the bifurcation line, the shorter the time for the system to return to the equilibrium state after the disturbance. This means that the further the points in the stability domain are from the bifurcation line, the better the stability of the system. Similarly, comparing S4 and S5, the larger the value of i K and the greater the distance ( from the bifurcation line, the shorter the time from divergence to the equal amplitude oscillation, and the larger the amplitude of the oscillation in the constant amplitude oscillation. In other words, the system becomes less stable when the points in the unstable domain are far from the bifurcation line. It can be seen that the improper selection of the control parameters of the system will cause the CSPSU to enter a state of constant amplitude oscillation. In order to protect the safe and stable operation of the system, it is necessary to guide the selection of the PI controller parameters according to the reference stability domain. At the same time, by substituting the pure imaginary characteristic root of the system into the calculation formula of the oscillation period and then taking the reciprocal, the oscillation frequency of the point on the bifurcation line can be obtained. The specific formula is shown in Equation (19).
Calculate the oscillation frequency corresponding to all the bifurcation points, and carry out the dynamic response simulation at the same time to obtain the amplitude value of the state variable x in the system during constant amplitude oscillation. The result is shown in Figure 5. The three equal-amplitude oscillation points P1, P2, P3 are selected to verify the results in Figure 5, and the control parameters are shown in Table 3, while the simulated curves are obtained as shown in Figure 6.    Figure 5, it can be clearly seen that with the change of p K from small to large, the oscillation frequency presents a monotonically increasing trend, becoming larger in a nearly linear manner. The amplitude of the oscillation changes in a more tortuous manner, while the overall trend is downward. There is a sharp drop in amplitude near the intersection of the bifurcation line with the horizontal axis. On the whole, the amplitude and frequency of the system's constant-amplitude oscillation have roughly the opposite trend; that is, the greater the frequency, the smaller the amplitude, and vice versa. This may indicate from the side that the energy of the system is generally conserved. The results shown in Figure 6 are the dynamic response of state variable x at points P1-P3. These response curves and the frequency marked on the right side of the figure verify the previous frequency-amplitude diagram.

The Effect of Active Power Step on the Stability Domain
The aforementioned stable region is drawn when the active power has a step drop of 10% of the rated value. In actual engineering, the magnitude and direction of the disturbances are unpredictable. For this reason, the influence of different disturbances on the system stable region is studied. We selected a total of six active power disturbances with different directions and magnitudes, and drew the stable region according to the same method as described above, and the results are shown in Figure 7. From Figure 5, it can be clearly seen that with the change of K p from small to large, the oscillation frequency presents a monotonically increasing trend, becoming larger in a nearly linear manner. The amplitude of the oscillation changes in a more tortuous manner, while the overall trend is downward. There is a sharp drop in amplitude near the intersection of the bifurcation line with the horizontal axis. On the whole, the amplitude and frequency of the system's constant-amplitude oscillation have roughly the opposite trend; that is, the greater the frequency, the smaller the amplitude, and vice versa. This may indicate from the side that the energy of the system is generally conserved. The results shown in Figure 6 are the dynamic response of state variable x at points P1-P3. These response curves and the frequency marked on the right side of the figure verify the previous frequency-amplitude diagram.

The Effect of Active Power Step on the Stability Domain
The aforementioned stable region is drawn when the active power has a step drop of 10% of the rated value. In actual engineering, the magnitude and direction of the disturbances are unpredictable. For this reason, the influence of different disturbances on the system stable region is studied. We selected a total of six active power disturbances with different directions and magnitudes, and drew the stable region according to the same method as described above, and the results are shown in Figure 7.
It can be seen from Figure 7 that as the active power disturbance changes from −0.3 to 0.3, the intersection points between the stable domain and the horizontal axis (i.e., K p ) gradually shifts to the left, and the intersection point with the vertical axis (i.e., K i ) gradually moves downward. The area of the stable region also gradually decreases, and the specific values are shown in Table 4. This means that the adjustable range of the PI controller parameters is also reduced for the stable operation of the system. This may be manifested in actual engineering operation: when the VSPSU is operating with 70% rated load under power generation conditions, the pressure to ensure its stable operation is lower than when operating with a rated load.   Table 4. This means that the adjustable ra controller parameters is also reduced for the stable operation of the system. manifested in actual engineering operation: when the VSPSU is operating wi load under power generation conditions, the pressure to ensure its stable lower than when operating with a rated load. To facilitate the calculation of the above indicators, the simulation co designed to subject the VSPSU to an active power disturbance, setting the act 10% of its rated value with a pulse width of 0.1 s. Under this disturbance, the will fluctuate for a short time from the initial equilibrium point and return to

System Parameter Trajectory Sensitivity Analysis Process
To facilitate the calculation of the above indicators, the simulation conditions are designed to subject the VSPSU to an active power disturbance, setting the active power to 10% of its rated value with a pulse width of 0.1 s. Under this disturbance, the system state will fluctuate for a short time from the initial equilibrium point and return to the original equilibrium point again under the regulation of the controller. The reason why the form of perturbation is not chosen from the aforementioned step perturbation used to calculate the stability domain is, on the one hand, to facilitate the calculation of the sensitivity index; on the other hand, the equilibrium point of the system will change after the step perturbation, which means that the parameter sensitivity will be affected by both the dynamic characteristics and the steady-state characteristics, while the study in the section is actually concerned with the effect of the parameter change on the dynamic characteristics of the system. According to the formula of trajectory sensitivity, it is only necessary to calculate the area of the dynamically changing part of the system variables. The sensitivity of the system state variables to the parameters can be obtained by comparing the difference between the areas before and after the parameter variation. The flow chart of parameter sensitivity analysis is shown in Figure 8.

System Parameter Trajectory Sensitivity Results
The parameters of the system that will vary within a reasonable range during the dynamic transients of the system are w , , , , , , , , , , According to Figure 7, a ± 5% perturbation is applied to each parameter on the basis of its initial value, and the simulation test of active power pulse disturbance is performed. The total duration of the simulation is 200 s, in which the disturbance occurs at  1s t . Due to the large number of parameters, only the parameter with the greatest influence on the system state (i.e., a T ) is selected to analyze the influence of the parameter perturbation on the system state variables. Taking the rotational speed of the VSPSU as an example, the simulation curve before and after the parameter ingestion is obtained, which is shown in Figure 9a, and the corresponding trajectory sensitivity calculation results are shown in Figure 9b (only the trajectory sensitivity of the first 100 s is selected to optimize the display effect).

System Parameter Trajectory Sensitivity Results
The parameters of the system that will vary within a reasonable range during the dynamic transients of the system are e x , e y , e h , e qx , e qy , e qh , T w , T a , T y , K p , K i , and the system variables to be analyzed are the five state variables that are selected when the model is built.
According to Figure 7, a ± 5% perturbation is applied to each parameter on the basis of its initial value, and the simulation test of active power pulse disturbance is performed. The total duration of the simulation is 200 s, in which the disturbance occurs at t = 1 s. Due to the large number of parameters, only the parameter with the greatest influence on the system state (i.e., T a ) is selected to analyze the influence of the parameter perturbation on the system state variables. Taking the rotational speed of the VSPSU as an example, the simulation curve before and after the parameter ingestion is obtained, which is shown in Figure 9a, and the corresponding trajectory sensitivity calculation results are shown in Figure 9b (only the trajectory sensitivity of the first 100 s is selected to optimize the display effect).
From Figure 9a, it can be seen that the positive perturbation of the T a slows down the decrease rate of the speed to a small extent and reduces the minimum and maximum during the dynamic process of the speed; i.e., at the same time when the speed decreases, its value subject to the positive perturbation of T a is larger than that subject to the negative regression of T a , and vice versa. As can be seen from Figure 9b, the sensitivity of different state variables to the T a varies, and the sensitivity of variables q, y, u to parameter T a is significantly greater than that of variable x, and the sensitivity of y, u to T a is basically the same, which can also be clearly seen from the derivative equations of the VSPSU. The sensitivity indexes of the complete system state variables to each parameter are listed in Table 5, and the depth of the color bars in the cells can visually reflect the size of the sensitivity indexes.  From Figure 9a, it can be seen that the positive perturbation of the a T slows down the decrease rate of the speed to a small extent and reduces the minimum and maximum during the dynamic process of the speed; i.e., at the same time when the speed decreases, its value subject to the positive perturbation of a T is larger than that subject to the negative regression of a T , and vice versa. As can be seen from Figure 9b, the sensitivity of different state variables to the a T varies, and the sensitivity of variables , , q y u to parameter a T is significantly greater than that of variable x, and the sensitivity of , y u to a T is basically the same, which can also be clearly seen from the derivative equations of the VSPSU. The sensitivity indexes of the complete system state variables to each parameter are listed in Table 5, and the depth of the color bars in the cells can visually reflect the size of the sensitivity indexes.  It can be seen from Table 5 that the sensitivity of each state variable to the e y , T a , K p is much higher than that of the other parameters; i.e., e y , T a , K p are the high sensitivity parameters, and the corresponding low sensitivity parameters are T y , e qh . In addition to the aforementioned highly sensitive parameters, the state variable q is also highly sensitive to the changes of the e qy , and the sensitivity of u, y is almost the same for different parameters. Moreover, the sensitivity of the state variable i dr is almost the same for different parameters.

Conclusions
As a long-term energy storage system, it is a very meaningful topic to study the stability and dynamic characteristics of the VSPSU at a time when the share of variable renewable energy access is continuously increasing.
The focus of this study was on the nonlinear bifurcation characteristics and parameter sensitivity of the VSPSU. Considering the nonlinear characteristics of head loss, a fifth-order nonlinear mathematical model of the VSPSU was developed. The numerical simulation verified the Hopf bifurcation theory derivation. When the system was subjected to active power step disturbance, the bifurcation in the system was supercritical. Considering the points located on the bifurcation line, the system perturbation response showed the characteristics of equal-amplitude oscillation, and when the PI controller proportional gain varied from small to large, the oscillation frequency was between 0.02 and 0.12, and showed a monotonic increasing trend. The amplitude of oscillation was opposite to the trend of frequency. As the active power step disturbance changed from −0.3 to 0.3, the area of the stability domain gradually decreased from 18.7364 to 17.3150.
Parameter sensitivity analysis of the established model yielded the correlation between the system full-state variables and subsystem parameter variations. The system state variables were more sensitive to changes in the transfer coefficient of turbine torque to guide vane opening, the unit inertia time constant, and the controller proportional gain, while they were less sensitive to changes in the receiver response time constant and the turbine flow-to-head transfer coefficient.
The stability domain can be used to guide the selection of unit control parameters for operation of the VSPSU, and the variation of the stability domain area with disturbance can provide a basis for rational load planning during the actual operation of the VSPSU, thus helping to exploit its ability to solve the intermittent problems of renewable energy sources.