# Sequential calibration of material constitutive model using mixed-effects calibration

Clément Laboulfie, Mathieu Balesdent, Loïc Brevault, François-Xavier Irisarri, Jean-François Maire, Sebastien da Veiga, Rodolphe Le Riche

► To cite this version:

Clément Laboulfie, Mathieu Balesdent, Loïc Brevault, François-Xavier Irisarri, Jean-François Maire, et al.. Sequential calibration of material constitutive model using mixed-effects calibration. Mechanics & Industry, 2023, 24, pp.32. 10.1051/meca/2023027. hal-04231804

HAL Id: hal-04231804

https://hal.science/hal-04231804

Submitted on 6 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sequential calibration of material constitutive model using mixed-effects calibration

Clément Laboulfie, Mathieu Balesdent, Loïc Brevault, François-Xavier Irisarri, Jean-François Maire, Sebastien da Veiga, Rodolphe Le Riche

► To cite this version:

Clément Laboulfie, Mathieu Balesdent, Loïc Brevault, François-Xavier Irisarri, Jean-François Maire, et al.. Sequential calibration of material constitutive model using mixed-effects calibration. Mechanics & Industry, 2023, 24, pp.32. 10.1051/meca/2023027. hal-04231804

HAL Id: hal-04231804

https://hal.science/hal-04231804v1

Submitted on 6 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Sequential calibration of material constitutive model using
mixed-effects calibration

Clément Laboulief<fnref n="1" />,*, Mathieu Balesdent<fnref n="2" />, Loïc Brevault<fnref n="2" />, François-Xavier Irisarri<fnref n="1" />, Jean-François Maire<fnref n="1" /> ID,
Sebastien Da Veiga<fnref n="3" />, and Rodolphe Le Riche<fnref n="4" /> ID

Received: 10 February 2023 / Accepted: 18 July 2023

**Abstract.** Identifying model parameters is nowadays intrinsically linked with quantifying the associated uncertainties. While classical methods allow to handle some types of uncertainties such as experimental noise, they are not designed to take into account the variability between the different test specimens, significant in particular for composites materials. The estimation of the impact of this intrinsic variability on the material properties can be achieved using population approaches where this variability is modeled by a probability distribution (e.g., a multivariate Gaussian distribution). The objective is to calibrate this distribution (or equivalently its parameters for a parametric distribution). Among the estimation methods can be found mixed-effects models where the parameters that characterize each replication are decomposed between the population averaged behavior (called fixed-effects) and the impact of material variability (called random-effects). Yet, when the number of model parameters or the computational time of a single run of the simulations increases (for multiaxial models for instance), the simultaneous, global identification of all the material parameters is difficult because of the number of unknown quantities to estimate and because of the required model evaluations. Furthermore, the parameters do not have the same influence on the material constitutive model depending for instance on the nature of the load (e.g., tension, compression). The method proposed in this paper enables to calibrate the model on multiple experiments. It decomposes the overall calibration problem into a sequence of calibrations, each subproblem allowing to calibrate the joint distribution of a subset of the model parameters. The calibration process is eased as the number as the number of unknown parameters is reduced compared to the full problem. The proposed calibration process is applied to an orthotropic elastic model with non linear longitudinal behavior, for a unidirectional composite ply made of carbon fibers and epoxy resin. The ability of the method to sequentially estimate the model parameters distribution is investigated. Its capability to ensure consistency throughout the calibration process is also discussed. Results show that the methodology allows to handle the calibration of complex material constitutive models in the mixed-effects framework.

**Keywords:** Mixed-effects models/model calibration/composite material/sequential calibration

# 1 Introduction

Following the increase of numerical facilities, virtual test-
ing is now widely used in order to partially substitute
numerical simulations to experimental campaigns because
of time and costs considerations. These numerical sim-
ulations 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 cali-
brated. Calibrating model parameters consists in finding

the parameters value that allow to best fit, in a specific
mathematical sense that should be specified, the exper-
imental 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 parame-
ters 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 com-
plexity of the manufacturing process. To characterize this
variability, international standards impose to perform test

<footnote n="1">¹ DMAS, ONERA, Université Paris Saclay, F-92322 Châtillon, France</footnote>

<footnote n="2">² DTIS, ONERA, Université Paris Saclay, F-91123 Palaiseau, France</footnote>

<footnote n="3">³ ENSAI (Campus de Ker Lann), 51 Rue Blaise Pascal 35172 Bruz</footnote>

<footnote n="4">⁴ LIMOS (CNRS, Mines Saint-Étienne & UCA), France</footnote>

<footnote>* e-mail: clement.laboulief@onera.fr</footnote>
Sequential calibration of material constitutive model using
mixed-effects calibration

Clément Laboulief<fnref n="1" />,*, Mathieu Balesdent<fnref n="2" />, Loïc Brevault<fnref n="2" />, François-Xavier Irisarri<fnref n="1" />, Jean-François Maire<fnref n="1" /> ID,
Sebastien Da Veiga<fnref n="3" />, and Rodolphe Le Riche<fnref n="4" /> ID

Received: 10 February 2023 / Accepted: 18 July 2023

**Abstract.** Identifying model parameters is nowadays intrinsically linked with quantifying the associated uncertainties. While classical methods allow to handle some types of uncertainties such as experimental noise, they are not designed to take into account the variability between the different test specimens, significant in particular for composites materials. The estimation of the impact of this intrinsic variability on the material properties can be achieved using population approaches where this variability is modeled by a probability distribution (e.g., a multivariate Gaussian distribution). The objective is to calibrate this distribution (or equivalently its parameters for a parametric distribution). Among the estimation methods can be found mixed-effects models where the parameters that characterize each replication are decomposed between the population averaged behavior (called fixed-effects) and the impact of material variability (called random-effects). Yet, when the number of model parameters or the computational time of a single run of the simulations increases (for multiaxial models for instance), the simultaneous, global identification of all the material parameters is difficult because of the number of unknown quantities to estimate and because of the required model evaluations. Furthermore, the parameters do not have the same influence on the material constitutive model depending for instance on the nature of the load (e.g., tension, compression). The method proposed in this paper enables to calibrate the model on multiple experiments. It decomposes the overall calibration problem into a sequence of calibrations, each subproblem allowing to calibrate the joint distribution of a subset of the model parameters. The calibration process is eased as the number as the number of unknown parameters is reduced compared to the full problem. The proposed calibration process is applied to an orthotropic elastic model with non linear longitudinal behavior, for a unidirectional composite ply made of carbon fibers and epoxy resin. The ability of the method to sequentially estimate the model parameters distribution is investigated. Its capability to ensure consistency throughout the calibration process is also discussed. Results show that the methodology allows to handle the calibration of complex material constitutive models in the mixed-effects framework.

**Keywords:** Mixed-effects models/model calibration/composite material/sequential calibration

# 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 exper-
imental 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 parame-
ters 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 com-
plexity of the manufacturing process. To characterize this
variability, international standards impose to perform test

<footnote n="1">¹ DMAS, ONERA, Université Paris Saclay, F-92322 Châtillon, France</footnote>

<footnote n="2">² DTIS, ONERA, Université Paris Saclay, F-91123 Palaiseau, France</footnote>

<footnote n="3">³ ENSAI (Campus de Ker Lann), 51 Rue Blaise Pascal 35172 Bruz</footnote>

<footnote n="4">⁴ LIMOS (CNRS, Mines Saint-Étienne & UCA), France</footnote>

<footnote>* e-mail: clement.laboulief@onera.fr</footnote>
repetitions [1] to determine the model parameters distri-
bution. In the classical approaches, model parameters are
calibrated independently on several specimens by mini-
mizing a fitting criterion [2] before inferring the model
parameters distribution by maximizing the joint likeli-
hood 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 observa-
tions and model predictions, the negative log-likelihood
of the model given the experiments [3,4]) does not con-
sider 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 vari-
ability 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 con-
trary to HBMs that are fully Bayesian, which usually
involves costly Markov Chain Monte Carlo (MCMC) sam-
pling. 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 analy-
ses 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 “ran-
dom” 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 com-
bination of a reference value (the average value given by
the manufacturer for a material batch) and of a devia-
tion due to the variability of the production process. The
parameters of this distribution can be determined by max-
imizing 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 elastic-
ity, viscosity, plasticity and damage [11]. Hence, various
types of modelings have been developed, for instance dam-
age [12], viscoelasticity [13,14], plasticity or viscoplasticity
[15], coupled viscoelasticity and viscoplasticity [16], cou-
pled 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 mod-
els for a single ply [19–22]). 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 like-
lihood 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 esti-
mation 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 pres-
ence of uncertainty, we will seek to take advantage of the
decomposition of the model to transform the overall cal-
ibration problem into a sequence of appropriate smaller
calibration subproblems, corresponding to the different
test configurations (e.g. different load profiles and stack-
ing sequences). In addition to decreasing the number of
parameters to be identified, this alleviates the compu-
tational costs as it avoids to take into account all the
phenomena at the same time.

To calibrate the material parameters, often, experi-
mental 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 decomposi-
tion 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 sequen-
tial procedure to calibrate material constitutive models
compliant with the mixed-effects framework. This proce-
dure 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 (trans-
formed 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 lim-
ited 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 intro-
duces 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.
repetitions [1] to determine the model parameters distri-
bution. In the classical approaches, model parameters are
calibrated independently on several specimens by mini-
mizing a fitting criterion [2] before inferring the model
parameters distribution by maximizing the joint likeli-
hood 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 observa-
tions and model predictions, the negative log-likelihood
of the model given the experiments [3,4]) does not con-
sider 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 vari-
ability 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 con-
trary to HBMs that are fully Bayesian, which usually
involves costly Markov Chain Monte Carlo (MCMC) sam-
pling. 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 analy-
ses 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 “ran-
dom” 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 com-
bination of a reference value (the average value given by
the manufacturer for a material batch) and of a devia-
tion due to the variability of the production process. The
parameters of this distribution can be determined by max-
imizing 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 elastic-
ity, viscosity, plasticity and damage [11]. Hence, various
types of modelings have been developed, for instance dam-
age [12], viscoelasticity [13,14], plasticity or viscoplasticity
[15], coupled viscoelasticity and viscoplasticity [16], cou-
pled 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 mod-
els for a single ply [19–22]). 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 like-
lihood 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 esti-
mation 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 pres-
ence of uncertainty, we will seek to take advantage of the
decomposition of the model to transform the overall cal-
ibration problem into a sequence of appropriate smaller
calibration subproblems, corresponding to the different
test configurations (e.g. different load profiles and stack-
ing sequences). In addition to decreasing the number of
parameters to be identified, this alleviates the compu-
tational costs as it avoids to take into account all the
phenomena at the same time.

To calibrate the material parameters, often, experi-
mental 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 decomposi-
tion 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 sequen-
tial procedure to calibrate material constitutive models
compliant with the mixed-effects framework. This proce-
dure 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 (trans-
formed 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 lim-
ited 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 intro-
duces 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 $\mathcal{F}(\cdot, \theta)$, the calibration aims at finding a set of model parameters labeled $\theta \in \mathbb{R}^d$ ($d$ standing for the number of parameters to be calibrated) which allows to mimic the experimental data. The random variable to $\theta$ is labeled $\Theta$ and $\theta$ 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 $\mathbf{Y}$ and $\mathbf{y} = (y_j)_{j \in [1, N]}$ its outcome, with $N$ the number of observations, considering one individual. Classically, the following decomposition applies:

$$ \mathbf{y} \sim \mathcal{F}(\cdot, \theta) + \xi, \quad (1) $$

where $\xi$ stands for the random vector of the errors which represents the experimental noise and the model bias. The outcomes of $\xi$ can be different from one measure to the other and are labeled $(\xi_j)_{j \in [1, N]}$. Equation (1) can be further detailed:

$$ y_j := \mathcal{F}(t_j, \theta) + \xi_j, \quad (2) $$

with $t_j$ the $j^{th}$ input measure, which is deterministic. A common practice of classical identification consists in minimizing, the sum of squares of the difference between the model prediction and the experimental data with respect to model parameters. This method is known as nonlinear least squares and the minimization is often conducted with the 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_{\Theta}$ whose

outcomes are the individual parameters $(\theta_i)_{i \in [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_{\Theta}$ whose outcomes are the individual parameters:

$$ \forall i \in [1, n], \quad \theta_i \sim f_{\Theta}. \quad (3) $$

Both the probability distribution $f_{\Theta}$ and its realizations $\theta_i$ are unknown and the aim is to determine them. In addition, if $f_{\Theta}$ is parametric (a Gaussian distribution for instance), let $\mathbf{\Pi}$ be its parameters and $f_{\Theta} = f_{\Theta|\mathbf{\Pi}}$. Identifying $f_{\Theta|\mathbf{\Pi}}$ is tantamount to determining $\mathbf{\Pi}$. Given $\mathbf{\Pi}$ and $\theta_i \sim f_{\Theta|\mathbf{\Pi}} \quad \forall i \in [1, n]$, the model output $\mathbf{y}_i$ can be written as

$$ \forall i \in [1, n], \quad \mathbf{y}_i \sim \mathcal{F}(\cdot, \theta_i) + \xi_i. \quad (4) $$

Without any other hypothesis, the outcomes of $\xi_i$ (labeled $(\xi_{ij})_{j \in [1, N_i]}$) are different for each individual and for each observation, with $N_i$ the number of observations of the $i^{th}$ individual. The global mixed-effects models for the $j^{th}$ output measure of the $i^{th}$ individual $y_{ij}$ reads as

$$ y_{ij} = \mathcal{F}(t_{ij}, \theta_i) + \xi_{ij}. \quad (5) $$

For the sake of simplicity and to ease the numerical derivations, $f_{\Theta|\mathbf{\Pi}}$ 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_{\Theta|\mathbf{\Pi}} := \mathcal{N}(\mu, \Sigma) \quad \text{and} \quad \mathbf{\Pi} := [\mu, \Sigma], \quad (6) $$

with $\mu \in \mathbb{R}^d$ the mean vector and $\Sigma \in \mathcal{M}_d(\mathbb{R})$ the covariance matrix. The second hypothesis is that for each individual and each measure, the discrepancy term is a Gaussian white noise (no bias, no correlation):

$$ \xi_{ij} \sim \mathcal{N}(0, \omega_{ij}^2), \quad (7) $$

with $\omega_{ij}$ the standard-deviation of the noise of the $j^{th}$ output measure of the $i^{th}$ individual. Furthermore, the noise is supposed to be homoscedastic, i.e. $\omega_{ij} = \omega_i \quad \forall j \in [1, N_i]$ [32]. This is a simplifying assumption but helps to limit the number of parameters to estimate. Finally, the vector of parameters to be calibrated is denoted $\Psi$:

$$ \Psi := (\mu, \Sigma, \Omega), $$
# 2 Population-based approaches for calibration

## 2.1 Formalization

For a given a material constitutive model $\mathcal{F}(\cdot, \theta)$, the calibration aims at finding a set of model parameters labeled $\theta \in \mathbb{R}^d$ ($d$ standing for the number of parameters to be calibrated) which allows to mimic the experimental data. The random variable to $\theta$ is labeled $\Theta$ and $\theta$ 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 $\mathbf{Y}$ and $\mathbf{y} = (y_j)_{j \in [1, N]}$ its outcome, with $N$ the number of observations, considering one individual. Classically, the following decomposition applies:

$$ \mathbf{y} \sim \mathcal{F}(\cdot, \theta) + \xi, \quad (1) $$

where $\xi$ stands for the random vector of the errors which represents the experimental noise and the model bias. The outcomes of $\xi$ can be different from one measure to the other and are labeled $(\xi_j)_{j \in [1, N]}$. Equation (1) can be further detailed:

$$ y_j := \mathcal{F}(t_j, \theta) + \xi_j, \quad (2) $$

with $t_j$ the $j$th input measure, which is deterministic. A common practice of classical identification consists in minimizing, the sum of squares of the difference between the model prediction and the experimental data with respect to model parameters. This method is known as nonlinear least squares and the minimization is often conducted with the 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_{\Theta}$ whose

outcomes are the individual parameters $(\theta_i)_{i \in [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_{\Theta}$ whose outcomes are the individual parameters:

$$ \forall i \in [1, n], \quad \theta_i \sim f_{\Theta}. \quad (3) $$

Both the probability distribution $f_{\Theta}$ and its realizations $\theta_i$ are unknown and the aim is to determine them. In addition, if $f_{\Theta}$ is parametric (a Gaussian distribution for instance), let $\Pi$ be its parameters and $f_{\Theta} = f_{\Theta|\Pi}$. Identifying $f_{\Theta|\Pi}$ is tantamount to determining $\Pi$. Given $\Pi$ and $\theta_i \sim f_{\Theta|\Pi}$ $\forall i \in [1, n]$, the model output $\mathbf{y}_i$ can be written as

$$ \forall i \in [1, n], \quad \mathbf{y}_i \sim \mathcal{F}(\cdot, \theta_i) + \xi_i. \quad (4) $$

Without any other hypothesis, the outcomes of $\xi_i$ (labeled $(\xi_{ij})_{j \in [1, N_i]}$) are different for each individual and for each observation, with $N_i$ the number of observations of the $i$th individual. The global mixed-effects models for the $j$th output measure of the $i$th individual $y_{ij}$ reads as

$$ y_{ij} = \mathcal{F}(t_{ij}, \theta_i) + \xi_{ij}. \quad (5) $$

For the sake of simplicity and to ease the numerical derivations, $f_{\Theta|\Pi}$ 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_{\Theta|\Pi} := \mathcal{N}(\mu, \Sigma) \text{ and } \Pi := [\mu, \Sigma], \quad (6) $$

with $\mu \in \mathbb{R}^d$ the mean vector and $\Sigma \in \mathcal{M}_d(\mathbb{R})$ the covariance matrix. The second hypothesis is that for each individual and each measure, the discrepancy term is a Gaussian white noise (no bias, no correlation):

$$ \xi_{ij} \sim \mathcal{N}(0, \omega_{ij}^2), \quad (7) $$

with $\omega_{ij}$ the standard-deviation of the noise of the $j$th output measure of the $i$th individual. Furthermore, the noise is supposed to be homoscedastic, i.e. $\omega_{ij} = \omega_i \forall j \in [1, N_i]$ [32]. This is a simplifying assumption but helps to limit the number of parameters to estimate. Finally, the vector of parameters to be calibrated is denoted $\Psi$:

$$ \Psi := (\mu, \Sigma, \Omega), $$
<figure><img src="image_2.png" /><figcaption>**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_{\Theta}$ the distribution.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

with $\Omega = \text{diag}(\omega_1^2, \dots, \omega_n^2)$. The population-based approaches seek $\Psi$ and provide an estimate of the individual parameters $(\theta_i)_{i \in [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 $\Psi$, 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 $\theta_i$ through fixed-effects $\beta$ and random effects $b_i$ as

$$ \theta_i = A_i\beta + B_i b_i, \qquad (8) $$

where $A_i, B_i$ are deterministic design matrices that belong to $M_d(\mathbb{R})$ with $d$ the number of model parameters. Design matrices other than $I_d$, the identity matrix of $M_d(\mathbb{R})$, 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 $\Pi$, but for a Gaussian $f_{\Theta|\Pi}$, the fixed-effects $\beta$ are equal to the population parameters average $\mu$ and the random-effects comply with $b_i \sim N(0_{\mathbb{R}^d}, \Sigma)$, $i \in [1, n]$. Furthermore, provided that the modeling proposed by the population exhibit a hierarchical structure, hierarchical-based techniques such as Hierarchical Bayesian models
<figure><img src="image_2.png" /><figcaption>**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_{\Theta}$ the distribution.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>

with $\Omega = \text{diag}(\omega_1^2, \dots, \omega_n^2)$. The population-based approaches seek $\Psi$ and provide an estimate of the individual parameters $(\theta_i)_{i \in [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 $\Psi$, 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 $\theta_i$ through fixed-effects $\beta$ and random effects $b_i$ as

$$ \theta_i = A_i\beta + B_i b_i, \qquad (8) $$

where $A_i, B_i$ are deterministic design matrices that belong to $M_d(\mathbb{R})$ with $d$ the number of model parameters. Design matrices other than $I_d$, the identity matrix of $M_d(\mathbb{R})$, 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 $\Pi$, but for a Gaussian $f_{\Theta|\Pi}$, the fixed-effects $\beta$ are equal to the population parameters average $\mu$ and the random-effects comply with $b_i \sim N(0_{\mathbb{R}^d}, \Sigma)$, $i \in [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 $\Pi$ (to which may be added $\Omega$, thus getting $\Psi$) 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 $\Pi$, in particular for the covariance matrix $\Sigma$. 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.

## 2.2 Likelihood of the mixed-effects

The calibration is often achieved by maximizing the likelihood of $\Psi$, $\mathcal{L}(\Psi) := f(y_1, \dots, y_n|\Psi)$ (where $f$ is a generic letter for probability density functions or PDFs), even if other methods can be found to estimate $\Psi$ [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 $\mathcal{L}_i(\Psi) := f(y_i|\Psi)$,

$$ \mathcal{L}(\Psi) := \prod_{i=1}^{n} \mathcal{L}_i(\Psi). \qquad (9) $$

Because the $\theta_i$ are not observed, the likelihood of the $i^{th}$ individual $\mathcal{L}_i(\Psi)$ is the integral of the marginal likelihood (i.e., the density of the output data $y_i$ given individual parameters $\theta_i$ and parameters $\Psi$) $f(y_i|\theta_i, \Pi, \omega_i)$ with

respect to all possible $\theta_i$ over $\mathbb{R}^d$:

$$ \mathcal{L}_i(\Psi) := \int_{\mathbb{R}^d} f(y_i|\theta_i, \Pi, \omega_i) f_{\Theta|\Pi}(\theta_i|\Pi) d\theta_i. \qquad (10) $$

In Equation (10), $f_{\Theta|\Pi}(\theta_i|\Pi)$ 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) $\hat{\Psi}$ is the result of the following maximization problem over the set of all the possible parameters $\Gamma_{\Psi}$:

$$ \hat{\Psi} := \arg\max_{\Psi \in \Gamma_{\Psi}} \mathcal{L}(\Psi). \qquad (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 $m_i$ the number of output measures of the $i^{th}$ individual and $k \in [1, m_i]$ the corresponding index. In the following, all individuals are supposed to have the same number of output measures, so $m_i = m$ $\forall i \in [1, n]$. Moreover, all measures are supposed to have the same number of observation points $N_i$. The main difference with Section 2.1 is that the model output of each observation point is no longer a scalar $y_{ij}$ but a vector $y_{ij} \in \mathbb{R}^m$. The same occurs for the discrepancy term, now labeled $\xi_{ij}$.

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

$$ \left\{ \begin{array}{l} \forall (i, j) \in [1, n] \times [1, N_i], y_{ij} = \mathcal{F}(t_{ij}, \theta_i) + \xi_{ij}, \quad (12a) \\ \forall i \in [1, n], \theta_i \sim f_{\Theta|\Pi}, \quad (12b) \\ \forall (i, j) \in [1, n] \times [1, N_i], \xi_{ij} \sim \mathcal{N}(0_{\mathbb{R}^m}, \Omega_{ij}), \quad (12c) \end{array} \right. $$

with $0_{\mathbb{R}^m}$ the zero vector of $\mathbb{R}^m$ and $\Omega_{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:

$$ \forall i \in [1, n], \forall (j_1, j_2) \in [1, N_i]^2, j_1 \neq j_2 $$
(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 frequentist and Bayesian frameworks. Note that in the Bayesian framework, both methods turn out to be equivalent. They will be here the chosen technique, mainly because of computational concerns in order to be able to perform the necessary computations. With 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.

## 2.2 Likelihood of the mixed-effects

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

$$L(Ψ) := \prod_{i=1}^{n} L_i(Ψ). \qquad (9)$$

Because the $θ_i$ are not observed, the likelihood of the $i^{th}$ individual $L_i(Ψ)$ is the integral of the marginal likelihood (i.e., the density of the output data $y_i$ given individual parameters $θ_i$ and parameters $Ψ$) $f(y_i|θ_i, Π, ω_i)$ with

respect to all possible $θ_i$ over $R^d$:

$$\mathcal{L}_i(Ψ) := \int_{\mathbb{R}^d} f(y_i|θ_i, Π, ω_i) f_{\Theta|Π}(θ_i|Π) dθ_i. \qquad (10)$$

In Equation (10), $f_{\Theta|Π}(θ_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) $\hat{Ψ}$ is the result of the following maximization problem over the set of all the possible parameters $Γ_Ψ$:

$$\hat{Ψ} := \arg\max_{Ψ \in Γ_Ψ} \mathcal{L}(Ψ). \qquad (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 $m_i$ the number of output measures of the $i^{th}$ individual and $k \in [1, m_i]$ the corresponding index. In the following, all individuals are supposed to have the same number of output measures, so $m_i = m$ $\forall i \in [1, n]$. Moreover, all measures are supposed to have the same number of observation points $N_i$. The main difference with Section 2.1 is that the model output of each observation point is no longer a scalar $y_{ij}$ but a vector $y_{ij} \in R^m$. The same occurs for the discrepancy term, now labeled $ξ_{ij}$.

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

$$\left\{ \begin{array}{l} \forall (i, j) \in [1, n] \times [1, N_i], y_{ij} = \mathcal{F}(t_{ij}, \theta_i) + \xi_{ij}, \quad (12a) \\ \forall i \in [1, n], \theta_i \sim f_{\Theta|\Pi}, \quad (12b) \\ \forall (i, j) \in [1, n] \times [1, N_i], \xi_{ij} \sim \mathcal{N}(0_{\mathbb{R}^m}, \Omega_{ij}), \quad (12c) \end{array} \right.$$

with $0_{\mathbb{R}^m}$ the zero vector of $\mathbb{R}^m$ and $\Omega_{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:

$$\forall i \in [1, n], \forall (j_1, j_2) \in [1, N_i]^2, j_1 \neq j_2$$
$$ \Rightarrow \mathrm{Cov}(\xi_{ij_1}, \xi_{ij_2}) = 0. \qquad (13) $$

Furthermore, for all individuals, the noise is supposed to be homoscedastic, that is to say: $\Omega_{ij} = \Omega_i \forall j \in [1, N_i]$. Note that the only difference between equations (5) and (12a) is that all terms are now vectors rather than scalars and that equation (12c) is the analogous of equation (7) for vectors instead of scalars.

Let us define for all individuals $\mathbf{Y}_i$, the matrix of the output data of the ith individual with $\mathbf{Y}_i := (\mathbf{y}_{i1}, \dots, \mathbf{y}_{iN_i}) \in \mathcal{M}_{m,N_i}(\mathbb{R})$. The rows of $\mathbf{Y}_i$ correspond to $\mathbf{y}_i^k = (y_{ij}^k)_{j \in [1, N_i]}$, that is to say the $N_i$ observations of the k-th measure of the ith individual. Similarly, it is possible to define the matrix of the errors as $\xi_i := (\xi_{i1}, \dots, \xi_{iN_i}) \in \mathcal{M}_{m,N_i}(\mathbb{R})$. The rows of the $\xi_i$ matrix correspond to $\xi_i^k = (\xi_{ij}^k)_{j \in [1, N_i]}$ and can be understood as the $N_i$ measurement errors of the k-th measure for the ith individual. Let us define $\tilde{\mathbf{Y}}_i$ the flattened array by row of the matrix of the output data $\mathbf{Y}_i$ and the same for the matrix of the discrepancy term $\tilde{\xi}_i$. Then, $\tilde{\mathbf{Y}}_i$ and $\tilde{\xi}_i$ belong to $\mathbb{R}^{mN_i}$. This allows to transform equation (12) into

$$ \left\{ \begin{array}{l} \forall i \in [1, n], \quad \tilde{\mathbf{Y}}_i = \mathcal{F}(\mathbf{t}_i, \boldsymbol{\theta}_i) + \tilde{\boldsymbol{\xi}}_i, \qquad (14a) \\ \forall i \in [1, n], \quad \boldsymbol{\theta}_i \sim_{\text{i.i.d.}} f_{\boldsymbol{\Theta}|\mathbf{\Pi}}, \qquad (14b) \\ \forall i \in [1, n], \quad \tilde{\boldsymbol{\xi}}_i \sim_{\text{i.i.d.}} \mathcal{N}(\mathbf{0}_{\mathbb{R}^{mN_i}}, \tilde{\boldsymbol{\Omega}}_i). \qquad (14c) \end{array} \right. $$

with $\tilde{\boldsymbol{\Omega}}_i = \boldsymbol{\Omega}_i \otimes \mathbf{I}_{N_i}$ where $\otimes$ denotes the Kronecker product, $\mathbf{I}_{N_i}$ the identity matrix of $\mathcal{M}_{N_i}(\mathbb{R})$ and $\mathbf{0}_{\mathbb{R}^{mN_i}}$ the zero vector of $\mathbb{R}^{mN_i}$. Hence, the set of the model parameters to be calibrated, labeled $\boldsymbol{\Psi}$ is defined as:

$$ \boldsymbol{\Psi} := (\boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\Omega}_1, \dots, \boldsymbol{\Omega}_m). $$

Finally, the likelihood of the mixed-effects becomes $f(\tilde{\mathbf{Y}}_1, \dots, \tilde{\mathbf{Y}}_n | \boldsymbol{\Psi})$. Equation (10) remains the same except for $\mathbf{y}_i$ that becomes $\tilde{\mathbf{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 := \int_{\mathbb{R}^d} \exp(-h(\mathbf{x})) d\mathbf{x}, $$

where $h(\cdot)$ is a function which must satisfy the following constraints:

1. $h(\cdot)$ admits a global minimum $\mathbf{x}_0$ that belongs to the integration interval,

2. $h(\cdot)$ is a twice-differentiable function,

3. its Hessian matrix $\mathcal{H}(h)$ computed at $\mathbf{x} = \mathbf{x}_0$ is a symmetric definite positive matrix.

The main idea is to assume that only regions close to $\mathbf{x}_0$ significantly contribute to the integral, leading to the following Laplace approximation of A:

$$ A \approx \exp(-h(\mathbf{x}_0)) \frac{(2\pi)^{\frac{d}{2}}}{\sqrt{|\mathcal{H}(h)(\mathbf{x}_0)|}}. \qquad (15) $$

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

$$ \begin{aligned} \mathcal{L}_i(\boldsymbol{\Psi}) &= \int_{\mathbb{R}^d} f(\tilde{\mathbf{Y}}_i | \boldsymbol{\theta}_i, \mathbf{\Pi}, \tilde{\boldsymbol{\Omega}}_i) f_{\boldsymbol{\Theta}|\mathbf{\Pi}}(\boldsymbol{\theta}_i | \mathbf{\Pi}) d\boldsymbol{\theta}_i \\ &= \int_{\mathbb{R}^d} \frac{1}{\sqrt{|\tilde{\boldsymbol{\Omega}}_i| |\boldsymbol{\Sigma}| (2\pi)^{d+mN_i}}} \qquad (16) \\ &\quad \times \exp \left( -\frac{g_i(\boldsymbol{\mu}, \mathbf{\Delta}, \mathbf{\Lambda}_i, \tilde{\mathbf{Y}}_i, \mathbf{b}_i)}{2} \right) d\mathbf{b}_i . \end{aligned} $$

Here, $\mathbf{\Delta}$ is the transpose of the result of the Cholesky decomposition of $\boldsymbol{\Sigma}^{-1}$ (so $\boldsymbol{\Sigma}^{-1} = \mathbf{\Delta}^T \mathbf{\Delta}$), $\mathbf{\Lambda}_i$ is the transpose of the Cholesky decomposition of $\tilde{\boldsymbol{\Omega}}_i^{-1}$ (so $\tilde{\boldsymbol{\Omega}}_i^{-1} = \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i$) and the function $g_i(\cdot)$ is defined by:

$$ g_i(\boldsymbol{\mu}, \mathbf{\Delta}, \mathbf{\Lambda}_i, \tilde{\mathbf{Y}}_i, \mathbf{b}_i) := \| \mathbf{\Lambda}_i(\tilde{\mathbf{Y}}_i - \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \mathbf{b}_i)) \|^2 \\ + \| \mathbf{\Delta} \mathbf{b}_i \|^2. \qquad (17) $$

The Laplace method is applied in two steps:

1. Search for the individual parameters (or rather the deviations), $\hat{\mathbf{b}}_i$, minimizing $g_i(\cdot)$ 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 $g_i(\cdot)$ at $\hat{\mathbf{b}}_i$ [32]:

$$ \mathcal{H}(g_i)(\hat{\mathbf{b}}_i) = \frac{\partial^2 \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i \partial \mathbf{b}_i^T} \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i (\tilde{\mathbf{Y}}_i - \mathcal{F}(\boldsymbol{\mu}, \mathbf{t}_i, \hat{\mathbf{b}}_i)) \\ + \frac{\partial \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i} \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i \frac{\partial \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i} \\ + \mathbf{\Delta}^T \mathbf{\Delta}. \qquad (18) $$

In practice, $\frac{\partial^2 \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i \partial \mathbf{b}_i^T} \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i (\tilde{\mathbf{Y}}_i - \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i))$ can be neglected if the model $\mathcal{F}(\cdot, \cdot)$ is close enough to the experiment $\tilde{\mathbf{Y}}_i$ [43]. The term $\frac{\partial \mathcal{F}(\boldsymbol{\mu}, \mathbf{t}_i, \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i}$ is evaluated using a finite difference scheme. As a result, the negative log-likelihood is finally expressed as

See equation (19) next page
$$ \Rightarrow \mathrm{Cov}(\xi_{ij_1}, \xi_{ij_2}) = 0. \qquad (13) $$

Furthermore, for all individuals, the noise is supposed to be homoscedastic, that is to say: $\Omega_{ij} = \Omega_i \forall j \in [1, N_i]$. Note that the only difference between equations (5) and (12a) is that all terms are now vectors rather than scalars and that equation (12c) is the analogous of equation (7) for vectors instead of scalars.

Let us define for all individuals $\mathbf{Y}_i$, the matrix of the output data of the ith individual with $\mathbf{Y}_i := (\mathbf{y}_{i1}, \dots, \mathbf{y}_{iN_i}) \in \mathcal{M}_{m,N_i}(\mathbb{R})$. The rows of $\mathbf{Y}_i$ correspond to $\mathbf{y}_i^k = (y_{ij}^k)_{j \in [1, N_i]}$, that is to say the $N_i$ observations of the k-th measure of the ith individual. Similarly, it is possible to define the matrix of the errors as $\xi_i := (\xi_{i1}, \dots, \xi_{iN_i}) \in \mathcal{M}_{m,N_i}(\mathbb{R})$. The rows of the $\xi_i$ matrix correspond to $\xi_i^k = (\xi_{ij}^k)_{j \in [1, N_i]}$ and can be understood as the $N_i$ measurement errors of the k-th measure for the ith individual. Let us define $\tilde{\mathbf{Y}}_i$ the flattened array by row of the matrix of the output data $\mathbf{Y}_i$ and the same for the matrix of the discrepancy term $\tilde{\xi}_i$. Then, $\tilde{\mathbf{Y}}_i$ and $\tilde{\xi}_i$ belong to $\mathbb{R}^{mN_i}$. This allows to transform equation (12) into

$$ \left\{ \begin{array}{l} \forall i \in [1, n], \quad \tilde{\mathbf{Y}}_i = \mathcal{F}(\mathbf{t}_i, \boldsymbol{\theta}_i) + \tilde{\boldsymbol{\xi}}_i, \qquad (14a) \\ \forall i \in [1, n], \quad \boldsymbol{\theta}_i \sim_{\text{i.i.d.}} f_{\boldsymbol{\Theta}|\mathbf{\Pi}}, \qquad (14b) \\ \forall i \in [1, n], \quad \tilde{\boldsymbol{\xi}}_i \sim_{\text{i.i.d.}} \mathcal{N}(\mathbf{0}_{\mathbb{R}^{mN_i}}, \tilde{\boldsymbol{\Omega}}_i). \qquad (14c) \end{array} \right. $$

with $\tilde{\boldsymbol{\Omega}}_i = \boldsymbol{\Omega}_i \otimes \mathbf{I}_{N_i}$ where $\otimes$ denotes the Kronecker product, $\mathbf{I}_{N_i}$ the identity matrix of $\mathcal{M}_{N_i}(\mathbb{R})$ and $\mathbf{0}_{\mathbb{R}^{mN_i}}$ the zero vector of $\mathbb{R}^{mN_i}$. Hence, the set of the model parameters to be calibrated, labeled $\mathbf{\Psi}$ is defined as:

$$ \mathbf{\Psi} := (\boldsymbol{\mu}, \boldsymbol{\Sigma}, \boldsymbol{\Omega}_1, \dots, \boldsymbol{\Omega}_m). $$

Finally, the likelihood of the mixed-effects becomes $f(\tilde{\mathbf{Y}}_1, \dots, \tilde{\mathbf{Y}}_n|\mathbf{\Psi})$. Equation (10) remains the same except for $\mathbf{y}_i$ that becomes $\tilde{\mathbf{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 := \int_{\mathbb{R}^d} \exp(-h(\mathbf{x})) d\mathbf{x}, $$

where $h(\cdot)$ is a function which must satisfy the following constraints:

1. $h(\cdot)$ admits a global minimum $\mathbf{x}_0$ that belongs to the integration interval,

2. $h(\cdot)$ is a twice-differentiable function,

3. its Hessian matrix $\mathcal{H}(h)$ computed at $\mathbf{x} = \mathbf{x}_0$ is a symmetric definite positive matrix.

The main idea is to assume that only regions close to $\mathbf{x}_0$ significantly contribute to the integral, leading to the following Laplace approximation of A:

$$ A \approx \exp(-h(\mathbf{x}_0)) \frac{(2\pi)^{\frac{d}{2}}}{\sqrt{|\mathcal{H}(h)(\mathbf{x}_0)|}}. \qquad (15) $$

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

$$ \begin{aligned} \mathcal{L}_i(\mathbf{\Psi}) &= \int_{\mathbb{R}^d} f(\tilde{\mathbf{Y}}_i | \boldsymbol{\theta}_i, \mathbf{\Pi}, \tilde{\boldsymbol{\Omega}}_i) f_{\boldsymbol{\Theta}|\mathbf{\Pi}}(\boldsymbol{\theta}_i | \mathbf{\Pi}) d\boldsymbol{\theta}_i \\ &= \int_{\mathbb{R}^d} \frac{1}{\sqrt{|\tilde{\boldsymbol{\Omega}}_i| |\boldsymbol{\Sigma}| (2\pi)^{d+mN_i}}} \qquad (16) \\ & \quad \times \exp \left( - \frac{g_i(\boldsymbol{\mu}, \mathbf{\Delta}, \mathbf{\Lambda}_i, \tilde{\mathbf{Y}}_i, \mathbf{b}_i)}{2} \right) d\mathbf{b}_i. \end{aligned} $$

Here, $\mathbf{\Delta}$ is the transpose of the result of the Cholesky decomposition of $\boldsymbol{\Sigma}^{-1}$ (so $\boldsymbol{\Sigma}^{-1} = \mathbf{\Delta}^T \mathbf{\Delta}$), $\mathbf{\Lambda}_i$ is the transpose of the Cholesky decomposition of $\tilde{\boldsymbol{\Omega}}_i^{-1}$ (so $\tilde{\boldsymbol{\Omega}}_i^{-1} = \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i$) and the function $g_i(\cdot)$ is defined by:

$$ g_i(\boldsymbol{\mu}, \mathbf{\Delta}, \mathbf{\Lambda}_i, \tilde{\mathbf{Y}}_i, \mathbf{b}_i) := \| \mathbf{\Lambda}_i(\tilde{\mathbf{Y}}_i - \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \mathbf{b}_i)) \|^2 \\ + \| \mathbf{\Delta} \mathbf{b}_i \|^2. \qquad (17) $$

The Laplace method is applied in two steps:

1. Search for the individual parameters (or rather the deviations), $\hat{\mathbf{b}}_i$, minimizing $g_i(\cdot)$ 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 $g_i(\cdot)$ at $\hat{\mathbf{b}}_i$ [32]:

$$ \mathcal{H}(g_i)(\hat{\mathbf{b}}_i) = \frac{\partial^2 \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i \partial \mathbf{b}_i^T} \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i (\tilde{\mathbf{Y}}_i - \mathcal{F}(\boldsymbol{\mu}, \mathbf{t}_i, \hat{\mathbf{b}}_i)) \\ + \frac{\partial \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i} \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i \frac{\partial \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i} \\ + \mathbf{\Delta}^T \mathbf{\Delta}. \qquad (18) $$

In practice, $\frac{\partial^2 \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i \partial \mathbf{b}_i^T} \mathbf{\Lambda}_i^T \mathbf{\Lambda}_i (\tilde{\mathbf{Y}}_i - \mathcal{F}(\mathbf{t}_i, \boldsymbol{\mu} + \hat{\mathbf{b}}_i))$ can be neglected if the model $\mathcal{F}(\cdot, \cdot)$ is close enough to the experiment $\tilde{\mathbf{Y}}_i$ [43]. The term $\frac{\partial \mathcal{F}(\boldsymbol{\mu}, \mathbf{t}_i, \hat{\mathbf{b}}_i)}{\partial \mathbf{b}_i}$ is evaluated using a finite difference scheme. As a result, the negative log-likelihood is finally expressed as

See equation (19) next page
$$ -\ln(\mathcal{L}(\Psi)) = -\sum_{i=1}^{n} \ln(\mathcal{L}_i(\Psi)) = -\sum_{i=1}^{n} \ln\left(\int_{\mathbb{R}^d} f(\tilde{Y}_i|\theta_i, \Pi, \tilde{\Omega}_i) f_{\Theta|\Pi}(\theta_i|\Pi) d\theta_i\right) \\ \approx \sum_{i=1}^{n} \left(\frac{1}{2}\ln(|\mathcal{H}(g_i)(\hat{b}_i)|) + \frac{g_i(\mu, \Delta, \Lambda_i, \tilde{Y}_i, \hat{b}_i)}{2} + \frac{1}{2}\ln(|\tilde{\Omega}_i|(2\pi)^{mN_i})\right) \\ + \frac{n}{2}\ln(|\Sigma|). \qquad (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 $\Psi$. Once $\Psi$ is chosen, the random effects $b_i$ are estimated by looking for the dominating contribution in the individual likelihood $\mathcal{L}_i$ (Eq. (16)) through the minimization of function $g_i(\cdot)$ (Eq. (17))

$$ \hat{b}_i = \arg\min_{b_i \in \Gamma_b} g_i(\mu, \Delta, \Lambda_i, \tilde{Y}_i, b_i), \qquad (20) $$

with $\Gamma_b$ the domain definition for the individual deviations. The individual deviations are then simply computed as $\theta_i = \mu + b_i$.

The determination of the individual parameters is a key step to ensure a proper likelihood maximization and should not be neglected. To explain its importance, it should be underlined that mixed-effects aim to tune the population parameters ($\Pi$) in order to find the best possible agreement between the assumed model parameters distribution ($f_{\Theta|\Pi}$) 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_{\Theta|\Pi}(\theta_i|\Pi)$ in equation (10) or $||\Delta b_i||^2$ in function $g_i(\cdot)$ 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_{\Theta|\Pi}$, 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,

$$ \hat{\Psi} = \arg\min_{\Psi \in \Gamma_{\Psi}} -\ln(\mathcal{L}(\Psi)), \qquad (21) $$

with $\Gamma_{\Psi}$ the research space for $\Psi$. 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(\mathcal{L}(\Psi))$ expressed with equation (19) with respect to $\Psi$ is difficult to compute analytically. Indeed, the differentiation with respect to the population parameters of function $g_i(\cdot)$ requires to differentiate the individual deviations $b_i$ defined in equation (20) with respect to the population parameters (i.e., $\mu$ and $\Sigma$ via $\Delta$). 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, $\hat{\Psi}$, despite the possible presence of local minima. The algorithm chosen here is the Covariance Matrix Adaptation 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 $\Psi$ labeled $p$ complies with $p \ge 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
$$
\begin{equation}
\begin{aligned}
-\ln(\mathcal{L}(\Psi)) = & -\sum_{i=1}^{n} \ln(\mathcal{L}_i(\Psi)) = -\sum_{i=1}^{n} \ln\left(\int_{\mathbb{R}^d} f(\tilde{Y}_i|\theta_i, \Pi, \tilde{\Omega}_i) f_{\Theta|\Pi}(\theta_i|\Pi) d\theta_i\right) \\
& \approx \sum_{i=1}^{n} \left(\frac{1}{2}\ln(|\mathcal{H}(g_i)(\hat{b}_i)|) + \frac{g_i(\mu, \Delta, \Lambda_i, \tilde{Y}_i, \hat{b}_i)}{2} + \frac{1}{2}\ln(|\tilde{\Omega}_i|(2\pi)^{mN_i})\right) \\
& + \frac{n}{2}\ln(|\Sigma|).
\end{aligned}
\tag{19}
\end{equation}
$$

## 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 $\Psi$. Once $\Psi$ is chosen, the random effects $b_i$ are estimated by looking for the dominating contribution in the individual likelihood $\mathcal{L}_i$ (Eq. (16)) through the minimization of function $g_i(\cdot)$ (Eq. (17))

$$ \hat{\mathbf{b}}_i = \arg\min_{\mathbf{b}_i \in \Gamma_b} g_i(\mu, \Delta, \Lambda_i, \tilde{Y}_i, \mathbf{b}_i), \quad (20) $$

with $\Gamma_b$ the domain definition for the individual deviations. The individual deviations are then simply computed as $\theta_i = \mu + \mathbf{b}_i$.

The determination of the individual parameters is a key step to ensure a proper likelihood maximization and should not be neglected. To explain its importance, it should be underlined that mixed-effects aim to tune the population parameters ($\Pi$) in order to find the best possible agreement between the assumed model parameters distribution ($f_{\Theta|\Pi}$) 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_{\Theta|\Pi}(\theta_i|\Pi)$ in equation (10) or $\|\Delta\mathbf{b}_i\|^2$ in function $g_i(\cdot)$ 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_{\Theta|\Pi}$, 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,

$$ \hat{\Psi} = \arg\min_{\Psi \in \Gamma_{\Psi}} -\ln(\mathcal{L}(\Psi)), \quad (21) $$

with $\Gamma_{\Psi}$ the research space for $\Psi$. 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(\mathcal{L}(\Psi))$ expressed with equation (19) with respect to $\Psi$ is difficult to compute analytically. Indeed, the differentiation with respect to the population parameters of function $g_i(\cdot)$ requires to differentiate the individual deviations $b_i$ defined in equation (20) with respect to the population parameters (i.e., $\mu$ and $\Sigma$ via $\Delta$). 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, $\hat{\Psi}$, despite the possible presence of local minima. The algorithm chosen here is the Covariance Matrix Adaptation 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 $\Psi$ labeled $p$ complies with $p \ge 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
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

into calibration subproblems that allow to identify dif-
ferent 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 meth-
ods 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 decomposi-
tion 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
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

into calibration subproblems that allow to identify dif-
ferent 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 meth-
ods 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 decomposi-
tion 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_{\Theta|\Pi}$, it gives :

$$f_{\Theta|\Pi}(\theta_i|\Pi) = c(F_1(\theta_{il}|\Pi_1), \dots, F_d(\theta_{id}|\Pi_d)|\Pi_{\text{cop}}) \times \prod_{l=1}^{d} f_l(\theta_{il}|\Pi_l), \quad (22)$$

with $\Pi_{\text{cop}}$ the copula parameters, $\Pi_l$ the parameters of the $l^{th}$ marginal distribution $f_{\Theta|\Pi}$ (so the marginal marginal distribution of $\theta_{il}$, $(i,l) \in [1, n] \times [1, d]$) with probability density function $f_l(\cdot)$ and cumulative density function $F_l(\cdot)$ and $c(\cdot)$, 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 fulfilled 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_{\Theta|\Pi}$. The sequential calibration process consists in defining a sequence of subsets of indices $S_q \subset [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 $\Pi_q$ to be calibrated at step $q$, leading to the definition of a specific likelihood function $L^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 $R_q$ is defined. If some of the parameters of the $q^{th}$ step are re-estimated, they are limited to their trust region $R_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 $[\alpha\Pi_{l,1}, \beta\Pi_{l,1}] \times \dots \times [\alpha\Pi_{l,q-1}, \beta\Pi_{l,q-1}]$ for the parameters $l \in [1, s_{q-1}]$, with $\alpha$ and $\beta$ indicating the size of the interval, $s_{q-1}$ the number of common parameters from step $q-1$ to $q$ and $[\Pi_{l,q-1}]_{l \in [1, s_{q-1}]}$ the corresponding coordinates. Note that the size of the intervals $\alpha$ and $\beta$ can also depend on the step $q$ and be different for each coordinate. Most of the time, $\alpha$ will be strictly smaller than $\beta$, meaning that the $[\Pi_{l,q-1}]_{l \in [1, s_{q-1}]}$ 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
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_{\Theta|\Pi}$, it gives :

$$f_{\Theta|\Pi}(\theta_i|\Pi) = c(F_1(\theta_{il}|\Pi_1), \dots, F_d(\theta_{id}|\Pi_d)|\Pi_{\text{cop}}) \times \prod_{l=1}^{d} f_l(\theta_{il}|\Pi_l), \quad (22)$$

with $\Pi_{\text{cop}}$ the copula parameters, $\Pi_l$ the parameters of the $l^{th}$ marginal distribution $f_{\Theta|\Pi}$ (so the marginal marginal distribution of $\theta_{il}$, $(i,l) \in [1, n] \times [1, d]$) with probability density function $f_l(\cdot)$ and cumulative density function $F_l(\cdot)$ and $c(\cdot)$, 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 fulfilled 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_{\Theta|\Pi}$. The sequential calibration process consists in defining a sequence of subsets of indices $S_q \subset [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 $\Pi_q$ to be calibrated at step $q$, leading to the definition of a specific likelihood function $L^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 $R_q$ is defined. If some of the parameters of the $q^{th}$ step are re-estimated, they are limited to their trust region $R_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 $[\alpha\Pi_{l,1}, \beta\Pi_{l,1}] \times \dots \times [\alpha\Pi_{l,q-1}, \beta\Pi_{l,q-1}]$ for the parameters $l \in [1, s_{q-1}]$, with $\alpha$ and $\beta$ indicating the size of the interval, $s_{q-1}$ the number of common parameters from step $q-1$ to $q$ and $[\Pi_{l,q-1}]_{l \in [1, s_{q-1}]}$ the corresponding coordinates. Note that the size of the intervals $\alpha$ and $\beta$ can also depend on the step $q$ and be different for each coordinate. Most of the time, $\alpha$ will be strictly smaller than $\beta$, meaning that the $[\Pi_{l,q-1}]_{l \in [1, s_{q-1}]}$ 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
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Different stages of the calibration process. Once the sequence of calibration problems to resolve is defined, the calibration process starts. At each iteration, the search space is defined depending on whether parameters have already been calibrated to preserve consistency between the different steps. Calibration stages are highlighted in red and those concerning the determination of the search space are highlighted in blue. At each iteration, $u_q$ refers to the integer labeling the subgroup of parameters that will be re-estimated in the current iteration. $N_G$ refers to the number of calibration subproblems. The process stops when all calibration subproblems are solved.</figcaption></figure>

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 $\Psi$ 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 from 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 $\alpha = 0.8$ and $\beta = 1.2$ (the boundaries correspond to $\pm 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.

# 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
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Different stages of the calibration process. Once the sequence of calibration problems to resolve is defined, the calibration process starts. At each iteration, the search space is defined depending on whether parameters have already been calibrated to preserve consistency between the different steps. Calibration stages are highlighted in red and those concerning the determination of the search space are highlighted in blue. At each iteration, $u_q$ refers to the integer labeling the subgroup of parameters that will be re-estimated in the current iteration. $N_G$ refers to the number of calibration subproblems. The process stops when all calibration subproblems are solved.</figcaption></figure>

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 $\Psi$ 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 from 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 $\alpha = 0.8$ and $\beta = 1.2$ (the boundaries correspond to $\pm 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.

# 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
<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

synthetic data first to study the consequences of such a methodology in Section 4.5.

## 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 $\sigma := (\sigma_{11}, \sigma_{22}, \sigma_{12})$ the Cauchy stress (in MPa) and $\varepsilon := (\varepsilon_{11}, \varepsilon_{22}, 2\varepsilon_{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:

$$ \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \end{pmatrix} = \underline{\underline{\mathbf{S}}} \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix} = \begin{pmatrix} S_{11} & S_{21} & 0 \\ S_{12} & S_{22} & 0 \\ 0 & 0 & S_{66} \end{pmatrix} \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix}, \quad (23) $$

with $\underline{\underline{\mathbf{S}}}$ the compliance tensor. $\underline{\underline{\mathbf{S}}}$ is a symmetric tensor so $S_{12} = S_{21}$. For practical interpretations, the coefficients of the compliance tensor can be related to the elastic moduli, that is to say the longitudinal modulus $E_{11}$, the transverse modulus $E_{22}$, the shear modulus $G_{12}$ and the Poisson's ratio $\nu_{12}$ [25].

The model presented above is elastic linear and can be used either in tension and compression. Yet, it has been observed that for some materials, in tension, 0° degree laminate tends to stiffen, while in compression, 0° degree laminate tend to soften [56]. This is the consequence of the fiber behavior. Indeed, when produced, the fibers are not always properly aligned and when loaded in tension, fibers get aligned increasing the stiffness, while in compression, the initial defaults are increased and tend to soften the material. Thus, the apparent modulus increases (non linearly) throughout the tension experiments and conversely on the compression experiments. Thus, the longitudinal compliance $S_{11}$ is changed and expresses as follows:

$$ S_{11}^{L} = S_{1}^{L} + \frac{(S_{11}^{0} - S_{1}^{L})\varepsilon_{0}^{L}}{(S_{11}^{0} - S_{1}^{L})\sigma_{11} + \varepsilon_{0}^{L}}, \quad (24) $$

with $S_{11}^0$ the elastic linear compliance (module at the origin), $S_1^L$ the asymptotic compliance and $\varepsilon_0^L$ 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 $S_{12}^L$ should also be updated and refers to $S_{12}^L := -\nu_{12}S_{11}^L$. 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 $S_1^T < S_{11}^0$ and $S_1^C > S_{11}^0$. With this parameterization, the model is not identifiable. In fact, it is over-parameterized, i.e. there are several combinations of both $S_1^L$ and $\varepsilon_0^L$ that give the same output as shown by Germain [57]. Thus, $\varepsilon_0^L$ fixed to a usual nominal value, here 0.005 [57]. Finally, there are 6 parameters to be calibrated $S_1^T, S_1^C, S_{11}^0, S_{22}, \nu_{12}$ and $S_{66}$.

In theory, the compliance matrix $\underline{\underline{\mathbf{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 $S_{66}$ 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-definiteness 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
<figure><img src="image_6.png" /><figcaption>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.</figcaption></figure>

synthetic data first to study the consequences of such a methodology in Section 4.5.

## 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 $\sigma := (\sigma_{11}, \sigma_{22}, \sigma_{12})$ the Cauchy stress (in MPa) and $\varepsilon := (\varepsilon_{11}, \varepsilon_{22}, 2\varepsilon_{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:

$$ \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \end{pmatrix} = \underline{\underline{\mathbf{S}}} \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix} = \begin{pmatrix} S_{11} & S_{21} & 0 \\ S_{12} & S_{22} & 0 \\ 0 & 0 & S_{66} \end{pmatrix} \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{pmatrix}, \quad (23) $$

with $\underline{\underline{\mathbf{S}}}$ the compliance tensor. $\underline{\underline{\mathbf{S}}}$ is a symmetric tensor so $S_{12} = S_{21}$. For practical interpretations, the coefficients of the compliance tensor can be related to the elastic moduli, that is to say the longitudinal modulus $E_{11}$, the transverse modulus $E_{22}$, the shear modulus $G_{12}$ and the Poisson's ratio $\nu_{12}$ [25].

The model presented above is elastic linear and can be used either in tension and compression. Yet, it has been observed that for some materials, in tension, 0° degree laminate tends to stiffen, while in compression, 0° degree laminate tend to soften [56]. This is the consequence of the fiber behavior. Indeed, when produced, the fibers are not always properly aligned and when loaded in tension, fibers get aligned increasing the stiffness, while in compression, the initial defaults are increased and tend to soften the material. Thus, the apparent modulus increases (non linearly) throughout the tension experiments and conversely on the compression experiments. Thus, the longitudinal compliance $S_{11}$ is changed and expresses as follows:

$$ S_{11}^{L} = S_{1}^{L} + \frac{(S_{11}^{0} - S_{1}^{L})\varepsilon_{0}^{L}}{(S_{11}^{0} - S_{1}^{L})\sigma_{11} + \varepsilon_{0}^{L}}, \quad (24) $$

with $S_{11}^0$ the elastic linear compliance (module at the origin), $S_1^L$ the asymptotic compliance and $\varepsilon_0^L$ 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 $S_{12}^L$ should also be updated and refers to $S_{12}^L := -\nu_{12}S_{11}^L$. 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 $S_1^T < S_{11}^0$ and $S_1^C > S_{11}^0$. With this parameterization, the model is not identifiable. In fact, it is over-parameterized, i.e. there are several combinations of both $S_1^L$ and $\varepsilon_0^L$ that give the same output as shown by Germain [57]. Thus, $\varepsilon_0^L$ fixed to a usual nominal value, here 0.005 [57]. Finally, there are 6 parameters to be calibrated $S_1^T, S_1^C, S_{11}^0, S_{22}, \nu_{12}$ and $S_{66}$.

In theory, the compliance matrix $\underline{\underline{\mathbf{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 $S_{66}$ 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-definiteness 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
$$ \Sigma = \begin{pmatrix} S_{11}^0 & S_1^T & S_1^C & \nu_{12} & S_{22} & S_{66} \\ S_1^0 & \mathbb{V}(S_{11}^0) & * & * & * & * \\ S_1^T & \mathrm{Cov}(S_{11}^0, S_1^T) & \mathbb{V}(S_1^T) & * & * & * \\ S_1^C & \mathrm{Cov}(S_{11}^0, S_1^C) & \mathrm{Cov}(S_1^T, S_1^C) & \mathbb{V}(S_1^C) & * & * \\ \nu_{12} & \mathrm{Cov}(S_{11}^0, \nu_{12}) & \mathrm{Cov}(S_1^T, \nu_{12}) & \mathrm{Cov}(S_1^C, \nu_{12}) & \mathbb{V}(\nu_{12}) & * \\ S_{22} & \mathrm{Cov}(S_{11}^0, S_{22}) & \mathrm{Cov}(S_1^T, S_{22}) & \mathrm{Cov}(S_1^C, S_{22}) & \mathrm{Cov}(S_{22}, \nu_{12}) & \mathbb{V}(S_{22}) \\ S_{66} & \mathrm{Cov}(S_{11}^0, S_{66}) & \mathrm{Cov}(S_1^T, S_{66}) & \mathrm{Cov}(S_1^C, S_{66}) & \mathrm{Cov}(\nu_{12}, S_{66}) & \mathrm{Cov}(S_{22}, S_{66}) & \mathbb{V}(S_{66}) \end{pmatrix} \quad (25) $$

can be summed up into a mean vector $\mu = (\mu_{S_{11}^0}, \mu_{S_1^T}, \mu_{S_1^C}, \mu_{\nu_{12}}, \mu_{S_{22}}, \mu_{S_{66}})^\top$ and covariance matrix

See equation (25) above

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, $S_{11}^0, S_1^T, S_1^C$ and $\nu_{12}$ characterize the fiber behavior and are almost exclusively monitor the 0° laminate. Similarly, the linear part of the longitudinal response of a tension/compression test on a 90° laminate is mainly controlled by $S_{22}$, the same applies for the linear part of the difference between the longitudinal and transverse responses for the ±45° laminate, mostly driven by $S_{66}$ [25,58,59]. Remember that, with 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 $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and $S_{22}$ cannot be estimated with the usual elementary experiments, and the same applies to those between $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and $S_{66}$ or between $S_{22}$ and $S_{66}$. In fact, from these elementary experiments, only the joint distribution of $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and the marginals of $S_{66}$ and $S_{22}$ can be estimated, which gives three distributions to calibrate (in addition to the residual models). This makes 18 population parameters (14 for the joint distribution of $S_{11}^0, S_1^T, S_1^C, \nu_{12}$, and 2 for both marginals of $S_{22}$ and $S_{66}$) to estimate.

From this assessment, two solutions arise. On the one hand, it is possible to calibrate all the experiments together by making sure that for each individual, only the relevant individual parameters are calibrated. This can be done through equation (8) by choosing relevant design matrices $A_i$ and $B_i$. For instance, with the mean vector defined earlier, one possible choice would be $A_i = B_i$ equal to the identity matrix of dimension 6 to which the two last rows would be deleted (and $A_i$ would have thus 4 rows with 6 columns) for the UD0° experiments, $A_i = B_i = (0, 0, 0, 0, 1, 0)^\top$ for the 90° experiments and $A_i = B_i = (0, 0, 0, 0, 0, 1)^\top$ 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 $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and the two other marginals independently from the others (which thus do not require to use design matrix). If the estimation of these distributions is carried out separately, this can be considered as a first application of the sequential strategy in an extreme case, given that between the three cases, there are not shared parameters, and the calibration of the three distributions can thus be achieved independently. The calibration of the marginals of $S_{66}$ and $S_{22}$ are carried out from experiments with univariate observations that exhibit a linear behavior. Thus, these cases only refer to univariate linear 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 $S_{11}^0, S_1^T, S_1^C$ and $\nu_{12}$.

One first remark that can be made is that $S_1^C$ only expresses when compressive loads are applied to the sample, and $S_1^T$ 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 $S_{11}^0, S_1^T, S_1^C$ and $\nu_{12}$ is parameterized by 13 and not 14 population parameters. $S_{11}^0$ and $\nu_{12}$ are activated both in tension and compression. Note that $\nu_{12}$ only expresses on the transverse responses, either in tension and compression. Here, there are shared parameters for two types of experiments, $S_{11}^0$ and $\nu_{12}$. Design matrices apply here. For instance, for the tension experiments, matrices $A_i$ and $B_i$ would correspond to the identity matrix of dimension 4 to which is deleted the third row (so both have 3 rows and 4 columns) and for the compression experiments, they could be taken as the identity matrix of dimension 4 to which is deleted the third row (so both have 3 rows and 4 columns). The interest of this test case is that it also allows to illustrate the sequential strategy. In particular, it enables to compare the accuracy of the calibration from the sequential strategy to the solution that implements design matrices. For 0° UD laminate, let consider a tension load $\sigma = (\sigma_{11}, 0, 0)$ with $\sigma_{11} > 0$. Then, following equation (23), the model output is determined by $S_{11}^0, S_1^T$ and $\nu_{12}$. It is thus possible to calibrate the parameters highlighted in red in equation (26):

See equation (26) next page
$$ \Sigma = \begin{pmatrix} S_{11}^0 & S_1^T & S_1^C & \nu_{12} & S_{22} & S_{66} \\ S_1^0 & \mathbb{V}(S_{11}^0) & * & * & * & * \\ S_1^T & \mathrm{Cov}(S_{11}^0, S_1^T) & \mathbb{V}(S_1^T) & * & * & * \\ S_1^C & \mathrm{Cov}(S_{11}^0, S_1^C) & \mathrm{Cov}(S_1^T, S_1^C) & \mathbb{V}(S_1^C) & * & * \\ \nu_{12} & \mathrm{Cov}(S_{11}^0, \nu_{12}) & \mathrm{Cov}(S_1^T, \nu_{12}) & \mathrm{Cov}(S_1^C, \nu_{12}) & \mathbb{V}(\nu_{12}) & * \\ S_{22} & \mathrm{Cov}(S_{11}^0, S_{22}) & \mathrm{Cov}(S_1^T, S_{22}) & \mathrm{Cov}(S_1^C, S_{22}) & \mathrm{Cov}(S_{22}, \nu_{12}) & \mathbb{V}(S_{22}) & * \\ S_{66} & \mathrm{Cov}(S_{11}^0, S_{66}) & \mathrm{Cov}(S_1^T, S_{66}) & \mathrm{Cov}(S_1^C, S_{66}) & \mathrm{Cov}(\nu_{12}, S_{66}) & \mathrm{Cov}(S_{22}, S_{66}) & \mathbb{V}(S_{66}) \end{pmatrix} \quad (25) $$

can be summed up into a mean vector $\mu = (\mu_{S_{11}^0}, \mu_{S_1^T}, \mu_{S_1^C}, \mu_{\nu_{12}}, \mu_{S_{22}}, \mu_{S_{66}})^\top$ and covariance matrix

See equation (25) above

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, $S_{11}^0, S_1^T, S_1^C$ and $\nu_{12}$ characterize the fiber behavior and are almost exclusively monitor the 0° laminate. Similarly, the linear part of the longitudinal response of a tension/compression test on a 90° laminate is mainly controlled by $S_{22}$, the same applies for the linear part of the difference between the longitudinal and transverse responses for the ±45° laminate, mostly driven by $S_{66}$ [25,58,59]. Remember that, with 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 $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and $S_{22}$ cannot be estimated with the usual elementary experiments, and the same applies to those between $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and $S_{66}$ or between $S_{22}$ and $S_{66}$. In fact, from these elementary experiments, only the joint distribution of $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and the marginals of $S_{66}$ and $S_{22}$ can be estimated, which gives three distributions to calibrate (in addition to the residual models). This makes 18 population parameters (14 for the joint distribution of $S_{11}^0, S_1^T, S_1^C, \nu_{12}$, and 2 for both marginals of $S_{22}$ and $S_{66}$) 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 $\mathbf{A}_i$ and $\mathbf{B}_i$. For instance, with the mean vector defined earlier, one possible choice would be $\mathbf{A}_i = \mathbf{B}_i$ equal to the identity matrix of dimension 6 to which the two last rows would be deleted (and $\mathbf{A}_i$ would have thus 4 rows with 6 columns) for the UD0° experiments, $\mathbf{A}_i = \mathbf{B}_i = (0, 0, 0, 0, 1)^\top$ for the 90° experiments and $\mathbf{A}_i = \mathbf{B}_i = (0, 0, 0, 0, 0, 1)^\top$ 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 $S_{11}^0, S_1^T, S_1^C, \nu_{12}$ and the two other marginals independently from the others (which thus do not require to use design matrix). If the estimation of these distributions is carried out separately, this can be considered as a first application of the sequential strategy in an extreme case, given that between the three cases, there are not shared parameters, and the calibration of the three distributions can thus be achieved independently. The calibration of the marginals of $S_{66}$ and $S_{22}$ are carried out from experiments with univariate observations that exhibit a linear behavior. Thus, these cases only refer to univariate linear 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 $S_{11}^0, S_1^T, S_1^C$ and $\nu_{12}$.

One first remark that can be made is that $S_1^C$ only expresses when compressive loads are applied to the sample, and $S_1^T$ 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 $S_{11}^0, S_1^T, S_1^C$ and $\nu_{12}$ is parameterized by 13 and not 14 population parameters. $S_{11}^0$ and $\nu_{12}$ are activated both in tension and compression. Note that $\nu_{12}$ only expresses on the transverse responses, either in tension and compression. Here, there are shared parameters for two types of experiments, $S_{11}^0$ and $\nu_{12}$. Design matrices apply here. For instance, for the tension experiments, matrices $\mathbf{A}_i$ and $\mathbf{B}_i$ would correspond to the identity matrix of dimension 4 to which is deleted the third row (so both have 3 rows and 4 columns) and for the compression experiments, they could be taken as the identity matrix of dimension 4 to which is deleted the third row (so both have 3 rows and 4 columns). The interest of this test case is that it also allows to illustrate the sequential strategy. In particular, it enables to compare the accuracy of the calibration from the sequential strategy to the solution that implements design matrices. For 0° UD laminate, let consider a tension load $\boldsymbol{\sigma} = (\sigma_{11}, 0, 0)$ with $\sigma_{11} > 0$. Then, following equation (23), the model output is determined by $S_{11}^0, S_1^T$ and $\nu_{12}$. It is thus possible to calibrate the parameters highlighted in red in equation (26):

See equation (26) next page
$$ \boldsymbol{\mu} = \begin{pmatrix} \mu_{S_{11}^0} \\ \mu_{S_1^C} \\ \mu_{S_1^T} \\ \mu_{\nu_{12}} \end{pmatrix} \quad \boldsymbol{\Sigma} = \begin{pmatrix} S_{11}^0 & S_1^C & S_1^T & \nu_{12} \\ S_1^0 & \mathbb{V}(S_{11}^0) & * & * \\ S_1^C & \mathbb{V}(S_1^C) & * & * \\ S_1^T & \mathbb{V}(S_1^T) & * & * \\ \nu_{12} & \mathbb{V}(\nu_{12}, S_1^C) & \mathbb{V}(\nu_{12}, S_1^T) & \mathbb{V}(\nu_{12}) \end{pmatrix}. \quad (26) $$

$$ \boldsymbol{\mu} = \begin{pmatrix} \mu_{S_{11}^0} \\ \mu_{S_1^C} \\ \mu_{S_1^T} \\ \mu_{\nu_{12}} \end{pmatrix} \quad \boldsymbol{\Sigma} = \begin{pmatrix} S_{11}^0 & S_1^C & S_1^T & \nu_{12} \\ S_1^0 & \mathbb{V}(S_{11}^0) & * & * \\ S_1^C & \mathbb{V}(S_1^C) & * & * \\ S_1^T & \mathbb{V}(S_1^T) & * & * \\ \nu_{12} & \mathbb{V}(\nu_{12}, S_1^C) & \mathbb{V}(\nu_{12}, S_1^T) & \mathbb{V}(\nu_{12}) \end{pmatrix}. \quad (27) $$

This allows to estimate a first trust region: $\mathcal{R}_1$. Let's now consider a compressive load $\sigma = (\sigma_{11}, 0, 0)$ with $\sigma_{11} < 0$. Then, following equation (23), the model output is determined by $S_{11}^0, S_1^C$ and $\nu_{12}$. The application of previous remarks shows that it is possible to calibrate the parameters highlighted in blue in equation (26):

See equation (27) above

Note that $\mu_{S_{11}^0}, \mu_{\nu_{12}}, \mathbb{V}(S_{11}^0), \mathbb{V}(S_{11}^0, \nu_{12})$ and $\mathbb{V}(\nu_{12})$ are limited to their trust region $\mathcal{R}_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 $S_{11}^0, S_1^T$ and $\nu_{12}$ by limiting the parameters of the joint distribution of $S_{11}^0, \nu_{12}$ to trust region $\mathcal{R}_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 $S_{11}^0, S_1^T, S_1^C, \nu_{12}, S_{22}$ and $S_{66}$ parameterized by 17 parameters with a dependence structure and 12 without) is divided with a sequential scheme into subproblems with at most 9 parameters with a dependence structure and 6 without. This reduction can illustrated further if another component is added, for instance the viscoelasticity behavior as in the ONERA Progressive Failure Model (OPFM) [56]. This component adds 6 new parameters and can be calibrated from both 90° and ±45° laminates (in which also express the linear compliances $S_{22}$ and $S_{66}$). 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 $O(d)$ while it is in $O(d^2)$ if all dependencies are considered, with $d$ the number of model parameters to estimate. Thus,
$$ \boldsymbol{\mu} = \begin{pmatrix} \mu_{S_{11}^0} \\ \mu_{S_1^C} \\ \mu_{S_1^T} \\ \mu_{\nu_{12}} \end{pmatrix} \quad \boldsymbol{\Sigma} = \begin{pmatrix} S_{11}^0 & S_1^C & S_1^T & \nu_{12} \\ S_1^0 & \mathbb{V}(S_{11}^0) & * & * \\ S_1^C & \mathbb{V}(S_1^C) & * & * \\ S_1^T & \mathbb{V}(S_1^T) & * & * \\ \nu_{12} & \mathbb{V}(\nu_{12}, S_1^C) & \mathbb{V}(\nu_{12}, S_1^T) & \mathbb{V}(\nu_{12}) \end{pmatrix}. \quad (26) $$

$$ \boldsymbol{\mu} = \begin{pmatrix} \mu_{S_{11}^0} \\ \mu_{S_1^C} \\ \mu_{S_1^T} \\ \mu_{\nu_{12}} \end{pmatrix} \quad \boldsymbol{\Sigma} = \begin{pmatrix} S_{11}^0 & S_1^C & S_1^T & \nu_{12} \\ S_1^0 & \mathbb{V}(S_{11}^0) & * & * \\ S_1^C & \mathbb{V}(S_1^C) & * & * \\ S_1^T & \mathbb{V}(S_1^T) & * & * \\ \nu_{12} & \mathbb{V}(\nu_{12}, S_1^C) & \mathbb{V}(\nu_{12}, S_1^T) & \mathbb{V}(\nu_{12}) \end{pmatrix}. \quad (27) $$

This allows to estimate a first trust region: $R_1$. Let's now consider a compressive load $\sigma = (\sigma_{11}, 0, 0)$ with $\sigma_{11} < 0$. Then, following equation (23), the model output is determined by $S_{11}^0, S_1^C$ and $\nu_{12}$. The application of previous remarks shows that it is possible to calibrate the parameters highlighted in blue in equation (26):

See equation (27) above

Note that $\mu_{S_{11}^0}, \mu_{\nu_{12}}, \mathbb{V}(S_{11}^0), \mathbb{V}(S_{11}^0, \nu_{12})$ and $\mathbb{V}(\nu_{12})$ are limited to their trust region $R_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 $S_{11}^0, S_1^T$ and $\nu_{12}$ by limiting the parameters of the joint distribution of $S_{11}^0, \nu_{12}$ to trust region $R_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 $S_{11}^0, S_1^T, S_1^C, \nu_{12}, S_{22}$ and $S_{66}$ parameterized by 17 parameters with a dependence structure and 12 without) is divided with a sequential scheme into subproblems with at most 9 parameters with a dependence structure and 6 without. This reduction can illustrated further if another component is added, for instance the viscoelasticity behavior as in the ONERA Progressive Failure Model (OPFM) [56]. This component adds 6 new parameters and can be calibrated from both 90° and ±45° laminates (in which also express the linear compliances $S_{22}$ and $S_{66}$). 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 $O(d)$ while it is in $O(d^2)$ if all dependencies are considered, with $d$ the number of model parameters to estimate. Thus,
<table><caption>**Table 1.** For the nonlinear elastic and viscoelastic behavioral 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.</caption><thead><tr><th>Calibration method</th><th>Maximum number of population parameters to estimate for the 2D nonlinear elasticity</th><th>Maximum number of population parameters to estimate for the 2D-viscoelasticity behavior</th></tr></thead><tbody><tr><td>Design matrices with dependence structure</td><td>17</td><td>42</td></tr><tr><td>Design matrices without dependence structure</td><td>12</td><td>16</td></tr><tr><td>Sequential approach with dependence structure</td><td>9</td><td>27</td></tr><tr><td>Sequential approach without dependence structure</td><td>6</td><td>12</td></tr></tbody></table>

because of the small number of available individuals, striking a balance between the complexity of the modeling and the number of available individuals 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 $\omega_i$s) is also to calibrate. Here again the implementation of the sequential approach allows to reduce the number of residual models (and thus of parameters) to estimate. Indeed, with the sequential approach, a focus is made on certain types of specimen (for instance compression experiments) while when with design matrices, different types of experiments are considered (tension and compression experiments). Thus, under the reasonable hypothesis of at least one residual model for each type of experiments, there are more residual models (so more parameters) to calibrate when design matrices are employed than in the different stages of the implementation of the sequential approach. Consequently, it appears that the sequential strategy is particularly suited to reduce the complexity of the minimization problem. It should be noticed that, in this application of the sequential strategy, the calibration subproblems have been defined from the analysis of the available experiments, allowing to select the parameters of model that would be used to describe the phenomena of interest. Furthermore, it should be emphasized that standard 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.

## 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 known to generate the synthetic data and are estimated with the proposed calibration approach.

A multivariate Gaussian distribution, defined by a mean vector $\mu_{\text{exact}}$ and a covariance matrix $\Sigma_{\text{exact}}$ is chosen to generate the synthetic data. This set of parameters will be denoted $\Pi_{\text{exact}}$ in the following. Individual parameters $\theta_{i,\text{exact}}$ are sampled from this distribution to compute, for each of these individuals the elastic model outputs $\varepsilon_{i,\text{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 $\mathcal{E}(\hat{\Pi})$, defined as:

$$ 
\begin{aligned}
\mathcal{E}(\hat{\Pi}) &:= (\mathcal{E}(\hat{\Pi}_1), \dots, \mathcal{E}(\hat{\Pi}_\ell)) \text{ with } \mathcal{E}(\hat{\Pi}_u) \\
&:= \frac{|\Pi_{\text{exact},u} - \hat{\Pi}_u|}{\Pi_{\text{exact},u}}, \quad \forall u \in [1, \ell],
\end{aligned} 
\quad (28) 
$$

noting $\ell$ the number of distribution parameters and the KL-divergence [60], KL($f_{\Theta|\Pi_{\text{exact}}}, f_{\Theta|\hat{\Pi}}$), between the exact $f_{\Theta|\Pi_{\text{exact}}}$ and the estimated distribution $f_{\Theta|\hat{\Pi}}$:

$$ 
\begin{aligned}
& \text{KL}(f_{\Theta|\Pi_{\text{exact}}}, f_{\Theta|\hat{\Pi}}) \\
& := \int_{\mathbb{R}^d} f_{\Theta|\Pi_{\text{exact}}}(\theta) \ln \left( \frac{f_{\Theta|\Pi_{\text{exact}}}(\theta)}{f_{\Theta|\hat{\Pi}}(\theta)} \right) d\theta.
\end{aligned} 
\quad (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 $\Pi_{\text{exact}} = [\mu_{\text{exact}}, \Sigma_{\text{exact}}]$ and $\hat{\Pi} = [\hat{\mu}, \hat{\Sigma}]$ respectively, the Kullback-Leibler (KL) divergence can be
<table><caption>**Table 1.** For the nonlinear elastic and viscoelastic behavioral 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.</caption><thead><tr><th>Calibration method</th><th>Maximum number of population parameters to estimate for the 2D nonlinear elasticity</th><th>Maximum number of population parameters to estimate for the 2D-viscoelasticity behavior</th></tr></thead><tbody><tr><td>Design matrices with dependence structure</td><td>17</td><td>42</td></tr><tr><td>Design matrices without dependence structure</td><td>12</td><td>16</td></tr><tr><td>Sequential approach with dependence structure</td><td>9</td><td>27</td></tr><tr><td>Sequential approach without dependence structure</td><td>6</td><td>12</td></tr></tbody></table>

because of the small number of available individuals, striking a balance between the complexity of the modeling and the number of available individuals 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 $\omega_i$s) is also to calibrate. Here again the implementation of the sequential approach allows to reduce the number of residual models (and thus of parameters) to estimate. Indeed, with the sequential approach, a focus is made on certain types of specimen (for instance compression experiments) while when with design matrices, different types of experiments are considered (tension and compression experiments). Thus, under the reasonable hypothesis of at least one residual model for each type of experiments, there are more residual models (so more parameters) to calibrate when design matrices are employed than in the different stages of the implementation of the sequential approach. Consequently, it appears that the sequential strategy is particularly suited to reduce the complexity of the minimization problem. It should be noticed that, in this application of the sequential strategy, the calibration subproblems have been defined from the analysis of the available experiments, allowing to select the parameters of model that would be used to describe the phenomena of interest. Furthermore, it should be emphasized that standard 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.

## 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 known to generate the synthetic data and are estimated with the proposed calibration approach.

A multivariate Gaussian distribution, defined by a mean vector $\mu_{\text{exact}}$ and a covariance matrix $\Sigma_{\text{exact}}$ is chosen to generate the synthetic data. This set of parameters will be denoted $\Pi_{\text{exact}}$ in the following. Individual parameters $\theta_{i,\text{exact}}$ are sampled from this distribution to compute, for each of these individuals the elastic model outputs $\varepsilon_{i,\text{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 $\mathcal{E}(\hat{\Pi})$, defined as:

$$ 
\begin{aligned}
\mathcal{E}(\hat{\Pi}) &:= (\mathcal{E}(\hat{\Pi}_1), \dots, \mathcal{E}(\hat{\Pi}_\ell)) \text{ with } \mathcal{E}(\hat{\Pi}_u) \\
&:= \frac{|\Pi_{\text{exact},u} - \hat{\Pi}_u|}{\Pi_{\text{exact},u}}, \quad \forall u \in [1, \ell],
\end{aligned} 
\quad (28) $$

noting $\ell$ the number of distribution parameters and the KL-divergence [60], KL($f_{\Theta|\Pi_{\text{exact}}}, f_{\Theta|\hat{\Pi}}$), between the exact $f_{\Theta|\Pi_{\text{exact}}}$ and the estimated distribution $f_{\Theta|\hat{\Pi}}$:

$$ 
\begin{aligned}
& \text{KL}(f_{\Theta|\Pi_{\text{exact}}}, f_{\Theta|\hat{\Pi}}) \\
& := \int_{\mathbb{R}^d} f_{\Theta|\Pi_{\text{exact}}}(\theta) \ln \left( \frac{f_{\Theta|\Pi_{\text{exact}}}(\theta)}{f_{\Theta|\hat{\Pi}}(\theta)} \right) d\theta.
\end{aligned} 
\quad (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 $\Pi_{\text{exact}} = [\mu_{\text{exact}}, \Sigma_{\text{exact}}]$ and $\hat{\Pi} = [\hat{\mu}, \hat{\Sigma}]$ respectively, the Kullback-Leibler (KL) divergence can be
<table><caption>**Table 2.** Exact values of the parameters distribution.</caption><thead><tr><th></th><th>S<sub>11</sub><sup>0</sup> [MPa<sup>-1</sup>]</th><th>S<sub>1</sub><sup>C</sup> [MPa<sup>-1</sup>]</th><th>S<sub>1</sub><sup>T</sup> [MPa<sup>-1</sup>]</th><th>ν<sub>12</sub> [×10<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Mean</td><td>8.52 × 10<sup>-6</sup></td><td>1.55 × 10<sup>-5</sup></td><td>6.47 × 10<sup>-6</sup></td><td>33.1</td></tr><tr><td>Standard-deviation</td><td>1.31 × 10<sup>-7</sup></td><td>1.73 × 10<sup>-6</sup></td><td>4.19 × 10<sup>-7</sup></td><td>0.660</td></tr></tbody></table>

estimated in closed form and refers to [61]

$$ \text{KL}(f_{\Theta|\Pi_{\text{exact}}}, f_{\Theta|\hat{\Pi}}) = \frac{1}{2} \left( \text{tr}(\hat{\Sigma}^{-1} \Sigma_{\text{exact}}) + (\boldsymbol{\mu}_{\text{exact}} - \hat{\boldsymbol{\mu}})^{\top} \hat{\Sigma}^{-1} (\boldsymbol{\mu}_{\text{exact}} - \hat{\boldsymbol{\mu}}) + \ln \left( \frac{|\hat{\Sigma}|}{|\Sigma_{\text{exact}}|} \right) - d \right). \quad (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 $\hat{\theta}_i$ or estimated output $\mathcal{F}(t_i, \hat{\theta}_i)$, the averaged relative error on the individual parameters $e(\hat{\theta}_1, \dots, \hat{\theta}_n)$ is defined as:

$$ \begin{gathered} e(\hat{\theta}_1, \dots, \hat{\theta}_n) := \frac{1}{n} \sum_{i=1}^{n} e_i(\hat{\theta}_i) \\ \text{with } e_i(\hat{\theta}_i) = (e_i(\hat{\theta}_i^1), \dots, e_i(\hat{\theta}_i^d)) \quad (31) \\ e_i(\hat{\theta}_i^a) = \frac{|\theta_{i,\text{exact}}^a - \hat{\theta}_i^a|}{\theta_{i,\text{exact}}^a} \quad \forall a \in [1, d], \end{gathered} $$

and the average error $d(\hat{\theta}_1, \dots, \hat{\theta}_n)$ between the exact and estimated strains for each coordinate

$$ \begin{gathered} d(\hat{\theta}_1, \dots, \hat{\theta}_n) := (d_1(\hat{\theta}_1, \dots, \hat{\theta}_n), \dots, d_m(\hat{\theta}_1, \dots, \hat{\theta}_n)) \\ \text{with } d_k(\hat{\theta}_1, \dots, \hat{\theta}_n) \\ := \frac{1}{n} \sum_{i=1}^{n} \frac{1}{N_i} \|\varepsilon_{i, \text{exact}}^k - \mathcal{F}^k(t_i, \hat{\theta}_i)\|_2^2 \quad \forall k \in [1, m], \quad (32) \end{gathered} $$

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_{\Theta|\Pi_{\text{exact}}}$. The above indicators are averaged over the repetitions and their dispersion is characterized by the coefficient of variation $\text{COV}(X)$ defined for any non-zero quantity $X$ as $\text{COV}(X) = \frac{\sqrt{\hat{V}(X)}}{\hat{M}(X)}$ with $\hat{M}(X)$ the empirical mean of $X$ and $\hat{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:

$$ \forall i \in [1, n], \quad \Omega_i = \Omega. $$

## 4.4 Generating virtual data

The exact means and standard-deviations of the parameters for the joint distribution $f_{\Theta|\Pi_{\text{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 $S_{11}^0$ and $\nu_{12}$ are chosen to fit with the values of dispersion given by Rollet [62]. Several correlations scenarios are considered: independence (between $S_1^C$ and $S_1^T$), positive correlation: $\rho(S_{11}^0, S_1^T) = 0.809$, $\rho(\nu_{12}, S_1^C) = 0.730$ and negative correlations $\rho(S_1^T, \nu_{12}) = -0.500$, $\rho(S_{11}^0, \nu_{12}) = -0.741$, $\rho(S_1^C, S_{11}^0) = -0.414$, noting $\rho$ the correlation. For each of the 20 repetitions, 50 independent identically distributed vectors of parameters corresponding to 50 virtual individuals are sampled from $f_{\Theta|\Pi_{\text{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 $(\varepsilon_{i,\text{exact}})_{i \in [1,n]}$ are computed. An heteroscedastic noise is added to the experimental data:

$$ \begin{gathered} \varepsilon_{i, \text{noisy}}^j = \varepsilon_{i, \text{exact}}^j \times (\mathbf{1}_{\mathbb{R}^2} + \tau \zeta) \\ \text{with } \zeta \sim \mathcal{N}(\mathbf{0}_{\mathbb{R}^2}, I_2) \text{ and } \tau = 0.02, \quad (33) \end{gathered} $$

with $\mathbf{1}_{\mathbb{R}^2}$ the unitary vector of $\mathbb{R}^2$, $\mathbf{0}_{\mathbb{R}^2}$ the zero vector of $\mathbb{R}^2$ and $I_2$ the identity matrix of $M_2(\mathbb{R})$. 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
<table><caption>**Table 2.** Exact values of the parameters distribution.</caption><thead><tr><th></th><th>S<sub>11</sub><sup>0</sup> [MPa<sup>-1</sup>]</th><th>S<sub>1</sub><sup>C</sup> [MPa<sup>-1</sup>]</th><th>S<sub>1</sub><sup>T</sup> [MPa<sup>-1</sup>]</th><th>ν<sub>12</sub> [×10<sup>-2</sup>]</th></tr></thead><tbody><tr><td>Mean</td><td>8.52 × 10<sup>-6</sup></td><td>1.55 × 10<sup>-5</sup></td><td>6.47 × 10<sup>-6</sup></td><td>33.1</td></tr><tr><td>Standard-deviation</td><td>1.31 × 10<sup>-7</sup></td><td>1.73 × 10<sup>-6</sup></td><td>4.19 × 10<sup>-7</sup></td><td>0.660</td></tr></tbody></table>

estimated in closed form and refers to [61]

$$ \text{KL}(f_{\Theta|\Pi_{\text{exact}}}, f_{\Theta|\hat{\Pi}}) = \frac{1}{2} \left( \text{tr}(\hat{\Sigma}^{-1} \Sigma_{\text{exact}}) + (\boldsymbol{\mu}_{\text{exact}} - \hat{\boldsymbol{\mu}})^{\top} \hat{\Sigma}^{-1} (\boldsymbol{\mu}_{\text{exact}} - \hat{\boldsymbol{\mu}}) + \ln \left( \frac{|\hat{\Sigma}|}{|\Sigma_{\text{exact}}|} \right) - d \right). \quad (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 $\hat{\theta}_i$ or estimated output $\mathcal{F}(t_i, \hat{\theta}_i)$, the averaged relative error on the individual parameters $e(\hat{\theta}_1, \dots, \hat{\theta}_n)$ is defined as:

$$ \begin{gathered} e(\hat{\theta}_1, \dots, \hat{\theta}_n) := \frac{1}{n} \sum_{i=1}^{n} e_i(\hat{\theta}_i) \\ \text{with } e_i(\hat{\theta}_i) = (e_i(\hat{\theta}_i^1), \dots, e_i(\hat{\theta}_i^d)) \quad (31) \\ e_i(\hat{\theta}_i^a) = \frac{|\theta_{i, \text{exact}}^a - \hat{\theta}_i^a|}{\theta_{i, \text{exact}}^a} \quad \forall a \in [1, d], \end{gathered} $$

and the average error $d(\hat{\theta}_1, \dots, \hat{\theta}_n)$ between the exact and estimated strains for each coordinate

$$ \begin{gathered} d(\hat{\theta}_1, \dots, \hat{\theta}_n) := (d_1(\hat{\theta}_1, \dots, \hat{\theta}_n), \dots, d_m(\hat{\theta}_1, \dots, \hat{\theta}_n)) \\ \text{with } d_k(\hat{\theta}_1, \dots, \hat{\theta}_n) \\ := \frac{1}{n} \sum_{i=1}^{n} \frac{1}{N_i} \|\varepsilon_{i, \text{exact}}^k - \mathcal{F}^k(t_i, \hat{\theta}_i)\|_2^2 \quad \forall k \in [1, m], \quad (32) \end{gathered} $$

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_{\Theta|\Pi_{\text{exact}}}$. The above indicators are averaged over the repetitions and their dispersion is characterized by the coefficient of variation $\text{COV}(X)$ defined for any non-zero quantity $X$ as $\text{COV}(X) = \frac{\sqrt{\hat{V}(X)}}{\hat{M}(X)}$ with $\hat{M}(X)$ the empirical mean of $X$ and $\hat{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:

$$ \forall i \in [1, n], \quad \Omega_i = \Omega. $$

## 4.4 Generating virtual data

The exact means and standard-deviations of the parameters for the joint distribution $f_{\Theta|\Pi_{\text{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 $S_{11}^0$ and $\nu_{12}$ are chosen to fit with the values of dispersion given by Rollet [62]. Several correlations scenarios are considered: independence (between $S_1^C$ and $S_1^T$), positive correlation: $\rho(S_{11}^0, S_1^T) = 0.809$, $\rho(\nu_{12}, S_1^C) = 0.730$ and negative correlations $\rho(S_1^T, \nu_{12}) = -0.500$, $\rho(S_{11}^0, \nu_{12}) = -0.741$, $\rho(S_1^C, S_{11}^0) = -0.414$, noting $\rho$ the correlation. For each of the 20 repetitions, 50 independent identically distributed vectors of parameters corresponding to 50 virtual individuals are sampled from $f_{\Theta|\Pi_{\text{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 $(\varepsilon_{i, \text{exact}})_{i \in [1, n]}$ are computed. An heteroscedastic noise is added to the experimental data:

$$ \begin{gathered} \varepsilon_{i, \text{noisy}}^j = \varepsilon_{i, \text{exact}}^j \times (1_{\mathbb{R}^2} + \tau\zeta) \\ \text{with } \zeta \sim \mathcal{N}(0_{\mathbb{R}^2}, I_2) \text{ and } \tau = 0.02, \end{gathered} \quad (33) $$

with $1_{\mathbb{R}^2}$ the unitary vector of $\mathbb{R}^2$, $0_{\mathbb{R}^2}$ the zero vector of $\mathbb{R}^2$ and $I_2$ the identity matrix of $M_2(\mathbb{R})$. 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
<figure><img src="image_7.png" /><figcaption>**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).</figcaption></figure>

<table><caption>**Table 3.** Search space for the population parameters.</caption><thead><tr><th>μ<sub>S<sub>11</sub></sub><sup>0</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</th><th>μ<sub>S<sub>1</sub></sub><sup>T</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</th><th>μ<sub>S<sub>1</sub></sub><sup>C</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</th><th>μ<sub>ν<sub>12</sub></sub> (10<sup>-2</sup>)</th></tr></thead><tbody><tr><td>[10.0, 110.0]</td><td>[50.0, 130]</td><td>[60.0, 250]</td><td>[15.0, 50.0]</td></tr><tr><td>σ<sub>S<sub>11</sub></sub><sup>0</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</td><td>σ<sub>S<sub>1</sub></sub><sup>T</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</td><td>σ<sub>S<sub>1</sub></sub><sup>C</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</td><td>σ<sub>ν<sub>12</sub></sub> (10<sup>-2</sup>)</td></tr><tr><td>[0.500, 4.50]</td><td>[0.500, 10.0]</td><td>[5.00, 30.0]</td><td>[0.100, 3.00]</td></tr><tr><td>ω<sub>long</sub><sup>T</sup> [×10<sup>-6</sup>]</td><td>ω<sub>trans</sub><sup>T</sup> [×10<sup>-5</sup>]</td><td>ω<sub>long</sub><sup>C</sup> [×10<sup>-5</sup>]</td><td>ω<sub>trans</sub><sup>C</sup> [×10<sup>-5</sup>]</td></tr><tr><td>[1.00, 10.0]</td><td>[3.16, 31.6]</td><td>[1.00, 10.0]</td><td>[3.16, 31.6]</td></tr></tbody></table>

calibration results with the tension experiments are investigated to check that the joint distribution of $S_{11}^0, \nu_{12}$ and $S_1^T$ 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.

### 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 $\pi$. 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 \times 1_{\mathbb{R}^p}$ with $p$ the dimension of $\Psi$. The population suggested by [44] corresponds to $4 + 3\ln(p)$ (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 minima, a multi-start version of the algorithm is chosen. 50 initialization points are considered, whose one corresponds 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
<figure><img src="image_7.png" /><figcaption>**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).</figcaption></figure>

<table><caption>**Table 3.** Search space for the population parameters.</caption><thead><tr><th>μ<sub>S<sub>11</sub></sub><sup>0</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</th><th>μ<sub>S<sub>1</sub></sub><sup>T</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</th><th>μ<sub>S<sub>1</sub></sub><sup>C</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</th><th>μ<sub>ν<sub>12</sub></sub> (10<sup>-2</sup>)</th></tr></thead><tbody><tr><td>[10.0, 110.0]</td><td>[50.0, 130]</td><td>[60.0, 250]</td><td>[15.0, 50.0]</td></tr><tr><td>σ<sub>S<sub>11</sub></sub><sup>0</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</td><td>σ<sub>S<sub>1</sub></sub><sup>T</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</td><td>σ<sub>S<sub>1</sub></sub><sup>C</sup> (10<sup>-7</sup>MPa<sup>-1</sup>)</td><td>σ<sub>ν<sub>12</sub></sub> (10<sup>-2</sup>)</td></tr><tr><td>[0.500, 4.50]</td><td>[0.500, 10.0]</td><td>[5.00, 30.0]</td><td>[0.100, 3.00]</td></tr><tr><td>ω<sub>long</sub><sup>T</sup> [×10<sup>-6</sup>]</td><td>ω<sub>trans</sub><sup>T</sup> [×10<sup>-5</sup>]</td><td>ω<sub>long</sub><sup>C</sup> [×10<sup>-5</sup>]</td><td>ω<sub>trans</sub><sup>C</sup> [×10<sup>-5</sup>]</td></tr><tr><td>[1.00, 10.0]</td><td>[3.16, 31.6]</td><td>[1.00, 10.0]</td><td>[3.16, 31.6]</td></tr></tbody></table>

calibration results with the tension experiments are investigated to check that the joint distribution of $S_{11}^0$, $\nu_{12}$ and $S_1^T$ 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.

### 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 $\pi$. 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 \times 1_{\mathbb{R}^p}$ with $p$ the dimension of $\Psi$. The population suggested by [44] corresponds to $4 + 3\ln(p)$ (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 minima, a multi-start version of the algorithm is chosen. 50 initialization points are considered, whose one corresponds 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
<table><caption>**Table 4.** Averaged calibrated mean, standard-deviations and correlations $\Pi_{1,mean}$ and averaged error $\mathcal{E}_{\text{mean}}(\hat{\Pi}_1)$ over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.</caption><thead><tr><th></th><th>S<sup>0</sup><sub>11</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th><th>ν<sub>12</sub> [×10<sup>-2</sup>]</th><th>S<sup>T</sup><sub>1</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Exact Means</td><td>85.2</td><td>33.1</td><td>64.7</td></tr><tr><td>Averaged Calibrated Means</td><td>85.3(0.305)</td><td>33.1(0.362)</td><td>64.8(0.872)</td></tr><tr><td>Averaged Errors</td><td>0.249(73.1)</td><td>0.306(64.9)</td><td>0.738(68.7)</td></tr><tr><td>On Calibrated Means (%)</td><td></td><td></td><td></td></tr><tr><td>Exact Standard-Deviations</td><td>1.31</td><td>0.660</td><td>4.19</td></tr><tr><td>Averaged Calibrated Standard Deviations</td><td>1.19(21.5)</td><td>0.672(11.7)</td><td>4.08(9.33)</td></tr><tr><td>Averaged Errors On Calibrated Standard-Deviations (%)</td><td>18.1(65.9)</td><td>8.31(105)</td><td>7.82(68.5)</td></tr><tr><td>Correlations</td><td>ρ(S<sup>0</sup><sub>11</sub>, ν<sub>12</sub>)</td><td>ρ(S<sup>0</sup><sub>11</sub>, S<sup>T</sup><sub>1</sub>)</td><td>ρ(S<sup>T</sup><sub>1</sub>, ν<sub>12</sub>)</td></tr><tr><td>Exact Correlations</td><td>-0.741</td><td>-0.500</td><td>0.809</td></tr><tr><td>Averaged Calibrated Correlations</td><td>-0.676(23.4)</td><td>-0.450(24.8)</td><td>0.851(16.3)</td></tr><tr><td>Averaged Errors On Calibrated Correlations (%)</td><td>17.9(81.4)</td><td>18.8(83.0)</td><td>15.9(53.3)</td></tr></tbody></table>

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μ].

### 4.5.2 Calibration of the joint distribution of $S_{11}^0, S_1^T$ and $\nu_{12}$

Given the modeling assumptions described in Sections 4.2 and 4.5, there are $p = 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 $\Sigma$ 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 $S_1^T$ and $\nu_{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 $S_{11}^0$ 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 $S_{11}^0$ mainly expresses at the beginning of the strain-stress curve, where the residuals have the lowest values, while $S_1^T$ expresses at the end of the experiment, for the highest residual values. Thus, function $g_i(\cdot)$ defined in equation (17) is necessarily more sensitive to the latter parameter than to $S_{11}^0$ 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 $S_{11}^0, S_1^T$ and $\nu_{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 $\zeta_{\text{exact}}$ (instead of the empirical COV), that is to say $\zeta_{\text{exact}} := \frac{\sqrt{\text{diag}(\Sigma_{\text{exact}})}}{\mu_{\text{exact}}}$ with $\mu_{\text{exact}}$ and $\Sigma_{\text{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 $\theta$.
<table><caption>**Table 4.** Averaged calibrated mean, standard-deviations and correlations $\Pi_{1,mean}$ and averaged error $\mathcal{E}_{\text{mean}}(\hat{\Pi}_1)$ over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.</caption><thead><tr><th></th><th>S<sup>0</sup><sub>11</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th><th>ν<sub>12</sub> [×10<sup>-2</sup>]</th><th>S<sup>T</sup><sub>1</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Exact Means</td><td>85.2</td><td>33.1</td><td>64.7</td></tr><tr><td>Averaged Calibrated Means</td><td>85.3(0.305)</td><td>33.1(0.362)</td><td>64.8(0.872)</td></tr><tr><td>Averaged Errors</td><td>0.249(73.1)</td><td>0.306(64.9)</td><td>0.738(68.7)</td></tr><tr><td>On Calibrated Means (%)</td><td></td><td></td><td></td></tr><tr><td>Exact Standard-Deviations</td><td>1.31</td><td>0.660</td><td>4.19</td></tr><tr><td>Averaged Calibrated Standard Deviations</td><td>1.19(21.5)</td><td>0.672(11.7)</td><td>4.08(9.33)</td></tr><tr><td>Averaged Errors On Calibrated Standard-Deviations (%)</td><td>18.1(65.9)</td><td>8.31(105)</td><td>7.82(68.5)</td></tr><tr><td>Correlations</td><td>ρ(S<sup>0</sup><sub>11</sub>, ν<sub>12</sub>)</td><td>ρ(S<sup>0</sup><sub>11</sub>, S<sup>T</sup><sub>1</sub>)</td><td>ρ(S<sup>T</sup><sub>1</sub>, ν<sub>12</sub>)</td></tr><tr><td>Exact Correlations</td><td>-0.741</td><td>-0.500</td><td>0.809</td></tr><tr><td>Averaged Calibrated Correlations</td><td>-0.676(23.4)</td><td>-0.450(24.8)</td><td>0.851(16.3)</td></tr><tr><td>Averaged Errors On Calibrated Correlations (%)</td><td>17.9(81.4)</td><td>18.8(83.0)</td><td>15.9(53.3)</td></tr></tbody></table>

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μ].

### 4.5.2 Calibration of the joint distribution of $S_{11}^0, S_1^T$ and $\nu_{12}$

Given the modeling assumptions described in Sections 4.2 and 4.5, there are $p = 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 $\Sigma$ 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 $S_1^T$ and $\nu_{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 $S_{11}^0$ 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 $S_{11}^0$ mainly expresses at the beginning of the strain-stress curve, where the residuals have the lowest values, while $S_1^T$ expresses at the end of the experiment, for the highest residual values. Thus, function $g_i(\cdot)$ defined in equation (17) is necessarily more sensitive to the latter parameter than to $S_{11}^0$ 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 $S_{11}^0, S_1^T$ and $\nu_{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 $\zeta_{\text{exact}}$ (instead of the empirical COV), that is to say $\zeta_{\text{exact}} := \frac{\sqrt{\text{diag}(\Sigma_{\text{exact}})}}{\mu_{\text{exact}}}$ with $\mu_{\text{exact}}$ and $\Sigma_{\text{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 $\theta$.
<table><caption>**Table 5.** Relative errors on the estimation of the individual parameters $e(\hat{\theta}_1, \dots, \hat{\theta}_n)$ averaged over 20 repetitions of the calibration process with 50 individual along with the spread $\zeta_{\text{exact}}$. The corresponding coefficients of variation in % are indicated between brackets.</caption><thead><tr><th></th><th>S<sup>0</sup><sub>11</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th><th>ν<sub>12</sub> [×10<sup>-2</sup>]</th><th>S<sup>T</sup><sub>1</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Relative error (%)</td><td>0.616(41.0)</td><td>0.613(57.4)</td><td>1.6(135)</td></tr><tr><td>Spread ζ<sub>exact</sub> (%)</td><td>1.53</td><td>2.00</td><td>6.48</td></tr></tbody></table>

<table><caption>**Table 6.** Error in model space $d(\hat{\theta}_1, \dots, \hat{\theta}_n)$ for each strain component.</caption><thead><tr><th>Strain Component</th><th>ε<sub>11</sub></th><th>ε<sub>22</sub></th></tr></thead><tbody><tr><td>Order of magnitude of the mean strain value</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Distance in model space</td><td>8.67 × 10<sup>-9</sup></td><td>8.66 × 10<sup>-10</sup></td></tr></tbody></table>

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 outweigh those of the linear behavior monitored by $S_{11}^0$. Yet, this statement seems to be in contradiction with the level of error. Meanwhile, a comparison to the exact spread values $\zeta_{\text{exact}}$ can bring some elements of justification. Indeed, this variable is the most concentrated around its mean value, especially compared to $S_1^T$. Thus, within $[\mu_{S_{11}^0} - 3\sigma_{S_{11}^0}, \mu_{S_{11}^0} + 3\sigma_{S_{11}^0}]$, (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 $S_{11}^0, S_1^T$ and $\nu_{12}$ either at the population or individual levels. Consequently, the first stage of the methodology is completed.

### 4.5.3 Calibration of the joint distribution of $S_{11}^0, S_1^C$ and $\nu_{12}$

In the second stage of the methodology, the objective is to identify the joint distribution of $S_{11}^0, S_1^C$ and $\nu_{12}$. Furthermore, the parameters of the joint distribution of $\nu_{12}$ and $S_{11}^0: \mu_{\nu_{12}}, \mu_{S_{11}^0}, \mathbb{V}(\nu_{12}), \mathbb{V}(S_{11}^0)$ and $\rho(\nu_{12}, S_{11}^0)$, already estimated in the previous step, will be limited to a trust

region corresponding to $\mathcal{R}_1 = [0.8\Pi_1, 1.2\Pi_1]$. Considering 20% of variation both enables a small exploration around the calibrated value of $\Pi_1$ and to remain concentrated around them. To bound the correlation parameter, $\rho(\nu_{12}, S_{11}^0)$, 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 $\mathcal{R}_1$.

First, one can notice that the distribution parameters involving $S_1^C$ (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 $S_1^C$ 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 $S_{11}^0$ and $\nu_{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 $S_{11}^0$ (the averaged error increases by 30%). Indeed, both in terms of averaged values or averaged error on the population parameters. The correlation between $S_{11}^0$ and $\nu_{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 (KL) divergence between the exact joint distribution of $S_{11}^0, \nu_{12}$ ($f_{\text{exact}}(S_{11}^0, \nu_{12})$) and the calibrated distribution ($f_{\text{calibrated}}(S_{11}^0, \nu_{12})$) in the two stages as depicts Figure 7 for the 20 repetitions. The main interest of such a quantity is that it takes into account all aspects within the same quantity, which sounds more appropriate as all parameters were calibrated jointly.

First, one can notice that the distribution parameters involving $S_1^C$ (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
<table><caption>**Table 5.** Relative errors on the estimation of the individual parameters $e(\hat{\theta}_1, \dots, \hat{\theta}_n)$ averaged over 20 repetitions of the calibration process with 50 individual along with the spread $\zeta_{\text{exact}}$. The corresponding coefficients of variation in % are indicated between brackets.</caption><thead><tr><th></th><th>S<sup>0</sup><sub>11</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th><th>ν<sub>12</sub> [×10<sup>-2</sup>]</th><th>S<sup>T</sup><sub>1</sub> [×10<sup>-6</sup>MPa<sup>-1</sup>]</th></tr></thead><tbody><tr><td>Relative error (%)</td><td>0.616(41.0)</td><td>0.613(57.4)</td><td>1.6(135)</td></tr><tr><td>Spread ζ<sub>exact</sub> (%)</td><td>1.53</td><td>2.00</td><td>6.48</td></tr></tbody></table>

<table><caption>**Table 6.** Error in model space $d(\hat{\theta}_1, \dots, \hat{\theta}_n)$ for each strain component.</caption><thead><tr><th>Strain Component</th><th>ε<sub>11</sub></th><th>ε<sub>22</sub></th></tr></thead><tbody><tr><td>Order of magnitude of the mean strain value</td><td>10<sup>-3</sup></td><td>10<sup>-3</sup></td></tr><tr><td>Distance in model space</td><td>8.67 × 10<sup>-9</sup></td><td>8.66 × 10<sup>-10</sup></td></tr></tbody></table>

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 outweigh those of the linear behavior monitored by $S_{11}^0$. Yet, this statement seems to be in contradiction with the level of error. Meanwhile, a comparison to the exact spread values $\zeta_{\text{exact}}$ can bring some elements of justification. Indeed, this variable is the most concentrated around its mean value, especially compared to $S_1^T$. Thus, within $[\mu_{S_{11}^0} - 3\sigma_{S_{11}^0}, \mu_{S_{11}^0} + 3\sigma_{S_{11}^0}]$, (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 $S_{11}^0, S_1^T$ and $\nu_{12}$ either at the population or individual levels. Consequently, the first stage of the methodology is completed.

### 4.5.3 Calibration of the joint distribution of $S_{11}^0, S_1^C$ and $\nu_{12}$

In the second stage of the methodology, the objective is to identify the joint distribution of $S_{11}^0, S_1^C$ and $\nu_{12}$. Furthermore, the parameters of the joint distribution of $\nu_{12}$ and $S_{11}^0: \mu_{\nu_{12}}, \mu_{S_{11}^0}, \mathbb{V}(\nu_{12}), \mathbb{V}(S_{11}^0)$ and $\rho(\nu_{12}, S_{11}^0)$, already estimated in the previous step, will be limited to a trust

region corresponding to $\mathcal{R}_1 = [0.8\Pi_1, 1.2\Pi_1]$. Considering 20% of variation both enables a small exploration around the calibrated value of $\Pi_1$ and to remain concentrated around them. To bound the correlation parameter, $\rho(\nu_{12}, S_{11}^0)$, 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 $\mathcal{R}_1$.

First, one can notice that the distribution parameters involving $S_1^C$ (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 $S_1^C$ 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 $S_{11}^0$ and $\nu_{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 $S_{11}^0$ (the averaged error increases by 30%). Indeed, both in terms of averaged values or averaged error on the population parameters. The correlation between $S_{11}^0$ and $\nu_{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 (KL) divergence between the exact joint distribution of $S_{11}^0, \nu_{12}$ ($f_{\text{exact}}(S_{11}^0, \nu_{12})$) and the calibrated distribution ($f_{\text{calibrated}}(S_{11}^0, \nu_{12})$) in the two stages as depicts Figure 7 for the 20 repetitions. The main interest of such a quantity is that it takes into account all aspects within the same quantity, which sounds more appropriate as all parameters were calibrated jointly.

First, one can notice that the distribution parameters involving $S_1^C$ (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
<table><caption>**Table 7.** Averaged calibrated mean, standard-deviations and correlations $\Pi_{2,mean}$ and averaged error $\mathcal{E}_{\text{mean}}(\hat{\Pi}_2)$ over 20 repetitions of the calibration process with samples of 50 individuals each. The coefficients of variation COV in % are indicated between brackets.</caption><thead><tr><th></th><th>$S_{11}^0$ [$\times 10^{-6}\text{MPa}^{-1}$]</th><th>$\nu_{12}$ [$\times 10^{-2}$]</th><th>$S_1^C$ [$\times 10^{-6}\text{MPa}^{-1}$]</th></tr></thead><tbody><tr><td>Exact Means</td><td>85.2</td><td>33.1</td><td>155.2</td></tr><tr><td>Averaged Calibrated Means</td><td>85.4(0.424)</td><td>33.1(0.339)</td><td>155.4(1.89)</td></tr><tr><td>Averaged Errors On Calibrated Means (%)</td><td>0.376(76.5)</td><td>0.273(74.6)</td><td>1.52(74.0)</td></tr><tr><td>Exact Standard-Deviations</td><td>1.31</td><td>0.660</td><td>17.3</td></tr><tr><td>Averaged Calibrated Standard Deviations</td><td>1.01(28.6)</td><td>0.669(10.5)</td><td>17.9(7.93)</td></tr><tr><td>Averaged Errors On Calibrated Standard-Deviations (%)</td><td>26.3(67.5)</td><td>9.41(54.8)</td><td>6.68(84.9)</td></tr><tr><td>Correlations</td><td>$\rho(S_{11}^0, \nu_{12})$</td><td>$\rho(S_1^C, \nu_{12})$</td><td>$\rho(S_{11}^0, S_1^C)$</td></tr><tr><td>Exact Correlations</td><td>-0.741</td><td>0.729</td><td>-0.414</td></tr><tr><td>Averaged Calibrated Correlations</td><td>-0.726(18.6)</td><td>0.713(6.40)</td><td>-0.354(56.0)</td></tr><tr><td>Averaged Errors On Calibrated Correlations (%)</td><td>14.7(74.6)</td><td>4.91(90.4)</td><td>41.1(69.0)</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** KL-divergences of the calibrated and exact joint distribution $S_{11}^0$ and $\nu_{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 $f_{\Theta|\Pi}$ is.</figcaption></figure>

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 $S_1^C$ 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 $S_{11}^0$ and $\nu_{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 $S_{11}^0$ (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 $S_{11}^0$ and $\nu_{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 $S_{11}^0, \nu_{12}$ ($f_{\text{exact}}(S_{11}^0, \nu_{12})$) and the calibrated distribution ($f_{\text{calibrated}}(S_{11}^0, \nu_{12})$) in the two stages as depicted Figure 7 for the 20 repetitions.

In Figure 7, the deterioration of the joint distribution of $S_{11}^0$ and $\nu_{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