Abstract
Co-Kriging (CK) modeling provides an efficient way to predict responses of complicated engineering problems based on a set of sample data obtained by methods with varying degree of accuracy and computation cost. In this work, the Gaussian random process (GRP) is introduced to construct a novel combination CK model (CK-GRP) to improve the prediction accuracy of the conventional CK model, in which all the sample information provided by different correlation models is well utilized. The features of the new model are demonstrated and evaluated for a numerical case and an engineering application. It is shown that the CK-GRP model proposed in this work is effective and can be used to improve the prediction accuracy and robustness of the CK model.
1. Introduction
Co-Kriging (CK) model is widely used in the optimization design for complex products [1, 2], for its advantage on reducing the computation cost in high-dimensional and strong nonlinear problem. In order to improve the prediction power of the CK model, the ordinary CK model needs to be improved, and the accuracy of the optimization design will be increased as well. Thus far, the majority of studies on improving the fitting performance of the CK model have focused on the following three aspects.
The first objective is to improve the prediction accuracy of the CK model. For example, Han et al. improved the computation accuracy and efficiency of the variable-fidelity surrogate model, combining the gradient-enhanced Kriging and a generalized hybrid bridge function [3]. To solve the variance inflation problem of the collocated CK model, a novel merged secondary variable approach is developed to construct the intrinsic collocated CK model [4, 5]; the prediction accuracy of the collocated CKSM is improved. In order to solve the problems of insufficient samples or too much computation expense of obtaining samples, Wang et al. constructed the indicator CK model, in which the training sample is extracted from the coarse samples using deconvolution process. And the prediction accuracy and efficiency of the indicator CK model are both improved [6].
The second objective is to improve the prediction efficiency of the CK model. Rumpfkeil et al. presented a kind of enhanced Kriging and CK model based on gradient and Hessian information; the computation efficiency of the CK model is improved [7]. Laurent et al. improved the modeling efficiency of the gradient enhanced CK model using a novel multiparametric strategy based on the Latin hyper cubic sampling method [8, 9]. For the purpose of solving the problem of the covariance, cross covariance cannot be calculated in the modeling of a CK model; Furrer and Genton proposed an aggregation CK model, which uses linear unbiased predictor to obtain the variances, and the fitting efficiency of the CK model is enhanced [10]. To improve the computation efficiency of the analysis of static, free vibration, and buckling for laminated composite plates, the moving Kriging (MK) interpolation method with Gaussian correlation model is chosen to interpolate kinematic variables, but the MK model is the original one and is not improved [11–13].
The third objective is to improve the fitting robustness of the CK model. Hoef J. M. et al. constructed a more flexible correlation model (CM) using moving average function, and the Fast Fourier Transform method is applied to solve the problems of large numbers of parameters and difficulties with integration. The simulation results demonstrate that the approach could improve the fitting robustness of the CK model [14]. Pardo-Igúzquiza et al. adopted convolution and deconvolution method to compute the CM and cross-CM, and the fitting robustness of the CK model is improved [15]. Han et al. established a kind of variable-fidelity CK model, and the effectiveness and fitting robustness are verified by concrete aerodynamics cases [16].
The precious studies improved the fitting accuracy, fitting efficiency, and the robustness of fitting results of the CK model using different approaches. The effectiveness of the corresponding modified CK model is demonstrated by different engineering applications or numerical simulations. However, the influence of CM on fitting performance is not considered by most of the researches, and the differences of the sample information obtained by different CM are not realized, too. The standard model of the CM is still determined by artificial selection, and the sample information is provided by the selected “optimal” CM, which restrict the fitting accuracy of a CK model; the optimization accuracy of the CK model-based complex product optimization cannot be guaranteed, too.
Compared with other CK models, this work proposes a combination CK model based on Gaussian random process (CK-GRP), in which the Gaussian random process (GRP) model is adopted as a multiple model fusion method to make full use of sample information obtained by different CMs. The effectiveness of the proposed CK-GRP model is demonstrated using one numerical case and an engineering application. It is shown that the proposed model has the advantages of higher fitting accuracy and strong robustness compared with other CK models in precious studies.
The rest of this work is organized as follows. In Section 2, the basic form of a CK model and the fundament types of CMs are presented. In Section 3, the GRP method is introduced to construct the CK-GRP, and the construction of the CK-GRP model is proposed. In Section 4, the effectiveness of the presented CK-GRP model is verified by two examples: one is a numerical case and the other is an engineering application. In Section 4, a brief conclusion of this work is provided.
2. CK Model and CMs
2.1. Basic Form of a CK Model
For an dimensional problem, suppose both the high accuracy function (expensive to evaluate) and low accuracy function (cheap to evaluate) are used to predict the response in any unknown point . Two datasets are obtained at the sampling sites:the corresponding response iswhere and are the number of sampling sites for the high and low accuracy prediction function, respectively, and in general.
The prediction value of the CK model iswhere and are the weighted factors of the high and low accuracy prediction function, respectively. Suppose there are 2 stationary random processes as shown in the following equation corresponding to and :
Then, the covariance and cross covariance of the random variables in different sites of the design space arewhere and are the process variances of and and the prediction value of the CK model could be deduced in the similar approach to the Kriging model. The prediction value of the CK model is (see [17] for the specific derivation)where
2.2. Fundament Types of CMs
CM is an important function to describe the spatial correlations of sample variables, and the prediction accuracy of the CK model is affected by it. In the modeling of a conventional CK model, the spatial correlation information of sample data is obtained using the selected CM based on experiences, and the semivariable function is calculated. The most commonly used CMs include exponential function (Exp), Gaussian function (Gaussian), spline function (Spline), linear function (Linear), and spherical function (Spherical). The specific forms of these CMs are illustrated in Table 1, in which represents the distance between any two points and and is correlation parameter.
Five CK models could be established based on one sample, according to the five kinds of CMs illustrated in Table 1. In the conventional modeling process of a CK model, a most suitable CM is selected as a variation function to extract sample information in general; then a CK model with best fitting accuracy is constructed. In this modeling process, the useful information provided by other CMs is ignored, which results in a loss of sample information, and the fitting accuracy of the CK model is relatively low. Therefore, in order to improve the accuracy of the CK model, the information provided by all the CMs needs to be extracted to establish a CK model.
3. Combination CK Model (CK-GRP) Based on Gaussian Random Process (GRP)
3.1. GRP Method
Gaussian random process has the characteristic that the combination of finite variables obeys multivariate normal distribution; it could be used to simulate the nonlinear characteristic of various responses flexibly [18]. Therefore, the Gaussian random process is adopted to construct the combination model of different CK models.
Assume is a p-dimensional input variable; : is the corresponding response function; then could be represented as a Gaussian random process:where is a mean value function, expressed as ; is a default polynomial function column vectors; is the corresponding unknown regression coefficients column vectors; is covariance function, to describe the spatial covariance of any two input variables and in a GRP, in which is a prior standard deviation, and is the spatial correlation function, which could be calculated using Gaussian correlation function illustrated in where is an unknown regression coefficients column vectors, used to describe the correlation decrement of the corresponding response values and , with the separation of and . Then, , , and together form a hyperparameters space , which could be calculated using maximum likelihood estimation method [19].
For multidimensional response data , in which the input variable is , the response data of every model obeys a multidimensional Gaussian distribution as shown in the following based on (8).
The response value and mean square error (MSE) of any input data are calculated, as shown in the following equations, respectively, after the hyperparameters space of GRP is calculated.
3.2. Construction of the CK-GRP
In order to improve the integrity of sample information for a CK model, the GRP model is introduced to improve and perfect the sample information of a CK model. And the concrete construction of the CK-GRP model is shown in Figure 1.

The basic steps of the construction are as follows:
Five different types of CMs are used to describe one same sample data structure (cheap data and expensive data contained), and five types of sample information are obtained
Five CK models are established, based on five types of single CM correspondingly, according to the sample information
The Latin hypercubic sampling method is used to obtain the input variables in the domain of definition , and the corresponding response values are calculated using five CK models, combined with the expensive data of original sample ; then the training sample of the combination CK model is formed
Taking as a training sample, five CK models based on one single CM correspondingly are fused using GRP method, and the CK-GRP model is established
Test the fitting accuracy of the CK-GRP model, and the expected improvement (EI) algorithm will be used to search a new filling point, if the response value do not satisfy the convergence condition. Then the response values in the new filling point are calculated in a different type of CK model, and a new training sample Ψ is obtained. On this basis, the CK-GRP model is reconstructed, and by testing the fitting accuracy of the CK-GRP model again until the response value satisfies the convergence condition, then the final CK-GRP model is obtained.
In the modeling of the CK-GRP model using GRP method, the weighted-sum method is adopted to construct the correlation model of the real responses and the prediction values of every model.where is the experiment response value (e.g., the calculated value of simulation software); ε is random observation error; is the unknown regression coefficients column vectors, each regression coefficient to one CK model; is the whole residual deviation of and real response .
As the combination CK model is constructed with several CK models, the multiresponse Gaussian process (MRGP) is adopted to modeling the responses of all the CK models:where is a vector of polynomial regression basis function; is the unknown regression coefficient matrix, and ; is the nonspatial covariance; is the spatial correlation function. Therefore, the covariance of any two CK models in the two points of the input space could be represented as . Similarly, the whole residual deviation could be expressed as a Gaussian random process with one single response:where the unknown hyperparameters could be calculated with the formula shown in [19]. According to (14) and (15), the Gaussian random process model of the experiment response value is
The covariance of and the response value of the th CK model are presented in the following from (13):where is a unit column vector; all the values of the elements are 0 except the th element which is 1. Similarly, the covariance of response values of the th CK model and th CK model iswhere is the covariance of the element for the multiresponse Gaussian process. And the hyperparameters could be calculated using maximum likelihood estimation method, which means that the covariance of the response value in the unknown point and other sample point could be calculated in the following:Similarly, the autocovariance of is
According to the superposition principle of the multiresponse Gaussian response, the experiment response value in the prediction point and the sample data d could be described as a multidimensional Gaussian distribution:where
According to (21), taking as the sample data, the prediction value in any unknown point could be calculated in the following:
The MSE of the corresponding prediction is
4. Numerical Examples
4.1. One-Dimensional Numerical Analytical Case
In order to verify the effectiveness and obtain the features of the presented CK-GRP model, a numerical analytical case [17] is adopted to analyze the effect on improving the prediction accuracy of a CK model. The function represents the high accuracy and large computation cost model. The low accuracy and small computation cost model is given by . The locations of the high and low accuracy samples are and, respectively. CK models with different CMs are constructed based on samples , and the corresponding response values. The prediction response values and their MSE are shown in Figure 2.

(a)

(b)
In Figure 2(a), all the CK models are able to pass through the point of the high prediction accuracy model, but there is large difference of the prediction accuracy between different CK models. The CK model based on Gaussian CM model has the highest prediction accuracy, which basically consists of the original function. According to Figure 2(b), to the one-dimensional test function, the fitting accuracy ordering from high to low of the CK models is CK-Gaussian > CK-Exp > CK-Linear > CK-Spherical > CK-Spline. In order to verify the prediction accuracy of the CK-GRP model further, the CK-GRP model is used to predict the test function, and the obtained prediction response results and their MSE values are shown in Figure 3.

(a)

(b)

(c)
From the analysis of Figure 3(a), the CK-GRP model has higher prediction accuracy compared with the CK-Gaussian model, and it almost coincides with the original test function. According to Figure 3(b), the maximum MSE of the CK-Gaussian model is 0.085. Otherwise, the maximum MSE of the CK-GRP model is only 5.231 × 10−9, which means that the prediction accuracy of the CK-GRP model has been improved significantly. The derivative of the MSE is illustrated in Figure 3(c); the maximum derivative value of the CK-Gaussian model MSE is 1.05 × 104, but the one of the CK-GRP model is only 4.51 × 10−4, which indicates that not only the prediction accuracy of the CK-GRP model is improved obviously, but also the prediction robustness is optimized.
The sample information obtained by different CMs needs to be utilized using GPR first in the modeling of the CK-GRP model. This will cause the increase of the computation time of the presented model. In order to analyze the computational cost of the CK-GRP model, the computation time of the CK-Gaussian model on 1000 test sample points is contrasted. The simulation results show that the computation time of the CK-Gaussian and the CK-GRP model are 0.85 seconds and 2.78 seconds separately. It means that the computational cost of the CK-GPR model is higher than the CK-Gaussian model significantly.
4.2. Modeling of Traveling Stability Indicator for a Rail Vehicle Based on CK-GRP
Traveling stability is an important indicator for the dynamic design of a rail vehicle. The full-scale vehicle test will be complicated and time consuming; as there is a strong nonlinear relation between the rail vehicle traveling stability indicator and design parameters, the number of design variables is large as well. However, the computation time will be decreased if the CK model is introduced to predict the values of the traveling stability indicator in different design parameters combination, and the design efficiency will be raised too. Taking the Chinese electric multiple units CRH2 as an object, the prediction accuracy of the CK-GRP model is verified using a full-scale vehicle test data and the rail vehicle dynamical model [20]. The traveling stability analysis model is shown in Figure 4.

The Sperling index is adopted as the traveling stability index of a rail vehicle; the concrete computational formula [21] iswhere is the vertical vibration acceleration (cm/s2); is the vertical vibration frequency (Hz); is the modifying coefficient of frequency.
Suspension system is an important buffering member of a rail vehicle to reduce the impact force from the irregularity of the track, and the primary vertical stiffness kpz, primary vertical damping Cpz, secondary vertical stiffness ksz, and secondary vertical damping Csz are the parameters that influence the traveling stability of a rail vehicle directly. And the secondary vertical stiffness ksz and secondary vertical damping Csz are the most important parameters that influence the traveling stability [22]. Therefore, the parameters of ksz and Csz are selected as the optimization variables of the vertical stability
Suppose the rail vehicle is traveling on the tracks at 200km/h, and the full-scale test and the dynamic simulation method are used to obtain the Sperling index values ; ; in different combination of suspension parameters combination , ; ; ; ; ; ;
The CK-Gaussian model and CK-GRP model are established, respectively, taking S1, S2 and the corresponding traveling stability index as the training sample. Then, the prediction response is obtained as shown in Figures 5(a) and 5(b), and the corresponding contour map is shown in Figures 5(c) and 5(d).

(a)

(b)

(c)

(d)
From the analysis of the prediction response in Figures 5(a) and 5(b), all the sample points are on the response surface, which indicates that both the CK-Gaussian model and the CK-GRP model could predict the Sperling index value of a rail vehicle in different suspension parameters combinations. The location of the maximum and minimum value is not much different, according to the contour map of Figures 5(c) and 5(d), which indicates that the predict trend of the CK-Gaussian model and the CK-GRP model are basically consistent.
In order to obtain the prediction feature of the CK-GRP model further, taking the MSE as the evaluate index and the CK-Gaussian model served as a contrast, the MSE of the prediction results of the CK-Gaussian model and the CK-GRP model is shown in Figures 6(a) and 6(b). According to Figure 6(a), the maximum MSE of the CK-Gaussian model is 6.950 × 10−4 and mean MSE is 3.153 × 10−4; while the maximum MSE of the CK-GRP model is 4.828 × 10−4, mean MSE is 3.611 × 10−4. It indicates that the prediction error of the CK-GRP model is comparatively low, and the effectiveness of the CK-GRP model of improving the fitting accuracy is verified.

(a)

(b)

(c)

(d)
5. Conclusion
In this work, a combination CK model (CK-GRP) to improve the fitting accuracy of the ordinary CK is proposed. This model is obtained by introducing the Gaussian random process to take full advantage of the sample information provided by different CK models, which are established based on one single different correlation model. The proposed CK-GRP model was demonstrated for a numerical analytical case and a rail vehicle example of predicting the traveling stability in different values of suspension parameters combinations. The examples show that the CK-GRP model has higher fitting accuracy and strong robustness compared with other CK models. But the computational cost is higher significantly, too. It means that the proposed CK-GRP model can potentially be applied to efficient rail vehicle dynamics data prediction as well as suspension parameters optimization or any other research areas where the CK model is in use.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The research was funded by RG Petro-Machinery (Group) Co., Ltd. Neither the funding agency nor any outside organization has participated in study design or has any conflicts of interest. RG Petro-Machinery (Group) Co., Ltd., had final approval of the manuscript.
Acknowledgments
The authors gratefully acknowledge the support from the doctoral researchers boosting program of Xi’an Shiyou University and the support from the RG Petro-Machinery (Group) Co., Ltd.