Hamiltonian conserved Crank-Nicolson schemes for a semi-linear wave equation based on the exponential scalar auxiliary variables approach

1.   Introduction

In this paper, we are concerned with the study of a two-dimensional multi-class traffic model. This work is placed in the constantly evolving framework of mathematical models for traffic flow. The goal of traffic models is to provide tools capable of helping traffic management, in order to optimize transport and obtain economic and environmental benefits, such as the reduction of vehicles queues and pollution.

Traffic models are divided into three main categories, which depend on the scale of observation: microscopic, macroscopic and kinetic models. Microscopic models follow the dynamics of each vehicle and are described by ordinary differential equations (ODEs), see e. g. [2,16,19,33,34,38]. Macroscopic models, based on fluid dynamics, consider aggregated quantities such as the density of vehicles and are governed by partial differential equations (PDEs), see e. g. [3,7,12,17,32,41,46]. Kinetic models [24,28,33,36,39,40] are between the previous two classes since they can be derived by microscopic models while macroscopic models can be derived by kinetic descriptions. We refer to [1,15] and references therein for a more complete review on traffic models.

In recent years, the ever-increasing amount of real data, due to new technologies, has widely influenced the research on mathematical models for traffic flow. The common goal of researchers is to exploit the available real data to build ad hoc traffic models, capable of simulating increasingly realistic scenarios. We refer to [8,9,13,37] for some inspiring examples of data-fitted traffic models. The common feature of these models is the application of vehicles trajectory data collected in datasets such as [25,44]. Datasets of this type generally contain data on multi-lane highways and are able to distinguish the type of vehicle.

The focus in this work is on macroscopic traffic models. In particular, we propose a multi-class generalization in two space dimensions of the well-known first order Lighthill-Whitham-Richards (LWR) model [32,41]. First order models such as the LWR are described by a single conservation law $\rho_{t}+f(\rho)_{x} = 0,$ where $\rho$ is the density of vehicles and $f(\rho)$ is the flux function. The shortcomings of first order models are well-known in literature, for instance the infeasible solutions with unbounded acceleration [30] or the inability to reproduce complex traffic phenomena like stop-and-go waves [29,43]. However, the extension of first order models to multi-class [4,20,45], multi-lane [23] or multi-dimensional [21] models has proven to be suitable to improve the deficiencies of the LWR model and to be able to describe also complex traffic phenomena.

As we have already mentioned, traffic datasets contain information related to multi-lane highways with different types of vehicles. Most traffic models refer to dynamics of single-lane traffic and therefore do not consider the movements related to lane changes. Our aim is to exploit now all the available data, including the line-changing behavior and the different vehicles classes. To this end, we propose an extension of the work by Herty, Fazekas and Visconti for a single-class traffic model [21] to a multi-class traffic model in two space dimensions. The most common approaches which include lane-changing are the two-dimensional models and the multi-lane models. The first approach is an emerging topic, and we refer to [6,21,22] for some examples. The second approach has been used for instance in [26,27], where the authors propose a microscopic, a kinetic and a fluid dynamic model with lane changing. Here, we stick to the two-dimensional approach and incorporate two types of vehicles interacting through the flux functions. The proposed model is then defined by the coupling of LWR-type models for two classes of vehicles in the $x$ and $y$ direction. The interaction between the two classes of vehicles is modeled by means of the flux functions which depend on the sum of vehicle densities as in [11,20]. With suitable assumptions on the flux functions, we study the two-dimensional Riemann problems and validate the model comparing the theoretical results with the solutions given by a numerical approximation of Lax-Friedrichs type. Then, we calibrate the flux and velocity functions with the German dataset [25] and compare the results of our model with real trajectories data. We also test the ability of the model of capturing vehicles overtaking.

The paper is organized as follows. In Section 2 we introduce the traffic model and study the Riemann problems. In Section 3, we describe the numerical scheme and validate the model via numerical tests. In Section 4, we calibrate the model with a German dataset and compare the results with the real trajectories of vehicles. In Section 5, we propose a modified version of the model calibrated with real data, and finally we investigate on the ability of the model to simulate vehicles overtaking compared to a multi-lane model.

2.   Two-dimensional multi-class model

In this section, we introduce the traffic model used throughout the paper. Let us consider two classes of vehicles, whose densities are denoted by $\rho$ and $\mu$, respectively. Our aim is to describe the dynamics of the two classes by means of a two-dimensional multi-class model. To this end, following [21], we introduce a LWR-type model in two dimensions for two classes of vehicles, i.e.,

where ${q^{x}}_{\rho,\mu}$ are the fluxes of $\rho$ and $\mu$ along the $x$-direction, and ${q^{y}}_{\rho,\mu}$ are the fluxes of $\rho$ and $\mu$ along the $y$-direction. Similarly to [20], we define the flux functions as

where ${c^{x}}$ and ${c^{y}}$ are parameters to be calibrated and ${r^{\mathrm{max}}}$ is the maximum density of vehicles. The velocity functions in $x$ and $y$ directions coincide for $\rho$ and $\mu$, and are defined by

Hence, ${c^{x}}$ and ${c^{y}}$ represent the maximum velocity in $x$ and $y$ direction. Note that we assume that the two classes of vehicles have the same velocity ${c^{x}}$ and ${c^{y}}$, and they have the same maximum density ${r^{\mathrm{max}}}$.

First of all, we present the properties of model (1). To simplify the notation, we normalize $\rho$ and $\mu$ in order to fix ${r^{\mathrm{max}}} = 1$. We introduce the following vectors

and matrices

Therefore, we can rewrite system (1) as

System (3) is hyperbolic if any linear combination of $A$ and $B$ is diagonalizable. Thus, for $({\kappa_{1}},{\kappa_{2}}) \in \mathbb R^{2}$, we define $C = {\kappa_{1}} A + {\kappa_{2}} B$. The eigenvalues of $C$ are

which are real for any couple $({\kappa_{1}},{\kappa_{2}})$, and they coincide if and only if $(\rho,\mu) = (0,0)$ or ${\kappa_{1}} = -{c^{y}}{\kappa_{1}}/{c^{x}}$. The associated eigenvectors are

The first eigenvalue is linearly degenerate, i.e., $\nabla\lambda_{1}\cdot \gamma_{1} = 0$, while the second one is genuinely nonlinear, i.e., $\nabla\lambda_{2}\cdot \gamma_{2} \neq 0$. The Riemann invariants are

2.1.   Two-dimensional Riemann problems

Next, we describe the two-dimensional Riemann problem [10,47] associated with (1). To simplify the computations, we introduce the variable $r = \rho+\mu$, so that problem (1) can be rewritten as

The Riemann problem in two dimensions is given by (4) with constant initial data in the four quadrants, i.e.

For convenience, we define ${\hat{f}}(r) = r{c^{x}}\left(1-r\right)$ and ${\hat{g}}(r) = r{c^{y}}\left(1-r\right)$.

Remark 1. The treatment of the two-dimensional Riemann problem proposed in [10] assumes convex flux functions $f$ and $g$. In order to recover this hypothesis in our case, it is sufficient to choose the parameters ${c^{x}}$ and ${c^{y}}$ negative. However, the concave case can be recovered from the following analysis through proper sign changes.

We look for self-similar solutions $r(t,x,y) = v(x/t,y/t)$ and therefore introduce $\xi = \tfrac{x}{t}$ and $\eta = \tfrac{y}{t}.$ We can rewrite the first equation of (4) as

which leads us to

where ${c^{x}}(1-2v) = {\hat{f}}'(v)$ and ${c^{y}}(1-2v) = {\hat{g}}'(v)$.

The set of singular points parametrized by $v$ is the straight line

Defining

the Rankine-Hugoniot jump condition is

Assuming that the normal vector $(d\eta,d\xi)$ is directed towards the positive side of the shock curves, the Oleinik's entropy condition is

for ${v_{0}}$ between ${v_{-}}$ and ${v_{+}}$.

The initial data in (5) for problem (4) in the variables $(\xi,\eta)$ is given by

The solution of problem (6) with initial data (10) is composed of elementary waves. There are five possible cases: (1) no shocks, (2) no rarefaction waves, (3) exactly one shock, (4) exactly one rarefaction wave, (5) two rarefaction waves and two shocks. In this work we skip the full details of the possible cases and refer to [47] for a detailed discussion. Let us highlight the five cases now:

(1) No shocks: This case is verified when $v_{3}<v_{2}< v_{4}<v_{1}$. Each couple $(v_{2},v_{1}), \,(v_{3},v_{4})\, (v_{1},v_{4})$ and $(v_{2},v_{3})$ is connected by rarefaction waves and the straight line ${\mathscr{S}}$ defines the points of connection between them. The solution is represented in Figure 1(A).

(2) No rarefaction waves: This case is verified when $v_{3}>v_{4}> v_{2}>v_{1}$. The couples $(v_{2},v_{1})$ and $(v_{2},v_{3})$ are connected by two shocks which collide in $A = (\gamma(v_{1},v_{2}),\nu(v_{2},v_{3}))$ while the couples $(v_{3},v_{4})$ and $(v_{4},v_{1})$ are connected by two shocks colliding in $B = (\gamma(v_{3},v_{4}),\nu(v_{1},v_{4}))$. Then, we have two shocks which connect $v_{1}$ and $v_{3}$. They start from the point $O = (\gamma(v_{1},v_{3}),\nu(v_{1},v_{3}))$ and terminate either in $A$ or in $B$. The solution is represented in Figure 1(B).

(3) Exactly one shock: This case is verified when

The first sub-case is represented in Figure 2(A). Using the Rankine-Hugoniot condition (8) it can be shown that the shock curve is concave, monotonically increasing in $(v_{1},v_{4})$, bounded by the base curve ${\mathscr{S}}$, tangentially intersects ${\mathscr{S}}$ and satisfies the entropy condition (9). A similar analysis holds for the second sub-case which is represented in Figure 2(B).

(4) Exactly one rarefaction wave: This case is verified when

The first possibility of initial data gives results similar to the previous case of exactly one shock wave. For the second initial datum there exist several sub-cases but we omit the details. In Figure 3 we show two examples for the two initial data configurations.

(5) Two shocks and two rarefaction waves: This case is verified when

The main difference between the two options of initial data is that in the first case the shock curves are not neighbors while in the second case they are neighbors. There are again several sub-cases, we only show an example of the two possible initial data sets in Figure 4.

In the next section, we introduce a suitable discretization for the multi-class model. For validation purposes of the proposed scheme, we aim to recover the theoretical results from above.

3.   Numerical discretization

The numerical analysis of system (3) can be done using the dimensional splitting approach [31] which means that the two-dimensional problem is split into two one-dimensional problems. Thus, equation (3) is split into

We uniformly divide the two-dimensional space $[0,{L^{x}}]\times[0,{L^{y}}]$ into a grid $\Omega = [0,{N_{x}}]\times[0,{N_{y}}]$ with $x$-step ${\Delta x}$ and $y$-step ${\Delta y}$. We refer to the cell of the grid as ${\Omega_{ij}}$. Defining $\lambda_{1,2}$ and $\gamma_{1,2}$ the eigenvalues of $A$ and $B$ respectively, the time step ${\Delta t}$ is determined by

Then, the time interval $[0,T]$ is divided into time steps of length ${\Delta t}$.

Starting from a given initial datum $U^{0}_{ij}$, the numerical scheme is defined by the Strang splitting as

We use the Local Lax-Friedrichs flux (also known as Rusanov flux) [5,Chapter 3] for $F$ and $G$, i.e.,

where $\alpha_{i+1/2,j}$ is the maximum modulus of the eigenvalues of the Jacobian matrix in the interval $(U_{i,j},U_{i+1,j})$.

3.1.   Validation

We now test the discretization method for the two-dimensional multi-class model (3) while comparing the numerical results to the theoretical solutions of the Riemann problems introduced in Section 2.1.

Our test setting is given by $\Omega = [-5,5]\times[-5,5]$ with ${\Delta x} = {\Delta y} = 0.02$. We fix the parameters of (2) to ${c^{x}} = {c^{y}} = -1$. As already observed in Remark 1, we fix negative parameters to recover convex flux functions. The time interval $[0,T] = [0,1]$ is divided into time steps of length ${\Delta t}$ satisfying condition (11). The initial datum for the two classes $\rho$ and $\mu$ is taken as in (10)

and $\mu_{0}(x,y) = \rho_{0}(x,y)/2$. For simplicity of notation we take the values $\rho_{i}\in\{1,2,3,4\}$ and then normalize $\rho$ and $\mu$ dividing by ${r^{\mathrm{max}}} = \rho^{\mathrm{max}}+\mu^{\mathrm{max}} = 6$.

Since we aim to recover the results of the analysis done in Section 2.1, where the plots are defined for the plane $(\xi,\eta)$ with $\xi = x/t$ and $\eta = y/t$, we note that for $t = 1$ the variables $\xi$ and $\eta$ coincide with $x$ and $y$. Therefore, we plot the contours of the numerical solution $U^{n}_{ij}$ at time $t^{n} = 1$ in order to identify the plane $(x,y)$ with the plane $(\xi,\eta)$ and for a better comparison.

As we have explained in Section 2.1, there are only five possible configurations of the solution, which are determined by the initial values $\rho_{i}$, $i = 1,\dots,4$.

(1) No shocks: We fix $\rho_{1} = 4,\, \rho_{2} = 2,\, \rho_{3} = 1$ and $\rho_{4} = 3$. As shown in Figure 5(A), we have only rarefaction waves connected by the straight line ${\mathscr{S}} = \{(x,y) \,|\, y = x\}$. The results in Figure 5(A) coincide with the theoretical solution shown in Figure 1(A).

(2) No rarefaction waves: We fix $\rho_{1} = 1,\, \rho_{2} = 2, \,\rho_{3} = 4$ and $\rho_{4} = 3$. In Figure 5(B), the points of connection between the shocks are $A = (\gamma(\rho_{1}+\mu_{1},\rho_{2}+\mu_{2}),\nu(\rho_{2}+\mu_{2},\rho_{3}+\mu_{3})) = (0.25,0.5), B = (\gamma(\rho_{3}+\mu_{3},\rho_{4}+\mu_{4}),\nu(\rho_{4}+\mu_{4},\rho_{1}+\mu_{1})) = (0.75,0)$ and $O = (\gamma(\rho_{1}+\mu_{1},\rho_{3}+\mu_{3}),\nu(\rho_{1}+\mu_{1},\rho_{3}+\mu_{3})) = (0.25,0.25)$ with $\gamma$ and $\nu$ defined in (7). The results in Figure 5(B) coincide with the theoretical solution shown in Figure 1(B).

(3) Exactly one shock: We fix $\rho_{1} = 3,\, \rho_{2} = 2,\, \rho_{3} = 1$ and $\rho_{4} = 4$. As shown in Figure 5(C), we consider the first sub-case described in Section 2.1, and the shock wave is below the straight line ${\mathscr{S}} = \{(x,y) \,|\, y = x\}$. The results in Figure 5(C) coincide with the theoretical solution shown in Figure 2(A).

(4) Exactly one rarefaction wave: We fix $\rho_{1} = 1,\, \rho_{2} = 2, \,\rho_{3} = 3$ and $\rho_{4} = 4$. The results in Figure 5(D) coincide with the theoretical solution shown in Figure 2(A), and similarly to the previous case we have that the only rarefaction wave is below the straight line ${\mathscr{S}} = \{(x,y) \,|\, y = x\}$.

(5) Two shocks and two rarefaction waves: We fix $\rho_{1} = 3,\, \rho_{2} = 1, \,\rho_{3} = 2$ and $\rho_{4} = 4$. As shown in Figure 5(E), the shock waves are not neighbors, but they are separated by the rarefaction waves and the straight line ${\mathscr{S}} = \{(x,y) \,|\, y = x\}$. The results in Figure 5(E) coincide with the theoretical solution shown in Figure 4(A).

4.   Data-driven multi-class model in 2D

In this section, we calibrate the two-dimensional multi-class model with a dataset of real trajectories data. We employ the public German dataset [25] which contains vehicle trajectories data recorded on the German motorway A3, nearby Frankfurt am Main. The analyzed area is a three lanes highway of about 900 meters in length and 12 meters in width, depicted in Figure 6. A system of five video cameras recorded the vehicles passing through the study area, collecting trajectory data for 20 minutes with a sampling period of about 0.2 seconds. We refer to [25] for a detailed description of the dataset and of the data collection method. We observe that the dataset distinguishes several types of vehicles, and particularly in this work we focus on the dynamics of cars and trucks.

4.1.   Fundamental diagrams

In order to calibrate the two-dimensional model with the German dataset, we need first to derive macroscopic quantities from the microscopic information provided by the dataset. Following [21], we describe how we derive the density of cars and trucks, $\rho$ and $\mu$, the speed in the two directions, ${u^{x}}$ and ${u^{y}}$, and the flux in the two directions, $q^{x}_{\rho,\mu}$ and $q^{y}_{\rho,\mu}$. Note that, as we have already observed in Section 2, the velocity functions coincide for $\rho$ and $\mu$ since ${c^{x}}$ and ${c^{y}}$ do not distinguish the class of vehicles.

We consider the data from the second camera from the right of Figure 6, thus we work with 20 minutes of real data. We introduce the time interval ($t_{0},t_{M})$, with $t_{0} = 0$ and $t_{M} = 20\,{\mathrm{min}}$, and uniformly divide it with a time step $dt$. Note that $dt$ is used to derive the macroscopic quantities from the microscopic ones, and is independent of the time step $\Delta t$ of the numerical scheme. We call $N_{\rho,\mu}(t_{k})$ the total number of cars and trucks at time $t_{k}$ and ${L^{x}}$ the length of the road along the main direction of travel. Then, we define

The German dataset only provides the position of vehicles with respect to the two directions, thus we need to derive the speed of vehicles from their positions. We assume that each vehicle travels at constant speed which corresponds to the slope of a linear approximation in the least square sense of the vehicle positions. We denote by $v^{x,y}_{i}$ the resulting microscopic speed of car $i$ and by $w^{x,y}_{i}$ the analogous speed of truck $i$. Since we assume that the two classes have the same speed function, we define the average speed as a function of the two classes

Finally, we combine equations (12) and (13) to define the flux functions as

Once we have the density, speed and flux data as functions of time, we aggregate them with respect to a certain time period $\widetilde T = \kappa dt$. In particular, we fix $dt = 1\,{\mathrm{s}}$ and $\kappa = 60$. In Figure 7, we show the speed-density and flux-density diagrams for the two classes of vehicles in the $x$ and $y$ directions.

The graphs show that the main direction of the flow is along the $x$-axis, according to the structure of the analyzed road, while the movements along the $y$-axis represent the lane changes. Note that both the flux and the velocity along the $y$-direction show negative values, due to the lane change that can occur in both directions. The maximum density value reached in Figure 7 is $60\,{\mathrm{veh/km}}$ for cars and $12\,{\mathrm{veh/km}}$ for trucks, and both values are much smaller than the maximum density of the road given by ${r^{\mathrm{max}}} = 400\,{\mathrm{veh/km}}$. More specifically, as we observe from Figure 7, there are more cars than trucks along the road, thus we calibrate the maximum density fixing the length of vehicles as if there are only cars on the road. Hence, we assume that the length of vehicles plus the safety distance is $7.5\,{\mathrm{m}}$, thus ${r^{\mathrm{max}}}$ is defined as

However, it should be noted that the dataset only contains data in free-flow regimes without capturing congested traffic phase.

Now we need to compute the parameters ${c^{x}}$ and ${c^{y}}$. The parameters ${c^{x}}$ and ${c^{y}}$ are chosen in order to minimize the $L^{2}$-norm between the flux functions defined in (2) and the fluxes derived from data in (14), i.e., we consider

The computation is performed using the MATLAB $\mathtt{fminbnd}$ tool, which is a specific solver for minimization problems. We obtain ${c^{x}} = 97.04$ and ${c^{y}} = -0.41$. Note that ${c^{y}}$ is negative, since the lane changes occur mainly towards the rightmost lane. The resulting speed and flux functions are shown in Figure 8. Since the flux functions depend on both $\rho$ and $\mu$, we have a family of flux and velocity functions. In particular, in Figure 8(B) we show the family of flux functions ${q^{x}_{\rho}}$ as $\mu$ changes. This means that at fixed value of $\mu$ we can move only along one of the flux curves. For instance, if $\mu = 0$, i.e. there are no trucks, then the fundamental diagram ${q^{x}_{\rho}}$ corresponds to the maximum flux curve in Figure 8(B), if $\mu = {r^{\mathrm{max}}}$ then no car can enter into the road, and thus ${q^{x}_{\rho}}\equiv0$. A similar discussion holds for the other plots of Figure 8.

We observe that the advantage of the multi-class model is that we can cover quite well the clouds of real data by means of the family of flux and velocity functions. However, we note that, since the German dataset contains data which refer only to the not congested phase of traffic, we do not have enough data to better calibrate congested traffic situations. In particular, the choice of ${c^{x}}$ and ${c^{y}}$ equal for both of the classes seems to overestimate the flux for the class $\mu$, in both the directions.

4.2.   Reconstruction of density from data

In this section we describe how to treat the microscopic data to define the initial density for the numerical scheme and the reference solution for the comparison of the results. The German dataset gives information about the position of vehicles every 0.2 seconds, thus we work with pointwise data. In order to define a density function $\rho(t,x,y)$ on a domain $\mathcal{D},$ we use a kernel density estimation, the Parzan-Rosenblatt window method [35,42]. The idea of this method is to consider the data points as a density distribution and then recover the global density by summing these distributions.

Let $N(t)$ be the number of cars at time $t$ and $(x_{i}(t),y_{i}(t))$ their positions, we define

In order to recover the smooth function $\rho$, we introduce a two-dimensional Gaussian kernel

and then define

We follow a similar procedure to estimate the density of trucks. The parameters $h_{x}$ and $h_{y}$ are bandwidths chosen in order to obtain an almost constant density profile for equidistant vehicles [14]. These parameters depend on the dimensions of the road, i.e., on a road of dimensions ${L^{x}}\times{L^{y}}$ we fix $h_{x} = {L^{x}}/20$ and $h_{y} = {L^{y}}/20$, with ${L^{x}}$ being the length of the road along the $x$-axis and ${L^{y}}$ the length along the $y$-axis.

For each video camera of the German dataset we work with records data for about $80\,{\mathrm{m}}$ in length and $12\,{\mathrm{m}}$ in width. The average speed of vehicles is such that they exit from the recording area after a few seconds. In order to test longer simulations and compare them with real data, we assume that the trajectory of each vehicle can be approximated by a linear movement. Indeed, let us consider a vehicle $i$ which crosses the road between a time interval $[t_{0},t_{1}]$. We compute the coefficients $a_{i}^{x,y}$ and $b_{i}^{x,y}$ such that we can approximate the $x$ and $y$ position as $x(t) = a_{i}^{x}+b_{i}^{x}t$ and $y(t) = a_{i}^{y}+b_{i}^{y}t$ minimizing the $L^{2}$-norm of the difference with the real positions. In this way, we are able to compute the "real" position of vehicles even when they exit the supervised area. The computed positions also allow for a comparison to the numerical results.

4.3.   Numerical test

Now, we compare the numerical simulations of model (1) with the real data computed from equation (16). The simulation refers to the data recorded by the second video camera of the German dataset.

Let us consider the domain $[0,{L^{x}}]\times[0,{L^{y}}]$ uniformly divided into a numerical grid $\Omega = [0,{N_{x}}]\times[0,{N_{y}}]$ with $x$-steps of length ${\Delta x}$ and $y$-steps of length ${\Delta y}$ during a time interval $[0,T]$ divided into time steps of length ${\Delta t}$ satisfying (11). The numerical solutions are computed by means of the numerical scheme introduced in Section 3, and they are denoted by

for cars and trucks respectively, with $x_{i} = i{\Delta x}$, $y_{j} = j{\Delta y}$ and $t^{n} = n{\Delta t}$. The ground-truth data are estimated by (16) as explained in Section 4.2, and they are denoted by

for cars and trucks, respectively. The parameters used in the following test are ${L^{x}} = 450\,{\mathrm{m}}$, ${L^{y}} = 14\,{\mathrm{m}}$, ${\Delta x} = {\Delta y} = 0.5\,{\mathrm{m}}$, $T = 5\,{\mathrm{s}}$, $h_{x} = 22.5$ and $h_{y} = 0.7$. The initial configuration of densities is recovered by the ground-truth data (16) starting from the time $\hat t = 14\,{\mathrm{s}}$ of the German dataset. Therefore, at the beginning, the numerical solution coincides with the ground-truth solution. At time $\hat t$ there are three cars and one truck along the three-lanes highway, so we analyze their dynamics. In Figure 9, we compare the contours of the ground-truth density data with the contours of the reconstructed density by the numerical simulation at the final time $T$.

The dotted lines divide the road into three lanes. After $5\,{\mathrm{s}}$ of simulation the real and the numerical configurations of density are quite similar for both populations of vehicles.

In order to better compare the numerical results with the ground-truth data, we introduce the following errors

The errors at time $T$ of the previous simulation are $E_{\rho}(T) = 0.06$ and $E_{\mu}(T) = 0.02$ computed with (17) and (18), respectively. In Figure 10 we plot the numerical errors between the numerical density and real data during $10$ seconds of simulation computed every 0.5 seconds. We observe that the error related to the truck is lower than the error related to cars and that both errors increase in time. However, they remain of order $10^{-2}$.

5.   Further data analyses

In the previous section we have seen that the two dimensional multi-class LWR (1) with the flux functions defined in (2) and calibrated with real data is able to simulate the dynamics of vehicles. However, our main assumption on the flux functions is that the coefficients ${r^{\mathrm{max}}}$, ${c^{x}}$ and ${c^{y}}$ are equal for both the classes of vehicles. This is a strong assumption, since it implies that cars and trucks have the same length and velocity. Therefore, we modify now the definition of the flux functions to differentiate more clearly the dynamics of the two classes.

5.1.   Test with real data

We consider again the German dataset [25]. Our aim is to consider different maximum density of cars and trucks, due to the different length of vehicles, and different parameters ${c^{x}}_{\rho}$, ${c^{x}}_{\mu}$, ${c^{y}}_{\rho}$ and ${c^{y}}_{\mu}$ in order to take into account velocity functions which depend on the class of vehicles. The maximum density of cars ${r^{\mathrm{max}}} = {r^{\mathrm{max}}}_{\rho} = 400 \,\mathrm{veh/km}$ coincides with equation (15), and we assume that the length of trucks is twice that of cars, hence we have ${r^{\mathrm{max}}}_{\mu} = 200 \,\mathrm{veh/km}$. We slightly modify the flux functions of (1) as

Note that in (19) the different maximum densities between cars and trucks is expressed by the term $(\rho+2\mu)/{r^{\mathrm{max}}}$.

With the introduction of different coefficients ${c^{x}}_{\rho}$, ${c^{x}}_{\mu}$, ${c^{y}}_{\rho}$ and ${c^{y}}_{\mu}$ we are able to better distinguish the behavior of the two classes of vehicles, by means of different maximum velocities for the two classes in both the directions.

Next, we repeat a procedure analogous to the one proposed in Section 4.1 to estimate the velocity functions and the fundamental diagrams. In particular, we define different velocity functions for $\rho$ and $\mu$ as

from which we recover the flux functions as

We estimate the parameters ${c^{x}}_{\rho}$, ${c^{x}}_{\mu}$, ${c^{y}}_{\rho}$ and ${c^{y}}_{\mu}$ in order to minimize the $L^{2}$-norm between the flux functions defined in (19) and the fluxes derived from data in (20) and compute

using again the $\mathtt{fminbnd}$ MATLAB tool. We obtain ${c^{x}}_{\rho} = 99.61$, ${c^{y}}_{\rho} = -0.40$, ${c^{x}}_{\mu} = 74.86$ and ${c^{y}}_{\mu} = -0.49$. Hence, the cars have a faster velocity than the trucks along the main direction of travel, while the velocity of lane-changing is quite similar between the two classes.

In Figure 11, we show the family of speed and flux functions obtained with the above described procedure. Note that the speed and flux functions related to trucks, shown Figures 11(B), 11(D), 11(F) and 11(H), are defined for $\mu\in[0,{r^{\mathrm{max}}}_{\mu}]$, with ${r^{\mathrm{max}}}_{\mu}$ being the half of ${r^{\mathrm{max}}}$. Again, we are able to cover the clouds of real data, but in this case the plots of ${q^{x}_{\mu}}$ and ${q^{y}_{\mu}}$ (Figures 11(F) and 11(H)) reach lower flux values with respect to Figures 8(C) and 8(F), according to the lower density and velocity of trucks recorded by the dataset. Hence, the overestimation of flux values for the class of trucks is highly reduced with the introduction of ${c^{x}}_{\mu}$ and ${c^{y}}_{\mu}$ compared to the results obtained in Section 4.1.

We repeat the same numerical test proposed in Section 4.3 with the new flux functions (19), estimating again the resulting errors with (17) and (18). The density plots we obtain are similar to the plots shown in Figure 9, thus we omit the picture. However, as shown in Figure 12, we obtain a better estimate of the errors compared to the test done in Section 4.3.

5.2.   Vehicles overtaking

A further investigation of our model is the testing of the ability of capturing vehicles overtaking. We consider the following flux functions

where we have different parameters ${c^{x}}_{\rho}$, ${c^{x}}_{\mu}$, ${c^{y}}_{\rho}$ and ${c^{y}}_{\mu}$ but the same maximum density ${r^{\mathrm{max}}}$. The idea is to simulate traffic dynamics with different maximum velocities for the two classes and verify if the faster vehicles are able to overtake the slower ones. Indeed, the presence of the component transverse to the main direction of motion naturally lends itself to the modeling of vehicles overtaking.

We consider a numerical grid $\Omega = [0,{N_{x}}]\times[0,{N_{y}}]$ with $x$-steps ${\Delta x}$ and $y$-steps ${\Delta y}$ during a time interval $[0,T]$ divided into time steps ${\Delta t}$ satisfying (11). In particular we work on a road with two lanes, with two cars and one truck. We fix the following parameters: ${L^{x}} = 100\,{\mathrm{m}}$, ${L^{y}} = 6\,{\mathrm{m}}$, ${\Delta x} = {\Delta y} = 0.2\,{\mathrm{m}}$ and $T = 4\,{\mathrm{s}}$. Moreover, we assume that ${c^{x}}_{\rho} = 80$ and ${c^{y}}_{\rho} = -0.4$, while ${c^{x}}_{\mu} = {c^{y}}_{\mu} = 0$, thus the truck does not move.

As we can see in Figure 13, we consider a road with two lanes with two cars and a truck. At the beginning of the simulation there is a car in the top lane and a truck in front of the other car in the bottom lane. Since the truck does not move and the cars are free to move along the $y$-axis, in Figure 13(B) we see that both cars move towards the north-east direction. In particular, the car in the top lane is leaving the road and the other one starts to overtake the truck, which acts as an obstacle along the main travel direction. Finally, Figure 13(C) shows that the car has been able to overtake the truck since it is exiting the road while the truck is still inside the domain.

5.2.1.   Comparison with a multi-lane model

Finally, we intend to compare the proposed 2D multi-class approach to a first order multi-class multi-lane model. Indeed, the inclusion of lane change dynamics well fits with the simulation of vehicle overtaking and allows us to compare the results obtained with the 2D multi-class model. Specifically, we extend the multi-lane LWR model proposed in [23] to a multi-class model. Hence, let us consider a road with two lanes and two classes of vehicles $\rho$ and $\mu$. The dynamics on the two lanes is described by

where $\rho^{1}$, $\mu^{1}$ and $q^{1}_{\rho,\mu}$ are the densities and the flux function of the two classes along lane 1, and $\rho^{2}$, $\mu^{2}$ and $q^{2}_{\rho,\mu}$ along lane 2. The functions $S_{\rho,\mu}$ regulate the lane changing and are defined as

where $u^{1,2}_{\rho,\mu}$ are the velocity functions related to $\rho$ and $\mu$ respectively along the two lanes and $C$ is a constant.

The flux functions are chosen similar to the ones used for the two-dimensional multi-class model. Therefore, we define them as

where ${r^{\mathrm{max}}}$ is the maximum density of the two classes. In order to compare such a model with the results obtained with our multi-class two-dimensional model from Section 5.2, we replicate an analogous test. Indeed, we consider a road $[0,L]$ with two lanes along which there are two cars and a truck during a time interval $[0,T]$. The parameters of the test are $L = 100\,{\mathrm{m}}$, ${\Delta x} = 0.2\,{\mathrm{m}}$, $T = 4\,{\mathrm{s}}$, $C = 1$ in (22), $c_{\rho} = 80$ and $c_{\mu} = 0$. We use a Godunov scheme [18] to approximate problem (21).

In Figure 14, we show the density of cars and trucks on the two lanes at different times. The plots in the first row show Lane 1, with a car and no trucks, the plots in the second row show Lane 2, with a car and a truck which does not move at all during the simulation since its velocity is 0. The source terms in (21) allow lane changing even if the traffic is not congested and vehicles are free to move along their lane. Hence, in Figure 14(E) we see that the density of cars increases at the end of Lane 2, due to the source term, while in Figures 14(B) and 14(C) the density of cars related to Lane 1 increases since the truck does not move and the cars change lane. At the end of the simulation the density of cars is higher on Lane 1 while it is close to 0 on Lane 2.

From the two tests proposed in Sections 5.2 and 5.2.1 we observe that the 2D multi-class model seems to be more suitable for capturing the overtaking of vehicles. Indeed, the two-dimensional description seem to fit better to such a dynamics, as shown in Figures 13 and 14.

6.   Conclusions

In this work we have introduced a two-dimensional multi-class traffic model. We have analyzed the two-dimensional Riemann problems related to our model and provided numerical validations with a numerical scheme based on dimensional splitting. Then, we have analyzed the model with a dataset of real trajectories data, focusing on the dynamics of cars and trucks. The dataset has been used to calibrate the flux and velocity functions and to compare the numerical results with ground-truth data. The numerical tests have shown the good approximation of the trajectories with our model, obtaining a numerical error of $10^{-2}$. We have improved the results modifying the flux functions in order to consider different maximum velocity values for the two classes of vehicles. Finally, we have tested the ability of the model to simulate vehicles overtaking, also in comparison with a first order multi-class multi-lane model.

Future investigations will aim at improving the proposed approach with second order models and deriving the corresponding microscopic model as a system of ordinary differential equations.