X Claeys
email: claeys@ann.jussieu.fr

Non-local variant of the Optimised Schwarz Method for arbitrary non-overlapping subdomain partitions

Non-local variant of the Optimised Schwarz Method for arbitrary non-overlapping subdomain partitions

Xavier

Introduction

The eective solution to large scale wave propagation problems relates to a wide range of applications and yet remains a challenge, in particular when simulating highly oscillatory phenomena. With the growing importance of parallel computing, an intense research eort has been dedicated, in recent years, to the development of domain decomposition strategies that can be eciently applied to wave propagation problems.

There is now a vast literature and a rich arsenal of well established domain decomposition techniques to deal with symmetric positive problems see e.g. [START_REF] Toselli | Domain decomposition methods algorithms and theory[END_REF][START_REF] Pechstein | Finite and boundary element tearing and interconnecting solvers for multiscale problems[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF]. By essence though, wave propagation does not fall into this symmetric positive framework and domain decomposition is much less developed for waves, from the point of view of both theory and eective numerical computation. In the case of harmonic regime propagation, the Optimized Schwarz Method (OSM) appears to be one of the most eective available approaches for domain decomposition in a wave context.

Ω 3 Ω 0 Ω 1 Ω 2 (a) Ω 3 Ω 0 Ω 1 Ω 2 (b) (c)
A general overview of this method and its numerous variants is given in [START_REF] Gander | A class of iterative solvers for the Helmholtz equation: factorizations, sweeping preconditioners, source transfer, single layer potentials, polarized traces, and optimized Schwarz methods[END_REF]. In OSM, the coupling of subdomains is maintained through transmission conditions at interfaces, and these transmission conditions are formulated in terms of ingoing and outgoing trace operators involving impedance coecients. The eciency of OSM crucially depends on the choice of these impedances.

The Optimized Schwarz Method, originally proposed by Lions [START_REF] Lions | On the Schwarz alternating method III. A variant for nonoverlapping subdomains[END_REF], was adapted by Després for wave propagation in [START_REF] Després | Décomposition de domaine et problème de Helmholtz[END_REF][START_REF] Després | Domain decomposition method and the Helmholtz problem[END_REF][START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle[END_REF][START_REF] Després | Domain decomposition method and the Helmholtz problem[END_REF] considering general non-overlapping partition of the computational domain and constant scalar impedance coecients. Although, in such a general geometrical setting, OSM with scalar impedance was proved to converge, no assessment was provided as regards the rate of convergence. In practice, the convergence could be slow. This was improved by Collino and Joly in [START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF][START_REF] Lecouvez | Quasi-local transmission conditions for non-overlapping domain decomposition methods for the Helmholtz equation[END_REF][START_REF] Claeys | Integral equation based optimized schwarz method for electromagnetics[END_REF] where the authors proposed operator valued self-adjoint positive impedance coecients and could establish geometric convergence of the method assuming that the subdomain partition does not involve any cross point i.e. point of adjacency of three interfaces (or one interface meeting the boundary of the compuational domain), see Fig. 1 above. In another series of contributions Antoine, Geuzaine and their collaborators [START_REF] Antoine | Optimized Schwarz domain decomposition methods for scalar and vector Helmholtz equations[END_REF][START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF][START_REF] Bouajaji | Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations[END_REF][START_REF] Antoine | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF][START_REF] Vion | Improved sweeping preconditioners for domain decomposition algorithms applied to time-harmonic Helmholtz and Maxwell problems[END_REF] considered the case of impedance coecients approaching appropriate Dirichlet-to-Neumann maps and obtained fastly converging numerical methods. Here also, the numerical methods were observed to be of good quality only when the subdomain partition does not contain any cross point.

To deal with cross points, the literature dedicated to substructuring DDM already oers techniques such as dual-primal FETI (see e.g. [START_REF] Pechstein | Finite and boundary element tearing and interconnecting solvers for multiscale problems[END_REF]Chap.5]) which are purely discrete methods where so-called primal dofs located at cross points are kept unchanged, while the other dofs are teared apart in substructures. Each iteration of the DDM solver then requires inversion of the matrix coupling all primal dofs. Dual-primal FETI was adapted to wave propagation by Boubendir [START_REF] Boubendir | Techniques de Décomposition de Domaine et Méthodes d'Equations Intégrales[END_REF][START_REF] Bendali | Non-overlapping domain decomposition method for a nodal nite element method[END_REF] by further augmenting local matrices to enforce the matching conditions in a way similar to the approach of Després. Energy estimates were established and led to a convergence result of DDM solvers although Section 4 in [START_REF] Bendali | Non-overlapping domain decomposition method for a nodal nite element method[END_REF] pointed that the convergence rate could not be proved mesh-uniform.

Although Boubendir's work already provides an eective treatment, cross-points remained a thorny issue (see [START_REF] Gander | On the applicability of Lions' energy estimates in the analysis of discrete optimized Schwarz methods with cross points[END_REF] for a thorough explanation) which recently attracted a renewed attention [START_REF] Antoine | Modave Corner treatment for high-order local absorbing boundary conditions in high-frequency acoustic scattering[END_REF][START_REF] Antoine | A non-overlapping domain decomposition method with high-order transmission conditions and cross-point treatment for Helmholtz problems[END_REF][START_REF] Gander | Cross-points in domain decomposition methods with a nite element discretization[END_REF][START_REF] Després | Corners and stable optimized domain decomposition methods for the Helmholtz problem[END_REF]. Very similar diculties arise in a dierent context: the derivation of Boundary Integral Equations (BIE) adapted to multi-domain scattering. The Multi-Trace Formalism (MTF) was introduced in [START_REF] Claeys | Multi-trace boundary integral formulation for acoustic scattering by composite structures[END_REF][START_REF] Claeys | Quasi-local multitrace boundary integral formulations[END_REF][START_REF] Claeys | Multitrace boundary integral equations[END_REF][START_REF] Claeys | Electromagnetic scattering at composite objects: a novel multi-trace boundary integral formulation[END_REF] as a complete framework for dealing with multi-domain BIE. From the perspective of functional analysis, MTF oers a clean treatment of cross-points. It would thus appear natural to try using the techniques developed in the Multi-Trace framework for dealing properly with cross points in Optimized Schwarz domain decomposition. This is precisely the aim of the present contribution.

In the present article, we rigorously establish equivalence bewteen a classical scattering problem in heterogeneous medium and a novel formulation [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] posed on the skeleton of a subdomain partition. This formulation is proved strongly coercive. We interpret Richardson's linear solver applied to this novel skeleton formulation as a new variant of the Optimized Schwarz Method. This new variant can be applied with any non-overlapping partition of the propagtion medium into Lipschitz subdomains, no matter the presence of cross-points. The key ingredient in this formulation is a non-local exchange operator used to enforce transmission conditions.

Such exchange operator has always existed in previous versions of OSM, but it was so far systematically assumed to be a local operator consisting in swapping the traces from both sides of each interface of the subdomain partition. The exchange operator we consider here is more elaborate, which is the main novelty of our approach. The regularity assumptions that we formulate are rather loose regarding both material coecients (piecewise Lipschitz) and geometry of subdomains (Lipschitz), which is another novelty made possible by our special exchange operator.

The formulation described in the present contribution cannot be considered as an extension of other pre-existing OSM strategies including [START_REF] Boubendir | Techniques de Décomposition de Domaine et Méthodes d'Equations Intégrales[END_REF][START_REF] Bendali | Non-overlapping domain decomposition method for a nodal nite element method[END_REF][START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF][START_REF] Lecouvez | Quasi-local transmission conditions for non-overlapping domain decomposition methods for the Helmholtz equation[END_REF][START_REF] Claeys | Integral equation based optimized schwarz method for electromagnetics[END_REF]. The dierence lies in our new exchange operator that is non-local and does not coincide with the standard local operator that swaps traces from each side of each interface. Because of this non-locality, even subdomains that are not adjacent will be coupled, contrary to other OSM strategies.

It should be mentionned that the present contribution is purely analytical and that, in its present form, this new variant of OSM does not seem appropriate for actual numerical computations. This is why we do not report on numerical results. A discrete counterpart of the present strategy is presented in [START_REF] Claeys | Robust treatment of cross points in Optimized Schwarz Methods available as[END_REF] where concrete numerical results are available. We still believe that the formulation we present here is an interesting theoretical object. In particular, it yields a strongly coercive formulation of Helmholtz problem which is not trivial: the derivation of coercive formulations for Helmoltz equation has been, in itself, the subject of recent attention [START_REF] Moiola | Is the Helmholtz equation really sign-indenite?[END_REF].

In the case of piecewise constant material coecients, Formulation [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] can also be used as a multi-domain coupling scheme for the solution to scattering problems by means of boundary integral formulation. In the particular case of piecewise constant coecients, the new formulation presented here can be considered as an alternative to other multi-domain BIE such as Multi-Trace [START_REF] Claeys | Multitrace boundary integral equations[END_REF], Boundary Element Tearing and Interconnecting [START_REF] Langer | Boundary element tearing and interconnecting methods[END_REF], or Rumsey's reaction principle [START_REF] Petersdor | Boundary integral equations for mixed Dirichlet, Neumann and transmission problems[END_REF]. However we emphasize that the present contribution is not specically oriented toward the derivation of a coupled boundary integral equation system, and a salient feature of the subsequent analysis lies in its ability to deal with heterogeneous propagation media.

Geometry and problem under study

In the present article, we are interested in a classical wave propagation problem in harmonic regime set in an heterogeneous medium in R d for d = 1, 2 or 3. We consider two essentially bounded measurable functions µ : R d → R + and κ : R d → C + , and we assume that there exist constants κ 0 , ρ 0 > 0 such that i)

sup x∈R d (|µ(x)| + |µ -1 (x)| + |κ(x)|) < +∞ ii) e{κ(x)} ≥ 0, m{κ(x)} ≥ 0, κ(x) = 0 ∀x ∈ R d iii) κ(x) = κ 0 and µ(x) = 1 for |x| > ρ 0 (1)
These assumptions are rather general yet reasonable enough to make the scattering problem we wish to examine properly well posed. We insist that we do not assume κ, µ to be piecewise constant. For some f ∈ L 2 (R d ) with bounded support, we wish to solve the following wave propagation problem in heterogeneous medium

           u ∈ H 1 loc (R d ) such that -div(µ∇u) -κ 2 u = f in R d , lim ρ→∞ ∂Bρ
|∂ ρ u -ıκ 0 u| 2 dσ ρ = 0.

(

) 2
where B ρ refers to the ball of radius ρ centered at 0, σ ρ is the associated surface measure, and ∂ ρ is the partial derivative with respect to |x|.

R d = ∪ J j=0 Ω j with Ω j ∩ Ω k = ∅ if j = k, each Ω j
is a Lipschitz domain, and Ω j is bounded for j = 0. The "skeleton" will refer to the union of all interfaces between subdomains

Γ = ∂Ω 0 ∪ • • • ∪ ∂Ω J .
We emphasize that such geometrical conguration allows the presence of junction points i.e. points where three subdomains or more abut. Examples of such non-overlapping multi-domain congurations are given in Fig. 1.

For the sake of simplicity, we make further regularity assumptions on material coecients in each subdomain, assuming that µ is Lipschitz regular in each subdomain, ∇µ j ∈ L ∞ (Ω j ) ∀j = 0 . . . J, where µ j := µ| Ω j .

(3) Assumptions (1)-(3) allow the coecients µ, κ to jump across the interfaces ∂Ω j ∩ ∂Ω k , but discards jumps of µ inside each subdomain. In particular, this setting includes the case where µ, κ are piecewise constant with respect to the subdomain partition. Problem (2) can be decomposed according to the subdomain partition introduced above, leading to wave equations in each subdomain coupled by transmission conditions imposed through each interface

           u ∈ H 1 loc (Ω j ) such that -div(µ∇u) -κ 2 u = f in Ω j , lim ρ→∞ ∂Bρ |∂ ρ u -ıκ 0 u| 2 dσ ρ = 0, (4)
   u| int ∂Ω j -u| int ∂Ω k = 0 ∀j, k = 0 . . . n µ j ∂ n j u| int ∂Ω j + µ k ∂ n k u| int ∂Ω k = 0 on ∂Ω j ∩ ∂Ω k . ( 5
)
where n j refers to the normal vector eld on ∂Ω j directed toward the exterior of Ω j , and ∂ n j v := n j • ∇v. The boundary traces coming into play in the transmission conditions above are taken from the interior of the subdomains, which is the meaning of the "int" superscript.

The present contribution will consist in deriving a strongly coercive reformulation of Problem (4)-( 5) of the form "identity+contraction". This reformulation will be posed in a space of trace on the skeleton Γ.

From here until the end of Section 5, we are going to develop a functional toolkit for traces on the skeleton of our subdomain partition. Besides topological considerations (spaces, norms,. . . ), this framework involves two important ingredients: a global DtN map T (see 3.4) and a non-local exchange operator Π (see Corollary 5.1). This framework is based on the Yukawa equation i.e. a PDE associated to the operator -∆ + γ -2 for some parameter γ > 0. We emphasize that the analysis presented in Section 3-5 is independent of the scattering problem [START_REF] Allaire | Numerical linear algebra[END_REF] and that the parameter γ is not a priori connected to material coecients µ, κ.

We will come back to our wave propagation problem in Section 6-7 where, relying on this functional toolkit, we shall derive a reformulation of ( 2) and establish its well-posedness and strong coercivity.

Trace spaces and operators

The treatment of interfaces between subdomains is a crucial aspect of any domain decomposition strategy, both for constructing or analysing it. As a consequence we pay a special attention to trace spaces.

Volume based spaces

First of all we need to x a few notations related to classical volume based function spaces.

For any Lipschitz domain Ω ⊂ R d , the space L 2 (Ω) will refer to square integrable functions equipped with the norm

ϕ 2 L 2 (Ω) := Ω |ϕ| 2 dx. The Sobolev space H 1 (Ω) := {ϕ ∈ L 2 (Ω), ∇ϕ ∈ L 2 (Ω) d } will be equipped with the norm v 2 H 1 (Ω) := ∇v 2 L 2 (Ω) + γ -2 v 2 L 2 (Ω) (6)
In this denition γ > 0 refers to a parameter that will be xed all through this article.

Occasionally we shall consider H(div, Ω) :

= {ψ ∈ L 2 (Ω) d , div(ψ) ∈ L 2 (Ω)} and H 1 (∆, Ω) := {ϕ ∈ H 1 (Ω), ∆ϕ ∈ L 2 (Ω)} equipped with the norm given by ϕ 2 H 1 (∆,Ω) := ϕ 2 H 1 (Ω) + ∆ϕ 2 L 2 (Ω) .
Finally if H(Ω) refers to any of the spaces introduced above, then

H loc (Ω) shall refer to all functions v : Ω → C such that vϕ ∈ H(Ω) for all ϕ ∈ C ∞ comp (R d ) := {ψ ∈ C ∞ (R d ), supp(ψ) bounded}.

Traces on the boundary of a single subdomain

Let us consider a Lipschitz open set Ω ⊂ R d such that either Ω or R d \ Ω is bounded. We shall refer to the space of Dirichlet traces H 1/2 (∂Ω) := {v| ∂Ω , v ∈ H 1 (Ω)} equipped with the norm, see [START_REF] Mclean | Strongly elliptic systems and boundary integral equations[END_REF]Thm.3.40],

v H 1/2 (∂Ω) := min{ ϕ H 1 (Ω) , ϕ| ∂Ω = v}. (7)
The space of Neumann traces H -1/2 (∂Ω) will be dened as the dual to H 1/2 (∂Ω) equipped with the corresponding canonical dual norm p

H -1/2 (∂Ω) := sup v∈H 1/2 (∂Ω) | p, v ∂Ω |/ v H 1/2 (∂Ω) .
Here v → p, v ∂Ω := p(v) simply refers to the action of p on v, so that (p, v) → p, v ∂Ω is a bilinear (not sesquilinear) form. As regards duality pairing, we shall also equivalently write v, p ∂Ω := p, v ∂Ω and ∂Ω pvdσ = p, v ∂Ω .

We will also equip the space of pairs of Dirichlet/Neumann traces with its own duality pairing.

Although many choices are possible, we use a skew-symmetric pairing that appears naturally in energy conservation calculus, dened by

[(u, p), (v, q)] ∂Ω := u, q ∂Ω -v, p ∂Ω for (u, p) and (v, q) in H + 1 2 (∂Ω) × H -1 2 (∂Ω). (8)
Note that this pairing does not involve any complex conjugation. Let n Ω refer to the normal vector eld on ∂Ω directed toward the exterior of Ω.

d : H 1 loc (Ω) → H 1/2 (∂Ω), τ Ω n : H 1 loc (∆, Ω) → H -1/2 (∂Ω) and τ Ω : H 1 loc (∆, Ω) → H 1/2 (∂Ω) × H -1/2 (∂Ω) uniquely dened by τ Ω d (ϕ) := ϕ| ∂Ω and τ Ω n (ϕ) := n Ω • ∇ϕ| ∂Ω , τ Ω (ϕ) := (τ Ω d (ϕ), τ Ω n (ϕ)) ∀ϕ ∈ C ∞ (Ω). (9)
Remark 3.1. In the case where Ω is one of the subdomains Ω j , j = 0 . . . J, assume that some

function v ∈ H 1 loc (Ω j ) satises div(µ∇v) ∈ L 2 loc (Ω j ).
Then according to [START_REF] Antoine | A quasi-optimal non-overlapping domain decomposition algorithm for the Helmholtz equation[END_REF] and ( 3), we have ∆v + µ -1 j ∇µ j • ∇v ∈ L 2 loc (Ω j ) and thus ∆v ∈ L 2 loc (Ω j ). In particular if u ∈ H 1 loc (Ω j ) satises div(µ∇u) + κ 2 u = f in Ω j with f as above, then u ∈ H 1 loc (∆, Ω j ) and its Dirichlet and Neumann traces in [START_REF] Börm | Ecient numerical methods for non-local operators[END_REF] are properly dened. Remark 3.2. In the case where Ω is the exterior subdomain Ω 0 , according to previous paragraphs, the normal vector n 0 is directed toward the exterior of Ω 0 and the traces in [START_REF] Börm | Ecient numerical methods for non-local operators[END_REF] are taken from the interior of Ω 0 .

Scalar products and Dirichlet-to-Neumann maps

For any v ∈ H 1/2 (∂Ω) let φ d (v) ∈ H 1 (Ω) refer to the unique element that achieves the minimum in [START_REF] Boubendir | Techniques de Décomposition de Domaine et Méthodes d'Equations Intégrales[END_REF] i.e. such that v H 1/2 (∂Ω) = φ d (v) H 1 (Ω) . Writing Euler's identity for this minimisation problem, we see that

Ω ∇φ d (v) • ∇ϕ + γ -2 φ d (v)ϕ dx = 0 ∀ϕ ∈ H 1 0 (Ω), which re-writes -∆φ d (v) + γ -2 φ d (v) = 0 in Ω. Then we introduce a so-called Dirichlet-to-Neumann (DtN) map T Ω := τ Ω n • φ d : H 1/2 (∂Ω) → H -1/2 (
T Ω (v) := n Ω • ∇φ d (v)| ∂Ω where φ d (v) ∈ H 1 (Ω) satises ∆φ d (v) -γ -2 φ d (v) = 0 in Ω φ d (v)| ∂Ω = v on ∂Ω. ( 10
)
This DtN map actually induces the scalar product associated to the norm [START_REF] Boubendir | Techniques de Décomposition de Domaine et Méthodes d'Equations Intégrales[END_REF]. First of all observe that φ d (u) = φ d (u) obviously, which implies T Ω (v) = T Ω (v). Next, according to the PDE satised by φ d in [START_REF] Claeys | A single trace integral formulation of the second kind for acoustic scattering[END_REF], applying Green's formula we obtain

Ω ∇φ d (u) • ∇φ d (v) + γ -2 φ d (u)φ d (v)dx = ∂Ω φ d (v) n Ω • ∇φ d (u)dσ = T Ω (u), v ∂Ω .
From this calculation it is clear that T Ω (u), v ∂Ω = T Ω (v), u ∂Ω in other words T Ω is a self-adjoint operator. This property plays an important role in the denition of the scalar product we wish to adopt for the trace space H 1/2 (∂Ω). Since by the dention of φ d we have u H 1/2 (∂Ω) = φ d (u) H 1 (Ω) , we can take the following as scalar product on the Dirichlet trace spaces

(u, v) H 1/2 (∂Ω) := T Ω (u), v ∂Ω for u, v ∈ H 1/2 (∂Ω). (11)
According to Riesz representation theorem, for any p ∈ H -1/2 (∂Ω) there exists a unique

ϕ p ∈ H 1/2 (∂Ω) such that p, v ∂Ω = (ϕ p , v) H 1/2 (∂Ω) = T Ω (ϕ p ), v for all v ∈ H 1/2 (∂Ω). Hence ϕ p = (T Ω ) -1 (p) and p 2 H -1/2 (∂Ω) = ϕ p 2 H 1/2 (∂Ω) = T Ω (ϕ p ), ϕ p ∂Ω = p, T -1 Ω (p) ∂Ω .
As a consequence the norm on Neumann data is induced by the following scalar product

(p, q) H -1/2 (∂Ω) := p, T -1 Ω (q) ∂Ω for p, q ∈ H -1/2 (∂Ω). (12)

Traces in a multi-domain setting

We will also need to consider cartesian products of Dirichlet or Neumann trace spaces based on the boundary of each subdomain of the partition, which we call multi-trace spaces dened as follows

H d (Γ) := H + 1 2 (∂Ω 0 ) × • • • × H + 1 2 (∂Ω J ), H n (Γ) := H -1 2 (∂Ω 0 ) × • • • × H -1 2 (∂Ω J ), H(Γ) d := Π j=0...J H + 1 2 (∂Ω j ) × H -1 2 (∂Ω j ). (13)
equipped with p 2 Hn(Γ)

:= p 0 2 H -1/2 (∂Ω 0 ) + • • • + p J 2 H -1/2 (∂Ω J ) for p = (p j ) J j=0 ∈ H n (Γ)
u, p := J j=0 u j , p j ∂Ω j ∀u = (u 0 , . . . , u J ) ∈ H d (Γ), ∀p = (p 0 , . . . , p J ) ∈ H n (Γ). (14)
The bilinear pairing dened above does not involve any complex conjugation operation. We shall indierently write p, u

:= u, p for u ∈ H d (Γ), p ∈ H n (Γ).
For the sake of conciseness, we shall denote T j instead of T Ω j . The operator T := diag j=0...J (T j ) :

H d (Γ) → H n (Γ) induces a scalar product underlying the norm of H n (Γ) through (p, q) Hn(Γ) = T -1 (p), q = J j=0 T -1 j (p j ), q j ∂Ω j (p, q) Hd(Γ) = J j=0 (p j , q j ) H -1/2 (∂Ω j ) . ( 15
)
for any p = (p j ) j=0...J and any q = (q j ) j=0...J in H n (Γ). As regards H(Γ), we shall consider a duality pairing given by the following skew symetric bilinear form

u, v := [u 0 , v 0 ] ∂Ω 0 + • • • + [u J , v J ] ∂Ω J for u = (u j ) J j=0 and v = (v j ) J j=0 in H(Γ). (16)
As regards trace operators, for the sake of conciseness, we shall denote τ j := τ Ω j and adopt similar conventions for τ j d and τ j n . We also introduce global trace operators that map into multi-trace spaces τ α (u) := (τ 0 α (u), . . . , τ J α (u)) for α = d, n τ (u) α := (τ 0 (u), . . . , τ J (u)).

(

) 17
4 Transmission conditions

Since we are considering a problem involving transmission conditions [START_REF] Antoine | A non-overlapping domain decomposition method with high-order transmission conditions and cross-point treatment for Helmholtz problems[END_REF], it is natural to introduce the subspace of H(Γ) consisting in all tuples of traces agreeing with these conditions:

this is what shall be called single-trace spaces dened by

X d (Γ) := { (v j ) J j=0 ∈ H d (Γ) | ∃ϕ ∈ H 1 (R d ), v j = ϕ| ∂Ω j ∀j } X n (Γ) := { (q j ) J j=0 ∈ H n (Γ) | ∃ψ ∈ H(div, R d ), q j = n j • ψ| ∂Ω j ∀j } X(Γ) d := { u = (u d , u n ) ∈ H(Γ) | u d ∈ X d (Γ), u n ∈ X n (Γ) } (18)
By construction, for a function u ∈ L 2 loc (R d ) such that u| Ω j ∈ H 1 loc (∆, Ω j ) for all j = 0 . . . J, the transmission conditions ( 5) are equivalent to the statement τ (u) ∈ X(Γ). The singletrace space has been extensively studied in the context of multi-trace formulations [START_REF] Claeys | Multi-trace boundary integral formulation for acoustic scattering by composite structures[END_REF]. The following caracterisation of this space was proved in [13, Prop.6.3]. For any u ∈ H(Γ) we have u ∈ X(Γ) ⇐⇒ u, v = 0 ∀v ∈ X(Γ).

Proof:

From [START_REF] Claeys | Second-kind boundary integral equations for electromagnetic scattering at composite objects[END_REF], it is clear that any u = (u d , u n ) ∈ H(Γ) actually belongs to X(Γ) if and only if u d ∈ X d (Γ) and u n ∈ X n (Γ). As a consequence, to prove the lemma, it suces to show that for any u d ∈ H d (Γ) and any u n ∈ H n (Γ)

we have i) u d ∈ X d (Γ) ⇐⇒ u d , q = 0 ∀q ∈ X n (Γ) ii) u n ∈ X n (Γ) ⇐⇒ u n , v = 0 ∀v ∈ X d (Γ)
We will only present the proof of i) since the proof for ii) is very similar. Take an arbitrary

u d = (u j d ) J j=0 ∈ H d (Γ). If u d ∈ X d (Γ), there exists ϕ ∈ H 1 (R d
) such that ϕ| ∂Ω j = u j d ∀j = 0 . . . J. Then for any q = (q j ) J j=0 ∈ X n (Γ), there exists ψ ∈ H(div, R d ) such that n j • ψ| ∂Ω j = q j ∀j = 0 . . . J. Applying a Green formula in each Ω j on the one hand, and in R d on the other hand, we obtain

u d , q = J j=0 u j d , q j ∂Ω j = J j=0 ∂Ω j n j • ψϕ dσ = J j=0 Ω j ∇ϕ • ψ + ϕ divψ dx = R d ∇ϕ • ψ + ϕ divψ dx = 0. ( 19
) Now assume that u d = (u 0 d , . . . , u J d ) ∈ H d (Γ) satises u d , q = 0∀q ∈ X n (Γ). For each j = 0 . . . J, introduce a lifting v j ∈ H 1 (Ω j ) such that v j | ∂Ω j = u j d , and set v(x) = 1 Ω 0 (x)v 0 (x)+ • • • + 1 Ω J (x)v J (x). We have clearly v ∈ L 2 (R d ) and, to prove that u d ∈ X d (Γ), it suces to show that v ∈ H 1 (R d ). Dene p ∈ L 2 (R d ) by p(x) = 1 Ω 0 (x)∇v 0 (x) + • • • + 1 Ω J (x)∇v J (x).
Pick an arbitrary ψ ∈ H(div, R d ), and set q = (q j ) J j=0 where q j := n j •ψ| ∂Ω j . Since q ∈ X n (Γ),

we have

R d v div(ψ)dx = J j=0 Ω j v div(ψ)dx = u d , q -J j=0 Ω j ψ • ∇v j dx = - R d ψ • pdx (20)
Since the above identity holds for any ψ ∈ H(div, R d ), we conclude that v admits a weak gradient over R d as a whole with p = ∇v in R d and, as a consequence v ∈ H 1 (R d ) and u d ∈ X(Γ).

As underlined during its proof, the above caracterisation implies that u ∈ H d (Γ) We have the direct sum

H n (Γ) = X n (Γ) ⊕ T(X d (Γ))
and it is orthogonal with respect to the scalar product induced by T -1 .

Proof:

First, according to Proposition 4.1, we have (p, T(u))

Hn(Γ) = p, u = 0 whenever p ∈ X n (Γ) and u ∈ X d (Γ). This proves that X n (Γ) is orthogonal to T(X d (Γ)) hence X n (Γ) ∩ T(X d (Γ)) = {0}.
Next pick an arbitrary p ∈ H n (Γ) and, by Riesz representation theorem, dene u as the unique element of X d (Γ) satisfying T(u), v = p, v for all v ∈ X d (Γ). As a consequence q = p-T(u) satises q, v = 0∀v ∈ X d (Γ) and thus belongs to X n (Γ) according to Proposition 4.1. This shows that H n (Γ) = X n (Γ) + T(X d (Γ)).

Potential theory

The problem (2) primarily considered in the present manuscript does not a priori lend itself to boundary integral equation techniques simply because (2) is a problem of propagation in heterogeneous media i.e. the PDEs involve a priori varying coecients. However several aspects of the solution strategy we wish to describe involve nonlocal operators. In particular, we shall need such theoretical tools for treatment of junctions. As a consequence, we dedicate the present section to recalling a few facts about boundary integral operators.

Layer potentials in a single subdomain

We rst introduce the Green kernel G (x) of the Yukawa's equation i.e. we dene G as the unique function solving -∆G + γ -2 G = δ 0 in R d and lim |x|→∞ G (x) = 0, where δ 0 is the Dirac measure centered at x = 0, and γ > 0 is a parameter that we have xed once and for all in 3.1. This kernel admits an explicit expression in terms of special functions namely

G (x) := K 0 (|x|/γ), x ∈ R 2 \ {0} for d = 2, G (x) := exp(-|x|/γ) 4π|x| , x ∈ R 3 \ {0} for d = 3. ( 21
)
where K 0 refers to the modied Bessel function of the second kind of order 0 also known as MacDonald function, see [47, 10.25]. With this kernel, and for any Lipschitz domain Ω ⊂ R d with bounded boundary, we can dene single and double layer potentials as follows: for any

(v, q) ∈ H 1/2 (∂Ω) × H -1/2 (∂Ω) we set Ψ Ω (v, q)(x) := Ψ Ω d (v)(x) + Ψ Ω n (q)(x), where Ψ Ω d (v)(x) := ∂Ω n Ω (y) • (∇G )(x -y)v(y)dσ(y),
where Ψ Ω n (q)(x) := ∂Ω G (x -y)q(y)dσ(y),

for all x ∈ R d \ ∂Ω.

For any v ∈ H 1/2 (∂Ω) × H -1/2 (∂Ω), we have (γ -2 -∆)Ψ Ω (v) = 0 both in Ω and R d \ Ω. Besides Ψ Ω (v)| Ω ∈ H 1 (∆, O) for O = Ω or O = R d \ Ω. For any x, y ∈ R d , x = y, dene G x : R d \ {x} → R + by G x (y) := G (x -y). Elementary calculus shows that Ψ Ω (u)(x) = [τ Ω (G x ), u] ∂Ω for all u ∈ H 1/2 (∂Ω) × H -1/2 (
satisfying (γ -2 -∆)u = 0 in Ω, we have Ψ Ω (τ Ω (u)) = 1 Ω (x)u(x) ∀x ∈ R d .
Here 1 Ω (x) = 1 if x ∈ Ω and 1 Ω (x) = 0 otherwise. In the representation formula above, the traces of solutions to the homogeneous PDE play a pivotal role. The potential operators actually provide a Calderón projector that maps onto such a space and can thus be used to caracterise them. A proof of the next result can be found in [ The operator τ Ω • Ψ Ω : H 1/2 (∂Ω) × H -1/2 (∂Ω) → H 1/2 (∂Ω) × H -1/2 (∂Ω) is a continuous projector whose range is the space C in (Ω) := {τ Ω (u) | u ∈ H 1 (Ω), (γ -2 -∆)u = 0 in Ω }.

Layer potentials in a multi-domain setting

We now establish a few results about potential theory that are specic to the multi-domain context. Part of the present section is a variant of results already available in [START_REF] Claeys | A single trace integral formulation of the second kind for acoustic scattering[END_REF][START_REF] Claeys | Quasi-local multitrace boundary integral formulations[END_REF]. However there are novelties. In particular we consider here a strongly coercive equation as opposed to the heterogeneous indenite problem of [START_REF] Claeys | A single trace integral formulation of the second kind for acoustic scattering[END_REF], which leads to stronger results and a completely new formulation of transmission conditions, see Proposition 5.4.

Considering Ω = Ω j for j = 0 . . . J, the result of the previous paragraph can be used directly in the multi-domain context. For the sake of conciseness, in the following, we shall

write Ψ j d , Ψ j n , Ψ j instead of Ψ Ω j d , Ψ Ω j
n , Ψ Ω j . We rst show that an explicitl formula for the orthogonal projector onto X n (Γ), can be obtained. We rely on so-called multi-potential operators Ψ d :

H d (Γ) → Π J j=0 H 1 loc (∆, Ω j ) and Ψ n : H n (Γ) → Π J j=0 H 1 loc (∆, Ω j ) dened as follows: for any u = (u d , u n ) ∈ H(Γ) we set Ψ(u)(x) = Ψ d (u d )(x) + Ψ n (u n )(x) where Ψ d (u d )(x) := J j=0 Ψ j d (u j d )(x), where Ψ n (u n )(x) := J j=0 Ψ j n (u j n )(x). ( 23
)
for any x ∈ R d \ Γ. Such operators have been rst considered in the context of the integral formulation of the second kind introduced in [START_REF] Claeys | A single trace integral formulation of the second kind for acoustic scattering[END_REF], see also [START_REF] Greengard | Stable and accurate integral equation methods for scattering problems with multiple material interfaces in two dimensions[END_REF][START_REF] Peng | Computations of electromagnetic wave scattering from penetrable composite targets using a surface integral equation method with multiple traces[END_REF][START_REF] Peng | A boundary integral equation domain decomposition method for electromagnetic scattering from large and deep cavities[END_REF][START_REF] Claeys | A second-kind Galerkin boundary element method for scattering at composite objects[END_REF][START_REF] Claeys | Second-kind boundary integral equations for electromagnetic scattering at composite objects[END_REF][START_REF] Claeys | Second kind boundary integral equation for multi-subdomain diusion problems[END_REF][START_REF] Claeys | Quasi-local multitrace boundary integral formulations[END_REF]. The multipotential operators satisfy many non-trivial properties. To begin with, the next proposition shows that they are closely related to global Dirichlet-to-Neumann maps. For

any u = (u d , u n ) ∈ H(Γ) we have τ • Ψ(u) = u ⇐⇒ u n = T(u d ).
Proof:

Pick a u = (u d , u n ) = (u j d , u j n ) J j=0 ∈ H(Γ) with u n = T(u d ). We have u j := (u j d , u j n ) = (u j d , T j (u j d )) ∈ C in (Ω j )
for each j = 0 . . . J. As a consequence, applying Proposition 5.1, we obtain τ k Ψ j (u j ) = δ j,k u j for any j, k = 0 . . . J. Summing the latter identity over j yields τ k Ψ(u) = u k for all k = 0 . . . J.

Let us now assume that u = (u We have Ψ(u) = 0 ∀u ∈ X(Γ).

d , u n ) = (u j d , u j n ) J j=0 ∈ H(Γ) satises τ • Ψ(u) = u. By construction, the function ψ := Ψ(u) is well dened in R d \ Γ, it satises (γ -2 -∆)ψ = 0 in each Ω j , j = 0 . . .

Proof:

Denoting as before G x (y) := G (x -y), recall that we have Ψ Ω (u)(x) = [τ Ω (G x ), u] ∂Ω u ∈ H 1/2 (∂Ω) × H -1/2 (∂Ω) and all x ∈ R d \ ∂Ω. Plugging this expression into the denition of the multi-potential operator yields Ψ(u)(x) = τ (G x ), u ∀u ∈ H(Γ), ∀x ∈ R d \ Γ. Now observe that for any x ∈ R d \Γ we have τ (G x ) ∈ X(Γ) hence applying Proposition 4.1 concludes the proof.

A direct consequence of the lemma above is that Ψ d (u d ) = 0 for all u d ∈ X d (Γ), and Ψ n (u n ) = 0 for all u n ∈ X n (Γ). We deduce in particular that X n (Γ) ⊂ Ker(τ n • Ψ n ). The proof of the next result is similar to that of [START_REF] Claeys | Quasi-local multitrace boundary integral formulations[END_REF]Lem.6.4].

Lemma 5.3.

We have

p -τ n • Ψ n (p) ∈ X n (Γ) for any p ∈ H n (Γ).
p -τ n • Ψ n (p) = v n -τ n • Ψ n (v n ) + T(u d ) -τ n • Ψ(u). Applying Lemma 5.1 yields T(u d ) -τ n • Ψ(u) = 0. Besides we have Ψ n (v n ) = 0 according to Lemma 5.2 since v n ∈ X n (Γ). To summarise, we have just established p -τ n • Ψ n (p) = v n ∈ X n (Γ), which concludes the proof.
Combining the previous two lemmas, we see that

(τ n • Ψ n )(Id -τ n • Ψ n ) = 0.

We have

Ker(τ n • Ψ n ) = Range(Id -τ n • Ψ n ) = X n (Γ), and τ n • Ψ n : H n (Γ) → H n (Γ) is a continuous projector.
The next result is new and gives further details about the image of this projector. We have Range(τ n • Ψ n ) = T(X d (Γ)) so that τ n • Ψ n is an orthogonal projector with respect to the scalar product induced by T -1 over H n (Γ).

Proof:

Taking account of both Proposition 4.2 and 5.3, we see that it suces to prove τ n • Ψ n (T(u)) = T(u) for all u ∈ X d (Γ). Hence consider any u = (u j ) J j=0 ∈ X d (Γ). According to Lemma 5.2 we have τ n • Ψ d (u) = 0. As a consequence, applying Proposition 5.1, we obtain

τ k n • Ψ n (T(u)) = τ k n • ( Ψ d (u) + Ψ n (T(u)) ) = τ k n • J j=0 Ψ j d (u j ) + Ψ j n (T j (u j )) = J j=0 τ k n • Ψ j d (u j ) + τ k n • Ψ j n (T j (u j )) = T k (u k ) (24)
for any k = 0 . . . J. Since this holds for all k, we obtain that τ n • Ψ n (T(u)) = T(u), which concludes the proof.

From the previous results, we immediately obtain an estimate on the norm of the projection, which will be key in the analysis of Section 7. Although it shares similarities with [11, Prop.

the next result is new. Remark 5.1. The denition of Π stems from an orthogonal projection with respect to the scalar product induced by the operator T dened in 3.4. So the operator Π is tightly related to T. This is an important fact and we will come back to this later.

In the subsequent analysis, this projector will be the key tool for caracterising elements of X(Γ) and thus enforcing transmission conditions across interfaces. The next result indeed provides a caracterisation of the single trace space that is a core novelty of the present contribution.

Proposition 5.4.

Consider any ω > 0. With the notations of the previous corollary, for any

u = (u d , u n ) ∈ H(Γ), we have u ∈ X(Γ) if and only if u n -ıωT(u d ) = Π(u n + ıωT(u d )).
Proof:

According to Corollary 5.

1, for u = (u d , u n ) ∈ H(Γ), we have u n ∈ X n (Γ) ⇐⇒ (Id-Π)u n = 0 and u d ∈ X d (Γ) ⇐⇒ (Id + Π)T(u d ) = 0. On the other hand, Range(Id + Π) ∩ Range(Id - Π) = {0} since (Id + Π)/2 is a projector, which leads to u ∈ X(Γ) ⇐⇒ (Id -Π)u n = ıω(Id + Π)T(u d ).
Rearranging this latter identity yields the conclusion of the proof.

Reformulation of wave equations

In this section we focus on the wave equations ( 4) that we will reformulate in terms of traces only. We adopt the approach developed by Collino, Ghanemi and Joly in [START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF] and further studied and extended in [START_REF] Lecouvez | Quasi-local transmission conditions for non-overlapping domain decomposition methods for the Helmholtz equation[END_REF][START_REF] Lecouvez | Iterative methods for domain decomposition without overlap with exponential convergence for the Helmholtz equation[END_REF]. This approach generalises the original work of Després [START_REF] Després | Domain decomposition method and the Helmholtz problem[END_REF][START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle[END_REF][START_REF] Després | Domain decomposition method and the Helmholtz problem[END_REF][START_REF] Després | Décomposition de domaine et problème de Helmholtz[END_REF] on Optimised Schwarz Method for Helmholtz equation. In the present section, we will derive a convenient caracterisation of

C + (Γ) := C + (Ω 0 ) × • • • × C + (Ω J ) where C + (Ω j ) := { (τ j d (ϕ), µ j τ j n (ϕ)) ∈ H 1/2 (∂Ω j ) × H -1/2 (∂Ω j ), C + (Ω j ) := { div(µ∇ϕ) + κ 2 ϕ = 0 in Ω j and C + (Ω j ) := { ϕ κ 0 -outgoing if j = 0. }. (25)
The space C + (Ω j ) is closed in H 1/2 (∂Ω j )×H -1/2 (∂Ω j ) and we will use these spaces to reformulate the wave equation in each subdomain. The following proposition provides an important decomposition of the multi-trace space. Although it was already established in [START_REF] Claeys | Multi-trace boundary integral formulation for acoustic scattering by composite structures[END_REF]Prop.6.1] based on arguments of potential theory that assume piecewise homogeneity of the propagation medium, below is a new proof that also applies in the case of heterogeneous media.

Proposition 6.1.

We have the direct sum H(Γ) = X(Γ) ⊕ C + (Γ).

Proof:

Let us rst show that X(Γ)

∩ C + (Γ) = {0}. Pick some u ∈ X(Γ) ∩ C + (Γ) decomposed in Dirichlet/Neumann components u = (u d , u n ) with u d = (u j d ) J j=0 ∈ H d (Γ) and u n = (u j n ) J j=0 ∈ H n (Γ).
For each j = 0 . . . J, let φ j ∈ H 1 loc (Ω j ) refer to the unique functions satisfying div(µ∇φ j ) + κ 2 φ j = 0 in Ω j , φ 0 is κ 0 -outgoing,

(τ j d (φ j ), µ j τ j n (φ)) = (u j d , u j n ) on ∂Ω j . (26)
Set φ := 1 Ω 0 φ 0 + • • • + 1 Ω J φ J , so that div(µ∇φ) + κ 2 φ = 0 in each Ω j and, since u = (τ j d (φ), µ j τ j n (φ)) j=0...J ∈ X(Γ) the function φ satises transmission conditions across Γ, so that div(µ∇φ) + κ 2 φ = 0 in R d and φ is κ 0 -outgoing. Well-posedness of the Helmholtz equation with outgoing radiation condition leads to φ = 0, hence u = 0, which proves that

X(Γ) ∩ C + (Γ) = {0}. ( 27
)
Now let us consider the general case of an arbitrary u ∈ H(Γ). Consider any lifting function ψ ∈ L 2 (R d ) with compact support such that ψ | Ω j ∈ H 1 (Ω j ) and τ j d (ψ ) = u j d for all j = 0 . . . J.

Next dene ψ ∈ H 1 loc (R d ) as the unique element of H 1 loc (R d ) satisfying J j=0 Ω j µ∇(ψ + ψ ) • ∇ϕ -κ 2 (ψ + ψ )ϕ dx = u n , τ d (ϕ) ∀ϕ ∈ H 1 comp (R d )
and lim ρ→∞ ∂Bρ |∂ ρ ψ -ıκ 0 ψ| 2 dσ ρ = 0

where H 1 comp (R d ) refers to the elements of H 1 (R d ) that are boundedly supported. Existence and uniqueness of such a ψ stems from well posedness of Helmholtz problems in unbounded heterogeneous media, see e.g. [START_REF] Colton | Inverse acoustic and electromagnetic scattering theory[END_REF]Chap.3]. Applying a Green formula in each Ω j , we obtain div(µ∇(ψ + ψ )) + κ 2 (ψ + ψ ) = 0 in each Ω j , j = 0 . . . J lim ρ→∞ ∂Bρ |∂ ρ ψ -ıκ 0 ψ| 2 dσ ρ = 0.

Setting v = (τ j d (ψ + ψ ), µ j τ j n (ψ + ψ )) j=0,...,J , the equations above imply that v ∈ C + (Γ). De- composing in Dirichlet/Neumann contributions v = (v d , v n ), we have v d -u d = (τ j d (ψ)) J j=0 ∈ X d (Γ) since ψ ∈ H 1 loc (R d ).
Moreover, applying Green formulas once more in [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle[END_REF], we see v n , τ d (ϕ) = u n , τ d (ϕ) for all ϕ ∈ H 1 (R d ). Using the weak caracterisation of single trace spaces given by Proposition 4.1, we conclude that u d -v d ∈ X d (Γ) and u n -v n ∈ X n (Γ) hence, setting w := uv ∈ X(Γ), so that, with the decomposition u = v + w, we have established H(Γ) = X(Γ) + C + (Γ) which, together with [START_REF] Després | Domain decomposition method and the Helmholtz problem[END_REF], concludes the proof.

Let us underline that the previous proposition does not require unique solvability of local subproblems. It only assumes existence and uniqueness of the global scattering problem. This decomposition result can be regarded as analogous to Proposition 4.2 although, in the result above, the direct sum is a priori not orthogonal.

The next property relates to energy conservation considerations and will thus play a key role in the forthcoming convergence analysis. Although such a result is classical in scattering theory (see e.g [START_REF] Colton | Integral equation methods in scattering theory[END_REF]Thm.3.12]), the next lemma is rephrased so as to t our framework. We also provide a proof for the sake of completeness. Lemma 6.1.

We have ı[u, u] ∂Ω j ≤ 0, ∀u ∈ C + (Ω j ) ∀j = 0 . . . J, and thus ı u, u ≤ 0 ∀u ∈ C + (Γ).

Proof:

For any u ∈ H 1/2 (∂Ω j )×H -1/2 (∂Ω j ), let ϕ ∈ H 1 loc (Ω j ) satisfy div(µ∇ϕ)+κ 2 ϕ = 0 in Ω j and (τ j d (ϕ), µ j τ j n (ϕ)) = u on ∂Ω j . For all j = 0 . . . J, we have ı[u, u] ∂Ω j = 2 m{ ∂Ω j µ j τ j n (ϕ)τ j d (ϕ)dσ}. In the case where j = 0, the domain Ω j is bounded so that we can apply a simple Green formula on the later identity,

ı[u, u] ∂Ω j = 2 m{ ∂Ω j µ j τ j n (ϕ)τ j d (ϕ)dσ} = 2 m{ Ω j ϕdiv(µ∇ϕ) + µ|∇ϕ| 2 dx} = 2 m{ Ω j µ|∇ϕ| 2 -κ 2 |ϕ| 2 dx} = -2 Ω j m{κ 2 }|ϕ| 2 dx ≤ 0.
Note that, in the above inequality, we used that Assumptions (1) imply m{κ 2 } ≥ 0. In the case of Ω 0 take any radius ρ 0 > 0 large enough to guarantee R d \ Ω 0 ⊂ B ρ 0 . We can apply the same calculus as above, considering B ρ ∩ Ω 0 instead of Ω 0 . Taking account of the radiation condition satised by ϕ(x) for |x| → ∞, and the fact that m{κ 2 } is boundedly supported

(since κ(x) = κ 0 for |x| > ρ 0 ), we obtain ı[u, u] ∂Ω j = -2 Ω 0 m{κ 2 }|ϕ| 2 dx + 2 m{ ∂Bρ ϕ∂ ρ ϕ dσ} ≤ 2 m{ ∂Bρ ϕ∂ ρ ϕ dσ} = - 1 κ 0 ∂Bρ 2 e{ıκ 0 ϕ∂ ρ ϕ} dσ = 1 κ 0 ∂Bρ |∂ ρ ϕ -ıκ 0 ϕ| 2 dσ - 1 κ 0 ∂Bρ |∂ ρ ϕ| 2 dσ -κ 0 ∂Bρ |ϕ| 2 dσ ≤ 1 κ 0 ∂Bρ |∂ ρ ϕ -ıκ 0 ϕ| 2 dσ ∀ρ > ρ 0 ≤ lim inf ρ→∞ 1 κ 0 ∂Bρ |∂ ρ ϕ -ıκ 0 ϕ| 2 dσ = 0.
v n + ıαT(v d ) 2 Hn(Γ) = v n 2 Hn(Γ) + ω 2 v d 2 Hd(Γ) + 2α e{ı v d , v n } = v n 2 Hn(Γ) + ω 2 v d 2 Hd(Γ) -2α m{ v d , v n } = v n 2 Hn(Γ) + ω 2 v d 2 Hd(Γ) + ıα v, v for α = ±ω. (29)
We shall assume that the scalar coecient ω > 0, usually referred to as impedance, is xed until the end of this article. From the above identity we deduce an expression for the dierence between ingoing and outgoing traces. We have

v n + ıωT(v d ) 2 Hn(Γ) -v n -ıωT(v d ) 2 Hn(Γ) = 2ıω v, v for all v = (v d , v n ) ∈ H(Γ).
So-called ingoing/outgoing Robin trace operators also play an important role in scattering theory so, in the present paragraph, we study these trace operators in more detail. Dene

τ j ± : H 1 (∆, Ω j ) → H -1/2 (∂Ω j ) by τ j ± (φ) := µ j τ j n (φ) ± ıωT j (τ j d (φ)) for φ ∈ H 1 (∆, Ω j ), τ ± := diag j=0...J (τ j ± ). (30)
Remark 6.1. It is crucial for the subsequent analysis that the impedance operator involved in the denition of τ ± coincides with the DtN map T dened in 3.4 that induces scalar products on H d (Γ) and H n (Γ).

The Robin trace operators can be considered for prescribing boundary data for the solution of wave equations in each subdomain. Due to the positivity of the DtN maps T j , the associated boundary value problems are systematically well posed.

Well posedness of the problem above is a straightforward consequence of Lemma 6.2. We have now all the ingredients required for deriving a reformulation of our scattering problem (2). Assume that u ∈ H 1 loc (R d ) is the unique solution to (2), and denote f := Π(τ + (φ f )) with φ f dened by [START_REF] Bouajaji | Approximate local magnetic-to-electric surface operators for time-harmonic Maxwell's equations[END_REF]. Then the tuple of traces p = τ -(u) satises p ∈ H n (Γ) and p -(Π • S)p = f. [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] Proof:

If u ∈ H 1 loc (R d ) refers to the unique solution to (2) then (τ j d (u), µ j τ j n (u)) j=0...J ∈ X(Γ) so, according to Proposition 5.4, we have τ -(u) = Π(τ + (u)). In addition, the function u -φ f solves an homogenous Helmholtz equation in each subdomain i.e. (div(µ∇

• )+κ 2 )(u-φ f ) = 0 in Ω j for each j = . . . J and u-φ f is κ 0 -outgoing radiating, so (τ j d (u-φ f ), µ j τ j n (u-φ f )) j=0...J ∈ C + (Γ). As a consequence we have τ + (u -φ f ) = S • τ -(u -φ f ) = S • τ -(u). Thus we conclude that τ -(u) = ΠS(τ -(u)) + Πτ + (φ f ).
The structure of this new formulation is strikingly close to standard Optimised Schwarz Methods (OSM). This appears clearly when comparing [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] with 2 in [START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF], see in particular Formula ( 45) and ( 51) of this reference.

Here also [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] appears adapted to domain decomposition. In the operator Id -Π • S, the factor S is block-diagonal, each block being associated to a dierent subdomain, so that matrixvector product is trivially parallelisable. Of course, each block of S involves a DtN operator. Once Equation [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] has been solved numerically and that p = (p j ) j=0...J has been computed, the solution to the boundary value problem (2) can be recovered in each Ω j separately by solving the following local problems div(µ∇u) + κ 2 u = -f in Ω j , u is κ 0 -outgoing, τ j -(u) = p j ∀j = 0 . . . J. [START_REF] Gander | On the applicability of Lions' energy estimates in the analysis of discrete optimized Schwarz methods with cross points[END_REF] The main new feature of the formulation we present here is the transmission operator Π.

Contrary to the exchange operator traditionally used in OSM, see e.g. Formula (42) in [START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF], our transmission operator Π is not local anymore. But it only involves exponentially decaying kernels, with a damping factor γ that can be tuned, so that Π can nevertheless be considered quasi-local. In addition, various techniques (H-matrices [START_REF] Bebendorf | Hierarchical matrices[END_REF][START_REF] Börm | Ecient numerical methods for non-local operators[END_REF][START_REF] Hackbusch | Hierarchical matrices: algorithms and analysis[END_REF], Fast Multipole Method [START_REF] Darve | The fast multipole method: numerical implementation[END_REF][START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF]) can be used to sparsify this operator further. We have Π • S(p) Hn(Γ) ≤ p Hn(Γ) for all p ∈ H n (Γ).

Due to the continuity of Q, we obviously have u ω ≤ Q ω • g ω , where Q ω is dened with [START_REF] Gander | Cross-points in domain decomposition methods with a nite element discretization[END_REF]. On the other hand multiplying (37) on the left by P ± we obtain P + T(u d ) = ıω (

) 39
Since the projectors P ± are orthogonal for the scalar product (•, •) Hn(Γ) we obtain

Using the identity obtained in [START_REF] Greengard | A new version of the fast multipole method for the Laplace equation in three dimensions[END_REF] to replace P -(u n ) 2 Hn(Γ) + ω 2 P + T(u d ) 2

Hn(Γ) in the iden- tity above, using that P ± = (Id ± Π)/2, and observing that (T(v), p) Hn(Γ) = v, p , we obtain

Using Lemma 6.1, the real part of the previous identity is bounded from below by e{(p -ΠS(p), p) Hn(Γ) } ≥ 2 g 2 ω . We conclude by using (39).

Lax-Milgram lemma combined with the previous theorem yields bijectivity of Id -ΠS as an obvious outcome. The operator Id -ΠS : H n (Γ) → H n (Γ) is an isomorphism.

Solution strategy

Let us briey discuss how, in practice, to solve [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] i.e. an equation of the form p-Π•S(p) = f. First of all, since Π 2 = Id, this equation can be transformed into (Π -S)p = Π(f) = τ + (φ f )

which is practically more convenient as it avoids handling a product of operators. A general Krylov solver such as GMRes could be considered for solving this equation. We refer the reader to [51, chap.6] for more details on this solver.

Figure 1 :

 1 Figure 1: Examples of subdomain partitions in 2D (a & b) and 3D (c) with 4 subdomains (3 bounded + exterior). There is no cross point in (a), and cross points are red dots in (b) and red dashed lines in (c).

  , and analogous denitions for Hd(Γ) and H(Γ) . The multi-trace space H(Γ) coincides with H d (Γ) × H n (Γ) through a re-ordering of traces which is why, when considering an element u = (u j d , u j n ) J j=0 ∈ H(Γ), we will sometimes commit a slight abuse of notation writing "u = (u d , u n )" to refer to the Dirichlet components u d = (u j d ) J j=0 ∈ H d (Γ) on the one hand, and the Neumann components u n = (u j n ) J j=0 ∈ H n (Γ) on the other hand. There is a natural duality between Dirichlet and Neumann multi-trace spaces through the bilinear pairing

  J by denition of the potential operators Ψ j 's. It also satises τ j d (ψ) = u j d and τ j n (ψ) = u j n for all j = 0 . . . J since τ • Ψ(u) = u. By denition of the DtN maps given in 3.3 and 3.4, this can be re-written u j n = T j (u j d )∀j = 0 . . . J. Put in vector form, this yields u n = T(u d ).The proof of the following result closely follows [10, Lem.5.1] and[START_REF] Claeys | Quasi-local multitrace boundary integral formulations[END_REF] Lem.6.3].Lemma 5.2.

  p ∈ H n (Γ) and, applying Proposition 4.2, decompose it as p = v n +T(u d ) where u d ∈ X d (Γ) and v n ∈ X n (Γ). According to Lemma 5.2 we have Ψ d (u d ) = 0 so that, setting u := (u d , T(u d )), we have

  From this we deduce immediately the following proposition, which is a variant of [11, Cor.6.1]. Proposition 5.3.

Corollary 5. 1 .

 1 Dene Π := Id -2τ n • Ψ n . Then we have Π 2 = Id and the operators (Id ± Π)/2 are continuous projectors with X n (Γ) := Ker(Id -Π) and T(X d (Γ)) := Ker(Id + Π). Besides the following continuity estimate holds: Π(p) Hn(Γ) = p Hn(Γ) ∀p ∈ H n (Γ).

6 . 1

 61 Robin trace operators The caracterisation of X(Γ) provided by Proposition 5.4 involved specic combinations of Neumann and Dirichlet trace operators. Let us bring the attention of the reader to the following elementary identity: for any v = (v d , v n ) ∈ H(Γ), and any ω > 0 we have

7 . 2

 72 Well-posedness of the new formulationLet us examine the properties of the operator Id -Π • S in detail. First of all Π • S continuously maps H n (Γ) into H n (Γ). In addition, combining Corollary 5.1 and Proposition 6.2, we obtain a contractivity result.Lemma 7.2.

- 1 P 2 Hn 2

 122 + (f)/2 = T(g d ) P -(u n ) = P -(f)/2 = g n ⇒ g 2 ω = ω 2 P + T(u d ) 2 Hn(Γ) + P -(u n ) u ω ≤ Q ω (ω 2 P + T(u d ) 2 Hn(Γ) + P -(u n ) 2 Hn(Γ) ). Since u ∈ range(Q) = C + (Γ) we have u n + ıωT(u d ) = S(u n -ıωT(u d ))and, as a consequence, (37) implies (Id -ΠS)(u n -ıωT(u d )) = f hence, according to Proposition 7.1, p = u n -ıωT(u d ), which leads to the estimate p Hn(Γ) /2 ≤ u 2 ω ≤ Q 2 ω g 2 ω .

  ΠS(p), p)Hn(Γ) -ıωT(u d ) -Π(u n + ıωT(u d )), u n -ıωT(u d )) Hn(Γ) = (P -(u n ) -ıωP + T(u d ), u n -ıωT(u d )) Hn(Γ) = P -(u n ) 2 Hn(Γ) + ω 2 P + T(u d ) 2Hn(Γ)-ıω(P + T(u d ), u n ) Hn(Γ) + ıω(P -(u n ), T(u d )) Hn(Γ)

  ΠS(p), p)Hn(Γ) = g 2 ω -(ıω/2) u d , u n + (ıω/2) u n , u d -(ıω/2)(ΠT(u d ), u n ) Hn(Γ) -(ıω/2)(Π(u n ), T(u d )) Hn(Γ) = g 2 ω -(ıω/2) u, u -ıω e{(ΠT(u d ), u n ) Hn(Γ) }

  ∂Ω), see e.g. [53, 6.6.3], [45, Chap.4] or [48, 1.2.2.4]. To be more explicit T Ω is dened by

  belongs to X d (Γ) if and only if u, p = 0∀p ∈ X n (Γ) and that, similarly, p ∈ H n (Γ) belongs to X n (Γ) if and only if u, p = 0 ∀u ∈ X d (Γ). The next result provides a new variant of the decomposition established in [14, Prop.6.1].

	Proposition 4.2.

Acknowledgement

This work received support from the French National Research Agency (ANR) through the NonlocalDD project, grant ref. ANR-15-CE23-0017-01.

Lemma 6.2.

For any g ∈ L 2 (Ω j ) with bounded support, and any h ∈ H -1/2 (∂Ω j ), there exists a unique φ ∈ H 1 loc (Ω j ) such that div(µ∇φ) + κ 2 φ = g in Ω j , and τ j -(φ) = h on ∂Ω j (and φ is κ 0outgoing if j = 0).

The proof of the previous lemma is a basic exercise on variationnal formulations, so it is left to the reader. We need to introduce resolvent operators that solve Helmholtz equation in each subdomain with a prescribed outgoing Robin boundary trace, the operator S j : H -1/2 (∂Ω j ) → H -1/2 (∂Ω j ) dened by

The following result, very similar to [START_REF] Després | Méthodes de décomposition de domaine pour les problèmes de propagation d'ondes en régime harmonique. Le théorème de Borg pour l'équation de Hill vectorielle[END_REF]Lem.4.1], was established in [START_REF] Lecouvez | Iterative methods for domain decomposition without overlap with exponential convergence for the Helmholtz equation[END_REF]Lem.3,Chap.2]. We provide the proof for the sake of completeness. Proposition 6.2.

The operator S = diag j=0...J (S j ) continuously maps H n (Γ) into H n (Γ) and is non-expansive: for all p ∈ H n (Γ) we have S(p) Hn(Γ) ≤ p Hn(Γ) .

Proof:

Pick an arbitrary p = (p j ) J j=0 ∈ H n (Γ). Applying Lemma 6.2, there exist functions φ j ∈ H 1 loc (Ω j ) such that div(µ∇φ j ) + κ 2 φ j = 0 in Ω j , and τ j -(φ j ) = p j on ∂Ω j (and φ j is κ 0 -outgoing if j = 0). Set v = (v d , v n ) := (τ j d (φ j ), µ j τ j n (φ j )) j=0,...,J , we have v n -ıωT(v d ) = p and v n + ıωT(v d ) = S(p). Since v ∈ C + (Γ) by construction, combining Corollary 6.1 and Lemma 6.1 concludes the proof.

The previous result shows that the scattering operator S is a contraction but it is not a priori an isometry. In the context of Problem (2), this is due to energy loss through radiation of waves toward innity and absorption properties of the propagation medium (positive imaginary part of κ 2 ).

Reformulation of the scattering problem

In the present section we describe a reformulation of the scattering problem (2) as an equivalently well posed problem.

Derivation of the formulation

To take account of the right hand side f ∈ L 2 (R d ), we introduce the oset function φ f ∈ L 2 loc (R d ) whose restriction to each subdomain φ f | Ω j belongs to H 1 loc (Ω j ) and is the unique

Coming back to Remarks 5.1 and 6.1, we stress that the estimate above is made possible by the fact the very same operator T comes into play in the denition of both Π and S (through the impedance trace operators τ ± ).

A direct consequence of this property is that the numerical range of the operator Id -Π • S is located in the complex right-half plane C + := {z ∈ C, e{z} ≥ 0}. This is denitely an interesting feature from the perspective of linear solvers. Next this operator is also one-to-one. Consider a p = (p j ) J j=0 ∈ H n (Γ) satisfying p = ΠS(p). Consider the function v ∈ L 2 loc (R d ) such that, its restriction in each subdomain v| Ω j belongs to H 1 loc (Ω j ) and satises div(µ∇v) +

. Hence, applying Proposition 5.4, we deduce that v ∈ X(Γ). Since, on the other hand, we have v ∈ C + (Γ) by construction, we conclude that v ∈ C + (Γ) ∩ X(Γ) = {0} according to Proposition 6.1. Hence p = v n -ıωT(v d ) = 0.

The next theorem is the main novelty of the present contribution. It shows that [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] oers a new strongly coercive formulation of the scattering problem [START_REF] Allaire | Numerical linear algebra[END_REF]. As discussed in [START_REF] Moiola | Is the Helmholtz equation really sign-indenite?[END_REF], this directly contradicts the widespread belief that harmonic wave propagation systematically leads to sign indenite formulations. There exists α > 0 such that e{((Id -Π • S)p, p) Hn(Γ) } ≥ α p 2 Hn(Γ) for all p ∈ H n (Γ).

Proof:

We need rst to introduce a few notations that we shall use only for this proof. According to Proposition 6.1, there exists a bounded projection operator Q : H(Γ) → H(Γ) with range(Q) = C + (Γ) and ker(Q) = X(Γ). For convenience, we set

Because ω > 0 is a simple xed positive constant, ω and H(Γ) are equivalent norms, and continuity of the projection Q is exactly equivalent to the boundedness of Q ω . We shall also consider the bounded orthogonal projectors P ± : H n (Γ) → H n (Γ) dened by

Now pick an arbitrary p ∈ H n (Γ).

where β ∈ (0, 1) is a relaxation parameter. Following Theorem 7 and Remark 9 in [START_REF] Collino | Domain decomposition method for harmonic wave propagation: a general presentation[END_REF],

a rough estimate can be derived for the convergence of Richardson's linear solver in this case. Let p ∞ refer to the unique solution to [START_REF] Bouajaji | A quasi-optimal domain decomposition algorithm for the time-harmonic Maxwell's equations[END_REF] and set e (n) := p ∞ -p (n) so that e (n+1) = ((1 -β)Id + βΠS)e (n) . Recall the convexity identity

which holds for any x, y ∈ H n (Γ) and any β ∈ (0, 1).

Hn(Γ)

In this estimate, the convergence factor (1 -α 2 β(1 -β)) 1/2 < 1 is thus minimized for β = 1/2 and takes the value (1 -(α/2) 2 ) 1/2 in this case.


published or not. The documents may come   L'archive ouverte pluridisciplinaire


Khouloud Dorgham
email: khouloud.dorgham@univ-artois.fr

Issam Nouaouri

Jean-Christophe Nicolas

Gilles Goncalves

-Morocco "New advances and challenges for sustainable and smart industries" A COLLABORATIVE SUPPLY CHAIN NETWORK DESIGN WITHIN A TERRITORY HOSPITAL GROUP

Keywords: Optimization, horizontal collaboration, hospital supply chain, logistics pooling

The management of logistics functions is essential for the overall good functioning of a health-care establishment to meet the efficiency and the effectiveness requirements that hospitals are increasingly faced with. Minimizing spending is currently leading health-care establishments to reason and optimize the physical flows in terms of the overall performance of their supply chains. Therefore, logistics pooling could be seen as a solution for hospitals to reduce costs and enhance the quality of service. In this work, linear programming (LP) model was proposed to design a shared hospital supply chain among various establishments of a "Territory Hospital Group" based on a horizontal collaborative logistic strategy. The objective is to rationalize, pool, and optimize the storage and the distribution of products between suppliers, warehouses, and cross-docks. Instances with small and medium sizes were generated based on a real situation and several tests were developed on different hypotheses to reveal the impact of the proposed logistics collaboration on economic costs.

INTRODUCTION

In recent decades, the grouping of establishments or communities (university communities, etc.) is becoming of great importance, and several research studies aimed at promoting new models of territorial organization in different areas of real life. Since July 2016, hospitals in France were also affected by this grouping, cause medical sector has been faced new challenges and issues that obliges to reorganize the activity of its different functional units. Encountered with these deep changes, it is essential to revise the structure of the health-care system in order to optimize the resources mobilized and control spending [START_REF] Pillay | Managerial competencies of hospital managers in south africa: a survey of managers in the public and private sectors[END_REF].

Therefore, in the interest of enhancing hospital logistics (managing products and material flows and distribution circuits), France hospitals have met to form Territorial Hospital Groups (THG) that vary mainly according to their establishment's parties, their budget, and the territories served. This reform aimed at pooling progressively certain support functions provided by the establishment such as logistic management and increasing cooperation between hospitals.

Following other sectors, collaboration strategy will be advantageous through an overall cost reduction, a supplier integration, and an optimization of processes such as storage of products in warehouses and pharmacies and their distribution to care units. Hence, to model this collaboration, mastery of logistics functions is essential for the overall good functioning of the hospital supply chain.

Generally, in a classic supply chain, every single hospital must review, release, treat, and monitor its supply cycle individually. The supplier's organization responds to this operation and must meet the unit care specifications. In this work, a linear programming model was proposed, at the strategic level of decision making, to design a pooled hospital supply chain among various establishments of a territory hospital group based on a horizontal collaborative logistic strategy.

The remainder of this paper is organized as follows: Section 2 reviews the literature on collaborative logistics network design problems with different real-word applications domains and methodologies. Section 3 presents a description of the proposed problem and its mathematical formulation.

In the fourth section, several instances were generated to demonstrate the pooling performance. Finally, Section 5 presents a conclusion with some future research studies.

Literature review

According to [START_REF] Moutaoukil | Modeling a logistics pooling strategy for agri-food smes[END_REF], collaboration strategy consists to share logistics means and resources in order to minimize costs and increase profits, it could be either on a vertical or horizontal level. The first type concerns partners who belong to the same logistics chain that operates at different levels of the supply network. Unlike the second type, that concerns partners of the same level who do not belong to the same logistics chain. Several collaborative supply chain design problems and models have been developed in the literature dealt with the design of a shared logistics network [START_REF] Mrabti | Vers un modèle de simulation de la mutualisation logistique 4.0[END_REF] to increase the efficiency of the supply chain and achieve economies of scale in a different real-world domain (agronomy, automobile, etc.) and for all decision-making levels.

In this section, we identified different research studies belonging to the strategic level of decision making that can be classified according to their field of application and their resolution approach.

Simulation is the most cited techniques in the literature to solve this problem.

As it was presented in [START_REF] Giso | GPA et GPA Mutualisée dans la grande distribution[END_REF], authors used simulation to compare traditional supply chain against vertical and horizontal collaborative strategies and confirmed that horizontal collaboration is the best strategy for reducing logistics costs. In [START_REF] Pooley | Modeling and evaluating shipment consolidation in a logistics system[END_REF], authors proposed a simulation approach to evaluate a logistic transportation consolidation strategy within a food manufacturing business to reduce costs by combining several types of logistical pooling; (1) pooling of warehouses, (2) pooling of platforms, and (3) pooling of means of transport. [START_REF] Wanke | Consolidation effects: Whether and how inventories should be pooled[END_REF] developed a simulation tool to determine the impact of inventory centralization and regular transshipment inventory-pooling models, on logistic costs and indicated the best pooling scenario. [START_REF] Leitner | Structural concepts for horizontal cooperation to increase efficiency in logistics[END_REF] used simulation techniques to evaluate pooling strategy on two projects in the automobile sector in Romania and Spain and optimize cost structures. Also, [START_REF] Moutaoukil | Modeling a logistics pooling strategy for agri-food smes[END_REF] uses several scenarios to compare the performance of a traditional logistics network against a pooling supply chain with a horizontal collaborative strategy and they used simulation technique to manage the agri-food SMEs flows. In 2017, [START_REF] Makaci | Pooled warehouse management: An empirical study[END_REF] developed a generic simulation model to solve a multiple case study approach that dealt with the management of shared warehouses, a list of indicators was presented to assess the performance of logistical pooling. Recently, in health-care domain, a simulation tool was developed by [START_REF] Nicolas | Méthode et outil pour la rationalisation des magasins et des flux au sein des groupements hospitaliers de territoire. 9ème Conférence Francophone en gestion et ingénierie des systèmes hospitaliers[END_REF], the decision-maker can create, choose and compare different pooling scenarios in order to measure their impact on the local maintenance and/or pooling of product flows among hospitals.

Other studies proposed heuristic optimization approaches and exact methods to solve the collaborative supply chain design problem. Authors in [START_REF] Groothedde | Towards collaborative, intermodal hub networks: A case study in the fast moving consumer goods market[END_REF] used heuristics to search the best combination of hubs in a pooled logistics chain in order to minimize costs. [START_REF] Nataraj | Consolidation centers in city logistics: A cooperative approach based on the location routing problem[END_REF] proposed a meta-heuristic approach to solve different horizontal collaboration scenarios in the transport domain in order to increase the vehicle filling rate and decrease economic cost. In [START_REF] Cheong | Logistics network design with supplier consolidation hubs and multiple shipment options[END_REF], a LP model was proposed to solve collaborative network design problem and decide the number, location, and operation of consolidation hubs by minimizing upstream and downstream costs of warehouses. [START_REF] Tuzkaya | A holonic approach based integration methodology for transportation and warehousing functions of the supply network[END_REF] proposed a LP model to solve a two-echelon supply chain problem for distributing automotive industry products from the suppliers to the warehouse and from the warehouse to the manufacturers to maximize profit.

For the best of our knowledge, in the operational research field, none of the existing research works has modeled or studied the impact of horizontal collaboration and pooling strategy in hospital sector and notably within territorial hospital groups, which provides us a strong motivation to study. This paper could be an extension of the work presented by [START_REF] Nicolas | Méthode et outil pour la rationalisation des magasins et des flux au sein des groupements hospitaliers de territoire. 9ème Conférence Francophone en gestion et ingénierie des systèmes hospitaliers[END_REF] where authors allow decisionmakers to choose and compare the pooling scenarios of products below a THG based on a set of criteria chosen by hospital partners. However, an optimal product pooling scenario could not be generated. In this paper, we propose an optimization linear programming model that offers an optimal product pooling scenario for the decision-maker within a THG to optimize the economic performance of hospital sector among the management of several logistics costs. Therefore, a multi-supplier, multiwarehouse, and multi-product network is considered as a collaborative pooled warehousing structure. The problem is designed as a minimum-cost flow graph generalized to several products.

PROBLEM DESCRIPTION

In our study, the hospital logistics chain could be designed as a layered network with |S| suppliers, |W | warehouses, and |P | products subfamily (food, cleaning materials, textiles, medicines, etc.).

Knowing that every warehouse has its managing strategy that characterizes its purchasing, procurement, and storage activities and it is dedicated for serving only one establishment. The process flow consists of three steps: (1) making a procurement order from suppliers, (2) shipment of products to the warehouse of the considered hospital, and (3) the distribution of these products to the unit care. Knowing that each care unit can only be supplied by its warehouse and that each warehouse manages its procurement, warehousing, and transportation activity independently. Such a supply chain problem can generate several logistics costs such as transportation costs, inventory holding costs, and ordering costs, etc. Which explains the request to migrate for a horizontal collaborative pooling strategy and to form a territorials hospital groups.

Therefore, our work consists in designing a pooled supply chain within various establishments of the THG. The objective is to find an optimal allocation of product flows, that are distributed from suppliers to warehouses and from warehouses to establishment, and to set up a pooling scenario that groups one or more sub-families (materials, food, etc.) of products in the suitable stores. Thus, in the present hospital logistics network, there are storage warehouses that hold the stock of one or more sub-families of products, and cross-dock stores that represent a point of material handling and distribution, where products are not stored for an extended time period. Knowing that a store can be, simultaneously, as a storage warehouse for one or more subfamilies and as a crossdock for other subfamilies.

Our pooling scenario is carried out in two stages, (1) the placement of certain sub-families of products on one or more warehouses and (2) their distribution from these warehouses to one or more cross-docks. The proposed model allows to specify optimally for each product sub-family, its source (which supplier), its storage locations (warehouse), and its distribution to cross-docks stores. The objective is to minimize the overall economic costs defined as follows:

• Full-Time Equivalent cost (FTE): represents the workload of employees.

• Transportation cost: denotes expenses related to the distribution of products from storage warehouses to cross-docks.

• Purchasing cost: concerning product prices set by suppliers.

• Ordering cost: generated during the management of orders which vary according to their annual number (administrative and logistical monitoring, reception and handling charges, etc.).

• Holding cost: related to the inventory storage (insurance, depreciation of facilities, rental, and maintenance of premises, etc.).

Different constraints should be respected; each facility's demand for a product must be satisfied and the maximum storage capacity of warehouses should not be exceeded. In this setting, we make the following assumptions:

• Suppliers have unlimited delivery capacity.

• The supply strategy (i.e. procurement periods, unit costs, etc.) of the warehouse where the pooling of products takes place is maintained.

• The number and locations of warehouses and cross-docks are assumed to be fixed and known.

• A given product could be distributed for warehouses by one or more suppliers at different prices.

• The product price proposed by a given supplier is fixed for all warehouses.

The notation sets, parameters, and decision variables used in the model are presented below. (2)

Decision variables

Objective function

In the economic objective function, we introduced the different costs that affect the pooling performance. It is represented by the following formula:

M inimize P C + HC + T C + F C + OC (3)
P C = P S W P C p,w,s x p,s,w (4)
HC = P W HC p,w IL p,w (5)
T C = P W C T C p,w,c y p,w,c (6)
F C = P S W
(F C1 p,w x p,s,w + F C2 p,c y p,w,c ) ( 7)

OC = P S W
x p,s,w OC p,w ( t P P p,w )

The economic function aims to minimize the summation of five types of costs; the purchasing cost, inventory holding costs at the warehouse, transportation cost, FTE cost and finally total ordering cost.

Constraints

The constraints of our model are as follows: Constraints 9 ensures that the demand of the unit care for each subfamily of product is satisfied. Equation 10 ensure that the total products quantity at each warehouse should not exceed its storage capacity.

Constraint 11 represents the balance among supplying, inventory and deliveries at each warehouse. Constraints 13 and 12 represent the types of decision variables and guarantee the positivity of the flows.

EXPERIMENTAL RESULTS

In this section, we present computational tests and analysis results to validate the computational efficiency and effectiveness of the model and to see the impact of the horizontal collaboration on our economic objective function.

Two different configurations are used, firstly we considered the pre-pooling scenario where products are delivered directly from supplier to unit care, then we considered the polling scenario, where collaboration between functional units of the THG is authorized (i.e. products can be stored in the warehouse of another unit care before being shipped). The experiments have been complimented on Intel(R) Core(TM) Duo i3-5005U 2.00 GHz and 4 Gb RAM. They have been solved using IBM ILOG CPLEX solver.

Data generation

For lack of real data, we performed computational experiments on a set of randomly generated test instances based on realistic parameter value obtained partially from a real situation of an existing territory hospital group. We considered a set of 8 instances that vary according their size and their cost structure. According to the number of warehouses, we can form two groups of instances (small and medium). Table 1 summarizes the parameters of each instances. The size of an instance is given by the number of potential warehouses (|W |) and the number of products (|P |). Continuous uniform distributions, MOSIM'20 -November 12-14, 2020 -Agadir -Morocco independent from each other, were considered in the random number generation of all variables. The number of suppliers is fixed and it is equal to two, the warehouse storage capacity and the demand of the care unit are generated randomly. The unit costs structure is determined as follows; the unit ordering cost has been supposed fixed (20e) for all warehouses and all types of products. Next, for each product p ∈ P , a purchasing unit cost is chosen randomly in the interval [1,20]. Also, for each warehouse m ∈ M , a fixed unit holding cost is set between [20,30] for all products of this warehouse.

Results and discussion

A comparison between the pre-pooling and the pooling scenarios will be considered to evaluate the effect of collaboration (is it advantageous for hospitals?). The pre-pooling scenario (SC1) presents the direct shipment of products from the supplier to the distribution center, and the pooling scenario (SC2) is obtained by solving the model presented previously with the objective economic function. For all instances, the optimal logistic cost of the prepooling scenario is usually lower than the pooling scenario, which can confirm that the pooling strategy within the THG is advantageous and provides a remarkable improvement at the economic level. Indeed, there is a significant reduction of the total cost that varies, depending on the instance, between 13% and 17% for small ones and between 11% and 29% for medium instances. This is reflected by an improvement in FTE cost, purchasing cost, ordering cost, and inventory holding cost. The CPU time (s) required by CPLEX for each scenario SC1 and SC2 is given in the third and fifth column, respectively. It can be seen that the CPU time increases with the number of warehouses and products, and the pooling scenario is always time-consuming than the pre-pooling configuration, e.g. the CPU time of SC1 is 80% lower than SC2 in M4.

Figure 1 and Figure 2 demonstrate economically the details of costs saving realized, they represent the average of each cost among all instances. Despite the appearance of an additional transportation cost, we note that we make gains in terms of other costs; 13% for the FTE cost since the solver will choose warehouses with lowest labor costs, and 29% for purchasing cost. In addition, a reduction in inventory holding cost and ordering cost with respectively 19% and 14%. This gain results generally from the optimal pooling of multiple products between warehouses according to better flow allocation. Horizontal collaboration does not only have an impact on logistics costs, but also it could generate additional added value on other factors. Therefore, the resulting analysis may be projected on different important dimensions such as the warehouse filling rate. This makes it possible to quantitatively evaluate the initial situation comparing by the collaborative configuration. Furthermore, it is possible to deduct the amount of products exchanged and grouped after collaboration for each instances by presenting the percentage of pooling products (Table 3).

For each instance, a total filling rate is measured among the overall warehouses before and after collaboration. We note that for all instances, the filling rate in SC2 decreases with an improvement ranging from 1% to 2%. Therefore, we can confirm that horizontal collaboration ensures better stock management and allows to save more free space and hold even more products. The fourth column represents the number of storage warehouses (#SW) kept open after collaboration. Moreover, according to pooled product percentage, it is noticed that for the majority of instances, more than 50% of products have been pooled, which confirms that collaboration is usually more advantageous.

Sensitivity analysis

In this section, we determine how target variables are affected based on changes in input parameters. It is a way to predict how changes in coefficients of the model can affect the optimal solution. Considering the importance of transportation cost generated during pooling, we motivated the analysis by increasing the unit transportation cost and assuming that warehouses are far apart to show the sensitivity of the model to this parameter. This study was carried out only on instance S1. Figure 3 shows the changes in the overall costs following the unitary increase in transportation cost. As it is clearly shown in Figure 3, when the unit transportation cost increases up to 45%, the total cost (objective value) gradually increases, while the other logistics costs remain constant. Above 45%, not only the transportation cost is increased but also all other logistics costs are influenced. Therefore, we can note that if the unit transportation cost is raised to a maximum of 45%, the optimal procurement and distribution plan (optimal solution) remains unchanged, but among 46% the optimal solution is no longer maintained. On the other side, we evaluate if the decrease of transportation cost affects the optimal solution and it was revealed that the coefficient can be decreased without bound by influencing only the cost of the objective value and not the optimal solution.

Conclusions

In this paper, we consider a collaborative supply chain within a territory hospital group in a multiproduct, multi-supplier, and multi-warehouse setting.

Our main objective was to demonstrate the pooling performance and improve the economic level of the supply chain. A LP optimization model was designed to organize product pooling within hospitals in order to minimize the overall logistic costs, including ordering cost, transportation cost, inventory holding cost, FTE cost, and purchasing cost. The proposed linear programming model can help managers to decide not only which product sub-family to pool, but also the storage locations of those products and the quantity transported between facilities to satisfy demand. The network presented in this study deals with capacitated facilities (warehouses/cross-dock). We test our model on small and medium instances generated. Good results are obtained; as logistics cost and warehouse filling rate.

This study can be extended in several directions. It will be possible to incorporate other different aspects of sustainable development such as environmental objectives in the context of horizontal collaboration. Also, to specify the problem for a more realistic scenario, we should deal with uncertain demand over the decision horizon.

  s,w -C y p,w,c = 0, ∀w ∈ W, p ∈ P (11) x p,s,w ≥ 0, ∀w ∈ W, p ∈ P, s ∈ S (12) y p,w,c ≥ 0, ∀w ∈ W, p ∈ P, c ∈ C (13)

Figure 1 :

 1 Figure 1: Transportation, holding and ordering costs variation

Figure 2 :

 2 Figure 2: FTE and purchasing costs variation

Figure 3 :

 3 Figure 3: Sensitivity analysis when transportation cost is increasing

  x p,s,w : quantity of product p transported from supplier s to warehouse w. y p,w,c : quantity of product p transported from warehouse w to cross-dock c.

	I p,w : Average inventory level of product p in a
	warehouse w:
	I p,w =	x p,s,w 2 * ( t P Pp,w )	(1)
	IL p,w :		Inventory value of products p in
	warehouse w:
	IL p,w =	P C s,w,p x p,s,w 2 * ( t P Pp,w )
				a
				warehouse w;
				F C1 p,w : Full-time equivalent unit cost of
				product p in warehouse w;
				F C2 p,c : Full-time equivalent unit cost of product
				p in cross-dock store c;
				T C p,w,c : Unit transportation cost of product p
				from warehouse w to cross-docks c;
				C w : maximum storage capacity of warehouse w;
				d p,c : demand of product p by cross-dock c;
				a p,w : Unit surface occupied by product p in the
				warehouse w (m 2 );

Sets S: set of suppliers, |S| =1..s; W : set of warehouses, |W | =1..w; C: set of cross-docks, |C| =1..c; P : set of products, |P | = 1..p; Variables / expression P C: Total purchasing cost; OC: Total ordering cost; T C: Total transportation cost; HC: Total holding cost; F C: Total Full-time equivalent cost; P C p,w,s : Unit purchasing cost of products p by the warehouse w from the supplier s; HC p,w : Possession rate of product p in a warehouse w;

OC p,w :Unit ordering cost of product p for t: calendar days; P P p,w : procurement period of product p for warehouse w (P P p,w = 0);

Table 2

 2

	Inst	SC1	CPU	SC2	CPU Gain
	S1	2160 10 3	2.5	1780 10 3	4.2	17.6%
	S2	3200 10 3	2.8	2647 10 3	7.0	17.2%
	S3	3876 10 3	112.4 3356 10 3	135.2 15.5%
	S4	4652 10 3	94.0	4111 10 3	162.4 13.1%
	M1	1316 10 4	0.3	1166 10 4	0.9	11.3%
	M2	1744 104	0.7	1238 10 4	1.4	29.0%
	M3	2226 10 4	124.0 1732 10 4	342.0 22.2%
	M4	2447 10 4	114.6 2174 10 4	562.0 11.2%

summarized the optimal solutions of different instances for the two scenarios.

Tableau 2: Optimal total cost for pre-pooling and polling scenarios

Acknowledgements

This work is funded byçà "Agence Regionale de la Santé" of Hauts-de-France Region.