The Modified Optimal Velocity Model: Stability Analyses and Design Guidelines
11 The Modified Optimal Velocity Model:Stability Analyses and Design Guidelines
Gopal Krishna Kamath, Krishna Jagannathan and Gaurav RainaDepartment of Electrical Engineering, Indian Institute of Technology Madras, Chennai 600 036, IndiaEmail: { ee12d033, krishnaj, gaurav } @ee.iitm.ac.in Abstract
Reaction delays are important in determining the qualitative dynamical properties of a platoon of vehicles travelingon a straight road. In this paper, we investigate the impact of delayed feedback on the dynamics of the ModifiedOptimal Velocity Model (MOVM). Specifically, we analyze the MOVM in three regimes – no delay, small delayand arbitrary delay. In the absence of reaction delays, we show that the MOVM is locally stable. For small delays,we then derive a sufficient condition for the MOVM to be locally stable. Next, for an arbitrary delay, we derive thenecessary and sufficient condition for the local stability of the MOVM. We show that the traffic flow transits fromthe locally stable to the locally unstable regime via a Hopf bifurcation. We also derive the necessary and sufficientcondition for non-oscillatory convergence and characterize the rate of convergence of the MOVM. These conditionshelp ensure smooth traffic flow, good ride quality and quick equilibration to the uniform flow. Further, since a Hopfbifurcation results in the emergence of limit cycles, we provide an analytical framework to characterize the type of theHopf bifurcation and the asymptotic orbital stability of the resulting non-linear oscillations. Finally, we corroborateour analyses using stability charts, bifurcation diagrams, numerical computations and simulations conducted usingMATLAB.
Index Terms
Transportation networks, car-following models, time delays, stability, convergence, Hopf bifurcation.
I. I
NTRODUCTION
Intelligent transportation systems constitute a substantial theme of discussion on futuristic smart cities. In thiscontext, self-driven vehicles are a prospective solution to address traffic issues such as resource utilization andcommute delays; see [1, Section 5.2], [2]–[4] and references therein. To ensure that these objectives are met, inaddition to ensuring human safety, the design of control algorithms for these vehicles becomes important. To thatend, it is imperative to have an in-depth understanding of human behavior and vehicular dynamics. This has led tothe development and study of a class of dynamical models known as the car-following models [5]–[10].
A conference version of this paper appeared in Proceedings of the rd Annual Allerton Conference on Communication, Control andComputing, pp. 538-545, 2015. DOI: 10.1109/ALLERTON.2015.7447051 a r X i v : . [ c s . S Y ] J u l Feedback delays play an important role in determining the qualitative behavior of dynamical systems [11]. Inparticular, these delays are known to destabilize the system and induce oscillatory behavior [10], [12]. In thecontext of human-driven vehicles, predominant components of the reaction delay are psychological and mechanicalin nature [12]. In contrast, the delay in self-driven vehicles arise due to sensing, communication, signal processingand actuation, and are envisioned to be smaller than human reaction delays [13].In this paper, we investigate the impact of delayed feedback on the qualitative dynamical properties of a platoonof vehicles traveling on a straight road. Specifically, we consider each vehicle’s dynamics to be modeled by theModified Optimal Velocity Model (MOVM) [10]. Motivated by the wide range of values assumed by reactiondelays in various scenarios, we analyze the MOVM in three regimes; namely, ( i ) no delay, ( ii ) small delay and ( iii ) arbitrary delay. In the absence of delays, we show that the MOVM is locally stable. When the delays arerather small, as in the case of self-driven vehicles, we derive a sufficient condition for the local stability of theMOVM using a suitable approximation. For the arbitrary-delay regime, we analytically characterize the region oflocal stability for the MOVM.In the context transportation networks, two additional properties are of practical importance; namely, ride quality(lack of jerky vehicular motion) and the time taken by the platoon to attain the desired equilibrium when perturbed.Mathematically, these translate to studying the non-oscillatory property of the MOVM’s solutions and the rate oftheir convergence to the desired equilibrium. In this paper, we also characterize these properties for the MOVM.In the context of human-driven vehicles, model parameters generally correspond to human behavior, and hencecannot be “tuned” or “controlled.” However, our work enhances phenomenological insight into the emergence andevolution of traffic congestion. For example, a peculiar phenomenon known as the “phantom jam” is observed onhighways [7], [8]. Therein, a congestion wave emerges seemingly out of nowhere and propagates up the highwayfrom the point of its origin. Such an oscillatory behavior in the traffic flow has typically been attributed to achange in the driver’s sensitivity, such as a sudden deceleration; for details, see [7], [8]. In general, feedback delaysare known to induce oscillations in state variables of dynamical systems [10], [12]. Since the MOVM explicitlyincorporates feedback delays, and relative velocities and headways constitute state variables of the MOVM, ourwork provides a theoretical basis for understanding the emergence and evolution of oscillatory phenomena such as“phantom jams.” In particular, our work serves to highlight the possible role of reaction delays in the emergenceof oscillatory phenomena in traffic flows. More generally, our results reveal an important observation: the trafficflow may transit into instability due to an appropriate variation in any subset of model parameters. To capture thiscomplex dependence of stability on various parameters, we introduce an exogenous, non-dimensional parameter inour dynamical model. We then analyze the behavior of the resulting system as the exogenous parameter is pushedjust beyond the stability boundary. We show that non-linear oscillations, termed limit cycles , emerge in the trafficflow due to a Hopf bifurcation .In the context of self-driven vehicles, reaction delays are expected to be smaller than their human counterparts [13].Hence, it would be realistically possible to achieve smaller equilibrium headways [1, Section 5.2]. This would, inturn, vastly improve resource utilization without compromising safety [3]. In this paper, based on our theoreticalanalyses, we provide some design guidelines to appropriately tune the parameters of the so-called “upper longitudinal control algorithm” [1, Section 5.2]. Mathematically, our analytical findings highlight the quantitative impact ofdelayed feedback on the design of control algorithms for self-driven vehicles. Specifically, our design guidelinestake into consideration various aspects of the longitudinal control algorithm such as stability, good ride qualityand fast convergence of the traffic to the uniform flow. In the event that the traffic flow does lose stability, ourdesign guidelines help tune the model parameters with an aim of reducing the amplitude and angular velocity ofthe resultant limit cycles.
A. Related work on car-following models
The motivation for our paper comes from the key idea behind the Optimal Velocity Model (OVM) proposed byBando et al. in [14]. However, the model considered therein was devoid of reaction delays. Thus, a new modelwas proposed in [6] to account for the drivers’ delays. Therein, the authors also claimed that these delays werenot central to capturing the dynamics of the system. In response, Davis showed via numerical computations thatreaction delays indeed played an important part in determining the qualitative behavior of the OVM [15]. This ledto a further modification to the OVM in [16]. However, this too did not account for the delay arising due to avehicle’s own velocity.It was shown in [17] that the OVM without delays loses local stability via a Hopf bifurcation. For the OVM withdelays, [18] performed an initial numerical study of the bifurcation phenomenon before supplying an analyticalproof in [9]. The theme issue on “Traffic jams: dynamics and control” [19] highlights the growing interest ina dynamical systems viewpoint of transportation networks. A recent exposition of linear stability analysis in thecontext of car-following models can be found in [20].From a vehicular dynamics perspective, most upper longitudinal controllers in the literature assume the lowercontroller’s dynamics to be well modeled by a first-order control system, in order to capture the delay lag [1, Section5.3]. The upper longitudinal controllers are then designed to maintain either constant velocity, spacing or time gap;for details, see [21] and the references therein. Specifically, it was shown in [21] that synchronization with thelead vehicle is possible by using information only from the vehicle directly ahead. This reduces implementationcomplexity, and does not mandate vehicles to be installed with communication devices.However, in the context of autonomous vehicles, communication systems are required to exchange various systemstates required for the control algorithm. This information is used either for distributed control [21] or coordinatedcontrol [22]. Formation and platoon stabilities have also been studied considering information flow among thevehicles [23], [24]. For an extensive review, see [25]. Usage of communication systems is also known to mitigatephantom jams [26].At a microscopic level, Chen et al. proposed a behavioral car-following model based on empirical data thatcaptures phantom jams [27]. Therein, the authors showed statistical correlation in drivers’ behavior before andduring traffic oscillations. However, no suggestions to avoid phantom jams were offered. To that end, Nishi et al. developed a framework for “jam-absorbing” driving in [28]. A “jam-absorbing vehicle” appropriately varies itsheadway with the aim of mitigating phantom jams. This work was extended by Taniguchi et al. [29] to include car-following behavior. Therein, the authors also numerically constructed the region in parameter space that avoidsformation of secondary jams.In the context of platoon stability, it has been shown that well-placed, communicating autonomous vehicles maybe used to stabilize platoons of human-driven vehicles [30]. More generally, the platooning problem has beenstudied as a consensus problem with delays [31]. Such an approach aids the design of coupling protocols betweeninteracting agents (in this context, vehicles). In contrast, we provide design guidelines to appropriately chooseprotocol parameters, for a given coupling protocol.
B. Our contributions
Our contributions are as follows.(1) We derive a variant of the OVM for an infinitely-long road – called the MOVM – and analyze it in threeregimes; namely, ( i ) no delay, ( ii ) small delay and ( iii ) arbitrary delay. We prove that the ideal case ofinstantaneously-reacting drivers is locally stable for all practically significant parameter values. We then derivea stability condition for the small-delay regime by conducting a linearization on the time variable.(2) For the case of an arbitrary delay, we derive the necessary and sufficient condition for the local stability ofthe MOVM. We then prove that, upon violation of this condition, the MOVM loses local stability via a Hopfbifurcation.(3) We provide an analytical framework to characterize the type of the Hopf bifurcation and the asymptotic orbitalstability of the emergent limit cycles using Poincar´e normal forms and the center manifold theory.(4) In the case of human-driven vehicles, our work enhances phenomenological insight into the emergence andevolution of traffic congestion. For example, the Hopf bifurcation analysis provides a mathematical frameworkto offer a possible explanation for the observed “phantom jams” [10]. In the case of self-driven vehicles, ourwork offers suggestions for their design guidelines.(5) We derive a necessary and sufficient condition for non-oscillatory convergence of the MOVM. This is useful inthe context of a transportation network since oscillations lead to jerky vehicular movements, thereby degradingride quality and possibly causing collisions.(6) We characterize the rate of convergence of the MOVM, thereby gaining insight into the time required forthe platoon to equilibrate, when perturbed. Such perturbations occur, for instance, when a vehicle departsfrom a platoon. Therein, we also bring forth the trade-off between the rate of convergence and non-oscillatoryconvergence of the MOVM.(7) We corroborate the analytical results with the aid of stability charts, bifurcation diagrams, numerical compu-tations and simulations performed using MATLAB.The remainder of this paper is organized as follows. In Section II, we summarize the OVM and derive theMOVM. In Sections III, IV and V, we characterize the stable regions for the MOVM in no-delay, small-delayand arbitrary-delay regimes respectively. We then derive the necessary and sufficient condition for non-oscillatoryconvergence of the MOVM in Section VI, and characterize its rate of convergence in Section VII. In Section VIII, we present the local Hopf bifurcation analysis for the MOVM. In Section IX, we corroborate our analyses usingMATLAB simulations before concluding the paper in Section X.II. M ODELS
In this section, we first provide an overview of the setting of our work. We then briefly explain the OVM, beforeending the section by deriving the MOVM.
A. The setting
We consider N + 1 idealistic vehicles (with length) traveling on an infinitely long, single-lane road with noovertaking. The lead vehicle is indexed with , the vehicle following it with , and so on. The acceleration of eachvehicle is updated based on a combination of its position, velocity and acceleration as well as those correspondingto the vehicle directly ahead. We use x i ( t ) , ˙ x i ( t ) and ¨ x i ( t ) to denote the position, velocity and acceleration of thevehicle indexed i at time t respectively. We also assume that the lead vehicle’s acceleration and velocity profilesare known. Specifically, we only consider leader profiles that converge to ¨ x = 0 and < ˙ x < ∞ in finite time;that is, there exists T < ∞ such that ¨ x ( t ) = 0 , ˙ x ( t ) = ˙ x > , ∀ t ≥ T . We also use the terms “driver” and“vehicle” interchangeably throughout. B. The Optimal Velocity Model (OVM)
The OVM, proposed by Bando et al. in [14], is based on the key idea that each vehicle in a platoon tries to attainan “optimal” velocity, which a function of its headway. Hence, each vehicle updates its acceleration proportionalto the difference between this optimal velocity and its own velocity. This was modified in [6] to account for thereaction delay. For N vehicles traveling on a circular loop of length L units, the dynamics is captured by [6] ¨ x ( t ) = a ( V ( x N ( t − τ ) − x ( t − τ )) − ˙ x ( t − τ )) , ¨ x i ( t ) = a ( V ( x i − ( t − τ ) − x i ( t − τ )) − ˙ x i ( t − τ )) , (1)for i ∈ { , · · · , N } . Here, a > is the drivers’ sensitivity coefficient, τ is the common reaction delay and V : R + → R + is called the Optimal Velocity Function (OVF). As pointed out in [32], an OVF satisfies:(i) Monotonic increase,(ii) Bounded above, and,(iii) Continuous differentiability.Let V max = lim y →∞ V ( y ) . The limit exists as a consequence of (i) and (ii) above. Also, (iii) ensures that an OVFwill also be invertible.
C. The Modified Optimal Velocity Model (MOVM)
Next, we derive a variant of the OVM for our setting. To that end, we define the relative spacing (or headway)and the relative velocity as y i ( t ) = x i − ( t ) − x i ( t ) and v i ( t ) = ˙ y i ( t ) = ˙ x i − ( t ) − ˙ x i ( t ) respectively. Substitutingthese in (1), we obtain the following system after some algebraic manipulations ˙ v ( t ) = ¨ x ( t ) + a ( ˙ x ( t − τ ) − V ( y ( t − τ )) − v ( t − τ )) , ˙ v k ( t ) = a ( V ( y k − ( t − τ k − )) − V ( y k ( t − τ k )) − v k ( t − τ k )) , ˙ y i ( t ) = v i ( t ) , (2)for i ∈ { , , · · · , N } and for k ∈ { , , · · · , N } . We refer to system (2) as the Modified Optimal Velocity Model(MOVM). Notice that, in contrast to (1), the MOVM no longer possesses the circular structure resulting from theperiodic boundary condition. Indeed, the second vehicle (with index ) now follows the lead vehicle rather thanthe vehicle with index N, since we consider the platoon to be traversing an infinitely long road.The MOVM is described by a system of Delay Differential Equations (DDEs). Since such systems are hard toanalyze, we obtain conditions for their local stability by analyzing them in the neighborhood of their equilibria.Such an analysis technique is called local stability analysis . To obtain the equilibrium, we equate the N derivativesin (2) to zero. This yields v ∗ i = 0 , y ∗ i = V − ( ˙ x ) , i = 1 , , · · · , N as an equilibrium for system (2). Therefore, tolinearize (2) about this equilibrium, we first consider a small perturbation u i ( t ) about the equilibrium of the relativespacing pertaining to vehicle indexed i. That is, u i ( t ) = y i ( t ) - y ∗ i . Next, we consider the Taylor’s series expansionof u i ( t ) , and set the leader’s profile to zero, to obtain the linearized model, given by ˙ v ( t ) = − du ( t − τ ) − av ( t − τ ) , ˙ v k ( t ) = du k − ( t − τ k − ) − du k ( t − τ k ) − av k ( t − τ k ) , ˙ u i ( t ) = v i ( t ) , (3)for i ∈ { , , · · · , N } and for k ∈ { , , · · · , N } . Here, d = aV (cid:48) ( V − ( ˙ x )) is the equilibrium coefficient, where theprime indicates differentiation with respect to a state variable. Henceforth, we denote ˜ d = V (cid:48) ( V − ( ˙ x )) . Therefore, d = a ˜ d. The MOVM is completely specified by the relative velocities v i ’s and the headways y i ’s. Therefore, the state ofthe MOVM at time “ t ” is given by S ( t ) = [ v ( t ) v ( t ) · · · v N ( t ) u ( t ) u ( t ) · · · u N ( t )] T ∈ R N . Thus, system (3)can be succinctly written in matrix form as ˙ S ( t ) = N (cid:88) k =0 A k S ( t − τ k ) , (4)where the matrices A k ∈ R N × N for each k . Also, τ is introduced for notational brevity and set to zero. Thematrices A k , k = 1 , · · · , N, are defined as follows. A = N × N N × N I N × N N × N . Here, N × N and I N × N denote zero and identity matrices of order N × N respectively. For ≤ k ≤ N − , wehave ( A k ) ij = − a, i = j = k, − d, i = k, j = N + k,d, i = k + 1 , j = k, , elsewhere , and ( A N ) ij = − a, i = j = N, − d, i = N, j = 2 N, , elsewhere . D. Optimal Velocity Functions (OVFs)
There are several functions that satisfy the properties mentioned in Section II-B. We mention four widely-usedOVFs [32], obtained by fixing a functional form for V ( · ) .(a) Underwood OVF: V ( y ) = V e − ymy . (b) Bando OVF: V ( y ) = V (cid:18) tanh (cid:18) y − y m ˜ y (cid:19) + tanh (cid:18) y m ˜ y (cid:19)(cid:19) . (c) Trigonometric OVF: V ( y ) = V (cid:18) tan − (cid:18) y − y m ˜ y (cid:19) + tan − (cid:18) y m ˜ y (cid:19)(cid:19) . (d) Hyperbolic OVF: V ( y ) = , y ≤ y ,V (cid:16) ( y − y ) n (˜ y ) n +( y − y ) n (cid:17) , y ≥ y . Here, V , y , y m , ˜ y and n are model parameters.As captured by [33, Figure ], the aforementioned OVFs behave similarly with varying headway. The followingare noteworthy: ( i ) The values attained by these OVFs, in the vicinity of the equilibrium, are almost the same, ( ii ) their slopes, evaluated at the equilibrium, are different. The linearized version of the MOVM, captured bysystem (4), brings forth the dependence on the slope via the variable ˜ d, and ( iii ) we make use of the Bando OVFthroughout this paper, except in Section VIII. Therein, we consider both the Bando OVF and the Underwood OVF,consistent with [10].We now proceed to understand the dynamical behavior of a platoon of cars running the MOVM. III. T
HE NO - DELAY REGIME
We first consider the idealistic case of instantaneously-reacting drivers. This results in zero reactions delays.Therefore, the model described by system (4) boils down to the following system of Ordinary Differential Equations(ODEs): ˙ S ( t ) = (cid:32) N (cid:88) k =0 A k (cid:33) S ( t ) . (5)We denote by A, the sum of matrices A k , which is known as the dynamics matrix . To characterize the stabilityof system (5), we require the eigenvalues of A to be negative [34, Theorem 5.1.1]. To that end, we compute itscharacteristic function as f ( λ ) = det ( λI N × N − A ) = det ( λ + a ) I N × N ˆ DI N × N λI N × N = 0 , where ˆ D is derived from the dynamics matrix A. The diagonal matrices of the above block matrix are invertible,and the off-diagonal matrices commute with each other. Hence, from [35, Theorem 3], the characteristic equationcan be simplified to ( λ + aλ − d ) N = 0 , which holds true if and only if λ + aλ − d = 0 . (6)Solving the above quadratic, we notice that the poles corresponding to system (5) will be negative if a > and ˜ d = V (cid:48) ( V − ( ˙ x )) > . We note that, from physical constraints, a > . Also, since V ( · ) is an OVF, itis monotonically increasing. Therefore, ˜ d > . Hence, for all physically relevant values of the parameters, thecorresponding poles will lie in the open left-half of the Argand plane. Thus, the MOVM is locally stable for allphysically relevant values of the parameters, in the absence of delays.IV. T
HE SMALL - DELAY REGIME
Having studied the MOVM in the absence of reaction delays, we now analyze it in the small-delay regime. Away to obtain insight for the case of small delays is to conduct a linearization on time. This would yield a systemof ODEs, which serves as an approximation to the original infinite-dimensional system (4), valid for small delays.We derive the criterion for such a system of ODEs to be stable, thereby emphasizing the design trade-off inherentamong various system parameters and the reaction delay.We begin by applying the Taylor’s series approximation to the time-delayed state variables thus: v i ( t − τ i ) ≈ v i ( t ) − τ i ˙ v i ( t ) , and u i ( t − τ i ) ≈ u i ( t ) − τ i ˙ u i ( t ) . Using this approximation for terms in (3), substituting v i ( t ) for ˙ u i ( t ) and re-arranging the resulting equations, we obtain the matrix equation B ˙ S ( t ) = A S ( t ) . (7)where the matrix A is the dynamics matrix, as defined in Section III, and B is a block matrix of the form B = B s N × N N × N I N × N , , where ( B s ) ij = − aτ i , i = j, , elsewhere . We note that, since B s is a diagonal matrix, so is B. Also, B is invertible if and only if aτ i (cid:54) = 1 , for each i. Thus,when aτ i (cid:54) = 1 , for each i, we define ˜ C = B − A, which is of the form ˜ C = ˜ C s ˜ C c I N × N N × N , , where ( C s ) ij = − a + dτ i − aτ i , i = j, − dτ j − aτ i , j = i − , , elsewhere , and ( C c ) ij = − d − aτ i , i = j, d − aτ i , j = i − , , elsewhere . For system (7) to be stable, the real part of eigenvalues of ˜ C must be negative [34, Theorem 5.1.1]. To that end,we compute its characteristic function as f ( λ ) = det ( λI N × N − ˜ C ) = det λI N × N − ˜ C s − ˜ C c − I N × N λI N × N = 0 . The diagonal matrices of the aforementioned block matrix are invertible, and the matrices in the second row thereincommute with each other. Hence, the characteristic equation simplifies to [35, Theorem 3] f ( λ ) = det (cid:16) λ ( λI N × N − ˜ C s ) − ˜ C c (cid:17) = 0 . On further simplification, this yields f ( λ ) = N (cid:89) i =1 (cid:0) (1 − aτ i ) λ + ( a − dτ i ) λ + d (cid:1) = 0 . For multiple terms in the above product to equal zero, their respective reaction delays must be equal. Such apossibility is not realistic, hence we ignore it. Therefore, for some i ∈ { , , · · · , N } , we have (1 − aτ i ) λ + ( a − dτ i ) λ + d = 0 . The roots of this quadratic equation are given by λ , = − ( a − dτ i ) ± (cid:112) ( a − dτ i ) − d (1 − aτ i )2(1 − aτ i ) . We now consider the following (exhaustive) cases. (1) Let aτ i > . Since d > , it follows that d (1 − aτ i ) < . Then, the eigenvalues are real. Further, one of theseeigenvalues will be positive and the other negative. Hence, we require aτ i < for system (7) to be stable.(2) Let ( a − dτ i ) ≥ d (1 − aτ i ) . Then, the eigenvalues are real. They are negative if and only if a − dτ i > ,i.e., ˜ dτ i < . Hence, we require ˜ dτ i < for system (7) to be stable.(3) Let ( a − dτ i ) < d (1 − aτ i ) . Then, the eigenvalues are complex. The real part of the eigenvalues will benegative if and only if a − dτ i > , i.e., ˜ dτ i < . Hence, we require ˜ dτ i < for system (7) to be stable.From the above cases, it is clear that system (7) is stable if and only if max( a, ˜ d ) τ i < , (8)for each i ∈ { , , · · · , N } . Recall that we obtained system (7) by truncating the Taylor’s series to first order.Hence, (8) is a sufficient condition for the local stability of the MOVM described by system (2), valid for smallvalues of the reaction delay. V. T
HE ARBITRARY - DELAY REGIME
Having studied system (2) in the no-delay and the small-delay regimes, in this section, we focus on the arbitrary-delay regime. We first derive the necessary and sufficient condition for the local stability of the MOVM. We thenshow that, upon violation of this condition, the corresponding traffic flow transits via a Hopf bifurcation to thelocally unstable regime.
A. Transversality condition
Hopf bifurcation is a phenomenon wherein, on appropriate variation of system parameters, a dynamical systemeither loses or regains stability because of a pair of conjugate eigenvalues crossing the imaginary axis in the Argandplane [11, Chapter 11, Theorem 1.1]. Mathematically, Hopf bifurcation analysis is a rigorous way of proving theemergence of limit cycles (isolated closed trajectory in state space) in non-linear dynamical systems.To determine if system (2) undergoes a stability loss via a Hopf bifurcation, we follow [37] and introduce anexogenous, non-dimensional parameter κ > . A general system of DDEs ˙ x ( t ) = f ( x ( t ) , x ( t − τ ) , · · · , x ( t − τ n )) , (9)is modified to ˙ x ( t ) = κf ( x ( t ) , x ( t − τ ) , · · · , x ( t − τ n )) , (10)with the introduction of the exogenous parameter. Note that ( i ) κ has no effect on the equilibrium of system (9),and ( ii ) we obtain system (9) by setting κ = 1 in system (10). We first linearize system (10) about its non-trivialequilibrium and derive its characteristic equation. We then search for a pair of conjugate eigenvalues on the imaginaryaxis in the Argand plane. This yields the necessary and sufficient condition for the local stability of system (10).Setting the exogenous parameter to unity then yields the necessary and sufficient condition for system (9). Theexogenous parameter so introduced helps simplify the requisite algebra and capture any interdependence among thesystem parameters. For the MOVM, introducing κ in (2) yields ˙ v ( t ) = ¨ x ( t ) + κa ( ˙ x ( t − τ ) − V ( y ( t − τ )) − v ( t − τ )) , ˙ v k ( t ) = κa ( V ( y k − ( t − τ k − )) − V ( y k ( t − τ k )) − v k ( t − τ k )) , ˙ y i ( t ) = κv i ( t ) , (11)for i ∈ { , , · · · , N } and for k ∈ { , , · · · , N } . We linearize this about the equilibrium v ∗ i = 0 , y ∗ i = V (cid:48) ( V − ( ˙ x )) , i = 1 , , · · · , N, and write it in matrix form to obtain ˙ S ( t ) = N (cid:88) k =0 ˜ A k S ( t − τ k ) , (12)where the matrices ˜ A k = κA k , for k = 1 , , · · · , N, where the matrices A k are as defined in Section II.The characteristic equation corresponding to system (12) is obtained as [34, Section 5.1] f ( λ ) = det (cid:32) λI N × N − N (cid:88) k =0 e − λτ k ˜ A k (cid:33) = 0 . The matrix in consideration is a block matrix of the form λI N × N − N (cid:88) k =0 e − λτ k ˜ A k = ˜ A ˜ B ˜ C ˜ D , where ˜ C = − κI N × N and ˜ D = λI N × N . Further, ˜ A is a diagonal matrix with the i th diagonal entry being λ + κae − λτ i , and ˜ B is a sparse lower-triangular matrix. Clearly, ˜ A and ˜ D are invertible, and ˜ C commutes with ˜ D .Therefore, the characteristic equation simplifies to the form [35, Theorem 3] f ( λ ) = det ˜ A ˜ B ˜ C ˜ D = det (cid:16) ˜ A ˜ D − ˜ B ˜ C (cid:17) = 0 . Simplifying the above expression, we obtain the characteristic equation pertaining to (12) as f ( λ ) = N (cid:89) i =1 ( λ + κaλe − λτ i + κ de − λτ i ) = 0 . (13)For multiple terms in the above product to equal zero, their respective reaction delays must be equal. Such apossibility is not realistic, hence we ignore it. Therefore, for some i ∈ { , , · · · , N } , we have λ + κaλe − λτ i + κ de − λτ i = 0 . (14)System (11) will be locally stable if and only if all the roots of (14) lie in the open left-half of the Argand plane [34,Theorem 5.1.1]. Therefore, we search for a conjugate pair of eigenvalues of (14) that crosses the imaginary axis inthe Argand plane. To that end, we substitute λ = jω in (14), with j = √− . We then equate the real and imaginaryparts to zero and obtain κaω sin ( ωτ i ) + κ d cos ( ωτ i ) = ω , (15) κaω cos ( ωτ i ) − κ d sin ( ωτ i ) = 0 . (16) Squaring and adding (15) and (16) yields ω − κ a ω − κ d = 0 . Solving for ω , we obtain ω , = κ (cid:32) a ± √ a + 4 d (cid:33) . Since we are searching for a positive root, we discard the negative root. The positive root of the above expressionis given by ω = κ (cid:115) a ( a + (cid:112) a + 4 ˜ d )2 . (17)For convenience, we write the above equation as ω = κχ. Notice that, on re-arranging (16), we obtain κ ˜ d tan( ωτ i ) = ω. Substituting for ω in the above equation and simplifying yields ω = 1 τ i tan − (cid:18) χ ˜ d (cid:19) . (18)Substituting ω in (16) and simplifying, we obtain κ cr = 1 τ i χ tan − (cid:18) χ ˜ d (cid:19) . (19)Thus, (18) and (19) yield the angular frequency of the oscillatory solution and the value of κ at which such asolution exists respectively.We now show that the MOVM undergoes a Hopf bifurcation at κ = κ cr . To that end, we need to prove thetransversality condition of the Hopf spectrum. That is, we must show that [11, Chapter 11, Theorem 1.1]Real (cid:18) d λ d κ (cid:19) κ = κ cr (cid:54) = 0 . (20)To that end, we differentiate (14) with respect to κ and simplify it, to obtainReal (cid:32)(cid:18) d λ d κ (cid:19) − (cid:33) κ = κ cr = κ cr ω τ i ( κ cr ˜ d cos( ω τ i ) + ω )( κ cr ˜ d cos( ω τ i ) + ω ) + ( κ cr ˜ d sin( ω τ i )) > . (21)The positivity in (21) follows because cos( ω τ i ) = κ cr ˜ d/ ( κ cr ˜ d + ω ) is positive. This expression follows from (16)using trigonometric manipulations. Also, Real ( z ) > if and only if Real (1 /z ) > ∀ z ∈ C . Hence, from (21) wehave Real (cid:18) d λ d κ (cid:19) κ = κ cr > . This proves the transversality of the Hopf spectrum. Therefore, the MOVM transits from the locally stable to thelocally unstable regime via a Hopf bifurcation at κ = κ cr . It can be shown that for sufficiently small values of κ, system (11) is locally stable. Additionally, the above strict inequality implies that the eigenvalues move fromleft to right in the Argand plane as κ is increased in the neighborhood of κ cr . Therefore, κ < κ cr is the necessaryand sufficient condition for local stability of system (11). B. Discussion
A few comments are in order.(1) Note that the characteristic equation (14) is transcendental, hence there exist infinitely many roots. However,system (11) loses local stability when the first conjugate pair of eigenvalues crosses the imaginary axis as Sensitivity parameter, a R eac ti ond e l a y , τ . . . SCN&SC
Fig. 1:
Stability chart:
Illustrates the necessary and sufficient condition (N&SC) (22) and the sufficient condition(SC) (8) for the MOVM, for small delays. The plot serves to validate our analysis presented in Section IV.the exogenous parameter is varied. Due to the positivity of the derivative in (21), system stability cannot berestored by increasing κ. (2) The equation of the stability boundary pertaining to system (11) is κ = κ cr . It is also called the Hopf boundaryof the said system. To obtain the Hopf boundary corresponding to the MOVM described by system (2), wetune the system parameters such that κ cr = 1 in (19). In particular, the MOVM is locally stable if and onlyif, for each i ∈ { , , · · · , N } , we have τ i < χ tan − (cid:18) χ ˜ d (cid:19) . (22)It is clear from (22) that when the reaction delay increases, the MOVM loses local stability via a Hopfbifurcation. Also note that when τ = 0 , (22) is trivially satisfied for all physically relevant parameter values.This is in agreement with the result derived in Section III. To validate the analysis presented in Section IV,we plot the Right Hand Sides (RHSs) of (8) and (22) for small values of the reaction delay in Fig. 1. Clearly,we notice from Fig. 1 that (8) indeed represents a sufficient condition for the local stability of the MOVM forsmall delays.(3) Loss of local stability via a Hopf bifurcation results in the emergence of limit cycles. Since the dynamicalvariables for the MOVM correspond to relative velocities and headways, these non-linear oscillations physicallymanifest as back-propagating congestion wave on a highway. Thus, as mentioned in the Introduction, ouranalysis provides a mathematical basis to the commonly-observed “phantom jam.”(4) Note that the non-dimensional parameter κ is not a model parameter; rather, it is an exogenous mathematicalentity introduced to aid the analysis and capture any interdependence among model parameters. It also serves tosimplify the algebra required to obtain the necessary and sufficient condition for local stability of the MOVM.Further, since substituting κ = 1 yields the MOVM, it is useful in a neighborhood around , i.e., near thestability boundary.(5) Gain parameters are known to destabilize feedback systems [1, Section 3.7]. Thus, we need to verify that the bifurcation phenomenon proved in this section is not an artefact of the exogenous parameter. To that end, weneed to verify that the MOVM also undergoes a Hopf bifurcation when one of the model parameters is chosenas the bifurcation parameter. It was shown in [38] that the transversality condition of the Hopf spectrum holdstrue for the characteristic equation of the form (14) (with κ = 1 ) when τ is used as the bifurcation parameter,although in a different context. VI. N ON - OSCILLATORY CONVERGENCE
In the previous three sections, we derived conditions for the MOVM to be locally stable in three different regimes.In the next two sections, we explore two important properties of the MOVM; namely, non-oscillatory convergenceand the rate of convergence.In the context of transportation networks, ride quality is of utmost importance. This, in turn, mandates that thevehicles avoid jerky motion. Since relative velocities and headways constitute dynamical variables for the MOVM,it boils down to studying the non-oscillatory property of its solutions. In particular, we derive the necessary andsufficient condition for non-oscillatory convergence of the MOVM. Mathematically, this amounts to ensuring thatthe eigenvalues corresponding to system (4) are negative real numbers.To derive the necessary and sufficient condition for non-oscillatory convergence of the MOVM, we begin with thecharacteristic equation corresponding to system (2), obtained by setting κ = 1 in (14). We also drop the subscript“ i ” for convenience. Thus, we obtain f ( λ ) = λ + ( aλ + d ) e − λτ = 0 . (23)To ensure non-oscillatory convergence of the MOVM, we require the roots of (23) to be real and negative. Tothat end, we substitute λ = σ + jω in (23), where j = √− . This yields aω sin( ωτ ) + ( aσ + d ) cos( ωτ ) = ( ω − σ ) e στ , and (24) aω sin( ωτ ) + ( aσ + d ) cos( ωτ ) = ( − σω ) e στ . (25)Squaring and adding (24) and (25), we obtain ( aω ) + ( aσ + d ) = ( ω + σ ) e στ . (26)To ensure that the roots are real, we require a condition for ω = 0 to be the only solution of (26). Substituting ω = 0 in (26), we obtain ( aσ + d ) = σ e στ . (27)Thus, the above condition is necessary for ω = 0 to be a solution of (26). To check whether it is also a sufficientcondition, we first separate the terms containing ω in (26) from those without it. This yields e στ ω + (2 σ e στ − a ) ω = ( aσ + d ) − σ e στ . Assuming ( aσ + d ) = σ e στ , we solve the above quadratic in ω to obtain ω = 0 or ω = a − σ e στ e στ . R eac ti ond e l a y , τ Sensitivity coefficient, a . . . . τ noc τ cr Fig. 2: Illustration of region of non-oscillatory convergence for the MOVM. Here, τ cr and τ noc represent theboundary of the locally stable region and the region of non-oscillatory convergence of the MOVM respectively.Notice the stringent requirements on the reaction delay for the solutions of the MOVM to be non-oscillatory.Thus, for ω = 0 to be the unique solution of (26), we require a = 2 σ e στ in addition to the condition mentionedin (27). That is, (23) has real eigenvalues if and only if ( aσ + d ) = σ e στ , and a = 2 σ e στ . (28)Solving the above two equations for the eigenvalue, we obtain σ = ˜ dm ± , with m ± = − ± √ . (29)Notice from the foregoing analysis that the eigenvalues are guaranteed to be negative if they are real. Substituting (29)in (23) and re-arranging, we obtain the boundary for the region of non-oscillatory convergence as e − ˜ dτm ± = − m ± ˜ da ( m + 1) . Notice that the Left Hand Side (LHS) in the above equation is a non-negative quantity. The RHS is non-negativefor m − but not for m + . Hence, we set m = m − in the above equation, and re-arrange to obtain τ noc = 1 m ˜ d ln (cid:18) − a ( m + 1) m ˜ d (cid:19) , (30)where τ noc represents the boundary for the region of non-oscillatory convergence in the τ -domain. Therefore, τ < τ noc represents the necessary and sufficient condition for non-oscillatory convergence of the MOVM. We notethat the following inequalities must be satisfied: < τ noc < τ cr , where τ cr is the RHS of (22).In summary, the necessary and sufficient condition for non-oscillatory convergence of the MOVM is τ i < m ˜ d ln (cid:18) − a ( m + 1) m ˜ d (cid:19) , (31)for each i ∈ { , , · · · , N } , when the RHS is positive and less than τ cr . We now illustrate the boundary for the region of non-oscillatory convergence of the MOVM described by (30). Inorder to better-understand the stringent constraints on system parameters to achieve non-oscillatory convergence, we also plot the necessary and sufficient condition for local stability (22) of the MOVM. To that end, we make use ofthe Bando OVF. We let the equilibrium velocity of the lead vehicle to be ˙ x = 5 m/s, and the model parameters as y ∗ = 2 m, ˜ y = 5 m and y m = 1 m. We then compute the corresponding V and ˜ d. We vary the sensitivity coefficient a from and , and compute the requisite boundaries using the scientific computation software MATLAB.Fig. 2 portrays regions of local stability and non-oscillatory convergence for the MOVM in the ( a, τ ) -space. For afixed a, the reaction delay must not exceed τ cr (respectively, τ noc ) for the MOVM to be locally stable (respectively,possess non-oscillatory solutions). Clearly, the values of τ need to be much smaller for the solutions of the MOVMto be non oscillatory as opposed to the stability of the MOVM, for a fixed value of a. In fact, as the sensitivityparameter a increases, the corresponding value of reaction delays required to ensure non-oscillatory convergencedecreases rapidly.We end this section with two remarks. ( i ) To the best of our knowledge, the analysis presented in this sectionis the first to address non-oscillatory convergence of systems with characteristic equations of the form (23) usingspectral-domain techniques. ( ii ) We can obtain ω = 0 as the only solution to (26) by a geometrical method asfollows. Re-arranging (26) yields ( ω + σ ) = ( a e − στ ) ω + ( aσ + d ) e − στ . Notice that the LHS and RHS of the above equation represent a parabola and a line in ω -domain respectively. Sincea parabola is strictly convex, the tangent to a parabola at any point will intersect it only at that point. Therefore,the slope of the line (captured by the coefficient of the term ω ) must equal the derivative of the LHS with respectto ω , evaluated at zero. That is, a e − στ = 2 σ . Also, the intercept of the line should be equal the value of theparabola evaluated at zero, i.e., ( aσ + d ) e − στ = σ . These are same as the conditions in (28).VII. R
ATE OF CONVERGENCE
In this section, we characterize the time required to attain the uniform traffic flow, once the traffic flow isperturbed (by events such as the departure of a vehicle from the platoon). Mathematically, it is related to the rate ofconvergence of solutions of the MOVM to the desired equilibrium. To that end, we follow [39] and first characterizethe rate of convergence of the MOVM. Then, using the notion of settling time, we derive an expression for thetime a platoon takes to attain the desired equilibrium following a perturbation.We begin by recalling the characteristic equation pertaining to system (4) from Section V-A. Dropping thesubscript “ i ” for ease of exposition, and setting κ = 1 in (14), we obtain λ + aλe − λτ + de − λτ = 0 . Using the change of variables z = λτ, the above equation results in z e z + a ∗ z + d ∗ = 0 , (32)where a ∗ = aτ and d ∗ = dτ . Notice that (32) has the same form as [39, Equation (22) ]. Hence, following [39],we substitute z = ψ − σ, where σ is non-negative and real, in (32) to obtain ( ψ − σψ + σ ) e ψ + a ∗ e σ ψ + ( d ∗ − a ∗ σ ) e σ = 0 . . . . R eac ti ond e l a y , τ Sensitivity coefficient, a τ cr τ noc . . . . (a) . . R eac ti ond e l a y , τ Sensitivity coefficient, a τ cr τ noc . . . . (b) Fig. 3:
Contour plots:
Contour lines of the rate of convergence overlaying the boundaries of the locally stableregion and the region of non-oscillatory convergence of the MOVM. While ( a ) is for low to medium values ofrate of convergence, ( b ) is for high values. From ( a ) , observe: ( i ) the rapid change in the gradient of the rate ofconvergence, and ( ii ) for lower values of the rate of convergence, non-oscillatory convergence is also guaranteed. Incontrast, ( b ) shows that very high rates of convergence cannot be achieved if the solutions are to be non oscillatory.The characteristic equation corresponding to the above system is obtained by substituting ψ = τ λ as λ + (cid:18) − στ (cid:19) λ + ( ae σ ) λe − λτ + (cid:16) d − aστ (cid:17) e σ e − λτ + (cid:18) σ τ (cid:19) = 0 . (33)The rate of convergence is the largest σ ≥ such that the root of (33) with the largest real part is negative [39].As pointed out in [39], finding such a σ analytically is intractable. Hence, we illustrate the variation of the rateof convergence numerically with respect to both the sensitivity parameter a and the reaction delay τ , using thescientific computation software MATLAB.We consider the Bando OVF, and set the following parameters: y m = 1 m, ˜ y = 5 m, y ∗ = 2 m and ˙ x = 5 m/s.We then compute the corresponding values of V and ˜ d. Next, we vary the sensitivity coefficient a in the range [1 , , and for each of its values, we compute the critical value of the reaction delay τ cr using (22). We then vary thereaction delay τ in the range [0 , τ cr ] , for each a. For every pair ( a, τ ) in this range, σ is increased from , till theroot of (33) with the largest real part crosses the imaginary axis in the Argand plane. Since the resulting plot wouldbe three dimensional, we present the corresponding contour plots in Fig. 3. For clarity in presentation, the contourplots are segregated as follows: Fig. 3a is for low to medium values of the rate of convergence, whereas Fig. 3b isfor high values. It can be seen from Fig. 3a that small changes in a or τ causes the rate of convergence to changefrom . to . . However, it would require relatively larger changes in a or τ for the rate of convergence to changefrom . to . . That is, the gradient of the rate of convergence increases rather rapidly with an increase in the rateof convergence. Also, for low values of the rate of convergence, non-oscillatory convergence can be guaranteed. Incontrast, Fig. 3b brings forth the trade-off between the rate of convergence and non-oscillatory convergence; very high rates of convergence cannot be achieved if the solutions are to be non oscillatory.The rate of convergence determines the time taken by a platoon to reach an equilibrium (denoted by T eMOV M ).To characterize T eMOV M , we first define the time taken by the i th pair of vehicles in the platoon following thestandard control-theoretic notion of “settling time.” That is, by t ei ( (cid:15) ) , we denote the minimum time taken by thetime-domain trajectory of the MOVM to enter and subsequently remain within the (cid:15) -band around the equilibrium.For simplicity, we drop the explicit dependence on (cid:15). Then, T eMOV M = N (cid:88) i =1 t ei . (34)It is clear that (34) is an upper bound on the time taken by the platoon to equilibrate. However, the equality holdssince the i th pair cannot equilibrate till the ( i − th pair has reached its equilibrium.VIII. H OPF BIFURCATION ANALYSIS
In the previous sections, we have characterized the stable region for the MOVM, and studied two of its mostimportant properties; namely, non-oscillatory convergence and the rate of convergence. We have also proved thatsystem (2) loses stability via a Hopf bifurcation, thus resulting in limit cycles. In this section, we provide ananalytical framework to characterize the type of the bifurcation and the asymptotic orbital stability of the emergentlimit cycles. We closely follow the style of analysis presented in [36], which uses Poincar´e normal forms and thecenter manifold theory.We begin by denoting the RHS of (11) as f i . That is, for i ∈ { , , · · · , N } ,f i (cid:44) aκ ( V ( y i − ( t − τ i − )) − V ( y i ( t − τ i )) − v k ( t − τ i )) . (35)Define µ = κ − κ cr . Notice that the system undergoes a Hopf bifurcation at µ = 0 , where κ = κ cr . Henceforth,we shall consider µ as the bifurcation parameter. Also, it is clear that when µ > , the exogenous parameter κ changes from κ cr to κ cr + µ, thus pushing the system into its unstable regime.We now provide a step-by-step overview of the detailed local bifurcation analysis, before delving into its technicaldetails. Step 1 : We begin by applying Taylor’s series expansion to the RHS of (35). Next, we separate the linear termsfrom their non-linear counterparts. This allows us to cast the resulting equation into the standard form of an OperatorDifferential Equation (OpDE).
Step 2 : When µ = 0 , the system has exactly one pair of purely imaginary eigenvalues with non-zero angularvelocity, as seen from (21). We call the linear space spanned by the corresponding eigenvectors as the criticaleigenspace. For the purpose of our analysis, we also require a locally invariant manifold that is a tangent to thecritical eigenspace at the system’s equilibrium. The center manifold theorem [36] guarantees the existence of sucha manifold. Step 3 : Next, we project the system onto its critical eigenspace and its complement when µ = 0 . Thus, we maywrite the dynamics of the original system on the center manifold as an ODE in a single complex variable. Step 4 : Finally, using Poincar´e normal forms, we evaluate the Lyapunov coefficient and the Floquet exponent.These, in turn, help characterize the type of the Hopf bifurcation and the asymptotic orbital stability of the emergentlimit cycles.We begin the analysis by expanding (35) about the equilibrium v ∗ i = 0 , y ∗ i = V (cid:48) ( V − ( ˙ x )) , i = 1 , , · · · , N, using Taylor’s series, to obtain ˙ v i ( t ) =( − κa ) v i,t ( − τ i ) + Ω (1) i y i,t ( − τ i ) + Ω (2) i y i,t ( − τ i ) + Ω (3) i y i,t ( − τ i ) + ζ (1) i y ( i − ,t ( − τ i − )+ ζ (2) i y i − ,t ( − τ i − ) + ζ (3) i y i − ,t ( − τ i − ) + higher order terms , ˙ y i ( t ) = κx i ( t ) (36)where we use the shorthand v i,t ( − τ i ) and y i,t ( − τ i ) to represent v i ( t − τ i ) and y i ( t − τ i ) respectively. The coefficientsare given by Ω (1) i = − κaV (cid:48) ( y ∗ i ) , Ω (2) i = − κaV (cid:48)(cid:48) ( y ∗ i ) , and Ω (3) i = − κaV (cid:48)(cid:48)(cid:48) ( y ∗ i ) . Also, the coefficients ζ (1) i , ζ (2) i and ζ (3) i assume the values Ω i (1) , Ω i (2) and Ω i (3) respectively for i > , and are zero for i = 1 . The termswith the primes V (cid:48) , V (cid:48)(cid:48) and V (cid:48)(cid:48)(cid:48) denote the first, second and third derivatives of the OVF with respect to the statevariable respectively.In the following, we use C k ( A ; B ) to denote the linear space of all functions from A to B which are k timesdifferentiable, with each derivative being continuous. Also, we use C to denote C , for convenience.With the concatenated state S ( t ) , note that (11) is of the form:d S ( t ) d t = L µ S t ( θ ) + F ( S t ( θ ) , µ ) , (37)where t > , µ ∈ R , and where for τ = max i τ i > , S t ( θ ) = S ( t + θ ) , S : [ − τ, −→ R N , θ ∈ [ − τ, . Here, L µ : C (cid:0) [ − τ, R N (cid:1) −→ R N is a one-parameter family of continuous, bounded linear functionals, whereasthe operator F : C (cid:0) [ − τ, R N (cid:1) −→ R N is an aggregation of the non-linear terms. Further, we assume that F ( S t , µ ) is analytic, and that F and L µ depend analytically on the bifurcation parameter µ , for small | µ | . Theobjective now is to cast (37) in the standard form of an OpDE:d S t d t = A ( µ ) S t + R S t , (38)since the dependence here is on S t alone rather than both S t and S ( t ) . To that end, we begin by transforming thelinear problem d S ( t ) / d t = L µ S t ( θ ) . We note that, by the Riesz representation theorem [40, Theorem 6.19], thereexists a N × N matrix-valued measure η ( · , µ ) : B (cid:0) C (cid:0) [ − τ, R N (cid:1)(cid:1) −→ R N × N , wherein each component of η ( · ) has bounded variation, and for all φ ∈ C (cid:0) [ − τ, R N (cid:1) , we have L µ φ = (cid:90) − τ d η ( θ, µ ) φ ( θ ) . (39)In particular, L µ S t = (cid:90) − τ d η ( θ, µ ) S ( t + θ ) . Motivated by the linearized system (12), we defined η = ˜ A ˜ B ˜ C ˜ D d θ, where ( ˜ A ) ij = − κdδ ( θ + τ i ) , i = j,κdδ ( θ + τ j ) , j = i − , i > , , otherwise, ( ˜ B ) ij = − κaδ ( θ + τ i ) , i = j, , otherwise, ˜ C = κI N × N and ˜ D = 0 N × N . For instance, when N = 2,d η = − κdδ ( θ + τ ) 0 − κaδ ( θ + τ ) 0 κdδ ( θ + τ ) − κdδ ( θ + τ ) 0 − κaδ ( θ + τ ) κ κ d θ. For φ ∈ C (cid:0) [ − τ, C N (cid:1) , we define A ( µ ) φ ( θ ) = d φ ( θ ) d θ , θ ∈ [ − τ, , (cid:82) − τ d η ( s, µ ) φ ( s ) ≡ L µ , θ = 0 , (40)and R φ ( θ ) = , θ ∈ [ − τ, , F ( φ, µ ) , θ = 0 . With the above definitions, we observe that d S t / d θ ≡ d S t / d t. Hence, we have successfully cast (37) in the formof (38). To obtain the required coefficients, it is sufficient to evaluate various expressions for µ = 0 , which weuse henceforth. We start by finding the eigenvector of the operator A (0) with eigenvalue λ (0) = jω . That is, wewant an N × vector (to be denoted by q ( θ ) ) with the property that A (0) q ( θ ) = jω q ( θ ) . We assume the form: q ( θ ) = [1 φ φ · · · φ N − ] T e jω θ , and solve the eigenvalue equations. That is, we need to solve − κde − jω τ − aφ N e − jω τ κde − jω τ + κ ( − dφ − aφ N +1 ) e − jω τ κde − jω τ + κ ( − dφ − aφ N +2 ) e − jω τ ... κde − jω τ N − + κ ( − dφ N − − aφ N − ) e − jω τ N − κβκφ β ... κφ N − β = jω φ φ ... φ N − φ N − , where β = j ( − je − jω τ ) /ω . This, in turn, necessitates the following assumptions:(i) − jω e jω τ + κdκ a = − e jω τ ω , (ii) For each i ∈ { , , · · · N − } , the following matrix is invertible: κde − jω τ i +1 + jω κae − jω τ i +1 κβ − jω , where β = j ( − je − jω τ ) /ω . Then, for i ∈ { , , · · · N − } ,φ N = κβjω , φ i = − jκω de − jω τ i ∆ M i , and φ N + i = − βκ de − jω τ i ∆ M i , where ∆ M i = ω − κdω sin( ω τ i +1 ) − κ βa cos( ω τ i +1 ) + j (cid:0) κ βa sin( ω τ i +1 ) − κdω cos( ω τ i +1 ) (cid:1) . We define the adjoint operator as follows: A ∗ (0) φ ( θ ) = − d φ ( θ ) d θ , θ ∈ (0 , τ ] , (cid:82) − τ d η T ( s, φ ( − s ) , θ = 0 , where d η T is the transpose of d η. We note that the domains of A and A ∗ are C (cid:0) [ − τ, C N (cid:1) and C (cid:0) [0 , τ ]; C N (cid:1) respectively. Therefore, if jω is an eigenvalue of A , then − jω is an eigenvalue of A ∗ . Hence, to find the eigenvectorof A ∗ (0) corresponding to − jω , we assume the form: p ( θ ) = B [ ψ N − ψ N − ψ N − · · · T e jω θ , and solve A ∗ (0) p ( θ ) = − jω p ( θ ) . Simplifying this, we obtain − κde − jω τ ψ N − + κde − jω τ ψ N − − βκψ N − − κde − jω τ ψ N − + κde − jω τ ψ N − − βκψ N − − κde − jω τ ψ N − + κde − jω τ ψ N − − βκψ N − ... − κde − jω τ N ψ N − βκ − κae − jω τ − κae − jω τ ... − κae − jω τ N = − jω ψ N − ψ N − ψ n − ... ψ . We require the following assumptions:(i) κd − jω e jω τ N κ a = − e jω τ ω , (ii) For each i ∈ { , , · · · N − } , the following matrix is invertible: jω − κae − jω τ i κβ κde − jω τ i − jω . Then, for i ∈ { , , · · · N − } , we obtain ψ N = jω e jω τ N κa , ψ N + i = jω κdψ N + i − e − jω τ N − i ∆ ˜ M i , and ψ i = κ adψ N + i − e − jω τ N − i ∆ ˜ M i , where ∆ ˜ M i = ω + κdω sin( ω τ N − i ) + κ βa cos( ω τ N − i ) + j (cid:0) κdω cos( ω τ N − i ) − κ βa sin( ω τ N − i ) (cid:1) . The normalization condition for Hopf bifurcation requires that (cid:104) p, q (cid:105) = , thus yielding an expression for B. For any q ∈ C (cid:0) [ − τ, C N (cid:1) and p ∈ C (cid:0) [0 , τ ]; C N (cid:1) , the inner product is defined as (cid:104) p, q (cid:105) (cid:44) ¯ p · q − (cid:90) θ = − τ θ (cid:90) ζ =0 ¯ p T ( ζ − θ ) d ηq ( ζ ) d ζ, (41)where the overbar represents the complex conjugate and the “ · ” represents the regular dot product. The value of B such that the inner product between the eigenvectors of A and A ∗ is unity can be shown to be B = 1 ζ + ζ + ζ + ζ , where ζ = (cid:18) e jω τ − e j ω τ − (cid:19) N − (cid:88) i =0 κψ N − i − ¯ φ i , ζ = N − (cid:88) i =0 (cid:18) e jω τ i +1 − e j ω τ i +1 jω (cid:19) κψ N − − i ( a ¯ φ i + d ¯ φ N + i ) ,ζ = N − (cid:88) i =0 (cid:18) e j ω τ i +1 − e jω τ i +1 jω (cid:19) κd ¯ φ i ψ N − − i , and, ζ = N − (cid:88) i =0 ψ N − − i ¯ φ i . In the above, we define φ = ψ = 0 for notational brevity.For S t , a solution of (38) at µ = 0 , we define z ( t ) = (cid:104) p ( θ ) , S t (cid:105) , and w ( t, θ ) = S t ( θ ) − Real ( z ( t ) q ( θ )) . Then, on the center manifold C , we have w ( t, θ ) = w ( z ( t ) , ¯ z ( t ) , θ ) , where w ( z ( t ) , ¯ z ( t ) , θ ) = w ( θ ) z w ( θ ) ¯ z w ( θ ) z ¯ z + · · · . (42)Effectively, z and ¯ z are the local coordinates for C in C in the directions of p and ¯ p respectively. We note that w is real if S t is real, and we deal only with real solutions. The existence of the center manifold C enables thereduction of (38) to an ODE in a single complex variable on C . At µ = 0, the said ODE can be described as ˙ z ( t ) = (cid:104) p, A S t + R S t (cid:105) , = jω z ( t ) + ¯ p (0) . F ( w ( z, ¯ z, θ ) + 2 Real ( z ( t ) q ( θ ))) , = jω z ( t ) + ¯ p (0) . F ( z, ¯ z ) . (43)This is written in abbreviated form as ˙ z ( t ) = jω z ( t ) + g ( z, ¯ z ) . (44)The objective now is to expand g in powers of z and ¯ z . However, this requires w ij ( θ ) ’s from (42). Once these areevaluated, the ODE (43) for z would be explicit (as given by (44)), where g can be expanded in terms of z and ¯ z as g ( z, ¯ z ) = ¯ p (0) . F ( z, ¯ z ) = g z g ¯ z g z ¯ z + g z ¯ z · · · . (45) Next, we write ˙ w = ˙ S t − ˙ zq − ˙¯ z ¯ q. Using (38) and (44), we then obtain the following ODE: ˙ w = A w − Real (¯ p (0) . F q ( θ )) , θ ∈ [ − τ, , A w − Real (¯ p (0) . F q (0)) + F , θ = 0 . This can be re-written using (42) as ˙ w = A w + H ( z, ¯ z, θ ) , (46)where H can be expanded as H ( z, ¯ z, θ ) = H ( θ ) z H ( θ ) ¯ z H ( θ ) z ¯ z + H ( θ ) z ¯ z · · · . (47)Near the origin, on the manifold C , we have ˙ w = w z ˙ z + w ¯ z ˙¯ z. Using (42) and (44) to replace w z ˙ z (and theirconjugates, by their power series expansion) and equating with (46), we obtain the following operator equations: (2 jω − A ) w ( θ ) = H ( θ ) , (48) −A w = H ( θ ) , (49) − (2 jω + A ) w ( θ ) = H ( θ ) . (50)We start by observing that S t ( θ ) = w ( θ ) z w ( θ ) ¯ z w ( θ ) z ¯ z + zq ( θ ) + ¯ z ¯ q ( θ ) + · · · . From the Hopf bifurcation analysis [36], we know that the coefficients of z , ¯ z , z ¯ z , and z ¯ z terms are used toapproximate the system dynamics. Hence, we only retain these terms in the expansions.To obtain the effect of non-linearities, we substitute the aforementioned terms appropriately in the non-linearterms of (36) and separate the terms as required. Therefore, for each i ∈ { , , . . . , N } , we have the non-linearityterm to be F i = F i z F i ¯ z F i z ¯ z + F i z ¯ z , (51)where, for i ∈ { , , · · · , N } , the coefficients are given by F i = Ω (1) i w i ( − τ i ) + ζ (1) i − w i − ( − τ i − ) , F i = Ω (1) i w i ( − τ i ) + ζ (1) i − w i − ( − τ i − ) , F i = Ω (1) i w i ( − τ i ) + ζ (1) i − w i − ( − τ i − ) , F i = 2Ω (2) i (cid:0) w i ( − τ i ) e jω τ i + 2 w i ( − τ i ) e − jω τ i (cid:1) + 2 ζ (2) i − (cid:0) w i − ( − τ i − ) e jω τ i − + 2 w i − ( − τ i − ) e − jω τ i − (cid:1) . For i ∈ { N + 1 , N + 2 , · · · , N } , each of these coefficients are zero. This is so, since last N states correspond tothe headways who evolution equations are linear. We represent the vector of non-linearities used in (43) as F = [ F F · · · F N ] T . Next, we compute g using F as g ( z, ¯ z ) = ¯ p (0) . F = ¯ B N (cid:88) l =1 ¯ ψ N − l F l . (52)Substituting (51) in (52), and comparing with (45), we obtain g x = ¯ B N (cid:88) l =1 ¯ ψ N − l F xl , (53)where x ∈ { , , , } . Using (53), the corresponding coefficients can be computed. However, computing g requires w ( θ ) and w ( θ ) . Hence, we perform the requisite computation next. For θ ∈ [ − τ, , H can be simplifiedas H ( z, ¯ z, θ ) = − Real (¯ p (0) . F q ( θ )) , = − (cid:18) g z g ¯ z g z ¯ z + · · · (cid:19) q ( θ ) − (cid:18) ¯ g ¯ z g z g z ¯ z + · · · (cid:19) ¯ q ( θ ) , which, when compared with (47), yields H ( θ ) = − g q ( θ ) − ¯ g ¯ q ( θ ) , (54) H ( θ ) = − g q ( θ ) − ¯ g ¯ q ( θ ) . (55)From (40), (48) and (49), we obtain the following ODEs: ˙ w ( θ ) = 2 jω w ( θ ) + g q ( θ ) + ¯ g ¯ q ( θ ) , (56) ˙ w ( θ ) = g q ( θ ) + ¯ g ¯ q ( θ ) . (57)Solving (56) and (57), we obtain w ( θ ) = − g jω q (0) e jω θ − ¯ g jω ¯ q (0) e − jω θ + e e jωθ , (58) w ( θ ) = g jω q (0) e jω θ − ¯ g jω ¯ q (0) e − jω θ + f , (59)for some vectors e and f , to be determined.To that end, we begin by defining the following vector: ˜ F (cid:44) [ F F · · · F N ] T . Equating (48) and (54),and simplifying, yields the operator equation: jω e − A (cid:0) e e jω θ (cid:1) = ˜ F . To solve this, we assume that thefollowing matrices to be invertible for each i ∈ { , , · · · , N } , jω o + κ ( a + d ) e − jω τ i κ ( a + d ) e − jω τ i − κτ jω o . Then, the operator equation can be simplified to (cid:0) jω + κ ( a + d ) e − jω τ (cid:1) e + κ ( a + d ) e − jω τ e N +1 (cid:0) jω + κ ( a + d ) e − jω τ (cid:1) e + κ ( a + d ) e − jω τ e N +2 ... (cid:0) jω + κ ( a + d ) e − jω τ N (cid:1) e N + κ ( a + d ) e − jω τ N e N − κτ e + 2 jω e N +1 ... − κτ e N + 2 jω e N = ˜ F . This, in turn, yields e i = 2 jω F i ∆ M ∗ i , and, e N + i = κτ F i ∆ M ∗ i , (60)for i ∈ { , , · · · , N } . Here, ∆ M ∗ i = − ω + 2 ω κ ( a + d ) sin( ω τ i ) + τ κ ( a + d ) cos( ω τ i ) + j (cid:0) ω κ ( a + d ) cos( ω τ i ) − τ κ ( a + d ) sin( ω τ i ) (cid:1) . Next, equating (49) and (55), and simplifying, we obtain the operator equation A f = − ˜ F , with ˜ F (cid:44) [ F F · · · F N ] T . On simplification, this equation yields − κτ ( a + d )( f + f N +1 ) − κτ ( a + d )( f + f N +2 ) ... − κτ ( a + d )( f N + f N ) κτ f ... κτ f N = − ˜ F . On solving this, we obtain for i ∈ { , , · · · , N } , f i = 0 , and, f N + i = F i κτ i ( a + d ) . (61)Substituting for e and f from (60) and (61) in (58) and (59) respectively, we obtain w ( θ ) and w ( θ ) . This, inturn, facilitates the computation of g . We can then compute c (0) = j ω (cid:18) g g − | g | − | g | (cid:19) + g ,α (cid:48) (0) = Real (cid:18) d λ d κ (cid:19) κ = κ cr , µ = − Real ( c (0)) α (cid:48) (0) , and β = 2 Real ( c (0)) . Here, c (0) is known as the Lyapunov coefficient and β is the Floquet exponent. It is known from [36] thatthese quantities are useful since ( i ) If µ > , then the bifurcation is supercritical , whereas if µ < , then the bifurcation is subcritical . ( ii ) If β > , then the limit cycle is asymptotically orbitally unstable , whereas if β < , then the limit cycle is asymptotically orbitally stable .We now present numerically-constructed bifurcation diagrams to gain some insight into the effect of variousparameters on the amplitude of the limit cycle. . . A m p lit ud e (r e l a ti v e v e l o c it y ) Bifurcation parameter, k y ∗ i = y ∗ i = y ∗ i = (a) . . A m p lit ud e (r e l a ti v e v e l o c it y ) Bifurcation parameter, k y ∗ i = y ∗ i = y ∗ i = (b) Fig. 4:
Bifurcation diagram:
Amplitude of the emergent limit cycles in relative velocity variable as a function of theexogenous parameter κ. ( a ) is for the Bando OVF, while ( b ) is for the Underwood OVF. For a fixed κ ∈ [1 , . , the Underwood OVF results in limit cycles of smaller relative velocity than its Bando counterpart. Bifurcation diagrams
To obtain bifurcation diagrams, we make use of DDE-BIFTOOL [41], [42]. We first input system (11) and theirfirst-order derivatives with respect to the state and delayed state variables to DDE-BIFTOOL. We then set κ = 1 and initialize the model parameters appropriately. We also fix a range of variation for the bifurcation parameter.DDE-BIFTOOL varies the bifurcation parameter accordingly and finds its critical value. We then increase the valueof κ and record the amplitude of the resulting limit cycle, thus obtaining the bifurcation diagram. We use the SIunits throughout. That is, time will be expressed in “seconds,” distance in “meters,” velocity in “meters per second”and the sensitivity coefficient in “inverse second.” For our comparison, we consider two optimal velocity functions;namely, the Bando OVF and the Underwood OVF.For the Bando OVF, we initialize the parameters as follows: N = 4 , a = 1 . , τ = 0 . , τ = 0 . , τ = 0 . and τ = 0 . . We fix y m = 2 and ˜ y = 5 , and compute V for each of y ∗ i = 1 , and . The vehicle indexed isconsidered to undergo a Hopf bifurcation. For the case of the Underwood OVF, we set the following values forthe parameters. N = 3 , a = 1 . , τ = 0 . , τ = 0 . and τ = 0 . . We fix y m = 2 , and compute V foreach of y ∗ i = 1 , and . The vehicle indexed is then considered to undergo a Hopf bifurcation. We choose theequilibrium velocity of the lead vehicle, ˙ x = 5 . The bifurcation diagrams are shown in Figure 4. As seen from the figure, the amplitude of the relative velocityincreases with an increase in κ. However, for a fixed value of the exogenous parameter, the Underwood OVFyields limit cycles with smaller relative velocity that its Bando counterpart, which is desirable. Also, notice thatthe amplitude of the emergent limit cycles increases with an increase in the equilibrium headway. This is intuitivebecause larger equilibrium headways offer more space for the resulting limit cycles to oscillate in. ˜ v , ˜ y ( × − ) Time (in seconds) −
20 30 40˜ v ( t ) ˜ y ( t ) (a) ˜ v , ˜ y ( × − ) Time (in seconds)0 −
102 64 8 ˜ v ( t ) ˜ y ( t ) (b) Fig. 5:
Simulation results:
Shows the variations in relative velocity and headway around their respective equilibria. ( a ) portrays the limit cycles predicted by (21), while ( b ) presents an instance of non-oscillatory behavior whenparameters are chosen appropriately satisfying (31).IX. S IMULATIONS
Thus far, we have analyzed the MOVM in no-delay, small-delay and arbitrary-delay regimes. We also studied twoof its important properties – non-oscillatory convergence and the rate of convergence. In the previous section, wepresented an analytical framework to characterize the type of Hopf bifurcation and the asymptotic orbital stabilityof the limit cycles that emerge when the stability conditions are marginally violated.In this section, we present the simulation results of the MOVM that serve to corroborate our analytical findings.We make use of the scientific computation software MATLAB to implement a discrete version of system (2), thussimulating the MOVM. We use T s = 10 − s as the update time. Throughout, we use SI units.To corroborate the insight from Section V, we make use of the Bando OVF. We set the parameters values as: N = 4 , a = 1 . , ˜ y = 5 , y m = 1 and y ∗ i = 3 for i = 1 , , , . We then compute the corresponding value of V using the functional form for the Bando OVF and τ cr using (22). Further, we set τ = τ cr / , τ = τ cr / ,τ = τ cr and τ = τ cr / . This ensures that the vehicle indexed ˙ x = 3 undergoes a Hopf bifurcation. The leader’svelocity profile is considered to be − e − t ) , thus yielding an equilibrium velocity of . We plot the variationof the relative velocity and the headway about their respective equilibria for the vehicle indexed . That is, we plot ˜ v ( t ) = v ( t ) − v ∗ = v ( t ) and ˜ y ( t ) = y ( t ) − y ∗ = y ( t ) − . Fig. 5a shows the emergence of limit cycles, aspredicted by the transversality condition of the Hopf spectrum (21).Next, we corroborate the analysis presented in Section VI using the Bando OVF. The parameters are initializedas: N = 4 , a = 2 , ˜ y = 25 , y m = 15 and y ∗ i = 15 for i = 1 , , , . We then compute the corresponding value of V using the functional form for the Bando OVF. We fix the equilibrium velocity of the leader to be ˙ x = 25 bysetting the its velocity profile as − e − t ) . We then compute τ noc using (30). We set the reaction delays as τ = τ noc / , τ = τ noc / , τ = τ noc / and τ = τ noc / . Fig. 5b shows an instance of the relative velocity and headway variations around their respective equilibria for the vehicle indexed . The headway and relative velocitiespossess the non-oscillatory behavior, as predicted by the analysis in Section VI.X. C
ONCLUDING R EMARKS
In this paper, we highlighted the importance of delayed feedback in determining the qualitative dynamicalproperties of a platoon of vehicles traveling on a straight road. Specifically, we analyzed the Modified OptimalVelocity Model (MOVM) in three regimes – no delay, small delay and arbitrary delay. We proved that, in theabsence of reaction delays, the MOVM is locally stable for all practically relevant values of model parameters. Wethen obtained a sufficient condition for the local stability of the MOVM by analyzing it in the small-delay regime.We also characterized the local stability region of the MOVM in the arbitrary-delay regime.We then proved that the MOVM undergoes a loss of local stability via a Hopf bifurcation. The resulting limitcycles physically manifest as a back-propagating congestion wave. Thus, our work provides a mathematical basisto explain the observed “phantom jams.” For the said analysis, we used an exogenous parameter that captures anyinterdependence among the model parameters.We then derived the necessary and sufficient condition for non-oscillatory convergence of the MOVM, with theaim of avoiding jerky vehicular motions. This, in turn, guarantees smooth traffic flow and improves ride quality. Next,we characterized the rate of convergence of the MOVM, which affects the time taken by a platoon to equilibrate.We also brought forth the trade-off between the rate of convergence and non-oscillatory convergence of the MOVM.Finally, we provided an analytical framework to characterize the type of Hopf bifurcation and the asymptoticorbital stability of the limit cycles which emerge when the stability conditions are violated. Therein, we madeuse of Poincar´e normal forms and the center manifold theory. We corroborated our analyses using stability charts,bifurcation diagrams, numerical computations and simulations conducted using MATLAB.
Avenues for further research
There are numerous avenues that merit further investigation. In this work, we have derived the conditions forpairwise stability of vehicles in a platoon, whose dynamics are captured by the MOVM. However, the string stabilityof such a platoon remains to be studied.From a practical standpoint, the parameters of the MOVM may vary, for varied reasons. Hence, it becomes im-perative that the longitudinal control algorithm be robust to such parameter variations, and to unmodeled dynamics.A
CKNOWLEDGEMENTS
This work is undertaken as a part of an Information Technology Research Academy (ITRA), Media Lab Asia,project titled “De-congesting India’s transportation networks.” The authors are also thankful to Debayani Ghosh,Rakshith Jagannath and Sreelakshmi Manjunath for many helpful discussions. R EFERENCES[1] R. Rajamani, “Vehicle Dynamics and Control,”
Springer , Second Edition, 2012.[2] A. Vahidi and A. Eskandarian, “Research advances in intelligent collision avoidance and adaptive cruise control,”
IEEE Transactions onIntelligent Transportation Systems , vol. 4, pp. 143-153 , 2003.[3] S. Greengard, “Smart transportation networks drive gains,”
Communications of the ACM , vol. 58, pp. 25-27, 2015.[4] V.A.C. van den Berg and E.T. Verhoef, “Autonomous cars and dynamic bottleneck congestions: The effects on capacity, value of time andpreference heterogeneity,”
Transportation Research Part B , vol. 94, pp. 43-60 , 2016.[5] D.C. Gazis, R. Herman and R.W. Rothery, “Nonlinear follow-the-leader models of traffic flow,”
Operations Research , vol. 9, pp. 545-567,1961.[6] M. Bando, K. Hasebe, K. Nakanishi and A. Nakayama, “Analysis of optimal velocity model with explicit delay,”
Physical Review E , vol.58, pp. 5429-5435, 1998.[7] D. Chowdhury, L. Santen and A. Schadschneider, “Statistical physics of vehicular traffic and some related systems,”
Physical Reports ,vol. 329, pp. 199-329, 2000.[8] D. Helbing, “Traffic and related self-driven many-particle systems,”
Reviews of Modern Physics , vol. 73, pp. 1067-1141, 2001.[9] G. Orosz and G. St ´ ep ´ an, “Subcritical Hopf bifurcations in a car-following model with reaction-time delay,” Proceedings of the RoyalSociety A , vol. 642, pp. 2643-2670, 2006.[10] G.K. Kamath, K. Jagannathan and G. Raina, “Car-following models with delayed feedback: local stability and Hopf bifurcation,” in
Proceedings of the rd Annual Allerton Conference on Communication, Control and Computing , 2015.[11] J.K. Hale and S.M.V. Lunel, “Introduction to Functional Differential Equations,”
Springer-Verlag , 2011.[12] R. Sipahi and S.I. Niculescu, “Analytical stability study of a deterministic car following model under multiple delay interactions,” in
Proceedings of Mechanical and Industrial Engineering Faculty Publications , 2006.[13] A. Kesting and M. Treiber, “How reaction time, update time, and adaptation time influence the stability of traffic flow,”
Computer-AidedCivil and Infrastructure Engineering , vol. 23, pp. 125-137, 2008.[14] M. Bando, K. Hasebe, A. Nakayama, A. Shibata and Y. Sukiyama, “Dynamical model of traffic congestion and numerical simulation,”
Physical Review E , vol. 51, pp. 1035-1042, 1995.[15] L.C. Davis, “Comment on analysis of optimal velocity model with explicit delay,”
Physical Review E , vol. 66, pp. 038101-2, 2002.[16] L.C. Davis, “Modifications of the optimal velocity traffic model to include delay due to driver reaction time,”
Physica A , vol. 319, pp.557-567, 2003.[17] I. Gasser, G. Sirito and B. Werner, “Bifurcation analysis of a class of ‘car following’ traffic models,”
Physica D , vol. 197, pp. 222-241,2004.[18] G. Orosz, B. Krauskopf and R.E. Wilson, “Bifurcations and multiple traffic jams in a car-following model with reaction-time delay,”
Physica D , vol. 211, pp. 277-293, 2005.[19] G. Orosz, R.E. Wilson and G. St ´ ep ´ an, (Editors), “Traffic jams: dynamics and control,” Philosophical Transactions A , vol. 368, 2010.[20] R.E. Wilson and J.A. Ward, “Car-following models: fifty years of linear stability analysis - a mathematical perspective,”
TransportationPlanning and Technology , vol. 34, pp. 3-18, 2011.[21] R. Rajamani and C. Zhu, “Semi-autonomous adaptive cruise control systems,”
IEEE Transactions on Vehicular Technology , vol. 51, pp.1186-1192, 2002.[22] Z. Qu, J. Wang and R.A. Hull, “Cooperative control of dynamical systems with application to autonomous vehicles,”
IEEE Transactionson Automatic Control , vol. 53, pp. 894-911, 2008.[23] R.U. Chavan, M. Belur, D. Chakraborty and D. Manjunath, “On the stability and formations in ad hoc multilane vehicular traffic,” in
Proceedings of the th International Conference on Communication Systems and Networks (COMSNETS) , 2015.[24] T.H. Summers, C. Yu, S. Dasgupta and B.D.O. Anderson, “Control of minimally persistent leader-remote-follower and coleader formationsin the plane,”
IEEE Transactions on Automatic Control , vol. 56, pp. 2778-2792 , 2011.[25] K.C. Dey, L. Yan, X. Wang, Y. Wang, H. Shen, M. Chowdhury, L. Yu, C. Qiu and V. Soundararaj, “A review of communication, drivercharacteristics, and controls aspects of cooperative adaptive cruise control (CACC),”
IEEE Transactions on Intelligent TransportationSystems , vol. 17, pp. 491-509, 2016.[26] M. Won, T. Park and S.H. Son, “Toward mitigating phantom jam using vehicle-to-vehicle communication,”
IEEE Transactions on IntelligentTransportation Systems , Early Access, 2016. [27] D. Chen, J. Laval, Z. Zheng and S. Ahn, “A behavioral car-following model that captures traffic oscillations,” Transportation ResearchPart B , vol. 46, pp. 744-761, 2012.[28] R. Nishi, A. Tomoeda, K. Shimura and K. Nishinari, “Theory of jam-absorbing driving,”
Transportation Research Part B , vol. 50, pp.116-129, 2013.[29] Y. Taniguchi, R. Nishi, T. Ezaki and K. Nishinari, “Jam-absorption driving with a car-following model,”
Physica A , vol. 433, pp. 304-315,2015.[30] G. Orosz, “Connected cruise control: modelling, delay effects, and nonlinear behaviour,”
Vehicle System Dynamics , vol. 54, pp. 1147-1176,2016.[31] M. di Bernardo, A. Salvi and S. Santini, “Distributed consensus strategy for platooning of vehicles in the presence of time-varyingheterogeneous communication delays,”
IEEE Transactions on Intelligent Transportation Systems , vol. 16, pp. 102-112, 2015.[32] M. Batista and E. Twrdy, “Optimal velocity functions for car-following models,”
Journal of Zhejiang University - Science A (AppliedPhysics & Engineering) , vol. 11, pp. 520-529, 2010.[33] G.K. Kamath, K. Jagannathan and G. Raina, “A computational study of a variant of the Optimal Velocity Model with no collisions,” in
Proceedings of th International Conference on Communication Systems and Networks (COMSNETS) , 2016.[34] I. Gy ¨ ori and G. Ladas, “Oscillation Theory of Delay Differential Equations With Applications,” Clarendon Press , 1991.[35] J.R. Silvester, “Determinants of block matrices,”
The Mathematical Gazette , vol. 84, pp. 460-467, 2000.[36] B.D. Hassard, N.D. Kazarinoff and Y.-H. Wan, “Theory and Applications of Hopf Bifurcation.,”
Cambridge University Press , 1981.[37] G. Raina, “Local bifurcation analysis of some dual congestion control algorithms,”
IEEE Transactions on Automatic Control , vol. 50, pp.1135-1146, 2005.[38] S. Manjunath and G. Raina, “FAST TCP: some queueing models and stability,” in
Proceedings of International Conference on SignalProcessing and Communications (SPCOM) , 2014.[39] S. Chong, S. Lee and S. Kang, “A simple, scalable, and stable explicit rate allocation algorithm for max − min flow control with minimumrate guarantee,” IEEE/ACM Transactions on Networking , vol. 9, pp. 322-335, 2001.[40] W. Rudin, “Real & Complex Analysis,”
Tata McGraw Hill Publications , Third Edition, 1987.[41] K. Engelborghs, T. Luzyanina and D. Roose, “Numerical bifurcation analysis of delay differential equations DDE-BIFTOOL,”
ACMTransactions on Mathematical Software (TOMS) , vol. 28, pp. 1-21, 2002.[42] K. Engelborghs, T. Luzyanina and G. Samaey, “DDE-BIFTOOL v. 2.00: a Matlab package for bifurcation analysis of delay differentialequations,”