Liudmila Tumash

Carlos Canudas-De-Wit

Maria Laura

Delle Monache

Boundary Control for Multi-Directional Traffic on Urban Networks

This paper is devoted to boundary control design for urban traffic described on a macroscopic scale. The state corresponds to vehicle density that evolves on a continuum two-dimensional domain that represents a continuous approximation of a urban network. Its parameters are interpolated as a function of distance to physical roads. The dynamics are governed by a new macroscopic multi-directional traffic model that encompasses a system of four coupled partial differential equations (PDE) each describing density evolution in one direction layer: North, East, West and South (NEWS). We analyse the class of desired states that the density governed by NEWS model can achieve. Then a boundary control is designed to drive congested traffic to an equilibrium with the minimal congestion level. The result is validated numerically using the real structure of Grenoble downtown (a city in France).

I. Introduction

Ever growing urban areas face the problem of transportation efficiency drop during rush hours. This triggers challenges for researchers to develop realistic traffic models to predict congestion formation as well as to suggest efficient control measures to mitigate it. The first macroscopic traffic model was introduced by Lighthill and Whitgham [START_REF] Lighthill | On kinematic waves, II: A theory of traffic flow on long crowded roads[END_REF] and Richards [START_REF] Richards | Shock waves on the highway[END_REF] who elaborated the kinematic wave theory for traffic. The LWR model is a first-order hyperbolic PDE based on conservation of vehicles that describes traffic evolution on a single road. Its key assumption is the existence of a concave relation between traffic flow and density known as fundamental diagram. LWR framework was extended to networks by adding Riemann problems at intersections in [START_REF] Lebaque | First-order macroscopic traffic flow models: Intersection modeling, network modeling[END_REF]. However, computational cost may become too high if one considers a large network, which requires development of macroscopic approaches for urban traffic modelling.

Alternatively, one can use 2D continuum models to describe urban traffic. These share a lot of features with pedestrian models [START_REF] Hughes | A continuum theory for the flow of pedestrians[END_REF] the main difference being that vehicles are restricted to move on roads. A model including diffusion and a drift term depending on density and network geometry was considered in [START_REF] Della Rossa | A distributed model of traffic flows on extended regions[END_REF]. Another work [START_REF] Jiang | A dynamic traffic assignment model for a continuum transportation system[END_REF] defines the direction of traffic motion by solving Eikonal equations. Recently, [START_REF] Mollier | Twodimensional macroscopic model for large scale traffic networks[END_REF] extended the LWR model to two dimensions with a space-dependent FD that incorporates network infrastructure parameters. A general method to calculate steady states in 2D LWR was presented in [START_REF] Tumash | Equilibrium Manifolds in 2D Fluid Traffic Models[END_REF], while [START_REF] Tumash | Topology-based control design for congested areas in urban networks[END_REF] elaborated a boundary controller for congested traffic that was further extended to traffic being in any (mixed) regime in [START_REF] Tumash | Boundary and VSL Control for Large-Scale Urban Traffic Networks[END_REF]. L. Tumash and C. Canudas-de-Wit are with Univ. Grenoble Alpes, CNRS, Inria, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. M. L. Delle Monache is with Department of Civil and Environmental Engineering, UC Berkeley, USA.

The aforementioned references consider only unidirectional traffic. The first attempt to include multiple directions was made recently by [START_REF] Lin | A predictive continuum dynamic user-optimal model for a polycentric urban city[END_REF] who deployed dynamic user-optimal principle for the path choice. Its main drawback is that traffic density may become unbounded (it is not based on a fundamental diagram). There are also other works [START_REF] Mollier | A step towards a multidirectional 2D model for large scale traffic networks[END_REF], [START_REF] Aghamohammadi | A Continuum Model for Cities Based on the Macroscopic Fundamental Diagram: a Semi-Lagrangian Solution Method[END_REF] proposing 2D multi-layer models with bounded densities. However, these models do not include mixing between direction layers (vehicles can not change the direction). Then, they are also not necessarily hyperbolic, i.e., their type varies with parameters, which exaggerates its analysis and numerical simulation. Our recent work [START_REF] Tumash | Multi-Directional Continuous Traffic Model For Large-Scale Urban Networks[END_REF] fixes both of these aspects introducing a multi-directional model that uses network geometry to assign direction layers to corresponding density evolutions: North, East, West and South (NEWS model). The direction layers are coupled, i.e., the model captures vehicles that change the direction of movement. Information about traffic flow direction is obtained from turning ratios at intersections.

In this paper we design a boundary control for the NEWS model that drives an initially congested state to the best possible desired equilibrium corresponding to congestion minimization, which equivalently means throughput maximization of the network. Our main contribution is an extensive analysis of possible space-varying profiles that the system can achieve, which is far from being trivial for multidirectional traffic systems. We then use Lyapunov methods to show exponential convergence to the desired state.

II. NEWS Model

We use the NEWS model introduced in [START_REF] Tumash | Multi-Directional Continuous Traffic Model For Large-Scale Urban Networks[END_REF] to predict traffic evolution on a general urban network. The state is described in terms of vehicle density ρ = (ρ N , ρ S , ρ W , ρ E ) T ∈ R 4 evolving in 4 cardinal direction layers: North (N), East (E), West (W) and South (S). It evolves on a bounded 2D continuum plane Ω ⊂ R 2 that is a rectangular domain bounded by x min , x max , y min and y max . It approximates a urban network whose road parameters are interpolated in the domain as a function of distance to these roads. The NEWS model describes evolution of ρ(x, y, t) : Ω × R + → R + by the following system of PDEs:

                                     ∂ρ N ∂t = 1 L ψ in N -ψ out N - ∂(c N ψ N ) ∂x - ∂(s N ψ N ) ∂y , ∂ρ S ∂t = 1 L ψ in S -ψ out S - ∂(c S ψ S ) ∂x - ∂(s S ψ S ) ∂y , ∂ρ W ∂t = 1 L ψ in W -ψ out W - ∂(c W ψ W ) ∂x - ∂(s W ψ W ) ∂y , ∂ρ E ∂t = 1 L ψ in E -ψ out E - ∂(c E ψ E ) ∂x - ∂(s E ψ E ) ∂y . (1)
Flux function ψ(x, y, ρ) is a four-dimensional vector ψ = (ψ N , ψ S , ψ W , ψ E ) T . Let us consider the North direction to simplify definitions. Then, ψ N (x, y, ρ N ) : Ω × R + → R + is a concave function that achieves maximum ψ max N (x, y) ∀(x, y) ∈ Ω (capacity) at the critical density ρ c N (x, y), and minimum is achieved at ψ N (x, y, 0) = ψ N (x, y, ρ max N ) = 0 with ρ max N (x, y) being the space-dependent traffic jam density. This concave relation is known as fundamental diagram (FD). Vehicles move freely to the North with a positive kinematic wave speed v N (x, y) if density in this direction is below the critical value, i.e., free-flow regime occurs for ρ N (x, y) ∈ [0, ρ c N (x, y)) ∀(x, y) ∈ Ω. Otherwise vehicles move in congested regime with a negative kinematic wave speed

ω N (x, y) if ρ N (x, y) ∈ [ρ c N (x, y), ρ max N (x, y)] ∀(x, y) ∈ Ω.
Traffic flow in other directions ψ S , ψ W , and ψ E can be retrieved from FD of corresponding directions. Here we use the triangular FD [START_REF] Daganzo | The cell transmission model: A dynamic representation of highway traffic consistent with the hydrodynamic theory[END_REF], and its parameters are defined ∀(x, y) ∈ Ω as:

ρ c N = ω N v N + ω N ρ max N , ψ max N = v N ω N v N + ω N ρ max N ,
and we also assume that ρ c = ρ max /3 for all directions.

In (1), ψ in (x, y, ρ) -ψ out (x, y, ρ) is the net flow at roads:

              ψ in N -ψ out N ψ in S -ψ out S ψ in W -ψ out W ψ in E -ψ out E               =               ψ S N + ψ WN + ψ EN -ψ NS -ψ NW -ψ NE ψ NS + ψ WS + ψ ES -ψ S N -ψ S W -ψ S E ψ NW + ψ S W + ψ EW -ψ WN -ψ WS -ψ WE ψ NE + ψ S E + ψ WE -ψ EN -ψ ES -ψ EW               .
Thereby ψ S N , ψ WN , ψ EN and so on are partial flows. As an example, ψ S N (x, y, ρ S , ρ N ) is a flow of cars that were going to South and then turned to North. It is defined as a function of demand from the South D S and supply of the North S N :

ψ S N = min {α S N D S (ρ S ), β S N S N (ρ N )} , (2)
where α S N and β S N are the turning ratio from South to North and the supply ratio of the North for vehicles arriving from the South, correspondingly. The demand and supply functions in (2) are defined for triangular FD as follows:

D S (ρ S ) =      v S ρ S , if 0 ≤ ρ S ≤ ρ c S , ψ max S , if ρ c S < ρ S ≤ ρ max S , (3)
S N (ρ N ) =      ψ max N , if 0 ≤ ρ N ≤ ρ c N , ω N (ρ max N -ρ N ), if ρ c N < ρ N ≤ ρ max N , (4)
Notice that demand can not exceed the maximal flow, since it is the highest possible flow determined by the number of lanes and speed limit (transportation capacity). We can now insert (4), (3) into (2) and obtain:

ψ S N = min{α S N v S ρ S , α S N ψ max S , β S N ω N ρ max N -ρ N , β S N ψ max N }. (5)
In [START_REF] Tumash | Multi-Directional Continuous Traffic Model For Large-Scale Urban Networks[END_REF] the NEWS model was first derived to describe density propagation in a vicinity of any intersection in a unified way. It was then extended to networks by defining intersection parameters everywhere in domain Ω through interpolation (details are given below). Thus, if we consider (1) for a fixed point (x, y) ∈ Ω being a location of some intersection with n out outgoing roads each having length l j for j ∈ {1, . . . , n out }, then L(x, y) in (1) should be interpreted as the mean length of outgoing roads of this intersection, and it was chosen such that the number of cars is conserved:

L =          n out j=1 ρ max j l j          /          n out j=1 ρ max j          .
Finally, c N , c S , c W , c E and s N , s S , s W , s E from (1) describe a general orientation of roads, e.g., consider North:

c N = n out j=1 p N θ j cos θ j ψ max j n out j=1 p N θ j ψ max j , s N = n out j=1 p N θ j sin θ j ψ max j n out j=1 p N θ j ψ max j ,
where p N θ j ∈ [0, 1] is a coefficient indicating orientation of road j w.r.t. North (1 means that road j is pointing exactly to the North, 0 means that angle θ N j ≥ π/2 w.r.t. North). Any urban network is given as a set of roads and intersections with given speed limits and number of lanes. The NEWS model ( 1) is a PDE that requires all intersection and FD parameters to be defined ∀(x, y) ∈ Ω (that is, also in areas where there are no roads). First, we define the maximal density ρ max by placing vehicles at every road with a headway distance of 6 m. We assume that every vehicle contributes to global density with a Gaussian kernel with standard deviation 50 m centred at its position (see [START_REF] Mollier | Twodimensional macroscopic model for large scale traffic networks[END_REF]). Further, we project all intersection and FD parameters α, β, L, cos θ, sin θ, ρ max , v and ω into NEWS formulation using network geometry data (see [START_REF] Tumash | Multi-Directional Continuous Traffic Model For Large-Scale Urban Networks[END_REF] for a detailed description). Then, these variables are interpolated using Inverse Distance Weighting that assignes values to parameters as a function of distance to roads, see also [START_REF] Mollier | Twodimensional macroscopic model for large scale traffic networks[END_REF].

III. Desired Steady-State

A. Congested Traffic

We will now linearize [START_REF] Lighthill | On kinematic waves, II: A theory of traffic flow on long crowded roads[END_REF], which is possible if traffic is considered to be either in pure free-flow or pure congested regime. Our goal is to mitigate congestions, and hence it is feasible to consider it for congested regime. Thus, we assume a heavily congested initial state that should be driven to the state of minimal congestion (remaining in the same regime). For the congested regime, minimum in ( 2) is always resolved to the benefit of supply, which in turn implies that

ψ S N = β S N ω N ρ max N -ρ N , ∀(x, y) ∈ Ω.
Using this simplification for model (1), we fix ρ 0 (x, y) ∀(x, y) ∈ Ω as an initial condition (a function of bounded variation), and introduce the following Initial Boundary Value Problem that describes evolution of congested traffic on domain Ω with Γ ⊂ Ω being its boundary: where Γ out ⊂ Γ is a set of boundary points (x, y) associated with the domain exit:

                           ∂ρ ∂t = 1 L (I -B)W ρ max -ρ - ∂[C W (ρ max -ρ)] ∂x - ∂[S W (ρ max -ρ)] ∂y , ρ(x, y, t) = u(x, y, t), ∀(x, y) ∈ Γ out ρ(x, y, 0) = ρ 0 (x, y), (6)
Γ out = (y max , x max , x min , y min ) T .
Congested traffic is controlled at downstream boundary Γ out by specifying the control vector u = (u N , u E , u W , u S ) T . See Fig. 1, where arrows denote boundaries to be activated for control in each direction. Boundary control can be physically realized by setting e.g. traffic lights at roads' exits. Finally, C, S , W and B in ( 6) are all 4 × 4 matrices such that C and S are diagonal matrices, W is a positive-definite diagonal matrix, and B is a nonnegative matrix:

C = diag(c), S = diag(s), W = diag(ω), B =               β NN β NE β NW β NS β EN β EE β EW β ES β WN β WE β WW β WS β S N β S E β S W β S S               .

B. Desired Steady-State

Define error from a desired space-varying equilibrium as:

ρ(x, y, t) = ρ(x, y, t) -ρ d (x, y), ∀(x, y) ∈ Ω.
To simplify the mathematical analysis, we restrict our study to desired profiles ρ d having values only in the congested regime, i.e., ρ d (x, y) ≥ ρ c (x, y) ∀(x, y) ∈ Ω.

The time derivative of error coincides with that of state [START_REF] Lebaque | First-order macroscopic traffic flow models: Intersection modeling, network modeling[END_REF]. The desired density is constant in time and thus:

1 L (I -B)W ρ max -ρ d = ∂[C W ρ max -ρ d ] ∂x + ∂[S W ρ max -ρ d ] ∂y , (7)
which lets us write the equation for error dynamics as

∂ ρ ∂t = 1 L (B -I)W ρ + ∂[C W ρ] ∂x + ∂[S W ρ] ∂y . (8)
We seek to find a desired density distribution that corresponds to congestion minimization, and then to design a boundary control that achieves that desired equilibrium. The desired density must remain in the congested regime, and its boundary values should be proportional to maximal densities at corresponding coordinates, i.e., ∃γ ∈ [1/3, 1] such that ρ (x, y) = γρ max (x, y), ∀(x, y) ∈ Γ out .

The range of constant γ coincides with the range of densities in congested regime (recall ρ c = 1/3ρ max ).

Problem 1. Find the desired space-varying density ρ d (x, y) ∀(x, y) ∈ Ω that corresponds to state of minimal congestion under constraints: ρ d (x, y) ≥ ρ c (x, y) ∀(x, y) ∈ Ω, and boundary values being proportional to maximal densities (9).

Remark 1. Minimizing congestion means finding

ρ d (x, y) ≥ ρ c (x, y) ∀(x, y) ∈ Ω such that the L ∞ norm ρ d (•) -ρ c (•) ∞ is minimized, i.
e., we want density to be as small as possible.

Remark 2. Physically, ( 9) implies that boundaries are filled in a homogeneous way. For example, imagine vehicles concentrated in a city center that tend to leave it simultaneously, e.g., when people drive back home from their offices.

In order to find a desired profile satisfying Problem 1, we need to solve [START_REF] Della Rossa | A distributed model of traffic flows on extended regions[END_REF]. First, we need to introduce a change of variables ρ(x, y) ∀(x, y) ∈ Ω as

ρ(x, y) = ρ max (x, y) -ρ d (x, y), (10)
which being inserted into (7), yields

1 L (I -B)W ρ = ∂[C W ρ] ∂x + ∂[S W ρ] ∂y . (11)
Then, the following steps are performed: 1) Initial guess: set the desired density values at the downstream boundaries Γ out equal to corresponding critical values, i.e., pick the lowest γ = 1/3, which leads to ρ(x, y) = 2 3 ρ max (x, y), ∀(x, y) ∈ Γ out .

2) Define a numerical grid: divide Ω into equal cells (n x = n y ). The discretization steps are ∆x = (x maxx min ) /n x , ∆y = (y maxy min ) /n y .

Then, a grid point is given by (i∆x, j∆y) for i ∈ {1, . . . , n x }, j ∈ 1, . . . , n y .

3) Discretize PDE (11): now consider North, and the same steps should be done for other directions. In accordance with the upwind scheme, discretization of c N ω N ρN and s N ω N ρN depends on signs of c N and s N . Hence, we can define 4 × 4 diagonal matrices Q x , Q y , R x and R y that capture the upwind scheme: c N,i, j > 0 :

Q x NN,i, j = c N,i+1, j ω N,i+1, j ∆x , R x NN,i, j = 0, else : Q x NN,i, j = 0, R x NN,i, j = - c N,i-1, j ω N,i-1, j ∆x ,
and the same can be written for s N and y-direction, for which we fix i and vary j. Define also a 4 × 4 matrix P as:

P i, j = 1 L i, j (B i, j -I)W i, j - |C i, j |W i, j ∆x - |S i, j |W i, j ∆y .
Using the definition of matrices P, Q x , Q y , R x and R y , we can now write the PDE system for ρ given by ( 11) in a discretized form that reads ∀(i, j) ∈ {1, . . . , n x } × {1, . . . , n y }:

P i, j ρi, j + Q x i, j ρi+1, j + Q y i, j ρi, j+1 + R x i, j ρi-1, j + R y i, j ρi, j-1 = 0. (12)
4) Find solution to system (12): using dimensional splitting, i.e., at each iteration first x and then y steps are performed. At each x step, terms ρi, j-1 and ρi, j+1 take fixed values from the previous iteration, while ρi-1, j and ρi+1, j are fixed for each y step. At x step the system ( 12) is solved for every j by block tridiagonal matrix algorithm, while at y step this algorithm is applied for every column i.

5) Find optimal solution: a numerical solution ρ for ( 12) is not necessarily optimal. Since ( 12) is a linear system, α ρ for α ∈ [0, 1] is also its solution. Let us estimate parameter α * providing the optimal equilibrium as in Problem 1. Consider the desired state ρ d obtained from [START_REF] Mollier | Twodimensional macroscopic model for large scale traffic networks[END_REF]:

ρ d = ρ max -α ρ. ( 13
)
By choosing α = 0 we obtain ρ d = ρ max , while by choosing α = 1 we achieve ρ d = ρ c at the boundaries (see step 1 and use ρ c = ρ max /3). This implies that by taking an intermediate value of α we guarantee congested regime at the boundaries. Let us calculate α * that provides for ρ d to be as close as possible to ρ c while preserving congested regime in Ω (see Remark 1), for which in general we can write:

ρ d ρ c ≥ 1 ⇒ ρ max -α ρ 1/3ρ max ≥ 1 ⇒ α ≤ 2 3 ρ max ρ .
From the discussion above it follows that the optimal state is achieved if ∃(x * , y * ), for which

α * = min (x,y)∈Ω r∈{N,S ,W,E} 2 3 ρ max r (x, y) ρr (x, y) . ( 14
)
The optimal desired profile in the whole domain Ω can be obtained from [START_REF] Tumash | Equilibrium Manifolds in 2D Fluid Traffic Models[END_REF]

for optimal α = α * ρ d (x, y) = ρ max (x, y) -α * ρ(x, y), (x, y) ∈ Ω, (15)
with α * given by [START_REF] Tumash | Topology-based control design for congested areas in urban networks[END_REF]. To get an expression for ρ d (x, y) at the boundary ∀(x, y) ∈ Γ out , we take ρ from step 1 and insert it into [START_REF] Tumash | Boundary and VSL Control for Large-Scale Urban Traffic Networks[END_REF], which yields

ρ d (x, y) = γ * ρ max (x, y), with γ * = 1 - 2 3 α * ∀(x, y) ∈ Γ out . (16)
This expression directly determines boundary control variables u(x, y) from ( 6), see details below.

IV. Boundary Control Design

After we have analyzed the desired profile corresponding to the state of minimal congestion (see Remark 1), let us formulate the boundary control design problem as follows.

Problem 2. Find a boundary controller u(x, y) that drives congested traffic governed by [START_REF] Lebaque | First-order macroscopic traffic flow models: Intersection modeling, network modeling[END_REF] to the desired space-varying state ρ d (x, y) given by ( 15) ∀(x, y) ∈ Ω as t → ∞.

To prove convergence to the desired profile, we have to restrict to a special network structure, e.g. a Manhattan grid. Assumption 1. Matrices C and S from ( 6) are constant in space, e.g., they can be defined as:

c N = 0, c S = 0, c W = -1, c E = 1, s N = 1, s S = -1, s W = 0, s E = 0. (17)
In general, further analysis requires these variables to be just constant in space, but we choose (17) for simplicity. We also make an assumption on supply ratios: Assumption 2. Matrix B is constant in space, which implies that every intersection has the same turning ratio pattern.

Theorem 1. Under Assumptions 1 and 2, let the boundary controller be defined ∀(x, y) ∈ Γ out as

u(x, y) = ρ d N (x, y max ), ρ d S (x, y min ), ρ d W (x min , y), ρ d E (x max , y) T , (18)
then ∃K, k > 0 such that

ρ(t) -ρ d 2 L 2 e -kt K ρ(0) -ρ d 2 L 2 ,
i.e., the state ρ(x, y, t) exponentially converges to the desired equilibrium ρ d (x, y) ∀(x, y) ∈ Ω as t → ∞.

Remark 3. Although for simplicity we made Assumption 1, control (18) can be applied to a more general network, as we will show on a numerical example.

Proof of Theorem 1. Let us first analyse matrix B -I. Its non-diagonal elements are positive and its diagonal elements are negative. Moreover, it has one zero eigenvalue and all others are negative, see Appendix I. Hence, B-I is a negative singular M-matrix with one zero eigenvalue. Thus, there exists a positive-definite diagonal 4 × 4 matrix D such that

D(B -I) + (B T -I)D ≤ 0. ( 19
)
Let us also introduce a diagonal 4 × 4 matrix E composed by exponential functions. We then define the following Lyapunov function candidate:

V =
(W ρ) T DE ∂[C W ρ] ∂x + ∂[S W ρ] ∂y dydx. (21)
Let us now denote the first term of (21) as V1 and the second term as V2 . Term V1 is negative due to (19) and the fact that matrix E is non-negative. We further consider V2 by inserting the values of matrices C and S (17) from Assumption 1 V2 = 2

x max

x min y max y min

ω E ρE D E e x ∂(ω E ρE ) ∂x -ω W ρW D W e -x ∂(ω W ρW ) ∂x + ω N ρN D N e y ∂(ω N ρN ) ∂y -ω S ρS D S e -y ∂(ω S ρS ) ∂y dydx.
This expression is then integrated by parts, which yields:

By setting the boundary controller as in (18) ∀t ∈ R + , we achieve ρN (x, y max , t) = 0, ρS (x, y min , t) = 0, ∀x ∈ [x min , x max ], ρW (x min , y, t) = 0, ρE (x max , y, t) = 0, ∀y ∈ [y min , y max ],

and one ensures that the first four integrals in (22) go to zero. The last term in (22) can be bounded as follows where we have used the definition of the Lyapunov function (20) and that ω is positive. This means that by inserting (23) into (22) and by using the bound from (24), we can write:

V = V1 + V2 ≤ V2 ≤ -kV, where k ∈ R + is a positive constant k = min (x,y)∈Ω r∈{N,E,W,S } ω r (x, y).
One can also prove that error ρ converges to zero in L 2 norm exponentially. Indeed, note that the Lyapunov function V from (20) defines an equivalent norm on density space:

m ρ 2 L 2 V M ρ 2 L 2 , with m = min (x,y)∈Ω r∈{N,S ,W,E} ω r (x, y)D r E r (x, y), M = max (x,y)∈Ω r∈{N,S ,W,E} ω r (x, y)D r E r (x, y).
By exponential convergence of Lyapunov functions we have

V(t) e -kt V(0) ⇒ ρ(t) 2 L 2 e -kt M m ρ(0) 2 L 2 .
Remark 4. Assumption 2 can be relaxed, if one can find matrix D that satisfies inequality (19) and whose elements D E (y) and D W (y) may depend on y, while D N (x) and D S (x) may depend on x.

V. Numerical Example

We consider Grenoble downtown with a total surface of 1.4 × 1 km 2 as a network. We define a numerical grid for n x = 60 and n y = 60, and deploy the 2D Godunov scheme to simulate density governed by ( 6) with downstream boundary conditions set to the desired optimal density as in (18), while the upstream boundary conditions are initialised with φ max . We will demonstrate how the boundary controller (18) mitigates congestion given the initial state

ρ 0 (x, y) = ρ max (x, y), ∀(x, y) ∈ Ω.
The results of control performance for the network of Grenoble downtown are shown in Fig. 2. Fig. 2a) illustrates the initial vehicle density (traffic jam). The optimal desired equilibrium ρ d is illustrated in Fig. 2b). It was found by first solving PDE for ρ (12) and then using [START_REF] Tumash | Boundary and VSL Control for Large-Scale Urban Traffic Networks[END_REF], where we use α * = 0.51 obtained using [START_REF] Tumash | Topology-based control design for congested areas in urban networks[END_REF]. Further, we show the impact of boundary controller (18) on congested traffic after t = 5 min, t = 20 min and t = 50 min in Fig. 2c), 2d) and 2e), respectively. The controlled state at t = 50 min looks identical to the desired equilibrium.

We deploy the Structural Similarity Measure (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] to enable a quantitative comparison between the controlled state on Fig. 2c) -e) and the desired state on Fig. 2b). This index is a perception-based metric used to detect structural changes. The range of SSIM is [-1, 1], where S S I M = 1 means that two images are identical and S S I M = -1 indicates that the second image is inverse of the first. To refine the computations, we also divide the network of Grenoble into 9 equal zones. SSIM is calculated for every zone, and then its mean value S S I M is found as an arithmetic average. The result is shown on Fig. 2f), where S S I M converges to 1 indicating that the desired profile is achieved.

VI. Conclusions

We investigated the multi-directional NEWS model from the control perspective. In particular, we analysed the class of desired equilibria that must satisfy a certain system of PDEs. We have posed and solved the problem of finding an equilibrium state that provides congestion minimization in a network, under the constraint that its range must remain in the congested regime. The desired state was assumed to be proportional to the maximal densities at the boundaries. Further, we proved the exponential convergence of a congested state to this desired equilibrium using Lyapunov methods. Finally, we used the real network of Grenoble downtown to produce a numerical example that stays in a good agreement with the theoretical result. Thereby, at each time step we have calculated the structural similarity index to show the convergence of our density distribution to the desired one.

We see two major promising directions for future studies: finding equilibria that admit mixed traffic regimes, and elaborating boundary control that is set to points on roads rather than being a continuous line.

Fig. 1 :

 1 Fig. 1: Vehicle density in Grenoble downtown. Downstream boundaries for control are indicated by colorful arrows: North in blue (u N ), East in dark red (u E ), West in green (u W ) and South in orange (u S ).

  ω N D N e y + ρ2 E ω E D E e x + ρ2 W ω W D W e -x + ρ2 S ω S D S e -y ) dydx,(20)where D N , D E , D W and D S are diagonal elements of matrix D. Function (20) is obviously positive-definite, since matrix WDE > 0. Let us now take its time derivative, and then insert the error dynamics ∂ ρ/∂t from (8), which yields: ) T DE(B -I) + (B T -I)DE W ρ dydx

2 +

 2 ( D W ω W ρW ) 2e x ( D E ω E ρE ) 2 D E ω E ρE ) 2e -x ( D W ω W ρW ) 2 ( D S ω S ρS ) 2e y ( D N ω N ρN ) 2 D N ω N ρN ) 2e -y ( D S ω S ρS ) 2 E (ω E ρE ) 2 + e -x D W (ω W ρW ) e y D N (ω N ρN ) 2 + e -y D S (ω S ρS ) 2 dydx.

  E (ω E ρE ) 2 + e -x D W (ω W ρW ) 2 + e y D N (ω N ρN ) 2 + e -y D S (ω S ρS ) 2 dydx ≤min (x,y)∈Ω r∈{N,S ,W,E} ω r (x, y)V (24)

Fig. 2 :

 2 Fig. 2: Boundary control in Grenoble downtown: a) initial congested state ρ 0 , b) desired equilibrium ρ d ; controlled state after: c) t = 5 min, d) t = 20 min, e) t = 50 min; d) S S I M between state and desired density as a function of time.

ACKNOWLEDGEMENT Scale-FreeBack project received funding from the European Research Council (ERC) under European Union's Horizon 2020 research and innovation program (grant agreement N 694209).

Appendix I Eigenvalues of matrix β -I

Let us now analyse eigenvalues of matrix B-I from [START_REF] Lebaque | First-order macroscopic traffic flow models: Intersection modeling, network modeling[END_REF]. To simplify the notations, we introduce B = B-I. By Gershgorin circle theorem, every eigenvalue of B lies within at least one of the Gershgorin discs d( bii , R i ), where d is a closed disc centered at bii with radius

Consider the first row of matrix B. The Gershgorin disc is centred at β NN -1 and its radius is

The remaining rows of matrix B can be analysed in the same way. Due to Gershgorin theorem, in general, every result looks similar to:

which implies that Reλ( B) ≤ 0 ∀λ( B) and if Reλ( B) = 0, then λ( B) = 0.

Let us consider λ( B) = 0 with x being the corresponding eigenvector:

x T B = 0 = x T λ( B).

Using the definition of matrix B, we further get

Thus, it follows that x is also the eigenvector of matrix B associated with the eigenvalue λ(B) = 1. Note that matrix B is a positive matrix, i.e., β i j > 0 for 1 ≤ i, j ≤ 4 (assume we have no zero turning ratios). Then by Perron-Frobenius theorem λ(B) = 1 is a Perron root (since all columns of B sum to 1), and thus it is a simple root. It follows that all the eigenvalues of matrix B = (B -I) are strictly negative and only one eigenvalue is zero.