A Bayesian statistical method for quantifying model form uncertainty and two model combination methods

Reliability Engineering and System Safety 129 (2014) 46–56

Contents lists available at ScienceDirect

Reliability Engineering and System Safety journal homepage: www.elsevier.com/locate/ress

A Bayesian statistical method for quantifying model form uncertainty and two model combination methods Inseok Park a,n, Ramana V. Grandhi b a b

Honda R&D Americas, Inc. Ohio Center, 21001 State Route 739, Raymond, OH 43067, USA Department of Mechanical and Materials Engineering, Wright State University, 3640 Colonel Glenn Hwy Dayton, OH 45435, USA

art ic l e i nf o

a b s t r a c t

Article history: Received 27 December 2010 Received in revised form 24 April 2014 Accepted 27 April 2014 Available online 5 May 2014

Apart from parametric uncertainty, model form uncertainty as well as prediction error may be involved in the analysis of engineering system. Model form uncertainty, inherently existing in selecting the best approximation from a model set cannot be ignored, especially when the predictions by competing models show significant differences. In this research, a methodology based on maximum likelihood estimation is presented to quantify model form uncertainty using the measured differences of experimental and model outcomes, and is compared with a fully Bayesian estimation to demonstrate its effectiveness. While a method called the adjustment factor approach is utilized to propagate model form uncertainty alone into the prediction of a system response, a method called model averaging is utilized to incorporate both model form uncertainty and prediction error into it. A numerical problem of concrete creep is used to demonstrate the processes for quantifying model form uncertainty and implementing the adjustment factor approach and model averaging. Finally, the presented methodology is applied to characterize the engineering benefits of a laser peening process. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Model form uncertainty Model probability Maximum likelihood estimation Bayesian inference Model combination

1. Introduction Structural analysis of complex physical phenomena is becoming more dependent on computer simulation as nonlinear modeling methods advance. A simulation model can vary depending on the underlying physics and engineering and the manner in which a mathematical model is converted into a simulation model. This implies that we may have two or more different simulation models to analyze an identical engineering system. In statistics, a variety of model selection criteria have been developed from different perspectives with the aim to select the best model from a plausible regression model set such as Akaike Information Criterion (AIC) [1], Bayesian Information Criterion (BIC) [2], Mallows' Cp [3] and minimum description length [4]. Generally, a statistical model selection process is implemented with the concept that the best model balance goodness of fit with simplicity (usually measured by counting the number of regression parameters) better than the other considered models. In general, uncertainty exists in the process of selecting the best model from a model set. Uncertainty involved in model selection, called model form uncertainty in this paper because of

n

Corresponding author. Tel: þ 1 937 877 8354. E-mail address: [email protected] (I. Park).

http://dx.doi.org/10.1016/j.ress.2014.04.023 0951-8320/& 2014 Elsevier Ltd. All rights reserved.

its existence in model form, is due to lack of confidence in selecting the best model. Ignoring model form uncertainty is problematic because it may lead to underestimating the variability of predictions or making erroneous predictions [5]. It has been argued that a simple way to account for model form uncertainty is by model combination which takes all the predictions by a model set into account (detailed in Section 3) [6]. Model combination produces predictions that incorporate the (epistemic) variation inherent in the statistical model selection process as well as the (aleatory and/or epistemic) variation conditional on each model. Model combination aims to predict unknown responses more reliably than each model in a set rather than better represent the physics of a real system or update predictions of each model given measured experimental data. The model selection criteria to select the best model from a plausible regression model set can be adapted to the task of model combination. Unlike the statistical regression models that fit observed experimental data, mathematical or simulation models used in the engineering field (specifically functional forms in the models) are fundamentally created based on scientific and engineering knowledge. Some of analytical models such as semi-empirical models have their functional forms derived from both theoretical knowledge and empirical data. Recently, an increasing number of papers on the model form uncertainty quantification have been published in the engineering field while there exists plenty of literature on the big topic of V &

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

V, which broadly covers the model form UQ [7,8]. Alvin et al. [9] used BMA to estimate the model form uncertainty in the frequency prediction of a component mounting bracket resulting from the use of three stochastic simulation models having the parametric uncertainty in elastic material modulus. The three simulation models had different levels of simplifying assumptions in their PDE forms, as well as different spatial meshes and different discrete solution variables. Model probability, which is assigned to each model to quantify model form uncertainty (the detailed description of model probability is given in Section 2.1), was simply assumed to be uniformly distributed over the considered simulation models. Zio and Apostolakis [10] used the adjustment factor approach—which is being described in Section 3.2—to estimate the model form uncertainty in the response predictions regarding the cumulative release of a radionuclide to water table given by six different models. The six models differ by some fundamental hypotheses on the groundwater flow and transport mechnisms. Model probabilities were evaluated based on expert opinions. Zhang and Mahadevan [11] estimated the failure probabilities for the butt welds of a steel bridge using two competing crack growth models (the Foreman and the Weertman crack growth models). They made a reliability analysis of fatigue life by averaging the estimated failure probabilities weighted by the probabilities of the crack models. The uncertainty in crack size measurement was quantified to evaluate model probabilities using Bayes’ theorem. Zouaoui and Wilson [12] used BMA to quantify the model form uncertainty in prediction of a system response (a message delay in a computer communication network). Although simulation models were used to predict a message delay, model probabilities were not assigned to the simulation models but to three different types of distributions that represent the uncertainty in an input variable (i.e. a message length). The distributions for a message length were assumed to be of exponential, normal and lognormal forms. McFarland and Bichon [13] used BMA to incorporate probability distribution model form uncertainty into the estimation of failure probability of a bistable MEMS device. As in the work of Zouaoui and Wilson, model probabilities were assigned to the three types of distributions (normal, lognormal and Weibull) that represent the uncertainty in an input variable (i.e. edge bias on beam widths). In the above-mentioned research, evaluating model probability did not rely on a statistical analysis of the degree of agreement between the physically observed data on system responses and the simulated model predictions of the data. In recent years, attention of researchers also in the engineering field has increasingly been drawn to the statistical evaluation of the discrepancies between the test and the simulated data for the model probability quantification. To evaluate model probability using the measured differences of physically observed and simulated data in a practical and effective way, Park et al. [14] developed a statistical approach (detailed in Section 2) based on the fundamental idea behind BIC. The variance in prediction errors of each model is estimated using the Maximum Likelihood Estimation (MLE), and the best point estimate of the variance plays a critical role in the quantification of model probability. However, the composite prediction that only incorporates model form uncertainty was shown to underestimate the uncertainty in response predictions. The present research accounts for both model form uncertainty and (unknown) prediction error involved in each model to obtain highly reliable prediction of a system response. The two types of uncertainties are merged into response prediction using a model combination technique called model averaging. The MLE based method for quantifying model form uncertainty is presented with a fully Bayesian estimation for doing it in Section 2. Two model combination methods, model averaging and the

47

adjustment factor approach, are described and compared in Section 3. In Section 4.1, the presented uncertainty quantification method and the model combination methods are illustrated with a numerical problem of concrete creep. Model form uncertainty and prediction errors associated with the FE analysis of a laser peening process are quantified and are incorporated into composite prediction in Section 4.2.

2. Quantification of model probability 2.1. Bayes' theorem for quantification of model probability Model probability is defined as the degree of belief that a model is the best approximation among a set of possibilities; here, the best approximating model is defined as the model that predicts system responses of interest, usually unknown, more accurately than the other models in a model set. Consider a set of models denoted by M1, M2, …, MK and experimental data D. Given experimental data D, Bayes' theorem presents a way to update prior probability of model Mk into posterior probability of Mk by PrðM k jDÞ ¼

PrðM k Þ  LðM k jDÞ K

; k ¼ 1; …; K

ð1Þ

∑ PrðM l Þ  LðM l jDÞ

l¼1

Pr(Mk) is prior probability of model Mk, the degree of belief that model Mk is the best approximation assessed prior to observing experimental data D. Zio and Apostolakis investigated the formal process of eliciting and interpreting expert judgments to quantify prior model probability Pr(Mk) [10]. The quantification of Pr(Mk) using a corpus of knowledge is arbitrary in nature because logically rigorous relations do not exist between the knowledge about a model set and prior model probability. Prior model probability Pr(Mk) is often given a uniform value to avoid the difficulty of numerically specifying prior knowledge. L(Mk|D) is called the likelihood of model Mk given experimental data D. L(Mk| D) is deeply related to the goodness of fit of model Mk to experimental data D relative to the other models in the model set. Since uniform prior model probability is assumed for this research, the only concern is the evaluation of model likelihood L (Mk|D). A methodology to evaluate model likelihood given experimental data is presented in Section 2.2. 2.2. Evaluation of model likelihood 2.2.1. Unknown prediction error Because a (mathematical or simulation) model is just an approximation to a real physical system, it unavoidably involves an error in its prediction of a response; here prediction error is the difference between model prediction and experimental data. A prediction error is unknown unless the corresponding response is measured from the considered physical system. Probability distribution is generally used to describe an unknown prediction error; more specifically, (aleatory and epistemic) uncertainty in prediction error is mathematically characterized by a probability distribution. In statistics, prediction errors are usually assumed to be normally distributed with zero mean and a constant variance [15]. It is often stated that all variations in observed experimental data that cannot be explained by the considered model are included in the error term [16]. The statistical theory for dealing with random prediction error is well-understood and allows for constructing easily interpretable statistical intervals for predictions. However, this way of describing prediction error does not accommodate any framework to discern between the sources causing prediction error.

48

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

To differentiate between the causes of prediction error, Kennedy and O'Hagan [17] first introduced a model discrepancy term, which quantifies the systematic discrepancy between the true system response and model prediction, in addition to an experimental error term (similar to a random error in regression model). The model discrepancy term empirically describes model bias due to the incompleteness in the science or engineering used to construct a model, imperfect numerical implementation and any inaccurate input parameters. The experimental error term is usually assumed to be an independent normal variable with zero mean. Model discrepancy and experimental error comprising prediction error are not discerned in this research because a separate incorporation of the model discrepancy term shifts predictions by a model from the initially predicted values (although this can reduce the biases involved in model predictions). Also, increasing the number of unknown parameters to separately include the discrepancy term may cause problem such as predictive ability (overfitting) and statistical identifiability under the circumstances there are generally quite limited amounts of experimental data available in the engineering applications. Not separating model discrepancy and experimental error term in mathematically specifying prediction error is thought to be desirable for addressing model form uncertainty quantification although this perspective may be inappropriate for other statistical inferences such as model updating and model validation. A common formulation to involve an unknown error into a prediction by each deterministic model is represented by y ¼ f k þ εk ; εk  Nð0; s2k Þ

ð2Þ

where y is an unknown system response of interest, and fk is the prediction of y by a deterministic model Mk. εk is the error involved in prediction fk that is conceptually understood as resulting from bias associated with model Mk and measurement error on experimental data. An unknown error εk is assumed to be a normal variable with zero mean and a constant variance. The main purpose of introducing this assumption is to facilitate the development of a practical method to quantify model form uncertainty given experimental data. Use of εk with zero mean does not shift model prediction fk and reflects the fact that each model claims that its prediction fk is the most probable value. However, the assumption of zero mean and constant variance on prediction error may be impracticable if model predictions are significantly biased against experimental data due to the effect of model inaccuracy rather than due to the randomness in experimental data, and the differences between experimental data and model predictions tend to keep increasing or decreasing. Prediction errors are assumed to be independent variables in this research due to insufficient experimental data to estimate the correlation among errors. The covariance matrix of prediction errors with non-zero covariance and non-constant variance elements should be estimated if the prediction errors are highly correlated, the homoscedasticity assumption is seriously violated, and experimental data is enough to evaluate the covariance elements.

2.2.2. Evaluation of model likelihood using observed experimental data Eq. (2), which represents the probabilistic prediction of a response y by each model Mk, is restated in a probability distribution form as shown in 1 ðy  f k Þ gðyjsk ; M k Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiexp  2s2k 2πsk 2

2

! ð3Þ

where sk is the unknown standard deviation of the error involved in model prediction fk. g (y|sk, Mk) is called prior predictive distribution of response y under model Mk. The likelihood function of sk for each model Mk given a single data is represented by ! ε2 1 ð4Þ Lðsk ; M k jdn Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiexp  kn2 ; εkn ¼ dn  f kn 2 sk 2πs2 k

where dn is a single point of observed experimental data D¼ {d1, …, dN}. f kn is the deterministic prediction of experimental data dn by each model Mk. εkn is considered to be a random sample from a normal distribution with zero mean and standard deviation sk. The assumption that prediction errors of model Mk are regarded as independent variables implies that experimental data are also independent variables. Because likelihood function L(sk, Mk|D) is the multiplication of L(sk, Mk|dn) over data D, L(sk, Mk|D) is represented by 0 1 N 2 !N2 ∑ ε B kn C N 1 B C expB  n ¼ 1 2 C ð5Þ Lðsk ; M k jDÞ ¼ ∏ Lðsk ; M k jdn Þ ¼ 2 @ A 2 πs 2 s n¼1 k k where N is the number of observed experimental data D. Model likelihood L(Mk|D) is expressed by marginal likelihood integral [5,6] Z Lðsk ; M k jDÞgðsk jM k Þdsk ð6Þ LðM k jDÞ ¼ where g(sk|Mk) is the prior distribution of unknown variable sk conditional on Mk. A normalized variant of BIC based on maximum likelihood estimation (MLE) [18] is implemented to evaluate model likelihood L(Mk|D). This method has the advantages that it is efficient (the marginal likelihood integral does not need to be directly solved), and the solution of the integral equation does not depend on how to specify prior distribution g(sk|Mk) (model likelihood L(Mk|D) is not a single number but a function of prior distribution g(sk|Mk)). The maximum likelihood estimator of sk is the value of unknown variable sk most supported by observed data D, which maximizes likelihood function L(sk, Mk|D). The logarithm of a likelihood function, called the log-likelihood function, is often used for the convenience of operation; the two functions have their maximum points at the same value. Taking the derivative of the logarithm of Eq. (5) with respect to sk and setting it equal to zero yields 2 3 N 2 ε ∑ kn 7 d d 6 6N 7 n¼1 ln Lðsk ; M k jDÞ ¼ 6 ðln 1  ln 2π  2 ln sk Þ  7 dsk d sk 4 2 2s2k 5 N

¼

N

sk

∑ ε2kn

þ n ¼ 13

sk

¼ 0:

Solving this equation for s2k gives N

s ¼ ^ 2k

∑ ε2kn

n¼1

N

ð7Þ

Putting the maximum likelihood estimator of sk shown above into the exponential term in Eq. (5) leads to !N2   1 N ð8Þ LðM k jDÞ C exp  2 ^ k2 2π s

Introducing the i.i.d. normal variables to describe prediction errors leads to Eq. (8) that is easy to implement; however, the validity of the i.i.d. assumption needs be looked into before using Eq. (8).

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

Given that uniform prior model probability is assumed, Eqs. (1) and (8) leads to ^ k ÞN=2 ð1=2π s 2

PrðM k jDÞ ¼

K

^ l ÞN=2 ∑ ð1=2π s

ð9Þ

2

49

If g(y|Mk,D) has mean E(y|Mk,D) and variance Var(y|Mk,D), the mean and variance of composite distribution g(y|D) are calculated by K

EðyjDÞ ¼ ∑ PrðM k jDÞEðyM k ; DjÞ

ð11Þ

k¼1

l¼1

For a very large dataset, posterior model probability Pr(Mk|D) evaluated based on MLE is fundamentally the same as it evaluated by directly solving the marginal likelihood equation. On the other hand, given a small amount of data, whether or not the direct integration and MLE methods provide an equal result largely depends on what prior distribution g(sk|Mk) is since the influence of g(sk|Mk) is significant compared to observed data. In chapter 4 and 5 where numerical problems are dealt with, the maximum likelihood based estimation of model probability Pr (Mk|D) is compared with the fully Bayesian estimation of it, which is made by numerically solving the marginal likelihood integral L (Mk|D) that involves a certain distribution function specified for g (sk|Mk). Two uninformative priors, uniform distribution and Jeffreys prior, which is invariant under reparameterization, are assumed for g(sk|Mk): g(sk|Mk) p 1, and g(sk|Mk) p 1/sk over 0 o sk o1 (improper priors), respectively.

K

K

k¼1

k¼1

VarðyjDÞ ¼ ∑ PrðMk jDÞVarðyM k ; DÞ þ ∑ PrðMk jDÞðEðyjMk ; DÞ  EðyjDÞÞ2

ð12Þ In Eq. (12), the first term, called within-model variance, represents the average degree of uncertainty in each model prediction. The second term, called between-model variance, represents the degree of uncertainty in response prediction resulting from model form uncertainty. In this research, g(y|D) in Eq. (10) is normal because g(y|Mk,D) is a normal distribution as shown in Eq. (3), and a linear combination of normal distributions is also normally distributed. Using Eqs. (11) and (12), the composite predictive distribution of response y, which incorporates both model form uncertainty and prediction error is expressed by gðyjDÞ ¼ N

K

K

K

k¼1

k¼1

k¼1

ð13Þ

3. Model combination Given that model form uncertainty is quantified by evaluating model probability, the model form uncertainty can be propagated into the prediction of a system response using a model combination method. In this paper, two model combination methods, called model averaging and the adjustment factor approach, are used to combine deterministic and probabilistic predictions from a set of models into a single composite prediction. Both model averaging and the adjustment factor approach are founded on the idea of averaging model predictions weighted by model probabilities. While model averaging is applicable in the case that probabilistic prediction of a system response is made by each model, the adjustment factor approach is applicable in the case that deterministic prediction of a response is made by each model. Model averaging can account for not only model form uncertainty but also other types of uncertainty in response prediction by each model. The adjustment factor approach accounts for the model form uncertainty involved in the deterministic response predictions from different models by accommodating normal and lognormal distribution forms. 3.1. Model averaging Model averaging [5,6] is the method that can incorporate model form uncertainty into composite response prediction given that each model prediction is represented by a probability distribution. 3.1.1. Integration of predictive distributions estimated by a model set The integration of the probabilistic predictions of response y by a model set is represented by K

gðyjDÞ ¼ ∑ PrðM k jDÞgðyjM k ; DÞ k¼1

!

^ k 2 þ ∑ PrðMk jDÞðf k EðyjDÞÞ2 ∑ PrðMk jDÞf k ; ∑ PrðM k jDÞs

ð10Þ

where posterior predictive distribution g(y|Mk, D) under model Mk represents the uncertainty involved in the prediction of response y by model Mk, evaluated given experimental data D. g(y|D) is called composite predictive distribution of response y given experimental data D. Posterior predictive distribution g(y|Mk,D) can be simply constructed by plugging the maximum likelihood ^ k of sk into Eq. (3). estimator s

In the case of specifying the prior distribution of sk, the posterior distribution of sk given test data D is estimated by applying Bayes' theorem as shown in gðsk jM k ; DÞ ¼ R

Lðsk ; M k jDÞgðsk jM k Þ Lðsk ; M k jDÞgðsk jM k Þdsk

ð14Þ

Given posterior distribution g(sk|Mk,D), the posterior predictive distribution of unknown response y under model Mk is obtained by marginalizing over g(sk|Mk,D) as expressed by Z  gðysk ; M k Þgðsk jM k ; DÞdsk ð15Þ gðyjM k ; DÞ ¼ Finally, the composite predictive distribution of y, g(y|D), and the mean and the variance of it can be computed using Eqs. (10)– (12). 3.2. Adjustment factor approach Mosleh and Apostolakis [19] suggest the adjustment factor approach to combine experts' estimates according to Bayes' theorem. The application of this approach was extended to the model form uncertainty problem [10]. The adjustment factor approach characterizes model form uncertainty as the uncertainty associated with the model given the highest model probability. Model form uncertainty is accounted for by an adjustment factor represented by a probability distribution. Depending on whether the assumption that the distribution representing the model form uncertainty in prediction of a system response is normal or log-normal, an additive or a multiplicative adjustment factor is used, respectively [10,14,20]. An additive adjustment factor is added to the prediction by the model given the highest model probability to incorporate model form uncertainty into a composite prediction. Similarly, a multiplicative adjustment factor multiplies the prediction by the model given the highest model probability. Because the multiplicative adjustment factor approach is not used for this research, the additive adjustment factor approach is only described below. When an additive adjustment factor is used, the prediction of a response y is represented by [10,19] n

y ¼ f þ Ea n

ð16Þ

50

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

where fn is the prediction of response y by the model assigned the highest model probability among the considered model set. Ena is an additive adjustment factor represented as a normal variable. Supposing that the predictions and probabilities of a set of models are known, the means and variances of both Ena and y are computed by K

n

EðEna Þ ¼ ∑ PrðM k jDÞðf k  f Þ k¼1

K

VarðEna Þ ¼ ∑ PrðM k jDÞðf k EðyjDÞÞ2 k¼1

n

ð17Þ

ð18Þ

Eðy=DÞ ¼ f þ EðEna Þ

ð19Þ

VarðyjDÞ ¼ VarðEna Þ

ð20Þ

The composite prediction of response y is also a normal variable with mean E(y|D) and variance Var(y|D). The bigger the absolute value of E(Ena ) is, the more the composite prediction is shifted from prediction fn as shown in Fig. 1. The bigger the variance of an adjustment factor Var(Ena ) is, the larger the degree of uncertainty involved in the prediction of response y that results from model form uncertainty). The mean E(y|D) represented by Eqs. (17) and (19) is identical to the mean of the composite prediction in Eq. (13) as can be seen from Eq. (17) being expanded and plugged into Eq. (19). The variance Var(y|D) expressed by Eqs. (18) and (20) is the same as the second term in the variance of the composite prediction in Eq. (13). An interval estimate of a system response is made to indicate the reliability of a composite prediction g(y|D). Considering that a composite prediction is represented by a normal distribution in both model averaging and the additive adjustment factor approach, a 95% confidence interval, which represents the interval that is supposed to include responses to be observed with a probability of 0.95, is established by h pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii ð21Þ EðyjDÞ  1:96 VarðyjDÞ; EðyjDÞ þ 1:96 VarðyjDÞ Because the i.i.d. assumption is involved in deriving mean E(y| D) and variance Var(y|D) of composite prediction g(y|D), using Eq. (21) is limited to the cases in which the assumption is valid.

4. Demonstration problems

methods are applied to the problem of simulating a laser peening process. 4.1. A numerical problem of concrete creep Creep is defined as the slow deformation a solid material undergoes under the influence of sustained stresses generated by external loads over an extended period of time. Creep in concrete is not fully understood because it involves many factors that influence the total amount of creep a specimen experiences. The factors considerably influencing the concrete creep include volume and mechanical properties of aggregate, types and contents of cement, water–cement ratio, member size, curing condition, temperature, relative humidity, age of concrete at loading, and stress intensity. 4.1.1. Mathematical models to predict concrete creep Several mathematical models have been suggested to estimate the amount of creep in concrete. For this problem, four empirical models recommended by different professional societies are used to predict the amount of creep: ACI 209 [21], AASHTO [22], CEB-FIP [23] and JSCE [24] creep models for design codes. The description of the four models is given in Appendix A. Although the four creep-prediction models can be used to predict the creep behavior of a concrete specimen, those models may contain considerable errors. This is because the effects of concrete ingredients and environmental conditions on creep behavior are too complicated to be fully understood, and a limited amount of creep-test results were used when deriving those empirical creep-prediction models. 4.1.2. Creep test and model predictions Experimental data on concrete creep available in the literature [25] are used as information to quantify the model probabilities. The creep coefficients measured from the creep test at several concrete ages are shown in Table 1. By putting the numerical values specifying the creep test (presented in Appendix B) into Eqs. (A.1)–(A.4), each creep-prediction model is represented by M 1 ðtÞ ¼ 1:924 

t 0:6 ; ACI  209 model 10 þ t 0:6

ð22Þ

M 2 ðtÞ ¼ 0:995 

t ; ASSHTO model 24:65 þ t

ð23Þ

 A closed-form numerical problem of concrete creep is used to demonstrate the presented model-form uncertainty quantification method and the two model combination methods. Next, the

M 3 ðtÞ ¼ 2:265 

t t þ 364:31

0:3 ; CEB  FIP model

  M 4 ðtÞ ¼ 1:541  1  expð  0:097t 0:6 Þ ; JSCE model

ð24Þ ð25Þ

where t is the age of concrete (days) between time of loading and each time of interest. The creep coefficients predicted by the equations above are plotted on a scale of 0–360 days in Fig. 2.

Fig. 1. Additive adjustment factor approach.

4.1.3. Quantification of model probability Prediction errors, assumed to be the four distinct normal variables with zero mean and unknown variances are incorporated into the four creep-prediction models shown in Eqs. (22)–(25). Using the differences between the measured and predicted values of creep coefficients shown in Table 1, the unknown variances of the four normal variables are estimated using Eq. (7). Using the calculated variance values and the number of the measured differences, nine in this problem, the likelihood of each model given the test data is computed using Eq. (8). The prior probabilities of the considered models, assumed to be uniform, are updated into the posterior model probabilities by Eq. (9). The calculated posterior probabilities are shown in Table 1. Also, the

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

51

Table 1 Measured creep coefficients at specific time steps along with predicted coefficients, and evaluated posterior model probabilities. Age of concrete after loading

Measured creep coefficient

Predicted creep coefficient ACI 209 model

AASHTO model

CEB-FIP model

JSCE model

1 day 7 days 14 days 21 days 28 days 90 days 180 days 270 days 360 days

0.259 0.544 0.639 0.727 0.796 1.054 1.189 1.297 1.410

0.175 0.468 0.630 0.737 0.817 1.151 1.333 1.428 1.489

0.039 0.220 0.360 0.458 0.529 0.781 0.875 0.912 0.931

0.386 0.688 0.842 0.946 1.026 1.394 1.625 1.753 1.836

0.142 0.413 0.580 0.698 0.788 1.177 1.368 1.446 1.485

Posterior model probability

MLE Jeffreys Uniform

9.06  10  1 9.04  10  1 8.80  10  1

6.48  10  6 6.67  10  6 2.42  10  5

8.46  10  6 8.64  10  6 3.04  10  5

9.41  10  2 9.62  10  2 1.20  10  1

Fig. 2. Time-dependent predictions of creep coefficients by four mathematical models. along with measured creep coefficients.

posterior probabilities of the four models are calculated by implementing the fully Bayesian inference with a uniform and a Jeffreys prior distribution specified for unknown standard deviation sk for each model Mk. As can be seen from Table 1, the MLE based method shows high accuracy in evaluating the posterior model probabilities in comparison with the fully Bayesian method while it saves a large amount of time in calculating it. The calculated posterior probabilities imply that the AIC 209 creep-prediction model is most likely to give the closest predictions of unobserved creep coefficients among the considered models. The AASHTO and the CEB-FIP creep-prediction models have considerably low probabilities compared to the AIC 209 and the JSCE models due to their poor performance on the measured data.

4.1.4. Application of two model combination methods The additive adjustment factor approach is implemented to adjust the prediction of a creep coefficient from the ACI-209 model —which has the highest posterior probability—using the predictions from the other three models. The mean and variance of the adjusted prediction of creep coefficient at every time step are

calculated using Eqs. (17)–(20). Using the calculated mean and variance of adjusted prediction of each creep coefficient, a 95% confidence interval is established at each time step using Eq. (21). The width of the confidence band, which is bounded by the two curves interpolating the end points of the confidence intervals at all the depths (Fig. 3), indicates the degree of model form uncertainty in the composite prediction. Because there is no true model among the considered models, most of the measured creep coefficients are positioned outside the confidence band; here, a true model is defined as the model that mirrors the real system under consideration and generates the same data as observed responses of the system given no measurement error involved. The uncertainties in the composite predictions of creep coefficients are underestimated because prediction errors are not incorporated into the composite predictions. Given the predictions by the four models including normally distributed errors, composite prediction of creep coefficient at each time step is estimated using Eq. (13). The width of a 95% confidence band established using Eq. (21) indicates the degree of uncertainty in the composite prediction due to uncertainties in both model form and prediction error (Fig. 4). A 95% confidence interval of creep coefficient at each time step is also estimated in

52

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

Fig. 3. A 95% confidence band of the composite prediction of creep coefficient incorporating model form uncertainty alone.

Fig. 4. 95% confidence bands of the composite predictions of creep coefficient incorporating both model form uncertainty and prediction errors.

the fully Bayesian framework. As shown in Fig. 4, the confidence intervals established using the MLE based method are a little narrower than those using the fully Bayesian method due to not accounting for the uncertainty involved in estimating the standard deviation sk of each model’s prediction errors. Since all the measured creep coefficients are included within the established confidence band, it can be inferred that incorporating both model form uncertainty and prediction errors makes the predictions of creep coefficients more reliable.

4.2. Finite element simulation for a laser peening process 4.2.1. Problem description Laser Peening (LP) is an advanced surface enhancement technique that has been shown to increase the fatigue life of metallic components. During the LP process, laser energy is converted into shock waves at the surface that induce compressive residual stresses. A detailed description of the LP process is found in the literature [26,27].

In simulating a LP process, accurate description of material behavior is a challenging task because of the high strain rates experienced by the material. During a LP process, the strain rates experienced by a material can reach as high as 106/s. In such high strain-rate processes, different material models are available to describe the elastic-plastic behavior. For this research, three material models are considered to describe the unknown material behavior: the Johnson–Cook (JC) model, the Zerilli–Armstrong (ZA) model, and the Khan–Huang–Liang (KHL) model. Each material model results in a different finite element (FE) simulation to predict the residual stress field induced by the LP process. The simulation requires an extensive computer effort due to the modeling of material behavior under high pressure shock waves with time marching numerical procedures. A schematic illustration of the LP FE model is shown in Fig. 5. The peak value of the pressure pulse inducing a residual stress field is 5.5 GPa in this research. The simulation results of the three FE models based on different material theories are shown along with measured experimental data (Fig. 6). The outcomes from each FE simulation are considered to contain trivial numerical errors. The significant

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

53

Fig. 5. Representative axi-symmetric FE mesh for LP simulation model.

differences between the three simulation results shown in Fig. 6 imply that model form uncertainty in the predictions of residual stresses might be considerable. The experimental data measured at ten points are used to estimate the probabilities of the three FE models.

Table 2 Prior and posterior probabilities of three FE models for axi-symmetric LP component. Model

Prior probability Posterior probability MLE 1

4.2.2. Quantification of model probability For this problem, the prior probabilities of the three FE models are assumed to be uniform (Table 2) because of the unavailability of information to quantify them. The measured differences between the experimental and model outcomes in Fig. 6 are used to evaluate the likelihoods of the three models by Eq. (8). Given the calculated model likelihoods, the equal prior model probabilities are updated into the posterior model probabilities using Eq. (1) (Table 2). The calculated posterior probability (1.706  10  6) of the FE model that includes the ZA material model (the ZA-based FE model) is significantly smaller than the posterior probabilities (6.318  10  1 and 3.682  10  1) of the two FE models that include the JC and the KHL material models (the JC- and the KHL-based FE models, respectively). To demonstrate the effectiveness of the MLE based method for model probability computation, the posterior model probabilities are also computed by implementing the fully Bayesian method. As shown in Table 2, the MLE based method estimates the posterior model probabilities with good accuracy even though the amount of measured data is small. 4.2.3. Composite prediction of a residual stress field Using model averaging, both model form uncertainty and prediction error are incorporated into the composite prediction of the residual stress field. Using Eqs. (13) and (21), a 95% confidence band is established to indicate the degree of uncertainty in the composite prediction of a residual stress at each depth due to the uncertainties in both model form and prediction error (Fig. 7). The fully Bayesian method is also utilized to establish the 95% confidence band. The MLE based method provides a little narrower confidence band than the two bands established using the fully Bayesian method as shown in Fig. 7 mainly because it does not take into account the uncertainty in estimating unknown standard deviation sk.

JC-based FE model 3.333  10 KHL-based FE model 3.333  10  1 ZA-based FE model 3.333  10  1

Jeffreys 1

Uniform 1

6.32  10 6.32  10 6.19  10  1 3.68  10  1 3.69  10  1 3.81  10  1 1.71  10  6 1.71  10  6 6.03  10  6

two uninformative priors for an unknown parameter to verify the accuracy of the MLE based method. Because the composite prediction of a response is made depending on a model set, the model-to-model variance in the prediction is incorporated into the composite prediction. The composite prediction is reliable since it incorporates prediction error of each model as well as model form uncertainty. A numerical problem of concrete creep is utilized to demonstrate the MLE based model form uncertainty quantification and the two model combination methods. The fully Bayesian inference method is also implemented to prove the effectiveness of the MLE based method. Finally, the presented method is applied to the FE analysis of a laser peening process. In each problem, the MLE based method is shown to provide fairly good approximations of the posterior model probabilities and the confidence intervals. Given deterministic predictions by a model set and measured experimental data, the presented method can be applied to any problem of model form uncertainty quantification regardless of the numbers of considered models and the experimental data.

Acknowledgments The authors acknowledge the support of this research work through the Contract FA8650-04-D-3446, DO # 25 sponsored by the Wright Patterson Air Force Base, Ohio.

Appendix A. Four creep-prediction models

5. Conclusion

ACI 209 creep-prediction model [21]

Model form uncertainty is quantified by evaluating the marginal model likelihood integral using MLE. Also, the marginal likelihood integral is directly solved in the Bayesian setting with

The creep-prediction model recommended by ACI Committee 209 uses an ultimate creep coefficient that may be corrected to account for a variety of factors influencing the magnitude of creep

54

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

Fig. 6. Residual stress comparisons between the predictions from three FE models and measured experimental data for axi-symmetric LP component.

Fig. 7. 95% confidence bands of the composite predictions of a residual stress field.

in concrete and also introduces a time-dependent function to express the growth in creep over time. Creep coefficient for a time of interest, defined as the ratio of creep strain for that time to initial elastic strain occurred by external load is represented by M 1 ðtÞ ¼ 2:35ðγ la  γ λ  γ vs  γ ψ  γ s  γ a Þ  where

t 0:6 10 þ t 0:6

ðA:1Þ

γla: loading age correction factor;

γ la ¼ 1:25t la  0:118 ðf or non  accelerated  cured concreteÞ where tla: loading age (days), γλ: relative humidity correction factor; γ λ ¼ 1:27 0:0067  RH (for RH440%) where RH: relative humidity (%), γvs: volume-to-surface area ratio correction factor;      γ vs ¼ 2=3 1 þ 1:13exp  0:0213  v=s

correction factor;

γ s ¼ 0:82 þ 0:00264s where s: observed slump (mm),

γa: air content correction factor;

γ a ¼ 0:46 þ0:09a Z1:0 where a: air content (%), t: length of time after loading (days). AASHTO creep-prediction model [22] The AASHTO (2007) creep-prediction model accounts for the volume-to-surface area ratio, the relative humidity, the concrete strength and the development of strength with time. Creep coefficient for each time step is determined by M 2 ðt Þ ¼ 1:9  ks  khc  kf  ktd  t i  0:118

ðA:2Þ

where v/s: volume-to-surface area ratio (mm), γψ: fine aggregate percentage correction factor;

where ks: factor accounting for the effect of the volume-to-surface area ratio;  ks ¼ 1:45  0:0051 v=s Z1:0;

γ ψ ¼ 0:88 þ 0:0024ψ

khc: relative humidity factor;

where

ψ: ratio of fine to total aggregate by weight (%), γs: slump

khc ¼ 1:56  0:008H

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

where H: relative humidity (%), kf: concrete strength factor; 35 kf ¼ ' 7 þ f ci where f'ci: compressive strength of concrete at 28 days (MPa), ktd: time development factor; ktd ¼

t 61  0:58f

'

ci þ t

;

t: age of concrete (days), defined as age of concrete between time of loading and time being considered for analysis of creep effects, ti: age of concrete when load is initially applied for accelerated curing (days), minus 6 days for moist-curing. CEB-FIP creep-prediction model [23]

55

m), W/C: water-to-cement ratio, RH: relative humidity (%), V: volume (mm3), S: surface area in contact with outside air (mm2), V/S: volume-to-surface area ratio (mm), Ec: elastic modulus of concrete at time of loading (MPa), t0, t' and t: age (days) of concrete at the beginning of drying, at the beginning of loading, and during loading, respectively; these factors should be corrected by

n 4000 t 0 ; t 0 and t ¼ ∑ Δt i exp 13:65  273 þ TðΔt i Þ=T 0 i¼1 where 1 1C.

Δti: number of days when the temperature is T (1C), T0:

Appendix B. Numerical values input in the four creep-prediction models

The mathematical model for creep prediction contained in the European design code, CEB-FIP Model Code 1990, takes into account the cement type and curing temperature that the ACI 209 and the AASHTO creep-prediction models do not account for. Creep coefficient for each time step is determined considering the fresh and hardened properties of a concrete mixture and environmental conditions as shown in  0:3 ðt t 0 Þ=t 1 M 3 ðtÞ ¼ ΦRH  β ðf cm Þ  β ðt 0 Þ  ðA:3Þ βH þ ðt  t 0 Þ=t 1

The dimensions, the mixture, the mechanical properties of the concrete specimen, and the conditions of the creep test are listed below:

where Ф(t, t0): creep coefficient, t: age of concrete at the moment considered (days), t0: age of concrete at time of loading (days), where

         

ΦRH ¼ 1 þ βðf cm Þ ¼ βðt 0 Þ ¼

1  RH=RH 0 0:46ðh=h0 Þ1=3 5:3

ðf cm =f cm0 Þ0:5

;

1 0:1 þ ðt 0 =t 1 Þ0:2 "



βH ¼ 150 1 þ 1:2

;

;

RH RH 0

18 #

h þ 250 r 1500 h0

where h: 2Ac/u (mm), fcm: mean compressive strength of concrete at 28 days (MPa), fcm0: 10 MPa, RH: relative humidity of the ambient environment (%), RH0: 100%, Ac: cross-sectional area (mm2), u: perimeter of the member in contact with the atmosphere (mm), h0: 100 mm, t1: 1 day. JSCE creep-prediction model [24] The creep-prediction model adopted in the JSCE Specification in Japan accounts for the unit cement and water contents, the water-cement ratio, the relative humidity, the volume-to-surface area ratio and the environmental temperature. Creep coefficient for each time step is estimated by h

i M 4 ðt Þ ¼ 1  exp  0:09ðt t 0 Þ0:6  ε0cr  Ec ðA:4Þ where

ε0cr ¼ ε0bc þ ε0dc ; ε0bc ¼ 1:5ðC þWÞ2:0 ðW=CÞ2:4 ðlog e t 0 Þ  0:67  10  9 ; ε0dc ¼ 4:5ðC þWÞ1:4 ðW=CÞ4:2 ðlog e ðV =S=10ÞÞ  2:2 ð1  RH=100Þ0:36 t 0  0:30  10  7 ; where C: unit cement content (kg/m3), W: unit water content (kg/

    

Volume of the specimen: 5.56  106 mm3 Surface of the specimen: 182,415 mm2 Volume-to-surface ratio: 30.48 Cross-sectional area: 182.41 mm2 Perimeter of the specimen in contact with the atmosphere: 478.78 mm Unit cement content: 379.7 kg/m3 Unit water content: 160.2 kg/m3 Water-to-cement ratio: 0.42 Ratio of fine to total aggregate: 0.362 Slump flow: 184.2 mm Air content: 3.5% Compressive strength at time of loading: 62.7 MPa Temperature: 22.8 1C Relative humidity: 50.0% Age of concrete at time of loading: 28 days

References [1] Akaike H. A new look at the statistical model identification. IEEE Trans Autom Control 1974;19(6):716–23. [2] Schwarz G. Estimating the dimension of a model. Ann Stat 1978;6(2):461–4. [3] Mallows CL. Some comments on Cp. Technometrics 1973;15(4):661–75. [4] Rissanen J. Modeling by the shortest data description. Automatica 1978;14 (5):465–71. [5] Draper D. Assessment and propagation of model uncertainty. J Royal Stat Soc Ser B 1995;57(1):45–97. [6] Raftery AE, Madigan D, Hoeting JA. Bayesian model averaging for linear regression models. J Am Stat Assoc 1997;92(437):179–91. [7] Li W, Chen W, Jiang Z, Lu Z, Liu Y. New validation metrics for models with multiple correlated responses. Reliab Eng Syst Saf 2014;127:1–11. [8] Yuan J, Ng S. A sequential approach for stochastic computer model calibration and prediction. Reliab Eng Syst Saf 2013;111:273–86. [9] Alvin KF, Oberkampf WL, Diegert KV, Rutherford BM.Uncertainty quantification in computational structural dynamics: a new paradigm for model validation. In: Proceedings of the 16th international modal analysis conference; 1998. p. 1191–8. [10] Zio E, Apostolakis G. Two methods for the structured assessment of model uncertainty by experts in performance assessments of radioactive waste repositories. Reliab Engi Syst Saf 1996;54(2–3):225–41. [11] Zhang R, Mahadevan S. Model uncertainty and Bayesian updating in reliability-based inspection. Struct Saf 2000;22(2):145–60. [12] Zouaoui F, Wilson JR.Accounting for input model and parameter uncertainty in simulation. In: Proceedings of the 2001 winter simulation conference; 2001. p. 290–9. [13] McFarland JM, Bichon BJ. Bayesian model averaging for reliability analysis with probability distribution model form uncertainty. In: Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC structures, structural dynamics and materials conference. CA: Palm Springs; 2009. [14] Park I, Amarchinta HK, Grandhi RV. A Bayesian approach for quantification of model uncertainty. Reliab Eng Syst Saf 2010;95(7):777–85.

56

I. Park, R.V. Grandhi / Reliability Engineering and System Safety 129 (2014) 46–56

[15] Gelman A, Hill J. Data analysis using regression and multilevel/hierarchical models. New York: Cambridge University Press; 2007. [16] Trucano TG, Swilera LP, Igusab T, Oberkampf WL, Pilch M. Calibration, validation, and sensitivity analysis: what's what. Reliab Eng Syst Saf 2006;91 (10–11):1331–57. [17] Kennedy MC, O'Hagan A. Bayesian calibration of computer models. J Royal Stat Soc. Ser B 2001;63(3):425–64. [18] Edwards AWF. Likelihood. Cambridge: Cambridge University Press; 1972. [19] Mosleh A, Apostolakis G. The assessment of probability distributions from expert opinions with an application to seismic fragility curves. Risk Anal 1986;6(4):447–61. [20] Riley ME, Grandhi RV. Quantification of model-form and predictive uncertainty for multi-physics simulation. Comput Struct 2011;89(11–12):1206–13. [21] ACI Committee 209, ACI 209R-92: Prediction of Creep, Shrinkage, and Temperature Effects in Concrete Structures (Reapproved 1997), Farmington Hills. Michigan: American Concrete Institute; 1997.

[22] AASHTO, AASHTO LRFD Bridge Design Specifications: Customary U.S. Units. 3rd ed. Washington D.C.: American Association of State Highway and Transportation Officials (AASHTO); 2007. [23] CEB, Creep and Shrinkage, In CEB-FIP Model Code 1990, Lausanne. Switzerland: Comite Euro-International du Beton (CEB); 1990. [24] Sakata K, Shimomura T. Recent progress in research on and code evaluation of concrete creep and shrinkage in Japan. J Adv Concr Technol 2004;2(2):133–40. [25] Kavanaugh B. Creep behavior of self-consolidating concrete [M.S. thesis]. Alabama: Auburn University; 2008. [26] Singh G, Grandhi RV. Mixed-variable optimization strategy employing multifidelity simulation and surrogate models. AIAA J 2010;48(1):215–23. [27] Amarchinta HK, Grandhi RV, Clauer AH, Langer K, Stargel D. Simulation of residual stress induced by a laser peening process through inverse optimization of material models. J Mater Process Technol 2010;210(14):1997–2006.