Mathieu Aucejo
email: mathieu.aucejo@cnam.fr

Olivier De Smet

Bayesian source identification using local priors

Keywords: Inverse problem, Force identification, Group generalized Tikhonov regularization, Local priors

come

Introduction

In an industrial context, the identification of sources exciting a structure remains an important topic. It can serve for instance to control vibration at sources, identify excitation sources in unmeasured areas, estimate reaction forces at boundaries or establish excitation models to numerically predict the response of structures. However, the identification of mechanical sources from vibration measurements is known to be an ill-posed inverse problem, meaning that the existence of a unique stable solution is not guaranteed. A classical approach to bypass this difficulty consists in including in the formulation of the inverse problem some prior information on the noise and the spatial distribution of sources to identify to constrain the space of solutions.

For this purpose, the Bayesian formalism is generally adopted, since it offers a rigorous mathematical framework that allows combining both probabilistic and mechanical data. Among all the methods used for mechanical source identification and deriving from the Bayesian statistics , the Tikhonov regularization is certainly the most popular [START_REF] Busby | Optimal regularization of an inverse dynamics problem[END_REF][START_REF] Thite | The quantification of structure-borne transmission paths by inverse methods. Part 2 : Use of regularization techniques[END_REF][START_REF] Liu | Reducing the impact of measurement errors when reconstructing dynamic forces[END_REF][START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF]. Unfortunately, this regularization leads to a systematic smoothing of regularized solutions [START_REF] Renzi | Identification Expérimentale de Sources vibratoires par Résolution du probleme Inverse modelisé par un opérateur Eléments Finis local (Vibratory sources identification by solving the inverse problem modeled by a local Finite Element operator[END_REF], which is not a desirable effect when a localized source has to be identified.

To remedy this problem, several approaches also based on the Bayesian inference have been proposed. In particular, it is useful to cite the work of Renzi et al. [START_REF] Renzi | Identification Expérimentale de Sources vibratoires par Résolution du probleme Inverse modelisé par un opérateur Eléments Finis local (Vibratory sources identification by solving the inverse problem modeled by a local Finite Element operator[END_REF][START_REF] Renzi | Identification of vibration excitation using a regularized finite element operator and a deconvolution post-process[END_REF] in which they apply a Richardson-Lucy deconvolution.

From the Bayesian standpoint, this approach consists in assuming a Poisson noise model and a uniform distribution for the source vector to identify [START_REF] Molina | Prior models and the richardsonlucy restoration method[END_REF]. Practically, this leads to impose a positivity constraint on the excitation field to identify, which implies the loss of the phase relationships between the identified sources. Another technique, based on the use of generalized Gaussian laws to reflect prior information on the noise and the sources to identify, has been proposed in [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF]. This approach has the advantage to alleviate the positivity constraint and to be flexible enough to identify sources of various types. However, in the procedures described above, the a priori on the spatial distribution of the sources remains global. Incidentally, poor identification can be obtained, if the actual excitation field combines both localized and distributed sources, since the a priori has to reflect a compromise between two contradictory distributions. Finally, it is worth citing the work of Zhang et al. [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF] in which they propose to solve the Bayesian force reconstruction problem with an uncertain model from Monte Carlo Markov Chain (MCMC) methods. Such an approach has the merit of providing a credible interval on the reconstructed force, which allows assessing the robustness of the identified solution.

The present paper aims at showing that one's prior knowledge of the system can be incorporated to aid the identification process. More precisely, the main objective is to provide a flexible identification methodology, based on Bayesian statistics, able to take advantage of prior local information available on both the nature and the location of excitation sources. To this end, the structure is divided into several identification regions, in which it is assumed that local priors on the nature of the sources follow generalized Gaussian distributions. In doing so, prior information can be accurately adapted to a prescribed zone of the studied structure and, therefore, one can fully exploit the spatial information that is generally available in practical situations.

From a mathematical point of view, the solution of the problem is defined as the maximum a posteriori estimate of the posterior distribution. Practically, one seeks the solution of the dual minimization problem, which is solved in an iterative manner using a Generalized Iteratively Reweighted Least-Squares (GIRLS) algorithm [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Rodriguez | An iteratively weighted norm algorithm for total variation regularization[END_REF].

To highlight the main features of the proposed methodology, the present article is divided into four parts. In section 2, the formulation of the proposed source identification problem is introduced. In section 3, the resolution of the resulting source identification problem is presented. As mentioned previously, the resolution is performed using the GIRLS algorithm described in [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF]. In order to make this paper clearer and more self-contained, the GIRLS algorithm is detailed by introducing further insights regarding the choice of the tuning parameters involved in its definition. Finally, the validity of the proposed methodology is illustrated both numerically and experimentally in sections 4 and 5. In particular, the proposed validations reveal that using local instead of global priors leads to more accurate identifications when a structure is excited by sources having different spatial distributions.

Formulation of the Bayesian source identification problem

This section aims at presenting the derivation of the formulation of the identification problem and introducing the main tuning parameters involved in the proposed formulation.

Problem description

Let us consider the situation where a vibration field X is caused by an unknown excitation field F. If the structure is linear, its dynamic behavior is completely defined by the transfer functions matrix H, relating the vibration field X to the excitation field F. Practically, the vibration field X, measured over the surface of a structure, is biased by two independent contributions. The first one is related to measurement errors, since the transducers used to perform measurements are imperfect. This contribution is generally modeled by an additive noise. The second one is related to modeling errors, which reflect that the model used to determine the transfer functions matrix H is imperfect. If modeling errors are supposed small enough, it is reasonable to incorporate this bias in the additive noise [START_REF] Zhang | Bayesian force reconstruction with an uncertain model[END_REF].

As a result, the measured vibration field X is obtained from the following direct formulation:

X = HF + N, (1)
where the transfer functions matrix H, the excitation field F and the global noise N are fully determined.

Unlike to the direct problem, the structural source identification problem consists in estimating the unknown excitation field F acting on a structure from the vibration field X measured on its surface only. In other words, the noise N is another unknown of the inverse problem. To deal with such a problem efficiently, the Bayesian paradigm is well adapted and is adopted in this work.

The Bayesian paradigm consists in considering all the parameters of the problem as random variables. Consequently, it is an uncertainty model, in which the uncertainty on each parameter is modeled by a probability distribution, describing the state of knowledge or the a priori on this parameter.

From the mathematical standpoint, the formulation of the Bayesian source identification problem relies on the Bayes' rule:

p(F| X) ∝ p( X|F)p(F), (2)
where:

- From the Bayes' rule, one seeks the most probable excitation field F given the measured vibration field X. Such a solution is a point estimate of F corresponding to a mode of the posterior probability distribution. Consequently, the solution of the identification problem is sought as the maximum a posteriori estimate, that is:

F = argmax F p(F| X) = argmax F p( X|F)p(F). (3)
Practically, it is generally easier to find the solution of the following dual minimization problem:

F = argmin F -ln p(F| X) = argmin F -ln p( X|F) -ln p(F). (4)
From Eqs. ( 2)-( 4), it is clear that the more prior information reflects the nature of the noise and the sources to identify, the more the posterior probability is high. From this observation, it is expected that proper results can be obtained if the likelihood function and the prior probability distribution are properly chosen, since the identification process gives the most probable solution given prior information on the noise and the nature of the sources to identify.

Choice of the likelihood function and prior probability distribution

As suggested in the previous section, the quality of the source identification is closely related to the choice of the likelihood function and prior probability distribution. That is why, this section aims at explaining the choices made in this paper to have a closed-form expression of the identification problem, that allows properly reflecting prior information on the noise and the sources to identify.

Choice of the likelihood function

As mentioned in section 2.1, the likelihood function reflects the uncertainty related to the measurement of the vibration field X. By definition [see Eq. ( 1)], this uncertainty is related to the noise N. Consequently, the likelihood function can be rewritten under the following form:

p( X|F) = p( X -HF|N), (5)
meaning that the likelihood function can be seen as the probability of obtaining X -HF = 0 given the noise N. In other words, the likelihood function measures the fidelity of the vibration model, corresponding to the noiseless vibration field X = HF, to the measured data. As a result, the likelihood function is also written:

p( X|F) = p(N), (6)
where p(N) is the prior probability distribution of the noise, representing the a priori of the experimenter on the nature of the noise inducing a discrepancy between the measured vibration field and its noiseless estimate before any measurement of the vibration field.

If the noise is supposed spatially white (meaning that the components of N are independent and identically distributed random variables) and due to multiple independent causes, then the prior probability distribution of the noise can be represented by a multivariate normal distribution with zero mean and variance α 2 :

p(N) ∝ exp - 1 2 N 2 2 α 2 , (7)
where • 2 is the L 2 -norm, also known as the Euclidean norm.

It can be inferred from Eqs. (1), ( 6) and ( 7) that the likelihood function is written:

p( X|F) ∝ exp - 1 2 X -HF 2 2 α 2 . ( 8
)
Finally, it is worth mentioning that in practical situations the noise may not be strictly white or even Gaussian. However, assuming a spatially Gaussian white noise allows limiting the number of hyperparameters of the Bayesian model. Consequently, this hypothesis can seem coarse or indecorous, but it is actually relevant as highlighted by the experimental results presented in the next of the paper.

Choice of the prior probability distribution

The prior probability distribution reflects the uncertainty related to the unknown excitation field F. Consequently, it can be seen as the a priori of the experimenter on the sources to identify. To take advantage of local information available a priori on the nature and the spatial distribution of sources, it is supposed that the structure is excited in N different regions by uncorrelated excitations of various types (localized or distributed). This allows considering the local excitation fields F i as independent random vectors and thus writing the prior probability distribution as the product of local prior probability distributions p(F i ), that is:

p(F) = N i=1 p(F i ), (9)
where p(F i ) reflects the prior knowledge of the nature of the sources in the zone i.

Furthermore, it is assumed that the components of each local excitation vector F i are independent and identically distributed random variables (i.e.

uncorrelated but having the same probability distribution) following a generalized Gaussian distribution with zero mean. As a consequence, each local excitation field follows a multivariate generalized Gaussian distribution with zero mean. From the mathematical standpoint, the local prior probability distributions are thus written:

p(F i ) ∝ exp - 1 q i L i F i q i q i β q i i , (10) where:
q i is the shape parameter of the distribution in the zone i. Its value is defined in the interval ]0, +∞[;

-• q i is the L q i -norm or quasi-norm, if q i ≥ 1 and q i < 1 respectively;
β i is the scale parameter of the distribution, which defines the dispersion of the distribution around the mean. It is therefore a generalized measure of the variance of the distribution;

-L i is a differentiation matrix controlling the regularity of the solution in the zone i.

It should be noted that the choice of a multivariate generalized Gaussian distribution allows a great flexibility for describing prior knowledge of the sources to identify as it will be made clearer in section 3.2.1. As for the assumption of a Gaussian white noise, the hypothesis of independence used to derive the prior probability distribution can be discussed. Actually, this hypothesis has been introduced here because the correlation between each component of the force vector is difficult to assess a priori. However, this assumption is typically violated when one tries to identify excitation forces and forces at boundaries. In this particular situation, forces at boundaries are functionally dependent on the excitation forces. From a theoretical standpoint, Tikhonov-like approaches are known in Machine Learning and Data Mining as naive Bayes classifiers. In these fields, it has been shown that naive Bayes classifiers works well not only when the parameters involved in the formulation are completely independent, but also when the parame-ters are functionally dependent [START_REF] Rish | An empirical study of the naive bayes classifier[END_REF]. For all the aforementioned reasons, the hypothesis of independence is relevant.

Practical form of the identification problem

In order to derive the practical form of the identification problem, Eqs. ( 8) and ( 10) are introduced into Eq. ( 4). In doing so, one obtains the general formulation of the identification problem given below:

F = argmin F 1 2 X -HF 2 2 α 2 + N i=1 1 q i L i F i q i q i β q i i . (11)
At this stage, 2N +1 parameters, namely (q i , β i , α), are necessary to mathematically describe prior information on the noise and the sources to identify.

However, it can be emphasized that the parameters q i and β i play a similar role in this formulation. Therefore, only one set of parameters, q i or β i , can be used. In the present paper, it has been chosen to work with the shape parameter q i , while setting the values of the parameters β q i i to a unique constant value β. Taking this observation into account, one finally gets the final form of the identification problem, which is a group generalized Tikhonov regularization :

F = argmin F 1 2 X -HF 2 2 + λ N i=1 1 q i L i F i q i q i , (12)
where λ = α 2 /β is the so-called regularization parameter, that controls the trade-off between prior information on the noise and prior information on the solution.

Before going any further, it could be interesting to discuss the frequency limitations of the proposed identification procedure. Actually, for a fully determined problem, the limitations are not related to the formulation of the inverse problem itself, but to the measurement or the calculation of the transfer functions matrix H and the measurement of the vibration field X.

Regarding the calculation of H, the frequency limitation is related to the model of the structure. If the model is purely analytic, then the limit is related to the assumptions made to establish the mechanical model. If the model is obtained by FEM, then the limit is fixed by the assumptions made to derive the mechanical model and the mesh used to discretize the structural domain (generally built using 6 nodes per wavelength). Regarding the measurement of the vibration field, the frequency limitation is mainly related to the bandwidths and the resolution of the transducers.

It can also be noticed that the proposed formulation includes the standard Tikhonov regularization, since the latter is recovered if q i = 2 for all i ∈ [1, N ] in Eq. [START_REF] Bjorck | Numerical methods for Least-Squares problems[END_REF]. Furthermore, it is emphasized that the selected zones can cover only a part of the whole structure, implying that no a priori on the spatial distribution of sources is assumed in non-selected zones. In these particular regions, the identification is therefore only based on the fitting to the experimental data, i.e. on prior information on the noise. Finally, it is worth noting that the robustness of the proposed method against modeling and measurement uncertainties can not be assessed, since the proposed formulation provides a point estimate of the force vector to identify. Nevertheless, the proposed methodology allows obtaining consistent identifications as it will be shown in sections 4 and 5.

Resolution of the source identification problem

Although the proposed formulation is highly flexible to describe prior knowledge of the sources to identify, its solution has generally no closed-form expression. That is why, the identification problem is solved using an iterative algorithm, which is a generalized version of the Iteratively Reweighted Least-Squares algorithm [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Rodriguez | An iteratively weighted norm algorithm for total variation regularization[END_REF][START_REF] Bjorck | Numerical methods for Least-Squares problems[END_REF]. The proposed algorithm has the advantage to be easy to implement and general enough to cope with a wide range of configurations depending on the values of q i . In the next of the paper, it is referred to as Generalized Iteratively Reweighted Least-Squares (GIRLS) algorithm.

General principle

The core idea of the GIRLS algorithm consists in replacing the direct resolution of the minimization problem [see Eq. ( 12)] by an equivalent iterative process having an explicit solution at each iteration. For this purpose, one has to notice that the L q -norm to the power of q can be rewritten as a squared weighted L 2 -norm as it is shown in Eq. ( 13):

∀x n , ∀q, 1 q x q q = 1 q n |x n | q = 1 2 n w(x n ) |x n | 2 = 1 2 W 1/2 x 2 2 , (13)
where w(x n ) = 2 q |x n | q-2 is the weighting coefficient and W = diag[w(x 1 ), . . . , w(x n )] is the corresponding definite positive weighting matrix.

Since the weighting matrix depends on x, the equality given in Eq. ( 13) cannot be verified directly, which justifies the use of an iterative algorithm to obtain the solution of the identification problem. As part of an iterative process, the aim is thus to find the solution x (k+1) at iteration k+1 from the solution x (k) at iteration k by setting w(x n ) = w(x (k) n ), so as to recover the identity given by Eq. ( 13) when the iterative process has converged.

When applying this idea to the proposed formulation, it can be inferred that the excitation field

F (k+1)
at iteration k+1 is solution of the following minimization problem:

F (k+1) = argmin F 1 2 X -HF 2 2 + λ (k+1) 2 N i=1 W (k) i 1/2 L i F i 2 2 , (14)
where λ (k+1) is the regularization parameter at iteration k+1.

For numerical purposes, the definition of the weighting matrix W (k) i is slightly modified compared to that given in Eq. ( 13), since it is defined such that:

W (k) i = diag 2 q i T L i F (k) i , (15)
with:

T x (k) j =      x (k) j q i -2 if x (k) j > (k) i (k) i q i -2 if x (k) j ≤ (k) i , (16) where x
(k) j is the j-th component of the vector x (k) and

(k) i is a small real positive number acting as a damping parameter. It allows avoiding infinite weights when x (k) j → 0 and q i < 2.

For implementation purposes, it should be emphasized that Eq. ( 14) can be more concisely expressed by introducing two matrices, W (k) and L, defined as follows:

W (k) =      W (k) 1 0 . . . 0 W (k) N      (17)
and

L =      L 1 0 . . . 0 L N      , (18)
if the selected zones entirely cover the studied structure. In all other cases,

L =      L 1 0 0 . . . . . . 0 L N 0      . (19)
Introducing W (k) and L into Eq.( 14), it readily comes that:

F (k+1) = argmin F 1 2 X -HF 2 2 + λ (k+1) 2 W (k) 1/2 LF 2 2 . ( 20
)

Choice of the tuning parameters

According to Eqs. ( 14), ( 15) and ( 16), 2N + 1 tuning parameters are involved in the resolution of the identification problem, namely the shape parameters of the distribution q i , the regularization parameter λ (k+1) and the damping parameters

(k)
i . In the following subsections, the selection of each parameter is detailed.

Choice of the shape parameters

To properly choose the value of the shape parameter in each zone, one has to keep in mind that the shape of the distribution or equivalently the shape of the penalty function, i.e. x q i q i , has a major influence on the solution of the identification problem [see Fig. 1]. Indeed, as shown in Fig. 1b, it is clear that the smaller q i is, the larger is the weight on small values of x. On the contrary, for large values of x, the smaller q i is, the smaller is the weight on these values. This difference in the weighting behavior for small and large values of x can be directly observed in the solutions obtained, since for q i ≤ 1 the solution will tend to have only a few non-zero values, while for q i = 2 the solution will tend to have only a few very small values. Figure 1: Selection of the shape parameter q i -(a) Probability distribution exp -

x q i q i qiβ q i i (β qi i = 1
) and (b) Corresponding penalty function x qi qi for (-) q i = 2, (--) q i = 1 and (-• -) q i = 0.5

From this observation, it can be inferred that q i can be set to 2 if the solution sought is supposed to be distributed [START_REF] Boyd | Convex optimization[END_REF], while q i can be chosen equal to or less than 1 if the solution is supposed to be localized [START_REF] Chartrand | Nonconvex regularization for image segmentation[END_REF][START_REF] Nikolova | Efficient reconstruction of piecewise constant images using nonsmooth nonconvex minimization[END_REF]. This rule of thumb essentially holds if the differentiation matrices L i correspond to the identity matrix. Indeed, if ones tries to identify a piecewise continuous excitation field in a zone i, one has to promote a solution for which the first derivative contains only a few non-zero values. In this situation, the differentiation matrix differs from the identity matrix. As an example, for a one dimensional structure, it is then recommended to choose q i ≤ 1 and use a differentiation matrix of the form:

L i =         -1 1 -1 1 . . . -1 1         . ( 21
)

Choice of the regularization parameter

The value of the regularization parameter can be chosen automatically from various selection methods, such as the Morozov's discrepancy principle [START_REF] Scherzer | The use of Morozov's discrepancy principle for Tikhonov regularization for solving nonlinear ill-posed problems[END_REF], the Generalized Cross Validation [START_REF] Golub | Generalized cross-validation as a method for choosing a good ridge parameter[END_REF], the Bayesian estimator [START_REF] Antoni | A Bayesian approach to sound source reconstruction: Optimal basis, regularization, and focusing[END_REF] or the L-curve principle [START_REF] Hansen | Rank-Deficient and Discrete Ill-Posed Problems: Numerical Aspects of Linear Inversion[END_REF]. In the present paper, the regularization parameter λ (k) is obtained using the L-curve principle at each iteration.

For the sake of completeness, it can be reminded that the L-curve principle consists in plotting, in a log-log scale, the L 2 -norm of the regularization term (i.e.

W (k) 1/2
LF 2 ) versus the L 2 -norm of the data fidelity term (i.e.

X-HF 2 ) for various values of the regularization parameter. Generally, the resulting curve exhibits an L-shape. The optimal value of the regularization parameter is the value for which the curvature of the L-curve is maximal, i.e.

corresponding to the corner of the L-curve. To illustrate the selection of the regularization parameter from the L-curve principle, an example extracted from the experimental validation presented in section 5 is proposed in Fig. 2. However, the computation of the regularization parameter at each iteration is intensive. Accordingly, it is interesting to determine whether a unique value can be used throughout the iterative process. As pointed out in [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF],

a series of numerical experiments has shown that using λ (k) = λ (0) allows obtaining similar results to those obtained from updated λ (k) within a comparable number of iterations. The logical outcome of this numerical study is that the choice of the regularization parameter can be performed once for all during the initialization step of the iterative process.

Choice of the damping parameters

The damping parameter 3]. From our own numerical tests, it has been found that using p = 5% leads to satisfying results.

|L i F (k) i |
As for the regularization parameter, accurate solutions can be obtained using a fixed value of the damping parameter. Indeed, from a set of numerical experiments, it has been noticed that using the value computed during the initialization step throughout the iterative process (i.e.

(k) i = (0) ) has no significant influence on both the quality and the time-performance of the algorithm.

Choice of the initial guess and stopping criterion

Since the proposed resolution algorithm is iterative, another crucial issue is the choice of the initial guess as well as that of the stopping criterion.

Choice of the initial guess

Choosing a good initial guess is one of the keys for a successful identification, since the proposed formulation is non-convex when q i < 1. In such a situation, the existence of a unique minimizer is not guaranteed. Consequently, the identified solution F strongly depends on the initial guess F

.

A good initial guess can be defined as a coarse solution of the problem, easy to compute, but sufficiently close to the solution sought to ensure the convergence of the iterative process. A possible choice that fulfills this requirement is the solution of the standard Tikhonov regularization, that is:

F (0) = argmin F 1 2 X -HF 2 2 + λ (0) 2 F 2 2 . (22)

Choice of the stopping criterion

As an iterative process, the GIRLS algorithm requires the definition of a practical test that allows determining when to stop the algorithm. In the present paper, the stopping criterion is related to the relative variation of the functional J( F

(k) ) = 1 2 X -H F (k) 2 2 + λ (k) 2 W (k-1) 1/2 L F (k) 2
2 between two successive iterations. As classically done in the literature, the relative variation δ of the functional J is defined such that:

δ = J F (k+1) -J F (k) J F (k) . ( 23
)
From there, the iterative process is stopped when the relative variation δ is less than or equal to some tolerance. Experimentally, it has been found that setting the tolerance to 10 -8 allows obtaining consistent identified solutions, while preserving the time-performance of the GIRLS algorithm.

Summary

To clearly highlight each step of the resulting iterative process, a comprehensive overview of the proposed version of the GIRLS algorithm is given in table 1. from Eq.( 22)

Compute (0) from the cumulative histogram of | F (0) | Initialize δ to 1 Iteration: while δ > tol Compute W (k) i
from Eq. ( 15) using

(k) i = (0)
Construct W (k) from Eq. ( 17)

Construct L from Eq. ( 18) or Eq.( 19)

Compute

F (k+1)
from Eq. ( 14) using λ (k+1) = λ (0)

Update δ using Eq. ( 23)

k ← k + 1 end return F ← F (k)

Numerical validation

In general, forces and moments can be applied on a structure or exists at boundaries. However, the identification of moments requires a careful attention [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF], since rotations are difficult to measure directly. For this reason, the proposed numerical validation is restricted to the identification of forces only, in order to focus on the particular interest of including local prior information in the identification process, which is the originality of the paper.

Description of the test case

In the present numerical validation, one seeks to identify a mechanical point force acting on a thin simply supported steel plate. In such a configuration, the excitation field exhibits two types of spatial distribution over the structure, namely a smooth distribution of the reaction forces at boundaries and a singular distribution around the location of the point force. The main simulation parameters of the proposed test case are listed in table 2.

The exact vibration displacement field X exact is computed from a FE mesh of the plate made up of 567 linear quadrilateral shell elements, assuming that only bending motions are measurable. In practical situations, the exact vibration field is corrupted by noise. This phenomenon is modeled by adding to the exact displacement field an additive Gaussian white noise as suggested by Eq. ( 1), that is:

X = X exact + N, (24)
where the noise vector N is such that [START_REF] Leclere | Acoustic imaging using under-determined inverseapproaches: Frequency limitations and optimal regularization[END_REF]: Furthermore, as suggested in [START_REF] Renzi | Vibratory source identification by using the finite element model of a subdomain of a flexural beam[END_REF], a FE model of the plate with free boundary conditions is used to compute the transfer functions matrix H, which is dynamically condensed over the measurable dofs only [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF][START_REF] Ibrahim | Structural force identification at unmeasured locations[END_REF]. The main interest in using free boundary conditions to model dynamic behavior of the plate is to allow the identification of external excitations acting on the structure as well as reaction forces at boundaries. Finally, one can notice that the proposed FE model of the plate is theoretically valid up to 650 Hz, corresponding to the frequency for which the structure is discretized with 6 nodes per bending wavelength.

N = β a e jφa , (25)

Application

The reference force vector F ref is obtained from the transfer functions matrix H and the exact displacement field X exact thanks to the following equation:

F ref = H -1 X exact . (26)
Fig. 4 presents the reference force vector to identify at 450 Hz. As expected, a unit point force appears at (x 0 ,y 0 ) = (1.04 m,0.93 m) as well as smooth reaction forces at boundaries of the plate. 26), one obtains the identified excitation vector F presented in Fig. 5.

Obviously, the source identification fails, since the reconstructed excitation field is highly dominated by the noise. This result is easily explained by the ill-conditioning of the transfer function matrix H. In the present case, the condition number of H is equal to 2×10 4 . This indicates the presence of small singular values, which are responsible for the amplification of the terms related to the noise and consequently cause the identification to fail [START_REF] Hansen | Discrete Inverse Problems: Insight and Algorithms[END_REF]. As done in Ref. [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF], the first idea to stabilize the identification problem consists in using a global a priori on the spatial distribution of sources. In such a configuration, it is clear that a compromise has to be found between the smoothness of the reaction forces at boundaries and the singularity of the point force. Considering explanations given in section 3.2.1, a relevant compromise is obtained after 31 iterations for q = 1.1 and L = I as pre-sented in Fig. 6. With this value of the shape parameter, the identified force amplitude F 0 is equal to 0.72 N instead of 1 N. However, the a priori on the spatial distribution of sources being perfectly adapted neither to localized sources nor to distributed sources, the identified sources do not perfectly match with the actual ones. An even better result can be obtained by using properly our prior knowledge of the location and the nature of the sources acting on the structure. Indeed, as presented in Fig. 4, the force vector to identify is sparse except in the vicinity of the boundaries of the plate. This observation leads to divide the structure into two identification regions: (i) a central region associated to the shape parameter q 1 and containing the point force only, in which a sparsity-promoting prior is required, and (ii) a region associated to the shape parameter q 2 and corre-sponding to the boundaries of the plate, in which a prior promoting smooth solutions has to be employed [see Fig. 7]. Following the considerations described in section 3.2.1, one sets (q 1 , q 2 ) = (0.5, 2). When choosing these values, one obtains the identified force vector presented in Fig. 8. It is noteworthy that the identified force vector is very similar to the reference one, meaning that chosen values are well adapted to describe the spatial distribution of sources in each selected zones. From a numerical standpoint, the immediate consequence of the accordance of the local priors with the actual local spatial distribution of sources is the decrease in the number of iteration of the GIRLS algorithm [START_REF] Aucejo | Structural source identification using a generalized Tikhonov regularization[END_REF]. Indeed, only 11 it-erations have been required to reach the tolerance and the identified force amplitude F 0 is equal to 0.995 N instead of 1 N. From this particular point of view, it can be noticed from Eq. ( 2) that the more prior information on the noise and the nature of the sources is meaningful, the more the confidence in the reconstruction is high. Numerically, this implies that the convergence of the GIRLS algorithm to a reliable solution is faster when choosing proper values of the shape parameters q i . The analysis presented above has been performed for a frequency that do not correspond to a natural frequency of the plate. However, it should be noted that source identification at natural frequencies of a lightly damped structure is generally far more difficult, since at these frequencies the shape of the vibration response of the structure is close to that of the corresponding mode shapes. Consequently, the nearfield information is somewhat masked.

From the standpoint of the identification process, this can be regarded as an increase of the noise level. To assess the ability of the proposed approach to perform consistent identification at the natural frequencies of the structure, the identification procedure is applied on the test case described previously, but for resonances of the plate, namely 85 Hz, 278 Hz and 572 Hz. The corresponding identifications are presented in Fig. 9. Obtained results show that for each natural frequency the shape of the identified force vectors are consistent with the reference ones. However, one has to noticed that the amplitude of the identified point force F 0 is generally underestimated. Indeed, the identified amplitude is around 0.83 N instead of 1 N for the three cases. Nevertheless, these results show that the proposed method allows having credible information on the sources acting on the structure at resonance frequencies.

Experimental validation

The structure of interest is a steel parallelepiped box, excited on one of its faces by a shaker supplied by a white noise and equipped with a force sensor, as presented in Fig. 10. The main features of the proposed experimental validation are listed in table 3.

In this experiment, the source identification is performed on the excited As for the numerical validation, let us start by identifying the force vector F using a global prior on the nature of the sources to identify. In the present case, it can be noticed that the analysis of the experimental set-up suggests that the spatial distribution of sources acting on the studied part of the structure is similar to that of the numerical test case. Consequently, it is natural to choose q = 1.1 and L = I. After 85 iterations, one obtains the identified force vector presented in Fig. 11. It can be noticed that the spatial distribution of the point force is spread around its actual location, indicating that the global a priori is not perfectly adapted to identify localized sources.

On the contrary, the spatial distribution of reaction forces is consistent with what is expected when analyzing the experimental set-up. It should however be noted that the identified forces at boundaries are the sum of pure reaction forces and equivalent forces due to the potential existence of boundary moments, because a force-only model of the structure is adopted here. As for the point force, one can, nevertheless, suspect that the identified amplitudes are underestimated, since the global a priori tries to reflect the compromise between the singularity of the point force and the smoothness of the reaction forces at boundaries. As observed in the numerical validation, a better identification can be obtained by using properly our a priori on the nature and the location of sources acting on the structure. From a careful analysis of the experimental set-up, one can infer that the target force vector is very sparse except in the vicinity of the boundaries of the studied plate, where the spatial distribution of the reaction force is likely to be smooth. Accordingly, two identification regions can be distinguished as illustrated in Fig. 12. The selected regions have the same characteristics as those defined in the numerical validation. This observation allows assuming that (q 1 , q 2 ) = (0.5, 2). Similarly, the reaction forces at boundaries remain consistent with our expectations. Finally, it can be noticed in passing that only 26 iterations have been necessary to reach the tolerance, confirming that properly exploiting the spatial information on the nature and the location of sources allows improving both the quality of the identification and the performance of the As for the numerical validation, the previous analysis has been performed outside natural frequencies of the structure. Nevertheless, this analysis is not sufficient to validate the proposed methodology in an experimental context, since identifying sources outside resonance frequencies is not the more challenging case. To have better insights on the identification ability of the approach, identifications are performed at resonance frequencies of the structure, namely 86.25 Hz, 212.5 Hz and 453.75 Hz. The results presented in Fig. 14 clearly show that the corresponding identifications are qualitatively consistent with our expectations. Furthermore, it should be noted that the location of the point force is still well identified, while its normalized ampli-tude is underestimated since it is around 0.75 for the three cases presented here.

Finally, to give a comprehensive overview of the proposed identification strategy, the force spectra identified at the point force location is given in Figure 15. In an ideal situation, the identified force amplitude should be equal to 1 at each studied frequency. It is obviously not the case in the present validation. However, one can notice that the identified force spectra is reasonably close to the target value, as the discrepancy varies from -5 dB to 2 dB over the entire frequency range and the frequency-averaged force amplitude is approximately equal to 1.02. Consequently, all these results are encouraging since credible results are obtained even with a model that does not strictly correspond to the actual structure.

Conclusion

In the present paper, a general methodology has been presented for identifying mechanical sources acting on a structure. The initial motivation of this study was to develop a method able to fully exploit the spatial information on the nature and the location of sources. Indeed, in the vast majority of the methods available in the literature, such an information remains unused, whereas it is generally roughly known in practical situations. To remedy this problem, we have used the Bayesian paradigm to derive a group generalized Tikhonov regularization, that allows introducing local a priori on the spatial distribution of sources. In doing so, the identification procedure is expected to be more flexible and accurate, since the a priori can be precisely adapted

  p(F| X) is the posterior probability distribution, representing the probability of observing F given the measured vibration field X. In other words, it represents what is known on the excitation field after making vibration measurements; -p( X|F) is the likelihood function, representing the probability of measuring X given an excitation field F. It reflects the uncertainty related to the measurement of the vibration field, which is itself related to the noise; -p(F) is the prior probability distribution, representing the prior knowledge of the excitation field before measuring the vibration field.

Figure 2 :

 2 Figure 2: Illustration of the L-curve principle -(a) Plot of the L-curve and (b) Curvature of the L-curve versus regularization parameter λ (k) -(•) Location of the optimal regularization parameter

  automatically for each zone from the cumulative histogram of |L i F (k) i | [20], meaning that the value of the damping parameter is chosen so that a fixed percentage p of the values of 18 |L i F (k) i | are less than or equal to (k) i [see Fig.

Figure 3 :

 3 Figure 3: Illustration of the choice of the damping parameter -Cumulative histogram of

E = 2 . 1 ×

 21 10 11 Pa Density ρ = 7800 kg.m -3 Structural damping η = 0.01 Amplitude of the force F 0 = 1 N Location of the force (x 0 , y 0 ) = (1.04 m, 0.93 m) Studied frequency f 0 = 450 Hz where β a = p n M -1/2 X exact 2 Z and φ a = U, with p n the noise percentage, M the number of measurable dofs, Z a Gaussian random vector and U a random vector uniformly distributed over the interval [0, 2π]. In what follows, the noise percentage p n is set to 2%.

Figure 4 :

 4 Figure 4: Reference force vector F ref at 450 Hz

Figure 5 :

 5 Figure 5: Naive identification of the force vector F at 450 Hz

Figure 6 :

 6 Figure 6: Identification of the force vector F at 450 Hz from corrupted data when using a global prior on the spatial distribution of sourcesq = 1.1 and L = I

Figure 7 :

 7 Figure 7: Definition of the selected zones -(•) zone 1 (sparse), (×) zone 2 (smooth) and (♦) location of the point force

Figure 8 :

 8 Figure8: Identification of the force vector F at 450 Hz from corrupted data when using local priors on the spatial distribution of sources -(q 1 ,q 2 ) = (0.5,2) and L = I

Figure 9 :

 9 Figure 9: Identification of the force vector F at resonance frequencies of the structure -(a) Reference at 85 Hz, (b) Identification at 85 Hz, (c) Reference at 278 Hz, (d) Identification at 278 Hz, (e) Reference at 572 Hz and (f) Identification at 572 Hz -(q 1 ,q 2 ) = (0.5,2) and L = I

Figure 10 :

 10 Figure 10: Experimental set-up

Figure 11 :

 11 Figure 11: Experimental identification of the force vector F c at 600 Hz using a global prior on the spatial distribution of sourcesq = 1.1 and L = I

Figure 12 :

 12 Figure 12: Definition of the selected zones -(•) zone 1 (sparse), (×) zone 2 (smooth) and (♦) location of the point force

Figure 13 :

 13 Figure 13: Experimental identification of the force vector F c at 600 Hz using local priors on the spatial distribution of sources -(q 1 , q 2 ) = (0.5, 2) and L = I

Table 1 :

 1 Generic GIRLS algorithm

	Proposed GIRLS Algorithm

Inputs: Transfer functions matrix: H, Measured vibration field X, Selected zones {i}, Differentiation matrix for each selected zones: L i , Tolerance: tol Output: Estimated force vector: F Initialization: Compute λ (0) from the L-curve principle Compute F (0)

Two extensions of the identification strategy developed in this paper could be considered. First, the application of the proposed approach to non-linear structures by linearizing the problem around a certain state of the system and apply the GIRLS algorithm. Second, the implementation of a full Bayesian approach based on MCMC techniques for which the point estimates obtained with the proposed methodology could be considered as a valuable starting point.


Térence Bayen
email: terence.bayen@umontpellier.fr

Alain Rapaport
email: supagro.alain.rapaport@inra.fr

Fatima Zahra

Fatima-Zahra Tani

Weak resilience of the chemostat model to a species invasion with non-autonomous removal rates

Keywords: Non-autonomous Dynamics, Chemostat Model, Periodic systems, Asymptotically Periodic Solutions

Introduction

The chemostat model describes microbial ecosystems which are continuously fed by nutrients, as it can be found in natural environments, such as lakes, lagoons, wetlands... and in experimental or industrial bioreactors. Being open systems, chemostats are naturally subject to external perturbations such as species invasion.

In this paper, we focus on the classical chemostat model, for which the Competitive Exclusion Principle (CEP) holds [START_REF] Armstrong | Competitive exclusion[END_REF]. For a single limiting resource, the CEP states that no more than one species (generically) survives in the long term under constant fed conditions (input substrate concentration and flow rate of the incoming resource), see, e.g., [START_REF] Harmand | The chemostat, mathematical theory of continuous culture of micro-organisms[END_REF][START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF].

For dynamical systems, resilience is often described as the ability of a system to return to an original state (typically a steady state) after a transient perturbation [START_REF] Martin | Defining Resilience Mathematically: From Attractors To Viability, In "Viability and Resilience of Complex Systems, Understanding Complex Systems[END_REF]. In the present work, we wish to study the resilience of the chemostat model to invasions by other species considered as disturbances. In particular, we are interested in the possibility of extending the resilience domain using time-varying input conditions.

For the chemostat model, it has already been pointed out that non-constant removal rates could allow the coexistence of two species, under some precise integral conditions, see for instance [START_REF] Smith | Competitive coexistence in an oscillating chemostat[END_REF][START_REF] Butler | A mathematical model of the chemostat with periodic washout rate[END_REF][START_REF] Lenas | Coexistence of three competing microbial populations in a chemostat with periodically varying dilution rate[END_REF][START_REF] Peng | Global Attractivity in a Periodic Chemostat with General Uptake Functions[END_REF][START_REF] Nakaoka | How can three species coexist in a periodic chemostat? Mathematical and Numerical Study[END_REF][START_REF] Wang | Competition in a chemostat with Beddington DeAngelis growth rates and periodic pulsed nutrient[END_REF] for periodic removal rates or [START_REF] Lobry | Stability loss delay in the chemostat with a slowly varying washout rate[END_REF] for slow varying environments. The idea is to create a time-varying growth environment which alternates the favored species. Recently, the question of quantifying the excursions of the state variables in the chemostat model under periodic removal rates, has been studied in [START_REF] Bayen | Optimal periodic control for scalar dynamics under integral constraint on the input, Mathematical Control and Related Fields[END_REF][START_REF] Bayen | Optimal periodic control of the chemostat with Contois growth function[END_REF], but for the mono-specific case only. To our knowledge, it has not been investigated how to synthesize time-varying removal rates allowing resilience of a mono-specific chemostat system in presence of an invasion by a new species, in such a way that the resident species returns to the same density level than before invasion an infinite number of times. The design of such time-varying removal rate is precisely the matter of the present work.

The paper is structured as follows. In Section 2, we state the resilience problem in the context of the chemostat model with two species. In particular, we show the existence of a threshold on the level of the resident species above which resilience is lost that allows us to introduce a concept of weak resilience in a time varying context. In Section 3, we provide a construction of a time-varying removal rate which guarantees weak resilience in the sense given in Section 2. In Section 4, we show that there exist weakly resilient periodic solutions, and conjecture that there exists an unique periodic solution associated to the time-varying removal rate that we construct. Then, we show that any solution of the system associated with this time-varying removal rate converges to this periodic solution.

Assumptions and definition of weak resilience

We start by recalling the chemostat model with two microbial species of concentrations x 1 , x 2 , respectively, that compete for a single resource of concentration s:

   ṡ = -µ 1 (s)x 1 -µ 2 (s)x 2 + D(s in -s), ẋ1 = (µ 1 (s) -D)x 1 , ẋ2 = (µ 2 (s) -D)x 2 ,
(2.1) in which the yield coefficients are equal to one. Parameters D and s in are respectively the removal rate (imposed by the input flow rate) and the input concentration of the resource. Here, species 1 and 2 play the respective roles of the resident and invasive species, as described in the introduction: the chemostat is considered first with the species 1 (the resident) alone and at some time t 0 (chosen equal to 0 for simplicity), the invasive species 2 appears. In the sequel, we consider that the following assumption on the growth functions

µ i (•), i = 1, 2 in (2.1) is fulfilled. Assumption 2.1. The functions µ i (•) are of class C 1 , monotone increasing with µ i (0) = 0, i = 1, 2.
For i = 1, 2, the break-even concentration for species i (related to the parameter s in ) is defined as

λ i (D) := sup {s ∈ [0, s in ] ; µ i (s) < D} ∈ [0, +∞], i = 1, 2.
When λ i (D) < +∞, and because µ i are strictly monotone, one has

µ i (λ i (D)) = D ⇐⇒ λ i (D) = µ -1 i (D).
Consider now the competition between the two species. The CEP states that only the species that realizes the minimum of the numbers λ i (D), i = 1, 2 (with D > 0) has a non null concentration at steady state. For a given constant D with 0 < D < µ 1 (s in ), this means that species 1 is excluded by an invasion by a species 2 when λ 2 (D) < λ 1 (D), or equivalently that the state of (2.1) converges asymptotically to (λ 2 (D), 0, s in -λ 2 (D)). Let us now consider a situation for which the dominance of one growth function over the other one is alternated with respect to the level of the resource, in the following way.

Assumption 2.2. There exists s ∈ (0, s in ) such that

(µ 1 (s) -µ 2 (s))(s -s) > 0, s ∈ [0, s in ] and s = s. (2.2)
In the rest of the paper, we assume that this assumption holds true. Note that for Monod's growth functions [START_REF] Monod | La technique de culture continue: Théorie et applications[END_REF] (which are quite popular in microbiology and bio-processes), one has

µ i (s) = μi s K i + s , i = 1, 2.
Assumption 2.2 then amounts to have the following condition

0 < μ2 K 1 -μ1 K 2 < (μ 1 -μ2 )s in to be fulfilled. From (2.
2), one must have µ 1 (s)µ 2 (s) = 0, therefore we set

D := µ 1 (s) = µ 2 (s).
We formulate the problem of invasion of species 1 by species 2 as follows. Suppose that only species 1 is present in a bioreactor, at the initial time (i.e., x 2 = 0 in (2.1)). For a constant value of D, a straightforward analysis shows that a necessary condition for species 1 to set up at steady state is to have λ 1 (D) < s in . Then, the corresponding steady state is given by (λ 1 (D), s inλ 1 (D)), which is globally asymptotically stable in the (s, x 1 ) positive plane (see for instance [START_REF] Harmand | The chemostat, mathematical theory of continuous culture of micro-organisms[END_REF]). Moreover, if D is less than D, then the concentration of species 1 at steady state, x eq 1 , is necessarily such that x eq 1 ≥ x1 with

x1 := s in -µ -1 1 ( D) = s in -s > 0.
Suppose now that a new species (species 2) that fulfills Assumption 2.2, invades the growth vessel and that the removal rate D is less than D. Then, the new system (2.1) (with x 2 = 0) will not be able to return to the original state (from the CEP cited before, the state of (2.1) converges asymptotically to (λ 2 (D), 0, s in -λ 2 (D))), after this perturbation. This is why, we say that the dynamics is not resilient when the concentration x 1 at steady state, before the invasion, is above the threshold x1 . The question of interest is to investigate if, even so, the system can be resilient with x 1 above the threshold x1 , when considering time-varying removal rate D(•). First, one can easily check that the domain

{(s, x 1 , x 2 ) ∈ R 3 + ; s + x 1 + x 2 = s in },
is an invariant and attractive set for the dynamics (2.1) when D(•) is persistently exciting 1 . Assuming that system (2.1) is already at steady state before invasion, we shall consider in the sequel the reduced dynamics on this domain, that is,

ẋ = f (x, D) := (µ 1 (s in -x 1 -x 2 ) -D)x 1 (µ 2 (s in -x 1 -x 2 ) -D)x 2 , (2.3)
defined on the invariant set

S := {x ∈ R 2 + ; x 1 + x 2 ≤ s in }.
From now on, we consider D(•) as a control variable, i.e., as a measurable function of time taking values within some interval [D m , D M ] where the minimum and maximum dilution rates D m and D M satisfy the inequality:

0 < D m < D < µ 1 (s in ) < D M . (2.4)
Note that D M is large enough to have the possibility to drive solutions of (2.3) to the washout 2 of both species.

To introduce resilience, we consider a threshold

x r 1 ∈ (x 1 , s in ),
for species 1 aiming at keeping x 1 above x r 1 as much as possible. It is then natural to introduce the subset of S, K(x r 1 ), defined as: K(x r 1 ) := {x ∈ S ; x 1 ≥ x r 1 and x 2 > 0}, and to ask about weak invariance properties of K(x r 1 ) for the dynamics (2.3) in the context of viability theory [START_REF] Aubin | Viability Theory[END_REF]. Recall that given a controlled system ẋ = g(x, u) (with g : R n × R m → R n ) and given a closed subset K ⊂ R n , the viability kernel, denoted by V iab(K), is defined as the largest subset of initial states x 0 ∈ K for which there is an admissible control u(•) such that the unique solution x(•) of the dynamics associated with u and such that x(0) = x 0 , satisfies x(t) ∈ K for any time t ≥ 0 (see [START_REF] Aubin | Viability Theory[END_REF]). We then say that the viability kernel is weakly invariant. Going back to (2.3), we assume in the rest of the paper (in addition to Assumptions 2.1, 2.2) that the following assumption is fulfilled:

Assumption 2.3. The threshold x r 1 satisfies 0 < D m < µ 1 (s in -x r 1 ). (2.5)
Remark 2.1. The choice of the three parameters D m , D M , and x r 1 is crucial throughout this work. In this approach, note that we first chose D m , D M satisfying (2.4), and then we supposed that the threshold x r 1 satisfies (2.5). It is worth to mention that we could alternatively fix x r 1 ∈ (x 1 , s in ) and then choose the minimal and maximal dilution rates in such a way that (2.5) and the inequality µ 1 (s in ) < D M are verified.

One has the following property, in terms of viability analysis (hereafter cl S is the closure of a set S ⊂ R 2 ).

1 By persistently exciting, we mean that the non-negative function D(•) is such that +∞ 0 D(t)dt = +∞, see [START_REF] Bastin | On-line estimation and adaptive control of bioreactors[END_REF]. 2 This means that for D M sufficiently large (D M > µ 1 (s in )), solutions of (2.3) converge to the origin.

Lemma 2.4. The viability kernel V iab(cl K(x r 1 )) of the set K(x r 1 ) for

(2.3) satisfies V iab(cl K(x r 1 )) = [x r 1 , s in ] × {0}. Proof. Recall that x r 1 is such that x r 1 ∈ (x 1 , s in -λ 1 (D m
)) (by (2.5)). Assumptions 2.1, 2.2 and (2.4) also imply the inequality 0 < λ 2 (D m ) < λ 1 (D m ) < s. Take an initial condition (x 1,0 , x 2,0 ) in cl K(x r 1 ). If x 2,0 = 0, the solution of (2.3) verifies x 2 (t) = 0 for any t > 0 and any time-varying D(•). From (2.5), we can choose a constant D such that 0 < D m ≤ D < µ 1 (s inx r 1 ), and we observe that (2.3) with D(t) = D satisfies

x 1 = x r 1 ⇒ ẋ1 = (µ 1 (s in -x r 1 ) -D)x r 1 > 0 and x 1 = s in ⇒ ẋ1 = -Ds in < 0. Thus, this constant D prevents x 1 (•) to leave the interval [x r
1 , s in ] over [0, +∞). Assume now that x 2,0 > 0. Then, a solution x(•) of (2.3) associated with an admissible time-varying function D(•) verifies x 2 (t) > 0 for any time t ≥ 0. Suppose by contradiction that x(•) stays in K(x r 1 ) for any time t ≥ 0. Then one has s(t) < s inx r 1 < s for any time t ≥ 0. By Assumption 2.2, one has for any t ≥ 0, µ 1 (s(t))µ 2 (s(t)) < 0 and thus, we deduce the inequality

ṡ(t) > -µ 2 (s(t))x 1 (t) -µ 2 (s(t))x 2 (t) + D m (s in -s(t)) = (D m -µ 2 (s(t)))(s in -s(t)), t ≥ 0. Notice that any positive solution ζ(•) of ζ = (D m -µ 2 (ζ(t)))(s in -ζ(t)) converges to λ 2 (D m ) when t → +∞. From (2.5) one has λ 2 (D m ) < s in -x r
1 , hence there exist t 1 ≥ 0 and s ∈ (0, λ 2 (D m )) such that one has ζ(t) > s for any t ≥ t 1 . From the preceding inequality, we deduce that s(•) satisfies s(t) ≥ ζ(t) for any time t ≥ 0. We thus deduce that for any time t ≥ t 1 , one has

µ 1 (s(t)) -µ 2 (s(t)) ≤ c := min{µ 1 (σ) -µ 2 (σ) ; σ ∈ [s, s in -x r 1 ]} < 0. If we differentiate the function q 1 := x 1 /x 2 w.r.t. t, we find that q1 = µ 1 (s(t)) -µ 2 (s(t)) q 1 , with s(t) = s in -x 1 (t) -x 2 (t).
One then obtains q1 < c q 1 . Therefore q 1 decreases to zero and x 1 as well, leading to a contradiction. We conclude that the only solutions of (2.3) that stay in cl K(x r 1 ) for any time are the ones starting with x 2,0 = 0 as was to be proved. This lemma shows that for any given threshold x r 1 satisfying (2.5), i.e., such that x1 < x r 1 < s inλ 1 (D m ), (or equivalently λ 1 (D m ) < s inx r 1 < s), the dynamics (2.3) is not resilient for the domain K(x r 1 ) in presence of species 2. This is precisely our starting point to introduce the concept of weak resilience.

Definition 2.5. Let x r 1 ∈ (x 1 , s in -λ 1 (D m )). The system (2.
3) is said to be weakly resilient for the set K(x r 1 ) if for any initial condition in K(x r 1 ), there exists a time-varying function D(•) with values in [D m , D M ] such that the corresponding solution of (2.3) satisfies meas {t ≥ 0 ; x(t) ∈ K(x r 1 )} = +∞. Such a function D(•) will be called a weakly resilient removal rate.

This definition is related to the minimal time crisis of controlled dynamics, studied in [START_REF] Bayen | About Moreau-Yosida regularization of the minimal time crisis problem[END_REF][START_REF] Bayen | About the minimal time crisis[END_REF][START_REF] Bayen | Minimal time crisis versus minimum time to reach a viability kernel : a case study in the prey-predator model[END_REF][START_REF] Doyen | Scale of viability and minimal time of crisis[END_REF], although we do not look in this paper for control functions minimizing the time spent outside the set K(x r 1 ) over a given time period [0, T ].

Construction of a weakly resilient removal rate

The aim of this section is to propose a robust and systematic way to build a time-varying D(•) taking alternatively the values D m and D M , and that is weakly resilient for (2.3), without requiring a precise knowledge of the expressions of the growth functions µ i (•). Recall that we suppose Assumptions 2.1, 2.2, and 2.3 to be fulfilled and that the parameter x r 1 satisfies x r 1 ∈ (x 1 , s inλ 1 (D m )). We begin by giving the main result (in Section 3.1) which gives a construction of a weakly resilient D(•). Next, we provide some properties of the dynamical system(2.3) in the domain S for a constant D (D = D M and D = D m ). Finally, we give the proof of Proposition 3.1 at the end of this section.

D M D m D λ 2 (D m ) λ 1 (D m ) s s in µ 1 µ 2 Figure 1: Illustration of Assumption 2.3.

Synthesis of a weakly resilient removal rate

In the following Proposition, we propose a time-varying D(•) allowing the dynamics (2.3) to be weakly resilient for the set K(x r 1 ). For any ε > 0, we define the set E E := (0, x1 ] × (0, ε].

Proposition 3.1. For ε small enough and any initial condition

x 0 := (x 1,0 , x 2,0 ) ∈ K(x r 1 ) there exists a piecewise constant function D(•) which alternates the values D M , D m on time intervals [T i , T i+1 ), i ∈ N satisfying: T 0 = 0 < T 1 < • • • < T i < T i+1 < • • • and lim i→∞ T i = +∞, (3.1)
where x(•) is the unique solution of (2.3) associated with the time-varying D(•) such that

(i) if x(T i ) / ∈ E, one has D(t) = D M for t ∈ [T i , T i+1
) with T i+1 is defined as the first next entry time in E.

(ii) if x(T i ) ∈ E, one has D(t) = D m for t ∈ [T i , T i+1
), the trajectory x(•) enters to the set K(x r 1 ) in finite time and T i+1 is defined as the first next exit time from K(x r 1 ). Finally, the time-varying D(•) is a weakly resilient removal rate.

We begin by a lemma which describes the asymptotic behavior of (2.3) when D is constant. Lemma 3.1. Any solution of (2.3) in S with a constant removal rate D converges asymptotically to an equilibrium.

Proof. First, consider an initial condition on the axes that are invariant by (2.3). Then, the variable x i (i equal to 1 or 2) is solution of a scalar autonomous dynamics on the x i -axis. Therefore, either it converges to an equilibrium point on the axis, or it tends to infinity, which is not possible as the domain S is bounded.

Consider now a positive initial condition in the set S. The corresponding solution then remains in the positive orthant, and one can consider the variables

ξ i = ln(x i ), i = 1, 2, whose dynamics is ξi = F i (ξ) := µ i s in -e ξ1 -e ξ2 -D, i = 1, 2. For a constant D, one has divF (ξ) = i=1,2 ∂ ξi F i (ξ) = - i=1,2
µ i s ine ξ1e ξ2 e ξi < 0. By Dulac's criterion, the system has no closed orbit and by Poincaré-Bendixon Theorem (see [START_REF] Perko | Differential equations and dynamical systems[END_REF]). We can then conclude that solutions of (2.3) converge asymptotically to an equilibrium, since trajectories are bounded.

In the sequel, we denote by z(•, ζ, D) the unique solution of (2.3) (over R) for an initial condition z(0) = ζ ∈ S and a constant D ∈ {D m , D M }.

Properties of the reduced dynamics with constant D = D M

We now provide asymptotic properties of (2.3) when D is constant equal to D M , that are illustrated on Fig. Proof. From Assumption 2.3, the system (2.3) has a unique steady state (0, 0) in S. From Lemma 3.1, we deduce that it is globally asymptotically stable on S. Consider a positive initial condition. From the expression of the dynamics (2.3), the solution x(•) = z(•, ζ, D M ) is clearly positive for any t ≥ 0, and by Assumption 2.3, one gets ẋi (t) < 0 for any t ≥ 0, thus z i (•, ζ, D M ) is decreasing for i = 1, 2. Consider then the function q 2 := x 2 /x 1 . A straightforward computation of its derivative yields

q2 = µ 2 (s in -x 1 (t) -x 2 (t)) -µ 1 (s in -x 1 (t) -x 2 (t)) q 2 .
By Assumption 2.2, one has µ 2 (s in )µ 1 (s in ) < 0. Since x(•) converges to 0, we deduce that there exist η > 0 and t M > 0 such that q2 (t) < -ηq 2 (t) for any time t > t M . This proves that q 2 (•) converges asymptotically to 0 and that trajectories are tangent to the x 1 -axis at (0, 0).

We are now in a position to introduce the following notation that will be used hereafter (see Fig. 2):

• Consider the point P r := (x r 1 , s inx r 1 ) on the boundary of S and set

x(•) := z(•, P r , D M ).

• The forward semi-orbit of (2.3) with D = D M passing through P r is denoted by (see Fig. 2):

γ+ := {x(t) ; t ≥ 0}.

• In view of Lemma 3.2, x1 (•) is decreasing and thus reaches x1 in finite time. Hence, there are t > 0 and δ > 0 satisfying: t := inf{t > 0; x1 (t) < x1 } and δ := x2 ( t).

Properties of the reduced dynamics with constant D = D m

We now turn to asymptotic properties of (2.3) with D = D m . In the sequel, we shall use the notation (i) It admits exactly three equilibria in S : E 0 := (0, 0), E m 1 := (x m 1 , 0) and E m 2 := (0, x m 2 ) which are respectively an unstable node, a saddle point, and a stable node.

x m i := s in -λ i (D m ), i = 1, 2. The scalar product in R 2 is written a • b with a, b ∈ R
(ii) One has x m 1 < x m 2 and the "strip" S m 1,2 defined as S m 1,2 := {x ∈ S ;

x 1 + x 2 ∈ [x m 1 , x m 2 ]} , is invariant (by (2.3) with D = D m ).
(iii) The edge (0, s in ] × {0} is the stable manifold of (2.3)

with D = D m at E m 1 on S. The unstable manifold at E m 1 in the domain S is denoted by W u (E m 1 ): it connects E m 1 to E m 2 and satisfies W u (E m 1 ) ⊂ S m 1,2 .
Proof. Inequality (2.5) and the monotonicity of the functions µ i (•) (Assumption 2.1) imply that there are two equilibria of (2.3) in S that are distinct of the origin and on the axes. They are uniquely defined by E m

1 and E m 2 . The Jacobian matrix at E m 1 and E m 2 are respectively given by

J(E m 1 ) := -µ 1 (sin -x m 1 )x m 1 -µ 1 (sin -x m 1 )x m 1 0 µ2(sin -x m 1 ) -Dm , J(E m 2 ) := µ1(sin -x m 2 ) -Dm 0 -µ 2 (sin -x m 2 )x m 2 -µ 2 (sin -x m 2 )x m 2 .
The point E m 1 is a saddle point because the eigenvalues of J(E m 1 ) are of opposite sign, E m 2 is a stable node because the eigenvalues of J(E m 2 ) are negative, and E 0 is clearly an unstable node which proves (i). It is worth noting that the unstable manifold

W u (E m 1 ) necessarily connects E m 1 to E m 2 by Lemma 3.1. Let us now prove (ii). From Assumption 2.2, one has x m 1 < x m 2 . When x ∈ S m 1,2 is such that x 1 + x 2 =
x m 1 , one has ẋ1 = 0 and ẋ2 > 0 whereas if

x 1 + x 2 = x m
2 , one has ẋ1 < 0 and ẋ2 = 0. Hence S m 1,2 is invariant. Let us finally prove (iii). The positive half axis x 1 > 0 is clearly the stable manifold of (2.3) with D = D m at E m 1 . Consider a non-null eigenvector v + of J(E m 1 ) associated with the positive eigenvalue µ 2 (s in -x m 1 )-D m , and let n be an outward normal to S m 1,2 at E m 1 . A straightforward calculation yields

v + = -1 1 + µ2(λ1(Dm))-Dm x m 1 µ 1 (λ1(Dm)) , n = -1 -1 , implying that v + • n = - µ 2 (λ 1 (D m )) -D m x m 1 µ 1 (λ 1 (D m )) < 0.
One then concludes that the vector v + points inward S m 1,2 at E m 1 . On another hand, from the Theorem of the stable and unstable manifolds [START_REF] Perko | Differential equations and dynamical systems[END_REF], we know that

W u (E m 1 ) is tangent to v + at E m 1 . Therefore, there is a neighborhood V of E m 1 in S m 1,2 such that W u (E m 1 ) ∩ V ⊂ S m 1,2 , but, as S m 1,2 is invariant, we conclude that W u (E m 1 ) ⊂ S m 1
,2 as was to be proved. Recall that the unstable manifold W u (E m 1 ) is a trajectory of (2.3) with D = D m , and that (2.3) satisfies ẋ1 < 0 on int S m 1,2 . So, W u (E m 1 ) can be parametrized as a function

x 1 → w u (x 1 ), x 1 ∈ [0, x m 1 ]
. Hereafter, hyp(w u ) stands for the hypograph of w u and let D ⊂ S be defined as (see Fig. 3):

D := hyp(w u ) ∩ S.
Note that the domain D is necessarily forward and backward invariant (for (2.3) with D = D m ) as its boundary is a locus of trajectories. Similarly as with D = D M , let us introduce the following notation (see Fig. 3):

• Since there is a unique intersection point between W u (E m 1 ) and the line {x 1 = x r 1 }, we set: x2 = w u (x r 1 ),

one can also write (

x r 1 , x2 ) = W u (E m 1 ) ∩ {x 1 = x r 1 }.
• From Assumption 2.3, one has x r 1 < x m 1 , thus we can fix a point P r := (x r 1 , x2,0 ) ∈ D such that:

0 < x2,0 < min(δ, x m 1 -x r 1 ), (3.3)
where x2,0 is small enough to ensure P r ∈ D.

• The backward semi-orbit of (2.3) with D = D m passing though P r is denoted by: γ-:= {x(t) ; t ≤ 0}, where x(•) = z(•, P r , D m ).

• Finally, define a positive parameter η > 0 as η := x m 1x r 1 -x2,0 .

s in x 2 x 1 x m 1 x m 1 x m 2 W u D s in x r 1 x 2 x 1 x2 x m 1 x m 1 W u γ- x m 2 D P r
(x 1 , x 2 ) ∈ {x ∈ D ; x 1 > 0 and x 1 + x 2 < x m 1 } ⇒ ẋ1 > 0, (x 1 , x 2 ) ∈ {x ∈ D ; x 1 > 0 and x 1 + x 2 > x m 1 } ⇒ ẋ1 < 0. (ii) As W u (E m 1
) is a trajectory, any solution of (2.3) with D = D m crosses the line {x 1 = x r 1 } at some point

(x r 1 , x 2 ) such that x 2 < x2 . (iii) Since D is backward invariant by (2.3) with D = D m and P r ∈ D, the inclusion γ-⊂ D is fulfilled.
The next lemma will be useful to define a small ε > 0 and times T and Ť (see Remarks 3.2 and 3.3 below).

Lemma 3.4. The curves γ+ and γintersect in the domain (0, x1 ) × (0, x2,0 ).

Proof. Consider the variable q1 = x1 /x 2 on the positive orthant. As previously, one has

q1 = µ 1 (s in -x1 (t) -x2 (t)) -µ 2 (s in -x1 (t) -x2 (t)) q1 .
From Assumption 2.2, µ 1 (s in )µ 2 (s in ) > 0. As x(t) → (0, 0) when t → -∞, there exist ť < 0 and c > 0 such that

µ 1 (s in -x1 (t) -x2 (t)) -µ 2 (s in -x1 (t) -x2 (t)) > c , t ≤ ť,
which shows that q1 (t) tends to 0 when t tends to -∞. Therefore γ-is tangent to the x 2 -axis at E 0 . As γ+ is tangent to the x 1 -axis at E 0 (Lemma 3.2), we deduce that the curve γ-is above the curve γ+ , in a neighborhood of the point E 0 in S. Since x1 (0) = x r 1 > x1 and x1 (t) tends to 0 when t tends to -∞, there exists t < 0 such that x1 ( t) = x1 with t = sup{t < 0; x1 (t) = x1 }.

This means that x1 (t) > x1 for all t ∈ ( t, 0]. As x2 (•) is increasing on D (cf. Remark 3.1), one has x2 ( t) < x2 (0) and from the choice of x2,0 = x2 (0) < δ (see condition (3.3)), one gets x2 (t) < δ, for all t ≤ 0 and in particular at t = t. The point (x 1 , x2 ( t)) of γ-is thus below (x 1 , δ), which belongs to γ+ . Therefore, γ+ and γ-have to cross at some point x( Ť ) with Ť < t, which verifies 0 < x1 ( Ť ) < x1 and 0 < x2 ( Ť ) < x2,0 .

Remark 3.2. Since γ-⊂ D (cf. Remark 3.1), the intersection between γand γ+ is also contained in D.

Lemma 3.4 implies that for each choice of the point P r , there is an intersection point

P ε := (x 1 , ε) ∈ D,
between γ+ and γ-such that 0 < x 1 < x1 and 0 < ε < x2,0 , see Fig. 4. By construction, there are T > 0 and Ť > 0 such that

P ε = x( T ) = x(-Ť ). Recall that E is by definition E := (0, x1 ] × (0, ε].
Since x1 < x r 1 and ε < x2,0 , the corner point (x 1 , ε) of E is below the curve γ-. Thus, one has also the inclusion E ⊂ D (because γ-⊂ D), see Fig. 4.

Remark 3.3. Given D m , D M and x r
1 that fulfill Assumption 2.3, the parameter ε can be chosen arbitrarily small taking the parameter x2,0 in the definition of P r small enough (recall (3.3)).

x m 1 ε K(x r 1 ) x1 x 1 γ+ γ- W u x m 1 x r 1 x 1 x 2 s in s in

Proof of Proposition 3.1

To help the reader, we provide in Appendix 1 a list of the notations used in Section 3. We start by giving the following definition.

Definition 3.5. The "southeast" order in R 2 (denoted by ) is defined as

∀(x, y) ∈ R 2 , x y ⇐⇒ {x 1 ≤ y 1 , x 2 ≥ y 2 }.
Notice that the dynamics (2.3) is competitive, and therefore (2.3) preserves the order (see [START_REF] Smith | Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems[END_REF]): for any admissible time-varying function D(•), one has:

∀(ζ 1 , ζ 2 ) ∈ S, ζ 1 ζ 2 =⇒ z(t, ζ 1 , D(t)) z(t, ζ 2 , D(t)), ∀t ≥ 0.
(3.4)

We have now given the necessary definitions and properties so that we can give to the proof of Proposition 3.1.

Proof. For clarity, we present the proof in several steps.

Step 1. Fix an initial condition x 0 ∈ K(x r 1 ) with x 2,0 < s inx r 1 . As x 0 / ∈ E, we set D = D M . Accordingly to Lemma 3.2, x(•) converges asymptotically to E 0 and thus T 1 := inf{t > 0 ; x(t) ∈ E} is well defined. Notice that there are two ways to reach E:

(i) x 1 (T 1 ) = x1 and x 2 (T 1 ) < ε, (ii) x 1 (T 1 ) ≤ x1 and x 2 (T 1 ) = ε.
Let us show that in both cases one has x 1 (T 1 ) > x 1 . As x1 > x 1 (recall Lemma 3.4), we get x 1 (T 1 ) > x 1 in case (i). In case (ii), one can also write T 1 = inf{t > 0 ; z 2 (t, x 0 , D M ) ≤ ε}. The order property (3.4) then implies that z(t, P r , D M ) z(t, x 0 , D M ), t ≥ 0.

As z( T , P r , D M ) = x( T ) = 1 , ε), we deduce that

z 1 ( T , x 0 , D M ) ≥ x 1 , z 2 ( T , x 0 , D M ) ≤ ε.
From the definition of T 1 , one has T 1 < T (one gets T 1 = T if x 0 = P r ) and thus

x 1 (T 1 ) = z 1 (T 1 , x 0 , D M ) ≥ x 1 . If x 1 (T 1 ) = x 1 , one should then have x 2 (T 1 ) = ε, that is, z(•, x 0 , D M ) = z(•, P r , D M )
which is not true. Hence, we obtain as well x 1 (T 1 ) > x 1 as was to be proved. In addition, notice that P ε x(T 1 ). Since E ⊂ D, the point x(T 1 ) necessarily satisfies x(T 1 ) ∈ D.

Step 2. At t = T 1 , we set D = D m . We use again the order property (3.4):

P ε x(T 1 ) =⇒ z(t, P ε , D m ) z(t, x(T 1 ), D m ), ∀t ≥ 0, (3.5)
and, as we have shown that x 1 (T 1 ) > x 1 , we obtain the inequalities

z 1 ( Ť , x(T 1 ), D m ) > z 1 ( Ť , P ε , D m ) = x r 1 , z 2 ( Ť , x(T 1 ), D m ) ≤ z 2 ( Ť , P ε , D m ) = x2,0 .
Therefore one has x(T 1 + Ť ) = z( Ť , x(T 1 ), D m ) ∈ int K(x r 1 ). One can then define a time T1 as:

T1 := inf{t > T 1 ; x 1 (t) > x r 1 },
which is such that T1 ∈ (T 1 , T 1 + Ť ). From the monotonicity of x 2 (•) in the set D (cf. Remark 3.1), one obtains

x 2 ( T1 ) < x 2 (T 1 + Ť ) = z 2 ( Ť , x(T 1 ), D m ) ≤ x2,0 .
As x(T 1 ) belongs to the set D with x 1 (T 1 ) > 0, z(•, x(T 1 ), D m ) converges asymptotically to the equilibrium E m 2 that lies on the x 2 -axis (cf. Remark 3.1). The time

T 2 > T 1 + Ť such that x 1 (T 2 ) = x r 1 , where x(t) := z(t -T 1 , x(T 1 ), D m ), t ≥ T 1 , is thus well defined. Moreover one has x(T 2 ) ∈ K(x r 1 ) and x(T 2 ) ∈ D as D is invariant by (2.3) (for D = D m ).
Step 3. At time T 2 , we have shown that x(T 2 ) belongs to the set K(x r 1 ), and also to the set D which implies that x 2 (T 2 ) < s inx r 1 . Therefore, we can consider x(T 2 ) as a new initial condition and apply iteratively the results of steps 1 and 2, defining an increasing sequence of times (T i ) i∈N . For i = 2k + 1 (with k ∈ N), one has x(T 2k+1 ) ∈ E and, as in step 2, we can define T2k+1 ∈ (T 2k+1 , T 2k+1 + Ť ) such that

x 1 ( T2k+1 ) = x r 1 with T2k+1 := inf{t > T 2k+1 ; x 1 (t) > x r 1 }. (3.6)
As shown in step 2, we necessarily have

x 2 ( T2k+1 ) < x2,0 . (3.7)
Because T 2k+2 > T 2k+1 + Ť , we get that lim i→+∞ T i = +∞ which concludes that property (3.1) is fulfilled.

Note that if we chose x 0 = P r then one obtains T 1 = T (and then by the unicity of the solution x(T 1 ) = P ε ). Furthermore, in this case, the time T1 is such that T1 = T 1 + Ť with x( T ) ∈ K(x r 1 ), since x(T 1 ) = x(-Ť ). The time T 2 > T1 that is the first exit time from K(x r 1 ) is well defined with x 2 (T 2 ) < s inx r 1 .

Step 4. We now show that the time spent by x 1 (•) above the threshold x r 1 is of infinite measure. For each k ∈ N, one has, from the definition of

T 2k+2 : meas{t ∈ [T 2k+1 , T 2k+2 ] ; x 1 (t) > x r 1 } = T 2k+2 -T2k+1 .
From Remark 3.1, one has ẋ1 (t) > 0 when s(t) = s inx 1 (t)x 2 (t) > λ 1 (D m ). At time T2k+1 , the inequality (3.7) implies that s( T2k+1 ) > s inx r 1 -x2,0 . We deduce that

s( T2k+1 ) -λ 1 (D m ) = s( T2k+1 ) -(s in -x r 1 -x2,0 ) + η > η > 0. (3.8)
Next, let us define

τ k := sup{θ > 0 ; s( T2k+1 + θ) > λ 1 (D m )}.
Then, one has necessarily x 1 (t) > x r 1 for any t ∈ [ T2k+1 , T2k+1 + τ k ] and one obtains the inequality

T 2k+2 -T2k+1 > τ k .
Let us now give a lower bound on the value of τ k . From equations (2.3), the following properties hold true:

x 1 > x r 1 , s > λ 1 (D m ) ⇒    s = s in -x 1 -x 2 < s in -x r 1 , x 1 = s in -s -x 2 < s in -λ 1 (D m ) = x m 1 , x 2 = s in -s -x 1 < s in -λ 1 (D m ) -x r 1 = x m 1 -x r 1 ,
from which one can obtain a lower bound on the speed at which the variable s decreases (as long as s is above λ 1 (D m ) and x 1 above x r 1 ):

ṡ = -µ 1 (s)x 1 -µ 2 (x)x 2 + D m (s in -s) ≥ -µ 1 (s in -x r 1 )x m 1 -µ 2 (s in -x r 1 )(x m 1 -x r 1 ) := -c .
Notice that c > 0. By integrating the above inequality, one can conclude that s stays above λ 1 (D m ) for a duration larger than (s( T2k+1 )λ 1 (D m ))/c . Thanks to (3.8), we can thus write

τ k > M := η c > 0,
where M > 0 does not depend on k. Finally, we have shown that for each k ∈ N, one has

meas{t ∈ [0, T 2k+2 ] ; x 1 (t) > x r 1 } > kM,
which shows that the time-varying D(•) is a weaky resilient removal rate as was to be proved.

x 2 ( T2k+1 ) < x 2 (T 2k+2 ) < x2 < x m 2 -x r 1 .
(3.9)

Convergence to periodic solutions

The goal of this section is to show that the weakly resilient removal rate defined by Proposition 3.1 generates asymptotically positive periodic solutions of system (2.3). Before addressing this point, we shall first prove the existence of a periodic solution of (2.3) associated with the time-varying D(•) given by Proposition 3.1.

To this end, let us introduce an operator O as

O : (0, s in -x r 1 ] → (0, s in -x r 1 ] x 2,0 → x 2 (T 2 )
where x(•) = (x 1 (•), x 2 (•)) is the unique solution of (2.3) for the initial condition (x r 1 , x 02 ) and the timevarying D(•) given by Proposition 3.1, parameters D m , D M , and P r being fixed. Notice that times T i , i ≥ 1, introduced in Proposition 3.1 depend on x 2,0 (in particular T 2 ). Hence, this operator slightly differs from the Poincaré map used for instance in [START_REF] Mottoni | Competition systems with periodic coefficients: a geometric approach[END_REF][START_REF] Smith | The theory of chemostat, dynamics of microbial competition[END_REF] for finding periodic solutions of dynamical systems, for which the period is fixed beforehand. We shall next examine properties of the operator O. Doing so, let us introduce the following notation (see Fig. 5):

• Denote by γ and γ, the orbits of (2.3) with D = D m passing respectively by (x r 1 , x m 1x r 1 ) and Pr . • Observe that γ is tangent to the segment {x

1 = x r 1 } ∩ S at (x r 1 , x m 1 -x r 1 ). Because x(•) converges to E m 2 (Lemma 3.
3), there are exactly two intersection points between γ and {x 1 = x r 1 } ∩ S, namely Pr and Qr :

= (x r 1 , x 2 ) with x 2 > x m 1 -x r 1 . (4.1)
• In the sequel, we denote by J the interval J

:= [x 2 , x2 ] (recall (3.2)). γ γ x m 1 x 2 x r 1 x m 1 x 2 s in Qr x 1

Properties of the operator O

In this section, we prove that O is continuous and decreasing. The continuity of O will follow from the continuity property of the first entry time into a set, that is related to Petrov's condition (see Appendix 2). For initial condition (x r 1 , x 2,0 ) with x 2,0 ∈ (0, S inx r 1 ], we denote T 1 (x 2,0 ), T 2 (x 2,0 ) the times T 1 , T 2 given by Proposition 3.1. Proof. For a given x 2,0 ∈ (0, s inx r 1 ], denote by x(•) the unique solution of (2.3) for the initial condition (x r 1 , x 02 ) and the time-varying D(•) given by Proposition 3.1.

Let us prove (i). For t ∈ [0, T 1 ], one has D(t) = D M and T 1 is defined as the first time that x(•) reaches the set T 1 := {x ∈ S ; x 1 ≤ x1 and x 2 ≤ ε}.

For t ∈ (T 1 , T 2 ), one has D(t) = D m . We know that there exists T1 ∈ (T 1 , T 1 + Ť ) such that x(•) crosses the line {x 1 = x r 1 } at time T1 with x 2 ( T1 ) < x2,0 (see step 2 of the proof of Proposition 3.1). Therefore, the state x( T1 ) is below the point P r . As x(•) cannot cross the orbit of P r (denoted γ) on ( T1 , T 2 ), x(•) crosses the line

{x 1 = x r 1 } at some time T 2 such that x 2 (T 2 ) > x 2 . (4.2)
Hence, T 2 can be defined as the first time t > T 1 + Ť such that x(•) reaches the set

T 2 := {x ∈ S ; x 1 ≤ x r 1 and x 2 ≥ x 2 },
at time t. On the positive set S + := {x ∈ S ; x 1 > 0 ; x 2 > 0}, that is invariant by (2.3), we are in a position to introduce the first entry time functions:

R 1 (x 0 ) := inf {t ≥ 0 ; z(t, x 0 , D M ) ∈ T 1 } , R 2 (x 0 ) := inf {t ≥ 0 ; z(t, x 0 , D m ) ∈ T 2 } ,
with x 0 ∈ S + . Then, for x 2,0 ∈ (0, s inx r 1 ], Proposition 3.1 allows to write the composition

O(x 2,0 ) = z 2 (T 2 -T 1 -Ť , x(T 1 + Ť ), D m ), (4.3) with        T 1 := R 1 (x r 1 , x 2,0 ), x(T 1 ) := z(T 1 , (x r 1 , x 2,0 ), D M ), x(T 1 + Ť ) := z( Ť , x(T 1 ), D m ), T 2 := T 1 + Ť + R 2 (x(T 1 + Ť )),
thanks to the definitions of T 1 , T 2 , and Ť . From the continuous dependency of an ODE w.r.t. initial conditions, (see, e.g., [START_REF] Perko | Differential equations and dynamical systems[END_REF]), the maps x 0 → z(t, x 0 , D) (for a fixed t) and t → z(t, x 0 , D) (for a fixed Lipschitz continuity of R 1 . Observe first that T 1 is convex, hence the (convex) normal cone to T 1 at some point x ∈ S + of its boundary is given by the expression

N T1 (x) =    R + × {0}, x 1 = x1 , x 2 < ε, {0} × R + , x 1 < x1 , x 2 = ε, R + × R + , x 1 = x1 , x 2 = ε.
Then, we can easily check that

{x 1 = x1 , x 2 ≤ ε} ⇒ f 1 (x, D M ) < φ 1 := (µ 1 (s in ) -D M )x 1 < 0, {x 1 ≤ x1 , x 2 = ε} ⇒ f 2 (x, D M ) < φ 2 := (µ 2 (s in -ε) -D M )ε < 0.
From the preceding inequalities, we deduce that for any point x on the boundary of T 1 in S + , one has

f (x, D M ) • ν < min(φ 1 , φ 2 ) ν , ν ∈ N T1 (x) \ {0}.
This allows us to conclude that R 1 is Lipschitz continuous over S + , thanks to Theorem 4.1.

Let T i be the sequence of times given by Proposition 3.1 for the initial condition x 0 = (x r 1 , x 2,0 ) and posit

F (t, x) = f (x, D M ), t ∈ [T 2k , T 2k+1 ) f (x, D m ), t ∈ [T 2k+1 , T 2k+2 ) (k ∈ N)
Clearly, on has T i+2 = T i + T 2 i.e. F is a T 2 -periodic dynamics.

We define now the following functions

F (τ, x) = T i+1 -T i T i+1 -T i F (τ, x), τ ∈ [T i , T i+1 ) g(t) = T i + T i+1 -T i T i+1 -T i (t -T i ), t ∈ [t i , T i+1 )
Clearly, the solution x(•) of ẋ = F (t, x) satisfies x(t) = x(g(t)) for any t ≥ 0, where x(•) is solution of Then, one has

lim i→+∞ T i+1 -T i T i+1 -T i = 1
and deduce that F is an asymptotically periodic dynamics with F as limit (see Definition 4.2 in Appendix 3). Moreover, one has lim k→+∞ x(kT 2 ) = x(T 2k ) = x 0 .

Therefore, we can apply the Theorem 4.3 (from [START_REF] Zhao | Asymptotic behavior for asymptotic periodic semiflows with applications[END_REF], recalled in Appendix 3) which gives lim t→+∞ x(t)x (t) = 0 and we have obtained thus the following limit lim t→+∞ x(t)x (g(t)) = 0.

We show now that the time shift σ(t) := tg(t) admits a limit σ. Notice that one has δ(T i ) = T i -T i for any i. Therefore, it is enough to prove that the sequence u k = T 2k -T 2k converges.

As O is contractive on J, there exists α ∈ (0, 1) such that As α < 1, on obtains that |u n+ku n | tends to 0 when n tends to ∞ uniformly in k. We conclude that the Cauchy sequence u k admits a limit, which gives the asymptotic shift σ.

2 . 3 . 2 .

 232 Lemma Any solution z(•, ζ, D M ) of (2.3) with ζ ∈ S converges asymptotically to the origin. Moreover, for any positive initial condition ζ in S, z 1 (•, ζ, D M ) and z 2 (•, ζ, D M ) are decreasing functions and the trajectory converges to the origin tangentially to the x 1 -axis.

Figure 2 :

 2 Figure 2: Phase portrait of (2.3) with constant D = D M and plot of the points P r , (x 1 , δ), and the semi-orbit γ+ .

2 Lemma 3 . 3 .

 233 and a denotes the euclidean norm of a vector a ∈ R 2 . The system (2.3) with D = D m possesses the following properties.

Figure 3 :

 3 Figure 3: Left: phase portrait of (2.3) with constant D = D m in the domain D whose boundary is the unstable manifold W u (E m 1 ). Right: plot of the point P r and the semi-orbit γ-.

Figure 4 :

 4 Figure 4: Plot of the set E := (0, x1 ] × ε] (in red) and the intersection point P ε = (x 1 , ε) between γ+ and γ-. The set K(x r 1 ) is depicted in blue.

Remark 3 . 4 .

 34 (i) In the proof of Proposition 3.1 (step 1), we have seen that x 1 (T 2k+1 ) > x 1 for any k ∈ N. Therefore, Proposition 3.1 remains valid if the set E is replaced byẼ := (x 1 , x1 ] × (0, ε].(ii) Notice also that the trajectory given by Proposition 3.1 reaches the set D at some time t ≤ T 1 , and then remains in this set for any future time. Indeed, from Remark 3.1, the trajectory belongs to D on any time interval [T 2k+1 , T 2k+2 ] (with D = D m ). On a time interval [T 2k+2 , T 2k+3 ], we have set D = D M , and we have seen in Lemma 3.2 that x 1 (•) and x 2 (•) are decreasing (with D = D M ). So, the trajectory also remains in D for t ∈ [T 2k+2 , T 2k+3 ]. One then concludes that the trajectory remains in D as well as in {x 1 > x 1 } (thanks to point (i) above) over [T 1 , +∞) (iii) Finally, note that at times T2k+1 and T 2k+2 , one has x 1 ( T2k+1 ) = x 1 (T 2k+2 ) = x r 1 and

Figure 5 :

 5 Figure 5: Plot of the orbits γ and γ and the point Qr .

Proposition 4 . 1 .

 41 The time-varying D(•) contructed in Proposition 3.1 fulfills the following continuity properties:(i) Times T 1 and T 2 are Lipschitz continuous functions of initial x 2,0 .(ii) The operator O is Lipschitz continuous and takes values in J.

  x 0 ) are Lipschitz continuous, given a constant D ∈ [D m , D M ]. Therefore, proving the Lipschitz continuity of T 1 and T 2 w.r.t x 2,0 essentially requires to prove the Lipschitz continuity of R 1 and R 2 over the set S + . Notice first that for constant D = D M , resp. D = D m , any solution in S + converges asymptotically to a steady state that belongs to the interior of T 1 , resp. T 2 . Therefore, the targets T 1 and T 2 can be reached in a finite horizon from any initial condition in S + , and thus, R 1 , R 2 are well defined with finite values in S + . To prove their Lipschitz continuity, we shall use Theorem 4.1 recalled in Appendix 2, showing that the inward pointing condition (4.8) is fulfilled on the boundary of T 1 and T 2 in S + .

  = F (τ, x(τ )), x(0) = x(0).Thanks to the continuity property of T i (see Proposition 4.1 and Remark 4.1) and (4.5), one concludes thatlim k→+∞ T 2k+i -T 2k = lim k→+∞ T i (x 2 (T 2k )) = T i (x 2,0 ), i = 1, 2, (4.6) which gives lim k→+∞ T 2(k+1) -T 2k = T 2 (x 2,0 ) = T 2 ,(4.7)and lim k→+∞ T 2k+i+1 -T 2k+i = T i+1 (x 2,0 ) -T i (x 2,0 ), i = 0, 1.

|x 2 ( 1 i=0α

 21 T 2(k-1) )x 2,0 | = |O(x 2 (T 2(k-2) )x 2,0 | ≤ α k-1 |x 2,0x 2,0 |, k > 1As the map x 2,0 → T 2 is Lipschitz continuous (cf Proposition 4.1), say of rank L, one obtains|u ku k-1 | = |T 2 (x 2 (T 2(k-1) ) -T 2 | ≤ Lα k-1 |x 2,0x 2,0 |and let us finally show that u k is a Cauchy sequence. For any n > 1 and k > 1, one has|u n+ku n | ≤ L k-n+i |x 2,0x 2,0 | = Lα n 1α k | 1α |x 2,0x 2,0 |

Acknowledgments

The authors thank the AVERROES program between Algeria and France for the PhD grant of F.-Z. Tani. The first author would like to thank INRA Montpellier and the UMR MISTEA for providing a half year delegation during the academic year 2017-2018. This work was supported by the LabEx NUMEV incorporated into the I-Site MUSE (AAP2017-2-08).

Lipschitz continuity of R 2 . Observe that T 2 is also convex, hence the (convex) normal cone to T 2 at some point x ∈ S + of its boundary is given by the expression

One can easily that the following properties are fulfilled:

As a consequence of (3.9), (4.1) and (4.2), note that one has x m 1x r 1 < x 2 < x m 2x r 1 . This allows us to conclude that ψ 1 < 0 and that ψ 2 > 0 yielding the inequality

for any x on the boundary of the set T 2 in S + . This proves the Lipschitz continuity of R 2 on S + using again Theorem 4.1. We conclude that both T 1 and T 2 are Lipschitz continuous w.r.t x 2,0 .

Let us prove now (ii). Recall that O can be written as function of T 1 and T 2 (see (4.3)). As T 1 and T 2 are Lipschitz continuous w.r.t x 2,0 then O as well. Combining (3.9) and (4.2) gives that O(x 2,0 ) ∈ J as was to be proved. Proof. Take two points x ± 2,0 ∈ (0, s inx r 1 ] such that x - 2,0 < x + 2,0 , and let us show that O(x - 2,0 ) > O(x + 2,0 ). Denote by x + (•), x -(•) the solutions generated by the time-varying D(•) given by Proposition 3.1 and for the initial conditions x + (0) = (x r 1 , x + 2,0 ) and x -(0) = (x r 1 , x - 2,0 ) respectively. One can then write x + (0) x -(0). For convenience, we denote by T + 1 , T + 2 and T - 1 , T - 2 the times T 1 , T 2 (as in Proposition 3.1) associated with x + (•) and x -(•) respectively, and let us set

First, let us note that the time-varying removal rate given by Proposition 3.1 satisfies D(t) = D M for both trajectories in a neighborhood of t = 0. Using the order property (3.4), one can then write

Thanks to this property, we must have

implying a contradiction with (4.4)). Therefore one gets (recall that trajectories with D = D M decrease), we obtain the inequality

Since the orbits of (2.3) with D = D M do not intersect, we also obtain that

, the time needed by x + (•) to reach the liner

We now consider x + (T + 1 ) and x -(T - 1 ) as initial conditions for (2.3) with D = D m . Then one gets

In the same way as previously, we deduce the inequality:

Since the orbits of (2.3) with D = D m do not intersect, one necessarily has x + 2 ( T + 1 ) > x - 2 ( T -1 ) (and also

). Finally, D is constant equal to D m for both trajectories until the first instant at which one of the two trajectory leaves the set K(x r 1 ). Since

Existence and attractivity of periodic solutions

In this section, we study how for any initial condition, the time-varying removal rate D(•) given in Proposition 3.1 allow system (2.3) to synchronize with a periodic solution (i.e. any solution of (2.3) associated with D(•) converges asymptotically to a periodic solution).

Existence of periodic solutions

The existence of a weakly resilient periodic trajectory follows from the previous results about the operator O. given by Proposition 3.

From Propositions 4.1 and 4.2, ϕ(•) is continuous and decreasing. Moreover, it verifies ϕ(x 2 ) > 0 for x 2 < x 2 and ϕ(x 2 ) < 0 for x 2 > x2 (because O is with values in J). We can then conclude that ϕ(•) possesses a unique zero in the interval (0, s inx 1,r ], or equivalently that there exists a unique fixed point x 2,0 of O. The solution x (•) for the initial condition (x r 1 , x 2,0 ) verifies x (T 2 ) = x (0), where T -2 is equal to the time T 2 generated by the time-varying D(•) given in Proposition 3.1, and is thus T 2 -periodic. We conclude that x (•) is the unique periodic solution such that x 1 (0) = x r 1 and O(x 2 (0)) = x 2 (0).

Attractivity of the periodic solution

Due to the particular structure of the non-autonomous dynamics (the times T i are not known explicitly), it appears that determining explicitly a bound on the Lipschitz rank of O is quite difficult. However, in all the simulations we performed, the operator O appears to be contractive, providing ε to be sufficiently small. We thus posit the following conjecture.

Conjecture 4.1. For ε sufficiently small, the operator O is contractive on J.

Proposition 4.3. Under the conjecture, for ε sufficiently small an any initial condition in K(x r 1 ), the solution x(•) generated by the time-varying D(•) given in Proposition 3.1 converges asymptotically to the periodic solution x (•) up to a time shift σ: lim

Proof. Fix an initial condition in K(x r 1 ). From Proposition 3.1 we know that the solution x 1 (•) reaches x r 1 in finite time. Let t 0 be the first time t for which x 1 (t 0 ) = x r 1 and let x 2,0 = x 2 (t 0 ). We can then consider, without any loss of generality, (x r 1 , x 2,0 ) as initial condition. Let T i be the sequence of times given by Proposition 3.1. The trajectory x(•) is then solution of the non-autonomous dynamics ẋ = F (t, x) with

As O is contractive, ona has the following limit

Appendix 1: list of notations

We remind in the next table all the parameters used in Section 3 and 4.

Appendix 2: Petrov's condition

We recall here a result about the continuity of the first entry time function (see Theorem 8.25 in [START_REF] Cannarsa | Semi-concave Funcions, Hamilton Jacobi Equations and Optimal Control[END_REF]), that is stated here for a non-controlled dynamics. Let g : R n → R n be a mapping of class C 1 with linear growth, and denote by y(•, y 0 ) the unique solution of the Cauchy problem: ẏ = g(y), y(0) = y 0 , defined over R + . Hereafter, we are given a non-empty compact subset T of R n and for y ∈ T , the set N p T (y) stands for the proximal normal cone to the set T at the point y (see [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF]). The standard inner product is written a • b for a, b ∈ R n , and a denotes the euclidean norm of the vector a. Theorem 4.1. Suppose that the Petrov condition

is fulfilled. Then, the first entry time function

If T is convex, the set N p T (y) coincides with the convex normal cone (see, e.g., [START_REF] Clarke | Functional Analysis, Calculus of Variations and Optimal Control[END_REF]) defined for y ∈ T as:

Appendix 3: Asymptotically periodic systems

We recall first the definition of asymptotically periodic semi-flows in R n .

Definition 4.2. A non-autonomous semiflow Φ: {(t, s); 0 ≤ s ≤ t < ∞}×R n → R n is asymptotically periodic with limit ω-periodic semi-flow T (t): R n → R n , t ≥ 0 if Φ(t j + n j ω, n j ω, x j ) → T (t)x, j → ∞ for any sequences t j → t, n j → ∞, x j → x when j → ∞, with x, x j in R n .

The following result can be found in [START_REF] Zhao | Asymptotic behavior for asymptotic periodic semiflows with applications[END_REF] (Theorem 3.1).

Theorem 4.3. Let Φ: {(t, s); 0 ≤ s ≤ t < ∞} × R n → R n be an asymptotically periodic semi-flow with limit ω-periodic semi-flow T (t): R n → R n , t ≥ 0. Denote T n (x) = φ(nω, 0, x) and S(x) = T (ω)x, n ≥ 0, x ∈ R n . If A 0 is a compact subset of R n invariant by the semi-flow S and y ∈ R n is such that d(T n (t), A 0 ) → 0 when n → ∞ then lim t→+∞ d(Φ(t, 0, y), T (t)A 0 ) = 0.