Piecewise deterministic Markov processes driven by scalar conservation laws

TL;DR

This paper studies PDMPs driven by scalar conservation laws, establishing existence and bounded variation of sample paths, and applies the framework to production and traffic models with data-driven flux variations.

Contribution

It introduces a novel PDMP framework driven by scalar conservation laws and demonstrates its applicability to real-world production and traffic systems.

Findings

Existence of PDMPs under certain conditions

Bounded variation estimates for sample paths

Effective modeling of flux scattering in data

Abstract

We investigate piecewise deterministic Markov processes (PDMP), where the deterministic dynamics follows a scalar conservation law and random jumps in the system are characterized by changes in the flux function. We show under which assumptions we can guarantee the existence of a PDMP and conclude bounded variation estimates for sample paths. Finally, we apply this dynamics to a production and traffic model and use this framework to incorporate the well-known scattering of flux functions observed in data sets.

Equations

$$u_t(x,t)+f(u(x,t){)}_x=0,u(x,0)=u_0(x),$$

$$u_t(x,t)+f^{\alpha (t)}(u(x,t){)}_x=0,u(x,0)=u_0(x),$$

$$\Vert u(\cdot ,t){\Vert }_{\infty }$$

$$\mathrm{TV}(u(\cdot ,t))$$

$$\Vert u(\cdot ,t)-u(\cdot ,s){\Vert }_{L^1}$$

$$\Vert u(\cdot ,t)-v(\cdot ,t){\Vert }_{L^1}\le \Vert u_0-v_0{\Vert }_{L^1}.$$

$$S_{s+t}^{\alpha }$$

$$S_0^{\alpha }$$

$$t\mapsto S_t^{\alpha }$$

$$\Vert S_t^{\alpha }u-S_t^{\alpha }v{\Vert }_{L^1}$$

$$t\mapsto S_t^{\alpha }{u}_0\text{ is the }$$

$${\varphi }_t:E\to E,\left(\begin{array}{c}u\\ \alpha \end{array}\right)\mapsto \left(\begin{array}{c}{S}_t^{\alpha }u,\\ \alpha \end{array}\right).$$

$$\Vert (s,u,\alpha )-(t,v,\beta )\Vert =|s-t|+\Vert u-v{\Vert }_{L^1}+|\alpha -\beta |.$$

$$\Vert (S_s^{\alpha }u,\alpha )-(S_t^{\beta }v,\beta )\Vert =\Vert S_s^{\alpha }u-S_t^{\beta }v{\Vert }_{L^1}+|\alpha -\beta |.$$

$$\Vert S_s^{\alpha }u-S_t^{\beta }v{\Vert }_{L^1}$$

$$\Vert S_t^{\alpha }v-S_t^{\beta }v{\Vert }_{L^1}$$

$$\le 2\Vert v-v_n{\Vert }_{L^1}+\Vert S_t^{\alpha }{v}_n-S_t^{\beta }{v}_n{\Vert }_{L^1}$$

$$\le 2\Vert v-v_n{\Vert }_{L^1}+t\Vert f^{\alpha }-f^{\beta }{\Vert }_{C^{0,1}}\mathrm{TV}(v_n),$$

$$\Vert S_s^{\alpha }u-S_t^{\beta }v{\Vert }_{L^1}$$

$$+2\Vert v-v_n{\Vert }_{L^1}+T\Vert f^{\alpha }-f^{\beta }{\Vert }_{C^{0,1}}\mathrm{TV}(v_n).$$

$$\Vert S_s^{\alpha }u-S_t^{\alpha }u{\Vert }_{L^1}<\frac{ϵ}{6},\Vert u-v{\Vert }_{L^1}<\frac{ϵ}{6},\Vert f^{\alpha }-f^{\beta }{\Vert }_{C^{0,1}}<\frac{ϵ}{\mathrm{TV}(v_n)6T},|\alpha -\beta |<\frac{ϵ}{6}$$

$$\Vert (S_s^{\alpha }u,\alpha )-(S_t^{\beta }v,\beta )\Vert <ϵ$$

$$\lambda :I\times \mathcal{B}(I)\times [0,T]\times L^1(\mathbb{R})\to {\mathbb{R}}_{>0},$$

$${\psi }_t(y)$$

$${\eta }_t(y,B)$$

$$X(t)={\varphi }_{t-T_n}(Y_n)\iff t\in [T_n,T_{n+1}),$$

$$\mathrm{TV}(u(t))=\mathrm{TV}(S_{t-T_1(\omega )}^{\alpha (T_1(\omega ),\omega )}u(t,\omega ))\le \mathrm{TV}(u(T_1(\omega ),\omega ))\le \mathrm{TV}(u_0)$$

$$\mathrm{TV}(u(t,\omega ))\le \mathrm{TV}(u_0).$$

$$f(\rho )=\min \{v\rho ,\mu \}$$

$$f^{\alpha }(\rho )=\mu (\alpha )(1-e^{-\frac{v(\alpha )}{\mu (\alpha )}\rho })$$

$$\lambda (\alpha ,B,t,\rho )=\bar{\lambda }(\rho ){\int }_B\frac{1}{\sqrt{2\pi {\sigma }^2(\rho )}}{e}^{-\frac{(z-\bar{\alpha }{)}^2}{2{\sigma }^2(\rho )}}dz$$

$$f^{\alpha }(\rho )=\frac{\theta -1}{{\alpha }^{\theta }}\frac{1}{\Gamma (\frac{\theta -1}{\alpha })}{\rho }^{\theta -1}{e}^{-\frac{\theta -1}{\alpha }\rho }$$

$$\lambda (\alpha ,B,t,\rho )=\bar{\lambda }(\alpha ,\rho ){\int }_B\frac{1}{2a(\alpha ,\rho )}\text{mathds}{1}_{[{\alpha }_0-a(\alpha ,\rho ),{\alpha }_0+a(\alpha ,\rho )]}(z)dz$$

Full text

{AMS}

60J25, 35L65

1 Introduction

The simplicity of scalar conservation laws allows to understand general behaviors of underlying models but, on the other hand, they are based on qualified assumptions as for example steady state or expected values. One possibility to widen this class of models are systems of conservation laws, where fluctuations and higher order moments can be governed. Another possibility to extend scalar conservation laws are stochastic effects. More precisely, starting from deterministic scalar conservation laws and a corresponding initial value problem (IVP)

[TABLE]

a natural extension is the incorporation of uncertainties. There already exist extensions based on a reformulation as stochastic differential equation like in [12] and partial stochastic differential equation as in [5, 17] in the literature. Also uncertain initial data as for example in [6] and random chosen flux functions [18] have been considered. In the latter work, the flux function is random and does not change randomly in time.

In contrast to [18], our goal is a stochastic process, which “chooses” a new flux function at random times, where these times and the random choice of the next flux function may dependent on the actual solution of the whole system. This can be easily motivated by, e.g. production models with machine failures [8, 10], and also opinion formation, change of state (gas to liquid or vice versa) are reasonable applications. This idea directly transfers us into the theory of piecewise deterministic Markov processes, see [14]. In detail, given a parametrized family of Lipschitz continuous flux functions $f^{\alpha}\in C^{0,1}(\mathbb{R})$ for $\alpha\in I\subset\mathbb{R}$, we are interested in a “solution” to

[TABLE]

where $\alpha(t)\in I$ denotes the current and random chosen flux function at time $t\in[0,T]$ and $x\in\mathbb{R}$.

We define how (2) has to be understood and how $\alpha(t)$ is specified in the subsequent section 2. This section is followed by applications and numerical results in the case of a production and traffic model in section 3.

2 Modeling Equations

Let $u\colon\mathbb{R}\to\mathbb{R}$ be a function, then we denote by $\operatorname{TV}(u)$ its total variation and define $\operatorname{BV}(\mathbb{R})=\{u\colon\mathbb{R}\to\mathbb{R}\colon\operatorname{TV}(u)<\infty\}$ as the set of all functions from $\mathbb{R}$ to $\mathbb{R}$ with total bounded variation, see, e.g. [13]. With this notation, it is well known as a result of Krusckov, see [13], that the IVP (1) has a unique weak entropy solution if $u_{0}\in\operatorname{BV}(\mathbb{R})\cap L^{1}(\mathbb{R})$ and if $f\in C^{0,1}(\mathbb{R})$, i.e., is Lipschitz continuous. Furthermore, the solution $u$ satisfies

[TABLE]

and is $L^{1}$ stable with respect to initial data

[TABLE]

Deterministic dynamics between jump times

In this section, we define the dynamics to (2) based on theory of PDMPs. Unfortunately, we cannot apply the theory of PDMPs directly on solutions to (1) with corresponding flux functions $f^{\alpha}$ since $\operatorname{BV}(\mathbb{R})$ is not separable and hence no Borel space. Following [2], we use the extended solution operator to (1) on $L^{1}(\mathbb{R})$ and denote it by $S_{t}^{\alpha}\colon L^{1}(\mathbb{R})\to L^{1}(\mathbb{R})$, where $\alpha$ indicates that flux function $f^{\alpha}$ is used. We directly deduce the following properties of the family $(S^{\alpha}_{t},t\in[0,T])$ for every $\alpha\in I$:

[TABLE]

Up to now, we have no specification of $\alpha(t)$. We define the state space $E=L^{1}(\mathbb{R})\times I$ equipped with the Borel $\sigma$-algebra $\mathcal{E}$ generated by the open sets induced by $\|(u,\alpha)\|=\|u\|_{L^{1}}+|\alpha|$ for $(u,\alpha)\in E$. Then $(E,\mathcal{E})$ is a Borel space.

Our aim is to switch the flux function only at random times, which results in deterministic dynamics between the jumps in the form of

[TABLE]

Properties (7)-(10) of $S$ directly translate to $\phi$. If we can show that $\phi\colon[0,T]\times E\to E$ is measurable, the dynamics $\phi$ is a candidate for deterministic dynamics in between jump times of a PDMP, see [14]. The following lemma 2.1 tells us a sufficient condition to prove measurability of $\phi$.

Lemma 2.1**.**

*Let the mapping $\alpha\mapsto f^{\alpha}$ from $I\to C^{0,1}(\mathbb{R})$ be continuous with $I\subset\mathbb{R}$ an interval, then $(t,u,\alpha)\mapsto(S^{\alpha}_{t}u,\alpha)$ is continuous from $[0,T]\times L^{1}(\mathbb{R})\times I\to L^{1}(\mathbb{R})\times I$ and consequently measurable. *

Proof 2.2**.**

Let $(s,u,\alpha),(t,v,\beta)\in[0,T]\times L^{1}(\mathbb{R})\times I$, then we use the norm

[TABLE]

According to this norm, we use

[TABLE]

To show continuity, we estimate $\|S_{s}^{\alpha}u-S_{t}^{\beta}v\|_{L^{1}}$ as follows:

[TABLE]

and conclude that we can make $\|S_{s}^{\alpha}u-S_{t}^{\alpha}u\|_{L^{1}}$ and $\|S_{t}^{\alpha}u-S_{t}^{\alpha}v\|_{L^{1}}$ sufficiently small by shrinking $\|(S_{s}^{\alpha}u,\alpha)-(S_{t}^{\beta}v,\beta)\|$ due to properties (9)-(10). Let $(v_{n},n\in\mathbb{N})$ be a sequence in $\operatorname{BV}(\mathbb{R})\cap L^{1}(\mathbb{R})$ satisfying $\lim_{n\to\infty}\|v_{n}-v\|_{L^{1}}=0$. We estimate $\|S_{t}^{\alpha}v-S_{t}^{\beta}v\|_{L^{1}}$ as follows:

[TABLE]

where we used the result from [11, p. 53] in the last estimate. Altogether, we find that

[TABLE]

Now, let $(t,v,\beta)\in[0,T]\times L^{1}(\mathbb{R})\times I$, $\epsilon>0$ and choose $n\in\mathbb{N}$ such that $\|v-v_{n}\|_{L^{1}}<\frac{\epsilon}{6}$ as well as $\delta>0$ such that

[TABLE]

implying

[TABLE]

*for all $(s,u,\alpha)\in[0,T]\times L^{1}(\mathbb{R})\times I$ satisfying $\|(s,u,\alpha)-(t,v,\beta)\|<\delta$. *

One simple example for a family of flux functions, which satisfies the continuity with respect to the parameter $\alpha\in I$ is given by $f^{\alpha}=\alpha f$ for $f\in C^{0,1}(\mathbb{R})$. Then $\|f^{\alpha}-f^{\beta}\|_{C^{0,1}}=\|f\|_{C^{0,1}}|\alpha-\beta|$.

Jump and jump time distributions

Following [14], we specify the transition intensities $q_{t}(y,B)\geq 0$, i.e., the rate to jump from $y\in E$ in a state in $B\in\mathcal{E}$ at time $t\in[0,T]$. This can be decomposed into $q_{t}(y,B)=\eta_{t}(y,B)\psi_{t}(y)$, where $\psi_{t}(y)$ is the total intensity that a jump occurs a time $t$ and $\eta_{t}(y,B)$ is the probability of a jump from $y$ into a state in $B$ provided a jump occurs at time $t$.

In order to use these intensities, we assume $(y,t)\mapsto\psi_{t}(y)$ to be measurable and for all $(y,t)$ we need $\int_{t}^{t+h}\psi_{s}(y)ds<\infty$ for $h=h(y,t)$ sufficiently small. For all $t$ we additionally assume that $\eta_{t}$ is a Marovian kernel, see, e.g. [1], for a definition. A further and natural assumption is that $\eta_{t}(y,\{y\})=0$ holds for all $(y,t)\in E\times[0,T]$.

At this point almost everything can happen at jump times but we fix the specific idea that the flux function only changes at the jump times. In detail, there is no jump in the solution of the conservation law component to inherit mass conservation again. To do so, we restrict on rates

[TABLE]

satisfying

$\sup\{\lambda(\alpha,I,t,u)\colon\alpha\in I,\;t\in[0,T],\;u\in L^{1}(\mathbb{R})\}\leq\lambda^{\text{max}}<\infty$, 2. 2.

for every $t\in[0,T]$, $(\alpha,u)\in E$ the mapping $B\mapsto\lambda(\alpha,B,t,u)$ is a measure, 3. 3.

for every $t\in[0,T]$, $B\in\mathcal{B}(I)$ the mapping $(\alpha,u)\mapsto\lambda(\alpha,B,t,u)$ is measurable, 4. 4.

for every $t\in[0,T]$, $(\alpha,u)\in E$ we have $\lambda(\alpha,\{\alpha\},t,u)=0$.

Then we define for every $y=(\alpha,u)\in E$ and $B\in\mathcal{E}$ the total intensity and jump distribution by

[TABLE]

Existence

Due to the uniform bound on $\psi_{t}$, we can use a so-called thinning algorithm to build the jump times $T_{n}$ and after jump locations $Y_{n}$ for $n\in\mathbb{N}_{0}$ iteratively, see [10, 15]. Since the number of jumps is finite $P$-almost surely, again due to the uniform bound on the rates, we obtain a stable random counting measure and theorem 7.3.1 from [14] can be applied. We obtain the following result

Theorem 2.3**.**

For every initial data $x_{0}=(\alpha_{0},u_{0})\in E$ there exists a stochastic process $X=(X(t),t\in[0,T])$ on some probability space $(\Omega,\mathcal{A},P)$, which satisfies

$X(0)=x_{0}$, 2. 2.

$X$* is a Markov process with respect to its natural filtration $\mathcal{F}^{X}=(\mathcal{F}^{X}_{t},t\in[0,T])$ given by $\mathcal{F}^{X}_{t}=\sigma(X(s),0\leq s\leq t)$,* 3. 3.

$X$* is piecewise deterministic and piecewise continuous, i.e., there exist jump times $T_{n}\in[0,T]$ and post jump locations $Y_{n}\in E$ for $n\in\mathbb{N}_{0}$ such that*

[TABLE]

where for convenience $T_{0}=0$ and $Y_{0}=x_{0}$.

Total Variation bounds and BV solutions

The extension of the solution to $L^{1}$ allowed us to use classical results from the theory of piecewise deterministic Markov processes to obtain the existence of a stochastic process, which satisfies our requirements. We expect that if the initial condition $u_{0}\in L^{1}(\mathbb{R})\cap\operatorname{BV}(\mathbb{R})$, then we deduce $u(t)\in L^{1}(\mathbb{R})\cap\operatorname{BV}(\mathbb{R})$ again as the following lemma shows.

Lemma 2.4**.**

*Let $X=(X(t),t\in[0,T])$ be the stochastic process from theorem 2.3 with $X(t)=(\alpha(t),u(t))\in E$. If $u(0)=u_{0}\in L^{1}(\mathbb{R})\cap\operatorname{BV}(\mathbb{R})$, then $u(t)\in L^{1}(\mathbb{R})\cap\operatorname{BV}(\mathbb{R})$ and $\operatorname{TV}(u(t))\leq\operatorname{TV}(u_{0})$. *

Proof 2.5**.**

Let $\omega\in\Omega$, $T_{n}(\omega)$ the jump times and $Y_{n}(\omega)$ the post jump locations of $X(\omega)$ for $n\in\mathbb{N}_{0}$. For $t\in[0,T_{1}(\omega))$ we have $\operatorname{TV}(u(t,\omega))=\operatorname{TV}(S^{\alpha_{0}}_{t}u_{0})\leq\operatorname{TV}(u_{0})$ by classical results on scalar conservation laws, see, e.g. [13]. At time $t=T_{1}$ the flux function changes and for $t\in[T_{1}(\omega),T_{2}(\omega))$ it follows

[TABLE]

by continuity of $t\mapsto u(t,\omega)$. Iteratively, we deduce

[TABLE]

Remark 2.6**.**

*Lemma 2.4 is only valid because we have no jumps in the $u$ component at jump times by construction. Using the same arguments, the mass in the $u$ component is preserved. *

3 Applications and numerical results

Since we motivated PDMPs driven by scalar conservation law dynamics by the scattering of real data, we discuss simulation results of two examples in this section. The first example is a production and the second example is a traffic flow model.

Production model

Macroscopic production models have been widely studied in the literature, see [3] for an overview. Since in production capacity drops occur due to machine failures or human influences, deterministic models have been extended to stochastic production models, see [4, 8, 9, 10]. Therein, a random flux function in the form of

[TABLE]

has been chosen with a deterministic production velocity $v>0$, a stochastic capacity $\mu$ for a production density $\rho$. The latter corresponds to the variable $u$ in our context. In [4, 9] the capacity $\mu$ is a Continuous Time Markov Chain, in [8] a semi-Markov process and in [10] a PDMP construction has been developed.

In contrast to the mentioned works, we consider a single production step instead of a network and use our more general setting that allows for further flux functions motivated by data sets, see e.g. [7]. One possible choice is

[TABLE]

for a continuous bounded capacity $\mu>0$ and velocity $v\geq 0$. Some calculation shows $\|f^{\alpha}-f^{\beta}\|_{C^{0,1}}=\mathrm{O}(|v(\alpha)-v(\beta)|+|\frac{v(\alpha)}{\mu(\alpha)}-\frac{v(\beta)}{\mu(\beta)}|)$ and the flux function fulfills the requirements to obtain the existence of a suitable stochastic process $X$, see theorem 2.3.

In Figure 1(a) flux functions for $\mu(\alpha)=1+\tanh(\frac{\alpha}{2})$ and $v(\alpha)=1+\tanh(\alpha)$ and different $\alpha$ are drawn. So, we can capture different production velocities and capacities by varying $\alpha$.

It remains to introduce jump rates $\lambda$ in the production setting. We want the total jump intensity to be dependent on the Work In Progress (WIP) on some interval $[a,b]\subset\mathbb{R}$, which is defined as $\text{WIP}(\rho(t))=\int_{a}^{b}\rho(x,t)dx.$ In detail, we assume as WIP increases, the probability of a change of the flux function increases and vice versa. The distribution of the post jump location is assumed to be symmetrical around $\bar{\alpha}\in\mathbb{R}$ with variance $\sigma^{2}>0$ and we exemplary use

[TABLE]

for every $\alpha\in\mathbb{R}$, $B\in\mathcal{B}(\mathbb{R})$, $t\in[0,T]$ and $\rho\in L^{1}(\mathbb{R})$. One reasonable choice for $\bar{\lambda}(\rho)$ is $\bar{\lambda}(\rho)=\lambda_{0}(1-e^{-\lambda_{1}\text{WIP}(\rho)})$ for some $\lambda_{0},\lambda_{1}>0$. For the subsequent simulation results, we assume $a=0$, $b=1$, $\lambda_{0}=5$, $\lambda_{1}=1$, $\sigma^{2}=10^{-2}$, $\bar{\alpha}=0$. The time horizon is $T=50$ and the numerical spatial domain is taken as large that boundary conditions have no influence at $x=0$ on the solution. The deterministic dynamics is approximated by a Godunov scheme and in figure 1(b) we see the result of one sample of the density flux relation at position $x=0$ generated by the model with initial data $\rho(x,0)=\frac{3}{2}(\sin(x)+1)e^{-\frac{|x|}{100}}$. The black markers consider to the density and flux at times $t=0,0.2,\dots,50$ and the black solid line in figure 1(b) represents the flux function for $\alpha=0$. We observe in this stochastic macroscopic production model the typical scattering effect like it is the case for microscopic production models driven by discrete event simulations in [7].

Traffic flow model

The scattering effect in the density flux diagram obtained by real data, see, e.g. [20, 21], is a fundamental pattern and important for the development of second order, stochastic and phase transition traffic flow models. In the so-called free phase we observe small fluctuations and an almost linearly increasing flux with respect to the density. At a critical density, the flux decreases in the so-called congested phase. The critical density and congested phase are characterized by higher variances, i.e. sacttering effects in data. There exist already stochastic approaches like in [16, 19] and a comprehensive overview is given in [22]. We will show that the framework, which we introduced in section 2 is able to capture the scattering effects as well.

As family of flux functions, we use, motivated by the shape of the probability density function of the Gamma distribution,

[TABLE]

for some parameter $\theta\geq 1$, $\alpha>0$, $\rho\geq 0$ and $\Gamma$ the Gamma function. If $\theta\geq 2$, we also have $f^{\alpha}\in C^{0,1}(\mathbb{R}_{\geq 0})$ and the maximum is attained at $\rho^{\ast}=\alpha$. In figure 2(a), we see the shape of the flux function by varying $\alpha\in[0.3,0.5]$ and $\theta=2.1$. We set

[TABLE]

for every $\alpha>0$, $B\in\mathcal{B}(\mathbb{R}_{>0})$, $t\in[0,T]$ and $\rho\in L^{1}(\mathbb{R})$. Here, we choose $\bar{\lambda}(\alpha,\rho)=\lambda_{0}+(\lambda_{1}-\lambda_{0})V(\alpha,\rho)$ for $\lambda_{0}=3$ as the minimal and $\lambda_{1}=10$ as the maximal rate, $a(\alpha,\rho)=\sqrt{\frac{9}{2\cdot 10^{3}}(V(\alpha,\rho)+1)}$ with $V(\alpha,\rho)=\int_{0}^{1}\mathds{1}_{\rho(x)\geq\alpha}dx$ and $\alpha_{0}=0.4$. The functional $V(\alpha,\rho)$ describes the portion of $[0,1]$, which is above the actual critical density $\alpha$ and always lies in between zero and one. To study the free phase, we use an initial condition in the form of $\rho_{0}(x)=(0.05+0.4\max\{\sin(x),0\})e^{-\frac{|x|}{100}}$. A sample of the density flux relation at $x=0$ as well as at $x=1$ is shown in figure 2(b) given at the times $t=0,0.1,\dots,50$. We observe a low scattering as expected. Contrary, in figure 2(c) a sample with initial condition $\rho_{0}(x)=(0.4+\max\{\sin(x),0\})e^{-\frac{|x|}{100}}$, i.e. congested case, is shown resulting in high scattering. Finally, in figure 2(d) the time evolution of the density and flux at $x=0$ in the congested case is shown. The density is not severely affected by the variation in $\alpha$ compared to the flux.

4 Conclusions

We have successfully incorporated random flux functions for scalar conservation laws in the sense of PDMPs. Additionally, we derived a sufficient condition for an arbitrary family of Lipschitz continuous flux functions such that we can guarantee the existence of a PDMP. The motivation of scattering effects in macroscopic models has been recovered in numerical simulation results in the case of a production and traffic flow model.

To cover more complex dynamics, like space dependent flux functions, the theory can be extended in a suitable way as future research. This can be relevant to model traffic accidents and models, where spatial events can happen. Additionally, systems of conservation laws should be examined as deterministic dynamics for PDMPs since the extension to $L^{1}$ solutions is not straightforward anymore.