As the actuator and core part of the HTRS, the governor which includes a controller and a servo system can balance the active power of power systems by adjusting the active power output from hydropower stations. The controller model applied in this paper is a parallel PI control structure, which can be represented as follows: d y d t = − K p d x t d t − K i x t , (1) where y is the output signal of the HURS governor; x _t is the relative deviation of the speed; and K _p and K _i represent the proportional gain and the integral gain, respectively.

The rigid model and elastic water hammer model are often utilized for a mathematical model of the pressure pipeline. For hydropower stations with long pressure pipelines, the elastic water hammer model is preferable to describe the characteristics of water flow in the pipe. However, as to hydropower stations with short pressure pipelines, it can be simplified to the rigid model. In this paper, the non-linear characteristics of flow and head loss are considered to establish the dynamic equation of the elastic water hammer model of the pressure pipeline (Liu and Guo, 2021). d q t d t = − h − Z F − 2 h t q t H 0 T w t , (2) where q _t is the flow rate of the pressure pipe; h is the relative deviation of the water head; Z _F is the water level of the downstream surge chamber; h _t is the head loss; H ₀ is the initial water head; and T _wt is the inertial time constant of the water flow.

The dynamic equation of the downstream surge chamber is (Chaudhry, 2014) d Z F d t = q t − q y T F , (3) where Z is the water level change of the downstream surge chamber; gr represents the flow rate of the diversion tunnel; gy is the flow rate of the tailrace tunnel; TF is the time constant of the downstream surge chamber; and T _F is a time constant of the surge chamber, T F = F H 0 Q y 0 .

The dynamic equation of the sloping roof tailrace tunnel is (Chaudhry, 2014) d q y d t = Z F − 2 h y q y H 0 − z y T w y + T w x . (4)

According to the work of Guo et al. (2015), T _wx and z _y are derived from Eqs 5, 6, respectively, and then, the dynamic equations of the sloping roof tailrace tunnel can be obtained by substituting Eqs 5, 6 into Eq. 4 as shown in Eq. 7. T w x = λ Q y 0 V x q y g H 0 c W ⁡ tan ⁡ α , (5) z y = λ Q y 0 q y H 0 c W , (6) d q y d t = − 2 h y q y H 0 + Z F − q y Q y 0 λ H 0 c W T w y + q y Q y 0 λ V x ⁡ cot ⁡ α g H 0 c W , (7) where h _y is the water head loss of the tailrace tunnel; z _y is the water level change of the tailrace tunnel; T _wy and T _wx represent the time constant of steady-state flow inertia and transient flow inertia, respectively; λ is the section coefficient of the tailrace tunnel; c is the wave velocity at the interface of the open and full-flow manifolds; W is the width of the tailrace tunnel; V _x is the flow rate at the open and full-flow interfaces; and α is the inclination angle of the tailrace tunnel.

The moment equation and flow equation of the hydraulic turbine are (Yang et al., 2019b) m t = e x x t + e y y + e h h , q t = e q x x t + e q y y + e q h h . (8)

In previous research, the first-order model of a synchronous generator was commonly employed to describe its inertia moment. Given the strong mechanical inertia of the hydraulic relay and water guide system, when the load changes, it is hard for the actuators to track and adjust quickly in real time, which leads to a lag phenomenon (Ma et al., 2021). To better reveal the dynamic characteristics of the HTRS coupled with the PG, this paper employs a second-order model of the synchronous generator, which depicts not only the rotational inertia but also the interaction between the electromagnetic power and the power angle. The mathematical model of synchronous generator is as follows: d x t d t = m t − e g x t + K a ∫ x t − x s d t + D a x t − x s + m g T a . (9)

Since there is an integral term in Eq. 9, the state variable ξ ₁ is defined as ξ 1 = ∫ x t − x s d t . Equation 9 can be transformed into d ξ 1 d t = x t − x s d x t d t = m t − e g x t + K a ξ 1 + D a x t − x s + m g T a , (10) where ξ ₁ is the intermediate state variable; x _s is the relative deviation value of the grid frequency; e _g is the self-regulation coefficient of the load; K _a is the equivalent synchronization coefficient; D _a is the equivalent damping coefficient; m _g is the relative deviation of the resistance torque; and T _a is the inertial time constant of the unit.

The structural block diagram of the equivalent PG is shown in Figure 2. Load disturbance only considers mg, ignoring load disturbance in the PG, i.e., p _l = 0. The dynamic equations of the equivalent PG can be deduced from Figure 2 as shown in the following (Guo and Peng, 2020): d ξ 2 d t = x s − ξ 2 T g , d x s d t = B p t − D s x s − ξ 2 T g R g T s , (11) where ξ ₂ is the intermediate state variable; T _g is the inertia time constant of the servo motor in the grid model; B is the power conversion coefficient; p _t is the power output; D _s is the self-regulation coefficient of the equivalent load; difference coefficient; and R _g is the inertia time constant of the grid equivalent unit.

By integrating Eqs 1–4 and Eqs 9–11, the eighth-order non-linear state equations are obtained as Eq. 12, which can reflect the coupling effect of the HTRS and PG. d x t d t = 1 T a − m g − e g x t + e x x t − D a x t − x s + e y y + e h e q h q t − e q x x t − e q y y − K a ξ 1 d y d t = − K i x t − K p T a − m g − e g x t + e x x t − D a x t − x s + e y y + e h e q h q t − e q x x t − e q v y − K a ξ d ξ 1 d t = x t − x s d x s d t = 1 T s − D s − B D a x s + B D a x t + B K a ξ 1 − ξ 2 R g T g d ξ 2 d t = x s − ξ 2 T g d q t d t = 1 T w t − 2 h t q t H 0 − 1 e q h q t − e q x x t − e q v y − Z F d Z F d t = 1 F H 0 q t − q y Q y 0 d q y d t = g c W − 2 h y q y + H 0 Z F − q y Q y 0 λ c g H 0 T w y W + q y Q y 0 V x λ ⁡ cot ⁡ α . (12)

Hopf bifurcation (Strogatz, 2014; Wiggins, 2013) is a simple yet important dynamic bifurcation phenomenon in non-linear systems, which is a type of localized dynamic bifurcation; specifically, as the bifurcation parameter varies, the system bifurcates abruptly from the equilibrium at the non-hyperbolic equilibrium from the extremal limit loop phenomenon. Accordingly, Hopf bifurcation theory with theoretical applicability, simple program design, and high calculation accuracy is generally applied to conduct non-linear analysis. Meanwhile, Hopf bifurcation theory has been applied to the dynamics research of the HTRS in relevant works. Thus, this theory can be adopted in this paper to investigate the non-linear dynamical behavior of the coupled system.

Under external perturbations, the system will generate steady-state or unsteady-state limit-loop oscillations at Hopf bifurcation points, which correspond to supercritical and subcritical bifurcations, respectively. The Jacobian matrix of the coupled system at the equilibrium point X _E is J(μ) = Df _x (X _E ,μ). Accordingly, the characteristic equation of the Jacobian matrix can be derived by det(J(μ) - λI) = 0 as follows: λ n + a 1 λ n − 1 + a 2 λ n − 2 + ⋯ + a n − 1 λ + a n = 0 , (13) where λ and a _i (μ) (i = 1,2, ,n) are eigenvalues and coefficients of the polynomial det(J(μ) - λI) = 0, respectively.

The existence of Hopf bifurcation can be validated by the following well-known Hurwitz criterion (Hassard et al., 1981): a i > 0 i = 1,2 , ⋯ , n , (14) Δ 2 = det a 1 1 a 3 a 2 > 0 , (15) Δ 3 = det a 1 1 0 a 3 a 2 a 1 a 5 a 4 a 3 > 0 , (16) Δ 7 = det a 1 1 0 0 0 0 0 a 3 a 2 a 1 1 0 0 0 a 5 a 4 a 3 a 2 a 1 0 0 a 7 a 6 a 5 a 4 a 3 a 2 a 1 0 a 8 a 7 a 6 a 5 a 4 a 3 0 0 0 a 8 a 7 a 6 a 5 0 0 0 0 0 a 8 a 7 = 0 . (17)

If Eqs 14–17 are satisfied at μ = μ _c , Eq. 13 has a pair of pure virtual eigenvalues λ _1,2 = ±iω _c and the system will bifurcate at μ = μ _c. Accordingly, the system will occur with periodic oscillations and generate limit cycles in phase space. The limit cycle period at Hopf bifurcation is T _LC = 2π∕ω _c . Furthermore, the type of Hopf bifurcation can be determined based on the cross-sectional coefficient σ ′ μ c = Re d λ d μ μ = μ c .

According to Eq. 12, variables x _t , y, ξ ₁, x _s , ξ ₂, q _t , Z _F , and q _y are chosen as state variables of the coupled system, and state vector X = (x _t , y, ξ ₁, x _s , ξ ₂, q _t , Z _F , q _y ) ^T can be obtained. With the exception of state variables, any other parameter in Eq. 12 can be chosen as the bifurcation parameter, denoted as μ. The choice of bifurcation parameter μ is based on the purpose of system research. For example, K _p and K _i are chosen as bifurcation parameters when studying the effect of the governor on the stability of the whole system.

The equilibrium point X _E = (x _tE , y _E , ξ _1E , x _sE , ξ _2E , q _tE , Z _FE , q _yE )^T can be obtained by solving f(x, μ) = 0, i.e., by making Eq. 12 = 0: x t E = 0 y E = m g c H 0 + 2 e q h h t + h y W + e q h Q y 0 λ c − 2 e h e q y h t + h y + e y H 0 + 2 e q h h t + h y W + − e h e q y + e q h e y Q y 0 λ ξ 1 E = 0 x s E = 0 ξ 2 E = 0 q t E = − c e q v H 0 m g W 2 c e h e q y h t + h y W − c e y H 0 + 2 e q h h t + h y W + e h e q y − e q h e y Q y 0 λ Z F E = − e q y m g 2 c h y W + Q y 0 λ 2 c e h e q y h t + h y W − c e y H 0 + 2 e q h h t + h y W + e h e q y − e q h e y Q y 0 λ q y E = − c e q y H 0 m g W 2 c e h e q y h t + h y W − c e y H 0 + 2 e q h h t + h y W + e h e q y − e q h e y Q y 0 λ . (18)

The Jacobian matrix J(μ) for solving the non-linear coupled system x ˙ = f x , μ is J μ = D f x x , μ = ∂ x ˙ t ∂ x t ∂ x ˙ t ∂ y ∂ x ˙ t ∂ ξ 1 ∂ x ˙ t ∂ x s ∂ x ˙ t ∂ ξ 2 ∂ x ˙ t ∂ q t ∂ x ˙ t ∂ Z F ∂ x ˙ t ∂ q y ∂ y ˙ ∂ x t ∂ y ˙ ∂ y ∂ y ˙ ∂ ξ 1 ∂ y ˙ ∂ x s ∂ y ˙ ∂ ξ 2 ∂ y ˙ ∂ q t ∂ y ˙ ∂ Z F ∂ y ˙ ∂ q y ∂ ξ ˙ 1 ∂ x t ∂ ξ ˙ 1 ∂ y ∂ ξ ˙ 1 ∂ ξ 1 ∂ ξ ˙ 1 ∂ x s ∂ ξ ˙ 1 ∂ ξ 2 ∂ ξ ˙ 1 ∂ q t ∂ ξ ˙ 1 ∂ Z F ∂ ξ ˙ 1 ∂ q y ∂ x ˙ s ∂ x t ∂ x ˙ s ∂ y ∂ x ˙ s ∂ ξ 1 ∂ x ˙ s ∂ x s ∂ x ˙ s ∂ ξ 2 ∂ x ˙ s ∂ q t ∂ x ˙ s ∂ Z F ∂ x ˙ s ∂ q y ∂ ξ ˙ 2 ∂ x t ∂ ξ ˙ 2 ∂ y ∂ ξ ˙ 2 ∂ ξ 1 ∂ ξ ˙ 2 ∂ x s ∂ ξ ˙ 2 ∂ ξ 2 ∂ ξ ˙ 2 ∂ q t ∂ ξ ˙ 2 ∂ Z F ∂ ξ ˙ 2 ∂ q y ∂ q ˙ t ∂ x t ∂ q ˙ t ∂ y ∂ q ˙ t ∂ ξ 1 ∂ q ˙ t ∂ x s ∂ q ˙ t ∂ ξ 2 ∂ q ˙ t ∂ q t ∂ q ˙ t ∂ Z F ∂ q ˙ t ∂ q y ∂ Z ˙ F ∂ x t ∂ Z ˙ F ∂ y ∂ Z ˙ F ∂ ξ 1 ∂ Z ˙ F ∂ x s ∂ Z ˙ F ∂ ξ 2 ∂ Z ˙ F ∂ q t ∂ Z ˙ F ∂ Z F ∂ Z ˙ F ∂ q y ∂ q ˙ y ∂ x t ∂ q ˙ y ∂ y ∂ q ˙ y ∂ ξ 1 ∂ q ˙ y ∂ x s ∂ q ˙ y ∂ ξ 2 ∂ q ˙ y ∂ q t ∂ q ˙ y ∂ Z F ∂ q ˙ y ∂ q y , (19) where ∂ x ˙ t ∂ x t = − D a − e y − e h e q x e q h + e x T a , ∂ x ˙ t ∂ y = − e h e q y e q h + e y T a , ∂ x ˙ t ∂ ξ 1 = − K a T a , ∂ x ˙ t ∂ x s = D a T a , ∂ x ˙ t ∂ ξ 2 = 0 , ∂ x ˙ t ∂ q t = e h e q h T a , ∂ x ˙ t ∂ Z F = 0 , ∂ x ˙ t ∂ q y = 0 , ∂ y ˙ ∂ x t = e h e q x K p + e q h D a + e g − e x K p − e q h K i T a e q h T a , ∂ y ˙ ∂ y = e h e q y K p − e q h e y K p e q h T a , ∂ y ˙ ∂ ξ 1 = K a K p T a , ∂ y ˙ ∂ x s = − D a K p T a , ∂ y ˙ ∂ ξ 2 = 0 , ∂ y ˙ ∂ q t = − e h K p e q h T a , ∂ y ˙ ∂ Z F = 0 , ∂ y ˙ ∂ q y = 0 , ∂ ξ ˙ 1 ∂ x t = 1 , ∂ ξ ˙ 1 ∂ y = 0 , ∂ ξ ˙ 1 ∂ ξ 1 = 0 , ∂ ξ ˙ 1 ∂ x s = − 1 , ∂ ξ ˙ 1 ∂ ξ 2 = 0 , ∂ ξ ˙ 1 ∂ q t = 0 , ∂ ξ ˙ 1 ∂ Z F = 0 , ∂ ξ ˙ 1 ∂ q y = 0 , ∂ x ˙ s ∂ x t = B D a T s , ∂ x ˙ s ∂ y = 0 , ∂ x ˙ s ∂ ξ 1 = B K a T s , ∂ x ˙ s ∂ x s = B D a − D s T s , ∂ x ˙ s ∂ ξ 2 = − 1 R g T g T s , ∂ x ˙ s ∂ q t = 0 , ∂ x ˙ s ∂ Z F = 0 , ∂ x ˙ s ∂ q y = 0 , ∂ ξ ˙ 2 ∂ x t = 0 , ∂ ξ ˙ 2 ∂ y = 0 , ∂ ξ ˙ 2 ∂ ξ 1 = 0 , ∂ ξ ˙ 2 ∂ x s = 1 , ∂ ξ ˙ 2 ∂ ξ 2 = − 1 T g , ∂ ξ ˙ 2 ∂ q t = 0 , ∂ ξ ˙ 2 ∂ Z F = 0 , ∂ ξ ˙ 2 ∂ q y = 0 , ∂ q ˙ t ∂ x t = e q x e q h T w t , ∂ q ˙ t ∂ y = e q y e q h T w t , ∂ q ˙ t ∂ ξ 1 = 0 , ∂ q ˙ t ∂ x s = 0 , ∂ q ˙ t ∂ ξ 2 = 0 , ∂ q ˙ t ∂ q t = 1 T w t − 1 e q h − 2 h t H 0 , ∂ q ˙ t ∂ Z F = − 1 T w t , ∂ q ˙ t ∂ q y = 0 , ∂ Z ˙ F ∂ x t = 0 , ∂ Z ˙ F ∂ y = 0 , ∂ Z ˙ F ∂ ξ 1 = 0 , ∂ Z ˙ F ∂ x s = 0 , ∂ Z ˙ F ∂ ξ 2 = 0 , ∂ Z ˙ F ∂ q t = Q y 0 F H 0 , ∂ Z ˙ F ∂ Z F = 0 , ∂ Z ˙ F ∂ q y = − Q y 0 F H 0 , ∂ q ˙ y ∂ x t = 0 , ∂ q ˙ y ∂ y = 0 , ∂ q ˙ y ∂ ξ 1 = 0 , ∂ q ˙ y ∂ x s = 0 , ∂ q ˙ y ∂ ξ 2 = 0 , ∂ q ˙ y ∂ q t = 0 , ∂ q ˙ y ∂ Z F = c g H 0 W c g H 0 T w y W + q y Q y 0 V x λ ⁡ cot ⁡ α , ∂ q ˙ y ∂ q y = − g Q y 0 V x λ c W − 2 h y q y + H 0 Z F − q y Q y 0 λ cot ⁡ α c g H 0 T w y W + q y Q y 0 V x λ ⁡ cot ⁡ α 2 + g − 2 c h y W − Q y 0 λ c g H 0 T w y W + q y Q y 0 V x λ ⁡ cot ⁡ α .

The characteristic polynomial of the Jacobian matrix J(μ) is x 8 + a 1 x 7 + a 2 x 6 + a 3 x 5 + a 4 x 4 + a 5 x 3 + a 6 x 2 + a 7 x + a 8 = 0 . (20)

Based on the aforementioned stability analysis criterion, the stability domain can be drawn on the K _p –K _i plane. Parameter values of the HTRS and PG are T _a = 10 s, T _wt = 3.5 s, F = 1,500 m², T _wy = 2.0 s, W = 7.5 m, m _g = −0.1, e _g = 0, e _x = −1, e _y = 1, e _h = 1.5, h _t = 1.46 m, H ₀ = 80 m, e _qx = 0, e _qy = 1, e _qh = 0.5, h _y = 1.12 m, Q _y0 = 500 m³/s, H _x = 18.3 m, α = 0.0599282, D _a = 0.073, K _a = 2, T _s = 40, D _s = 0.4, R _g = 0.2, T _g = 40, and B = 0.1, and x denotes the eigenvalue of the polynomial J(μ) with value 3.

The bifurcation line is a crucial stability indicator comprising all Hopf bifurcation points on the K _p –K _i plane, which divides the whole parameter plane into stable and unstable domains. Thus, the location of bifurcation line determines stability margin in the parameter plane of the coupled system that also reflects dynamic characteristics.

With respect to the non-linear coupled system researched in this paper, K _i is chosen as the bifurcation parameter. Then, the stable domain and bifurcation line of the coupled system can be determined by solving Eq. 12, and they are shown in Figure 3A. In Figure 3A, three state points s ₁, s ₂, and s ₃ are chosen to investigate dynamic response of the coupled system under different parameter values. Coordinate values of the three selected state points and respective theoretical states of the dynamic response are shown in Table 1. In addition, from Figure 3B, it can be concluded that σ'(μ _c ) > 0, which indicates that the Hopf bifurcation is supercritical.

According to chaos theory and Hopf bifurcation theory, the discriminant conditions for Hopf bifurcation to occur in the system are Eqs 14–17. With the aim of verifying the correctness of the conclusions in Table 1, numerical simulation experiments are utilized in this section. Three state points s ₁, s ₂, and s ₃ are substituted into the state equation to solve the dynamic response of system state variables. Dynamic characteristics and the phase space trajectory of the state variables for the three corresponding state points are shown in Figure 4. The equilibrium point of the system is often chosen as the initial state of the system when simulation experiments of non-linear systems are carried out. Therefore, all the parameters of the system are substituted into the system of state equations and the equilibrium point of the state variables of the system is solved, and the state of the equilibrium point is obtained to be (−0.1335, 0). So, (−0.1335, 0) is chosen to be the initial state of the system.

The conclusions that can be obtained from Figure 4 are as follows:

(1) Figure 4A shows that Hopf bifurcation occurs at s ₁, which is coherent with the aforementioned theoretical analysis. From the phase space trajectory in Figure 4A, it can be concluded that at point s ₁, the system oscillates from the equilibrium point and quickly enters limit cycles. Furthermore, when the load disturbance of the coupled system is m _g = −0.1, state variables x _s and x _t oscillate with equal amplitude from the equilibrium point.

The highly complex dynamic behavior and significant impact on the stability region are consequences of the complex non-linear characteristics inherent to the HTRS. In order to reveal the coupling effect between the HTRS and PG, the coupled system was partitioned into two subsystems for non-linear dynamics analysis. This type of analysis enables a deeper understanding of the dynamic characteristics of the HTRS and PG, revealing the interrelationship between them and allowing for a more precise evaluation of system stability.

Subsystem 1 is the hydropower station considering a downstream surge chamber and sloping roof tailrace tunnel, while subsystem 2 is the PG. For subsystem 1, since the downstream surge chamber and sloping roof tailrace tunnel can decrease water hammer pressure during the transition process of the hydropower station, the influence of the downstream surge chamber and sloping roof tailrace tunnel on the operation of the hydropower station has to be considered. Although the tailrace tunnel can reduce the loss of outlet kinetic energy by using a turbine runner, water level fluctuation of the surge chamber and flow inertia variation of the tailrace tunnel will also interfere with the system. Moreover, the disturbance interaction not only directly affects stability of the hydropower station but also generates complex dynamic behavior during transients. Due to the fact that subsystem 2 is directly coupled to the hydropower station, changes of PG load and the connection and exit of power will affect safe operation of the coupled system. Hence, it is essential to study subsystem 1 and subsystem 2 on the dynamic behavior and stability of the coupled system.

The Jacobi matrix of the coupled system is obtained by differentiating state variables according to Eq. 12. According to matrix theory and Hopf bifurcation theory, a pair of purely imaginary eigenvalues will appear when Hopf bifurcation occurs, so this pair of eigenvalues is selected to analyze the stability of the coupled system, as shown in Figure 8.

Based on Figure 8, it can be concluded that two eigenvalues have the same change rule with controller parameters. When K _i is constant, with the increase of K _p , the real part of eigenvalue increases and the imaginary part decreases. The coupled system is unstable when the real part of eigenvalues is greater than zero, so it gradually transitions from the stable region to unstable region as K _p increases. Furthermore, as K _i increases, the area where the real part of the eigenvalue is greater than zero gradually increases, indicating that the instability region gradually increases. Therefore, by reasonably adjusting the controller parameters, the stability of the coupled system can be effectively improved.

The sensitivity of state variables is evaluated using standard deviation to verify the sensitivity of the coupled system to variations of HTRS parameters. To facilitate the analysis of the sensitivity of each parameter to the system state variables, the standard deviation of each parameter is quantified. s p = ∑ i = 1 n x i − x ¯ 2 n , (21) σ p = s p ∑ p = 1 m s p , (22) where x i = ∂ χ j ∂ ν k , χ = (q _y , Z _F , y, q _t , x _s , x _t ), and ν = (T _a , F, T _wt , T _wy ). x ¯ denotes the mean value of x _i for all data from i = 1 to i = n. s _p is the standard deviation of χ induced by HTRS parameters. σ _p is the ratio of each standard deviation. The sensitivity analysis was performed as in Figure 9, and the standard deviation was processed to obtain the ratio of s.d. in order to facilitate the analysis and comparison of the effect of each parameter on the system. Specific analysis results are shown in Figures 9A, 10.

Figure 9A describes the sensitivity analysis to changes of HTRS parameters at m _g = −0.1, K _p = 1, and K _i = 0.1. It can be seen from Figure 9A that the highest standard deviations of q _y and Z _F are correlated with T _wy while the highest standard deviations of y, q _t , x _s , and x _t are correlated with T _wt and T _a .

The sensitivity of HTRS parameters to state variables is shown in Figure 10. From Figure 10, it can be concluded that the state variables q _y and Z _F are more influenced by T _wy while effects of T _wt and T _a on the state variables x _s and x _t are more significant. Thus, q _y and Z _F are most sensitive to T _wy , but x _s and x _t are most sensitive to T _wt and T _a .

This section studies the sensitivity of PG parameters to state variables. The specific research method and normalization are the same as in the previous section, and the system operates under load condition, i.e., m _g = −0.1, K _p = 1, and K _i = 0.1. Sensitivity analysis to changes of PG parameters is shown in Figure 9B.

According to Figure 9B, state variables of the coupled system are most sensitive to parameters B and R _g , while the uncertainty sources of other parameters are less important. The results show that state variables are most sensitive to variability of B and R _g .

The sensitivity of PG parameters to state variables is shown in Figure 11. Since the sensitivity of PG parameters to state variables is not in the same order of magnitude, it is divided into two sub-graphs for better visualization. From Figure 11, it can be seen that q _y , x _s , x _t , and Z _F exhibited most sensitivity to the variation of B and R _g , while other parameters are less sensitive.

From Figures 10, 11, it can be concluded that the overall effect of the variation of HTRS parameters on x _t is greater than that of the variation of PG parameters. On the contrary, the overall effect of variable PG parameters on q _y , x _s , and Z _F is greater than that of variable HTRS parameters.