The BGK approximation of kinetic models for traffic

TL;DR

This paper introduces a modified BGK kinetic model for traffic flow that achieves conditional stability and captures stop-and-go waves, bridging microscopic and macroscopic traffic models.

Contribution

A new BGK-type kinetic model is derived that stabilizes solutions in congested traffic and reproduces realistic wave phenomena, enhancing traffic flow modeling.

Findings

The modified model reproduces stop-and-go waves as backward propagating signals.

It provides a mesoscopic description linking microscopic and macroscopic models.

The model achieves conditional stability in congested traffic regimes.

Abstract

We study spatially non-homogeneous kinetic models for vehicular traffic flow. Classical formulations, as for instance the BGK equation, lead to unconditionally unstable solutions in the congested regime of traffic. We address this issue by deriving a modified formulation of the BGK-type equation. The new kinetic model allows to reproduce conditionally stable non-equilibrium phenomena in traffic flow. In particular, stop and go waves appear as bounded backward propagating signals occurring in bounded regimes of the density where the model is unstable. The BGK-type model introduced here also offers the mesoscopic description between the microscopic follow-the-leader model and the macroscopic Aw-Rascle and Zhang model.

Figures

Tables

Table 1. Table 3.1 : Classification of flow regimes for the Boltzmann-type kinetic equation ( 2.1 ).

$ϵ$	Regime	Kinetic model	Continuum flow model
$0$	Equilibrium flow	Boltzmann	Mass conservation law (LWR model)
$\to 0$	Viscous flow		Diffusion equation
$\asymp 1$	Transitional		Extended hydrodynamic equations
$>1$	Rarefied		-
$\to \infty$	Free molecular flow	Collision-less Boltzmann	-

Equations

$${\partial }_{t}f(t,x,v)+v{\partial }_{x}f(t,x,v)=\frac{1}{ϵ}Q[f,f](t,x,v),$$

$$f(t,x,v):{\mathbb{R}}^{+}\times \mathbb{R}\times V\to {\mathbb{R}}^{+}$$

$$\rho (t,x)={\int }_{V}f(t,x,v)\mathrm{d}v,(\rho u)(t,x)={\int }_{V}vf(t,x,v)\mathrm{d}v,u(t,x)=\frac{1}{\rho }{\int }_{V}vf(t,x,v)\mathrm{d}v$$

$${\partial }_{t}f(t,v)=\frac{1}{ϵ}Q[f,f](t,v).$$

$$Q[f,f](t,x,v)={\int }_V\mathrm{d}{v}_{\ast }{\int }_V\mathrm{d}{v}^{\ast }f(t,v_{\ast },x)\mathcal{P}(v_{\ast }\to v|v^{\ast };\rho )f(t,v^{\ast },x)-f(t,v,x){\int }_V\mathrm{d}{v}^{\ast }f(t,v^{\ast },x),$$

$$\mathcal{P}(v_{\ast }\to v|v^{\ast };\rho )=\{\begin{array}{ll}P\delta (v-\min \{v_{\ast }+{\Delta }_a,V_M\})+(1-P)\delta (v-\max \{v_{\ast }-{\Delta }_b,0\})& v_{\ast }\le v^{\ast }\\ P\delta (v-\min \{v_{\ast }+{\Delta }_a,V_M\})+(1-P)\delta (v-\max \{v^{\ast }-{\Delta }_b,0\})& v_{\ast }>v^{\ast }\end{array}$$

$$\mathcal{P}(v_{\ast }\to v|v^{\ast };\rho )\ge 0,{\int }_V\mathcal{P}(v_{\ast }\to v|v^{\ast };\rho )\mathrm{d}v=1.$$

$$M_f(v;\rho )=\sum _{j=1}^N{f}_j^{\infty }(\rho ){\delta }_{v_j}(v),f_j^{\infty }>0\forall j=1,\dots ,N,$$

$$f_j^{\infty }=\{\begin{array}{ll}0& P\ge \frac{1}{2}\\ \frac{-2(1-P){\sum }_{k=1}^{j-1}{f}_k^{\infty }+(1-2P)\rho +\sqrt{{\left[(1-2P)\rho -2(1-P){\sum }_{k=1}^{j-1}{f}_k^{\infty }\right]}^2+4P(1-P)\rho f_{j-1}^{\infty }}}{2(1-P)}& \text{mbox}otherwise\end{array}$$

$$Q_{\text{eq}}(\rho )=\left(\rho U_{\text{eq}}(\rho )\right)={\int }_{V}v{M}_f(v;\rho )\mathrm{d}v,U_{\text{eq}}(\rho )=\frac{1}{\rho }{\int }_{V}v{M}_f(v;\rho )\mathrm{d}v.$$

$${\partial }_t\rho +{\partial }_x{Q}_{\text{eq}}(\rho )=0$$

$${\partial }_{t}f(t,x,v)+v{\partial }_{x}f(t,x,v)=\frac{1}{ϵ}\left(M_f(v;\rho )-f(t,x,v)\right),$$

$$f_{ij}^{(1)}$$

$$f_{ij}^{n+1}$$

$$f_{ij}^{n+1}=M_f(v_j;{\rho }_i^n)-\frac{\delta t}{\delta x}\left(\mathcal{F}\left(M_f(v_j;{\rho }_{i+1}^n),M_f(v_j;{\rho }_i^n)\right)-\mathcal{F}\left(M_f(v_j;{\rho }_i^n),M_f(v_j;{\rho }_{i-1}^n)\right)\right),$$

$${\rho }_i^{n+1}={\rho }_i^n-\frac{\delta t}{\delta x}\left(\mathcal{F}\left({\rho }_{i+1}^n,{\rho }_i^n\right)-\mathcal{F}\left({\rho }_i^n,{\rho }_{i-1}^n\right)\right).$$

$$\mathcal{F}\left(f_{i+1,j}^{(1)},f_{ij}^{(1)}\right)=\frac{1}{2}\left(v_j{f}_{i+1,j}^{(1)}+v_j{f}_{ij}^{(1)}-\alpha \left(f_{i+1,j}^{(1)}-f_{i+1,j}^{(1)}\right)\right).$$

$${\rho }_0(x)={\rho }_{\min }+\left({\rho }_{\max }-{\rho }_{\min }\right)\exp \left(-50(x-0.5{)}^2\right).$$

$$\underset{P\to {\frac{1}{2}}^{+}}{\lim }{\int }_{\mathcal{V}}v{M}_f(v;\rho )\mathrm{d}v=\underset{P\to {\frac{1}{2}}^{-}}{\lim }{\int }_{\mathcal{V}}v{M}_f(v;\rho )\mathrm{d}v={\rho }_{\text{crit}}{V}_M.$$

$${\partial }_t\rho +{\partial }_x{Q}_{\text{eq}}(\rho )=ϵ{\partial }_x(\mu (\rho ){\rho }_x),$$

$$\mu (\rho )={\int }_V{v}^2{\partial }_{\rho }{M}_f(v;\rho )\mathrm{d}v-{\left({\int }_{V}v{\partial }_{\rho }{M}_f(v;\rho )\mathrm{d}v\right)}^2={\int }_V{v}^2{\partial }_{\rho }{M}_f(v;\rho )\mathrm{d}v-Q_{\text{eq}}^{\prime }(\rho {)}^2.$$

$$f(t,x,v)=M_f(v;\rho )+ϵf_1(t,x,v),\text{with }{\int }_V{f}_1(t,x,v)\mathrm{d}v=0,$$

$$f_1(t,x,v)=-{\partial }_t{M}_f(v;\rho )-v{\partial }_x{M}_f(v;\rho )+O(ϵ),$$

$${\partial }_t\rho +{\partial }_x{Q}_{\text{eq}}(\rho )$$

$$=ϵ{\partial }_x{\int }_{V}v({\partial }_t{M}_f(v;\rho )+v{\partial }_x{M}_f(v;\rho ))\mathrm{d}v+O({ϵ}^2)$$

$$=ϵ{\partial }_x\left[{\rho }_t{\int }_{V}v{\partial }_{\rho }{M}_f(v;\rho )\mathrm{d}v+{\rho }_x{\int }_V{v}^2{\partial }_{\rho }{M}_f(v;\rho )\right]+O({ϵ}^2)$$

$$=ϵ{\partial }_x\left[-{\rho }_x{Q}_{\text{eq}}^{\prime }(\rho ){\int }_{V}v{\partial }_{\rho }{M}_f(v;\rho )\mathrm{d}v+{\rho }_x{\int }_V{v}^2{\partial }_{\rho }{M}_f(v;\rho )\right]+O({ϵ}^2).$$

$${\partial }_{\rho }\left({\int }_V{v}^2{M}_f(v;\rho )\mathrm{d}v\right)>Q_{\text{eq}}^{\prime }(\rho {)}^2$$

$$Q_{\text{eq}}^{\prime }(\rho )={\int }_{v}v{\partial }_{\rho }{M}_f(v;\rho )\mathrm{d}v<0,{\partial }_{\rho }\mathrm{Var}(v)={\partial }_{\rho }{\int }_V(v-U_{\text{eq}}(\rho ){)}^2{M}_f(v;\rho )\mathrm{d}v<0$$

$${\partial }_{\rho }{\int }_V{v}^2{M}_f(v;\rho )\mathrm{d}v$$

$${\partial }_t\rho +{\partial }_x(\rho u)$$

$${\partial }_t(\rho u)+{\partial }_x{\int }_V{v}^{2}f(t,x,v)\mathrm{d}v$$

$${\int }_V{v}^{2}f(t,x,v)\mathrm{d}v\approx \rho u^2+h(\rho )$$

Full text

MSC

90B20, 35Q20, 35Q70

Keywords

Vehicular traffic, BGK models, Chapman-Enskog expansion, multi-scale models, kinetic equations

1 Introduction

There are mainly three modeling scales in the mathematical description of vehicular traffic flow. The microscopic one is based on the prediction of the trajectory of each single vehicle via systems of ordinary differential equations. The macroscopic one is based on the assumption that traffic flow behaves like a fluid where individual vehicles cannot be identified. Therefore, the flow is represented by a density function and evolved in space and time via partial differential equations. The last one is the mesoscopic scale which offers an intermediate description of traffic between the microscopic and the macroscopic scale. In fact, mesoscopic or kinetic models are characterized by a statistical description of the microscopic states of vehicles but, at the same time, still provide the macroscopic aggregate representation of traffic flow e.g. linking collective dynamics to pairwise interactions among vehicles at a smaller microscopic scale.

Prototype examples of microscopic models are given by the optimal velocity model by Bando [3] and the follow-the-leader model [15]. The first and classical example of a macroscopic model for traffic was introduced by Lighthill and Whitham [30] and independently by Richards [44]. This model is based on the single conservation law for the density assuming that the dynamics is governed by a velocity function which is a function of the density only. An extension of this is provided by the second-order macroscopic model introduced by Aw and Rascle [2] and by Zhang [48]. It is based on a system of two partial differential equations describing also the evolution of the speed of the flow. There is a proved link between common microscopic and macroscopic models. More precisely, the follow-the-leader model approaches to the Lighthill-Wihtham-Richards model and to the Aw-Rascle-Zhang model in the limit of infinitely many vehicles [1, 10, 11, 20, 21].

This work is mainly concerned with the analysis of traffic flow models at kinetic level. However, we point-out that throughout the paper the study is performed by linking models at different scales. Classical kinetic models for vehicular traffic [26, 27, 37, 40, 41] are based on a Boltzmann-type equation that, in classical kinetic theory, describes the statistical behavior of a system of particles in a gas. Other type of kinetic formulations are available for modeling traffic flow and they also allow to simplify the complex structure of the integro-differential Boltzmann-type equation. Among them, we mention Vlasov-Fokker-Planck type models in which the integral structure of the collision kernel is replaced by differential operators, obtained also by means of suitable time scaling, e.g. see [18, 22], and discrete-velocity models, e.g. see [9, 12, 19].

We study spatially non-homogeneous Boltzmann-type kinetic models for traffic and in particular the BGK approximation, originally introduced by Bhatnagar, Gross and Krook [4] for mesoscopic models of gas particles. The BGK equation replaces the Boltzmann-type kernel with a relaxation towards the equilibrium distribution of the full kinetic equation. For this reason, the BGK model provides a good description of the Boltzmann-type equation in the limit of many interactions, since in this regime the kinetic distribution is close to equilibrium. However, when the interaction frequency is extremely high to immediately drive the kinetic distribution towards the equilibrium one, the kinetic equation converges to the solution of the Lighthill-Whitham-Richards model. Since the flow is at equilibrium, non-equilibrium phenomena are not described in this regime. If the interaction frequency is still large but there is some residual between the kinetic and the equilibrium distribution, then the kinetic model is able to reproduce non-equilibrium waves. We are interested in the study of the stability of the BGK model in the latter case. Using a Chapman-Enskog expansion we show that the BGK equation for traffic flow problems yields an advection-diffusion equation having a negative diffusion coefficient in dense traffic. This analysis is closely related to the investigation of the sub-characteristic condition for relaxation systems [6, 23]. This implies possible unbounded growth of non-equilibrium waves and the corresponding behavior is investigated at the particle description and also at the macroscopic level of the BGK model by looking at the system of the second-order moment equations.

Due to the possible instability of the kinetic model, here we aim to derive a novel BGK equation which is stable or weakly-stable in dense traffic. As in Section 3.2, we define the BGK model being stable if the sub-characteristic condition is always verified and weakly-stable if the sub-characteristic condition is violated in a proper subdomain of the admissible values of the density. As discussed in [46], a weakly-stable model for traffic is not necessarily unrealistic: the instability occurring in the regime can be regarded as model for the appearance of moving unstable but bounded waves, as e.g. stop-and-go waves. The boundedness of these instabilities is numerically observed in Section 3.3. In [46] the analysis of the sub-characteristic condition has been performed for the Aw-Rascle and Zhang model, showing that it could be either satisfied or weakly-satisfied, depending on the choices of the equilibrium speed or flux and of the pressure function. Based on this study, our goal is to develop a BGK-type equation which has as second-order macroscopic structure a Aw-Rascle and Zhang type model so that it may be at least weakly-stable. The new kinetic model for traffic is derived by starting from the follow-the-leader model that can be seen as a discretization of the Aw-Rascle and Zhang model. Therefore, we do not only introduce a well-posed or weakly well-posed kinetic equation but we also contribute in deriving the mesoscopic description of common models for traffic flow. The diagram in Figure 1.1 schematically synthesizes the contribution of this work.

The manuscript is organized as follows. In Section 2 we briefly review the homogeneous kinetic model for traffic flow introduced in [43]. In Section 3 we introduce the non-homogeneous version of the kinetic model and we study non-equilibrium phenomena by using the classical BGK approximation of the full kinetic equation. We show that model is unstable, namely it is able to produce unbounded growth of small perturbations, via Chapman-Enskog expansions and Grad’s moment method. In Section 4 we propose a way to derive an alternative formulation of a BGK-type model for traffic flow. This approach has also the advantage of prescribing the mesoscopic step between the microscopic follow-the-leader model and the macroscopic Aw-Rascle and Zhang model. Finally, we discuss the results and possible perspectives in Section 5.

2 Review of a homogeneous kinetic model

Recent homogeneous kinetic model for traffic flow have been discussed e.g. in [43, 45]. The model proposed therein is used as starting point for the discussion in the following. However, the particular choice below is not required for the shown results later.

General kinetic traffic models of the form read

[TABLE]

where

[TABLE]

is the mass distribution function of the flow, i.e. $f(t,x,v)\mathrm{d}x\mathrm{d}v$ gives the number of vehicles in $[x,x+\mathrm{d}x]$ with velocity in $[v,v+\mathrm{d}v]$ at time $t>0$. Moments of the kinetic distribution function $f$ allow to define the usual macroscopic quantities for traffic flow

[TABLE]

yielding density, flux and mean speed of vehicles, respectively, at time $t$ and position $x$. The space $V:=[0,V_{M}]$ is the space of microscopic speeds. We suppose that $V$ is limited by a maximum speed $V_{M}>0$, which may depend on several factors, and that the density is also limited by a maximum density $\rho_{M}$. Throughout the work we will take dimensionless quantities, thus $V_{M}=1$ and $\rho_{M}=1$.

In [43] the spatially homogeneous kinetic model corresponding to (2.1) has been discussed:

[TABLE]

Both in (2.1) and in (2.4), $\epsilon$ is a positive quantity. It is reasonable to assume that it may depend on the density $\rho$ and possibly its first derivative, see Remark 3.1. The equilibrium solution of (2.4) is called equilibrium distribution, or Maxwellian in analogy with kinetic models for rarefied gas dynamics, and we indicate it by $M_{f}$.

The equilibrium distribution $M_{f}$ depends on the modeling of the collision kernel $Q[f,f]$ that accounts for the kinetic interactions between vehicles. The formulation of $Q[f,f]$ draws inspiration from the classical Boltzmann collision kernel for gas particles. Thus $Q[f,f]$ is split in the difference of a gain term and a loss term. More specifically,

[TABLE]

where we consider binary interactions, in which a test vehicle with velocity $v_{*}$ interacts with a field vehicle with speed $v^{*}$, emerging with speed $v$ as a result of the interaction and with probability given by $\mathcal{P}$. The kinetic model for traffic flow studied in [43] can be generalized to the following probability distribution accounting for stochastic speed transitions, see [45]:

[TABLE]

where $P=P(\rho)\in[0,1]$ is a decreasing function of the density giving the probability of accelerating, $\Delta_{a}$ and $\Delta_{b}$ are the acceleration and the braking parameters, respectively. We can think of $\Delta_{a}$ as the instantaneous physical acceleration of a vehicle. The parameter $\Delta_{b}$ instead corresponds to an uncertainty in the estimate of the field vehicle speed, either one’s own, or the one of the vehicle ahead, as in the first and the second line respectively of (2.6). Note that the model is continuous across the line $v_{*}=v^{*}$ and that mass conservation holds true since

[TABLE]

In the following we will consider $\Delta_{b}=0$ and $\Delta_{a}=\Delta v$ constant, so that one recovers exactly the model in [43].

In [43] a result on the existence and uniqueness of a particular class of stationary solutions has been established:

Theorem 2.1** (Theorem 3.1 and Theorem 3.4 in [43]).**

Let $P=P(\rho)$ be a given function of the density $\rho\in[0,\rho_{M}]$ such that $P\in[0,1]$. Let $\{v_{j}\}_{j=1}^{N}$ be a set of velocities in $[0,V_{M}]$. The distribution function

[TABLE]

with $\sum_{j=1}^{N}f_{j}^{\infty}(\rho)=\rho$, is the unique stable weak stationary solution of the model (2.4)-(2.6) with $\Delta_{b}=0$ provided $v_{j}=v_{1}+j\Delta v$, $j=1,\dots,N$, and

[TABLE]

for $j=1,\dots,N-1$ and $f_{N}^{\infty}=\rho-\sum_{j=1}^{N-1}f_{j}^{\infty}$.

Modeling with kinetic equations such as (2.1) is driven by the interaction kernel $Q[f,f]$ and for this reason homogeneous kinetic models are widely studied for traffic flow problems. Using the steady-state of the kinetic model it is possible to define the flux and the mean speed of vehicles at equilibrium as

[TABLE]

Notice that Theorem 2.1 states that $M_{f}$ depends only on a finite number of speeds, that it is parameterized by the local density $\rho$ and that, for each fixed value of $\rho$, it is unique. Therefore, the maps $\rho\mapsto Q_{\text{eq}}(\rho)$ and $\rho\mapsto U_{\text{eq}}(\rho)$ are actually functions prescribing a uniquely defined correspondence between the density and the quantities (2.7). These relations provide the so-called simulated fundamental diagrams of traffic which are used in order to validate the model.

In Figure 2.1 we show the fundamental diagrams of homogeneous kinetic model (2.4)-(2.6) that are computed by using the equilibrium distribution prescribed by Theorem 2.1. The top row shows the correspondence $\rho\mapsto Q_{\text{eq}}(\rho)$ and the bottom row shows $\rho\mapsto U_{\text{eq}}(\rho)$. We take $P(\rho)=1-\rho$ and we consider the case of $N=2$ speeds (left panels) and $N=3$ speeds (right panels) with $\Delta v=\frac{1}{N-1}$. In Figure 2.2 we show that the fundamental diagrams can fit experimental flux diagrams. This is for instance the case of the experimental measurements provided by the Minnesota Department of Transportation in 2013 and tuning the parameters of the model as $N=3$ speeds, $\Delta v=\frac{1}{2}$ and $P(\rho)=1-\rho^{\frac{1}{4}}$. Observe that the flux at maximum density is not zero since we compute the diagram up to a value of the density $\widetilde{\rho}$ such that $\widetilde{\rho}<\rho_{M}$. This is due to the observation that experimental data contain a residual movement even in the congested phase. The model is able to describe the phase transition and to catch the two regimes of traffic compared with experimental data. For low values of $\rho$, more precisely for $\rho<\rho_{\mbox{\scriptsize crit}}$, where $\rho_{\mbox{\scriptsize crit}}$ is the solution of $P(\rho)=\frac{1}{2}$, the flux is linear in $\rho$. This is the phase of free flow, when $U_{\text{eq}}(\rho)=V_{M}$. For larger values of $\rho$, i.e. $\rho$ larger than the critical density $\rho_{\mbox{\scriptsize crit}}$, the role of the interactions increases, and the flux decreases. This corresponds to the congested phase of traffic flow, see also [24].

Remark 2.2**.**

Thanks to Theorem 2.1 we have that the equilibrium distribution $M_{f}$ has a fixed dependence on $v$, and depends on $x$ and $t$ only through the local density $\rho(x,t)$. However, for $\Delta_{b}=0$, it was shown in [43] that unstable equilibria can be found, for which $M_{f}$ depends not only on $\rho$, but also on the initial distribution $f(t=0,x,v)$. Those equilibria are unstable, in the sense that they disappear and drive to the solutions given by Theorem 2.1 when the initial data are perturbed. If the braking uncertainty $\Delta_{b}>0$, it was shown in [45] that only stable equilibria are found, they still depend on a finite number of speeds. In this work we still focus on the case $\Delta_{b}=0$ and consider initial distributions $f(t=0,x,v)$ such that only stable equilibria occur.

3 Uncontrolled instabilities in standard non-homogeneous kinetic models

3.1 Motivation: non-equilibrium regimes

The properties of the fundamental diagrams derived in Section 2 are the result of the modeling of the microscopic interactions. Although the model is able to reproduce typical features of the experimental diagrams, we point-out that the simulated diagrams are only single valued. This is a direct consequence of the uniqueness of the equilibrium distribution provided by Theorem 2.1. This means that the model cannot reproduce and explain the scattering behavior of the experimental measurements, at least not at equilibrium. Mathematical models being able to describe multivalued fundamental diagrams were also studied and they link the observed scattered data to several reasons: for instance the presence of heterogeneous components, such as different drivers [31, 33] or populations of vehicles [42], or the deviation of microscopic speeds at equilibrium [13], or the presence of non-equilibrium phenomena such as stop and go waves [25, 46, 49].

In this work we will focus on the last motivation and consider non-equilibrium phenomena that might lead to complex traffic phenomena as e.g. stop and go waves. Such phenomena cannot be described by spatially homogeneous models but there is the need to reconsider the space dependence and thus to take into account the full kinetic equation (2.1).

3.2 Non-homogeneous kinetic models

In the non-homogeneous case, the relaxation parameter $\epsilon$ plays a crucial role since it balances the weight between the convection and the source term. In analogy to gas dynamics, where the interaction frequency is given by the Knudsen number, different values of $\epsilon$ lead to different regimes. The situation is synthesized in Table 3.1. If we allow for $\epsilon=0$, i.e. we suppose that the interactions are so frequent to instantaneously relax $f$ to the local equilibrium distribution $M_{f}$, we are in the so-called equilibrium flow regime where (2.1) is solving the conservation law for the density

[TABLE]

with $Q_{\text{eq}}(\rho)$ defined in (2.7). Equation (3.1) is a first-order macroscopic model for traffic flow  [30, 44]. Instead, we expect that if $\epsilon$ is small, but not vanishing, then we are in a regime where the kinetic equation (2.1) is solving the continuity equation (3.1) with a diffusive perturbation of order $O(\epsilon)$. For $\epsilon\asymp 1$ we are in the transitional regime where the counterpart of the kinetic equation is given by extended continuum hydrodynamic equations. In traffic flow, this is the case of the Aw-Rascle [2] and Zhang [48] model. For $\epsilon>1$, but not too large, we are in the kinetic regime or rarefied gas regime and finally for $\epsilon\to\infty$ we obtain the regime of the collision-less kinetic equation.

Here, we will consider a type of approximation for the full non-linear Boltzmann-type equation which formally holds for small values of $\epsilon$, namely the BGK approximation. This technique is based on replacing the Boltzmann-type kernel in (2.1) by a relaxation term:

[TABLE]

where $M_{f}(v;\rho)$ is the distribution of the microscopic speeds $v$ at equilibrium, corresponding to the density $\rho$ at time $t$ and position $x$ derived from the spatially homogeneous equation (2.4). Equation (3.2) assumes that the equilibrium distribution is known. For small values of $\epsilon$, the BGK model provides the same behavior of the full kinetic equation (2.1). In contrast, (3.2) is only an approximation of (2.1) in the regime with $\epsilon$ large.

Remark 3.1** (Backward propagating information).**

Let us consider $\epsilon$ fixed and $\epsilon\gg 1$, i.e. the regime in which the transport term in (3.2) is stronger then the relaxation term. It is well known that in this regime, since the velocities in traffic are non-negative, the information propagate in the wrong direction for congested traffic situations. In this phase, in fact, we expect to see backward propagation of signals but the transport part move them forward. This phenomenon is shared both by a full kinetic model (2.1) and the BGK equation (3.2). For this reason, several models with the ability of overcoming this drawback were introduced in the mathematical literature [12, 27]. This effect could be linked to a fixed choice of $\epsilon$. Instead, in analogy with the Knudsen number in gas models, $\epsilon$ should be a decreasing function of the density. In this way, $\epsilon$ is large in the free flow phase since the interactions are less frequent and the convective term rules the dynamics. On the contrary, since $\epsilon$ decreases when density increases, i.e. interactions among vehicles are dominant, we expect that relaxation is faster when $\rho$ is high. It is clear from the fundamental diagrams that, in the regimes in which $f$ is close to the local equilibrium $M_{f}$, signals move backward in congested phases. See also the left panel of Figure 3.1 where we show the characteristic speeds of the kinetic model when $f$ is close to equilibrium for several values of the speed number $N$. In other words we claim that the regimes described in Table 3.1 depend on the phase of traffic.

Numerical scheme.

We briefly describe the numerical discretization of BGK using a first-order implicit-explicit (IMEX) scheme in time [34, 39] by treating the convective term explicitly and the collision term implicitly. This allows to treat the stiffness for small values of $\epsilon$. Moreover, observe that, due to the linearity of the source term with respect to $f$, the scheme remains still explicit.

Let us consider the computational domain $(t,x,v)\in[0,T_{M}]\times[0,L]\times[0,1]$ where $T_{M}$ is the final time and introduce a mesh $I_{ij}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]\times\{v_{j}\}$, where $x_{i}=(i-\frac{1}{2})\delta x$, $i=1,\dots,N_{x}$ and $v_{j}=(j-1)\delta v$, $j=1,\dots,N_{v}$. Here, the numerical parameter $\delta v$ is a sub-multiple of the physical parameter $\Delta v$ defined in Section 2. In this way, $N$, i.e. the number of discrete speeds defining the equilibrium distribution $M_{f}$, see Theorem 2.1, is a sub-multiple of $N_{v}$. Although in the transient speeds that are not spaced by $\Delta v$ give a contribution to the dynamics, while at equilibrium only the $N$ speeds spaced by $\Delta v$ give a non-zero contribution, here we always consider the simple case with $\delta v=\Delta v$. For piece-wise constant approximation of cell-averages in space we have $f_{ij}^{n}\approx f(n\delta t,x_{i},v_{j})$ and obtain

[TABLE]

for $i=1,\dots,N_{x}$, $j=1,\dots,N_{v}$ and $n\geq 0$. Observe that the explicit knowledge of the equilibrium distribution $M_{f}$, given by Theorem 2.1, is crucial for the performance of the numerical scheme. The quantity $\mathcal{F}$ can be one of the known Lipschitz continuous, conservative and consistent numerical flux function. The asymptotic preserving [14, 35] result holds for the scheme (3.3).

Proposition 3.2** (Asymptotic preserving).**

Scheme (3.3) is asymptotic preserving in the sense that it solves a global first-order discretization of the macroscopic equation (3.1) with the flux function $Q_{\text{eq}}(\rho)$ given by (2.7) when $\epsilon\to 0$.

Proof.

By equation (3.3a), we have that $f_{ij}^{(1)}\to M_{f}(v_{j};\rho_{i}^{n})$, if $\epsilon\to 0$. Substituting $f_{ij}^{(1)}$ in (3.3b) then, for $\epsilon\to 0$, we find

[TABLE]

and hence

[TABLE]

This scheme is an explicit conservative, consistent and stable discretization of the macroscopic equation (3.1) with the flux function $Q_{\text{eq}}(\rho)$ given by (2.7), provided $\mathcal{F}(\cdot,\cdot)$ is a conservative, consistent and stable numerical flux function for the non-linear conservation law (3.1). ∎

Hence, we choose the time-step $\delta t$ using the CFL condition applied to the non-linear conservation law (3.1) and in view of the nonlinear transport on macroscopic level we employ the Lax-Friedrichs flux:

[TABLE]

Non-equilibrium effects.

We show numerically that the multivaluedness in fundamental diagrams can be reproduced as a result of non-equilibrium effects.

The top rows of Figure 3.2, 3.3 and 3.4 contain the solution obtained with $\epsilon=0$, i.e. in the equilibrium regime, on the left, and the solution found with $\epsilon=10^{-1}$ on the right. The solution is shown in several instants of time, starting from the green curve and ending with the blue thick profile. The final solutions are compared in the bottom left part of the figure: the profile with $\epsilon=10^{-1}$ is depicted with the continuous red line, and the profile with $\epsilon=0$ with the dashed black line.

The dots in the right bottom panel contain the evolution of the flux $\rho u$ defined in (2.3) which is drawn as a function of $\rho$. Note that, at equilibrium, the fundamental diagram $Q_{\text{eq}}(\rho)$ is a function of $\rho$. When the system is not in equilibrium, the flux $\rho u$ depends not only on $\rho$, but also on the distribution $f$. Again, the green dots are obtained in the initial stages of the simulation, and they fade into blue as the solution evolves. If the flow were in equilibrium, as in the case $\epsilon=0$, then $f$ coincides with the equilibrium distribution $M_{f}$ (uniquely determined by $\rho$). Thus the dispersion in these plots is a measure of deviations from equilibrium. Hence, the dispersion of the points in blue may reflect the multivaluedness of the fundamental diagram

In all simulations we consider $N=N_{v}=3$ speeds and so the physical acceleration is $\Delta v=\frac{1}{2}$. The spatial domain is $[0,1]$ and discretized by $N_{x}=400$ cells. The final time is $T_{M}=0.2$. We consider an initial smooth perturbation in the density

[TABLE]

with periodic boundary conditions. The initial kinetic distribution $f(t=0,x,v)$ is the uniform distribution over the velocity space with density $\rho_{0}(x)$.

In Figure 3.2, $\rho_{\min}=0.1$ and $\rho_{\max}=0.3$. The perturbation in the density is below the critical density $\rho_{\text{crit}}=\frac{1}{2}$. Thus here the whole solution moves towards the right for $\epsilon=0$ because in this regime the flux is linear with respect to the density. In the case $\epsilon=10^{-1}$ the solution also moves to the right but with very little deformation. As consequence, the flux diagram shows a very small dispersion, which means that the flow is almost at equilibrium.

If we choose $\rho_{\min}=0.6$ and $\rho_{\max}=0.9$, the initial perturbation occurs on a dense traffic flow, see Figure 3.3. Now, the expected propagation speed in the fundamental diagram is negative, and in fact the solution for $\epsilon=0$ moves towards the left. The solution with $\epsilon=10^{-1}$ also moves towards the left, but new maxima and minima appear which means that the maximum principle is no longer satisfied. These are purely non-equilibrium effects, as shown also in the flux diagram. The equilibrium fundamental diagram is shown with the red line, and it can be seen bisecting the $\rho u$ values computed during the evolution of the case $\epsilon=10^{-1}$.

The last smooth test we show has a large initial pulse, with $\rho_{\min}=0.1$ and $\rho_{\max}=0.9$. The results are in Figure 3.4. The low density part of the wave on the right accelerates and propagates towards the right. The high density part of the wave propagates backward, increasing the slope of the profile. Thus a shock forms, propagating towards the left. This structure appears quickly in the case $\epsilon=0$ while it is slower for $\epsilon=10^{-1}$ where however we see the appearance of a new maximum. Again the solution with $\epsilon=10^{-1}$ exhibits non equilibrium waves which are represented by the scattered values in the flux diagram.

Remark 3.3** (Solutions of Riemann problems).**

As already observed, in the relaxation limit the full kinetic equation (2.1) and the BGK approximation (3.2) solve the conservation law (3.1) with $Q_{\text{eq}}$ computed using the equilibrium distribution of the homogeneous kinetic model. We notice that the fundamental relation $Q_{\text{eq}}$ defined in (2.7) has two main properties.

It is continuous: in fact, from Theorem 2.1 it is clear that

[TABLE]

Thus, there is no need to treat the macroscopic model (3.1) with special techniques for conservation laws with discontinuous fluxes. 2. 2.

It is not convex: in fact, it is linear in the free flow regime (i.e. when $\rho\leq\rho_{\text{crit}}$) and it is convex in the congested regime (i.e. when $\rho>\rho_{\text{crit}}$). See the right panel in Figure 3.1 which shows the behavior of the second derivative of $Q_{\text{eq}}$ for different values of the speed number.

The non-convexity of the flux-density diagram represents the most interesting property since it implies that the wave structure at equilibrium is non-trivial, because even a simple Riemann problem initial data will result in more than a single wave. For a more detailed discussion we refer to [5, 29].

Unbounded instabilities.

The results of the previous simulations show that non-homogeneous models are important to describe non-equilibrium phenomena. Note that, in the case of the perturbation in dense traffic, the non-equilibrium effects provides new maxima and minima that apparently grow in time. These results suggest that the appearance of stop and go waves on highways could be interpreted and reproduced as non-equilibrium effects.

Here, we prove that the BGK model (3.2) has the drawback that it could lead to unbounded unstable waves in the congested regime of traffic. In fact, let us consider the same setting of the simulations provided in Figure 3.3, namely a perturbation of the density in the congested regime. In Figure 3.5 we show the solutions for a larger final time with $\epsilon=10^{-1}$ in the left panel and $\epsilon=10^{-2}$ in the right panel. We stop the simulations as soon as the density profile exceeds the maximum density. We observe that the BGK model propagates the initial profile backward producing unstable waves that do not remain bounded. This phenomenon will be investigated below using the Chapman-Enskog expansion.

Consider the Chapman-Enskog expansion of the BGK equation (3.2) as first order perturbation in $\epsilon$ of $f$ around the equilibrium distribution $M_{f}(v;\rho)$. This leads to the advection-diffusion equation

[TABLE]

where the diffusion coefficient $\mu(\rho)$ is given by

[TABLE]

In fact, considering the expansion

[TABLE]

observing that

[TABLE]

plugging the approximation of $f$ into the BGK equation (3.2) and integrating with respect to the velocity we obtain

[TABLE]

If $\mu(\rho)<0$ then the advection-diffusion equation is ill-posed and therefore has solutions with unbounded growth even starting from small perturbations. In the case of the kinetic model (3.2), the sign of the diffusion coefficient depends on the equilibrium distribution $M_{f}$. The request $\mu(\rho)>0$ is

[TABLE]

and, since $Q_{\text{eq}}(\rho)$ is the fundamental diagram at equilibrium, this condition requires that the square of the characteristic velocities is bounded by the variation of the kinetic energy in each regime. The condition (3.7) depends only on the equilibrium distribution $M_{f}$ which is given by the homogeneous kinetic model and is the result of the modeling of the microscopic interactions.

In Figure 3.6 we show the sign of the diffusion coefficient (3.6) in the case of the equilibrium distribution computed in [43], and recalled in Section 2, with two and three speeds corresponding to $\Delta v=1$ and $\Delta v=\frac{1}{2}$, respectively. We observe that $\mu(\rho)\geq 0$ if $\rho\leq\rho_{\text{crit}}=\frac{1}{2}$, while $\mu(\rho)<0$ if $\rho>\rho_{\text{crit}}=\frac{1}{2}$. This result could explain the unbounded growth of perturbations in the density, numerically observed in Figure 3.5.

In the following result, we prove that the instability of the model does not depend on the choice of the equilibrium distribution but it occurs for any suitable equilibrium of kinetic traffic models.

Proposition 3.4**.**

Assume that $\exists\,\widetilde{\rho}\in(0,1)$ such that

[TABLE]

for all $\rho\in(\widetilde{\rho},1)$. Then the quantity $\mu(\rho)$ given in (3.6) is negative $\forall\,\rho\in(\widetilde{\rho},1)$.

Proof.

Conditions (3.8) are sufficient to $\int_{V}v^{2}\partial_{\rho}M_{f}(v;\rho)\mathrm{d}v<0$. In fact

[TABLE]

$\forall\,\rho\in(\widetilde{\rho},1)$. ∎

More specifically, for the case of two microscopic speeds and in view of the previous proposition, the model can be seen as the kinetic relaxation system of Jin and Xin [23] that cannot satisfy the sub-characteristic condition since $v>0$ and $Q_{\text{eq}}^{\prime}(\rho)<0$ in the congested regime.

Remark 3.5**.**

We observe that conditions (3.8) do occur in traffic flow. In particular, the first one, $Q^{\prime}_{\text{eq}}(\rho)<0$, is always verified in the congested phase of traffic. The second condition instead requires that the variance of the microscopic speeds decreases as the density increases. This is a reasonable assumption in traffic since we expect that the freedom in choosing a microscopic velocity reduces when traffic becomes dense.

In view of the results provided by the Chapman-Enskog analysis, we state the following definition which establishes that the BGK model is unstable.

Definition 3.6**.**

A kinetic model is said to be stable if $\mu(\rho)\geq 0$, $\forall\,\rho\in[0,\rho_{M}]$, weakly-unstable if $\mu(\rho)<0$ on an interval $(\rho_{1},\rho_{2})$ properly contained in $[0,\rho_{M}]$ and unstable if $\mu(\rho)<0$ on an interval $(\rho_{1},\rho_{2})$ in which either $\rho_{1}=0$ or $\rho_{2}=\rho_{M}$.

In order to further investigate the instability of the BGK model when applied to traffic flow problems, we must go beyond its hydrodynamic limit. For this reason, we consider higher order macroscopic moments, following the philosophy of Grad’s moment method in gas dynamics. Therefore, the second-order macroscopic model derived from the BGK equation (3.2) is obtained without closing the continuity equation with the distribution at equilibrium and computing the evolution equation for the first moment. Then we have

[TABLE]

System (3.9) is closed and thus solvable when a closure law for the second moment of the kinetic distribution $f$ is prescribed.

If we consider the following approximation

[TABLE]

with $h(\rho)$ being an increasing function of the density, then from (3.9) it is possible to recover the Payne and Whitham model [38] second-order model for traffic flow:

[TABLE]

It is well-know that model (3.10) suffers of the following drawback (see for instance the works by Daganzo [8]): the characteristic velocities of (3.10) are $\lambda_{1}=u-\sqrt{h^{\prime}(\rho)}$ and $\lambda_{2}=u+\sqrt{h^{\prime}(\rho)}$. Thus one characteristic speed is always larger than the speed $u$ of vehicles. 2. 2.

If we instead approximate

[TABLE]

then (3.9) reads

[TABLE]

that is again the Payne and Whitham model without the anticipation term. However, the model still suffers of the following drawbacks. First, the system is not strictly hyperbolic since the eigenvalues are $\lambda_{1}\equiv\lambda_{2}=u$ leading to a lack of strict hyperbolicity. The possible unbounded growth of small perturbations is highlighted by the fact that the sub-characteristic condition for the model is $-\rho^{2}U_{\text{eq}}^{\prime}(\rho)>0$ which is never satisfied. 3. 3.

It is easy to show that if we approximate around the equilibrium distribution, i.e.

[TABLE]

the relaxation system (3.9) has a sub-characteristic condition $\mu(\rho)<0$ with $\mu(\rho)$ given by (3.6) and therefore in the whole congested regime.

Remark 3.7**.**

The particle dynamics resulting from the BGK model (3.2) are described by the ODE system

[TABLE]

where $i$ is the label of a vehicle. System (3.12) represents the microscopic model by Bando [3] which is know to be stable provided $\epsilon<\frac{1}{2|U_{\text{eq}}^{\prime}(\rho)|}$, i.e. if the relaxation is large.

The last results of this paragraph show that properties of the BGK model studied via Chapman-Enskog expansion occurs also in its second-order systems of moments and in its particle formulations.

Closing BGK at second order moments one provides the structure of common second-order models in the literature on traffic flow. Due to this consideration we aim to understand what is the kinetic structure representing the mesoscopic scale of suitable second-order macroscopic models.

In the following section we investigate properties of the second-order macroscopic model [2, 48, 46].

3.3 The case of the Aw-Rascle and Zhang model

The prototype example of second-order macroscopic traffic models is given by the Aw and Rascle [2] and Zhang [48] (ARZ) model where the conservation law for the density $\rho(t,x)$ of vehicles at time $t\geq 0$ and position $x$ is coupled with an additional equation for the evolution of the macroscopic speed of the flow $u(t,x)$. We focus on the non-homogeneous version of the ARZ model which writes in conservation form

[TABLE]

where $q$ is the second conservative variable and it is defined as $q:=\rho(u+h(\rho))$. The function $h=h(\rho)$ is a strictly increasing function of the density and it is called hesitation function or traffic pressure although it is homogeneous to a velocity. The quantity $\epsilon$ is a time which rules the relaxation speed of the flux $\rho u$ to the equilibrium flux $Q_{\text{eq}}(\rho)$ which is a given function of the density, i.e. the fundamental diagram. Here $Q_{\text{eq}}$ is not necessarily the one given in (2.7).

The ARZ model was introduced in order to overcome the drawbacks of the second-order macroscopic model introduced by Payne and Whitham and recalled in equation (3.10). The characteristic velocity of the ARZ model are $\lambda_{1}=u-\rho h^{\prime}(\rho)\leq u=\lambda_{2}$ and thus information propagates at most with the vehicles’ velocity.

In primitive variables the system reads

[TABLE]

In (3.14) $U_{\text{eq}}(\rho)=\frac{1}{\rho}Q_{\text{eq}}(\rho)$ is the equilibrium speed.

Both (3.13) and (3.14) can be understood as relaxation systems [23]. In the limit $\epsilon\to 0$, model (3.14) converges towards solutions of the conservation law (3.1). In contrast to second-order macroscopic models, equation (3.1) has only one characteristic velocity given by

[TABLE]

If $\epsilon$ is small, but not vanishing, (3.14) approaches the advection-diffusion equation (3.5) where the diffusion coefficient $\mu(\rho)$ is given by

[TABLE]

In fact, using a Chapman-Enskog expansion, the coefficient (3.15) is obtained by considering a first-order expansion of the speed $u=U_{\text{eq}}(\rho)+\epsilon u_{1}$ around the equilibrium velocity function $U_{\text{eq}}(\rho)$

[TABLE]

Finally, plugging the expansion of $u$ just obtained into the first equation of (3.14) we get the advection-diffusion equation (3.5) with coefficient $\mu(\rho)$ given by (3.15). The condition $\mu(\rho)>0$ provides the so-called sub-characteristic condition [6, 23]. For the ARZ model $\mu(\rho)>0$ is satisfied if

[TABLE]

The sub-characteristic condition guarantees the stability of solutions of the ARZ model based on initial conditions where $\rho(t,x)=\overline{\rho}$ and $\overline{u}=U_{\text{eq}}(\overline{\rho})$. In fact, the LWR model (3.1) and the ARZ model (3.14) share solutions of the type $(\overline{\rho},\overline{u})$ and in this case the sub-characteristic condition (3.16) can be rewritten in terms of the characteristic velocities of the two models computed in $(\overline{\rho},\overline{u})$ as

[TABLE]

The above condition implies that the characteristic velocity of the LWR model must lie between the two characteristic velocities of the ARZ model. Whitham [47] proved that this inequality is verified if and only if solutions of the second-order macroscopic model with initial condition $(\overline{\rho},\overline{u})$ are linearly stable, i.e. small perturbations to $(\overline{\rho},\overline{u})$ decay in time [6].

We stress the fact that condition (3.16) strongly restricts the possible choice of $U_{\text{eq}}$ and $h$. For instance, the Greenshields’ closure $U_{\text{eq}}(\rho)=1-\rho$ and the hesitation function $h(\rho)=C\rho^{\gamma}$, $C,\gamma>0$ given in [2] violate (3.16) if $\rho^{\gamma-1}\leq\frac{1}{C\gamma}$.

In Figure 3.7 we numerically analyze the sign of the diffusion coefficient $\mu(\rho)$ of the ARZ model given in (3.15) for several choices of the equilibrium velocity and the hesitation function. We observe that standard choices lead to regimes where the diffusion coefficient is negative in the congested regime and consequently the sub-characteristic condition is not verified. This produces instabilities of uniform flow as discussed above. However, for suitable choices of $U_{\text{eq}}$ and $h$, regimes where $\mu(\rho)<0$ can be properly contained in $[0,\rho_{M}]$, and in this case numerical evidence shows the existence of a bounded density regime.

This result was already mentioned and analyzed in [46] where the authors link the appearance of such regimes to the presence of nonlinear traveling wave solutions, called jamitons, and identified as the stop and go waves observed in real traffic situations.

Bounded waves.

Numerical evidence suggests that if the instability of the model is restricted on a set $(\rho_{1},\rho_{2})$ properly contained in $[0,\rho_{M}]$ the amplitude of unstable waves remains bounded. In fact, the growth of a wave is smoothed when the solution enters in the regime $[0,\rho_{M}]\setminus(\rho_{1},\rho_{2})$ where the diffusion is positive.

We point-out that this consideration is purely heuristic and based on numerical investigations. To the best of our knowledge there are few contributions on the analysis of advection-diffusion equations with a changing sign diffusion coefficient, see e.g. [7, 16, 28, 32].

In Figure 3.8 we numerically study the evolution of a density perturbation with the ARZ model as we already did for the BGK model in the previous section. Here we consider the case in which the bump in the density intersects a regime where the diffusion is negative. In particular, we consider the equilibrium flux $Q_{\text{eq}}$ and the equilibrium velocity $U_{\text{eq}}$ prescribed by the three-velocity kinetic distribution at equilibrium, see (2.7). The pressure or hesitation function is chosen as $h(\rho)=2\rho$. This is the situation described in Figure 7(f) where $\mu(\rho)<0$ for $\rho\in(0.5,0.66)$. The initial condition is given by (3.4) with $\rho_{\min}=0.55$ and $\rho_{\max}=0.9$. The space domain is $[-1,1]$ and discretized with $N_{x}=400$ cells. The top-left panel shows the density profile at different instants in time. We observe that unstable waves appear in the regime where the diffusion coefficient is negative and they move backward. The top-right panel shows a zoom in the region where the unstable waves are rising and the bottom-left panels shows the evolution of the macroscopic speed. The growth is bounded as confirmed in the bottom-right panel where we show the density profile for a longer final time and using periodic boundary conditions.

In the following, we consider models that allow for a Chapman-Enskog expansion as advection-diffusion equation with a negative diffusion coefficient in a subset $(\rho_{1},\rho_{2})$ properly contained in $[0,\rho_{M}]$. The violation of the sub-characteristic condition lead to instabilities that in turn can be regarded as models for stop and go waves. These waves are however unbounded if the regime of negative diffusion is not properly contained in $[0,\rho_{M}]$ as in the case of the classical BGK model.

4 Derivation of a modified BGK-type model for traffic flow

In this section we address the issue of deriving a modified BGK-type equation for traffic flow satisfying the following properties.

Property 1**.**

A kinetic model for traffic flow should be weakly-unstable in the sense of Definition 3.6.

Property 2**.**

The second-order system of moments of a kinetic model for traffic flow should recover an ARZ-type macroscopic model.

The starting point for the derivation of the BGK-type model are microscopic follow-the-leader (FTL) models. The FTL model is proved to converge to the ARZ model in the macroscopic limit [1].

The kinetic model derived with this approach will also automatically guarantee Property 2. Moreover, since we have shown that an ARZ-type model may have a negative diffusion coefficient in a small density regime, see for instance Figure 7(d)-7(f), we expect that the BGK-type equation will also satisfy Property 1.

We will also observe that the ARZ-type second-order macroscopic model derived by the new mesoscopic structure is rather general since it keeps the details that kinetic provides, and it allows to recover the classic ARZ model as special case.

4.1 Microscopic follow-the-leader model

In [1] a first-order semi-discretization in space of the Lagrangian formulation of the ARZ model is shown to be equivalent to the microscopic follow-the-leader model:

[TABLE]

with

[TABLE]

and where $x_{i}$ and $v_{i}$ are the microscopic position and velocity of the vehicle $i$ at time $t\in\mathbb{R}^{+}$, $U_{\text{eq}}(\rho_{i})$ is a target speed depending on the value of the local density $\rho_{i}=\frac{\Delta X}{x_{i+1}-x_{i}}$. The quantity $\Delta X$ is the typical length of a vehicle. In the following computations we assume $\Delta X$ normalized to $1$. The constants $C_{\gamma}>0$, $\gamma>0$ and the relaxation time $\epsilon$ are given parameters. The link between the microscopic and the macroscopic model is also studied in [10, 11, 17].

Following [1], we consider an alternative but equivalent formulation of the model (4.1) based on the quantity $w_{i}:=v_{i}+p(\rho_{i})$. The function $p=p(\rho_{i})$ is the so-called traffic pressure and satisfies $p(\rho)\geq 0$, $p^{\prime}(\rho)>0$. In terms of the new variable $w_{i}$, model (4.1) reads

[TABLE]

provided $p(\rho)$ is chosen as

[TABLE]

with $K$ given by (4.2). Model (4.3) is identical to (4.1) but the introduction of $w_{i}$ allows us to rewrite the acceleration equation as a relaxation step. Note that the choice of the function $p(\rho)$ is determined by the modeling of $K$. Since a different interaction term leads to different pressure functions, we consider the case of a general function $K$.

Next we introduce stochastic particle interactions in terms of the desired speed $w_{i}$ that is an approximation of (4.3). Then, we derive the corresponding kinetic equation for the evolution of the distribution of the desired velocity via the relaxation algorithm, see [36, Section 4.2.2]. Finally, we rewrite the kinetic equation in terms of the actual microscopic speeds of vehicles and we verify that Property 1 and 2 hold true for the new kinetic model for traffic.

4.2 A BGK-type model for traffic

Let $g=g(t,x,w):\mathbb{R}^{+}\times\mathbb{R}\times W\to\mathbb{R}^{+}$ be the kinetic distribution function such that $g(t,x,w)\mathrm{d}x\mathrm{d}w$ gives the number of vehicles at time $t\in\mathbb{R}^{+}$ with position in $[x,x+\mathrm{d}x]\subset\mathbb{R}$ and desired speed in $[w,w+\mathrm{d}w]\subset W$. The desired speed $w$ is assumed to be the speed that drivers want to keep in “optimal” situations. We define $W:=[w_{\min},+\infty)$ the space of the microscopic desired speeds where $w_{\min}>0$ may be interpreted as the minimum speed limit in free-flow conditions.

The macroscopic density, i.e. the number of vehicles per unit length, at time $t$ and position $x$ is defined by

[TABLE]

and we define the macroscopic quantity

[TABLE]

In order to derive the evolution equation for the kinetic distribution $g=g(t,x,w)$, we first prescribe the update of the microscopic states in a probabilistic interpretation that has a close link with the microscopic particle model (4.3). Let us assume a set of $N$ particles $(x_{i}^{0},w_{i}^{0})$ is given at time $t_{0}$ characterized by the position $x_{i}^{0}$ and the desired speed $w_{i}^{0}$. We can think of $(x_{i}^{0},w_{i}^{0})$ as samples from the distribution $g(t=0,x,w)$. The update of the microscopic states $(x_{i},w_{i})$ at time $t_{0}+\Delta t$ is generated by the following algorithm.

Step 1.

Reconstruct the density $\rho_{i}=\frac{1}{x_{i+1}-x_{i}}$ at time $t$

Step 2.

Compute $x_{i}=x_{i}^{0}+(w_{i}^{0}-p(\rho_{i}))\Delta t$

Step 3.

For each $i$, with probability $\frac{\Delta t}{\epsilon}$ replace $w_{i}^{0}$ with a velocity sample from a distribution $M_{g}$ with density $\rho_{i}$; otherwise set $w_{i}=w_{i}^{0}$.

So far, the distribution $M_{g}=M_{g}(w;\rho)$ only has to fulfill the requirement

[TABLE]

However, we additionally require that

[TABLE]

In fact, we observe that under this assumption the probabilistic microscopic algorithm prescribes an average speed which is equal to the output velocity prescribed by a first-order time discretization of (4.3).

According to [36, Section 4.2.2], we interpret the interactions introduced above as a stochastic particle method solving

[TABLE]

This equation is still a BGK-type equation since the collision kernel is linear and describes the relaxation of $g$ towards a given distribution $M_{g}$ parameterized by the density $\rho$.

It is important to point-out that, compared to classic kinetic theory, our approach is different in the sense that $M_{g}$ is an “equilibrium distribution” with a modified idea of microscopic velocity. Thanks to (4.6), we impose a-priori $M_{g}$ but it is still based on the knowledge of the classical Maxwellian $M_{f}$, which is related to the classical concept of microscopic velocity, by means of $U_{\text{eq}}(\rho):=\frac{1}{\rho}\int vM_{f}\mathrm{d}v$. We do not impose $M_{f}$ a-priori (and so $U_{\text{eq}}(\rho)$ and consequently $M_{g}$), but we use the equilibrium distribution $M_{f}$ that comes out from the modeling of microscopic interactions of the spatially homogeneous kinetic model in Section 2. Therefore, the added details provided by a kinetic approach are still present.

Moreover, we observe that Maxwellian $M_{f}$ of other kinetic models can be employed, allowing to define $M_{g}$ and the BGK model (4.7) with the modified concept of microscopic speeds. Before concluding this section, we further investigate relation between $M_{f}$ and $M_{g}$.

Let $f=f(t,x,v)$ be the kinetic distribution on the microscopic speed $v$ of vehicles, such that $f(t,x,v)\mathrm{d}x\mathrm{d}v$ gives the number of vehicles at time $t$ with position in $[x,x+\mathrm{d}x]\subset\mathbb{R}$ and speed in $[v,v+\mathrm{d}v]\subset V:=[0,+\infty)$. From (4.7) it is possible to derive the equation for $f$ via a nonlinear change of variables.

[TABLE]

The deviation is controlled by $p(\rho)$ in such a way that $p(0)=0$ so that $v=w$, i.e. drivers can behave as they wish in the vacuum regime, and $p(\rho_{M})\leq w_{\min}$ in order to ensure the bound $v\geq 0$. In terms of $f=f(t,x,v)$ we define the density as

[TABLE]

and the flux of vehicles as

[TABLE]

Further, $M_{g}(w;\rho)=M_{g}(v+p(\rho);\rho),\ w\in W,\ v\in V$ and $M_{f}=M_{f}(v;\rho)$, and setting

[TABLE]

relation (4.6) is satisfied if $U_{\text{eq}}(\rho)$ is the equilibrium velocity provided by the Maxwellian $M_{f}$ of a kinetic model for traffic. In fact, for each fixed density $\rho$, we have that

[TABLE]

and

[TABLE]

Hence, $\int_{V}vM_{f}(v;\rho)\mathrm{d}v=\rho U_{\text{eq}}(\rho).$

4.3 Properties of the BGK-type model

We verify that Property 1 and Property 2 hold true for the kinetic equation (4.7).

Property 2 - Macroscopic limits.

We show that Property 2 holds by computing the system of moment equations of (4.7), up to second order.

The continuity equation is obtained by integration in $w$ space as

[TABLE]

where $q$ is defined by (4.5). If one assumes that $g$ is at equilibrium , $g=M_{g}$, we obtain the classical first-order model for traffic flow

[TABLE]

Instead, the second-order moment equation gives the system

[TABLE]

System (4.13) is not closed since the knowledge of the second moment of $g$ is required. In particular, the ARZ model (3.13) is recovered by closing as

[TABLE]

and taking the function $p(\rho)=h(\rho)$.

Property 1 - Chapman-Enskog expansion.

In order to verify that Property 1 holds true for the kinetic equation (4.7) we use Chapman-Enskog expansion approximating

[TABLE]

and defining $Q_{\text{eq}}(\rho)=\rho U_{\text{eq}}(\rho)$. Following similar computations given in Section 3.1 for the classical BGK equation, we obtain

[TABLE]

Note that, using the definition (4.11) of the equilibrium distribution for the velocity $v$, we compute

[TABLE]

Thus, the kinetic equation (4.7) solves the advection-diffusion equation

[TABLE]

with

[TABLE]

It is easy to verify that (4.15) is exactly (3.15) if the closure (4.14) is considered. Observe that, compared to (3.6), the diffusion coefficient (4.15) contains two additional terms which depend on the function $p(\rho)$. Therefore, it is possible, for a given distribution $M_{f}$, to find a suitable $p(\rho)$ such that $\mu(\rho)>0$ also in the congested regime. Recall that $\mu(\rho)$ given in (3.6) was unconditionally negative in the congested phase of traffic for the classical BGK model (3.2). In particular, it is possible to find $p(\rho)$ in order to guarantee that Property 1 holds for the equation (4.7) and thus have a small regime where instabilities can occur but they are bounded.

Setting $Q_{\text{eq}}(\rho)=\rho U_{\text{eq}}(\rho)$ the second two terms of the diffusion coefficient (4.15) can be written in terms of the equilibrium speed function as

[TABLE]

Therefore, $\mu(\rho)=\mu_{\text{BGK}}(\rho)+C(\rho)$ where $C(\rho)=-\rho^{2}p^{\prime}(\rho)U_{\text{eq}}^{\prime}(\rho)\geq 0$ since $p$ and $U_{\text{eq}}$ are an increasing and a non-increasing function of the density, respectively. This means that, for enough large $C(\rho)$, the additional term helps to make the diffusion coefficient (3.6) positive. In Figure 4.1 we numerically show this result for the case of the homogeneous kinetic model reviewed in Section 2 with $2$ and $3$ microscopic speeds and two different pressure functions, i.e. $p(\rho)=\frac{3}{2}\rho^{2}$ and $p(\rho)=\rho^{3}$. Observe that, in particular, there is a small density regime where the diffusion is negative and thus equation (4.7) unstable in the viscous limit. Compare Figure 4.1 with Figure 3.6. The given formulations of the pressure functions provide exactly the interaction term (4.2) in the microscopic follow-the-leader model for suitable choices of the parameters $C_{\gamma}$ and $\gamma$.

5 Conclusions

In this paper we studied the stability of the classical BGK formulation of a kinetic model for traffic flow. Via Chapman-Enskog expansion we have shown that the BGK model leads to an advection-diffusion equation with a negative diffusion coefficient for a density greater than the critical one. This leads to unbounded growth of waves.

The issue has been addressed by deriving a modified BGK-type equation as mesoscopic description of the classical microscopic follow-the-leader model. The new model has a second-order moment system of ARZ–type and, as direct consequence, it result to be conditionally well-posed. The instabilities occurring in the regime of high densities have now a bounded growth and can be regarded as model for stop and go waves. Further, with the new kinetic model, we have also introduced the mesoscopic description between the microscopic follow-the-leader model and the macroscopic ARZ model.

Acknowledgments

This work has been supported by HE5386/13-15.

M. Herty and G. Visconti would like to thank the German Research Foundation (DFG) for the kind support within the Cluster of Excellence “Internet of Production” (IoP ID390621612).

G. Puppo and G. Visconti are members of the “National Group for Scientific Computation (GNCS-INDAM)”.