1. science
2. physics

Journal of Information & Computational Science 12:2 (2015) 657–672
Available at http://www.joics.com
January 20, 2015
Stochastic Sensitivity Analysis for Models with Parameter
Correlation and Model Approximation Method ?
Xian Dong a , Zhan Wang a,b , Jianrong Pan a,b,∗
a Department
b State
of Civil Engineering, South China University of Technology, Guangzhou 510640, China
Key Laboratory of Subtropical Building Science, South China University of Technology
Guangzhou 510640, China
Abstract
In the stochastic sensitivity analysis, a large number of simulation models lead to low computational
efficiency. With the dimension increasing in assigned problems, the accurate is difficult to achieve
by popular regression methodologies. Without considering the correlations of parameters, inaccurate
sensitivity coefficients would be calculated by analyzing the effect of parametric variables on structures,
moreover, the existing methods only calculate the local gradient as sensitivity. According to these
problems, approximate model and global sensitivity method are employed for design sensitivity analysis
of structures. The approximate model is constructed by the hybrid neural network which possesses
significant learning capacity and generalization capability with a small amount of information. The
improved chaotic particle swarm is applied to optimize the neural network so as to improve the
computational efficiency and accuracy. The proposed sensitivity analysis method values the global
response of the outputs by varying all the input parameters at a time with the correlations of parameters.
Uniform design and Latin hypercube sampling are used to sample points. Numerical analysis show that
the proposed method can successfully measure the actual sensitivity.
Keywords: Approximate Model; Sensitivity Coefficient; Correlation between Parameters; Improved
Chaotic Particle Swarm; Stochastic
1
Introduction
System optimization or design reinforcement require estimating the influence of various parameters on structures, that is sensitivity analysis. The actual effect of engineering structures caused
by mechanical properties of material, stochastic external loads and stochastic section errors would
be inadvertent through deterministic parameters models, therefore the uncertainty and stochastic
of parameters [1, 2] need to be considered in complicated structures. In recent years, the methods
?
Project supported by the National Nature Science Foundation of China (No. 51178192, 51378009, 51378219),
and Key Laboratory of Subtropical Building Science Foundation, South China University of Technology
(No. 2012ZA05).
∗
Corresponding author.
Email address: [email protected] (Jianrong Pan).
1548–7741 / Copyright © 2015 Binary Information Press
DOI: 10.12733/jics20105305
658
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
of uncertainty and sensitivity analysis for structures have been developed by domestic and foreign scholars including Monte Carlo simulation [3], stochastic finite element method [4], response
surface method [5, 6] and so on.
The major advantage of Monte Carlo simulation is simplicity in formulation, however there
is a consequence of this method that too many complex models cause high computational cost.
The polynomial functions for inputs and outputs would be established based on response surface
method which is convenient for compiling calculation programs, however it is nearly impossible
to achieve the satisfactory precision when the parameters of function having highly nonlinear
relation. Stochastic finite element method which including perturbation stochastic finite element
method and spectral stochastic finite element method has complex random theories. Variation
coefficients would be limited in a special range conducting analysis by the method of perturbation
stochastic finite element.
Analysis of structure stochastic sensitivity needs for a large amount of simulations, computational cost is higher for the structure is more complex, so far, the most effective way is to replace
the physical models by establishing the approximate models [7, 8]. The sensitivity techniques have
been classified into two groups: local sensitivity analysis methods and global sensitivity analysis
methods. The local sensitivity analysis methods reflect the local response of the output(s) by
varying input parameters one at a time then holding other parameters at determined values, the
global sensitivity analysis methods reflect the global response of the output(s) in a finite region
[9]. The local sensitivity analysis technique has been widely used for easily implement, but the
effect of global parameters would be covered for physical parameters with different dimension
magnitudes. The correlation of parameters as an essential component effects the output results
in complex structure systems, in addition, the inaccurate results would be obtained by neglecting
the parameter relativity. By now the sensitivity analysis methods have evolved as a major area
of research, but only a few studies have considered the correlations between input parameters in
the sensitivity analysis of structures.
In this paper, the analysis of stochastic parameters sensitivity can be divided into two procedures: approximate model representation and stochastic sensitivity analysis. In Section 2, we
use hybrid neural network (optimized by the method of improved Chaotic particle swarm) to
approximate the structure models so as to improve the computational efficiency, avoid computing
large finite element models. As an application of improved chaotic particle swarm to optimize the
neural network, the high magnitudes nonlinear function of complex structure would be built by
limited samples. Section 3 transforms the correlated variables to independent normal variables
equivalently by nataf model, Section 4 develops a stochastic sensitivity method to estimate the
global effect of numerous parameters based on varying input parameters at the same time. Section
5 tests some representative examples and shows that the proposed procedure is useful in solving
the problems that neural network is easily to fall into local minimum and slow convergence etc.
2
2.1
Method of Establishing Approximate Model
Algorithm of Feed-forward Neural Network with Back-propagation
Learning
Neural network has the behaviors of self-learning, self-organization and simulate nonlinear
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
659
relationship, it will simulate the complicated nonlinearly relation between inputs and outputs
based on the sample space divided by continuous hypersurface which is constituted of bounded
continuous differentiable functions. The fundamental of feed-forward neural network is it adjusts
the weights and threshold value by error back propagation, which takes use of errors caused by
the forward propagation of input layer, hidden layer and output layer feeding back to neurons of
each layer.
(1) The mean square error is the sum of each neuron error when the output has been derived
from the network, one sample error E is given by:
n1 −1
1X
E=
(dl − yl )2
2 l=0
(1)
The sum error of all the samples is given by:
n
Ec =
n −1
1
1XX
(dl − yl )2
2 n=1 l=0
(2)
The network weight is adjusted by gradient descent method
∂Ec
4wik = −
∂wik
(3)
where dl is the sample output, yl is the network output.
(2) When the output layers have iterated t times the equations are given by:
∂Ec
0
∂wkl
0
p
p
p
X
X
X
0
0
∂Ec
∂E p
∂E p ∂ylp ∂ulp
0
=
=
=−
(dpl − ylp )f (ulp )xkp
0
0
0
0p
p
∂yl ∂ul ∂wkl
∂wkl
∂wkl
p=1
p=1
p=1
0
0
wkl (t + 1) = wkl (t) − η
=−
p
X
0
(dpl − ylp )ylp (1 − ylp )xkp
p=1
=
− ylp )ylp (1
t
X
0
0
0
∂ulp =
wkl xkp
k=0
0
p
yl = f (ulp )
p
δkl
(dpl
(4)
(5)
− ylp )
(6)
(7)
(8)
0
only affects yl .
where p is the number of samples, E p is a function about y0 , y1 , . . . , ym−1 , wkl
(3) Hidden layers neurons are trained by incremental learning rules such as:
p
X
∂Ec
p p
δik
xi
=−
∂wik
p=1
wik (t + 1) = wik (t) − η
p
=
δik
m−1
X
l=0
0
0
(9)
∂Ec
∂wik
0
p
wik xip (1 − xip )
δik
(10)
(11)
660
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
wik (t + 1) = wik (t) + η
p
X
p p
δik
xi
(12)
p=0
The input layers, hidden layers and output layers are connected by the neural network with the
map of n-dimension space to m-dimension space. The number of hidden neurons is n1 , wik is the
weight of input layers to hidden layers, bk is the threshold value of input layers to hidden layers,
0
0
wkl is the weight of hidden layers to output layers, bl is the threshold value of hidden layers to
output layers, f is a transfer function. The calculation process is expressed as:

n−1
P

 x0k = f (
wik xi − bk )
k = 0, 1, 2, . . . , n1 − 1

i=0
(13)
n−1
P 0 0
0


l = 0, 1, 2, . . . , m − 1
wkl xk − bl )
 xl = f (
k=0
The learning rate, momentum coefficient and other parameters have a serious impact on the
feedforward network. The network will fall into local extremum and could not obtain the optimal
weight distribution for the improper values of these parameters, therefore, in this paper the
network will be optimized by chaos particle swarm algorithm.
2.2
The Algorithm of Traditional Chaotic Particle Optimization
Particle swarm optimization originated from the behavior of biotic population and has widely used
to solve optimization problems. In this algorithm, any candidate solution of optimization can be
analogy to a particle in v-dimensional search space. These particles moved in a v-dimensional
search space with a certain speed which adjusted by the positions and velocities of its own
and companions’. The qualities of all the particles are evaluated by a fitness function which is
determined by the optimization goal. Each particle is influenced by the fitness function moving
toward to its best local position, therefore the swarm moves toward to the best global position in
an identified space.
A biotic population modal is represented as X = (x1 , x2 , ..., xk )T consisting of k particles in vdimensional search space, the ith particle’s position is xi = (xi,1 , xi,2 , ..., xi,v )T and the velocity is
vi = (vi,1 , xi,2 , ..., xi,v )T . The equations describing the process of the particle swarm optimization
are given by:
j+1
j
vi,v
= wvi,v
+ m1 rand()(pbestji,v − xji,v )m2 rand()(gbestji,v − xji,v )
(14)
j
j+1
xj+1
i,v = xi,v + vi,v
(15)
where rand() is a random number in interval (0, 1), m1 , m2 are acceleration constants. Where
v
denotes the ith0 s velocity in v-dimensional space which iterated in the jth time, xji,v is the
vi,j
ith0 s position in the vth dimensional space which is iterated in the jth time; pbestji,v denotes the
best position of the particle i in the vth dimensional space, gbestji,v is the global best position of
the swarm in the vth dimensional space. w is a decline inertial weight decreasing from 0.9 to 0.5
[10] with the increasing iteration times and w can be expressed as:



 m1 rand() ∈ [0, 2]
m2 rand() ∈ [0, 2]


k
 w(k) = −0.5
+ 0.9
maxnumber
(16)
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
661
Chaos is a motion state, that it is seemingly irregular but determined and unpredictable ubiquitous in certainty nonlinear dynamical systems. The characteristic of chaotic system is extremely
sensitive to infinitely-small initial changes and perturbation. It will deviate from the original
direction of evolution after prolonged exercise even by a condition of tiny disturbance. The
unpredictability of the chaotic system derives from unstable state belonging to the system.
One of typical chaos system is logistic mapping model. Suppose that xn is a parameter related
to the individual number of an insect in the nth year and xn+1 is the number of the (n + 1)th
year, both xn and xn+1 are all integers. The equation about the parameters xn and xn+1 can be
described as v = f (xn ), n = 1, 2, . . ., where n is a nonnegative integer. Logistic model can be
given as:
xn+1 = µxn (1 − xn ), n = 1, 2, . . . ,
(17)
where µ is a control parameter deciding the state of the system, when the value rang of µ is [3, 4],
the dynamic action of the system including two states period-doubling and chaos. If µ takes the
value range [3.56994567, 4], the system will get into chaotic state, if µ = 4, the system will be
totally in chaotic state, the algorithm which making the swarm chaos with µ = 4 can be carried
out as following steps:
(1) Map xi,v to the interval [0, 1.0], xi,v is the location of particle xi in the vth dimension, [av , bv ]
is the domain of xi in the vth dimension.
hxi,v =
xi,v − av
b v − av
(18)
= 4hxji (1−hxji )
(2) Generate a chaotic sequence h1x,v , h2x,v , . . . , hJx,v by iterating an equation hxj+1
i
for j times;
(3) Map particles of the chaotic sequence to the primary searching space following the function
xji,v = av + hji,v (bv − av ), chaos column points xji = (hji,1 , hji,2 , . . . , hji,n ) in each dimensional
chaotic sequence can be derived from Tent mapping.
2.3
Improved Chaotic Particle Swarm Optimization
The phenomenon of particle swarm premature convergence is increasingly outstanding when solving complex optimization problems. So far the main method to solve the problem is increasing
the size of the particle swarm, but this approach could not solve this problem fundamentally.
Chaotic variables which are stochastic motion state variables obtained by deterministic equations. Using chaotic variables to optimize the searching process can jump out of local minima,
keep the diversity of population, improve global optimization performance [11].
There will be no significance in the process of searching the optimal solution if the new particles
which generated by chaotification always in the state of chaos. The algorithm can be constantly
close to the optimal solution only under the alternate states of chaos and stability [12].
Fitness functions according to the objective function are established to estimate whether the
particles are in premature convergence or not. In this paper, we set the relative errors of the
network outputs and the actual response as the fitness functions based on the algorithm of back
662
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
propagation. To insure the functions are nonnegative, the expressions are represented as follows:
v
u n
n
.X
uX
2
t
u(x) =
(yi − f (xi ))
f (xi )
i=1
yi = f
n
³X
0
0
0
´
wki xk − bi
k=0
0
xk = f
m
³X
(19)
i=1
(20)
´
wlk xl − bk
(21)
l=0
where u(x) is the judgment function of premature, yi is the network output, xi is the network
input, x0k is hidden layer neuron.
The status of particle swarms that it is in steady state or chaos is determined by the movement
distance of current particle, the equation about it can be given as:
smove = abs(xi (t) − xi (t − 1))
(22)
where t is a number of the particle evolutionary generations. It is considered that algorithm has
reached to the local optimization when smove or u(x) is less than a limit value, and the variables
will get into chaotization by chaotic searching. The concrete steps are represented in Section 2.2.
2.4
Optimization of Neural Network by Chaotic Particle Swarm Algorithm
Encode matrixes of every particle by coding strategy, the matrixes are consisted of two kinds
particles, the first kind particles are the weights of input layers to the hidden layers and hidden
layers to the output layers, the second kind particles are the threshold values of input layers to
hidden layers and hidden layers to the output layers. These matrixes are expressed as:
P article(1) = [W1 , W2 ]


w1,1
w1,2
. . . w1,n

 w2,1 w2,2 . . .

W1 =  .
..
..
 ..
.
.

wm,1 wm,2 . . .

w
wn+1,2
 n+1,1
 wn+2,1 wn+2,2

W2 = 
..
..

.
.

wn+k,1 wn+k,2
w2,n
..
.






wm,n
. . . wn+1,m
. . . wn+2,m
..
..
.
.
. . . wn+k,m
(23)







(24)
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
P article(2) = [Z1 , Z2 ]

663

z1,1

 z2,1

Z1 =  .
 ..

zm,1

z
 1,k+1
 z2,k+1

Z2 = 
..

.

zm,k+1
z1,2
. . . z1,k
z2,2
..
.
. . . z2,k
..
..
.
.






zm,2 . . . zm,k
z1,k+2
. . . z1,k+g
z2,k+2
..
.
. . . z2,k+g
..
..
.
.
(25)







(26)
zm,k+2 . . . zm,k+g
where n is the dimension of input layer variables, m is the number of hidden layer neurons, k is the
dimension of output layers, g is the population size. The computational process of optimization
hybrid neural network are described as follows:
(1) Make parameters initialized. Assuming that the size of particle population is G, inertia weight
is w, learning factors are m1 and m2 , the maximum value of particle swarm evolutionary
generations is T , the energy function is E, the iteration number of chaotic map is Q, the time
of function initial changed is f = 0, the maximum number of energy changing times of PSO
is fmax , the judgment of premature convergence is u, the judgment of steady state is s. Set
the value of current evolutionary generation is t = 1, the weights and the threshold values of
the network are separately represented by two kinds of particles which randomly generated
in feasible region.
(2) The energy functions are constructed by the fitness function and the movement distance
function of particles. The global best position of initial particles swarm which is expressed
as gbest would be selected by the initial fitness value. The best position of current particles
swarm is expressed as pbest.
(3) The position and velocity of the new particle are updated according to the equations14 and
15. After estimating the energy values of the new generated particles, go straight to step 4
for the ith particle value is higher than value u and s, if the energy value of the ith particle
is lower than value u or s, selects the position of the particles’ lowest energy value to be the
new position that the particles are q chaotic column points which have generated by logistic
map or tent map.
(4) Evaluate the energy values of the new particles are estimated, if the energy values are lower
than the values which are in before the update, the value of pbest and the corresponding
position will be updated, if the lowest energy value is lower than gbest, the value of gbest and
the corresponding position will be updated and then f = 0, otherwise, f = f + 1.
(5) If f ≥ fmax , the error back propagation is used to have a local search near the gbest. The
searching result will instead of pbest and update the corresponding position for it is better
than gbest, otherwise, the searching result will instead of the worst particle of pbest, f = 0,
t = t + 1. The process is finished when t ≥ T , otherwise, go back to step 3.
(6) Save the results and optimize the network.
664
3
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
Transformation between Correlated Variables and Independent Variables
In order to achieve the equivalent transformation between correlated variables and independent
variables, there are three problems to be solved, the first is Joint probability density functions of
correlated random variables could not be determined; the second is correlation coefficients would
have changed with the normalization of random variables; the third is changing the correlation
coefficients to zero is not equivalent to transform correlated variables to independent variables,
in order to conquer these problems, Section 3.1 and 3.2 would discuss in detail.
3.1
Transformation between Correlated Non Normal Variables and
Correlated Normal Variables
That normal variables are independent when the correlation coefficients are zero is a mathematical
theorem, so firstly, converts non-normal random variables to normal random variables by the
functions of equal probability transformation. Consider a random vector X = [X1 X2 . . . X3] and
a standard normal distribution vector Y = [Y1 Y2 . . . Y3 ]. Then the relationship between random
distribution and normal distribution are given by:
(
Fi (xi ) = Φ(yi )
i = 1, 2, . . . , n
(27)
xi = Fi−1 [Φ(yi )]
Joint probability density functions deriving from the rules of implicit function derivation are
expressed as follows:
f1 (x)f2 (x) . . . fn (x)
fx (x) =
φn (y, ρ0 )
(28)
φ(y1 )φ(y2 ) . . . φ(yn )
1 T −1
1
e( 2 y ρ0 y)
φn (y, ρ0 ) = p
(29)
n
(2π) det(ρ0 )
where Fi (·) is a marginal probability distribution function, Fi−1 [·] is the inverse function of Fi (·),
Φ(·) is a standard normal distribution function, ρ0 is the equivalent correlation coefficient matrix
of Y , fi (·), i = 1, 2, . . . , n is the marginal probability density function corresponding with Fi (·),
φ(·) is the probability density function of standard normal variables, det(·) represents the matrix
determinant, (·)−1 is the function of matrix inversion. The correlation coefficient ρ0xy of standard
normal random vector y can be derived from the following equation:
Z +∞ Z +∞ ³
xi − µxi ´³ xj − µxj ´
fXi ,Xj (xi , xj ) dxi dxj
ρij =
σxi
σxj
−∞
−∞
(30)
Z +∞ Z +∞ ³ −1
−1
Fi (Φ(yi )) − µxi ´³ Fj (Φ(yj )) − µxj ´
=
φ2 (yi , yj , ρ0ij ) dxi dxj
σxi
σxj
−∞
−∞
where ρij is the correlation coefficient between xi and xj of random vector X, ρ0ij is the correlation
coefficient matrix component of standard normal random vector y, fXi Xj (xi , xj ) is the joint probability density function between the ith and jth variables. The matrix ρ0 could be determined
by the nonlinear function (30) for the marginal probability density function of random variables
and ρ have been calculated.
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
3.2
665
Transformation between Correlated Normal Variables and Independent Normal Variables
Suppose that the correlation coefficient matrix ρ0 have been obtained by Section 3.1, for an
intermediate variable M is represented as:
√
M =V · Λ
(31)
where Λ is a diagonal matrix composed by the eigenvalues λ of the correlation coefficient matrix
ρ0 , V is a matrix composed by the eigenvectors v corresponding with the eigenvalues λ. The
relationship between independent standard normal variables U and correlated standard normal
variables Y can be represented as follows:
U = M −1 · Y
4
(32)
Stochastic Sensitivity Estimation for Structure Systems
Using Approximate Model and Monte Carlo
To analyze stochastic response of structure by Monte Carlo method needs large number of simulations. In this paper, hybrid neural network is used to establish a reasonable network between
random variables and structure response with selecting small-scale certain samples, in addition,
the hybrid neural network is optimized by the algorithm of improved chaos particle swarm, moreover enormous workload would be avoided by the approximate model.
4.1
4.1.1
The Sampling Strategy of the Approximate Model
Uniform Design Sampling Method for Neural Network
Uniform design [13] is a method that setting the experiment points uniformly, it’s algorithm is
the combination between number theory and multivariate statistical. More information about
the experiment model can be obtained from the fewer experiments by uniform design, so it is
suitable for the multifactorial experiments and the model system parameters which are completely
unknown. The regression model is given by:
J = E(x1 , x2 , . . . , xn ) + ε
(33)
where E is a known function class, the value intervals of xn are C n ⊂ [0, 1]n .
The uniformity of the design points is measured by deviation. Assume that γ = {xi | i =
1, 2, . . ., n} are the n test points of C n , g is an arbitrary function of E satisfy the conditions of
g ∈ C n , N (g, ) is the number of points which are in the range g ∈ C n , the deviation of point set
y in C n is
¯
¯
¯ N (g, γ)
¯
e(n, γ) = sup ¯¯
− (ν[0, g])¯¯
(34)
n
(g∈C n )
where v([o, g]) denotes the rectangular volume which consists of n points, sup denotes the upper
bound. In this section, we adopt good lattice point method to construct an algorithm which
666
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
making the points distributing uniformly in the range of experiment. (1) Make an assumption
that n is the experiment times, c is a number that it is less than n, C is a vector consist of
positive integers c that the greatest common divisor between n and c is 1. (2) The jth column
of uniform design table is consist of uij = icj [modn], if icj is more than n, subtracting h which
is an appropriate time of n from icj and the difference between icj and h is in interval [1, n], the
uij are given by:
 j

u1 = c j


(


uij + cj
uij + cj ≤ n
ui+1,j =

uij + cj − n
uij + cj ≥ n




i = 1, . . . , n − 1
4.1.2
(35)
Process of Random Sampling
Latin Hypercube Sampling (LHS) is a method for generating a sequence of points from multidimensional distributions [14]. The distinction between Latin hypercube sampling and direct
sampling is that taking advantage of LHS to avoid sampling repeatedly. For an assumption that
the range of each variables is divided into n equal probability intervals, F (xi ) is a probability
distribution function of random variables. The n + 1 boundary points which have been generated
in domain are given by:
³k´
(k)
k = 1, 2, . . . , n + 1
(36)
bi = Fi−1
n
The samples in each equal probability interval are expressed as:
(k)
xi
=
Fi−1
³ k − 0.5 ´
n
k = 1, 2, . . . , n
(37)
The sampling process of considering the correlation between parameters are listed as below
(1) Generate n-dimensional random variables which follow normal distribution, assume a computational accuracy θ, n-dimensional random variables X;
(2) Calculate the correlation coefficient matrix of normal distribution random variables;
(3) Transform correlated random normal variables to random independent variables by the algorithm of Section 3;
(4) Calculate the correlation coefficient matrix ρ of the original model, obtain the equivalent
correlation coefficient matrix ρ0 according to the literature [15];
(5) Generate random variables following normal distribution with the correlation coefficient ρ0
according to the inverse transformation of Section 3, calculate the absolute values of the
differences between correlation coefficient matrix of the generated variables and ρ0 , if all the
differences are less than θ, go to step 6, otherwise, go back to step 2;
(6) Transform the correlated random normal variables correlation coefficient matrix ρ0 to the
correlated random non-normal variables based on the equal probability transformation.
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
4.2
667
Calculation of Structure Sensitivity Coefficients
The sensitivity of each parameter is determined by the output variations corresponding to the
range of each input, the outputs are the predicted results based on neural network and the inputs
are the random samples according to probability density function. Considering the impact of each
random variable on structure response, correlation coefficients are used to measure the correlation
between random variables and structure response, the coefficients of variation are brought in to
eliminate the effect of the random variable dimensions and measurement scales. The equations
are given by:
n
P
(xi − x)(f (xi ) − f (xi ))
i=1
Rr = r n
(38)
P
2
2
(xi − x) (f (xi ) − f (xi ))
i=1
v
u³ u
X Rr ´
Rr u
Sr =
/t
γr
γ
r=1 r
(39)
where Rr is the Pearson correlation coefficients, xi is random variable, n is Monte Carlo sampling
times, f (xi ) is the network output according to xi , x is the mean value of xi , f (xi ) is the mean
value of f (xi ); γr is the coefficient of variation according to Monte Carlo random samples, u is a
number represented the random variables total dimensions, Sr is the sensitivity value of random
variables, the absolute value of coefficient reflects the influence of random variables on structure
response, positive and negative reflect the degree of correlation between random variables and
structure response. The flow chart of structure stochastic sensitivity analysis based on Section 2,
Section 3 and Section 4 is shown in Fig. 1.
5
Examples and Analysis
In this section, comparing approximate models which are established by the method of Section
2 with Monte Carlo to verify the accuracy and then two engineering examples are selected to
calculate sensitivity values by stochastic analysis. To assess the accuracy of the meta-model,
the values of 10000 additional test samples which obtained by Monte Carlo are employed approximately as exact solutions. It will demonstrate the superiority of this algorithm since the
performance functions of engineering structures are complex implicit functions. The Monte Carlo
results considering relativity are abbreviated to MC, the results of hybrid neural network are
abbreviated to HNN.
5.1
5.1.1
Verify the Approximate Model
Example 1
Considering a non-linear function with three uniform random variables, f (x) = x21 + x22 + x23 +
x1 x2 x3 , x1 ∈ [10, 20], x2 ∈ [15, 25], x3 ∈ [5, 15]. The obtained results are presented in Table 1,
the frequency histogram are shown in Fig. 4.
668
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
Determine
sampling number
Particle swarm
initialization
Not
satisfy
Sample according to
the variable distributions
Not
satisfy
Random analysis of the
structure by Monte Carlo
Save the network
structure
Satisfy
The operation of
particle swarm (a)
Judge the accuracy
Stochastic sensitivity
analysis
Judge premature
and steady
Output the stucture
response
Y
Generate chaotic
sequences
Construct the network
of structure
Calculate the structure
response
The operation of
particle swarm (b)
Optimize the neural
network
Sample according to
the sampling strategy
Output the final
particles
Determine the statistical
parameters and distributions
Judge limit
value of iteration
N
Satisfy
Fig. 1: Stochastic sensitivity analysis flow char
Table 1: Analysis results of the structure response (example 1)
Statistics
Sample number
Mean value (mm)
Standard deviation
Compute time (min)
MC
10000
3778.9
1287.6
1.67
HNN
50
3780.9
1286.9
0.18
Relative error
/
0.053%
0.054%
/
5.1.2
Example 2
A frame structure whose beam and column are all I shaped cross section is acted upon by horizontal load P (Fig. 2). Assuming that the parameters of I shaped cross section and load are
random input variables with uniform distributions and the distributions of these parameters are
shown in Table 2, transverse displacement of the joint 3 is the output variable. The structure
stochastic is analyzed based on these parameters. The analysis results of the structure response
are showed in Table 3 and the frequency histogram is shown in Fig. 5.
Table 2: Distribution parameters of random variables (example 2)
—
Value interval
Beam height
Beam width
Beam flange
Beam web
Load
(mm)
(mm)
thickness (mm)
thickness (mm)
(KN)
200-300
120-180
8-12
8-12
14-20
669
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
1-1
P
2
3
10
10
250
4m
150
1-1
1
4
Fig. 2: Schematic of truss structure
Table 3: Analysis results of the structure response (example 2)
Statistics
Sample number
Mean value (mm)
Standard deviation
Compute time (min)
MC
10000
23.7
9.8
16.67
HNN
50
23.51
9.61
0.1
Relative error
/
0.8%
1.94%
/
5.1.3
Example 3
Joint is the basic component of frame structure. According to the mechanics of joints, all kinds
of joints have the characteristic of semi-rigid. The prerequisite of researching the semi-rigid joints
is to explore the impact of each component of joint. So a welded steel joint is utilized as the
model. The connection of joint is fully welded, the flange connection is slop butt weld and the
web connection is twin fillet welding seam (Fig. 3). Assume that the beam height, beam width,
beam flange thickness, beam web thickness and column flange thickness are input variables with
uniform distribution and the distributions of these parameters are shown in Table 4, the joint
rotation is output variable. The analysis results of the structure response are shown in Table 5
and the frequency histogram is shown in Fig. 6. By comparing the results of HNN with the results
of MC in Table 1, 3 and 5, it is evident that the results of HNN meet requirement of the accuracy,
while requiring much fewer model runs. Such as in example1, the relative error of the mean
value between MC and HNN is 0.053%, the relative error of the standard deviation is 0.054%; in
example 3, the compute time of HNN is 0.92 hour while the compute time of MC is 166.67 hour.
The frequency histograms of Fig. 4, Fig. 5 and Fig. 6 indicate that the approximations are close
to the true distributions, moreover HNN proves useful for the sensitivity analysis of complicated
models, the more complex the model, the more obvious advantages of the algorithm.
Table 4: Distribution parameters of random variables (example 3)
—
Value interval
Beam height
Beam width
Beam flange
Beam web
Load
(mm)
(mm)
thickness (mm)
thickness (mm)
(KN)
200-300
120-180
8-14
6-14
14-20
670
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
Table 5: Analysis results of the structure response (example 3)
Statistics
Sample number
Mean value (mm)
Standard deviation
Compute time (min)
MC
10000
5.06
0.68
166.7
HNN
50
5.08
0.72
0.92
Relative error
/
0.8%
1.94%
/
16
300
16
800
30 140 30
MC calculations
Probability frequency
700
5
Beam
5
8
Joint region
234
8
Column
HNN calculations
600
500
400
300
200
100
0
1000 2000 3000 4000 5000 6000 7000 8000 9000
V
Fig. 3: Schematic diagram of joint
Fig. 4: Frequency histogram (example 1)
700
800
MC calculations
MC calculations
600
HNN calculations
Probability frequency
Probability frequency
700
600
500
400
300
200
10
20
30
40
50
60
70
Displacements of the third joint (mm)
Fig. 5: Frequency histogram (example 2)
5.2
500
400
300
200
100
100
0
0
HNN calculations
80
0
3.5
4.0
4.5
5.0
5.5 6.0
Rotation (1e-4)
6.5
7.0
7.5
Fig. 6: Frequency histogram (example 3)
Stochastic Sensitivity
Calculate the sensitivity coefficients of example 2 and example 3 by the algorithm of section
the correlation coefficients between beam height and beam width are both 0.8 in example 2 and
example 3. The sensitivity results are shown in Fig. 7 and Fig. 8. BW is the beam width, BH is
the beam height, BFT is the beam flange thickness, BWT is the beam web thickness, CFT is the
column flange thickness, P is the load.
The sensitivity coefficients in Fig. 7 reflect the positive and negative correlation between pa-
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
0.4
0.4
0.2
0.2
0
0
−0.2
−0.2
−0.4
−0.4
−0.6
Sensitivity without considering
parameters correlation
Sensitivity considering parameters
correlation
−0.8
−1.0
BW
BH
BFT
BWT
P
Fig. 7: Sensitivity coefficients (example 2)
−0.6
Sensitivity without considering
parameters correlation
Sensitivity considering parameters
correlation
−0.8
−1.0
671
BH
BW
BFT
BWT
CFT
Fig. 8: Sensitivity coefficients (example 3)
rameters and structure response. When the parameters correlation is not considered, beam height
and beam web thickness are positively correlated with joint 3 displacement, that is to say, the
displacement of joint 3 is increased with an increase of beam height and web thickness. The
sensitivity coefficients considering parameters correlation are −0.7215, −0.5315, −0.334, −0.106
and 0.2722 for BW, BH, BFT, BWT and P in turn. However this conclusion is obviously consistent with the actual situation that all the parameters are negatively correlated with displacement
of joint 3 except the external load with considering parameters correlation. Fig. 8 reflects the
differences between considering correlation and without considering correlation. In the condition
of considering parameters correlation, BH is the most important sensitive parameter for rotation,
other sensitive parameters are BW, CFT, BWT and BFT in order of importance. The sensitivity
coefficient of BW is 0.2 under the condition of without considering parameters correlation, that
is mean the rotation is increased with increasing the beam width and beam flange thickness. It
is obvious false for without considering parameters correlation. However all the parameters are
negatively correlated with the rotation for considering parameters correlation is accordance with
the actual form mechanics.
6
Conclusions
This study suggests hybrid neural network as the approximate model and global sensitivity technique to stochastic sensitivity problems. The conclusions are as follows:
(1) The algorithm of chaos-particle swarm has been improved by optimizing the particles alternately between steady and chaos. The proposed technique that hybrid neural network was
optimized by the improved chaos-particle swarm provided a satisfactory approach for highly
non-linear implicit functions.
(2) Agreement appears much closer between these results which were calculated by hybrid neural
network and the Monte Carlo results. In addition, these results of approximate models could
be used for sensitivity analysis in probability design.
(3) Compared with Monte Carlo which would obtain complex structure responses by finite element, numerical analysis has demonstrated that the algorithm in this paper shows great
672
X. Dong et al. / Journal of Information & Computational Science 12:2 (2015) 657–672
superiority in computational efficiency.
(4) The effect of the random variable dimensions and measurement scales on structure sensitivity
can be removed by the combination of global sensitivity algorithm and stochastic analysis.
The parameter sensitivities are then obtained numerically using the global method. The calculations of the sensitivity coefficients considering parameters correlation are more consistent
with the actual.
References
[1]
J. Choi, J. W. Harvey, M. H. Conklin, Use of multi-parameter sensitivity analysis to determine
relative importance of factors influencing natural attenuation of mining contaminants [C], Proceedings of the Technical Meeting, Charleston, South Carolina, March, 1999, 8-12
[2]
M. Horsen, J. C. Refsgaard, S. Hansen, Assessment of uncertainty in simulation of nit rate leaching
to aquifers at catchment scale [J], Journal of Hydrology, 242(3), 2001, 210-227
[3]
Phelim Boyle, Mark Broadieb, Paul Glasserman, Monte Carlo methods for security pricing [J],
Journal of Economic Dynamics and Control, 21(8), 1997, 1267-1321
[4]
Daojin Lin, Quan Qin, Reliability stochastic finite element analysis of in-plane buckling of an
existing arch bridge [J], Engineering Mechanics, 22(6), 2005, 122-126
[5]
Ran Hu, Dianqing Li, Chuangbing Zhou et al., Structural reliability analysis using stochastic
response surface method [J], Engineering Mechanics, 27(9), 2010, 192-200
[6]
Lufeng Yang, Xianfeng Yang, Bo Yu et al., Improved stochastic response surface method based on
stepwise regression analysis [J], Chinese Journal of Computational Mechanics, 30(1), 2013, 88-93
[7]
Long Tang, Guangyao Li, Hu Wang, Kriging-HDMR meramodeling technique for nonlinear problems [J], Chinese Journal of Theoretical and Applied Mechanics, 4, 2011, 780-784
[8]
Hu Wang, Long Tang, G. Y. Li, Adaptive MLS-HDMR metamodeling techniques for high dimensional problems [J], Expert Systems with Applications, 38(11), 2011, 14117-14126
[9]
Chonggang Xu, Gertner, George Zdzislaw, Uncertainty and sensitivity analysis for models with
correlated parameters [J], Reliability Engineering and System Safety, 93(10), 2008, 1563-1573
[10] Y. Shi, R. C. Eberhart, Empirical study of particle swarm optimization [C], Proceedings of the
1999 IEEE Congress on Evolutionary Computation, CEC 99, 1999
[11] Fei Gao, Hengqing Tong, Parameter estimation for system based on particle swarm optimization
[J], Acta Phys., 55(2), 2006, 577-582
[12] Xiaobo Xu, Kangfeng Zheng, Dan Li et al., New chaos-particle swarm optimization algorithm [J],
Journal on Communication, 33(1), 2012, 24-29
[13] Kaitai Fang, Changxing Ma, Orthogonal Design and Uniform Design [M], Beijing: Science Press,
2001
[14] W. Clifford Hansen, C. Jon Helton, J. C´edric Sallaberry, Use of replicated Latin hypercube sampling to estimate sampling variance in uncertainty and sensitivity analysis results for the geologic
disposal of radioactive waste [J], Procedia-Social and Behavioral Sciences, 2(6), 2010, 7674-7675
[15] A. Der Kiureghian, P. L. Liu, Structural reliability under incomplete probability information [J],
Journal of Engineering Mechanics, 112(1), 1986, 85-104