Keywords: AMS Subject classifications: 62F10, 62E17, 62P12 Skew distributions, Generalized Extreme Value distribution, Extreme Value Theory, Probability Weighted Moments, Return Level Estimation

) on the skew normal distribution, we propose an extension of the Generalized Extreme Value (GEV) distribution, the SGEV. This new distribution allows for a better fit of maxima and can be interpreted as both the distribution of maxima when maxima are taken on dependent data and when maxima are taken over a random block size. We propose to estimate the parameters of the SGEV distribution via the Probability Weighted Moments method. A simulation study is presented to provide an application of the SGEV on block maxima procedure and return level estimation. The proposed method is also implemented on a real-life data.

Skew distributions and Extreme Value Theory

The idea of modifying a symmetric probability density function of a random variable (r.v.) in a multiplicative fashion, introducing skewness, has been implemented for the univariate normal distribution by Azzalini ([2], [START_REF] Azzalini | Further results on a class of distributions which includes the normal ones[END_REF]) yielding the so-called skew-normal distribution having the probability density function of the form 2φ(x)Φ(λx), x, λ ∈ R (where φ(•) and Φ(•) denote respectively the probability distribution function (pdf) and cumulative distribution function (cdf) of the standard normal distribution). This distribution presents the advantage to look like a normal one in the core data but allowing some asymmetry. The model has since been extended to the gaussian related distributions (multivariate normal or Student distribution) and its extensions are provided in the book edited by Genton [START_REF] Genton | Skew-Elliptical Distributions and Their Applications: A Journey Beyond Normality[END_REF] and in the review paper [START_REF] Azzalini | The skew-normal distribution and related multivariate families[END_REF] and the monograph by Azzalini [START_REF] Azzalini | The Skew-Normal and Related Families[END_REF] who regularly updates list of references on his web page. The skew-normal distribution has been also applied in many fields like climatology (Flecher et al. [START_REF] Flecher | A Stochastic Daily Weather Generator for Skewed Data[END_REF]) or insurance (Vernic [START_REF] Vernic | Multivariate skew-normal distributions with applications in insurance[END_REF]) and adapted to classical statistical tools like Space-state models (Naveau et al. [START_REF] Naveau | A skewed Kalman filter[END_REF]) or random fields (Allard and Naveau [START_REF] Allard | A New Spatial Skew-Normal Random Field Model[END_REF]). Since then, several authors have tried to generalize these results to skewing arbitrary symmetric (or not) pdf's. See Gupta et al. [START_REF] Gupta | Some skewsymmetric models[END_REF] for a most detailed discussion of skew-symmetric models based on the normal, Student's t, Cauchy, Laplace, logistic and uniform distributions. In that case, the main feature of these models is that a new parameter λ is introduced to control skewness and kurtosis. Nadarajah and Kotz ( [START_REF] Nadarajah | Skew distributions generated from different families[END_REF], [START_REF] Nadarajah | Skew Models I[END_REF] and [START_REF] Nadarajah | Skew Models II[END_REF]) proposed some skew models based on Cauchy, Laplace, Logistic, Bessel and Pearson (among others). Some relations between Extreme Value Theory and Skew Distribution have already been studied by Chang and Genton [START_REF] Chang | Extreme Value Distributions for the Skew-Symmetric Family of Distributions[END_REF] in the case of skewsymmetric family of distributions and by Lysenko et al. [START_REF] Lysenko | Multivariate extremes of generalized skew-normal distributions[END_REF] in the case of multivariate generalized skew-normal distribution.

In the same spirit, we propose here an extension of the Generalized Extreme Value (GEV) distribution. This distribution, introduced by Jenkinson [START_REF] Jenkinson | The frequency distribution of the annual maximum (or minimum) values of meteorological elements[END_REF], is the limiting distribution for extremes, both normalized maxima and minima of a random sample. In a more formal way, it is known that the only possible limiting form of a normalized maximum of a random sample when a non-degenerate limit exists is captured by the GEV

G(x; σ, γ, µ) =        exp -(1 + γ x-µ σ ) -1/γ , if 1 + γ x-µ σ > 0, γ = 0, exp -exp(-x-µ σ ) , if x ∈ R, γ = 0,
where µ ∈ R, σ > 0 and γ are called the location, scale and shape parameters, respectively.

The GEV model has been applied to environmental sciences to fit, for example, yearly maximal temperatures or yearly maximal river discharges. However, the limiting distribution does not always provide a good fit in the finite samples that arise in practice, since the extreme value theory is based on asymptotic arguments.

Our purpose, here, is to generalize the GEV distribution by addition of an extra parameter, λ, introduced to model flexibility in the data and so better adapt to real situations. The main features of this distribution will be discussed together with two methods of estimation of its parameters.

The paper is organized as follows. In Section 2, we will define the distribution, study some of its characteristics and present some possible interpretations. In Section 3, the parameters of the distribution will be estimated by the classical methods of maximum likelihood and probability-weighted moments. Finally, a simulation study and an application to real data will be performed (Section 4 and 5) to investigate the behavior of the estimator in small and moderate samples.

The skew generalized extreme value distribution

Definition 1. A random variable X is said to have a Skew Generalized Extreme Value (SGEV) distribution if its probability density function is given by f

(x) = (λ + 1)g(x)G λ (x) (1)
where λ > -1, g(.) and G(.) are respectively the density function and the cumulative distribution function of a GEV distribution with location parameter µ, scale parameter σ and shape parameter γ.

The cumulative distribution fucntion is given by

F (x) = G λ+1 (x).
The distribution will be denoted for short by SGEV (µ, σ, γ, λ). It follows from (1) that the pdf of X can be written as

f (x) = λ+1 σ 1 + γ x-µ σ -1/γ-1 exp -(1 + λ) 1 + γ x-µ σ -1/γ if 1 + γ x-µ σ > 0, γ = 0, ( λ+1 σ ) exp -(λ + 1) exp -x-µ σ exp -x-µ σ if x ∈ R, γ = 0.
When λ = 0, equation ( 1) reduces to the GEV pdf. The parameter λ allows for a greater degree of flexibility and we can expect this to be useful in many more practical situations.

The quantile function of the cdf F is given by

Q(u) =      µ + σ γ -log u λ+1 -γ -1 if γ = 0, µ -σ log -log u λ+1 if γ = 0
and allows to compute quantiles in an easy way.

We also find that the expectation, the second-order moment and the variance of X are respectively given by

E(X) = µ - σ γ [1 -(λ + 1) γ Γ(1 -γ)] , E(X 2 ) = µ - σ γ + σ γ(λ + 1) -γ Γ(1 -γ) 2 , V (X) = σ 2 γ 2 (λ + 1) -2γ Γ(1 -2γ) -Γ 2 (1 -γ) , if λ > -1
, where Γ denotes the Gamma function Γ(t) = ∞ 0 e -u u t-1 du, t > 0. Simple modifications of Chang and Genton [START_REF] Chang | Extreme Value Distributions for the Skew-Symmetric Family of Distributions[END_REF] leads to the conclusion that the extreme value index of the SGEV distribution is γ.

The SGEV distribution offers some interesting practical interpretations which are going to be presented in the following remarks.

Remark 1. Let X 1 , X 2 , . . . be i.i.d. r.v.'s with cdf F and let N be a Poisson r.v. with intensity λ * > 0. Moreover, let N be independent of {X i } i≥1 . Define the following maximum with a random block size

M = max(X 1 , . . . , X N ) if N ≥ 1, 0 if N = 0.
Then, we have

P (M ≤ x) = +∞ n=0 F n (x) λ * n e -λ * n! = exp {-λ * [1 -F (x)]} = exp -F (x) λ * . Now, denote G(x) = exp[-F (x)]. It follows P (M ≤ x) = G λ * (x) (2)
which corresponds to SGEVs with strictly positive parameters λ * = λ + 1.

Define the cdf H(x; σ, γ, µ) by

H(x; σ, γ, µ) =      1 -1 + γ( x-µ σ ) -1/γ , if γ = 0, 1 -exp -( x-µ σ ) , if γ = 0,
where x ≥ µ if γ ≥ 0, µ ≤ x ≤ µ -σ/γ if γ < 0 and µ ∈ R is the location parameter, σ > 0 the scale parameter and γ ∈ R the shape parameter. This distribution, which appears as the limit distribution of scaled excesses over high thresholds, is called a Generalized Pareto Distribution (GPD) (Pickands [25]).

If F is a GPD tail, then the cdf of M is the cdf of a SGEV on the interval [µ, ∞[, as it is shown in Figure 1. If x ≤ µ, then P(M ≤ x) = P(N = 0) = exp(-λ * ). The equality (2) offers then a first interpretation for the SGEV distribution: the cdf G λ * (x), λ * = λ + 1, on the interval [µ, ∞[ can be viewed as the cdf of a maximum taken over exceedances in a random Poisson block size.

Remark 2. The SGEV is also the distribution of the maxima when they are taken on dependent data. In that case, λ is between -1 and 0. The extremal index θ (see, for example, Beirlant et al. [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]) is then equal to λ + 1. The SGEV is more general since the value of λ needs only to be greater than -1. So, the family of SGEV includes the classical GEV one, but also the GEV one based on dependent data. This is quite interesting because this property could be used to define a new plug-in estimator of the extremal index, θ = λ + 1, where λ is the estimator of λ obtained by the PWM method.

Parameter estimation

Let (X 1 , . . . , X n ) be a sample of SGEV r.v.s. In this section, the parameters (µ, σ, γ, λ) of the SGEV distribution will be estimated by both the method of maximum likelihood (ML) and the probability weighted moments (PWM) method.

Maximum likelihood

In case γ = 0, the log-likelihood function for a sample (X 1 , . . . , X n ) of SGEV r.v.s is given by

logL(µ, σ, γ, λ) = n log(λ + 1) -n log σ - 1 γ + 1 n i=1 log 1 + γ x i -µ σ -(λ + 1) n i=1 1 + γ x i -µ σ -1/γ . If γ = 0, the log-likelihood function reduces to logL(µ, σ, λ) = n log(λ + 1) -n log σ -(λ + 1) n i=1 exp - x i -µ σ - n i=1 x i -µ σ .
The ML estimator (μ, σ, γ, λ) is obtained by maximizing the log-likelihood function. Clearly, no explicit solution exists to these equations so that numerical procedures are called for. For a discussion about the properties of the MLE see Smith [START_REF] Smith | Maximum likelihood estimation in a class of nonregular cases[END_REF]. In the paper by Smith [START_REF] Smith | Maximum likelihood estimation in a class of nonregular cases[END_REF], the ML estimators of the parameters of the GEV distribution are shown to be asymptotically normal whenever γ > -1/2. Since the SGEV distribution is based on the GEV distribution, it is easy to check that the classical properties of the ML estimator of (µ, σ, γ, λ) hold whenever γ > -1/2. Remark 3. The ML method usually gives poor results when it has to estimate a skew parameter: for example, for moderate sample sizes, the ML estimator is infinite with positive probability when dealing with skew normal distribution (Sartori [26]). The issue mentioned here has a close connection with the Fisher singularity problem discussed in Hallin andLey (2012a, 2012b). In Section 4.1 and in Tables 1 to 4, the results obtained by estimating the parameters by Maximum Likelihood are disappointing. This is the reason why, when estimating the return level in Section 4.2, we will avoid the ML method and prefer the Probability Weighted Moments one, which will be described in the next section.

The method of probability weighted moments

The probability weighted moments, introduced by Greenwood et al. [START_REF] Greenwood | Probability weighted moments: definition and relation to parameters of several distributions expressable in inverse form[END_REF], are a simple generalization of the usual moments of a probability distribution. They are defined in terms of the cdf H of a r.v. Y as

E{Y p [H(Y )] r [1 -H(Y )] s }
with p, r and s real numbers. Particularly useful special cases are the probability-weighted moments

β r = E{Y [H(Y )] r }
which have been used by [START_REF] Hosking | Estimation of the Generalized Extreme-Value distribution by the method of probabilityweighted moments[END_REF] to estimate the parameters of the GEV distribution. For the GEV distribution, they equal

β r = (r + 1) -1 µ - σ γ 1 - Γ(1 -γ) (1 + r) -γ γ < 1, γ = 0. ( 3
)
This method of moments has been also applied to the truncated skewnormal distribution (Flecher et al. [START_REF] Flecher | Estimating the Closed Skew-Normal distributions parameters using weighted moments[END_REF]). It can be equally applied to estimate the parameters of the SGEV distribution. Let X be a SGEV (µ, σ, γ, λ) r.v. with cdf F and let Y be a GEV (µ, σ, γ) r.v. with cdf G. It is easy to check that the probability weighted moments of a SGEV distribution, denoted by ν a , are special probability weighted moments of a GEV distribution:

ν a = E[XF a (X)] = xF a (x)f (x)dx = xG (λ+1)a (x)(λ + 1)g(x)G λ (x)dx = (λ + 1) xG (λ+1)a+λ (x)g(x)dx = (λ + 1)E[Y G λ(a+1)+a (Y )] = (λ + 1)β λ(a+1)+a
By taking into account (3), we obtain

ν a = E[XF a (X)] = 1 a + 1 µ - σ γ 1 - Γ(1 -γ) ((λ + 1)(a + 1)) -γ .
In practice, we just need to choose four different values of a and so set up the necessary number of moment equations. Here, we choose a = 0, 1, 2, 3, because they are trivial to estimate. The estimators of γ, σ, µ and λ are then obtained by solving the following system of equations:

                         ν 0 = µ -σ γ [1 -(λ + 1) γ Γ(1 -γ)] , 2ν 1 = µ -σ γ [1 -2 γ (λ + 1) γ Γ(1 -γ)] , 3ν 2 = µ -σ γ [1 -3 γ (λ + 1) γ Γ(1 -γ)] , 4ν 3 = µ -σ γ [1 -4 γ (λ + 1) γ Γ(1 -γ)] , (4)
after replacing ν a by its empirical version νa :

νa = 1 n n i=1 X (i) i n a where X (1) ≤ X (2) ≤ ... ≤ X (n)
represent the ordered sample. This system is quiet hard to solve numerically.

A promising avenue consists of first estimating γ by solving the following equation 3ν

2 -ν0 2ν 1 -ν0 = 3 γ -1 2 γ -1 . (5)
Once γ given, find the estimations of λ and σ by solving the following system of equations

         2ν 1 -ν0 = σ γ Γ(1 -γ)( λ + 1) γ (2 γ -1) , 4ν 3 -3ν 2 = σ γ Γ(1 -γ)( λ + 1) γ (4 γ -3 γ ) , (6)
and finally, σ, γ and λ given, obtain the estimation of µ, that is

μ = ν0 + σ γ 1 -( λ + 1) γ Γ(1 -γ) .
The case γ = 0 can also be covered by this method: a system of three equations with three unknown parameters will then be obtained by taking the first three equations in the system (4) and letting γ go to 0. Nevertheless, in the case γ = 0, all moments exist and then, we do not need a new method based on moments to get good results since the classical moments method will be satisfactory.

To derive the limiting distribution of (μ, σ, γ, λ), we need to study the limiting behavior of νa . By direct application of the central limit theorem, one can easily prove that νa is asymptotically normally distributed with mean ν a and asymptotic variance

(λ + 1) 2 n σ 2 γ 2 [(a + 1)(λ + 1)] 2γ Γ(1 -2γ)K a(λ + 1) + λ (a + 1)(λ + 1) -Γ 2 (1 -γ)
where K denotes the hypergeometric function

K(x) = 1 F 2 (k, 2k; 1 + k; -x).
A similar argument applies to any linear combination of the νa (a = 0, 1, 2, . . .) and it follows that the νa are asymptotically jointly normal.

Theorem 1. Let (X 1 , . . . , X n ) be a sample of SGEV. If γ < 0.5, then the rescaled PWM estimator vector defined by (4) and denoted by

√ n     ν0 ν1 ν2 ν3     converges in distribution, as n → ∞, towards N         ν 0 ν 1 ν 2 ν 3     , Ω    
where Ω denotes the covariance matrix.

The elements of the covariance matrix Ω are given in Appendix A.

Corollary 2. The estimators of µ, σ, γ and λ obtained by the probability weighted moments method are asymptotically normally distributed and are asymptotically unbiased if γ < 0.5.

The proof of the corollary is a direct application of Theorem 1 and the Delta Method.

Remark 4. Moments ν 0 , ..., ν 3 have been chosen in the estimation procedure of the parameters of the SGEV distribution. The same moments could be used to obtain different moment methods. Equations similar to (5) could be written down with any of the 4 sets of 3 moments, among ν 0 , ..., ν 3 . Similarly, any pair of moments among the 6 possible pairs can lead to a system of two equations similar to that in [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. Finally, µ can also be estimated using any of the 4 moments. The estimation variance of the selected combination of moments is the lowest among all possible combinations, since the moments of lowest degree have been selected. This procedure can then be improved using all combinations simultaneously. Consider for example the estimation of γ. Four equations similar to (5) can be written using each triplet of moments. An estimator γ can be defined as the value minimizing the mean squared errors between the left hand side and the right hand side of these equations. A similar approach can then be implemented for estimating σ and λ given γ, and finally for estimating µ given all other estimates.

A simulation study 4.1 Parameters estimation

We perform a simulation study in order to compare the performances of the ML and PWM estimators of the parameters of the SGEV distribution.

Simulations were performed for sample size n = 100 with γ taking values -0.5, 0, 0.5, 1 and 1.5. Without loss of generality, the location and scale parameters were set to µ = 0 and σ = 1. We set also the skew parameter λ to 0.2, 0.5, 1 and 1.5. For each combination of values of γ and λ, we generated 100 random samples from the SGEV distribution and for each sample, we estimate the parameters µ, σ, γ and λ by ML and PWM method. The results are listed in Tables 1 to 4.

Remark that some chosen values for γ are not in the range of validity of the estimation methods. In practice, on real data, the value of γ is unknown, so values which are not in the validity range were used to discuss the robustness of the method. For example, when γ is greater than 0.5 but lesser than 1, moments still exist but the estimators are no more asymptotically normal. In this case, the PWM estimator gives acceptable results when estimating µ and σ (sometimes λ). On the other side, when γ ≥ 1 the PWM method is not satisfactory anymore.

As previously mentioned, the Maximum Likelihood method does not perform well especially when estimating the skew parameter λ. Indeed, in the case γ < 0.5, we can notice the superiority of the PWM method for estimating the skew parameter λ, while ML and PWM provides good and comparable estimations for the three other parameters. In particular, because of the quite large sample size (n = 100) in our simulation study, the ML estimator is preferable to estimate γ (when γ < 0.5). On the other side, both methods do not perform well when γ gets greater (if γ > 0.5), especially for the skew parameter. For the PWM method, this is due to the usual regularity conditions underlying the asymptotic properties of the estimators (γ < 0.5) while for the ML method this may be due to some algorithmic problems. This shows that both methods are not really robust when used outside their range of validity. For a more detailed discussion about the performance of both ML and PWM estimators, see Hosking and Wallis [START_REF] Hosking | Estimation of the Generalized Extreme-Value distribution by the method of probabilityweighted moments[END_REF]. Remark 5. When looking at the results of the simulation study, it appears that the parameter γ is always fairly well estimated by ML. This could suggest to consider a mixed estimation procedure which will consist in estimating γ with the ML method and the other parameters with the PWM method. This kind of procedure would lead to a vector of estimators whose asymptotic law could be quite hard to obtain. An alternative could be to consider the estimation obtained by the PWM method as a starting value for the ML method.

Return level estimation

Let Z 1 ,...,Z n be a sample of i.i.d. random variables from a distribution function H. We are now interested in estimating z p := H -1 (1 -p). The method of Block Maxima consists in replacing the original sample Z 1 , ..., Z n by Y 1 , ..., Y m where Y i is the maximum on the i th block of size k i.e. Y i = max{Z (i-1)k+1 , ...., Z ik } with n = m × k. Since the r.v.'s Y i are asymptotically GEV distributed, the return level can be estimated using (see Beirlant et al. [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF])

ẑ * p =      µ Y + σ Y γ Y -log(1 -p) k -γ Y -1 , if γ Y = 0, µ Y -σ Y log -log(1 -p) k , if γ Y = 0, (7)
leading to the GEV estimator ẑGEV If Z 1 ,...,Z n are i.i.d. with gaussian distribution N (0, 1) we also have (see Embrechts et al. [START_REF] Embrechts | Modelling Extremal Events for Insurance and Finance[END_REF])

P Z (n,n) -d N n c N n ≤ z → GEV (µ = 0, σ = 1, γ = 0)
where Z (n,n) is the maximum of the sample,

c N n = (2 log n) -1/2 ,
and

d N n = 2 log n - log(4π) + log log n 2(2 log n) 1/2 .
The distribution of Y i can be approximated by a Obviously, the theoretical estimation z BM p will get closer to z p as the block size increases. In practice, we need to choose a block size large enough to ensure a good approximation of the GEV distribution but allowing also enough maxima observations to obtain a good parameter estimation. In that way, the flexibility of the SGEV distribution can give a good compromise. Once the sample of maxima obtained, the return level can be estimated by

GEV (µ = d N k , σ = c N k , γ = 0)
ẑSGEV p =      µ s + σs γs -log(1 -p) k/(λs+1) -γs -1 , if γ s = 0, µ s -σ s log -log(1 -p) k/(λs+1) , if γ s = 0, (8)
where µ s , σ s , γ s et λ s are the parameter estimators obtained after fitting a SGEV on the sample of maxima Y 1 , . . . , Y m .

Similar results stand for the Log-gamma distribution LΓ(α, β) with

c LΓ n = (Γ(β)) -1 (log n) β-1 n 1/α and d n = 0.
We want to answer the two following questions : in the case of slow convergence of the maximum towards the GEV distribution (e.g. the Gaussian case), how the SGEV behaves and how does the SGEV for heavy tailed distributions? Figures 2 and3 represent the estimation of the quantile as a function of the block size. For each chosen distribution (N (0, 1) in Figure 2 and Loggamma(2,2) in Figure 3), we simulated 500 samples of size n = 500 or n = 5000. For each combination, we compute the theoretical value z BM p (red line) to show how good the approximation by a GEV distribution is, the estimation of the quantile based on the maximum by block using both GEV (ẑ GEV p , grey line) and SGEV (ẑ SGEV p , black line) and the 95% empirical confidence intervals (dashed lines). The true value of the quantile is represented by the horizontal line. We also represent, on the bottom row of each figure, the value of the skew parameter estimation λ s , together with a 95% empirical confidence interval. By a first look at these simulations, we observe that the return level estimation is not so good when dealing with Gaussian distributions (cf Figure 2). This is not surprising since it is well known that the normal distribution converges very slowly to the GEV distribution, but it is also an indicator of the limitations of the Extreme Value Theory. Although the SGEV approximation is not great, it seems to propose an improvement in terms of bias and variance.

When dealing with a heavy tailed distribution (cf Figure 3), the quality of the return level estimation using the SGEV approximation gets better even if the approximation of the sample of maxima seems good (see the theoretical z BM p ).

In view of these simulations, the SGEV proposes a complementary method to the classical Extreme Value Theory by improving the return level estimation both in the case of a slow convergence or in the heavy tailed case.

An application to real data

We apply the proposed procedure to the Fort Collins precipitation data. The data consist of the annual daily precipitation maxima (in mm) recorded by a weather station from 1948 to 2001 in Fort Collins (Colorado, USA) (cf Guillou et al. [START_REF] Guillou | Return level bounds for discrete and continuous random variables[END_REF] or Diebolt et al. [START_REF] Diebolt | Improving Probability-Weighted Moments Methods for the Generalized Extreme Value Distribution[END_REF] for more details and explanation on this data set).

The values of the estimated parameters, obtained by ML and PWM methods, are listed in Table 5 when fitting a SGEV distribution or a GEV distribution. For the GEV distribution, the estimation of the return level associated to a period of 100 years is equal to 573.72 (when adopting the PWM method) and 653.33 (when adopting the ML method). By fitting the SGEV distribution, we obtained an estimation of the return level period equal to 608.83 (with the PWM method) and 653.21 (with the ML method). This application was selected to illustrate a number of practical points. First of all, we remark that if we look at the same estimation method (either ML or PWM method) the return level estimation does not vary so much if we take a GEV or a SGEV tail approximation. This is probably due to the fact that the annual maximum is taken on a sample of important size (equal to 365) and so the GEV approximation is good enough to provide interesting results. Indeed, the values of λ are near 0 indicating that the GEV distribution may be a good choice of fit for this sample. In addition, we remark that there is no loss using SGEV instead of GEV distribution. So, in our opinion, the first one should be preferred since it englobes more cases than the classical GEV distribution. Indeed, the SGEV distribution could be used to define the right block size. One can try different block size and choose the one which gives the value of the parameter Authors would like to thank two anonymous referees for their helpful comments that led to an improvement of the present paper (we especially thank for Remarks 4 and 5).

A Asymptotic covariance matrix

The asymptotic covariance matrix of the estimator vector obtained by the probability weighted moments method (see Theorem 1) can be written as

Ω = (λ + 1) 2 σ 2 γ 2     Ω 0,0 Ω 0,1 Ω 0,2 Ω 0,3 Ω 0,1 Ω 1,1 Ω 1,2 Ω 1,3 Ω 0,2 Ω 1,2 Ω 2,2 Ω 2,3 Ω 0,3 Ω 1,3 Ω 2,3 Ω 3,3    
and the elements on the main diagonal of Ω can be specified as follows -(λ(a + 2) + a + 1) γ Γ(1 -2γ)K (λ + 1)(a + 1) λ(a + 2) + a + 1 +2 ((λ + 1)(a + 1)) γ × (λ(a + 2) + a + 1) γ ((λ + 1)(a + 2)) γ Γ 2 (1 -γ) , Ω a,a+2 = 1 2 ((a + 3)(λ + 1)) 2γ Γ(1 -2γ)K λ(a + 1) + a (λ + 1)(a + 3)

Ω 0,0 = (λ + 1) 2γ Γ(1 -2γ)K λ λ + 1 -Γ 2 (1 -γ) ,
-(λ(a + 3) + a + 2) γ Γ(1 -2γ)K (λ + 1)(a + 1) λ(a + 3) + a + 2 +2 ((λ + 1)(a + 1)) γ × (λ(a + 3) + a + 2) γ ((λ + 1)(a + 3)) γ Γ 2 (1 -γ) , Ω a,a+3 = 1 2 ((a + 4)(λ + 1)) 2γ Γ(1 -2γ)K λ(a + 1) + a (λ + 1)(a + 4)

-(λ(a + 4) + a + 3) γ Γ(1 -2γ)K (λ + 1)(a + 1) λ(a + 4) + a + 3 +2 ((λ + 1)(a + 1)) γ × (λ(a + 4) + a + 3) γ ((λ + 1)(a + 4)) γ Γ 2 (1 -γ) .

p

  after replacing µ Y , σ Y and γ Y by μY , σY and γY , the ML or PWM estimators of the GEV parameters based on the sample Y 1 , ..., Y m .

  and the return level of probability p can be approximated by z BM p after replacing µ Y and σ Y by d N k and c N k , and γ Y by 0 in Equation (7) (in the case γ Y = 0).

Figure 1 :

 1 Figure 1: Cumulative distribution function of M if X i ∼ GP D(µ = 10, σ = 1, γ = 0) and λ = 1 (full line) of the Remark 1 and cumulative distribution function of the GEV distribution (dashed line).

Ω 1 , 1 =

 11 (2λ + 2) 2γ Γ(1 -2γ)K 2λ + 1 2λ + 2 -Γ 2 (1 -γ) , Ω 2,2 = (3λ + 3) 2γ Γ(1 -2γ)K 3λ + 2 3λ + 3 -Γ 2 (1 -γ) , Ω 3,3 = (4λ + 4) 2γ Γ(1 -2γ)K 4λ + 3 4λ + 4 -Γ 2 (1 -γ) ,where K and Γ are respectively the hypergeometric and gamma functions. The remaining elements of the matrix areΩ a,a+1 = 1 2 ((a + 2)(λ + 1)) 2γ Γ(1 -2γ)K λ(a + 1) + a (λ + 1)(a + 2)

Table 1 :

 1 For each combination of µ, σ and γ, we generate 100 samples of SGEV(µ,σ,γ,λ = 0.2) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Parameter True value Mean for ML Mean for PWM sd for ML sd for PWM
	µ	0	0.184	-0.014	0.452	0.103
	σ	1	0.908	0.953	0.288	0.080
	γ	-0.5	-0.517	-0.433	0.065	0.075
	λ	0.2	0.038	0.213	0.620	0.028
	µ	0	0.097	0.009	0.075	0.105
	σ	1	1.020	0.955	0.086	0.089
	γ	0	0.0168	0.039	0.051	0.072
	λ	0.2	0.066	0.196	0.086	0.007
	µ	0	-0.131	0.018	0.086	0.093
	σ	1	0.921	1.002	0.103	0.126
	γ	0.5	0.511	0.467	0.107	0.104
	λ	0.2	0.383	0.192	0.128	0.069
	µ	0	-0.317	-0.024	0.068	0.264
	σ	1	0.666	1.248	0.096	0.331
	γ	1	1.047	0.772	0.138	0.109
	λ	0.2	0.773	0.364	0.191	0.316
	µ	0	-0.399	-0.528	0.058	0.805
	σ	1	0.388	1.780	0.087	0.877
	γ	1.5	1.492	0.897	0.159	0.085
	λ	0.2	1.317	0.916	0.253	0.937
	λ nearest to zero.
	6 Conclusions
	A new class of distributions, called SGEV distributions, has been proposed
	as a generalization of the classical GEV distribution, thanks to the intro-
	duction of a new parameter λ. Estimators of the parameters and quantiles
	of the SGEV distribution have been derived by ML and PWM methods.
	As mentioned in the above sections, the SGEV distribution is more flexible
	than the classical GEV distribution and so, can fit better maxima when deal-
	ing with small maxima samples or when the maxima come from small size

Table 2 :

 2 For each combination of µ, σ and γ, we generate 100 samples of SGEV(µ,σ,γ,λ = 0.5) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Parameter True value Mean for ML Mean for PWM sd for ML sd for PWM
	µ	0	0.022	0.006	0.445	0.094
	σ	1	1.060	0.951	0.353	0.128
	γ	-0.5	-0.349	-0.439	0.786	0.094
	λ	1	0.124	1.005	1.330	0.027
	µ	0	0.241	0.009	0.249	0.089
	σ	1	1.048	0.965	0.177	0.100
	γ	0	0.196	0.031	0.424	0.069
	λ	1	0.383	0.997	0.331	0.004
	µ	0	0.157	-0.014	0.107	0.105
	σ	1	1.081	1.043	0.116	0.173
	γ	0.5	0.513	0.454	0.128	0.117
	λ	1	0.702	1.004	0.186	0.044
	µ	0	0.011	-0.273	0.074	0.429
	σ	1	0.988	1.326	0.106	0.444
	γ	1	1.028	0.763	0.132	0.119
	λ	1	0.958	1.408	0.162	0.830
	µ	0	-0.100	-2.445	0.085	1.914
	σ	1	0.827	2.334	0.166	1.825
	γ	1.5	1.514	0.894	0.169	0.077
	λ	1	1.262	3.123	0.219	2.416

Table 3 :

 3 For each combination of µ, σ and γ, we generate 100 samples of SGEV(µ,σ,γ,λ = 1) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Parameter True value Mean for ML Mean for PWM sd for ML sd for PWM
	µ	0	0.137	-0.006	0.272	0.098
	σ	1	0.936	0.935	0.222	0.119
	γ	-0.5	-0.523	-0.421	1.005	0.085
	λ	1.5	0.185	1.507	0.808	0.0198
	µ	0	0.179	0.009	0.345	0.100
	σ	1	1.029	0.956	0.192	0.126
	γ	0	0.591	0.027	0.746	0.069
	λ	1.5	0.596	1.497	0.755	0.005
	µ	0	0.185	0.011	0.262	0.124
	σ	1	1.114	1.021	0.180	0.206
	γ	0.5	0.694	0.471	0.407	0.123
	λ	1.5	1.089	1.508	1.991	0.058
	µ	0	0.138	-0.461	0.103	0.676
	σ	1	1.104	1.389	0.154	0.522
	γ	1	1.041	0.770	0.158	0.121
	λ	1.5	1.222	2.219	0.188	1.473
	µ	0	0.006	-2.905	0.070	2.740
	σ	1	0.969	2.732	0.134	1.403
	γ	1.5	1.509	0.904	0.171	0.079
	λ	1.5	1.501	3.127	0.193	2.492

Table 4 :

 4 For each combination of µ, σ and γ, we compute 100 samples of SGEV(µ,σ,γ,λ = 1.5) and calculate the PWM and ML estimations, together with the respective mean and standard deviation (sd).

	Model Method	µ	σ	γ	λ
	SGEV	ML	119.41 50.22 0.35 -0.15
		PWM 107.23 49.99 0.30 0.05
	GEV	ML	111.31 47.39 0.35	-
		PWM 112.47 50.57 0.27	-

Table 5 :

 5 Estimations of µ, σ, γ and λ by ML and PWM methods for the Fort Collins data set, obtained by fitting a SGEV or a GEV distribution.

Acknowledgements

This work has been partially supported by the MIRACCLE-GICC project, the French Research National Agency (ANR) under the reference ANR-08-BLAN-0314-01, the EU-FP7 "ACQWA" Project (www.acqwa.ch) under contract Nr 212250, the ANR-MOPERA project, the PEPER-GIS project, the ANR-McSim project and the research chair Actuariat responsable sponsored by Generali.

blocks. Two interpretations of the SGEV distribution can be given. First, the cdf can be viewed as the cdf of a maximum over a random block size when the maximum is taken over exceedances (see Remark 1). Secondly, the SGEV distribution can be seen as the limiting distribution of the maxima when the original data are not independent. In that case, the parameter λ is just θ -1, where θ is the extremal index. This means that the estimator of λ, obtained by the PWM method, can be looked at as a new and efficient estimator of the extremal index θ.

Moreover, the SGEV could be used to find the right block size, the one such that λ is nearest to zero.


Hang Wan

Michael David

William Derigent

Modelling Digital Twins as a recursive Multi-Agent Architecture: application to energy management of communicating materials

Keywords: Communicating materials, digital twin, recursive multi-agent system, architecture modelling

In recent years, there has been an increasing interest in intelligent product especially for industrial applications. A special intelligent product was proposed ten years ago, and this paradigm was named "communicating material", i.e. materials equipped with micro components able to sense, store and communicate data. Communicating material needs a recursive architecture to be managed with a digital twin. A literature analysis is proposed concerning recursive multi-agent systems and architectures regarding communicating material applicability. The most promising reference architecture was adapted for modelling digital twin of a communicating material under energy constraint. The composition and decomposition mechanisms are detailed to improve energy performance by generating an aggregated product with a global view. An applicative scenario relative to the different aggregation modes between two communicating materials, illustrates the proposed architecture ability.

INTRODUCTION

The rapid development of Internet of Things (IoT) leads to the intelligent product paradigm [START_REF] Wong | The intelligent product driven supply chain[END_REF]. In this vision, manufactured products can complete tasks independently and can also cooperate with other products to respond to changes in environment. This adaptive ability makes it widely used in industry (from manufacturing to logistics) [START_REF] Meyer | Intelligent Products: A survey[END_REF].

In 2009, the "communicating material" (CM) concept was proposed by [START_REF] Kubler | Problem definition methodology for the "communicating material" paradigm[END_REF]. It is a special type of intelligent product. This material can communicate with environment. More specially, they can process, store monitoring data in their own structure and exchange information with other CMs or users. Compared with traditional intelligent product, CM has two more special abilities: the capability of being intrinsically and wholly communicating which support physical transformation (decomposition/composition); and the capability of managing its own data. This concept has been applied in different works. A recent example is the communicating concrete in the McBIM project [START_REF] Derigent | Materials communicating with the BIM: Aims and first results of the MCBIM project[END_REF]. A Wireless Sensor Network (WSN) is embedded into a concrete component to monitor, process and store structural health information for construction industry. CM works in a resource constrained environment which brings challenges in data and energy management [START_REF] Wan | Energy-efficient chain-based data gathering applied to communicating concrete[END_REF]. Indeed, in this project, the communicating concrete should work from manufacturing to exploitation (at least ten years).

To meet this long-term service requirement, solutions are required to monitor and forecast remaining energy of the WSN embedded in the materials. Distributed solutions over the network are often suboptimal. Using a virtual representation could allow to achieve a global optimization. Indeed, applying the concept of Digital Twin (DT) seems to be a suitable solution, since it would bring a WSN virtual model to store all information concerning the physical part and provide optimal behavior decisions (with a centralized solution) to extend its lifetime. DT concept was first introduced by Michael Grieves in 2002 and formalized in 2016: a digital twin is a set of virtual information constructs that describes the physical product. It includes all information that could be obtained from the physical part [START_REF] Grieves | Digital Twin: Mitigating Unpredictable, Undesirable Emergent Behavior in Complex Systems[END_REF]. Many other DT definitions have been proposed. Form Rosen, a DT is realistic models of the current state of process and behaviors which can interact with real environment [START_REF] Rosen | About the importance of autonomy and digital twins for the future of manufacturing[END_REF]. [START_REF] Gabor | A Simulation-Based Architecture for Smart Cyber-Physical Systems[END_REF] define a DT as the simulation of physical object to predict future state of the system. A comparison of these different views is provided by Negri [START_REF] Negri | A Review of the Roles of Digital Twin in CPS-based Production Systems[END_REF]. In our definition, a DT is a digital representation of real CM products, which uses the input data to mirror, predict and provides energy efficient solutions to extend CM lifetime.

Fig. 1. Digital Twin for communicating concrete

The DT example for communicating concrete is shown in Fig. 1. Real communicating nodes (CNs) are represented by agents in concrete DT. Colors of agent represent energy level of CN from 100% (blue) to 0 (black). Node agents are grouped in a concrete agent which can store all product information and monitor data of the real concrete. Concrete agent can provide optimal orders with a global view to meet environment changes. These orders could be a monitoring frequency change, a new aggregation technique application at node or a different routing strategy for data collection.

During all lifecycle phases, each concrete element can manage its own data and make own decision to accommodate environmental changes. Besides, it can aggregate with other elements into an aggregated product with comprehensive knowledge of energy information. This flexible aggregation ability brings challenges to their dedicated DTs: How to construct a DT with many agents (multi-agent) and with many abstraction levels? How to realize the recursion that aggregate multiple independent DTs into an "aggregated DT" to improve performance and extend lifetime?

To address the mentioned issues, this paper aims to design a recursive multi-agent architecture to support this kind of management. In next section, the recent recursive holarchy and multi-agent architecture are reviewed. The CM recursive agent model is detailed in section 3. A global view of the recursive architecture and the composition/decomposition mechanisms are presented in section 4. CM application scenario is given in section 5. Finally, section 6 concludes this paper and presents some perspectives for future work.

RELATED WORK

The concept of intelligent product has been widely applied in our daily life. The recursion of multi products has attracted a lot of attention, especially for the design of control system. [START_REF] Mortellec | Embedded holonic fault diagnosis of complex transportation systems[END_REF] present SURFER, an embedded holonic fault diagnosis architecture for complex transportation systems (CTS) where train, vehicles and doors are considered as intelligent products. The generic holonic architecture is illustrated in Fig. 2. All elements are represented by Holons. The term Holon was proposed by [START_REF] Koestler | The ghost in the machine[END_REF]. A holon consists of an information processing part combined with a physical part. In this architecture, H1 represents a train. It has two children (vehicles H11 and H12). Besides, each vehicle also has two doors. The system S1 is decomposed into subsystem S11 and S12, S112 is a sub-system of S11. In this hierarchical structure, each level takes results from lower layer into account and report results to upper layer. This approach presents information exchanges on the same level and the recursion for different levels. The authors describe a modeling framework but does not detail the composition and decomposition mechanisms needed to manage a dynamic system. [START_REF] Sallez | On the activeness of physical internet containers[END_REF] introduce a recursive holarchy in logistics for Physical-Internet or (PI)-containers. These containers have some communicational and decisional capabilities. According to their size and function, three types are defined from the largest to the smallest: transport container, handling container and packaging container. Encapsulation and composition concepts are defined for the composition of containers. Indeed, containers of the same type can be composed as a "composite" PI-container. Small containers (handling and packaging) can be encapsulated in a bigger one. However, the compositional mechanisms are not formally specified. [START_REF] Esmaeili | A socially-based distributed self-organizing algorithm for holonic multi-agent systems: Case study in a task environment[END_REF] use a location information-based algorithm to construct the holonic multi-agent system (HMAS) for large scaled applications (LSA). In the proposed four-level architecture, a holon is selected as the head at each level except for the fourth layer. Holons can communicate with those at same level, or with the holons in higher or lower levels through the heads. This location information-based system can work well with static applications. However, the energy cost for reconstruction with this model is very high, it is hard to meet the high flexibility requirement of CM application.

A recursive agent meta model (RAMM) was proposed by [START_REF] Vo | An operational meta-model for handling multiple scales in agentbased simulations[END_REF] for multi-scale modeling (MSM). This metamodel has four main parts: species, agent, environment, and scheduler. Agents of different species can form a hierarchical structure. For example, a root agent located at level 0 consists of two agent species at level 1 (species 1 and species 2). At the lower level, species 3 is a micro species of species 1. The authors define the capture and release methods to add or remove micro species. However, the interactions between agents are not detailed. [START_REF] Occello | A multilevel generic multi-agent architecture for supervision of collective cyber-physical systems[END_REF] The hierarchical MAS-R structure is shown in Fig. 3. At the level 0, agents a0 i are "Elementary" agents which have a physical applicative part and a recursive agent part. The agent ai j represents the agent j on level i. These agents on high levels (i is greater than zero) are called "Partial" agent (no physical part). A "Complete Elementary" agent is composed by its elementary agent and all theirs directly related "Partial" agents (from abstraction levels). "Composed" agents are the gathering of few "Partial" agents from the same abstraction level. Besides, the authors also detailed the composition mechanism. Once composition conditions of an agent are met, it creates parent agent and invites others to join. This agent then adds accepted agents to the aggregated agent list. At the end, it sends "update" messages to all agents in the list to ensure that they have the same parent. The concept of recursion has been discussed or applied in all mentioned works. A summary is given in Table 1. We try to apply a dynamic recursive architecture to CM applications. Four characteristics are used to analyze these presented works: application domain, architecture, formalization of agents and aggregation mechanism. Application domain presents their applied filed. The architecture records whether the approach is fixed or dynamic. Whether they proposed an agent model or composition mechanisms are recorded in formalization of agent and aggregation mechanisms respectively. From Table 1, it can see that the needs of recursion are in different areas. The method given by (Le [START_REF] Mortellec | Embedded holonic fault diagnosis of complex transportation systems[END_REF] and [START_REF] Esmaeili | A socially-based distributed self-organizing algorithm for holonic multi-agent systems: Case study in a task environment[END_REF] are adapted to fixed systems, however, it is hard to apply for dynamic application like CM.

Although the aggregation notions (encapsulation and composition) and the interactions between agents were given by [START_REF] Sallez | On the activeness of physical internet containers[END_REF], the authors did not detail the aggregation mechanism. [START_REF] Vo | An operational meta-model for handling multiple scales in agentbased simulations[END_REF] provide a recursive agent metal-model, in which the aggregation mechanisms (capture and release) are discussed. However, the interactions between agents located at the same level or different levels are not detailed. MAS-R [START_REF] Occello | A multilevel generic multi-agent architecture for supervision of collective cyber-physical systems[END_REF] presents a recursive agent model, details interactions and aggregation mechanisms for a dynamic architecture. This method seems to be suitable for CM applications. In following sections, we try to apply this approach to CM context. The adapted agent model, recursive architecture and aggregation mechanisms are detailed for the CM application.

CM RECURSIVE AGENT MODEL

In this section, we adapt the MAS-R agent model for CM application where we add energy and some other information to the model. The CM agent model consists in five parts:

• The Knowledge part of agent (K): contains the agent's ID, energy level, recursion state, abstraction level, dynamic functions (it is used to calculate recursion state), energy model (it records energy consumption of its components). Cooperative Agent List (CAL) and other basic information;

• The Composition part (C): contains the composition function to construct a higher level;

• The Decomposition part (D): contains the decomposition function which asks parent agent to destruct the abstraction level;

• The Observer (O): observes the recursion state in D to call the recursive (composition or decomposition) functions and the communication functions to update neighbor information if necessary;

• The Recursive interaction part (RI): contains the interaction functions with other agents in the same group or with agents on other levels (parents or children), data processing and state updating.

The O part of agent periodically applies dynamic functions on the K part to analyze the recursion state. Once it finds that aggregation or decomposition conditions are met, it runs the composition/decomposition functions to create a higher abstract level or remove the existed abstract level (composition and decomposition mechanisms are detailed in next section).

Once the composition (or decomposition) process is done, the updating status process in the RI part will be performed to update the information stored in K.

Agents can be regrouped as a composed CM (see Fig. 4). The real node coupled with a node agent is an elementary material holon. These node agents belong to the same composed agent as shown by link 1 in Fig. 4. They can discuss with agents located at the same level by classical message as the interaction link 2. The abstraction link which present the link between an agent with its parent is illustrated as the link 3. Node agents use the transformation functions in RI part to send information to the composed agent.

CM RECURSIVE ARCHITECTURE

In different phases of the CM product lifecycle, one composed agents may cooperate with others. The recursive architecture is presented in this section to improve energy efficiency and extend service time.

The CM recursive architecture is shown in Fig. 5. In this example, there are six CNs in the physical part. Each node is represented by an elementary agent (Level 0) in the digital part. For an agent CMi j , i and j represent the abstraction level and ID respectively. Elementary agents (from CM0 1 to CM0 6 ) are aggregated into two CM agents (CM1 1 and CM1 2 ) on level 1 and CM2 1 on level 2. The agent located at the highest level has the global view of all CM products.

Fig. 5. The CM recursive architecture

There are three operators for aggregation, information exchange between agents in this recursive architecture:

• the PA operator (Process Agent): is responsible for the composition and decomposition processes.

• the PI operator (Process Interaction): is defined as the transformation of an interaction between two abstract levels, for example, the level n -1 and n.

• the PO operator (Process Organization): is defined as the transposition of relationship between levels. It can transfer the relationship from lower levels to higher levels.

When the observer part of an agent (for example CM0 1 in Fig. 5) finds that composition conditions are met, it applies the PA operator to begin the composition process. The sequence diagram of composition process for generating the parent agent CM1 1 is illustrated in Fig. 6. There are three elementary agents on level 0. Composition process exchanges can be explained in five steps:

1. The CM0 1 creates a parent agent CM1 1 and adds itself to the CAL of CM1 1 (message 1 in Fig. 6).

2. The request_compose messages are then sent to all cooperative agents in its CAL (message 2 and 3).

3. The invited agents process the composition request and reply an accept_compose message to accept invitation (or a refuse_compose message to refuse). CAL is also included in the accept_compose message (message 4 and 5).

4. The agent CM0 1 adds accepted agents to CAL and sends the energy information to its parent node (message 6). It then broadcasts update_compose messages to its child agents (message 7 and 8).

5. Once received the update_compose message, the accepted agents update their own CAL to ensure that they have the same parent.

Fig. 6. The sequence diagram of composition

Once CM1 1 is created and the structure is stable, it begins the process of energy model construction. This composed agent CM1 1 sends request_energy messages to all children to get energy model information.

During the CM1 1 lifecycle, its child agents check periodically whether they have relationship with another group (such as CM0 3 can cooperate with CM0 4 belonging to another CM). In that case, the agent CM0 3 requests the parent information of agent CM0 4 . This information is then transferred to its parent agent CM1 1 to update its CAL. By this way, the lower level is transferred to a higher level (the PO operator). The composition process will be repeated at CM1 1 or CM1 2 when composition conditions are met.

Once all composition process and energy model construction are completed as in Fig. 5. The agent at the highest level (the agent CM 2 1 ) can ask current energy level of child agents to estimate CM lifetime. The estimation process example is shown in Fig. 7.

The agent CM2 1 first send request_energy messages to its child agents (CM1 1 and CM1 2 ). These agents then check whether they have current energy level information. They can reply it directly if they have (such as the aggregation agent CM1 2 in this example). Otherwise, the agent (CM1 1 ) asks all child agents for their energy information thanks to the PI operator.

Once received all replies, the composed agent (CM1 1 ) can aggregate information and reply an aggregated ack_energy message (minimal or average energy level) to parent agent CM2 1 . At the end, the agent CM2 1 can then estimate CM lifetime based on the recent energy levels and energy consumption model.

Fig. 7. CM lifetime estimation process

If the agent CM 0 3 has run out its energy. Once the observer of a cooperative agent (for example CM0 2 ) notices the disappearance of this agent, the decomposition process begins as following:

1. The agent CM0 2 sends a request_reduction to its parent agent CM1 1 .

2. The agent CM1 1 processes the request and updates CAL, update_reduction messages are then sent to the agents in CAL.

3. The child agents update their own CAL based on the received update_reduction message.

To resume this section, the MAS-R adaptation for CM application is presented. The composition and decomposition mechanisms for constructing or destroying abstraction levels are detailed. In addition, the residual energy estimation process is also given for recursive architecture. In next section, the usefulness of this architecture is demonstrated on an applicative example.

APPLICATIVE SCENARIO

In some phases such as the construction for concrete element, two aggregated CM products may cooperate to improve energy performance as shown in Fig. 8. Each aggregated element (green and red agents) is composed by three elementary agents. The cooperation between the two independent DTs is discussed hereafter.

Once the observer part of an agent (CM0 3 in the green part) detected the red part, it sends request_coop to warm its parent. The parent agent (CM1 1 ) decides to cooperate or not and sends the ack_coop to all children. There are three cooperation possibilities (refuse, association and aggregation). The first case is no relation. Once the elementary agents receive a refuse message, they will maintain the original working state.

If the parent agent finds that the red part is related, it then replies an association message to all children. Only the connected agents (CM 0 3 in green and CM 0 4 in red) can exchange messages with each other. The PO operator can transfer this relation between the connect agents to its parent agents (CM1 1 in green and CM1 2 ). Therefore, these two parts can always know each other's states and cooperate to accomplish specific tasks if necessary.

Once the recursive states are met, all elementary agents will stop activity and wait orders from their parent node. The composed agent (CM1 1 ) then begins the composition process thanks to PA operator as in Fig. 6. Once the composition process and energy model construction are completed, the agent CM2 1 can then estimate CM lifetime as that in Fig. 7.

CONCLUSION AND PERSPECTIVES

This work presents the communicating material concept which is a special intelligent product with composition and decomposition abilities. Recent works applying MAS recursion are reviewed and compared with four characteristics to find an appropriate approach for CM application. The recursive architecture MAS-R proposed by [START_REF] Occello | A multilevel generic multi-agent architecture for supervision of collective cyber-physical systems[END_REF] is quite proper for the communicating material. The recursive architecture can meet the dynamic requirement of CM application. However, the intelligence distribution of agents on different levels are not discussed in MAS-R. We apply MAS-R approach and develop the missing parts needed for CM energy management. CM agent model and the recursive MAS architecture with composition/decomposition mechanisms are presented. For CM application, the agent located at highest level may estimate the residual lifetime and evaluate different strategies to provide an optimal solution to extend network lifetime. Applying the MAS-R method to CM energy management seems to be promising as illustrated on the CMs interaction scenario.

In the following works, the composition/decomposition conditions as well as the dynamic functions will be specified.

The interactions between agents on different levels will be developed. Besides, the intelligence at agents like data aggregation, residual energy estimation, selection of different data collection strategies will be specified. We will then apply the design method to the McBIM project, and develop a recursive digital twin to manage the energy and monitoring data of sensor nodes in real concretes.

Fig. 2 .

 2 Fig. 2. Holonic architectural model for embedded diagnosis (Le Mortellec et al., 2013) present SURFER, an embedded holonic fault diagnosis architecture for complex transportation systems (CTS) where train, vehicles and doors are considered as intelligent products. The generic holonic architecture is illustrated in Fig. 2. All elements are represented by Holons.The term Holon was proposed by[START_REF] Koestler | The ghost in the machine[END_REF]. A holon consists of an information processing part combined with a physical part. In this architecture, H1 represents a train. It has two children(vehicles H11 and H12). Besides, each vehicle also has two doors. The system S1 is decomposed into subsystem

  present a generic Recursive Multi-Agent System (MAS-R) for supervision and observation of largescale complex systems. A MAS-R agent has four parts: Knowledge (K), Observation (O), Composition (C) and Recursive Interaction (RI). Basic knowledge and abstraction level of agent are stored in K part. O is responsible for observing the state of composition. Composition tasks are performed in C. The RI part of an agent updates its own data, processes, and transfers data with other agents.

Fig. 3 .

 3 Fig. 3. The structure of MAS-R

Fig. 4 .

 4 Fig. 4. Composed CM agent

Fig. 8 .

 8 Fig. 8. Cooperation of two aggregated agents

Table 1 .

 1 Summary table of recursive MAS

	Reference	Application	domain	Architecture	Formalization	of agent	Aggregation	mechanisms
	SURFER (Mortellec et al., 2013)	CTS	Fixed	No	No
	PI-Container (Sallez et al., 2015)	Logistics	Dynamic	No	No
	HMAS (Esmaeili et al., 2017)	LSA	Fixed	No	Yes
	RAMM (Vo et al., 2012)	MSM	Dynamic	Yes	Yes
	MAS-R (Occello et al., 2019)	WSN	Dynamic	Yes	Yes

ACKNOWLEDGMENTS

Authors thank financial support from the French National Research Agency (ANR) under the McBIM project, grant number ANR-17-CE10-0014.