Issue 
Mechanics & Industry
Volume 24, 2023
History of matter: from its raw state to its end of life



Article Number  32  
Number of page(s)  27  
DOI  https://doi.org/10.1051/meca/2023027  
Published online  08 September 2023 
Sequential calibration of material constitutive model using mixedeffects calibration
^{1}
DMAS, ONERA, Université Paris Saclay,
F92322
Châtillon, France
^{2}
DTIS, ONERA, Université Paris Saclay,
F91123
Palaiseau, France
^{3}
ENSAI (Campus de Ker Lann),
51 Rue Blaise Pascal
35172
Bruz
^{4}
LIMOS (CNRS, Mines SaintÉtienne & UCA),
France
^{*} email: clement.laboulfie@onera.fr
Received:
10
February
2023
Accepted:
18
July
2023
Identifying model parameters is nowadays intrinsically linked with quantifying the associated uncertainties. While classical methods allow to handle some types of uncertainties such as experimental noise, they are not designed to take into account the variability between the different test specimens, significant in particular for composites materials. The estimation of the impact of this intrinsic variability on the material properties can be achieved using population approaches where this variability is modeled by a probability distribution (e.g., a multivariate Gaussian distribution). The objective is to calibrate this distribution (or equivalently its parameters for a parametric distribution). Among the estimation methods can be found mixedeffects models where the parameters that characterize each replication are decomposed between the population averaged behavior (called fixedeffects) and the impact of material variability (called randomeffects). Yet, when the number of model parameters or the computational time of a single run of the simulations increases (for multiaxial models for instance), the simultaneous, global identification of all the material parameters is difficult because of the number of unknown quantities to estimate and because of the required model evaluations. Furthermore, the parameters do not have the same influence on the material constitutive model depending for instance on the nature of the load (e.g., tension, compression). The method proposed in this paper enables to calibrate the model on multiple experiments. It decomposes the overall calibration problem into a sequence of calibrations, each subproblem allowing to calibrate the joint distribution of a subset of the model parameters. The calibration process is eased as the number as the number of unknown parameters is reduced compared to the full problem. The proposed calibration process is applied to an orthotropic elastic model with non linear longitudinal behavior, for a unidirectional composite ply made of carbon fibers and epoxy resin. The ability of the method to sequentially estimate the model parameters distribution is investigated. Its capability to ensure consistency throughout the calibration process is also discussed. Results show that the methodology allows to handle the calibration of complex material constitutive models in the mixedeffects framework.
Key words: Mixedeffects models / model calibration / composite material / sequential calibration
© C. Laboulfie et al., Published by EDP Sciences 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Following the increase of numerical facilities, virtual testing is now widely used in order to partially substitute numerical simulations to experimental campaigns because of time and costs considerations. These numerical simulations require the use of material constitutive models which describe different phenomena such as elasticity, viscoelasticity, plasticity, damage, etc. These models, in turns, involve several parameters that need to be calibrated. Calibrating model parameters consists in finding the parameters value that allow to best fit, in a specific mathematical sense that should be specified, the experimental responses of interest. This task is challenging because the model may not be able to fully reproduce the observations, the experimental data can be noisy and the material properties subject to inherent variability. Consequently, a faithful calibration of the model parameters requires to take into account the different sources of uncertainty in the identification process.
The properties of composite materials are known to be subject to a significant variability because of the complexity of the manufacturing process. To characterize this variability, international standards impose to perform test repetitions [1] to determine the model parameters distribution. In the classical approaches, model parameters are calibrated independently on several specimens by minimizing a fitting criterion [2] before inferring the model parameters distribution by maximizing the joint likelihood of the estimated parameters given a probabilistic model (e.g., multivariate Gaussian) [3]. This procedure is not fully satisfying because the fitting criterion (e.g., the sum of squares of the differences between observations and model predictions, the negative loglikelihood of the model given the experiments [3,4]) does not consider material variability in its definition. To overcome this difficulty, it has been proposed [5,6] to use population approaches [7,8] to characterize material variability. In this framework, the different repetitions are described by a specific set of model parameters known as the individual parameters, distributed according to a joint probability distribution (e.g., a multivariate Gaussian distribution). This distribution represents the impact of material variability on the model parameters. The mixedeffects models and the Hierarchical Bayesian models (HBMs) [9] belong to the populationbased approaches. The main difference between them is that mixedeffects models can be applied either in the frequentist or the Bayesian paradigm contrary to HBMs that are fully Bayesian, which usually involves costly Markov Chain Monte Carlo (MCMC) sampling. Given the computational costs of a single run of a contemporary mechanical models, the focus is set on the mixedeffects approach to be able to carry out the analyses more easily. The mixedeffects notion comes from the fact that there are “fixed” effects that are shared by the entire population of individuals (i.e., specimen) and “random” effects that are specific to each individual of the population. For instance, the Young’s modulus measured on a tension test individual can be considered as the combination of a reference value (the average value given by the manufacturer for a material batch) and of a deviation due to the variability of the production process. The parameters of this distribution can be determined by maximizing the appropriate likelihood of mixedeffects. This kind of approach has already been applied successfully to calibrate simple and uniaxial material constitutive models [10].
Because of their heterogeneous and anisotropic nature, the response of composite materials to external loads is a combination of several phenomena, mainly elasticity, viscosity, plasticity and damage [11]. Hence, various types of modelings have been developed, for instance damage [12], viscoelasticity [13,14], plasticity or viscoplasticity [15], coupled viscoelasticity and viscoplasticity [16], coupled viscoelasticity, viscoplasticity and damage [17,18], etc. The modeling can be carried out at the laminate or at the ply scale (see for instance damageplasticity models for a single ply [19–22]). In addition to the inplane features, complex constitutive models can also describe outofplane damage such as delamination [23]. Thus, the constitutive models to calibrate usually involves a large number of parameters, increasing the size of the search space and making the optimization of the likelihood more difficult. Furthermore, the computation of the likelihood of the mixedeffects is challenging as it is expressed as a product of multivariate integrals. Given the computational time of a single run of those models (especially for multiscale nonlinear models), the estimation of the likelihood requires to set up appropriate strategies (e.g., metamodeling [24]). Therefore, in order to decrease the complexity of the calibration in the presence of uncertainty, we will seek to take advantage of the decomposition of the model to transform the overall calibration problem into a sequence of appropriate smaller calibration subproblems, corresponding to the different test configurations (e.g. different load profiles and stacking sequences). In addition to decreasing the number of parameters to be identified, this alleviates the computational costs as it avoids to take into account all the phenomena at the same time.
To calibrate the material parameters, often, experimental tests of different natures are available, either in terms of loads (e.g., tension, compression or shear tests) or in terms of stacking sequences (e.g., unidirectional 0° laminates or quasiisotropic laminates). Nevertheless, the sensitivity of the model output with respect to the model parameters depends on the nature of the test [25] and all the model parameters characterizing a given phenomenon (elasticity for instance) cannot always be estimated with the same type of test. Consequently, in the decomposition of the calibration process, the choice of the relevant parameters must be performed carefully to avoid illposed problems.
To solve such problems, this paper introduces a sequential procedure to calibrate material constitutive models compliant with the mixedeffects framework. This procedure aims to calibrate the same orthotropic elastic model on multiple experiments with different stacking sequences. Prior to calibration, for each available individual, the parameters to which the model output is sensitive are determined from expert knowledge. This allows to define a sequence of calibration subproblems with a separation of the calibrated parameters. Note that the separation may not be strict as some parameters may be sensitive to several tests. For a given calibration subproblem (transformed into an optimization problem) attached to the appropriate experimental data, we only need to consider the joint distribution of the parameters relevant to this specific experiment. To ensure consistency between the different steps of the sequential calibration process, for a given step, the search space of the distribution parameters already estimated with previous subcalibrations is limited to a trust region. The distribution parameters of the model coefficients that are not sensitive are fixed during the resolution of this calibration subproblem.
This paper is organized as follows. Section 2 introduces the populationbased approaches and in particular mixedeffects models, with a focus on the treatment of multivariate data. In Section 3, the proposed sequential calibration process is presented in details. In order to illustrate the performance of the calibration strategy, the methodology is applied to the calibration of an orthotropic elastic model with nonlinear longitudinal behavior. In Section 4, a numerical study is carried out to investigate the performance of this methodology. Concluding remarks can be found in Section 5.
2 Populationbased approaches for calibration
2.1 Formalization
For a given a material constitutive model ℱ(·, θ), the calibration aims at finding a set of model parameters labeled θ ∈ ℝ^{d} (d standing for the number of parameters to be calibrated) which allows to mimic the experimental data. The random variable to θ is labeled Θ and θ stand for its outcome. The calibration problem in the presence of uncertainty can be formalized as follows [26,27]. Let us note the random vector of the output data Y and y = (y_{j})_{j∈〚1,N〛} its outcome, with N the number of observations, considering one individual. Classically, the following decomposition applies:
where ξ stands for the random vector of the errors which represents the experimental noise and the model bias. The outcomes of ξ can be different from one measure to the other and are labeled (ξ_{j})_{j∈〚1,N〛}. Equation (1) can be further detailed:
with t_{j} the j^{th} input measure, which is deterministic. A common practice of classical identification consists in minimizing, the sum of squares of the difference between the model prediction and the experimental data with respect to model parameters. This method is known as nonlinear least squares and the minimization is often conducted with the LevenbergMarquardt algorithm [28,29]. It should be noted that other criteria can be chosen to perform the calibration of the constitutive model [30]. To consolidate the identification of the model parameters, the sensitivity matrix that corresponds to the gradient of the constitutive model with respect to the model parameters can be used as proposed in [31]. Taking into account such information allows to better adapt the minimization algorithm to the local variations of the fitting criterion. Yet, such methods do not account for the intrinsic variability of the material properties. On the contrary, populationbased approaches explicitly model the variability between individuals within a population. Figure 1 illustrates this fundamental difference with respect to the classical identifications.
Populationbased approaches explicitly model the variability between individuals within a population. Figure 1 illustrates this fundamental difference with respect to the classical approaches summarized previously.
In the classical approaches (Fig. 1.a), individual variability is neglected and combined with other types of uncertainties, which allows to describe all the individuals by the same value of model parameters. In the populationbased approaches (Fig. 1.b), each sample (experiment replication) is assigned a specific parameter vector value called the individual parameters. It is assumed that there exists an underlying probability distribution f_{Θ} whose outcomes are the individual parameters (θ_{i})_{i∈〚1,n〛} with n the number of individuals (corresponding to the number of experiment repetitions) [7,8]. Both the underlying probability distribution of the model parameters and the value of the individual parameters have to be determined. In addition, it remains possible to take into account other sources of uncertainty, such as measurement noise. Thus, populationbased approaches are particularly well suited to situations in which individual variability cannot be neglected with respect to other sources of uncertainty
The populationbased approaches framework [7,8] assumes that there exists a probability distribution f_{Θ} whose outcomes are the individual parameters:
Both the probability distribution f_{Θ} and its realizations θ_{i} are unknown and the aim is to determine them. In addition, if f_{Θ} is parametric (a Gaussian distribution for instance), let Π be its parameters and f_{Θ} = f_{ΘΠ} Identifying f_{ΘΠ} is tantamount to determining Π. Given Π and θ_{i} ~ f_{ΘΠ} ∀_{i} ∈ _{〚1,n〛}, the model output y_{i} can be written as
Without any other hypothesis, the outcomes of ξ_{i} (labeled ) are different for each individual and for each observation, with N_{i} the number of observations of the i th individual . The global mixedeffects models for the j^{th} output measure of the i^{th} individual y_{ij} reads as
For the sake of simplicity and to ease the numerical derivations, f_{ΘΠ} is chosen in this paper to be a multivariate Gaussian distribution of dimension d (remember that d corresponds to the number of parameters to calibrate) but the principles of the approach remain valid for other distributions:
with µ ∈ ℝ^{d} the mean vector and Σ∈ ℳ_{d}(ℝ) the covariance matrix The second hypothesis is that for each individual and each measure, the discrepancy term is a Gaussian white noise (no bias, no correlation):
with ω_{ij} the standarddeviation of the noise of the j^{th} output measure of the i^{th} individual Furthermore, the noise is supposed to be homoscedastic, i.e. ω_{ij} = ω_{i} ∀ j ∈ 〚1, N_{i}〛 [32] This is a simplifying assumption but helps to limit the number of parameters to estimate Finally, the vector of parameters to be calibrated is denoted Ψ
with . The populationbased approaches seek Ψ and provide an estimate of the individual parameters (θ_{i})_{i∈[1,n]} as a byproduct. The complete hierarchical modeling, that combines both the population level (the repetition of experiments) and the individual level (the specimens themselves), is summarized in Figure 2.
Figure 2 highlights the hierarchical link between the population level (that is to say the repetitions of specimens) and the individual level for which standard calibration methods operate. In fact, population approaches appear as a generalization of the standard calibration methods with the additional feature that they statistically account for the variability that arises from repetitions of specimens.
To identify Ψ, that is to say the population and residual model, different methods can be set up in either frequentist or Bayesian fashions. Among the frequentist methods can be found mixedeffects models [7,8,34]. The notion of mixedeffects comes from the decomposition of the individual parameters θ_{i} through fixedeffects β and random effects b_{i} as
where A_{i}, B_{i} are deterministic design matrices that belong to ℳ_{d} (ℝ) with d the number of model parameters. Design matrices other than I_{d}, the identity matrix of ℳ_{d} (ℝ), can be used for example by removing them through null coefficients in the matrices if some of the model parameters are deterministic or do not express themselves on all experiments [32] (e.g., compressive properties on the tension behavior). Generally, there is no direct interpretation between the fixedeffects, random effects and Π, but for a Gaussian f_{ΘΠ}, the fixedeffects β are equal to the population parameters average µ and the randomeffects comply with ,i ∈ 〚1, n〛. Furthermore, provided that the modeling proposed by the population exhibit a hierarchical structure, hierarchicalbased techniques such as Hierarchical Bayesian models (HBMs) [6,9] can be implemented. HBMs are a fully Bayesian approach that estimates the population parameters with a Bayesian paradigm. As in all Bayesian problems, it requires the definition of a prior on the population parameters Π (to which may be added Ω, thus getting Ψ) that allows to derive the expression of the posterior distribution of the population parameters given the available data thanks to the Bayes rule. HBMs techniques may lead to two challenges. The first is the choice of the prior distribution of the population parameters Π, in particular for the covariance matrix Σ. The second resides in the estimation of the posterior distribution from Bayes rule. Except for specific simple linear models to calibrate, the posterior distribution cannot be expressed analytically (in particular because of the normalization constant [35]) and should be sampled with, for instance, MarkovChainMonteCarlo techniques (e.g., the MetropolisHastings algorithm [36] or the Gibbs sampler [37]). Moreover, because of the hierarchical nature of HBMs, it is necessary to use nested MCMC. This makes the sampling procedure harder than when no population modeling is considered, for which it is already complicated to achieve. Given that MCMC techniques are computationally intensive, this may not be applicable in applications involving computationally costly simulation codes. On the contrary, mixedeffects models provide more flexibility as they can be applied in either the frequentist and Bayesian frameworks. Note that in the Bayesian framework, both methods turn out to be equivalent. They will be here the chosen technique, mainly because of computational concerns in order to be able to perform the necessary computations. With mixedeffects, the population parameters are most of the time estimated by likelihood maximization [7,8,32], that will be derived in the next section for first models with univariate output and then for multivariate output. Note that all the derivations presented in the following assume that the design matrices are set equal to the identity matrix for the sake of clarity. If it is not the case, the relevant linear change of variables should be applied.
Fig. 1 Difference between the populationbased and the classical approaches. On the left, a), usual methods with a single parameter vector for all individuals. On the right, b), populationbased approach in the mixedmodels effects framework with f_{Θ} the distribution. 
Fig. 2 Different levels of modeling in the population approach (adapted from [8] and [33]). The dashed lines indicate random sampling and the ellipsis the distributions. 
2.2 Likelihood of the mixedeffects
The calibration is often achieved by maximizing the likelihood of Ψ, ℒ(Ψ) := f(y_{1},…, y_{n}Ψ) (where f is a generic letter for probability density functions or PDFs), even if other methods can be found to estimate Ψ [38]. The stepbystep derivation of the likelihood function can be found in [10]. Under the assumption of independent individuals (corresponding to independent experiment replications, for instance with different samples), it reads as the product of all the individual likelihoods ℒ_{i}(Ψ) := f (y_{i}Ψ),
Because the θ_{i} is are not observed, the likelihood of the i^{th} individual ℒ_{i}(Ψ) is the integral of the marginal likelihood (i.e., the density of the output data y_{i} given individual parameters θ_{i} and parameters Ψ) f (y_{i}θ_{i}, Π, ω_{i}) with respect to all possible θ_{i} over ℝ^{d}:
In Equation (10), f_{ΘΠ}(θ_{i}Π) refers to the PDF of the model parameters distribution computed for the individual parameters. Note that, as a multivariate integrals, the individual likelihoods are challenging to compute. In particular, it should be ensured that none of them collapse to 0, as it would otherwise flaw the estimation of the complete likelihood function. The maximum likelihood estimator (MLE) is the result of the following maximization problem over the set of all the possible parameters Γ_{Ψ} :
In practice, the logarithm of the likelihood is computed to ease the numerical optimization [39].
2.3 Mixedeffects for multivariate models
The calibration of mixedeffects for univariate models has been extensively discussed [7,8,32]. This happens when only one measure is available on each individual. Yet, sometimes, one has access to several measures on the same individual (with longitudinal and transverse gauges or digital image correlation for instance). Then, the model to be calibrated becomes multivariate in the sense that the output is no longer a scalar but rather a vector, requiring to adapt the mixedeffects equations [40].
Let’s note m_{i} the number of output measures of the ith individual and k ∈ 〚1, m_{i}〛 the corresponding index. In the following, all individuals are supposed to have the same number of output measures, so m_{i} = m ∀i ∈ 〚1, n〛. Moreover, all measures are supposed to have the same number of observation points N_{i}. The main difference with Section 2.1 is that the model output of each observation point is no longer a scalar y_{ij} but a vector y_{ij}·∈ ℝ^{m}. The same occurs for the discrepancy term, now labeled ξ_{ij}.
The equations governing the mixedeffects are [7,8,32]:
with the zero vector of ℝ^{m} and Ω_{ij}· the covariance matrix of the discrepancy term which describes both the errors of the different sensors and the possible correlations between the errors on two different sensors. The observations are supposed serially uncorrelated as in earlier studies [41], meaning that the measurement errors at two different observation times are independent:
Furthermore, for all individuals, the noise is supposed to be homoscedastic, that is to say: Ω_{ij} = Ω _{i} ∀ j ∈ 〚1, N_{i}〛Note that the only difference between equations (5) and (12a) is that all terms are now vectors rather than scalars and that equation (12c) is the analogous of equation (7) for vectors instead of scalars
Let us define for all individuals Y_{i}, the matrix of the output data of the ith individual with Y_{i} := () . The rows of Y_{i} correspond to , that is to say the N_{i} observations of the kth measure of the ith individual. Similarly, it is possible to define the matrix of the errors as ξ_{i} := () . The rows of the ξ_{i} matrix correspond to and can be understood as the N_{i} measurement errors of the kth measure for the ith individual. Let us define the flattened array by row of the matrix of the output data Y_{i} and the same for the matrix of the discrepancy term ξ_{i}. Then, and belong to . This allows to transform equation (12) into
with where ⊗ denotes the Kronecker product, the identity matrix of (ℝ) and ℝ the zero vector of ℝ. Hence, the set of the model parameters to be calibrated, labeled Ψ is defined as:
Finally, the likelihood of the mixedeffects becomes . Equation (10) remains the same except for y_{i} that becomes .
2.4 Likelihood estimation with Laplace approximation
The evaluation of the likelihood function requires to compute the individual likelihoods, which implies to estimate the multidimensional integral of equation (10). A fundamental method to compute the integral is the MonteCarlo method [42]. However, such an approach requires too many model evaluations to keep the computational time reasonable and may lead to numerical noise in the estimation of the multivariate integral resulting in challenge in the likelihood maximization. An alternative approach, named Laplace approximation, can be used [32]. The Laplace approximation [32] applies to integrals of the type
where h(·) is a function which must satisfy the following constraints:
h(·) admits a global minimum x_{0} that belongs to the integration interval,
h(·) is a twicedifferentiable function,
its Hessian matrix ℋ(h) computed at x = x_{0} is a symmetric definite positive matrix.
The main idea is to assume that only regions close to x_{0} significantly contribute to the integral, leading to the following Laplace approximation of A:
Given the modeling choices (Eq. (14)), the individual likelihoods read as:
Here, Δ is the transpose of the result of the Cholesky decomposition of Σ^{−1} (so Σ^{−1} = Δ^{T} Δ), Λ_{i} is the transpose of the Cholesky decomposition of and the function ɡ_{i}(·) is defined by:
The Laplace method is applied in two steps:
Search for the individual parameters (or rather the deviations), , minimizing ɡ_{i}(·) equation (17),
Computation of the Laplace approximation with equation (15).
To finalize the approximation of the likelihood, it remains to estimate the Hessian matrix of ɡ_{i}(·) at [32]:
In practice, can be neglected if the model ℱ( · , · ) is close enough to the experiment [43]. The term is evaluated using a finite difference scheme. As a result, the negative loglikelihood is finally expressed as
2.5 Individual parameters estimation
The mixedeffects model also allows to infer the individual parameters. With the Laplace method, individual parameters are byproducts of the likelihood function calculation given Ψ. Once Ψ is chosen, the random effects b_{i} are estimated by looking for the dominating contribution in the individual likelihood ℒ_{i} (Eq. (16)) through the minimization of function ɡ_{i}(·) (Eq. (17))
with Γ_{b} the domain definition for the individual deviations. The individual deviations are then simply computed as θ_{i} = µ + b_{i}.
The determination of the individual parameters is a key step to ensure a proper likelihood maximization and should not be neglected. To explain its importance, it should be underlined that mixedeffects aim to tune the population parameters (Π) in order to find the best possible agreement between the assumed model parameters distribution (f_{ΘΠ}) and the empirical distribution of the individual parameters. Thus, the observation data should provide enough information to accurately estimate the individual parameters. In particular, the material constitutive model should be appropriately chosen with respect to the observations and not present imperfections (e.g., not overparametrized). Otherwise, if the individual parameters cannot be learnt with sufficient precision, the regularization term (f_{ΘΠ}(θ_{i}Π) in equation (10) or ‖Δb_{i}‖^{2} in function ɡ_{i}(·) in equation (17)) will coerce them to the mean value and thus flaw the estimation of the standarddeviations and correlations parameters. In fact, generally, it appears that to ensure a proper estimation of a given marginal of f_{ΘΠ}, the corresponding parameter should be properly estimated over the different individuals, and the same occurs for the parameters involved in a dependency relation.
2.6 Likelihood maximization estimate
Now that the expression of the likelihood for a mixedeffects model has been established, the model parameters can be identified by minimizing the negative loglikelihood,
with Γ_{Ψ} the research space for Ψ. Remember that, as computing the likelihood is numerically challenging, the Laplace approximation of the likelihood is used instead in the optimization process. To solve the minimization problem (Eq. (21)), the first possibility would be to rely on gradientbased algorithms, which need either the gradient or an approximation of the gradient of the objective function. Here, the gradient of — ln(ℒ(Ψ)) expressed with equation (19) with respect to Ψ is difficult to compute analytically. Indeed, the differentiation with respect to the population parameters of function ɡ_{i}(·) requires to differentiate the individual deviations b_{i} defined in equation (20) with respect to the population parameters (i.e., µ and Σ via Δ). Moreover, the loglikelihood function often presents numerous local minima that also makes the optimization task by gradientbased algorithms complicated.
A second possibility is to use gradientfree optimization algorithms, and in particular evolutionary algorithms. Though they require many likelihood estimations, they exhibit a greater ability to reach the global maximum of the likelihood, , despite the possible presence of local minima. The algorithm chosen here is the Covariance Matrix Adapatation Evolution Strategy (CMAES) algorithm [44]. Given that a multivariate Gaussian distribution of dimension d ( d stands for the number of parameters to be calibrated) is parametrized by at least d means and d variance parameters, and that mixedeffects models allow to estimate the variance of the discrepancy term, there is always at least 2d +1 parameters to estimate in equation (21) and thus, the dimension of Ψ labeled Ƥ complies with Ƥ ⩾ 2d + 1. The complete principle of mixedeffects calibration that starts with the population modeling and ends with the individual estimation step is depicted in Figure 3.
3 Sequential calibration process
Advanced mechanical constitutive models characterize different aspects of the behavior of the material they describe, among which elasticity, viscoelasticity behavior, damage, etc. Each of these phenomena may be described by its own set of model parameters, which finally makes the number of parameters to be calibrated large ( 21 for the 3 — d anisotropic elastic model for instance [25]). Consequently, a common practice [45,46] in classical calibration consists in dividing the full calibration problem into calibration subproblems that allow to identify different subsets of the complete parameters vector [47,48]. In these works, the different subproblems are associated with a specific part of the model parameters. The methods consist in solving a sequence of calibration problems taking advantage in each step of the knowledge acquired so far. To be more precise, for a given step, in a frequentist framework, the coefficients calibrated in the previous steps are either taken as fixed inputs or severely bounded. In a Bayesian fashion, the identified posterior probability density of one step becomes a prior density for the next stages [49,50].
The sequential approach requires to define a sequence of calibration subproblems, which must be performed carefully. Indeed, as each step depends on the previous ones, an error at one step may flaw the whole calibration process. In some cases, the decomposition in subproblems can be rather easy to perform with the help of expert knowledge [47,48]. For example, Jung et al. [47] choose to separate the model parameters between those that describe the electrical and mechanical behaviors of an energy harvester. Yet, in some cases, such a decomposition is not available. Consequently, to avoid an arbitrary choice, another option consists in relying on a statistical criterion. For instance, given the available experiences, it can be decided to select the parameters to which the model output is sensitive following the analysis of the rank of the Fisher Information Matrix (FIM, the Hessian matrix of the loglikelihood) [51]. This takes advantage of the fact that, for classical calibration problems, assuming a Gaussian discrepancy ξ in equation (1) makes the negative loglikelihood equivalent to least squares, allowing to compute the Fisher Information Matrix from the derivative of the material constitutive model with respect to the model parameters.
The principal limitation of the previous studies is that the statistical modeling developed do not consider the inherent variability of the individuals on which the experiments are performed (and thus do not require the population modeling). On the contrary, in this paper, the calibration is carried out in the context of the population approaches, which complicates the modeling and the estimation procedure. This paper aims to adapt the sequential calibration strategy to the mixedeffects models. Remember that the objective of the mixedeffects framework is to characterize the material variability observed on a population of individuals modeled by the model parameters distribution. According to Sklar’s theorem [52], there exists a unique decomposition of any probability density function (PDF) between marginals and copula. Applied to the model parameters distribution f_{ΘΠ}, it gives :
with Π_{cop} the copula parameters, Π_{l}; the parameters of the l^{th} marginal of distribution f_{ΘΠ} (so the marginal marginal distribution of θ_{il} (i, l) ∈ 〚1,n〛 × 〚1,d〛) with probability density function f_{l}(·) and cumulative density function F_{l}; (·) and c(·), the density function of the copula (Gaussian in this paper). In a case where all parameters are independent (the copula density is constant equal to one), following equation (22), it is possible to calibrate separately the marginals of the relevant parameters for each type of experiments. Note that in the presented procedure, the parameters whose dependencies cannot be estimated (because they do not act simultaneously on the same test) are set as independent. Furthermore, the Gaussian copula (chosen in this paper) is compatible with a sequential scheme. Indeed, remember a Gaussian copula is parameterized by a correlation matrix which is symmetric positivedefinite. Sylvester’s criterion [53] states that if the determinant of the submatrices extracted from the full correlation matrix by deleting a row and the symmetric column is strictly positive, then the full correlation matrix remains positivedefinite. This condition is fullfilled if these submatrices are positivedefinite. Then, it is possible to sequentially estimate a Gaussian copula, allowing, along with the estimation of the corresponding marginals, a sequential calibration of f_{ΘΠ}. The sequential calibration process consists in defining a sequence of subsets of indices S_{q} ⊂ 〚1, d〛 corresponding to the parameters to which the model parameters is sensitive (determined for instance from the FIM) defining the marginals and correlations parameters Π_{q} to be calibrated at step q, leading to the definition of a specific likelihood function ℒ^{q}. Note that these subsets are not necessarily disjoint. In practice, it is not likely that all parameters will be sensitive enough to be propely estimated with standard experiments. Thus, it can be considered that in each step, the joint distribution of a part of the model parameters will be estimated, reducing the number of population parameters to calibrate at each step. Once these parameters are estimated, a trust region ℛ_{q} is defined. If some of the parameters of the q^{th} step are reestimated, they are limited to their trust region ℛ_{q} (this choice is discussed further below). This allows to ease the overall likelihood optimization as it focuses on a specific regions of the parameter space, but also enables to make sure that the calibrated distribution remains consistent with the data used in the previous stages. The definition of the sequence of calibration relies on the observation that the identification of a given marginal requires the model output to be sensitive enough to the corresponding parameters. This allows to design relevant experiments to estimate properly the PDF of interest. Figure 4 sums up the different stages of the method.
The choice of the trust region to which the parameters are limited is a key feature of the process. In theory, it should illustrate the trust we have in the estimated parameters and allow to keep some consistency between the population parameters estimates throughout the different steps. From a numerical point of view, the trust regions ease the negative loglikelihood minimization as the search space in which the population parameters are refined at step q should be close to the Maximum likelihood Estimate (MLE) of step q — 1. Nevertheless, as likelihood functions may present local optima, one drawback of this method is that the search space can be reduced around a local optimum from which it can be hard to escape, making the minimization of the negative loglikelihood more prone to local convergence than in usual cases.
Several options can be explored to choose the trust region. The easiest one is an hypercube centered on the population parameters that have already been estimated. For instance, from step q — 1 to q, it would correspond to [αΠ_{l,1}, βΠ_{l,1}] ×·… × [αΠ_{l,q–1}, βΠ_{l,q–1}] for the parameters l ∈ 〚1, Sq–1〛, with α and ß indicating the size of the interval, s_{q–1} the number of common parameters from step q – 1 to q and the corresponding coordinates. Note that the size of the intervals a and β can also depend on the step q and be different for each coordinate. Most of the time, α will be strictly smaller than β, meaning that the can be updated. Still, they can be set equal, which is tantamount to fixing the population parameters to their current estimate. If this trust regions management is easy to implement in practice, its main drawback is that it is arbitrary as it defines a range of variation which does not necessarily reflect the level of uncertainty of the population estimates.
A more appropriate choice for tuning this trust region would be a confidence interval, given either by the asymptotic theory or by bootstrap [8,32]. These proposals are relevant as they quantify the uncertainty on the population parameters that comes from the lack of data in addition to other sources of errors. Thus, they seem to be the most appropriate candidates for the trust region. However, these confidence interval techniques suffer from several limitations that make them difficult to use in practice. Indeed, aside from their respective issues discussed later, both of them are based on asymptotic assumptions, that is to say when the number of individuals increases to infinity. Of course, the asymptotic regime is not achievable because for times and costs considerations, and in practice, few (below 100) individuals are available. Thus, the intervals that reflect the uncertainty could be wider than expected, limiting the interest of a sequential approach based on a limitation of the search space. Furthermore, both techniques suffer from their own drawbacks. For instance, the asymptotic theory is based on the Fisher Information Matrix that is the expectation with respect to the population parameters of the secondorder derivative of the loglikelihood. Its evaluation is difficult and demands appropriate estimation techniques [54,55]. The observed FIM, in fact the secondorder derivative of the loglikelihood with respect to Ψ is easier to estimate but less accurate, which is a key feature given that the trust region is described by a covariance matrix computed as the inverse of the FIM, leading to possible numerical issues. Bootstrap does not suffer form this limitation as, from the bootstrap samples, confidence intervals can be estimated straightforwardly, without the need of any approximation. Yet, the accuracy of the estimations requires as much samples as possible, which is synonymous a heavy computational burden because the each sample is the result of a likelihood maximization. Finally, note that both of these methods remain arbitrary to some extent as it is necessary to define the levels of confidence that defines the size of the trust region.
Consequently, for this first implementation, the chosen trust region is the first proposed, that is to say defined from arbitrary intervals centered on the population estimates with α = 0.8 and ß =1.2 (the boundaries correspond to ±20% around the calibrated values). Such a range should allow to make a compromise between exploration and intensification for the maximization of the likelihood. It was determined here by trial and error. Such settings should be considered as hyperparameters of the method and ought to be carefully chosen.
Fig. 3 Recap of the mixedeffects calibration procedure. Modeling stages are highlighted in red, those that concern the calibration itself are highlighted in blue. 
Fig. 4 Different stages of the calibration process. Once the sequence of calibration problems to resolve is defined, the calibration process starts. At each iteration, the search space is defined depending on whether parameters have already been calibrated to preserve consistency between the different steps. Calibration stages are highlighted in red and those concerning the determination of the search space are highlighted in blue. At each iteration, u_{q} refers to the integer labeling the subgroup of parameters that will be reestimated in the current iteration. N_{G} refers to the number of calibration subproblems. The process stops when all calibration subproblems are solved. 
4 Applications
This section presents the application of the considered sequential calibration to the identification of an orthotropic non linear elastic model. After a brief description of the model in Section 4.1, the applied sequential strategy is exposed in Section 4.2 before describing the assessment protocol of the proposed method in Section 4.3. The objective is to validate the model on synthetic data first to study the consequences of such a methodology in Section 4.5.
Fig. 5 Illustration of the role of the model parameters. The solid line indicates the model output, the blue dashed lines the asymptotes and the green dashed line the tangent at the origin. Exponent T stands for tension load and C for compressive load. 
4.1 Presenting the elastic model
This section aims to present the orthotropic elastic behavior of an elementary unidirectional (UD) ply in plane stress [25]. Let us note σ := (σ_{11}, σ_{22}, σ_{12}) the Cauchy stress (in MPa) and ε := (ε_{11}, ε_{22}, 2ε_{12}) the observed strain dimensionless within the material frame characterized by the orthogonal frame (1, 2, 3). Axis 1 represents the fiber direction, axis 2 represents the matrix direction. The orthotropic elastic model can be written (using the Voigt notation [25]) in the material axes:
with the compliance tensor. is a symmetric tensor so S_{12} = S_{21}. For practical interpretations, the coefficients of the compliance tensor can be related to the elastic moduli, that is to say the longitudinal modulus E_{11}, the transverse modulus E_{22}, the shear modulus G_{12} and the Poisson’s ratio ν_{12} [25].
The model presented above is elastic linear and can be used either in tension and compression. Yet, it has been observed that for some materials, in tension, 0° degree laminate tends to stiffen, while in compression, 0° degree laminate tend to soften [56]. This is the consequence of the fiber behavior. Indeed, when produced, the fibers are not always properly aligned and when loaded in tension, fibers get aligned increasing the stiffness, while in compression, the initial defaults are increased and tend to soften the material. Thus, the apparent modulus increases (non linearly) throughout the tension experiments and conversely on the compression experiments. Thus, the longitudinal compliance S_{11} is changed and expresses as follows:
with the elastic linear compliance (module at the origin), the asymptotic compliance and the abscissa of the intercept between the abscissa axis of the asymptotic curve. The symbol L indicates whether tension (T) or compression (C) is considered. The transverse compliance should also be updated and refers to . Figure 5 details the role of the different parameters.
As the compliance of the material decreases in tension and increases in compression, the parameters should comply with and . With this parameterization, the model is not identifiable. In fact, it is overparameterized, i.e. there are several combinations of both and that give the same output as shown by Germain [57]. Thus, fixed to a usual nominal value, here 0.005 [57]. Finally, there are 6 parameters to be calibrated , , , S_{22},v_{12} and S_{66}.
In theory, the compliance matrix should be positivedefinite [25] To enforce this property on the individual compliance tensor, two requirements have to be met: the strict positivity of S_{66} and the positivedefiniteness of the submatrix containing all other compliances (this is in fact the application of the Sylvester criterion [53]). The first point is easy to check but the second one deserves more attention. Indeed, it can only be guaranteed if all the involved compliances can be calibrated on the same specimen. In that case, it is possible to enforce such property by implementing constraints between the compliances or by decomposing the compliance matrix with a relevant factorization (e.g., Cholesky, spherical decomposition). However, none of these solutions is fully satisfactory as they either increase the difficulty of the auxiliary minimizations or use parameters that are difficult to interpret from a physical point of view. In any other case, for most calibration methods, it is not possible to guarantee positivedefinitness of the compliance matrix in a simple way, because, on a given specimen, all parameters of interest are not available. A more convenient way to proceed, chosen in this paper, consists in calibrating the elastic behavior without accounting to this constraint and to discard the realizations of matrices that are not definite positive in a subsequent.
4.2 Application of the proposed sequential strategy
Given the modeling choices, the model parameters distribution is described by 27 parameters, namely 6 mean parameters, 6 variance parameters and 15 covariance parameters. The distribution parameters can be summed up into a mean vector and covariance matrix
where the * stands for the symmetric coefficient of the covariance matrix.
To identify the elastic constants, a usual practice is to rely on tests on unidirectional 0°, 90° and ±45° laminates [25,58,59]. Indeed, from a physical point of view, , , and v_{12} characterize the fiber behavior and are almost exclusively monitor the 0° laminate. Similarly, the linear part of the longitudinal response of a tension/compression test on a 90° laminate is mainly controlled by S_{22}, the same applies for the linear part of the difference between the longitudinal and transverse responses for the ±45° laminate, mostly driven by S_{66} [25,58,59]. Remember that, with mixedeffects, to estimate properly a given marginal, the corresponding parameters should be accurately estimated over the different identified individuals, and the same occurs for the parameters involved in statistical dependency relation. Thus, from the analyses made earlier, the correlations between , , , v_{12} and S_{22} cannot be estimated with the usual elementary experiments, and the same applies to those between , , v_{12} and S_{66} or between S_{22} and See. In fact, from these elementary experiments, only the joint distribution of , , , v_{12} and the marginals of S_{66} and S_{22} can be estimated, which gives three distributions to calibrate (in addition to the residual models). This makes 18 population parameters (14 for the joint distribution of , , v_{12} and 2 for both marginals of S_{22} and See) to estimate.
From this assessment, two solutions arise. On the one hand, it is possible to calibrate all the experiments together by making sure that for each individual, only the relevant individual parameters are calibrated. This can be done through equation (8) by choosing relevant design matrices A_{i} and B_{i}. For instance, with the mean vector defined earlier, one possible choice would be A_{i} = B_{i} equal to the identity matrix of dimension 6 to which the two last rows would be deleted (and A_{i} would have thus 4 rows with 6 columns) for the UD0° experiments, A_{i} = B_{i} = (0, 0, 0, 0, 1, 0)^{⊤} for the 90° experiments and A_{i} = B_{i} = (0, 0, 0, 0, 0, 1)^{⊤} for the ±45° experiments. Yet, to ease the calibration process, as the parameters that describe the behavior of the 0° laminates do not express on the behavior of the 90° laminates and similarly between the other laminates, it is possible to identify the joint distribution of , , , v_{12} and the two other marginals independently from the others (which thus do not require to use design matrix). If the estimation of these distributions is carried out separately, this can be considered as a first application of the sequential strategy in an extreme case, given that between the three cases, there are not shared parameters, and the calibration of the three distributions can thus be achieved independently. The calibration of the marginals of S_{66} and S_{22} are carried out from experiments with univariate observations that exhibit a linear behavior. Thus, these cases only refer to univariate linear mixedeffects models, that is to say the simplest case that can be found and extensively discussed in [7,8,32]. Consequently, the study is here focused on the calibration of the joint distribution of , , and v_{12}.
One first remark that can be made is that only expresses when compressive loads are applied to the sample, and only expresses when tension loads are applied. Thus, from the elementary tests, there is no configuration in which they are activated simultaneously, and thus the correlation between them cannot be estimated. Therefore, contrary to the above comments, the joint distribution of , , and ν_{12} is parameterized by 13 and not 14 population parameters. and v_{12} are activated both in tension and compression. Note that v_{12} only expresses on the transverse responses, either in tension and compression. Here, there are shared parameters for two types of experiments, and v_{12}. Design matrices apply here. For instance, for the tension experiments, matrices A_{i} and B_{i} would correspond to the identity matrix of dimension 4 to which is deleted the third row (so both have 3 rows and 4 columns) and for the compression experiments, they could be taken as the identity matrix of dimension 4 to which is deleted the third row (so both have 3 rows and 4 columns). The interest of this test case is that it also allows to illustrate the sequential strategy. In particular, it enables to compare the accuracy of the calibration from the sequential strategy to the solution that implements design matrices. For 0° UD laminate, let consider a tension load σ = (σ_{11}, 0, 0) with σ_{11} > 0. Then, following equation (23), the model output is determined by , and v_{12}. It is thus possible to calibrate the parameters highlighted in red in equation (26):
This allows to estimate a first trust region: ℛ_{1}. Let’s now consider a compressive load σ = (σ_{11}, 0, 0) with σ_{11} < 0. Then, following equation (23), the model output is determined by , and v_{12}. The application of previous remarks shows that it is possible to calibrate the parameters highlighted in blue in equation (26):
Note that , , , , v_{12}) and 𝕍(v_{12}) are limited to their trust region ℛ_{1} in the second step. All the distribution parameters of the model are now estimated through these two steps. From this point, it could be decided to restart the calibration with the tension experiments to refine the estimation of the joint distribution of , and v_{12} by limiting the parameters of the joint distribution of ,v_{12} to trust region ℛ_{2} (determined from parameters calibrated in the second phase) and the others to another that would be based on the calibration in the first phase. This would be repeated upon reaching convergence (that should be defined), as in a fixed point problem. This first difficulty is to define convergence (at the level of the population parameters, of the calibrated distributions, etc.). This is of prime importance because a stationary behavior could only be met for some convergence criteria, and nothing guarantees that only few iterations would be enough, increasing the computational costs. Still, note that this paper is focused on the sequential strategy, leaving the description and the generalization of the sequential strategy with a fixedpoint approach to future investigations. The application of the sequential strategy has only been presented in one way, starting with the tension experiments and ending with the compression experiments. Yet, there is no particular justification to such a choice. Consequently, it is here necessary to carry out the sequential strategy in both ways, that is to say starting with the tension experiments and ending with the compression experiments and then starting with the compression experiments and ending with the tension experiments.
One of the main advantage of the sequential approach is to reduce the maximum number of parameters to estimate compared to the classical method using design matrices. For instance, the calibration of the nonlinear elasticity (that is to say the joint distribution of , , , v_{12}, S_{22} and S_{66} parameterized by 17 parameters with a dependence structure and 12 without) is divided with a sequential scheme into subproblems with at most 9 parameters with a dependence structure and 6 without. This reduction can illustrated further if another component is added, for instance the viscoelasticity behavior as in the ONERA Progressive Failure Model (OPFM) [56]. This component adds 6 new parameters and can be calibrated from both 90° and ±45° laminates (in which also express the linear compliances S_{22} and See). With the same modeling assumptions (Gaussian marginals and copula), the joint distribution is parameterized with 42 parameters with a dependence structure (all correlations but two can be estimated) and 16 without. Implementing a sequential scheme allows to substitute to these problems two calibration of distributions with respectively 27 parameters with a dependence structure and 12 without. Note that this analysis does not take into account the residual models which is discussed later. Table 1 recaps the number of population parameters to estimate with the sequential approach or using design matrices for both elasticity and viscoelasticity behavior considering a dependence structure or not.
A first remark that can be made from Table 1 is the interest to separate the identification of the parameters between the different phenomena, which can be considered as an application of the sequential approach. Furthermore, it shows that the sequential approach allows to limit the number of population parameters to calibrate for each phenomena. This is particularly true for the viscoelasticity phase, in which the implementation of the sequential approach with a dependence structure cuts almost by half the number of population parameters, even if this remains difficult to carry out. Table 1 also illustrates that a dependence modeling inflates the number of population parameters. This is expected because without dependence structure, the number of population parameters is in 𝒪(d) while it is in 𝒪(d^{2}) if all dependencies are considered, with d the number of model parameters to estimate. Thus, because of the small number of available individuals, striking a balance between the complexity of the modeling and the number of available individuals appears appears as necessary, provided that the correlation coefficients demand a large amount of data to be accurately estimated. So far, only the number of distribution parameters to estimate has been discussed. Yet, with mixedeffects, the residual model (here the ω_{i}s) is also to calibrate. Here again the implementation of the sequential approach allows to reduce the number of residual models (and thus of parameters) to estimate. Indeed, with the sequential approach, a focus is made on certain types of specimen (for instance compression experiments) while when with design matrices, different types of experiments are considered (tension and compression experiments). Thus, under the reasonable hypothesis of at least one residual model for each type of experiments, there are more residual models (so more parameters) to calibrate when design matrices are employed than in the different stages of the implementation of the sequential approach. Consequently, it appears that the sequential strategy is particularly suited to reduce the complexity of the minimization problem. It should be noticed that, in this application of the sequential strategy, the calibration subproblems have been defined from the analysis of the available experiments, allowing to select the parameters of model that would be used to describe the phenomena of interest.. Furthermore, it should be emphasized that standard mixedeffects and the sequential strategy enable to characterize the same distributions. The only difference between the two methods is to know whether the full distribution is estimated straightforwardly or sequentially. The interest of this test case is that it allows to test the sequential strategy in a case where it can be compared to its direct counterpart. The next section describes the tools used to analyze the results.
For the nonlinear elastic and viscoelastic behavorial law of the OPFM model, number of population parameters to be calibrated with design matrices or using the sequential approach considering a dependence structure or not. The residual models are not taken into account. Note that for the sequential approach, the number of population parameters refers to highest number of parameters to be estimated in the different stages of the method.
4.3 Application of the sequential strategy in the mixedeffects framework
To test the ability of the proposed methodology to calibrate the elastic model with mixedeffects, it is chosen to calibrate the model with synthetic data. Both the distribution parameters and the individual parameters are are known to generate the synthetic data and are estimated with the proposed calibration approach.
A multivariate Gaussian distribution, defined by a mean vector µ_{exact} and a covariance matrix Σ_{exact} is chosen to generate the synthetic data. This set of parameters will be denoted Π_{exact} in the following. Individual parameters θ_{i,exact} are sampled from this distribution to compute, for each of these individuals the elastic model outputs ε_{i,exact}. The mixedeffects approach is applied to the synthetic data and the results obtained are confronted to their exact counterpart.
The results are first compared at the population level to ensure that the model parameters distribution is properly estimated, using the relative error on the distribution parameters , defined as:
noting 𝓁 the number of distribution parameters and the KLdivergence [60], , ), between the exact and the estimated distribution :
In general, it is computed with coarse MonteCarlo. Yet, if both the exact and calibrated model parameters distribution are Gaussian multivariate distribution of dimension d parameterized by Π_{exact} = [μ_{exact}, Σ_{exact}] and respectively, the KullbackLeibler (KL) divergence can be estimated in closed form and refers to [61]
The KLdivergence expresses a notion of distance between two distributions. One interest of such measure is that it accounts simultaneously for all aspects (marginals and correlations) when the error on the population parameters only studies the calibration of one population parameter at the time. In addition, it is here specifically relevant as the population parameters are estimated jointly and the KLdivergence can be understood as a joint measure of the errors. The adequation between the calibrated and the synthetic data can also be evaluated for each individual. To verify that the different individuals are correctly identified either in terms of calibrated individual parameters or estimated output ℱ(t_{i}, ), the averaged relative error on the individual parameters e(,… ) is defined as:
and the average error d(,…, ) between the exact and estimated strains for each coordinate
are computed. Furthermore, to account for the specificity of the sampled individuals, the calibration is repeated 20 times with different individuals generated from the same joint distribution . The above indicators are averaged over the repetitions and their dispersion is characterized by the coefficient of variation COV(X) defined for any nonzero quantity X as with (X) the empirical mean of X and (X) its empirical variance.
Finally, the maximum likelihood estimates (MLE) also enables to determine the variance of the residual model. To keep the number of parameters reasonable and ease the minimization of the loglikelihood, it is decided here to assign to each individual and measure the same covariance matrix of the discrepancy term:
Exact values of the parameters distribution.
4.4 Generating virtual data
The exact means and standarddeviations of the parameters for the joint distribution are described in Table 2.
Mean values are consistent with material properties of composite material with carbon fibers and epoxy resin such as those described in Laurin [56]. Standarddeviations of and v_{12} are chosen to fit with the values of dispersion given by Rollet [62]. Several correlations scenarios are considered: independence (between and ), positive correlation : ρ(, ) = 0.809, ρ(v_{12}, ) = 0.730 and negative correlations ρ(,v_{12}) = −0.500, ρ(, v_{12}) = −0.741, ρ(, ) = −0.414, noting ρ the correlation. For each of the 20 repetitions, 50 independent identically distributed vectors of parameters corresponding to 50 virtual individuals are sampled from to achieve statistical consistency of the estimators. Note that 50 repetitions in an experimental context can be considered as much, but with synthetic data, it eases the analysis of the calibration results.
For each of these samples, given a tension stress that goes from 0 to 1500 MPa with 32 points regularly spaced and a compressive load that goes from 0 to − 1000 MPa with 28 points regularly spaced, the corresponding model outputs (ε_{i,exact})_{iϵ〚1, n〛} are computed. An heteroscedastic noise is added to the experimental data:
with the unitary vector of ℝ^{2}, the zero vector of ℝ^{2} and I_{2} the identity matrix of ℳ_{2}(ℝ). Yet, in the calibration process, the noise is considered as homoscedastic for all measures (cf. Eq.(7)) as otherwise, this would make too much residual models to estimate (up to 2800 in compression and 3200 in tension). Exact and noisy data are depicted in Figure 6.
4.5 Calibration with synthetic data
This section is dedicated to the presentation of the calibration results in the mixedeffects framework with the sequential strategy using synthetic data. First, the calibration results with the tension experiments are investigated to check that the joint distribution of ,v_{12} and is properly estimated in Section 4.5.2. Then, the results of the second stage, that is to say with the compression experiments are analyzed to assess the ability of the sequential procedure to estimate the full joint probability distribution in Section 4.5.3. Both orders are considered (tension to compression and then conversely), before being compared to the straightforward calibration of the full joint distribution using design matrices in Section 4.5.4.
Fig. 6 Example of a set of 50 synthetic stressstrain curves, without and with added noise (lines and dots, respectively). The longitudinal (on the left) and transverse responses (on the right) are depicted for both the tension (upper line) and compression experiments (lower line). 
Search space for the population parameters.
4.5.1 Numerical settings
The implementation of the code has been carried out in Python with the NumPy library [63] for matrix calculation and OpenTURNS library [64] for the probabilistic modeling.
The global optimization is carried out with the CovarianceMatrixAdaptationEvolutionStrategy [44]. The optimization variables are normalized between the 0 and 1. The bounds on the population parameters are reported in Table 3. They are determined in two stages. First, the individual parameters are estimated, from which empirical means and standarddeviations are computed. They are used to deduce upper and lower bounds. Note that these bounds can be adapted if deemed necessary when optimization variables hit either upper or lower bounds. The parameters for the correlations are bounded between 0 and π. These mathematical bounds cannot be changed to ensure the positivedefinite character of the correlation matrix. Those on the variance of the residual models are determined by trialerror. The initialization point is taken as with Ƥ the dimension of Ψ. The population suggested by [44] corresponds to 4 + 3ln(Ƥ) (the upper rounded value is chosen). The maximum number of iterations is set to 400 to carry out the global optimization and appears to be sufficient in the numerical simulations. The tolerance criterion for the negative loglikelihood is 10^{−6} and the same for the criterion is chosen for the optimization variable.
The auxiliary optimization to identify the individual parameters is the NelderMead algorithm in NLopt library [65]. One of its main advantage is that it is derivative free and thus more robust than gradient descent as it does not require any gradient or Hessian approximation, and it does not require the same number of model evaluations than evolutionary that are computationally intensive. To avoid local minimas, a multistart version of the algorithm is chosen. 50 initialization points are considered, whose one correspond to the mean vector which is the most probable value for the individual parameters. The 49 others are determined from a Latin Hypercube Sampling (without random shuffle) available in OpenTURNS [64], and like for the global optimization, the variables are normalized between 0 and 1. The individual parameters are bounded between [0.3µ, 1.7µ].
Averaged calibrated mean, standarddeviations and correlations Π_{1,mean} and averaged error ε_{mean}() over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.
4.5.2 Calibration of the joint distribution of , and ν_{12}
Given the modeling assumptions described in Sections 4.2 and 4.5, there are Ƥ = 11 parameters to be calibrated in this first step corresponding to the tension/compression step: 3 mean parameters, 3 standarddeviation parameters, 3 correlation parameters and 2 parameters for the residual model. The calibration of the dependence structure is achieved implementing a spherical parameterization of the Cholesky decomposition of the covariance matrix [66]. This technique has multiple advantages, notably the guarantee to provide a positivedefinite matrix for Σ without adding constraints to the optimization algorithm. Another interest is that the optimization variables are easy to bound, contrary to those that appear in the Cholesky decomposition, which is helpful when defining the trust regions.
To start the analysis result, the estimation of the calibrated means, standarddeviations and correlations reported in Table 4 can be discussed.
These results indicate that all the mean parameters are properly calibrated as the average relative error remains below 1% in average over all the repetitions. The standarddeviation parameters and the correlations are not estimated with the same accuracy, which can be understood as they require more samples to converge to their exact values. Still, the standarddeviations of both and ν_{12} similarly to the correlation between them are estimated with acceptable level of errors (at most 15% in average). It is interesting to notice that even if 50 repetitions of experiments can be considered as a lot in a mechanical context, it remains relatively few in a statistical context as some repetitions reach high levels of error for the population parameters, as indicate the COVs.Remark that all population parameters linked to except its mean exhibit higher level of errors than the others. This suggests a poorer calibration of the corresponding individual parameters. To explain this point, remember that mainly expresses at the beginning of the strainstress curve, where the residuals have the lowest values, while expresses at the end of the experiment, for the highest residual values. Thus, function ɡ_{i}(·) defined in equation (17) is necessarily more sensitive to the latter parameter than to that can be flawed with more error. Studying the correlations, the higher (in absolute value) the exact correlation, the more accurate the estimation. This is not surprising if remembering that correlations can be understood as constraints between the model parameters. The higher (in absolute value) it is, the more it coerces the individual parameters to follow a given pattern which makes the likelihood function sensitive to it, allowing to provide an accurate estimation. Globally, the calibration gives estimation of the joint distribution of , and ν_{12} that can be considered as rather accurate provided the difficulties coming from the concentration of the linear compliance around its mean.
To check the consistency of the calibration, it is mandatory to study the calibration of the individual parameters as they influence the determination of the population parameters. Table 5 reports the error on the individual parameters as well as the error on strain outputs. It also contains the measure of the spread of each model parameter labeled ς_{exact} (instead of the empirical COV), that is to say with µ_{exact} and Σ_{exact} the exact mean vector and covariance matrix. It measures how much each parameter is concentrated around its mean values. In fact, it is a coefficient of variation that is here renamed, to prevent confusions with the empirical coefficient of the averaged relative errors (denoted COV), that also measures the concentration around their mean value for several quantities that are different from θ.
In Tables 5 and 6, the individual parameters are, in average, always properly estimated with a maximum error of at most 1.6% for each parameter. This a key feature as a poor calibration of the individual parameters would prevent a proper estimation of the joint distribution of the model parameters. Furthermore, in model space, the square root of the model error reaches about 10^{−5} for both strain component, that can be considered low compared to the strain values that are above 10^{−3} for most of the experiments. This confirms the good calibration of the experiments in average for all repetitions. It has been stated earlier that the linear compliance was not estimated as accurately as it could be because the residual of the nonlinear part outweight those of the linear behavior monitored by . Yet, this statement seems to be in contradiction with the level of error. Meanwhile, a comparison to the exact spread values ς_{exact} can bring some elements of justification. Indeed, this variable is the most concentrated around its mean value, especially compared to . Thus, within [ ,], (the 99% credible interval for Gaussian distributed variables), most of the parameters will give the same model error, which makes them difficult to identify exactly. This effect has a small impact on the mean parameter because it is a position parameter and the same interval of linear compliances is always targeted. However, this gets different for the standarddeviation which indicates the dispersion of the model parameters and depends on the position of the individual parameters with respect to the mean value, making it more sensitive to estimation errors. Still, despite of the limitations highlighted here, these different indicators allow to validate the calibration of joint distribution of , and v_{12} either at the population or individual levels. Consequently, the first stage of the methodology is completed.
Relative errors on the estimation of the individual parameters e(,…, ) averaged over 20 repetitions of the calibration process with 50 individual along with the spread ς_{exact}. The corresponding coefficients of variation in % are indicated between brackets.
Error in model space d(,…, ) for each strain component.
4.5.3 Calibration of the joint distribution of , and v_{12}
In the second stage of the methodology, the objective is to identify the joint distribution of , and v_{12}. Furthermore, the parameters of the joint distribution of v_{12} and , , 𝕍(v_{12}), 𝕍() and ρ(v_{12}, ), already estimated in the previous step, will be limited to a trust region corresponding to ℛ_{1} = [0.8Π_{1}, 1.2Π_{1}]. Considering 20% of variation both enables a small exploration around the calibrated value of Π_{1} and to remain concentrated around them. To bound the correlation parameter, ρ(v_{12}, ) , a spherical decomposition of the covariance matrix is used [66]. The number of population parameters to estimate remains the same ( 11 ), with six of them limited to 7 ℛ_{1} .
First, one can notice that the distribution parameters involving (not calibrated in step 1) are well estimated. Indeed, the average relative error on the marginal parameters is, over all the repetitions, below 10%. The correlation coefficients show larger level of errors, but it can be understood as it is rather weak (−0.41) which adds to the limited number of individuals. Otherwise, the calibration of the population parameters that involve can be considered as satisfactory, provided that this parameter is the most dispersed around its mean value (its exact coefficient of variation is about 10% which makes it more subjected to estimation errors). It is now possible to focus on the joint distribution of and v_{12}. From Table 7, it seems that the estimation on this joint distribution is downgraded between the two steps. Indeed, except for the mean parameters, the two standarddeviations appear to be flawed from their initial estimation (Tab. 4), in particular for the standarddeviation of (the averaged error increases by 30%). Indeed, both in terms of averaged values or averaged error on the population parameters. The correlation between and v_{12} does not seem to be much affected and even exhibits a small improvement. Rather than focusing independently on each parameter, a more suited way to quantify this alleged degradation is the KullbackLeibler (KL) divergence between the exact joint distribution of , v_{12} (ƒ_{exact}(, v_{12})) and the calibrated distribution (ƒ_{calibrated}(, v_{12} )) in the two stages as depicts Figure 7 for the 20 repetitions. The main interest of such a quantity is that it takes into account all aspects within the same quantity, which sounds more appropriate as all parameters were calibrated jointly.
First, one can notice that the distribution parameters involving (not calibrated in step 1) are well estimated. Indeed, the average relative error on the marginal parameters is, over all the repetitions, below 10%. The correlation coefficients show larger level of errors, but it can be understood as it is rather weak (−0.41) which adds to the limited number of individuals. Otherwise, the calibration of the population parameters that involve can be considered as satisfactory, provided that this parameter is the most dispersed around its mean value (its exact coefficient of variation is about 10% which makes it more subjected to estimation errors). It is now possible to focus on the joint distribution of and ν_{12}. From Table 7, it seems that the estimation on this joint distribution is downgraded between the two steps. Indeed, except for the mean parameters, the two standarddeviations appear to be flawed from their initial estimation (Tab. 4), in particular for the standarddeviation of (the averaged error increases by 30%). This downgrade can be noticed both the averaged values or averaged errors on the population parameters. The correlation between and ν_{12} does not seem to be much affected and even exhibits a small improvement. Rather than focusing independently on each parameter, a more suited way to quantify this alleged degradation is the KullbackLeibler divergence between the exact joint distribution of , ν_{12} (ƒ_{exact}(, ν_{12})) and the calibrated distribution (ƒ_{calibrated}(, ν_{12})) in the two stages as depict Figure 7 for the 20 repetitions.
In Figure 7, the deterioration of the joint distribution of and ν_{12} suggested by Table 7 seems to be confirmed. For instance, the median of the KLdivergences increases (multiplied by 3), and the ranges of variation of the KLdivergences are also increased by two (without consideration of the outliers). This seems mainly due to poorer estimation in the compression phase of the standarddeviation of . Consequently, a deterioration of the estimation of joint distribution of and ν_{12} can be observed in the second phase. Nevertheless, this deterioration remains limited for both the relative errors and the KLdivergences. In particular, even if they are larger in the second phase than in the first phase, the KLdivergences have the same orders of magnitude. Indeed, except for one repetition, the KLdivergences range from 0 to 0.7 in the first step and from 0 to 1 in the second step, which can be considered as a clear but small downgrade. This does not necessarily invalidate the sequential method, that would require a comparison with the direct estimation method proposed in Section 4.5.4 but rather hints that the calibration is harder to perform using the compression experiments instead of tension experiments. One possible explanation could be a poorer calibration of the individual linear compliances (to be investigated later). To finish the analysis, it remains to check the calibration of the individuals, either in the parameters and model space. The suited indicators are displayed in Tables 8 and 9 that also contain, as earlier, the measure of the spread of each model parameter labeled Çexact (instead of the empirical COV), that is to say with µ_{exact} and Σ_{exact} the exact mean vector and covariance matrix.
From Tables 8 and 9, the same remarks made on Tables 5 and 6 arise. Indeed, in average, for all repetitions, all experiments seem properly calibrated either in model or parameters space. However, with respect to their exact spread, the nonlinear compliance can be considered as better calibrated than the others. The only difference is that the elastic linear parameters , ν_{12} are flawed with higher errors, partly explaining the poorer calibration of their joint distribution highlighted above. This could come from the nonlinearity of the strainstress curve more pronounced in compression than in tension, making the sum of the residuals of the nonlinear part greater than those from the linear part, here reduced to few points, but this remains to be confirmed by a detailed study.
The results from this section and the previous one show that by conducting the calibration starting with the tension experiments and then continuing with the compression experiments, it is possible to estimate accurately the full joint distribution of the longitudinal nonlinear elasticity provided the limitations of the model and of population parameters values. Still, the sequential calibration has been performed in an arbitrarily chosen order and there is no particular reason to favor this one to others. Thus, it is necessary to check to what the extent the choice of the order and more generally the calibration strategy influences the calibration results. In the next section, the sequential strategy is achieved in the reverse order to check in particular whether this choice downgrades the estimation of the model parameters distribution, before being compared to the method that implements design matrices.
Averaged calibrated mean, standarddeviations and correlations Π_{2,mean} and averaged error ε_{mean} () over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.
Fig. 7 KLdivergences of the calibrated and exact joint distribution and ν_{12} for the 20 repetitions in the two stages. Remember that the KLdivergence expresses a notion of distance between two distributions and the lower is its value, the more accurate the calibration of ƒ_{ΘΠ} is. 
Relative errors on the estimation of the individual parameters e(,…, ) averaged over 20 repetitions of the calibration process with the spread ς_{exact}. The corresponding coefficients of variation in % are indicated between brackets.
Error in model space d(,…, ) for each strain component.
4.5.4 Comparison of the different calibration strategies
In this section, the objective is to present the results provided by the other solutions that allow to estimate the model parameters distribution. Remember that among these solutions can be found either the implementation of the sequential strategy in the reverse order: compression first, then traction, and with the use of design matrices. Results from the sequential calibration in the reverse order are presented, before a comparison with the calibration with design matrices. Note that all the optimization settings for the sequential scheme are kept the same as in Section 4.5.2. The main issue here is to know whether one of the orders provides better results (in terms of errors). To investigate this point, the result of main interest is directly studied, i.e., how does evolve the calibration of the joint distribution of , v_{12} ? In particular, is the estimation downgraded or improved between the two stages ? To answer this question, the KLdivergence between the exact joint distribution of , ν_{12} (ƒ_{exact}(, ν_{12})) and the calibrated distribution (ƒ_{calibrated}(, v_{12})) in the two stages can be estimated.
Here, contrary to the previous case, there is a clear improvement of the calibration results between the two phases. Indeed, the median of the KLdivergences is almost divided by 7. This is interesting because, here, the tension phase corrected the errors made in the compression phase that appears harder to conduct as shown by all the indicators studied so far. In fact, the observed decrease mainly comes from the fact that, with the compression experiments, the estimation of the population parameters that concern the joint distribution of and v_{12} shows higher levels of error compared the to estimate of the population parameters from the tension experiments, as it arises from Tables 4 and 7. Thus, choosing this order of calibration implies an improvement when the other was synonymous of a deterioration. It is also worth noticing that, with the compression experiments, the joint distribution of the and v_{12} is estimated with increased accuracy when its calibration is carried out in the second phase rather than in the first phase, as show the KLdivergences in Figures 7 and 8. Indeed, though their medians are comparable in both situations, the range of variation of the KLdivergences (without consideration of the outliers) in the second phase of the calibration process is about half of their range of variation in the first phase of the calibration process when starting with the compression experiments. This illustrates the interest to use a sequential scheme. Indeed, the MLE calibrated in the first phase (conducted with the tension experiments and described in Sect. 4.5.2) is used to define a trust region, which reduces the size of the search space by focusing the optimization within a region that concentrates proposals of population parameters that seems to be close to the value of MLE of the second phase. This eases the optimization and seems to limit the deterioration of the population parameters estimates, contrary to the situation in which the population parameters are learnt from scratch (that is to say when the compression experiments are used first). Still, this assumption would demand to be consolidated with uncertainty quantification techniques and confirmed by appropriate convergence studies. It can be interesting to know whether this is a global improvement (i.e., the estimation of all parameters was improved) or if it is concentrated on one specific parameter. To study this particular point, the averaged error on the population parameters of the joint distribution of and ν_{12} whose values are reported in Tables 10 and 11 can be investigated.
From these tables, it follows that the estimation of all the population parameters is slightly improved from the first to the second phase except , whose averaged estimation error is reduced by about 30%. This mainly comes from the fact that the linear compliances are identified with better accuracy with the tension experiments than with the compression experiments, as show Tables 5 and 8. Still, it can be noticed that the calibration of the standarddeviation of v_{12} and of the correlation is not improved as for the other parameters given that the coefficient of variation of their corresponding averaged error is increased, indicating that this improvement is not uniform among the repetitions. However, globally, with this particular model, provided the modeling assumptions and the values chosen for population parameters, it seems better to identify first the joint distribution of the population parameters in compression, before a refinement with the tension experiments. Of course, this conclusion does not apply in all cases and is only specific to our model of interest, parameters value and modeling choices.
So far, the results provided by a sequential strategy have only been discussed. However, in this particular case, it is possible to use classical mixedeffects to estimate, with the help of design matrices, the 13 parameters of the joint distribution of , , and ν_{12} (recall that Cov(, ) cannot be estimated with the actual data). Furthermore, with the two different kind of experiments considered here, the tension and compression experiments, a reasonable assumption is to consider different residual models for each kind of experiment. Provided that there are two different curves to identify (longitudinal and transverse measurements), that makes 4 parameters for the two residual models. In total, p = 17 parameters have to be estimated. This illustrates a first interest to use a sequential scheme as with such techniques, only p = 11 parameters have to be calibrated in the different phases. The three calibration techniques are studied here. They respectively stand for the calibrated distribution in the second phase of the first implementation of the sequential strategy (the one obtained after the calibration with the compression experiments) labeled C_{1}, the calibrated distribution in the second phase of the second implementation (the one obtained after the calibration with the tension experiments) C2 and finally, the one that comes from the complete calibration the joint distribution of , , and ν_{12}, labeled C_{3}. Several items can be investigated among which the error on all population parameters in the three cases, as reported in Table 12.
In Table 12, it is noticeable that the estimation of the marginals of the asymptotic compliance in tension and compression are slightly downgraded in case C_{3} compared to cases C_{1} and C_{2} . This could come from the fact in case C_{3}, all variables are considered together and may be more difficult to tune than when the population parameters corresponding to a single type of experiment are estimated. Still, the differences remain small and may also originate from numerical issues. Furthermore, it appears that, surprisingly, the calibration that considers both compression and tension does not necessarily provide the best results for the calibration of the joint distribution of and ν_{12}. Indeed, as these parameters express on both types of experiments with 50 individuals each (so 100 individuals in total), we could expect it to be the best solution among the three possibilities given that this distribution is estimated with the highest number of individual possible. Obviously, it is not true in this particular case. In fact, it seems that the calibration process makes tradeoffs between the tension and compression experiments, and provided that the linear compliances are identified with less accuracy on the compression than the tension experiments. The phenomenon of tradeoff can be illustrated in the parameters space by depicting the PDFs of the elastic linear parameters, estimated with either the tension or compression experiments in Figure 9.
In Figure 9, on the left, it can be noticed that, if the calibration is carried out only from the compression experiments, because of the higher levels of error, the marginal parameters of are estimated less accurately than when they are estimated from only the tension experiment or when both types of experiments are jointly considered, which is consistent with previous results. For this repetition, adding the compression experiments to the tension one slightly decreases the error on the marginal estimation. This seems to be linked to the level of errors on the individual parameters as discussed earlier. The calibration of the other marginals and correlations do not exhibit specific problems provided the limitations of the population parameters values and the model issues. The KLdivergence is computed for two distributions: the joint distribution of and v_{12} on the one hand and on the other hand for the full joint distribution, and illustrated in Figure 10.
In Figure 10, it is interesting to notice that the KLdivergences are similar in terms of values and range of variation, either for the joint distribution of and v_{12} or for the complete joint distribution. Still, it is possible to observe the KLdivergences have significant range of variation (from 10^{−1} to over 10^{8}). This can be partly explained by the fact that the complete distribution is parameterized by 13 parameters (against 5 for the joint distribution of and v_{12}) whose 2 correlations exhibit large level of errors (in particular ρ(,)). These results show that the three different cases provide similar results in terms of accuracy of the calibrated distribution, suggesting that implementing a sequential approach does not downgrade much the estimation of the population parameters. A final comparison can be made to study the computational behavior of the sequential scheme with respect to the method that uses design matrices. Two indicators are studied here: the number of (material) model evaluations and the number of likelihood evaluations, all reported in Figure 11 for the three cases over the twenty repetitions.
Figure 11 shows that with a sequential scheme, in total over the two phases, the likelihood function is evaluated a comparable number of times in each case, even if in Case C_{2} (when ending with the tension experiments), a slight increase can be noticed. Still, the increase remains an within acceptable range. On the contrary, the interest of using a sequential approach arises with the study of number of model evaluations, which is of prime interest. Indeed, the costly part of the likelihood computation refers to the resolution of the minimization problem defined in equation (20). Focusing on the number of model evaluations demonstrate the interest of using a sequential scheme as over the two phases, the number of model evaluations is approximately divided by two from case C_{3} to case C_{1} and by 40% from case C_{3} to case C_{2}, which can be considered as a significant gain in terms of computational cost. Still, the number of model evaluations remains high, about 10^{8}. This number can be partly explained by the choices and settings of the optimization algorithms used for the minimization of both the opposite loglikelihood and function g_{i}(·). For more complex, it may be necessary to combine the sequential strategy with other approaches such as surrogate modeling to guarantee the appliability of the calibration method. If such solution is chosen, the surrogate model would be easier to train compared to the standard case. Indeed, throughout the sequential process, the search spaces over which population parameters evolve are limited to preserve consistency between the different steps. Therefore, the range of validity of the chosen surrogates can be smaller, making them easier to build, contrary to the standard case in which the entire search space is considered. Note that, as in most procedures implementing surrogate models, it should be checked that the trained surrogate model does not introduce errors that could impair the calibration process. The observed reduction can be explained by different reasons. First, notice in Figure 11, for the three cases, the number of likelihood evaluations are (in total for cases C_{1} and C_{2} ) similar. In fact, the number of likelihood evaluations in each stage of the sequential approaches is about half of the number of likelihood evaluations when design matrices are used. This can be understood as there more parameters to estimate (17 instead of 11). However, for case C_{3}, for each likelihood evaluation, the models for both the compression and tension experiments are identified. This adds to the fact that, in the second step of the methodology, the individual calibrations are eased.Indeed, the mean parameters around which the individual deviations b_{i} are searched (see Eq. (20)) are already estimated with few error (smaller than 1% in average). This is even more important here because of the concentration of the linear elastic parameters around their average values. This is illustrated in Figure 12, which compares for the twenty repetitions, the number of model evaluations for the second steps of the sequential strategy, i.e., case C_{1} with the compression experiments and case C_{2} with the tension experiments to their counterpart of case C_{3}.
Figure 12 confirms that thanks to the focus realized by the first identification for the mean parameters, it seems easier to identify the individual parameters. Indeed, the number of model evaluations is almost cut by half between the case C3 and the second stages of both cases C_{1} and C_{2} . This significant improvement partly explains the gain observed in terms of model evaluations. A further study on the convergence behavior should be carried out to demonstrate rigorously the interest of the sequential approach compared to classical mixedeffects methods with design matrices.
This study of the sequential scheme allows to illustrate some advantages and drawbacks. Among the benefits of such a strategy can be cited a proper estimation of the population parameters, at least as accurate as the one that can be obtained with the use of design matrices in this testcase. Furthermore, it allows to significantly reduce the number of model evaluations necessary to carry out the calibration of the population parameters. Finally, it also reduces the number of population parameters of estimated in each calibration problem, allowing to consider more complex phenomena in future works. However, the process exhibits its own limitations. For instance, it could be interesting to test different solutions for the trust regions instead of a simple hypercube. In addition, the routine could have been continued until reaching convergence as in a fixedpoint approach at the expense of a larger computational costs. Still, this testcase is interesting as the results can be compared to the standard calibration process that sometimes cannot be carried out because of the computational costs and the dimension of the optimization problem.
Fig. 8 KLdivergences of the calibrated and exact joint distribution and ν_{12} for the 20 repetitions in the two stages. Remember that the KLdivergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒ_{Θ∣Π} is. 
Averaged marginal parameters (∏_{1,mean} from the first phase and ∏_{2,mean} from the second phase), alongside the averaged error (ε_{mean}() and ε_{mean}()) over 20 repetitions of the calibration process with samples of 50 individuals each. The first two columns refer to the error made on the population parameters estimated in the first stage ∏_{1} (with the compression experiment) and the two left to parameters ∏_{2} estimated in the second phase (with the tension experiment). The coefficients of variation COV in % are indicated between brackets.
Averaged marginal parameters correlation parameter in both phases and averaged error ε_{mean}()_{mean,1/2} over 20 repetitions of the calibration process with samples of 50 individuals each. The first column refers to the error made on the population parameters in the second stages of the sequential methodology (with the compression experiment on the left and with the tension experiment on the right). The coefficients of variation COV in % are indicated between brackets.
Averaged error ε_{mean}(Π) on the population parameters calibrated in the three cases over 20 repetitions of the calibration process with samples of 50 individuals each. For cases C_{1} and C_{2}, the error on the parameters of the joint distribution of and v_{12} correspond to those obtained after completion of the sequential scheme, that is to say after the second step. The coefficients of variation COV in % are indicated between brackets.
Fig. 9 For the 17^{th} repetition, PDF of and ν_{12} and the corresponding individual parameters for the calibration from case C_{1} (with the compression experiments in the second phase), from case C_{2} (with the tension experiments in the second phase) and considering both types of experiments (case C_{3}). 
Fig. 10 On the left, KLdivergences between the calibrated and exact joint distribution of and v_{12} and on the right, KLdivergences between the calibrated and exact joint full distributions for the 20 repetitions and 3 cases. Remember that the KLdivergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒ_{ΘΠ} is. 
Fig. 11 On the left, number of model evaluations for the three cases and on the right, number of likelihood evaluations for the three cases. In cases C_{1} and C_{2}, the costs of both stages are combined together. 
Fig. 12 On the left, total number of model evaluations for the estimation of the individual parameters with the compression experiments in case C_{3} and on the second stage of case C_{1}. On the right, total number of model evaluations for the estimation of the individual parameters with the tension experiments in case C_{3} and on the second stage of case C_{2}. 
5 Conclusion
This article presents a sequential calibration procedure that aims to calibrate complex material constitutive models compliant with the mixedeffects framework using multivariate data. The method consists in defining a sequence of calibration subproblems, each one corresponding to a subset of the available experimental data. The subdivision of the complete problem into subproblems can be based on the choice of the expert and on sensitivity analyses. Indeed, wellchosen tests can activate different parts of the model, i.e. different subsets of the model parameters. Thus, each type of test makes it possible to define a subproblem whose objective is to identify the joint distribution of a subset of the model parameters. A parameter can be involved in several subproblems. In this case, it is identified once and the search space of the corresponding distribution parameters is reduced to a trust region in the subsequent calibrations, to ensure consistency between the different subproblems.
The method is applied to a virtual test case to the calibration of the orthotropic model in plane stress with nonlinear longitudinal behavior. This testcase is interesting as it allows to compare the sequential strategy to its usual counterpart using design matrices that is not necessary available in all configurations. A first implementation allows to split the identification of the joint distribution of the compliances with the usual elementary UD 0°, UD 90° and ±45° laminates. The second allows to divide the calibration of the joint distribution of the longitudinal and transverse nonlinear elasticity of UD0^{°} laminates between the tension and compression experiments. The consistency of the results throughout the calibration process is ensured by bounding the parameters estimated twice. The results of the methodology applied on synthetic data generated with the exact model demonstrate the ability of the proposed procedure to estimate properly the model parameters distribution in the presence of significant material variability. The results obtained considering the different orders of calibration are investigated and compared to the straightforward calibration of the complete joint distribution using design matrices. It shows that in particular this calibration method is not necessarily the best in this test case, the one step calibration method is not necessarily the best, in particular because of the tradeoffs made between the different type of experiments. An application (not shown in this paper) with experiments carried on a material with carbon and epoxy resin confirm this point and shows in particular that with a sequential strategy, the number of model evaluations to provide an estimation of the full joint distribution is significantly reduced. This is a key feature, and combined to the preservation of the calibrated distributions throughout the calibration process, it shows that sequential calibration is not only compatible with mixedeffects models but scales better to more complex model with a high number of parameters to be calibrated.
Though the proposed method seems promising, the proposed application case is simple and very fast to evaluate. More complex models usually require to solve differential equations which can be much more time consuming, making the overall computational cost of the method intractable with direct calls to the model. In that regard, combining the proposed approach with surrogate models falls as a direct perspective of this work.
References
 U.D.O. Defense, Technomic Publishing Company, Materials Sciences Corporation, American Society for Testing, and Materials, Composite Materials HandbookMIL 17, Vol. III: Materials Usage, Design, and Analysis, The Composite Materials HandbookMIL 17. Taylor & Francis, 1999 [Google Scholar]
 A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, 2005. [CrossRef] [Google Scholar]
 L. Wasserman, All of Statistics: A Concise Course in Statistical Inference, Springer Publishing Company, Incorporated, 2010 [Google Scholar]
 A. Antoniadis, R. Carmona, Regression non lineaire et applications. Economie et statistiques avan¸c´ees, Economica, 1992 [Google Scholar]
 G.C. Ballesteros, P. Angelikopoulos, C. Papadimitriou, P. Koumoutsakos, Bayesian hierarchical models for uncertainty quantification in structural dynamics, Vulnerabil. Uncertain. Risk: Quantif. Mitig. Manag. 162, 1615–1624 (2014) [CrossRef] [Google Scholar]
 M. Song, I. Behmanesh, B. Moaveni, C. Papadimitriou, Accounting for modeling errors and inherent structural variability through a hierarchical bayesian model updating approach: an overview, Sensors 20, 3874–3901 (2020) [CrossRef] [PubMed] [Google Scholar]
 E. Demidenko, Mixed Models. Theory and Applications with R, John Wiley & Sons, 2013 [Google Scholar]
 M. Lavielle, Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools, CRC press, 2014 [CrossRef] [Google Scholar]
 P. Congdon, Bayesian Statistical Modelling, John Wiley & Sons, 2007 [Google Scholar]
 C. Laboulfie, M. Balesdent, L. Brevault, F.X. Irisarri, J.F. Maire, S. Da Veiga, R. Le Riche, Calibration of material model parameters using mixedeffects booktitle = 25ème Congres Français de Mecanique, model, 2021, pp. 258–295 [Google Scholar]
 J. Lemaitre, J.L. Chaboche, Mechanics of solid materials, Cambridge University Press, 1994 [Google Scholar]
 J. Lemaitre, R. Desmorat. Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures. Springer Science & Business Media, 2006 [Google Scholar]
 S. Moreau, A. Chrysochoos, J.M. Muracciole, B. Wattrisse, Analysis of thermoelastic effects accompanying the deformation of pmma and pc polymers, Comptes Rendus Mec. 333, 648–653 (2005) [CrossRef] [Google Scholar]
 R. Schapery, On the characterization of nonlinear viscoelastic materials, Polymer Eng. Sci. 9, 295–310 (1969) [CrossRef] [Google Scholar]
 J.L. Chaboche, A review of some plasticity and viscoplasticity constitutive Theories, Int. J. Plasticity 24, 1642–1693 (2008) [CrossRef] [Google Scholar]
 B. Miled, I. Doghri, L. Delannay, Coupled viscoelasticviscoplastic modeling of homogeneous and isotropic polymers: numerical algorithm and analytical solutions. Comput. Methods Appl. Mech. Eng. 200, 3381–3394 (2011). [Google Scholar]
 A. Krairi, I. Doghri, A thermodynamicallybased constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage, Int. J. Plasticity 60, 163–181 (2014) [Google Scholar]
 F. Praud, G. Chatzigeorgiou, J. Bikard, F. Meraghni, Phenomenological multimechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation, Mech. Mater. 114, 9–29 (2017) [Google Scholar]
 C. Hochard, P.A. Aubourg, J.P. Charles, Modelling of the mechanical behaviour of wovenfabric cfrp laminates up to failure, Composites Sci. Technol. 61, 221–230 (2001) [CrossRef] [Google Scholar]
 P. Ladeveze, E. LeDantec, Damage modelling of the elementary ply for laminated composites, Composites Sci. Technol. 43, 257–267 (1992) [CrossRef] [Google Scholar]
 A. Launay, M.H. Maitournam, Y. Marco, I. Raoult, F. Szmytka, Cyclic behaviour of short glass fibre reinforced polyamide: experimental study and constitutive equations, Int. J. Plasticity 27, 1267–1293 (2011) [CrossRef] [Google Scholar]
 F. Praud, G. Chatzigeorgiou, Y. Chemisky, F. Meraghni, Hybrid micromechanicalphenomenological modelling of anisotropic damage and anelasticity induced by microcracks in unidirectional composites, Composite Struct. 182, 223–236 (2017) [CrossRef] [Google Scholar]
 F. Lachaud, Delaminage de materiaux composites a fibres de carbone et a matrices organiques: etude numerique et experimentale, suivi par emission acoustique. PhD thesis, 1997 [Google Scholar]
 A. Forrester and A. Sobester, A.and Keane. Engineering Design via Surrogate Modelling – A Practical Guide, Wiley, 2008 [CrossRef] [Google Scholar]
 J.M. Berthelot, Materiaux CompositesComportement Mecanique et Analyse des Structures, 2012. [Google Scholar]
 M. Gallagher, J. Doherty, Parameter estimation and uncertainty analysis for a watershed model, Environ. Model. Softw. 22, 1000–1020 (2007) [CrossRef] [Google Scholar]
 M. Kennedy, A. O’Hagan, Bayesian calibration of computer models, J. Roy. Stat. Soc. B (Stat. Methodol.) 63, 425–464 (2001) [CrossRef] [Google Scholar]
 K. Levenberg, A method for the solution of certain nonlinear problems in least Squares, Q. Appl. Math. 2, 164–168 (1944) [CrossRef] [Google Scholar]
 D.W. Marquardt, An algorithm for leastsquares estimation of nonlinear parameter, J. Soc. Ind. Appl. Math. 11, 431–441 (1963) [CrossRef] [Google Scholar]
 R. Pintelon, J. Schoukens, System Identification, John Wiley & Sons, Ltd, 2012 [CrossRef] [Google Scholar]
 F. Meraghni, Y. Chemisky, B. Piotrowski, R. Echchorfi, N. Bourgeois, E. Patoor, Parameter identification of a thermodynamic model for superelastic shape memory alloys using analytical calculation of the sensitivity matrix, Eur. J. Mech. A/Solids 45, 226–237 (2014) [CrossRef] [Google Scholar]
 J. Pinheiro, D. Bates, MixedEffect Models in S and Splus, Vol. 96, 2002 [Google Scholar]
 J.B. Nagel, B. Sudret, A unified framework for multilevel uncertainty quantification in Bayesian inverse problems, Probab. Eng. Mech. 43 68–84 (2016) [CrossRef] [Google Scholar]
 S.L. Beal, L.B. Sheiner, Estimating population kinetics, Crit. Rev. Biomed. Eng. 8, 195–222 (1982) [Google Scholar]
 D. Blei, A. Kucukelbir, J. McAuliff, Variational inference: a review for statisticians, J. Am. Stat. Assoc. 112, 859–877 (2017) [CrossRef] [Google Scholar]
 J. Wakefield, The Bayesian analysis of population pharmacokinetic models, J. Am. Stat. Assoc. 91, 62–75 (1996) [CrossRef] [Google Scholar]
 P. Damlen, J. Wakefield, S. Walker, Gibbs sampling for Bayesian nonconjugate and hierarchical models by using auxiliary variables, J. Roy. Stat. Soc. B (Stat. Methodol.) 61, 331–344 (1999) [CrossRef] [Google Scholar]
 M. Davidian, D.M. Giltinan, Nonlinear Models for Repeated Measurement Data, Chapman & Hall/CRC Monographs on Statistics & Applied Probability, Taylor & Francis, 1995 [Google Scholar]
 R. Fisher, K. Pearson, On an Absolute Criterion for Fitting Frequency Curves. Messenger of Mathematics, 1911 [Google Scholar]
 D. Hall, M. Clutter, Multivariate multilevel nonlinear mixedeffects models for timber yield predictions, Biometrics 60, 16–24 (2004) [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
 J. McFarland, S. Mahadevan, L. Romeroand, V. Swiler, Calibration and uncertainty analysis for computer simulations with multivariate output, AIAA J. 46, 1253–1265 (2008) [CrossRef] [Google Scholar]
 N. Metropolis, S. Ulam, The monte carlo method, J. Am. Stat. Assoc. 44, 335–341 (1949) [CrossRef] [PubMed] [Google Scholar]
 D. Bates, D. Hamilton, D. Watts, Calculation of intrinsic and parametereffects curvatures for nonlinear regression models. Commun. Stat. Simul. Comput. 12, 469–477 (1983) [CrossRef] [Google Scholar]
 N. Hansen, S. Muller, P. Koumoutsakos, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMAES), Evol. Comput. 11, 1–18 (2003) [CrossRef] [PubMed] [Google Scholar]
 W. Oberkampf, C. Roy, Verification and Validation in Scientific Computing, Cambridge University Press, 2010 [CrossRef] [Google Scholar]
 L. Schwer, Guide for Verification and Validation in Computational Solid Mechanics, ASME, 2006 [Google Scholar]
 B. Jung, H. Yoon, H. Oh, G. Lee, M. Yoo, B. Youn, Y. Hu. Hierarchical model calibration for designing piezoelectric energy harvester in the presence of variability in material properties and geometry, Struct. Multidiscipl. Optim. 53, 161–173 (2016) [CrossRef] [MathSciNet] [Google Scholar]
 B.D. Youn, B. Jung, Z. Xi, S.B. Kim, W.R. Lee, A hierarchical framework for statistical model calibration in engineering product development, Comput. Methods Appl. Mech. Eng. 200, 1421–1431 (2011) [CrossRef] [Google Scholar]
 A. Urbina, T. Paez, T. Consulting, A bayes network approach to uncertainty quantification in hierarchicallydeveloped computational models, Int. J. Uncertain. Quantif. 2, 173–193 (2012) [CrossRef] [MathSciNet] [Google Scholar]
 J. Ye, M. Mahmoudi, K. Karayagiz, L. Johnson, R. Seede, Y. Karaman, R. Arroyave, A. Elwany, Bayesian calibration of multiple coupled simulation models for metal additive manufacturing: a Bayesian network approach. ASCEASME J. Risk Uncertain. Eng. Syst. B: Mech. Eng. 8, 011111 (2022) [Google Scholar]
 T. Kimand, B. Youn, Identifiabilitybased model decomposition for hierarchical Calibration, Struct. Multidiscipl. Optim. 60, 1801–1811 (2019) [Google Scholar]
 M. Sklar, Fonctions de repartition a n dimensions et leurs marges. in Annales de l’ISUP, Vol. 8, 1959, pp. 229–231 [Google Scholar]
 G. Gilbert, Positive definite matrices and sylvester’s criterion, Am. Math. Monthly 98, 44–46 (1991) [CrossRef] [Google Scholar]
 F. Mentre, A. Mallet, D. Baccar, Optimal design in randomeffects regression models, Biometrika 84, 429–442 (1997) [CrossRef] [MathSciNet] [Google Scholar]
 C. Vong, S. Ueckert, J. Nyberg, A.C. Hooker, Handling below limit of quantification data in optimal trial design, J. Pharmacokinet. Pharmacodyn. (2014) [Google Scholar]
 F. Laurin, Approche Multiechelle des Mecanismes de Ruine Progressive des Materiaux Stratifies et Analyse de la Tenue de Structures Composites (in French) PhD thesis, 2005. [Google Scholar]
 J. Germain, Evaluation des capacites prédictives d’un modele avance pour la prevision de la tenue de plaques stratifiees perforées (in french), Theses, Universite ParisSaclay, January 2020 [Google Scholar]
 I. Daniel, O. Ishai, Engineering Mechanics of Composite Materials, Oxford university press New York, 2006 [Google Scholar]
 J. Tomblinand, Y. Ng, S. Keshavanarayana, Material qualification and equivalency for polymer matrix composite material systems, Technical report, 2001 [Google Scholar]
 S. Kullback, R. A. Leibler, On information and sufficiency, Ann. Math. Stat. 22, 79–86 (1951) [CrossRef] [Google Scholar]
 L. Pardo, Statistical Inference based on Divergence Measures, CRC Press, 2018 [CrossRef] [Google Scholar]
 Y. Rollet, Vers une Maîtrise des Incertitudes en Calcul des Structures Composites (in French), PhD thesis, 2007 [Google Scholar]
 C.R. Harris, K. Jarrod Millman, S. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N.J. Smith, R. Kern, M. Picus, S. Hoyer, M.H. van Kerkwijk, M. Brett, A. Haldane, J. del Rio Fernandez, M. Wiebe, P. Peterson, P. GérardMarchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, T.E. Oliphant, Array programming with NumPy, Nature 585, 357–362 (2020) [CrossRef] [PubMed] [Google Scholar]
 M. Baudin, A. Dutfoy, B. Iooss, A.L. Popelin, OpenTURNS: an industrial software for uncertainty quantification in simulation. in R. Ghanem, D. Higdon, H. Owhadi (Eds.), Handbook of Uncertainty Quantification, Springer International Publishing, 2016, pp. 1–38 [Google Scholar]
 S.J. Johnson, The nlopt nonlinearoptimization package [Google Scholar]
 J. Pinheiro, D. Bates, Unconstrained parametrizations for variance=covariance matrices. Stat. Comput. 6, 289–296 (1996) [CrossRef] [Google Scholar]
Cite this article as: Clément Laboulfie, Mathieu Balesdent, Loïc Brevault, FrançoisXavier Irisarri, JeanFrançois Maire, Sebastien Da Veiga, Rodolphe Le Riche, Sequential calibration of material constitutive model using mixedeffects calibration, Mechanics & Industry 24, 32 (2023)
All Tables
For the nonlinear elastic and viscoelastic behavorial law of the OPFM model, number of population parameters to be calibrated with design matrices or using the sequential approach considering a dependence structure or not. The residual models are not taken into account. Note that for the sequential approach, the number of population parameters refers to highest number of parameters to be estimated in the different stages of the method.
Averaged calibrated mean, standarddeviations and correlations Π_{1,mean} and averaged error ε_{mean}() over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.
Relative errors on the estimation of the individual parameters e(,…, ) averaged over 20 repetitions of the calibration process with 50 individual along with the spread ς_{exact}. The corresponding coefficients of variation in % are indicated between brackets.
Averaged calibrated mean, standarddeviations and correlations Π_{2,mean} and averaged error ε_{mean} () over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.
Relative errors on the estimation of the individual parameters e(,…, ) averaged over 20 repetitions of the calibration process with the spread ς_{exact}. The corresponding coefficients of variation in % are indicated between brackets.
Averaged marginal parameters (∏_{1,mean} from the first phase and ∏_{2,mean} from the second phase), alongside the averaged error (ε_{mean}() and ε_{mean}()) over 20 repetitions of the calibration process with samples of 50 individuals each. The first two columns refer to the error made on the population parameters estimated in the first stage ∏_{1} (with the compression experiment) and the two left to parameters ∏_{2} estimated in the second phase (with the tension experiment). The coefficients of variation COV in % are indicated between brackets.
Averaged marginal parameters correlation parameter in both phases and averaged error ε_{mean}()_{mean,1/2} over 20 repetitions of the calibration process with samples of 50 individuals each. The first column refers to the error made on the population parameters in the second stages of the sequential methodology (with the compression experiment on the left and with the tension experiment on the right). The coefficients of variation COV in % are indicated between brackets.
Averaged error ε_{mean}(Π) on the population parameters calibrated in the three cases over 20 repetitions of the calibration process with samples of 50 individuals each. For cases C_{1} and C_{2}, the error on the parameters of the joint distribution of and v_{12} correspond to those obtained after completion of the sequential scheme, that is to say after the second step. The coefficients of variation COV in % are indicated between brackets.
All Figures
Fig. 1 Difference between the populationbased and the classical approaches. On the left, a), usual methods with a single parameter vector for all individuals. On the right, b), populationbased approach in the mixedmodels effects framework with f_{Θ} the distribution. 

In the text 
Fig. 2 Different levels of modeling in the population approach (adapted from [8] and [33]). The dashed lines indicate random sampling and the ellipsis the distributions. 

In the text 
Fig. 3 Recap of the mixedeffects calibration procedure. Modeling stages are highlighted in red, those that concern the calibration itself are highlighted in blue. 

In the text 
Fig. 4 Different stages of the calibration process. Once the sequence of calibration problems to resolve is defined, the calibration process starts. At each iteration, the search space is defined depending on whether parameters have already been calibrated to preserve consistency between the different steps. Calibration stages are highlighted in red and those concerning the determination of the search space are highlighted in blue. At each iteration, u_{q} refers to the integer labeling the subgroup of parameters that will be reestimated in the current iteration. N_{G} refers to the number of calibration subproblems. The process stops when all calibration subproblems are solved. 

In the text 
Fig. 5 Illustration of the role of the model parameters. The solid line indicates the model output, the blue dashed lines the asymptotes and the green dashed line the tangent at the origin. Exponent T stands for tension load and C for compressive load. 

In the text 
Fig. 6 Example of a set of 50 synthetic stressstrain curves, without and with added noise (lines and dots, respectively). The longitudinal (on the left) and transverse responses (on the right) are depicted for both the tension (upper line) and compression experiments (lower line). 

In the text 
Fig. 7 KLdivergences of the calibrated and exact joint distribution and ν_{12} for the 20 repetitions in the two stages. Remember that the KLdivergence expresses a notion of distance between two distributions and the lower is its value, the more accurate the calibration of ƒ_{ΘΠ} is. 

In the text 
Fig. 8 KLdivergences of the calibrated and exact joint distribution and ν_{12} for the 20 repetitions in the two stages. Remember that the KLdivergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒ_{Θ∣Π} is. 

In the text 
Fig. 9 For the 17^{th} repetition, PDF of and ν_{12} and the corresponding individual parameters for the calibration from case C_{1} (with the compression experiments in the second phase), from case C_{2} (with the tension experiments in the second phase) and considering both types of experiments (case C_{3}). 

In the text 
Fig. 10 On the left, KLdivergences between the calibrated and exact joint distribution of and v_{12} and on the right, KLdivergences between the calibrated and exact joint full distributions for the 20 repetitions and 3 cases. Remember that the KLdivergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒ_{ΘΠ} is. 

In the text 
Fig. 11 On the left, number of model evaluations for the three cases and on the right, number of likelihood evaluations for the three cases. In cases C_{1} and C_{2}, the costs of both stages are combined together. 

In the text 
Fig. 12 On the left, total number of model evaluations for the estimation of the individual parameters with the compression experiments in case C_{3} and on the second stage of case C_{1}. On the right, total number of model evaluations for the estimation of the individual parameters with the tension experiments in case C_{3} and on the second stage of case C_{2}. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.