Kinetic models for traffic flow resulting in a reduced space of microscopic velocities

The purpose of this paper is to study the properties of kinetic models for traffic flow described by a Boltzmann-type approach and based on a continuous space of microscopic velocities. In our models, the particular structure of the collision kernel allows one to find the analytical expression of a class of steady-state distributions, which are characterized by being supported on a quantized space of microscopic speeds. The number of these velocities is determined by a physical parameter describing the typical acceleration of a vehicle and the uniqueness of this class of solutions is supported by numerical investigations. This shows that it is possible to have the full richness of a kinetic approach with the simplicity of a space of microscopic velocities characterized by a small number of modes. Moreover, the explicit expression of the asymptotic distribution paves the way to deriving new macroscopic equations using the closure provided by the kinetic model.


(Communicated by Lorenzo Pareschi)
Abstract. The purpose of this paper is to study the properties of kinetic models for traffic flow described by a Boltzmann-type approach and based on a continuous space of microscopic velocities. In our models, the particular structure of the collision kernel allows one to find the analytical expression of a class of steady-state distributions, which are characterized by being supported on a quantized space of microscopic speeds. The number of these velocities is determined by a physical parameter describing the typical acceleration of a vehicle and the uniqueness of this class of solutions is supported by numerical investigations. This shows that it is possible to have the full richness of a kinetic approach with the simplicity of a space of microscopic velocities characterized by a small number of modes. Moreover, the explicit expression of the asymptotic distribution paves the way to deriving new macroscopic equations using the closure provided by the kinetic model.

1.
Introduction. The purpose of kinetic theory for traffic flow is to provide an aggregate representation of the distribution of vehicles on the road, thanks to a detailed characterization of the microscopic interactions among the vehicles, which play an important role in the macroscopic trend of the flow. The goal is to obtain Existence of quantized equilibria: All equilibria are quantized:  This framework permits to take the stochasticity of the drivers' behavior into account, thanks to the probability distribution which assigns a weight to the possible driver's decisions, while maintaining the general kinetic setting, based on the deterministic evolution of the distribution function. We propose two models that follow the framework just described: one is based on quantized velocity jumps, i.e., if acceleration occurs the new speed is obtained by increasing the pre-interaction velocity of a fixed quantity ∆v (δ model); the other one is based on a continuous uniform distribution defined on a bounded interval parametrized by ∆v (χ model), see also [14,15].
In this paper, ∆v is a finite parameter which models the physical velocity jump performed by vehicles when they increase their speeds as a result of an interaction. Clearly, this parameter may depend on the mechanical characteristics of vehicles, see [24], but in this paper we will assume that ∆v is fixed. In the section on macroscopic properties, §5, we show that ∆v is related to the maximum acceleration and we discuss how this parameter can be chosen through experimental data used in [16].
This paper also shows that the fundamental diagrams (or closure laws) obtained with the perhaps more natural, but more complex and more computationally demanding, χ model are very close to those provided by the simpler δ model. We thus investigate the equilibria of the δ model, both from an analytical and a numerical point of view. Analytically we find that velocity distributions formed by a linear combination of Dirac δ's may be equilibria only if the δ's are centered at velocities spaced by multiples of ∆v. Next, we compute equilibria using a numerical scheme capable of converging also to possible absolutely continuous equilibria. Here we find only the quantized equilibria described above, independently of the discretization parameters. This fact suggests that the class of discrete-velocity equilibria is the only one that the continuous-velocity Boltzmann-type δ model possesses. This situation is summarized graphically in Figure 1.
The paper is organized as follows. In Section 2 we briefly recall the Boltzmanntype kinetic equation and we specialize it by giving two sets of interaction rules.
The resulting δ and χ models are discussed in depth in Section 3 and in Section 4, respectively. In particular, in Section 3 we prove the existence of a class of quantized steady-state distributions. Since we are unable to prove their uniqueness, we discretize the model by approximating the kinetic distribution with a piecewise constant function and we then show by numerical evidence that the class of stationary solutions of the δ model is only the one we have already studied analytically. In Section 4, we show the somewhat surprising result that the equilibrium distributions of the χ model yield a macroscopic flow that is extremely well approximated by the discrete-velocity-based closure law resulting from the δ model. This is illustrated in the final Section 5, where we show that the fundamental diagrams, that is the flux-density relationships, obtained from the two models tend to coincide under grid refinement. Next we compare these diagrams with experimental data, finding that our models reproduce well experimental fundamental diagrams and thus they capture the characteristics of macroscopic traffic flow. Further, we compute the macroscopic acceleration induced by the model, proving in particular its link with ∆v. We end the paper with a section summarizing the main results of this work, and proposing possible applications and further developments. Finally, the details of the matrix elements resulting from the discretization of the χ model are written in an Appendix.
2. The general form of the kinetic model. In this section we briefly recall the general structure of a Boltzmann-type kinetic traffic model, which we will then specialize by prescribing a set of binary interaction rules in order to derive two models which differ only in the modeling of the acceleration interaction. Both models are defined on a continuous velocity space and they are characterized by a parameter ∆v related to the typical acceleration of a vehicle.
We will focus on the space homogeneous case, because we want to investigate the structure of the collision term and of the resulting equilibrium distributions. In particular, we show that the simplified model (δ model) permits to describe the complexity of the equilibrium solutions with a very small number of discrete velocities.
Let f = f (t, v) : R + × V → R + be the kinetic distribution function, where V = [0, V max ] is the domain of the microscopic speeds and V max is the maximum speed, which may depend on the mechanical characteristics of the vehicles, on imposed speed limits, environmental conditions (such as the quality of the road, the weather conditions, etc). The statistical distribution f is such that f (t, v)dv gives the number of vehicles with velocity in [v, v + dv] at time t.
As usual, macroscopic quantities are obtained as moments of the distribution function f with respect to the velocity v: where ρ is the density, i.e. the number of vehicles per unit length (typically, kilometers), u is the macroscopic speed and ρu is the flux of vehicles. Note that ρ can also be interpreted as the reciprocal of the average distance between cars, see [2].
In the homogeneous case, the Boltzmann-type equation can be written as where Q[f, f ](t, v) is the collisional operator which describes the relaxation to equilibrium due to the microscopic binary interactions among vehicles. For mass conservation to hold, the collision term must satisfy In fact, this ensures that, in the space homogeneous case, the density remains constant in time.
The collisional operator is usually split into a gain term G[f, f ] and a loss term L[f, f ], that model statistically the interactions which lead to gain or to loose the test speed v. Denoting with A(v * → v|v * ) the probability that the velocity v ∈ V results from a microscopic interaction between candidate vehicles with velocity v * and field vehicles with speed v * , the model writes as an integro-differential equation in which η(v * , v * ) is the interaction rate possibly depending on the relative speed of the interacting vehicles, e.g. η(v * , v * ) = |v * − v * | as in [5,14]. Although such a choice would make the model richer, in [25] we found that a constant interaction rate is already sufficient to account for many aspects of the complexity of traffic. Another possibility is to consider η as dependent on the local congestion of the road, that is η = η(ρ). However this is not relevant in the homogeneous case, where ρ is constant, for then η would affect only the relaxation time towards equilibrium. Thus in this paper we will set η =constant.
Notation. In the whole paper, in order to shorten formulas, we adopt the following traditional shorthand f (t, v * ) = f * , f (t, v * ) = f * , etc. Note in particular that in the space homogeneous case f * and f * are not different distribution functions, but the evaluation of the same f (t, v) at two different points v * and v * .
We will suppose that A depends also on the macroscopic density ρ in order to account for the influence of the macroscopic traffic conditions (local road congestion) on the microscopic interactions among vehicles, see [23,14,11,25]. Thus, we suppose that A fulfills Assumption 1.
where ρ max is the maximum density of vehicles, for instance the maximum number of vehicles per unit length in bumper-to-bumper conditions. Remark 1. Any transition probability density A that satisfies Assumption 1 guarantees mass conservation since 2.1. Choice of the probability density A. The probability density A assigns a post-interaction speed in a non-deterministic way, consistently with the intrinsic stochasticity of the drivers' behavior. The construction of A is at the core of a kinetic model. Here, it is obtained with a very small set of rules.
• If v * ≤ v * , i.e. the candidate vehicle is slower than the field vehicle, the post-interaction rules are: Do nothing: the candidate vehicle keeps its pre-interaction speed with probability 1 − P 1 , thus v = v * ; Accelerate: the candidate vehicle accelerates to a velocity v > v * with probability P 1 . • If v * > v * , i.e. the candidate vehicle is faster than the field vehicle, the post-interaction rules are: Accelerate: in order to overtake the leading vehicle, the candidate vehicle accelerates to a velocity v > v * with probability P 2 ; Brake: the candidate vehicle decelerates to the velocity v = v * with probability 1 − P 2 , thus following the leading vehicle. From the previous rules, we observe that the probability density A has a term which will be proportional to a Dirac delta function at v = v * , due to the interaction which preserves the pre-interaction microscopic speed (the "Do nothing" alternative). Note that this is a "false gain" for the distribution f , because the number of vehicles with speed v is not altered by this interaction.
In the following, we assign the speed after braking as proposed in [22] and used also in [7,8] in the context of a discrete velocity model. Namely, we suppose that if a vehicle brakes, interacting with a slower vehicle, it slows down to the speed v * of the leading vehicle. Thus, after the interaction it gets the speed v = v * without overtaking the leading field vehicle. Instead, for the post-interaction speed due to acceleration we propose two different models.
Quantized acceleration (δ model): the output velocity v is obtained by accelerating instantaneously from v * to the velocity min {v * + ∆v, V max }. Considering all possible outcomes, the resulting probability distribution, in this case, is Uniformly distributed acceleration (χ model): the new velocity v is uniformly distributed between v * and min{v * + ∆v, V max }. On the whole, the resulting probability distribution becomes Note that the acceleration of a vehicle in (3) is similar to the one assumed in [7,8], which however were based on a discrete velocity space. In [7,8] the acceleration parameter ∆v is chosen as the distance between two adjacent discrete velocities, thus ∆v depends on the number of elements in the speed lattice. In this work, ∆v is a physical parameter that represents the ability of a vehicle to change its preinteraction speed v * . With this choice, ∆v does not depend on the discretization of the velocity space and the maximum acceleration is bounded, as in [16]. In contrast, deceleration can be larger than ∆v, and this fact reflects the hypothesis that drivers tend to brake immediately if the flow becomes more congested, while they react more slowly when they can accelerate (see the concept of traffic hysteresis in [29] and references therein).
The acceleration performed in the χ model has some points of contact with the microscopic rules prescribed in [12,14]. In [12], however, the post-interaction speed is selected through a random process in the interval [v * , V max ]. Here instead, the post-interaction speed is deterministic. In [14] instead, the velocity after acceleration is uniformly distributed over a range of speeds between v * and v * +α(V max −v * ), where α is supposed to depend on the local density; in a similar way, the output velocity from a braking interaction is assumed to be uniformly distributed in [βv * , v * ], with β ∈ [0, 1].
In the following, the probabilities P 1 and P 2 are taken as P 1 = P 2 =: P and P will be a function of the local density only, as assumed for instance in [23] where P = 1 − ρ/ρ max . More generally, from a modeling point of view, P should be a decreasing function of ρ, see also [10] or [26]. For instance in [25] we have considered where γ ∈ (0, 1) can be chosen to better fit experimental data. In [7] the choice is a parameter describing environmental conditions, for instance road or weather conditions. In more sophisticated models, one could also choose P as a function of the relative speed of interacting vehicles, but we will not explore this possibility in the present work. The simplified choice P 1 = P 2 and the interaction rules described at the beginning of this section guarantee the continuity of the transition probability (3) and (4) along v * = v * . (3) and (4) for A include terms of the form δ v * (v), which actually describe false gains mentioned above, because the velocity of the candidate vehicle does not change. They are automatically compensated by false losses, as it can be seen by rewriting the classical kinetic loss term of equation (2) in the form

Remark 2. Both choices
3. The δ velocity model. Now, we focus on the steady states of model (3). We start with the existence of a particular set of equilibrium solutions of the continuous model, which are computed analytically. Next, we consider a finite volume discretization of the model, and we show that the discrete equilibria have precisely the structure found before analytically, thus suggesting that the particular set of equilibria found analytically are the only equilibria of the system.
It can be proven [9] that the Cauchy problem associated to (2) is well posed provided the probability density A is Lipschitz continuous with respect to v * and v * in a suitable Wasserstein metric. This is indeed the case of the A defined in (3), with P 1 = P 2 .
Using the expression (3) for A, we rewrite the gain term in (2) as Proof. Without loss of generality suppose that {v j } N j=1 is an ordered set of velocities is a weak stationary solution of the δ model, it satisfies the following steady weak form of equation (2): where φ ∈ C c (V ) is a test function, with C c (V ) the space of continuous functions having compact support contained in V , and the probability density A is given in (3). Substituting the expression of f ∞ in the above equation we obtain The proof will be organized as follows: in order to determine an equation for the f ∞ j 's, we consider a particular family of test functions φ j defined as piecewise linear functions such that φ j (v j ) = 1 and φ j (v i ) = 0, ∀ i = j; in this way, first we find an equation for f ∞ 1 , then we show that f ∞ 2 = 0 if v 2 = v 1 +∆v and finally by induction we prove that if f ∞ j = 0, then v j = v 1 + j∆v, for some j ∈ {1, . . . , N }. Let j = 1, equation (6) Due to the particular construction of φ 1 and using N j=1 f ∞ j = ρ, the above expression reduces to If P > 1/2, only f ∞ 1 = 0 is acceptable, because the other root is negative. If instead P < 1/2, both roots can be accepted, but only f ∞ 1 = ρ 1−2P 1−P > 0 is stable. This argument will be used for selecting a single root throughout the proof. Now, let j = 2 and φ = φ 2 . Equation (6) writes as Note that φ 2 (min{v h +∆v, V max }) = 1 if h = 1 and v 1 +∆v = v 2 . While φ 2 (min{v h + ∆v, V max }) = 0 otherwise, due to the particular choice of φ 2 which is centered in v 2 and can be taken with support smaller than 2 |v 2 − v 1 − ∆v| around this point. Then, if v 2 = v 1 + ∆v, the constant coefficient of (7) is zero for all P ∈ [0, 1] and the solutions of the equation Instead, if v 2 = v 1 + ∆v, the third term in equation (7) is P ρf ∞ 1 ≥ 0, for all P ∈ [0, 1]. More precisely, if P ≥ 1/2 then f ∞ 1 = 0, thus the constant coefficient vanishes and again one concludes that (7) is positive and the leading coefficient is negative, the equation has two real roots with opposite signs. Therefore and this is the only case in which f ∞ 2 can be non-zero. We now proceed by induction. Suppose that v k − v 1 is an integer multiple of ∆v and Taking the test function φ = φ j , the equation for f ∞ j writes as

∆v and the above equation becomes
which has two real roots. If P ≥ 1/2, using the inductive hypothesis f ∞ k = 0 for all k ≤ j − 1. Thus f ∞ j = ρ 1−2P 1−P (which is again not acceptable since it is negative) or f ∞ j = 0, confirming the induction. If P < 1/2, we have two real roots with opposite signs, so one proves that f ∞ j can be chosen strictly positive. If instead v j = v 1 + j∆v then the constant term of equation (8) is zero and the two roots are f ∞ which is negative for all values of P ∈ [0, 1] using the Lemma 3.2 below.
In the previous proof we use the following technical fact.
Proof. We proceed by induction on k. Let k = 2 To this end, the explicit formulation of the gain term is now useful. Notice that the Dirac delta function at v = min {v * + ∆v, V max } can be split as because the velocity jump of size ∆v, leading to the output velocity , the post-interaction velocity will be v = V max . Thus the gain term of the δ model can be written as where H α (x) denotes the Heaviside step function with jump located in α. The last term in the expression of G means that, as a result of the microscopic interactions, the mass P ρ

Vmax
Vmax−∆v f * dv * is allocated entirely to the speed V max . Note that, in space-nonhomogeneous models, f * and f * may refer to distributions evaluated at different locations in space, see for instance [14] and [15]. For this reason we keep the integrals over field and candidate particles separate.
Suppose for simplicity that the acceleration parameter ∆v satisfies ∆v = V max /T with T ∈ N. We consider a discretization of the velocity space defining the velocity Figure 2. Structure of the probability matrices of the δ model, with ∆v = V max /2. The shaded areas correspond to the nonzero elements of the matrices. For the meaning of the different hatchings, please see Table 1 cells Note that all cells have amplitude δv = V max /(N − 1) except I 1 and I N which have amplitude δv/2.
We consider a piecewise constant approximation of the kinetic distribution so that where f j represents the number of vehicles traveling with velocity v ∈ I j . By integrating the kinetic equation (1) over the cells I j and using f N (t, v) in place of f (t, v) we obtain the following system of ordinary differential equations whose initial conditions and ρ is the initial density, which remains constant during the time evolution in the spatially homogeneous case. Pattern Entries of the matrices Proportional to P 1 − P P 1 Table 1. Table describing the patterns in the matrices of Figure 2, 3 and 9.
As it can be checked using (12), these matrices are stochastic with respect to This property comes from Assumption 1, and it guarantees mass conservation.
Recall that the elements of the matrix A j δ hk are the probabilities that the candidate vehicle with velocity in I h , interacting with a field vehicle with velocity in I k , acquires a velocity in I j . The fact that these matrices are sparse means that a velocity in I j can be acquired only for special values of the velocity of candidate and field vehicles. In particular, the j-th row of the matrix A j δ contains the probability of what we called "false gains" in Remark 2, that is the probability that the candidate vehicle does not change its speed. The non zero elements of the j-th column are the probabilities that a candidate vehicle acquires a speed in I j by braking down to the speed of the leading vehicle. The non zero rows, located at h = j − r + δ r , r + , h − 1 and h + 1, contain the probabilities that the candidate vehicle accelerates by ∆v, acquiring therefore a velocity in v * + ∆v ∈ I j starting from a velocity v * in I h−1 , I h or I h+1 . The band between the rows j − r + δ r , r + + 2 and j − 1 is filled with zeros, because in the δ model the acceleration is quantized. As we will see in Section 4.1, this band will be filled by non zero elements in the χ model, where the acceleration is distributed uniformly between [0, ∆v].
In Table 2 we summarize the numerical parameters introduced in order to discretize the continuous-velocity model.
3.1.1. δv as integer sub-multiple of ∆v. For the special choice of the velocity grid which ensures that r = ∆v/δv ∈ N, the formulae (12) simplify and the resulting interaction matrices are given in Figure 3. Notice that the rows j − r ± 1 are filled with zeros. In fact, since δv is an integer sub-multiple of ∆v, a velocity in I j can be obtained as result of an acceleration only if the pre-interaction speed is a velocity in I j−r .
The structure of the matrices A j δ determines the equilibrium of the discrete model (13). In Figure 4 we show the function f ∞ Figure 3. Structure of the probability matrices of the δ model with δv integer sub-multiple of ∆v. The shaded areas correspond to the non-zero elements of the matrices. For the meaning of the different hatchings, please see Table 1.  We mark with red circles on the x-axes the center of the T + 1 cells obtained with r = 1. integrating numerically the system of equations up to steady state, for a few typical cases. In all numerical tests we take P = 1 − ρ, with V max = ρ max = 1. As initial macroscopic densities, we choose ρ = 0.3, 0.6 (plots to the left and right of This means that, as δv → 0, only a finite number T + 1 of velocities carry a non-zero mass of vehicles at equilibrium. More precisely, the discrete asymptotic function f ∞ N (v) is different from zero only in the T + 1 cells I 1 , I r , I 2r , . . . , I N . Therefore, as time goes to infinity the number of nonzero values of the f j 's appearing in (10) is univocally determined by the acceleration term ∆v = V max /T . The previous considerations can be supported by the numerical results in Figure  5, in which we show the evolution towards equilibrium of the f j 's, j = 1, . . . , N . In this figure, ∆v = V max /3, and the different plots are obtained starting from a uniform initial distribution, namely f j (t = 0) = ρ/N, ∀j = 1, . . . , N , with r = 1 (green), r = 2 (blue) and r = 3 (red), which correspond to N = 4, 7 and 10 velocity cells respectively, and ρ = 0.6. It is clear that under grid refinement the number of nonzero steady values does not change. In fact, note that a different dynamics towards equilibrium is observed, for different values of the number N of cells, but as equilibrium is approached, the values of the f j 's go to zero except for the velocities corresponding to integer multiples of ∆v. Moreover, the non zero values of the steady-state distribution f ∞ do not depend on the discretization parameter δv. This fact can be also deduced by looking at equations (12) of the discrete collision operator. These expressions are not functions of δv. Thus all exact values of the equilibria can be obtained using a coarse grid. Theorem 3.4 below confirms the structure of the equilibria that we have just observed in the numerical results and it states that the steady-state solution of the δ model, prescribed by Theorem 3.1, can be reconstructed numerically on the grid with δv = ∆v. To this end, we first recall a result from [6], where the existence and well posedness of the solution of such systems is proved and then we show that the discrete equilibria of additional equations resulting from the choice r > 1 give actually no contribution.
where the matrices A j are stochastic matrices with respect to the index j, i.e. j A j hk = 1 for all h, k. Then the system admits a unique non-negative global solution f ≥ 0 satisfying the a priori estimate The following result, together with Remark 3, shows that all equilibria of the discrete model are of the quantized form described by Theorem 3.1. Thus it establishes the correspondences symbolized by the right vertical and the middle horizontal arrows in Figure 1. T +1 k=1 (f 1 ) k = ρ. Proof. We already know from Theorem 3.3 that the solution of (13) exists, is nonnegative and is uniquely determined by the initial condition.
To prove the statement, we compute explicitly the equilibrium solutions of (13), using the explicit expression of the collision kernel given in (12a), (12b), (12c) and (12d), with r ∈ N and N = rT + 1. Since here we are interested in the solutions of the homogeneous problem, we will take identical distributions for the candidate and the field vehicles, i.e. f j = f j .
For j = 1, using the expression (12a) and the fact that This is a quadratic equation for f 1 , which has the two roots f 1 = 0 and f 1 = ρ(1 − 2P )/(1 − P ). It is easy to see that one solution is stable, and the other one unstable, depending on the value of P . Here we are interested only in the stable root, so we find Thus, no vehicle is in the lowest speed cell I 1 if P ≥ 1 2 , which, for the simple case P = 1 − ρ means that all cars are moving if ρ ≤ 1 2 . The case j = 1 we just computed is typical. Also for larger values of j, we find a quadratic equation for the unknown f j , which involves only previously computed values of f k , k < j. Thus, we can easily compute iteratively all components of f r .
For 2 ≤ j ≤ r, the equilibrium equation is obtained by using again the expressions (12a) with r ∈ N. Then Start from j = 2. Clearly, for P ≥ 1 2 , substituting equation (14), we again have (f r ) 2 = 0. For P < 1 2 , the equation for f 2 , with f 1 given by (14), becomes Comparing with the equation for f 1 , we see that now the stable root is f 2 = 0. Thus, at equilibrium, we have (f r ) 2 = 0, for all values of P . Analogously, it is easy to see that (f r ) j = 0, ∀j = 3, . . . , r.
For r + 1 ≤ j < N , in place of (12a), we use (12b) and (12c) in order to obtain the equilibrium equation. Then The equation has a positive discriminant thus it always admits two real roots. To fix ideas, let us consider r + 1 ≤ j ≤ 2r. If j = r + 1, since (f r ) k = 0, ∀k = 2, . . . , r, equation (15) becomes Thus, we find (f r ) r+1 = 0 for P ≥ 1/2, because the equation for f r+1 becomes identical to (3.1.1). If instead P < 1/2, substituting the expression for f 1 , the equation for f r+1 becomes This equation has a negative and a positive real root, which is stable. Thus otherwise. Now, let r + 1 < j ≤ 2r. Since f j−r = 0, the constant term of equation (15) is zero. Then, as seen for j = 2, . . . , r, it is easy to prove that, for each j = r + 2, . . . , 2r and for all values of P , one solution is negative and thus (f r ) j = 0.
Clearly, this procedure can be repeated and we find that otherwise if j = lr + 1, l = 0, . . . , T − 1, while (f r ) j = 0 otherwise. From these considerations, the thesis easily follows. Finally, just note that for the last value, f N , we can use mass conservation Notice that the result above proves that equilibria are determined by the initial density ρ, which is constant in time in the spatially homogeneous case, but they do not depend on the number of cells N used to approximate the kinetic distribution. In fact, the stable equilibria just computed correspond exactly to the values f ∞ Remark 3 (The case of a generic δv). In order to further investigate the existence of stable equilibria, we can seek more general ones numerically. In fact, the finite volume discretization (10) is capable of converging to absolutely continuous equilibria, but so far we proved that it converges to sums of Dirac masses if δv is an integer sub-multiple of ∆v. Here we apply our discretization scheme with noninteger ratios r = ∆v/δv and show that the resulting equilibria converge, in the sense of distributions, to the equilibria already described in Theorems 3.1 and 3.4.
When r = ∆v/δv is not integer, Theorem 3.4 cannot be applied, but numerical integration of equation (13) shows that, for large time, the solution approaches equilibria that have masses concentrated at points spaced by ∆v. More precisely, in this more general case, only a small finite number of f ∞ j 's is nonzero and the cumulative distribution function induced by the discrete f ∞ N (v), cf. (10), approximates the cumulative distribution of a sum of Dirac masses centered at multiples of ∆v. That is, F N (v) converges to a piecewise constant function with jump discontinuities at multiples of ∆v. See Figure 6. Furthermore, Figure 6 shows that the jump discontinuities of the cumulative distribution in the limit δv → 0 are located exactly in the points computed analytically by Theorem 3.4 in the case r ∈ N. In particular, the figure shows the cumulative distribution at equilibrium computed by solving numerically (13) for several values of the discretization parameter δv with non-integer ratio ∆v/δv. In both panels the density is ρ = 0.6, while ∆v = 1/3 in the left plot and ∆v = 1/5 in the right one. The continuous red curve represents the cumulative distribution of the stationary solution for δv = ∆v, computed by Theorem 3.4. The black dotted curve is computed with N = 15 velocity cells, the green dot-dashed line with N = 30 and finally the blue dashed one with N = 60. We observe that as δv → 0 the cumulative distribution tends to that of the analytic solution of Theorem 3.1, where only the velocities centered in multiples of ∆v give rise to a jump discontinuity.   N − 3)). The thick lines highlight the components f j and the blue ones are for those that appear in stable equilibria, i.e. with j = kr + 1 for k = 0, . . . , T .
We conclude the section with a few remarks on the structure of the equilibria of (13).
Remark 4 (Reduced velocity space). The δ model is characterized by a small number of non-trivial values for the microscopic velocities, which allows one to compute analytically the equilibrium of the system. This fact provides an analytic closure of the macroscopic equations. This fact can also be exploited from a numerical point of view, justifying the use of a small number of discretization points in simulations aimed at capturing equilibrium effects. In this sense, the model seems to suggest that kinetic corrections can be accounted for with a computational cost which is not much higher than the one needed for a macroscopic model.
Remark 5 (Unstable equilibria). Theorem 3.4 gives the uniqueness of the stable equilibrium of the model with ∆v/δv = r ∈ N. Unstable ones may occur if the initial condition is such that f 1 (0) = 0. In fact, the interaction rules related to the case v * > v * do not generate a post-interaction velocity v which is lower than v * . Thus if f 1 (0) = 0, i.e. if there are no vehicles with velocity v 1 = 0 at the initial time, there will not be interactions leading to an increase of f 1 . This consideration can be generalized: if f j (0) = 0 for j = 1, . . . ,j < r, then the computed equilibria will be f j = (f r ) j−j , where f r is the vector containing the stable equilibria. In this sense, the equilibrium solution of the δ model does not only depend on ρ but also on the initial condition f (0, v). These solutions however are unstable: a small perturbation on f 1 (t = 0) is enough to trigger the evolution towards the stable equilibrium, which depends only on ρ.
This is illustrated in Figure 7: in the left panel we show the evolution towards equilibrium when f j (t = 0) = 0 for all classes (this is the stable equilibrium), while in the middle we show the case when f 1 = f 2 = f 3 = 0. In the rightmost panel we show a perturbation of the previous case, where f 1 takes a very small but nonzero value. It is clear that the evolution goes at first towards the unstable equilibrium of the middle panel, but then, in the long run, the stable equilibrium of Theorem 3.4 emerges.
Remark 6 (Convergence rate to equilibrium). In Figure 8 we study the rate of convergence towards the stable steady-states. For the set of densities ρ ∈ {0.2, 0.3, 0.4, 0.6, 0.7, 0.8} we integrate numerically the system (13) for large times   Figure 8. Speed of convergence towards the stable equilibria of the δ model. The initial condition is a small random perturbation of the steady-states.
starting from a small random perturbation of the stable equilibrium. In both panels, we use a linear scale for the x-axis (time) and a logarithmic scale for the y-axis (error). The error is computed at each step as The figure suggests that the rate of convergence towards the stable equilibrium depends on the density. In fact, in the left panel we consider different values of the ratio ∆v/δv = r, in particular r = 1 (solid lines) and r = 2 (dashed lines), and we note that the slopes of the curves corresponding to the same value of ρ are the same. A similar behavior is observed in the right plot in which two different ∆v are taken, ∆v = 1/5 (solid lines) and ∆v = 1/3 (dashed lines). We can conjecture that for large enough times the distance from equilibrium behaves as e(t) Ce −M (ρ)t f (t = 0) − f ∞ 2 , where C = C(r, ∆v), M (ρ) > 0. 4. The χ velocity model. The structure of the steady-state distribution of the δ model clearly depends on the particular choice of the acceleration interaction made in (3), in which a vehicle accelerates by jumping from its pre-interaction velocity v * to the new velocity v * + ∆v. Thus it could seem quite natural that only velocities 0, ∆v, 2∆v, . . . , V max give a non zero contribution at equilibrium.
Here we study the χ model, already introduced in Section 2.1 (see equation (4)), in which vehicles can assume a post-interaction velocity uniformly distributed over a range of speeds when the acceleration interaction occurs. We will show that, although this model is more refined than the δ model, at equilibrium the essential information is already caught by the simpler δ model.
Using the formulation (4) for the transition probability density A, we rewrite the gain term in (2) as Notice that the χ function can be split as Figure 9. Structure of the probability matrices of the χ model with δv integer sub-multiple of ∆v.
hence substituting in the above equation and evaluating explicitly the integrals, we find Observe that (16) differs from (9) only in the terms proportional to P .

Discretization of the model.
To compute the steady-state solution of the χ model, we need to integrate the equations numerically. We use the same discretization of the velocity space V introduced to discretize the δ model and therefore the kinetic distribution is approximated as in (10). Integrating the kinetic equation (1) over each cell, we find the system of ODEs (11), but now the gain term is given by (16). Although in this case the integrals are laborious, they can be computed recalling that ∆v = V max /T with T ∈ N and assuming that δv is an integer sub-multiple of ∆v. Thus we will take N − 1 ≡ 0 (mod T ) and r = N −1 T . The details are reported in Appendix A. Here we just point out that, as in the case of the δ model, the ODE system can be conveniently rewritten in vector form as where f = [f 1 , . . . , f N ] T ∈ R N is the vector of the unknown functions, e j ∈ R N denotes the vector with a 1 in the j-th coordinate and 0's elsewhere, 1 T N = [1, . . . , 1] ∈ R N and A j χ is the j-th interaction matrix such that A j χ hk contains the probabilities that a candidate vehicle with velocity in I h interacting with a field vehicle with velocity in I k acquires a velocity in I j .
Since the matrices appearing in (17) are again stochastic, we can apply Theorem 3.3 to guarantee the well posedness of the associated Cauchy problem.
In Figure 9 we show the structure of the χ matrices: they are less sparse than the δ matrices because of the uniformly distributed acceleration in [v * , v * + ∆v], where v * is the pre-interaction speed. In fact the matrix A j χ contains non-zero elements also in the rows from the (j − r + 1)-th to the (j − 1)-th, see the shaded areas in Figure 9, which represent non-zero probabilities of accelerating to a speed in I j . Since ∆v = rδv, exactly r rows fill up. Instead, the area drawn using hatching contains the same probabilities already shown in Figure 3 for the case of the δ model, with r ∈ N. Note that the elements of the matrices depend on δv (see equations (A.1)). Thus, in contrast to the δ model, steady solutions of the ODE system (17) depend on the number of velocity cells N chosen to approximate the kinetic distribution (see (10)). In other words, although δv is an integer sub-multiple of ∆v, this model does not converge, as time goes to infinity, to the asymptotic distribution on the coarse grid as in the δ model. However, we do recover the asymptotic distribution on the coarse grid in the limit δv → 0.
Notice from equations (A.1) that all the elements in the rows j −r, . . . , j −1 of the matrices A j χ , j = 1, . . . , N , tend to 0 as 1/r when the grid is refined. In particular, for j = 1, . . . , r, A j χ → A j δ . This consideration is not true for the matrices A j χ , for j = r + 1, . . . , N .

4.2.
Expected speed of the δ and the χ model. Despite their differences the χ and the δ model are deeply related. This can be seen by computing the expected output speed in each model resulting from a fixed pre-interaction speed. We define the expected value v of the post-interaction velocity as For brevity we indicate with A δ (v) and A χ (v) the probability densities given in (3) and (4) respectively. Again we assume that P 1 = P 2 = P as function of the density ρ. For the δ model we obtain In contrast, if we consider the χ model we have By comparing the last lines of (19) and (20), it is clear that Remark 7 (Uniformly distributed deceleration). The computation of the expected speed shows a link between the δ model and the χ model. In fact, under the constraint ∆v δ = 1 2 ∆v χ , the two models provide the same expected speed. Similarly, if we consider a braking scenario in which the candidate vehicle brakes to a speed uniformly distributed in an interval centered on v * , then the expected speed results to be again (1 − P )v * .
Remark 8. Let us compare the A j χ matrices (Figure 9) for r < j ≤ N − r with a given ∆v and the corresponding A j δ matrices (Figure 3) with ∆v 2 . For the case r ∈ N, the isolated nonzero row of A j δ is at j − r 2 , which corresponds to the middle of the green shaded area in A j χ . Moreover, for any fixed ∆v, it can be proved that, for δv → 0, the sum of the quantities located in the shaded area of A j χ is equal to the total contribution provided by the j − r 2 -th row of the δ matrices obtained with the acceleration parameter ∆v 2 . In other words, as δv → 0, the total effect of the j − r 2 -th row of A j δ with ∆v 2 is spread over r + 1 rows in the matrices A j χ with ∆v, which are the rows shaded in green in Figure 9.

Macroscopic acceleration.
In order to explain the relation of ∆v with the acceleration of the vehicles in the model, we compute the rate of change of the macroscopic velocity: where v is a function of v * and v * as defined in (18). In the case of the δ model, v δ is given by (19) and thus Given an initial distribution f (t = 0, v), the equation above yields the evolution of the macroscopic acceleration in time. It is easy to study analytically this quantity at the initial time. In particular, we compute the initial acceleration in the case in which all vehicles are still but the density is below the value for which P = 1/2 (that is ρ = 1/2 when taking P = 1 − ρ).   The above equation shows that the acceleration of the vehicles in the δ model depends linearly on ∆v. Analogously, for the χ model, using (20), one obtains which reinforces the remark made in (21) about the similarities of the χ and the δ model when ∆v χ = 2∆v δ .
Remark 9 (Acceleration). Recall that (ηρ) −1 is a time. Thus ηρ∆v is the builtin acceleration of the model, which, not surprisingly, is linked to ∆v. We can use dimensional arguments to estimate the order of magnitude of ∆v. According to [16], the maximum acceleration of cars is approximately 2.5 m/s 2 . The maximum speed is approximately V max 28m/s, and we expect the maximum acceleration when P = 1. Thus ηρ∆v 2.5 m/s 2 .
The estimates above provide the trend of the macroscopic acceleration starting from rest. For the general case, we now study the evolution of the macroscopic velocity u in time, up to steady state. These data are shown in Fig. 10 and Fig. 11, for various combinations of the model parameters. The results shown are obtained integrating the equations for the δ and the χ model found in (13) and (17), respectively, with r ∈ N, and computing at each time u(t) = 1 ρ V vf (t, v) dv. Figure 10 shows a typical case in which we expect acceleration. The density is ρ = 0.15, well below the value corresponding to P = 1/2 when P = 1 − ρ, and we start with an initial distribution in which f 1 (t = 0) = ρ, while f j (t = 0) = 0, for all j ≥ 2. Thus initially all vehicles are still, and, since the density is low, they will accelerate to reach the maximum speed. The duration of the transient depends on the product η∆v = η/T , for a fixed density, as it is apparent in the right panel of the figure color. As expected, the macroscopic velocity for the χ and the δ model behave very similarly, provided the parameter ∆v is chosen correctly.
Next, in Figure 11, we show the evolution of the macroscopic velocity in two cases when we expect deceleration for the δ model. Namely, we consider ρ = 0.65 in the left panel and ρ = 0.9 in the right panel. The initial distribution is f N (t = 0) = ρ− , f 1 (t = 0) = , and f j (t = 0) = 0, j = 2, . . . , N − 1. The value is introduced to ensure convergence to the stable equilibrium, see Remark 5. Here = 0.01. In other words, we start with a congested traffic, in which initially almost all vehicles are traveling at the fastest speed available. Clearly, this situation is somewhat artificial, but surely we expect the vehicles to brake. Since braking does not depend on ∆v, we expect that the relaxation time towards equilibrium depends mainly on η and only weakly on T . This is clearly seen in both pictures. The macroscopic speed to which the model relaxes on the other hand will depend on ∆v and on ρ, but not on η. Note that when ρ = 0.9, in all cases considered here, the equilibrium speed is nearly zero: in fact the traffic is extremely jammed. For ρ = 0.65 instead, we expect that the traffic will have a residual speed, because we are well below the value P = 1/2, but cars are not "bumper to bumper", and this residual speed does depend on ∆v.

Fundamental diagrams.
As already discussed in the previous section, the nonzero elements of the matrix A j χ can be lumped in the matrix A j δ for N sufficiently large, with the only exception of the elements in the rectangle r × N (see Figure  9) of the matrices for j = N − r + 1, . . . , N . This is shown, for instance, in the evaluation of the expected value of the resulting speeds due to acceleration interactions in (19) and (20), which are comparable, except again for high speed values close to V max (and different from it by at most ∆v). Thus, although the χ model is apparently more refined then the δ model, we expect both models to provide similar macroscopic information, for large N . This is usually analyzed by computing the density and the flux as moments of the asymptotic kinetic distribution f ∞ (v):  Notice that, for the δ model, Theorem 3.4 ensures that only few velocities, obtained with δv = ∆v, are necessary to describe completely the exact asymptotic kinetic distribution. We expect therefore that the macroscopic behavior of the δ model will be apparent even on the coarse velocity grids, i.e. for r = 1. Figure 12 shows the fundamental diagrams provided by the δ model (blue curves) and the χ model (red curve), computed with ∆v δ = 1 2 ∆v χ , for two different values of ∆v δ . In the left panel, r = 1, while r = 20 on the right. The figure shows that the diagram of the χ model is very similar to the diagram of the δ model when N → ∞ and this result is in agreement with the fact that the expected output speed of the two models is mostly the same (i.e., the same in a large range of preinteraction speeds) when choosing the acceleration parameter of the δ model as a half of the acceleration parameter of the χ model. The only difference is provided by the maximum speeds which, as already noted, are slightly different. Note that the similarity of the fundamental diagrams does not mean that the asymptotic equilibrium functions of the χ and of the δ model converge to the same function as N tends to ∞.
Observe that the fundamental diagrams given by the δ model in both plots in each line of Figure 12 use the same information. In fact, following the results of Theorem 3.4, the macroscopic flux is given by  Figure 13. Fundamental diagrams resulting from the δ model with acceleration parameter ∆v δ = 1 4 . The probability P is taken as in (5) with γ = 1 (blue data), γ = 3/4 (green data) and γ = 1/4 (cyan data). The dashed lines are the fluxes in the limit r → ∞.
where v j denotes the center of the cell I j and f r is the vector containing the equilibria of the system (13) with ∆v/δv = r ∈ N. Recalling the definition of I j , we have that v 1 = ∆v/4r, v N = V max − ∆v/4r and v (l−1)r+1 = (l − 1)∆v. Thus, in order to compute the fundamental diagram of the δ model with any value of r, it is enough to compute the equilibria f 1 using r = 1 and then compute the flux with the formula above. In particular, using only f 1 , one may also compute the fundamental diagram of the δ model also in the limit r → ∞ with the formula The dashed blue line in all panels of Figure 12 shows the quantity Flux δ (∞) just defined. Note that in the case of the χ model, for each value of r, one has instead to compute the full equilibrium distribution with N = rT + 1 velocities.
When increasing r, we observe that the flux at ρ max approaches zero. This is because for ρ = ρ max , (f 1 ) 1 is the only non zero component at equilibrium, all vehicles travel at a velocity in the lowest speed class I 1 and the flux is therefore ∆v 4r (f 1 ) 1 . Similarly, in the free phase all vehicles travel at a velocity in the highest speed class I N and the flux is therefore (V max − ∆v 4r )(f 1 ) T +1 = (V max − ∆v 4r )ρ. The free-phase flux is therefore linear in ρ and its slope approaches V max when r → ∞.
In Figure 12, we observe that both models provide a sharp decrease in the flux, beyond the critical density, namely the value of the density marking the transition from free to congested flow. This phenomenon is well known in traffic modeling, and it is called capacity drop, see [29] and references therein. From Theorem 3.4 it is apparent that, for the δ model, the critical density corresponds to a bifurcation of the equilibrium solutions. In fact, one deduces that for P ≥ 1 2 the equilibrium distribution is f ∞ (v) = ρδ Vmax (v), which means that all vehicles travel at maximum speed. Only when P < 1 2 the lower speed classes begin to fill up. Thus, the physical concept of phase transition in traffic flow theory has a rigorous mathematical counterpart in the present model. Using the law given in (5) the value of ρ for which P = 1/2 is ρ c := 1 2 1 γ and then we may act on γ in order to change the critical density. For instance, see Figure 13 in which we plot three fundamental diagrams of the δ model with P as in (5) and γ = 1 ( * -markers), γ = 3/4 (×), γ = 1/4 ( ). Remark 10. Observe from Figure 12 that the critical density of the χ model approaches the critical density of the δ model when N → ∞. In fact, since the matrix A 1 χ → A 1 δ for r → ∞, we also have that (f χ r ) 1 → f δ 1 1 . More precisely, the analogous of (3.1.1) for the χ model is and the stable equilibrium is thus Figure 14 we show the fundamental diagrams of the χ model for r = 1 and r = 20, together with a few representative f j 's at equilibrium, as functions of ρ. In the left part, for r = 1, two phase transitions appear in the fundamental diagram (top left). Comparing with the bottom left plot, the origin of this phenomenon can be appreciated. A first transition occurs when the density becomes large enough to force a few drivers to brake: thus the second largest speed class I N −1 starts being populated (green dashed curve), while the fastest speed class begins to be depleted (red curve). A second transition occurs when some vehicles enter the lowest speed class (blue curve). This latter transition is the one that, when increasing r, moves towards the critical density ρ = 1/2, see Remark 10. The first phase transition is not observable for large r, because f N −1 is related to the velocity v N −1 → V max , as δv → 0, so that the transition of vehicles from I N to I N −1 is not enough to determine an abrupt change in the flux.

5.3.
Comparison with experimental data. Figure 15 shows the comparisons of the results produced by the δ model with experimental data published in [27]. In  Figure 15. Comparison between experimental data and the diagram resulting from the δ model, with ∆v = 1/3, P = 1 − ρ 3 /4 (left panel) and ∆v = 1/4, P = 1 − ρ 1 /4 (right panel). The experimental diagram is reproduced by kind permission of Seibold et al. [27].
the left figure we have tuned the critical density to reproduce the correct position of the phase transition. The experimental data are normalized and the fundamental diagram computed by the model is provided for all values of the density between 0 and ρ max , which corresponds to a situation in which all vehicles are bumper-to-bumper and still. The figure on the right stems from the observation that experimental data contain a residual movement even in the congested phase. Thus the bumper-to-bumper situation is never actually observed. Therefore in the figure on the right we also tune the maximum densityρ max actually observed, with ρ max < ρ max . In this case we obtain a very good agreement with experimental data. With the present model we do not reproduce the scattering of the data, but this can be explained keeping into account a mixture of two different populations of drivers and/or vehicles, as we have proposed in [25]. 6. Conclusions and perspectives. In this work we have studied two kinetic models for vehicular traffic based on a Boltzmann-like term describing binary interactions. We have assumed a continuous space of microscopic speeds and we have analyzed the space homogeneous case to study the asymptotic behavior of the distribution function together with the resulting flux-density diagrams.
Our models are characterized by a parameter ∆v, that has physical relevance and is related to the maximum speed variation in a unit of time. The models are defined by the transition probability of gaining a given velocity and they differ only in the modeling of the acceleration interaction.
First of all we have studied the case in which the resulting speed after an acceleration is obtained by a velocity jump from v * to v * + ∆v, where v * is the preinteraction speed. We have referred to this model as δ model and we have found a class of asymptotic distributions which is atomic with respect to the velocity variable. In other words it is a linear combination of Dirac delta functions centered in a finite number of velocities. The number T of delta functions contributing to the stable equilibrium distribution is controlled by the acceleration parameter through the relation ∆v = V max /T . This result means that the number of discrete velocities necessary to completely describe the equilibrium distribution function is implicitly determined by the acceleration parameter ∆v and therefore is small.
Instead, in the χ model we have prescribed the acceleration interaction in a way that is closer to the modeling given in [14], but again respecting the physical relevance of ∆v. In fact, we have assumed that the output speed after acceleration is uniformly distributed over the range [v * , v * + ∆v]. We have shown that the χ model with acceleration parameter ∆v gives a macroscopic behavior similar to the one provided by the simpler δ model with acceleration parameter ∆v/2, as it can be seen by studying the macroscopic properties of the two models and comparing their fundamental diagrams, notwithstanding the fact that the respective asymptotic distribution functions do not approach each other. Thus the χ model, despite the more sophisticated description of interactions, gives the same macroscopic information of the simpler and computationally much cheaper δ model, at least at equilibrium. Moreover we have proved that both models provide a bounded macroscopic acceleration, studying the evolution in time of the macroscopic velocity, and its relation with the parameter ∆v.
The results obtained in this work suggest that a small number of velocities is sufficient for the kinetic modeling of traffic. This is crucial to make kinetic modeling of complex traffic flows amenable to computations. Note also that here the acceleration remains controlled by ∆v, in contrast to models based on a lattice of velocities, in which the possible outcomes of an interaction and the acceleration of vehicles depend on the particular lattice chosen [6]. Moreover the complete knowledge of the equilibrium distribution is crucial to derive macroscopic models with a rich enough closure law resulting naturally as consequence of the microscopic interactions. Thus, without the need of prescribing heuristic speed-density relations, we obtain fundamental diagrams with a phase transition and a capacity drop as those occurring in experimental data.
Finally, the particular structure of the equilibrium distribution provided by the δ model allows one to generalize this framework to the case of multipopulation models, as in [24], in order to recover multivalued fundamental diagrams, see also [25].
Appendix A. Matrix elements for the discretization of the χ model. In order to compute the A j χ matrices, we observe that the gain operator of the δ model given in (9) differs from the gain operator of the χ model only in the terms proportional to P appearing in the equation (16). Therefore, we just show the terms resulting from When the terms above are integrated over the cells I 1 , we get For j = r + 2, . . . , N − r − 1, For j = N − r, For j = N − r + 1, . . . , N − 1, Finally, for j = N , The matrices A j χ can be formed by removing the underlined terms in (12) with r ∈ N and adding the contributions given in (A.1).