Robust Multiscale Identification of Apparent Elastic Properties at Mesoscale for Random Heterogeneous Materials with Multiscale Field Measurements

The aim of this work is to efficiently and robustly solve the statistical inverse problem related to the identification of the elastic properties at both macroscopic and mesoscopic scales of heterogeneous anisotropic materials with a complex microstructure that usually cannot be properly described in terms of their mechanical constituents at microscale. Within the context of linear elasticity theory, the apparent elasticity tensor field at a given mesoscale is modeled by a prior non-Gaussian tensor-valued random field. A general methodology using multiscale displacement field measurements simultaneously made at both macroscale and mesoscale has been recently proposed for the identification the hyperparameters of such a prior stochastic model by solving a multiscale statistical inverse problem using a stochastic computational model and some information from displacement fields at both macroscale and mesoscale. This paper contributes to the improvement of the computational efficiency, accuracy and robustness of such a method by introducing (i) a mesoscopic numerical indicator related to the spatial correlation length(s) of kinematic fields, allowing the time-consuming global optimization algorithm (genetic algorithm) used in a previous work to be replaced with a more efficient algorithm and (ii) an ad hoc stochastic representation of the hyperparameters involved in the prior stochastic model in order to enhance both the robustness and the precision of the statistical inverse identification method. Finally, the proposed improved method is first validated on in silico materials within the framework of 2D plane stress and 3D linear elasticity (using multiscale simulated data obtained through numerical computations) and then exemplified on a real heterogeneous biological material (beef cortical bone) within the framework of 2D plane stress linear elasticity (using multiscale experimental data obtained through mechanical testing monitored by digital image correlation).


Introduction
Within the framework of linear elasticity theory, the numerical modeling and simulation of heterogeneous materials with hierarchical complex random microstructure give rise to many scientific challenges.Their modeling is a topical issue with numerous applications in diverse material sciences, including for instance sedimentary rocks, natural composites, fiber-or nano-reinforced composites, some concretes and cementitious materials, some porous media, some living biological tissues, among many others [1].Although such materials are often considered and modeled as deterministic and homogeneous elastic media at macroscale in most practical applications, they are not only random and heterogeneous at microscale but they also usually cannot be explicitly described by any local morphological and mechanical properties of their constituents and easily reconstructed in a computational framework in the presence of multiple interfaces.The modeling and identification of their elastic properties at meso-or microscales have been the subject of many research works in recent decades.Nowadays, with the recent developments achieved around the construction of stochastic models for tensor-valued random elasticity fields and their experimental inverse identification using field imaging techniques, one of the most promising ways consists in introducing a prior stochastic model of the apparent elasticity tensor field of heterogeneous materials of the considered microstructure at a given mesoscale.Note that this mesoscopic scale allows the introduction of the spatial correlation length(s) of the microstructure, and that for materials with a hierarchical structure, such as cortical bone or tendon, different mesoscopic scales can be defined.Such a mesoscopic stochastic modeling of random heterogeneous elastic media can further be used to characterize the macroscopic mechanical properties in the context of the stochastic homogenization over a representative volume element (RVE) subdomain.This representative volume element should be, provided that it exists, sufficiently large compared to the microscale and sufficiently small compared to the macroscale.In the present probabilistic context, a major question concerns the statistical inverse identification of a prior stochastic model parameterized by a small or moderate number of hyperparameters using only partial and limited experimental data.

Overview of Inverse Methods for the Mechanical Characterization of Micro/Meso-Structural Properties
The inverse methods for the experimental identification of elastic properties of homogeneous or heterogeneous materials at macroscale and/or mesoscale have been the subject of numerous research works over the three past decades.The first methods related to the experimental characterization and description of random microstructural morphologies by using image analysis techniques have been introduced and developed by the end of the 1980s [2][3][4][5][6] for the numerical modeling and simulation of random microstructures made up with heterogeneous materials.Since the early 1990s, significant technological advances in the field of optical measuring instruments, such as digital cameras equipped with Charge-Coupled Device (CCD) or Complementary MetalOxideSemiconductor (CMOS) image sensors and microscope objectives, have widely contributed to the emergence of imaging techniques such as two-dimensional (2D) or three-dimensional (3D) digital image correlation (DIC) for identification purposes.DIC techniques [7][8][9] are now commonly used in solid mechanics and material sciences for experimental measurements of elastic displacement fields of samples under external loading [10][11][12][13][14][15][16] in order to identify mechanical properties of complex microstructures for heterogeneous materials [13,[17][18][19][20][21][22][23][24] with different classes of material symmetries.The recent milestones achieved around data acquisition systems and processing softwares for 3D images obtained for example by X-ray computed microtomography (µCT) [25][26][27][28][29][30], magnetic resonance imaging (MRI) [31][32][33][34], optical coherence tomography (OCT) [35][36][37][38][39] or any other non-invasive and non-destructive testing technique for the reconstruction of 3D images in high resolution, have allowed the development of three-dimensional measurements of displacement fields by digital volume correlation (DVC) [9,15,[40][41][42][43][44][45][46][47][48][49][50].Such 3D full-field measurements offer the potential of identifying stochastic models of 3D tensor-valued random elasticity fields at different scales for the mechanical characterization of 3D real microstructures made up of heterogeneous materials.

Multiscale Statistical Identification Method
In keeping with the aforementioned works, an innovative methodology has been recently proposed in Reference [140] for the multiscale statistical inverse identification of a prior stochastic model of the random apparent elasticity field at mesoscale for a heterogeneous anisotropic elastic microstructure.This multiscale identification procedure has been formulated within the framework of 3D linear elasticity theory under the following assumptions: (i) at macroscale, the elasticity tensor is deterministic and homogeneous and therefore independent of the spatial coordinates; (ii) at a given mesoscale, the tensor-valued random elasticity field is the restriction to a mesoscopic subdomain of a statistically homogeneous random field indexed by R 3 , allowing to be consistent with the assumption for the existence of a representative volume element in the framework of stochastic homogenization [68,128].
The proposed method allows for the multiscale inverse identification of (i) the tensor-valued random field that models the apparent elasticity tensor field at a given mesoscale, and (ii) the effective elasticity tensor at macroscale, for a heterogeneous anisotropic elastic material with a random microstructure whose morphological and mechanical properties cannot be properly described and reconstructed in a computational framework from the local topology and mechanical behavior of its constitutive phases.The prior stochastic model of the random elasticity field is constructed by using the MaxEnt principle [68,[72][73][74][75][76][77][78], initially derived within the general framework of information theory [141][142][143].We then obtain a second-order mean-square continuous non-Gaussian positive-definite symmetric real matrix-valued random field.In addition, an explicit algebraic representation has been established in Reference [144].Such a prior stochastic model of random elasticity field has been used, in particular, for stochastic boundary value problems, such as static linear elasticity problems [68,128,144].It is classically parameterized by a small or moderate number of scalar-, vector-and/or tensor-valued hyperparameters, namely the mean function of the random elasticity field, a dispersion coefficient controlling the level of statistical fluctuations of the random elasticity field around its mean function and spatial correlation lengths characterizing the spatial correlation structure of the random elasticity field.The statistical inverse problem for the identification of this prior stochastic model is formulated as a multi-objective optimization problem for which the optimal parameters are the optimal values of the hyperparameters of the stochastic model.However, within the framework of this identification methodology, it can be shown that the mean function of the random elasticity field cannot directly be identified using only the available experimental kinematic field measurements at mesoscale.The experimental values of the stress fields associated with the kinematic fields observed experimentally at mesoscale should also be known, but these values are not available in practice.Conversely, it can also be shown that the other hyperparameters (dispersion coefficient and spatial correlation lengths) controlling the statistical fluctuations of the random elasticity field cannot directly be identified using only the available experimental kinematic field measurements at macroscale.Consequently, such a statistical inverse identification procedure requires multiscale experimental field measurements that must be made simultaneously at both macroscopic and mesoscopic scales, since by assumption only a single specimen submitted to a given external loading at macroscale is experimentally tested.A stochastic homogenization method is then used to propagate the uncertainties at mesoscale towards the macroscale under the classical assumption of scale separation between macroscale and mesoscale, so that a sufficiently large mesoscopic subdomain can be defined within the macroscopic domain and considered as a representative volume element.However, it should be noted that it is not necessary for this representative volume element to be the same size as the mesoscopic domain(s) of observation on which the experimental measurements are performed.Thus, the multiscale statistical inverse problem is formulated as a multi-objective optimization problem that consists in minimizing a (vector-valued) multi-objective cost function defined by three numerical indicators corresponding to single-objective cost functions [140], namely (i) a macroscopic numerical indicator allowing the distance between the measured experimental fields and the computed numerical fields to be quantified at macroscale, (ii) a mesoscopic numerical indicator allowing the distance between the statistical fluctuations exhibited by the measured experimental fields and the ones exhibited by the computed numerical fields to be quantified at mesoscale, and (iii) a multiscale numerical indicator allowing the distance between the elasticity tensor at macroscale and the effective elasticity tensor constructed by computational stochastic homogenization of the random apparent elasticity field in a representative volume element at mesoscale.

Drawbacks and Limitations of the Multiscale Identification Method
The multiscale identification method proposed in Reference [140] has been first validated by numerical simulations on in silico materials and then successfully applied to the experimental characterization of the elastic properties of a biological tissue (beef cortical bone) within the framework of 2D plane stress linear elasticity from multiscale optical measurements of displacement fields performed at both macroscopic and mesoscopic scales on a single cortical bone specimen under static external loading at macroscale [145].Nevertheless, the proposed identification method has some drawbacks that limit its use.First, it should be noted that the cost functions introduced for the multi-objective optimization problem are not dedicated to a particular hyperparameter of the prior stochastic model of the random field to be identified.Therefore, the only approach considered for solving the multi-objective optimization problem was to use a global optimization algorithm (genetic algorithm) that belongs to the class of random search, genetic and evolutionary algorithms [146][147][148][149][150][151][152][153][154][155][156] to randomly explore the admissible set of hyperparameters.Despite a suitable parameterization (population size at each new generation, random generation of initial population, selection procedure for reproduction including crossover and mutation operators, elite count, stopping criteria, etc.) of the genetic algorithm used in Reference [140] and the use of parallel processing and computing, the computational cost for solving the multi-objective optimization problem is high.This is due in particular to the large stochastic dimension of the tensor-valued random elasticity field.
Secondly, during the validation and implementation of the multiscale identification method proposed in Reference [140], it was found that, for different mesoscopic domains of observation within the same macroscopic domain, the resolution of the multi-objective optimization problem led to different optimal values of hyperparameters from one domain to another.Indeed, the experimental field measurements over each mesoscopic domain of observation can be modeled as different random fields, and therefore the multi-objective cost function on each mesoscopic domain of observation is a deterministic function of these random fields.This explains why the statistics of the multi-objective cost function are different from one mesoscopic domain of observation to another.In Reference [140], the multi-objective cost function has been replaced by the statistical average of the multi-objective cost functions calculated over each of the mesoscopic domains of observation.

Improvements of the Multiscale Identification Method and Novelty of the Paper
In order to overcome the issues outlined above, this research work aims to present two major improvements of the methodology initially proposed in Reference [140] allowing the statistical inverse identification of the tensor-valued random elasticity field at mesoscale to be performed with a better computational efficiency, higher accuracy and improved robustness.First, we introduce an additional mesoscopic numerical indicator allowing the distance between the spatial correlation length(s) of the measured experimental kinematic fields and the one(s) of the computed numerical kinematic fields to be quantified at mesoscale, so that each hyperparameter of the prior stochastic model has its own dedicated single-objective cost function, thus allowing the time-consuming global optimization algorithm (genetic algorithm) used in Reference [140] to be avoided and replaced with a more efficient algorithm, such as a fixed-point iterative algorithm, for solving the underlying multi-objective optimization problem.Secondly, in the case where experimental field measurements are available on several mesoscopic domains of observation, we propose to not replace "naively" the multi-objective cost function by its empirical mean over all the mesoscopic domains of observation, but to consider a multi-objective optimization problem for each mesoscopic domain of observation.Thus, each mesoscopic domain of observation leads to a possible solution of the values of the hyperparameters.Each of these values is then considered as a realization of a random vector of hyperparameters whose prior stochastic model is constructed by using the MaxEnt principle, and whose hyperparameters can be determined by using the MLE method, in order to improve both the robustness and the accuracy of the inverse identification method of the prior stochastic model.

Outline of the Paper
The paper is organized as follows.Following this introduction, Section 2 presents the general assumptions for solving the underlying multiscale statistical inverse problem.Then, Section 3 is dedicated to the description of the multiscale experimental test configuration for obtaining experimental data at both macroscale and mesoscale.Section 4 describes the prior stochastic model of the fourth-order tensor-valued random elasticity field and its parameterization.Section 5 focuses on the objectives of the multiscale statistical inverse problem and the multiscale identification strategy.Next, Section 6 presents the construction of the macroscopic, mesoscopic and multiscale numerical indicators that are used for solving the multiscale statistical inverse problem as a multi-objective optimization problem.In this section, a focus is made on the improvements proposed by this paper in the definition of these numerical indicators with respect to the previous work presented in Reference [140].The multi-objective optimization problem is then set in Section 7 and some numerical methods for solving such a multi-objective problem are presented in Section 8. Section 9 discusses an improvement proposed in this paper for a robust identification when some experimental field measurements are available on several mesoscopic domains of observation.Section 10 presents a numerical validation of the proposed multiscale identification methodology on in silico test specimens within the framework of 3D linear elasticity under 2D plane stress assumption and in the general 3D case, for which the multiscale experimental data have been numerically simulated.Finally, Section 11 presents an experimental application to a real heterogeneous biological material constituted of beef cortical bone within the framework of linear elasticity under 2D plane stress assumption, for which the multiscale experimental data have been obtained from a single static uniaxial compression test performed on a specimen of beef femoral cortical bone and monitored by 2D digital image correlation at both macroscale and mesoscale.Lastly, Section 12 gives some conclusions and potential perspectives of this work.

Assumptions for Solving the Multiscale Statistical Inverse Problem
In the present work, we address the statistical inverse identification of the elastic properties for a complex microstructure made up of a heterogeneous anisotropic material and considered as a random linear elastic medium.In this section, we first state suitable assumptions for solving this multiscale statistical inverse problem.Within the framework of linear elasticity theory, probability theory and computational stochastic homogenization in micromechanics and multiscale mechanics of heterogeneous materials, the following assumptions related to scale separation, stationarity and ergodicity properties are introduced: • there exists a scale separation between macroscale and mesoscale, so that a mesoscopic subdomain can be defined and for which the dimensions are sufficiently large with respect to the size of the heterogeneities and sufficiently small with respect to the size of the macroscopic domain.Such a mesoscopic subdomain can then be considered as a representative volume element; • the random apparent elasticity tensor field at mesoscale is the restriction to one or more bounded mesoscopic subdomain(s) of a second-order stationary random field indexed by R 3 , and consequently the mean function of the random elasticity field at mesoscale is independent of the spatial coordinates; • the random apparent elasticity tensor field at mesoscale is ergodic in average in the mean-square sense, so that the homogenized elasticity tensor at macroscale calculated by stochastic homogenization of the random apparent elasticity field in a mesoscopic subdomain corresponding to a representative volume element can be considered as almost deterministic, in the sense that (i) its spatial average reaches an asymptotic convergence with a very high level of probability for a sufficiently large mesoscopic subdomain size, and therefore (ii) its level of statistical fluctuations around its mean function at macroscale can be considered as negligible, thus yielding a deterministic homogenized elasticity tensor at macroscale.
In this work, we focus on the class of heterogeneous materials that can be considered as random elastic media and for which the hypothesis stated on the scale separation between macroscale and mesoscale is verified.It should be noted that, if such a scale separation assumption was not satisfied, then the multiscale statistical inverse problem under consideration would be an ill-posed problem if only a single experimental field measurement at macroscale was available, because in this case the macroscopic elasticity (or compliance) tensor must be modeled by a random tensor and a single experimental measurement is not sufficient to identify its stochastic model.The proposed identification methodology is therefore not adapted to this case and would require several experimental field measurements at macroscale as well as modifications of the macroscopic and multiscale indicators introduced in Section 6, and also the introduction of additional numerical indicators at macroscale.Hereinafter, since the present identification methodology is developed within the framework of linear elasticity theory, we will use the terminology "strain field" to make reference to the "linearized strain field" for the sake of conciseness.

Multiscale Experimental Test Configuration
The difficulties related to the acquisition of the experimental measurements for the inverse identification procedure to be carried out are induced not only by the complex nature of the heterogeneous anisotropic elastic microstructure but also by the need to obtain multiscale kinematic field measurements at two different scales (macroscale and mesoscale) for a single test specimen under given static loading conditions through a multiscale DIC performed simultaneously at both macroscale and mesoscale.To overcome such difficulties, a suitable experimental protocol, including the preparation of the test specimen, the development of a measuring bench, the acquisition system of digital images and the DIC method, has been set up in Reference [145] for the acquisition of 2D multiscale optical measurements of displacement fields performed at both macroscale and mesoscale on a single beef cortical bone specimen submitted to a static vertical uniaxial compression test.Such a living biological tissue with a complex hierarchical microstructure is of particular interest in the present context of multiscale modeling and identification for random heterogeneous materials.The multiscale experimental test configuration is briefly recalled here.A sketch of the multiscale experimental configuration of the specimen at macroscale and mesoscale is represented in Figure 1.
∂Ω meso exp,q u meso exp,q and displacement field u meso exp,q measured in each mesoscopic domain of observation Ω meso exp,q , for q = 1, . . ., Q.
The test specimen has a cubic shape and is submitted to a simple external load.On the upper side of the specimen, a surface force field is applied, while the opposite side of the specimen is clamped.Then, during the same and unique experimental loading, the displacement fields at both macroscale and mesoscale are simultaneously measured, for instance in using two optical digital cameras equipped with CCD imaging sensors with different spatial resolutions for the simultaneous acquisition of displacement field optical measurements at both macroscopic and mesoscopic scales.The measurements are performed on the domain Ω macro exp at macroscale and on the domain Ω meso exp at mesoscale that are 2D or 3D parts of the specimen at macroscale and mesoscale, respectively.These domains can be 3D in the case of microtomography techniques for the acquisition of 3D experimental data, or they can be 2D in the case of digital camera techniques for the acquisition of 2D experimental data.Note that in case the dimensions of the mesoscopic domain of observation Ω meso exp are very small with respect to the dimensions of the macroscopic domain of observation Ω macro exp , then more information can be used by collecting additional experimental field measurements at mesoscale on Q non-overlapping mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q for which the relative mutual locations into the test specimen are not necessarily recorded.The experimental database is then constituted of the vector-valued experimental displacement fields u macro exp and u meso exp,1 , . . ., u meso exp,Q , respectively, at macroscale on Ω macro exp and at mesoscale on Ω meso exp,1 , . . ., Ω meso exp,Q .The experimental tensor-valued strain fields ε macro exp and ε meso exp,1 , . . ., ε meso exp,Q , respectively associated to the experimental displacement fields u macro exp and u meso exp,1 , . . ., u meso exp,Q , can be calculated by post-processing through interpolation techniques.

Prior Multiscale Stochastic Model and Its Hyperparameters
At the macroscale, the specimen under test is modeled as a deterministic homogeneous linear elastic medium for which the effective mechanical properties are represented by a deterministic model of the fourth-order elasticity tensor C macro (a) that is independent of spatial position x and parameterized by a vector a belonging to an admissible set A macro .The vector-valued parameter a is constituted of the algebraically independent coefficients spanning the macroscopic elasticity tensor C macro (a) having a given symmetry class induced by linear elastic material symmetries.At the mesoscale, the specimen under test is modeled as a random heterogeneous linear elastic medium for which the apparent mechanical properties are represented by a prior stochastic model of the fourth-order tensor-valued random elasticity field.In Reference [144], the ensemble SFE + of non-Gaussian second-order stationary random fields has been introduced and constructed in using the theory of information, the MaxEnt principle and the theory of random matrices.Such a family of tensor-valued random fields is completely parameterized by the values of their mean function, a dispersion coefficient usually denoted as δ, and d n(n + 1)/2 = (d 3 (d + 1) 2 + 2 d 2 (d + 1))/8 = 63 possibly different spatial correlation lengths, with d = 3 and n = d(d + 1)/2 = 6 in 3D linear elasticity (see References [128,144] for a definition of the spatial correlation lengths of a random field).All these parameters are independent of the spatial position x since every tensor-valued random field in SFE + is second-order stationary on R 3 by construction.In addition, the dispersion coefficient δ introduced in Reference [144] is such that where n = d(d + 1)/2 = 6 with d = 3 in 3D linear elasticity.Hence, any tensor-valued random field in SFE + has no statistical fluctuations when δ = 0 and consequently its values are almost surely (a.s.) equal to its mean function.In addition, the level of statistical fluctuations of any tensor-valued random field in SFE + increases with the value of δ.Consequently, the highest statistical fluctuations are obtained when δ = δ sup .Ensemble SFE + has been especially constructed in Reference [144] for offering a prior stochastic model that can be used for modeling the tensor-valued apparent elasticity (or compliance) fields at mesoscale.Consequently, in this paper, we will use the same approach and the prior stochastic model of the elasticity tensor field C meso (resp.the compliance tensor field S meso ) will be defined as the restriction to a given bounded subdomain in R 3 of a random tensor field belonging to SFE + and indexed by R 3 .
The prior stochastic model of C meso or S meso can then be deduced from each other by inverse of each other.In this work, we will only consider the special case for which the spatial correlation structure of C meso (resp.S meso ) is defined by only 3 (instead of 63) different values 1 , 2 , 3 for the spatial correlation lengths and consequently some of the 63 spatial correlation lengths are mutually equal to each other.Furthermore, the mean function of C meso (resp.S meso ) can be represented by a set of n sym n(n + 1)/2 parameters h 1 , . . ., h n sym that might have or not physical meaning in mechanical engineering such as Young's moduli, Poisson's ratios, bulk and shear moduli, and so forth (see for instance Section 10).Finally, the hyperparameters of the prior stochastic model of C meso (resp.S meso ) are δ, 1 , 2 , 3 and h 1 , . . ., h n sym that can be gathered into the vector-valued hyperparameter b = (δ, , h) in which = ( 1 , 2 , 3 ) and h = (h 1 , . . ., h n sym ).Hereinafter, the set of all the admissible values of vector h is denoted by H meso and the admissible set of vector b is denoted by B meso .

Objectives of the Multiscale Statistical Inverse Problem
The deterministic model of C macro (a) at macroscale and the prior stochastic model of C meso (b) at mesoscale have to be identified by calculating the optimal values a macro and b meso of the vector-valued parameter a ∈ A macro and the vector-valued hyperparameter b ∈ B meso , respectively, according to the experimental kinematic field measurements available at both macroscale and mesoscale.While the vector-valued parameter a can completely be identified by solving a usual deterministic inverse problem using only the available experimental field measurements at macroscale, the vector-valued hyperparameter b = (δ, , h) cannot directly be identified by solving a statistical inverse problem using only the available experimental field measurements at mesoscale.More precisely, the dispersion parameter δ and the vector of spatial correlation lengths require only experimental field measurements at mesoscale to be identified, whereas the vector h requires additional experimental field measurements at macroscale to be identified.Indeed, the hyperparameters δ and controlling respectively the level of statistical fluctuations and the spatial correlation structure of the random elasticity field require experimental field measurements with a sufficiently fine spatial resolution to be identified, while the hyperparameters h representing the mean elasticity field would require the experimental values of the stress fields associated with the kinematic (displacement or strain) fields observed experimentally at mesoscale to be identified, but these values are not available in practice.The complete statistical information on random field C meso (b) must then be transferred to the macroscale in order to identify its mean function C meso using the available experimental field measurements at macroscale.A natural choice for such a transfer of information consists in computing the effective elasticity tensor C eff (b) by a computational stochastic homogenization method and in comparing it with the previously identified elasticity tensor C macro (a).Thus, unlike the vector-valued parameter a, the vector-valued hyperparameter b requires multiscale experimental field measurements (at macroscale and mesoscale) to be completely identified, thus leading to a challenging multiscale statistical inverse problem to be solved.Since by assumption only a single specimen is experimentally tested under a given static external loading applied at macroscale, the experimental field measurements must be performed simultaneously at both macroscale and mesoscale on the single test specimen, but they do not need to be performed on the whole domain of the specimen.

Strategy for Solving the Multiscale Statistical Inverse Problem
Due to the major difficulties stated above and induced by the complexity of the challenging multiscale statistical inverse problem to be solved, a first complete methodology concerning such a multiscale identification has been recently proposed in Reference [140], in which a multiscale statistical inverse identification strategy is introduced and developed for an elastic microstructure with heterogeneous anisotropic statistical fluctuations within the framework of 3D linear elasticity theory.The proposed strategy allows for the identification of (i) the optimal value a macro of vector-valued parameter a, and (ii) the optimal value b meso of vector-valued hyperparameter b, by using the experimental displacement field measurements at both macroscale and mesoscale.The multiscale experimental identification methodology originally developed in Reference [140] consists in introducing and constructing three different numerical indicators allowing the multiscale statistical inverse problem to be formulated as a multi-objective optimization problem.In the present work, we develop an improved multiscale experimental identification methodology involving four numerical indicators that are sensitive to the variation of the parameters and hyperparameters to be identified, which are: 1.A macroscopic numerical indicator J macro (a), dedicated to the identification of parameter a, that allows for quantifying the distance between the experimental strain field ε macro exp associated to the experimental displacement field u macro exp measured at macroscale in the macroscopic domain Ω macro exp and the strain field ε macro (a) associated to the displacement field u macro (a) computed from a deterministic homogeneous linear elasticity boundary value problem (with both Dirichlet and Neumann boundary conditions) that models the experimental test configuration at macroscale and involves the unknown deterministic elasticity tensor C macro (a); 2. A mesoscopic numerical indicator J meso δ (b), dedicated to the identification of hyperparameter δ, that allows for quantifying the distance between a pseudo-dispersion coefficient δ ε exp modeling the level of spatial fluctuations of the experimental strain field ε meso exp associated to the experimental displacement field u meso exp measured at mesoscale in a mesoscopic domain of observation Ω meso exp , and a random pseudo-dispersion coefficient D E (b) representing the level of statistical fluctuations of the random strain field E meso (b) associated to the random displacement field U meso (b) computed from a stochastic heterogeneous linear elasticity boundary value problem (with only Dirichlet boundary conditions) that models the experimental test configuration at mesoscale and involves the random elasticity tensor field C meso (b) with an unknown level of statistical fluctuations δ that must be identified; 3. Another mesoscopic numerical indicator J meso (b), dedicated to the identification of hyperparameter = ( 1 , 2 , 3 ), that allows for quantifying the distance between the 3 different pseudo-spatial correlation lengths ε exp,1 , ε exp,2 , ε exp,3 of the experimental strain field ε meso exp in each spatial direction, measured at mesoscale in a mesoscopic domain of observation Ω meso exp , and the 3 pseudo-spatial correlation lengths of the random strain field E meso (b) in each spatial direction, computed from the same mesoscopic stochastic boundary value problem as for J meso δ (b) for which the random elasticity tensor field C meso (b) has a spatial correlation structure induced and characterized by an unknown vector of spatial correlation lengths = ( 1 , 2 , 3 ) that must be identified; 4. A multiscale numerical indicator J multi h (a, b), dedicated to the identification of hyperparameter h, that allows for quantifying the distance between the homogeneous deterministic elasticity tensor C macro (a) at macroscale and the effective elasticity tensor C eff (b) resulting from a computational stochastic homogenization in a representative volume element Ω RVE at mesoscale of the random elasticity tensor field C meso (b) whose mean function C meso is unknown and must be identify.
The multiscale statistical inverse problem then consists in identifying the optimal values a macro and b meso of the parameters a in A macro and hyperparameters b in B meso , respectively, by solving a multi-objective optimization problem that consists in minimizing the (vector-valued) multi-objective cost function J (a, b) = J macro (a), J meso δ (b), J meso (b), J multi h (a, b) involving the four aforementioned numerical indicators.However, for further computational savings, the multi-objective optimization problem can be decomposed into (i) a single-objective optimization problem that consists in minimizing J macro (a) for identifying the optimal vector-valued parameter a macro using only the experimental field measurements at macroscale, and (ii) a multi-objective optimization problem that consists in minimizing J meso (b) = J meso δ (b), J meso (b), J multi h (a macro , b) for identifying the optimal vector-valued hyperparameter b meso using the experimental field measurements at mesoscale and exploiting the optimal vector-valued parameter a macro previously identified at step (i).

Construction of the Numerical Indicators for Solving the Multiscale Statistical Inverse Problem
In this section, the construction of the macroscopic, mesoscopic and multiscale numerical indicators for solving the multiscale statistical inverse problem is presented., so that there is no rigid body motion during the test.Within the context of linear elasticity theory, the deterministic boundary value problem at macroscale consists in finding the vector-valued displacement field u macro and the associated tensor-valued Cauchy stress field σ macro satisfying the following equilibrium equations, stress-strain constitutive equation and Neumann and Dirichlet boundary conditions in which div denotes the divergence operator of a second-order tensor-valued field with respect to x, the colon symbol : denotes the classical twice contracted tensor product, n macro is the unit normal vector to ∂Ω macro pointing outward Ω macro and ε macro is the classical tensor-valued strain field associated to displacement field u macro and defined by in which ε denotes the deterministic linear operator mapping the displacement field to the corresponding linearized strain field, the superscript T denotes the transpose operator and ∇ denotes the gradient operator of a vector-valued field with respect to x. Recall that, as the material is assumed to be deterministic and homogeneous at macroscale, the unknown fourth-order deterministic elasticity tensor C macro (a) involved in constitutive Equation ( 3) is independent of x and parameterized by a parameter a belonging to an admissible set A macro depending on the considered material symmetry class.A sketch of the deterministic boundary value problem at macroscale is represented in Figure 2a.

Stochastic Mesoscopic Boundary Value Problem for the Mesoscopic Indicators
At mesoscale, the stochastic boundary value problem modeling the experimental test configuration described in Section 3 is written over an open bounded domain Ω meso ⊂ R 3 with mesoscopic dimensions.A given domain of observation Ω meso exp corresponds to one given 2D or 3D part Ω meso obs of Ω meso .Within the context of linear elasticity theory, the stochastic boundary value problem at mesoscale consists in finding the vector-valued random displacement field U meso and the associated tensor-valued random Cauchy stress field Σ meso satisfying the following equilibrium equations, stress-strain constitutive equation and Dirichlet boundary conditions where E meso is the tensor-valued random strain field associated to random displacement field U meso and defined by Note that non-homogeneous Dirichlet boundary conditions (9) are prescribed on the whole boundary ∂Ω meso of Ω meso , which correspond to the displacement field u meso exp that is experimentally measured over a given domain of observation Ω meso exp on the test specimen at mesoscale.Note also that (8) can equivalently be rewritten as where S meso (b) = (C meso (b)) −1 is the random compliance tensor field of the considered material at mesoscale.For some linear elasticity problems, such as with 2D plane stress assumption, constitutive Equation ( 11) is more appropriate than (8).A sketch of the stochastic boundary value problem at mesoscale is represented in Figure 2b.

Macroscopic Numerical Indicator
Within the context of inverse identification, the optimal identified value a macro of parameter a can be determined by exploiting the sensitivity of the model strain field ε macro with respect to a and using the experimental strain field ε macro exp , which is obtained in Ω macro exp but can be rewritten in Ω macro obs , through the introduction of a macroscopic numerical indicator J macro (a) defined for any vector a ∈ A macro by where |Ω macro obs | denotes the measure of domain Ω macro obs and • F denotes the Frobenius (or Hilbert-Schmidt) norm.Macroscopic numerical indicator J macro (a) allows for quantifying the spatial average over the macroscopic domain Ω macro obs of the distance between the model strain field ε macro (a) and the experimental strain field ε macro exp at macroscale.The optimal vector-valued parameter a macro can then be identified by minimizing J macro (a) over all vector-valued parameter a in A macro , provided that the model strain field ε macro (a) computed by solving the deterministic boundary value problem ( 2)-( 6) is sufficiently sensitive to parameter a.

Mesoscopic and Multiscale Numerical Indicators
Within the context of statistical inverse identification, the optimal identified values b meso = (δ meso , meso , h meso ) of b = (δ, , h) can be determined by exploiting the sensitivity of some quantities of interest of the stochastic boundary value problem ( 7)- (10) with respect to δ, = ( 1 , 2 , 3 ) and h, respectively, and using their counterparts coming from the experimental measurements through the introduction of two mesoscopic numerical indicators J meso δ (b) and J meso (b) and one multiscale numerical indicator J multi h (a, b).

Mesoscopic Numerical Indicator Associated to the Dispersion Parameter
A first mesoscopic numerical indicator J meso δ (b) is introduced to identify the dispersion parameter δ controlling the level of statistical fluctuations of random elasticity field C meso (b) at mesoscale and defined for any vector b ∈ B meso by where E denotes the mathematical expectation, D E (b) is a positive-valued random variable that models the random level of spatial fluctuations of the random solution obtained by solving the stochastic boundary value problem ( 7)-( 10) at mesoscale and where δ ε exp is its counterpart for the experimental test specimen at mesoscale, such that where E meso (b) and ε meso exp are the spatial averages of random strain field E meso (b) and experimental strain field ε meso exp , respectively, and where where |Ω meso obs | denotes the measure of domain Ω meso obs .Note that it can easily be shown that E meso (b) = ε meso exp for all b ∈ B meso a.s. and consequently E meso (b) is a deterministic tensor.Also, since random strain field E meso (b) is a priori nor statistically homogeneous neither ergodic in average, E meso (b) does not correspond to the statistical mean function of E meso (b) and therefore V E (b) (resp.D E (b)) does not correspond to the variance (resp.dispersion coefficient) of E meso (b).The mesoscopic numerical indicator J meso δ (b) defined by ( 13) allows for quantifying the relative distance between the statistical mean value of D E (b) and its experimental observation δ ε exp .It should also be noted that a mesoscopic numerical indicator similar to this one was introduced in Reference [140], but with different expressions than that of ( 13), ( 15) and ( 16) for the definitions of J meso δ (b) and V E (b), respectively.

Mesoscopic Numerical Indicator Associated to the Spatial Correlation Lengths
A second mesoscopic numerical indicator J meso (b) is introduced to identify the vector of spatial correlation lengths = ( 1 , 2 , 3 ) characterizing the spatial correlation structure of random elasticity field C meso (b) (or random compliance field S meso (b)) and defined for any vector b ∈ B meso by where L E α (b) is a positive-valued random variable that models the spatial correlation length along the α-th spatial direction (relative to the spatial coordinate x α ) characterizing the spatial correlation structure of the statistical fluctuations of random strain field E meso (b) and where ε exp,α is its observation for the experimental test specimen at mesoscale.Usual signal processing methods (such as the periodogram method) are used for estimating L E α (b) and ε exp,α by considering the approximation that they are independent of x which is not the case since E meso (b) and ε meso exp are usually not statistically homogeneous because of the non-homogeneous Dirichlet boundary conditions (9) involving the experimental displacement field u meso exp on ∂Ω meso .The mesoscopic numerical indicator J meso (b) defined by (17) allows for quantifying the relative distance between the statistical mean values of

Multiscale Numerical Indicator Associated to Computational Stochastic Homogenization
A multiscale numerical indicator J multi h (a, b) is introduced to identify the mean function C meso (b) of the random elasticity field C meso (b) at mesoscale and defined for any vector a ∈ A macro and any vector b ∈ B meso by where C eff (b) is the effective elasticity tensor constructed by computational stochastic homogenization of C meso (b) in an open bounded mesoscopic domain Ω RVE , which is assumed to be a representative volume element.It should be noted that, under scale separation assumption, C eff (b) is actually a random tensor for which the level of statistical fluctuations tends to zero when the size of domain Ω RVE tends to infinity [68,128,131].This is the reason why the statistical mean value E{C eff (b)} has been considered in the definition (18) of J multi h (a, b) instead of the effective elasticity tensor C eff (b) itself.The multiscale indicator J multi h (a, b) defined by (18) allows for quantifying the relative distance between (i) the macroscopic elasticity tensor C macro (a) involved in the deterministic boundary value problem ( 2)-( 6) at macroscale, and (ii) the statistical mean value of the effective elasticity tensor C eff (b) calculated by a computational stochastic homogenization method in the mesoscopic subdomain Ω RVE of the random elasticity field C meso (b) involved in the stochastic boundary value problem ( 7)-( 10) at mesoscale.

Comments
It should be noted that in the original formulation initially proposed [140], the numerical indicator J meso (b) was not introduced.The improved formulation proposed in the present work is more advanced than the original formulation initially proposed in Reference [140] to the extent that it involves an additional mesoscopic numerical indicator, namely J meso (b), so that the parameter a and the three components δ, and h of the hyperparameter b each have their own dedicated numerical indicator.Thus, the number of single-objective cost functions being equal to the number of parameters to optimize, it is possible to substitute the computationally expensive global search algorithm used in Reference [140], which belongs to the class of random search, genetic and evolutionary algorithms [146][147][148][149][150][151][152][153][154][155][156], with a more computationally efficient optimization algorithm, such as the fixed-point iterative algorithm considered in the present work (see Section 8).Indeed, even using parallel processing and computing tools, the computational cost incurred by the global optimization algorithm (genetic algorithm) used in Reference [140] remains high due to the large stochastic dimension of the tensor-valued random elasticity field C meso (b), so that the multi-objective optimization problem can be numerically intractable, with the current available computer resources, in very high stochastic dimension for large-scale (non-)linear computational models of three-dimensional random microstructures.The computational cost of the genetic algorithm is compared to the one of the fixed-point iterative algorithm in terms of the number of evaluations of the stochastic computational model in the 2D validation example presented in Section 10.1.It provides a measure of the computational efficiency that is independent of the computer hardware used to perform the numerical simulations.Lastly, it should be noted that an alternative mesoscopic numerical indicator J meso δ (b) is used compared to the previous work in Reference [140] without degrading the performance in terms of accuracy.

Multiscale Statistical Inverse Problem Formulated as a Multi-Objective Optimization Problem
The multiscale statistical inverse identification of parameter a and hyperparameter b can be performed simultaneously by formulating the multiscale statistical inverse problem as a multi-objective optimization problem, that is (a macro , b meso ) = arg min where J (a, b) is the (vector-valued) multi-objective cost function consisting of the four aforementioned numerical indicators as single-objective cost functions and defined for any vector a ∈ A macro and any vector b ∈ B meso by In accordance with the strategy for solving the multiscale statistical inverse problem (see Section 5.2), for a better computational efficiency, the multiscale statistical inverse identification of a and b is performed sequentially by splitting the multi-objective optimization problem into two subproblems solved one after the other: 1. a macroscale inverse problem formulated as a single-objective optimization problem that consists in calculating the optimal value a macro of parameter a in A macro that minimizes the macroscopic numerical indicator J macro (a), that is 2. a mesoscale statistical inverse problem formulated as a multi-objective optimization problem that consists in calculating the optimal value b meso of hyperparameter b in B meso that minimizes the two mesoscopic numerical indicators J meso δ (b) and J meso (b) as well as the multiscale numerical indicator J multi h (a macro , b) simultaneously, that is where J meso (b) is the (vector-valued) multi-objective cost function defined for any vector b ∈ B meso by

Numerical Methods for Solving the Multi-Objective Optimization Problem
The deterministic boundary value problem (2)-( 6) defined on domain Ω macro at macroscale and the stochastic boundary value problem ( 7)-( 10) defined on a subdomain Ω meso ⊂ Ω macro at mesoscale are both discretized using a classical displacement-based finite element method (FEM) [157,158].The mathematical expectations of the quantities of interest of the stochastic boundary value problem ( 7)-(10) involved in the three numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) are estimated using the Monte Carlo numerical simulation method [106][107][108]159,160] with N s independent realizations {C meso (θ r )} 1 r N s of C meso .For the computation of the optimal value a macro , the classical single-objective optimization problem (21) is solved using the Nelder-Mead simplex algorithm [161][162][163][164][165]. For the computation of the optimal value b meso , the non-trivial multi-objective optimization problem (22) does not admit a single global optimal solution, but inherently gives rise to a set of optimal solutions (called Pareto optima) resulting from a trade-off among the three components J meso δ (b), J meso (b) and J multi h (a macro , b) of the multi-objective cost function J meso (b) which are competing and a priori conflicting.Based on the concept of noninferiority [166] (also called Pareto optimality) for characterizing the components of a multi-objective function, a noninferior (or Pareto optimal) solution is such that an improvement in any objective function requires a degradation of some of the other objective functions, whereas an inferior solution is such that an improvement can be attained in all the objective functions.The set of all the noninferior solutions in the parameter space is called the Pareto optimal set and the corresponding objective function values in the multidimensional objective function space is called the Pareto optimal front.The interested reader can refer to References [151][152][153][154][155][156]167] and the references therein for an overview of nonlinear multi-objective optimization methods including the fundamental principles, some Pareto (near-)optimality conditions and a number of traditional and evolutionary optimization algorithms.In Reference [140], the multi-objective optimization problem under consideration has been successfully solved by using the genetic algorithm [151,156] that allows for constructing and finding a set of local Pareto optimal solutions that should be sufficiently representative of the whole Pareto optimal set and as many and diverse as possible for further selection [153,167].The best compromise optimal solution is selected among all the potential Pareto optimal solutions as the one that minimizes the distance to a utopian solution that is constituted by the individual optimal solutions of the conflicting components of the multi-objective function, which corresponds to the origin of the Pareto front.
In the present work, a dedicated numerical indicator has been set up specifically for each component of hyperparameter b = (δ, , h), allowing for the use of a simpler and more efficient multi-objective optimization algorithm, namely a fixed-point iterative algorithm.Starting from an ad hoc initial guess, it consists in sequentially minimizing J meso δ (b), J meso (b) and J multi h (a macro , b) respectively with respect to δ, and h in their sets of admissible values that are such that b = (δ, , h) belongs to B meso .The iterative process is stopped when the residual norm between two iterates becomes lower than a user-specified prescribed tolerance for each of the three single-objective optimization problems.Numerical results have shown that, for the problem under consideration, such a fixed-point iterative algorithm can achieve the same precision as the genetic algorithm in terms of convergence but with a lower overall computational cost (see the numerical examples in Sections 10 and 11).The main drawback of such a numerical optimization algorithm lies in the choice of the initial values used to start the algorithm that may be critical for the localization of the final global convergence region.Besides, note that the fixed-point iterative method introduced in this work could a priori be applied to the original formulation proposed in Reference [140], but it would lead to minimize the objective function J meso δ (b) with respect to δ and simultaneously given the other hyperparameters h.Although it is possible, the problem is that J meso δ (b) is very sensitive to δ but less sensitive with respect to , since it has been tailored to perform the identification of the optimal value δ meso of δ and not the one meso of .Consequently, using such a fixed-point iterative strategy would yield uncertainties on the identified value meso of .It is the reason why the additional objective function J meso (b) has been introduced and for which the sensitivity is of first order with respect to and of second order with respect to δ.

Probabilistic Model for a Robust Identification of the Hyperparameters
When several non-overlapping mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q are available for experimental measurements for the same test specimen instead of a unique observation domain Ω meso exp , then the solution of the multi-objective optimization problem presented in Section 7 can yield different optimal values b meso 1 , . . ., b meso Q of hyperparameter b meso when experimental data comes from one mesoscopic domain of observation to another since mesoscopic indicators J meso δ (b) and J meso (b) depend on the values of experimental displacement fields u meso exp,1 , . . ., u meso exp,Q that are measured on each of them.Consequently, the optimal value b meso of hyperparameter b should be considered as uncertain and should be modeled as a vector-valued random variable B = (D, L, H) for which b meso 1 , . . ., b meso Q are assumed to be Q independent realizations.Thus, in Reference [140], a robust identification of the optimal value b opt is proposed by averaging the identified values b meso 1 , . . ., b meso Q .Nevertheless, in the present work, an improved strategy is proposed that consists in constructing a prior stochastic model of the vector-valued hyperparameter B by using the MaxEnt principle [68,72,73,77]  (n + 1)/(n + 5) = √ 7/11 ≈ 0.7977 < 1 (with n = 6 in linear elasticity), (iii) the random components of random vector L are (statistically independent) positive-valued random variables a.s.for which the mean value is given in ]0 , +∞[ and the values are unlikely near zero by construction of the mesoscale modeling, otherwise it would mean the current scale of the computational model is not correct and too large, (iv) the random components of random vector H take their values a.s. in the admissible set H meso .We then have for all b = (δ, , h) ∈ B meso , where p D (δ) = 1 ) are positive parameters to be identified.We refer the reader to Reference [168] for a detailed construction of the prior stochastic model of H and a rigorous characterization of the statistical dependence between the components of random elasticity tensors exhibiting a.s.some given material symmetry properties for the six highest levels of linear elastic symmetries.For the special case of isotropic materials, we have H meso = ]0 , +∞[×]0 , +∞[ and the prior probability density function p H of random vector H is written as for all h = (h 1 , h 2 ) ∈ H meso , in which where The mean values of H 1 and H 2 are respectively equal to (1 − λ)/λ 1 and (1 − 5λ)/λ 2 , and the dispersion coefficients of H 1 and H 2 are respectively equal to 1/ √ 1 − λ and 1/ √ 1 − 5λ.Note that the probability density functions of H 1 and H 2 both involve the same hyperparameter λ < 1/5 that controls the level of statistical fluctuations of both H 1 and H 2 .In addition, H 1 and H 2 cannot be deterministic variables, since their dispersion coefficients are non zero whatever the value of λ < 1/5.Finally, the probabilistic model of B = (D, L, H) involves the unknown vector-valued hyperparameter s = (s 1 , s 2 , The optimal value s opt of s is determined using the MLE method with the available data that are the Q independent realizations b meso 1 , . . ., b meso Q of random vector B. The MLE method consists in computing s opt by solving the following optimization problem where s → L(s; b meso The accuracy of the identified optimal value s opt is then all the higher as the number Q of mesoscopic domains of observation is large but at the expense of a higher computational cost.Lastly, the optimal value b opt of vector-valued hyperparameter b ∈ B meso is computed by solving the following optimization problem b opt = arg max b∈B meso p B (b; s opt ).
Hence, optimal value b opt corresponds to the most probable value of random vector B according to the identified probability distribution represented by its probability density function p B (•; s opt ) parameterized by s opt .Note that the averaging approach presented in Reference [140] is a particular case of the MLE method presented in this section if the prior stochastic models of D, L and H are uniform random variables.It is the reason why a better robust identification is expected since the prior stochastic model of B has been improved in this work.In the present work, since D is modeled as a uniform random variable on ]0 , δ sup [, the optimal value δ opt of δ is simply obtained by averaging the Q independent realizations δ meso 1 , . . ., δ meso Q of D. A more advanced prior stochastic model for D could have been considered, for instance by adding as available information that its mean value is given and its values are unlikely near zero, thus leading to a unimodal probability density function p D with support ]0 , δ sup [ and with a higher parameterization than the simple uniform probability density function considered here.

Numerical Validation of the Multiscale Identification Method on In Silico Materials in 2D Plane Stress and 3D Linear Elasticity
In this section, we present a numerical application of the improved multiscale identification methodology proposed in the present work within the framework of 2D plane stress and 3D linear elasticity theories by using in silico materials for which the macroscopic and mesoscopic mechanical properties are known.The required multiscale "experimental" kinematic fields have been obtained through numerical simulations using one random realization of the random elasticity field in SFE + (see Section 4) not restricted from R 3 to some mesoscopic domain Ω meso but restricted to the whole macroscopic domain Ω macro for a given experimental value b meso exp of hyperparameter b ∈ B meso .The solution of a deterministic boundary value problem over this macroscopic domain Ω macro is then computed for a heterogeneous random elasticity field whose spatial correlation lengths correspond to the characteristic sizes of the heterogeneities at microscale.This deterministic boundary value problem is solved using a classical numerical method (FEM) whose computational cost is high and potentially prohibitive in 3D, what can be avoided by computational homogenization methods, but it is required to completely simulate the multiscale "experimental" measurements.

Validation on an In Silico Specimen in Compression Test in 2D Plane Stress Linear Elasticity
For this first numerical validation example, a 2D plane stress assumption is considered.Macroscopic domain of observation Ω macro exp is a 2D square domain and it exactly corresponds to the cross-section of macroscopic domain Ω macro and such that Ω macro obs = Ω macro exp since the test specimen is in silico.The dimensions of 2D macroscopic domain of observation Ω macro exp are 1×1 cm 2 in a fixed Cartesian frame (O, x 1 , x 2 ) of R 2 .It is possible to introduce a set of Q = 16 non-overlapping 2D square mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q ⊂ Ω macro exp for which the mesoscale dimensions are 1×1 mm 2 (see Figure 1 for a schematic representation of domains of observation Ω macro exp and Ω meso exp,1 , . . ., Ω meso exp,Q ).Consequently, observation domain Ω meso obs , for which the dimensions are also 1×1 mm 2 , is defined as the 2D square cross-section of mesoscopic domain Ω meso .Deterministic surface force field f macro is uniformly distributed on the top boundary of macroscopic domain Ω macro and applied along the (downward vertical) −x 2 direction with an intensity of 5 kN such that f macro = 5 kN/cm 2 = 5×10 7 N/m 2 , while the bottom boundary of macroscopic domain Ω macro is clamped.

Parameterization of the Macroscopic and Mesoscopic Models
At macroscale, the solution of deterministic boundary value problem ( 2)-( 6) with 2D plane stress assumption depends only on 6 components {S macro (a)} ijkh of deterministic compliance tensor S macro (a) with i, j, k, h ∈ {1, 2}.Consequently, the solution at macroscale depends only on the components of a 2D fourth-order compliance tensor S macro 2D (a) that is defined by {S macro 2D (a)} ijkh = {S macro (a)} ijkh for all i, j, k, h ∈ {1, 2}.Then, a 2D fourth-order elasticity tensor at macroscale can be introduced and defined by C macro 2D (a) = (S macro 2D (a)) −1 .Since within the framework of linear elasticity theory, any isotropic material is completely characterized by a bulk modulus κ and a shear modulus µ at macroscale, then we have the vector-valued parameter a = (κ, µ).In particular, we have chosen the experimental value  7)-( 10) with 2D plane stress assumption depends only on 6 components {S meso (b)} ijkh of random compliance tensor field S meso (b) with i, j, k, h ∈ {1, 2} or equivalently on every 21 components {C meso (b)} ijkh of random elasticity tensor field C meso (b) with i, j, k, h ∈ {1, 2, 3}.It is the reason why we have chosen to construct the prior stochastic model of the random compliance tensor field S meso (b) as presented in Section 4 and the stochastic boundary value problem ( 7)-( 10) is solved in using (11) rather than (8).Furthermore, its mean function S meso is spatially constant and models an isotropic elastic medium that is completely characterized by a mean bulk modulus κ and a mean shear modulus µ at mesoscale.Consequently, the vector-valued hyperparameter b = (δ, , h) involves only (i) a dispersion parameter δ, (ii) a spatial correlation length that is such that 1 = 2 = in order to be consistent with the effective model at macroscale for which the material is assumed to be isotropic and with 3 = +∞ in order to be consistent with the 2D plane stress assumption, and (iii) a vector-valued hyperparameter h = (κ, µ) gathering the mean bulk modulus κ and the mean shear modulus µ at mesoscale.In particular, we have chosen the experimental value b meso exp = (δ meso exp , meso exp , κ meso exp , µ meso exp ) with δ meso exp = 0.40, meso exp = 125 µm, κ meso exp = 13.75GPa and µ meso exp = 3.587 GPa, corresponding to a mean Young's modulus E meso exp = 9.900 GPa and a mean Poisson's ratio ν meso exp = 0.380 GPa.For identification purposes and further computational savings, we consider a reduced admissible set B meso ad ⊂ B meso for the vector-valued hyperparameter b = (δ, , κ, µ) such that . ., Ω meso exp,Q are then searched on this multidimensional grid of n V ×n V ×n V ×n V points in the hypercube B meso ad .Within the framework of linear elasticity under 2D plane stress assumption, both the deterministic boundary value problem ( 2)-( 6) and the stochastic boundary value problem ( 7)-( 10) are solved by discretizing the 2D macroscopic and mesoscopic domains of observation Ω macro obs and Ω meso obs in space using the FEM.The finite element meshes of 2D square domains Ω macro obs and Ω meso obs are structured meshes made up with 4-nodes linear quadrangular elements with Gauss-Legendre quadrature rule.The stochastic boundary value problem ( 7)-( 10) at mesoscale is solved using the Monte Carlo numerical method.Mesh convergence analyses of the numerical solutions of the deterministic boundary value problem ( 2)-( 6) at macroscale and of the stochastic boundary value problem ( 7)- (10) at mesoscale have been performed in order to define accurate finite element approximations at both macroscopic and mesoscopic scales.The finite element mesh of 2D macroscopic domain Ω macro obs is a regular grid containing 25×25 quadrangular elements with uniform element size h macro = 0.4 mm = 4×10 −4 m in each spatial direction.It thus comprises 676 nodes and 625 elements, with 1300 unknown degrees of freedom (dofs).The finite element mesh of 2D mesoscopic domain Ω meso obs is a regular grid containing 100×100 quadrangular elements with uniform element size h meso = 10 µm = 10 −5 m in each spatial direction.It thus comprises 10, 201 nodes and 10, 000 elements, with 20, 000 unknown dofs.The number of Gauss integration points per spatial correlation length used for numerical quadrature over 2D macroscopic domain of observation Ω macro obs and 2D mesoscopic domain of observation Ω meso obs is n G = 4 in each spatial direction.Concerning the computational stochastic homogenization with 2D plane stress assumption, we consider a 2D square domain Ω RVE of side length B RVE defined in a Cartesian frame (O, x 1 , x 2 ) and we use the homogenization method with static uniform boundary conditions (i.e., with homogeneous stresses) which is appropriate for linear elasticity under 2D plane stress assumption.Note that only the components {S eff (b)} ijkh with i, j, k, h ∈ {1, 2} can be calculated.We then obtain a 2D fourth-order effective compliance tensor As the mathematical expectations involved in each of the numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) are evaluated using the Monte Carlo numerical method, statistical convergence analyses of their statistical estimators with respect to the number of independent realizations N s have been carried out and a convergence has been reached for N s = 500.Sensitivity analyses of each of the three numerical indicators have been performed with respect to each of the hyperparameters δ, , h = (κ, µ), respectively, in the reduced admissible set B meso ad = [0.25 , 0.50]× [20 , 250] µm× [8.5 , 17] GPa× [2.15 , 4.50] GPa.Hence, it can be shown that each numerical indicator is sufficiently sensitive to the variation of its dedicated hyperparameter and that the multi-objective optimization problem (22) to be solved is well-posed.
Recall the multiscale statistical inverse problem has been formulated into two decoupled optimization problems in a and b, respectively, to be solved sequentially (see Section 7), namely (i) a macroscale single-objective optimization problem (21) for the inverse identification of the optimal value a macro of parameter a in its admissible set A macro , and (ii) a mesoscale multi-objective optimization problem (22) for the statistical inverse identification of the global optimal value b opt of hyperparameter b in its reduced admissible set B meso ad .

Resolution of the Single-Objective Optimization Problem at Macroscale
In this paragraph, we present the results of the first single-objective optimization problem (21) at macroscale which consists in minimizing the macroscopic numerical indicator J macro (a) constructed in the macroscopic domain of observation Ω macro exp for identifying the optimal value a macro of a at macroscale.The single-objective optimization problem (21) at macroscale has been solved using the Nelder-Mead simplex algorithm.The identification results are reported in Table 1 and show that the relative error between the identified optimal value a macro = (13.901,3.685) in [GPa] and the reference experimental value a macro exp = (14.328,3.670) in [GPa] used for the construction of the numerically simulated "experimental" database remains small (less than 3% and 0.5% for κ macro and µ macro , respectively), allowing to validate the proposed identification methodology in 2D plane stress linear elasticity for the resolution of the single-objective optimization problem (21) at macroscale.In this paragraph, we present the results of the second multi-objective optimization problem ( 22) at mesoscale which consists in simultaneously minimizing the three numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) constructed in each of the Q = 16 mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q using the optimal parameter a macro = (13.901,3.685) in [GPa] previously identified at macroscale (see the previous paragraph) for identifying the global optimal value b opt of b at mesoscale.The multi-objective optimization problem (22) at mesoscale has been solved using the fixed-point iterative algorithm on the one hand and the genetic algorithm on the other hand for comparison purposes.In order to analyze the numerical efficiency of these two resolution approaches, instead of evaluating the computing time which strongly depends on the computer hardware used, we choose in this work to compare the number of evaluations of the random solution of the stochastic boundary value problem ( 7)-( 10) at mesoscale (i.e., the number of calls to the deterministic numerical model at mesoscale) required by each algorithm to achieve the desired convergence.
The identification results obtained with the fixed-point iterative algorithm are summarized in Table 2 for the set of Q = 16 mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q , namely the set of Q identified values b meso 1 , . . ., b meso Q and numbers of iterations n 1 , . . ., n Q required to reach the desired convergence, with a convergence criterion on the residual norm between two iterations that must be less than a prescribed tolerance set to 10 −9 , and the global optimal value b opt computed by using the MLE method.On the one hand, there are greater variations between the identified values meso 1 , . . ., meso Q and δ meso 1 , . . ., δ meso Q , reflecting the fact that the two associated mesoscopic numerical indicators J meso δ (b) and J meso (b) depend directly on the experimental field measurements on each mesoscopic domain of observation.On the other hand, the lower variability between the identified values κ meso 1 , . . ., κ meso Q and µ meso 1 , . . ., µ meso Q can be explained by the fact that the associated multiscale numerical indicator J multi h (a macro , b) does not depend directly on the experimental field measurements on each mesoscopic domain of observation but is rather conditioned by the identified values meso 1 , . . ., meso Q and δ meso 1 , . . ., δ meso Q .Thus, the relative errors calculated on these two hyperparameters are essentially due to the quality of the discretization of the reduced admissible set B meso ad .In particular, the fixed-point iterative algorithm has selected the same identified value µ meso 1 717 GPa (among the n V = 10 test points in [2.15 , 4.50] GPa) for the Q = 16 mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q .Clearly, a finer grid (with n V > 10) might yield different values for the identified hyperparameter µ meso selected by the optimization algorithm.It is the reason why a prior probabilistic model for the identified hyperparameters has been introduced.The number of evaluations of the stochastic computational model needed by the fixed-point iterative algorithm is given by n FP tot = 3 n V N s ∑ 16 q=1 n q , where the superscript FP refers to "Fixed-Point" and n V is the number of evaluations of a numerical indicator to search for the minimum with respect to the associated hyperparameter.Figure 3 shows the probability density functions p D , p L , p K and p M of random variables D, L, K and M, respectively, which are defined in Section 9 with the two components H 1 = K and H 2 = M of random vector H = (K, M).As suggested by the identification results shown in Table 2 and as already mentioned in Section 9, a more advanced prior stochastic model for D would have been preferable to obtain a unimodal probability density function p D with support ]0 , δ sup [ and which would be concentrated around the reference experimental value δ meso exp = 0.4.Besides, although all the independent realizations µ meso 1 , . . ., µ meso Q of M given in Table 2 are equal to the same identified value 3.717 GPa, the probability density function p M does not correspond to the Dirac measure on R at point 3.717 GPa but to a gamma distribution with a very small dispersion around this value, since for the prior probabilistic model of H = (K, M) considered here, K and M cannot be deterministic variables (see Section 9).We finally obtain the global optimal value b opt = (0.391, 135.328, 12.273, 3.717) ) with relative errors less than 3%, 9%, 11% and 4% for δ opt , opt , κ opt and µ opt , respectively, with respect to the reference experimental value b meso exp = (0.40, 125, 13.75, 3.587) in ([−], [µm], [GPa], [GPa]) used to construct the numerically simulated "experimental" database, allowing to validate the proposed identification methodology in 2D plane stress linear elasticity for the resolution of the multi-objective optimization problem ( 22) at mesoscale.Figure 4 shows the evolution of the global optimal values δ opt , opt , κ opt µ opt estimated by the MLE method as a function of the number Q of independent realizations b meso 1 , . . ., b meso Q of random vector B = (D, L, K, M).Although the number Q remains low (less than or equal to 16), we observe that each of the global optimal values tends to converge towards an objective value when Q increases, which demonstrates that the use of the MLE method with the prior probabilistic model of B proposed in this work allows a robust identification of the vector-valued hyperparameter b = (δ, , κ, µ).
In terms of computational efficiency, we can see in Table 2 that the numbers of iterations n 1 , . . ., n Q required to achieve the desired convergence are relatively low (less than or equal to 4) for each of the Q = 16 mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q , leading to a number of calls to the deterministic numerical model at mesoscale of 855, 000.Table 3 contains the global optimal values b opt and the corresponding relative errors (with respect to the reference experimental value b meso exp ) obtained for different values N s ∈ {5, 50, 500} of the number of independent realizations generated for the statistical estimation of the mathematical expectations involved in the different numerical indicators.It can be seen that a strong decrease in the value of N s allows a considerable gain in computing time while maintaining similar results for the identified global optimal values, which can be explained by the use of the MLE method which makes the resolution of the statistical inverse identification problem more robust with respect to the convergence of the statistical estimators used in the numerical indicators of the multi-objective optimization problem (22).The identification results obtained with the genetic algorithm are summarized in Table 4 for the set of Q = 16 mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q , namely the set of Q identified values b meso 1 , . . ., b meso Q and numbers of generations n 1 , . . ., n Q required to reach the desired convergence, and the global optimal value b opt computed by using the MLE method.The initial population used to start the genetic algorithm contains n I = 40 individuals.Figure 5 shows an example of different 2D cross-sections of the Pareto front for the first mesoscopic domain of observation Ω meso exp,1 .The best comprise optimal solution corresponds to the point marked with a green circle on the different 2D cross-sections of the Pareto front, because according to the explanations given in Section 8, it is chosen among all the noninferior (Pareto optimal) solutions generated and selected in the optimal Pareto set (represented by red stars in Figure 5) as the one that minimizes the distance at the origin of the Pareto front in the multidimensional space (of dimension 3) of the multi-objective cost function J meso (b).For reasons of limitation in terms of calculation cost, the number N s of independent realizations used for the statistical estimation of the mathematical expectations involved in the numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) is reduced to N s = 50.Although the statistical convergence of the three numerical indicators is not achieved, the results of Table 3 show that, thanks to the probabilistic modeling of the hyperparameters and the maximum likelihood estimation, the results of the statistical inverse identification method are not significantly affected by a decrease in the value of N s and are therefore robust with respect to the statistical fluctuations of the different numerical indicators.The number of evaluations of the stochastic computational model needed by the genetic algorithm is given by n GA tot = 3 n I N s ∑ 16 q=1 n q , where the superscript GA refers to "Genetic Algorithm". Figure 6 shows the probability density functions p D , p L , p K and p M of random variables D, L, K and M, respectively.We finally deduce the global optimal value b opt = (0.372, 128.401, 11.656, 3.306) in ([−], [µm], [GPa], [GPa]) with relative errors less than 8%, 3%, 16% and 8% for δ opt , opt , κ opt and µ opt , respectively, with respect to the reference experimental value b meso exp = (0.40, 125, 13.75, 3.587) in ), which are acceptable (reasonably good) and similar to the errors obtained by the fixed-point iterative algorithm.There are still some fluctuations in the values κ meso 1 , . . ., κ meso Q and µ meso 1 , . . ., µ meso Q identified on each of the Q = 16 mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q , which was not the case for the fixed-point iterative algorithm.This underlies the numerical resolution of the Pareto front, which depends on the number n V of values in each dimension of the parameter search space.In terms of computational efficiency, we can see that the number n GA tot = 19, 176, 000 of evaluations of the stochastic computational model (resulting from the number of individuals n I = 40 in the initial population and the number of population generations n 1 , . . ., n Q ) is much higher than that n FP tot = 87, 000 required by the fixed-point iterative algorithm with N s = 50 (see Table 3).Finally, the fixed-point iterative algorithm allows significant computational savings (in terms of computational cost) compared to the genetic algorithm.

Validation on an In Silico Specimen in Compression Test in 3D Linear Elasticity
In this section, we present a second validation example in 3D linear elasticity.We assume there are Q = 3 test specimens on which are applied exactly the same external loads at macroscale.Recall that for the validation, the "experimental" tests are actually performed in silico.Macroscopic domain of observation Ω macro exp is exactly the same 3D cubic domain for each test specimen and corresponds to 3D experimental field measurements on the full volume of each test specimen.As for the previous 2D validation example, since the experimental field measurements are actually performed in silico, we also have Ω macro obs = Ω macro exp .The dimensions of each 3D macroscopic domain of observation Ω macro exp are 2×2×2 mm 3 .For each test specimen, the mesoscale dimensions of 3D mesoscopic domain of observation Ω meso exp are 0.5×0.5×0.5 mm 3 (see Figure 7).Deterministic surface force field f macro is uniformly distributed on the top boundary of macroscopic domain Ω macro exp and applied along the (downward vertical) −x 3 direction with an intensity of 2 kN such that f macro = 50 kN/cm 2 = 5×10 8 N/m 2 , while the bottom boundary of macroscopic domain = Ω macro (in green) which contains a 3D cubic mesoscopic domain of observation Ω meso exp = Ω meso (in red) for the numerical validation in 3D linear elasticity.

Parameterization of the Macroscopic and Mesoscopic Models
Within the framework of linear elasticity theory, any material that is isotropic at macroscale can be completely characterized by a bulk modulus κ and a shear modulus µ.Consequently, we have chosen the parameterization a = (κ, µ).In particular, we have chosen the experimental value a macro At mesoscale, we have chosen to construct the prior stochastic model of the random elasticity tensor field C meso as presented in Section 4 and the stochastic boundary value problem ( 7)-( 10) is solved in using (8) rather than (11).Furthermore, its mean function C meso is spatially constant and models an isotropic elastic medium that is completely characterized by a mean bulk modulus κ and a mean shear modulus µ at mesoscale.Consequently, the vector-valued hyperparameter b = (δ, , h) involves only (i) a dispersion parameter δ, (ii) a spatial correlation length that is such that 1 = 2 = 3 = in order to be consistent with the effective model at macroscale for which the material is assumed to be isotropic, and (iii) a vector-valued hyperparameter h = (κ, µ) gathering the mean bulk modulus κ and the mean shear modulus µ at mesoscale.In particular, we have chosen the experimental value b meso exp = (δ meso exp , meso exp , κ meso exp , µ meso exp ) with δ meso exp = 0.32, meso exp = 80 µm, κ meso exp = 145 GPa and µ meso exp = 67.3GPa, corresponding to a mean Young's modulus E meso exp = 174.85GPa and a mean Poisson's ratio ν meso exp = 0.2990 GPa.As already mentioned in Section 10. of hyperparameters b for each of the 3 in silico test specimens are then searched on this multidimensional grid of n V ×n V ×n V ×n V points in the hypercube B meso ad .The classical displacement-based FEM is used for the spatial discretization of (i) the deterministic boundary value problems defined by ( 2)-( 6) in replacing C macro by Q = 3 independent realizations of the random apparent elasticity tensor field C meso on Ω macro instead of Ω meso to simulate both the "experimental" macroscopic displacement field u macro exp in Ω macro exp = Ω macro and the mesoscopic displacement field u meso exp in Ω meso exp = Ω meso , (ii) the deterministic boundary value problem defined by ( 2)-( 6) to calculate the macroscopic displacement field u macro in domain Ω macro obs = Ω macro , and (iii) the stochastic boundary value problems defined by ( 7)- (10) to calculate the random mesoscopic displacement fields U meso in using experimental data obtained by solving (i) that are the experimental displacement fields u meso exp measured on the boundary of domain Ω meso exp = Ω meso obs for each realization of C meso .The stochastic solver used for solving the stochastic boundary value problem ( 7)- (10) at mesoscale is the Monte Carlo numerical method.As 3D macroscopic and mesoscopic domains Ω macro and Ω meso are cubic domains, we consider for each of them a spatial discretization with a structured mesh made up with 8-nodes linear hexahedral elements with Gauss-Legendre quadrature rule.The finite element mesh of 3D macroscopic domain Ω macro is made with the same spatial discretization as the one used for the 2D validation example at macroscale, that is a structured mesh of 25×25×25 = 15, 625 hexahedral elements with uniform element size h macro = 80 µm = 8×10 −5 m in each spatial direction.The finite element mesh of 3D mesoscopic domain Ω meso is made with the same spatial discretization as the one used for the 2D validation example at mesoscale and whose element size depends on the smallest spatial correlation length considered, that is a structured mesh of 20×20×20 = 8000 hexahedral elements with uniform element size h meso = /(n G /2) = (50 µm)/2 = 25 µm = 2.5×10 −5 m in each spatial direction, with n G = 4 Gauss integration points per spatial correlation length.
Concerning the computational stochastic homogenization, as for the 2D validation example, the size B RVE of representative volume element Ω RVE is defined as a function of the spatial correlation length such that B RVE = 20× = 20×50 µm = 1 mm = 10 −3 m.
Recall that the multiscale statistical inverse problem has been formulated into two decoupled optimization problems in a and b, respectively, to be solved sequentially (see Section 7), namely (i) a macroscale single-objective optimization problem (21) for the inverse identification of the optimal value a macro of parameter a in its admissible set A macro , and (ii) a mesoscale multi-objective optimization problem (22) for the statistical inverse identification of the global optimal value b opt of hyperparameter b in its reduced admissible set B meso ad .

Resolution of the Single-Objective Optimization Problem at Macroscale
In this paragraph, we present the results of the first single-objective optimization problem (21) at macroscale which consists in minimizing the macroscopic numerical indicator J macro (a) constructed in each of the Q = 3 in silico test specimens for identifying the optimal value a macro of a at macroscale.The single-objective optimization problem (21) at macroscale has been solved using the Nelder-Mead simplex algorithm.The identification results are reported in Table 5 and show that the relative error between the identified optimal value a macro = (138.783,64.355) in [GPa] and the reference experimental value a macro exp = (138.758,64.377) in [GPa] used for the construction of the numerically simulated "experimental" database remains very small (less than 0.02% and 0.04% for κ macro and µ macro , respectively), allowing to validate the proposed identification methodology in 3D linear elasticity for the resolution of the single-objective optimization problem (21) at macroscale.In this paragraph, we present the results of the second multi-objective optimization problem (22) at mesoscale which consists in simultaneously minimizing the three numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) constructed in each of the Q = 3 in silico tests specimens using the optimal parameter a macro = (138.783, 64.355) in [GPa] previously identified at macroscale (see the previous paragraph) for identifying the global optimal value b opt of b at mesoscale.
In contrast, unlike the 2D validation example, the multi-objective optimization problem (22) has been solved only with the fixed-point iterative algorithm using the same convergence criterion on the residual norm between two iterations that must be less than a tolerance set to 10 −9 and by searching for the solution of the multi-objective optimization problem (22) in a multidimensional grid of n V ×n V ×n V ×n V points in the reduced admissible set B meso ad ⊂ R 4 .The genetic algorithm has not been used because the resulting computational cost was too high with the available computational resources.The number of independent realizations for the statistical estimation of the mathematical expectations involved in the different numerical indicators is set to N s = 500.The number of evaluations of the stochastic computational model needed by the fixed-point iterative algorithm is given by n FP tot = 3 n V N s ∑ 3 q=1 n q .Table 6 reports the identification results obtained with the fixed-point iterative algorithm for the set of Q = 3 in silico tests specimens, namely the set of identified values b meso 1 , b meso 2 , b meso 3 and numbers of iterations n 1 , n 2 , n 3 required to reach the desired convergence (with a tolerance set to 10 −9 ), and the global optimal value b opt computed by using the MLE method.As for the 2D validation example, there are more significant variations between the identified values meso being almost identical for each in silico test specimen, we directly identify the global optimal values κ opt and µ opt without using the MLE method for the random variables K and M. Figure 8 shows the probability density functions p D and p L defined in Section 9 and associated to random variables D and L, respectively.We finally obtain the global optimal value b opt = (0.330, 91.236, 150.000, 64.722) in ([−], [µm], [GPa], [GPa]) with relative errors less than 4% for δ opt , opt , κ opt and µ opt with respect to the reference experimental value b meso exp = (0.32, 80, 145, 67.3) in ([−], [µm], [GPa], [GPa]) used to construct the numerically simulated "experimental" database, allowing to validate the proposed identification methodology in 3D linear elasticity for the resolution of the multi-objective optimization problem (22) at mesoscale.In terms of computational efficiency, we can see in Table 6 that the numbers of iterations n 1 , n 2 , n 3 required to achieve the desired convergence are relatively low (less than or equal to 4) for each of the 3 in silico test specimens yielding a number of calls to the deterministic numerical model at mesoscale of 150, 000.
Finally, the results obtained for the identification of the parameters of the deterministic model at macroscale and of the hyperparameters of the prior stochastic model at mesoscale for both validation examples in 2D plane stress and 3D linear elasticity, for which the reference experimental values are known a priori, demonstrate the efficiency, accuracy and robustness of the improved identification methodology, thereby allowing to apply it in the next section to a real biological material (beef femur cortical bone) with real experimental field measurements.Lastly, let us mention that the fixed-point iterative algorithm introduced in the present work to solve the multi-objective optimization problem allows a considerable gain in terms of computational cost compared to the genetic algorithm used in Reference [140].

Numerical Application of the Multiscale Identification Method to Real Beef Cortical Bone in Plane Stress Linear Elasticity
In this section, we present a numerical application of the proposed multiscale identification methodology within the framework of 3D linear elasticity with 2D plane stress assumption by using a real experimental database made up of 2D multiscale optical measurements of displacement fields (obtained by DIC method) for only a single test specimen of cortical bone coming from a beef femur.The multiscale experimental test configuration corresponds to the one described in Section 3 and already considered in the 2D and 3D numerical validation examples presented in Section 10.Technical details concerning the experimental protocol (specimen preparation, measuring bench, optical image acquisition system and DIC method) for obtaining the multiscale field measurements (performed simultaneously at both macroscale and mesoscale) can be found in Reference [145].The unique test specimen at macroscale is a cubic shaped sample with dimensions 1×1×1 cm 3 prepared from bovine cortical bone.Even though such a biological tissue is often considered and modeled as a deterministic homogeneous medium with a transversely isotropic linear elastic behavior at macroscale ( 10 mm), its microstructure at mesoscale (from 500 µm to 5 mm) contains randomly arranged osteons with some resorption cavities (lacuna), that are the principal types of inclusions/inhomogeneities, embedded in a matrix constituted by circumferential interstitial lamella surrounding Haversian canals.As a consequence, it is an anisotropic (heterogeneous) composite material with a complex hierarchical structure, which can be considered and modeled as a random linear elastic medium at mesoscale, and is therefore well adapted to the experimental application of the multiscale identification methodology developed in the present work.The single specimen is clamped on its lower face and loaded under vertical uniaxial compression onto its upper face with a maximal resultant force of 9 kN so as to preserve a linear elastic material behavior.In order to reduce the measurement noises (induced by the speckled pattern technique, the lighting of the observed 2D face, the optical image acquisition system, etc.), a Gaussian spatial filter classically used in image processing has been applied to smooth the experimental displacement fields u macro exp = (u macro exp,1 , u macro exp,2 ) and u meso exp = (u meso exp,1 , u meso exp,2 ) measured at macroscale and at mesoscale, respectively.The images of experimental displacement fields at macroscale and at mesoscale have been filtered with a 2D Gaussian smoothing kernel with standard deviation 3.5.This value has been chosen as a qualitative compromise allowing to regularize/smooth the experimental kinematic fields without removing the spatial fluctuations that are of the same order of magnitude as the lower bound of the search interval for the spatial correlation length .Such a spatial filter is also necessary to prevent the optimization algorithms from converging to a zero spatial correlation length.Figures 9 and 10 represent the two components of macroscopic experimental displacement field u macro exp over the 2D macroscopic domain Ω macro exp and the ones of mesoscopic experimental displacement field u meso exp over the 2D mesoscopic domain Ω meso exp , respectively, before and after application of the Gaussian spatial filter.

Parameterization of the Macroscopic and Mesoscopic Models
In accordance with the experimental configuration and associated multiscale measurements, 2D plane stresses are assumed and consequently, the deterministic and stochastic computational models at macroscale and mesoscale are the same as those used for the 2D validation example presented in Section 10.1.Thus, the modeling at macroscale and at mesoscale for the prior stochastic model, the hyperparameters and the parameterization are also exactly the same as in Section 10.1, namely defining S meso in SFE + and introducing vector-valued parameter a = (κ, µ) and vector-valued hyperparameter b = (δ, , κ, µ).Optimal values of the latter are assumed to be in the reduced admissible set B meso ad ⊂ B meso constructed from information available in the literature such that δ ∈ [0.30 , 0.65], ∈ [50 , 100] µm, κ ∈ [9.5 , 11] GPa and µ ∈ [3.5 , 5.0] GPa, instead of the full admissible space B meso = ]0 , δ sup [×]0 , +∞[× ]0 , +∞[ 2 with δ sup = (n + 1)/(n + 5) = √ 7/11 ≈ 0.7977 < 1 (with n = 6 in linear elasticity).As in the 2D validation example, this reduced admissible set B meso ad is discretized into n V = 10 points evenly spaced in each dimension for which the three numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) defined in Section 6.4 are evaluated and compared.
As for Section 10.1, both the deterministic boundary value problem (2)-( 6) set on the macroscopic domain Ω macro and the stochastic boundary value problem ( 7)- (10) set on the mesoscopic domain Ω meso with 2D plane stress assumption are solved by discretizing the 2D domains of observation Ω macro obs and Ω meso obs in space using the FEM.As 2D macroscopic and mesoscopic domains of observation Ω macro obs and Ω meso obs are square domains, we consider for both a spatial discretization with a structured mesh made up with 4-nodes linear quadrangular elements with Gauss-Legendre quadrature rule, in order to be consistent with the regular grids used for the acquisition and discretization of experimental data.The 2D macroscopic domain Ω macro obs with macroscale dimensions 1×1 cm 2 is discretized with a structured mesh of 9×9 = 81 quadrangular elements with uniform element size h macro = 1.111 mm = 1.111×10 −3 m in each spatial direction.The 2D mesoscopic domain Ω meso obs with mesoscale dimensions 1×1 mm 2 is discretized with a structured mesh of 99×99 = 9801 quadrangular elements with uniform element size h meso = 10.10 µm = 1.010×10 −5 m in each spatial direction.As for the 2D validation example, the size B RVE of representative volume element Ω RVE is defined with respect to the spatial correlation length such that B RVE = 20× .The stochastic boundary value problem ( 7)- (10) at mesoscale is solved using the Monte Carlo numerical method and statistical convergence analyses have been systematically performed to set the number of independent realizations for the statistical estimation of the mathematical expectations involved in the different numerical indicators to the value N s = 500.In this paragraph, we present the results of the first single-objective optimization problem (21) at macroscale which consists in minimizing the macroscopic numerical indicator J macro (a) constructed in the macroscopic domain of observation Ω macro obs for identifying the optimal value a macro of a at macroscale.The single-objective optimization problem (21) at macroscale has been solved using the Nelder-Mead simplex algorithm.Table 7 gives the identified optimal value a macro = (11.335,4.781) in [GPa], corresponding to a macroscopic transverse bulk modulus κ macro = 11.335GPa and a macroscopic transverse shear modulus µ macro = 4.781 GPa, or equivalently to a macroscopic transverse Young's modulus E macro = 12.575 GPa and a macroscopic transverse Poisson's ratio ν macro = 0.3151, which are in coherence with the values already published and available in the literature for this type of biological material.In this paragraph, we present the results of the second multi-objective optimization problem (22) at mesoscale which consists in simultaneously minimizing the three numerical indicators J meso δ (b), J meso (b) and J multi h (a macro , b) constructed in the mesoscopic domain of observation Ω meso obs using the optimal parameter a macro = (11.335,4.781) in [GPa] previously identified at macroscale (see the last paragraph) for identifying the global optimal value b opt of b at mesoscale.The multi-objective optimization problem (22) has been solved only by using the fixed-point iterative algorithm (with a convergence criterion on the residual norm between two iterations that must be less than a tolerance set to 10 −9 ) and by searching for the solution of the multi-objective optimization problem (22) in a multidimensional grid of n V ×n V ×n V ×n V points in the reduced admissible set B meso ad ⊂ R 4 .The number of evaluations of the stochastic computational model needed by the fixed-point iterative algorithm is given by n FP tot = 3 n V N s n FP , where n FP is the number of iterations required to reach the desired convergence for the considered mesoscopic domain of observation Ω meso obs .Table 8 gives the identified optimal value b meso = (0.533, 61.111, 10.500, 4.667) in ([−], [µm], [GPa], [GPa]) obtained with the fixed-point iterative algorithm, corresponding to a dispersion parameter δ meso = 0.533, a spatial correlation length meso = 61.111µm, a mesoscopic mean transverse bulk modulus κ meso = 10.500GPa and a mesoscopic mean transverse shear modulus µ meso = 4.667 GPa, or equivalently to a mesoscopic mean transverse Young's modulus E meso = 12.194 GPa and a mesoscopic mean transverse Poisson's ratio ν meso = 0.3064.The number of iterations n FP required to achieve the desired convergence with the fixed-point iterative algorithm over the mesoscopic subdomain Ω meso is n FP = 5, leading to a number of evaluations of the stochastic computational model equal to n FP tot = 7500.The identification results obtained at mesoscale are also in agreement with the information provided in the literature for this type of biological material.Indeed, from a physical standpoint, the identified spatial correlation length meso = 61.111µm turns out to be of the same order of magnitude as the distance between two adjacent lamellae of an osteon in bovine (beef femur) cortical bone.Moreover, such a value of spatial correlation length is in accordance with the assumption of scale separation between macroscale and mesoscale.

Conclusions
In the present work, we have revisited the multiscale identification methodology recently proposed in Reference [140] for the mechanical characterization of the apparent elastic properties of a complex microstructure made up of a heterogeneous anisotropic material that can be considered as a random linear elastic medium within the framework of 3D linear elasticity theory.Such a multiscale identification has been performed by solving a challenging multiscale statistical inverse problem (requiring multiscale experimental field measurements) formulated as a multi-objective optimization problem.This latter can be decomposed into a first single-objective optimization problem defined at macroscale and a second multi-objective optimization problem defined at mesoscale, to be solved sequentially and involving cost functions (numerical indicators) sufficiently sensitive to the variation of the parameters and hyperparameters to be identified.These numerical indicators allow for quantifying and minimizing the distance between some relevant quantities of interest resulting from the multiscale experimental field measurements at macroscale and mesoscale on the one hand, and their counterparts obtained through forward numerical simulations of a deterministic computational model at macroscale and of a stochastic computational model at mesoscale corresponding to the experimental configuration on the other hand.We consider an ad hoc prior stochastic model introduced in Reference [144] for the numerical modeling and simulation of the random elasticity field, which is parameterized by a small number of hyperparameters.We also employ a stochastic computational homogenization method for the transfer of statistical information from mesoscale to macroscale.The multiscale identification methodology leads to the identification of the optimal values of (i) the parameters involved in the deterministic model of the effective (deterministic and homogeneous) elasticity tensor at macroscale and (ii) the hyperparameters involved in the prior stochastic model of the apparent (random and heterogeneous) elasticity tensor field at mesoscale In the present paper, we have proposed two main improvements of the multiscale statistical inverse identification methodology of the prior stochastic model.First, we have introduced an additional single-objective cost function (numerical indicator) at mesoscale dedicated to the identification of the spatial correlation length(s) involved in the prior stochastic model, allowing the newly formulated multi-objective optimization to be solved with a better computational efficiency by using a (computationally cheap) fixed-point iterative algorithm instead the (costly) global optimization algorithm (genetic algorithm) used in Reference [140].The identification results obtained with the fixed-point iterative algorithm are promising and comparable to that obtained with the genetic algorithm in terms of accuracy.Second, an ad hoc probabilistic modeling of the hyperparameters involved in the prior stochastic model and identified on different mesoscopic domains of observation has been proposed in order to improve both the robustness and the precision of the statistical inverse identification method of the prior stochastic model.Finally, the improved identification methodology has been first validated on in silico materials within the framework of 2D plane stress and 3D linear elasticity with numerically simulated multiscale experimental data, and then successfully applied to real heterogeneous biological material within the framework of 2D plane stress linear elasticity with real multiscale experimental measurements of 2D displacement fields obtained from a static uniaxial compression test performed on a single specimen made of bovine cortical bone and monitored by 2D digital image correlation at both macroscale and mesoscale.In line with this work, several perspectives could be addressed: (i) the multi-objective optimization problem could be solved by using machine learning based on artificial neural networks with a numerical database generated from the stochastic computational model to train an artificial neural network in an (offline) preliminary phase and to use the trained neural network to perform the statistical inverse identification in a computationally cheap (online) computing phase for further reducing the computational cost; (ii) the proposed methodology could be applied to real multiscale experimental measurements of full 3D displacement fields obtained for example by X-ray computed microtomography and digital volume correlation, and also to other types of random heterogeneous materials; (iii) the proposed methodology could be improved by identifying a posterior stochastic model of the non-Gaussian random elasticity (or compliance) field in high stochastic dimension at the mesoscale of an anisotropic heterogeneous linear elastic microstructure using the identified prior stochastic model.

Figure 1 .
Figure 1.Multiscale experimental configuration: displacement field u macro exp measured in the macroscopic domain of observation Ω macro exp

6. 1 .|
Deterministic Macroscopic Boundary Value Problem for the Macroscopic Indicator At macroscale, the deterministic boundary value problem modeling the experimental test configuration described in Section 3 is written over an open bounded domain Ω macro ⊂ R 3 with macroscopic dimensions of the specimen.The experimental domain of observation Ω macro exp is simulated as one given 2D or 3D part Ω macro obs of Ω macro .The boundary ∂Ω macro of Ω macro consists of two disjoint and complementary parts Γ macro N , on which Neumann boundary conditions are applied, and Γ macro D , on which Dirichlet boundary conditions are applied, such that ∂Ω macro = Γ macro N = 0, where |Γ macro D | denotes the 2D measure of Γ macro D .A given deterministic surface force field f macro is applied on Γ macro N , while homogeneous Dirichlet conditions are applied on Γ macro D

Figure 2 .
Boundary value problems at (a) macroscale and (b) mesoscale.(a) Deterministic boundary value problem characterized by deterministic elasticity tensor C macro (a) at macroscale: deterministic displacement field u macro (a) computed at macroscale in Ω macro ; (b) Stochastic boundary value problem characterized by random elasticity tensor field C meso (b) at mesoscale: random displacement field U meso (b) computed at mesoscale in Ω meso .
and the available information allowing for the explicit construction and parametric representation of the probability density function p B : b → p B (b) of random vector B. A robust identified value b opt is finally obtained using the MLE method [68-71] with the independent realizations b meso 1 , . . ., b meso Q .The available information for constructing the prior stochastic model of B is as follows: (i) random variables D, L and H are mutually statistically independent, (ii) random variable D takes its values a.s. in ]0 , δ sup [ with δ sup =

1 , 1 ,
. . ., b meso Q ) is the log-likelihood function for the Q independent realizations b meso . . ., b meso Q of B which is defined for all s ∈ S by L(s; b meso 1 , . . ., b meso Q exp , µ macro exp ) with κ macro exp = 13.901GPa and µ macro exp = 3.685 GPa, corresponding to a Young's modulus E macro exp = 10.158GPa and and a Poisson's ratio ν macro exp = 0.3782.At mesoscale, the solution of stochastic boundary value problem (
. Then, a 2D fourth-order effective elasticity tensor can be defined as C eff 2D (b) = (S eff 2D (b)) −1 .A convergence analysis of the statistical estimator of its statistical fluctuations with respect to the representative volume element size B RVE has been performed.A representative volume element size B RVE = 20× = 400 µm = 4×10 −4 m has been found to be sufficient to reach negligible statistical fluctuations for the construction of the multiscale numerical indicator J multi h (a macro , b) that is calculated by replacing C macro (a) and C eff (b) with C macro 2D (a) and C eff 2D (b), respectively, in (18).

Figure 5 .
Different 2D cross-sections of the Pareto front with the noninferior (Pareto optimal) solutions represented by red stars and the best compromise optimal solution surrounded by a green circle for the mesoscopic domain of observation Ω meso exp,1 .(a) cross-section (J meso δ (b), J meso (b)); (b) cross-section (J meso δ (b), J multi h (a macro , b)); (c) cross-section (J meso (b), J multi h (a macro , b)).

Figure 7 .
Figure 7. Illustration of the test specimen occupying the 3D cubic macroscopic domain of observation Ω macro exp exp = (κ macro exp , µ macro exp ) with κ macro exp = 138.783GPa and µ macro exp = 64.355GPa, corresponding to a Young's modulus E macro exp = 167.218GPa and a Poisson's ratio ν macro exp = 0.2992.

3 ,meso 1 , κ meso 2 , κ meso 3 and µ meso 1 , µ meso 2 , µ meso 3
again reflecting the fact that the two associated mesoscopic numerical indicators J meso δ (b) and J meso (b) depend directly on the experimental field measurements on each in silico test specimen.The identified values κ

Figure 8 .
Fixed-point iterative algorithm: probability density functions p D and p L of random variables D and L, respectively.(a) p D (δ); (b) p L ( ).

Figure 9 .Figure 10 .
Components u macro exp,1 and u macro exp,2 of macroscopic experimental displacement field u macro exp over the 2D macroscopic domain Ω macro exp before and after application of the Gaussian spatial filter.(a) u macro exp,1 unfiltered; (b) u macro exp,1 filtered; (c) u macro exp,2 unfiltered; (d) u macro exp,2 filtered.Components u meso exp,1 and u meso exp,2 of macroscopic experimental displacement field u meso exp over the 2D mesoscopic domain Ω meso exp before and after application of the Gaussian spatial filter.(a) u meso exp,1 unfiltered; (b) u meso exp,1 filtered; (c) u meso exp,2 unfiltered; (d) u meso exp,2 filtered.

11. 2 .
Numerical Results of the Multiscale Statistical Inverse Identification 11.2.1.Resolution of the Single-Objective Optimization Problem at Macroscale

Table 1 .
Comparison between the identified optimal value a macro and the reference experimental value a macro exp .

Table 2 .
Fixed-point iterative algorithm: comparison between the global optimal value b opt obtained from the Q = 16 identified values b meso δ [µm] κ [GPa] µ [GPa] n q

Table 3 .
Fixed-point iterative algorithm: comparison between the global optimal value b opt and the reference experimental value b meso exp for different values of the number N s of independent realizations generated for the statistical estimation of the mathematical expectations involved in the different numerical indicators.

Table 4 .
Genetic algorithm: comparison between the global optimal value b opt obtained from the Q = 16 identified values b meso 1 , . . ., b meso Q for each of the Q mesoscopic domains of observation Ω meso exp,1 , . . ., Ω meso exp,Q , and the reference experimental value b meso exp .δ [µm] κ [GPa] µ [GPa]

Table 5 .
Comparison between the identified optimal value a macro and the reference experimental value a macro exp .

Table 6 .
Fixed-point iterative algorithm: comparison between the global optimal value b opt obtained from the 3 identified values b meso