Uncertainty damping in kinetic traffic models by driver-assist controls

TL;DR

This paper introduces a kinetic traffic model incorporating microscopic uncertainties and demonstrates control strategies via driver-assist vehicles that reduce macroscopic variability, enhancing traffic predictability.

Contribution

It presents a novel kinetic model with uncertainty and develops control methods at the microscopic level to dampen uncertainty propagation in traffic flow.

Findings

Control strategies can significantly reduce traffic flow variability.

A relationship between driver-assist penetration rate and uncertainty damping is established.

Traffic predictability improves with increased implementation of control protocols.

Abstract

In this paper, we propose a kinetic model of traffic flow with uncertain binary interactions, which explains the scattering of the fundamental diagram in terms of the macroscopic variability of aggregate quantities, such as the mean speed and the flux of the vehicles, produced by the microscopic uncertainty. Moreover, we design control strategies at the level of the microscopic interactions among the vehicles, by which we prove that it is possible to dampen the propagation of such an uncertainty across the scales. Our analytical and numerical results suggest that the aggregate traffic flow may be made more ordered, hence predictable, by implementing such control protocols in driver-assist vehicles. Remarkably, they also provide a precise relationship between a measure of the macroscopic damping of the uncertainty and the penetration rate of the driver-assist technology in the traffic stream.

Tables

Table 1. Table 1: Choices of the Φ k subscript Φ 𝑘 \Phi_{k} ’s depending on the distribution Ψ Ψ \Psi of the uncertain parameter

Distribution $\Psi$	Polynomials ${\Phi }_k$	Support of the ${\Phi }_k$’s
Uniform	Legendre	Compact interval
Beta	Jacobi	Compact interval
Gamma	Laguerre	${ℝ}_{+}$
Binomial	Krawtchouk	$\mathbb{N}$

Equations

$$\begin{array}{rl}{v}^{\prime }& =v+\gamma I(v,v_{\ast };z)+D(v)\eta \\ v_{\ast }^{\prime }& =v_{\ast }.\end{array}$$

$$\mathrm{Prob}(z\le \bar{z})={\int }_{-\infty }^{\bar{z}}\Psi (z)dz.$$

$$I(v,v_{\ast };z)=P(\rho ;z)(1-v)+(1-P(\rho ;z))(P(\rho ;z)v_{\ast }-v),$$

$$P(\rho ;z)=(1-\rho {)}^z,z>0.$$

$$⟨\eta ⟩=0,⟨{\eta }^2⟩={\sigma }^2,$$

$$\begin{array}{rl}{v}^{\prime }& =v+\gamma [P(\rho ;z)(1-v)+(1-P(\rho ;z))(P(\rho ;z)v_{\ast }-v)]+D(v)\eta \\ v_{\ast }^{\prime }& =v_{\ast }\end{array}$$

$$(1-\gamma )v+D(v)\eta \ge 0,$$

$$\eta \ge c(\gamma -1),cD(v)\le v.$$

$$(\gamma -1)(1-v)+D(v)\eta \le 0,$$

$$\eta \le c(1-\gamma ),cD(v)\le 1-v.$$

$$|\eta |\le c(1-\gamma ),cD(v)\le \min \{v,1-v\},$$

$$\frac{d}{dt}{\int }_0^1\phi (v)f(t,v;z)dv=\frac{1}{2\tau }{\int }_0^1{\int }_0^1⟨\phi (v^{\prime })-\phi (v)⟩f(t,v;z)f(t,v_{\ast };z)dvd{v}_{\ast }$$

$$\frac{d}{dt}{\int }_0^{1}f(t,v;z)dv=0,$$

$$V(t;z):={\int }_0^{1}vf(t,v;z)dv,$$

$$\displaystyle\begin{aligned} \frac{dV}{dt}&=\frac{1}{2\tau}\int_{0}^{1}\int_{0}^{1}I(v,\,v_{\ast};\,z)f(t,\,v;\,z)f(t,\,v_{\ast};\,z)\,dv\,dv_{\ast}\\ &=\frac{\gamma}{2\tau}\left[P(\rho;\,z)(1-V)-\bigl{(}1-P(\rho;\,z)\bigr{)}^{2}V\right].\end{aligned}$$

$$V_{\infty}(\rho;\,z):=\frac{P(\rho;\,z)}{P(\rho;\,z)+{\bigl{(}1-P(\rho;\,z)\bigr{)}}^{2}}.$$

$${\bar{V}}_{\infty }(\rho ):={\mathbb{E}}_z(V_{\infty }(\rho ;z))={\int }_{{\mathbb{R}}_{+}}{V}_{\infty }(\rho ;z)\Psi (z)dz$$

$${\varsigma }_{\infty }^2(\rho ):={\mathrm{Var}}_z(V_{\infty }(\rho ;z))={\int }_{{\mathbb{R}}_{+}}{V}_{\infty }^2(\rho ;z)\Psi (z)dz-{\bar{V}}_{\infty }^2(\rho ).$$

$$\left\{(\rho ,q)\in [0,1]\times {\mathbb{R}}_{+}:q\in \left[\rho {\bar{V}}_{\infty }(\rho )-\rho {\varsigma }_{\infty }(\rho ),\rho {\bar{V}}_{\infty }(\rho )+\rho {\varsigma }_{\infty }(\rho )\right]\right\}$$

$$\mathrm{Prob}(z=z_k)={\alpha }_k\in [0,1],\sum _{k=1}^n{\alpha }_k=1,$$

$$\Psi (z)=\sum _{k=1}^n{\alpha }_k\delta (z-z_k),$$

$${\bar{V}}_{\infty }(\rho )=\sum _{k=1}^n{\alpha }_k{V}_{\infty }(\rho ;z_k),{\varsigma }_{\infty }^2(\rho )=\sum _{k=1}^n{\alpha }_k{V}_{\infty }^2(\rho ;z_k)-{\left(\sum _{k=1}^n{\alpha }_k{V}_{\infty }(\rho ;z_k)\right)}^2.$$

$$\Psi (z)=\frac{1}{b-a}\chi (a\le z\le b),$$

$${\bar{V}}_{\infty }(\rho )=\frac{2}{\sqrt{3}(b-a)\log (1-\rho )}{\left[\arctan \left(\frac{2x-1}{\sqrt{3}}\right)\right|}_{x={(1-\rho )}^a}^{x={(1-\rho )}^b}.$$

$${\varsigma }_{\infty }^2(\rho )=\frac{1}{3(b-a)\log (1-\rho )}{\left[\frac{\sqrt{x}-2}{x-\sqrt{x}+1}+\frac{2}{\sqrt{3}}\arctan \left(\frac{2\sqrt{x}-1}{\sqrt{3}}\right)\right|}_{x={(1-\rho )}^{2a}}^{x={(1-\rho )}^{2b}}-{\bar{V}}_{\infty }^2(\rho ).$$

$${\int }_0^1{\int }_0^1⟨\phi (v^{\prime })-\phi (v)⟩f_{\infty }(v;z)f_{\infty }(v_{\ast };z)dvd{v}_{\ast }=0$$

$$\gamma =ϵ,{\sigma }^2=\lambda ϵ,\tau =\frac{ϵ}{2},$$

$$\frac{d}{dt}{\int }_0^1\phi (v)f_{ϵ}(t,v;z)dv=\frac{1}{ϵ}{\int }_0^1{\int }_0^1⟨\phi (v^{\prime })-\phi (v)⟩f_{ϵ}(t,v;z)f_{ϵ}(t,v_{\ast };z)dvd{v}_{\ast },$$

$$\frac{d}{dt}{\int }_0^1\phi (v)f_{ϵ}(t,v;z)dv$$

$$+\frac{\lambda }{2}{\int }_0^1{\phi }^{\prime \prime }(v)D^2(v)f_{ϵ}(t,v;z)dv+R_{ϵ}(f_{ϵ},f_{ϵ})[\phi ](t;z),$$

$$\frac{d}{dt}{\int }_0^1\phi (v)f(t,v;z)dv$$

Full text

**Keywords: kinetic traffic model, Boltzmann-type equation with cut-off, Fokker-Planck limit, Uncertainty Quantification

**

Mathematics Subject Classification: 35Q20, 35Q70, 35Q84, 35Q93, 90B20

1 Introduction

In recent times, the challenge of vehicle automation has inspired new paradigms for the management and the governance of traffic and has imposed a deep reflection about the behavioural aspects involved in vehicle dynamics, see e.g. [3, 20]. On one hand, the emerging technologies constitute a potentially high innovation in the realm of the smart cities, for example with clear benefits in terms of emission reduction and of mitigation of road risk factors. On the other hand, the time horizon of their actual full implementation is still uncertain and questionable. For this reason, recent advances in the modelling of driver-assist/autonomus vehicles have focused, in particular, on the assessment of the effectiveness of Advanced Driver-Assistance Systems (ADAS) in a mixed scenario, with few automated vehicles embedded in a stream of human-manned vehicles [44]. In this context, it is particularly interesting to address the design of suitable control protocols, which are sufficiently robust to cope with some behavioural aspects of human drivers, and then to study their aggregate impact on the traffic stream, see e.g. [40].

In this work, we study a kinetic traffic model with structural uncertainty in the vehicle-to-vehicle interactions. Such an uncertainty, represented by a parameter whose value is not known deterministically, accounts, for example, for the heterogeneous response of different types of vehicles to the speed variations. Kinetic models for traffic flow have a quite long history, going back to the pioneering works [34, 36] up to more recent advances, such as e.g. [16, 23]. The kinetic approach has the clear advantage of linking, in a self-consistent way, the different scales of the problem: from the microscopic one of the interactions among the vehicles to the mesoscopic one of the aggregate distribution of observable quantities, such as the fundamental diagrams, and further also to the macroscopic one of the density waves flowing along a road [44]. In [37, 38] a kinetic approach has been developed to explain the scattering of the fundamental diagram of traffic. The latter is an experimentally measured relation between the traffic congestion and the vehicle flux, which typically features consistent macroscopic fluctuations at high vehicle density, see e.g. [22]. The mentioned works are based on multi-population/mixture models, which take into account the heterogeneous composition of the traffic stream by means of several kinetic equations, one for each class of vehicles. It is worth mentioning that also different approaches have been considered in the literature to model heterogeneous traffic flows. We recall, for instance, the car-following microscopic models [25, 27] and the systems of conservation laws with a flux function depending simultaneously on all the class densities [4, 9]. In the present work, we show that a kinetic explanation of the scattering of the fundamental diagram in the congested traffic regime can be effectively obtained also by introducing, in the microscopic interactions of a single-population model, an uncertain parameter accounting for the different reaction strengths of the various classes of vehicles. On one hand, this method leads to an increased dimensionality of the kinetic problem, owing to the necessity to handle the uncertain parameter. On the other hand, however, it avoids the necessity to deal with systems of kinetic equations, which would require to define several microscopic interaction rules and mesoscopic collision operators for the various classes of vehicles, thereby making the whole approach more amenable to additional modelling and analytical developments.

Along this line, in this work we exploit the potential offered by the uncertain kinetic setting to further address a control problem, aimed at reducing the scattering of the fundamental diagram. Specifically, we design two possible control protocols, which, consistently with the aforementioned ADAS technologies, are applied at the level of the microscopic interactions among the vehicles and, in particular, act in such a way to align the speed of the vehicles to a prescribed congestion-dependent value. Since, in a realistic scenario, all vehicles are not equipped with the ADAS technology, such controls are active only on a possibly small percentage of vehicles, corresponding to the so-called penetration rate. Our main result is that such microscopic control strategies are able to dampen the structural uncertainty present in the vehicle-to-vehicle interactions, thereby reducing the macroscopic scattering of the fundamental diagram in a way precisely linkable to the penetration rate of the ADAS technology. An immediate consequence of this is a more ordered, hence predictable, macroscopic flow of the vehicles, with clear advantages for the traffic governance at large scale. It is worth anticipating that this theoretical result holds for any possible statistical distribution of the uncertain parameter. Sticking to a purely model-based approach, in this paper we assume that such a distribution is given. However, we mention that, in the realm of Uncertainty Quantification (UQ), several tools from Bayesian inference have been developed to estimate the distribution of uncertain parameters in physical models, see e.g. [24]. As far as the traffic dynamics considered here are concerned, in the recent work [18] the statistical distribution of the uncertain parameter entering our interaction rules has been estimated from raw traffic data so as to fit the observable speed distributions.

From the methodological point of view, the proposed control setting takes advantage of recently introduced methods for the binary control of Boltzmann-type kinetic equations, see [1, 2] and also [43, 44] for the specific application to vehicular traffic. First, we derive an uncertain control, namely one obtained from a microscopic cost functional depending pointwise on the uncertain parameter of the interactions. As a second case, we derive a deterministic control, namely one obtained from a microscopic cost functional averaged with respect to the uncertain parameter. These two microscopic controls lead to two different Boltzmann-type models. The one associated with the uncertain control is a Maxwellian model, since the corresponding Boltzmann equation features a unitary interaction kernel. Conversely, the one associated with the deterministic control is a kinetic model with cut-off, because the corresponding Boltzmann equation needs a non-constant interaction kernel which discards some microscopic interactions possibly violating the physical bounds of the microscopic variable. We discuss the analytical properties of both models and, in particular, we show that, in the asymptotic regime of the quasi-invariant interactions [10, 41], they can be described by the same Fokker-Planck equation, which leads to explicitly computable and realistic steady states highlighting the reduction of the scattering of the fundamental diagram discussed above. We stress that the adopted kinetic approach makes possible a multiscale investigation of the propagation of both the structural uncertainty and the effect of the control across the different scales of the problem. We also propose several numerical tests, which, by appealing to specific numerical methods for Uncertainty Quantification for kinetic and mean field equations, cf. e.g. [6, 13, 19, 21, 49, 50], visualise the theoretical results and allow us to consider various probability distributions of the uncertain parameter, spanning also interesting cases in which explicit computations are not possible.

In more detail, the paper is structured as follows: in Section 2, we discuss a basic microscopic model of the interactions among the vehicles without control and we stress, in particular, the role of the uncertain parameter. Then we pass to a Boltzmann-type kinetic description, whence we obtain explicitly the asymptotic trend of the mean speed, which we use to define the fundamental diagram together with its scattering induced by the uncertain parameter. Finally, in the quasi-invariant interaction limit, we recover a Fokker-Planck-type description, whence we deduce an explicit form of the equilibrium speed distribution along with its uncertainty. In Section 3, we introduce the two microscopic control strategies to be applied to the interactions among the vehicles. Specifically, in Section 3.1, we discuss the uncertain control, which leads to a Boltzmann-type kinetic description for Maxwellian-like particles and allows us to prove a precise result concerning the reduction of the scattering of the fundamental diagram in terms of the penetration rate of the control strategy. In Section 3.2, we discuss instead the deterministic control, which is more realistically implementable in practice but requires to deal with a more difficult Boltzmann-type kinetic description with cut-off. However, we prove that, in the quasi-invariant interaction regime, such a description yields the same limit equations for both the mean speed and the statistical speed distribution as the previous description, whence we recover the validity of all the results proved before. In Section 4, we present several numerical tests supporting the theoretical findings of the previous sections. Finally, in Section 5 we summarise the results of the paper and we briefly outline possible research developments.

2 Scattering of the fundamental diagram from uncertain binary interactions

2.1 Description of the microscopic interactions with uncertainty

In accordance with the general approach of the kinetic theory, a kinetic description of traffic flow is based on the identification of microscopic interaction rules for pairs of vehicles. We assume, in particular, that these interactions modify the speed of the vehicles in consequence of accelerations and decelerations. Therefore, we characterise the microscopic state of a generic vehicle by means of a variable $v\in[0,\,1]$ representing its (dimensionless) speed. Denoting by $v_{\ast}\in[0,\,1]$ the speed of the leading vehicle, we describe the speed variation in a binary interaction as

[TABLE]

In (1), $\gamma>0$ is a proportionality parameter and $I$ is the interaction function, which describes the acceleration/deceleration dynamics mentioned above. This function depends on the pre-interaction speeds $v$, $v_{\ast}$ of the interacting vehicles and on an uncertain parameter, viz. random variable, $z\in\mathbb{R}$ with known probability distribution $\Psi=\Psi(z):\mathbb{R}\to\mathbb{R}_{+}$, i.e.

[TABLE]

Such an uncertain parameter represents a structural uncertainty in the binary rule modelled by $I$, due to the fact that the physics of the interactions among the vehicles is inevitably partly heuristic.

Getting inspiration from [44], we consider the following interaction function:

[TABLE]

where $P(\rho;\,z)\in[0,\,1]$ is the probability of acceleration. The function (2) expresses the fact that with probability $P$ the $v$-vehicle accelerates towards the maximum speed, while with probability $1-P$ it adapts to the fraction $P$ of the speed of the $v_{\ast}$-vehicle. The probability of acceleration depends on the (dimensionless) density of the vehicles $\rho\in[0,\,1]$, in such a way that the higher $\rho$ the lower $P$, because a dense traffic hinders accelerations. Moreover, $P$ depends also on the uncertain parameter $z$. A possible form is:

[TABLE]

From Figure 1, we observe that the mapping $\rho\mapsto P(\rho;\,z)$ expresses the same qualitative trend ($P$ decreasing with $\rho$) for all $z>0$. Nevertheless, it may induce quantitatively different dynamics, because such a trend is either concave, linear or convex depending on the specific value of $z$. It is therefore reasonable to regard the exponent $z\in\mathbb{R}_{+}$ as an uncertain parameter, considering that there are apparently no a priori motivations to opt for a specific value of $z$ in the heuristic model (3). Moreover, since $\rho\in[0,\,1]$, we notice that the higher $z$ the lower the probability of acceleration at all densities $\rho>0$. This suggests that the values of $z$ may be associated with different classes of vehicles characterised by a different motility in the traffic stream. We will come back to this interpretation in Section 2.2.

Returning to (1), the term $D(v)\eta$ describes stochastic fluctuations caused by the intrinsic randomness of the driver behaviour. Specifically, $\eta\in\mathbb{R}$ is a centred random variable, i.e. one with zero mean:

[TABLE]

where $\langle\cdot\rangle$ denotes the average with respect to the probability distribution of $\eta$ and $\sigma>0$ is the standard deviation of $\eta$. Moreover, $D:[0,\,1]\to\mathbb{R}_{+}$ is a speed-dependent non-negative diffusion coefficient modulating the amplitude of the stochastic fluctuation. We stress that, unlike $z$, the random variable $\eta$ does not describe a structural uncertainty of the model but rather a stochastic completion of the deterministic dynamics expressed by the interaction function $I$.

Finally, we notice that (1) prescribes that the $v_{\ast}$-vehicle does not change speed during an interaction. This is because vehicle interactions are normally anisotropic, the leading vehicle being unperturbed by the rear vehicle.

2.1.1 Physical admissibility of the interaction rules

In order to be physically admissible, the interaction rules

[TABLE]

have to guarantee $v^{\prime},\,v_{\ast}^{\prime}\in[0,\,1]$ for every $v,\,v_{\ast}\in[0,\,1]$. This is actually obvious for $v_{\ast}^{\prime}$, while it is more delicate for $v^{\prime}$.

Let us assume $\gamma\in[0,\,1]$. We begin by observing that, since $0\leq P(\rho;\,z)\leq 1$ and $v_{\ast}\geq 0$, a sufficient condition for $v^{\prime}\geq 0$ is

[TABLE]

which is certainly satisfied if there exists a constant $c>0$ such that

[TABLE]

Likewise, since $v_{\ast}\leq 1$, a sufficient condition for $v^{\prime}\leq 1$ is

[TABLE]

which is satisfied if

[TABLE]

On the whole, the physical admissibility of (5) is guaranteed by the following sufficient conditions:

[TABLE]

where $c>0$ is an arbitrary constant. These conditions imply that the stochastic fluctuation $\eta$ is bounded and that $D(0)=D(1)=0$.

2.2 Boltzmann-type aggregate description and traffic diagrams

A statistical description of the aggregate dynamics resulting from the superposition of many binary interactions (5) may be obtained by introducing the distribution function $f=f(t,\,v;\,z)$, where $t\geq 0$ is the time, see [8]. Specifically, $f$ is such that $f(t,\,v;\,z)dv$ gives the probability that, at time $t$, a vehicle travels with a speed comprised between $v$ and $v+dv$, given the uncertain parameter $z$. Since, under (6), all interactions (5) are physically admissible, the distribution function $f$ evolves according to the following Boltzmann-type equation for Maxwellian-like particles, here written in weak form, see [32]:

[TABLE]

for every observable quantity $\varphi:[0,\,1]\to\mathbb{R}$, namely any quantity which may be expressed as a function of the microscopic state $v$ of the vehicles. This equation states that the time variation of the mean of $\varphi$ (left-hand side) is due to the mean of the average variation of $\varphi$ in a representative binary interaction (right-hand side). The parameter $\tau>0$ at the right-hand side is a relaxation time or, in other words, $\frac{1}{\tau}$ is the interaction frequency.

We have stressed that $f$ is parametrised by the uncertain parameter $z$, because its evolution depends, among other things, on the uncertain interaction function $I$ contained in $v^{\prime}$. Consequently, (7) is a stochastic kinetic equation, whose solution may be regarded, for fixed $t$, $v$, as a random variable given as a function of the uncertain parameter $z$.

Letting $\varphi(v)=1$ in (7), we discover

[TABLE]

therefore our kinetic model conserves in time the total mass of the system. We notice that the interpretation of the mapping $v\mapsto f(t,\,v;\,z)$ as a probability density function on $[0,\,1]$ for all $t$, $z$ is compatible with this property of equation (7). With $\varphi(v)=v$ and recalling (4), we obtain instead from (7) the evolution of the mean speed

[TABLE]

i.e.:

[TABLE]

This equation can be solved explicitly to find the mapping $t\mapsto V(t;\,z)$ parametrised by $\rho$. Here, we are in particular interested in the asymptotic value reached by $V$ for $t\to+\infty$, say $V_{\infty}$, which describes the mean speed emerging when interactions are in equilibrium:

[TABLE]

Notice that, due to the uncertain parameter $z$, this is in turn a stochastic quantity. We may compute the expectation of $V_{\infty}$ with respect to $z$ as

[TABLE]

and its variance as

[TABLE]

Next, we observe that we may use the mapping $\rho\mapsto\rho\bar{V}_{\infty}(\rho)$ to define the fundamental diagram of traffic, namely the equilibrium relationship between the traffic density and the macroscopic flux of the vehicles. Together with its $z$-standard deviation $\rho\mapsto\rho\varsigma_{\infty}(\rho)$, it produces the following set

[TABLE]

in the density-flux plane, where most of the random flux values $\rho V_{\infty}(\rho;\,z)$ lie. The set (12) explains the scattering of the fundamental diagram, typically found in empirical measurements of the flow of vehicles, as the result of the superposition of different microscopic dynamics produced by different values of $z$ and weighted by the corresponding probability measure $\Psi(z)dz$.

For a general probability distribution $\Psi$, the exact computation of (10), (11) may be non-trivial and one often needs to rely on numerical quadrature formulas. However, for special classes of probability distributions $\Psi$, such as those considered below, analytical results can be obtained.

Let us consider the case in which $z\in\{z_{1},\,z_{2},\,\dots,\,z_{n}\}\subset\mathbb{R}_{+}$ is a discrete random variable with law

[TABLE]

so that

[TABLE]

where $\delta(z-z_{k})$ is the Dirac delta distribution centred in $z_{k}$. Then

[TABLE]

We observe that $V_{\infty}(\rho;\,z_{k})$, $k=1,\,\dots,\,n$, is the result of the microscopic dynamics (5) with $z=z_{k}$. If we interpret the $z_{k}$’s as characteristic values of certain classes of vehicles (for instance, cars, lorries, motorcycles, …) which may travel along the road, the formulas above show that $\bar{V}_{\infty}(\rho)$, $\varsigma_{\infty}^{2}(\rho)$ originate from the average superposition of the macroscopic dynamics produced by each of such classes. Within this interpretation, the $\alpha_{k}$’s can be understood as the proportions of the various classes of vehicles present in the traffic stream. This provides a point of view on multi-class traffic models based on structural uncertainties in the microscopic composition of the traffic “mixture”. We mention that other multi-class traffic models are already present in the literature, see e.g. [4, 37], which also explain the scattering of the fundamental diagram by appealing to similar physical motivations but different mathematical formalisations.

Let us assume, for instance, that $z$ takes only the two values $z_{1}=1$, with probability $\alpha_{1}$, and $z_{2}=3$, with probability $\alpha_{2}$. The largest value of $z$, i.e. $z_{2}$, may represent e.g., lorries, which accelerate less especially at high traffic density, cf. (3). Conversely, the smallest value of $z$, i.e. $z_{1}$, may represent e.g., cars. We consider two different compositions of the traffic stream: first, one with $\alpha_{1}=70\%$ of cars and $\alpha_{2}=30\%$ of lorries; then, the opposite one with $\alpha_{1}=30\%$ of cars and $\alpha_{2}=70\%$ of lorries. In Figure 2, we see that both the fundamental diagram and the region (12) change realistically according to the composition of the traffic stream. In both cases, we notice that the scattering of the fundamental diagram is quite limited at low density, i.e. in the so-called free flow regime, when the flux grows almost linearly with $\rho$. It becomes more marked at higher density, i.e. in the so-called congested flow regime, when the flux decreases non-linearly with $\rho$. This is very nicely in agreement with the typical experimental observations, see e.g. [22, Chapter 2] and also [15, 39].

As a second example, let us consider $z\sim\mathscr{U}([a,\,b])$, i.e. the case in which $z$ is a continuous random variable uniformly distributed in the interval $[a,\,b]$ with $0\leq a<b$. Hence

[TABLE]

where $\chi$ denotes the characteristic function (specifically, $\chi(a\leq z\leq b)=1$ if $a\leq z\leq b$ while $\chi(a\leq z\leq b)=0$ otherwise). Using (3) and (9), from (10) we compute:

[TABLE]

Also the variance $\varsigma_{\infty}^{2}(\rho)$ can be given an explicit representation, indeed from (11) we obtain:

[TABLE]

Figure 3 shows the area defined by (12) for $z\sim\mathscr{U}([1,\,3])$, i.e. $a=1$ and $b=3$ in the formulas above. Also in this case, we notice that the computed scattering of the fundamental diagram reproduces correctly the qualitative empirical features mentioned in the previous example.

Remark 2.1*.*

In the two examples discussed above, the statistical distribution $\Psi$ has been postulated. As already mentioned in the Introduction, the problem of estimating $\Psi$ is beyond the scope of this paper. Nevertheless, recent efforts in this direction based on the analysis of real traffic data, see e.g. [17, 18], show that the types of distribution and the range of $z$ considered here are consistent with the results obtained from experimental measurements.

2.3 Fokker-Planck description and equilibria

Besides the fundamental diagram of traffic, from the kinetic description of the microscopic dynamics (5) we may also obtain information on the statistical distribution of the speed of the vehicles emerging when interactions are in equilibrium, i.e. for large times ($t\to+\infty$ in the limit). Such a distribution corresponds to the Maxwellian distribution of the classical kinetic theory of gases.

By definition, the equilibrium distribution, say $f_{\infty}=f_{\infty}(v;\,z)$, makes the right-hand side of (7) vanish:

[TABLE]

for all observable quantities $\varphi$. Nevertheless, it is in general difficult to determine $f_{\infty}$ by tackling directly this equation, at least from the analytical point of view. If, however, the interactions (5) are quasi-invariant, i.e. if they produce each time a small change of speed, then it is possible to approximate asymptotically the Boltzmann-type integro-differential equation (7) with a Fokker-Planck partial differential equation, which may allow for an explicit computation of the steady states.

Quasi-invariant interactions are reminiscent of the grazing collisions introduced in the classical kinetic theory [45, 46]. In our context, they amount to assuming $\gamma,\,\sigma^{2}\ll 1$ in (5), so that both the deterministic part of the interactions and the stochastic fluctuations are small. To balance such a smallness and observe asymptotic trends, the interaction frequency has to be increased accordingly. We obtain such an effect by means of the following scaling:

[TABLE]

where $0<\epsilon\ll 1$ is a scaling parameter and $\lambda>0$ a proportionality constant. From (7) we obtain therefore

[TABLE]

where we have stressed that the solution $f$ depends now also on the scaling parameter $\epsilon$.

From now on, we take $\varphi$ smooth, specifically $\varphi\in C^{\infty}_{c}(0,\,1)$. We assume that, for every $\epsilon>0$ and every $z\in\mathbb{R}_{+}$, (14) admits a unique solution $f_{\epsilon}(\cdot,\,\cdot;\,z)\in C(\mathbb{R}_{+};\,L^{1}(0,\,1))$ which is a probability density function in the variable $v$, cf. [14, Appendix A]. Moreover, we assume that there exists a probability density function $f(\cdot,\,\cdot;\,z)\in C(\mathbb{R}_{+};\,L^{1}(0,\,1))$ such that, up to subsequences, $f_{\epsilon}(\cdot,\,\cdot;\,z)\to f(\cdot,\,\cdot;\,z)$ in $C(\mathbb{R}_{+};\,L^{1}(0,\,1))$ when $\epsilon\to 0^{+}$.

Expanding $\varphi(v^{\prime})$ in Taylor’s series up to the third order about $v$ (because $v^{\prime}\approx v$ in the quasi-invariant regime), we get

[TABLE]

where $I$ is given by (2) and $R_{\epsilon}(f_{\epsilon},\,f_{\epsilon})[\varphi]$ is a remainder which depends on $\varphi^{\prime\prime\prime}$ and on $(v^{\prime}-v)^{3}$, cf. [44]. In the limit $\epsilon\to 0^{+}$, it can be proved that $R_{\epsilon}(f_{\epsilon},\,f_{\epsilon})[\varphi]\to 0$ for all $\varphi\in C^{\infty}_{c}([0,\,1])$, see again [44] for details. Consequently, in view of the assumptions stated above, we find that the asymptotic equation satisfied by $f$ is

[TABLE]

i.e., integrating by parts the right-hand side and invoking the compactness of the support and the arbitrariness of $\varphi$:

[TABLE]

where we have substituted the expression (2) of $I$.

Equation (15) is the announced Fokker-Planck asymptotic model, which describes the trend of the system in the quasi-invariant regime. Notice that (15) is a Fokker-Planck equation with non-constant coefficients, due to the fact that the mean speed $V(t;\,z)$ is not conserved in time by the interactions (5). The stationary solution $f_{\infty}$ to (15), which provides an asymptotic approximation of the Maxwellian of (7), solves

[TABLE]

where $V_{\infty}$ is the asymptotic mean speed (9). In particular, if we choose111The function (16) does not actually comply with the requirement in (6), due to the vertical tangents at $v=0,\,1$. Nevertheless, such a $D$ may be obtained as the uniform limit for $\epsilon\to 0^{+}$ of e.g., the sequence of functions $D_{\epsilon}(v):=\sqrt{[(1+\epsilon)v(1-v)-\frac{\epsilon}{4}]_{+}}$, where $[\cdot]_{+}$ denotes the positive part. It can be checked that $D_{\epsilon}$ complies with (6) with $c=\sqrt{\frac{\epsilon}{1+\epsilon}}$, see [41] for details. Therefore, the use of (16) in (15) is justified after performing the quasi-invariant limit.

[TABLE]

we obtain that the unique stationary solution with unitary mass is

[TABLE]

where $\operatorname{Beta}(\cdot,\,\cdot)$ is the beta function. We notice that (17) is a beta probability density function with mean $V_{\infty}(\rho;\,z)$, as expected, and variance $\frac{\lambda V_{\infty}(\rho;\,z)(1-V_{\infty}(\rho;\,z))}{2+\lambda}$. Interestingly, beta probability densities have been found to fit particularly well the experimental data on the speed distribution of the vehicles [26, 28].

Actually, (17) is a family of equilibrium distributions parametrised by $z$. In order to quantify the intrinsic uncertainty in this description of the equilibrium, we may refer to the mean equilibrium distribution:

[TABLE]

and to its variance

[TABLE]

Notice, from (10), that

[TABLE]

Figures 4, 5 show $\bar{f}_{\infty}$, together with the uncertainty bands determined by its standard deviation $\pm\sqrt{\operatorname{Var}_{z}(f_{\infty}(\cdot;\,z))}$, for various discrete and continuous probability distributions of $z$. The curves have been computed numerically by means of suitable quadrature formulas in $z$, starting from the analytical expression (17). We observe that, depending on the distribution of $z$, in the mean distribution $\bar{f}_{\infty}$ several local peaks may appear, which suggest a slightly multi-modal average trend of the equilibrium speed distribution. Such a trend is actually documented in experimental data [18, 26, 28].

3 Uncertainty damping by driver-assist control

In this section, we tackle the problem of damping the scattering of the fundamental diagram, i.e., owing to (12), of reducing the variance $\varsigma_{\infty}^{2}(\rho)$. Although the scattering is a macroscopic manifestation of the traffic, its origin relies in the structural uncertainty affecting the microscopic interactions (5). Therefore, we will proceed to control a priori (some of) the microscopic interactions (5), assessing then the aggregate impact of such a control, rather than trying to constrain a posteriori the macroscopic flow of the vehicles as a whole. Indeed, in our view, actually feasible ways of controlling a multi-agent system, such as vehicular traffic, make use of bottom-up strategies, which target the simple microscopic interactions, rather than of strategies pointing directly to the aggregate behaviour. In the context of vehicular traffic, such bottom-up strategies refer to the control of the local responses given by the driver-vehicle system to the external stimuli coming from other vehicles by means of Advanced Driver-Assistance Systems (ADAS).

Let us consider, therefore, the following controlled binary interaction:

[TABLE]

where $I$ is given by (2), $u$ is a control to be determined from the optimisation of a suitable functional and $\Theta\in\{0,\,1\}$ is a Bernoulli random variable taking into account whether the randomly selected $v$-vehicle is or is not equipped with the ADAS technology. In particular, the law of $\Theta$ is

[TABLE]

where $p\in[0,\,1]$ is the so-called penetration rate, i.e. the percentage of ADAS vehicles in the traffic stream.

We propose, in the following, two different functionals, whose minimisation yields an instantaneous optimal control $u^{\ast}$ to be used in (18). These functionals correspond, in principle, to different control strategies, which however, as we will see, lead to the same macroscopic effect. Remarkably, one of such strategies has an eminently theoretical interest, whereas the other is closer to an actual implementability.

3.1 Pointwise uncertainty control

Let $v_{d}(\rho)\in[0,\,1]$ be a desired (or recommended) speed determined as a function of the vehicle density $\rho$ and suggested to the ADAS vehicles by e.g., some sensors along the road, which are able to detect the level of traffic congestion. We assume $0\leq v_{d}(\rho)\leq 1$ for all $\rho\in[0,\,1]$. Further physically meaningful characteristics of the desired speed are that the mapping $\rho\mapsto v_{d}(\rho)$ be non-increasing and that $v_{d}(0)=1$, $v_{d}(1)=0$. A prototypical example may be $v_{d}(\rho)=1-\rho$.

Automated traffic flow dynamics have been modelled recently as cruise controls implemented either in the vehicle-to-vehicle interactions, see e.g. [15, 29, 35, 44], or directly in the aggregate traffic dynamics [11, 12]. In general, the effect of such controls is the stabilisation of the traffic flow depending on the penetration rate of the automated vehicles, as also testified by the field experiment reported in [40]. Taking inspiration from this idea, an optimal control $u^{\ast}$ may be thought of as aimed at keeping the post-interaction speed $v^{\prime}$ of an ADAS vehicle as close as possible to $v_{d}(\rho)$. This is obtained by minimising e.g. the following functional:

[TABLE]

subject to the constraint (18). Notice that $J_{1}$ depends pointwise on the uncertain parameter $z$ through $v^{\prime}$. Moreover, $\nu u^{2}$, with $\nu>0$, is a quadratic penalisation term, which expresses the cost of the control.

The minimisation of $J_{1}$ may be carried out by forming the Lagrangian:

[TABLE]

where $\mu\in\mathbb{R}$ is the Lagrange multiplier, and then by considering the optimality conditions:

[TABLE]

whence we deduce the optimal control

[TABLE]

Notice that, for each $z\in\mathbb{R}_{+}$, this is a feedback control, because it is expressed as a function of the pre-interaction speeds $v$, $v_{\ast}$. Plugged into (18), it produces the controlled binary interactions

[TABLE]

3.1.1 Physical admissibility of the interactions (20)

The physical admissibility of the interactions (20) requires $v^{\prime}\in[0,\,1]$ for all $v,\,v_{\ast}\in[0,\,1]$. Considering the expression (2) for $I$, together with the observations made in Section 2.1.1 plus $v_{d}(\rho)\geq 0$ and $\Theta\leq 1$, we see that a sufficient condition for $v^{\prime}\geq 0$ is

[TABLE]

which is certainly satisfied if there exists a constant $c>0$ such that

[TABLE]

Likewise, since $v_{d}(\rho)\leq 1$, a sufficient condition for $v^{\prime}\leq 1$ is

[TABLE]

which is certainly satisfied if

[TABLE]

On the whole, (20) is physically admissible if there exists $c>0$ such that

[TABLE]

Again, this implies that the random variable $\eta$ has to be bounded and that the diffusion coefficient $D$ vanishes at $v=0,\,1$. Notice that if $\nu\to+\infty$, i.e. if the control is so penalised that the optimum is not to control the interactions ($u^{\ast}=0$, cf. (19)), then we recover exactly (6).

3.1.2 Boltzmann-type kinetic description

Since, under (21), all the interactions (20) are physically admissible, the evolution of the distribution function $f=f(t,\,v;\,z)$ is described by a Boltzmann-type equation for Maxwellian-like particles analogous to (7). The only difference is that now, in the “collisional” term appearing on the right-hand side of (7), we need to take into account a further average with respect to the random variable $\Theta$. Thus we write:

[TABLE]

where $v^{\prime}$ is given by (20).

In particular, we are interested in the trend of the mean speed $V=V(t;\,z)$, which we obtain with $\varphi(v)=v$ considering furthermore that $\Theta^{2}$ in (20) is a Bernoulli random variable with parameter $p$, because so is $\Theta$ by definition. Hence:

[TABLE]

In order to better characterise the role of the various parameters in the large time trend of the mean speed, it is useful to refer to the quasi-invariant interaction regime introduced in (13). In this case, to observe in the limit also a balanced contribution of the controlled part of the interactions, we scale also the control penalisation $\nu$ as

[TABLE]

where $\kappa>0$ is a proportionality parameter which expresses the strength of the control in the quasi-invariant scale. Then, we obtain that the limit equation satisfied by $V$ for $\epsilon\to 0^{+}$ is

[TABLE]

where

[TABLE]

is what we call the effective penetration rate. In particular, the equilibrium mean speed is

[TABLE]

Since $0\leq P(\rho;\,z),\,v_{d}(\rho)\leq 1$, we have222We use the notation $a\lesssim b$ to mean that there exists a constant $C>0$, whose exact value is unimportant, such that $a\leq Cb$.

[TABLE]

which shows that the instantaneous control (19) has the effect of keeping the average traffic stream close to the suggested optimal speed. From this, we easily deduce the important consequence:

[TABLE]

meaning that the control reduces indeed the scattering of the fundamental diagram caused by the uncertain parameter $z$. The magnitude of such a reduction is estimated by the inverse of ${(p^{\ast})}^{2}$, hence it is higher for either a higher penetration rate $p$ or a lower control penalisation $\kappa$.

Interestingly, this result does not actually depend on the specific form (2) of the interaction function $I$. A fruitful generalisation is enlightened by the following theorem:

Theorem 3.1**.**

In (20), let the interaction function $I$ be bounded, i.e. let a constant $I_{\max}>0$ exist such that

[TABLE]

Assume, moreover, that $I$ is such that the interactions (20) are all physically admissible in the sense discussed in Section 3.1.1. Then

[TABLE]

Proof.

We observe that it suffices to prove an estimate like (26), for then the thesis follows immediately from (27), which holds independently of the specific form of $I$.

Since all interactions (20) are physically admissible by assumption, we may use the Boltzmann-type equation (22). For $\varphi(v)=v$ we find then

[TABLE]

which, in the quasi-invariant limit, becomes

[TABLE]

Since $v_{d}(\rho)$ is time-independent, we have $\frac{dV}{dt}=\frac{d}{dt}\left(V-v_{d}(\rho)\right)$, therefore we may write

[TABLE]

whence, integrating in time,

[TABLE]

with $V_{0}(z):=V(0;\,z)$. From the boundedness of $I$ we further obtain

[TABLE]

thus we discover:

[TABLE]

whence finally

[TABLE]

and we are done. ∎

In Figure 6, we show the fundamental diagram $\rho\bar{V}_{\infty}(\rho)$, cf. (25), and its standard deviation $\rho\varsigma_{\infty}(\rho)$, cf. (11), for both a discrete (binomial) and a continuous (uniform) probability distribution of the uncertain parameter $z$ and increasing values of the effective penetration rate $p^{\ast}$. The curves have been obtained numerically by means of suitable quadrature formulas in $z$, starting from the analytical expression (25). As predicted by (27), and more in general by Theorem 3.1, the results display a clear damping of the uncertainty at the macroscopic level, thanks to the action of the microscopic control. This is particularly evident in Figure 6(c)-(d), which may be compared with Figure 3 illustrating the related uncontrolled case.

3.1.3 Fokker-Planck description and equilibria

Similarly to Section 2.3, we may obtain more details on the large time trend emerging from the controlled interactions (20) by studying the quasi-invariant limit of the full Boltzmann-type equation (22). We recall that, in this case, this amounts to considering the quasi-invariant scaling (13)-(23). Proceeding like in Section 2.3, we determine the following Fokker-Planck equation with non-constant coefficients:

[TABLE]

where $V$ is the solution to (24), whence, at the steady state,

[TABLE]

with $V_{\infty}$ given by (25). The unique distribution with unitary mass solving this equation for $D$ like in (16) is

[TABLE]

namely again a beta probability density function parametrised by $z$ and incorporating the action of the microscopic control through the parameter $p^{\ast}$.

In Figure 7, we show the mean equilibrium distribution $\bar{f}_{\infty}$, together with the uncertainty bands produced by its standard deviation $\pm\sqrt{\operatorname{Var}_{z}(f_{\infty}(v;\,z))}$, for a discrete (binomial) and a continuous (uniform) probability distribution of $z$ and increasing values of the effective penetration rate $p^{\ast}$. The curves have been obtained numerically, like in the previous cases, by means of suitable quadrature formulas in $z$, starting from the analytical expression (29). A comparison of Figure 7(a)-(b) with Figure 4 and of Figure 7(c)-(d) with Figure 5(a) clearly shows that, by increasing $p^{\ast}$, the mean value of the expected equilibrium distribution $\bar{f}_{\infty}$ tends to the desired speed $v_{d}(\rho)$, cf. (26), and, moreover, that the uncertainty affecting $\bar{f}_{\infty}$ is globally dampened with respect to the uncontrolled case.

3.2 Uncertainty control on average

The instantaneous optimal control (19) has a theoretical interest, as confirmed by Theorem 3.1, but it may not be straightforwardly implementable in practice. Indeed, it is an uncertain control, because it depends on the specific value of the uncertain parameter $z$. In this section, we focus on the construction of a deterministic control, i.e. one independent of $z$, and we prove that it actually performs analogously to the uncertain one as far as the damping of the $z$-uncertainty is concerned.

To this purpose, we consider the following new functional to be minimised:

[TABLE]

and we form the corresponding Lagrangian:

[TABLE]

The main difference with respect to the case considered in Section 3.1 is that here we perform an optimisation on average, rather than pointwise, in $z$. Consequently, we also require that the constraint given by the interaction dynamics (18) be satisfied on average, as it is clear from the expression of $\mathscr{L}_{2}$ above.

The optimality conditions take now the form:

[TABLE]

whence the instantaneous optimal control turns out to be

[TABLE]

We observe that this control coincides with the $z$-average of the uncertain control (19). Moreover, also (30) is a feedback control.

Plugging into (18), we obtain the interaction rules

[TABLE]

In this case, it is not as straightforward as before to individuate sufficient conditions on $\eta$ and $D(v)$, like (6) or (21), which guarantee the physical admissibility of these interactions. The reason is that, when attempting to ensure $v^{\prime}\in[0,\,1]$ by bounding the random fluctuation and the diffusion coefficient, it is not easy to compare effectively $I(v,\,v_{\ast};\,z)$ and $\mathbb{E}_{z}(I(v,\,v_{\ast};\,z))$ for a possibly generic probability distribution of $z$, not even if we refer to the specific interaction function (2). Consequently, we have to admit that, in general, some interactions (31) may not be physically admissible (i.e., they may produce either $v^{\prime}<0$ or $v^{\prime}>1$) and need therefore to be discarded from the statistical description of the system.

3.2.1 Boltzmann-type description with cut-off

In order to discard the possibly unphysical microscopic interactions (31), we consider the following Boltzmann-type equation with cut-off:

[TABLE]

which differs from (22) because of the non-constant interaction kernel

[TABLE]

This makes the study of (32) more challenging than in the previous cases, because it introduces additional non-linearities. For instance, the trend of the mean speed $V=V(t;\,z)$ is given, for $\varphi(v)=v$, by the equation

[TABLE]

whose right-hand side can no longer be easily written in terms of $V$ itself, like in (8) or (24). To tackle this case, we may fruitfully work in the quasi-invariant regime, taking inspiration from the technique presented in [10]. In particular, we will show in the next Section 3.2.2 that the Fokker-Planck equation originating from (32) in such a regime is still (28), which implies in particular that the properties of the equilibria remain unchanged with respect to the case of the pointwise uncertain control discussed in Section 3.1.

3.2.2 Fokker-Planck asymptotics

Similarly to (14), under the quasi-invariant scaling (13)-(23) equation (32) may be rewritten as

[TABLE]

We make the same assumptions as in Section 2.3 concerning the existence and uniqueness of the solution to this equation for every $\epsilon>0$ and the convergence, up to subsequences, of $f_{\epsilon}(\cdot,\,\cdot;\,z)$ to some $f(\cdot,\,\cdot;\,z)\in C(\mathbb{R}_{+};\,L^{1}(0,\,1))$ as $\epsilon\to 0^{+}$. Moreover, from now on we take $\varphi\in C^{\infty}_{c}(0,\,1)$ and, inspired by (6) and (21), we assume $cD(v)\leq\min\{v,\,1-v\}$ for some $c>0$ along with $D(v)>0$ for all $v\in(0,\,1)$.

In the quasi-invariant limit, the term $A_{\epsilon}(f_{\epsilon},\,f_{\epsilon})[\varphi]$ produces the right-hand side, in weak form, of the Fokker-Planck equation (28). Indeed, the additional term with $\mathbb{E}_{z}(I(v,\,v_{\ast};\,z))$ contained in the interaction rules (31), which is absent in (20), can be easily shown to be of order $\epsilon$, hence infinitesimal in the limit.

Therefore, the goal is to show that the remainder $R_{\epsilon}(f_{\epsilon},\,f_{\epsilon})[\varphi]$ vanishes for $\epsilon\to 0^{+}$. To this purposes, we observe preliminarily that

[TABLE]

and that, in view of the scaling (13) of $\sigma^{2}$, we may represent $\eta=\sqrt{\lambda\epsilon}Y$, where $Y$ is a random variable with $\langle Y\rangle=0$ and $\langle Y^{2}\rangle=1$. We assume moreover that $\eta$, hence also $Y$, has finite third order moment, i.e. $\langle\lvert Y\rvert^{3}\rangle<+\infty$. Then, upon scaling the other parameters of the interactions (31) according to (13)-(23) and expanding $\varphi(v^{\prime})$ in Taylor’s series about $v$, we estimate333Here and henceforth we use the notation $a\lesssim b$ to mean that there exists a constant $C>0$, independent of $\epsilon$ and whose specific value is unimportant, such that $a\leq Cb$.

[TABLE]

where $\overline{D}:=\sup_{v\in[0,\,1]}D(v)$ and moreover

[TABLE]

with

[TABLE]

Plugging (34) into (33) and considering that $1-\chi(0\leq v^{\prime}\leq 1)\leq 1$ and that, for $\epsilon$ sufficiently small, it results $\alpha\lesssim\epsilon$, we compute:

[TABLE]

The last six terms on the right-hand side are bounded and tend to zero as $\epsilon\to 0^{+}$. Conversely, the first three terms need to be estimated more carefully.

From the interactions (31) scaled according to (13)-(23), we see that a sufficient condition for $v^{\prime}\geq 0$ is

[TABLE]

Let us take $\epsilon<\frac{\kappa}{\kappa+1}$, so that $\epsilon\left(1+\frac{1}{\kappa}\right)-1<0$. Then, invoking the assumptions (6), (21) on $D$, a further sufficient condition for $v^{\prime}\geq 0$ is

[TABLE]

If we restrict $v$ to any compact subset $U$ of $(0,\,1)$ then there exists a constant $\underline{D}>0$ such that $D(v)\geq\underline{D}$ for all $v\in U$. Consequently, the right-hand side of (36) is in $U$ a continuous function of $\epsilon$, which tends to $-\infty$ when $\epsilon\to 0^{+}$. Hence, it is negative for $\epsilon$ small enough, say $\epsilon<\epsilon_{0}$. On the whole, we may define

[TABLE]

and state that, for $v\in U$, a sufficient condition for $v^{\prime}\geq 0$ is $Y\geq-b_{\epsilon}(v)$ provided $\epsilon<\min\{\frac{\kappa}{\kappa+1},\,\epsilon_{0}\}$. Arguing analogously, we deduce that, for $v\in U$, a sufficient condition for $v^{\prime}\leq 1$ is $Y\leq b_{\epsilon}(v)$ with the same restriction on $\epsilon$ as above, so that finally

[TABLE]

It follows

[TABLE]

whence, from (35) with $U:=\operatorname{supp}{\varphi^{\prime}}\cup\operatorname{supp}{\varphi^{\prime\prime}}\subset\subset(0,\,1)$,

[TABLE]

Let now $p,\,q\in[1,\,+\infty]$ be such that $\frac{1}{p}+\frac{1}{q}=1$. Hölder’s and Chebyshev’s inequalities imply

[TABLE]

where the last one holds for $\epsilon<\min\{\frac{\kappa}{\kappa+1},\,\epsilon_{0}\}$. Therefore, we may deduce:

[TABLE]

and finally

[TABLE]

which implies $R_{\epsilon}(f_{\epsilon},\,f_{\epsilon})[\varphi]\to 0$ for $\epsilon\to 0^{+}$. Notice, in particular, that $\langle\lvert Y\rvert^{\frac{5}{2}}\rangle^{\frac{2}{5}}\leq\langle\lvert Y\rvert^{3}\rangle^{\frac{1}{3}}<+\infty$ by Hölder’s inequality again.

Consequently, as already anticipated, the large time aggregate effect of the deterministic control (30) is actually the same as that of the uncertain control (19). Indeed, thanks to the result just obtained, the quasi-invariant limit for $\varphi(v)=v$ yields equation (24), whence we deduce that the asymptotic mean speed is again (25) and that all the conclusions drawn in Section 3.1, in particular (26), (27), (29) and Theorem 3.1, hold true also in this case.

In Figure 8, we exemplify the validity of the theoretical result just obtained. We show the numerical approximation of the large time average solution $\bar{f}_{\infty}(v)$ to the Boltzmann-type equation (32), computed by means of a standard Monte Carlo scheme with rejection in the variable $v$ and a collocation method for uncertainty quantification in the variable $z$ in the case $z\sim\mathscr{U}([1,\,3])$, cf. [31, 32] and [13, 42]. By comparing, for decreasing $\epsilon$, such a solution with the $\bar{f}_{\infty}(v)$ obtained from (29) by means of the Gauss-Legendre quadrature formula in $z$, we observe the expected consistency of the Boltzmann-type model (32) with the Fokker-Planck model (28) in the quasi-invariant limit.

4 Numerical tests

In this section, we present several numerical tests, which confirm the theoretical findings of the previous sections. First, we focus on the Fokker-Planck-type modelling derived in Section 2.3 in the case of uncontrolled microscopic interactions. Next, we show the effectiveness of the control approaches proposed in Sections 3.1, 3.2.

4.1 Generalised polynomial chaos for uncertainty quantification

Numerical methods for Uncertainty Quantification (UQ) gained a lot of momentum in recent times, thanks to their spectral convergence on the random field obtained under suitable regularity assumptions [47, 48]. The development of UQ methods for kinetic equations encoding uncertain quantities presents several difficulties, which are essentially due to the increased dimensionality and the intrinsic structural properties of the solution. Recently, a lot of efforts have been devoted to the study of kinetic equations with uncertainties, see e.g. [13, 19, 50]. Despite their spectral accuracy, Stochastic Galerking-generalised Polynomial Chaos (SG-gPC) methods typically lead to the loss of physical properties of the solution, such as the non-negativity of the distribution function, the entropy dissipation and the conservation of the steady states. We refer the interested reader to [7, 6, 49] for further insights into the preservation of the structural properties of the solutions to kinetic equations for the collective behaviour of multi-agent systems.

For completeness, in the following we report some fundamental notions on SG-gPC methods, in order to introduce the numerical scheme which we use for our numerical tests.

Let $(\Omega,\,\mathscr{F},\,\operatorname{Prob})$ be a probability space, where $\Omega$ is an abstract sample space, $\mathscr{F}$ a $\sigma$-algebra of subsets of $\Omega$ and $\operatorname{Prob}$ a probability measure on $\mathscr{F}$. Let us introduce, furthermore, a random variable

[TABLE]

where $I_{Z}\subseteq\mathbb{R}$ is the range of $z$ and $\mathscr{B}(\mathbb{R})$ is the Borel $\sigma$-algebra of subsets of $\mathbb{R}$. We assume that the distribution of $z$ is expressed by a known probability density function $\Psi:I_{Z}\to\mathbb{R}_{+}$.

Given $M\in\mathbb{N}$, the SG-gPC method is based on introducing the linear space of polynomials of degree at most $M$, generated by an orthonormal set $\{\Phi_{k}\}_{k=0}^{M}$ satisfying

[TABLE]

where

[TABLE]

is the Kronecker delta. The choice of the $\Phi_{k}$’s depends clearly on the density function $\Psi$. This aspect has been extensively investigated in the UQ literature in relation with the so-called Wiener-Askey scheme for orthogonal polynomials, see [47, 48]. In the following, we present some examples where the law of $z$ and the related set of orthonormal polynomials are among those reported in Table 1.

Assume that the kinetic distribution function $f$ solving either Fokker-Planck equation (15), (28) is sufficiently smooth, in particular $f(t,\,v;\,\cdot)\in L^{2}(I_{Z})$. Then, we can introduce the expansion

[TABLE]

where \hat{f}_{k}(t,\,v):=\mathbb{E}\bigl{(}f(t,\,v;\,z)\Phi_{k}(z)\bigr{)} for every $k=0,\,\dots,\,M$. Plugging the approximation $f^{M}$ of $f$ into the Fokker-Planck equation (15), we obtain

[TABLE]

Representing $f^{M}$ by means of the expansion above, multiplying the equation by $\Phi_{h}$ and using the orthonormality of the set $\{\Phi_{k}\}_{k=0}^{M}$ yields the following coupled system of $M+1$ deterministic partial differential equations ruling the evolution of the projections $\{\hat{f}_{h}\}_{h=0}^{M}$ of $f$ onto the polynomial space:

[TABLE]

where

[TABLE]

By introducing the vector $\hat{\mathbf{f}}(t,\,v):=(\hat{f}_{0}(t,\,v),\,\dots,\,\hat{f}_{M}(t,\,v))$ and the $(M+1)\times(M+1)$-matrix $\mathbf{E}[\hat{\mathbf{f}}](t,\,v):=\{E_{hk}[f^{M}](t,\,v)\}_{h,k=0}^{M}$, whose entries are given by (39), problem (38) may be recast in vector form. In particular, we obtain the following deterministic system of $M+1$ Fokker-Planck equations:

[TABLE]

In the case of controlled microscopic interactions, we apply the same procedure to the Fokker-Planck equation (28) and, in place of (38), we get:

[TABLE]

where now

[TABLE]

Consequently, problem (41) can be recast in vector form as

[TABLE]

where $\mathbf{E}[\hat{\mathbf{f}}](t,\,v):=\{E_{hk}[f^{M}](t,\,v)\}_{h,k=0}^{M}$ is the $(M+1)\times(M+1)$ matrix whose entries are given by (42) and $\mathbf{I}$ is the $(M+1)\times(M+1)$ identity matrix.

In particular, for our purposes, it is of paramount importance to compute accurately the large time trend of the kinetic distribution function $f$. Aiming at preserving the steady state, we solve numerically (40), (43) by means of a micro-macro gPC scheme, see [13, 30], together with a second order semi-implicit time integration, see [5, 33].

4.2 The uncontrolled case

We begin by considering the Fokker-Planck model (40), i.e., in particular, the one in which microscopic interactions are not controlled, cf. (15), complemented with a deterministic initial condition independent of the traffic density $\rho$:

[TABLE]

Concerning the distribution of the uncertain parameter $z$, we explore mainly the following two cases:

1. (i)

uniform: $z\sim\mathscr{U}([1,\,3])$, corresponding to a Legendre polynomial chaos expansion; 2. (ii)

binomial: $z$ is such that $z-1\sim\operatorname{B}\!\left(50,\,\frac{1}{50}\right)$, corresponding to a Krawtchouk polynomial chaos expansion,

cf. Table 1. Notice that, in both cases, the mean value of the exponent of the probability of acceleration (3) is $\mathbb{E}_{z}(z)=2$.

In Figure 9, we show the trend of the $L^{2}$-error produced by the numerical scheme with respect to the analytical steady state $f_{\infty}(v;\,z)$ given by (17) for an increasing number $M$ of modes of the gPC expansion (37). For this convergence test, we considered a uniform discretisation of the variable $v\in[0,\,1]$ by $41$ gridpoints and a uniform discretisation of the time interval $[0,\,T]$, with final time $T=60$, by a time step $\Delta{t}=1$. Moreover, we fixed the traffic density to $\rho=0.4$. From Figure 9(a), which refers to the uniformly distributed $z$ with Legendre polynomial chaos expansion, we observe that we reach essentially the machine precision with a relatively small number of modes. The same spectral convergence to zero of the error may be appreciated in Figure 9(b), which illustrates the case of the binomial distribution with Krawtchouck polynomial chaos expansion and, for completeness, also the case in which $z$ has beta distribution on the interval $I_{Z}=[1,\,3]$ with zero mean and variance equal to $\frac{1}{3}$. Notice that this case corresponds actually to $z\sim\mathscr{U}([1,\,3])$, using however the Jacobi polynomial chaos expansion.

Figures 10, 11 show the contours of the mean $\bar{f}(t,\,v)$ and of the variance $\operatorname{Var}_{z}(f(t,\,v;\,z))$ of the solution $f$ to the Fokker-Planck equation (15) for $\rho=0.2,\,0.4,\,0.6$, computed by means of the polynomial chaos expansion as

[TABLE]

We used $M=20$ modes and the same discretisations of $v$ and $t$ described above. It is evident that the distribution $\Psi$ of $z$ strongly affects the results. In particular, we observe substantial differences in the instantaneous variance, which highlights the regions where the solution to (15) is more sensitive to the uncertainty introduced by $z$.

4.3 The controlled case: uncertainty damping by control methods

We consider now the Fokker-Planck model (41), in which the microscopic interactions (20) are controlled by the ADAS technology through the uncertain control (19). Owing to the results of Section 3.2.1, we know that the observed large time trends are also produced by the microscopic interactions (31) implementing the deterministic control (30).

In particular, we prescribe the same deterministic initial distribution (44) as in the uncontrolled case and, as far as the distribution of the uncertainty is concerned, we explore the same cases $z\sim\mathscr{U}([1,\,3])$ and $z-1\sim\operatorname{B}\!\left(50,\,\frac{1}{50}\right)$ already considered in the uncontrolled case. Moreover, we fix the penetration rate to $p=0.1$, meaning that $10\%$ of the vehicles in the traffic stream are equipped with the ADAS technology, and we consider various penalisation coefficients $\kappa>0$ of the control. We remind that, in (28), these two factors are linked through the parameter $p^{\ast}$. Finally, we set the optimal speed to $v_{d}(\rho)=1-\rho$, like in the cases illustrated in Figures 6, 7.

With the aim of investigating the ability and the effectiveness of the control method to dampen the structural uncertainty of the model, in Figures 12, 13 we show the variance $\operatorname{Var}_{z}(f(t,\,v;\,z))$ of the solution $f$ to the Fokker-Planck equation (28) for $\rho=0.2,\,0.4,\,0.6$ and for a decreasing penalisation $\kappa=10^{-1},\,10^{-2}$ of the control, which implies an effective penetration rate increasing from $p^{\ast}=1$ to $p^{\ast}=10$. We observe that the control produces in time a confinement of the $z$-variance of $g$ around the optimal speed $v_{d}(\rho)$. In particular, such a confinement is tighter and is reached more quickly with a high effective penetration rate $p^{\ast}$.

To further stress this fact, in Figure 14 we show a numerical approximation of the asymptotic variance $\operatorname{Var}_{z}(f_{\infty}(v;\,z))$, where $f_{\infty}(v;\,z)$ is given by (29), computed at a sufficiently large time ensuring that the equilibrium has been numerically reached. We fix $\rho=0.4$ and we consider a decreasing penalisation coefficient $\kappa=10^{5},\,10^{-1},\,10^{-2}$. We see that the action of the control confines progressively the variability of $g_{\infty}$ around the optimal speed $v=v_{d}(\rho)$, namely around $v=0.6$ in the present case. At the macroscopic level, this produces the uncertainty damping observed in the fundamental diagrams, cf. Figure 6.

5 Conclusions

In this paper, we have shown that the scattering of the fundamental diagram of traffic, usually observed from experimental data at high vehicle density, can be explained organically in terms of the propagation of some uncertainty from the microscopic scale of the interactions among the vehicles to the macroscopic scale of the hydrodynamical flow of the vehicles. Such an uncertainty may be associated with different types of vehicles present in a heterogeneous traffic stream. We have shown that the classical methods of the collisional kinetic theory, suitably coupled with arguments typical of the uncertainty quantification, provide a sound mathematical-physical background to investigate such a propagation across the scales and to design possible countermeasures to dampen it. In particular, we have suggested that these countermeasures may take advantage of automatic driver-assistance systems, namely of feedback controls implemented on a certain percentage of vehicles in the traffic stream. By means of the theoretical paradigm recalled above, we have studied the large scale aggregate impact of some simple microscopic control protocols. Specifically, we have proved that it is actually possible to reduce the macroscopic scattering of the fundamental diagram by inducing some vehicles in the stream to align locally their speed to a congestion-dependent recommended speed. We have proposed that this can be done in two alternative ways, namely by either a stochastic or a deterministic control, and we have proved that, in the long run, these two strategies are actually equivalent. The analysis of the two strategies has required to deal with two different types of kinetic equations: a simpler Boltzmann-type equation for Maxwellian-like particles in the case of the stochastic control; a more difficult Boltzmann-type equation with cut-off in the case of the deterministic control, which imposes a non-constant interaction kernel. In this respect, this paper gives also a contribution to the investigation of this second case, which is less standard, hence to date still under-explored, in the literature of the kinetic models of multi-agent systems.

The possibility to take advantage of the emerging technologies in the realm of vehicle automation to control the fundamental diagram, and in particular to reduce its scattering, is, in our opinion, a non-trivial fact with potential beneficial impacts on the governance of traffic at large scale. Furthermore, from the point of view of the mathematical research, we believe that this work paves the way for the study of macroscopic models of traffic flow with uncertainty and uncertainty control through driver-assist vehicles, cf. [11, 12], grounded on the multiscale framework of the kinetic theory.

Acknowledgements

This research was partially supported by the Italian Ministry for Education, University and Research (MIUR) through the “Dipartimenti di Eccellenza” Programme (2018-2022), Department of Mathematical Sciences “G. L. Lagrange”, Politecnico di Torino (CUP: E11G18000350001) and Department of Mathematics “F. Casorati”, University of Pavia; and through the PRIN 2017 project (No. 2017KKJP4X) “Innovative numerical methods for evolutionary partial differential equations and applications”.

This work is also part of the activities of the Starting Grant “Attracting Excellent Professors” funded by “Compagnia di San Paolo” (Torino) and promoted by Politecnico di Torino.

Both authors are members of GNFM (Gruppo Nazionale per la Fisica Matematica) of INdAM (Istituto Nazionale di Alta Matematica), Italy.