Stability and instability in a class of car following model on a closed loop

Physica A 388 (2009) 2476–2482

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Stability and instability in a class of car following model on a closed loop Alan McKee a , Mark McCartney b,∗ a

Pembroke College, Cambridge, CB2 1RF, United Kingdom

b

School of Computing and Mathematics, University of Ulster, Newtownabbey, BT37 0QB, United Kingdom

article

info

Article history: Received 9 October 2008 Received in revised form 8 January 2009 Available online 26 February 2009 PACS: 45.30.+s 45.50.Jf Keywords: Traffic flow Time delay differential equations

a b s t r a c t A velocity-matching car following model is modified to represent the motion of n vehicles travelling on a closed loop. Each vehicle is given a preferred velocity profile, which it attempts to achieve while also attempting to maintain a zero relative velocity between itself and the vehicle in front. The crucial distinctive of the looped model, as opposed to ‘non-looped’ models, is that the last vehicle in the stream is itself being followed by the lead (first) vehicle. The model gives rise to a system of n coupled time delay differential equations which are solved approximately (using a Taylor series expansion in time delay) and numerically using a fourth-order Runge–Kutta routine. The stability of the model is considered and an analytic form of the stable region in parameter space is found in the limit as n approaches infinity. © 2009 Elsevier B.V. All rights reserved.

1. Introduction Historically the mathematical modelling of road traffic flow has been approached in broadly two ways: macroscopic and microscopic models. Macroscopic models ignore the behaviour of individual vehicles and concentrate instead on global properties of traffic flow, such as traffic density (vehicles per metre) and flow rate past a point on the road (vehicles per second). In classic examples of macroscopic modelling, traffic is modelled as a continuous compressible fluid or a kinetic gas. Microscopic models on the other hand consider the behaviour of individual vehicles and their motion in relation to their nearest neighbours. Overviews of the modelling approaches used in traffic dynamics can be found elsewhere [1–5]. Recent work, see for example Refs. [6–11], has investigated the behaviour of traffic moving on a closed circuit or loop using both microscopic and macroscopic models. The work discussed in this paper applies this idea to a traditional linear stimulus–response car following model and explores what impact placing the vehicles in a closed loop has on the resulting behaviour of the system. 2. The Gazis–Herman–Rothery class of car following model Consider a column of n vehicles moving along a stretch of road where passing is not possible. If we label the vehicle at the head of the column i = 1, the next vehicle i = 2, etc., then a well established way [5] to model such a system is via a set of coupled differential equations of the form x¨ i (t ) = λ

∗

(˙xi (t ))m (˙xi−1 (t − T ) − x˙ i (t − T )) , (xi−1 (t − T ) − xi (t − T ))l

i = 2, . . . n

Corresponding author. E-mail address: [email protected] (M. McCartney).

0378-4371/$ – see front matter © 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.physa.2009.02.019

(1)

A. McKee, M. McCartney / Physica A 388 (2009) 2476–2482

2477

where xi (t ) is the position of the ith vehicle at time t, relative to a fixed position upstream, λ is a (positive) sensitivity coefficient measuring how strongly the following driver responds to the stimuli of the relative separation between the vehicles, (xi−1 (t − T ) − xi (t − T )), their relative velocities, (˙xi−1 (t − T ) − x˙ i (t − T )), and the current driving speed, x˙ i (t ). The indices m and l are positive constants (usually chosen to be integers) and T is the time it takes a driver (and vehicle) to respond to observed changes in relative velocity and position with respect to the (i − 1)th vehicle in front. For obvious reasons such models are described as car following models. Typically the motion of the lead vehicle, i = 1, is given by some prescribed function of velocity and the subsequent behaviour of the following vehicles is analysed. The presence of chaotic behaviour in (1) has been sought without success [12]. However it has been found that increasing the complexity of (1) by the addition of a nonlinear inter-vehicle spacing term [13], or simplifying (1) by setting m = l = 0, but allowing vehicles to overtake each other [9] have been found to give rise to chaotic motion. The simplest form of (1), the so-called velocity matching car following model, was first introduced in the late 1950’s by Chandler, Herman and Montroll [14], and can be expressed in terms of vehicle velocities as dui (t )

(2) = λ [ui−1 (t − T ) − ui (t − T )] dt where ui (t ) = x˙ i (t ) is the velocity of the ith vehicle at time t. Investigating (2) on a loop or ring road naturally leads us to allow the first (lead) vehicle in the stream to have a preferred velocity profile, w1 (t ), which the driver tries to achieve while also seeking to follow the vehicle in front. Thus he now accelerates at a rate which varies with both the relative velocity between his car and the vehicle in front and the difference between his current and preferred velocities at time t: du1 (t ) dt dui (t )

= λ [un (t − T ) − u1 (t − T )] + α [w1 (t − T ) − u1 (t − T )] (3)

= λ [ui−1 (t − T ) − ui (t − T )]

∀ i = 2, . . . , n

dt where α > 0. The behaviour of this model has been investigated in Ref. [10] and the work discussed here extends this model so that all n cars have a preferred velocity. 3. A looped model where all vehicles have a preferred velocity profile Generalising (3) so that all vehicles have a preferred velocity profile yields du1 (t ) dt dui (t )

= λ [un (t − T ) − u1 (t − T )] + α [w1 (t − T ) − u1 (t − T )] (4)

= λ [ui−1 (t − T ) − ui (t − T )] + α [wi (t − T ) − ui (t − T )]

dt or in matrix format

∀ i = 2, . . . , n

˙ (t ) = Au (t − T ) + α w (t − T ) u

(5)

where

− (λ + α)  λ    A= 0  ..  .

···

0

− (λ + α)

0

..

λ .. .

− (λ + α) .. .

···

0

0

0

.

λ



0

   .   

.

.. .

. λ

− (λ + α)

.. ..

0

(6)

Although the velocity matching model described by (2) and its looped variants (3) and (4) are extremely simplified models of the car following process as described by (1), it has the advantage of allowing us to isolate the effects of velocity matching behaviour on a closed loop, and, as is shown later in this paper, investigate some of the properties of the model analytically. In general, (4) will be solved numerically. However, as an alternative to the full numerical solution, we can obtain an approximate solution to the system for small values of T by using Taylor series expansions of ui (t − T ) and wi (t − T ) about t. By retaining the first two terms of the Taylor series expansions we obtain a system of differential equations of the form

  dw 1 ( t ) + α w1 ( t ) − T dt dt dt dt   dui (t ) dui−1 (t ) dui (t ) dw i ( t ) = λui−1 (t ) − λT − (λ + α) ui (t )+(λ + α) T +α wi (t )− T du1 (t )

dt

= λun (t ) − λT

dun (t )

dt

− (λ + α) u1 (t )+(λ + α) T

du1 (t )

dt

dt

(7)

∀ i = 2, . . . , n

2478

A. McKee, M. McCartney / Physica A 388 (2009) 2476–2482

Fig. 1. Results for Taylor series expansion (7) and Runge–Kutta (RK4) solution of (4) for a two-car ring model in which w1 (t ) = 20 ms−1 , w2 (t ) = 10 ms−1 , λ = 0.3 s−1 , α = 0.8 s−1 and T = 0.7 s, the time step in the fourth-order Runge–Kutta routine, ∆t = 0.01 s. Dashed line: Taylor series expansion solution for u1 (t ). Dashed–dotted line: RK4 solution for u1 (t ). Solid line: Taylor series expansion solution for u2 (t ). Dotted line: RK4 solution for u2 (t ). Note the initial unphysical backwards motion of the Taylor series expansion solution for u2 (t ).

or alternatively in matrix form as

˙ (t ) = (I + T A)−1 Au (t ) + (I + T A)−1 w ˜ (t ) u

(8)

  where A is as defined by (6) and w ˜ i (t ) = α wi (t ) − T dwdti (t ) . We expect a solution of the form u (t ) =

n X

ci ki eβi t + up (t )

(9)

i =1

where the βi are the solutions of the characteristic equation det (I + T A)−1 A − β I = 0




(10)

and ki are the corresponding eigenvectors. For a system of two vehicles it is easy to show that

α αT − 1 2λ + α β2 = (2λ + α) T − 1 β1 =

(11)

and hence the complementary function admits exponentially growing solutions when T > Tc =

1 2λ + α

.

When n = 2, w (t ) = w1 u1 ( t ) = − u2 ( t ) = −

w1 + w2 2

w1 + w2 2

(12)


T w2 is constant, and the initial condition u (0) = 0, the solution of (9) is α

t

2λ+α α (w1 − w2 ) (2λ+α) (λ + α) w1 + λw2 t T −1 + e 2 (2λ + α) 2λ + α 2λ+α α (w1 − w2 ) (2λ+α) λw 1 + (λ + α) w2 t T −1 + + e . 2 (2λ + α) 2λ + α

e αT −1 − α t e αT −1

(13)

Fig. 1 shows the results for a system of two vehicles in a ring for the fourth-order Runge–Kutta solution of (5) and the analytic solution of (8). The time delay is taken to be T = 0.7 s with λ = 0.3 s−1 , α = 0.8 s−1 , which from (12) gives Tc = 5/7 s. The preferred velocities are w1 (t ) = 20 ms−1 , w2 (t ) = 10 ms−1 . As Fig. 1 illustrates, even for values of T quite close to Tc the Taylor series approach can provide a reasonably good approximation to the RK4 solution for larger times. Indeed, at a modelling level we may be tempted to prefer it because the solutions appear to approach the post-transient behaviour monotonically. However, closer observation of Fig. 1 reveals that

A. McKee, M. McCartney / Physica A 388 (2009) 2476–2482

2479

the Taylor series approximation predicts unphysical backwards movement of the second vehicle for times close to t = 0. That such unphysical initial backwards motion can be a general feature of the approximate time delay model can be seen by considering (8) at time t = 0 with the initial condition u (0) = 0. In such a case, (8) can be written as

˙ (0) = (I + T A)−1 w ˜ (0) u

(14)

which has the solution

α u˙ i (0) = n a − bn

n −i X

a

n−j j−1

b

n X

˙ i+j (0) + wi+j (0) − T w


a

n −j j −1

b

˙ i+j−n (0) wi+j−n (0) − T w

!


(15)

j=n−i+1

j=1

where a = λT

(16)

b = (λ + α) T − 1. If we then consider the simple case where

w1 = w wi = 0 ∀ i = 2 . . . n

(17)

then (15) becomes u˙ i (0) = α

ai − 1 b n − i an − bn

.

(18)

Thus in (18) if b < 0 and |b| > |a| we will have initial negative acceleration when n − i is even, and if b < 0 and |b| < |a| we will have initial negative acceleration when n − i is odd. Given the initial condition u (0) = 0, this negative acceleration will lead to initial backward motion. Apart from the Taylor series approximation to (5), a further simplification is to set the time delay T = 0, yielding

˙ (t ) = Au (t ) + α w (t ) . u

(19)

If we now choose each vehicle to have a constant preferred velocity, w becomes a constant vector and (19) has a solution of the form (9), where the βi are the solutions of the characteristic equation det (A − β I) = 0

(20)

and ki are the corresponding eigenvectors. The solutions of (20) are

  2π ij βj = λ e n − 1 − α
where we note that R βj < 0 ∀j = 1 . . . n.

(21)

If w is a constant vector then the particular integral of (19) is given by up = −α A−1 w

(22)

and the ith entry of up is found to be given by upi =

"

α  (α + λ) 1 −


λ n α+λ

i  X

 j =1

λ α+λ

i−j

wj +

n  X j=i+1

λ α+λ

n+i−j

# wj .

(23)

Thus the post-transient behaviour of the ith vehicle for the zero time delay system (19) with constant preferred velocities is simply the weighted mean of the preferred velocities with the weighting factor on the velocity of the jth vehicle in front λ of the ith being α+λ .


j

4. Stability of the model We investigate the stability of (5) by investigating solutions of the matrix equation

˙ (t + T ) = Au (t ) . u

(24)

Using a trial complementary function of the form ucf (t ) = keβ t

(25)

2480

A. McKee, M. McCartney / Physica A 388 (2009) 2476–2482

then (24) can be expressed as

β eβ T keβ t = Akeβ t ;

(26)

βi T

thus ki and βi = βi e are the eigenvectors and eigenvalues of A and we expect a solution of the form (9). The eigenvalues of A are given by the solution of the characteristic equation ∗

det A − β ∗ I = 0,




(27)

which reduces to

µ + γ + zez


n

− γn = 0

(28)

where z = βT

γ = λT µ = αT .

(29)

For the system to be defined as locally stable the eigenvalues, βi , must be real and negative ∀i = 1, . . . , n. Eq. (28) has the solution zez = γ ei

2π k n

− γ − µ for k = 0, 1, . . . , n − 1,

and, noting that z ∈ R ⇒ zez ∈ R and for n ≥ 3, γ e n ≥ 3. For n = 2, the system is locally stable when

µ<

i 2nπ

(30)

− γ − µ 6∈ R, we see that the system is never locally stable for

1

− 2γ . (31) e Since local stability can only occur in a system of two vehicles on a closed loop, we extend our stability analysis to investigate the conditions under which the system possesses what may be described as weak local stability. In this case the eigenvalues, βi , must have negative real parts. In order to analytically investigate the stability region boundary in (γ , µ) space then we use a trial version for the eigenvalue of the form β = iy

(32)

where y ∈ R. Substituting (32) into (30) and equating real and imaginary parts gives Y sin Y = µ − γ Y cos Y = γ sin





 cos

2π k

2π k n



 −1



(33)

n

where Y = yT ∈ R. Eqs. (33) then give a parameterisation of the boundary curves for the stability region in (γ , µ) space. We note that the boundary given by k = j is the same as that given by k = n − j, and hence the weak local stability region is   bounded by at most 2n + 1 curves along with the two axes. Fig. 2 shows the weak local stability region in (γ , µ) space for a system of 11 vehicles; it plots the curves parameterised by (33), and the region of weak local stability is shaded in grey. We note that, when n = 2, (33) yields

µ<

π

− 2γ . (34) 2 Comparing (34) with (31) we note that, as expected, the region of weak local stability is a superset of the region of local stability. 5. System stability as n → ∞ We next investigate the stability of the model as n → ∞. To do this we let θ = 2πn k , and in the limit n → ∞ we treat θ ∈ [0, π ] as a continuous variable and find the envelope of the set of curves parameterised by θ . Changing from (γ , µ) co-ordinates to polar co-ordinates

γ = r cos ϕ µ = r sin ϕ

(35)

then (33) can be written as Y sin Y = r sin ϕ − r cos ϕ (cos θ − 1) Y cos Y = r cos ϕ sin θ .

(36)

A. McKee, M. McCartney / Physica A 388 (2009) 2476–2482

2481

Fig. 2. Region of weak local stability (grey) for the system (4). The curves are the solutions to (33) for a system of 11 vehicles.

Fig. 3. Boundary of the weak local stability region for the system (4) in the limit n → ∞ as given by (37).

To find the boundary of the weak local stability region in the limit n → ∞, we then minimise r with respect to θ for fixed φ . Differentiating (36) with respect to θ and setting ddrθ = 0 leads to the following parameterisation for the boundary curve:

q γ =

1 + Y2




1 + tan2 Y

Y + tan Y


Y cos Y

 µ = Y sin Y − γ 1 − q

 1 − Y tan Y 1 + Y2




1 + tan2 Y




(37)

.

The form of (37) is shown in Fig. 3. Note that for a fixed value of γ (e.g. in Fig. 3 consider γ ' 0.55) the system is unstable for very low values of µ and then exhibits weak stability as µ increases. However, as µ continues to increase the system becomes unstable again. This phenomenon results from the competition between the sensitivity, λ = γ /T , of the vehicle towards the vehicle in front of him and his sensitivity α = µ/T towards his own preferred velocity. Thus for a range of values this competition is constructive, resulting is stable behaviour. 6. Concluding remarks In this paper we have investigated the effect of taking a velocity matching car following model, placing the vehicles on a closed loop and giving each car a preferred velocity. A range of analytic and numerical results have been presented to

2482

A. McKee, M. McCartney / Physica A 388 (2009) 2476–2482

analyse the resulting mathematical model. We find that the approximate solution technique of expanding the time delay as a Taylor series, and then solving the resulting coupled ordinary differential equations analytically, can lead to unphysical initial backwards motion for certain parameter choices. Analysis of the local stability of the model shows that the system only exhibits local stability under restricted conditions when there are two vehicles on the loop. If there are three or more vehicles on the loop then the system is never locally stable. Extending the stability analysis to allow for disturbances which are damped, but may be oscillatory (which we describe as weak local stability), we find that the model will admit weak locally stable solutions for any number of vehicles and that there is a limit for the stability region as the number of vehicles n → ∞. In this latter case, the result is pleasing both mathematically and physically. It is pleasing mathematically, because we are able to give a result about the stability of an infinite set of coupled time delay differential equations. It is pleasing physically because as the number of vehicles on the looped network increases we approach a well defined region of parameter space for which the model is well behaved. As noted earlier, the velocity matching model studied in this paper is the simplest example of the class of car following model as described by (1). In particular, no reference has been made to the influence of relative separations of vehicles, or to the possibility of collisions. It is intended that future work will address these issues by investigating more complex models. Acknowledgements We acknowledge the advice of three anonymous referees in the writing of the final version of this paper. The work carried out in this paper was partially funded by the Nuffield Foundation under grant number URB/35452. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14]

B.S. Kerner, The Physics of Traffic, Springer, Berlin, New York, 2004. D.C. Gazis, Traffic Theory, Kluwer Academic Publishers, Boston, 2002. D. Helbing, Reviews of Modern Physics 73 (2001) 1067–1141. D. Chowdhury, A. Schadschneider, K. Nishinari, Physics of Life Reviews 2 (2005) 318–352. M. Brackstone, M. McDonald, Transportation Research F 2 (2000) 181–196. M. Bando, K. Hasebe, A. Nakayama, A. Shibata, Y. Sugiyama, Physical Review E 51 (1995) 1035–1042. L.C. Davis, Physica A 319 (2003) 557–567. W. Dahui, W. Ziqiang, F. Ying, Physical Review E 76 (2007) 016105. S. Jamison, M. McCartney, Chaos 17 (2007) 033116. S. Jamison, M. McCartney, Investigating a class of car following model on a ring, in: Benjamin Heydecker (Ed.), Mathematics in Transport, Elsevier, 2007, pp. 97–110. M. McCartney, Transportation Research B (2007) (in press). D. Jarrett, Z. Xiaoyan, in: A.J. Crilly, R.A. Earnshaw, H. Jones (Eds.), The Dynamic Behaviour of Road Traffic Flow: Stability or Chaos? in: Applications of Fractals and Chaos, Springer-Verlag, Berlin, Heidelberg, New York, 1993. D.J. Low, P.S. Addison, Nonlinear Dynamics 16 (1998) 127–151. F.E. Chandler, R. Herman, E.W. Montroll, Operations Research 6 (1958) 165–184.