A Grasshopper Optimization Algorithm-Based Response Surface Method for Non-Probabilistic Structural Reliability Analysis with an Implicit Performance Function

: Non-probabilistic reliability analysis has great developmental potential in the ﬁeld of structural reliability analysis, as it is often difﬁcult to obtain enough samples to construct an accurate probability distribution function of random variables based on probabilistic theory. In practical engineering cases, the performance function (PF) is commonly implicit. Monte Carlo simulation (MCS) is commonly used for structural reliability analysis with implicit PFs. However, MCS requires the calculation of thousands of PF values. Such calculation could be time-consuming when the structural systems are complicated, and numerical analysis procedures such as the ﬁnite element method have to be adopted to obtain the PF values. To address this issue, this paper presents a grasshopper optimization algorithm-based response surface method (RSM). First, the method employs a quadratic polynomial to approximate the implicit PF with a small set of the actual values of the implicit PF. Second, the grasshopper optimization algorithm (GOA) is used to search for the global optimal solution of the scaling factor of the convex set since the problem of solving the reliability index is transformed into an unconstrained optimal problem. During the search process in the GOA, a dynamic response surface updating strategy is used to improve the approximate accuracy near the current optimal point to improve the computing efﬁciency. Two mathematical examples and two engineering structure examples that use the proposed method are given to verify its feasibility. The results compare favorably with those of MCS. The proposed method can be non-invasively combined with ﬁnite element analysis software to solve non-probabilistic reliability analysis problems of structures with implicit PF with high efﬁciency and high accuracy.


Introduction
To address the uncertain factors in a structure, it is necessary to reasonably select a mathematical model that quantifies the uncertain factors and structural reliability. Many analytical-and simulation-based methods have been developed to calculate the failure probability and reliability index of structures [1,2]. Among them, Monte Carlo simulation (MCS) is probably the most widely used high-precision structural reliability analysis method. However, when the probability of structural failure is small, MCS often needs to simulate a large number of samples to obtain accurate probability values, which will increase a lot of computational costs. Some improved methods, such as the subset simulation method [3], the line sampling method [4], and the importance sampling method (IS) [5], still have a high computation when facing complex problems. Analytical methods, such as first-order 2 of 21 reliability methods (FORMs) and second-order reliability methods (SORMs), are computationally efficient but are often not stable enough in the face of structures with implicit PFs. For some improved algorithms, such as the directional stability transformation method (DSTM) [6] and finite-step adaptive length (FAL) [7], it is difficult to solve the convex set reliability problem [8]. On the other hand, by the type of random variable in the structures, Cremona [9] divided structural reliability models into random reliability models, fuzzy reliability models, non-probabilistic reliability models, and mixtures of the three models according to the mechanism of uncertainty factors. Among the various uncertainty models, the most widely used is the probabilistic reliability model based on probability theory and stochastic processes [10][11][12]. However, it is difficult to obtain sufficient statistical sample data in practical engineering structures to construct the probability distribution function of a probabilistic model. Additionally, a small fitting error of the probability distribution function of a random reliability model can lead to a large structural reliability calculation error [13]. This is especially true when there are few samples with which the model can calculate [14].
To address the condition of insufficient samples, in 1993, Ben-Haim [15] and Elishakoff [16] proposed a non-probabilistic reliability model for structural analysis and applied the maximum uncertainty allowed in the system to evaluate the structural reliability. Nonprobabilistic reliability theory uses an artificially defined set model to describe the uncertain variables of the structure. The non-probabilistic reliability model only needs the range of the uncertain variables and then constructs the mapping relationship between the uncertain variables and the corresponding feasible domain of the structure under the combined action of multiple variables. It can be seen from this that it is particularly important to use the non-probabilistic reliability method in the reliability analysis of complex engineering structures. Due to factors such as safety and cost, it is difficult to carry out enough field experiments to obtain enough actual samples to establish reliable probability distribution functions in complex engineering structures [17]. Furthermore, small deviations in the probability density function due to data noise are likely to lead to large result errors [13]. However, in the same situation, using a small number of field experiments, the range of random variables can be quickly obtained, and then a structural reliability analysis with high reliability can be constructed through the non-probabilistic reliability method. Then, by comparing the set model and the limit state surface defined by the performance function (PF), the structural reliability index can be evaluated. Subsequently, scholars have conducted many studies that can be summarized into two categories: proposing different non-probabilistic reliability models aimed at better describing multivariable problems and studying how to quantify non-probabilistic reliability indicators according to different models. In 1998, Qiu and Elishakoff [18] used interval theory to establish an interval model to optimize a six-bar truss. In 2011, Wang et al. [19,20] proposed a non-probabilistic reliability model measurement index based on an interval model, where the index was defined by the ratio of the volume of the safety domain of the structure to the total volume of the basic interval variable domain. The interval model can only deal with independent variables [21]; therefore, scholars further proposed ellipsoid models to deal with dependent variables. In 1998, Pantelides et al. [22] proposed an ellipsoid model and compared the structural analysis results of different working conditions based on the ellipsoid model. In 2011, Jiang et al. [23] proposed a modeling method to establish a multidimensional ellipsoid model. To make the dependence of uncertain variables explicit, they introduced a covariance matrix to evaluate the correlations of dependent variables. In 2016, Kang and Zhang [24] constructed a mathematical formula system to develop a compact ellipsoid model. Recently, to address the problem of including independent and dependent variables, scholars have proposed more general models. In 2015, Jiang et al. [21] proposed a multidimensional parallelepiped model that can efficiently deal with complex 'multi-source uncertainty' problems. Guo et al. [25] proposed a mixed model composed of a hypercube interval model and a hyper-ellipsoid model. In 2021, Cao et al. [26] proposed a quantifying strategy for non-probabilistic uncertainties by combining the general polygonal convex set model and cluster polygonal convex set model. However, with the development of non-probabilistic reliability models and quantitative indices of reliability, solving nonprobabilistic reliability problems quickly and accurately has become a restriction on related field development.
To address the above problems, scholars have explored many methods. Jiang et al. [23] applied the first-order reliability method (FORM) to a non-probabilistic convex model based on multidimensional ellipsoids for uncertainty. Wang et al. [27] proposed an improved second-order reliability method (SORM) for reliability analysis in the presence of both random and interval variables. These methods perform well in reliability analysis with explicit PFs; however, in practical engineering cases, the PF is commonly implicit. Therefore, it is difficult to obtain the exact boundaries of the structural failure domain and the reliable domain. For more complex structures, the reliability can only be calculated by numerical simulation methods such as finite element analysis (FEA), which are very time-consuming in most practical engineering applications. Therefore, the use of an explicit surrogate model of PFs for structural reliability analysis has become mainstream. One of the most commonly used methods is the response surface method (RSM). In 1989, Faravelli [28] first proposed a quadratic polynomial regression model to approximate a function, and this method is also called the classical RSM. The response surface can quickly and accurately fit an implicit PF and, at the same time, can eliminate the interference caused by sample noise through probabilistic methods. RSMs can be simply divided into two types: regression-based models and interpolation-based models [29]. Regression-based models, such as Gaussian process regression [12,30,31], use the probabilistic properties of known samples to perform regression analysis to fit functions in the feasible domain. An alternative model is based on general polynomial chaos expansion (gPCE) [32][33][34]. In this surrogate model, the random variables are placed in the random space expanded by orthogonal polynomials. The coefficients of each dimension are multiplied by the corresponding orthogonal polynomials, and then added to obtain the explicit surrogate model of the implicit PF. Interpolation-based models, such as support vector machines [35,36] and artificial neural networks [37][38][39], correct the values of the weight parameters of the model by learning the characteristics of the known samples to interpolate the function in the feasible domain. However, the RSMs mentioned above are all aimed at probabilistic reliability, and research on RSMs aimed at non-probabilistic reliability is still very rare. This paper proposes a grasshopper optimization algorithm (GOA)-based RSM for the non-probabilistic reliability analysis of structures with implicit PFs. In this method, a quadratic polynomial response surface is presented to fit the PF. Notably, for general engineering structures, although the newly proposed RSM has high accuracy, the fitting speed is slow. Additionally, the hyper-parameters in models such as support vector machine (SVM) are difficult to determine without manual parameter adjustment [40]. In contrast, the application of the quadratic polynomial RSM in the field of probabilistic reliability is very mature [41] and can be quickly fitted while largely guaranteeing the fitting accuracy in general structures. Then, the proposed method uses a dynamic fitting strategy to find the approximate location of the generalized design point and then performs local fitting of the PF around the generalized design point. To find the globally optimal solution faster in the process of calculating the non-probabilistic reliability index, the GOA is used to search for the reliability index and optimal design point of the structure. This paper is organized as follows: Section 2 introduces the convex model theory to quantify the nonprobabilistic reliability with the reliability index. Section 3 presents the framework of the proposed method. According to the framework, Section 4 introduces the key technologies of the proposed method. Section 5 introduces the detailed implementation of the proposed method based on the key technologies. Section 6 gives four examples-two mathematical examples and two engineering structure examples-to verify the feasibility of the proposed method. Finally, Section 7 summarizes the conclusion for the proposed method according to the examples.

Reliability Index of Non-Probabilistic Reliability
A non-probabilistic reliability model is an effective way to deal with structures with insufficient samples. The question arises of how to characterize the safety of structures in a non-probabilistic setting. Many scholars have proposed multiple reliability indices to quantify the non-probabilistic reliability [42][43][44]. The use of a scaling factor for convex sets is one of the most common methods.
According to the properties of the convex model, the total uncertainty of a structure is the intersection set of m convex sets of all uncertainties [37]. This intersection set can also be considered a feasible domain for the structure to remain stable under the impact of multiple uncertain variables. It can be represented as: where α = {α 1 , α 2 , . . . , α n } is the size parameter of the convex set family of total uncertain variables, X = {x 1 , x 2 , . . . , x n } is the position parameter of the convex set family of total uncertain variables, and X m = (x n m−1 +1 , . . . , x n m ) represents a subset of uncertain variables with the same convex set. The stability of the structure under a certain working condition can be determined by Equation (1). The scale factor λ can be used to freely scale the intersection set in Equation (1). The scaled intersection set is as follows: The factor acts on the size parameter of each convex set. When the above intersection set in Equation (2) is increased by the factor λ to intersect the failure domain divided by the PF, this indicates that the structure is in a failure state under some combination of uncertain variables. Therefore, λ must be reduced until the intersection is tangent to the failure domain. This tangent point represents the generalized design point of the structure on the PF [45]. Figure 1 shows the non-probabilistic reliability index of a PF considering two uncertain variables in x-space. The failure domain is divided by PF g(X) = 0, and the blue solid area denotes the normalized convex model. It is defined by the intersection of an ellipsoid convex set X 1 and an interval convex set X 2 . When the convex set expands to be tangent to the PF, the value of λ is the non-probabilistic reliability index η, and the tangent point A is the generalized design point. Therefore, the minimum λ value can be defined as the reliability index η of structural reliability as follows: The index η can be used to simply measure the degree of non-probabilistic reliability to judge whether the structure is in a safe or failure state. The structural reliability problem is then transformed into the problem of finding the optimal solution of η.
The index η can be used to simply measure the degree of non-probabilistic re to judge whether the structure is in a safe or failure state. The structural reliability p is then transformed into the problem of finding the optimal solution of η.
The scaling factor of the convex set is used to obtain the optimization formu The optimization formula is obtained by transforming Equation (3) into an uncons optimization problem through the penalty method: ρ x x is the functional relation of the convex set model and P is the The scaling factor of the convex set is used to obtain the optimization formula of η. The optimization formula is obtained by transforming Equation (3) into an unconstrained optimization problem through the penalty method: where ρ(x, x) is the functional relation of the convex set model and P is the penalty coefficient, which is generally a large constant. Finally, the non-probabilistic reliability problem is transformed into the problem of finding the optimal solution of Equation (4), which is based on the convex set model X and PF g(X). In other words, Equation (4) can be considered the fitness function of some optimization algorithm, and the optimal solution can be sought.

Framework of the Proposed Method
To solve non-probabilistic reliability problems quickly and accurately, there are two possible ways: (1) while ensuring the accuracy of calculation, the time of performing FEA should be reduced as much as possible; that is, the number of calls to the structural PF should be reduced; (2) while ensuring accuracy, the number of samples required by the surrogate model should be reduced as much as possible; that is, the fitting time of the surrogate model should be reduced. As it is also difficult to obtain sufficient samples in practical engineering cases, a dynamic updating mechanism can be used to correct the fitting errors of the surrogate model in local areas on the basis of small samples. The framework of the dynamic surrogate model for non-probabilistic reliability analysis is shown in Figure 2. x* and η are the optimal design point and reliability index searched by the optimization process. ness function. The third part uses the searched optimal design point as the new ce point to construct new sample points and refit the PF if the output does not reach con gence. Finally, the iteration reaches convergence and outputs the optimal solution. T three parts make the computational cost of structural reliability very low. The key t nologies of each part of the proposed method are presented in Section 4.

Flexible Construction of the Response Surface with a Small Dataset
To construct the response surface function, a sample dataset related to the uncer variables must be obtained. The initial sample datasets are constituted by a series of po xi (i = 1, 2, …, n) and a series of corresponding PF values zi calculated by FEA. The qu of the sample points directly affects whether the quadratic polynomial response sur can approximate the real PF well.
In traditional structure reliability analysis, the methods of sample point design c monly used in response surface construction are the random sampling design, two-l factorial design [46,47], central composite factorial design [28], Bucher design [48], L hypercube sampling design [49], IS design [5,50], and directional simulation [51]. sampling design method this paper uses is the Bucher design. This is because the num of sample points constructed by the Bucher design grows more slowly than that of o design methods when the number of undetermined variables increases. When the st ture is complex and the uncertain variables have a complex high-dimensional form increasing number of sample points will reduce the efficiency of fitting the respo To establish the framework, the dynamic surrogate model needs to include the following three parts: (1) the response surface model; (2) the optimization process; and (3) the dynamic updating process. The first part constructs the response surface model to fit the PF and constructs the fitness function in Equation (4) for non-probabilistic reliability analysis. The second part uses optimization methods to find the optimal value of the fitness function. The third part uses the searched optimal design point as the new center point to construct new sample points and refit the PF if the output does not reach convergence. Finally, the iteration reaches convergence and outputs the optimal solution. These three parts make the computational cost of structural reliability very low. The key technologies of each part of the proposed method are presented in Section 4.

Flexible Construction of the Response Surface with a Small Dataset
To construct the response surface function, a sample dataset related to the uncertain variables must be obtained. The initial sample datasets are constituted by a series of points x i (i = 1, 2, . . . , n) and a series of corresponding PF values z i calculated by FEA. The quality of the sample points directly affects whether the quadratic polynomial response surface can approximate the real PF well.
In traditional structure reliability analysis, the methods of sample point design commonly used in response surface construction are the random sampling design, two-level factorial design [46,47], central composite factorial design [28], Bucher design [48], Latin hypercube sampling design [49], IS design [5,50], and directional simulation [51]. The sampling design method this paper uses is the Bucher design. This is because the number of sample points constructed by the Bucher design grows more slowly than that of other design methods when the number of undetermined variables increases. When the structure is complex and the uncertain variables have a complex high-dimensional form, an increasing number of sample points will reduce the efficiency of fitting the response surface. For example, when there are n factors, the two-level factorial design will construct 2 n sample points, and the central composite factorial design will construct 2 n + 2n + 1 sample points. Both grow exponentially in number. In contrast, the Bucher design uses the difference set between the central composite factorial design and two-level factorial design, constructing only 2n + 1 sample points, which is linearly increasing. The sampling method in the two-dimensional condition is shown in Figure 3: surface. For example, when there are n factors, the two-level factorial design will construct 2 n sample points, and the central composite factorial design will construct 2 n + 2n + 1 sample points. Both grow exponentially in number. In contrast, the Bucher design uses the difference set between the central composite factorial design and two-level factorial design, constructing only 2n + 1 sample points, which is linearly increasing. The sampling method in the two-dimensional condition is shown in Figure 3: In Figure 3, the center of the initial sample points is the position parameter x* = (x1, x2, …, xn). The vertical and horizontal offset distances are kf times the size parameter αi After constructing the initial sample points, the FEM is used to calculate the PF values of the sample points. The obtained function values g(x1, x2, …, xn) and g(x1, x2, …, xi ± kfαi, …, xn) and the sample points will be combined by a one-to-one correspondence into 2n + 1 samples of the response surface model.
In this paper, a polynomial response surface is used to fit the PF. When a dynamic fitting strategy is employed, the quadratic polynomial response surface does not need to represent the response of a nonlinear model to multiple uncertainty parameters over the entire feasible region. Instead, it is only necessary to respond to the nonlinear model of the uncertainty parameters near the generalized design point, that is, to realize the local fitting of the response surface near the generalized design point. In local fitting, since the local PF itself is not too complicated, the use of high-order polynomials will not greatly improve the accuracy but will increase a lot of computational costs. So, the use of quadratic polynomials to fit nonlinear PFs is sufficient to achieve high-precision fitting. Simple quadratic polynomials are usually used, as follows: where X are the basic random variables of the structure; a, bi, ci, and dij are the undetermined coefficients of the polynomial; ( ) g X is the real PF; and ˆ( ) g X is the quadratic polynomial response surface function.
Given the variables X, the corresponding value ˆ( ) g X is calculated by the FEM.
Then, the undetermined coefficients can be solved with the following equation: where     In Figure 3, the center of the initial sample points is the position parameter x* = (x 1 , x 2 , . . . , x n ). The vertical and horizontal offset distances are kf times the size parameter α i (interval model k = 1, ellipsoid model k = e i ).
After constructing the initial sample points, the FEM is used to calculate the PF values of the sample points. The obtained function values g(x 1 , x 2 , . . . , x n ) and g(x 1 , x 2 , . . . , x i ± kfα i , . . . , x n ) and the sample points will be combined by a one-to-one correspondence into 2n + 1 samples of the response surface model.
In this paper, a polynomial response surface is used to fit the PF. When a dynamic fitting strategy is employed, the quadratic polynomial response surface does not need to represent the response of a nonlinear model to multiple uncertainty parameters over the entire feasible region. Instead, it is only necessary to respond to the nonlinear model of the uncertainty parameters near the generalized design point, that is, to realize the local fitting of the response surface near the generalized design point. In local fitting, since the local PF itself is not too complicated, the use of high-order polynomials will not greatly improve the accuracy but will increase a lot of computational costs. So, the use of quadratic polynomials to fit nonlinear PFs is sufficient to achieve high-precision fitting. Simple quadratic polynomials are usually used, as follows: where X are the basic random variables of the structure; a, b i , c i , and d ij are the undetermined coefficients of the polynomial; g(X) is the real PF; andĝ(X) is the quadratic polynomial response surface function. Given the variables X, the corresponding valueĝ(X) is calculated by the FEM. Then, the undetermined coefficients can be solved with the following equation: where λ =(A, b 1 , b 2 , . . . , b n , c 1 , c 2 , . . . , c n , d 1 , . . . , d n ) T ,ĝ(X) = (ĝ 1 ,ĝ 2 , . . . ,ĝ 2n+1 ) T , and A is the coefficient matrix.
To make the response surface fit the real PF better near the checking point, a quadratic polynomial without a cross-product term is usually selected, as follows: Note that in the non-probabilistic reliability model, if X is a non-probabilistic random variable, the quadratic polynomial response surface function is selected by the above formula. In the case of two random variables X and Y, the response surface model expressed in Equation (7) does not consider the cross product term of random variables; that is, it is equivalent to two variables being independent variables. Therefore, the quadratic polynomial response surface in this case can be obtained by adding the quadratic polynomial constructed by X and the quadratic polynomial constructed by Y, and can be expressed as the following Equation (8): The explicit quadratic polynomial response surface model established above is obtained by fitting the implicit PF with sample points. Therefore, the response surface model can be used to replace the PF to divide the structural reliability domain and failure domain and replace the PF to participate in the non-probabilistic reliability analysis of the structure, thus avoiding the impediment of the implicit PF to the structural reliability analysis. Based on the response surface function Equation (8), the optimization problem Equation (4) can be established to determine the non-probabilistic reliability.

Quickly Searching for the Global Optimal Solution Using the GOA
Due to the diversity of random variables in the structure, the PFs used for structural reliability analysis in engineering structures are often implicit multimodal functions, which means Equation (4) has multiple local optimal solutions. Traditional optimization methods are based on gradient descent to find the global optimum, which often leads to local optimization [52]. Nature-inspired algorithms, such as the genetic algorithm (GA) [53], particle swarm algorithm (PSO) [54], ant colony algorithm (AC) [55], etc., create random input values as candidate optimization schemes. Then, iterative search for the global optimum is achieved by evaluating each scheme, observing the target value, and changing and combining schemes according to the output. However, these algorithms still experience many difficulties in solving complex optimization problems [56]. For example, PSO lacks dynamic adjustment of the search speed and falls into local optimization. GA has a strong global optimization ability but weak local search ability, which makes it difficult to converge quickly. The grasshopper algorithm (GOA) [57] has an adaptively adjusted search speed and can be well combined with non-probabilistic structural reliability solution methods to speed up the solution process [58]. Therefore, GOA becomes the optimization method used in this paper to find the optimal solution of the fitness function Equation (4).
By observing and simulating the collective behavior of grasshopper swarms as they search for food, in 2017, Saremi et al. [57] proposed the grasshopper optimization algorithm (GOA). The basic idea of the GOA is as follows.
The current position of the grasshopper operator X i can be described as: where S i is social interaction, G i is the force of gravity on the i-th grasshopper, and A i is the wind advection. These parameters can be defined as follows: where d ij is the distance between the i-th and j-th grasshoppers, s is the function shown in Equation (13) to define the strength of social force, andd ij is a unit vector of the i-th and j-th grasshoppers. f is the intensity of attraction, and l is the attractive length scale. g is the gravitational constant, andê g is a unit vector towards the center of the Earth. u is a constant drift, andê w is a unit vector in the direction of the wind. Then, the next position of the grasshopper swarm can be defined as follows: wherex k i indicates the next position of the grasshopper operator x i in dimension k. ul k and f l k represent the upper and lower bounds of the model in dimension k. If the new positions of some grasshopper operators exceed the upper and lower bounds of the model space, they will be pulled back to the model space. d k ij represents the distance between grasshopper operator x k i and grasshopper operator x k j in dimension k. T k g represents the component of the currently found optimal solution in dimension k. c 1 and c 2 are adaptive shrinkage parameters, which can be calculated as: where t is the current iteration time. c max and c min are the maximum and minimum values of the adaptive shrinkage parameters. In this paper, they are set as 1 and 1.0 × 10 −5 . Through the GOA, the proposed method can efficiently search for the optimal solution to Equation (4). The proposed method will gradually approach the generalized design point and the reliability index of non-probabilistic structural reliability.

Dynamically Updating the Response Surface
After obtaining the initial quadratic polynomial response surface model, the optimal design point can be searched by the GOA. Notably, compared with the initial sample points of the uncertain variables, the optimal design point is much closer to the generalized design point. This indicates that when the structure is impacted by the combination of uncertain variables represented by this optimal point, it is the most likely situation to cause structural failure in the current. This means that this optimal design point is also an effective sample point for the response surface model to learn the characteristics of structural failure. Therefore, we construct an iterative strategy called the dynamic fitting strategy of the response surface: we initialize the response surface model, find the optimal design point, use the searched optimal design point as the new center point to construct new sample points if the outputs do not reach convergence, retrain the response surface model, and find the new optimal design point with the GOA. In addition, since more samples near the generalized design point are found in the dynamic fitting strategy, the response surface model will be fitted near the generalized design point. Therefore, the response surface does not need to be fitted in the entire feasible domain of the PF, thus reducing the cost of fitting. Figure 4 shows the two-dimensional dynamic fitting strategy of the response surface. After two iterations, the response surface model is very close to the PF near the generalized design point.
new sample points if the outputs do not reach convergence, retrain the response surface model, and find the new optimal design point with the GOA. In addition, since more samples near the generalized design point are found in the dynamic fitting strategy, the response surface model will be fitted near the generalized design point. Therefore, the response surface does not need to be fitted in the entire feasible domain of the PF, thus reducing the cost of fitting. Figure 4 shows the two-dimensional dynamic fitting strategy of the response surface. After two iterations, the response surface model is very close to the PF near the generalized design point.

Implementation of the Proposed Method
The proposed method can be implemented through the following steps, as shown in Figure 5.
Step 1: Create the initial sample points using the Bucher design method. Set x* = (x1, x2, …, xn) as the center of the initial sample points and kfαi as the offset distance.
Step 3: Solve the undetermined coefficients a, bi, ci, dij, and eij of the quadratic polynomials.
Step 4: Construct the quadratic polynomial response surface function in the form of Equation (8).
Step 5: Establish the optimization problem of Equation (4) using the scaling factor of the convex set.
Step 6: Set the hyper-parameters of the GOA. Search for the global optimal solution of Equation (4) by the GOA to obtain the optimal design point x* and reliability index value η of the current iteration step.
Step 7: Compare the reliability index value of the current iteration step with that of the previous iteration step. If the convergence conditions are satisfied (|η(i) − η(i-1)| < 10 −3 , where i is the iteration step), stop the iteration and obtain the reliability index. Otherwise, go to Step 8.
Step 8: Generate new sample points by using the optimal design point x* as the center of Bucher design and calculate the sample FEA value. Then, return to Step 3.

Implementation of the Proposed Method
The proposed method can be implemented through the following steps, as shown in Figure 5.

Examples and Discussion
To verify the feasibility of this method, four examples-two mathematical examples and two engineering structure examples-are given.

Example 1
The nonlinear PF of a structure is [37]: The variation range of the uncertain variables can be expressed by the following nonprobabilistic interval model:  Step 1: Create the initial sample points using the Bucher design method. Set x* = (x 1 , x 2 , . . . , x n ) as the center of the initial sample points and kfα i as the offset distance.
Step 3: Solve the undetermined coefficients a, b i , c i , d ij , and e ij of the quadratic polynomials.
Step 4: Construct the quadratic polynomial response surface function in the form of Equation (8).
Step 5: Establish the optimization problem of Equation (4) using the scaling factor of the convex set.
Step 6: Set the hyper-parameters of the GOA. Search for the global optimal solution of Equation (4) by the GOA to obtain the optimal design point x* and reliability index value η of the current iteration step. Step 7: Compare the reliability index value of the current iteration step with that of the previous iteration step. If the convergence conditions are satisfied (|η(i) − η(i-1)| < 10 −3 , where i is the iteration step), stop the iteration and obtain the reliability index. Otherwise, go to Step 8.
Step 8: Generate new sample points by using the optimal design point x* as the center of Bucher design and calculate the sample FEA value. Then, return to Step 3.

Examples and Discussion
To verify the feasibility of this method, four examples-two mathematical examples and two engineering structure examples-are given.

Example 1
The nonlinear PF of a structure is [37]: The variation range of the uncertain variables can be expressed by the following non-probabilistic interval model: where the position parameters are x 1 = x 2 = 0 and the size parameters are α 1 = α 2 = 1. The Bucher experimental design method is used to construct the initial sample points. The position parameter of each variable is taken as the center point. Twice the value of the size parameter is taken as the offset distance of the model. After three calculation steps, the iteration reaches convergence. The explicit expression of the quadratic polynomial response surface at each iteration step is: where y(X), α, and K in each iteration step are: Step 1:ĝ (X) ( Step 2:ĝ (X) ( Step 3:ĝ (X) (2) Through the above calculation, the response surface image of each iteration step can be obtained, as shown in Figure 6. Comparing the response surface with the PF and comparing the optimal design point of each iteration step with the generalized design point, it can be seen that the method in this paper can quickly fit the PF and can quickly find the generalized design point. Notably, image (c) is the part in the red box in image (d). This indicates that the proposed method performs local dynamic fitting to the PF in the vicinity of the domain where the generalized design point is located. The proposed method does not focus too greatly on areas far from the generalized design point. Therefore, the fitting process for the PF can be further accelerated on the premise that the generalized design point can be accurately found.  Table 1 gives the optimal design points and reliability indices searched by the GOA in each iteration step. Table 2 shows the calculation results of different methods.   Table 1 gives the optimal design points and reliability indices searched by the GOA in each iteration step. Table 2 shows the calculation results of different methods. In Table 2, the number of function calls f c can be calculated from the number of sampling points (2n + 1) and the number of iteration steps s. The calculation formula is as follows: The results of the Monte Carlo simulation (MCS) method with a large number of samples are used as relatively accurate solutions. The result of directly using the GOA to solve the optimization problem of Equation (4) is used as a comparison. Tables 1 and 2 show two valuable pieces of information. First, the GOA method has higher precision, and the relative error is only 0.015%. Second, the relative error of the proposed method is close to that of the GOA method; however, the number of fitness function calls is significantly lower, only 1/34 that of the GOA method. We place tighter constraints on the convergence conditions of dynamic iterations in our method (|η(i) − η(i − 1)| < 10 −6 ). The obtained reliability index value is 2.5844, the relative error is reduced to 0.008%, and the number of function calls is 226. It shows that the method in this paper can further increase the accuracy, with a moderate increase in the amount of computation, while still ensuring the computational efficiency. It can be proven that (1) the GOA has a strong optimization capability and can deal with non-probabilistic reliability optimization problems with high precision; (2) the proposed method combines the RSM with the GOA and can improve the efficiency of the solution by reducing the number of calls to the fitness function. The method proposed in this paper has some applicability in solving the non-probabilistic reliability index with the whole interval model of the nonlinear PF of structural reliability.

Example 2
The PF of a structure is Ref. [37]: The variation range of the uncertain variables satisfies the following mixed model, which combines the non-probabilistic interval model and super ellipsoid model: where the position parameters are X = (x 1 , x 3 , x 3 ) = (10, 5, 7) and the size parameters are α 1 = α 2 = 1. It is worth noting that GOA has a total of six hyperparameters: the population size N, the maximum number of iterations t, the upper and lower bounds ul and fl of the search, and the maximum and minimum adaptive shrinkage parameters c max and c min . The upper and lower search bounds ul and fl represent the search space of the model, which can be reasonably estimated according to the variation range of random variables. The population size N and number of iterations t have similar effects on the operation of the algorithm: when N or t increases, the search operator searches more thoroughly in the entire model space, but if they are too large, it will cause a lot of unnecessary searches, which greatly increases the running time of the algorithm and reduces the efficiency. If N or t is too small, the global optimal solution cannot be searched, which affects the accuracy of the understanding. Therefore, N and the maximum t are two empirical parameters. In general, the more complex the optimization problem to be solved, the greater the N, and the greater the t. The adaptive shrinkage parameter c gradually decreases with the increase in the number of iterations, so that the grasshopper algorithm will not quickly converge to the target. In the early stage of the algorithm, a larger c is set, so that the grasshopper operator is encouraged to perform global optimization, and a small c is set in the later stage of the optimization, so that the grasshopper operator can perform local optimization as much as possible and move toward the target. In addition, the problem dimension D is equal to the number of uncertain variables contained in the optimization problem. Subsequent examples follow the above principles.
The Bucher experimental design method is used to construct the initial sample points. The iteration results of Example 2 are shown in Table 3. The comparison results of using different methods are shown in Table 4, and Figure 7 shows the final response surface g(x) = 0 and the PFĝ(>x) = 0. The comparison between the final searched optimal point and the generalized design point is also shown. After four calculation steps, the iteration reaches convergence. In Table 3, with the increase in the iteration steps, the optimal design points and reliability index values tend to become stable gradually. In Table 4, the MCS method is used to simulate 10 6 calculation results as a relatively accurate solution, and the result of directly using the GOA to solve the optimization problem of Equation (4) is used as a comparison. The error of the proposed method is only 0.353%, and the calculation accuracy is satisfactory. The calculation accuracy of the proposed method is similar to that of the GOA, and the number of function calls is only 29 while the number of function calls of the GOA is 506. The calculation efficiency of the proposed method is obviously higher. Furthermore, Figure 7 shows the dynamic fitting process of the response surface model and the optimal design point of each step. Notably, the range shown in the image in step 4 is smaller than those of the other three steps, and the last step shows that this dynamic fitting strategy makes the response surface model fit the PF very well around the generalized design point. The GOA can quickly determine the optimal design point and makes the iteration quickly reach the convergence condition in four steps. To check the fitting ability of this quadratic polynomial response surface model to the actual PF, after each iteration step, the error between the PF and the response surface model at the optimal design point is calculated (see Figure 8). It can be seen from Figure  8 that the sample points set at the beginning have a large error in the positive direction, and then the dynamic fitting strategy of the proposed method quickly adjusts the parameters, and even a negative error occurs; however, the absolute error is reduced. Finally, according to the errors of the first two steps, the fitting error of the response surface is quickly corrected and reduced, which shows that the surrogate model can fit the implicit PF well. To check the fitting ability of this quadratic polynomial response surface model to the actual PF, after each iteration step, the error between the PF and the response surface model at the optimal design point is calculated (see Figure 8). It can be seen from Figure 8 that the sample points set at the beginning have a large error in the positive direction, and then the dynamic fitting strategy of the proposed method quickly adjusts the parameters, and even a negative error occurs; however, the absolute error is reduced. Finally, according to the errors of the first two steps, the fitting error of the response surface is quickly corrected and reduced, which shows that the surrogate model can fit the implicit PF well.

Example 3: A 10-Truss Girder
Example 3 is a 10-truss girder, as shown in Figure 9. An interval model is used to solve the reliability index with an implicit PF. The length of the horizontal and vertical bars is L, the cross-sectional area of each bar is A i (i = 1, 2, . . . , 10), the elastic modulus is E, and the acting positions of the external loads P 1 , P 2 , and P 3 are shown in the figure. Suppose L, A i (i = 1, 2, . . . , 10), E, P 1 , P 2 , and P 3 are 15 relatively independent interval variables, the position and size parameters of which are shown in Table 5.

Example 3: A 10-Truss Girder
Example 3 is a 10-truss girder, as shown in Figure 9. An interval model is used to solve the reliability index with an implicit PF. The length of the horizontal and vertical bars is L, the cross-sectional area of each bar is Ai (i = 1, 2, …, 10), the elastic modulus is E, and the acting positions of the external loads P1, P2, and P3 are shown in the figure. Suppose L, Ai (i = 1, 2, …, 10), E, P1, P2, and P3 are 15 relatively independent interval variables, the position and size parameters of which are shown in Table 5.

Example 3: A 10-Truss Girder
Example 3 is a 10-truss girder, as shown in Figure 9. An interval model is used to solve the reliability index with an implicit PF. The length of the horizontal and vertical bars is L, the cross-sectional area of each bar is Ai (i = 1, 2, …, 10), the elastic modulus is E, and the acting positions of the external loads P1, P2, and P3 are shown in the figure. Suppose L, Ai (i = 1, 2, …, 10), E, P1, P2, and P3 are 15 relatively independent interval variables, the position and size parameters of which are shown in Table 5.    The vertical displacement D(x) of node 2 is considered the intuitive stiffness index in this structure. Therefore, a PF can be established, with D(x) not exceeding 3% of L: The vertical displacement D(x) of node 2 is calculated by the FEA software ANSYS. The parameters of the GOA are set as follows: the dimension is D = 15, the population size is N = 60, the maximum and minimum adaptive shrinkage parameters are c max = 1 and c min = 1.0 × 10 −5 , the number of iterations is t = 5000, the upper and lower bounds of the model space are ul = 10 and fl = 0.0, and the search range of the scaling factor λ is [0, 10].
The calculation results are shown in Table 6. In the table, the search result of the genetic algorithm is taken as the relatively accurate solution [37]. The calculation results of the weighted linear RSM [37] and the ANN [37] are used as a comparison. The result of the weighted linear RSM is the highest, and the accuracy of the proposed method is similar to that of the ANN. However, the number of function calls of the proposed method is significantly lower than those of the above two methods, indicating that the number of FEA calls used in the calculation process of the proposed method is less, so the calculation speed is higher, and the method is more efficient. In another paper by the author of this paper [59], the use of SVM to build a response surface model was proposed and applied to the same structural example. Note that the quadratic polynomial requires fewer samples than SVM, the regression effect is more stable, and the support vector fitting effect is too dependent on the hyperparameters.

Example 4: A 23-Truss Girder
Example 4 is a 23-truss girder, as shown in Figure 10 [60]. A mixed model is used to solve the reliability index with an implicit PF. E is the elastic modulus, A is the crosssectional area, the load P is the uncertain variable of the convex model, and the parameter distributions of the uncertain variables are shown in Table 7.
The vertical displacement D(x) of node 2 is calculated by the FEA software ANSY The parameters of the GOA are set as follows: the dimension is D = 15, the population siz is N = 60, the maximum and minimum adaptive shrinkage parameters are cmax = 1 and cm = 1.0 × 10 −5 , the number of iterations is t = 5000, the upper and lower bounds of the mod space are ul = 10 and fl = 0.0, and the search range of the scaling factor λ is [0, 10].
The calculation results are shown in Table 6. In the table, the search result of the g netic algorithm is taken as the relatively accurate solution [37]. The calculation results the weighted linear RSM [37] and the ANN [37] are used as a comparison. The result the weighted linear RSM is the highest, and the accuracy of the proposed method is sim lar to that of the ANN. However, the number of function calls of the proposed method significantly lower than those of the above two methods, indicating that the number FEA calls used in the calculation process of the proposed method is less, so the calculatio speed is higher, and the method is more efficient. In another paper by the author of th paper [59], the use of SVM to build a response surface model was proposed and applie to the same structural example. Note that the quadratic polynomial requires fewer sam ples than SVM, the regression effect is more stable, and the support vector fitting effect too dependent on the hyperparameters.

Example 4: A 23-Truss Girder
Example 4 is a 23-truss girder, as shown in Figure 10 [60]. A mixed model is used solve the reliability index with an implicit PF. E is the elastic modulus, A is the cros sectional area, the load P is the uncertain variable of the convex model, and the paramet distributions of the uncertain variables are shown in Table 7.   Figure 10. Twenty-three-truss girder. Interval model kN 49 7.5 The vertical displacement D(x) of the midpoint is considered the intuitive stiffness index in this structure. Therefore, the performance can be established, with D(x) not exceeding 0.15 m: The vertical displacement D(x) of the midpoint is calculated by the FEA software ANSYS. The parameters of the GOA are set as follows: the dimension is D = 10, the population size is N = 80, the maximum and minimum adaptive shrinkage parameters are c max = 1 and c min = 1.0 × 10 −5 , the number of iterations is t = 6000, the upper and lower bounds of the model space are ul = 20.0 and fl = 0.0, and the search range of the scaling factor λ is [0, 10].
In Table 8, the MCS method is used to simulate 10 6 calculation results as a relatively accurate solution, and the result of directly using the quadratic polynomial RSM based on sequential quadratic programming (RSM-SQP, implemented using the MATLAB toolbox) to solve the optimization problem of Equation (4) is used as a comparison. It can be seen that in terms of both the number of function calls and the accuracy of the reliability index, the proposed method is better than RSM-SQP.

Conclusions
This paper presents a grasshopper optimization algorithm-based RSM for non-probabilistic structural reliability analysis with an implicit PF. The main conclusions are as follows:

1.
By conducting the proposed method, the reliability index of a structure can be quickly found with fewer function calls. The error self-correction fitting of the response surface model near the generalized design point is achieved. This significantly reduces the number of calls to the PF in the calculation process and avoids the decrease in the solution speed caused by the global fitting of the response surface model. Compared with MCS and other methods, the proposed method is more suitable for performing non-probabilistic reliability analysis with an implicit PF, especially for non-probabilistic reliability analysis of complex engineering structures with high computational accuracy requirements; however, it has lower time consumption because there are fewer calls to the FEA.

2.
The proposed method is nonintrusive. An existing structural analysis program or FEA software can be directly invoked to carry out the proposed method. It can be easily applied by structural engineers without any modification to structural analysis programs or FEA software.
The efficiency and convenience of the proposed method make it a new way of quickly solving reliability problems for structures with implicit PFs. However, this research is suitable only for reliability analysis with a unique optimal solution. The case of multiple solutions needs further research. The method in this paper is not suitable for structures with multiple optimal solutions or structures with discontinuous PFs. In the above cases, the proposed method will not be able to obtain all possible results. In addition, in the reliability calculation of more complex engineering structures, the quadratic polynomial response surface model has great limitations, and it is impossible to fit the high-order features of the structural limit state surface from the sample points. Therefore, the following work will focus on applying different response surface models with higher fitting accuracy and more efficient optimization methods to the dynamic fitting framework proposed in this paper.