Subcritical Hopf bifurcations in a car-following model with reaction-time delay

Gábor Orosz , Gábor Stépán

Abstract

A nonlinear car-following model of highway traffic is considered, which includes the reaction-time delay of drivers. Linear stability analysis shows that the uniform flow equilibrium of the system loses its stability via Hopf bifurcations and thus oscillations can appear. The stability and amplitudes of the oscillations are determined with the help of normal-form calculations of the Hopf bifurcation that also handles the essential translational symmetry of the system. We show that the subcritical case of the Hopf bifurcation occurs robustly, which indicates the possibility of bistability. We also show how these oscillations lead to spatial wave formation as can be observed in real-world traffic flows.

Keywords:

1. Introduction

The so-called uniform flow equilibrium of vehicles following each other on a road is a kind of steady state, where equidistant vehicles travel with the same constant velocity. Ideally, this state is stable. Indeed, it is the goal of traffic management that drivers choose their traffic parameters to keep this state stable and also to reach their goal, that is, to travel with a speed close to their desired speed. Still, traffic jams often appear as congestion waves travelling opposite to the flow of vehicles (Kerner 1999). The formation of these traffic jams (waves) is often associated with the linear instability of the uniform flow equilibrium, which should be a rare occurrence. However, it is also well known among traffic engineers that certain events, such as a truck pulling out of its lane, may trigger traffic jams even when the uniform flow is stable. We investigate a delayed car-following model and provide a thorough examination of the subcriticality of Hopf bifurcations related to the drivers' reaction-time delay. This explains how a stable uniform flow can coexist with a stable traffic wave.

The car-following model analysed in this paper was first introduced in Bando et al. (1995) without the reaction-time delay of drivers, and that was investigated by numerical simulation. Then, numerical continuation techniques were used in Gasser et al. (2004) and Berg & Wilson (2005) by applying the package Auto (Doedel et al. 1997). Recently, Hopf calculations have been carried out in Gasser et al. (2004) for arbitrary numbers of cars following each other on a ring.

The reaction-time delay of drivers was first introduced in Bando et al. (1998), and its importance was then emphasized by the study of Davis (2003). In these papers, numerical simulation was used to explore the nonlinear dynamics of the system. The first systematic global bifurcation analysis of the delayed model (Davis 2003) was presented in Orosz et al. (2004), where numerical continuation techniques, namely the package Dde-Biftool (Engelborghs et al. 2001), were used. The continuation results were extended to large numbers of cars in Orosz et al. (2005), where the dynamics of oscillations, belonging to different traffic patterns, were analysed as well. In this paper, we perform an analytical Hopf bifurcation calculation and determine the criticality of the bifurcation as a function of parameters for arbitrary numbers of vehicles in the presence of the drivers' reaction delay.

While the models without delay are described by ordinary differential equations (ODEs), presenting the dynamics in finite-dimensional phase spaces, the appearance of the delay leads to delay differential equations (DDEs) and to infinite-dimensional phase spaces. The finite-dimensional bifurcation theory that is available in basic textbooks (Guckenheimer & Holmes 1997; Kuznetsov 1998) have been extended to DDEs in Hale & Verduyn Lunel (1993), Diekmann et al. (1995), Kolmanovskii & Myshkis (1999) and Hale et al. (2002). The infinite-dimensional dynamics make the bifurcation analysis more abstract. In particular, the Hopf normal form calculations require complicated algebraic formalism and algorithms, as is shown in Hassard et al. (1981), Stépán (1986, 1989), Campbell & Bélair (1995), Orosz (2004) and Stone & Campbell (2004). Recently, these Hopf calculations have been extended for systems with translational symmetry in Orosz & Stépán (2004), which is an essential property of car-following models. This situation is similar to the S1 symmetry that occurs in laser systems with delay, see Verduyn Lunel & Krauskopf (2000) and Rottschäfer & Krauskopf (2004). In Orosz & Stépán (2004), the method was demonstrated in the over-simplified case of two cars on a ring. Here, these calculations are extended to arbitrarily many cars, providing general conclusions for the subcriticality of the bifurcations and its consequences for flow patterns. Our results are generalization of those in Gasser et al. (2004) for the case without reaction-time delay. In particular, we prove that this delay makes the subcriticality of Hopf bifurcations robust.

2. Modelling issues

The mathematical form of the car-following model in question was introduced and non-dimensionalized in Orosz et al. (2004). Here, we recall the basic features of this model. Periodic boundary conditions are considered, i.e. n vehicles are distributed along a circular road of overall length L; see figure 1. (This could be interpreted as traffic on a circular road around a large city, e.g. the M25 around London, even though such roads usually possess higher complexity.) As the number of cars is increased, the significance of the periodic boundary conditions usually tends to become smaller.

Figure 1

Vehicles flowing clockwise on a circular road.

We assume that drivers have identical characteristics. Considering that the ith vehicle follows the (i+1)th vehicle and the nth car follows the 1st one, the equations of motion can be expressed asEmbedded Image(2.1)where the dot stands for time derivative. The position, the velocity, and the acceleration of the ith car are denoted by xi, Embedded Image and Embedded Image, respectively. The so-called optimal velocity function Embedded Image depends on the distance of the cars Embedded Image, which is usually called the headway. The argument of the headway contains the reaction-time delay of drivers which now is rescaled to 1. The parameter α>0 is known as the sensitivity and 1/α>0 is often called the relaxation time. Due to the rescaling of the time with respect to the delay, the delay parameter is hidden in the sensitivity α and the magnitude of the function V(h); see details in Orosz et al. (2004).

Equation (2.1) expresses that each driver approaches an optimal velocity, given by V(h), with a characteristic relaxation time of 1/α, but reacts to its headway via a reaction-time delay 1. The general features of the optimal velocity function V(h) can be summarized as follows:

  1. V(h) is continuous, non-negative and monotone increasing, since drivers want to travel faster as their headway increases. Note that in the vicinity of the Hopf bifurcation points, V(h) is required to be three times differentiable for the application of the Hopf theorem (Guckenheimer & Holmes 1997; Kuznetsov 1998).

  2. Embedded Image as h→∞, where Embedded Image is known as the desired speed, which corresponds, for example, to the speed limit of the given highway. Drivers approach this speed with the relaxation time of 1/α when the traffic is sparse.

  3. V(h)≡0 for h∈[0,1], where the so-called jam headway is rescaled to 1, see Orosz et al. (2004). This means that drivers attempt to come to a full stop if their headways become less than the jam headway.

One might, for example, consider the optimal velocity function to take the formEmbedded Image(2.2)as already used in Orosz et al. (2004, 2005). This function is shown together with its derivatives in figure 2. Functions with shapes similar to equation (2.2) were used in Bando et al. (1995, 1998) and Davis (2003).

Figure 2

(a) The optimal velocity function (2.2) and (b–d) its derivatives.

Note that the analytical calculations presented here are valid for any optimal velocity function V(h): it is not necessary to restrict ourself to a concrete function in contrast to the numerical simulations in Bando et al. (1998) and Davis (2003) and even to the numerical continuations in Orosz et al. (2004, 2005).

The dimensionless parameters α and v0 can be obtained from their dimensional counterparts as shown in Orosz et al. (2004). If we consider the reaction-time delay 1–2 s and the relaxation time 1–50 s, we obtain the sensitivity as their ratio: α∈[0.02,2]. Also, if we assume desired speeds in the range 10–40 m s−1 and jam headways 5–10 m, their ratio times the reaction time provides the dimensionless desired speed v0∈[1,20]. The linear stability diagram does not change qualitatively when v0 is varied in this range as shown in Orosz et al. (2005).

3. Translational symmetry and Hopf bifurcations

The stationary motion of the vehicles, the so-called uniform flow equilibrium, is described byEmbedded Image(3.1)whereEmbedded Image(3.2)andEmbedded Image(3.3)Note that one of the constants Embedded Image can be chosen arbitrarily due to the translational symmetry along the ring. Henceforward, we consider the average headway h*=L/n as a bifurcation parameter. Increasing h* increases the length L of the ring, which involves scaling all headways hi accordingly.

Let us define the perturbation of the uniform flow equilibrium byEmbedded Image(3.4)Using Taylor series expansion of the optimal velocity function V(h) about Embedded Image up to third order of Embedded Image, we can eliminate the zero-order termsEmbedded Image(3.5)whereEmbedded Image(3.6)At a critical/bifurcation point Embedded Image, the derivatives take the values Embedded Image, Embedded Image and Embedded Image. Now, and further on, prime denotes differentiation with respect to the headway.

Introducing the notationEmbedded Image(3.7)equation (3.5) can be rewritten asEmbedded Image(3.8)where Embedded Image. The matrices Embedded Image and the near-zero analytic function Embedded Image are defined asEmbedded Image(3.9)Here, Embedded Image stands for the n-dimensional identity matrix, while the matrix Embedded Image and the functions Embedded Image are defined asEmbedded Image(3.10)

System (3.8) possesses a translational symmetry, therefore, the matrices Embedded Image satisfyEmbedded Image(3.11)that is, the Jacobian Embedded Image has a zero eigenvalueEmbedded Image(3.12)for any value of parameter h*. Furthermore, the near-zero analytic function Embedded Image preserves this translational symmetry, that is,Embedded Image(3.13)for all c≠0 satisfying Embedded Image; see details in Orosz & Stépán (2004).

The steady state y(t)≡0 of equation (3.8) corresponds to the uniform flow equilibrium (3.1) of the original system (2.1). Considering the linear part of (3.8) and using the trial solution Embedded Image with Embedded Image and Embedded Image, the characteristic equation becomesEmbedded Image(3.14)According to equation (3.11), the relevant zero eigenvalue (3.12) is one of the infinitely many characteristic exponents that satisfy equation (3.14). This exponent exists for any value of the parameter b1, that is, for any value of the bifurcation parameter h*.

At a bifurcation point defined by Embedded Image, i.e. by Embedded Image, Hopf bifurcations may occur in the complementary part of the phase space spanned by the eigenspace of the zero exponent (3.12). Then, there exists a complex conjugate pair of pure imaginary characteristic exponentsEmbedded Image(3.15)which satisfies equation (3.14). To find the Hopf boundaries in the parameter space, we substitute Embedded Image into equation (3.14). Separation of the real and imaginary parts givesEmbedded Image(3.16)The so-called wave number k can take the values Embedded Image (for even n) and Embedded Image (for odd n). The wave numbers Embedded Image are not considered since they belong to conjugated waves producing the same spatial pattern. Note that k=1 belongs to the stability boundary of the uniform flow equilibrium, while the cases k>1 result in further oscillation modes around the already unstable equilibrium. Furthermore, for each k the resulting frequency is bounded so that Embedded Image; see Orosz et al. (2004, 2005).

Note also that the function Embedded Image shown in figure 2b is non-monotonous, and so a b1cr boundary typically leads either to two or to zero Embedded Image boundaries. For a fixed wave number k, these boundaries tend to finite values of h* as Embedded Image, as is shown in Orosz et al. (2005). Using trigonometric identities, equation (3.16) can be transformed toEmbedded Image(3.17)which is a useful form used later in the Hopf calculation together with the resulting formEmbedded Image(3.18)

With the help of the identityEmbedded Image(3.19)we can calculate the following necessary condition for Hopf bifurcation as the parameter b1 is varied asEmbedded Image(3.20)whereEmbedded Image(3.21)Since equation (3.20) is always positive, this Hopf condition is always satisfied. Now, using the chain rule and definition (3.6), condition (3.20) can be calculated further as the average headway h* is varied to giveEmbedded Image(3.22)This condition is fulfilled if and only if Embedded Image, which is usually satisfied except at some special points. For example, the function Embedded Image shown in figure 2c becomes zero at a single point over the interval Embedded Image. Notice that Embedded Image is zero for Embedded Image and for Embedded Image, but the critical headway hcr never takes these values for α>0. Conditions (3.20) and (3.22) ensure that the characteristic roots (3.15) cross the imaginary axis with a non-zero speed when the parameters b1 and h* are varied.

4. Operator differential equations and related eigenvectors

The DDE (3.8) can be rewritten in the form of an operator differential equation (OpDE) which reflects its infinite-dimensional dynamics. For the critical bifurcation parameter Embedded Image, we obtainEmbedded Image(4.1)where the dot still refers to differentiation with respect to the time t and Embedded Image is defined by the shift Embedded Image, Embedded Image on the function space Embedded Image of continuous functions mapping Embedded Image. The linear and nonlinear operators Embedded Image, Embedded Image are defined asEmbedded Image(4.2)Embedded Image(4.3)where the matrices Embedded Image and the nonlinear function Embedded Image are given byEmbedded Image(4.4)

The translational symmetry conditions (3.11) and (3.13) are inherited, that is,Embedded Image(4.5)andEmbedded Image(4.6)for all c≠0 satisfying (L+R)c=0.

In order to avoid singularities in Hopf calculations, we follow the methodology and algorithm of Orosz & Stépán (2004) and eliminate the eigendirection belonging to the relevant zero eigenvalue λ0=0 (equation (3.12)). The corresponding eigenvector Embedded Image satisfiesEmbedded Image(4.7)which, applying the definition (4.2) of operator Embedded Image, leads to a boundary value problem with the constant solutionEmbedded Image(4.8)One finds thatEmbedded Image(4.9)where each component of the vector Embedded Image is equal to 1. Here, Embedded Image is a scalar that can be chosen freely, in particular, we chooseEmbedded Image(4.10)

In order to project the system to s0 and to its complementary space, we also need the adjoint operatorEmbedded Image(4.11)where asterisk denotes either adjoint operator or transposed conjugate vector and matrix. The eigenvector Embedded Image of Embedded Image associated with the eigenvalue Embedded Image satisfiesEmbedded Image(4.12)This gives another boundary value problem, which has the constant solutionEmbedded Image(4.13)Here, we obtainEmbedded Image(4.14)However, Embedded Image is not free, but is determined by the normality conditionEmbedded Image(4.15)Defining the inner productEmbedded Image(4.16)condition (4.15) gives the scalar equationEmbedded Image(4.17)from which we obtainEmbedded Image(4.18)

Note that, as the vectors s0 and n0 are the right and left eigenvectors of the operator Embedded Image belonging to the eigenvalues Embedded Image and Embedded Image, similarly, the vectors S0 and N0 are the right and left eigenvectors of the matrix Embedded Image, belonging to the same zero eigenvalues. For the non-delayed model (Bando et al. 1995), the vectors S0 and N0 are the left and right eigenvectors of the Jacobian (L+R), that is, their structure is related the spatial structure of the system along the circular road.

We are now able to eliminate the singular dynamics generated by the translational symmetry so that we project the system to s0 and to its complementary space. With the help of the eigenvectors s0 and n0, the new state variables Embedded Image and Embedded Image are defined asEmbedded Image(4.19)Using the derivation given in Orosz & Stépán (2004), we split the OpDE (4.1) intoEmbedded Image(4.20)Its second part is already fully decoupled, and can be redefined asEmbedded Image(4.21)where the new nonlinear operator Embedded Image assumes the formEmbedded Image(4.22)Considering Embedded Image given by equation (4.3), and using the expressions (3.9), (3.10) and (4.4), and the eigenvectors (4.9) and (4.14), we obtainEmbedded Image(4.23)where the definition Embedded Image is applied. Note that the system reduction related to the translational symmetry changes the nonlinear operator of the system while the linear operator remains the same. This change has an essential role in the centre-manifold reduction presented below.

Let us consider a Hopf bifurcation at a critical point Embedded Image. First, we determine the real and imaginary parts Embedded Image of the eigenvector of the linear operator Embedded Image, which belongs to the critical eigenvalue Embedded Image (3.15), that is,Embedded Image(4.24)After the substitution of definition (4.2) of Embedded Image, the solution of the resulting boundary value problem can be written asEmbedded Image(4.25)with constant vectors Embedded Image satisfying the homogeneous equationEmbedded Image(4.26)

Using formula (3.17) and following the steps (A 1)–(A 3) of appendix A, one may solve (4.26) and obtainEmbedded Image(4.27)where the scalar parameters u and υ can be chosen arbitrarily and the vectors C, Embedded Image areEmbedded Image(4.28)with the wave number k used in equations (3.16) and (3.17). The cyclic permutation of the components in C and S results in further vectors S1 and S2 that still satisfy equation (4.26). This result corresponds to the Embedded Image symmetry of the system, that is, all the cars have the same dynamic characteristics. Choosing u=1 and υ=0 yieldsEmbedded Image(4.29)

The real and imaginary parts Embedded Image of the eigenvector of the adjoint operator Embedded Image associated with Embedded Image are determined byEmbedded Image(4.30)

The use of definition (4.11) of Embedded Image leads to another boundary value problem having the solutionEmbedded Image(4.31)where the constant vectors Embedded Image satisfyEmbedded Image(4.32)

Applying equation (3.17) and following the steps (A 4)–(A 7) of appendix A, the solution of equation (4.32) is obtained asEmbedded Image(4.33)The scalar parameters Embedded Image, Embedded Image are determined by the orthonormality conditionsEmbedded Image(4.34)which using the inner product definition (4.16), results in the two scalar equationsEmbedded Image(4.35)Substituting equations (3.17), (4.29) and (4.33) into (4.35) and using the second-order trigonometric identities (B 3)–(B 5) of appendix B, we obtainEmbedded Image(4.36)with the solutionEmbedded Image(4.37)where Embedded Image is defined by equation (3.21). The substitution of equation (4.37) into (4.33) leads toEmbedded Image(4.38)

As Embedded Image and Embedded Image are the right and left eigenvectors of the operator Embedded Image belonging to the eigenvalues Embedded Image and Embedded Image, the vectors Embedded Image and Embedded Image are similarly the right and left eigenvectors of the matrix Embedded Image belonging to the same eigenvalues. Note that for the non-delayed model (Bando et al. 1995) the vectors Embedded Image and Embedded Image are the left and right eigenvectors of the Jacobian Embedded Image, that is, their structure is again related the spatial structure of the system.

5. Centre-manifold reduction

Now, we investigate the essential dynamics on the two-dimensional centre manifold embedded in the infinite-dimensional phase space. Since the eigenvectors s1 and s2 span a plane tangent to the centre manifold at the origin, we use s1, s2 and n1, n2 to introduce the new state variablesEmbedded Image(5.1)where Embedded Image and Embedded Image. Using the derivation presented in Orosz & Stépán (2004), we can reduce OpDE (4.21) to the formEmbedded Image(5.2)The scalar parameters Embedded Image are induced by the translational symmetry and their expressions are determined in Orosz & Stépán (2004) asEmbedded Image(5.3)Since in our case Embedded Image, we obtainEmbedded Image(5.4)

The power series form of (5.2) is also given in Orosz & Stépán (2004) asEmbedded Image(5.5)where the subscripts of the constant coefficients Embedded Image and the vector ones Embedded Image refer to the corresponding jth and kth orders of z1 and z2, respectively. The terms with the coefficients Embedded Image and Embedded Image come from the linear combinations of s1(ϑ) and s2(ϑ). The translational symmetry only enters through the terms with coefficients Embedded Image, so that the terms with coefficients Embedded Image and Embedded Image refer to the structure of the modified nonlinear operator Embedded Image (equation (4.22)). Using the third- and fourth-order trigonometric identities (B 6)–(B 11) of appendix B, we can calculate these coefficients for wave numbers kn/2, kn/3 and kn/4 as given by equations (A 8) and (A 9) in appendix A.

Note that the cases Embedded Image, Embedded Image and Embedded Image result in different formulae for the above coefficients, but the final Poincaré–Lyapunov constant will have the same formula as in the case of general wave number k. The detailed calculation of these ‘resonant’ cases is not presented here.

Approximate the centre manifold locally as a truncated power series of w depending on the coordinates z1 and z2 asEmbedded Image(5.6)There are no linear terms since the plane spanned by the eigenvectors s1 and s2 is tangent to the centre manifold at the origin. Third and higher order terms are dropped. The unknown coefficients Embedded Image can be determined by calculating the derivative of w and substituting that into the third equation of (5.5). The solution of the resulting linear boundary value problem, given in detail in Orosz & Stépán (2004), isEmbedded Image(5.7)where the vectors Embedded Image satisfyEmbedded Image(5.8)Here, we also used thatEmbedded Image(5.9)in accordance with equation (A 8).

One can find that the 2n-dimensional equation for H0 is decoupled from the 4n-dimensional equation for H1, H2 in equation (5.8). Since (L+R) is singular due to the translational symmetry (4.5), the non-homogeneous equation for H0 in (5.8) may seem not to be solvable. However, its right-hand side belongs to the image space of the coefficient matrix (L+R) due to the translational symmetry induced terms Embedded Image. We obtain the solutionEmbedded Image(5.10)with the undetermined parameter κ that has no role in the following calculations.

At the same time, the non-homogeneous equation for H1, H2 in equation (5.8) is not effected by the vectors Embedded Image. Using equation (3.17) and following the steps (A 10)–(A 14) of appendix A one can find the solutionEmbedded Image(5.11)whereEmbedded Image(5.12)

Now, using equation (5.11) in (5.6) and (5.7), we can calculate Rw(−1) which appears in the first two equations of (5.5). In this way, we obtain the flow restricted onto the two-dimensional centre manifold described by the ODEsEmbedded Image(5.13)where the coefficients Embedded Image have already been determined by equation (A 8) in (5.5), while the coefficientsEmbedded Image(5.14)originate in the terms involving Rw(−1). To determine equation (5.14), the trigonometric identities (B 3)–(B 11) of appendix B have been used. The coefficients Embedded Image are not displayed, since they have no role in the following calculations.

We note that the coefficients Embedded Image of the second-order terms are not changed by the centre-manifold reduction. The so-called Poincaré–Lyapunov constant in the Poincaré normal form of equation (5.13) can be determined by the Bautin formula (Stépán 1989)Embedded Image(5.15)The bifurcation is supercritical for negative and subcritical for positive values of Δ. We found that Δ>0 is always true when Embedded Image which is the case for real traffic situations (many vehicles n with a few waves k). This can be proven as is detailed below.

The first part of the expression (5.15) of Δ in front of the parenthesis is always positive, since Embedded Image. Within the parenthesis in equation (5.15), the first term is positive since equation (3.16) implies Embedded Image, which yields critical headway values Embedded Image such that Embedded Image; see figure 2b,d. The second term in the parenthesis in equation (5.15) contains the ratio of two complicated expressions, which by using equations (3.16) and (5.12), can be rearranged in the formEmbedded Image(5.16)andEmbedded Image(5.17)Since equation (5.17) is always positive, the sign of equation (5.16) is crucial for deciding the overall sign of Δ. According to equation (3.16) Embedded Image, that is, the realistic case Embedded Image implies the oscillation frequency Embedded Image.

Figure 3a shows the numerator Embedded Image for some particular values of α demonstrating that Embedded Image for Embedded Image. Note that if α→0 then Embedded Image may become negative (see figure 3a for Embedded Image), but this is a physically unrealistic case where drivers intend to reach their desired speed Embedded Image extremely slowly.

Figure 3

Quantities defined by (5.16) and (5.17) as a function of the frequency ω, for representative values of parameter α. (a) Numerator Embedded Image and (b) ratio Embedded Image.

Moreover, the ratio of equations (5.16) and (5.17), Embedded Image, is not only positive for Embedded Image but also Embedded Image when ω→0 (i.e. when n→∞) as shown in figure 3b. This feature provides robustness for subcriticality. Note that subcriticality also occurs for optimal velocity functions different from equation (2.2), e.g. for those that are considered in Orosz et al. (2004).

Using definition (3.6), formulae (3.18) and (3.22) and expressions (5.15)–(5.17), the amplitude A of the unstable oscillations is obtained in the formEmbedded Image(5.18)Thus, the first Fourier term of the oscillation restricted onto the centre manifold isEmbedded Image(5.19)Since close to the critical bifurcation parameter Embedded Image, we have Embedded Image, equation (5.19) yieldsEmbedded Image(5.20)where the vectors S1, S2 are given in equation (4.29).

The parameter Embedded Image enters Δ and A via the derivatives Embedded Image, Embedded Image and Embedded Image which are all proportional to v0. Consequently, Δ depends linearly on v0 causing no sign change and v0 disappears from A. Furthermore, v0 is embedded in the critical parameter Embedded Image which is determined from equation (3.16) by inverting Embedded Image. However, one may check that this is relevant for Embedded Image only, when small v0 may result in supercriticality as demonstrated in Orosz et al. (2004). In contrast, the realistic case Embedded Image leads to robust subcriticality as explained in §6 and also demonstrated in Orosz et al. (2005).

Note that zero reaction time delay results in Embedded Image as shown in Gasser et al. (2004). In that case, subcriticality appears only for extremely high values of the desired speed v0 when the term 6b3cr becomes greater than Embedded Image at the critical points (of the non-delayed model). Consequently, the presence of the drivers' reaction-time delay has an essential role in the robustness of the subcritical nature of the Hopf bifurcation. This subcriticality explains how traffic waves can be formed when the uniform flow equilibrium is stable, as is detailed in the subsequent section.

6. Physical interpretation of results

The unstable periodic motion given in equation (5.20) corresponds to a spatial wave formation in the traffic flow, which is actually unstable. Substituting (4.29) into (5.20) and using definition (3.7), one can determine the velocity perturbation asEmbedded Image(6.1)The interpretation of this perturbation mode is a wave travelling opposite to the car flow with spatial wave number k (i.e. with spatial wavelength Embedded Image). The related wave speed isEmbedded Image(6.2)where the elimination of the frequency ω with the help of equation (3.18) leads toEmbedded Image(6.3)Since the uniform flow equilibrium (3.1) travels with speed Embedded Image, the speed of the arising wave isEmbedded Image(6.4)

By considering the optimal velocity function (2.2), we obtain cwave<0, that is, the resulting wave propagates in the opposite direction to the flow of vehicles. Note that the non-delayed model introduced in Bando et al. (1995) exhibits the same wave speed apart from some differences in the coefficient of the correction term Embedded Image. If one neglects this correction term, the wave speed becomes independent of n and k, which corresponds to the results obtained from continuum models; see e.g. Whitham (1999).

In order to check the reliability of the Poincaré–Lyapunov constant (5.15) and the amplitude estimation (5.18), we compare these analytical results with those obtained by numerical continuation techniques with the package Dde-Biftool (Engelborghs et al. 2001). In figure 4a, we demonstrate the subcriticality for n=9 cars and k=1 wave. The horizontal axis corresponds to the uniform flow equilibrium, that is, stable for small and large values of h* (shown by green solid line) but unstable for intermediate values of h* (shown by red dashed line) in accordance with formula (3.16) and figure 2b. The Hopf bifurcations, where the equilibrium loses its stability, are marked by blue stars. The branches of the arising unstable periodic motions given by equation (5.18) are shown as red dashed curves.

Figure 4

The amplitude A of velocity oscillations as a function of the average headway parameter h* (a) and the corresponding velocity profiles at h*=2.9 (b) and (c) for n=9 cars, k=1 wave and parameters α=1.0, v0=1.0. (a) Horizontal axis (A≡0) represents the uniform flow equilibrium and the analytical results are coloured: green solid and red dashed curves represent stable and unstable branches, respectively, and blue stars stand for Hopf bifurcations. Grey curves correspond to numerical continuation results: solid and dashed curves refer to stable and unstable states and grey crosses represent fold bifurcations. The points marked by Embedded Image and Embedded Image refer to the velocity profiles shown in panels (b) and (c), respectively. (b) Stable stop-and-go oscillations are shown for point P1 in green and for point Embedded Image in grey and (c) unstable oscillations are displayed for point P2 in red and for point Embedded Image in grey.

In §2, conditions (i)–(iii) require that the optimal velocity function is bounded so that Embedded Image; see figure 2a. Maximum principles show that for t→∞ the velocity of any solution is contained by the interval [0, v0]. This suggests that ‘outside’ the unstable oscillating state there exists an attractive oscillating state which may include stopping, accelerating, travelling with the desired speed v0, and decelerating. The amplitude of this solution is defined as Embedded Image. In figure 4a the horizontal green line at Embedded Image represents this stable stop-and-go oscillation. The corresponding stop-and-go wave propagates against the traffic flow, since vehicles leave the traffic jam at the front and enter it at the rear; see figure 1. The above heuristic construction reveals the existence of bistability on each side of the unstable equilibrium between the Hopf bifurcation point and the point where the branches of unstable and stable oscillations intersect each other (outside the frame in figure 4a). For such parameters, depending on the initial condition, the system either tends to the uniform flow equilibrium or to the stop-and-go wave.

In figure 4a, we also displayed the results of numerical continuation carried out with the package Dde-Biftool (Engelborghs et al. 2001). Grey solid curves represent stable oscillations while grey dashed curves represent unstable ones. The fold bifurcation points, where the branches of stable and unstable oscillations meet, are marked by grey crosses. The comparison of the results shows that the analytical approximation of the unstable oscillations is quantitatively reliable in the vicinity of the Hopf bifurcation points. The heuristic amplitude v0/2 of the stop-and-go oscillations is slightly larger than the numerically computed ones. The analytically suggested bistable region is larger than the computed one (between the Hopf point and the fold point), since the third degree approximation is not able to predict fold bifurcations of periodic solutions. In order to find these fold bifurcation points, it is necessary to use numerical continuation techniques as presented in Orosz et al. (2005). Nevertheless, qualitatively the same structure is obtained by the two different techniques.

As was already mentioned in §3, the wave numbers k>1 are related to Hopf bifurcations in the parameter region, where the uniform flow equilibrium is already unstable. This also means that the corresponding oscillations for k>1 are unstable independently of the criticality of these Hopf bifurcations. Still, we found that these Hopf bifurcations are all robustly subcritical for any wave number k (except for large Embedded Image). Consequently, the only stable oscillating state is the stop-and-go motion for k=1. On the other hand, several unstable solutions may coexist as is explained in Orosz et al. (2005). Note that analytical and numerical results agree better as the wave number k is increased because the oscillating solution becomes more harmonic.

To represent the features of vehicles' motions, the velocity oscillation profiles of the first vehicle are shown in figure 4b,c for the points Embedded Image and Embedded Image marked in figure 4a, for headway h*=2.9. Again, the coloured curves correspond to the analytical results, while the grey curves are obtained by numerical continuation. In figure 4b,c, the time window of each panel is chosen to be the period of the first Fourier approximation given by equation (3.18) (red curve in figure 4c) and the dashed vertical lines indicate oscillation periods computed numerically with Dde-Biftool.

Figure 4b shows the stop-and-go oscillations. The heuristic construction (green curve) is obtained by assuming that the stopping and flowing states are connected with states of constant acceleration/deceleration, which is qualitatively a good approximation of the numerical result (grey curve). Figure 4c compares the unstable periodic motions computed analytically (red curve) from the Hopf calculation with those from numerical continuation (grey curve).

For a perturbation ‘smaller’ than the unstable oscillation, the system approaches the uniform flow equilibrium. If a larger perturbation is applied then the system develops stop-and-go oscillations and a spatial stop-and-go travelling wave appears as demonstrated in figure 1. Since the period of the stable and unstable oscillations are close to each other, the stable stop-and-go wave travels approximately with the speed of the unstable travelling wave; see equation (6.4) for k=1.

7. Conclusion

A nonlinear car-following model has been investigated with special attention paid to the reaction-time delay of drivers. By considering the average headway as a bifurcation parameter, Hopf bifurcations were identified. In order to investigate the resulting periodic motions, the singularities related to the essential translational symmetry had to be eliminated. Then the Hopf bifurcations were found to be robustly subcritical leading to bistability between the uniform flow equilibrium and a stop-and-go wave. The appearing oscillations manifest themselves as spatial waves propagating backward along the circular road.

In the non-delayed model of Bando et al. (1995), subcriticality and bistability occur only for extremely high values of the desired speed v0, as it is demonstrated in Gasser et al. (2004). We proved that subcriticality and bistability are robust features of the system due to the drivers' reaction-time delay, even for moderate values of the desired speed. This delay, which is smaller than the macroscopic time-scales of traffic flow, plays an essential role in this complex system because it changes the qualitative nonlinear dynamics of traffic.

Due to the subcriticality, stop-and-go traffic jams can develop for large enough perturbations even when the desired uniform flow is linearly stable. These perturbations can be caused, for example, by a slower vehicle (such us a lorry) joining the inner lane flow for a short-time interval via changing lanes. It is essential to limit these unwanted events, for example, by introducing temporary regulations provided by overhead gantries. Still, if a backward travelling wave shows up without stoppings, it either dies out by itself or gets worse ending up as a persistent stop-and-go travelling wave. In order to dissolve this undesired situation, an appropriate control can be applied using temporary speed limits given by overhead gantries that can lead the traffic back ‘inside’ the unstable travelling wave and then to reach the desired uniform flow. For example, the MIDAS system (Lunt & Wilson 2003) installed on the M25 motorway around London is able to provide the necessary instructions for drivers.

Acknowledgments

The authors acknowledge with thanks discussions with and comments of Eddie Wilson and Bernd Krauskopf on traffic dynamics and on numerical bifurcation analysis. This research was supported by the University of Bristol under a Postgraduate Research Scholarship and by the Hungarian National Science Foundation under grant no. OTKA T043368.

Footnotes

  • Alternative address: Research Group on Dynamics of Vehicles and Machines, Hungarian Academy of Science, PO Box 91, Budapest 1521, Hungary.

    • Received September 1, 2005.
    • Accepted January 6, 2006.

References

View Abstract