# Priority-based Riemann solver for traffic flow on networks

Maria Laura Delle Monache, Paola Goatin, Benedetto Piccoli

► To cite this version:

Maria Laura Delle Monache, Paola Goatin, Benedetto Piccoli. Priority-based Riemann solver for traffic flow on networks. Communications in Mathematical Sciences, 2018, 16 (1), pp.185-211. 10.4310/CMS.2018.v16.n1.a9. [https://doi.org/10.4310/CMS.2018.v16.n1.a9](https://doi.org/10.4310/CMS.2018.v16.n1.a9). [arXiv:1013.3682](https://arxiv.org/abs/1013.3682)

HAL Id: hal-01336823

https://inria.hal.science/hal-01336823v1

Submitted on 23 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Priority-based Riemann solver for traffic flow on networks*

Maria Laura Delle Monache† Paola Goatin‡

Benedetto Piccoli§

## Abstract

In this article we introduce a new Riemann solver for traffic flow on networks. The Priority Riemann solver (PRS) provides a solution at junctions by taking into consideration priorities for the incoming roads and maximization of through flux. We prove existence of solutions for the solver for junctions with up to two incoming and two outgoing roads and show numerically the comparison with previous Riemann solvers. Additionally, we introduce a second version of the solver that considers the priorities as softer constraints and illustrate numerically the differences between the two solvers.

**Keywords:** Scalar conservation laws, Traffic flow, Riemann solver

**AMS subject classifications:** 90B20, 35L65

## 1 Introduction

Conservation law on network is now a mature field with an increasing number of contributions in recent years. The theory for the scalar case is quite

<footnote>*This research was supported by the NSF grant CNS #1446715, by KI-Net "Kinetic description of emerging challenges in multiscale problems of natural sciences" - NSF grant # 1107444 and by the INRIA associated team Optimal REroute Strategies for Traffic management (ORESTE).</footnote>

<footnote>†Department of Mathematical Sciences, Rutgers University - Camden, Camden, NJ, USA (email: ml.dellemonache@rutgers.edu).</footnote>

<footnote>‡Inria Sophia Antipolis - Méditerranée, France (email: paola.goatin@inria.fr).</footnote>

<footnote>§Department of Mathematical Sciences and CCIB, Rutgers University - Camden, Camden, NJ, USA (email: piccoli@camden.rutgers.edu).</footnote>
developed (see [9, 15, 18]), with most results based on the concept of Riemann
solver. The latter is the network equivalent to the classical Riemann solvers
for conservation laws on the real line and provide a solution to Riemann
problems at junctions, i.e., Cauchy problems with constant initial data on
each road.

This theory was applied to different domains, including vehicular traffic [14],
supply chains [1], irrigation channels [3] and others. For a complete account
of recent results and references we refer the reader to the survey [6].

For vehicular traffic, authors considered many different traffic situations
to be modeled, thus proposing a rich set of alternative junction models even
for the scalar case, see [8, 9, 12, 13, 14, 18, 22, 24]. Here, we first propose a
new model which considers priorities among the incoming roads as the first
criterion and maximization of flux as the second. The main idea is that the
road with the highest priority will use the maximal flow taking into account
also outgoing roads constraints. If some room is left for additional flow then
the road with the second highest priority will use the left space and so on.
A precise definition of the new Riemann solver, called *Priority Riemann Solver*, is based on a traffic distribution matrix *A* (Definition 11), a priority
vector *P* = (*p*₁, ..., *p*ₙ) (with *p*ᵢ ≥ 0 and ∑ᵢ *p*ᵢ = 1) and requires a recursion
method, which is described in Algorithm 1.

We also model special situations in which some outgoing roads do not absorb
traffic from some incoming ones and propose an alternative solver with softer
priorities, see Algorithm 2.

During the writing of this manuscript we discovered that our priority-based
Riemann solver may be obtained as limit of solvers defined by Dynamic
Traffic Assignment based on junctions with queues [5].

The general existence theorem of [15] can be applied to every Riemann
solver satisfying three general properties, called (P1)-(P3), but can not be
applied in the present case. Indeed, the proof is based on estimates on the
flow total variation in space on the network in terms of the total variation in
time of the flow through the junction Γ, see definition (5). In turn the latter
is bounded thanks to the general property (P3), which ensures that waves
bringing flux decrease to the junction provoke a decrease of Γ. Such property
(P3) is not satisfied by the Priority Riemann Solver (see the Appendix: case
A2 with flux increase corresponding to Figure 6a).

Therefore, we achieve existence via a new set of general properties. Property
(P1) is the same as that of [15], while we modify (P2) and (P3) by using
estimates involving not only Γ but also the maximal flow along the priority
vector P in the set of admissible flows, see definition (13).

Then we apply the general theory to the Priority Riemann Solver by proving that the new (P1)-(P3) are satisfied for junctions with at most two incoming and two outgoing roads.

Then, to illustrate the Priority Riemann Solver, the one with soft priorities and compare with existing ones, we implement numerical simulations via the Godunov scheme.

The paper is organized as follows. In Section 2 we introduce the basic definitions of the theory of conservation laws on networks, then in Section 3 we define our Priority Riemann Solver and prove existence of solutions to Cauchy problems in Section 4. In Section 5, an alternative definition of the Riemann Solver with softer priorities is described and lastly, in Section 6, we propose a numerical discretization and show some numerical simulations comparing our Solvers to existing ones. The Appendix 7 collects the proof of the main theorem of the paper.

## 2 Basics

In this section we recall the basic definitions and results of the theory of conservation laws on networks, based on the concept of Riemann solver at junctions. Due to finite propagation speed of waves, to achieve existence results for Cauchy problems it is not restrictive to focus on a single junction. For details on how to extend the results to a general network, we refer the reader to [14, 15].

Fix a junction J with n incoming roads $I_1, \dots, I_n$ and m outgoing roads $I_{n+1}, \dots, I_{n+m}$, where $I_i =] -\infty, 0] \ (i \in \{1, \dots, n\})$ and $I_j = [0, +\infty[ \ (j \in \{n+1, \dots, n+m\})$. The traffic on each road $I_l$ ($l \in \{1, \dots, n+m\}$) is modeled using the celebrated Lighthill-Whitham-Richards model (briefly LWR, see [21, 23]):

$$ \partial_t \rho_l + \partial_x f(\rho_l) = 0, \quad t \ge 0, x \in I_l, \qquad (1) $$

where $\rho_l(t, x) \in [0, \rho_{\max}]$, is the car density, $v_l = v_l(\rho_l)$ is the average velocity and $f_l = f_l(\rho_l) = \rho_l v_l(\rho_l)$ is the flux. For simplicity, throughout the paper we assume $\rho_{\max} = 1$ and $f_l = f$ for all $l = 1, \dots, n+m$.

We make the following assumptions on the flux function f:

(H) $f: [0, 1] \to \mathbb{R}$ is a Lipschitz continuous and concave function satisfying

$$ 1. f(0) = f(1) = 0; $$
2. there exists a unique $\rho_{cr} \in ]0, 1[$ such that $f$ is strictly increasing in $[0, \rho_{cr}$[ and strictly decreasing in $] \rho_{cr}, 1]$.

As usual, entropic solutions and weak solutions at junctions are given by:

**Definition 2.1** A function $\rho_l \in C([0, +\infty[; L_{loc}^1(I_l))$ is an entropy-admissible solution to (1) in the arc $I_l$ if, for every $k \in [0, 1]$ and every $\varphi : [0, +\infty[ \times I_l \to \mathbb{R}$ smooth, positive and with compact support in $]0, +\infty[ \times (I_l \setminus \{0\})$, it holds

$$ \int_{0}^{+\infty} \int_{I_l} \left( |\rho_l - k| \partial_t \varphi + \operatorname{sgn}(\rho_l - k) (f(\rho_l) - f(k)) \partial_x \varphi \right) dxdt \geq 0. \quad (2) $$

**Definition 2.2** A collection of functions $\rho_l \in C([0, +\infty[; L_{loc}^1(I_l))$, ($l \in \{1, \dots, n+m\}$) is a weak solution at $J$ if

1. for every $l \in \{1, \dots, n+m\}$, the function $\rho_l$ is an entropy-admissible solution to (1) in the road $I_l$;

2. for every $l \in \{1, \dots, n+m\}$ and for a.e. $t > 0$, the function $x \mapsto \rho_l(t, x)$ has a version with bounded total variation;

3. for a.e. $t > 0$, it holds

$$ \sum_{i=1}^{n} f(\rho_i(t, 0-)) = \sum_{j=n+1}^{n+m} f(\rho_j(t, 0+)), \quad (3) $$

where $\rho_l$ stands for the version with bounded total variation of 2.

A Riemann problem at the junction $J$ is a Cauchy problem with constant initial data on each road. More precisely, given $\rho_{1,0}, \dots, \rho_{n+m,0} \in [0, 1]$, the corresponding Riemann problem at $J$ is given by

$$ \left\{ \begin{array}{l} \partial_t \rho_l + \partial_x f(\rho_l) = 0, \\ \rho_l(0, \cdot) = \rho_{0,l}, \end{array} \right. \quad l \in \{1, \dots, n+m\}. \quad (4) $$

For a collection of functions $\rho_l \in C([0, +\infty[; L_{loc}^1(I_l))$ ($l \in \{1, \dots, n+m\}$) such that, for every $l \in \{1, \dots, n+m\}$ and a.e. $t > 0$, the map $x \mapsto \rho_l(t, x)$ has a version with bounded total variation, we define the functionals

$$ \Gamma(t) := \sum_{i=1}^{n} f(\rho_i(t, 0-)) \quad (5) $$
and

$$TV_f(t) := \sum_{l=1}^{n+m} TV(f(\rho_l(t, \cdot))). \qquad (6)$$

Notice that $\Gamma$ is the flux through the junction, i.e. the total number of cars crossing the junction $J$ per unit of time, while $TV_f$ is the total variation of the flux on the whole network. From the flux bounds we easily derive:

$$0 \leq \Gamma(t) \leq nf(\rho_{cr}). \qquad (7)$$

A Riemann solver at $J$ is defined by:

**Definition 2.3** A Riemann solver $\mathcal{RS}$ is a function

$$\mathcal{RS}: \begin{array}{rcl} [0, 1]^{n+m} & \to & [0, 1]^{n+m} \\ (\rho_{1,0}, \dots, \rho_{n+m,0}) & \mapsto & (\bar{\rho}_{1}, \dots, \bar{\rho}_{n+m}) \end{array}$$

satisfying the following properties

$$1. \sum_{i=1}^{n} f(\bar{\rho}_i) = \sum_{j=n+1}^{n+m} f(\bar{\rho}_j);$$

2. for every $i \in \{1, \dots, n\}$, the classical Riemann problem

$$\left\{ \begin{array}{ll} \rho_t + f(\rho)_x = 0, & x \in \mathbb{R}, t > 0, \\ \rho(0, x) = \begin{cases} \rho_{i,0}, & \text{if } x < 0, \\ \bar{\rho}_i, & \text{if } x > 0, \end{cases} \end{array} \right.$$

is solved with waves with negative speed;

3. for every $j \in \{n+1, \dots, n+m\}$, the classical Riemann problem

$$\left\{ \begin{array}{ll} \rho_t + f(\rho)_x = 0, & x \in \mathbb{R}, t > 0, \\ \rho(0, x) = \begin{cases} \bar{\rho}_j, & \text{if } x < 0, \\ \rho_{j,0}, & \text{if } x > 0, \end{cases} \end{array} \right.$$

is solved with waves with positive speed.

Moreover, the Riemann solver $\mathcal{RS}$ must satisfy the consistency condition if

$$\mathcal{RS}(\mathcal{RS}(\rho_{1,0}, \dots, \rho_{n+m,0})) = \mathcal{RS}(\rho_{1,0}, \dots, \rho_{n+m,0})$$

for every $(\rho_{1,0}, \dots, \rho_{n+m,0}) \in [0, 1]^{n+m}$.
For future use, we now provide some definitions for the LWR model and for Riemann problems at junctions, for more details see [14].

**Definition 2.4** We say that $(\rho_{1,0}, \dots, \rho_{n+m,0})$ is an equilibrium for the Riemann solver $RS$ if

$$RS(\rho_{1,0}, \dots, \rho_{n+m,0}) = (\rho_{1,0}, \dots, \rho_{n+m,0}).$$

**Definition 2.5** We say that a datum $\rho_i \in [0, 1]$ in an incoming road is a good datum if $\rho_i \in [\rho_{cr}, 1]$ and a bad datum otherwise.

We say that a datum $\rho_j \in [0, 1]$ in an outgoing road is a good datum if $\rho_j \in [0, \rho_{cr}]$ and a bad datum otherwise.

We also define the following function:

**Definition 2.6** Let $\tau: [0, 1] \to [0, 1]$ be the map such that:

1. $f(\tau(\rho)) = f(\rho)$ for every $\rho \in [0, 1]$;

2. $\tau(\rho) \neq \rho$ for every $\rho \in [0, 1] \setminus \{\rho_{cr}\}$.

Clearly, the function $\tau$ is well defined and satisfies

$$0 \leq \rho \leq \rho_{cr} \iff \rho_{cr} \leq \tau(\rho) \leq 1, \quad \rho_{cr} \leq \rho \leq 1 \iff 0 \leq \tau(\rho) \leq \rho_{cr}.$$

Given initial data (of Riemann type) $\rho_{1,0}, \dots, \rho_{n+m,0} \in [0, 1]$ we define:

1. for every $i \in \{1, \dots, n\}$

$$\gamma_i^{\max} = \begin{cases} f(\rho_{i,0}), & \text{if } 0 \leq \rho_{i,0} \leq \rho_{cr}, \\ f(\rho_{cr}), & \text{if } \rho_{cr} \leq \rho_{i,0} \leq 1; \end{cases} \qquad (8)$$

2. for every $j \in \{n+1, \dots, n+m\}$

$$\gamma_j^{\max} = \begin{cases} f(\rho_{cr}), & \text{if } 0 \leq \rho_{j,0} \leq \rho_{cr}, \\ f(\rho_{j,0}), & \text{if } \rho_{cr} \leq \rho_{j,0} \leq 1; \end{cases} \qquad (9)$$

3. for every $l \in \{1, \dots, n+m\}$

$$\Omega_l = [0, \gamma_l^{\max}]. \qquad (10)$$
Moreover, we have the following result (see [15]):

**Proposition 2.1** It holds:

1. For every $i \in \{1, \dots, n\}$, an element $\bar{\gamma}$ belongs to $\Omega_i$ if and only if there exists $\bar{\rho}_i \in [0, 1]$ such that $f(\bar{\rho}_i) = \bar{\gamma}$ and point 2 of Definition 2.3 is satisfied.

2. For every $j \in \{n+1, \dots, n+m\}$, an element $\bar{\gamma}$ belongs to $\Omega_j$ if and only if there exists $\bar{\rho}_j \in [0, 1]$ such that $f(\bar{\rho}_j) = \bar{\gamma}$ and point 3 of Definition 2.3 is satisfied.

## 3 Definition of the Priority Riemann Solver

In this section we define a new Riemann solver based on priorities. For this purpose, we first fix a matrix $A$ belonging to the set of matrices:

$$ A := \left\{ A = \{a_{ji}\}_{\substack{i=1,\dots,n \\ j=n+1,\dots,n+m}} : \sum_{j=n+1}^{n+m} a_{ji} = 1 \quad \forall i \right\} \quad (11) $$

and a priority vector $P = (p_1, \dots, p_n) \in \mathbb{R}^n$, with $p_i > 0$, $\sum_i p_i = 1$, indicating priorities among incoming roads.

Consider the closed, convex and non-empty set

$$ \Omega = \left\{ (\gamma_1, \dots, \gamma_n) \in \prod_{i=1}^{n} \Omega_i : A \cdot (\gamma_1, \dots, \gamma_n)^T \in \prod_{j=n+1}^{n+m} \Omega_j \right\}, \quad (12) $$

and define:

$$ \bar{h} = \sup\{h \in \mathbb{R}^+ : hP \in \Omega\}. \quad (13) $$

Given Riemann data $(\rho_{1,0}, \dots, \rho_{m+n,0})$, we define a vector $Q = (\bar{\gamma}_1, \dots, \bar{\gamma}_n)$ of incoming fluxes by a recursive procedure. First we explain the procedure in steps and then provide a pseudo-code in Algorithm 1.

* **STEP 1.** For every $i \in \{1, \dots, n\}$ define

$$ h_i = \max\{h : h p_i \le \gamma_i^{\max}\} = \frac{\gamma_i^{\max}}{p_i}, $$
and for every $j \in \{n+1..., n+m\}$ define

$$h_j = \max \left\{ h : (A \cdot (hP))_j = h \left( \sum_i a_{ji} p_i \right) \le \gamma_j^{\max} \right\} = \frac{\gamma_j^{\max}}{\sum_i a_{ji} p_i}.$$

In other words, $h_i$ is the maximal $t$ so that $hP$ verifies the flux constraint for the $i$-th road, similarly for $h_j$.

$$\text{Set } \bar{h} = \min_{ij} \{h_i, h_j\}.$$

We distinguish two cases:

- CASE 1. If there exists $j$ such that $h_j = \bar{h}$, then we set $Q = \bar{h}P$ and we are done.

- CASE 2. Otherwise, let $I_1 = \{i \in \{1, ..., n\} : h_i = \bar{h}\}$ (by assumption $I_1 \neq \emptyset$). We set $Q_i = \bar{h}p_i$ for $i \in I_1$ and we go to next step.

• **STEP S.** In step $S-1$ we defined a set $I_{S-1}$ and, by induction, all components of $Q$ are fixed for $i \in J_S = I_1 \cup \dots \cup I_{S-1}$. We let $|J_S| < n$ denote the cardinality of $J_S$ and denote by $J_S^c$ the complement of $J_S$ in $\{1, ..., n\}$. We now define $h_i$ for $i \in J_S^c$ by:

$$h_i = \max\{h : h p_i \le \gamma_i^{\max}\} = \frac{\gamma_i^{\max}}{p_i},$$

and for every $j \in \{n+1..., n+m\}$ define

$$h_j = \max \left\{ h : \sum_{i \in J_S} a_{ji} Q_i + h \left( \sum_{i \in J_S^c} a_{ji} p_i \right) \le \gamma_j^{\max} \right\}.$$

We then proceed similarly to STEP 1, setting $\bar{h} = \min_{ij} \{h_i, h_j\}$ and distinguishing two cases:

- CASE 1. If there exists $j$ such that $h_j = \bar{h}$, then we set $Q_i = \bar{h}P_i$ for $i \in J_S^c$ and we are done.

- CASE 2. Otherwise, let $I_S = \{i \in J_S^c : h_i = \bar{h}\}$ (by assumption $I_S \neq \emptyset$). We set $Q_i = \bar{h}p_i$ for $i \in I_S$. If $J_S \cup I_S = \{1, ..., n\}$ then we stop, otherwise we go to next step.
**Algorithm 1 Recursive definition of PRS**

Set $J = \emptyset$ and $J^c = \{1, \dots, n\} \setminus J$.

while $|J| < n$ do

$$ \forall i \in J^c \rightarrow h_i = \max\{h : h p_i \le \gamma_i^{max}\} = \frac{\gamma_i^{max}}{p_i}, $$

$$ \forall j \in \{n+1, \dots, n+m\} \rightarrow h_j = \sup\{h : \sum_{i \in J} a_{ji} Q_i + h(\sum_{i \in J^c} a_{ji} p_i) \le \gamma_j^{max}\}. $$

Set $\hbar = \min_{ij}\{h_i, h_j\}$.

if $\exists j$ s.t. $h_j = \hbar$ then

Set $Q = \hbar P$ and $J = \{1, \dots, n\}$.

else

Set $I = \{i \in J^c : h_i = \hbar\}$ and $Q_i = \hbar p_i$ for $i \in I$.

Set $J = J \cup I$.

end if

end while

We are now ready to define the Priority Riemann Solver (briefly PRS).

**Definition 3.1** Let $Q = (\bar{\gamma}_1, \dots, \bar{\gamma}_n)$ be the vector of incoming fluxes defined by Algorithm 1, then the vector of outgoing fluxes is given by $A \cdot Q^T = (\bar{\gamma}_{n+1}, \dots, \bar{\gamma}_{n+m})^T$.

For every $i \in \{1, \dots, n\}$, set $\bar{\rho}_i$ equal either to $\rho_{i,0}$ if $f(\rho_{i,0}) = \bar{\gamma}_i$, or to the solution to $f(\rho) = \bar{\gamma}_i$ such that $\bar{\rho}_i \ge \rho_{cr}$. For every $j \in \{n+1, \dots, n+m\}$, set $\bar{\rho}_j$ equal either to $\rho_{j,0}$ if $f(\rho_{j,0}) = \bar{\gamma}_j$, or to the solution to $f(\rho) = \bar{\gamma}_j$ such that $\bar{\rho}_j \le \rho_{cr}$. Finally, $PRS: [0, 1]^{n+m} \rightarrow [0, 1]^{n+m}$ is given by

$$ PRS(\rho_{1,0}, \dots, \rho_{n+m,0}) = (\bar{\rho}_1, \dots, \bar{\rho}_n, \bar{\rho}_{n+1}, \dots, \bar{\rho}_{n+m}). \quad (14) $$

# 4 Existence result for Cauchy problems

Given initial data of bounded variation $\rho_{0,1}(\cdot), \dots, \rho_{0,n}(\cdot) : [0, +\infty[ \rightarrow [0, 1]$ and $\rho_{0,n+1}(\cdot), \dots, \rho_{0,n+m}(\cdot) : ]-\infty, 0] \rightarrow [0, 1]$ the corresponding Cauchy problem is defined by:

$$ \left\{ \begin{array}{l} \partial_t \rho_l + \partial_x f(\rho_l) = 0, \\ \rho_l(0, x) = \rho_{0,l}(x), \end{array} \right. \quad l \in \{1, \dots, n+m\}. \quad (15) $$

To solve Cauchy problems one can construct approximate solutions via Wave Front Tracking (WFT). In simple words, one first approximate the initial
data by piecewise constant functions, then solve the corresponding Riemann problems within roads and at junctions approximating rarefaction waves by a fan of rarefaction shocks and solve new Riemann problems when waves interact with each other or with the junction. We refer the reader to [14] for details. Notice that all waves in a WFT approximate solution are shocks, i.e. traveling discontinuities. For every wave we will usually indicate by $\rho_l$, respectively $\rho_r$, the left limit, respectively right limit, of the approximate solution at the discontinuity point. To prove convergence of WFT approximations, one needs to estimate the number of waves, the number of wave interactions and provide estimates on the total variation of approximate solutions. The general theory of [15] is based on three properties which guarantee such estimates. Along the same idea we define three general properties (P1)-(P3) which will ensure existence of solutions.

The first property requires that equilibria are determined only by bad data values (and coincides with (P1) in [15]), more precisely:

**Definition 4.1** We say that a Riemann solver $\mathcal{RS}$ has the property (P1) if the following condition holds. Given $(\rho_{1,0}, \dots, \rho_{n+m,0})$ and $(\rho'_{1,0}, \dots, \rho'_{n+m,0})$ two initial data such that $\rho_{l,0} = \rho'_{l,0}$ whenever either $\rho_{l,0}$ or $\rho'_{l,0}$ is a bad datum, then

$$ \mathcal{RS}(\rho_{1,0}, \dots, \rho_{n+m,0}) = \mathcal{RS}(\rho'_{1,0}, \dots, \rho'_{n+m,0}). \quad (16) $$

The second property requires for bounds in the increase of the flux variation for waves interacting with $J$. More precisely the latter is bounded in terms of the strength of the interacting wave as well as the sum of the changes in the incoming fluxes and in $\bar{h}$ (see (13)). Moreover, the increase in $\bar{h}$ is bounded by the strength of the interacting wave.

**Definition 4.2** We say that a Riemann solver $\mathcal{RS}$ has the property (P2) if there exists a constant $C \ge 1$ such that the following condition holds. For every equilibrium $(\rho_{1,0}, \dots, \rho_{n+m,0})$ of $\mathcal{RS}$ and for every wave $(\rho_i, \rho_{i,0})$ for $i = 1, \dots, n$ (respectively $(\rho_{j,0}, \rho_j)$ for $j = n+1, \dots, n+m$) interacting with $J$ at time $\bar{t} > 0$ and producing waves in the arcs according to $\mathcal{RS}$, we have

$$ \begin{aligned} & \mathrm{TV}_f(\bar{t}+) - \mathrm{TV}_f(\bar{t}-) \\ & \le C \min \left\{ |f(\rho_{l,0}) - f(\rho_l)|, |\Gamma(\bar{t}+) - \Gamma(\bar{t}-)| + |\bar{h}(\bar{t}+) - \bar{h}(\bar{t}-)| \right\} \end{aligned} \quad (17) $$

and

$$ \bar{h}(\bar{t}+) - \bar{h}(\bar{t}-) \le C |f(\rho_{l,0}) - f(\rho_l)|. \quad (18) $$
Finally, we state the third property: if a wave interacts with $J$ and provokes a flux decrease then $\bar{h}$ decreases and the increase of $\Gamma$ is bounded by the change in $\bar{h}$.

**Definition 4.3** We say that a Riemann solver $\mathcal{RS}$ has the property (P3) if the following holds. For every equilibrium $(\rho_{1,0}, \dots, \rho_{n+m,0})$ of $\mathcal{RS}$ and for every wave $(\rho_i, \rho_{i,0})$ with $f(\rho_i) < f(\rho_{i,0})$ for $i = 1, \dots, n$ (respectively $(\rho_{j,0}, \rho_j)$ with $f(\rho_j) < f(\rho_{j,0})$ for $j = n+1, \dots, n+m$) interacting with $J$ at time $\bar{t} > 0$ and producing waves in the arcs according to $\mathcal{RS}$, we have

$$ \Gamma(\bar{t}+) - \Gamma(\bar{t}-) \le C |\bar{h}(\bar{t}+) - \bar{h}(\bar{t}-)|, \quad (19) $$

$$ \bar{h}(\bar{t}+) \le \bar{h}(\bar{t}-). \quad (20) $$

**Theorem 4.1** If a Riemann solver satisfies (P1)-(P3), then every Cauchy problem with initial data of bounded variation admits a weak solution.

In order to prove Theorem 4.1, we need first to provide some definition and results. We start by giving the following:

**Definition 4.4** A wave along a WFT approximate solution generated at time $t = 0$ inside a road is called original, while the ones generated by $J$ are called not original. If two original waves interact, then the resulting wave is still called original, while if an original wave interacts with a not original wave then the resulting wave is not original.

We now have the following result:

**Proposition 4.1** Let $(\rho_l, \rho_r)$ be a wave generated on an incoming road $I$ from the junction at time $\bar{s}$. Assume that there exists a time $\bar{t} > \bar{s}$ at which the wave interacts with $J$ (after interacting with waves inside $I$) and call $\bar{\rho}_l$, respectively $\bar{\rho}_r$, its left, respectively, right limit at $\bar{t}-$. If $I$ is an incoming road then we have $\bar{\rho}_r \ge \rho_{cr}$ and $f(\bar{\rho}_l) < f(\bar{\rho}_r)$. If $I$ is an outgoing road then we have $\bar{\rho}_l \le \rho_{cr}$ and $f(\bar{\rho}_r) < f(\bar{\rho}_l)$.

**Proof.** We prove the result for an incoming road, the other case being similar. First notice that $(\rho_l, \rho_r)$ must have negative speed, thus if $\rho_l < \rho_{cr}$ then $\rho_r > \tau(\rho_r) \ge \rho_{cr}$, while if $\rho_l \ge \rho_{cr}$ then $\rho_r \ge \rho_{cr}$. Therefore in both cases we have $\rho_r \ge \rho_{cr}$. If the wave interacts with waves coming from the left then the value of $\rho_r$ does not change. If the wave interacts with a wave $(\rho_r, \hat{\rho})$
coming from the right, then the wave was generated from the junction J (or obtained by interactions of waves generated from J) and thus must satisfy $\hat{\rho} \ge \rho_{cr}$. Finally, $\bar{\rho}_r \ge \rho_{cr}$. Then, since $(\bar{\rho}_l, \bar{\rho}_r)$ must have positive speed, we deduce that $\bar{\rho}_l < \rho_{cr}$ and $f(\bar{\rho}_l) < f(\bar{\rho}_r)$, thus we conclude. $\square$

From Proposition 4.1 we have the following:

**Corollary 4.1** If a wave $(\rho_l, \rho_r)$ interacts with J from an incoming road and satisfy $f(\rho_l) > f(\rho_r)$ then it is an original wave. If a wave $(\rho_l, \rho_r)$ interacts with J from an outgoing road and satisfy $f(\rho_r) > f(\rho_l)$ then it is an original wave.

**Proof of Theorem 4.1.** WFT approximate solutions can be constructed because (P1) holds true (see [15]).

We now prove that the total variation of the flux $TV(f)$ remains uniformly bounded in time along WFT approximate solutions. The main idea is to first bound the total variation in time of $\bar{h}$ and then of $\Gamma$. This in turn will provide the desired estimate.

Let us indicate with PV the positive variation of a function and with NV the negative one. Then:

$$TV(\bar{h}) = PV(\bar{h}) + NV(\bar{h})$$

and

$$PV(\bar{h}) = PV^O(\bar{h}) + PV^R(\bar{h}),$$

where $PV^O$ is the variation due to interactions of original waves with the junction and $PV^R$ the one due to returning waves (i.e. not original).

From (P2) we get:

$$PV^O(\bar{h}) \le CTV(f(\rho_0)) \le C \max_{\rho \in [0,1]} |f'(\rho)| TV(\rho_0)$$

and from (P3) and Corollary 4.1 it follows $PV^R(\bar{h}) = 0$. Then $PV(\bar{h})$ is bounded and, since $\bar{h} \le f^{\max}/\max_i p_i$, also $TV(\bar{h})$ is bounded. Similarly, for $\Gamma$ we can write

$$TV(\Gamma) = PV(\Gamma) + NV(\Gamma)$$

and

$$PV(\Gamma) = PV^O(\Gamma) + PV^R(\Gamma).$$

Following the proof of [15, Lemma 12], (P2) implies

$$PV^O(\Gamma) \le (C + 2)TV(f(\rho_0)) \le (C + 2) \max_{\rho \in [0,1]} |f'(\rho)| TV(\rho_0).$$
From (P3) we have that

$$PV^R(\Gamma) \leq CTV(\bar{h}),$$

which we just proved to be bounded. Therefore TV(Γ) is also bounded.

Now, define Int the set of times at which a wave interacts with the junction J and, for $s \in Int$ let us indicate by $\Delta TV_f(s)$ the change due to the interaction.

By (P2) we have

$$TV(f(t)) \leq TV(f(\rho_0)) + \sum_{s \in \text{Int}, s \leq t} \Delta TV_f(s) \\ \leq \max_{\rho \in [0,1]} |f'(\rho)| TV(\rho_0) + C(TV(\Gamma) + TV(\bar{h})).$$

Once $TV(f)$ is bounded, one can obtain a bound on $TV(\rho)$ as in [15] and conclude by passing to the limit in WFT approximate solutions. □

**Proposition 4.2** *The Priority Riemann Solver PRS satisfies (P1)-(P3) for junctions with $n \leq 2$, $m \leq 2$ and $0 < a_{ji} < 1$ for all i, j.*

The technical proof is deferred to Section 7.

# 5 Solver with softer priorities

In this section we define a different version of the Riemann solver that uses priorities as softer constraints. In particular, this solver will differ from the solver PRS defined in Section 3 when one of the entries of the matrix A, defined in (11), vanishes, see Figure 1. Notice that the softer priority of the SPRS will allow some flow from road 2 to pass through the junction, when the maximal flow from road 1 is already reached. This reflects the situation where the physical geometry of the junction allows for traffic from road 2 to road 4 (no traffic goes from 2 to 3) even if the traffic from road 1 to road 3 is maximal and has higher priority. For this purpose, we consider a matrix A that may have $a_{ji} = 0$ for some i, j and a priority vector $P = (p_1, \dots, p_n) \in \mathbb{R}^n$, with $p_i > 0$, $\sum_i p_i = 1$.

Then the Riemann solver with softer priorities (briefly SPRS) can be defined by the following recursive algorithm:
<figure><img src="image_2.png" /><figcaption>Figure 1: Different solution for the two solvers: *PRS* and *SPRS*. The distribution matrix *A* satisfies $a_{13} = 0$, i.e. no cars enter road 3 from road 1. Correspondingly the constraint on the flux $\gamma_3$ is represented by a vertical line. The *PRS* selects the point *Q* in (a), while *SPRS* the point *Q* in (b). Since priorities are softer the flux through the junction of the *SPRS* solution is higher than that of the *SPRS*.</figcaption></figure>

**Algorithm 2 Recursive definition of SPRS**

Set $J = \emptyset$ and $J^c = \{1, \dots, n\} \setminus J$.

while $|J| < n$ do

$\forall i \in J^c \rightarrow h_i = \max\{h : hp_i \le \gamma_i^{max}\} = \frac{\gamma_i^{max}}{p_i}$,

$\forall j \in \{n+1, \dots, n+m\} \rightarrow h_j = \sup\{h : \sum_{i \in J} a_{ji}Q_i + h(\sum_{i \in J^c} a_{ji}p_i) \le \gamma_j^{max}\}$.

Set $\hbar = \min_{ij}\{h_i, h_j\}$.

if $\exists j$ s.t. $h_j = \hbar$ then

Set $I = \{i \in J^c : a_{ji} \neq 0\}$ and $Q_i = hp_i$ for $i \in I$.

else

Set $I = \{i \in J^c : h_i = \hbar\}$ and $Q_i = hp_i$ for $i \in I$.

end if

Set $J = J \cup I$.

end while

We are now ready to define the Softer Priority Riemann Solver.
**Definition 5.1** Let $Q = (\bar{\gamma}_1, \dots, \bar{\gamma}_n)$ be the vector of incoming fluxes defined by Algorithm 2, then the vector of outgoing fluxes is given by $A \cdot Q^T = (\bar{\gamma}_{n+1, \dots, n+m})^T$.

For every $i \in \{1, \dots, n\}$, set $\bar{\rho}_i$ equal either to $\rho_{i,0}$ if $f(\rho_{i,0}) = \bar{\gamma}_i$, or to the solution to $f(\rho) = \bar{\gamma}_i$ such that $\bar{\rho}_i \ge \rho_{cr}$. For every $j \in \{n+1, \dots, n+m\}$, set $\bar{\rho}_j$ equal either to $\rho_{j,0}$ if $f(\rho_{j,0}) = \bar{\gamma}_j$, or to the solution to $f(\rho) = \bar{\gamma}_j$ such that $\bar{\rho}_j \le \rho_{cr}$. Finally, $SPRS: [0, 1]^{n+m} \to [0, 1]^{n+m}$ is given by

$$SPRS(\rho_{1,0}, \dots, \rho_{n+m,0}) = (\bar{\rho}_1, \dots, \bar{\rho}_n, \bar{\rho}_{n+1}, \dots, \bar{\rho}_{n+m}). \quad (21)$$

# 6 Numerical scheme and numerical simulations

To illustrate the $PRS$ and $SPRS$ dynamics we provide some simulations based on the well-known Godunov scheme [16] on networks (see [14]), which is based on solutions to Riemann problems.

Define a numerical grid on $[0, T] \times \mathbb{R}$ given by:

* $\Delta x$ is the fixed space grid size;

* $\Delta t^\nu, \nu \in \mathbb{N}$, is the time grid size satisfying the CFL condition [10]:

$$\Delta t^\nu \max_{j \in \mathbb{Z}} |f'(u_j^\nu)| \le \frac{1}{2} \Delta x \quad (22)$$

* $(t^\nu, x_j) = (t^{\nu-1} + \Delta t^\nu, j \Delta x)$ for $\nu \in \mathbb{N}$ and $j \in \mathbb{Z}$ are the grid points.

Consider a scalar conservation laws equipped with initial data:

$$\begin{aligned} \partial_t u + \partial_x f(u) &= 0, & x \in \mathbb{R}, t \in [0, T], \\ u(0, x) &= u_0(x), & x \in \mathbb{R}. \end{aligned} \quad (23)$$

An approximate solution of the problem is constructed first by taking a piecewise constant approximation of the initial data

$$u_j^0 = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} u_0(x) dx, \quad j \in \mathbb{Z}, \quad (24)$$
and then defining $u_j^\nu$ recursively from $u_j^0$ as follows. Under the CFL (22) the waves generated by different Riemann problem at the cell interfaces do not interact and the scheme can be written as follows

$$u_j^{\nu+1} = u_j^{\nu} - \frac{\Delta t^{\nu}}{\Delta x} (G(u_j^{\nu}, u_{j+1}^{\nu}) - G(u_{j-1}^{\nu}, u_j^{\nu})), \quad (25)$$

where the numerical flux $G$ is given by

$$G(u, v) = \begin{cases} \min_{z \in [u, v]} f(z) & \text{if } u \le v \\ \max_{z \in [v, u]} f(z) & \text{if } v \le u. \end{cases} \quad (26)$$

To impose boundary conditions and conditions at junctions we use the classical approach introduced in [7].

### Boundary conditions.

Each road is divided into $M$ cells, numbered from 1 to $M$. Boundary conditions are imposed using ghost cells. For an incoming road $I_i$ we define:

$$u_{i,1}^{\nu+1} = u_{i,1}^{\nu} - \frac{\Delta t^{\nu}}{\Delta x} (G(u_{i,1}^{\nu}, u_{i,2}^{\nu}) - G(u_{i,0}^{\nu}, u_{i,1}^{\nu})) \quad (27)$$

where $u_{i,0}^\nu$ is the value of the density at the boundary.

The outgoing boundary for $I_j$ is treated in the same way by defining:

$$u_{j,M}^{\nu+1} = u_{j,M}^{\nu} - \frac{\Delta t^{\nu}}{\Delta x} (G(u_{j,M}^{\nu}, u_{j,M+1}^{\nu}) - G(u_{j,M-1}^{\nu}, u_{j,M}^{\nu})) \quad (28)$$

with $u_{i,M+1}^\nu$ the value of the density at the outgoing boundary.

### Conditions at the junction.

For $I_i$ with $i \in \{1, \dots, n\}$ that is connected at the junction at the right endpoint we set:

$$u_{i,M}^{\nu+1} = u_{i,M}^{\nu} - \frac{\Delta t^{\nu}}{\Delta x} (Q_i^{\nu} - G(u_{i,M-1}^{\nu}, u_{i,M}^{\nu})) \quad \text{for } i \in \{1, \dots, n\}, \quad (29)$$

while for the outgoing roads, connected at the junction with the left endpoint we have:

$$u_{j,1}^{\nu+1} = u_{j,1}^{\nu} - \frac{\Delta t^{\nu}}{\Delta x} (G(u_{j,1}^{\nu}, u_{j,2}^{\nu}) - Q_j^{\nu}) \quad \text{for } j \in \{n+1, \dots, n+m\}, \quad (30)$$

where $Q_i^\nu, Q_j^\nu$ are the incoming and outgoing fluxes given by the Riemann solvers at junction corresponding to the initial data

$$(u_{1,M}^{\nu}, \dots, u_{n,M}^{\nu}, u_{n+1,1}^{\nu}, \dots, u_{n+m,1}^{\nu})$$

(see Algorithms 1, 2).
## 6.1 Numerical results

For the simulations, we set the length of each road equal to 1 and incoming roads are parametrized by the interval $I_i = [-1, 0]$ while outgoing roads are given by $I_j = [0, 1]$, with the junction placed at $x = 0$. Moreover, we fix $f(\rho) = \rho(1 - \rho)$, thus $\rho_{cr} = 0.5$.

1. **Case I:** Comparison *PRS* vs. *SPRS*.

This case illustrates the different dynamics given by the two Riemann solvers proposed in this article.

We consider a junction with 2 incoming roads ($I_1, I_2$) and 2 outgoing roads ($I_3, I_4$). We fix the matrix $A$ and the priority vector $P$ as follows:

$$A = \begin{bmatrix} 0.6 & 0 \\ 0.4 & 1 \end{bmatrix} \quad P = \begin{bmatrix} 0.7 & 0.3 \end{bmatrix}. \qquad (31)$$

We consider the following initial data:

$$\rho_{1,0} = 0.6, \quad \rho_{2,0} = 0.2, \quad \rho_{3,0} = 0.85, \quad \rho_{4,0} = 0.2. \qquad (32)$$

The different results of the simulations (see Figures 2) can be seen in particular in road 2 ad 4. We observe that *SPRS* allows more flux through the junction than *PRS*, for which we observe the formation of a big shock moving backwards on road 2.

2. **Case II:** Comparison *PRS* vs. *RS*_{CGP}.

We propose here a comparison between the *PRS* with the Riemann solver proposed by Coclite, Garavello and Piccoli in [9] and briefly referred to as *RS*_{CGP}.

We consider a 2 × 2 junction and we fix the matrix $A$ and the priority vector $P$ as follows:

$$A = \begin{bmatrix} 0.5 & 0.6 \\ 0.5 & 0.4 \end{bmatrix} \quad P = \begin{bmatrix} 0.7 & 0.3 \end{bmatrix}. \qquad (33)$$

We consider the following initial data:

$$\rho_{1,0} = 0.2, \quad \rho_{2,0} = 0.6, \quad \rho_{3,0} = 0.3, \quad \rho_{4,0} = 0.8. \qquad (34)$$

The simulations (see Figures 3) show clearly the different solutions of the Riemann solvers. In particular, *RS*_{CGP} creates a big shock in the incoming road 1 decreasing its flux. This wave does not appear in our Riemann solver *PRS*.
<figure><img src="image_3.png" /><figcaption>Figure 2: Case I: Solution of the problem using *PRS* on the left and *SPRS* on the right.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3: Case II: Solution of the problem using *PRS* on the left and *RS<sub>CGP</sub>* on the right.</figcaption></figure>
### 3. Case III: 3 × 2 junction.

We fix the matrix A and the priority vector P as follows:

$$A = \begin{bmatrix} 0.5 & 0.6 & 0.2 \\ 0.5 & 0.4 & 0.8 \end{bmatrix} \quad P = \begin{bmatrix} 0.5 & 0.3 & 0.2 \end{bmatrix}. \qquad (35)$$

We consider the following initial data:

$$\rho_{1,0} = 0.2, \quad \rho_{2,0} = 0.6, \quad \rho_{3,0} = 0.3, \quad \rho_{4,0} = 0.8, \quad \rho_{5,0} = 0.2. \qquad (36)$$

<figure><img src="image_5.png" /><figcaption>Figure 4: Case III: Solution of the problem using PRS</figcaption></figure>
Due to the lower priorities given to roads 2 and 3 we can see that queues are created in the two incoming roads, see Figure 4. Note also that this case cannot be handled by $RS_{CGP}$ since $n > m$.

# 7 Appendix: Proof of Proposition 4.2

The construction of $PRS$ depends only on the matrix $A$, the priority vector $P$ and the sets $\Omega_l$. The latter, in turn, depends only on bad data, thus property (P1) holds true.

We prove (P2) and (P3) for the case $n = m = 2$ and distinguish between three different generic situations for the initial equilibrium: demand constrained, demand/supply constrained and supply constrained (where demand indicates flow from incoming road and supply flow to outgoing ones). In the first situation the incoming roads act as constraint in the definition of the set $\Omega$ (see (12)) and the equilibrium corresponds to the point $Q_0$ as in Figure 5a. The second case corresponds to one incoming and one outgoing road acting as constraint and to the point $Q_0$ as in Figure 8a. Finally, the third case corresponds to outgoing roads acting as constraint and to the point $Q_0$ as in Figure 12a.

* **Case A: Demand constrained.** By symmetries, it is not restrictive to assume that the priority line $hP$, $h > 0$, intersects the constraint $\gamma_2^{\text{inc}} = \gamma_{2,0}$. We have to distinguish several subcases:

**Case A1:** The incoming wave is $(\rho_1, \rho_{1,0})$ (on road 1). Since $\gamma_{1,0}$ is an active constraint, $\rho_{1,0} \le \rho_{cr}$ and $\rho_1 \le \rho_{cr}$. We distinguish the two situations:

- If $f(\rho_1) > f(\rho_{1,0})$ we define $\gamma_{1,1} = \frac{\gamma_3 - a_{32}\gamma_2}{a_{31}}$ and $\gamma_{1,0} \le \gamma_{1,1} \le \gamma_1$ (see Figure 5a). We get:

$$ 
\begin{aligned}
\mathrm{TV}(f)^+ &= |\gamma_1 - \gamma_{1,1}| + a_{31} |\gamma_{1,1} - \gamma_{1,0}| + a_{41} |\gamma_{1,1} - \gamma_{1,0}|, \\
\mathrm{\Delta TV}(f) &= (a_{31} + a_{41} - 1) |\gamma_{1,1} - \gamma_{1,0}|, \\
\mathrm{\Delta \Gamma} &= \gamma_{1,1} - \gamma_{1,0} > 0, \\
\mathrm{\Delta \bar{h}} &= 0.
\end{aligned}
$$

Hence, (P2) holds and (P3) doesn't need to be verified.
If $f(\rho_1) > f(\rho_{1,0})$, we define $\gamma_{1,1} = \frac{p_1}{p_2}\gamma_2$, hence $\gamma_{1,0} \ge \gamma_{1,1} \ge \gamma_1$ (see Figure 5b). We have:

$$
\begin{align*}
\mathrm{TV}(f)^- &= |\gamma_1 - \gamma_{1,0}|, \\
\mathrm{TV}(f)^+ &= a_{31} |\gamma_1 - \gamma_{1,0}| + a_{41} |\gamma_1 - \gamma_{1,0}|, \\
\Delta\mathrm{TV}(f) &= (a_{31} + a_{41} - 1) |\gamma_1 - \gamma_{1,0}|, \\
\Delta\Gamma &= (\gamma_1 - \gamma_{1,0}) < 0, \\
\Delta\bar{h} &= \frac{1}{p_1}(\gamma_1 - \gamma_{1,1}) < 0.
\end{align*}
$$

Hence, (P2) and (P3) hold.

<figure><img src="image_6.png" /><figcaption>Figure 5: Case A1</figcaption></figure>

**Case A2:** The incoming wave is $(\rho_2, \rho_{2,0})$ (on road 2). Since $\gamma_{2,0}$ is an active constraint, $\rho_{2,0} \le \rho_{cr}$ and $\rho_2 \le \rho_{cr}$.

If $f(\rho_2) > f(\rho_{2,0})$, we define

$$
\gamma_{2,1} = \frac{\gamma_3 - a_{31}\gamma_1}{a_{32}} \quad \text{and} \quad \gamma_{2,2} = \frac{p_2}{a_{31}p_1 + a_{32}p_2} \gamma_3^{\max}
$$

so that $\gamma_{2,0} \le \gamma_{2,1} \le \gamma_{2,2} \le \gamma_2$, see Figure 6a. Note that this case is the same as in A1 except for the case in the drawing. In this case we have:

$$
\mathrm{TV}(f)^- = |\gamma_2 - \gamma_{2,0}| = |\gamma_2 - \gamma_{2,2}| + |\gamma_{2,2} - \gamma_{2,1}| + |\gamma_{2,1} - \gamma_{2,0}|,
$$
$$ \begin{aligned} \mathrm{TV}(f)^+ = & \frac{a_{32}}{a_{31}} |\gamma_{2,2} - \gamma_{2,1}| + |\gamma_2 - \gamma_{2,2}| + a_{32} |\gamma_{2,1} - \gamma_{2,0}| \\ & + \left| a_{42} |\gamma_{2,1} - \gamma_{2,0}| + \left(a_{42} - a_{41} \frac{a_{32}}{a_{31}}\right) |\gamma_{2,2} - \gamma_{2,1}| \right|, \end{aligned} $$

$$ \begin{aligned} \Delta \mathrm{TV}(f) = & \left( \frac{a_{32}}{a_{31}} - 1 \right) |\gamma_{2,2} - \gamma_{2,1}| + (a_{32} - 1) |\gamma_{2,1} - \gamma_{2,0}| \\ & + \left| a_{42} |\gamma_{2,1} - \gamma_{2,0}| + \left(a_{42} - a_{41} \frac{a_{32}}{a_{31}}\right) |\gamma_{2,2} - \gamma_{2,1}| \right|, \end{aligned} $$

$$ \Delta\Gamma = (\gamma_{2,1} - \gamma_{2,0}) + \left(1 - \frac{a_{32}}{a_{31}}\right) (\gamma_{2,2} - \gamma_{2,1}), $$

$$ \Delta \bar{h} = \frac{1}{p_2}(\gamma_{2,2} - \gamma_{2,0}) $$

Hence, (P2) holds while (P3) doesn't need to be checked.

— If $f(\rho_2) < f(\rho_{2,0})$ one has $\gamma_2 < \gamma_{2,0}$, see Figure 6b. Therefore:

$$ \mathrm{TV}(f)^- = |\gamma_2 - \gamma_{2,0}|, $$

$$ \mathrm{TV}(f)^+ = (a_{32} + a_{42}) |\gamma_2 - \gamma_{2,0}|, $$

$$ \Delta \mathrm{TV}(f) = 0 \quad (\text{recall that } a_{32} + a_{42} = 1), $$

$$ \Delta\Gamma = (\gamma_2 - \gamma_{2,0}) < 0, $$

$$ \Delta \bar{h} = \frac{1}{p_2}(\gamma_2 - \gamma_{2,0}) < 0. $$

Hence, (P2) and (P3) hold.

**Case A3:** The incoming wave is $(\rho_{3,0}, \rho_3)$ (on road 3, the case of road 4 being similar).

— If $f(\rho_3) < f(\rho_{3,0})$ We define $\gamma_{3,1} = \left(a_{31}\frac{p_1}{p_2} + a_{32}\right)\gamma_2$ so that $\gamma_3 \le \gamma_{3,1} \le \gamma_{3,0}$ (see Figure 7):

$$ \mathrm{TV}(f)^- = |\gamma_3 - \gamma_{3,0}|, $$

$$ \mathrm{TV}(f)^+ = \frac{a_{41} + 1}{a_{31}} |\gamma_{3,1} - \gamma_{3,0}| + \frac{(a_{41} + 1)p_1 + (a_{42} + 1)p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_3 - \gamma_{3,1}|, $$

$$ \Delta \mathrm{TV}(f) = \frac{a_{41}}{a_{31}} |\gamma_{3,1} - \gamma_{3,0}| + \left( \frac{(a_{41} + 1)p_1 + (a_{42} + 1)p_2}{a_{31}p_1 + a_{32}p_2} - 1 \right) |\gamma_3 - \gamma_{3,1}|, $$

$$ \Delta\Gamma = \frac{1}{a_{31}}(\gamma_{3,1} - \gamma_{3,0}) + \frac{a_{41}p_1 + a_{42}p_2}{a_{31}p_1 + a_{32}p_2}(\gamma_3 - \gamma_{3,1}) < 0, $$
<figure><img src="image_7.png" /><figcaption>Figure 6: Case A2</figcaption></figure>

$$ \Delta \bar{h} = \frac{1}{a_{31}p_1 + a_{32}p_2} (\gamma_3 - \gamma_{3,1}) < 0. $$

Hence, (P2) and (P3) hold.

— If $f(\rho_3) > f(\rho_{3,0})$ and $\gamma_3 > \gamma_{3,0}$ then we stay demand constrained and nothing happens.

* **Case B: Supply constrained, priority line intersects a demand constraint.** Even in this case, it is not restrictive to assume that the priority line $hP$, $h > 0$, intersects the constraint $\gamma_1^{\text{inc}} = \gamma_{1,0}$. We spit the proof in several subcases depending on the origin of the incoming wave:

**Case B1:** The incoming wave is $(\rho_1, \rho_{1,0})$ (on road 1). We define $\gamma_{1,1} = \frac{a_{42} - a_{32}}{a_{31}a_{42} - a_{32}a_{41}} \gamma_1$ (see Figure 8a and 8b). We distinguish the two situations:

— If $f(\rho_1) > f(\rho_{1,0})$ we get $\gamma_{1,0} \le \gamma_{1,1} \le \gamma_{1,2} \le \gamma_1$, where $\gamma_{1,2} = \frac{p_1}{a_{41}p_1 + a_{42}p_2} \gamma_4$ (see Figure 8a). Then we have:

$$ \text{TV}(f)^- = |\gamma_1 - \gamma_{1,0}|, $$
<figure><img src="image_8.png" /><figcaption>Figure 7: Case a3 - Decreasing $\gamma_3$</figcaption></figure>

$$TV(f)^+ = |\gamma_1 - \gamma_{1,2}| + \left( \frac{a_{41}}{a_{42}} + \left| a_{31} - a_{32} \frac{a_{41}}{a_{42}} \right| \right) |\gamma_{1,2} - \gamma_{1,1}| \\ + \left( \frac{a_{31}}{a_{32}} + \left| a_{41} - a_{42} \frac{a_{31}}{a_{32}} \right| \right) |\gamma_{1,1} - \gamma_{1,0}|,$$

$$\Delta TV(f) = \left( \frac{a_{41}}{a_{42}} + \left| a_{31} - a_{32} \frac{a_{41}}{a_{42}} \right| - 1 \right) |\gamma_{1,2} - \gamma_{1,1}| \\ + \left( \frac{a_{31}}{a_{32}} + \left| a_{41} - a_{42} \frac{a_{31}}{a_{32}} \right| - 1 \right) |\gamma_{1,1} - \gamma_{1,0}|,$$

$$\Delta \Gamma = \left(1 - \frac{a_{31}}{a_{32}}\right) (\gamma_{1,1} - \gamma_{1,0}) + \left(1 - \frac{a_{41}}{a_{42}}\right) (\gamma_{1,2} - \gamma_{1,1}),$$

$$\Delta \bar{h} = \frac{1}{p_1}(\gamma_{1,2} - \gamma_{1,0}).$$

Hence, (P2) holds and (P3) doesn't need to be checked.

If $f(\rho_1) < f(\rho_{1,0})$ and $\gamma_{1,0} \ge \gamma_{1,1} \ge \gamma_{1,2} \ge \gamma_1$ (see Figure 8b), we define $\gamma_{1,2} = \frac{\gamma_3 - a_{32}\gamma_2}{a_{31}}$ and we get:

$$TV(f)^- = |\gamma_1 - \gamma_{1,0}|,$$
$$TV(f)^{+} = (a_{31} + a_{41}) |\gamma_{1} - \gamma_{1,2}| + \left( \frac{a_{31}}{a_{32}} + \left| a_{41} - a_{42} \frac{a_{31}}{a_{32}} \right| \right) |\gamma_{1,2} - \gamma_{1,1}| \\ + \left( \frac{a_{41}}{a_{42}} + \left| a_{31} - a_{32} \frac{a_{41}}{a_{42}} \right| \right) |\gamma_{1,1} - \gamma_{1,0}|,$$

$$\Delta TV(f) = \left( \frac{a_{31}}{a_{32}} + \left| a_{41} - a_{42} \frac{a_{31}}{a_{32}} \right| - 1 \right) |\gamma_{1,2} - \gamma_{1,1}| \\ + \left( \frac{a_{41}}{a_{42}} + \left| a_{31} - a_{32} \frac{a_{41}}{a_{42}} \right| - 1 \right) |\gamma_{1,1} - \gamma_{1,0}|,$$

$$\Delta \Gamma = (\gamma_1 - \gamma_{1,2}) + \left(1 - \frac{a_{31}}{a_{32}}\right) (\gamma_{1,2} - \gamma_{1,1}) + \left(1 - \frac{a_{41}}{a_{42}}\right) (\gamma_{1,1} - \gamma_{1,0}),$$

$$\Delta \bar{h} = \frac{1}{p_1}(\gamma_1 - \gamma_{1,0}) < 0.$$

Hence, (P2) and (P3) hold.

**Case B2:** The incoming wave is $(\rho_2, \rho_{2,0})$ (on road 2).

<figure><img src="image_9.png" /><figcaption>Figure 8: Case B1</figcaption></figure>

– If $f(\rho_2) > f(\rho_{2,0})$: nothing happens (see Figure 9a).

– If $f(\rho_2) < f(\rho_{2,0})$ (see Figure 9b) we define $\gamma_{2,1} = \frac{p_1}{p_2}\gamma_1$ and we
compute:

$$TV(f)^- = |\gamma_2 - \gamma_{2,0}|,$$

$$TV(f)^+ = (a_{32} + a_{42}) |\gamma_2 - \gamma_{2,0}|,$$

$$\Delta TV(f) = (a_{32} + a_{42} - 1) |\gamma_2 - \gamma_{2,0}| = 0,$$

$$\Delta \Gamma = (\gamma_2 - \gamma_{2,0}) < 0,$$

$$\Delta \bar{h} = \frac{1}{p_2}(\gamma_2 - \gamma_{2,1}) < 0.$$

Hence, (P2) and (P3) hold.

**Case B3:** The incoming wave is $(\rho_{3,0}, \rho_3)$ (on road 3).

<figure><img src="image_10.png" /><figcaption>Figure 9: Case B2</figcaption></figure>

– If $f(\rho_3) > f(\rho_{3,0})$ (see Figure 10a), we define

$$\gamma_{3,1} = \min \begin{cases} a_{31}\gamma_1 + a_{32}\gamma_2^{\max}, \\ a_{31}\gamma_1 + \frac{a_{32}}{a_{42}}(\gamma_4^{\max} - a_{41}\gamma_1). \end{cases}$$

We get:

$$TV(f)^- = |\gamma_3 - \gamma_{3,0}|,$$
$$TV(f)^+ = \left( \frac{1}{a_{32}} + \frac{a_{42}}{a_{32}} \right) |\gamma_{3,1} - \gamma_{3,0}| + |\gamma_3 - \gamma_{3,1}|,$$

$$\Delta TV(f) = \left( \frac{1}{a_{32}} + \frac{a_{42}}{a_{32}} - 1 \right) |\gamma_{3,1} - \gamma_{3,0}|,$$

$$\Delta \Gamma = \frac{1}{a_{32}}(\gamma_{3,1} - \gamma_{3,0}) > 0,$$

$$\Delta \bar{h} = 0.$$

Hence, (P2) holds and (P3) does not need to be checked.

— If $f(\rho_3) < f(\rho_{3,0})$ (see Figure 10b), we define $\gamma_{3,1} = \left(a_{31} + a_{32}\frac{p_2}{p_1}\right)\gamma_1$ and we compute:

$$TV(f)^- = |\gamma_3 - \gamma_{3,0}|,$$

$$TV(f)^+ = \frac{p_1}{a_{31}p_1 + a_{32}p_2} |\gamma_3 - \gamma_{3,1}| + \frac{1}{a_{32}} |\gamma_{3,1} - \gamma_{3,0}| + \frac{p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_3 - \gamma_{3,1}| \\ + \frac{a_{42}}{a_{32}} |\gamma_{3,1} - \gamma_{3,0}| + \frac{a_{41}p_1 + a_{42}p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_3 - \gamma_{3,1}|,$$

$$\Delta TV(f) = \left( \frac{(1 + a_{41})p_1 + (1 + a_{42})p_2}{a_{31}p_1 + a_{32}p_2} - 1 \right) |\gamma_3 - \gamma_{3,1}| + \left( \frac{1 + a_{42}}{a_{32}} - 1 \right) |\gamma_{3,1} - \gamma_{3,0}|,$$

$$\Delta \Gamma = \frac{p_1 + p_2}{a_{31}p_1 + a_{32}p_2}(\gamma_3 - \gamma_{3,1}) + \frac{1}{a_{32}}(\gamma_{3,1} - \gamma_{3,0}) < 0,$$

$$\Delta \bar{h} = \frac{1}{a_{31}p_1 + a_{32}p_2}(\gamma_3 - \gamma_{3,1}) < 0.$$

Hence, (P2) and (P3) hold.

**Case B4:** The incoming wave is $(\rho_{4,0}, \rho_4)$ (on road 4).

— If $f(\rho_4) > f(\rho_{4,0})$ nothing changes.

— If $f(\rho_4) < f(\rho_{4,0})$: the same as case **B3** decreasing exchanging the roles of $\gamma_3$ and $\gamma_4$.

* **Case C: Supply constrained, priority line intersects a supply constraint.** It is not restrictive to assume that the priority line $hP$, $h > 0$, intersects the constraint $a_{31}\gamma_1^{\text{inc}} + a_{32}\gamma_2^{\text{inc}} = \gamma_{3,0}$. We distinguish the following subcases, depending on the origin of the incoming wave:
<figure><img src="image_11.png" /><figcaption>Figure 10: Case B3</figcaption></figure>

**Case C1:** The incoming wave is $(\rho_1, \rho_{1,0})$ (on road 1).

– If $f(\rho_1) > f(\rho_{1,0})$, nothing changes.

– If $f(\rho_1) < f(\rho_{1,0})$, the analysis is similar to case **B1** decreasing (see Figure 11).

**Case C2:** The incoming wave is $(\rho_2, \rho_{2,0})$ (on road 2). This case is symmetric to **C1**.

**Case C3:** The incoming wave is $(\rho_{3,0}, \rho_3)$ (on road 3).

– If $f(\rho_3) > f(\rho_{3,0})$ (see Figure 12a), we define

$$ \gamma_{3,1} = \min \left\{ \begin{array}{l} \left( a_{31} + a_{32} \frac{p_2}{p_1} \right) \gamma_1^{\max}, \\ \left( a_{31} \frac{p_1}{p_2} + a_{32} \right) \gamma_2^{\max}, \\ \frac{a_{31} p_1 + a_{32} p_2}{a_{41} p_1 + a_{42} p_2} \gamma_4^{\max}, \end{array} \right. $$
<figure><img src="image_12.png" /><figcaption>Figure 11: Case C1 - decreasing $\gamma_1$</figcaption></figure>

and we compute:

$$TV(f)^- = |\gamma_3 - \gamma_{3,0}|,$$

$$TV(f)^+ = \frac{p_1 + p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_{3,1} - \gamma_{3,0}| + |\gamma_3 - \gamma_{3,1}| + \frac{a_{41}p_1 + a_{42}p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_{3,1} - \gamma_{3,0}|,$$

$$\Delta TV(f) = \left( \frac{(a_{41} + 1)p_1 + (a_{42} + 1)p_2}{a_{31}p_1 + a_{32}p_2} - 1 \right) |\gamma_{3,1} - \gamma_{3,0}|,$$

$$\Delta \Gamma = \frac{p_1 + p_2}{a_{31}p_1 + a_{32}p_2} (\gamma_{3,1} - \gamma_{3,0}) > 0,$$

$$\Delta \bar{h} = \frac{1}{a_{31}p_1 + a_{32}p_2} (\gamma_{3,1} - \gamma_{3,0}) > 0.$$

Hence, (P2) holds and (P3) does not count since we are increasing fluxes.

– If $f(\rho_3) < f(\rho_{3,0})$ (see Figure 12b) we get:

$$TV(f)^- = |\gamma_3 - \gamma_{3,0}|,$$

$$TV(f)^+ = \frac{p_1 + p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_3 - \gamma_{3,0}| + \frac{a_{41}p_1 + a_{42}p_2}{a_{31}p_1 + a_{32}p_2} |\gamma_3 - \gamma_{3,0}|,$$

$$\Delta TV(f) = \left( \frac{(a_{41} + 1)p_1 + (a_{42} + 1)p_2}{a_{31}p_1 + a_{32}p_2} - 1 \right) |\gamma_3 - \gamma_{3,0}|,$$

$$\Delta \Gamma = \frac{p_1 + p_2}{a_{31}p_1 + a_{32}p_2} (\gamma_3 - \gamma_{3,0}) < 0,$$
$$ \Delta \bar{h} = \frac{1}{a_{31}p_1 + a_{32}p_2} (\gamma_3 - \gamma_{3,0}) < 0. $$

Hence, (P2) and (P3) hold.

**Case C4:** The incoming wave is $(\rho_{4,0}, \rho_4)$ (on road 4).

<figure><img src="image_13.png" /><figcaption>Figure 12: Case C3</figcaption></figure>

- If $f(\rho_4) > f(\rho_{4,0})$, nothing happens.

- If $f(\rho_4) < f(\rho_{4,0})$, the situation is similar to case C3 with the roles of $\gamma_3$ and $\gamma_4$ reversed.

## References

[1] D. Armbruster, P. Degond, and C. Ringhofer. A model for the dynamics of large queuing networks and supply chains. *SIAM J. Appl. Math.*, 66(3):896–920, 2006.

[2] C. Bardos, A. Y. le Roux, and J.-C. Nédélec. First order quasilinear equations with boundary conditions. *Comm. Partial Differential Equations*, 4(9):1017–1034, 1979.
[3] Georges Bastin, Alexandre M. Bayen, Ciro D'Apice, Xavier Litrico, and Benedetto Piccoli. Open problems and research perspectives for irrigation channels. *Netw. Heterog. Media*, 4(2):i-v, 2009.

[4] A. Bressan. *Hyperbolic systems of conservation laws: the one-dimensional Cauchy problem*. Oxford university press, 2000.

[5] A. Bressan and A. Nordli. The Riemann solver for traffic flow at an intersection with buffer of vanishing size. to appear., 2016.

[6] Alberto Bressan, Sunčica Čanić, Mauro Garavello, Michael Herty, and Benedetto Piccoli. Flows on networks: recent results and perspectives. *EMS Surv. Math. Sci.*, 1(1):47-111, 2014.

[7] G. Bretti, R. Natalini, and B. Piccoli. Fast algorithms for a traffic flow model on networks. *Discrete and Continuous Dynamical Systems - Series B*, 6(3):427-448, 2006.

[8] Y. Chitour and B. Piccoli. Traffic circles and timing of traffic lights for cars flow. *Discrete Contin. Dyn. Syst. Ser. B*, 5(3):599-630, August 2005.

[9] G.M. Coclite, M. Garavello, and B. Piccoli. Traffic flow on a road network. *SIAM J. Math. Anal.*, 36(6):1862-1886, 2005.

[10] R. Courant, K. Friedrichs, and H. Lewy. On the partial differential equations of mathematical physics. *IBM journal*, 11(2):215-234, 1967.

[11] C.F. Daganzo. The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory. *Transportation Research Part B*, 28:269-287, 1994.

[12] Ciro D'Apice and Benedetto Piccoli. Vertex flow models for vehicular traffic on networks. *Math. Models Methods Appl. Sci.*, 18(suppl.):1299-1315, 2008.

[13] M. L. Delle Monache, J. Reilly, S. Samaranayake, W. Krichene, P. Goatin, and A. M. Bayen. A PDE-ODE model for a junction with ramp buffer. *SIAM Journal on Applied Mathematics*, 74(1):22-39, 2014.
[14] M. Garavello and B. Piccoli. *Traffic flow on networks*, volume 1 of *AIMS Series on Applied Mathematics*. American Institute of Mathematical Sciences (AIMS), Springfield, MO, 2006. Conservation laws models.

[15] M. Garavello and B. Piccoli. Conservation laws on complex networks. *Ann. I. H. Poincaré*, 26:1925–1951, 2009.

[16] S. K. Godunov. A finite difference method for the numerical computation of discontinuous solutions of the equations of fluid dynamics. *Mathematiccheckii Sbornik*, 47:271–290, 1959.

[17] M. Herty, J. Lebacque, and S. Moutari. A novel model for intersections of vehicular traffic flow. *Netw. Heterog. Media*, 4:813–826, 2009.

[18] Helge Holden and Nils Henrik Risebro. A mathematical model of traffic flow on a network of unidirectional roads. *SIAM J. Math. Anal.*, 26(4):999–1017, 1995.

[19] S. N. Kružhkov. First order quasilinear equations with several independent variables. *Mathematiccheckii Sbornik*, 81(123):228–255, 1970.

[20] R.J. LeVeque. *Finite Volume Methods for Hyperbolic Problems*. Cambridge Texts in Applied Mathematics. Cambridge University Press, 2002.

[21] M. J. Lighthill and G. B. Whitham. On kinematic waves. II. A theory of traffic flow on long crowded roads. *Proc. Roy. Soc. London Ser. A*, 229:317–346, 1955.

[22] Alessia Marigo and Benedetto Piccoli. A fluid dynamic model for T-junctions. *SIAM J. Math. Anal.*, 39(6):2016–2032, 2008.

[23] P. I. Richards. Shock waves on the highway. *Operations Research*, 4:42–51, 1956.

[24] S. Samaranayake, J. Reilly, W. Krichene, M. L. Delle Monache, P. Goatin, and A. Bayen. Discrete-time system optimal dynamic traffic assignment (SO-DTA) with partial control for horizontal queuing networks. preprint, https://hal.inria.fr/hal-01095707, 2014.