Inverse Uncertainty Quantification in Model Calibration Using Probabilistic and Interval Approaches

1. Introduction

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,3,4]. A recent overview of optimization techniques to solve the calibration task is given in [5,6]. 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 [7].

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. In [1,8] a straight-forward linearization based on the maximum likelihood approach was introduced, where the parameter scatter is assumed 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 can be obtained only implicitly from the model input-output dependencies in many cases. A good overview and recent developments of this technique are given e.g. in [9,10,11,12,13]. 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 multi-dimensional parameter distribution [14]. Most of these methods rely on the original Metropolis-Hastings algorithm [15]. However, the MCMC methods require usually a large number of model evaluations and meta-models are often applied to accelerate the identification procedure [2,16]. Additionally, the Bayesian updating approach requires an accurate estimate of the measurement covariance matrix which is often not available. In the following paper we want to investigate this problem in detail.

Additionally to probabilistic methods, interval approaches have been developed which model the uncertainty of the model inputs and outputs as interval numbers [17,18] or confidence sets [19]. Based on this assumption, interval predictor models have been used for model parameter identification [20]. The unknown model parameter bounds can be obtained for inverse problems from the given model output bounds by interval optimization methods [21,22]. In [23] 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 results usually in minimum and maximum values for each parameter considered in the model calibration. In [24,25] bounding boxes have been used to represent the multi-dimensional domain of the feasible input parameter space. This idea has been extended quite recently to consider correlations between the model parameters [26,27].

Let us introduce the following uniaxial bilinear elasto-plasticity model as an illustrative example

$$\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.$$

(2)

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. In Figure 2 the reference curve is given 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$.

If we assume, that the measurement points $y_i^{*}$ contain some scatter or uncertainty, the calibrated optimal parameter set could show significant deviation from the reference parameters as shown in Figure 2, 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),$$

(3)

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 3 for a of noise term of $\alpha =0.1N/m{m}^2$ by considering 20 measurement points as shown in Figure 2. The figure indicates, that all three identified parameter show significant scatter and the yield stress and the hardening modulus show a significant correlation to each other.

In the following paper, we will investigate different methods to estimate the uncertainty of the calibrated parameters from given measurement uncertainty. Based on classical probabilistic approaches as the maximum likelihood estimator or the Bayesian updating procedure, the scatter of the identified parameters can be estimated as a multi-variate probability density function. Both procedures require a sufficient accurate knowledge or estimate of the scatter of the measurement points, often modeled by a Gaussian covariance matrix. Furthermore, we present an additional 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 can be obtained for this assumption using a constrained optimization procedure. However, the identification of the whole possible parameter domain for the given measurement bounds is not straight forward. Therefore, we introduce an efficient line-search method, where not only the bounds itself but also the shape of the feasible parameter domain can be identified. This enables the quantification of the interaction and the uniqueness of the input parameters. 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 a numerical example, we identify the fracture parameters of plain concrete with respect to a wedge splitting test, where six unknown material parameters should be identified.

2. Sensitivity Analysis of the Model Parameters

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, this parameter should be kept fixed and not considered in the calibration process.

A quite common approach for a sensitivity analysis is a global variance-based method using Sobol indicices [28,29]. 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 applied for optimization problems as discussed in [30].

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$$

(4)

we can assume a uniform distribution for each parameter

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

(5)

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 [28]

$$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)}},$$

(6)

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 [29] 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)}},$$

(7)

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 [31]. In [32,33] 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).$$

(8)

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) [33] 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 [34,35], Kriging [36], Moving Least Squares [37] as well as artificial neural networks [38] 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 [39]. Further information on the overall framework and the CoP measure can be found in [40,41]. 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) [42], which reduces the spurious correlations between the inputs parameters.

In Figure 4 the sensitivity indices from the squared Pearson correlations are shown for the bilinear example by evaluating 20 points on the stress-strain curve. 1000 Latin Hypercube samples have been generated by considering the following parameter bounds: $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 figure clearly indicates, that the Young’s modulus shows 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 indeces 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 results of the sensitivity analysis by using the MOP approach are given additionally in the figure and 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 indeces are more reliable measures to judge on the sensitivity of the model parameters for this example.

3. Probabilistic Approaches

3.1. Maximum Likelihood Approach

In the maximum likelihood approach [1,8] 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}).$$

(9)

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 [8]

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

(10)

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.$$

(11)

If the measurement errors $y_i^{*}-y_i$ are assumed to 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}},$$

(12)

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}$ with the maximum likelihood approach can be obtained by standard optimization procedures without exact knowledge of the measurement covariance matrix.

3.2. Markov estimator

The objective function in Eq. 11 can be linearized as follows [8]

$$\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},$$

(13)

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}}}.$$

(14)

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,8]

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

(15)

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. 15 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.

The evaluation of Eq. 15 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 small number of measurement repetitions 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}.$$

(16)

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}.$$

(17)

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},$$

(18)

where n is the number of measurement points.

In Figure 5 the estimated scatter of the parameters of the bilinear model is shown. Again 20 measurement points with ${\sigma }_y=0.1N/m{m}^2$ are assumed, which results in a covariance matrix according to Eq. 16. The figure indicates a similar scatter and correlation of the parameters as obtained with the 1000 calibration runs in Figure 3. For the calculation of the sensitivity matrix in Eq. 14 only 6 model evaluations are necessary.

The sensitivity matrix in Eq. 14 can be visualized as shown in Figure 6 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).$$

(19)

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)}}.$$

(20)

In Figure 7 the estimated parameter scatter normalized with the parameter ranges is shown depending on the size of an independent measurement error according to Eq. 16 for 20 measurement points. The figure indicates a linear relation of the standard deviation of all parameters as expected from the formulation in Eq. 17. 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 6. Additionally to the influence of the measurement error, the influence of the number of measurement points is shown in Figure 7 for a constant ${\sigma }_y=0.1N/m{m}^2$. The figure clearly indicates, that with an increasing number of independent measurement points, the estimated scatter and thus the uncertainty decreases. This problem is discussed in [43] and could be resolved only, if a realistic correlation between the measurement points is assumed in the covariance matrix.

3.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)}}.$$

(21)

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

(22)

The so-called prior distribution $P\left(\mathbf{p}\right)$ can be assumed either from previous information or as an 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 [14]. For further details on Bayesian updating the interested reader is referred to [9]. In our study the well-known Metropolis-Hastings algorithm [15] is applied, which considers only the forward evaluation of the likelihood function to generate samples of the posterior distribution. Further details on this algorithm can be found in [2].

In Figure 8 the estimated parameter scatter of the bilinear model is shown by assuming again independent measurement errors of ${\sigma }_y=0.1N/m{m}^2$. 20.000 samples have been generated by the Metropolis-Hastings algorithm by neglecting the first 2.000 samples, the so-called burn-in phase. The pair-wise anthill plots and the correlation matrix indicate a good agreement with the results of the Markov estimator in Figure 5 for this example. This can be explained by the approximately linear behavior of the model for small measurement errors.

In Figure 9 the corresponding variation of the response values are given. 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.

4. Interval Approaches

4.1. 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 10. 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.$$

(23)

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 [23], is shown in principle in Figure 10, 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. 23. 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}$$

(24)

where the specific $\alpha$-level defines the lower and upper bounds $y_i^l$ and $y_i^u$ of each response value.

In Figure 11 the obtained minimum and maximum parameter values are shown exemplarily for the bilinear model by considering a constant measurement interval size of $\Delta y_i=\pm 0.2N/m{m}^2$. Every optimization run was carried out by using an evolutionary algorithm from the Ansys optiSLang software package [39] with 25 generations and a population size of 20. The obtained best parameter sets of each optimization run are shown in the figure, which span a three-dimensional box around the reference parameter set. However, the figure illustrates that the detection of the lower and upper parameter values alone is not sufficient to investigate the interactions between the input parameters. Therefore, we present a new radial line-search approach in the next section, which investigates the whole multi-dimensional domain of the defined input parameters.

4.2. Radial Line-Search Approach

Assuming interval constrains for each response value $y_i$ according to Eq. 23, 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\}.$$

(25)

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 12: 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},$$

(26)

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

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

(27)

The direction vectors ${\mathbf{a}}_l$ can be chosen randomly or uniformly distributed as shown in Figure 12 for two input parameters. For each direction a line-search is performed in order to obtained 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 known from the reliability analysis [44,45].

In Figure 13 the initial radial search points and the obtained boundary points are shown for the 2D subspaces of the three parameter bilinear model by assuming an interval of $\Delta y_i=\pm 0.2N/m{m}^2$. The line-search procedure is performed for each subspace separately while keeping the third parameter constant at the reference values. 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.

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

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

(28)

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. 28 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.

In Figure 14 the initial search points are shown for the bilinear example on the 3D sphere by considering the parameter bounds of the 2D analysis. The figure indicates a very uniform distribution of the points. Additionally, the figure shows the final line-search points for these 500 directions, where each point has been obtained by 10 bisection steps along the specific direction. The figure shows a Delaunay triangulation of the convex hull of these points utilizing the MATLAB plotting library [48]. 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 du to the pure elastic behavior and the resulting perfect correlation of the stress-strain curves.

5. Application Example

5.1. Sensitivity Analysis and Deterministic Parameter Optimization

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

Six unknown parameters were 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}$ of the bi-linear softening law. In Table 1 the reference parameters according to [49] 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 was investigated by generating 200 samples within the defined bounds by using an improved Latin-Hypercube Sampling according to [42]. The corresponding force-displacement curves of these samples are given in Figure 17 additionally to the experimental values according to [49]. The total effect sensitivity indices of the model parameters were estimated for each force value on the load-displacement curve by using the Metamodel of Optimal Prognosis [33] approach as discussed in Section 2. 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 did not consider it in the following parameter calibration and uncertainty estimation. The other parameter show a certain influence whereby the Young’s modulus, the tensile strength and the fracture energy are most important.

The optimal parameter set was determined by a least-squares minimization using Eq. 1, whereby the deviation of the response values were 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 [39] was 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.

5.2. Probabilistic Uncertainty Estimation

In a next step the Markov estimator and the Bayesian updating approaches were applied to estimate the parameter uncertainties. The covariance matrix of the 13 measurement points was defined according to Eq. 16 as a diagonal matrix, whereby the constant standard deviation was 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 were obtained from the estimated covariance matrix by using Eq. 17. The partial derivatives required in Eq. 14 have been estimated using a central differences approach with an interval of $0.01\%$ of the initial parameter ranges. Thus, only 10 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, which is the fundamental assumption of the Markov estimator.

The Bayesian updating has been performed with the same diagonal measurement covariance matrix used for the Markov estimator. 50.000 samples have been 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 and 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. For the fracture energy and the shape parameters similar statistical properties have been obtained. Even the estimated parameter correlations are quite close to the Markov estimator results.

5.3. Radial Line-Search

Finally, the proposed radial line-search algorithm is applied for 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 and 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 [39] 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.

6. Discussion

In the presented paper probabilistic and interval approaches for the uncertainty quantification of calibrated simulation models have been investigated. The Markov estimator and the Bayesian updating require the knowledge of the measurement covariances to estimate the multi-variate parameter distribution functions. If the covariance is estimated by a diagonal matrix assuming independent measurement errors, the number of considered measurement points directly influences the estimated parameter uncertainty and the results might be difficult to interpret.

The presented interval optimization approach assumes only maximum and minimum values of the measurement points and determines the corresponding boundary of the parameter values. From the viewpoint of the author, these results are easier to interpret. The introduced novel radial line-search procedure will detect a given number of boundary points in an efficient manner. Fekete point sets are utilized for cases with more than two input parameters to ensure an uniform representation of the boundary. However, the computation requires a large number of model evaluations but can be computed in parallel for each independent search direction. Similar to Bayesian updating, surrogate models could be applied to increase the efficiency.

In all investigated methods in this paper, the modeling error is neglected so far. The consideration of model uncertainty and model bias would be the scope of future research.