1 Introduction

The hydraulic turbine governing system (HTGS) is a complex, nonlinear, multivariable, time-varying and non-minimum phase system^[1-2], whose parameters change significantly as operating point changes, especially in large perturbation. An approximate linearization^[3-4] model is often employed to describe the HTGS in small perturbation, while some nonlinear factors, such as dead-zone, output saturation and clearance in feedback system of the servomotor, would have obviously effect on HTGS^[5-7] and should be taken into consideration according to research needs. Taking the output saturation of distributing valve for example, the position of main distributing valve may reach saturation when the perturbation approaches 10% of total output, which in turn affects dynamic behaviors of HTGS. Bifurcation analysis has been used to evaluate how the control system alters the steady-state dynamics in Refs.^[8-10]. In addition, existence of such bifurcation phenomena and other rich dynamics, such as synchronization, symmetric bifurcation and chaos that appear in many nonlinear dynamical systems^[11-12]are also studied in detail based on nonlinear models^[13-15], and have made a lot of meaningful research results^[16-17].

Bifurcation phenomena in HTGS have been considered as a cause of decreased availability and reduced stability^[18-19] for a long time. These studies are instrumental in identifying and documenting the oscillations in hydraulic generating units, however, the reports are scarce on the possible existence, the direction of Hopf bifurcation in HTGS and its effects on the oscillatory behaviors. In fact, with the head and capacity increasing of unit set, self-excited oscillations and unstable oscillation problems, caused by hydraulic, mechanical, electrical and other factors, have become prominent increasingly. Nonlinear dynamic analysis and stability control have been one of the issues that people are most concerned about in the hydropower station operation already^[20-22].

Motivated by the above discussions, the existence of Hopf bifurcation and its direction for a PID controller type speed governing system with saturation nonlinearity are investigated by using a small perturbation model in this paper. Simulations on both sides of bifurcation are presented by numerical calculations. Theory analysis and simulation results suggest the PID parameters located in limit cycle part should be avoided when a PID or a gain-scheduling PID scheme is used in the system. Otherwise, Hopf bifurcation or undesirable oscillation that is detrimental for the safe and stable operation of unit set may occur in the governing system.

2 Nonlinear Mathematical Model of HTGS

The closed-loop system is often regarded as a serial one with a linear part and a nonlinear part^[23-24]when studying the systems with saturation nonlinearity. A PID governing system with saturation is employed in this paper, as shown in Fig. 1, where x_r is the reference input, x and y are output of the speed and position of the servomotor, respectively, and m_r is the load perturbation. Transfer function of the PID type governor is described as follows:

${{G}_{PID}}\left( s \right)={{K}_{p}}+({{K}_{i}}/s)+{{K}_{d}}s$

(1)

where K_p, K_i and K_d are the proportional, integral and differential gains, respectively. An approximate linearization model G_t(s) is adopted for the hydro turbine in small perturbation in this paper, that:

${{G}_{t}}\left( s \right)={{e}_{y}}(1-e{{T}_{w}}s)/(1+{{e}_{qh}}{{T}_{w}}s)$

(2)

where G_t(s) denotes the transfer function from servomotor of the guide vane to hydro turbine moment in which a rigid water column model and a linear model of piping systems are employed^[25-26], e=e_qye_h/e_y-e_qh.

The dynamic equation for turbine output torque m_t can be described as:

${{{\dot{m}}}_{t}}=\frac{1}{{{e}_{qh}}{{T}_{w}}}\left[ -{{m}_{t}}+{{e}_{y}}y-\frac{e{{e}_{y}}}{{{T}_{w}}{{T}_{y}}}\left( u-y \right) \right]$

(3)

where e_y is a partial derivative of turbine torque with respect to the guide vane; e_qy is a partial derivative of the flow with respect to the guide vane; e_qh is a partial derivative of the flow with respect to the head; e_h is a partial derivative of turbine torque with respect to the head; T_w is water inertia time constant of a pressure guide-water system.

A synchronous generator G_g(s), connected to an infinite bus through a transmission line, is regarded as the target system. The second-order nonlinear model can be expressed by:

$~\left\{ \begin{align} & \delta ={{\omega }_{0}}\omega \\ & \dot{\omega }=({{m}_{t}}-{{m}_{e}}-K\omega )/{{T}_{ab}} \\ \end{align} \right.$

(4)

where δ, ω, K and T_ab denote the rotor angle, relative speed deviation, damping coefficient, sum of the machine starting time and load time constants, respectively. Electromagnetic torque of generator m_e is equal to its electromagnetic power P_e:

${{m}_{e}}~=\text{ }{{P}_{e}}$

Electromagnetic power P_e can be described with the following equation:

${{P}_{e}}=\frac{E{{\prime }_{q}}{{V}_{s}}}{x\prime dx\Sigma }\text{sin }\!\!~\!\!\text{ }\delta +\frac{{{V}_{s}}^{2}}{2}\frac{x\prime dx\Sigma -xqx\Sigma }{x\prime dx\Sigma xqx\Sigma }\text{sin}~2\delta $

(5)

Where the effects of speed deviations, damping coefficient andtorque variations are all included in the analysis of the generator dynamic characteristics^[27]. And

$\left\{ \begin{align} & x\prime dx\Sigma =x{{\prime }_{d}}+{{x}_{T}}+0.5{{x}_{L}} \\ & xqx\Sigma ={{x}_{q}}+{{x}_{T}}+0.5{{x}_{L}} \\ \end{align} \right.$

(6)

where E′_q and x′_qxΣ are the transient electric potential and synchronous reactance of q-axis, respectively; V_s is the infinite bus voltage and x′_dxΣ is the transient reactance of d-axis; x′_d and x_q are the direct axis transient reactance and quadrature axis reactance; x_T and x_L are the transformer short circuit reactance and the transmission line reactance, respectively.

Eqs.(4) -(6) are the simplified second-order generator model based on approximate linearization turbine model, which has been considered to be sufficient and widely applied in stability and dynamic quality analysis of HTGS from a power system operation view point. Higher order nonlinear generator model can be employed according to research needs.

Saturation in actuators is inevitable in many feedback control systems. Therefore, output saturation of the servomotor position in a HTGS will be taken into consideration in this paper. The nonlinear block of saturation can be expressed by u=N(z), N(z) represents the nonlinear function of the saturation part. To get the nonlinear equation of a closed-loop system, a new variable z is defined as follows:

$z={{K}_{p}}+({{K}_{i}}/s)+{{K}_{d}}s$

Then u=N(K_p +K_i /s+K_ds). The nonlinear block in Fig. 2 is given by:

$u=N\left( z \right)=({{y}_{m}}/2{{z}_{0}})(\left| z+{{z}_{0}} \right|-\left| z-{{z}_{0}} \right|)$

(7)

For convenience, set y_m=1 and z₀=1. Then Eq.(7) can be expressed by:

$\begin{align} & u=N\left( z \right)=\frac{1}{2}(|{{K}_{p}}+\frac{{{K}_{i}}}{s}+{{K}_{d}}s+1|- \\ & |{{K}_{p}}+\frac{{{K}_{i}}}{s}+{{K}_{d}}s-1|) \\ \end{align}$

(8)

Based on the discussions above, the four-dimensional differential equations that coupled each individual part of HTGS with saturation nonlinearity can be written as:

$\left\{ \begin{align} & \dot{\delta }={{\omega }_{0}}\omega \\ & \dot{\omega }=\frac{1}{{{T}_{ab}}}({{m}_{t}}-K\omega -\frac{E{{\prime }_{q}}{{V}_{s}}}{x\prime dx\Sigma }\text{sin}~\delta -\frac{{{V}_{s}}^{2}}{2}\frac{x\prime dx\Sigma -xqx\Sigma }{x\prime dx\Sigma xqx\Sigma }\text{sin}~2\delta ) \\ & {{{\dot{m}}}_{t}}=\frac{1}{{{e}_{qh}}{{T}_{w}}}(-{{m}_{t}}+{{e}_{y}}y-\frac{e{{e}_{y}}{{T}_{w}}}{{{T}_{y}}}\dot{y}) \\ & \begin{array}{*{35}{l}} \dot{y}=\frac{1}{{{T}_{y}}}(\frac{1}{2}(\left| {{K}_{p}}+\frac{{{K}_{i}}}{s}+{{K}_{d}}s+1 \right|- \\ \left| {{K}_{p}}+{{K}_{i}}s+{{K}_{d}}s-1 \right|)-y)\left( 9 \right) \\ \end{array} \\ \end{align} \right.$

(9)

Eq.(9) is the water-electro-mechanical coupled nonlinear model of HTGS that integrates turbine model of rigid water hammer, water pipes linear model^[28], PID type governor model and second-order generator model. The complex nonlinear nature problem within the system could be revealed much better by Eq.(9) compared to a linear one.

Numerical calculations of Eq.(9) can be carried out by using the Function ode45 in MATLAB to investigate the nonlinear dynamical behaviors of HTGS, e.g., Hopf bifurcation point, direction of Hopf bifurcation, bifurcation surface of PID parameters and time domain response waveforms of state variables. Time step is set to be 0.01 in the study. In addition, the Matcont software is also used to calculate the supercritical Hopf bifurcation points in the study.

3 The Existence of Dynamic Hopf Bifurcation and Its Direction 3.1 Existence Determination of Hopf Bifurcation

The existence conditions of the Hopf bifurcation for a four-dimensional nonlinear system, such as Eq.(9) above, should satisfy the following theorem according to the Hopf bifurcation theory of nonlinear dynamical systems. Meanwhile, the values collection of Hopf bifurcation can be also determined when Hopf occurs^[29].

Theorem: the Jacobi matrix characteristic polynomial for a nonlinear system y=F(x, v) at equilibrium point (x=0, x∈R⁴) can be defined as:

$f(\lambda ,\text{ }v)=~{{\lambda }^{4}}+{{p}_{3}}\left( v \right){{\lambda }^{3}}+{{p}_{2}}\left( v \right){{\lambda }^{2}}+{{p}_{1}}\left( v \right)\lambda +{{p}_{0}}\left( v \right)$

(10)

where v∈R¹ is the bifurcation parameter. The following conditions (a) and (b) hold if v=0:

(a) The coefficient p_i > 0 (i = 0, 1, 2, 3) , and

${{p}_{3}}{{p}_{2}}{{p}_{1}}={{p}_{3}}^{2}{{p}_{0}}+{{p}_{1}}^{2}$

(11)

$\left( b \right)\text{ }p{{\prime }_{0}}\ne (\frac{{{p}_{3}}}{{{p}_{0}}{{p}_{1}}}-\frac{{{p}_{1}}}{{{p}_{3}}})\frac{{{p}_{3}}p{{\prime }_{1}}-p{{\prime }_{3}}{{p}_{1}}}{{{p}_{3}}^{2}}+\frac{p{{\prime }_{2}}{{p}_{1}}}{{{p}_{3}}}$

(12)

where *′= d* /dv, p_i(i=1, 2, 3) is the coefficient of the characteristic equation. When | v | is sufficiently small, Hopf bifurcation exists in one side of v=0, and the period of limit cycle can be described as follows:

$T=2\pi \sqrt{({{p}_{3}}/{{p}_{1}})}$

(13)

Eqs.(10) -(13) are the direct conditions for existence of Hopf bifurcation for a four-dimensional nonlinear system. All the points on the surface that determined by Eq.(11) are Hopf bifurcation critical points that would lead to the amplitude oscillation of the unit, as shown in Fig. 3. For example, the system will suffer a periodic oscillation phenomenon, namely Hopf bifurcation, if conditions (a) and (b) are satisfied, i.e., the PID parameters are on the critical surface of Fig. 3. Structure stability of the system will change dramatically, leading to the system structure instability and complex nonlinear oscillations when the PID parameters located in the unstable area.

PID parameters instability region that lies beneath the surface can guarantee the HTGS working in a safe and stable state; correspondingly, the instability region that called limit cycle is located in certain area above surface. PID parameters on the surface or in limit cycle should be avoided when a PID or a gain-scheduling PID scheme is used. Otherwise, Hopf bifurcation or undesirable oscillations may occur. Fig. 3 can be applied to investigate the relationship between the PID parameters and system stability.

Taking v_d=K_d-K_d* or v_p=K_p-K_p* as the bifurcation parameter, when |v| is sufficiently small, the sustained periodic oscillation phenomenon will take place according to the theorem. For example, let K_i*=1, K_p*=2, the bifurcation parameter K_d can be calculated to be K_d1*=0.313 and K_d2*=2.036, respectively, as shown in Fig. 4. K_d1* and K_d2* can be also obtained through the intersection values between the ordinate line (K_p=2) and the boundary curve (K_i=1) , as shown in Fig. 5. All the points that in limit cycle, for example, K_i=1, K_d=2, K_p=2.5 (v_p>0) , would lead to the undesirable diffused oscillation, as shown in Figs. 6 and 7. Furthermore, the stable and unstable PID parameters can be obtained much easier through Fig. 5 when compared with Fig. 3.

Fig. 4 shows the bifurcation diagram of speed ω calculated by the Matcont software. Two bifurcation points of K_d are presented in Fig. 4 when K_p=2, K_i=1, i.e. K_d1*=0.313 and K_d2*=2.036. Speed ω converge to equilibrium point 0 when K_d1＜K_d＜K_d2 and to stable limit cycle when K_d>K_d2 or K_d＜K_d1. What’s more, ω will converge to stable limit cycle when 0 ＜K_d＜K_d1* or K_d> K_d2*.

The stability of system can also be analyzed from perspective of mathematics. The properties and Jacobi matrix roots for different values of K_d are presented in Table 1. The Jacobi matrix has two conjugate complex roots with negative real parts when the system is in the stable domain (K_d1*＜K_d＜K_d2*), and then the system tends to be stable (equilibrium point); when the parameters are the critical bifurcation point (K_d1*, K_d2*), Jacobi matrix has two negative real roots and two pure imaginary conjugate complex roots, the system is in center of bifurcation; When K_d＜K_d1* or K_d>K_d2*, the matrix has two conjugate complex roots, a negative real part one and a positive real part one, and system will converge to stable limit cycle at this time.

Theoretical analysis of Table 1 can be considered as an explanation for the phenomenon that occurs in Fig. 4. Dynamic behaviors of nonlinear systems in critical bifurcation point can be seen from Fig. 4 and Table 1 clearly. Moreover, it can be obtained that when K_d is far away from two bifurcation points and near to middle place between K_d1* and K_d2*, the system is more stable and the convergence rate is much quicker. In practical applications, values of K_d should be between two critical places and the farther the better. Excessive and tedious calculations can be avoided by using the theorem above. However, the theorem does not give the direction of bifurcation; further study will be discussed in Section 3.2.

Point 1 (K_i=1, K_d=0.313, K_p=2.32, v_p>0) , point 2 (K_i=1, K_d=0.313, K_p=2.0, v_p=0) and point 3 (K_i=1, K_d=0.313, K_p=1.6, v_p＜0) are three typical points in Fig. 8, respectively. Direction of Hopf bifurcation can be told from these points, which will be verified by the simulation results in Figs. 6-7 and Figs. 9-12.

3.2 Direction of Hopf Bifurcation

Researches in Section 3.1 can be applied to calculate the existence conditions of Hopf bifurcation, but direction and stability of bifurcation need to be further studied. Eq.(10) will have two negative real roots and two pure imaginary conjugate complex roots if condition (a) holds:

${{\lambda }_{k}}=\pm i{{\omega }_{0}}=\pm \text{i}\sqrt{({{p}_{1}}/{{p}_{3}})},\text{ }k=1,\text{ }2$

(14)

Taking v_p=K_p-K_p* as bifurcation parameter; then the application of Eq.(10) yields:

$\frac{d{{\lambda }_{k}}}{dv}{{|}_{v=0}}=\frac{1.56{{\omega }_{0}}^{2}\mp 2.78i{{\omega }_{0}}}{(3{{p}_{3}}{{\omega }_{0}}^{2}-{{p}_{1}})\pm 2i{{\omega }_{0}}(2{{\omega }_{0}}^{2}-{{p}_{2}})}$

(15)

The real part of Eq.(15) is:

$\begin{align} & a\prime \left( {{v}_{0}} \right)=\text{Re}\frac{d{{\lambda }_{k}}}{dv}{{|}_{v=0}}= \\ & {{\omega }_{0}}^{2}\frac{1.56(3{{p}_{3}}{{\omega }_{0}}^{2}-{{p}_{1}})\pm 5.56{{\omega }_{0}}^{2}(2{{\omega }_{0}}^{2}-{{p}_{2}})}{{{(3{{p}_{3}}{{\omega }_{0}}^{2}-{{p}_{1}})}^{2}}+4{{\omega }_{0}}^{2}{{(2{{\omega }_{0}}^{2}-{{p}_{2}})}^{2}}}, \\ & k=1,\text{ }2\text{ } \\ \end{align}$

(16)

Where the denominator and numerator of the right side in Eq.(16) are positive values according to condition (a) and Eq.(17) , i.e., Re (dλ_k/dv) |_v=0>0, (k=1, 2) , the Hopf bifurcation points that presented in Table 1 is supercritical. The equilibrium point of the system is a stable focus if v＜0; Correspondingly, a periodic movement will be generated from the equilibrium point when v>0 (taking point 1 v_p of Fig. 8 for example), in this case, a stable limit cycle will appear in the system and a diffused sustained oscillation of output can be observed, as shown in Figs. 6 and 7.

The corresponding sustained oscillation curves of δ, ω and the system stable focuses at point 2, where v_p=0, are presented in Figs. 9 and 10, respectively, in which the oscillation period is T≈1.053 s, that is consistent with the calculation result from Eq.(13) . The governing system makes a periodic oscillation instead of converging to equilibrium point 0 and further reveals that the system is not able to be stable.

The curves for δ, ω and stable focus at point 3, where v_p＜0, are shown in Figs. 11 and 12. The system state variables will finally converge to the equilibrium point 0. Moreover, the system is more likely to converge to an equilibrium point 0 quickly by making the PID parameter apart from the Hopf bifurcation points, the further the better.

Analysis and simulation above show that a supercritical Hopf bifurcation exists in HTGS. Figs. 6-7 and Figs. 9-12 illustrate the direction of Hopf bifurcation and its characteristic: Hopf bifurcation is supercritical when the real part Re(dλ_k/dv)|_v=0>0, (i=1, 2) , and the equilibrium point of the system is a stable focus if v＜0, and when v>0 a periodic movement will occur.

4 Effect of Different Parameters on Stability

Curves in Fig. 5 represent K_p=K_p* for different values of K_d when K_i=0, 1, 3, respectively. A stable limit cycle will appear in HTGS according to conditions (a) and (b). When K_p and K_d locate in the limit cycle part of Fig. 5, a sustained oscillation with certain period will appear in the system, which is not generally expected in the flied operation. Consequently, suitable parameters should be selected carefully so that the PID parameters can be located in stable region in Fig. 5 to guarantee the unit set operating in a stable state. Similarly, curves in Figs. 13-16 represent K_p* for different values of K_d when T_w, T_ab, T_y and K changes, respectively.

The stability regions have decreased radically when T_w and T_y increase, indicating that T_w and T_y are the detrimental factors that lead to instability of system, consistent with the conclusions of other scholars. Meanwhile, stability regions expand with T_ab and K increasing, which can help the system converge to the equilibrium point quickly.

5 Example and Numerical Simulations

Parameters of the HTGS that works under certain operating conditions are setting as follows: ω₀=314, T_w=0.8 s, T_ab=9.0 s, T_y=0.1 s, K=2, E′_q=1.35, V_s=1, x′_dxΣ=1.15, x′_qxΣ=1.474, e_qh =0.5, e_y =1.0, e=0.7, respectively. Jacobi matrix characteristic polynomial coefficients p_i of the system Jacobi matrix characteristic polynomial coefficients p_i at equilibrium point (0, 0, 0, 0) are obtained as follows:

$\left\{ \begin{array}{*{35}{l}} {{p}_{3}}=12.72-1.56*{{K}_{d}} \\ {{p}_{2}}=2.78*{{K}_{d}}-1.56*{{K}_{p}}+62.07 \\ {{p}_{1}}=2.78*{{K}_{p}}-1.56*{{K}_{i}}+434.15 \\ {{p}_{0}}=2.78*{{K}_{i}}+857.2 \\ \end{array} \right.$

(17)

Three-parameters implicit polynomial is obtained after substituting Eq.(17) into Eq.(11) and a space bifurcation surface can be carried out, as shown in Fig. 3. Bifurcation diagram for ω is illustrated in Fig. 4. Curve of bifurcation points is shown in Fig. 8, in which three typical points for direction of Hopf bifurcation are included. The limit cycles for the three typical points of Fig. 8 are shown in Figs. 7, 10 and 12, and the corresponding oscillation curves of δ and ω are shown in Figs. 6, 9 and 11. Figs. 13-16 in part 4 present the curves that represent K_p^* for different values of K_d when K_i=1, where the system parameters T_w, T_ab, T_y, K change, respectively.

6 Conclusions

A fully coupled nonlinear mathematical model for a PID governing system with saturation has been developed in this paper. Extensive investigations of the Hopf bifurcation behaviors of an example HTGS under various conditions were conducted using the proposed model. Unexpected limit cycle oscillation in HTGS was clearly illustrated from various perspectives when PID parameters meet certain conditions. The study also revealed that stability of the system was closely related to system parameters such as T_w, T_ab, indicating the parameters should be carefully selected to guarantee safe and stable operation of hydro turbine unit. The work in this paper suggests that Hopf bifurcation can be considered as one of the possible reasons for sustained oscillations phenomenon commonly seen in HTGS. Simulation results validate the theoretical analysis. This paper presents a novel analysis approach to research the stability of HTGS.