Research on LQR control of nonlinear active suspension cooperatively optimized by differential geometry and genetic algorithms

Abstract

With the continuous advancement of the automotive industry’s demand for driving comfort and handling stability, active suspension systems have become a research hotspot due to their ability to achieve effective balancing among multiple performance metrics. However, traditional active suspension control algorithms face challenges such as insufficient modeling accuracy and excessive reliance on empirical parameter optimization when dealing with nonlinear systems. To address these issues, this study proposes an improved linear quadratic regulator (LQR) control strategy that integrates differential geometry theory and genetic algorithm (GA) optimization to enhance the control performance of nonlinear active suspension systems. First, a quarter-car nonlinear active suspension model is established using differential geometry theory, effectively preserving the system’s nonlinear characteristics and improving modeling precision. Building on this foundation, an enhanced LQR controller is designed, combined with the global optimization capability of the genetic algorithm, to construct a fitness function based on multiple suspension performance metrics. This approach achieves automated optimization of the weight coefficient matrix and multi-objective balancing. To validate the effectiveness of the proposed control strategy, comparative simulations under discrete impact road and Class C random road conditions are conducted using the MATLAB-Simulink platform. Results demonstrate that, compared to PID control, the GA-LQR strategy reduces the root mean square (RMS) values of body vertical acceleration, suspension dynamic deflection, and tire dynamic displacement by 6.94%, 57.72%, and 4.4%, respectively, under discrete impact road conditions. Under Class C road conditions, these reductions further increase to 22.1%, 63.6%, and 16.25%. This study provides a high-performance and robust control strategy for nonlinear suspension systems, laying a theoretical foundation for improving overall vehicle performance and driving handling stability.

Introduction

With the continuous development of the automotive industry and the gradual improvement of people’s living standards, ride comfort has become one of the key indicators for measuring vehicle quality. As a crucial component of the vehicle chassis, the suspension system primarily consists of elastic elements, guidance mechanisms, and shock absorbers. Its main function is to mitigate shocks experienced by the vehicle when driving on uneven roads, thereby enhancing ride comfort. Based on working principles and adjustment methods, vehicle suspension systems can be categorized into three types: passive suspension, semi-active suspension, and active suspension. Passive suspension systems have fixed damping coefficients, which imposes certain limitations on balancing vehicle handling performance and ride smoothness. To improve vehicle maneuverability while effectively suppressing vertical, roll, and pitch motions caused by road impacts, researchers have developed semi-active suspension systems with adjustable damping characteristics and active suspension systems capable of proactively generating actuation forces to control body movements. Compared to passive and semi-active systems, active suspensions demonstrate superior ability to balance conflicting performance metrics such as ride comfort, suspension travel, and road grip, garnering significant attention in the industry.^1–5

The performance of active suspension systems is directly determined by the quality of their control algorithms. ⁶ Current mainstream control algorithms can generally be divided into two categories: model-based control algorithms and model-free control algorithms. Model-free control algorithms, which do not rely on precise system modeling but instead achieve feedback control through input-output data, demonstrate unique advantages when handling nonlinear controlled systems. Common model-free control methods include fuzzy control^7–9 and PID control.^10,11 Researchers like Sun and Yang developed an adjustable fuzzy control algorithm for suspension systems, effectively improving vehicle comfort and driving stability. ⁷ Rashid et al. proposed a quarter-car semi-active suspension control method incorporating fuzzy logic controllers, leveraging the broad adjustment capabilities of magnetorheological (MR) dampers to reduce body vibrations. ⁸ Laoufi and Eskandarian established a seven-degree-of-freedom nonlinear full-vehicle dynamics model and implemented fuzzy controllers, with simulations confirming enhanced ride smoothness. ⁹ Han et al. introduced an adaptive fuzzy PID control strategy enabling self-adjustment of control parameters under complex road conditions. ¹⁰ However, both fuzzy control and PID control often require extensive trial-and-error processes and heavy reliance on operator expertise when formulating control rules or tuning gain coefficients. Such experience-driven parameter adjustment methods struggle to meet modern requirements for high precision and automation in control algorithms.

To achieve optimal suspension control performance, model-based suspension control algorithms have attracted extensive attention. Representative approaches include the Linear Quadratic Regulator (LQR), ¹² adaptive control,^13–15 and Model Predictive Control (MPC).^16,17 Among these, the LQR algorithm has been widely adopted in industrial applications for stabilizing linear time-invariant systems due to its superior stability and robustness. ¹⁸ To address system parameter uncertainties and external stochastic disturbances, Bai and Wang ¹⁹ proposed a robust optimal integral sliding mode control strategy by integrating LQR optimal control with sliding mode control. For tractor systems with significant unsprung mass, Ahn et al. ²⁰ designed a Linear Quadratic Gaussian (LQG) optimal control algorithm. Jin and Luo ²¹ employed an equivalent linearization method to transform nonlinear suspension systems into linear ones by minimizing the mean square error to determine equivalent damping and stiffness matrices. However, existing studies often oversimplify suspension system modeling by treating inherently nonlinear parameters such as spring stiffness coefficients and damping coefficients as fixed values.^22–24 In practice, suspension systems exhibit significant complexity, and neglecting the nonlinear characteristics of stiffness and damping leads to substantial deviations between predicted and actual results. ¹⁴ Therefore, comprehensive consideration of these nonlinear factors is essential during suspension system modeling to ensure model accuracy and control strategy. ⁷

In recent years, research on transforming complex nonlinear systems into linear models through differential geometry theory has attracted significant attention. Compared with conventional linearization methods, the differential geometry approach accounts for the geometric structure of nonlinear systems, thereby better preserving their inherent nonlinear characteristics during transformation and yielding linearized models with enhanced accuracy.^25,26 Chen et al. ²⁷ decoupled a nonlinear suspension system using differential geometry methods, isolating the wheel and sprung mass as independent linear subsystems. To achieve precise control of a magnetically suspended double-gimbal control moment gyroscope (MSDGCMG), Chen and Chen ²⁸ proposed a novel linearization methodology based on differential geometry theory, significantly improving the linearization precision of MSDGCMG. Pan et al. ²⁹ employed differential geometric linearization theory to convert a Static Var Generator (SVG) nonlinear system into a linear equivalent, subsequently developing an SVG control algorithm that demonstrated superior robustness and stability in simulations. Rezapour et al. ³⁰ derived robotic motion equations using differential geometric techniques and implemented partial feedback linearization to linearize the direct-drive portion of system dynamics, with simulation results validating their theoretical framework. These research advancements collectively demonstrate the distinctive advantages of differential geometry theory in addressing system nonlinearities.

The integration of differential geometry theory enhances the precision of suspension system models, which in turn improves the performance of Linear Quadratic Regulator (LQR) controllers. However, this advancement simultaneously increases controller complexity, particularly in determining the weighting matrices Q and R within the LQR framework. Conventionally, these weighting matrices have been predominantly determined through empirical engineering practices—an approach plagued by inherent subjectivity and inefficiency. ¹⁸ Furthermore, the lack of direct correspondence between control performance optimization and the configuration of LQR weighting matrices exacerbates the complexity and difficulty of parameter tuning.

The emergence of intelligent optimization algorithms has introduced novel solutions to controller parameter calibration. Kumar ³¹ proposed an Adaptive Particle Swarm Optimization (APSO) algorithm to optimize weighting matrix elements in LQR controllers, demonstrating its viability as an effective alternative to traditional weight selection methods. Meng et al. ³² developed a Simulated Annealing (SA)-based methodology for LQR weight matrix optimization, significantly enhancing calibration efficiency. Nevertheless, both APSO and SA algorithms exhibit susceptibility to local optima, leading to diminished optimization capabilities. In contrast, Genetic Algorithms (GA) have gained widespread adoption due to their inherent advantages in population diversity, crossover mechanisms, and parallel search capabilities, which confer superior robustness and adaptability in addressing global optimization, complex constraints, and multimodal problems.

This study introduces differential geometry theory to construct an improved LQR controller with high model accuracy and optimizes the weighting coefficient matrices of the controller using a Genetic Algorithm (GA). The research aims to overcome the control performance degradation caused by treating nonlinear characteristic parameters as constants in traditional LQR controllers when handling suspension models, as well as the problem in configuring weighting coefficient matrices. First, a 1/4-scale nonlinear active suspension model is established using differential geometry theory, explicitly incorporating the nonlinear characteristics of suspension parameters. Based on this model, an improved LQR controller considering these nonlinear features is developed. Subsequently, a performance evaluation function, that is, a fitness function, is designed to assess the performance of the LQR controller. The weighting coefficient matrices of the LQR controller are optimized via the Genetic Algorithm to determine the optimal feedback control strategy. Finally, MATLAB-Simulink simulations are conducted to replicate vehicle driving conditions under discrete impact road models and random road models. The effectiveness of the proposed control strategy is verified by comparing the vertical body acceleration, suspension dynamic deflection, and tire dynamic displacement between passive suspension systems and PID-controlled active suspension systems.

Construction of nonlinear active suspension model

According to Newton’s second law, a 1/4 nonlinear active suspension model is established. As shown in Figure 1, the 1/4 suspension model is a two-degree-of-freedom model, including the vertical displacements of the vehicle body and unsprung mass. The spring and damping of the suspension are located between the vehicle body and unsprung mass. The stiffness of the vehicle tire is equivalent to a linear spring. The equations of motion of the two-degree-of-freedom nonlinear 1/4 suspension model are shown in equation (1):

$$\{\begin{array}{c}{x_b}^{\prime \prime }{m}_b=-F_{\mathrm{ks}}-F_c-u\hfill \\ {x_w}^{\prime \prime }{m}_w=+F_{\mathrm{ks}}+F_c+u-K_t(x_w-x_r)\hfill \end{array}$$

(1)

Figure 1.

One-fourth nonlinear active suspension model.

In the equation, ${x_b}^{″}$ and ${x_w}^{″}$ denote the acceleration of the car body and the unsprung mass, respectively. $F_{\mathrm{ks}}$ and $F_c$ represent the suspension spring force and damping force, respectively. $x_w$ and $x_r$ refer to the displacement of the unsprung mass and the road surface input displacement, respectively. $K_t$ stands for the vertical dynamic stiffness of the tire. $m_b$ and $m_w$ indicate the car body and unsprung mass, respectively. $u$ represents the suspension control force.

In order to obtain the accurate changes of spring output force and damping output force, referring to a certain class B car in CarSim software, the relationship between spring output force and suspension stroke and the relationship between damping output force and suspension movement speed are obtained respectively. Through the MATLAB curve fitting toolbox, the fitting curve is obtained. The relationship between the compression stroke of the suspension and the spring force is shown in Figure 2. It can be seen that in the range of −0.005 to 0 m of the suspension stroke, the suspension stroke and the spring output force are approximately linearly related, and the rest show a nonlinear relationship. When the suspension stroke is zero, that is, in a stationary state, the spring output force is not zero. The reason is that the spring is affected by the gravity of the vehicle body. The spring has been deformed when the vehicle is stationary. The nonlinear spring fitting curve is in the form of a cubic equation, as shown in the following equation (2):

$$F_{\mathrm{ks}}=F_{k0}+k_0(x_b-x_w{)}^3+k_1(x_b-x_w{)}^2+k_2(x_b-x_w)$$

(2)

Figure 2.

Suspension spring force fitting curve.

In the equation, $F_{k0}$ represents the spring output force when the vehicle is at rest; $k_0$ , $k_1$ , and $k_2$ are the fitting coefficients of the third-order polynomial, with values of 95,350,000, 4,038,000, and 251,500, respectively. In the calculation process of the two-degree-of-freedom nonlinear suspension model, the compression of the spring by the car body when the vehicle is at rest is not considered. The fitting curve is written as shown in equation (3):

$$F_{\mathrm{ks}}=k_0(x_b-x_w{)}^3+k_1(x_b-x_w{)}^2+k_2(x_b-x_w)$$

(3)

Similarly, the relationship between the expansion and contraction rate of the damper and the damping force is shown in Figure 3. The nonlinear damping fitting curve can be written in the form of a quadratic polynomial, as shown in equation (4):

$$F_c=b_0({x_b}^{\prime }-{x_w}^{\prime }{)}^2+b_1({x_b}^{\prime }-{x_w}^{\prime })$$

(4)

Figure 3.

Suspension damping force fitting curve.

In the equation, $b_0$ and $b_1$ are the coefficients of the second-order polynomial, valued at −322.6 and 3193, respectively.

Control methods

Linearization of nonlinear 1/4 active suspension model

Theoretical basis of differential geometry

For nonlinear problems, the differential geometry method can retain the high-order terms of nonlinearity, and realize the linear decoupling of nonlinear systems through coordinate transformation and nonlinear retention methods. Write the state equation and output equation of the nonlinear suspension as:

$$\{\begin{array}{c}{x}^{\prime }=f(x)+g(x)u+\eta \omega \hfill \\ Y=h(x)+d(x)u\hfill \end{array}$$

(5)

In the equation, $x\in R^n$ denotes the system state vector; $y\in R^1$ is the system output; $u\in R^1$ is the controllable system input; $f$ and $g$ are the nonlinear function vectors.

There is a positive integer r such that the system meets the condition in the vicinity of a specific operating point $X=X_0$ :

$$\{\begin{array}{c}{L}_g{L}_f^{i}h(x)\equiv 0\mathrm{i}=1,2···,\mathrm{r}-2\hfill \\ L_g{L}_f^{r-1}h(x)\ne 0\hfill \end{array}$$

(6)

In this context, the relative order of the system in the neighborhood of $X_0$ is r. In equation (6): The operation of $L_g{L}_f^{i}h(x)$ is called Lie derivative operation, and the operation rule is as shown in equation (7).

$$\{\begin{array}{c}{L}_f^0\lambda (x)=\lambda (x)\hfill \\ L_f\lambda (x)=\sum _{i=1}^{i=n}\frac{\partial \lambda (x)}{\partial x_i}{f}_i(x)\hfill \\ L_f^k\lambda (x)=L_f(L_f^{k-1}\lambda (x))\hfill \\ L_g{L}_f^k\lambda (x)=\frac{\partial (L_f^k(\lambda (x)))}{\partial x}g(x)\hfill \end{array}$$

(7)

When the relative order of the nonlinear system expression is r and satisfies $\mathrm{\Omega }(x)·\eta =0$ , the linear feedback expression can be solved:

$$u=\frac{1}{L_g{L}_f^{r-1}h(x)}[-L_f^{r}h(x)+v]$$

(8)

Where:

$$\mathrm{\Omega }(x)=\left[\begin{array}{c}\frac{\partial h(x)}{\partial (x)}\hfill \\ \multicolumn{1}{c}{\frac{\partial (L_{f}h(x))}{\partial (x)}}\\ \multicolumn{1}{c}{···}\\ \multicolumn{1}{c}{\frac{\partial (L_f^{r-1}h(x))}{\partial (x)}}\end{array}\right]$$

(9)

Differential geometric linear suspension model

Select suspension dynamic deflection $(x_b-x_w)$ , vehicle body vertical velocity ${x_b}^{\prime }$ , wheel dynamic displacement $(x_w-x_r)$ and unsprung mass vertical velocity ${x_w}^{\prime }$ as state variables; select vehicle body vertical acceleration ${x_b}^{″}$ , suspension dynamic deflection $(x_b-x_w)$ and wheel dynamic displacement $(x_w-x_r)$ as output observation values. The state quantity matrix expression is:

$$x={\left[\begin{array}{cccc}{x}_b-x_w\hfill & {x_b}^{\prime }\hfill & x_w-x_r\hfill & {x_w}^{\prime }\hfill \end{array}\right]}^T={\left[\begin{array}{cccc}{x}_1\hfill & x_2\hfill & x_3\hfill & x_4\hfill \end{array}\right]}^T$$

The matrix expression for the output observation values is:

$$Y={\left[\begin{array}{ccc}{x_b}^{\prime \prime }\hfill & x_b-x_w\hfill & x_w-x_r\hfill \end{array}\right]}^T={\left[\begin{array}{ccc}{y}_1\hfill & y_2\hfill & y_3\hfill \end{array}\right]}^T$$

The state space equation and output equation for the nonlinear active suspension are defined as:

$$\{\begin{array}{c}{x}^{\prime }=f(x)+g(x)u+\eta \omega \hfill \\ Y=h(x)+d(x)u\hfill \end{array}$$

(10)

Where:

$$\begin{array}{c}f(x)={\left[\begin{array}{cccc}{x}_2-x_4\hfill & f_1{(x)}^{\prime }\hfill & x_4\hfill & f_2{(x)}^{\prime }\hfill \end{array}\right]}^T\hfill \\ g(x)={\left[\begin{array}{cccc}0\hfill & \frac{1}{m_b}\hfill & 0\hfill & -\frac{1}{m_w}\hfill \end{array}\right]}^T\hfill \\ \eta ={\left[\begin{array}{cccc}0\hfill & 0\hfill & -1\hfill & 0\hfill \end{array}\right]}^T\hfill \\ h(x)={\left[\begin{array}{ccc}{f}_1{(x)}^{\prime }\hfill & x_1\hfill & x_3\hfill \end{array}\right]}^T\hfill \\ d(x)={\left[\begin{array}{ccc}\frac{1}{m_b}\hfill & 0\hfill & 0\hfill \end{array}\right]}^T\hfill \end{array}$$

$$\begin{array}{c}{f}_1(x{)}^{\prime }=[-95350000\times (x_b-x_w{)}^3-4038000\times (x_b-x_w{)}^2\hfill \\ -251500\times (x_b-x_w)+322.6\times ({x_b}^{\prime }-{x_w}^{\prime }{)}^2\hfill \\ -3193\times ({x_b}^{\prime }-{x_w}^{\prime })]/m_b\hfill \end{array}$$

$$\begin{array}{c}{f}_2(x{)}^{\prime }=[95350000\times (x_b-x_w{)}^3+4038000\times (x_b-x_w{)}^2\hfill \\ +251500\times (x_b-x_w)-322.6\times ({x_b}^{\prime }-{x_w}^{\prime }{)}^2\hfill \\ +3193\times ({x_b}^{\prime }-{x_w}^{\prime })]/m_w\hfill \end{array}$$

Let $h_1(x)=x_1$ , when $k=0$ , have:

$$\begin{array}{c}{L}_g{L}_f^0{h}_1(x)=\frac{\partial L_f^0{h}_1(x)}{\partial x^T}g(x)\hfill \\ =\frac{\partial h_1(x)}{\partial x^T}g(x)=\left[\begin{array}{cccc}1\hfill & 0\hfill & 0\hfill & 0\hfill \end{array}\right]·\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\frac{1}{m_b}}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-\frac{1}{m_w}}\end{array}\right]=0\hfill \end{array}$$

When $k=1$ , have

$$\begin{array}{c}{L}_g{L}_f^0{h}_1(x)=\frac{\partial L_f^1{h}_1(x)}{\partial x^T}g(x)\hfill \\ =\frac{\partial [\frac{\partial h_1(x)}{\partial x^T}f(x)]}{\partial x^T}g(x)=\left[\begin{array}{cccc}0\hfill & 1\hfill & 0\hfill & -1\hfill \end{array}\right]·\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\frac{1}{m_b}}\\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-\frac{1}{m_w}}\end{array}\right]\hfill \\ =\frac{1}{m_b}+\frac{1}{m_w}\ne 0\hfill \end{array}$$

It can be concluded that the relative order $r=2$

$$\begin{array}{c}\mathrm{\Omega }(x)=\left[\begin{array}{cccc}\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{-1}\end{array}\hfill \end{array}\right]\hfill \\ \mathrm{\Omega }(x)·\eta =\left[\begin{array}{cccc}\begin{array}{c}1\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{-1}\end{array}\hfill \end{array}\right]·\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{-1}\\ \multicolumn{1}{c}{0}\end{array}\right]=\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]\hfill \end{array}$$

We can derive:

$$\begin{array}{l}{u}_1=95350000\times {(x_b-x_w)}^3+4038000\times {(x_b-x_w)}^2\\ +251500\times (x_b-x_w)-322.6\times {({x^{\prime }}_b-{x^{\prime }}_w)}^2\\ +3193\times ({x^{\prime }}_b-{x^{\prime }}_w)-\frac{m_b{k}_t}{m_b+m_w}{x}_3+\frac{m_b{m}_w}{m_b+m_w}{v}_1\end{array}$$

Similarly let $h_2(x)=x_2$ , the relative order $r=1$ .

$$\begin{array}{c}{u}_2=-[95350000\times (x_b-x_w{)}^3+4038000\times (x_b-x_w{)}^2\hfill \\ +251500\times (x_b-x_w)-322.6\times ({x_b}^{\prime }-{x_w}^{\prime }{)}^2\hfill \\ +3193\times ({x_b}^{\prime }-{x_w}^{\prime })]+m_b{v}_2\hfill \end{array}$$

Let $h_3(x)=x_3$ , the relative order $r=2$ .

$$\begin{array}{c}{u}_3=-[95350000\times (x_b-x_w{)}^3+4038000\times (x_b-x_w{)}^2\hfill \\ +251500\times (x_b-x_w)-322.6\times ({x_b}^{\prime }-{x_w}^{\prime }{)}^2\hfill \\ +3193\times ({x_b}^{\prime }-{x_w}^{\prime })]-k_t{x}_3-m_w{v}_3\hfill \end{array}$$

Let $h_4(x)=x_4$ , the relative order $r=1$ .

$$\begin{array}{c}{u}_4=95350000\times (x_b-x_w{)}^3+4038000\times (x_b-x_w{)}^2\hfill \\ +251500\times (x_b-x_w)-322.6\times ({x_b}^{\prime }-{x_w}^{\prime }{)}^2\hfill \\ +3193\times ({x_b}^{\prime }-{x_w}^{\prime })-k_t{x}_3-m_w{v}_4\hfill \end{array}$$

Contrasting with the previous equation, the control rate expression for the nonlinear active suspension model is chosen as:

$$\begin{array}{c}u=95350000\times (x_b-x_w{)}^3+4038000\times (x_b-x_w{)}^2\hfill \\ +251500\times (x_b-x_w)-322.6\times ({x_b}^{\prime }-{x_w}^{\prime }{)}^2\hfill \\ +3193\times ({x_b}^{\prime }-{x_w}^{\prime })-v\hfill \end{array}$$

(11)

Substituting equation (11) into (5), we obtain the linear mathematical model of the nonlinear active suspension:

$$\{\begin{array}{c}{x}^{\prime }=\mathrm{Ax}+\mathrm{Bv}+\Gamma \omega \hfill \\ Y=\mathrm{Cx}+\mathrm{Dv}\hfill \end{array}$$

(12)

Where:

$$\begin{array}{c}A=\left[\begin{array}{cccc}0\hfill & 1\hfill & 0\hfill & -1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 1\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -\frac{k_t}{m_w}\hfill & 0\hfill \end{array}\right]\hfill \\ B={\left[\begin{array}{cccc}0\hfill & -\frac{1}{m_b}\hfill & 0\hfill & \frac{1}{m_w}\hfill \end{array}\right]}^T\hfill \\ \Gamma ={\left[\begin{array}{cccc}0\hfill & 0\hfill & -1\hfill & 0\hfill \end{array}\right]}^T\hfill \\ C=\left[\begin{array}{cccc}\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{1}\end{array}\hfill & \begin{array}{c}0\hfill \\ \multicolumn{1}{c}{0}\\ \multicolumn{1}{c}{0}\end{array}\hfill \end{array}\right]\hfill \\ D=\left[\begin{array}{ccc}-\frac{1}{m_b}\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \end{array}$$

Design of LQR controller optimized by genetic algorithm

In order to improve the ride comfort and operational stability of the vehicle, this paper aims to reduce the vertical acceleration of the vehicle body, suspension dynamic deflection and tire dynamic displacement. Therefore, the performance index of the LQR controller designer is:

$$\begin{array}{c}J={\int }_0^{\infty }[q_1{z_b}^{\prime \prime 2}+q_2(x_b-x_w{)}^2+q_3(x_w-x_r{)}^2]\mathrm{dt}\hfill \\ ={\int }_0^{\infty }[Y^T\mathrm{QY}+u^T\mathrm{Ru}]\mathrm{dt}\hfill \end{array}$$

(13)

Let the weight coefficient matrix be:

$$\begin{array}{c}Q=\left[\begin{array}{ccc}{q}_1\hfill & \hfill & \hfill \\ \hfill & q_2\hfill & \hfill \\ \hfill & \hfill & q_3\hfill \end{array}\right]\hfill \\ R=\left[R\right]\hfill \end{array}$$

Equation (13) is further arranged as:

$$\begin{array}{c}J={\int }_0^{\infty }[(\mathrm{Cx}+\mathrm{Dv}{)}^{T}Q(\mathrm{Cx}+\mathrm{Dv})+v^T\mathrm{Rv}]\mathrm{dt}\hfill \\ ={\int }_0^{\infty }[x^T{Q}_{d}x+2x^T{N}_{d}v+v^T\mathrm{Rv}]\mathrm{dt}\hfill \end{array}$$

(14)

Where:

$$\begin{array}{c}{Q}_d=C^T\mathrm{QC}\hfill \\ N_b=C^T\mathrm{QD}\hfill \\ R_d=R+D^T{\mathrm{QD}}_{°}\hfill \end{array}$$

Enter the equation (14):

$$\begin{array}{c}{Q}_d=\left[\begin{array}{cccc}{q}_2\hfill & 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& q_3\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & 0\hfill & 0\hfill \end{array}\right]\hfill \\ N_d={\left[\begin{array}{cccc}0& 0& 0& 0\end{array}\right]}^T\hfill \\ R=\frac{q_1}{m_b^2}\hfill \end{array}$$

The optimal control force of the actuator output is:

$$u=-\mathrm{kx}(t)$$

(15)

In the equation, $k$ is the optimal control feedback matrix; $x(t)$ is the state feedback matrix at any time.

The optimal feedback gain matrix $k$ can be determined using MATLAB’s linear quadratic optimal controller, as shown in equation (16).

$$[k,s,e]=\mathrm{LQR}(A,B,Q,R,N)$$

(16)

As can be seen from equations (13) to (16), the optimal control force is greatly related to the weight coefficient. However, the weight coefficient is often selected based on experience. Although the subjective adjustment of the weight coefficient can obtain the corresponding optimal control law, there is great subjectivity. To solve this problem, this paper uses a genetic algorithm to optimize the weight coefficients of the optimal linear controller (LQR) controller. Considering that the three performance indicators of the suspension system are both related and contradictory, the weighted summation method used in this paper can turn multi-objective optimization into single-objective optimization, simplifying the complexity of the problem. At the same time, it is necessary to ensure that the suspension stroke cannot exceed the limit stroke of the suspension. Therefore, the fitness function designed in this study is shown in equation (17). The optimal weight coefficient is determined by finding the minimum value of the fitness function.

$$\mathrm{minimize}\mathrm{L}=\frac{\mathrm{BA}(x)}{B{A}_{\mathrm{pas}}}+\frac{\mathrm{SWS}(x)}{\mathrm{SW}{S}_{\mathrm{pas}}}+\frac{\mathrm{DTP}(x)}{\mathrm{DT}{P}_{\mathrm{pas}}}$$

(17)

In the equation, $\mathrm{BA}$ , $\mathrm{SWS}$ , $\mathrm{DTP}$ represent the body vertical acceleration, suspension dynamic travel, and tire dynamic displacement under the corresponding control laws. $B{A}_{\mathrm{pas}}$ , $\mathrm{SW}{S}_{\mathrm{pas}}$ , $\mathrm{DT}{P}_{\mathrm{pas}}$ are the performance indicators for the passive suspension.

The operation steps of genetic algorithm are as follows: ① Generate a population and assign it to the weight coefficient $q_1$ , $q_2$ , $q_3$ . The optimal feedback gain matrix is obtained through equation (16). The optimal control force can be obtained from equation (15) and brought into the nonlinear suspension model to obtain the suspension performance index. ② Calculate the fitness value of each group in various populations by equation (17) and judge whether the termination condition of genetic algorithm is satisfied. If satisfied, exit the genetic algorithm and obtain the optimal individual. If not satisfied, jump to step ③. ③ The genetic algorithm performs selection, crossover, mutation and generates a new population, and jumps to step ②. The flow chart of optimizing weighting coefficients of LQR controller is shown in Figure 4.

Figure 4.

LQR controller optimization weighting coefficient flow chart.

Taking the discrete impact road surface as an example, in the optimization process of genetic algorithm, the change of fitness function value of the optimal individual is shown in Figure 5. As the population is continuously optimized, the fitness value of the optimal individual continuously decreases and finally converges to 0.00495. At this time, the corresponding optimal individual $(q_1,q_2,q_3)=(206041,239891,111158)$ .

Figure 5.

Optimal individual and fitness function change.

The optimal feedback matrix determined through the genetic algorithm is:

$$k=[-135916,-7846.71,-95124.6,-3205.98]$$

The optimal control force is then determined as:

$$u=135916\times x_1+7846.71\times x_2+95124.6\times x_3+3205.98\times x_4$$

By the same token, the optimal controlling force of Class C random pavement is:

$$\begin{array}{c}u=262720.21\times x_1+9394.54\times x_2+92914.48\times x_3\hfill \\ -3076.31\times x_4\hfill \end{array}$$

Simulation results and analysis

The relevant parameters of a certain B-class car are: vehicle body mass $m_b=277.5\mathrm{kg}$ , unsprung mass $m_w=54.5\mathrm{kg}$ , and tire vertical dynamic stiffness $K_t=232000N/m$ . The road excitation adopts the discrete impact road surface model and the C-class road model. The formula of the discrete impact road surface model is as shown in equation (18). The schematic diagram of the discrete impact road surface model is shown in Figure 6.

$$\{\begin{array}{c}y=4*t-102.5\le \mathrm{t}<2.51\hfill \\ y=0.042.51\le \mathrm{t}<2.53\hfill \\ y=-4*t+10.162.53\le \mathrm{t}<2.54\hfill \\ y=0\mathrm{else}\hfill \end{array}$$

(18)

Figure 6.

Diagram of discrete impact pavement model.

The time-domain model of Class C road roughness is as shown in equation (19). The schematic diagram of the Class C road model is shown in Figure 7.

$$z_g(t{)}^{\prime }=-2\pi n_{1}u{z}_g(t)+2\pi n_0\sqrt{G_q(n_0)u}\omega (t)$$

(19)

Figure 7.

Class C road height information.

In the equation, $z_g$ is the road height; $n_1$ is the lower cutoff spatial frequency, taking a value of $0.01m^{-1}$ ; $u$ is the vehicle speed, taking a value of $20m/s$ ; $n_0$ is the reference spatial frequency, taking a value of $0.1m^{-1}$ ; $\sqrt{G_q(n_0)}$ is the road roughness coefficient, with a value of $256\times {10}^{-6}{m}^3$ for a C-class road; $\omega (t)$ is road white noise, with noise power set to 0.5 and sampling time set to 0.005.

According to the 1/4 suspension dynamics formula, a simulation model is built in Simulink software, and the step size is selected as 0.005. After converting the nonlinear suspension model into a linear model by using the differential geometry method, the GA-LQR control algorithm is used to obtain the optimal feedback control law, thereby obtaining the feedback control force of the nonlinear suspension system. The control flow chart of the nonlinear suspension system is shown in Figure 8. The working conditions are discrete impact road surface and Class C random road surface respectively. The observed quantities are body vertical acceleration, suspension dynamic deflection and tire dynamic displacement. The GA-LQR control algorithm is compared with passive suspension and traditional PID algorithm. The Simulink active suspension simulation diagram is shown in Figure 9.

Figure 8.

Nonlinear suspension control flow chart.

Figure 9.

Simulink active suspension simulation.

As can be seen from Figures 10 to 12, under the condition of discrete impact road surface, compared with passive suspension and active suspension using traditional PID control algorithm, the GA-LQR control algorithm can not only significantly reduce the vertical acceleration of the vehicle body, suspension dynamic deflection and tire dynamic displacement, but also reduce the time for each observed quantity to reach a stable state.

Figure 10.

Body vertical acceleration response under discrete impact.

Figure 11.

Suspension dynamic disturbance response under discrete impact.

Figure 12.

Tire dynamic displacement response under discrete impact.

It can be concluded from Table 1 that compared with passive suspension, the root mean square values of vehicle body vertical acceleration, suspension dynamic disturbance and tire dynamic displacement by using GA-LQR control algorithm are reduced by 17.3%, 45.3% and 10.4% respectively. Compared with the traditional PID control algorithm, the root mean square values of vehicle body vertical acceleration, suspension dynamic disturbance and tire dynamic displacement by using GA-LQR control algorithm are reduced by 6.94%, 57.72% and 4.4% respectively. The traditional PID control algorithm cannot reduce the suspension dynamic disturbance, and compared with passive suspension, it is increased by 29.4%.

Table 1.

Comparison of performance indexes in discrete impact conditions.

Performance index	Passive suspension	PID	GA-LQR
Body vertical acceleration $(m/s^2)$	3.57	3.17	2.95
Suspension dynamic disturbance $(m)$	0.0034	0.0044	0.00186
Tire dynamic displacement $(m)$	0.0048	0.0045	0.0043

As can be seen from Figures 13 to 15, under the condition of Class C random road surface, compared with passive suspension and active suspension using traditional PID control algorithm, the GA-LQR control algorithm can reduce the vertical acceleration of the vehicle body, suspension dynamic deflection and tire dynamic displacement, among which the suspension dynamic deflection can be significantly reduced.

Figure 13.

Body vertical acceleration response under Class C road surface.

Figure 14.

Suspension dynamic deflection response under Class C road surface.

Figure 15.

Tire dynamic displacement response under Class C road surface.

It can be concluded from Table 2 that compared with passive suspension, the root mean square values of vehicle body vertical acceleration, suspension dynamic disturbance and tire dynamic displacement by using GA-LQR control algorithm are reduced by 28.2%, 63.6% and 23% respectively. Compared with the traditional PID control algorithm, the root mean square values of vehicle body vertical acceleration, suspension dynamic disturbance and tire dynamic displacement by using GA-LQR control algorithm are reduced by 22.1%, 63.6% and 16.25% respectively. The traditional PID control algorithm cannot reduce suspension dynamic disturbance and tire dynamic displacement.

Table 2.

Comparison of performance indexes under Class C pavement conditions.

Performance index	Passive suspension	PID	GA-LQR
Body vertical acceleration $(m/s^2)$	6.38	5.88	4.58
Suspension dynamic disturbance $(m)$	0.0066	0.0066	0.0024
Tire dynamic displacement $(m)$	0.0087	0.008	0.0067

Conclusion

This paper conducts research on the control strategy of nonlinear active suspension systems in view of the significant nonlinear characteristics of automobile suspension systems. First, the nonlinear active suspension model is transformed into a linear model by using differential geometry theory. Then, the linear quadratic regulator (LQR) is optimized by genetic algorithm (GA) to obtain the feedback control force of the active suspension, and comparative studies are carried out under Class C random road surface and discrete impact road surface conditions respectively. The simulation results show that:

• (1) The fitness function designed in this study can effectively take into account each evaluation index of the suspension system. At the same time, combined with genetic algorithm, it can realize the rapid optimization of the weight coefficient of the linear quadratic regulator controller. The setting relationship between the active suspension performance and the weight coefficient matrix in the LQR controller is clarified, avoiding the empiricism and subjectivity in the traditional weight coefficient setting process and significantly improving the coefficient setting efficiency.

• (2) The genetic algorithm optimized linear quadratic regulator (GA-LQR) control algorithm proposed in this study can effectively improve the comprehensive vibration reduction performance of the suspension system. Compared with the traditional proportional-integral-derivative (PID) control, under the discrete impact road surface, the GA-LQR algorithm reduces the root mean square values of vehicle body vertical acceleration, suspension dynamic deflection and tire dynamic displacement by 6.94%, 57.72% and 4.4% respectively; under Class C road surface, the above three indicators are reduced by 22.1%, 63.6% and 16.25% respectively.