Gabriel R Barrenechea

Franz Chouly
email: franz.chouly@inria.fr

A finite element method for the resolution of the Reduced Navier-Stokes/Prandtl equations

Keywords: incompressible flow, RNS/P equations, finite elements

A finite element method to solve the Reduced Navier-Stokes Prandtl (RNS/P) equations is described. These equations are an asymptotic simplification of the full Navier-Stokes equations, obtained when one dimension of the domain is of one order smaller than the others. These are therefore of particular interest to describe flows in channels or pipes of small diameter. A low order finite element discretization, based on a piecewise constant approximation of the pressure, is proposed and analyzed. Numerical experiments which consist in fluid flow simulations within a constricted pipe are provided. Comparisons with Navier-Stokes simulations allow to evaluate the performance of prediction of the finite element method, and of the model itself.

Introduction

For some kinds of flow problems, it is physically relevant to simplify the full Navier-Stokes equations assuming that one or two characteristic lengths are predominent. Among the most classical examples are the Prandtl's boundary layer equations (see e.g. [START_REF] Schlichting | Boundary-layer theory[END_REF]) or the hydrostatic approximation for shallow water flows (see e.g. [START_REF] Azérad | Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics[END_REF]). The interest of these simplifications is a refined analysis of the fluid flow problem, through a better understanding of its relevant scalings. It can be also motivated by practical considerations, as it might reduce the computation cost far above those of the most efficient Navier-Stokes solvers, despite the considerable effort that has been made in this direction (see e.g. [START_REF] Turek | Efficient Solvers for Incompressible Flow Problems[END_REF]).

For flows in long pipes of small diameter, such a simplification, justified by an asymptotic analysis, has been derived and called Reduced Navier-Stokes/Prandtl (RNS/P) equations [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF]. These equations have been employed successfully to model various problems, specifically in the domain of biological flows: in a stenosis [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF], in the laryngeal glottis [START_REF] Lagrée | Characterization of the pressure drop in a 2D symmetrical pipe: some asymptotical, numerical and experimental comparisons[END_REF], in the pharyngeal duct [START_REF] Van Hirtum | In vitro validation of some flow assumptions for the prediction of the pressure distribution during obstructive sleep apnea[END_REF][START_REF] Lagrée | Asymmetrical effects in a 2D stenosis[END_REF]. In particular, the prediction of pressure and wall shear stress distributions has been compared to references such as Navier-Stokes simulations [START_REF] Lagrée | Asymmetrical effects in a 2D stenosis[END_REF] or experimental measurements [START_REF] Van Hirtum | In vitro validation of some flow assumptions for the prediction of the pressure distribution during obstructive sleep apnea[END_REF]. It results from these comparisons that the RNS/P predictions are quite close to the references, with a relative error of a few percents.

For these equations, a numerical method based on finite differences has been proposed and tested in [START_REF] Lagrée | Characterization of the pressure drop in a 2D symmetrical pipe: some asymptotical, numerical and experimental comparisons[END_REF]. It is a streamwise marching algorithm, inspired from the classical methods designed to solve the heat equation [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. Even if the proposed method is cheap and adaptable to some different geometries, it has some important drawbacks. Among them, we can quote:

-first, it lacks of robustness. Specifically, for some categories of geometries, such as constrictions, numerical problems occur after the separation of the flow. This is mostly due to recirculation effects, which cannot be easily taken into account in the finite differences framework. The standard method is to use the "FLARE" approximation [START_REF] Reyhner | The interaction of a shockwave with a laminar boundary layer[END_REF] which consists in removing the u∂ x u term when the longitudinal velocity is negative [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF]. However, this is an ad-hoc approximation and does not ensure a correct computation in the whole domain.

-Some care has to be taken when adapting a finite differences scheme to complicated geometries, or to a three-dimensional problem [START_REF] Briley | Numerical method for predicting three-dimensional steady viscous flow in ducts[END_REF].

-If we are interested in fluid-structure interaction problems (such as in the upper airways [START_REF] Chouly | Numerical and experimental study of expiratory flow in the case of major upper airway obstructions with fluid-structure interaction[END_REF]), with the finite element method for solving the motion of the structure, the transmission of the forces at the interface can not be done in a simple and natural way (see [START_REF] Chouly | Numerical and experimental study of expiratory flow in the case of major upper airway obstructions with fluid-structure interaction[END_REF] for the details).

As a result of the previous considerations, in this work we are interested in the first steps towards a finite element method for the resolution of the RNS/P equations, which will avoid some of these disadvantages. Even if the presentation and the results are given for the bidimensional case, the method can be easily extended to the tridimensional case.

The plan of the paper is as follows. First, the complete boundary value problem is given in Section 2. The finite element method is described and analyzed in Section 3; in particular, since finite elements of common use for the Navier-Stokes equations -namely Taylor-Hood element [START_REF] Zienkiewicz | The Finite Element Method for Fluid Dynamics -Sixth Edition[END_REF] and the Mini element [START_REF] Ern | Theory and Practice of Finite Elements[END_REF] -do not provide a correct approximation, we use a specific finite element, originally proposed in [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF] for the Stokes equation. It is shown that with this method, the discrete problem admits a solution. Numerical experiments, presented in Section 4, have been carried out to confirm the analysis, and to test the precision of the method through comparison with Navier-Stokes simulations, taken as a reference. Finally some concluding remarks and perspectives are drawn.

The boundary value problem

The Reduced Navier-Stokes/Prandtl (RNS/P) equations are derived from the Navier-Stokes (1) equation. For the sake of simplicity, one can assume a newtonian, steady, incompressible, laminar and bidimensional flow:

(u • ∇) u = -1 ρ ∇p + ν△u + g, ∇ • u = 0, ( 1
)
where u is the velocity, p is the pressure, ρ is the density, ν is the kinematic viscosity and g is the external force field; g is in a great amount of applications the gravity field but may also stand for any kind of other external influence (e.g. a magnetic field). To derive the RNS/P equations, we need two assumptions, namely: 1. if we note D 2 the transversal dimension of the domain and D 1 the longitudinal dimension, then D 1 /D 2 ≫ 1 (Fig. 1). 2. if the Reynolds number is defined as Re = U 0 D 2 /ν, where U 0 stands for the maximal velocity at the entry, then Re ≫ 1 .

Then, the Navier-Stokes equation ( 1) can be simplified in order to obtain Reduced Navier-Stokes / Prandtl (RNS/P) equations (see [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF] for the derivation):

   u 1 ∂ x 1 u 1 + u 2 ∂ x 2 u 1 = -1 ρ ∂ x 1 p + ν ∂ 2 x 2 2 u 1 + g 1 , ∂ x 2 p = 0, ∂ x 1 u 1 + ∂ x 2 u 2 = 0. (2)
Here, (u 1 , u 2 ) are respectively the longitudinal and the transversal components of the velocity u, and g 1 is the longitudinal component of the external force field g. In the case of a gravity field, it means of course that the gravity is taken into account only if the duct is not horizontal. Boundary conditions consist of no slip on the lower and upper walls as well as an inlet flow at the entrance. The exit of the domain is considered as free. As a result, the RNS/P equations are the Prandtl boundary layer equations [START_REF] Chorin | A Mathematical Introduction to Fluid Mechanics -Third Edition[END_REF] with two major differences:

1. the domain in the RNS/P formulation is bounded in the transverse direction and there is no more fitting at the infinity with the inviscid flow; 2. the pressure distribution in the domain is an unknown.

Let us consider Ω which is a polygonal domain in R 2 with boundary ∂Ω; Γ i ⊂ ∂Ω is the entry (inlet flow), Γ w ⊂ ∂Ω is the rigid wall (with no-slip boundary conditions) and Γ o ⊂ ∂Ω is the exit (outlet flow) (Fig. 1). We give now the full boundary value problem that aims to be solved, in a non-dimensional form:

                   u 1 ∂ x 1 u 1 + u 2 ∂ x 2 u 1 + ∂ x 1 p -1 Re ∂ 2 x 2 2 u 1 = g 1 in Ω, ∂ x 2 p = 0 in Ω, ∂ x 1 u 1 + ∂ x 2 u 2 = 0 in Ω, u 1 = u 0 1 on Γ i , u 1 = 0 on Γ w , u 2 n 2 = 0 on Γ i ∪ Γ w , σ RNSP n = 0 on Γ o . (3)
The velocity profile u 0 1 at the entry may be arbitrary, usually a flat profile or a Poiseuille profile. Note that the condition u 2 n 2 = 0 comes from the fact that u 2 is expected to be less regular than u 1 (as in [START_REF] Azérad | Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics[END_REF]), so that the classical trace theorems (see e.g. [START_REF] Raviart | Introduction à l'Analyse Numérique des Équations aux Dérivées Partielles[END_REF]) do not allow to define u 2 on the boundary (in next section is introduced the space in which u 2 is defined). The condition on Γ o is a "do-nothing" condition, where σ RNSP = -p I + 1

Re ∂ x 2 u 1 e 1 ⊗ e 2 is a degenerated Cauchy stress tensor associated to the RNS/P equations, and n the normal unit vector oriented outward the domain.

The finite element method

As in the case of Stokes or Navier-Stokes equation, a finite element method can be proposed to solve the RNS/P equations. Moreover, some of the techniques already known to obtain a discrete approximation and to analyze it can be adapted to this case. Nevertheless, as in the case of the primitive equations of the ocean, in which the pressure is also constant in one direction of the space, finite elements such as Taylor-Hood or Mini element are inappropriate for discretization [START_REF] Azérad | Analyse des équations de Navier-Stokes en bassin peu profond et de l'équation de transport[END_REF]. Hence, we propose a discretization using an element which was first studied in [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF].

Weak formulation

To avoid technical difficulties, the problem (3) is rewritten with homogeneous Dirichlet boundary conditions on ∂Ω (Γ w = ∂Ω, Γ i = ∅, Γ o = ∅) 3 . Let us first present some notations. By L 2 (Ω) we denote the space of square integrable scalar functions on Ω, (•, •) Ω stands for the inner product in

L 2 (Ω) (in L 2 (Ω) 2 or in L 2 (Ω) 2×2 , if necessary); ||.|| 0,Ω stands for the norm in L 2 (Ω) associated to (•, •) Ω . L 2 0 (Ω)
is the subspace of functions in L 2 (Ω) with zero mean value on Ω. H 1 (Ω) is the space of square integrable scalar functions on Ω, with square integrable first derivatives. In the sequel, we will need the following space

H 1 (∂ x 2 , Ω) = {v ∈ L 2 (Ω) | ∂ x 2 v ∈ L 2 (Ω)}. (4)
H 1 0 (Ω) stands for the closed subspace of H 1 (Ω) with vanishing trace on ∂Ω. Similarly, we note H 1 0 (∂ x 2 , Ω) the following space:

H 1 0 (∂ x 2 , Ω) = {v ∈ H 1 (∂ x 2 , Ω) | vn 2 | ∂Ω = 0}. (5)
H 1 (Ω), H(Ω) and H 0 (Ω) are the following spaces of vector-valued functions:

H 1 (Ω) = H 1 (Ω) × H 1 (Ω), H(Ω) = H 1 (Ω) × H 1 (∂ x 2 , Ω), H 0 (Ω) = H 1 0 (Ω) × H 1 0 (∂ x 2 , Ω). (6)
H 1 (Ω) is a Hilbert space with the following scalar product:

(u, v) 1,Ω = (u, v) Ω + (∇u, ∇v) Ω . (7)
H(Ω) is also a Hilbert space with the scalar product:

(u, v) H(Ω) = (u, v) Ω + (∇u 1 , ∇v 1 ) Ω + (∂ x 2 u 2 , ∂ x 2 v 2 ) Ω . (8)
This implies the following property on the norms:

∀v ∈ H 1 (Ω), ||v|| H(Ω) ≤ ||v|| 1,Ω , (9)
with ||.||

c : (u, v, w) -→ (u 1 ∂ x 1 v 1 + u 2 ∂ x 2 v 1 , w 1 ) Ω , a Re : (u, v) -→ 1 Re (∂ x 2 u 1 , ∂ x 2 v 1 ) Ω , a λ : (u, v) -→ (λ∇ • u, ∇ • v) Ω , a : (u, v) -→ a Re (u, v) + a λ (u, v). (10)
Here, λ is a non-negative (possibly equal to zero) scalar field over Ω. The least-squares term a λ can be added into the variational formulation without affecting the solution. Although irrelevant for the continuous problem, this least-squares term has been introduced as a standard stabilization for convection-dominated flow problems (see, e.g., [START_REF] Braack | Stabilized finite element methods for the generalized Oseen problem[END_REF] for a recent review on this issue) an has an effect on the solution of the discrete problem (see Section 4 for a discussion). We also introduce the continuous bilinear form b(•,

•) : L 2 0 (Ω) × H 0 (Ω) → R defined by: b : (p, v) -→ -(p, ∇ • v) Ω . (11)
Using these forms, we present the following weak formulation for (3):

Find (u, p) ∈ H 0 (Ω) × L 2 0 (Ω) such that:    ∀v ∈ H 0 (Ω), c(u, u, v) + a(u, v) + b(p, v) = (g, v) Ω , ∀q ∈ L 2 0 (Ω), b(q, u) = 0. ( 12
)
Note that from the boundary value problem, we have g = (g 1 , 0). Nevertheless, we will consider the general case g 2 = 0 in the rest of the text.

in H(Ω) for which trace vanishes on Γ w . The inlet flow is imposed by setting u 1 equal to u 0 1 on Γ i . The test function v belongs to H Γw ∪Γ i (Ω) since the "do-nothing" condition is chosen for the outlet flow on Γ o . Due to this condition, the functions p and q are not restricted to L 2 0 (Ω) and are now in L 2 (Ω).

We now return to the case of homogeneous Dirichlet condition for the analysis of the problem, and we introduce the following space:

H b (Ω) DEF = ker b = {v ∈ H 0 (Ω) | ∀q ∈ L 2 0 (Ω), b(q, v) = 0} = {v ∈ H 0 (Ω) | ∇ • v = 0}, ( 13
)
where the last equality arises from the fact that ∇ • v ∈ L 2 0 (Ω). As for the Navier-Stokes and Stokes equations, we reformulate this mixed weak formulation into two dependent problems. For stability reasons (see Lemma 2 below), the convective term c(•, •, •) is transformed using the following well-known relationship:

∀u, v, w ∈ H 0 (Ω), c s (u, v, w) = - 1 2 (∇ • u, v 1 w 1 ) Ω . ( 14
)
As a result, c s (•, •, •), the symmetric part of c(•, •, •), vanishes on H b (Ω) and the trilinear form is equal to its antisymmetric part:

∀u, v, w ∈ H b (Ω), c(u, v, w) = c a (u, v, w) = 1 2 ((u 1 ∂ x 1 v 1 + u 2 ∂ x 2 v 1 , w 1 ) Ω -(u 1 ∂ x 1 w 1 + u 2 ∂ x 2 w 1 , v 1 ) Ω ) . (15)
Using ( 14)-( 15) we can propose the following equivalent formulation for the problem [START_REF] Franca | Stabilized finite element methods. II. The incompressible Navier-Stokes equations[END_REF]:

-Find u ∈ H b (Ω) such that : ∀v ∈ H b (Ω) , c a (u, u, v) + a(u, v) = (g, v) Ω , (16)
-Find p ∈ L 2 0 (Ω) such that :

∀v ∈ H 0 (Ω) , b(p, v) = (g, v) Ω -c a (u, u, v) -a(u, v). (17)

Finite element spaces

For the continuous problem ( 16)-( 17), a discretization first proposed in [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF] for the Stokes problem has been chosen: a discontinuous approximation of the pressure is preferred as it provides local conservation of the mass (cf. [START_REF] Gresho | Incompressible Flow and the Finite Element Method[END_REF]). Let {T h } h>0 be a regular family of admissible triangulations of Ω (cf. [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]). For each K ∈ T h , h K stands for the element diameter and h = max (K∈T h ) h K . The following space has been chosen for the velocity field:

H h = (H h,2 ∩ H 1 0 (Ω)) × (H h,1 ∩ H 1 0 (∂ x 2 , Ω)), (18)
where, for k = 1, 2:

H h,k = {v ∈ C 0 ( Ω) | ∀K ∈ T h , v| K ∈ P k (K)}. ( 19
)
We also need to introduce the space H h,b , defined as follows:

H h,b = {v h ∈ H h | ∀q h ∈ Π h , b(q h , v h ) = 0}. ( 20
)
We note that H h,b is not necessarily a subspace of H b (Ω). The pressure is approximated using the following space:

Π h = {q ∈ L 2 0 (Ω) | ∀K ∈ T h , q| K ∈ P 0 (K)}. ( 21
)
The finite element associated to this choice is called P 2 /P 1 /P 0 . Using this pair of spaces, we propose the following finite element method for ( 16)-( 17):

-Find u h ∈ H h,b such that:

∀v h ∈ H h,b , c a (u h , u h , v h ) + a(u h , v h ) = (g, v h ) Ω . (22)
-Find p h ∈ Π h such that:

∀v h ∈ H h , b(p h , v h ) = (g, v h ) Ω -c a (u h , u h , v h ) -a(u h , v h ). ( 23
)

Analysis of the discrete problem

The aim of this section is to analyze the discrete problem ( 22)-( 23). It will be proved that it admits at least one solution. We start with the following technical, but fundamental result:

Lemma 1. The mapping

||.|| x 2 : H h,b -→ R v h = (v h,1 , v h,2 ) -→ ||v h || x 2 = ||∂ x 2 v h,1 || 0,Ω , (24)
defines a norm on H h,b .

Proof. The only property to check is that

||v h || x 2 = 0 implies v h = 0 in Ω. The other properties arise directly from the fact that ||.|| 0,Ω is a norm on L 2 (Ω). Let us consider v h ∈ H h,b such that ||v h || x 2 = 0.
The Poincaré inequality (see [START_REF] Raviart | Introduction à l'Analyse Numérique des Équations aux Dérivées Partielles[END_REF])

||v h,1 || 0,Ω ≤ C(Ω) ||∂ x 2 v h,1 || 0,Ω , (25)
implies that v h,1 = 0 in Ω. As a result, v h,2 satisfies ∀q h ∈ Π h , (q h , ∂ x 2 v h,2 ) Ω = 0. ( 26
)
Let us now remark that the function v h,2 can be written as follows on each triangle K i of the mesh:

v h,2 | K i (x 1 , x 2 ) = α i + β i x 1 + γ i x 2 . ( 27
)
For two given elements K i and K j of the mesh, the function

qh = 1 |K i | 1 1 K i - 1 |K j | 1 1 K j , (28)
belongs to Π h , and then, using [START_REF] Van Hirtum | In vitro validation of some flow assumptions for the prediction of the pressure distribution during obstructive sleep apnea[END_REF] we easily see that

γ i = γ j . ( 29
) Let D = {(x 0 1 , x 2 ) | x 2
∈ R} be a vertical line, for any arbitrary x 0 1 such that D ∩ Ω = ∅. This line intersects the mesh in a sequence of adjacent triangles (K i ) i∈{0,...,n} . The property (29) and the continuity of v h,2 imply that:

∀i ∈ {0, . . . , n}, α i = α 0 , β i = β 0 , γ i = γ 0 . ( 30
)
The boundary conditions that satisfies v h,2 are then such that: α 0 = β 0 = γ 0 = 0. Since the same argument may be used for every (or almost every)

x 0 1 such that D ∩ Ω = ∅, then v h = 0 in Ω. Remark. A consequence of this lemma is that (H h,b , (•, •) x 2 ) is a Hilbert space with the scalar product : (u h , v h ) x 2 = (∂ x 2 u h,1 , ∂ x 2 v h,1 ) Ω .
Lemma 2. For λ ≥ 0, the problem ( 22) admits at least one solution u h ∈ H h,b .

Proof. We follow an approach similar to the one presented in [START_REF] Temam | Navier-Stokes Equations[END_REF] for the full Navier-Stokes equations. In (H h,b , (•, •) x 2 ), finite dimensional Hilbert space, we introduce the mapping f from H h,b into itself as follows. For v h ∈ H h,b , f (v h ) is the unique vector such that:

∀w h ∈ H h,b , (f (v h ), w h ) x 2 = c a (v h , v h , w h ) + a(v h , w h ) -(g, w h ) Ω . (31)
It is easy to check that f is a continuous mapping. Now, for λ ≥ 0, a(•, •) satisfies

∀v h ∈ H h,b , a(v h , v h ) ≥ 1 Re ||v h || 2 x 2 , (32)
and, using the Cauchy-Schwarz inequality and the equivalence of norms in a finite dimensional space there exists a positive constant C 1,h such that:

∀v h ∈ H h,b , |(g, v h ) Ω | ≤ C 1,h ||g|| 0,Ω ||v h || x 2 . ( 33
)
Using the previous results, f satisfies

(f (v h ), v h ) x 2 = c a (v h , v h , v h ) + a(v h , v h ) -(g, v h ) Ω = a(v h , v h ) -(g, v h ) Ω ≥ 1 Re ||v h || 2 x 2 -C 1,h ||g|| 0,Ω ||v h || x 2 ≥ ||v h || x 2 ( 1 Re ||v h || x 2 -C 1,h ||g|| 0,Ω ). If we choose k > Re C 1,h ||g|| 0,Ω , then for ||v h || x 2 = k, (f (v h ), v h ) x 2 > 0.
As a result, the lemma 1.4 p.164 in [START_REF] Temam | Navier-Stokes Equations[END_REF] ensures the existence of a solution u h of the equation f (u h ) = 0, in other words, a solution of the discrete problem [START_REF] Shome | Modeling of airflow in the pharynx with application to sleep apnea[END_REF].

Remark. For the inequality (33), a better majoration can be given in the case g 2 = 0. Indeed, from the Cauchy-Schwarz inequality and the Poincaré inequality, we have:

∀v h ∈ H h,b , |(g, v h ) Ω | ≤ C(Ω)||g 1 || 0,Ω ||v h || x 2 . (34)
Hence, if we define k 0 = (Re C(Ω) ||g 1 || 0,Ω ), the Lemma 1.4 in [START_REF] Temam | Navier-Stokes Equations[END_REF] also ensures that ||u h || x 2 ≤ k 0 , thus ||u h,1 || 0,Ω ≤ C(Ω) k 0 , using again the Poincaré inequality. In other words, the set {||u h,1 || 0,Ω } (h>0) is bounded.

For the problem [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF], we now have the following lemma:

Lemma 3. The pair P 2 /P 1 /P 0 is inf-sup stable, i.e., there exists a constant β > 0, independent of h, such that:

inf q h ∈Π h sup v h ∈H h b(q h , v h ) ||v h || H(Ω) ||q h || 0,Ω ≥ β. (35)
Then, for a given u h , the problem ( 23) admits one unique solution p h ∈ Π h .

Proof. In [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF] it is proved that inf

q h ∈Π h sup v h ∈H h b(q h , v h ) ||v h || 1,Ω ||q h || 0,Ω ≥ β, (36)
and the result arises from [START_REF] Davis | Algorithm 832: UMFPACK V4.3-an unsymmetric-pattern multifrontal method[END_REF].

Collecting the previous results, we can state the main theorem of this section:

Theorem 1. The problem ( 22)-( 23) admits at least one solution (u h , p h ). Furthermore, in the case of g 2 = 0, the set {||u h,1 || 0,Ω } (h>0) is bounded by (C(Ω) 2 Re ||g 1 || 0,Ω ), where C(Ω) is the constant from the Poincaré inequality.

Remark :

(1) For the Taylor-Hood element and the Mini element, the inf-sup condition is also valid, which ensures that the discrete problem (23) has a unique solution for a given u h . Nevertheless, for these elements, and λ = 0, the problem [START_REF] Shome | Modeling of airflow in the pharynx with application to sleep apnea[END_REF] might have no solution. Indeed, Lemma 2 might not be valid since for these elements, ||.|| x 2 might not be a norm on H h,b (note that the specific properties of the P 2 /P 1 /P 0 element have been used in the proof of Lemma 1). This has been confirmed by the numerical experiments that fail for these elements.

(2) Note that no majoration of the transverse velocity u h,2 has been provided. This is due to the very particular nature of the RNS/P equations, that allow a weak control on this variable.

Description of the algorithm of resolution

For numerical simulations, boundary conditions that are not homogeneous have been considered (see equation [START_REF] Berger | Flows in stenotic vessels[END_REF] and remark below equation ( 12)). For the inlet flow, a Poiseuille profile has been chosen:

u 0 1 (x 2 ) = 4(1 -x 2 )x 2 . ( 37
)
The Newton method has been used to deal with the non-linearity that arises from the convective term. At each step of the Newton loop, the linearized discrete problem is solved using a multi-frontal Gauss LU factorization (cf. [START_REF] Davis | A combined unifrontal/multifrontal method for unsymmetric sparse matrices[END_REF]) implemented in the package UMFPACK4 (cf. [START_REF] Davis | Algorithm 832: UMFPACK V4.3-an unsymmetric-pattern multifrontal method[END_REF]). The complete scheme of the numerical resolution is given in Fig. 2.

The first numerical parameter is λ, which is the coefficient for the least-squares term a λ (•, •).

It has been chosen for the numerical experiments as a constant and not as a scalar function.

The two other parameters are n Re and ε N which are respectively the number of steps in the continuation loop 1 and the convergence criterion for the Newton loop. This latter has been fixed to 10 -7 for all the simulations. The convergence in the Newton loop has been measured through computation of:

||(du h , dp h )|| max = max T h ||(du h , dp h )|| 2 , (38)
where (du h , dp h ) is the increment in each iteration, defined in Fig. 2. All the numerical results have been obtained using FreeFEM++ software [START_REF] Hecht | Freefem++ 2.16-1 documentation[END_REF].

Numerical results and discussion

The problem consists of computing the fluid flow in a constricted pipe, a type of geometry which corresponds to a great variety of situations: flow in a Venturi pipe [START_REF] Turek | Efficient Solvers for Incompressible Flow Problems[END_REF], in a collapsible tube [START_REF] Carpenter | Flow in collapsible tubes and past other highly compliant boundaries, chapter 2. Flows in Deformable Tubes and Channels[END_REF], in a stenosis [START_REF] Berger | Flows in stenotic vessels[END_REF], in the vocal folds [START_REF] Scherer | Intraglottal pressure profiles for a symmetric and oblique glottis with a divergence angle of 10 degrees[END_REF] or in the human pharynx [START_REF] Shome | Modeling of airflow in the pharynx with application to sleep apnea[END_REF], etc. The geometry can either be symmetric (for instance in a stenosis) or asymmetrical (for instance in the human pharynx or at the base of the tongue). Here, we have considered the asymmetrical problem. The characteristics of one representative simplified geometry are given in Fig. 3: it is a straight pipe which is constricted because of a bump in the upper border. The pertinent parameters for this type of problem are the width δ and the height h b of the bump, as well as the Reynolds number Re.

The simulations have been carried out for three types of geometries:

-(geometry 1) A long pipe with a slightly curved upper wall (δ = 5, h b = 0.2). It corresponds to an ideal case in which the assumptions of validity of the RNS/P equations should be encountered. -(geometry 2) A pipe with a small obstacle (δ = 0.1, h b = 0.2), a case described in [START_REF] Lagrée | The RNS/Prandtl equations and their link with other asymptotic descriptions: application to the wall shear stress scaling in a constricted pipe[END_REF]. It permits to test the method in a more realistic situation, with separation of the flow above a given Reynolds number. -(geometry 3) A severe constriction (δ = 0.5, h b = 0.5), a case described in [START_REF] Lagrée | Asymmetrical effects in a 2D stenosis[END_REF]. The interest is to test the method and the model itself in a situation corresponding to the limit of validity of the RNS/P equations.

The meshes for each case are depicted Fig. 4. The range for the Reynolds number Re is 1 -1000 for the geometry 1, 1 -500 for the geometry 2 and 1 -100 for the geometry 3. In a first set of experiments, the value for λ has been fixed to 0. For comparison, the complete Navier-Stokes equations have been solved, on the same geometry and with the same mesh.

The numerical strategy is the same as for RNS/P equations: a continuation method with a Newton loop to treat the non-linearity. The number of steps of the continuation method is always the same for Navier-Stokes and RNS/P. The convergence criterion ε N for the Newton loop is fixed to the same value as for the RNS/P equations: 10 -7 . For discretization, P 2 /P 1 Taylor-Hood elements have been considered (quadratic interpolation on the velocity and linear interpolation on the pressure) [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. No stabilization of the finite element approximation has been used in the convection-dominated regime, since no pure oscillations have been observed in the numerical solutions. The pressure drop ∆P between the inlet and the outlet, which is an output of the simulations, has been compared. Moreover, the force F sup exerted by the fluid on the upper wall has been computed since it is of particular interest in the case of fluid-structure interaction. This force is defined as

F sup = Γsup σ f n dΓ, ( 39
)
where Γ sup is the upper part of the boundary, σ f is the tensor of fluid constraints and n is the inner unit vector normal to the boundary. It is of interest to decompose F sup as:

F sup = F p sup + 1 Re F τ sup , (40)
where F p sup is the contribution of the pressure:

F p sup = Γsup (-p n) dΓ, (41)
and F τ sup is the contribution of the shear stress:

F τ sup = Γsup (∇u + ∇u T )n dΓ. ( 42
)
In practice, for incompressible flows, the contribution from the shear stress is negligible with respect to the contribution from the pressure (see e.g. [START_REF] Chouly | Numerical and experimental study of expiratory flow in the case of major upper airway obstructions with fluid-structure interaction[END_REF]). The results of the computations and of the comparisons (pressure drop ∆P and quadratic norms of F sup , F p sup , F τ sup ) are presented in Table 1.

For the geometry 1 and for a Reynolds number of 1000, the pressure distribution p and the horizontal velocity u 1 are depicted Fig. 5. The predictions of ∆P and F sup are in good adequation with those from the Navier-Stokes simulations, taken as a reference. The maximal error is found for Re = 1000, and is of 2 % for ∆P and of 8 % for F sup . This error corresponds to the quantity

|∆P NS -∆P RNS/P | ∆P NS , ( 43
)
where NS indicates the prediction from Navier-Stokes simulations and RNS/P the prediction from RNS/P simulations (the same computation is done for F sup ). As a result, the simulations for this geometry have permitted to validate the numerical method.

For the geometry 2 and Re = 100, the pressure distribution p and the horizontal velocity u 1 are depicted in Fig. 6. For Re = 1, the prediction of ∆P and F sup corresponds to the prediction from Navier-Stokes simulations, with errors of 6 % and 1.5 %, respectively. When the Reynolds number Re is 100, and the convection such that recirculation effects are observed behind the obstacle in Navier-Stokes simulations, the prediction of ∆P and F sup remains satisfying (errors of 9 % and 2 % respectively), though the RNS/P equations are in principle not adapted for the simulation of recirculation effects, because of the assumption ∂ x 2 p = 0. The reason is that in this case, the global effect of recirculation is weak. Note however that the pressure distribution at the level of the bump differs between RNS/P and Navier-Stokes. Though the drop is approximatively of the same magnitude, it is anticipated in the RNS/P simulation, with an abrupt pressure recovery that is not present in the Navier-Stokes simulation. Note also that the effects of recirculation affect the values of the velocity and therefore the value of F τ sup , which is different between RNS/P and Navier-Stokes simulations (error of 37.5 %). For a Reynolds number Re of 500, the recirculation is stronger. This does not prevent the RNS/P simulation to converge but the results are quite different from those of the Navier-Stokes simulation: an error of 24 % for ∆P and an error of 36 % for F sup .

For the geometry 3 and Re = 1, p and u 1 depicted Fig. 7. The adequation between RNS/P and Navier-Stokes simulations is still satisfying for Re = 1, with an error of 13.5 % for ∆P and an error of 8 % for F sup . This is slightly higher than in the precedent cases, but this fact is somehow expected since for this kind of geometry the RNS/P equations are in their limit of validity. Due to the recirculation effects, the error is more important when Re is increased to 50 or 100, and is up to 23 % for ∆P and 21 % for F sup (for Re = 100).

Computing times5 for each simulation are indicated in Table 1 for λ = 0. It appears from the results that for Re ≫ 1, the RNS/P simulations are faster than the NS simulations, of approximatively 30 %. This gain comes from the fact that less Newton iterations are required with RNS/P to reach the same convergence criterion. In the case Re = O(1), the computing time for RNS/P is comparable to the computing time for NS, except for the geometry 3 in which it is higher for RNS/P. This might be due to the linear system which should be ill-conditioned for small Reynolds. This case is however not physically relevant, as the RNS/P equations have been designed for flows with high Reynolds number.

The influence of the least-squares term a λ (•, •) has finally been assessed. Simulations have been carried out with λ = 1 for the geometries 1 and 2. Simulations with the geometry 1 reveal the first interest of this least-squares term: it decreases significantly the computing time. Indeed, a speed-up of approximatively 2 is achieved in comparison to the reference simulations with Navier-Stokes (see Table 2). This term has in particular two effects: it makes easier the resolution of the linear system at each iteration, and it decreases the number of Newton iterations at each continuation step. Of course, for this geometry, this term does not improve the accuracy since results were already very close to those obtained from Navier-Stokes when λ was fixed to 0. Simulations for the geometry 2 also confirm these considerations on the computing time (see Table 3). The role of simulations on this geometry is then to show the other interest of the term a λ (•, •): it may improve the performance of prediction. For this geometry and Re ≤ 100, the results are very close to the ones obtained with λ = 0, so that the difference with the Navier-Stokes simulations is nearly the same. Concerning the prediction of the shear stress component F τ sup at Re = 100, the results are better with λ = 1, as the error is of 15 % instead of the 37.5 % mentioned previously for λ = 0. For Re > 100, the results presented in Table 3 show clearly that the performance of prediction is improved with λ = 1. Concerning ∆P , the error is of 7.5 % for λ = 1 instead of 14 % for λ = 0 when Re = 200 (respectively, 6 % and 24 % when Re = 500). For F sup , the error decreases from 5 % to 2 % when λ changes from 0 to 1, when Re = 200 (respectively from 36 % to 15 % when Re = 500). Indeed, when the Reynolds number is high and the convection effects become predominant in the fluid, the least-square term allows an ad-hoc reproduction of the recirculation behind the obstacle, though underestimated if compared to those observed in the Navier-Stokes simulations (see Fig. 8). This explains the better estimation of F τ sup (see Table 3) and F sup in particular.

As the choice λ = 1 is of course not the only possible, a sensitivity study has been carried out, which results are presented in Table 4. First, we noted that the value of this parameter has not a significant impact on the computing time when it is chosen different from 0: for λ ≥ 1, this time remained to be 1'21 (1'44 for λ = 0). Nevertheless, concerning the performance of prediction, the value of λ = 2 is the optimal for F sup (0.6 % of error) and λ = 3 is the optimal for ∆P (0.7 % of error), so a value of λ between 2 and 3 should be the best for carrying out simulations. Though increasing λ reduces the error on F τ sup , it increases the error on ∆P and F sup if λ > 3, so values of λ greater than 3 might not be a good choice.

Concluding remarks

A finite element method to solve the Reduced Navier-Stokes/Prandtl (RNS/P) equations has been described and tested. The discretization is based on an element originally defined in [START_REF] Stenberg | Error analysis of some finite element methods for the Stokes problem[END_REF] for the Stokes equations. With this element, called P 2 /P 1 /P 0 , the velocity is approximated with continuous piecewise quadratic and linear functions, while the pressure is approximated with a piecewise constant function. This element, as Taylor-Hood and the Mini element, verifies the inf-sup condition for the RNS/P equations. However, the subproblem of computing the discrete velocity u h might be ill-posed for the Taylor-Hood and the Mini element, which systematically leaded to failures in the numerical experiments. In opposite, we show that with the proposed P 2 /P 1 /P 0 element, this subproblem admits a solution. Furthermore, it was shown that the longitudinal component of u h is bounded by a constant independent of h, the parameter of discretization.

Numerical experiments permitted to test the method for some particular geometries, and to show that the prediction of the pressure drop and of the constraints on the surrounding walls is comparable to the prediction from Navier-Stokes simulation. Moreover, the current method avoids limitations of the precedent finite differences method used for instance in [START_REF] Lagrée | Characterization of the pressure drop in a 2D symmetrical pipe: some asymptotical, numerical and experimental comparisons[END_REF], which lacks of adaptivity for some geometries and for coupling with other physical entities, such as a moving wall. Also, a bidimensional problem has been chosen for the simplicity of the presentation, but the extension to three-dimensional geometries is straightforward. Even in the case of moderate recirculation effects (small obstacle in a duct for instance), the method still gives a correct approximation of the predicted variables, whereas in the finite difference context, the computation had to be stopped, and only a prediction of the pressure and of the velocity on the upstream part of the domain was provided [START_REF] Lagrée | Asymmetrical effects in a 2D stenosis[END_REF]. In the case of strong recirculation, for instance in a severe constriction, both finite differences and finite elements are unsatisfying. However, in this case, the RNS/P equations are out of their domain of validity, and a Navier-Stokes solver should be used instead for an accurate simulation.

Concerning the performance of the method, a gain has been observed systematically in comparison with Navier-Stokes simulation. The reason is that less Newton iterations are required to reach convergence with the RNS/P equations. Moreover, the least-squares stabilization term has a positive impact on the speed-up, by its effect on the resolution of the linear system at each Newton iteration, and as it decreases even more the number of Newton iterations. In addition, this gain should be increased significantly, more specifically in the three-dimensional case, by improving the method for the resolution of the linear system which is built at each iteration. For instance, it would be interesting to design preconditioners adapted to this system. Another possibility would consist in implementing and studying the effect of stabilization techniques much more sophisticated than the simple one that has been used (see e.g. [START_REF] Franca | Stabilized finite element methods. II. The incompressible Navier-Stokes equations[END_REF] for the Navier-Stokes case). These two points will be object of future research.

Another perspective concerns shallow water flows, for which the Navier-Stokes equations are simplified in a very similar manner [START_REF] Azérad | Analyse des équations de Navier-Stokes en bassin peu profond et de l'équation de transport[END_REF][START_REF] Azérad | Mathematical justification of the hydrostatic approximation in the primitive equations of geophysical fluid dynamics[END_REF]. Nevertheless, the term ∂ 2

x 2 1 u 1 is conserved in the laplacian, though it is not justified asymptotically. Thus, it would be interesting to study the applicability of the proposed method to shallow water equations. The last perspective is the extension of this finite element method to a coupled fluid/structure interaction problem and to three-dimensional geometries. Init:

• Read the simulation parameters and the mesh.

• Init the velocity and the pressure: (u h , p h ).

• Init the Reynolds number Re.

Loop 1: Continuation strategy.

• Loop 2: Newton iteration.

• Solve the linearized problem:

Find (du h , dp h ) such that

8 > > > > > > > > > > < > > > > > > > > > > : ∀v h , ca(u h , du h , v h ) + ca(du h , u h , v h ) + a(du h , v h ) -(dp h , ∇ • v h )Ω = (g, v h )Ω -ca(u h , u h , v h ) -a(u h , v h ) + (p h , ∇ • v h )Ω, ∀q h , -(q h , ∇ • du h )Ω = (q h , ∇ • u h )Ω.
• Update the velocity and the pressure:

u h ← u h + du h p h ← p h + dp h • End Loop 2.
• Increase the Reynolds number Re.

End Loop 1.

End.

Fig. 1 .

 1 Fig. 1. The domain Ω for the resolution of the RNS/P equations, with the notations for the different parts of the boundary. Note that x1 and x2 are non-dimensional coordinates.

Fig. 2 .Fig. 3 .Fig. 4 . 18 -Fig. 5 . 19 -Fig. 6 .Fig. 7 .

 2341851967 Fig. 2. The complete algorithm of numerical solving of the RNS/P equations.

Fig. 8 .

 8 Fig. 8. Fluid flow in a constricted pipe. Geometry 2. Effect of the least-squares term a λ (•, •). Re = 500. λ = 1. Is depicted: the horizontal velocity u1. (a) RNS/P simulation. (b) Navier-Stokes simulation.

  1,Ω the norm on H 1 (Ω) associated to (•, •) 1,Ω . Moreover, let us note c(•, •, •), a Re (•, •), a λ (•, •) and a(•, •) the continuous trilinear and bilinear forms on H 0 (Ω) defined by:

Table 1 .

 1 Results of the computations for the RNS/P equations (left) and the full Navier-Stokes equations (right). The pressure difference and the force on the superior wall are indicated, as well as the computing time (in minutes/seconds) of the simulations.

	geometry Re ∆P ||F sup|| ||F p sup || ||F τ sup || time	geometry Re ∆P ||F sup|| ||F p sup || ||F τ sup || time
		1 215.4 2218.0 2216.8 96.6 0'36		1 216.1 2174.8 2173.6 97.1 0'41
	(1)	100 2.16 20.80 20.79	96.6 0'45	(1)	100 2.16 20.36 20.35	97.0 1'12
		500 0.46	3.51	3.51	97.6 2'07		500 0.46	3.38	3.37	98.0 2'53
		1000 0.26	1.62	1.62	100.1 3'29		1000 0.25	1.51	1.50	100.3 5'10
		1 40.98 105.14 103.03 21.92 0'45		1 43.52 106.66 105.66 18.55 0'40
	(2)	100 0.42	1.03	0.99	21.87 1'03	(2)	100 0.46	1.01	1.00	15.91 1'23
		500 0.10	0.17	0.17	20.81 4'16		500 0.13	0.13	0.12	9.41 6'15
		1 71.33 157.66 150.72 31.96 1'16		1 82.47 171.20 165.16 26.58 0'38
	(3)	50 1.53	2.53	2.48	31.38 1'03	(3)	50 1.85	2.92	2.90	24.35 1'30
		100 0.85	0.94	0.94	30.95 4'32		100 1.10	1.19	1.19	21.16 6'20
	equations	∆P ||F sup|| ||F p sup || ||F τ sup || time	equations	∆P ||F sup|| ||F p sup || ||F τ sup || time
	Navier-Stokes 1.09 9.58	9.56	97.04 1'32	Navier-Stokes 0.46 3.38	3.37	98.0 2'53
	RNS/P (λ = 0) 1.09 9.81	9.79	96.61 1'10	RNS/P (λ = 0) 0.46 3.51	3.51	97.6 2'07
	RNS/P (λ = 1) 1.09 9.79	9.77	96.95 0'40	RNS/P (λ = 1) 0.46 3.46	3.45	97.8 1'22

Table 2 .

 2 Influence of the parameter λ for simulations with the geometry 1. The Reynolds number Re is 200 (left) and 500 (right).

	equations	∆P ||F sup|| ||F p sup || ||F τ sup || time	equations	∆P ||F sup|| ||F p sup || ||F τ sup || time
	Navier-Stokes 0.25 0.46	0.45	13.79 2'44	Navier-Stokes 0.13 0.13	0.12	9.41 6'15
	RNS/P (λ = 0) 0.21 0.48	0.48	21.56 1'44	RNS/P (λ = 0) 0.10 0.17	0.17	20.81 4'16
	RNS/P (λ = 1) 0.23 0.45	0.45	16.61 1'23	RNS/P (λ = 1) 0.12 0.15	0.15	13.09 3'16

Table 3 .

 3 Influence of the parameter λ for simulations with the geometry 2. The Reynolds number Re is 200 (left) and 500 (right).

	λ ∆P ||F sup|| ||F p sup || ||F τ sup || time
	0 14.1	5.3	5.5	56.3 1'44
	1 7.5	1.7	1.1	20.4 1'23
	2 3.6	0.6	0.1	13.7 1'21
	3 0.7	1.3	2.2	9.8	1'20
	4 1.8	3.4	4.5	7.0	1'21
	5 4.3	5.5	6.9	4.9	1'21

Table 4 .

 4 Study of the parameter λ for simulations with the geometry 2. The Reynolds number Re is 200. The relative error in % in comparison to Navier-Stokes simulations is given for each variable.

However, the weak formulation in the general case will be given in the next remark.

Remark. In the case of non-homogeneous boundary conditions, the weak formulation (12) is still valid with the following modifications: u belongs to H Γw (Ω), the subspace of functions

http://www.netlib.org/linalg

The computer on which simulations have been carried out is a Power PC G4 1.2 GHz, with RAM of 256 MO.

Acknowledgements

A part of this work was performed during Franz Chouly's postdoctoral stay at the Departamento de Ingeniería Matemática, Universidad de Concepción, Chile, funded by Alfa Project and Fundación Andes, Chile (grant C-14055). The authors would like to thank R. Rodríguez, C. Pérez, P.-Y. Lagrée, M. Fernández and J.-F. Gerbeau for many helpful discussions and comments.


A Eubanks

M Mimi

Dembélé Keita

C Anoma

T T E Dah

E Mensah

G Maradan

M Bourrelly

M Mora

L Riegel

Understanding how peer relationships influence peerdelivered HIV prevention interventions among Ugandan female sex workers: a case study

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

with index-reported test results (n = 105), 82% (23/28) who tested positive initiated treatment (80% HIVST vs 92% clinic-based; p = 0.60) and 20% (15/76) who tested negative initiated PrEP (16% HIVST vs 23% clinic-based; p = 0.57). Conclusions: HIVST kit distribution to newly diagnosed HIV-positive individuals effectively increased sexual partner testing, had high positivity yield and linkage to treatment. One-fifth of HIV-uninfected partners initiated PrEPthus innovations to link to prevention services are urgently needed.

PE16.15

Understanding how peer relationships influence peerdelivered HIV prevention interventions among Ugandan female sex workers: a case study from HIV self-testing Background: Female sex workers (FSWs) have tightly connected peer networks and are at high risk of HIV infection. Studies demonstrate that uptake and continuation of HIV prevention interventions increase when they are peer-delivered. We analyzed qualitative data from a peer-delivered HIV self-testing (HIVST) intervention among FSWs in urban Uganda to understand how underlying peer dynamics can facilitate or hinder the uptake of HIVST delivered by peers. Methods: Between October and November 2016, 960 FSWs were enrolled in a four-month randomized trial testing different peerdelivered HIVST models in Kampala. An additional 120 FSWs were trained as peer educators (PEs) to deliver HIVST. FSWs were >= 18 years old, self-reported exchanging sex for money or goods (past month) and had not recently tested for HIV (past 3 months). PEs either distributed HIVST kits (arm 1) or coupons to obtain HIVST kits from local clinics (arm 2) to FSWs. FSWs (n = 30) and PEs (n = 28) from these intervention arms were randomly sampled to complete indepth interviews and focus group discussions, respectively. Verbatim transcripts were analyzed thematically using inductive and deductive codes.

Results: The median age of FSWs was 28 years (IQR: 24 to 32) and PEs was 33 years (IQR: 29 to 37). We found that peer dynamics, including social hierarchies, communication networks, economic power, and HIV status disclosure, acted both as facilitators and barriers to HIVST uptake, Figure 1. For example, communication networks facilitated HIVST uptake by effectively relaying accurate information about HIVST procedures, but these networks were vulnerable to misinformation and rumors, thus hindering HIVST uptake. Conclusions: In Uganda, FSW peer networks have complex existing dynamics that act as facilitators and barriers to peer-delivered HIV prevention interventions. Future peer-delivered HIV prevention interventions should be carefully designed around existing dynamics within peer networks to maximize initial and repeat intervention uptake.

PE16.16

Loss to follow-up among MSM on PrEP in West Africa (CohMSM-PrEP ANRS12369-Expertise France)

  A. Eubanks 1 ; M. Mimi 1 ; B. Demb el e Keita 2 ; C. Anoma 3 ; T. T.E. Dah 4 ; E. Mensah 5 ; G. Maradan 6 ; M. Bourrelly 1 ; M. Mora 1 ; L. Riegel 7 ; D. Rojas Castro 7 ; I. Yaya 8 ; B. Spire 1 ; C. Laurent 8 and L. Sagaon-Teyssier 1 1 Aix-Marseille University, INSERM, IRD, SESSTIM, Marseille, France, 2 ARCAD-SIDA, Bamako, Mali, 3 Espace Confiance, Cote D'Ivoire, 4 Centre Muraz, Bobo-Dioulasso, Burkina Faso, 5 Espoir Vie Togo, Togo,

  University of Washington, Department of Global Health, Seattle, United States, 3 International Research Consortium, Kampala, Uganda, 4 University of California, Francis I. Proctor Foundation, San Francisco, United States

	M. McGowan 1 ; S. Roche 2 ; A. Nakitende 3 ; J. Wachinger 1 ;
	D.K. Musoke 3 ; C.E. Oldenburg 4 ; T. Bӓrnighausen 1 and K. Ortblad 2
	1 Heidelberg University, Heidelberg Institute of Global Health (HIGH),
	Heidelberg, Germany, 2

confronted by unknowns concerning the feasibility in this context, due to the highly vulnerable nature of MSM (stigma, precarity, high-risk sex). We aimed to estimate the attrition rate and identify the factors associated with loss to follow-up (LTFU) in a cohort of MSM on PrEP in West Africa. Methods: Since 2017, CohMSM-PrEP has offered a comprehensive prevention package for MSM in Mali, Cote d'Ivoire, Burkina Faso, and Togo. Quarterly follow-up includes PrEP (daily or event-driven) and socio-behavioral data collection. Participants from a previous MSM cohort and new participants were enrolled. LTFU was defined as not attending the last two scheduled follow-up visits. The Kaplan-Meier technique and log-rank test were used to estimate time to LTFU and to test for significance between groups. The Cox proportional hazards regression model was used to determine predictors of LTFU and adjusted by confounders. Results: 585 participants were recruited from November 2017-January 2020. The median follow-up time was 15.6 months. During this period, 119 participants were LTFU (20%). The median follow-up time for LTFU participants was 3 months. The attrition rate was 1.8/ 100 person-years. Newly enrolled participants left the cohort at a higher rate than former CohMSM participants (p-value:<0.001). Factors associated with LTFU can be found in Table 1 PE16.17

How men can support women in taking PrEP: perspectives from young women, male partners and peers in Siaya County, western Kenya 1) YW taking PrEP aged 20 to 24 who were enrolled in the DREAMS program, and 2) male peers/partners aged ≥ 18 years who were either a friend or sexual partner of a YW.

Participants completed several photo assignments focused on YW's PrEP adherence and persistence, including influential factors among males and how men could support YW. Photographs were presented and discussed in same-gendered groups using the SHOWeD method. Discussions were audio recorded, translated and transcribed into English, and analyzed using applied thematic analysis. Here we focus on overall themes that describe males' influences on YW's use of PrEP, identified across all photo assignments on adherence and persistence.

Results: Twenty-two YW (average age 22.5) and 17 men (average age 28.1) participated in the photo assignments. The majority were married and living with their partners (YW: n = 14; men: n = 10). YW's photographs and discussions typically depicted negative male influences on YW's PrEP use-i.e., men were more often viewed as barriers than as supporters. YW also described having little autonomy over their sexual lives, explaining that men tend to dictate when and how sexual encounters occur. Men's photographs and discussions suggested that men would support YW's PrEP use if PrEP was better promoted in the community and if men were more knowledgeable about PrEP in general. Men also explained that YW's PrEP use is hindered by stigmatizing peer and community attitudes.

Conclusions: Currently, men appear to hinder more than help YW's PrEP use. However, with greater male and community engagement about PrEP, men could play an important role in facilitating and normalizing PrEP use among YW.

PE16.18

Integrating oral pre-exposure prophylaxis services to public HIV care clinics in Kenya: Results from a pragmatic stepped-wedge randomized trial


Student Member, IEEE Francesco Lisi
email: francesco.lisi@phd.unipi.it

Stefano Fortunati
email: ste-fano.fortunati@centralesupelec.fr

Maria Sabrina Sabrina Greco
email: maria.greco@unipi.it

Fulvio Gini
email: fulvio.gini@unipi.it

Enhancement of a state-ofthe-art RL-based detection algorithm for Massive MIMO radars

Keywords: Adaptive selection, beamforming, constant false alarm rate, massive MIMO radar, Reinforcement Learning, SARSA, target detection

In the present work, a reinforcement learning (RL) based adaptive algorithm to optimise the transmit beampattern for a co-located massive MIMO radar is presented. Under the massive MIMO regime, a robust Wald-type detector, able to guarantee certain detection performances under a wide range of practical disturbance models, has been recently proposed. Furthermore, an RL/cognitive methodology has been exploited to improve the detection performance by learning and interacting with the surrounding unknown environment. Building upon previous findings, we develop here a fully adaptive and data-driven scheme for the selection of the hyper-parameters involved in the RL algorithm. Such an adaptive selection makes the Wald-RL-based detector independent of any ad-hoc, and potentially sub-optimal, manual tuning of the hyper-parameters. Simulation results show the effectiveness of the proposed scheme in harsh scenarios with strong clutter and low SNR values.

I. INTRODUCTION

The concept of Cognitive Radar (CR) has been firstly introduced by Haykin [START_REF] Haykin | Cognitive radar: a way of the future[END_REF]. The need of adaptation to the changes in the environment is intrinsic in the radar detection problem due to the presence of multiple sources of non-stationarity, such as variations in the clutter statistics over time or changes in the target scenario. Unlike communication systems where the transmitter and receiver are physically separated, in a monostatic radar system the transmitter and the receiver are in the same position, allowing the latter to easily transmit information to the former. MIMO radar systems can be divided into two fundamentals categories: widely separated radars and co-located radars. As [START_REF] Haimovich | MIMO radar with widely separated antennas[END_REF] suggests, a radar with widely separated antennas can exploit the spatial diversity of the target radar cross section to obtain a diversity gain similar to the one in communications. In a co-located MIMO radar, the antennas are closely spaced and each element of the array transmits a different probing signal, contrary to standard phased array where all the elements transmit the same waveform with variable amplitude and phase. In [START_REF] Li | MIMO radar with colocated antennas[END_REF] the authors discuss the advantages of parameter identifiability led by co-located MIMO radars. Following the recent works [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF], [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF], in this correspondence, we focus our attention on co-located massive MIMO radars exploiting RL-based techniques [START_REF] Sutton | Reinforcement learning: An introduction[END_REF], [START_REF] Mohri | Foundations of machine learning[END_REF] to implement the cognition loop [START_REF] Greco | Cognitive radars: On the road to reality: Progress thus far and possibilities for the future[END_REF]. Specifically, in [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF] the authors proved that, if the number of virtual spatial antenna channels is high enough (massive MIMO regime), a robust Waldtype test can be derived to guarantee the constant false alarm rate (CFAR) property under a wide variety of disturbance models using a single snapshot. Then, in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF], a reinforcement learning (RL) approach is proposed to optimise the transmission beampattern of a co-located massive MIMO radar exploiting the robust Wald-type detector in [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF].

Two original contributions are proposed in this correspondence. The first one is the introduction of two new policies, called quasi ε-greedy policy and quasi ε-greedy policy with target recovery. The former improves the performance of the system by forcing the RL algorithm to focus its power in a number of angular bins greater or equal to the number of detected targets; in addition, the latter exploits a mechanism that allows the system to recover a missed target more quickly. The second original contribution consists of an adaptive algorithm that selects the SARSA algorithm hyper-parameters (ε and α) based on the received signal. This allows the system to adapt to the variations of the surrounding environment and to achieve better detection capability both in stationary and non-stationary scenarios. Afterwards, the RL-based algorithm proposed in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] with the quasi ε-greedy policy with target recovery and the adaptive selection of ε and α is tested in various simulation scenarios. The results confirmed that the updated algorithm achieves higher detection performance, compared to the one in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF]. Both our algorithm and the one proposed in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] exploit the robust Wald-type detector described in [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF].

In Section II we describe the detection problem and the main properties of the robust Wald-type detector derived in [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF]. Section III provides the background on the RL-based detector proposed in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF]. In Section IV two variations of the ε-greedy policy are proposed, while an adaptive algorithm to select the ε and α parameters is described in Section V. In Section VI we validate the proposed contributions via simulation. The simulated scenarios are described in Appendix A. Finally, conclusions are drawn in Section VII.

NOTATIONS: In this paper we use upper case letters A and lower case ones a to denote matrices and vectors respectively. (•) T and (•) H denote a matrix transpose and conjugate transpose respectively, while (•) * denotes the conjugate operator. I N denotes the N × N identity matrix, while 0 N denotes an all zeros N × N matrix.

II. THE CFAR DETECTION ALGORITHM

Consider a co-located MIMO radar with N T transmit and N R receive antennas. Both transmit and receive arrays are uniform linear arrays (ULA) with inter-element distance d = λ/2. The transmitted signal vector can be expressed as s(t) = W Φ(t) where W ∈ C N T ×N T is a weighting matrix and Φ(t) ∈ C N T is a vector containing a set of orthonormal signals. After sampling the correct range-doppler bin at the output of the matched filter the received signal can be expressed as [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF]

y = αh + c ∈ C N ×1 , (1)
with

N = N T N R and h = W T a T (ν 0 ) ⊗ a R (ν 0 ).
Finally, a T (ν 0 ) and a R (ν 0 ) are the transmit and receive steering vectors that depend on the spatial frequency ν 0 = (d/λ) sin(θ 0 ) where θ 0 is the target angle of arrival (in a given reference frame).

Assuming that the received signal is processed at each time instant k by a bank of L spatial filters tuned to a fixed grid of angular bins Θ = {θ l } L l=1 , the detection problem can be formulated as the following hypothesis testing problem [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] H 0 : y k,l = c k,l ,

H 1 : y k,l = α k,l h k,l + c k,l . (2)
In the following, we assume that the grid is chosen in order to uniformly span the spatial frequency interval [-0.5, 0.5]. As a test statistic, the robust Wald-type detector Λ k,l is adopted [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF], [START_REF] Fortunati | Erratum to "Massive MIMO radar for target detection[END_REF]:

Λ k,l = 2|h H k,l y k,l | 2 h H k,l Γk,l h k,l , (3)
where Γk,l is an estimate of the covariance matrix of the noise vector c k,l . The detector discriminates between H 0 and H 1 by comparing the statistic with a threshold λ in each angular bin. When N → ∞, the statistic satisfies the CFAR property, i.e. the Probability of False Alarm (P F A ) is constant and can be selected by choosing the threshold as λ = -2 ln (P F A ), under a wide rage of practical disturbance processes with unknown statistical characterisation. For a P F A equal to 10 -4 , if N ≥ 10 4 the MIMO radar is assumed to operate in the massive MIMO regime and the previous property is satisfied [START_REF] Fortunati | Massive MIMO radar for target detection[END_REF].

III. RL-BASED ALGORITHM

After having ensured the CFARness of the detection scheme through the Wald-type detector in (3), the work [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] focused on how to exploit the large degrees of freedom, offered by a massive MIMO radar, to maximise the Probability of Detection (P D ). This can be achieved by properly selecting the weighting matrix W k to shape the transmit beampattern a T T (ν)W k W H k a * T (ν) [START_REF] Fuhrmann | Transmit beamforming for MIMO radar systems using signal cross-correlation[END_REF]. In order to do so, an RL-based algorithm has been proposed to allow the radar focusing the power in the angular bins associated to the targets' angular position [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF]. In the following, a concise summary of the main definitions and results obtained in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] will be provided since they represent the starting point of the original developments proposed in this correspondence. Let us start by introducing the Markov Decision Process (MDP) characterising the learning/detection task at hand. For the sake of clarity, it is worth mentioning that the following definitions of state, action and reward, associated to the considered MDP, are slightly different with respect to the ones introduced in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF]. Even though the essence remains exactly the same, this new formulation is more precise and self-consistent.

A. The set of the states

The state space of our MDP is denoted as

S ≜ ! s (i) " K i=0 , where K < ∞ is the maximum number of detectable target.
The state of the system at time instant k is defined as

s k = s (i k ) with [5] i k ≜ min # L-1 $ l=0 Λk,l , K % , (4)
where Λk,l ≜ u(Λ k,l -λ) and u(•) corresponds to the Heaviside step function1 .

B. The set of the actions

The action set of the MDP is indicated as

A ≜ ! a (j) " K j=0 .
If the action a (j k ) is selected at time instant k, then the beamforming algorithm focuses its power in the j k angular bins with the highest decision statistic.

Let ! l

(n) k " L n=1 be the set containing the angular bin's indexes sorted in descending order with respect to the decision statistic at time instant

k (Λ k,l (1) k ≥ ... ≥ Λ k,l (L) k
), we define the set

Ω k ≜ & ∅ , j k = 0 {l (n) k } j k n=1 , j k ∕ = 0 , (5)
containing the indexes of the angular bins associated to the j k highest decision statistics.

The system chooses the matrix W k according to [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF], [START_REF] Wang | Reinforcement learning-based waveform optimization for MIMO multi-target detection[END_REF]

W k = ' ( ) ( * W ort ≜ + Pmax N T • I N T , Ω k = ∅ # arg max W min l∈Ω k a T T (ν l )W W H a * T (ν l ) subject to tr ! W W H " ≤ Pmax , Ω k ∕ = ∅ (6)
where P max is the maximum transmitted power.

C. The reward function

Let's define the set

Ψ k ≜ ! l (n) k " K n=1 , where ! l (n) k " L n=1
is the same set as in Section B, and the sets

Φ k ≜ & ∅ , i k = 0 ! l (n) k " i k n=1 , i k ∕ = 0 , (7)
and Φk ≜ Ψ k \Φ k , where i k is defined in (4). The reward is defined as [5]

r k ≜ $ l∈Φ k PD,k,l - $ l∈ Φk PD,k,l (8)
where PD,k,l is an estimate of the P D of a target located in the l th angular bin at the k th time instant.

D. SARSA algorithm

The goal of RL algorithms is to find the optimal policy π * , i.e. the one that maximises the state value

function V π (s) ≜ E π # +∞ , h=0 γ h r k+h+1 - -s k = s % , ∀s ∈ S,
where γ ∈ [0, 1] is a damping factor. By defining the state-action value function associated to the policy π

as Q(s, a) ≜ E π # +∞ , h=0 γ h r k+h+1 - -(s k = s) ∩ (a k = a) % , it
can be proved that the greedy action associated to the state [START_REF] Sutton | Reinforcement learning: An introduction[END_REF], [START_REF] Mohri | Foundations of machine learning[END_REF]. The SARSA algorithm, named after the update rule sequence State-Action-Reward-State-Action, allows the system to recursively compute the Q matrix associated to the optimal policy when the dynamics of the environment are unknown.

s is equal to π * (s) = arg max a ′ ∈A Q π * (s, a ′ )
The algorithm proposed in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] consists of setting the initial values s 0 = s (0) , a 0 = a (0) and Q 0 = 0 K+1 and proceeding by computing the new state s k+1 and reward r k+1 , selecting a new action a k+1 following the ε-greedy policy and then updating Q as

Q k+1 (s k , a k ) = Q k (s k , a k )+ + α k+1 (r k+1 + γQ k (s k+1 , a k+1 ) -Q k (s k , a k )). (9)
Once the action has been selected the algorithm computes the W k+1 matrix following (6) and then the system transmits the new pulses. The learning rate α k ∈ (0, 1) is the weight given to the new information with respect to the old one.

E. ε-greedy policy

To ensure the convergence of the SARSA algorithm, the new action a k+1 must be selected according to a policy that guarantees that each state-action pair is visited infinitely many times. The ε-greedy policy is one of the most used ones in the RL-literature [START_REF] Sutton | Reinforcement learning: An introduction[END_REF], [START_REF] Mohri | Foundations of machine learning[END_REF]. The new action is selected as

π (1) k (s k ) = & a (greedy) k , w.p. 1 -ε U . A -{a (greedy) k } / , w.p. ε , (10)
where "w.p." stands for with probability, a

= arg max a ′ ∈A Q k (s k , a ′ ) (greedy) k
and U{C} denotes a function that selects randomly one of the elements in the set C with uniform probability. The ε parameter controls the exploration-exploitation tradeoff : low ε values correspond to a system that chooses the greedy action most of the time (exploitation), while high values to one that selects a random action more frequently (exploration).

IV. Policy improvement

This section and the following Section V present the original contribution of this correspondence. We start by providing two improved versions of the ε-greedy policy.

A. Quasi ε-greedy policy

The problem arising from choosing the ε-greedy policy in our specific application can be easily explained with an example.

Example 1. Consider a static scenario with two targets as in Scenario 1 (see Appendix A). Suppose that the system performs the optimal action at the k th iteration, i.e. selects action a k = a (2) and the set Ω k , defined in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF], contains the indices associated to the actual position of the two targets. If the system chooses the random action a k+1 = a (1) at the next time instant and focuses all the power in only one of the two targets, with high probability it will end up in s k+2 = s (1) . At this point, even if the system selects the action a k+2 = a (2) , the set Ω k+2 isn't guaranteed to contain the indexes of the two targets. The system will stay in state s (1) until the missed target will be in the set Ω k+m .

To overcome this issue a new policy is proposed here. We define the quasi ε-greedy policy as:

π (2) k (s k ) = & U ! A ′ (s k ) -{a (greedy) k } " , w.p. ε a (greedy) k , w.p. 1 -ε (11
) where A ′ (s (i) ) ≜ {a (j) , j = i, ..., K}. If the quasi εgreedy policy is combined with Q 0 = I K+1 then, when the system is in state s (i) , it can't focus its power in a number of angular bins less than i, that corresponds to the number of detected targets if these aren't more than K.

B. Quasi ε-greedy policy with target recovery

Although the quasi ε-greedy policy solved one of the main issues encountered with the ε-greedy policy, it doesn't solve the target loss problem described in the following example.

Example 2. Consider a scenario with two targets in the radar scene, like in Scenario 1 (see Appendix A). Suppose that, at a given time instant k, the system is in state s (2) , having detected both targets, and randomly selects action a [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] . Although two of the five angular bins in which the system focuses its power correspond to the position of the targets, the system might lose one of them if its SNR is low and end up in state s (1) . Now, with high probability, the two highest decision statistics still correspond to the angular position of the targets. Consequently, if the system selects action a (2) , it recovers the lost target and gets back in state s (2) . On the other hand, since the system is in state s (1) , it might happen that it selects action a (1) randomly or because this is the greedy action associated to state s [START_REF] Haykin | Cognitive radar: a way of the future[END_REF] , especially in the early phases when the system hasn't figured out the scene yet. In this case the system focuses all its power in the direction of the detected target and the decision statistic of the missed one drops. Now it will take a long time for the system to recover it.

To solve this "missed target" problem, we define a new policy as:

π (3) k (s k , s k-1 ) ≜ & arg max a∈A Q k (s k-1 , a) , i k < i k-1 π (2) k (s k ) , i k ≥ i k-1 .
(12) In words, this policy can be motivated as follows: when the system detects an higher or equal number of targets, it chooses the new action according to the quasi εgreedy policy. On the contrary, if the number of detected targets is smaller than the one at the previous iteration the system tries to recover them as soon as possible by choosing the greedy action associated to the state at the previous time instant s k-1 .

V. Adaptive selection of ε and α

The ε and α hyper-parameters are both kept constant in the algorithm proposed in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF]. Although this is one of the most popular choices for non-stationary environments, it has two major drawbacks. Firstly, in order to use the algorithm, the user must select them in advance and this operation requires a certain knowledge of both the problem and algorithm, thus limiting its usability. Secondly, the parameters' optimal value changes over time, both in stationary and non-stationary environments. To overcome these problems we propose an adaptive algorithm that selects both the parameters based on the collected data. Let r k be the reward at time instant k and d k the sequence defined as

d k ≜ & r 1 , k = 1 r k -r k-1 , k ∕ = 1 . ( 13
)
Then, ε and α can be updated, at the k th iteration, according to the following strategy:

x k+1 = ' ( ) ( * max{c 1 • x k , x min }, |d k | < η 1 min{c 2 • x k , x max }, η 1 < |d k | < η 2 x max , |d k | > η 2 (14)
where x corresponds to ε or α, c 1 ∈ (0, 1) and c 2 ∈ (1, +∞). The initial value of x is set to x 0 = x max . c 1 , c 2 , η 1 and η 2 are constants, but differ for the ε and α algorithm. Table I lists all the constants' values. x k+1 is not updated if the system was in exploration mode in the two previous time instants, i.e. k -1 and k. If the system selects a random action at time instant k, the reward r k+1 may drop due to a bad choice of the action causing |d k+1 | = |r k+1 -r k | to surpass η 1 even though the scenario hasn't changed. If the system then chooses the correct action at time k + 1, the reward r k+2 rises back to a value around r k , but |d k+2 | = |r k+2 -r k+1 | is likely to be over η 1 due to the low value of r k+1 .

The thresholds η 1 and η 2 depend on the definition of the reward. For the one defined in (8), if the system misses one target or detects a new one the reward variation will be around 1 depending on the target SNR. A good choice for η 1 is 0.5. Moreover, η 2 should be chosen high enough to guarantee that the value is reset only when a sudden change in the scenario happens: some empirical tests suggested that η 2 ≥ 1.8 meets this specification.

VI. SIMULATION RESULTS

In this section the performances of the two policies and the adaptive algorithm previously described are validated via simulations in Scenarios 1 and 2 (see Appendix A). Then the RL-based cognitive algorithm (RL-C) with TABLE II: Simulation parameters. L is the number of spatial filters in the receiver and K is the maximum number of detectable targets; γ is the damping factor of the cumulative reward and Q 0 is the initial value of the state-action value matrix.

P F A N T N R L K Pmax γ Q 0 10 -4 100 100 20 5 1 0.8 I K+1
the adaptive selection of ε and α and the quasi ε-greedy policy with target recovery is compared with the orthogonal algorithm and a "non RL-based" cognitive algorithm (NRL-C) in two non-stationary scenarios (3 and 4). The orthogonal algorithm is the one with an omnidirectional beampattern, while the NRL-C algorithm focuses the power in the angular bins where a detection occurred.

As performance benchmark, we plot an upper bound on the P D obtained by a clairvoyant beamformer that focuses its power in the exact (and generally unknown) direction of the targets. Finally, the proposed RL-C algorithm is compared with the one in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] in a stationary scenario. Due to a lack of space only the most relevant figures are shown, but the interested reader can obtain all the simulations' results using the code available via the link in the first page of the correspondence. Table II lists all the parameters' values used in the simulations.

Figure 1 shows the P D of target 2 in Scenario 1 for the ε-greedy policy and the two proposed policies. The figure shows that both the quasi ε-greedy policy and quasi ε-greedy policy with target recovery have better detection performances than the classical ε-greedy one. In particular, the quasi ε-greedy policy with target recovery has the best performances among the three policies.

Figure 2 compares the performance of the adaptive ε algorithm with the two static cases with ε = ε min and ε = ε max while keeping α constant (α = 0.5); Figure 3 does the same for α when ε is kept constant (ε = 0.5). All the graphs show that the adaptive algorithm combines the positive effects of having high values of both the parameters in the initial transitory phase, when the system has to gather information about the surrounding environment, and low values in the following phase when the system approaches a steady state.

Scenario 3 consists of two targets with fixed angular position and variable SNR, which is shown in the inset in Figure 4a. Figure 4 shows the P D of both targets. Since the SNR increases in the interval [START_REF] Haykin | Cognitive radar: a way of the future[END_REF]100] and decreases in [101,200], the P D varies accordingly. It is interesting to note that there is a delay between the drop of the P D and the time instant when the SNR starts to go down (k = 100) for both the RL and NRL cognitive algorithms. For example, considering the performance of the RL-C algorithm in Figure 4b, the P D goes down after k = 110; in the interval [101,110] the positive effect of the algorithm's learning capability prevails over the negative effect of the decreasing of the target's SNR. The RL-based algorithm shows better detection performance TABLE III: P D of the four targets in Scenario 3 in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] for a fixed P F A value of 10 -4 . The values in the second row are obtained using the proposed RL-C algorithm, while the ones in the third row are extracted from Figure 6 in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF].

Target 1 (ν = -0.2) 2 (ν = 0) 3 (ν = 0.2) 4 (ν = 0.3)
RL-C 1.00 0.98 0.99 0.97 [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] 1.00 0.74 0.91 0.73 than the non RL-based cognitive one, confirming that the system is able to exploit the information gathered from the surrounding environment. Since the SNR of both targets is low, the orthogonal algorithm can't detect them.

Even though the RL-based algorithm shows far better performance than the other algorithms, there's still a gap with the upper bound on the P D . Scenario 4 consists of three stationary targets. Target 1 and 2 are in the scene at time k = 1 and disappear at time k = 101 and k = 301 respectively. Target 3 appears at time k = 201 and stays in the scene until the end of the simulation. Figure 5 confirms that the RL-based algorithm is the one with the best detection performance among the three algorithms. Figure 6 and Figure 7 show the temporal evolution of the ε and α parameters: the adaptive algorithm is able to track the variations in the scene and adjust the parameters accordingly.

Finally, Table III compares the performance of the proposed RL-C algorithm with one obtained in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF]. Thanks to the developed quasi ε-greedy policy with target recovery and the adaptive selection of ε and α, a significant improvement has been achieved.

VII. CONCLUSION

In the present paper we introduced two variations of the ε-greedy policy. In addition, we proposed an adaptive algorithm to select the SARSA ε and α parameters, which increases both the RL-based algorithm performance and usability, since the user doesn't have to set them manually. Then the enhanced version of the algorithm proposed in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] was tested in various stationary and non-stationary scenarios. The results confirmed that the updated algorithm leads to a significant improvement of the detection capability of the massive MIMO radar system. On the other hand, the new policies and the adaptive selection of the hyper-parameters slightly increase the computational overhead, which must be taken into account when the algorithm is implemented by a real-time system. Furthermore, the time complexity of the algorithm grows exponentially with N T , due to the growth of the dimension of the weighting matrix W [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF], limiting its scalability. Future works will investigate the possibility to fill the gap between the performance of the improved RLbased algorithm and the upper bound. More specifically, we will try to fuse the position information (i.e. target tracking) with the RL-based detection algorithm.

In all the simulated scenarios the noise process is an AR(6) process with t-distributed independent and identically distributed (i.i.d.) innovations with the same parameters as the ones in [START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF].

Table IV contains all the details that characterise the four simulated scenarios.

  for the years 2015-2020, and AESS VP Publications (2018-2020) and EiC of the IEEE Aerospace and Electronic Systems Magazine. She's now IEEE SPS Director-at-Large for Region 8 and President-Elect of AESS (2022-23). Her general interests are in the areas of statistical signal processing, estimation and detection theory. In particular, her research interests include clutter models, coherent and incoherent detection in non-Gaussian clutter, CFAR techniques, radar waveform diversity, bistatic/mustistatic active and passive radars, cognitive radars and integration of sensing and communications. She co-authored many book chapters and more than 220 journal and conference papers. Fulvio Gini (Fellow IEEE) received the Doctor Engineer (cum laude) and the Research Doctor degrees in electronic engineering from the University of Pisa, Italy, in 1990 and 1995 respectively. In 1993 he joined the Department of Ingegneria dell'Informazione of the University of Pisa, where he become Associate Professor in 2000 and he is Full Professor since 2006. Prof. Gini is the Deputy Head of the Department since November 2016. From July 1996 through January 1997, he was a visiting researcher at the Department of Electrical Engineering, University of Virginia, Charlottesville. He is an Associate Editor for the IEEE Transactions on Aerospace and Electronic Systems (AES) since Jan. 2007 and for the Elsevier Signal Processing journal since Dec. 2006. He has been AE for the Transactions on Signal Processing (2000-06) and a Senior AE of the same Transaction (2016-17). He was a Member of the EURASIP JASP Editorial Board (2003-09). He was the Area Editor for the Special issues of the IEEE Signal Processing Magazine (2012-14). He was co-recipient of the 2001 and 2012 IEEE AES Society's Barry Carlton Award for Best Paper published in the IEEE Transactions on AES and co-recipient of the 2020 EURASIP JASP Best Paper Award. He was recipient of the 2022 IEEE AESS Warren White Award for Excellence in Radar Engineering, the 2020 EURASIP Meritorious Service Award, the 2003 IEE Achievement Award, and the 2003 IEEE AES Society Nathanson Award to the Young Engineer of the Year. He is the IEEE AES Society Awards Chair. He was Member of the IEEE AES Society Radar System Panel (2008-2020), member of the IEEE AES Society Board of Governors (BoG) (2017-2022) and member of the IEEE Signal Processing Society BoG (2021-2023). He was a member of the IEEE SPS Awards Board (2016-2018) and of the IEEE SPS Conference Board (2019-2020). He has been a Member of the Signal Processing Theory and Methods (SPTM) Technical Committee (TC) of the IEEE Signal Processing Society and of the Sensor Array and Multichannel (SAM) TC for many years. He is a member of the IEEE TAB Awards and Recognition Committee (TABARC) (2020-2021). He was a Member of the Board of Directors (BoD) of the EURASIP Society, the EURASIP Award Chair (2006-2012) and the EURASIP President (2013-2016). He was the General co-Chair of the 2020 IEEE Radar Conference, Florence (Italy). He was the Technical co-Chair of the 2006 EURASIP Signal and Image Processing Conference (EUSIPCO 2006), Florence (Italy), of the 2008 Radar Conference, Rome (Italy), and of the 2015 IEEE CAMSAP workshop, Cancun (Mexico). He was the General co-Chair of the 2nd Workshop on Cognitive Information Processing (CIP2010), of the 2014 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP 2014), and of the 2nd, 3rd and 4th editions of the workshop on Compressive Sensing in Radar (CoSeRa). He authored or co-authored 11 book chapters, 150 journal papers, 175 conference papers, and 2 national patents (under approval). H-index: 49, 9415 citations.

  quasi ε-greedy quasi ε-greedy with target recovery

Fig. 1 :Fig. 2 :Fig. 3 :

 123 Fig. 1: Policy comparison: P D of target 2 (Scenario 1).

Fig. 4 :

 4 Fig. 4: P D of the targets in Scenario 3. The inset figure in (a) shows the SNR of both targets expressed in dB.

Fig. 5 :

 5 Fig. 5: P D of the targets in Scenario 4.

Fig. 6 :Fig. 7 :

 67 Fig. 6: ε k sequence (Scenario 4)

  E{•} denotes the statistical expectation. The Kronecker product is represented by ⊗. A closed interval of numbers between a and b is denoted by [a, b], while a set containing only a and b is denoted by {a, b}. The absolute value is represented by |•|. The relative complement of set A with respect to set B is denoted as B\A.

TABLE I :

 I Adaptive ε and α algorithm parameters. ε min (α min ) and ε max (α max ) are the endpoints of the interval to which ε (α) belongs; c 1 and c 2 correspond to the multiplicative decrease and increase factors; η 1 and η 2 are the lower and upper threshold respectively.

	x	x min	xmax	c 1	c 2	η 1	η 2
	ε	0.1	0.8	0.8	2	0.5	1.8
	α	0.2	0.6	0.9	2.5	0.5	1.8

TABLE IV :

 IV Target scenarios.

	Scenario	Time Interval	Target	Angular Bin	ν	SN R dB
	1	[1,300]	1 2	7 16	-0.20 0.25	-20 -20
	2	[1,100]	1	17	0.30	-20
	3	[1,200]	1 2	7 16	-0.20 0.25	variable (Inset Fig.4a)
		[1,100]	1 2	5 13	-0.30 0.10	-18 -21
	4	[101,200] [201,300]	2 2 3	13 13 17	0.10 0.10 0.30	-21 -21 -20
		[301,400]

In this correspondence the notation (•) k,l is used, unlike the one (•) k l used in[START_REF] Ahmed | A reinforcement learning based approach for multitarget detection in Massive MIMO radar[END_REF] .

CORRESPONDENCE:

The work of Francesco Lisi was supported by the Laboratoire des signaux et systèmes (L2S). The work of M.S. Greco and F. Gini has been partially supported by the Italian Ministry of Education and Research (MIUR) in the framework of the CrossLab project (Departments of Excellence) of the University of Pisa, laboratory of Industrial Internet of Things (IIoT).