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

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 oscillationmode will emerge. Next, a method to determine OEB bifurcation in the timedelayed 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.


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][9][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][11][12][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][19][20][21][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  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  function to propose a mechanism to the occurrence of this bifurcation.As an available tool, Lambert  function [24,25] has been utilized to analyze some complex behaviors in nonlinear system including stability analysis of time-delayed system [26][27][28][29].It is applied to derive a general method to determine OEB bifurcation in this paper.
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  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  function.In Section 5, the presented approach is demonstrated by numerical simulations.Finally, Section 6 concludes and summarizes this paper.

Definition of OEB Bifurcation
In this paper, all studies are based on the following timedelayed LTI system with single lag: where A 0 , A  ∈  × ; x ∈   is the vector of the state variables;  ≥ 0 is the constant of time delay; () is the system trajectory in the range of [−, 0].It is well known that the following equation is the corresponding characteristic equation of system Σ: where  ∈  is the eigenvalue to be determined.The eigenvalue spectrum of system Σ is denoted by Γ: Theoretically, it can be calculated from (2) directly when  is given.Further, denote the set of the real eigenvalue by Γ  and the set of the complex eigenvalue with nonzero imaginary by Γ  : where   ,   are the eigenvalue numbers of Γ  and Γ  , respectively.Then the following (5) holds: where   is the eigenvalue number of Γ.
It can be found that when  = 0, system Σ degenerates into the following LTI differential equation: where A = A 0 + A  .And, the number of eigenvalues is ; that is,  =0 = .When  increases, the values of   ,   , and   will change accordingly.Therefore, they can be considered as functions of  as follows: (),   (), and () 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   > 0 for system Σ of (1).Its Γ  has    real eigenvalues and Γ  has    complex eigenvalues with  ≤   .And, after  >   , the numbers of Γ  and Γ  turn to    and    .If the following equations hold: there is an oscillation emergence bifurcation at  =   .
According to Definition 1, suppose there is a real eigenvalue  ∈  in Γ at  ≤   .It will split into a pair of conjugated eigenvalues  + ,  − as shown in Figure 2. It is manifested that the number of Γ  will decrease by one, the number of Γ  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  + ,  − will emerge after OEB bifurcation.

Lambert 𝑊 Function and OEB Bifurcation
Our aim is to derive a general method based on Lambert  function to determine OEB bifurcation in system Σ defined by (1).A brief introduction to Lambert  function is firstly given.
The following function is called Lambert  function: where ,  ∈ .Its inverse function is expressed as Lambert  function is not injective and the relation  is multivalued (except at 0).It has many solution branches in complex plane [25], which can be denoted by   ,  = 0, ±1, ±2, . ... If the attention is restricted to real-valued , then the relation is defined only for  ≥ −1/.For −1/ ≤  < 0, there are two possible real values of () as shown in Figure 3.We denote the branch satisfying −1 ≤ () by  0 (solid line in Figure 3), and the branch satisfying () ≤ −1 by  −1 (dashed line in Figure 3). 0 is referred to the principal branch of Lambert  function.It increases from  0 (−1/) = −1 to  0 (+∞) = +∞, while branch encounter at point  as shown in Figure 3, and the value of point  is given as follows: Supposed that  ∈  in ( 9) and ( 10) is a system parameter and it varies with another parameter  ∈ : Substitute ( 12) into (10).It can be found that  is also a function of : Then, the - curve can be plotted.An example is shown in Figure 4, where - curve crosses the dotted line of  = −1/ at some points which can be divided into two categories: If  is on the interval of (  ,   ), that is, the intervals above the line with  > −1/, (13) has real solutions.While if  is on the interval of (  ,  +1 ), that is, the grey intervals with  < −1/, (13) only has complex solutions with nonzero imaginary.
Recalling the definition of OEB, we can find that   is just a set of OEB bifurcations of the system expressed by (13).That is, at any point of   , a real eigenvalue will turn into a pair of conjugated complex eigenvalues with the increment of .
Here the following typical transcendental item in onedimensional (1-dim) delayed-system is treated as an example for further discussion to Figure 4: where  0 ,  ∈ , and  ≥ 0. The value of  can be derived from (15) based on Lambert  function as follows: where  From Figures 4 and 5, it is asserted that the following formula holds at the OEB bifurcation points when   (⋅) ∈ : And, if   (⋅) has complex values on one side of the OEB bifurcation point, the above condition will be modified as follows: where Re() is the real part of .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.

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

1-Dim System.
The following 1-dim system is considered in this subsection: where  ≥ 0, , ,  ∈ , and  ̸ = 0.The characteristic equation of ( 20) is where  =  − .It can be rewritten as According to (15),  can be solved from (22) based on Lambert  function: where Similar to Figure 5, the OEB bifurcation can be determined from (24) via condition (18).

2-Dim System.
Let us consider the following 2-dim system: where  ≥ 0,   ,   ∈ , ,  = 1, 2; It means that there are two negative and real eigenvalues at  = 0, so system given by ( 25) is stable without time delay.Supposed that  is an eigenvalue of (25) at  > 0 and the following equation holds: System characteristic equation can be derived as follows: With  1 = −( In order to determine the OEB bifurcation in the system, we reorganize (28) into the following form: where , , , V,  ∈ .And, it can be expanded as follows: Comparing ( 30) with ( 28), the following equations hold: And, the following result can be obtained: when  2 3 − 4 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  2  3 − 4 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): where  ∈  is a temporary variable.
Via the first formula in (35), the following equation holds: According to (15) and ( 16), with  0,1 = −1/,  1 = V − ,  1 = ,  has the following form: where Via the second formula in (35), the following equation holds: Similarly, having defined  0,2 = −1/,  2 = 1/−,  2 = , the following result can be obtained: where Equations ( 37) and ( 40) can be utilized to solve the values of  and  when  is given.They can be further taken into (38) and (41) to get the values of  1 and  2 .And if any one of the following equations holds, there is an OEB bifurcation: Obviously, if  is a complex eigenvalue,  defined by (34) and  1 ,  2 given by ( 38) and (41) are all complex values.They are all real values if  is a real eigenvalue.Hence, condition (19) is suitable for this case.
(ii) Consider  2 3 − 4 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: where , , , V,, ∈ .Also, its expansion is given as follows: Comparing it with (28), the following equations hold: Hence, To guarantee that , , , V,  are real numbers, the following relationship must hold: Since  can be selected arbitrarily, it is set as follows: Substituting ( 53) into (51) yields Since  2 3 − 4 5 < 0, it must be true that  5 > 0. So  and  are two real numbers.According to (49)-(51), , V,  are all real numbers.In order to use condition (19), (44) is written as 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: Suppose that there exists a real variable  ∈  such that Inserting it into (59) yields Defining  0,3 = sign(  ()),  3 = V − ,  3 = , eigenvalue  is expressed as where Defining  0,4 = sign(  ())⋅,  4 = /,  4 = (−1), the similar result from (61) is obtained: where Equations ( 40), (62), and (64) can be used to solve the values of , , and  when  is given.They can be further taken into (41), (63), and (65) to get the values of  2 ,  3 , and  4 .As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system: (66)

High-Dimensional Systems.
Consider system Σ given by (1) with  ≥ 3. Its characteristic equation has the following general form: where  =  − and  − () is a polynomial of : (i) 3-Dim Case ( = 3).Take 3-dim system as an example to explain the principle of proposed method.General form of the characteristic equation for  = 3 is given as below: It can be converted into the following form so as to utilize the previous results: where CE 2 (, ) is described by (28).And,   ∈ ,  = 1, 2, . . ., 9, are variables to be determined.Setting CE 2 (, ) = , (70) can be written as Via the first formula in (71), It is obvious that CE  2 (, ) 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), where ã = ( 9 −  6 )/( −  8 ), b = − 7 /( −  8 ).The result of 1-dim system in Section 4.1 can be applied to (73).And the following result can also be derived: 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 As a conclusion, any of the following conditions hold, and there is an OEB bifurcation in the system: where ( 76) and (77) are concluded from (72) with conditions expressed by (42)-(43) when  2 3 − 4 5 ≥ 0 (or conditions expressed by (66) when  2  3 − 4 5 < 0).
(ii) Case of High Dimension with  > 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  = 4, (80)  The first formula is a 3-dim case, and the second is a 2-dim case.Hence, 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 highdimensional systems with  ≥ 3 can be obtained: Then, a higher dimension system can be converted into two lower dimension systems so as to determine the occurrence of OEB bifurcation accordingly.

1-Dim System.
Consider the following 1-dim system: According to (24), the - 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.
Similar to Figure 5, point  with (, ) = (0.0173, −1/) 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.
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 and converges oscillatedly in the range of  > 0.0173.It is manifested that there is OEB bifurcation at  = 0.0173, which leads to a new oscillation mode and increases the number of system eigenvalues by one.

2-Dim System.
The following two cases are employed to validate the proposed method.
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  with  = 0.0476.According to (38) and (41), the curves of - 1 and - 2 can be depicted in Figure 10.Solid line represents that the corresponding   ,  = 1,2, is real, while dotted line represents that it is complex.Some results nearby point  are listed in Table 2.It can be found that condition (42) holds at point  with 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  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 13 plots the curves of - 2 , - 3 , and - 4 .Similarly, solid and dotted lines are used to represent the real and complex values of   ,  = 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.
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.

3-Dim
System.Consider a 3-dim system defined by (1) with A 0 , A  as follows: It can be found that the system have three real eigenvalues with  = 0: where A = A 0 + A  .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  with  = 0.0431 so as to yield an OEB bifurcation.
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 - 1 , - 2 , and - 5 curves are illustrated, respectively, by Figures 16 and 17.It is notable that in Figure 16 or Figure 17, the - 1 curve (black line) collides with the horizontal line of −1/ at  = 0.0431.The condition of (76) holds at this point so as to indicate occurrence of the OEB bifurcation.
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.4.

Conclusions
In this paper, a method based on Lambert  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, 3dim, 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.

1 Figure 1 :
Figure 1: A real eigenvalue splits into a pair of conjugated eigenvalues.

Figure 5
Figure 5 depicts the curve of - with  0 = −0.5,  = 1.0.It is evident that point  with  = 0.23196 in the figure is an OEB bifurcation of system (15) according to Definition 1 and Figure4.From Figures4 and 5, it is asserted that the following formula holds at the OEB bifurcation points when   (⋅) ∈ :

Figure 12 :Figure 13 :
Figure 12: Real and imaginary parts of two eigenvalues.

Figure 14 :
Figure 14: Eigenvalue locus with time delay increasing.