Reliability assessment based on an adaptive response surface method considering correlation among random variables

Taejun Cho1 , Changhwan Jang2 , Jeongbae Lee3 , Seongsoo Kim4

1, 3, 4Daejin University, Pocheon 487-711, Korea

2Korean Intellectual Property Office, Daejeon 302-701, Korea

4Corresponding author

Journal of Vibroengineering, Vol. 15, Issue 4, 2013, p. 1961-1971.
Received 27 September 2012; accepted 1 August 2013; published 31 December 2013

Copyright © 2013 Vibroengineering This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License

Although the Monte-Carlo Simulation (MCS) technique can evaluate a reliability of most structural systems, its processing time equals, approximately, the reciprocal of the probability of failure. While the Stochastic Finite Element (SFE) method could help to solve such a drawback, it is limited to specific computer programs, in which the mean and the coefficient of random variables are estimated by a perturbation, or by a weighted integral method. Therefore, SFE may not be easily applicable when using commercial software or systems that are not prepared with the prerequisite programming. To overcome these limitations, the RSM can be applied, because its accuracy depends on both the distance of axial points, and the linearity of the Limit State Functions (LSFs). The correlation among random variables and the response of a system is evaluated by composing a Bayesian belief nets (BBN). Consequently, the proposed Linear Adaptive Weighted Response Surface Method (LAW-RSM) with BBN modeling produces improved converged reliability indices than conventional RSMs and detail observation for the uncertainties in structural components.

Keywords: risk assessment, response surface method, adaptive weighted, Bayesian belief networks.

1. Introduction

The research presented in this paper was motivated not just by the limitations in conventional probabilistic simulation methods, but the limitation of correlation modeling for systems as well. Monte-Carlo Simulation (MCS) is usually employed to any system for evaluating the spatial and time-dependent variability of uncertainties, regardless of the strength of the correlation among its components. However, the processing time of MCS is approximately proportional to the reciprocal of the probability of failure. To overcome this disadvantage, stratified sampling in survey engineering specifically applied to MCS, Latin hyper cube method, and Markov chain modeling have been developed. Essentially, they reduce the processing time by means of algorithms specially designed to do so.

Other promising approach to reduce computational demands is the combination of the Monte Carlo based SFE with analytical procedures, such as the Response Surface Method (RSM). However, the more random variables the problem has, the more time spent to perform the FE-computations to determine points on the limit state, time that may even be higher than the one required by direct MCS.

So far, some applications of the RSM have been mentioned, but no more details of such a method have been provided. Basically, the RSM was developed by Box and Wilson (1951), and offers an alternative way for simulating systems indirectly, in less time and with less effort than traditional approaches such as MCS. Although the RSM was first used for maximizing the efficiency of a chemical factory’s operational area, now it is widely employed in many civil engineering applications (Bucher et al., 1990; Rajaschekhar et al., 1996). Such a method has the following main two advantages:

1) The implicit Limit State Function (LSF) consists of user-selected input variables, and

2) The stochastic probability can be analyzed by selecting time or space-dependent input variables for the LSF.

The LSF (also called ‘performance function’), is an expression in which each random variable must be transformed into a standard normal Gaussian random variable. To do so, a mapping process is carried out, where the physical value of a vector called ‘x’, for each random variable, should match a standard normal vector called ‘u’, by associating the probability levels of the corresponding cumulative distribution functions. This process leads to the conversion of the probability of failure into a standard-normal space.

After defining the performance function, the probability of its occurrence is easily obtained by differentiating the resultant equation. Therefore, unlike MCS, with the RSM it is possible to calculate such a probability, even for complicated structures with numerous degrees of freedom. Moreover, since the RSM reproduces the response of a complex structural system in an explicit form, it can be applied for predicting the behavior of, for instance, deteriorated concrete structures by cyclic freeze-thaw (Cho, 2007; Cho et al., 2013), or estimation the probability of bridge’s failure modes, or a technique for fatigue reliability evaluation of a steel welding member (Nowak et al., 2007; Park and Kang 2013).

Having approximated the structure’s response to a fitted explicit LSF, the RSM is next utilized to calculate the probability of its occurrence. Using the fitted function of important selected variables, the LSF is then employed for the safety evaluations, with a significantly reduced amount of structural analyses, than that of MCS (note that a complete numerical example will be detailed later, showing the RSM application).

Based on the selection of axis points, the resultant reliability index can be changed. In fact, this dependency might be improved by a response surface augmented moment method (RSMM) (Lee et al., 2006) using 3n full factorial experimental design, with 5 levels and weights, resulting from the estimation of the respective moment.

The main objective and scope of this piece of research are to review the benefits and limitations of conventional RSM techniques, and to compare their convergence while changing the linearity of LSFs, by means of a numerical example.

It is important to mention that both cases only consider the use of quadratic polynomial response surfaces, without cross product. This particular type has been used, because it can help to fill the existing gap between an exact solution (theoretically obtained), and the iteratively fitted one.

Since the latter one is a linear regression determined by an iterative method, the aim of this piece of research is to compare both quadratic and linear approaches and, precisely, identify the exact solution. In the end, it is intended to develop an improved RSM that converges more rapidly than traditional ones, regardless of the evaluation of linear or nonlinear LSFs.

Regarding correlation modeling for structural systems bayesian belief nets model is introduced in section 2.2, applied to an example in section 3.2.

2. Problems of basic response surface method

2.1. Diverging example

To analyze this type of problem, a face slab of a Concrete Face Rock-fill Dam (CFRD) has been selected (see Fig. 1). Such a section of the continuous slab is considered as a cantilever beam for the analysis, being the support the lower part. With the setting of a combination of fixed and seismic load equivalent to 0.154 g, the maximum moment external loads at the critical location of the cantilever beam were calculated with Strand-7, a commercial finite element program.

The continuous face slab of the target CFRD dam is supported by uniformly compacted foundation. Hence the settlement of structure will be ignorable. The material model is mainly composed of 24 MPa compressive strength concrete and reinforcements of 300 MPa ultimate strength, which limits surface crack in terms of over 0.4 % reinforcing volume ratio. The material, structural and geometric information for the selected elements are like the following Fig. 1. As can be seen, the maximum moment was 8.51 kN-m.

Fig. 1. a) FE Model, b) resultant factored external moment of the example structure

 a) FE Model, b) resultant factored external moment of the example structure

a) Mesh for FE Model

 a) FE Model, b) resultant factored external moment of the example structure

b)Mmax=51 kN-m

Next, three axial points were chosen to construct the response function with a gap of ± (standard deviation) between center points (the average of random variables), and axial points. The performance of the target beam was then calculated using the LSFs for estimating the probability of flexural failure, when the value of the function is less than or equal to zero.

A maximum crack width is now calculated in order to predict the probability for this width to exceed the design criteria. To do so, a response surface based on a quadratic RSM has been utilized, and compared with the results of both MCS and FOSM. Under an external service moment of 7.30∙10kN-mm at the critical location of the cantilever beam, the serviceability LSF for the crack width was evaluated as:

g = w a l l o w - w m a x ,
w a l l o w = 0.0035 × t c = 0.037  mm ,   (for a general corrosive environment) ,
w m a x = 1.08 β f s d c , min A s 100000 = 0.009  mm ,


t c = d i a m e t e r   o f   r e i n f o r c e m e n t 2 = 10.50 mm ,

d c , min = t h e   m i n i m u m   c o v e r   d e p t h = 20 mm ,

β = ( H - k D ) / ( D - k D ) = 1.068 ,

H = t h e   h e i g h t   o f   t h e   b e a m ,

k = - n p + n p 2 + 2 n p = 0.229 ,

p = 0.003775 ,

n = r a t i o   o f   e l a s t i c i t y   o f   s t e e l   a n d   c o n c r e t e = E s / E c = 200000 / 4700 f c k = 9 ,

f c k = t h e   28 t h   d a y   c o m p r e s s i v e   c o n c r e t e   s t r e n g t h ,

D = d e p t h   o f   t h e   b e a m = 400 mm ,

f s = t h e   s t e e l   s t r e s s = M e / A s D - x 3 = 14.518   MPa ,

x = k D = 86.930 mm ,

A s = t h e   c o n f i n i n g   a r e a   f o r   r e i n f o r c e m e n t s = 2 d y B N u m b e r   o f   r e i n f o r c e m e n t s =   8000 m m 2 ,

d y = d c , min = 20 mm .

Table 1 shows the random variables used in each case (where Me is the unfactored external moment), as well as their mean values and coefficients of variation (C.O.V.) assuming a normal distribution.

Table 1. Random variables and statistical values

Random variables
Notation for the Random variables
Mean value
M e (N-mm)
X 1
A s ( mm2)
X 2
D (mm)
X 3
f c k (N/mm2)
X 4
d c , m i n (mm)
X 5

The performance of the target beam has been calculated using the LSFs presented in Eqs. (4) and (5), composed by a quadratic RSM as follows:

g = w a l l o w - a 0 + i = 1 n a i x i + i = 1 n a i x i 2 ,

and the LSF for MCS or FOSM is:

g = w a l l o w - 1.08 β f s d c , min A s 100000 ,

where g() is again the LSF and ai represents the coefficients of the response function. Making use of the values calculated above, and assuming a 0.1 coefficient of variation for the wallow, the statistical values of the selected random variables are estimated in Table 5. The load term in the LSF is evaluated by RSM+FOSM, FOSM, and MCS, equivalent to the expected maximum crack width calculated as X6 in the same Table, following the Korean Concrete Institute recommendations.

The coefficients of the response surface constructed by a least square method, while varying the area of reinforcements in Table 2, are as follows:

g = w a l l o w - 0.734 + 1.899 10 - 8 x 1 - 9.564 10 - 3 x 2 - 1.151 10 - 3 x 3 + 1.012 10 - 10 x 4 + 3.869 10 - 3 x 5 + 6.262 10 - 22 x 1 2 + 4.521 10 - 5 x 2 2 + 9.769 10 - 7 x 3 2 - 2.108 10 - 12 x 4 2 - 3.869 10 - 5 x 5 2 .

The response surface functions (expressions in parenthesis in Eq. (6) are used as the demand terms in the LSFs. After constructing the LSFs of such Eq. (6), the reliability indices can be calculated by FOSM. In contrast, reliability indices for Eq. (6) were evaluated by MCS and FOSM. Three kinds of reinforcement areas are considered for varying demand terms in the LSFs. The compared results of RSM (+FOSM) against MCS and FOSM, show increased differences than those in the previous example (only 6 %), being now 53.4 % the maximum difference, when the reinforcement area is 351.9 mm2 (see Table 3).

Table 2. Coefficients of RSM with varying area of reinforcements (unit = see Table 1)

Analysis cases
Random variables
A s = 87.972
A s = 175.944
A s = 351.9
X 1 = Me
X 2 = As
X 3 = D
X 4 = fck
X 5 = dc,min
X 6 = wmax
X 6 = wmax
X 6 = wmax

Again, this value compared with the results obtained earlier, related to the differences among the simulation methods, reveals bigger discrepancies between RSM (+FOSM) and FOSM (or MCS), mainly as a result of the increased nonlinearity in the LSF – Eq. (2) and Eq. (5).

In Table 3, the three different methods are used and compared for the evaluation of the same LSFs. This example shows the diverging issues already mentioned, which lead directly to the proposition of an improved RSM.

In Table 3, the negative reliability indices indicate that there are severe discrepancies between ultimate limit state and serviceability limit state. Note that the safety margin for the serviceability limit state is relatively low, due to the risk of the two events. Having analyzed these examples with traditional approaches, now one that overcomes the limitations of the previous two will be presented.

Table 3. Reliability index and Probability of failure by three different evaluation methods for the three ultimate resistance cases of investigation with increasing area of reinforcements

Amount of reinforcement
Reliability index
Quadratic RSM
FOSM [25]
Comparison by RSM
A s = 87.972 mm2
A s = 175.94 mm2
A s = 351.9 mm2
NA*: If the probability of failure (Pf) is too small, it is almost not feasible to obtain the Pf by Monte-Carlo Simulation.

2.2. Bayesian belief network model (BBN)

The correlation among random variables for infra structures has been studied in a various way. The most well-known general modeling technique is event tree analysis (ETA) or fault tree analysis (FTA), which are convenient for qualitative analysis but hard to model for various sources. Markov process modeling overcame the weakness of ETA and FTA although it is not appropriate for modeling complex system.

The correlation modeling has been specifically dealt by bayesian update modeling technique. Since 1985, BBN could model a highly correlated causal system in an inductive way, which cannot determine the distribution of effects from results to cause direction. Therefore, the inverse way of BBN is proposed in this research in order to identify the role of random variables to consequential result considering correlation among random variables of each member.

BBN are a graphical structure known as a directed acyclic graph (DAG). A DAG is a set of nodes and arcs, or directed edges, between nodes, such that there is no directed cycle.

BBN is typically built by first finding an appropriate structure either by interviewing an expert, or by selecting a good model from training data, and then using a training sample to estimate the parameters (Darwiche, 2000). The resulting belief net is then used to answer queries based on joint probability.

Compared with event tree analysis (ETA) or fault tree analysis(FTA), BBN could evaluate and detect a partial failure in a probabilistic quantification, which is not allowed in ETA or FTA. It additionally overcome the exponentially increased complexity in stochastic modeling of Markov Model by employing zipping techniques. This feature is shown by showing superiority to Markov Chain Monte Carlo Simulation as well, i.e., Metropolis method or Gibbs Sampling in terms of simplification by zipping.

The steps of uncertainty analysis in a BBN model are as follows:

1) defines a set of random variables;

2) defines a joint distribution for these random variables;

3) defines other variables which are functions of the random variables;

4) Monte Carlo samples the entire joint distributions (probabilistic and functional variables);

5) communicates the probabilistic prediction results,

where, the random variables are assigned marginal distributions, followed by a specification of a DAG to capture conditional relations. Probabilistic influence between parent and child is represented as conditional rank correlation. A joint probability density for the probabilistic nodes is built using the joint normal copula to realize the dependence relations. This density can be updated analytically, or sampled and analyzed with Monte Carlo Simulations.

3. The linear adaptive weighted response surface method

3.1. Improvement techniques

In spite of the two issues discussed before, the merits of the basic RSM are remarkable, and perhaps that is the reason why it has been widely applied for carrying out reliability analysis. The approximation of structural responses, however, would show relatively large errors depending on the nonlinearity form of the LSFs. Irfan et al. (2005) proposed a weighted regression method in which the response surface function was formulated by assigning higher weights to the variable that was closer to the limit state. To do so, they adopted an n×n diagonal matrix of weights and the best design, amongst the responses from the performance function corresponding to the design matrix, was selected. Consequently, the distance among input random variables is optimally determined. This occurs when the LSFs approaches to zero. Thus, the best response can be given as:

f b e s t = m i n g x i ,

where fbest is the determined best value amongst the axis points, and g(x)i is the fitted response on the response surface, composed by center and axis points at the iteration step.

The component of the weighted matrix in each iterative calculation can be expressed in the following manner:

w i = exp - g x i - f b e s t f b e s t ,

where wi is the component of the weighted matrix, and needs to be multiplied by the center or axis points, for composing an improved response surface for the next iterative step. Consequently, the longer arm length between axis points is reduced, providing closer design points than the old axis points.

To improve the convergence and accuracy when dealing with nonlinear LSFs, a first approximation could be to integrate the weighted matrix with the adaptive iteration developed by Rackwitz-Fiessler, leading to the establishment of the adaptive weighted response surface (AWRS) function. However, having combined the improved techniques (both linear and nonlinear response surface functions using basic, adaptive, and adaptive weighted method) the authors concluded that the linear adaptive weighted response surface method (LAW-RSM), was the best alternative. The following examples are aimed at demonstrating this statement.

3.2. Validation of the LAW-RSM

The random variables in the example of the simply supported beam selected for analysis are: distributed load (W), span length (L), modulus of elasticity (E), and the moment of inertia (see Fig. 2).

Fig. 2. Scheme of the beam considered in Example 1

 Scheme of the beam considered in Example 1

The limit state being considered is deflection, and the allowable deflection is specified to be L/360. The maximum deflection is 0.0069WL4/EI, and it occurs at 0.446L from either end of the beam.

The limit state function is:

g W , L , E , I = L 360 - 0.0069 W L 4 E I .

The mean and coefficient of variation for each of the variables under study are presented in Table 4.

Table 4. Parameters of random variables in Example 1

X 1
8E-4 m4
Moment of Inertia
X 2
2.00E+07 kN/m2
Elasticity Modulus
X 3
10 kN/m
Distributed Load
5 m
Span length

In the example considered, the linear RSM was accurate with the weighted regression method. As can be seen in Table 5, while the application of the quadratic LAW-RSM resulted in a difference of 42.23 %, compared with the results of the classical MCS, the LAW-RSM showed less than 3 % of difference with such a simulation approach. However, the convergence and accuracy will be changed depending on the linearity order of limit state function.

Table 5. Comparison of converged reliability indices from applied RS methods

Analysis method
Reliability Index, β
Probability of Failure, Pf
Error (%)
Number of Simulations = 1,000,000
Linear RSM
1. Basic
2. Adaptive
3. Adaptive Weighted
Rackwitz-Fiessler with Exponential weighting
Quadratic RSM
4. Basic
5. Adaptive
6. Adaptive Weighted
Rackwitz-Fiessler with Exponential weighting

3.3. Correlation model for random variables in BBN

BBN were a convenient alternative to develop the model. They are graphs, whose nodes are connected by arcs. While nodes represent variables, arcs represent the strength of the relationship between them (Hanea, 2008). The employed bayesian belief net is answer queries for the following probability density function based on conditional probability:

P D e f l e c t i o n I n e r t i a , E l a s t i c i t y , L o a d I n t e n s i t y , S p a n l e n g t h ,

Following Markov condition for BBN (chain rule) the above joint probability is then calculated as:

P X 1 , X 2 , , , X n = i = 1 n P X i p a r e n t s X i ,

where, Pxy is a coefficient of correlation between random variables, xi and uidenotes each random variable, n represents the total number of random variables, Sx,Sy are standard deviation of variables, and x-,y- are mean value of random variables.

Four variables were considered as important for the model, they are: (i) moment of inertia, (ii) elasticity, (iii) distributed load intensity and (iv) length of a span (see Table 4). The analytical conditioning in BBNs on four probabilistic nodes can be conditionalized, whereby the conditional distribution is computed in a few minutes. This is possible with the joint normal copula. Random variables are treated as transformations of joint normal variables.

If Z=(Zi,...,Zn) denotes a joint normal distribution of normal variables Zi with mean zero and unit variance, then random variables Xi,...,Xn are written as

X i , . . . , X = F 1 - 1 Φ Z 1 , F n - 1 Φ Z n ,

where Fi is the cumulative distribution function of random variable Xi, and Fn is the cumulative distribution function of the standard normal variable. If Fi is continuous invertible, this transformation is rank preserving, the rank correlations of Xi, Xj will be the same as that of Zi, Zj.

The rank correlation of normal variables is computed from the product moment correlation with the Pearson transformation. For details see (Hanea et al., 2006)

Conditioning on the value Xi=x entails conditioning Z on Zi=Φ-1(Fi(x)). Letting Z(Xi=x) denote this conditional distribution, the conditional distribution of Xi,...,Xn, when Xi,...,Xn are joined by the joint normal copula, is given by:

X 1 , X n X i = x = F 1 - 1 Φ Z 1 X i = x , F n - 1 Φ Z n X i = x .

The results of conditioning among random variables could be presented in graphical format in various ways. As showin in Fig. 3, if used with the results of Monte Carlo simulation, it may be hard to understand due to the complex correlation.

Fig. 3. Correlation among random variables shown in Cobwebs

 Correlation among random variables shown in Cobwebs

Fig. 4. a) BBN model among random variables, b) probability density function of each varibale

 a) BBN model among random variables, b) probability density function of each varibale


 a) BBN model among random variables, b) probability density function of each varibale


Therefore we have selected quantified regression analysis between two variables. One of the most important variables in structural analysis would be deflection. The correlation between deflection and other variables are investigated presented in Fig. 4 and Fig. 5. By the results, it turns out that due to the contribution of span length to the load intensity the correlation between deflection and span length show the largest sensitivity while modulus of elasticity show the lowest correlation, shown in Figure 5. This is quite different with simply supported beam case, where deflection is directly proportional to Moment load while it is reversely proportional to elasticity of modulus times inertia of secondary moment of the sections.

Consequently, while the proposed improved RSM method could present faster converged and improved evaluation of reliability of structures, Bayesian belief net model additionally provide contribution of each random variables in terms of correlation among components.

Fig. 5. Correlation and conditional expectation of deflection on other random variables

 Correlation and conditional expectation of deflection on other random variables

a) Correlation of defl ection on random variables

 Correlation and conditional expectation of deflection on other random variables

b) Conditional expectation of deflection on random variables

4. Conclusions

A few stage analyses have shown that the RSM could construct a regression model, from the responses of any complicated structural system, which may be very efficient for carrying out a reliability analysis, or a statistical prediction from the constructed response surface. This in turn helps to save time and effort significantly when performing millions of iterative calculations. However, the accuracy of RSM varies depending on the distances of axial points, as well as the linearity of the LSFs. Selecting more efficient techniques from previous improved RSM models, and having discarded inexact or slowly converged ones, an improved linear adaptive weighted response surface method (LAW-RSM) was developed.

The weighting matrix was determined based on the closeness to LSF as a zero value, which additionally increases the convergence with closeness to the failure point and optimal determination of the distance among input random variables. To find the failure point more effectively, new centre points were iteratively calculated by the Rackwitz-Fiessler method, rather than the Bucher’s approach.

Based on the findings, a computer program has been developed and verified, by means of two examples, with the LSFs in the linear and quadratic forms of response surface functions. The developed LAW-RSM shows better converged results, regardless of the linearity of LSFs, than the other five RSM methods considered in this piece of research. By updating design points to the failure points, with weighting on important random variables in LAW-RSM, the authors believe that the two main features of engineering analysis, the accuracy and the stability of calculations, have been obtained by constructing the LSF.

In addition, conditioning on random variables is modeled in a Bayesian belief net model (BBN). BBN provides the quantification of correlation among structural components, which is essential in more economic and safer structural design.


This research was a part of the project titled “Development of the Advanced Technology of Nuclear Power Plant Structures Quality on Performance Improvement and Density Reinforcement (201016101004L)” in the Nuclear Power Technology Development Project funded by Korea Institute of Energy Technology Evaluation and Planning (KETEP) and the Ministry of Knowledge Economy, Korea.


  1. Box G. E. P., Wilson K. B. On the experimental attainment of optimum conditions. Journal of Royal Statistical Society, Vol. 13, 1951, p. 1-45, 507-516. [Search CrossRef]
  2. Bucher C. G., Bourgund U. A fast and efficient response surface approach for structural reliability problems. Structural Safety, Vol. 7, Issue 1, 1990, p. 57-66, 507-516. [Search CrossRef]
  3. Cho T. Prediction of Cyclic Freeze-Thaw Damage in Concrete Structures Based upon Response Surface Method. Construction and Building Materials, Vol. 21, 2007, p. 2031-2040. [Search CrossRef]
  4. Cho T., Kim S., Lee J., Jang C. An improved Bayesian inference model for auto healing of concrete specimen due to a cyclic freeze-thaw experiment. Journal of Vibroengineering, Vol. 15, Issue 2, 2013, p. 884-891. [Search CrossRef]
  5. Hanea A. M., Kurowicka D., Cooke R. M. Hybrid Method for Quantifying and Analyzing Bayesian Belief Nets. Quality and Reliability Engineering International, Vol. 22, 2006, p. 709-729. [Search CrossRef]
  6. Irfan K., McMahonb C. A. A response surface method based on weighted regression for structural reliability analysis. Probabilistic Engineering Mechanics, Vol. 20, 2005, p. 11-17. [Search CrossRef]
  7. Lee S. H., Kwak B. M. Response surface augmented moment method for efficient reliability analysis. Structural Safety, Vol. 28, 2006, p. 261-272. [Search CrossRef]
  8. Nowak A. S., Cho T. Prediction of the Combination of Failure Modes for an Arch Bridge System. Journal of Constructional Steel Research, Vol. 63, Issue 12, 2007, p. 1561-1569. [Search CrossRef]
  9. Rajaschekhar M. R., Ellingwood B. R. A new look at the response surface approach for reliability analysis. Structural Safety, Vol. 12, 1993, p. 205-216. [Search CrossRef]
  10. Park Y., Kang D. Fatigue reliability evaluation technique using probabilistic stress-life method for stress range frequency distribution of a steel welding member. Journal of Vibroengineering, Vol. 15, Issue 1, 2013, p. 77-88. [Search CrossRef]