Tom Druet

Arnaud Recoquillay

Bastien Chapuis

Emmanuel Moulin

Passive guided waves tomography for structural health monitoring

In this paper is presented a baseline-free quantitative method for imaging corrosion flaws in thin plates. It only requires an embedded guided wave sensor network used in a fully passive way, i.e. without active emission of waves. This method is called passive guided waves tomography. The aim of this development is the use of this method for the Structural Health Monitoring (SHM) of critical structures with heavy limitations on both sensor's intrusiveness and diagnostic's reliability because it allows to use sensors that cannot emit elastic waves such as Fiber Bragg Gratings (FBGs) which are less intrusive than piezoelectric transducers (PZTs). The idea consists in using passive methods in order to retrieve the impulse response from elastic diffuse fields -naturally present in structures -measured simultaneously between the sensors. In this paper, two passive methods are studied: the ambient noise cross-correlation and the passive inverse filter. Once all the impulse responses between the sensors are retrieved, they are used as input data to perform guided wave tomography.

I. INTRODUCTION

In most industries such as aeronautic, space, railway, shipbuilding, petrochemical or nuclear industry a failure can create a disaster, and thus asserting the structural integrity is critical. Hence, maintenance operations have a great impact on both costs and security.

The Structural Health Monitoring (SHM) consists in acquiring and analyzing data from sensors embedded in a structure in order to evaluate its health (SAE International, 2013). It allows to perform either scheduled maintenance, and thus to anticipate potential failures, or, in the long term, SHM-triggered maintenance. To that end, a SHM system must be able to provide a technical diagnostic with several levels of information on the flaw -detection, localization, classification, estimates of size -and then a technical prognostic in order to estimate the remaining lifetime of the structure.

Various defects needing to be taken care of exists in industrial applications. These defects can appear during the fabrication or during the life of the structure. The SHM generally addresses the issue of the inspection during the operational life of the structure. The flaws can either appear instantaneously or slowly progress over time. This paper deals with the latter kind of flaws such as corrosion or erosion.

A physical way to interrogate large and thin structures such as aircraft fuselages, pipes or ship hulls can be offered by guided elastic waves emitted and received by a network of piezoelectric transducers (PZTs) (Yu and a) tom.druet@cea.fr [START_REF] Yu | In situ 2-D piezoelectric wafer active sensors arrays for guided wave damage detection[END_REF]. Generally, the SHM methods based on elastic guided waves use a baseline. This baseline is compared to the current state by performing, for instance, a subtraction between signals in order to highlight the contribution of the defect. If the algorithms do not take into account external phenomena, these methods with a baseline are usually not robust to changes of environmental conditions. This may cause false alarms, which is very penalizing for SHM systems [START_REF] Croxford | Strategies for guided-wave structural health monitoring[END_REF]. Many methods aim at correcting the effects of the changes of environmental conditions [START_REF] Croxford | Efficient temperature compensation strategies for guided wave structural health monitoring[END_REF][START_REF] Liu | Robust ultrasonic damage detection under complex environmental conditions using singular value decomposition[END_REF][START_REF] Putkis | Continuous baseline growth and monitoring for guided wave shm[END_REF] but, for now, they suffer from performance problems or seem to be impossible to implement in real structures.

In order to avoid the use of a baseline, it has been decided to increase the number of sensors used compared to classical solutions [START_REF] Michaels | Detection, localization and characterization of damage in plates with an in situ array of spatially distributed ultrasonic sensors[END_REF][START_REF] Quaegebeur | Correlationbased imaging technique for fatigue monitoring of riveted lapjoint structure[END_REF][START_REF] Zhao | Active health monitoring of an aircraft wing with embedded piezoelectric sensor/actuator network: I. defect detection, localization, and growth monitoring[END_REF]. This allows to get more information about the structure inspected, and thus to avoid the need of a baseline. Guided wave tomography algorithms are used here to that end. Moreover, this kind of algorithms are able to not only localize but also quantify the severity of defects consisting in a loss of thickness such as corrosion or erosion. However, increasing the number of sensors has a cost, particularly in terms of intrusiveness. A response to that issue could be to use undersampled data [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF][START_REF] Shi | Ultrasonic wave-speed diffraction tomography with undersampled data using virtual transducers[END_REF] but, for now, the number of sensors' reduction achieved is of a factor two to four with respect to the optimal number. This reduction seems to be not enough for some industrial sectors such as aeronautic industry which has a tremendous added mass constrain.

Another interesting possibility is to design a guided wave based SHM system with minimal intrusiveness thanks to the use of fiber Bragg gratings (FBGs). However, these sensors cannot emit waves, they can only work in a passive way. This means that passive methods are needed to retrieve the response between two sensors. Indeed, in the past years researchers have shown that it is possible to use ambient structure-borne elastic noise instead of actively emit the waves in the structure in order to retrieve the impulse response of the medium between two sensors (Sabra et al., 2007a). Moreover, it has been shown in [START_REF] Druet | Passive guided waves measurements using fiber Bragg gratings sensors[END_REF] that passive guided wave measurements are possible thanks to the coupling between FBGs and passive methods. In this paper is shown how to perform passive tomography by coupling guided wave tomography adapted to SHM framework [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF] with passive methods. Finally, the passive methods have several other advantages in the context of SHM [START_REF] Chehami | Detection and localization of a defect in a reverberant plate using acoustic field correlation[END_REF][START_REF] Sabra | Structural health monitoring by extraction of coherent guided waves from diffuse fields[END_REF]. First, since it is not necessary to inject energy in the structure for wave emission, the power consumption of the SHM system can be significantly reduced. This could be of crucial importance for embedded systems with limited or no access to power sources like aircraft or underwater monitoring systems. Then, the complexity of both the electrical cabling and the embedded electronics are reduced, since no emission circuitry is necessary.

This paper shows in a first section the methodology to perform passive guided wave tomography. Two passive methods are presented: the first one is called "ambient noise cross-correlation" and the second one is called "passive inverse filter". Then some keys presented in [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF] to adapt guided waves tomography to the SHM framework are recalled. Finally, comparisons between active and passive signals are presented and passive guided wave tomographies for both passive methods are studied for an ambient noise respecting more or less the equipartition hypothesis.

II. METHODS

In this section is presented how to perform a passive guided wave tomography. First the basis of the two passive methods used in this paper are given and then the main steps of guided wave tomography are recalled.

A. Passive methods: from cross-correlation to passive inverse filter

The first step of any imaging algorithm is the data acquisition. This step usually consists in successively emitting a wave with a sensor and measuring the signal with all others. These signals are called in the following active signals. Another way to obtain equivalent data is to measure the ambient noise at all sensors simultaneously and then process it to reconstruct the data used in the imaging algorithm. This method is called in the following passive acquisition and the corresponding signals are called passive signals.

A few algorithms exist to retrieve the data from the ambient noise. The most known is the ambient noise cross-correlation [START_REF] Lobkis | On the emergence of the Green's function in the correlations of a diffuse field[END_REF], which has been applied in various domains (Campillo and Paul, 2003; Davy et al., 2016; Duvall Jr et al., 1993; Lobkis and Weaver, 2001; Roux et al., 2004; Sabra et al., 2007a , b). This method is fully justified in many cases, but some hypotheses are needed concerning the noise spatial repartition, and it is well known that the quality of the reconstruction is strongly linked to the fulfilment of these hypotheses. To increase the quality of the reconstruction, it is then possible to use the correlation of the coda of the noise cross-correlation [START_REF] Stehly | Reconstructing green's function by correlation of the coda of the correlation (C3) of ambient seismic noise[END_REF]: when many sensors are involved in the acquisition, this method enhance the spatial distribution of the signal by using the coda of all reconstructed signals.

Another possibility is the passive inverse filter [START_REF] Gallot | A passive inverse filter for green's function retrieval[END_REF], which is well suited when many receivers are used as in the case of the guided-wave tomography. As the passive inverse filter is based on the inversion, after regularization, of the noise response matrix, it is a global process, which reconstructs all passive signals from all measured ambient noises simultaneously. In this section will be described the passive inverse filter. To the authors' knowledge, a proof of this method is still an open subject which would give interesting insights on its necessary hypotheses. The description below gives the idea of the method.

In all the following, we consider linear acoustic waves. For a given source f , the acoustic field ϕ in the domain

Ω ⊂ R d , d ∈ {2, 3}, is solution to ∂ 2 ϕ ∂t 2 - 1 c 2 ∆ϕ = f, in Ω × R, (1)
to which is added relevant initial and boundary conditions.

First of all, let us start with a short recall of the ambient noise cross-correlation. Consider a set of point sensors placed at (x j ), j = 0, . . . , N -1. The empirical cross-correlation of the signals u measured over the time interval [0, T ], T > 0, at x i and x j is

C(x i , x j , t; T ) = 1 T T 0 u(x i , τ )u(x j , τ + t) dτ . (2)
Suppose now that the source term f is an unknown noise source, which is a space-time stationary random field that is also delta correlated in space, which means that, for all t, τ > 0 and points x and y, the following relation holds

f (x, t)f (y, τ ) = K(x)δ(x -y)F (t -τ ), (3)
where • denotes the ensemble average, K is a function characterizing the spatial distribution of the noise sources and F is the normalized time correlation function, its Fourier transform being linked to the spectral energy density of the ambient noise. A common further assumption to be able to recover active signals from passive acquisitions using the cross-correlation is then that the ambient noise is well distributed in space, that is K ≡ 1 over a suitable subdomain of Ω, typically all Ω. This last assumption combined with ( 3) is often refered to as the equipartition hypothesis. It can then be proved [START_REF] Garnier | Passive imaging with ambient noise[END_REF]) that lim

T →+∞ dC dt (x i , x j , t; T ) ∝ F * [G(x i , x j , •) -G(x i , x j , -•)] (t), (4)
where * denotes the convolution operation in time and G(•, x j , •) is the Green function of the problem, that is the solution of ( 1) with f = δ(•-x j ). An interesting fact about the cross-correlation is its link to the time reversal operator [START_REF] Gallot | Acoustical imaging in reverberating medias[END_REF]: consider the Fourier transform of the limit of the empirical cross-correlation (2) as T goes to infinity,

C(x i , x j , ω) = u(x i , ω) u(x j , ω). (5)
For a given frequency ω, define the cross-correlation matrix C as

C ij = C(x i , x j , ω). ( 6
)
This matrix is then equal to the time-reversal operator, which would equal to uu * , where u is the signal vector, the coefficients of which are the Fourier transform of the measurements at each sensor. Consider now that P noise acquisitions are done or equivalently if T is big enough, that the acquired ambient noise can be cut into P parts, P > N . Instead of u, consider U which is a N ×P matrix containing all noise acquisitions. The link between the cross-correlation matrix and the time-reversal operator holds:

C = UU * . ( 7
)
This enables the introduction of the passive inverse filter: instead of considering UU * , the passive inverse filter consists in right-multiplying U by a regularized inverse of U, namely Ũ-1 . More precisely, the regularization is done thanks to a singular value decomposition, in which only the biggest singular values are kept, which will be explained at the end of this section. If we denote U = VDW the singular value decomposition of U, where V and W are unitary matrices containing the left and right singular vectors and D is a diagonal matrix containing the singular values of U σ i sorted in descending order, then

U Ũ-1 = V I M 0 0 0 V * , (8)
where M < N is the number of singular values kept after regularization. Again, this operator is close to the timereversal one:

UU * = VD N V * , (9)
where D N is a N × N diagonal matrix of coefficients |σ i | 2 . Hence, the passive inverse filter can be seen as a regularized version of the cross-correlation [START_REF] Seydoux | Preprocessing ambient noise cross-correlations with equalizing the covariance matrix eigenspectrum[END_REF].

The difference between the two operators is the relative weight given to each singular direction. Whereas the time-reversal operator gives a bigger weight to the first singular directions, with a quick decrease due to the squared singular values, the passive inverse filter gives the same weight to all selected singular directions, which seems to enable the extraction of pieces of information from the ambient noise which would be neglected otherwise. In the same way, the normalization at each frequency implies that they all have the same weight in the final reconstructed signal, which induces a whitening of the reconstructed signal. The difficulty is then to correctly choose the singular directions of interest, which will be kept in the regularized inverse of U. To do so, a classic technique is to consider the variations of the modulus of the singular values: a usual behavior consists in a fast decay followed by a slower decay. The point of the change of slope is a good truncation candidate [START_REF] Gallot | A passive inverse filter for green's function retrieval[END_REF] and will be used in our experimental study.

B. Guided wave tomography

In this article, the used guided wave tomography algorithm is based on the original iterative Hybrid Algorithm for Robust Breast Ultrasound Tomography (HAR-BUT) (Huthwaite and Simonetti, 2013) but adapted to the heavy constraints imposed by the SHM framework. This adaptation has been detailed in a previous work [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF]. The same algorithm is used here to produce images of corrosion-like defects but this time using passive data. In this section are recalled the main elements of the imaging algorithm and a regularization method called "variable relaxation" which help to improve significantly the quality of the image when undersampled data are used.

In all the following, a single mode is considered. The wavenumbers and velocities used in this section are implicitly those of this mode. The iterative version of HAR-BUT includes the following principal steps. In order to initialize the algorithm, a time-of-flight tomography image of the plate thickness is computed in order to fulfill the Born approximation in the subsequent diffraction tomography step. For the time-of-flight tomography, the structure to be imaged is represented at the position x by the scalar slowness field s(x) = 1 /vG(x), where v G is the group velocity of the considered mode. To represent the slowness s(x), we project it on a finite basis. The projection of the scalar slowness s(x) thus reads

s(x) = i,j S ij β ij (x), ( 10
)
where β is a 2D B-spline basis and S ij is the coefficient associated to the B-spline β ij . The slowness is then linked to the time-of-flight τ k→l for a trajectory γ k→l between the emitter k and the receiver l:

τ k→l = i,j S ij A ijkl , (11)
where A ijkl = γ k→l β ij (γ) dγ is the contribution of the basis element β ij to the time-of-flight τ k→l . The image is then updated thanks to an absolute estimator based on the difference between the data τ defect k→l and the computed time-of-flight τ (n-1) k→l on the (n -1) th map. The n th iteration of the time-of-flight tomography step reads:

S (n) ij = S (n-1) ij + δ k S (n) ij , where δ k S (n)
ij is the correction applied to the coefficient S ij and is computed according to the following formulation:

δ k S (n) ij = 1 N ijk l W ijkl τ defect k→l -τ (n-1) k→l λ k→l , (12)
where λ k→l is the trajectory's length, W ijkl is a weighting of the path, consisting of A ijkl in the simplest case but usually including a Hamming window to avoid side effects, and the normalization factor N ijk reads:

N ijk = l A ijkl .
Once the time-of-flight tomography image is computed, it becomes the first guess of an iterative diffraction tomography process. The image is then represented thanks to an object function O(x), which reads

O(x) = k(x) 2 -k 2 0 = k 2 0 v 0 v(x) 2 -1 , ( 13
)
where k 0 is the background wavenumber, v 0 the background phase velocity and v(x) the phase velocity at the position x, which is linked to the thickness thanks to the dispersion curves when the algorithm has converged. One can show (Huthwaite, 2016a) that O(x) reads

O(x) = π -π π -π ϕ s k0 (x l ) G 0 (x l ; x)G 0 (x; x k ) W(θ k , θ l ) dθ k dθ l , (14)
where ϕ s is the scattered field, G 0 the Green's function of the unperturbed domain, W a weighting function coming from a variable change, x k , respectively x l , the position of the emitter k, respectively the receiver l, and θ i , i ∈ {k, l}, the angles of the emitter and receiver. Note that it is possible (Huthwaite, 2016a) to replace scattered field ϕ s by the total field ϕ. Iterative HARBUT consists in computing iteratively the correction O δ (x), which is added to the previous estimation O

b (x) so that O(x) = O b (x) + O δ (x).
More details about original HARBUT can be found in (Huthwaite and Simonetti, 2013) but also in (Huthwaite, 2014 , 2016a , b , c; Huthwaite et al., 2013) for many following developments. Concerning the adaptations to the constraints imposed by the SHM framework which are used also in this paper, they can be found in [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF]. We recall here one of these adaptations which concern the regularization called "variable relaxation" and which allows to limit the amplification of spatial fluctuations at each iteration due to the noise and the lack of sensors, while still allowing us to accurately reconstruct localized defects above the noise threshold. Let e ij and ẽ be two parameters depending on the used algorithm. More precisely, let e ij be the difference between the computed thickness and the pristine one denoted ẽ for the time-of-flight tomography, and let e ij be the object function O and ẽ = -min O for the diffraction tomography. Let e ij be the new considered value at the position ij. The "variable relaxation" is defined as follows:

e ij =      e ij 1 + 1 z 2 ij α 2 if |e ij | < γβẽ, e ij otherwise, (15)
with β a threshold depending on the amplitude fluctuation of the noise level of the image, γ a constant determining the relaxation limit, α a constant defining the strength of the regularization and z ij reads

z ij =          |e ij | βẽ if |e ij | < βẽ, |e ij | 1 2 βẽ 1 -cos π |eij |-γβẽ (1-γ)βẽ if βẽ ≤ |e ij | < γβẽ. (16)
As explained in [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF], we chose to take γβ = 1, β = 0.1, α = 4 for the time-of-flight tomography and α = 2 for the diffraction tomography.

III. RESULTS

In this section, we present experimental results on passive tomography. First, the experimental setup is presented. Then the differences between the crosscorrelation and the passive inverse filter presented in the section II A are illustrated with classical signal processing tools. In order to confirm the good convergence of passive methods, passive and active signals are compared for several distances in the following section. Finally, we present passive guided wave tomography images for both cross-correlation signals and passive inverse filter signals.

A. Experimental setup

The experimental setup used to obtain all the results of this paper is shown in figure 1. We consider an aluminum plate of density ρ = 2.7 g.cm -3 , longitudinal velocity V L = 6360 m.s -1 , transversal velocity V T = 3140 m.s -1 and dimensions 1400 mm × 1250 mm × 2 mm. On this aluminum plate are glued 30 piezoelectric transducers of diameter ∅ = 18 mm using epoxy glue. These sensors are equally distributed on a circle of radius r = 300 mm.

The ambient noise required to perform passive tomography has been generated in the plate by manually spraying compressed air randomly on its surface [START_REF] Druet | Passive guided waves measurements using fiber Bragg gratings sensors[END_REF][START_REF] Larose | Reconstruction of Rayleigh-Lamb dispersion spectrum based on noise obtained from an air-jet forcing[END_REF]. The distance between the nozzle and the plate was of approximately 5 cm. Ten seconds sampled at 2.5 MHz of this ambient noise were recorded and then post-processed in order to obtain passive signals for all the couples of sensors.

B. Cross-correlation versus Passive inverse filter

In this section are illustrated the differences between the cross-correlation and the passive inverse filter on typical passive signals.

First, we present raw passive signals in figure 2 for a distance between sensors of 52 cm. The case of the crosscorrelation is presented in figure 2(a) whereas the case of the passive inverse filter is presented in 2(b). We also present in figure 3 there spectral content.

Concerning the cross-correlation, we can identify wave packets around 400 µs for the causal part and -400 µs for the anticausal one. They correspond approximately to the theoretical time-of-flight of the A 0 mode for low frequencies, around 20 kHz -30 kHz. This seems to be consistent because the spectral content of the crosscorrelation is centered around 25 kHz as can be seen figure 3. Moreover, the cross-correlation of figure 2(a) is globally symmetrical, which is expected from the passive reconstruction and is an indication of a good reconstruction.

Concerning, the passive inverse filter, we would like to recall that the 30 sensors where used altogether to obtain the passive signal between the two sensors of figure 2(b) as explained in section II A. It is not possible to identify specific wave packets on the passive signal of figure 2(b) because of the whitening effect of the passive inverse filter: indeed, as can be seen figure 3, the spectral density of the passive inverse filter has a wider frequency range and the spectral density is somewhat constant over this frequency range, whereas the signal reconstructed thanks to the cross-correlation has the same energy repartition as the elastic noise. Thus high-frequency components, including electromagnetic noise, will contribute to the passive signal resulting from the inverse filter process. This noise induces significant autocorrelation peak at zero time as can be observed in figure 2(b). Since this peak, as well as the whole part of the signal be- fore the first modal wavepacket arrival (at approximately 100 µs), does not contain useful physical information, it is set to zero between -65 µs and 65 µs. This operation eliminates undesirable electromagnetic noise, while keeping the information concerning the guided modes. The figure 3 shows that the electromagnetic coupling has its energy at frequencies higher than 150 kHz / 200 kHz.

In order to compare the physical information of the passive signals, the signals are projected in a time-frequency plan thanks to a reassigned spectrogram [START_REF] Auger | Improving the readability of time-frequency and time-scale representations by the reassignment method[END_REF] presented in figure 4. The same experiment as before has been conducted but with a larger bandwidth in order to highlight the whitening effect of the passive inverse filter. To this end, we used two piezoelectric transducers of diameter 13 mm spaced 60 cm apart from each other. The reassigned spectrogram corresponding to the cross-correlation signal is plotted figure 4(a) and the one corresponding to the passive inverse filter is plotted figure 4(b). Once again it can be seen that the passive inverse filter has a whitening effect compared to the cross-correlation. Moreover, the power spectral density in both cases corresponds to the theoretical time-of-flight (red dash-line) of the A 0 mode.

In both cases, the reconstructed signals are coherent with the theoretical ones. Signals reconstructed thanks to the passive inverse filter have a wider frequency range, which may be of interest for applications. It should be noted that the influence of the spatial repartition of the noise could not be fully studied in this case as the used plate was very reverberant, inducing a good spatial repartition of the ambient noise even for a localized noise source.

C. Passive and active signals comparisons

In this section are presented comparisons between active and passive signals for the frequency range of interest of this paper and for both the cross-correlation and the passive inverse filter.

Active signals s(A, B, t) are obtained by emission of a 5 cycles toneburst at 20.4 kHz from a piezoelectric transducer placed in A. The receiver is placed in B. The elastic guided waves are measured after their propagation for several distances between the emitter and the receiver (see figure 5). In order to get the same spectral density for passive signals named s P (A, B, t) and active signals s(A, B, t), the following normalization is performed [START_REF] Druet | Passive guided waves measurements using fiber Bragg gratings sensors[END_REF]:

s P (A, B, ω) =        C(A, B, ω) C(A, B, ω) | s(A, B, ω)| ω ∈ Ω ε , C(A, B, ω) otherwise, (17)
where • denotes the Fourier transform of a variable and the set Ω For both cases, a very good match is observed. The theoretical time-of-flight is also plotted and corresponds to the maximum of the envelop of the wave packet corresponding to the A 0 mode. Note that in both cases (figure 5(a) and 5(b)) the electromagnetic coupling around the time t = 0 µs is observed.

ε reads Ω ε = ω ∈ R | C(A, B, ω)| > ε max ω | C(A, B, ω)| ,
The study in this section allows us to conclude that the passive signals seems to be of adequate quality to perform passive imaging and more precisely passive guided wave tomography. This is the subject of the following section.

D. Passive images with guided wave tomography

In order to test passive tomography a defect was created into the zone of interest (see the figure 6) by machining locally the aluminum plate. A reference of the defect has been obtained thanks to an immersion ultrasonic scan, such a scan being significantly more precise than a global method such as guided wave tomography. The reference is plotted on figure 7.

During the acquisition of the ambient noise, an effort has been made to respect, as a first step, the equipartition hypothesis as much as possible. This means that the compressed air jet was moved on the whole surface in order to spend an equivalent time on each surface areas. Passive tomography is first computed without the use of the "variable relaxation". The result is visible fig- ure 8(a) for the cross-correlation and figure 8(b) for the passive inverse filter. We can see on both images typical aliasing due to the fact that we use undersampled data (see [START_REF] Druet | Autocalibration method for guided wave tomography with undersampled data[END_REF] for more details). Moreover, both passive tomography images seem to be of the same quality.

In order to compensate the aliasing effect, the "variable relaxation" is added in the guided wave tomography process, as explained section II B. The results are presented figure 9, image (a) for the cross-correlation and image for the passive inverse filter. Once again, both images are of same quality but this time the aliasing effect is significantly reduced. The background of the image is very clean so that the defect clearly appears in the inspected area. Moreover, compared to the reference, the localization of the defect is well identified and the geometry is almost fully reconstructed. In order to compare more quantitatively the results, we propose cross-sections of the images on figure 10. The position of the crosssection is visible on the reference (figure 7). The small variations within the defect visible in the reference of the figure 7 are the only thing that is not reconstructed by the passive tomography algorithm. This is because the working wavelength is too large to be able to reconstruct such a variation.

Passive tomography gives really good results with an ambient noise which respects as much as possible the equipartition condition but this condition is not always respected in real industrial applications. That is why we now study a case which is closer to reality. Indeed, this time the compressed air jet is only moved onto left-half of the plate area. This is supposed to throw the energy repartition off balance. The passive guided wave tomographies corresponding to that configuration are presented figure 11(a) for the cross-correlation and figure 11(b) for the passive inverse filter. The results are of the same quality as the ones for the equipartition case (figure 9). This can be explained by the fact that the aluminum plate is a very reverberant structure compensating for the bad equipartition of the ambient noise.

To conclude, passive guided wave tomography gives a really good image of corrosion flaws whether it be with the cross-correlation or with the passive inverse filter. However, as explained in section II A, the passive inverse filter should give better results than the cross-correlation in cases where the noise does not fill the equipartition hypothesis. The structure studied in this paper did not permit to break this hypothesis and to illustrate this fact. Future studies in this direction seem interesting.

IV. CONCLUSIONS

In this paper were presented imaging results of corrosion-like flaws using a passive guided wave tomography algorithm. To do so, two passive data reconstruction were compared: the usual ambient noise cross-correlation and the passive inverse filter. Even though the results are very similar for both methods, the passive inverse filter seem to be a promising method as its working hypotheses seem to be weaker. Indeed, the passive inverse filter can be seen as a regularized version of the cross-correlation and thus, it should be more robust. A thorough study is yet to be done. Furthermore, the passive inverse filter is naturally well-suited for configurations using many sensors, which is the case of tomography algorithms.

This study was done using classical piezoelectric transducers, but the use of passive algorithms allow the use of sensors not able to emit waves, for example the fiber Bragg gratings. These sensors are really interesting because of their low intrusiveness in the context of SHM and our current work focus on the use of fiber Bragg gratings.

The application presented in this paper is a rather well-controlled laboratoy configuration. We plan on applying the method to real structures, using the ambient noise coming from their environment to image it.

FIG. 1

 1 FIG. 1. (Color online) Aluminum plate and piezoelectric transducers positions.

  FIG. 2. Raw signals of (a) cross-correlation and (b) passive inverse filter.

  FIG. 4. (Color online) Reassigned spectrogram of (a) crosscorrelation raw signal and (b) Passive inverse filter raw signal.

  FIG. 5. (Color online) Comparison between active signals and (a) cross-correlation (b) passive inverse filter.

FIG

  FIG. 6. (Color online) Defect position.

  FIG. 8. (Color online) Passive guided wave tomography without regularization for (a) the cross-correlation and (b) the inverse filter.

  FIG. 9. (Color online) Passive guided wave tomography with regularization for (a) the cross-correlation and (b) the passive inverse filter.

  FIG. 10. (Color online) Cross-sections of the passive guided wave tomographies with regularization.

  FIG. 11. (Color online) Passive guided wave tomography with regularization for (a) Correlation (b) Passive inverse filter. Non equipartition case.

J. Acoust. Soc. Am. / 23 October 2020

J. Acoust. Soc. Am. / 23 October 2020


Guillaume Fraysse
email: guillaume.fraysse@orange.com

Jonathan Lejeune

Julien Sopena

Pierre Sens

A resource usage efficient distributed allocation algorithm for 5G Service Function Chains

Keywords: Computer network management, Distributed algorithms, Network Slicing, Distributed Systems, k-mutex, drinking philosophers, deadlock

Recent evolution of networks introduce new challenges for the allocation of resources. Slicing in 5G networks allows multiple users to share a common infrastructure and the chaining of Network Function (NF)s introduces constraints on the order in which NFs are allocated. We first model the allocation of resources for Chains of NFs in 5G Slices. Then we introduce a distributed mutual exclusion algorithm to address the problem of the allocation of resources. We show with selected metrics that choosing an order of allocation of the resources that differs from the order in which resources are used can give better performances. We then show experimental results where we improve the usage rate of resources by more than 20% compared to the baseline algorithm in some cases. The experiments run on our own simulator based on SimGrid.

Introduction

The flexibility of 5G networks allows the apparition of new services. Complex services rely on Slices split across multiple Network Service Provider (NSP)s. The allocation of a service is now not only the allocation of a single NF but the chaining of multiple NFs. Chains of NFs introduce a constraint on the order of the allocation.

We argue that in such use cases the system can be modeled as a distributed system. The various NFs from the different NSPs can be abstracted as resources. One or more instances of these resources can be available for each of these resources. We define the allocation of a chain of NFs as a request of a set of resources with an associated request order. The allocation of a request means allocating all the resources in the set while respecting this request order.

In distributed systems, allocation of resources is seen as a Mutual Exclusion problem. Several variants of the problem have been defined considering a single resource [START_REF] Lamport | Time, Clocks, and the Ordering of Events in a Distributed System[END_REF][START_REF] Ricart | An Optimal Algorithm for Mutual Exclusion in Computer Networks[END_REF] or multiple instances of a single type of resource [START_REF] Raymond | A Tree-Based Algorithm for Distributed Mutual Exclusion[END_REF][START_REF] Raynal | A Distributed Solution to the K-out of-M Resources Allocation Problem[END_REF]. Very few solutions have been proposed for the problem of Mutual Exclusion for systems with multiple instances of multiple resources [START_REF] Bouabdallah | A Distributed Token-Based Algorithm for the Dynamic Resource Allocation Problem[END_REF][START_REF] Lejeune | Reducing Synchronization Cost in Distributed Multi-Resource Allocation Problem[END_REF].

We propose here a distributed algorithm to solve the problem of the allocation of resources to create chains of NFs in 5G Slices for Multi-Domain use cases. We address this as a distributed Mutual Exclusion problem. We show that for systems with multiple instances of resources the selection of the instances to allocate has an influence on performances. Our algorithm extends the LASS algorithm [START_REF] Lejeune | Reducing Synchronization Cost in Distributed Multi-Resource Allocation Problem[END_REF] for systems with N instances of M resources. We introduce a subroutine to select the instance of the resource based on the orders of the requests as the LASS algorithm does not address this constraint. The algorithm is based on the transmission of a token that contains the permissions to use the resources. In a Network Functions Virtualization (NFV) network the resources are the nodes themselves, and they can't be transferred from one node to another. In such system each node is the manager of its own resource. We propose an extension to manage systems where the decisions to allocate the resources are made locally by each node.

We introduce a broader description of these use cases as well as some related work and background in Section 2. We describe our problem in Section 3. In Section 4 we introduce our algorithm, and define the allocation order as distinct from the order in the requests. We then introduce the methodology used to evaluate the algorithm with the SimGrid [START_REF] Casanova | Versatile, Scalable, and Accurate Simulation of Distributed Applications and Platforms[END_REF] simulator and show the experimental results in Section 5. We finally present our conclusions and future works in Section 6.

Related work

The architecture of Telecom Networks is rapidly evolving. Operators have launched the NFV [START_REF][END_REF] initiative at the European Telecommunications Standards Institute (ETSI) to move the NFs from dedicated hardware to virtualized infrastructures based on Cloud and Software-Defined Networking (SDN) technologies.

Allocation of a single NF is rarely sufficient. More complex services require multiple NFs to inter-operate while respecting an order. To this end ETSI defined the VNF Forwarding Graph (VNF FG)s [START_REF]ETSI GS NFV 001: Network Functions Virtualisation (NFV) Use Cases[END_REF] as the analogue to the connection of physical appliances with cables. Following the description of the VNF FG use case in 2013, RFC7665 [START_REF] Pignataro | Service Function Chaining (SFC) Architecture[END_REF] introduced in 2015 Service Function Chaining (SFC) to allow the allocation of a chain of NFs. In 2017, RFC8402 [START_REF] Previdi | Segment Routing Architecture[END_REF] introduced Segment Routing to allow a node to send a list of instructions to be run by subsequent nodes. problem of resources allocation for multiple NFs in the correct order.

5G is the first generation of mobile networks to include NFV as an enabler for new types of services. 3GPP has introduced Network Slices [29,[START_REF] Rost | Network Slicing to Enable Scalability and Flexibility in 5G Mobile Networks[END_REF] that enables multiple virtual end-to-end networks to share the same NFV infrastructures. A service offered by a Slice can rely on the infrastructures of multiple NSPs, it is then called Multi-Domain. This can be the case when a large op-erator has split its network in multiple subdomains (e.g. datacenters) or when a use case may require the infrastructures of multiple operators.The European 5G Public Private Partnership (5G PPP) initiative launched projects that defined several use cases based on Slices and Multi-Domain. The SLICENET project's eHealth use case [29] requires multiple NSPs to provide Slices that are chained together to provide the service.

Multiple centralized solutions exist for the allocation of resources such as Virtual Network Function (VNF) in networks or the placement of Virtual Machines in Cloud infrastructure [START_REF] Mills | Comparing VM-Placement Algorithms for On-Demand Clouds[END_REF][START_REF] Rai | Generalized Resource Allocation for the Cloud[END_REF][START_REF] Jennings | Resource Management in Clouds: Survey and Research Challenges[END_REF][START_REF] Widjajarto | Cloud Computing Reference Model: The Modelling of Service Availability Based on Application Profile and Resource Allocation[END_REF]. These approaches focus on finding an optimal placement depending on a set of constraints. Some papers focus on finding heuristics to respect Service Level Agreement (SLA)s [START_REF] Machida | Redundant Virtual Machine Placement for Fault-Tolerant Consolidated Server Clusters[END_REF] or security rules [START_REF] Jhawar | Supporting Security Requirements for Resource Management in Cloud Computing[END_REF]. These problems are mostly addressed with Integer Linear Programming (ILP) formulations. This centralized method may not be adequate for Multi-Domain use cases when it is not be possible to have a centralized manager or when the cost of building a global view of the system has a high cost. A centralized method also often requires an a priori knowledge of the requests.

We propose to address such systems as distributed systems, and propose a solution where there is no centralized manager. In such systems resources all execute the algorithms locally and get (resp. send) information from (resp. to) other resources by the passing of messages.

The allocation of resources in distributed systems can be handled as a Mutual Exclusion problem on these resources. The Mutual Exclusion is a fundamental problem in distributed systems and was first described by E. W. Dijkstra in 1965 [START_REF] Dijkstra | Solution of a Problem in Concurrent Programming Control[END_REF] for systems where multiple processes try to allocate concurrently a single shared resource. Allocating this resource allows them to execute a portion of code known as Critical Section (CS) allowing processes to use the resource exclusively. Multiple solutions [START_REF] Lamport | Time, Clocks, and the Ordering of Events in a Distributed System[END_REF][START_REF] Ricart | An Optimal Algorithm for Mutual Exclusion in Computer Networks[END_REF][START_REF] Suzuki | A Distributed Mutual Exclusion Algorithm[END_REF][START_REF] Naimi | A Log (N) Distributed Mutual Exclusion Algorithm Based on Path Reversal[END_REF] have been proposed.

The mutual exclusion problem was later generalized in two ways: -for systems with one instance of M resources known as the dining philosophers problem, when the requests are static, and drinking philosophers problem, when the requests are dynamic. It was defined by K. M. Chandy and J. Misra in 1984 [START_REF] Chandy | The Drinking Philosophers Problem[END_REF]. -for systems with k instances of a single resource, known as the k-mutex problem [START_REF] Raymond | A Tree-Based Algorithm for Distributed Mutual Exclusion[END_REF]. A variant of this problem is known as the k-out of-M resources allocation problem [START_REF] Raynal | A Distributed Solution to the K-out of-M Resources Allocation Problem[END_REF] when one process tries to allocate multiple instances of a single type of resource. Algorithms to solve drinking philosophers problems, need to address potential conflicts between two requests. A conflict occurs when two requests try to allocate a common resource. If two requests don't conflict, they are allowed to enter their CS simultaneously.

The Dining/Drinking philosophers problem was generalized in 1990 by Awerbuch and Saks [1] as the Dynamic Job Scheduling problem where processes require resources that can be used by a single process at a time. A job can only be executed when all its required resources are available for exclusive use by the process. The problem is also related to the job-shop scheduling optimiza-tion problem, in which n jobs have to be scheduled on m machines with the shortest-length schedule. This problem is NP-hard [START_REF] Garey | The Complexity of Flowshop and Jobshop Scheduling -Mathematics of Operations Research[END_REF].

Algorithms addressing the mutual exclusion for systems with M resources can be divided into two groups: incremental algorithms and simultaneous algorithms. Algorithms in the first group incrementally allocate resources according to a static total order on the resources. They are using consecutive mutexes on each of the M resources. E. W. Dijkstra's algorithm from this group [START_REF] Dijkstra | Hierarchical Ordering of Sequential Processes[END_REF] is the baseline algorithm for our comparison and is detailed in Section 5.2. Algorithms in the second group do not set a predefined total order on the resources but try to simultaneously allocate resources for multiple requests. To achieve this multiple mechanisms have been proposed. Some require a knowledge of the conflict graph [START_REF] Chandy | The Drinking Philosophers Problem[END_REF][START_REF] Ginat | An Efficient Solution to the Drinking Philosophers Problem and Its Extensions[END_REF]. Others rely on a broadcast mechanism with high messages complexity [START_REF] Maddi | Token Based Solutions to M Resources Allocation Problem[END_REF] or a global lock that is not efficient when the concurrency of requests is high [START_REF] Bouabdallah | A Distributed Token-Based Algorithm for the Dynamic Resource Allocation Problem[END_REF]. All the simultaneous algorithms have in common to build a total order of the requests to schedule them.

Finally, it is possible to extend drinking philosopher and k-mutex problems by considering systems with N instances of M types of resources and requests for k instances of 1 or more types, we call it the k -M -N problem.

In a system with N instances of M resources, it is necessary to decide which instance to allocate for a given request. Once an instance has been selected, we have a simplification of the k -M -N problem to the drinking philosophers problem, where each instance is uniquely identified by its location. To the best of our knowledge this problem has not been specifically addressed. Some papers address the drinking philosophers problems and mentioned possible extension of their algorithms to systems with N instances [START_REF] Bouabdallah | A Distributed Token-Based Algorithm for the Dynamic Resource Allocation Problem[END_REF][START_REF] Lejeune | Reducing Synchronization Cost in Distributed Multi-Resource Allocation Problem[END_REF] but did not consider the selection of the instances as a specific constraint.

Algorithms from the state of the art don't consider the latency of the network links. They also do not address that the nodes selected for a chain of NFs need to respect a specific request order. In our model we add a weight to the edges of the graph to take this latency into consideration and be able to compute a path that respects the order in which resources are used. They also do not take into consideration that network links are not First In First Out (FIFO) channels. Our algorithm makes no assumption on the order in which messages are received.

The LASS algorithm [START_REF] Lejeune | Reducing Synchronization Cost in Distributed Multi-Resource Allocation Problem[END_REF] is a simultaneous algorithm that addresses systems with a single instance of M resources. It has been shown that its performance are better than those of incremental algorithms. It builds allocation vectors for all requests. These vectors are then used to compute a total order on requests, as detailed in section 3. Our algorithm extends it and includes a preemption mechanism that is used when messages are received in an order that is different from the total order of the requests.

N instances of M resources. In an example of VNF FG described in [START_REF]ETSI GS NFV 001: Network Functions Virtualisation (NFV) Use Cases[END_REF], packets need to traverse an Intrusion Detection Systems (IDS), a Firewall (FW) and a Load-Balancer (LB). The left part of Figure 1 shows this example with three NSPs.The figure shows the cases where there are three instances of each resource distributed across the three NSPs. We model our system as a non-directed connected static graph G=(N , E) where N is the set of nodes and E the set of edges. A node contains at most one resource. Edges have positive weights that allow to model the latency of the links between nodes. A weight of 0 on an edge allows to model a system with multiple resources on a single node. A node with 2 resources can be modeled in the graph as two nodes holding one resource each and connected by a zero-weight edge. We note C the set of types of resources.

Each node in the graph can issue an allocation request. A request is modeled as a couple Req(n, {c 1 , . . ., c s }) where :

n is the requesting node, -{c 1 , . . . , c s } where c r ∈ C, ∀r is an ordered set of types of resources needed. The request order gives the order of the resources in the request. The order of resources can be different across requests. The right side of Figure 1 shows a model of the use case introduced above. We pose a request Req

1 = Req(n 1 , {c 3 , c 2 , c 1 }) in this system.
n 1 is the requesting node, 3 types of resources c 1 , c 2 and c 3 are requested. The request order is c 3 ≤ c 2 ≤ c 1 , i.e. first c 3 , then c 2 and finally c 1 .

Our algorithm does not require a knowledge of the conflict graph like [START_REF] Chandy | The Drinking Philosophers Problem[END_REF][START_REF] Ginat | An Efficient Solution to the Drinking Philosophers Problem and Its Extensions[END_REF]. It requires that each node has knowledge of its neighbors and knows where to find each type of resource so that each node can send messages to others.

The first subroutine of the algorithm presented in Section 4 computes a path in the graph. A path contains all the nodes from, and including, the requesting node, to the last resource requested, respecting the request order. The path contains the nodes holding the requested resources as well as the nodes connecting those, e.g. a valid path for Req 1 is the ordered set of nodes (n 1 , n 5 , n 7 , n 8 , n 7 , n 5 ). The originating node is n 1 , the first type of resource requested is c 3 and n 8 holds an instance of it. It is necessary to go through n 5 and n 7 to reach n 8 from n 1 . Then n 7 and n 5 hold the two other requested resources.

A request is satisfied when the path contains nodes that hold all the types of resources in the correct request order and all the requested resources are allocated to it, allowing the requesting node to enter its CS.

The algorithm builds an allocation vector for each request. Each entry of an allocation vector is a pair (n, value), called a counter, where n ∈ N is the node having the requested resource and value is a positive integer incremented by the node upon reception of the request, as detailed in Section 4. It is not possible for two requests to get the same counter value for a node n in their allocation vectors. The allocation vector for request Req r is noted V Reqr = ((n r 1 , counter r n1 ), . . . , (n r s , counter r ns )) where r is the identifier of the request and s is the size of request Req r , e.g. a valid allocation vector for Req 1 is V Req1 = ((n 8 , 3), (n 7 , 2), (n 5 , 4)). The allocation vectors allow the algorithm to compute a total order of the requests.

To sort the requests according to the global order we define the precedence of a request Req as its rank in the total order. Req i precedes Req j if the average of the counters of the allocation vector V Reqi is less than the average of the counters of the allocation vector V Reqj . For instance, if we consider

Req 2 = Req(n 4 , {c 3 , c 2 , c 1 }) and Req 3 = Req(n 6 , {c 1 , c 2 , c 3 }) with V Req2 = ((n 6 , 3), (n 3 , 3), (n 4 , 1)
) and V Req3 = ((n 4 , 2), (n 3 , 2), (n 6 , 1)) their respective allocation vectors, Req 3 will be allocated before Req 2 since 3+3+1 3

> 2+2+1 3 .
If two allocation vectors have the same average value, it is necessary to add a way to break the tie. No generic method is proposed here since it can be implementation-specific, but the method used in our experimental platform is detailed in section 5.3.

Algorithms

Our algorithm consists of two consecutive subroutines:

the path computation subroutine, in which the algorithm selects the resources instances to be allocated and computes a routing path between them that respects the allocation order present in the request, -the allocation subroutine in which the algorithm allocates the resources selected during the path computation subroutine.

Path computation subroutine

This subroutine assumes that each node of the system has some local knowledge on how to reach each type of resource. Each node keeps an up-to-date local routing table containing the name of the neighbor node that is the closest to each type of resource as well as the distance to this node. The entry in the routing table can be the node itself for the type of resource it is holding. How these routing tables are built is out of scope of this paper. The solution used in our simulation is described in Section 5. Table 1 gives some routing tables for the system in Figure 1.

Type Node D c1 n2 1 c2 n1 0 c3 n5 9 (a) n1 Type Node D c1 n5 0 c2 n1 1 c3 n8 8 (b) n5 Type Node D c1 n5 7 c2 n7 0 c3 n7 1
(c) n7 Table 1. Routing tables for n1, n5 and n7. D is for Distance.

In this example, node n 1 's shortest path to the type of resource c 3 is of length 9 and starts at n 5 : (n 1 , n 5 , n 7 , n 8 ), and node n 5 's shortest path to the type of resource c 3 is of length 8 and starts at n 7 : (n 5 , n 7 , n 8 ).

n1 n5 n7 n8 t R O U T I N G R O U T I N G R O U T IN G R O U T I N G { ( n 8 , 3 ) } R O U T I N G { ( n 8 , 3 ) , ( n
7 , 2 ) } A L L O C A L L O C A LL O C A C K E N D C S E N D C S E N D C S

Fig. 2. Algorithm execution for Req1

Path computation relies on the ROUTING message. A first ROUTING message is sent by the requesting node to the node holding the first resource requested, according to its routing table. This node then sends a ROUTING message to the next node on the path to the node holding the next resource. The operation is repeated on each node until a complete path with all the requested resources in the correct order has been found. Figure 2 shows the messages sent during the execution of the algorithm in the sample system.

This subroutine computes a valid path with the selected instances of all the requested resources in the order given in the request. It does not check that the resources are available nor does it guarantee that the path computed is the shortest path containing all the resources.

Once the algorithm has reached the node holding the last resource requested, the allocation subroutine described below starts from the last node of the path.

Building the allocation vector We compute a total order of the requests to preserve the liveness property, i.e. it guarantees that all requests are satisfied in a finite time. We use the method from the LASS algorithm to compute this order based on vectors build for each request.

The first step is to build the allocation vector for the request. The vector is built by each of the node on the computed path. Each node has a local counter that is initialized to 0. This counter acts as a logical clock, but contrary to Lamport's logical clocks [START_REF] Lamport | Time, Clocks, and the Ordering of Events in a Distributed System[END_REF] it is local and only incremented by the node when it receives a new request. As such it is then not possible for two requests to get the same counter value for a node in their allocation vector. The first request receives the value 1, the second receives the value 2, and so on. Upon reception of a request a node increments its local counter. It then updates the allocation vector with the value of its counter. The updated allocation vector is inserted in the request message when it is forwarded to the next node in the path. As a simplification this has been merged in our implementation with the ROUTING messages. Thus, the allocation vector is built during the forwarding of the ROUTING messages as shown in Figure 2.

The allocation subroutine

The allocation subroutine allocates the instances of resources selected during the path computation subroutine. Nodes receive allocations messages in an arbitrary order, and it can be completely different from the order computed on the requests. To deal with this situation, the subroutine has a preemption mechanism to enforce the total order of the requests. Simulations show that the allocation order within a request has a strong impact on the performance of the algorithm, we compare here two allocation orders. All this is detailed below.

The allocation The core of the allocation subroutine is based on the ALLOC messages. In a system where all the resources are initially in the IDLE state, this message is sent to the first node in the allocation order, detailed below. This node enters the ALLOCATED state and sends an ALLOC message to the next node. The operation is then repeated until the last node in the allocation order is reached. Then the last node send an ALLOC ACK message to the requesting node to inform it that the allocation of all resources has been made. It then enters its CS and starts using the resources. Upon leaving its CS it sends a END CS message that is forwarded along the path to all the nodes holding the requested resources. The messages sent for the allocation of Req 1 are shown on Figure 2.

The allocation order Each request has an associated partial request order for the resources within the request, i.e. the order in which the resources are used, cf. Section 3. We define the allocation order as the order used by the algorithm to allocate the resources. There is no relation between the request order and the allocation order and they can be different for a same request. As a mean of comparison, we introduce two different allocation orders that are then evaluated in section 5.

The by values allocation order sends the ALLOC messages according to the values of the counters in the allocation vector. The last node of the path sends a ROUTING ACK message to the requesting node. Upon reception of this message the requesting node sends an ALLOC message to the node with the highest counter value in the allocation vector. The allocation then follows the order of the counters in the allocation vector. As an allocation based on this order starts by allocating the node with the highest counter value, it reduces the probability that requests with higher precedence will arrive during the rest of the allocation of the request.

The reverse order allocation order sends the ALLOC messages in the reverse order of the routing path. The first subroutine follows the path as it selects it. With this reverse order, the allocation follows the path backwards to go back to the requesting node.

For Req 1 with the computed path (n 1 , n 5 , n 7 , n 8 , n 7 , n 5 ) if the allocation vector is V Req1 = ((n 8 , 3), (n 7 , 2), (n 5 , 4)), the allocation for reverse order follows the order n 5 , n 7 , n 8 . For by values the order is n 5 , n 8 ,n 7 .

Preemption of resources Since the system is distributed a request can arrive on a node already in the ALLOCATED state for a request that has a lower precedence. This can lead to deadlocks. To manage these situations the algorithm preempts the resources to enforce the global order on the requests. This requires an additional state for the nodes, PREEMPTING, and two additional messages: PREEMPT and PREEMPT ACK . The state diagram in Figure 3a shows how ALLOC messages are handled, Figure 3b shows how PREEMPT messages are handled. Apart from the case where the node is IDLE described above there are two other cases to consider. If a node is already PREEMPTING, or if it is ALLOCATED and the new request has a lower precedence than the request that currently holds the resource, then the request is stored in a local WAITING queue. Otherwise, i.e., when the node is ALLOCATED and the new request has a higher precedence than the request that currently holds the resource, the algorithm performs a preemption of the resource on the request that currently holds it, named current. To perform a preemption a node sends a PREEMPT message to the node that received its resource, i.e., the node to which it previously sent an ALLOC message for current. If the node that received the PREEMPT is not the last node in the path of current, it continues the preemption to put current on hold. For this it sends a PREEMPT message to the next node in the allocation order of current. current resumes later when it becomes the request with the highest precedence in the WAITING queue.

Figure 4 shows how preemption works in the case of Req 2 and Req 3 . When n 3 receives the ALLOC message for Req 3 , it is already ALLOCATED because it has previously received the ALLOC message for Req 2 . It then sends a PRE-EMPT message along the path of Req 2 to n 6 . n 6 accepts the preemption, stores Req 2 in its local WAITING queue and sends back a PREEMPT ACK message. Req 3 then sends an ALLOC message for its last required resource to n 4 . Req 3 is now satisfied and enters its CS. When it leaves its CS, n 6 resumes Req 2 .

Further considerations on preemptions When a preemption occurs, the algorithm ensures that the resource is always released. Either the node that receives the PREEMPT message decides to release it immediately, either it waits for the current request to leave its CS. As CS have finite durations the algorithm ensures that the resource is released in finite time.

If multiple requests have a higher precedence than the one currently holding its resource, a node can receive a PREEMPT message when it is already PREEMPTING. In this case, priority is given to the request with the highest precedence.

As the communications channels are not FIFO, even if PREEMPT messages are sent along the same path than the ALLOC messages, a node may receive a PREEMPT before it has received the corresponding ALLOC . To deal with such situations, each node maintains an IGNORE list. When it receives a PREEMPT for a request and it has not yet received the corresponding ALLOC , it stores the request in the IGNORE list. When it finally receives the ALLOC for a request that is present in the IGNORE list, it ignores the ALLOC and removes the request from the list.

Experimental results

We present an evaluation of the two allocation orders reverse path and by values of the algorithm detailed above. We compare their performance through metrics that are introduced below to Dijkstra's Incremental algorithm detailed in Section 5.2 as the baseline.

Metrics and evaluation

Several metrics are used to compare the results of the algorithms:

the average usage rate is the average time during which resources are used. It is the sum of the times during which each resource is used divided by the overall duration of the experiment for all resources. 100% means that all resources are used all the time. 50% means that 50% of the resources are used on average. The objective is to maximize this metric, -the average waiting time is the average time spent by requests between the moment at which they are emitted and the moment they are satisfied.

The objective is to minimize this metric, -The average number of messages per request is the ratio between the total number of messages sent in the system for the duration of the test and the number of requests. The objective is to minimize this metric.

Dijkstra's Incremental algorithm

The baseline algorithm for our experiments is Dijkstra's Incremental algorithm [START_REF] Dijkstra | Hierarchical Ordering of Sequential Processes[END_REF]. It does not require any additional assumption on the system which allows to evaluate it against the same systems as our algorithm described above. This algorithm is selected because it gives the best average usage rate among all the state of the art algorithms evaluated. Other algorithms are not included here for the sake of space. Our implementation of the algorithm relies on the same first subroutine described above for the selection of the path, but does not require the building of the allocation vectors. The second subroutine of the Incremental algorithm relies on a static global order of the nodes. For the test system in Figure 1, our implementation considers that the global order of nodes is n 1 < n 2 < . . . < n 9 according to the subscript value. This algorithm has a drawback we call the domino effect. It is possible when the number of conflicts is high that nodes wait for each other's requests to be finished. Until these resources become available all the resources from nodes that come after in the order are unavailable. The probability of occurrence of the domino effect increases with the size of the requests.

Simulating the system with SimGrid

The algorithms have been tested on topologies from the Internet Topology Zoo [START_REF] Knight | The Internet Topology Zoo[END_REF], a dataset of topologies provided by Operators. The topologies are enriched with resources distributed across their nodes. Weights are attributed to the network links between the nodes to simulate the latency. To limit the number of parameters, all experiments in this paper use the same constant weight for all the links. The results presented here are based on the Cesnet200706 topology which consists of N = 44 nodes. This topology was selected because it has a sufficient enough size to avoid the possible bias in the results from topologies that have few nodes. Larger topologies lead to longer simulation times to get similar results. In this topology, the degrees of the nodes vary from 1 to 18, with an average degree of 2.

The simulator is based on SimGrid 3.23 [START_REF] Casanova | Versatile, Scalable, and Accurate Simulation of Distributed Applications and Platforms[END_REF]. SimGrid provides multiple Routing Models for the routing of packets between nodes. The experiments use only the Full routing model. This model requires all routes to be explicitly described in the configuration file. It ensures that no routing algorithm is used to connect two nodes that are not explicitly connected by a Link to each other. It allows the routing of the packets to be made at the application level and thus allows to simulate the first subroutine of the algorithm.

The routing tables necessary for the first subroutine are built statically by the simulator during the initialization using Dijkstra's Shortest Path algorithm [START_REF] Dijkstra | A Note on Two Problems in Connexion with Graphs[END_REF]. Routing tables are based on the Link and Hosts of the configuration.

The system is under maximum load during the tests. To achieve this, the simulation starts as many requests as there are nodes in the system, and a new request is sent as soon as a request ends. The contents of the requests are generated by the node following a linear random function. The size of the requests in a single experiment is constant. As shown by other experimental results not included here for sake of space, using a constant size does not affect the results significantly.

The duration of all experiments is the same. The time spent in CS by the requests is also constant. Empirically we settled for a duration of CS of 300, 000 and a duration of experiment of 500, 000, 000 so that time spent in CS is orders of magnitude longer than the time spent to send a message between nodes. These durations are in simulator time unit. This approximately simulates CS of 30s and a total simulation duration of 14 hours. Experiments show that this duration is long enough to make the impact of randomness negligible and the results representative.

In this experimental platform, if two allocation vectors have the same average value, the id of the nodes are used to break the tie. Since requests all have the same size, we first compare the id of the first nodes in each allocation vector. If it is the same node for both requests, we compare the id of second nodes and so on. It this comparison also results in a tie, when both requests are for the same nodes, then the internal identifiers of the requests are compared.

Figure 5 shows the results of the evaluation of the metrics defined in Section 3 for the Incremental and the algorithm detailed in Section 4 using two allocation orders : reverse path and by values in two different system configurations detailed below. The x-axis show the size of requests for a given simulation. Both axes of all the figures use a logarithmic scale. This section shows the evaluation on a system with one instance of C = 44 types of resources.

For requests of size 1, Figures 5a and5c show that the algorithm itself has no influence on the Average Usage Rate or the Average Waiting Time. In such configuration, no preemption is performed.

For requests of size 2 to 7, we can observe that Average Usage Rate decreases for all the algorithms. We reach a minimum between 7 and 9 depending on the algorithms. This can be deduced from Maekawa's proof of his quorum algorithm [START_REF] Maekawa | A N Algorithm for Mutual Exclusion in Decentralized Systems[END_REF], based on Finite Projective Planes properties, that √ N ≈ 6.63 is the minimum size so that all requests have at least one intersection. Usage rate rises up again after this lowest point because requests become larger and even if there are more conflicts between requests, each request leads to the allocation of a larger number of resources.

If the size of requests is greater than half the size of the system, then it is not possible to allocate concurrently two requests. In such situations the only possibility is to allocate the requests sequentially. The Average Usage Rate of the system then grows linearly with the size of the requests. The figure has been truncated at size 30 for the sake of readability. It is almost 100% for requests of size N as each consecutive request allocates all the resources. It is not exactly 100% due to the cost of the algorithm.

Our algorithms show an improvement on the Incremental of the Average Usage Rate with both allocation orders. The by values allocation order gives the best results for all metrics. For the Average Usage Rate we can see improvement of up to 20% from the Incremental. As shown in Figure 5b, by values does not generate more messages than the Incremental whereas reverse path shows a larger number of messages for requests of size 4 and more. This is due to the high number of preemptions taking place. The by values order limits the number of preemptions by starting with the node that has the highest counter value. This node is the most likely to receive requests with a higher precedence. Once its resource is allocated, the probability that the other nodes are preempted by requests with higher precedence gets lower. Experimental results show in Figure 5c that the by values order does not impact negatively the Average Waiting Time. The Average Waiting Time for the Incremental algorithm is significantly worse than for any of the variants. This is due to the domino effect. For the sake of space a full comparison with a near-optimal allocation is not included here, but results show that even the best results included here are 10 to 20 points lower than a near-optimal solution until requests of size 11, after which the difference starts to decrease.

System with N instances of M resources

Figure 6 shows the Average Usage Rate of the same algorithms on the same topology but with a different placement of the resources. Instead of a single instance of each resource, the system holds 4 instances for each of C = 11 different types. The figure includes three additional algorithms. Each of this additional algorithm is a variant of one of the three presented above: it uses the same allocation subroutine but a different path computation subroutine. As the path computation subroutine detailed in section 4.1 uses a static routing table, a node always selects the same node for a type of resource. The result is that the load is not well balanced across all the instances, which leads to lower Average Usage Rate. For requests of size 4, the lowest Average Usage Rate observed, the algorithm with the best result, the by values, reaches around 17%. This is lower than the worst result for the configuration with a single instance of M resources in Figure 5a where the Incremental reaches 27%.

As shown in the three additional algorithms, the load balancing improves with a simple round-robin on the different instances of each type of resource during the path computation subroutine. For example with 4 instances of c 1 , the first request selects the first instance, the second request the second one, and so on. It starts over with the fifth request that selects the first one. The selection of the route has a significant impact when there are more than one instance of the resources and improves the Average Usage Rate.

Conclusion and future works

We introduced in this paper a new algorithm for distributed Mutual Exclusion for systems with N instances of M resources, known as the k -M -N problem. This is applicable to the allocation of Multi-Domain Chains of NFs in 5G Slices. We show an improvement of up to 20% of the Average Usage Rate of resources from the baseline Incremental algorithm, with an Average Waiting Time that can be several orders of magnitude lower and no degradation of the number of messages. The results show the impact of the allocation order in which allocation is performed.The presented results focus on a few key parameters.

For the path allocation subroutine of the algorithm, the results showed the importance of the selection of the instance when N > 1. We plan to study how the performance can be further improved.

Our approach to the allocation of resources in the allocation subroutine is pessimistic, i.e., it considers that deadlocks are going to happen. But their probability of occurrence can be low in some situations. We plan to study an optimistic approach that let deadlocks occur and fix them a posteriori.

We also plan to implement the algorithm in the scheduler of a NFV platform, if possible in the live network of an Operator.

Fig. 1 .

 1 Fig. 1. A system with 3 types of resources: c1 (IDS), c2 (FW), and c3 (LB)

Fig. 3 .

 3 Fig. 3. State diagrams for ALLOC and PREEMPT messages

Fig. 4 .

 4 Fig. 4. Allocations for Req2 Req3, preemption of n6 by Req2

5. 4 Fig. 5 .

 45 Fig. 5. Evaluation for one instance of 44 types of resources

Fig. 6 .

 6 Fig. 6. Evaluation for 4 instances of 11 types of resources

The allocation of resources for VNF Forwarding Graph (VNF FG) in Multi-Domain 5G slices is addressed as a Mutual Exclusion problem for systems with