Emissions minimization on road networks via Generic Second Order Models

1.   Introduction

The aim of this paper is to build a model to estimate and minimize traffic emissions by regulating traffic dynamics. Such regulation corresponds to the choice of suitable model parameters, which in turn represent traffic signals and traffic light timing. Specifically, we extend the Generic Second Order Model (GSOM), introduced in ^[3,29], to road networks, pair it to an emission model and then minimize a functional comprising ${\mathrm{NO_{x}}}$ emissions and travel time.

Estimating traffic emissions is an important and challenging problem. First, most emission models are based on the knowledge of vehicle speed and acceleration. Thus, at macroscopic level, a first-order system based only on conservation of cars, such as the Lighthill-Whitham-Richards (LWR) model ^[31,35], is not sufficient to feed an emission model. In particular, the LWR model admits solutions with unbounded acceleration ^[30] which would lead to unrealistic emissions estimates. It is then necessary to consider a so-called second-order model, i.e., a model with two equations: a first equation for the conservation of mass and a second for the conservation or balance of a modified momentum, which may model drivers' property. The first second-order model goes back to Payne and Whitham ^[33,38]. After criticisms to the model, see ^[12], a new line of research originated starting with the Aw-Rascle-Zhang (ARZ) model ^[5,39], which successfully addressed criticisms to the Payne-Whitham approach. More recently, various second-order models were proposed ranging from generalizations of the ARZ, such as in ^[14,17], to phase transition models as in ^[8,10] and GSOM in ^[3,27,29]. Such models are characterized by a family of fundamental diagrams (density-flow functions) and, due to their multi-faceted nature, are particularly appropriate to fit real traffic data. We refer to ^[15,34] for more details on data-fitted second order models.

Traffic models on networks have been widely studied in last two decades and authors have considered many different traffic scenarios proposing a rich amount of alternative models at junctions. The LWR model has been extended to road networks in several papers, see for example ^{[13,18,20,25]}. The ARZ model on networks was considered in ^[19,23,24] and phase-transition models in ^[11,21]. In this paper we consider a road network with merge (two incoming and one outgoing roads) and diverge (one incoming and two outgoing roads) junctions. On each road, we assume that the traffic flow evolution is described by the GSOM

where $\rho$ is the density of vehicles, $v = V(\rho, w)$ is the velocity function, and $w$ is a property of drivers. Notice that the first equation in (1.1) models the conservation of cars, while the second is the passive advection of the variable $w$, which gives rise to different fundamental diagrams. To define the solution on the whole network we follow the approach proposed in ^[18] based on the concept of Riemann Problem at a junction, which is a Cauchy problem with constant initial data on each road. Solutions to Riemann Problems are required to maximize the flux while conserving the density $\rho$ and total property $y = \rho w$ through the junction. To determine a unique solution to Riemann Problems, we need to introduce additional criteria, which depend on the type of junction. For diverge junctions, a traffic distribution parameter is assigned to outgoing roads as done in ^[24] for the ARZ model. For merge junctions, a priority rule between incoming roads is considered, as it was done for the LWR model in ^[9]. More precisely, for a fixed priority parameter $\beta\in[0, 1]$, given the two incoming fluxes ${\hat q}_{1}$, ${\hat q}_{2}$, we require:

Equation (1.2) establishes a proportional relationship between the two incoming fluxes. For instance, if $\beta = 0$ only traffic from the first road is allowed and vice versa for $\beta = 1$. Therefore, traffic lights can be easily represented by time-varying priority parameters. This rule, together with the maximization of flux and conservation of $\rho$ and $y$, determines unique values of the variable $w$ on each road. In fact, the value ${\hat w}_3$ on the outgoing road is given by a convex combination of the values ${\hat w}_1$ and ${\hat w}_2$ of the two incoming roads, i.e.,

As a result, the maximal flux that can be received by the outgoing road, i.e. the supply, depends on the priority rule. The final solution is determined by maximizing the flow through the junction respecting the priority rule, but relaxing the latter in case the supply exceeds the demand from the road with higher priority. In rough words, the supply is given to incoming roads according to the priority rule and redistributed in case of surplus. The complete procedure to build the solution for a merge junction is explained in details in Definition 3.4. The extension of the proposed algorithms to a generic junction with $n$ incoming and $m$ outgoing roads can be found in ^[6] together with some preliminary results on the bounds on the total variation of junction waves.

The solution on networks to GSOM is then used to feed an emission model, focusing on the emission of nitrogen oxides (${\mathrm{NO_{x}}}$). Several studies deal with estimating emissions from dynamic traffic models, see for instance ^{[1,2,7,22,26,36,37]} and references therein. In particular, in ^[2] the authors deal with minimizing emissions by acting on the parameters of the model, while in ^[22] the authors analyze the possible benefits on emissions deriving from the limitation of traffic. The interest on ${\mathrm{NO_{x}}}$ gases in our work is due to their negative effects on health ^[40] and to their connection with ozone ^[4]. Minimizing only emissions would result in extreme solutions blocking traffic, thus we consider a cost function including a term measuring travel times. Therefore, we express the cost of emissions and travel time over the whole network as:

where $E^\gamma_r$, respectively $\mathcal{V}^{\gamma}_{r}$, is the emission rate, respectively velocity, along the road $r$, while $c_1$ and $c_2$ are weights. The functional ${\mathcal{F}}$ depends on the parameter vector $\gamma$ governing the traffic dynamic, which is comprised of the traffic distribution and priority parameters. Our interest is in minimizing ${\mathcal{F}}(\gamma)$ and compare different type of intersections, such as traffic lights and roundabouts. Note that, we consider a control $\gamma$ which varies in a compact set and therefore the minimum exists. Due the the high nonlinearity of ${\mathcal{F}}(\gamma)$, explicit analytical solutions can not be found in general. Therefore, we resort to numerical optimization to compute the optimal vectors $\gamma$. First, we focus on a merge junction and compare a priority-based junction with one regulated by a traffic light. The latter corresponds to alternating the values $\beta = 0$ and $\beta = 1$ for the green and red phases. These cycles are parameterized by the green-phase duration $t_g$ and the red-phase duration $t_r$. The numerical results show that it is possible to find an optimal $\beta$ and an optimal couple $(t_g, t_r)$, and that the two types of junctions perform similarly when minimizing emissions and travel time.

Next, we analyze how the solution to the minimization problem depends on the initial traffic state $(\rho, w)$. Here we interpret $w$ as drivers' preferred speed: low values of $w$ correspond to slow drivers, and high values of $w$ to fast drivers. For the priority-ruled junction, the minimum of the functional is achieved by giving high priority to the incoming road with higher density and fast drivers. Similarly, for the traffic light, the road with higher density must have a longer green-phase, except for high congestion when the opposite happens. In the latter situation, the sensitivity with respect to $w$ is greater.

We then focus on a more complex situation of a roundabout with two incoming and two outgoing roads. The roundabout has four additional stretch of roads to connect incoming to outgoing roads and form a circle. As before we compare priority-based junctions with traffic lights, by choosing optimally the priorities and the traffic light timing. Further, we compare the results obtained by the optimal traffic light timing with a periodic one showing that the former is convenient especially in reducing travel times. The numerical tests show that, when few vehicles enter the network, traffic lights produce lower emissions and travel times compared to the priority-based case. In congested situations, instead, the use of priorities produces higher levels of emissions but with shorter travel times w.r.t. traffic lights dynamics. It is worth to notice that traffic light timing can be easily adjusted in time, while changing priority-based rule would be more challenging. Overall, traffic lights outperform traffic signals in terms of emissions for roundabouts and perform better also taking into account travel times for low densities. Moreover, the optimal traffic light timing are more robust for variation of the functional weights. Interestingly, there is an increasing diffusion of roundabouts in Europe and US given the expected better performance in terms of output. This study shows that traffic signals should be added to roundabouts if one aims also at lowering emissions. This is a first example of how the model can be used to support decision makers for sustainable traffic management.

The paper is organized as follows. In Section 2 we define the GSOM and the Riemann problem at junctions. In Section 3 we describe the solution to the Riemann problem for diverge and merge junctions. In Section 4 a functional is formulated to estimate emission rate and travel time, while in Section 5 we provide details for the numerical approach. Sections 6 and 7 are devoted to the numerical tests for optimal controls and estimation of ${\mathrm{NO_{x}}}$ emissions. In Section 8 we draw our conclusions. Finally, in Appendix A we report some additional numerical tests for the roundabout.

2.   The Riemann Problem for GSOM at a junction

In order to extend the GSOM model to networks, one has to analyze the Riemann problem at a junction, i.e., the Cauchy problem with constant initial data on each road incident to the junction.

Recall the GSOM model equations (1.1). The variable $w$ parametrizes a family of fundamental diagrams $Q(\rho, w) = \rho V(\rho, w)$. The usual assumptions on $Q$ and $V$ are:

(H1) $Q(0, w) = 0$ and $Q({\rho^{\mathrm{max}}}(w), w) = 0$ for each $w$, where ${\rho^{\mathrm{max}}}(w)$ is the maximum density of vehicles for $Q(\cdot, w)$.

(H2) $Q(\rho, w)$ is strictly concave with respect to $\rho$, i.e., $\frac{ \partial^{2}Q}{ \partial\rho^{2}} < 0$.

(H3) $Q(\rho, w)$ is non-decreasing with respect to $w$, i.e., $Q_{w}\geq0$.

(H4) $V(\rho, w)\geq0$ for each $\rho$ and $w$.

(H5) $V(\rho, w)$ is strictly decreasing with respect to $\rho$, i.e., $V_{\rho} < 0$ for each $w$.

(H6) $V(\rho, w)$ is non-decreasing with respect to $w$, i.e., $V_{w}\geq0$.

From (H2) and (H3), for every $w$ the curve $\rho\to Q(\cdot, w)$ has a unique point of maximum, denoted by $\sigma(w)$, and we set ${Q^{\mathrm{max}}}(w) = Q(\sigma(w), w)$. Moreover, when $\rho = 0$ there is not a unique maximum velocity. For every $w$ we set ${V^{\mathrm{max}}}(w) = V(0, w)$.

The eigenvalues of Eq (1.1) are

The concavity of the flux implies $\lambda_1\leq\lambda_2$ and $\lambda_1 = \lambda_2$ if and only if $\rho = 0$, thus for $\rho\ne 0$ the system is strictly hyperbolic. The eigenvectors associated with the eigenvalues are

The first eigenvalue is genuinely nonlinear, i.e., $\nabla\lambda_1\cdot\gamma_1 \neq0$, while the second one is linearly degenerate, i.e., $\nabla\lambda_2\cdot\gamma_2 = 0$. Hence, the curves of the first family are 1-shocks or 1-rarefaction waves, while the curves of the second family are 2-contact discontinuities. Finally the Riemann invariants are

The first Riemann invariant ${z_{1}}$ is constant along 1-shock and 1-rarefaction waves, while the second Riemann invariant ${z_{2}}$ is constant along the 2-contact discontinuities.

By defining the total property $y = \rho w$, Eq (1.1) can be rewritten in conservative form as

where $v = V\Big(\rho, \frac{y}{\rho}\Big)$.

We recall now the main definitions concerning traffic models on road networks and we refer to ^{[13,18,20,25]} for further details. A road is modeled by an interval $I = (a, b)\subset \mathbb R$, with possibly $a = -\infty$ or $b = +\infty$. A junction $J$ is a collection of roads $((I_1, \ldots, I_n), (I_{n+1}, \ldots, I_{n+m}))$ where $I_1, \ldots, I_n$ are the incoming roads and $I_{n+1}, \ldots, I_{n+m}$ are the outgoing ones. We define a network as a couple $({\mathcal{I}}, {\mathcal{J}})$ where ${\mathcal{I}}$ is a finite collection of roads $I_{r}$, and ${\mathcal{J}}$ is a finite collection of junctions $J$.

On each road $I_{r}$, the traffic dynamic is described by a GSOM as

with $v_{r} = V\left(\rho_{r}, \frac{y_{r}}{\rho_{r}}\right)$, for $x\in I_{r}$ and $t\geq0$. The construction of a solution on the whole network is obtained via wave-front tracking starting from solutions to Riemann problems to Eq (2.3) at each junction. More precisely, given constant initial data on each road, we look for possible waves with negative speed for incoming roads and positive ones on outgoing roads. This is necessary to have conservation of mass through the junction, see ^[18]. In the literature this construction is usually called Riemann Solver, which is a map that assigns a solution to each constant initial data on the roads of the network, see [18,Definition 4.2.2] for further details. To isolate the admissible waves, we study the sign of the eigenvalues in Eq (2.1) and Eq (2.2). By the concavity of the flux function, the first eigenvalue ${\lambda_{1}}(\rho, w) = \rho+\rho{V_\rho}(\rho, w) = Q_{\rho}(\rho, w)$ satisfies ${\lambda_{1}}\geq0$ for $\rho\leq\sigma(w)$ and ${\lambda_{1}} < 0$ for $\rho > \sigma(w)$. The second eigenvalue is given by ${\lambda_{2}}(\rho, w) = V(\rho, w)$, thus by (H4) the speed of the 2-contact discontinuity is always non-negative.

In order to describe the flux maximization, let us consider the supply and demand functions, see ^[17] for details and discussion. The supply function $s({\rho}, {w})$ is defined as

and the demand function $d({\rho}, {w})$ as

2.1.   Incoming roads

Let us consider an incoming road at a junction. Only waves with negative speed are admissible. Since ${\lambda_{2}}\geq0$, we can have only 1-shock or 1-rarefaction waves.

We fix a left state ${U_{-}} = ({\rho^{-}}, {w^{-}})$ and look for the set of all admissible right states ${\hat U} = ({\hat\rho}, {\hat w})$ that can be connected to ${U_{-}}$ with waves with negative speed. Along the 1-waves the variable $w$ is conserved, therefore only the density $\rho$ changes. This case is analogous to the definition of admissible solutions on incoming roads for first order traffic models, see for instance ^[18]. We refer to ^[24] where the ARZ model is treated in a similar way.

Proposition 1. Let $V$ be a velocity function that verifies properties (H4)–(H6) and let ${U_{-}} = ({\rho^{-}}, {w^{-}})$ be a left state on an incoming road.

If ${\rho^{-}} = 0$, then the only admissible right state is ${\hat U} = {U_{-}}$.

If ${\rho^{-}}\neq0$, then the set of possible right states ${\hat U} = ({\hat\rho}, {\hat w})$ verifies ${\hat w} = {w^{-}}$ and:

1. If ${\rho^{-}}\leq\sigma({w^{-}})$, then ${\hat\rho}\in{\mathcal{N}}({U_{-}}) = \{{\rho^{-}}\}\cup({\widetilde\rho^{-}}({w^{-}}), {\rho^{\mathrm{max}}}({w^{-}})]$, where ${\widetilde\rho^{-}}({w^{-}})$ is the density such that $Q({\widetilde\rho^{-}}({w^{-}}), {w^{-}}) = Q({\rho^{-}}, {w^{-}})$.

2. If ${\rho^{-}} > \sigma({w^{-}})$, then ${\hat\rho}\in{\mathcal{N}}({U_{-}}) = [\sigma({w^{-}}), {\rho^{\mathrm{max}}}({w^{-}})]$.

Moreover, denoting by $d$ the demand function defined in Eq (2.5), it holds

Proof. First assume ${\rho^{-}}\neq0$. If ${\rho^{-}}\leq\sigma({w^{-}})$ (Figure 1 top-left) to have $\lambda_1\leq 0$ there are two possibilities: either ${\hat U} = {U_{-}}$, or moving above the density value ${\widetilde\rho^{-}}({w^{-}}) > \sigma({w^{-}})$ by a jump with zero speed. Indeed, since $Q({\widetilde\rho^{-}}({w^{-}}), {w^{-}}) = Q({\rho^{-}}, {w^{-}})$, the Rankine-Hugoniot condition $s({\widetilde U^{-}}-{U_{-}}) = Q({\widetilde\rho^{-}}({w^{-}}), {w^{-}})-Q({\rho^{-}}, {w^{-}})$ implies that the speed of the discontinuity $s$ is zero. In this case we can move with a 1-shock with negative speed towards any right state ${\hat U}$ with ${\hat w} = {w^{-}}$ and ${\widetilde\rho^{-}}({w^{-}}) < {\hat\rho}\leq{\rho^{\mathrm{max}}}({w^{-}})$. If ${\rho^{-}} = 0$ then ${\widetilde\rho^{-}}({w^{-}}) = {\rho^{\mathrm{max}}}({w^{-}})$, therefore the solution is ${\hat U} = {U_{-}}$.

If ${\rho^{-}} > \sigma({w^{-}})$, every state ${\hat U}$ with ${\hat w} = {w^{-}}$ and ${\hat\rho}\in[\sigma({w^{-}}), {\rho^{\mathrm{max}}}({w^{-}})]$ is connected to ${U_{-}}$ with waves with negative speed (Figure 1 top-right). In particular, we have a 1-rarefaction wave if ${\hat\rho}\leq{\rho^{-}}$ and a 1-shock if ${\hat\rho} > {\rho^{-}}$.

2.2.   Outgoing roads

Let us consider an outgoing road at a junction. We are interested in the waves with positive speed, thus we can have a 1-shock or 1-rarefaction wave and a 2-contact discontinuity.

We fix a right state ${U_{+}} = ({\rho^{+}}, {w^{+}})$ and look for the set of all admissible left states ${\hat U} = ({\hat\rho}, {\hat w})$ that can be connected to ${U_{+}}$ with waves with positive speed. We emphasize that along the 1-waves the $w$ is conserved and only the density $\rho$ changes. We therefore assume that the value ${w^{-}}$ is given and depends on the states on the incoming roads (see Section 3). On the other hand, along the 2-wave the velocity $V(\rho, w)$ is conserved. Then, the definition of the admissible states ${\hat U}$ depends on the existence of an intermediate point ${U^{\dagger}} = ({\rho^{\dagger}}, {w^{\dagger}})$ such that ${w^{\dagger}} = {w^{-}}$ and $V({\rho^{\dagger}}, {w^{\dagger}}) = V({\rho^{+}}, {w^{+}})$.

Proposition 2. Let $V$ be a velocity function that verifies properties (H4)–(H6). For a given value ${w^{-}}$ and a given right state ${U_{+}} = ({\rho^{+}}, {w^{+}})$ with associated velocity ${v^{+}} = V({\rho^{+}}, {w^{+}})$, if ${v^{+}}\leq{V^{\mathrm{max}}}({w^{-}})$ then there exists a unique point ${U^{\dagger}} = ({\rho^{\dagger}}, {w^{\dagger}})$ such that ${w^{\dagger}} = {w^{-}}$ and $V({\rho^{\dagger}}, {w^{\dagger}}) = {v^{+}}$.

Proof. If ${v^{+}}\leq{V^{\mathrm{max}}}({w^{-}})$ then the equation $V(\rho, {w^{-}}) = {v^{+}}$ admits a solution. By (H5), $\partial_\rho V < 0$ and, by the implicit function theorem, there exists $\rho(w; {v^{+}})$ such that $V(\rho(w; {v^{+}}), w) = {v^{+}}$. Moreover, (H5)–(H6) imply

We then have ${w^{\dagger}} = {w^{-}}$ and ${\rho^{\dagger}} = \rho({w^{-}}; {v^{+}})$.

Proposition 3. Let $V$ be a velocity function that verifies properties (H4)–(H6), ${U_{+}} = ({\rho^{+}}, {w^{+}})$ a right state on an outgoing road, and ${v^{+}} = V({\rho^{+}}, {w^{+}})$ the associated velocity. For a given value ${w^{-}}$, a left state ${\hat U} = ({\hat\rho}, {\hat w})$, which can be connected to ${U_{+}}$ with positive speed waves, satisfies ${\hat w} = {w^{-}}$ and the following.

(i) If $v^+ \leq {V^{\mathrm{max}}}({w^{-}})$, let ${U^{\dagger}} = ({\rho^{\dagger}}, {w^{\dagger}})$ be the intersection point between the level curves $\{{z_{2}} = {v^{+}}\}$ and $\{{z_{1}} = {w^{-}}\}$, then ${w^{\dagger}} = {w^{-}}$ and

1. if ${\rho^{\dagger}}\leq\sigma({w^{\dagger}})$, then ${\hat\rho}\in{\mathcal{P}}({U_{+}}) = [0, \sigma({w^{\dagger}})]$;

2. if ${\rho^{\dagger}} > \sigma({w^{\dagger}})$, then ${\hat\rho}\in{\mathcal{P}}({U_{+}}) = [0, {\widetilde\rho^{\dagger}}({w^{\dagger}}))\cup\{{\rho^{\dagger}}\}$, where ${\widetilde\rho^{\dagger}}({w^{\dagger}})$ is the density such that $Q({\widetilde\rho^{\dagger}}({w^{\dagger}}), {w^{\dagger}}) = Q({\rho^{\dagger}}, {w^{\dagger}})$.

(ii) If $v^+ > {V^{\mathrm{max}}}({w^{-}})$ then ${\hat\rho}\in{\mathcal{P}}({U_{+}}) = [0, \sigma({w^{-}})]$.

Moreover, denoting by $s$ the supply function defined in Eq (2.4), it holds

Proof. If ${v^{+}}\leq{V^{\mathrm{max}}}({w^{-}})$, by Proposition 2 there exists a unique point ${U^{\dagger}}$ such that ${w^{\dagger}} = {w^{-}}$ and $V(\rho, {w^{\dagger}}) = {v^{+}}$. Thus, if ${\rho^{\dagger}}\leq\sigma({w^{\dagger}})$, then every state ${\hat U}$ with ${\hat w} = {w^{-}}$ and ${\hat\rho}\in[0, \sigma({w^{\dagger}})]$ can be connected to ${U^{\dagger}}$ by waves with positive speed (Figure 1 bottom-left). In particular we have a 1-rarefaction wave if ${\rho^{\dagger}}\leq{\hat\rho}$ and a 1-shock if ${\rho^{\dagger}} > {\hat\rho}$. Then, ${U^{\dagger}}$ is connected to ${U_{+}}$ by a 2-contact discontinuity which has positive speed.

If ${\rho^{\dagger}} > \sigma({w^{\dagger}})$ (Figure 1 bottom-right), we have two possibilities: no wave, then ${\hat U} = {U^{\dagger}}$, or moving below the density value ${\widetilde\rho^{\dagger}}({w^{\dagger}}) < \sigma({w^{\dagger}})$ by a jump with positive speed. In this case, a 1-rarefaction connects to an intermediate state ${\hat U}$ with ${\hat w} = {w^{\dagger}}$ and $0\leq{\hat\rho}\leq{\widetilde\rho^{\dagger}}({w^{\dagger}})$, then a 2-contact discontinuity connects to ${U_{+}}$.

Otherwise, if ${v^{+}} > {V^{\mathrm{max}}}({w^{-}})$ then the equality $V(\rho, {w^{-}}) = {v^{+}}$ can not hold. It holds ${\rho^{\dagger}} = 0$ and the admissible left state ${\hat\rho}$ has to be in $[0, \sigma({w^{-}})]$.

To summarize, we denote

where $\rho(\cdot; {v^{+}})$ is the implicit function given by the equation $V(\rho, w) = {v^{+}}$, which is well defined as stated in Proposition 2.

Remark 1. For numerical purposes, we use the Collapsed Generalized Aw-Rascle-Zhang (CGARZ) model, see ^[16] and Section 5. This model is characterized by a maximum velocity ${V^{\mathrm{max}}}$ common to any $w$. Hence, the case of ${v^{+}} > {V^{\mathrm{max}}}(w)$ never holds for the CGARZ model.

3.   The GSOM on networks

In this section we apply Propositions 1 and 3 to define the Riemann Solver for merge and diverge junctions. To identify a unique solution we assume the maximization of the flux and the conservation of $\rho$ and $y = \rho w$ across the junction. Moreover, we assume that a distribution parameter on outgoing roads and a priority rule on incoming ones are given.

3.1.   Diverge junction

We consider the case of a junction with one incoming and two outgoing roads. Given a left state ${U_{-}}_{1}$ for the incoming road and two right states ${U_{+}}_{2}$ and ${U_{+}}_{3}$ for the outgoing roads, our aim is to determine the junction values ${\hat U}_{i} = (\hat\rho_i, \hat w_i)$, $i = 1, 2, 3$, giving rise to a boundary-value problem on each road. The solutions to the latter pieced together provide a solution to the Riemann problem at the junction. First, introduce a traffic distribution parameter $\alpha\in(0, 1)$: vehicles are distributed in proportion $\alpha$ and $1-\alpha$ on the roads 2 and 3, respectively. Note that the cases $\alpha = 0$ or $\alpha = 1$ reduce the problem to a simple 1 to 1 junction, thus in this analysis we exclude the two extreme values. Set $\hat q_i = \hat \rho_i \hat v_i$, $\hat v_i = V(\hat\rho_i, \hat w_i)$, $i = 1, 2, 3$, then the conservation of $\rho$ and $y$ across the junction reads:

By Proposition 1 we have ${\hat w}_{1} = {w^{-}}_{1}$, and by Eqs (3.1)–(3.4) we deduce ${\hat w}_{1} = {\hat w}_{2}$ and ${\hat w}_{1} = {\hat w}_{3}$, hence ${\hat w}_{2} = {\hat w}_{3} = {w^{-}}_{1}$. Now the states ${\hat U}_{i}$ correspond to six unknowns for which we have five equations. Using the free parameter $q = \hat q_1$ and, by Eqs (2.6) and (2.7) we get the constraints

where, by Proposition 3, ${w^{\dagger}}_{2} = {w^{\dagger}}_{3} = {w^{-}}_{1}$ and ${\rho^{\dagger}}_{2}$, ${\rho^{\dagger}}_{3}$ are given by Eq (2.8) with ${w^{-}} = {w^{-}}_i$ and ${v^{+}} = {v^{+}}_i$, $i = 2, 3$, respectively. To satisfy Eq (3.5) and maximize the outgoing flux, it holds

and

Then, the junction density values are ${\hat\rho}_{1}\in{\mathcal{N}}({U_{-}}_{1})$ such that $Q({\hat\rho}_{1}, {w^{-}}_1) = {\hat q}_{1}$ and ${\hat\rho}_{j}\in{\mathcal{P}}({U_{+}}_{j})$ such that $Q({\hat\rho}_{j}, {w^{-}}_{j}) = {\hat q}_{j}$, $j = 2, 3$. In ^[24,28], the authors obtain the same solution for the ARZ model.

3.2.   Merge junction

We consider the case of a junction with two incoming and one outgoing roads. Given two states ${U_{-}}_{1}$ and ${U_{-}}_2$ for the incoming roads and a state ${U_{+}}_{3}$ for the outgoing road, we look for the junction values ${\hat U}_{1}$, ${\hat U}_{2}$ and ${\hat U}_{3}$. As done before, we set $\hat q_i = \hat\rho_i \hat v_i$, $i = 1, 2, 3$, and we assume that vehicles from roads 1 and 2 enter into the road 3 with the following priority rule

where $\beta\in[0, 1]$. Note that for $\beta = 0$ or $\beta = 1$, one of the two incoming roads is completely stopped at the junction, and the problem reduces to the 1 to 1 case.

The conservation of $\rho$ and $y$ across the junction yields:

By Proposition 1, we have that ${\hat w}_{1} = {w^{-}}_{1}$ and ${\hat w}_{2} = {w^{-}}_{2}$. Equation (3.7) combined with Eqs (3.6) and (3.8), implies

Hence, ${\hat w}_1$, ${\hat w}_2$ and ${\hat w}_3$ are defined and ${\hat q}_1$, ${\hat q}_2$ and ${\hat q}_3$ have to satisfy Eqs (3.6) and (3.7). It remains a free parameter and, in order to define a unique solution, we impose the maximization of the flux on the outgoing road. By Eqs (2.6) and (2.7), we get the constraints

where ${\rho^{\dagger}}_{3}$ is given by Eq (2.8) with ${w^{-}} = {\hat w}_3$ and ${v^{+}} = V({\rho^{+}}_3, {w^{+}}_3)$. From now on, we set

We assume that both ${d_{1}}$ and ${d_{2}}$ are greater than 0. Indeed, the trivial case of ${d_{1}} = {d_{2}} = 0$ means that no vehicles cross the intersection, and the case of ${d_{1}} = 0$ or ${d_{2}} = 0$ reduces the junction to the 1 to 1 type. In order maximize the flux on the outgoing road we set in Eq (3.10)

To summarize, the couple $(\hat{q_{1}}, \hat{q_{2}})$ is given by the intersection point $P$ between the following two lines

where the first one represents the priority rule Eq (3.6), while the second one represent the conservation equation (3.7) coupled with Eq (3.11). In Eq (3.12), ${\mathrm{r}}$ coincides with the axis $x = 0$ when $\beta = 1$. Note that, since ${\rho^{\dagger}}_{3} = {\rho^{\dagger}}_{3}({\hat w}_{3};{v^{+}}_3)$ and ${\hat w}_{3}$ depends on $\beta$, the maximum flux that can be received by the outgoing road is a function of the priority rule, i.e., $s({{\rho^{\dagger}}_{3}}, {{\hat w}_{3}}) = {s_{3}}(\beta)$.

The intersection point between $\mathrm{r}$ and $\mathrm{s}$ is

If $P\in\Omega = [0, d_{1}]\times[0, d_{2}]$, we can set ${\hat q}_{1} = (1-\beta){s_{3}}(\beta)$ and ${\hat q}_{2} = \beta{s_{3}}(\beta)$. Otherwise, if $P\notin\Omega$, then the point does not satisfy the constraints (3.10), and we need to relax one of our constraints. We propose two possible approaches:

(RP) The relation Eq (3.6) is satisfied with $\beta$ fixed a priori, while the outgoing flow in Eq (3.11) is not maximized. This is the case for instance of a stop sign or a traffic policeman regulating the junction.

(AP) The priority parameter $\beta$ is modified, thus allowing to maximize the outgoing flux. This is the case, for instance, of unsupervised junction.

To detail the procedure to compute the junction densities $\hat \rho_i$, $i = 1, 2, 3$, first recall that $\hat w_1 = \hat w_2 = {w^{-}}_1$ and $\hat w_3 = (1-\beta){w^{-}}_1+\beta{w^{-}}_2$ as stated in Eq (3.9). We introduce the parameter

which identifies the priority line in Eq (3.12) that passes through the point $({d_{1}}, {d_{2}})$. If $P\notin\Omega$, we distinguish two cases:

(i) $\beta\geq{\beta_{d}}$ then the $y$-coordinate of $P$, $\beta{s_{3}}(\beta)$ is greater than the upper bound ${d_{2}}$. Then we fix ${\hat q}_2 = {d_{2}}$ and look for an admissible value ${\hat q}_1$;

(ii) $\beta < {\beta_{d}}$ then the $x$-coordinate $(1-\beta){s_{3}}(\beta) > {d_{1}}$. Then we fix ${\hat q}_1 = {d_{1}}$ and look for an admissible value ${\hat q}_2$.

We first describe the RP algorithm. For a given priority parameter ${\bar\beta}\in[0, 1]$, to satisfy the priority rule, the solution must lie on the line $(1-{\bar\beta}){\hat q}_2 = {\bar\beta}{\hat q}_1$. For this reason, when $P\notin\Omega$ the couple $({\hat q}_{1}, {\hat q}_{2})$ will be defined by the intersection point between the priority line and the boundary $\partial\Omega$, see for instance the point $Q$ in Figures 2b and 2c.

Definition 3.1. Riemann Solver Algorithm RP. Let ${\bar\beta}\in[0, 1]$ and let ${\hat w}_{3}$ and $P$ be as in Eq (3.9) and Eq (3.14) with $\beta = {\bar\beta}$, respectively. Assume that Eq (3.7) holds. Define $(\hat q_1, \hat q_2)\in\Omega$ as follows:

1. If $P\in\Omega$, then ${\hat q}_{1} = (1-{\bar\beta}){s_{3}}({\bar\beta})$ and ${\hat q}_{2} = {\bar\beta}{s_{3}}({\bar\beta})$.

2. If $P\notin\Omega$ and ${\bar\beta}\geq{\beta_{d}}$, then ${\hat q}_{1} = (1-{\bar\beta}){d_{2}}/{\bar\beta}$ and ${\hat q}_{2} = {d_{2}}$.

3. If $P\notin\Omega$ and ${\bar\beta} < {\beta_{d}}$, then ${\hat q}_{1} = {d_{1}}$ and ${\hat q}_{2} = {\bar\beta}{d_{1}}/(1-{\bar\beta})$.

The density value ${\hat\rho}_{i}\in{\mathcal{N}}({U_{-}}_{i})$ is determined by the equality $Q({\hat\rho}_{i}, {\hat w}_{i}) = {\hat q}_{i}$, $i = 1, 2$, while ${\hat\rho}_{3}\in{\mathcal{P}}({U_{+}}_{3})$ is determined by $Q({\hat\rho}_{3}, {\hat w}_{3}) = {\hat q}_{3}$.

To describe the AP algorithm we need some preliminary results.

Lemma 3.2. The supply function $s({\rho^{\dagger}}_3(w; \cdot), w)$ is non-decreasing in $w$.

Proof. By Eq (2.4) we can have

In the first case, assumption (H3) applies; in the second one, by Proposition 2 we have $V({\rho^{\dagger}}_3(w; \cdot), w) = {v^{+}}_3$ and ${\rho^{\dagger}}_3(w; \cdot)$ is non-decreasing in $w$.

To study the function ${s_{3}}(\beta) = s({{\rho^{\dagger}}_{3}({\hat w}_{3}(\beta); \cdot)}, {{\hat w}_{3}(\beta)})$ with respect to $\beta$, we distinguish two cases:

(a) ${w^{-}}_1 \leq {w^{-}}_2$ then both ${\hat w}_3$ and ${s_{3}}$ are increasing in $\beta$;

(b) ${w^{-}}_1 > {w^{-}}_2$ then both ${\hat w}_3$ and ${s_{3}}$ are decreasing in $\beta$.

Lemma 3.3. Let ${\bar\beta}\in[0, 1]$ and ${\beta_{d}}$ given in Eq (3.15).

1. If ${\bar\beta}\geq{\beta_{d}}$ and ${\bar\beta} {s_{3}}({\bar\beta}) > {d_{2}}$, then there exists at least a $\beta\in[0, {\bar\beta})$ such that $\beta{s_{3}}(\beta) = {d_{2}}$.

2. If ${\bar\beta} < {\beta_{d}}$ and $(1-{\bar\beta}){s_{3}}({\bar\beta}) > {d_{1}}$, then there exists at least a $\beta\in({\bar\beta}, 1]$ such that $(1-\beta){s_{3}}(\beta) = {d_{1}}$.

Proof. We first prove point $1$. Consider the two cases (a) and (b), i.e. ${w^{-}}_1 < {w^{-}}_2$ and ${w^{-}}_1 > {w^{-}}_2$, respectively.

If ${w^{-}}_1 < {w^{-}}_2$ then the function $f(\beta) = \beta{s_{3}}(\beta)$ is increasing in $[0, 1]$ and such that $f(0) = 0$ and $f({\bar\beta}) > {d_{2}}$ by hypothesis; therefore, there exists a unique $\beta^* < {\bar\beta}$ such that $f(\beta^*) = {d_{2}}$.

If ${w^{-}}_1 > {w^{-}}_2$ then ${s_{3}}(\beta)$ is decreasing in $\beta$ and the behavior of the function $f(\beta) = \beta{s_{3}}(\beta)$ is not known a priori. However the function $f$ is continuous and such that $f(0) = 0$ and $f({\bar\beta}) > {d_{2}}$ by hypothesis; therefore there exists at least a $\beta < {\bar\beta}$ such that $f(\beta) = {d_{2}}$.

The proof of point $2$ is entirely similar, so we skip the details.

The AP algorithm is described in the following definition. As mentioned above, the algorithm adapts the priority parameter to maximize the outgoing flux while keeping the parameter as close as possible to its initial value.

Definition 3.4. Riemann Solver Algorithm AP. Let ${\bar\beta}\in[0, 1]$ and $P$ be as in Eq (3.14) with $\beta = {\bar\beta}$. Assume that Eq (3.7) holds. Define $(\hat q_1, \hat q_2)\in\Omega$ as follows:

1. If $P\in\Omega$ then ${\hat q}_{1} = (1-{\bar\beta}){s_{3}}({\bar\beta})$, ${\hat q}_{2} = {\bar\beta}{s_{3}}({\bar\beta})$ and ${\hat w}_{3} = (1-{\hat\beta}){w^{-}}_{1}+{\hat\beta}{w^{-}}_{2}$ with ${\hat\beta} = {\bar\beta}$.

2. If $P\notin\Omega$ and ${\bar\beta}\geq{\beta_{d}}$, then for $\beta^{*} = \max\{\beta\in[0, {\bar\beta})\, :\, \beta{s_{3}}(\beta) = {d_{2}}\}$, we set ${\hat\beta} = \max\{\beta^{*}, {\beta_{d}}\}$, ${\hat q}_{1} = \min\{(1-{\hat\beta}) {s_{3}}({\hat\beta}), {d_{1}}\}$, ${\hat q}_{2} = {d_{2}}$ and ${\hat w}_{3} = (1-{\hat\beta}){w^{-}}_{1}+{\hat\beta}{w^{-}}_{2}$.

3. If $P\notin\Omega$ and ${\bar\beta} < {\beta_{d}}$, then for $\beta^{*} = \min\{\beta\in({\bar\beta}, 1] \, :\, (1-\beta){s_{3}}(\beta) = {d_{1}}\}$, we set ${\hat\beta} = \min\{\beta^{*}, {\beta_{d}}\}$, ${\hat q}_{1} = {d_{1}}$, ${\hat q}_{2} = \min\{{\hat\beta}{s_{3}}({\hat\beta}), {d_{2}}\}$ and ${\hat w}_{3} = (1-{\hat\beta}){w^{-}}_{1}+{\hat\beta}{w^{-}}_{2}$.

The density value ${\hat\rho}_{i}\in{\mathcal{N}}({U_{-}}_{i})$ is determined by the equality $Q({\hat\rho}_{i}, {\hat w}_{i}) = {\hat q}_{i}$, $i = 1, 2$, while ${\hat\rho}_{3}\in{\mathcal{P}}({U_{+}}_{3})$ is determined by $Q({\hat\rho}_{3}, {\hat w}_{3}) = {\hat q}_{3}$.

Proposition 4. The couple $({\hat q}_{1}, {\hat q}_{2})$ in Definition 3.4 satisfies the constraints in Eq (3.10).

Proof. The given couple $({\hat q}_1, {\hat q}_2)$ verifies the first two constraints in Eq (3.10) by construction. Therefore, it remains to prove that ${\hat q}_3 = {\hat q}_{1}+{\hat q}_{2}\leq{s_{3}}({\hat\beta}) = {s_{3}}({\rho^{\dagger}}_{3}({\hat w}_3;\cdot), {\hat w}_{3})$.

We start from the case 2 of Definition 3.4. In light of Lemma 3.3 case 1, the value $\beta^{*}$ is well defined. Moreover, since the slope of r increases with $\beta$, the point $((1-\beta^*){s_{3}}(\beta^*), {d_{2}})$ is such that: if $\beta^*\geq {\beta_{d}}$ then $(1-\beta^*){s_{3}}(\beta^*)\leq d_1$ and if $\beta^* > {\beta_{d}}$ then $(1-\beta^*){s_{3}}(\beta^*) > d_1$. Therefore, we focus on these two possibilities:

● If $\beta^*\geq{\beta_{d}}$ then ${\hat\beta} = \beta^*$ and ${\hat q}_1 = (1-{\hat\beta}){s_{3}}({\hat\beta})$. Hence, ${\hat q}_1+{\hat q}_2 = (1-{\hat\beta}){s_{3}}({\hat\beta}) + {d_{2}} = (1-{\hat\beta}){s_{3}}({\hat\beta}) + {\hat\beta}{s_{3}}({\hat\beta}) = {s_{3}}({\hat\beta})$ and the thesis follows. This is the case, for instance, of point $R$ in Figure 2b.

● If $\beta^* < {\beta_{d}}$ then ${\hat\beta} = {\beta_{d}}$ and ${\hat q}_1 = {d_{1}}$. The couple $({\hat q}_{1}, {\hat q}_{2}) = ({d_{1}}, {d_{2}})$ is admissible if ${d_{1}}+{d_{2}}\leq{s_{3}}({\beta_{d}})$. Since ${d_{1}} = (1-{\beta_{d}}){d_{2}}/{\beta_{d}}$ we have ${d_{1}}+{d_{2}} = {d_{2}}/{\beta_{d}}$. From the definition of $\beta^*$, for each $\beta\in(\beta^{*}, {\bar\beta}]$ it holds $\beta{s_{3}}(\beta) > {d_{2}}$, and we get the thesis:

This is the case of point $S$ in Figure 2b.

The proof of case 3 follows similarly, see Figure 2c for an example of possible configuration.

Remark 2. Let us consider the particular case of ${w^{-}}_1 = {w^{-}}_2 = {w^{+}}_3$, i.e., the variable $w$ is constant on the roads network. The diverge junction can be treated exactly as the LWR model at junctions, as done in ^[18]. For the merge junction we observe that the assumption of $w$ constant implies that the straight line ${\mathrm{s}}$ defined in Eq (3.13) coincides for all $\beta$, therefore the solution is limited to the points $P$ or $Q$ in Figure 2, excluding the points $R$ and $S$. Thus, we recover again the LWR model on networks, as treated in ^[18].

4.   Minimize emissions and travel time

The emission of pollutants is strictly connected to speed and acceleration of vehicles. In this section we set up an optimization problem to minimize the ${\mathrm{NO_{x}}}$ emission rates due to vehicular traffic.

Consider Eq (2.3) on a network with roads $I_{r}$, $r = 1, \ldots, N_r$, during a time interval $[0, T]$. Following ^[7], we use the microscopic emission model proposed in ^[32] which estimates the emission rate $E_{i}$ of vehicle $i$ at time $t$ using the instantaneous speed $v_i(t)$ and acceleration $a_i(t)$. We then define

where $E_0$ is a lower-bound of emission and $f_1$ to $f_6$ are emission constants associated with ${\mathrm{NO_{x}}}$, see Table 1.

Let $\Gamma\subset \mathbb R^{k}$ be the set of $k$ control parameters $\gamma = (\gamma_1, \ldots, \gamma_k)$ governing the traffic dynamic. These are given by the traffic distribution and priority parameters $\alpha$ and $\beta$ of Section 3. Let ${N_{r}}$ be the number of roads and $E^{\gamma}_{r}(x, t)$ be the emission rate in $x$ at time $t$ related to $\gamma$ and road $r$. Note that $E^{\gamma}_{r}(x, t)$ is estimated through Eq (4.1) by summing all the vehicle contributions in $x$ at time $t$. We introduce the following operator to estimate the total emission rate on a road network as a function of $\gamma\in\Gamma$,

To guarantee acceptable travel times, we include a velocity term thus getting the objective function

where $c_1$ and $c_2$ are two proper weights and $\mathcal{V}^{\gamma}_{r} = \max\{V^{\gamma}_{r}(x, t), \varepsilon\}$, $\varepsilon > 0$, with $V^\gamma_r$ velocity function of the traffic model, related to control parameter $\gamma$ and to road $r$. The parameter $\varepsilon$ allows to exclude the null speeds in the calculation. Our goal is to solve the minimization problem

The minimum exists since $\Gamma$ is a compact set in $\mathbb R^{k}$, but we do not expect uniqueness. Due to the complexity and the strictly nonlinear dependence of the functional ${\mathcal{F}}$ on the control $\gamma$, we treat the problem numerically using global search.

5.   Numerical setup

From the GSOM family, we choose the CGARZ model ^[16] for simulations. The CGARZ model assumes that there is a unique maximum density ${\rho^{\mathrm{max}}}$ independent of $w$ at which the vehicles stop, i.e. $V({\rho^{\mathrm{max}}}, w) = 0$ for all $w$. Furthermore, it assumes a given free-flow threshold density $\rho_f$ such that the flux of vehicles is not influenced by $w$ when $\rho\leq\rho_f$ (free-flow regime). Thus, the flux is described by a single-valued fundamental diagram in free-flow regimes and by a multi-valued function in congestion. For $\rho\in[0, {\rho^{\mathrm{max}}}]$, we have

Following ^[7], we assume a lower and upper bound for $w$, i.e., $0\leq {w_{L}}\leq w\leq {w_{R}}$, a Greenshields flux function in the free-flow phase, i.e.,

and a flux in congested phase given by

where ${w_{L}} = Q_{f}(\rho_{f})$, ${w_{R}} = Q_{f}({\rho^{\mathrm{max}}}/2)$ and ${\rho^{\mathrm{max}}}/2$ is the critical density of $Q_{f}(\cdot)$. The velocity function is then given by

With these choices, the property $w$ describes drivers attitude with respect to speed. Low values of $w$ describe slow drivers, and high values of $w$ fast drivers.

We consider now the traffic model in Eq (2.3) with flux function given in Eq (5.1), and we divide each road into ${N_{x}}$ cells $[x_{j-1/2}, x_{j+1/2})$ of length ${\Delta x}$ centered in $x_{j}$, and the time interval into ${N_{t}}+1$ steps $t^{n} = n{\Delta t}$.

The model is then solved numerically with the 2CTM scheme ^[7] with suitable boundary conditions at the extremes of the network. We use the theory given in Sections 3.1 and 3.2 to build the numerical solution at junctions.

The 2CTM numerical scheme is described by the two equations

where $F^{\rho, n}_{r, j\pm1/2}$ and $F^{y, n}_{r, j\pm1/2}$ are the numerical fluxes, $\rho^{n}_{j}$ and $y^{n}_{j}$ are the $j$-cell average

for any time $t^{n}$, respectively. For the CGARZ model, the numerical fluxes are defined as

where $d({\cdot}, {\cdot})$ and $s({\cdot}, {\cdot})$ are the demand and supply functions defined in Eqs (2.5) and (2.4) respectively. Since $y = \rho w$ the numerical fluxes $F^{y, n}_{r, j\pm1/2}$ are such that

Moreover, the stability of the scheme is guaranteed by the CFL condition

where for the CGARZ model ${V^{\mathrm{max}}} = V(0, w)$ for all $w$. See ^[7] and references therein for more details.

The emission model (4.1) is based on vehicles speed and acceleration. The latter is obtained by computing the total derivative of $V(\rho, w)$, i.e.,

where

By simple computations, for the GSOM we have

Once $\rho^{n}_{r, j}$, $w^n_{r, j}$ and $V^{n}_{r, j}$ are known, from Eq (5.6) we get the discrete acceleration

Then, we compute the emission rate by adapting equation (4.1) to the numerical framework. Hence, for each time $t^{n}$ and cell $x_j$, $j = 1, \ldots, {N_{x}}$, of road $r$ we define

where $E_0 = 0$ and the coefficients $f_1$ to $f_6$ are collected in Table 1.

The functional ${\mathcal{F}}(\gamma)$ in Eq (4.3) is then discretized as

where $E^{\max}$ is the maximum emission rate, $\varepsilon$ is the rounded minimum velocity, and, in order to have comparable quantities for the emission and travel time functional, the weights $c_{1}$ and $c_{2}$ are given by

From now on we assume $\varepsilon = 1\, {\mathrm{km/h}}$. As shown in Appendix A, this choice of weights does not substantially affect the numerical results described in the following sections. Thus ${\mathcal{F}}$ in Eqs (4.3) and (5.9) is an appropriate functional to analyze the cost in emission and travel time.

6.   Case study of a merge junction

Let us consider the merge junction depicted in Figure 3, where we assume road 1 to be a ramp merging to roads 2 and 3. We assume the junction to be governed first by a priority rule and then by a traffic light. The latter is modeled by alternating $\beta = 0$ and $\beta = 1$ in time.

The model parameters in Eq (5.1) and those for the numerical tests are fixed in Table 2. The initial data is assumed to be constant on all the three roads and is chosen according to Table 3.

First we focus on the emission functional

We look for the parameter $\beta\in[0, 1]$ which minimizes ${\mathcal{F}_{E}}(\beta)$, and analyze the RP algorithm given in Definition 3.1. In Figure 4 on the left, we show ${\mathcal{F}_{E}}(\beta)$ for $\beta$ varying in $[0, 1]$. The optimal priority rule is given by $\beta^{opt} = 0$, i.e., no vehicle enters the junction from road 2. This result is unrealistic and thus motivates the use of the extended functional (4.3) including travel times. In Figure 4 on the right, we then show the test result for functional ${\mathcal{F}}$ with $\varepsilon = 1\, {\mathrm{km/h}}$. The optimal parameter is $\beta^{opt} = 0.64$ and ${\mathcal{F}}(\beta^{opt}) = 8.10$ which is an admissible and realist solution.

Optimal traffic light We model a traffic light placed at the end of roads 1 and 2 (see Figure 3b) by alternating $\beta = 0$ and $\beta = 1$ in time. Specifically, for $\beta = 0$ the traffic light is green for road 1 and red for road 2, on the contrary for $\beta = 1$ it is red for road 1 and green for road 2. The controls are given by the green phase duration $t_{g}$ (when $\beta = 0$) and red phase duration $t_{r}$ (when $\beta = 1$). The problem (4.4) is studied for $\Gamma = G \times R\subset \mathbb R^{2}$, where $G$ and $R$ are the intervals where $t_g$ and $t_r$ vary, and the cost functional ${\mathcal{F}}(\gamma) = {\mathcal{F}}(t_{g}, t_{r})$. Fixing $G = R = [0, 90\, {\mathrm{s}}]$, in Figure 5 we plot ${\mathcal{F}}(t_{g}, t_{r})$ with initial traffic data given in Table 3. The optimal times are $t^{opt}_{g} = 5\, {\mathrm{s}}$ and $t^{opt}_{r} = 10\, {\mathrm{s}}$ and ${\mathcal{F}}(t^{opt}_{g}, t^{opt}_{r}) = 8.18$. We observe that the region bounded by dark-blue lines identifies the points with functional values close to the minimum one. Therefore, many couples $(t_g, t_r)$ allow to have low emissions and travel time.

In summary, in Table 4 we compare the minimum values of ${\mathcal{F}_{E}}(\gamma)$, ${\mathcal{F}_{T}}(\gamma)$ and ${\mathcal{F}}(\gamma)$ obtained with $\gamma = \beta^{opt}$ and $\gamma = (t^{opt}_{g}, t^{opt}_{r})$. The optimal values are very close. The numerical tests show that ${\mathcal{F}}$ has a convex shape, both with respect to $\beta^{opt}$ and $(t^{opt}_{g}, t^{opt}_{r})$.

6.1.   Sensitivity to initial data

Here we investigate numerically the sensitivity of the minimization problem (4.4) with respect to the initial traffic states for constant initial data on all three roads. We consider two different traffic scenarios:

(i) $\rho^0_{2, 3} < {\rho_{f}}$, i.e. free flow traffic conditions on roads 2 and 3. Specifically, we fix $\rho^0_{2} = \rho^0_{3} = 15\, {\mathrm{veh/km}}$ and $w^0_{2} = w^0_{3} = ({w_{L}}+{w_{R}})/2$ along the roads;

(ii) $\rho^0_{2, 3} > {\rho_{f}}$, i.e. congested traffic conditions on roads 2 and 3. Specifically, we fix $\rho^0_{2} = \rho^0_{3} = 60\, {\mathrm{veh/km}}$ and $w^0_{2} = w^0_{3} = ({w_{L}}+{w_{R}})/2$ along the roads.

The optimal control is computed as function of the initial datum on road 1: $(\rho^0_1, w^0_1)\in [0, {\rho^{\mathrm{max}}}]\times [{w_{L}}, {w_{R}}]$.

Priority rule We focus on RP algorithm. Recall that values of $\beta < 0.5$ give the priority to road 1, while values of $\beta > 0.5$ give the priority to road 2. In Figure 6 we highlight the level curve related to $\beta = 0.5$ using a bold line. In Figure 6a we show the result for the free-flow case (i). The optimal priority $\beta^{opt}$ decreases as $\rho_1^0$ increases. Specifically, if $\rho^0_1 < \rho^0_2 = \rho^0_3 = 15\, {\mathrm{veh/km}}$ then road 2 should have the priority and $\beta^{opt}$ is independent of the speed attitude of drivers $w^{0}_{1}$. On the other hand, if $\rho^0_1 > 15\, {\mathrm{veh/km}}$ then road 1 should have the priority. In this case, $\beta^{opt}$ depends on $w^{0}_{1}$. In fact, it decreases more rapidly for high values of $w^{0}_{1}$. Hence, vehicles with fast drivers should cross the junction in a higher percentage ($1-\beta$) than vehicles with slow drivers. In Figure 6b we show the result for the congested case (ii). As before $\beta^{opt}$ decreases as $\rho_1^0$ increases. We observe that road 2 should always have the priority when slow drivers ($w^0_1 = w_L$) arrive from road 1. On the other hand, road 1 should have the priority for high values of $\rho^0_1$ and $w^0_1$ (region to the right of the curve $\beta^{opt} = 0.5$).

Traffic light Here we analyze how the ratio between the optimal green and red duration $t^{opt}_g/t^{opt}_r$ varies with respect to $(\rho^0_1, w^0_1)$ for the two traffic scenarios (i) and (ii). In Figure 7 we show the result for the free-flow case (i). The left plot represents the level curves of ${\mathcal{F}}$ computed with the optimal couple $(t^{opt}_g, t^{opt}_r)$: the minimum value is increasing in $\rho^0_1$ independently of $w^0_1$; the dependence on $w^0_1$ only occurs when many slow vehicles arrive from road 1 (bottom right of the figure). The right plot shows the level curves of the ratio $t^{opt}_g/t^{opt}_r$, where the bold line identifies the curve with $t^{opt}_g/t^{opt}_r = 1$. We observe that for small values of $\rho^0_1$, the red phase should be longer than the green one. On the other hand, when $\rho^0_1$ increases, the ratio becomes greater than one, and thus vehicles coming from road 1 should have a longer green phase. Again, the solution is not very sensitive to the variations of the speed attitude of drivers $w^0_1$. In Figure 8 we show the result for the congested case (ii). The behavior of ${\mathcal{F}}(t_g^{opt}, t_r^{opt})$ on the left plot is analogous to case (i), while the trend of the ratio $t^{opt}_g/t^{opt}_r$ changes. Indeed, the green phase should be longer than the red one only for high values of $\rho^0_1$ and low values of $w^{0}_{1}$. Finally we observe that in both cases, the minimum of the functional is not very sensitive to small perturbations of optimal $(t^{opt}_g, t^{opt}_r)$. We can summarize the results as follows. For the priority-ruled junction, we obtain the minimum of the functional ${\mathcal{F}}$ by giving the priority to the incoming road with higher density and favoring fast drivers. For the traffic light too, the road with higher density should have a longer green phase. However, when the three roads are congested, vehicles with slow drivers should have a longer green phase. As expected, the sensitivity with respect to $w$ is greater when traffic is congested, that is when it is more influenced by $w$.

7.   Emissions at roundabouts

In this section we study emissions and travel times for a roundabout, modeled combining merge and diverge junctions as depicted in Figure 9. There are four junctions: $J_{1}$ and $J_{3}$ of type $2\to 1$ (merge); $J_{2}$ and $J_{4}$ of type $1\to 2$ (diverge). We focus on the AP algorithm to compute the minimum of problem (4.4), obtaining the priority parameters $\gamma = (\beta_{J_{1}}, \beta_{J_{3}})\in[0, 1]^{2}$. We also compute the optimal timing $\gamma = ((t_{g}, t_{r})_{J_{1}}, (t_{g}, t_{r})_{J_{3}})\in [G\times R]^{2}$ for the roundabout with traffic lights placed at the two merge junctions $J_{1}$ and $J_{3}$, with $G = R = [25\, {\mathrm{s}}, 90\, {\mathrm{s}}]$. We exclude traffic light phases smaller than $25\, {\mathrm{s}}$, and compare the roundabout with priorities with that with traffic lights. Furthermore, we consider the case of periodic traffic lights with $(t_{g}, t_{r})_{J_{1}} = (t_{g}, t_{r})_{J_{3}} = (45\, {\mathrm{s}}, 45\, {\mathrm{s}})$ to compare the optimal traffic lights with naive ones.

The two diverging junctions $J_2$ and $J_4$ have a fixed distribution parameter $\alpha = 0.6$. The model parameters ${\rho_{f}}$, ${\rho^{\mathrm{max}}}$, ${\rho_{c}}$ and ${V^{\mathrm{max}}}$, the length of the roads $L$ and the space step ${\Delta x}$ are fixed as in Table 2. The length of the simulations is $T = 1\, {\mathrm{h}}$ and the time step ${\Delta t} = 2.57\, {\mathrm{s}}$. The initial density is assumed to be null for each road. We analyze three traffic scenarios determined by the density of vehicles which enter into the network from roads 1 and 5. On the latter, we used Dirichlet boundary conditions:

and $w^{n}_{0, r} = ({w_{L}}+{w_{R}})/2$ for $r = 1, 5$. We use Neumann boundary conditions for roads 3 and 7, thus allowing all vehicles to exit the roundabout. The initially empty network is filled up for the first 20 min of simulation, then no more vehicles access the network until the final time $T = 1\, {\mathrm{h}}$. In this way, the emissions are measures both for loading and unloading of the roundabout.

In Table 5 we show the optimal controls and the corresponding functionals values. We observe that ${\mathcal{F}_{E}}$, ${\mathcal{F}_{T}}$ and ${\mathcal{F}}$ grow as the number of vehicles entering the network increases, both for priorities and traffic lights dynamics. In particular, in the case of $\bar\rho = 15\, {\mathrm{veh/km}}$ in Eq (7.1), the traffic lights dynamics produce 20% lower emissions and 2% lower travel times with respect to priorities. In congested situations, instead, the emissions are reduced by about 11% in presence of traffic lights, while the travel times are 6% longer compared to priority-ruled dynamics. Moreover, for all the boundary data $\bar\rho$ the dynamics with periodic traffic lights produce higher emissions and longer travel times compared to the dynamics with optimized traffic lights; in congested situations the increase of ${\mathcal{F}}$ is about 30%. On the other hand, periodic traffic lights produce lower emissions than optimized priorities dynamics, but significantly higher travel times. The higher levels of emissions associated with priorities can be observed also in Figure 10, where we plot the emissions on each road of the network at different times. The emissions associated with traffic lights dynamics show an oscillating behavior which is not observed in the priorities case, see plots in Figures 10a–e. At the final time of the simulation, plots Figures 10c and 10f, the emissions are close to 0 as nearly all vehicles have left the network.

Finally, in Figure 11, we show the change in time of the total emission rates in the whole network. The trend in emission rates is the same for the three cases: emissions rise as vehicles enter the network and then decrease to 0. The peak value grows as $\bar\rho$ increases. Again, periodic traffic lights correspond to higher emissions than optimized ones. In Table 6 we report the total number of vehicles that enter the network for the three tests and the corresponding total amount of emissions produced with the two optimized traffic dynamics. We observe that emissions are more than double when $\bar\rho = 40\, {\mathrm{veh/km}}$ compared to $\bar\rho = 15\, {\mathrm{veh/km}}$ and almost triple when $\bar\rho = 80\, {\mathrm{veh/km}}$ with respect to $\bar\rho = 15\, {\mathrm{veh/km}}$, while the difference between the case of $\bar\rho = 40\, {\mathrm{veh/km}}$ and the one of $\bar\rho = 80\, {\mathrm{veh/km}}$ is smaller.

To check the robustness of our results, we computed the minima of the functional ${\mathcal{F}}$ for different values of the weights $c_{1}$ and $c_{2}$ in Appendix A. The specific values of the functional obviously varies as we change the weights, but not the qualitative and quantitative comparison of priorities with traffic lights. Moreover, the optimal traffic light timing appears to be more robust than the optimal priorities.

8.   Conclusions

In this work, we have extended the Generic Second Order Model to a road network with merge and diverge junctions and proposed a tool to estimate and minimize traffic emissions by regulating traffic dynamics. Such regulation corresponds to the choice of suitable model parameter $\gamma$ that governs the distribution of traffic in a diverge and priorities in a merge.

Different scenarios have been considered, such as: a traffic policeman who strictly enforces the priority rule (RP algorithm), an uncontrolled intersection where drivers tend to maximize the flow (AP algorithm), and the presence of a traffic light. A functional measuring emissions and travel times was tested numerically on a single merge junction, showing that the minimum is achieved by giving the priority and a longer green traffic light to the incoming road with higher density and fast drivers. On the other hand, the test performed on a roundabout has pointed out that traffic lights appear to be convenient with respect to priorities for emissions, especially at low densities. This indicates that the increasingly common roundabouts may benefit from the installation of traffic lights at entrances. We conclude by stating that our approach is very flexible and can easily be used as a decision support for traffic management.

Acknowledgments

C. B. and M. B. acknowledge the Italian Ministry of University and Research (MUR) for supporting this research with funds coming from the PRIN Project 2017 "Innovative numerical methods for evolutionary partial differential equations and applications" (No. 2017KKJP4X). This work was carried out within the research project "SMARTOUR: Intelligent Platform for Tourism" (No. SCN 00166) funded by the MUR with the Regional Development Fund of European Union (PON Research and Competitiveness 2007–2013). C. B. and M. B. are members of the INdAM Research group GNCS. The research of B. P. was partially supported by the NSF CPS Synergy project "Smoothing Traffic via Energy-efficient Autonomous Driving" (STEAD) CNS 1837481. The research of B. P. is based upon work supported by the U.S. Department of Energy's Office of Energy Efficiency and Renewable Energy (EERE) under the Vehicle Technologies Office award number CID DE-EE0008872. The views expressed herein do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

Conflict of interest

The authors declare that there is no conflict of interest.

A.   Senstivity of ${\mathcal{F}}$ to weights $c_{1}$ and $c_{2}$

In this appendix we investigate the sensitivity of the functional ${\mathcal{F}}$ with respect to the weights $c_{1}$ and $c_{2}$ in Eq (5.9) for the roundabout. Our aim is to compare the optimal controls obtained by giving more importance once to emissions and once to the travel time. Therefore, we define ${\mathcal{F}}_{c_{1}} = \kappa c_{1}{\mathcal{F}_{E}}+c_{2}{\mathcal{F}_{T}}$ and ${\mathcal{F}}_{c_{2}} = c_{1}{\mathcal{F}_{E}}+\kappa c_{2}{\mathcal{F}_{T}}$ with $\kappa = 10, 100$.

In Tables A1 and A2 we report the optimal controls computed for ${\mathcal{F}}_{c_{1}}$ and ${\mathcal{F}}_{c_{2}}$, using the Dirichlet boundary conditions in Eq (7.1) for different $\bar\rho$ as in Section 7. First, we observe that the values of the functional ${\mathcal{F}}_{c_{1}}$ are lower than those of the functional ${\mathcal{F}}_{c_{2}}$. Therefore, giving more importance to emissions rather than to travel time allows to reduce the total cost. Analogously to the case of functional ${\mathcal{F}}$ studied in Section 7, in all cases traffic lights dynamics are convenient in terms of emissions production, while the travel time is shorter when traffic is ruled by priorities. Finally, note that the optimal priorities are influenced by the choice of the functional, while the optimal traffic light timing is always the same for all the tests.