Modeling random traffic accidents by conservation laws

TL;DR

This paper presents a stochastic traffic flow model incorporating accidents, based on conservation laws, with theoretical analysis and numerical methods to understand the impact of accidents on traffic dynamics.

Contribution

It introduces a novel stochastic model for traffic accidents using conservation laws, with bounded variation analysis and a numerical solution algorithm.

Findings

Total variation remains bounded over finite time

The model effectively captures random accidents in traffic flow

A numerical algorithm for simulating the model is developed

Abstract

We introduce a stochastic traffic flow model to describe random traffic accidents on a single road. The model is a piecewise deterministic process incorporating traffic accidents and is based on a scalar conservation law with space-dependent flux function. Using a Lax-Friedrichs discretization, we show that the total variation is bounded in finite time and provide a theoretical framework to embed the stochastic process. Additionally, a solution algorithm is introduced to also investigate the model numerically.

Equations

$$F^{p,s,c}(x,\rho )=c_a(x,p,s,c)c_{road}(x)f(\rho ).$$

$${\rho }_t+(F^{p,s,c}(x,\rho ){)}_x=0,\rho (x,0)={\rho }_0(x),$$

$$\mathrm{TV}(\Psi (\rho (\cdot ,t)))<\infty \text{ for all }t\in [0,T],\Psi (\rho )=\mathrm{sgn}(\rho -{\rho }^{\ast })\frac{f({\rho }^{\ast })-f(\rho )}{f({\rho }^{\ast })}.$$

$$a,a^{\prime },f,f^{\prime }\in L^{\infty }(\mathbb{R}),a^{\prime },a^{\prime \prime }\in L^1(\mathbb{R}),{\rho }_0\in \mathrm{BV}(\mathbb{R}).$$

$${\rho }_i^{j+1}={\rho }_i^j-$$

$$\displaystyle-(\frac{1}{2\lambda}(\rho_{i-1}^{j}-\rho_{i}^{j})+\frac{1}{2}(a_{i-1}f(\rho_{i-1}^{j})+a_{i}f(\rho_{i}^{j}))\Big{)}.$$

$$\lambda \Vert a{\Vert }_{\infty }\Vert f^{\prime }{\Vert }_{\infty }\le 1,$$

$$|{\rho }_i^{j+1}|$$

$$=\frac{1}{2}|{\rho }_{i+1}^j+{\rho }_{i-1}^j-\lambda (a_{i+1}{f}^{\prime }({\xi }_i)({\rho }_{i+1}^j-{\rho }_{i-1}^j)+f({\rho }_{i-1}^j)(a_{i+1}-a_{i-1}))|$$

$$\le \frac{1}{2}\left(|{\rho }_{i+1}^j|(1-\lambda a_{i+1}{f}^{\prime }({\xi }_i))+|{\rho }_{i-1}^j|(1+\lambda a_{i+1}{f}^{\prime }({\xi }_i))\right)+\frac{\lambda }{2}f({\rho }_{i-1})|a^{\prime }({\eta }_i)|2\Delta x$$

$$\le \Vert {\rho }^j{\Vert }_{\infty }+\Delta t\Vert f{\Vert }_{\infty }\Vert a^{\prime }{\Vert }_{\infty }.$$

$$\mathrm{TV}({\rho }^{j+1})$$

$$+({\rho }_{i+1}^j-{\rho }_i^j)-\lambda (a_{i+1}f({\rho }_{i+1}^j)-a_{i}f({\rho }_i^j))|$$

$$=\frac{1}{2}\sum _{i\in \mathbb{Z}}|({\rho }_{i-1}^j-{\rho }_{i-2}^j)(1+\lambda a_{i-1}{f}^{\prime }({\xi }_{i-\frac{3}{2}}^j))+\lambda f({\rho }_{i-2}^j)\Delta x{a}^{\prime }({\eta }_{i-\frac{3}{2}})$$

$$+({\rho }_{i+1}^j-{\rho }_i^j)(1-\lambda a_{i+1}{f}^{\prime }({\xi }_{i+\frac{1}{2}}^j))-\lambda f({\rho }_i^j)\Delta x{a}^{\prime }({\eta }_{i+\frac{1}{2}})|$$

$$\le \frac{1}{2}\sum _{i\in \mathbb{Z}}|{\rho }_{i-1}^j-{\rho }_{i-2}^j||1+\lambda a_{i-1}{f}^{\prime }({\xi }_{i-\frac{3}{2}}^j)|$$

$$+\frac{1}{2}\sum _{i\in \mathbb{Z}}|{\rho }_{i+1}^j-{\rho }_i^j||1-\lambda a_{i+1}{f}^{\prime }({\xi }_{i+\frac{1}{2}}^j)|$$

$$+\frac{\lambda \Delta x}{2}\sum _{i\in \mathbb{Z}}|a^{\prime }({\eta }_{i+\frac{1}{2}})f({\rho }_i^j)-a^{\prime }({\eta }_{i-\frac{3}{2}})f({\rho }_{i-2}^j)|.$$

$$|a^{\prime }({\eta }_{i+\frac{1}{2}})f({\rho }_i^j)-a^{\prime }({\eta }_{i-\frac{3}{2}})f({\rho }_{i-2}^j)|$$

$$\le$$

$$\mathrm{TV}({\rho }^{j+1})\le (1+\Delta t\Vert a^{\prime }{\Vert }_{\infty }\Vert f^{\prime }{\Vert }_{\infty })\mathrm{TV}({\rho }^j)+\Delta t\frac{3}{2}\Vert f{\Vert }_{\infty }\Vert a^{\prime \prime }{\Vert }_1.$$

$$\mathrm{TV}({\rho }^{j+1})$$

$$=:C_1.$$

$$\mathrm{TV}({\rho }^{j+m})-\mathrm{TV}({\rho }^j)$$

$$\le \Delta t\sum _{k=0}^{m-1}(\Vert a^{\prime }{\Vert }_{\infty }\Vert f^{\prime }{\Vert }_{\infty }\mathrm{TV}({\rho }^{j+k})+\frac{3}{2}\Vert f{\Vert }_{\infty }\Vert a^{\prime \prime }{\Vert }_1)$$

$$\le m\Delta t(C_1\Vert a^{\prime }{\Vert }_{\infty }\Vert f^{\prime }{\Vert }_{\infty }+\frac{3}{2}\Vert f{\Vert }_{\infty }\Vert a^{\prime \prime }{\Vert }_1).$$

$$|\mathrm{TV}({\rho }^{j+m})-\mathrm{TV}({\rho }^j)|\le \widetilde{C_1}|t-\tilde{t}|.$$

$${\mathrm{TV}}_{\mathbb{R}\times [0,T]}(\rho )=\sum _{j=0}^{\frac{T}{\Delta t}}\Delta t\mathrm{TV}({\rho }^j)+\sum _{i\in \mathbb{Z}}\Delta x\sum _{j=0}^{\frac{T}{\Delta t}-1}|{\rho }_i^{j+1}-{\rho }_i^j|.$$

$$\sum _{j=0}^{\frac{T}{\Delta t}}\Delta t\mathrm{TV}({\rho }^j)\le T{C}_1.$$

$$|{\rho }_i^{j+1}-{\rho }_i^j|$$

$$=\frac{1}{2}|({\rho }_{i+1}^j-{\rho }_i^j)(1-\lambda a_{i+1}{f}^{\prime }({\xi }_{i+\frac{1}{2}}))-({\rho }_i^j-{\rho }_{i-1}^j)(1+\lambda a_i{f}^{\prime }({\rho }_{i-\frac{1}{2}}))$$

$$-\lambda \Delta x(f({\rho }_i^j)a^{\prime }({\eta }_{i+\frac{1}{2}})+f({\rho }_{i-1}^j)a^{\prime }({\eta }_{i-\frac{1}{2}}))|$$

$$\le \frac{1}{2}|({\rho }_{i+1}^j-{\rho }_i^j)|(1-\lambda a_{i+1}{f}^{\prime }({\xi }_{i+\frac{1}{2}}))+\frac{1}{2}|{\rho }_i^j-{\rho }_{i-1}^j|(1+\lambda a_i{f}^{\prime }({\xi }_{i-\frac{1}{2}}))$$

$$+\frac{1}{2}\lambda \Delta x(f({\rho }_i^j)|a^{\prime }({\eta }_{i+\frac{1}{2}})|+f({\rho }_{i-1}^j)|a^{\prime }({\eta }_{i-\frac{1}{2}}))|),$$

$$\sum _{i\in \mathbb{Z}}|{\rho }_i^{j+1}-{\rho }_i^j|$$

$$\le C_1+\lambda \Vert f^{\prime }{\Vert }_{\infty }\Vert a^{\prime }{\Vert }_1$$

$$\sum _{i\in \mathbb{Z}}\Delta x\sum _{j=0}^{\frac{T}{\Delta t}-1}|{\rho }_i^{j+1}-{\rho }_i^j|$$

$$=:C_2.$$

$$C_F:={\int }_{\mathbb{R}}{F}^{p,s,c}(x,\rho (x))dx\le \Vert a{\Vert }_{\infty }\Vert v{\Vert }_{\infty }{\int }_{\mathbb{R}}\rho (x)dx$$

$${\mu }_{p,s,c,\rho }^F(B)={\int }_B\frac{1}{C_F}{F}^{p,s,c}(x,\rho (x))dx$$

Full text

AMS Classification: 35L65, 60J25, 90B20

Keywords: conservation laws, traffic flow, random accidents, piecewise deterministic processes

1 Introduction

Macroscopic traffic flow models based on hyperbolic conservation laws have been intensively investigated during the last decades, see [13, 14] for an overview. The various research directions include theoretical and numerical investigations such for instance well-posedness [4], coupled models [6], network extensions [14, 19], optimal control [16], or more recently, data-driven approaches [10] while stochastic traffic models have been less considered [20, 30].

Typically, macroscopic traffic flow equations are either characterized by first-order models for the evolution of the traffic density or second-order models, where an additional equation for the velocity is considered. So far, the modeling of traffic accidents (or incidents) has been considered in a deterministic setting [18, 27, 28], queueing theory approaches [3, 22] or kinetic models [12]. There are only a few contributions, where the presence of accidents is described by a stochastic process [22].

Therefore, the aim of this paper is to combine the stochastic modeling of accidents with the Lighthill-Whitham-Richards (LWR) model [25] of first-order type and to provide a framework that allows for theoretical and numerical studies. The idea is to include random effects directly in the flux function such that failures depend on the current traffic density.

We assume that accidents happen at random times and have an impact on the road capacity around the accident. Based on the LWR model, we incorporate these accidents by a space-dependent flux function determining the deterministic structure between the random accidents. Obviously, the profile of the traffic density has an impact on the probability of an accident. For instance, fluctuations in the density lead to different velocities of the cars and an accident is more likely as it is the case for stationary traffic situations. The traffic density does not only influence the probability of an accident. It also indicates where an accident could happen as for example at the end of a traffic jam. In order to capture these ideas, we face two building blocks, i.e. the deterministic dynamics between accidents and the stochastic nature, which interrupts the deterministic flow at random times. This directly leads to the well-known piecewise deterministic processes (PDPs), see [7, 21]. In [15], the latter idea has been used to incorporate random machine failures of machines based on hyperbolic dynamics, where the product density influences machine failures and vice versa. Compared to [15], we face different challenges here: First, we deal with a nonlinear dynamics with a space-dependent flux function, which does not admit total variation bounds in general and we prove under which conditions we can guarantee these bounds. Second, the position of an accident depends on the current density, which makes the modeling more involved. Additionally, the classical thinning algorithm, see [24], to sample the times of an accident might lead to large computational costs.

There are different works about hyperbolic equation based dynamics connected to randomness as for example random velocity fields [1, 26] and propagation of uncertainty [11]. However, in these works, there is no influence of the conserved quantity on the stochastic nature, i.e. no bi-directional relation between the deterministic and stochastic ideas.

The paper is organized as follows: in Section 2, we present the modeling of accidents within the LWR model and show that the total variation of the new model is bounded. Furthermore, the stochastic process is characterized such that accident probabilities can be embedded. In Section 3, a stochastic solution algorithm based on a Lax-Friedrichs discretization is introduced to analyze the occurrence of traffic accidents from a numerical point of view.

2 Modeling of accidents

We introduce how accidents as capacity drops can be incorporated into the LWR model. As we will see, this leads to a conservation law with space-dependent flux function. The latter equation is then extended to the possibility of a single (or multiple) random accidents.

2.1 General setting

Let $f\colon[0,1]\to[0,\infty)$ be a function of LWR type, i.e. $f(\rho)=\rho(1-\rho)$ with $f(0)=f(1)=0$, $f^{\prime\prime}\leq\underline{c}<0$ for some $\underline{c}<0$ and a unique $\rho^{\ast}\in(0,1)$ such that $f^{\prime}(\rho^{\ast})=0$. To describe the capacities of the road, we assume a function $c_{road}\colon\mathbb{R}\to\mathbb{R}_{>0}$ and use $c_{road}(x)f(\rho)$ as space-dependent flux. An appropriate choice for $c_{road}$ might be piecewise constant, describing the dependency of speed limits or the number of lanes.

We interpret an accident on a road as capacity reduction within an interval $I(p,s)\subset(p-s,p+s)$ of length $s$, where $p\in\mathbb{R}$ denotes the position and $s\in\mathbb{R}$ the size of the accident. The amount of capacity reduction is denoted by $c\in[0,c_{\text{max}}]$ with $0\leq c_{\text{max}}<1$ such that the road capacity at $p$ is given by $(1-c)c_{road}(p)$. We denote by $x\mapsto c_{a}(x,p,s,c)$ the capacity function of the accident. Then, it is natural to define the space-dependent flux function

[TABLE]

Altogether, we end up with the following Cauchy problem

[TABLE]

which admits a unique entropy solution, see [5] if $\operatorname{TV}(c_{a}(\cdot,p,s,c)c_{road}(\cdot))<\infty$, $\operatorname{TV}(\rho_{0})<\infty$ and if $c_{a}(\cdot,p,s,c)c_{road}(\cdot)$ is differentiable with except of finitely many points. Additionally, we need that

[TABLE]

This does not imply $\operatorname{TV}(\rho)<\infty$, which we will need in the modeling of stochastic accidents later. However, the following lemma provides conditions on the data such that the solution to the scalar conservation law (1) remains in $\operatorname{BV}(\mathbb{R})$.

Lemma 2.1**.**

Let $a(x):=c_{a}(x,p,s,c)c_{road}(x)$ satisfy $a\in C^{2}(\mathbb{R})\cap\operatorname{TV}(\mathbb{R})$ and let $f$ be an LWR flux. Furthermore, we assume

[TABLE]

Then there exists a constant $C=C(T,\|a\|_{\infty},\|a^{\prime}\|_{\infty},\|a^{\prime\prime}\|_{1},\|f\|_{\infty},\|f^{\prime}\|_{\infty},\operatorname{TV}(\rho_{0}))$ such that the solution to (1) satisfies $\operatorname{TV}(\rho(t))\leq C$ for all $t\in[0,T]$ and $\|\rho(t)\|_{\infty}\leq\|\rho_{0}\|_{\infty}+T\|a^{\prime}\|_{\infty}\|f\|_{\infty}$. Additionally, the mapping $t\mapsto\operatorname{TV}(\rho(t))$ is Lipschitz continuous on $[0,T]$.

Proof.

We prove the lemma by using the Lax-Friedrichs scheme given by

[TABLE]

The convergence of the Lax-Friedrichs scheme has been studied in [23], whereas in [31, 32] the Godunov scheme has been examined. For our purpose, the Lax-Friedrichs scheme is a suitable choice avoiding the study of various cases as needed for the Godunov scheme. We start with the $L^{\infty}$ estimate followed by the $\operatorname{BV}$-estimate and conclude that the numerical scheme converges to the unique solution of the Cauchy problem.

$\bullet$ $L^{\infty}$ estimate. Using the CFL condition

[TABLE]

we deduce that

[TABLE]

The latter implies $\|\rho^{j}\|_{\infty}\leq\|\rho_{0}\|_{\infty}+T\|f\|_{\infty}\|a^{\prime}\|_{\infty}$.

$\bullet$ $\operatorname{BV}$ estimates. Using the same arguments as in the $L^{\infty}$ estimates, we can estimate the spatial $\operatorname{BV}$ bound as follows:

[TABLE]

Using the CFL condition and

[TABLE]

yields

[TABLE]

Hence, we have

[TABLE]

Furthermore, we deduce the following bound on the time difference of the total variation

[TABLE]

If $t=j\Delta t$ and $\tilde{t}=(j+m)\Delta t$, then

[TABLE]

In order to use a compactness argument for the numerical scheme to converge, we need the total variation in space and time. For piecewise constant function $\rho$ it holds

[TABLE]

We can directly estimate the first expression by

[TABLE]

To analyze the second expression we start with

[TABLE]

where we use the CFL condition and $f\geq 0$. This leads to

[TABLE]

and therefore

[TABLE]

Let $(\Delta t_{n})_{n\in\mathbb{N}}$ be a sequence, which converges to zero and $\Delta x_{n}=\frac{\Delta t_{n}}{\lambda}$ be the corresponding spatial discretization, satisfying the CFL condition. The constructed sequence of piecewise constant functions $(\bar{\rho}_{n})_{n\in\mathbb{N}}$ has a subsequence $(\bar{\rho}_{n_{l}})_{l\in\mathbb{N}}$, which converges to some $\rho\in\operatorname{BV}(\mathbb{R}\times[0,T])$ in $L^{1}_{loc}(\mathbb{R})$ by Helly’s theorem. A Kruzkov type inequality, see [23], and a Lax-Wendroff type argument show that $(\bar{\rho}_{n})_{n\in\mathbb{N}}$ converges to a weak entropy solution, which is unique by [5]. Consequently, the limiting solution is the solution to the IVP satisfying claimed properties of the lemma. ∎

Hence, we are now able to mathematically introduce traffic accidents as partial road capacity drops via the function $a$.

2.2 Random traffic accidents

The parameters to incorporate a traffic accident in equation (1) are the position $p$, the size $s$ and the capacity drop $c$. From the modeling perspective the position is the first parameter to consider since there exists a dependency on the current traffic situation: if there are no cars, or cars are fully stopped by a traffic jam, we expect no accident, whereas if cars drive with high speed and the density is high at the same time, we expect a higher probability of an accident. Also, we observe accidents at the end of traffic jams. To summarize, the following modeling ideas should be included:

a higher distance between cars at lower speed implies a lower accident probability and vice versa, 2. 2.

a higher accident probability at increasing density (as for example tailbacks).

Regarding 1. The flow exactly describes the combination of density, i.e. car distances, and velocities such that at places where $\rho=\rho^{\ast}$ the probability of an accident can be assumed to be the most highest. This idea corresponds to a probability capturing random accidents caused by human failures solely (i.e. excluding tailbacks). If $v$ is uniformly bounded on $[0,1]$, the normalizing constant

[TABLE]

is finite and we can define the family of probability measures

[TABLE]

for $\rho\in\operatorname{BV}(\mathbb{R})=\{\rho\in L^{1}(\mathbb{R})\colon\operatorname{TV}(\rho)<\infty,\rho\in[0,1]\}$ and $B\in\mathcal{B}(\mathbb{R})$, where the latter denotes the Borel $\sigma$-algebra on $\mathbb{R}$. Here, we assume $\|\rho_{0}\|_{1}>0$ then it follows $C_{F}\neq 0$ by assumptions on $F^{p,s,c}$. The probability measure $\mu^{F}_{p,s,c,\rho}$ exactly describes the probability distribution of the position of an accident caused by the flows.

Regarding 2.: In 1. only the information of the flow is used to specify the probability of the position of an accident. Here, we incorporate the fact that at ends of tailbacks the probability of an accident is much higher, i.e. if the derivative of $\rho$ is positive. Generally, for $\rho\in L^{1}(\mathbb{R})$ we can not assign a proper derivative $D\rho$ but if $\rho\in\operatorname{BV}(\mathbb{R})$ we can argue as follows: on the one hand, a classical derivative of $\rho\in\operatorname{BV}(\mathbb{R})$ does not exist but on the other hand, the derivative of $\rho$ corresponds to a signed Radon measure $D\rho$ by a consequence of Riesz representation theorem. Furthermore, it holds for $\rho\in L^{1}(\mathbb{R})$ that

[TABLE]

where $|D\rho|$ is the total variation of the measure $D\rho$ and is given by

[TABLE]

In the latter equation we used the Hahn decomposition, i.e. there exists a measurable set $\tilde{B}\in\mathcal{B}(\mathbb{R})$ such that $\rho^{+}(B)=D\rho(B\cap E)\geq 0$ and $\rho^{-}(B)=-D\rho(B\cap(\mathbb{R}\setminus E))\geq 0$ satisfy $D\rho(B)=D\rho^{+}(B)-D\rho^{-}(B)$ for every $B\in\mathcal{B}(\mathbb{R})$. For further details, we refer the reader to [2, 9, 17, 29].

A natural probability measure for $\rho\in\operatorname{BV}(\mathbb{R})$ to describe positions of potential accidents caused by increasing densities is then given by

[TABLE]

for every $B\in\mathcal{B}(\mathbb{R})$ provided $D\rho^{+}(\mathbb{R})>0$.

Summarizing, we define

[TABLE]

for some fixed $\beta\in[0,1]$. That means, if $\beta=1$, the influence of increasing densities is neglected (end of tailbacks) and if $\beta=0$, only the latter effect is incorporated. The case $D\rho^{+}(\mathbb{R})=0$ means that there is no increasing part in the function $\rho$, which implies together with $\rho\in L^{1}(\mathbb{R})$ and $\operatorname{TV}(\rho)<\infty$ that only $\rho=0$ can fulfill $D\rho^{+}(\mathbb{R})=0$.

We only have discussed the probability distribution for the position $p$ of the accidents so far. We assume that the size $s$ follows the probability distribution $\mu^{size}$ on $(\mathbb{R},\mathcal{B}(\mathbb{R}))$ and the capacity reduction $c$ follows $\mu^{cap}$ on $([0,1),\mathcal{B}([0,1)))$. In a natural way, we collect the details using the product space

[TABLE]

with norm

[TABLE]

for $y=(p,s,c,\rho)\in E$ to define a Banach space $E$. Furthermore, we denote by $\mathcal{E}=\sigma(E)$ the smallest $\sigma$-algebra generated by the open sets induced by the norm $\|\cdot\|_{E}$. Finally, we define for every $y\in E$ and every $B\in\mathcal{E}$ the product measure

[TABLE]

where $\epsilon_{z}$ is the Dirac measure with unit mass in $z$. Since $\eta(y,B)$ describes the transition from no accident to one accident, we expect $\eta$ to be a kernel as the following lemma shows.

Lemma 2.2**.**

Let $(p,s,c)\mapsto\int_{\mathbb{R}}c_{a}(x,p,s,c)dx$ be continuous. Then $\eta$ defines a Markovian kernel on $(E,\mathcal{E})$, which additionally satisfies $\eta(y,\{y\})=0$ for every $y\in E$ if either $\mu^{size}(\{s\})=0$ for all $s$ or $\mu^{cap}(\{c\})=0$ for all $c\in[0,1)$.

Proof.

Let $y\in E$, then $B\mapsto\eta(y,B)\geq 0$ is a measure and also $\eta(y,E)=1$ by construction. Given a set $B\in\mathcal{E}$, the mapping $y\mapsto\eta(y,B)$ is measurable if $y\mapsto\mu_{y}^{pos}$ is measurable since $\epsilon_{\rho}$ is measurable in $\rho$. We have $\eta(y,\{y\})=\mu_{\rho}^{p}(\{p\})\mu^{s}(\{s\})\mu^{c}(\{c\})\epsilon_{\rho}(\{\rho\})=0$.

It remains to show that $\rho\mapsto\mu^{p}_{\rho}$ is measurable. For every $B\in\mathcal{B}(\mathbb{R})$ one verifies for $\rho\neq 0$ that

[TABLE]

Take $y=(p,s,c,\rho)$, $\tilde{y}=(\tilde{p},\tilde{s},\tilde{c},\tilde{\rho})\in E$, satisfying $\rho$, $\tilde{\rho}\neq 0$. We deduce

[TABLE]

We also have

[TABLE]

Hence, the mapping $y\mapsto\mu^{pos}_{y}(B)$ is continuous and therefore measurable. ∎

So far, we only have specified the probability distribution of a jump in the case that a jump occurs. To construct the time of a jump, or accident, we additionally need information about how likely a jump at time $t$ is. This can be done with rate functions and is based on the ideas of a marked point process, or, deterministic Markov processes, see [7, 21].

A possible choice for a rate function $\psi\colon E\to(0,\infty)$ is given by

[TABLE]

where $\lambda^{F},\lambda^{D}>0$ scale the influence of accidents caused by high fluxes and ends of tailbacks, respectively. For fixed $y=(p,s,c,\rho)\in E$, the rate $\psi(y)$ is finite. More precisely, if $\bar{\rho}(x,t)$ is a weak entropy solution to the IVP (1), then for $a(x)=c_{a}(x,p,s,c)c_{road}(x)$ it holds that

[TABLE]

We have to keep in mind that for $y=(p,s,c,\rho)\in E$ the values $\|a\|_{\infty},\|a^{\prime}\|_{\infty},\|a^{\prime\prime}\|_{1}$ might differ. We know that $\|a\|_{\infty}=1$ and $a^{\prime\prime}=0$ for all $x\in\mathbb{R}\setminus I(p,s)$ by assumption. Hence, $\|a^{\prime\prime}\|_{1}\leq|I(p,s)|\||a^{\prime\prime}\|_{\infty}$. Therefore, we assume $a\in C^{2}(\mathbb{R})$, cf. Lemma 2.1.

Let $\phi\colon E\to E$ be the deterministic evolution, i.e.

[TABLE]

where $\rho(t)$ is the unique weak entropy solution to the IVP (1) with initial datum $\rho_{0}$ and the parameters $p_{0},s_{0},c^{\max}=c_{0}$.

Let $(U_{i},i\in\mathbb{N})$ be a sequence of independent and identically distributed (i.i.d) random variables on some probability space $(\Omega,\mathcal{A},P)$ each having a uniform distribution on $[0,1]$. Furthermore, let $(\xi_{i},i\in\mathbb{N})$ be a sequence of i.i.d exponentially distributed random variables on the same probability space $(\Omega,\mathcal{A},P)$ and independent of $(U_{i},i\in\mathbb{N})$ and choose $t_{n}\in[0,T]$, $y_{n}\in E$. The following thinning algorithm produces the next jump time $T_{n+1}$ and corresponding post jump location $Y_{n+1}$.

One can show, see [15], that

[TABLE]

for $t\geq t_{n}$ and $B\in\mathcal{E}$.

We set $T_{0}=0$ and $Y_{0}=(p_{0},s_{0},c_{0},\rho_{0})\in E$ and apply the thinning algorithm iteratively. In every iteration we obtain a new upper bound $\bar{\lambda}$ on the rates, which might increase but stays finite for finitely many iterations. Let denote $((T_{n},Y_{n}),n\in\mathbb{N}_{0})$ the constructed jump times and post-jump locations, then we define the piecewise deterministic process (PDP) $(X(t),t\in[0,T])$ as

[TABLE]

Remark 2.3*.*

The total variation bound on the solution is quite pessimistic for reasonable initial datum. 2. 2.

The total variation bound can be very large in small time intervals and the Algorithm 1 can not be used efficiently to simulate the model. 3. 3.

We expect $X$ being a Markov process but standard results, see [21] can not be applied since $\operatorname{BV}$ is no Borel space and the existence of regular conditional distributions is not guaranteed.

Multiple accidents on roads. In order to implement multiple accidents in the model, we label accidents and extend the state space as follows:

• •

positions are now given by $\vec{p}\in\mathbb{R}^{\mathbb{N}}$,

• •

sizes of the accidents are $\vec{s}\in\mathbb{R}^{\mathbb{N}}$,

• •

capacity reductions $\vec{c}\in[0,1)^{\mathbb{N}}$

and set

[TABLE]

with the norm

[TABLE]

Let $\lambda_{A}>0$ be the rate of an accident and $\lambda_{R}>0$ be the rate of resolving an accident. We define $m(\vec{c})=\min\{i\colon c_{i}=0\}$ and $\pi_{i}(z,\vec{v})=(v_{1},\dots,v_{i-1},z,v_{i+1},\dots)\in\mathbb{R}^{\mathbb{N}}$. A natural choice for the jump distribution is then given by

[TABLE]

Here, $\mu^{pos}_{y}=\beta\mu^{F}_{y}+(1-\beta)\mu^{D}_{\rho}$, where $\mu^{F}_{y}(B)=\int_{B}\frac{1}{C_{F}}F^{\vec{p},\vec{s},\vec{c}}(x,\rho(x))dx$ and $F^{\vec{p},\vec{s},\vec{c}}(x,\rho)=c_{road}(x)f(\rho)\prod_{i\in\mathbb{N}}c_{a}(x,p_{i},s_{i},c_{i})$. The sum $N(\vec{c})=\sum_{i\in\mathbb{N}}\mathds{1}_{c_{i}>0}$ corresponds to the number of accidents and we see that $B\mapsto\eta(y,B)$ is a probability measure. Since $\pi_{i}$ and $m$ are measurable functions, the mapping $y\mapsto\eta(y,B)$ is measurable if again $y\mapsto\mu^{pos}_{y}$ is measurable, see Lemma 2.2. Since $\lambda_{A}$ corresponds to the rate of an accident, we choose again

[TABLE]

and

[TABLE]

The upper bound on the rate function is now given by

[TABLE]

where $a(x)=c_{road}(x)\prod_{i\in\mathbb{N}}c_{a}(x,p_{i},s_{i},c_{i})$, $y=(\vec{p},\vec{s},\vec{c},\rho(t))$ and $\rho(t)$ is the unique weak entropy solution to (1).

We explain the choice of (6) by the following example. We consider two accidents with capacity reduction $\vec{c}=(0.5,0,0.5,0\dots)$, i.e. $N(\vec{c})=2$ and $m(\vec{c})=2$. We set $B_{\vec{p}}=\mathbb{R}\times B_{p}\times\mathbb{R}\times\cdots\in\sigma(\mathbb{R}^{\mathbb{N}})$, $B_{\vec{s}}=\mathbb{R}\times B_{s}\times\mathbb{R}\times\cdots\in\sigma(\mathbb{R}^{\mathbb{N}})$ and $B_{\vec{c}}=B_{c_{1}}\times B_{c_{2}}\times B_{c_{3}}\times\mathbb{R}\times\cdots\in\sigma([0,1)^{\mathbb{N}})$. Then, we set $B=B_{\vec{p}}\times B_{\vec{s}}\times B_{\vec{c}}\times\operatorname{BV}(\mathbb{R})$ and obtain

[TABLE]

This implies that the probability of resolving the first accident and no new accident, i.e. $B_{c_{1}}=\{0\}$, $B_{c_{3}}=B_{c_{2}}=\emptyset$, is given by

[TABLE]

In the same manner we obtain the probability of having a new accident somewhere with some size and no repairs, i.e. $B_{c_{1}}=B_{c_{3}}=\emptyset$, $B_{c_{2}}=B_{p}=\mathbb{R}$ and $B_{s}=[0,1)$,

[TABLE]

Hence, if $\lambda_{A}=\lambda_{R}$, the probabilities are equal with value $\frac{1}{3}$.

3 Numerical treatment and computational results

The Cauchy problem (1) is numerically solved using the Lax-Friedrichs scheme with a temporal step size $\Delta t>0$ and a fixed relation $\frac{\Delta t}{\Delta x}$ such that the scheme converges to the weak entropy solution $\rho$ of the Cauchy problem, cf. Lemma 2.1. We denote by

[TABLE]

the cell means of the initial datum $\rho_{0}$ for $X_{i}=i\Delta x$ and $i\in\mathbb{Z}$.

Since the position, size and capacity reduction stays constant between the jumps, we define the discrete deterministic dynamics as

[TABLE]

where ${\rho}_{0}$ is a piecewise constant function on $[x_{i-\nicefrac{{1}}{{2}}},x_{i+\nicefrac{{1}}{{2}}})$ given by the cell means ${\rho}_{i}^{0}$. Further, ${\rho}(t)$ is the piecewise constant function given by the numerical scheme with step size $\Delta t$ and a possibly smaller last step size to reach exactly $t$.

Then, we then approximate $\mu^{F}_{y}$ by

[TABLE]

and

[TABLE]

Thanks to the piecewise constant cell averages, we enjoy an explicit representation of $D\rho^{+}$ as

[TABLE]

The discretized version of the rate function $\psi(y)$ is then given by

[TABLE]

In order to use Algorithm 1, we need a uniform upper bound on $\bar{\psi}$ which will depend on the number of accidents and grows exponentially due to the total variation bound in Lemma 2.1. In [24], less restrictive bounds have been used to define an appropriate algorithm but the bounds propsed will also depend on the exponential growth of the estimation of the total variation. We will introduce an approximate scheme, where the jump times are not simulated exactly in the following. The idea is based on the simulation algorithm introduced in [8], where an algorithm has been proposed to approximate a continuous-time Markov Chain.

The probability that an accident occurs at a time $T_{n+1}$, which is before $T_{n}+\Delta t$ is given by

[TABLE]

as $\Delta t\to 0$. This is true since $t\mapsto\psi(\phi_{t}(Y_{n}))$ is Lipschitz continuous by using Lemma 2.1 and the properties of $C_{F}$, i.e.

[TABLE]

Equation (7) motivates the following algorithm to approximate the next jump time $T^{a}_{n+1}$.

The parameters $\varrho\in(0,1]$, $\Delta t_{ref}>0$ are user-defined and $(U_{i},i\in\mathbb{N})$ is a sequence of i.i.d. uniformly distributed random variables. The parameter $\Delta t_{ref}$ allows to control the accuracy of the algorithm as the reference step size and $\varrho$ is the acceptance ratio in the case that $\Delta t_{ref}$ and $T$ are large. We see that Algorithm 2 uses an adaptive step size, where the adaptivity is incorporated by the current value of the rate function $\psi(y)$. We do not need any uniform bound, which is the obvious advantage and reduces the computational costs. Note that the exact solution operator $\phi$ has to be replaced by the discrete one in numerical implementations.

It remains to introduce the simulation procedure in the case that an accident happens or an accident does not cause capacity drop anymore, i.e. the simulation of $\eta(y,\cdot)$. The highest index $i$, where $c_{i}>0$ and $c_{i}=0$ corresponds exactly to $N(\vec{c})$ by construction if we start with $c_{j}>0$ for $j=1,\dots N(\vec{c})$ and $c_{j}=0$ for $j>N(\vec{c})$. One can use the well-known composition method, i.e. the distribution is a weighted sum of distributions, and we obtain the following procedure:

Choose whether an accident happens $Z_{1}=1$ or an accident is resolved $Z_{1}=0$ by a Bernoulli distributed random variable with $P(Z_{1}=1)=\frac{\lambda_{A}}{\lambda_{R}N(\vec{c})+\lambda_{A}}$. 2. 2.

• •

Case $Z_{1}=1$: Choose independently a position $p_{N(\vec{c})+1}$ according to the law $\mu_{y}^{pos}$, a size $s_{N(\vec{c})+1}$ according to $\mu^{size}$ and $c_{N(\vec{c})+1}\sim\mu^{cap}$ the corresponding capacity drop.

• •

Case $Z_{1}=0$: Choose a uniformly distributed index on $\{1,\ldots,N(\vec{c})\}$ to indicate which accident got removed.

Simulating the new position is straightforward since $i$ is picked according to

[TABLE]

and then the position within cell $i$ as a uniform distribution on $[x_{i-\nicefrac{{1}}{{2}}},x_{i+\nicefrac{{1}}{{2}}})$.

3.1 Simulation results

We assume a bounded road $[-L,L]\subset\mathbb{R}$ in the following with periodic boundary conditions $\rho(-L,t)=\rho(L,t)$ for (1) to avoid difficulties with boundary treatment. We assume possibly different road capacities on $[-L,L]$, i.e. let

[TABLE]

for $-L=x_{0}<x_{1}<\cdots x_{M}=L$ with $c_{m}\geq\underline{c}$ for $m=0,\dots M-1$ and $c_{0}=c_{M-1}$. The latter condition avoids a discontinuity for the periodic boundary conditions and implies that cars leaving at $x=L$ enter in the same manner at $x=-L$ again. Since we need enough regularity on $c_{road}$ to apply the total variation bound on the solution of (1), we use a mollifier $M_{\epsilon}$ with support $[-\epsilon,\epsilon]$ and $\int_{\mathbb{R}}M_{\epsilon}(x)dx=1$. Then, $c_{road}(x)=\tilde{c}_{road}\ast M_{\epsilon}(x)=\int_{\mathbb{R}}\tilde{c}_{road}(y)M_{\epsilon}(x-y)dy\in C^{\infty}$ and $|\operatorname{supp}(c_{road}^{\prime\prime})|\leq 2\epsilon M$.

We use the same ideas for the capacity reduction and define $\tilde{c}_{a}(x,p,s,c)=1-c\mathds{1}_{(p-\frac{s}{2},p+\frac{s}{2})}(x)$ for $p\in[-L,L]$, $s\in(-L,L)$ and $c\in[0,1-c_{min}]$. By defining $c_{a}(x,p,s,c)=\tilde{c}_{a}(\cdot,p,s,c)\ast M_{\epsilon}(x)$, and using $a(x)=c_{road}(x)\prod_{i\in\mathbb{N}}c_{a}(x,p_{i},s_{i},c_{i})$, we deduce

[TABLE]

as required.

Remark 3.1*.*

We face only finitely many accidents $P$-a.s. such that the infinite product in $a$ can be represented by a finite product. Therefore, the differentiation of $a$ can be understood in the classical sense.

The first example is devoted to the understanding of the dynamics of the LWR model with accidents derived in the previous sections. We are interested whether the modeling ideas can be also observed in computational experiments. The data we use is as follows: a time horizon $T=60$, a spatial discretization $\Delta x=\frac{1}{50}$ of $[-10,10]$ and $\Delta t_{ref}=\frac{1}{20}$. The initial density is chosen constant as $\rho_{0}(x)=0.4$ and the LWR flux is given by $f(\rho)=\rho(1-\rho)$. We assume a road capacity given by the non-smooth version as

[TABLE]

which implies a capacity reduction on $[0,5]$ caused by e.g. roads under constructions. To incorporate capacity drops caused by accidents, we use the function

[TABLE]

In numerical investigations, we have recovered that smoothing the latter functions does not significantly change the results for a fixed spatial step size $\Delta x$ and $\epsilon<\frac{\Delta x}{2}$, which reduces the computational costs significantly. For the stochastic part, we use $\lambda_{R}=\frac{1}{2}$, $\lambda_{D}=\frac{1}{10}$, $\lambda_{F}=\frac{1}{105}$ and assume

[TABLE]

as well as $\varrho=1$.

A first insight into the behavior of the model. Having all the parameters at hand, except $\beta$ from equation (4), we can get first insights into the behavior of the model using numerical simulations for varying $\beta.$ The latter parameter describes the influence of the current flux on the position of possible accidents, see (3).

Figure 1 shows the traffic density (black bold line) for different points in time and using only the information of $D\rho^{+}$ to determine the position of an accident, i.e. $\beta=0$. The rectangles in the figures indicate the range of the road affected by an accident, where a bright color corresponds to a capacity drop of $0.99$ and the other color of $0.5$, see $\mu^{c}$ in (8). Since the initial distribution is constant with a value of $0.4$, we draw the density at the first time at which an accident happens in Figure 1(a). Due to a spatial inhomogeneous road capacity $c_{road}(x)$, the initial density profile changed to a non-constant equilibrium traffic density. As we would expect, the accident happens at the incresaing part of the density, i.e. at the end of the traffic jam, with a road capacity reduction of $0.99$. At this position a traffic jam occurs until the accident is removed, see Figure 1(b). The traffic density relaxes to an equilibrium density again and the second accident happens at the end of the traffic jam as Figure 1(c) indicates. Again a capacity reduction of $0.99$ has been randomly chosen and a third accident occurs right after the second accident. The latter can be seen in Figure 1(d), which shows the traffic density at the time, where the second accident gets resolved.

At the time, where both accidents are resolved, see Figure 1(e), we see the high impact of the previous accidents on the density, which does not reach the equilibrium state until the next accident occurs as Figure 1(f) shows. Again, the position of the accident is at an increasing part of the density. Altogether, we see that our model is able to map the ideas of accidents at places with an increasing density and the numerical solutions look very confident using the CFL condition with equality.

In the following, we discuss simulation results using the parameter $\beta=0.5$ shown in Figure 2. We face an approximately equilibrium density at the time of a first accident again, see Figure 2(a). Here, the accident occurs close to the position zero, which is not an increasing part of the density. The accident is therefore created by the flux, which is uniform on the interval [-10,10] while the density is close to equilibrium.

As Figure 2(b) shows, the second accident happens at the traffic jam end. After the first accident has been removed, a third accident occurs and Figure 2(c) shows the traffic density at the time right before the fourth accident occurs. The fourth accident is inside the area of the second accident and has a small size of impact, see Figure 2(d). The latter accident occurred at this position since the flux around $\rho=0.5$ is the most highest and we are not in a stationary state.

Numerical verification of the approximate scheme. In order to verify numerically that the approximate algorithm works well, we study the distribution of the first jump time, i.e. the first time of an accident. Formula (5) exactly describes the cumulative distribution function (CDF), which can be approximated using the Lax-Friedrichs scheme to approximate $\phi$. Using the left-sided rectangular rule to approximate $\int_{0}^{t}\psi(\phi_{\tau}(y_{0}))d\tau$ and the Matlab function ecdf***Documentation: https://de.mathworks.com/help/stats/ecdf.html, 2019 to compute the empirical cumulative distribution function (ECDF) yields the results shown in Figure 3 computed by using $10^{4}$ samples of the first accident time $T_{1}^{a}$. First of all, we observe a very good fitting of the CDF by the ECDF computed with the approximation Algorithm 2. This implies that the corresponding probability distributions are close (in the weak sense). Furthermore, we observe that the parameter $\beta$ has no significant influence on the shape or values of the CDF as Figures 3(a) and 3(b) show.

In order to compare a histogram generated by the approximation procedure with the exact probability density function (pdf) $g(t)$, we can differentiate (5) and obtain

[TABLE]

Figure 4 shows a histogram of samples of $T_{1}^{a}$ and the theoretical result $g(t)$. We observe a good agreement between both quantities again, also independent of the choice of $\beta$.

Finally, we discuss the distribution of the first accident’s position. Figure 4 shows the histogram of samples of the first accident’s position, where we distinguish the cases $\beta=0$ and $\beta=0.5$ again. In both cases, the probability having an accident at position $x=-4$ is the most highest, which corresponds to the congestion end in the stationary traffic profile, see Figure 1(c) for example. One significant difference between $\beta=0$ and $\beta=0.5$ can be observed for $x\in[0,5]$, where in the case of $\beta=0$, i.e. no flux information, no accident happens.

In contrast, for $\beta=0.5$, there is a strictly positive probability having an accident in this interval, which is clear since the stationary value of $\rho$ is approximately at the maximal flow, i.e. at $0.5$.

To conclude, the numerical simulations inherit the ideas for the stochastic traffic flow model and the numerical results are convincing.

4 Conclusion

We successfully have derived a stochastic traffic flow model capturing random traffic accidents. Furthermore, a tailored numerical approximation scheme has been introduced, which also has been validated in numerical simulation examples.

The stochastic traffic flow model allows for road capacity planning and controlling variable speed limit systems in such a way that traffic accidents are rarely events, which might be future research. Additionally, the extension to a second order traffic models and networks can be considered.

Acknowledgments

This work was supported by the BMBF project ENets (05M18VMA) and the DFG project GO1920/10-1.