# A Method to Determine Oscillation Emergence Bifurcation in Time-Delayed LTI System with Single Lag

Authors: Yu Xiaodan; Jia Hongjie; Wang Chengshan; Jiang Yilang
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

One type of bifurcation named oscillation emergence bifurcation (OEB) found in time-delayed linear time invariant (abbr. LTI) systems is fully studied. The definition of OEB is initially put forward according to the eigenvalue variation. It is revealed that a real eigenvalue splits into a pair of conjugated complex eigenvalues when an OEB occurs, which means the number of the system eigenvalues will increase by one and a new oscillation mode will emerge. Next, a method to determine OEB bifurcation in the time-delayed LTI system with single lag is developed based on Lambert W function. A one-dimensional (1-dim) time-delayed system is firstly employed to explain the mechanism of OEB bifurcation. Then, methods to determine the OEB bifurcation in 1-dim, 2-dim, and high-dimension time-delayed LTI systems are derived. Finally, simulation results validate the correctness and effectiveness of the presented method. Since OEB bifurcation occurs with a new oscillation mode emerging, work of this paper is useful to explore the complex phenomena and the stability of time-delayed dynamic systems.


---

## Body

## 1. Introduction

In the nature, system trend is not only dependent on the current operating point but also subject to its previous conditions. This phenomenon is called time delay [1, 2]. As a main cause to the control deterioration and system instability, time delays exist widespread in the industry of communication [3], biology [4], mechanism [5], chemistry [6], society [7], and power utility [8–10]. Therefore, it is an important issue and frequent topic to evaluate the impact of time delays on the stability of dynamic systems.Many methods to analyze the stability of time-delayed dynamic system have been proposed in the past literatures. For example, [10–13] applied Lyapunov theory to derive various stability criteria for time-delayed systems. To seek the suitable Lyapunov function (or Lyapunov functional) and to construct the corresponding criteria were the two key points of such methods. Linear matrix inequality (LMI) technique, Lyapunov-Razumikhin theory, and Lyapunov-Krasovskii theory were usually used in the derivation. Just as we know, characteristic equation of time-delayed system has transcendental element, which is very difficult to handle in the analysis. In order to surmount this problem, [14, 15] utilized Rekasius substitution method and [16] used Pade approximation method to approximate the transcendental elements and transfer the characteristic equation into a polynomial, which was used to calculate the system critical eigenvalues. Reference [9, 17] treated the stability region method as a tool to determine the small signal stability region of time-delayed dynamic system. All points in the stability region were small signal stable, and the boundary of the stability region was composed by some bifurcation points, such as saddle node bifurcation, Hopf bifurcation, and singularity induced bifurcation. Reference [18–22] applied bifurcation theory to analyze some complex phenomena existing in the time-delayed systems, such as Hopf bifurcation, period doubling bifurcation, and chaos.Time-delayed system is a complex dynamic system, in which many interesting and intricate dynamic behaviors exist. For example, Figure1 from [23] demonstrated that a real eigenvalue split into a pair of conjugated complex eigenvalues at point P with increase of time delay in a time-delayed power system. In this paper, this phenomenon is called oscillation emergence bifurcation (OEB). And, we apply the theory of bifurcation analysis and Lambert W function to propose a mechanism to the occurrence of this bifurcation. As an available tool, Lambert W function [24, 25] has been utilized to analyze some complex behaviors in nonlinear system including stability analysis of time-delayed system [26–29]. It is applied to derive a general method to determine OEB bifurcation in this paper.Figure 1
A real eigenvalue splits into a pair of conjugated eigenvalues.The contents of this paper are organized into six sections. Following this introduction, definition of OEB bifurcation is put forward in Section2 and some of its features are also discussed in this section. In Section 3, Lambert W function is introduced and a mechanism of the OEB’s occurrence is proposed. Section 4 presents a method to determine OEB bifurcation in time-delayed LTI system with single lag based on the Lambert W function. In Section 5, the presented approach is demonstrated by numerical simulations. Finally, Section 6 concludes and summarizes this paper.


## 2. Definition of OEB Bifurcation

In this paper, all studies are based on the following time-delayed LTI system with single lag:(1)
Σ
:
{
x
˙
=
A
0
·
x
+
A
i
·
x
(
t
-
τ
)
x
(
t
)
=
ϕ
(
t
)
,
t
∈
[
-
τ
,
0
]
,
where A
0
,
A
i
∈
R
n
×
n; x
∈
R
n is the vector of the state variables; τ
≥
0 is the constant of time delay; ϕ
(
t
) is the system trajectory in the range of [
-
τ
,
0
].It is well known that the following equation is the corresponding characteristic equation of systemΣ:
(2)
det
⁡
(
s
-
A
0
-
A
i
·
e
-
s
·
τ
)
=
0
,
where s
∈
C is the eigenvalue to be determined. The eigenvalue spectrum of system Σ is denoted by Γ:
(3)
Γ
=
[
s
1
,
s
2
,
s
3
,
…
]
.
Theoretically, it can be calculated from (2) directly when τ is given. Further, denote the set of the real eigenvalue by Γ
R and the set of the complex eigenvalue with nonzero imaginary by Γ
I:
(4)
Γ
R
=
[
s
1
,
s
2
,
…
,
s
n
R
]
,
Γ
I
=
[
s
~
1
,
s
~
2
,
…
,
s
~
n
I
]
,
where n
R, n
I are the eigenvalue numbers of Γ
R and Γ
I, respectively. Then the following (5) holds:
(5)
Γ
=
[
Γ
R
,
Γ
I
]
,
n
τ
=
n
R
+
n
I
,
where n
τ is the eigenvalue number of Γ.It can be found that whenτ
=
0, system Σ degenerates into the following LTI differential equation:
(6)
x
˙
=
A
·
x
,
where A
=
A
0
+
A
i. And, the number of eigenvalues is n; that is, n
τ
=
0
=
n. When τ increases, the values of n
R, n
I, and n
τ will change accordingly. Therefore, they can be considered as functions of τ as follows:
(7)
n
R
=
N
R
(
τ
)
,
n
I
=
N
I
(
τ
)
,
n
τ
=
N
(
τ
)
=
N
R
(
τ
)
+
N
I
(
τ
)
.
N
R
(
τ
), N
I
(
τ
), and N
(
τ
) are three functions of τ to yield the numbers of real eigenvalues, complex eigenvalues, and whole eigenvalues of system Σ.Definition 1 (definition of oscillation emergence bifurcation).
Presume that there is a value of time delayτ
c
>
0 for system Σ of (1). Its Γ
R has n
R
o real eigenvalues and Γ
I has n
I
o complex eigenvalues with τ
≤
τ
c. And, after τ
>
τ
c, the numbers of Γ
R and Γ
I turn to n
R
n and n
I
n. If the following equations hold:
(8)
n
R
n
=
n
R
o
-
1
,
n
I
n
=
n
I
o
+
2
,
n
τ
n
=
n
τ
o
+
1
,
n
τ
n
=
n
R
n
+
n
I
n
,
n
τ
o
=
n
R
o
+
n
I
o
,there is an oscillation emergence bifurcation at τ
=
τ
c.According to Definition1, suppose there is a real eigenvalue s
∈
R in Γ at τ
≤
τ
c. It will split into a pair of conjugated eigenvalues s
+, s
- as shown in Figure 2. It is manifested that the number of Γ
R will decrease by one, the number of Γ
I will increase by two, and number of Γ will increase by one after an OEB bifurcation occurs. At the same time, a new oscillation mode subject to s
+, s
- will emerge after OEB bifurcation.Figure 2
Illustration of OEB bifurcation.


## 3. LambertW Function and OEB Bifurcation

Our aim is to derive a general method based on LambertW function to determine OEB bifurcation in system Σ defined by (1). A brief introduction to Lambert W function is firstly given.The following function is called LambertW function:
(9)
w
e
w
=
z
,
where w
,
z
∈
C. Its inverse function is expressed as
(10)
w
=
W
(
z
)
,
that
is
,
W
(
z
)
e
W
(
z
)
=
z
.LambertW function is not injective and the relation W is multivalued (except at 0). It has many solution branches in complex plane [25], which can be denoted by W
k, k
=
0
,
±
1
,
±
2
,
…. If the attention is restricted to real-valued W, then the relation is defined only for z
≥
-
1
/
e. For -
1
/
e
≤
z
<
0, there are two possible real values of W
(
z
) as shown in Figure 3. We denote the branch satisfying -
1
≤
W
(
z
) by W
0 (solid line in Figure 3), and the branch satisfying W
(
z
)
≤
-
1 by W
-
1 (dashed line in Figure 3). W
0 is referred to the principal branch of Lambert W function. It increases from W
0
(
-
1
/
e
)
=
-
1 to W
0
(
+
∞
)
=
+
∞, while branch W
-
1 decreases from W
-
1
(
-
1
/
e
)
=
-
1 to W
-
1
(
0
-
)
=
-
∞. They encounter at point P as shown in Figure 3, and the value of point P is given as follows:
(11)
P
:
(
z
,
w
)
=
(
-
1
e
,
-
1
)
.Figure 3
W
0 and W
-
1 branches of Lambert W function.Supposed thatz
∈
R in (9) and (10) is a system parameter and it varies with another parameter u
∈
R:
(12)
z
=
f
z
(
u
)
,
z
∈
R
.Substitute (12) into (10). It can be found that w is also a function of u:
(13)
w
=
W
(
z
)
=
W
(
f
z
(
u
)
)
=
W
z
(
u
)
.Then, theu
-
z curve can be plotted. An example is shown in Figure 4, where u
-
z curve crosses the dotted line of z
=
-
1
/
e at some points which can be divided into two categories:
(14)
M
s
=
(
M
1
,
M
2
,
…
,
M
i
,
…
)
,
N
s
=
(
N
1
,
N
2
,
…
,
N
i
,
…
)
.Figure 4
Curve ofz
=
f
z
(
u
).Ifu is on the interval of (
M
i
,
N
i
), that is, the intervals above the line with z
>
-
1
/
e, (13) has real solutions. While if u is on the interval of (
N
i
,
M
i
+
1
), that is, the grey intervals with z
<
-
1
/
e, (13) only has complex solutions with nonzero imaginary.Recalling the definition of OEB, we can find thatN
s is just a set of OEB bifurcations of the system expressed by (13). That is, at any point of N
s, a real eigenvalue will turn into a pair of conjugated complex eigenvalues with the increment of u.Here the following typical transcendental item in one-dimensional (1-dim) delayed-system is treated as an example for further discussion to Figure4:
(15)
e
-
τ
s
=
a
0
(
s
-
r
)
,
where a
0
,
r
∈
R, and τ
≥
0. The value of s can be derived from (15) based on Lambert W function as follows:
(16)
s
=
r
+
1
τ
W
(
τ
e
-
τ
r
a
0
)
=
r
+
1
τ
W
(
z
)
,
where
(17)
z
=
f
z
(
τ
)
=
τ
e
-
τ
r
a
0
.Figure5 depicts the curve of τ
-
z with a
0
=
-
0.5, r
=
1.0. It is evident that point N with τ
=
0.23196 in the figure is an OEB bifurcation of system (15) according to Definition 1 and Figure 4.Figure 5
Curve ofz
=
f
z
(
τ
).From Figures4 and 5, it is asserted that the following formula holds at the OEB bifurcation points when f
z
(
·
)
∈
R:
(18)
z
=
-
1
e
,
d
z
d
τ
=
d
f
z
d
τ
<
0
.
And, if f
z
(
·
) has complex values on one side of the OEB bifurcation point, the above condition will be modified as follows:
(19)
z
=
-
1
e
,
d
Re
(
z
)
d
τ
=
d
Re
(
f
z
)
d
τ
<
0
,
where Re
(
z
) is the real part of z. For a time-delayed system, if its characteristic equation can be transformed into a suitable form, conditions of (18) or (19) may be adopted as a criterion to determine the OEB bifurcation.


## 4. An Approach to Determine OEB Bifurcation

In this section, we aim to derive a general method to transform the characteristic equation of time-delayed system into a form that can be handled by LambertW function so as to determine the occurrence of OEB bifurcation. Methods suitable for one-dimensional (1-dim) and two-dimensional (2-dim) time-delayed LTI systems are first derived and then they are expanded to the high-dimensional systems.


### 4.1. 1-Dim System

The following 1-dim system is considered in this subsection:(20)
x
˙
=
a
x
+
b
x
(
t
-
τ
)
,
where τ
≥
0, a
,
b
,
x
∈
R, and b
≠
0. The characteristic equation of (20) is
(21)
CE
1
(
s
,
D
)
=
s
-
a
-
b
·
D
=
0
,
where D
=
e
-
τ
s. It can be rewritten as
(22)
D
=
e
-
τ
s
=
1
b
(
s
-
a
)
.
According to (15), s can be solved from (22) based on Lambert W function:
(23)
s
=
a
+
W
(
z
)
τ
,
where
(24)
z
=
f
z
(
τ
)
=
b
·
τ
·
e
-
τ
a
.Similar to Figure5, the OEB bifurcation can be determined from (24) via condition (18).


### 4.2. 2-Dim System

Let us consider the following 2-dim system:(25)
[
x
˙
1
x
˙
2
]
=
[
a
11
a
12
a
21
a
22
]
[
x
1
x
2
]
+
[
b
11
b
12
b
21
b
22
]
[
x
1
(
t
-
τ
)
x
2
(
t
-
τ
)
]
,
where τ
≥
0, a
i
j
,
b
i
j
∈
R, i
,
j
=
1,2; x
=
[
x
1
,
x
2
]
T
∈
R
2. Assume s
1
<
s
2
<
0 are two eigenvalues of the following matrix:
(26)
A
~
=
[
a
11
+
b
11
a
12
+
b
12
a
21
+
b
21
a
22
+
b
22
]
.It means that there are two negative and real eigenvalues atτ
=
0, so system given by (25) is stable without time delay. Supposed that s is an eigenvalue of (25) at τ
>
0 and the following equation holds:
(27)
det
⁡
{
[
s
-
a
11
-
a
12
-
a
21
s
-
a
22
]
-
[
b
11
b
12
b
21
b
22
]
e
-
τ
s
}
=
0
.
System characteristic equation can be derived as follows:
(28)
CE
2
(
s
,
D
)
=
s
2
+
c
1
s
+
c
2
+
c
3
D
s
+
c
4
D
+
c
5
D
2
=
0
.
With c
1
=
-
(
a
11
+
a
22
), c
2
=
a
11
a
22
-
a
12
a
21, c
3
=
-
(
b
11
+
b
22
), c
4
=
a
11
b
22
+
a
22
b
11
-
(
a
12
b
21
+
a
21
b
12
), c
5
=
b
11
b
22
-
b
12
b
21, and D
=
e
-
τ
s.In order to determine the OEB bifurcation in the system, we reorganize (28) into the following form:
(29)
(
s
+
ζ
+
w
·
D
)
(
s
+
z
+
y
·
D
)
-
v
=
0
,
where ζ
,
y
,
w
,
v
,
z
∈
R. And, it can be expanded as follows:
(30)
s
2
+
(
ζ
+
z
)
s
+
(
ζ
z
-
v
)
+
(
w
+
y
)
D
s
+
(
w
z
+
ζ
y
)
D
+
w
y
D
2
=
0
.Comparing (30) with (28), the following equations hold:
(31)
ζ
+
z
=
c
1
,
ζ
z
-
v
=
c
2
,
w
+
y
=
c
3
,
w
z
+
ζ
y
=
c
4
,
w
y
=
c
5
.
And, the following result can be obtained:
(32)
w
=
c
3
±
c
3
2
-
4
c
5
2
,
y
=
c
3
-
w
,
z
=
y
c
1
-
c
4
y
-
w
,
ζ
=
c
1
-
z
,
v
=
ζ
z
-
c
2
when c
3
2
-
4
c
5
<
0 and there is no real solution for (31). That is (28) cannot be converted into the form of (29) in real domain. In order to use the criterion given by (18) or (19), the following two scenarios are considered, respectively.(i) Consider
c
3
2
-
4
c
5
≥
0. Equation (28) can be converted into the form of (29) in real domain. Equation (29) can be rewritten as follows in order to use the condition given by (18) or (19):
(33)
(
s
+
ζ
+
w
·
D
)
(
s
+
z
+
y
·
D
)
=
v
,
(34)
⟹
s
+
ζ
+
w
D
v
=
1
s
+
z
+
y
D
:
=
K
,
(35)
⟹
{
s
+
ζ
+
w
D
=
K
v
s
+
z
+
y
D
=
1
K
,
where K
∈
C is a temporary variable.Via the first formula in (35), the following equation holds:
(36)
D
=
e
-
τ
s
=
-
1
w
(
s
+
ζ
-
K
v
)
.
According to (15) and (16), with a
0,1
=
-
1
/
w, r
1
=
K
v
-
ζ, c
1
=
τ, s has the following form:
(37)
s
=
S
1
(
s
)
=
r
1
+
W
(
z
1
)
c
1
,
where
(38)
z
1
=
f
1
(
τ
)
=
c
1
e
-
c
1
r
1
a
0,1
.Via the second formula in (35), the following equation holds:
(39)
D
=
e
-
τ
s
=
-
1
y
(
s
+
z
-
1
K
)
.Similarly, having defineda
0,2
=
-
1
/
y, r
2
=
1
/
K
-
z, c
2
=
τ, the following result can be obtained:
(40)
s
=
S
2
(
s
)
=
r
2
+
W
(
z
2
)
c
2
,
where
(41)
z
2
=
f
2
(
τ
)
=
c
2
e
-
c
2
r
2
a
0,2
.Equations (37) and (40) can be utilized to solve the values of s and K when τ is given. They can be further taken into (38) and (41) to get the values of z
1 and z
2. And if any one of the following equations holds, there is an OEB bifurcation:
(42)
z
1
=
-
1
e
,
d
Re
(
z
1
)
d
τ
=
d
Re
(
f
1
)
d
τ
<
0
,
(43)
z
2
=
-
1
e
,
d
Re
(
z
2
)
d
τ
=
d
Re
(
f
2
)
d
τ
<
0
.
Obviously, if s is a complex eigenvalue, K defined by (34) and z
1, z
2 given by (38) and (41) are all complex values. They are all real values if s is a real eigenvalue. Hence, condition (19) is suitable for this case.(ii) Consider
c
3
2
-
4
c
5
<
0. Since (28) cannot be converted into the form of (29) in real domain at this condition, it is needed to find a new method for converting (28) into the following form:
(44)
(
s
+
ζ
+
w
·
D
)
(
s
+
z
+
y
·
D
)
-
q
·
D
·
s
-
v
=
0
,
where ζ
,
y
,
w
,
v
,
z
,
q
∈
R. Also, its expansion is given as follows:
(45)
s
2
+
(
ζ
+
z
)
s
+
(
ζ
z
-
v
)
+
(
w
+
y
-
q
)
D
s
+
(
w
z
+
ζ
y
)
D
+
w
y
D
2
=
0
.
Comparing it with (28), the following equations hold:
(46)
ζ
+
z
=
c
1
,
ζ
z
-
v
=
c
2
,
w
+
y
=
c
3
+
q
,
w
z
+
ζ
y
=
c
4
,
w
y
=
c
5
.Hence,(47)
w
=
c
3
+
q
±
(
c
3
+
q
)
2
-
4
c
5
2
,
(48)
y
=
c
3
+
q
-
w
=
c
3
+
q
∓
(
c
3
+
q
)
2
-
4
c
5
2
,
(49)
z
=
y
c
1
-
c
4
y
-
w
,
(50)
ζ
=
c
1
-
z
,
(51)
v
=
ζ
z
-
c
2
.
To guarantee that ζ, y, w, v, z are real numbers, the following relationship must hold:
(52)
(
c
3
+
q
)
2
-
4
c
5
≥
0
.Sinceq can be selected arbitrarily, it is set as follows:
(53)
q
=
-
c
3
+
5
c
5
.Substituting (53) into (51) yields
(54)
w
=
c
3
+
q
±
(
c
3
+
q
)
2
-
4
c
5
2
=
(
5
±
1
)
c
5
2
,
y
=
c
3
+
q
-
w
=
(
5
∓
1
)
c
5
2
.Sincec
3
2
-
4
c
5
<
0, it must be true that c
5
>
0. So w and y are two real numbers. According to (49)–(51), ζ, v, z are all real numbers. In order to use condition (19), (44) is written as
(55)
(
s
+
ζ
+
w
D
)
(
s
+
z
+
y
D
)
=
v
+
q
D
s
,
(56)
⟹
s
+
ζ
+
w
D
q
D
s
+
v
=
1
s
+
z
+
y
D
:
=
K
,
(57)
⟹
{
s
+
ζ
+
w
D
=
K
(
v
+
q
D
s
)
s
+
z
+
y
D
=
1
K
.
0000000It can be indicated that the second formula in (57) is equivalent to the one in (35). Thus, it can generate the same results given by (40), (41), and condition (43). The new conditions will be derived from the first formula in (57) as follows:
(58)
s
+
ζ
+
w
D
=
K
(
v
+
q
D
s
)
⟹
D
(
K
q
s
-
w
)
=
s
+
ζ
-
K
v
(59)
⟹
D
=
e
-
τ
s
=
s
+
ζ
-
K
v
K
q
s
-
w
∶
=
P
n
(
s
)
P
d
(
s
)
.Suppose that there exists a real variableQ
∈
R such that
(60)
e
-
Q
τ
s
=
sign
⁡
(
P
n
(
s
)
)
·
(
s
+
ζ
-
K
v
)
.Inserting it into (59) yields
(61)
e
-
(
1
-
Q
)
τ
s
=
sign
⁡
(
P
d
(
s
)
)
K
q
s
-
w
.Defininga
0,3
=
sign
⁡
(
P
n
(
s
)
), r
3
=
K
v
-
ζ, c
3
=
Q
τ, eigenvalue s is expressed as
(62)
s
=
S
3
(
s
)
=
r
3
+
W
(
z
3
)
c
3
,
where
(63)
z
3
=
f
3
(
τ
)
=
c
3
e
-
c
3
r
3
a
0,3
.Defininga
0,4
=
sign
⁡
(
P
d
(
s
)
)
·
K
q, r
4
=
w
/
K
q, c
4
=
(
Q
-
1
)
τ, the similar result from (61) is obtained:
(64)
s
=
S
4
(
s
)
=
r
4
+
W
(
z
4
)
c
4
,
where
(65)
z
4
=
f
4
(
τ
)
=
c
4
e
-
c
4
r
4
a
0,4
.Equations (40), (62), and (64) can be used to solve the values of s, K, and Q when τ is given. They can be further taken into (41), (63), and (65) to get the values of z
2, z
3, and z
4. As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system:
(66)
z
2
=
-
1
e
,
d
Re
(
z
2
)
d
τ
=
d
Re
(
f
2
)
d
τ
<
0
,
z
3
=
-
1
e
,
d
Re
(
z
3
)
d
τ
=
d
Re
(
f
3
)
d
τ
<
0
,
z
4
=
-
1
e
,
d
Re
(
z
4
)
d
τ
=
d
Re
(
f
4
)
d
τ
<
0
.


### 4.3. High-Dimensional Systems

Consider systemΣ given by (1) with n
≥
3. Its characteristic equation has the following general form:
(67)
CE
(
s
,
D
)
=
∑
i
=
0
n
[
p
n
-
i
(
s
)
D
i
]
=
0
,
where D
=
e
-
τ
s and p
n
-
i
(
s
) is a polynomial of s:
(68)
p
n
-
i
(
s
)
=
∑
j
=
0
n
-
i
a
j
s
j
.
Some general forms of the characteristic equation for system Σ are listed in Table 1. N
c is the number of coefficients in the general form. For example, there are five coefficients which are a
1
~
a
5 in the characteristic equation for n
=
2, and there are nine coefficients a
1
~
a
9 in the characteristic equation with n
=
3.Table 1
General forms of the characteristic equation for time-delayed system.
Dim
General forms of the characteristic equation
N
c
1
s
+
a
1
+
a
2
D
=
0
2
2
s
2
+
a
1
s
+
a
2
+
(
a
3
s
+
a
4
)
D
+
a
5
D
2
=
0
5
3
s
3
+
a
1
s
2
+
a
2
s
+
a
3
+
(
a
4
s
2
+
a
5
s
+
a
6
)
D
+
(
a
7
s
+
a
8
)
D
2
+
a
9
D
3
=
0
9
4
s
4
+
a
1
s
3
+
a
2
s
2
+
a
3
s
+
a
4
+
(
a
5
s
3
+
a
6
s
2
+
a
7
s
+
a
8
)
D
+
(
a
9
s
2
+
a
10
s
+
a
11
)
D
2
+
(
a
12
s
+
a
13
)
D
3
+
a
14
D
3
=
0
14
5
s
5
+
a
1
s
4
+
a
2
s
3
+
a
3
s
2
+
a
4
s
+
a
5
+
(
a
6
s
4
+
a
7
s
3
+
a
8
s
2
+
a
9
s
+
a
10
)
D
+
(
a
11
s
3
+
a
12
s
2
+
a
13
s
+
a
14
)
D
2
+
(
a
15
s
2
+
c
16
s
+
a
17
)
D
3
+
(
a
18
s
+
a
19
)
D
4
+
a
20
D
5
=
0
20(i) 3-Dim Case (
n
=
3
). Take 3-dim system as an example to explain the principle of proposed method. General form of the characteristic equation for n
=
3 is given as below:
(69)
CE
3
(
s
,
D
)
=
s
3
+
a
1
s
2
+
a
2
s
+
a
3
+
(
a
4
s
2
+
a
5
s
+
a
6
)
D
+
(
a
7
s
+
a
8
)
D
2
+
a
9
D
3
=
0
.It can be converted into the following form so as to utilize the previous results:(70)
CE
3
(
s
,
D
)
=
CE
2
(
s
,
D
)
·
(
s
+
c
6
+
c
7
D
)
-
(
c
8
s
+
c
9
)
=
0
,
where CE
2
(
s
,
D
) is described by (28). And, c
i
∈
R, i
=
1,2
,
…
,
9, are variables to be determined. Setting CE
2
(
s
,
D
)
=
K, (70) can be written as
(71)
CE
2
(
s
,
D
)
-
K
=
0
,
K
·
(
s
+
c
6
+
c
7
D
)
-
(
c
8
s
+
c
9
)
=
0
.
Via the first formula in (71),
(72)
CE
2
′
(
s
,
D
)
≔
CE
2
(
s
,
D
)
-
K
=
s
2
+
c
1
s
+
(
c
2
-
K
)
+
c
3
D
s
+
c
4
D
+
c
5
D
2
=
0
.
It is obvious that CE
2
′
(
s
,
D
) is a typical 2-dim case discussed above. Hence, the result of Section 4.2 can be directly applied to (72).Further, via the second formula of (71),
(73)
K
·
(
s
+
c
6
+
c
7
D
)
-
(
c
8
s
+
c
9
)
=
0
⟹
s
-
a
~
-
b
~
·
D
=
0
,
where a
~
=
(
c
9
-
K
c
6
)
/
(
K
-
c
8
), b
~
=
-
K
c
7
/
(
K
-
c
8
). The result of 1-dim system in Section 4.1 can be applied to (73). And the following result can also be derived:
(74)
z
5
=
f
5
(
τ
)
=
b
~
·
τ
·
e
-
τ
a
~
.From the above discussion, it can be asserted that a 3-dim system can be transformed to a combination of a 1-dim system and a 2-dim system. Therefore, conclusions in Sections4.1 and 4.2 can be applied. For simplification, it is denoted as
(75)
CE
3
(
s
,
D
)
=
CE
2
(
s
,
D
)
⊕
CE
1
(
s
,
D
)
.As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system:(76)
z
1
=
-
1
e
,
d
Re
(
z
1
)
d
τ
=
d
Re
(
f
1
)
d
τ
<
0
,
(77)
z
2
=
-
1
e
,
d
Re
(
z
2
)
d
τ
=
d
Re
(
f
2
)
d
τ
<
0
,
(78)
z
5
=
-
1
e
,
d
Re
(
z
5
)
d
τ
=
d
Re
(
f
5
)
d
τ
<
0
,
where (76) and (77) are concluded from (72) with conditions expressed by (42)-(43) when c
3
2
-
4
c
5
≥
0 (or conditions expressed by (66) when c
3
2
-
4
c
5
<
0).(ii) Case of High Dimension with
n
>
3. Similar to 3-dim case, the characteristic equation of high-dimensional system can also be converted into lower ones so as to apply the previous derived results. For example, for n
=
4,
(79)
CE
4
(
s
,
D
)
=
CE
3
(
s
,
D
)
·
(
s
+
c
10
+
c
11
D
)
-
(
c
12
s
2
+
c
13
s
+
c
14
)
=
0
.SettingCE
3
(
s
,
D
)
=
K, (79) can be rewritten as
(80)
CE
3
(
s
,
D
)
-
K
=
0
,
CE
4
(
s
,
D
)
=
K
·
(
s
+
c
10
+
c
11
D
)
-
(
c
12
s
2
+
c
13
s
+
c
14
)
=
0
.The first formula is a 3-dim case, and the second is a 2-dim case. Hence,(81)
CE
4
(
s
,
D
)
=
CE
3
(
s
,
D
)
⊕
CE
2
(
s
,
D
)
=
[
CE
2
(
s
,
D
)
⊕
CE
1
(
s
,
D
)
]
⊕
CE
2
(
s
,
D
)
.
And, the cases discussed for 2-dim and 1-dim cases can be applied to determine OEB bifurcation in 4-dim system.Further, the following recursion formula for high-dimensional systems withn
≥
3 can be obtained:
(82)
CE
i
(
s
,
D
)
=
CE
i
-
1
(
s
,
D
)
⊕
CE
i
-
2
(
s
,
D
)
,
i
=
3,4
,
…
.Then, a higher dimension system can be converted into two lower dimension systems so as to determine the occurrence of OEB bifurcation accordingly.


## 4.1. 1-Dim System

The following 1-dim system is considered in this subsection:(20)
x
˙
=
a
x
+
b
x
(
t
-
τ
)
,
where τ
≥
0, a
,
b
,
x
∈
R, and b
≠
0. The characteristic equation of (20) is
(21)
CE
1
(
s
,
D
)
=
s
-
a
-
b
·
D
=
0
,
where D
=
e
-
τ
s. It can be rewritten as
(22)
D
=
e
-
τ
s
=
1
b
(
s
-
a
)
.
According to (15), s can be solved from (22) based on Lambert W function:
(23)
s
=
a
+
W
(
z
)
τ
,
where
(24)
z
=
f
z
(
τ
)
=
b
·
τ
·
e
-
τ
a
.Similar to Figure5, the OEB bifurcation can be determined from (24) via condition (18).


## 4.2. 2-Dim System

Let us consider the following 2-dim system:(25)
[
x
˙
1
x
˙
2
]
=
[
a
11
a
12
a
21
a
22
]
[
x
1
x
2
]
+
[
b
11
b
12
b
21
b
22
]
[
x
1
(
t
-
τ
)
x
2
(
t
-
τ
)
]
,
where τ
≥
0, a
i
j
,
b
i
j
∈
R, i
,
j
=
1,2; x
=
[
x
1
,
x
2
]
T
∈
R
2. Assume s
1
<
s
2
<
0 are two eigenvalues of the following matrix:
(26)
A
~
=
[
a
11
+
b
11
a
12
+
b
12
a
21
+
b
21
a
22
+
b
22
]
.It means that there are two negative and real eigenvalues atτ
=
0, so system given by (25) is stable without time delay. Supposed that s is an eigenvalue of (25) at τ
>
0 and the following equation holds:
(27)
det
⁡
{
[
s
-
a
11
-
a
12
-
a
21
s
-
a
22
]
-
[
b
11
b
12
b
21
b
22
]
e
-
τ
s
}
=
0
.
System characteristic equation can be derived as follows:
(28)
CE
2
(
s
,
D
)
=
s
2
+
c
1
s
+
c
2
+
c
3
D
s
+
c
4
D
+
c
5
D
2
=
0
.
With c
1
=
-
(
a
11
+
a
22
), c
2
=
a
11
a
22
-
a
12
a
21, c
3
=
-
(
b
11
+
b
22
), c
4
=
a
11
b
22
+
a
22
b
11
-
(
a
12
b
21
+
a
21
b
12
), c
5
=
b
11
b
22
-
b
12
b
21, and D
=
e
-
τ
s.In order to determine the OEB bifurcation in the system, we reorganize (28) into the following form:
(29)
(
s
+
ζ
+
w
·
D
)
(
s
+
z
+
y
·
D
)
-
v
=
0
,
where ζ
,
y
,
w
,
v
,
z
∈
R. And, it can be expanded as follows:
(30)
s
2
+
(
ζ
+
z
)
s
+
(
ζ
z
-
v
)
+
(
w
+
y
)
D
s
+
(
w
z
+
ζ
y
)
D
+
w
y
D
2
=
0
.Comparing (30) with (28), the following equations hold:
(31)
ζ
+
z
=
c
1
,
ζ
z
-
v
=
c
2
,
w
+
y
=
c
3
,
w
z
+
ζ
y
=
c
4
,
w
y
=
c
5
.
And, the following result can be obtained:
(32)
w
=
c
3
±
c
3
2
-
4
c
5
2
,
y
=
c
3
-
w
,
z
=
y
c
1
-
c
4
y
-
w
,
ζ
=
c
1
-
z
,
v
=
ζ
z
-
c
2
when c
3
2
-
4
c
5
<
0 and there is no real solution for (31). That is (28) cannot be converted into the form of (29) in real domain. In order to use the criterion given by (18) or (19), the following two scenarios are considered, respectively.(i) Consider
c
3
2
-
4
c
5
≥
0. Equation (28) can be converted into the form of (29) in real domain. Equation (29) can be rewritten as follows in order to use the condition given by (18) or (19):
(33)
(
s
+
ζ
+
w
·
D
)
(
s
+
z
+
y
·
D
)
=
v
,
(34)
⟹
s
+
ζ
+
w
D
v
=
1
s
+
z
+
y
D
:
=
K
,
(35)
⟹
{
s
+
ζ
+
w
D
=
K
v
s
+
z
+
y
D
=
1
K
,
where K
∈
C is a temporary variable.Via the first formula in (35), the following equation holds:
(36)
D
=
e
-
τ
s
=
-
1
w
(
s
+
ζ
-
K
v
)
.
According to (15) and (16), with a
0,1
=
-
1
/
w, r
1
=
K
v
-
ζ, c
1
=
τ, s has the following form:
(37)
s
=
S
1
(
s
)
=
r
1
+
W
(
z
1
)
c
1
,
where
(38)
z
1
=
f
1
(
τ
)
=
c
1
e
-
c
1
r
1
a
0,1
.Via the second formula in (35), the following equation holds:
(39)
D
=
e
-
τ
s
=
-
1
y
(
s
+
z
-
1
K
)
.Similarly, having defineda
0,2
=
-
1
/
y, r
2
=
1
/
K
-
z, c
2
=
τ, the following result can be obtained:
(40)
s
=
S
2
(
s
)
=
r
2
+
W
(
z
2
)
c
2
,
where
(41)
z
2
=
f
2
(
τ
)
=
c
2
e
-
c
2
r
2
a
0,2
.Equations (37) and (40) can be utilized to solve the values of s and K when τ is given. They can be further taken into (38) and (41) to get the values of z
1 and z
2. And if any one of the following equations holds, there is an OEB bifurcation:
(42)
z
1
=
-
1
e
,
d
Re
(
z
1
)
d
τ
=
d
Re
(
f
1
)
d
τ
<
0
,
(43)
z
2
=
-
1
e
,
d
Re
(
z
2
)
d
τ
=
d
Re
(
f
2
)
d
τ
<
0
.
Obviously, if s is a complex eigenvalue, K defined by (34) and z
1, z
2 given by (38) and (41) are all complex values. They are all real values if s is a real eigenvalue. Hence, condition (19) is suitable for this case.(ii) Consider
c
3
2
-
4
c
5
<
0. Since (28) cannot be converted into the form of (29) in real domain at this condition, it is needed to find a new method for converting (28) into the following form:
(44)
(
s
+
ζ
+
w
·
D
)
(
s
+
z
+
y
·
D
)
-
q
·
D
·
s
-
v
=
0
,
where ζ
,
y
,
w
,
v
,
z
,
q
∈
R. Also, its expansion is given as follows:
(45)
s
2
+
(
ζ
+
z
)
s
+
(
ζ
z
-
v
)
+
(
w
+
y
-
q
)
D
s
+
(
w
z
+
ζ
y
)
D
+
w
y
D
2
=
0
.
Comparing it with (28), the following equations hold:
(46)
ζ
+
z
=
c
1
,
ζ
z
-
v
=
c
2
,
w
+
y
=
c
3
+
q
,
w
z
+
ζ
y
=
c
4
,
w
y
=
c
5
.Hence,(47)
w
=
c
3
+
q
±
(
c
3
+
q
)
2
-
4
c
5
2
,
(48)
y
=
c
3
+
q
-
w
=
c
3
+
q
∓
(
c
3
+
q
)
2
-
4
c
5
2
,
(49)
z
=
y
c
1
-
c
4
y
-
w
,
(50)
ζ
=
c
1
-
z
,
(51)
v
=
ζ
z
-
c
2
.
To guarantee that ζ, y, w, v, z are real numbers, the following relationship must hold:
(52)
(
c
3
+
q
)
2
-
4
c
5
≥
0
.Sinceq can be selected arbitrarily, it is set as follows:
(53)
q
=
-
c
3
+
5
c
5
.Substituting (53) into (51) yields
(54)
w
=
c
3
+
q
±
(
c
3
+
q
)
2
-
4
c
5
2
=
(
5
±
1
)
c
5
2
,
y
=
c
3
+
q
-
w
=
(
5
∓
1
)
c
5
2
.Sincec
3
2
-
4
c
5
<
0, it must be true that c
5
>
0. So w and y are two real numbers. According to (49)–(51), ζ, v, z are all real numbers. In order to use condition (19), (44) is written as
(55)
(
s
+
ζ
+
w
D
)
(
s
+
z
+
y
D
)
=
v
+
q
D
s
,
(56)
⟹
s
+
ζ
+
w
D
q
D
s
+
v
=
1
s
+
z
+
y
D
:
=
K
,
(57)
⟹
{
s
+
ζ
+
w
D
=
K
(
v
+
q
D
s
)
s
+
z
+
y
D
=
1
K
.
0000000It can be indicated that the second formula in (57) is equivalent to the one in (35). Thus, it can generate the same results given by (40), (41), and condition (43). The new conditions will be derived from the first formula in (57) as follows:
(58)
s
+
ζ
+
w
D
=
K
(
v
+
q
D
s
)
⟹
D
(
K
q
s
-
w
)
=
s
+
ζ
-
K
v
(59)
⟹
D
=
e
-
τ
s
=
s
+
ζ
-
K
v
K
q
s
-
w
∶
=
P
n
(
s
)
P
d
(
s
)
.Suppose that there exists a real variableQ
∈
R such that
(60)
e
-
Q
τ
s
=
sign
⁡
(
P
n
(
s
)
)
·
(
s
+
ζ
-
K
v
)
.Inserting it into (59) yields
(61)
e
-
(
1
-
Q
)
τ
s
=
sign
⁡
(
P
d
(
s
)
)
K
q
s
-
w
.Defininga
0,3
=
sign
⁡
(
P
n
(
s
)
), r
3
=
K
v
-
ζ, c
3
=
Q
τ, eigenvalue s is expressed as
(62)
s
=
S
3
(
s
)
=
r
3
+
W
(
z
3
)
c
3
,
where
(63)
z
3
=
f
3
(
τ
)
=
c
3
e
-
c
3
r
3
a
0,3
.Defininga
0,4
=
sign
⁡
(
P
d
(
s
)
)
·
K
q, r
4
=
w
/
K
q, c
4
=
(
Q
-
1
)
τ, the similar result from (61) is obtained:
(64)
s
=
S
4
(
s
)
=
r
4
+
W
(
z
4
)
c
4
,
where
(65)
z
4
=
f
4
(
τ
)
=
c
4
e
-
c
4
r
4
a
0,4
.Equations (40), (62), and (64) can be used to solve the values of s, K, and Q when τ is given. They can be further taken into (41), (63), and (65) to get the values of z
2, z
3, and z
4. As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system:
(66)
z
2
=
-
1
e
,
d
Re
(
z
2
)
d
τ
=
d
Re
(
f
2
)
d
τ
<
0
,
z
3
=
-
1
e
,
d
Re
(
z
3
)
d
τ
=
d
Re
(
f
3
)
d
τ
<
0
,
z
4
=
-
1
e
,
d
Re
(
z
4
)
d
τ
=
d
Re
(
f
4
)
d
τ
<
0
.


## 4.3. High-Dimensional Systems

Consider systemΣ given by (1) with n
≥
3. Its characteristic equation has the following general form:
(67)
CE
(
s
,
D
)
=
∑
i
=
0
n
[
p
n
-
i
(
s
)
D
i
]
=
0
,
where D
=
e
-
τ
s and p
n
-
i
(
s
) is a polynomial of s:
(68)
p
n
-
i
(
s
)
=
∑
j
=
0
n
-
i
a
j
s
j
.
Some general forms of the characteristic equation for system Σ are listed in Table 1. N
c is the number of coefficients in the general form. For example, there are five coefficients which are a
1
~
a
5 in the characteristic equation for n
=
2, and there are nine coefficients a
1
~
a
9 in the characteristic equation with n
=
3.Table 1
General forms of the characteristic equation for time-delayed system.
Dim
General forms of the characteristic equation
N
c
1
s
+
a
1
+
a
2
D
=
0
2
2
s
2
+
a
1
s
+
a
2
+
(
a
3
s
+
a
4
)
D
+
a
5
D
2
=
0
5
3
s
3
+
a
1
s
2
+
a
2
s
+
a
3
+
(
a
4
s
2
+
a
5
s
+
a
6
)
D
+
(
a
7
s
+
a
8
)
D
2
+
a
9
D
3
=
0
9
4
s
4
+
a
1
s
3
+
a
2
s
2
+
a
3
s
+
a
4
+
(
a
5
s
3
+
a
6
s
2
+
a
7
s
+
a
8
)
D
+
(
a
9
s
2
+
a
10
s
+
a
11
)
D
2
+
(
a
12
s
+
a
13
)
D
3
+
a
14
D
3
=
0
14
5
s
5
+
a
1
s
4
+
a
2
s
3
+
a
3
s
2
+
a
4
s
+
a
5
+
(
a
6
s
4
+
a
7
s
3
+
a
8
s
2
+
a
9
s
+
a
10
)
D
+
(
a
11
s
3
+
a
12
s
2
+
a
13
s
+
a
14
)
D
2
+
(
a
15
s
2
+
c
16
s
+
a
17
)
D
3
+
(
a
18
s
+
a
19
)
D
4
+
a
20
D
5
=
0
20(i) 3-Dim Case (
n
=
3
). Take 3-dim system as an example to explain the principle of proposed method. General form of the characteristic equation for n
=
3 is given as below:
(69)
CE
3
(
s
,
D
)
=
s
3
+
a
1
s
2
+
a
2
s
+
a
3
+
(
a
4
s
2
+
a
5
s
+
a
6
)
D
+
(
a
7
s
+
a
8
)
D
2
+
a
9
D
3
=
0
.It can be converted into the following form so as to utilize the previous results:(70)
CE
3
(
s
,
D
)
=
CE
2
(
s
,
D
)
·
(
s
+
c
6
+
c
7
D
)
-
(
c
8
s
+
c
9
)
=
0
,
where CE
2
(
s
,
D
) is described by (28). And, c
i
∈
R, i
=
1,2
,
…
,
9, are variables to be determined. Setting CE
2
(
s
,
D
)
=
K, (70) can be written as
(71)
CE
2
(
s
,
D
)
-
K
=
0
,
K
·
(
s
+
c
6
+
c
7
D
)
-
(
c
8
s
+
c
9
)
=
0
.
Via the first formula in (71),
(72)
CE
2
′
(
s
,
D
)
≔
CE
2
(
s
,
D
)
-
K
=
s
2
+
c
1
s
+
(
c
2
-
K
)
+
c
3
D
s
+
c
4
D
+
c
5
D
2
=
0
.
It is obvious that CE
2
′
(
s
,
D
) is a typical 2-dim case discussed above. Hence, the result of Section 4.2 can be directly applied to (72).Further, via the second formula of (71),
(73)
K
·
(
s
+
c
6
+
c
7
D
)
-
(
c
8
s
+
c
9
)
=
0
⟹
s
-
a
~
-
b
~
·
D
=
0
,
where a
~
=
(
c
9
-
K
c
6
)
/
(
K
-
c
8
), b
~
=
-
K
c
7
/
(
K
-
c
8
). The result of 1-dim system in Section 4.1 can be applied to (73). And the following result can also be derived:
(74)
z
5
=
f
5
(
τ
)
=
b
~
·
τ
·
e
-
τ
a
~
.From the above discussion, it can be asserted that a 3-dim system can be transformed to a combination of a 1-dim system and a 2-dim system. Therefore, conclusions in Sections4.1 and 4.2 can be applied. For simplification, it is denoted as
(75)
CE
3
(
s
,
D
)
=
CE
2
(
s
,
D
)
⊕
CE
1
(
s
,
D
)
.As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system:(76)
z
1
=
-
1
e
,
d
Re
(
z
1
)
d
τ
=
d
Re
(
f
1
)
d
τ
<
0
,
(77)
z
2
=
-
1
e
,
d
Re
(
z
2
)
d
τ
=
d
Re
(
f
2
)
d
τ
<
0
,
(78)
z
5
=
-
1
e
,
d
Re
(
z
5
)
d
τ
=
d
Re
(
f
5
)
d
τ
<
0
,
where (76) and (77) are concluded from (72) with conditions expressed by (42)-(43) when c
3
2
-
4
c
5
≥
0 (or conditions expressed by (66) when c
3
2
-
4
c
5
<
0).(ii) Case of High Dimension with
n
>
3. Similar to 3-dim case, the characteristic equation of high-dimensional system can also be converted into lower ones so as to apply the previous derived results. For example, for n
=
4,
(79)
CE
4
(
s
,
D
)
=
CE
3
(
s
,
D
)
·
(
s
+
c
10
+
c
11
D
)
-
(
c
12
s
2
+
c
13
s
+
c
14
)
=
0
.SettingCE
3
(
s
,
D
)
=
K, (79) can be rewritten as
(80)
CE
3
(
s
,
D
)
-
K
=
0
,
CE
4
(
s
,
D
)
=
K
·
(
s
+
c
10
+
c
11
D
)
-
(
c
12
s
2
+
c
13
s
+
c
14
)
=
0
.The first formula is a 3-dim case, and the second is a 2-dim case. Hence,(81)
CE
4
(
s
,
D
)
=
CE
3
(
s
,
D
)
⊕
CE
2
(
s
,
D
)
=
[
CE
2
(
s
,
D
)
⊕
CE
1
(
s
,
D
)
]
⊕
CE
2
(
s
,
D
)
.
And, the cases discussed for 2-dim and 1-dim cases can be applied to determine OEB bifurcation in 4-dim system.Further, the following recursion formula for high-dimensional systems withn
≥
3 can be obtained:
(82)
CE
i
(
s
,
D
)
=
CE
i
-
1
(
s
,
D
)
⊕
CE
i
-
2
(
s
,
D
)
,
i
=
3,4
,
…
.Then, a higher dimension system can be converted into two lower dimension systems so as to determine the occurrence of OEB bifurcation accordingly.


## 5. Case Studies

### 5.1. 1-Dim System

Consider the following 1-dim system:(83)
x
˙
=
20
x
-
30
x
(
t
-
τ
)
.According to (24), the τ
-
z curve can be depicted in Figure 6 and the corresponding eigenvalue locus is depicted in Figure 7. For a pair of conjugated eigenvalues, only the one above the real axis is showed in the figures in this paper.Figure 6
Curve ofτ
-
z.Figure 7
Curve of an eigenvalue.Similar to Figure5, point N with (
τ
,
z
)
=
(
0.0173
s
,
-
1
/
e
) in Figure 6 is an OEB bifurcation of system described by (83). The variation of the system eigenvalue in Figure 7 verifies the analysis above. That is, the real eigenvalue splits into a pair of conjugated complex eigenvalues at the OEB bifurcation with τ
=
0.0173
s.In order to demonstrate the impact of OEB bifurcation to the system defined by (83), its unit step responses with different τ values are shown in Figure 8. It can be found that the system trajectory converges monotonously in the range of τ
<
0.0173
s and converges oscillatedly in the range of τ
>
0.0173
s. It is manifested that there is OEB bifurcation at τ
=
0.0173
s, which leads to a new oscillation mode and increases the number of system eigenvalues by one.Simulation results with variousτ values.
(a)
Simulation results withτ
=
0.01
s and τ
=
0.0173
s
(b)
Simulation results withτ
=
0.02
s and τ
=
0.03
s


### 5.2. 2-Dim System

The following two cases are employed to validate the proposed method.(i) Case A.Consider
(84)
A
0
=
[
a
11
a
12
a
21
a
22
]
=
[
1
-
3
5
-
6
]
,
A
i
=
[
b
11
b
12
b
21
b
22
]
=
[
-
4
4
2
-
3
]
,
eig
(
A
-
)
=
[
-
2.0
,
-
10.0
]
,
c
3
2
-
4
c
5
=
33.0
>
0
,
where A
- is defined by (26).It is obvious that the conditions described by (42)-(43) are suitable for this case.Via the tracing method proposed by [23], we can obtain the system eigenvalue locus. Figure 9 depicts the curve of the second eigenvalue with τ increasing (the first one is always on the real axis, so it is not shown in the figure). There is an OEB bifurcation at point Q with τ
=
0.0476
s. According to (38) and (41), the curves of τ
-
z
1 and τ
-
z
2 can be depicted in Figure 10. Solid line represents that the corresponding z
i, i
=
1,2, is real, while dotted line represents that it is complex. Some results nearby point Q are listed in Table 2. It can be found that condition (42) holds at point Q with (
τ
,
z
1
)
R
=
(
0.0476
s
,
-
1
/
e
) and it indicates the occurrence of OEB bifurcation precisely.Table 2
z
1
,
z
2 values near point Q in Figure 10.
τ
Re
(
z
1
)
Re
(
z
1
)
+
1
/
e
Re
(
z
2
)
Re
(
z
2
)
+
1
/
e
0.04758
−0.3675333
0.0003462
−0.0858096
0.2820698
0.04759
−0.3676635
0.0002159
−0.0865306
0.2813488
0.04760
−0.3678304
0.0000491
−0.0879570
0.2799224
0.04761
−0.3679437
−0.0000643
−0.0883511
0.2795284
0.04762
−0.3680362
−0.0001567
−0.0883415
0.2795379Figure 9
Locus of the second eigenvalue withτ increasing.Figure 10
Curves ofτ
-
z
1 and τ
-
z
2.(ii) Case B.Consider
(85)
A
0
=
[
a
11
a
12
a
21
a
22
]
=
[
1
2
5
-
6
]
,
A
i
=
[
b
11
b
12
b
21
b
22
]
=
[
-
4
-
1
2
-
3
]
,
eig
(
A
-
)
=
[
-
2.0
,
-
10.0
]
,
c
3
2
-
4
c
5
=
-
7.0
<
0
,
where A
- is defined by (26).It can be asserted that the conditions described by (66) are appropriate for this case. Similar to Case A, the eigenvalue locus with τ increasing can be obtained via the tracing method expressed by [23]. The eigenvalue locus is shown in Figure 11, while real and imaginary parts of the eigenvalues are shown in Figure 12. Note that two real eigenvalues at the beginning move close to each other on the real axis. They encounter at point A and then turn to a pair of conjugated complex eigenvalues. According to Definition 1, there is no OEB bifurcation in the system since the number of eigenvalues does not increase here.Figure 11
Eigenvalue locus of Case B.Real and imaginary parts of two eigenvalues.
(a)
Real parts of two eigenvalues withτ increasing
(b)
Imaginary parts of two eigenvalues withτ increasingFigure13 plots the curves of τ
-
z
2, τ
-
z
3, and τ
-
z
4. Similarly, solid and dotted lines are used to represent the real and complex values of z
i, i
=
2,3
,
4. It is obvious that any condition identified by (66) does not hold here, so there is no OEB bifurcation in the system.Figure 13
Curves ofτ
-
z
2, τ
-
z
3, and τ
-
z
4.Case studies above demonstrate that conditions expressed by (42)-(43) or (66) are correct. And, the method proposed can correctly indicate the occurrence of OEB bifurcation in the 2-dim time-delayed system.


### 5.3. 3-Dim System

Consider a 3-dim system defined by (1) with A
0, A
i as follows:
(86)
A
0
=
[
-
5
1
-
1
1
-
4
2
0
0
-
2
]
,
A
i
=
[
2
1
0
1
-
1
-
1
1
3
-
8
]
.
It can be found that the system have three real eigenvalues with τ
=
0:
(87)
eig
(
A
-
)
=
[
-
1.8692
,
-
5.5749
,
-
10.5559
]
,
where A
-
=
A
0
+
A
i.Similar to Section5.2, we can get the eigenvalue locus shown in Figure 14, where the first and second eigenvalues always move on the real axis, while the third eigenvalue splits into a pair of conjugated complex eigenvalues at point B with τ
=
0.0431
s so as to yield an OEB bifurcation.Figure 14
Eigenvalue locus with time delay increasing.Table3 gives the coefficients in (69) and (70). According to the definition of z
1
,
z
2
,
z
5 by (76)–(78), we can depict the curves of τ
-
z
1, τ
-
z
2, and τ
-
z
5 in Figure 15. The τ
-
z
5 curve (red line) hits the horizontal line of -
1
/
e at τ
=
0.0431
s and condition (78) holds which indicate the occurrence of the OEB bifurcation.Table 3
Coefficients ofC
E
3
(
s
,
D
) given by (69) and (70).
a
1
a
2
a
3
a
4
a
5
a
6
a
7
a
8
a
9
11
37
38
7
60
117
−8
−16
−29
c
1
c
2
c
3
c
4
c
5
c
6
c
7
c
8
c
9
8.2422
15.7565
−0.5508
−0.7162
−3.8406
2.7578
7.5508
1.4870
5.4538Figure 15
Curves ofτ
-
z
1, τ
-
z
2, and τ
-
z
5.It should be mentioned that the transformation from (69) to (70) is not unique. Table 4 shows two sets of coefficients in (70). The corresponding τ
-
z
1, τ
-
z
2, and τ
-
z
5 curves are illustrated, respectively, by Figures 16 and 17. It is notable that in Figure 16 or Figure 17, the τ
-
z
1 curve (black line) collides with the horizontal line of -
1
/
e at τ
=
0.0431
s. The condition of (76) holds at this point so as to indicate occurrence of the OEB bifurcation.Table 4
Two sets of coefficients in (70).
c
1
c
2
c
3
c
4
c
5
c
6
c
7
c
8
c
9
Set I
6.4863
13.8753
9.2544
32.8516
12.8635
4.5137
−2.2544
6.1525
24.6287
Set II
7.2715
7.6935
5.2964
27.8647
−17.0229
3.7285
1.7036
−2.1947
−9.3149Figure 16
Curves ofτ
-
z
1, τ
-
z
2, and τ
-
z
5 with set I of Table 4.Figure 17
Curves ofτ
-
z
1, τ
-
z
2, and τ
-
z
5 with set II of Table 4.From the case study above, it is verified that conditions expressed by (76)–(78) can indicate the occurrence of OEB bifurcation in the 3-dim time-delayed LTI system.


## 5.1. 1-Dim System

Consider the following 1-dim system:(83)
x
˙
=
20
x
-
30
x
(
t
-
τ
)
.According to (24), the τ
-
z curve can be depicted in Figure 6 and the corresponding eigenvalue locus is depicted in Figure 7. For a pair of conjugated eigenvalues, only the one above the real axis is showed in the figures in this paper.Figure 6
Curve ofτ
-
z.Figure 7
Curve of an eigenvalue.Similar to Figure5, point N with (
τ
,
z
)
=
(
0.0173
s
,
-
1
/
e
) in Figure 6 is an OEB bifurcation of system described by (83). The variation of the system eigenvalue in Figure 7 verifies the analysis above. That is, the real eigenvalue splits into a pair of conjugated complex eigenvalues at the OEB bifurcation with τ
=
0.0173
s.In order to demonstrate the impact of OEB bifurcation to the system defined by (83), its unit step responses with different τ values are shown in Figure 8. It can be found that the system trajectory converges monotonously in the range of τ
<
0.0173
s and converges oscillatedly in the range of τ
>
0.0173
s. It is manifested that there is OEB bifurcation at τ
=
0.0173
s, which leads to a new oscillation mode and increases the number of system eigenvalues by one.Simulation results with variousτ values.
(a)
Simulation results withτ
=
0.01
s and τ
=
0.0173
s
(b)
Simulation results withτ
=
0.02
s and τ
=
0.03
s


## 5.2. 2-Dim System

The following two cases are employed to validate the proposed method.(i) Case A.Consider
(84)
A
0
=
[
a
11
a
12
a
21
a
22
]
=
[
1
-
3
5
-
6
]
,
A
i
=
[
b
11
b
12
b
21
b
22
]
=
[
-
4
4
2
-
3
]
,
eig
(
A
-
)
=
[
-
2.0
,
-
10.0
]
,
c
3
2
-
4
c
5
=
33.0
>
0
,
where A
- is defined by (26).It is obvious that the conditions described by (42)-(43) are suitable for this case.Via the tracing method proposed by [23], we can obtain the system eigenvalue locus. Figure 9 depicts the curve of the second eigenvalue with τ increasing (the first one is always on the real axis, so it is not shown in the figure). There is an OEB bifurcation at point Q with τ
=
0.0476
s. According to (38) and (41), the curves of τ
-
z
1 and τ
-
z
2 can be depicted in Figure 10. Solid line represents that the corresponding z
i, i
=
1,2, is real, while dotted line represents that it is complex. Some results nearby point Q are listed in Table 2. It can be found that condition (42) holds at point Q with (
τ
,
z
1
)
R
=
(
0.0476
s
,
-
1
/
e
) and it indicates the occurrence of OEB bifurcation precisely.Table 2
z
1
,
z
2 values near point Q in Figure 10.
τ
Re
(
z
1
)
Re
(
z
1
)
+
1
/
e
Re
(
z
2
)
Re
(
z
2
)
+
1
/
e
0.04758
−0.3675333
0.0003462
−0.0858096
0.2820698
0.04759
−0.3676635
0.0002159
−0.0865306
0.2813488
0.04760
−0.3678304
0.0000491
−0.0879570
0.2799224
0.04761
−0.3679437
−0.0000643
−0.0883511
0.2795284
0.04762
−0.3680362
−0.0001567
−0.0883415
0.2795379Figure 9
Locus of the second eigenvalue withτ increasing.Figure 10
Curves ofτ
-
z
1 and τ
-
z
2.(ii) Case B.Consider
(85)
A
0
=
[
a
11
a
12
a
21
a
22
]
=
[
1
2
5
-
6
]
,
A
i
=
[
b
11
b
12
b
21
b
22
]
=
[
-
4
-
1
2
-
3
]
,
eig
(
A
-
)
=
[
-
2.0
,
-
10.0
]
,
c
3
2
-
4
c
5
=
-
7.0
<
0
,
where A
- is defined by (26).It can be asserted that the conditions described by (66) are appropriate for this case. Similar to Case A, the eigenvalue locus with τ increasing can be obtained via the tracing method expressed by [23]. The eigenvalue locus is shown in Figure 11, while real and imaginary parts of the eigenvalues are shown in Figure 12. Note that two real eigenvalues at the beginning move close to each other on the real axis. They encounter at point A and then turn to a pair of conjugated complex eigenvalues. According to Definition 1, there is no OEB bifurcation in the system since the number of eigenvalues does not increase here.Figure 11
Eigenvalue locus of Case B.Real and imaginary parts of two eigenvalues.
(a)
Real parts of two eigenvalues withτ increasing
(b)
Imaginary parts of two eigenvalues withτ increasingFigure13 plots the curves of τ
-
z
2, τ
-
z
3, and τ
-
z
4. Similarly, solid and dotted lines are used to represent the real and complex values of z
i, i
=
2,3
,
4. It is obvious that any condition identified by (66) does not hold here, so there is no OEB bifurcation in the system.Figure 13
Curves ofτ
-
z
2, τ
-
z
3, and τ
-
z
4.Case studies above demonstrate that conditions expressed by (42)-(43) or (66) are correct. And, the method proposed can correctly indicate the occurrence of OEB bifurcation in the 2-dim time-delayed system.


## 5.3. 3-Dim System

Consider a 3-dim system defined by (1) with A
0, A
i as follows:
(86)
A
0
=
[
-
5
1
-
1
1
-
4
2
0
0
-
2
]
,
A
i
=
[
2
1
0
1
-
1
-
1
1
3
-
8
]
.
It can be found that the system have three real eigenvalues with τ
=
0:
(87)
eig
(
A
-
)
=
[
-
1.8692
,
-
5.5749
,
-
10.5559
]
,
where A
-
=
A
0
+
A
i.Similar to Section5.2, we can get the eigenvalue locus shown in Figure 14, where the first and second eigenvalues always move on the real axis, while the third eigenvalue splits into a pair of conjugated complex eigenvalues at point B with τ
=
0.0431
s so as to yield an OEB bifurcation.Figure 14
Eigenvalue locus with time delay increasing.Table3 gives the coefficients in (69) and (70). According to the definition of z
1
,
z
2
,
z
5 by (76)–(78), we can depict the curves of τ
-
z
1, τ
-
z
2, and τ
-
z
5 in Figure 15. The τ
-
z
5 curve (red line) hits the horizontal line of -
1
/
e at τ
=
0.0431
s and condition (78) holds which indicate the occurrence of the OEB bifurcation.Table 3
Coefficients ofC
E
3
(
s
,
D
) given by (69) and (70).
a
1
a
2
a
3
a
4
a
5
a
6
a
7
a
8
a
9
11
37
38
7
60
117
−8
−16
−29
c
1
c
2
c
3
c
4
c
5
c
6
c
7
c
8
c
9
8.2422
15.7565
−0.5508
−0.7162
−3.8406
2.7578
7.5508
1.4870
5.4538Figure 15
Curves ofτ
-
z
1, τ
-
z
2, and τ
-
z
5.It should be mentioned that the transformation from (69) to (70) is not unique. Table 4 shows two sets of coefficients in (70). The corresponding τ
-
z
1, τ
-
z
2, and τ
-
z
5 curves are illustrated, respectively, by Figures 16 and 17. It is notable that in Figure 16 or Figure 17, the τ
-
z
1 curve (black line) collides with the horizontal line of -
1
/
e at τ
=
0.0431
s. The condition of (76) holds at this point so as to indicate occurrence of the OEB bifurcation.Table 4
Two sets of coefficients in (70).
c
1
c
2
c
3
c
4
c
5
c
6
c
7
c
8
c
9
Set I
6.4863
13.8753
9.2544
32.8516
12.8635
4.5137
−2.2544
6.1525
24.6287
Set II
7.2715
7.6935
5.2964
27.8647
−17.0229
3.7285
1.7036
−2.1947
−9.3149Figure 16
Curves ofτ
-
z
1, τ
-
z
2, and τ
-
z
5 with set I of Table 4.Figure 17
Curves ofτ
-
z
1, τ
-
z
2, and τ
-
z
5 with set II of Table 4.From the case study above, it is verified that conditions expressed by (76)–(78) can indicate the occurrence of OEB bifurcation in the 3-dim time-delayed LTI system.


## 6. Conclusions

In this paper, a method based on LambertW function theory is presented to determine the oscillation emergence bifurcation (OEB) in time-delayed LTI system with single lag. Conditions of OEB’s occurrence in 1-dim, 2-dim, 3-dim, and higher dimension systems are derived, respectively. And, case studies confirm the correctness and preciseness of the proposed method. Work of this paper is beneficial to explore the complex dynamic phenomena existing in the time-delayed system. However, the proposed method is only valid for the time-delayed system with single lag. How to expand it to the multiple time-delayed dynamic systems will be further discussed in the future.



---

*Source: 823937-2014-07-03.xml*