Interacting moving bottlenecks in traffic flow

1.   Introduction

Moving bottleneck models have been introduced in the engineering literature starting from the end of the last century to describe the interaction of slow moving vehicles, such as buses or trucks, with bulk traffic, see ^[9,28,30]. Further approaches have then be proposed in the applied mathematics community ^{[1,15,27,39,41]}, always following the underlying idea of using multi-scale models consisting of Partial Differential Equations (PDEs) to model the overall traffic flow evolution, coupled with Ordinary Differential Equations (ODEs) accounting for the moving bottlenecks trajectories. Despite the intrinsic difficulties arising as a result of the flux constraint induced by the presence of the slow vehicles, corresponding existence and stability results have been provided, see also ^{[17,20,31,32]}, together with numerical schemes to compute approximate solutions capturing the emerging non-classical phenomena ^{[2,3,11,14,16]}.

More recently, moving bottleneck models have been applied to traffic flow regulation by means of Autonomous Vehicles (AVs), see e.g., ^{[4,5,6,7,8,34]}. Indeed, AVs can act as endogenous actuators of moving variable speed limits, improving throughput by hindering upstream traffic. This possibility has already been the object of several theoretical studies (see e.g., ^{[13,24,26,35,40,42]} and [19,Chapter 11]), but also artificial ^[25] and real world ^[38] experiments. In particular, it has been demonstrated that even low AV penetration rates can have a significant impact on reducing congestion phenomena and the associated energy consumption. In this context, macroscopic approaches are needed to address the curse of dimensionality in control design for microscopic models ^[18].

Even if some of the above works account for the presence of multiple AVs on the road, very few of them allow for AV interactions, and only in very specific modeling frameworks ^[21,37]. In this paper, we aim at generalizing the control framework designed in ^[20,34] for the Lighthill-Whitham-Richards model ^[33,36] to multiple, possibly interacting, moving bottlenecks, as if they were distributed on several lanes and they were allowed both to merge (if on the same lane) and to overtake (if on different lanes). This can be achieved by a careful study of the corresponding Riemann-type problems, and a natural generalization of the numerical reconstruction technique presented in ^[3].

The proposed model can be used to develop a general simulation environment accounting for several AVs interacting with the surrounding traffic flow and among themselves, which in turns allows to design optimization-based control strategies having fleets of AVs as actuators.

The paper is organized as follows: Section 2 details the model framework, providing a mathematical description of the fully coupled PDE-ODE model and the interactions among AVs at the level of Riemann problems. Section 3 describes the numerical scheme and Section 4 shows the results of the numerical experiments.

2.   PDE-ODE model with interacting moving bottlenecks

Following ^[15,20], we consider the following strongly coupled PDE-ODEs system

for $\ell = 1, \ldots, N$ and $t > 0, \, x \in \mathbb R$, where $\rho = \rho(t, x) \in [0, R]$ denotes the macroscopic traffic density at time $t \ge 0$ and at position $x \in \mathbb R$, $f = f(\rho) = \rho v(\rho)$ is the (strictly convex) flux and $v \in \mathbf{C^{2}}\left([0, R]; [0, V]\right)$ is a strictly decreasing function such that $v(0) = V$ and $v(R) = 0$, which represents the average speed of cars ($R$ being the maximal vehicle density attainable on the considered road section). Above, $\rho_0\in [0, R]$ and $y_\ell^0\in \mathbb R$, $\ell = 1, \ldots, N$, are respectively the initial traffic density and AV positions, while the function $F_\alpha$ in Eq (2.1c), $\alpha \in\ [0, 1]$, represents the road capacity reduction rate due to the presence of the AV at $x = y_\ell(t)$, $\ell = 1, \ldots, N$, acting as a moving bottleneck which imposes a unilateral flux constraint at its position. For simplicity, we assume that this parameter is the same for all AVs, and it is a function of the number of lanes: $\alpha = (M-1)/M$, where $M\in {\mathbb N}$ denotes the number of lanes. (The extension to a more general framework would follow the same lines detailed in this work.) To determine the function $F_\alpha$, we consider the rescaled flux function

which is strictly concave function and such that $f_{\alpha}(0) = f_{\alpha}(\alpha R) = 0$. For every $u \in [0, V]$, we define the point $\tilde \rho_{u}$ as the unique solution to the equation $f'_{\alpha}(\rho) = u$. Moreover, for every $u \in [0, V]$, we define the function

Hence, if $\dot y(t) = u$, the function $F_\alpha$ in Eq (2.1c) is defined by

If $\dot y(t) = v(\rho(t, y(t)+))$, the inequality (2.1c) is trivially satisfied since the left-hand side is zero. Finally, the points $0 \le \check \rho_u \le \tilde \rho_u \le \hat \rho_u \le \rho_u^* \le R$ are uniquely defined by

see ^[15,20] and Figure 1. We remark that $\check \rho_V = \tilde \rho_V = \hat \rho_V = \rho^*_V = 0$.

Solutions to Eq (2.1) are defined as follows (see [20,Definition 3.1]):

Definition 2.1. The $N+1$-tuple $\left(\rho, y_1, \ldots, y_N\right)$ provides a solution to (2.1) if the following conditions hold.

1. $\rho \in \mathbf{C^{{0}}} \left(\mathbb R^+; \mathbf{L^{1}_{loc}} \left(\mathbb R; [0, R]\right)\right)$ and ${ \rm{TV}}\left(\rho(t)\right) < + \infty$ for all $t \in \mathbb R^+$;

2. $y_\ell \in \mathbf{W^{1, 1}_{loc}}(\mathbb R^+; \mathbb R)$ for $\ell = 1, \ldots, N$;

3. For every $\kappa\in \mathbb R$ and for all $\varphi\in\mathbf{C_c^{1}} (\mathbb R^2; \mathbb R^+)$ it holds

4. For a.e. $t > 0$, $f \left(\rho\left(t, y_\ell(t)\pm\right)\right) - \dot y_\ell(t) \rho \left(t, y_\ell(t)\pm\right) \le F_\alpha\left(\dot y_\ell(t)\right)$ for $\ell = 1, \ldots, N$;

5. For a.e. $t > 0$, $\dot y_\ell(t) = \min\left\{u_\ell(t), v\left(\rho\left(t, y_\ell(t)+\right)\right)\right\}$ for $\ell = 1, \ldots, N$.

As long as AV trajectories do not intersect, one can apply the theory and the numerical schemes developed in ^{[3,14,15,16,20]}. In particular, we recall that the Constrained Riemann Solver is defined as follows.

Definition 2.2. The Constrained Riemann Solver $\mathcal{R}^\alpha$ for

is the map $\mathcal{R}^\alpha:[0, R]^2\times[0, V]\to\mathbf{L^{1}_{loc}}(\mathbb R, [0, R])$ defined as follows.

1. If $f(\mathcal{R}(\rho_L, \rho_R)(u)) > F_\alpha(u)+u\, \mathcal{R}(\rho_L, \rho_R)(u)$, then

2. If $f(\mathcal{R}(\rho_L, \rho_R)(u))\leq F_\alpha(u)+u\, \mathcal{R}(\rho_L, \rho_R)(u)$, then

Above, $\mathcal{R}$ denotes the standard (i.e., without the constraint (2.2c)) Riemann solver for Eqs (2.2a)–(2.2d), i.e., the (right continuous) map $(t, x)\mapsto \mathcal{R}(\rho_L, \rho_R)(x/t)$ given by the standard weak entropy solution to Eqs (2.2a)–(2.2d).

In this paper, we are interested in describing what happens when two vehicles $AV_i$ and $AV_j$, with $i, j\in\{1, \ldots, N\}$, $i\not = j$, interact. Let us assume $y_i(t) < y_j(t)$ with constant speeds $u_i(t) = u_i > u_j(t) = u_j$ for $t\in\, ]\bar t-\varepsilon , \bar t[$, so that $y_i(\bar t) = y_j(\bar t)$ for some $\bar t > 0$, see Figure 2. Depending on the surrounding traffic density conditions, only three situations may occur in a Riemann-like setting (piece-wise constant density with at most one jump discontinuity):

● The constraint (2.1c) is enforced for $\ell = i$ (we will say that $AV_i$ is active), so that $\rho(t, x) = \hat\rho_{u_i}$ for $x < y_i(t)$ and $\rho(t, x) = \check\rho_{u_i}$ for $x > y_i(t)$. Since, at $x = y_j(t)$, it holds $f \left(\check\rho_{u_i}\right) - u_j \check\rho_{u_i} \le F_\alpha\left(u_j\right)$ (see Figure 2), the constraint is not enforced for $\ell = j$ ($AV_j$ is then inactive).

● If $AV_j$ is active, we have $\rho(t, x) = \hat\rho_{u_j}$ for $x < y_j(t)$ and $\rho(t, x) = \check\rho_{u_j}$ for $x > y_j(t)$. At $x = y_i(t)$, it therefore holds $f \left(\hat\rho_{u_j}\right) - u_i \hat\rho_{u_j} \le F_\alpha\left(u_i\right),$ thus $AV_i$ is inactive.

● If both $AV_i$ and $AV_j$ are inactive, we may have either a constant density or a classical shock coinciding with one of the AV trajectories. In any case, the density must satisfy $\rho(t, x)\in [0, \check\rho_{u_i}]\cup [\hat\rho_{u_j}, R]$.

We remark that $AV_i$ and $AV_j$ cannot be both active at the same time, unless other waves are present between them.

To detail the problem evolution after the interaction (at $t > \bar t$), we distinguish if AVs are moving in different lanes or if they are located on the same lane.

Same lane interactions.

Let us assume that the two vehicles are in the same lane, so for $t > \bar t$ the upstream vehicle must adapt to the downstream vehicle speed: we will have $y_i(t) = y_j(t)$ with $u_i(t) = u_j(t)$ for $t\geq \bar t$. Three situations may occur, see Figure 3:

● If, before the interaction, $AV_i$ was active and $AV_j$ inactive, at $t > \bar t$ the solution will consists of a classical shock joining $\hat\rho_{u_i}$ to $\hat\rho_{u_j}$, followed by a non-classical shock at $x = y_j(t) = y_i(t)$ between $\hat\rho_{u_j}$ and $\check\rho_{u_j}$ and a rarefaction wave from $\check\rho_{u_j}$ to $\check\rho_{u_i}$, see Figure 3a. In particular, after the interaction, both AVs will be active.

● If, before the interaction, $AV_i$ was inactive and $AV_j$ active, at $t > \bar t$ the solution will consists only of a non-classical shock joining $\hat\rho_{u_j}$ to $\check\rho_{u_j}$, see Figure 3b. Again, after the interaction, both AVs will be active.

● If, before the interaction, $AV_i$ and $AV_j$ were both inactive, at $t > \bar t$ the solution will consists at most of the classical shock that was already present before the interaction, see Figure 3c. After the interaction, both AVs will still be inactive.

In any case, after the interaction, the two vehicles behave as $AV_j$.

Different lane interactions.

If the two interacting vehicles travel in different lanes, for $t > \bar t$ they will follow their own trajectory and $AV_i$ will overtake $AV_j$: we will have $y_i(t) \geq y_j(t)$ with $u_i(t)\geq u_j(t)$ for $t\geq \bar t$. Three situations may occur, see Figure 4:

● If, before the interaction, $AV_i$ was active and $AV_j$ inactive, at $t > \bar t$ the solution will consists of a classical shock joining $\hat\rho_{u_i}$ to $\hat\rho_{u_j}$, followed by a non-classical shock at $x = y_j(t)$ from $\hat\rho_{u_j}$ to $\check\rho_{u_j}$, a classical shock from $\check\rho_{u_j}$ to $\hat\rho_{u_i}$ and a non-classical shock at $x = y_i(t)$ from $\hat\rho_{u_i}$ to $\check\rho_{u_i}$, see Figure 4a. In particular, after the interaction, both AVs will be active.

● If, before the interaction, $AV_i$ was inactive and $AV_j$ active, at $t > \bar t$ the solution will consists of a non-classical shock at $x = y_j(t)$ from $\hat\rho_{u_j}$ to $\check\rho_{u_j}$, followed by a classical shock from $\check\rho_{u_j}$ to $\hat\rho_{u_i}$, a non-classical shock at $x = y_i(t)$ from $\hat\rho_{u_i}$ to $\check\rho_{u_i}$ and a classical shock from $\check\rho_{u_i}$ to $\check\rho_{u_j}$, see Figure 4b. In particular, after the interaction, both AVs will be active.

● If, before the interaction, $AV_i$ and $AV_j$ were both inactive, at $t > \bar t$ the solution will consists at most of the classical shock that was already present before the interaction, see Figure 4c. After the interaction, both AVs will still be inactive.

3.   Numerical scheme

We extend the reconstruction strategy introduced in ^[3] to include moving bottleneck interactions. The scheme is composed of two parts:

Numerical approximation of Eqs (2.1a)–(2.1c).

To approximate the conservation equation (2.1a), we use a conservative finite volume scheme for the constrained hyperbolic PDE using a flux reconstruction technique at the constraint locations, which allows to capture the non-classical shocks sharply. Let $\Delta x$ and $\Delta t$ be the fixed space and time steps satisfying the standard Courant-Friedrichs-Lewy (CFL) condition ^[10]:

and set $x_{j-1/2} = j\Delta x$, $x_j = (j+1/2)\Delta x$ for $j\in \mathbb{Z}$, and $t^n = n\Delta t$ for $n\in \mathbb{N}$. The initial data $\rho_0$ is approximated by the piece-wise constant function obtained by averaging it on the discretization cells $C_j = [x_{j-1/2}, x_{j+1/2}]$, namely

Away from AV positions, Eq (2.1a) is approximated by the standard Godunov scheme ^[22], whose numerical fluxes at cell interfaces $F_{{{j+}}\frac{{1}}{2}}^{{n}} = F (\rho_j^n, \rho_{j+1}^n)$ can in this case be derived using the supply-demand formula ^[29]

where

$\rho_{cr} = \mathop {{\rm{argmin}}}\limits_{\rho \in [0,R]} f(\rho)$ being the point of maximum of the flux function $f$.

Let now the approximate $\ell$-th AV position be $y_\ell^n\in C_{m_{\ell}}$ for some $m_{\ell}\in \mathbb{Z}$, $\ell = 1, \ldots, N$. If

we assume a moving bottleneck at $\bar{x}_{m_\ell} = x_{m_\ell-1/2}+d_{m_{\ell}}^n \Delta x$ with $d_{m_{\ell}}^n = \dfrac{\check{\rho}_{u_\ell^n}-\rho_{m_{\ell}}^n}{\check{\rho}_{u_\ell^n}-\hat{\rho}_{u_\ell^n}}$. If $0\leq d_{m_{\ell}}^n\leq 1$, then $\bar{x}_{m_\ell}\in C_{m_{\ell}}$ and we set

We can then update the density by means of the conservative formula

Above, we set $u_\ell^n = u_\ell(t^n)$. If two or more AVs are located in the same cell, the moving bottlenecks are treated sequentially one after the other, starting from those which are not active (i.e., satisfy Eq (2.1c)), then processing those which are active (violating Eq (2.1c)).

Remark 1. For better resolution, we also apply the above reconstruction technique to classical shocks, as described in [3,Section 3.1].

Numerical approximation of Eq (2.1b).

To track the AV trajectories, at each time step, we update the positions $y_{\ell}^n$ using an explicit Euler scheme

If two AVs are in the same lane and in the same cell, and the upstream $AV_{\ell_1}$ moves faster than the downstream $AV_{\ell_2}$, i.e., $u_{\ell_1}^n > u_{\ell_2}^n$, then we set

(A more accurate computation of the interaction point could be implemented to improve simulation accuracy.)

Remark 2. The above procedure seems not coherent with the observation that two vehicles travelling at close speeds $u_i\approx u_j$ would act as a single bottleneck with $\alpha_{ij} = \alpha_i+\alpha_j-1$. In this perspective, another option could be to reconstruct a single front moving at speed $u_{ij} = (u_i+u_j)/2$ with $\alpha = \alpha_{ij}$ whenever two vehicles are in the same cell on two different lanes. In this case, the procedure should be the following:

If $f(\mathcal{R}(\rho_{m-1}^n, \rho_{m+1}^n)(u_{ij})) > F_{\alpha_{ij}}(u_{ij}) + u_{ij} \mathcal{R}(\rho_{m-1}^n, \rho_{m+1}^n)(u_{ij})$,

we replace $\rho_m^n$ by $\hat\rho_{u_{ij}}$ and $\check\rho_{u_{ij}}$, with the jump located at

where $d_m^n\in[0, 1]$ is given by

To reconstruct the numerical flux at the cell interface $x_{m+\frac{1}{2}}$, we compute

and we set

4.   Numerical tests

In this section we illustrate the performances of the numerical scheme described in the previous Section 3 in capturing AV interactions correctly. We first analyze the Riemann-like cases described in Section 2, to conclude with a more general test case.

We consider a three lane stretch of road of length $50$ km, parametrized by the interval $[0, 50]$ with absorbing boundary conditions. We consider the quadratic flux function proposed by Greenshields ^[23] that assumes a linear decreasing dependence of the speed on the traffic density

where $V = 140$ km/h denotes the maximal speed and $R = 400$ veh/km the maximal (bump-to-bump) density on the road. The capacity reduction ratio is set to $\alpha = 0.6$. In all simulations, we set $\Delta x = 0.2$ and $\Delta t = 0.9\Delta x/V$.

For the Riemann-like cases, we consider the following initial data:

The results are depicted in Figures 5 and 6. Except small spurious oscillations emerging from AV interactions in some cases, the solutions, and in particular non-classical shocks, are accurately reproduced.

To give a more general example, we consider the following initial data: $\rho_0(x)\equiv 200$ and

Moreover, $AV_1$ and $AV_3$ are on the same lane.

Figure 7 displays the solution, accounting for the different interaction types (queuing and overtaking). We observe that $AV_1$ and $AV_3$ are initially inactive while $AV_2$ and $AV_4$ are active. Moreover, $AV_1$ and $AV_3$ overtake the preceding vehicles and finally merge. Also, $AV_1$ is always inactive until it merges with $AV_3$, which becomes active after the interaction with $AV_4$, while $AV_4$ is always active. On the contrary, $AV_2$ becomes inactive when it reaches an high traffic density region.

5.   Conclusion

We have presented a strongly coupled PDE-ODE model describing the interaction of a small number of autonomous vehicles with bulk traffic on a multi-lane highway. The model allows for AV queuing and overtaking, expanding the range of traffic control applications of previous works ^{[4,5,6,7,8,34]}. A first study on traffic management opportunities offered by this framework is presented in ^[12], showing that low penetration rates are sufficient to reach nearly optimal improvements of a selected performance index.

Acknowledgements

The second author was partially supported by ERASMUS+/KA1 "NORTH SOUTH TRAINEESHIP" academic year 2020/21.

Conflict of interest

The authors declare there is no conflict of interest.