Traffic flow with multiple quenched disorders

Traffic flow ^††NOTICE: this is the author’s version of a work that was accepted for publication in Physical Review E. Changes resulting from the publishing process, such as peer review, editing, corrections, structural formatting, and other quality control mechanisms may not be reflected in this document. Changes may have been made to this work since it was submitted for publication. A definitive version was subsequently published in Physical Review E DOI: 10.1103/PhysRevE.00.002100. is a classic problem in non-equilibrium statistical physics that has attracted the attention of physicists and engineers for decades. Various phenomena pertaining to traffic have been studied using a variety of models Chowdhury et al. (2000); Helbing (2001); Nagel et al. (2000). However, the phenomena that the system exhibits are so complex that there are many aspects of it which are not yet completely understood. One interesting aspect which drew the attention of physicists is the emergence of collective behavior in the traffic system in the presence of a quenched disorder in a physical parameter which led to power-laws and phase transitions.

Ben-Naim, Krapivsky and Redner (BKR) Ben-Naim et al. (1994) analytically derived the power-law behavior of various quantities like average speed, average platoon size (or length) of vehicles behind the slowest vehicle in the traffic stream when there is a quenched disorder in the free-flow speeds of the vehicles, which are the speeds of vehicles when unobstructed by other vehicles (also sometimes referred to as maximum speeds and desired speeds). They assumed that the vehicles (or the clusters of vehicles) move ballistically before interacting with their leader vehicle (cluster) and instantaneously merging with it. The assumed process is similar to the well-studied ballistic aggregation except that the momentum is not conserved in the present case. The system evolves into a stationary state in which there is a highly dense platoon of vehicles all moving at the same speed and a large gap ahead of the leader. BKR showed that power-laws for time evolution of the average platoon size behind the slowest vehicle and the average speed of the system have exponents which depend on the exponent of the free-flow speed distribution i.e., if $p(v)\sim v^{\mu}$ as $v\sim 0$,

where $N_{\mathrm{c}}$ is the cluster size and

Ben-Naim and Krapivsky went on further to develop a kinetic model for traffic flow Ben-Naim and Krapivsky (1998, 1999) which also includes vehicle overtaking.

Traffic flow has been extensively studied as a totally asymmetric simple exclusion process (TASEP) as originally introduced by Spitzer Spitzer (1970) and its variants. Krug and Ferrari Krug and Ferrari (1996); Krug (2000) studied the TASEP with quenched disorder in particle hopping rates and observed the formation of clusters with their size increasing with time. They conjectured that the power law for the cluster-size and speed should be similar to (Eq. (1) and Eq. (2)). Simulation studies by Bengrine et al. Bengrine et al. (1999) supported their conjecture. A further confirmation came from the studies of Ktitarev et al. Ktitarev et al. (1997) using a cellular automaton (CA) that is closely related to the TASEP with parallel updating of particle positions. Evans Evans (1996) independently studied the TASEP with particle-wise quenched disorder in the hopping rates and obtained exact solutions for the steady state of the system. Evans showed that the clustering of the vehicles is analogous to Bose-Einstein condensation. Krug and Ferrari Krug and Ferrari (1996) also studied the phase transition from the low density jammed phase, where the platoon formation happens, to the high density laminar phase with no such platoon formation. They observed that the type of phase transition depends upon the distribution of the quenched disorder. In both works it is assumed that the quenched disorder in the hopping rates has the form

The transition was observed to be of first order for $\mu>1$ and second order for $0<\mu\leq 1$. It was also observed that this phase transition need not match the phase transition from the free-flow phase to the congested phase.

Car-following models are another popular way of simulating the traffic flow. While CA models of traffic are inspired by that of a lattice gas, car-following models are continuum models in space and time. A variety of car following models exist in the literature Newell (2002); Helbing (2001); Treiber and Kesting (2013). These models, albeit time consuming to simulate, mimic real life traffic more closely than the TASEP based models and can have physically meaningful model parameters. This advantage made the car-following models highly popular in engineering applications and there are commercial software packages which can simulate various complex traffic flow conditions using the car-following models Fellendorf and Vortisch (2010); Behrisch et al. (2011).

There are limited studies on heterogeneity in car-following and there exist studies on the effect of slow vehicles Masukura et al. (2007, 2009); Li et al. (2016). However, to our best of knowledge there exist no studies on the effects of quenched disorders on probabilistic car-following. Further, the effect of having quenched disorders in multiple parameters is not studied even in the TASEP and CA models, and remains unknown. Viewing the quenched disorder as a heterogeneity in the vehicles, the actual traffic system can have heterogeneity in various parameters in addition to free-flow speed. For example, there can be heterogeneity in the minimum gap chosen by each driver, there can be heterogeneity in the reaction time of each driver and so on. The behavior of the traffic system would depend on the combined effect of various such parameters. Thus, in this work, we study the power-law behavior in various dynamical quantities and the phase transitions in the Newell’s car following model Newell (2002) by introducing quenched disorder into its parameters.

Newell’s car-following model, quenched disorders and the simulation details are briefly discussed in Sec. II. In Sec. III, the power-laws in various quantities observed in the simulations of the Newell’s model are discussed and analyzed. In Sec. IV, we study the stationary state gap distribution determined from the simulations and also develop a mean field theory to calculate it. In Sec. V, we study the phase transitions in the system using the simulations and analyze the results. The paper is concluded in Sec. VI.

Car-following models generally involve conjectured relationships between a vehicle’s speed (or equivalently their acceleration) and and the spatial gap between them and their leader and (in some models) the speeds of their leaders. Newell’s car-following model conjectures that the speed of vehicle $i$ at time $t$ is a function of the spatial gap from their leader $i-1$ at time $t-\Delta t_{\mathrm{adapt}}$, where $\Delta t_{\mathrm{adapt}}$ is a speed adaptation time (or a reaction time). The equation of motion is

where $s_{i}(t)\equiv x_{i-1}(t)-x_{i}(t)$ is the gap between vehicle $i$ and their leader $i-1$ at time $t$ and $v(\cdot)$ is a conjectured speed-spacing relation. Newell’s model basically assumes the following: (i) The path taken by the follower is a spatio-temporal translation of the path taken by the leader. (ii) The speed ($v$) of a vehicle depends on the gap ($s$) ahead of it as given below

We chose Newell’s car following model for our study as the model has only three parameters all of which are physically meaningful; the jam density ($\rho_{\mathrm{j}}$), the free-flow (maximum) speed of the vehicle ($v_{\mathrm{f}}$) and the backward wave speed ($w$). It is the simplest car-following model with interpretable parameters and it has been empirically validated by others Ahn et al. (2004); Chiabaut et al. (2009, 2010).

At a microscopic level, the density $\rho$ perceived by a driver may be interpreted as the reciprocal of the gap $s$ ahead of his/her vehicle. Thus, $S_{\mathrm{c}}=\rho_{\mathrm{c}}^{-1}$ is the critical gap,

is the critical density, and $S_{\mathrm{j}}=\rho_{\mathrm{j}}^{-1}$ is the reciprocal of jam density and is the minimum gap chosen by the driver. Eqs.(5,6) essentially mean that a triangular fundamental diagram is assumed for each driver. The equations allow for a simple interpretation that the vehicles do not interact as long as $s>S_{\mathrm{c}}$. The interactions are present only when $s\leq S_{\mathrm{c}}$ and are totally asymmetric i.e., the interaction effects only the follower.

The inputs to the simulation are (i) a track/circuit of chosen length, (ii) boundary conditions for the track, (iii) the initial spatial and speed distributions of the vehicles, and (iv) the appropriate choices of the driver-vehicle related parameters $\rho_{\mathrm{j}}$, $v_{\mathrm{f}}$, and $w$. If all the driver-vehicle units were identical, $\rho_{\mathrm{j}}$, $v_{\mathrm{f}}$ and $w$ would be the same for all the vehicles. We employ a forward Euler scheme to solve (4):

The time step length used to simulate the system, $\Delta t$, is chosen in such a way that the maximum displacement is less than the minimum of $S_{\mathrm{j}}$ so that no accidents happen. The position of each vehicle is simultaneously updated using the simple kinematic equations following the speed rule of Eq. (5). For the present work, we used periodic boundary conditions for the track. The initial spatial distribution of the vehicles is taken to be uniform and each vehicle is given its maximum speed so that the system is initialized in free-flow. Heterogeneity is introduced via a quenched disorder in each parameter of the model Jabari et al. (2014). Thus, each vehicle has a different set of parameters drawn from a chosen distribution. This means that each vehicle has a different fundamental diagram depending on its parameters. It has been shown in previous work by one of us Jabari et al. (2014, 2018); Zheng et al. (2018) that the quenched disorders in the parameters of the model can be accurately described by a distribution which is bounded from above and below. The traffic flow problem naturally appeals for a heterogeneity of that kind. For example, the free-flow (maximum) speed of vehicles is expected to vary between a minimum non-zero value and a maximum value. Similarly, the minimum gaps chosen by the drivers and the backward wave speed have strictly positive minimum values and finite upper bounds. In our simulations, we use the beta distribution for the parameters. It is a bounded domain distribution that takes as special cases other bounded domain distributions. It can be tuned to be symmetric, left, or right-skewed, i.e., it offers flexibility in simulation. While it is not the only distribution with these particular features, it is the most widely-known distribution with these features, and it is a conventional choice for traffic variables Haight (1963); Jabari et al. (2020). The general form of the (standard) beta distribution is

where $B(n,l)$ is the beta function with parameters $n$ and $l$, and $I_{\mathcal{X}}(x)$ is the indicator function, which returns 1 if $x\in\mathcal{X}$ and returns 0 otherwise. Let $A$ be any of the three parameters $v_{\mathrm{f}}$, $\rho_{\mathrm{j}}$ or $w$ with minimum and maximum values denoted by $A^{\min}$ and $A^{\max}$. The distribution of $A$, a generalized beta distribution is

In this work, we do not consider a quenched disorder in the speed adaptation time, $\Delta t_{\mathrm{adapt}}$. This has subtle impacts that we will consider in a future paper. In this paper, the delays are treated deterministically via the time step parameter $\Delta t$.

In this section, we determine from simulations the probability distribution of gap $p(s)$ for stationary states of the system in the regime where the BKR power-laws hold. We also develop a mean field theory to calculate $p(s)$.

To determine $p(s)$ from simulations, quenched disorder averaging is done over an ensemble of a few hundred copies of the stationary states of the system, depending on track length. A typical $p(s)$ we obtain is plotted in Fig. 7. The distribution $p(s)$ has two distinct components that require separate treatments, the first component corresponds to gaps among non-leading vehicles ($p_{\mathrm{p}}(s)$), i.e., within the platoon; the other corresponds to the gap ahead of the slowest vehicle, which is greater than the maximum gap within the platoon ($S_{\mathrm{p}}^{\max}$). The two components, $p_{\mathrm{p}}(s)$ and $p_{\mathrm{g}}(s)$ are not probabilities but are supported on disjoint sets. We may hence write

Consider a (following) vehicle ($f$) moving at a speed $v=v_{\mathrm{f}}>v_{\mathrm{f},l}$ approaching a platoon, where $v_{\mathrm{f},l}$ is the free-flow speed of the leader ($l$) of the platoon. The follower $f$ starts interacting with the last vehicle in the platoon in accordance with its gap choice characteristics and thus starts slowing down so as to equalize its speed to that of $l$. From Eq. (5), assuming $s>S_{\mathrm{j}}$, the gap from their leader upon joining the platoon (i.e., assuming the speed of the platoon leader) is

This gap is equal to a critical gap of the platoon leader $l$: $s=S_{\mathrm{c},l}$. Thus, each vehicle in the platoon would have a gap that is randomly chosen from the critical gap distribution of the leader of the platoon. We denote this gap by distribution $p_{\mathrm{c},l}(s)$. The above arguments are true even if the distributions of the quenched disorder are correlated. However, the present work assumes the beta distributions of the $w$ and the $S_{\mathrm{j}}$ are independent and hence $p_{\mathrm{c},l}(s)$ is just the product of the beta distributions of $w$ and $S_{\mathrm{j}}$ given a $v_{\mathrm{f},l}$.

We now derive the form of the stationary state gap distribution applying mean field assumptions. As the number of vehicles $N$ gets large, two things happen, all platoons coalesce into a single platoon and the speed of the leader tends to the minimum free-flow speed, i.e., $v_{\mathrm{f},l}\rightarrow v_{\mathrm{f}}^{\min}$ in the stationary limit and as $N$ gets large. Consequently, in this regime, the speed distribution tends to a point mass at $p(v)\rightarrow\delta(v-v_{\mathrm{f}}^{\min})$. Thus the distribution of gap is

We dropped the dependence on $v=v_{\mathrm{f}}^{\min}$ in the last equality to simplify notation and as this should be understood from context. For any vehicle selected at random from the $N$ vehicles in the system (i.e., chosen with probability $N^{-1}$), there are $N-1$ ways that a follower is chosen and one way that a leader is chosen. Let $p_{\mathrm{p}}(s_{1},\ldots,s_{N-1})$ denote the joint probability of the $N-1$ gaps in the platoon and $p_{l}(s)$ denote the probability density associated with the leader. Then

Note that the indices in the above equation do not represent the order of the vehicles; under stationary conditions, the gaps between the vehicles in a platoon are governed by their individual quench disorders and are independent of one another. The first term in the sum is $p_{\mathrm{p}}(s)$ and the second term is $p_{\mathrm{g}}(s)$. Since vehicles chose their gaps independently (under stationary conditions), the former simplifies to

Hence, the gaps within the platoon are dependent entirely on the leader’s critical gap choice. That is, $S_{\mathrm{p}}^{\max}=S_{\mathrm{c},l}^{\max}$ (in the limit) and from Eq. (6),

which has a bounded domain as $S_{\mathrm{j}}$ and $w$ both have bounded domains. We next analyze the second component

by analyzing the limiting behavior as $N$ and $L$ get large but maintain the relation $NL^{-1}=\overline{\rho}$. To do so, consider the distribution of

where $S_{\mathrm{p},i}$ is the $i$th gap within the platoon. For large $N$ (and $L$), this approaches a normally distributed random variable with mean $\langle S_{l}\rangle=L-(N-1)\langle S_{\mathrm{p}}\rangle$ and variance $(N-1)\sigma_{\mathrm{p}}^{2}\equiv(N-1)(\langle S_{\mathrm{p}}^{2}\rangle-\langle S_{\mathrm{p}}\rangle^{2})$. Hence, for large $N$

That the variance of $S_{l}$ is proportional to $N$ implies that the variability in the gap ahead of the platoon leader grows as the circuit gets longer and the number of vehicles is increased. However, the variability in this gap in proportion to $L$ shrinks with $N$. In other words, $S_{l}L^{-1}$ is normally distributed with mean $1-(N-1)L^{-1}\langle S_{\mathrm{p}}\rangle\approx 1-\overline{\rho}\langle S_{\mathrm{p}}\rangle$ and variance $(N-1)L^{-2}\sigma_{\mathrm{p}}^{2}=(N-1)N^{-2}\overline{\rho}\sigma_{\mathrm{p}}^{2}$.

It can be observed from the Fig. 7 that the stationary state distribution $p_{\mathrm{p}}(s)$ determined from our simulations is similar to $p_{\mathrm{c},l}(s)$ as determined above and that $p_{\mathrm{g}}(s)$ resembles a normal distribution (with small variance for large $N$). The above arguments breakdown when the speed is no more $v_{\mathrm{f}}^{\min}$.

Thus, the average queue size and the average relative speed variations show power-law behavior as derived by BKR, which is an artifact of the quenched disorder in $v_{\mathrm{f}}$. The effects of quenched disorders in $\rho_{\mathrm{j}}$ and $w$ are that $p_{\mathrm{p}}(s)$ and $S_{\mathrm{c}}$ are bounded. We checked validity of the results for different values of the beta distribution parameters $n$ and $l$ corresponding to PDF of $v_{\mathrm{f}}$.

As explained in the introduction, there are studies of phase transitions from the platoon formation phase to the laminar phase (where there is no platooning) in the TASEP and in the CA models. The TASEP simulations of Krug and Ferrari Krug and Ferrari (1996) and Bengrine et al. Bengrine et al. (1999) showed that the type of transition depends on the exponent of the distribution of hopping rates (Eq. (3)). They also observed that the transition point for the phase transition from the platoon formation phase to the laminar phase differs from the transition point for free-flow to congestion transition. Further, Krug and Ferrari showed that the gap distribution close to the transition is a power-law with the exponent depending on that of the hopping rate. Results of Ktitarev et al. Ktitarev et al. (1997) corroborate those observations even though the quenched disorder is in deceleration rates. In this section, we discuss the results of our simulations of the stationary states of the system and note the similarities and differences between our results and that of TASEP and CA models.

First we obtain the expressions for the average gap $\langle\bar{S}\rangle$ and the variance in the gap $\Delta_{S}^{2}$. We have that

as expected. Similarly,

and for large $N$, we have that

where we have used $\langle S_{\mathrm{p}}\rangle=\langle S_{\mathrm{c},l}\rangle$ in the limit based on the analysis above. From the last term, we see that the variance depends on the size of the system. This is a result of the formation of a single platoon of vehicles behind the slowest vehicle and the mere existence of a gap in front of the platoon leader, which becomes unbounded in size as $L$ gets large. As the density reaches the critical density of the platoon leader $\overline{\rho}\rightarrow S_{{\mathrm{c}},l}^{-1}$, the dependence on $L$ vanishes and $\Delta_{S}^{2}$ in this density regime is governed by intra-platoon gaps. Indeed, as $\overline{\rho}\rightarrow S_{\mathrm{c},l}^{-1}$, the gap in front of the platoon leader disappears and the platoon leader ceases to be the platoon leader.

Thus we take $\overline{\rho}_{\mathrm{c}}=\langle S_{\mathrm{c},l}\rangle^{-1}$ as the phase transition point and note that when $\overline{\rho}=\overline{\rho}_{\mathrm{c}}$ (i.e., at the phase transition), we have that

in the thermodynamic limit ($N\rightarrow\infty$, $L\rightarrow\infty$ and $NL^{-1}=\overline{\rho}_{\mathrm{c}}$). Krug and Ferrari characterized the phase transition to be of second order if the variance diverges as $\overline{\rho}_{\mathrm{c}}$ is approached. If the variance is finite, the transition would be of first order. In the present case, the variance $\Delta_{S}^{2}$ becomes finite as the phase transition is approached as can be seen from Eqs. (35) and  (36) implying that the transition is first order. As the above derivation doesn’t use any information regarding the quenched disorder in $v_{\mathrm{f}}$ except for the fact that there exists a $v^{\mathrm{min}}_{\mathrm{f}}$, the transition is expected to be of first order for any value of the parameters of the quenched disorder of $v_{\mathrm{f}}$.

Our simulations further confirm these results. The simulations are done for various $\overline{\rho}$ close to $\overline{\rho}_{\mathrm{c}}$. In order to characterize the transition, we determined the following quantities as a function of $\overline{\rho}$: the average maximum gap $\langle S^{\max}\rangle$ in the system, the probability of finding a gap greater than $S_{\mathrm{c},l}^{\max}$, i.e., $\mathbb{P}(S>S_{\mathrm{c},l}^{\max})$, $\Delta_{S}^{2}$, and the average flow $\langle J\rangle=\overline{\rho}\langle V\rangle$.

The left pane of Fig. 8 shows $\mathbb{P}(S>S_{\mathrm{c},l}^{\max})$ for various track lengths. The probability decreases sharply as the transition density $\overline{\rho}_{\mathrm{c}}$ is approached. The ensemble averaged maximum gap $\langle S^{\max}\rangle$ also decreases sharply as $\overline{\rho}_{\mathrm{c}}$ is approached. See the right pane of Fig. 8. The sharpness increases as the track length $L$ is increased indicating a discontinuity in the thermodynamic limit. In the inset of Fig. 8, the variances calculated from simulation and the theory (Eq. (35)) for a track of length 50 km are shown. An excellent match between them solidifies our theoretical findings. Thus, all the quantities we determined classify the transition as a first order transition. The flow-density diagram is plotted in Fig. 9.

The sharpness of the free-flow to congestion transition also increased with track length. The transition point in the thermodynamic limit is obtained by linearly extrapolating the free-flow curve and the congestion curves away from the transition point and by determining their intersection point. It can be seen that the transition point in the thermodynamic limit exactly matches with the predicted value of $\overline{\rho}_{\mathrm{c}}$. Thus, both the phase transitions happen at the same traffic density. Eq. (35) suggests that it would be first order for any values of parameters in quenched disorder for $v_{\mathrm{f}}^{\min}$. We confirmed it by performing simulations with other values of the parameters. See Fig. 10.

All the above results are in sharp contrast to those observed in the TASEP Krug and Ferrari (1996); Bengrine et al. (1999) and CA Ktitarev et al. (1997) simulations. Another interesting point to observe is that the triangular form of the flow-density diagram is preserved even after heterogeneity is included which says that the heterogeneous system may be mapped onto an effective homogeneous system with a suitable choice of parameters.

We have simulated the traffic system in its full generality within Newell’s model by incorporating quenched disorders in the free-flow speed, the jam density and the kinematic wave speed. The quenched disorders are modeled using beta distributions. We observed that system follows power-laws for average platoon size and the average speed as derived by Ben-Naim et al. We found that the gap distribution is a bounded domain distribution. To see if there is any equivalence between the sticky gas and the present system, we calculated the second moment of local traffic density, the density distribution and the spatial auto-correlation in local traffic densities for a particular set of parameters. We observed that the second moment of density and the density distribution do not follow a power-law indicating no equivalence between the present system and the sticky gas. As opposed to TASEP, there was no power-law behavior in the gap distribution at the critical density corresponding to a phase transition from the platoon formation phase to the laminar phase. Using mean field assumptions we derived the form of the gap distribution in the system in the platoon formation phase. It turns out that the gap distribution of vehicles in the platoon in the thermodynamic limit is simply the critical gap distribution of the leader. It is shown that the variance becomes finite as the phase transition point is approached and thus the transition is always of first order independent of the form of the quenched disorder, which is again in contrast to TASEP results. The transition density is shown to be nothing but the reciprocal of the average gap in the platoon. Simulation results support our theory.

Overall, the present results do not exactly conform with those of TASEP with quenched disorder. However, we expect that the present simulations closely mimic the real picture than the TASEP and related models. Thus the results of this work may be viewed as a benchmark for coarse grained models which use hydrodynamics or kinetic theory and incorporate heterogeneity.