Louis-Brahim Beaufort

Pierre-Yves Massé

Antonin Reboulet

Laurent Oudre

Network Reconstruction Problem for an Epidemic Reaction-Diffusion

We study the network reconstruction problem for an epidemic reactiondiffusion. These models are an extension of deterministic, compartmental models to a graph setting, where the reactions within the nodes are coupled by a diffusion. We study the influence of the diffusion rate, and the network topology, on the reconstruction and prediction problems, both from a theoretical and experimental standpoint. Results first show that for almost every network, the reconstruction problem is identifiable. Then, we show that the faster the diffusion, the harder the reconstruction, but that increasing the sampling rate may help in this respect. Second, we demonstrate that it is possible to classify symmetrical networks generating the same trajectories, and that the prediction problem can still be solved satisfyingly, even when the network topology makes exact reconstruction difficult.

Introduction

Network reconstruction problems, in which one aims at reconstructing a network structure from the observation of a signal evolving on it, is an important topic of current research, spanning over numerous domains (Marc Timme and Casadiego 2014; Shandilya and Marc Timme 2011; [START_REF] Dong | Laplacian matrix learning for smooth graph signal representation[END_REF]; Le Bars et al. 2019; Sardellitti, Barbarossa, and Lorenzo 2019). Indeed, the widespread use of networks as a modelling tool in fields as diverse as telecommunications (Pastor-Satorras and Vespignani 2004; Newman, Watts, and Strogatz 2002), genetics [START_REF] Gardner | Inferring Genetic Networks and Identifying Compound Mode of Action via Expression Profiling[END_REF]; [START_REF] Karlebach | Modelling and analysis of gene regulatory networks[END_REF], ecology [START_REF] Hanski | Metapopulation Biology: Ecology, Generics, and Evolution[END_REF]; Tamburello, O. Ma, and M. Côté 2019), or transportation of goods or humans [START_REF] Youn | Price of Anarchy in Transportation Networks: Efficiency and Optimality Control[END_REF][START_REF] Perfido | Towards Sustainable Water Networks: Automated Fault Detection and Diagnosis[END_REF], to name but a few, makes understanding the connections between their structure, or internal properties, and the phenomena which happen over them, a crucial issue.

Recently, Prasse and Piet Van Mieghem (2020) have addressed the reconstruction problem for a wide class of epidemiological models (Sahneh, C. Scoglio, and P. Van Mieghem 2013). These models have gained considerable attention since the early 20 th century, following notably the classic works of [START_REF] Ogilvy | A contribution to the mathematical theory of epidemics[END_REF]. In those, individuals are categorized in compartments which describe their status with respect to an infectious disease, and the models describe the way they transition from compartments to compartments as the disease spreads through contacts, and they react (heal) to it (Odo Diekmann, H. Heesterbeek, and Britton 2012). Quickly, the early scalar models have been enhanced, by embedding them into networks (Pastor-Satorras and Vespignani 2001; Pastor-Satorra et al. 2015; Nowzari, Preciado, and Pappas 2016), in order to refine the analysis of the influence of contacts between individuals, on the spread of the disease.

In their work, Prasse and Piet Van Mieghem (2020) ask two questions: first, can the network structure be retrieved from the observation of the dynamics? Secondly, even in the case of a negative answer, is it possible to approximate the structure well enough to predict the future evolution of the disease? Even if the problem they study has a linear structure, they show the answers to these questions are not straightforward. We propose to address the same questions on another very important class of network-epidemiological models, the epidemic reaction-diffusion models, also known as metapopulation models with explicit movement [START_REF] Arino | Diseases in Metapopulations[END_REF]. Just like the model studied in [START_REF] Prasse | Network Reconstruction and Prediction of Epidemic Outbreaks for General Group-Based Compartmental Epidemic Models[END_REF], the nodes of the graph represent sub-populations: for instance, the cities in the transportation network of a country. However, the interactions between populations is no longer described by a static contact structure, but by a diffusion. Accordingly, the internal dynamics of each sub-population follow a standard deterministic epidemiological model (SIS, SEIR, ...) (Odo Diekmann, H. Heesterbeek, and Britton 2012) while flows of individuals go from node to node through a diffusion. Following their apparition in population dynamics in the 1970's, these models have since gained considerable attention in the field of mathematical epidemiology (Brauer and [START_REF] Arino | Spatio-temporal spread of infectious pathogens of humans[END_REF]). However, we are not aware the inverse problem has been studied for these models, up to now.

The standard network reconstruction procedure (Marc Timme and Casadiego 2014) is the optimisation of some regression error. However, lack of identifiability, and bad conditioning, may prevent it from being highly efficient. Structural and qualitative analysis of the model is therefore of much importance to better understand the dynamics at hand, and guide the reconstruction work. Moreover, such analysis may provide insights for other models incorporating diffusion as well [START_REF] Haehne | Detecting Hidden Units and Network Size from Perceptible Dynamics[END_REF]. Our contributions are therefore the following. On the one hand, we conduct theoretical analysis on the influence of the diffusion rate on the reconstruction, and illustrate our results by experiments. On the other hand, we study the influence of network topology, both theoretically, using notably the notion of graph automorphisms (Simon, Taylor, and We first present background material in Section 2. Next, we present the problem we address, conduct some initial identifiability analysis, and describe the experimental setup, in Section 3. Then, in Section 4, we study the influence of the diffusion rate. Finally, in Section 5, we study the influence of the network topology. The proofs are deferred to the appendices. The code for the experiments, implemented in Python, is available on the git repository: https://reine.cmla.ens-cachan.fr/ masse/network_reconstruction_reaction_diffusion.

Background

We first present the classical epidemiological models (Section 2.1), before giving a short overview on network reconstruction techniques (Section 2.2). Then, we present our contributions (Section 2.3). Finally, we introduce our notations (Section 2.4).

Deterministic, Compartmental Epidemiological Models

Deterministic, compartmental epidemiological models represent the propagation of a disease within a population by first segmenting the population in compartments, describing the status with respect to the disease (Odo Diekmann, H. Heesterbeek, and Britton 2012). Classical compartments include the "susceptible" (S), which gathers people which may contract the disease when confronted to "infected" (I) people, who later will have "recovered" (R). Transitions from compartments to com-partments are governed by differential equations. One simple and generic model, which we use for simplicity throughout our study, is the SIR model. Three scalar functions s, i and r track the numbers of people in each compartment, and they evolve according to, for all t ≥ 0,

             ds dt = -βsi di dt = βsi -δi dr dt = δi. (1)
Here, β and δ are positive real numbers. The parameter β is often called the infection rate, and δ is the curing rate. The quantity δ -1 may be interpreted as the average time an individual remains infected, before healing (O. Diekmann, J. A. P. Heesterbeek, and Roberts 2009). The fact it is positive means people heal in finite time. It is well-known that the system of Equation (1) has a global solution for every initial condition (s 0 , i 0 , r 0 ) with only nonnegative coordinates, and that solutions tend to equilibria of the form (s ∞ , 0, r ∞ ) (Odo Diekmann, H. Heesterbeek, and Britton 2012).

Works have extended these models to graphs, in order to increase their representative power ( Definition 1 (Diffusion Matrix). A diffusion matrix M is first Metzler, that is for i = j, we have M ij ≥ 0. Then, it is irreducible2 . Thirdly, its columns have vanishing sums.

The resulting reaction-diffusion dynamics is given by

             dS dt = -βS I + MS dI dt = βS I -δI + MI dR dt = δI + MR, (2)
where, as we explain below in Section 2.4, represents the coordinate-wise product3 . Standard results guarantee that the solution to Equation (2) is global, and converges to a fix point of the form (S, 0, R), as t → ∞ (Arino 2009). Moreover, the total population is preserved, that is n S n (t) + I n (t) + R n (t) is constant. Finally, standard Perron-Frobenius theory [START_REF] Meyer | Matrix Analysis and Applied Linear Algebra[END_REF] shows a diffusion matrix admits a unique stationary distribution, that is a positive vector μM summing to 1 such that Mμ M = 0. Moreover, solutions of dX/dt = MX with initial condition X 0 having a nonzero coordinate along μM converge to μM , as t → ∞. In particular, since the total population S + I + R satisfies this equation, and its initial condition has a nonzero coordinate along μM (one quickly checks it equals n S n (0)+I n (0)+R n (0), which is nonzero as S(0), I(0) and R(0) have nonnegative coordinates), it converges to this stationary distribution, as t → ∞.

Background on Issues in Network Reconstruction

The network reconstruction problem from observations, where one aims at expliciting the topology of a network of N nodes, by observing the values taken by some dynamical system which evolves on it, has been extensively studied in the literature (see for instance the review (Marc Timme and Casadiego 2014)). The network is described by some matrix M (typically, the adjacency matrix, possibly weighted). Observations are often gathered in two matrices, Ŷ and Ô, which typically belong to ∈ M N ×K (R), where K is the number of measurements. Often, Ŷ gathers estimates of the time derivatives of the state of the dynamical system in each node, at the different measurement times, and Ô is the so-called observation matrix gathering the values in each node, at the same times. Then, one knows the relation Ŷ = M Ô must be satisfied. Therefore, the problem amounts to solving this regression equation. We first describe the different observations possible, then address the solving of the regression.

Observations may first consist in measurements of the answer the system gives to some user-driven perturbation of its dynamics [START_REF] Gardner | Inferring Genetic Networks and Identifying Compound Mode of Action via Expression Profiling[END_REF]; Yeung, Tegnér, and Collins 2002; Yu and Parlitz 2010). In the case of non linear dynamics, these perturbations may occur near a fixed point, the interest being that the first-order expansion of this system then depends linearly on the network [START_REF] Gardner | Inferring Genetic Networks and Identifying Compound Mode of Action via Expression Profiling[END_REF], so that the observations Ŷ , Ô, and the network matrix M, indeed satisfy the Ŷ = M Ô equation. Alternatively, observations may be obtained through mere observation of the system (Shandilya and Marc For instance, in Bussel, Kriener, and M. Timme (2011), the authors use detailed knowledge of the evolution of a synthetic model of a biological synaptic network between spiking times, to obtain the matrices Ŷ and Ô verifying the Ŷ = M Ô equation. On the other hand, [START_REF] Casadiego | Model-free inference of direct network interactions from nonlinear collective dynamics[END_REF]) only assume some very general relation between the first order derivatives of the dynamical system, and the values it takes, in order to obtain similar relations.

Once obtained the observations such that the equation Ŷ = M Ô holds, one must then solve the regression problem. It may be over-determined, if K > N , or under determined, if K < N [START_REF] Stoer | Introduction to Numerical Analysis[END_REF]. Even if X has rank N , it may be ill-conditioned, thus preventing efficient solving by mere matrix inversion. To address these issues, a standard choice is to minimise the regression error with respect to some norm. One choice is then between L 1 or L 2 (least-squares) optimisation. The former induces sparsity, which may be desirable. For instance, [START_REF] Mangan | Inferring biological networks by sparse identification of nonlinear dynamics[END_REF] assume the dynamics decompose in some well-chosen basis, and that most of the coefficients in the expansion vanish. They then identify a subspace to which the vector of coefficients belongs, and finally use standard algorithms to find the sparsest vector in this subspace. In [START_REF] Yeung | Reverse engineering gene networks using singular value decomposition and robust regression[END_REF], the authors use an SVD decomposition of some observation matrix to parametrize the set of networks consistent with the data, and then use sparse regression to find the sparsest such network. In W.-X. [START_REF] Wang | Time-series based prediction of complex oscillator networks via compressive sensing[END_REF], the authors decompose the dynamics over some infinite basis, then use compressed sensing to evaluate the coefficients, only few of them are then nonzero.

Least-square optimization is on the other hand less costly, and better suited for over-determined systems. In Prasse and Piet Van Mieghem (2020), the authors use a least-square optimisation, but add a L 1 penalty in order to enforce some degree of sparsity, thus solving:

min M∈M N (R N ) Ŷ -M Ô 2 + λ M 1 ,
where λ is selected by cross-validation.

Finally, let us also mention [START_REF] Tyrcha | Network Inference with Hidden Units[END_REF], which in another vein differentiates the dynamics of the model, in order to train it to reproduce the observations, as is usual for Recurrent Neural Networks.

Contributions of the article

In our work, we study the reconstruction, and prediction, problems, for an epidemic reaction-diffusion. We assume known a model, and we consider that observations are a given, standalone time-series, which is arguably the harder setting observations wise, and which seems more relevant in the case of epidemic dynamics. We first show that for almost every network, the reconstruction problem is identifiable (Lemma 6). Then, we show that the quicker the diffusion, the lower the numerical rank of the observation matrix (Corollary 8), and the harder the reconstruction, but that increasing sampling helps reconstruct the network. Then, we classify symmetrical networks generating the same trajectories (Lemma 9). Finally, we show experimentally, on synthetic data constructed with random graph generators exhibiting different topologies, that reconstruction is easier for more "constrained" topologies, and that the prediction problem can still be solved satisfyingly even when the network topology makes exact reconstruction difficult. We use least-squares under constraints to solve numerically the reconstruction problem (see Section 3.3), in the experiments.

Notations and main definitions

We use the capital letter X to designate vectors on R 3N , for some integer N ≥ 1, which write X = (S, I, R), with S, I, R ∈ R N . Lower case x designates vectors on R 3 , with x = (s, i, r), and s, i and r real numbers. Whenever we consider some function f defined over R N , we choose to extend the notation in a straightforward way to R 3N , by writing, for X = (S, I, R) as above, f (X) := (f (S), f (I), f (R)). Whenever X is a function describing a trajectory of a dynamical system, X : R + → R 3N , we write Y = f (X) the function defined by, for all t ≥ 0, Y (t) := f (X(t)). As a result, combining with what precedes, for X = (S, I, R), for all t ≥ 0, we have Y (t) = (S(t), I(t), R(t)). For a set S ⊂ R N , for X = (S, I, R) ∈ R 3N , we write X ∈ S to mean that S, I and R belong to S. Then, for two vectors u and v of equal dimensions, we write u v their coordinate wise product: u v is the vector (u i v i ) i . Finally, 1 N designates the vector of R N will all coordinates equal to 1. Let us finally introduce the following three definitions, which help us formalise the setting.

Definition 2 (Model). We call model, and write M = (M, (β, δ) , X 0 ), a tuple consisting of a diffusion matrix M, epidemiological parameters gathered in β and δ, and an initial condition

X 0 = (S 0 , I 0 , R 0 ) ∈ R 3N .
Definition 3 (Flow on a model). A flow Φ on the set of models M is a mapping:

Φ : M → C R + , R 3N M → Φ (M)
defined by, for each model M = (M, (β, δ) , X 0 ), for all t ≥ 0, Φ t (M) is the value at time t of the solution of the differential Equation (2), with initial condition X 0 , that is

Φ t (M) = (S(t), I(t), R(t)).
Slightly abusing notations, in the following, we sometimes write Φ(M) when (β, δ) and X 0 are fixed, so that the model only depends on the choice of the diffusion matrix.

Let us finally introduce the observation matrix (Tyrcha and Hertz 2014), and the vectors of estimates of the reaction terms, and the derivatives. We do not observe the whole trajectories, but only some samples of them. For some integer K ≥ 1, let us then consider the sampling times 0 = t 0 < t 1 < t 2 < . . . < t K . Let for each node n Ŝn (t k ) be the (possibly noisy) observation of compartment S in node n at time t k , (and likewise for the other compartments). We can also estimate the vectors of derivatives, and of reaction terms, of Equation (2), from the observations, as is done in Shandilya and Marc Timme (2011). For every

1 ≤ k ≤ K, we define ρS (t k ) = β Ŝ(t k ) • Î(t k )
the vector of reaction terms on S at time t k , and

DS (t k ) = Ŝ(t k ) -Ŝ(t k-1 ) (t k -t k-1
) -1 the estimate of the derivative on S at time t k . We do likewise for the other compartments.

Definition 4 (Observation matrix, derivatives and reaction terms). Let the observation matrix on S be

ÔS =    Ŝ1 (t 1 ) . . . Ŝ1 (t K ) . . . . . . . . . ŜN (t 1 ) . . . ŜN (t K )    ∈ M N ×K (R) .
Note likewise ÔI the observations on I, and ÔR those on R. Define finally the matrix by block

Ô ((t k ) , Φ(M)) = ÔS , ÔI , ÔR ∈ M N ×3K (R)
. This is the observation matrix associated with the sampling times (t k ), and the flow Φ(M).

Likewise, we write ρS , DS ∈ M N ×K (R) the matrices of reaction terms (resp. derivatives) on S, and likewise for the other compartments. We finally define ρ = (ρ S , ρI , ρR ) ∈ M N ×3K (R) and D = DS , DI , DR ∈ M N ×3K (R).

Problem Studied, Identifiability and Experimental Setup

We present the problem we address (Section 3.1), then initate the study of its identifiability (Section 3.2), and finally present the setting in which we conduct the experiments (Section 3.3).

Problematic: Reconstruction and Prediction

We assume known the initial condition X 0 , and the epidemiological parameters β n 's and δ n 's. Some unknown diffusion matrix M * then generates a flow (S, I, R) = Φ(M * ), and given as observations the trajectories S, I and R, we address the following two questions.

1. Question 1. Can we estimate M * from the observations? In other words, do they first uniquely define M * ? And in so, is it possible to estimate M * from them?

2. Question 2. Can we predict the future evolution of the system if, for some T train < ∞, we can only observe the trajectories in some initial phase 0 ≤ t ≤ T train of the system, that is the observations only consist in (Φ t (M * ), t ≤ T train )?

Identifiability of the Diffusion

We now conduct some preliminary analysis on the identifiability of the diffusion.

We have the following characterisation of the set of diffusion matrices M which produce the same trajectories as M * (see Appendix A for a proof, and likewise for future results).

Lemma 5 (Diffusions Generating the Same Trajectories). Let M * be a diffusion matrix, and write (S, I, R) = Φ (M * ). Then, every matrix M = M * + H, such that first M is a diffusion matrix, and secondly such that for all t ≥ 0, we have4 Φ t (M) ∈ ker H, produces the same trajectories as M * .

As a result, provided the vector space generated by the trajectories, that is by the vectors S(t), I(t) and R(t), for t ≥ 0, is the whole space R N , then the answer to our first question is affirmative (as the only H possible vanishes over the whole space, therefore vanishes). Therefore, a fundamental question governing the issue of the identifiability of the diffusion matrix is the existence of strict subspaces of R N in which the trajectories evolve. This moreover gives us a practical criterion to evaluate if the diffusion matrix generating a given trajectory is unique: we check if the observation matrix has rank N , which is sufficient to guarantee the uniqueness. Now, often the trajectories generate the whole space, as the next result shows.

Lemma 6 (Almost Everywhere Identifiability). Let 0 ≤ t 1 < ... < t N < ∞ be a subdivision of the nonnegative real half-axis. Then, for almost every M, X 0 , for all β, δ, writing M = (M, (β, δ) , X 0 ), the space generated by the samples of the trajectories at instants t 1 , ..., t N (that is

Φ t1 (M), . . . , Φ t N (M)), is equal to R N .
This might give the impression the problem is solved, for almost every M and X 0 . Indeed, assume the space generated by the trajectories is the whole of R N . Then, the observation matrix Ô ((t k ) , Φ(M)) has rank N . We therefore know the image of M on a basis, which fully determines it. However, the conditioning of the observation matrix is often poor in practise, so that reconstruction of M * by extracting a basis is inefficient. In the next two sections, we investigate two reasons why this is the case. Firstly, we study the influence of the speed of diffusion in Section 4. Secondly, we consider the topology of the graph in Section 5. Nonetheless, the fact the diffusion is often unique means that, when running a reconstruction algorithm, we can have good hope it will succeed in finding a good fit, which we show is the case in the experiments.

Experimental Set-Up

For each experiment, we start by generating a graph from a random graph generator. We use four random graphs, with different topologies: the Erdős-Rényi and the Waxman graph, which are quite connected graphs, and therefore quite "close" to a complete graph, and the Relaxed Caveman and extended Barabási-Albert graphs, which are less connected, and exhibit a more clustered structure.

In that sense, the diffusion is more constrained by these graphs, and we expect the reconstruction problems to be easier in that case. We study graphs of sizes ranging from N = 20 to N = 140. Then, we draw the coefficients of the diffusion matrix M * uniformly on [0, 1]. Finally, when needed, we renormalise the diffusion matrix by the typical time of diffusion τ > 0. The epidemiological coefficients β n and δ n are drawn at random, such that for each node n, we have

β n ∼ |N (1, 2)| and δ n ∼ N 3 × 10 -2 , 6 × 10 -2 .
Next, we simulate the ground truth trajectories on the time interval [0, T ∞ ], with T ∞ = 10. We use a uniform time discretisation step of (∆t) groundtruth = 10 -3 , and a Runge-Kutta discretisation scheme of order 4. For each setting, we repeat the experiments at least 5 times, so as to control the stochastic fluctuations. Finally, the train set in which we sample the observations is [0, T train ], with T train = 2. We use various sample steps (∆t) sample , ranging between 2 × 10 -3 and 10 -2 . Recall from Section 2.1 that we write μM the stationary distribution of the diffusion matrix M. In each experiment, we use as initial condition X 0 = (S(0), I(0), R(0)) = (s 0 μM , i 0 μM , r 0 μM ), where s 0 , i 0 and r 0 are nonnegative real numbers, and i 0 > 0. As a result, the vector of initial susceptibles S(0) is proportional to the stationary distribution, and likewise for I(0) and R(0). We compute the reconstructed diffusion matrix by solving

min M∈M N (R) D -R -M Ô 2 , such that      M i,j ≥ 0, i = j, j M i,j = 0 for all node i.
This is a convex optimisation problem. We solved it using the Python package CVXPY (Diamond and Boyd 2016; [START_REF] Agrawal | A rewriting system for convex optimization problems[END_REF]. We write M rec the matrix obtained. Moreover, to truly enforce the fact M rec is a diffusion matrix, we postprocessed the matrix obtained by enforcing that column sums vanish: for every node n, we replaced the diagonal coefficient

M rec (n, n) by - N i=1 M rec (i, n).
To assess the reconstruction, we use two metrics. First, we use the AUC [START_REF] Fawcett | An introduction to ROC analysis[END_REF]) on the presence of edges, as Prasse and Piet Van Mieghem (2020). It is computed thanks to the corresponding fonction in Scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine Learning in Python[END_REF]. Secondly, we evaluate the prediction error, that is the norm of the difference between the trajectories computed with the true model, and those computed with the reconstructed diffusion M rec , by

1 N 1 T ∞ -T train pmax p=1 Φ tp (M * ) -Φ tp (M rec ) 2 (t p -t p-1 ) , (3)
where p max = T∞-Ttrain (∆t) grountruth and (t p ) is the discretisation scheme used for the simulations, whose beginning has been removed, so that t 0 = T train .

Computations with the Barabási-Albert graph were proner to numerical instabilities. We believe this is due to its topology being more constrained. As a result, we modified a bit the experimental setting for this specific graph, increasing sampling to (∆t) sample = 4 × 10 -3 , and increasing the number of repetitions to 30.

Influence of the Diffusion Rate

We now study the influence of the diffusion rate on the feasibility of the network reconstruction, first theoretically (Section 4.1), then experimentally (Section 4.2).

Analysis

One difficulty of the network reconstruction problem is the conditioning of the observation matrix (Definition 4), which may be poor. In particular, its numerical rank may be significantly lower than N , as observed also in (Prasse and Piet Van Mieghem 2020). In our case, this may be partly due to the homogenisation performed by the diffusion. Indeed, given different epidemiological parameters, and different population sizes, the internal dynamics of the different nodes evolve differently. However, the diffusion tends to homogeneise each compartment, so that S(t) tends to a vector proportional to the stationary distribution, μM , and likewise for I(t) and R(t). As a result, the diffusion tends to worsen the conditioning of a basis. This effect depends on the time-scale at which diffusion occurs, with respect to that at which the reactions in each node occur. We first show, in the following Lemma 7, that when the typical time of evolution of the diffusion, τ , goes to 0 (equivalently, the diffusion rate 1/τ goes to infinity), and in the presence of fixed epidemiological parameters, the trajectories tend to those of a scalar SIR systems, which coefficients we express in terms of the β n 's, the δ n 's and the stationary 2). Lemma 7 (Limit Trajectories for Diffusion Rate going to Infinity). Let M = (M, (β, δ) , X 0 ) be a model, and assume the initial condition X 0 is such that S 0 , I 0 and R 0 are proportional to the stationary distribution μM . Write (s, i, r) the solutions of the scalar system

             ds dt = -βsi di dt = βsi -δi dr dt = δi,
with s(0) = n S n (0), and likewise for i and r, and with

β = n β n μM (n) 2 , and δ = n δ n μM (n).
Then, for any T > 0, Φ τ (M) → (sμ M , iμ M , r μM ), as → 0, uniformly on [0, T ].

We illustrate Lemma 7 on Figure 1. We ran experiments according to the protocol described in Section 3.3, using Erdős-Rényi and Relaxed Caveman graphs, for a range of values of τ . For (t p ) the discretisation scheme used for the simulations, and p max the number of t p 's, we plot the error

1 N 1 T ∞ pmax p=1 Φ tp ( M τ ) -(s(t p )μ M , i(t p )μ M , r(t p )μ M ) 2 2 (t p -t p-1 )
between the trajectories obtained with the vector model, and those computed from the scalar model. We indeed see it goes to 0, as τ → 0. Moreover, the discrepancy is bigger for the Relaxed Caveman graph, than for the Erdős-Rényi one: indeed, the latter is more connected, and therefore there are much more exchanges between the nodes, so that it is closer to a kind of "average" model, which the scalar limit is. Then, from Lemma 7, we immediately have the following corollary which describes its consequences for the numerical rank of the observation matrix.

Corollary 8 (Numerical Rank of the Observation Matrix for Diffusion Rate going to Infinity). We make the same assumptions as in Lemma 7. Let, for some integer K ≥ 1, (t k ) 1≤k≤K be a family of sample times. Let τ > 0, and let us We illustrate this convergence on Figure 2. We use the same protocol as for Figure 1, but this time display the numerical rank. We see it gets lower and lower, as τ → 0. It is lower for the Erdős-Rényi graph, probably for the same reasons given above.

Experiments: Diffusion Rate, Sampling Frequency

Let us now investigate the consequences of this phenomenon, for the practical reconstruction problem. We first show on Figure 3 the AUC as a function of the typical time of diffusion τ , for a fixed sampling rate. The AUC increases as the typical time of diffusion decreases, as we expected. It is bigger for the Relaxed Caveman graph, which has "more structure" than the Erdős-Rényi one. Then, we study how increased sampling may help prediction for high diffusion rates. We therefore ran experiments for different values of τ , and different sampling rates.

On Figure 4, we show a heatmap of the AUC, with different sampling steps, and diffusion rates, for an Erdős-Rényi graph of 30 nodes. The darker the color, the smaller the AUC is. On each row, we see colors get darker as we go to the right: this means that, for each fixed diffusion rate, the AUC deteriorates as the sampling step increases. On each column, we see colors get darker as we move to the top: this means that, for each sampling step, the AUC worsens as the diffusion rate increases. Overall, we see that the bottom left triangle is lighter (sampling is high enough with respect to the diffusion rate, AUCs are big), while the top right triangle is darker (sampling is low with respect to the diffusion rate, AUCs are lower).

Influence of the Network Topology

We now study the influence of the network topology, first theoretically, from an algebraic standpoint (Section 5.1), then experimentally (Section 5.2).

Symmetries

Thanks to Lemma 5 we know that identifiability of M is linked to the dimension of the vector space spanned by the flow Φ(M). Now, symmetries of the model may cause the trajectories to live in low dimensional spaces. Indeed, they often lower the dimensions of the studied spaces by eliminating the dependencies of equations in some variables. For instance, a 2 dimensional problem in physics which is invariant under rotations around the origin will have a solution which will only depend on the distance to the origin. These principles have been applied successfully to numerous fields, and have been used in the context of mathematical epidemiology (Simon, Taylor, and Kiss 2011; Ward and López-García 2019) to reduce the number of calculations needed to simulate the propagation of diseases. We now investigate the influence of symmetries on the inverse problem we study.

Let us first define precisely symmetries. We write S N the symmetric group of order N , and σ its elements, which are called permutations. We write P (σ) the permutation matrix associated to the permutation σ. Then, for any V ∈ R N , P (σ)V = V if, and only if, for every orbit of σ, for every i, j in this orbit, we have V i = V j . A vector X = (S, I, R) is symmetric with respect to σ if P (σ)S = S, and likewise for I and R. This notion extends to groups of permutations, as follows.

Let H be a subgroup of S N , and define Fix(H) (Lang 2012) the space of vectors stable by H, that is:

Fix(H) = V ∈ R N | ∀σ ∈ H, P (σ)V = V .
Then, X is said to be symmetric with respect to H if X ∈ Fix(H). This extends to flows by saying that a flow (S,

I, R) = Φ(M) is symmetric with respect to H if Φ(M) ∈ Fix(H).
Finally, we say that a model M = (M, (β, δ) , X 0 ) is symmetric with respect to some permutation σ if, writing P = P (σ), we have M = P MP -1 , P β = β, P δ = δ and P X 0 = X 0 . We then say σ is an automorphism of M, extending in a straigthforward way the notion of graph automorphism [START_REF] Hell | Graphs and Homomorphisms[END_REF]). Indeed, if σ is an automorphism of M, then it is in particular an automorphism of the underlying weighted graph, meaning that for all nodes i, j ∈ N , the edges i j and σ(i) σ(j) have the same weight: M i,j = M σ(i),σ(j) . We write Aut(M) the group of model automorphisms of M.

We first establish, in Lemma 9, that trajectories generated by a diffusion M are symmetrical with respect to some group H first if, and only if, there exists a diffusion admitting all permutations in H as automorphisms which generates the same trajectories and secondly if, and only if, the diffusion M stabilizes Fix(H).

Lemma 9 (Networks Generating Symmetrical Trajectories). Let M be a diffusion matrix, β, δ be the vectors of epidemiological coefficients, and H be a subgroup of S N . Assume that β and δ are symmetric with respect to H. Then, the following conditions are equivalent.

1. Symmetries of the Trajectories. For all X 0 ∈ Fix(H), the flow of (M, (β, δ) , X 0 ) is symmetric with respect to H.

Symmetrical Generating Diffusion.

There exists a diffusion matrix M such that H ⊂ Aut M, (β, δ), X 0 and such that for all X 0 ∈ Fix(H), the flow of M, (β, δ), X 0 equals the flow of (M, (β, δ), X 0 ).

Stabilization by the Diffusion. M stabilizes Fix(H), that is MFix(H) ⊂ Fix(H).

We now show that, in the spirit of Lemma 5, trajectories symmetrical with respect to H are generated by diffusions which differ by a matrix Z vanishing on Fix(H). These matrices represent the fact flows between nodes with identical S, I and R values may be redirected freely within themselves, provided the outgoing flows are modified accordingly5 . The nodes where the flows are identical are those in the same orbits under H (Lang 2012), that is the nodes i and j such that, for some σ ∈ H, we have j = σ(i). Define therefore, for all 1 ≤ i < j ≤ N , and for all 1 ≤ k ≤ N -1, the redirection matrix

Z i,j,k = E k,i -E k,j -E N,i + E N,j ,
where the E r,s 's matrices are the vectors of the canonical basis of M N (R). This matrix removes one unit of rate from the edge j k, and adds one unit of rate on the edge i k. It does the reverse with respect to the node N , taking one unit of rate from i N and adding it to j N , in order to enforce the fact that the sums of Z i,j,k vanish, that is as much rate goes to each node than goes out.

Lemma 10 (Flow Redirection within the Orbits). Under the same assumptions as in Lemma 9, let H be the biggest group of symmetries letting invariant the trajectories. Then, the affine space of matrices producing the same trajectories as M for every initial condition X 0 ∈ Fix(H) is exactly the subspace generated by the Z i,j,k 's, for all i and j which are in the same orbit under H. This space has dimension at least

(N -1)(N -#{different trajectories}).
The dimension of this space is a lower bound on the dimension of the affine space of matrices generating the same trajectories as M. To summarize, given a diffusion matrix M, we have given an explicit description of a set of matrices giving the same trajectories as M. As a result, if the diffusion matrix we try to reconstruct gives symmetrical trajectories, and if we have an algorithm which gives us one solution of the reconstruction problem, then we are able to find many such matrices explicitly, though we cannot single the original M out. Note that this has consequences on the conditioning of the observation matrix. Indeed, its rank is then necessarily bounded by #{different trajectories}. As such, if the model presents symmetries, then several singular values of the observation matrix will be zero, and in a neighbourhood of M as well, the numerical rank will be bounded by N -#{different trajectories}. This proves that the nearest a model is to a symmetrical model, the most difficult it is to reconstruct the diffusion matrix.

Network Topology Experiments

We now study experimentally the influence of the network topology on the estimation and prediction problems. As explained in Section 3.3, we study two metrics: the AUC on the presence or absence of edges, and the prediction error, and we present results for various sizes of graphs, and various types of random networks, exhibiting different topologies. The plots are box plots, where the solid lines are the medians of values, and the shaded areas gather the [10%, 90%] intervals of values.

We present, on Figure 5, the AUC as a function of the number of nodes, for various types of random graphs. The full set of parameters used is available in the code. The AUC is quite good for small graphs, more than 0.8, but decreases as the number of nodes increases. As expected, the more constrained the topology, the better the AUC. Indeed, it is in general best for the Barabási-Albert graph, and the second best is often the Relaxed Caveman graph. The Waxman gaph, and above all the Erdős-Rényi one, exhibit the worse AUCs.

On Figure 6, we show the prediction error. We see the prediction errors in general are quite low, less than 2×10 -4 , and diminish with the number of nodes. Moreover, the Erdős-Rényi graph consistently exhibits the lowest error. These results are consistent with each other, in the sense that it seems the more the graph has connections, the easiest it is to predict the future behaviour of the system (more edges, either through more nodes, or through the topology, in the case of the Erdős-Rényi graph). However, they are opposite to the results for the AUCs. Therefore, they tend to suggest that the more constrained the topology, the easier it is to reconstruct the network, but the more mixing there is, the easiest it is to predict the future evolution of the system. We did not display the prediction errors for the Barabási-Albert graph, as it was about 5 times greater than for the other graphs, and exhibited also high variance. We believe it comes probably first from the fact it is proner to numerical instabilities, as we said in Section 3.3. Secondly, it is also probably due to its topology being more constrained: as a result, small errors on the reconstruction lead to high errors on the prediction. Finally, on Figure 7, we see the numerical rank of the observation matrix for several graphs. It tends to stagnate or decrease as the number of nodes increases, which is not surprising, as large matrices tend to have small singular values, which therefore do not contribute to the numerical rank. It is consistently higher for the Barabási-Albert graph, which structure is more constrained.

Conclusions, future works

In this article, we have studied the reconstruction, and prediction, problems, for an epidemic reaction-diffusion. We have proved that for almost every network, the reconstruction problem is identifiable. Then, we have shown that the quicker the diffusion, the lower the numerical rank of the observation matrix, and the harder the reconstruction, but that increasing sampling helped reconstruct the network. Then, we have classified symmetrical networks generating the same trajectories. Finally, we showed experimentally, on synthetic data constructed with random graph generators exhibiting different topologies, that reconstruction was easier for more "constrained" topologies, and that the prediction problem could still be solved satisfyingly even when the network topology makes exact reconstruction difficult.

We have studied the case when the observations we have consist of all the trajectories in all the nodes. Other studies have considered partial observations, as in the presence of missing nodes [START_REF] Tyrcha | Network Inference with Hidden Units[END_REF][START_REF] Haehne | Detecting Hidden Units and Network Size from Perceptible Dynamics[END_REF], or more generally partial observations (Nitzan, Casadiego, and TImme 2017; Ioannidis, Romero, and Giannakis 2018; Ioannidis, Shen, and Giannakis 2019). This would be an interesting extension to further our work.

Then, measures often come with a level of noise, and studying the robustness of reconstruction in the presence of noise should be another direction of study.

Finally, studies on real-world data would allow to attempt at separating the specific influence of diffusion, through transportation networks for instance, on the spread of epidemics, from that of internal (within cities, or countries) reactions.

A Proofs for Section 3. Likewise, I and R satisfy the corresponding equations, with diffusion matrix M instead of M * . Now, (S, I, R) satisfies the differential equation Equation (2) with diffusion matrix M, and starts at the initial condition X 0 . Therefore, by unicity of the solutions of this system, we have (S, I, R) = Φ(M).

To prove Lemma 6, we first prove that the trajectories generated by a linear system often generate the whole space. Before doing this, we need the following technical lemmas.

Lemma 11 (Technical Result). For every integer N ≥ 1, let (λ n ) 1≤n≤N and (t k ) 1≤k≤N be two families of distinct (within each family) real numbers. Then, for every N ≥ 1, the following property holds true:

P(N ) : ∀1 ≤ n ≤ N, N k=1 µ k e λnt k = 0 ⇒ (∀1 ≤ k ≤ N, µ k = 0).
Proof. We prove it by induction on N . The case N = 1 is immediate. Let N ≥ 2, and let us show that if P(N -1) holds, then P(N ) holds as well. Assume that for all 1 ≤ n ≤ N , we have N k=1 µ k e λnt k = 0. Then, for all 1 ≤ n ≤ N , we have

N -1 k=1 µ k e λn(t k -t N ) + µ N = 0.
Let us consider the mapping:

g : R λ → N -1 k=1 µ k e λ(t k -t N ) + µ N .
Then, g has N distinct roots (the λ n 's), therefore by Rolle's theorem, its derivative admits N -1 distinct roots. We therefore obtain N -1 values λ n such that, for all 1 ≤ n ≤ N -1, λ n < λ n < λ n+1 and dg dλ (λ n ) = 0. As a result, for all 1 ≤ n ≤ N -1, we have

0 = dg dλ (λ n ) = N -1 k=1 µ k (t k -t N )e λn(t k -t N ) .
By induction, for all 1 ≤ k ≤ N -1, we have µ k (t k -t N ) = 0, and therefore µ k = 0, as the t k 's are distinct. Then, µ N = 0 as well, hence the result.

We need the following lemma. Though well-known, we could not locate a reference, so added it here for completeness.

Lemma 12 (Almost Every Diffusion Matrices has Distinct Eigenvalues). Almost every diffusion matrix has N distinct eigenvalues.

Proof. First, for all τ > 0, and t ≥ 0, we have

S τ (t) = n S τ,n (t) μM + ν τ (t), (5)
where ν τ (t) belongs to the set H = ν ∈ R N | n ν n = 06 . Moreover, ν τ (t) is bounded uniformly in τ > 0 and t ≥ 0, as all the S τ (t)'s are bounded by the total population, that is the sum of the coordinates of the initial condition X 0 .

Then, for all τ > 0, t → ν τ (t) is differentiable thanks to Equation (5) and, by differentiating Equation (5), we see ν τ satisfies a differential equation of the form:

dν τ dt (t) = M τ ν τ (t) + γ τ (t),
where γ τ is a quantity depending on many things, but which is uniformly bounded in τ > 0 and t ≥ 0, again thanks to the fact that the S τ (t)'s, the I τ (t)'s and the R τ (t)'s are bounded by the total size of the population. Moreover, for all ν ∈ R N , we have Mν ∈ H, since the columns of M sum to 0. As a result, for all τ > 0, and all t ≥ 0, we have dντ dt (t) ∈ H (since H is a finite dimensional vector space, so derivatives of functions living on it stay in it), and M τ ν τ (t) ∈ H (by what precedes), so that γ τ (t) also belongs to H.

Finally, for all τ > 0, and all t ≥ 0, we have

ν τ (t) = exp M t τ ν τ (0) + t 0 exp M t -s τ γ τ (s) ds = t 0 exp M t -s τ γ τ (s) ds,
since ν τ (0) = 0, as the initial condition is proportional to μM by assumption. Let us now fix ε > 0. Let B ⊂ H be a ball such that that, for all τ > 0 and t ≥ 0, we have γ τ (t) ∈ B. Since M only has eigenvalues with (strictly) negative eigenvalues on H, we may find some threshold u min > 0 such that, for all u ≥ u min , for all ν ∈ B, we have exp (Mu) ν < ε.

Moreover, there exists a constant κ ≥ 1 such that, for all u ≥ 0, for all ν ∈ B, we have exp (Mu) ν ≤ κ. Let us now consider τ ≤ 1 umin κ ε (it is chosen so that, for t -s ≥ ε κ , we have t-s τ ≥ u min ). Therefore, for all t ≤ T , we have

ν τ (t) ≤ t 0 exp M t -s τ γ τ (s) ds = ε κ 0 exp M t -s τ γ τ (s) ds + T ε κ exp M t -s τ γ τ (s) ds ≤ ε κ κ + T ε κ ε ds ≤ ε + ε T - ε κ ≤ ε (1 + T ) .
As a result, for all τ ≤ 1 umin κ ε , with u min and κ chosen independently of τ , we have sup 0≤t≤T ν τ (t) ≤ ε (1 + T ). We have therefore proven that d (S τ (t), Rμ M ) → 0, as τ → 0. We would prove likewise the result for I τ and R τ , which concludes the proof.

We can now prove Lemma 7.

Proof. Let T > 0. Now, let us write, for all t ≥ 0, using the notations of the proof of Lemma 14,

S τ (t) = s τ (t)μ M + ν τ (t), (6)
and likewise for I τ and R τ . Since S 0 , I 0 and R 0 are proportional to the stationary distribution, we know, thanks to Lemma 14, that ν τ (t) tends towards 0, uniformly on each [0, T ], with T ≥ 0. Since for all τ , S τ , I τ and R τ are nonnegative, and bounded by the total population, the family of functions (s τ , i τ , r τ ), defined on [0, T ], has values in a bounded set of the continuous functions from R + to R 3 , endowed with the infinity norm on each of s, i and r, that is (s, i, r) ∞ = max ( s ∞ , i ∞ , r ∞ ). Moreover since, for all τ ≥ 0, for all 0 ≤ t ≤ T , we have

S τ (t) = S τ (0) - t 0 βS τ (s) I τ (s) + M τ S τ (s) ds,
we obtain, summing along the coordinates,

s τ (t) = s τ (0) - n β n μM (n) 2 t 0 s τ (s)i τ (s) ds.
As a result, s τ is differential on [0, T ], its derivative satisfies

ds τ dt = - n β n μM (n) 2 s τ (t)i τ (t),
and its derivative is threfore bounded on [0, T ], uniformly on τ . The same holds for i τ and r τ . Therefore, (s τ , i τ , r τ ) is equi-continuous. As a result, the family ((s τ , i τ , r τ )) τ ≥0 is pre-compact (Sutherland 2004) in the Banach space of functions from [0, T ] to R 3 , endowed with the infinity norm defined above, so that, provided it admits an unique adherence value, it converges towards this one.

Let us consider a converging subsequence, and still index it by τ , to simplify notations. As a result, the limit s satisfies

s(t) = s(0) + n β n μM (n) 2 s 0 s(s)i(s) ds,
and likewise for i and s. Moreover, for all t ≥ 0, we know that s τ (0) = n S τ (0) = n S n (0) which does not depend on τ , therefore s(0) = n S n (0), and likewise for i 0 and r 0 . Therefore, (s, i, r) is solution of the scalar system described in the statement of the Lemma, and by uniqueness of the solutions of this system, satisfying the initial condition (s 0 , i 0 , r 0 ), the tuple is uniquely defined. Therefore, the family of (s τ , i τ , r τ )'s admits a unique adherence value, and converges towards this one. Plugging back into Equation (6), we see that S τ → sμ M , as τ → 0, uniformy on [0, T ], and likewise for I τ and R τ , which concludes the proof.

Let us prove Corollary 8.

Proof. Applying the results of Lemma 7 with T = t K , we know that the observation matrix writes Ô ((t k ) , Φ τ (M)) = μM ⊗ (s(t 1 ), . . . , s(t K ), i(t 1 ), . . . i(t K ), r(t 1 ) . . . r(t K )) + O (ε(τ )) , (7) as τ → 0, where ε(τ ) → 0, when τ → 0. Indeed, using the notations of the proof of Lemma 7, we know that for each 1 ≤ k ≤ K, the column S(t k ) of the observation matrix (for instance), writes S(t k ) = s(t k )μ M + ν τ (t k ), and ν τ tends to 0, as τ → 0, uniformly on [0, t K ]. The first term of Equation ( 7) is of rank 1, as s + i + r = 1, identically. The conclusion follows from the continuity of the numerical rank of a matrix (the numerical rank is the number of singular values greater than some threshold, and these values depend continuously on the matrix).

C Proofs for Section 5.1, symmetries

Let us first show the effect of node-renumbering on the trajectories.

Lemma 15 (Node re-numbering). Let M = (M, (β, δ) , X 0 ) be a model and P a permutation matrix. Then, for all t ≥ 0, we have Φ t (P • M) = P • Φ t (M) = (P S(t), P I(t), P R(t)) , writing (S(t), I(t), R(t)) = Φ t (M), and using the conventions of Section 2.4. This implies immediately that if P is the matrix of an automorphism of our model, then for all t ≥ 0, we have Φ t (P • M) = Φ t (M), and therefore P • Φ t (M) = Φ t (M). In other words, if i and j are in a same orbit of P , then the trajectories at nodes i and j are the same: for all t ≥ 0, S i (t) = S j (t), and likewise for I and R.

Proof. Let us show that (P S, P I, P R) is a solution of the differential equation also satisfied by the flow Φ P MP -1 , (P β, P δ) , P X 0 , which is enough to conclude by unicity of the solutions sharing the same initial condition. Let 1 ≤ i ≤ N . Then, we have So P S satisfies the required equation, and I, R do as well, which we show using the same method, and which allows us to conclude.

We can now prove Lemma 9.

Proof. Let us first prove 3) ⇒ 2). Let σ ∈ H. For any x ∈ H, we have P (σ)MP (σ) -1 x = P (σ)Mx as x ∈ Fix(H) = Mx as Mx ∈ Fix(H).

We now average the nodes which give the same trajectories, which is the standard method of averaging under a group action. Let We therefore obtain by construction that, for all X 0 ∈ Fix(H), we have H ⊂ Aut( M, (β, δ), X 0 ). The fact M and M agree on Fix(H) is a direct consequence of the assumptions of the lemma. The Z i,j,k 's are clearly linearly independent. Let us show that they generate the space of all Z's.

Let Z vanish on F ix(H). We will make all of the coefficients of Z vanish by substracting multiples of Z i,j,k s, which will prove that Z is indeed a linear combination of the Z i,j,k s. Let O = {i 1 < ... < i m } be an orbit of {1, ..., n} of cardinal m under the action of H.

Let us remark that for 1 ≤ k ≤ N -1, 1 ≤ l < m, Z i l ,i l+1 ,k satisfies the conditions of the lemma and has its (k, i l ) coefficient equal to 1, its (k, i l+1 ) coefficient equal to -1.

Therefore, Z -Z k,i1 Z i1,i2,k has its (k, i 1 ) coefficient equaling zero, and its (k, i 2 ) coefficient equal to Z k,i2 + Z k,i1 , and aside from the last line (which we will ignore for the moment) these are the only coefficients changing.

Then, if m ≥ 3, we can reiterate this by considering Z -Z k,i1 Z i1,i2,k -(Z k,i2 + Z i1,k )Z i2,i3,k and the obtained matrix will have the (k, i 2 ) coefficient vanishing and the (k, i 3 ) coefficient changing to Z k,i3 + Z k,i2 + Z k,i1 , and the (k, i 1 ) coefficient is still 0.

We iterate this method exactly m -1 times to obtain Z . By construction, Z has each of the (k, i l ), 1 ≤ l ≤ m coefficients vanishing except maybe the l = m one, equaling Z k,i1 + ... + Z k,im .

But then this one is also zero. Indeed, if (e j ) j is the canonical basis, as m l=1 e i l ∈ F ix(H), we have Z m l=1 e i l = 0, and by looking the k-th coefficient, we obtain Z k,im = Z k,i1 + ... + Z k,im = 0.

We can then iterate this construction on every line except the last (meaning for 1 ≤ k ≤ N -1) and every orbit to obtain Z . By construction, every line of Z is zero, except maybe the last (k = N ), but then as the columns of Z have a vanishing sum (as Z is a linear combination of Z and the Z i,j,k ), Z = 0. Thus, Z is in the space generated by the Z i,j,k .

  Kiss 2011; Ward and López-García 2019), and experimentally. Similar questions have been asked for other epidemic models (Ganesh, Massoulié, and Towsley 2005; Durrett 2010; Vajdi and Caterina Scoglio 2018; Prasse, Devriendt, and Piet Van Mieghem 2021).

Figure 1 :

 1 Figure 1: Discrepancy between the Scalar Model and the True Model as a Function of the Typical Time of Diffusion

Figure 2 :

 2 Figure 2: Numerical rank, of the observation matrix, for a fixed family of sample times, and two random graphs of 30 nodes, for various diffusion rates τ .

Figure 3 :

 3 Figure 3: AUCs for Various Diffusion Rates

Figure 4 :

 4 Figure 4: Heatmap of AUCs wrt. Typical Time of Diffusion and Sampling Step

Figure 5 :

 5 Figure 5: AUC on the Adjacency Matrix for Various Graphs and Sizes of Graphs

Figure 6 :

 6 Figure 6: Prediction Error on the Trajectories for Various Sizes of Graphs

=

  dS σ(i) dt = -β σ(i) S σ(i) I σ(i) + [M S] σ(i) = -[P β] i [P S] i [P I] i + [P M S] i = -[P β] i [P S] i [P I] i + [(P M P -1 )(P S)] i .Therefore, we have dP S dt = -(P β) ((P S) (P I)) + P MP -1 (P S) .

P

  (σ)MP (σ) -1 .

  Driessche 2001; Van den Driessche and Watmough 2002; W. Wang and Zhao 2005; Allen et al. 2007; Tien et al. 2015;

  2, Identifiability of the DiffusionLemma 5. Let M = M * + H be a diffusion matrix, such that for all t ≥ 0, we have S(t) ∈ ker(H), I(t) ∈ ker(H) and R(t) ∈ ker(H). Let us check (S, I, R) = Φ(M * ) is then a solution of Equation (2) with diffusion matrix equal to M. Indeed, for all t ≥ 0, we have

	dS dt	(t) = -βS(t) I(t) + M * S(t) by assumption
		= -βS(t) I(t) + M * S(t) + HS(t) as S(t) ∈ ker H
		= -βS(t) I(t) + MS(t).

We adopt this terminology, for lack of a universally agreed term for these matrices.

This is possible if the graph is directed because we ask it is then strongly connected.

For instance, for a node n, the equation on Sn reads: dSn/dt = -βnSn(t)In(t) + N i=1 M n,i S i (t). Our notation is not standard, but we think it makes clearer the link of the graph system of Equation (2) with the scalar system of Equation (1).

As explained in Section 2.4, Φt(M) ∈ ker H means that S(t), I(t) and R(t) belong to ker H.

In fact, Z matrices, like the diffusion matrices, describe rates. However, as long as nodes have equal values, modifying the rates, or the flows going out of them, becomes equivalent.

This is a consequence of the decomposition R N = Rμ M ⊕ H.

Proof. Let E be the subset of diffusion matrices which do not have N distinct eigenvalues, and let us prove it has null measure. We know that E is the zero set of the discriminant ∆ applied to the characteristic polynomial [START_REF] Lang | Algebra. Graduate Texts in Mathematics[END_REF], defined for every matrix A of order N by

The set of diffusion matrices is contained in the linear subspace of matrices A satisfying 1 T N A = 0, which is of dimension N 2 -N . We equip it with the standard Lebesgue measure over R N 2 -N .

Then, for any diffusion matrix D ∈ E, and any diffusion matrix M, D + λM is in E only for a finite number of λ ∈ R. Indeed, λ → F (D+λM) = ∆•det(XId-D-λM) is a polynomial function, so it is either zero or has a finite number of zeroes. But then it is nonzero at λ = 0, because D ∈ E. Therefore, the indicator function χ of E is zero almost everywhere on {D + λM, λ ∈ R}, so it is also zero almost everywhere on {M + λD, λ ∈ R}.

Let us finally fix some D / ∈ E, and some vector space H such that H ⊕ Vect(D) = {A, 1 T N A = 0}. Thanks to the Fubini-Tonelli theorem, we have

Therefore, E has zero Lebesgue measure.

Lemma 13 (Almost Everywhere Identifiability, linear case). Let 0 ≤ t 1 < ... < t N < ∞ be a subdivision of the nonnegative real half-axis, M ∈ M N (R), and y 0 ∈ R N . Let y be the solution of the following differential equation:

Then, (y(t 1 ), ..., y(t N )) is a basis of R N for almost every M, y 0 .

Proof. Let us prove the result for M and y 0 satisfying the additional assumptions that all the eigenvalues of M have multiplicity 1, and that every coordinate of y 0 in an eigenbasis of M is nonzero (note that as M has N distinct eigenvalues, the associated subspaces are 1-dimensional, so the eigenbasis is unique up to permutation or scaling of the vectors). Thanks to Lemma 12, we will then have proved the result as stated, that is for almost every M, and also almost every y 0 (since the y 0 's with at least one 0 coordinate live in a union of N hyperplanes x i = 0 where x i is i-th the coordinate in the eigenbasis, which has zero Lebesgue measure.).

Since y satisfies the linear equation dy/dt = My, we know that, for all t ≥ 0, we have y(t) = exp(tM)S 0 . Let us write λ 1 , ..., λ N the distinct eigenvalues of M, and (e n ) a corresponding eigenbasis. For all t ≥ 0, we can decompose y(t) along this eigenbasis. Let us write y 1 (t), . . . , y N (t) the corresponding coefficients so that, for all t ≥ 0, we have S(t) = n y n (t)e n . By assumption, for every n, we have y n (0) = 0. As a result, for every t ≥ 0, we have

We want to show that (y(t k )) 1≤k≤N is a basis of R N . Let k µ k y(t k ) = 0 be a linear dependence relation. Since we have

µ k e t k λn , we know that, for all n ∈ N , we have

µ k e t k λn = 0, using the unicity of coordinates in the basis (e n ) n and the fact that for all node n, we have y n (0) = 0. We use Lemma 11 to conclude.

We can now prove Lemma 6.

Proof. Let K = S + I + R be the total population irrespective of infection status (for each node n, K n = S n + I n + R n is the population of node n). By definition, for all t ≥ 0, K(t) is in the space generated by the trajectories. Then, K follows the differential equation

Using Lemma 13, we see that (K(t 1 ), ..., K(t N )) generates R N for almost every M, and K(0) = S 0 + I 0 + R 0 , so for almost every M, S 0 , I 0 , R 0 . As a result, for almost every M, X 0 = (S 0 , I 0 , R 0 ), the space generated by the trajectories contains a family which generate R N . This proves our claim.

B Proofs for 4.1, analysis of the influence of the diffusion rate

Let us first prove the following result.

Lemma 14 (Trajectories Close to a Line for Infinitely Quick Diffusion). Let us assume the initial condition X 0 is such that S 0 , I 0 and R 0 are proportional to the stationary distribution.

For every τ > 0, let us write Φ τ the solution of the system of Equation (2) where M is replaced by M/τ , that is, for every t ≥ 0, Φ τ t = (S τ (t), I τ (t), R τ (t)). Then, for all T > 0, sup

when τ → ∞, and likewise for I τ and R τ .

The results extends to cases when X 0 is not proportional to the stationary distribution, only taking the supremum over some interval [t(τ ), T ], where t(τ ) tends to 0, when τ → 0, and t(τ ) represents the time it takes for the system to converge to the stationary distribution.

the averaging. Now, since M is symmetric with respect to H, thanks to Lemma 15, we know that the trajectories it generates are also symmetric with respect to H. As a result, they belong to Fix(H). Therefore, they also satisfy the differential equation with M, as we have just proven that M and M agree on Fix(H).

Then, 2) ⇒ 1) is a direct consequence of Lemma 15.

To prove 1) ⇒ 3), let S 0 ∈ Fix(H), and let us show that MS 0 ∈ Fix(H). Choose I 0 and R 0 such that I 0 ∈ Fix(H), and define X 0 = (S 0 , I 0 , R 0 ). For every σ ∈ H, for every node i, we have

by assumption on the coefficients, and S σ(i) (0)I σ(i) (0) = S i (0)I i (0) by assumption. Moreover, we have

where the second equality is a consequence of the fact that trajectories remain in Fix(H). As a result, we have

Therefore, for every σ ∈ H, we have P (σ)MS 0 = MS 0 so that, by definition, we have MS 0 ∈ Fix(H).

Corollary 10. Let us first prove that Z = M -M vanishes on Fix(H). For all X 0 ∈ Fix(H), since M and M produce the same trajectories, we have, for all t ≥ 0, dS dt = -βS(t) I(t) + MS(t) = -βS(t) I(t) + M S(t), so that ZS(t) = (M -M ) S(t) = 0. As a result, ZS(0) = 0. This is true for all S(0) ∈ Fix(H) (as X 0 = (S 0 , I 0 , R 0 ) is arbitrary provided S 0 , I 0 and R 0 all belong to Fix(H)), so Z vanishes on Fix(H).

Corollary 10. This result is a particular case of the following Lemma 16, when we let H be the set of permutations under which the trajectories are invariant. In that case, the number of orbits of H is the number of different trajectories.

Lemma 16 (Matrices vanishing on Fix(H)). For any group of permutations H, a basis of the space of matrices Z vanishing on Fix(H) is given by the Z i,j,k 's introduced before Corollary 10, where i and j are in the same orbit under H. Consequently, the dimension of this space is (N -1)#{orbits under H}.

Proof. Let Z be such a matrix. As Z = M -M with M, M diffusion matrices, the columns of Z have vanishing sums. Moreover, Z has to vanish on any vector fixed by H. These vectors are precisely the x ∈ R N such that x i = x σ(i) for all i ∈ N , and for all σ ∈ H. Thus, they are the x's such that E k,i x = E k,σ(i) x for every i, k ∈ N , and σ ∈ H. Therefore, Z i,j,k x = 0 for all x ∈ Fix(H) when i, j, k satisfy


Pascal Robert
email: pascal.robert@grenoble-inp.fr

Nicolas Beraud

Matthieu Museau

Henri Paris

A guide to select strain sensors to be embedded in smart metal parts built by WAAM

Keywords: Smart parts, Wire and Arc Additive Manufacturing, Strain sensors, Core integration, Sensor selection guide

Resume: Additive Manufacturing gives designers freedom to explore new way of thinking smart parts. Smart parts can be defined as parts that react, adjust themselves and/or communicate. Such tasks can be performed through embedded devices or in an independent way when a change in their environment occurs, decreasing or replacing a human intervention. They are an important step toward the establishment of the 4.0 industry. This study is focused on smart part with strain monitoring. Embedding sensors enable core measurement while ensuring protection of the device from outside impairment. Additive manufacturing makes it possible to embed such kind of sensor inside a part during the manufacturing process. These newly produced parts would be able to monitor their use and optimize their lifetime. Wire Arc Additive Manufacturing is a good candidate to produce smart parts. However, this process uses gas metal arc welding technology to build metallic parts in an additiveminded way and thus imply high current and high temperatures. These consequences make it difficult for some kind of strain sensing technologies to be embedded soundly. This paper benchmark various sensor technologies for monitoring part strain and groups them into three approaches: embedding electronic devices, embedding fibre optic sensors and embedding passive exogenous indicators. Each approach is presented through one of its major sensor technology, and their pro and cons are discussed. Finally, a guide is proposed and tested on three examples, to select the strain-sensing technology adapted to restrictions imposed by the specifications of each example and by the WAAM process.

Introduction

The age of internet of thing is coming. Every day, a new connected refrigerator, a new connected car are sold. This phenomenon is enabled by the irruption of smart parts in the day-to-day life. The arrival of this kind of object in the industry is an important step in the establishment of the 4.0 industry [START_REF] Porter | How smart, connected products are transforming companies[END_REF] and is notably favoured by the freedom given by additive manufacturing.

The purpose of smart parts is to replace or decrease human intervention in the performance of a task by reacting and or communicating internal or external data when changes occur in them or their environment [START_REF] Meyer | Intelligent Products: A survey[END_REF]. The most commonly used method of creating smart parts is to embed electronic devices (sensors &/or actuator) during the manufacturing process.

Wire and Arc Additive Manufacturing (WAAM), is a process using gas metal arc welding technology in an additive manner to create metal parts [START_REF] Williams | Wire + Arc Additive Manufacturing[END_REF]. Its high deposition rate and ability to create parts without size restrictions are assets that interest industrial investors in this process.

By embedding a strain sensor in a metallic part, it can be possible to monitor its uses and enhance its lifespan. However, welding technologies develop high current and tension that are harmful to any kind of electronic devices.

This paper investigates several strain sensing technologies that can be embedded in a wire and arc additive manufacturing process to create a smart part. They can be clustered in three approach: embedding of an electronic device, embedding of an optical fibre sensor, and monitoring of an embedded indicator (often through non-destructive testing technology). A selection guide based on different criteria that can be adapted to different specifications is finally presented.

Sensor functional specifications

The research scope shall be reminded: only strain sensors are investigated; they must be integrated during the WAAM manufacturing process; the information sought must stem from core measurement of the strain in the part. Requirements of the sensor functional specifications are listed in the following table as "Main Function" or Constraints" and are intended to assist in the selection of a sensor technology.

C8

To be robust

C9

To provide additional information easily Numerous strain-monitoring technologies are reviewed in the literature. The technologies suitable to monitor core strain in a part made by WAAM process can be grouped into three approach. First, strain monitoring by embedding an electronic device such as a piezoresistive or a piezoelectric strain gauge. Second, using fibre optic sensor such as Bragg grating sensors embedded in the WAAM part. Finally, observing the strain of a passive physical indicator (a void or a sheet of a different material) embedded inside the part with a sensing technology. Many technologies used in non-destructive testing such as ultrasonic phased array detection or induction testing can be used to control the embedded indicator.

Some of the technologies studied such as laser telemetry, digital image correlation, vibrating wire sensor or coordinate measuring machine are discarded. They can monitor strain but not in a "smart-part" perspective.

Benchmark of Sensors

Embedding an electronic device

The piezoresistive sensor technology is a classic and economical way of controlling strain. It works by measuring the tension between the ends of an electronic circuit that is stick or embedded on a part which the strain has to be measured. The circuit is piezoresistive and thus, it resistance change when it is strained. Measuring its tension variation hence enable to measure its resistance variation and thus, its elongation. These sensors can be produced as thin film and be embedded in a part with not much volume hindrance. Thin film sensors are successfully embedded in a Stainless Steel Part during a Selective Laser Melting process by [START_REF] Li | Embedded sensors in layered manufacturing[END_REF] , [START_REF] Li | Embedding and characterization of fiber-optic and thin-film sensors in metallic structures[END_REF]. and [START_REF] Choi | Microfabrication and characterization of metal-embedded thin-film thermomechanical microsensors for applications in hostile manufacturing environments[END_REF].

A piezoelectric sensor produce a current when it is strained. By measuring this current, the strain of the part hosting the sensor can be calculated. (Hossain et al., 2016) embedded successfully this kind of sensor in a titanium part produced by Electron Beam Melting process.

It must be noted that these two kinds of sensors needed a thermo-electric protective case. Indeed, these sensors can be easily damaged by high current and temperature. Piezo-electric sensors can moreover lose definitely their ability to produce current under strain when they are heated above a specific temperature known as "Curie temperature" (350°C for the one used by Hossain). Furthermore, these sensors must be well bonded with their host part in order to measure strain correctly and thus, their protections must assure the bonding quality. (Figure 1) show how Li and Choi produced these cases. However, the protection thickness must be adapted to WAAM specificities (melting depth, temperature, current) that are different from the powder bed process used by Li, Choi and Hossain and may reach several millimetres. [START_REF] Li | Embedded sensors in layered manufacturing[END_REF] The electric signal used as support of the information is easy to treat by the user to compute the part strain. Furthermore, this kind of device does not require heavy instrumentation to be monitored continuously. They require to be connected to wire and sockets that must be equally protected from WAAM conditions. A basic gauge measure strain in one direction but more complex sensors shaped as "rosette" can measure strain with the same accuracy on a whole plane. Easy access to terminals going out from the part has to be assured to plug it to an acquisition unit or to a wireless sensor node as a part of a wireless sensor network [START_REF] Gupta | A Review on Wireless Sensor Networks[END_REF]. The sensor node may communicate with a transceiver using different technologies such as radio frequencies, wifi or Bluetooth [START_REF] Karray | A comprehensive survey on wireless sensor node hardware platforms[END_REF].

Embedding an optical fibre sensor

Fibre Bragg Grating (FBG) sensors is a family of optical fibre sensor sensitive to temperature and strain. Those sensors are widely investigated in literature of smart part manufacturing. They are obtained by treating locally an optical fibre by UV to create an active part (a Bragg grating). It exists different type of Bragg grating. The classic one reacts as a notch filter to light going through the fibre reflecting only the spectrum centred on a specific "Bragg" wavelength determined by the Bragg grating pitch. When the active part is stretched, its Bragg wavelength shift and this value can be used to calculate the strain of the host part, see figure 2. A more thorough explanation on the production of these sensors and on their underlying mechanisms is presented by [START_REF] Ferdinand | Capteurs à fibres optiques à réseaux de Bragg -Fabrication et caractéristiques[END_REF].

Figure 2 : Bragg Wavelength shift in the reflected spectrum of an FBG [START_REF] Luyckx | The in situ measurement of the effect of plastic deformation on the magnetic properties of steel: Part II -Permeability curves[END_REF] With fibre sensor, light convey the measure and interact hardly with the metallic host part where an electric signal would need special sheeted electrical wire to be conveyed throughout the part without being hindered by magnetic effect. Furthermore, silica is a good insolent and cannot melt because of the welding current used. Hence, optical fibre does not require any electrical insulation. However, the grating produced on the fibre can be erased if it is heated at more than 700°C. Furthermore, Silica fibre does not present a good bonding ability with metals. Hence, they need a protection for their integration during an additive manufacturing process. This protection is often bimetallic, the first one is sputtered on the fibre and is used as a seed to grow the second one as a thick nickel or copper layer by electroplating.

FBG sensors are widely used for health monitoring of big structure such as bridge, buildings, airplane wings and ships mast. Multiple Bragg Gratings can be written on a single optical fibre, it is ideal to control lengthy structures on multiple points. [START_REF] Alemohammad | Fabrication of smart cutting tools with embedded optical fiber sensors using combined laser solid freeform fabrication and moulding techniques[END_REF]) created a smart cutting tool for turning by embedding an FBG on a steel part by electroplating nickel. [START_REF] Li | Embedded sensors in layered manufacturing[END_REF], (Li and Prinz, 2003) [START_REF] Havermann | Study on fibre optic sensors embedded into metallic structures by selective laser melting[END_REF], has embedded FBG in a stainless steel part during a SLM process. They used Nickel as a protective layer to assure the bonding with steel and have a protection presenting a mechanical behaviour close from the host part material. However, Havermann explains that slippage between the fibre and the protective layer could occur during the manufacturing process and results in pre-compressing the FBG.

Grandal embed in her various studies FBG in tin during Laser and Gas Tungsten Arc welding process, [START_REF] Grandal | Technique for embedding fiber optics in metallic structures for smart material applications[END_REF], [START_REF] Grandal | Laser Cladding-Based Metallic Embedding Technique for Fiber Optic Sensors[END_REF]. She used copper and nickel protection on the fibre and advised to embed it with three welding beads parallel to the fibre in order to improve the bonding and reduce localised micro-bending. [START_REF] Feng | Temperature Sensing of Metal-Coated Fiber Bragg Grating[END_REF], [START_REF] Zhang | Sensitivity alteration of fiber Bragg grating sensors with additive micro-scale bi-material coatings[END_REF], [START_REF] Mathew | Stainless steel component with compressed fiber Bragg grating for high temperature sensing applications[END_REF], [START_REF] Kim | Characteristics of strain transfer and the reflected spectrum of a metal-coated fiber Bragg grating sensor[END_REF], [START_REF] Li | Strain Transfer Coefficient Analyses for Embedded Fiber Bragg Grating Sensors in Different Host Materials[END_REF] worked on numerical simulations of different behaviour of FBG sensor. They studied its temperature and strain sensitivity with or without different kind of protection. They study how the sheath material, thickness and bonding affect strain transfer to the fibre.

Work remain to be done in order to embed successfully an FBG during a WAAM process with high melting temperature material such as aluminium in a first step or steel in a second step.

It must be reminded that optical fibre sensors are brittle and dispendious. Moreover, spectrum analyser necessary to read the signal provided by fibres are expensive. The accuracy of these sensor is high (µε) but is diminished by the protecting sheath. Hence, these sensors are interesting to control long structure when multiples grating can be written on one long fibre and provide multiple sensing points. They need however, a protective case to be integrated to a metallic part. Numerous studies and simulations are done in order to understand how these cases impact sensitivity of the FBG to temperature and strain. Solutions exist to dissociate strain influence from temperature influence in the FBG behaviour. Continuous monitoring require that the fibre is always connected to a spectrum analyser. It is thus, suitable for part where the fibre ingress location is static and easy to protect from environment hazard.

Embedding and monitoring of a passive exogenous indicator

By embedding an exogenous indicator inside a part and monitoring it with a selected technology, it is possible to monitor locally strain in a part. The main requirement is that the indicator must be differentiated from the host part by the sensing technology used.

Various technology used for non-destructive testing or structural health monitoring can be used in this way such as X-ray technology, Thermographic imaging, Ultrasonic testing or Induction testing. Non-destructive testing technology suitable for WAAM part are reviewed by (Chauveau, 2018), (Bacelar, 2017) and [START_REF] Lopez | Mapping of non-destructive techniques for inspection of wire and arc additive manufacturing[END_REF].

Embedding an exogenous indicator give freedom to select properties of the indicator such as its size or its material. Hence, its hindrance on the part can be minimized by adjusting it shape and size and its bonding can be improved by choosing its material in respect with the sensitivity of the monitoring technology used. However, welding on a very thin sheet can make it melt easily without additional heat and current exchangers.

Various sensing technology are imaging techniques that scan the part through acoustic wave, x-ray or light spectrum emitted. They monitor the strain in the part by controlling the size of the indicator. These technologies need deep knowledge to calibrate them and expensive equipment in order to compute an image from signals sent and received. Many data are generated and sent to the acquisition unit in order to achieve the imaging and it might be challenging to control continuously strain from these technologies.

Other technologies can monitor strain in a part by controlling other indicator properties that change with strain. Potential drop, Magnetic Permeability Induction Testing (MPIT) and Eddy Current Testing are three of this kind of approach. The first one is sensitive to conductivity, the second one to magnetic permeability and the last one to both. All these technologies have inspection depth of few millimetres and are influenced by the microstructure of the host part. MPIT appears as the less limited and is presented in the next pages. [START_REF] Langevin | Variation de la perméabilité magnétique des aciers ordinaires sous l'action des efforts de traction mécanique[END_REF], observed that certain grade of steel had magnetostrictive properties. It means that under strain, their magnetic properties change and more exactly their relative magnetic permeability. Hence, it is possible to measure the strain of a part by measuring the shift occurring in its relative magnetic permeability. The relative magnetic permeability (µr) of a material define how it is interacting with a magnetic field. Some materials does not interact or hardly with a magnetic field and are called paramagnetic, (µ r ≈1, e.g. aluminium). Some attenuate it and are called diamagnetic (µr <1, e.g. copper). Finally, some enhance it greatly and are known as ferromagnetic (µr >>1, e.g. steel).

It is possible to create a magnetic field with a coil alimented by an alternative current. The coil impedance is influenced by the behaviour of the magnetic field it creates. It means that, if a ferromagnetic material is seen by the field, the coil impedance shifts according to the volume and permeability of the material put through the field. Thus, if the magnetic permeability of the material changes due to strain, this strain can be calculated from the shift in the coil impedance. The impedance of a coil is measured through its phase and its magnitude, i.e. its delay between tension and intensity, and the attenuation of the tension generator at its terminals. By plotting its Bode diagram, it can be observed that a coil has three behaviour in function of signal frequency: resistive, inductive and capacitive. By monitoring the impedance at the resonance frequency, i.e. at the inductive-capacitive regimen interface, a higher accuracy is reached in the impedance shift deduction and thus, in the permeability shift.

(Makar and Tanner, 1998), [START_REF] Gorkunov | Correlation between the stress-strain state parameters and magnetic characteristics of carbon steels[END_REF], [START_REF] Stupakov | Investigation of magnetic response to plastic deformation of low-carbon steel[END_REF], [START_REF] Swartzendruber | Effect of plastic strain on magnetic and mechanical properties of ultralow carbon sheet steel[END_REF] [START_REF] Gorkunov | Correlation between the stress-strain state parameters and magnetic characteristics of carbon steels[END_REF] (Arbenz, 2016) studied several systems for measuring the magnetic permeability of a part and developed her own system adapted to complex geometry.

To monitor a part strain of characteristic length (L) by using the permeability shift of an indicator, some requirements are needed. The part and the indicator must have very different permeability (µ). The indicator must be ferromagnetic (µd>>1), and the host part must be paramagnetic (µp≈1, aluminium or austenitic stainless steel), otherwise the host part limit the coil interrogating depth greatly. Their respective permeability behaviour under strain must be known and the insertion depth (δ) and the coil-part distance (d) must be accurate. An alternative voltage (u) with a magnitude (U) and a pulsation (𝜔) is applied to the coil. With these data, the indicator permeability influence on the coil impedance (Z) can be isolated with an LCR meter that measures its magnitude and its phase and linked to the part strain (ε), see figure 4. A sensing technology of each approach is evaluated table 2 in its ability to perform requirements deduced from the functional design presented table 1. Five performance levels (0, eliminatory; 1, poor; 3 mediocre; 7 good and 9 excellent) are assigned to each technology to compute a score. These performance levels can be assessed by a panel of experts and are based in this paper on the characteristics of sensors technologies depicted previously. No additional information 1

Next, weights are assigned to each function to give more or less importance in their fulfilment in regards to specifications of a particular subject case. These weights can be assessed by an expert group. Three different examples buildable by WAAM are selected: an airplane wing spar [START_REF] Williams | Wire + Arc Additive Manufacturing[END_REF], a lifting hook [START_REF] Jackson | Huisman's 3D printed WAAM hook lifts 80,000 kg load, 3D Printing Industry[END_REF] and a frame of a Computerized Numerical Control (CNC) Machine. Their own requirements lead to different weight distributions in functions identified, presented Table 3. can be seen that strain gauge is advised to monitor strain in a frame of a CNC machine, Fibre Bragg Grating sensor is suitable for monitoring strain in an airplane wing spar and magnetic permeability induction testing seems appropriate to monitor strain in lifting hook.

Conclusion

WAAM appear to be a suitable technology to produce smart parts. The research scope of this paper focuses on strain sensors. Different technologies and sensors were classified into three approach: embedding an electronic sensor, an optical fibre sensor or a passive exogenous strain indicator. They were analysed based on their ability to be embedded during WAAM process and their suitability for the host part requirements. Thanks to this study, the sensor selection grid is proposed and a guide is suggested to help the user to select the most suitable technology. Nevertheless, this study is still only based on the literature and future experiments need to be conducted to ensure that the advised technologies are suitable for measuring core strain in WAAM parts. Furthermore, the sensors scope can be broadened to temperature or acceleration.

Figure 1 :

 1 Figure 1 : Thin Film Strain Gauge encapsulation[START_REF] Li | Embedded sensors in layered manufacturing[END_REF]

  studied the magnetostrictive properties of various grade of steel. Gorkunov observe Figure 3, two behaviours matching with the elastic and plastic domains of a) St3, b) C45 and c) C45 quenched and tempered steel under strain.

Figure 3 :

 3 Figure 3: Magnetic permeability of three grade of steel under strain[START_REF] Gorkunov | Correlation between the stress-strain state parameters and magnetic characteristics of carbon steels[END_REF]

Figure 4 :

 4 Figure 4: Inductive behaviour of a coil in font of: a) the host part not strained b) the host part and the indicator not strained c) the host part and the indicator strained 4 Sensor Selection Grid

Table 1 :

 1

		Sensor Functional Design
	MF1	To measure core strain accurately
	C1	To be inserted during WAAM process with few requirements
	C2	To be compact
	C3	To obtain the sensing signal with few requirements
	C4	To compute strain from signal easily
	C5	To be economical
	C6	To monitor strain continuously with few requirement
	C7	To be adaptable to all kind of part (geometry and metal grade)

Table 2 :

 2

	Sensor Selection Grid
	PF-C	STRAIN	FBG	MPIT
		GAUGES
	MF1: To measure	Through thick	Through protection	Dependant to coil
	core strain	protection	9	position
	accurately	3		3
	C1: To be inserted	Requires heat	Requires heat	Requires heat and
	during WAAM	and current	protection	current evacuation
	process with few	protection	3	7
	requirements	1
	C2: To be compact	Very Invasive	Invasive long sensor	minimally invasive
		thick protection	with protection	9
		1	3
	C3: To obtain the	Sockets, DC	Ingress, Laser and	Positioning area,
	sensing signal with	generator and	Embedded and	AC generator and
	few requirements	embedded	protected non sensing	coil initial
		protected wires	optical fibre is	calibration are
		are mandatory	mandatory	mandatory
		3	3	7
	C4: To compute	Linear law with	Spectrum analyser	RLC meter needed
	strain from signal	tension	needed	7
	easily	9	1
	C5: To be	Very Economical	Expensive sensor and	Economical sensor,
	economical	sensor,	spectrum analyser;	RLC meter needed
		protection	protection needed	7
		needed	1
		9
	C6: To monitor	Need a Wifi or	Need the spectrum	Need a Wifi or
	strain continuously	Bluetooth	analyser and the laser	Bluetooth high rate
	with few	transfer device	to be connected to	transfer device and
	requirements	and a DC	the part and	an AC generator
		generator with a	acquisition unit,	with a battery
		battery linked to	ingress area need to	linked to the coil
		the part	be static	fixed on the part
		9	1	3
	C7: To be	Need access to	Best with long part,	Better with thin
	adaptable to all	sockets	need access to ingress	and paramagnetic
	kind of part	9	1	part, need
	(geometry and			positioning surface
	metal grade)			3
	C8: To be robust	Sockets and	Ingress and optical	Can be recalibrated
		electrical wire	fibre, sensible to	Can be monitored
		1	temperature	with another
			3	technique

Table 3

 3 Functions weight distributions for three examples Finally, scores can be computed in order to select the sensor or technology the most suited for each case. Table3and 4 are multiplied together as matrix to produce scores listed

	WEIGHT\MF-C	MF1 C1 C2 C3 C4 C5 C6 C7 C8 C9
	CNC Machine's Frame	9	3		1	3	3	9	9	3	9
	Lifting Hook	7	9	7	7	9	9	7	7	9	1
	Airplane Wing Spar	9	7	3	3	3	3	3	1	7	9

Table 4 .

 4 Table 4 Score obtained by each technology for three examples

	SCORE \ TECHNOLOGY	Strain gauge	FBG	MPIT
	CNC Machine's Frame	286	228	243
	Lifting Hook	358	200	397
	Airplane Wing Spar	170	232	229
	Hence, it

Acknowledgments

Authors want to acknowledge Thibault RICART, Libor RUFER and Skandar BASROUR from TIMA Laboratory for their help in sensors investigations. S.MART, Grenoble INP and the Auvergne-Rhône-Alpes region shall be thanked for their investment in WAAM technologies and research thesis.