# Active Suspension Control Based on Estimated Road Class for Off-Road Vehicle

Authors: Mingde Gong; Haohao Wang; Xin Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Road input can be provided for a vehicle in advance by using an optical sensor to preview the front terrain and suspension parameters can be adjusted before a corresponding moment to keep the body as smooth as possible and thus improve ride comfort and handling stability. However, few studies have described this phenomenon in detail. In this study, a LiDAR coupled with global positioning and inertial navigation systems was used to obtain the digital terrain in front of a vehicle in the form of a 3D point cloud, which was processed by a statistical filter in the Point Cloud Library for the acquisition of accurate data. Next, the inverse distance weighting interpolation method and fractal interpolation were adopted to extract the road height profile from the 3D point cloud and improve its accuracy. The roughness grade of the road height profile was utilised as the input of active suspension. Then, the active suspension, which was based on an LQG controller, used the analytic hierarchy process method to select proper weight coefficients of performance indicators according to the previously calculated road grade. Finally, the road experiment verified that reasonable selection of active suspension’s LQG controller weightings based on estimated road profile and road class through fractal interpolation can improve the ride comfort and handling stability of the vehicle more than passive suspension did.


---

## Body

## 1. Introduction

In the running course of a vehicle, the road excitation characteristics considerably influence its ride comfort and handling stability [1, 2]. Suspension is an important component that plays a critical function in the transmission of force between vehicle’s tire and body and alleviates the impact of the road surface. Consequently, obtaining the pavement roughness in front of a vehicle to control its suspension and improve ride comfort and handling stability has been attracting the attention of researchers.Existing active suspension regulates its parameters only when the road interferes with the vehicle [1, 2]. The use of the sensor system to detect the road on which a vehicle will be running and of the built-in controllers to process the related data can maximise the potential of active suspension. Therefore, a suspension system with a preview function is vital to future intelligent chassis.Stereo cameras [3, 4] and LiDAR sensors [5–7] can be adopted for this function by measuring the road height information in front of a running vehicle with certain accuracy. In [2], the use of a low-cost ultrasonic sensor was examined to compute road surface height estimation in real time by using simple algorithms due to the drawbacks of sensors (i.e., stereo vision or LiDAR), such as costliness and computation complexity. However, the ultrasonic sensor does not have the high precision of a laser radar. In [7], a LiDAR coupled with a global positioning system (GPS) and an inertial navigation system (INS) was used to measure two off-road surfaces, but details of the measurement and data processing were not provided. In [8], the author proposed the use of the overall response of the preceding vehicles to generate preview controller information for follower vehicles; no sensor is used to measure the terrain online.Although many scholars have proposed the adoption of preview information for active suspension, detailed descriptions have not been provided. Instead, researchers have merely introduced concepts or assumed the road height profile to be known and obtained by statistical methods. This study proposes an approach that uses LiDAR coupled with INS and GPS to obtain a preview of the 3D point cloud of the terrain in a geodetic coordinate system and extract the road height profile in front of a vehicle. This work also elaborates the methods and steps of data processing. The road roughness grade is identified on the basis of the aforementioned information. Then, the Linear-Quadratic-Gaussian (LQG) controller’s weight coefficients of performance indicators are optimised in different road roughness grades. Finally, we select proper weight coefficients that correspond to previously calculated road grades to improve the vehicle’s ride comfort and handling stability.This study is organised as follows. Section2 provides the mobile mapping system. Section 3 presents the data processing of sensors. Section 4 uses the methods of the two preceding sections to experiment and obtain the 3D point cloud of the test site. Section 5 adopts the road height profile information utilised in fractal interpolation theory in Section 6 to improve accuracy. This profile information is extracted from the 3D point cloud by using the inverse distance weighting (IDW) interpolation method. Section 7 solves the grade of the road height profile’s roughness. Section 8 establishes the suspension vertical dynamic model and designs the LQG controller for active suspension. Section 9 calculates the controller’s weight coefficients of performance indicators for different road grades using the analytic hierarchy process (AHP). In Section 10, a road experiment verifies that reasonable selection of LQG controller weightings based on estimated road profile and road class through fractal interpolation can improve the ride comfort and handling stability of the vehicle more than passive suspension does.


## 2. Mobile Mapping System

To control active suspension with preview information, a device that maps the front terrain of the running vehicle must be developed. Laser radar scanning is a rapid 3D measurement method and has existed for a while. But it is seldom used for the realization of active suspension. The mobile LiDAR system is mainly composed of a LiDAR sensor, GPS, and INS (Figure1). In the running course, the LiDAR is mounted in front of the vehicle and returns the 3D point cloud of the road on the basis of its own coordinate system. One vital challenge in accurate mapping is estimating the attitude of LiDAR, which affects the precision of the point cloud. To solve this problem, the INS, which is installed near the LiDAR, collects the attitude angle that is used to compensate for motion during scan acquisition to improve the precision of data in real time. We need to realise the accurate location of the LiDAR to spread the 3D point cloud on the earth surface. Although GPS has higher positioning accuracy than other approaches, it is vulnerable to a wide variety of interferences, such as the multipath effect from radars, electromagnetic interference, and block of signals [9, 10]. INS is a self-contained system that consists of inertial measurement units, which can provide dynamic measurement for a short period with high frequent update. It is a complete autonomous navigation system with remarkable concealment, strong anti-interference capability, and immunity to meteorological conditions [11]. However, its measurement error may accumulate over time due to drifting effects [12]. The integration of GPS and INS has been proposed and widely implemented for vehicle applications to leverage the strengths of the two systems and to offset their individual drawbacks because they are more advantageous than single navigation systems [13, 14]. In these integrated systems, the Kalman filter is a popular fusion method that has been used in recent years due to its practicability and suitability [15]. The centimeter-level positioning accuracy can be obtained in real time by using real-time kinematic (RTK) technology.Figure 1
Mobile mapping system.


## 3. Point Cloud Data Processing

The multilinear LiDAR can launch more one laser lines whilst rotating around the laser emission centre at a certain frequency and returns the point cloud of the obstacle in the form of relative distanceLi and angle θ, which depend on the LiDAR’s coordinate system. The point cloud is stacked in front of the vehicle instead of being spread on the earth surface because the LiDAR coordinate system moves in real time with the vehicle. To resolve this problem, in the first step, we should coordinate the transformation of the point cloud. Relative distance Li and angle θ of the point cloud are transformed into the Cartesian coordinate. In the second step, data fusion algorithm is then implemented. GPS/INS integrated navigation system is used to transform the data obtained from first step into the WGS-84 geodetic coordinate system, which is stationary.


### 3.1. Coordinate Transformation of Point Cloud

The LiDAR’s vertical angle resolution is marked asα. Figure 2 shows that OL-XLYLZL is a LiDAR coordinate system whose origin OL coincides with the laser emission centre. The angle between the instantaneous laser line and the XL axis is defined as θ, and Li is the distance between OL and the obstacle. According to geometric relations, the 3D point cloud is represented in the Cartesian coordinate system as follows:(1)xLyLzL=Li·cos⁡αi∗·cos⁡θcos⁡αi∗·sin⁡θsin⁡αi∗αi∗=2α,i=1α,i=2-α,i=3-2α,i=4θ∈-60∘,50∘Figure 2
Operation mode of ibeo LUX 2010.


### 3.2. Data Fusion

The coordinate systems of the mobile mapping system include the LiDAR coordinate systemOL-XLYLZL, INS coordinate system OI-XIYIZI, local horizontal coordinate system OG-XGYGZG, and WGS-84 geodetic coordinate system O84-X84Y84Z84 (Figure 3).Figure 3
Coordinate systems of mobile mapping system.The coordinate values of the point cloud[xL,yL,zL] based on the LiDAR coordinate system are transformed to the WGS-84 geodetic coordinate system as follows:(2)xW84yW84zW84=xGPSyGPSzGPS+RWRNRMxLyLzL+ΔxLIΔyLIΔzLI-ΔxIGΔyIGΔzIGwhere RM is the rotation matrix between LiDAR and INS, whose element is calculated from the installation angle error of two sensors. The error angles around x, y, and z are β, α, and γ, which represent roll, pitch, and course angle errors, respectively, during the installation process.(3)RM=Rγz·Rαy·RβxRγz=cos⁡γsin⁡γ0-sin⁡γcos⁡γ0001Rαy=cos⁡α0-sin⁡α010sin⁡α0cos⁡αRβx=1000cosβsinβ0-sinβcosβwhere RN is the rotation matrix between INS and the local horizontal coordinate system and the matrix element is calculated from the measured values of INS, where φ, ω, and k represent pitch, roll, and course angles, respectively.(4)RN=Rκ·Rφ·RωRκ=cos⁡κsin⁡κ0-sin⁡κcos⁡κ0001Rφ=cos⁡φ0-sin⁡φ010sin⁡φ0cos⁡φRω=1000cos⁡ωsin⁡ω0-sin⁡ωcos⁡ωwhere RW is the rotation matrix between the local horizontal coordinate system and the WGS-84 geodetic coordinate system and the matrix element is calculated from the measured values of the loosely coupled GPS/INS system, where B, L, and H represent latitude, longitude, and altitude, respectively.(5)RW=-cos⁡Lsin⁡B-sin⁡L-cos⁡Lcos⁡B-sin⁡Lsin⁡Bcos⁡L-sin⁡Lcos⁡Bcos⁡B0-sin⁡Bwhere ∆xLI,∆yLI,∆zLIT are the offset values between the LiDAR and the INS coordinate system. ∆xIG,∆yIG,∆zIGT represent the offset values between the INS and the GPS coordinate system.The location of GPS in the WGS-84 geodetic coordinate system can be expressed as follows:(6)xGPSyGPSzGPS=a1-e2sin2B+H·cos⁡B·cos⁡La1-e2sin2B+H·cos⁡B·sin⁡La1-e2sin2B·1-e2+H·sin⁡Bwhere α is a semimajor axis of WGS-84 and e is the first eccentricity.(7)a=6378137me2=0.006694379


## 3.1. Coordinate Transformation of Point Cloud

The LiDAR’s vertical angle resolution is marked asα. Figure 2 shows that OL-XLYLZL is a LiDAR coordinate system whose origin OL coincides with the laser emission centre. The angle between the instantaneous laser line and the XL axis is defined as θ, and Li is the distance between OL and the obstacle. According to geometric relations, the 3D point cloud is represented in the Cartesian coordinate system as follows:(1)xLyLzL=Li·cos⁡αi∗·cos⁡θcos⁡αi∗·sin⁡θsin⁡αi∗αi∗=2α,i=1α,i=2-α,i=3-2α,i=4θ∈-60∘,50∘Figure 2
Operation mode of ibeo LUX 2010.


## 3.2. Data Fusion

The coordinate systems of the mobile mapping system include the LiDAR coordinate systemOL-XLYLZL, INS coordinate system OI-XIYIZI, local horizontal coordinate system OG-XGYGZG, and WGS-84 geodetic coordinate system O84-X84Y84Z84 (Figure 3).Figure 3
Coordinate systems of mobile mapping system.The coordinate values of the point cloud[xL,yL,zL] based on the LiDAR coordinate system are transformed to the WGS-84 geodetic coordinate system as follows:(2)xW84yW84zW84=xGPSyGPSzGPS+RWRNRMxLyLzL+ΔxLIΔyLIΔzLI-ΔxIGΔyIGΔzIGwhere RM is the rotation matrix between LiDAR and INS, whose element is calculated from the installation angle error of two sensors. The error angles around x, y, and z are β, α, and γ, which represent roll, pitch, and course angle errors, respectively, during the installation process.(3)RM=Rγz·Rαy·RβxRγz=cos⁡γsin⁡γ0-sin⁡γcos⁡γ0001Rαy=cos⁡α0-sin⁡α010sin⁡α0cos⁡αRβx=1000cosβsinβ0-sinβcosβwhere RN is the rotation matrix between INS and the local horizontal coordinate system and the matrix element is calculated from the measured values of INS, where φ, ω, and k represent pitch, roll, and course angles, respectively.(4)RN=Rκ·Rφ·RωRκ=cos⁡κsin⁡κ0-sin⁡κcos⁡κ0001Rφ=cos⁡φ0-sin⁡φ010sin⁡φ0cos⁡φRω=1000cos⁡ωsin⁡ω0-sin⁡ωcos⁡ωwhere RW is the rotation matrix between the local horizontal coordinate system and the WGS-84 geodetic coordinate system and the matrix element is calculated from the measured values of the loosely coupled GPS/INS system, where B, L, and H represent latitude, longitude, and altitude, respectively.(5)RW=-cos⁡Lsin⁡B-sin⁡L-cos⁡Lcos⁡B-sin⁡Lsin⁡Bcos⁡L-sin⁡Lcos⁡Bcos⁡B0-sin⁡Bwhere ∆xLI,∆yLI,∆zLIT are the offset values between the LiDAR and the INS coordinate system. ∆xIG,∆yIG,∆zIGT represent the offset values between the INS and the GPS coordinate system.The location of GPS in the WGS-84 geodetic coordinate system can be expressed as follows:(6)xGPSyGPSzGPS=a1-e2sin2B+H·cos⁡B·cos⁡La1-e2sin2B+H·cos⁡B·sin⁡La1-e2sin2B·1-e2+H·sin⁡Bwhere α is a semimajor axis of WGS-84 and e is the first eccentricity.(7)a=6378137me2=0.006694379


## 4. Terrain Measurement

In this study, integrated navigation system of MEMS is adopted, of which core components are the gyroscope, accelerometer, and the high-performance Beidou/GPS receiver. The position, velocity, azimuth angle, attitude angle, three-axis acceleration, and angular velocity can be output by the Kalman filter integrated navigation algorithm that can solve the problem of INS error drift over time and RTK technology. It enters the pure inertial navigation mode without satellite navigation signal and can maintain high navigation accuracy in a short time.The ibeo LUX 2010® selected for project implementation detects objects and their distance by means of laser beams. It scans the surroundings with several rotating laser beams, receives the echoes with a photo diode receiver, processes the data by means of a time of flight calculation, and issues the processed data via the interfaces Ethernet and/or CAN.The specifications of integrated navigation system and LiDAR are shown in Tables1 and 2.Table 1
Specification of integrated navigation system.
Azimuth accuracy Attitude accuracy Positioning accuracy Speed accuracy ≤ 0 . 1 ∘ ≤ 0 . 1 ∘ ≤ 2 c m ≤ 0.1 m / s  Interface Data update rate Voltage Power waste  R S 232 ≤ 500 H z D C : 9 V ~ 36 V ≤ 7 WTable 2
Specification of LiDAR.
Scanning frequency Angular resolution Number of laser lines Laser wavelength 50 H z 0 . 5 ∘ 4 905 n m  Laser grade Vertical angular resolution Voltage Interface  C l a s s    I α = 0 . 8 ∘ D C : 9 V ~ 27 V C A N / E t h e r n e tThe calibration of the integrated navigation system and the LiDAR needs to be carried out. The main work is to measure the relative position and angular error between the sensors.The sensor position center is given by the manufacturer, so the relative position error between sensors is measured by a ruler, a triangle ruler, and the like.Sensors are mounted on an adjustable platform. Integrated navigation system’s pitch angle and roll angle are adjusted to 0 according to its output value, and the heading angle is made to coincide with the longitudinal symmetry axis of the vehicle by using a laser scale. Based on the adjusted integrated navigation system, the relative angle error between sensors is measured by a level meter.The data rate of LiDAR is 50Hz and that of integrated navigation system is 500Hz. In the process of data solution, a unified time stamp and Kalman filter prediction model are used to solve the problem of inconsistent working frequency of multisensor. The vehicle computer time is regarded as the unified time stamp of the system. According to Kalman filter prediction model, the sensor data are converted to the working frequency of the LiDAR, thus forming the multisensor data at the same frequency.The dirt terrain around a university was selected as the test site (Figure4). Multithreading technology was used to ensure simultaneous data acquisition from integrated navigation system and LiDAR. The contact surface between the tire and the ground is the reference zero surface. Figure 5 shows the road map of the selected dirt terrain after data processing in MATLAB®.Figure 4
Test site.Figure 5
A road map of dirt terrain.Table3 shows the point cloud properties of the dirt terrain.Table 3
3D point cloud properties of dirt terrain.
Point numbers 346036 Longitudinal position range Minimum 4.0002 (m) Maximum 20.9994 (m) Lateral position range −4.9935 (m) 2.9998 (m) Vertical height range −0.1477 (m) 0.1997 (m)


## 5. Extraction of Road Height Profile

The 3D point cloud of the dirt terrain acquired with LiDAR inevitably suffers from noise contamination and contains outliers [16–18] due to the limitations of sensors [19], the inherent noise of the acquisition device [20], and the lighting or reflective nature of the surface or artefact in the scene [21].Therefore, the raw point cloud should be filtered to obtain accurate data that are suitable for further processing. For this project, we only consider the road that the vehicle is going to pass through and other road elevation information is ignored to simplify the calculation. Thus, the use of the road height profile extracted from the 3D point cloud of the terrain is beneficial. This paper assumes that the road ahead is empty. In the follow-up study, vehicle trajectory planning should be carried out in a short time according to the steering wheel angle.


### 5.1. Filtering of Point Cloud

The Point Cloud Library (PCL) is a large-scale open project [22] for 2D/3D image and point cloud processing, the framework of which contains numerous state-of-the-art algorithms, including filtering, feature estimation, surface reconstruction, registration, model fitting, and segmentation. For example, these algorithms can be used to filter outliers from noisy data, stitch 3D point clouds together, segment relevant parts of a scene, extract key points, and compute descriptors to recognise objects in the world on the basis of their geometric appearance and create and visualise surfaces from point clouds.The StatisticalOutlierRemoval filter in PCL is used to filter outliers. Compared with the unused filter, the denoised 3D point cloud of the dirt terrain decreased from 346036 to 322628 (Figure6).Figure 6
Denoised road map of the dirt terrain.


### 5.2. Extraction of Road Height Profile

We can predict the trajectory of a vehicle in advance according to the vehicle’s driving direction. We assume that the vehicle drives in a straight line, and the equation for the hypothetical driving track is shown in the following formula:(8)longitudinalrange∈4m,21mlateralrange=-2mThe 2D road height profile, whose accuracy should correspond with that of the 3D point cloud, is discrete because the 3D point cloud data is discrete. Therefore, we should calculate the accuracy and then extract the road height profile according to the calculated precision.Since the resolution of point cloud is affected by many factors, such as the number of lasers in LiDAR and vehicle speed, the average resolution of point cloud is used in the subsequent calculation.The volume of the denoised 3D point cloudV and the single point occupancy volume v can be written as follows:(9)V=rangex×rangey×rangez(10)rangex=20.9994-4.0003=16.9991rangey=2.9998+4.9935=7.99330rangez=0.1977+0.1477=0.34540(11)v=Vthenumberof3DpointcloudWe regardv as the cube box. Thus, the average resolution d of the 3D point cloud can be considered as the length of the cube. The calculation shows that the value of d is approximately5cm.IDW interpolation is a frequently used technique for LiDAR data [23]. The core idea of this technique is that nearby points are more alike than distant ones are. The weights assigned to the points closer to the prediction location are greater than those far away. It is calculated by using the following equation:(12)h^=∑i=1n1/dihi∑i=1n1/diwhere hi is the known height value and di is the distance to the known points.According to the precision of the denoised point cloud, a point to be interpolated is selected every5 cm from the vehicle trajectory. Six nearest neighbours to that point are selected to interpolate using IDW with the weight of 1/di. Thereafter, the road height profile is obtained (Figure 7). The key road points of the vehicle’s trajectory, which not only effectively reduce the number of point clouds but also improve the calculation speed, are thus obtained. These points are vital for the next step.Figure 7
Road height profile using IDW algorithm.


## 5.1. Filtering of Point Cloud

The Point Cloud Library (PCL) is a large-scale open project [22] for 2D/3D image and point cloud processing, the framework of which contains numerous state-of-the-art algorithms, including filtering, feature estimation, surface reconstruction, registration, model fitting, and segmentation. For example, these algorithms can be used to filter outliers from noisy data, stitch 3D point clouds together, segment relevant parts of a scene, extract key points, and compute descriptors to recognise objects in the world on the basis of their geometric appearance and create and visualise surfaces from point clouds.The StatisticalOutlierRemoval filter in PCL is used to filter outliers. Compared with the unused filter, the denoised 3D point cloud of the dirt terrain decreased from 346036 to 322628 (Figure6).Figure 6
Denoised road map of the dirt terrain.


## 5.2. Extraction of Road Height Profile

We can predict the trajectory of a vehicle in advance according to the vehicle’s driving direction. We assume that the vehicle drives in a straight line, and the equation for the hypothetical driving track is shown in the following formula:(8)longitudinalrange∈4m,21mlateralrange=-2mThe 2D road height profile, whose accuracy should correspond with that of the 3D point cloud, is discrete because the 3D point cloud data is discrete. Therefore, we should calculate the accuracy and then extract the road height profile according to the calculated precision.Since the resolution of point cloud is affected by many factors, such as the number of lasers in LiDAR and vehicle speed, the average resolution of point cloud is used in the subsequent calculation.The volume of the denoised 3D point cloudV and the single point occupancy volume v can be written as follows:(9)V=rangex×rangey×rangez(10)rangex=20.9994-4.0003=16.9991rangey=2.9998+4.9935=7.99330rangez=0.1977+0.1477=0.34540(11)v=Vthenumberof3DpointcloudWe regardv as the cube box. Thus, the average resolution d of the 3D point cloud can be considered as the length of the cube. The calculation shows that the value of d is approximately5cm.IDW interpolation is a frequently used technique for LiDAR data [23]. The core idea of this technique is that nearby points are more alike than distant ones are. The weights assigned to the points closer to the prediction location are greater than those far away. It is calculated by using the following equation:(12)h^=∑i=1n1/dihi∑i=1n1/diwhere hi is the known height value and di is the distance to the known points.According to the precision of the denoised point cloud, a point to be interpolated is selected every5 cm from the vehicle trajectory. Six nearest neighbours to that point are selected to interpolate using IDW with the weight of 1/di. Thereafter, the road height profile is obtained (Figure 7). The key road points of the vehicle’s trajectory, which not only effectively reduce the number of point clouds but also improve the calculation speed, are thus obtained. These points are vital for the next step.Figure 7
Road height profile using IDW algorithm.


## 6. Interpolation of Road Height Profile Based on Fractal Theory

The preview sensor must be able to image sudden obstacles such as the high curb stone of standard EU 1340 and DIN 483. For this, the preview sensor has to scan the road with a grid width of2cm. So an interpolation method close to the nature geometry should be selected to improve the accuracy of the road height profile, thereby enabling precise control of the active suspension.Fractal interpolation, which is based on fractal theory, irregularity, roughness, and self-similarity of the dataset itself, uses the iterated function system (IFS) to focus on the retention of the global characteristics of a dataset [24, 25]. Recently, several studies have been including fractal compression and fractal interpolation methods in many digitised terrain reconstruction applications [26, 27]. Fractal interpolation can describe natural phenomena more realistically than can classic interpolation methods [25, 28]. This theory utilises its similarity principle to observe detailed structures hidden in a chaotic appearance as a thinking tool for the understanding of global appearance from local characteristics [28]. Figure 8 shows the difference between traditional and fractal geometries on data simulation. Fractal geometry can express the details of a simulated object, whereas traditional geometry cannot [29]. The calculation of fractal interpolation is described in Section 10 “Methodology.”Figure 8
Data interpolation. (a) Spline interpolation. (b) Fractal interpolation.
(a) (b)


### 6.1. Quantitative Analysis of Road Height Profile’s Fractal Characteristics

The quantitative analysis of fractal characteristics uses the “box covering method” to obtain the dual logarithmic coordinate diagram and to analyse its properties further. The calculation steps are as follows:( 1 ) The square with a length r that covers the fractal curve is used (Figure 7).( 2 ) For different r, counting the number of squares Nr is to cover the fractal curve.( 3 ) The dual logarithmic coordinate diagram of ln⁡Nr-ln⁡1/r is plotted.( 4 ) The correlation coefficient of discrete points in the dual logarithmic coordinate diagram is calculated.Table4 lists the values of ri, Nri, ln⁡1/ri, and ln⁡Nri.Table 4
Coverage statistics.
r i(m) 0.005 0.01 0.02 0.04 0.08 N ( r i ) 71400 18700 5100 1275 426  ln ⁡ 1 / r i 5.2983 4.6052 3.9120 3.2189 2.5257  ln ⁡ 1 / N r i 11.1761 9.8363 8.5370 7.1507 6.0544The data ofln1/ri and ln1/Nri are fitted by using the least squares method. Thus, a dual logarithmic regression line graph is obtained (Figure 9). The regression linear equation is as follows:(13)ln⁡Nri=1.8652·ln⁡1ri+1.2541Figure 9
Double logarithmic map of road height profile based on box covering method.The correlation coefficientR of discrete points ln⁡1/Nri and ln⁡1/ri is as follows:(14)R=99.94%where R=99.94% indicates that the road height profile has strong self-similarity. It also reflects the fractal characteristics of the curve quantitatively. Therefore, the theory and method of fractal interpolation can be used to study Figure 7.


### 6.2. Fractal Interpolation of Road Height Profile

The data obtained from Section5.2 are interpolated by fractal interpolation. Four data points are inserted between each interval, which is equivalent to obtaining a road height profile every1 cm. Figure 10 shows the road height information after fractal interpolation.Figure 10
Road height information after fractal interpolation.


## 6.1. Quantitative Analysis of Road Height Profile’s Fractal Characteristics

The quantitative analysis of fractal characteristics uses the “box covering method” to obtain the dual logarithmic coordinate diagram and to analyse its properties further. The calculation steps are as follows:( 1 ) The square with a length r that covers the fractal curve is used (Figure 7).( 2 ) For different r, counting the number of squares Nr is to cover the fractal curve.( 3 ) The dual logarithmic coordinate diagram of ln⁡Nr-ln⁡1/r is plotted.( 4 ) The correlation coefficient of discrete points in the dual logarithmic coordinate diagram is calculated.Table4 lists the values of ri, Nri, ln⁡1/ri, and ln⁡Nri.Table 4
Coverage statistics.
r i(m) 0.005 0.01 0.02 0.04 0.08 N ( r i ) 71400 18700 5100 1275 426  ln ⁡ 1 / r i 5.2983 4.6052 3.9120 3.2189 2.5257  ln ⁡ 1 / N r i 11.1761 9.8363 8.5370 7.1507 6.0544The data ofln1/ri and ln1/Nri are fitted by using the least squares method. Thus, a dual logarithmic regression line graph is obtained (Figure 9). The regression linear equation is as follows:(13)ln⁡Nri=1.8652·ln⁡1ri+1.2541Figure 9
Double logarithmic map of road height profile based on box covering method.The correlation coefficientR of discrete points ln⁡1/Nri and ln⁡1/ri is as follows:(14)R=99.94%where R=99.94% indicates that the road height profile has strong self-similarity. It also reflects the fractal characteristics of the curve quantitatively. Therefore, the theory and method of fractal interpolation can be used to study Figure 7.


## 6.2. Fractal Interpolation of Road Height Profile

The data obtained from Section5.2 are interpolated by fractal interpolation. Four data points are inserted between each interval, which is equivalent to obtaining a road height profile every1 cm. Figure 10 shows the road height information after fractal interpolation.Figure 10
Road height information after fractal interpolation.


## 7. Roughness Grade Division of Road Height Information

Vehicle suspension is actively controlled using the LQG controller to enable the vehicle to obtain improved ride comfort and handling stability on the basis of road height information. Therefore, suitable LQG controller’s weight coefficients of performance indicators that vary in different road roughness grades should be selected. The International Organisation for Standardisation (ISO 8608, 2016) proposed a methodology for classifying road profiles on the basis of power spectral density [30], which divides road roughness into eight grades (A–H).The road height information must be initially divided. The road surface roughness index unifies the traditional statistical parameter and the fractal dimension and is thus suitable for the grade of the road surface [31]. In [31], the expression of road roughness index is as follows:(15)R∗=Rq1/D-1whereD and Rq are the fractal dimension and root mean square (RMS) error of the road surface, respectively. The road height information in this study is as follows:(16)Rq=15.7037×10-3mD=1.8652R∗=24.1180×10-3maccording to the data provided in [31], which belongs to B grade. Next, the LQG controller of the square active suspension is designed to solve the controller’s weight coefficients under different road roughness grades.


## 8. Optimal Design of LQG Controller

Before designing the LQG controller of the active suspension, a dynamic model of passive and active suspensions of the vehicle, which is the basis of follow-up calculation, must be established.


### 8.1. Suspension System Dynamics

Figure11 establishes a quarter vehicle model with passive and active suspensions, which have two degrees of freedom and heave dynamics. Unlike passive suspension, active suspension has only one input of road excitation and more active controlling force; moreover, a force actuator is installed in place of the damper variable. In this implementation, the tire is represented by a linear spring. In Table 5, the parameters of the quarter vehicle model are shown. Since the active suspension is a parallel structure of actuator and spring, the passive suspension replaces actuator with damper, so the suspension stiffness of both is set to different values.Table 5
Quarter vehicle model parameters.
Symbol Quantity Value Unit m s Sprung mass 320 k g  m u Unsprung mass 40 k g  Ks1 Passive suspension stiffness 22000 N / m  K s 2 Active suspension stiffness 20000 N / m  K t Tire stiffness 200000 N / m  C s Passive suspension damping 1000 N ⋅ s / m  Z s 1 Sprung mass’s displacements of passive suspension - m  Z s 2 Sprung mass’s displacements of active suspension - m  Z u 1 Unsprung mass’s displacements of passive suspension - m  Z u 2 Unsprung mass’s displacements of active suspension - m  q Displacements of road excitation - m  u Active force - NFigure 11
Quarter vehicle model. (a) Passive suspension. (b) Active suspension.
(a) (b)According to the quarter vehicle model based on Newton’s second law, the equations of motion of the system are as follows [32]:( 1 ) Passive suspension:(17)muZ¨u1=Ktq-Zu1-Ks1Zu1-Zs1-CSZ˙u1-Z˙s1msZ¨s1=Ks1Zu1-Zs1+CSZ˙u1-Z˙s1( 2 ) Active suspension:(18)muZ¨u2=Ktq-Zu2-Ks1Zu2-Zs2-umsZ¨s2=Ks1Zu2-Zs2+uThe state equation is established, and the system state variableX is defined as follows:(19)X1=Z˙s1Z˙u1Zs1Zu1TX2=Z˙s2Z˙u2Zs2Zu2TThe suspension system can be expressed in the following state equation:(20)X˙1=A1·X1+B1·qX˙2=A2·X2+B2·UU=uqwhere(21)
A
1
=
-
C
s
m
s
C
s
m
s
-
K
s
1
m
s
K
s
1
m
s
C
s
m
u
-
C
s
m
u
K
s
1
m
u
-
K
s
1
+
K
t
m
u
1
0
0
0
0
1
0
0
B
1
=
0
K
t
m
u
0
0
A
2
=
0
0
-
K
s
2
m
s
K
s
2
m
s
0
0
K
s
2
m
u
-
K
s
2
+
K
t
m
u
1
0
0
0
0
1
0
0
B
2
=
1
m
s
0
-
1
m
u
K
t
m
u
0
0
0
0The system output Y is defined as follows:(22)Y1=Z¨s1Zu1-Zs1q-Zu1TY2=Z¨s2Zu2-Zs2q-Zu2TThe output equation of the system is as follows:(23)Y1=C1·X+D1·qY2=C2·X+D2·Uwhere (24)C1=-CsmsCsms-Ks1msKs1ms00-11000-1D1=001C2=00-Ks2msKs1ms00-11000-1D2=1ms00001


### 8.2. LQG Controller

The main performance evaluation indicators in designing the controller for the quarter active suspension are as follows:(1) vehicle body acceleration, which evaluates ride comfort and handing stability; (2) suspension dynamic deflection, which affects body posture; and (3) tire dynamic deflection on behalf of the tire road holding [33].Here, the comprehensive performance evaluation indicator denoted byJ is defined as follows [34, 35]:(25)J=limT→∞⁡1Tq1·Zu2t-q2+q2·Zs2t-Zu2t2+q3·Z¨s22tdtwhere q1,q2,q3 are the weight coefficients of tire dynamic deflection, suspension dynamic deflection, and vehicle body acceleration, respectively. The weight coefficient q3 of the vehicle body acceleration, which is used as a base value for the other coefficients to simplify the calculation, is set to 1 because it is important for vehicle body acceleration.Then, formula (25) can be rewritten as an integral type of quadratic function by the following:(26)J=limT→∞1TX2TQX2+UTRU+2X2TNUdtwhere(27)Q=0000000000q2+q3Ks22ms2-q2-q3Ks22ms200-q2-q3Ks22ms2q1+q2+q3Ks22ms2N=0000-Ks2ms20Ks2ms2-q1R=q3ms200q1When the vehicle parameter values and weight coefficientsq1,q2,q3 are determined, the optimal control feedback gain matrix K can be obtained from the following Riccati equation:(28)PA2+A2TP-PB2+NR-1B2TP+NT+Q=0(29)K=R-1B2TP+NTThen, according to the state variablesX2 at any time, the optimal control force of the actuator denoted by u can be described as follows:(30)u=-K·X2The solving method of the weight coefficients in different road grades is discussed in the succeeding section.


## 8.1. Suspension System Dynamics

Figure11 establishes a quarter vehicle model with passive and active suspensions, which have two degrees of freedom and heave dynamics. Unlike passive suspension, active suspension has only one input of road excitation and more active controlling force; moreover, a force actuator is installed in place of the damper variable. In this implementation, the tire is represented by a linear spring. In Table 5, the parameters of the quarter vehicle model are shown. Since the active suspension is a parallel structure of actuator and spring, the passive suspension replaces actuator with damper, so the suspension stiffness of both is set to different values.Table 5
Quarter vehicle model parameters.
Symbol Quantity Value Unit m s Sprung mass 320 k g  m u Unsprung mass 40 k g  Ks1 Passive suspension stiffness 22000 N / m  K s 2 Active suspension stiffness 20000 N / m  K t Tire stiffness 200000 N / m  C s Passive suspension damping 1000 N ⋅ s / m  Z s 1 Sprung mass’s displacements of passive suspension - m  Z s 2 Sprung mass’s displacements of active suspension - m  Z u 1 Unsprung mass’s displacements of passive suspension - m  Z u 2 Unsprung mass’s displacements of active suspension - m  q Displacements of road excitation - m  u Active force - NFigure 11
Quarter vehicle model. (a) Passive suspension. (b) Active suspension.
(a) (b)According to the quarter vehicle model based on Newton’s second law, the equations of motion of the system are as follows [32]:( 1 ) Passive suspension:(17)muZ¨u1=Ktq-Zu1-Ks1Zu1-Zs1-CSZ˙u1-Z˙s1msZ¨s1=Ks1Zu1-Zs1+CSZ˙u1-Z˙s1( 2 ) Active suspension:(18)muZ¨u2=Ktq-Zu2-Ks1Zu2-Zs2-umsZ¨s2=Ks1Zu2-Zs2+uThe state equation is established, and the system state variableX is defined as follows:(19)X1=Z˙s1Z˙u1Zs1Zu1TX2=Z˙s2Z˙u2Zs2Zu2TThe suspension system can be expressed in the following state equation:(20)X˙1=A1·X1+B1·qX˙2=A2·X2+B2·UU=uqwhere(21)
A
1
=
-
C
s
m
s
C
s
m
s
-
K
s
1
m
s
K
s
1
m
s
C
s
m
u
-
C
s
m
u
K
s
1
m
u
-
K
s
1
+
K
t
m
u
1
0
0
0
0
1
0
0
B
1
=
0
K
t
m
u
0
0
A
2
=
0
0
-
K
s
2
m
s
K
s
2
m
s
0
0
K
s
2
m
u
-
K
s
2
+
K
t
m
u
1
0
0
0
0
1
0
0
B
2
=
1
m
s
0
-
1
m
u
K
t
m
u
0
0
0
0The system output Y is defined as follows:(22)Y1=Z¨s1Zu1-Zs1q-Zu1TY2=Z¨s2Zu2-Zs2q-Zu2TThe output equation of the system is as follows:(23)Y1=C1·X+D1·qY2=C2·X+D2·Uwhere (24)C1=-CsmsCsms-Ks1msKs1ms00-11000-1D1=001C2=00-Ks2msKs1ms00-11000-1D2=1ms00001


## 8.2. LQG Controller

The main performance evaluation indicators in designing the controller for the quarter active suspension are as follows:(1) vehicle body acceleration, which evaluates ride comfort and handing stability; (2) suspension dynamic deflection, which affects body posture; and (3) tire dynamic deflection on behalf of the tire road holding [33].Here, the comprehensive performance evaluation indicator denoted byJ is defined as follows [34, 35]:(25)J=limT→∞⁡1Tq1·Zu2t-q2+q2·Zs2t-Zu2t2+q3·Z¨s22tdtwhere q1,q2,q3 are the weight coefficients of tire dynamic deflection, suspension dynamic deflection, and vehicle body acceleration, respectively. The weight coefficient q3 of the vehicle body acceleration, which is used as a base value for the other coefficients to simplify the calculation, is set to 1 because it is important for vehicle body acceleration.Then, formula (25) can be rewritten as an integral type of quadratic function by the following:(26)J=limT→∞1TX2TQX2+UTRU+2X2TNUdtwhere(27)Q=0000000000q2+q3Ks22ms2-q2-q3Ks22ms200-q2-q3Ks22ms2q1+q2+q3Ks22ms2N=0000-Ks2ms20Ks2ms2-q1R=q3ms200q1When the vehicle parameter values and weight coefficientsq1,q2,q3 are determined, the optimal control feedback gain matrix K can be obtained from the following Riccati equation:(28)PA2+A2TP-PB2+NR-1B2TP+NT+Q=0(29)K=R-1B2TP+NTThen, according to the state variablesX2 at any time, the optimal control force of the actuator denoted by u can be described as follows:(30)u=-K·X2The solving method of the weight coefficients in different road grades is discussed in the succeeding section.


## 9. Solution of Weight Coefficients

We need to solve eight sets of weight coefficients to enable the active suspension to adapt to different road grades. According to [36], random road excitation can be written as follows:(31)q˙t=-2πn00vqt+2πn0Sqn0vWtwhere n00=0.011m-1 is the spatial lower cut-off frequency. Wt stands for Gaussian white noise with zero mean. v represents the speed of the vehicle.According to formula (31), under the speed 40 km/h, we obtain eight grades of road surface roughness, as shown in Figure 12.Figure 12
Road surface roughness of different grades.AHP is a multicriteria decision-aiding method that is based on a solid axiomatic foundation. It is a systematic procedure for dealing with complex decision-making problems, in which many competing alternatives exist [37–39]. Thus, we adopt this method to solve the weight coefficients of the LQG controller. The process of calculating LQG controller coefficients using AHP method is described in Section 10 “Methodology.”Table6 lists the passive suspension performance statistics under different road grades.Table 6
Performance statistics of passive suspension under different road grades.
Road grades σ B A 1 m / s 2 σ D T D 1 m σ S W S 1 m A 0.3418 0.0011 0.0034  B 0.6836 0.0023 0.0069  C 1.3671 0.0046 0.0137  D 2.7342 0.0091 0.0275  E 5.4684 0.0182 0.0550  F 10.9368 0.0365 0.1099  G 21.8737 0.0729 0.2198  H 43.7473 0.1458 0.4396On the high-grade road surface, the ride comfort is mainly optimized, and on the poor grade road surface, the handling stability is mainly optimized. The subjective judgment matrices of different road grades are as follows:(32)HA=129121319131HB=139131319131HC=1348431418141HD=1346431416141HE=1131111311HF=1121111211HG=1111115151HH=1111116161Table7 shows the calculated LQG controller’s weight coefficients in different road roughness grades.Table 7
LQG controller weighting factor in different road grades.
Road grades q 1 q 2 A (qA) 39323.4009 1259.5536  B (qB) 30009.3714 1100.3214  C (qC) 86562.0060 1716.5706  D (qD) 95273.8155 2079.4772  E (qE) 62422.0080 4760.8143  F (qF) 71455.3616 6238.4316  G (qG) 52648.7606 16933.7086  H (qH) 49544.3806 17994.7505Because the road height profile calculated in Section7 belongs to B grade, the weight coefficients in the LQG controller are selected as follows:(33)q1=30009.3714q2=1100.3214


## 10. Methodology

### 10.1. Fractal Interpolation

#### 10.1.1. Affine Transformation

A two-dimensional affine transformationW is formed as follows:(34)Wx′y′=abcdxy+efwhere a,b,c,d,e,f are real numbers, px,y is a point on the plane, and W is a linear affine transformation.The transformation matrixA can be decomposed into rotation R, expansion L, and distortion T.(35)A=abcd=R·T·LR=cos⁡θ-sin⁡θsin⁡θcos⁡θT=1tan⁡αtan⁡α1L=l100l2Figure13 shows a decomposition diagram of an affine transformation.Figure 13
Affine transform decomposition diagram.


#### 10.1.2. IFS and Fractal Interpolation

A set of data pointsPi=xi,yi∈R2∣i=0,1,2,⋯,N satisfies x0<x1<x2<⋯<xN. The attractor of IFS is a continuous function of interpolating data. An interpolation curve across a series of data points Pi and map Po and PN to Pn-1 and Pn is formed.Considering IFSR2;Wn,n=0,1,2,⋯,N, Wn is an affine transformation with the following form:(36)Wnxy=an0cndn·xy+enfnThe aforementioned equation satisfies the following:(37)Wnx0y0=xn-1yn-1WnxNyN=xnynwhere bn=0 guarantees that the interpolation function of each cell does not overlap. Assume that dn is constant; then we can solve the parameters as follows [38]:(38)an=xn-xn-1xN-x0en=xNxn-1-x0xnxN-x0cn=yn-yn-1-dnyN-y0xN-x0fn=xNyn-1-x0yn-dnxNy0-x0yNxN-x0For a set of data pointsPi=xi,yi∈R2∣i=0,1,2,⋯,N, N numbers of affine transformations exist. The fractal interpolation function of the affine transformation of nth is as follows:(39)xni=anxi+enyni=cnxi+dnyi+fni=0,1,···,N;n=1,2,⋯,N


### 10.2. Analytic Hierarchy Process

#### 10.2.1. Calculating Objective Weight Coefficient

According to the parameters of a passive suspension, its dynamic simulation model is operated under eight road grades to obtain the performance statisticsσBA1,σDTD1,σSWS1, which are the RMS of vehicle’s vertical acceleration, tire dynamic deflection, and suspension dynamic deflection, respectively. The objective weight coefficient of vertical acceleration is supposed to be 1. According to [38], the equation of calculating other objective weight coefficients is as follows:(40)σBA12·1=β1·σDTD12=β2·σSWS12


#### 10.2.2. Determining Subjective Weight Coefficients

( 1 ) A subjective judgement matrix is created as follows.h i j is the comparison value of the importance of the indicators i and j, and a fundamental 1-to-9 scale can be used to rank the judgments (Table 8).Table 8
Comparison values of importance between indicators.
i / j Equally Moderately Strongly Very Extremely h i j 1 3 5 7 9The subjective judgement matrixH is built through pairwise comparison of each decision factor.(41)H=1h12⋯⋯h1n1h121h23⋯h2n⋯1h23⋯⋯⋯⋮⋮⋱⋱⋮1h1n1h2n⋯⋯1( 2 ) Matrix H is calculated by multiplying the vector of every row as follows:(42)M=M1,M2,⋯,MnTMi=∏j=1nhiji,j=1,2,⋯,n( 3 )   M n is calculated as follows:(43)W¯=W1¯,W2¯,⋯,Wn¯TWi¯=Mini=1,2,⋯,n( 4 ) Positive vector W is calculated as follows:(44)W=W¯∑i=1nWi¯( 5 ) The maximum eigenvalue of matrix H is calculated and checked for consistency as follows.If every element in matrixH satisfies the equation hij=1/hji and hij=hik∙hkj, then the matrix is the consistency matrix. The subjective judgement matrix is often not perfectly consistent due to the randomness of people’s judgment. These judgment errors can be detected by the consistency ratio CR.(45)λmax=∑i=1nHWinWii=1,2,⋯,n(46)CR=λmax-nRIn-1where RI is a random consistency index of matrix H.Whenn is 10, RI is equal to 1.4851. The solution is correct when CR is less than 1. If CR is more than 1, then the evaluation matrix H needs to be revised on the basis of [39].( 6 ) Subjective weight coefficients.The subjective weight coefficient of the vehicle’s vertical acceleration is 1. The subjective weight coefficients of the tire dynamic deflectionγ1 and suspension dynamic deflection γ2 are determined according to the following formula:(47)W1=W2γ1=W3γ2


#### 10.2.3. Calculating LQG Controller’s Weight Coefficientsq1 and q2

The general weights of the evaluating indexes related to ride comfort can be obtained on the basis of the objective weights and subjective weight coefficients.(48)q1=β1·γ1q2=β2·γ2


## 10.1. Fractal Interpolation

### 10.1.1. Affine Transformation

A two-dimensional affine transformationW is formed as follows:(34)Wx′y′=abcdxy+efwhere a,b,c,d,e,f are real numbers, px,y is a point on the plane, and W is a linear affine transformation.The transformation matrixA can be decomposed into rotation R, expansion L, and distortion T.(35)A=abcd=R·T·LR=cos⁡θ-sin⁡θsin⁡θcos⁡θT=1tan⁡αtan⁡α1L=l100l2Figure13 shows a decomposition diagram of an affine transformation.Figure 13
Affine transform decomposition diagram.


### 10.1.2. IFS and Fractal Interpolation

A set of data pointsPi=xi,yi∈R2∣i=0,1,2,⋯,N satisfies x0<x1<x2<⋯<xN. The attractor of IFS is a continuous function of interpolating data. An interpolation curve across a series of data points Pi and map Po and PN to Pn-1 and Pn is formed.Considering IFSR2;Wn,n=0,1,2,⋯,N, Wn is an affine transformation with the following form:(36)Wnxy=an0cndn·xy+enfnThe aforementioned equation satisfies the following:(37)Wnx0y0=xn-1yn-1WnxNyN=xnynwhere bn=0 guarantees that the interpolation function of each cell does not overlap. Assume that dn is constant; then we can solve the parameters as follows [38]:(38)an=xn-xn-1xN-x0en=xNxn-1-x0xnxN-x0cn=yn-yn-1-dnyN-y0xN-x0fn=xNyn-1-x0yn-dnxNy0-x0yNxN-x0For a set of data pointsPi=xi,yi∈R2∣i=0,1,2,⋯,N, N numbers of affine transformations exist. The fractal interpolation function of the affine transformation of nth is as follows:(39)xni=anxi+enyni=cnxi+dnyi+fni=0,1,···,N;n=1,2,⋯,N


## 10.1.1. Affine Transformation

A two-dimensional affine transformationW is formed as follows:(34)Wx′y′=abcdxy+efwhere a,b,c,d,e,f are real numbers, px,y is a point on the plane, and W is a linear affine transformation.The transformation matrixA can be decomposed into rotation R, expansion L, and distortion T.(35)A=abcd=R·T·LR=cos⁡θ-sin⁡θsin⁡θcos⁡θT=1tan⁡αtan⁡α1L=l100l2Figure13 shows a decomposition diagram of an affine transformation.Figure 13
Affine transform decomposition diagram.


## 10.1.2. IFS and Fractal Interpolation

A set of data pointsPi=xi,yi∈R2∣i=0,1,2,⋯,N satisfies x0<x1<x2<⋯<xN. The attractor of IFS is a continuous function of interpolating data. An interpolation curve across a series of data points Pi and map Po and PN to Pn-1 and Pn is formed.Considering IFSR2;Wn,n=0,1,2,⋯,N, Wn is an affine transformation with the following form:(36)Wnxy=an0cndn·xy+enfnThe aforementioned equation satisfies the following:(37)Wnx0y0=xn-1yn-1WnxNyN=xnynwhere bn=0 guarantees that the interpolation function of each cell does not overlap. Assume that dn is constant; then we can solve the parameters as follows [38]:(38)an=xn-xn-1xN-x0en=xNxn-1-x0xnxN-x0cn=yn-yn-1-dnyN-y0xN-x0fn=xNyn-1-x0yn-dnxNy0-x0yNxN-x0For a set of data pointsPi=xi,yi∈R2∣i=0,1,2,⋯,N, N numbers of affine transformations exist. The fractal interpolation function of the affine transformation of nth is as follows:(39)xni=anxi+enyni=cnxi+dnyi+fni=0,1,···,N;n=1,2,⋯,N


## 10.2. Analytic Hierarchy Process

### 10.2.1. Calculating Objective Weight Coefficient

According to the parameters of a passive suspension, its dynamic simulation model is operated under eight road grades to obtain the performance statisticsσBA1,σDTD1,σSWS1, which are the RMS of vehicle’s vertical acceleration, tire dynamic deflection, and suspension dynamic deflection, respectively. The objective weight coefficient of vertical acceleration is supposed to be 1. According to [38], the equation of calculating other objective weight coefficients is as follows:(40)σBA12·1=β1·σDTD12=β2·σSWS12


### 10.2.2. Determining Subjective Weight Coefficients

( 1 ) A subjective judgement matrix is created as follows.h i j is the comparison value of the importance of the indicators i and j, and a fundamental 1-to-9 scale can be used to rank the judgments (Table 8).Table 8
Comparison values of importance between indicators.
i / j Equally Moderately Strongly Very Extremely h i j 1 3 5 7 9The subjective judgement matrixH is built through pairwise comparison of each decision factor.(41)H=1h12⋯⋯h1n1h121h23⋯h2n⋯1h23⋯⋯⋯⋮⋮⋱⋱⋮1h1n1h2n⋯⋯1( 2 ) Matrix H is calculated by multiplying the vector of every row as follows:(42)M=M1,M2,⋯,MnTMi=∏j=1nhiji,j=1,2,⋯,n( 3 )   M n is calculated as follows:(43)W¯=W1¯,W2¯,⋯,Wn¯TWi¯=Mini=1,2,⋯,n( 4 ) Positive vector W is calculated as follows:(44)W=W¯∑i=1nWi¯( 5 ) The maximum eigenvalue of matrix H is calculated and checked for consistency as follows.If every element in matrixH satisfies the equation hij=1/hji and hij=hik∙hkj, then the matrix is the consistency matrix. The subjective judgement matrix is often not perfectly consistent due to the randomness of people’s judgment. These judgment errors can be detected by the consistency ratio CR.(45)λmax=∑i=1nHWinWii=1,2,⋯,n(46)CR=λmax-nRIn-1where RI is a random consistency index of matrix H.Whenn is 10, RI is equal to 1.4851. The solution is correct when CR is less than 1. If CR is more than 1, then the evaluation matrix H needs to be revised on the basis of [39].( 6 ) Subjective weight coefficients.The subjective weight coefficient of the vehicle’s vertical acceleration is 1. The subjective weight coefficients of the tire dynamic deflectionγ1 and suspension dynamic deflection γ2 are determined according to the following formula:(47)W1=W2γ1=W3γ2


### 10.2.3. Calculating LQG Controller’s Weight Coefficientsq1 and q2

The general weights of the evaluating indexes related to ride comfort can be obtained on the basis of the objective weights and subjective weight coefficients.(48)q1=β1·γ1q2=β2·γ2


## 10.2.1. Calculating Objective Weight Coefficient

According to the parameters of a passive suspension, its dynamic simulation model is operated under eight road grades to obtain the performance statisticsσBA1,σDTD1,σSWS1, which are the RMS of vehicle’s vertical acceleration, tire dynamic deflection, and suspension dynamic deflection, respectively. The objective weight coefficient of vertical acceleration is supposed to be 1. According to [38], the equation of calculating other objective weight coefficients is as follows:(40)σBA12·1=β1·σDTD12=β2·σSWS12


## 10.2.2. Determining Subjective Weight Coefficients

( 1 ) A subjective judgement matrix is created as follows.h i j is the comparison value of the importance of the indicators i and j, and a fundamental 1-to-9 scale can be used to rank the judgments (Table 8).Table 8
Comparison values of importance between indicators.
i / j Equally Moderately Strongly Very Extremely h i j 1 3 5 7 9The subjective judgement matrixH is built through pairwise comparison of each decision factor.(41)H=1h12⋯⋯h1n1h121h23⋯h2n⋯1h23⋯⋯⋯⋮⋮⋱⋱⋮1h1n1h2n⋯⋯1( 2 ) Matrix H is calculated by multiplying the vector of every row as follows:(42)M=M1,M2,⋯,MnTMi=∏j=1nhiji,j=1,2,⋯,n( 3 )   M n is calculated as follows:(43)W¯=W1¯,W2¯,⋯,Wn¯TWi¯=Mini=1,2,⋯,n( 4 ) Positive vector W is calculated as follows:(44)W=W¯∑i=1nWi¯( 5 ) The maximum eigenvalue of matrix H is calculated and checked for consistency as follows.If every element in matrixH satisfies the equation hij=1/hji and hij=hik∙hkj, then the matrix is the consistency matrix. The subjective judgement matrix is often not perfectly consistent due to the randomness of people’s judgment. These judgment errors can be detected by the consistency ratio CR.(45)λmax=∑i=1nHWinWii=1,2,⋯,n(46)CR=λmax-nRIn-1where RI is a random consistency index of matrix H.Whenn is 10, RI is equal to 1.4851. The solution is correct when CR is less than 1. If CR is more than 1, then the evaluation matrix H needs to be revised on the basis of [39].( 6 ) Subjective weight coefficients.The subjective weight coefficient of the vehicle’s vertical acceleration is 1. The subjective weight coefficients of the tire dynamic deflectionγ1 and suspension dynamic deflection γ2 are determined according to the following formula:(47)W1=W2γ1=W3γ2


## 10.2.3. Calculating LQG Controller’s Weight Coefficientsq1 and q2

The general weights of the evaluating indexes related to ride comfort can be obtained on the basis of the objective weights and subjective weight coefficients.(48)q1=β1·γ1q2=β2·γ2


## 11. Experiment

When the vehicle’s speed is constant, the more point cloud data are collected and the higher data fusion rate is required. Therefore, in order to meet the needs of high-speed data processing of active suspension, firstly, only collect and store the topographic data in front of the vehicle; secondly, collect the data by using the suitable number of lines of LiDAR; thirdly, the algorithm adopted in this paper not only effectively collects enough point cloud data on the ground but also satisfies the accuracy of the solution after data fusion and obtains the elevation information of the road surface for effective control of suspension.Road experiments of the prototype vehicle are performed to verify the accuracy and efficiency of the active suspension control algorithm that is based on LiDAR-previewed terrain information. The main body of the prototype vehicle is assembled and constructed using the FAW Besturn car chassis (Figure14). It includes the mechanical, hydraulic servo, and the electric control parts. The single active suspension actuator unit is shown in Figure 15, and the left side is the servo control actuator cylinder. Figure 4 shows the test site.Figure 14
Prototype vehicle with active suspension.Figure 15
Single active suspension actuator.Figures16–22 show the experiment data. In order to properly validate selection of control weightings based on estimated road profile and road class, the performance evaluation indicators of active control for different weight factors q1 and q2 that are under different road class are reported in Figures 16–18. Figures 19–21 depict the performance evaluation indicator comparison between active suspension with LQG controller whose controller coefficient is selected as qB and passive suspension. Figure 22 shows changes in the control force of the active suspension actuators.Figure 16
Vertical acceleration of active suspension.Figure 17
Tire dynamic deflection of active suspension.Figure 18
Active suspension dynamic deflection.Figure 19
Vertical acceleration.Figure 20
Tire dynamic deflection.Figure 21
Suspension dynamic deflection.Figure 22
Active control force.Table9 shows the RMS of performance statistics of LQG-active suspension system based on different weight coefficients q1 and q2 and passive suspension system. Draw the data in Table 9 as a bar chart as shown in Figures 23–25. As can be seen from Figure 23, the RMS of vertical acceleration of LQG active suspension controller based on grade B pavement coefficients qB is about 19.46% lower than that of passive suspension and about 10.33%, 10.09%, 9.77%, 8.32%, 7.09%, 3.54%, and 0.78% lower than the LQG active suspension controller based on other pavement grade coefficients qF, qG, qE, qH, qD, qA, and qC.Table 9
RMS values of active and passive suspension performance indicators.
Active suspension Passive suspension q A q B q C q D q E q F q G q H Vertical acceleration (m/s2) 0.8581 0.8277 0.8342 0.8909 0.9173 0.9231 0.9206 0.9028 1.0277  Tire deflection (m) 0.0246 0.0242 0.0258 0.0247 0.0251 0.0256 0.0248 0.0245 0.0270  Suspension deflection (m) 0.0344 0.0335 0.0337 0.0342 0.0338 0.0339 0.0336 0.0343 0.0368Figure 23
Bar-chart of suspension’s vertical acceleration.Figure 24
Bar-chart of suspension’s tire dynamic deflection.Figure 25
Bar-chart of suspension’s dynamic deflection.In Figure24, the RMS of tire dynamic deflection of LQG active suspension controller based on grade B pavement coefficients qB is about 15.9% lower than that of passive suspension and about 10.33%, 10.09%, 9.77%, 8.32%, 7.09%, 3.54%, and 0.78% lower than the LQG active suspension controller based on other pavement grade coefficients qF, qG, qE, qH, qD, qA, and qC.In Figure25, the RMS of suspension dynamic deflection of LQG active suspension controller based on grade B pavement coefficients qB is about 8.97% lower than that of passive suspension and about 2.62%, 2.33%, 1.46%, 1.18%, 0.89%, 0.59%, and 0.29% lower than the LQG active suspension controller based on other pavement grade coefficients qA, qH, qD, qF, qE, qC, and qG.The experimental data show that reasonable selection of LQG controller weightings based on estimated road profile and road class through fractal interpolation can improve vehicle ride comfort and handling stability effectively more than passive suspension does.


## 12. Conclusions

This study focuses on extracting a road height profile from the 3D point cloud data of the terrain in front of a vehicle as input to the suspension control model. Three innovations are mainly presented. Firstly, a LiDAR coupled with INS and GPS is used to measure the 3D point cloud of the terrain in front of the vehicle. Although many researchers have proposed the use of preview information for active suspension, they did not provide detailed descriptions and instead only presented concepts. Secondly, the road height profile is extracted using the IDW interpolation method. The main difference is that this method has been used for 3D digital terrain obtained by airborne LiDAR systems previously and is rarely used in a vehicle’s LiDAR system to obtain preview information of a road. Thirdly, fractal interpolation is used to improve the accuracy of road height information. This interpolation method is close to natural geometry. Thus, the interpolated data in this work are more precise than traditionally interpolated data.Road experiment verifies that reasonable selection of LQG controller weightings based on estimated road profile and road class through fractal interpolation can improve the ride comfort and handling stability of the vehicle more than passive suspension does.



---

*Source: 3483710-2019-06-27.xml*