A Transformation-Based Improved Kriging Method for the Black Box Problem in Reliability-Based Design Optimization

Abstract

In order to overcome the drawbacks of expensive function evaluation in the practical reliability-based design optimization (RBDO) problem, researchers have proposed the black box-based RBDO method. The algorithm flow of the commonly employed RBDO method for the black box problem consists of the outer construction loop of the surrogate model of the constraint function and the inner surrogate model-based solving loop. To improve the solving ability of the black box RBDO problem, this paper proposes a transformation-based improved kriging method to increase the effectiveness of the two loops identified above. For the outer loop, a sample distribution-based learning function is suggested to improve the construction efficiency of the surrogate model of the constraint function. For the inner loop, a paired incremental sample-based limit reliability boundary construction approach is suggested to transform the RBDO problem into an equivalent deterministic design optimization problem that can be efficiently solved by classical optimization algorithms. The test results of five cases demonstrate that the proposed method can accurately construct the surrogate model of the constraint function and efficiently solve the black box RBDO problem.

1. Introduction

In traditional engineering, researchers overlook uncertainty during design and optimization, also called deterministic design optimization (DDO), which occasionally leads to unexpected failures. To address the uncertainty in the design, researchers have proposed reliability-based design and optimization (RBDO).

A standard RBDO is a double-loop optimization algorithm. The inner loop is a reliability analysis loop, while the outer loop is a criteria-based optimization loop [1]. In order to further improve the solution efficiency of the RBDO algorithm, some researchers have proposed many decoupling methods, such as the single-loop approach (SLA), the safety factor approach (SFA), sequential approximate programming (SAP), and sequential optimization and reliability assessment (SORA) [2,3,4]. Furthermore, researchers have further improved the above decoupling method to significantly enhance the practicality and generalizability of the RBDO method. Meng et al. [5,6] modified the reliability analysis process to solve nonlinear problems and proposed an efficient and robust AC-SLA. Wu et al. [7] designed a safety factor for the constraint function and optimized it continuously in iteration. Yi et al. [8,9] proposed the performance measure approach-based SAP and compared it to the reliability index approach-based SAP to verify the efficiency of the algorithm in different situations. However, among the above four decoupling methods, the most widely used is SORA, which has had a significant amount of improvements [10,11,12,13]. Other researchers have improved the solving capability of the RBDO method by improving the inner loop reliability analysis. In cases where the expression of the RBDO problem is known, researchers have proposed a large number of criteria for reliability analysis [14,15,16,17].

However, in practical engineering, information on the constraint function of the RBDO problem is usually unclear and the cost of evaluation is high. Therefore, researchers have considered the function as a black box and have introduced the surrogate model technique, as well as proposed the corresponding RBDO method for the black box problem. The commonly employed surrogate model techniques in the black box RBDO problem are the kriging [18,19], radial basis function (RBF) [20], support vector machine (SVM) [21,22], and artificial neural network (ANN) [23,24]. Meanwhile, some sampling strategies have been introduced into the construction of the surrogate model, significantly improving the construction efficiency of the surrogate model during the solution of the RBDO problem; such strategies include the directional method [25], constraint boundary sampling (CBS) [26,27], and importance sampling (IS) [28]. Furthermore, the introduction of the surrogate model improves the solving ability of the RBDO method and solves many problems that are difficult to solve using the traditional RBDO method, such as curse of dimensionality [29], highly nonlinear functions [30], and problems with time variants [31].

The surrogate model-based RBDO method usually consists of an outer surrogate model construction loop and an inner surrogate model-based solving loop. For the outer surrogate model construction loop, the model construction strategy is always a hot topic. However, it has been suggested that the existing commonly employed construction strategy, also called learning function (LF), ignores the sample distribution characteristics, which may lead to insufficient model construction efficiency [32]. For the inner surrogate model-based solving loop, it has been found that the existing methods are not effective in dealing with designs with complex failure modes, which may lead to inefficiency of the iteration [33]. In order to solve the above two problems, a transformation-based improved kriging method is proposed in this paper. For the outer surrogate model construction loop, a distribution-based LF is suggested to adjust the sampling process of the sample of the surrogate model, which effectively improves the construction efficiency and increases the model accuracy. For the inner surrogate model-based solving loop, a paired incremental sample-based transformation approach is suggested, which transforms the inner RBDO problem into an equivalent DDO problem for fast solution. This equivalent transformation effectively avoids the effect of complex failure modes on RBDO. In summary, this paper provides an efficient learning function for the surrogate model construction of the constraint function. Meanwhile, a decoupled solution framework for the black box RBDO problem based on equivalent transformation is proposed.

The remainder of this paper is arranged as follows: In Section 2, the black box RBDO problem and the factors affecting the efficiency of its solving are introduced. In Section 3, sample distribution-based LF for constraint function surrogate model construction is described. Meanwhile, several test functions are employed to examine the efficiency of the LF in solving the black box problem. In Section 4, a paired incremental sample-based transformation approach for the inner RBDO problem is presented. In Section 5, the algorithm flow of the proposed transformation-based improved kriging method for the black box problem in RBDO is elaborated. In Section 6, five test cases are employed to verify whether the proposed method is effective. All test cases are treated as the black box problem. In Section 7, some achievements of this paper are summarized.

2. The Black Box Problem in RBDO

In general, a reliability-based design optimization problem can be expressed via the following mathematical form:

$$\begin{array}{ll}find& X\\ min& f\left(X\right)\\ s.t.& P\left\{G\left(X\right)\le 0\right\}\le \mathsf{\Phi }\left(-{\beta }_0\right)\end{array}$$

(1)

where $X$ represents the design variable, $f\left(\cdot \right)$ represents the objective function, $G\left(\cdot \right)$ represents the constraint function, and ${\beta }_0$ represents the objective reliability index.

As mentioned in the Introduction, the problem in Equation (1) can be solved efficiently using the traditional RBDO method when the expression of the constraint function is known. However, in practical engineering, the expression of the constraint function may be unknown, or evaluation may be expensive. Under such circumstances, these practical problems are treated as black box RBDO problems [21] and are solved according to the framework shown in Figure 1.

In the above framework, the solution efficiency of the algorithm is mainly influenced by the following factors:

(1) The efficiency of constructing the surrogate model of the constraint function

The sampling strategy of the training sample directly influences the efficiency of the surrogate model construction. The way of obtaining additional samples based on the optimal design in the RBDO loop is the core of the strategy formulation. Under the appropriate strategy, only a small number of training samples are needed to construct an accurate model where the final reliability design is located.

(2) The efficiency of solving the surrogate model-based RBDO problem

The selection of the initial design and optimization criteria of the design are the main factors affecting the efficiency of solving the surrogate model-based RBDO problem. An inappropriate initial design may increase the iteration step of the optimization. The simple criterion may cause the iteration of the design to be misled by the nonlinear constraint function model, resulting in the design oscillating in the local region [33].

Therefore, in this paper, the two key factors affecting the efficiency of solving the black box RBDO problem are studied and corresponding improvement approaches are proposed.

3. Surrogate Model Construction with the Distribution-Based Sampling Strategy

In this section, the sampling strategy for the surrogate model construction is described, which introduces the distribution factor into the classical strategy and decreases the negative impact of uneven sample distribution on the efficiency of model construction.

The efficiency global optimization algorithm (EGO) [34] treats the prediction results of the kriging-based surrogate model as though they obey a normal distribution. This idea is widely used in RBDO and has been improved, forming a series of adaptive kriging methods [23,35,36,37]. As mentioned above, the core part of the adaptive kriging method is the sampling strategy, which is also called the learning function (LF) in research on the RBDO [38,39]. The most commonly employed LF is expressed as follows:

$$\mathrm{LF}\left(x\right)=\frac{\left|u\left(x\right)\right|}{mse\left(x\right)}$$

(2)

where $u\left(\cdot \right)$ is the prediction result of the kriging model and $mse\left(\cdot \right)$ is the prediction error evaluated by the kriging model.

Employing the LF can effectively search for an additional sample, which has significant prediction errors near the limit state boundary of the surrogate model of the constraint function. However, due to the different nonlinearities in the location of the samples, the level of prediction error near the samples is different. For example, taking the estimation error of 90% as the threshold, the trust regions of different samples in the normal space of one case are shown in Figure 2, where $U_x$ represents the projection coordinate of the sample in the normal space.

Ignoring the level of prediction error of different samples may lead to a concentration of training samples in the local nonlinear region, which further leads to a reduction in the efficiency of model construction [32]. Therefore, based on the authors’ previous research [40], a sample distribution balancing approach is introduced to form the distribution-based sampling strategy. The mathematical expression of the improved learning function LF’ proposed in this paper is as follows:

$${\mathrm{LF}}^{\prime }\left(x\right)=\frac{\left|u\left(x\right)\right|}{mse\left(x\right)}·e^{\left(-\frac{{‖x-T_i‖}^2}{R_i{}^2}\right)}$$

(3)

where $i$ represents the i-th existing training samples, $T_i$ represents the closest existing training sample, and $R_i$ represents the trust radius of the closest existing training sample.

The trust radius of each sample is initialized by the following formula:

$$R_i=‖p-T_i‖$$

(4)

where $p$ represents the nearest location from the training samples, $T_i$, where the level of the prediction error is equal to the objective reliability index. The location, $p$, is obtained by solving the following optimization problem:

$$\begin{array}{ll}find& p\\ min& \left|\left|\frac{mse\left(p\right)}{u\left(T_i\right)}\right|-\mathsf{\Phi }\left(-{\beta }_0\right)\right|\end{array}$$

(5)

where ${\beta }_0$ represents the objective reliability index of the RBDO problem.

With an increase in the number of training samples, the accuracy of the model near the samples gradually improves and the trust radius of the samples also changes adaptively.

To demonstrate that adding sample balancing measures into the LF can improve the efficiency of model construction, several test functions are used here as black box functions for verification. The test functions are listed in

Table 1. Information of the test functions.

Test Function	Expression
Branin	$f\left(x_1,x_2\right)=-0.3979+{\left(x_2-\frac{51}{40{\pi }^2}{x}_1^2+\frac{5}{\pi }{x}_1-6\right)}^2+10\left(1-\frac{1}{8\pi }\right)cos\left(x_1\right)+10$
Shubert	$f\left(x_1,x_2\right)=186.7309+\left({\displaystyle \sum _{i=1}^{5}icos\left(\left(i+1\right)x_1+i\right)}\right)\left({\displaystyle \sum _{i=1}^{5}icos\left(\left(i+1\right)x_2+i\right)}\right)$
Griewank	$f\left(x\right)={\displaystyle \sum _{i=1}^d\frac{x_i^2}{4000}}-{\displaystyle \prod _{i=1}^{d}cos\left(\frac{x_i}{\sqrt{i}}\right)}+1$
Bohachevsky	$f\left(x_1,x_2\right)=x_1^2+2x_2^2-0.3cos\left(3\pi x_1\right)-0.4cos\left(4\pi x_2\right)+0.7$
Levy	$\begin{array}{ll}f\left(x\right)=& sin{\left(\pi {\omega }_1\right)}^2+{\displaystyle \sum _{i=1}^{d-1}{\left({\omega }_i-1\right)}^2\left[1+10sin{\left(\pi {\omega }_i+1\right)}^2\right]}\\ & +{\left({\omega }_d-1\right)}^2\left[1+sin{\left(2\pi {\omega }_d\right)}^2\right]\end{array}$ ${\omega }_i=1+\frac{x_i-1}{4}$

.

The test results with the LF and the LF’ are summarized in

Table 2. Test result comparison between the LF and the LF’.

Test Function	LF		LF’
	Best	Iteration	Mean	Variation	Worst	Iteration
Branin	0.0001	30	0.0001	3.3121 × 10−6	0.0063	27
Shubert	−21.9672	17	−0.0105	6.2978 × 10−5	−0.0410	18
Griewank	0.0021	52	0.0012	1.1359 × 10−6	0.0017	47
Bohachevsky	0.0209	8	0.0010	2.7831 × 10−6	0.0043	5
Levy	0.0152	24	0.0020	5.0631 × 10−6	0.0058	22

. The following conclusions can be drawn from analyzing the data in the table: (1) Compared to the LF strategy, employing the LF’ strategy reduces the number of samples that are required to construct a surrogate model of the constraint function; (2) compared to the LF strategy, the design obtained with the LF’ strategy is more accurate; (3) upon repeated testing, the designs obtained with the LF strategy show slight variation and the results are robust. In summary, the LF’ strategy can improve the efficiency of the construction of the surrogate model of the constraint function. The first factor affecting the efficiency of solving the black box RBDO problem mentioned in the previous section is addressed.

4. A Paired Incremental Sample-Based Transformation Method for Surrogate Model-Based RBDO

The LF strategy proposed in Section 2 effectively improves the efficiency of the construction of the surrogate model of the constraint function. In this section, the aim is to improve the efficiency of solving the surrogate model-based RBDO problem.

As mentioned above, most of the current methods for the surrogate model-based RBDO problem employ the nested solution framework, consisting of reliability analysis and criteria-based design optimization, resulting in the efficiency being limited by the selection of the initial design and the optimization criteria of the design. To remove this limitation, a paired incremental sample-based RBDO method is proposed in this paper, which decouples the nested reliability analysis and design optimization by directly constructing a limit reliability boundary, where the reliability index satisfies the optimization requirement. Furthermore, with the help of the constructed boundary, the original RBDO problem is transformed into an equivalent DDO problem.

The limit reliability boundary is constructed based on the samples of the surrogate model of the constraint function and the paired incremental samples. Moreover, paired incremental samples are obtained based on the reliability analysis results of the sample of the surrogate model of the constraint function. According to the surrogate model of the constraint function and following the reliability analysis formula, the most probable point (MPP) of the sample $x$ can be obtained:

$$\begin{array}{ll}find& U_{MPP}\\ min& ‖U_{MPP}‖\\ s.t.& \hat{G}\left(U_{MPP}\right)=0\end{array}$$

(6)

where $U_{MPP}$ represents the MPP coordinate in normal space.

The modal length of the MPP in Equation (6) is also known as the reliability index of the sample $x$ and labeled “${\beta }_x$”. The MPP represents the minimum distance point on the constraint function to the sample $x$ in normal space. According to the geometry, the minimum distance point from any point located in the vector direction from the sample $x$ to the MPP to the constraint function is the MPP. Therefore, it is possible to intercept any reliability index in the vector direction of the sample $x$ to the MPP.

Interception of the paired incremental samples proposed in this paper is shown in Figure 3. After obtaining the MPP of the sample $x$, in the vector direction of the MPP, paired incremental samples with a reliability index $\left({\beta }_0+\epsilon \right)$ and $\left({\beta }_0-\epsilon \right)$ are obtained equally spaced, with $\epsilon$ as the spacing. By obtaining several sets of incremental samples, a boundary model satisfying the objective reliability index can be constructed. However, Figure 3 only shows a case in which the sample reliability index is larger than the objective reliability index; the situation is more complicated in practical application. Therefore, it is necessary to develop an appropriate paired-sample interception strategy to meet different reliability analysis results.

According to the relationship between the reliability index of the incremental samples and the objective reliability index, the paired samples are divided into an upper incremental sample, $P_{+}$, and a lower incremental sample, $P_{-}$. Based on the reliability index of the sample $x$, the interception of upper incremental sample $P_{+}$ is classified into three cases, as shown in Equation (7):

$$P_{+}=\left\{\begin{array}{ll}{U}_{MPP}\left(\frac{{\beta }_0+\epsilon }{{\beta }_X}\right)& {\beta }_X\ge {\beta }_0+\epsilon \\ U_{MPP}& {\beta }_0\le {\beta }_X<{\beta }_0+\epsilon \\ \varnothing & {\beta }_X<{\beta }_0\end{array}\right.$$

(7)

where $U_{MPP}$ represents the MPP in normal space, $\epsilon$ represents the interception spacing, and ${\beta }_X$ represents the reliability index of the sample $x$.

Based on the reliability index of the sample $x$, the interception of the lower incremental sample $P_{-}$ is classified into three cases, as shown in Equation (8):

$$P_{-}=\left\{\begin{array}{ll}{U}_{MPP}& {\beta }_{\mathbf{X}}\le {\beta }_0-\epsilon \\ U_{MPP}\left(\frac{{\beta }_0-\epsilon }{{\beta }_{\mathbf{X}}}\right)& else\end{array}\right.$$

(8)

From Equations (7) and (8), it can be seen that the generation of incremental samples depends on the reliability index of the sample $\mathbf{x}$ and interception spacing. In order to appropriately set the interception spacing, two cases were employed to examine the effect of different interception spacings on the accuracy of the limit reliability boundary. The examination process was as follows: (1) random sampling in the feasible domain; (2) constructing the corresponding limit reliability boundary based on random samples and generated paired incremental samples based on different interception spacings, $\epsilon =\left[\sigma ,2\sigma ,3\sigma ,4\sigma ,5\sigma \right]$; (3) resampling randomly and obtaining the true reliability index and the limit reliability boundary-based reliability index; (4) calculating the mean error and the variance of the error between the true reliability index and the model-based reliability index, which are displayed in Figure 4. In addition, $\sigma$ represents the variance of the design variable.

Figure 4 illustrates that different interception spacings have an effect on the accuracy of the limiting reliability boundary. The interception spacing showed consistently trending effects on the mean and variance of the errors in different test cases. With an increase in the interception spacing, the mean and variance of the errors first decreased and subsequently increased. When the interception spacing was two or three times the variance of the design variable, the accuracy of the limit reliability boundary was higher compared to other cases.

This section transforms the constraint function surrogate model-based RBDO problem into a model-based equivalent DDO problem by constructing a limit reliability boundary. This transformation decouples the traditional nested RBDO process and avoids the selection of an initial design, which can effectively improve the solving efficiency of the RBDO loop.

5. Optimization Flow of the Proposed Method

In Section 3 and Section 4, a distribution-based sampling strategy for refinement of the kriging model of the constraint function and a paired incremental sampling-based method for transformation of the RBDO problem were described. These approaches can effectively address the two main factors that decrease the efficiency of solving the black box RBDO problem, respectively. In this section, the Algorithm 1 flow of the proposed method for the black box RBDO problem is elaborated as follows, and the invocation of the distribution-based sampling strategy and the paired incremental sample-based transformation method in the flow is shown in Figure 5.

Algorithm 1 Algorithm flow for black box RBDO problem

Step 1: Activate Latin hypercube sampling (LHS) to obtain the initial samples. The number of the initial sample is $5N$, where $N$ represents the dimension number of the design variable. Step 2: Calculate the response value of the samples and construct the kriging surrogate model of the constraint function. Step 3: Based on the constructed kriging surrogate model of the constraint function, a paired incremental sample-based transformation method is generated to transform the RBDO problem into an equivalent DDO problem. The whole process is divided into four substeps. Step 3.1: Analyze the reliability of the samples employed in the construction of the kriging surrogate model of the constraint function and obtain the reliability index and MPP of each sample. Step 3.2: According to the generation strategy of the paired incremental sample, obtain the paired incremental sample of each sample. Step 3.3: Employ the paired incremental samples and the samples of the constraint function model to construct a limit reliability boundary. Step 3.4: Transform the RBDO problem into an equivalent DDO problem using the limit reliability boundary. Step 4: Obtain the optimal design of the equivalent DDO problem with particle swarm optimization (PSO). Step 5: Comparing the current optimal value, $Y_t$, and the last step optimum value, $Y_{t-1}$, if the differential value is smaller than the threshold, $\epsilon$, as per Equation (9), go to Step 9; otherwise, go to Step 6. $$\frac{\left|Y_t-Y_{t-1}\right|}{\left|Y_{t-1}\right|}<\epsilon$$ (9) Step 6: Obtain an additional sample according to the improved learning function LF’. Step 7: Output/replace the current optimal design as the additional sample. Step 8: Add the additional samples from Steps 6 and 7 into the sample set and go back to Step 2. Step 9: Output the optimal result and end the optimization process.

6. Test Examples

In this section, several practical examples are used to verify the availability of the proposed method, which involves nonlinear functions and multidimensional and engineering problems. During the example testing, all constraint functions were considered as black box functions. The function expressions were only used to check the accuracy of the kriging surrogate model constructed by the proposed method. From analysis of the verification results, the proposed method can solve the above problems more effectively than the method cited in the references. The kriging model constructed in this paper employs the Dace model [41]. The regression function and correlation function of the model are set as zero-order polynomial function and Gaussian function, respectively.

6.1. Haupt Example

In this part, the Haupt example [26] is introduced for testing. This numerical example includes two uncertain design variables and is constrained by two functions. One constraint is the trigonometric function, which can verify the ability to solve nonlinear problems. The mathematical expression of the Haupt example is described as follows:

$$\begin{array}{cc}find\hfill & \mathbf{x}={\left[x_1,x_2\right]}^{\mathrm{T}}\hfill \\ min\hfill & f={\left(x_1-3.7\right)}^2+{\left(x_2-4\right)}^2\hfill \\ s.t.\hfill & P\left[G_i\left(\mathbf{x}\right)\le 0\right]\le \mathsf{\Phi }\left(-{\beta }_i\right)\hspace{1em}i=1,2\hfill \\ & 0\le x_1\le 3.5,0\le x_2\le 3.5\hfill \\ & G_1\left(\mathbf{x}\right)=-x_{1}sin\left(4x_1\right)-1.1x_{2}sin\left(2x_2\right)\hfill \\ & G_2\left(\mathbf{x}\right)=x_1+x_2-3\hfill \\ & x_i\sim \left(x_i,{0.1}^2\right)fori=1,2\hfill \\ & {\beta }_1={\beta }_2=2\hfill \end{array}$$

(10)

Through 11 iterations, we obtained the optimization results of this example. In Figure 6, these four figures show the construction of the kriging model of the constraint function and limit reliability boundary. In the figure, the actual boundary of the constraint function is shown by the solid green line, the surrogate model of the constraint function is shown by the solid red line, and the surrogate model of the limit reliability boundary of the constraint function is shown by the dashed blue line. In order to distinguish the surrogate model of the two constraint functions, we used the dashed line to denote the surrogate model of the first constraint function and the dotted dashed line to denote the surrogate model of the second constraint function.

Because one of the constraint functions is linear, kriging models of the constraint function and limit reliability boundary were constructed at the beginning of the iteration. As the iteration progressed, the kriging model of another constraint function was gradually constructed. As can be seen from the figure, the early stage of the algorithm was mainly to improve the fitness of the constraint function kriging model, and the construction of the limit reliability boundary was relatively rough. From the third to the fourth pictures, the reliability boundary near the optima was well constructed.

In order to make the verification process more convincing, this paper selected the standard SORA, SORA with Quantiles [42], and IPFR [43] methods for comparison. We compared the final calculated value, the optimal point location, the call times of the constraint functions, and the probability performance of the optimum on each constraint function. The comparison results are listed in

Table 3. Comparison results summary of the Haupt example.

Method	Result	Optimal Point	Function Call	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$
Proposed method	1.3055	2.8196, 3.2717	44	2.0053	Inf
SORA	1.3188	2.8210, 3.2652	509	2.0366	Inf
SORA with Quantiles	1.3039	2.8241, 3.2674	52	1.9987	Inf
IPFR	1.3070	2.8215, 3.2681	497	2.0133	Inf

.

Table 3 shows that the optimal values obtained by the testing methods were similar. In the specific value, 1.3055 was the minimum, except for the SORA with Quantiles method. However, the reliability index of the best calculation result obtained by the SORA with Quantiles method was only 1.9987 for the first constraint function. Although this value is within the allowable tolerance of error, it is not as good as the calculation result of the proposed method in a practical situation, where the reliability index was 2.0053. IPFR’s performance was similar to that of SORA.

6.2. Choi Example

In the same reference [26], researchers employed the Choi example, which has two uncertain design variables and three constraint functions. The mathematical expression of the Choi example is described as follows:

$$\begin{array}{cc}find\hfill & \mathbf{x}={\left[x_1,x_2\right]}^{\mathrm{T}}\hfill \\ min\hfill & f=x_2-0.4x_1\hfill \\ s.t.\hfill & P\left[G_i\left(\mathbf{x}\right)\le 0\right]\le \mathsf{\Phi }\left(-{\beta }_i\right)i=1,2,3\hfill \\ & 0\le x_1\le 10,0\le x_2\le 10\hfill \\ & G_1\left(\mathbf{x}\right)=\frac{x_1^2{x}_2}{20}-1\hfill \\ & G_2\left(\mathbf{x}\right)={\left(x_1-6.5\right)}^2+8\left(x_2-3\right)\hfill \\ & G_3\left(\mathbf{x}\right)=\frac{80}{x_1^2+8x_2+5}-1\hfill \\ & x_i\sim \left(x_i,{0.6}^2\right)fori=1,2,3\hfill \\ & {\beta }_1={\beta }_2={\beta }_3=1\hfill \end{array}$$

(11)

Figure 7 shows the construction process of the surrogate model in this case. Similar to the previous case, the limit reliability boundary surrogate model with different constraint functions is used with different labels. The first constraint function is denoted by the dashed line, the second by the dotted line, and the third by the dotted line. In this case, the nonlinearity of the three constraint functions is weaker than that of the first constraint function in the previous case. Therefore, compared to case 6.1, the number of samples used in the surrogate model construction is smaller and the algorithm converges faster.

The same evaluation method was employed as that in the Haupt example to evaluate the proposed method’s ability in the Choi example. In

Table 4. Comparison results summary of the Choi example.

Method	Result	Optimal Point	Function Call	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{\beta }}_3$
Proposed method	1.1399	6.1085, 3.5833	68	4.9710	1.0000	1.0000
SORA	1.1390	6.1105, 3.5832	480	4.9717	0.9995	0.9972
SORA with Quantiles	1.1451	6.1064, 3.5877	87	4.9772	1.0075	0.9992
IPFR	1.1403	6.1114, 3.5849	492	4.9748	1.0022	0.9945

, the comparison results are shown. In terms of optimal values, the proposed method (1.1399) and SORA (1.1390) were better than the other two methods, with values of 1.1451 and 1.1403, respectively. At the same time, it still demonstrated considerable advantage in constraining function calls.

6.3. Hock and Schittkowski Problem

In this subsection, a test example from Hock and Schittkowski was selected. This problem has 10 uncertain design variables and eight constraint functions. The mathematical expression of the Reducer design example is described as follows:

$$\begin{array}{cc}find\hfill & \mathbf{x}={\left[x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9,x_{10}\right]}^{\mathrm{T}}\hfill \\ min\hfill & f=x_1^2+x_2^2+x_1{x}_2-14x_1-16x_2+{\left(x_3-10\right)}^2+\hfill \\ & 4{\left(x_4-5\right)}^2+{\left(x_5-3\right)}^2+2{\left(x_6-1\right)}^2+5x_7+7{\left(x_8-11\right)}^2+\hfill \\ & 2{\left(x_9-10\right)}^2+{\left(x_{10}-7\right)}^2+45\hfill \\ s.t.\hfill & P\left[G_i\left(\mathbf{x}\right)\ge 0\right]\le \mathsf{\Phi }\left(-{\beta }_i\right)i=1,2,3,4,5,6,7,8,9,10\hfill \\ & x_i\ge 0i=1,2,3,4,5,6,7,8,9,10\hfill \\ & G_1\left(\mathbf{x}\right)=\frac{4x_1+5x_2-3x_7+9x_8}{105}-1\hfill \\ & G_2\left(\mathbf{x}\right)=10x_1-8x_2-17x_7+2x_8\hfill \\ & G_3\left(\mathbf{x}\right)=\frac{-8x_1+2x_2+5x_9-2x_{10}}{12}-1\hfill \\ & G_4\left(\mathbf{x}\right)=\frac{3{\left(x_1-2\right)}^2+4{\left(x_2-3\right)}^2+2x_3{}^2-7x_4}{120}-1\hfill \\ & G_5\left(\mathbf{x}\right)=\frac{5x_1^2+8x_2+{\left(x_3-6\right)}^2-2x_4}{40}-1\hfill \\ & G_6\left(\mathbf{x}\right)=\frac{0.5{\left(x_1-8\right)}^2+2{\left(x_2-4\right)}^2+3x_5^2-x_6}{30}-1\hfill \\ & G_7\left(\mathbf{x}\right)=x_1^2+2{\left(x_2-4\right)}^2+2x_1{x}_2+14x_5-6x_6\hfill \\ & G_8\left(\mathbf{x}\right)=-3x_1+6x_2+12{\left(x_9-8\right)}^2-7x_{10}\hfill \\ & x_i\sim \left(x_i,{0.02}^2\right)fori=1,2,3,4,5,6,7,8,9,10\hfill \\ & {\beta }_1={\beta }_2={\beta }_3={\beta }_4={\beta }_5={\beta }_6={\beta }_7={\beta }_8={\beta }_9={\beta }_{10}=3\hfill \end{array}$$

(12)

In

Table 5. Comparison results of the Hock and Schittkowski problem.

Method	Result	Optimal Point	Function Call
Proposed method	27.7692	2.1327, 2.3363, 8.7063, 5.0909, 0.9284, 1.4581, 1.3842, 9.8050, 8.1495, 8.4750	409
SORA	27.7438	2.1350, 2.3299, 8.7088, 5.1011, 0.9302, 1.4640, 1.3889, 9.8101, 8.1517, 8.4628	1862
SORA with Quantiles	27.7455	2.1380, 2.3231, 8.7054, 5.0945, 0.9221, 1.4489, 1.3950, 9.8152, 8.1542, 8.4525	480
IPFR	27.7548	2.1290, 2.3460, 8.7100, 5.0970, 0.9290, 1.4530, 1.3770, 9.8000, 8.1360, 8.4680	1769

, the optimization comparison results are shown, while in

Table 6. Probabilistic constraints comparison of the Hock and Schittkowski problem.

Method	${\mathit{\beta }}_{6,8}$	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	${\mathit{\beta }}_3$	${\mathit{\beta }}_4$	${\mathit{\beta }}_5$	${\mathit{\beta }}_7$
Proposed method	Inf	3.0374	3.0050	3.0029	3.0823	3.0008	3.0340
SORA	Inf	2.9981	2.9945	2.9816	3.0094	3.0167	3.0083
SORA with Quantiles	Inf	2.9737	3.0158	3.0044	3.0556	3.0041	3.0192
IPFR	Inf	2.9924	3.0102	3.0257	3.0377	2.9910	2.9955

, the reliability index comparison is shown. According to the provided theoretical results, the optimal point is located at the intersection area of six constraint boundaries, which tests the searching ability of the proposed method. As per Table 5, all methods obtained the best solution in terms of error tolerance. Regarding the function call times, SORA called the constraint function 1862 times because of the large number of related constraint boundaries. In contrast, the proposed method called 409 times, while SORA with Quantiles called 480 times. According to the reliability index comparing results in Table 6, the proposed method performed well in the second (3.0050), third (3.0029), and fifth (3.0008) constraint functions. The other two methods only met the constraint requirements within the error tolerance.

6.4. Reducer Design

In order to increase the sufficiency of demonstration, an engineering example of reducer design optimization was selected [4]. This example is a multidimensional example with seven uncertain design variables and 10 constraint functions. Simultaneously, the range scale of the design variables vastly different from one another, which increases the difficulty of constructing an alternative model. The structural diagram of the reducer is shown in Figure 8.

This example can be expressed in mathematical form as follows:

$$\begin{array}{cc}find\hfill & \mathbf{x}={\left[x_1,x_2,x_3,x_4,x_5,x_6,x_7\right]}^{\mathrm{T}}\hfill \\ min\hfill & f=0.7854x_1{x}_2^2\left(3.3333x_3^2+14.9334x_3-43.0934\right)-\hfill \\ & 1.508x_1\left(x_6^2+x_7^2\right)+7.477\left(x_6^3+x_7^3\right)+0.7854\left(x_4{x}_6^2+x_5{x}_7^2\right)\hfill \\ s.t.\hfill & P\left[G_i\left(\mathbf{x}\right)\le 0\right]\le \mathsf{\Phi }\left(-{\beta }_i\right)i=1,2,3,4,5,6,7\hfill \\ & 2.6\le x_1\le 3.6,0.7\le x_2\le 0.8,17\le x_3\le 28,7.3\le x_4\le 8.3,7.3\le x_5\le 8.3,\hfill \\ & 2.9\le x_6\le 3.9,5\le x_7\le 5.5\hfill \\ & G_1\left(\mathbf{x}\right)=\frac{27}{x_1{x}_2^2{x}_3}-1,G_2\left(\mathbf{x}\right)=\frac{397.5}{x_1{x}_2^2{x}_3^2}-1,G_3\left(\mathbf{x}\right)=\frac{1.93x_4^3}{x_2{x}_3{x}_6^4},G_4\left(\mathbf{x}\right)=\frac{1.93x_5^3}{x_2{x}_3{x}_7^4}\hfill \\ & G_5\left(\mathbf{x}\right)=\frac{\sqrt{{\left(745x_4/x_2{x}_3\right)}^2+1.69\times {10}^6}}{0.1x_6^3}-1100\hfill \\ & G_6\left(\mathbf{x}\right)=\frac{\sqrt{{\left(745x_5/x_2{x}_3\right)}^2+1.575\times {10}^8}}{0.1x_7^3}-850\hfill \\ & G_7\left(\mathbf{x}\right)=x_2{x}_3-40,G_8\left(\mathbf{x}\right)=5-\frac{x_1}{x_2},G_9\left(\mathbf{x}\right)=\frac{x_1}{x_2}-12\hfill \\ & G_{10}\left(\mathbf{x}\right)=\frac{1.5x_6+1.9}{x_4}-1,G_{11}\left(\mathbf{x}\right)=\frac{1.1x_7+1.9}{x_5}-1\hfill \\ & x_i\sim \left(x_i,{0.005}^2\right)fori=1,2,3,4,5,6,7\hfill \\ & {\beta }_1={\beta }_2={\beta }_3={\beta }_4={\beta }_5={\beta }_6={\beta }_7=3\hfill \end{array}$$

(13)

The engineering example involved seven dimensions, so we could not use the graph method to show the search process of the algorithm. Therefore, in

Table 7. Comparison results of the Hock and Schittkowski problem.

Method	Result	Optimal Point	Function Call
Proposed method	3040.24	3.5771, 0.7000, 17.0000, 7.3128, 7.7571, 3.3661, 5.3032	396
SORA	3040.69	3.5789, 0.7000, 17.0000, 7.3000, 7.7641, 3.3662, 5.3027	1058
SORA with Quantiles	3039.21	3.5770, 0.7000, 17.0000, 7.3114, 7.7602, 3.3660, 5.3016	451
IPFR	3041.29	3.5785, 0.7000, 17.0000, 7.3092, 7.7639, 3.3671, 5.3034	994

, the comparison results of the optimization are shown. The call times of the constraint function and the constraint index of the optimal value of each constraint function are shown in

Table 8. Probabilistic constraints comparison of the Hock and Schittkowski problem.

Method	${\mathit{\beta }}_{2,3,4,7,9,10}$	${\mathit{\beta }}_1$	${\mathit{\beta }}_5$	${\mathit{\beta }}_6$	${\mathit{\beta }}_8$	${\mathit{\beta }}_{11}$
Proposed method	Inf	6.7029	3.1717	3.0063	3.0241	3.1723
SORA	Inf	6.7363	3.1965	3.2058	3.0947	4.1881
SORA with Quantiles	Inf	6.7010	3.1522	2.9860	3.0202	3.8262
IPFR	Inf	6.7289	3.2730	3.2458	3.0790	4.0576

. First, the optimal results can be obtained by several methods within the error tolerance. Like the previous two examples, the SORA method called the constraint functions up to 1058 times, while the proposed method only called them 396 times. Thus, the proposed method demonstrated good performance on four constraint functions that reached the limit value in terms of the reliability index. The SORA with Quantities method performed the best in the 5th and 8th constraints, but the deviation was larger in the 11th constraint.

From the comparison of the four testing examples above, we can infer the following: The proposed method can solve the RBDO problem effectively. Among the methods compared, the function call times were the lowest. Meanwhile, the obtained optimization results of the proposed method were equivalent to those of the other methods within the error tolerance. In terms of the reliability index, the optimal result was slightly conservative compared to that of the other methods. However, we could achieve the same effect as the other methods by setting the deviation range for the target value of the reliability index.

6.5. Fan of the Automotive Cooling System

To further test the practicability of the algorithm in engineering problems, we used a simulation model from previous work for RBDO. The object of this calculation example is an automobile cooling system fan, and the design parameters are shown in Figure 9. The objective and constraint functions of this example have no mathematical expressions, and all data are directly output from the simulation model.

According to the actual situation, the initial calculation domain inlet was set as the pressure stagnation inlet with a value of –60 Pa. Meanwhile, the calculation domain outlet was set as the atmospheric pressure outlet with a value of 0 Pa. Lastly, the rotation speed of the rotating calculation domain was set as 1200 rpm. Information of the random design variables of the cooling fan is listed in

Table 9. Information of the random design variables of the cooling fan.

Variable	PDF.	Standard Deviation σ
x1	Normal	0.01
x2	Normal	0.01
w	Normal	0.03
R	Normal	0.5
L	Normal	0.02

.

The efficiency of the cooling fan was treated as a constraint condition, and the threshold was set to 30%. By monitoring the parameters of airflow pressure and fan blade torque, the efficiency of the axial fan can be calculated by Equation (14):

$$\eta =\frac{\Delta P\times M\times 9550}{T\times n}$$

(14)

where $\Delta P$ (pa) is the pressure decrease between the inlet and outlet faces, $M$ (kg/s) is the mass flow rate of the fan domain, $T$ (N/m) is the torque of the fan blade, and $n$ (r/s) is the rotation speed of the fan blade. This RBDO problem can be expressed simply as follows:

$$\begin{array}{ll}find& X=[x_1,x_2,w,R,L{]}^{\mathrm{T}}\\ min& -Massflowrat{e}_{inlet}\left(X\right)\\ s.t.& P[G_i\left(X\right)\le 0]\le \mathsf{\Phi }\left(-{\beta }_i\right)\hspace{1em}i=1,2\\ & 13\le x_1\le 30,15\le x_2\le 30,20\le w\le 50,50\le R\le 80,300\le L\le 500\\ & G_1(X)=\eta -0.3\\ & G_2(X)=Torqu{e}_{rotation}-0.8625\\ & X\sim \left(X,{\sigma }_i{}^2\right)\hspace{1em}for\hspace{1em}i=1,2,3,4,5\\ & {\beta }_1={\beta }_2=2\end{array}$$

(15)

To make the computational fluid dynamics simulation converge quickly, the simulation domain was divided into three parts: The inlet domain, the outlet domain, and the fan rotation domain. At the same time, for the inlet and outlet domains in which the flow field is relatively simple, the mesh size was set to 15 mm, while for the rotation domain, the mesh size was set to 5 mm. This measure ensures that the model CFD results can be obtained quickly without distortion in the simulation. The situation of the simulation domain mesh generation is shown in Figure 10.

After 34 iteration steps, the termination condition was satisfied. The model was modified according to the obtained optimal solution and imported into a CFD software for simulation, and the mass flow rate was read. The optimization results of DDO and RBDO are shown in

Table 10. Results comparison of DDO and RBDO.

Result Type	Optimal Result	Mass Flow Rate of the Fan Inlet Domain (Objective)	Reliability Index (MCS)
			${\mathit{\beta }}_1$	${\mathit{\beta }}_2$
DDO	27.89, 15.23, 41.67, 52.74, 437.67	0.6237	1.6356	1.6188
RBDO	27.99, 14.99, 41.65, 52.57, 437.71	0.6194	2.0094	2.0131

, and reliability indexes were obtained and reliability indexes were obtained by the MCS method. The velocity vector diagram of the cross-section in the simulation domain is shown in Figure 11. Through comprehensive analysis of the figure and table, it can be seen that the design point that meets the reliability requirements was successfully found through the method proposed in this paper.

7. Conclusions

In this paper, a transformation-based improved kriging method was proposed for the black box RBDO problem, which commonly occurs in practical engineering, and the solution efficiency of the problem was effectively improved. The improvement was presented for two main factors affecting the solving efficiency of the black box RBDO problem:

(1) A distribution-based improved learning function was proposed for the sample update strategy of the surrogate model construction of the constraint function in the outer loop of the problem solution, which can balance the distribution of additional samples and improve the efficiency of the surrogate model construction.

(2) A paired incremental sample-based transformation method was proposed for the inner RBDO solution. This method constructs a limit reliability boundary based on incremental samples and transforms the problem into an equivalent DDO problem for solution, which can avoid the initial design selection and improve the solution capability for the nonlinear RBDO problem.

The test results showed that the proposed method can reduce the size of the sample for the surrogate model of the constraint function, improve the accuracy of the surrogate model, and obtain the optimal design of the black box RBDO problem more efficiently and accurately.

The method proposed in this paper can solve the general black box RBDO problem. However, for the extreme nonlinear or high-dimensional cases, the efficiency of the surrogate model construction in the proposed method may be reduced, and further improvement of the proposed method is needed to extend its generalizability continuously.