# Singular Perturbation Theory-Based Qualitative Dynamics Investigation of Flywheel Energy Storage System in Discharge Mode

Authors: Weiya Zhang; Yongli Li; Xiaoyong Chang; Nan Wang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

An investigation on qualitative dynamics in a voltage-current dual-loop controlled flywheel energy storage system (FESS) operating in discharge mode is presented in this paper, providing novel insights into the effect of two-timescale characteristics on the safety and stability of energy transmission of FESS. Based on singular perturbation theory, a two-timescale approach is proposed to separate the FESS into the fast and slow subsystems. Stability analysis of the transient fixed points confirms the effects of systemic parameters on FESS’s dynamics and indicates that the FESS shifts from the spiking state to the quiescent state when the slow variable crosses the bifurcation point of the fast subsystem. Mechanism analysis reveals that the root cause of the qualitative dynamics is the voltage instability of the FESS. Moreover, the feasibility boundaries of key parameters are derived, and application requirements of the proposed approach are also discussed, guiding the extension of the approach to engineering applications and solving the dynamics analysis problem to some extent at a theoretical analysis level. Constant voltage discharge experiment is performed based on the FESS test bench built in Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, which validates the theoretical results.


---

## Body

## 1. Introduction

Permanent magnet (PM) brushless dc motor (BLDCM) controlled FESS, with the advantages of high density, low maintenance, long lifetime, and good compactness, has become a new trend for energy storage, applied more and more in the uninterruptible power supply, rail transportation, and smart grids [1, 2]. As shown in Figure 1, the FESS mainly consists of two self-contained parts, that is, the mechanical part (flywheel motor system) and the electrical part (power drive system). In general, the flywheel motor is controlled by power electronic circuits. Due to the existence of the intrinsic nonlinearity, various nonlinear dynamics occur during the operation of FESS when the system state changes, which has an influence on the safety and stability of energy transmission.FESS test bench. (a) Flywheel motor system, (b) power drive system.
(a)
(b)So far, rotor dynamic problems of FESS’s mechanical part have been seen as a cause of decreased rotor dynamic performance and reduced stability [3–9]. However, few researchers have revealed that, as a strongly coupled system, FESS has more complex nonlinear dynamics due to the interaction and difference between the mechanical part and the electrical part. This complex nonlinear dynamics has much more direct influence on the safety and stability of energy transmission and thus affects the safety and stability of FESS. Zhang et al. have investigated the nonlinear dynamics of FESS from the viewpoint of the interaction between the mechanical part and the electrical part [10]. Turning to the difference of the two parts we find that the electrical variables have significantly faster dynamics than the mechanical variables. As such, the dynamics of the voltage and current are faster than that of the rotate speed of flywheel, making the FESS be a typical two-timescale system [11]. It is well known that bursting phenomenon is observed when a slow variable controls the fast dynamics in some two-timescale systems such as neuronal systems and biological systems [12, 13]. Bursting is a state of switching between the spiking state (SP) and the quiescent state (QS). Generally, QS indicates all the variables are at rest or exhibit small amplitude oscillations, while SP indicates variables may behave in large amplitude oscillations.In this paper, we aim to present evidences that as FESS operating in discharge mode, a small change in parameter values around the bifurcation points of FESS’s fast system will lead to qualitative dynamics of the full-system, and investigate the effect of two-timescale characteristics on such dynamics, which is similar to nonrecurrent bursting. The analysis of dynamical systems with two timescales is a subject whose history interweaves three different viewpoints: nonstandard analysis [14], classical asymptotics [15], and geometric singular perturbation theory [16]. The first two methods lead to relatively large errors, and the geometric singular perturbation method is used to get the analytical solution of simple multiple timescale systems. The two-timescale approach proposed by Rinzel [17] is the classic approach to deeply investigate the two-timescale bifurcation dynamics, which gives a full description of the steady state, and periodic solution set of the fast subsystem, reflecting the global bifurcation structure of the fast subsystem with the slow variables treated as parameters [18]. First, stability analysis of the transient fixed points is proposed to study the bifurcation set of the fast subsystem, showing that as the slow variable varies, the fast subsystem loses stability from the originally stable state to Hopf bifurcation, and the dynamical evolution of the full-system is close in accordance with that of the fast subsystem. Not only reveling that the operation state of the FESS shifts when the slow variable crosses the bifurcation point of the fast subsystem, but also giving a way to predict the occurrence and evolution of qualitative dynamics of FESS in discharge mode. Then, the bifurcation mechanism analysis of the fast subsystem is proposed, offering an intuitive explanation of the origin of the nonrecurrent dynamics of the full-system. Furthermore, the feasibility regions are shown and provide instructions to parameters setting of FESS. Finally, the application requirements of the proposed approach are also discussed, guiding the extension of the approach to dynamics analysis of other electromechanical coupling systems.This paper is organized as follows. In Section2, the normalized dynamic model of FESS is established and numerical simulations have been taken. In Section 3, the two-timescale approach based on singular perturbation theory is proposed and applied. A brief analysis of the application of the proposed approach is shown in Section 4. Also, the observed instability phenomena are observed experimentally, as presented in Section 5. Finally, Section 6 concludes this paper.


## 2. Modelling and Two-Timescale Characteristics of FESS

As shown in Figure2, the modelling of FESS includes the flywheel motor system (the flywheel rotor driven by BLDCM) and the power drive system (the electrical subsystem and feedback control subsystem). The flywheel motor system is designed with a full bridge (IGBTs) at its output electrical terminals and DC-DC converter at the dc link. While diodes perform uncontrolled rectification, the DC-DC converter adjusts the voltage of the BLDCM in order to make it suitable for the load. Electronic commutation is achieved using a microprocessor-based controller with a Hall-effect position and a current sensor as input to generate gating signals for IGBTs.Figure 2
Structure of the FESS.


### 2.1. Modelling of FESS

Before modeling the FESS, five assumptions of the flywheel motor system are described as follows: (a) the saturation of the core is neglected; (b) the losses of eddy and hysteresis are ignored; (c) the distribution of air gap is uniform; (d) the self-inductance and mutual inductance among the windings are independent of the position of the rotor; (e) ignore the commutation process. The physical structure of the flywheel motor system is shown in Figure3(a), while the schematic diagram is shown in Figure 3(b). It has been assumed that the phase resistance r
m, the self-inductance L, and mutual inductance M of all the windings are equal. Assuming further that there is no change in the rotor reluctance with angle, hence, the circuit equations of the three windings in phase variables are
(1)
[
e
a
e
b
e
c
]
=
[
r
m
0
0
0
r
m
0
0
0
r
m
]
[
i
A
i
B
i
C
]
+
[
L
M
M
M
L
M
M
M
L
]
d
d
t
[
i
A
i
B
i
C
]
+
[
u
A
u
B
u
C
]
,
where u
A, u
B, and u
C are phase voltages, i
A, i
B, and i
C are phase currents, and e
a, e
b, and e
c are the induced back EMFs. For i
A
+
i
B
+
i
C
=
0 and denoting L
-
M as l
m, then the coupled circuit of the stator windings in terms of the machine electrical constants can be derived in Figure 3(c). The circuit equation can be written as
(2)
[
e
a
e
b
e
c
]
=
[
r
m
0
0
0
r
m
0
0
0
r
m
]
[
i
A
i
B
i
C
]
+
[
l
m
0
0
0
l
m
0
0
0
l
m
]
d
d
t
[
i
A
i
B
i
C
]
+
[
u
A
u
B
u
C
]
.Flywheel motor system: (a) physical structure, (b) schematic diagram, (c) three-phase equivalent circuit, (d) two-phase conducting equivalent circuit (the current direction corresponds to the generator operation).
(a)
(b)
(c)
(d)According to the switching pattern (commutation function), only two phases are active at the same time, while the third one is silent (Figure3(d)) [19]. For example, the voltage equation during the operation of phase A
-
B can be written as
(3)
u
2
=
-
2
r
m
i
a
-
2
l
m
d
i
a
d
t
+
(
e
A
-
e
B
)
=
-
R
m
i
a
-
L
m
d
i
a
d
t
+
K
e
ω
,
where K
e, B
v, J, R
m
=
2
r
m, and L
m
=
2
l
m denote line back-EMF constant, friction coefficient, moment of inertia, line resistance, and line inductance, respectively. Considering the equation of motion for BLDCM and treating the susceptibility and flux as constant, the differential equations for flywheel motor system are derived as
(4)
d
i
a
d
t
=
K
e
ω
L
m
-
R
m
i
a
L
m
-
u
2
L
m
,
d
ω
d
t
=
-
K
m
i
a
J
-
B
v
ω
J
,
where K
m denotes torque constant and the angular velocity ω, the phase current i
a, and the motor voltage u
2 are state variables.As shown in Figure4, the electrical subsystem include a bidirectional DC-DC converter, which consists of a transfer inductor L, capacitors C
1 and C
2, IGBT V
1 with diode V
D
1, and IGBT V
2 with diode V
D
2 and a three-phase full bridge converter which consists of IGBT-diode V
T
1
-
6. Mechanical energy is converted into electrical energy through the diodes of the three-phase full bridge converter when electrical energy is transmitted to the load via bidirectional DC-DC converter and regulated as load demands. During the discharge mode, the switching of V and V
D determines the two different modes, which are Mode 1 (n
T
<
t
<
n
T
+
d
T) when V
2 is on and V
D
1 is off and Mode 2 (n
T
+
d
T
<
t
<
(
n
+
1
)
T) when V
2 is off and V
D
1 is on. A dual-loop proportional-integral scheme is applied to regulate the duty cycle of switches, of which the PI coefficients are (K
p
v, K
i
v) and (K
p
i, K
i
i). The sampling coefficients of voltage and current are h
v
=
0.02424 and h
i
=
0.08.Figure 4
Schematic diagram of power drive system.Dimensionless parameters are introduced to put the system into standard form. Define the nominal output voltage asV
ref, the nominal angular velocity as ω
ref, and the load resistance as R. To be clear, the nomenclature for all variables, abbreviations, and parameters are listed in Tables 1, 2, and 3. With these definitions, the other variables can be normalized on the basis of t
=
L
C
1
τ, X
1
(
i
L
)
=
i
L
R
/
V
ref, X
2
(
u
1
)
=
u
1
/
V
ref, X
3
(
u
2
)
=
u
2
/
V
ref, X
4
(
i
a
)
=
i
a
R
/
V
ref, X
5
(
ω
)
=
ω
/
ω
ref, X
6
(
x
1
c
)
=
x
1
c, X
7
(
x
2
c
)
=
x
2
c, a
1
=
R
L
C
1
/
L, a
2
=
-
L
C
1
/
(
R
C
1
), a
3
=
L
C
1
/
(
R
C
2
), a
4
=
R
K
e
ω
ref
L
C
1
/
(
V
ref
L
m
), a
5
=
-
R
m
L
C
1
/
L
m, a
6
=
-
R
L
C
1
/
L
m, a
7
=
-
K
m
V
ref
L
C
1
/
(
J
R
ω
ref
), a
8
=
-
B
v
L
C
1
/
J, a
9
=
h
v
V
ref
L
C
1, a
10
=
-
h
i
V
ref
L
C
1
/
R, a
11
=
-
K
p
v
h
v
h
i
V
ref
L
C
1, a
12
=
K
i
v
h
i
L
C
1, b
1
=
a
9, and b
2
=
-
a
11.Table 1
Nomenclature for state variables.
State variables (with dimensions)
State variables (dimensionless)
Inductor currenti
L
X
1
(
i
L
)
Load voltageu
1
X
2
(
u
1
)
Motor voltageu
2
X
3
(
u
2
)
Phase currenti
a
X
4
(
i
a
)
Angular velocityω
X
5
(
ω
)
Outer corrector outputx
1
c
X
6
(
x
1
c
)
Inner corrector outputx
2
c
X
7
(
x
2
c
)Table 2
Nomenclature for abbreviations.
Systematic matrix coefficient
Input matrix coefficient
a
1, a
2, a
3, a
4, a
5, a
6, a
7, a
8, a
9, a
10, a
11
b
1, b
2Table 3
Nomenclature for parameters with values.
Parameter
Value
Voltage-loopK
p
v
60
Voltage-loopK
i
v
49
Current-loopK
p
i
40
Current-loopK
i
i
36
Nominal angular velocityω
ref
800 rads−1
Transfer inductanceL
0.7 mH
Filter capacitanceC
1
1 mF
Filter capacitanceC
2
1 mF
Line resistanceR
m
0.5 Ω
Line inductanceL
m
0.36 mH
Torque constantK
m
0.0673 Nm/A
Line back-EMF constantK
e
0.0625 V/rads−1
Friction coefficientB
v
0.0001 Nm/rads−1
Moment of inertiaJ
0.123 kgm2
Load resistanceR
5 ΩThen, the dimensionless switched dynamical equations of the present FESS are(5)
X
•
(
τ
)
=
A
1
X
(
τ
)
+
B
1
,
Mode
1
,
X
•
(
τ
)
=
A
2
X
(
τ
)
+
B
2
,
Mode
2
,
where
(6)
A
1
=
[
0
0
a
1
0
0
0
0
0
a
2
0
0
0
0
0
-
a
3
0
0
a
3
0
0
0
0
0
a
6
a
5
a
4
0
0
0
0
0
a
7
a
8
0
0
0
-
a
9
0
0
0
0
0
a
10
a
11
0
0
0
a
12
0
]
,
A
2
=
[
0
-
a
1
0
0
0
0
0
-
a
2
a
2
0
0
0
0
0
0
0
0
a
3
0
0
0
0
0
a
6
a
5
a
4
0
0
0
0
0
a
7
a
8
0
0
0
-
a
9
0
0
0
0
0
a
10
a
11
0
0
0
a
12
0
]
,
B
1
=
Β
2
=
[
0
0
0
0
0
b
1
b
2
]
,
X
(
τ
)
=
[
X
1
(
i
L
)
X
2
(
u
1
)
X
3
(
u
2
)
X
4
(
i
a
)
X
5
(
ω
)
X
6
(
x
1
c
)
X
7
(
x
2
c
)
]
T
.


### 2.2. Dynamic Characteristics

The operation principle of FESS is concerned with two processes: charge and discharge. In the charge mode, the flywheel motor system is driven by power grid and the electric energy is stored in the form of mechanical energy. Once the unit receives a signal of discharge, the flywheel rotor starts to decelerate and drives the BLDCM to generate electricity. This paper reports that the parametric regions have a significant effect on the two-timescale characteristics of FESS and the discharge performance degrades due to nonlinear dynamics when FESS’s systemic parameters fall into a certain area.Here, based on the exact state (5), a series of numerical simulations are carried out to make an initial evaluation of the possible dynamics. As shown in Table 3, the physical system parameters are set according to the actual FESS, and the PI control parameters of the feedback control subsystem are carefully designed to maintain closed-loop performance of the power drive system in spite of varying conditions. For energy storage unit as FESS, the stability in conjunction with the discharge targets is the primary consideration, so we restrict our attention to the case when V
ref varies, while the other parameters are fixed. Essentially, for each set of parameters, cycle-by-cycle time domain waveforms are generated by solving the appropriate linear equation in a subinterval of time, according to different switching states. The waveforms for different V
ref are shown in Figures 5–7.Time histories atV
ref
=
5
V, for (a) X
1
(
i
L
), (b) oscillatory state of (a), (c) steady state of (a), (d) X
2
(
u
1
) & X
3
(
u
2
).
(a)
(b)
(c)
(d)Time histories atV
ref
=
10
V, for (a) X
1
(
i
L
), (b) oscillatory state of (a), (c) steady state of (a), (d) X
2
(
u
1
) & X
3
(
u
2
).
(a)
(b)
(c)
(d)Time histories atV
ref
=
15
V, for (a) X
1
(
i
L
), (b) oscillatory state of (a), (c) steady state of (a), (d) X
2
(
u
1
) & X
3
(
u
2
).
(a)
(b)
(c)
(d)Figure5(a) shows that at V
ref
=
5
V, X
1
(
i
L
) takes a combined oscillation of large and small amplitudes as the flywheel rotor slows down. At the beginning of discharge, X
1
(
i
L
) shows a quasiperiodic oscillation with a high-value magnitude. Such oscillations decay rapidly when X
5
(
ω
) drops to 0.5; then X
1
(
i
L
) becomes steady. The close-up views are shown in Figures 5(b)-5(c). Dynamics of X
3
(
u
2
) and X
2
(
u
1
) as X
5
(
ω
) decreases are shown in Figure 5(d). For V
ref
=
10
V and V
ref
=
15
V shown in Figures 6-7, similar oscillations with different amplitudes and frequencies happen to X
1
(
i
L
). Although the control parameters are carefully designed according to the closed-loop transfer function of the power drive system, the interaction and difference between the mechanical part and the electrical part still cause complex nonlinear dynamics. It can be seen that larger V
ref determines stronger oscillations and X
5
(
ω
) changes much slower than other state variables. This indicates the existence of two-timescale characteristics in the discharge mode of FESS. Moreover, the oscillations of the electrical variables weaken along with the decrease of the mechanical variables, which will be proved related to the qualitative dynamics in the following section.


## 2.1. Modelling of FESS

Before modeling the FESS, five assumptions of the flywheel motor system are described as follows: (a) the saturation of the core is neglected; (b) the losses of eddy and hysteresis are ignored; (c) the distribution of air gap is uniform; (d) the self-inductance and mutual inductance among the windings are independent of the position of the rotor; (e) ignore the commutation process. The physical structure of the flywheel motor system is shown in Figure3(a), while the schematic diagram is shown in Figure 3(b). It has been assumed that the phase resistance r
m, the self-inductance L, and mutual inductance M of all the windings are equal. Assuming further that there is no change in the rotor reluctance with angle, hence, the circuit equations of the three windings in phase variables are
(1)
[
e
a
e
b
e
c
]
=
[
r
m
0
0
0
r
m
0
0
0
r
m
]
[
i
A
i
B
i
C
]
+
[
L
M
M
M
L
M
M
M
L
]
d
d
t
[
i
A
i
B
i
C
]
+
[
u
A
u
B
u
C
]
,
where u
A, u
B, and u
C are phase voltages, i
A, i
B, and i
C are phase currents, and e
a, e
b, and e
c are the induced back EMFs. For i
A
+
i
B
+
i
C
=
0 and denoting L
-
M as l
m, then the coupled circuit of the stator windings in terms of the machine electrical constants can be derived in Figure 3(c). The circuit equation can be written as
(2)
[
e
a
e
b
e
c
]
=
[
r
m
0
0
0
r
m
0
0
0
r
m
]
[
i
A
i
B
i
C
]
+
[
l
m
0
0
0
l
m
0
0
0
l
m
]
d
d
t
[
i
A
i
B
i
C
]
+
[
u
A
u
B
u
C
]
.Flywheel motor system: (a) physical structure, (b) schematic diagram, (c) three-phase equivalent circuit, (d) two-phase conducting equivalent circuit (the current direction corresponds to the generator operation).
(a)
(b)
(c)
(d)According to the switching pattern (commutation function), only two phases are active at the same time, while the third one is silent (Figure3(d)) [19]. For example, the voltage equation during the operation of phase A
-
B can be written as
(3)
u
2
=
-
2
r
m
i
a
-
2
l
m
d
i
a
d
t
+
(
e
A
-
e
B
)
=
-
R
m
i
a
-
L
m
d
i
a
d
t
+
K
e
ω
,
where K
e, B
v, J, R
m
=
2
r
m, and L
m
=
2
l
m denote line back-EMF constant, friction coefficient, moment of inertia, line resistance, and line inductance, respectively. Considering the equation of motion for BLDCM and treating the susceptibility and flux as constant, the differential equations for flywheel motor system are derived as
(4)
d
i
a
d
t
=
K
e
ω
L
m
-
R
m
i
a
L
m
-
u
2
L
m
,
d
ω
d
t
=
-
K
m
i
a
J
-
B
v
ω
J
,
where K
m denotes torque constant and the angular velocity ω, the phase current i
a, and the motor voltage u
2 are state variables.As shown in Figure4, the electrical subsystem include a bidirectional DC-DC converter, which consists of a transfer inductor L, capacitors C
1 and C
2, IGBT V
1 with diode V
D
1, and IGBT V
2 with diode V
D
2 and a three-phase full bridge converter which consists of IGBT-diode V
T
1
-
6. Mechanical energy is converted into electrical energy through the diodes of the three-phase full bridge converter when electrical energy is transmitted to the load via bidirectional DC-DC converter and regulated as load demands. During the discharge mode, the switching of V and V
D determines the two different modes, which are Mode 1 (n
T
<
t
<
n
T
+
d
T) when V
2 is on and V
D
1 is off and Mode 2 (n
T
+
d
T
<
t
<
(
n
+
1
)
T) when V
2 is off and V
D
1 is on. A dual-loop proportional-integral scheme is applied to regulate the duty cycle of switches, of which the PI coefficients are (K
p
v, K
i
v) and (K
p
i, K
i
i). The sampling coefficients of voltage and current are h
v
=
0.02424 and h
i
=
0.08.Figure 4
Schematic diagram of power drive system.Dimensionless parameters are introduced to put the system into standard form. Define the nominal output voltage asV
ref, the nominal angular velocity as ω
ref, and the load resistance as R. To be clear, the nomenclature for all variables, abbreviations, and parameters are listed in Tables 1, 2, and 3. With these definitions, the other variables can be normalized on the basis of t
=
L
C
1
τ, X
1
(
i
L
)
=
i
L
R
/
V
ref, X
2
(
u
1
)
=
u
1
/
V
ref, X
3
(
u
2
)
=
u
2
/
V
ref, X
4
(
i
a
)
=
i
a
R
/
V
ref, X
5
(
ω
)
=
ω
/
ω
ref, X
6
(
x
1
c
)
=
x
1
c, X
7
(
x
2
c
)
=
x
2
c, a
1
=
R
L
C
1
/
L, a
2
=
-
L
C
1
/
(
R
C
1
), a
3
=
L
C
1
/
(
R
C
2
), a
4
=
R
K
e
ω
ref
L
C
1
/
(
V
ref
L
m
), a
5
=
-
R
m
L
C
1
/
L
m, a
6
=
-
R
L
C
1
/
L
m, a
7
=
-
K
m
V
ref
L
C
1
/
(
J
R
ω
ref
), a
8
=
-
B
v
L
C
1
/
J, a
9
=
h
v
V
ref
L
C
1, a
10
=
-
h
i
V
ref
L
C
1
/
R, a
11
=
-
K
p
v
h
v
h
i
V
ref
L
C
1, a
12
=
K
i
v
h
i
L
C
1, b
1
=
a
9, and b
2
=
-
a
11.Table 1
Nomenclature for state variables.
State variables (with dimensions)
State variables (dimensionless)
Inductor currenti
L
X
1
(
i
L
)
Load voltageu
1
X
2
(
u
1
)
Motor voltageu
2
X
3
(
u
2
)
Phase currenti
a
X
4
(
i
a
)
Angular velocityω
X
5
(
ω
)
Outer corrector outputx
1
c
X
6
(
x
1
c
)
Inner corrector outputx
2
c
X
7
(
x
2
c
)Table 2
Nomenclature for abbreviations.
Systematic matrix coefficient
Input matrix coefficient
a
1, a
2, a
3, a
4, a
5, a
6, a
7, a
8, a
9, a
10, a
11
b
1, b
2Table 3
Nomenclature for parameters with values.
Parameter
Value
Voltage-loopK
p
v
60
Voltage-loopK
i
v
49
Current-loopK
p
i
40
Current-loopK
i
i
36
Nominal angular velocityω
ref
800 rads−1
Transfer inductanceL
0.7 mH
Filter capacitanceC
1
1 mF
Filter capacitanceC
2
1 mF
Line resistanceR
m
0.5 Ω
Line inductanceL
m
0.36 mH
Torque constantK
m
0.0673 Nm/A
Line back-EMF constantK
e
0.0625 V/rads−1
Friction coefficientB
v
0.0001 Nm/rads−1
Moment of inertiaJ
0.123 kgm2
Load resistanceR
5 ΩThen, the dimensionless switched dynamical equations of the present FESS are(5)
X
•
(
τ
)
=
A
1
X
(
τ
)
+
B
1
,
Mode
1
,
X
•
(
τ
)
=
A
2
X
(
τ
)
+
B
2
,
Mode
2
,
where
(6)
A
1
=
[
0
0
a
1
0
0
0
0
0
a
2
0
0
0
0
0
-
a
3
0
0
a
3
0
0
0
0
0
a
6
a
5
a
4
0
0
0
0
0
a
7
a
8
0
0
0
-
a
9
0
0
0
0
0
a
10
a
11
0
0
0
a
12
0
]
,
A
2
=
[
0
-
a
1
0
0
0
0
0
-
a
2
a
2
0
0
0
0
0
0
0
0
a
3
0
0
0
0
0
a
6
a
5
a
4
0
0
0
0
0
a
7
a
8
0
0
0
-
a
9
0
0
0
0
0
a
10
a
11
0
0
0
a
12
0
]
,
B
1
=
Β
2
=
[
0
0
0
0
0
b
1
b
2
]
,
X
(
τ
)
=
[
X
1
(
i
L
)
X
2
(
u
1
)
X
3
(
u
2
)
X
4
(
i
a
)
X
5
(
ω
)
X
6
(
x
1
c
)
X
7
(
x
2
c
)
]
T
.


## 2.2. Dynamic Characteristics

The operation principle of FESS is concerned with two processes: charge and discharge. In the charge mode, the flywheel motor system is driven by power grid and the electric energy is stored in the form of mechanical energy. Once the unit receives a signal of discharge, the flywheel rotor starts to decelerate and drives the BLDCM to generate electricity. This paper reports that the parametric regions have a significant effect on the two-timescale characteristics of FESS and the discharge performance degrades due to nonlinear dynamics when FESS’s systemic parameters fall into a certain area.Here, based on the exact state (5), a series of numerical simulations are carried out to make an initial evaluation of the possible dynamics. As shown in Table 3, the physical system parameters are set according to the actual FESS, and the PI control parameters of the feedback control subsystem are carefully designed to maintain closed-loop performance of the power drive system in spite of varying conditions. For energy storage unit as FESS, the stability in conjunction with the discharge targets is the primary consideration, so we restrict our attention to the case when V
ref varies, while the other parameters are fixed. Essentially, for each set of parameters, cycle-by-cycle time domain waveforms are generated by solving the appropriate linear equation in a subinterval of time, according to different switching states. The waveforms for different V
ref are shown in Figures 5–7.Time histories atV
ref
=
5
V, for (a) X
1
(
i
L
), (b) oscillatory state of (a), (c) steady state of (a), (d) X
2
(
u
1
) & X
3
(
u
2
).
(a)
(b)
(c)
(d)Time histories atV
ref
=
10
V, for (a) X
1
(
i
L
), (b) oscillatory state of (a), (c) steady state of (a), (d) X
2
(
u
1
) & X
3
(
u
2
).
(a)
(b)
(c)
(d)Time histories atV
ref
=
15
V, for (a) X
1
(
i
L
), (b) oscillatory state of (a), (c) steady state of (a), (d) X
2
(
u
1
) & X
3
(
u
2
).
(a)
(b)
(c)
(d)Figure5(a) shows that at V
ref
=
5
V, X
1
(
i
L
) takes a combined oscillation of large and small amplitudes as the flywheel rotor slows down. At the beginning of discharge, X
1
(
i
L
) shows a quasiperiodic oscillation with a high-value magnitude. Such oscillations decay rapidly when X
5
(
ω
) drops to 0.5; then X
1
(
i
L
) becomes steady. The close-up views are shown in Figures 5(b)-5(c). Dynamics of X
3
(
u
2
) and X
2
(
u
1
) as X
5
(
ω
) decreases are shown in Figure 5(d). For V
ref
=
10
V and V
ref
=
15
V shown in Figures 6-7, similar oscillations with different amplitudes and frequencies happen to X
1
(
i
L
). Although the control parameters are carefully designed according to the closed-loop transfer function of the power drive system, the interaction and difference between the mechanical part and the electrical part still cause complex nonlinear dynamics. It can be seen that larger V
ref determines stronger oscillations and X
5
(
ω
) changes much slower than other state variables. This indicates the existence of two-timescale characteristics in the discharge mode of FESS. Moreover, the oscillations of the electrical variables weaken along with the decrease of the mechanical variables, which will be proved related to the qualitative dynamics in the following section.


## 3. Singular Perturbation Theory-Based Qualitative Dynamics Analysis

In this section, based on the singular perturbation theory, the two-timescale approach is proposed to separate the full-system into the fast and slow subsystems, providing a way for analyzing the interaction of the two-timescale dynamics. Treating the slow variable as constant, stability analysis of the transient fixed points of the full-system is proposed to describe the evolution of dynamics of the full process.


### 3.1. Two-Timescale Model

Averaged model approach is a common method to analyze the physical mechanism of converters which neglects the switching details but focuses on the envelope of the dynamical motion. The following analysis will adopt this approach. Set the duty ratio asd. Set all the derivatives to zero, and we get
(7)
X
1
•
(
I
L
)
=
(
d
-
1
)
a
1
X
2
(
U
1
)
+
d
a
1
X
3
(
U
2
)
=
0
,
X
2
•
(
U
1
)
=
(
d
-
1
)
a
2
X
1
(
I
L
)
+
a
2
X
2
(
U
1
)
=
0
,
X
3
•
(
U
2
)
=
a
3
X
4
(
I
a
)
-
d
a
3
X
1
(
I
L
)
=
0
,
X
4
•
(
I
a
)
=
a
4
X
5
(
Ω
)
+
a
5
X
4
(
I
a
)
+
a
6
X
3
(
U
2
)
=
0
,
X
5
•
(
Ω
)
=
a
7
X
4
(
I
a
)
+
a
8
X
5
(
Ω
)
=
0
,
X
6
•
(
X
1
c
)
=
-
a
9
X
2
(
U
1
)
+
a
9
=
0
,
X
7
•
(
X
2
c
)
=
a
10
X
1
(
I
L
)
+
a
11
X
2
(
U
1
)
+
a
12
X
6
(
X
1
c
)
-
a
11
=
0
,
where X
1
(
I
L
), X
2
(
U
1
), X
3
(
U
2
), X
4
(
I
a
), X
5
(
Ω
), X
6
(
X
1
c
), and X
7
(
X
2
c
) are average values of X
1
(
i
L
), X
2
(
u
1
), X
3
(
u
2
), X
4
(
i
a
), X
5
(
ω
), X
6
(
x
1
c
), and X
7
(
x
2
c
) during a duty cycle, U
ramp is the independent sawtooth peak voltage of the control subsystem, a
13
=
K
i
i
/
U
ramp, a
14
=
K
p
i
K
p
v
h
i
/
U
ramp, a
15
=
-
K
p
i
K
p
v
h
v
h
i
V
ref
/
U
ramp, a
16
=
-
K
p
i
h
v
h
i
V
ref
/
(
R
U
ramp
), and d
=
a
13
X
7
(
X
2
c
)
+
a
14
X
6
(
X
1
c
)
+
a
15
X
2
(
U
1
)
+
a
16
X
1
(
I
L
)
-
a
15. a
i
(
i
=
1,2
,
…
,
12
) are real constants with different magnitudes, among which a
7 and a
8 are about 100 times smaller than the others. Thus X
5
(
ω
) changes much slower than other variables, which proves the existence of state variables with two timescales. The full-system can be divided into fast subsystem and slow subsystem, that is, the fast subsystem α which contains (
X
1
(
i
L
)
,
X
2
(
u
1
)
,
X
3
(
u
2
)
,
X
4
(
i
a
)
,
X
6
(
x
1
c
)
,
X
7
(
x
2
c
)
) and the slow subsystem β which contains X
5
(
ω
). From (7) we can see a
4 reflects the coupling effect between X
5
(
ω
) and X
4
(
i
a
) and also the coupling effect between the fast and slow subsystems. Define a
4 as coupling coefficient; the influence on dynamics of FESS from a
4 will be studied later.


### 3.2. Hopf Bifurcation Set for Fast Subsystem

To study the mechanism of the full-system’s dynamics, the concept of transient fixed point is proposed; for general multitimescale system, define the transient fixed point of the full-system as the fixed point of the fast subsystem under slow variables with different fixed values. During a sufficiently short time period[
τ
k
,
τ
k
+
Δ
τ
], the movement trend of the fast subsystem during [
τ
k
,
τ
k
+
Δ
τ
] can be predicted by the eigenvalues on the basis at X
5
¯
=
X
5
(
ω
)
|
τ
k. Thus the microstructure of the full trajectory can be described. The transient fixed point at various X
5
¯ is(8)
[
X
1
X
2
X
3
X
4
X
6
X
7
]
T
=
[
1
1
-
D
1
1
-
D
D
1
-
D
a
10
(
D
-
1
)
a
12
1
a
13
(
D
+
a
16
D
-
1
+
a
10
a
14
(
1
-
D
)
a
12
)
]
T
,where D
=
(
2
a
6
+
a
4
X
5
¯
-
a
5
+
(
a
5
-
a
4
X
5
¯
-
2
a
6
)
2
-
4
a
6
(
a
4
X
5
¯
+
a
6
)
)
/
2
a
6.By studying the movement of the eigenvalues of the Jacobian under varyingX
5
¯, stability information such as the occurrence of bifurcations can be obtained, which reveals the bifurcation of the transient fixed points follows: (9)
det
⁡
[
λ
I
-
J
(
X
)
]
=
0
,
J
(
X
)
=
[
a
1
a
16
(
X
2
+
X
3
)
D
a
1
-
a
1
+
a
1
a
15
X
2
+
a
1
a
16
X
3
D
a
1
0
a
1
a
14
(
X
2
+
X
3
)
a
1
a
13
(
X
2
+
X
3
)
D
a
2
-
a
2
+
a
2
a
16
X
1
a
2
+
a
2
a
15
X
1
0
0
a
2
a
14
X
1
a
2
a
13
X
1
-
D
a
3
-
a
3
a
16
X
1
-
a
3
a
15
X
1
0
a
3
-
a
3
a
14
X
1
-
a
3
a
13
X
1
0
0
a
6
a
5
0
0
0
-
a
9
0
0
0
0
a
10
a
11
0
0
a
12
0
]
,where λ is the eigenvalue, I is the identity matrix, and J
(
X
) is the Jacobian matrix at [X
1, X
2, X
3, X
4, X
6, X
7]. Parameters are set as Table 3, and loci of the eigenvalues are shown in Figure 8. When all the eigenvalues are in the left half-side of complex plane, the system is stable. When a couple of complex conjugate eigenvalues simultaneously cross the imaginary axis, Hopf bifurcation occurs [20, 21].Eigenvalue loci asX
5
¯ decreases from 0.9 to 0.1 at (a) V
ref
=
5
V, (b) V
ref
=
10
V, and (c) V
ref
=
15
V.
(a)
(b)
(c)From Figures8(a)–8(c), as X
5
¯ changes from 0.9 to 0.1, the eigenvalues (λ
3
,
λ
4
,
λ
5
,
λ
6) stay on the left half-side of complex plane, while a pair of conjugate complex eigenvalues (λ
1
,
λ
2) are firstly on the right half-side and cross the imaginary axis as X
5
¯ reaches 0.4903, 0.4269, and 0.3774, respectively.


### 3.3. Studies on Qualitative Dynamics

It would be imperative to know how the influence on the stability of the transient fixed points from the slow variable is reflected in the full-system. The study of the internal relations between the properties of the transient fixed points and the full-system’s dynamics can predict the occurrence and evolution of nonlinear dynamics of FESS in the discharge mode. Among all the fast variables, the inductor currenti
L serves as a link between the access system and the flywheel motor system; thus, the stability and dynamic characteristics of i
L influence the energy transmission process a lot. In this part we focus mainly on the dynamics of X
1
(
i
L
) against X
5
(
ω
). Parameters are set as Table 3, and obviously in Figure 9, the transient fixed points curves (ES1, ES2, and ES3) are L-form curves with various X
5
(
ω
), which is divided in a stable part (dashed line) and an unstable part (solid line). Hopf bifurcation points (H
1,H
2, andH
3) are the joints of the two parts.Figure 9
Transient fixed points curves (ES1, ES2, and ES3) with the variation of X
5
(
ω
) at V
ref
=
5
V, 10
V, 15
V.Considering the different kinds of equilibria, the stable node represents the quiescent state (QS), which indicates all the variables are at rest or exhibit small amplitude oscillations. The stable limit cycle surrounding the unstable focus represents the spiking state (SP), which indicates variables may behave in large amplitude oscillations [22]. Figures 10(a)–10(f) show the stroboscopic phase trajectory [23] of the full-system as well as the transient fixed points curves under different reference voltages in two and three dimensions. The phase trajectory for V
ref
=
5
V is plotted in Figures 10(a)-10(b), from which we can see thatH
1 divides the full-system’s trajectory into two qualitatively different parts. Hopf bifurcation occurs atH
1. The trajectory starting atA
1 moves along ES1 toH
1, where the trajectory tends to a limit cycle oscillation, the direction of which can be demonstrated by X
5
•
(
ω
) (see the expression in (7)). BetweenH
1 andA
1, the difference between the fast and slow subsystems on timescales causes large amplitude oscillations around ES1 from the beginning of the discharge process, leading the system to SP (fromA
1 toH
1). The repetitive oscillation stays until the trajectory meetsH
1, at which Hopf bifurcation of the transient fixed point takes place. Using the center manifold theory, the curvature coefficient at the Hopf bifurcation point is less than 0, which proved the existence of supercritical bifurcation. Therefore, the amplitudes of the oscillations decrease gradually after the trajectory passing byH
1 and SP settles down to QS (fromH
1 toB
1). The coupling strength of the fast and slow subsystems causes another type of small amplitude oscillations around ES1 in QS. Then the phase trajectory reachesB
1, at which the full-system becomes stable. When the FESS is in QS, the slow subsystem only influences the position of the transient fixed points but does not affect the dynamics of the full-system. For the reason thatH
1 is the unique Hopf bifurcation point to join SP and QS, there is only one state conversion in a discharge cycle. The above process completes one period of the nonrecurrent qualitative dynamics.Stroboscopic phase trajectory of the full-system as well as the transient fixed points curves under differentV
ref in two and three dimensions at (a)-(b) V
ref
=
5
V, (c)-(d) V
ref
=
10
V, and (e)-(f) V
ref
=
15
V.
(a)
(b)
(c)
(d)
(e)
(f)Similar phenomena occur when the reference voltages are 10 V and 15 V. ForV
ref
=
10
V and 15
V in Figures 10(c)–10(f), theH
2 andH
3 still join SP and QS. Comparing with qualitative dynamics at V
ref
=
5
V, the real part of the pair of complex conjugate eigenvalues of the Jacobian at the transient fixed point corresponding toA
2 andA
3 are larger than that corresponding toA
1, which causes much more intense oscillations in SP. Above all, the slow variable modulates the qualitative dynamics by acting essentially as parameters to the full-system. The FESS shifts from SP to QS with the change of the slow variable. These shifts occur when the slow variable crosses the bifurcation point on the transient fixed points curve. Furthermore, from the expression of a
4 we can see V
ref determines the value of a
4 largely, which represents the strength of the coupling between the two subsystems. When V
ref increases, the coupling strength weakens and QS lasts shorter time. The qualitative dynamics with the variation of the coupling strength are shown in Figure 11. It can be found that the duration T of QS decreases quickly as V
ref increases. As is shown in Figure 11, when V
ref
=
5
V, T
=
87.41; when V
ref
=
10
V, T
=
62.74; and when V
ref
=
15
V, T
=
48.14.Time domain waveforms of QS at (a)V
ref
=
5
V, (b) V
ref
=
10
V, and (c) V
ref
=
15
V.
(a)
(b)
(c)


### 3.4. Mechanism Analysis Based on Homotopy Method

According to the analysis above, the transient fixed points curves obtained possesses Hopf bifurcation points at which the full-system can be divided into a stable part and an unstable part, and qualitative dynamics is closely bound up with the properties of the transient fixed points. Therefore, the mechanism analysis of the bifurcation of the transient fixed points can give an intuitive explanation of the origin of complex oscillations of the full-system. A state-to-eigenvalue correspondence can be set up to reveal the physical mechanism of the qualitative dynamics by tracing the changing trend of the eigenvalues. The homotopy method [24–26] can be applied to link eigenvalues of the Jacobian matrix J
(
X
) of a dynamic model to the corresponding state variables through the following homotopy relation:
(10)
H
(
r
)
=
(
1
-
r
)
F
(
X
)
+
r
J
(
X
)
(
0
≤
r
≤
1
)
,
where F
(
X
)
=
diag
⁡
[
J
11
,
J
22
,
J
33
,
J
44
,
J
55
,
J
66
] is the diagonal Jacobian and r is the homotopy parameter. When r varies in the interval [
0
,
1
] and the difference between each two adjacent values is sufficiently small, homotopy method takes the trajectory of eigenvalues of H
(
r
) as a continuous path. Following the paths from r
=
0 to r
=
1, the correspondence between eigenvalues and state variables of the fast subsystem is established. Parameters are set as Table 3, and Figure 12 shows the trace of the sorted eigenvalues by making r as abscissa and the real part of eigenvalues of H
(
r
) as vertical coordinate at V
ref
=
10
V, X
5
¯
=
0.9. More traces are calculated and show the same correspondence, which is shown in Table 4.Table 4
The correspondence between eigenvalues and state variables.
State variables
eigenvalue
X
2
(
u
1
), X
3
(
u
2
)
λ
1, λ
2
X
1
(
i
L
), X
4
(
i
a
)
λ
3, λ
4
X
6
(
x
1
c
), X
7
(
x
2
c
)
λ
5, λ
6Figure 12
Real part of eigenvalues trace through homotopy method for FESS.From the analysis above, Hopf bifurcation occurs whenλ
1 and λ
2 cross the imaginary axis. Thus λ
1 and λ
2 are the key factors which dominate the stability of the FESS. Considering the correspondence between the eigenvalues and the state variables, we can conclude that the system instability derives from the voltage instability; that is, the voltage instability is the root inducement of the qualitative dynamics. Therefore, dynamic voltage stability control scheme plays a key role in stabilizing the FESS, which will be the focus of our further work.A detailed look into the influencing factors on the qualitative dynamics is taken further by presenting the boundaries of SP and QS in terms ofu
2 and u
1. It shows that the region of QS gets smaller with the increase of u
1, which represents the load level (Figure 13).Figure 13
Critical boundary between SP and QS withu
2
-
u
1.


## 3.1. Two-Timescale Model

Averaged model approach is a common method to analyze the physical mechanism of converters which neglects the switching details but focuses on the envelope of the dynamical motion. The following analysis will adopt this approach. Set the duty ratio asd. Set all the derivatives to zero, and we get
(7)
X
1
•
(
I
L
)
=
(
d
-
1
)
a
1
X
2
(
U
1
)
+
d
a
1
X
3
(
U
2
)
=
0
,
X
2
•
(
U
1
)
=
(
d
-
1
)
a
2
X
1
(
I
L
)
+
a
2
X
2
(
U
1
)
=
0
,
X
3
•
(
U
2
)
=
a
3
X
4
(
I
a
)
-
d
a
3
X
1
(
I
L
)
=
0
,
X
4
•
(
I
a
)
=
a
4
X
5
(
Ω
)
+
a
5
X
4
(
I
a
)
+
a
6
X
3
(
U
2
)
=
0
,
X
5
•
(
Ω
)
=
a
7
X
4
(
I
a
)
+
a
8
X
5
(
Ω
)
=
0
,
X
6
•
(
X
1
c
)
=
-
a
9
X
2
(
U
1
)
+
a
9
=
0
,
X
7
•
(
X
2
c
)
=
a
10
X
1
(
I
L
)
+
a
11
X
2
(
U
1
)
+
a
12
X
6
(
X
1
c
)
-
a
11
=
0
,
where X
1
(
I
L
), X
2
(
U
1
), X
3
(
U
2
), X
4
(
I
a
), X
5
(
Ω
), X
6
(
X
1
c
), and X
7
(
X
2
c
) are average values of X
1
(
i
L
), X
2
(
u
1
), X
3
(
u
2
), X
4
(
i
a
), X
5
(
ω
), X
6
(
x
1
c
), and X
7
(
x
2
c
) during a duty cycle, U
ramp is the independent sawtooth peak voltage of the control subsystem, a
13
=
K
i
i
/
U
ramp, a
14
=
K
p
i
K
p
v
h
i
/
U
ramp, a
15
=
-
K
p
i
K
p
v
h
v
h
i
V
ref
/
U
ramp, a
16
=
-
K
p
i
h
v
h
i
V
ref
/
(
R
U
ramp
), and d
=
a
13
X
7
(
X
2
c
)
+
a
14
X
6
(
X
1
c
)
+
a
15
X
2
(
U
1
)
+
a
16
X
1
(
I
L
)
-
a
15. a
i
(
i
=
1,2
,
…
,
12
) are real constants with different magnitudes, among which a
7 and a
8 are about 100 times smaller than the others. Thus X
5
(
ω
) changes much slower than other variables, which proves the existence of state variables with two timescales. The full-system can be divided into fast subsystem and slow subsystem, that is, the fast subsystem α which contains (
X
1
(
i
L
)
,
X
2
(
u
1
)
,
X
3
(
u
2
)
,
X
4
(
i
a
)
,
X
6
(
x
1
c
)
,
X
7
(
x
2
c
)
) and the slow subsystem β which contains X
5
(
ω
). From (7) we can see a
4 reflects the coupling effect between X
5
(
ω
) and X
4
(
i
a
) and also the coupling effect between the fast and slow subsystems. Define a
4 as coupling coefficient; the influence on dynamics of FESS from a
4 will be studied later.


## 3.2. Hopf Bifurcation Set for Fast Subsystem

To study the mechanism of the full-system’s dynamics, the concept of transient fixed point is proposed; for general multitimescale system, define the transient fixed point of the full-system as the fixed point of the fast subsystem under slow variables with different fixed values. During a sufficiently short time period[
τ
k
,
τ
k
+
Δ
τ
], the movement trend of the fast subsystem during [
τ
k
,
τ
k
+
Δ
τ
] can be predicted by the eigenvalues on the basis at X
5
¯
=
X
5
(
ω
)
|
τ
k. Thus the microstructure of the full trajectory can be described. The transient fixed point at various X
5
¯ is(8)
[
X
1
X
2
X
3
X
4
X
6
X
7
]
T
=
[
1
1
-
D
1
1
-
D
D
1
-
D
a
10
(
D
-
1
)
a
12
1
a
13
(
D
+
a
16
D
-
1
+
a
10
a
14
(
1
-
D
)
a
12
)
]
T
,where D
=
(
2
a
6
+
a
4
X
5
¯
-
a
5
+
(
a
5
-
a
4
X
5
¯
-
2
a
6
)
2
-
4
a
6
(
a
4
X
5
¯
+
a
6
)
)
/
2
a
6.By studying the movement of the eigenvalues of the Jacobian under varyingX
5
¯, stability information such as the occurrence of bifurcations can be obtained, which reveals the bifurcation of the transient fixed points follows: (9)
det
⁡
[
λ
I
-
J
(
X
)
]
=
0
,
J
(
X
)
=
[
a
1
a
16
(
X
2
+
X
3
)
D
a
1
-
a
1
+
a
1
a
15
X
2
+
a
1
a
16
X
3
D
a
1
0
a
1
a
14
(
X
2
+
X
3
)
a
1
a
13
(
X
2
+
X
3
)
D
a
2
-
a
2
+
a
2
a
16
X
1
a
2
+
a
2
a
15
X
1
0
0
a
2
a
14
X
1
a
2
a
13
X
1
-
D
a
3
-
a
3
a
16
X
1
-
a
3
a
15
X
1
0
a
3
-
a
3
a
14
X
1
-
a
3
a
13
X
1
0
0
a
6
a
5
0
0
0
-
a
9
0
0
0
0
a
10
a
11
0
0
a
12
0
]
,where λ is the eigenvalue, I is the identity matrix, and J
(
X
) is the Jacobian matrix at [X
1, X
2, X
3, X
4, X
6, X
7]. Parameters are set as Table 3, and loci of the eigenvalues are shown in Figure 8. When all the eigenvalues are in the left half-side of complex plane, the system is stable. When a couple of complex conjugate eigenvalues simultaneously cross the imaginary axis, Hopf bifurcation occurs [20, 21].Eigenvalue loci asX
5
¯ decreases from 0.9 to 0.1 at (a) V
ref
=
5
V, (b) V
ref
=
10
V, and (c) V
ref
=
15
V.
(a)
(b)
(c)From Figures8(a)–8(c), as X
5
¯ changes from 0.9 to 0.1, the eigenvalues (λ
3
,
λ
4
,
λ
5
,
λ
6) stay on the left half-side of complex plane, while a pair of conjugate complex eigenvalues (λ
1
,
λ
2) are firstly on the right half-side and cross the imaginary axis as X
5
¯ reaches 0.4903, 0.4269, and 0.3774, respectively.


## 3.3. Studies on Qualitative Dynamics

It would be imperative to know how the influence on the stability of the transient fixed points from the slow variable is reflected in the full-system. The study of the internal relations between the properties of the transient fixed points and the full-system’s dynamics can predict the occurrence and evolution of nonlinear dynamics of FESS in the discharge mode. Among all the fast variables, the inductor currenti
L serves as a link between the access system and the flywheel motor system; thus, the stability and dynamic characteristics of i
L influence the energy transmission process a lot. In this part we focus mainly on the dynamics of X
1
(
i
L
) against X
5
(
ω
). Parameters are set as Table 3, and obviously in Figure 9, the transient fixed points curves (ES1, ES2, and ES3) are L-form curves with various X
5
(
ω
), which is divided in a stable part (dashed line) and an unstable part (solid line). Hopf bifurcation points (H
1,H
2, andH
3) are the joints of the two parts.Figure 9
Transient fixed points curves (ES1, ES2, and ES3) with the variation of X
5
(
ω
) at V
ref
=
5
V, 10
V, 15
V.Considering the different kinds of equilibria, the stable node represents the quiescent state (QS), which indicates all the variables are at rest or exhibit small amplitude oscillations. The stable limit cycle surrounding the unstable focus represents the spiking state (SP), which indicates variables may behave in large amplitude oscillations [22]. Figures 10(a)–10(f) show the stroboscopic phase trajectory [23] of the full-system as well as the transient fixed points curves under different reference voltages in two and three dimensions. The phase trajectory for V
ref
=
5
V is plotted in Figures 10(a)-10(b), from which we can see thatH
1 divides the full-system’s trajectory into two qualitatively different parts. Hopf bifurcation occurs atH
1. The trajectory starting atA
1 moves along ES1 toH
1, where the trajectory tends to a limit cycle oscillation, the direction of which can be demonstrated by X
5
•
(
ω
) (see the expression in (7)). BetweenH
1 andA
1, the difference between the fast and slow subsystems on timescales causes large amplitude oscillations around ES1 from the beginning of the discharge process, leading the system to SP (fromA
1 toH
1). The repetitive oscillation stays until the trajectory meetsH
1, at which Hopf bifurcation of the transient fixed point takes place. Using the center manifold theory, the curvature coefficient at the Hopf bifurcation point is less than 0, which proved the existence of supercritical bifurcation. Therefore, the amplitudes of the oscillations decrease gradually after the trajectory passing byH
1 and SP settles down to QS (fromH
1 toB
1). The coupling strength of the fast and slow subsystems causes another type of small amplitude oscillations around ES1 in QS. Then the phase trajectory reachesB
1, at which the full-system becomes stable. When the FESS is in QS, the slow subsystem only influences the position of the transient fixed points but does not affect the dynamics of the full-system. For the reason thatH
1 is the unique Hopf bifurcation point to join SP and QS, there is only one state conversion in a discharge cycle. The above process completes one period of the nonrecurrent qualitative dynamics.Stroboscopic phase trajectory of the full-system as well as the transient fixed points curves under differentV
ref in two and three dimensions at (a)-(b) V
ref
=
5
V, (c)-(d) V
ref
=
10
V, and (e)-(f) V
ref
=
15
V.
(a)
(b)
(c)
(d)
(e)
(f)Similar phenomena occur when the reference voltages are 10 V and 15 V. ForV
ref
=
10
V and 15
V in Figures 10(c)–10(f), theH
2 andH
3 still join SP and QS. Comparing with qualitative dynamics at V
ref
=
5
V, the real part of the pair of complex conjugate eigenvalues of the Jacobian at the transient fixed point corresponding toA
2 andA
3 are larger than that corresponding toA
1, which causes much more intense oscillations in SP. Above all, the slow variable modulates the qualitative dynamics by acting essentially as parameters to the full-system. The FESS shifts from SP to QS with the change of the slow variable. These shifts occur when the slow variable crosses the bifurcation point on the transient fixed points curve. Furthermore, from the expression of a
4 we can see V
ref determines the value of a
4 largely, which represents the strength of the coupling between the two subsystems. When V
ref increases, the coupling strength weakens and QS lasts shorter time. The qualitative dynamics with the variation of the coupling strength are shown in Figure 11. It can be found that the duration T of QS decreases quickly as V
ref increases. As is shown in Figure 11, when V
ref
=
5
V, T
=
87.41; when V
ref
=
10
V, T
=
62.74; and when V
ref
=
15
V, T
=
48.14.Time domain waveforms of QS at (a)V
ref
=
5
V, (b) V
ref
=
10
V, and (c) V
ref
=
15
V.
(a)
(b)
(c)


## 3.4. Mechanism Analysis Based on Homotopy Method

According to the analysis above, the transient fixed points curves obtained possesses Hopf bifurcation points at which the full-system can be divided into a stable part and an unstable part, and qualitative dynamics is closely bound up with the properties of the transient fixed points. Therefore, the mechanism analysis of the bifurcation of the transient fixed points can give an intuitive explanation of the origin of complex oscillations of the full-system. A state-to-eigenvalue correspondence can be set up to reveal the physical mechanism of the qualitative dynamics by tracing the changing trend of the eigenvalues. The homotopy method [24–26] can be applied to link eigenvalues of the Jacobian matrix J
(
X
) of a dynamic model to the corresponding state variables through the following homotopy relation:
(10)
H
(
r
)
=
(
1
-
r
)
F
(
X
)
+
r
J
(
X
)
(
0
≤
r
≤
1
)
,
where F
(
X
)
=
diag
⁡
[
J
11
,
J
22
,
J
33
,
J
44
,
J
55
,
J
66
] is the diagonal Jacobian and r is the homotopy parameter. When r varies in the interval [
0
,
1
] and the difference between each two adjacent values is sufficiently small, homotopy method takes the trajectory of eigenvalues of H
(
r
) as a continuous path. Following the paths from r
=
0 to r
=
1, the correspondence between eigenvalues and state variables of the fast subsystem is established. Parameters are set as Table 3, and Figure 12 shows the trace of the sorted eigenvalues by making r as abscissa and the real part of eigenvalues of H
(
r
) as vertical coordinate at V
ref
=
10
V, X
5
¯
=
0.9. More traces are calculated and show the same correspondence, which is shown in Table 4.Table 4
The correspondence between eigenvalues and state variables.
State variables
eigenvalue
X
2
(
u
1
), X
3
(
u
2
)
λ
1, λ
2
X
1
(
i
L
), X
4
(
i
a
)
λ
3, λ
4
X
6
(
x
1
c
), X
7
(
x
2
c
)
λ
5, λ
6Figure 12
Real part of eigenvalues trace through homotopy method for FESS.From the analysis above, Hopf bifurcation occurs whenλ
1 and λ
2 cross the imaginary axis. Thus λ
1 and λ
2 are the key factors which dominate the stability of the FESS. Considering the correspondence between the eigenvalues and the state variables, we can conclude that the system instability derives from the voltage instability; that is, the voltage instability is the root inducement of the qualitative dynamics. Therefore, dynamic voltage stability control scheme plays a key role in stabilizing the FESS, which will be the focus of our further work.A detailed look into the influencing factors on the qualitative dynamics is taken further by presenting the boundaries of SP and QS in terms ofu
2 and u
1. It shows that the region of QS gets smaller with the increase of u
1, which represents the load level (Figure 13).Figure 13
Critical boundary between SP and QS withu
2
-
u
1.


## 4. Application

### 4.1. Feasibility Regions Analysis

This part will apply the two-timescale approach to derive the feasibility regions of FESS in the discharge mode and then provide instructions to parameters setting of FESS. From the analysis above we know that the feasibility regions of the fast subsystem dominate that of the full-system. Therefore, the FESS is stable when all the eigenvalues of the fast subsystem are in the left half-side of complex plane withX
5
(
ω
) decreasing from 1 to 0.The feasibility regions of key parameters are shown in Figure14. Figure 14(a) shows the feasibility boundary in the parameter space of V
ref versus the outer loop control parameters K
p
v and K
i
v and Figure 14(b) shows the feasibility boundary in the parameter space of V
ref versus the inner loop control parameters K
p
i and K
i
i, all of which clearly illustrate the effect of those sensitive parameters on the feasibility regions. The space in front of the critical surface corresponds to stable operation and the space behind it corresponds to unstable operation. The results can be used as instructions to the parameters setting of the access unit of FESS itself and constraints to improve the safety and stability of FESS and the power system.Feasibility regions of key parameters, (a)V
ref versus K
p
v and K
i
v at K
p
i
=
40, K
i
i
=
36, (b) V
ref versus K
p
i and K
i
i at K
p
v
=
60, K
i
v
=
49.
(a)
(b)


### 4.2. Application Requirements

Obviously, the applicability and rationality of the proposed two-timescale approach with transient fixed points analysis mainly depend on the existence of state variables with two timescales, which is not to be considered as an exact criterion but as a guideline [27]. From model (7) we can see a
7
=
-
K
m
V
ref
L
C
1
/
(
J
R
ω
ref
) and a
8
=
-
B
v
L
C
1
/
J represent the change rates of X
5
(
ω
); thus, J dominates the difference between the fast and slow variables on timescales. For typical electromechanical coupling systems, the two-timescale characteristics are ubiquitous but in degree. As is shown in the foregoing analysis, the proposed approach does well in predicting the qualitative dynamics when the magnitude difference of fast and slow variables is about 100 times. Set J to 0.5, 0.25, and 0.1 times of its original value; other parameters are set as Table 3; then the corresponding magnitude difference will be 50, 25, and 10 times. Figure 15 shows the stroboscopic phase trajectory of the full-system in the phase plane X
5
(
ω
)
-
X
1
(
i
L
) as well as the transient fixed points bifurcation diagram at V
ref
=
10
V. At J
′
=
0.5
·
J in Figure 15(a),H
4 divides the phase trajectory of the full-system into two qualitatively different parts, which agrees very well with ES4. At J
′
=
0.25
·
J in Figure 15(b), ES5 basically corresponds to the dynamics of the full-system. At J
′
=
0.1
·
J in Figure 15(c), ES6 has deviated from the correct equilibrium position.Stroboscopic phase trajectory of the full-system inX
5
(
ω
)
-
X
1
(
i
L
) as well as the transient fixed points curves at (a) J
′
=
0.5
·
J, (b) J
′
=
0.25
·
J, and (c) J
′
=
0.1
·
J.H
4,H
5, andH
6-Hopf bifurcation point.
(a)
(b)
(c)Considering that the two-timescale approach based on singular perturbation theory is a kind of model reduction method, Ghorbel and Spong [28] have given out the condition for the reduction of multitimescales system model; the equilibrium of the fast system must be close to that of the full-system. We can see that though the SP and QS caused by Hopf bifurcation are especially apparent in qualitative dynamics, they can still become weakening or even disappear as the difference of state variables on timescales diminishes. And meanwhile, the deviation between the transient fixed points curves and the stroboscopic phase trajectory of the full-system becomes larger, which reflects the inapplicability of the proposed approach. Therefore, the applicability and rationality of the proposed approach we concern here mainly refers to the bifurcation characteristics of the fast subsystem can actually reflect what extent of the full-system’s dynamical evolution. From the numerical simulations in Figure 15 we can see the proposed approach is applicable when the magnitude difference of state variables is bigger than 25 times; more simulations have been done and gave out the same conclusion that when the magnitude difference is not big enough, the proposed approach is not applicable.


## 4.1. Feasibility Regions Analysis

This part will apply the two-timescale approach to derive the feasibility regions of FESS in the discharge mode and then provide instructions to parameters setting of FESS. From the analysis above we know that the feasibility regions of the fast subsystem dominate that of the full-system. Therefore, the FESS is stable when all the eigenvalues of the fast subsystem are in the left half-side of complex plane withX
5
(
ω
) decreasing from 1 to 0.The feasibility regions of key parameters are shown in Figure14. Figure 14(a) shows the feasibility boundary in the parameter space of V
ref versus the outer loop control parameters K
p
v and K
i
v and Figure 14(b) shows the feasibility boundary in the parameter space of V
ref versus the inner loop control parameters K
p
i and K
i
i, all of which clearly illustrate the effect of those sensitive parameters on the feasibility regions. The space in front of the critical surface corresponds to stable operation and the space behind it corresponds to unstable operation. The results can be used as instructions to the parameters setting of the access unit of FESS itself and constraints to improve the safety and stability of FESS and the power system.Feasibility regions of key parameters, (a)V
ref versus K
p
v and K
i
v at K
p
i
=
40, K
i
i
=
36, (b) V
ref versus K
p
i and K
i
i at K
p
v
=
60, K
i
v
=
49.
(a)
(b)


## 4.2. Application Requirements

Obviously, the applicability and rationality of the proposed two-timescale approach with transient fixed points analysis mainly depend on the existence of state variables with two timescales, which is not to be considered as an exact criterion but as a guideline [27]. From model (7) we can see a
7
=
-
K
m
V
ref
L
C
1
/
(
J
R
ω
ref
) and a
8
=
-
B
v
L
C
1
/
J represent the change rates of X
5
(
ω
); thus, J dominates the difference between the fast and slow variables on timescales. For typical electromechanical coupling systems, the two-timescale characteristics are ubiquitous but in degree. As is shown in the foregoing analysis, the proposed approach does well in predicting the qualitative dynamics when the magnitude difference of fast and slow variables is about 100 times. Set J to 0.5, 0.25, and 0.1 times of its original value; other parameters are set as Table 3; then the corresponding magnitude difference will be 50, 25, and 10 times. Figure 15 shows the stroboscopic phase trajectory of the full-system in the phase plane X
5
(
ω
)
-
X
1
(
i
L
) as well as the transient fixed points bifurcation diagram at V
ref
=
10
V. At J
′
=
0.5
·
J in Figure 15(a),H
4 divides the phase trajectory of the full-system into two qualitatively different parts, which agrees very well with ES4. At J
′
=
0.25
·
J in Figure 15(b), ES5 basically corresponds to the dynamics of the full-system. At J
′
=
0.1
·
J in Figure 15(c), ES6 has deviated from the correct equilibrium position.Stroboscopic phase trajectory of the full-system inX
5
(
ω
)
-
X
1
(
i
L
) as well as the transient fixed points curves at (a) J
′
=
0.5
·
J, (b) J
′
=
0.25
·
J, and (c) J
′
=
0.1
·
J.H
4,H
5, andH
6-Hopf bifurcation point.
(a)
(b)
(c)Considering that the two-timescale approach based on singular perturbation theory is a kind of model reduction method, Ghorbel and Spong [28] have given out the condition for the reduction of multitimescales system model; the equilibrium of the fast system must be close to that of the full-system. We can see that though the SP and QS caused by Hopf bifurcation are especially apparent in qualitative dynamics, they can still become weakening or even disappear as the difference of state variables on timescales diminishes. And meanwhile, the deviation between the transient fixed points curves and the stroboscopic phase trajectory of the full-system becomes larger, which reflects the inapplicability of the proposed approach. Therefore, the applicability and rationality of the proposed approach we concern here mainly refers to the bifurcation characteristics of the fast subsystem can actually reflect what extent of the full-system’s dynamical evolution. From the numerical simulations in Figure 15 we can see the proposed approach is applicable when the magnitude difference of state variables is bigger than 25 times; more simulations have been done and gave out the same conclusion that when the magnitude difference is not big enough, the proposed approach is not applicable.


## 5. Experimental Verification

To verify the analysis in this paper, constant voltage discharge experiment is carried out; the parameter values are set as those in Table3. The driving motor of the flywheel was solved by BLDCM, which indicates high reliability and the rotational speed being up to 8000 r/min. First charge up the FESS to 60% of its rated speed; then catch the time-variant dynamics under discharge mode at different V
ref. Set V
ref
=
15
V; the trajectory of i
L evolves with ω; when X
5
(
ω
)
=
0.60 (60% of rated speed) in Figure 16(a), i
L shows a large-scale oscillation with a frequency of 666.7 Hz and an amplitude of 4 A, implying the FESS is in SP; when X
5
(
ω
)
=
0.37 in Figure 16(b), the amplitudes of quasiperiodic oscillation behaviors decrease gradually and approach a nearly stable limit cycle, as the property by Hopf bifurcation of the fast subsystem. The waveform of i
L is superposed with a sinusoidal oscillation with a frequency of 625 Hz and an amplitude of 2.4 A; then it settles down to QS; when X
5
(
ω
)
=
0.25 and X
5
(
ω
)
=
0.20 in Figures 16(c)-16(d), the phase trajectory of i
L
-
u
1 is a point and the FESS operates steady.Discharge experiment forV
ref
=
15
V, (a) X
5
(
ω
)
=
0.60, (b) X
5
(
ω
)
=
0.37, (c) X
5
(
ω
)
=
0.25, and (d) X
5
(
ω
)
=
0.20.
(a)
(b)
(c)
(d)The structures of qualitative dynamics change withV
ref. Figures 17-18 present phase portraits for V
ref
=
10
V and V
ref
=
5
V, respectively. Set V
ref
=
10
V; when X
5
(
ω
)
=
0.60 in Figure 17(a), i
L shows a large-scale oscillation with a frequency of 625 Hz and an amplitude of 4 A; when X
5
(
ω
)
=
0.42 in Figure 17(b), the waveform of i
L mixes with a sinusoidal oscillation with a frequency of 588 Hz and an amplitude of 2 A; when X
5
(
ω
)
=
0.25 and X
5
(
ω
)
=
0.20 in Figures 17(c)-17(d), the phase trajectory of i
L
-
u
1 is a point and the FESS operates steady. Set V
ref
=
5
V; when X
5
(
ω
)
=
0.60 in Figure 18(a), i
L shows a large-scale oscillation with a frequency of 600 Hz and an amplitude of 3.8 A; when X
5
(
ω
)
=
0.49 in Figure 18(b), the waveform of i
L mixes with a sinusoidal oscillation with a frequency of 90 Hz and an amplitude of 0.67 A; when X
5
(
ω
)
=
0.25 and X
5
(
ω
)
=
0.20 in Figures 18(c)-18(d), the phase trajectory of i
L
-
u
1 is a point and the FESS operates steady.Discharge experiment forV
ref
=
10
V, (a) X
5
(
ω
)
=
0.60, (b) X
5
(
ω
)
=
0.42, (c) X
5
(
ω
)
=
0.25, and (d) X
5
(
ω
)
=
0.20.
(a)
(b)
(c)
(d)Discharge experiment forV
ref
=
5
V, (a) X
5
(
ω
)
=
0.60, (b) X
5
(
ω
)
=
0.49, (c) X
5
(
ω
)
=
0.25, and (d) X
5
(
ω
)
=
0.20.
(a)
(b)
(c)
(d)


## 6. Conclusion

This paper investigates the qualitative dynamics of the voltage-current dual-loop controlled FESS, which is mainly shown as the fast oscillations of the inductor current and the motor voltage weakens along with the slowdown of the flywheel rotor. By the proposed two-timescale approach based on singular perturbation theory, the state variables are separated into fast and slow variables. First, it is shown that the stability of FESS is closely bound up with that of the transient fixed points. The FESS shifts from SP to QS when the slow variable crosses the bifurcation point on the transient fixed points curve. Larger eigenvalues’ real parts determine more intense oscillations in SP, and larger coupling coefficient leads to longer duration of QS. Further analysis shows that the evolution of the full-system’s dynamics is dominated by the difference between the slow and fast variables on timescales, whereas the qualitative dynamics are mainly caused by the voltage instability. Moreover, the feasibility regions of the main system parameters are derived, in which stability operation and power transmission quality of FESS can be guaranteed. Finally, an applicability investigation shows that when the difference of state variables on timescales is not big enough, the proposed approach is not applicable. This paper provides insights into the effect of two-timescale characteristics on the safety and stability of energy transmission of FESS. The results can be used as instructions to the parameters setting of the FESS itself and constraints to improve the safety and stability of FESS and the smart grids.



---

*Source: 487596-2014-08-03.xml*