Abstract
A nonlinear mathematical model for hydroturbine governing system (HTGS) has been proposed. All essential components of HTGS, that is, conduit system, turbine, generator, and hydraulic servo system, are considered in the model. Using the proposed model, the existence and stability of Hopf bifurcation of an example HTGS are investigated. In addition, chaotic characteristics of the system with different system parameters are studied extensively and presented in the form of bifurcation diagrams, time waveforms, phase space trajectories, Lyapunov exponent, chaotic attractors, and Poincare maps. Good correlation can be found between the model predictions and theoretical analysis. The simulation results provide a reasonable explanation for the sustained oscillation phenomenon commonly seen in operation of hydroelectric generating set.
1. Introduction
Many hydropower plants have been built [1–3] worldwide to harness the energy of falling or running water for electricity purpose. An important part of hydropower plant is the hydraulic turbine governing system (HTGS), which serves to maintain safe, stable, and economical operation of hydropower generating unit [4]. The HTGS is in nature a complex nonlinear, multivariable, time-varying, and nonminimum phase system, which involves the interactions between hydraulic system, mechanical system, and electrical system [5]. The complex dynamic behaviors of the HTGS significantly influence the operation conditions of hydroelectric generating set. For instance, oscillatory problems in hydroelectric generating units were reported to be closely related to the possible Hopf bifurcation and chaotic oscillatory behaviors in HTGS [6–9]. In the absence of a model to conveniently predict the dynamic behaviors of the system, means to address the practical operational issues of hydroelectric generating units, that is, sustained oscillation phenomenon [6, 10], have to be limited.
The literature review reveals that considerable research efforts have been devoted to the modeling of each individual part of the HTGS [11]. For instance, elastic model and nonelastic model [12, 13] for the conduit system have been adopted in long pipeline and short pipeline systems, respectively. In addition, various hydroturbine models were also available in [10, 14–26]. While Sanathanan [14] claimed that the output of hydroturbine was proportional to hydraulic head and volume flow, Hannett et al. [15] approximated this term by its first-order Taylor formula and Kishor et al. [16, 17] proposed six detailed expressions about these six transfer coefficients with turbine speed and head, which were shown to provide reasonably accurate predictions for HTGS with small disturbance. X. Liu and C. Liu [10] studied small disturbance stability of hydropower plant with complex conduit with a linear turbine model, but they fell short of including elastic water-hammer effects. Bakka et al. made significant contributions to modeling, simulation [18–22], and control [23–26] of wind turbine system, which stimulates modeling of hydroturbine in this paper.
In regards to the power generator, models with different orders could also be found in [27–29]. While higher order power generator models typically represent better accuracy, the overall computational efficiency would be greatly sacrificed especially in case of the modeling of a fully coupled governing system. To compromise with the computational efficiency, the power generator model has to be selected carefully such that it provides accurate approximation with acceptable computational effort.
Aside from the abovementioned models for individual part of the HTGS, investigations on bifurcation and chaos for the HTGS at the system level were reported only in a few studies [6, 30, 31]. For instance, Konidaris and Tegopoulos [31] investigated the oscillatory problems in hydraulic generating units; Mansoor et al. [6] successfully reproduced an oscillatory phenomenon whose causes were difficult to be identified with limited recorded data. He also proposed a methodology to improve the stability of the control system. While these studies were instrumental in identifying and documenting the oscillations in hydraulic generating units, neither of them noticed the effects of Hopf bifurcation on the oscillatory behaviors. Ling and Tao [30] first reported the influence of bifurcation phenomenon on the sustained oscillations in his study of Hopf bifurcation behavior of HTGS with saturation. However, his model was oversimplified by employing first-order model for the power generator and the PI governing system. Chen et al. [32] developed a novel nonlinear dynamical model for hydroturbine governing system with a surge tank and studied exhaustively the influences of different parameters for the first time. However, they failed to include adequate theoretical analysis and description on bifurcation. Determination of the existence and detailed calculation of Hopf bifurcation were also not considered in their work.
The above literature review shows that numerous well-developed models for each individual part of HTGS are available, but investigations of bifurcation and chaotic oscillations in HTGS with a fully coupled model of the system are rarely seen. As such, this paper aims to develop a fully coupled nonlinear dynamical model for HTGS and investigate the bifurcation and chaotic oscillatory behaviors of the HTGS. In addition, a theorem for existence determination of Hopf bifurcation for four-dimensional nonlinear system has been proposed for convenient prediction of the bifurcation of critical points, which would otherwise be impractical especially for high-dimensional systems due to tremendous computational demand by conventional analysis methods, that is, Lyapunov-Schmidt (L-S) method, center manifold method [33, 34], or normal form theory [35].
There are three main contributions of this paper compared with prior works. First, a new four-dimensional fully coupled nonlinear mathematical model of HTGS was presented and the parameters were from a practical power station, which made the work more consistent with actual project compared with [10, 30] work. Second, the theorem for stability, Hopf bifurcation, and dynamic quality analysis of four-dimensional system that can avoid excessive and tedious calculations was firstly introduced in the paper, providing a new approach for HTGS analysis and computation. Third, nonlinear dynamical behaviors of the above system with different parameters were studied in detail and necessary numerical simulation results were presented.
The paper is outlined as follows. First, the formulations of a fully coupled nonlinear dynamical model will be presented. Next, the investigations of the Hopf bifurcation and chaotic behaviors of the system will be described. Finally, a brief conclusion will be given.
2. Nonlinear Mathematical Model of HTGS
HTGS consists of five parts, that is, conduit system, hydroturbine, governor, electrohydraulic servo system, and power generator. Model for each individual part has been well developed. Water from reservoir enters tunnel first and then flows through penstock before reaching turbine gate. Next, it flows into scroll casing to promote the hydroturbine to rotate. The power generator and hydroturbine are connected by a shaft coupling. Water that flows into the hydroturbine can be regulated by wicket gates, which are controlled by the governor system. The governor system operates accordingly given the deviation between electric demand and developed torque [35].
2.1. Conduit System Model
A no-elastic model [28] is employed for the conduit system in this study. The unsteady flow partial differential equations in pressure pipes can be described as where , , are parameters of the penstock. They denote the diameter, head loss, and area of the pipeline, respectively. , are hydraulic head and turbine flow in penstock in operating condition, is pressure wave velocity, and is the length from upstream.
The head and flow equation between two sections of penstock can be deduced from (1). It can be described as [28]
Subscripts , are symbols of upstream and downstream section of pipeline, respectively. and are the composite equations of parameters of the penstock, , .
, , can be written as
Providing that the head loss is negligible, and can be rewritten as follows:
With the hydraulic friction losses being trivial and (tunnel connects with reservoir directly), the head and flow function is simplified as follows:
It can be seen from (5) that the water hammer transfer function is a nonlinear hyperbolic tangent function, which was inconvenient to use and should be expanded by series. Substituting (4) into (5), the equation for the water strike transfer function at point is obtained as
Typically, there are the two models, , as follows: where , . Equations (6) to (8) represent three kinds of water hammer models: the first two models are called elastic water hammer model and the last one is rigid water hammer model that is employed in the paper. In the above equations, the pipeline is assumed to have constant cross-sectional area over the full length, which is impossible in actual engineering. Relevant parameters in the field are obtained by the following formulas: where , , denote the length, cross-sectional area, and the velocity of pipeline , respectively, and there are pipes in total.
2.2. Turbine Model of Rigid Water Hammer
For small perturbation around the rated operating point, the equation of the turbine can be represented as below:
The six constants of hydroturbine , , , , , are the partial derivatives of the torque and flow with respect to turbine speed, guide vane, and head, respectively. These constants may vary as the operating point changes.
In the dynamic models of the turbine and conduit system [12, 14] where the relative deviation is used to represent the state variables, the relationship between turbine torque and its output power is
As the unit speed changes little, the speed deviation , leading to . Then the transfer function of the turbine and conduit system is where is transfer function of conduit system defined in (6).
In case of rigid water hammer, the coefficient in (6) is 0 and the dynamic equation for the turbine output torque, as shown in Figure 1, becomes
2.3. Generator Model
A synchronous generator [36], which connected to an infinite bus through a transmission line, is considered as the target system. The second-order nonlinear dynamical model, after making standard considerations, can be written as where , , , and denote the rotor angle, relative speed deviation, damping coefficient, and mechanical starting time, respectively. The electromagnetic torque of the generator is equal to its electromagnetic power :
The electromagnetic power can be calculated with the following formula: where the effects of speed deviations, damping coefficient, and torque variations are all included in the analysis of the generator dynamic characteristics
Equations (14) to (17) are the simplified second-order nonlinear generator model based on the turbine model of rigid water hammer, which has been widely applied in nonlinear controller design and stability analysis of HTGS. It is often used for the stability characteristics and dynamic quality analysis of HTGS from the perspective of power system. Higher order nonlinear generator model could be employed according to research needs.
2.4. Hydraulic Servo System Model
The servomotor, which acts as the actuator, is used to amplify the control signals and provide power to operate the guide vane. Its transfer function can be written as where is the engager relay time constant.
2.5. Governor Model
At present, a parallel PID controller [18, 20, 22–26] is widely used in hydraulic turbine governors [27, 30] in filed. Its transfer function is given as
Substituting (19) into (18), the results can be written in state-space form. At last, the mixed function can be obtained as follows:
Based on the discussions above, the differential equations that coupled each individual part of the turbine nonlinear control system can be written as
Equation (21) is the four-dimensional water-electromechanical coupled model of HTGS that integrates the turbine model of rigid water hammer, the water pipes linear model, and the nonlinear dynamic generator model. Compared to the linear model, it could reflect the complex nonlinear nature problem within the system much better. Equation (21) could be applied to analyze and simulate the dynamic characteristics of the HTGS.
At the equilibrium point , the following condition has to be satisfied:
Equations (21) and (22) can be numerically solved to investigate the nonlinear behaviors of HTGS such as Hopf bifurcation points, bifurcation surface of PID adjustment coefficients, time domain response waveforms of state variables, and Lyapunov exponent with MATLAB by using a variable-step continuous solver based on the four-order Runge-Kutta formula. Time step is 0.01 in the study.
3. The Existence of Dynamic Hopf Bifurcation
3.1. Existence Determination of Hopf Bifurcation
For a four-dimensional nonlinear system, the criteria for the existence of Hopf bifurcation are given in the following theorem. Meanwhile, the bifurcation value collection can also be determined when Hopf bifurcation occurs.
Theorem 1. For a nonlinear system , , is the bifurcation parameter, is the equilibrium point, and the Jacobi matrix characteristic polynomial at the equilibrium point is
If , the following conditions hold.
(1) The coefficient (),
(2)
where , when is sufficiently small, Hopf bifurcation exists in one side of , and the period of limit cycle can be described as follows:
This theorem is the direct algebraic criterion to determine the existence of Hopf bifurcation for a four-dimensional system, which can avoid excessive and tedious calculations. However, the theorem does not give bifurcation direction. Each point on the surface that is determined by (24) is bifurcation critical point when Hopf bifurcation occurs (the whole bifurcation critical points constitute the entire bifurcation collection). For example, periodic oscillation phenomenon, namely, Hopf bifurcation, will occur, if the system parameters and the PID parameters satisfy conditions 1 and 2 of the above theorem in a certain operating condition. The surface determines the scope of the system’s stability region also. In fact, from the perspective of nonlinear dynamical systems structure stability analysis, when the PID parameters can meet with some certain conditions, the structure stability of the system will change dramatically and then comes to the system structure instability, which leads to complex nonlinear oscillations of the system.
3.2. Analysis and Simulation for Hopf Bifurcation
Taking a practical power plant, for example, setting parameters , s, s, s, , , , , , , , and , respectively, analyzes theoretical calculation and simulation in the paper based on the data. The followings are the system Jacobi matrix characteristic polynomial coefficients at equilibrium point :
It requires all the characteristic polynomial coefficients () by condition 1 when using Theorem 1 to determine the existence of Hopf bifurcation. For example, ; then , which means that the Hopf bifurcation area is limited to the range of . Similarly, the ranges of and can be obtained. After substituting into (24), a three-parameter implicit polynomial is obtained, and a space bifurcation surface can be made according to the polynomial, as shown in Figure 2.
All the points on the surface are Hopf bifurcation critical points that would lead to the amplitude oscillation of the unit; stability region lies beneath the surface; PID parameters in stability region ensure the stable and safe operation of units; the instability region is located above surface in certain area. Values of the PID parameters that are above the surface or on the surface should be avoided when a PID or a gain-scheduling PID scheme is used in the system. Otherwise, the Hopf bifurcation or undesirable oscillations may occur in the governing system. Figure 2 can be used to study the relationship between the stability of system and PID parameters, and stable PID parameters could be chosen from it directly. When taking or as the bifurcation parameter, it can be calculated that the left hand side of (25) is equal to zero, while the right hand side is not; thus condition 2 of the theorem holds.
According to Theorem 1, when is sufficiently small, Hopf bifurcation phenomenon of the system will take place. For example, when and , the bifurcation parameters and can be calculated from (24). and can also be got by the intersection values between the curve () and the line (), as shown in Figure 3. The curve () is the boundary between stability region and limit cycle; point that located in limit cycle, for example, , , and , would lead to the undesirable oscillation. Figure 3 is more intuitive to tell the stale and unstable PID parameters compared with Figure 2.
Figure 4 is the bifurcation diagram of unit speed . When , , and is the bifurcation parameter, there are two bifurcation points in the figure; that is, and as presented calculated by the Matcont software. It is observed that the system is stable and system state variable will converge to the equilibrium point 0 when ; it will converge to a stable limit cycle when or , and the amplitude of the limit cycle oscillation corresponds to the ordinate value of limit cycle curve. The bifurcation thresholds and in Figure 4 are consistent with the intersection values in Figure 3 and the results calculated from (24) as shown in Figure 2. Table 1 shows the properties of the system and the changes of Jacobi matrix eigenvalues with different . When the system is in stable domain (), all eigenvalues of Jacobi matrix have negative real parts and the system tends to be stable. As the parameter changes to the critical bifurcation point , the eigenvalues of Jacobi matrix have complex conjugate eigenvalues with zero real part; that is, the system is in a supercritical state. When or , the matrix has two complex conjugate eigenvalues; one has a negative real part and the other has a positive real part, system state variables cannot converge to the equilibrium point at this time, and the system is in stable amplitude oscillation (limit cycle) in this region.
Table 1 provides theoretical explanation for the phenomenon that occurs in Figure 4. Dynamic behaviors of nonlinear systems in critical bifurcation point on both sides can be seen from Figure 4 and Table 1 clearly. Furthermore, it can be obtained from simulation and theoretical analysis that the farther from the two bifurcation points ( and ) and the nearer to middle place the is, the quicker the convergence rate is and the more stable the system is. In practical applications, value of should be between two critical places and the farther the better, which will be verified by the simulation results in Figures 5, 6, 7, and 8.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Time domain response waveforms of and after stabilization can be seen in Figure 5 where the PID parameters are , , and , which are in limit cycle area in Figure 3, indicating that system state variables and will converge to a stable limit cycle. The motion of the system at this value of adjustment coefficients is sustained periodic oscillation and the oscillations amplitude of is the limit cycle curve ordinate value at point , as can be seen from Figure 4. In addition, the oscillation cycle of and at bifurcation point is about 1.053 s. Meanwhile, it can be got that , easily from (27), so the period of limit cycle s can be calculated from (26) at the Hopf bifurcation point, which is consistent with the simulation result. These results indicate that the governing system makes a periodic vibration and further reveal that the system is not able to be stable. Similarly, the waveforms of and for the value of , which are beneath the limit cycle curve (), namely, the stable parameters, are presented in Figures 6 to 8. After a period of time, each state ( and ) will converge to equilibrium point 0 and the system keeps steady at these values of adjustment coefficients; the biggest and value are 0.5, 0.01, respectively. It is observed that when value of is farther from the two critical points ( and ) and nearer to middle place (1.175), the convergence rate is quicker and the system is more stable.
Figures 9 and 10 show the phase space orbits of the system variables. It can be observed that the closed orbits limit cycles for these parameters formed in a certain region in the phase space orbit when ; thus , , and . The system shows a diffused and nonlinear growth oscillations curve. Hopf bifurcation occurs and the system tends to be a stable oscillation state where , , , , indicating that the system is in sustained periodic oscillation state; this value of the differential adjustment coefficient cannot be used in practical applications. Therefore, when the turbine regulating system uses PID (or gain scheduling PID) strategy, the values of adjustment coefficients should be beneath the limit cycle curve. Otherwise, there may be a sustained and stable oscillations, instead of tending to be stable, and the governing system would lose stability finally, leading to an unstable control.
4. Analysis and Simulation for Chaos
When , with the continually increasing of , unit speed on bifurcation diagram remains as a single rising curve until . The bifurcation diagram gets an intricate pattern when , the system state parameters converge neither to the equilibrium point nor to a stable limit cycle, and it undergoes a random motion without any rules, namely, the chaotic motion. There are rich and complex nonlinear dynamical behaviors in this area, as shown in Figure 11.
Lyapunov exponent is a quantitative indicator to measure the system dynamic behavior, which indicates the average rate of convergence or divergence of the system in phase space among different adjacent tracks. The existence of chaotic dynamics for the system can be judged intuitively by the largest Lyapunov exponent depending on whether it is greater than 0 or not. The system chaotic motion can be obtained by Lyapunov exponent with different , as presented in Figure 12. When , the largest Lyapunov exponent LEl is zero, indicating that the system is in periodic motion state (limit cycle). In the vicinity of , LE1 changes between 0 and 1, and thus, the system movement state alternates between periodic motion and chaotic motion; LE1 and LE2 are both greater than zero till and means chaotic motion that leads to the unstable state of the control system exists.
Figure 13 shows the system stability domain, limit cycles, and chaotic region of the system when , There are three different areas in the Fig, namely, the stable region, limit cycle, and chaos region, respectively; boundary between periodic motion and chaotic zone is determined by the three PID parameters especially when the maximum Lyapunov exponent LEl is greater than zero.
Chaotic motion of nonlinear systems can generate attractors in the phase space with unique nature. The original stable periodic motion turns to instability after entering into the chaotic region. Chaotic attractor has a complex motion internal with countless unstable periodic orbits studded in it. Overall, the system trajectories are always stretching and varying within a certain range and quite disordered in some part, leading to the chaotic strange attractor. When , , and , chaotic attractors that are similar to the one with scroll structure in Chua’s circuit [37] generate in the and phase space, as shown in Figures 14 and 15. It is observed that they have the characteristic of being globally bounded, but being local unstable, which would do great harm to the HTGS in operation, further reveals that the system would lose stability finally. Figure 16 is the time domain response waveforms of and , result shows that the system makes a disordered and aperiodic oscillation, and behaviors of the system are hardly to be predicted. The rotor angle and turbine speed have left their original position largely, oscillating between −50 and 50 −0.1 to 0.1, respectively, with a random period, and the system is in a chaos, which means that this value of the differential adjustment coefficient should be avoided.
(a)
(b)
Poincare maps are shown in Figures 17 and 18. It can be seen that numerous intersections of the continuous trajectory turn up in the Poincare maps on section of and . However, they neither fill with the entire phase space nor distribute discretely in the and phase space but occur in a certain region frequently, constituting a piece of dense point just like cloud. These results indicate that the system is in chaotic state and the system cannot tend to be stable. Figures 11, 12, 13, 14, 15, 16, 17, and 18 verify the chaotic and aperiodic movement of the system from mathematical perspective.
5. Conclusion
In this paper, a fully coupled nonlinear mathematical model for HTGS has been developed. Extensive investigations of the nonlinear behaviors of an example HTGS were conducted using the proposed model. The existence and stability of Hopf bifurcation were studied first. Both the simulation results and theoretical analysis predicted that Hopf bifurcation would occur when the governor’s parameters satisfy certain criteria. It was also found that the convergence rate and stability of the system were closely related to the bifurcation parameters, indicating that the PID adjustment coefficient parameters have to be carefully selected to guarantee the safe and stable operation of hydroelectric generating unit. Relevant stability ranges of PID adjustment coefficients for the example HTGS were presented in the paper. In addition to the study of Hopf bifurcation, chaotic behaviors of the system were also investigated extensively using different methods. It was found that in certain region of PID parameter space, chaos phenomenon would occur and lead to the chaotic oscillatory and unstable control of the system, which has been identified as a cause of decreased availability and reduced stability of the system. The above investigations on Hopf bifurcation and chaos phenomenon of the example HTGS clearly illustrated the complex nonlinear nature of the system. The simulation results based on the proposed model were shown to have good correlation with the theoretical analysis predictions. The model presented in this study for the analysis of the dynamic behaviors of HTGS can be further extended to a more sophisticated system.
Some related topics, such as data-driven framework [38], chaos control, and Hopf Bifurcation in HTGS with nonlinear saturation, would be considered in the future work to achieve more practical oriented results.
Nomenclature
| : | Laplace transform of , p.u. |
| : | Laplace transform of , p.u. |
| : | Initial turbine flow of turbine, m3/s |
| : | Initial hydraulic head of turbine, m |
| : | Gravitational acceleration |
| : | Penstock length, m |
| : | Penstock area, m2 |
| : | Penstock diameter, m |
| : | Head loss coefficient, p.u. |
| : | Pressure wave velocity, m/s |
| : | Elastic time constant, s |
| : | Water starting time, s |
| : | Incremental deviation of guide vane/wicket gate position, p.u. |
| : | Partial derivatives of turbine torque with respect to head, guide vane, and speed, p.u. |
| : | Partial derivatives of the flow with respect to head, guide vane, and turbine speed, p.u. |
| : | Rotor angle, rad |
| : | Turbine/rotor speed, rad/s |
| : | Base angular speed, rad/s |
| : | Speed deviation |
| : | Transient electric potential of -axis, p.u. |
| : | Infinite bus voltage, p.u. |
| : | Transient reactance of -axis, p.u. |
| : | Synchronous reactance of -axis, p.u. |
| : | Direct axis transient reactance, p.u. |
| : | Quadrature axis reactance, p.u. |
| : | Transformer short circuit reactance, p.u. |
| : | Transmission line reactance, p.u. |
| : | Mechanical starting time, s |
| : | Mechanical torque of turbine, N.m |
| : | Electrical torque, N.m |
| : | Damping factor, p.u. |
| : | Terminal active power |
| : | Output power |
| : | Engager relay time constant |
| : | Proportional adjustment coefficient |
| : | Integral adjustment coefficient |
| : | Differential adjustment coefficient. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.