Estimate of traffic emissions through multiscale second order models with heterogeneous data

1.   Introduction

In this paper, we focus on the development of models specifically designed to take advantage of the availability of heterogeneous data. By heterogeneous data we mean not only data coming from different sources, but especially data coming from different scales of observation and different modes of monitoring. We refer in particular to Lagrangian data, which provide information on the trajectories followed by vehicles, and to Eulerian data, which measure the transit of cars from fixed locations. In the case of vehicles, the Lagrangian data are typically GPS data (i.e. trajectory data, from which the instantaneous speed is derived), while the Eulerian data come from fixed sensors placed along the road, capable of counting cars and measuring their speed.

Alongside this analysis we consider the problem of estimating emissions from vehicular traffic on complex networks. The continuous traffic growth is in fact associated with negative environmental effects, which are related to both air quality and climate change. In order to assess the impact of traffic emissions on the environment and human health, an accurate estimate of their rates and location is required.

1.1.   Related work

In this article we mainly follow the approach of [11], where the authors propose a new traffic model that integrates position and velocity information of a single tracked vehicle into the velocity function of the Lighthill-Whitham-Richards (LWR) model [25,28]. LWR is a first order (i.e. a single equation) model that describes traffic dynamics in a road through the density of vehicles and their average speed. We extend the ideas of [11] to the case of several vehicle trajectory data and to the family of Generic Second Order Models (GSOM) introduced in [23]. GSOM encompasses the majority of second order (i.e. two equations) traffic models that, unlike first order ones, are able to reproduce bounded traffic accelerations [22]. Our choice leads to a multiscale type model that exploits the best features of the macroscopic and microscopic approaches. Multiscale models has been already considered in several papers. We refer, for instance, to [10,16] where the macroscopic LWR model is merged with the classical microscopic follow the leader model. In [12] the authors propose a multiscale approach obtained by coupling a first order macroscopic model with a second order microscopic one that is used only under specific traffic conditions. The interested reader can find other examples of multiscale models in [7,15,18,20,24]. Hereafter we cite some works closer to our goals. In [9], the LWR model is combined with an ordinary differential equation representing the trajectory of a slow vehicle acting as a moving bottleneck. This approach can be considered as a way to include real trajectory data in a macroscopic traffic model and is therefore comparable with our scopes. We also refer to [8] and the references therein, to mention a robust method of traffic estimation involving mixed fixed and mobile sensor data using the Hamilton-Jacobi equations.

Once the traffic state variables have been estimated, they can be used as input for the so-called emission models, which evaluate the mass of pollutants emitted. In this respect, in [19,36] the authors provide general frameworks to integrate macroscopic traffic flow models and microscopic emission models. In [2,3] the traffic modeling relies on the LWR model and a reaction-diffusion model describes the spread of carbon monoxide in the air with a source term associated with traffic dynamics. In the recent work [17] the authors analyze a reaction-diffusion model, based on LWR traffic dynamics, to control nitrogen oxides emissions. In [30] the authors suggest a new methodology to estimate in real-time the emission rates of pollutants and describe their diffusion in air. Furthermore, in [4] the authors propose a computational tool to estimate pollutant emissions due to vehicular traffic using second order traffic models. This approach approximates emissions well when a large amount of data is available to feed the traffic model. The present work extends this method by including microscopic data in the second order macroscopic traffic model, and provides a good estimate of pollutant emissions even when few data is available.

1.2.   Background and motivations

In order to exploit as much information as possible, we consider traffic models that take into account different and heterogeneous traffic data available along a road. A single source of data is generally not sufficient to estimate and forecast traffic volumes on a road. In Figure 1 on the left we provide an example of the partial information coming from GPS data; speed is observed for a limited number of vehicles and only at specific points in time and space. At the same time, in Figure 1 on the right we show an example of flux data coming from a fixed sensor; this type of data is not enough to correctly calculate traffic quantities such as densities, as it only provides average speed values or noisy, time-sampled flow information.

To deal with these two types of data, we have followed the approach proposed in [11], which is an effective and efficient way of coupling macroscopic and microscopic variables to estimate traffic volumes. This approach allows us to perform data fusion directly at the model level: Eulerian flow data measured by stationary sensors is used as boundary condition for the differential equations describing the traffic dynamics, while the Lagrangian data from the GPS sensor is used to correct macroscopic quantities in real time. More precisely, let $p = p(t)$ be the position of a single vehicle along the road at time $t$. Starting from the conservation law of vehicles density $\rho = \rho(x,t)$,

the measured trajectory $p = p(t)$ is directly encoded in the time- and space-dependent speed law ${\mathcal{U}}$ as

if $\dot{p}(t)+u(\rho)\neq0$, otherwise ${\mathcal{U}} = 0$. Here, $u = u(\rho)$ is a given speed function depending only on the density and $\chi(\xi)$ is a smooth, non-negative function such that $\chi(\xi) = 1$ for $|\xi|\leq \ell$ and $\chi(\xi) = 0$ when $|\xi|\geq L$, for two fixed constant $\ell,L$ with $\ell<L$. Thus, when the point $(x,t)$ is sufficiently far from the position of the vehicle $p = p(t)$ ($\chi = 0$), then ${\mathcal{U}}$ is defined by the velocity function $u(\rho)$ as in the LWR model, otherwise it is given by the harmonic mean between $u$ and the velocity $\dot p(t)$ of the vehicle.

We extend this idea to the case of more trajectory data available. We propose two ways to incorporate the Lagrangian data of $N$ vehicles. In the first one, the velocity function ${\mathcal{U}}$ takes into account only the speed $\dot p_k(t)$ of the $\kappa$-vehicle closest to the point $x$ at time $t$; in the second one, it depends on the mean speed of the set of vehicles closest to point $x$. The difference between the two approaches only shows up when several vehicles with different speeds are close to each other. In both cases, the resulting flow function ${\mathcal{F}} = \rho{\mathcal{U}}$ is strictly concave with respect to $\rho$. We do not investigate the model from an analytical point of view. We need ${\mathcal{U}}$ to inherit reasonable properties of $u$ and $\dot p$, such as positivity and boundedness, and we still refer to [11] for analytical aspects.

In absence of real trajectory data, the macroscopic traffic model introduced above can be coupled with a microscopic one. Lagrangian data, in fact, can be generated from a microscopic model and then included in the velocity function ${\mathcal{U}}$ as previously described. One therefore obtains a multiscale type model capable of describing some traffic phenomena generally captured only by a detailed, microscopic description of the dynamics.

The proposed approaches can be easily applied to other traffic models. We can thus incorporate the information from $N$ Lagrangian data into the family of Generic Second Order Model (GSOM) [23] described by

where $\rho(x,t)$ is the density of vehicles, $w(x,t)$ is a vehicle/driver property or invariant, which is conserved along trajectories, and ${\mathcal{V}}$ is the velocity field. By taking into account the Lagrangian information, ${\mathcal{V}} = {\mathcal{V}}(x,t,\rho,w; \textbf{p})$ is defined as

if $\dot{p}_{\kappa}(t)+v(\rho,w)\neq0$, otherwise ${\mathcal{V}} = 0$. Here $\textbf{p} = (p_1(t),\ldots,p_N(t))$ is the vector of positions of the $N$ vehicles and $\kappa = \kappa(x,t)$ is the index of the closest vehicle to point $x$ at time $t$. The function $v = v(\rho,w)$ is a given analytical function which describes the macroscopic velocity field in the native (without tracking vehicles) second order model.

1.3.   Emissions estimate

The use of a second order model leads to good approximations of the acceleration of vehicles and consequently to the estimate of traffic emissions at ground level, that is one of the main goal of our work. Most emission models are based on both vehicle speed and acceleration, see for instance [6,32] and references therein. Here we explore the use of two types of formula which compute the emissions associated with the motion of vehicles. The models we consider have been introduced in [26] and [1], respectively. In both formulas the contribution of the Lagrangian data is incorporated in the terms of speed and acceleration evaluated by the traffic model. Specifically, by computing the acceleration function $a = a(x,t)$ as the time material derivative of the new function ${\mathcal{V}}$ introduced above it directly depends on the acceleration $\ddot{p}_{\kappa}(t)$ of the nearest vehicle.

We compare the two emission models and we show how the integration of real data affects their results. We observe that the integration of trajectory data into the macroscopic traffic model increases the order of accuracy of the emission estimate, even when there are few data available. In particular, the formula proposed in [26] gives better results.

We conclude our study with a real-life application using trajectory and fixed sensors data provided by Autovie Venete S.p.A. on the Italian A4 (Trieste-Venice) highway. With this test we link heterogeneous traffic source data to emission estimates along a road network, and at the same time we provide an approximation of the source term that feeds air pollutant diffusion and chemical reaction models. The numerical results show how the GPS data influences the solution of the traffic model and gives good reproductions of the emission peaks at the macroscopic scale.

1.4.   Main goal

In summary, we propose a second order traffic model that returns macroscopic traffic quantities by incorporating microscopic information. Microscopic trajectories are included in the definition of the velocity field in order to perturb the velocity and acceleration values at the macroscopic level. This methodology combines the computational efficiency of a macroscopic model with the accuracy of a microscopic representation. This makes it particularly suitable as an input for estimating the mass of emitted pollutants when an aggregate description is required. With a few Lagrangian trajectories, it is in fact possible to reproduce significant emission variations at the macroscopic scale. The procedure is very flexible and can be used with real measurements or with vehicle trajectories generated by Lagrangian models.

1.5.   Paper organization

In Section 2 we propose two possible extensions of the first order LWR model to integrate Lagrangian data from $N$ vehicles. A numerical test is proposed to show the differences between the two approaches. In Section 3, we apply the ideas given in the previous section to the family of GSOM. Then, we compute the acceleration as the material derivative of the velocity function, making explicit its dependence on the acceleration of the single vehicles. In Section 4 we describe two models to estimate the emissions produced by vehicular traffic. In Section 5.1 and 5.2, we propose numerical tests to show how the integration of Lagrangian data into the GSOM impacts on the traffic dynamic and on the estimate of emissions, respectively. We conclude with Section 5.3, which describes the use of GPS and fixed sensors data provided by Autovie Venete S.p.A. on a portion of the Italian highway network. Finally, in Appendix A we report two technical proofs.

2.   A macroscopic first order model embedding Lagrangian data

Assume to know the trajectory of $N$ vehicles, and let $\textbf{p} = (p_1(t),\ldots,p_N(t))$ be the vector of their positions at time $t$. Following [11], we assume $p_i = p_i(t)$, to be continuous and smooth functions such that $\dot{p}_i\geq 0$ for $i = 1,\cdots,N$ and a.e. $t\in{\mathbb{R}}^+$. We propose two approaches to incorporate information from the vehicles into the velocity function of the conservation law (1).

(CV) Closest Vehicle. We define ${\mathcal{U}}$ in (1) as

where $p_\kappa(t)$ is the position of the $\kappa$-vehicle closest to point $x$. With a slight abuse of notation, we have

(ACVs) Average on the Closest Vehicles. For each vehicle position $p_i(t)$, $i = 1,\ldots,N$, we set

We introduce the function $\phi(x,t)$ that counts the number of vehicles whose position at time $t$ is close to point $x$, i.e. such that $\chi(x-p_i(t))>0$,

We then define the function ${\mathcal{U}}$ in (1) as the average value of the speeds of the closest vehicles, that is

The two approaches are different only when two or more trajectories are very close to each other, otherwise they coincide.

We use the same assumptions made in [11], namely: $\rho\in[0,{\rho^{\mathrm{max}}}]$, $u = u(\rho)$ is a regular and non increasing function with $u({\rho^{\mathrm{max}}}) = 0$, $f(\rho) = \rho u(\rho)$ is strictly concave in $\rho$ and, for each $z>0$ the function ${\rho z u}/{(z+u)}$ is strictly concave in $\rho$; the flux function ${\mathcal{F}}(x,t,\rho) = \rho{\mathcal{U}}(x,t,\rho; \textbf{p})$, with ${\mathcal{U}}$ in (2) or (5), is strictly concave with respect to $\rho$. Thus, there exists a unique density value $\sigma(x,t)$ where the flux reaches its maximum ${\mathcal{F}}(x,t,\sigma) = {\mathcal{F}}^{\max}(x,t)$ and we can define the sending $S = S(x,t,\rho)$ and receiving $R = R(x,t,\rho)$ functions in the standard way [21],

The critical density $\sigma(x,t)$ is defined by ${\partial}_{\rho}{\mathcal{F}}(x,t,\sigma) = 0$. In the (CV) case, for $(x,t)$ such that $\chi(x-p(t))\neq 0$ or $\chi(x-p(t))\neq 1$, it is implicitly defined by the non-linear relation

Also in the (ACVs) case the computation of $\sigma$ requires the solution of a non-linear problem,

In both cases, for each point $(x,t)$, the calculation of $\sigma = \sigma(x,t)$ must be done numerically.

Remark 2.1. Numerical solvers for nonlinear equations such as (7) or (8) have a high computational cost. Since the critical density $\sigma$ is the maximum point of the flow function ${\mathcal{F}}$, one can reduce the complexity of the computation by sampling the values of ${\mathcal{F}} = {\mathcal{F}}(\cdot,\cdot,\rho)$ into a vector and applying a search for the maximum value of the vector that returns the corresponding density index.

To highlight the differences between the two proposed models (CV) and (ACVs), in Figure 2 we plot different curves ${\mathcal{F}}(x,t,\rho)$ with one or more trajectory data. We set $u(\rho) = {u^{\mathrm{max}}}({\rho^{\mathrm{max}}}-\rho)/{\rho^{\mathrm{max}}}$ and

In Figure 2 on the left we fix the time $t$ and we plot the curves ${\mathcal{F}}(x,t,\rho)$ in the case of a single trajectory data. In this case the two approaches (CV) and (ACVs) coincide, and we observe that the critical density $\sigma(x,t)$ increases with $\chi(x-p(t))$ while ${\mathcal{F}}^{\max}(x,t)$ decreases. In Figure 2 on the right, we compare the flux function ${\mathcal{F}}(x,t,\rho)$ given by the (CV) and (ACVs) approaches assuming to have three vehicles travelling with different speed and close enough to influence the velocity function (5). The solid line represents the flux function ${\mathcal{F}}(x,t,\rho)$ with the (ACVs) approach on the cell $x$ where all the three vehicles contribute to the velocity computation, i.e. $\phi(x,t) = 3$. The dotted lines, instead, represent the flux function with the (CV) model in the cells $x$ where $\chi(x-p_{\kappa}(t)) = 1$, $\kappa\in\{1,2,3\}$. We observe that $\sigma(x,t)$ increases with the decrease of vehicle velocity. Indeed, for small velocities the vehicles need more time to fill the road, thus we obtain higher values of the critical density.

Let us consider now a numerical grid on a road $[a,b]$. Let $N_{x}$ be the number of cells $[x_{j-1/2},x_{j+1/2})$ of size ${\Delta x}$, and $N_t+1$ the number of intervals of length ${\Delta t}$ into which we divide the period of time $[0,T]$. Let $\rho^{n}_{j} = \rho(x_{j},t^{n})$ be the density of vehicles into the cell $x_{j}$ at time $t^{n}$, defined as the cell average

To approximate the model (1) we use the Cell Transmission Model (CTM) [13]. The numerical scheme has the following structure

with the numerical flux $F^{n}_{j+1/2}$ defined as

where $S$ and $R$ are the sending and receiving functions introduced in (6). At each time step $n$ and for each cell centered in $x_j$, the critical densities $\sigma(x_j,t^n)$ must be computed numerically. As observed in Remark 2.1, the computational costs are reduced by calculating the values ${\mathcal{F}}^n_j(\rho) = {\mathcal{F}}(x_j,t^n,\rho)$ as $\rho\in[0,{\rho^{\mathrm{max}}}]$ varies and applying a search for the maximum of these values to compute the index $l$ such that $\sigma(x_j,t^n)\approx \rho_l$. For instance, on the MATLAB platform the ${\texttt{max}}$ function allows to efficiently compute the index of the maximum of the sampled flux values and is significantly faster than the function ${\texttt{fsolve}}$ required to solve non-linear equations as (7) or (8).

With standard arguments it is possible to prove the following properties of the scheme (10)-(11). Both results given in the next Proposition 2.2 and 2.3 apply to the (CV) and (ACVs) models. For completeness, in Appendix A we give the details of the proofs for the (CV) approach.

Proposition 2.2. The approximate solution $\boldsymbol\rho^n = (\ldots, \rho^n_{j-1},\rho^n_j,\rho^n_{j+1},\ldots)$ constructed via the scheme $(10)-(11)$ satisfies the bounds

if the following Courant-Friedrichs-Levy (CFL) condition holds:

Proposition 2.3. The stability of the scheme is guaranteed by the condition

where $\dot{p}_\kappa$ is the velocity of the $\kappa$-th vehicle, with $\kappa$ defined in $(3)$, and ${u^{\mathrm{max}}}$ is the maximum velocity value given by the analytical model.

Remark 2.4. If all the trajectories on the road have a speed always lower than the parameter ${u^{\mathrm{max}}} = u(0)$ of the analytical model, then we use the standard CFL condition ${\Delta t}\leq {\Delta x}/{u^{\mathrm{max}}}$.

Now we propose an example to show the differences between the two approaches (CV) and (ACVs). Recall that $u(\rho) = {u^{\mathrm{max}}}({\rho^{\mathrm{max}}}-\rho)/{\rho^{\mathrm{max}}}$. We fix the following parameters: ${\rho^{\mathrm{max}}} = 100\,{\mathrm{veh/km}}$, ${u^{\mathrm{max}}} = 90\,{\mathrm{km/h}}$, $a = 0$, $b = 3\,{\mathrm{km}}$, ${\Delta x} = 0.1\,{\mathrm{km}}$, $T = 1\,{\mathrm{min}}$. The time step ${\Delta t} = 0.2\,{\mathrm{s}}$ is chosen to satisfy the CFL condition in (13). The function $\chi$ in (2) is defined in (9) with $\ell = 2{\Delta x}$ and $L = 6{\Delta x}$. We test two different initial data, namely a rarefaction wave and a shock one, given by

Since the two definitions of ${\mathcal{U}}$ in (2) and (5) differ for $\phi(x,t)>1$ with $\phi$ in (4), we fix the trajectories of three vehicles as

so that $\phi(x_{k},t)>1$ for some cells $x_{k}$ and time $t$. In Figure 3 we compare the results at different times starting from the rarefaction wave on the top and the shock one on the bottom. Let us consider the rarefaction case: the main differences between the solutions to (1) with ${\mathcal{U}}$ in (2) and (5) can be observed in plots (A), (B) and (C), i.e. when the three vehicles are really close to each other ($\phi(x,t)>1$). Indeed, the solutions related to ${\mathcal{U}}$ in (2) (red-solid lines) show oscillations with higher picks than those related to (5) (blue-dotted lines), where the velocity value is averaged. At the final time of the simulation, plot (D), the two solutions are quite similar, since the three vehicles are far enough from each other ($\phi(x,t)\leq1$ for each $x$) and the two definitions of ${\mathcal{U}}$ in (2) and (5) coincide. A similar analysis holds for the shock wave initial data. Therefore, the main difference between the two dynamics (red-solid and blue-dotted lines) concerns the height of the peaks of the oscillations caused by the individual behavior of the tracked vehicles.

We conclude this section by highlighting another important feature of the model described above. The proposed approach also allows the coupling of a macroscopic model with a microscopic one. Put simply, trajectory data can be generated from a microscopic model, and then included in the velocity function ${\mathcal{U}}$ following the (CV) or (ACVs) methodology. By coupling, for example, the first order macroscopic model (1) with a second order microscopic one of the follow-the-leader type, the resulting multiscale model would have the following form

where the $N$-th vehicle is the leader which has a special dynamics, and the function $A$ represents the acceleration to be defined, see for instance [38].

We propose a numerical test, to show how the model (14) can reproduce typical traffic phenomena. In Figure 4 on the left, we generate trajectories by the second order microscopic model used in [12] (equations (1)-(6)-(7)), which is specifically designed to reproduce stop & go waves. Starting from a constant initial density along a road, in Figure 4 on the right we observe how the microscopic dynamics cause a variation on the macroscopic density, leading to the reproduction of the stop & go phenomenon at the macroscopic level. In this simulation, we have used the (CV) approach. For completeness, in Figure 5 we superpose the density profile obtained with the (CV) method to the (ACVs) one. The two profiles are very similar because the not-perturbed vehicles generated by the microscopic dynamics are all moving at the same speed.

The advantage of using a model of the type proposed in (14) is clear from Figure 6, where a small sample of microscopic data is sufficient to describe a stop & go wave that is difficult to simulate through a macroscopic model (particularly of the first order).

3.   A macroscopic second order model embedding Lagrangian data

In this section we apply the ideas introduced in Section 2 to the Generic Second Order Models (GSOM) [23], described by

where $\rho(x,t)$ is the density of vehicles and $w(x,t)$ is a property of vehicles advected by the velocity function ${\mathcal{V}}$. Similarly to the first order (CV) case, given $N$ vehicles with known trajectory $\textbf{p} = (p_1(t),\ldots,p_N(t))$, the velocity function ${\mathcal{V}} = {\mathcal{V}}(x,t,\rho,w; \textbf{p})$ is defined as previously given in the introduction by

if $\dot{p}_{\kappa}(t)+v(\rho,w)\neq0$ and ${\mathcal{V}}(x,t,\rho,w; \textbf{p}) = 0$ if $\dot{p}_{\kappa}(t)+v(\rho,w) = 0$, with $\kappa = \kappa(x,t)$ given in (3). Following [14], we assume that the flux function $Q(\rho,w) = \rho v(\rho,w)$ satisfies these properties: $Q(\rho,w)\in {C^{1}}$ and $Q({\rho^{\mathrm{max}}},w) = 0$ for each $w$; the flux is strictly concave with respect to $\rho$ and it is non-decreasing with respect to $w$ for each $\rho$. The flux function $Q$ defines the velocity function as $v(\rho,w) = Q(\rho,w)/\rho$. The assumptions on $Q$ imply that the function $v$ is regular and it is non-increasing with respect to $\rho$.

As a consequence of the properties of $Q$ and $v$, the flux function ${\mathcal{Q}}(x,t,\rho,w) = \rho{\mathcal{V}}(x,t,\rho,w; \textbf{p})$ with ${\mathcal{V}}$ in (16) is strictly concave with respect to the density $\rho$. Thus, as in the first order case, there exists a unique density value $\sigma = \sigma(x,t,w)$ where ${\mathcal{Q}}$ reaches its maximum ${\mathcal{Q}}^{\max}(x,t,w) = {\mathcal{Q}}(x,t,\sigma,w)$ with respect to $\rho$ and we can extend the definition of sending and receiving functions in (6) as

with $\sigma(x,t,w)$ critical density obtained by numerically solving the non-linear equation ${\partial}_{\rho}{\mathcal{Q}}(x,t,\sigma,w) = 0$, for each triple of values $(x,t,w)$, see Remark 2.1.

Under the same numerical setting introduced in Section 2, for $j = 0,\ldots,N_x-1$, $n = 0,\ldots,N_t$, we extend the scheme (10) to

with ${\mathcal{V}}^n_j = {\mathcal{V}}(x_j,t^n,\rho^n_j,w^n_j; \textbf{p})$ and

Since the variable $w$ is advected forward by the motion of vehicles travelling with positive velocity, the second equation has been approximated with an up-wind finite difference scheme.

With computations similar to the first order case given in Proposition 2.2 and 2.3, the following CFL condition

guarantees the stability of the scheme (17). Specifically, we get

where $\dot{p}_\kappa$ is the velocity of the $\kappa$-th vehicle, with $\kappa$ defined in (3), and ${V^{\mathrm{max}}}$ is the maximum velocity value given by the analytical model, i.e. ${V^{\mathrm{max}}} = \max_w v(0,w)$.

The finite volume formulation given in (17) allows for easy handling of flow and velocity data from fixed sensors. In order to exploit this information, the sensor data are used as boundary conditions on the incoming side of a road. Usually, sensor data are aggregated per minute and need to be interpolated to be available at each numerical time step $t^n$. Let ${\mathcal{Q}}^{n}_{{\mathrm{sens}}}$ and ${\mathcal{V}}^n_{{\mathrm{sens}}}$ be the sensor flux and speed data after interpolation, respectively. There are two main ways to include the real data into the numerical scheme (17): i) replace $Q^{n}_{-1/2} = {\mathcal{Q}}^{n}_{{\mathrm{sens}}}$ in the first cell $[x_{-1/2},x_{1/2}]$ making sure that the sensor datum is an admissible value with respect to the analytical bounds; ii) fix the density value $\rho^n_{0} = {\mathcal{Q}}^{n}_{{\mathrm{sens}}}/{\mathcal{V}}^n_{{\mathrm{sens}}}$ and use it in the calculation of $Q^{n}_{-1/2}$. In both cases, once the density at the left boundary has been computed, we estimate $w^{n}_{0}$ such that $v(\rho^{n}_{0},w^{n}_{0}) = {\mathcal{V}}^{n}_{{\mathrm{sens}}}$.

In the following, we employ the first approach i) that uses the sensor data in an almost pure form and avoids the additional approximation needed to calculate the density.

3.1.   Acceleration estimate

One of the main advantages of second order models is their greater accuracy in approximating velocity, which allows improvements in the estimate of vehicle acceleration.

Let us consider the model (15), with the velocity function ${\mathcal{V}}$ in (16). We compute the acceleration function $a = a(x,t)$ as the material derivative ($D/Dt = d/dt + v \cdot d/dx$) of ${\mathcal{V}} = {\mathcal{V}}(x,t,\rho(x,t),w(x,t); \textbf{p})$ with respect to the time $t$:

Therefore, from (15),

For $(\dot{p}(t),v(\rho,w))\neq(0,0)$ we have

where $\chi^\prime(\cdot)$ is the first derivative of $\chi$ with respect to its argument. Let us observe that the acceleration formula directly depends on the acceleration $\ddot{p}(t)$ of the vehicle that is influencing the speed function of the model.

4.   Emissions estimate from traffic quantities

Here we focus on the estimation of pollutants production at ground level due to vehicular traffic, whose impact on air quality is a long-standing and complex problem. Following [4], we propose a computational approach that combines the traffic simulation model with an emission one. We focus our attention on ${\mathrm{NO_{x}}}$ emissions, which are dangerous to human health and are precursors of the ozone, that also has negative effects on health [31,35].

Once approximated the density, velocity and acceleration of vehicles, we use them to estimate the emissions produced by vehicular traffic. Here we consider two microscopic emission models and we follow [4] to extend them to macroscopic variables.

Let us consider a vehicle $i$ moving at speed $v_{i}$ and subject to a certain acceleration $a_{i}$ at time $t$. The first microscopic emission model considered (see [26]), estimates the emissions associated with the motion of the vehicle as

where $f_{1}$–$f_{6}$ are coefficients that depend on the type of vehicle, fuel and pollutant considered. In Table 1 we collect the coefficients associated with ${\mathrm{NO_{x}}}$ emissions for a petrol car, and we refer to [26,Table 2] for an exhaustive list of coefficients related to other types of vehicles and pollutants.

The second emission model is based on [1]. For

the emissions associated with the vehicle $i$ are estimated by

where $P$ is the following $4\times4$ matrix (see [37,Appendix A1])

To extend these two models to the macroscopic case, we consider $M$ vehicles in a stretch of road moving at the same speed $\bar v(t)$ and subject to the same acceleration $\bar a(t)$. The emission rates associated with traffic in (20) are modified in

while from (21) we obtain

with

Hereafter we refer to the emission model (22) as the E-max-formula and to (23) as the E-exp-formula. In Figure 7 we show a comparison between these two formulations at a microscopic level, observing that in this example their results are quite similar. The details of this numerical test are described in Section 5.2.

5.   Numerical tests

This section is devoted to the numerical tests. First, we focus on the second order model (15), then we analyze emission models and finally, we propose an application on a road network representing a portion of the Italian A4 motorway, combining real GPS data and fixed sensors.

5.1.   Traffic dynamics

Let us begin with a test to illustrate the features of the second order traffic model (15)-(16). First of all, we choose the Collapsed-Generalized-Aw-Rascle-Zhang (CGARZ) model [14] among the family of GSOM. In the CGARZ model the definition of the flow function is characterized by the distinction between the free and congested flow traffic regime. Hence, we define the flux function as

for a given density threshold ${\rho_{f}}$ separating the two phases and

where, following [4], we set

for ${w_{L}}$ and ${w_{R}}$ minimum and maximum given value of $w$, respectively. With these choices, the functions $Q(\rho,w)$ and $v(\rho,w) = Q(\rho,w)/\rho$ satisfy the hypotheses required by model (15).

The numerical simulations are performed with the CTM scheme (17). We set ${\rho^{\mathrm{max}}} = 100\,{\mathrm{veh/km}}$, ${\rho_{f}} = 10\,{\mathrm{veh/km}}$, ${V^{\mathrm{max}}} = 90\,{\mathrm{km/h}}$, ${w_{L}} = g({\rho_{f}})$ and ${w_{R}} = g({\rho^{\mathrm{max}}}/2)$. The parameters $a$, $b$, ${\Delta x}$, $T$, ${\Delta t}$, the function $\chi$ and its parameters $\ell$ and $L$, are the same used for the test proposed in Section 2.

We consider three vehicles and simulate their trajectory as follows

We then fix the initial data

Figure 8 shows the final time of the simulation. We observe that the discontinuity in the variable $w$ generates a rarefaction wave on the right half of the road. The presence of the GPS trajectories, instead, produces non-classical shocks which are well captured by the CTM scheme.

5.2.   Application to emissions

In this section we propose a test to show how the integration of GPS data impacts the estimate of emissions due to vehicular traffic. To this end, we consider $N$ vehicles moving according to the following equations during the time interval $[0,T]$:

with $c = 0.3$, $k_{i} = 20+5(i-1)/(N-1)$ and $x_{0,i} = 1+0.05(i-1)$. With this choice, the vehicles are initially equally spaced every $50\,{\mathrm{m}}$ and then are subjected to different accelerations and decelerations, which we expect will affect the emissions. We fix $T = 20\,{\mathrm{min}}$ and the length of the road equal to $10\,{\mathrm{km}}$. Our analysis is then focused on a stretch of the road, specifically the interval $[4\,{\mathrm{km}}, 7\,{\mathrm{km}}]$, which is empty both at the beginning and at the end of the simulation, see Figure 9.

In order to compare microscopic quantities with macroscopic ones, we first derive the macroscopic density and velocity from vehicle trajectories through a kernel density estimation (KDE). We use the Parzen-Rosenblatt window method [27,29], which associates a density distribution with the position of the vehicles and then derives the global density by adding these distributions. More precisely, let $x_{i}(t)$ and $\nu_{i}(t)$ be the position and velocity of the $N$ vehicles, respectively. We define

where $\delta$ is the Dirac delta function, so that

In order to recover the smooth density and velocity $\rho$ and $\nu$, we introduce the Gaussian kernel

where $h$ is a smoothing parameters, which is chosen to obtain an almost constant density profile for equidistant vehicles. We then define

With this methodology we are able to reconstruct the initial density, $\rho_{0} = \rho(x,0)$, and velocity, $\nu_{0} = \nu(x,0)$, of the macroscopic model (15). Once $\rho_{0}$ and $\nu_{0}$ are known, we recover $w_{0}$ such that $v(\rho_{0},w_{0}) = \nu_{0}$.

The aim of this test is to compare the macroscopic emissions associated with the traffic model (15) with the microscopic ones given by vehicle trajectories. Indeed, we estimate the traffic quantities $\rho^{n}_{j}$ and $v^{n}_{j}$ by means of the CTM scheme introduced in Section 3, and the acceleration $a^{n}_{j}$ as described in Section 3.1. We then use these quantities in (22) and (23) to obtain the macroscopic emissions with the two proposed E-max-formula and E-exp-formula, denoted by $E^{n}_{j}$ and $\widehat E^{n}_{j}$, respectively. At the same time, the microscopic velocity and acceleration given in (25) are used in (20) and (21) to recover the microscopic emissions $e^{n}_{i}$ and $\widehat e^{n}_{i}$ for each vehicle $i$. We then compute the total amount of emissions along the road at time $t^{n}$ as

similarly for $\widehat E^{n}_{{\mathrm{tot}}}$ and $\widehat e^{n}_{{\mathrm{tot}}}$. Furthermore, we are interested in analyzing the effects on emissions of the number of tracked vehicles in model (15). To this end, let us introduce ${\widetilde N}$ as the number of vehicles which influence the macroscopic traffic model. We perform three simulations: in the first one ${\widetilde N}$ coincides with $N$, hence all the vehicles used for the microscopic dynamic influence the macroscopic model; then we consider only one vehicle out of two for model (15), hence ${\widetilde N} = N/2$, and finally one in four, ${\widetilde N} = N/4$.

For our simulations, we fix a $10\,{\mathrm{km}}$ long road and a time horizon of 20 minutes. The model parameters are the same of the first numerical test of Section 5.1. The bandwidth $h$ of the KDE is fixed as $h = 100 \,{\mathrm{m}}$. We use homogeneous Dirichlet boundary conditions on the left and homogeneous Neumann conditions on the right, which correspond to allow vehicles to leave the road. We then focus on the interval $[4\,{\mathrm{km}}, 7\,{\mathrm{km}}]$.

On the top plots of Figure 10 we compare the density of vehicles at $t = 10\,{\mathrm{min}}$ computed with model (15) and the one estimated through the KDE techniques. The black and fuchsia points represent vehicles position; in particular, the fuchsia points identify the vehicles which are tracked in model (15). From the plots we observe that the macroscopic model (15) (blue-solid line) produces some peaks in correspondence of vehicles position, but the profile of the KDE density (red-dotted line) is quite well-captured. The accuracy of the results increases with the number of monitored vehicles ${\widetilde N}$. The yellow line with circles represents the density associated with (15) without considering the trajectory data, obtained by setting ${\mathcal{V}} = v(\rho,w)$. We observe that the lack of GPS data produces a density profile completely different from the one recovered by microscopic data.

In the central (bottom) plots of Figure 10 we compare $E^{n}_{{\mathrm{tot}}}$ ($\widehat E^{n}_{{\mathrm{tot}}}$), with and without trajectory data, and $e^{n}_{{\mathrm{tot}}}$ ($\widehat e^{n}_{{\mathrm{tot}}}$) during the whole simulation. We observe that the trend and the absolute values of emission rates using the E-max-formula (22) are well captured by model (15), even when the number of tracked vehicles decreases. The total emission profile obtained from the model without GPS data is smoother than the other two since the dynamics do not take into account any vehicle acceleration or deceleration. On the other hand, the macroscopic E-exp-formula (23) well catches the behavior of emission rates but overestimates the absolute values.

Finally, in Table 2 we estimate the absolute error in $L^{1}$ norm at the final time of simulation, $\big\lVert{E^{{N_{t}}}_{{\mathrm{tot}}}-e^{{N_{t}}}_{{\mathrm{tot}}}}\big\rVert_{1}$ and $\big\lVert{\widehat E^{{N_{t}}}_{{\mathrm{tot}}}-\widehat e^{{N_{t}}}_{{\mathrm{tot}}}}\big\rVert_{1}$, with and without real trajectory data for ${\widetilde N} = N,N/2,N/4$. Note that the error associated with the model without real data is not influenced by ${\widetilde N}$, thus we report a unique value in the third and last column of the table. For the E-max-formula, the error slightly grows with the decrease of ${\widetilde N}$ and we gain an order of magnitude using the model with real data. For the E-exp-formula the error seems not to be affected by the number of tracked vehicles and is higher than the one obtained with the first emission formula.

To sum up, we can give a good approximation of emissions even when few real trajectory data are available. Moreover, the E-max-formula (22) allows for a better approximation of macroscopic emissions than the E-exp-formula (23) compared to the ground-truth emissions obtained using (20) and (21). Therefore, in the rest of the paper we will use only the E-max-formula.

Remark 5.1. The results of the tests described above are not in contrast with those proposed in [4,Section 3.1], where the CGARZ model without GPS data is used with the NGSIM dataset [34]. Indeed, the latter contains data for more than 5000 vehicles in 500 meters of road during 45 minutes of data recording. Hence, the large amount of real data allows quite accurate approximations of the emissions just using the CGARZ model with initial data and boundary conditions recovered from the real trajectories of vehicles.

5.3.   A real-life application with fixed sensors and GPS data

To conclude, we propose an application of traffic model (15) using real trajectory and fixed sensors data provided by Autovie Venete S.p.A. on the Italian A4 (Trieste-Venice) highway. The latter is a two lane motorway network travelled by light and heavy vehicles (cars and trucks). In Figure 11 we draw a sketch of the network with three diverge and three merge junctions connecting six roads. The triangles represent the fixed sensors which record the flux of vehicles that enter the network. Since we have access to more real truck trajectories than car ones, here we focus only on the dynamic of heavy vehicles and the goal is to apply the methodology described above to estimate the ${\mathrm{NO_{x}}}$ emission rates in a real-life scenario. The model can be used in the same way with real trajectory data of cars or other vehicle classes. As already discussed at the end of Section 5.2, we use the E-max-formula (22) for a more precise approximation of macroscopic emission rates. We use the coefficients $f_{1}-f_{6}$ in (22) calibrated for trucks, see [26,Table 2].

First of all we describe the main features of the network. We have three incoming roads (1, 3 and 5) and three outgoing roads (2, 4 and 6). The dynamic of vehicles along each road is described by the model (15). We use the notation $\rho_{r,j}^{n}$ to identify the density on road $r$ in cell $x_{j}$ at time $t^{n}$; similarly for the variable $w$. A special treatment is required at junctions, which are divided into two types: the diverge junctions, that divide a road in two, and the merge ones, which join two roads into one. Diverge junctions are ruled by a distribution coefficient that specifies the percentage of vehicles which prefer one road rather than the other one. For instance, junction $\mathrm{D1}$ has a certain distribution parameter $\alpha_{\mathrm{D1}}$ defining the percentage of vehicles which continue their path from road 1 to road 2, while $1-\alpha_{\mathrm{D1}}$ represents vehicles which continue on road 6. Similarly for the other diverge junctions. Following [5], merge junctions are characterized by a priority rule that establishes which of the two roads sends more flow of vehicles. We denote by $\beta_{\mathrm{M1}}$, $\beta_{\mathrm{M2}}$ and $\beta_{\mathrm{M3}}$ these parameters in $[0,1]$, and we calibrate them as a property of the network.

Now we describe the treatment of real data and we begin with data from GPS devices. In Figure 12 we show an example of real trajectory data related to heavy vehicles on the six roads of the network. This data provides information on the position and velocity of vehicles at various times. The time interval of the recorded positions is not constant and can vary from a few seconds to many minutes. In order to integrate the trajectory data into the numerical scheme, once a vehicle is located into a monitored road, we need to know its position with respect to the time step ${\Delta t}$. To this end, we linearly interpolate the given data to reconstruct missing information when necessary. Furthermore, this data is used to estimate the distribution parameters $\alpha_{\mathrm{D1}}$, $\alpha_{\mathrm{D2}}$ and $\alpha_{\mathrm{D3}}$ of diverge junctions. Indeed, we use the trajectory data to reconstruct the main paths along the network to compute the percentage of vehicles that at the junction continue towards one road rather than the other one. The analysis of the paths of heavy vehicles during one month of data allowed us to estimate $\alpha_{\mathrm{D1}} = 0.78$, $\alpha_{\mathrm{D2}} = 0.78$ and $\alpha_{\mathrm{D3}} = 0.48$.

In addition to GPS data, there are fixed sensors along the highway network that record the flow and speed of vehicles crossing it every minute. We denote the sensor data on road $r$ by ${\mathcal{Q}}_{r,{\mathrm{sens}}}$ and ${\mathcal{V}}_{r,{\mathrm{sens}}}$. As mentioned ad the end of Section 3, we use sensors data for the boundary conditions on the incoming side of the roads entering the network, i.e. roads 1, 3 and 5. Indeed, we modify the incoming numerical flux of scheme (17) in the first cell of these roads. Specifically, since the sensor data is given every minute, we interpolate it in a piecewise constant way, thus we define ${\mathcal{Q}}^{n}_{r,{\mathrm{sens}}} = Q_{r,{\mathrm{sens}}}(\hat t)$, where $\hat t = \lfloor 60n{\Delta t}\rfloor$ is the minute corresponding to $t^{n}$. Then we compute the density $\rho^{n}_{r,0}$ by applying (17) with $Q^{n}_{r,-1/2} = {\mathcal{Q}}^{n}_{r,{\mathrm{sens}}}$. Once computed the density, we estimate $w^{n}_{r,0}$ such that $v(\rho^{n}_{r,0},w^{n}_{r,0}) = {\mathcal{V}}^{n}_{r,{\mathrm{sens}}}$, with $V^{n}_{r,{\mathrm{sens}}} = {\mathcal{V}}_{r,{\mathrm{sens}}}(\lfloor 60 n{\Delta t}\rfloor)$ velocity value captured by the sensor. In Figure 13 we show the variation in time of the flux of heavy vehicles for the three sensors.

To simulate the traffic dynamics we proceed as follows: first, we estimate the initial condition starting from an empty network; then, we approximate the dynamics of vehicles by means of model (15). To estimate the initial traffic state, we start from an empty network with $w = ({w_{L}}+{w_{R}})/2$. As explained above, the traffic state on the incoming boundary of roads entering the network, $(\rho^{n}_{r,0},w^{n}_{r,0},y^{n}_{r,0})$ for $r\in\{1,3,5\}$, is estimated from sensors data. Moreover we assume homogeneous Neumann conditions on the outgoing side of roads 2, 4 and 6, hence vehicles are free to leave the network. The roads are then filled for half an hour through model (15) with both trajectory and sensors data, using the CTM scheme introduced in Section 2. The traffic state after half an hour of simulation is the initial state for model (15), and we use it to approximate the dynamics and then estimate vehicular emissions via (22).

In our test we consider the network depicted in Figure 11. We use again the CGARZ model with the flux function defined at the beginning of Section 5. Let us denote by $L_{r}$ the length of road $r$, $r = 1,\dots,6$. According to the properties of the highway, we fix $L_{1} = L_{4} = 41\,{\mathrm{km}}$, $L_{2} = L_{3} = L_{4} = L_{5} = 36\,{\mathrm{km}}$, $\beta_{M1} = \beta_{M3} = 0.2$, $\beta_{M2} = 0.5$, ${V^{\mathrm{max}}} = 90\,{\mathrm{km/h}}$ and ${\rho^{\mathrm{max}}} = 56\,{\mathrm{veh/km}}$. In particular, by assuming that the length of each heavy vehicle is $18\,{\mathrm{m}}$ (length + safety distance), the maximum density ${\rho^{\mathrm{max}}}$ is computed by dividing the total employable space (two lanes) by the space occupied by the vehicles, thus ${\rho^{\mathrm{max}}} = 56\,{\mathrm{veh/km}}$. Moreover, we set ${\rho_{f}} = 10\,{\mathrm{veh/km}}$, ${w_{L}} = g({\rho_{f}}) = 739$ and ${w_{R}} = g({\rho^{\mathrm{max}}}/2) = 1260 ,$ with $g$ defined in (24).

Once estimated the initial state, we consider a simulation with $T = 1.5\,{\mathrm{h}}$ influenced by the trajectory data represented in Figure 12. In Figure 14 we show the density of vehicles along the network at different times, excluding the six junctions and focusing only on the main roads. We observe that the inclusion of real trajectory data allows model (15) to capture the formation of small congestions, mainly on the horizontal roads, which could not be observed otherwise.

In Figure 15 we draw the variation in time of the total emissions produced along roads 2 and 3, where the dynamics are mainly influenced by the presence of GPS data. To explain the different trend of the curves, we need to look at the traffic dynamics after $60\,\min$ of simulation, see Figure 14(C). Here we observe a slowing vehicle at the end of road 1 that forces the macroscopic dynamics to slow down abruptly. Therefore, few trucks enter road 2, causing the sharp decrease in emissions shown in Figure 15 (blue line). This phenomenon is not captured by the model without GPS data (red line). Then, the congestion at the end of road 1 slowly melts away and the trucks cross the intersection towards road 2, causing the slowdown still visible at the end of the simulation, see Figure 14(D). The trend of the two emission curves, calculated with and without GPS data (blue and red line respectively), differ from then on, since the first model is able to detect significant changes in traffic speed and acceleration. On road 3, instead, the presence of real data produces higher emission values than the model without vehicle trajectories.

In Figure 16, we further investigate the dynamics on the last 15 km of road 3, drawing the density, speed, acceleration and ${\mathrm{NO_{x}}}$ emissions estimated at different times. Note that the direction of motion is from right to left, according to the network depicted in Figure 11. We observe that the presence of real trajectory data is responsible for the sudden changes in speed and acceleration that strongly influence the density and the emissions estimates.

Finally, we analyze the diffusion of ${\mathrm{NO_{x}}}$ pollutants into the air. Let us denote by $\psi$ the concentration of ${\mathrm{NO_{x}}}$ in unit of weight per unit of volume. First of all we assume that $\psi$ is constant along the $z$-axis, in order to reduce to a two-dimensional problem on a domain $\Omega$ in which the network is located. More precisely, we set $\Omega = [-L_{x},L_{x}]\times[-L_{y},L_{y}]$ where $L_{x} = 15\,{\mathrm{km}}$ and $L_{y} = 5\,{\mathrm{km}}$; hence we focus our attention on a rectangle around the junctions that involves only the four horizontal roads, see Figure 17. The domain $\Omega$ is discretized through a grid of steps ${\Delta x} = {\Delta y} = 10\,{\mathrm{m}}$; to do this, we divide each cell of the traffic dynamics ($100\,{\mathrm{m}}$ long) into 10 smaller cells that inherit the density, velocity, acceleration and emissions previously computed. Then, we use the emissions obtained from the traffic dynamics as source term of the following diffusion problem

where $\mu$ is the diffusion coefficient, fixed as $\mu = 10^{-8}\,{\mathrm{km}}^{2}/{\mathrm{h}}$ for aerosols [33]. To define the source term $S$ we exploit the emission rates $E_{r}(x,t)$ computed through (22) on the horizontal roads of the network. Since the traffic model is one-dimensional in space, we trivially extend the emissions into the $y$-axis as

where $[a_{r},b_{r}]$ and $[c_{r},d_{r}]$ are the horizontal and vertical length of the road $r$ in $\Omega$, respectively (in our one-dimensional configuration $d_r-c_r = {\Delta y}$). Then, the source term $S$ is defined as $S(x,y,t) = E(x,y,t)/{\Delta x}^{3}$.

The diffusion problem (26) is numerically solved with an explicit finite difference scheme. In Figure 18 we show the source term of ${\mathrm{NO_{x}}}$ emissions and their diffusion in air. In the top plots we observe several red peaks corresponding to high amounts of ${\mathrm{NO_{x}}}$ emissions due to the traffic dynamics. The red peaks on road 3 are connected to the high emission values shown in the last column of Figure 16. The bottom plots show the spread of ${\mathrm{NO_{x}}}$ concentration in $\Omega$. We observe that the diffusion of pollutants is slower than the source term, and is strongly influenced by microscopic information. Indeed, the presence of real data allows us to take into account variations in speed and acceleration that we would not be aware of without the trajectory data. Therefore, the use of GPS data enables us to approximate the source term well and to reproduce at macroscopic level the emission peaks due to unpredictable traffic dynamics.

6.   Conclusions

In this paper, we proposed a second order macroscopic traffic model that integrates multiple trajectory data into the velocity function as a tool to compute traffic quantities for estimating emissions. The combination of a macroscopic model with microscopic data is suggested by the computational efficiency of the former and the high accuracy of the latter. The proposed numerical tests show that, even when trajectory data are sparse, they make it possible to reproduce variations in speed and acceleration that otherwise would not be observable. As a result, we obtained more accurate approximations of emissions.

In the near future we plan to improve the study by distinguishing between light and heavy vehicles. We will consider macroscopic multi-class models and appropriate emission formulas that are calibrated with respect to the vehicle types.

Appendix A.   Technical proofs

Proof of Proposition 2.2. It is sufficient to prove that under condition (12) the method is monotone. Let

Below we check that

We have

By construction, the sending and receiving function in (6) are monotone in $\rho$, i.e. $\partial S(\cdot,\cdot,\rho)/\partial\rho \geq 0$ and $\partial R(\cdot,\cdot,\rho)/\partial \rho \leq 0$.

In the following we set $S^n_j = S(x_j,t^n,\rho^n_j)$ and $R^n_j = R(x_j,t^n,\rho^n_j)$. For $l = j-1$, we have

Similarly, for $l = j+1$

Therefore, in both cases $l = j\pm 1$, $\partial H(j,n,\boldsymbol\rho^n)/\partial \rho^n_l \geq 0$. When $l = j$, we have to check that

We analyze the following four cases:

ⅰ) If $S^n_{j-1} \leq R^n_j$ and $S^n_{j} \leq R^n_{j+1}$, then

Thus, to obtain the inequality (27), we assume

where $S^n_j = {\mathcal{F}}(x_j,t^n,\rho^n_j)$ or $S^n_j = {\mathcal{F}}^{\max}(x_j,t^n)$. This leads to the condition

ⅱ) If $S^n_{j-1} \leq R^n_j$ and $S^n_{j} > R^n_{j+1}$, then

Hence, (27) follows.

ⅲ) If $S^n_{j-1} > R^n_j$ and $S^n_{j} > R^n_{j+1}$, then

Therefore, we assume

where $R^n_j = {\mathcal{F}}(x_j,t^n,\rho^n_j)$ or $R^n_j = {\mathcal{F}}^{\max}(x_j,t^n)$. This leads again to

ⅳ) If $S^n_{j-1} > R^n_j$ and $S^n_{j} \leq R^n_{j+1}$, then

Moreover,

Hence, we need again

Summing up, if

then the operator $H(j,n,\boldsymbol\rho)$ is non-decreasing with respect to all the components of $\boldsymbol\rho$. Thus,

from which the thesis follows.

Proof of Proposition 2.3. The flux function ${\mathcal{F}}(x,t,\rho) = \rho{\mathcal{U}}(x,t,\rho)$, with ${\mathcal{U}}$ as in (2), satisfies $\partial_\rho{\mathcal{F}} = {\mathcal{U}}+\rho\partial_\rho{\mathcal{U}}$. Since $\partial_\rho{\mathcal{U}} \leq 0$ by construction, we have

Let us consider the single trajectory $p_\kappa(t)$, with $\kappa = \kappa(x,t)$ defined in (3), then

where

By simple computations we have

which depends on the selected trajectory $p_\kappa(t)$ at each point $(x,t)$. Therefore, for all $(x,t,\rho)\in{\mathbb{R}}\times{\mathbb{R}}^+\times[0,{\rho^{\mathrm{max}}}]$

and thus we can assume the CFL condition

Acknowledgments

This work was partially funded by the company Autovie Venete S.p.A. The authors acknowledge the Italian Ministry of University and Research (MUR) to support this research with funds coming from the project "Innovative numerical methods for evolutionary partial differential equations and applications" (PRIN Project 2017, No. 2017KKJP4X). The work was also carried out within the research project "SMARTOUR: Intelligent Platform for Tourism" (No. SCN_00166) funded by the MUR with the Regional Development Fund of European Union (PON Research and Competitiveness 2007-2013). The authors are members of the INdAM Research group GNCS.