Macroscopic Manifestations of Traffic Waves in Microscopic Models

Microscopic car-following traffic models are systems of ordinary differential equations (ODEs), where the dynamics of each vehicle are described by a dynamic equation of motion. To derive the trajectories of individual vehicles, the ODEs could prescribe the vehicle velocity (first-order models), or the velocity and the acceleration (second-order models). Second-order car-following models are of the form

$$\dot{v}(t)=f(s(t),v(t),\Delta v(t))\;.$$

(1)

Here $s$ is the gap (measured in metres) to the vehicle ahead (the “lead vehicle”), $v$ is the vehicle’s velocity (measured in $\text{m}/\text{s}$), and $\Delta v$ is the relative velocity between the lead vehicle and the vehicle itself, defined for vehicle $i$ as $\Delta v_{i}=v_{i-1}(t)-v_{i}(t)$, where $v_{i}$ is its own speed and $v_{i-1}$ the lead vehicle’s speed. Many second order models have been proposed; below we discuss two important examples of them.

Macroscopic traffic models come in the form of partial differential equations (PDEs) that describe aggregate quantities such as traffic density or flow rate, instead of resolving the individual vehicles. A fundamental example is the Lighthill-Whitham-Richards (LWR) model. It is a first-order hyperbolic conservation law, describing the density of vehicles on a road with no entries or exits. It was first proposed in [16, 24] as an equation describing the relation between density $\rho$ and flow rate $q$ as

$$\frac{\partial{\rho}}{\partial t}+\frac{\partial{q}}{\partial x}=0\;.$$

(9)

Equation (9) is often referred to as the continuity equation, and will be used in later sections to denote a FD relationship $q=Q(\rho)$, via a bulk velocity $u$, where $q=\rho u$. The scalar conservation law (9) satisfies a maximum principle, hence it is devoid of dynamic instabilities and cannot be used as a model for phantom traffic jams. It is important to note that other macroscopic models exist that allow the study of phantom traffic jams like the inhomogeneous Aw-Rascle-Zhang (ARZ) model [2, 35, 9], which provides a second evolution equation for the bulk velocity field and can be seen as a special case of generic macroscopic second-order models [13, 6]. Note that second-order macroscopic models are not considered here as this work only studies the manifestations of microscopic models, but not their macroscopic limits.

In this work, the setups used for the traffic simulations can be categorised into three main types: ring roads, infinite roads, and bottlenecks. For all setups, we consider vehicles on single-lane roads with positions $x_{i}(t)$ and speeds $v_{i}(t)$, where vehicle $i$ follows vehicle $(i-1)$ and the equations of motion of all vehicles follow (1). In a ring road setup, $N$ vehicles circulate in a loop, where the $1^{\text{st}}$ vehicle follows the $N^{\text{th}}$. Even though a ring road is not a realistic representation of traffic on highways, this setup creates a controlled environment to study congestion, traffic waves, and car-following behaviour, as for example employed in real-world experiments [28, 29, 27]. The infinite road setup simulates traffic on an unbounded road providing a realistic representation of traffic on real-world highways when a lead vehicle with a large gap ahead is present. Whereas, the bottleneck road setup simulates traffic on a bounded road segment with prescribed inflow and outflow conditions, where the outflow condition often causes localised disruption in traffic and growing traffic jams. In real life, a bottleneck can be caused by the geometry of the road or reductions in speed limits.

In our simulations, the infinite road setup consists of a platoon of vehicles. The lead vehicle in that platoon is moving always at a constant speed and is not affected by any added noise or perturbations. The road segment with a bottleneck is implemented via (i) an inflow layer where at each time step vehicles are seeded onto the road segment at a fixed inflow speed and spacing (satisfying an equilibrium relationship) whenever the vehicles already on the road have moved to create enough space, and (ii) an outflow segment where the speeds of every vehicle that has left the domain is set to a fixed value (lower than the inflow speed). Those vehicles remain being tracked until they can be removed because they have ceased to affect any other vehicles. Note that, because below we employ the bottleneck scenario for congested flow, traffic flow theory dictates that it is the outflow conditions that will transport information into the domain. We use those different setups in the scenarios of §4 as each offers unique insights into traffic behaviour.

The stability of a given equilibrium state can be studied through a linear stability analysis [33]. We linearise the equations of motion around an equilibrium state by choosing the position of vehicle $i$ to be $x_{i}(t)=(s_{\text{eq}}+\ell)i+v_{\text{eq}}t+y_{i}$, where $y_{i}$ is an infinitesimal perturbation. Then, we substitute these vehicle trajectories into (1), Taylor-expand around the equilibrium state, and ignore all higher order terms. The perturbation equation obtained is

	$\displaystyle\ddot{y}_{i}(t)$	$\displaystyle=\alpha_{1}(y_{i-1}-y_{i})-\alpha_{2}\dot{y}_{i}+\alpha_{3}\dot{y}_{i-1}\;,$		(10)
	$\displaystyle\text{where}\qquad\alpha_{1}=\frac{\partial f}{\partial s}\;,$	$\displaystyle\qquad\alpha_{2}=\frac{\partial f}{\partial(\Delta v)}-\frac{\partial f}{\partial v}\;,\qquad\alpha_{3}=\frac{\partial f}{\partial(\Delta v)}\;,\qquad$		(11)

and all the partial derivatives are evaluated at the chosen equilibrium state. Then, the growth or decay of solutions to (10) is characterised by performing a Laplace transform ansatz $y_{i}(t)=c_{i}e^{\omega t}$, where $c_{i},\omega\in\mathbb{C}$. This yields the interpretation of (10) as an input/output (I/O) system, $c_{i}=F(\omega)c_{i-1}$, where the transfer function is

$$F(\omega)=\frac{\alpha_{1}+\alpha_{3}\omega}{\alpha_{1}+\alpha_{2}\omega+\omega^{2}}\;.$$

(12)

The temporal growth or decay of the lead vehicle’s velocity profile is captured by $\mathrm{Re}(\omega)$ in equation (12). The frequency of oscillation of the lead vehicle’s velocity profile is represented by $\mathrm{Im}(\omega)$ in equation (12). The transfer function’s modulus $|F|$ is the growth or decay of the perturbation amplitude from one vehicle to the next. An equilibrium state is defined to be string unstable if speed and spacing fluctuations of each vehicle are greater in amplitude than those preceding it [33]. String instability is a property that can eventually (when growing out of the linear regime) give rise to stop-and-go wave structures in car-following model simulations. With the setup in (12), stability means that $|F(\omega)|\leq 1\;\forall\omega\in i\mathbb{R}$ (considering $\mathrm{Re}(\omega)=0$ is the only way to ensure that $e^{\omega t}$ remains infinitesimally small for all $t\in\mathbb{R}$). Evaluating $|F(\omega)|^{2}$ for $\omega\in i\mathbb{R}$ yields that this stability criterion can be written as a condition on the partial derivatives of $f$, as follows:

$$\alpha_{2}^{2}-\alpha_{3}^{2}-2\alpha_{1}\geq 0\;.$$

(13)

Hence, an equilibrium state’s stability is determined by the partial derivatives of $f(s,v,\Delta v)$ with respect to the state variables at that equilibrium state and (13).

For example, some straightforward calculations yield that for the IDM, stability vs. instability of an equilibrium background density depends only on the choices of the IDM parameters $a$ and $b$. The other parameters determine the shape of the fundamental diagram and are chosen accordingly, where $s_{0}$ affects the jam density or the maximum traffic road density (reasonable values around $2\text{m}$), $v_{0}$ affects the maximum flow rate (reasonable values between $30\text{m}/\text{s}$ and $45\text{m}/\text{s}$, but is often set to the maximum speed limit value), $T$ affects the maximum flow rate (reasonable values around $1\text{s}^{-1}$), and $\delta$ affects the curvature of the fundamental diagram (usually set to $4$). Figure 1 shows for different choices of $a$ and $b$ the density value at which the IDM pivots from stable to unstable. It also shows an example of the stability vs. instability regions of the IDM on the level of the fundamental diagram for specific parameter values $a=1.3\text{m}/\text{s}^{2}$ and $b=2\text{m}/\text{s}^{2}$. It should be noted that when analyzing real-world traffic data the maximum density of vehicles that a road can accommodate while maintaining a relatively smooth and free-flowing traffic condition is often referred to as the “critical density” beyond which traffic flow begins to degrade significantly. However, the present analysis of traffic models highlights that the maximum-capacity density and the onset-of-instability density are actually different quantities and may have quite different values. Here, the term “critical density” refers to the former.

Simulating the ODE-based microscopic models is carried out using the explicit time-stepping numerical scheme Forward Euler (FE) method. In the absence of noise (see below), one could also use higher order Runge-Kutta methods (e.g., RK4) to gain higher accuracy. In all simulations in this work we choose a time step of size $\Delta t=0.1\text{s}$. This choice ensures that approximation errors induced by the numerical schemes will not create any undesirable effects like vehicle crossings, non-convergence in velocity/position trajectories, or vehicles travelling too far in one time step.

To obtain a simulation of a model that properly captures unstable behaviour, we augment the velocity solution of the model in (1) by a scaled white noise term (chosen to be a Gaussian with a mean of zero). This means that vehicle trajectories are the solution to a stochastic differential equation arising by augmenting (1) with an additive Brownian motion. The solution to such stochastic differential equations is a velocity with noise scaled so that the impact of the noise over a time interval is a normally distributed random perturbation with mean zero and variance equal to $\sigma^{2}$ times the interval’s length, where $\sigma$ denotes the magnitude of the noise. An Euler-Maruyama discretisation of the stochastic differential equation for a discretisation time step $\Delta t$ for the $i^{\text{th}}$ vehicle is given as

	$\displaystyle x_{i}(t+\Delta t)$	$\displaystyle=x_{i}(t)+\Delta tv_{i}(t)\;,$		(14)
	$\displaystyle v_{i}(t+\Delta t)$	$\displaystyle=v_{i}(t)+f(s_{i}(t),v_{i}(t),\Delta v_{i}(t))+\sqrt{\Delta t}~{}\mu(0,\sigma)\;.$		(15)

The noise in the simulations serves as a simplified model of the inaccuracies that human drivers introduce while trying to maintain a uniform flow driving on a road, or other fast perturbations like gusts of wind. Figure 2 shows the importance of including perturbations/noise in simulations for wave development. Two simulations of $60$ vehicles on a $1500\text{m}$ ring-road are initialised at equilibrium, following the IDM with $a=1.3\text{m}/\text{s}^{2}$ and $b=2\text{m}/\text{s}^{2}$. The background density of this simulation is $40\text{~{}veh}/\text{km}$, so it lies in the unstable region for the given $a$ and $b$ values (see fundamental diagram in Figure 1), and with enough perturbations, the model dynamics develop waves. The first simulation is carried without noise, and we can see that the speed profile of a vehicle in this simulation remains constant. On the other hand, the second simulation is carried out with the instability triggered by adding noise during the first $500\text{s}$ only, and we see the subsequent development and growth of waves from the variations in the speed profile.

To establish a connection between the micro and macro scales, we derive macroscopic quantities from microscopic models. There are several ways established to do so (like in [14]); one is to use Gaussian kernels of the form

$$G(x)=Z^{-1}e^{-(\frac{x}{h})^{2}}\;,$$

(16)

where $h$ is the width of the kernel, chosen to be wider than the minimum spacing between vehicles but small enough to capture localised information about the density and flow rate. For $N$ vehicles on the road, we reconstruct via standard techniques the density as the superposition of these Gaussian kernels

$$\rho(x,t)=\sum_{j=1}^{N}G(x-x_{j}(t))\;,$$

(17)

and the flow rate as also the superposition of kernels weighted by velocities

$$q(x,t)=\sum_{j=1}^{N}\dot{x}_{j}(t)G(x-x_{j}(t))\;.$$

(18)

Consequently we can define the reconstruction of bulk traffic velocity as

$$u(x,t)=\frac{q(x,t)}{\rho(x,t)}\;.$$

(19)

The task of choosing the suitable kernel width $h$ is an art by itself, but one should keep in mind that very large $h$ values would smear out features that might be interesting, and very small choices of $h$ would introduce unwanted oscillations.

A crucial advantage of using the above technique to reconstruct macroscopic densities and flow rates is that the reconstructed $(\rho,q)$-pairs exactly satisfy the continuity equation (9), and such solutions form lines on the fundamental diagram [12]. To prove that latter fact, we consider a solution of the continuity equation (9) that is a single profile moving with speed $s$ along the road. Then we can write the quantities $\rho$ and $q$ as functions of a single variable $\eta=x-st$. This implies that $\rho_{t}=-s\frac{d\rho}{d\eta}$ and $q_{x}=\frac{dq}{d\eta}$, and thus (9) becomes $\frac{d}{d\eta}(-s\rho+q)=0$. Integration yields $-s\rho+q=m$, where the integration constant $m$ is the mass flux of vehicles relative to the wave. Thus one obtains the relationship

$$q=m+s\rho\;,$$

meaning that pairs of density and flow rate, $(\rho,q)$, that satisfy (9) and are travelling waves, form straight lines in the fundamental diagram plane, where the line’s slope $s$ equals the speed of the travelling wave.

Therefore, if we plot point-wise all the reconstructed pairs $(\rho,q)$ at a specific instance in time when waves are fully developed they will fall on a line of a slope equal to the speed of the travelling wave (see §3.5) and the average of those point-wise (in both time and space) reconstructed $(\rho,q)$-pairs we call the effective or average state (see Figure 3 as an example with the blue dot representing the average state). For a given simulation setup, the curve formed by all the effective states (emerging from many different initial equilibrium states) in the unstable regime is denoted as the reduced fundamental diagram curve. Given an initial equilibrium state in the unstable regime, we show interest in determining how the corresponding average state moves in the $(\rho,q)$-space as waves develop in different setups where certain quantities like density, flow rate, and bulk velocity are held fixed (see §4).

Waves are travelling perturbations in the distribution of vehicles and travelling backwards with respect to vehicles on the road. The state of traffic flow can be determined by the variation of the average vehicle speeds on the road. If this average is high and close to the road’s speed limit, then we are in a free flow state. Otherwise, we are in a traffic wave state ranging in intensity from mild (slower vehicles in congestion but with no stoppage) to stop-and-go traffic (experiencing complete stoppage). Quantifying the intensity of traffic waves is affected by the margin of variation of the speeds of vehicles on the road. In this work, we show specific interest in traffic waves and how to interpret them as structures in the fundamental diagram.

Macroscopically, we consider two points in the ($\rho$,$q$)-plane representing the different traffic states: the right state $(\rho_{R},q_{R})$ and left state $(\rho_{L},q_{L})$. We define the speed of the shock based on the Rankine–Hugoniot conditions [15] as:

$$v_{s}=\frac{q_{R}-q_{L}}{\rho_{R}-\rho_{L}}\;.$$

(20)

For a macroscopic interpretation of microscopic waves, we define a ‘jamiton line’ as the nonlinear travelling wave solution on the fundamental diagram [25]. The word jamiton was introduced by the authors of [7] in analogy to the nonlinear travelling waves in physics called solitons. In the ($\rho$,$q$)-plane, the jamiton line is determined as the best fit line of macroscopic density and flow rate data extracted using Gaussian kernels at a specific instance in time from microscopic simulations. It should be stressed that since we are postprocessing numerical results, the fit is needed because (a) a travelling wave solution may not be fully established yet; and due to (b) numerical approximation errors; and (c) noise effects. The fit will be a straight line segment with slope $v_{s}$ and two end points $(\rho_{R},q_{R})$ and $(\rho_{L},q_{L})$ representing the two traffic states across the wave (see Figure 3 where the endpoints of the line across the fundamental diagram represents the left and right states).

With a careful choice of parameters in both the IDM and OVM, the two models can be reasonably compared. For all the comparisons conducted below we use the following parameters: in (8) we set $\alpha=1.085\frac{1}{\text{s}}$, $\beta=22.0779\frac{\text{m}^{2}}{\text{s}}$, $\nu=2$, $a_{\text{m}}=1.3\frac{\text{m}}{\text{s}^{2}}$, $b_{\text{m}}=5\frac{\text{m}}{\text{s}^{2}}$, $s_{0}=2\text{m}$, $v_{0}=30\frac{\text{m}}{\text{s}}$, whereas in (2) we set $a=1.3\frac{\text{m}}{\text{s}^{2}}$, $b=2\frac{\text{m}}{\text{s}^{2}}$, $v_{0}=30\frac{\text{m}}{\text{s}}$, $s_{0}=2\text{m}$, $T=\frac{1}{\text{s}}$, and $\delta=2$. Those parameters are chosen based on the principled methododology explained below.

To be able to compare jamiton lines of two different models, the two models should first have the same fundamental diagrams, i.e., they should exhibit the same equilibrium $(s,v)$-relationships from $f(s,v,0)_{\text{IDM}}=f(s,v,0)_{\text{OVM-Modified}}=0$.

To match the fundamental diagram of the IDM and OVM we start with an IDM parameter $\delta=2$ (rather than the popular choice $\delta=4$, because for $\delta=2$ the resulting mathematical expressions become simpler) so the IDM equilibrium spacing function can be easily inverted and is given by

$$S(v)=\frac{s_{0}+Tv}{\sqrt{1-(\frac{v}{v_{0}})^{2}}}\;.$$

(21)

From (21) we derive a formula for the OVM optimal velocity function that reads as

$$V(s)=\frac{-s_{0}+\sqrt{s_{0}^{2}-(s_{0}^{2}-s^{2})(\frac{s^{2}}{T^{2}v_{0}^{2}}+1)}}{T(\frac{s^{2}}{T^{2}v_{0}^{2}}+1)}\;.$$

(22)

Note that $V(s)$ is zero for the minimum gap $s_{0}$ and asymptotes at a maximum velocity $v_{0}$ for $s\to\infty$.

After achieving the same fundamental diagram equilibrium curves for both models, we now additionally ensure that the stability vs. instability regions on those curves are also comparable for both models. For that, we start with IDM parameter choices $a$ and $b$ that set the stability regions. We then ensure that the OVM model parameters, namely $\alpha$ and $\beta$, are uniquely determined such that the maximum growth rate of waves and the range of unstable background densities of interest are comparable to those of the IDM. To establish that, we enforce two constraints:

1. (1)

   The background density at which the corresponding equilibrium solution for both models pivot from stable to unstable is the same. This imposes a linear relation between $\alpha$ and $\beta$ that reads as $\beta=(V^{\prime}(s_{s})-\frac{\alpha}{2})s_{s}^{2}$ , where $s_{s}$ is the equilibrium spacing corresponding to the background density at which the IDM model we started with pivots from stable to unstable.

2. (2)

   The OVM parameter $\alpha$ (and consequently $\beta$) is determined by solving the minimisation problem

   $$\min_{\alpha}\int_{\Omega_{\rho}}\lvert\max_{\omega}F_{IDM}(\rho,\omega)-\max_{\omega}F_{OVM}(\rho,\omega)\rvert^{2}q(\rho)d\rho\;,$$

   which minimises the difference of the maximum growth rate of waves $F(w)$ defined in (12) of both models over all possible densities weighted by the flow rate.

After these cross-calibrations of the these two important microscopic models, OVM vs. IDM, have been conducted, we can now compare them in terms of their jamiton lines generated for different background densities.

Applying the methodologies described above to compare the two models, we find a striking fundamental difference in the waves properties.

Considering the fundamental diagram of a second order microscopic car-following model, and starting with a specific equilibrium state on the fundamental diagram where vehicles are equispaced and moving at the same speed, how does the effective (or average) density and flow rate state evolve in the ($\rho$,$q$)-space as we perturb the system?
To answer that we consider three different scenarios where we fix certain quantities (density, velocity, flow rate) and study the effect of fixing such quantities on the evolution of the effective state.

Note that the presented scenarios are crucial situations that apply in real-world experiments and/or configurations, thus they are also important building blocks for simulations. As highlighted earlier, to isolate the car-following behaviour in this work we only consider a simple case where heterogeneity in driving behaviour and lane switching are not modelled. One can draw an analogy of car-following behaviour in the presence of waves to gas dynamics in higher temperature where more local oscillations rise, and as a consequence the same amount of gas/vehicles “wiggle” more, thus occupying more space/volume. In all three scenarios, we start with equi-spaced vehicles initialised at equilibrium speed and we add perturbations to the system at all times to check how the effective state evolves. Notice that for all three scenarios corresponding to Figures 8, 9, and 10 (top), as time evolves, the velocities of all vehicles on the road pass through a transient phase before being suitably close to the travelling wave state limit $t\to\infty$ when traffic waves are well developed. Therefore, in all scenarios it is important to consider the the regions in space-time where/when the waves are fully developed, this includes (a) for all scenarios since we start with equi-spaced configuration, we need to remove the initial and transient layers in time, in addition (b) for scenarios 2 and 3 we also need to remove spacial boundary layers where waves are not fully developed.

In this subsection, we (i) ensure the robustness of the process of finding the effective/average state and its evolution in the density-flow rate plane for each scenario, (ii) determine the curves that the various average states form, and (iii) compare the effective states curves resulting in the reduced fundamental diagram for the IDM vs. OVM in all three scenarios. We run an ensemble of simulations at each background density of interest with the initial equilibrium states chosen to be in the unstable regime and the density is not extremely high (i.e., not close to the maximum density). This ensures that waves will develop, and eliminates non-interesting cases of fully occupied roads. For each background density, the mean of the calculated average states from an ensemble of 25 simulations is plotted for the IDM and OVM in Figures 11 and 12 respectively, with error bars in both the density and flow rate direction representing the standard deviations over the respective ensembles. The same methodology is carried for all three scenarios at each initial equilibrium state.

The results verify that for both the IDM and OVM, the effective states start at the initial equilibrium state and move away from the fundamental diagram: (i) downwards on a vertical line in the ring road scenario, (ii) towards the origin on a line connecting the initial equilibrium state with the origin in the infinite road scenario, and (iii) to the left on a horizontal line in the bottleneck scenario, for all unstable background densities. This confirms that although each simulation has a random component, in all three scenarios we have the same fundamental, systematic principle: waves correspond to the same kind of car-following behaviour that the vehicles exhibit on the road, independent of the scenario set up; it is only the manifestations where the corresponding effective states end up in the ($\rho$,$q$)-space that is different depending on the setup.

A further important observation can be made about the shape of the reduced fundamental diagram for the IDM vs. OVM. Note that for the IDM the averaged effective states for all unstable background densities and from all three different scenarios are reasonably close to one single line (a minor exception is for lower background densities that are close to the boundary of instability), verifying that the shape of the reduced fundamental diagram is scenario-independent as all effective states for all scenarios in the unstable regime lie on the same line (see Figure 11). For the OVM, the averaged effective states for all unstable background densities and from all three different scenarios are reasonably close to a fitted jamiton curve rather than a line (see Figure 12). This also verifies that the shape of the reduced OVM fundamental diagram is scenario-independent as all effective states for all scenarios in the unstable regime lie on the same curve. The results in this study show that the presence of traffic waves always yields the same effective velocity-spacing relationship, thus microscopic models in the unstable flow regime exhibit many of the fundamental structural properties that have been proven for second-order macroscopic traffic models [25, 23].