Carlotta Donadello

Nicolas Forcadel

Philippe Helluy

Laurent Navoret

Yannick Co-Encadrant

Ulrich Privat

Razafison

Emmanuel Franck

Thibault Liard

Optimal control of nonlinear hyperbolic systems on networks: gradient-based and deep Learning approaches

CONTENTS

CONTENTS

Strasbourg University

Introduction en français

Les équations dérivées de la mécanique des fluides peuvent être utilisées pour modéliser une large variété de phénomènes physiques intervenant dans des situations très différentes. Au-delà même des études de comportements de liquides ou de gaz qui viennent immédiatement à l'esprit avec, par exemple, des applications en météorologie, climatologie ou en aérodynamisme, il est pertinent de noter qu'une très large classe de systèmes présente un comportement similaire lorsqu'ils sont étudiés à l'échelle macroscopique. L'un de leurs principaux points communs est certainement que, contrairement à ce que l'on peut observer par exemple avec des modèles paraboliques en mécanique du solide, ils ne propagent l'information qu'à une vitesse finie. Les exemples incluent bien sûr la propagation d'ondes électromagnétiques modélisées par les équations de Maxwell, mais également le cheminement de requêtes http à travers un réseau internet de routeurs interconnectés, ou l'évolution d'une épidémie entre différentes villes reliées par des routes. L'objectif principal de cette thèse est la simulation numérique et le contrôle optimal de ce type de système, en utilisant dans un premier temps diverses méthodes déterministes basées sur le gradient, puis l'apprentissage profond avec des réseaux de neurones pour réduire efficacement la taille du problème.

La première partie de ce manuscrit s'intéresse à un modèle hydrodynamique de trafic routier qui est une extension par Coclite et al. [START_REF] Coclite | Traffic flow on a road network[END_REF] du modèle de Lighthill-Whitham-Richards [START_REF] Lighthill | On kinematic waves -ii -a theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF] à un graphe représentant un réseau routier. La gestion des intersections d'un tel réseau constitue un défi autant du point de vue de la modélisation que du contrôle, à cause notamment du couplage fortement non linéaire qui relie des solutions très peu régulières (L 1 ). Par ailleurs, pour des raisons de persistance des données ainsi que de reproductibilité des résultats, le code associé aux simulations numérique est libre. Sa documentation fait l'objet du chapitre 3.

Dans un deuxième temps, nous nous intéresserons à un moyen d'accélérer les calculs, largement utilisé dans la recherche et l'industrie, les bases réduites [START_REF] Ito | A reduced-order method for simulation and control of fluid flows[END_REF]. Il s'agit d'une question cruciale, car les algorithmes d'optimisation et de contrôle peuvent avoir besoin d'usage intensif de simulations non linéaires, et il est impératif de rendre ces appels rapides tout en conservant une bonne fiabilité des résultats. Nous verrons que l'approche classique de réduction de modèle se heurte à la forte non-linéarité du système considéré, et proposerons une approche de réduction robuste utilisant un réseau de neurones pour évaluer le terme fortement non linéaire.

Le manuscrit se conclut par un projet réalisé en colaboration avec Leo Meyers et Florent Noisette sous la direction de Bastien Polizzi, Sébastian Minjeaud, Olivier Bernard et Thierry Goudon lors d'une école d'été. Ces travaux traitent eux-aussi de non-linéarité dans un modèle de mécanique des fluides, cette fois-ci appliqué à la biologie. Il s'agit en effet de simuler de façon précise l'évolution de micro-algues et de leur matrice extracellulaire dans un milieu aqueux. L'apport principal de ce chapitre est le développment d'un prédicteur-correcteur sur grille décalée permettant de restituer correctement le champ des vitesses malgré des couplages fortement non linéaires entre les phases ainsi que plusieurs forces et phénomènes biologiques comme la gravitation ou la

∂ t ρ(t, x) + ∂ x f (ρ(t, x)) = 0, (t, x) ∈ (0, T ) × [a, b]. (1)
à laquelle on adjoint des conditions aux limites appropriées faisant intervenir des hypothèses de modélisation quant à la répartition des flux des véhicules aux jonctions à chaque instant.

Un point important est la régularité et la compacité des solutions à variation bornée, cette notion, illustrée défavorable dont nous chercherons à nous prémunir. Nous nous intéressons à cette question que nous reformulons comme un problème de contrôle optimal. Nous considérons un modèle macroscopique de trafic routier sur réseau, semi-discrétisé en espace et décidons de nous donner la possibilité de contrôler le flux à certaines jonctions.

Notre objectif est de lisser et diminuer le trafic le long d'un chemin donné dans un temps donné. Une contrainte de parcimonie est imposée sur les contrôles, afin de s'assurer que les stratégies optimales sont réalisables en pratique.

On est donc ammenés, après discrétisation de [START_REF] Ahmed | On closures for reduced order models $-$ A spectrum of first-principle to machine-learned avenues[END_REF] avec N c ∈ N mailles sur chacune des N r ∈ N routes :

ρ(t) = ρ i,j (t) 1≤i≤Nr 1≤j≤Nc =   1 △x i,j x i,j+ 1 2 x i,j-1 2 ρ(t, x) dx   1≤i≤Nr 1≤j≤Nc , (2)
CONTENTS CONTENTS Strasbourg University à minimiser une fonctionnelle du type :

J (u) = i∈roads j∈mesh(i) ρ i,j (T ; u), (3)
sous la contrainte dynamique donnée par le schéma aux volumes finis contrôlé suivant :

         dρ dt t = f F V (ρ(t), γ(t)) , t ∈ (0, T ) γ(t) = ϕ LP (ρ(t), u(t)) , t ∈ (0, T ) ρ(0) = ρ 0 , (4)
γ 1 (t) = γ 2 (t) = min γ demand 1 (t), (1 -u 2 (t))γ supply 2 (t) , (5)
où γ(t) = (γ 1 (t), γ 2 (t)) et γ demand

1 (t), γ supply 2
(t) sont respectivement les flux maximaux en entrée et en sortie de la jonction. Ils sont calculés à partir des densitées à l'instant t.

Remark 1. Il est important de noter ici que la façon de définir la dépendance de ϕ LP par rapport au contrôle u est un choix de modélisation crucial auquel il faut faire attention pour conserver le caractère bien posé du problème démontré dans [START_REF] Coclite | Traffic flow on a road network[END_REF], en particulier quant à la conservation de la densité de voitures à travers les jonctions.

Nous effectuons ensuite une analyse du problème de contrôle optimal, prouvant l'existence d'un contrôle optimal et dérivant des conditions d'optimalité, que nous réécrivons sous la forme d'une unique équation fonctionnelle. Nous utilisons ensuite cette formulation pour dériver un algorithme mixte que nous interprétons comme un mélange de deux méthodes : une méthode de descente combinée avec une méthode de point fixe.

Nous vérifions finalement par des expériences numériques l'efficacité de cette méthode sur des exemples de graphes, d'abord simples, puis plus complexes, avant de mettre en évidence l'efficacité de notre approche en la comparant aux méthodes standard.

Chapitre 3 : Un code Julia pour le contrôle optimal d'un réseau routier une structure similaire. La fin du chapitre entre un peu plus dans les détails propres au langage de programmation utilisé pour présenter des perspectives d'amélioration du code du point de vue de son empreinte mémoire et de sa rapidité, deux sujets particulièrement reliés en Julia.

Chapitre 4 : Outils pour les modèles réduits et le machine learning

La décomposition en modes propres permet de capturer la structure sous-jacente à un jeu de données et de filtrer les informations redondantes. Cette technique de compression est efficace pour reformuler un système d'équations de grande dimension N >> 1 tel que :

dρ(t) dt = F (ρ(t)), ρ(t) ∈ R N , (6)
en un système tronqué de taille réduite m << N :

d ρ(t) dt = F ( ρ(t)), ρ(t) ∈ R m , (7)
où ρ est obtenu par projection (compression) de ρ sur une base judicieusement choisie pour garantir la fiabilité des données réduites. Typiquement, nous travaillerons avec N = 800 et m ≤ 20 et noterons Φ ∈ R N ×n de sorte que :

ρ = Φ T ρ. (8)
Cette méthode dite d'interpolation empirique (EIM), qui existe aussi dans sa version dynamique (DEIM) pour les cas faiblement non linéaires, souffre néanmoins grandement de la forte non-linéarité du second membre du système.

Nous choisissons d'y remédier à l'aide de réseaux de neurones. Après en avoir exposé brièvement les CONTENTS CONTENTS Strasbourg University principaux aspects, nous rappelons une version de la descente de gradient qui permet de calculer de façon efficace l'erreur de prédiction du réseau par rapport à un grand jeu de données et qui est la base des méthodes utilisées dans les librairies standards de machine learning comme Adam ou RMS prop [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]: la descente de gradient stochastique [START_REF] Robbins | A stochastic approximation method[END_REF].

Chapitre 5 : Approches comparées d'hyper-réduction par apprentissage profond

Nous utilisons les bases réduites pour approcher une famille d'équations différentielles ordinaires non linéaires, indexée par un paramètre du modèle que l'on notera µ :

dρ µ (t) dt = F (ρ µ (t)), ρ µ (t) ∈ R N . (9)
Dans le chapitre 5, nous traiterons spécifiquement du cas où le modèle est donné par (4) sur une route à une entrée et une sortie, avec pour paramètre la densité maximale de voiture autorisée sur la route sortante. Le comportement de la solution numérique pour différentes valeurs du paramètre est illustré Figure 5, et l'on remarque que les basses valeurs du paramètre favorisent l'apparition de chocs tandis que les hautes valeurs génèrent des ondes de raréfactions.

En supposant que l'on puisse déduire ρ µ d'une représentation réduite ρµ ∈ R m de sorte que ρ µ ≈ Φ ρµ , la non-linéarité annulerait l'intérêt de l'utilisation d'un modèle réduit. Nous serions en effet ammenés à résoudre :

Φ T dρ µ (t) dt = d ρµ (t) dt = Φ T F (Φ ρµ (t)), ρ µ (t) ∈ R N , ρµ (t) ∈ R m , (10)
impliquant une évaluation du flux F en grande dimension.

Nous formons alors un réseau de neurones, paramétré par θ, dont le but sera d'apprendre à évaluer le flux

F / / F (ρ) Φ T ρ Φ O O Φ T •F •Φ
. . fortement non linéaire en n'utilisant que la base réduite (flèche bleue sur la Figure 6), c'est-à-dire de sorte que le calcul ne fasse pas intervenir la grande dimension du problème initial (flèche rouge sur la Figure 6).

Les deux approches utilisent le même réseau de neurones Fθ , un jeu de données ρµ (t) pour certaines valeurs de (t, µ) et ne diffèrent que par leur fonction de perte. Dans le premier cas, appelé "Neural Closure" (NC), nous cherchons à minimiser l'erreur entre le second membre prédit par le réseau de neurones Fθ ( ρµ , µ) et la compression Φ T (•) du second membre provenant du modèle de grande dimension F ( ρµ ; µ) :

L N C (θ) = error Fθ ( ρµ , µ), Φ T F ( ρµ ; µ) . (11)
Dans la seconde approche, appelée "Differentiable Physics" (DP), nous notons :

S 1 θ ( ρµ (t n )) := ρµ (t n ) + ∆t Fθ ( ρµ (t n ), µ) (12)
l'intégrateur en temps construit en utilisant le réseau de neurones comme flux d'intégration, et écrivons S K θ les composées itérées de S 1 θ nous transportant l'état réduit d'un instant t n à un instant t n+K . Nous souhaitons ainsi minimiser l'erreur commise par cet intégrateur sur K pas de temps en le comparant à la compression de la trajectoire de référence

Φ T ρ µ (t n ), • • • , Φ T ρ µ (t n+K ) : L DP K (θ) = error(S K θ ( ρµ (•)), ρµ (• + K)). (13)
C'est avec cette dernière approche que nous obtenons les résultats les plus prometteurs.

Contributions et communications

Cette thèse a donné lieu aux contributions et communications suivantes :

• Articles :

-Optimal scenario for road evacuation in an urban environment, preprint HAL (2023)

-Comparative approaches to hyper-reduction using deep learning, article en cours de rédaction -Numerical schemes for mixture theory models with filling constraint: application to biofilm ecosystems. Ce papier est soumis, et sera publié dans le journal ESAIM: Proceedings and Surveys.

• Orateur invité :

-

Tools for traffic modeling and analysis

This chapter first outlines road traffic model categories, then delves into macroscopic modeling, particularly nonlinear intersection coupling. It also previews the tools for analyzing the control problem and its numerical approximation. Emphasizing the significance of solutions' regularity and bounded variation, which plays a vital role in solution existence and iterative algorithm convergence toward optimal control.

Overview of existing traffic models

Depending on the specific scenario being studied, various approaches exist for modeling the phenomenon under consideration. These approaches vary in scale and the level of accuracy they aim to achieve.

At the microscopic scale, the behavior of individual vehicles can be described using Lagrangian microscopic models. These models focus on the trajectories, including speed and position, of each vehicle. Examples of such models include car-following models, as referenced in works such as [START_REF]Simulation des strassenverkehrsflusses[END_REF][START_REF] Gipps | A behavioural car-following model for computer simulation[END_REF][START_REF]Congested traffic states in empirical observations and microscopic simulations[END_REF], as well as cellular automaton models, as discussed in [START_REF] Krug | Universality classes for deterministic surface growth[END_REF][START_REF] Biham | Self-organization and a dynamical transition in traffic-flow models[END_REF][START_REF]A cellular automaton model for freeway traffic[END_REF]. Additionally, kinetic traffic models, introduced in the 1960s, take a mesoscopic approach, representing vehicles as distributions in position-velocity space. These models have been explored in works such as [START_REF] Prigogine | On a generalized boltzmann-like approach for traffic flow[END_REF][START_REF] Coscia | On the mathematical theory of vehicular traffic flow ii discrete velocity kinetic models[END_REF][START_REF] Delitala | Mathematical modeling of vehicular traffic: a discrete kinetic theory approach[END_REF][START_REF] Fermo | Fundamental diagrams for kinetic equations of traffic flow[END_REF][START_REF] Materne | An explicit kinetic model for traffic flow[END_REF][START_REF] Paveri-Fontana | On boltzmann-like treatments for traffic flow: A critical review of the basic model and an alternative proposal for dilute traffic analysis[END_REF][START_REF] Wegener | A kinetic model for vehicular traffic derived from a stochastic microscopic model[END_REF].

Moving up to a macroscopic scale, the flow of vehicles is treated as a continuous medium. One notable macroscopic model is the one introduced by Lighthill, Whitham, and Richards (commonly denoted as LWR) in the 1950s, as documented in [START_REF] Lighthill | On kinematic waves -ii -a theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF]. In the subsequent sections of this work, we will primarily use the LWR model, as its continuous representation of flow lends itself well to formulating dynamic programming problems.

It's worth noting that the original first-order model from the 1950s has been extended to second-order models resembling Euler-type systems with viscosity, as discussed in [START_REF] Payne | Model of freeway traffic and control[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF][START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF][START_REF]A non-equilibrium traffic model devoid of gas-like behavior[END_REF]. These extensions provide a more comprehensive understanding of traffic dynamics.

For a comprehensive introduction to this topic, see [START_REF] Garavello | Traffic flow on networks: conservation law models[END_REF].

Macroscopic modeling of road traffic on a network

To study road traffic on a macroscopic scale, we use a continuous model introduced by Lighthill, Whitham, and rise to numerous studies. The junction problem is of central importance, since it is necessary to make sense of solutions that are a priori multivalued, as well as to properly model the distribution of vehicles from one road to another. This section starts with a formal derivation of the model to give an intuition of the underlying physics.

Then, classical results of the domain are presented.

We write ρ(t, x) the number of cars per unit of length at position x ∈ R and time t ≥ 0. This quantity allows us to define the mass on any interval I := [a, b] as being:

m(t) = b a ρ(t, x) dx. (1.1)
This is a conserved quantity of the system, which means that:

m ′ (t) = b a ∂ t ρ(t, x) dx = 0. (1.2)
The flow of vehicles, denoted as f (ρ(t, x)), is the mass of vehicles that can pass through x at time t. Assuming no source or loss of vehicles, its variations at the boundaries of a road give us the variation of the mass within it, i.e.

m ′ (t) = b a ∂ t ρ(t, x) dx = ingoing flux -outoing flux = f (ρ(t, a)) -f (ρ(t, b)) = - b a ∂ x f (ρ(t, x)) dx.
We therefore obtain for some T ≥ 0 the following scalar conservation law:

∂ t ρ(t, x) + ∂ x f (ρ(t, x)) = 0, (t, x) ∈ (0, T ) × [a, b]. (1.3)
In [Lighthill, Whitham, Richards, 1955], the flux is given more specifically by:

f (ρ) = ρv(ρ) = v max ρ(1 - ρ ρ max ), (1.4)
with v max and ρ max the maximum speed and density allowed on the road. This choice is motivated by the fact that v(0) = v max and v(ρ max ) = 0. We drive indeed faster on an empty road than during a traffic jam, as illustrated in Figure 1.1.

Since this type of model is well known to produce discontinuous solutions in finite time, even for smooth initial data, we therefore need to weaken the notion of solution.

Definition 1 (Weak solution on a road). We say that ρ is a weak solution on road I if for any φ : R ≥0 × I → R smooth with compact support on R >0 × I, we have:

+∞ 0 I (ρ∂ t φ + f (ρ)∂ x φ) dxdt = 0. (1.5)
For the associated Cauchy problem to be well-posed, the search for weak solutions is not yet sufficient to guarantee uniqueness. The family of functions we need is the weak entropic ones, in other words, those that are Tools for traffic Strasbourg University physically relevant. These considerations are very standard and are discussed in detail in Coclite et al. [START_REF] Gabriellle | A fluid-dynamic traffic model on road networks[END_REF].

On a network

A road network of N r roads and N v junctions is modelled by a directed graph G = (V, E) whose edges

E = {I i := [a i , b i ] ⊂ R | 1 ≤ i ≤ N r },
represent the roads and the vertices

V = {J k = (i k 1 , • • • , i k n k ), (j k 1 , • • • , j k m k ) | 1 ≤ k ≤ N v },
the junctions, each junction J k having n k incoming and m k outgoing roads.

The LWR model is straightforwardly extended to describe the behaviour of the vehicles densities (ρ i ) i on roads I up to a time T ≥ 0, which readily yields the following Cauchy problem on the network:

∂ t ρ i (t, x) + ∂ x f i (ρ i (t, x)) = 0, (t, x) ∈ [0, T ] × [a i , b i ], ρ i (0, x) = ρ i,0 (x), x ∈ [a i , b i ], (1.6)
where flow and speed are given by:

f i (ρ i ) = ρ i v i (ρ i ), v i (ρ i ) = v max i 1 - ρ i ρ max i , (1.7

Vehicles distribution at junctions

Defining the coupling of the roads at junctions in an appropriate way is however a non-trivial problem, since the density becomes multivaluated with as many crossing characteristics as the arity of the intersection. This has given rise to numerous works [REF], and we will follow the approach of Coclite et al. [START_REF] Coclite | Traffic flow on a road network[END_REF] by linking the solution at junctions from its flow traces.

We write for each time t ∈ (0, T ) the flow traces γ(t) = γ L (t), γ R (t) of the flux at boundaries, such that:

Tools for traffic Strasbourg University • > • > • [a 1 γ L 1 b 1 ] γ R 1 [a 2 γ L 2 b 2 ] γ R 2 Ψ R 0 Ψ L 3 Figure 1.
2: Junction with one ingoing and one outoing road.

• γ R i (t) := f i ρ i (t, b - i )
is the flow going in the junction from the road i, located at the right1 of the ingoing road.

• γ L j (t) := f j ρ j (t, a + j ) is the flow going from the junction to the road j, located at the left of the outgoing road.

• (Ψ k (t)) k denote the flows at the junctions located next to graph's leaves, using ghost cells convention.

We provide an illustration of these notations in Figure 1.2, on a junction between one ingoing and one outoing road.

Considering a generic junction with n ingoing and m outgoing roads, Coclite et al. introduce the relation

n i=1 f i (ρ i (t, b i )) = n+m j=n+1 f j (ρ j (t, a j )), (1.8)
to ensure the conservation of the flux. This means that the number of cars entering the intersection is the same as the number leaving it. But it does not determine uniquely where they come from or where they are going at if there are several incoming or outgoing roads (i.e. for m, n ≥ 2). They therefore add constaints to the model such as assuming that there is a raw-stochastic distribution matrix A ∈ R m×n that links the input of a junction (γ R i ) 1≤i≤n to its output (γ L j ) n+1≤j≤n+m = A(γ R i ) 1≤i≤n . We then write

A := (α ji ) n+1≤j≤n+m 1≤i≤n , 0 < α ji < 1, n+m j=n+1 α ji = 1, (1.9)
with α ji being the percentage of drivers arriving from the i-th road that take the j-th outgoing road.

Finally, the main modelling assumption made to achieve coupling between roads is that drivers behave in such a way as to maximize flow, while respecting the intrinsic restrictions imposed, for example, by road size or speed limits. As a result, the flow at junctions will be an optimal compromise between maximum demand satisfaction and available supplies. These restrictions are computed as follows: by f i concavity, there indeed exists a unique

σ i ∈ [0, ρ max i ] such that f ′ i (σ i ) = 0.
In this point, f i reaches its maximum so we define (γ R,max i ) 1≤i≤n the maximum flows that can be obtained on incoming roads, i.e. the demand flux, Figure 1.3a, as:

γ R,max i = f i (ρ i (t, b i )), if ρ i (t, b i ) ∈ [0, σ i ], f i (σ i ), if ρ i (t, b i ) ∈ [σ i , ρ max i ], = f i (min{σ i , ρ i (t, b i )}). (1.10)
We define the same way the supply flux (γ L,max j ) n+1≤j≤n+m for the outgoing roads, see Figure 1.3b:

γ L,max j = f j (σ j ), if ρ j (t, a j ) ∈ [0, σ j ], f j (ρ j (t, a j )), if ρ j (t, a j ) ∈ [σ j , ρ max j ].
= f j (max{σ j , ρ j (t, a j )}).

(1.11) In order to write the optimal flow repartition, let us introduce the pavements:

Tools for traffic

Ω i := [0, γ R,max i ], 1 ≤ i ≤ n, Ω j := [0, γ L,max j ], n + 1 ≤ j ≤ n + m, (1.12)
which are used to define the following polytope:

Ω := {γ := (γ 1 , • • • , γ n ) ∈ Ω 1 × • • • × Ω n | Aγ ∈ Ω n+1 × • • • × Ω n+m }. (1.13)
The optimal trade-off principle stated earlier then leads to solve, at every junction and for each time, the following linear programming:

max γ∈Ω 1 • γ, (1.14)
where

1 := (1, • • • , 1) ∈ R n .
See again Coclite et al. [START_REF] Coclite | Traffic flow on a road network[END_REF] for the definition of weak entropic solutions at junctions.

Having defined the notion of solution on each roads and junctions, we are now able to define the solution on an entire network.

Definition 2 (Admissible solution on a network). Given initial datum ρi :

I i → R and boudary conditions ψ i : [0, +∞[→ R functions of L ∞ , we say that a collection of functions ρ = (ρ 1 , • • • , ρ Nr ) with ρ i : [0, +∞[×I i → R continuous as functions from [0, +∞[ into L 1 loc
, is an admissible solution if ρ i is a weak entropic solution to (1.6) on I i , and such that at each junction ρ is a weak solution and is an admissible weak solution in case of bounded variations.

Note that we're working in the context of bounded variations. This notion has several implications which will be discussed in section 1.3. For now, let's quickly recall the notion of Riemann problem, and how to solve it in this context of networks. This is interesting both for theoretical reasons with the existence of a solution to the Cauchy problem on a network by the so-called wavefront tracking algorithm [START_REF] Coclite | Traffic flow on a road network[END_REF], and for numerical reasons since it allows to construct finite volume schemes.

Tools for traffic

Strasbourg University

Analysis of the Riemann problem

The existence and uniqueness for the Cauchy problem on network are established in [START_REF] Coclite | Traffic flow on a road network[END_REF] using the wavefront tracking algorithm, which relies on approximating the solution by piecewise constant functions and solving each local problem using the fact that in hyberbolic dynamics, informations travels at finite speed. We will first recall Riemann problems on single roads before giving the extensions on junctions and networks.

Riemann problem at roads

We aim to solve for each road the following kind of Cauchy problem, with Heavyside-like initial condition:

         ∂ t ρ(t, x) + ∂ x f (ρ(t, x)) = 0, ρ(0, x) = ρ 0 (x) =    ρ L , x < 0 ρ R , x > 0 (1.15)
This is a very well known problem, see [START_REF] Bressan | Hyperbolic systems of conservation laws[END_REF] for instance, and the solution is given in the following theorem.

Theorem 1. If ρ L < ρ R then f ′ (ρ L ) > f ′ (ρ R
) and we have an entropic shock propagating at the speed σ given by the Rankine-Hugoniot relation:

σ = f (ρ L ) -f (ρ R ) ρ L -ρ R = v max 1 - ρ L + ρ R ρ max , (1.16)
ρ(t, x) = r(ξ) =    ρ L , ξ < σ, ρ R , ξ > σ.
(1.17)

Conversely, if ρ L > ρ R then f ′ (ρ L ) < f ′ (ρ R
) and we are in the case of a rarefaction wave, given by:

ρ(t, x) = r(ξ) =          ρ L , ξ < f ′ (ρ L ), g(ξ), f ′ ρ L < ξ < f ′ (ρ R ), ρ R , f ′ (ρ R ) < ξ, (1.18)
where

g(ξ) = ρmax 2 1 -ξ vmax is such that f ′ (g(ξ)) = ξ.

Riemann problem at junction

The Riemann problem at a junction J := ((

I 1 , • • • , I n ), (I n+1 , • • • , I n+m )
) is given by Equation 1.6 on each road (I k ) 1≤k≤n+m with respect to the optimal coupling constraints from Equation 1.14, where the initial condition is a constant density on each road. According to the aformentioned modeling hypothesis, there exists a unique weak solution of this problem in the sense given by [ [START_REF] Coclite | Traffic flow on a road network[END_REF], Theorem 3.1].

Cauchy problem on network

Using the wavefront tracking strategy developed in [START_REF] Dafermos | Polygonal approximations of solutions of the initial value problem for a conservation law[END_REF], knowledge of the Riemann problem is sufficient to obtain unique solutions to generic Cauchy problems. In fact, the idea of the wavefront tracking algorithm is to approximate the solution by a family of piecewise constant functions, solve the Riemann problems at each interface and then pass to the limit with the appropriate topology. It is important to note that uniqueness is only given in the following sense:

The solution given by the convergence of the wavefront tracking algorithm does not depend on the choice of the piecewise constant approximations.

However, there is no guarantee that another type of strategy will give the same solution. As far as we known, the latter problem remains open.

Definition 3 (Wavefront tracking approximation). The wavefront tracking algorithm reads as follows:

• let ρ be a piecewise constant map defined on the network

• we solve riemann problems on each road at the discontinuity points of ρ as well as the riemann problem at junctions from the traces of ρ

• for each riemann problem, if the obtained solution is a shockwave, do nothing if the obtained solution is a rarefaction, we spatially discretize it in order to obtain discontinuities travelling at Rankine-Hugoniot speed. We thus deal with shockwaves only. An illustration of this step is reproduced Figure 1.4 and comes from [START_REF] Liard | Analysis and numerical solvability of backward-forward conservation laws[END_REF] • we obtained a new piecewise constant solution, to wich we apply the sames steps as above.

Theorem 2 (Wavefront tracking approximation property). Let ρ ν be a wavefront tracking approximation of Equation 1.3 with grid size given by ν. If this squence converges in L 1 loc when ν → 0, then the limit is a weak entropic solution of Equation 1.3.

Theorem 3 [START_REF] Coclite | Traffic flow on a road network[END_REF]. Let f be the flow from the LWR model, and consider a network in which each junction has at most two ingoing and two outoing roads. Let ρ ∈ L 1 loc be an initial datum. Then, for any time T > 0, we have the following properties:

Tools for traffic

Strasbourg University

• All approximating sequences generated by the wave-front tracking algorithm yields a unique solution ρ to System (1.15) on [0, T ], such that ρ(0, •) = ρ, obtained as the limit of piecewise-constant solutions whose flows derivatives never cancel.

• If the initial datum ρ belongs to L 1 (0, T ) then ρ(t, •) belongs to L 1 .

• Let ρ(t, •), ρ ′ (t, •) ∈ L 1 for a.e. time t ≥ 0, both obtained by wavefront tracking approximation. Then:

ρ(t, •) -ρ ′ (t, •) L 1 ≤ ρ(0, •) -ρ ′ (0, •) L 1 .

About BV spaces

We introduce now the space of functions with bounded variation, which is important both for the existence of the Cauchy problem for the LWR model and the existence of an optimal control. Given a function f belonging to L 1 ([0, T ]), its total variation on [0, T ] is defined as

TV(f ) = sup T 0 f (t)ϕ ′ (t) dt, ϕ ∈ C 1 c ([0, T ]), ∥ϕ∥ L ∞ ([0,T ]) ≤ 1 .
We denote by BV (0, T ) the space of functions of bounded variations,

BV (0, T ) := f ∈ L 1 (0, T ), T V (f ) < +∞ ,
and, endowed with the norm

∥•∥ BV (0,T ) = T V (•) + ∥•∥ L 1 (0,T ) ,
the space BV (0, T ), ∥•∥ BV (0,T ) is Banach. It possesses the following usefull properties:

• monotonicity: any function f in BV (0, T ) can be rewritten as f = g -h, where g and h are increasing functions, and conversely. In particular, any monotone function is in BV (0, T ).

• density: BV (0, T ) = C ∞ c (0, T ), for the topology associated to the ∥•∥ BV (0,T ) norm.

• compactness: This property is a key ingredient of the proof of Theorem 5 in order to obtain convergence of the control from the convergence of a minimizing sequence of the cost function.

Theorem 4. (Helly selection theorem). for M > 0, the function set {f ∈ BV (0, T ) | ∥f ∥ BV ≤ M } is a relatively compact subset of L 1 (0, T ), i.e. for every sequence

(f n ) n in BV (0, T ) such that ∃M ∈ R, ∀n, T V (f n ) ≤ M , then there exists f ∈ L 1 (0, T ) such that f n L 1
-→ f , up to a subsequence.

• lower semi-continuity:

BV (0, T ) ∋ f n L 1 loc -→ f ⇒ lim inf n→+∞ T V (f n ) ≥ T V (f ).
This property is also used in the proof mentionned above, to conclude that the limit control from the minimizing sequence solves the optimal control problem.

From what we've stated, BV functions are essentially an algebra of discontinuous functions that are almost everywhere differentiable. This makes this space very convenient for working with discontinuous solutions of PDEs, such as the shocks generated by the nonlinear scalar conservation laws we deal with in this thesis. Not only does this give us the existence of weak solutions for the Cauchy problem on network with BV initial data, but regularization properties also help to avoid the so-called "chattering" phenomenon in control theory. Indeed, having a solution in BV is morally having its gradient bounded for the L 1 -norm. See chapter 2 for more on this subject.

We here give the following well known inclusions, true on any compact interval, to provide a global view:

C 1 ⊂ Lip ⊂ AC ⊂ BV ⊂ D -a.e.,
where AC and D -a.e. are the sets of absolutely continuous and almost everywhere differentiable functions respectively.

Discrete framework: while the above considerations are of great help in proving most of the theoretical guarantees of our work, it is also important to highlight the behavior of the total variation on the discretizations of our problems. The discretization of a function

f : t ∈ [0, T ] → R on a mesh 0 = t 0 < t 1 < • • • < t N = T
leads us to consider functions of the type:

f disc (t) := N -1 n=0 f n 1 [tn,t n+1 [ (t),
where we write f n := f (t n ), and a direct computation yields:

T V (f disc ) = N -1 n=0 | f n+1 -f n | .
Since we will often need to manipulate differentiable expressions in order to perform, for instance, a gradient Tools for traffic descent, we introduce for some ν > 0 the regularized total variation, more extensively studied in subsection 2.4.1:

T V ν (f disc ) = N -1 n=0 a ν (f n+1 -f n ).
It uses the following well-known smoothed absolute value, illustrated for different values of ν in Figure 1.6:

a ν (x) := x 2 + ν 2 .
Chapter 2

Optimal scenario for road evacuation in an urban environment

The following chapter relates an article written in collaboration with Emmanuel Franck (INRIA TONUS, Strasbourg), Laurent Navoret (INRIA TONUS, Strasbourg), and Yannick Privat (IECL, Nancy).

Abstract: How to free a road from vehicle traffic as efficiently as possible and in a given time, in order to allow for example the passage of emergency vehicles? We are interested in this question which we reformulate as an optimal control problem. We consider a macroscopic road traffic model on networks, semi-discretized in space and decide to give ourselves the possibility to control the flow at junctions. Our target is to smooth the traffic along a given path within a fixed time. A parsimony constraint is imposed on the controls, in order to ensure that the optimal strategies are feasible in practice.

We perform an analysis of the resulting optimal control problem, proving the existence of an optimal control and deriving optimality conditions, which we rewrite as a single functional equation. We then use this formulation to derive a mixed algorithm interpreting it as a mix between two methods: a descent method combined with a fixed point method. We verify with numerical experiments the efficiency of this method on examples of graphs, first simple, then more complex. We highlight the efficiency of our approach by comparing it to standard methods.

We propose a code implementing this approach in the Julia language, available to any user at the following link: Keywords: Traffic network, Optimal Control, Fluid model, Hyperbolic PDE, optimization methods AMS classification: 35Q49, 65M08, 49K15, 49K30.

Introduction

Road traffic modeling. With the concentration of populations in cities where transport flows are constantly increasing, road traffic modeling has become a central issue in urban planning. Whether it is to configure traffic lights, to adapt public transport offers or to manage in an optimal way situations of high congestion. The knowledge of the behavior of a road network is an important issue for the management of a crisis in an urban environment. It is indeed necessary to quickly predict traffic, especially traffic jams, in order to organize rescue operations. Before creating decision support tools for crisis management involving road traffic, a first step is to develop a method for monitoring road traffic.

Optimal scenario for road evacuation in an urban environment

Strasbourg University

In this work, we will focus on the evacuation of a traffic lane in finite time. It can be, for example, an axis that we wish to free in order to provide access to first aid. We approach this question in the form of an optimal control problem of a graph where each edge corresponds to a traffic lane. The main objective of this work is the determination of a prototype algorithm answering this question.

State of the art. Depending on the situation of interest, several types of approaches co-exist to model this phenomenon, depending on the scale and the desired accuracy. At the microscopic scale, trajectories (speed and position) of each vehicle can be described through Lagrangian microscopic models such as car-following models [START_REF] Gipps | A behavioural car-following model for computer simulation[END_REF][START_REF]Simulation des strassenverkehrsflusses[END_REF][START_REF]Congested traffic states in empirical observations and microscopic simulations[END_REF] or cellular automaton models [START_REF] Biham | Self-organization and a dynamical transition in traffic-flow models[END_REF][START_REF] Krug | Universality classes for deterministic surface growth[END_REF][START_REF]A cellular automaton model for freeway traffic[END_REF]. Kinetic traffic models were also introduced in the 60's [START_REF] Coscia | On the mathematical theory of vehicular traffic flow ii discrete velocity kinetic models[END_REF][START_REF] Delitala | Mathematical modeling of vehicular traffic: a discrete kinetic theory approach[END_REF][START_REF] Fermo | Fundamental diagrams for kinetic equations of traffic flow[END_REF][START_REF] Materne | An explicit kinetic model for traffic flow[END_REF][START_REF] Paveri-Fontana | On boltzmann-like treatments for traffic flow: A critical review of the basic model and an alternative proposal for dilute traffic analysis[END_REF][START_REF] Prigogine | On a generalized boltzmann-like approach for traffic flow[END_REF][START_REF] Wegener | A kinetic model for vehicular traffic derived from a stochastic microscopic model[END_REF] and deal with vehicles on a mesoscopic scale in the form of distributions in position-velocity space.

Finally, on a macroscopic scale, the vehicle flow is then considered as a continuous medium. This is the case of the famous model introduced by Lighthill, Whitham and Richards (denoted LWR in what follows) [START_REF] Lighthill | On kinematic waves -ii -a theory of traffic flow on long crowded roads[END_REF][START_REF] Richards | Shock waves on the highway[END_REF].

We will use it in the rest of this article, since the continuous description of the flow is particularly well adapted to the writing of a dynamic programming problem. It is notable that this first order model introduced in the 50's was later extended to Euler-type systems with viscosity, of second order [START_REF] Payne | Model of freeway traffic and control[END_REF][START_REF] Whitham | Linear and nonlinear waves[END_REF][START_REF] Aw | Resurrection of "second order" models of traffic flow[END_REF][START_REF]A non-equilibrium traffic model devoid of gas-like behavior[END_REF].

In this paper, we introduce and analyze an optimal control problem modeling an evacuation scenario for an axis belonging to a road network. A close control model has already been introduced in [START_REF]Optimal control for traffic flow networks[END_REF], but under restrictive assumptions about the flow regime in order to avoid congestion phenomena, which is what we want to avoid. It has been in particular highlighted that the flow may be not invertible with respect to control parameters, meaning that density at the nodes may not be reconstructed from boundary conditions. In [START_REF] Goatin | Speed limit and ramp meter control for traffic flow networks[END_REF], a method to manage variable speed limits combined with coordinated ramp metering within the framework of the Lighthill-Whitham-Richards

(LWR) network model has been introduced. They consider a "first-discretize-then-optimize" numerical approach to solve it numerically. Note that the optimal control problem we deal with in the following has some similarities with the one proposed by the authors of this article, but differs in several respects: our control variables aim to control not just the ramps metering, but an entire axis. Furthermore, our approach to solving the underlying optimal control problem is quite different, in particular our use of optimality conditions leading to an efficient solution algorithm.

In [START_REF] Reilly | Adjoint-based optimization on a network of discretized scalar conservation law PDEs with applications to coordinated ramp metering[END_REF], a similar control problem is discussed, which is addressed using a piecewise-linear flux framework driving to a much simplified adjoint system involving piecewise constant Jacobian matrices in time and space.

This approach enables faster gradient computation, allowing them to tackle real-world simulations related to ramp-metering configuration. In our study, we opted to use the model introduced in [START_REF] Coclite | Traffic flow on a road network[END_REF] and dealt with the non-linear flux using automatic differentiation techniques to compute the gradient. This kind of approach has also been used in [START_REF] Bayen | Adjoint-based control of a new eulerian network model of air traffic flow[END_REF] for air traffic flow with a modified LWR-based network model.

We also mention several contributions aiming at dealing with concrete or real time cases. The paper [START_REF] Fügenschuh | Combinatorial and Continuous Models for the Optimization of Traffic Flows on Networks[END_REF] deals with this kind of control problem for simplified nonlinear (essentially without considering congestions) and linear formulations of the model, considering large scale networks. In [START_REF] Kurzhanskiy | Active traffic management on road networks: a macroscopic approach[END_REF], an algorithm for real-time traffic state estimation and short-term prediction is introduced.

Several other close control problems have been investigated in the past: in [START_REF] Göttlich | Modeling and optimizing traffic light settings in road networks[END_REF], the authors use switching controls to deal with traffic lights at an 8×4 junction in a piecewise linear flow framework. The paper [START_REF] Ancona | On the optimization of conservation law models at a junction with inflow and flow distribution controls[END_REF] deals with an optimal control problem of the same kind, formulated from the continuous LWR model, solved

Optimal scenario for road evacuation in an urban environment Strasbourg University numerically using an heuristic random parameters search. We refer to [START_REF] Ancona | Analysis and control on networks: Trends and perspectives[END_REF][START_REF] Bayen | Control Problems for Conservation Laws with Traffic Applications: Modeling, Analysis, and Numerical Methods[END_REF] for a broad view on this topic.

Regarding now controlled microscopic models, let us mention [START_REF] Malena | Traffic estimation and mpc-based traffic light system control in realistic real-time traffic environments[END_REF] where an approach based on model predictive control (MPD) is developed and [START_REF] Baumgart | Optimal control of traffic flow based on reinforcement learning[END_REF] using a reinforcement learning algorithm. In the book [START_REF] Martin Treiber | Traffic Flow Dynamics[END_REF], the last chapter is devoted to the modeling of various optimization problems related to road traffic. There has also been some research into the optimal switching of the traffic lights to maximize the traffic flow using a mixed-integer model [START_REF] Göttlich | Modeling and optimizing traffic light settings in road networks[END_REF]. Other works have looked at traffic regulation by comparing and discussing the use of lights and circles [START_REF] Chitour | Traffic circles and timing of traffic lights for cars flow[END_REF].

Concerning now the control of mesoscopic and macroscopic models, there are software contributions seeking in particular to bring real time answers and guidance tools in case of accident. Let us also mention [START_REF] Jayakrishnan | An evaluation tool for advanced traffic information and management systems in urban networks[END_REF][START_REF] Contram | Structure of the model[END_REF] dealing in part with traffic prediction and control problems.

This article uses a controlled model close to the one studied in [START_REF]Optimal control for traffic flow networks[END_REF] which is devoted to the control of the macroscopic LWR model on a directed graph using the framework introduced in [START_REF] Coclite | Traffic flow on a road network[END_REF]. However, we are interested in taking into account some properties of the model that were not addressed in [START_REF]Optimal control for traffic flow networks[END_REF], such as congestion phenomena.

Hence, we will not focus on individual paths since we are interested in the overall dynamics with a focus on congestion. Particular attention is paid to the modeling of junctions as they introduce a nonlinear coupling between roads.

Main objectives. Our goal is to design a numerical method to control the flow of vehicles in all fluid regimes, including saturated regimes/congestions, using dynamical barriers at each road in the graph.

We thus propose a road traffic model in the form of a controlled graph at each junction. This models an urban area structured by roads. For practical reasons, even if we wish to take into account possible congestion/shocks in the system, we position ourselves in a differentiable framework by adding weak diffusion terms in the system, through a semi-discretization in space. Hence, the resulting problem consists in minimizing the vehicle density at the final time on a given route under ODE constraint, obtained by finite volume semi-discretization of the LWR system. We then enrich this problem with a constraint on the number of simultaneously active roadblocks in order to take into account staffing issues.

Plan of this article. While the macroscopic model describing the flow is very standard (Lighthill-Whitham-Richards, 1955) the problem of flow distribution at junctions, inspired by more recent works [START_REF] Gabriellle | A fluid-dynamic traffic model on road networks[END_REF], [START_REF] Coclite | Traffic flow on a road network[END_REF], realizes a nonlinear coupling between the different edges. The distribution of vehicles at junctions is thus modeled by an optimal allocation process that depends on the maximum possible flows at the junctions, through a linear programming (LP) problem targeting the maximization of incoming flows. In order to influence the traffic flow, we introduce control functions that are defined at each road entrance and acts as barriers by weighting the capacity of a road leaving a junction to admit new vehicles.

After investigating the existence of optimal control, we then derive necessary optimality conditions that we reformulate in an exploitable way and introduce an optimization algorithm based on a hybrid combination of a variable step gradient method and a well-chosen fixed point method translating the optimality conditions. We then illustrate the efficiency of the introduced method using examples of various graphs modeling, in particular, a traffic circle or a main road surrounded by satellite roads.

Our contribution can be summarized in the form of a computational code allowing to process complex graphs, Optimal scenario for road evacuation in an urban environment Strasbourg University but nevertheless of rather small dimension. This code, written in the Julia language, is available and usable at the following link:

https://github.com/mickaelbestard/TRoN.jl Notations. Throughout this article, we will use the following notations:

• For x ∈ R n , x + will denote the positive part of x, namely x + = max{x, 0 R n }, the max being understood component by component;

• BV(0, T ) denotes the space of all functions of bounded variation on (0, T );

• W 1,∞ (0, T ) denotes the space of all functions f in L ∞ (0, T ) whose gradient in the sense of distributions also belongs to L ∞ (0, T );

• ∥ • ∥ R d (resp. ⟨•, •⟩ R d ) stands for the standard Euclidean norm (resp. inner product) in R d ; • Let F : R p × R q → R m and (x 0 , y 0 ) ∈ R p × R q .
We will denote by ∂ x F (x 0 , y 0 ), the Jacobian matrix of the mapping x → F (x, y 0 ). When no ambiguity is possible, we will simply denote it by ∂ x F ;

• N r : number of roads in the model;

• N c number of mesh cells per road for the first-order Finite Volumes scheme;

• N T : number of time discretization steps in the numerical schemes;

A controlled model of trafic flow

In this section, we first present the traffic flow model on a road network and its semi-discretization. Then we introduce the precise control model: the control at junctions will translates regulation action using traffic signs or traffic lights.

Traffic dynamics on network without control

The road network is a directed graph of N r roads, with N r ∈ N * , where each edge corresponds to a road and each vertex to a junction. The roads will all be denoted as real intervals

[a i , b i ] for some i ∈ 1, N r . Of course,
this writing does not determine the topology of the graph. It is necessary to define, for each junction, the indices of the incoming and outgoing roads in order to characterize the directed graph completely.

We are interested in the time evolution of the densities on each road. On the i-th road, the evolution of the density ρ i (t, x) is provided by the so-called LWR model

               ∂ t ρ i (t, x) + ∂ x f i (ρ i (t, x)) = 0, (t, x) ∈ (0, T ) × [a i , b i ], ρ i (0, x) = ρ 0 i (x), x ∈ [a i , b i ], f i (ρ i (t, a i )) = γ L i (t), t ∈ (0, T ), f i (ρ i (t, b i )) = γ R i (t), t ∈ (0, T ), (LWR)
Optimal scenario for road evacuation in an urban environment

Strasbourg University

where the flux f i (ρ) is given by

f i (ρ) = ρ v max i 1 - ρ ρ max i , (2.1)
where v max i and ρ max i denotes respectively the maximal velocity and density allowed on the i-th road. We will consider without loss of generality that

ρ max i = v max i = 1.
The initial density is denoted ρ 0 i (x) and γ L i and γ R i are functions allowing to prescribe the flux at the left and right boundaries of the interval.

The fluxes γ at junctions are determined using the model introduced in [START_REF] Coclite | Traffic flow on a road network[END_REF]. Considering a junction J, the sets of indices corresponding to the in-going and out-going roads are denoted J in and J out . We assume that there is a statistical behavior matrix A J defined by

A J := (α ji ) (j,i)∈Jout×J in , 0 < α ji < 1, j∈Jout α ji = 1,
with α ji the proportion of vehicles going to the j-th outgoing road among those coming from the i-th ingoing road. Outgoing fluxes have to satisfy the relation

(γ L j ) j∈Jout = A J (γ R i ) i∈J in , (2.2)
and we have the balance between in-going and out-going fluxes:

i∈J in γ R i = j∈Jout γ L j .
Let us note that we have natural constraints on the fluxes. Indeed the LWR flux function f i reaches its maximum at ρ max i /2. Therefore, ingoing fluxes γ R i have to be smaller than the following upper bound:

γ R,max i = f i (min{ρ i (t, b i ), ρ max i /2}),
which takes into account reduced demands when the density is lower than ρ max i /2 in the ingoing road. Similarly outgoing fluxes γ L j have to be smaller than the upper bound:

γ L,max j = f j (max{ρ j (t, a j ), ρ max j /2}),
which takes into account reduced capacities when the density is larger than ρ max i /2 in the outgoing road. We refer to [START_REF] Coclite | Traffic flow on a road network[END_REF] for details. Consequently, the fluxes belong to following set:

Ω J =    (γ R i ) j∈J in , (γ L j ) i∈Jout ∈ i∈J in 0, γ R,max i j∈Jout 0, γ L,max j with (γ L j ) j∈Jout = A J (γ R i ) i∈J in   
Then we assume that drivers succeed in maximizing the total flow. Consequently the fluxes are solution to the following Linear Programming (LP) problem:

max ((γ R i ) j∈J in ,(γ L j ) j∈J out )∈ΩJ i∈J in γ R i , (2.3)
Optimal scenario for road evacuation in an urban environment

Strasbourg University

This definition is valid at any junction of the network. Note that this problem does not have necessarily a unique solution in the case where there are more incoming roads than outgoing roads. In that case, a priority modeling assumption should be added to select one solution. We refer to [START_REF] Bretti | Numerical approximations of a traffic flow model on networks[END_REF] for more details. In the sequel, we suppose that this linear programming problem have a unique solution solution and we formally write:

γ(t) = ϕ LP (ρ(t)),
for the simultaneous resolution of the linear programming problems at all the junctions of the network. Note that the dependency on ρ(t) results from the definition of the upper bounds involved in the definition of the sets Ω J .

In the following, we will only consider networks with at most two ingoing roads and two outgoing roads, i.e.

either 1 × 1, 1 × 2, 2 × 1 or 2 × 2 junctions.
This gives existence and stability of the LWR Cauchy problem [START_REF] Coclite | Traffic flow on a road network[END_REF],

and has the advantage that the function ϕ LP can be explicitly provided [START_REF] Shi | A maximum-principle-satisfying finite volume compact-weno scheme for traffic flow model on networks[END_REF] while it already enables to model a huge variety of road networks.

Control at junctions

We introduce the vector of controls u = (u i (•)) 1≤i≤Nr ∈ L ∞ ([0, T ], R Nr ) at every road entrance. Note that we implicitly assume that every road junction is controlled, but our model allows without any difficulty to neutralize some controls in order to model the fact that only some junctions are controlled.

We interpret each control u j as a rate, assuming that at each time t ∈ (0, T ), the maximum flow out a junction and going into road j is weighted by a coefficient u j (t) ∈ [0, 1]. This allow in particular to keep valid all well-posedness considerations mentioned in [START_REF] Coclite | Traffic flow on a road network[END_REF]. At junction J, we thus define the following polytope of constraints

Ω J (u) = (γ R i ) j∈J in , (γ L j ) j∈Jout ∈ i∈J in 0, γ R,max i j∈Jout 0, (1 -u j )γ L,max j with (γ L j ) j∈Jout = A J (γ R i ) i∈J in
This initial model is not completely satisfactory since a full control on one outgoing road would result on zero outgoing fluxes for all the ingoing roads. Indeed relation (2.2) implies that ingoing fluxes are linear combination of outgoing fluxes with positive weights. This is definitely not the desired behavior as we would expect that the trafic flow would be deviated to the uncontrolled outgoing roads. To solve this problem, we choose to make the statistical behavior matrix A J also dependent on the control u. More precisely, we ask that, as soon as a road entry is fully controlled and the other are not controlled, the proportion of vehicles entering the roas is set to 0.

To illustrate this point, let us describe the control of 1 × 2 junction, i.e. with one incoming and 2 outgoing roads.

Then the traffic distribution matrix A writes:

A(u) = α(u) 1 -α(u) ,
where α(u) denotes the proportion of vehicles going in the first outgoing road. Then, denoting u = (u 1 , u 2 ) the controls of the two outgoing roads, the function α(u 1 , u 2 ) has to be chosen such that:

(a) if the first outgoing road is fully controlled and the second is not controlled, then α(u) has to vanish, resulting

Optimal scenario for road evacuation in an urban environment Strasbourg University in the condition: α(1, 0) = 0, and then the entire vehicle flow goes in the second outgoing road.

(b) inversely, if the second outgoing road is fully controlled and the first one is not controlled, then 1 -α(u)

has to vanish, resulting in the condition: 1 -α(0, 1) = 0, and then the entire vehicle flow goes in the first outgoing road.

We also ask for the following additional property:

(c) the parameters are not modified if the outgoing roads are equally controlled, which results in the condition:

∀u ∈ [0, 1], α(u, u) = ᾱ,
where ᾱ is a given distribution coefficient.

In Appendix 2.6, we provide an explicit construction of such a function. We also treat the case of 2 × 2 junctions.

Note that junctions with only one outgoing roads (1 × 1 and 2 × 1) do not require such modification.

Like in the previous section, the simultaneous resolution of the Linear Porgramming problems at each junctions junction is now denoted:

γ = ϕ LP (ρ, u)
The explicit expressions of ϕ LP in the set of cases we deal with are provided in Appendix 2.6. Beyond the explicit expression of ϕ LP (ρ, u), what matters is that ϕ LP is a Lipschitz function with respect to (ρ, u). In the following, in order to use tools of differentiable optimization, we will consider a C 1 approximation of this function. This issue is commented at the end of Appendix 2.6. According to these comments, we will assume from now on:

The function ϕ LP is Lipschitz, and C 1 with respect to its second variable u.

(H ϕ LP )

Semi-discretized model

For algorithmic efficiency reasons, we prefer to be able to define the sensitivity of the different data of the problem with respect to the control. Since the model (LWR) is known to generate possible irregularities in the form of shocks, we have decided to introduce regularity through a semi-discretization of the model in space. The relevance of this choice will be discussed in the concluding section.

Hence, we discretize the model (LWR) with a first-order Finite Volume (FV) scheme. We consider N c ∈ N * mesh cells per road: the discretization points on the i-th roads are denoted

a i = x i,1/2 < x i,3/2 < . . . < x i,Nc-1/2 = b i and the space steps △x i,j = x i,j+1/2 -x i,j-1/2 .
Then the discrete densities and the boundary fluxes are denoted:

ρ(t) = ρ i,j (t) 1≤i≤Nr 1≤j≤Nc =   1 △x i,j x i,j+ 1 2 x i,j-1 2 ρ(t, x) dx   1≤i≤Nr 1≤j≤Nc , γ(t) = γ L 1 (t), γ R 1 (t), . . . , γ L Nr (t), γ R Nr (t) .
With a slight abuse of notation, we have written similarly the discrete variable as the continuous one, since we will essentially work on the semi-discretized model. Then the semi-discretized dynamics is given by the Optimal scenario for road evacuation in an urban environment Strasbourg University

differential equation          dρ dt t = f F V (ρ(t), γ(t)) , t ∈ (0, T ) γ(t) = ϕ LP (ρ(t), u(t))) , t ∈ (0, T ) ρ(0) = ρ 0 , (LWR-sd)
where the finite volume flow at the j-th mesh cell of the i-th road reads

f F V (ρ, γ) ij =          -1 △x ij F i (ρ i,j , ρ i,j+1 ) -γ L i , if j = 1, -1 △x ij (F i (ρ i,j , ρ i,j+1 ) -F i (ρ i,j-1 , ρ i,j )) , if 2 ≤ j ≤ N c -1, -1 △x ij γ R i -F i (ρ i,j-1 , ρ i,j ) , if j = N c , (2.4)
with F i (u, v) the so-called local Lax-Friedrich numerical flux given by

F i (u, v) = f i (u) + f i (v) 2 -max |f ′ i (u)|, |f ′ i (v)| (v -u) 2 .
We refer to [START_REF] Leveque | Numerical methods for conservation laws[END_REF] for an introduction to Finite Volume approximations.

Conclusion: an optimal control problem

We are interested in an approximate controllability problem which consists in emptying a given route as much as possible for a fixed end time T > 0. Let us introduce χ path ⊂ 1, N r , a set of indices corresponding to the route we wish to empty in a time T .

We would like to minimize a functional with respect to the control u, representing the sum of all the densities on this path, in other words

C T (u) = i∈χ path Nc j=1 ρ i,j (T ; u) = c • ρ(T ; u),
where ρ(t; u) denotes the solution to the controlled system (LWR-sd) and c = (c i ) 1≤i≤Nr the vector defined by

c i = 1 if i ∈ χ path and c i = 0 else.
Of course, it is necessary to introduce a certain number of constraints on the sought controls, in accordance with the model under consideration, and to model that the obtained control is feasible in practice. We are thus driven to consider the following constraints:

(i) 0 ≤ u i (•) ≤ 1, meaning that at each time, the control is a vehicle acceptance rate on a road;

(ii) Nr i=1 u i (•) ≤ N max : we impose at each time a maximum number of active controls in order to take into account the staff required for roadblocking.

(iii) Regularity: we will assume that u is of bounded variation, and write u ∈ BV([0, T ]; R Nr ). This constraint Optimal scenario for road evacuation in an urban environment Strasbourg University involves the total variation 1 of the control and models that the roadblock is supposed not to "blink" over time. From a control point of view, we aim at avoiding the so-called chattering phenomenon.

The first constraint above will be included in the set of admissible controls. Concerning the other two constraints, we have chosen to include them as penalty/regularization terms in the functional. Of course, other choices would be quite possible and relevant.

Let N max ∈ N * be an integer standing for the maximal number of active controls for this problem and θ = (θ S , θ B ) ∈ R 2 be two non-negative parameters. According to all the considerations above, the optimal control problem we will investigate reads: inf

u∈U ad J (u), (P θ )
where the admissible set of controls is defined by

U ad = L ∞ ([0, T ], [0, 1] Nr ), (2.5)
and the regularized cost functional J writes

J (u) = C T (u) + θ S S(u) + θ B B(u), (2.6)
where S(u) denotes the regularizing term modeling the limitation on the number of active controls and B(u) is the total variation of u in time, namely

S(u) = 1 2 T 0 Nr i=1 u i (t) -N max 2 + dt and B(u) = Nr i=1
TV(u i ).

(2.7)

2.3 Analysis of the optimal control problem (P θ )

In this section, we will investigate the well-posedness and derive the optimality conditions of Problem (P θ ).

These conditions will form the basis of the algorithms used in the rest of this study.

Well-posedness of Problem (P θ )

We follow the so-called direct method of calculus of variations. The key point of the following result is the establishment of uniform estimates of ρ with respect to the control variable u, in the W 1,∞ (0, T ) norm.

Theorem 5. Let us assume that θ B is positive. Then, Problem (P θ ) has a solution.

Proof. Let (u n ) n∈N be a minimizing sequence for Problem (P θ ). Observe first that, by minimality, the sequence (J (u n )) n∈N is bounded, and therefore, so is the total variation (B(u n )) n∈N . We infer that, up to a subsequence, (u n ) n∈N converges in L 1 (0, T ) and almost everywhere toward an element u * ∈ BV(0, T ). Since the class U ad is closed for this convergence, we moreover get that u * ∈ L ∞ ([0, T ], [0, 1] Nr ).

1 Given a function f belonging to L 1 ([0, T ]), the total variation of f in [0, T ] is defined as

TV(f ) = sup T 0 f (t)ϕ ′ (t) dt, ϕ ∈ C 1 c ([0, T ]), ∥ϕ∥ L ∞ ([0,T ]) ≤ 1 .
Optimal scenario for road evacuation in an urban environment

Strasbourg University

For the sake of readability, we will denote similarly a sequence and any subsequence with a slight abuse of notation.

In what follows, it is convenient to introduce the function g : R Nc×Nr × [0, 1] Nr defined by g(ρ, u) = f F V (ρ, ϕ LP (ρ, u)), so that ρ solves the ODE system ρ ′ = g(ρ, u). Let us set ρ n := ρ(• ; u n ).

Step 1: the function g is continuous. Indeed, note that the function ϕ LP is continuous (and even Lipschitz), according to (H ϕ LP ). We also refer to the comments of modeling issues at the end of Section 2.2.2. According to (2.4), f F V has therefore the same regularity as the numerical flow ρ → F(ρ), defined by

F j,j+1 i (ρ) = f j i + f j+1 i 2 - max{|df j i |, |df j+1 i |} 2 ρ j+1 i -ρ j i , with f j i = ρ j i (1 -ρ j i
) and df j i = 1 -2ρ j i , which is obviously continuous.

Step 2: the sequence (ρ n ) n∈N is uniformly bounded in W 1,∞ (0, T ). For u ∈ U ad given, let us introduce the function L : t → 1 2 ∥ρ(t)∥ 2 R Nc×Nr . Up to standard renormalizations, we will suppose in this proof that v max i = ρ max i = 1 and that ∆x i does not depend on the index i. According to the chain rule, L is differentiable, and its derivative reads 2

L ′ (t) = ⟨ρ ′ (t), ρ(t)⟩ R Nc×Nr = ⟨g(ρ(t), u(t)), ρ(t)⟩ R Nc×Nr = - 1 ∆x Nr i=1 Nc j=1 F j,j+1 i (ρ) -F j-1,j i (ρ) ρ j i = - 1 ∆x Nr i=1   γ R i ρ Nc i + Nc-1 j=1 F j,j+1 i (ρ)ρ j i -γ L i ρ 1 i - Nc j=2 F j-1,j i (ρ)ρ j i   = 1 ∆x Nr i=1   γ L i ρ 1 i -γ R i ρ Nc i + Nc-1 j=1 F j,j+1 i (ρ) ρ j+1 i -ρ j i   .
For a given i ∈ 1, N r , and j ∈ 1, N c -1 , we have

F j,j+1 i (ρ) ρ j+1 i -ρ j i = f j i + f j+1 i 2 - max{|df j i |, |df j+1 i |} 2 ρ j+1 i -ρ j i ρ j+1 i -ρ j i ≤ f i (σ i ) + (ρ j+1 i -ρ j i ) ρ j+1 i -ρ j i , since f j i ≤ f i (σ i ) = 1/4
according to the explicit expression (2.1) of f i , and

|df j i | ≤ v max i = 1.
Since the local Lax-Friedrich numerical scheme is monotone, we infer that ρ j i is lower than ρ max = 1. We then have

F j,j+1 i (ρ) ρ j+1 i -ρ j i ≤ 2f i (σ i ) + ρ j+1 i -ρ j i 2 = 1 2 + ρ j+1 i -ρ j i 2 ≤ 1 2 + (ρ j i ) 2 + (ρ j+1 i ) 2 -2ρ j i ρ j+1 i ≤ 1 2 + (ρ j i ) 2 + (ρ j+1 i ) 2 ,
2 we drop the time dependancy notation for readability.

Optimal scenario for road evacuation in an urban environment Strasbourg University by positivity of the ρ j i . Using the majoration

γ L i ρ 1 i -γ R i ρ Nc i ≤ f i (σ i ) × 2 = 1/2, all the calculations above yield L ′ (t) ≤ 1 ∆x Nr i=1   1 2 + Nc-1 j=1 1 2 + (ρ j i ) 2 + (ρ j+1 i ) 2   ≤ 1 ∆x N r N c 2 + ||ρ(t)|| 2 .
We infer the existence of two positive numbers ᾱ, β that do not depend on u, such that L ′ (t) ≤ ᾱL(t) + β for a.e. t ∈ (0, T ). By using a Grönwall-type inequality, we get

L(t) ≤ L(0)e ᾱt + t 0 βe ᾱ(t-s) ds = L(0) + β ᾱ (1 -e -ᾱt ) e ᾱt , (2.8)
Step 3: conclusion. According to (2.8), the sequence

(ρ n ) n∈N is bounded in L ∞ (0, T ). Since ρ ′ n = g(ρ n , u n ) a.
e. in (0, T ) and g is continuous, it follows that (ρ n ) n∈N is uniformly bounded in W 1,∞ (0, T ). Therefore, by using the Arzela-Ascoli theorem, this sequence converges up to a subsequence in C 0 ([0, T ]) toward an element

ρ * ∈ W 1,∞ (0, T ).
Let us recast System (LWR-sd) as a fixed point equation, as

for every t ∈ (0, T ), ρ n (t) -ρ 0 = t 0 g(ρ n (s), u n (s)) ds,
we obtain by letting n go to +∞,

ρ * (t) -ρ 0 = t 0 g(ρ * (s), u * (s)) ds,
meaning that ρ * solves System (LWR-sd) with u * as control. Finally, according to the aforementioned convergences and since the functionals S and B are convex, it is standard that one has

lim n→+∞ C T (u n ) = C T (u * ), lim inf n→+∞ S(u n ) ≥ S(u * ), lim inf n→+∞ B(u n ) ≥ B(u * ).
We get that J (u * ) ≤ lim inf n→+∞ J (u n ), and we infer that u * solves Problem (P θ ). ■

Optimality conditions

We have established the existence of an optimal control in Theorem 5. In order to derive a numerical solution algorithm, we will now state the necessary optimality conditions on which the algorithm we will build is based.

One of the difficulties is that the functional we use involves non-differentiable quantities. We will therefore first use the notion of subdifferential to write the optimality conditions. In Section 2.4 dedicated to numerical experiments, we will explain how we approximate these quantities.

Let us first compute the differential of the functional C T . To this aim, we introduce the tangent cone to the set U ad .

Definition 4. Let u ∈ U ad . A function h in L ∞ (0, T ) is said to be an admissible perturbation of u in U ad if,
for every sequence of positive real numbers (ε n ) n∈N decreasing to 0, there exists a sequence of functions h n converging to h for the weak-star topology of L ∞ (0, T ) as n → +∞, and such that u + ε n h n ∈ U ad for every n ∈ N.

Proposition 6. Let u ∈ U ad and (ρ, γ) the associated solution to (LWR-sd). We introduce the two matrices M and N defined from the Jacobian matrices of f F V and ϕ LP as

M (ρ, γ, u) = (∂ ρ f F V )(ρ, γ) + (∂ γ f F V )(ρ, γ) (∂ ρ ϕ LP )(ρ, u), (2.9
)

N (ρ, γ, u) = (∂ γ f F V )(ρ, γ)(∂ u ϕ LP )(ρ, u), (2.10)
where we use the notational conventions introduced in Section 2.1.

The functional C T is differentiable in the sense of Gâteaux and its differential reads

dC T (u)h = T 0 (N (ρ, γ, u) ⊤ p) • h dt, (2.11)
for every admissible perturbation h, where p is the so-called adjoint state, defined as the unique solution to the

Cauchy system    p ′ + M (ρ, γ, u) ⊤ p = 0 in (0, T ), p(T ) = c.
(2.12)

Remark 2. In Proposition 6 above, the matrices M and N express respectively the way by which ρ interacts with γ and γ interacts with u.

Proof of Proposition 6. Let u ∈ U ad . The Gâteaux differentiability of C T , U ad ∋ u → ρ and U ad ∋ u → γ is standard, and follows for instance directly of the proof of the Pontryagin Maximum Principle (PMP, see e.g. [START_REF] Lee | Foundations of optimal control theory[END_REF]). Although the expression of the differential of C T could also be obtained by using the PMP, we provide a short proof hereafter to make this article self-contained.

Let h ∈ L ∞ (0, T, U) be an admissible perturbation of u in U ad . One has

dC T (u) h = c • . ρ(T ), (2.13)
where

. ρ = d(u → ρ)h solves the system          ρ′ = (∂ ρ f F V )(ρ, γ) ρ + (∂ γ f F V )(ρ, γ) γ, γ = (∂ ρ ϕ LP )(ρ, u) ρ + (∂ u ϕ LP )(ρ, u)h, ρ(0) = 0.
We infer that ρ solves

   ρ′ = M (ρ, γ, u) ρ + N (ρ, γ, u)h, ρ(0) = 0, (2.14)
with M and N as defined in (2.9). Let us multiply the main equation of (2.12) by ρ in the sense of the inner product, and integrate over (0, T ). We obtain:

T 0 ρ • dp dt dt + T 0 ρ • (M (ρ, γ, u) ⊤ p) dt = 0.
Optimal scenario for road evacuation in an urban environment

Strasbourg University

Similarly, let us multiply the main equation of (2.14) by p in the sense of the inner product, and integrate over (0, T ). We obtain:

T 0 p • d ρ dt dt - T 0 M (ρ, γ, u) ρ • p dt = T 0 N (ρ, γ, u)h • p dt.
Summing the two last equalities above yields

ρ(T ) • p(T ) -ρ(0) • p(0) = T 0 h • (N (ρ, γ, u) ⊤ p) dt.
Using this identity with p(T ) = c and ρ(0) = 0 results in Expression (2.11). ■

From this result, we will now state the optimality conditions for Problem (P θ ). Let us first recall that, according to [36, Proposition I.5.1], the subdifferential of the total variation is well-known, given by

∂ TV(u * ) = {η ∈ C 0 ([0, T ]; R Nr ) | ∥η∥ ∞ ≤ 1 and ηdu * = TV(u * )}.
Let us denote by e i the i-th vector of the canonical basis of R Nr Theorem 7. Let u * = (u * i ) 1≤i≤Nr , denote a solution to Problem (P θ ), (ρ, γ) the associated solution to (LWR-sd). and let i ∈ 1, N r . There exists η * ∈ ∂ TV(u * ) such that

• on {u * i = 0}, one has Ψ • e i ≥ 0, • on {u * i = 1}, one has Ψ • e i ≤ 0, • on {0 < u * i < 1}, one has Ψ • e i = 0,
where the function Ψ : [0, T ] → R Nr is given by

Ψ(t) = N (ρ(t), γ(t), u * (t)) ⊤ p * (t) + θ S Nr i=1 u * i (t) -N max + + θ B η * (t)
and where p * denotes the adjoint state introduced in Proposition 6, associated to the control choice u * .

Remark 3. Written in this way, the first order optimality conditions are difficult to use. In the next section, we will introduce an approximation of the total variation of u * leading to optimality conditions more easily usable within an algorithm.

Proof of Theorem 7. To derive the first order optimality conditions for this problem, it is convenient to introduce the so-called indicator function ι U ad given by

ι U ad (u) = 0 if u ∈ U ad +∞ else.
Observing that the optimization problem we want to deal with can be recast as

min u∈L ∞ ((0,T );R Nr ) J θ (u) + ι U ad (u),
Optimal scenario for road evacuation in an urban environment Strasbourg University it is standard in nonsmooth analysis to write the first order optimality conditions as:

0 ∈ ∂ (J θ (u * ) + ι U ad (u * )) ,
or similarly, by using standard computational rules [START_REF] Ekeland | Convex analysis and variational problems[END_REF],

-∂C T (u * ) -θ S ∂S(u * ) ∈ θ B ∂B(u * ) + ∂ι U ad (u * ), Let u ∈ U ad . The condition above yields the existence of η * = (η * i ) 1≤i≤Nr ∈ ∂ TV(u * ) such that dC T (u * )(u -u * ) + θ S dS(u * )(u -u * ) + θ B Nr i=1 ⟨η * i , u i -u * i ⟩ L 2 (0,T ) ≥ 0.
Since u is arbitrary, we infer that for any admissible perturbation h of u * (see Definition 4), one has

dC T (u * )h + θ S dS(u * )h + θ B Nr i=1 ⟨η * i , h i ⟩ L 2 (0,T ) ≥ 0, or similarly T 0 h •   N (ρ, γ, u) ⊤ p * + θ S Nr i=1 u * i -N max + + θ B η *   dt ≥ 0. (2.15)
To analyze this optimality condition, let us introduce the function Ψ : [0, T ] → R Nr defined by

Ψ(t) = N (ρ, γ, u) ⊤ p * (t) + θ S Nr i=1 u * i -N max + + θ B η * (t). (2.16)
In what follows, we will write the optimality conditions holding for the i-th component of u * , where i ∈ 1, N r is given.

Let us assume that the set I = {0 < u * i < 1} is of positive Lebesgue measure. Let x 0 be a Lebesgue point of u * i in I and let (G n ) n∈N be a sequence of measurable subsets with G n included in I and containing x 0 . Let us consider h = (h j ) 1≤j≤Nr such that h j = 0 for all j ∈ 1, N r \{i} and h i = 1 Gn . Notice that u * ± ηh belongs to U ad whenever η is small enough. According to (2.15), one has

± Gn Ψ(t) • e i dt ≥ 0.
Dividing this inequality by |G n | and letting G n shrink to {x 0 } as n → +∞ shows that one has Ψ(t) • e i = 0, a.e. in I.

Let us now assume that the set I 1 = {u * i = 1} is of positive Lebesgue measure. Then, by mimicking the reasoning above, we consider x 1 , a Lebesgue point of u * i in I 1 , and h = (h j ) 1≤j≤Nr such that h j = 0 for all j ∈ 1, N r \{i} and h i = -1 Gn , where (G n ) n∈N is a sequence of measurable subsets with G n included in I 1

Optimal scenario for road evacuation in an urban environment Strasbourg University and containing x 1 . According to (2.15), one has

- Gn Ψ(t) • e i dt ≥ 0.
As above, we divide this inequality by |G n | and let G n shrink to {x 1 } as n → +∞. We recover that Ψ(t) • e i ≤ 0.

Regarding now the set I 0 = {u * i = 0}, the reasoning is a direct adaptation of the case above, which concludes the proof. ■

Towards a numerical algorithm

In this section, we introduce an exploitable approximation of the problem we are dealing with and describe the algorithm implemented in the numerical part.

An approximate version of Problem (P θ )

The fact that Problem (P θ ) involves the total variation of control makes the problem non-differentiable. Of course, dedicated algorithms exist to take into account such a term in the solution, for instance Chambolle's projection algorithm [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. Nevertheless, in order to avoid too costly numerical approaches, we have chosen to consider a simple differentiable approximation of the term B(u), namely

B ν (u) = Nr i=1 TV ν (u i ), (2.17)
where ν > 0 is a small parameter and the total variation TV(u i ) is approximated by a differentiable functional in L 2 (0, T ), denoted TV ν (u i ), where ν > 0 stands for a regularization parameter. The concrete choice of the differentiable approximation of the TV standard will be discussed in the rest of this section. We will also give some elements on its practical implementation.

We are thus led to consider the following approximate version of Problem (P θ ):

inf u∈U ad J θ,ν (u) . (P θ,ν )
where U ad is given by (2.5), and J θ,ν is given by

J θ,ν (u) = C T (u) + θ S S(u) + θ B B ν (u), (2.18)
We will see that this approximation is in fact well adapted to a practical use. Indeed, the first order optimality conditions for this approximated problem can be rewritten in a very concise and workable way, unlike the result stated in Theorem 7. This is the purpose of the following result.

Theorem 8. Let u * ∈ U ad denote a local minimizer for Problem (P θ,ν ) and (ρ, γ) the associated solution to (LWR-sd). Then, u * satisfies the first order necessary condition

Λ(•) = 0 a.e. on [0, T ],
Optimal scenario for road evacuation in an urban environment

Strasbourg University

where Λ : [0, T ] → R Nr is defined by

Λ(t) = min {u * (t), max {u * (t) -1, ∇ u J θ,ν (u * )(t)}} ,
and

∇ u J θ,ν (u * ) : [0, T ] ∋ t → N (ρ(t), γ(t), u(t)) ⊤ p * (t) + θ S Nr i=1 u * i (t) -N max + + θ B Nr i=1 ∇ u TV ν (u i )(t),
where p * has been introduced in Theorem 7, the min, max operations being understood componentwise, and the term ∇ u denoting the gradient with respect to u in L 2 (0, T ).

Proof. The proof is similar to the proof of Theorem 7. Indeed, let u * be a local minimizer for Problem (P θ,ν ).

The first order optimality conditions are given by the so-called Euler inequation and read dJ θ,ν (u * )h ≥ 0, or similarly

T 0 ∇ u J θ,ν (u * ) • h dt ≥ 0.
for every admissible perturbation h, as defined in Definition 4. Let us fix i ∈ 1, N r . By mimicking the reasoning involving Lebesgue points in the proof of Theorem 7, we get

• on {u * i = 0}, one has ∇ u J θ,ν (u * ) • e i ≥ 0, • on {u * i = 1}, one has ∇ u J θ,ν (u * ) • e i ≤ 0, • on {0 < u * i < 1}, one has ∇ u J θ,ν (u * ) • e i = 0.
Note that, on {u * i = 0}, one has ∇ u J θ,ν (u * ) • e i ≥ 0 and then Λ(t)

• e i = min{0, max{-1, ∇ u J θ,ν (u * )(t) • e i }} = 0. On {u * i = 1}, one has ∇ u J θ,ν (u * ) • e i ≤ 0 and therefore Λ(t) • e i = min{1, max{0, ∇ u J θ,ν (u * )(t) • e i }} = 0. On {0 < u * i < 1}, one has ∇ u J θ,ν (u * ) • e i = 0 and therefore Λ(t) • e i = min{u * i , max{u * i (t) - 1, 0}} = 0.
Conversely, let us assume that Λ(•) = 0. On {u * i = 0}, one has .

0 = Λ(t) • e i = min{0, max{-1, ∇ u J θ,ν (u * )(t) • e i }}, so that max{-1, ∇ u J θ,ν (u * )(t) • e i } ≥ 0 and finally, ∇ u J θ,ν (u * )(t) • e i ≥ 0
Hence, the term u n i corresponds to the control at road i and time n. The simplest discrete version of the TV semi-norm is given by

TV(u) ≃ Nr i=1 N T n=2 |u n i -u n-1 i |.
Optimal scenario for road evacuation in an urban environment

Strasbourg University

In order to manipulate differentiable expressions, we introduce for some ν > 0, the smoothed discrete TV operator is given by

TV ν : R Nr×N T ∋ (u n i ) 1≤i≤Nr 1≤n≤N T → Nr i=1 N T n=2 a ν (u n i -u n-1 i ), (2.19)
where

a ν : R ∋ x → √ x 2 + ν 2 . Let (i 0 , n 0 ) ∈ 1, N r × 1, N T .
In what follows, we will use a discrete equivalent of Theorem 8, involving the gradient of TV ν , obtained from the expression

∂ u n 0 i 0 TV ν (u) =      a ′ ν (u n 0 i 0 -u n 0 -1 i 0 ) -a ′ ν (u n 0 +1 i 0 -u n 0 i 0 ) if 2 ≤ n 0 ≤ N T -1 a ′ ν (u N T i 0 -u N T -1 i 0 ) if n 0 = N T -a ′ ν (u 2 i 0 -u 1 i 0 ) if n 0 = 1, where u = (u n i ) 1≤i≤Nr 1≤n≤N T is given

Numerical solving of the primal and dual problems

The models we use are already discretized in space. We now explain how we discretize them in time. We recall that, at the end of section 2.4.1, we have already considered that the controls are assimilated to piecewise constant functions on each cell of the considered mesh and on each time step.

Finite volumes scheme for the primal problem. The main transport equation on network is solved by integrating for each road the finite volume flow given by (LWR-sd) with an explicit Euler scheme:

ρ n+1 i,j = ρ n i,j - △t n △x i,j F n i,j+ 1 2 -F n i,j-1 2 , 1 ≤ i ≤ N r , 1 ≤ j ≤ N c ,
using the local-Lax numerical flux

F n i,j+ 1 2 = F(ρ n i,j , ρ n i,j+1 ), 2 ≤ j ≤ N c -1,
and the Neumann boundary conditions

F n i, 1 2 := γ n i,L , F n i,Nc+ 1 2 := γ n i,R , 1 ≤ i ≤ N r ,
where γ n i,L and γ n i,R are obtained as the solution to the linear programming system

γ n = ϕ LP (ρ n , u n ).
Euler scheme for adjoint problem. To solve the backward ODE (2.12), it is convenient to introduce Z(t) := (M (T -t)) T , and q(t) := p(T -t) such that we are now dealing with the Cauchy system:

   q ′ (t) = Z(t) q(t), t ∈ (0, T ), q(0) = c,
Optimal scenario for road evacuation in an urban environment Strasbourg University that we integrate with a classical explicit Euler scheme:

q 1 = c, q n+1 = (I + △t n Z n ) q n , n = 1, • • • , N T -1.
The solution is finally recovered using that p n = q N T -n+1 .

Optimization algorithms

The starting point of the algorithm we implement is based on a standard primal-dual approach, in which the state and the adjoint are computed in order to deduce the gradient of the considered functional. We combine it with a projection method in order to guarantee the respect of the L ∞ constraints on the control. This method has the advantage of being robust, as it generally allows a significant decrease of the cost functional. On the other hand, it often has the disadvantage of being very local, which results in an important dependence on the initialization.

Moreover, one can expect that there are many local minima, since the targeted problem is intrinsically of infinite dimension, which makes the search difficult.

We will propose a modification of this well-known method, using a fixed point method inspired by the optimality condition stated in Theorem 8.

Projected gradient descent

A direct approach is to consider the gradient algorithm, in which we deal with the condition 0 ≤ u k i ≤ 1 by projection, according to Algorithm 1. The specific difficulty of this approach is to find a suitable descent-step δ k , which must be small enough to ensure descent but large enough for the algorithm to converge in a reasonable number of iterations. We hereby combine this algorithm with a scheduler Equation 2.20 inspired from classical learning rate scheduler in deep-learning [START_REF] Wu | Demystifying Learning Rate Policies for High Accuracy Training of Deep Neural Networks[END_REF] to select an acceptable descent step, where δ 0 and decay are given real numbers.

δ k := δ 0 1 + decay ×k . ( 2

.20)

Algorithm 1 Optimal control by projected gradient descent method (GD)

Require: ρ 0 , u 0 , tol > 0, itermax > 0 Initialization: k = 0 while ||u k+1 -u k || > tol and k ≤ itermax do (ρ, γ) ← solution to      ρ ′ = f FV (ρ, γ) γ = ϕ LP (ρ, u k ) ρ(0) = ρ 0 p ← solution to p ′ + M (ρ, γ, u k ) T p = 0 p(T ) = c u k+1 ← proj [0,1] (u k -δ k ∇J θ,ν (u k )), with δ k > 0 such that J θ,ν (u k+1 ) ≤ J θ,ν (u k ) k ← k + 1 end while
Fixed point method A fixed-point (FP) algorithm is derived using the first-order optimality conditions stated in Theorems 7 and 8, by rewriting them in a fixed-point formulation. We have seen that they write under the form

u i ∈ I µ ⇒ ∂ u i J θ,ν (u) ∈ E µ ,
where µ ∈ {0, * , 1}, I 0 = {0}, I * = (0, 1),

I 1 = {1}, E 0 = R + , E * = {0}, E 1 = R -. This rewrites as u i ∈ I µ ⇒ u i ∈ F µ := {u i |∂ u i J θ,ν (u) ∈ E µ }.
This leads us to compute u i by using the following fixed-point relationship

u k+1 i = 1 × 1 {∂u i J θ,ν (u k )<-κ} + u k i 1 {-κ≤∂u i J θ,ν (u k )≤κ} + 0 × 1 {∂u i J θ,ν (u k )>κ} , (2.21)
where κ > 0 is a small threshold preventing too small gradients from modifying the control. This finally leads us to Algorithm 2.

Remark 4. The advantage of the fixed-point formulation can be illustrated by the fact that some common control configurations can lead to a flat gradient while optimisation is still possible, i.e. gradient descent is too local to detect certain appropriate (and existing) descent direction. Let's construct such a case where the fixed-point algorithm has a major advantage over gradient descent. We recall the constraints of the linear programming for the 2 × 2 case:

0 ≤ γ R 1 ≤ γ R,max 1 , (2.22
)

0 ≤ γ R 2 ≤ γ R,max 2
,
(2.23)

0 ≤ α u γ R 1 + β u γ R 2 ≤ (1 -u 3 )γ L,max 3 , (2.24
)

0 ≤ (1 -α u )γ R 1 + (1 -β u )γ R 2 ≤ (1 -u 4 )γ L,max 4
.

(2.25)

Considering for instance a case where u 3 = 1, then Equation 2.24 yields (2.28) in Appendix for the definition of these functions). It remains that necessarly, γ R 1 = γ R 2 = 0 and the junction is blocked. The point here is that the only way to influence the network -therefore acting on the value of J θ,ν -is by having a control such that γ R 1 or γ R 2 becomes positive, and this can be obtained only when u 4 is set to 0. Indeed, we would now have α u = P ᾱ ε (1) = P β ε (1) = β u = 0, releasing the constraint given by Equation 2.24 on the values of γ R 1 , γ R 2 .

α u γ R 1 + β u γ R 2 = 0. By positivity, this implies that α u γ R 1 = β u γ R 2 = 0. If furthermore u 4 > 0, then it follows that α u = P ᾱ ε (1 -u 4 ) ̸ = 0 and β u = P β ε (1 -u 4 ) ̸ = 0 (see
However, if the algorithm finds a perturbation δu 4 > 0 small enough such that u 4 -δu 4 ̸ = 0, it would yields J θ,ν (u) = J θ,ν (u -δu 4 e 4 ), thus ∂ u 4 J θ,ν (u) = 0. As a result, the gradient descent step is stationnary for u 4 even though the perturbation δu 4 was in the good direction in order to decrease J θ,ν after several more iterations.

With the fixed-point algorithm instead, we find the suitable control u 4 in one iteration from the previous configuration thanks to the indicator functions, since the sign of the gradient contains all the information that we needed. This allows us to circumvent the threshold phenomenon of the gradient described above.

Optimal scenario for road evacuation in an urban environment

Strasbourg University

Algorithm 2 Optimal control with Fixed Point method (FP)

Require: ρ 0 , u 0 , tol > 0, itermax > 0

Initialization: k = 0 while ||u k+1 -u k || > tol and k ≤ itermax do (ρ, γ) ← solution to          ρ ′ = f FV (ρ, γ) γ = ϕ LP (ρ, u k ) ρ(0) = ρ 0 p ← solution to    p ′ + M (ρ, γ, u k ) T p = 0 p(T ) = c u k+1 ← 1 {∇J θ,ν (u k )<-κ} + u k 1 {-κ≤∇J θ,ν (u k )≤κ} + 0 × 1 {∂u i J θ,ν (u k )>κ} k ← k + 1 end while

Hybrid GDFP methods

Since gradient descent theoretically always guarantees a direction of descent (provided we choose a small enough step) but can be slow and/or remain trapped in a local extremum (see Remark 4), and since the fixed point appears more exploratory but does not guarantee descent, it seems worthwhile to investigate the hybridization of both methods.

The chosen algorithm is implemented by computing most of the iterations by gradient descent and using the fixed point method every K ∈ N iterations in the expectation of escaping from possible basins of attraction of local minimizers. Note that K is a hyper-parameter of the method. This is summarized in Algorithm 3.

Instead of performing FP steps regularly, we can choose to space them in order to have more more iterations for the gradient descent to converge. This second version is given in Algorithm 4 where the time between two FP steps increases by a factor of τ .

Algorithm 3 Optimal control with hybrid GD & FP algorithm (GDFP)

Require: ρ 0 , u 0 , tol > 0, itermax > 0

Initialization: k = 0 while ||u k+1 -u k || > tol and k ≤ itermax do (ρ, γ) ← solution to      ρ ′ = f FV (ρ, γ) γ = ϕ LP (ρ, u k ) ρ(0) = ρ 0 p ← solution to p ′ + M (ρ, γ, u k ) T p = 0 p(T ) = c if k%K == 0 then u k+1 ← 1 {∇J θ,ν (u k )<-κ} + u k 1 {-κ≤∇J θ,ν (u k )≤κ} + 0 × 1 {∂u i J θ,ν (u k )>κ} else u k+1 ← proj [0,1] (u k -δ k ∇J θ,ν (u k )), with δ k > 0 such that J θ,ν (u k+1 ) ≤ J θ,ν (u k ) end if k ← k + 1 end while
Optimal scenario for road evacuation in an urban environment

Strasbourg University

Algorithm 4 Optimal control with GDFP with spaced FP steps

Require: ρ 0 , u 0 , tol > 0, itermax > 0, K 0 , τ Initialization: k = 0, K = K 0 while ||u k+1 -u k || > tol and k ≤ itermax do (ρ, γ) ← solution of      ρ ′ = f FV (ρ, γ) γ = ϕ LP (ρ, u k ) ρ(0) = ρ 0 p ← solution of p ′ (t) + M T p = 0 p(T ) = c if k%K == 0 then u k+1 ← 1 {∇J (u k )<-κ} + u k 1 {-κ≤∇J (u k )≤κ} + 0 × 1 {∂u i J θ,ν (u k )>κ} K ← τ K else u k+1 ← proj [0,1] (u k -δ k ∇J (u k )), with δ k > 0 such that J (u k+1 ) ≤ J (u k ) end if k ← k + 1 end while

Numerical Results

This section presents some results obtained in various situations using the methods presented above. In all that follows, we assume that ρ max = 1, v max = 1 and L = 1. We will first validate our approach on single junctions before studying more complex road networks. These are namely: a traffic circle and a three lanes network of intermediate size with a configuration unfavorable to our objective.

Single junctions

We start by considering single junctions of type 1×1, 1×2, 2×1 and 2×2 as they will be the building blocks of the more complex networks.

We consider 50 mesh cells per road and the initial density equals 0.66. The route to empty is always composed of an incoming and an outgoing road. The time interval of the simulations is adjusted in order to allow the route to be completely emptied: the final time T thus equals respectively 6, 3.5, 10 and 5 for the four junctions. We start the optimization algorithms with initial controls equal to 0 and set the convergence threshold to having ∥Λ∥ less than 10 -1 , with a prescribed maximum of 100 iterations. Neither constraints on the number of controls (θ S = 0) nor BV regularization (θ B = 0) are considered for these test cases. Figure 2.1 shows us the comparison between the cost functional history when using the gradient descent (GD), the fixed-point (FP) and the hybrid (GDFP) methods. Specific numerical parameters of the algorithm are given in Table 2.1. We observe fundamental differences in behaviour between the cost functionals obtained by GD and the one obtained by FP. Indeed, GD allows a regular descent whereas FP generates jumps and oscillations, which allows a better exploration of the parameters and avoids certain unsatisfactory local minima. On Figure 2.2, we also plotted the evolution of the optimality conditions Λ for all methods: as expected, we observe that this quantity reaches values close to 0 at the optimal point. The obtained controls are shown in Figure 2.3 and seem relevant. For instance, in the 2×1 case, the control of the outgoing road 3 is fully activated only after time t ≈ 8 to enable the emptying of the incoming road 1 first.

Optimal scenario for road evacuation in an urban environment

Strasbourg University

We further note that the controls are essentially sparse and non-oscillating. They are almost bang-bang, i.e. take only values 0 and 1.

These results seem to suggest that the hybrid method is the most likely to generalize to larger graphs because of its ability to explore and find critical points while still being able to provide convergence. It is therefore the one we will use in the following. Optimal scenario for road evacuation in an urban environment

Symbol

Traffic circle

We consider the traffic circle test case as proposed in [START_REF] Qiu | Runge-kutta discontinuous galerkin method for traffic flow model on networks[END_REF]. It is composed of 8 roads as depicted in Fig. ?? The route to evacuate is chosen to be (1, 2, 3, 5, 8) (see Figure 2.4a) and there is no constraints on the maximal number of controls (θ S = 0) or BV regularization (θ B = 0). The parameters are given in Table 2.2.

including
Results are gathered on Figure 2.5. On the upper left panel, the cost functional history is depicted. The first four iterations of the gradient descent make the cost functional decrease very slowly, except for the second iteration, and then a FP step triggers the jump seen at iteration 5. A few GD steps are then observed and then another FP step results in a second jump at iteration 10. Then the gradient descent is able to reach a satisfactory local minimum by the 12th iteration. We note that the optimality function Λ (Fig. 2.5) follows essentially the same behaviour and reaches a small value at the final iteration.

The control obtained by the algorithm is given in Figure 2.5c, bottom left. We observe first that roads 1 and 6 are allways controlled since they are entry roads on the network and would add new vehicles, and then that road 3 entrance is always controlled to drive the flow to the outgoing road 4. Furthermore, while road 5 entrance is mostly controlled after time t = 8 to let cars from road 3 leave the route before, control on road 7 entrance is

Three lanes network

We consider the three lanes network composed of 23 roads as shown in Fig. In the sequel, we would like to determine an optimal control in order to empty the central lane, made of roads (4, 6, 7, 8, 9) (see Figure 2.6a), with N max = 5. Note that, contrary to the previous test cases, controlling this route may not prevent the flow to circulate from the ingoing to the outgoing road.

First we consider the optimal problem with essentially no limitation on the number of active controls (θ S = 10 -8 ) nor BV regularization (θ B = 10 -8 ). When ||Λ|| 2 becomes lower than 1 we consider that we are close to a basin of attraction, and we give increased effect to regularization by setting θ S = 10 -4 and θ B = 10 -6 .

We use the numerical parameters of Table 2.2: in particular the FP steps are used at iterations 5 and 15. We then observe a slight increasing of the loss J when the regularization coefficients θ S and θ B are increased, mainly due to the high values of S(u) and B ν (u). After the fixed-point step at iteration 15, all the coefficients of the loss are decreasing and reach the convergence tolerance by the 17-th iteration, as shown in Figure 2.7b. A zoom on this phenomena is given Figure 2.8.

The cost C T (u) remains equal to 0.23 whereas the FP step at iteration 15 clearly reduced the staffing constraint S(u) by more than an order of magnitude, at the cost of a slight increase in the total variation. Note that the GD steps were stationnary and that only the FP step led to convergence. The penalization algorithm has clearly solved the two issues: the number of active controls is now less than 7, which is still higher than 5 but far preferable than the earlier 16 in a approached control point of view, and there is no time oscillations of the controls even if it would have been a mathematically correct way to reduce the staffing constraint.

Symbol Name Value

Conclusion

This article has therefore introduced a new network traffic control problem where the control acts not only on the flux constraints at the entrance of the road but also on the statistical behaviour matrices at the junctions.

The theoretical analysis of the optimal control problem led to an original numerical method combining the gradient descent method and a fixed-point strategy. The numerical results obtained on the traffic circle and on the Optimal scenario for road evacuation in an urban environment

where P is a polynomial of degree at most 2, that satisfies P (-1) = α(0, 1) = 1, P (0) = α(u, u) = ᾱ, P (1) = α(1, 0) = 0.

The standard Lagrange interpolation formula yields

P (x) = x(x -1) 2 + ᾱ(1 -x 2 ).
(2.26)

However, due to the division by α in the expression of γ R 1 , the case α(u 2 = 1, u 3 = 0) = 1 -u 2 = 0 is degenerate. To overcome this definition problem, we introduce a small perturbation parameter ε > 0 and a modification of the polynomial P denoted P ε , such that α ε := P ε (u 2 -u 3 ) and

P ε (-1) = 1 -ε 2 , P ε (0) = α, P ε (1) = ε 2 .
The same reasoning as above yields P ε (x) = P (x) + ε 2 x. Since the interpolation procedure does not guarantee that the range of the function be contained in [ε 2 , 1 -ε 2 ], we compose the obtained expression with a projection onto the set of admissible values, to get at the end

α ε (u 2 , u 3 ) := proj [ε 2 ,1-ε 2 ] (P ε (u 2 -u 3 )) . (2.27)
Similarly, we also replace the expression (1-u j ) by (1-u j +ε) to avoid the case where (1-u 2 ) = α(u 2 , u 3 ) = 0.

We finally obtain the following regularized solution at 1 × 2 junctions

γ R 1 = min γ R,max 1 , min (1 -u 2 + ε) γ L,max 2 α ε (u 2 , u 3 ) , (1 -u 3 + ε) γ L,max 3 1 -α ε (u 2 , u 3 ) , γ L 2 = α ε γ R 1 , γ L 3 = (1 -α ε )γ R 1 .
Two ingoing and one outgoing roads (n = 2, m = 1). The ingoing roads are indexed by 1, 2 and the outgoing one is indexed by 3. In any case, we have γ L 3 = γ R 1 + γ R 2 , and γ R 1 , γ R 2 are computed as follows:

• we set γ L,max

3,u := (1 -u 3 ) γ L,max 3 ; • if γ R,max 1 + γ R,max 2 ≤ γ L,max 3,u , then γ R 1 = γ R,max 1 and γ R 2 = γ R,max 2 ;
• else,

-if γ R,max 1 ≥ q γ L,max 3 and γ R,max 2 ≥ (1 -q)γ L,max 3,u , then γ R 1 = q γ L,max 3,u and γ R 2 = (1 -q) γ L,max 3,u . -if γ R,max 1 < qγ L,max 3,u and γ R,max 2 ≥ (1 -q)γ L,max 3,u , then γ R 1 = γ R,max 1 and γ R 2 = γ L,max 3 -γ R,max 1 . -if γ R,max 1 ≥ qγ L,max 3,u and γ R,max 2 < (1 -q)γ L,max 3,u , then γ R 1 = γ L,max 3,u -γ R,max 2 and γ R 2 = γ R,max 2 .
Two ingoing and two outgoing roads (n = m = 2). The ingoing roads are indexed by 1, 2 and the outgoing ones are indexed by 3, 4. We have

γ 3 γ 4 ⊤ = A(u) γ 1 γ 2 ⊤ with A(u) = α(u) β(u) 1 -α(u) 1 -β(u) ,
Optimal scenario for road evacuation in an urban environment

where, given two real number ᾱ, β in (0, 1), the coefficients are defined in the same way as in the case 1 × 2, as

α ε (u) = proj [ε 2 ,1-ε 2 ] P ᾱ ε (u 3 -u 4 ) , β ε (u) = proj [ε 2 ,1-ε 2 ] P β ε (u 3 -u 4 ) .
Here, P ξ ε is defined by

P ξ ε (x) = x(x -1) 2 + ξ(1 -x 2 ) + ε 2 x, ξ ∈ {ᾱ, β}. (2.28)
Let us set γ L,max i,u

:= (1 -u i ) γ L,max i , for i = 3, 4.
To solve the LP 2×2 problem, we first analyze the simplex standing for the polytope of constraints, represented on Figure 2.10: the constraints related on the incoming routes form the rectangle

Ω in := [0, γ R,max 1 ] × [0, γ R,max 2 ]
while those related on the outgoing roads correspond to the regions below the curves C 3 and C 4 , respective graphs of the functions

γ 1 → (1 -u 3 )γ L,max 3 -α u γ 1 β u and γ 1 → (1 -u 4 )γ L,max 4 -(1 -α u )γ 1 1 -β u .
A convexity argument, standard in linear optimization, yields that the exists (at least) a solution lying on a vertex of the set of constraints. To discuss on optimality of each vertex, we first compute Γ u , the intersection of C 3 and C 4 . The point Γ u = (Γ 1u , Γ 2u ) is the intersection point of the lines with cartesian equations

0 0 0 0 γ 1 γ 2 C 3 C 4 Γu
α u γ 1 +β u γ 2 = γ max 3,u and (1 -α u )γ 1 + (1 -β u )γ 2 = γ max 4,u . We obtains    Γ 1u = (1 -β u )γ max 3,u -β u γ max 4,u /△ u , Γ 2u = -(1 -α u )γ max 3,u -α u γ max 4,u /△ u , (2.29)
with if Γ 1u >= 0 and Γ 2u < 0, then γ 2 = 0;

△ u = α u (1 -β u ) -β u (1 -α u ). A
* if α u < β u , then γ 1 = γ max 3,u /α u ; * else, γ 1 = γ max 4,u /(1 -α u ). -if Γ 1u < 0 and Γ 2u ≥ 0, then γ 1 = 0; * if α u < β u , then γ 2 = γ max 3u /β u ; * else, γ 2 = γ max 4u /(1 -β u ). -if Γ 1u < 0 and Γ 2u < 0, then γ 1 = γ 2 = 0; -else, (γ 1 , γ 2 ) = (Γ 1u , Γ 2u ). • Second case: Γ 1u > γ max 1 and Γ 2u > γ max 2 . Then, (γ 1 , γ 2 ) = (γ max 1 , γ max 2
).

• Third case:

Γ 1u > γ max 1 and Γ 2u ≤ γ max 2 ; -if α u < β u , then γ 2 = proj [0,γ max 2 ] γ max 3u -α u γ max 1 β u , γ 1 = γ max 3 /α u if γ 2 = 0 γ max 1 if γ 2 ∈ (0, γ max 2 ];
else,

γ 2 = proj [0,γ max 2 ] γ max 4u -(1 -α u )γ max 1 1 -β u , γ 1 = γ max 4 /(1 -α u ) if γ 2 = 0 γ max 1 if γ 2 ∈ (0, γ max 2
];

• Fourth case: Γ 1u ≤ γ max 1 and Γ 2u > γ max 2 ;

if α u > β u , then

γ 1 = proj [0,γ max 1 ] γ max 3u -β u γ max 1 α u , γ 2 = γ max 3 /β u if γ 1 = 0 γ max 2 if γ 1 ∈ (0, γ max 1 ];
else,

γ 1 = proj [0,γ max 1 ] γ max 4u -(1 -β u )γ max 1 1 -α u , γ 2 = γ max 4 /(1 -β u ) if γ 1 = 0 γ max 2 if γ 1 ∈ (0, γ max 1 ].
Regular approximation of ϕ LP . In the following, in order to use differentiable optimization techniques, we will replace the function ϕ LP , which is only Lipschitz, by an approximation C 1 , replacing all operations consisting in taking the minimum or the maximum of two quantities by regular approximations:

min{x, y} ← x + y -(x + y) 2 + η 2 2 and max{x, y} ← x + y + (x + y) 2 + η 2 2
for a given η > 0.

Optimal scenario for road evacuation in an urban environment

Strasbourg University

Remark 5 (Some remarks on the LP problem). It is notable that In [START_REF] Qiu | Runge-kutta discontinuous galerkin method for traffic flow model on networks[END_REF], a close LP problem has been solved. However, we found that a case corresponding to a particular constraint configuration was forgotten in their analysis. This error could lead to erroneous vertices provifding negative values to some flows. The configuration that was not taken into account in this publication is illustrated on A Julia code for optimal control of a road network

This chapter aims at valorizing and making exploitable the code used for the numerical simulations of the our article presented in chapter 2, that can be found at https://github.com/mickaelbestard/TRoN.jl, as well as retracing some programming notions learned during this thesis. To do so, we will first present the programming interface allowing to reproduce identically the results presented previously, as well as to create new test cases from the existing graphs (front-end: familiarization with the interface and customization of existing function calls). In a second more technical part, we will go into more detail on the design of the interface itself in order to give the interested reader the possibility to tailor the code to his needs, with for example the customization of the cost function, the optimization algorithm or the graph representing the road network. The final part will deal with the limitations and the ways to improve the code in a perspective of scaling up on large problems, like for example on a whole city. (back-end: reappropriation of the code while being aware of some key aspects of the Julia language).

1 Why Julia? Overview of the code An efficient implementation of the algorithms presented above implies two objectives that may be contradictory at first sight. First, we need a language that is sufficiently advanced to allow us to easily build our model, such as the graph structure or the automatic gradient calculations. On the other hand, the complexity of our model requires high numerical performance, especially on larger graphs. While the usual workflow is often to first write a small proof of concept in Python or Matlab before rewriting and optimizing it in a compiled language such as C/C++/Fortran, we chose to develop and optimize (in the numerical sense) in the same language. Julia seems indeed well suited to this purpose.

However, if one of the main ambitions of Julia is to offer the syntactic simplicity of Python with the performance of C, it can be very helpful to understand the internal mechanisms of Julia, and in particular the management of memory allocations, in order to take full advantage of this language and manage to write code that is efficient, generic and easily readable. This is why we have grouped most of the code complexity in a simplified user interface that we will explain in the first sections.

Files organization: We use a common project organization, with a separation of folders into 'sources', 'testscases' and 'unit-tests'. Packages are managed via a 'Pkg' environment which allows to create a local module, which we call TRON (TRaffic On Networks). This also avoids compatibility problems between different possible versions of the libraries, making the project easier to maintain.

The informations about the problem are distributed along different data structures1 :

• Junction, containing the local data of the nodes such as the adjacent routes, static and dynamic distribution matrices, priorities parameters...

• Network, which centralizes all data relating to the network, i.e. the directed graph with local metadata such as maximum local speed and density, the set of junctions, the LWR model flow and its discretization to a numerical flow, ...

• Mesh, the spatial discretization of a road, i.e. the cells with their area, the lengths of each road2 , ...

• Model, that unifies all previously introduced physical quantities and stores the computations related to the direct problem as well as the jacobians of f F V and ϕ LP .

• The structure Optim connects the road graph with our optimizers in order to perform the optimal control algorithm.

• Finally, Diagnostics is a placeholder that concisely gathers all the data to be saved/plotted. which interact with the following keys external libraries:

• Graph.jl which implements for instance the directed graph structure with useful functions to easily manipulate it (vertices(graph::SimpleDiGraph), edges(graph::SimpleDiGraph), source or destination of an edge, ...)

• ForwardDiff.jl which is used for automatic differentiation in the structure Model.

The arborescence is summarized in the UML diagram Figure 3.1. The convention used is that custom structures are in bold, existing libraries are in italics, and the "+" and "-" signs denote public and private fields, respectively.

Presentation of the interface

This part serves as a general presentation of the code and gives the necessary indications for its use as a standalone tool.

You must have cloned the following git directory https://github.com/mickaelbestard/TRoN.jl and launched the Tron module (cf readme.md).

Pre-configured test cases are located in the folder Tron/example, containing the subfolders:

Strasbourg University

• junctions : 1x1, 1x2, 2x1, 2x2

• circle

• threeways

As an example, let us consider the test case of a junction with one input and one output, located in the folder junctions/1x1/. This folder contains the file Test_1x1.jl, which we import via § ¤ 1 include("Test_1x1.jl")

¦ ¥

In order to trigger the compilation, it is usual to run the code a first time with small parameters. This can be done here by using the default ones with the following function call3 : § ¤

1 test_Optimizer(); ¦ ¥
Once this is done, we launch the proper optimization loop and store the result afterwards: § ¤

1 opt = test_Optimizer(final_time=6, Nc=50, algo="GD") ¦ ¥
The possible values of the parameter algo, corresponding to the different choices of optimization algorithms in chapter 2 are:

• "GD" : projected gradient descent

• "FP" : fixed point algorithm

• "GDFP" : hybridization of GD and FP algorithms Notice that other parameters such as penalized constraints and convergence accuracy can also be prescribed, but the default values are those of the paper and can be left unchanged for now.

We then retrieve the results via § ¤

1 diagnostics(opt, test_case="junction") ¦ ¥
which produces the convergence graphs as well as the optimal snapshots and places them in a sub-folder diagnostics of the current folder. The test_case parameter is optional, and is used to obtain a heatmap instead of a standard 1D graphs, see Figure 3.2. The value "junction" is common to the four test cases of the junction folder, and can be replaced respectively by "circle" or "threeways" in the corresponding test cases. The graph structure is coded via the library Graph.jl which allows, once a graph is initialized, to link node i to node j via add_edge!(graph, i, j).

We then indicate the nodes of the network for which we will have to impose conditions on the edges, i.e. the leaves of the graph. Note that for now this code assumes constant boundary conditions at leaves, it has to be modified to make them time-dependant.

Lines 21 to 24 are used to create and assign the default distribution matrix at each junction. It is given with constant values but will later depend on the value of the control, see chapter 2 for more explanations on the mathematical model.

It is important to note that the size of each route is assumed to be constant equal to 1 throughout the code, and the same is true for the number of cells per route.

Edit optimizer: If we now want to change the way the descent direction is calculated, we need to implement a function with the following signature: § ¤

1 function my_descent(opt::Optim) 2 # some julia code 3 end

¦ ¥

The Optim structure opt::Optim notably contains the following data:

• control at the current iteration Optim.uk = u k and at the previous iteration Optim.ukm = u k-1

• the gradient of the cost function to be minimized, at the current iteration: Optim.dJ = ∇J(u k )

• tolerance threshold for gradient cancellation Optim.κ = κ (in the fixed-point algo, a gradient coordinate is considered as null if its absolute value is below this threshold)

The only constraint that this function must respect is to return a descent direction ∆u k , so that the iteration of the control loop is always written in the form :

u k+1 ← proj(u k -∆u k ),
where proj is the projection of the coordinates between 0 and 1. This translates in particular for the fixed-point algorithm by the following expression:

A Julia code for optimal control of a road network

¦ ¥

where the returned value is

∆u k := u k -F P = u k -(1 ∇J(u k )<-κ + u k 1 |∇J(u k )|≤κ
) so that we have

u k+1 ← proj(u k -∆u k ) = u k -(u k -F P ) = F P,
with all the coordinates of F P being in [0, 1] by construction.

Minimize another functional: Our optimal control software can easily address other problems by modifying the cost functional. We may want to minimize the average speed on some roads, the pollution, etc... The function to be called to calculate the cost function is noted ComputeCost!, and can be very generic under the only conditions that it :

• takes opt::Optim as input

• modifies the attribute opt.J

• returns opt.J, given the current implementation of the linear search algorithm (see line 6): § ¤

Limitations of the code and perspectives

Although this code is somewhat efficient, a huge bottleneck still needs to be adressed for the gradient computation.

Indeed, we chose to be the most generic possible and use automatic differentiation in order to quickly be able to use a new model and/or a new control problem, but derivating through the junctions operators is graph-dependant which makes this task more complex.

While this approach works well on our test-cases, it still relies too much on the hardward capacity of the user with several GiB allocated through the process. The main reason is the way Julia and automatic differentiation work, and more specifically the conflict between generic user code and specialised (and therefore fast) machine code. This is the topic of the next paragraph.

Compiler v/s Interpreter

It is important to understand the concepts of compilation and interpretation. A compiler translates human-readable code, the source code, into machine-readable code (assembler/binary) that can be executed later without having that is language-independent and statically typed. After a series of analyses and optimizations, LLVM generates native machine code that is hardware-optimized and efficient because it runs without further interpretation.

Because each function call is specialized with respect to its signature, it is possible to write generic callers that are then statically typed at runtime. This is a key feature of the language, called Multiple Dispatch. Source code written in a way that allows LLVM to generate this kind of optimized code is said to be Type Stable and ensures high performance. Let us illustrate this with some examples. § ¤

A Julia code for optimal control of a road network Strasbourg University significant. For instance, a number like 0.1 is expressed internally in the following way:

0.1 ⇝ + sign 1.00000000 mantissa × 10 -1 exponent ,
and inaccuracy occurs when an operation is performed between two numbers of very different orders of magnitude, as illustrated Table 3.2, where we write the successive operations used to compute 10 5 + 0.066567896.

Operation Remark 1.00000000 × 10 5 + 0.066567896 = 1.00000000 × 10 5 + 0.00000066 × 10 5 Convertion to the same exponent = 1.00000066 × 10 5

Loss of information Table 3.2: Loss of information using floating-points arithmetic.

The consequence on differentiation is that, for ε equal to the machine precision,

x 0 + ε -x 0 =    0, Float, ε, Maths.
The solution of the forward automatic differentiation implemented in the Julia package ForwardDiff is to separate x 0 and ε along two dimensions, using the so-called dual Numbers.

Dual Numbers. The dual numbers are a two-dimensional commutative unitary associative algebra over the real numbers, arising from the reals by addition of a new element ε with the property ε 2 = 0 (ε is a nilpotent element).

They were introduced by William Clifford in 1873. Such a number writes uniquely as z = a + εb =: (a, b), for some a, b ∈ R. Their use in automatic differentiation is straightforward, because for these numbers Taylor expansions yields:

f (x 0 + εh) = f (x 0 ) + f ′ (x 0 )hε + O(ε 2 ) = f (x 0 ) + f ′ (x 0 )hε = f (x 0 ), f ′ (x 0 )h ,
meaning that we compute both the value and the derivative of the function in one evaluation, in an exact way.

This framework readily extends in higher dimensions into multi-duals:

∀i, j ε i ε j = 0, f (x 0 + ε) = f (x 0 ) + ∇f (a) • ε = f (x 0 ), ∇f (x 0 )
We can again compute ∇f in a single (vector) evaluation instead of computing ∂ i f separately.

Towards type-stable automatic differentiation in Julia, conclusion and perspectives

The use of ForwardDiff has many advantages in terms of implementation (one line only to compute a Jacobian of composed vector functions) but introduces a high execution cost, as we are in a case where type-stability and thus pre-allocations are more difficult to obtain simultaneously. We gain in simplicity of use but the optimization of the code requires a real understanding of some internal mechanisms of Julia.

Chapter 4

Tools for ROMs and ML

Eigenmode decomposition simplifies data by preserving the essential structure and removing redundancies.

This allows complex equations to be transformed into a compact, reliable system by means of a well-chosen projection. However, in the presence of strong non-linearity, conventional methods such as EIM and DEIM run into difficulties. We propose to address this issue using neural networks. After a brief introduction to the latter, the chapter ends with a presentation of stochastic gradient descent. This is an efficient approach to high-dimensional optimization that lies behind methods such as Adam and RMS prop, which are commonly found in standard machine learning libraries.

Reduced Order Modeling

The numerical solution of optimal control problems involves a large number of calls to the solver of the physical system under consideration, as with the direct adjoint looping (DAL) method [START_REF] Ru | The application of the method of quasi-reversibility to the sideways heat equation[END_REF] presented in chapter 2. In order to reduce the computation time associated with these methods, we are interested in the possibility of capturing most of the dynamics with as few degrees of freedom as possible. A common approach to achieving this goal is to use the Reduced Order Model (ROM) framework. This framework is designed to create a compact yet precise model capable of accurately representing a range of parameter variations. The ultimate goal is to enable real-time simulation of parameterized models. This approach is rooted in the recognition that natural phenomena typically exhibit an inherent structure that allows them to be effectively characterized by the superimposition of their most significant patterns. Our approach is based on a priori data from full-order models (FOMs), such as finite elements or finite volumes, which we use to construct an appropriate low-rank subspace onto which we project the dynamics. However, this method quickly runs into difficulties due to significant nonlinearities. To solve this problem, we will develop an adapted closure mechanism to efficiently handle projection errors. Our approach is to exploit the remarkable approximation capabilities of a neural network through a data-driven strategy.

See [START_REF] Ahmed | On closures for reduced order models $-$ A spectrum of first-principle to machine-learned avenues[END_REF][START_REF] Snyder | Reduced Order Model Closures: A Brief Tutorial[END_REF] for a comprehensive review of this topic. We're now going to introduce some of the tools that we will need in the following.

Best low-rank approximation

Any matrix A ∈ M n,m (R) admits a proper orthogonal decomposition given via singular values decomposition (SVD) due to the fact that we can always apply the spectral theorem to the symetric positive definite matrices A T A and AA T . This yields the following theorem:

Theorem 9. (SVD) Let A ∈ M n,m (
R) of rank p ≤ min {n, m}. Then there exists positive real numbers

σ 1 ≥ • • • ≥ σ p > 0, and orthogonal matrices U ∈ O n (R), V ∈ O m (R), such that: A = U ΣV T , Σ = D 0 0 0 , D = diag(σ 1 , • • • , σ p ). (4.1)
Moreover, this theorem is immediately related to our compression concern, since it allows us to approximate any matrix by one of lower rank:

Theorem 10. Eckart-Young (1936) Let A ∈ M n,m (R) of rank p.
Then, for any r ≤ p:

min rg(X)=r ∥A -X∥ F = ∥A -A r ∥ F , (4.2)
where A r = U r Σ r V T r , with U r , Σ r and V T r being the r-th firsts columns of U, Σ and V T respectively.

POD based Reduced Order Modeling (POD-ROM)

Let us consider a system of ODEs in high dimension N x >> 1, obtained for instance from the discretization of a PDE and indexed by a family of parameters µ ∈ R:

dx dt (t; µ) = F (x(t; µ); µ), x(t; µ) ∈ R Nx . (4.3)
In order to reduce the dimension, we make the fundamental hypothesis that the high dimensionnal trajectory x(t; µ) can be represented by a trajectory x(t; µ) that lies on a low-dimensionnal hyperplane, i.e. there exists a vectorial subspace V := span{ϕ j } r j=1 with r << N x , and x ref ∈ R Nx such that:

x(t; µ) ≈ x ref + Φx(t; µ), Φ =    | | ϕ 1 • • • ϕ r | |    ∈ R Nx×r . (4.4)
We will also call Φ a decoder since it maps a low-rank solution x ∈ R r onto the full dimensional space R Nx .

From now on, we will suppose without loss of generality that x ref ≡ 0, for the sake of readability.

Remark 6. Note that the fundamental hypothesis above assumes some linearity in the behavior of the solution, thus will not hold in non-linear cases. This is the topic of subsection 1.3.

We apply the SVD, which is named the proper orthogonal decomposition (POD) in this context, on snapshots of the full model (4.3) to find an orthogonal basis on which to express the solution, such that most of the useful information is stored in the first low-frequency coefficients. The snapshots are data considered to be of high-fidelity since they are computed using the full order model on a fine mesh.

The purpose of the POD is then to determinate a low-rank approximation basis, a decoder, Φ from the datum x(t; µ). To this end, we use the approach of [START_REF] Ito | A reduced-order method for simulation and control of fluid flows[END_REF] that is to form the snapshots matrix:

X := x(t 1 , µ 1 ), • • • , x(t Nt , µ 1 ), • • • , x(t 1 , µ Nµ ), • • • , x(t Nt , µ Nµ ) ,

Tools for ROMs and ML

Strasbourg University [START_REF] Michael | Sparse mri: The application of compressed sensing for rapid mr imaging[END_REF][START_REF] Michael | Compressed sensing mri[END_REF][START_REF] William | Energy preserved sampling for compressed sensing mri[END_REF] and we use it here to sparsely sample the nonlinarity in the ODE. This is achieved thanks to the low-rank structure of the flux that we capture via a second SVD. We first reduce the flux separately, which is called hyper-reduction [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF]. We build the flux snapshots:

N := F (x(t 1 , µ 1 ); µ 1 ), • • • , F (x(t Nt , µ 1 ); µ 1 ), • • • , F (x(t 1 , µ Nµ ); µ Nµ ), • • • , F (x(t Nt , µ Nµ ); µ Nµ ) ,
that we also decompose using SVD:

N = ΞΛW T ,
where Ξ = (ξ i ) and W = (w j ) are orthogonal, and Λ = (λ k ) the matrix of corresponding singular values. The information given by this new modal basis allows us to construct a sparse sampling matrix P ∈ R p×Nx thanks to a greedy algorithm that selects the p sampling locations such as maximizing recursively the projection error between

ξ k and [ξ 1 , • • • , ξ k-1 ], for k ≤ p.
More details about the construction of P and further considerations can be found in the work of [START_REF] Petros | Machine learning for fluid mechanics[END_REF]. The sparse sampling matrix P is then used to construct the oblic projector of rank p << N x :

Π ob := Ξ p P T Ξ p -1 P T ,
such that we now replace the costly resolution of Equation 4.7 by:

dx dt (t; µ) = Φ T Π ob F (Φx(t; µ)). (4.10)
The advantage of this formulation is that the flow evaluation now requires only p << N x evaluation of the nonlinarity, at the coordinates selected by the projector.

Deep learning tools

The last section of this chapter introduces classical considerations in deep learning, which is a tool that is more and more mixed with usual numerical analysis [START_REF] Nathan | Workshop report on basic research needs for scientific machine learning: Core technologies for artificial intelligence[END_REF] and is now extensively used in various fields such as image processing, language analysis or fast PDE resolution. We first introduce the overall context of neural networks with a focus on multilinear perceptrons (MLP). Then, we give some insights about the backbone of the numerical feasibility of those methods, the so-called stochastic gradient descent (SGD).

Neural Networks

A neural network is a parametric non-linear function, that is obtained as the composition of linear weighted functions with nonlinear (activation) functions denoted σ. The classic way of representing this is by stacking layers together:

Input layer z 0 , Hidden layers

z k = σ W k z k-1 + b k , 1 ≤ k ≤ l -1 Output layer z l = W l z l-1 + b l ,
Tools for ROMs and ML

Strasbourg University

where we call θ := {W k , b k } 1≤k≤l the trainable parameters, and activation functions σ(•) are usually choosen amongst tanh, ReLU, softmax, ...

Neural networks provide a very effective tool for approximating functions, and its use can be justified theoretically by numerous universal approximation theorems of which we give here one of the first formulations by Cybenko:

Theorem 11. Universal approximation, [28, Cybenko, 1989]. The set of feed-forward neural networks is dense in the set of functions that are continuous on a given compact, for the uniform topology:

Writing K a compact subset of R n and defining

N N = {x → C • (σ • (W • x + b)); A ∈ M k,n (R), b ∈ R k , C ∈ M m,k (R)}, then N N = C(K; R m ), for ∥•∥ ∞,K .
There are refinements of this theorem which give the same kind of results for less regular functions, such as Bochner-Lebesgue integrable functions, and for bounded networks with a control of the dimension, see [START_REF] Funahashi | On the approximate realization of continuous mappings by neural networks, neural networks[END_REF] for a review. Although other architectures have proved highly effective for certain classes of problems (CNN in image processing, RNN for time series, transformers for language processing, ...), we will restrict ourselves to the class introduced above.

To find the optimal trainable parameters θ = {W k , b k }, we perform a gradient-based optimization of a loss function L that mainly quantifies the data attachment of the parameterized model, as well as regularization criteria. A prototypical example of such methods is explained in the next subsection.

Stochastic Gradient Decent

A central element in the implementation of neural networks is the ability to automatically and accurately compute high-dimensional gradients. Automatic differentiation being the subject of subsection 5.2, we introduce here the essential tool for the practical feasibility of high-dimensional computations: the stochastic gradient descent [START_REF] Robbins | A stochastic approximation method[END_REF].

We begin by discussing its theoretical construction before presenting the algorithms that enable it to be used in practice.

Continuous problem: Let (Ω, A, P) be a probability space, (U, < •, • >) be a Hilbert space whose norm is denoted by ∥•∥, and U ad ⊂ U be a closed convex subset. Let W : Ω → W be a random variable, and j : U × W → R a function that is differentiable in its first variable and integrable in the second:

∀w ∈ W, j(•, w) is differentiable, ∀u ∈ U ad , j(u, •) ∈ L 1 (Ω, A, P) .
The generic problem of interest is the following minimization problem:

(P) : inf u∈U ad J(u), where J : U ad → R u → E W (j(u, W )) . (4.11)
Tools for ROMs and ML

Strasbourg University

Writting µ := P • W -1 the pushforward measure of P by W , we have:

E W (j(u, W )) = Ω j(u, W (ω)) dP(ω) = W j(u, w) dµ(w).
By differentiability of j and derivation under the integral sign, we also have :

∀u ∈ U ad , ∇J(u) = E W (∇ u j(u, W )). (4.12)
Approximated problem: The numerical approximation of the previous problem leads us to consider the following approximate Monte-Carlo problem:

(P app ) : inf u∈U ad J app , where J app : U ad → R u → 1 N N i=1 j(u, w i ) , (4.13)
where w 1 , • • • , w m is the realization of an i.i.d. m-sample of W , i.e. an i.i.d. sequence W 1 , • • • , W m where ∀i, W i has the same distribution as W .

We therefore obtain the following estimate:

∇J(u) ≃ ∇J app (u) = 1 N N i=1 ∇ u j(u, w i ). (4.14)
This formulation is not numerically advantageous, as it involves calculating a large number of gradients at each step, N being in machine learning the size of the dataset. It is indeed conventional [START_REF] Combettes | Proximal Splitting Methods in Signal Processing[END_REF] to consider a projected gradient descent algorithm:

u 0 given, u k+1 ← Π U ad u k -δ k ∇J app (u) , (4.15)
involving the calculation of m gradients, possibly in high dimension, at each iteration k.

To address this problem, we will avoid calculating the expectation along W by randomly picking a sample following the law of W at each iteration, forming a sequence of estimators converging to the solution of the problem:

STOCHASTIC GRADIENT DESCENT ALGORITHM (SGD) 1. Initialization. k = 0 : chose u 0 ∈ U and ρ 0 ∈ R * + .
2. Iteration k. Pick a sample w k+1 from the random variable W and update:

u k+1 = Π U ad u k -δ k ∇ u j(u k , w k+1 ) .
We notice that the stochastic gradient ∇ u j(u k , w) is an unbiased estimator of the deterministic gradient ∇J(u)

according to (4.12).

To improve convergence, we use the stochastic gradient "minibatch", which consists of randomly drawing a "batch", i.e. an n-arrangement in 1, N n , calculating the gradient using the n corresponding realizations of the Chapter 5

Comparative approaches to hyper-reduction using deep learning 1 General framework

In order to use large-scale modeling tools, for example in a control context where the resolution of a large nonlinear system is repeatedly required, it is relevant to consider order reduction. Such methods aim to simulate complex models in real time by appropriately compressing their degrees of freedom. Morally, the approach is close to what is classically done with Fourier series, where the unknown is decomposed on a modal basis so as to preserve the essential dynamics once the high frequencies have been cut off:

ρ(t) ≈ Nm k=0 a k (t)ϕ k .
In our work, the expansion basis will be created from data generated by the model, called snapshots.

For our purpose, we want to solve congestion problems in a large city, involving solving the LWR model several times on a large road network in order to use the optimal control tools introduced in chapter 2. The system can be parameterized with respect to quantities such as maximum velocities (v max i ) i , maximum densities accurate reconstruction of the full state. With only 5 modes, the FV state is somewhat diffuse, but increasing the number of modes allows the FV wavefront to be tracked very well, with a few spurious oscillations due to modal decomposition.

Φ T ρ Φ O O Φ T •F •Φ . . Fθ 0 0 F (ρ), Fθ ( ρ)
The same comparison with the 20 modes is made in Figure 5.9. We see more oscillations in the reduced solution of reference, and a lot of small oscillations for the DEIM state. The NC yields a too small amplitude and the same wavefront as before, whereas the DP's amplitude is truly above the reference but still with a reasonable wavefront location. averaging greatly improves the results. Indeed, the full solution is very well fitted both with respect to amplitude and wavefront locations. Furthermore, the error in time is flattened, suggesting a better stability.

Conclusion and perspectives

In this work, we have introduced a new way to compute a closed hyper-reduction based on a deep learning approach. Two methods are investigated, inspired from the trajectory fitting and the model fitting frameworks.

The numerical results obtained on the single-junction network have shown that the trajectory fitting method seems effective in providing an accurate reduced dynamics for all kind of flows, even the most congested ones.

Introduction

There are many physical cases of flows composed of different gas or liquids interacting together. For example, tissue bodies and tumors can be described as a set of interacting viscoelastic materials. Powder-snow avalanches can be described as a mixture of fluid phases. Similarly, the rheology of the gut microbiota and its interactions with chyme (a mixture of partially digested food and water) and the host can be modeled using mixture theory [START_REF] Labarthe | A mathematical model to investigate the key drivers of the biogeography of the colon microbiota[END_REF]. Complex flows can also be found in many engineering applications involving multiphase systems such as boiling water in nuclear reactors. Therefore, the framework of mixture theory is a common tool to model and study complex flows.

Mathematical models based on mixture theory take the form of systems of partial derivative equations, coupled with algebraic constraints. The theoretical analysis of such systems and the characterization of the qualitative properties of the solutions are extremely complicated [START_REF] Beirão | Diffusion on viscous fluids. existence and asymptotic properties of solutions[END_REF][START_REF] Goudon | On a model for mixture flows: Derivation, dissipation and stability properties[END_REF][START_REF] Kazhikhov | The correctness of boundary-value problems in a diffusion model of an inhomogeneous liquid[END_REF][START_REF] Secchi | On the motion of viscous fluids in the presence of diffusion[END_REF]. Thus, it is important to develop efficient numerical methods able to accurately capture the solutions [START_REF] Calgaro | Modeling and simulation of mixture flows: Application to powdersnow avalanches[END_REF][START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Dutykh | Mathematical modeling of powder-snow avalanche flows[END_REF].

In this article, we are interested in applying mixture models to describe biofilm dynamics. Indeed, mixture theory revealed a powerful approach to represent microbial biofilms where a consortium of cells is embedded in a polymeric structure [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF].

In mixture theory, the unknowns of the model are requested to satisfy certain constraints. As far as the continuous equations are considered, several equivalent formulations of these constraints can be derived and used to bring out the properties of the model. However, the preservation of these constraints by a numerical scheme is a challenge and, once a discretization setup has been adopted, it is not clear that all the formulations of the constraints remain equivalent. This issue can induce a loss of stability and accuracy, and eventually a Numerical schemes for mixture theory models Strasbourg University dramatic loss of key physical properties of the simulated flows. Thus, we adapt and extend the numerical scheme proposed in [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] in order to preserve these constraints. The numerical scheme will be tested and illustrated with a multiphasic model representing the development of a photosynthetic biofilm, with the application for biofuel, protein, or drug production.

The paper is organized as follows. The first section is dedicated to the mixture theory framework with a presentation of the simplified model used to test our numerical scheme. The second section details the numerical scheme and its properties. The third section presents the results and comparison with standard numerical schemes.

2 Mixture theory framework: application to biofilms

Mixture theory framework

The mixture theory framework [START_REF] Rajagopal | Mechanics of mixtures[END_REF], also known as mixture mechanics or continuum mechanics for fluid dynamics, enables describing multi-phasic systems at the mesoscopic scale which is an intermediary scale between microscopic and macroscopic scales. It was introduced in the 1960s by Truesdell [91,[START_REF] Truesdell | The classical field theories[END_REF][START_REF] Truesdell | Rational thermodynamics[END_REF] and generalizes Navier-Stokes equations to multi-phasic systems. The mixture theory framework assumes that each component of the mixture might be present at every point in space and at any time. Moreover, the system's physical properties (ex. viscosity, incompressibility, ...) are naturally included.

Consider a mixture of k components indexed by α ∈ 1, k . Each component is locally described by its volumetric mass density ρ α , its volume fraction ϕ α , and its local velocity v α . The volume fraction represents the relative volume occupied by a component in an elementary normalized piece of volume. Thus, assuming that there is no vacuum they satisfy the algebraic constraint

k α=1 ϕ α = 1. (6.1)
The mixture dynamic depends on mass transfers which are modeled through mass balance equations (6.2a) and the local forces applied to the system which are accounted for through momentum balance equations (6.2b).

Thus, for each component the state variables satisfy the equations:

∂ t (ρ α ϕ α ) + ∇ x • (ρ α ϕ α v α ) = Γ α , (6.2a)
∂ t (ρ α ϕ α v α ) + ∇ x • (ρ α ϕ α v α ⊗ v α ) + ∇ x π α + ϕ α ∇ x P = ∇ x • (ϕ α τ α ) + F α + ϕ α ρ α g + Γ α v α , (6.2b)
where Γ i is the mass exchange term, π α is the elastic tensor, P is the common pressure, τ α the viscous stress tensor, F α the friction forces, and g the gravity force. Depending on the considered application some forces can be neglected and some others might be added.

Depending on the targeted application one can add for each component an extra equation for the evolution of the density ρ α . Nevertheless, liquids are weakly compressible, especially when pressure variations are small. Therefore, in most cases, for liquids the component densities ρ α can be assumed constant. When all the component volumetric mass densities are assumed constant, the mass balance equations (6.2a) are

equivalent to ∂ t (ϕ α ) + ∇ x • (ϕ α v α ) = Γ α /ρ α .
Then summing these equations for each phase leads to the pseudo Numerical schemes for mixture theory models Strasbourg University rate D p . Thus, the mass balance equations for a dissolved component within the phase α writes:

∂ t (ρ α ϕ α θ p ) + ∇ x • (ρ α ϕ α θ p v α ) -∇ x • (ρ α ϕ α D p ∇ x θ p ) = Γ p . (6.9)
where again the source term Γ p represents the mass exchange associated to component p.

Mixture model for biofilm

We focus on a simplified 1D model for biofilms. Biofilms are made of microorganisms A (microalgae, bacteria, or a consortium of both) and an extra-cellular matrix E. The biofilm is usually immersed in water L. Therefore, according to mixture theory framework, see section 2.1, each component α ∈ {A, E, L} is described through three macroscopic variables: the mass density ρ α , the volume fraction ϕ α , and the velocity v α . By definition, the volume fractions satisfy at any time the algebraic volume-filling constraint (6.1) which reads in this case:

ϕ A + ϕ E + ϕ L = 1.
In the one-dimensional case, the mass balance equations (6.2a) writes:

∂ t (ρ α ϕ α ) + ∂ x (ρ α ϕ α v α ) = Γ α , α ∈ {A, E, L}. (6.10)
In this context, the volumetric mass densities ρ α can be assumed to be constant. Thus, the mixture averaged velocity satisfies the pseudo incompressibility constraint (6.11) which writes here:

∂ x ϕ A v A + ϕ E v E + ϕ L v L = Γ A ρ A + Γ E ρ E + Γ L ρ L . (6.11)
For biofilms, there are various biological processes to be taken into account. The main processes are growth, extra-cellular matrix excretion, and death. These reactions are schematically represented in Table 6.1. The parameters η α are pseudo-stoichiometric coefficients that quantify how much a reactant (ex. liquid, algae, substrate, ...) or a product (ex. algae, extra-cellular matrix, ...) is consumed or produced when a reaction occurs.

The functions ψ i are the reaction rates. They describe the speed at which reactions take place as a function of the local composition of the mixture. The source terms read as follows: The growth is mainly induced by substrate (S) assimilation and liquid (L) absorption. However, as a first approximation, we assume that the substrate is in excess. Thus, the growth rate ψ g takes the form ψ g = µ g ρ A ϕ A ϕ L , where µ g is the maximal growth rate. The extra-cellular matrix excretion ψ e and the death rate ψ d are assumed to be proportional to the quantity of microalgae, thus ψ e = µ e ρ A ϕ A and ψ d = µ d ρ A ϕ A respectively.

Γ A = ψ g -ψ e -ψ d , Γ E = ψ e + η E ψ d , Γ L = (1 -η E )ψ d -η L ψ g . Biological reaction representation Name Reactant(s) Rate Product(s) Growth η L L + η S S ψg -----→ A Excretion A ψe -----→ E Death A ψ d -----→ η E E + (1 -η E )L
Numerical schemes for mixture theory models Strasbourg University Nevertheless, biofilms are very complex ecosystems and the biological processes are very simplified here. Thus, a model extension accounting for substrate and oxygen is presented in section 5.3.

In the one-dimensional case and neglecting the gravity, for α ∈ {A, E, L} the momentum balance equations simplify into:

∂ t (ρ α ϕ α v α ) + ∂ x ρ α ϕ α v 2 α + ∂ x π α = -ϕ α ∂ x P + 4 3 ∂ x (µ α ϕ α ∂ x v α ) + F α + Γ α v α .
To keep the model as simple as possible, let us assume that the elastic tensor takes the form of a pressure law, see equation (6.5), for the tissues (ie. algae and extra-cellular matrix). Since the liquid phase is not elastic this term is null for the liquid, namely π L = 0. Similarly, let us assume that the friction forces are constant and symmetric.

Thus, in the expression (6.8) for F α , the term f α,α ′ for (α, α ′ ) ∈ {A, E, L} 2 and α ̸ = α ′ are constant and such

that f α,α ′ = f α ′ ,α .
The model is supplemented by boundary conditions. Let Ω = [0, L] be the domain and ∂Ω its boundary.

In 1D, the domain should correspond to a biofilm core drilling in the orthogonal axis of the support where the biofilm develops. The velocities at the bottom of the domain, which corresponds to the surface on which the biofilm develops, vanish v α (0) = 0, α ∈ {A, E, L}. However, the velocity on the top must satisfy a constraint induced by the incompressibility constraint (6.3). Indeed, the integration over the whole domain of equation ( 6.3) combined with the null velocity at the bottom leads to

ϕ A v A + ϕ E v E + ϕ L v L (x = L) = L 0 Γ A ρ A + Γ E ρ E + Γ L ρ L dx
To enforce this condition, let assume that on the top, the velocities are given by v

x=L = L 0 Γ A ρ A + Γ E ρ E + Γ L ρ L dx.
Remark 10. Although there is no biophysical reason to impose the equality between the top velocities, this assumption remains acceptable in this context. Indeed, our focus concerns the biofilm development and the final time considered prevents the biofilm to reach the top of the domain. Therefore, in our context, the hypothesis that all top velocity are equals should not affect the dynamics of the biofilm growth.

Synthesis of model equations

According to the previous section the PDE system under consideration writes:

ϕ A + ϕ E + ϕ L = 1, (6.12a)
∂ t ϕ A + ∂ x (ϕ A v A ) = Γ A ρ A , (6.12b)
∂ t ϕ E + ∂ x (ϕ E v E ) = Γ E ρ E , (6.12c)
∂ t ϕ L + ∂ x (ϕ L v L ) = Γ L ρ L , (6.12d)
∂ t (ρ A ϕ A v A ) + ∂ x ρ A ϕ A v 2 A + ∂ x π A = -ϕ A ∂ x P + 4 3 ∂ x µ A ϕ A ∂ x v A + F A + Γ A v A , (6.12e)
∂ t (ρ E ϕ E v E ) + ∂ x ρ E ϕ E v 2 E + ∂ x π E = -ϕ E ∂ x P + 4 3 ∂ x µ E ϕ E ∂ x v E + F E + Γ E v E , (6.12f) ∂ t (ρ L ϕ L v L ) + ∂ x ρ L ϕ L v 2 L = -ϕ L ∂ x P + 4 3 ∂ x µ L ϕ L ∂ x v L + F L + Γ L v L , (6.12g)
where the sources terms (Γ α ) α , the elastic tensors (π α ) α and the drag forces (F α ) α are given by:

Γ A = ψ g -ψ e -ψ d , Γ E = ψ e + η E ψ d , Γ L = (1 -η E )ψ d -η L ψ g , (6.13a
)

ψ g = µ g ρ A ϕ A ϕ L , ψ e = µ e ρ A ϕ A , ψ d = µ d ρ A ϕ A , (6.13b
)

π α = γ α ϕ α ϕ ⋆ α βα , α ∈ {A, E}, (6.13c
)

F α = α ′ ̸ =α f α,α ′ (v α ′ -v α ), f α,α ′ = f α ′ ,α α ∈ {A, E, L}. (6.13d)
The system (6.12) is supplemented with the boundary conditions

v α (x = 0) = 0, and v α (x = L) = L 0 Γ A ρ A + Γ E ρ E + Γ L ρ L dx,
for all α ∈ {A, E, L}.

The initial data for the volume faction can be chosen arbitrarily provided they are biologically relevant.

However, to enforce the algebraic constraint on the sum over all the volume fractions (6.1), the velocities have to satisfy the incompressibility constraint (6.11) at all times and therefore the initial velocities must verify this constraint as well. Thus, the initial velocities are defined through a pressure P computed using the incompressibility constraint, see section 3.3, by v 0 α = ṽ0 α -∂xP ρα where ṽ0 α is the initial desired velocity. Here, the system is assumed to be initially at rest so ṽ0 α = 0 for all the phases.

Most of the parameters come from [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] or [START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF]. The viscosity coefficient for microalgae and the extra-cellular-matrix are taken from [START_REF] Peterson | Viscoelasticity of biofilms and their recalcitrance to mechanical and chemical challenges[END_REF]. All the parameter values are gathered in table 6.2.

Numerical scheme

In this section, we are interested in the numerical approximation of the PDE system (6.12). Nevertheless, the general principles and in particular the treatment of the pseudo incompressibility constraint remain valid in a more general context. In such PDE systems, the pressure is defined through the volume filling constraint (6.1), namely ϕ A + ϕ E + ϕ L = 1 for the considered model. The treatment of this constraint and thus the definition of the pressure is always an issue and requires specific treatment. To this end, the momentum equations are treated using a projection correction method inspired by the numerical method introduced by Chorin [START_REF] Chorin | The numerical solution of the navier-stokes equations for an incompressible fluid[END_REF][START_REF] Chorin | Numerical solution of the navier-stokes equations[END_REF][START_REF] Chorin | On the convergence of discrete approximations to the navier-stokes equations[END_REF] and Temam [START_REF] Temam | Sur l'approximation de la solution des équations de navier-stokes par la méthode des pas fractionnaires (i)[END_REF] for incompressible viscous flows. In a nutshell, the momentum equation is decomposed using a time splitting to separate the contribution of the pressure as follows: Friction coefficient between E and L 20 kg m -3 day -1

∂ t (ρ α ϕ α v α ) + ∂ x ρ α ϕ α v 2 α + ∂ x π α = f rac43∂ x µ α ϕ α ∂ x v α + F α + Γ α v α , (6.14a)
∂ t (ρ α ϕ α v α ) + ϕ α ∂ x P = 0. ( 6

Projection correction method

Let us start with the presentation of the time discretization. Let T ∈ R + be the final time and (t n ) n⩾0 a subdivision of [0, T ] such that t n = n k=0 ∆t k . Consider α ∈ {A, E, L} a phase and its associated volume fraction ϕ α and velocity v α . Then, ϕ n α (x) and v n α (x) denote, respectively, their approximation at time t n . To shorten the notations, let us drop the space variable x and denote δt = ∆t n+1 . Assuming that all the quantities are known at time t n , the approximated solution at time t n+1 = t n + δt is computed using the following steps:

1. Update the volume fractions according to the mass balance equations (6.12b)-(6.12d):

ϕ n+1 α = ϕ n α - δt ρ α ∂ x (ϕ n α v n α ) + δt ρ α Γ α .
2. Update the momentum equations without the contribution of the pressure term by solving the following system:

ϕ n+1 α v n+ 1 2 α -ϕ n α v n α = δt ρ α -∂ x ϕ n α (v n α ) 2 -∂ x π n α + 4 3 ∂ x µ α ϕ n+1 α ∂ x v n+ 1 2 α + F α ϕ n+1 , v n+ 1 2 + Γ n α v n α . (6.15)
3. Compute the pressure using the incompressibility constraint (6.3). This step is detailed in subsection 3.3.

4. Update the velocity using the pressure with:

v n+1 α = v n+ 1 2 α - δt ρ α ∂ x P n+1 .

1D space discretization

Following [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF], the space is discretized using staggered grids. This enables avoidance of any odd/even decoupling in the stencil of the discrete version of the system. Moreover, the use of staggered grids also allows to have or deduce naturally the quantity of interest (e. g. deduce the pressure gradient on the velocity mesh grid). Let (x i ) i∈ 0,I be a regular subdivision of the domain Ω such that x i = i∆x with ∆x = L I the mesh step. Let also define the mesh cell centers:

x i+ 1 2 = i + 1 2 ∆x for i ∈ 0, I -1 .
The model variables are located:

• at the mesh cell centers for the volume fraction and the pressure:

ϕ α,i+ 1 2 , P i+ 1 2 for 0 ⩽ i ⩽ I -1
• at the mesh cell edges for the velocities: v α,i for 0 ⩽ i ⩽ I. Model unknowns are discretized using a finite volume scheme. The transport terms in the mass balance equations (6.12b)-(6.12d) are written:

ϕ n+1 α,i+ 1 2 = ϕ n α,i+ 1 2 - δt ρ α ∆x F i+1 (ϕ n α , v n α ) -F i (ϕ n α , v n α ) + δt ρ α Γ α,i+ 1 2 (6.16)
where F i represents the numerical mass flux at the interface x i , which is a function of the neighboring cells.

There are multiple relevant choices for the definition of the numerical flux. For the sake of simplicity, to ensure Numerical schemes for mixture theory models Strasbourg University .

x 0 v 0 x 1 v 1 x 2 v 2 x 3 v 3 x 4 v 4 x 5 v 5 x 1
stability and since it is well adapted to staggered grids, it is convenient to use upwind numerical flux. Thus, the discrete mass flux is defined by

F i (ϕ, v) = F + ϕ i-1 2 , v i + F - i ϕ i+ 1 2 , v i with F + (ϕ, v) =    0 if v ⩽ 0, ϕv if v > 0, and
F -(ϕ, v) =    ϕv if v < 0, 0 if v ⩾ 0.
All the volume fractions are updated using equation (6.16). Thus, the volume-filling constraint enforcement is not guaranteed and depends on the strategy used to compute the pressure, see sections 3.3 and 4.2.

Remark 11. To update the volume fractions and ensure volume-filling constraint enforcement another strategy consists to use equation (6.16) for all the components except one (usually the liquid) which is computed using the algebraic volume-filling constraint (6.1): ϕ α ′ = 1 -α̸ =α ′ ϕ α as done in [START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF].

For the momentum balance equation, following [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF], the transport term is also discretized using an upwind strategy based on the material velocity v, that is the momentum flux is defined by

G i+ 1 2 = v n α,i 2 F + ϕ i-1 2 , v i + F + ϕ i+ 1 2 , v i+1 + v n α,i+1 2
F -ϕ i+ 1 2 , v i + F -ϕ i+ 3 2 , v i+1 .
The other terms of equation (6.15) are discretized using standard approximations. Remark that interpolation on the dual mesh is required only for the zeroth order terms like the momentum supply induced by mass exchanges or friction forces. For these terms, the approximation of the volume fraction on the dual mesh is obtained by approximating the volume fractions using the values in the neighboring cells:

ϕ i = 1 2 ϕ i-1 2 + ϕ i+ 1 2
. Therefore, dropping the α for readability, equation (6.15) is discretized as follows:

ϕ n+1 i v n+ 1 2 i - 4δt 3∆x 2 µ ϕ n+1 i+ 1 2 v n+ 1 2 i+1 -ϕ n+1 i+ 1 2 + ϕ n+1 i-1 2 v n+ 1 2 i + ϕ n+1 i-1 2 v n+ 1 2 i-1 -δtF ϕ n+1 i , v n+ 1 2 i = ϕ n i v n i - δt ∆xρ G i+ 1 2 + π ϕ n i+ 1 2 -G i-1 2 -π ϕ n i-1 2 + δtΓ(ϕ n i )v n i .
Remark 12. In this projection step, the viscosity and the friction are treated implicitly. For the viscosity, this treatment enables the relaxation of the CFL constraint and avoids numerical instabilities.

Remark 13. Like in [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF] the computation of friction forces requires a specific treatment. Indeed, the friction forces depend on the difference between the phase velocities, and when a phase vanishes the velocity can not be deduced from the momentum (ie. ϕv). In the considered applications, areas of pure liquid or biofilm are important so the adaptation of the initial data to avoid phase vanishing is irrelevant. To overcome this difficulty, a strategy consists to treat these terms implicitly so the velocity can be directly computed using the above equation. However, this is costly because it imposes to solve at each time step a linear system of size: number of phases×mesh grid size.

Finally, the space discretization of the correction step is given by: v n+1

Pressure approximation

Let us detail the third step of the projection correction method. This is the key step to enforce the algebraic constraint on the sum over all the volume fractions (6.1). The standard strategy consists in plugging the time discrete version of equation (6.14b):

ϕ n+1 α v n+1 α = ϕ n+1 α v n+ 1 2 α
δt ρα ϕ n+1 α ∂ x P n+1 into the incompressibility constraint (6.11) to obtain the following equation on the pressure:

∂ x α ϕ n+1 α v n+ 1 2 α - δt ρ α ϕ n+1 α ∂ x P n+1 = α Γ α ρ α (6.17)
Thus the pressure can be obtained by solving a non-linear and inhomogeneous Poisson equation. As mentioned above, this strategy relies on the use of the continuous version of the incompressibility constraint. Therefore, there is no guarantee that the algebraic volume-filling constraint will be fulfilled at the discrete level.

To enforce the algebraic volume filling constraint, we adapt the strategy proposed in [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF], which consists in using the fully discretized mass balance equations to deduce the appropriate discrete incompressibility constraint.

To this end, let us assume that the constraint α ϕ n α,i+ 1 2 = 1 is satisfied for all times (t n ) n⩾0 and in all the grid mesh cells. Thus, the sum of the equations (6.16) over the phases leads to Consequently, to ensure that the algebraic volume-filling constraint is met, the pressure must be the solution of the non-linear equation (6.19). The solution can be approximated using Newton's methods. In practice, although this method is more expensive than the standard approach its cost remains reasonable. Indeed, the Jacobian matrix is explicitly known and the solution at the previous time step reveals to be a good initial guess so only very few iterations are necessary to converge. Both strategies are compared in subsection 4.2.

Numerical results

The aim of the paper is to present and test a numerical method able to simulate mixture models for biofilms by guaranteeing the preservation of the algebraic volume filling constraint. Another challenge when one wants to go towards the applications, relies on the difficulty to calibrate the parameters of the model. Many parameters are, up to our knowledge, not available in the current literature and very difficult to extrapolate from experimental data. For example, in [START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF] the elastic tensor (ie. π α ) settings are calibrated so that the biofilm front velocity matches observations. Consequently, any modification of the model requires recalibration. To avoid such difficulties, subsection 4.1 presents numerical simulations based on the numerical scheme presented in section 3, but assumes that the viscosity can be neglected, which enables reusing parameters from [START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF] for the elastic tensors. Secondly, subsection 4.2 presents comparisons between the two strategies to approximate the pressure, still neglecting the viscosity. Finally, subsection 5.1 presents the dynamic of the full model including viscosity and recalibration of the elastic tensors.

Initially, the mixture is only made of microalgae and liquid and the volume fractions are set by ϕ 0 A (x) = max{0, 0.05(x -0.1)(x + 0.1)}, ϕ 0 E = 0, and ϕ 0 L = 1 -ϕ A . (6.20)

As mentioned in subsection 2.3, the system is assumed to be at rest. Thus, the initial velocities are defined by v 0 α = -1 ρα ∂ x P where the pressure P is determined according to the strategy presented in subsection 3.3 to enforce the algebraic constraint on the sum over all the volume fractions (6.1). (6.12) where the viscous terms are neglected, namely µ α = 0 for α ∈ {A, E, L}. The simulation is made using the numerical scheme presented in section 3 and using the strategy based on the adaptation of [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] for the computation of the pressure, see subsection 3.3. In these figures, the left side corresponds to the surface where the biofilm sticks and develops and the right side corresponds to the side covered by the liquid, where nutrients are brought.

Biofilm dynamic without viscosity

According to Figures 6.2a, 6.2b and 6.2c, there is front propagation corresponding to the biofilm (dashed orange curve) development within the liquid. As in [START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF], two areas can be distinguished within the biofilm. For example in Figure 6.2c, on the left side, namely for x ∈ [0, ∼ 6]µm, the biofilm is mainly made of extra-cellular matrix (ie. E), whereas on the biofilm front, namely for x ∈ [∼ 6, 7.6]µm, the biofilm is mainly made of microalgae. On the opposite, the right side, namely for x > 7.6, is made of pure water.

For the biofilm components, the velocities are positive near the front, which is expected and explains the biofilm expansion. Otherwise, the liquid velocity is negative in the biofilm region, which means that the liquid is drained into the biofilm due to its consumption for the biofilm growth.

Volume filling constraint validation

Let us compare the two strategies presented in subsection 3.3 to enforce the algebraic volume-filling constraint, that is computing the pressure P either as the solution of discretization the linear equation (6.17) (standard strategy), or as the solution of the non-linear equation (6.19) (adapted strategy). To this aim, as in subsection 4.1, the system (6.12) without the viscous terms is simulated, but using the standard strategy to enforce the volume The rescaled light intensity is the ratio between the received light and the optimal light intensity I opt , namely:

I(t, x) = I surf I opt exp - L x τ 1 -ϕ L (t, y) dy , (6.22)
where I surf is the light intensity at the surface of the tank (ie. x = L) and τ the attenuation coefficient of the biofilm, assuming that microalgae and extra-cellular matrix have the same attenuation rate. The parameter values associated to the light are gathered in Table 6.2. The global dynamic is comparable to the dynamic observed in Figure 6.2. In particular, the biofilm front position travels at the same speed, and, again, there are two areas within the biofilm: the back which is mainly made of an extra-cellular matrix (for x ∈ [0, 6] in subfigure 6.5c) and the front which is mainly made of microalgae (for x ∈ [6, 7.6]µm in subfigure 6.5c). However, as expected, the volume fraction of biofilm is lower. Indeed, taking into account the effect of light reduces the growth in the shadowed areas and thus the biomass of microalgae.

The extra-cellular matrix is also reduced since it is made from microalgae excretion and dead organisms. The parameter values associated with the inclusion of the dissolved components are gathered in Table 6.2.

Numerically, the transport and reaction terms in the mass balance equations for the solutes are treated similarly to the other components. The diffusion terms are treated implicitly to ensure stability without constraining the CFL condition.

Figure 6.6 represents the time dynamic of mixture components when accounting for the contribution of light and solutes. The global dynamic is comparable to the dynamic observed in Figure 6.5. In particular, the biofilm front position travels at a comparable speed. Again there are two areas within the biofilm: the back which is mainly made of an extra-cellular matrix (for x ∈ [0, 6.5] in subfigure 6.6c), and the front which is mainly made of microalgae (for x ∈ [5.5, 6.6]µm in subfigure 6.6c). However, the biofilm front velocity is slightly slower here. This can be explained by the fact that the lack or excess of solutes in the active part of the biofilm slightly reduces its growth. Indeed, for example, at t = 360h, within the biofilm area, the concentration of substrate is reduced by 12.4% and the concentration of inorganic carbon is reduced by 5.1% relatively to the input values (ie. θ in ). Besides, the concentration of oxygen is increased by 46.3% relatively to θ in,C . These discrepancies are larger at the beginning and tend to decrease over time, see Figure 6.9 in the supplementary material. These results are in good agreement with the results presented in [START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF][START_REF] Polizzi | Understanding photosynthetic biofilm productivity and structure through 2d simulation[END_REF].

Conclusions and perspectives

This article proposes an adaptation of the numerical scheme presented in [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] able to enforce the volume filling constraint in mixture models including mass exchanges. As in [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] the strategy consists in deducing the discrete version of the incompressibility constraint from the discretized mass balance equations. Numerical simulations show that this method enables the enforcement of the total volume filling constraint at the discrete level.

In addition, on the modeling side, previous models from the literature are enriched by the inclusion of viscous terms. These terms are essential to properly model biofilms in their fluidic environment especially when there is a mixing of the surrounding fluid. In this context, this work has allowed us to highlight the importance of designing well-balanced numerical scheme able to efficiently capture the biofilm dynamic when including the viscosity.

Indeed, including the viscosity requires to recalibrate model parameters; in particular the elastic tensors need to be strongly rescaled in order to recover realistic front features. However, with these parameters, the numerical set up is more demanding to reach convergence. This difficulty leads to consider further the design of a specific well-balanced scheme for the problem. To this end, the use of well-balanced numerical schemes able to preserve the equilibrium at the biofilm front can be considered.

Finally, in further works, it would be interesting to include additional biological features. Among others, biofilms are generally multi-species. The framework of mixture theory is well adapted to incorporate different species and such extensions are affordable if the interaction between the species and their metabolisms is known.

To make the model even more realistic and predictive its calibration on experimental data is also particularly interesting. In conclusion, real-life biofilms are 3D and therefore the extension and implementation of the numerical method in 2D and 3D should be considered.

Comment optimiser l'évacuation en cas de crise et l'intervention des services d'urgence sur les réseaux routiers ? Comment simuler rapidement la dynamique du trafic routier avec des modèles réduits ? Cette thèse aborde ces questions en utilisant des approches de contrôle optimal et d'intelligence artificielle. Nous modélisons la première question à l'aide d'un problème de contrôle visant à libérer un axe routier de façon optimale, en utilisant la signalisation et des barrages. Notre méthode utilise une semi-discrétisation de l'équation de Lighthill, Whitham et Richards (LWR) sur un graphe de réseau routier. Nous introduisons des fonctions de commande à chaque jonctions pour traduire une action de délestage de l'axe de circulation. Dans un souci d'implémentation pratique, nous prenons également en compte la contrainte d'obtenir des contrôles parcimonieux. De façon pratique, un tel problème est complexe à mettre en oeuvre numériquement, en raison de sa dimension, du caractère fortement non-linéaire du critère et du nombre important de minima locaux. Nous introduisons et testons un algorithme numérique reposant sur l'analyse des conditions d'optimalité d'optimalité de ce problème. Un code de résolution dans le langage Julia est mis à disposition en open source. Dans une seconde partie, nous cherchons à proposer des modèles réduits de trafic afin de faciliter les simulations numériques de réseaux routiers. Nous cherchons à préserver la fiabilité du modèle tout en réduisant son ordre en utilisant des bases réduites associées à des techniques d'apprentissage profond.

How to optimize crisis evacuation and emergency service response on road networks? How can road traffic dynamics be rapidly simulated using scale models? This thesis addresses these questions using optimal control and artificial intelligence approaches. We model the first question using a control problem aimed at optimally clearing a road axis, using signalling and barricades. Our method uses a semi-discretization of the Lighthill, Whitham and Richards (LWR) equation on a road network graph. We introduce control functions at each junction to reflect a relieving action on the traffic axis. For the sake of practical implementation, we also take into account the constraint of obtaining sparse controls. In practice, such a problem is difficult to implement numerically, due to its size, the highly non-linear nature of the criterion and the large number of local minima. We introduce and test a numerical algorithm based on the analysis of optimality conditions for this problem. A solver code in the Julia language is made available as open source. In a second part, we propose reduced traffic models to ease numerical simulations of road networks. We aim to preserve the reliability of the model while reducing its order by using reduced bases associated with deep learning techniques.

Figure 1 :Chapitre 1 :

 11 Figure 1: Macroscopic view of Paris' road network. https://www.data.gouv.fr/fr/reuses/circulation-a-paris/

Figure 2 ,Chapitre 2 :

 22 Figure 2, est prépondérante dans l'existence de solutions ou la convergence d'un algorithme itératif vers un contrôle minimisant effectivement le critère considéré.

Figure 2 :

 2 Figure 2: Exemples de fonctions classiques bv et non bv.

( a )

 a Road network and route to empty in red. T = 10.0, J(u) = 55.741 0 max (b) Density at final time without control.

Figure 3 :

 3 Figure 3: (Rond-point) Gauche : Configuration du rond-point. Les routes (2, 3, 5, 7) du cercle sont dans le sens inverse des aiguilles d'une montre, les routes (1, 6) sont entrantes et les routes (4, 8) sont sortantes. Le trajet à vider (1, 2, 3, 5, 8) est en rouge. A droite : Simulations numériques sans contrôle au temps T = 10.

  où f F V est le flux volumes finis et γ est le flux numérique aux jonctions. La fonction ϕ LP introduite par Coclite et. al est modifiée dans ce manuscrit pour tenir compte du contrôle. On calcule la répartition optimale des flux aux jonctions à l'aide d'un problème de programmation linéaire sur la somme des traces des flux entrants. À titre d'exemple, voici son expression pour une jonction très simple ayant une route entrante et une route sortante :

Figure 4 :

 4 Figure 4: Global structure of the Julia code for optimal control of traffic flow.

Figure 5 :

 5 Figure 5: Densités de véhicules au temps final pour différentes valeurs de la densité maximale sur la deuxième route. Discrétisation spatiale : N c = 400 cellules par route.

Fθ 0 0 FFigure 6 :

 06 Figure 6: Objectif schématisé de la méthode d'hyper-réduction par apprentissage profond.

Chapitre 6 :

 6 Schémas numériques pour les modèles de théorie des mélanges avec contrainte de remplissage : application aux écosystèmes de biofilms Le dernier chapitre traite de l'évolution de micro-algues dans un milieu aqueux, modélisée par un mélange de trois phases. La complexité du modèle provient de la non-linéarité des termes d'advection, des couplages, ainsi que des phénomènes physiques (gravitation, pression) et biologiques (photosynthèse) pris en compte. La stratégie numérique que nous avons utilisée repose sur l'utilisation de grilles décalées pour améliorer l'évaluation des vitesses aux interfaces des mailles, ainsi que d'un prédicteur-correcteur utilisant le gradient de la pression pour CONTENTS CONTENTS Strasbourg University corriger l'erreur commise sur le calcul des vitesses. Il s'agit d'une prépublication issue d'un projet réalisé lors du Centre d'Eté Mathématique de Recherche Avancée en Calcul Scientifique (CEMRACS) en 2022 avec Léo Meyer et Florent Noisette, sous la direction de Bastien Polizzi, Sébastian Minjeaud, Olivier Bernard et Thierry Goudon. Ce projet s'insère naturellement dans le contexte général de ce manuscrit, car il vise à trouver des solutions numériques pour résoudre le problème de la forte non-linéarité dans un modèle fluide.

Figure 7 :

 7 Figure 7: Vue macroscopique de micro-algues. Azeredo, Joana et al. "Critical review on biofilm methods." Critical reviews in microbiology vol. 43,3 (2017): 313-351.

-----

  Simulation numérique avec Julia, Séminaire INFOMATHS. Sorbonne Université et Université Paris-Dauphine (2022) Optimal scenario for road evacuation in an urban environment, Workshop on Mathematical modeling, Analysis and Approximation of Vehicular and Pedestrian Traffic. Tour (2023) • Conférences Nationales et Internationales : -Simulation numérique avec Julia, Congrès SMAI, La Grande Motte (2021) -CANUM, Evian (2022) -International Conference on Hyperbolic Problems, Malaga (2022) Séminaire EDP, IRMA Strasbourg (2021) Séminaire doctorants, IRMA Strasbourg (2021, 2022, 2023) Rencontres doctorants "PhD Aways Days", Strasbourg, Luxembourg (2022) Groupe de travail Machine Learning, Stochastic gradient descent, IRMA Strasbourg (2022) • Vulgarisation scientifique : -Finaliste Alsace Ma thèse en 180 secondes, Strasbourg (2022) Part I Optimal Control Applied to Road Networks Chapter 1

) with v max i ≥ 0

 0 the maximum speed allowed by the model on the road i and ρ i being bounded by a given ρ max i ≥ 0. See Figure 1.1.

Figure 1 . 1 :

 11 Figure 1.1: Local vehicles speed and flow in the LWR model, with ρ max = v max = 1.

Figure 1 . 3 :

 13 Figure 1.3: Demand and supply fluxes for ρ max = v max = 1.

Figure 1 . 4 :

 14 Figure 1.4: Rarefaction wave and its spatial discretization on a uniform mesh. T.Liard and E. Zuazua, 2023[START_REF] Liard | Analysis and numerical solvability of backward-forward conservation laws[END_REF]

(

  a) A non-continuous function that have bounded variations: the stairs function, often encountered when modeling shockwaves. (b) A continuous function (by continuation) which does not have bounded variations: y(x) = x cos(1/x), with instabilities in the vicinity of zero.

Figure 1 . 5 :

 15 Figure 1.5: Examples of classical bv and non-bv functions.

Figure 1 . 6 :

 16 Figure 1.6: Smoothing of the absolute value with the function a ν .

Figure 2 . 1 :

 21 Figure 2.1: (Single junctions) Cost functional as function of iterations for the gradient descent (GD), the fixed point (FP) and the hybrid (GDFP) methods.

Figure 2 . 2 :

 22 Figure 2.2: (Single junctions) Iterations of the convergence criterion ||Λ||.

Figure 2 . 3 :

 23 Figure 2.3: (Single junctions) Optimal controls obtained with the GDFP method. The routes to empty are respectively (1, 2), (1, 3), (1, 3), (1, 3).

( a )

 a Road network and route to empty in red. T = 10.0, J(u) = 55.741 0 max (b) Density at final time without control.

Figure 2 . 4 :

 24 Figure 2.4: (Traffic circle) Left: Traffic circle configuration. The roads (2, 3, 5, 7) in the circle are counter-clockwise, roads (1, 6) are incoming and roads (4, 8) are outgoing. The route to empty (1, 2, 3, 5, 8) is in red. Right: Numerical simulations without control at time T = 10.

quite the opposite: up to time t = 8 ,

 8 it is activated to let the flow circulate as road 8 entrance is open, then from time t = 8 to t = 10, it is deactivated as outgoing road 8 entrance is now closed. Finally, on Figure 2.5d, we can check that the route is effectively empty at final time. 1 20 10 -6 10 -2 5 2 10 5 × 10 -2 10 -1Three lanes network 0.66 1 5 10 -3 10 -2

  Cost functional history (b) ||Λ||2 history (c) Control (d) Density at final time

Figure 2 . 5 :

 25 Figure 2.5: (Traffic circle). Results obtained with GDFP algorithm with spaced FP steps.

2 . 6 :

 26 all the roads are directed from left to right, with road 1 and road 11 being respectively the incoming and outgoig leaves. Initial density is sill constant equal to 0.66 and a corresponding incoming flux is imposed at the entrance of the network on road 1. Neumann boundary conditions are used at the end of the outgoing road 11. The traffic distribution at each junction are uniform: there are no preferred trajectories. Fig. ?? depicts the result of the simulation without control: the density in the central lane tends to be saturated as several roads lead onto it.

Figure 2 .

 2 Figure2.7a represents the cost functional history during the optimization process. We observe first a rapid decrease of the cost functional during the gradient descent steps, which reaches a basin of attraction at iteration 9 with J = 0.23 and ||Λ|| 2 = 1.001 × 10 -2 . Figure2.7c shows the control provided by the algorithm for this iteration. The time evolution of the road densities with this control is represented in Figure2.9a. The control performs rather well regarding the final cost functional and the final road densities on the central line. However, as might be expected, this remains rather unsatisfactory given the relatively high number of active controls and

Figure 2 . 6 :

 26 Figure 2.6: (Three lanes network) Left: Three lanes network configuration. Road 1 is incoming and road 11 is outgoing and all the others are directed from left to right. The route to empty (4, 6, 7, 8, 9) is in red. Right: Numerical simulations without control at time T = 30.

  Cost functional. (b) ||Λ||2 history.(c) Control (low constraints).(d) Control (increased constraints).

Figure 2 . 7 :

 27 Figure 2.7: (Three lanes network with scheduled constraints) Controls obtained with the GDFP algorithm.

Figure 2 . 8 :

 28 Figure 2.8: Cost functional between iterations 9 and 17.

( a )

 a Mean density (low constraints). T = 30.0, J(u) = 0.234 0 max (b) Final time snapshot (low constraints). (c) Mean density (increased constraints). T = 30.0, J(u) = 0.238 0 max (d) Final time snapshot (increased constraints).

Figure 2 . 9 :

 29 Figure 2.9: (Three lanes network with scheduled constraints) Densities obtained with the GDFP algorithm.

Figure 2 . 10 :

 210 Figure 2.10: Polytope of constraints for the 2x2 LP. The constraints on the input flow are represented by the colored area.

Figure 2 50 Table 2 . 4 : 3 C 4 SFFigure 2 . 11 :

 2502434211 Figure 2.11: Graphical resolution of the 2 × 2 LP problem. polytope of constraint (left) and a zoom (right). F is the erroneous solution and S is the correct one.

( a )

 a Density at initial time on road 1 (local coordinates) (b) Density at initial time on road 2 (local coordinates) (c) Density at initial time on the 1x1 network.

Figure 3 . 2 :

 32 Figure 3.2: Standard 1D graphics and heatmap display of the 1x1 junction. Plotted density: 0.3(sin(5x) + 1), for 0 ≤ x ≤ 2 and ∆x = 1/50.

  to re-read the source code. The interpreter does this translation line-by-line during the execution of the program, which allows greater flexibility since the types can depend on runtime values, but is much slower because less optimisation is possible. There are actually intermediate methods of compilation/interpretation, and Julia uses one of them, called Just-In-Time (JIT) compilation. It is a kind of compromise between the first two methods: the source is dynamically compiled into optimised machine code while the program is running. Specifically, Julia uses LLVM (Low Level Virtual Machine), a compilation platform designed for compiletime and run-time optimization. It first translates the source code into an intermediate representation (LLVM-IR)

Figure 5 . 1 :Figure 5 . 2 :F/

 5152 Figure 5.1: By Slevin48an209 -CC BY-SA 4.0, https://commons.wikimedia.org/w/index.php?curid=95887101

Figure 5 . 3 :

 53 Figure 5.3: Objective of the deep hyper-reduction method.

Figure 5 . 5 :

 55 Figure 5.5: Loss convergence over 500 epochs for the Neural Closure strategy.

Figure 5 .

 5 Figure 5.10 and Figure 5.11, demonstrate that successively adding time iterations (K = 50) and time

Figure 5 . 6 :Figure 5 . 7 :Figure 5 . 8 :Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :

 56575859510511 Figure 5.6: Loss convergence over 2000 epochs for the Differentiable Physics strategy using 5 modes.

Figure 6 .

 6 Figure 6.1 gives an example of the staggered grids with the localization of model variables.

Figure 6 . 1 :

 61 Figure 6.1: Regular staggered grid in one dimension for 5 mesh cells with the volume fractions and the velocities locations. The pressure P and the phase volume fractions (ϕ α ) α are located at the mesh cell centers x i+ 1 2 0⩽i⩽4

  in the standard strategy, an equation on the pressure or its gradient can be deduced by using the time discrete version of equation (6.14b). Since in the correction step, the volume fractions remain unchanged, the time discrete version of equation (6.14b) simplifies into v n+1 α ∂ x P n+1 . Injecting this relation into equation(6.18) gives:

Figure 6 .

 6 Figure 6.2 presents the numerical results for different times of system (6.12) where the viscous terms are

( a )Figure 6 . 4 :

 a64 Figure 6.4: Mixture components volume fractions (first row) and velocities (second row) for different times. In this simulation, the viscosity is included and the elastic tensors for the biological phases are multiplied 9 • 10 7 so the biofilm front velocity matches experimental measurements. The simulation is made using 2048 mesh cells for the space grid.

Figure 6 .

 6 Figure 6.5 represents the time dynamic of mixture components when accounting for the contribution of light.

  Numerical schemes for mixture theory modelsStrasbourg UniversityThe external supply for the dissolved components is modeled through Dirichlet boundary conditions at the top of the bioreactor, namely at x = L. Otherwise, the no flux boundary condition at the bottom of the bioreactor is modeled using the Neumann boundary condition: ∂ x θ p | x=0 = 0 for p ∈ {S, C, O}.

( a )Figure 6 . 6 :

 a66 Figure 6.6: Mixture components volume fractions (first row) and velocities (second row) for different times.

  Dans un souci de reproductibilité des résultats et pour rendre les travaux du précédent chapitre réutilisables, nous mettons le code source à disposition et présentons son utilisation. Ce programme est codé en Julia et peut être trouvé en accès libre en suivant ce lien :

	CONTENTS	CONTENTS	Strasbourg University
		Optim
		+ model:Model
		+ uk:Vector{Vector{Float64}}
		+ J:Real
		+ ?J:Vector{Vector{Float64}} + p_t:Vector{Vector{Float64}	Model
	Diagnostics	+ redlane:Vector{Integer} + diagnostics:Diagnostics +...	+ net:Network + mesh:Mesh + ?_t:Vector{Vector{Real}}
	+ losshistory:Vector{Real} + residualhistory:Vector{Real} + optimalityhistory:Vector{Real} + u_optimal:Vector{Vector{Real}}	+ ComputeDual!() + ComputeCost!() + ComputeGradientCost!() + self() #Optimisation algorithm	+ ?_t:Vector{Vector{Real}} + u_t:Vector{Vector{Real}} + ?f_FV:Matrix{Real} + ?? _LP:Matrix{Real}
	+ ?_optimal:Vector{Vector{Real}} +...		-f_FV!(?, ?):Vector{Real} -? _LP!(?, u):Vector{Real}
			+ ComputePrimal!()
	+ diagnostics(opt:Optim):plot		+ AutoJac!():Matrix{Real}
		Mesh	Network
	ForwardDiff.jl Dual Numbers, Automatic differentiation	+ NumberOfCenters:Integer} + centers:Vector{Real} + nodes:Vector{Real} + dx:Real + Length:Real	+ graph:SimpleDiGraph + Junctions:Vector{Junction} + ?max_t:Vector{Real} + vmax_t:Vector{Real} + adjacencyMatrix:Matrix{Integer}
		+ Mesh(Length, NumberOfCenters, NumberOfNodes)	-flow(?, ?max, vmax):Real -Fnum(?L, ?R, side):Real
			Junction
		SimpleDiGraph.jl
			+ Jmat:Matrix{Real}
		Oriented graph structure	+ priorities:Vector{Real}
		(topological)
			-LP(?in, ?out, u):Vector{Real}

https://github.com/mickaelbestard/TRoN.jl Le chapitre commence par la présentation de la structure générale du code ainsi que l'interface permettant d'interagir avec, que ce soit pour reproduire les résultats ou le réutiliser pour d'autres problèmes faisant intervenir

  . A similar reasoning yields the optimality conditions on {u * i = 1} and {0 < u * i < 1}. The conclusion follows. ■ Practical computation of the TV operator gradient in a discrete framework. From a practical point of view, we discretize the space of controls, which leads us to consider a time discretization denoted (t n ) 1≤n≤N T with N T ∈ N * fixed, as well as piecewise constant controls in time, denoted (u n

i ) 1≤i≤Nr 1≤n≤N T

Table 2 . 1 :

 21 Numerical parameters for single junctions.

Table 2 .

 2

2: Numerical parameters for the traffic circle and the three lanes networks.

Table 2 . 3 :

 23 Numerical parameters for the full model.

3

  Editing an existing test caseLet's now dig a little deeper into the code to modify the test parameters. To do this, let's first unpack the function Function call parameters, cf Table3.1. lines 7-8: A controlled model is initialized from the network. The latter can be chosen among the other pre-configured ones (1x2, circle, ...) or modified via the parameters (given here with their default values) rmax_t=[ 1., 1.], vmax_t=[1., 1.], density_at_leaves_t=[0.66, 0.]. These parameters are vectors whose i-th coordinate applies to the i-th route.

	A Julia code for optimal control of a road network	Strasbourg University
	symbol		meaning	type
	final_time		Final time (dimensionless) for the direct problem	Float
	test_Optimizer. Nc Number of cells per route	Int
	Nmax			Maximum number of simultaneous controls allowed	Int
	§ coef_pen		Penalization coefficient θ S of the constraint i u i (t) ≤ N max	Float	¤
	1 function test_Optimizer(;final_time=0.3, Nc=5, Nmax=2, coef_pen=0, coef_reg=0, 2 itermax=3, CV_tol=1e-2, algo="GD", lane=Bool.([1,1]), u_init=0.0, coef_reg Penalization coefficient θ B of the BV regularization Float
	3 itermax	initstate=Constant(value=0.66), fp_tol=0.0, fp_trigger=10, delay_factor=1.0, Maximum number of iterations for the optimization algorithm Int
	lane	4	ismaskcontrol=false, network = graph_1x1(1.0), lr_init=1e-1, decay=1e-1 "One-hot vector" indicating the routes to be evacuated Vector{Bool}
	u_init	5 6	)	Control initial value	Float
	7 initialstate	model = Model(network, Nc, Tf=final_time) Initial condition of the direct problem	InitFunction (custom)
	fp_tol	8 9	init(model) Tolerance on the sign of the gradient for the fixed point algorithm Float
	10 fp_trigger	optim = Optim(model, Nmax, coef_pen, coef_reg) FP triggering delay in GDFP	Int
	11 12 delay_factor	set_scheduler(optim, lr_init=lr_init, decay=decay) set_redlane!(optim, lane) spacing of fp_trigger in GDFP	Float
	13 CV_tol	optim.fp_tol = fp_tol Convergence tolerance for the optimization algo	Float
	14 ismaskcontrol Indicates if we restraint some roads to be controlled optim.trigger = fp_trigger	Bool
	15 16 network	optim.delay_factor = delay_factor optim.ismaskcontrol = ismaskcontrol Road network to optimize	Network (custom)
	lr_init	17	[u .= u_init for u in optim.ukm] Initial learning rate (can be scheduled)	Float
	decay	18 19	[u .= u_init for u in optim.uk] [u .= u_init for u in optim.ukp] Decay in the lr scheduler	Float
		20
		21	optim(itermax=itermax, CV_tol=CV_tol, algo=algo, initstate=initstate)
		22
		23	return optim
	24 end ¦				¥
	Line by line:
	lines 1-4: line 10: Creation of a structure containing the data of the control problem
	lines 11-19: Assigning parameters in the optimization structure

line 21: The structure is also callable as a function. Launching the optimization loop.

The structure optim is returned to the user and contains all the data of the problem.

Table 3 . 1 :

 31 Call parameters of the function test_Optimizer.Lines 3 to 7 are used to create the structure Network which contains the macroscopic LWR model as well as what could be called a topological directed graph, the nodes having no coordinates.

	function graph1x1(;rmax_t=[1.,1.], vmax_t=[1.,1.], density_at_leaves_t=[0.66,0])
	Nr = 2
	Nv = 3
	Nl = 2
	Net = Network(Nr, Nv, Nl, rmax_t, vmax_t)
	# graph topology
	add_edge!(Net.graph, 1,2)
	add_edge!(Net.graph, 2,3)
	# placement of "leafs" nodes and corresponding density values
	# Graph is: 1 (enter node) --> 2 (junction) --> 3 (exit node)
	set_list_leaves(Net, [1,3])
	set_ρ_leaves_t(Net, density_at_leaves_t)
	init_Network(Net)

4 Customize the code

The purpose of this section is to show what to modify in the code to solve other similar problems.

Create a graph: Let's first study how to implement its own graph through the example of the design of the 1x1 graph, whose code is given by § ¤

  As an example, here is the implementation of the cost function used in the article4 : § ¤

	12		end
	13	end
	14	if opt.isregularized
	15		TVU = TotalVariations(u)
	16	end
	17	ρT_sum = dot(opt.redlane, opt.model.ρ_t[end])
	18
	19	opt.J = ρT_sum + θs * Nmax_deviation + θb * TVU
	20 end ¦		¥
		1	### linesearch
		2	iter=0
		3	while J1 >= J0 && iter < maxiter
		4	descent_step * = decay_coef
		5	ukp_guess = proj.(opt.uk .-descent_step * opt.dJ)
		6	J1 = ComputeCost!(opt, u=ukp_guess)
		7	iter += 1
	¦	8	end	¥
	function ComputeCost!(opt::Optim; u=opt.model.u_t,
			initialfunc=Constant(value=0.66))
		θs, θb = opt.regularizations
		set_control!(opt.model, u)
		ComputePrimal!(opt.model, initialfunc)
		Nmax_deviation = 0.
		TVU = 0.

if opt.ispenalized for n in eachindex(u) Nmax_deviation += max(0, sum(u[n]) -opt.Nmax)ˆ2 * opt.model.∆t_t[n]

Table 6 . 1 :

 61 Schematic representation of the biochemical reactions considered in the model.

the graph being oriented, we define "left" and "right" canonically from the direction of each road.

Julia does not use objects, but structures. We present as "methods" functions that call and modify or use the attributes of such structures.

In this version of the code all roads are assumed to have the same lengths, but it is possible to generalize the implementation readily.

The semicolon is used to prevent the returned object from being displayed, and is not necessary as it would be in C, for example.

IMPORTANT: in Julia, the last computed value is automatically returned by the function, hence the absence here of the return instruction.

function unstable(x)

if x<rand()

return "negative x"

end

function stable(x)

if x<0

Table 6.2: Model parameters. The parameters come from [10, 74, 73].

Acknowledgement

• This work was supported by Conseil Regional de Bourgogne Franche-Comté (France).

Optimal scenario for road evacuation in an urban environment Strasbourg University three-lane network have shown that the method is relatively effective in providing relevant controls. It is also possible to impose constraints on the number of active controls. This work can be pursued in several directions. Firstly, it would be interesting to show theoretically that the optimal controls are indeed bang-bang, as suggested by the numerical results. This could provide new improvements of the optimization algorithm. Secondly, the indicator and projection functions make our algorithms very sensitive to threshold or scaling effects. Further analysis could be of great help in determining hyperparameters such as θ or κ generically, since they are set by hand for the time being. Finally, the problem has been tackled here using a discretize (in space) then optimize strategy. Another possibility would have been to first optimize and then discretize. Such strategy may result in another kind of algorithms.

Obviously, it would be interesting to test the algorithm on larger road networks. For this, it might be interesting to couple this method with interesting initialisation methods where an initial control would be obtained with a macroscopic graph. Such an approach has, for example, been developed in a different context in epidemiology [START_REF] Courtès | Reduced modelling and optimal control of epidemiological individual-based models with contact heterogeneity[END_REF].

Appendix

Expression of the function ϕ LP

We provide hereafter the explicit expression of ϕ LP for junctions with at most 2 ingoing and 2 outgoing roads.

We will write (γ R , γ L ) = ϕ LP (ρ, u).

One ingoing and one outgoing roads (n = m = 1). In that case, the solution of the LP problem (2.3) is given by γ

One ingoing and two outgoing roads (n = 1, m = 2). Setting α = α 21 , we obtain

However, this formulation raises a modeling problem. Indeed, this solution does not distinguish between a blocked road (γ R,max j = 0) and a controlled road (1 -u j = 0). Indeed, the first case corresponds to a congestion involving cars on both sides of the intersection, preventing vehicles from entering the second -potentially emptyroad, while in the second case, drivers will not try to enter the controlled road and we revert to the 1 × 1 configuration with the remaining roads.

To capture this difference in driver behavior, we will adjust the parameter α with respect to the control u.

Therefore, α will be considered as a function of (u 2 , u 3 ). More precisely, we introduce a function denoted P of the variable x = u 2 -u 3 , taking the constant value ᾱ ∈ (0, 1) at x = 0 and yielding the relevant 1 × 1 case whenever one road is fully controlled. For the sake of simplicity, we look for a function (u 2 , u 3 ) → P (u 2 -u 3 ),

¦ ¥

This function illustrates a typical case of type instability. In fact, the return type of the function is a number or a string depending on the sign of x-rand(), which is only known at runtime. This prevents the compiler from writing specialized code, resulting in a more complicated intermediate representation (LLVM-IR). We see in the llvm code that some instructions are needed to convert the output to the correct type during execution. § ¤

1 julia> @code_llvm unstable(-1) 2 ; @ REPL [START_REF] Biham | Self-organization and a dynamical transition in traffic-flow models[END_REF]:

¦ ¥

The corresponding LLVM translation is now really short and also fairly human readable: § ¤ 1 julia> @code_llvm stable(-1) 2 ; @ REPL [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF]:1 within `stable3 define i64 @julia_stable_532(i64 signext %0) #0 { 4 top: 5 %1 = icmp sgt i64 %0, 0 6 %. = select i1 %1, i64 %0, i64 0 7 ; @ REPL [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] within `stable8 ret i64 %. 9 }

¦ ¥

Essentially, it means that:

1. The register %1 stores the value of the comparison between the register %0 of type i64 (integer on 64 bits) and 0.

2. Then, the register %. is filled with the value stored in %0 or 0, both of type i64, depending on the value of %1.

3. The last line shows that the returned value is of type i64.

The key point to remember herewith is that it is advantageous to have type-stable code in order to reduce the size of the LLVM code and hence its execution time.

Accurate numerical differentiation within floating-point numbers arithmetic

Before actually talking about automatic differentiation, it is primordial to explain the drawbacks of the floatingpoints arithmetic when it comes to compute accurate local growth rate of a function.

The first approach that comes to mind is, for a function x → f (x), to consider its derivative in some x 0 to be:

for a small ε > 0. A Taylor expansion tells us that this expression yields f ′ (x 0 ) + O(ε) wherever f is C 1 , so we want ε to be as small as possible, especially since errors add up when composing the derivatives. But if this reasonning is relevant in a mathematical framework, it simply causes the derivative to vanish at some point, because there is no guarantee that x 0 + ε and x 0 will be encoded as two differents numbers.

Floating-points arithmetic. The fact that there are infinitely many numbers and only a finite number of transistors implies that not all real numbers can be represented, and that only a certain number of decimals are A Julia code for optimal control of a road network

Strasbourg University

Although the current version of the code is satisfactory from a research point of view (a few minutes to run test cases of one or two dozens of roads), it is still very insufficient if we envisage a scale-up to an industrial environment, with test cases of an entire city with several hundred or even thousands of roads to simulate.

After profiling the code, it is clear that most of the computation time is spent computing the gradient of the linear programming problem, which is not surprising since the graph structure is used extensively.

Beyond the necessary considerations of parallelization, it seems that large gains in execution time and memory are achievable in a first step by reaching a suitable mix between static pre-allocations and dynamic allocations. In the first case the initial launch cost is increased for the benefit of a faster execution, contrary to the second case where the initial cost is lower at the cost of an overload of the execution time. In this case, it is not guaranteed that the hardware supports too many simultaneous allocations, but a too slow execution of the optimization algorithm is also undesirable, more specifically in a crisis management context.

The natural perspective that emerges from this is to achieve real-time optimal control over large networks, such as thousands of edges for a city like Paris. Amongst all tools that can be used to achieve this, graph reduction and reduced order models seems to be of great interest and the later is investigated in chapter 5 which is devoted to the hyper-reduction of the finite volumes scheme on a network.

Part II Numerical Schemes for Strongly Nonlinear Fluid Dynamics

Tools for ROMs and ML Strasbourg University which is of size N x × N µ N t . We then look for any K ≤ r for a decoder Φ K that solves:

By Theorem 10, the solution is given by the K firsts columns of the eigen vectors of the SVD decomposition:

Remark 7. The numbers σ i , called the singular values of X , are equal to the eigenvalues of both X T X and X X T , and are ordered in decreasing order:

To make an analogy with Fourier series, note that the eigenvectors associated with the first eigenvalues are called low-frequency modes. Ignoring the information provided by the smallest eigenvalues is then similar to using a low-pass filter.

With this construction we now are able to project (4.3) on the hyperplane generated by the snapshots. Writing x := Φ T x, we now consider:

that is approximated, using our fundamental hypothesis, by:

In the case where F is linear, we are then led to solve:

We perform a single off-line calculation of the full-dimensional matrix, denoted as L := Φ T F Φ. Subsequent numerical integration exclusively operates within the reduced dimension. During the online phase, for example using explicit Euler integration on a regular grid (t n ) 1≤n≤Nt and writing xn µ := x(t n ; µ), we compute:

A nonlinear function F leads to the obligation to compute the high dimensional expression F (Φx n µ ) at each time step. This is obviously not efficient in a reduction framework. An approach to deal with this case is presented in the next section.

Nonlinear case: DEIM

The nonlinarity of the flux makes the reduced basis from snapshots inneficient to compress the trajectory of the system because the evaluation of the nonlinearity costs as mush as the resolution of the high fidelity model. The construction of an actual low rank evaluation of the non linear flux that we present here is inspired by what is often known as sparse sampling [START_REF] Candes | Stable signal recovery from incomplete and inaccurate measurements[END_REF] and relies on the statement that meaningful data are sparse in the set of all possibles data, i.e. natural data are compressible. From that, one can infer that a suitable configuration of sparse interpolations points (sometimes refered to as sensors) allows an accurate reconstruction of the signal without needing the entier sampling of the high fidelity data of reference. This principle has a lot of applications Tools for ROMs and ML Strasbourg University w i , and repeating until all the training data have been run through. Then start again.

2. Iteration k.

• We draw w k+1 from W , which we partition into N/n batches (b t ) t of size n.

• batches iterations:

Comparative approaches to hyper-reduction using deep learning Strasbourg University (ρ max i ) i , or initial data (ρ i,0 ) i . For the sake of generality, we will refer to such parameters as µ in all that follows.

For more details on this model, please refer to chapter 2.

After a spatial discretization using N c cells per road, thus writing N x := N r N c the total degrees of freedom of the parameterized numerical unknown ρ(• ; µ) ∈ R Nx on the network, we obtain the following system of ODEs:

(5.1)

Remark 8. We need a fairly large number of cells per route to obtain a high-fidelity model, which will serve as a baseline in our investigations.

The goal now is to build a surrogate model with an approached solution ρ ∈ R Nm associated to a low rank flux F( ρ), with N m << N x , on the form:

We aim to identify an appropriate low-rank subspace that can provide an accurate solution representation.

In line with the POD-ROM framework, our goal is to develop an algorithm that divides the process into offline and online phases. During the offline phase, we create this subspace using prior data, while the online phase exclusively involves low-dimensional computations.

To achieve this, we employ the snapshot approach, wherein we calculate the full-order solution for a set comprising N µ parameters

Writing N s := N µ ×N t the number of snapshots and ρ j,n µ k = ρ j (t n ; µ k ), we have:

to which we apply a Proper Orthogonal Decomposition (POD):

and Σ ∈ M Nx,Ns (R) of rank p such that:

Comparative approaches to hyper-reduction using deep learning

Strasbourg University

The reduced basis with N m modes is given by the first N m columns of the matrix U , i.e.

We use it to apply a Galerkin projection on the high-fidelity system, which yields the following reduced ODE:

that we rewrite, using the approximation from the fundamental hypothesis Φ ρ ≈ ρ:

(5.6)

However, the non-linearity of F raises important questions. Indeed, we can see that evaluating Φ T F(Φ •) has essentially the same cost as F and moreover needs to be done at each time step, thus obliterating any gain in order reduction. A classical way to deal with this problem is to add an hyper-reduction such as the Dynamic Empirical Interpolation Method (DEIM) presented in subsection 1.3:

However, strong nonlinearity in the flux F makes it necessary to use an increasing number of modes to get an accurate dynamic. Therefore, we usually add a closure term, which amounts to consider:

and determine an appropriate C µ that corrects out of the retained modes the error made on the truncated modes.

This point of view is widely discussed in [START_REF] Ahmed | On closures for reduced order models $-$ A spectrum of first-principle to machine-learned avenues[END_REF][START_REF] Snyder | Reduced Order Model Closures: A Brief Tutorial[END_REF].

We highlight the limitations of the DEIM approach for genuine nonlinear flow scenarios in Figure 5.2. In this context, we employ the LWR model on two roads as the governing equation, with specific parameters and model details provided in subsection 3.1.

We chose to simultaneously compute a hyper-reduction flux and a closure term to ensure both reliability and stability of the reduction, all within a model-agnostic framework that derives the flux solely from low-dimensional information, as schematically illustrated Figure 5.3. In the following section, we describe our approach to address this challenge.

Presentation of our strategy

The approach that we propose in this work is to use a neural network to learn the hyper-reduction and the closure in a same function. This leads us to consider the system:

Comparative approaches to hyper-reduction using deep learning

where Fθ is a neural network with weights θ = (θ 1 , • • • , θ p ) that takes ( ρ(t n ; µ), µ) as input and returns a prediction of the compressed flux Fθ ( ρ(t n ; µ), µ) ≈ Φ T F(ρ(t n ; µ); µ). We see that this network is precisely designed to recover the benefits of a reduction, since the online stage uses the reduced dimension only, see Figure 5.3. Note that learning the hyper-reduction only would amount to learn Φ t F(Φ ρ) instead, where ρ is the compression of the full order trajectories.

Remark 9. To simplify the notations in what follows, we will always consider a constant time step ∆t such that, for each integer n, we write t n = (n -1)∆t. Moreover, we write ρ n µ := ρ(t n ; µ) for all (n, µ).

In everything that follows, we consider an Euler explicit time integration of Equation 5.9:

i.e. the neural network's objective is to calculate a time integrator for a single iteration.

We will compare two approaches that produce the same type of network, but differ in their training strategies.

In either case, the network will be a multi-linear perceptron (MLP), with hyper-parameters described in detail below.

The dataset is assembled as follows:

• Snapshot construction: For each µ ∈ P µ we solve the high-fidelity Equation 5.1 from which we save the pairs ρ µ , F (ρ µ ; µ) ; µ ∈ P µ .

(5.11)

• Construction of the reduced space: An SVD is performed on the snapshots to form and store the encoding matrix Φ T as well as the projections of the states and the fluxes

• Neural network I/O: The network will therefore have inputs ρproj µ , µ to predict the targets Fproj µ .

Our neural network is specifically designed to acquire knowledge about an integration flux, which in other words means it grasps the underlying physics of our problem or the model. During training, we employ two well-established methods commonly used in this field: model fitting and trajectory fitting [START_REF] Ahmed | On closures for reduced order models $-$ A spectrum of first-principle to machine-learned avenues[END_REF]. The former involves straightforward computation based on the neural network's output, while the latter is conceptually more intuitive as it remains independent of the model and directly compares solutions.

Neural closure (NC):

Our first approach of the deep-hyper-reduction strategy is to compare directly the output of the neural network to its prediction target, which leads to consider the minimization of the model error:

(5.13)

Comparative approaches to hyper-reduction using deep learning Strasbourg University

Differentiable programming (DP):

In the second approach, we attempt to learn the whole integrator over a given number K ∈ N of time steps. More precisely, writing the explicit time discretization of Equation 5.9 as:

we define by induction the integrator S K θ yielding the reduced state ρn+K µ :

We will hereafter denote the reduced state at time t n+k predicted by the neural network from a low-dimension state ρn µ at time t n by: ρDP,n+k µ,θ

Our baseline will be the compression of the high-fidelity state ρproj µ := Φ T ρ µ , thus we search for optimal parameters θ opt such that:

We therefore minimize the following loss:

(5.17)

Another loss that we will consider is the time average of the previous one:

This formulation has the advantage of using information not only from ρn+K µ but also from the entire trajectory { ρn+k µ , 0 ≤ k ≤ K}. In this way, we hope to obtain a more generalizable neural network. We observe that for K = 1, this approach boils down to the model learning introduced before, indeed: normalization [START_REF] Santurkar | How Does Batch Normalization Help Optimization?[END_REF]. We will therefore conduct a comparison between model training without this renormalization and trajectory learning specifically for the case when K equals 1. This comparison will help emphasize the impact of time rescaling.

3 Numerical results

Test-case presentation

We want to describe different types of dynamics in the two-road model presented in Figure 1.2, which are influenced by the maximum density allowed on each road and denoted by µ := ρ max ∈ R 2 . We initialize the road network with constant values that are respectively 0.6 on the road going in the junction and 0.1 for the outgoing one. We impose constant density of 0.6 at the network entrance and a free exit condition. Setting Comparative approaches to hyper-reduction using deep learning Architecture and hyperparameters: All approaches are trained on the same fully connected neural network of 4 layers with 32 neurons each and tanh activation functions. We trained the networks with Adam optimizer, batches of size 16 and a 30% dropout rate. The learning rate was scheduled as follows:

iterations 1-10 : lr = 10 -3 , iterations k>10 :

We uniformly sampled 9 values of µ 2 between 0.1 and 0.9 as validation data for both approaches.

Results

We base our comparisons on the DEIM method. Figure 5.2 shows the limitations of this approach for our problem, since the state is not well captured at all with oscillations instead of shock, and no useful information about its height. The DEIM state is essentially equal to the initial condition.

Training: In In contrast, the DP method shows great reconstruction properties since the modes are well learned, allowing

Comparative approaches to hyper-reduction using deep learning

Strasbourg University

This work is still preliminary and can be pursued in several directions. Firstly, it would be interesting to extend the method for more complex networks as those introduced in chapter 2. A suitable interface to the optimal control algorithm would enable the problem to be extended to larger networks, which could lead to real applications using standard computers, and therefore more easily deployable.

Secondly, the training process can be improved with better architecture and/or parameters. It would be interesting to use reccurent neural networks or long-short-term memory, as they are suited for time evolving problems. Investigations on suitable activation functions and normalizations should help improve the overall process.

It would also be interesting to test numerically the conjecture made by Noack et al. [START_REF] Noack | The need for a pressure-term representation in empirical Galerkin models of incompressible shear flows[END_REF] and mentioned in [START_REF] Ahmed | On closures for reduced order models $-$ A spectrum of first-principle to machine-learned avenues[END_REF] that the network trained with model regression (Neural Closure) extrapolates better outside the training interval than the one trained with trajectory regression (Differentiable Physics). However, the latter is claimed to be more accurate in the reconstructive regime, i.e. inside the training interval.

Numerical schemes for mixture theory models Strasbourg University incompressibility constraint:

This means that the local divergence of the averaged mixture velocity is equal to the local volume variation induced by mass exchanges.

The elastic tensor π α can be interpreted as the internal pressure of the component. There are several ways to model this term depending on the nature of the component. When the component α represents particles, as in [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF], there is a close-packing limit. This property can be enforced by using an appropriate expression for π α as

, with γ α > 0, and

where 0 < ϕ ⋆ α < 1 is the so-called close-packing volume fraction limit. When the component α represents softer material like living tissues it can take the form of standard pressure law:

, with γ α > 0, and

where 0 < ϕ ⋆ α < 1 is a threshold, see [START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Preziosi | Multiphase modelling of tumour growth and extracellular matrix interaction: mathematical tools and applications[END_REF]. More complex laws, based on the Flory-Huggins theory:

enable accounting for colligative properties at low concentrations, see [START_REF] Cogan | The role of the biofilm matrix in structural development[END_REF].

The viscous stress tensor τ α is defined by

where the constant µ > 0 stands for the component dynamic viscosity and t ∇v α stands for the transpose of the velocity differential matrix.

The friction force F α is induced by the difference in the relative speed of the mixture components:

with f α,α ′ the friction force law between the components pair α and α ′ . As a first approximation, it can be assumed that f α,α ′ is a strictly positive constant. However, the friction between two components should vanish when one of them disappears. Thus, a more realistic alternative is to consider that friction depends on the local composition and use instead f α,α ′ (ϕ α ϕ α ′ ) r α,α ′ . Nevertheless, the total momentum conservation principle enforces that

Dissolved components, like substrate, can be included. A dissolved component p within a phase α is described through its concentration θ p . In addition to the transport by the phase, it can also diffuse within the phase at a Numerical schemes for mixture theory models filling constraint. The results for the mixture components volume fractions at t = 360h are presented in Figure 6.3a and can be compared to Figure 6.2c. According to these plots, the results are comparable. Similarly, the shape of the pressure gradients curves are also similar, see Figure 6.3b. Nevertheless, according to the dotted purple curve in Figure 6.3b, there is a significant discrepancy close to the biofilm front (ie. at x = 7.6mm)

in the pressure gradients between the two correction strategies. Note that the front is the active part of the biofilm. Namely, it is in this area that the source terms are the largest and lead to significant changes in mixture composition. Therefore, it is expected that the effect of the pressure gradient is observable notably there.

Besides, the pressure can be interpreted as the Lagrange multiplier associated with the volume-filling constraint. Thus, it is important to compare how these strategies enable enforcing at the discrete level the volume filling constraint (6.1). To this end, Figure 6.3c represents the sum of the volume fractions within the domain a time t = 360h. According to this plot, the strategy adapted from [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] enables ensuring the volume filling constraint, whereas the standard strategy does not. Numerically, the maximal error on the volume-filling constraint for the standard strategy is 1.007 • 10 -3 , whereas with the adapted strategy, it is 5.107 • 10 -15 , namely the order of magnitude of the precision used in Newton's method. Moreover, with this adapted method, the error remains negligible throughout the simulation whereas, with the standard strategy, it varies over time, see Figure 6.7 in Appendix 6.

Model extensions

Following insights coming from [START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF], this section presents various relevant extensions of the model and their numerical simulations.

Numerical schemes for mixture theory models (right) at time t = 360h. In figures 6.3b and 6.3c the blue curve represents correspond to a simulation made using the standard strategy for the pressure gradient computation (ie. solving a Poisson equation) and the dashed orange curve represents the results obtained using the strategy adapted from [START_REF] Berthelin | Multifluid flows: A kinetic approach[END_REF] which require the resolution of a non-linear equation (ie. Newthon's method).

Including the viscosity

Adding the viscous terms for the components requires recalibrating the model parameters. Indeed, the viscosity is a measure of the component's resistance to deformation. Therefore, when accounting for the viscosity, the parameters associated with the component's ability to deform must be adapted. In particular, the elastic tensors for the microalgae and the extra-cellular matrix must be recalibrated. Moreover, up to our knowledge, there is no direct measurement of the parameters and they are calibrated, see [START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF], such that the biofilm front velocity matches experimental measurements, see [START_REF] Schnurr | Factors affecting algae biofilm growth and lipid production: A review[END_REF]. However, such calibration is extremely complex because the biofilm front velocity depends also on many other parameters like the growth or death rate. Nevertheless, to get the right order of magnitude of the biofilm front velocity the elastic tensor coefficients must be significantly increased: multiplied by 9 • 10 7 so is set to γ A = γ E = 4.5 • 10 -3 kg/m/day. The global dynamic remains comparable to the dynamic observed in Figure 6.2. In particular, there is still a biofilm traveling front. Again, there are two areas within the biofilm: the back which is mainly made of an extra-cellular matrix (for x ∈ [0, 5] in subfigure 6.4c), and the front which is mainly made of microalgae (for x ∈ [3.97, 6.7]µm in subfigure 6.6c). Nevertheless, a major discrepancy is that the microalgae remain more located at the front when including the viscosity. This is particularly visible at t = 240h when comparing Figure 6.2b and Figure 6.4b. In addition, at t = 240h, the velocities order of magnitude close to the front is larger when including the viscosity. However, the interpretation of this observation is tricky. Indeed, the shift in the elastic tensors for the biological phases imposes the use of very refined mesh grids to properly capture the biofilm dynamic. Thus, it would be of particular interest to design and use well-balanced numerical scheme able to preserve the biofilm front structure. For more details about the numerical convergence of the scheme, see appendix 6.

Including light intensity

A microalga is a photosynthetic organism. Thus, microalgae require light to grow. When microalgae develop within a biofilm, the upper layers overshadow the lower layers. Following [START_REF] Clarelli | Mathematical models for biofilms on the surface of monuments[END_REF][START_REF] Clarelli | A fluid dynamics model of the growth of phototrophic biofilms[END_REF][START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF] to account for these mechanisms, the microalgae growth rate becomes

Including light intensity and solutes

Following [START_REF] Polizzi | A time-space model for the growth of microalgae biofilms for biofuel production[END_REF], let us now include three dissolved components: the substrate (S), the inorganic carbon (C), and the oxygen (O). As mentioned in subsection 2.1, the dynamic for dissolved components is modeled using a convection-diffusion reaction equation (6.9).

In a nutshell, the substrate represents the nitrate which is a nutrient of primary importance for the growth of autotrophic organisms like microalgae. Besides, roughly speaking, photosynthesis is the assimilation of inorganic carbon using light energy by autotrophic organisms. Photosynthesis releases oxygen. Thus, including these components is of primary interest. Taking into account these compounds also allows us to include the process of respiration. Basically, respiration is the opposite mechanism of photosynthesis and its consideration allows us to better describe the dynamic of thick biofilms. Indeed, the process of respiration becomes non-negligible in the absence of light, namely in the biofilm's inner layers.

As for the light, the contributions of the dissolved components to the photosynthesis process are accounted for in the growth through the multiplication by functions f p , p ∈ {S, C, O} which represent how the growth is modified by the local concentration of these components. On the one hand, limited contribution in high concentration regimes of the substrate and the inorganic carbon is modeled using Monod's law: f p = θp Kp+θp . On the other hand, the inhibition induced by high oxygen concentration is modeled by the sigmoidal function

Thus, including the contribution of the dissolved components and the light intensity, the algae growth rate becomes:

The respiration process is modeled by ψ r = µ r ϕ A θ O Kr+θ O where µ r is the maximal respiration rate and K r the half-saturation constant for the oxygen.

The modification of the microalgae growth rate and the inclusion of the respiration process requires to adapt the source terms for the phases as follows:

As for a phase, the source terms for a dissolved component is the sum of the pseudo-stoichiometric coefficients multiplied by the reaction rates. Thus, for the dissolved components, the source terms are

Appendix

Spatiotemporal equilibrium

The spatiotemporal equilibrium states for the system (6.12) correspond to the state solution where the source terms of all phases vanish, namely: Γ i = 0 for i ∈ {A, E, L}. In particular, Γ E = 0 induces ψ e + η E ψ d = 0 which lead to ϕ A = 0. Thus, the only spatiotemporal stationary state is the null state, namely ϕ i = 0 for i ∈ {A, E, L}.

Time dynamic of the volume fraction sum

Numerical convergence analysis for the model including viscosity

Numerical experiments have shown that numerical parameters need to be significantly reduced to reach convergence when the viscosity is included, and re-estimating the elastic tensor accordingly. Indeed, as mentioned in section 5.1, when including the viscosity, the elastic tensor coefficients must be rescaled and multiplied by 9 • 10 7

Numerical schemes for mixture theory models to obtain realistic front velocities for the biofilm. Figure 6.8 shows the convergences of the numerical scheme in both cases: with and without the viscous term. As expected, numerical convergence is obtained in both cases.

Nevertheless, as presented in Figure 6.8 the convergences rate is lower when the viscous term is included. This explains at least partially why the numerical parameters need to be significantly reduced to reach acceptable precision for the application considered when the viscosity is included

Relative variation of solutes concentration

Conclusion

What about the continuous problem?

Our control problem is formulated as a spatially discretized model, which is classically called a "discretize then optimize" approach. It would be interesting to compare it with the "optimize then discretize" approach, which consists in using a continuous optimization algorithm that is then applied to a discretized model. This would raise the question of the existence and uniqueness of optimal control, as well as its regularity, and it would also be interesting to study first-order optimality conditions in this framework. Furthermore, we have shown that for any spatial discretization step ∆x, the associated optimal control problem P ∆x has a solution u ∆x . This leads to the following natural questions, related to what happens when this step shrinks to zero:

• Does u ∆x converge to some u 0 when ∆x → 0?

• What is the nature of this convergence?

• Is there a topological dependence of the control on ∆x?

A more general question concerns the Γ-convergence of solutions (see e.g. [START_REF] Attouch | Variational analysis in Sobolev and BV spaces: applications to PDEs and optimization[END_REF][START_REF] Braides | Gamma-convergence for Beginners[END_REF]), which would imply at the same time the convergence of minimizers as ∆x → 0, and the one of the value function. This question would help justify the use of a discrete model as a good approximation to the continuous model derived from Coclite et.

al. [START_REF] Coclite | Traffic flow on a road network[END_REF]. It is also expected that a better understanding of the continuous problem will not only provide the basis for new algorithms, but also guarantee the convergence and stability of solutions.

Scaling the problem up

An intriguing research objective lies in the exploration of strategies for optimizing traffic flow within sprawling metropolitan road networks. However, a formidable obstacle currently looms large, primarily related to the substantial computational resources demanded by the control algorithm discussed in chapter 2 when dealing with extensive network configurations.

Starting with the extension of the deep hyper-reduction approach to 1x2, 2x1 and 2x2 junctions, the synthesis of the main elements of this thesis into a coherent methodology could be an interesting way to tackle this numerical challenge. In this approach, we would construct a simplified model tailored to accommodate the inherent complexities of significantly larger networks, thereby mitigating the computational demands. Within the framework of this lightweight but robust model, we would conduct controlled experiments to discern effective traffic control strategies.

123 Numerical schemes for mixture theory models Afterwards, we hope that the findings and lessons learned from these controlled experiments using the scaled-down model can potentially be extended and used to address the practical issues of traffic management in metropolitan areas. This approach aims to bridge the gap between theoretical research and real-world applications, potentially leading to enhanced efficiency in managing complex urban traffic networks.

It is paramount to emphasize that the successful implementation of this methodology relies on a thorough exploration of computational methodologies, modeling techniques and validation protocols, guaranteeing both its feasibility and its effectiveness in the field of real-world metropolitan traffic optimization.

A complementary way of scaling would be to design a meta-graph of a complex network to reduce the number of paths to be considered in the control loop. We can imagine considering a precise network only in a neighborhood of the axis to be evacuated and factorizing the rest of the graph to obtain non-physical meta-routes, but computing the appropriate boundary conditions to solve the control problem. This type of partitioning was studied in particular in S. V. Dongen's thesis [START_REF] Dongen | Graph Clustering by Flow Simulation[END_REF], where he introduced the Markov Clustering Process (MCL process). This consists of empirically averaging random flows traversing a graph and grouping the most frequently connected nodes into a single class. In fact, such classes appear to be sufficiently isolated from the rest of the graph that they can be treated as hypersummits, which would greatly reduce the size of the problem.

The Figure 6.10 show us the reduction of a 150-node graph to 14 simple subgraphs, which we could imagine treating as meta-summits or independent subgraphs, at least with greatly reduced coupling to allow efficient parallelization of computations.

Beyond road traffic

As pointed out by Ciro D'Apice, Rosanna Manzo, and B. Piccoli in [START_REF] D'apice | Packet flow on telecommunication networks[END_REF], the way we model traffic on road networks shares similarities with how we model telecommunication networks or logistical supply chains.

Therefore, it could be worthwhile to explore the possibility of applying the tools developed in this thesis to these contexts. This could help in investigating issues related to network congestion and finding optimal ways to efficiently utilize the network links, all within a resource-efficient approach that minimizes the resources needed for their operation.