Enhanced dung beetle optimizer for Kriging-assisted time-varying reliability analysis

1.   Introduction

To improve the effective solution of complex engineering problems, obtaining more accurate reliability analysis results is a long-term task ^[1,2,3,4,5]. A variety of reliability analysis methods and optimization methods are used to solve the complexity of real-world problems. Among them, time-varying reliability analysis (TRA), as a method to evaluate the failure probability of the system during the design life, is widely studied in the research and practical application of engineering systems ^[6,7]. There are two methods to solve the TRA problem: the outcrossing rate method and extreme value method ^[8]. However, these TRA methods face the problems of complexity and high computational cost ^[9,10].

Many scholars have proposed a TRA method based on the surrogate model to solve the computational difficulty faced by traditional analysis methods ^[11,12,13]. A substitute model is used to replace complex performance evaluation, which reduces the computational pressure of the TRA method that requires a lot of computational cost ^{[14,15,16,17]}. Common surrogate models include the physical information neural network (PINN) ^[18,20], the Kriging model ^[21,22], and the Gaussian model ^[23,24]. Dang et al. proposed a TRA method that combines PINN and Bayesian theory. Compared with Monte Carlo simulation (MCS), which only uses finite element analysis, the reliability analysis method combined with a surrogate model greatly reduces the computational difficulty ^[25]. Chen et al. proposed a TRA method assisted by the Kriging model. Due to the adoption of the Kriging model building method which actively tracks failure samples, the efficiency of the new TRA method has a high improvement effect ^[26]. Wang et al. proposed a TRA method based on Gaussian processes. This study improves the efficiency of surrogate modeling by combining machine learning algorithms and the Gaussian regression process ^[27]. Relatively, the Kriging model can provide both the predicted value and the predicted variance, which makes the adaptive Kriging model receive higher attention ^[28,29,30]. With the aid of the surrogate model, the efficiency of TRA is improved by higher efficiency. However, in the face of the problem of determining the specific time of prediction error, the efficiency of TRA assisted by the surrogate model has a lot of room for improvement. This is an optimization problem, so more excellent optimization algorithms are needed to improve the accuracy and efficiency of TRA methods ^[31,32]. Compared with the gradient-based optimization algorithm, the heuristic algorithm has higher applicability and global convergence ^[33]. Therefore, it is of great significance to develop a better heuristic optimization algorithm for TRA in engineering applications.

The group optimization technique is an important heuristic algorithm, which develops excellent algorithms from the coordinated behavior observed in natural organisms as an inspiration ^[34,35,36]. The optimization process of group optimization technology is divided into two stages ^[37,38,39]. The first stage is to generate a set of random individuals in the search space. The second stage is to combine, change, and evolve these random individuals repeatedly. Common swarm optimization algorithms use biological principles to push the population to the optimal or near-optimal scheme, which is the main difference between different optimization algorithms ^[40,41]. These algorithms include genetic algorithms (GA) ^[42], whale optimization (WOA) ^[43], particle swarm optimization (PSO) ^[44,45,46], differential evolution (DE) ^[35,47], and grey wolf optimizer (GWO) ^[48] to mention some. In the face of complex optimization problems, traditional optimization algorithms often still face the problem of low computational efficiency and accuracy ^[37,49]. Although these optimization algorithms have since been extensively studied and optimized, the practicability of new optimization algorithms is still challenged due to the complexity of the problem ^{[50,51,52,53,54,55,56]}. By observing how dung beetles navigate and move dung balls over long distances, the developed dung beetle optimization (DBO) algorithm has a good balance between exploration and development ^[31]. Therefore, the DBO algorithm is widely studied and applied in research and engineering. However, DBO algorithms can fall into local optimization traps when dealing with some optimization problems, so there is potential for further improvement.

To solve this problem, this study proposes an enhanced DBO (EDBO) algorithm. Different from references ^{[31,42,43,44,45,46,47,48]}, the EDBO algorithm uses chaos to generate more uniform population samples. This makes the new heuristic algorithm easier to get out of the dilemma in the face of local optimization traps. In addition, this study also uses the Gaussian random walk strategy to update the position of the overall population, which improves the global convergence ability of the algorithm. Different from references ^[25,26,27], this study introduces the heuristic algorithm into the TRA method based on the adaptive Kriging model. Different from reference ^[57], the EDBO algorithm is used to linearize the limit state function in this study. With the help of the EDBO algorithm, the new TRA method has higher accuracy and computational efficiency. Different from reference ^[9], this study uses the EDBO algorithm for TRA. Compared with the TRA based on the differential optimization DE algorithm, the proposed EDBO algorithm-assisted reliability method has higher computational efficiency. Compared with the previous TRA methods, the proposed TRA framework has higher computational efficiency and analysis performance with the help of the EDBO algorithm and adaptive Kriging model establishment strategy. In the classical CEC 2005 (IEEE Congress on Evolutionary Computation) and CEC 2021 benchmark tests, the EDBO algorithm is compared with five popular heuristic optimization algorithms, and the superiority of the EDBO algorithm is proved. In addition, EDBO is used to maximize the prediction error of the adaptive Kriging method in TRA, which further improves the efficiency and accuracy of TRA. Finally, the effectiveness and superiority of the proposed TRA method are verified by four classical engineering problems.

The remainder of the paper is organized as follows. The general formulation of the adaptive Kriging-based TRA and the EDBO is thoroughly reviewed in Section 2. Then, in Section 3, the proposed EDBO algorithm is explained in detail. In Section 4 the standard real-world benchmark functions are used to validate the performance of the proposed method. In Sections 5 and 6, EDBO-based TRA and four widely used TRA test cases are used to validate the performance of the proposed method. Finally, Section 7 provides a summary and prospect of the research findings.

2.   Review on TRA and DBO

This section briefly reviews some preliminaries, basic definitions, and procedures of adaptive Kriging-based TRA and DBO.

2.1.   Adaptive Kriging-based TRA

One of the primary steps in this approach is the implementation of accurate surrogate models ^[58]. This is particularly crucial because the most probable point trajectory (MPPT), denoted as ${\widehat{u}}_{Mpp}\left(t\right)$, is solely varying on a single input variable, which is the time parameter $t$. This characteristic simplifies the process of constructing the one-dimensional response surface using gaussian process regression (GPR) ^{[57,59,60,61,62,63]}. By leveraging the capabilities of GPR, we can efficiently model the behavior of the MPPT over time ^[64].

Another key aspect of this methodology is the strategic use of the adaptive sampling approach. This approach is designed to focus exclusively on conducting MPP searches at critical time instances that have been carefully identified ^{[65,66,67,68]}. These instances are chosen to maximize the mean square error (MSE) of the GPR model ${\widehat{u}}_{Mpp}\left(t\right)$. By concentrating our computational resources on these critical time points, we substantially reduce in the overall number of MPP searches compared to the more traditional time-dependent reliability analysis (TDTRA) method ^[69,70]. This reduction in computational effort translates to significant gains in efficiency without compromising the accuracy of our results.

Utilize a Gaussian process regression to estimate the trajectory of the MPP in a time-varying scenario. Let ${t}_{i}$ be an arbitrary time moment within the interval $\left[{t}_{s}, {t}_{e}\right]$. Represent the MPP of the instantaneous limit-state function $g\left(X, Y\left({t}_{i}\right), {t}_{i}\right) = 0$ as ${\widehat{u}}_{Mpp}\left(t\right)$. As ${t}_{i}$ varies from ${t}_{s}$ to ${t}_{e}$, ${\widehat{u}}_{Mpp}\left(t\right)$ transitions from ${\widehat{u}}_{Mpp}\left(t\right)$ to${\widehat{u}}_{Mpp}\left(t\right)$. When these MPPs are interconnected, they form a trajectory in the U-space known as the MPPT, denoted as ${u}_{Mpp}\left(t\right)$ for $t\in \left[{t}_{s}, {t}_{e}\right]$. Figure 1 illustrates the MPPT with solid curves representing the limit-state boundaries at different time instances and a dashed curve representing the MPPT.

However, for a general time-varying performance function, obtaining ${u}_{Mpp}\left(t\right)$ through MPP searches at multiple discretized time instances, as in the TDTRA method, becomes computationally impractical ^[71]. To alleviate computational expenses, surrogate modeling is used to approximate the MPPT. The Gaussian process regression is chosen due to its significant advantage in quantifying prediction variance and its successful applications in the field of reliability analysis.

To determine the precise time ${t}^{*}$ at which the most probable point occurs, it becomes imperative to employ heuristic optimization techniques, as the underlying problem exhibits a non-convex nature. Non-convex optimization involves intricate and mathematically undefined relationships. This research introduces and scrutinizes various non-convex optimization approaches. In this thesis, we leverage a recent hybrid method that combines an enhanced DBO method to enhance the accuracy of ${t}^{*}$, consequently exerting an influence on both the MPP and the probability of failure. Additionally, we provide an exposition of other optimization methods that have been previously utilized for the purpose of comprehensive comparison.

The Kriging meta-modeling technique, also called GPR, is a widely employed statistical model in engineering design and optimization ^[72]. It proves especially valuable in scenarios where computational resources are limited. When dealing with costly simulation-based functions, Kriging is proficient at constructing a mathematical model based on a finite set of observed points, serving as an effective approximation. This cost-effective Kriging model aids in gaining insights into the target function by enabling exploration of the design space.

An efficient adaptive sampling approach should encompass three key components:

(1) Local Exploitation: This aspect entails modifying the sampling procedure to concentrate on regions with significant prediction errors. The local exploitation aspect can be represented in several ways, including by using cross-validation errors, prediction errors at validation sites, local optima locations, and local response variation analysis, among other techniques.

(2) Global Exploration: This element ensures that the global meta-modeling objective is maintained, preventing the oversight of potentially significant regions. It achieves this by employing criteria based on distances.

(3) Balancing Local and Global Considerations: This final component is pivotal in striking a balance between local exploitation and global exploration, significantly influencing sampling performance. However, many existing adaptive sampling approaches for Kriging tend to rely on fixed balancing rules, which may not be optimal for sampling performance.

Therefore, the maximizing error prediction estimation (MEPE) method for adaptive Kriging ^[57] is used for constructing surrogate models for the performance functions of mechanical structures.

2.2.   DBO

There exist numerous species of dung beetles, including Copris ochus Motschulsky, Onthophagus gibbulus, and Caccobius jessoensis Harold, among others ^[31]. These insects, widely distributed in nature, are recognized for their diet primarily consisting of animal dung. Dung beetles play a crucial role as decomposers in ecosystems worldwide. Notably, research indicates an intriguing behavior among dung beetles, wherein they form dung into spherical balls and roll them, as depicted in Figure 2. Notably, dung beetles aim to transport their dung balls swiftly and efficiently, thus avoiding competition from other dung beetles.

As illustrated in Figure 2, a dung beetle is observed rolling a dung ball backward, larger than its size. Conversely, an intriguing behavior of dung beetles involves utilizing celestial cues, particularly from sources like the sun, the moon, and polarized light, for navigation, ensuring the dung ball's movement along a straight trajectory.

However, in conditions of complete darkness where no light sources are available, the path of the dung beetle becomes curved and occasionally circular instead of straight. It's important to note that various natural factors, such as wind and uneven terrain, can cause dung beetles to deviate from their intended course. Moreover, dung beetles may encounter obstacles hindering their progress in the rolling process. Consequently, dung beetles often ascend onto the dung ball to perform a series of movements, including rotations and pauses, known as a "dance" which determines their subsequent direction of movement.

To mimic the behavior of rolling a ball, dung beetles need to maintain a consistent direction of movement across the entire search area. The trajectory of a dung beetle's movement is depicted in Figure 3. In this illustration, the dung beetle utilizes the sun for navigation, with the red arrow denoting the rolling direction. As the rolling progresses, the position of the dung beetle, engaged in rolling the ball, is continuously adjusted and can be formulated as follows:

where in the context of an iterative process denotes the current iteration. The position of the ${i}^{th}$dung beetle at the ${t}^{th}$ iteration is represented by ${x}_{i}\left(t\right)$. Here, the constant $k$ falls within the range $\left(0, \left.0.2\right]\right.$ and signifies the deflection coefficient. Additionally, b represents another constant within the range $\left(\mathrm{0, 1}\right)$. The natural coefficient, denoted by $\alpha$, takes values of either -1 or 1. ${x}^{w}$ signifies the global worst position, while $\mathit{\Delta} x$ mimics changes in light intensity.

Alpha denotes natural factors such as wind and uneven ground that can deflect dung beetles from their original path. A value of 1 for alpha means no deviation, while -1 indicates a deviation from the original direction. Alpha is randomly set to 1 or -1 to mimic real-world complexities. Similarly, a higher $\mathit{\Delta} x$ value signifies a weaker light source. Additionally, $k$ and $b$ are fixed at 0.1 and 0.3, respectively, according to ^[31].

The dung beetle, when obstructed, adjusts its direction by engaging in a dance-like behavior crucial for ball-rolling dung beetles, wherein the tangent function within the range $\left[0, \pi \right]$ is employed to derive the new rolling direction, after which the beetle rolls the ball backward, updating its position accordingly and defined as Eq (2).

where $\theta$ is the deflection angle $\left[0, \pi \right]$.

Dung beetles roll dung balls to secure locations for hiding, vital for safeguarding their offspring, thus motivating the development of a boundary selection strategy to mimic suitable egg-laying areas based on natural behavior, which is shown by Eq (3).

where ${X}^{*}$ represents the current best position within a limited spawning area defined by $L{b}^{*}$and $U{b}^{*}$ while $R = 1-\frac{t}{{T}_{max}}$ with ${T}_{max}$ as the maximum iteration number progresses. $Lb$ and $Ub$ signify the lower and upper bounds of the optimization problem, respectively.

Female dung beetles select brood balls within the identified spawning area for egg laying, with the boundary range dynamically changing based on the $R$ value, as indicated by Eq (3), resulting in the dynamic positioning of the brood ball throughout the iteration process.

where the position of the ${i}^{th}$ brood ball at the ${t}^{th}$ iteration is denoted as ${B}_{i}\left(t\right)$, with ${b}_{1}$ and ${b}_{2}$ representing two separate random vectors sized $1\times D$, where $D$ signifies the dimension of the optimization problem; it's important to highlight that the brood ball's position is confined strictly within a defined range, namely, the spawning area.

Upon reaching adulthood, certain dung beetles surface from the ground in search of food, prompting the establishment of an optimal foraging area to direct their foraging activities, mimicking their natural behavior, with the boundaries of this area precisely defined by the Eq (5) below.

where ${X}^{b}$ represents the best global position, while $L{b}^{b}$ indicates the lower and upper bounds of the optimal foraging area, with other parameters defined in Eq (3), leading to the subsequent update of the small dung beetle's position as shown in Eq (6).

where ${x}_{i}\left(t\right)$ represents the position of the ${i}^{th}$ small dung beetle at the ${t}^{th}$ iteration, with ${C}_{1}$ being a normally distributed random number, and ${C}_{2}$ being a random vector within the range of $\left(\mathrm{0, 1}\right)$.

Thieving dung beetles, a prevalent natural occurrence depicted in Figure 4, target dung balls from other beetles; this behavior is highlighted in Eq (5), where ${X}^{b}$ represents the optimal food source, implying that the vicinity of ${X}^{b}$ is prime territory for food competition, with the thief's position information being updated accordingly throughout the iteration process as shown in the Eq (7) below.

where ${x}_{i}\left(t\right)$ represents the positional data of the ${i}^{th}$ thief at the ${t}^{th}$ iteration, while $g$ is a normally distributed random vector with a size of $1\times D$. $S$ denotes a constant value.

3.   The proposed enhanced DBO

This section presents proposed enhancements to the DBO algorithm, alongside the corresponding framework in adaptive Kriging-based TRA problems.

3.1.   Improved tent chaotic mapping (initialization)

One crucial aspect of optimization is the initialization of the population, which significantly influences the algorithm's convergence and ability to find optimal solutions efficiently. Traditional methods often rely on random or uniform initialization techniques, which may not always lead to effective exploration of the search space.

In recent years, chaotic dynamics have emerged as a promising avenue for improving population initialization in optimization algorithms ^[73]. Chaos theory deals with complex, deterministic systems that exhibit seemingly random behavior. One particular chaotic mapping function, known as the tent map, has garnered attention for its ability to introduce controlled randomness into optimization processes.

The tent map operates on the principle of iteratively transforming an initial condition through a nonlinear equation, resulting in chaotic behavior characterized by sensitivity to initial conditions ^[74]. This property makes the tent map suitable for generating diverse and evenly distributed initial populations, enhancing the exploration capabilities of optimization algorithms. Using this property, chaos can be harnessed to design population initialization processes that promote a more even distribution, thereby creating ample opportunities for generating potential optimal solutions. The tent chaotic mapping is represented by the Eq (8).

Examination of the chaotic iterative sequence derived from the tent chaotic mapping function reveals the presence of minor cycles and unstable periodic points. To mitigate the emergence of such cycles and points throughout the iteration process, an element of randomness can be incorporated into the initial tent chaotic mapping formula. This enhanced formulation of the tent chaotic mapping is depicted in Eq (9).

where rand is a random variable from 0 to 1.

The improved tent chaotic mapping performance is tested with the simulation, as shown in Figure 5.

The illustration indicates that incorporating the random variable results in a more evenly spread-out set of initial values. This modification preserves the randomness, ergodicity, and regularity inherent in the original tent chaotic mapping. The improved tent chaotic mapping is implemented to replace random initialization in the DBO algorithm, thereby improving the distribution quality of the initial population within the search space and bolstering its global search effectiveness.

3.2.   Gaussian random walk (exploration)

The Gaussian function calculates the probability density value of a Gaussian distribution at a given point $x$ with mean $\mu$ and standard deviation $\sigma$. It follows the Eq (10) for a Gaussian distribution:

The function returns the probability density value $y$ at the given point $x$ within the distribution defined by $\mu$ and $\sigma$. $\sigma$ is a parameter that controls the standard deviation of the Gaussian distribution. It determines the magnitude of the perturbations applied to each solution.

The weight or probability associated with the current solution based on its similarity to the best solution is calculated by the Gaussian random walk ^[75]. The weight is computed using the Gaussian function. The step size is then calculated by multiplying the weight with the difference between the current solution and best solution. The direction and magnitude of the perturbation to be applied to the current solution are determined by this step size.

Finally, the perturbation is added to the current solution. Random numbers are generated from a standard normal distribution (rand), which are then scaled by the step size. The resulting scaled values are added element-wise to the current solution, resulting in an updated solution, as shown in Eq (11).

3.3.   The Procedure of the Proposed EDBO Algorithm

In this study, two strategies are adopted to improve the exploration and exploitation ability of the DBO algorithm. To solve optimization problems, the initialization of the population adopts an improved tent chaotic mapping method, enhancing the randomness, ergodicity, and regularity of the population. Then, the position of the thief beetle, which is explained in Eq (7), is improved using the Gaussian random walk exploration.

The main procedures of the proposed EDBO are as follows:

(1) Initialize the population of blood mass randomly and update with the tent chaotic mapping function.

(2) Calculate the fitness values and their corresponding position by Eq (1).

(3) Update the position using Eq (2) up to Eq (6).

(4) Then, update the thief position using the Gaussian random walk-by Eq (11)

(5) If the current optimal position has better fitness than the previous one, replace the previous one with the current position. Then, return to step (2) if the iteration count is not reached.

(6) If no, check the iteration count and stop the iteration if the iteration count is reached.

The flowchart of the procedures is shown in Figure 6.

4.   Case study and discussion

The standard real-world benchmark functions and three widely used TRA test cases are used to validate the performance of the proposed method in this section, such as the accuracy and efficiency. The proposed EDBO is compared with three notable metaheuristic algorithms: WOA, PSO, DE, GWO, and the original DBO. Table 1 meticulously outlines the relevant parameters for each algorithm.

CEC 2005 and CEC 2021 are widely used standard real-world benchmark test function sets. In this section, they are used to illustrate the performance of the EDBO. ${F}_{min}$represents a theoretical minimum. The maximum iterative number ${T}_{max}$ and population size are 100 and 30, respectively.

Ten benchmark functions from CEC 2005 and ten benchmark functions (F1–F10) from CEC 2021 were used; the proposed EDBO method has the lowest error among all comparative algorithms. Among the 10 test functions, the standard deviation of the proposed algorithm in calculating the 10 test functions is the smallest. This indicates that while EDBO has high accuracy, it also has high robustness. Figures 7 and 8 show the iterative convergence curves of different algorithms and the corresponding two-dimensional views of the test functions. The figures indicate that in the majority of instances, the slope of the convergence curve of the proposed method is the largest. At the same time, it can also quickly jump out of the local optimum and reach the global optimum.

In Figures 7 and 8, the proposed EDBO method is presented with a red color having a triangle mark. The overall performance of the proposed method is the best for all the benchmark functions, which illustrates the applicability of the proposed method when targeting most functions.

From Figures 7–9, it can be seen that the proposed method can reach the global optimum with the fewest number of function calls in all of the test functions. In different convergence curves, the proposed method has the fastest convergence speed in many of the test functions. The overall convergence speed is fast, which illustrates the less computational cost of the proposed method compared with other algorithms.

5.   The proposed EDBO method for TRA

The conventional TRA problem is to figure out how likely it is that a structure will stay safe within a certain amount of time, given several dynamic uncertainties. The adaptive Kriging-based TRA algorithm uses surrogate modeling to approximate the performance function. Since it is an adaptive process, it will use a certain kind of learning function that finds the best time point that will maximize the prediction error. This results in an optimization problem that maximizes the prediction error which is based on the learning function. In this section, MSE is maximized to find the next time point in the adaptive Kriging process. As a result, the approximate MPPT (AMPPT) coupled with the EDBO is proposed. Furthermore, the procedure of the EDBO-based AMPPT algorithm is outlined as follows:

(1) Discretize time into N initial samples.

(2) Apply the MPP search at the time point ${t}^{*}$.

(3) Build a Kriging model for the performance function.

(4) Using the prediction variance as a prediction error, find the next ${t}^{*}$ that maximizes it.

(5) Check for the target error and update the Kriging model.

(6) If the target error is met, calculate the time-varying reliability by spectral decomposition and the MCS method ^[76].

The flowchart of the EDBO-based TRA is shown in Figure 10.

6.   The proposed EDBO-based TRA examples

6.1.   Cantilever tube structure

The tube structure has a diameter and a wall thickness of $d$ and $h$, respectively ^[77]. The component is shown in Figure 11, and the structure is subjected to three random forces $\left({F}_{1}\left(t\right), {F}_{2}, P\right)$ as well as one random torque $T\left(t\right)$. The force ${F}_{1}\left(t\right)$ and torque $T\left(t\right)$ are represented as stochastic processes. The initial yield strength of the component is ${R}_{0}$. The structure's yield strength, $R\left(t\right) = \left(1+0.01t\right){R}_{0}$, decreases monotonically due to material deterioration. The structure will collapse if the maximal stress $\mathit{max}\left(t\right)$ is greater than the yield strength $R\left(t\right)$.

At elevated temperatures, the time-varying limit state function is given by:

where ${\sigma \left(t\right)}_{max}$ is calculated by

The reliability of this structure has been monitored over five years. The proposed AMPPT-EDBO techniques are compared with previously done MCS, TDTRA, and the adaptive Kriging-based AMPPT using differential evolution (AMPPT-DE).

For MCS, 300 equally spaced time points are discretized within the interval [0, 5] years, generating 106 samples. Two cases with different discrete step lengths have been presented for TDTRA, AMPPT-DE, and AMPPT-EDBO. Case 1 is set to 1/2 year, and case 2 is set to 1/6 year, respectively. In the AMPPT-DE and AMPPT-EDBO, three-time instants (0, 2.5, and 5) are employed to generate the initial MPP samples. Table 2 provides the probability characteristics of each variable that affects the cantilever tube structure. It's worth noting that some entries are marked as "NA" (Not Available), indicating that specific data points were not applicable in those particular cases.

In Figure 12, the MCS method is represented by the yellow line, while the AMPPT-DE method is presented by deep orange and light orange colors for cases 1 and 2, respectively. The TDTRA method is shown by purple and green colors for cases 1 and 2, respectively. Finally, the proposed AMPPT-EDBO method is presented with blue and red colors for cases 1 and 2, respectively. The probability of failure increases gradually for all the methods. In addition, both in cases 1 and 2, the three methods have nearly comparable failure probability with variation of ±0.0001. The details of the comparison are given in Table 3.

Tables 3 and 4 present the failure probability and reliability index, respectively. As inferred from Tables 3 and 4, the failure probability of the TDTRA, AMPPT-DE, and AMPPT-EDBO increases with the number of discrete step lengths, represented as cases 1 and 2.

Table 5 presents the accuracy and efficiency of the methods. Since MCS generates huge samples, its prediction is considered to be the more reliable and accurate. Therefore, the accuracy of other methods is calculated by comparing the failure probability with that of MCS. The efficiency is the number of calls that are required to calculate the failure probability. Table 5 shows that the accuracy of all three discretization-based methods is similar with slight deviation. In case 2, AMPPT-EDBO has the accuracy of $2.97\%$, AMPPT-DE has $2.97\%$, while TDTRA's is $3.9\%$. However, the efficiency comparison between the methods shows a significant difference. MCS has by far the largest number of calls, which is 3×10⁸. TDTRA has a relatively lower number of calls: 1541 for case 1 and 3514 for case 2; however, it is still computationally costly. AMPPT-DE needs 356 calls for case 1 and 671 for case 2. The proposed AMPPT-EDBO needs the lowest number of calls, which is 263 for case 1 and 263 for case 2.

Therefore, the proposed AMPPT-EDBO method decreases computational cost without compromising the accuracy of the reliability analysis. Since the adaptive Kriging simplifies the limit state function, the high computational cost needed for the highly nonlinear limit state function is significantly decreased. In addition, the enhanced DBO further improves the adaptive Kriging by quickly finding the best time point that minimizes fitness error in the adaptive Kriging process.

6.2.   Corroded steel beam

This example is based on the simply supported beam, as shown in Figure 13. The steel beam has an $L$ span and an $A-A$ rectangular cross-section ^[78]. The corrosion effect causes the initial width and height of the cross-section, ${w}_{0}$and ${h}_{0}$, to linearly decrease with time. Both a load brought on by the self-weight and a concentrated dynamic load $F\left(t\right)$ are imposed on the beam. Failure will occur if the maximum stress on the beam exceeds the yield strength of the material.

The limit state function is described as follows:

In this context, the symbols ${\rho }_{st}$ and ${\sigma }_{y}$ represent the yield strength and density of the steel material, respectively, while $k$ indicates the corrosion rate. For this particular situation, the variables $y, {w}_{0}$, and ${h}_{0}$ are treated as stochastic variables. The concentrated load $F\left(t\right)$ is assumed to follow a stationary Gaussian process. The probabilistic properties of all these variables are outlined in Table 6, provided below.

The reliability of this structure has been monitored over 30-years period. The proposed AMPPT-EDBO techniques are compared with previously done MCS, TDTRA, and AMPPT-DE.

For MCS, 500 equally spaced time points are discretized within the interval [0, 30] years, generating 2×10⁶ samples. Two cases with different discrete step lengths have been presented for TDTRA, AMPPT-DE, and AMPPT-EDBO. Case 1 is set to 1/3th year, and case 2 is set to 1/10th year, respectively. In the AMPPT-DE and AMPPT-EDBO, three-time instants (0, 15, and 30) are employed to generate the initial MPP samples.

In Figure 14, the MCS method is represented by the yellow line, while the AMPPT-DE method is presented by deep orange and light orange colors for cases 1 and 2, respectively. For cases 1 and 2, the TDTRA method is represented by the colors purple and green, respectively. In cases 1 and 2, the proposed AMPPT-EDBO method is shown in blue and red, respectively. The probability of failure rises steadily for every method with an increase in time. Furthermore, it appears that the three methods in Cases 1 and 2 have comparable failure probabilities, which affirms that the accuracy of the prediction methods is fairly comparable. Table 7 provides the comparison's details.

Tables 7 and 8 present the failure probability and reliability index, respectively. As inferred from Tables 7 and 8, the probability of failures using the TDTRA, AMPPT-DE, and AMPPT-EDBO increases with an increase in discretization interval, which is proved in cases 1 and 2. Both failure probability and reliability index are less for smaller discretization intervals, as in case 1, and gradually increase as the discretization interval increases, evidenced in case 2.

Table 9 displays the results for the accuracy and efficiency of the four methods. According to Table 9, TDTRA, AMPPT-DE, and the proposed AMPPT-EDBO methods have comparable accuracy. In case 2, AMPPT-EDBO achieves an accuracy of $0.37\%$, while AMPPT-DE has $14.13\%$, and TDTRA's accuracy is $12.3\%$. However, when comparing efficiency, there is a significant contrast between the methods. MCS requires the highest number of calls, which is 1×10⁹. TDTRA has a relatively lower number of calls, 1817 for case 1 and 4885 for case 2. As shown in Table 9, the number of calls increases with a larger discretization interval, represented by cases 1 and 2. AMPPT-DE necessitates 335 number of calls for case 1 and 335 for case 2. In comparison, the proposed AMPPT-EDBO method requires the fewest number of calls, 304 for case 1 and 299 for case 2. Therefore, the proposed AMPPT-EDBO method significantly improves efficiency compared to the rest of the presented methods.

6.3.   Solid rocket engine shell

The shell structure of the solid rocket engine's time-variant reliability is examined in this illustration ^[79]. Due to their simple architecture and low number of moving parts, solid rocket engines are widely used to boost satellite launch vehicles. We give a schematic of a solid rocket engine, where we suppose that ${D}_{t}$ stands for the nozzle throat's diameter and ${D}_{s}$ stands for the inner shell's diameter. The propellant has a length of $L$ and a starting thickness of ${H}_{0}$. The solid propellant is housed inside the cylindrical metal casing. A hole in the cylindrical shell is used as the combustion chamber. A flame front develops on the surface of the propellant once it has caught fire. The solid rocket engine is presented in Figure 15.

The system of differential equations can be computed by the Runge-Kutta method, which will yield the time-variant pressure $P\left(t\right)$. We use the parameters $\left({C}^{*}, \rho, {D}_{t}\right)$ as random variables for calculating $P\left(t\right)$. As shown in Figure 15, solid rocket engines are used on air-to-air and air-to-ground missiles, on model rockets, and as boosters for satellite launchers. In a solid rocket, the fuel and oxidizer are mixed into a solid propellant packed into a solid cylinder. A hole through the cylinder serves as a combustion chamber.

The system of differential equations can be computed by:

where ρ is the density of the propellant. ${C}^{*}$ is a coefficient of the chamber characteristic velocity with uncertainty. a, n, and Γ are uniform values given by $a = 7.9752\times 1{0}^{-6}$, $n = 0.4356$, $\varGamma = 0.6528$.

The ${A}_{t}$ and ${A}_{f}\left(t\right)$ are the areas of the nozzle throat and the flame front, respectively, which can be computed by:

The beginning circumstances for the system of differential equations are given by

To avoid excessively computing the differential equations, a polynomial response surface is built for $P\left(t\right)$ to formulate it as a function of the time parameter and random variables$\left({C}^{*}, \rho, {D}_{t}\right)$. This is represented by:

Therefore, the time-dependent limit state function of the solid rocket engine shell is given by:

The time period analysis is in the interval of ^[7,10] seconds. For the MCS, 300 discretized equal distant time instants are created, and 10⁶ samples are generated at each instant. Similar to the previous examples, two cases were set, with discrete step lengths of 1/3rd and 1/10th of a second for cases 1 and 2, respectively. Once again, the proposed AMPPT-EDBO method is compared with previously done MCS, TDTRA, and AMPPT-DE. The MCS is set similarly to the previous example 6.2.

Table 10 lists the random variables for the components with their probabilistic features, including the type of distribution, mean value, coefficient of variation, and autocorrelation coefficient function, if available.

In Figure 16, the MCS method is represented by the yellow line, while the AMPPT-DE method is presented by deep orange and light orange colors for cases 1 and 2, respectively. For cases 1 and 2, the TDTRA method is represented by the colors purple and green, respectively. In cases 1 and 2, the proposed AMPPT-EDBO method is shown in blue and red, respectively. In a similar trend to the previous examples 5.1 and 5.2, the probability of failure increases with an increase in time for every method. Table 11 provides a detailed comparison of the probability of failure for each method.

Tables 11 and 12 present the failure probability and reliability index, respectively. With a similar trend with previous examples 6.1 and 6.2, the accuracy of the TDTRA, AMPPT-DE, and AMPPT-EDBO increase with an increase in discretization interval, which is represented in cases 1 and 2.

Table 13 presents the accuracy and efficiency of the methods. In case 2, the relative error of AMPPT-EDBO is 2.43%, while AMPPT-DE's is 14.63% and TDTRA's is 14.63%. In this example, the accuracy of the proposed AMPPT-EDBO method is much better than the TDTRA and AMPPT-DE methods. In addition, the efficiency comparison between the methods shows a significant difference. MCS has by far the largest number of calls, which is 3×10⁸. TDTRA has a relatively lower number of calls: 1490 for case 1 and 4854 for case 2. AMPPT-DE needs $664$ calls for case 1 and 557 for case 2. The proposed AMPPT-EDBO needs the lowest number of calls, which is 435 for case 1 and 456 for case 2.

The proposed AMPPT-EDBO method not only has efficiency improvement by decreasing computational cost but is also proven to have better accuracy compared to traditional TDTRA and Kriging-based AMPPT-DE.

6.4.   Bar truss

This example considers the 23-bar truss ^[61] presented in Figure 17.

The cross-section areas and Young's modulus of the bars are considered random variables, denoted as ${A}_{1}$ and ${E}_{1}$ for horizontal bars, and ${A}_{2}$and ${E}_{2}$ for sloping bars, respectively. The truss is subjected to six time-variant stochastic loads ${F}_{i}\left(t\right)\left(i = \mathrm{1, 2}, ..., 6\right)$. The maximum defection of the truss, denoted as ${d}_{m}\left(t\right)$ in Figure 17, is required to be less than the threshold $d\left(t\right) = {d}_{0}\left(1-0.01t\right),$ where ${d}_{0}$ is the initial value of $d\left(t\right)$. Therefore, the time-variant performance function is given by

where ${d}_{m}\left(t\right)$ is obtained by the finite element method.

Eleven stochastic process/random variables are involved in this example, and the probabilistic characteristics of all involved variables are listed in Table 14. The 23-bar truss is assessed over a time span of [0, 10] years. The proposed AMPPT-EDBO techniques are compared with MCS, TDTRA, and AMPPT-DE.

The simulation steps for the 23-bar truss are presented as follows. Since the force is a Gaussian process, 10 sample points are generated by using the mean, standard deviation, and autocorrelation values, which are already given. Then, since the area of the truss bar is also a normal distribution, it is set as an input variable in the Ansys software, with the directional deformation at the desired point set as an output variable. The input data points are then generated according to the given normal distribution information and set as input in the Ansys software. MATLAB is used to call the Ansys software to solve the directional deformation for the given input points, and the output is saved in an Excel file. The distribution of the directional deformation is then estimated from the recorded Excel file, and the failure probability of the truss bar structure is solved using the given limit state function.

In Figure 18, the ANSYS simulation result of the 23-bar truss is presented for total deformation in (a) and direction deformation along the y-axis in (b). The total deformation plot in part (a) provides an overall view of how the truss structure is deforming under the applied loads or boundary conditions. The y-direction deformation plot in part (b) isolates the deformation in a specific axis, which can be useful for understanding the directionality of the structural response and identifying potential problem areas or stress concentrations along the y-axis.

In Figure 19, the MCS method is represented by the yellow line, while the AMPPT-DE method is presented by deep orange and light orange colors for cases 1 and 2, respectively. For cases 1 and 2, the TDTRA method is represented by the colors purple and green, respectively. In cases 1 and 2, the proposed AMPPT-EDBO method is shown in blue and red, respectively. As in the examples, the probability of failure increases gradually with time for all the methods. The details of the comparison are given in Table 15.

Tables 15 and 16 present the probability of failure and reliability index, respectively. As is already the trend in previous examples, the accuracy of TDTRA, AMPPT-DE, and AMPPT-EDBO increases with an increase in the discretization interval, which is seen in cases 1 and 2.

Table 17 presents the accuracy and efficiency of the methods. In case 2, AMPPT-EDBO has an accuracy of 21.2%, AMPPT-DE has 22.6%, and TDTRA's is 22.4%. In this example, the accuracy of the proposed AMPPT-EDBO method is much better than the TDTRA and AMPPT-DE methods. In addition, the efficiency comparison between the methods shows a significant difference. MCS has by far the largest number of calls, which is 5×10⁸. TDTRA has a relatively lower number calls: 315 for case 1 and 630 for case 2. AMPPT-DE needs 84 calls for case 1 and 84 for case 2. The proposed AMPPT-EDBO needs the lowest number of calls, which is 57 for both cases 1 and 2.

The proposed AMPPT-EDBO method not only has efficiency improvement by decreasing computational cost but is also proven to have better accuracy compared to traditional TDTRA and Kriging-based AMPPT-DE.

7.   Conclusions

In this study, a novel EDBO algorithm is proposed to solve the optimization problem in Kriging model-assisted adaptive TRA. The algorithm combines tent chaos and Gaussian random walk to improve the performance of optimization techniques. The proposed optimization method is applied to solve the optimization problem in the process of adaptive Kriging model-assisted reliability analysis. Finally, the effectiveness of this method is illustrated by four examples of mechanical structure models. Compared with the traditional method, this method reduces the computational cost while maintaining a low relative error rate, which reflects the high efficiency, accuracy, and reliability of the proposed method. In the future, it is a direction with great potential to develop better optimization algorithms to solve complex reliability analysis problems. In addition, we will also develop a new reliability analysis framework based on excellent optimization algorithms to evaluate the reliability of more complex engineering problems in future work.

Author contributions

Yunhan Ling: Conceptualization, methodology, software, validation, formal analysis, investigation, writing-original draft, writing-review & editing, funding acquisition; Yiqing Shi: Formal analysis, investigation, software, validation, funding acquisition, writing-review & editing; Huimin Hou: Formal analysis, validation, investigation, writing-review & editing; Lidong Pan: Methodology, software, validation, writing-review & editing; Hao Chen: Software, validation, formal analysis, investigation, writing-review & editing; Peixin Liang: Conceptualization, methodology, software; Shiyuan Yang: Conceptualization, methodology, writing-original draft, writing-review & editing, funding acquisition; Peng Nie: Formal analysis, investigation, writing-review & editing; Jiahao Han: Formal analysis; Debiao Meng: Conceptualization, methodology, writing-review & editing, funding acquisition. All authors have read and agreed to the published version of the manuscript.

Acknowledgments

This research was funded by the National Key Research and Development Program (Grant No. 2022YFB3706904); the Guangdong Basic and Applied Basic Research Foundation (No. 2022A1515240011); China Scholarship Council (No. 202406070025 and 202406070043).

Conflict of interest

The authors declare there is no conflict of interest.