Inverse Uncertainty Quantification in Material Parameter Calibration Using Probabilistic and Interval Approaches

1. Introduction

Within the calibration of material models, the numerical results of a simulation model are compared with the experimental measurements. One common approach is to minimize the squared errors between the simulation response and the measurements, the so-called least squares approach [1]. This minimization task can be solved by a single-objective optimization algorithm to obtain the optimal parameter set. The optimization procedure requires special attention, if the optimization task is not convex and contains several local minima, as discussed in [2]. A recent overview of optimization techniques applied to solve a calibration task is given in [3,4]. Often a maximum likelihood approach [1] is used to formulate the optimization goal. To ensure a successful identification with optimization techniques, the sensitivity of the model parameters with respect to the considered model outputs should be investigated initially, as discussed in [5].

Additionally to the parameter values itself, the estimation of the parameter uncertainty due to scatter in the experimental measurements and inherent randomness of the material properties requires further investigations. Generally, probabilistic and non-probabilistic techniques have been developed to estimate the parameter uncertainty [6]. A well-known probabilistic method is the maximum likelihood approach [1,7], which uses a straight-forward linearization to represent the parameter scatter as a multi-variate normal distribution. A generalization of this method is the Bayesian updating technique [1], where the joint distribution function of the model parameters could also be non-normally, but can be obtained only implicitly from the model input-output dependencies in the general case. A good overview and recent developments of this technique are given e.g. in [8,9,10,11]. Since the parameter distribution function can be obtained in closed form only for special cases, inverse sampling procedures such as Markov Chain Monte Carlo methods (MCMC) are usually applied to generate discrete samples of the unknown multidimensional parameter distribution [12]. Most of these methods rely on the original Metropolis-Hastings algorithm [13]. However, the MCMC methods usually require a large number of model evaluations and surrogate models are often applied to reduce the numerical burden of the identification procedure [2,14]. Recent developments utilize Bayesian neural networks to represent the forward simulation model including the corresponding parameter and responses uncertainty in a probabilistic approximation model [15,16,17]. Alternatively to the ordinary least squares formulation, Kalman filter identification has been introduced [18] and recently applied for material parameter identification in [19,20]. This technique is a Bayesian estimator using prior information to identify the parameters. However, all the mentioned probabilistic methods require an accurate estimate of the measurement covariance matrix, which is often not available.

Additionally to probabilistic methods, interval approaches have been developed which model the uncertainty of the model inputs and outputs as interval numbers [21,22] or confidence sets [23]. Based on this assumption, interval predictor models have been used for model parameter identification [24]. The unknown model parameter bounds can be obtained for inverse problems from the given model output bounds by interval optimization methods [25,26]. In [27] interval optimization was introduced for fuzzy numbers to obtain a quasi-density of unknown model outputs in forward problems. However, this method can also be applied to inverse problems, which usually results in minimum and maximum values for each parameter considered in the model calibration. In [28,29] bounding boxes have been used to represent the multidimensional domain of the feasible input parameter space. This idea has been extended quite recently to consider correlations between the model parameters [30,31].

In the following paper, we will investigate different methods to estimate the uncertainty of the calibrated parameters from given measurement uncertainty. Using the classical probabilistic approaches as the maximum likelihood estimator and the Bayesian updating procedure as benchmark methods, the scatter of the identified parameters will be estimated as a multi-variate probability density function. Furthermore, we present an interval-based approach by assuming the scatter of the measurements not as correlated random numbers but as individual intervals with known minimum and maximum values. The corresponding possible minimum and maximum values for each input parameter will be obtained using a constrained optimization procedure. Furthermore, the identification of the whole feasible parameter domain for the given measurement bounds is obtained by a novel radial line-search method, which is extended for more than two model parameters using Fekete point sets on a unit hypersphere. With help of this approach, a uniform discretization of the feasible parameter domain boundary is obtained. Furthermore, we put emphasis on the initial sensitivity analysis, which should be considered for all approaches to neglect unimportant parameters in the calibration procedure. As numerical examples, we identify first the material parameters of a bilinear elasto-plasticity model and second the fracture parameters of plain concrete with respect to a wedge splitting test.

2. Materials and Methods

Within the calibration of material models, the numerical results of a simulation model $\mathbf{y}\left(\mathbf{p}\right)$ are compared with the experimental measurements ${\mathbf{y}}^{*}$ as illustrated in Figure 1. One common approach is to minimize the squared errors between the simulation response and the measurements, the so-called least squares approach [1]:

$${\left[{\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right)\right]}^T\left[{\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right)\right]\to min.$$

(1)

This minimization task can be solved by a single-objective optimization algorithm to obtain the optimal parameter set ${\mathbf{p}}_{opt}$. The optimization procedure requires special attention, if the optimization task is not convex and contains several local minima, as discussed in [2].

2.1. Sensitivity Analysis of the Forward Problem

In order to obtain a sufficient accuracy of the calibrated model parameters, a certain sensitivity of these parameters with respect to the model responses is necessary. If all model responses are not sensitive with respect to a specific parameter, this parameter can not be identified within an inverse approach. In such a case, a non-sensitive parameter should be kept fixed and not considered in the calibration process.

A quite common approach for a sensitivity analysis is the global variance-based method using Sobol indicices [32,33]. In this method, the contribution of the variation of the input parameters with respect to the variation of a certain model output can be analyzed. Originally, this method was introduced for sensitivity analysis in the context of uncertainty quantification, but it can similarly apply for optimization problems as discussed in [34]. If we define lower and upper bounds for the m unknown model parameters $p_j$

$$p_i^l\le p_i\le p_i^u,i=1,\dots ,m,$$

(2)

we can assume a uniform distribution for each parameter

$$P_i\sim \mathcal{U}[p_i^l,p_i^u].$$

(3)

The variance contribution of the parameter $P_i$ with respect to a model output $Y_j$ can be quantified with the first order sensitivity measure introduced by [32]

$$S_{P_i}\left(Y_j\right)={\displaystyle \frac{V_{P_i}\left(E_{{\mathbf{P}}_{\sim i}}\left(Y_j\right|P_i)\right)}{V\left(Y_j\right)}},$$

(4)

where $V\left(Y_j\right)$ is the unconditional variance of $Y_j$ and $V_{P_i}\left(E_{{\mathbf{P}}_{\sim i}}\left(Y_j\right|P_i)\right)$ is called the variance of conditional expectation with ${\mathbf{P}}_{\sim i}$ denoting the set of all parameters but $P_i$. $V_{P_i}\left(E_{{\mathbf{P}}_{\sim i}}\left(Y_i\right|P_i)\right)$ measures the first order effect of $P_i$ on the model output. Since first order sensitivity indices quantify only the decoupled influence of each parameter, the total effect sensitivity indices have been introduced by [33] as follows

$$S_{P_i}^T\left(Y_j\right)=1-{\displaystyle \frac{V_{{\mathbf{P}}_{\sim i}}\left(E_{P_i}\left(Y_j\right|{\mathbf{P}}_{\sim i})\right)}{V\left(Y_j\right)}},$$

(5)

where $V_{{\mathbf{P}}_{\sim i}}\left(E_{P_i}\left(Y_j\right|{\mathbf{P}}_{\sim i})\right)$ measures the first order effect of ${\mathbf{P}}_{\sim i}$ on the model output which does not contain any effect corresponding to $P_i$.

Generally, the computation of $S_i$ and $S_i^T$ requires special estimators with a large number of model evaluations [34]. In [35] regression based estimators have been introduced, which require only a single sample set of the model input parameters and the corresponding model outputs. One very simply estimator for the first order index is the squared Pearson correlation coefficient

$$S_{P_i}^{corr}\left(Y_j\right)={\rho }^2(P_i,Y_j).$$

(6)

The Pearson correlation assumes a linear dependence between model input and output. If this is fulfilled, the sum of the first order indices should be close to one. However, if the model output is a strongly nonlinear function of the parameters, the sum is smaller than one and a reliable estimate is not possible with this approach.

With the Metamodel of Optimal Prognosis (MOP) [35] a more accurate meta-model based approach for variance-based sensitivity analysis was introduced, where the quality of the meta-model is estimated with the Coefficient of Prognosis (CoP) by means of a sample cross-validation. If the CoP is close to one, the meta-model can represent the main part of the model output variation and thus the sensitivity estimates are sufficiently accurate. In the MOP approach, different types of meta-models such as classical polynomial response surface models [36,37], Kriging [38], Moving Least Squares [39] as well as artificial neural networks [40] are tested automatically for a specific model output and evaluated by means of the CoP measure. This approach is available with the Ansys optiSLang software package [41]. Further information on the CoP measure can be found in [42]. The samples for the regression based estimators can be generated pure randomly within the defined parameter bounds. However, in order to increase the accuracy of the sensitivity estimates, we apply an improved Latin Hypercube sampling (LHS) [43], which reduces the spurious correlations between the inputs parameters.

2.2. Maximum Likelihood Approach

In the maximum likelihood approach [1,7] the likelihood of the parameters is defined to be proportional to the conditional probability of the measurements ${\mathbf{y}}^{*}$ for a given parameter set $\mathbf{p}$

$$L=k·f\left({\mathbf{y}}^{*}\right|\mathbf{p}).$$

(7)

By assuming, that the chosen model can represent the material behavior perfectly, the remaining gap between model responses and experimental observations $({\mathbf{y}}^{*}-\mathbf{y})$ can be explained only by measurement errors. In this case, $P\left({\mathbf{y}}^{*}\right|\mathbf{p})$ is equivalent to the probability of reproducing the measurement errors. Assuming a multivariate normal distribution for the measurement errors yields [7]

$$P\left({\mathbf{y}}^{*}\right|\mathbf{p})=P({\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right))={\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^m\left|{\mathbf{C}}_{\mathbf{yy}}\right|}}}exp\left[-{\displaystyle \frac{1}{2}}{[{\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right)]}^T{\mathbf{C}}_{\mathbf{yy}}^{-1}[{\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right)]\right],$$

(8)

where ${\mathbf{C}}_{\mathbf{yy}}$ is the covariance matrix of the measurement uncertainty. Maximizing the likelihood L is equivalent to minimize $S=-2lnL$ and yields to the following objective function

$$J\left(\mathbf{p}\right)={\left[{\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right)\right]}^T{\mathbf{C}}_{\mathbf{yy}}^{-1}\left[{\mathbf{y}}^{*}-\mathbf{y}\left(\mathbf{p}\right)\right]\to min.$$

(9)

If the measurement errors $y_i^{*}-y_i$ are assumed to be independent, the objective function equals a weighted least squares formulation

$$J=\sum _{i=1}^n{w}_i{(y_i^{*}-y_i)}^2,w_i={\displaystyle \frac{1}{{\sigma }_{y,i}^2}},$$

(10)

where the weight factor $w_i$ for each measurement value is the inverse of its variance ${\sigma }_{y,i}^2$. For constant measurement errors this objective simplifies to the well-known least squares formulation introduced in Eq. 1. In this case the optimal parameter set ${\mathbf{p}}_{opt}$ can be obtained by the maximum likelihood approach with standard optimization procedures without exact knowledge of the measurement covariance matrix.

The objective function in Eq. 9 can be linearized as follows [7]

$$\Delta \mathbf{p}={\left({\mathbf{A}}^T{\mathbf{C}}_{\mathbf{yy}}^{-1}\mathbf{A}\right)}^{-1}{\mathbf{A}}^T{\mathbf{C}}_{\mathbf{yy}}^{-1}\Delta \mathbf{y},$$

(11)

where $\mathbf{A}$ is a sensitivity matrix, which contains the partial derivatives of the model responses with respect to the parameters

$$\mathbf{A}\left(\mathbf{p}\right)={\displaystyle \frac{\partial \mathbf{y}\left(\mathbf{p}\right)}{\partial \mathbf{p}}}.$$

(12)

Assuming, that the optimal parameter set ${\mathbf{p}}_{opt}$ has been found by minimizing the objective function, we can estimate the sensitivity matrix at the optimal parameter set e.g. by a central differences method.

Assuming a linear relation between the model parameters and responses, the covariance matrix of the parameters can be estimated with the so-called Markov estimator [1,7]

$${\mathbf{C}}_{\mathbf{pp}}={\left({\mathbf{A}}_{opt}^T{\mathbf{C}}_{\mathbf{yy}}^{-1}{\mathbf{A}}_{opt}\right)}^{-1}.$$

(13)

This covariance matrix ${\mathbf{C}}_{\mathbf{pp}}$ contains the necessary information about the scatter of the individual parameters, as well as possible pairwise correlations. The distribution of the parameters is Gaussian as assumed for the measurements. The estimated scatter and correlations are a very useful information to judge about the accuracy and uniqueness of the identified parameters in a calibration process. This method can be interpreted as an inverse first order second moment method: assuming a linear relation between parameters and measurements and a joint normal distribution of the measurements, we can directly estimate the joint normal covariance of the parameters. However, the application of this linearization scheme requires an accurate estimate of the optimal parameter set ${\mathbf{p}}_{opt}$ by a previous least squares minimization. Furthermore, the inversion of the matrix in Eq. 13 requires, that every parameter is sensitive w.r.t. the model output. This should be investigated with an initial sensitivity analysis, as explained in Section 2.1.

The linearization in Eq. 11 is of course a simplification of the general nonlinear case and contains only the first term of the Taylor series expansion of the non-linear model responses. We will investigate the influence of this linearization by means of the second example in this paper. Further discussion on applicability of this approach for nonlinear models can be found in [7] and [2].

The evaluation of Eq. 13 requires an estimate of the measurement covariance matrix ${\mathbf{C}}_{\mathbf{yy}}$. If a measurement procedure is repeated very often, this matrix can be estimated directly from measurement statistics. However, this is often not possible. A valid estimate could be to assume independent measurement errors. The resulting covariance matrix is then diagonal, where at least the scatter of the individual measurement points has to be estimated. From a few measurement repetitions only, this can be estimated straight forward. If even such information is not available, the measurement scatter can be assumed to be constant ${\sigma }_y$ for each output value

$${\mathbf{C}}_{\mathbf{yy}}={\sigma }_y^2\mathbf{I}.$$

(14)

In this case, the estimated covariances of the parameters are directly scaled by the assumed measurement scatter ${\sigma }_y$

$${\mathbf{C}}_{\mathbf{pp}}={\sigma }_y^2{\left({\mathbf{A}}_{opt}^T{\mathbf{A}}_{opt}\right)}^{-1}.$$

(15)

If only a single measurement curve is available for the investigated responses, the misfit of the calibration process can be used for a rough estimate. Assuming the measurement errors to be constant and independent, we obtained the following estimate

$${\sigma }_y\approx RMSE=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-y_i^{*})}^2},$$

(16)

where n is the number of measurement points.

The sensitivity matrix in Eq. 12 can be visualized as shown in the result section by normalizing the partial derivatives with respect to the parameter ranges

$${\tilde{A}}_{y_j}\left(p_i\right)={\displaystyle \frac{\partial y_j\left(\mathbf{p}\right)}{\partial p_i}}·(p_i^u-p_i^l).$$

(17)

By using the variance formula for a sum of linear factors, we can estimate local variance-based first order sensitivity indices from this partial derivatives by assuming again a uniform distribution of the parameters within the defined bounds

$$S_{P_i}^{local}\left(Y_j\right)={\displaystyle \frac{{\left(\frac{\partial y_j\left(\mathbf{p}\right)}{\partial p_i}\right)}^2·V\left(P_i\right)}{{\sum }_{k=1}^m{\left(\frac{\partial y_j\left(\mathbf{p}\right)}{\partial p_k}\right)}^2·V\left(P_k\right)}}={\displaystyle \frac{{\tilde{A}}_{y_j}^2\left(p_i\right)}{{\sum }_{k=1}^m{\tilde{A}}_{y_j}^2\left(p_k\right)}}.$$

(18)

2.3. Bayesian Updating

Another probabilistic method is the Bayesian updating approach [1], which can represent non-linear relations between the parameters and the responses but requires much more numerical effort to obtain the statistical estimates. In this approach, the unknown conditional distribution of the parameters is formulated with respect to the given measurements by using the Bayes’ theorem as follows

$$P\left(\mathbf{p}\right|{\mathbf{y}}^{*})={\displaystyle \frac{P\left({\mathbf{y}}^{*}\right|\mathbf{p})P\left(\mathbf{p}\right)}{P\left({\mathbf{y}}^{*}\right)}}.$$

(19)

The likelihood function $P\left({\mathbf{y}}^{*}\right|\mathbf{p})$ can be formulated as probability density function, whereby often a multi-normal distribution function is utilized similarly as by Markov estimator

$$P\left({\mathbf{y}}^{*}\right|\mathbf{p})={\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^m\left|{\mathbf{C}}_{\mathbf{yy}}\right|}}}exp\left[-{\displaystyle \frac{1}{2}}{({\mathbf{y}}^{*}-\mathbf{y})}^T{\mathbf{C}}_{\mathbf{yy}}^{-1}({\mathbf{y}}^{*}-\mathbf{y})\right].$$

(20)

The so-called prior distribution $P\left(\mathbf{p}\right)$ can be assumed either from previous information or as a uniform distribution within initially defined parameter bounds. Since the normalization term $P\left({\mathbf{y}}^{*}\right)$ is usually not known, the posterior parameter distribution $P\left(\mathbf{p}\right|{\mathbf{y}}^{*})$ could not be obtained in closed form. However, in such a case, Markov Chain Monte Carlo methods could be used to generate discrete samples of the posterior distribution. A recent overview of these methods is given in [12]. For further details on Bayesian updating, the interested reader is referred to [6]. In our study, the well-known Metropolis-Hastings algorithm [13] is applied, which considers only the forward evaluation of the likelihood function to generate samples of the posterior distribution. Detailed benchmark investigations on the Metropolis-Hastings algorithm as well as an efficient approximation of the likelihood function can be found in [2].

2.4. Interval Optimization

Instead of modeling the measurement points as correlated random numbers, which requires the estimate of the full covariance matrix, we can assume, that each measurement point is an interval number with known minimum and maximum bounds as shown in Figure 2. For such a given interval of each measurement point, we want to obtain the corresponding range of the model parameters which fulfill

$$y_i^l\le y_i(p_1,\dots ,p_m)\le y_i^u,i=1,\dots ,n.$$

(21)

In the classical interval optimization approach, usually the interval of an uncertain response is obtained from the known interval of the input parameters by several forward optimization runs. A specific approach for this purpose, the $\alpha$-level optimization [27], as shown in principle in Figure 2, where the most probable value is assumed for each input parameter with the maximum $\alpha$-level and a linear descent to zero is defined between the maximum and the lower and upper bounds. The $\alpha$-level optimization will detect the shape of the corresponding $\alpha$-levels of the unknown outputs. In an inverse approach, we need to search for the possible minimum and maximum values of the input parameters $p_j$ for a given level, while each response value $y_i$ is constrained to the corresponding interval as defined in Eq. 21. The corresponding constraint optimization task can be formulated for each input parameter $p_j$ as an individual minimization and maximization task as follows

$$\begin{array}{cc}\hfill p_j& \to min,max\hfill \\ \hfill with& p_k^l\le p_k\le p_k^u,\hfill & \hfill k=& 1,\dots ,m\hfill \\ \hfill & y_i^l\le y_i(p_1,\dots ,p_m)\le y_i^u,\hfill & \hfill i=& 1,\dots ,n,\hfill \end{array}$$

(22)

where the specific $\alpha$-level defines the lower and upper bounds $y_i^l$ and $y_i^u$ of each response value. This approach detects the lower and upper parameter bounds only and will not investigate possible interactions in the multidimensional parameters space. Therefore, we present a new radial line-search approach in the next section, which investigates the whole feasible domain of the possible model parameter values.

2.5. Radial Line-Search Approach

Assuming interval constrains for each response value $y_i$ according to Eq. 21, the corresponding domain of the feasible input parameter domain $\mathbb{A}$ is defined as

$$\mathbb{A}=\left\{{\mathbf{p}}_A\right|y_i^l\le y_i\left({\mathbf{p}}_A\right)\le y_i^u\forall 1\le i\le n\}.$$

(23)

This domain of feasible parameter sets ${\mathbf{p}}_A$ could be analyzed in case of two input parameters with a line-search approach as illustrated in Figure 3: starting from the optimal parameter set ${\mathbf{p}}_{opt}\in \mathbb{A}$, we search in every direction for the maximum radius of the feasible parameter domain

$${\mathbf{p}}_l=r_l·{\mathbf{a}}_l+{\mathbf{p}}_{opt},r_l\to max\mathrm{subjected}\mathrm{to}{\mathbf{p}}_l\in \mathbb{A},$$

(24)

where ${\mathbf{a}}_l$ is a direction vector of unit length

$$\parallel {\mathbf{a}}_l\parallel =1.$$

(25)

The direction vectors ${\mathbf{a}}_l$ can be chosen randomly or uniformly distributed as shown in Figure 3 for two input parameters. For each direction a line-search is performed in order to obtain the maximum possible radius $r_l$ for this direction. This search can be realized by bisection or other methods. The obtained boundary points ${\mathbf{p}}_l$ result in a uniform discretization of the feasible parameter domain $\mathbb{A}$. This approach is similar to the directional sampling method developed for probability integration in reliability analysis [44,45].

In order to extend this method to more than two input parameters, the search directions ${\mathbf{a}}_l$ have to be discretized on a 3D- or higher dimensional (hyper-) sphere. One very efficient discretization method to obtain almost equally distributed points on the hyper-sphere are the Fekete point sets [46], which are generated using a minimum potential energy criterion

$$\mathsf{\Pi }=\sum _{1\le i<j\le n_f}{\parallel {\mathbf{a}}_i-{\mathbf{a}}_j\parallel }^{-2}\to min,$$

(26)

where $n_f$ is the number of points on the unit hyper-sphere. The minimization of the potential energy results in equally distributed point sets and will require a significant smaller number of model evaluations to obtain a suitable representation of the parameter domain boundary. An efficient optimization procedure to solve Eq. 26 is proposed in [46]. Once the initial points have been generated, a line-search is performed for each of the search directions ${\mathbf{a}}_l$ similarly as in two dimensions. This search procedure is performed in our study by using a simple bisection method [47]. However, since every direction can be evaluated independently, a parallel computing of the required model evaluations is possible. The obtained points on the boundary can be used to visualize a discretized shape of the feasible parameter domain. If more than two or three parameters are investigated, this visualization could be displayed in 2D or 3D subspaces.

3. Results

3.1. Bilinear Elasto-Plasticity Model

As a first example, we will investigate a uniaxial bilinear elasto-plasticity model defined as follows

$$\sigma \left(ϵ\right)=\left\{\begin{array}{cc}E·ϵ& ϵ\le {\displaystyle \frac{E}{{\sigma }_Y}}\\ {\sigma }_Y+H·\left(ϵ-{\displaystyle \frac{{\sigma }_Y}{E}}\right)& ϵ>{\displaystyle \frac{E}{{\sigma }_Y}}\end{array}\right.$$

(27)

where $\sigma \left(ϵ\right)$ is the uniaxial stress depending on the uniaxial strain $ϵ$. The model parameters are the Young’s modulus E, the yield stress ${\sigma }_Y$ and the hardening modulus H. Figure 4 shows a reference curve for a selected parameter set $E_{ref}=1000N/m{m}^2$, ${\sigma }_{Y,ref}=4.0N/m{m}^2$ and $H_{ref}=100N/m{m}^2$, which is used in the following investigations.

3.1.1. Deterministic Calibration with Noisy Measurements

By assuming additional uncertainties in the measurement points $y_i^{*}$, the obtained optimal parameter set could significantly deviate from the reference parameters as shown in Figure 4, where an additional Gaussian noise term is used to generate synthetic measurement points

$$y_i^{*}=\sigma ({ϵ}_i,E_{ref},{\sigma }_{Y,ref},{\sigma }_{Y,ref})+\alpha ·\mathcal{N}(0,1),$$

(28)

where $\mathcal{N}(0,1)$ is a normally distributed random number with zero mean and unit variance. If we repeat the calibration of the model with different samples of the noise term, the calibrated optimal parameters will show a certain scatter, which depends on the size of the noise term. 1000 samples of the calibrated parameters are shown exemplarily in Figure 5 for a noise term of $\alpha =0.1N/m{m}^2$ by considering 20 measurement points as shown in Figure 4. The figure indicates, that all three identified parameters show significant scatter and the yield stress and the hardening modulus show a significant correlation to each other.

3.1.2. Sensitivity Analysis

By means of the bilinear example, we will apply the presented probabilistic approaches and the interval optimization methods. Initially, we perform a sensitivity analysis to investigate the different estimators discussed in Section 2.1. For this purpose, we calculate the sensitivity indices from the squared Pearson correlations by evaluating 20 points on the stress-strain curve. 1000 Latin Hypercube samples have been generated by considering the parameter bounds of $E\in [800,1200]N/m{m}^2$, ${\sigma }_Y\in [2.0,6.0]N/m{m}^2$ and $H\in [50,500]N/m{m}^2$. The obtained sensitivity indices shown in Figure 6 indicate, that the Young’s modulus has a significant influence for smaller strain values whereas the yield stress and the hardening modulus become more important with increasing strains. However, the sum of the individual sensitivity indices is not close to one for strain values between 0.2 and 0.6 % which is caused be the non-linear relation between these stress values with respect to the model parameters. The sensitivity indices by using the MOP approach are evaluated in a second step using the same 1000 samples. As indicated in figure, these estimates show a similar behavior as the correlation based estimates. However, the approximation quality of the underlying meta-model is almost excellent for all stress values. Thus, the estimated total effect indices are more reliable measures to judge on the sensitivity of the model parameters for this example.

3.1.3. Maximum Likelihood Approach

In a next step, the Markov estimator is applied using Eq. 13 to estimate the parameter covariance matrix from the measurement covariance matrix. Again, 20 measurement points are considered on the hardening curve with a constant measurement error of ${\sigma }_y=0.1N/m{m}^2$. The covariance matrix is assembled according to Eq. 14. The required sensitivity matrix in Eq. 12 is estimated by using central difference at the optimal parameter set, which requires only 6 additional model evaluations. In Figure 7 the estimated scatter of the parameters is visualized by generating 1000 samples of the estimated parameter covariances. The figure indicates a similar scatter and correlation of the parameters as obtained with the 1000 calibration runs in Figure 5. Figure 8 shows the normalized derivatives of the sensitivity matrix.

In a next step, the influence of the assumed measurement error ${\sigma }_y$ and the number of measurement points $y_i$ is investigated. Figure 9 shows the estimated parameter scatter normalized with the parameter ranges depending on the size of an independent measurement error according to Eq. 14 for 20 measurement points. The figure indicates a linear relation of the standard deviation for all three parameters, as expected from the formulation in Eq. 15. The normalized parameter scatter of the yield stress is significantly smaller compared to the scatter of the Young’s and hardening moduli. This is caused by the fact, that sensitivity of the yield stress is on average larger as observed in the normalized partial derivatives in Figure 8. Additionally, the influence of the number of measurement points is shown in Figure 9 for a constant ${\sigma }_y=0.1N/m{m}^2$. This figure clearly indicates, that with an increasing number of independent measurement points, the estimated scatter and thus the uncertainty decreases. This problem is investigated in detail in [48] and could be resolved only, if a realistic correlation between the measurement points is assumed in the covariance matrix.

3.1.4. Bayesian Updating

The Bayesian updating approach is investigated next by using a constant measurement error ${\sigma }_y=0.1N/m{m}^2$ and a diagonal covariance matrix as within the previous investigations with Markov estimator. The posterior distribution in Eq. 19 is evaluated for the bilinear model using the Metropolis-Hastings algorithm with 20.000 samples, whereby the first 2.000 samples of the so-called burn-in phase are neglected in the statistical analysis. In Figure 10 the estimated parameter scatter is shown. The pair-wise anthill plots and the correlation matrix indicate a good agreement with the results of the Markov estimator shown in Figure 7. This may be explained by an approximately linear behavior of the model for small parameter variations.

In Figure 11 the estimated mean values and the variation of the model response values are visualized for the generated samples. The figure indicates a good agreement of the mean response values with the deterministic reference curve. Around a deterministic yield strain of $0.4\%$ a larger deviation could be observed. More interesting is the standard deviation of the individual stress responses. In contrast to the assumption of a constant measurement error of ${\sigma }_y=0.1N/m{m}^2$, the variation of all stress values is smaller caused by the inherent correlations of the model responses.

3.1.5. Interval Optimization and Radial Line-Search

Next, the presented interval approaches are applied for the bilinear model by considering a constant measurement interval size of $\Delta y_i=\pm 0.2N/m{m}^2$. First, the minimum and maximum parameter values are obtained with 2 optimization runs for each model parameter by using an evolutionary algorithm from the Ansys optiSLang software package [41] with 25 generations and a population size of 20. The obtained best parameter sets of each optimization run are shown in Figure 12, which span a three-dimensional box around the reference parameter set.

Second, the radial line-search approach is applied for different 2D subspaces by considering two model parameters as optimization variables while the third is kept at the reference value. In Figure 13 the initial radial search points and the finally obtained boundary points are shown for all 2D subspaces. The initial parameter bounds have been chosen from the previous optimization results as $E\in [900,1100]N/m{m}^2$, ${\sigma }_Y\in [3.5,4.5]N/m{m}^2$, and $E\in [0,200]N/m{m}^2$ and 50 uniformly distributed search directions are considered for each subspace. The subspace plots in Figure 13 indicate an almost rectangular boundary of the feasible Young’modulus-yield stress subspace and of the Young’s modulus - hardening modulus subspace, meanwhile in the third subspace, the yield stress - hardening modulus subspace, the boundary indicates a significant dependence of the feasible parameter values.

Finally, the bilinear example is investigated with the radial line-search with all three parameters. For this purpose, 500 uniformly distributed Fekete points are generated on the unit sphere by considering the parameter bounds of the 2D analysis. The line-search is performed for each direction individually with a bisection with 10 refinement steps. In Figure 14 the initial search points are shown together with the obtained boundary points for all 500 directions. In the figure, a Delaunay triangulation of the convex hull of boundary points is used for better visualization purposes by utilizing the MATLAB plotting library [49]. The obtained feasible parameter domain indicates a dependence between yield stress and hardening modulus similar to the 2D analyses. The ranges of the response values of the obtained boundary points are given in Figure 15, where the most values approach to the $\pm 0.2N/m{m}^2$ interval. For strain values below $0.4\%$ the maximum and minimum response values decrease linearly due to the pure elastic behavior and the resulting perfect correlation of the stress-strain curves.

3.2. Tension Softening Model for Concrete Cracking

3.2.1. Sensitivity Analysis and Deterministic Parameter Optimization

By means of the second numerical example, the material parameters of plain concrete are identified using the measurement results of a wedge splitting test published in [50]. The geometry of the investigated concrete specimen is shown in Figure 16. The softening curve is obtained by a displacement-controlled simulation. The simulation model is discretized by 2D finite elements and considers an elastic base material and a predefined crack with bilinear softening law. The displacements are measured as the relative displacements between the load application points. Further details about the simulation model can be found in [51].

Six unknown parameters are identified with the presented methods: the Young’s modulus E, the Poisson’s ratio $\nu$, the tensile strength $f_t$, the Mode-I fracture energy $G_f$ and the two shape parameters ${\alpha }_{ft}$ and ${\alpha }_{wc}$, which define the kink of the bilinear softening law in terms of the tensile stress ${\sigma }_1={\alpha }_{ft}·f_t$ and the relative displacements normal to the crack surface $\Delta u_{N1}={\alpha }_{wc}·\Delta u_{Nc}$. In Table 1 the reference parameters according to [50] and the defined parameter bounds for the calibration are given.

In a first step, the sensitivity of the material parameters with respect to the simulated force values is investigated by generating 200 samples within the defined bounds by using an improved Latin-Hypercube Sampling according to [43]. The corresponding force-displacement curves of these samples are given in Figure 17 additionally to the experimental values according to [50]. The total effect sensitivity indices of the model parameters are estimated for each force value on the load-displacement curve by using the Metamodel of Optimal Prognosis [35] approach, as discussed in Section 2.1. These sensitivity indices are shown in Figure 18 depending on the displacement values. The figure indicates, that the Poisson’ ratio is not sensitive to any of the simulation results. Because of this finding, we do not consider it in the following parameter calibration and uncertainty estimation. The other parameters show a certain influence whereby the Young’s modulus, the tensile strength and the fracture energy are most important.

The optimal parameter set is determined by a least-squares minimization using Eq. 1, whereby the deviation of the response values is evaluated at the 13 displacement values of the experimental measurement points. As optimization approach, an Adaptive Response Surface Method of the Ansys optiSLang software package [41] is utilized, which required 220 model evaluations. The obtained optimal parameter values are given in Table 1 and the corresponding load-displacement curve is shown in Figure 19, which indicates almost perfect agreement with the experimental values.

3.2.2. Probabilistic Uncertainty Estimation

In a next step, the Markov estimator and the Bayesian updating approaches are applied to estimate the parameter uncertainties. The covariance matrix of the 13 measurement points are defined accordingly to Eq. 14 as a diagonal matrix, whereby a constant standard deviation is assumed as the root mean squared error of the deterministic calibration as

$${\sigma }_y=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-y_i^{*})}^2}=151.4N.$$

The estimated parameter uncertainty using the Markov estimator are given in Table 2. The mean values are the optimal parameter values from the previous deterministic optimization. The standard deviation and the parameter correlation matrix are obtained from the estimated covariance matrix by using Eq. 15. The partial derivatives required in Eq. 12 are estimated using a central differences approach with an interval of $0.01\%$ of the initial parameter ranges. Thus, only 10 additional model evaluations were necessary to perform the Markov estimator for this example.

The results in Table 2 indicate, that the shape parameters and the Young’s modulus have a significantly larger coefficient of variation (CoV) as the tensile strength and the fracture energy. This can be explained, if we compare the sensitivity indices of the parameters given in Figure 18, where the fracture energy and the tensile strength have in average the largest values. The estimated parameter correlation matrix is given in Figure 20, which indicate a significant correlation between the Young’s modulus and the tensile strength. Further dependencies could be observed between the fracture energy and one of the shape parameters. Additionally, the figure shows some interesting pair-wise anthill plots with 10.000 samples generated with the multi-variate normal distribution considering the estimated parameter covariance matrix of the Markov estimator.

Next, the Bayesian updating is performed with the same diagonal measurement covariance matrix used for the Markov estimator. 50.000 samples are generated with the Metropolis-Hastings algorithm, whereby 5.000 samples of the burn-in phase are neglected in the statistical analysis. The remaining 45.000 samples are shown additionally in Figure 20. In Table 2 the corresponding statistical estimates are given. The table and the pair-wise anthill plots indicate, that the mean values differ significantly from the Markov estimator results for the Young’s modulus and the tensile strength. Good agreement of the statistical properties is obtained for the fracture energy and the shape parameters. Even the estimated parameter correlations are quite close to the Markov estimator results.

3.2.3. Interval Approaches

Finally, the proposed radial line-search algorithm is demonstrated by means of this example. In a first step, only the Young’s modulus and the tensile strength are considered as uncertain parameters, meanwhile the remaining parameters are kept constant. The resulting two-dimensional search points are shown in Figure 21 for an assumed measurement interval of $\Delta y_i=\pm 0.5kN$, which corresponds to 3.3 times the RMSE. This interval is equivalent to the $99.9\%$ confidence interval of the previously assumed measurement uncertainty. The figure clearly indicates, that the feasible 2D parameter domain is non-symmetric and bounded by different linear and non-linear segments. Different $\alpha$-levels could be investigated, if the measurement interval is modified in the analysis. Figure 21 shows additionally the domain boundary for different values of $\Delta y_i$, where a non-linear shrinkage of the domain boundary could be observed. The figure clearly indicates, that the obtained parameter domain boundary is not symmetric as assumed by the Markov estimator.

In a second step, all five sensitive parameters are considered in the radial line-search. 5.000 directions have been investigated in this analysis, where the ranges of the hyper-sphere are given in Table 3. The table shows additionally the minimum and maximum possible parameters bounds obtained by the interval optimization approach by assuming the measurement interval as $\Delta y_i=\pm 0.5kN$. For each minimization and maximization task, again the evolutionary algorithm from the Ansys optiSLang software package [41] with 25 generations and a population size of 20 was applied. The minimum and maximum parameter values are shown together with the obtained 5.000 boundary points of the radial line-search approach in Figure 22 for some selected parameter subspaces. The figure indicates for the Young’s modulus and tensile strength a non-symmetric dependence. The dependence between the other parameters seems to be similar as the results obtained with the Bayesian updating.

Since the evaluation of the radial line search requires a huge number of model evaluations, we will investigate the accuracy of different approximation methods by means of this example. First a linear approximation model from the 10 samples of the central differences approach within an interval of $\pm 0.05\%$ of the initial parameter ranges is investigated. We apply the radial line-search on this 5D approximation model as investigated previously only in the 2D subspace of the Young’s modulus and the tensile strength as uncertain parameters, meanwhile the remaining parameters are kept constant. In Figure 23 the results by using the original FEM model are compared to results obtained with the linear approximation model. The figure indicates already a sufficient accuracy, which represents the main characteristics of the feasible parameter domain. The 2D analysis has been chosen to enable a more clear visualization of the approximation errors.

Additionally, we apply a full quadratic approximation model based on 100 Latin-Hypercube samples generated within the $\pm 0.05\%$ domain of all 5 parameters around the optimal parameter values ${\mathbf{p}}_{opt}$. The results of the 2D line-search shown in Figure 23 are quite similar as obtained with the linear model, but still show a certain deviation from the FEM solutions. In order to get a better representation of the non-linear dependencies of the model responses and the model parameters in the range of the parameter scatter, we apply the Metamodel of Optimal Prognosis as a more sophisticated surrogate on the response values obtained from the 100 LHS samples within the $\pm 0.05\%$ parameter ranges. As indicated additionally in Figure 23, the straight forward application of the MOP already yields to quite accurate results. An additional improvement of the surrogate model could be achieved by considering the obtained 100 boundary points of the line-search procedure as additional support points within the surrogate model, which results in an almost perfect agreement of the line-search results.

4. Conclusions

In the presented paper we have investigated well-known probabilistic approaches as benchmark methods for the uncertainty quantification of calibrated material parameters. The applied Markov estimator and the Bayesian updating require the knowledge of the measurement covariances to estimate the multi-variate probability density function of the calibrated model parameters. Often this matrix can be estimated quite roughly based on few measurement curves. If the covariance matrix is assumed as a diagonal matrix this implies independent measurement errors and the estimated parameter scatter directly depends on the number of considered points on the measurement curve. As a consequence, the obtained results might be difficult to interpret and a successful application of the probabilistic approaches would require a more accurate estimate of the measurement covariance matrix.

In order to overcome this limitation, an interval optimization approach which assumes only the maximum and minimum values of the measurement points was extended in the presented paper. With help of this approach, the full feasible parameter domain fulfilling the defined measurement intervals could be determined by a novel radial line-search procedure, which will search for the boundary of the feasible parameter space based on the underlying simulation model. The boundary domain is discretized with uniformly distributed Fekete points, which define the radial directions based on the optimal parameter set obtained in the previous deterministic calibration. For each direction, the boundary point was obtained by a one-dimensional line-search procedure. Compared to previous interval search procedures as published in [26,30,31], the presented approach does not require any assumption on the parameter correlation.

However, the computation requires numerous model evaluations but can be computed in parallel for each independent search direction. The computational burden could be reduced significantly by using either low order approximation schemes or by applying more sophisticated surrogate models on support points around the optimal parameter set. Here is place for further research in order to determine a fully automatic adaptive scheme based on surrogate models.

As a result of the presented approach, we obtain the boundary of the feasible parameter domain by means of a given number of discretization points in the parameter input space. The transfer to a closed description of the parameter uncertainty e.g. with correlated interval numbers is an open question which would motivate further research.