Research on nonlinear model and fuzzy fractional order PI λ D μ control of air suspension system

Abstract

To improve the ride comfort of wheeled armored vehicles, air springs are used. To describe the vehicle motion more accurately, a nine-degree-of-freedom air suspension system for the whole vehicle was established, and its equations of motion were derived. Through theoretical analysis of the stiffness characteristics and forces on the air springs, the nonlinear restoring force was obtained as a cubic polynomial of the air spring displacement. The simulation results obtained by fitting the polynomial and radial basis function curves with MATLAB/Simulink software are consistent with the actual test results, thus verifying the correctness of the nonlinear air spring polynomial model. Finally, a fuzzy fractional order PI^λD^μ controller is designed and simulated for the vehicle-seat-body model in terms of wheel dynamic load, suspension dynamic deflection, body acceleration, and other indicators. The simulation results show that the fuzzy fractional order PI^λD^μ Proportion Integral Differential (PID) control strategy has better overall performance than the PID control strategy, fuzzy control strategy, and fuzzy PID control strategy.

1. Introduction

1.1. Current status of air suspension research

The control strategy of air suspension is one of the key technologies of the vehicle suspension system, which is also related to the ride comfort and handling stability of vehicles with air springs. In the 1960s, scholars have proposed the concept of active suspension. With the development of industry in the 1970s, active ^1,2 and semi-active ^3,4 suspension systems have been developed. Scholars have applied optimal control ,^5–8 adaptive control ,^9–11 fuzzy control ,^12,13 and artificial neural networks ^14–16 to vibration control of vehicle suspension systems. Kim proposed asynchronous and synchronous methods to control the height of air suspension . ¹⁷ Fuzzy control was used to solve the uncertainty problem of time-varying mode, and the correctness of the two methods was verified by simulation. Gao Zepeng et al. ¹⁸ used fuzzy control to model the electronically controlled air suspension and solved the overshoot phenomenon in the suspension static regulation under unbalanced load . New system control and fault handling method are proposed for the four-point closed-loop air suspension system . ¹⁹ The sliding mode control with phase compensation feedback signal is used to improve the control accuracy and reduce the vibration of the car body. The fault handling method provides fault detection, diagnosis, and management for the closed-loop air suspension. Then Kim used sliding mode control to adjust the height of the suspension and improve the roll and pitch of the vehicle body. ²⁰ Xiao-Qiang Sun et al. ²¹ proposed a fuzzy-sliding control for the electronically controlled air suspension to improve the accuracy of body height adjustment and attitude leveling, for that the car acquired better control characteristics. Prabu established a nonlinear mathematical model of two-degree-of-freedom air suspension system. ²² Based on the neural-slide control theory, Niao-na Zhang et al. ²³ control the vibration of a two-degree-of-freedom semi-active air suspension. The simulation results show that this control method can reduce chattering. W.N. Bao uses fuzzy and adaptive sliding control to simulate the indexes of suspension ride comfort. Compared with the traditional sliding control, the FASM control effect is more obvious. ²⁴ By comparing three tuning methods of Ziegler Nicholas, Refined Ziegler Nicholas, and Optimal control, the controller is designed by setting the appropriate gain value. Three gain values of PID controller are obtained by three different adjusting methods to control the air suspension system. The fractional order sliding mode control for semi-active suspension is proposed by Sy Dzuyen Nguyen. The system obtained high performance control under uncertainty and external disturbance. ²⁵ Sathishkumar used a fuzzy control strategy to actively control the 3-DOF air suspension system model of 1/4 vehicle. ²⁶ In order to analyze the performance of the active suspension system, the seat and body acceleration were selected to judge the state of the control system. At the same time, the semi-active air suspension system of light vehicle was analyzed by fuzzy and PID control strategies. ²⁷ Eskandary developed a mathematical model of an air suspension system with independent height and stiffness regulation. ²⁸ By controlling the air pressure of the air chamber, the stiffness and height of the vehicle can be changed simultaneously under different driving conditions to keep the performance of the vehicle in optimal condition. Xubin Dong et al. ²⁹ proposed a fractional air disaster eradication method based on the adaptive clonal algorithm for active suspension, which reduces the requirement of model accuracy. The optimum mode and self-adaptive PID control for active suspension are presented by Nonlinear Dyn et al. The optimum mode and self-adaptive PID control for active suspension are presented by M. J. Mahmoodabadi et al. ³⁰ Based on the design of self-adaptive rules, the parameters of the system are adjusted by means of the mode and system, so that the system has a small overshoot and greatly reduces the vibration.

In recent years, more and more researches have been done on active suspension, and the combined control effect has been optimized. Many control system strategies have been discussed by researchers to enhance the performance of air suspension systems. This paper focuses on the combination of fitting response curves for air springs and comparison of performance of different controllers for active suspension. Based on the theoretical analysis of the stiffness characteristics and force of air spring, the air spring model was selected according to the vehicle type, and then the air spring characteristic polynomial fitting was carried out to establish the air spring polynomial model. The fuzzy fractional order PI^λD^μ control theory was proposed, and the controller was designed by analysis. The positive results of the fuzzy fractional order PI^λD^μ control were verified by Simulink simulation of the dynamic load, dynamic deflection, and vertical acceleration of the suspension.

1.2. Mathematical methods and the application in other fields

Most physical problems in life can be solved by differential equations, and there are many methods to obtain exact solutions of equations: differential constraint method, similar transformation method, etc. Lie symmetry analysis is the easiest to implement and the most effective method to obtain exact solutions, which was initially used only for local symmetry analysis, ³¹ and later Bluman et al. ³² expanded the symmetry category to appear non-local symmetry. Subhankar Sil et al. ³³ obtained the exact solution of the equation as well as the exact solution of the productive soliton equation by non-local symmetry. Adomian decomposition method (ADM) is an effective method for linear and nonlinear partial differential equations, which can provide analytic approximation solutions for differential equations, but there are shortcomings in solving the Riemann problem. Dia Zeidan et al. improved the approximation method for hyperbolic partial differential equations in the selection of characteristic curves to compensate for the shortcomings of the small radius of convergence of AMD. And in addition, the smooth solution of hyperbolic partial differential will be in finite time, breaking within a finite time gives rise to the concept of discontinuous waves. Dia Zeidan et al.^34–35 studied the evolution of the plane and cylindrical symmetric motion of the gas dynamics equations and discontinuity surfaces using the van der Waals model for magnetic gases. Rahul Kumar Chaturvedi et al. ³⁶ also studied the Riemann and non-classical Riemann solutions of the non-simultaneous partial differential equations for the fluid MCG model, obtaining that the fundamental in the solution waveforms can change due to external forces. For a two-phase fluid (at least one phase is a fluid) as in a suspension model with kinematic coupling, where the flow of one phase always involves the relative motion of the other phase, the numerical solution and characteristics of the two-phase flow equations and the Riemann numerical solution of the compressible two-phase flow were studied by D. Zeidan. ³⁷ The Riemann analytic solution of the drift flux problem was determined in a paper by S. Kuila ³⁸ in a constant form and described its complete solution strategy. The Lie symmetry problem in the differential equations for polyatomic gases was analyzed and the evolutionary behavior of the excitation structures, weak discontinuity surfaces, characteristic excitations, and their interactions in gases under dynamic pressure was studied by Dia Zeidan et al. ³⁹

S. Ahmad et al. ⁴⁰ proposed the use of fuzzy fractional order differential equations in the Caputo sense to describe novel coronavirus infection systems, leading to a better understanding of the dynamical structure of their behavior. Ullah A ⁴¹ used fuzzy Caputo fractional derivatives to study complex population dynamics and fuzzy Laplace transform and Adomian decomposition while to the general numerical solution of the model. Ahmad N et al. ⁴² proved the existence and uniqueness of the fuzzy fractional Volterra–Fredholm integer difference equation .Ahmad S et al. ⁴³ derived a semi-analytic solution of the time-fractional Fisher equation under fuzziness and analyzed the convergence and error of the Laplace–Adomian decomposition method (LADM). Later, Ahmad S ⁴⁴ studied the fuzzy dispersion PDE under Caputo, Caputo–Fabrizio, and Atangana–Baleanu fractional operators and verified the values of different orders and uncertainties.

2. Stiffness and force characteristics of air spring

The load acting on the air spring depends on the compressed gas in the spring. According to the gas equation of state, the volume and pressure satisfy the following equation in the working process of the air spring ⁴⁵

$$\left(p-p_a\right)V^m=const$$

(1)

in which p—Absolute pressure of gas in a gasbag at static equilibrium position, $p_a$ —atmospheric pressure, V—Air spring volume, m—Variable index, and const—Constant. Variable index $m$ is determined by the deformation speed of air spring. When the spring vibration speed is small or static characteristics are studied, the gas state in the spring is approximately isothermal, and $m$ takes 1.0. When the vibration speed of the spring is large or the dynamic characteristics are studied, the gas state in the spring is approximately adiabatic, and $m$ takes 1.38. Considering the atmospheric pressure, the air spring is filled with gas and has loads on the air springs,^45,46 then

$$F=\left(p-p_{\text{a}}\right)A$$

(2)

in which $A$ is the effective area of the air spring at the specified working status and its value varies with the height of the air spring.^46,47

$$A=A\left(h\left(t\right)\right)$$

(3)

in which h is the height of air spring.

The stiffness of air springs can be obtained by derivation of air spring displacement. $x$ ¹⁸

$$K=\frac{dF}{dx}=A\hspace{0.17em}\frac{dp}{dx}+\left(p-p_{\text{a}}\right)\hspace{0.17em}\frac{dA}{dx}$$

(4)

According to the analysis of the thermodynamic process inside the air spring and the equation (1) as well as (4), the vertical stiffness of the air spring can be given by the following formula.

$$K=m{A}^2\frac{p-p_{\text{a}}}{V}+\left(p-p_{\text{a}}\right)\hspace{0.17em}\frac{dA}{dx}$$

(5)

Among them, the payload area satisfies ⁴⁸

$$A=-\frac{dV}{dx}$$

(6)

At the same time, the effective area and the displacement are linear in the change rate of the payload area of the air spring ⁴⁸

$$A=Cx+A_0$$

(7)

Substituting equation (7) into equation (4) yields

$$K=m\left(p-p_{\text{a}}\right)\frac{{\left(Cx+A_0\right)}^2}{V}+C\left(p-p_{\text{a}}\right)=\frac{m{C}^2\left(p-p_{\text{a}}\right)}{V}{x}^2+\frac{2mC{A}_0\left(p-p_{\text{a}}\right)}{V}x+\frac{m{A}_0^2\left(p-p_{\text{a}}\right)}{V}+C\left(p-p_{\text{a}}\right)$$

(8)

It can be concluded from equation (8) that the air spring suspension stiffness is a quadratic polynomial of displacement, so it can be derived that the load force is a cubic polynomial of displacement. This conclusion can provide a theoretical basis for the polynomial fitting of air spring force.

3. Mathematical model of air spring force

3.1. Polynomial fitting of air spring force

In the paper, the parameters of wheeled armored vehicle suspension system in reference are selected as shown in Table 1, ⁴⁹ in which the no-load mass is 7400 Kg. According to the manual of air spring, the Qianjin brand 1884N rubber air spring is selected. The specific data are shown in Tables 2 and 3. The performance curve of the air spring is shown in Figure 1. ⁵⁰

Table 1.

Suspension system parameters of “vehicle-seat-body” model of vehicle.

Parameters	Name	Numerical value	Unit
m h	Human body mass	65	kg
m c	Seat mass	30	kg
m s	Sprung mass (no load and full load)	7400/9400	kg
m u1, m u2 , m u3, m u4	Unsprung mass of suspension	400	kg
k h	Seat cushion equivalent stiffness	8228.7	N/m
k c	Seat spring stiffness	24×103	N/m
k u1, k u2, k u3, k u4	Tire equivalent stiffness of tire	539×103	N/m
c h	Seat cushion equivalent damping	152.8	Ns/m
c c	Seat equivalent damping	875	Ns/m
c u1, c u2, c u3, c u4	Damping of front and rear tire	5000	Ns/m
A	Distance between front axle center and center of mass	1.45	m
B	Distance between rear axle center and center of mass	1.75	m
W	1/2 wheelbase	0.9	m
D	Vertical distance from seat to y axis	0.5	m
E	Vertical distance from seat to x axis	0.4	m

Table 2.

Basic data of 1884N rubber air spring.

Name	Numerical value
Load capacity/Kg	1080–3046
Working stroke/mm	240
Recommended assembly design height/mm	280
Assembly stretch maximum height/mm	430
Assembly compressed minimum height/mm	190

Table 3.

1884N rubber air spring load capacity data sheet.

Assembly height/mm	Load capacity/Kg
	0.3 MPa	0.5 MPa	0.8 MPa
430	580	980	1550
350	900	1560	2600
280	1080	1850	3046
240	1450	2400	3700
190	2250	3300	5200

Figure 1.

Performance curve of 1884N rubber air spring.

According to the load capacity data table, the cubic polynomial fitting is carried out by using MATLAB software to get the curve of load capacity and displacement of air spring under different inflating pressure, as shown in Figure 2. In general, the height of the assembly is zero. The displacement is positive when compression, and the displacement is negative when stretching.

Figure 2.

Air spring load force-displacement curve.

It can be seen that the load force of the air spring is always positive and the direction is upward under different inflation pressures. However, compared with the load capacity of static balance, the size of air spring is positive and negative, that the load force is positive and upward when compressed below the static balance position; the load force is negative when stretched above the static balance position and downward. The relative spring negative load is a restoring force for the vibration of the car body, that is, $-\mathrm{\Delta }F=F-F_0$ . $\mathrm{\Delta }F$ —Nonlinear restoring force of air spring, $F$ —Air spring load force, and $F_0$ — Load acting on static equilibrium. $F_0=mg$ , $m$ —sprung mass. It should be addressed that the validity of these formulas requires suspension ratio (ratio between spring velocity and wheel vertical velocity) equal to 1.0.

The cubic polynomial function of the air spring nonlinear restoring force according to the polynomial fitting is

$$\mathrm{\Delta }F=a_{1}x+a_2{x}^2+a_3{x}^3$$

(9)

When the air pressure is 0.3 MPa, the polynomial fitting coefficient is

$$\left(\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right)=\left(\begin{array}{ccc}54.221& 545.814& 2739.960\end{array}\right)$$

When the air pressure is 0.5 MPa, the polynomial fitting coefficient is

$$\left(\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right)=\left(\begin{array}{ccc}77.646& 561.740& 2993.145\end{array}\right)$$

When the air pressure is 0.8 MPa, the polynomial fitting coefficient is

$$\left(\begin{array}{ccc}{a}_1& a_2& a_3\end{array}\right)=\left(\begin{array}{ccc}102.111& 911.286& 5951.281\end{array}\right)$$

At this time, the curves of nonlinear restoring force and displacement of air spring are shown in Figure 3.

Figure 3

Nonlinear restoring force-displacement curve of air spring.

3.2. Verification of polynomial model

For various types of air springs, manufacturers will give the performance curve of rubber air springs through test data. Guizhou Qianjin Rubber Co., Ltd. has given the performance curve of Qianjin brand 1884N rubber air spring as shown in Figure 1. ⁵⁰ Because Figure 1 is a performance curve of mass and displacement, the performance curve of load force and displacement of air spring can be obtained. The polynomial curve of the air spring is compared with the experimental curve as shown in Figure 4. The residual values are calculated, as shown in Table 4. From the calculated residual values, we can see that the residual values are low and the polynomial model curve can well predict the actual performance curve of the air spring, which verifies the correctness of the polynomial model.

Figure 4

Comparison diagram of fitting curves. Note: PFC: polynomial fitting curve; RBFFC: radial basis function fitting curve; TDC: test data curve.

Table 4

Polynomial fitting residuals.

Pressure/MPa	0.3	0.5	0.8
Residual/KN	0.30	0.89	0.36

3.3. Comparison with radial basis function fitting

Radial basis function (RBF) has good nonlinear approximation performance and convenient training. It has been widely used in the signal process, structure analysis, control decision system, and other disciplines and engineering fields and is especially suitable for solving curve-fitting problems. The main function used in the design of radial basis function neural network in MATLAB neural network toolbox is newrb. The paper will use this function to fit the function of air spring force, as shown in Figure 5. At the same time, compared with the polynomial fitting model curve, as shown in Figure 4 and calculated the residual value of radial basis function fitting, as shown in Table 5. It is known that the residual value of radial basis function fitting is larger than that of polynomial fitting in Tables 4 and 5. So polynomial fitting is better than radial basis function fitting. The polynomial fitting function will be used in the next modeling process.

Figure 5

Radial basis function fitting curve.

Table 5

Radial basis function fitting residuals.

Pressure/MPa	0.3	0.5	0.8
Residual/KN	0.72	2.05	0.87

4. Dynamic modeling of nonlinear air suspension system

In the paper, the active air suspension system of nine degrees vehicle-seat-body model of vehicle is shown in Figure 6.

Figure 6

Active air suspension system of nine degrees vehicle-seat-body model of vehicle.

In the picture,

$m_h$ 、 $m_c$ 、 $m_s$ —Human body 、 seat and vehicle body mass,

$m_{u1}$ 、 $m_{u2}$ 、 $m_{u3}$ 、 $m_{u4}$ —Front and rear unsprung mass,

$k_h$ 、 $k_c$ --—Seat cushion and seat equivalent stiffness,

$k_{s1}$ 、 $k_{s2}$ 、 $k_{s3}$ 、 $k_{s4}$ -—Stiffness coefficient of front and rear suspension air springs,

$k_{u1}$ 、 $k_{u2}$ 、 $k_{u3}$ 、 $k_{u4}$ —Stiffness coefficient of front and rear tires,

$c_h$ 、 $c_c$ —Seat cushion and seat equivalent damping,

$c_{s1}$ 、 $c_{s2}$ 、 $c_{s3}$ 、 $c_{s4}$ —Damping coefficient of front and rear suspension,

$c_{u1}$ 、 $c_{u2}$ 、 $c_{u3}$ 、 $c_{u4}$ —Damping coefficient of front and rear tire,

$z_h$ 、 $z_c$ —Vertical displacement of body and seat,

$a$ 、 $b$ _{—Distance between front and rear axle to mass center of vehicle body,}

$e$ 、 $d$ — The vertical distance from the seat to the X and Y axes.

$w$ — 1/2 wheelbase;

$q_1$ 、 $q_2$ 、 $q_3$ 、 $q_4$ —Road roughness input of front and rear wheel,

$z_s$ —Vertical displacement at mass center of vehicle body.

$\theta$ —Pitch angle at body centroid,

$\phi$ —Roll angle at body centroid,.

$z_{s1}$ 、 $z_{s2}$ 、 $z_{s3}$ 、 $z_{s4}$ —Vertical displacement of front and rear sprung mass,

$z_{u1}$ 、 $z_{u2}$ 、 $z_{u3}$ 、 $z_{u4}$ —Vertical displacement of front and rear unsprung mass,

Assuming that the body is a rigid body, when the vehicle roll angle and pitch angle are small, body displacement above the suspension can be expressed as

$$\{\begin{array}{l}{z}_{s1}=z_s-a\theta +w\phi \hfill \\ z_{s2}=z_s-a\theta -w\phi \hfill \\ z_{s3}=z_s+b\theta +w\phi \hfill \\ z_{s2}=z_s+b\theta -w\phi \hfill \\ z_e=z_s-e\theta +d\phi .\hfill \end{array}$$

(10)

According to Newton’s law of motion, the differential equation of dynamics are as follows

$$m_h{\ddot{z}}_h+f_{kh}+f_{ch}=0$$

(11)

$$m_c{\ddot{z}}_c+f_{kc}+f_{cc}-f_{kh}-f_{ch}=0$$

(12)

$$m_s{\ddot{z}}_s+f_{ks1}+f_{cs1}-u_1+f_{ks2}+f_{cs2}-u_2+f_{ks3}+f_{cs3}-u_3+f_{ks4}+f_{cs4}-u_4-f_{kc}-f_{cc}=0$$

(13)

$$I_{\theta }\ddot{\theta }+a\left(f_{ks1}+f_{cs1}-u_1+f_{ks2}+f_{cs2}-u_2\right)+e\left(f_{kc}+f_{cc}\right)-b\left(f_{ks3}+f_{cs3}-u_3+f_{ks4}+f_{cs4}-u_4\right)=0$$

(14)

$$I_{\phi }\ddot{\phi }+w\left(f_{ks1}+f_{cs1}-u_1+f_{ks3}+f_{cs3}-u_3\right)-d\left(f_{kc}+f_{cc}\right)-w\left(f_{ks2}+f_{cs2}-u_2+f_{ks4}+f_{cs4}-u_4\right)=0$$

(15)

$$m_{u1}{\ddot{z}}_{u1}+f_{ku1}+f_{cu1}+u_1-f_{ks1}-f_{cs1}=0$$

(16)

$$m_{u2}{\ddot{z}}_{u2}+f_{ku2}+f_{cu2}+u_2-f_{ks2}-f_{cs2}=0$$

(17)

$$m_{u3}{\ddot{z}}_{u3}+f_{ku3}+f_{cu3}+u_3-f_{ks3}-f_{cs3}=0$$

(18)

$$m_{u4}{\ddot{z}}_{u4}+f_{ku4}+f_{cu4}+u_4-f_{ks4}-f_{cs4}=0$$

(19)

in which $f_{kh}=k_h\left(z_h-z_c\right)$ , $f_{ch}=c_h\left({\dot{z}}_h-{\dot{z}}_c\right)$ —Elastic force and damping force of cushion; $f_{kc}=k_c\left(z_c-z_e\right)$ , $f_{cc}=c_c\left({\dot{z}}_c-{\dot{z}}_e\right)$ —Seat spring force and damping force;

• $f_{ks1}=a_1\left(z_{s1}-z_{u1}\right)+a_2{\left(z_{s1}-z_{u1}\right)}^2+a_3{\left(z_{s1}-z_{u1}\right)}^3$ _, $f_{ks2}=a_1\left(z_{s2}-z_{u2}\right)+a_2{\left(z_{s2}-z_{u2}\right)}^2+a_3{\left(z_{s2}-z_{u2}\right)}^3$ _, $f_{ks3}=a_1\left(z_{s3}-z_{u3}\right)+a_2{\left(z_{s3}-z_{u3}\right)}^2+a_3{\left(z_{s3}-z_{u3}\right)}^3$ _, $f_{ks4}=a_1\left(z_{s4}-z_{u4}\right)+a_2{\left(z_{s4}-z_{u4}\right)}^2+a_3{\left(z_{s4}-z_{u4}\right)}^3$ —Front and rear nonlinear air spring force;

• $f_{cs1}=c_{s1}\left({\dot{z}}_{s1}-{\dot{z}}_{u1}\right)$ _, $f_{cs2}=c_{s2}\left({\dot{z}}_{s2}-{\dot{z}}_{s2}\right)$ _, $f_{cs3}=c_{s3}\left({\dot{z}}_{s3}-{\dot{z}}_{u3}\right)$ _, $f_{cs4}=c_{s4}\left({\dot{z}}_{s4}-{\dot{z}}_{s4}\right)$ —Front and rear damping force of shock absorber;

• $f_{ku1}=k_{u1}\left(z_{u1}-q_1\right)$ _, $f_{ku2}=k_{u2}\left(z_{u2}-q_2\right)$ _, $f_{ku3}=k_{u3}\left(z_{u3}-q_3\right)$ _, $f_{ku4}=k_{u4}\left(z_{u4}-q_4\right)$ —Front and rear elastic force of tire;

• $f_{cu1}=c_{u1}\left({\dot{z}}_{u1}-{\dot{q}}_1\right)$ _, $f_{cu2}=c_{u2}\left({\dot{z}}_{u2}-{\dot{q}}_2\right)$ _, $f_{cu3}=c_{u3}\left({\dot{z}}_{u3}-{\dot{q}}_3\right)$ _, $f_{cu4}=c_{u4}\left({\dot{z}}_{u4}-{\dot{q}}_4\right)$ —Front and rear damping force of tire;

• $u_1$ 、 $u_2$ 、 $u_3$ 、 $u_4$ —Control force generated by the front and rear suspension active controllers.

5. Design of fuzzy fractional order PI^λD^μ control system

5.1. Structure of fuzzy fractional order PI^λD^μ controller

The paper selects the vertical displacement of front and rear sprung mass $z_{sf}$ and $z_{sr}$ as the feedback input. Fuzzy controller takes its error $e$ and its change rate $ec$ as input, and correction parameters $\mathrm{\Delta }{k}_p$ 、 $\mathrm{\Delta }{k}_i$ 、 $\mathrm{\Delta }{k}_d$ as output. The fractional order PI^λD^μ controller output parameters $k_p$ 、 $k_i$ 、 $k_d$ are shown in equations 19–21. The schematic diagram of fuzzy fractional order PI^λD^μ control structure is shown in Figure 7.

$$k_p={k^{\prime }}_p+\mathrm{\Delta }{k}_p$$

(20)

$$k_i={k^{\prime }}_i+\mathrm{\Delta }{k}_i$$

(21)

$$k_d={k^{\prime }}_d+\mathrm{\Delta }{k}_d$$

(22)

in which ${k^{\prime }}_p$ 、 ${k^{\prime }}_i$ 、 ${k^{\prime }}_d$ are the initial value.

Figure 7

Fuzzy adaptive fractional order PI^λD^μ control principle diagram.

The output of the fractional order PI^λD^μ controller is a combination of input multiple, fractional derivative, and fractional integral. The formula is

$$u\left(t\right)=K_{p}e\left(t\right)+K_i{D}_t^{-\lambda }e\left(t\right)+K_d{D}_t^{\mu }e\left(t\right).$$

(23)

in which $\lambda$ 、 $\mu$ —Integral and differential power of fraction.

5.2. Fuzzy subset of input and output variables

For input and output variables, fuzzy state is used by seven fuzzy sets, both for the $\left\{NB,NM,NS,O,PS,PM,PB\right\}$ . The fuzzy domain of error $e$ and error change rate $ec$ are in (−3,3). The fuzzy domain of $\mathrm{\Delta }{K}_p$ 、 $\mathrm{\Delta }{K}_i$ 、 $\mathrm{\Delta }{K}_d$ is in (−1,1). Fuzzy subsets of input variable and output variable use triangular membership functions.^51,52

5.3. Fuzzy control rule

The fuzzy control rules is an important component of fuzzy controller. It is described by means of the input and output of a controller relationship between them, namely, the fuzzy relationship between them. Although the method of fuzzy reasoning is given, the fuzzy control rule is based on people’s Cognition. In the paper, the fuzzy controller with two inputs is described by two language fuzzy sets, forming 49 control rules. The fuzzy control rule table of three parameters is shown in Table 6.^53,54

Table 6.

Fuzzy control rules of different parameter.

$e$ $ec$	NB	NM	NS	O	PS	PM	PB
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ NB $\mathrm{\Delta }{K}_d$	PB NB PS	PB NB PS	PM NB O	PM NM O	PS NM PS	PS O PM	O O PB
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ NM $\mathrm{\Delta }{K}_d$	PB NB NS	PB NB NS	PM NM NS	PM NM NS	PS NS O	O O PM	O O PM
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ NS $\mathrm{\Delta }{K}_d$	PM NM NB	PM NM NM	PM NS NM	PS NS NS	O O O	NS PS PS	NM PS PM
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ O $\mathrm{\Delta }{K}_d$	PM NM NB	PS NS NM	PS NS NM	O O NS	NS PS O	NM PS PS	NM PM PM
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ PS $\mathrm{\Delta }{K}_d$	PS NS NB	PS NS NM	PS O NS	NS PS NS	NS PS O	NM PM PS	NM PM PS
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ PM $\mathrm{\Delta }{K}_d$	O O NM	O O NS	NS PS NS	NM PM NS	NM PM O	NM PB PS	NB PB PS
$\mathrm{\Delta }{K}_p$ $\mathrm{\Delta }{K}_i$ PB $\mathrm{\Delta }{K}_d$	O O PS	NS O O	NS PS O	NM PM O	NM PM O	NB PB PB	NB PB PB

5.4. Comparative analysis of simulation

In order to show the superiority of the controller designed in this paper, PID controller, fuzzy controller, and fuzzy PID controller are also designed. We take transfer function as

$$G\left(s\right)=\frac{1}{s^2+s+1}$$

(24)

Calculating the unit step response and getting the simulation curve as shown in Figure 8. USS—unit step signal curve, USP—unit step response curve, PIDC—PID control curve, Fuzzy—fuzzy control curve, Fuzzy-PID—Fuzzy-PID control curve, and Fuzzy-PI^λD^μ—Fuzzy fractional order PI^λD^μ control curve.

Figure 8

Unit step response curve.

From Figure 8, it can be seen that the fuzzy fractional order PI^λD^μ control system has the best control effect by combining the overshoot, reaction time, stability time, delay time, rise time, steady-state error, and other factors.

6. Simulation analysis

When the vehicle travels on grade C road and $u=20\hspace{0.17em}m/s$ , the formula for follows road roughness coefficient is as

$$S_q\left(n_0\right)=256\times {10}^{-6}\hspace{0.17em}{m}^2/m^{-1}\hspace{0.17em}\left(n_0=0.1\hspace{0.17em}{m}^{-1}\right)$$

(25)

Velocity spectrum density is a white noise

$$S_{\dot{q}}\left(f\right)={\left(2\pi f\right)}^2{S}_q\left(f\right)=4{\pi }^2{n}_0^2{S}_q\left(n_0\right)u$$

(26)

The white noise input is regarded as random road displacement excitation. The road condition of the same side tire is basically the same and the rear wheel input is only a little delayed. The time difference is related to the wheelbase and speed, then the time difference formula is $\mathrm{\Delta }t=L/u$ , $L$ —the wheelbase, and $u$ —the speed of the vehicle. The front and rear wheel input response can be obtained, as shown in Figure 9.

Figure 9.

Input response of front and rear wheel road. Note: (a) front left wheel, (b) frond right wheel, (c) rear left wheel, (d) rear right wheel.

Wheel dynamic load, suspension dynamic deflection, vertical acceleration of vehicle body, vertical acceleration of seat, vertical acceleration of human body, pitch angle acceleration, and roll angle acceleration are taken as control objects. Through trial and error, we can get the λ=0.1, μ=0.5 in the fuzzy fractional order PI^λD^μ control strategy. From the evaluation index curves in Figures 10–22 and the root means square values in Table 7, we can see that the root mean square values of wheel dynamic load of active nonlinear air suspension controlled by fuzzy fractional order PI^λD^μ on grade C pavement is better than that of passive suspension, PID control suspension, fuzzy control suspension, and fuzzy PID control suspension. The root means square values of vertical acceleration of vehicle body, seat, human body, and pitch angle acceleration are also improved than that of another control suspension. Although the front and rear suspensions have larger dynamic deflection than the fuzzy PID control, they are still greatly improved compared with the passive air suspension. Therefore, the comprehensive performance of the active air suspension controlled by fuzzy fractional order PI^λD^μ is better and the ride comfort and handling stability of the vehicle are effectively improved.

Figure 10

Curve of the front left wheel dynamic load under full load.

Figure 11.

Curve of the front right wheel dynamic load under full load.

Figure 12.

Curve of the rear left wheel dynamic load under full load.

Figure 13.

Curve of the rear right wheel dynamic load under full load.

Figure 14.

Curve of the front left suspension dynamic deflection under full load.

Figure 15.

Curve of the front right suspension dynamic deflection under full load.

Figure 16.

Curve of the rear left suspension dynamic deflection under full load.

Figure 17.

Curve of the rear right suspension dynamic deflection under full load.

Figure 18.

Curve of the vehicle body vertical acceleration under full load.

Figure 19.

Curve of the pitch angle acceleration under full load.

Figure 20.

Curve of the roll angle acceleration under full load.

Figure 21.

Curve of the seat vertical acceleration under full load.

Figure 22.

Curve of the human body vertical acceleration under full load.

Table 7

Performance comparison of vehicle under full load and different control.

Performance evaluation index (RMS)	Passive suspension	PID	Performance improvement, %	Fuzzy	Performance improvement, %	Fuzzy-PID	Performance improvement, %	Fuzzy-PIλDμ	Performance improvement, %
Front left wheel dynamic load(N)	1865.2854	1549.9896	16.90	1691.8529	9.30	1429.3445	23.37	1391.8978	25.38
Front right wheel dynamic load(N)	2196.6896	1886.9756	14.10	2038.0913	7.22	1772.1874	19.32	1691.6869	22.99
Rear left wheel dynamic load(N)	2259.6805	1931.4612	14.53	2089.8354	7.52	1804.0534	20.16	1692.9434	25.08
Rear right wheel dynamic load(N)	2135.6268	1778.8837	16.70	1954.9036	8.46	1641.5605	23.13	1588.4174	25.62
Front left suspension dynamic deflection(m)	0.0197	0.0163	17.26	0.0182	7.61	0.0148	24.87	0.0153	22.34
Front right suspension dynamic deflection(m)	0.0226	0.0191	15.49	0.0210	7.08	0.0180	20.35	0.0175	22.57
Rear left suspension dynamic deflection(m)	0.0226	0.0191	15.49	0.0210	7.08	0.0180	20.35	0.0175	22.57
Rear right suspension dynamic deflection(m)	0.0228	0.0190	16.67	0.0214	6.14	0.0173	24.12	0.0179	21.49
Vertical acceleration of vehicle body (m/s2)	0.4121	0.2974	27.83	0.3268	20.70	0.2621	36.40	0.2021	50.96
Pitch angle acceleration (rad/s2)	0.1379	0.1086	21.25	0.1136	17.62	0.0998	27.63	0.0792	42.57
Roll angle acceleration (rad/s2)	0.1425	0.0962	32.49	0.1032	27.58	0.0869	39.02	0.0729	48.84
Vertical acceleration of seat (m/s2)	0.5326	0.4007	24.77	0.4347	18.38	0.3576	32.86	0.2806	47.32
Vertical acceleration of human body (m/s2)	0.8283	0.6564	20.75	0.6899	16.71	0.5942	28.26	0.4517	45.47

7. Conclusion

The analysis of air spring stiffness characteristics and restoring force show that air spring stiffness as well as restoring force are quadratic and cubic polynomials of the displacement. The air spring polynomial model used in this paper can better predict the actual performance curve of the air spring through the residual analysis calculation. Through modeling analysis, the pavement excitation curves were obtained by using Matlab/Simulink at a vehicle speed of 20 km/h on a Class C road surface. By designing the fuzzy fractional order PI^λD^μ control and simulating the dynamic load, dynamic deflection, and pitch angle of each wheel and suspension using Simulink, we can see that the fuzzy fractional order PI^λD^μ control strategy largely attenuates the vertical acceleration of the body, seat, and human body quite high in terms of ride, compared with the fuzzy PID control strategy, PID control strategy, and fuzzy control strategy. In addition, it is also outstanding in reducing the pitch angle and roll angle, which improves the steering stability and smoothness of the car.

In this paper, only a common input of white noise class C road is considered, which is not good enough to illustrate the comprehensive performance of fuzzy fractional order PI^λD^μ control. In addition, when considering the vertical acceleration of the human body, only the single rider case is considered, and the results may deviate slightly from the actual situation. We will make up for the previously considered situation in the future. With the advancement of intelligence, the air suspension will also develop in the direction of active “intelligence."