# A Novel Combination Co-Kriging Model Based on Gaussian Random Process

Authors: Huan Xie; Wei Zeng; Hong Song; Wen Sun; Tao Ren
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 Section2, 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 ann dimensional problem, suppose both the high accuracy function y1 (expensive to evaluate) and low accuracy function y2 (cheap to evaluate) are used to predict the response in any unknown point x. Two datasets are obtained at the sampling sites:(1)S1=x1,1,…,x1,m1T∈Rm1×nS2=x2,1,…,x2,m2T∈Rm2×n;the corresponding response is(2)y1=y11,…,y1m1T∈Rm1y2=y21,…,y2m2T∈Rm2,where m1 and m2 are the number of sampling sites for the high and low accuracy prediction function, respectively, and m2>m1 in general.The prediction value of the CK model is(3)y^1x=λTyS=λ1Ty1+λ2Ty2,where λ1T and λ2T 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 y1 and y2:(4)Y1x=β1+Z1xY2x=β2+Z2x.Then, the covariance and cross covariance of the random variables in different sites of the design space are(5)covZ1x1,Z1x1=σ12R11x1,x1covZ2x2,Z2x2=σ22R22x2,x2covZ1x1,Z2x2=σ1σ2R12x1,x2,where σ12 and σ22 are the process variances of Y1(x) and Y2(x) 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)(6)y^1x=φTβ+rTxR-1yS-Fβ,where(7)φ=01,β=β1β2=FTR-1F-1FTR-1yS,r=r1xr2x,R=R11R12R21R22,yS=y1σ1σ2y2,F=1001∈Rm1+m2×2.


### 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 dk represents the distance between any two points xi and xj and θk is correlation parameter.Table 1
Fundamental types of CMs in a CK model.
CMs R k ( θ , d k ) Gaussian exp ⁡ - θ k d k 2 Exp exp ⁡ - θ k d k Spline 1 - 15 ξ k 2 + 30 ξ k 3  0 ≤ ξ k ≤ 0.2 1.25 1 - ξ k 3  0.2 ≤ ξ k ≤ 1 0  ξ k ≥ 1 Linear m a x 0,1 - θ k d k Spherical 1 - 1.5 ξ k + 0.5 ξ k 3 ,   ξ k = m i n 1 , θ k d kFive CK models could be established based on one sample, according to the five kinds of CMs illustrated in Table1. 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.


## 2.1. Basic Form of a CK Model

For ann dimensional problem, suppose both the high accuracy function y1 (expensive to evaluate) and low accuracy function y2 (cheap to evaluate) are used to predict the response in any unknown point x. Two datasets are obtained at the sampling sites:(1)S1=x1,1,…,x1,m1T∈Rm1×nS2=x2,1,…,x2,m2T∈Rm2×n;the corresponding response is(2)y1=y11,…,y1m1T∈Rm1y2=y21,…,y2m2T∈Rm2,where m1 and m2 are the number of sampling sites for the high and low accuracy prediction function, respectively, and m2>m1 in general.The prediction value of the CK model is(3)y^1x=λTyS=λ1Ty1+λ2Ty2,where λ1T and λ2T 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 y1 and y2:(4)Y1x=β1+Z1xY2x=β2+Z2x.Then, the covariance and cross covariance of the random variables in different sites of the design space are(5)covZ1x1,Z1x1=σ12R11x1,x1covZ2x2,Z2x2=σ22R22x2,x2covZ1x1,Z2x2=σ1σ2R12x1,x2,where σ12 and σ22 are the process variances of Y1(x) and Y2(x) 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)(6)y^1x=φTβ+rTxR-1yS-Fβ,where(7)φ=01,β=β1β2=FTR-1F-1FTR-1yS,r=r1xr2x,R=R11R12R21R22,yS=y1σ1σ2y2,F=1001∈Rm1+m2×2.


## 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 dk represents the distance between any two points xi and xj and θk is correlation parameter.Table 1
Fundamental types of CMs in a CK model.
CMs R k ( θ , d k ) Gaussian exp ⁡ - θ k d k 2 Exp exp ⁡ - θ k d k Spline 1 - 15 ξ k 2 + 30 ξ k 3  0 ≤ ξ k ≤ 0.2 1.25 1 - ξ k 3  0.2 ≤ ξ k ≤ 1 0  ξ k ≥ 1 Linear m a x 0,1 - θ k d k Spherical 1 - 1.5 ξ k + 0.5 ξ k 3 ,   ξ k = m i n 1 , θ k d kFive CK models could be established based on one sample, according to the five kinds of CMs illustrated in Table1. 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.Assumex=[x1,x2,…,xp]∈Rp is ap-dimensional input variable; y: Rp→R is the corresponding response function; then y(x) could be represented as a Gaussian random process:(8)yx~gpmx,Vx,x′,where m(x) is a mean value function, expressed as hxTβ; h(x) is a default polynomial function column vectors; β is the corresponding unknown regression coefficients column vectors; V(x,x′)=σ2R(x,x′) is covariance function, to describe the spatial covariance of any two input variables x and x′ in a GRP, in which σ is a prior standard deviation, and Rx,x′ is the spatial correlation function, which could be calculated using Gaussian correlation function illustrated in (9)Rx,x′=exp-∑k=1pωkxk-xk′2,where ω=ω1,ω2,…,ωpT is an unknown regression coefficients column vectors, used to describe the correlation decrement of the corresponding response valuesy(x) and y(x′), with the separation of x and x′. Then, ω, β, and σ2 together form a hyperparameters space Φ=(β,σ2,ω), which could be calculated using maximum likelihood estimation method [19].For multidimensional response dataΨ=yx1,yx2,…,yxnT, in which the input variable is X=x1T,…,xnTT, the response data of every model obeys a multidimensional Gaussian distribution as shown in the following based on (8).(10)yx1⋮yxn~NHβ,σ2R=Nhx1T⋮hxnTβ,σ2Rx1,x1…Rx1,xn⋮⋱⋮Rx1,xn⋯Rxn,xn.The response valuey^(xp) and mean square error (MSE) MSE(y^(xp)) of any input data xp are calculated, as shown in the following equations, respectively, after the hyperparameters space Φ of GRP is calculated.(11)y^xp=hxpTβ+Rxp,XR-1d-Hβ,(12)MSEy^xp=σ21-Rxp,XR-1Rxp,XT+hxpT-ΗR-1Rxp,XTT×ΗR-1H-1×hxpT-ΗR-1Rxp,XT.


### 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 Figure1.Figure 1
Construction of the CK-GRP model.The basic steps of the construction are as follows:( 1 ) 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 φi(x)(i=1,2,…,5) are obtained( 2 ) Five CK models y^i(x,CMi),i=1,2,…,5 are established, based on five types of single CM correspondingly, according to the sample information φi(x)(i=1,2,…,5)( 3 ) The Latin hypercubic sampling method is used to obtain the input variables in the domain of definition X=x1T,…,x5TT, and the corresponding response values are calculated using five CK models, combined with the expensive data of original sample x6T,yx6T; then the training sample Ψ=yx1,yx2,…,yx6T of the combination CK model is formed( 4 ) Taking Ψ=yx1,yx2,…,yx6T as a training sample, five CK models based on one single CM correspondingly are fused using GRP method, and the CK-GRP model y(x)=gph·Tβ,σ2R·,· is established( 5 ) 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 responsesy(x) and the prediction values ym(x) of every model.(13)yx=yex-ε=ymxρ+δx,where ye(x) is the experiment response value (e.g., the calculated value of simulation software); ε is random observation error; ρ=[ρ1,ρ2,…,ρi] is the unknown regression coefficients column vectors, each regression coefficient to one CK model; δ(x) is the whole residual deviation of ym(x)ρ and real response y(x).As the combination CK model is constructed with several CK models, the multiresponse Gaussian process (MRGP) is adopted to modeling the responsesym(x) of all the CK models:(14)ymx~gphm·TBm,∑mRm·,·,where hm(x) is a vector of polynomial regression basis function; Bm=βm1,βm2,…,βmi is the unknown regression coefficient matrix, and βmi=[β1m[i],β2m[i],…,βnm[i]]; ∑m is the nonspatial covariance; Rmx,x′ is the spatial correlation function. Therefore, the covariance of any two CK models in the two points x,x′ of the input space could be represented as covymix,ymjx′=∑i,jmRmx,x′). Similarly, the whole residual deviation δ(x) could be expressed as a Gaussian random process with one single response:(15)δ·~gphδ·Tβδ,σδ2Rδ·,·,where the unknown hyperparameters Θ1δ={βδ,σδ,ωδ} could be calculated with the formula shown in [19]. According to (14) and (15), the Gaussian random process model ye(x) of the experiment response value is(16)ye·~gphm·TBmρ+hδ·Tβδ,∑mρRm·,·+σδ2Rδ·,·.The covariance ofye(x) and the response value of the ith CK model are presented in the following from (13):(17)covyex,ymix′=covymxρ+σx+ε,ymx′ei12=covymxρ+σx+ε,ymx′ei=ρTei∑mRmx,x′,where ei is a unit column vector; all the values of the elements are 0 except the ith element which is 1. Similarly, the covariance of response values of the ith CK model and jth CK model is(18)covymix,ymjx′=covymxei,ymx′ej=eiT∑mejRmx,x′,where eiT∑mej is the covariance ∑m of the (i,j) 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 ye(xq) in the unknown point xq and other sample point ye(x) could be calculated in the following:(19)covyexq,yex=covymxqρ+δxq+ε,ymxρ+δx+ε=ρT∑mρRmxq,x+σδ2Rδxq,x.Similarly, the autocovariance of ye(xq) is(20)covyexq,yexq=covymxqρ+δxq+ε,ymxqρ+δxq+ε=ρT∑mρRmxq,xq+σδ2Rδxq,xq.According to the superposition principle of the multiresponse Gaussian response, the experiment response valueye(xq) in the prediction point and the sample datad could be described as a multidimensional Gaussian distribution:(21)yexqd~gpHqHβ,TqVqVdTqT,where(22)Hq=ρ1HmxqT⋯ρiHmxqT∣HδxqT=ρT⊗HmxqT∣HδxqT,(23)Tq=ρT∑me1Rmxq,xm1⋯ρT∑meiRmxq,xm1∣ρT∑mρRmxq,xe+σδ2Rδxq,xeT,(24)Vq=ρT∑mρRmxq,xq+σδ2Rδxq,xq.According to (21), taking d as the sample data, the prediction value ye(xq) in any unknown point xq could be calculated in the following:(25)y^exq∣d=Hqβ+TqVd-1d-Hβ.The MSE of the corresponding prediction is(26)MSEy^exq∣d=Vq-TqVd-1TqT+HqT-HTVd-1TqTTωHqT-HTVd-1TqT.


## 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.Assumex=[x1,x2,…,xp]∈Rp is ap-dimensional input variable; y: Rp→R is the corresponding response function; then y(x) could be represented as a Gaussian random process:(8)yx~gpmx,Vx,x′,where m(x) is a mean value function, expressed as hxTβ; h(x) is a default polynomial function column vectors; β is the corresponding unknown regression coefficients column vectors; V(x,x′)=σ2R(x,x′) is covariance function, to describe the spatial covariance of any two input variables x and x′ in a GRP, in which σ is a prior standard deviation, and Rx,x′ is the spatial correlation function, which could be calculated using Gaussian correlation function illustrated in (9)Rx,x′=exp-∑k=1pωkxk-xk′2,where ω=ω1,ω2,…,ωpT is an unknown regression coefficients column vectors, used to describe the correlation decrement of the corresponding response valuesy(x) and y(x′), with the separation of x and x′. Then, ω, β, and σ2 together form a hyperparameters space Φ=(β,σ2,ω), which could be calculated using maximum likelihood estimation method [19].For multidimensional response dataΨ=yx1,yx2,…,yxnT, in which the input variable is X=x1T,…,xnTT, the response data of every model obeys a multidimensional Gaussian distribution as shown in the following based on (8).(10)yx1⋮yxn~NHβ,σ2R=Nhx1T⋮hxnTβ,σ2Rx1,x1…Rx1,xn⋮⋱⋮Rx1,xn⋯Rxn,xn.The response valuey^(xp) and mean square error (MSE) MSE(y^(xp)) of any input data xp are calculated, as shown in the following equations, respectively, after the hyperparameters space Φ of GRP is calculated.(11)y^xp=hxpTβ+Rxp,XR-1d-Hβ,(12)MSEy^xp=σ21-Rxp,XR-1Rxp,XT+hxpT-ΗR-1Rxp,XTT×ΗR-1H-1×hxpT-ΗR-1Rxp,XT.


## 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 Figure1.Figure 1
Construction of the CK-GRP model.The basic steps of the construction are as follows:( 1 ) 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 φi(x)(i=1,2,…,5) are obtained( 2 ) Five CK models y^i(x,CMi),i=1,2,…,5 are established, based on five types of single CM correspondingly, according to the sample information φi(x)(i=1,2,…,5)( 3 ) The Latin hypercubic sampling method is used to obtain the input variables in the domain of definition X=x1T,…,x5TT, and the corresponding response values are calculated using five CK models, combined with the expensive data of original sample x6T,yx6T; then the training sample Ψ=yx1,yx2,…,yx6T of the combination CK model is formed( 4 ) Taking Ψ=yx1,yx2,…,yx6T as a training sample, five CK models based on one single CM correspondingly are fused using GRP method, and the CK-GRP model y(x)=gph·Tβ,σ2R·,· is established( 5 ) 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 responsesy(x) and the prediction values ym(x) of every model.(13)yx=yex-ε=ymxρ+δx,where ye(x) is the experiment response value (e.g., the calculated value of simulation software); ε is random observation error; ρ=[ρ1,ρ2,…,ρi] is the unknown regression coefficients column vectors, each regression coefficient to one CK model; δ(x) is the whole residual deviation of ym(x)ρ and real response y(x).As the combination CK model is constructed with several CK models, the multiresponse Gaussian process (MRGP) is adopted to modeling the responsesym(x) of all the CK models:(14)ymx~gphm·TBm,∑mRm·,·,where hm(x) is a vector of polynomial regression basis function; Bm=βm1,βm2,…,βmi is the unknown regression coefficient matrix, and βmi=[β1m[i],β2m[i],…,βnm[i]]; ∑m is the nonspatial covariance; Rmx,x′ is the spatial correlation function. Therefore, the covariance of any two CK models in the two points x,x′ of the input space could be represented as covymix,ymjx′=∑i,jmRmx,x′). Similarly, the whole residual deviation δ(x) could be expressed as a Gaussian random process with one single response:(15)δ·~gphδ·Tβδ,σδ2Rδ·,·,where the unknown hyperparameters Θ1δ={βδ,σδ,ωδ} could be calculated with the formula shown in [19]. According to (14) and (15), the Gaussian random process model ye(x) of the experiment response value is(16)ye·~gphm·TBmρ+hδ·Tβδ,∑mρRm·,·+σδ2Rδ·,·.The covariance ofye(x) and the response value of the ith CK model are presented in the following from (13):(17)covyex,ymix′=covymxρ+σx+ε,ymx′ei12=covymxρ+σx+ε,ymx′ei=ρTei∑mRmx,x′,where ei is a unit column vector; all the values of the elements are 0 except the ith element which is 1. Similarly, the covariance of response values of the ith CK model and jth CK model is(18)covymix,ymjx′=covymxei,ymx′ej=eiT∑mejRmx,x′,where eiT∑mej is the covariance ∑m of the (i,j) 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 ye(xq) in the unknown point xq and other sample point ye(x) could be calculated in the following:(19)covyexq,yex=covymxqρ+δxq+ε,ymxρ+δx+ε=ρT∑mρRmxq,x+σδ2Rδxq,x.Similarly, the autocovariance of ye(xq) is(20)covyexq,yexq=covymxqρ+δxq+ε,ymxqρ+δxq+ε=ρT∑mρRmxq,xq+σδ2Rδxq,xq.According to the superposition principle of the multiresponse Gaussian response, the experiment response valueye(xq) in the prediction point and the sample datad could be described as a multidimensional Gaussian distribution:(21)yexqd~gpHqHβ,TqVqVdTqT,where(22)Hq=ρ1HmxqT⋯ρiHmxqT∣HδxqT=ρT⊗HmxqT∣HδxqT,(23)Tq=ρT∑me1Rmxq,xm1⋯ρT∑meiRmxq,xm1∣ρT∑mρRmxq,xe+σδ2Rδxq,xeT,(24)Vq=ρT∑mρRmxq,xq+σδ2Rδxq,xq.According to (21), taking d as the sample data, the prediction value ye(xq) in any unknown point xq could be calculated in the following:(25)y^exq∣d=Hqβ+TqVd-1d-Hβ.The MSE of the corresponding prediction is(26)MSEy^exq∣d=Vq-TqVd-1TqT+HqT-HTVd-1TqTTωHqT-HTVd-1TqT.


## 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 y1=6x-22sin⁡(12x-4),x∈[0,1] represents the high accuracy and large computation cost model. The low accuracy and small computation cost model is given by y2=0.5y1+10(x-0.5)-5. The x locations of the high and low accuracy samples are S1={0,0.6,1.0} andS2={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}, respectively. CK models with different CMs are constructed based on samples S1, S2 and the corresponding response values. The prediction response values and their MSE are shown in Figure 2.Figure 2
Comparison of different CK models for one-dimensional numerical analytical case: (a) approximation response; (b) approximation mean square error.
(a) (b)In Figure2(a), all the CK models are able to pass through the point y1 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.Figure 3
Comparison of the CK-GRP model and the CK-Gaussian for the numerical analytical case: (a) approximation response; (b) approximation mean square error; (c) derivative of approximation mean square error.
(a) (b) (c)From the analysis of Figure3(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.Figure 4
Traveling stability analysis model of a rail vehicle.The Sperling index is adopted as the traveling stability index of a rail vehicle; the concrete computational formula [21] is(27)W=0.896a3f·Ff10,where a is the vertical vibration acceleration (cm/s2); f is the vertical vibration frequency (Hz); F(f) 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 stiffnesskpz, primary vertical dampingCpz, secondary vertical stiffnessksz, and secondary vertical dampingCsz are the parameters that influence the traveling stability of a rail vehicle directly. And the secondary vertical stiffnessksz and secondary vertical dampingCsz are the most important parameters that influence the traveling stability [22]. Therefore, the parameters ofksz andCsz are selected as the optimization variables X=(x1,x2) of the vertical stabilitySuppose 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 valuesW1={2.151;2.203;2.218}; W2={2.201;2.184;2.176;2.160;2.211;2.215;2.192;2.202;2.171;2.213} in different combination of suspension parameters combination S1={(98.758.8);(106.622.1);(22547.4)}, S2={(82.945.3);(90.851.6); (114.530.1); (122.326.3)  (138.138.9); (146.160);(185.549.5); (201.324.2); (209.228.4); (218.921.7)}.The CK-Gaussian model and CK-GRP model are established, respectively, taking S1, S2 and the corresponding traveling stability index W1,W2 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).Figure 5
Prediction response and contour map of the CK-Gaussian model and the CK-GRP model: (a) prediction response of the CK-Gaussian model; (b) prediction response of the CK-GRP model; (c) contour map of the CK-Gaussian model; (d) contour map of the CK-GRP model.
(a) (b) (c) (d)From the analysis of the prediction response in Figures5(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.Figure 6
Prediction response MSE and its contour map of the CK-Gaussian model and the CK-GRP model: (a) MSE of the CK-Gaussian model; (b) MSE of the CK-GRP model; (c) MSE contour map of the CK-Gaussian model; (d) MSE contour map of the CK-GRP model.
(a) (b) (c) (d)


## 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 y1=6x-22sin⁡(12x-4),x∈[0,1] represents the high accuracy and large computation cost model. The low accuracy and small computation cost model is given by y2=0.5y1+10(x-0.5)-5. The x locations of the high and low accuracy samples are S1={0,0.6,1.0} andS2={0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0}, respectively. CK models with different CMs are constructed based on samples S1, S2 and the corresponding response values. The prediction response values and their MSE are shown in Figure 2.Figure 2
Comparison of different CK models for one-dimensional numerical analytical case: (a) approximation response; (b) approximation mean square error.
(a) (b)In Figure2(a), all the CK models are able to pass through the point y1 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.Figure 3
Comparison of the CK-GRP model and the CK-Gaussian for the numerical analytical case: (a) approximation response; (b) approximation mean square error; (c) derivative of approximation mean square error.
(a) (b) (c)From the analysis of Figure3(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.Figure 4
Traveling stability analysis model of a rail vehicle.The Sperling index is adopted as the traveling stability index of a rail vehicle; the concrete computational formula [21] is(27)W=0.896a3f·Ff10,where a is the vertical vibration acceleration (cm/s2); f is the vertical vibration frequency (Hz); F(f) 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 stiffnesskpz, primary vertical dampingCpz, secondary vertical stiffnessksz, and secondary vertical dampingCsz are the parameters that influence the traveling stability of a rail vehicle directly. And the secondary vertical stiffnessksz and secondary vertical dampingCsz are the most important parameters that influence the traveling stability [22]. Therefore, the parameters ofksz andCsz are selected as the optimization variables X=(x1,x2) of the vertical stabilitySuppose 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 valuesW1={2.151;2.203;2.218}; W2={2.201;2.184;2.176;2.160;2.211;2.215;2.192;2.202;2.171;2.213} in different combination of suspension parameters combination S1={(98.758.8);(106.622.1);(22547.4)}, S2={(82.945.3);(90.851.6); (114.530.1); (122.326.3)  (138.138.9); (146.160);(185.549.5); (201.324.2); (209.228.4); (218.921.7)}.The CK-Gaussian model and CK-GRP model are established, respectively, taking S1, S2 and the corresponding traveling stability index W1,W2 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).Figure 5
Prediction response and contour map of the CK-Gaussian model and the CK-GRP model: (a) prediction response of the CK-Gaussian model; (b) prediction response of the CK-GRP model; (c) contour map of the CK-Gaussian model; (d) contour map of the CK-GRP model.
(a) (b) (c) (d)From the analysis of the prediction response in Figures5(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.Figure 6
Prediction response MSE and its contour map of the CK-Gaussian model and the CK-GRP model: (a) MSE of the CK-Gaussian model; (b) MSE of the CK-GRP model; (c) MSE contour map of the CK-Gaussian model; (d) MSE contour map of the CK-GRP model.
(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.



---

*Source: 6372572-2018-09-03.xml*