Modelling the spatial spread of COVID-19 in a German district using a diffusion model

1.   Introduction

Macroscopic traffic flow models based on hyperbolic conservation laws have been intensively investigated during the last decades, see ^[1,2] for an overview. The various research directions include theoretical and numerical investigations such for instance well-posedness ^[3], coupled models ^[4], network extensions ^[2,5], optimal control ^[6], or more recently, data-driven approaches ^[7] while stochastic traffic models have been less considered ^[8,9].

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 ^[10,11,12], queueing theory approaches ^[13,14] or kinetic models ^[15]. There are only a few contributions, where the presence of accidents is described by a stochastic process ^[13].

Therefore, the aim of this paper is to combine the stochastic modeling of accidents with the Lighthill-Whitham-Richards (LWR) model ^[16] 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 ^[17,18]. In ^[19], 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 ^[19], we face different challenges here: since accidents happen at various spatial positions, we use a space-dependent flux function capturing traffic accidents as spatial capacity drops in the LWR model. This is the deterministic part only and we need an appropriate model representing the stochasticity of traffic accidents. Hereby, we model the position, size and the capacity reduction caused by the accident. The modeling of random positions is done by using two relevant effects leading to traffic accidents: high flux and tailbacks. High fluxes imply that cars drive at an appropriate speed but are also comparably dense such that inattention leads to rear-end collision. The same holds true at tailbacks where we have a jump from low to high density, i.e., the derivative of the traffic density is strongly positive. From the mathematical point of view, we need a notion of derivative of the traffic density and we use functions of so-called bounded variation. Unluckily, the LWR model with space-dependent flux does not admit solutions which are of bounded variation in space. Therefore, we pose additional assumptions to the classical ones guaranteeing that the solution is of bounded variation. The sampling of the times at which accidents happen is also more involved because of a pessimistic bound on the traffic accident rate. We therefore introduce a numerical scheme to get rid-off this bound and show that it approximates the model in numerical examples.

There are different works about hyperbolic equation based dynamics connected to randomness as for example random velocity fields ^[20,21] and propagation of uncertainty ^[22]. 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

Before we start with the mathematical modeling of accidents, we motivate basic ideas first. The Lighthill-Whitham-Richards (LWR) model is described by the solution to the following hyperbolic conservation law

where $\rho(x, t) \in [0, \rho_{max}]$ is the traffic density at $x \in \mathbb{R}$ and $t > 0$ bounded by $\rho_{max} > 0$ and at time $t = 0$ an initial profile $\rho_0$ is prescribed. The function $v \colon [0, \rho_{max}] \to [0, v_{max}]$ is the velocities of cars, where $v(0) = v_{max} > 0$, $v(\rho_{max}) = 0$ and $v^\prime < 0$ imply that it decreases with increasing densities. It is also common to write $f(\rho) = \rho v(\rho)$ as the so-called LWR flux function if $v$ satisfies the assumptions mentioned before. If we fix the maximal density $\rho_{max}$ but choose a different flux function $\tilde{v} < v$, cars drive slower using $\tilde{v}$. Alternatively, $\rho\tilde{v}(\rho) < \rho v(\rho)$ means a lower flux and hence a lower capacity of a road. Using a space-dependent velocity $v(x, \rho)$, where $x \mapsto v(x, \rho)$ is a piecewise constant velocity function, i.e., on road segments are used different $v$, we can model capacity differences on roads caused by different number of lanes, speed limits, or, as we will use, accidents in the form of capacity drops. Without loss of generality, we assume $\rho_{max} = 1$ in the following.

2.1.   General setting

Let $f \colon [0, 1]\to [0, \infty)$ be a function of LWR type, e.g., $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_r\colon \mathbb{R} \to \mathbb{R}_{ > 0}$ and use $c_r(x)f(\rho)$ as space-dependent flux. An appropriate choice for $c_r$ 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_r(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

To illustrate the function $c_a$, we state the following examples satisfying the assumptions from before:

Example 2.1.

● Discontinuous: $c_a(x, p, s, c) = 1-c \mathbb{1}_{(p-s, p+s}(x)$

● Continuous: $c_a(x, p, s, c) = 1-c \mathbb{1}_{(p-s, p+s}(x)(1-\frac{(x-p)^2}{s^2})$

Altogether, we end up with the following Cauchy problem

In order to state existence and uniqueness results for (2.1), we introduce the notion of total variation $\operatorname{TV}(\rho)$ for a function $\rho$ which we define as

The space of functions with bounded total variation is denoted by $\operatorname{BV}(\mathbb{R}) = \{\rho \in L^1(\mathbb{R}) \colon \operatorname{TV}(\rho) < \infty\}$ and we call a function a BV function if it is a function in $\operatorname{BV}(\mathbb{R})$. The Cauchy problem (2.1) admits a unique entropy solution, see ^[23] if $\operatorname{TV}(c_a(\cdot, p, s, c)c_r(\cdot)) < \infty$, $\rho_0 \in \operatorname{BV}(\mathbb{R})$ and if $c_a(\cdot, p, s, c)c_r(\cdot)$ is differentiable with except of finitely many points. Additionally, we need that

The latter conditions arise from the technique used in ^[23] which is the so-called front tracking technique. The idea is to approximate the initial condition by a simple step function leading to finitely many Riemann problems. Using the singular mapping $\Psi$, which is a one-to-one mapping, the existence of a solution is shown in the same manner as it is the case with space-independent fluxes, see for example ^[24]. Due to $\Psi$, the existence of $\rho$ is also ensured. Within the proof, the total variation bounds are essential and require the assumptions on $\rho_0$ and $c_a(x, p, s, c)c_r(x)$. However, this only shows that $\operatorname{TV}(\Psi(\rho(\cdot, t))$ is finite for all $t \in [0, T]$ but 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 (2.1) remains in $\operatorname{BV}(\mathbb{R})$.

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

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 (2.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

Here, $\rho_i^{j}$ is the approximated value of $\rho(x_i, t^j)$, where $x_i = i \Delta x$, $t^j = j\Delta t$ for $i \in \mathbb{Z}$, $j \in \mathbb{N}_0$ are the spatial and temporal grids. The quotient $\lambda = \frac{\Delta t}{\Delta x}$ is an important constant used to prescribe the domain of dependence of the numerical scheme. Note that $\lambda$ is the reciprocal of a speed. To obtain a numerical reasonable scheme for the Cauchy problem (2.1), we have to ensure that the analytic domain of dependence is contained in the numerical domain of dependence since otherwise information gets lost. The largest velocity at which information may propagate in our model can be recovered from the slope of the characteristic curves and is given by $\sup_{x, \rho} |a(x)f^\prime(\rho)| \leq \|a\|_\infty\|f^\prime\|_\infty$. The distance information is then bounded by $\Delta t \|a\|_\infty\|f^\prime\|_\infty \leq \Delta x$, which is equivalent to

the so-called Courant-Friedrichs-Lewy (CFL) condition. The convergence of the Lax-Friedrichs scheme has been studied in ^[25], whereas in ^[26,27] 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.

● $L^\infty$ estimate. Using the CFL condition

we deduce that

The latter implies $\|\rho^j\|_\infty \leq \|\rho_0\|_\infty+T \|f\|_\infty \|a^\prime\|_\infty$.

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

Using the CFL condition and

yields

Hence, we have

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

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

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

We can directly estimate the first expression by

To analyze the second expression we start with

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

and therefore

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 ^[25], 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 ^[23]. 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

In the following, we introduce a LWR model which incorporates traffic accidents. For this reason, we start with a simple example to strengthen our motivation.

Example 2.3. In Figure 1(a), we observe a traffic situation given as the traffic density (solid line) and the corresponding velocity (dashed line) in the case of the classical LWR model with $v(\rho) = 1-\rho$. In order to model traffic accidents, we should be able to prescribe at which positions $x$ a traffic accident is more likely. If we choose the traffic density and interpret it as probability (modulo normalizing), we obtain a bad estimation since at positions, where cars are in a traffic jam ($\rho \approx 1$), we assume the highest probability of an accident. Exactly the opposite holds true if we choose the velocity as a measure of likeliness of an accident, where at places with no cars is the probability of an accident the highest. Therefore, we choose the flux as a representative of the probability of an accident, see Figure 1(b). This is reasonable since the flux describes the amount of cars traveling at this place.

A very common reason for accidents are significant changes in the traffic density, i.e., end of traffic jams or distinctive traffic waves. This is not captured by the flux and needs an additional modeling approach, which is discussed later on in this section.

The parameters to incorporate a traffic accident in Eq (2.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, see example 2.3. Also, we observe accidents at the end of traffic jams. To summarize, the following modeling ideas should be included:

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

2. a higher accident probability at places, where we observe an increase in the density (as for example tailbacks).

Regarding 1. The flow $F^{p, s, c}(x, \rho)$ exactly describes the combination of density, i.e., car distances, and space-dependent velocities such that at places with high capacity and where $\rho = \rho^\ast$ the probability of an accident can be assumed to be the highest. This idea corresponds to a probability capturing random accidents caused by human failures solely (i.e., excluding tailbacks). If the velocity function $v(\rho)$ is uniformly bounded for $\rho \in [0, 1]$, the normalizing constant

is finite and we can define the family of probability measures

for $\rho \in \operatorname{BV}(\mathbb{R})$ 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, see the motivating example 2.3.

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

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

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 ^{[28,29,30,31]}.

To clarify this approach, we consider the density situation given in Figure 1(b). In this case, we obtain

where $\epsilon_z$ denotes the Dirac measure with unit mass in $z$. We see directly that $\operatorname{TV}(\rho) = D\rho^+(\mathbb{R})+ D\rho^-(\mathbb{R})$. If we normalize $D \rho^+$ by $D \rho^+(\mathbb{R})$, we obtain a discrete probability distribution which support consists of the places, where a jump from low to high density is observed and the positions are related to a probability proportional to the height of the jump. This probability distribution incorporates both, the position and the fact that a higher speed difference from high to low determines the probability of an accident. The jump height in the density is directly correlated to the speed since $v^\prime(\rho) < 0$. Therefore, 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

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

Summarizing, we define

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)))$.

Remark 2.4. One could also include a dependence of the distributions $\mu^{size}$, $\mu^{cap}$ on $\rho$. We keep these distributions independent of $\rho$ here, since it is not obvious how a functional dependence may look like and if it is really the case in traffic flow.

In a natural way, we collect the details using the product space

with norm

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

where we recover that $\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.5. 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 by definition of $\eta(y, \cdot),$ the mapping $B \mapsto \eta(y, B) \geq 0$ is a measure. Since $E = \mathbb{R} \times \mathbb{R} \times [0, 1) \times \operatorname{BV}(\mathbb{R}) \in \mathcal{E}$ and by utilizing the definition of a product measure, we have

In the same manner, we have

since either $\mu^{size}(\{s\}) = 0$ for all $s$ or $\mu^{cap}(\{c\}) = 0$ for all $c \in [0, 1)$ by assumption.

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$. 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

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

We also have

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 ^[17,18]. Let $\psi(y) > 0$ be the traffic accident rate for a given state $y$ of the system. The normalization coefficient $C_F = C_F(\rho)$ maps the state of the system to a specific value. From the definition of $C_F$, we know that it is the total flux, which is in fact positive correlated to the accident probability, i.e., a higher flux implies a higher accident probability and vice versa. Analogously, $D\rho^+(\mathbb{R})$ is positive correlated to the accident probability since a higher number of tailbacks and higher increases in jumps increase the probability of an accident at tailbacks. Hence, a possible choice for a rate function $\psi \colon E \to (0, \infty)$ is given by

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 (2.1), then for $a(x) = c_a(x, p, s, c)c_r(x)$ it holds that

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.2.

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

where $\rho(t)$ is the unique weak entropy solution to the IVP (2.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 ^[19], that

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

Remark 2.6.

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

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. We expect $X$ being a Markov process but standard results, see ^[18] 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

with the norm

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

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_r(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}} \mathbb{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.5. Since $\lambda_A$ corresponds to the rate of an accident, we choose again

and

The upper bound on the rate function is now given by

where $a(x) = c_r(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 (2.1).

We explain the choice of (2.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

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

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)$,

Hence, if $\lambda_A = \lambda_R$, the probabilities are equal with value $\frac{1}{3}$.

3.   Numerical treatment and computational results

The Cauchy problem (2.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.2. We denote by

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

where ${\rho}_0$ is a piecewise constant function on $[x_{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}, x_{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 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 approximate $\mu^F_{y}$ by

and

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

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

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.2. In ^[32], 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 ^[33], 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

as $\Delta t \to 0$. This is true since $t \mapsto \psi(\phi_t(Y_n))$ is Lipschitz continuous by using Lemma 2.2 and the properties of $C_F$, i.e.,

Equation (3.1) 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:

1. 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. ● 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

and then the position within cell $i$ as a uniform distribution on $[x_{i - {1 \mathord{\left/ {\vphantom {1 2}} \right. } 2}}, x_{i + {1 \mathord{\left/ {\vphantom {1 2}} \right. } 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 (2.1) to avoid difficulties with boundary treatment. We assume possibly different road capacities on $[-L, L]$, i.e., let

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_r$ to apply the total variation bound on the solution of (2.1), we use a mollifier $M_\epsilon$ with support $[- \epsilon, \epsilon]$ and $\int_ \mathbb{R} M_ \epsilon(x) dx = 1$. Then, $c_r(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_r^{\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 \mathbb{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_r(x)\prod_{i \in \mathbb{N}}c_a(x, p_i, s_i, c_i)$, we deduce

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

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

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

as well as $\varrho = 1$.

A first insight into the behavior of the model. Having all the parameters at hand, except $\beta$ from Eq (2.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 (2.3).

Figure 6 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^{cap}$ in (3.2). 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 6(a). Due to a spatial inhomogeneous road capacity $c_r(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 6(b). The traffic density relaxes to an equilibrium density again and the second accident happens at the end of the traffic jam as Figure 6(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(f), 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(g), 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(h) 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.

Figure 3 shows a sample path of the model with $\beta = 1$, i.e., only flux information is used to determine the position of an accident. We observe in Figure 3(a)-3(b) an accident at densities, which are close to $\rho^\ast = 0.5$ here, which is exactly what we would expect as the flux $f$ without spatial capacity restrictions is maximal there. The third and fourth accident occur at places with the maximal road capacity, i.e., not in $[0, 5]$.

In summary, we recover the modeling ideas in the simulation results again. More precisely, for $\beta = 0$ we only have accidents at positions with an increasing density profile and in the case $\beta = 1$, the flux determines the position of the traffic accidents.

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 (2.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 $\texttt{ecdf}$^* to compute the empirical cumulative distribution function (ECDF) yields the results shown in Figure 5 computed by using $10^4$ samples of the first accident time $T_1^a$.

^* Documentation: https://de.mathworks.com/help/stats/ecdf.html, 2019

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 4(a) and 4(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 (2.5) and obtain $g(t) = \psi(\phi_t(y_0))e^{-\int_0^t \psi(\phi_\tau(y_0))d\tau}.$ Figure 5 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 5 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 6(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$.

Impact of boundary conditions. In the following, we pose a different type of boundary conditions, i.e., we assume an inflow condition at the left boundary and a homogeneous Neumann boundary condition at the right boundary. We have to keep in mind that the inflow $G_{in}(t)$ at the left boundary $x_l$ can be posed only at times where $f^\prime(\rho(x_l), t) > 0$, i.e., information can enter into the road segment. Otherwise, we also use homogeneous Neumann boundary conditions at $x_l$. We use the same parameters and functions as in the latter numerical investigations and discuss the influence of $G_{in}(t)$ on the model. If we have no accident, the road capacity is given by the bottleneck $\min_{x}c_r(x)f(\rho^\ast) = \frac{5}{4}$ in our setting.

Figure 6 shows a simulation result with an inflow $G_{in}(t) = \frac{15}{16} < \frac{5}{4}$ being below the capacity of the road segment. Since $\beta = 0$, we have accidents at tailbacks only. In Figure 7, an inflow $G_{in}(t) = \frac{5}{4}$ is assumed leading to more distinctive traffic jams than the inflow before. Accidents are consequently more likely at positions which are further to the left of the bottleneck compared to the previous case. In Figure 7(c), we have an accident right left close to the boundary such that the inflow has not been assigned due to $f^\prime \leq 0$ at this place. After this accident has been resolved, the inflow is again obtained as one can see in Figure 7(d).

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 derivation from a microscopic model and the extension to a second order traffic models and networks can be considered.

Acknowledgments

This work was supported by the BMBF project ENets (05M18VMA), the DFG project GO1920/10-1 and the DAAD project "Stochastic dynamics for complex networks and systems" (Project-ID 57444394).

Conflict of interest

All authors declare no conflicts of interest in this paper.