Delay-dependent stability analysis by using delay-independent integral evaluation

Automatica 70 (2016) 153–157

Contents lists available at ScienceDirect

Automatica journal homepage: www.elsevier.com/locate/automatica

Technical communique

Delay-dependent stability analysis by using delay-independent integral evaluation✩ Qi Xu a , Gabor Stepan b , Zaihua Wang a,1 a

State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, 210016 Nanjing, China

b

Department of Applied Mechanics, Budapest University of Technology and Economics, H-1521 Budapest, Hungary

article

info

Article history: Received 13 May 2015 Received in revised form 9 November 2015 Accepted 10 March 2016

Keywords: Time delay Stability Retarded Neutral Definite integral Delay-dependent Delay-independent

abstract For the stability analysis of time-delay systems, the available methods usually require the exact evaluation of some quantities. The definite integral stability method, originated from the Argument Principle or the Cauchy Theorem, is effective because it only requires a rough estimation of the testing integral over a finite interval to judge stability. However, no general rule is given in the literature for properly choosing the upper limit of the testing integral. In this paper, two simple algorithms are presented for finding the parameter-dependent critical upper limit and a parameter-independent upper limit without any restriction on the number of time delays. These results improve and complete the definite integral stability method. As illustrated by the numerical examples, the proposed algorithms work effectively and accurately. © 2016 Elsevier Ltd. All rights reserved.

1. Introduction Time-delay systems are described by delay differential equations (DDEs) (Erneux, 2009; Michiels & Niculescu, 2007; Niculescu, 2001; Stepan, 1989). The present study is restricted to autonomous Retarded DDEs (RDDEs) and Neutral DDEs (NDDEs) only. Time delays come usually from controllers, filters, actuators, or the contact problems of pure mechanical systems. In control applications, we mention here the car following problem modeled by RDDEs, where the delay is the human driver’s response time or the signal communication and processing time in autonomous vehicles (Brackstone & McDonald, 1999; Campbell, Egerstedt, How, & Murray, 2010; Orosz, Wilson, & Stepan, 2010); NDDE is used in modeling and reducing the sway of container cranes with a delayed

✩ This work was supported by National Natural Science Foundation of China under Grant 11372354, Hungarian National Development Agency under Grant TÉT12-CN-1-2012-0012, Hungarian Scientific Research Foundation OTKA under grant number K101714, Funding of Jiangsu Innovation Program for Graduate Education under Grant CXZZ13-0145. The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Fen Wu under the direction of Editor André L. Tits. E-mail addresses: [email protected] (Q. Xu), [email protected] (G. Stepan), [email protected] (Z. H. Wang). 1 Corresponding author.

http://dx.doi.org/10.1016/j.automatica.2016.03.028 0005-1098/© 2016 Elsevier Ltd. All rights reserved.

controller (Masoud & Nayfeh, 2003; Nayfeh, Masoud, & Nayfeh, 2011; Zhang, Wang, & Hu, 2012). As an example for a pure mechanical system with delay originated in the contact of tool and workpiece, RDDE is used to model the machine tool vibration of the cutting process (Long & Balachandran, 2007; Stepan, 2001). The possible negative effect of the delay on stability and performance has been a key issue in engineering applications. Lots of methods, criteria, or algorithms are available for the stability analysis of DDEs (Erneux, 2009; Kuang, 1993; Michiels & Niculescu, 2007; Niculescu, 2001; Olgac & Sipahi, 2002; Stepan, 1989). Lyapunov functionals and Lyapunov–Krasovskii functionals have been used for stability analysis of DDEs applicable also for global stability analysis (Gu, 2010; Gu & Liu, 2009; Niculescu, 2001). By rewriting the discrete-delay DDEs into distributed-delay DDEs, Lyapunov–Krasovskii functionals provide sufficient stability conditions of DDEs (Gu & Niculescu, 2000; Ivanescu, Niculescu, Dugard, Dion, & Verriest, 2003; Niculescu, 1999). In order to obtain more dedicated results, methods based on the analysis of characteristic functions are preferred. Among the first ones is the Pontryagin method that goes back to 1942 (Pontryagin, 1942). The method of stability switch is effective in finding the stable intervals for a given parameter (Kuang, 1993). A different version of the stability switch method is based on introducing Rekasius substitution, where the critical conditions can be studied with polynomials (Olgac, 2004; Olgac & Sipahi, 2002; Sipahi &

154

Q. Xu et al. / Automatica 70 (2016) 153–157

Olgac, 2006). Although the Nyquist criterion works effectively in some applications (Fu, Olbrot, & Polis, 1989, 1991), the Nyquist plot may not be good in judging stability because it may intersect itself in an intricate way, especially in case of oscillatory systems with multiple delays (Abdallah, Dorato, Benites-Read, & Byrne, 1993). The Argument Principle based stability criteria proposed and developed in Hassard (1997), Kolmanovskii and Myshkis (1999), Stepan (1989) and Xu and Wang (2014) enable one to calculate analytically the number of characteristic roots in the right-half complex plane (unstable roots, for short) for RDDEs and most of NDDEs. To do this, Hassard (1997) and Stepan (1989) suggest the calculation of all the finite number of roots with positive real parts of a transcendental real function associated with the complex characteristic function, instead of calculating the corresponding improper integral. The definite integral method works effectively for RDDEs (Kolmanovskii & Myshkis, 1999) by changing the corresponding improper integral into a proper one, and it is proved to be applicable also for most of NDDEs (Xu & Wang, 2014). A special and useful feature of the definite integral method is that it only requires a rough estimation of the testing integral. In this way, it is of less complexity and its computational cost is low with a properly chosen upper limit of the integral. However, no general rules are given for estimating this upper limit. This paper aims at proposing two general algorithms for choosing such an upper limit. With these, the number of unstable roots is calculated easily, so that the stability can be judged. The rest of the paper is organized as follows. In Section 2, a brief introduction of the stability criteria in definite integral form is given for DDEs. Then in Section 3, the algorithm for finding the critical upper limit of the definite integral is presented and illustrated with an example. In Section 4, the proposed algorithm is generalized to parameter-independent stability tests, and it is used to plot stability charts of NDDEs in an example with multiple delays. Finally in Section 5, some concluding remarks are drawn.

Let us consider a linear time-delay system described by m 

Ni x˙ (t − τi ) = Ax(t ) +

i =1

m 

Bi x(t − τi )

(1)

where x ∈ Rn , A, Bi , Ni ∈ Rn×n . The characteristic function of Eq. (1) is given in the form n 

αi (e−λτ1 , . . . , e−λτm )λn−i

(2)

i=0

where αi (z1 , . . . , zm ), i = 0, 1, . . . , n are real polynomials with respect to z1 = e−λτ1 , . . . , zm = e−λτm . Eq. (1) is a RDDE when Nk = 0 for all k = 1, 2, . . . , m, that is, α0 (z1 , . . . , zm ) ≡ 0, while it is a NDDE when at least one Nk ̸= 0 for some k = 1, 2, . . . , m, that is, α0 (z1 , . . . , zm ) ̸= 0. A trivial solution of a RDDE is asymptotically stable in Lyapunov sense if and only if all the characteristic roots stay in the open left half of the complex plane (Kuang, 1993). For a NDDE, further condition is needed that these roots are uniformly bounded away from the imaginary axis, because they may have accumulation points on the imaginary axis (Kuang, 1993). Assume that f (λ) has no roots on the imaginary axis. Let N be the number of all the unstable characteristic roots of Eq. (2) located in the right half complex plane, and let ∆C denote the argument change over the contour C shown in Fig. 1, which consists of {C1 : λ = Reθ i | θ ∈ (−π /2, π /2)} and {C2 : λ = ω i| ω ∈ (−R, R)}. Then the Argument Principle or Cauchy Theorem gives

N = lim

R→+∞

def

F (T1 , T2 ) =

T2



 ℜ

f ′ (ω i)



f (ω i)

T1

dω.

(4)

For a RDDE, Eq. (3) leads to Kolmanovskii and Myshkis (1999)

N =

n 2

−

1

π

lim F (0, T ).

(5)

T →+∞

Hence simplified from Eq. (5) using a definite integral, N = 0 if and only if there is a sufficient large T > 0 such that F (0, T ) > (n − 1)π /2, see Kolmanovskii and Myshkis (1999). For a NDDE, the limit in (5) does not converge. However, under the following assumption (Hale & Lunel, 2002) sup

ℜ(λ)>0, |λ|→∞

  α0 (e−λτ1 , . . . , e−λτm ) < 1,

(6)

1 ∆C arg(f (λ)) = lim R →+∞ 2π 2π i

Lemma 1. Assume that f (λ) has no roots on the imaginary axis, and condition (6) holds, then there exists a sufficiently large T0 > 0, such that for all T ≥ T0



i =1

f (λ) = λn +

Let ℜ(z ) denote the real part of complex number z, and define a testing integral F (T1 , T2 ) as

a lemma is derived in Xu and Wang (2014) from Eq. (3):

2. Problem statement

x˙ (t ) +

Fig. 1. The contour C of the line integral in Eq. (3).

f ′ (λ)

 C

f (λ)

dλ.

(3)

N ∈

−

F (0, T )

π

+

n−1 2

,−

F (0, T )

π

+

n+1



2

.

(7)

Because the length of the interval given in (7) is less than 1, the exact number N can be located in the interval (7) using an arbitrary T that is slightly larger than T0 . With the use of standard mathematical softwares, N can be calculated simply by using the round off command

 N = round

n 2

−

F (0, T )



π

.

(8)

When N = 0, system (1) is asymptotically stable. Thus, under assumption (6), the key issue for the stability analysis is to develop some effective algorithms for finding an upper limit T in the testing integral F (0, T ). 3. Parameter-dependent critical upper limit This section is devoted to finding a critical upper limit T0 of the testing integral, with which Eq. (8) gives the correct integer N for any T > T0 . The key idea is to introduce two real functions by using the real and imaginary parts of i−n f (ω i), given by R(ω) + S (ω) i = i−n f (ω i).

(9)

Q. Xu et al. / Automatica 70 (2016) 153–157

This was first used in Hassard (1997) for RDDEs to present the criterion in a unified form for f (ω i) independently whether the highest degree n in f (ω i) is even or odd. Then R(ω) = (1 + b0 (ω))ωn + b1 (ω)ωn−1 + · · · , S (ω) = c0 (ω)ωn + c1 (ω)ωn−1 + · · · , where bi (ω), ci (ω), i = 0, 1, . . . , n, are polynomials with respect to cos(ωτj ), sin(ωτj ), j = 1, . . . , m. To understand the key features of R(ω) and S (ω) graphically, let us consider the plot of i−n f (ω i) defined by

{(R(ω), S (ω)) : 0 6 ω < +∞}.

(10)

It starts at the point (R(0), S (0)), which is on the real/imaginary axis when n is even/odd, and it ends in the open right half complex plane due to 1 + b0 (ω) > 0 that is a consequence of condition (6). From Eq. (9), straightforward calculation shows that d dω

arctan

S (ω) R(ω)

 =ℜ

f ′ (ω i) f (ω i)



.

(11)

Hence, as ω varies from 0 to T , the value of F (0, T ) in Eq. (8) is the argument change of the curve (10) from (R(0), S (0)) to (R(T ), S (T )). Theorem 1. Assume that f (λ) has no roots on the imaginary axis, and condition (6) holds. Let T0 (τ1 , . . . , τm ) be the maximal positive root of R(ω) (if exists), then for all T > T0 (τ1 , . . . , τm ), the integer N satisfies Eq. (7).

155

Hence, according to (14), (18), for any T > T0 , one has F (0, T ) ∈



(n − 2N − 1)π (n − 2N + 1)π , 2



2

.

(19)

Case (ii). The crossing of the plot (10) at ω = T0 is on the negative half of the imaginary axis, namely R(T0 ) = 0, S (T0 ) < 0, and Eqs. (4) and (11) yield 0 < F (T0 , T ) =



T T0

d dω

arctan

S (ω) R(ω)

dω < π .

(20)

The same inequality (19) can be proved similarly from Eq. (20) as done for Case (i). Eq. (19) is equivalent to Eq. (7). This completes the proof.  Remark 1. If R(ω) has no positive root, one can choose T0 = 0 and prove N = round(n/2) as in Theorem 1. Remark 2. For any T that is slightly smaller than the maximal positive root T0 of R(ω), Eq. (8) generates an integer that is different from N . Therefore, T0 is critical indeed for choosing the upper limit. For example, Case (i) implies F (T0 , T ) > 0 for T = T0 −ε, 0 < ε ≪ 1, then Eqs. (12) and (18) yield F (0, T ) > (n − 2N + 1)π /2, which contradicts to Eq. (7), and the further calculation would yield an integer smaller than N .

Proof. Due to condition (6), the function R(ω) has finite number of parameter-dependent real zeros. Denote the maximal positive one (if exists) by T0 , then for all T > T0 , from Eq. (4) one has

Algorithm 1. Based on Theorem 1, the stability test of Eq. (1) with the parameter-dependent critical upper limit can be carried out in three steps:

F (0, T ) = F (0, T0 ) + F (T0 , T )

Step 1. Calculate the characteristic function f (λ) in Eq. (2), and the real part R(ω) of i−n f (ω i). Step 2. Find the largest positive root T0 of R(ω) = 0, choose an arbitrary number T > T0 , and calculate N by rounding off as in Eq. (8). Step 3. If N = 0, the trivial solution of the system is asymptotically stable, otherwise it is unstable.

(12)

Case (i). The crossing of the plot (10) at ω = T0 is on the positive half of the imaginary axis. In this case, R(T0 ) = 0, S (T0 ) > 0. Then Eqs. (4) and (11) generates

− π < F (T0 , T ) =

T



d dω

T0

arctan

S (ω) R(ω)

dω < 0

(13)

and

Example 1. Considering a second order NDDE

− π + F (0, T0 ) < F (0, T ) < F (0, T0 )

(14)

x¨ (t ) − 0.5x¨ (t − τ ) + 1000x˙ + 2000x˙ (t − τ ) + 500x(t ) = 0,

follows. Next, Lemma 1 implies that there must exist a large enough T ′ > T0 such that

whose characteristic function is

(n − 2N − 1)π

f (λ) = (1 − 0.5e−τ λ )λ2 + (1000 + 2000e−τ λ )λ + 500.

2

< F (0, T ′ ) <

(n − 2N + 1)π 2

.

(15)

From Eq. (13), T ′ > T0 gives

− π < F (T0 , T ) = ′



T′ T0

d dω

arctan

S (ω) R(ω)

dω < 0.

(16)

Now, the real part of i−2 f (ω i) takes the form of

Eqs. (15) and (16) lead to

(n − 2N − 1)π 2

< F (0, T0 ) <

(n − 2N + 3)π 2

.

(17)

Because the point (R(T0 ), S (T0 )) is located on the imaginary axis, while (R(0), S (0)) is on the real/imaginary axis when n is even/odd, F (0, T0 ) must be k1 π + π2 when n is even, or k2 π when n is odd, for some integers k1 , k2 . Moreover, for even n, the only number of the form k1 π + π2 located in ((n − 2N − 1)π /2, (n − 2N + 3)π /2) is (n − 2N + 1)π /2. For n is odd, (n − 2N + 1)π /2 is also the only number of the form k2 π . Thus Eq. (17) gives F (0, T0 ) =

(n − 2N + 1)π 2

.

When τ = 0.01, the above NDDE has 3 pairs of unstable characteristic roots: 18.260 ± 1440.546i, 40.900 ± 852.668i and 58.956 ± 281.235i, calculated by the numerical method proposed in Xu and Wang (2014).

(18)

R(ω) = (1 − 0.5 cos(0.01ω))ω2 − 2000 sin(0.01ω)ω − 500. The largest positive root of R(ω) is about 2038.166. Choose T = 2038.2, one has F (0, T ) = −17.273, hence Eq. (8) gives N = 6. With T = 2038.1, however, Eq. (8) leads to a wrong result N = 7. Moreover, by applying Algorithm 1 to each node of τ in a given interval τ ∈ [0, 0.011], N can be evaluated as shown in Fig. 2. The stable delay interval is about τ ∈ [0, 0.00079). In many applications, a guess 100 or 200 of the upper limit T0 can be large enough (Xu & Wang, 2014). But Example 1 shows that T0 can be very large, and a correct guess of T0 in judging the stability is needed indeed.

156

Q. Xu et al. / Automatica 70 (2016) 153–157

Fig. 2. The number N of the unstable roots for Example 1 over τ ∈ [0, 0.011].

4. Parameter-independent upper limit The critical upper limit used in Theorem 1 and Algorithm 1 depends on the delays and other system parameters, and it needs to be repeatedly calculated like in plotting Fig. 2. Thus, parameterindependent upper limit is preferable when some parameters, like the delays, are not fixed. Such an upper limit is also useful as an initial guess in finding the T0 in Step 2 of Algorithm 1 by using an iteration method like Newton–Raphson method. Theorem 2. Assume that f (λ) has no roots on the imaginary axis, and condition (6) holds, then there exists a constant T0 independent of the delays τ1 , . . . , τm such that for all T > T0 and all delays, N satisfies Eq. (7).

Fig. 3. The stability chart of Eq. (22) with respect to τ1 and τ2 with a = −0.8, b = −0.15, c = −1.5, d = 2, f = −5, g = 1, by choosing the delay independent upper limit T = 170.1 in Theorem 2. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Based on this theorem, Algorithm 2 can be proposed analogously to Algorithm 1 by replacing the maximal positive root of R(ω) with the maximal positive root of Rinf (ω) to be T0 .

For specified delay values, say τ1 = τ2 = 0.2, the upper limit T does not need to be as large as 170.1. In fact, Theorem 1 gives the critical upper limit T0 = 19.207, so T = 19.3 is a suitable upper limit for this case. Fig. 3 represents the stability chart of Eq. (22) in (τ1 , τ2 ) ∈ [0, 0.6] × [0, 0.65], repeatedly using Algorithm 2 to each combination of the meshed (τ1 , τ2 ). The stable region with N = 0 is colored in red, the unstable regions with N = 2 are colored in blue, and the unstable regions with N > 4 are in black. The results fit perfectly to the stability chart presented in Sipahi and Olgac (2006) using the cluster treatment of characteristic roots (CTCR). While CTCR is computationally more efficient for constructing stability charts in the space of delay parameters, the definite integral method requires less coding effort because it avoids the automatic selection of the intricate stable regions by means of root tendency calculations. In addition, the proposed method is efficient computationally for a system with fixed parameters independently from the number of time delays.

Example 2. Consider a first-order linear NDDE (Sipahi & Olgac, 2006)

5. Conclusion

Proof. Due to the harmonic periodicity of sin(ωτj ) and sin(ωτj ) in the coefficients of R(ω) in Eq. (9), each coefficient bi (ω) has its infimum inf bi (ω) with respect to τ1 , τ2 , . . . , τm , which becomes independent of all delay parameters. Thus for ω > 0, def

R(ω) > Rinf (ω) =(1 + inf b0 (ω))ωn + inf b1 (ω)ωn−1 + · · ·

(21)

Let T0 be the maximal positive root of Rinf (ω), then T0 is also independent of the delays and not less than the maximal positive root of R(ω). From Theorem 1, all T > T0 is suitable for locating N in Eq. (7), which completes the proof. 

x˙ (t ) − ax˙ (t − τ1 ) − bx˙ (t − τ2 )

= cx(t − τ1 ) + dx(t − τ2 ) + fx(t − τ1 − τ2 ) + gx(t ).

(22)

The characteristic equation is f (λ) = 0, where f (λ) = (1 − ae−τ1 λ − be−τ2 λ )λ

− ce−τ1 λ − de−τ1 λ − f e−(τ1 +τ2 )λ − g .

(23)

Here n = 1, and one has R(ω) = 1 − a cos(τ1 ω) − b cos(τ2 ω) ω + c sin(τ1 ω)





+ d sin(τ2 ω) − f sin(τ1 ω + τ2 ω),

(24)

Rinf (ω) = (1 − |a| − |b|)ω − |c | − |d| − |f | = 0.

(25)

and

With fixed a = −0.8, b = −0.15, c = −1.5, d = 2, f = −5, g = 1 as used in Sipahi and Olgac (2006), Eq. (25) becomes 0.05 ω = 8.5, which yields T0 = 170. Thus, we can choose T = 170.1 as the upper limit for any combination of τ1 and τ2 .

This paper presents two theorems and two algorithms for finding a sharp estimation of the upper limit of the definite integral method for the stability analysis of DDEs. The stability analysis of relevant classes of NDDEs is addressed. The upper limit of the testing integrals can be chosen as any number that is larger than the maximal root of a real function associated with the characteristic function under critical stability condition. It can be either parameter-dependent or parameter-independent. One advantage of the proposed algorithm is its easy implementation in computer algebra languages: the number of unstable characteristic roots can be determined simply by using the round off commands of the testing integral in any commercial mathematical software. In this way, the stability or instability can be determined directly, and also the kinds of possible bifurcations can be identified like Hopf, double Hopf, saddle–node, and so on, in the corresponding nonlinear DDEs. Another advantage of the proposed algorithm is that its implementation does not get more complex in the presence of many discrete delay parameters, even in neutral DDE systems.

Q. Xu et al. / Automatica 70 (2016) 153–157

Numerical examples show that the proposed schemes work efficiently and accurately. References Abdallah, C.T., Dorato, P., Benites-Read, J., & Byrne, R. (1993). Delayed positive feedback can stabilize oscillatory systems. In American control conference (pp. 3106–3107). Brackstone, M., & McDonald, M. (1999). Car-following: a historical review. Transportation Research Part F: Traffic Psychology and Behaviour, 2(4), 181–196. Campbell, M., Egerstedt, M., How, J. P., & Murray, R. (2010). Autonomous driving in urban environments: approaches, lessons and challenges. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 368(1928), 4649–4672. Erneux, T. (2009). Applied delay differential equations. London: Springer-Verlag. Fu, M., Olbrot, A. W., & Polis, M. (1989). Robust stability for time-delay systems: the edge theorem and graphical tests. IEEE Transactions on Automatic Control, 34(8), 813–820. Fu, M., Olbrot, A. W., & Polis, M. (1991). The edge theorem and graphical tests for robust stability of neutral time-delay systems. Automatica, 27(4), 739–741. Gu, K. (2010). Stability problem of systems with multiple delay channels. Automatica, 46(4), 743–751. Gu, K., & Liu, Y. (2009). Lyapunov–Krasovskii functional for uniform stability of coupled differential-functional equations. Automatica, 45(3), 798–804. Gu, K., & Niculescu, S.-I. (2000). Additional dynamics in transformed time-delay systems. IEEE Transactions on Automatic Control, 45(3), 572–575. Hale, J. K., & Lunel, S. M. V. (2002). Strong stabilization of neutral functional differential equations. IMA Journal of Mathematical Control and Information, 19(1 and 2), 5–23. Hassard, B. (1997). Counting roots of the characteristic equation for linear delaydifferential systems. Journal of Differential Equations, 136(2), 222–235. Ivanescu, D., Niculescu, S.-I., Dugard, L., Dion, J.-M., & Verriest, E. I. (2003). On delaydependent stability for linear neutral systems. Automatica, 39(2), 255–261. Kolmanovskii, V., & Myshkis, A. (1999). Introduction to the theory and applications of functional differential equations. Dordrecht: Kluwer Academic Publishers. Kuang, Y. (1993). Delay differential equations: with applications in population dynamics. London: Academic Press.

157

Long, X. H., & Balachandran, B. (2007). Stability analysis for milling process. Nonlinear Dynamics, 49(3), 349–359. Masoud, Z. N., & Nayfeh, A. H. (2003). Sway reduction on container cranes using delayed feedback controller. Nonlinear dynamics, 34(3–4), 347–358. Michiels, W., & Niculescu, S.-I. (2007). Advances in design and control, Stability and stabilization of time-delay systems. Philadelphia: Society for Industrial and Applied Mathematics. Nayfeh, A., Masoud, Z., & Nayfeh, N. (2011). A smart sway controller for cranes from theory to laboratory to industry. In Vibration problems ICoVP (pp. 14–29). Niculescu, S.-I. (1999). A model transformation class for delay-dependent stability analysis. In American control conference (pp. 314–318). Niculescu, S.-I. (2001). Delay effects on stability: a robust control approach. London: Springer-Verlag. Olgac, N. (2004). A practical method for analyzing the stability of neutral type LTItime delayed systems. Automatica, 40(5), 847–853. Olgac, N., & Sipahi, R. (2002). An exact method for the stability analysis of timedelayed linear time-invariant (LTI) systems. IEEE Transactions on Automatic Control, 47(5), 793–797. Orosz, G., Wilson, R. E., & Stepan, G. (2010). Traffic jams: dynamics and control. Philosophical Transactions of the Royal Society of London A (Mathematical and Physical Sciences), 368(1928), 4455–4479. Pontryagin, L. S. (1942). On the zeros of some elementary transcendental functions. Izvestiya Rossiiskoi Akademii Nauk. Seriya Matematicheskaya, 6(3), 115–134. (Russian) English translation (1955) in American Mathematical Society Translations, Vol. 2, pp. 95–110. Sipahi, R., & Olgac, N. (2006). Complete stability analysis of neutral-type first order two-time-delay systems with cross-talking delays. SIAM Journal on Control and Optimization, 45(3), 957–971. Stepan, G. (1989). Retarded dynamical systems: stability and characteristic functions. Longman Scientific & Technical Essex. Stepan, G. (2001). Modelling nonlinear regenerative effects in metal cutting. Philosophical Transactions of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 359(1781), 739–757. Xu, Q., & Wang, Z. H. (2014). Exact stability test of neutral delay differential equations via a rough estimation of the testing integral. International Journal of Dynamics and Control, 2(2), 154–163. Zhang, L., Wang, H. L., & Hu, H. Y. (2012). Symbolic computation of normal form for Hopf bifurcation in a neutral delay differential equation and an application to a controlled crane. Nonlinear Dynamics, 70(1), 463–473.