An Efficient and Time-Saving Reliability Analysis Strategy for Complex Mechanical Structure

With the development of production technology, the mechanical structure has become more and more complicated, which makes the simulation process of the corresponding computer model very time-consuming. As a result, the reliability analysis needs to consume huge time resources. To deal with this problem, an improved method, which balances the accuracy of failure probability and the efficiency of computation, is proposed. The novelty of the proposed method is the use of an efficient point selection strategy, the $k$ -means algorithm, and the constraint of correlation among the training points. The $k$ -means algorithm can divide the candidate points into a few groups. Therefore, we can update the Kriging model by selecting several sample points which have large contributions to improve the accuracy of failure probability in each iteration. Meanwhile, a constraint is applied to control the relative location among points of DoE to avoid redundant information. The efficiency and accuracy of the proposed method are verified through two numerical examples. Finally, a type of artillery coordination system as a representative of the complex mechanisms is mentioned. And the proposed method is applied to calculate the reliability of the position accuracy of the coordination process, which proves the significance of the proposed method in engineering practice.


I. INTRODUCTION
For the reliability analysis of an engineering structure, the m-dimension random variables x and the performance function G(x) are represented as x = [x 1 , x 2 , . . . , x m ] T and G(x) = G(x 1 , x 2 , . . . , x m ) respectively. The performance function G(x) = 0 divides the variable domain into two parts: the failure domain G(x) ≤ 0 and the safe domain G(x) > 0. The failure probability can be defined as follows: where f (·) is the joint probability density function of x.
To obtain the value of failure probability, multiple methods have been proposed. The first order and second order reliability methods (FORM and SORM) [1] evaluate the failure probability through the most probable point (MPP).
The associate editor coordinating the review of this manuscript and approving it for publication was Yu Liu . However, the accuracy can hardly be ensured while dealing with non-linear and high dimensional problems. Monte Carlo Simulation (MCS) [2], [3], as a sampling method, is accurate and widely used. But it is extremely time-consuming which needs to call a large number of the performance function, especially for the problems with low failure probability. To improve the efficiency, scholars propose some methods to reduce the variance of MCS, such as line sampling [4], [5], importance sampling [6]- [8], and subset simulation [9], [10]. However, these improved methods based on MCS also cannot satisfy the computational requirements.
Nowadays, surrogate models are widely used in the structural reliability analysis. Such methods replace the performance function with a surrogate model which can improve the computational efficiency evidently. The surrogate models include Support Vector Machine [11], [12], Response Surface [13], [14] {Liu, 2016 #1105;Rusek, 2016 #1106}, Neural Networks [15], Polynomial Chaos Expansions [16], Kriging [17]- [19], and so on. VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ As is known, the Kriging model, which has been used in optimization problems and structural reliability analysis, is constructed based on the design of experiments (DoE).
To build the Kriging model in an efficient way, numbers of adaptive methods have been proposed in the last few years [20]- [22]. The efficient global reliability analysis (EGRA) method [20], [23], which is inspired by efficient global optimization (EGO) method [24], uses the expected feasibility function (EFF) to select the new points. According to EFF, the points which have high variance near the region of the limit state are selected to update the Kriging model. But sometimes the points with low probability density may be select, which have little influence on the failure probability. Basing on the Kriging model and MCS, Echard et al. [21] propose an active learning method (AK-MCS). Zhifu Zhu et al. propose a new criterion of selecting points based on the full information of the Kriging model [25]. Sun et al. [26] update the Kriging model through the points selected by function LIF. To solve the problems with low failure probability, the AK-IS method which replaces MCS with the importance sampling method is proposed [27]. Lv et al. [28] propose an efficient learning method (AK-LS) with a new function H which is introduced from the information entropy theory. Inspired by the AK-IS method, ALK-KDE-IS is proposed to solve the low failure probability problems by using kerneldensity-estimation [29]. Morató et al. [30] apply the Kriging model on turbine support structures. Zhang et al. [31] propose a new projection-outline-based active learning strategy to improve the efficiency of analysis. Chai et al. [17] propose a new learning function L f and apply it to fatigue crack reliability. Yun et al. [32] propose the step-wise AK-MCS method to solve the problems with fuzzy failure probability. Linxiong et al. [33] use the improved EGO method and secondary point selection strategy to improve the efficiency of reliability analysis. Xu et al. [34] improve the AK-MCS method based on the modified subset simulation to solve the problems with low failure probability. Wang and Sun [35] propose a stepwise accuracy improvement strategy of Kriging (SAIS). The results in the literature show that the SAIS method can obtain an accurate result with fewer samples.
Nevertheless, there are still some drawbacks in the SAIS method. First, it can only add one sample points to the DoE to update the current Kriging model in each iteration. It means that in one iteration, the accuracy of the Kriging model can be improved only in a small region near the added point which leads to inefficiency. Second, there is no measure to limit the relative location among the added points to avoid redundant information and useless performance function evaluation. Aiming at solving the engineering problem, this article proposes an efficient and time-saving Kriging reliability analysis strategy by using parallel computation and considering the correlation of the added points in DoE. The new method is suitable for solving implicit high-dimensional problems with long simulation time and has great engineering value.
The paper is organized as follows: Section II provides an introduction to the Kriging model and the SAIS method.
Section III proposes a new method to reduce the iterations and improve the efficiency of computation. Then two examples are studied to illustrate the high efficiency of the new method. In Section IV, a type of artillery coordination system with flexible structure is analyzed. Section V presents the conclusion.

A. THE KRIGING MODEL
The Kriging model, a linear unbiased predictor model, includes two parts, one is a deterministic part and the other is a stationary Gaussian process. The performance function is presented as follows: where f(x) is the basis function, and β is the corresponding regression coefficient. The covariance of z(x j ) and z(x k ) is expressed as follows: Cov where σ 2 is the variance of the random process, R(x j , x k ; θ) is the correlation function of z(x j ) and z(x k ), θ is the parameter vector of R(x j , x k ; θ). In this article, the correlation function is represented as.
where x i j and x i k are the i th components of x j and x k . θ i is the i th components of θ. The parameter vector θ can be calculated by maximum likelihood estimation. More details are introduced in [21].
According to the theory of the Gauss process, the value of performance function at point x follows a normal distribution.

B. THE SAIS METHOD 1) THE ERROR OF THE FAILURE PROBABILITY
For the Kriging model, the uncertainty of the sign of G(x) can be estimated as follows: where U (x) is proposed in [21].
The SAIS method calculates the absolute error betweenP f and P f by Eq. (7).
The expectation of |P f − P f | is presented as where E N is proposed in [26] and used as an accuracy measure in [35].
C. THE METHOD OF SELECTING THE BEST NEXT POINT The current Kriging model is constructed with [x 1 , x 2 , · · ·, x N ] and [y 1 , y 2 , · · ·, y N ] T . If a given point x N +1 with y N +1 is added into the Kriging model, the accuracy measure E N will be improved to E N +1 , (10) and the improvement and its expectation are presented as, (12) A large value of (x N +1 ) means that the point x N +1 has a high contribution to improving the accuracy of failure probability. Therefore, the best new point can be select as follows:

III. THE IMPROVED NEW METHOD
Currently, the two main challenges of engineering problems are the accuracy of estimated failure probability and the efficiency of calculation. Thus, how to improve the SAIS method to reduce iterations and save the time-consuming, is considered in this section.

A. THE STRATEGY OF SELECTING BEST POINTS 1) APPLICATION OF CORRELATION
Now, an example with two-dimensional variables proposed in [21] is studied. Using the SAIS method, when theP f is accurate enough, the number of iterations is 48. Fig. 1 shows the training points of Kriging. We can find that some points are too close to others, presented as red squares. This phenomenon is not helpful to improve the accuracy of the Kriging model. To deal with this issue, a hypothesis is formed: When we add x N +1 to the current DoE with y N +1 and construct a new Kriging model, the accuracy of the Kriging model in a region near the x N +1 is high enough and it is unnecessary to select another training point in this region.
Thus, how to determine the size of the region is a key. A simple method is to control the distance between added points (||x j − x k || 2 ) directly. But different variables have different influences on the accuracy of the Kriging model, that is, the above constraint is rough. Considering Eq. (4), the information of relative position among points can be presented as the correlation function R(x j , x k ; θ), which weights the effects of different variables. The correlation condition is defined as where x new is the new added point,  [36].

2) APPLICATION OF PARALLEL COMPUTATION
Now, consider the problem of parallelizability. The black dotted line represents the actual function and the red solid line means the fitting curve, the black stars represent the best points selected by SAIS. In Fig. 2, we can found that from the 18 th iteration to the 22 nd iteration, four best new points are added and the Kriging model is updated as well. However, these points are all located in a small region represented by the ellipse. According to the correlation function R(x i , x j ; θ), only a small part of the Kriging model which belongs to the vicinity of the added points can be improved. And the added points have little effect on the accuracy of other regions. It is noteworthy that when the accuracy of the Kriging model is rough, the best new points for successive iterations may lie in a small region. Until the accuracy of this region has VOLUME 8, 2020 been improved and is higher than other regions, the new best point can be selected from other more inaccuracy regions. For efficiency, this situation needs to be avoided. Aiming at updating the Kriging model from multiple regions simultaneously and improving the efficiency, parallel computation is used. For complex engineering problems, the parallel operation can greatly reduce the simulation time by performing multiple simulations at the same time.
As shown in Fig. 3, yellow triangles mean the candidate points, which are used to update the Kriging model. The candidate points are generated by Monte Carlo Markov Chain (MCMC). As an example, in order to improve the efficiency, we select 4 best next points directly in one iteration based on the value of (x N +1 ). The defect can be found in Fig. 3. The black stars B and C, which are close to A, have little influence on improving the accuracy of the Kriging model. This condition is not wanted and how to avoid it is a key.
A direct method is dividing the variable space into some sub-spaces, and then choose one best point from each subspace. However, it is difficult to find a suitable space divided strategy for this issue and the candidate points generated by MCMC does not cover the whole space. In this article, the points divided method will be substituted for space divided method. A widely used point divided method is the k-means algorithm. Firstly, generate candidate points by the MCMC method. Secondly, the k-means algorithm is applied to divide the points into k groups. Then select the best point from each group. In this way, k best points (x 1 best , x 2 best , . . . , x k best ) can be selected in one iteration. In Fig. 4, the candidate points are divided into four groups by the k-means algorithm, presented by different colors. Meanwhile, we can consider that the fitting curve of the Kriging model and the x-space consist of four parts corresponding four colors, named as P 1 , P 2 , P 3 , P 4 . We can select the best next points from each part. Now, discuss how to select the sample points correctly. Some methods also use parallel computation to improve the efficiency [29,37]. These methods select the best points by using learning functions at the same time without considering the correlation condition between each other, which may lead to the situation shown in Fig. 5. To avoid this flaw, the constraint mentioned above is used again. The best points will be selected one by one instead of at the same time.

B. APPROXIMATION OF THE BEST POINTS
As is known, the further the points distance from the Kriging model, the lower effect on the accuracy of the Kriging model. In other words, the points of one part cannot remarkably influence the accuracy of the Kriging model belonging to other parts. Aiming at improving the efficiency of computation and simplifying calculation, the weak influences of points from other parts are ignored. So the new strategy of selecting the best points is defined as follows, The number of parts divided by the k-means method is k, P i means the i th part (i ∈ {1, 2, . . . , k}). For the i th part, the candidate points which cannot satisfy Eq. (15) are eliminated, and the best point will be select from the remaining points. Assume that the current Kriging model is constructed with [x 1 , x 2 , · · ·, x N ] and [y 1 , y 2 , · · ·, y N ] T . For a given point x i N +1 with y i N +1 is addition to the Kriging model, the accuracy of the pseudo Kriging model which belongs i th part is represented as E i N +1 .
The improvement is represented as follows. (17) Without calling the performance function, the true value of y i N +1 is unknown, and the value of E i (x i N +1 , y i N +1 ) cannot be calculated. But according to Eq.TABLE 3, the distribution information of y i N +1 is known. Basing on the distribution of y i N +1 , the expectation (x i N +1 ) can be represented as Eq. (18). And the best point is defined as Eq. (19).
To simplify calculation, the sign of G(x) with U N (x) > 2 can be considered as correctly predicted. A new symbol Ẽ i (x i N +1 , y i N +1 ) is defined and the conditional PDF is used (20), as shown at the bottom of the next page. The random samples with f (x|U N (x) < 2) are generated by MCMC. The integration of E i (x i N +1 , y i N +1 ) can be expressed as Eq. (21).
where x i M ,1 , x i M ,2 , . . . , x i M ,S i are the random samples belongs to P i . S i represents the number of samples in P i . S is the total number of samples generated by MCMC. The value of S is discussed in [35], the increasing of S cannot remarkably influence the iterations. In this article, According to Eq.(19)∼(21), the best point x i best which has large contribution to improving the accuracy of Kriging model can be represented as.
The integral of Eq. (22) can be calculated by Gauss-Hermite quadrature.
where N G is the number of quadrature points, v and w are defined as quadrature points and weights.

C. THE PROCEDURES OF THE IMPROVED METHOD
Step 1: Generate the initial DoE and y. Select the initial DoE [x 1 , x 2 , · · ·, x N ] by using Latin hypercube sampling (LHS), and get the values of y = [y 1 , y 2 , · · ·, y N ] T by calling the performance function.
Step 2: Generate N MCS random points with the MCS.
Step 3: Construct the initial Kriging model with the DoE and y. Use the DACE, a MATLAB tool box, to construct the initial Kriging modelĜ N (x), (t = 0).
Step 4: Generate the candidate points by MCMC in the region with U (x) ≤ 2. And divide the candidate points into k groups by using k-means algorithm.
Step 5: Select the best next points and update the current DoE. Select the best next points from each group x best = [x 1 best , x 2 best , . . . , x k best ] by the method referred above. And get the corresponding performance values y = [y 1 best , y 2 best , . . . , y k best ] T with the performance function. Finally add the k points into the DoE.
Step 7: Calculate the value ofP f ,N +t and E N +t with the N MCS random points. If the stopping criterion is satisfied, turn to step 8, else turn back to step 4.
The stopping criterion of the improved method is presented as Eq.(25) where Step 8: Calculate the coefficient of variation ofP f According to Eq. (27), the N MCS random samples generated by MCS must be large enough to ensure δ P f ≤ 0.03. If not, we should generate more random samples, and turn to step 4. The procedure of the proposed method is also illustrated in Fig. 6. VOLUME 8, 2020 The performance function of this example, a series system with four parts, is expressed as follows [21], [38], where x 1 and x 2 are mutually independent variables: x 1 ∼ N (0, 1), x 2 ∼ N (0, 1). Analyze this example by the improved method. Fig. 7 shows that the use of correlation [r] can avoid the phenomenon showed in Fig. 1 successfully. To reduce the influence of different DoE, Monte Carlo samples, and MCMC samples on the results, the average results of running 10 times for different methods are summarized in Table 1. And Fig.8 and Fig.9 show the results of different methods for running one time. The number of points in the initial DoE is 6. N call means the number of calling the perform function. N it is the number of iterations. ε is the relative error between P f andP f . SAIS-R means improving the SAIS method only considering the correlation of added points in DoE. The black line in Fig. 8(a) expresses the failure probability gained by MCS. And the black line in Fig. 8(b) is the threshold of error (where, [e ε ] = 0.01).
From Fig. 8 and Table 1, it can be seen that, with the increasing of N it , the proposed method is more efficient than SAIS, and the value ofP f calculated by the improved method converges faster than other methods. WhenP f is convergent, the improved method only needs less than 20 iterations. Comparing different values of k, it needs 20.6 iterations to be convergent when k = 2 and when k = 6, it needs only 6.7 iterations. Now, discuss the effect of different values of k. Fig. 9(a) and Fig. 9(b) present the estimate of failure probabilityP f and the e ε corresponding different k respectively. It can be found that the increasing of k does not remarkably influences the number of calling the performance function.

2) A MODIFIED RASTRIGIN FUNCTION
The perform function, a complicated non-linear function with two variables, is presented as follows [35], where x 1 , x 2 are standard normal distributed random variables.
In this example N 0 =6, [e ε ] = 0.05. The criterion reliability result is calculated by MCS with N MC =10 5 . The average results of running 10 times for different methods are summarized in Table 2. Fig.10 and Fig.12 show the results of different methods for running one time.
In Table 2, Compared with other methods, we can find that the new method is more efficient. In Fig. 11, with the same   N call , the SAIS-R method fits the curve of G(x) = 0 better than SAIS, which confirms the feasibility of the constraint based on correlation function. In Fig. 12, the influence of different values of k is discussed. It is found that the value of k cannot significantly influence the convergence of failure probability and relative error.

IV. RELIABILITY ANALYSIS OF A TYPE OF ARTILLERY COORDINATOR
With the development of technology and industry, the mechanical mechanisms have become more and more complex. Scholars found it hard to apply the traditional classical kinematics and dynamics theory to perform the structural design and optimization analysis of these complicated mechanisms, then they put forward the multi-body dynamics theory. The multi-body system contains three types: the rigid multibody system, rigid-flexible coupling multi-body system, and multi-flexible body system. In recent years, scholars have begun to study the multi-body dynamics of the mechanism with flexible bodies. Zhang et al. [39] establish the dynamic analysis model of the parallel mechanism with flexible body components. Olivier et al. [40] uses the ANSYS software to establish the dynamic model of the multi-flexible body and applies it to the analysis of the semi-active suspension of the automobile and the manipulator. Taking the port crane as an example, Zhang et al. [41] put forward a rigid-flexible coupling modeling method, which simplifies the flexible model and improves the efficiency of simulation significantly. He et al. [42] build a rigid-flexible coupling model to analyze the vibration of train bridge. In this section, we establish a rigid-flexible coupling model of a type of artillery coordinator, and analyze the reliability of the position accuracy.
The structure of the coordinator mainly includes: trunnion, reducer, coordination arm, hydraulic cylinder, pendulum mechanism, and balancing machine. The coordination arm is welded by section steel and mounted on the trunnion to rotate around the trunnion. The coordination arm is welded with rotary support of the shell tray and lower support of the oil cylinder. The shell tray is mounted on the upper left of the coordination arm and is connected to the coordination arm by a pin shaft. The shell tray is welded by a semicircular plate and a trapezoidal steel plate in which the shell is stored.
When it works, the coordinator catches the shell pushed out by the ammunition supply silo and coordinates it to a specified angle. Next, under the action of the overturning cylinder, the shell is quickly flipped to the ammunition ramming line.
Then the shell is transported to the gun bore. The diagram of the coordination process is shown in Fig. 13. The coordinator is considered to be failed if the coordinator cannot transport the shell to the specified location within a given time. The main factors that influence the motion reliability of the coordination mechanism are as follows: (1) The size error and assembly error will influence the output precision of the coordination and further influence the position accuracy of the coordination.
(2) The output error of the driving components (motor, hydraulic, etc.). If the hydraulic cylinder cannot drive the telescopic rod to the theoretical position at the specified time point, the output of the mechanism will deviate from the theoretical value, which may influence the position accuracy of the mechanism.
(3) The deformation of components. Under the condition of high speed and heavy load, the structure will be deformed, and the position accuracy of the structure will inevitably be reduced.
Considering the above factors, a rigid-flexible coupling parametric model of artillery coordinator is established to fit the actual situation. Fig. 14 is the main flow of establishing the rigid-flexible coupling model of the coordinator based on ADAMS and ANSYS software.   15 shows the parametric model of rigid-flexible coupling for the coordinator. The pendulum mechanism shown in red is the flexible body part. Fig. 16 is the displacement curve of the center of shell in the standard state. The whole motion process is 1.0 seconds. The displacement of the center of the shell at t = 1.0s is collected and used as the standard value ξ . The performance function of the position accuracy is defined as,  where, γ is permissible error, which can be defined according to the position accuracy requirement (in this article γ = 3.0mm); r 1 -r 4 represent the radius of the hole; l x and l y represent the error in the x, y direction of the center of shell caused by the installation and size error. l 1 , l 2 represent the displacement of balancing machine and hydraulic cylinder; ξ (r 1 , r 2 , r 3 , r 4 , l x , l y , l 1 , l 2 ) is the true displacement at 1.0s obtained by simulation based on the error factors. Simulink-ADAMS joint simulation (Fig. 17) is used for reliability analysis. AK-MCS, LIF, SAIS, and the improved method with k = 4 are used to perform the reliability analysis. According to the complexity of the model and  the computer performance, each simulation time is 28min. The number of initial DoE is 20. The result can be found in Fig. 18 and Table 4.
From Fig. 18 and Table 4, we can found that the result obtained by the four methods are basically the same, which proves the accuracy of the improved method. Comparing the N it and total time, the improved method can greatly reduce the reliability analysis time from more than 100 hours to only 23.8 hours, which has important application value.

V. CONCLUSION
This article aims at improving the efficiency of evaluating the value of failure probability by reducing the time of calculation. An improved method, based on the SAIS method, is proposed. The advantages of the new method are as follows: (1) To improve the efficiency of computation, the proposed method can select multiple new points, which have high contributions on the accuracy of the failure probability, in each iteration to update the current Kriging model from multiple regions.
(2) The simulations of sample points are carried out by parallel computation which can perform multiple simulations in one simulation time.
(3) Aiming at simplifying calculation and saving time, we improve the equations for selecting the best next points in [35].
(4) To avoid redundant information, the correlation function is used to limit the relative location among the points in DoE.
The new method is verified by two examples. Comparing with other methods, the improved method needs fewer iterations to be convergent. And increasing the values of k can reduce the number of iterations significantly. Finally, an artillery coordination system is introduced. And the new method is used for the reliability analysis of position accuracy. It can be found that for complex and time-consuming problems, this method has a significant reduction of the reliability analysis time, which has great engineering value.
Because the proposed method calculates the failure probability by MCS, a large number of sample points are required when solving the reliability problems with small failure probability. According to Eq. (28), if P f = 10 −7 ∼ 10 −5 , it needs about 10 8 ∼10 10 sample points, which expends a large amount of computer memory. In this case, it will cost long time in each iteration. To avoid above problem, we can replace MCS with the IS method. In this way, we can greatly reduce the number of samples and improve the efficiency of calculation.