# A Reliability-Based Multidisciplinary Design Optimization Strategy considering Interval Uncertainty Based on the Point-Infilled Kriging Model

Authors: Chao Yuan; Hao Zhang; Yong Sun; Yunhan Ling; Lidong Pan; Yue Sun; Dianyu Fu; Zhengguo Hu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Generally, traditional uncertainty design optimization (UDO) methods are based on probability density distribution function or fuzzy membership function. In this situation, a large amount of uncertain information is necessary to construct the UDO model accurately. While, the interval UDO methods require less design information. Only the upper and lower bounds of interval uncertainty are utilized to construct the optimization model. In this study, to enhance the efficiency and accuracy of UBDO considering interval uncertainty, a reliability-based multidisciplinary design optimization (RBMDO) strategy using the point-infilled Kriging model is proposed. In the given method, a double-nested RBMDO model considering interval uncertainty is established. The collaborative optimization is utilized to deal with coupling relationships among complex systems. Then, the point-infilled Kriging response surface strategy is introduced to approximate the RBMDO model. The procedure of the interval multidisciplinary collaborative optimization method based on the Kriging model is discussed. Two examples are given to illustrate the application of the proposed method.


---

## Body

## 1. Introduction

The traditional equipment and systems are usually constructed by deterministic optimization model and solved by traditional method such as the probability method. However, in the practice process, there are various uncertain factors such as material characteristics, manufacturing errors, and operating environment [1–4]. Dealing with the needs of modern equipment and systems, it is urgent to further research the theory and method of uncertain design optimization. In the future, uncertain design optimization is one of the important research directions in the field of advanced manufacturing.When analyzing the uncertain factors of large and complex equipment and systems, it is difficult to obtain the real distribution of uncertain variables because of the large amount of data, but it is relatively convenient to obtain the interval bounds of uncertain variables, and the amount of data required is relatively small [5]. Many researchers are looking for effective methods to reduce the computational burden and improve computational accuracy [6–8], while the interval number method has been proposed in the field of mathematics. Interval design optimization is an uncertainty optimization method that does not need to calculate the probability distribution function of uncertain design variables, which uses interval range to characterize the uncertainty range of variables. It requires a relatively low number and accuracy of data, and the calculation is relatively simple, which leads to its obvious advantages in engineering. Wu et al. [9] used the Taylor expansion method to build a polynomial approximation model and combined the interval algorithm to solve the uncertainty problem of the vehicle suspension. Liu et al. [10] proposed a novel interval uncertainty formulation for exploring the impact of epistemic uncertainty on reliability-constrained design performance in a composite stiffened panel and a cantilever tube. Wang et al. [11] proposed a dimension-by-dimension interval uncertainty modeling and analysis method, which was solved by volume ratio theory.Approximate models can improve the design optimization efficiency of complex systems, which lays the foundation for efficient design of complex systems [12, 13]. Rafiee and Faiz [14] used the Box–Behnken response surface method and particle swarm optimization (PSO) algorithm to establish an approximate model of a high-dimensional system. Qian et al. [15] proposed a general sequential constraint update method based on approximate model confidence intervals, which introduced target switching and sequential updating strategies based on the interpolation uncertainty problem of the Kriging model. Meng et al. [16] solved the multidisciplinary problem with uncertainty based on the CO method and approximate model, which improved the computational efficiency of design optimization. Meanwhile, a six-sigma method is used to design a permanent magnet motor for hybrid vehicles. Lieu et al. [17] proposed an adaptive surrogate model using neural networks for reliability analysis by randomly selecting initial experimental design points from a given population of Monte Carlo simulations and to build an approximate model of the performance function, which improved the accuracy of the surrogate model.Design optimization considering uncertainty is one of the factors affecting the whole equipment performance. In addition, modern engineering systems usually involve overlapping and fusion of multidiscipline, which often consists of multiple subsystems and components. So, the design method of multidisciplinary integration will be the trend of future research [18, 19]. In this context, multidisciplinary design optimization (MDO) theory has attracted the attention of researchers, which is a classical method for the design optimization of complex system.The principle is to deal with the coupling relationship between disciplines through some decoupling methods [20–23]. At the same time, to optimize the performance of the entire system through a single-layer or multilayer optimization strategy, the multidisciplinary coupling problem is solved hierarchically [24–26]. Since the uncertainty in the subsystem is propagated in the multidisciplinary coupling, to quantify the impact of uncertain factors accurately and take advantage of MDO simultaneously, the RBMDO has received extensive attention and has become a research hotspot [27–29]. Zaman and Mahadevan [30] proposed an algorithm for reliability-based design optimization; meanwhile, they have proposed a specialized formulation to deal with the epistemic uncertainty caused by interval data. Wang et al. [31] developed a reliability analysis method based on sequence optimization for multidisciplinary systems under mixed uncertainty, which decoupled reliability analysis and solution of deterministic MDO problems.In the design optimization process of engineering equipment and system, it is often encountered that the double-nested optimization model is difficult to solve due to the high complexity of the system. According to the numerical analysis method to solve the double-nested interval optimization problem, there are two problems as follows. (1) In general, the inner layer conducts uncertainty analysis while the outer layer conducts design variable solving. The outer layer will call the inner layer, solving cycle dozens or even hundreds of times at a time, resulting in an exponential increase in computation, which often takes a long time and affects efficiency. (2) In practical engineering, it is often difficult to obtain the real optimization model. Aiming at the above problems, this study proposes an interval optimization method of the point-filled Kriging approximation model. On the basis of the improved PSO, the double-nested optimization model is simplified to improve efficiency. At the same time, the accuracy of Kriging approximation model is improved by adding the distribution of test points to meet the engineering requirements.In the process of design optimization of complex equipment and system, it is inevitable that there will be coupling between different disciplines and the uncertainty between different disciplines will also be coupled, which makes it a challenge to solve the uncertainty optimization of complex systems involving multiple disciplines. So, this study proposes an optimization method of interval management design based on the Kriging model. In this model, the RBMDO problem considering interval uncertainty is transformed into a deterministic MDO problem, which is analyzed by a collaborative optimization decoupling framework.The content of this study is as follows. Section2 includes a double-nested uncertainty optimization model and its solution strategy. A Kriging model based on the point-infilled criterion is constructed in Section 3. Section 4 presents an RBMDO method problem considering interval uncertainty based on the point-infilled Kriging model. A mathematical example and an engineering example are used to verify the efficiency and accuracy of the proposed method in Section 5.


## 2. Double-Nested Optimization Strategy under Interval Uncertainty

The constructions of interval uncertainty model and interval optimization model analysis have been completed in previous parts. The uncertainty optimization problem has been transformed into a nested optimization problem. The inner layer deals with uncertain factors and the outer layer solves the optimization model. However, it is difficult to solve nested problems directly; thus, the algorithm requires certain adaptability and convergence. In this study, a two-layer nested optimization strategy based on improved PSO is used to solve the problem. The inner layer analyzes the uncertainty of the model and transforms the uncertainty model. The outer layer uses the evaluation function method to construct a new optimization model to solve according to the results from the inner layer.


### 2.1. Interval Optimization Model Based on Evaluation Function

For a multidisciplinary optimization model, it is not satisfying to find an effective solution or a weakly effective solution in the calculation process. Generally, a decision-making method is introduced to make decisions for the solution to obtain the optimal solution that can satisfy the engineering application. In the process of multidisciplinary optimization, an evaluation function is constructed to convert the multidisciplinary optimization problem into a single-objective optimization problem. As shown in (1), it is a multidisciplinary optimization model:(1)minDVFX=minf1x,f2x,…,fnxs.t.hX=0gX≥0DV=X,where f• represents the objective function, h• express equality constraints, g• depicts an inequality constraint, DV is the design variable space, and X is design variables. In this study, f•, h•, and g• are defaulted to be real-valued continuous functions that can obtain second-order continuous partial derivatives.The optimization model constructed by using the evaluation function for the above model is shown in the following equation:(2)minDVμFX=minμf1x,f2x,…,fnxs.t.hX=0gX≥0DV=X.Evaluation function methods that are commonly used include ideal point method, maximum and minimum value method, and linear weighting method [32]. This study will use the ideal point method to calculate and solve the transformed multi-objective function. The principle of the ideal point method is to find the minimum value of each objective function, then comprehensively consider all objective functions. Finally, we make each objective function as close to its respective minimum value as possible. Compared with the linear weighting method and the maximum and minimum method, the ideal point method has great advantages in calculation accuracy. The weight can be set according to the requirements of the parameters in the actual project.For each subobjective function of the objective functionFX=minf1x,f2x,…,fnx, the optimal solution after minimization is shown in the following equation:(3)fi∗x=minfixi=1,2,…,n.Assuming thatxii=1,2,…,n corresponding to the optimal solution of each objective function is the same, then xi is the optimal solution of FX. In general, the optimal values of each objective function do not correspond to exactly the same. The minimum value of each objective function is shown in (4). The minimum value point of each objective function is called the ideal point of the model:(4)F∗X=f1∗x,f2∗x,…,fn∗x.We introduce another point in the target spaceRn. We take the distance between the introduced point f and the ideal point f∗ as the evaluation function. We find the value as close as possible to the ideal point within the feasible range X. The evaluation function can be constructed as follows:(5)minx∈XμF=minμf−f∗2=min∑i=1nμifi−fi∗2,where μi indicates the importance factor of the ith objective function. For the convenience of calculation, it can be converted to (6) for solving in the actual calculation process, which improves the solving efficiency:(6)minx∈XμF=minx∈X∑i=1nμifi−fi∗,where∑i=1nμi=1 for convenience of calculation.The calculation steps of this method are as follows:Step 1: we seek the ideal point. We find the minimum value of each objective function,fi∗x=minfix, i=1,2,…,n. The ideal point is shown in the following equation:(7)F∗X=f1∗x,f2∗x,…,fn∗x.Step 2: we check the ideal point. If the ideal points are all equal,x1=x2=⋯=xn, the optimal solution x∗=xi is output. If x1,x2,…,xn is not exactly the same, the minimization distance is solved by (6). The optimal solution x∗ is obtained by solving the minimization problem.Through a series of model processing, the optimization model can be expressed as follows:(8)minDVfMX=minμfLX−fL∗+1−μfRX−fR∗s.t.PhIX=aI≥λPgIX≥bI≥θDV=X,X∈Rn.For the converted optimization model of (8), the optimization model is processed by the penalty function method. The optimization model of the above equation can be transformed into the following:(9)minDVf=fMX+σPXPX=∑i=1mΘPgIX≥bI−θ+∑i=1nΓPhIX=aI−λDV=X,X∈Rn,where σ is the penalty factor, PX is the penalty function, and Θ and Γ are the continuous function. The typical method is as follows:(10)Θ=max0,−PgIX≥bI−θαΓ=max0,−PhIX=aI−λβ,where α≥1andβ≥1; it can be set α=β=2 generally.


### 2.2. Double-Nested Optimization Strategy Based on Improved PSO

The above section has completed the construction of interval uncertainty model and the analysis of interval optimization model. For such models, people usually transform the uncertain optimization problem into a nested optimization problem to solve it. The inner layer solves the uncertain factors, and the outer layer solves the optimization model. It is difficult to solve the nested problem directly, and the algorithm is required to have certain adaptability and convergence. PSO algorithm is a widely used optimization algorithm. This section mainly solves the optimization model through a double-layer nested optimization solution strategy based on improved PSO. The solution strategy under the intelligent optimization algorithm can improve the solution efficiency and avoid falling into the local optimal solution in a certain range.The difference between PSO and other algorithms is that there is no evolutionary operator. In anN-dimensional space, there are m particle swarms, where the ith particle is represented as an N-dimensional vector x⟶i=xi1,xi2,…,xiN, i=1,2,…,m. The position of each particle is the current solution. According to the x⟶i position, the fitness value of the objective function is obtained, and then, the quality of the solution is judged.The velocity of theith particle is recorded as V⟶=Vi1,Vi2,…ViN. The current optimal position of the ith particle is recorded as P⟶iN=Pi1,Pi2,…,PiN. The current optimal position of the entire particle swarm is P⟶O=PO1,PO2,…,PON.Assumingfx is the objective function to be minimized, the current best position of ith particle is shown in the following equation:(11)pit+1=pit⟶fxit+1≥fpitXit+1⟶fxit+1<fpit.The evolution running equation for particles is defined as follows:(12)vint+1=vint+c1PSOr1ntpint−xint+c2PSOr2ntpont−xintxint+1=xint+vint+1,where i=1,m,n=1,N, c1PSO and c2PSO represent the step size, r1PSO and r2PSO are random numbers in the range of 0,1, and viN∈−vmax,vmax and vmax are the maximum speed. If the search space is in −vmax,vmax, it can be set as vmax=kxmax0.1≤k≤1.0.The particle velocity is improved on the basis of (12), which is shown in the following equation:(13)vint+1=ωspeed⋅vint+c1r1npint−xint+c2r2ntpont−xontxint+1=xint+vint+1,where ωspeed is the inertia weight, which is used to control the particle speed. When ωspeed is large, it is suitable for global search optimization. When ωspeed is small, it is suitable for local search optimization. We change its particle search ability by resizing ωspeed. To better search for the global optimal solution, ωspeed = 1.1 in this study.According to [33], if the system is not sensitive with the overall reliability, the constraint possibility value can be small; otherwise, a larger constraint probability value needs to be set. The setting of the objective function importance factor depends on the needs of actual engineering problems. In the subsequent examples of this study, appropriate parameter values are selected according to different situations.In this study, both the inner and outer optimization algorithms use the maximum number of iterations as the termination condition. The solution flowchart is shown in Figure1.Figure 1
Double-nested optimization strategy.When solving the interval optimization problem according to the numerical analysis method, there are two levels of nesting problems: the inner layer is called upon to solve the loop dozens or even hundreds of times for one cycle of the outer layer and computational costs would increase exponentially, which also affects the calculation efficiency tremendously. Therefore, the Kriging model based on the point-infilled criterion is proposed to overcome the efficiency problem.


## 2.1. Interval Optimization Model Based on Evaluation Function

For a multidisciplinary optimization model, it is not satisfying to find an effective solution or a weakly effective solution in the calculation process. Generally, a decision-making method is introduced to make decisions for the solution to obtain the optimal solution that can satisfy the engineering application. In the process of multidisciplinary optimization, an evaluation function is constructed to convert the multidisciplinary optimization problem into a single-objective optimization problem. As shown in (1), it is a multidisciplinary optimization model:(1)minDVFX=minf1x,f2x,…,fnxs.t.hX=0gX≥0DV=X,where f• represents the objective function, h• express equality constraints, g• depicts an inequality constraint, DV is the design variable space, and X is design variables. In this study, f•, h•, and g• are defaulted to be real-valued continuous functions that can obtain second-order continuous partial derivatives.The optimization model constructed by using the evaluation function for the above model is shown in the following equation:(2)minDVμFX=minμf1x,f2x,…,fnxs.t.hX=0gX≥0DV=X.Evaluation function methods that are commonly used include ideal point method, maximum and minimum value method, and linear weighting method [32]. This study will use the ideal point method to calculate and solve the transformed multi-objective function. The principle of the ideal point method is to find the minimum value of each objective function, then comprehensively consider all objective functions. Finally, we make each objective function as close to its respective minimum value as possible. Compared with the linear weighting method and the maximum and minimum method, the ideal point method has great advantages in calculation accuracy. The weight can be set according to the requirements of the parameters in the actual project.For each subobjective function of the objective functionFX=minf1x,f2x,…,fnx, the optimal solution after minimization is shown in the following equation:(3)fi∗x=minfixi=1,2,…,n.Assuming thatxii=1,2,…,n corresponding to the optimal solution of each objective function is the same, then xi is the optimal solution of FX. In general, the optimal values of each objective function do not correspond to exactly the same. The minimum value of each objective function is shown in (4). The minimum value point of each objective function is called the ideal point of the model:(4)F∗X=f1∗x,f2∗x,…,fn∗x.We introduce another point in the target spaceRn. We take the distance between the introduced point f and the ideal point f∗ as the evaluation function. We find the value as close as possible to the ideal point within the feasible range X. The evaluation function can be constructed as follows:(5)minx∈XμF=minμf−f∗2=min∑i=1nμifi−fi∗2,where μi indicates the importance factor of the ith objective function. For the convenience of calculation, it can be converted to (6) for solving in the actual calculation process, which improves the solving efficiency:(6)minx∈XμF=minx∈X∑i=1nμifi−fi∗,where∑i=1nμi=1 for convenience of calculation.The calculation steps of this method are as follows:Step 1: we seek the ideal point. We find the minimum value of each objective function,fi∗x=minfix, i=1,2,…,n. The ideal point is shown in the following equation:(7)F∗X=f1∗x,f2∗x,…,fn∗x.Step 2: we check the ideal point. If the ideal points are all equal,x1=x2=⋯=xn, the optimal solution x∗=xi is output. If x1,x2,…,xn is not exactly the same, the minimization distance is solved by (6). The optimal solution x∗ is obtained by solving the minimization problem.Through a series of model processing, the optimization model can be expressed as follows:(8)minDVfMX=minμfLX−fL∗+1−μfRX−fR∗s.t.PhIX=aI≥λPgIX≥bI≥θDV=X,X∈Rn.For the converted optimization model of (8), the optimization model is processed by the penalty function method. The optimization model of the above equation can be transformed into the following:(9)minDVf=fMX+σPXPX=∑i=1mΘPgIX≥bI−θ+∑i=1nΓPhIX=aI−λDV=X,X∈Rn,where σ is the penalty factor, PX is the penalty function, and Θ and Γ are the continuous function. The typical method is as follows:(10)Θ=max0,−PgIX≥bI−θαΓ=max0,−PhIX=aI−λβ,where α≥1andβ≥1; it can be set α=β=2 generally.


## 2.2. Double-Nested Optimization Strategy Based on Improved PSO

The above section has completed the construction of interval uncertainty model and the analysis of interval optimization model. For such models, people usually transform the uncertain optimization problem into a nested optimization problem to solve it. The inner layer solves the uncertain factors, and the outer layer solves the optimization model. It is difficult to solve the nested problem directly, and the algorithm is required to have certain adaptability and convergence. PSO algorithm is a widely used optimization algorithm. This section mainly solves the optimization model through a double-layer nested optimization solution strategy based on improved PSO. The solution strategy under the intelligent optimization algorithm can improve the solution efficiency and avoid falling into the local optimal solution in a certain range.The difference between PSO and other algorithms is that there is no evolutionary operator. In anN-dimensional space, there are m particle swarms, where the ith particle is represented as an N-dimensional vector x⟶i=xi1,xi2,…,xiN, i=1,2,…,m. The position of each particle is the current solution. According to the x⟶i position, the fitness value of the objective function is obtained, and then, the quality of the solution is judged.The velocity of theith particle is recorded as V⟶=Vi1,Vi2,…ViN. The current optimal position of the ith particle is recorded as P⟶iN=Pi1,Pi2,…,PiN. The current optimal position of the entire particle swarm is P⟶O=PO1,PO2,…,PON.Assumingfx is the objective function to be minimized, the current best position of ith particle is shown in the following equation:(11)pit+1=pit⟶fxit+1≥fpitXit+1⟶fxit+1<fpit.The evolution running equation for particles is defined as follows:(12)vint+1=vint+c1PSOr1ntpint−xint+c2PSOr2ntpont−xintxint+1=xint+vint+1,where i=1,m,n=1,N, c1PSO and c2PSO represent the step size, r1PSO and r2PSO are random numbers in the range of 0,1, and viN∈−vmax,vmax and vmax are the maximum speed. If the search space is in −vmax,vmax, it can be set as vmax=kxmax0.1≤k≤1.0.The particle velocity is improved on the basis of (12), which is shown in the following equation:(13)vint+1=ωspeed⋅vint+c1r1npint−xint+c2r2ntpont−xontxint+1=xint+vint+1,where ωspeed is the inertia weight, which is used to control the particle speed. When ωspeed is large, it is suitable for global search optimization. When ωspeed is small, it is suitable for local search optimization. We change its particle search ability by resizing ωspeed. To better search for the global optimal solution, ωspeed = 1.1 in this study.According to [33], if the system is not sensitive with the overall reliability, the constraint possibility value can be small; otherwise, a larger constraint probability value needs to be set. The setting of the objective function importance factor depends on the needs of actual engineering problems. In the subsequent examples of this study, appropriate parameter values are selected according to different situations.In this study, both the inner and outer optimization algorithms use the maximum number of iterations as the termination condition. The solution flowchart is shown in Figure1.Figure 1
Double-nested optimization strategy.When solving the interval optimization problem according to the numerical analysis method, there are two levels of nesting problems: the inner layer is called upon to solve the loop dozens or even hundreds of times for one cycle of the outer layer and computational costs would increase exponentially, which also affects the calculation efficiency tremendously. Therefore, the Kriging model based on the point-infilled criterion is proposed to overcome the efficiency problem.


## 3. Interval Optimization Method Based on the Point Strategy of the Kriging Model

### 3.1. Point-Infilled Kriging Model

The first step in building an approximate model is to conduct a design of experiment for collecting sample points. The method and uniformity of sample point’s collection will affect the construction accuracy directly of the approximate model. Since LHS method has a good balance between high sampling efficiency and convenience, which is widely used to construct approximate models in practical engineering [34], the LHS method is selected for sampling in this study.The Kriging model consists of an approximate regression and a nonparametric part [35]:(14)f^x=∑i=1nβiqix+Zx,where x represents the sample point, f^• represents the approximate objective function, βi represents the polynomial coefficients of the approximate objective function, and qi• represents the polynomial basis function in the polynomial function. The combination of polynomial basis functions forms the basis function vector, which can be expressed as qx=q1x,q2x,…,qpxT. p represents the number of regression parameters, related to the number of design variables. Z• represents a random process that provides an approximation of the local deviation of the analyzed object.In the optimization based on the approximate model, the selection of update points will affect the calculation efficiency and the accuracy of optimization results directly. LHS is used as the sample collection method for the constructed Kriging model, which has stable performance and high fitting accuracy. The Kriging model has the predictive ability, which can improve computational efficiency to a certain extent [36]. However, the Kriging model has sparse sample points and insufficient model construction accuracy near the optimal solution. It is not easy to find the optimal solution that satisfies the accuracy requirements in practical engineering applications. From there, many scholars use local exploration and global search to improve the function fitting accuracy of approximate models. The reliable methods are the maximum probability of improvement (PI) criterion and the maximum improvement expectation criterion (EI criterion). In this study, the distribution of sample points is strengthened by the EI function. The maximum value is obtained by analyzing the probability density of the existing samples at the extreme value of the current predicted value. The EI function is defined as follows [37]:(15)EIx=fmin−f^xΦfmin−f^xs^2x+sϕfmin−f^xs^2xs>00s=0,where Φ• represents the normal distribution function, ϕ• represents the probability density function, fmin represents the minimum value in the current sample point, and f^• and s^• are the predicted value and standard deviation. By calculating the maximum value of the EI function and infill sampling criterion, the distance between different sample points in the Kriging model can be narrowed and the sample points can be distributed adaptively. This section uses the EI function to predict and infill points to the initially constructed Kriging model, enhance the distribution of sample points in the feasible region, enhance the fitting accuracy, and solve the global optimal solution accurately. The model construction and calculation principle flow of this method is shown in Figure 2. The implementation steps are as follows:Step 1: in a given design space, the LHS method is used for uniform sampling to obtain the initial test sample point set and the initial sample library is established.Step 2: we analyze the system model. We train the sample point set through the Kriging approximation model method of the initial sample library established in the first step. We construct the Kriging approximation modelf^x of the optimized model fx.Step 3: we use the improved PSO as a solver to obtain the current minimum point off^x. We obtain the maximum extreme value point of the EI function by solving the EI function near the obtained current minimum point.Step 4: we use the currently calculated maximum extremum point of the EI function as a newly infill sampling point; meanwhile, we update the test sample point set.Step 5: we use the updated sample point set as a new test sample to correct the established approximate model. In this way, the fitting accuracy of the approximate model is improved until the newly infill sampling points change within a range that meets the required accuracy. We stop infilling points to get an approximate model with the required accuracy.Figure 2
Flowchart of the Kriging approximation model for the infilling criterion.


### 3.2. Model Test

In this section, the point-infilled Kriging approximation model strategy is verified by testing the function. The mathematical model of the test function is shown as follows:(16)minfx1,x2=0.01x1−x222+1−x2+24−x12+7sin0.5x2sin0.5x1x2s.t.2≤x1≤7,0≤x2≤5,where f• is the objective function and x1 and x2 are the design variables in the test example. The space of the design variables is normalized. The true response surface of the test function is shown in Figure 3, and the corresponding contour of the test function is shown in Figure 4. Figures 3 and 4 show that the test function has a strong degree of nonlinearity, and it can better reflect the effectiveness of the established approximate model for testing.Figure 3
The true response surface of the test function.Figure 4
The true contour of the test function.The fitting verification of the test function shown in (16) is carried out by using the approximate model based on the point-infilled strategy. First, Latin hypercube sampling is used for sampling. Twenty and 30 sampling points are selected, and the sampling results are shown in Figure 5. It shows that the sampling points are uniform.Figure 5
Schematic diagram of LHS sampling. (a) 20 sampling points (b) 30 sampling points.
(a)(b)The Kriging approximation model is used to approximate the fitting test function. The corresponding contour of 20 and 30 sampling points with fitting are shown in Figure6. Figure 6 shows the accuracy of the approximate model varies with the number of initial sample points. The more the initial sample points, the higher the accuracy of the constructed approximate model, which is more similar to the real model.Figure 6
Kriging approximation model. (a) 20 sampling points. (b) 30 sampling points.
(a)(b)On the basis of the Kriging approximation model, the EI function point-infilled strategy is used to modify the approximation model, and the region range of the optimal point is quantified. The number of points of EI function is set to 10, and the approximate model is shown in Figure7. It can be seen from the figure that the EI function has a good effect on the revision of the approximate model, and the built model is more similar to the real model. In addition, the infilled points are concentrated is the area where the best points are located.Figure 7
Kriging approximation model with EI point-infilled strategy. (a) 20 sampling points. (b) 30 sampling points.
(a)(b)In order to more clearly express the process of constructing Kriging approximate model with EI point-infilled strategy, the initial sampling of 30 test points of LHS is shown in Figure8.Figure 8
Kriging approximate model adding points process: (a) add 1 test point, (b) add 3 test points, (c) add 5 test points, (d) add 7 test points, (e) add 9 test points, and (f) add 10 test points.
(a)(b)(c)(d)(e)(f)Compared with the Kriging approximation model, the point-infilled Kriging approximation model has a higher fitting degree to the test function and is closer to the actual model. In addition, the point range of EI function is consistent with the range of the optimal solution, so it can be better applied in practical projects.


### 3.3. Interval Optimization Method Based on the Kriging Model with Point Strategy

In complex engineering problems considering interval uncertainty, the approximate model method is an important method to solve the complex optimization model. Interval uncertainty optimization problem can be transformed into approximate interval uncertainty optimization problem as shown below:(17)minDVf^X,XUIs.t.Ph^X,XUI=aI≥λPg^X,XUI≥bI≥θDV=X,XUI,XUI∈XUL,XUR,X∈Rn,where f^• represents the approximate objective function, h^• represents the approximate model with equality constraints, g^• represents the approximate model with inequality constraints, and P• represents the constraint possibility degree function. By applying the evaluation function method and penalty function method in Section 2.1, it can be transformed into the approximate model as shown below:(18)minDVf^=f^MX+σP^XP^X=∑i=1mϕPg^IX≥bI−θ+∑i=1nφPh^IX=aI−λDV=X,X∈Rn,where f^ represents the approximate model of penalty function under the approximate model of objective function and constraint conditions and f^M• represents the multiobjective evaluation function.The optimization method process based on the point-infilled Kriging approximation model considering interval uncertainty is shown in Figure9.Figure 9
Flowchart of interval optimization method based on point-infilled Kriging model strategy.The Kriging model based on the point-infilled criterion can replace the real model to solve the problem in the case of an unknown optimization model. This approximation model considering interval uncertainty solves the efficiency problem caused by double-nested optimization of interval uncertainty problem. Meanwhile, the approximation model also has high accuracy.


## 3.1. Point-Infilled Kriging Model

The first step in building an approximate model is to conduct a design of experiment for collecting sample points. The method and uniformity of sample point’s collection will affect the construction accuracy directly of the approximate model. Since LHS method has a good balance between high sampling efficiency and convenience, which is widely used to construct approximate models in practical engineering [34], the LHS method is selected for sampling in this study.The Kriging model consists of an approximate regression and a nonparametric part [35]:(14)f^x=∑i=1nβiqix+Zx,where x represents the sample point, f^• represents the approximate objective function, βi represents the polynomial coefficients of the approximate objective function, and qi• represents the polynomial basis function in the polynomial function. The combination of polynomial basis functions forms the basis function vector, which can be expressed as qx=q1x,q2x,…,qpxT. p represents the number of regression parameters, related to the number of design variables. Z• represents a random process that provides an approximation of the local deviation of the analyzed object.In the optimization based on the approximate model, the selection of update points will affect the calculation efficiency and the accuracy of optimization results directly. LHS is used as the sample collection method for the constructed Kriging model, which has stable performance and high fitting accuracy. The Kriging model has the predictive ability, which can improve computational efficiency to a certain extent [36]. However, the Kriging model has sparse sample points and insufficient model construction accuracy near the optimal solution. It is not easy to find the optimal solution that satisfies the accuracy requirements in practical engineering applications. From there, many scholars use local exploration and global search to improve the function fitting accuracy of approximate models. The reliable methods are the maximum probability of improvement (PI) criterion and the maximum improvement expectation criterion (EI criterion). In this study, the distribution of sample points is strengthened by the EI function. The maximum value is obtained by analyzing the probability density of the existing samples at the extreme value of the current predicted value. The EI function is defined as follows [37]:(15)EIx=fmin−f^xΦfmin−f^xs^2x+sϕfmin−f^xs^2xs>00s=0,where Φ• represents the normal distribution function, ϕ• represents the probability density function, fmin represents the minimum value in the current sample point, and f^• and s^• are the predicted value and standard deviation. By calculating the maximum value of the EI function and infill sampling criterion, the distance between different sample points in the Kriging model can be narrowed and the sample points can be distributed adaptively. This section uses the EI function to predict and infill points to the initially constructed Kriging model, enhance the distribution of sample points in the feasible region, enhance the fitting accuracy, and solve the global optimal solution accurately. The model construction and calculation principle flow of this method is shown in Figure 2. The implementation steps are as follows:Step 1: in a given design space, the LHS method is used for uniform sampling to obtain the initial test sample point set and the initial sample library is established.Step 2: we analyze the system model. We train the sample point set through the Kriging approximation model method of the initial sample library established in the first step. We construct the Kriging approximation modelf^x of the optimized model fx.Step 3: we use the improved PSO as a solver to obtain the current minimum point off^x. We obtain the maximum extreme value point of the EI function by solving the EI function near the obtained current minimum point.Step 4: we use the currently calculated maximum extremum point of the EI function as a newly infill sampling point; meanwhile, we update the test sample point set.Step 5: we use the updated sample point set as a new test sample to correct the established approximate model. In this way, the fitting accuracy of the approximate model is improved until the newly infill sampling points change within a range that meets the required accuracy. We stop infilling points to get an approximate model with the required accuracy.Figure 2
Flowchart of the Kriging approximation model for the infilling criterion.


## 3.2. Model Test

In this section, the point-infilled Kriging approximation model strategy is verified by testing the function. The mathematical model of the test function is shown as follows:(16)minfx1,x2=0.01x1−x222+1−x2+24−x12+7sin0.5x2sin0.5x1x2s.t.2≤x1≤7,0≤x2≤5,where f• is the objective function and x1 and x2 are the design variables in the test example. The space of the design variables is normalized. The true response surface of the test function is shown in Figure 3, and the corresponding contour of the test function is shown in Figure 4. Figures 3 and 4 show that the test function has a strong degree of nonlinearity, and it can better reflect the effectiveness of the established approximate model for testing.Figure 3
The true response surface of the test function.Figure 4
The true contour of the test function.The fitting verification of the test function shown in (16) is carried out by using the approximate model based on the point-infilled strategy. First, Latin hypercube sampling is used for sampling. Twenty and 30 sampling points are selected, and the sampling results are shown in Figure 5. It shows that the sampling points are uniform.Figure 5
Schematic diagram of LHS sampling. (a) 20 sampling points (b) 30 sampling points.
(a)(b)The Kriging approximation model is used to approximate the fitting test function. The corresponding contour of 20 and 30 sampling points with fitting are shown in Figure6. Figure 6 shows the accuracy of the approximate model varies with the number of initial sample points. The more the initial sample points, the higher the accuracy of the constructed approximate model, which is more similar to the real model.Figure 6
Kriging approximation model. (a) 20 sampling points. (b) 30 sampling points.
(a)(b)On the basis of the Kriging approximation model, the EI function point-infilled strategy is used to modify the approximation model, and the region range of the optimal point is quantified. The number of points of EI function is set to 10, and the approximate model is shown in Figure7. It can be seen from the figure that the EI function has a good effect on the revision of the approximate model, and the built model is more similar to the real model. In addition, the infilled points are concentrated is the area where the best points are located.Figure 7
Kriging approximation model with EI point-infilled strategy. (a) 20 sampling points. (b) 30 sampling points.
(a)(b)In order to more clearly express the process of constructing Kriging approximate model with EI point-infilled strategy, the initial sampling of 30 test points of LHS is shown in Figure8.Figure 8
Kriging approximate model adding points process: (a) add 1 test point, (b) add 3 test points, (c) add 5 test points, (d) add 7 test points, (e) add 9 test points, and (f) add 10 test points.
(a)(b)(c)(d)(e)(f)Compared with the Kriging approximation model, the point-infilled Kriging approximation model has a higher fitting degree to the test function and is closer to the actual model. In addition, the point range of EI function is consistent with the range of the optimal solution, so it can be better applied in practical projects.


## 3.3. Interval Optimization Method Based on the Kriging Model with Point Strategy

In complex engineering problems considering interval uncertainty, the approximate model method is an important method to solve the complex optimization model. Interval uncertainty optimization problem can be transformed into approximate interval uncertainty optimization problem as shown below:(17)minDVf^X,XUIs.t.Ph^X,XUI=aI≥λPg^X,XUI≥bI≥θDV=X,XUI,XUI∈XUL,XUR,X∈Rn,where f^• represents the approximate objective function, h^• represents the approximate model with equality constraints, g^• represents the approximate model with inequality constraints, and P• represents the constraint possibility degree function. By applying the evaluation function method and penalty function method in Section 2.1, it can be transformed into the approximate model as shown below:(18)minDVf^=f^MX+σP^XP^X=∑i=1mϕPg^IX≥bI−θ+∑i=1nφPh^IX=aI−λDV=X,X∈Rn,where f^ represents the approximate model of penalty function under the approximate model of objective function and constraint conditions and f^M• represents the multiobjective evaluation function.The optimization method process based on the point-infilled Kriging approximation model considering interval uncertainty is shown in Figure9.Figure 9
Flowchart of interval optimization method based on point-infilled Kriging model strategy.The Kriging model based on the point-infilled criterion can replace the real model to solve the problem in the case of an unknown optimization model. This approximation model considering interval uncertainty solves the efficiency problem caused by double-nested optimization of interval uncertainty problem. Meanwhile, the approximation model also has high accuracy.


## 4. The Optimization Method of Interval Multidisciplinary Design Based on the Kriging Model

### 4.1. RBMDO Method considering Interval Uncertainty

In this study, RBMDO turns uncertainty into reliability constraints, which are added to the MDO process. In the previous sections, the RBMDO problem considering interval uncertainty has been transformed into deterministic MDO problem. Based on the traditional deterministic MDO strategy, this section studies multidisciplinary problems considering interval uncertainty.The MDO process considering interval uncertainty is similar to the traditional single-discipline interval optimization process. The difference is that the solution process of the interval uncertainty multidisciplinary optimization model is three-layer nested iterative optimization. The inner layer needs to decouple multidisciplinary problems, the middle layer analyzes the decoupled interval uncertainty model, and the outer layer solves the design variables through intelligent optimization algorithms. The multidisciplinary optimization model under interval uncertainty can be expressed as(19)FindXs,XminDVfXs,X,XsU,XU,Ys.t.hXs,X,XsU,XU,Y=aIgXs,X,XsU,XU,Y≤bIDV=Xs,X,XsU,XUXsU∈XsUL,XsUR,XU∈XUL,XUR,where XsU represents the interval uncertainty shared design variable, XsUL and XsUR are the lower and upper bounds, respectively, XU represent the interval uncertainty local design variables, XUL and XUR are the lower and upper bounds, respectively, aI represents the equality constraint interval, and bI represents the inequality constraint interval.The CO method is used to decouple the optimization model of (19). The system layer optimization model can be expressed as(20)System:FindZs,ZjminDVf=fZs,Zj,ZsU,ZUj,Y^s.t.JiZs,Zj,ZsU,ZUj,Y^≤εDV=Zs,Zj,ZsU,ZUj,Y^,where Y^ represents the coupling design variable, ZsU represents the uncertainty shared design variable at the system layer, ZUj represents the jth uncertainty design variable at the system layer, and Ji•=0 represents the discipline consistency constraint of discipline i, which can be expressed as(21)JiZs,Zij,Y^=Xsi−Zs22+Xij−Zj22+XsUi−ZsU22+XUij−ZUj22+Y•i−Y^•i22+Yi•Xsi,Xij,Y•i−Y^i•22,where Y•i represents the coupling variable to discipline i from other disciplines and Y^i·• represents the coupling variable input from discipline i. The discipline optimization model of the ith discipline in the discipline layer can be expressed as(22)Discipline−i:minDVJi=JiXsi,Xij,XsUi,XUij,Y•is.t.hiXsi,Xij,XsUi,XUij,Y•i=aiIgiXsi,Xij,XsUi,XUij,Y•i≤biIDV=Xsi,Xij,XsUi,XUij,Y•i,where XsUi represents the interval shared design variable of discipline i and XUij represents the jth interval local design variable of discipline i.


### 4.2. Interval Multidisciplinary Collaborative Optimization Method Based on the Kriging Model

In the previous section, the multidisciplinary collaborative optimization method with interval uncertainty was studied. Since the multidisciplinary systems usually have complex optimization models or difficulty in building optimization models. In this section, the approximate model technique is introduced to solve the MDO model under interval uncertainty. By approximating the system layer and discipline layer objective functions and constraints under the CO decoupling strategy, the difficulty of solving is reduced. The multidisciplinary collaborative optimization method framework under interval uncertainty based on the approximate model is shown in Figure10.Figure 10
Frame diagram of interval multidisciplinary collaborative optimization method based on the approximate model.


#### 4.2.1. Constructing the MDO Model under Interval Uncertainty Based on the Kriging Approximation Model

First of all, for the system layer optimization model, the Kriging model is fitted, which is shown as follows:(23)System:minDVf^=f^Zs,Zj,ZsU,ZUj,Y^s.t.JiZs,Zj,ZsU,ZUj,Y^≤εDV=Zs,Zj,ZsU,ZUj,Y^,where f^• represents the approximate model of the objective function. The fitting function can be expressed as(24)f^Zs,Zj,ZsU,ZUj,Y^=qZs,Zj,ZsU,ZUj,Y^Tβ^+rTZs,Zj,ZsU,ZUj,Y^R−1F−Qβ^,where qi• is the basis function vector, β^ is the regression parameter matrix obtained by the optimal linear unbiased estimation, r• is the correlation vector between the sample value and the predicted value, R• is the spatial correlation function matrix, F is the response function, and Q is the polynomial matrix. The parameter representation is shown below:(25)qZs,Zj,ZsU,ZUj,Y^=q1Zs,Zj,ZsU,ZUj,Y^⋯qpZs,Zj,ZsU,ZUj,Y^β^=QTR−1Q−1QTR−1FrZs,Zj,ZsU,ZUj,Y^=RZs,Zj,ZsU,ZUj,Y^;Zs1,Z1,ZsU1,ZU1,Y^1⋯RZs,Zj,ZsU,ZUj,Y^;Zsn,Zn,ZsUn,ZUn,Y^nR=RZs1,Z1,ZsU1,ZU1,Y^1;Zs1,Z1,ZsU1,ZU1,Y^1RZs1,Z1,ZsU1,ZU1,Y^1;Zsn,Zn,ZsUn,ZUn,Y^nRZsn,Zn,ZsUn,ZUn,Y^n;Zs1,Z1,ZsU1,ZU1,Y^1RZsn,Zn,ZsUn,ZUn,Y^n;Zsn,Zn,ZsUn,ZUn,Y^nF=f^1,f^2,…f^nTQ=qTZs1,Z1,ZsU1,ZU1,Y^1qTZsn,Zn,ZsUn,ZUnT.After the Kriging model of the system layer has been constructed, new sample points are infilled to the approximate model constructed by the system layer through the EI function criterion during the optimization and solution process. We rebuild the Kriging model and enhance its fitting accuracy until the solution optimization criteria are satisfied. The approximate model of the system layer is constructed.In the process of multidisciplinary optimization, the discipline layer conducts generally discipline analysis on the constraints of the discipline itself. For the discipline layer, the Kriging approximation model is used to fit its constraints under the CO strategy. The fitted optimization model is shown as follows:(26)Discipline−i:minDVJi=JiXsi,Xij,XsUi,XUij,Y•is.t.h^iXsi,Xij,XsUi,XUij,Y•i=aiIg^iXsi,Xij,XsUi,XUij,Y•i≤biIDV=Xsi,Xij,XsUi,XUij,Y•i,where h^i• represents the equality constraints of approximate fitting of discipline i and g^i• represents the inequality constraints of approximate fitting of discipline i. The fitting process can be expressed as follows:(27)h^iXsi,Xij,XsUi,XUij,Y•i=qXsi,Xij,XsUi,XUij,Y•iTλ^+rTXsi,Xij,XsUi,XUij,Y•iR−1H−Dλ^,where qi• is the basis function vector, λ^ the regression parameter matrix obtained by the optimal linear unbiased estimation, r• is the correlation vector, R• is the spatial correlation function matrix, H is the response function, and D is the polynomial matrix.(28)g^Xsi,Xij,XsUi,XUij,Y•i=qXsi,Xij,XsUi,XUij,Y•iTω^+rTXsi,Xij,XsUi,XUij,Y•iR−1G−Pω^,where qi• is the basis function vector, ω^ is the regression parameter matrix obtained by the optimal linear unbiased estimation, r• represents the correlation vector between the sample value and the predicted value,R• represents the spatial correlation function matrix, G is the response function, and P represents the polynomial matrix.After the preliminary approximate model is constructed for the constraints at the discipline layer, new sample points are infilled to the constructed approximate model through the EI function criterion. The Kriging approximate model is reconstructed, which can enhance the construction accuracy of the approximate model. The MDO problem with the interval uncertainty optimization approximation model under the CO framework is completed.


#### 4.2.2. MDO Problem considering Interval Uncertainty Solution Method Based on the Kriging Model

Through the analysis of the interval uncertainty multidisciplinary problem, the solution process can be divided into three layers of nested solution: multidisciplinary problem decoupling, interval uncertainty analysis, and objective function solution. However, due to the complexity of the multidisciplinary problem model, the use of approximate model technology to approximate the optimization model can reduce the computational difficulty in the solution process. It can also solve the problem that the model cannot converge due to the complexity of the model. From the above analysis, this section solves the RBMDO problem considering interval uncertainty based on the point-infilled Kriging model. The solution idea is shown in Figure11 and the solution steps are as follows:Step 1: we establish an interval uncertainty MDO optimization model. The hierarchical processing is performed through the CO decoupling strategy, the processed system layer, and discipline layer optimization models are obtained.Step 2: in a given design space, LHS sampling method is adopted to get initial uniform test sample points’ set. Then, the initial sample library is set up by the initial samples. We set the maximum number of iterationsn, stopping criterion for adding points ε and stopping criterion for iteration μ.Step 3: the sample point set is trained by the Kriging approximation model method. The Kriging approximation model is constructed according to the system layer objective function, constraints, and discipline layer constraints.Step 4: we perform uncertainty analysis on the system layer and discipline layer optimization approximation models, respectively. Then, the improved PSO optimization algorithm is used as the solver to obtain the current optimal point of the optimization model. The EI function value near the current optimal solution is judged. If it converges, we infill a point to terminate, and the current solution is the optimal solution; if it does not converge, we jump to the next step.Step 5: we use the currently calculated maximum extreme point of the EI function as the newly infill sampling point. Meantime, we update the test sample point set.Step 6: we infill the updated sample point set as a new test sample to the initial sample point set. We go to step 3 to rebuild the approximate model and judge the distance between the optimal solution of thenth step and the optimal solution of the n − 1 step. If it can converge, we stop the iteration and obtain the optimal solution.Figure 11
The solution idea of the RBMDO problem considering interval uncertainty based on the point-infilled Kriging model.


## 4.1. RBMDO Method considering Interval Uncertainty

In this study, RBMDO turns uncertainty into reliability constraints, which are added to the MDO process. In the previous sections, the RBMDO problem considering interval uncertainty has been transformed into deterministic MDO problem. Based on the traditional deterministic MDO strategy, this section studies multidisciplinary problems considering interval uncertainty.The MDO process considering interval uncertainty is similar to the traditional single-discipline interval optimization process. The difference is that the solution process of the interval uncertainty multidisciplinary optimization model is three-layer nested iterative optimization. The inner layer needs to decouple multidisciplinary problems, the middle layer analyzes the decoupled interval uncertainty model, and the outer layer solves the design variables through intelligent optimization algorithms. The multidisciplinary optimization model under interval uncertainty can be expressed as(19)FindXs,XminDVfXs,X,XsU,XU,Ys.t.hXs,X,XsU,XU,Y=aIgXs,X,XsU,XU,Y≤bIDV=Xs,X,XsU,XUXsU∈XsUL,XsUR,XU∈XUL,XUR,where XsU represents the interval uncertainty shared design variable, XsUL and XsUR are the lower and upper bounds, respectively, XU represent the interval uncertainty local design variables, XUL and XUR are the lower and upper bounds, respectively, aI represents the equality constraint interval, and bI represents the inequality constraint interval.The CO method is used to decouple the optimization model of (19). The system layer optimization model can be expressed as(20)System:FindZs,ZjminDVf=fZs,Zj,ZsU,ZUj,Y^s.t.JiZs,Zj,ZsU,ZUj,Y^≤εDV=Zs,Zj,ZsU,ZUj,Y^,where Y^ represents the coupling design variable, ZsU represents the uncertainty shared design variable at the system layer, ZUj represents the jth uncertainty design variable at the system layer, and Ji•=0 represents the discipline consistency constraint of discipline i, which can be expressed as(21)JiZs,Zij,Y^=Xsi−Zs22+Xij−Zj22+XsUi−ZsU22+XUij−ZUj22+Y•i−Y^•i22+Yi•Xsi,Xij,Y•i−Y^i•22,where Y•i represents the coupling variable to discipline i from other disciplines and Y^i·• represents the coupling variable input from discipline i. The discipline optimization model of the ith discipline in the discipline layer can be expressed as(22)Discipline−i:minDVJi=JiXsi,Xij,XsUi,XUij,Y•is.t.hiXsi,Xij,XsUi,XUij,Y•i=aiIgiXsi,Xij,XsUi,XUij,Y•i≤biIDV=Xsi,Xij,XsUi,XUij,Y•i,where XsUi represents the interval shared design variable of discipline i and XUij represents the jth interval local design variable of discipline i.


## 4.2. Interval Multidisciplinary Collaborative Optimization Method Based on the Kriging Model

In the previous section, the multidisciplinary collaborative optimization method with interval uncertainty was studied. Since the multidisciplinary systems usually have complex optimization models or difficulty in building optimization models. In this section, the approximate model technique is introduced to solve the MDO model under interval uncertainty. By approximating the system layer and discipline layer objective functions and constraints under the CO decoupling strategy, the difficulty of solving is reduced. The multidisciplinary collaborative optimization method framework under interval uncertainty based on the approximate model is shown in Figure10.Figure 10
Frame diagram of interval multidisciplinary collaborative optimization method based on the approximate model.


### 4.2.1. Constructing the MDO Model under Interval Uncertainty Based on the Kriging Approximation Model

First of all, for the system layer optimization model, the Kriging model is fitted, which is shown as follows:(23)System:minDVf^=f^Zs,Zj,ZsU,ZUj,Y^s.t.JiZs,Zj,ZsU,ZUj,Y^≤εDV=Zs,Zj,ZsU,ZUj,Y^,where f^• represents the approximate model of the objective function. The fitting function can be expressed as(24)f^Zs,Zj,ZsU,ZUj,Y^=qZs,Zj,ZsU,ZUj,Y^Tβ^+rTZs,Zj,ZsU,ZUj,Y^R−1F−Qβ^,where qi• is the basis function vector, β^ is the regression parameter matrix obtained by the optimal linear unbiased estimation, r• is the correlation vector between the sample value and the predicted value, R• is the spatial correlation function matrix, F is the response function, and Q is the polynomial matrix. The parameter representation is shown below:(25)qZs,Zj,ZsU,ZUj,Y^=q1Zs,Zj,ZsU,ZUj,Y^⋯qpZs,Zj,ZsU,ZUj,Y^β^=QTR−1Q−1QTR−1FrZs,Zj,ZsU,ZUj,Y^=RZs,Zj,ZsU,ZUj,Y^;Zs1,Z1,ZsU1,ZU1,Y^1⋯RZs,Zj,ZsU,ZUj,Y^;Zsn,Zn,ZsUn,ZUn,Y^nR=RZs1,Z1,ZsU1,ZU1,Y^1;Zs1,Z1,ZsU1,ZU1,Y^1RZs1,Z1,ZsU1,ZU1,Y^1;Zsn,Zn,ZsUn,ZUn,Y^nRZsn,Zn,ZsUn,ZUn,Y^n;Zs1,Z1,ZsU1,ZU1,Y^1RZsn,Zn,ZsUn,ZUn,Y^n;Zsn,Zn,ZsUn,ZUn,Y^nF=f^1,f^2,…f^nTQ=qTZs1,Z1,ZsU1,ZU1,Y^1qTZsn,Zn,ZsUn,ZUnT.After the Kriging model of the system layer has been constructed, new sample points are infilled to the approximate model constructed by the system layer through the EI function criterion during the optimization and solution process. We rebuild the Kriging model and enhance its fitting accuracy until the solution optimization criteria are satisfied. The approximate model of the system layer is constructed.In the process of multidisciplinary optimization, the discipline layer conducts generally discipline analysis on the constraints of the discipline itself. For the discipline layer, the Kriging approximation model is used to fit its constraints under the CO strategy. The fitted optimization model is shown as follows:(26)Discipline−i:minDVJi=JiXsi,Xij,XsUi,XUij,Y•is.t.h^iXsi,Xij,XsUi,XUij,Y•i=aiIg^iXsi,Xij,XsUi,XUij,Y•i≤biIDV=Xsi,Xij,XsUi,XUij,Y•i,where h^i• represents the equality constraints of approximate fitting of discipline i and g^i• represents the inequality constraints of approximate fitting of discipline i. The fitting process can be expressed as follows:(27)h^iXsi,Xij,XsUi,XUij,Y•i=qXsi,Xij,XsUi,XUij,Y•iTλ^+rTXsi,Xij,XsUi,XUij,Y•iR−1H−Dλ^,where qi• is the basis function vector, λ^ the regression parameter matrix obtained by the optimal linear unbiased estimation, r• is the correlation vector, R• is the spatial correlation function matrix, H is the response function, and D is the polynomial matrix.(28)g^Xsi,Xij,XsUi,XUij,Y•i=qXsi,Xij,XsUi,XUij,Y•iTω^+rTXsi,Xij,XsUi,XUij,Y•iR−1G−Pω^,where qi• is the basis function vector, ω^ is the regression parameter matrix obtained by the optimal linear unbiased estimation, r• represents the correlation vector between the sample value and the predicted value,R• represents the spatial correlation function matrix, G is the response function, and P represents the polynomial matrix.After the preliminary approximate model is constructed for the constraints at the discipline layer, new sample points are infilled to the constructed approximate model through the EI function criterion. The Kriging approximate model is reconstructed, which can enhance the construction accuracy of the approximate model. The MDO problem with the interval uncertainty optimization approximation model under the CO framework is completed.


### 4.2.2. MDO Problem considering Interval Uncertainty Solution Method Based on the Kriging Model

Through the analysis of the interval uncertainty multidisciplinary problem, the solution process can be divided into three layers of nested solution: multidisciplinary problem decoupling, interval uncertainty analysis, and objective function solution. However, due to the complexity of the multidisciplinary problem model, the use of approximate model technology to approximate the optimization model can reduce the computational difficulty in the solution process. It can also solve the problem that the model cannot converge due to the complexity of the model. From the above analysis, this section solves the RBMDO problem considering interval uncertainty based on the point-infilled Kriging model. The solution idea is shown in Figure11 and the solution steps are as follows:Step 1: we establish an interval uncertainty MDO optimization model. The hierarchical processing is performed through the CO decoupling strategy, the processed system layer, and discipline layer optimization models are obtained.Step 2: in a given design space, LHS sampling method is adopted to get initial uniform test sample points’ set. Then, the initial sample library is set up by the initial samples. We set the maximum number of iterationsn, stopping criterion for adding points ε and stopping criterion for iteration μ.Step 3: the sample point set is trained by the Kriging approximation model method. The Kriging approximation model is constructed according to the system layer objective function, constraints, and discipline layer constraints.Step 4: we perform uncertainty analysis on the system layer and discipline layer optimization approximation models, respectively. Then, the improved PSO optimization algorithm is used as the solver to obtain the current optimal point of the optimization model. The EI function value near the current optimal solution is judged. If it converges, we infill a point to terminate, and the current solution is the optimal solution; if it does not converge, we jump to the next step.Step 5: we use the currently calculated maximum extreme point of the EI function as the newly infill sampling point. Meantime, we update the test sample point set.Step 6: we infill the updated sample point set as a new test sample to the initial sample point set. We go to step 3 to rebuild the approximate model and judge the distance between the optimal solution of thenth step and the optimal solution of the n − 1 step. If it can converge, we stop the iteration and obtain the optimal solution.Figure 11
The solution idea of the RBMDO problem considering interval uncertainty based on the point-infilled Kriging model.


## 4.2.1. Constructing the MDO Model under Interval Uncertainty Based on the Kriging Approximation Model

First of all, for the system layer optimization model, the Kriging model is fitted, which is shown as follows:(23)System:minDVf^=f^Zs,Zj,ZsU,ZUj,Y^s.t.JiZs,Zj,ZsU,ZUj,Y^≤εDV=Zs,Zj,ZsU,ZUj,Y^,where f^• represents the approximate model of the objective function. The fitting function can be expressed as(24)f^Zs,Zj,ZsU,ZUj,Y^=qZs,Zj,ZsU,ZUj,Y^Tβ^+rTZs,Zj,ZsU,ZUj,Y^R−1F−Qβ^,where qi• is the basis function vector, β^ is the regression parameter matrix obtained by the optimal linear unbiased estimation, r• is the correlation vector between the sample value and the predicted value, R• is the spatial correlation function matrix, F is the response function, and Q is the polynomial matrix. The parameter representation is shown below:(25)qZs,Zj,ZsU,ZUj,Y^=q1Zs,Zj,ZsU,ZUj,Y^⋯qpZs,Zj,ZsU,ZUj,Y^β^=QTR−1Q−1QTR−1FrZs,Zj,ZsU,ZUj,Y^=RZs,Zj,ZsU,ZUj,Y^;Zs1,Z1,ZsU1,ZU1,Y^1⋯RZs,Zj,ZsU,ZUj,Y^;Zsn,Zn,ZsUn,ZUn,Y^nR=RZs1,Z1,ZsU1,ZU1,Y^1;Zs1,Z1,ZsU1,ZU1,Y^1RZs1,Z1,ZsU1,ZU1,Y^1;Zsn,Zn,ZsUn,ZUn,Y^nRZsn,Zn,ZsUn,ZUn,Y^n;Zs1,Z1,ZsU1,ZU1,Y^1RZsn,Zn,ZsUn,ZUn,Y^n;Zsn,Zn,ZsUn,ZUn,Y^nF=f^1,f^2,…f^nTQ=qTZs1,Z1,ZsU1,ZU1,Y^1qTZsn,Zn,ZsUn,ZUnT.After the Kriging model of the system layer has been constructed, new sample points are infilled to the approximate model constructed by the system layer through the EI function criterion during the optimization and solution process. We rebuild the Kriging model and enhance its fitting accuracy until the solution optimization criteria are satisfied. The approximate model of the system layer is constructed.In the process of multidisciplinary optimization, the discipline layer conducts generally discipline analysis on the constraints of the discipline itself. For the discipline layer, the Kriging approximation model is used to fit its constraints under the CO strategy. The fitted optimization model is shown as follows:(26)Discipline−i:minDVJi=JiXsi,Xij,XsUi,XUij,Y•is.t.h^iXsi,Xij,XsUi,XUij,Y•i=aiIg^iXsi,Xij,XsUi,XUij,Y•i≤biIDV=Xsi,Xij,XsUi,XUij,Y•i,where h^i• represents the equality constraints of approximate fitting of discipline i and g^i• represents the inequality constraints of approximate fitting of discipline i. The fitting process can be expressed as follows:(27)h^iXsi,Xij,XsUi,XUij,Y•i=qXsi,Xij,XsUi,XUij,Y•iTλ^+rTXsi,Xij,XsUi,XUij,Y•iR−1H−Dλ^,where qi• is the basis function vector, λ^ the regression parameter matrix obtained by the optimal linear unbiased estimation, r• is the correlation vector, R• is the spatial correlation function matrix, H is the response function, and D is the polynomial matrix.(28)g^Xsi,Xij,XsUi,XUij,Y•i=qXsi,Xij,XsUi,XUij,Y•iTω^+rTXsi,Xij,XsUi,XUij,Y•iR−1G−Pω^,where qi• is the basis function vector, ω^ is the regression parameter matrix obtained by the optimal linear unbiased estimation, r• represents the correlation vector between the sample value and the predicted value,R• represents the spatial correlation function matrix, G is the response function, and P represents the polynomial matrix.After the preliminary approximate model is constructed for the constraints at the discipline layer, new sample points are infilled to the constructed approximate model through the EI function criterion. The Kriging approximate model is reconstructed, which can enhance the construction accuracy of the approximate model. The MDO problem with the interval uncertainty optimization approximation model under the CO framework is completed.


## 4.2.2. MDO Problem considering Interval Uncertainty Solution Method Based on the Kriging Model

Through the analysis of the interval uncertainty multidisciplinary problem, the solution process can be divided into three layers of nested solution: multidisciplinary problem decoupling, interval uncertainty analysis, and objective function solution. However, due to the complexity of the multidisciplinary problem model, the use of approximate model technology to approximate the optimization model can reduce the computational difficulty in the solution process. It can also solve the problem that the model cannot converge due to the complexity of the model. From the above analysis, this section solves the RBMDO problem considering interval uncertainty based on the point-infilled Kriging model. The solution idea is shown in Figure11 and the solution steps are as follows:Step 1: we establish an interval uncertainty MDO optimization model. The hierarchical processing is performed through the CO decoupling strategy, the processed system layer, and discipline layer optimization models are obtained.Step 2: in a given design space, LHS sampling method is adopted to get initial uniform test sample points’ set. Then, the initial sample library is set up by the initial samples. We set the maximum number of iterationsn, stopping criterion for adding points ε and stopping criterion for iteration μ.Step 3: the sample point set is trained by the Kriging approximation model method. The Kriging approximation model is constructed according to the system layer objective function, constraints, and discipline layer constraints.Step 4: we perform uncertainty analysis on the system layer and discipline layer optimization approximation models, respectively. Then, the improved PSO optimization algorithm is used as the solver to obtain the current optimal point of the optimization model. The EI function value near the current optimal solution is judged. If it converges, we infill a point to terminate, and the current solution is the optimal solution; if it does not converge, we jump to the next step.Step 5: we use the currently calculated maximum extreme point of the EI function as the newly infill sampling point. Meantime, we update the test sample point set.Step 6: we infill the updated sample point set as a new test sample to the initial sample point set. We go to step 3 to rebuild the approximate model and judge the distance between the optimal solution of thenth step and the optimal solution of the n − 1 step. If it can converge, we stop the iteration and obtain the optimal solution.Figure 11
The solution idea of the RBMDO problem considering interval uncertainty based on the point-infilled Kriging model.


## 5. Example Analysis

### 5.1. Mathematical Example Analysis

This example contains two discipline coupling problems. The mathematical example model contains 8 design variables and 4 uncertain design variables. The optimization problem can be expressed as(29)minDVf=f1+f2f1=x12+x2+xu1−2,f2=x32+x42xu2x12=x52+x62,x32=x12+x62+x72+x82s.t.g1x,xu=x2+xu12−x1−2≤0,1g2x,xu=x3−2−x42xu2≤0,1g3x,xu=x5+xu32−x62≤0,0.5g4x,xu=x1+x7−2xu42−x62≤0,1.5g5x,xu=x62+x72xu42−x82≤0,20≤x≤10,xu1∈−1,1,xu2∈0.8,1.2,xu3∈−0.5,1,xu4∈0.9,1.1,DV=x,xu,where f1 and f2 represent the objective functions of discipline 1 and discipline 2, respectively, x=x1,x2,…,x7,x8 represent the design vector, and xu=xu1,xu2,xu3,xu4 represent the uncertainty design variable. For the optimization model shown in (29), the above models are decoupled from disciplines based on the CO decoupling framework. During the solution process, each iteration of the system layer and the discipline layer will be solved once. Figure 12 shows the optimization model after decoupling the above optimization model CO.Figure 12
Decoupled optimization model under the CO framework.In Figure12, z represents the system-layer design variables, zu represents the system-layer uncertainty design variables, x represents the discipline-layer design variables, and xu represents the discipline-layer uncertainty design variables. The discipline layer design vector x1 is the coupling design variable.The coupling design variables between discipline 1 and discipline 2 are decoupled by the CO method into the system layer and the discipline layer for optimization, respectively. The system layer coordinates the inconsistency among various disciplines by controlling the compatibility and consistency constraints of each discipline. The discipline layer optimizes with the compatibility and consistency constraints of various disciplines as the optimization goal. In the respective solution process of the discipline layer and the system layer because each discipline has a compatible consistency constraint function, the solution process between each discipline is independent relatively. There is no need to consider the solutions of other disciplines, just pass the respective solution results to the system layer for unified analysis and solution. This section solves the above model by constructing an approximate model, thereby proving the effectiveness of the approximate model in interval multidisciplinary problems. Before solving by approximate model optimization, we set the constraint probability to 0.98, the multiobjective function importance factor to 0.5, the maximum number of iterations to 200, and the compatibility consistency constraint toε1=ε2=0.05. In the given space of design variables, the LHS experimental design method is used to sample to obtain the initial sample points. In order to compare the optimization efficiency and optimization accuracy, 20, 30, 40, and 50 experimental design points were sampled at the system layer for fitting and 20 experimental design points were sampled at the discipline layer for approximate fitting. The Kriging approximation model is used to fit the system-layer objective function, constraints and discipline-layer constraints to build an approximate model. The EI function point-infilled criterion is used to improve the accuracy of the approximate model. The constructed approximate models at the system layer and discipline layer are shown in equations (30)–(32):System layer:(30)System:minDVF=f^z1,z2,z3,z4,zu1,zu2s.t.G1=g^1z1,z2,zu1≤0,1,G2=g^2z3,z4,zu2≤0,1J1=z1−x122+z6−x6122≤ε1,J2=z3−x322+z6−x6222≤ε20≤z≤10,zu1∈−1,1,zu2∈0.8,1.2,DV=z,zu.Discipline layer 1:(31)Discipline−1:minDVJ1=z1−x122+z6−x6122wherex1=x52+x620.5s.t.G3=g^3x5,x6,xu3≤0,0.50≤x≤10,xu3∈−0.5,1,DV=x,xu.Discipline layer 2:(32)Discipline−2:minDVJ2=z3−x322+z6−x6222wherex3=x12+x62+x72+x820.5s.t.G4=g^4x1,x6,x7,xu4≤0,1.5G5=g^5x6,x7,x8,xu4≤0,20≤x≤10,xu4∈0.9,1.1,DV=x,xu.To verify the accuracy of the method proposed in this section, the results obtained by solving the real mathematical model are shown in Table1. The approximate model is called through the multidisciplinary system solution strategy under interval uncertainty, and the approximate model is solved by the improved PSO algorithm. The solution results are shown in Table 2. It can be seen from the results that the fewer the initial sampling points, the more iterations. For the interval multidisciplinary system, the approximate model will be reconstructed multiple times at the discipline layer for each iteration. Efficiency decreases as the number of uncertainty analyses increases. For the solution results, the final solution results with different number of sampling points is small, the error from the real solution range is also small.Table 1
Solving results of the real mathematical model.
Optimal design vectorObjective function intervalx = [0.77, 1.42, 1.54, 0.71, 0.01, 0.77, 1.01, 0.40][3.493, 5.164]Table 2
Results of interval solution based on the Kriging approximation model with points infilled.
Sample timesIteration timesOptimal design vectorObjective function interval2078[0.79, 0.87, 1.60, 0.67, 0.00, 0.79, 1.14, 0.04][3.617, 5.375]3052[0.75, 1.06, 1.48, 0.71, 0.01, 0.75, 0.97, 0.35][3.645, 5.256]4036[0.76, 1.44, 1.47, 0.72, 0.00, 0.76, 0.97, 0.30][3.529, 5.224]5027[0.80, 1.29, 1.63, 0.78, 0.01, 0.79, 1.15, 0.29][3.546, 5.201]Through a mathematical example, the number of different test points is sampled to compare and solve. The solution results are compared with the real solution results, which verify the RBMDO strategy considering interval uncertainty based on the point-infilled Kriging model proposed in this section.


### 5.2. Design Optimization of Engine Piston

As the core component of the engine, the piston is one of the parts with the worst working environment in the engine. In the design process, it is necessary to ensure that it can withstand various needs such as thermal load, mechanical load, and high-speed operation reliability. Figure13(a) shows a simplified three-dimensional model of a certain type of piston and Figure 13(b) shows the corresponding parameters of the piston [38].Figure 13
Piston model diagram. (a) 3D model of piston and (b) main parameters of piston.
(a)(b)The piston material is aluminum alloy. It is optimized with the design goal of minimizing piston mass. The design variables in lightweight design usually depend on the selection of the dimensions of each part in the structural design process. Design variables and uncertainty design variables are shown in Tables3 and 4.Table 3
Piston design variables.
Design variable namesDesign variablesInitial design variable valuesDesign variable rangeOverall cylinder height,h1x1 (mm)6860≤x1≤76Fire bank height,h3x2 (mm)1210≤x2≤151st, 2nd, and 3rd ring groove heights,h4x3 (mm)4.53≤x3≤61st and 2nd bank heights,h5x4 (mm)52≤x4≤7Piston pin seat height,h8x5 (mm)2520≤x5≤30Piston pin seat width,lx6 (mm)2015≤x6≤25Top thickness,h2x7 (mm)129≤x7≤14Piston top outer diameter,d1x8 (mm)5248≤x8≤55Table 4
Uncertainty variables of piston interval.
Uncertain variable nameUncertainty variableInitial uncertainty variable valueInterval rangePiston top inner diameter,d2xu1 (mm)42[38, 46]Skirt ring groove height,h6xu2 (mm)8[6, 10]Skirt height,h7xu3 (mm)26[25, 30]In Tables3 and 4, x1,x2,…,x7,x8 is the design variable and xu1,xu2,xu3 is the interval uncertainty design variable. The piston model is analyzed by introducing a multidisciplinary strategy; this study divides the disciplines into thermal and strength disciplines. The design variables that have a greater impact on the thermal discipline are X1=x1,x2,x7,x8,xu1,xu3, which are set as the design variables of the thermal discipline. The design variable that has a greater impact on the strength discipline is X2=x3,x5,x6,x7,x8,xu1,xu2, which are set as the design variables of the strength discipline. The coupling relationship between the multiple disciplines of the piston is shown in Figure 14.Figure 14
Coupling relationship between multiple disciplines of piston.According to the construction process of the Kriging approximation model, this example uses the LHS method to sample 55 sets of data uniformly within the design space. Then, parametric modeling and finite element analysis were used to obtain sample points, 50 sets of data were used to fit the piston optimization model to construct an approximate model. The other 5 sets of data are used to test the accuracy of the constructed approximate model. After points infilled through the EI function, the average error of the approximate model is controlled within 10%. Part of the sampling data is shown in Table5. Taking the initial design variables of the piston as an example, the boundary conditions are given and the finite element analysis is performed according to the working conditions of the piston. The model is divided into finite element mesh firstly as shown in Figure 15. The mesh contains 471,307 nodes and 117,681 elements.Table 5
Partial sampling points of LHS test design.
Sample points12⋯2728⋯5455Experimental design variables (mm)x162.368.7⋯72.861.5⋯67.970.4x210.213.4⋯12.913.7⋯11.812.7x34.23.6⋯5.44.1⋯5.73.9x42.84.6⋯5.65.4⋯6.83.6x522.324.6⋯29.926.3⋯27.521.6x616.318.9⋯24.622.1⋯24.217.4x79.610.8⋯12.112.8⋯12.910.3x848.652.1⋯54.155.0⋯52.949.3xu140.242.5⋯42.144.3⋯43.838.7xu27.86.5⋯8.77.5⋯9.76.9xu326.428.3⋯28.525.9⋯29.125.7Temperature (°C)391408⋯409383⋯395403Stress (MPa)308295⋯315299⋯306304Quality (kg)0.440.50⋯0.530.42⋯0.430.51Figure 15
Finite element mesh of the piston model.The heat transfer process of the piston is very complicated during the working process. In this study, three boundary conditions of heat transfer at the top of the piston, external heat transfer, and inner cavity heat transfer are set in the heat transfer analysis. The convection heat transfer coefficient is determined by a combination of experience and theory generally. The convection heat transfer coefficientαh between the gas in the cylinder liner and the top of the piston, the convection heat transfer coefficient βh between the outside of the piston and the cooling water, and the convection heat transfer coefficient γh in the piston cavity are shown below:(33)αh=1b2/0.84+1/1/2π∫02π2.47Cm3⋅PgTgdθpβh=1.163300+1800GcScγh=T1−T2κT1−Toiδp,where Tg represents the instantaneous temperature in the cylinder, Pg represents the instantaneous pressure in the cylinder, Cm represents the average speed of the piston, and b2 represents the thickness of the thermal insulation coating. The thickness of the thermal insulation coating of the piston is 0.1 mm, Gc represents the cooling water flow, Sc represents the cross-sectional area of the coolant passage, Toi represents the crankcase temperature, T1 is the top temperature, T2 is the cavity temperature, κ is the material thermal conductivity, and δp is the thickness of the top of the piston. According to the set boundary conditions, the piston model is subjected to finite element analysis. The obtained temperature field distribution cloud diagram is shown in Figure 16. The temperature at the center of the top of the piston is the highest, with the highest point being 329.05°C. The temperature at the bottom of the piston skirt is the lowest, with the lowest point being 118.36°C.Figure 16
Cloud diagram of temperature field distribution.After the thermal analysis is completed, the force analysis of the model is performed. We apply air pressure and inertial force to the model based on the above thermal analysis model. According to [38], it can be determined that the air pressure at the top of the piston and the side of the fire bank height is 1.2–1.6 MPa, generally. The pressure at the first piston ring port is about 0.76 P, the pressure at the second ring port is about 0.25 P, and other places are ignored. When the piston moves to the top dead center and bottom dead center, it will be affected by inertial force, which can be calculated by(34)Fp=mpωp2Rpcosαp+RpLpcos2αp,where mp represents the piston mass, Rp represents the piston radius, Lp represents the connecting rod length, ωp represents the crank rotation angular velocity, and αp represents the crank angle. After the above boundary conditions are applied, the stress and deformation analysis is performed. The cloud diagram of the analysis results is shown in Figure 17. It can be seen that the stress at the center of the top of the piston is the largest, with a maximum of 353.47 MPa.Figure 17
Coupling stress cloud diagram of piston.According to the RBMDO strategy considering interval uncertainty based on the point-infilled Kriging model proposed in this study, taking the lightweight of the piston as the optimization goal, discipline 1 is the thermal discipline and discipline 2 is the strength discipline. According to the above experimental design points, the approximate fitting method based on the point-infilled Kriging criterion is used under the CO decoupling strategy to construct approximate models for the system layer and the discipline layer, respectively. The constructed approximate mathematical optimization model can be expressed as equations (35)–(37).System layer:(35)System:minFZ,Zu=f^Z,Zus.t.JiZ,Zu≤εZ=z1,z2,z3,z4,z5,z6,z7,z8Zu=zu1,zu1,zu3.Discipline layer 1:(36)Discipline1:minJ1=x1−z12+x2−z22+x7−z72+x8−z82+xu1−zu12+xu3−zu32s.t.G1=T=T^x1,x2,x7,x8,xu1,xu3≤350,37060≤x1≤76,10≤x2≤159≤x7≤14,48≤x8≤55xu1∈38,46,xu3∈25,30.Discipline layer 2:(37)Discipline2:minJ2=x3−z32+x5−z52+x6−z62+x7−z72+x8−z82+xu1−zu12+xu2−zu22s.t.G2=σ=σ^x3,x5,x6,x7,x8,xu1,xu2,T≤400,4403≤x3≤6,20≤x5≤30,15≤x6≤259≤x7≤14,48≤x8≤55xu1∈38,46,xu2∈6,10,where F• represents the objective function, Z represents the system layer design variable, Zu represents the system layer uncertainty design variable., x=x1,x2,…,x7,x8 represents the discipline layer design variable, xu=xu1,xu2,xu3 represents the discipline layer uncertainty design variable, andJi• represents the compatibility consistency constraint.Safety and reliability need to be ensured during the working process of the piston. In the calculation process, both the system layer and the discipline layer constraint probability are set to 0.98. The objective function importance factor is set to 0. The compatibility consistency constraint is set toε1=ε2=0.05. The system layer iteration number is set to 300. Under the CO framework, the uncertainty analysis is performed by using the Kriging approximation model and the optimization algorithm is used to solve the problem. Subsequently, the iterative optimization solution is performed through the EI point-infilled criterion. The results obtained by the two models are shown in Table 6.Table 6
Optimization results of the piston model.
Namex1x2x3x4x5x6Initial value68124.552520Kriging67.1213.264.335.4623.9219.69Point-infilled Kriging66.8413.184.265.3423.7819.16Namex7x8TmaxσmaxmInitial value1252329.05353.470.481Kriging11.8950.97[309, 320][368, 385][0.458, 0.473]Point-infilled kriging11.6750.49[315, 322][379, 390][0.449, 0.461]It can be seen that the quality of the optimized results of the Kriging model based on the point-infilled criterion is lower than the traditional Kriging model. The maximum mass value is reduced by 9% compared with the initial value and the maximum temperature is reduced by 7°C. The maximum stress value is increased compared with the initial value, but the maximum allowable stress value is not exceeded, which can satisfy the actual engineering.


## 5.1. Mathematical Example Analysis

This example contains two discipline coupling problems. The mathematical example model contains 8 design variables and 4 uncertain design variables. The optimization problem can be expressed as(29)minDVf=f1+f2f1=x12+x2+xu1−2,f2=x32+x42xu2x12=x52+x62,x32=x12+x62+x72+x82s.t.g1x,xu=x2+xu12−x1−2≤0,1g2x,xu=x3−2−x42xu2≤0,1g3x,xu=x5+xu32−x62≤0,0.5g4x,xu=x1+x7−2xu42−x62≤0,1.5g5x,xu=x62+x72xu42−x82≤0,20≤x≤10,xu1∈−1,1,xu2∈0.8,1.2,xu3∈−0.5,1,xu4∈0.9,1.1,DV=x,xu,where f1 and f2 represent the objective functions of discipline 1 and discipline 2, respectively, x=x1,x2,…,x7,x8 represent the design vector, and xu=xu1,xu2,xu3,xu4 represent the uncertainty design variable. For the optimization model shown in (29), the above models are decoupled from disciplines based on the CO decoupling framework. During the solution process, each iteration of the system layer and the discipline layer will be solved once. Figure 12 shows the optimization model after decoupling the above optimization model CO.Figure 12
Decoupled optimization model under the CO framework.In Figure12, z represents the system-layer design variables, zu represents the system-layer uncertainty design variables, x represents the discipline-layer design variables, and xu represents the discipline-layer uncertainty design variables. The discipline layer design vector x1 is the coupling design variable.The coupling design variables between discipline 1 and discipline 2 are decoupled by the CO method into the system layer and the discipline layer for optimization, respectively. The system layer coordinates the inconsistency among various disciplines by controlling the compatibility and consistency constraints of each discipline. The discipline layer optimizes with the compatibility and consistency constraints of various disciplines as the optimization goal. In the respective solution process of the discipline layer and the system layer because each discipline has a compatible consistency constraint function, the solution process between each discipline is independent relatively. There is no need to consider the solutions of other disciplines, just pass the respective solution results to the system layer for unified analysis and solution. This section solves the above model by constructing an approximate model, thereby proving the effectiveness of the approximate model in interval multidisciplinary problems. Before solving by approximate model optimization, we set the constraint probability to 0.98, the multiobjective function importance factor to 0.5, the maximum number of iterations to 200, and the compatibility consistency constraint toε1=ε2=0.05. In the given space of design variables, the LHS experimental design method is used to sample to obtain the initial sample points. In order to compare the optimization efficiency and optimization accuracy, 20, 30, 40, and 50 experimental design points were sampled at the system layer for fitting and 20 experimental design points were sampled at the discipline layer for approximate fitting. The Kriging approximation model is used to fit the system-layer objective function, constraints and discipline-layer constraints to build an approximate model. The EI function point-infilled criterion is used to improve the accuracy of the approximate model. The constructed approximate models at the system layer and discipline layer are shown in equations (30)–(32):System layer:(30)System:minDVF=f^z1,z2,z3,z4,zu1,zu2s.t.G1=g^1z1,z2,zu1≤0,1,G2=g^2z3,z4,zu2≤0,1J1=z1−x122+z6−x6122≤ε1,J2=z3−x322+z6−x6222≤ε20≤z≤10,zu1∈−1,1,zu2∈0.8,1.2,DV=z,zu.Discipline layer 1:(31)Discipline−1:minDVJ1=z1−x122+z6−x6122wherex1=x52+x620.5s.t.G3=g^3x5,x6,xu3≤0,0.50≤x≤10,xu3∈−0.5,1,DV=x,xu.Discipline layer 2:(32)Discipline−2:minDVJ2=z3−x322+z6−x6222wherex3=x12+x62+x72+x820.5s.t.G4=g^4x1,x6,x7,xu4≤0,1.5G5=g^5x6,x7,x8,xu4≤0,20≤x≤10,xu4∈0.9,1.1,DV=x,xu.To verify the accuracy of the method proposed in this section, the results obtained by solving the real mathematical model are shown in Table1. The approximate model is called through the multidisciplinary system solution strategy under interval uncertainty, and the approximate model is solved by the improved PSO algorithm. The solution results are shown in Table 2. It can be seen from the results that the fewer the initial sampling points, the more iterations. For the interval multidisciplinary system, the approximate model will be reconstructed multiple times at the discipline layer for each iteration. Efficiency decreases as the number of uncertainty analyses increases. For the solution results, the final solution results with different number of sampling points is small, the error from the real solution range is also small.Table 1
Solving results of the real mathematical model.
Optimal design vectorObjective function intervalx = [0.77, 1.42, 1.54, 0.71, 0.01, 0.77, 1.01, 0.40][3.493, 5.164]Table 2
Results of interval solution based on the Kriging approximation model with points infilled.
Sample timesIteration timesOptimal design vectorObjective function interval2078[0.79, 0.87, 1.60, 0.67, 0.00, 0.79, 1.14, 0.04][3.617, 5.375]3052[0.75, 1.06, 1.48, 0.71, 0.01, 0.75, 0.97, 0.35][3.645, 5.256]4036[0.76, 1.44, 1.47, 0.72, 0.00, 0.76, 0.97, 0.30][3.529, 5.224]5027[0.80, 1.29, 1.63, 0.78, 0.01, 0.79, 1.15, 0.29][3.546, 5.201]Through a mathematical example, the number of different test points is sampled to compare and solve. The solution results are compared with the real solution results, which verify the RBMDO strategy considering interval uncertainty based on the point-infilled Kriging model proposed in this section.


## 5.2. Design Optimization of Engine Piston

As the core component of the engine, the piston is one of the parts with the worst working environment in the engine. In the design process, it is necessary to ensure that it can withstand various needs such as thermal load, mechanical load, and high-speed operation reliability. Figure13(a) shows a simplified three-dimensional model of a certain type of piston and Figure 13(b) shows the corresponding parameters of the piston [38].Figure 13
Piston model diagram. (a) 3D model of piston and (b) main parameters of piston.
(a)(b)The piston material is aluminum alloy. It is optimized with the design goal of minimizing piston mass. The design variables in lightweight design usually depend on the selection of the dimensions of each part in the structural design process. Design variables and uncertainty design variables are shown in Tables3 and 4.Table 3
Piston design variables.
Design variable namesDesign variablesInitial design variable valuesDesign variable rangeOverall cylinder height,h1x1 (mm)6860≤x1≤76Fire bank height,h3x2 (mm)1210≤x2≤151st, 2nd, and 3rd ring groove heights,h4x3 (mm)4.53≤x3≤61st and 2nd bank heights,h5x4 (mm)52≤x4≤7Piston pin seat height,h8x5 (mm)2520≤x5≤30Piston pin seat width,lx6 (mm)2015≤x6≤25Top thickness,h2x7 (mm)129≤x7≤14Piston top outer diameter,d1x8 (mm)5248≤x8≤55Table 4
Uncertainty variables of piston interval.
Uncertain variable nameUncertainty variableInitial uncertainty variable valueInterval rangePiston top inner diameter,d2xu1 (mm)42[38, 46]Skirt ring groove height,h6xu2 (mm)8[6, 10]Skirt height,h7xu3 (mm)26[25, 30]In Tables3 and 4, x1,x2,…,x7,x8 is the design variable and xu1,xu2,xu3 is the interval uncertainty design variable. The piston model is analyzed by introducing a multidisciplinary strategy; this study divides the disciplines into thermal and strength disciplines. The design variables that have a greater impact on the thermal discipline are X1=x1,x2,x7,x8,xu1,xu3, which are set as the design variables of the thermal discipline. The design variable that has a greater impact on the strength discipline is X2=x3,x5,x6,x7,x8,xu1,xu2, which are set as the design variables of the strength discipline. The coupling relationship between the multiple disciplines of the piston is shown in Figure 14.Figure 14
Coupling relationship between multiple disciplines of piston.According to the construction process of the Kriging approximation model, this example uses the LHS method to sample 55 sets of data uniformly within the design space. Then, parametric modeling and finite element analysis were used to obtain sample points, 50 sets of data were used to fit the piston optimization model to construct an approximate model. The other 5 sets of data are used to test the accuracy of the constructed approximate model. After points infilled through the EI function, the average error of the approximate model is controlled within 10%. Part of the sampling data is shown in Table5. Taking the initial design variables of the piston as an example, the boundary conditions are given and the finite element analysis is performed according to the working conditions of the piston. The model is divided into finite element mesh firstly as shown in Figure 15. The mesh contains 471,307 nodes and 117,681 elements.Table 5
Partial sampling points of LHS test design.
Sample points12⋯2728⋯5455Experimental design variables (mm)x162.368.7⋯72.861.5⋯67.970.4x210.213.4⋯12.913.7⋯11.812.7x34.23.6⋯5.44.1⋯5.73.9x42.84.6⋯5.65.4⋯6.83.6x522.324.6⋯29.926.3⋯27.521.6x616.318.9⋯24.622.1⋯24.217.4x79.610.8⋯12.112.8⋯12.910.3x848.652.1⋯54.155.0⋯52.949.3xu140.242.5⋯42.144.3⋯43.838.7xu27.86.5⋯8.77.5⋯9.76.9xu326.428.3⋯28.525.9⋯29.125.7Temperature (°C)391408⋯409383⋯395403Stress (MPa)308295⋯315299⋯306304Quality (kg)0.440.50⋯0.530.42⋯0.430.51Figure 15
Finite element mesh of the piston model.The heat transfer process of the piston is very complicated during the working process. In this study, three boundary conditions of heat transfer at the top of the piston, external heat transfer, and inner cavity heat transfer are set in the heat transfer analysis. The convection heat transfer coefficient is determined by a combination of experience and theory generally. The convection heat transfer coefficientαh between the gas in the cylinder liner and the top of the piston, the convection heat transfer coefficient βh between the outside of the piston and the cooling water, and the convection heat transfer coefficient γh in the piston cavity are shown below:(33)αh=1b2/0.84+1/1/2π∫02π2.47Cm3⋅PgTgdθpβh=1.163300+1800GcScγh=T1−T2κT1−Toiδp,where Tg represents the instantaneous temperature in the cylinder, Pg represents the instantaneous pressure in the cylinder, Cm represents the average speed of the piston, and b2 represents the thickness of the thermal insulation coating. The thickness of the thermal insulation coating of the piston is 0.1 mm, Gc represents the cooling water flow, Sc represents the cross-sectional area of the coolant passage, Toi represents the crankcase temperature, T1 is the top temperature, T2 is the cavity temperature, κ is the material thermal conductivity, and δp is the thickness of the top of the piston. According to the set boundary conditions, the piston model is subjected to finite element analysis. The obtained temperature field distribution cloud diagram is shown in Figure 16. The temperature at the center of the top of the piston is the highest, with the highest point being 329.05°C. The temperature at the bottom of the piston skirt is the lowest, with the lowest point being 118.36°C.Figure 16
Cloud diagram of temperature field distribution.After the thermal analysis is completed, the force analysis of the model is performed. We apply air pressure and inertial force to the model based on the above thermal analysis model. According to [38], it can be determined that the air pressure at the top of the piston and the side of the fire bank height is 1.2–1.6 MPa, generally. The pressure at the first piston ring port is about 0.76 P, the pressure at the second ring port is about 0.25 P, and other places are ignored. When the piston moves to the top dead center and bottom dead center, it will be affected by inertial force, which can be calculated by(34)Fp=mpωp2Rpcosαp+RpLpcos2αp,where mp represents the piston mass, Rp represents the piston radius, Lp represents the connecting rod length, ωp represents the crank rotation angular velocity, and αp represents the crank angle. After the above boundary conditions are applied, the stress and deformation analysis is performed. The cloud diagram of the analysis results is shown in Figure 17. It can be seen that the stress at the center of the top of the piston is the largest, with a maximum of 353.47 MPa.Figure 17
Coupling stress cloud diagram of piston.According to the RBMDO strategy considering interval uncertainty based on the point-infilled Kriging model proposed in this study, taking the lightweight of the piston as the optimization goal, discipline 1 is the thermal discipline and discipline 2 is the strength discipline. According to the above experimental design points, the approximate fitting method based on the point-infilled Kriging criterion is used under the CO decoupling strategy to construct approximate models for the system layer and the discipline layer, respectively. The constructed approximate mathematical optimization model can be expressed as equations (35)–(37).System layer:(35)System:minFZ,Zu=f^Z,Zus.t.JiZ,Zu≤εZ=z1,z2,z3,z4,z5,z6,z7,z8Zu=zu1,zu1,zu3.Discipline layer 1:(36)Discipline1:minJ1=x1−z12+x2−z22+x7−z72+x8−z82+xu1−zu12+xu3−zu32s.t.G1=T=T^x1,x2,x7,x8,xu1,xu3≤350,37060≤x1≤76,10≤x2≤159≤x7≤14,48≤x8≤55xu1∈38,46,xu3∈25,30.Discipline layer 2:(37)Discipline2:minJ2=x3−z32+x5−z52+x6−z62+x7−z72+x8−z82+xu1−zu12+xu2−zu22s.t.G2=σ=σ^x3,x5,x6,x7,x8,xu1,xu2,T≤400,4403≤x3≤6,20≤x5≤30,15≤x6≤259≤x7≤14,48≤x8≤55xu1∈38,46,xu2∈6,10,where F• represents the objective function, Z represents the system layer design variable, Zu represents the system layer uncertainty design variable., x=x1,x2,…,x7,x8 represents the discipline layer design variable, xu=xu1,xu2,xu3 represents the discipline layer uncertainty design variable, andJi• represents the compatibility consistency constraint.Safety and reliability need to be ensured during the working process of the piston. In the calculation process, both the system layer and the discipline layer constraint probability are set to 0.98. The objective function importance factor is set to 0. The compatibility consistency constraint is set toε1=ε2=0.05. The system layer iteration number is set to 300. Under the CO framework, the uncertainty analysis is performed by using the Kriging approximation model and the optimization algorithm is used to solve the problem. Subsequently, the iterative optimization solution is performed through the EI point-infilled criterion. The results obtained by the two models are shown in Table 6.Table 6
Optimization results of the piston model.
Namex1x2x3x4x5x6Initial value68124.552520Kriging67.1213.264.335.4623.9219.69Point-infilled Kriging66.8413.184.265.3423.7819.16Namex7x8TmaxσmaxmInitial value1252329.05353.470.481Kriging11.8950.97[309, 320][368, 385][0.458, 0.473]Point-infilled kriging11.6750.49[315, 322][379, 390][0.449, 0.461]It can be seen that the quality of the optimized results of the Kriging model based on the point-infilled criterion is lower than the traditional Kriging model. The maximum mass value is reduced by 9% compared with the initial value and the maximum temperature is reduced by 7°C. The maximum stress value is increased compared with the initial value, but the maximum allowable stress value is not exceeded, which can satisfy the actual engineering.


## 6. Discussion

With the increase in structure and system complexity in engineering, there are more and more uncertain factors affecting the stable operation of the system. In the design process, the uncertainty of system should be considered in many aspects. As an important part of nonprobabilistic uncertainty design optimization method, the interval number method has shown a good application prospect in many practical projects. In this study, the interval uncertainty design optimization method based on the approximate model and the interval multidisciplinary design optimization method based on the approximate model are studied. The main work is as follows:(1)
The optimization model considering interval uncertainty and its solution method are studied. In the process of solving the optimization model, the ideal point method and penalty function method in the evaluation function method are combined to calculate the optimization model considering uncertainty, and the interval uncertainty model is solved by the double-nested optimization algorithm based on improved PSO.(2)
The Kriging approximation model based on point-infilled criterion is studied. The optimization method based on the Kriging approximation model can simplify the complex mathematical optimization model. The approximate model method is applied to the interval uncertainty optimization problem to solve the low efficiency of the double-nested optimization in the decoupling process of the interval uncertainty problem.(3)
The MDO decoupling strategy based on CO considering interval uncertainty and its optimization solution method are studied. Considering the MDO model and the decoupling of discipline relationship, the coupling relationship between various disciplines based on the framework of CO solution are studied, and an optimization model considering interval uncertainty based on the point-infilled Kriging approximation model is constructed. Finally, a mathematical example and an engineering example of piston design optimization are used to verify the effectiveness of the interval MDO method based on the Kriging approximation model.


### 6.1. Research Contribution

This thesis mainly studies three aspects based on interval uncertainty theory. First, the construction of the interval uncertainty model and the transformation and solution of the uncertainty model are studied. Second, a new idea of solving the interval uncertainty optimization model is further proposed based on the approximate model technology. Third, the approximate model and the interval uncertainty theory are combined with the multidisciplinary optimization method. By combining the interval uncertainty analysis and optimization method with the approximate model and MDO problem, it provides a new way for the design optimization of large-scale engineering equipment and system and also provides a new method for the theory of interval uncertainty design optimization.


### 6.2. Limitations and Future Research Directions

(1)
At present, in interval uncertainty optimization, to use the deterministic optimization solution method, the uncertainty optimization model needs to be converted into a deterministic optimization model. However, it translates to a two-level nested optimization problem or even three-level nested optimization problem in the RBMDO problem. With the increasing complexity of the problem, the number of computer iterations increases exponentially during the solution process, which will affect the computational efficiency greatly. Therefore, it is necessary to study some new methods to solve the interval uncertainty problems, so as to ensure the generality of interval optimization methods in more complex problems.(2)
This paper studies the problem of nonprobabilistic uncertainty under interval uncertainty. However, many systems also have uncertainty factors based on probability distribution. Therefore, it needs to be studied in the future that how to quantify both the uncertainty factors based on probability distribution and based on nonprobability distribution in the process of optimization.


## 6.1. Research Contribution

This thesis mainly studies three aspects based on interval uncertainty theory. First, the construction of the interval uncertainty model and the transformation and solution of the uncertainty model are studied. Second, a new idea of solving the interval uncertainty optimization model is further proposed based on the approximate model technology. Third, the approximate model and the interval uncertainty theory are combined with the multidisciplinary optimization method. By combining the interval uncertainty analysis and optimization method with the approximate model and MDO problem, it provides a new way for the design optimization of large-scale engineering equipment and system and also provides a new method for the theory of interval uncertainty design optimization.


## 6.2. Limitations and Future Research Directions

(1)
At present, in interval uncertainty optimization, to use the deterministic optimization solution method, the uncertainty optimization model needs to be converted into a deterministic optimization model. However, it translates to a two-level nested optimization problem or even three-level nested optimization problem in the RBMDO problem. With the increasing complexity of the problem, the number of computer iterations increases exponentially during the solution process, which will affect the computational efficiency greatly. Therefore, it is necessary to study some new methods to solve the interval uncertainty problems, so as to ensure the generality of interval optimization methods in more complex problems.(2)
This paper studies the problem of nonprobabilistic uncertainty under interval uncertainty. However, many systems also have uncertainty factors based on probability distribution. Therefore, it needs to be studied in the future that how to quantify both the uncertainty factors based on probability distribution and based on nonprobability distribution in the process of optimization.


## 7. Conclusion

This study mainly focuses on three aspects: interval uncertainty optimization model and its solution method, interval uncertainty design optimization method based on the approximate model, and interval multidisciplinary design optimization method. The main conclusions of this study are as follows:(1)
The interval uncertainty optimization model is established and the solution method is provided. The ideal point method of the evaluation function method is used to deal with the multidisciplinary optimization model transformed from interval uncertainty. Then, an improved PSO algorithm is introduced to solve the double-nested optimization model.(2)
The Kriging model based on the point-infilled strategy is established. The test results show that the proposed model not only improves the computation efficiency but also can guarantee high accuracy of the approximation model.(3)
A RBMDO strategy considering interval uncertainty based on the point-infilled Kriging model is established. Firstly, MDO is introduced into the process of interval uncertainty optimization, and an optimization method based on interval multidisciplinary collaborative is proposed. Then, an interval multidisciplinary design optimization method based on the Kriging approximation model is built by combining the interval optimization method based on point-filled Kriging strategy and the interval multidisciplinary collaborative optimization method. Thirdly, in order to effectively solve the interval multidisciplinary approximate optimization model, a three-level nested interval multidisciplinary design optimization solution method is proposed. In the end, a mathematical example is used to verify the effectiveness of the proposed method. According to the example of the engine piston, the quality of the optimized results of the Kriging model based on the point-infilled criterion is lower than the traditional Kriging model. The maximum mass value is reduced by 9% compared with the initial value and the maximum temperature is reduced by 7°C. The maximum stress value is increased compared with the initial value, but the maximum allowable stress value is not exceeded, which can satisfy the actual engineering. The example verifies that the method proposed in this study is superior compared with the traditional method.



---

*Source: 8582511-2022-11-17.xml*