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

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.

(ProQuest: ... denotes non-US-ASCII text omitted.)

Academic Editor:Xinghuo Yu

Key Laboratory of Smart Grid of Ministry of Education, Tianjin University, Tianjin 300072, China

Received 23 January 2014; Accepted 21 May 2014; 3 July 2014

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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, Figure 1 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.

[figure omitted; refer to PDF]

The contents of this paper are organized into six sections. Following this introduction, definition of OEB bifurcation is put forward in Section 2 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: [figure omitted; refer to PDF] where _{A 0} , _{A i} ∈ ^{R n × n} ; x ∈ ^{R n} is the vector of the state variables; τ ...5; 0 is the constant of time delay; [varphi] ( 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 Σ : [figure omitted; refer to PDF] where s ∈ C is the eigenvalue to be determined. The eigenvalue spectrum of system Σ is denoted by Γ : [figure omitted; refer to PDF] 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} : [figure omitted; refer to PDF] where _{n R} , _{n I} are the eigenvalue numbers of _{Γ R} and _{Γ I} , respectively. Then the following (5) holds: [figure omitted; refer to PDF] where _{n τ} is the eigenvalue number of Γ .

It can be found that when τ = 0 , system Σ degenerates into the following LTI differential equation: [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] _{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 τ ...4; _{τ c} . And, after τ > _{τ c} , the numbers of _{Γ R} and _{Γ I} turn to ^{n R n} and ^{n I n} . If the following equations hold: [figure omitted; refer to PDF] there is an oscillation emergence bifurcation at τ = _{τ c} .

According to Definition 1, suppose there is a real eigenvalue s ∈ R in Γ at τ ...4; _{τ 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.

[figure omitted; refer to PDF]

3. Lambert W Function and OEB Bifurcation

Our aim is to derive a general method based on Lambert W 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 Lambert W function: [figure omitted; refer to PDF] where w , z ∈ C . Its inverse function is expressed as [figure omitted; refer to PDF]

Lambert W 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 ...5; - 1 / e . For - 1 / e ...4; z < 0 , there are two possible real values of W ( z ) as shown in Figure 3. We denote the branch satisfying - 1 ...4; W ( z ) by _{W 0} (solid line in Figure 3), and the branch satisfying W ( z ) ...4; - 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: [figure omitted; refer to PDF]

Figure 3: _{W 0} and _{W - 1} branches of Lambert W function.

[figure omitted; refer to PDF]

Supposed that z ∈ R in (9) and (10) is a system parameter and it varies with another parameter u ∈ R : [figure omitted; refer to PDF]

Substitute (12) into (10). It can be found that w is also a function of u : [figure omitted; refer to PDF]

Then, the u - 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: [figure omitted; refer to PDF]

Figure 4: Curve of z = _{f z} ( u ) .

[figure omitted; refer to PDF]

If u 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 that _{N 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 Figure 4: [figure omitted; refer to PDF] where _{a 0} , r ∈ R , and τ ...5; 0 . The value of s can be derived from (15) based on Lambert W function as follows: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Figure 5 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 of z = _{f z} ( τ ) .

[figure omitted; refer to PDF]

From Figures 4 and 5, it is asserted that the following formula holds at the OEB bifurcation points when _{f z} ( · ) ∈ R : [figure omitted; refer to PDF] And, if _{f z} ( · ) has complex values on one side of the OEB bifurcation point, the above condition will be modified as follows: [figure omitted; refer to PDF] 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 Lambert W 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: [figure omitted; refer to PDF] where τ ...5; 0 , a , b , x ∈ R , and b ...0; 0 . The characteristic equation of (20) is [figure omitted; refer to PDF] where D = ^{e - τ s} . It can be rewritten as [figure omitted; refer to PDF] According to (15), s can be solved from (22) based on Lambert W function: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Similar to Figure 5, the OEB bifurcation can be determined from (24) via condition (18).

4.2. 2-Dim System

Let us consider the following 2-dim system: [figure omitted; refer to PDF] where τ ...5; 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: [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] System characteristic equation can be derived as follows: [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] where ζ , y , w , v , z ∈ R . And, it can be expanded as follows: [figure omitted; refer to PDF]

Comparing (30) with (28), the following equations hold: [figure omitted; refer to PDF] And, the following result can be obtained: [figure omitted; refer to PDF] 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} ...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): [figure omitted; refer to PDF] where K ∈ C is a temporary variable.

Via the first formula in (35), the following equation holds: [figure omitted; refer to PDF] According to (15) and (16), with _{a 0,1} = - 1 / w , _{r 1} = K v - ζ , _{c 1} = τ , s has the following form: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Via the second formula in (35), the following equation holds: [figure omitted; refer to PDF]

Similarly, having defined _{a 0,2} = - 1 / y , _{r 2} = 1 / K - z , _{c 2} = τ , the following result can be obtained: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] where ζ , y , w , v , z , q ∈ R . Also, its expansion is given as follows: [figure omitted; refer to PDF] Comparing it with (28), the following equations hold: [figure omitted; refer to PDF]

Hence, [figure omitted; refer to PDF] To guarantee that ζ , y , w , v , z are real numbers, the following relationship must hold: [figure omitted; refer to PDF]

Since q can be selected arbitrarily, it is set as follows: [figure omitted; refer to PDF]

Substituting (53) into (51) yields [figure omitted; refer to PDF]

Since ^{c 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 [figure omitted; refer to PDF]

It 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: [figure omitted; refer to PDF]

Suppose that there exists a real variable Q ∈ R such that [figure omitted; refer to PDF]

Inserting it into (59) yields [figure omitted; refer to PDF]

Defining _{a 0,3} = sign ... ( _{P n} ( s ) ) , _{r 3} = K v - ζ , _{c 3} = Q τ , eigenvalue s is expressed as [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

Defining _{a 0,4} = sign ... ( _{P d} ( s ) ) · K q , _{r 4} = w / K q , _{c 4} = ( Q - 1 ) τ , the similar result from (61) is obtained: [figure omitted; refer to PDF] where [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF]

4.3. High-Dimensional Systems

Consider system Σ given by (1) with n ...5; 3 . Its characteristic equation has the following general form: [figure omitted; refer to PDF] where D = ^{e - τ s} and _{p n - i} ( s ) is a polynomial of s : [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF]

It can be converted into the following form so as to utilize the previous results: [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] Via the first formula in (71), [figure omitted; refer to PDF] It is obvious that ^{CE 2 [variant prime]} ( 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), [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF]

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 Sections 4.1 and 4.2 can be applied. For simplification, it is denoted as [figure omitted; refer to PDF]

As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system: [figure omitted; refer to PDF] where (76) and (77) are concluded from (72) with conditions expressed by (42)-(43) when ^{c 3 2} - 4 _{c 5} ...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 , [figure omitted; refer to PDF]

Setting _{CE 3} ( s , D ) = K , (79) can be rewritten as [figure omitted; refer to PDF]

The first formula is a 3-dim case, and the second is a 2-dim case. Hence, [figure omitted; refer to PDF] 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 with n ...5; 3 can be obtained: [figure omitted; refer to PDF]

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: [figure omitted; refer to PDF]

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 omitted; refer to PDF]

Figure 7: Curve of an eigenvalue.

[figure omitted; refer to PDF]

Similar to Figure 5, 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

[figure omitted; refer to PDF]

(b) Simulation results with τ = 0.02 s and τ = 0.03 s

[figure omitted; refer to PDF]

5.2. 2-Dim System

The following two cases are employed to validate the proposed method.

(i) Case A. Consider [figure omitted; refer to PDF] 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.2795379

Figure 9: Locus of the second eigenvalue with τ increasing.

[figure omitted; refer to PDF]

Figure 10: Curves of τ - _{z 1} and τ - _{z 2} .

[figure omitted; refer to PDF]

(ii) Case B. Consider [figure omitted; refer to PDF] 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.

[figure omitted; refer to PDF]

Real and imaginary parts of two eigenvalues.

(a) Real parts of two eigenvalues with τ increasing

[figure omitted; refer to PDF]

(b) Imaginary parts of two eigenvalues with τ increasing

[figure omitted; refer to PDF]

Figure 13 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} .

[figure omitted; refer to PDF]

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: [figure omitted; refer to PDF] It can be found that the system have three real eigenvalues with τ = 0 : [figure omitted; refer to PDF] where A - = _{A 0} + _{A i} .

Similar to Section 5.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.

[figure omitted; refer to PDF]

Table 3 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 of C _{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.4538

Figure 15: Curves of τ - _{z 1} , τ - _{z 2} , and τ - _{z 5} .

[figure omitted; refer to PDF]

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

Figure 16: Curves of τ - _{z 1} , τ - _{z 2} , and τ - _{z 5} with set I of Table 4.

[figure omitted; refer to PDF]

Figure 17: Curves of τ - _{z 1} , τ - _{z 2} , and τ - _{z 5} with set II of Table 4.

[figure omitted; refer to PDF]

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 Lambert W 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.

Acknowledgment

This work is supported by the Natural Science Foundation of China under Grant nos. 51277128 and 51377117.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

References

[1] J. K. Hale Theory of Functional Differential Equations , World Publishing Corporation, Beijing, China, 2003.

[2] L. Dugard, E. I. Verriest Stability and Robust Control of Time Delay Systems , Springer, New York, NY, USA, 1997.

[3] F. S. Gentile, J. L. Moiola, E. E. Paolini, "Nonlinear dynamics of internet congestion control: a frequency-domain approach," Communications in Nonlinear Science and Numerical Simulation , vol. 19, no. 4, pp. 1113-1127, 2014.

[4] T. Ito, H. Ohashi, T. Tamaki, T. Yamaguchi, "Mathematical modeling of molecular recognition by an ion-gating membrane oscillator," Journal of Membrane Science , vol. 448, pp. 231-239, 2013.

[5] Y. Yu, H. Tang, X. Han, Q. Bi, "Bursting mechanism in a time-delayed oscillator with slowly varying external forcing," Communications in Nonlinear Science and Numerical Simulation , vol. 19, no. 4, pp. 1175-1184, 2014.

[6] S. Xu, J. Bao, "Distributed control of plant-wide chemical processes with uncertain time-delays," Chemical Engineering Science , vol. 84, pp. 512-532, 2012.

[7] S. N. Rader, M. L. Reagan, B. Janoiko, J. E. Johnson, "Human- in-the-loop operations over time delay: NASA analog missions lessons learned," in Proceedings of the 43rd International Conference on Environmental Systems, pp. 1-6, Vail, Colo, USA, 2013.

[8] H. Jia, Y. Mu, Y. Qi, "A statistical model to determine the capacity of battery-supercapacitor hybrid energy storage system in autonomous microgrid," International Journal of Electrical Power & Energy Systems , vol. 54, pp. 516-524, 2014.

[9] H. Jia, X. Yu, Y. Yu, C. Wang, "Power system small signal stability region with time delay," International Journal of Electrical Power and Energy Systems , vol. 30, no. 1, pp. 16-22, 2008.

[10] X. Yu, H. Jia, C. Wang, "CTDAE & CTODE models and their applications to power system stability analysis with time delays," Science China Technological Sciences , vol. 56, no. 5, pp. 1213-1223, 2013.

[11] Y. He, M. Wu, J. She, "Delay-dependent exponential stability of delayed neural networks with time-varying delay," IEEE Transactions on Circuits and Systems II: Express Briefs , vol. 53, no. 7, pp. 553-557, 2006.

[12] J. YiLang, J. Tao, J. HongJie, D. ChaoYu, "A novel LMI criterion for power system stability with multiple time-delays," Science China Technological Sciences , vol. 57, no. 7, pp. 1392-1400, 2014.

[13] R. H. Gielen, M. Lazar, S. V. Rakovic, "Necessary and sufficient Razumikhin-type conditions for stability of delay difference equations," IEEE Transactions on Automatic Control , vol. 58, no. 10, pp. 2637-2642, 2013.

[14] N. Olgac, R. Sipahi, "An exact method for the stability analysis of time-delayed linear time-invariant (LTI) systems," IEEE Transactions on Automatic Control , vol. 47, no. 5, pp. 793-797, 2002.

[15] Q. Sun, H. An, H. Jia, X. Yu, "An improved power system stability criterion with multiple time delays," in Proceedings of the IEEE General Meeting, vol. 1-8, pp. 786-792, Calgary Alberta, Calif, USA, July 2009.

[16] F. Rossi, P. Colaneri, R. Shorten, "Padé discretization for linear systems with polyhedral Lyapunov functions," IEEE Transactions on Automatic Control , vol. 56, no. 11, pp. 2717-2722, 2011.

[17] J. Hongjie, Y. Xiaodan, W. Chengshan, W. Wei, Z. Yuan, Z. Jinli, L. Peng, "Study on power system extended small signal stability region (DE-SSSR) in time delay space," in Proceedings of the 15th IEEE Mediterranean Electrotechnical Conference (MELECON '10), vol. 1, pp. 1569-1574, Valletta, Malta, April 2010.

[18] Y. Ding, W. Jiang, P. Yu, "Double Hopf bifurcation in delayed van der Pol-duffing equation," International Journal of Bifurcation and Chaos in Applied Sciences and Engineering , vol. 23, no. 1, 2013.

[19] X. Xu, H. Y. Hu, H. L. Wang, "Stability, bifurcation and chaos of a delayed oscillator with negative damping and delayed feedback control," Nonlinear Dynamics: An International Journal of Nonlinear Dynamics and Chaos in Engineering Systems , vol. 49, no. 1-2, pp. 117-129, 2007.

[20] F.-X. Wu, "Stability and bifurcation of ring-structured genetic regulatory networks with time delays," IEEE Transactions on Circuits and Systems I: Regular Papers , vol. 59, no. 6, pp. 1312-1320, 2012.

[21] T. Co, "Relay-stabilization and bifurcations of unstable SISO processes with time delay," IEEE Transactions on Automatic Control , vol. 55, no. 5, pp. 1131-1141, 2010.

[22] S. Guo, W. Jiang, "Global stability and Hopf bifurcation for Gause-type predator-prey system," Journal of Applied Mathematics , vol. 2012, 2012.

[23] X. Yu, H. Jia, C. Wang, "An eigenvalue spectrum tracing algorithm and its application in time delayed power systems," Automation of Electric Power Systems , vol. 36, no. 24, pp. 10-38, 2012.

[24] D. J. Jeffrey, D. E. G. Hare, R. M. Corless, "Unwinding the branches of the Lambert W function," The Mathematical Scientist , vol. 21, no. 1, pp. 1-7, 1996.

[25] R. M. Corless, G. H. Gonnet, D. E. G. Hare, D. E. Knuth, "On the Lambert W function," Advances in Computational Mathematics , vol. 5, no. 1, pp. 329-359, 1996.

[26] E. Jarlebring, T. Damm, "The Lambert W function and the spectrum of some multidimensional time-delay systems," Automatica , vol. 43, no. 12, pp. 2124-2128, 2007.

[27] C. Hwang, Y. Cheng, "A note on the use of the Lambert W function in the stability analysis of time-delay systems," Automatica , vol. 41, no. 11, pp. 1979-1985, 2005.

[28] Y. Sun, A. G. Ulsoy, "Solution of a system of linear delay differential equations using the matrix lambert function," in Proceedings of the 26th American Control Conference, pp. 2433-2438, Minneapolis, Minn, USA, June 2006.

[29] S. Yi, S. Yu, J. H. Kim, "Analysis of neural networks with time-delays using the Lambert W function," in Proceedings of the American Control Conference (ACC '11), pp. 3221-3226, San Francisco, Calif, USA, July 2011.