Structural Safety 32 (2010) 183–190

Contents lists available at ScienceDirect

Structural Safety
journal homepage: www.elsevier.com/locate/strusafe

Interval Monte Carlo methods for structural reliability

Hao Zhang a,*, Robert L. Mullen b,*, Rafi L. Muhanna c,*
a
School of Civil Engineering, University of Sydney, NSW 2006, Australia
b
Department of Civil Engineering, Case Western Reserve University, Cleveland, OH 44106, USA
c
School of Civil and Environmental Engineering, Georgia Institute of Technology, Savannah, GA 31407, USA

a r t i c l e i n f o a b s t r a c t

Article history: This paper considers structural reliability assessment when statistical parameters of distribution func-
Received 31 August 2009 tions can not be determined precisely due to epistemic uncertainty. Uncertainties in parameter estimates
Received in revised form 4 January 2010 are modeled by interval bounds constructed from confidence intervals. Reliability analysis needs to con-
Accepted 6 January 2010
sider families of distributions whose parameters are within the intervals. Consequently, the probability of
Available online 2 February 2010
failure will vary in an interval itself. To estimate the interval failure probability, an interval Monte Carlo
method has been developed which combines simulation process with the interval analysis. In this
Keywords:
method, epistemic uncertainty and aleatory uncertainty are propagated separately through finite ele-
Bayesian approach
Epistemic uncertainty
ment-based reliability analysis. Interval finite element method is utilized to model the ranges of struc-
Interval analysis tural responses accurately. Examples are presented to compare the interval estimates of limit state
Interval finite element probability obtained from the proposed method and the Bayesian approach.
Monte Carlo simulation Ó 2010 Elsevier Ltd. All rights reserved.
Parameter uncertainty

1. Introduction the scope of this work, the statistical uncertainty will be consid-
ered here. In particular, we are interested in reliability assessment
A major step in structural reliability analysis is the modeling when the parameters of distribution functions can not be deter-
and quantification of various sources of uncertainty. It is common mined precisely due to small sample size.
in engineering practice to distinguish between aleatory uncer- The Bayesian approach is routinely used to consider the uncer-
tainty and epistemic uncertainty [1,2]. Aleatory uncertainty is tainties associated with estimation of parameters of a probability
due to the inherent random nature of physical quantities (e.g., vari- distribution. The unknown parameters are assumed to be (Bayes-
abilities in yield strength of steel). Aleatory uncertainties are gen- ian) random variables [3]. The epistemic uncertainty and aleatory
erally modeled by random variables. In contrast to aleatory uncertainty are combined through the total probability theorem.
uncertainty, epistemic uncertainty is knowledge-based, and arises With the Bayesian approach, subjective judgments are required
from imperfect modeling, simplification and limited availability of to estimate the Bayesian random variables. The estimate of the
database. Possible sources of epistemic uncertainty include model Bayesian random variables can be improved when additional data
uncertainty and statistical uncertainty. Model uncertainty is re- become available. Before receiving more data, however, the Bayes-
lated to the discrepancy between real structural behavior and its ian approach remains a subjective representation of uncertainty.
simplified representation in mathematical models such as finite In this paper, incomplete knowledge of the distribution param-
element (FE) models. Statistical uncertainty is another important eters is modeled by the interval bounds constructed from confi-
source of epistemic uncertainty. The probability distribution to de- dence intervals. Based on the observational data, a confidence
scribe a random phenomenon is generally not precisely known. interval is established over which the parameter is located at a
The statistical parameters (e.g., mean and standard deviation) are specified level of confidence [3]. The epistemic uncertainties in
usually estimated by statistical inference from sampled observa- estimating the parameters are reflected in the widths of the inter-
tional data and a point estimator is used to approximate the ‘true’ vals. With the interval approach, reliability assessment needs to
parameter. Thus the distribution is itself subject to some uncer- consider families of distributions whose parameters are within
tainty. Statistical uncertainty may be significant if only a limited the intervals. One practical way to describe the ensemble of distri-
sample of data is available. While the model uncertainty is beyond butions is to specify its lower and upper bounds. The mathematical
frameworks using this methodology include Dempster-Shafer evi-
* Corresponding authors. Tel.: +61 2 93513923; fax: +61 2 93513343 (H. Zhang).
dence theory [4], random set theory [5] and probability boxes [6].
E-mail addresses: [email protected] (H. Zhang), [email protected] (R.L. Mul- The computation procedures are typically a combination of inter-
len), rafi[email protected] (R.L. Muhanna). val analysis and the Cartesian product method [7]. Lower and

0167-4730/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved.
doi:10.1016/j.strusafe.2010.01.001
184 H. Zhang et al. / Structural Safety 32 (2010) 183–190

upper bounds on the limit state probability are computed. Pen- the framework of the Bayesian approach, one can compute the fre-
metsa and Grandhi [8] considered random variables and interval quency distribution of pf , based on which a Bayesian confidence
variables simultaneously in structural reliability analysis. Function interval on pf can be estimated. Interval estimate of failure proba-
approximation was used to reduce the number of simulations. In bility can provide useful information to decision-makers about
[6,9], probability boxes and Dempster-Shafer structures were used the variability in reliability or risk [14]. When applying the Bayesian
to bound imprecisely specified probability distributions. The meth- approach, subjective judgment is needed to estimate fh ðhÞ. The esti-
od was applied to environmental risk assessment. Tonon et al. [10] mate of fh ðhÞ can be improved by using the Bayesian updating rule
considered the reliability analysis for an aircraft wing at the early when more data become available. Before receiving additional data,
stage of design process where the observational data is not however, the Bayesian approach remains a subjective representa-
point-valued but set-valued. Random set theory was employed to tion of epistemic uncertainty.
represent the envelope of all probability distributions compatible
with the available information. In [11], random sets were con-
2.1. Interval approach
structed from limited observational data by applying inequalities
of Tchebycheff to the sample mean. The method was utilized to
This paper adopts the confidence interval approach to represent
bound the statistics of the displacement response of a cantilever
the unknown parameters h. Let H denote the confidence intervals,
sheet pile wall. In [12], interval analysis was combined with
and h is a generic (arbitrary) element h 2 H. Under this assumption
first-order reliability method (FORM). The unknown means and
one needs to consider families of distributions whose parameters
standard deviations of random variables were modeled as interval
are in the intervals. Conceivably, the probability of failure pf will
numbers. Interval analysis was applied to the closed-form solu-
not be unique and vary in an interval. We are interested in estimat-
tions of FORM, and the bounds on structural system reliability in-
ing a lower bound and an upper bound of pf .
dex were evaluated.
A visualization of all possible distributions with h 2 H can be
Despite the research progress, the computing effort, especially
obtained by means of upper and lower distribution functions. Let
when Cartesian product method is used, is a barrier to the practical
FðxÞ denote the cumulative distribution function (CDF) for the ran-
application of non-traditional uncertainty models. The issue of
dom variable X. For every x, an interval ½FðxÞ; FðxÞ generally can be
computational cost becomes more serious when the reliability
readily found to bound the possible values of FðxÞ, i.e.,
analysis is FE-based, i.e., the structural responses are obtained
FðxÞ 6 FðxÞ 6 FðxÞ, for h 2 H. Such a pair of two CDFs FðxÞ and
through FE analyses.
FðxÞ construct a so-called probability box or probability bounds [6].
This paper proposes an interval Monte Carlo method to propa-
Fig. 1 shows the probability box for a normal distribution with
gate interval parameters through FE-based reliability assessment.
an interval mean of [2.0, 3.0] and a standard deviation of 0.5. In this
The FE portion of the analysis is carried out by the authors’ recently
simple example, it is easy to verify that FðxÞ is the CDF of the nor-
developed interval finite element method (FEM). An interval esti-
mal variable with a mean of 3 and FðxÞ is the one with a mean of 2.
mate on the failure probability is computed. Two plane structures
Probability box represents a general framework to represent
are analyzed to illustrate the proposed method. The interval failure
imprecisely specified distributions. It can represent not only distri-
probabilities obtained from the proposed method and the Bayesian
butions with unknown parameters, but also distributions with un-
approach are compared through the examples.
known type or even unknown dependencies [6].
Cartesian product method is routinely used for computing with
2. Reliability analysis under parameter uncertainty probability boxes [6,7]. In this method, a probability box is discret-
ized into a list of pairs fðA1 ; m1 Þ; . . . ; ðAi ; mi Þ; . . .g, in which Ai are
The basic reliability problem is defined by the multiple intervals and mi are their associated probability masses. Ai can be
integration termed as focal elements and mi can be viewed as the probability
Z Z that Ai is the range of x [17]. Thus a probability box is analogous
pf ¼ PðGðXÞ 6 0Þ ¼  fX ðxÞdx: ð1Þ to a discrete probability distribution except that the probability
GðXÞ60
mass is assigned to an interval rather than to a precise point. Dif-
Here pf represents the probability of failure of the structure. ferent discretization methods have been proposed, such as the
X ¼ ðX 1 ; . . . ; X n ÞT is the n-dimensional vector of the basic random Outer Discretization Method and the Averaging Discretization
variables representing uncertain quantities such as applied loads, Method [16]. The two discretization methods are graphically dem-
material strength and stiffness. fX ðxÞ is the joint probability density
function for X. GðXÞ is the limit state function and GðXÞ 6 0 defines
the failure state.
Let h denote the (unknown) statistical parameters of the distri-
bution function fX ðxÞ. In the Bayesian approach h are modeled to be
random variables, thus fX ðxÞ becomes a conditional distribution
function fXjh ðxjhÞ. Clearly the presence of h implies that the proba-
bility pf is random itself. The expectation of the conditional failure
CDF

probability is often computed to characterize the total uncertainty

[1,13]
Z Z Z F
p~f ¼ pf ðhÞfh ðhÞdh ¼  fXjh ðxjhÞfh ðhÞdxdh ð2Þ
h GðX;hÞ60

where fh ðhÞ is the joint probability distribution function of the

parameters h. However, there is a perception that the mean (or
other point estimate) of pf does not fully characterize the epistemic
x
uncertainty in the failure probability [14,15]. Alternatively, episte-
mic and aleatory uncertainties can be propagated through reliabil- Fig. 1. A probability box defined by a normal distribution with a mean of [2.0, 3.0]
ity analysis separately to obtain an interval estimate of pf . Within and a standard deviation of 0.5.
 H. Zhang et al. / Structural Safety 32 (2010) 183–190 185

onstrated in Fig. 2a and b. In practice, unbounded distributions are The basic Monte Carlo simulation formula can be extended to
truncated to a finite range. Clearly, the accuracy of the approxima- the case when fX ðxÞ is a probability box with h 2 H. When h varies
tion depends on how fine the discretization is [7]. in intervals, the randomly simulated basic variables x ^ j vary in
Consider the limit state function GðXÞ in Eq. (1). Suppose that intervals accordingly. The limit state function Gðx ^ j Þ becomes a
the basic variables X ¼ ðX 1 ; . . . ; X n ÞT are represented by n probabil- function of h as well, i.e., Gðx^j ; hÞ. If the minimum and the maxi-
ity boxes. If the basic variables are statistically independent, the mum values of Gðx ^ j ; hÞ can be determined
random relation pertaining to ðX 1 ; . . . ; X n Þ can be defined by a prob-
^j ; hÞÞ 6 Gðx
MinðGðx ^j ; hÞ 6 MaxðGðx
^j ; hÞÞ; for h 2 H ð6Þ
ability box fðR1 ; q1 Þ; . . . ; ðRl ; ql Þg on the Cartesian product of the fo-
cal elements of X 1      X n [17,18]. The lower and upper bound then
on the probability of GðXÞ 6 0 can be evaluated as [4]
^ j ; hÞÞ 6 0 6 I½Gðx
I½MaxðGðx ^j ; hÞ 6 0 6 I½MinðGðx
^ j ; hÞÞ 6 0: ð7Þ
X
f ¼
p qi ð3Þ
Applying Eq. (7) in (5) gives
Ri :0PinfðGðRi ÞÞ
X
pf ¼ qi ð4Þ 1 XN
1 XN
^ j ; hÞÞ 6 0 6
I½MaxðGðx ^ j ; hÞ 6 0
I½Gðx
Ri :0PsupðGðRi ÞÞ N j¼1 N j¼1
in which inf( ) and sup( ) denote infimum and supremum of the 1X N
function, respectively. p f and pf represent an upper and lower 6 ^ j ; hÞÞ 6 0:
I½MinðGðx ð8Þ
N j¼1
bound for pf respectively.
The above computation procedure requires that the image of Thus, Eq. (8) provides an interval estimate for pf
every focal element Ri through the limit state function Gð Þ be cal-
culated. If each probability box X i has k focal elements, the total 1X N
pf  ^j ; hÞÞ 6 0;
I½MaxðGðx
number of focal elements from the Cartesian products N j¼1
n
X 1      X n is k . For realistic engineering problems with large ð9Þ
n
number of n and/or k, the computing effort of performing k struc- 1X N
f 
p ^ j ; hÞÞ 6 0;
I½MinðGðx for h 2 H:
tural analyses is prohibitive. Williamson [7] introduced a conden- N j¼1
sation strategy of constructing coarser discretization for
2
probability boxes to reduce the calculation number to ðn  1Þk .
However, there is a trade off between computational cost and 3.2. Computational aspects
accuracy of results. To overcome the difficulty associated with
the Cartesian product method, an interval Monte Carlo simulation The first step in the implementation of interval Monte Carlo
procedure has been developed. simulation is the generation of intervals in accordance with the
prescribed probability boxes. The inverse transform method [3] is
often used to generate random numbers. Consider a random vari-
3. Interval Monte Carlo simulation
able X with CDF FðxÞ. If ðu1 ; u2 ; . . . ; um Þ is a set of values from the
standard uniform variate, then the set of values
3.1. Basic formula
xi ¼ F 1
X ðui Þ; i ¼ 1; 2; . . . ; m ð10Þ
In Monte Carlo simulation, the probability of failure is approxi-
mated as [1] will have the desired CDF FðxÞ. The inverse transform method can
be extended to perform random sampling from a probability box.
1 XN
Suppose that an imprecise CDF FðxÞ is bounded by FðxÞ and FðxÞ,
pf  ^ j Þ 6 0
I½Gðx ð5Þ
N j¼1 as shown in Fig. 3. For each ui in Eq. (10), two random numbers
are generated
where N is the total number of simulations conducted, x ^j represents
xi ¼ F 1 ðui Þ; xi ¼ F 1 ðui Þ: ð11Þ
the jth randomly simulated vector of basic variables, and I[ ] is the
indicator function, having the value 1 if [ ] is ‘true’ and the value 0 if Such a pair of xi and xi form an interval ½xi ; xi  which contains all pos-
[ ] is ‘false’. sible simulated numbers from the ensemble of distributions for a
CDF

CDF

x x

Fig. 2. Discretization of a probability box. (a) Outer Discretization Method; (b) Averaging Discretization Method [16].
186 H. Zhang et al. / Structural Safety 32 (2010) 183–190

connections between elements, avoiding element coupling. The

penalty method is then employed to impose constraints to recover
the connections between elements, and to ensure the compatibility
of the displacements. The system equation in the interval FEM
takes the following form
ðK þ QÞu ¼ p ð14Þ
where K is the interval stiffness matrix, u is the interval displace-
ment vector, p is the interval load vector, and Q is the deterministic
penalty matrix. Eq. (14) can be transformed into a fixed point
equation
Rp  RðK þ Q Þu0 þ ðI  RðK þ QÞÞu ¼ u ð15Þ
in which R is a nonsingular preconditioning matrix, and u0 is an
approximate deterministic solution. It can be verified that
Fig. 3. Generation of random number from distribution with interval parameters.
u ¼ u þ u0 . Based on Eq. (15), interval fixed point iteration is con-
structed [33,36]
particular ui . The method is graphically demonstrated in Fig. 3 for
Z þ Cu ðlÞ ¼ u ðlþ1Þ ð16Þ
one-dimensional case.
with Z ¼ Rp  RðK þ Q Þu0 , and the iterative matrix C ¼ IRðK þ Q Þ.
3.3. Computing the ranges of structure responses – interval FEM The iteration converges when

The computational effort of the interval Monte Carlo simulation u ðlþ1Þ # u ðlÞ : ð17Þ
is contingent on the efficiency of computing the range (max. and Then u ðlþ1Þ þ u0 guarantees to contain the exact solution set of Eq.
min.) of structural responses through FE analyses when the simu- (14). The original interval fixed point iteration implicitly assumes
lated basic variables vary in intervals. This task can be performed that the coefficients of K vary independently between their bounds.
by using the interval FEM. A variety of solution techniques have This assumption is not valid for the system equations that arise in
been proposed for the interval FEM, including the combinatorial the interval FEM. Special formulation has to be developed to remove
method [19,20], perturbation method [21,22], sensitivity analysis coefficient-dependence in the algorithm. By using the EBE tech-
method [23], optimization method [24,25], and interval analysis nique, it is possible to decompose the interval stiffness matrix K
method [26–28]. In this paper, the interval FE analysis is formu- into two parts
lated as an interval analysis problem. The interval analysis and
interval FE formulation is briefly introduced below. Further details K ¼ SD ð18Þ
are provided in the authors’ previous work [29–31]. in which S is a deterministic matrix and D is an interval diagonal
Interval analysis concerns the numerical computations involv- matrix whose diagonal entries are the interval variables associated
ing interval numbers. The four elementary operations of real arith- with each element (e.g., interval modulus of elasticity). The term Z
metic, namely addition (+), subtraction (), multiplication () and in Eq. (16) can then be reintroduced as
division () can be extended to intervals. Operations
 2 fþ; ; ; g over interval numbers x and y are defined by the Z ¼ Rp  RðK þ QÞu0 ¼ Rp  RQu0  RSDu0
general rule [32,33] ¼ Rp  RQu0  RSMd: ð19Þ
x  y ¼ ½minðx  yÞ; maxðx  yÞ for  2 fþ; ; ; g ð12Þ
It must be noted that in Eq. (19) Du0 is introduced as Md in which M
in which x and y denote generic elements x 2 x and y 2 y. Software is a deterministic matrix and d is an interval vector [26]. The com-
and hardware support for interval computation are available (e.g., ponents of d are the diagonal entries of D with the difference that
[34,35]). every interval variable occurs only once in d. This treatment elimi-
For a real-valued function f ðx1 ; . . . ; xn Þ, the interval extension of nates the coefficient-dependence in Z, which is critical for obtaining
f( ) is obtained by replacing each real variable xi by an interval var- a tight bound.
iable xi and each real operation by its corresponding interval arith- The interval fixed point iteration converges if and only if
metic operation. From the fundamental property of inclusion qðjCjÞ < 1 [37], where qðjCjÞ is the spectral radius of the absolute
isotonicity [32], the range of the function f ðx1 ; . . . ; xn Þ can be rigor- value of the iterative matrix C. To achieve a small qðjCjÞ, the choice
ously bounded by its interval extension function R ¼ ðQ þ SÞ1 is made such that

ff ðx1 ; . . . ; xn Þ j x1 2 x1 ; . . . ; xn 2 xn g # f ðx1 ; . . . ; xn Þ: ð13Þ C ¼ I  RðQ þ SDÞ ¼ I  RðQ þ SÞ  RSðD  IÞ ¼ RSðD  IÞ: ð20Þ

Eq. (13) indicates that f ðx1 ; . . . ; xn Þ contains the range of f ðx1 ; . . . ; xn Þ Numerical tests have shown that fast convergence (within 10 iter-
for all xi 2 xi . ations) generally can be achieved by using the above modified iter-
A natural idea to implement interval FEM is to apply the inter- ative algorithm. The developed linear elastic interval FEM has been
val extension to the deterministic FE formulation. Unfortunately, successfully applied to plane frame structures, as well as continuum
such a naïve use of interval analysis in FEM yields meaningless mechanics problems [29–31]. The structural responses can be accu-
and overly wide results [27,28]. The reason is that in interval arith- rately and efficiently computed.
metic each occurrence of an interval variable is treated as a differ-
ent, independent variable. It is critical to the formulation of the 4. Example 1: truss structure
interval FEM to identify the dependence between the interval vari-
ables and prevent the widening of results. In this paper, an ele- Fig. 4 shows a planar truss. The serviceability limit state of
ment-by-element (EBE) technique is utilized for element deflection is considered. The deflection limit at the midspan is
assembly [27,30]. The elements are detached so that there are no set to be 7.5 cm. Linear elastic analyses were performed.
 H. Zhang et al. / Structural Safety 32 (2010) 183–190 187

Table 2
Bayesian confidence intervals for pf , Case 1 (truss in Fig. 4).

Confidence 90% 95% 99%

pf [0.08%, 0.20%] [0.063%, 0.22%] [0.014%, 0.27%]

normal distribution with a mean of l and a standard deviation of

r. Assuming a uniform prior, for samples of size 30 the mean
cross-sectional p ffiffiffiffiffiffi of A1  A6 and A7 
area A15 are modeled by
pffiffiffiffiffiffi
Nð10:32; 0:516= 30Þ and Nð6:45; 0:323= 30Þ respectively [3].
Fig. 4. Truss structure.
Sample statistics for the loads are used in this case.
The distribution of pf was numerically determined through a
Table 1 double-loop simulation procedure. In the outer loop, the means
Sample statistics for the basic random variables (truss in Fig. 4). of cross-sectional areas were sampled. With the simulated mean
Variables Sample mean Sample standard deviation No. samples
values, the system failure probability was computed in the inner
loop. 10,000 simulations were made for the mean cross-sectional
A1  A6 (cm2) 10.32 0.516 30
A7  A15 (cm2) 6.45 0.323 30
areas. The relative frequency distribution of pf is shown in Fig. 5
Ln P 1 4.483 0.09975 20 with a mean of 0.14% and a COV (coefficient of variation) of 0.3.
Ln P 2 5.582 0.09975 20 It can be seen that the distribution of pf is skewed. A lognormal
Ln P 3 4.483 0.09975 20 is a good fit to the simulated data of pf . Based on the simulated dis-
Unit of P: kN. tribution of pf , one can determine the Bayesian confidence interval
for pf corresponding to a specified probability. For example, one
The cross-sectional areas for the 15 members and the three can state the probability that pf will be in the interval from
loads are identified as the basic random variables. All the 18 ran- 0.08% to 0.2% is 0.9. Table 2 summarizes the intervals for pf with
dom variables are assumed to be mutually statistically indepen- different confidence. The Bayesian confidence intervals were
dent. Assume that based on experience, the cross-sectional areas determined such that the central point of the interval equals to
can be modeled by normals, and the loads modeled by lognormals. the mean.
Suppose the statistics for the random variables were estimated
from limited samples of data. Table 1 gives the available sample 4.1.2. Interval approach
statistics for the cross-sectional areas and the logarithm of the Table 3 summarizes the confidence intervals for the mean
loads (Ln P). The Young’s modulus is assumed deterministic cross-sectional areas at different confidence levels. The random
(200 GPa). The limit state probability is found to be 0.14% if the cross-sectional areas with interval means are bounded by probabil-
sample statistics are used. We now take into account the effects ity boxes. For each confidence level, 500,000 interval Monte Carlo
of parameter uncertainty on reliability assessment. Two cases are simulations were performed with randomly generated intervals
considered: (1) the uncertain means for the cross-sectional areas, for the cross-sectional areas according to their probability bounds.
and (2) the uncertain means for the (logarithm of) loads. The inter- Interval FEM was utilized to find the minimum and the maximum
val Monte Carlo method and the Bayesian approach are used to ob- of the deflection at the midspan. The last row of Table 3 gives the
tain interval estimates for the failure probability. The results from interval estimates of pf , corresponding to the mean cross-sectional
the two approaches are compared. areas at a particular confidence level. For example, if the 95% con-
fidence intervals are used to bound the mean cross-sectional areas,
the failure probability is found to be between 0.048% and 0.30%. As
4.1. Case 1 expected, the width of the interval estimate for pf becomes wider
when a higher confidence level is used for the unknown
4.1.1. Bayesian approach parameters.
Assume that the standard deviations of the cross-sectional areas The computational efficiency of the proposed interval Monte
are equal to those obtained from the samples. Let Nðl; rÞ denote a Carlo simulation is evident in this example when compared with

0.04

0.035
mean = 0.0014
COV = 0.3
0.03
Relative frequency

0.025

0.02

0.015

0.01

0.005

0
0.0004 0.0010 0.0016 0.0022 0.0028 0.0034
Probability of failure
Fig. 5. Relative frequency distribution for pf obtained from the Bayesian approach, Case 1 (truss in Fig. 4).
188 H. Zhang et al. / Structural Safety 32 (2010) 183–190

Table 3 Table 4
Interval estimates for pf obtained from interval Monte Carlo simulations, Case 1 (truss Interval estimates for pf obtained from the interval Monte Carlo simulations, Case 2
in Fig. 4). (truss in Fig. 4).

Variables Confidence level for the mean cross-sectional area Variables Confidence level for ki
90% 95% 99% 90% 95% 99%
A1  A6 (cm2) [10.17, 10.48] [10.14, 10.51] [10.08, 10.56] k1 ; k3 [4.4465, 4.5199] [4.4395, 4.5269] [4.4258, 4.5407]
A7  A15 (cm2) [6.35, 6.55] [6.34, 6.57] [6.30, 6.60] k2 [5.5452, 5.6186] [5.5381, 5.6256] [5.5244, 5.6393]
pf [0.062%, 0.26%] [0.048%, 0.30%] [0.031%, 0.43%] pf [0.025%, 0.67%] [0.019%, 0.72%] [0.01%, 1.18%]

interval pf changes from [0.019%, 0.72%] to [0.01%, 1.18%] when the

the Cartesian product method. The Cartesian product method
confidence level for ki increases from 95% to 99%. The high upper
would require 1015 analyses if each of the 15 probability boxes is
bound of the interval pf corresponds to the long upper tail of the
discretized into 10 focal elements. Such a computing effort makes
distribution of pf obtained in the Bayesian approach.
the Cartesian product method impractical in this example.

4.3. Comparison of Bayesian approach and interval approach

4.2. Case 2
In this example the Bayesian approach and the interval ap-
4.2.1. Bayesian approach proach both imply that the mean values of (logarithm of) loads
The natural logarithm of loads (Ln Pi ; i ¼ 1; 2; 3) are normal dis- and member cross-sectional areas have significant effects on the
tributions. Assume that the standard deviations for Ln P i are equal structural reliability. The system reliability is particularly sensitive
to their sample standard deviations. Let ki denote the mean of Ln P i . to the means of (logarithm of) loads. Both approaches suggest that
According to the Bayesian approach, with uniform priors the pos-
pffiffiffiffiffiffi additional data, particularly for loads, should be collected if more
terior distributions of p and k3 are Nð4:483; 0:09975= 20Þ, and
k1ffiffiffiffiffiffi confidence in the reliability estimate is needed.
k2 is Nð5:582; 0:09975= 20Þ. 10,000 samples were drawn for each Although the Bayesian approach and the interval approach pro-
ki . The system failure probability associated with the sampled ki vide comparable interval estimates for the failure probability,
was computed and the frequency distribution for pf is shown in these two approaches are conceptually different. The results have
Fig. 6. It is evident that the frequency distribution is broad and different meanings. The interval estimate of pf obtained from the
skewed to the right with a long upper tail. The mean value of pf Bayesian approach has an explicit statistical implication. Compar-
is 0.17% with a very high COV of 0.76. The long upper tail and high ing Tables 2 and 3 indicates that the interval failure probabilities
COV of pf suggests that the parameters ki have a significant effect obtained from the interval approach tend to be wider than those
on the failure reliability. The results shown in Fig. 6 also suggest of the Bayesian approach. This is mainly due to the fact that prior
that the system failure probability is between 0.017% and 0.34% information on the unknown parameters was assumed and incor-
with 90% confidence. porated in the Bayesian approach. On the other hand, the interval
approach only uses the information of the bounds on unknown
4.2.2. Interval approach parameters, thus does not provide a statistical statement on the
Table 4 presents the confidence intervals on ki and the corre- resulting interval failure probability. Although some subjective
sponding interval estimates on the failure probability. 500,000 judgment is still needed in the interval approach, such as selecting
interval Monte Carlo simulations were performed. Comparing Ta- appropriate confidence levels for the unknown parameters, it gen-
bles 3 and 4 shows that the widths of the interval failure probabil- erally requires less subjective information than the Bayesian ap-
ities are noticeably wider in Case 2 than in Case 1. This suggests proach. The authors make no judgment as to the relative validity
that the failure probability is more sensitive to the errors in ki than of the Bayesian approach and the interval approach. The purpose
the mean cross-sectional areas. As evidenced in Table 4, the inter- here is to demonstrate the feasibility of propagating interval infor-
val failure probability only becomes slightly wider when the confi- mation on unknown parameters through reliability assessment by
dence level for ki increases from 90% to 95%. On the other hand, the using the proposed interval Monte Carlo methods.

0.06

0.05 mean = 0.0017

COV = 0.76
Relative frequency

0.04

0.03

0.02

0.01

0
0.000 0.002 0.004 0.006 0.008 0.010
Probability of failure
Fig. 6. Relative frequency histogram for pf obtained from the Bayesian approach, Case 2 (truss in Fig. 4).
 H. Zhang et al. / Structural Safety 32 (2010) 183–190 189

between column-to-column, and between beam-to-beam. No cor-

relation exists between column-to-beam.
The serviceability limit state considered is H=500 (21 mm) for
roof drift, where H is the total height of the frame. Linear elastic
analyses were performed. The bounds of the cumulative frequency
distribution of the roof drift were obtained after 10,000 interval
Monte Carlo simulations, and is as shown in Fig. 8. The probability
of failure is found to be between 0.11% and 1.29%.

6. Conclusion

An interval Monte Carlo method has been developed for reli-

ability assessment under parameter uncertainties represented by
confidence intervals. The interval information of unknown param-
Fig. 7. Two-bay two-story frame. eters and the inherent uncertainties are propagated separately
through reliability analysis. Interval FEM is utilized to model the
ranges of structural responses accurately and efficiently. Interval
5. Example 2: frame structure estimates of failure probability are computed which can provide
a statement of confidence in the results of the reliability estimate.
The second example demonstrates the application of the devel- A wide interval pf implies that epistemic uncertainties are large,
oped method on frame structures. A two-bay two-story steel frame thus additional data should be collected. The developed method
shown in Fig. 7 is considered. The frame geometry and member can also be used to study the sensitivities of failure probability
sizes are based on those of [38]. In the figure the column is denoted with respect to the changes in distribution parameters.
as ‘C’ and the beam as ‘B’. Subscripts indicate member number. The A truss structure and a steel frame have been analyzed to illus-
frame is subjected to lateral loads at each floor and vertical loads at trate the proposed method. In the truss example, the results pro-
the top. The loads are assumed deterministic. The Young’s modulus duced by using the proposed method and the Bayesian approach
is 200 GPa. The sectional properties (moment of inertia I and cross- were compared. Both approaches can obtain interval estimates
sectional area A) of each member are considered as basic random on the limit state probability and provide decision-makers and
variables. Suppose only interval bound information is available designers useful information about the variabilities in reliability
for the means of the random variables. Table 5 summarizes the estimates. The interval pf obtained from the proposed method
interval means and the standard deviations for the basic random tends to be wider than that from the Bayesian approach. This is
variables. All random variables are assumed normals. For simplic- to be expected since the interval approach only uses the informa-
ity, it is assumed that the moment of inertia and cross-sectional tion of the confidence bounds on the unknown parameters, while
area are statistically independent. Perfectly correlation is assumed in the Bayesian approach prior information is assumed and incor-
porated in analysis.
The introduced interval Monte Carlo method is based on direct
Table 5 sampling. Future work can explore the possibility of extending the
Random sectional properties for the frame of Fig. 7. interval analysis with more advanced sampling methods such as
the importance sampling to improve the computational efficiency.
Member lI (cm4) rI (cm4) lA (cm2) rA (cm2)
C1 [5327.2, 5494.9] 270.55 [35.4,36.5] 1.80
C2 [62696.6, 64670.2] 3184.2 [246.4,254.2] 12.5
References
C4 [2204.6, 2274.0] 111.97 [22.5,23.2] 1.14
C3, C5, C6 [50813.0,52412.4] 2580.6 [203.3, 209.7] 10.32 [1] Melchers RE. Structural reliability analysis and prediction. 2nd ed. West
Sussex, England: John Wiley & Sons; 1999.
B1 [116787.9, 120464.1] 5931.3 [157.5,162.5] 8.0
[2] Kiureghian AD, Ditlevsen O. Aleatory or epistemic? Does it matter? Struct Saf
B2 [319629.9, 329691.1] 16233.0 [252.2,260.1] 12.81
2009;31:105–12.
B3 [25078.7, 25868.1] 1273.7 [75.0,77.3] 3.81
[3] Ang AH-S, Tang W. Probability concepts in engineering planning and design.
B4 [133998.7, 138216.6] 6805.4 [176.0,181.5] 8.94 Basic principles, vol. 1. John Wiley; 1975.
[4] Dempster AP. Upper and lower probabilities induced by a multi-valued
l: mean; r: standard deviation.
mapping. Ann Mater Stat 1967;38:325–39.
[5] Kendall DG. Foundations of a theory of random sets. In: Harding E, Kendall D,
editors. Stochastic geometry. New York: Wiley; 1974. p. 322–76.
[6] Ferson S, Kreinovich V, Ginzburg L, Myers DS, Sentz K, Constructing probability
1 boxes and Dempster-Shafer structures, Tech. Rep. SAND2002-4015, Sandia
National Laboratories, 2003.
[7] Williamson R. Probabilistic arithmetic. PhD thesis, University of Queensland,
Cumulative frequency

0.8 Austrialia, 1989.

[8] Penmetsa RC, Grandhi RV. Efficient estimation of structural reliability for
problems with uncertain intervals. Comput Struct 2002;80:1103–12.
0.6
[9] Tucker WT, Ferson S. Probability bounds analysis in environmental risk
assessment. In: Applied biomathematics, Setauket, 2003.
0.4 [10] Tonon F, Bae H-R, Grandhi R, Pettit C. Using random set theory to calculate
reliability bounds for a wing structure. Struct Infrastruct Eng 2006;2:191–200.
[11] Oberguggenberger M, Fellin W. Reliability bounds through random sets: non-
0.2 parametric methods and geotechnical applications. Comput Struct
2008;86(10):1093–101.
[12] Qiu Z, Yang D, Elishakoff I. Probabilistic interval reliability of structural
0
systems. Int J Solids Struct 2008;45:2850–60.
16 17 18 19 20 21 22 23
[13] Der Kiureghian A. Analysis of structural reliability under parameter
Roof drift (mm) uncertainties. Probab Eng Mech 2008;23:351–8.
[14] Ellingwood BR, Kinali K. Quantifying and communicating uncertainty in
Fig. 8. Bounds of the cumulative frequency distribution of the roof drift (frame in seismic risk assessment. Struct Saf 2009;31:179–87.
Fig. 7). [15] Reid SG. Confidence and risk. Struct Saf 2009;31:98–104.
190 H. Zhang et al. / Structural Safety 32 (2010) 183–190

[16] Tonon F. Using random set theory to propagate epistemic uncertainty through [27] Muhanna RL, Mullen RL. Uncertainty in mechanics problems-interval-based
a mechanical system. Reliab Eng Syst Saf 2004;85:169–81. approach. J Eng Mech, ASCE 2001;127(6):557–66.
[17] Dubois D, Prade H. Random sets and fuzzy interval analysis. Fuzzy Sets Syst [28] Dessombz O, Thouverez F, Laıˆné J-P, Jézéquel L. Analysis of mechanical
1991;42:87–101. systems using interval computations applied to finite elements methods. J
[18] Hall J, Rubio E, Anderson MG. Random sets of probability measures in slope Sound Vib 2001;238(5):949–68.
hydrology and stability analysis. J Appl Math Mech 2004;84(10–11):710–20. [29] Muhanna RL, Mullen RL, Zhang H. Penalty-based solution for the interval finite
[19] Rao SS, Berke L. Analysis of uncertain structural systems using interval element methods. J Struct Eng, ASCE 2005;131(10):1102–11.
analysis. AIAA J 1997;35(4):727–35. [30] Zhang H. Nondeterministic linear static finite element analysis: an interval
[20] Ganzerli S, Pantelides CP. Load and resistance convex models for optimum approach. PhD thesis, Georgia Institute of Technology, USA, 2005.
design. Struct Optim 1999;17:259–68. [31] Muhanna RL, Zhang H, Mullen RL. Combined axial and bending stiffness in
[21] Qiu Z, Elishakoff I. Antioptimization of structures with large uncertain-but- interval finite-element methods. J Struct Eng, ASCE 2007;133(12):1700–9.
non-rand parameters via interval analysis. Compt Meth Appl Mech Eng [32] Moore RE. Interval analysis. Englewood Cliffs (NJ): Prentice-Hall Inc.; 1966.
1998;152:361–72. [33] Neumaier A. Interval methods for systems of equations. Cambridge University
[22] McWilliam S. Anti-optimisation of uncertain structures using interval analysis. Press; 1990.
Comput Struct 2000;79:421–30. [34] Sun microsystems. Interval arithmetic in high performance technical
[23] Pownuk A. Efficient method of solution of large scale engineering problems computing. Sun microsystems (a white paper), September 2002.
with interval parameters. In: Muhanna RL, Mullen RL (Eds.), Proceedings of [35] Knüppel O. Profil/bias-a fast interval library. Computing 1994;53:277–88.
NSF workshop on reliable engineering computing, Savannah, GA, 2004. [36] Rump SM. Rigorous sensitivity analysis for systems of linear and nonlinear
[24] Koyluoglu U, Cakmak S, Ahmet N, Soren RK. Interval algebra to deal with equations. Math Comput 1990;54(190):721–36.
pattern loading and structural uncertainty. J Eng Mech, ASCE [37] Rohn J, Rex G. Enclosing solutions of linear equations. SIAM J Numer Anal
1995;121(11):1149–57. 1998;35(2):524–39.
[25] Möller B, Beer M. Fuzzy-randomness-uncertainty in civil engineering and [38] Ziemian RD. Advanced methods of inelastic analysis in the limit states design
computational mechanics. Springer-Verlag; 2004. of steel structures. PhD thesis, Cornell University, Ithaca, NY, 1990.
[26] Mullen RL, Muhanna RL. Bounds of structural response for all possible
loadings. J Struct Eng, ASCE 1999;125(1):98–106.