Uncertainty quantification in a macroscopic traffic flow model calibrated on GPS data

1.   Introduction

Macroscopic traffic flow models consisting in (systems of) partial differential equations are used to simulate traffic flows on road networks since decades ^[1]. Yet, these models are usually fully deterministic, and the coupling with real data poses severe difficulties, which require advanced data assimilation techniques (see e.g., ^[2] and references therein) and may result in poor prediction outcomes.

In this paper, we focus on the basic Lightwill-Whitham-Richards (LWR) first order model ^[3,4], augmented with random variables in the velocity function and the initial condition to account for real data uncertainty. This model is specifically designed to cope with the traffic data set we were provided, which consists of floating car data coming from embedded GPS devices. For alternative stochastic traffic flow models, we refer the reader to ^{[5,6,7,8,9,10,11]}.

Several stochastic methods have been proposed in literature to evaluate uncertainty propagation in stochastic PDE models. For hyperbolic conservation laws, the so-called non-intrusive methods, like multilevel Monte Carlo finite volumes ^[12,13] or stochastic collocation ^[12,14], allow to use the underlying deterministic code but suffer of slow convergence rate and curse of dimensionality. On the other side, intrusive methods, like polynomial chaos expansion ^[15], require deep modifications of the deterministic simulation code and are not suitable for discontinuous solutions. The main aim of this work is to test qualitatively the semi-intrusive approach proposed by Abetaall and Congedo ^[16] in the context of traffic flow modeling. The underlying idea is to extend the spatial computational domain to the probabilistic components and to compute conditional expectations of the flux in the probabilistic direction. To evaluate the expectations of the flux, we use a piecewise polynomial approximation computed using a reconstruction method. This polynomial reconstruction is the same used for space finite volume methods, except that the measure is no longer the standard Lebesgue measure but the probabilistic one. Compared to the above mentioned approaches, the Abetaall-Congedo method requires minor modification of the deterministic code and ensure rapid convergence. Our study being qualitative, we rely on some simplifying assumptions in the choice of the stochastic parameters and the corresponding distribution functions. We observe that the approach can be easily extended to any probability distribution function, in case of correlated random variables and it is suitable to discontinuous solutions. A deeper statistical analysis of the data set would be necessary for a more rigorous study and more significant quantitative results.

The traffic data available for this research were provided by the company AUTOROUTES TRAFIC ^[17] and they are presented and treated in [18,Section 4]. They correspond to a sector of the French A8 highway, also called la Provençale, between the exit no. 45 (Antibes) and the exit no. 49 (Nice St Isidore), for a total length of 17, 5 km, see Figure 1a. In this study, we will consider the direction from Antibes to Nice St Isidore, denoted as Direction 1. Data were collected on four Tuesdays (March 19 and 26 and April 2 and 9, 2013) from 6 a.m. to 11 a.m. and are divided into two categories: GPS data and magnetic loop detector data. GPS data, supplied by COYOTE embedded systems, include the device ID, the position and the velocity of the car, sent every minute. In order to show the behaviour of data, in Figure 1b we report in a space-time plot the speeds registered on March 19.

Loop detector data are supplied by the highway operator ESCOTA ^[19] at 16 locations and consist of hourly and 6 min flux averages. They will be used for the estimation of the real amount of cars travelling on the considered section, see Section 3.2.

The rest of the paper is organized as follows. In Section 2 we recall the deterministic LWR model and the corresponding numerical scheme used later. In Section 3 we introduce the stochastic setting and the random variables of interest. The semi-intrusive approach is described and tested in Section 4. Section 5 is devoted to validation against real data and conclusion and perspectives are presented in Section 6.

2.   The deterministic LWR model

Macroscopic traffic flow models describe the evolution of the position of vehicles on a (infinite) road identified with the real line in terms of averaged quantities, such as the density $\rho = \rho(t, x)$, $t\in {\mathbb R}_+, x\in {\mathbb R}$ and the average speed of cars $v = v(t, x)$. The first model was introduces in the mid '50 by Lighthill, Whitham and Richards ^[3,4] and it is based on the conservation of the number of cars, which is expressed by the following scalar conservation law:

where $q = \rho v$ is the traffic flow. To close the equation, the LWR model assumes that $v = v(\rho)$ is a non-increasing function of the density. In this work, we will use a modified Newell-Daganzo velocity function ^[20,21], which is characterized by a linear decreasing free-flow speed and a hyperbolic velocity in congested regime, and we add a downward jump at $\rho = \rho_c$ to model the capacity drop observed in real traffic (cf. Figure 2):

where $v_{max}$ is the maximal free-flow speed, $\omega_f$ is the backward propagating wave-speed, $\rho_c$ is the critical density (the limit density between the fluid and congested phases), $\rho_{max}$ is the maximal density corresponding to a bumper-to-bumper situation and $\rho_a$ is a further parameter accounting for the capacity drop, so that

This choice has shown to best fit our data, as results from ^[18].

Approximate solutions to (2.1) can be computed via Godunov finite volume method ^[22], which is equivalent to the supply-demand method (or Cell Transmission Model) used by the transportation engineers ^[20,23,24].

Let $\Delta x$ and $\Delta t$ be respectively the space and time grid sizes, $x_i = i\Delta x$, $i\in\mathbb{Z}$, be the grid points and $C_i = [x_{i-1/2}, x_{i+1/2}[$ the space mesh cells. We aim at constructing approximate solutions

In order to iteratively define $\rho_\Delta$, the demand and supply functions are usually defined respectively as

In the case of a discontinuous flux function, the above definitions are modified as follow ^[25]:

where $q(\rho_c\pm)$ denote the left and right traces of $q$ at $\rho_c$.

The Godunov numerical flux at the interface $x_{i+1/2}$ is then defined by

and the recursive numerical scheme is given by

under the classical CFL condition

where $\lambda^n$ is the maximum of the absolute values of wave speeds at time $t^n$, see [25,Section 5] for more details. The scheme (2.6) is initialized taking

for a given initial datum $\rho_0:\mathbb{R}\to [0, \rho_{max}]$.

3.   Stochastic scalar conservation laws

Following ^[12,13,26], we will consider $N = 2$ independent random variables $X: = \{X_1, X_2\}$, which will account for uncertainty in the maximal speed and the initial datum respectively, defined on a probability space ${\mathcal P} = (\Omega, {\mathcal B}_{\Omega}, \mu)$, where $\omega: = (\omega_1, \ldots, \omega_s)\in\Omega: = \Pi_{i = 1}^s [a_i, b_i]$ is the sample space of random parameters, ${\mathcal B}_{\Omega}$ is the $\sigma$-algebra of the Borel sets on $\Omega$ and $\mu$ the Lebesgue measure. We assume that $X_i:\Omega \to {\mathbb R}$ is measurable and we denote by $f_{X_i}: {\mathbb R} \to {\mathbb R}_+$ its probability density function. Since the variable are independent, their cumulative distribution is given by $f_X(x_1, x_2) = f_{X_1} (x_1) f_{X_2} (x_2)$. Hence, for any measurable real valued function $g$, its expected value is given by

We are interested in the stochastic Cauchy problem

where $F = F(\cdot; X_1(\cdot)): (\Omega, {\mathcal B}_{\Omega}) \to \left(L^{\infty}({\mathbb R}; {\mathbb R}); {\mathcal B}(L^{\infty}({\mathbb R}; {\mathbb R}))\right)$ is the stochastic flux function and $\rho_0 = \rho_0(\cdot; X_2(\cdot)): (\Omega, {\mathcal B}_{\Omega}) \to \left(L^1({\mathbb R}; {\mathbb R}); {\mathcal B}(L^1({\mathbb R}, {\mathbb R}))\right)$ the stochastic initial condition. Since we are led to consider flux functions with jump discontinuity in $\rho$, we refer to the corresponding theory ^[27,28,29]. Let us assume that $F(\rho; X): = F_-(\rho; X) +(F_+(\rho; X)-F_-(\rho; X))H(\rho-\rho_c)$, where $H$ denotes the Heaviside function and $\tilde H$ denotes the multivalued extension of $H$ ($\tilde H (0)\in [0, 1]$). We say that $F$ is jump continuous at $\rho_c$ if the left and right limits at $\rho = \rho_c$ exist and are finite.

Definition 3.1. (Adapted from [13,Definition 3.2]) A measurable mapping $\rho : (\Omega, {\mathcal B}_{\Omega}) \to C\left({\mathbb R}_{t > t_0}; L^1({\mathbb R}; {\mathbb R})\right)$ is a random entropy solution of (3.1) if

● For $\mu$-a.e. $\omega\in\Omega$, it satisfies

for all test functions $\phi\in C^1_c ([t_0, +\infty[\times {\mathbb R}; {\mathbb R})$.

● For $\mu$-a.e. $\omega\in\Omega$ and for each convex entropy $\eta\in C^1(\mathbb{R}; \mathbb{R})$, there exists a function $w\in L^{\infty}(\mathbb{R}^+ \times \mathbb{R}; [0, 1])$ such that $w(t, x)\in \tilde H(\rho(t, x, \omega))$ a. e., it holds

for all test functions $\phi\in C^1_c ([t_0, +\infty[\times {\mathbb R}; {\mathbb R}_+)$, where

Well-posedness results for problem (3.1) in the case $F\in W^{1, \infty}({\mathbb R}; {\mathbb R})$ can be found in [13,Theorem 3.3] and [12,Theorem 3.11].

3.1.   Random velocity

To account for vehicle speed variability, we consider a stochastic velocity function in the form

expressing the perturbation from the equilibrium velocity $v$. The distribution of the perturbation depends on several factors as, for instance, the drivers behavior and the weather conditions. We use GPS real data to fit the distribution function. To this end, we plot speed values against densities and we find a spline curve interpolating the medians, separating the density domain in cells and computing the median speed of each spatial cell, see Figure 3.

We then compute the $y$-distance of each point from the spline and we analyze the normality of the distributions in congested and free flow phases separately in order to understand if we can model both distributions with a Gaussian. We remove outliers to reduce the error due to rescaling and we perform a Kolmogorov-Smirnov test on normalized data. We found p-values for free and congested flow of 0.015 and 0.170 respectively, allowing us to conclude that we have no statistical significance to reject the hypothesis of normality. The QQ-plots in Figure 4 and the empirical and normal CDFs in Figure 5 support the conclusion.

Our goal is to apply a suitable perturbation to the velocity in order to reproduce this deviation from a median value. Hence, we have to adjust the parameters of the distribution of these perturbations. Their standard deviations are respectively 0.12 and 0.23. For simplicity, we will use the same perturbation for free-flow and congested conditions. We take $\sigma = 0.2$, which corresponds to the normal distribution law $N_{(\mu = 0, \sigma = 0.2)}$. Since we need a distribution with bounded domain, we replace the normal distribution by a triangular law $T_{[-0.5, 0, 0.5]}$, which still has standard deviation equal to 0.2.

3.2.   Random initial datum

The second random variable we consider is a perturbation applied to the initial condition to compensate for the lack of information on the penetration rate of GPS data. In our specific case, the amount of equipped vehicles represents only a little percentage of the total volume of traffic. Its variation depends on the time of the day and, indirectly, on traffic density.

To estimate the percentage of equipped vehicles compared to the total traffic density, we use the corresponding real flux measurements obtained by some magnetic loop detectors. In particular, we consider measurements taken at km 180, in the direction west-est (Antibes to Nice), and we compare them with the GPS data at the same time and location, see Figure 6.

Percentages per hour and per day are collected in Table 1. We can observe that the penetration rate of GPS-equipped vehicles has the same pattern every day: early in the morning the percentage of cars equipped with a COYOTE system is very low (0.9–1.4%), then it starts growing until it reaches the maximum between 9 and 10 a.m. (2.3–3%). Finally it starts decreasing again. This is probably due to the fact that COYOTE systems are mainly used by some specific categories of people such as taxi or truck drivers.

Taking a least square interpolation of the percentages of equipped vehicles (Figure 7a), we obtain rescaled values (Figure 7b) and we can analyze the error of the GPS flux with respect to the real one. We can consequently derive the density from the flux and model its uncertainty. To compute the density from the flux, which is not a monotone function, we need to know whether the traffic is in congested or free flow phase (cf. Figure 2).

Since we know when traffic at km 180 is in free or congested flow condition, we can distinguish the flux values and compute the corresponding density through the inverse of the corresponding branch of the flux function.

Given the rescaled GPS densities and the real ones, we analyze the distribution of the relative error of the former with respect to latter. We then analyze the error distribution as a function of density in order to model a perturbation on the initial datum. To fit the data, we use as an exponential function of the density multiplied by a random variable $X_2$:

where $\alpha$ and $\beta$ are positive parameters. We set $\alpha$ and $\beta$ to model the maximal perturbation. We look for parameters so that the maximum of the perturbation contains at least the 99% of the relative errors of the density (Figure 8). We set the first constraint as a perturbation of 100% in $\rho = 0$, which implies $\beta = 1$, and the second one as a perturbation of 60% in $\rho_c$. Therefore we use

We consider a uniform distribution $f_{X_2} = U_{[-1, 1]}$. Remark that this choice ensures that the density does not attain negative values.

3.3.   A stochastic LWR model for uncertainty quatification in GPS data

We are thus led to consider problem (3.1)

and

where $\rho(t_0, x)$ is the density reconstructed from data at time $t = t_0$, see Section 5.1. In the next sections, we will explain how to estimate the space-time evolution of statistical moments for the above model, and compare the results to our data set.

4.   Numerical scheme

We aim at computing approximate solutions $\rho_\Delta (\omega)$ of (3.1) by an iterative procedure $\rho^{n+1}_\Delta (\omega) = {\mathcal H}(\rho^n_\Delta (\omega))$. Following ^[16], besides the space and time discretization introduced in Section 2, we introduce a partition of the probability space $\Omega\subset\mathbb{R}^s$, i.e., a set of $\Omega_j$, $j = 1, ..., N$, of mutually independent subsets covering $\Omega$:

We look for approximate solutions of the form

For every $\Omega_j$, its probability measure is

Therefore, we want to approximate the solution of (3.1) by the conditional expectation

Once the conditional expectation is computed, the mean and the variance of $\rho_\Delta (\omega)$ on the cell $C_i$ at time $t^n$ are given by

Since the operator ${\mathcal H}$ is nonlinear, we need to compute the conditional expectation $E[g(X)]$ of a function evaluation $g(X)$ in $\Omega_j$, given the conditional expectations $E[X|\Omega_j]$. For each $\Omega_j$, we wish to find a polynomial $P_j \in \mathbb{P}^p(\mathbb{R}^s)$ of degree $p$, defined on a stencil $S_j$, i.e., a set $S_j = \{\Omega_k \}_{k \in I_k}$ with $\Omega_j \in S_j$, such that

This means that the conditional expectation of the reconstructed polynomial $P_j$ is equal to the conditional expectation of $X$ in each $\Omega_k$ of the stencil. Once this polynomial $P_j$ is known, we can estimate

by using a quadrature formula in each $\Omega_k$.

We focus now on the case of a one-dimensional probability space $\Omega\subset \mathbb{ {\mathbb R}}$, the two-dimensional framework being a straightforward generalization.

The Godunov scheme (2.5)–(2.6) applied to (3.1) gives

As we said before, the expectation of the numerical flux can be approximated as

where $P_{i, j}$ is a piecewise polynomial reconstruction of the density in the geometric cell $C_i$ and probabilistic cell $\Omega_j$ and computed on the stencil $S_j$.

In order to compute the integral (4.4), we use the third order Gaussian quadrature:

where

Therefore, the integral (4.4) can be computed as:

The most important step is computing the polynomial reconstruction, which defines the order of the method. For any cell $\Omega_j$, we define a polynomial $P_j \in \mathbb{P}^p (\mathbb{R})$ of degree $p$, described by a stencil $S_j = \{\Omega_j, \Omega_{j_1}, ..., \Omega_{j_p}\}$, where $j_1, ..., j_p \not = j$. Since in the space variable we use a first order finite volume method, we expect that the convergence order will not change increasing the order of the polynomial in the probability space. We will test the approach with a piecewise constant and piecewise linear interpolations (using a MUSCL scheme), analyzing the convergence orders.

● $p = 0$: first order reconstruction. We take $S_j = \{ \Omega_j \}$ and constant polynomials.

● $p = 1$: ENO second order reconstruction. We evaluate two linear polynomials, and take the least oscillatory one. We introduce the average mid-points

For the cell $\Omega_j$, we define two polynomials of degree 1: $p_j^-$ is constructed using the cells $\Omega_{j-1}$ and $\Omega_j$ and $p_j^+$ is defined on $\Omega_j$ and $\Omega_{j+1}$. For $p_j^+$ we have:

such that

Since $E[x-\omega_j|\Omega_j] = 0$ by definition of $\omega_j$, we obtain the $2 \times 2$ system

and similarly forl $p_{i, j}^-(\omega) = a_{i, j}^- \dfrac{\omega-\omega_j}{\omega_{j+1}-\omega_j} + b_{i, j}^-$. From (4.7) and (4.8), we recover

We choose $P_{i, j}^n$ equal to the one that realizes the least oscillation: $\min\left\{{\left|{E(\rho_{i, j+1}^n|\Omega_{j+1}) - E(\rho_{i, j}^n|\Omega_{j})}\right|}, {\left|{E(\rho_{i, j}^n|\Omega_{j}) - E(\rho_{i, j-1}^n | \Omega_{j-1}) }\right|}\right\}$.

The procedure is summarized in Algorithm 1.

4.1.   Numerical tests

To validate the approach, we show some results concerning the stochastic conservation law (3.1) with random flux function (3.4) without capacity drop and piece-wise constant (deterministic) initial datum $\rho_0(x) = \rho_L = 10$ for $x < x_0$ and $\rho_0(x) = \rho_R = 80$ for $x > x_0$, $x_0 = 0.5$. In this case the mean and the variance can be computed analytically: for each realization of the random variable, the solution is given by a shock wave moving with speed

and we deduce

We consider both uniform and triangular probability distributions and the piece-wise constant and ENO polynomial reconstructions described above, which give similar results. Figures 9 and 10 show the results obtained with the ENO reconstruction for $\Delta x = 0.01$ and $\Delta x = 0.002$ respectively, and $N = 40$ probabilistic cells in $[-1, 1]$.

To analyze the $L_1$-convergence for the two polynomial reconstructions, we denote by $\bar\mu_{\Delta}$ and $\bar\mu$ respectively the means computed by the semi-intrusive method and analytically. We then keep fixed the space mesh, while refining the probabilistic discretization, and compute the error

Figure 11a shows that the errors corresponding to the first and second order interpolations are very close, since they are conditioned by the first order spatial approximation. Figure 11b displays the probabilistic convergence curves for $N = 2, \ldots, 2^{10}$, corresponding to different space mesh sizes $\Delta x = 2^{-n} 10^{-2}$, $n = 1, \ldots, 4$. This confirms that if the probabilistic discretization is fine enough, the main contribution comes from the spatial error.

To evaluate the computing time, we compare the Abetaall-Congedo approach using a piecewise constant polynomial reconstruction of the density with the Monte Carlo method. Figure 12 plots the time as a function of the $L^1$-error, for $\Delta x = 0.001$ and a triangular distribution $T_{[-0.5, 0, 0.5]}$ for the speed perturbation. We remark that the Abetaall-Congedo approach is significantly faster than Monte Carlo.

5.   Uncertainty quantification in traffic flow simulations

5.1.   Model calibration

We calibrate the fundamental diagram following [18,Section 5]. In particular, for the section around km 180 in direction 1, we take $v_{max} = 125$ km/h, $\omega_f = 17$ km/h, $\rho_{max} = 614$ vehicles/km, $\rho_c = 120$ vehicles/km and $\rho_a = 300$ vehicles/km.

To estimate the initial conditions at some fixed time $t_0 \geq 6:00$ am, for every space cell $C_i$, we consider a weighted average of the data measured in the hour preceding time $t_0$, with an exponentially decreasing weight (for more details see ^[18,30]):

Then, as initial condition, we compute $\rho_i^{t_0}$ as

where $\rho_i^{data}(t_k)$ is the density measured in cell $C_i$ at time $t_k$. The coefficient $a$ depends on the reliability of the data and on the variability of the traffic in the hour before $t_0$. After some tests, we decided to fix $a\in[1.5, 2]$ to give more importance at data close to $t_0$ and we increase it if we need to consider a longer period.

To avoid dealing with boundary conditions, we run simulations on a larger space domain so that information coming from upstream and downstream boundaries do not reach the domain of interest.

5.2.   Uncertainty on the mean velocity

For the speed uncertainty, we use the Abetaall-Congedo approach using piece-wise constant polynomials. The triangular distribution $T_{[-0.5, 0, 0.5]}$ is discretized in 20 cells and we set $t_0 =$ 9:00 am. We remark that when a constant interpolation is used, the approach is not very different from Monte Carlo methods, but the computational time is drastically reduced, see Section 4.1. We plot the speed evolution for $X_1 = 0$ in Figure 13a and the mean and standard deviations in Figure 13b. In both cases, we compare the simulation results with the real values. We observe that the results of the simulation fit the behavior of the actual data and most of the values fall in the standard deviation range (59% and 91% respectively at $t_0+15$ min and $t_0+30$ min).

5.3.   Uncertainty on the initial condition

We analyze here the influence of a perturbation on the initial density. The goal is to find the variance corresponding to the actual variability due to the lack of data. Unfortunately, we saw in Section 3.2 that our data cover only a very low percentage of the actual values. Consequently, when we fit the variability, we found a non-linear perturbation which has high values. This implies a loss of significance in our model. Actually, when we apply the perturbation to a real case (we use the above-mentioned example), we can notice that the standard deviation dramatically increases (Figure 14). In this case, 86% of the data falls in the standard deviation range both at $t_0+15$ min and $t_0+30$ min.

5.4.   Complete uncertainty estimation

Considering both speed and initial data uncertainties in the model, we have to deal with a two dimensional probability space. Following Section 3.3, we note $\omega_1\in\Omega^1 = \cup_j\Omega^1_j$ the speed random parameter and $\omega_2\in\Omega^2 = \cup_l\Omega^2_l$ the initial density one, and assume they are independent. Using Algorithm 1, we define

and

respectively

Remark that, by the law of total variance, we can compute

We report in Figure 15 the approximated mean and standard deviation corresponding to the section between Antibes and Saint-Laurent-du-Var on April 2, 2013, 9–9:30 am. We observe that the standard deviation is very close to the previous case in Section 5.3, where only the initial density perturbation was considered. (Now, 95% of the data falls in the standard deviation range both at $t_0+15$ min and $t_0+30$ min.) This can be explained by the fact that, in our specific case, the information on density is poor and the variability of the initial datum is very large.

5.5.   Travel time estimation

Uncertainty quantification can be extended to travel time predictions. Since each car in the data set has an identification label, we know the time needed to cover a certain distance. Besides, numerical simulations allow us to compute travel times from mean velocities, and their mean and standard deviation. We report in the following figures the travel times in a space-time plot: in red the trajectory of a selected vehicle in the data set, in blue the computed mean and in green the computed mean $\pm$ its standard deviation. We observe that the forecasts are good when the traffic is in free flow or congested, see Figures 16 and 17. On the contrary, when traffic conditions change, like in Figure 18, initial conditions are not sufficient to have good results and information about boundary inflow and outflow become essential.

6.   Conclusion

In this paper, we have proposed a stochastic macroscopic traffic flow model accounting for GPS data uncertainty. We used an efficient semi-intrusive numerical method to evaluate uncertainty propagation in traffic simulations and we tested the results against real data collected on a busy sector of the French A8 highway close to Nice. The results are qualitatively satisfactory, in particular when traffic does not change phase. From the application point of view, these preliminary results can be improved considering higher order finite volume methods for the spatial discretization, and taking into account boundary conditions and their uncertainty. Moreover, an improved analysis of the statistical distribution of the data set is mandatory for more relevant results. Besides, well-posedness of stochastic conservation laws with discontinuous flux function needs to be investigated.

Acknowledgments

This research was mainly carried out during the first author's internship at Inria Sophia Antipolis Méditerranée.

Conflict of interest

The authors declare there is no conflicts of interest.

Supplementary

Semi-intrusive deterministic algorithm