# Nonlinear Dynamical Analysis of Hydraulic Turbine Governing Systems with Nonelastic Water Hammer Effect

Authors: Junyi Li; Qijuan Chen
Journal: Journal of Applied Mathematics
 (2014)

---

## 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.


---

## Body

## 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
(1)
Momentum equation:
∂
H
∂
x
+
1
g
A
∂
Q
∂
t
+
f
Q
2
2
g
D
A
2
=
0
,
Continuity equation:
∂
Q
∂
x
+
g
A
a
2
∂
H
∂
t
=
0
,
where D, f, A are parameters of the penstock. They denote the diameter, head loss, and area of the pipeline, respectively. H, Q are hydraulic head and turbine flow in penstock in operating condition, a is pressure wave velocity, and x 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]
(2)
[
H
A
(
s
)
Q
A
(
s
)
]
=
[
c
h
(
r
Δ
x
)
-
Z
c
s
h
(
r
Δ
x
)
-
s
h
(
r
Δ
x
)
Z
c
c
h
(
r
Δ
x
)
]
[
H
B
(
s
)
Q
B
(
s
)
]
.SubscriptsB, A are symbols of upstream and downstream section of pipeline, respectively. r and Z
c are the composite equations of parameters of the penstock, r
=
L
C
s
2
+
R
C
s, Z
c
=
r
/
C
s.L, R, C can be written as
(3)
L
=
Q
0
g
A
H
0
,
C
=
g
A
H
0
(
a
2
Q
0
)
,
R
=
(
f
Q
0
2
)
(
g
D
A
2
H
0
)
.Providing that the head loss is negligible,r and Z
c can be rewritten as follows:
(4)
r
=
1
a
s
,
Z
c
=
2
h
w
,
h
w
=
a
Q
0
(
2
g
A
H
0
)
.With the hydraulic friction losses being trivial andH
B
(
s
)
=
0 (tunnel connects with reservoir directly), the head and flow function is simplified as follows:
(5)
H
A
(
s
)
Q
A
(
s
)
=
-
Z
c
s
h
(
r
Δ
x
)
c
h
(
r
Δ
x
)
=
-
Z
c
t
h
(
r
Δ
x
)
.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 G
h
(
s
) at point A is obtained as
(6)
G
h
(
s
)
=
-
2
h
w
t
h
(
0.5
T
r
s
)
=
-
2
h
w
∑
i
=
0
n
(
(
0.5
T
r
s
)
2
i
+
1
/
(
2
i
+
1
)
!
)
∑
i
=
0
n
(
(
0.5
T
r
s
)
2
i
/
(
2
i
)
!
)
.Typically, there are the two models,i
=
0,1, as follows:
(7)
G
h
(
s
)
=
-
2
h
w
t
h
(
0.5
T
r
s
)
=
-
2
h
w
(
1
/
48
)
T
r
3
s
3
+
(
1
/
2
)
T
r
s
(
1
/
8
)
T
r
2
s
2
+
1
,
(8)
G
h
(
s
)
=
-
2
h
w
t
h
(
0.5
T
r
s
)
=
-
T
w
s
,
where T
r
=
2
Δ
x
/
a, T
w
=
L
Q
0
/
A
g
H
0. 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:
(9)
L
=
∑
i
=
1
n
L
i
,
A
=
∑
i
=
1
n
L
i
∑
i
=
1
n
(
L
i
/
A
i
)
,
T
r
=
∑
i
=
1
n
2
L
i
a
i
,
T
w
=
∑
i
=
1
n
(
L
i
/
A
i
)
Q
0
g
H
0
where L
i, A
i, a
i denote the length, cross-sectional area, and the velocity of pipeline i, respectively, and there are n 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:(10)
m
t
=
e
x
x
+
e
y
y
+
e
h
h
,
q
=
e
q
x
x
+
e
q
y
y
+
e
q
h
h
.The six constants of hydroturbinee
x, e
y, e
h, e
q
x, e
q
y, e
q
h 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
(11)
P
m
=
m
t
+
Δ
ω
.As the unit speed changes little, the speed deviationΔ
ω
=
0, leading to P
m
=
m
t. Then the transfer function of the turbine and conduit system is
(12)
G
t
(
s
)
=
e
y
1
+
e
G
h
(
s
)
1
-
e
q
h
G
h
(
s
)
,
where G
h
(
s
) is transfer function of conduit system defined in (6).In case of rigid water hammer, the coefficientn in (6) is 0 and the dynamic equation for the turbine output torque, as shown in Figure 1, becomes
(13)
m
˙
t
=
1
e
q
h
T
w
[
-
m
t
+
e
y
y
-
e
e
y
T
w
T
y
(
u
-
y
)
]
.Figure 1
Dynamic model of turbine and conduit system.


### 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
(14)
δ
˙
=
ω
0
ω
,
ω
˙
=
1
T
a
b
(
m
t
-
m
e
-
K
ω
)
,
where δ, ω, K, and T
a
b denote the rotor angle, relative speed deviation, damping coefficient, and mechanical starting time, respectively. The electromagnetic torque of the generator m
e is equal to its electromagnetic power P
e:
(15)
m
e
=
P
e
.The electromagnetic power can be calculated with the following formula:(16)
P
e
=
E
q
′
V
s
x
d
x
∑
′
sin
δ
+
V
s
2
2
x
d
x
∑
′
-
x
q
x
∑
x
d
x
∑
′
x
q
x
∑
sin
2
δ
,
where the effects of speed deviations, damping coefficient, and torque variations are all included in the analysis of the generator dynamic characteristics
(17)
x
d
∑
′
=
x
d
′
+
x
T
+
1
2
x
L
,
x
q
∑
=
x
q
+
x
T
+
1
2
x
L
.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(18)
G
s
(
s
)
=
1
1
+
T
y
s
,
where T
y 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
(19)
G
2
(
s
)
=
(
k
p
+
k
i
s
+
k
d
s
)
.Substituting (19) into (18), the results can be written in state-space form. At last, the mixed function can be obtained as follows:
(20)
d
y
d
t
=
1
T
y
(
-
k
p
ω
-
k
i
∫
Δ
ω
-
k
d
ω
˙
-
y
)
.Based on the discussions above, the differential equations that coupled each individual part of the turbine nonlinear control system can be written as(21)
δ
˙
=
ω
0
ω
,
ω
˙
=
1
T
a
b
(
m
t
-
D
ω
-
E
q
′
V
s
x
d
x
∑
′
sin
δ
m
m
m
m
m
m
-
V
s
2
2
x
d
x
∑
′
-
x
q
x
∑
x
d
x
∑
′
x
q
x
∑
sin
2
δ
)
,
m
˙
t
=
1
e
q
h
T
w
(
e
e
y
T
w
T
y

-
m
t
+
e
y
y
mmmmm
-
e
e
y
T
w
T
y
(
-
k
p
ω
-
k
i
ω
0
δ
-
k
d
ω
˙
-
y
)
)
,
y
˙
=
1
T
y
(
-
k
p
ω
-
k
i
ω
0
δ
-
k
d
ω
˙
-
y
)
.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(
0,0
,
0,0
), the following condition has to be satisfied:
(22)
ω
=
0
,
m
t
=
e
y
y
,
y
=
-
k
i
δ
ω
0
,
e
y
k
i
δ
ω
0
+
E
q
′
V
s
x
d
x
∑
′
sin
δ
+
V
s
2
2
x
d
x
∑
′
-
x
q
x
∑
x
d
x
∑
′
x
q
x
∑
sin
2
δ
=
0
.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.


## 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
(1)
Momentum equation:
∂
H
∂
x
+
1
g
A
∂
Q
∂
t
+
f
Q
2
2
g
D
A
2
=
0
,
Continuity equation:
∂
Q
∂
x
+
g
A
a
2
∂
H
∂
t
=
0
,
where D, f, A are parameters of the penstock. They denote the diameter, head loss, and area of the pipeline, respectively. H, Q are hydraulic head and turbine flow in penstock in operating condition, a is pressure wave velocity, and x 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]
(2)
[
H
A
(
s
)
Q
A
(
s
)
]
=
[
c
h
(
r
Δ
x
)
-
Z
c
s
h
(
r
Δ
x
)
-
s
h
(
r
Δ
x
)
Z
c
c
h
(
r
Δ
x
)
]
[
H
B
(
s
)
Q
B
(
s
)
]
.SubscriptsB, A are symbols of upstream and downstream section of pipeline, respectively. r and Z
c are the composite equations of parameters of the penstock, r
=
L
C
s
2
+
R
C
s, Z
c
=
r
/
C
s.L, R, C can be written as
(3)
L
=
Q
0
g
A
H
0
,
C
=
g
A
H
0
(
a
2
Q
0
)
,
R
=
(
f
Q
0
2
)
(
g
D
A
2
H
0
)
.Providing that the head loss is negligible,r and Z
c can be rewritten as follows:
(4)
r
=
1
a
s
,
Z
c
=
2
h
w
,
h
w
=
a
Q
0
(
2
g
A
H
0
)
.With the hydraulic friction losses being trivial andH
B
(
s
)
=
0 (tunnel connects with reservoir directly), the head and flow function is simplified as follows:
(5)
H
A
(
s
)
Q
A
(
s
)
=
-
Z
c
s
h
(
r
Δ
x
)
c
h
(
r
Δ
x
)
=
-
Z
c
t
h
(
r
Δ
x
)
.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 G
h
(
s
) at point A is obtained as
(6)
G
h
(
s
)
=
-
2
h
w
t
h
(
0.5
T
r
s
)
=
-
2
h
w
∑
i
=
0
n
(
(
0.5
T
r
s
)
2
i
+
1
/
(
2
i
+
1
)
!
)
∑
i
=
0
n
(
(
0.5
T
r
s
)
2
i
/
(
2
i
)
!
)
.Typically, there are the two models,i
=
0,1, as follows:
(7)
G
h
(
s
)
=
-
2
h
w
t
h
(
0.5
T
r
s
)
=
-
2
h
w
(
1
/
48
)
T
r
3
s
3
+
(
1
/
2
)
T
r
s
(
1
/
8
)
T
r
2
s
2
+
1
,
(8)
G
h
(
s
)
=
-
2
h
w
t
h
(
0.5
T
r
s
)
=
-
T
w
s
,
where T
r
=
2
Δ
x
/
a, T
w
=
L
Q
0
/
A
g
H
0. 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:
(9)
L
=
∑
i
=
1
n
L
i
,
A
=
∑
i
=
1
n
L
i
∑
i
=
1
n
(
L
i
/
A
i
)
,
T
r
=
∑
i
=
1
n
2
L
i
a
i
,
T
w
=
∑
i
=
1
n
(
L
i
/
A
i
)
Q
0
g
H
0
where L
i, A
i, a
i denote the length, cross-sectional area, and the velocity of pipeline i, respectively, and there are n 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:(10)
m
t
=
e
x
x
+
e
y
y
+
e
h
h
,
q
=
e
q
x
x
+
e
q
y
y
+
e
q
h
h
.The six constants of hydroturbinee
x, e
y, e
h, e
q
x, e
q
y, e
q
h 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
(11)
P
m
=
m
t
+
Δ
ω
.As the unit speed changes little, the speed deviationΔ
ω
=
0, leading to P
m
=
m
t. Then the transfer function of the turbine and conduit system is
(12)
G
t
(
s
)
=
e
y
1
+
e
G
h
(
s
)
1
-
e
q
h
G
h
(
s
)
,
where G
h
(
s
) is transfer function of conduit system defined in (6).In case of rigid water hammer, the coefficientn in (6) is 0 and the dynamic equation for the turbine output torque, as shown in Figure 1, becomes
(13)
m
˙
t
=
1
e
q
h
T
w
[
-
m
t
+
e
y
y
-
e
e
y
T
w
T
y
(
u
-
y
)
]
.Figure 1
Dynamic model of turbine and conduit system.


## 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
(14)
δ
˙
=
ω
0
ω
,
ω
˙
=
1
T
a
b
(
m
t
-
m
e
-
K
ω
)
,
where δ, ω, K, and T
a
b denote the rotor angle, relative speed deviation, damping coefficient, and mechanical starting time, respectively. The electromagnetic torque of the generator m
e is equal to its electromagnetic power P
e:
(15)
m
e
=
P
e
.The electromagnetic power can be calculated with the following formula:(16)
P
e
=
E
q
′
V
s
x
d
x
∑
′
sin
δ
+
V
s
2
2
x
d
x
∑
′
-
x
q
x
∑
x
d
x
∑
′
x
q
x
∑
sin
2
δ
,
where the effects of speed deviations, damping coefficient, and torque variations are all included in the analysis of the generator dynamic characteristics
(17)
x
d
∑
′
=
x
d
′
+
x
T
+
1
2
x
L
,
x
q
∑
=
x
q
+
x
T
+
1
2
x
L
.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(18)
G
s
(
s
)
=
1
1
+
T
y
s
,
where T
y 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
(19)
G
2
(
s
)
=
(
k
p
+
k
i
s
+
k
d
s
)
.Substituting (19) into (18), the results can be written in state-space form. At last, the mixed function can be obtained as follows:
(20)
d
y
d
t
=
1
T
y
(
-
k
p
ω
-
k
i
∫
Δ
ω
-
k
d
ω
˙
-
y
)
.Based on the discussions above, the differential equations that coupled each individual part of the turbine nonlinear control system can be written as(21)
δ
˙
=
ω
0
ω
,
ω
˙
=
1
T
a
b
(
m
t
-
D
ω
-
E
q
′
V
s
x
d
x
∑
′
sin
δ
m
m
m
m
m
m
-
V
s
2
2
x
d
x
∑
′
-
x
q
x
∑
x
d
x
∑
′
x
q
x
∑
sin
2
δ
)
,
m
˙
t
=
1
e
q
h
T
w
(
e
e
y
T
w
T
y

-
m
t
+
e
y
y
mmmmm
-
e
e
y
T
w
T
y
(
-
k
p
ω
-
k
i
ω
0
δ
-
k
d
ω
˙
-
y
)
)
,
y
˙
=
1
T
y
(
-
k
p
ω
-
k
i
ω
0
δ
-
k
d
ω
˙
-
y
)
.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(
0,0
,
0,0
), the following condition has to be satisfied:
(22)
ω
=
0
,
m
t
=
e
y
y
,
y
=
-
k
i
δ
ω
0
,
e
y
k
i
δ
ω
0
+
E
q
′
V
s
x
d
x
∑
′
sin
δ
+
V
s
2
2
x
d
x
∑
′
-
x
q
x
∑
x
d
x
∑
′
x
q
x
∑
sin
2
δ
=
0
.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 systemy
=
F
(
x
,
μ
), x
∈
R
4, μ
∈
R
1 is the bifurcation parameter, x
=
0 is the equilibrium point, and the Jacobi matrix characteristic polynomial at the equilibrium point is
(23)
f
(
λ
,
μ
)
=
λ
4
+
p
3
(
λ
)
λ
3
+
p
2
(
λ
)
λ
2
+
p
1
(
λ
)
λ
+
p
0
(
λ
)
.
Ifμ
=
0, the following conditions hold.
(1) The coefficientp
i
>
0 (i
=
0,1
,
2,3),
(24)
p
3
p
2
p
1
=
p
3
2
p
0
+
p
1
2
.
(2)(25)
p
0
′
≠
(
p
3
p
0
p
1
-
p
1
p
3
)
p
3
p
1
′
-
p
3
′
p
1
p
3
2
+
p
2
′
p
1
p
3
,
where *
′
=
d
*
/
d
μ, when |
μ
| is sufficiently small, Hopf bifurcation exists in one side of μ
=
0, and the period of limit cycle can be described as follows:
(26)
T
=
2
π
p
3
p
1
.
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ω
0
=
314, T
w
=
0.8 s, T
a
b
=
9.0 s, T
y
=
0.1 s, D
=
2, E
q
′
=
1.35, V
s
=
1, x
d
x
∑
′
=
1.15, x
q
x
∑
′
=
1.474, e
q
h
=
0.5, e
y
=
1.0, and e
=
0.7, 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 (
0,0
,
0,0
):
(27)
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
.It requires all the characteristic polynomial coefficientsp
i
>
0 (i
=
0,1
,
2,3) by condition 1 when using Theorem 1 to determine the existence of Hopf bifurcation. For example, p
3
>
0; then K
d
<
8.17857, which means that the Hopf bifurcation area is limited to the range of K
d
<
8.17857. Similarly, the ranges of K
i and K
p can be obtained. After substituting p
i 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.Figure 2
Bifurcation surface ofK
p, K
i, K
d.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. Figure2 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 μ
d
=
K
d
-
K
d
* or μ
p
=
K
p
-
K
p
* 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 Theorem1, when |
μ
| is sufficiently small, Hopf bifurcation phenomenon of the system will take place. For example, when K
i
*
=
1 and K
p
*
=
2, the bifurcation parameters K
d
1
*
=
0.313 and K
d
2
*
=
2.036 can be calculated from (24). K
d
1
* and K
d
2
* can also be got by the intersection values between the curve (K
i
=
1) and the line (K
p
=
2), as shown in Figure 3. The curve (K
i
=
0,1
,
3) is the boundary between stability region and limit cycle; point that located in limit cycle, for example, K
i
=
1, K
d
=
2.5, and K
p
=
2, would lead to the undesirable oscillation. Figure 3 is more intuitive to tell the stale and unstable PID parameters compared with Figure 2.Figure 3
Curves of bifurcation points withK
i
=
0,1
,
3.Figure4 is the bifurcation diagram of unit speed ω. When K
p
=
2, K
i
=
1, and K
d is the bifurcation parameter, there are two bifurcation points in the figure; that is, K
d
1
=
0.313 and K
d
2
=
2.036 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 K
d
1
<
K
d
<
K
d
2; it will converge to a stable limit cycle when K
d
<
K
d
1 or K
d
>
K
d
2, and the amplitude of the limit cycle oscillation corresponds to the ordinate value of limit cycle curve. The bifurcation thresholds K
d
1 and K
d
2 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 K
d. When the system is in stable domain (K
d
1
*
<
K
d
<
K
d
2
*), 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 (
K
d
1
*
,
K
d
2
*
), the eigenvalues of Jacobi matrix have complex conjugate eigenvalues with zero (
-
4.5
e
-
5
,
-
2.7
e
-
6
) real part; that is, the system is in a supercritical state. When K
d
<
K
d
1
* or K
d
>
K
d
2
*, 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
Changes of system properties with the bifurcation parameterK
d.
Bifurcation parameterK
d
K
d < K
d
1
*
K
d = K
d
1
* = 0.313
K
d
1
* < K
d < K
d
2
*
K
d = K
d
2
* = 2.036
K
d > K
d
2
*
Jacobi matrix eigenvalues
Two complex conjugate eigenvalues, negative real part one and positive real part one
-
2.46
-
9.77
-
4.5
e
-
5
±
5.98
i
Two complex conjugate eigenvalues with negative real parts
-
2.76
-
6.80
-
2.7
e
-
6
±
6.77
i
Two complex conjugate eigenvalues, negative real part one and positive real part one
Property
Stable limit cycle
center of bifurcation  in a supercritical state
Stable focus, stable domain
Center of bifurcation in a supercritical state
Stable limit cycleFigure 4
Bifurcation diagram of rotor speedω.Table1 provides theoretical e 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 (K
d
1
*
=
0.313 and K
d
2
*
=
2.036) and the nearer to middle place the K
d is, the quicker the convergence rate is and the more stable the system is. In practical applications, value of K
d 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.Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
0.309.
(a)
(b)Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, K
d
=
0.6.
(a)
(b)Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
0.8.
(a)
(b)Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
1.15.
(a)
(b)Time domain response waveforms ofδ and ω after stabilization can be seen in Figure 5 where the PID parameters are K
i
=
1, K
p
=
2, and K
d
=
0.309, 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 K
d
=
0.309, 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 p
3
=
12.24, p
1
=
438.15 easily from (27), so the period of limit cycle T
=
1.050 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 K
d
=
0.6,0.8,1.15, which are beneath the limit cycle curve (K
i
=
1), 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 K
d is farther from the two critical points (K
d
1
=
0.313 and K
d
2
=
2.036) and nearer to middle place (1.175), the convergence rate is quicker and the system is more stable.Figures9 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 μ
d
=
0.002; thus K
i
=
1, K
p
=
2, and K
d
=
2.038. The system shows a diffused and nonlinear growth oscillations curve. Hopf bifurcation occurs and the system tends to be a stable oscillation state where -
0.06
<
m
t
<
0.06, -
0.04
<
ω
<
0.04, -
1
<
δ
<
1, -
0.04
<
y
<
0.04, 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.Figure 9
Phase orbit ofδ, ω, and m
t.Figure 10
Phase orbit ofm
t, ω, and y.


## 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 systemy
=
F
(
x
,
μ
), x
∈
R
4, μ
∈
R
1 is the bifurcation parameter, x
=
0 is the equilibrium point, and the Jacobi matrix characteristic polynomial at the equilibrium point is
(23)
f
(
λ
,
μ
)
=
λ
4
+
p
3
(
λ
)
λ
3
+
p
2
(
λ
)
λ
2
+
p
1
(
λ
)
λ
+
p
0
(
λ
)
.
Ifμ
=
0, the following conditions hold.
(1) The coefficientp
i
>
0 (i
=
0,1
,
2,3),
(24)
p
3
p
2
p
1
=
p
3
2
p
0
+
p
1
2
.
(2)(25)
p
0
′
≠
(
p
3
p
0
p
1
-
p
1
p
3
)
p
3
p
1
′
-
p
3
′
p
1
p
3
2
+
p
2
′
p
1
p
3
,
where *
′
=
d
*
/
d
μ, when |
μ
| is sufficiently small, Hopf bifurcation exists in one side of μ
=
0, and the period of limit cycle can be described as follows:
(26)
T
=
2
π
p
3
p
1
.
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ω
0
=
314, T
w
=
0.8 s, T
a
b
=
9.0 s, T
y
=
0.1 s, D
=
2, E
q
′
=
1.35, V
s
=
1, x
d
x
∑
′
=
1.15, x
q
x
∑
′
=
1.474, e
q
h
=
0.5, e
y
=
1.0, and e
=
0.7, 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 (
0,0
,
0,0
):
(27)
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
.It requires all the characteristic polynomial coefficientsp
i
>
0 (i
=
0,1
,
2,3) by condition 1 when using Theorem 1 to determine the existence of Hopf bifurcation. For example, p
3
>
0; then K
d
<
8.17857, which means that the Hopf bifurcation area is limited to the range of K
d
<
8.17857. Similarly, the ranges of K
i and K
p can be obtained. After substituting p
i 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.Figure 2
Bifurcation surface ofK
p, K
i, K
d.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. Figure2 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 μ
d
=
K
d
-
K
d
* or μ
p
=
K
p
-
K
p
* 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 Theorem1, when |
μ
| is sufficiently small, Hopf bifurcation phenomenon of the system will take place. For example, when K
i
*
=
1 and K
p
*
=
2, the bifurcation parameters K
d
1
*
=
0.313 and K
d
2
*
=
2.036 can be calculated from (24). K
d
1
* and K
d
2
* can also be got by the intersection values between the curve (K
i
=
1) and the line (K
p
=
2), as shown in Figure 3. The curve (K
i
=
0,1
,
3) is the boundary between stability region and limit cycle; point that located in limit cycle, for example, K
i
=
1, K
d
=
2.5, and K
p
=
2, would lead to the undesirable oscillation. Figure 3 is more intuitive to tell the stale and unstable PID parameters compared with Figure 2.Figure 3
Curves of bifurcation points withK
i
=
0,1
,
3.Figure4 is the bifurcation diagram of unit speed ω. When K
p
=
2, K
i
=
1, and K
d is the bifurcation parameter, there are two bifurcation points in the figure; that is, K
d
1
=
0.313 and K
d
2
=
2.036 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 K
d
1
<
K
d
<
K
d
2; it will converge to a stable limit cycle when K
d
<
K
d
1 or K
d
>
K
d
2, and the amplitude of the limit cycle oscillation corresponds to the ordinate value of limit cycle curve. The bifurcation thresholds K
d
1 and K
d
2 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 K
d. When the system is in stable domain (K
d
1
*
<
K
d
<
K
d
2
*), 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 (
K
d
1
*
,
K
d
2
*
), the eigenvalues of Jacobi matrix have complex conjugate eigenvalues with zero (
-
4.5
e
-
5
,
-
2.7
e
-
6
) real part; that is, the system is in a supercritical state. When K
d
<
K
d
1
* or K
d
>
K
d
2
*, 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
Changes of system properties with the bifurcation parameterK
d.
Bifurcation parameterK
d
K
d < K
d
1
*
K
d = K
d
1
* = 0.313
K
d
1
* < K
d < K
d
2
*
K
d = K
d
2
* = 2.036
K
d > K
d
2
*
Jacobi matrix eigenvalues
Two complex conjugate eigenvalues, negative real part one and positive real part one
-
2.46
-
9.77
-
4.5
e
-
5
±
5.98
i
Two complex conjugate eigenvalues with negative real parts
-
2.76
-
6.80
-
2.7
e
-
6
±
6.77
i
Two complex conjugate eigenvalues, negative real part one and positive real part one
Property
Stable limit cycle
center of bifurcation  in a supercritical state
Stable focus, stable domain
Center of bifurcation in a supercritical state
Stable limit cycleFigure 4
Bifurcation diagram of rotor speedω.Table1 provides theoretical e 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 (K
d
1
*
=
0.313 and K
d
2
*
=
2.036) and the nearer to middle place the K
d is, the quicker the convergence rate is and the more stable the system is. In practical applications, value of K
d 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.Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
0.309.
(a)
(b)Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, K
d
=
0.6.
(a)
(b)Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
0.8.
(a)
(b)Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
1.15.
(a)
(b)Time domain response waveforms ofδ and ω after stabilization can be seen in Figure 5 where the PID parameters are K
i
=
1, K
p
=
2, and K
d
=
0.309, 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 K
d
=
0.309, 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 p
3
=
12.24, p
1
=
438.15 easily from (27), so the period of limit cycle T
=
1.050 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 K
d
=
0.6,0.8,1.15, which are beneath the limit cycle curve (K
i
=
1), 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 K
d is farther from the two critical points (K
d
1
=
0.313 and K
d
2
=
2.036) and nearer to middle place (1.175), the convergence rate is quicker and the system is more stable.Figures9 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 μ
d
=
0.002; thus K
i
=
1, K
p
=
2, and K
d
=
2.038. The system shows a diffused and nonlinear growth oscillations curve. Hopf bifurcation occurs and the system tends to be a stable oscillation state where -
0.06
<
m
t
<
0.06, -
0.04
<
ω
<
0.04, -
1
<
δ
<
1, -
0.04
<
y
<
0.04, 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.Figure 9
Phase orbit ofδ, ω, and m
t.Figure 10
Phase orbit ofm
t, ω, and y.


## 4. Analysis and Simulation for Chaos

WhenK
d
>
2.5, with the continually increasing of K
d, unit speed ω on bifurcation diagram remains as a single rising curve until K
d
=
5.72. The bifurcation diagram gets an intricate pattern when K
d
>
5.72, 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.Figure 11
Bifurcation diagram ofω with K
p
=
2, K
i
=
1.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 differentK
d, as presented in Figure 12. When K
d
≦
5.72, the largest Lyapunov exponent LEl is zero, indicating that the system is in periodic motion state (limit cycle). In the vicinity of 5.72
≦
K
d
≦
6.0, 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 K
d
=
8.0 and means chaotic motion that leads to the unstable state of the control system exists.Figure 12
Lyapunov exponent withK
d.Figure13 shows the system stability domain, limit cycles, and chaotic region of the system when K
i
=
1.0, 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 K
d when the maximum Lyapunov exponent LEl is greater than zero.Figure 13
Stability domain, limit cycles, and chaotic region ofK
p, K
d with K
i
=
1.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. WhenK
p
=
2.0, K
i
=
1.0, and K
d
=
6.5, chaotic attractors that are similar to the one with scroll structure in Chua’s circuit [37] generate in the δ
,
ω
,
m
t and δ
,
ω
,
y 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.Figure 14
The chaotic attractor withK
i
=
1, K
p
=
2, and K
d
=
6.5.Figure 15
The chaotic attractor withK
i
=
1, K
p
=
2, and K
d
=
6.5.Time waveforms ofδ and ω with K
i
=
1, K
p
=
2, and K
d
=
6.5.
(a)
(b)Poincare maps are shown in Figures17 and 18. It can be seen that numerous intersections of the continuous trajectory turn up in the Poincare maps on section of ω
=
0 and m
t
=
0. However, they neither fill with the entire phase space nor distribute discretely in the δ
,
y
,
m
t and δ
,
ω
,
y 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.Figure 17
Poincare map withω
=
0.Figure 18
Poincare map withm
t
=
0.


## 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.



---

*Source: 412578-2014-06-25.xml*