Improved Vehicle Vibration Control through Optimization of Suspension Parameters Using the Response Surface Method and a Non-Linear Programming with a Quadratic Lagrangian Algorithm

Abstract

Vibration-control techniques generally fall into two categories: passive methods that optimize the structure of the suspension to absorb any impact from the ground, and active methods that directly control the vertical force of the suspension by hydraulic or electric actuators when the vehicle traverses a bumpy road. In this study, a vibration-control method is described that employs both an optimal controller and suspension parameter optimization. Continuous speed bumps are implemented to simulate more complex and realistic driving conditions. First, a vehicle system is modeled using a semi-recursive multibody formulation, which allows for a more precise description of the longitudinal–vertical dynamics. Then, an optimal control method for vehicle vibration control is introduced. Second, the Latin hypercube design is utilized to analyze the response surface methodology (RSM) model. For suspension optimization, the RSM model and the non-linear programming with a quadratic Lagrangian (NLPQL) algorithm are employed. Thirdly, both passive suspension optimization and active motion control are employed for vibration control. The results indicate that the presented method can effectively control vehicle vibration, decreasing the average vibration by 30.8%. The results suggest that the novel approach can also enhance the ride comfort in autonomous vehicles traversing, e.g., a series of speed bumps.

1. Introduction

Vibration control is a common method for enhancing travel comfort in vehicles and is one of the most important evaluation criteria. The vibration affects the driver’s confidence and the passengers’ well-being and safety. The irregularity of the road is typically correlated with the vibration of a moving vehicle. In order to increase safety, speed bumps are frequently placed near essential buildings and specific areas, such as intersections with heavy traffic, tunnel entrances and exits and school and community entrances. Besides slowing down vehicles, they also tend to increase vehicle vibrations and significantly reduce ride comfort, even at safe/legal (slower) velocities. Consequently, the comfort of vehicle drivers and passengers could be enhanced by investigating how to reduce the vibration of, e.g., “smart” vehicles on roads with frequent speed bumps [1,2,3,4,5].

Commonly used vibration-control methods tend to primarily optimize the suspension parameters, such as the spring and damping characteristics. Song presented a new fuzzy sliding mode controller for vibration control of semi-active automotive suspension systems with magnetorheological dampers [6]. Xu presented a hybrid controller with hybrid acceleration driven damping algorithms to reduce the suspension spring mass and the vibration of the electric vehicle motor [7]. To capture the vibration response needed to optimize suspension parameters, accurate full-vehicle dynamics models are essential. The multibody dynamics model of a vehicle (with a detailed description of all suspensions and road/tire interactions) is well-suited to analyze the problem [8,9,10]. However, the majority of multibody-dynamics models are computationally expensive for this type of optimization [11,12,13]. Hence, the computational cost for real-life problems could be reduced by developing real-time and efficient multibody-dynamics models (e.g., recursive or semi-recursive models that use a system-tree-topology). The semi-recursive multibody model is typically more efficient than the multibody models used by the majority of commercial software applications. Pan developed a vehicle model using semi-recursive dynamics equations and validated it against commercial software solutions [14,15]. In this study, an efficient semi-recursive multibody formulation is introduced to model the vehicle dynamics [16,17,18]. This model takes each component’s dynamic properties into account and enables precise dynamic responses.

To further reduce the computational cost and increase the efficiency of the suspension optimization, several approximation- or meta-modeling methods were developed (based on limited samples of vehicle dynamics). In addition, a semi-recursive multibody model was devised to ensure that the samples represent the vehicle dynamics accurately. The widely-used approximation methods include the response surface methodology (RSM) model, the orthogonal polynomial model, and the Kriging model. They are widely used for suspension optimization to improve the comfort level in vehicles [19,20,21,22]. The RSM model approximates the function more precisely over a brief (local) range due to its use of fewer experiments, simple algebraic expressions, and design optimization simplification. In addition, the RSM model accommodates complex response relationships by choosing appropriate regression models. Due to its high efficiency and minimal cost, the RSM model was chosen as the approximation model for suspension optimization in this study.

There are two kinds of methods for engineering optimizations: (1) Stochastic methods, like genetic algorithms [23,24,25,26]. Typically, a large initial population is required to discover a (sufficiently) viable solution, unfortunately leading to make stochastic methods computationally costly [27,28]. (2) Gradient-based optimization. For an appropriate iterative refinement of the design space, these methods significantly rely on the initial starting values and demand accurate gradient calculations [29,30,31,32]. Baumal, McPhee, and Calamai compared the computational efficiency of a genetic algorithm with a gradient-based optimization algorithm for a pitch-plane vehicle model. The genetic algorithm converged on an optimal solution that offered only a 4% performance advantage over the gradient-based algorithm [33]. In addition, thousands more function evaluations are required. A gradient-based optimization method that uses optimization criteria and the RSM model is incorporated to find the optimal suspension parameters in this study.

On solving optimization problems with non-linear objective functions, Sun proposed a system-level design optimization method based on the actual operating environment and multi-objective optimization based on sequential subspace optimization to improve the performance and efficiency of electric vehicle motors [34,35]. While Shi proposed a Fuzzy-based sequential Taguchi robust optimization method to improve the motor’s comprehensive performance and optimization efficiency [36]. The non-linear programming with a quadratic Lagrangian (NLPQL) method is regarded as more suitable for the vertical vibration problems (the single-objective optimization) caused by a rough road surface addressed in this study than the previously mentioned optimization techniques. Accordingly, the global optimal solution for the optimization parameters can be efficiently and effectively found. In this study, an embedded active speed-control scheme and the (passive) optimization of the suspension parameters are combined. To the best of the author’s knowledge, the method of this study, which employs both a response surface methodology and non-linear programming with a quadratic Lagrangian (NLPQL) algorithm (to suppress the vehicle vibration), has rarely (if ever) been reported.

The highlights of this study can be summarized as follows:

–

An efficient vehicle dynamics model was developed using a semi-recursive multi-body dynamics approach, enabling an accurate description of the longitudinal and vertical dynamic response of the vehicle. Building upon this model, an optimal control algorithm for reducing vibrations when driving over frequent speed bumps was devised.

–

A Latin hypercube experimental design was utilized to increase the efficacy of simulation data collection and reduce the number of required simulations. This method assures efficient sampling across the input parameter space, enabling a thorough investigation of suspension design possibilities. In addition, the Response Surface Methodology (RSM) was applied to simulation data to construct a surrogate model for the design optimization (vibration control) problem.

–

Leveraging the quadratic RSM model developed in the previous step, a non-linear programming using a quadratic Lagrangian algorithm was applied to determine an optimal solution for the suspension parameters. The identification of the most suitable suspension parameters that meet the intended performance criteria is allowed by taking into account both the goal of minimizing vibrations and the constraints imposed by the vehicle’s maneuverability in this algorithm.

Overall, an effective and systematic approach for optimizing suspension parameters and achieving superior vibration control in vehicles is provided by the combination of the efficient vehicle dynamics model, the Latin hypercube design of experiments, the quadratic RSM model, and the non-linear programming algorithm.

The rest of the manuscript is organized as follows. In Section 2, a semi-recursive multibody method and an optimal vehicle controller is described. In Section 3, a metamodel based on the optimal Latin hypercube experimental design and surface response methodology is introduced. The accuracy of the metamodel is also investigated. In Section 4, the parameters for the suspensions via non-linear programming with a quadratic Lagrangian algorithm are optimized, in accordance with better control of the vehicle vibration. Finally, in Section 5, the outcomes of this study are summarized.

2. Vehicle-Vibration Suppression Using an Optimal Control Algorithm

2.1. Multibody Model of the Vehicle

This study modeled vehicle dynamics with an efficient semi-recursive multibody method. The virtual power principle was used at each body’s center of gravity (C.G.) to formulate the dynamic equations of the whole system. Consider an open-loop multibody system in which the first velocity transformation is used to characterize the Cartesian velocities and accelerations relative to one another in order to reduce dimensionality. It is important to note that the tree-topology of an open-loop system was used, including its connectivity and path matrix. The equations of motion of the open-loop system, in turn, can be expressed as [37,38]:

$${\mathbf{R}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{T}}^{\mathrm{T}}\overline{\mathbf{M}}\mathbf{T}{\mathbf{R}}_{\mathrm{d}}\ddot{\mathbf{z}}={\mathbf{R}}_{\mathrm{d}}^{\mathrm{T}}{\mathbf{T}}^{\mathrm{T}}(\overline{\mathbf{Q}}-\overline{\mathbf{M}}\mathbf{T}{\dot{\mathbf{R}}}_{\mathrm{d}}\dot{\mathbf{z}})$$

(1)

where, $\mathbf{T}$ denotes the path matrix, which is an upper triangular matrix that represents the system connectivity, and ${\mathbf{R}}_{\mathrm{d}}$ contains the first velocity-transformation matrix. $\overline{\mathbf{M}}$ and $\overline{\mathbf{Q}}$, respectively, denote the generalized mass matrix and external forces. $\dot{\mathbf{z}}$ and $\ddot{\mathbf{z}}$, respectively, denote the relative (joint) velocities and accelerations.

Typically, the vehicle system is a medium-scale closed-loop multibody system comprising numerous closed-loop structures, such as suspension systems. Therefore, direct modeling of this system using a recursive formulation is not possible, and preprocessing procedures are necessary to convert it into open-loop structures. Several mechanical joints were severed or a few thin rods were removed to break open the closed chains in the first stage. Subsequently, an open-loop system can be established to fully use Equation (1). However, consideration must be given to the loop-closure constrained equations that result from the cut joints and removed rods. A closed-loop multibody system’s differential-algebraic equations of motion can be obtained by introducing Lagrange multipliers. A second velocity transformation employing the Jacobian of the loop-closure constraints is introduced to eradicate the Lagrange multipliers and further reduce the dimension of the equations of motion. It is capable of describing the relative coordinates using a minimal set of independent relative coordinates. The equations of motion of closed-loop multibody systems, in turn, can be formulated as [15,39,40,41]

$${\mathbf{R}}_{\mathbf{z}}^{\mathrm{T}}{\mathbf{R}}_{\mathrm{d}}^{\mathrm{T}}{\overline{\mathbf{M}}}^{\sum }{\mathbf{R}}_{\mathrm{d}}{\mathbf{R}}_{\mathbf{z}}{\ddot{\mathbf{z}}}^{\mathrm{i}}={\mathbf{R}}_{\mathbf{z}}^{\mathrm{T}}{\mathbf{R}}_{\mathrm{d}}^{\mathrm{T}}\left[{\overline{\mathbf{Q}}}^{\sum }-{\mathbf{T}}^{\mathrm{T}}\overline{\mathbf{M}}\frac{d\left(\mathbf{T}{\mathbf{R}}_{\mathrm{d}}{\mathbf{R}}_{\mathbf{z}}\right)}{dt}{\dot{\mathbf{z}}}^{\mathrm{i}}\right]$$

(2)

where, ${\mathbf{R}}_{\mathbf{z}}$ contains the second velocity-transformation matrix. ${\overline{\mathbf{M}}}^{\sum }$ and ${\overline{\mathbf{Q}}}^{\sum }$ denote the accumulated mass matrix and external force, respectively. ${\dot{\mathbf{z}}}^{\mathrm{i}}$ and ${\ddot{\mathbf{z}}}^{\mathrm{i}}$ denote the independent velocities and accelerations, respectively.

In vehicle multibody systems, Equation (2) can be used for real-time simulation because the length of the independent relative acceleration vector ${\ddot{\mathbf{z}}}^{\mathrm{i}}$ equals the degrees of freedom of the vehicle system. The variables (independent relative accelerations) include the accelerations along the X-, Y-, and Z-axis, and the angular accelerations around the X-, Y-, and Z-axis. Thus, the equations of motion can be used to accurately describe the lateral–longitudinal–vertical coupling dynamics. This vehicle dynamics model was verified and validated by multibody system dynamics theory and commercial software applications, respectively [14,42]. In addition, a standard differential form is adopted through the equations of motion while allowing stable integration for closed-loop vehicle multibody systems by utilizing a variety of numerical integrators. The explicit 4th-order Runge–Kutta method, a good trade-off between computational efficiency and solution accuracy, is widely employed for long-term simulations.

The parameter information of the vehicle, which is studied in this paper, is shown in

Table 1. Important parameters of the vehicle model.

Parameter	Value
Degrees of freedom	16
Vehicle mass	2519 kg
Front suspension mass	38.78 kg
Rear suspension mass	188.53 kg
Tire rolling radius	0.3495 m
Distance from center of gravity (C.G.) to front axle	1767 mm
C.G. height	1200 mm
Distance from C. G. to rear axle	2333 mm
Wheelbase	4100 mm
Stiffness of front absorber	(4400, 5300) Nm/rad
Damping of front absorber	(3400, 5400) N/(m/s)
Damping of rear absorber	(5800, 9000) N/(m/s)

. The essential characteristics of the vehicle, including vertical acceleration and displacement and longitudinal velocity and displacement, can be obtained in real-time based on Equation (2). In the following section, an active speed-control scheme was incorporated into the optimized controller to reduce vehicle vibration.

2.2. Optimal Control Algorithm

In this study, speed bumps or humps were used to imitate a continuous bumpy road, according to Chinese national and industry standards. The speed bumps were usually 300 mm to 600 mm in width and 30 mm to 60 mm in height. The arc-shaped bumps with a height of 50 mm or 60 mm have better speed-control effects [43]. In this work, the road profile consists of five continuous speed bumps (with a height of 60 mm and a width of 400 mm). The bumpy road is illustrated in Figure 1. Due to the fact that the vehicle speed is below the speed threshold of the speed bumps, this investigation was conducted to ensure safe driving. In order to enhance riding comfort, vertical vibration control was subsequently implemented.

The optimal control methods include the dynamic programming method, the variational method, and the linear quadratic programming method. Due to its simplicity and reliability, the linear quadratic programming method was utilized for controller design in this study. It is critical to formulating the state-space equation between the inputs and outputs of the system to enable optimal control. Because the bumpy road results in the variation of longitudinal velocity and vertical vibration, a longitudinal–vertical coupling model is constructed based on the kinematics of actual road bumps, as described in Figure 2 [44].

In Figure 2, $r_1$ and $r_2$ represent, respectively, the radius of the tires and speed bumps. The speed bumps, $\delta r$, and $\delta r_c$ represent the vertical variation of the C.G. of the tire and the vehicle when passing the speed bumps, and $\delta \theta$ represents the rotation angle for the tire at each time-step during the numerical integration process. The vertical displacement of the tire created by the wheel passing through the arc profile of the speed bumps is converted into the vertical vibration displacement of the vehicle, thereby linking the longitudinal velocity and vertical vibration. In the modeling process, the deformation of the tire is negligible since it has little effect on the vehicle vibration. The vehicle suspension is represented by the rigidity and damping coefficients. The system state equation is described as follows:

$$\begin{array}{cc}\hfill \dot{\mathbf{x}}& =\mathbf{Ax}+\mathbf{Bu}+\Gamma \mathbf{w}\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \mathbf{y}& =\mathbf{Cx}\hfill \end{array}$$

(4)

where, $\mathbf{x}$ represents the state space vector, which is defined as $\mathbf{x}={[v_x,x_z,v_z]}^{\mathrm{T}}$. $\mathbf{y}$ represents the control output, which is defined as $\mathbf{y}=v_x$. $\mathbf{u}$ represents the control inputs, which is defined as $\mathbf{u}=a_x$. $v_x$ and $a_x$ denote the longitudinal velocity and acceleration of the vehicle, respectively. $x_z$ and $v_z$ denote the vertical displacement and velocity of the vehicle, respectively. $\mathbf{w}$ represents the system disturbance. $\mathbf{A}$, $\mathbf{B}$, $\mathbf{C}$, and $\Gamma$, respectively, represent the coefficient matrices or vectors of the controller. They are not full rank in this work.

The performance index of the optimal controller presented for vehicle control can be described as follows:

$$J=\frac{1}{2}{\int }_0^{\infty }\left[{\mathbf{x}}^{\mathrm{T}}\left(t\right)\mathbf{Q}\left(t\right)\mathbf{x}\left(t\right)+{\mathbf{u}}^{\mathrm{T}}\left(t\right)\mathbf{R}\left(t\right)\mathbf{u}\left(t\right)\right]dt$$

(5)

where, $\mathbf{Q}$ represents a positive semi-definite diagonal coefficient matrix. $\mathbf{R}$ represents a symmetric positive definite matrix. Furthermore, the structure of the optimal controller was developed—see Figure 3. The structure of the optimal controller was meticulously developed and designed. The optimal controller is resulted through a thorough and systematic approach that takes into account the controlled object (longitudinal–vertical vehicle multibody system), system dynamics and control strategies (Linear Quadratic Regulator, LQR). Under the condition that the system state equations (Equations (3) and (4)) of the controlled object are satisfied, the optimal control gain K is determined by minimizing the value of the performance index J (Equation (5)). It aims to provide a robust and efficient control output for the multibody system, ensuring that the requirement for suppressing vertical vibration is met.

In general, this optimal control will not increase the change in longitudinal acceleration rather than suppressing vertical vibration because the driving torques are controlled and updated every 100 ms smoothly in this work. The vehicle speed based on the real-time vehicle states obtained using the multibody dynamics model is modified via the optimal control. The time phase of the optimal controller can be customized or updated based on specific needs. Thus, the vehicle’s pace can be prevented from fluctuating too abruptly or frequently, allowing for more stable regulation.

In addition, the longitudinal acceleration control will no longer be implemented. The vertical vibration amplitude and peaks determine whether or not longitudinal acceleration should be regulated. The control process and the threshold are defined with precision by the optimal control strategy. After traversing the speed bumps, the vehicle’s ultimate speed remains nearly constant, thereby reducing the effect of speed on vibration suppression. As a consequence, a clear and effective vibration suppression effect is obtained.

3. Metamodels Based on the Response Surface Method

For extra suppression of the vehicle vibration, a suspension-parameter optimization was carried out based on the presented vehicle multibody model and optimal control strategy. Within the fluctuation range of the initial suspension parameters (rigidity and damping), the design was optimized. The vehicle was equipped with front double-wishbone suspensions and rear leaf-spring suspensions. The front double-wishbone suspension is suspended by a torsional spring, which is always used in vehicles due to ample assembly space. Due to their rigidity, the leaf springs’ deformation excitation is minimal when the vehicle travels over a rocky road. Consequently, the spring’s rigidity is not defined as a design variable during optimization. The optimization procedure is illustrated in Figure 4.

3.1. Optimal Latin Hypercube for the Metamodels

Design of experiment (DOE) is a systematic and efficient method for examining the relationship between the factors that affect a system and the response of that system. Fewer experiments, less time, and less expensive calculations are necessitated. In engineering and scientific investigation, more precise characteristics are typically employed. The principal experimental design methods include complete factorial design, partial factorial design, orthogonal experimental design, uniform experimental design, Latin hypercube design, and optimal Latin hypercube design.

The experimental design of the Latin hypercube has a sufficient capacity for occupying space. A lower experimental cost and a looser level classification for the impact factor are reached compared to a complete factorial design. By choosing suitable sampling points, this method can be adapted to various design-optimization problems. In actual situations, however, there is a great deal of sampling uncertainty regarding spatial uniformity. This issue is addressed by developing the optimized Latin hypercube experimental procedure. It makes it possible that all experimental points can be distributed as evenly as possible in the design space [45]. This is due to the fact that the space-filling capacity and equilibrium of the optimized Latin hypercube experimental method are superior. In the vehicle suspension, the main design parameters are the front suspension stiffness $k_f$, the front suspension damping $c_f$, and the rear suspension damping $c_r$. The optimization factors and levels are represented in

Table 2. Factors and levels of the design variables.

Factor	Number of Levels	Levels (mm)	Original Level (mm)
$k_f$ (Nm/rad)	40	4400, 4422.5, 4445, …, 5255, 5277.5, 5300	4900
$c_f$ (N/(m/s))	40	3400, 3450, 3500, …, 5300, 5350, 5400	4400
$c_r$ (N/(m/s))	40	5800, 5880, 5960, …, 8840, 8920, 9000	7400

. Suspension’s design parameters, which were all divided into 40 different levels and each level appeared once in different experiment, in order to ensure a balanced and uniform distribution of each factor level, are regarded as the design variables of the Latin hypercube experimental. Since the large variation interval of the suspension parameters can lead to numerical instability during the integration process, the values of the levels chosen here range within a stable set that guarantees a stable integration for computing and controlling vehicle dynamics based on the numerical experiments, i.e., primarily using percentage boundary constraints, centered on the initial values with upper and lower limits floating by 10% ($k_f$) and 20% ($c_f$ and $c_r$). The experimental design matrix and the dynamic response results for each sample are shown in

Table 3. Experimental design matrix.

Run	Variables			Responses (Vibration)
	${\mathit{k}}_{\mathit{f}}$ (Nm/rad)	${\mathit{c}}_{\mathit{f}}$ (N/(m/s))	${\mathit{c}}_{\mathit{r}}$ (N/(m/s))	Acceleration (m/s2)	Displacement (m)
1	5046.154	3502.564	6866.667	20.967	0.043
2	4492.308	4374.359	7523.077	20.788	0.05
3	4746.154	4938.462	8917.949	21.124	0.044
4	5207.692	4476.923	6292.308	20.921	0.042
5	4815.385	3451.282	7605.128	20.93	0.045
6	4700	4015.385	6948.718	20.771	0.048
7	4907.692	4425.641	5882.051	20.68	0.047
8	4792.308	4066.667	7933.333	20.981	0.045
9	5184.615	5143.59	6374.359	20.891	0.043
10	4446.154	4835.897	8425.641	20.944	0.049
…	…	…	…	…	…
40	4538.462	4169.231	7523.077	20.988	0.047

. A Latin hypercube matrix with dimensions is generated basing on the number of factors and the desired number of samples. Each row in the matrix represents a sample and each column represents a factor. In this experiment, there were 40 experimental samples and each with three factors. In the Latin hypercube matrix, it is ensured that each level of each factor occurs exactly once in each column and row, resulting in a balanced and stratified sampling design. The semi-recursive vehicle multi-body model constructed in the preceding section was used to simulate each sample, and the corresponding response variable values (displacement and acceleration) for each sample were recorded.

When a vehicle travels over a continuous series of speed bumps, road vertical vibrations can be produced owing to the excitation of the speed bump, which can be characterized by vibration displacements and accelerations. The peaks and averages of the vibration displacements and accelerations are therefore defined as “responses”. They are also treated as the design optimization problem’s optimization objectives. Using the optimized Latin hypercube method, experimental samples and corresponding responses are collected in order to develop a metamodel, which will be discussed next.

3.2. RSM Model for the Vibration Control

The metamodel model drastically improves the optimization process’s efficacy while maintaining sufficient precision. In the past decade, several metamodels, such as the RSM, Kring, and radial basis function models, were utilized extensively in engineering optimizations. With fewer experiments, the RSM model accurately approximates the function within a local range, thereby facilitating the optimization process [46,47]. Furthermore, the RSM model can be applied in complex response relationships by using suitable regression models. In this work, the linear polynomial and quadratic polynomial expressions of the RSM are used and compared in terms of accuracy, which can be described as:

$$\begin{array}{cc}\hfill & y(x)={\beta }_0+\sum _{i=1}^k{\beta }_i{x}_i\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & y(x)={\beta }_0+\sum _{i=1}^k{\beta }_i{x}_i+\sum _{i=1}^k{\beta }_{ii}{x}_i^2+\sum _{1\le i<j\le k}^k{\beta }_{ij}{x}_i{x}_j\hfill \end{array}$$

(7)

where, the factors $x_i$ (i = $1,2,3$) are suspension’s design parameters $k_f$, $c_f$ and $c_r$, therefore the number of variables k is 3. The response variable y is the predicted vibration acceleration, ${\beta }_0$ is the intercept term, ${\beta }_i$ (i = $1,2,3$) are the linear polynomial regression coefficient associated with the factors $x_i$. Equation (7) includes not only the linear terms, but also the interaction and quadratic terms. The second term on the right-hand side of the equation represents the quadratic terms, while the third term represents the interaction terms. ${\beta }_{ij}$ ($1\le i<j\le 3$) represents the interaction terms coefficient. These regression coefficients provide information regarding the effect of each factor and their interactions on the response variable, enabling prediction and optimization of the response within the experimental range.

The RSM model is founded on the “two replacements at a time” method to determine whether all items can be replaced, beginning with a constant term. This is carried out to identify coefficients that result in the smallest error. Using this method, the linear and quadratic RSM models for vibration displacements and accelerations can be developed. According to the Pareto analysis, both the front suspension stiffness $k_f$ and the rear suspension damping $c_r$ are the most significant parameters for vehicle vibration. Thus, the RSM models are described making use of these two parameters—see Figure 5. The effect of the main factors (suspension parameters) on the vibration responses is illustrated.

The overall accuracy of the RSM models can be evaluated and tested by error functions. $W_T$, $W_R$ and $W_E$ are the sum of squared deviations of the samples, the sum of squared regression deviations and the sum of squared residual deviations, respectively. The corresponding expressions are written as follows:

$$\begin{array}{cc}\hfill & W_T=\sum _{i=1}^n{(y_i-\overline{y})}^2\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & W_R=\sum _{i=1}^n{({\tilde{y}}_i-\overline{y})}^2\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & W_E=\sum _{i=1}^n{(y_i-{\tilde{y}}_i)}^2\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \overline{y}=\frac{1}{n}\sum _{i=1}^n{y}_i,\hfill \end{array}$$

(11)

where, $y_i$ represents the value of the responses of the i-th sample, ${\tilde{y}}_i$ represents the approximated value of the responses of the i-th sample, n represents the total number of samples, and $\overline{y}$ represents the mean value of the responses of the n samples. The accuracy of the RSM models is evaluated using the ${\mathrm{R}}^2$ determination, which is expressed by the following equation:

$${\mathrm{R}}^2=1-\frac{W_E}{W_T}$$

(12)

where, the numerator component represents the sum of the squared differences between the actual and predicted values, and the denominator component represents the sum of the squared differences between the actual and mean values. The closer the value of ${\mathrm{R}}^2$ is to 1, the more accurate the RSM model is. The closer the value of ${\mathrm{R}}^2$ is to 0, the worse the model’s fit is. Because the ${\mathrm{R}}^2$ determination is affected by the number of samples and cannot truly and quantitatively explain the accuracy of the RSM model. The ${\mathrm{R}}_{adj}^2$ is introduced to balance the effect of the number of samples. It can be expressed as

$${\mathrm{R}}_{adj}^2=1-\frac{(1-R^2)(n-1)}{(n-p-1)}$$

(13)

where, p is the number of features, and the value range of the ${\mathrm{R}}_{adj}^2$ is from 0 to 1.

Table 4. Summary of RSM model statistics.

	Linear		Quadratic
	${\mathbf{R}}^2$	${\mathbf{R}}_{\mathit{adj}}^2$	${\mathbf{R}}^2$	${\mathbf{R}}_{\mathit{adj}}^2$
First peak	0.999	0.999	0.999	0.998
Second peak	0.994	0.994	0.999	0.998
Third peak	0.973	0.970	0.999	0.998
Fourth peak	0.537	0.484	0.960	0.942
Fifth peak	0.890	0.878	0.995	0.992
Average peak	0.981	0.979	0.999	0.998
Average displacement	0.999	0.999	0.999	0.999

describes the summary statistics of the ${\mathrm{R}}^2$ and ${\mathrm{R}}_{adj}^2$ determination for the linear and the quadratic RSM model. For the fourth and fifth peaks, the linear RSM model has a relatively poor fit accuracy. This is because the characteristics of the vertical acceleration responses for the fourth and fifth peaks are highly non-linear, and the vehicle dynamics cannot be adequately described by the linear RMS model. In addition, the quadratic RSM model has a higher overall fit accuracy, resulting in a superior fit for each acceleration peak. This also indicates that the quadratic RSM model’s predicted responses are closer to the actual responses. It can be seen from Table 4 that the minimum value for ${\mathrm{R}}_{adj}^2$ is 0.942, which approaches 1. This suggests that the quadratic RSM model is capable of producing accurate outcomes. Consequently, it could be substituted for the actual physical model for design optimization. However, the linear RSM model could be useful if focusing on the minimum computational cost. There is a trade-off between computational cost and approximation accuracy.

4. Optimization of the Suspension

4.1. Non-Linear Programming Using the Quadratic Lagrangian (NLPQL) Algorithm

The vehicle system’s capacity to transport body mass and absorb the impact of rough roads is determined through the suspension’s stiffness. Moreover, caused by suspension’s damping, the vehicle vibration is reduced. The greater the damping coefficient, the greater the effect on vibration suppression. Nonetheless, a negative effect on system rigidity is realized when it is very large. In other terms, damping can diminish the effectiveness of elastic components and damage vehicle components. Therefore, an optimization of the suspension’s parameters is required to determine the optimal suspension rigidity and damping. Given the constraints imposed by the vertical displacement of the center of mass, the vibration acceleration is minimized in this manner when the vehicle traverses the rough road.

Several optimization methods (i.e., algorithms) can be used: numerical optimization, direct search, global search, gradient optimization, etc. In this investigation, the NLPQL gradient optimization algorithm is utilized. Based on the current solution, an updated solution in the design space is determined by the gradient algorithm using a forward approach and an appropriate search step. The concept underlying the NLPQL algorithm is executed to express the objective function as the second-order Taylor equation and to linearly process the constraints [48,49]. Thus, the problem of design optimization can be transformed into a quadratic programming problem. Consistent with the prior quadratic RSM model, a better convergence is achieved. The design optimization problem is now defined, and the related pseudocode is described as Algorithm 1:

Algorithm 1 Design optimization.

Require: Minimize vibration acceleration ($\min f\left(\mathbf{x}\right)$, $\mathbf{x}=(k_f,c_f,c_r)$) Ensure: 4400 Nm/rad ≤ $k_f$ ≤ 5300 Nm/rad  3400 N/(m/s) ≤ $c_f$ ≤ 5400 N/(m/s)  5800 N/(m/s) ≤ $c_r$ ≤ 9000 N/(m/s)  $z_c$ ≤ $z_0$ + $\epsilon$ = 0.05 + 0.003 = 0.053 m  (Note: $z_c$ represents the average peak of the vibration displacements,  $z_0$ represents the threshold of the vibration displacements.)

In the above algorithm, the design variables are the front suspension stiffness $k_f$, front suspension damping $c_f$, and rear suspension damping $c_r$, which are defined for a specific range using our numerical experiences to assure numerical stability. The objective is optimization to minimize the vibration acceleration when the vehicle traverses a road with speed bumps. Minimum ground clearance was also used as a design constraint in consideration of the driving manoeuvrability of the vehicle. This represents that the constraint stipulates that the vehicle’s vertical displacement must not exceed a predetermined limit, which can be defined by a threshold value and a variable interval $\epsilon$.

For the purpose of analyzing the efficiency of the NLPQL algorithm and comparing its optimal results, the large-cardinality full-factor method was utilized for verification and comparison. Transversal analysis for all possible factor combinations at all levels is necessitated in full-factor design. It is presumed that the effects of various factors on the results of quantitative observation are equivalent. Each experiment has a distinct level for each variable. Thus, a vast quantity of information covering the entire design space was obtained. Nevertheless, while it can accurately assess the effects of each factor at all levels, it is cumbersome and time-consuming for large-scale issues. In the NLPQL algorithm, each factor is divided into 10 levels to ensure that the sampling density has sufficient spatial coverage accuracy. The total number of samples is 1000 (groups), and the interval sampling values for different levels are 100, 200, and 320, respectively. The optimization results (using the RSM models and NLPQL algorithm) are illustrated in Figure 6 and Figure 7. Due to space constraints, it is essential to note that only the optimal stiffness and damping of the front suspension are shown in Figure 6 and Figure 7.

The optimization results for both the NLPQL and the full factor methods are shown in

Table 5. Suspension parameters before and after the optimization.

Method	Stiffness of Front Absorber (Nm/rad)	Damping of Front Absorber (N/(m/s))	Damping of Rear Absorber (N/(m/s))
Initial parameters	4901	4400	7400
NLPQL method	4607	4270	5800
Full-factor method	4550	4200	5800

. Optimized parameters that are consistent with the results of the full-factor method are yielded through the NLPQL method. Any deviations occur within an allowable interval that satisfies the design specifications. Nevertheless, compared to the time-consuming full-factor method, the NLPQL method is significantly more efficient while maintaining adequate precision.

4.2. Vibration-Control Results

To visualize the effects of vibration control utilizing the optimal suspension parameters and the optimal controller, a vehicle traversing a series of speed obstacles was subjected to a dynamic maneuver. The vehicle’s initial velocities were set to 25 km/h and 45 km/h. The vibration control results are shown in Figure 8. The first five values denote the peak value for vibration acceleration for the speed bump series (as shown in Figure 1). The last peak value denotes the average vibration acceleration on the bumpy road. Two different initial vehicle speeds were used. The control effects were visualized using four distinct control methods: no optimal controller and optimal parameters (method A), optimal controller but without optimal parameters (method B), PID controller with optimal parameters (method C), and optimal controller with optimal parameters (method D). Figure 8 illustrates that when the vehicle traverses a bumpy road, its vertical acceleration magnitude increases with the initial speed (from 25 km/h to 45 km/h). In addition, the maximal acceleration of vibration occurs at the third peak when the initial velocity is 45 km/h (relatively high). This is due to the fact that the vehicle’s speed does not diminish significantly over the bumps. Nonetheless, the vibration builds up over time. This indicates that the vehicle is unable to assimilate sufficient energy and will likely experience greater vibrations at a higher initial speed on a rough road. In addition, the maximal acceleration of vibration decelerates as the initial speed increases.

Table 6. The peaks of vertical acceleration (25 km/h).

25 km/h	Control Method
	A	B	C	D
First peak	14.60	13.92 (↓ 4.61%)	13.56 (↓ 7.09%)	13.39 (↓ 8.29%)
Second peak	15.72	13.44 (↓ 14.49%)	14.10 (↓ 10.31%)	13.36 (↓ 14.99%)
Third peak	15.68	9.73 (↓ 37.93%)	11.75 (↓ 25.04%)	10.44 (↓ 33.43%)
Fourth peak	15.65	9.73 (↓ 37.78%)	8.52 (↓ 45.42%)	8.54 (↓ 45.42%)
Fifth peak	15.24	7.53 (↓ 50.55%)	8.92 (↓ 41.45%)	7.46 (↓ 51.03%)
Mean value	15.38	10.87 (↓ 29.28%)	11.33 (↓ 26.33%)	10.64 (↓ 30.82%)

and

Table 7. The peaks of vertical acceleration (45 km/h).

45 km/h	Control Method
	A	B	C	D
First peak	17.93	18.46 (↑ 2.92%)	17.20 (↓ 4.11%)	17.51 (↓ 2.35%)
Second peak	21.48	19.63 (↓ 8.62%)	19.54 (↓ 9.00%)	18.59 (↓13.44%)
Third peak	24.49	23.85 (↓ 2.63%)	26.44 (↓ 7.95%)	24.58 (↑ 0.37%)
Fourth peak	22.00	22.34 (↑ 1.53%)	20.72 (↓ 5.80%)	21.55 (↓ 2.03%)
Fifth peak	22.58	20.38 (↓ 9.75%)	21.94 (↓ 2.84%)	20.10 (↓ 11.00%)
Mean value	21.70	20.93 (↓ 3.54%)	21.17 (↓ 2.43%)	20.47 (↓ 5.67%)

present the vertical vibration suppression numerical results for the four different control strategies at speeds of 25 km/h and 45 km/h. The efficacy of the other three control methods is determined through a benchmark served by the numerical results derived from the simulation of the dynamics without any control (method A). It can be seen from the tables that, at the lower speed (25 km/h), the average suppression for the five peaks is 29.28% for the optimal controller, and the maximum suppression can reach 50.55%. By optimizing the suspension parameters, the average suppression can be improved to reach 30.8%. A more significant suppression effect is realized than the PID controller with optimal parameters (the average suppression 26.33%). For an initial speed of 45 km/h, the average suppression for the five peaks is 3.54% for the optimal controller, and the maximum suppression reached 9.75%. Taking into consideration the optimal suspension parameters, the optimal controller suppresses peaks by an average of 5.7%, whereas the PID controller only suppresses peaks by 2.43%, which is consistent with the findings at modest vehicle speeds. Finally, the vehicle vibration control method, which uses the RSM model, the NLPQL algorithm, and the optimal controller, is the most effective. Passive suspension-parameter optimization is integrated with an active speed-control scheme from the new vibration control method.

Based on the results shown in Figure 6 and Figure 7 and Table 6 and Table 7, the following conclusions can be drawn:

–

Intricate coupling effects are exhibited by the three suspension parameters when the vehicle traverses a series of speed bumps.

–

One peak value is insufficient to characterize the vibration performance. Instead, the average acceleration of the five shocks is a more appropriate metric for describing vibration control.

–

Consistent with real-world observations, opposing trends are displayed via the results for vehicle accelerations and displacements. To optimize the design, vehicle displacements are therefore converted into constraints.

The NLPQL method is used as a computationally efficient and scalable mathematical optimization technique. However, it is worth noting that the initial point selection has a significant impact on the algorithm’s convergence and ultimate results. In the research, this technique is used to optimize the suspension parameters in order to minimize the vibration acceleration of a vehicle traveling over a rough road. As the objective function, the previous quadratic RSM model is utilized, offering superior convergence compared to the conventional suspension dynamics model. To avoid falling into a local optimum solution that does not satisfy the actual engineering requirements, a specific range of vehicle suspension parameter variables is defined based on numerical experience. In future applications, suspension optimization parameters may be used directly in the control of the vehicle’s active suspension, i.e., the suspension stiffness and damping are adjusted to the current optimum parameters based on the road conditions, and then the vehicle speed is adjusted to further suppress vehicle vibrations using a controller such as the optimal controller.

5. Conclusions

The vibration problem that occurs when a vehicle traverses a series of speed obstacles was partially addressed by developing a suspension-optimization method employing an optimal control algorithm. In addition, a multibody vehicle model was created using a semi-recursive multibody formulation that accurately depicts the critical vehicle characteristics. For the vehicle dynamics model, an optimal controller was used to suppress all vertical vibrations. And a quadratic RSM model was proposed with an optimal Latin hypercube approach to improve the vehicle vibration controller. The RSM model was analyzed and its accuracy was confirmed. The suspension parameters were optimized using an NLPQL algorithm and the RSM model, with the full-factor method used to validate the optimization results.

The results indicate that the new method of vibration control, which combines a meta-model-based suspension optimization with an optimal control algorithm, enables more effective suppression of vibrations on rough roads. In addition, the average vibration acceleration decreased by more than 30% at an initial speed of 25 km/h. The average and maximum vibration acceleration decreased by more than 5% and 10% at 45 km/h, respectively. Overall, these results suggest that the new method could be used to control the vibration of autonomous vehicles to improve passenger comfort and safety. Future research could focus on the parameter-uncertainty analysis for vehicle vibration control.