The linear stability properties of car-following models of highway traffic are analysed. A general family of models is introduced and the subsequent analysis developed in terms of its partial derivatives. Two measures of wave propagation, namely (i) the group velocity and (ii) the signal velocity, are introduced and computed. These measures are used to classify how instability propagates disturbances, measured relative to the frame of the road along which the vehicles drive. Detector data suggest that disturbances should propagate only in an upstream direction (convective upstream instability), and it is shown how to parametrize models to agree with data and avoid unrealistic downstream propagation (absolute and convective downstream instability).
Since the mid-1990s, there has been an intense and sometimes argumentative discussion about the spatio-temporal patterns observed in highway traffic flow, see Kerner & Rehborn (1996, 1997), Kerner (2002a,b), Treiber et al. (2000, 2010), Helbing et al. (2009) and Schönhof & Helbing (2009). A typical spatio-temporal pattern is displayed in figure 1. The features of interest are a region of congestion developing at an on-ramp, out of which nucleate stop-and-go waves, which propagate upstream at about 15–20 km h−1—a figure that seems almost universal in that it has been observed on highways in many different countries across the Western world (Zielke et al. 2008). One explanation for stop-and-go waves is that uniform traffic flow is unstable (either linearly or nonlinearly) under certain conditions. The inevitable fluctuations (owing to lane-changing, noise, etc.) at the level of individual vehicles then become magnified so as to result in large-amplitude patterns with macroscopic spatial scales (i.e. much longer than the range of the individual vehicle interactions).
To understand the ‘up-scaling’ process, we analyse the stability properties of car-following models that describe vehicles as discrete entities moving in continuous space and time, whose motions solve simple ordinary (or delay) differential equations coupling the motion of nearest neighbours. In terms of the linear stability properties that we focus on here, this analysis began in the late 1950s (Chandler et al. 1958). In particular, Herman et al. (1959) developed the concept, now known as string instability, in which a growing wave envelope may propagate up a column of vehicles, even when the ‘local’ dynamics of each individual vehicle are stable. In mathematical language, string instability is a special case of convective spatio-temporal instability, in that growth is via a localized travelling envelope, in contrast to absolute instability where the linearized model experiences blow-up at all points in the spatial domain.
In realistic traffic models, driver behaviour is governed by vehicles in front, rather than behind, so instability is necessarily convective when measured in the frame of the vehicles. The complication is that data such as figure 1 are presented in the frame of the road, rather than in the frame of the vehicles, each of which is driving forward relative to the road. So an instability that is convective in the vehicle frame can be either absolute or convective in the road frame, depending on whether vehicles drive forward faster than the wave envelope propagates backwards relative to them (figure 2). The importance of this distinction between convective and absolute instability in the vehicle and road frames has been recognized by Treiber et al. (2000, 2010), but so far a full mathematical analysis has been lacking.
In this paper, we present rigorous mathematical criteria that distinguish between convective string instability and absolute string instability in the road frame. The calculations are developed for a general family of car-following models introduced in §2. The basic structure of the linearization is then explained in §3. The chief idea is to calculate ranges of wave speeds that correspond to the propagation of information. The calculations that then follow are of two types:
— Firstly (§4), we compute the group velocity that can be supported by growing modes. This work is thus a generalization of Mitarai & Nakanishi (2000), who computed group velocity for a particular car-following model known as the optimal velocity model. More recent work by Helbing & Johansson (2009) revisited the optimal velocity model and in addition computed the group velocity for a family of related macroscopic (partial differential equation) models. All of these group velocity calculations are attractive because they are simple both conceptually and in terms of their mathematical details. However, this approach is based on the flawed assumption that group velocity is equivalent to the velocity of information—when, in fact, this assumption is true only for non-dissipative media (Brillouin 1960).
— The correct approach (§5) employs the calculation of signal velocity, which may be defined loosely as the velocity at which a new signal can penetrate a medium that is at rest. This calculation employs the steepest descent method to determine the large-time asymptotics of the vehicle trajectories. This leads to an intriguing mathematical problem concerning the saddle selection procedure that is described in detail in appendix B.
Unfortunately (from the point of view of the practitioner), the extra complication of the signal velocity calculation is required, because the group velocity and the signal velocity are only the same for non-dissipative media (Brillouin 1960). For example, in light propagation, it is known that dissipative media may exhibit superluminal group velocity (Garrett & McCumber 1970; Wang et al. 2000).
Section 6 gives numerical results for two particular examples of car-following models. In particular, figure 1 suggests that in congested conditions, information should propagate only upstream in the road frame, and hence we seek models and parametrizations that replicate this behaviour. In fact, the absence of downstream growth in macroscopic data is not incompatible with convective downstream or absolute linear instability, if it is combined with more complicated model features. For example, one might try to construct set-ups in which the downstream growth is trapped at junctions by the spatial heterogeneity of the road, and where the consequent upstream propagation results from nonlinear effects.
However, in the parsimonious framework that we consider, it is a sensible starting point to require that models display only convective upstream linear instability in the road frame. As we shall see, this is a rather strong criterion in model parametrization, which not all models obey for all parameter values. This practical theme is taken up further in §7, where conclusions and opportunities for further work are presented.
2. Modelling details
We follow the general car-following model framework developed in Wilson (2008). Our starting point is the standard situation depicted in figure 3. We consider a single lane of traffic with identical vehicles labelled 1, 2, etc., in the upstream direction. Displacements and velocities are denoted xn(t) and vn(t)≥0, respectively, and our models shall also involve the front-to-front spacing hn(t):=xn−1(t)−xn(t)>0 of consecutive vehicles, sometimes referred to as the headway. Note that overtaking is neglected in our framework in return for analytical tractability, and our approach is to view lane changes (and other imperfections, such as noise, heterogeneous road, differences between drivers, etc.) as external perturbations to a deterministic single-lane model whose stability should then be analysed.
In their simplest form, car-following models consist of a set of coupled differential equations for the trajectory of each vehicle, which supplements the kinematic relations 2.1with a behavioural model 2.2Here, dot denotes differentiation with respect to time, and equation (2.2) mimics how drivers accelerate or decelerate in response to the distance to the vehicle in front, the relative velocity 2.3and their own velocity. The point of this paper is to derive results for models in the general form of equation (2.2) with the minimum number of additional assumptions. However, to demonstrate the theory in operation, we illustrate the paper with two representative examples:
— The optimal velocity with relative velocity (OVRV) model, see Jiang et al. (2001) and Ward (2008). We have 2.4where V is the so-called optimal velocity function defining the maximum safe speed for a given headway, which drivers relax towards at rate α>0. The term (β>0) models the tendency for drivers to brake when closing in on their predecessor, and to accelerate when the gap is increasing. The OVRV model is not intended to model quantitative details of driver behaviour, but rather to capture the correct qualitative features in the simplest possible functional form. When β=0, the OVRV model reduces to the much-studied optimal velocity model due to Bando et al. (1995). We will follow Bando et al. (1995) and work in dimensionless variables with the choice 2.5
— The intelligent driver model (IDM), see Treiber et al. (2000). We have 2.6where 2.7and the notation and standard parameter values are given in table 1. In contrast to the OVRV model, the IDM is an attempt to model driver behaviour quantitatively in dimensional terms, and we may view equations (2.6) and (2.7) as a proxy for the complicated schemes that one usually finds in commercial microsimulation packages.
Note in practice there are many possible refinements to equation (2.2), such as the inclusion of reaction-time delay, multi-anticipation (where the motion of more than one vehicle ahead is considered), lane-changing effects and heterogeneity of the vehicle fleet and driver population—but these are beyond the scope of this paper: see Brackstone & McDonald (1999) and Helbing (2001) for reviews.
3. Uniform flow and linear stability preliminaries
Motivated by data, we focus on models described by equation (2.2) that have an equilibrium speed-headway function V such that 3.1In some models, such as the OVRV (equation (2.3)), the function V is provided explicitly as a parameter, whereas in others, such as the IDM (equations (2.6) and (2.7)), it is derived by isolating v* from the relation f(h*,0,v*)=0. Whichever is the case, we require that V satisfies: (i) V (0)=0, (ii) V ′≥0, and (iii) as .
The consequence of equation (3.1) is a one-parameter family of steady driving solutions known as uniform flows, in which xn(t)=x0−nh*+tV (h*). Thus, in a uniform flow solution, all vehicles drive at the same speed V (h*) (so that the relative velocity of any pair of vehicles is zero) with the same time-independent headway h*. In the t–x plane, vehicle trajectories are thus equally spaced parallel straight lines.
Our approach is to consider small perturbations to the uniform flow equilibria by setting and , where and are small. Assuming f is sufficiently smooth, this linearization yields 3.2where the partial derivatives are evaluated at the constant equilibrium arguments (h*,0,V (h*)), and necessary constraints for rational driver behaviour are 3.3These conditions simply state that drivers increase their acceleration in response to an increase in headway or relative velocity but tend to decrease their acceleration as their own velocity increases. At very large headway, the sensitivity of the dynamics to the vehicle in front is very slight, so we should also allow . However, since the focus of our stability calculations is on close-following situations, we shall maintain strict inequality in equation (3.3).
Since equation (2.3) implies and , we may apply the difference operator to equation (3.2) to obtain 3.4which is the key equation in this paper. Fundamentally, equation (3.4) describes the second-order dynamics of vehicle n driven by those of vehicle n−1 ahead of it. We stress that in this modelling framework, perturbations can only propagate upstream relative to the vehicle frame; however, it remains an open question as to which direction disturbances propagate relative to the road frame.
Let us consider the left-hand side differential operator of equation (3.4). Its characteristic equation (cf. Kesting & Treiber 2008) is 3.5Since all coefficients are positive (by equation (3.3)), the roots (that we call platoon eigenvalues) are either complex conjugate with negative real parts, or both real and negative. Hence, the differential operator is stable. The consequence is that if an instantaneous perturbation is applied to a finite platoon of vehicles whose leader’s velocity is held fixed, then as , every vehicle returns to equilibrium. Hence, the models that we consider do not exhibit platoon instability (sometimes referred to as local instability in the classical traffic literature).
However, the key point is that a platoon-stable model can nevertheless display instability when observed in a moving frame of reference. This is the essential point of convective instability, also known as string instability in the classical traffic literature. To understand this point, consider a semi-infinite column of vehicles in uniform flow whose leader drives at a fixed speed. Let us suppose that the second vehicle in the column is instantaneously ‘kicked’ out of equilibrium, and let us study the ensuing fluctuations that propagate back up the column of vehicles. Each vehicle in turn will be disturbed momentarily from equilibrium, but then return to steady driving owing to platoon stability. However, if the column is string-unstable, the maximum departure from equilibrium that each vehicle experiences will increase as one goes further and further up the column of vehicles. See Wilson & Ward (2010) for a more detailed discussion of this process.
4. Group velocity calculations
Our first attempt at understanding string instability is a development of the general Fourier analysis presented in Wilson (2008). The point is to solve equation (3.4) via the ansatz 4.1which yields the quadratic characteristic equation 4.2to solve for the (generally complex) growth rate λ in terms of the discrete wavenumber θ, 0<θ≤π. Here, a small (positive) value of θ corresponds to very long wavelength fluctuations, and θ=π gives the shortest possible wavelength in this discrete setting, i.e. a perturbation that is period 2 in the vehicle index. The solutions to equation (4.2) are given exactly by 4.3where stability of a given mode is governed by the root λ+ with greatest real part. In Wilson (2008), it is proved that if any given wavenumber θ is unstable, i.e. if λr(θ):=Re(λ+(θ))>0, then all smaller wavenumbers (i.e. longer wavelengths) are unstable too—in the sense that for . Hence, stability of a given uniform flow solution is controlled by long wavelengths in that one only needs to test the sign of λr(0+).
From equation (4.2), we have and hence that the real and imaginary parts of λ are even and odd functions of θ, respectively. Thus, we seek the expansion 4.4with real coefficients λ1, λ2, …, etc. Here, the sign of the real part 4.5governs whether a given mode grows or decays, and we define 4.6so that ω(θ)/θ gives the phase velocity measured relative to the discrete lattice coordinate, with a positive value denoting upstream propagation.
— For λ2<0, the given uniform flow is linearly stable in that all wavenumbers are linearly stable.
— For λ2>0, the given uniform flow is linearly unstable to small wavenumber (long wavelength) modes, and we may speak of a window of instability, i.e. a range of unstable wavenumbers for which λr(θ)>0.
Plots of λr(θ) for the stable and unstable cases are displayed in figure 4.
When a uniform flow is only marginally unstable, the window of instability is confined to relatively small wavenumbers. In this case, λr(θ) is well approximated by the first two terms in its series, see equation (4.5). Since it may be shown that λ4<0 when λ2>0 (appendix A), it follows that 4.9However, the exact value of may also be calculated numerically from equation (4.3).
The point now is to determine the speed of packets of waves (focusing later on unstable modes for which λr(θ)>0). In the lattice frame (i.e. relative to the column of vehicles), the upstream speed of a packet of wavenumber θ is given by the well-known group velocity ω′(θ)=dω/dθ, which must be multiplied by the underlying lattice spacing h* to give a dimensional speed. However, since vehicles (i.e. the lattice points) are moving downstream with the underlying velocity V (h*), the velocity of the wave packet relative to the road frame is given by 4.10where a positive (resp. negative) value denotes downstream (resp. upstream) propagation. Henceforth, we refer to cg(θ) defined by equation (4.10) as simply the group velocity.
To calculate cg(θ) exactly, we need to differentiate the quadratic formula in equation (4.3) directly. However, for small θ, we may alternatively use equation (4.7) and the first two terms of equation (4.6) to obtain 4.11Since it may be shown that λ3>0 when λ2>0 (appendix A), it follows that cg(θ) is an increasing function for small θ. In fact, cg(θ) is increasing for all θ, but the proof requires the direct differentiation of equation (4.3) (details omitted here).
We now return to the fundamental question, namely how to classify the instability of uniform flow. The approach that we take here is to view the group velocity cg(θ) as a proxy for the velocity of propagation of information. In fact, as we have discussed, this approach is flawed, because the traffic models that we consider are not non-dissipative and so instead the signal velocity should be used. (This is a more difficult calculation that follows in §5.)
However, for the time being, we identify the range of group velocities exhibited by unstable wavenumbers, and consequently the range of velocities with which wave packets can be propagated without decay. Specifically, we take an unstable uniform flow (for which λ2>0) and define 4.12In some sense, c−g and c+g thus capture, respectively, the velocities of the upstream and downstream edges of the wedge of instability (figure 2). The three types of instability are then characterized as follows:
— c−g<c+g<0. Convective upstream instability (wave packets only propagate upstream).
— c−g<0<c+g. Absolute instability (wave packets propagate in both directions).
— 0<c−g<c+g. Convective downstream instability (wave packets propagate only downstream).
This classification is depicted in the numerical examples in figure 5. Because cg(θ) is an increasing function, it follows that 4.13Thus, we have the exact formula 4.14but the exact calculation of c+g requires the application of numerical methods to equation (4.3), to calculate θmax and the derivative ω′. However, for cases that are only marginally unstable, in that is small, we may use equation (4.9) in combination with equation (4.11) to obtain 4.15
The analytical details of the group velocity calculation are now complete. Later, in §6, we classify whole ranges of uniform flows for the exemplar models introduced in §2, and in particular we chart how the stability changes with parameters. This process involves calculating the transitions between absolute and convective upstream (resp. downstream) instability defined by the locus of c+g=0 (resp. c−g=0).
5. Signal velocity calculations
So far we have used the group velocity (which characterizes the propagation of wave packets) in order to classify instability in the linearized general car-following model defined in equation (2.2). We now use asymptotic methods to calculate the signal velocity, at which information penetrates a medium that is initially at rest (Brillouin 1960). As we have discussed, the group velocity and the signal velocity are usually different. The group velocity is found by a simpler calculation, but it is the signal velocity that correctly represents the propagation of information.
Our set-up is a notional experiment in which a semi-infinite column of vehicles n≥0 are given equilibrium initial data hn(0)=h* and vn(0)=V (h*). We then prescribe a small disturbance to the trajectory of the lead vehicle n=0, which is localized in time, and we examine how this perturbation drives the motion of the vehicles upstream. Thus, for n≥1, we take the Laplace transform of equation (3.4) and obtain 5.1where Fn(z) denotes the Laplace transform of . This recursion relation can be solved to yield 5.2where 5.3For n≥1, we apply the Laplace inversion formula 5.4where γ is chosen to the right of the poles of g(z), given by the platoon eigenvalues μ that solve equation (3.5). In fact, this integral can be solved explicitly by calculating the residue at the poles, but since these are degree n, we must apply the Leibnitz formula to compute the (n−1)th derivative of the regular part of the integrand. This procedure is tractable but yields a complicated double sum that cannot be analysed further.
Hence, instead we analyse equation (5.4) by asymptotic methods. Our approach is to apply the method of steepest descents to equation (5.4) to extract the behaviour: see Bender & Orszag (1999, ch. 6) for details. The general scope of the method is the integral 5.5where f(z) and ρ(z) are analytic except at isolated singularities (poles, branch cuts, etc.). The idea is then to consider an alternative contour 𝒞′ that is obtained by deforming 𝒞. If this can be achieved without crossing the integrand’s singularities, then I is unchanged by Cauchy’s theorem. We seek such a 𝒞′ that passes through a saddle point z* at which ρ′(z*)=0, and at which is attained the local and global maxima of Re(ρ(z)) along the contour. In consequence, as , the local properties of the integrand at the saddle point dominate the integral, and the leading-order asymptotics may be computed in the form 5.6where .
In its present form, the integral (5.4) gives ρ(z)=z, and consequently there is no saddle point. The trick is thus to write the integrand in the form 5.7so that we may proceed with f(z)=F0(z) and 5.8where κ=n/t. For the method of steepest descents to work, we require ρ(z) and f(z) to be independent of t, so we fix κ while we apply the asymptotics. In essence, we thus analyse wave propagation along a ray with speed κ upstream in the lattice frame.
To clarify the structure of the calculation: we have , where I is given asymptotically (for large t) by equation (5.6). In equation (5.6), the saddle point z* is a function of κ=n/t, which is found by the procedure that we now develop.
The saddle points are found by computing 5.9and thus we seek solutions of 5.10yielding the cubic equation 5.11parametrized by κ.
In general, equation (5.11) must be solved numerically. The question is then for a given κ, which of the three roots of equation (5.11) should be taken as z* in formula (5.6). The requirements are that (i) the saddle point in question is a local maximum (rather than a local minimum) of Re(ρ(z)) in the direction of the contour, (ii) that it is also the global maximum of Re(ρ(z)) on the contour, and (iii) that the contour may be obtained by deformation of the line without crossing the poles of g(z) (i.e. the platoon eigenvalues μ). These requirements concern the level-set geometry of Re(ρ(z)) that is considered in detail in appendix B.
The selection procedure developed in appendix B has some delicate features for small κ that nevertheless may be automated by computer. However, for , the following particularly simple prescription holds:
Case 1. Equation (5.11) has one real and two complex conjugate roots. The two complex conjugate roots are the saddles z* of interest and the contour should be deformed over both, to incorporate two terms of the form of equation (5.6). However, each complex conjugate contributes the same exponential growth factor, so that either may be used in the criterion that follows.
Case 2. Equation (5.11) has three real roots. The greatest, i.e. is the one of interest.
Once the correct saddle z*(κ) has been selected, formula (5.6) may be used to derive the asymptotics of as we have described above. Note that this procedure gives rise automatically to a real answer for . In case 2, this is because the saddle in question is real, whereas in case 1, it is due to cancellation of the conjugate components from the pair of saddles that is used.
Figure 6a shows an example of just how close the asymptotic agreement is. However, more importantly, we may extract the dominating exponential growth factor from equation (5.6), in the form 5.12where 5.13If Φ(κ)>0, then perturbations grow along the ray in question, whereas if Φ(κ)<0, they decay.
The next task is to consider how Φ depends on κ. We have no proofs as to the general structure, but rather we use a numerical procedure that loops over κ, solving the cubic equation (5.11) and selecting the correct saddle for each individual value of κ. A typical computation is illustrated in figure 6b. In general, for linearly unstable uniform flows, we find κ1, κ2>0 such that Φ(κ)>0 for κ1<κ<κ2, and Φ(κ)<0 elsewhere. Thus, κ1 and κ2 identify critical speeds bounding a wedge in which perturbations grow. Thus, we may establish signal velocities 5.14measured relative to the road frame, which bound the envelope in which perturbations grow. The signs of the signal velocities c−s and c+s may then be analysed in the same way as those of the group velocities c−g and c+g (§4) in order to classify the type of stability.
6. Examples and illustrations
So far, in §§4 and 5, we have developed general analytical methods in order to establish a recipe that can be used to classify the stability of any car-following model that complies with the general formulation (2.2). To illustrate what may be achieved, we now apply these techniques to the exemplar models introduced in §2.
Note that it is no longer sufficient to study the stability of a single uniform flow solution. Rather, a car-following model should be characterized by analysing the stability class of all of its uniform flow solutions, as the equilibrium headway (or velocity) and the model parameters are varied.
To begin with, we fix the model parameters by choosing (i) α=0.6 and β=0.2 for the OVRV model and (ii) by selecting the standard published parameters for the IDM model (table 1). We then scan through a range of equilibrium velocities v* and compute the corresponding bounding group velocities c±g (§4) and signal velocities c±s (§5).
See figure 7a,b for the results. In each plot, we observe ‘bubbles’ formed by the pairs of curves (c−g,c+g) and (c−s,c+s), which mark the extent of the unstable regime. Note that in general the group velocities and signal velocities do not agree, except at the left- and right-hand edges of the bubbles that correspond to the onset of instability—at which points the underlying linear models are non-dissipative.
Furthermore, in figure 7a,b, we classify ranges of equilibrium velocity v* according to the bounding signal velocities, so that
— Cu (convective upstream instability) corresponds to ;
— A (absolute instability) corresponds to ; and
— Cd (convective downstream instability) corresponds to .
In the unshaded ranges labelled S, the dynamics are linearly stable because λ2<0, see equation (4.8).
When traffic is unstable, we find 6.1and we conjecture that this is a general result. In consequence, the range (c−s,c+s) of possible signal velocities overlaps but is downstream of the range (c−g,c+g) of possible group velocities. Extreme caution is thus required if group velocity is used as a proxy for the propagation of information. In particular, a model may be considered ‘good’ from the point of view of group velocity, in that and thus downstream propagation in the road frame appears to be avoided. However, one may have c+s>0 so that instability is in fact absolute.
To investigate the classification of instability further, we now compute two-parameter diagrams in which both the equilibrium velocity and a single-model parameter are varied. For the latter (and for the purposes of illustration), we choose the sensitivity parameters, namely α in the OVRV model (equations (2.4) and (2.5)) and a in the IDM model (equations (2.6) and (2.7)).
See figure 8 for the results, where the regimes are computed and labelled according to the signal velocity, in the same way as for the one-parameter pictures in figure 7. The boundaries between the different instability regions (denoted by solid lines) are thus computed by tracing out the loci of c−s=0 and . In addition, dashed lines show the loci of c−g=0 and c+g=0, so that it is clear that using the group velocity will not achieve the same classification as using the signal velocity. The boundary of the stable regime S is found by computing the locus of λ2=0.
Note that figure 8 establishes that the stability charts for the OVRV model and the IDM are qualitatively the same. We have found that most ‘reasonable’ car-following models seem to generate similar pictures, but this assertion can only be supported by a much more detailed numerical study that is beyond the scope of this paper.
As we have discussed, we consider it appropriate that practical car-following models are parametrized so that linear stability is only ever convective upstream (type Cu) in the road frame. However, in the IDM model, the standard choice a=0.73 m s−2, indicated by the horizontal dotted line in figure 8b, attains all three types of instability. Thus, our analysis suggests that a larger value, e.g. a≃1.2 m s−2, is better.
In fact, in figure 8, we may observe that as the sensitivity parameter (α or a) is reduced, the onset of instability is convective upstream, i.e. the desirable kind of instability. Only when the sensitivity parameter is reduced much further do we also obtain (undesirable) absolute and convective downstream instability. This suggests that when parametrizing microsimulation models, one should generally choose parameters that are only slightly unstable at the linear level, so that only type Cu instability is obtained.
Why is the onset of instability (as α or a is reduced) convective? Recall (figure 7) that at onset the bounding signal and group velocities coincide so that the group velocity classifies the type of instability correctly. Further, at onset, , so that c+g=cg(0+)=c−g from equation (4.13). Thus, there is a single value of group velocity (rather than an admissible window), and instability is thus convective upstream (resp. downstream) if the group velocity is negative (resp. positive). To distinguish the two possible types, equation (4.14) gives 6.2at onset, where ρ*=1/h*, and 6.3with . Relation (6.3) is the so-called fundamental diagram relating macroscopic flow and density. Equation (6.2) then captures the well-known property of hydrodynamic theory that the local wave-propagation speed is equal to the gradient of the fundamental diagram.
To guarantee that the onset of instability is convective upstream, it thus suffices to parametrize one’s car-following model so that the first instability (as sensitivity is reduced) occurs to the right of the fundamental diagram’s maximum value, so that Q′(ρ*)<0.
Finally, note that the locus of c−g=0 corresponds in general to the maximum of the fundamental diagram at which Q′(ρ*)=0—a condition that is independent of the sensitivity (α or a), thus giving rise to the vertical sections of the dashed line in figure 8a,b.
The message of this paper is that there are important distinctions in the type of instability that highway traffic models display. The concepts presented here are standard in the mathematical wave propagation literature (Huerre & Rossi 1998, §3), but are not appreciated by traffic modellers and in particular the microsimulation community. However, we view the stability properties as an essential part of the design and parametrization of car-following models. Firstly, no good simulation model should display platoon instability (equation (3.5)). Secondly, if a model displays string instability, then one should consider how it propagates disturbances relative to the frame of the road. Our view is that only convective upstream (type Cu) instability should be allowed. This is because spatio-temporal patterns in detector data that display growth only do so in the upstream direction. However, we should emphasize that the absence of downstream growth in data does not rule out downstream linear instability (of types Cd and A), if it is combined with more complicated model features (e.g. nonlinearity, or spatial heterogeneity at junctions). Thus, the arguments presented here should be viewed only as the first step in a more detailed modelling debate.
We have seen that there are two quantities that can be used to describe the overall wave propagation properties and classify the direction of instability. These are the group velocity (§4) and the signal velocity (§5). The group velocity describes the propagation of wave packets, whereas the signal velocity describes how a new disturbance penetrates a medium that is initially at rest. Our view is that traffic instability is triggered by instantaneous microscale events in the vicinity of junctions, and so the signal velocity is the correct measure to determine the envelope of the resulting disturbance and so classify the type of instability.
Unfortunately, the calculation of the signal velocity is more complicated than that for the group velocity, but in principle, all of the steps shown here can be automated using automatic differentiation and root-finding methods, if a car-following model were to be supplied as a black-box routine that complies with the structure described in §2. However, more work is required to generalize the model structure still further, for example, to incorporate reaction-time delay and multi-anticipation effects.
Finally, an example of the kind of practical lesson that can be achieved by our analysis is displayed in figure 8—which although it is model-specific, suggests more broadly that the undesirable types of string instability (convective downstream Cd and absolute A) may be eliminated if models are parametrized so that linear instability is only marginal. Recent analysis of empirical stop-and-go wave data (Zielke et al. 2008) has indicated surprisingly small growth rates that are consistent with the suggestion that linear instability is a marginal effect. However, an alternative perspective (Wilson 2008) for further investigation is that nonlinear stability properties are also crucial in determining the overall wave propagation properties.
We thank Martin Treiber who brought this problem area to our attention at a meeting in Hamburg in October 2007, and for interesting conversations with him since, in Vienna (May 2008) and in Dresden (April 2009). J.A.W. acknowledges the support of the Science Foundation Ireland (grant ref. 06/MI/005). R.E.W. acknowledges the support of EPSRC Advanced Research Fellowship EP/E055567/1 and the Highways Agency for access to MIDAS inductance loop data.
Appendix A. Details of the small wavenumber expansion
In §4, we introduced the perturbation expansion for λ+ in terms of the real coefficients , which we solved for up to second order. Furthermore, we showed that λ1<0 and the sign of λ2 controls stability. In this section, we prove that λ3>0 and λ4<0 when λ2>0, i.e. when the system is unstable.
Collecting the third-order terms in the small θ perturbation expansion of the dispersion relation given by equation (4.2), we find A1Using equation (3.3), we observe that A2since A3Thus, when λ2>0, all the grouped terms in parenthesis in equation (A1) have negative sign and therefore λ3>0. Collecting the fourth-order terms, we find A4Using the same argument as before, when λ2>0, we find A5and A6since all the grouped terms in square brackets are negative. Thus, when λ2>0, all the grouped terms in parenthesis in equation (A4) have positive sign and therefore λ4<0.
Appendix B. Details of the saddle selection procedure
We now give further details of the selection of the saddle point z* required in the signal velocity calculation as described in §5. This is a rather technical argument which the end-user may skip. Recall that we seek saddle points of ρ(z) defined by equation (5.8), which are thus roots of the cubic equation (5.11) parametrized by fh, , fv and κ. Ultimately, we are interested in calculating the correct choice of saddle as κ is varied. However, here we develop the prescription for fixed (but general) values of all of the parameters including κ.
Because the cubic equation (5.11) has real coefficients, ρ(z) has generically either three real saddles or a single real saddle and a complex conjugate pair. Our tactic is to analyse the local behaviour of ρ(z) at the real saddle(s) and to determine whether the direction of their curvature is suitable for employment in the steepest descent method. If no real saddle can be used, we infer that the complex conjugates must be chosen instead.
When a candidate saddle has thus been identified, the global geometry must be checked to ensure that the original contour can be deformed to pass through the saddle without crossing singularities. Moreover, we require that ϕ(z):=Re(ρ(z)) attains not only a local maximum but also its global maximum on the deformed contour 𝒞′ at the selected saddle point. In practice, this last point can be checked by inspecting the level-set geometry of ϕ(z) and ensuring that 𝒞′ does not re-cross the level sets that emanate from the saddle.
We shall seek deformed contours 𝒞′ that are symmetric about the real axis in the complex plane, and which thus cross the real axis orthogonally, i.e. parallel to the imaginary axis at the point of crossing. If ϕ has a local maximum in the imaginary direction at real z*, then by the Cauchy–Riemann equations, ϕ has a local minimum in the real direction. Real saddles of ϕ thus correspond to turning points of the real function B1and candidates for the steepest descent method correspond to local minima. Here x is real, and we define B2and B3The key simplification is that we may analyse the real geometry of equation (B1) without solving the cubic equation (5.11), but rather by studying its (real) singularities, which are of two types:
— P: the single zero of p(x). Note as .
— Q: the real zeroes of q(x), if there are any. These correspond to the platoon eigenvalues that solve equation (3.5). Note as .
Note also that as and as . The geometry of ϕ(x) may thus be classified according to the order of the singularities on the real line, which we may describe intuitively via the notation QQP, QPQ, PQQ or P. Here, case P denotes when the platoon eigenvalues are complex.
We now analyse how the singularity structure depends on the choice of parameters fh, , fv and κ. To simplify the study, we observe that a non-dimensional time can be introduced, for example, in equation (3.4). This eliminates the parameter fh and introduces rescaled parameters and and subsequently . Hence, in what follows, we set fh=1 and analyse how the problem depends on , fv and κ only, where the ‘tildes’ are dropped for notational convenience.
The order of the singularities of ϕ(x) is independent of κ, and so may be described entirely in terms of the parameter plane (figure 9). Here, the boundaries between the different classes can be described analytically, and in particular, it may be proved that only cases QQP and P occur for parameters corresponding to instability. For these cases, the putative geometries for ϕ and a notation to describe them are shown in figure 10.
The next task is to investigate which of the geometries in figure 10 are actually achieved, as κ is varied. This is a problem with intriguing mathematical structure, but to simplify matters, we have addressed it by numerical search. The results are summarized in figure 11. In each of the sub-plots, κ is held fixed, and the geometry in the plane is analysed. To focus attention on parameters corresponding to instability, the horizontal axis in these pictures is rescaled so that the unstable regime is a rectangular strip that fills the scope of each sub-plot.
— The QQP cases (0,1,0,0) and (0,1,0,2) and the P cases (1,0) and (1,2) dominate for most parameters. For these, the simple saddle selection rules as presented in §5 apply.
— There is a thin slice of parameters confined to 0<κ<0.1 (or in the dimensional units used in the main body of the paper) where cases (0,1,2,0) and (3,0) occur. For these, the middle saddle should actually be selected. To clarify: if the real saddles are given by z1<z2<z3, then z2 should be selected. These rare cases are contrary to the usual prescription given in §5.
— The cases (0,3,0,0) and (2,1,0,0), which have no admissible saddles, do not occur. In fact, this result can be proved rigorously by considering the small and large κ limits of the cubic equation (5.11) and its derivative with respect to κ; however, this analysis is beyond the present discussion.
Finally, figure 12 shows the level-set geometry of ϕ(z) and the saddle selection procedure and deformed contours for each of the six geometries that occur. These plots establish that the rules described above provide saddles and contours that satisfy the global geometric constraints of the steepest descent method.
- Received August 20, 2010.
- Accepted January 17, 2011.
- This journal is © 2011 The Royal Society