High-Precision Kriging Modeling Method Based on Hybrid Sampling Criteria

: Finding new valuable sampling points and making these points better distributed in the design space is the key to determining the approximate effect of Kriging. To this end, a high-precision Kriging modeling method based on hybrid sampling criteria (HKM-HS) is proposed to solve this problem. In the HKM-HS method, two inﬁlling sampling strategies based on MSE (Mean Square Error) are optimized to obtain new candidate points. By maximizing MSE (MMSE) of Kriging model, it can generate the ﬁrst candidate point that is likely to appear in a sparse area. To avoid the ill-conditioned correlation matrix caused by the too close distance between any two sampling points, the MC (MSE and Correlation function) criterion formed by combining the MSE and the correlation function through multiplication and division is minimized to generate the second candidate point. Furthermore, a new screening method is used to select the ﬁnal expensive evaluation point from the two candidate points. Finally, the test results of sixteen benchmark functions and a house heating case show that the HKM-HS method can effectively enhance the modeling accuracy and stability of Kriging in contrast with other approximate modeling methods.


Introduction
In the field of engineering design and scientific practice, simulation models are widely used to imitate and analyze expensive black box problems. The rapid development and advancement of technology makes some black box problems have greater nonlinearity and computational complexity. Therefore, the simulation models will undoubtedly become more complex and the computing cost will be more expensive. All of these bring additional challenges to the black box simulation problems. Under these circumstances, approximate models, especially for the surrogate models [1,2] (also known as meta-models or response surfaces), are widely applied to better simulate complex black box problems. Using surrogate models to approximately replace expensive simulations can not only significantly reduce the time spent, but also effectively improve the modeling accuracy.
Due to the actual demand of engineering and science, different types of surrogate models have sprung up. The classics mainly include Support Vector Regression (SVR) [3,4], Polynomial Response Surface (PRS) [5,6], Radial Basis Function (RBF) [7,8], Multiple Adaptive Regression Spline (MARS) [9] and Kriging [10][11][12]. The PRS model can solve loworder, low-dimensional and easy-to-access engineering design problems. However, it is not suitable for solving high-dimensional nonlinear multi-modal design problems. RBF is an interpolation method that uses the sum of the weights of the basic functions to approximate complex black box problems. The RBF method can more accurately approximate nonlinear and high-dimensional black box problems, but the approximate requires it to use more function evaluations, which in turn increases the time cost. The computational complexity of the SVR model depends on the number of support vectors rather than the dimension of the sample space. In this way, the dimension problem can be avoided. However, the This strategy assumes that every unobserved point in the region may be the next update point. During the sampling process, the expected improvement function is transformed into a cumulative probability distribution function, and an appropriate new sample is intelligently drawn within a certain probability region. In addition, Li proposed a Krigingbased multi-point sequential sampling optimization (KMSSO) method [28]. So far, there are few relevant literatures and research results to solve this problem. The reasons for this phenomenon are as follows. When adding multiple points at the same time, the correlation between each point needs to be considered, otherwise it may lead to an ill-conditioned matrix, or find the same point and cause modeling failure. Adding multiple points at the same time may take too much time, resulting in low modeling efficiency. Adding multiple points at the same time does not guarantee that each point is promising and may fail to improve the accuracy.
To overcome this bottleneck problem, a high-precision Kriging modeling method based on hybrid sampling criteria (HKM-HS) is proposed. Two sampling strategies, MMSE and MC, are included in HKM-HS. First, the MMSE sampling strategy is used to generate the first candidate point. If the second candidate point is directly generated by MMSE, the distance between the two candidate points may be very close, resulting in the illconditioned correlation matrix. Therefore, the correlation function is added as a selection condition when selecting the second candidate point. MMSE is to ensure that the candidate points have the maximum amount of information, while the minimum correlation function is to ensure that the distance between the candidate points and the existing sample points is not too close, thereby ensuring that the new correlation matrix will not be ill-conditioned. The MC criterion is obtained by combining the maximum mean square error and the minimum correlation function through the multiplication and division method, and the second candidate point is obtained by MC sampling criterion. Finally, the evaluation points are selected from two candidate points.
The rest of this work is organized as follows. The second part discusses the basic background of the Kriging model. The third part proposes two sampling criteria and gives specific steps for determining new sampling points. The fourth part tests 16 benchmark functions and two actual cases and gives the test results. The fifth part is the application of HKM-HS method in house heating. The sixth part is the conclusion.

Kriging
In statistics, the interpolation process of Kriging is finished by a Gaussian process controlled by the priori covariance. Under appropriate a priori assumptions, the Kriging model can provide the best linear unbiased prediction. The specific details are stated as follows.
Suppose there are m sample points X = [x 1 , . . . , x m ] T , X ∈ m×n and object response Y = [y 1 , . . . , y m ] T . The Kriging model is expressed as follows.
The Equation consists of two parts. The first part µ(x) is the trend function. When µ(x) is equal to 0, µ and ∑ p i=1 β i f i (x) are called simple Kriging, ordinary Kriging, and standard Kriging models, respectively. Where f i (x) and β i are the regression function and the corresponding regression function coefficients, respectively. Standard Kriging uses the regression function to determine the change of process mean [29].
The second part of Equation (1) is a stochastic process model established by observing the data and quantifying the correlation of data. z(x) is considered as the realization of the stochastic process Z(x) with mean of 0 and variance of σ 2 . The covariance between sampling points is shown in Equation (2). where process variance σ 2 is a scalar and θ is a key parameter of the kernel function. By optimizing θ, the correlation between the design points can be adjusted adaptively. The expression of the kernel function is as follows.
There are many options for spatial correlation functions. However, the most widely used correlation function is the Gaussian function. The expression of the Gaussian model is as follows.
Under the above assumptions, the best linear unbiased estimate of y(x) iŝ r and R in Equation (5) are Among them, the least square estimationβ of β can be obtained by calculating Equation (8).β In addition, the predicted mean square errorŝ(x) ofŷ(x) can be calculated by Equation (9).
It can be seen from Equations (3), (8) and (9) that the matrix R, vectorβ, and the estimated varianceŝ(x) of the unknown point all depend on the value of the parameter θ. Based on the maximum likelihood estimation theory, the unconstrained optimization algorithm is used to maximize Equation (11) to obtain the optimal θ value [30].
In fact, the best approximation curve is not necessarily obtained from the optimal value, as long as the value is close to the optimal value, there will be a good approximation result.

AME
The AME (adaptive maximum entropy) [21] sampling method determines the new sampling point by maximizing the determinant of the adjusted covariance matrix. The adjusted covariance matrix expression is as follows. where η i and η j are the adjust factors for points x i and x j , respectively. They are determined by the errors at x i and x j . The parameter γ can be adjusted to η i to balance local exploitation and global exploration. The adjusted covariance matrix not only considers the distance factor between two points, but also considers the error of the two points. Sampling more points in the area with a larger error can improve the accuracy of the meta-model. When x i is a known sample point,

HKM-HS Method
The purpose of this paper is to add one or two new sampling points in an optimization iteration of the Kriging model. The HKM-HS method with two infilling sampling criteria can achieve this purpose. The specific details of the two infilling sampling criteria are as follows.

Maximizing Mean Square Error (MMSE)
A good sampling criterion can not only improve the accuracy and efficiency of modeling, but also reduce the cost of modeling. The criterion for evaluating the merits and demerits of the sampling methods are as follows. Whether the newly discovered sampling points have more information in the design field and whether the sampling points can be more evenly distributed in the entire design space under the premise of maintaining certain independence has yet to be determined. The common sampling criteria are maximum mean square error (MMSE), integrated mean square error (IMSE) [22], an adaptive Bayesian sequential sampling method [21] and maximum entropy sampling method (ME) [24]. The ME method tends to add candidate points far away from the current sample points, without considering the objective information of the design point. The adaptive Bayesian sequential sampling method considers not only the distance factor but also the error factor when adding design points, but the running time is slightly longer than the MMSE and IMSE methods. Both MMSE and IMSE methods can obtain the next sampling point by combining the known sample data with the Kriging model. The ideas of the two methods are similar, but the MMSE method is simpler and easier. The MMSE method can obtain sampling points with potentially useful information from the current Kriging model. In the sampling process, the distance between sampling points is considered to ensure that the new sampling data is evenly distributed in the entire design space. Therefore, MMSE sampling criterion is selected as the final sampling criterion. The established model is as follows.
max f 1 (x) =ŝ(x) (13) whereŝ is a function about the unknown point x. The purpose of Equation (13) is to find the candidate point x when the mean square error is maximized.

MC Criterion
The spatial correlation function (SCF) can affect the smoothness of the Kriging model. This correlation applies to not only the known sample points, but also to the quantified observations. There are eight commonly used correlation function models in SCF [31][32][33]. See Table 1 for details.

Spatial Correlation Function Model
The following is an example to prove the advantages of the Gaussian correlation function. A simple one-dimensional problem is expressed as [25]: Different correlation functions are used in Kriging modeling process. Figure 1 shows the trend and error of the approximate function.  Figure 1b is the absolute error between the actual value and the predicted value. The Figure 1a shows the relationship between the original function and the seven Kriging approximate functions based on different correlation functions. The approximate results show that the selection of correlation function does have a certain influence on the accuracy of the Kriging modeling and the selection of new samples. When the kernel function of the Kriging is the Gaussian function, the approximate function is closest to the original function. Figure 1b indicates that the Kriging approximation function with Gaussian kernel function has the smallest error, which shows that the Kriging approximation function has the highest prediction accuracy.

AE (absolute error) in
To further verify the effectiveness of the Gaussian kernel function, a three-dimensional function and a five-dimensional function are used as examples to verify it. Figure 2a,b are graphs of the average absolute error of each experiment run 20 times. Obviously, the Kriging model based on the Gaussian kernel function has the smallest average absolute error, which proves the effectiveness of the Gaussian kernel function. When comparing the effectiveness of multiple kernel functions, the Exponential Gaussian kernel function does not participate in the comparison. There are two main reasons. One is that the specific expression varies according to the parameter θ n+1 , so there is no fixed expression. The other is that the Exponential kernel function and the Gaussian kernel function are both special forms of the Exponential Gaussian kernel function; that is, the exponential kernel and the Gaussian kernel belong to the Exponential Gaussian kernel.
According to these analyses, the Gaussian kernel function is finally selected as the correlation function. Its specific expression is shown in Equation (4).
The smaller the correlation between the new sampling point and the known sampling points, the better. If the correlation value is large, the matrix R may be ill-conditioned. In view of this, the requirement of minimizing the correlation function in Equation (13) need be met to minimize the correlation between new sampling points and all design points.
where parameter r i = R(θ, x, x i ) i = 1, . . . , m is the correlation between new sampling point and the ith sample point. Multiplication and division can be used to deal with multi-objective optimization problems [34,35]. The multiplication or division of multiple objective optimization functions can make multi-objective optimization problems have a clearer meaning and can transform the multi-objective optimization problem into a single-objective optimization problem, so that the single-objective contains the optimization intent of each single objective in the multi-objective. If multiple objective functions need to be maximized, we can multiply them together to form a new objective function f (x) = f 1 (x)· f 2 (x) · · · f k (x); or we can use a weighted geometric mean to construct a new objective function, such as If some of the multiple objective functions need to be maximized and others need to be minimized, the product of all objective functions to be maximized can be divided by the product of all objective functions that are minimized and then maximized, or the product of all the objective functions that need to be minimized is divided by the product of all the objective functions to be maximized and then minimized, and so on. In this way, the problem can be transformed into a problem of maximizing or minimizing the value in the form of "multiplication and division." For the problems raised in this paper, we can get the following single-objective optimization model: where f 1 (x) represents the mean square error function, as shown in Equation (13); f 2 (x) represents the correlation function, as shown in Equation (14). Equation (15) is obtained by combining the two objective functions of Equations (13) and (14) through multiplication and division. The significance of building model (15) is to transform two objective functions into a single objective function. The range of each The range of f 1 cannot be determined because the MSE of different functions is different, but the value of f 1 is always greater than 0.
In addition, the following situations may exist. When the range of f 1 is much smaller than that of f 2 , a small change of the function f 1 will greatly affect the optimization results. In this situation, the function f 1 plays a leading role, and the selection of point pays more attention to the correlation between candidate point and all sample points. Similarly, when the range of f 1 is much bigger than that of f 2 , the function f 2 dominates the optimization results, and the selection of point pays more attention to the amount of information of the candidate points. In either case, the point found can be used as a candidate point. Therefore, Equation (15) can be optimized directly. The optimal value obtained at this time is the second candidate sampling point.

Screening Method
Two candidate sampling points were generated through two infilling sampling criteria. In this paper, the PSO (particle swarm optimization) algorithm of MATLAB was used to optimize the sampling criteria. The first candidate point x i was obtained by optimizing the MMSE function (Equation (13)) through the PSO algorithm. The second candidate point x j was obtained by optimizing the G(x) function (Equation (15)) with the PSO algorithm. Then, one or two candidate sampling points were selected as new sampling points. The method of selecting sampling points is described as follows.
First, set the correlation threshold to 0.001. The range of correlation is [0, 1]. The closer the correlation value is to 0, the higher the mutual independence between the two points. Similarly, the closer the correlation value is to 1, the lower the mutual independence between the two sample points. In addition, the distance between two points with low independence may be too close, resulting in an ill-conditioned correlation matrix. Set the correlation threshold to 0.001, which is a value close enough to 0. This setting can ensure that the two candidate points have independent information and avoids the appearance of an ill-conditioned correlation matrix, which may lead to the consequences of modeling failure.
Second, determine the correlation value between two candidate points. The correlation value between two candidates x i and x j is determined by Gaussian correlation function.
If the correlation value is less than 0.001, it means that the correlation between the two candidates x i and x j is very small. At this time, both x i and x j are selected as new sampling points. If the correlation value is equal to or greater than 0.001, it means that the correlation between the two candidates x i and x j is large. At this time, if both x i and x j are selected as new sampling points, the correlation matrix may be ill-conditioned. Therefore, in this case, select the point with large mean square error as the new sampling point.

Implementation of the HKM-HS Method
This section introduces the specific implementation process of the HKM-HS model. Table 2 shows the specific implementation steps of the HKM-HS model. Table 3 is the pseudocode of the HM-HS algorithm. Figure 3 is the flowchart of the HKM-HS model.

HKM-HS Method
Step 1. Generate initial design points. The LHD (Latin Hypercube Design) method is used to generate m = 2n initial design point x i (i= 1, . . . , m). Then, set the initial sample set to X.
Step 2. Determine the sample sets X and Y. Estimate the expensive objective function of each initial sample point x i (i= 1, . . . , m) to obtain y(x i ). Then, set Y to be the set of all y(x i ). For one or two update points obtained by the infilling sampling criterion, their expensive function evaluation values can be determined. Then the update point and its function evaluation value are added to the sample sets X and Y, respectively. If there is only one new sampling point, then m = m + 1; if there are two new sampling points, then m = m + 2.
Step 3. Build Kriging model. The Kriging model will be established by the new data sets X and Y formed by Step 2, and the construction of Kriging is realized by DACE toolbox in MATLAB.
Step 4. Infill sampling criteria. The new infilling sampling method has two criteria, one is MMSE and the other is the MC sampling strategy based on mean square error and correlation function. See Sections 3.1 and 3.2 for details. The PSO algorithm of MATLAB is used to optimize the sampling criteria. The first candidate point x i is obtained by optimizing the MMSE function (Equation (12)) through PSO algorithm. The second candidate point x j is obtained by optimizing the G(x) function (Equation (14)) with PSO algorithm.
Step 5. Determine new sampling point. First, two candidate sampling points x i and x j are generated with the new infilling sampling criteria. Then, one or two candidate sampling points are selected as new sampling points. See Section 3.3 for details.
Step 6. Stopping criterion. The maximum number of function evaluations is set as N max = 20n. Determine the relationship between m and N max . If the number of sampling points m is greater than N max , stop adding points and go to step 7. If the condition is not met, go back to step 2.
Step 7. Stop. This is the output global approximate Kriging model.

Numerical Experiment
To test the modeling accuracy and stability of the HKM-HS method, the 16 benchmark functions were tested with HKM-HS, MMSE, LHD and AME [21] methods, respectively, and the test results of the four methods were compared. Then the HKM-HS algorithm was applied to the problems of spring design, simple pendulum system and a house heating case.

Benchmark Function Test
In this paper, a total of 16 benchmark functions were tested, and the dimensional range of these functions was 2 to 16 dimensions. Among them, there were 6 two-dimensional benchmark functions. Descriptions of these functions can be found in Table 4. Table 4. Benchmark functions.

Dimension Expression
These benchmark functions were tested by the HKM-HS, MMSE, LHD and AME algorithms, respectively. The number of design points to build these four models is the same (step 6), and the size of initial design point was set as 2n (step 1). There are two main motivations. (a) There is little difference between the initial point setting in the following papers and that in this paper. In the papers [36,37], 2n sample points were selected as the initial points in some experiments. This shows that it is appropriate to select 2n as the number of initial sample points. (b) The accuracy of the initial model is relatively low when the size of the initial design point is fixed to 2n. Experiments performed under this condition can better show the improvement of model accuracy, which can verify the effectiveness of the proposed method. In the modeling process, new sample points are obtained (step 4-5) and the Kriging model is updated constantly (step 2-3).
The six two-dimensional functions were classified according to their similarity in significant physical properties and shapes, so as to more accurately determine the influence of the algorithm on different types of functions. The characteristic of the first type of function is that it has many local minimums. Bukin and Schwefel functions belong to the first type. The results of the MSE obtained by the test are shown in Figures 4 and 5.  The characteristic of the second type of benchmark function is that it is the plateshaped function. The Mccormick function belongs to the second type. Figure 6 indicates the test results of MSE. The third type of benchmark function is a valley-shaped function. The sixHump function belongs to the third type. The results of its MSE test are shown in Figure 7. The remaining two functions are divided into the fourth type. Figures 8 and 9 are their MSE results.    The mean square error results of different types of two-dimensional benchmark functions tested by the HKM-HS, MMSE, LHD and AME algorithms are shown in Figures 4-9. These figures show, based on the test results of these six benchmark functions, that the HKM-HS algorithm has the smallest mean square error, which indicates that the HKM-HS algorithm has good test results for different types of two-dimensional benchmark functions. Compared with other three modeling methods, the HKM-HS modeling method has higher accuracy.
There are ten high-dimensional benchmark functions. The benchmark functions are divided into four types. The characteristic of the first type of function is that it has many local minima. The Levy function and Rastrigin function belong to the first type of function. The property of the second type of function is that it is the plate-shaped function. The DixonPrice and Rosenbrock functions belong to the second type. The third type of function has a steep ridge. The Michalewicz function and Michalewicz 10 function belong to the third type of function. The remaining four functions are divided into the fourth type. Each benchmark function is tested 30 times. The average value and standard deviation of the 30 test results is taken as the data in Table 5. The RMSE results in Table 5 are obtained by leaving a cross-validation method. The calculation formula of RSME is as follows.
where k is the maximum expensive evaluation number. Assuming that x i is one of the sample points, the variance at x i is evaluated asŝ 2 i .ŝ 2 i is the predicted variance obtained from the Kriging model constructed by the remaining k − 1 sample data.      The RMSE results of 10 benchmark functions with different dimensions tested using HKM-HS, AME, MMSE and LHD algorithms are listed in Table 5. Table 5 concluded that in the test of 10 benchmark functions, the test result of 8 functions showed that the HKM-HS algorithm has the smallest average value of RMSE. For the Levy function, the LHD method has the smallest average value and standard deviation of RMSE. The AME algorithm has a better result in the DixonPrice function. The test results of these four methods were compared, and the comparison results show that the accuracy and stability of the HKM-HS method is slightly higher than the other three methods.
The Levy and DixonPrice benchmark functions belong to the first and second types of functions, instead of focusing on one type of function. This shows that the method proposed in this paper is applicable to various types of functions. Figures 10-14 also show that the results obtained by the HKM-HS algorithm have better stability and accuracy. In short, the test results of 10 benchmark functions show that compared with AME, MMSE and LHD, the HKM-HS method has higher modeling accuracy and stability. Therefore, it is feasible to apply the HKM-HS method to practical cases. Figure 15 is a schematic diagram of the spring design. Arora was the first to propose the problem of spring design [38]. The purpose of designing the spring is to make the weight of the spring as small as possible. The minimum deflection, shear stress, and oscillation frequency are used as the constraints of the model.

Spring Design Problem
proposed in this paper is applicable to various types of functions. Figures 10-14 show that the results obtained by the HKM-HS algorithm have better stability and a racy. In short, the test results of 10 benchmark functions show that compared with A MMSE and LHD, the HKM-HS method has higher modeling accuracy and stabi Therefore, it is feasible to apply the HKM-HS method to practical cases. Figure 15 is a schematic diagram of the spring design. Arora was the first to prop the problem of spring design [38]. The purpose of designing the spring is to make weight of the spring as small as possible. The minimum deflection, shear stress, and cillation frequency are used as the constraints of the model. The model expression formula is:

Spring Design Problem
The constraints are: d P P D Figure 15. Spring design.
The model expression formula is: The constraints are: The HKM-HS method studies unconstrained problems, so it is necessary to transform the constraint problem into an unconstrained problem. If any variable in the variable interval satisfies all constraints, then the problem in the interval is an unconstrained problem. The nature of the constraints determines that the three design variables are closely related and influence each other, so it is impossible to obtain an exact variable interval. Therefore, we only need to find the subinterval which satisfies all the constraints. By simplifying and calculating the constraints, the subinterval that satisfies all the constraints obtained in this paper is x 1 ∈ [−3, −2.6], x 2 ∈ −5 × 10 4 , −4 × 10 4 and x 3 ∈ −1.1 × 10 −7 , −9.1 × 10 −8 .
HKM-HS, MMSE, LHD and AME methods are used to model the spring design problem. The modeling results are shown in Table 6 and Figure 16.  The average value and standard deviation of RMSE obtained by running each method 30 times are listed in Table 6. In Table 6, the average RMSE value obtained by the HKM-HS method is the smallest. This shows that in the spring design problem, the accuracy of modeling with the HKM-HS method is higher. Figure 16 and the standard deviation of RMSE reflect the higher stability of the HKM-HS method. Figure 17 is a schematic diagram of a ten-bar planar truss. Barthelemyy was the first to study the ten-bar planar truss [39,40]. The purpose of the study is to minimize the weight of the truss system under the condition that each steel bar stress σ i meets the constraint. The ten design variables x i (i = 1, . . . , 10) are the cross-sectional areas of each steel bar. The problem is expressed as

Ten-Bar Planar Truss Problem
subject to 0.645cm 2 ≤ x i ≤ 64.5cm 2 i = 1, · · · , 10 − 172375 kpa ≤ σ i ≤ 172375 kpa i = 1, · · · , 10 Table 7 shows the RMSE results of 30 runs of each method. Figure 18 is the box plot of RMSE. From Table 7 and Figure 18, it can be seen that compared with AME, MMSE and LHD, the HKM-HS has the smallest average value and standard deviation of RMSE.
This shows that the HKM-HS method has the highest modeling accuracy and stability. Therefore, compared with the other three methods, the HKM-HS method is more suitable for solving this problem.

House Heating
With the continuous progress of society, the problem of an energy shortage cannot be ignored. Energy conservation has received widespread attention worldwide, and China is no exception. Every winter, house heating is a way to consume energy that cannot be underestimated. In order to reduce energy consumption and reduce pollution to the external environment as much as possible, it is necessary to discuss the issue of house heating [41,42]. There are many factors that affect the heating cost of a house, such as the size of the house, the thermal properties of the house materials, the thermal resistance of the house, the characteristics of the heater, the cost of electricity, and the indoor and outdoor temperatures. Different influencing factors have different degrees of influence on the heating cost of houses. Therefore, predicting the heating cost under different conditions has great practical significance for reducing energy consumption, reducing economic costs and reducing environmental pollution. Figure 19 is a schematic diagram of house heating. In order to study the relationship between house heating cost and influencing factors, we need to set some influencing factors as variables. Define the length, width and height of the house as variables x 1 , x 2 and x 3 , respectively. Suppose the number of windows in the house is six and set the length and width of the windows as variables x 4 and x 5 . Set the initial room temperature (outdoor temperature) as variable x 6 . Set the hot air temperature of the heater as variable x 7 . Set the thickness of the glass wool on the wall as variable x 8 . In addition, it is assumed that the flow rate of the heater is 3600 kg/hr, the electricity cost is 0.09 $/kWhr and the constant pressure specific heat capacity of air is 1005.4 J/kg. The house heating simulation model was established by MATLAB and SIMULINK.
The variables x 1 , · · · , x 8 are used as the input parameters of the simulation, and the heating cost of the house is used as the output variable. The heating cost of the house is related to time. Here we take the day as the unit and the heating cost of the house per day as the output variable. The HKM-HS method was used to determine the new valuation point, that is, to determine the different values of 8 variables (influencing factors), and then estimate the heating cost of a day through simulation. From this, the degree of influence of different influencing factors on heating costs can be explored.
In the process of using the HKM-HS method to achieve house heating, first randomly select 8 initial sample points x i (i = 1, · · · , 8) through the LHD method and simulate the initial sample points to obtain y i (i = 1, · · · , 8). Secondly, select new estimation points by the method proposed in this paper and perform simulation until the total number of sample points is 24. In order to ensure the accuracy of the experiment, the average RMSE value of the ten results is taken as the error. RMSE is obtained by the cross-validation method. Table 8 shows the average value and standard deviation of ten RMSE values obtained by applying HKM-HS, MMSE, LHD and AME methods to house heating. The error results show that the modeling accuracy of the HKM-HS method is very good. Figures 20-22 is an analysis of the housing heating problem.    Figure 22, the changing trend of x 6 and x 7 is almost the same as that of the objective value y. This shows that the hot air temperature of the heater and the initial room temperature of the house have a great influence on the heating cost of the house. This may be because the temperature setting of the heater and initial room temperature determine the use time and start-up frequency of the heater. The higher the temperature setting of the heater, the greater the electricity consumption, which will lead to increased electricity cost. On the contrary, the lower the initial room temperature of the house, the higher the electricity cost. (2) In Figure 20, the changing trends of x 1 , x 2 , and x 3 have some similarities with the changing trends of the objective value y, and the similarity of x 1 is more obvious. This indicates that the geometric structure of the house also has a certain influence on the heating cost. The length of the house has the most obvious influence on the heating cost, because the length occupies the largest weight of the house size. The larger the house, the more electricity the heater consumes for heating, and the higher the heating cost.
(3) It can be seen in Figure 21 that the changing trends of x 4 , x 5 , and x 8 are less consistent with that of the objective value y. This implies that these factors have little influence on the heating cost of houses. In short, the HKM-HS method can predict the heating cost of influencing factors under different conditions. Setting a suitable heater temperature according to the size of the house and the initial room temperature can effectively reduce heating cost.

Conclusions
The paper proposes a high-precision Kriging modeling method based on hybrid sampling criteria (HKM-HS). First, LHD was used to select the initial sampling point in this method. Secondly, MMSE and MC criteria were optimized to find two candidate points, which were further screened to obtain the final sampling points. Then the HKM-HS modeling method was used to test 16 benchmark functions and three engineering cases. The test results show that the HKM-HS is almost better than other three Kriging-based modeling methods in terms of accuracy and stability. Furthermore, as an actual simulation problem, a housing heating case was used by the proposed method to show that the HKM-HS has a certain degree of engineering application.
However, it is not appropriate for this method to deal with higher dimensional problems. It is the main reason why the increase in dimensions not only causes the construction cost of the Kriging model to rise rapidly, but also cannot guarantee the final modeling accuracy. In view of this, the exploration of high-dimensional problems based on the Kriging model has attracted the attention of some researchers. The key point is how to achieve dimensionality reduction of the high-dimensional Kriging model while ensuring certain accuracy requirements. At present, there are two main ideas to deal with this key issue. One is to transform high-dimensional modeling problems into lowdimensional modeling problems through effective dimensionality reduction methods, thereby increasing the modeling efficiency of Kriging; the other core idea is to reduce the amount of calculation of hyper-parameters in the Kriging modeling process because optimizing hyper-parameters is a very time-consuming task. These research ideas and directions are yet to be further explored and developed by researchers.