A Least Squares Ensemble Model Based on Regularization and Augmentation Strategy

Surrogate models are often used as alternatives to considerably reduce the computational burden of the expensive computer simulations that are required for engineering designs. The development of surrogate models for complex relationships between the parameters often requires the modeling of high-dimensional functions with limited information, and it is challenging to choose an effective surrogate model over the unknown design space. To this end, the ensemble models—combined with different surrogate models—offer effective solutions. This paper presents a new ensemble model based on the least squares method, which is a regularization strategy and an augmentation strategy; we call the model the regularized least squares ensemble model (RLS-EM). Three individual surrogate models—Kriging, radial basis function, and support vector regression—are used to compose the RLS-EM. Further, the weight factors are estimated by the least squares method without using the global or local error metrics, which are used in most existing methods. To solve the collinearity in the least squares calculation process, a regularization strategy and an augmentation strategy are developed. The two strategies help explore the unknown regions and improve the accuracy on one hand; on the other hand, the collinearity can be reduced, and the overfitting phenomenon that may occur can be avoided. Six numerical functions, from two-dimensional to 12-dimensional, and a computer numerical control (CNC) milling machine bed design problem are used to verify the proposed method. The results of the numerical examples show that RLS-EM saves a considerable amount of computation time while ensuring the same level of robustness and accuracy compared with other ensemble models. The RLS-EM used for the CNC milling machine bed design problem also shows good accuracy characteristics compared with other ensemble methods.


Introduction
Computational simulations, such as finite element analysis (FEA) or computational fluid dynamics, have been displaying steady progress in describing engineering systems, and these simulations play a key role in optimizing the design of complex engineering equipment.However, computer simulations may consume a considerable amount of time for complicated simulations in engineering design.Therefore, a surrogate modeling method has been developed rapidly over the last three decades as an alternative for computationally expensive simulations that consumes less time [1].A wide variety of surrogate models have been used in engineering design, such as polynomial response surface (PRS) [2,3], Kriging (KRG) [3][4][5][6], radial basis function (RBF) [7,8], and support vector regression (SVR) [8][9][10][11][12][13].The PRS and SVR models can identify global trends for a given input data set; whereas, owing to the interpolation characteristics, KRG and RBF have higher local accuracy around the training points.Reviews about the surrogate models can be found in [14][15][16][17].
The rapid development of various surrogate models provides researchers with a lot more flexibility while they are selecting models for engineering design problems.However, it is challenging to choose the optimal model for a specific application before all the different surrogate models are constructed [18].Since practical engineering applications often exhibit different linear and nonlinear properties, no single surrogate model can exhibit high performance in all scenarios; each surrogate model has advantages and disadvantages [19][20][21][22].The PRS model fits the linear relationship between the inputs and the outputs well, while the KRG and RBF model are more suited for complex nonlinear relationships between the input and output datasets.The SVR model is suitable for both linear and nonlinear relationships between the inputs and the outputs, as it can choose different kernel functions and hyper-parameters [23][24][25][26][27][28][29][30].To make better use of the advantages of each model, as well as avoid wasting existing individual models, researchers have combined different surrogate models into an ensemble model to develop weighted average surrogate (WAS) models [31,32].
Several ensemble models have been developed in the literature, and studies have shown that the ensemble model combines the predictive power of each individual surrogate model to improve accuracy and robustness.Existing ensemble models are most commonly based on error correlation or prediction variance, and they can be classified according to global and local error metrics [31][32][33][34].Zerpa et al. [33] constructed a WAS using several individual surrogate models for the optimization of alkaline-surfactant-polymer flooding processes, and found that the WAS exhibited better performance than individual surrogates.Goel et al. [31] proposed different approaches in which the weights of the ensemble were determined based on the generalized mean square cross-validation error (GMSE).Acar et al. [32] proposed an optimization method to calculate the weight factors by minimizing the GMSE.Zhou et al. [34] used a recursive process to obtain the weight factors, updating them in each iteration until the convergence goal was reached.The study described above shows that the weight factors are evaluated as a global metric.Unlike the global error metric, the local error metric method assesses weight factors in relatively small spaces, even point-by-point [35][36][37][38][39]. Acar [35] used various local error measures to construct an ensemble model and presented a local error measure of pointwise cross-validation error.Lee et al. [36] presented a v-nearest points cross-validation method to calculate the weight factors in a local region.Hierarchical design space reduction [37,40], hybrid, and adaptive meta-modeling [38] are also effective methods based on local error metrics.
The essence of an ensemble model is to assign weight factors to the known individual surrogate models and then sum the results of each model.The accuracy of each model affects the accuracy of the ensemble model; thus far, individual surrogate models are expected to have relatively high accuracy to meet the robustness and the predictive performance requirements.When high-precision individual models are obtained, the weight factors can be calculated by regression methods to improve computational efficiency and save computational cost, instead of other optimization algorithms or error metrics.
Motivated by the regression idea, this paper proposes a novel ensemble modeling technique named the regularized least squares ensemble model (RLS-EM).Three individual surrogate models, KRG, RBF, and SVR, are used to develop the RLS-EM; the least squares algorithm with a regularization strategy and an augmentation strategy is used to calculate the weight factors.The regularization strategy and augmentation strategy help to solve the collinearity problem caused by the inherent interpolation properties of the KRG and RBF models, and the similar prediction values at some sample points.The augmentation strategy is carried out on the unexplored regions, which can help improve the accuracy of the surrogate models, while the regularization strategy helps to avoid the potential overfitting phenomenon.The RLS-EM aims to take advantage of the well-performing ensemble surrogate model to guarantee the robustness and accuracy for different problems from low to relatively high dimensions.The RLS-EM aims to take advantage of the well-performing ensemble surrogate model to guarantee the robustness and accuracy for different problems from low to relatively high dimensions.
The remainder of this paper is organized as follows.In the next section, a brief introduction to the ensemble methods is presented.Then, the development of the proposed RLS-EM is described.
Several numerical functions and an engineering application are tested in the following section.Finally, the conclusions are presented.

Background of Ensemble Methods
Usually, the surrogate model technique is utilized to construct several different surrogates and select the best one.However, this scenario has two major shortcomings [41].It is wasteful to discard the so-called inaccurate models, and the accuracy of the surrogate model is affected by the sample points.This is because the surrogate model may exhibit different precisions for different data sets.To overcome these drawbacks, ensemble methods are proposed.
An ensemble surrogate model is a weighted combination of several individual surrogate models [42], which is defined as: where f e (x) is the prediction of the ensemble, M is the number of surrogate models used, and w i is the weight factor for the ith surrogate f i (x).Evidently, the larger weights are assigned to the more accurate surrogate models, and vice versa.Zerpa et al. [33] proposed the evaluation of the weight factors w i in a linear ensemble as: where V i is the prediction variance of the ith surrogate model.Goel et al. [31] considered PRS, KRG, and RBF, and proposed an ensemble scheme to estimate the weight factors in a WAS, including the BestPRESS (BP), the PRESS weighted surrogate (PWS), and the non-parametric PRESS weighted surrogate (NPWS).
Taking the prediction sum of squares (PRESS) as the error measure, the NPWS is given as: where E j is the GMSE of the ith surrogate model calculated from: where y(x i ) is the true response at the ith data point x i , and ŷ(−i) j (x i ) is the corresponding prediction from the surrogate model constructed using all except the ith data point x i , and N is the number of sample points.
The model with the least PRESS error is assigned a weight factor of one, and all the other models are assigned zero weight factors; this strategy is called the BP model [31].
The PWS uses the GMSE as a global error metric to select the weight factors using a heuristic formulation, which is formulated as follows: The weighting scheme requires the user to specify parameters α and β, which control the contribution of the individual surrogates; α and β are assumed to be 0.05 and −1, respectively [31].
Acar and Rais-Rohani [32] used GMSE as the global error metric and proposed an optimization algorithm to calculate the weight factors; the algorithm is expressed as follows: Viana et al. [43] proposed an ensemble surrogate model called optimal weighted surrogate (OWS); the OWS is represented as follows: The correlation matrix of the error from the individual surrogate models that are used to constitute the ensemble surrogate model is expressed as follows: where e i and e j are the vectors of cross-validation errors (i.e., PRESS) for the ith and jth surrogate models, respectively.The application of the ensemble models can be found in [44][45][46][47].

Basic Formulation of the Least Squares Method
A general linear regression model can be represented as follows [47]: where p i (x) represents any function about the variable x or simply the variable x, M is the number of regression terms, and ε is the approximation error.For convenience, p i (x) is characterized with X i , for N samples with M dimensions, X = [p 1 (x i ), p 2 (x i ), . . ., p M (x i )], and the corresponding responses Y = [y 1 , y 2 , . . ., y N ] T .Then, the matrix form of linear regression is represented as: where w = [w 1 , w 2 , . . ., w M ] T , and the error term ε = [ε 1 , ε 2 , . . ., ε N ] T , supposing that the errors are normally and independently distributed, with zero mean and finite variance, that is ε i ∼ N(0, σ 2 ).
Based on the Gauss-Markov theorem [47], the weight factors ŵ = [ ŵ1 , ŵ2 , • • • , ŵM ] calculated by the OLS method form the best linear unbiased estimator, which can be represented as: and ŵ satisfies the following equations:

Samples Adding by the Augmentation Strategy
The RLS-EM proposed in this paper seeks to simultaneously capture the global and local accuracy.Since the PRS may exhibit lower accuracy in some nonlinear applications, RLS-EM only combines KRG, RBF, and SVR to meet the local accuracy and global trend requirements.The predicted values of the KRG and RBF models at the training points are equal to the actual function values, so the collinearity is unavoidable.To solve the collinearity problem, an augmentation strategy and a regularization strategy were developed.The augmentation strategy is used to reduce the influence of the collinearity on one hand; on the other hand, the augmentation strategy helps to improve the accuracy of the model in the unexplored area.
The N samples obtained by Latin hypercube sampling (LHS) technique are used to construct the KRG, RBF, and SVR surrogate models, and the corresponding prediction values at the samples are

P kr (x
where the P kr (x i ) is the absolute value of the difference between the KRG and RBF model at the ith adding points.The prediction values of the KRG, RBF, and SVR surrogate models on X add are represented as Ŷadd = [ ŷ1 (x iadd ), ŷ2 (x iadd ), ŷ3 (x iadd )], iadd = 1, 2, . . ., N add , and the actual functions values at the N add points are Y = [y N+1 , y N+2 , . . ., y N+Nadd ] T .Then, the augmentation matrix system of Ŷexpand and Y are represented as: Appl.Sci.2019, 9, 1845

The Regularization Strategy in the Least Squares System
A regularization term is added to further reduce the impact of collinearity.Due to the interpolation properties of the KRG and RBF models, ŷ1 (x i ) and ŷ2 (x i ) are equal to the actual function values at the N samples.The relatively high precision of KRG, RBF, and SVR surrogate models may also predict approximately equal values at some samples in the set X add .By adding a regularization term multiplying an identity matrix, the matrix coefficients ŵexpand can be estimated from the augmented matrix inversion system as follows: where I 3×3 is an identity matrix, and λ is the regularization parameter.Since the linear correlation in ŶT expand is expected to be lower than that in Ŷ, and the regularization item further reduces the linear correlation, the accuracy and robustness on evaluating w by means of ŵexpand is expected to be better.The weight factors for the ensemble model are calculated as: The regularized least squares ensemble method can be expressed as follows: Calculate the inverse of the augmented matrix system for ŵexpand by Equation ( 15), and the standardized weight factors by Equation ( 16).
However, the regularization parameter λ should be confirmed before using Equation (15); a search algorithm was developed to obtain the optimal regularization parameter value λ * , and the detailed pseudo codes are summarized in Algorithm 2.
Algorithm 2 Search for the optimal regularization parameter λ* Begin: ) is set for λ, and l = min = 1, r = max = q.2: While λ min < λ max , mid = (l + r)/2 , go to step 3, else go to step 8. 3: Randomly divide the predicted values of the KRG, RBF, and SVR surrogate models of the N + N add samples into k (we use k = 5 in this paper) equal parts.4: The Ŷtrain matrix is made up by the predicted values of three individual surrogate models in the k − 1 group, and by singular value decomposition (SVD), which can be expressed as where Ŷtrain = [ Ŷ1 , Ŷ2 , Ŷ3 ], U = (U 1 , U 2 , U 3 ) is a m t × 3 orthogonal matrix, and m t is the number of the k−1 parts of the samples.
) is a 3 × 3 orthogonal matrix.5: After the SVD, ŵe is calculated for λ l and λ r by: ŵe = Vdiag( ŵ1 is the initial weight factors for λ l , and ŵ2 is the initial weight factors for λ r .6: Calculate the weight factors with Equations ( 15) and ( 16) for ŵ1 and ŵ2 , separately, and construct the f e (λ l ) and f e (λ r ) with Equation (1).7: Calculate the RMSE of f e (λ l ) and f e (λ r ) with Equation ( 20), and the current optimal λ c values are calculated as: Back to 2 8: The optimal λ* is equal to the λ c after iteration, and it can be used to construct the RLS-EM.Output: λ*.

Case Studies
In this section, we compare the performance of the RLS-EM with that of the individual models KRG, RBF, and SVR (the detailed construction can be seen in Appendix A) and the ensemble models BP, PWS, NPWS, and OWS described in Section 2. Three types of error metrics were used to evaluate the performances of different surrogate models: root mean squared error (RMSE), which provides a global error measure over the design space; average absolute error (AAE), which ensures that the positive and negative errors will not counteract; and the coefficient of determination (R 2 ), which is a statistical measure of how close the data are to the fitted regression line.
where y is the mean of the observed responses, y i denotes the observed response for x i , ŷi denotes the corresponding prediction, and N t is the number of evaluation points.
We implement the RLS-EM with MATLAB routines, the KRG model was based on a design and analysis of computer experiment toolbox named DACE [48], the RBF model was developed by Sarra [49], and the SVR model was based on the LIBSVM, a library for support vector machines, which was developed by Chang and Lin [50].Four ensemble models including BP, PWS, NPWS, and OWS were implemented in the MATLAB toolbox developed by Viana [51].The cases have been executed with MATLAB R2018a on a computer Intel (R) Core (TM) i7-8700K, CPU @3.7 GHz, 32.0 Gb RAM, 64 bits, and Windows 10.

Numerical Examples
Six numerical examples varying from two-dimensional (2-D) to 12-dimensional (12-D) [42,44] were chosen to test the performance of RLS-EM: LHS was used to generate the training and testing sets, the MATLAB routine "lhsdesign" with "maximin" criterion and 100 iterations were used to generate the (N + N add ) samples and N t tests.The summary of the sampling in the numerical cases is provided in Table 1.Table 2 lists the setup details of the individual models, which were used to develop the ensemble model based on different variable dimensions and nonlinearities.Each individual model has significant differences between variables with different dimensions and different degrees of nonlinearity, e.g., for the low-dimensional variables such as variable with numbers two and four, constant regression can satisfy the accuracy requirements of a KRG model, while the high-dimensional variables require quadratic regression to obtain a more accurate model.Similarly, the kernel parameters and regularization parameters of different dimensional variables with different degrees of nonlinearity are different for the SVR model.Thus, the KRG, RBF, and SVR model setting information for different dimensional variables are listed in detail, as shown in Table 2.
To validate the performance of the different surrogate models, 100 runs were executed for each of the numerical examples.The MATLAB routine "boxplot" was used for easy visualization and comparison.The three-dimensional surface plots of the Branin-Hoo and Camelback functions are shown in Figures 1 and 2, respectively.The nine surface plots in Figure 1 show that each surrogate model fits the Branin-Hoo function well.However, Figure 2 shows that the different surrogate models have considerable differences in the Camelback function fitting.Despite the two functions being highly nonlinear, the RLS-EM can accurately approximate the actual functions.The boxplots of RMSE, AAE, and R 2 for the different test functions are shown in Figures 3-5; the mean and standard deviations of the different surrogate models for the performances are listed in Table 3.After 100 runs were executed, the mean and standard deviation of the RMSE, AAE, and R 2 metrics for the numerical examples are shown in Table 3.For each metric of the numerical examples, the values to the left of the symbol "/" are the mean of the different models, and the values below are the standard deviations corresponding to the models.For the RMSE and AAE metrics, the smaller mean values indicate the better model accuracy, and the smaller standard deviation values show the better robustness.The R 2 metric with a mean value is closer to one and a smaller standard deviation indicate a more accurate and more robust model.  1 The above parameters were set according to experience and can be fine-tuned by many optimization algorithms, which is not covered in this study.

KRG RBF SVR
Linear regression, Gaussian correlation, θ0 = ndv (1/9) , 0.01 < θi < 20 Gaussian basis functions, kernel parameter γ = 1, polynomial term = 1 Gaussian kernel γ = 0.5, regularization parameter C = 100, quadratic loss ε = 0.001 12 KRG RBF SVR Quadratic regression, Gaussian correlation, θ0 = ndv (1/12) , 0.01 < θi < 20 Gaussian basis functions, kernel parameter γ = 2, polynomial term = 1 Gaussian kernel γ = 0.5, regularization parameter C = 100, quadratic loss ε = 0.0001 To validate the performance of the different surrogate models, 100 runs were executed for each of the numerical examples.The MATLAB routine "boxplot" was used for easy visualization and comparison.The three-dimensional surface plots of the Branin-Hoo and Camelback functions are shown in Figures 1 and 2, respectively.The nine surface plots in Figure 1 show that each surrogate model fits the Branin-Hoo function well.However, Figure 2 shows that the different surrogate models have considerable differences in the Camelback function fitting.Despite the two functions being highly nonlinear, the RLS-EM can accurately approximate the actual functions.The boxplots of RMSE, AAE, and R 2 for the different test functions are shown in Figures 3-5; the mean and standard deviations of the different surrogate models for the performances are listed in Table 3.After 100 runs were executed, the mean and standard deviation of the RMSE, AAE, and R 2 metrics for the numerical examples are shown in Table 3.    From Table 3, and Figures 3-5, we can see that no individual surrogate model is always accurate for all test cases, KRG fits the Branin-Hoo function well, while RBF shows a better fitting precision than KRG for the Camelback function.The superiority of ensemble models is not evident for low-dimensional variables functions; however, as the variable dimension and the degree of nonlinearity increase, the ensemble models perform better than most of the individual surrogate models.RLS-EM outperforms all the models in most of the error metrics for the six numerical problems.It shows good fitting performance lower RMSE values on the Camelback, Hartmann-3, Hartmann-6, and Dixon-Price functions.The RMSE and AAE values in Table 3 and their boxplots in Figures 3-5 also show that the RLS-EM is robust.than KRG for the Camelback function.The superiority of ensemble models is not evident for low-dimensional variables functions; however, as the variable dimension and the degree of nonlinearity increase, the ensemble models perform better than most of the individual surrogate models.RLS-EM outperforms all the models in most of the error metrics for the six numerical problems.It shows good fitting performance and lower RMSE values on the Camelback, Hartmann-3, Hartmann-6, and Dixon-Price functions.The RMSE and AAE values in Table 3 and their boxplots in  The BP, PWS, NPWS, and OWS use GMSE as error metric, and they require more computation time than the individual surrogate models, especially for high-dimensional problems.The GMSE error metric takes a relatively longer time to repeatedly construct the individual surrogate models, and the computation time is also affected by the number of divisions.In the RLS-EM, the individual surrogate models are constructed based on the initial samples, the weight factors are obtained by the regularization least squares method, which helps avoid the time spent on repetitively constructing the individual surrogate models.Figure 6 shows the computational cost of Hartmann-6, Extended-Rosenbrock, and Dixon-Price problems, which are represented by the subscript numbers of 1, 2 and 3, respectively.Further, Figure 6 shows that, as the variable dimension increases, BP, PWS, NPWS, and OWS are considerably more time-consuming than RLS-EM.The BP, PWS, NPWS, and OWS use GMSE as error metric, and they require more computation time than the individual surrogate models, especially for high-dimensional problems.The GMSE error metric takes a relatively longer time to repeatedly construct the individual surrogate models, and the computation time is also affected by the number of divisions.In the RLS-EM, the individual surrogate models are constructed based on the initial samples, the weight factors are obtained by the regularization least squares method, which helps avoid the time spent on repetitively constructing the individual surrogate models.Figure 6 shows the computational cost of Hartmann-6, Extended-Rosenbrock, and Dixon-Price problems, which are represented by the subscript numbers of 1, 2 and 3, respectively.Further, Figure 6 shows that, as the variable dimension increases, BP, PWS, NPWS, and OWS are considerably more time-consuming than RLS-EM.

Deformation Prediction for the CNC Milling Machine Bed
A CNC milling machine is mainly composed of a bed, column, slider, and toolbox among other components.The column and slider, under static conditions, exert a large force on the bed, which is expressed by the red arrows in Figure 7.When the milling machine is being operated, the bed is also  The BP, PWS, NPWS, and OWS use GMSE as error metric, and they require more computation time than the individual surrogate models, especially for high-dimensional problems.The GMSE error metric takes a relatively longer time to repeatedly construct the individual surrogate models, and the computation time is also affected by the number of divisions.In the RLS-EM, the individual surrogate models are constructed based on the initial samples, the weight factors are obtained by the regularization least squares method, which helps avoid the time spent on repetitively constructing the individual surrogate models.Figure 6 shows the computational cost of Hartmann-6, Extended-Rosenbrock, and Dixon-Price problems, which are represented by the subscript numbers of 1, 2 and 3, respectively.Further, Figure 6 shows that, as the variable dimension increases, BP, PWS, NPWS, and OWS are considerably more time-consuming than RLS-EM.

Deformation Prediction for the CNC Milling Machine Bed
A CNC milling machine is mainly composed of a bed, column, slider, and toolbox among other components.The column and slider, under static conditions, exert a large force on the bed, which is expressed by the red arrows in Figure 7.When the milling machine is being operated, the bed is also

Deformation Prediction for the CNC Milling Machine Bed
A CNC milling machine is mainly composed of a bed, column, slider, and toolbox among other components.The column and slider, under static conditions, exert a large force on the bed, which is expressed by the red arrows in Figure 7.When the milling machine is being operated, the bed is also affected by the milling impact from the toolbox.Since the deformation has a great influence on the milling precision, the design of the milling machine bed needs good resistance to the deformation; thus, it is very important to accurately predict the deformation during design.
As the milling force is small, we only considered the column and slider weights applied to the milling machine bed, and we predicted the static deformation.The simplified structure of the bed is mainly controlled by eight variables, which are shown in Figure 8.The variables' design space is set as As the milling force is small, we only considered the column and slider weights applied to the milling machine bed, and we predicted the static deformation.The simplified structure of the bed is mainly controlled by eight variables, which are shown in Figure 8   Owing to the heavy computation time, the data set for the milling machine bed design was fixed; thus, it was not possible to generate 100 different designs of experiments cyclically.To solve this problem, in each of the 100 runs, the points for the data sets (N, Nadd, Nt) were chosen randomly at respective ratios from the 200 sampling points.The results of the test are listed in Table 4.

Table 4.
Comparison for the design of milling machine bed.thus, it is very important to accurately predict the deformation during design.
As the milling force is small, we only considered the column and slider weights applied to the milling machine bed, and we predicted the static deformation.The simplified structure of the bed is mainly controlled by eight variables, which are shown in Figure 8   Owing to the heavy computation time, the data set for the milling machine bed design was fixed; thus, it was not possible to generate 100 different designs of experiments cyclically.To solve this problem, in each of the 100 runs, the points for the data sets (N, Nadd, Nt) were chosen randomly at respective ratios from the 200 sampling points.The results of the test are listed in Table 4.

Table 4.
Comparison for the design of milling machine bed.Owing to the heavy computation time, the data set for the milling machine bed design was fixed; thus, it was not possible to generate 100 different designs of experiments cyclically.To solve this problem, in each of the 100 runs, the points for the data sets (N, N add , N t ) were chosen randomly at respective ratios from the 200 sampling points.The results of the test are listed in Table 4. From Table 4, RLS-EM has the best RMSE and R 2 values; further, it has the second-best AAE value.The performance of BP is better than that of KRG, SVR, and other ensembles in AAE; however, the performances of RMSEs of OWS, NPWS, and OWS are better than those of the individual surrogate models.The results reveal that because the linear or nonlinear relationships inside are unknown, when encountering a black-box engineering problem, using an individual surrogate model to approximate the relationship between the design variables and the responses may yield inaccurate results.However, the inaccuracies of each individual surrogate model do not considerably affect the approximate accuracy of RLS-EM; the RLS-EM performs well for the deformation prediction of the milling machine bed.When all the samples are obtained by a time-consuming FEA analysis process, there is no significant increase in the amount of computation caused by the search of the regularization parameter compared to the time-consuming error metrics used in other ensemble methods.Therefore, RLS-EM is an effective engineering problem modeling method that effectively improves the computational efficiency while keeping the modeling accuracy.

Results
In this work, a new method that combines the advantages of least squares method, the regularization, and the augmentation is developed to construct a better and time-saving ensemble model in the cases that only a small number of sample points are available.The weight factors are calculated by the least squares method with the regularization strategy and the augmentation strategy.The augmentation strategy helps to obtain the augmented samples in the unexplored regions by a sample exploration method.On one hand, it helps to improve the accuracy of the individual surrogate models; on the other hand, the augmentation strategy helps to reduce the collinearity problem caused by the intrinsic properties of KRG and RBF and the approximate prediction values on some densely distributed regions.The regularization strategy with an optimal search method to find the best regularization parameter helps to further reduce the collinearity and avoid the potential overfitting problem.
Six numerical functions and a 9-D CNC milling machine bed deformation prediction problem were used to test the proposed RLS-EM method.Four other ensemble models and KRG, RBF, and SVR were adopted for comparison with RLS-EM.The results show that for the numerical functions, the RLS-EM model can provide satisfactory robustness and accuracy, with better or equivalent levels compared to other ensemble methods, while saving a considerable amount of computational cost.The results of the CNC milling machine bed deformation prediction problem also show that the RLS-EM has a good accuracy and robust performance.
In the future work, the hyperparametric optimization will be studied in the RLS-EM, which will further help improve the accuracy and robustness of the RLS-EM, and RLS-EM-based optimization will be studied too.

Algorithm 1 5 :
the KRG, RBF, and SVR prediction values at x i , respectively, the corresponding actual function values Y = [y 1 , y 2 , . . ., y N ] T .The augmentation strategy is implemented to add additional samples in the exploration regions that are far from the N original samples.The number of N add points (the set is X add ) is obtained from Algorithm 1. Pseudo code of augmentation strategy for adding samples Input: X = [ x 1 , x 2 , . . ., x N ].1: Set S pri add empty, S = X.2: Obtain 3 × N add samples by LHS, put them in X lhs .3: For i = 1: 3N add do 4: Calculate the distance of all the members in X lhs to the samples in S. Move the sample with the largest distance from X lhs to S pri add and S. 6: End for 7: Construct KRG, RBF by S, calculate the uncertainties with (13) at the sample set S pri add , storage the difference values in P kr .8: Sort P kr from the largest to the least, choose the top N add values of the corresponding samples, and put them into X add .Output: X add .

Figure 2 .
Figure 2. Surface plots of the Camelback function.Figure 2. Surface plots of the Camelback function.

Figure 3 .
Figure 3. Boxplots of RMSE for the six numerical examples.Figure 3. Boxplots of RMSE for the six numerical examples.

Figure 3 .
Figure 3. Boxplots of RMSE for the six numerical examples.Figure 3. Boxplots of RMSE for the six numerical examples.

Figure 4 .
Figure 4. Boxplots of AAE for the six numerical examples.

Figure 5 .
Figure 5. Boxplots of R 2 for the six numerical examples.

Figure 6 .
Figure 6.Comparison of computational time of ensemble models.

Figure 5 .
Figure 5. Boxplots of R 2 for the six numerical examples.

Figure 5 .
Figure 5. Boxplots of R 2 for the six numerical examples.

Figure 6 .
Figure 6.Comparison of computational time of ensemble models.

Figure 6 .
Figure 6.Comparison of computational time of ensemble models.
x 1 ∈ [40,60], x 2 ∈ [40,60], x 3 ∈ [50,80], x 4 ∈ [40,60], x 5 ∈ [20,40], x 6 ∈ [20,45], x 7 ∈ [15,30], x 8 ∈ [50,80], and x 9 ∈ [40,60]; all the variables are in millimeters.The force of the beam is 56.5 kN and that of the slider is 23.68 kN.The slider is positioned at the initial position of the bed.FEA simulations were carried out to obtain the sample set and the corresponding deformation values.An RLS-EM model was constructed to evaluate the deformation of the bed under the two forces, which are based on the variables with different size values.A total number of 200 sample points were selected for the construction and verification of the proposed ensemble model.thus, it is very important to accurately predict the deformation during design.
. The variables' design space is set as x1 [40,60], x2 [40,60], x3 [50,80], x4 [40,60], x5 [20,40], x6 [20,45], x7 [15,30], x8 [50,80], and x9 [40,60]; all the variables are in millimeters.The force of the beam is 56.5 kN and that of the slider is 23.68 kN.The slider is positioned at the initial position of the bed.FEA simulations were carried out to obtain the sample set and the corresponding deformation values.An RLS-EM model was constructed to evaluate the deformation of the bed under the two forces, which are based on the variables with different size values.A total number of 200 sample points were selected for the construction and verification of the proposed ensemble model.

Figure 7 .
Figure 7. Milling machine tool and the sketch of the bed.

3 Figure 8 .
Figure 8. Sectional sketch for the design variables.

Figure 7 .
Figure 7. Milling machine tool and the sketch of the bed.
. The variables' design space is set as x1 [40,60], x2 [40,60], x3 [50,80], x4 [40,60], x5 [20,40], x6 [20,45], x7 [15,30], x8 [50,80], and x9 [40,60]; all the variables are in millimeters.The force of the beam is 56.5 kN and that of the slider is 23.68 kN.The slider is positioned at the initial position of the bed.FEA simulations were carried out to obtain the sample set and the corresponding deformation values.An RLS-EM model was constructed to evaluate the deformation of the bed under the two forces, which are based on the variables with different size values.A total number of 200 sample points were selected for the construction and verification of the proposed ensemble model.

Figure 7 .
Figure 7. Milling machine tool and the sketch of the bed.

3 Figure 8 .
Figure 8. Sectional sketch for the design variables.

Figure 8 .
Figure 8. Sectional sketch for the design variables.

Table 1 .
Summary specifications for numerical cases.

Table 4 .
Comparison for the design of milling machine bed.