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

© C. Laboulfie et al., Published by EDP Sciences 2023

Licence Creative CommonsThis 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 log-likelihood 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 mixed-effects models and the Hierarchical Bayesian models (HBMs) [9] belong to the population-based approaches. The main difference between them is that mixed-effects 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 mixed-effects approach to be able to carry out the analyses more easily. The mixed-effects 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 mixed-effects. 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 damage-plasticity models for a single ply [1922]). In addition to the in-plane features, complex constitutive models can also describe out-of-plane 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 mixed-effects 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 multi-scale 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 quasi-isotropic 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 ill-posed problems.

To solve such problems, this paper introduces a sequential procedure to calibrate material constitutive models compliant with the mixed-effects 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 sub-calibrations 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 population-based approaches and in particular mixed-effects 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 Population-based 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 = (yj)j∈〚1,N its outcome, with N the number of observations, considering one individual. Classically, the following decomposition applies:

y(,θ)+ξ,(1)

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:

yj:=(tj,θ)+ξj,(2)

with tj the jth 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 Levenberg-Marquardt 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, population-based approaches explicitly model the variability between individuals within a population. Figure 1 illustrates this fundamental difference with respect to the classical identifications.

Population-based 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 population-based 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, population-based approaches are particularly well suited to situations in which individual variability cannot be neglected with respect to other sources of uncertainty

The population-based approaches framework [7,8] assumes that there exists a probability distribution fΘ whose outcomes are the individual parameters:

i1,n,  θifΘi.i.d..(3)

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 yi can be written as

i1,n,yii.i.d.(,θi)+ξi.(4)

Without any other hypothesis, the outcomes of ξi (labeled (ξij)j1,Ni) are different for each individual and for each observation, with Ni the number of observations of the i th individual . The global mixed-effects models for the jth output measure of the ith individual yij reads as

yij=(tij,θi)+ξij.(5)

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:

fΘ|Π:=𝒩(μ,Σ) and Π:=[ μ,Σ ],(6)

with µ ∈ ℝd the mean vector and Σ∈ ℳd(ℝ) the covari-ance matrix The second hypothesis is that for each individual and each measure, the discrepancy term is a Gaussian white noise (no bias, no correlation):

ξij𝒩i.i.d.(0,wij2),(7)

with ωij the standard-deviation of the noise of the jth output measure of the ith individual Furthermore, the noise is supposed to be homoscedastic, i.e. ωij = ωij ∈ 〚1, Ni〛 [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 Ω=diag(ω12,,ωn2). The population-based approaches seek Ψ and provide an estimate of the individual parameters (θi)i∈[1,n] as a by-product. 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 mixed-effects models [7,8,34]. The notion of mixed-effects comes from the decomposition of the individual parameters θi through fixed-effects β and random effects bi as

θi=Aiβ+Bibi,(8)

where Ai, Bi are deterministic design matrices that belong to ℳd (ℝ) with d the number of model parameters. Design matrices other than Id, 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 fixed-effects, random effects and Π, but for a Gaussian fΘ|Π, the fixed-effects β are equal to the population parameters average µ and the random-effects comply with biN(0d,Σ),i ∈ 〚1, n〛. Furthermore, provided that the modeling proposed by the population exhibit a hierarchical structure, hierarchical-based 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, Markov-Chain-Monte-Carlo techniques (e.g., the Metropolis-Hastings 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, mixed-effects models provide more flexibility as they can be applied in either the frequen-tist 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 mixed-effects, 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.

thumbnail Fig. 1

Difference between the population-based and the classical approaches. On the left, a), usual methods with a single parameter vector for all individuals. On the right, b), population-based approach in the mixed-models effects framework with fΘ the distribution.

thumbnail 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 mixed-effects

The calibration is often achieved by maximizing the likelihood of Ψ, ℒ(Ψ) := f(y1,…, yn|Ψ) (where f is a generic letter for probability density functions or PDFs), even if other methods can be found to estimate Ψ [38]. The step-by-step 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 (yi|Ψ),

(Ψ):=i=1ni(Ψ).(9)

Because the θi is are not observed, the likelihood of the ith individual ℒi(Ψ) is the integral of the marginal likelihood (i.e., the density of the output data yi given individual parameters θi and parameters Ψ) f (yi|θi, Π, ωi) with respect to all possible θi over ℝd:

i(Ψ):=df(yi|θi,Π,wi)fΘ|Π(θi|Π)dθi.(10)

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 ΓΨ :

Ψ^:=argmaxΨ Γ Ψ(Ψ).(11)

In practice, the logarithm of the likelihood is computed to ease the numerical optimization [39].

2.3 Mixed-effects for multivariate models

The calibration of mixed-effects 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 mixed-effects equations [40].

Let’s note mi the number of output measures of the ith individual and k ∈ 〚1, mi〛 the corresponding index. In the following, all individuals are supposed to have the same number of output measures, so mi = mi ∈ 〚1, n〛. Moreover, all measures are supposed to have the same number of observation points Ni. The main difference with Section 2.1 is that the model output of each observation point is no longer a scalar yij but a vector yij·∈ ℝm. The same occurs for the discrepancy term, now labeled ξij.

The equations governing the mixed-effects are [7,8,32]:

{ (i,j)1,n×1,Ni,yij=(tij,θi)+ξij,i1,n,θifΘ|Π,(i,j)1,n×1,Ni,ξij𝒩(0m,Ωij),i.i.d. (12a) (12b) (12c)

with 0m 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:

i1,n,(j1,j2)1,Ni2,j1j2Cov(ξij1,ξij2)=0.(13)

Furthermore, for all individuals, the noise is supposed to be homoscedastic, that is to say: Ωij = Ω ij ∈ 〚1, Ni〛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 Yi, the matrix of the output data of the ith individual with Yi := (yi1,,yiNi) Mm,Ni(). The rows of Yi correspond to yik=(yijk)j1,Ni, that is to say the Ni observations of the k-th measure of the ith individual. Similarly, it is possible to define the matrix of the errors as ξi := (ξi1,ξiNi) Mm,Ni(). The rows of the ξi matrix correspond to ξik=(ξijk)j1,Ni and can be understood as the Ni measurement errors of the k-th measure for the ith individual. Let us define Y˜i the flattened array by row of the matrix of the output data Yi and the same for the matrix of the discrepancy term ξi. Then, Y˜i and ξ˜i belong to mNi. This allows to transform equation (12) into

{ i1,n,Y˜i=(ti,θi)+ξ˜i,i1,n,θii.i.d.fΘ|Πi1,n,ξ˜iNi.i.d.(0mNi,Ω˜i). (14a) (14b) (14c)

with where ⊗ denotes the Kronecker product, INi the identity matrix of Ni (ℝ) and ℝ the zero vector of ℝ. Hence, the set of the model parameters to be calibrated, labeled Ψ is defined as:

Ψ:=(μ,    Σ,    Ω1,   ,    Ωm)

Finally, the likelihood of the mixed-effects becomes f(Y˜1,,Y˜n|Ψ). Equation (10) remains the same except for yi that becomes Y˜i.

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 Monte-Carlo 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

A:dexp(h(x)) dx,

where h(·) is a function which must satisfy the following constraints:

  1. h(·) admits a global minimum x0 that belongs to the integration interval,

  2. h(·) is a twice-differentiable function,

  3. its Hessian matrix ℋ(h) computed at x = x0 is a symmetric definite positive matrix.

The main idea is to assume that only regions close to x0 significantly contribute to the integral, leading to the following Laplace approximation of A:

Aexp(h(x0)) )(2π)d2| (h)(x0) |.(15)

Given the modeling choices (Eq. (14)), the individual likelihoods read as:

i(Ψ)=df(Y˜i|θi,Π,Ω˜i)fΘ|Π(θi|Π)dθi               =d1| Ω˜i || Σ |(2π)d+mNi                 ×exp(gi(μ,Δ,Λi,Y˜i,bi)2)dbi.(16)

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 Ω˜i1(so Ω˜i1=ΛiTΛi) and the function ɡi(·) is defined by:

gi(μ,Δ,Λi,Y˜i,bi):= Λi(Y˜i(ti,μ+bi)) 2              + Δbi 2.(17)

The Laplace method is applied in two steps:

  1. Search for the individual parameters (or rather the deviations), b^i, minimizing ɡi(·) equation (17),

  2. 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 b^i [32]:

(gi)(b^i)=2(ti,μ+b^i)bibiTΛiTΛi(Y˜i(μ,ti,b^i))                        +(ti,μ+b^)TbiΛiTΛi(ti,μ+b^i)bi                         +ΔTΔ.(18)

In practice, 2(ti,μ+b^i)bibiTΛiTΛi(Y˜i(ti,μ+b^i)) can be neglected if the model ℱ( · , · ) is close enough to the experiment Y˜i [43]. The term (μ,ti+b^i)bi is evaluated using a finite difference scheme. As a result, the negative log-likelihood is finally expressed as

ln((Ψ))=i=1nln(i(Ψ))=i=1nln(df(Y˜i|θi,Π,Ω˜i)fΘ|Π(θi|Π)dθi)                             i=1n(12ln(| (gi)(b^i) |)+gi(μ,Δ,Λi,Y˜i,b^i)2+12ln(| Ω˜i |(2π)mNi))                              +n2ln(| Σ |).(19)

2.5 Individual parameters estimation

The mixed-effects model also allows to infer the individual parameters. With the Laplace method, individual parameters are by-products of the likelihood function calculation given Ψ. Once Ψ is chosen, the random effects bi are estimated by looking for the dominating contribution in the individual likelihood ℒi (Eq. (16)) through the minimization of function ɡi(·) (Eq. (17))

b^i=argminbiΓbgi(μ,Δi,Λi,Y˜i,bi),(20)

with Γb the domain definition for the individual deviations. The individual deviations are then simply computed as θi = µ + bi.

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 mixed-effects 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 over-parametrized). Otherwise, if the individual parameters cannot be learnt with sufficient precision, the regularization term (fΘ|Π(θi|Π) in equation (10) or ‖Δbi2 in function ɡi(·) in equation (17)) will coerce them to the mean value and thus flaw the estimation of the standard-deviations 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 mixed-effects model has been established, the model parameters can be identified by minimizing the negative log-likelihood,

Ψ^=argminΨΓΨln((Ψ)),(21)

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 gradient-based 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 bi defined in equation (20) with respect to the population parameters (i.e., µ and Σ via Δ). Moreover, the log-likelihood function often presents numerous local minima that also makes the optimization task by gradient-based algorithms complicated.

A second possibility is to use gradient-free 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 Covari-ance Matrix Adapatation Evolution Strategy (CMA-ES) 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 mixed-effects 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 mixed-effects 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 frequen-tist 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 log-likelihood) [51]. This takes advantage of the fact that, for classical calibration problems, assuming a Gaussian discrepancy ξ in equation (1) makes the negative log-likelihood 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 mixed-effects models. Remember that the objective of the mixed-effects 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 :

fΘ|Π(θi|Π)=c(F1(θi1|Π1),,Fd(θid|Πd)|Πcop)                               ×l=1dfl(θil|Πl),(22)

with Πcop the copula parameters, Πl; the parameters of the lth marginal of distribution fΘ|Π (so the marginal marginal distribution of θil (i, l) ∈ 〚1,n〛 × 〚1,d〛) with probability density function fl(·) and cumulative density function Fl; (·) 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 positive-definite. Sylvester’s criterion [53] states that if the determinant of the sub-matrices extracted from the full correlation matrix by deleting a row and the symmetric column is strictly positive, then the full correlation matrix remains positive-definite. This condition is full-filled if these sub-matrices are positive-definite. 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 Sq ⊂ 〚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 qth step are re-estimated, 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 log-likelihood 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 log-likelihood 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, sq–1 the number of common parameters from step q – 1 to q and [ Πl,q1 ]l1,sq1 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 [ Πl,q1 ]l1,sq1 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 second-order derivative of the log-likelihood. Its evaluation is difficult and demands appropriate estimation techniques [54,55]. The observed FIM, in fact the second-order derivative of the log-likelihood 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 hyper-parameters of the method and ought to be carefully chosen.

thumbnail Fig. 3

Recap of the mixed-effects calibration procedure. Modeling stages are highlighted in red, those that concern the calibration itself are highlighted in blue.

thumbnail 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, uq refers to the integer labeling the subgroup of parameters that will be re-estimated in the current iteration. NG 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.

thumbnail 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:

(ε11ε222ε12)=S__(σ11σ22σ12)=(S11S210S12S22000S66)   (σ11σ22σ12),(23)

with S__ the compliance tensor. S__ is a symmetric tensor so S12 = S21. For practical interpretations, the coefficients of the compliance tensor can be related to the elastic moduli, that is to say the longitudinal modulus E11, the transverse modulus E22, the shear modulus G12 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 S11 is changed and expresses as follows:

S11L=S1L+(S110S1L)ε0L(S110S1L)σ11+ε0L,(24)

with S110 the elastic linear compliance (module at the origin), S1L the asymptotic compliance and ε0L 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 S12L should also be updated and refers to S12L:=v12S11L. 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 S1T<S110 and S1C<S110. With this parameterization, the model is not identifiable. In fact, it is over-parameterized, i.e. there are several combinations of both S1L and ε0L that give the same output as shown by Germain [57]. Thus, ε0L fixed to a usual nominal value, here 0.005 [57]. Finally, there are 6 parameters to be calibrated S1T, S1C, S110, S22,v12 and S66.

In theory, the compliance matrix S__ should be positive-definite [25] To enforce this property on the individual compliance tensor, two requirements have to be met: the strict positivity of S66 and the positive-definiteness 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 positive-definitness 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 μ=(μS110   μS1T   μS1C   μv12   μS22   μS66) and covariance matrix

Σ=S110S1TS1Cv12S22S66(S110S1TS1Cv12S22S66𝕧(S110)*****Cov(S110,S1T)𝕧(S1T)****Cov(S110,S1C)Cov(S1T,S1C)𝕧(S1C)***Cov(S110,v12)Cov(S1T,v12)Cov(S1C,v12)𝕧(v12)**Cov(S110,S22)Cov(S1T,S22)Cov(S1C,S22)Cov(S22,v12)𝕧(S22)*Cov(S110,S66)Cov(S1T,S66)Cov(S1C,S66)Cov(v12,S66)Cov(S22,S66)𝕧(S66)) (25)

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, S110, S1T, S1C and v12 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 S22, the same applies for the linear part of the difference between the longitudinal and transverse responses for the ±45° laminate, mostly driven by S66 [25,58,59]. Remember that, with mixed-effects, 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 S110, S1T, S1C, v12 and S22 cannot be estimated with the usual elementary experiments, and the same applies to those between S110, S1T, S1C v12 and S66 or between S22 and See. In fact, from these elementary experiments, only the joint distribution of S110, S1T, S1C, v12 and the marginals of S66 and S22 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 S110, S1T, S1C v12 and 2 for both marginals of S22 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 Ai and Bi. For instance, with the mean vector defined earlier, one possible choice would be Ai = Bi equal to the identity matrix of dimension 6 to which the two last rows would be deleted (and Ai would have thus 4 rows with 6 columns) for the UD0° experiments, Ai = Bi = (0, 0, 0, 0, 1, 0) for the 90° experiments and Ai = Bi = (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 S110, S1T, S1C, v12 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 S66 and S22 are carried out from experiments with univariate observations that exhibit a linear behavior. Thus, these cases only refer to univariate linear mixed-effects 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 S110, S1T, S1C and v12.

One first remark that can be made is that S1C only expresses when compressive loads are applied to the sample, and S1T 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 S110, S1T, S1C and ν12 is parameterized by 13 and not 14 population parameters. S110 and v12 are activated both in tension and compression. Note that v12 only expresses on the transverse responses, either in tension and compression. Here, there are shared parameters for two types of experiments, S110 and v12. Design matrices apply here. For instance, for the tension experiments, matrices Ai and Bi 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 S110, S1T and v12. It is thus possible to calibrate the parameters highlighted in red in equation (26):

μ=(μS110μS1CμS1Tμv12) =S110S1CS1Tv12(S110S1CS1Tv12𝕍(S110)***Cov(S1C,S110)𝕍(S1C)**Cov(S1T,S110)0𝕍(S1T)*Cov(v12,S110)Cov(v12,S1C)Cov(v12,S1T)𝕍(v12)).(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 S110, S1C and v12. The application of previous remarks shows that it is possible to calibrate the parameters highlighted in blue in equation (26):

μ=(μS110μS1CμS1Tμv12) =S110S1CS1Tv12(S110S1CS1Tv12𝕍(S110)***Cov(S1C,S110)𝕍(S1C)**Cov(S1T,S110)0𝕍(S1T)*Cov(v12,S110)Cov(v12,S1C)Cov(v12,S1T)𝕍(v12)).(27)

Note that μS110, μυ12, 𝕍(S110), Cov( S110, v12) and 𝕍(v12) 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 S110, S1T and v12 by limiting the parameters of the joint distribution of S110,v12 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 fixed-point 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 S110, S1T, S1C, v12, S22 and S66 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 S22 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 𝒪(d2) 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 mixed-effects, the residual model (here the ωis) 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 mixed-effects 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.

Table 1

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 mixed-effects framework

To test the ability of the proposed methodology to calibrate the elastic model with mixed-effects, 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 mixed-effects 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:

ε(Π^)   :=(ε(Π^1),,ε(Π^)) with ε(Π^u)              :=| Πexact,uΠ^u |Πexact,u,u1,,(28)

noting 𝓁 the number of distribution parameters and the KL-divergence [60], KL(fΘ|Πexact, fΘ|Π^), between the exact fΘ|Πexact and the estimated distribution fΘ|Π^ :

KL (fΘ|Πexact,fΘ|Π^):=dfΘ|Πexact(θ)In (fΘ|Πexact(θ)fΘ|Π^(θ)) dθ.(29)

In general, it is computed with coarse Monte-Carlo. 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 Kullback-Leibler (KL) divergence can be estimated in closed form and refers to [61]

KL (fΘ|Πexact,fΘ|Π^)=12( tr(Σ^1Σexact)+ (μexactμ^)Σ^1(μexactμ^)+In (| Σ^ || Σexact |)d ).(30)

The KL-divergence 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 KL-divergence 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 θ^i or estimated output (ti, θ^i), the averaged relative error on the individual parameters e(θ^1,… θ^n) is defined as:

e(θ^1,,θ^n):=1ni=1nei(θ^i) with  ei(θ^i)=(ei(θ^i1),   ,   ei(θ^id))ei(θ^ia)=| θi,exactaθ^ia |θi,exactaa1,d,(31)

and the average error d(θ^1,…, θ^n) between the exact and estimated strains for each coordinate

d(θ^1,,θ^n):=(d1(θ^1,,θ^n),   ,   dm(θ^1,,θ^n)) with  dk(θ^1,,θ^n):=1ni=1n1Ni εi,exactkk(ti,θ^i) 22k1,m,(32)

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 fΘ|Πexact. The above indicators are averaged over the repetitions and their dispersion is characterized by the coefficient of variation COV(X) defined for any non-zero quantity X as COV(X)=V^(X)M^(X) with M^(X) the empirical mean of X and V^ (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 log-likelihood, it is decided here to assign to each individual and measure the same covariance matrix of the discrepancy term:

i1,n,  Ωi=Ω.

Table 2

Exact values of the parameters distribution.

4.4 Generating virtual data

The exact means and standard-deviations of the parameters for the joint distribution fΘ|Πexact 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]. Standard-deviations of S110 and v12 are chosen to fit with the values of dispersion given by Rollet [62]. Several correlations scenarios are considered: independence (between S1C and S1T), positive correlation : ρ(S110, S1T) = 0.809, ρ(v12, S1C) = 0.730 and negative correlations ρ(S1T,v12) = −0.500, ρ(S110, v12) = −0.741, ρ(S1C, S110) = −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 fΘ|Πexact 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:

εi,noisyj=εi,exactj×(12+τζ) with ζ𝒩(02,I2) and τ=0.02,(33)

with 12 the unitary vector of ℝ2, 02 the zero vector of ℝ2 and I2 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 mixed-effects 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 S110,v12 and S1T 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.

thumbnail Fig. 6

Example of a set of 50 synthetic stress-strain 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).

Table 3

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 Covariance-Matrix-Adaptation-Evolution-Strategy [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 standard-deviations 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 positive-definite character of the correlation matrix. Those on the variance of the residual models are determined by trial-error. The initialization point is taken as 0.65×1 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 log-likelihood 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 Nelder-Mead 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 multi-start 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µ].

Table 4

Averaged calibrated mean, standard-deviations and correlations Π1,mean and averaged error εmean(Π^1) 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 S110, S1T 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 standard-deviation 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 positive-definite 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, standard-deviations 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 standard-deviation 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 standard-deviations of both S1T 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 S110 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 S110 mainly expresses at the beginning of the strain-stress curve, where the residuals have the lowest values, while S1T 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 S110 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 S110, S1T 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 ςexact:=diag(Σexact)μexact 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 S110 . 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 S1T. Thus, within [μS1103σS110 ,μS110+3σS110], (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 standard-deviation 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 S110, S1T and v12 either at the population or individual levels. Consequently, the first stage of the methodology is completed.

Table 5

Relative errors on the estimation of the individual parameters e(θ^1,…, θ^n) 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.

Table 6

Error in model space d(θ^1,…, θ^n) for each strain component.

4.5.3 Calibration of the joint distribution of S110, S1C and v12

In the second stage of the methodology, the objective is to identify the joint distribution of S110, S1C and v12. Furthermore, the parameters of the joint distribution of v12 and S110:μv12, μS110, 𝕍(v12), 𝕍(S110) and ρ(v12, S110), 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, ρ(v12, S110 ) , 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 S1C (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 S1C 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 S110 and v12. 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 standard-deviations appear to be flawed from their initial estimation (Tab. 4), in particular for the standard-deviation of S110 (the averaged error increases by 30%). Indeed, both in terms of averaged values or averaged error on the population parameters. The correlation between S110 and v12 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 Kullback-Leibler (KL) divergence between the exact joint distribution of S110, v12exact(S110, v12)) and the calibrated distribution (ƒcalibrated(S110, v12 )) 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 S1C (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 S1C 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 S110 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 standard-deviations appear to be flawed from their initial estimation (Tab. 4), in particular for the standard-deviation of S110 (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 S110 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 Kullback-Leibler divergence between the exact joint distribution of S110 , ν12exact(S110, ν12)) and the calibrated distribution (ƒcalibrated(S110, ν12)) in the two stages as depict Figure 7 for the 20 repetitions.

In Figure 7, the deterioration of the joint distribution of S110 and ν12 suggested by Table 7 seems to be confirmed. For instance, the median of the KL-divergences increases (multiplied by 3), and the ranges of variation of the KL-divergences are also increased by two (without consideration of the outliers). This seems mainly due to poorer estimation in the compression phase of the standard-deviation of S110 . Consequently, a deterioration of the estimation of joint distribution of S110 and ν12 can be observed in the second phase. Nevertheless, this deterioration remains limited for both the relative errors and the KL-divergences. In particular, even if they are larger in the second phase than in the first phase, the KL-divergences have the same orders of magnitude. Indeed, except for one repetition, the KL-divergences 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 S110 (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 ςexact:=diag(Σexact)μexact 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 S110, ν12 are flawed with higher errors, partly explaining the poorer calibration of their joint distribution highlighted above. This could come from the non-linearity of the strain-stress 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.

Table 7

Averaged calibrated mean, standard-deviations and correlations Π2,mean and averaged error εmean (Π^2) over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.

thumbnail Fig. 7

KL-divergences of the calibrated and exact joint distribution S110 and ν12 for the 20 repetitions in the two stages. Remember that the KL-divergence expresses a notion of distance between two distributions and the lower is its value, the more accurate the calibration of ƒΘ|Π is.

Table 8

Relative errors on the estimation of the individual parameters e(θ^1,…, θ^n) averaged over 20 repetitions of the calibration process with the spread ςexact. The corresponding coefficients of variation in % are indicated between brackets.

Table 9

Error in model space d(θ^1,…, θ^n) 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 S110, v12 ? In particular, is the estimation downgraded or improved between the two stages ? To answer this question, the KL-divergence between the exact joint distribution of S110, ν12exact(S110, ν12)) and the calibrated distribution (ƒcalibrated(S110, v12)) 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 KL-divergences 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 S110 and v12 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 S110 and v12 is estimated with increased accuracy when its calibration is carried out in the second phase rather than in the first phase, as show the KL-divergences in Figures 7 and 8. Indeed, though their medians are comparable in both situations, the range of variation of the KL-divergences (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 S110 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 σS110, 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 standard-deviation of v12 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 mixed-effects to estimate, with the help of design matrices, the 13 parameters of the joint distribution of S1T, S1C, S110 and ν12 (recall that Cov(S1T, S1C) 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 C1, 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 S1T, S1C, S110 and ν12, labeled C3. 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 C3 compared to cases C1 and C2 . This could come from the fact in case C3, 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 S110 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 trade-off 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 S110 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 KL-divergence is computed for two distributions: the joint distribution of S110 and v12 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 KL-divergences are similar in terms of values and range of variation, either for the joint distribution of S110 and v12 or for the complete joint distribution. Still, it is possible to observe the KL-divergences have significant range of variation (from 10−1 to over 108). This can be partly explained by the fact that the complete distribution is parameterized by 13 parameters (against 5 for the joint distribution of S110 and v12) whose 2 correlations exhibit large level of errors (in particular ρ(S1C,S110)). 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 C2 (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 C3 to case C1 and by 40% from case C3 to case C2, which can be considered as a significant gain in terms of computational cost. Still, the number of model evaluations remains high, about 108. This number can be partly explained by the choices and settings of the optimization algorithms used for the minimization of both the opposite log-likelihood and function gi(·). 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 C1 and C2 ) 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 C3, 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 bi 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 C1 with the compression experiments and case C2 with the tension experiments to their counterpart of case C3.

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 C1 and C2 . 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 mixed-effects 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 test-case. 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 fixed-point approach at the expense of a larger computational costs. Still, this test-case 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.

thumbnail Fig. 8

KL-divergences of the calibrated and exact joint distribution S110 and ν12 for the 20 repetitions in the two stages. Remember that the KL-divergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒΘ∣Π is.

Table 10

Averaged marginal parameters (1,mean from the first phase and 2,mean from the second phase), alongside the averaged error (εmean(Π^2) and εmean(Π^2)) 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.

Table 11

Averaged marginal parameters correlation parameter ρ^(S110,v12)mean,1/2 in both phases and averaged error εmean(ρ^(S110,v12))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.

Table 12

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 C1 and C2, the error on the parameters of the joint distribution of S110 and v12 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.

thumbnail Fig. 9

For the 17th repetition, PDF of S110 and ν12 and the corresponding individual parameters for the calibration from case C1 (with the compression experiments in the second phase), from case C2 (with the tension experiments in the second phase) and considering both types of experiments (case C3).

thumbnail Fig. 10

On the left, KL-divergences between the calibrated and exact joint distribution of S110 and v12 and on the right, KL-divergences between the calibrated and exact joint full distributions for the 20 repetitions and 3 cases. Remember that the KL-divergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒΘ|Π is.

thumbnail 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 C1 and C2, the costs of both stages are combined together.

thumbnail Fig. 12

On the left, total number of model evaluations for the estimation of the individual parameters with the compression experiments in case C3 and on the second stage of case C1. On the right, total number of model evaluations for the estimation of the individual parameters with the tension experiments in case C3 and on the second stage of case C2.

5 Conclusion

This article presents a sequential calibration procedure that aims to calibrate complex material constitutive models compliant with the mixed-effects 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, well-chosen 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 test-case 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 trade-offs 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 mixed-effects 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

  1. U.D.O. Defense, Technomic Publishing Company, Materials Sciences Corporation, American Society for Testing, and Materials, Composite Materials Handbook-MIL 17, Vol. III: Materials Usage, Design, and Analysis, The Composite Materials Handbook-MIL 17. Taylor & Francis, 1999 [Google Scholar]
  2. A. Tarantola, Inverse Problem Theory and Methods for Model Parameter Estimation, Society for Industrial and Applied Mathematics, 2005. [CrossRef] [Google Scholar]
  3. L. Wasserman, All of Statistics: A Concise Course in Statistical Inference, Springer Publishing Company, Incorporated, 2010 [Google Scholar]
  4. A. Antoniadis, R. Carmona, Regression non lineaire et applications. Economie et statistiques avan¸c´ees, Economica, 1992 [Google Scholar]
  5. 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]
  6. 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]
  7. E. Demidenko, Mixed Models. Theory and Applications with R, John Wiley & Sons, 2013 [Google Scholar]
  8. M. Lavielle, Mixed Effects Models for the Population Approach: Models, Tasks, Methods and Tools, CRC press, 2014 [CrossRef] [Google Scholar]
  9. P. Congdon, Bayesian Statistical Modelling, John Wiley & Sons, 2007 [Google Scholar]
  10. C. Laboulfie, M. Balesdent, L. Brevault, F.-X. Irisarri, J.-F. Maire, S. Da Veiga, R. Le Riche, Calibration of material model parameters using mixed-effects booktitle = 25ème Congres Français de Mecanique, model, 2021, pp. 258–295 [Google Scholar]
  11. J. Lemaitre, J.-L. Chaboche, Mechanics of solid materials, Cambridge University Press, 1994 [Google Scholar]
  12. J. Lemaitre, R. Desmorat. Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures. Springer Science & Business Media, 2006 [Google Scholar]
  13. 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]
  14. R. Schapery, On the characterization of nonlinear viscoelastic materials, Polymer Eng. Sci. 9, 295–310 (1969) [CrossRef] [Google Scholar]
  15. J.-L. Chaboche, A review of some plasticity and viscoplasticity constitutive Theories, Int. J. Plasticity 24, 1642–1693 (2008) [CrossRef] [Google Scholar]
  16. B. Miled, I. Doghri, L. Delannay, Coupled viscoelastic-viscoplastic modeling of homogeneous and isotropic polymers: numerical algorithm and analytical solutions. Comput. Methods Appl. Mech. Eng. 200, 3381–3394 (2011). [Google Scholar]
  17. A. Krairi, I. Doghri, A thermodynamically-based constitutive model for thermoplastic polymers coupling viscoelasticity, viscoplasticity and ductile damage, Int. J. Plasticity 60, 163–181 (2014) [Google Scholar]
  18. F. Praud, G. Chatzigeorgiou, J. Bikard, F. Meraghni, Phenomenological multi-mechanisms constitutive modelling for thermoplastic polymers, implicit implementation and experimental validation, Mech. Mater. 114, 9–29 (2017) [Google Scholar]
  19. C. Hochard, P.-A. Aubourg, J.-P. Charles, Modelling of the mechanical behaviour of woven-fabric cfrp laminates up to failure, Composites Sci. Technol. 61, 221–230 (2001) [CrossRef] [Google Scholar]
  20. P. Ladeveze, E. LeDantec, Damage modelling of the elementary ply for laminated composites, Composites Sci. Technol. 43, 257–267 (1992) [CrossRef] [Google Scholar]
  21. 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]
  22. F. Praud, G. Chatzigeorgiou, Y. Chemisky, F. Meraghni, Hybrid micromechanical-phenomenological modelling of anisotropic damage and anelasticity induced by microcracks in unidirectional composites, Composite Struct. 182, 223–236 (2017) [CrossRef] [Google Scholar]
  23. 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]
  24. A. Forrester and A. Sobester, A.and Keane. Engineering Design via Surrogate Modelling – A Practical Guide, Wiley, 2008 [CrossRef] [Google Scholar]
  25. J.-M. Berthelot, Materiaux Composites-Comportement Mecanique et Analyse des Structures, 2012. [Google Scholar]
  26. M. Gallagher, J. Doherty, Parameter estimation and uncertainty analysis for a watershed model, Environ. Model. Softw. 22, 1000–1020 (2007) [CrossRef] [Google Scholar]
  27. M. Kennedy, A. O’Hagan, Bayesian calibration of computer models, J. Roy. Stat. Soc. B (Stat. Methodol.) 63, 425–464 (2001) [CrossRef] [Google Scholar]
  28. K. Levenberg, A method for the solution of certain nonlinear problems in least Squares, Q. Appl. Math. 2, 164–168 (1944) [CrossRef] [Google Scholar]
  29. D.W. Marquardt, An algorithm for least-squares estimation of nonlinear parameter, J. Soc. Ind. Appl. Math. 11, 431–441 (1963) [CrossRef] [Google Scholar]
  30. R. Pintelon, J. Schoukens, System Identification, John Wiley & Sons, Ltd, 2012 [CrossRef] [Google Scholar]
  31. 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]
  32. J. Pinheiro, D. Bates, Mixed-Effect Models in S and S-plus, Vol. 96, 2002 [Google Scholar]
  33. 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]
  34. S.L. Beal, L.B. Sheiner, Estimating population kinetics, Crit. Rev. Biomed. Eng. 8, 195–222 (1982) [Google Scholar]
  35. D. Blei, A. Kucukelbir, J. McAuliff, Variational inference: a review for statisticians, J. Am. Stat. Assoc. 112, 859–877 (2017) [CrossRef] [Google Scholar]
  36. J. Wakefield, The Bayesian analysis of population pharmacokinetic models, J. Am. Stat. Assoc. 91, 62–75 (1996) [CrossRef] [Google Scholar]
  37. P. Damlen, J. Wakefield, S. Walker, Gibbs sampling for Bayesian non-conjugate and hierarchical models by using auxiliary variables, J. Roy. Stat. Soc. B (Stat. Methodol.) 61, 331–344 (1999) [CrossRef] [Google Scholar]
  38. M. Davidian, D.M. Giltinan, Nonlinear Models for Repeated Measurement Data, Chapman & Hall/CRC Monographs on Statistics & Applied Probability, Taylor & Francis, 1995 [Google Scholar]
  39. R. Fisher, K. Pearson, On an Absolute Criterion for Fitting Frequency Curves. Messenger of Mathematics, 1911 [Google Scholar]
  40. D. Hall, M. Clutter, Multivariate multilevel nonlinear mixed-effects models for timber yield predictions, Biometrics 60, 16–24 (2004) [CrossRef] [MathSciNet] [PubMed] [Google Scholar]
  41. 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]
  42. N. Metropolis, S. Ulam, The monte carlo method, J. Am. Stat. Assoc. 44, 335–341 (1949) [CrossRef] [PubMed] [Google Scholar]
  43. D. Bates, D. Hamilton, D. Watts, Calculation of intrinsic and parameter-effects curvatures for nonlinear regression models. Commun. Stat. Simul. Comput. 12, 469–477 (1983) [CrossRef] [Google Scholar]
  44. N. Hansen, S. Muller, P. Koumoutsakos, Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES), Evol. Comput. 11, 1–18 (2003) [CrossRef] [PubMed] [Google Scholar]
  45. W. Oberkampf, C. Roy, Verification and Validation in Scientific Computing, Cambridge University Press, 2010 [CrossRef] [Google Scholar]
  46. L. Schwer, Guide for Verification and Validation in Computational Solid Mechanics, ASME, 2006 [Google Scholar]
  47. 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]
  48. 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]
  49. A. Urbina, T. Paez, T. Consulting, A bayes network approach to uncertainty quantification in hierarchically-developed computational models, Int. J. Uncertain. Quantif. 2, 173–193 (2012) [CrossRef] [MathSciNet] [Google Scholar]
  50. 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. ASCE-ASME J. Risk Uncertain. Eng. Syst. B: Mech. Eng. 8, 011111 (2022) [Google Scholar]
  51. T. Kimand, B. Youn, Identifiability-based model decomposition for hierarchical Calibration, Struct. Multidiscipl. Optim. 60, 1801–1811 (2019) [Google Scholar]
  52. M. Sklar, Fonctions de repartition a n dimensions et leurs marges. in Annales de l’ISUP, Vol. 8, 1959, pp. 229–231 [Google Scholar]
  53. G. Gilbert, Positive definite matrices and sylvester’s criterion, Am. Math. Monthly 98, 44–46 (1991) [CrossRef] [Google Scholar]
  54. F. Mentre, A. Mallet, D. Baccar, Optimal design in random-effects regression models, Biometrika 84, 429–442 (1997) [CrossRef] [MathSciNet] [Google Scholar]
  55. 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]
  56. 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]
  57. 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 Paris-Saclay, January 2020 [Google Scholar]
  58. I. Daniel, O. Ishai, Engineering Mechanics of Composite Materials, Oxford university press New York, 2006 [Google Scholar]
  59. J. Tomblinand, Y. Ng, S. Keshavanarayana, Material qualification and equivalency for polymer matrix composite material systems, Technical report, 2001 [Google Scholar]
  60. S. Kullback, R. A. Leibler, On information and sufficiency, Ann. Math. Stat. 22, 79–86 (1951) [CrossRef] [Google Scholar]
  61. L. Pardo, Statistical Inference based on Divergence Measures, CRC Press, 2018 [CrossRef] [Google Scholar]
  62. Y. Rollet, Vers une Maîtrise des Incertitudes en Calcul des Structures Composites (in French), PhD thesis, 2007 [Google Scholar]
  63. 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érard-Marchant, 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]
  64. M. Baudin, A. Dutfoy, B. Iooss, A.-L. Popelin, Open-TURNS: 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]
  65. S.J. Johnson, The nlopt nonlinear-optimization package [Google Scholar]
  66. 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çois-Xavier Irisarri, Jean-François Maire, Sebastien Da Veiga, Rodolphe Le Riche, Sequential calibration of material constitutive model using mixed-effects calibration, Mechanics & Industry 24, 32 (2023)

All Tables

Table 1

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.

Table 2

Exact values of the parameters distribution.

Table 3

Search space for the population parameters.

Table 4

Averaged calibrated mean, standard-deviations and correlations Π1,mean and averaged error εmean(Π^1) over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.

Table 5

Relative errors on the estimation of the individual parameters e(θ^1,…, θ^n) 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.

Table 6

Error in model space d(θ^1,…, θ^n) for each strain component.

Table 7

Averaged calibrated mean, standard-deviations and correlations Π2,mean and averaged error εmean (Π^2) over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.

Table 8

Relative errors on the estimation of the individual parameters e(θ^1,…, θ^n) averaged over 20 repetitions of the calibration process with the spread ςexact. The corresponding coefficients of variation in % are indicated between brackets.

Table 9

Error in model space d(θ^1,…, θ^n) for each strain component.

Table 10

Averaged marginal parameters (1,mean from the first phase and 2,mean from the second phase), alongside the averaged error (εmean(Π^2) and εmean(Π^2)) 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.

Table 11

Averaged marginal parameters correlation parameter ρ^(S110,v12)mean,1/2 in both phases and averaged error εmean(ρ^(S110,v12))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.

Table 12

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 C1 and C2, the error on the parameters of the joint distribution of S110 and v12 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

thumbnail Fig. 1

Difference between the population-based and the classical approaches. On the left, a), usual methods with a single parameter vector for all individuals. On the right, b), population-based approach in the mixed-models effects framework with fΘ the distribution.

In the text
thumbnail 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
thumbnail Fig. 3

Recap of the mixed-effects calibration procedure. Modeling stages are highlighted in red, those that concern the calibration itself are highlighted in blue.

In the text
thumbnail 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, uq refers to the integer labeling the subgroup of parameters that will be re-estimated in the current iteration. NG refers to the number of calibration subproblems. The process stops when all calibration subproblems are solved.

In the text
thumbnail 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
thumbnail Fig. 6

Example of a set of 50 synthetic stress-strain 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
thumbnail Fig. 7

KL-divergences of the calibrated and exact joint distribution S110 and ν12 for the 20 repetitions in the two stages. Remember that the KL-divergence expresses a notion of distance between two distributions and the lower is its value, the more accurate the calibration of ƒΘ|Π is.

In the text
thumbnail Fig. 8

KL-divergences of the calibrated and exact joint distribution S110 and ν12 for the 20 repetitions in the two stages. Remember that the KL-divergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒΘ∣Π is.

In the text
thumbnail Fig. 9

For the 17th repetition, PDF of S110 and ν12 and the corresponding individual parameters for the calibration from case C1 (with the compression experiments in the second phase), from case C2 (with the tension experiments in the second phase) and considering both types of experiments (case C3).

In the text
thumbnail Fig. 10

On the left, KL-divergences between the calibrated and exact joint distribution of S110 and v12 and on the right, KL-divergences between the calibrated and exact joint full distributions for the 20 repetitions and 3 cases. Remember that the KL-divergence expresses a notion of distance between two distributions and the lower the more accurate the calibration of ƒΘ|Π is.

In the text
thumbnail 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 C1 and C2, the costs of both stages are combined together.

In the text
thumbnail Fig. 12

On the left, total number of model evaluations for the estimation of the individual parameters with the compression experiments in case C3 and on the second stage of case C1. On the right, total number of model evaluations for the estimation of the individual parameters with the tension experiments in case C3 and on the second stage of case C2.

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.