Ride comfort optimization of a multi-axle heavy motorized wheel dump truck based on virtual and real prototype experiment integrated Kriging model

Abstract

The optimization of hydro-pneumatic suspension parameters of a multi-axle heavy motorized wheel dump truck is carried out based on virtual and real prototype experiment integrated Kriging model in this article. The root mean square of vertical vibration acceleration, in the center of sprung mass, is assigned as the optimization objective. The constraints are the natural frequency, the working stroke, and the dynamic load of wheels. The suspension structure for the truck is the adjustable hydro-pneumatic suspension with ideal vehicle nonlinear characteristics, integrated with elastic and damping elements. Also, the hydraulic systems of two adjacent hydro-pneumatic suspension are interconnected. Considering the high complexity of the engineering model, a novel kind of meta-model called virtual and real prototype experiment integrated Kriging is proposed in this article. The interpolation principle and the construction of virtual and real prototype experiment integrated Kriging model were elucidated. Being different from traditional Kriging, virtual and real prototype experiment integrated Kriging combines the respective advantages of actual test and Computer Aided Engineering simulation. Based on the virtual and real prototype experiment integrated Kriging model, the optimization results, obtained by experimental verification, showed significant improvement in the ride comfort by 12.48% for front suspension and 11.79% for rear suspension. Compared with traditional Kriging, the optimization effect was improved by 3.05% and 3.38% respectively. Virtual and real prototype experiment integrated Kriging provides an effective way to approach the optimal solution for the optimization of high-complexity engineering problems.

Introduction

The ride comfort of a vehicle is significantly influenced by its suspension system, the road surface roughness, and the vehicle speed. Many research works have been done to optimize the suspension to improve the ride comfort. Jayachandran and Krishnapillai ¹ optimized the values of stiffness and damping coefficient for a passive suspension system and introduced a semi-active control system to improve the ride comfort for a passenger car. Kaldas et al. ² proposed a new optimization methodology for optimizing the damper top mount characteristics to improve the vehicle ride comfort. Eriksson and Friberg ³ calculated the response of the bus by using a finite element model describing the complete vehicle. The body of the bus was condensed using a substructuring technique to reduce the optimization time. Finally, the optimized stiffness and related damping of each power unit rubber insulator were found to achieve the optimization of ride comfort. Yang et al. ⁴ proposed a systematic methodology to evaluate and improve the ride comfort. In the study, a new index called annoyance rate was presented to indicate the quantitative correlation between objective method and subjective comment. But the precondition of further investigation and optimization is an accurate model. He et al. ⁵ established a four-axle heavy truck and road coupling system and carried out a multi-objective optimization based on modified skyhook damping control to improve the ride comfort, while the traditional optimization method could not reduce the computational burden.

Even in the field of construction machinery such as earth moving equipment, ride comfort has also become an important competitive parameter for manufacturers of mobile machines. It means the challenge to improve the ride comfort of the machinery, namely, the components for suspension system should be chosen properly. Often the design would be confined by commercially available components. ⁶ Evaluation and improvement of ride comfort are made efficient by model-based prototyping. ⁷ In order to find the proper parameters, most engineers carried out the optimization based on the computer model^8,9 because it is possible to change parameters and evaluate the performance in a cost-efficient way compared to changing and testing a full-scale prototype. ¹⁰ Just as in the research of Yong Yang et al., ⁴ the high-precision simulation model is needed.

In the practical complex engineering application, the high precision of the model more likely means it will cost more time to make the optimization. Therefore, the computational expense would be increased. In order to alleviate the computational burden of engineering analyses, the Kriging surrogate model has been used for a variety of applications, especially in the design of complex mechanical system and aerospace engineering.

Kriging was inspired by the work of Danie G Krige, a South African geologist, for mining estimation, but he never formalized the method. Then, Georges Matheron ¹¹ developed the theory of regionalized variables based on the work done by Krige and called the method Kriging in his honor. Noel Cressie ¹² summarized the historical origins of Kriging in detail. The standard approach, first applied to deterministic computer simulations by Sacks et al. ¹³ called ordinary Kriging (OK), is a variant of Design and Analysis of Computer Experiments (DACE). Since then, it has been widely used in approximate computation models. Along with the development of Computer Aided Engineering (CAE) technique, the modeling and optimization of Kriging based on CAE analysis have offered perfect performance. Martin and Simpson^14,15 conducted a study on using Kriging models to approximate deterministic computer models and discussed the applicability of various Kriging variants. Duchaine et al. ¹⁶ introduced the Kriging meta-model with a design of experiments (DOE) approach to limit the impact of many evaluations of flow computations which are CPU- and time-consuming. The application of this method was proved to be feasible with available computing power and yields manageable response time. Based on the simulated samples, Lim et al. ¹⁷ developed a surrogate model using Kriging for the fast computation of hybrid sandwich panel in the design optimization. Li et al. ¹⁸ proposed an approach Kriging-assisted multi-objective genetic algorithm (K-MOGA) for multi-objective design optimization. Kriging-based meta-modeling was embedded within a MOGA. Examples were used to illustrate that K-MOGA reduced the number of simulation calls by evaluating some individuals in the population by Kriging meta-models instead of the simulation. In the work of Shinkyu Jeong et al., ¹⁹ the genetic algorithm using the Kriging model for objective function evaluation was introduced. A criterion “expected improvement (EI)” was used for the selection of the additional sample point. This made it possible not only to improve the accuracy of the response surface but also to explore the global optimum efficiently. Yung-Chang Cheng and Cheng-Kang Lee ²⁰ integrated the uniform design and Kriging interpolation, took the critical hunting speed of the vehicle as the optimization function, and created an input–output segregate model based on the Kriging interpolation method to optimize the suspension parameters of a high speed railway vehicle, and obtained an ideal result. Wang Hao et al. ²¹ put forward a modified sequential Kriging optimization to change the twice optimization problem into once by adding more than one point at the same time. In their work, the proposed method could approach the global optimization quickly and accurately in the application of active suspension.

However, in some particularly complex project, especially with strong nonlinearity, the order of Kriging surrogate model obtained from the pure CAE analysis could be lowered. In that case, it would lead to the loss of the extremum or the optimal solution.

In this article, the hydro-pneumatic suspension design and optimization for a multi-axle heavy motorized wheel dump truck (Figure 1) is presented based on approximate model.

Figure 1.

The multi-axle heavy motorized wheel dump truck.

To find the proper design, the actual test results are introduced to modify the meta-model to remedy the above-mentioned shortage. Hence, a novel meta-model called virtual and real prototype experiment integrated Kriging (VRPEI-Kriging) is proposed to perform the ride comfort optimization of the multi-axle heavy motorized wheel dump truck.

Modeling

The vehicle dynamic model

With the continuous development of national economy construction, heavy motorized wheel dump truck has been used in national road construction, mining and other kinds of engineering, and construction of essential transport infrastructure. The cost of the design research and development is huge, and the test cost also accounts for a large proportion in the process of design optimization. To make the design and optimization convenient, the dynamic model is built up in modules using batch files and the model can easily be parameterized, as shown in Figure 2.

Figure 2.

Dynamic model of dump truck.

The suspension structure for the vehicle is adjustable hydro-pneumatic suspension with ideal vehicle nonlinear characteristics, which integrated elastic element and damping element. It can improve not only the ride comfort of the vehicle but also a comparatively perfect regulating function. For instance, no matter in no load or full load, the suspension can keep the vehicle at a constant height from the ground and then make the suspension stiffness adapt to the changes of no load or full load. In addition to these, the adjustable hydro-pneumatic suspension can realize the switch of elastic suspension to adjust the attitude of the vehicle, such as leaning forward and backward, rising and falling, and so on.

There are eight axles in the dump truck. Each axle with the same structure is equipped with a hydro-pneumatic spring. Because the ride comfort of the vehicle is mainly considered here, the steering system is not taken into account. In the suspension system, the rubber bushing between connectors is given approximately limited to the lack of model parameters. The tire response is crucial to evaluate ride comfort realistically on off-road conditions. A tire model developed for an off-road vehicle with big tires is applied in the dynamic model. The tire model needs only few modeling parameters to handle short-wave terrain with obstacles. ²²

The characteristics of the hydraulic system play a pivotal role in the vehicle ride comfort. The hydraulic system of two adjacent hydro-pneumatic suspension is interconnected. In the half vehicle model, as shown in Figure 3, the pressure of the same hydraulic system in I and II is always the same, so does III and IV, as shown in Figure 4. This structure helps keep the traveling vehicle’s body balanced and avoid inhomogeneous stress.

Figure 3.

The vehicle plane diagram.

Figure 4.

Hydro-pneumatic spring structure diagram.

For the vehicle suspension structure, the corresponding simulation model of the hydro-pneumatic spring was established in AMESim, as shown in Figure 5.

Figure 5.

The simulation model of hydro-pneumatic spring.

From the hydraulic simulation model, the variational piston speed was obtained and force changed with the piston displacement. Then, the hydro-pneumatic spring characteristic was introduced into ADAMS to get the complete CAE simulation model, as shown in Figure 6.

Figure 6.

The simulation flow chart in modeFRONTIER.

The working environment of the mine car is bad and the road traffic is poor. Considering the large mass of the vehicle, the mine road was constructed referring to the standard three-dimensional (3D) random road and maintained regularly. The 3D pavement model was built based on sinusoidal superposition principle, as shown in Figure 7.

Figure 7.

The 3D pavement model.

Design variables and constraints

Design variables

For hydro-pneumatic suspension, the initial volume and initial pressure of accumulator, and the diameter of damping orifice are the main factors that affect the stiffness and damping of hydro-pneumatic spring. These parameters were adjusted to achieve the adjustable attitude of the vehicle. For the convenience of adjusting the parameters of hydro-pneumatic spring in actual experiments, based on the known hydraulic cylinder structure (Figure 4) and the mass distribution coefficient, select the initial pressure P0 and initial volume V0 of accumulator, and the diameter of damping orifice d0 in front hydro-pneumatic spring and rear hydro-pneumatic spring as the design variables to be optimized. The range of variables value is shown in Table 1.

Table 1.

Range of variables value.

		Initial	Minimum	Maximum
Front hydro-pneumatic spring	Initial pressure of accumulator P0_f (Pa)	5.5 × 106	3 × 106	1 × 107
	Volume of accumulator V0_f (L)	5	3	10
	Damping orifice d0_f (m)	1.6 × 10−2	1.2 × 10−2	3.5 × 10−2
Rear hydro-pneumatic spring	Initial pressure of accumulator P0_r (Pa)	1 × 107	5 × 106	1.5 × 107
	Volume of accumulator V0_r (L)	4	3	10
	Damping orifice d0_r (m)	9 × 10−3	6 × 10−3	1.2 × 10−2

According to the requirement of the mine car and characteristics of the structure, the design variable values should be constrained by the offset frequency of the front and rear suspension, the working stroke of the hydraulic cylinder, and the dynamic load of the wheel.

Offset frequency

The sprung mass is very large in full load, so the front and rear suspension offset frequency is relatively low. But low offset frequency may reduce the operation stability of the vehicle; therefore, the offset frequency of front and rear suspension should be constrained. Offset frequency is defined as

$$n=\frac{1}{2\pi }\sqrt{\frac{K}{M_b}}$$

(1)

where $n$ stands for the suspension offset frequency, $K$ for the suspension stiffness, and $M_b$ for the suspension sprung mass.

Working stroke

Hydro-pneumatic spring is installed between the frame and the axle, so the working stroke of hydro-pneumatic spring, from no load to full load, is subjected to the conditions of the vehicle structure. Thus, it is necessary to constrain the dynamic displacement of the suspension. The working stroke is defined as

$$S_{k-m}=\frac{(P_m-P_k)V_k}{P_m{A}_1}$$

(2)

where $S_{k-m}$ stands for the amount of compression for hydro-pneumatic spring from no load to full load, $P_m$ for the pressure of the gas chamber in full load, $P_k$ for the pressure of the gas chamber in no load, $V_k$ for the volume of the gas chamber in no load, and $A_1$ for the cross-sectional area of the piston rod.

Dynamic load of wheels

According to the characteristics of hydro-pneumatic spring’s variable stiffness, the stiffness of the suspension is larger when the vehicle is at full load. Therefore, the working stroke and the dynamic load of the wheel are smaller than those of other engineering vehicles. The dynamic load of the wheel impacts the adhesion between the wheel and the road, associated with the grip of the tires, and plays a crucial role in handling stability and driving safety of the vehicle. Thus, the dynamic load of front and rear wheels should be constrained in the optimization. The dynamic load coefficient is defined as

$$\xi =\frac{{\left[\frac{1}{T}\underset{0}{\overset{T}{\int }}{F}_d^2\mathrm{dt}\right]}^{\frac{1}{2}}}{G}$$

(3)

where $\xi$ stands for the dynamic load coefficient, $T$ for the analysis of vibration time, $F_d$ for the dynamic loads of the wheels, and $G$ for the static load of the wheels.

Optimization objective

In the research of vehicle ride comfort, the sprung mass’ centroid of the whole vehicle model includes three directions of vibration acceleration—vertical, longitudinal, and lateral. In this article, the root mean square (RMS) value of vertical vibration acceleration in the center of the sprung mass is taken as the objective function $J$ . The expression is

$$J=a_w={[\underset{0.5}{\overset{80}{\int }}{W}^2(f)G_a(f)\mathrm{df}]}^{\frac{1}{2}}$$

(4)

where $a_w$ stands for the RMS value of vertical vibration acceleration in the center of the sprung mass, $G_a(f)$ for the acceleration power spectral density function, and $W(f)$ for the frequency weighting function which can be expressed as

$$W(f)=\{\begin{array}{cc}\hfill \begin{array}{c}\hfill 0.5\hfill \\ \multicolumn{1}{c}{f/4}\\ \multicolumn{1}{c}{1}\\ \multicolumn{1}{c}{12.5/f}\end{array}\hfill & \hfill \begin{array}{c}\hfill (0.5<f<2)\hfill \\ \multicolumn{1}{c}{(2<f<4)}\\ \multicolumn{1}{c}{(4<f<12.5)}\\ \multicolumn{1}{c}{(12.5<f<80)}\end{array}\hfill \end{array}$$

(5)

In all, the optimization model can be obtained

$$\begin{array}{cc}\hfill \mathrm{find}\hfill & \hfill \mathbf{X}=[P0\_f,V0\_f,d0\_f,P0\_r,V0\_r,d0\_r]\hfill \\ \multicolumn{1}{c}{\mathrm{minimize}}& \hfill \mathbf{F}(\mathbf{X})=J\hfill \\ \multicolumn{1}{c}{\mathrm{subject}\mathrm{to}}& \hfill \{\begin{array}{cc}\multicolumn{1}{c}{}& {\mathbf{X}}_{min}\le \mathbf{X}\le {\mathbf{X}}_{max}\hfill \\ \multicolumn{1}{c}{}& n_{min}\le n_f\le n_{max}\hfill \\ \multicolumn{1}{c}{}& n_{min}\le n_r\le n_{max}\hfill \\ \multicolumn{1}{c}{}& \begin{array}{c}\hfill S_{k-m}\le 300\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\hfill {\xi }_f\le \frac{1}{3},\hfill & \hfill {\xi }_r\le \frac{1}{3}\hfill \end{array}}\end{array}\hfill \end{array}\hfill \end{array}$$

where ${\mathbf{X}}_{min}$ and ${\mathbf{X}}_{max}$ represent the minimum and maximum values of design variables, respectively; $n_f$ and $n_r$ stand for the offset frequency of front and rear suspension, respectively; and ${\xi }_f$ and ${\xi }_r$ for the pre and post tires’ relative dynamic load, respectively.

VRPEI-Kriging meta-model

The VRPEI-Kriging model is based on the VRPEI-Kriging interpolation as the core. In the case of virtual simulation results, it combines the actual test results to construct the approximate meta-model. The main difference with the traditional Kriging model lies in that VRPEI-Kriging considers the analysis results with high reliability from the actual tests, which have cost thousands of dollars and countless hours.

Principle of VRPEI-Kriging interpolation

Take ${\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,\dots ,{\mathbf{x}}_n\in \Omega$ as the known points from the computer, where Ω is the domain of the computer model. The corresponding responses are $\mathbf{Z}=[Z({\mathbf{x}}_1),Z({\mathbf{x}}_2),\dots ,Z({\mathbf{x}}_n)]$ . And m points from the actual test model be ${\mathbf{t}}_1,{\mathbf{t}}_2,\dots ,{\mathbf{t}}_m$ . The resulting outputs are $\mathbf{T}=[T({\mathbf{t}}_1),T({\mathbf{t}}_2),\dots ,T({\mathbf{t}}_m)]$ . Note that ${\mathbf{x}}_i\ne {\mathbf{t}}_k$ . Mingling Z with T as the final response domain, consider a linear predictor of the output; thus, the predictor is

$$Z^{*}(\mathbf{x})={\displaystyle \sum _{i=1}^n}{\alpha }_{i}Z({\mathbf{x}}_i)+{\displaystyle \sum _{k=1}^m}{\beta }_{k}T({\mathbf{t}}_k)$$

(6)

${\alpha }_i$ and ${\beta }_k$ as the weight coefficient, respectively, describe the contribution degree of sample value and experimental value on the estimate value. In statistics

$$E[Z({\mathbf{x}}_i)]=E[T({\mathbf{t}}_k)]=u$$

(7)

In accordance with the unbiased estimation

$$E[Z^{*}(\mathbf{x})]=E[{\displaystyle \sum _{i=1}^n}{\alpha }_{i}Z({\mathbf{x}}_i)+{\displaystyle \sum _{k=1}^m}{\beta }_{k}T({\mathbf{t}}_k)]=u$$

(8)

$${\displaystyle \sum _{i=1}^n}{\alpha }_i+{\displaystyle \sum _{k=1}^m}{\beta }_k=1$$

(9)

Considering the proportion of the number of test points, the environment of the test, and so on, the test weight ratio T_w is introduced

$${\displaystyle \sum _{i=1}^n}{\alpha }_i=1-T_w$$

(10)

$${\displaystyle \sum _{k=1}^m}{\beta }_k=T_w$$

(11)

The RMS of estimation error is

$$\begin{array}{cc}\multicolumn{1}{c}{{\sigma }_E^2}& =D[Z(\mathbf{x})-Z^{*}(\mathbf{x})]=D[Z(\mathbf{x})]\hfill \\ \multicolumn{1}{c}{}& +D[Z^{*}(\mathbf{x})]-2C(Z(\mathbf{x}),Z^{*}(\mathbf{x}))\hfill \\ \multicolumn{1}{c}{}& =C(\mathbf{x},\mathbf{x})+{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{j=1}^n}{\alpha }_i{\alpha }_{j}C({\mathbf{x}}_i,{\mathbf{x}}_j)\hfill \\ \multicolumn{1}{c}{}& +{\displaystyle \sum _{k=1}^m}{\displaystyle \sum _{l=1}^m}{\beta }_k{\beta }_{{}_l}C({\mathbf{t}}_k,{\mathbf{t}}_l)\hfill \\ \multicolumn{1}{c}{}& +2{\displaystyle \sum _{i=1}^n}{\displaystyle \sum _{k=1}^m}{\alpha }_i{\beta }_{k}C({\mathbf{x}}_i,{\mathbf{t}}_k)\hfill \\ \multicolumn{1}{c}{}& -2({\displaystyle \sum _{i=1}^n}{\alpha }_{i}C({\mathbf{x}}_i,\mathbf{x})+{\displaystyle \sum _{k=1}^m}{\beta }_{k}C({\mathbf{t}}_k,\mathbf{x}))\hfill \end{array}$$

(12)

So we want to get $\mathbf{\alpha }=[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n{]}^T$ and $\mathbf{\beta }=[{\beta }_1,{\beta }_2,\dots ,{\beta }_m{]}^T$ to obtain the minimum ${\sigma }_E^2$ .

When the random function does not meet the second order smoothly, and satisfy the intrinsic hypothesis, $\gamma (h)=C(0)-C(h)$ .

Using the Lagrange method, we obtained that

$$\left(\begin{array}{ccc}\hfill \mathbf{A}\hfill & \hfill \mathbf{B}\hfill & \hfill \mathbf{G}\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}^T}& \hfill \mathbf{D}\hfill & \hfill \mathbf{H}\hfill \\ \multicolumn{1}{c}{{\mathbf{G}}^T}& \hfill {\mathbf{H}}^T\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill \mathbf{\alpha }\hfill \\ \multicolumn{1}{c}{\mathbf{\beta }}\\ \multicolumn{1}{c}{\mathbf{\chi }}\end{array}\right)=\left(\begin{array}{c}\hfill {\mathbf{K}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_2}\\ \multicolumn{1}{c}{{\mathbf{K}}_3}\end{array}\right)$$

(13)

where

$$\mathbf{A}={\left(\begin{array}{cccc}\hfill \gamma ({\mathbf{x}}_1,{\mathbf{x}}_1)\hfill & \hfill \gamma ({\mathbf{x}}_1,{\mathbf{x}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{x}}_1,{\mathbf{x}}_n)\hfill \\ \multicolumn{1}{c}{\gamma ({\mathbf{x}}_2,{\mathbf{x}}_1)}& \hfill \gamma ({\mathbf{x}}_2,{\mathbf{x}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{x}}_2,{\mathbf{x}}_n)\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{\gamma ({\mathbf{x}}_n,{\mathbf{x}}_1)}& \hfill \gamma ({\mathbf{x}}_n,{\mathbf{x}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{x}}_n,{\mathbf{x}}_n)\hfill \end{array}\right)}_{n\times n}$$

$$\mathbf{B}={\left(\begin{array}{cccc}\hfill \gamma ({\mathbf{x}}_1,{\mathbf{t}}_1)\hfill & \hfill \gamma ({\mathbf{x}}_1,{\mathbf{t}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{x}}_1,{\mathbf{t}}_m)\hfill \\ \multicolumn{1}{c}{\gamma ({\mathbf{x}}_2,{\mathbf{t}}_1)}& \hfill \gamma ({\mathbf{x}}_2,{\mathbf{t}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{x}}_2,{\mathbf{t}}_m)\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{\gamma ({\mathbf{x}}_n,{\mathbf{t}}_1)}& \hfill \gamma ({\mathbf{x}}_n,{\mathbf{t}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{x}}_n,{\mathbf{t}}_m)\hfill \end{array}\right)}_{n\times m}$$

$$\begin{array}{cc}\multicolumn{1}{c}{}& \mathbf{D}={\left(\begin{array}{cccc}\hfill \gamma ({\mathbf{t}}_1,{\mathbf{t}}_1)\hfill & \hfill \gamma ({\mathbf{t}}_1,{\mathbf{t}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{t}}_1,{\mathbf{t}}_m)\hfill \\ \multicolumn{1}{c}{\gamma ({\mathbf{t}}_2,{\mathbf{t}}_1)}& \hfill \gamma ({\mathbf{t}}_2,{\mathbf{t}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{t}}_2,{\mathbf{t}}_m)\hfill \\ \multicolumn{1}{c}{⋮}& \hfill ⋮\hfill & \hfill \ddots \hfill & \hfill ⋮\hfill \\ \multicolumn{1}{c}{\gamma ({\mathbf{t}}_m,{\mathbf{t}}_1)}& \hfill \gamma ({\mathbf{t}}_m,{\mathbf{t}}_2)\hfill & \hfill \cdots \hfill & \hfill \gamma ({\mathbf{t}}_m,{\mathbf{t}}_m)\hfill \end{array}\right)}_{m\times m},\hfill \\ \multicolumn{1}{c}{}& \mathbf{G}={\left(\begin{array}{cc}\hfill \begin{array}{c}\hfill 1\hfill \\ \multicolumn{1}{c}{1}\\ \multicolumn{1}{c}{⋮}\end{array}\hfill & \hfill \begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{⋮}\end{array}\hfill \\ \multicolumn{1}{c}{1}& \hfill 0\hfill \end{array}\right)}_{n\times 2},\mathbf{H}={\left(\begin{array}{cc}\hfill \begin{array}{c}\hfill 0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{⋮}\end{array}\hfill & \hfill \begin{array}{c}\hfill 1\hfill \\ \multicolumn{1}{c}{1}\\ \multicolumn{1}{c}{⋮}\end{array}\hfill \\ \multicolumn{1}{c}{0}& \hfill 1\hfill \end{array}\right)}_{m\times 2}\hfill \end{array}$$

$${\mathbf{K}}_1={[\gamma ({\mathbf{x}}_1,\mathbf{x}),\gamma ({\mathbf{x}}_2,\mathbf{x}),\dots ,\gamma ({\mathbf{x}}_n,\mathbf{x})]}^T$$

$${\mathbf{K}}_2={[\gamma ({\mathbf{t}}_1,\mathbf{x}),\gamma ({\mathbf{t}}_2,\mathbf{x}),\dots ,\gamma ({\mathbf{t}}_m,\mathbf{x})]}^T$$

$${\mathbf{K}}_3=[1-T_w,T_w{]}^T$$

$$\mathbf{\alpha }=[{\alpha }_1,{\alpha }_2,\dots ,{\alpha }_n{]}^T$$

$$\mathbf{\beta }=[{\beta }_1,{\beta }_2,\dots ,{\beta }_m{]}^T$$

$\mathbf{\chi }=[\epsilon ,\lambda {]}^T$ for the Lagrange factor.

The obtained minimum RMS is expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{{\sigma }_{Emin}^2=}& {\displaystyle \sum _{i=1}^n}{\alpha }_i\gamma ({\mathbf{x}}_i,\mathbf{x})+{\displaystyle \sum _{k=1}^m}{\beta }_k\gamma ({\mathbf{t}}_k,\mathbf{x})\hfill \\ \multicolumn{1}{c}{}& +(1-T_w)\epsilon +T_w·\lambda -\gamma (\mathbf{x},\mathbf{x})\hfill \end{array}$$

(14)

In the VRPEI-Kriging interpolation model, the introduction of test weight ratio T_w is the key to distinguish from other interpolations, such as Simple Kriging (SK), OK, Universal Kriging (UK), CoKriging (CK), Bayesian Kriging (BK), indicator Kriging (IK), and so on. T_w , a complex variable, is determined by the proportion of experimental date, the condition of actual test, the distribution of the test date, and so on.

In the actual construction of the VRPEI-Kriging approximation model, the choice of experimental factor depends on

$$\begin{array}{cc}\hfill \mathrm{find}\hfill & \hfill T_w\in (0,1)\hfill \\ \multicolumn{1}{c}{s.t.}& \hfill \begin{array}{cc}\hfill min\hfill & \hfill {\sigma }_{Emin}^2\hfill \end{array}\hfill \end{array}$$

From (13) and (14), it can be obtained that

$$\begin{array}{cc}\multicolumn{1}{c}{{\sigma }_{Emin}^2}& ={\displaystyle \sum _{i=1}^n}{\alpha }_i\gamma ({\mathbf{x}}_i,\mathbf{x})+{\displaystyle \sum _{k=1}^m}{\beta }_k\gamma ({\mathbf{t}}_k,\mathbf{x})\hfill \\ \multicolumn{1}{c}{}& +(1-T_w)\epsilon +T_w·\lambda -\gamma (\mathbf{x},\mathbf{x})\hfill \\ \multicolumn{1}{c}{}& =\left(\begin{array}{ccc}\hfill {\mathbf{\alpha }}^T\hfill & \hfill {\mathbf{\beta }}^T\hfill & \hfill {\mathbf{\chi }}^T\hfill \end{array}\right)\left(\begin{array}{c}\hfill {\mathbf{K}}_1\hfill \\ \multicolumn{1}{c}{{\mathbf{K}}_2}\\ \multicolumn{1}{c}{{\mathbf{K}}_3}\end{array}\right)-\gamma (\mathbf{x},\mathbf{x})\hfill \\ \multicolumn{1}{c}{}& =\left(\begin{array}{ccc}\hfill {\mathbf{\alpha }}^T\hfill & \hfill {\mathbf{\beta }}^T\hfill & \hfill {\mathbf{\chi }}^T\hfill \end{array}\right)\left(\begin{array}{ccc}\hfill \mathbf{A}\hfill & \hfill \mathbf{B}\hfill & \hfill \mathbf{G}\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}^T}& \hfill \mathbf{D}\hfill & \hfill \mathbf{H}\hfill \\ \multicolumn{1}{c}{{\mathbf{G}}^T}& \hfill {\mathbf{H}}^T\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill \mathbf{\alpha }\hfill \\ \multicolumn{1}{c}{\mathbf{\beta }}\\ \multicolumn{1}{c}{\mathbf{\chi }}\end{array}\right)-\gamma (\mathbf{x},\mathbf{x})\hfill \end{array}$$

(15)

The function is

$$f=\left(\begin{array}{ccc}\hfill {\mathbf{\alpha }}^T\hfill & \hfill {\mathbf{\beta }}^T\hfill & \hfill {\mathbf{\chi }}^T\hfill \end{array}\right)\left(\begin{array}{ccc}\hfill \mathbf{A}\hfill & \hfill \mathbf{B}\hfill & \hfill \mathbf{G}\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}^T}& \hfill \mathbf{D}\hfill & \hfill \mathbf{H}\hfill \\ \multicolumn{1}{c}{{\mathbf{G}}^T}& \hfill {\mathbf{H}}^T\hfill & \hfill 0\hfill \end{array}\right)\left(\begin{array}{c}\hfill \mathbf{\alpha }\hfill \\ \multicolumn{1}{c}{\mathbf{\beta }}\\ \multicolumn{1}{c}{\mathbf{\chi }}\end{array}\right)$$

(16)

The optimum T_w is found to minimize $f$ .

Construction of VRPEI-Kriging approximation model

A VRPEI-Kriging model could be given in equation (17). The model consists of two parts. The first part, $\mathbf{f}{(x)}^T\mathbf{\mu }$ , is a simple linear regression of the data. It indicates a global model of the design space. The second part models the systematic lack-of-fit or deviations from the linear model so that the whole model interpolates the experimental points generated according to one of the designs of experiment approach

$$\hat{y}(\mathbf{x})=\mathbf{f}{(\mathbf{x})}^T\mathbf{\mu }+e(\mathbf{x})$$

(17)

The random process $e(\mathbf{x})$ is assumed to be a Gaussian stationary process with zero mean and covariance

$$V(x_i,x_j)={\sigma }^2\gamma (x_i,x_j)$$

(18)

The spatial correlation function, $\gamma (•)$ , affects the smoothness of the model and differentiability of the response surface.

For the convenience of description, we named the known points from the computer model as ${\mathbf{x}}_1,{\mathbf{x}}_2,{\mathbf{x}}_3,\dots ,{\mathbf{x}}_n$ and m points from the actual test model as ${\mathbf{x}}_{n+1},{\mathbf{x}}_{n+2},{\mathbf{x}}_{n+3},\dots ,{\mathbf{x}}_{n+m}$ ; the resulting outputs are $\mathbf{Y}=[y({\mathbf{x}}_1),y({\mathbf{x}}_2),\dots ,y({\mathbf{x}}_n),y({\mathbf{x}}_{n+1}),y({\mathbf{x}}_{n+2}),\dots ,y({\mathbf{x}}_{n+m})]$ . Thus, at any point x, a linear predictor of the output can be expressed

$$\hat{y}(\mathbf{x})={\mathbf{\lambda }}^T(\mathbf{x})\mathbf{Y}$$

(19)

The VRPEI-Kriging approach treats $\hat{y}(\mathbf{x})$ as a random variable and finds the best linear unbiased predictor, ${\mathbf{\lambda }}^T(\mathbf{x})\mathbf{Y}$ , which minimizes the mean square error (MSE) of the prediction

$$\mathrm{MSE}[\hat{y}(\mathbf{x})]=E[{\mathbf{\lambda }}^T(\mathbf{x})\mathbf{Y}-y(\mathbf{x}{)]}^2$$

(20)

subject to the unbiasedness constraint

$$E[{\mathbf{\lambda }}^T(\mathbf{x})\mathbf{Y}]=E[y(\mathbf{x})]$$

(21)

UK is defined with a set of regression functions

$$\mathbf{f}(x)=[f_1(x),f_2(x),...,f_p(x)]$$

(22)

Define a vector $\mathbf{F}$ that is the value of $\mathbf{f}(x)$ evaluated at each of the known locations

$$\mathbf{F}=[\mathbf{f}(x_1),\mathbf{f}(x_2),...,\mathbf{f}(x_n),\mathbf{f}(x_{n+1}),\mathbf{f}(x_{n+2}),...,\mathbf{f}(x_{n+m}{)]}^T$$

(23)

The correlation matrix, $\mathbf{R}$ , is composed of the spatial correlation function evaluated at each possible combination of known points

$$\mathbf{R}={\left(\begin{array}{cc}\hfill \mathbf{A}\hfill & \hfill \mathbf{B}\hfill \\ \multicolumn{1}{c}{{\mathbf{B}}^T}& \hfill \mathbf{D}\hfill \end{array}\right)}_{(n+m)\times (n+m)}$$

(24)

This matrix, $\mathbf{R}$ , by definition of the spatial correlation function, $\gamma (•)$ , being positive semi-definite, is a positive semi-definite matrix. It is also symmetric since $\gamma ({\mathbf{x}}_1,{\mathbf{x}}_2)=\gamma ({\mathbf{x}}_2,{\mathbf{x}}_1)$ , and the diagonal consists of all ones since $\gamma ({\mathbf{x}}_1,{\mathbf{x}}_1)=1$ . Set a vector that represents the correlation between an unknown point and all known points

$$\begin{array}{c}\multicolumn{1}{c}{\mathbf{r}(\mathbf{x})={[\gamma ({\mathbf{x}}_1,\mathbf{x}),\gamma ({\mathbf{x}}_2,\mathbf{x}),\dots ,\gamma ({\mathbf{x}}_n,\mathbf{x}),\gamma ({\mathbf{x}}_{n+1},\mathbf{x}),\gamma ({\mathbf{x}}_{n+2},\mathbf{x}),\dots ,\gamma ({\mathbf{x}}_{n+m},\mathbf{x})]}^T}\end{array}$$

(25)

If ${\mathbf{\lambda }}^T(\mathbf{x})$ solves this constrained minimization problem, then ${\mathbf{\lambda }}^T(\mathbf{x})\mathbf{Y}$ is called a best linear unbiased predictor for $\hat{y}(\mathbf{x})$ . Solving ${\mathbf{\lambda }}^T(\mathbf{x})$ and substituting into (17), we obtain

$$\hat{y}(\mathbf{x})={\mathbf{f}}^T(\mathbf{x})\hat{\mathbf{\mu }}+{\mathbf{r}}^T(\mathbf{x}){\mathbf{R}}^{-1}(\mathbf{Y}-\mathbf{F}\hat{\mathbf{\mu }})$$

(26)

where the regression coefficient $\hat{\mathbf{\mu }}$ is

$$\hat{\mathbf{\mu }}=({\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{F}{)}^{-1}{\mathbf{F}}^T{\mathbf{R}}^{-1}\mathbf{Y}$$

(27)

The selection of the correlation function is chosen by the user when generating the VRPEI-Kriging model. In the statistical and engineering literature, the Gaussian function ¹² is by far the most popular and is also used in this work. It is defined as

$$\gamma ({\mathbf{x}}_i,{\mathbf{x}}_j)=\exp (-{\displaystyle \sum _{k=1}^p}{\theta }_k{|x_i^k-x_j^k|}^2)$$

(28)

where ${\theta }_k$ $(0\le {\theta }_k\le \infty )$ is the vector of unknown correlation parameters which is of length p, the number of design variables. Also $x_i^k$ , $x_j^k$ are the kth element of the sample points ${\mathbf{x}}_i,{\mathbf{x}}_j$ .

Under the assumption that $e(\mathbf{x})$ in equation (17) follows a Gaussian process, the likelihood function for it is directly related to its probability distribution function given by

$$L[\widehat{\mu },{\sigma }^2,R|Y]=\frac{1}{\sqrt{{(2\pi {\sigma }^2)}^n\left|R\right|}}\exp \left[-\frac{1}{2{\sigma }^2}(Y-F\widehat{\mu })R^{-1}(Y-F\widehat{\mu })\right]$$

(29)

The negative of the Napierian Logarithm (LN)-likelihood is given by

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{LN}=\frac{n}{2}\ln \left[2\pi {\sigma }^2\right]+\frac{1}{2}\ln \left[\left|\mathbf{R}\right|\right]}\\ \multicolumn{1}{c}{+\frac{1}{2{\sigma }^2}(\mathbf{Y}-\mathbf{F}\hat{\mathbf{\mu }}){\mathbf{R}}^{-1}(\mathbf{Y}-\mathbf{F}\hat{\mathbf{\mu }})}\end{array}$$

(30)

From equation (30), the maximum likelihood estimation of ${\sigma }^2$ could be obtained

$${\hat{\sigma }}^2=\frac{(\mathbf{Y}-\mathbf{F}\hat{\mathbf{\mu }}){\mathbf{R}}^{-1}(\mathbf{Y}-\mathbf{F}\hat{\mathbf{\mu }})}{n}$$

(31)

By plugging it into equation (30), we obtain the concentrated log-likelihood function which is only related to the parameter $\mathbf{\theta }$

$$\mathrm{LN}(\mathbf{\theta })=\frac{n}{2}\ln \left[2\pi {\sigma }^2\right]+\frac{1}{2}\ln \left[\left|\mathbf{R}\right|\right]+\frac{n}{2}$$

(32)

So the value of $\mathbf{\theta }$ could be obtained by solving the minimization problem

$$\begin{array}{cc}\multicolumn{1}{c}{\mathrm{find}}& {\theta }_k>0\hfill \\ \multicolumn{1}{c}{\min }& \frac{n}{2}\ln \left[2\pi {\sigma }^2\right]+\frac{1}{2}\ln \left[\left|\mathbf{R}\right|\right]+\frac{n}{2}\hfill \end{array}$$

(33)

Therefore, if $\mathbf{R}$ is known, the maximum likelihood estimates of $\hat{\mathbf{\mu }}$ and ${\sigma }^2$ are easily calculated.

Assessment of VRPEI-Kriging meta-model

Due to the fact that the interpolation model goes exactly through the observations, it is meaningless to assess the capability of the meta-model to reproduce the observed data. Thus, the approach of cross-validation becomes the most prevalent method to examine the accuracy of the interpolation model. ²³

Apparently, based on unbiased estimation, the smaller the root mean square error (RMSE) between the approximate model and actual model, the higher the accuracy. Therefore, the RMSE over the domain of the model is a crucial measurement. The maximum absolute error represents the maximum deviation between the estimated value and the actual. So, it should be taken into consideration to quantify the quality of the VRPEI-Kriging model. In the construction of VRPEI-Kriging, the minimized RMSE is acquired during the solution of optimum T_w .

As an interpolation model based on OK, the assessment of VRPEI-Kriging should contain the evaluation of OK. The first is to leave out one or several sample data and to fit meta-models alternatively based on the rest of the sample data. The next step is to evaluate the predicted values of surrogate models on those left-out data by means of actual test. The accuracy of the approximate model could be ensured based on the above steps. Then, the part of the real prototype experiments results should be considered. Without a large number of actual test data, the intuitive way is to evaluate the VRPEI-Kriging with the optimization results. The positive effect on actual optimization results could be used to prove the superiority of VRPEI-Kriging.

Optimization and result analysis

The optimization was carried out following the procedure in Figure 8. First, the engineering problem was optimized based on the traditional Kriging model. By integrating the actual test results from prototype experiments, the VRPEI-Kriging meta-model was constructed for further optimization.

Figure 8.

The frame of optimization process.

In the numerical simulation of VRPEI-Kriging, a certain amount of information should be known through the deterministic virtual test. Therefore, an appropriate number of sample points should be chosen in the design of space within a certain range to perform the numerical experiments.

The choice of sample points in the design space is very vital. If the selected range of sampling space is too large, it will reduce the structural response of the overall simulation accuracy. On the contrary, the greatest point of impact on reliability may be excluded from the sample space and then affect the calculation accuracy. Considering the need of actual engineering, sampling is performed in the sampling space by means of Latin Hypercube Sampling (LHS): ²³

1. Sample 100 times for the six selected design variables by means of LHS.

2. Calculate the sampling values based on the combined simulation model to obtain the RMS value of vertical vibration acceleration of front and rear suspension, which could be taken as the corresponding response of the Kriging model.

3. Use the calculated value obtained above, in modeFRONTIER, to establish the approximation model of the vertical vibration of weighted acceleration of the front and rear axle sprung mass.

4. Select five sampling points in the design space and complete the calculation using the simulation model and the approximation model, respectively, and thus record the error between them. Set the model accuracy as 0.05. If the error meets the requirement, then stop the iteration to obtain the specific parameters of the Kriging approximation model; otherwise, add these computed values to the response domain and then go to step (3) to reconstruct the approximate model.

According to the above steps, the correlation parameters $\mathbf{\theta }$ can be estimated by minimizing $\mathrm{LN}(\mathbf{\theta })$ in equation (32). Based on the known $\hat{\mathbf{\theta }}$ , the regression coefficients $\mathbf{\mu }$ could also be estimated according to equation (27). $\hat{\mathbf{\theta }}$ and $\hat{\mathbf{\mu }}$ are obtained as $\hat{\mathbf{\theta }}=(6.431,2.776,1.877,0.682,7.116,5.242)$ and $\hat{\mathbf{\mu }}=(0.269,0.121,0.092,0.035,0.327,0.241)$ . The Kriging approximation model was obtained as follows (Figures 9 and 10).

Figure 9.

Approximation model of the RMS value of vertical vibration acceleration of the front axle sprung mass: (a) the changed a_{w_f} with V0_f and P0_f, (b) the changed a_{w_f} with V0_f and d0_f, and (c) the changed a_{w_f} with P0_f and d0_f.

Figure 10.

Approximation model of the RMS value of vertical vibration acceleration of the rear axle sprung mass: (a) the changed a_{w_r} with V0_r and P0_r, (b) the changed a_{w_r} with V0_r and d0_r, and (c) the changed a_{w_r} with P0_r and d0_r.

Different from the building of the traditional Kriging meta-model, the construction of VRPEI-Kriging cannot be separated from the sample data of actual test. Besides the above results data obtained through the virtual simulation model, in reference to the constraints and under the premise of considering actual test conditions, 10 samples were randomly selected in the sampling interval. According to the sample data, the parameters of hydro-pneumatic suspension were adjusted and the actual test was performed. (Figure 11). The corresponding experimental results were recorded and analyzed, as shown in Table 2.

Figure 11.

The layout of sensors and data acquisition (left: the layout of front acceleration sensor; middle: the layout of rear acceleration sensor; right: data acquisition system).

Table 2.

Actual experimental data.

	Front hydro-pneumatic spring			Rear hydro-pneumatic spring			The test results
	P0_f (Pa)	V0_f (L)	d_f (m)	P0_r (Pa)	V0_r (L)	d_r (m)	aw_f (m/s2)	aw_r (m/s2)
1	8.50E+06	10.0	0.018	6.50E+06	3.0	0.008	1.4971	1.1805
2	1.00E+07	6.0	0.020	1.50E+07	5.0	0.012	2.4307	1.2892
3	7.50E+06	8.0	0.014	1.20E+07	3.5	0.009	1.7531	1.3027
4	8.00E+06	3.0	0.015	1.20E+07	6.5	0.007	2.6844	1.1542
5	6.50E+06	7.0	0.018	7.50E+06	4.5	0.006	1.4686	1.3161
6	5.00E+06	3.5	0.016	8.00E+06	6.5	0.007	2.3126	1.1291
7	6.00E+06	6.0	0.020	9.50E+06	4.5	0.008	1.4182	1.4425
8	7.00E+06	4.0	0.019	1.00E+07	8.0	0.006	2.5234	1.5374
9	6.00E+06	5.0	0.020	7.50E+06	6.0	0.007	2.1326	1.2557
10	7.00E+06	6.5	0.016	1.50E+07	6.5	0.012	1.9307	1.2090

Taking the Kriging meta-model as the basis and introducing the actual experimental date based on virtual simulation, the minimized RMS of the error between the estimated values and test values was taken as the optimization target to solve the optimal test weight factor T_w . Thus, the approximation model of VRPEI-Kriging was constructed in the same way as the Kriging meta-model. Genetic algorithm, after 200 generation algorithm iterations, was used to solve the optimal solution based on Kriging and VRPEI-Kriging meta-model, respectively. The results are shown in Table 3.

Table 3.

Results of optimization.

	Front hydro-pneumatic spring			Rear hydro-pneumatic spring			The acceleration RMS
	P0_f (Pa)	V0_f (L)	d_f (m)	P0_r (Pa)	V0_r (L)	d_r (m)	aw_f (m/s2)	aw_r (m/s2)
Before optimization	5.5E+06	6.0	0.016	1.0E+07	4.0	0.009	1.3555	1.2285
Kriging	6.5E+06	9.0	0.012	9.5E+06	6.0	0.006	0.8931	0.7457
VRPEI-Kriging	6.5E+06	8.0	0.012	9.0E+06	8.0	0.006	1.1294	0.8968

RMS: root mean square; VRPEI-Kriging: virtual and real prototype experiment integrated Kriging.

According to the obtained optimal parameters, actual experiments were carried out. We compare the test results as follows (Table 4).

Table 4.

Comparison of optimization.

	Meta-model	Front suspension (m/s2)	Rear suspension (m/s2)
Before optimization	–	1.3555	1.2285
Virtual optimization	Kriging	0.8931	0.7457
	VRPEI-Kriging	1.1294	0.8968
Actual test	Kriging	1.2236	1.1217
	VRPEI-Kriging	1.1863	1.0837

VRPEI-Kriging: virtual and real prototype experiment integrated Kriging.

From the correlation tables, as shown in Figure 12, it is obvious that the optimization results based on Kriging and VRPEI-Kriging meta-model improved the ride comfort of the vehicle to the same extent.

Figure 12.

Comparison of optimization.

Just from the virtual optimization, the results, obtained based on the Kriging surrogate model, showed greater degree of improvement than on VRPEI-Kriging. For example, the acceleration RMS of front suspension was improved by 34.12% (from 1.3555 to 0.893 m/s²) according to the Kriging model, while it only improved by 16.68% (from 1.3555 to 1.1294 m/s²) based on VRPEI-Kriging. A similar conclusion can be obtained from the curve of rear suspension.

Through the verification by the actual test, we found that the acceleration RMS of front suspension was just 1.2236 m/s², actually improved by only 9.73% based on Kriging. While based on VRPEI-Kriging, the result was improved by 12.48%. It is similar for the rear suspension. The result, based on VRPEI-Kriging, was improved by 11.79%, better than 8.69% from the optimization based on the Kriging model. Compared with traditional Kriging, the data analysis showed that the optimization effect was improved by 3.05% and 3.38%, respectively, based on VRPEI-Kriging.

What is more, there is every reason to think that the optimization based on VRPEI-Kriging is more efficient by mutual comparison of the test data.

Conclusion

The ride comfort optimization of a multi-axle heavy motorized wheel dump truck was done completely based on the VRPEI-Kriging model. The complete CAE simulation model was established, including the dynamics model, the hydro-pneumatic spring model, and 3D pavement model. Under the constraints of natural frequency, the working stroke, and the dynamic load of the wheel, the vibration acceleration RMS of the sprung mass center was optimized. Considering the complexity of the engineering and the optimization cost, the optimization was executed based on the VRPEI-Kriging meta-model to find the optimal parameters for the front and rear hydro-pneumatic suspension. The results showed that the acceleration RMS was actually improved by 12.48% for front suspension and 11.79% for rear suspension. It is meaningful to relieve the fatigue of the frame and the suspension.

The VRPEI-Kriging meta-model proposed in this article was built based on traditional Kriging. The interpolation principle was deduced in detail. The VRPEI-Kriging model considers the actual test results with high reliability, while the traditional did not. However, in some particularly complex project with strong nonlinearity, the order of Kriging surrogate model obtained from the pure CAE analysis could be lowered. In that case, it would lead to the loss of the extremum or the optimal solution. Then, we need the actual test results to modify the approximate model to remedy the shortage. VRPEI-Kriging precisely contributes to cover the potential defect, especially in an intricate engineering problem whose CAE simulation model is almost impossible to be established with high accuracy. Maybe there are still some deficiencies in the VRPEI-Kriging, such as the arithmetic solution of the test weight ratio T_w . It will be further studied in later work. Nonetheless, the optimization in this article illustrates the efficiency of the VRPEI-Kriging meta-model. The ride comfort of the multi-axle heavy motorized wheel dump truck was improved effectively based on the VRPEI-Kriging model.