Research on air suspension with novel dampers based on glowworm swarm optimization proportional–integral–derivative algorithm

Abstract

In order to enhance automotive ride performance, permanent magnets and magnetic valves were introduced in magnetorheological dampers and hence permanent magnets and magnetic valves introduced in magnetorheological damper’s air suspension was studied. First, on the basic principle of providing magnetic field by permanent magnets and adjusting damping force by magnetic valves, a novel permanent magnets and magnetic valves introduced in magnetorheological damper was designed. On the basis of equivalent surface current model, magnetic induction intensity activated by permanent magnets in damping channel was calculated, and damping force model was built based on hydromechanics theory. Damping force–displacement experiments and damping force–velocity experiments were carried out, and the theoretical model of permanent magnets and magnetic valves introduced in magnetorheological damper was verified. Then, physical and mathematical models of 4 degree-of-freedom permanent magnets and magnetic valves introduced in magnetorheological damper air suspension model were built. In order to enhance the controlling effects of the studied air suspension, optional moving operation, moving variation operation and operation of step size changing with fluorescein proportionality were introduced into glowworm swarm optimization algorithm, and glowworm swarm optimization proportional–integral–derivative controlling algorithm was designed. Finally, three experiments with input of impact road surface, sine wave road surface, and white noise road surface were carried out, and the experimental results verified that the working performance of permanent magnets and magnetic valves introduced in magnetorheological damper air suspension based on glowworm swarm optimization proportional–integral–derivative controlling algorithm was good.

Introduction

Automotive suspension systems, with a great impact on riding performance and handling performance, have always been one of the studying focuses, so that suspension systems have been studied from different aspects by scholars. At present, many scholars focus their research on air suspension systems. For example, Nieto et al. ¹ developed an air suspension with two connecting pipes of different sizes: the riding comfort of automobiles could be improved if the shorter pipe was used, while the handling performance could be improved if the longer one was used. Eskandary et al. ² made two rubber airbags for an air suspension. By controlling the pressure in the chambers of the two airbags, automotive ride height and stiffness can be adjusted, which could enhance automotive suspension performance. Sun et al. ³ developed a newly designed controller for an air suspension. Through the controller, the automotive height could be adjusted so that the roll and pitch angle could be regulated. Li et al. ⁴ established a model car test bench with interconnected air suspension, and research on the torsion elimination performance was conducted through the test bench. Abid et al. ⁵ investigated an air spring suspension equivalent system and optimization methods of obtaining good parameters of air spring system. Chen et al. ⁶ developed a height adjusting system for the air suspension of agriculture transport vehicles. With the help of designed single neuron-adaptive proportional–integral–derivative (PID) control algorithm, the adjusting system had many advantages, such as shorter response time and faster response speed. Li and Li ⁷ studied an effective vertical stiffness formula for air spring used by automotive air suspensions, and four experimental methods were employed to study static and quasi-static vertical stiffness. Wong et al. ⁸ developed a numerical model of the rolling lobe air spring, and in the model, finite element method was used. The accuracy of the numerical model was proved to be very good. Zargar et al. ⁹ developed a nonlinear mathematical model for air springs. Experiments were carried out to verify the model and the experimental results indicated that the model was effective. Meanwhile, many scholars focus their study on other aspects of air suspensions and have abstained fruitful academic achievements. ¹⁰

Moreover, in order to enhance suspension working performance, many scholars focus their study on magnetorheological damper (MRD) suspensions. For instance, Fakhraei et al. ¹¹ conducted chaotic vibration analysis of an articulated vehicle with MRDs, and the analysis results could help to design heavy vehicles with MRDs. Hu et al. ¹² built a quarter-car model with an MRD and designed a novel hybrid fuzzy controller for the suspension which could enhance the riding comfort of automobiles. Tang et al. ¹³ developed a state-observer-based Takagi–Sugeno fuzzy controller for a semi-active suspension with MRDs, and the controller was verified by a designed bench testing system. Balamurugan et al. ¹⁴ studied a novel modified algebraic model for MRDs, through which a semi-active quarter suspension model was built. Then two nested controllers with good controlling effects were studied for the suspension system. Ata and Salem ¹⁵ investigated the sky-hook, hybrid, and fuzzy-hybrid controllers for tracked vehicle with MRDs, with their advantages and disadvantages being discussed through experimental analysis. Sulaiman et al. ¹⁶ developed a semi-active controller for a light-heavy vehicle suspension with MRDs, and the controlling effect on tire force control was explored. Tyan et al. ¹⁷ built a controller for analyzing the working performance of semi-active suspension with MRDs, and through simulation experiments, the controller’s working effect was verified. At the same time, many scholars have researched MRD suspensions from other aspects. ¹⁸

Based on the literature mentioned above, a conclusion can be drawn that the working performance of automotive suspensions can be enhanced when either air springs or MRDs is used in suspensions. Meanwhile, the literature also indicates in MRD suspensions study, the suspension springs are regarded as normal ones, while in air suspensions study, the suspension dampers are regarded as normal ones. In order to enhance comprehensive performance, some scholars introduced permanent magnets into MRDs. Dong ¹⁹ investigated an MRD with permanent magnet used by automobiles, and numerical simulation results showed that automotive suspensions with the newly designed MRD had good ride comfort. Kim et al. ²⁰ designed a new MRD, whose damping force could be tuned by the permanent magnet. Experiments were conducted, and it was demonstrated that the damping force of the proposed MRD could be tuned rapidly. Hu et al. ²¹ proposed a self-sensing MRD composed of permanent magnets and coils, and experimental results verified that the self-sensing MRD had good self-sensing ability and controllable damping capability. Zhang et al. ²² proposed a novel MRD with permanent magnets which could solve fail-safe problem of normal dampers. Du et al. ²³ investigated an MRD with permanent magnets which could solve the problem of magnetorheological fluid sedimentation of normal MRDs. The above literature indicates that comprehensive performance can be enhanced greatly through introducing permanent magnets into MRDs. Therefore, suspensions with MRDs embedded with permanent magnets will be the developing trend of automotive suspensions. Xiao et al. ²⁴ investigated a suspension with MRDs embedded with permanent magnets, but the suspension is traditional one with normal springs. However, the research on air suspensions with MRDs embedded with permanent magnets is rare to be seen, and there is a lot of research room. In order to enhance the working performance of air suspensions, a permanent magnets and magnetic valves introduced magnetorheological damper (PMMVMRD) was designed, and its force models were built. Then, on this basis, PMMVMRD air suspension models were built. Meanwhile, in order to enhance controlling effects of air suspensions, glowworm swarm optimization (GSO) algorithm was developed, and on this basis, glowworm swarm optimization proportional–integral–derivative (GSO-PID) controlling algorithm suitable for PMMVMRD air suspension was designed. Finally, simulations were performed for PMMVMRD air suspension based on GSO-PID controlling algorithm in three road inputs. In addition, a new type of PMMVMRD, which can be adjusted through a magnetic valve and magnetic field generated by solenoid coils and a permanent magnet, will be investigated in future. When the difference between actual damping force and desired damping force is large, adjusting method of magnetic valve can be used. When their deference is small, adjusting method of controlled magnetic field generated by solenoid coils and a permanent magnet can be used so that accurate control can be realized. If the new type of PMMVMRD can be used in air suspensions, riding comfort and handling performance can be possibly further enhanced. This article is organized as follows. In the coming five sections, mathematical models of proposed PMMVMRD are built, and its experiments are carried out as well. Next, GSO algorithm is developed, and GSO-PID controlling algorithm is designed. Subsequently, models of PMMVMRD air suspension are built. Then, simulation experiments are carried out by taking impact road surface and sine wave road surface as inputs. Finally, conclusions are drawn in the last section.

Mathematical models and experiments of PMMVMRD

Mechanic structure and working principle of PMMVMRD

By referring to MRD mechanical structure in literature 24, the piston head shown in Figure 1 is I shaped, and a circular permanent magnet (black end is N pole and white pole is S pole) is mounted in the recess. Damping channel forms between the outer surface of the piston and the inner surface of the cylinder barrel of damper. A stainless-steel resistance-magnet plate is mounted on each side of the piston so that magnetic field direction can be perpendicular to damping channel, whose magnetic field is caused only by the circular permanent magnet. There are some components in the piston head, such as a magnetic coil, a valve plug, a back spring, a screwed pipe, and magnetic valve pipes (a bigger pipe is on the right side and three symmetric circular pipes are on the left sides), which make up a magnetic valve. When no electricity is applied in the magnetic coil, magnetic valve pipe will be shut off by the back spring, while electricity is applied, the attractive force generated by the coil will overcome the elastic force to open the magnetic valve. The opening degree will be bigger with the electricity current increasing. Thus, the working principle of PMMVMRD is as follows. When electricity in the magnetic coil is 0, the valve will shut off, and all the magnetorheological liquid will flow through the damping channel, and hence strong damping force will be created. When electricity with certain strength is applied in the coil, valve will be opened to corresponding degree, some magnetorheological liquid will flow through magnetic valve, and the remaining will flow through the damping channel. Damping force will not be created by the liquid flowing through magnetic valve, but only by that through damping channel. Therefore, damping force can be adjusted through regulating the open degree of magnetic valve.

Figure 1.

Mechanical-structure piston head of PMMVMRD.

Prediction model of damping force

According to the mechanical structure of PMMVMRD, the diagram of magnetic circuit can be achieved and shown in Figure 2, which indicates that the line of magnetic force is along the route A₁→A₂→A₃→A₄→A₃→A₂→A₁ in a circle.

Figure 2.

Diagram of magnetic circuit.

Suppose that total flow rate is Q in the piston movement, so equation (1) can be achieved as

$$Q=Q_1+Q_2$$

(1)

where Q₁ and Q₂ are flow rates through magnetic valve and damping channel, respectively.

Total flow rate Q can be calculated by equation (2)

$$Q=A_p·v$$

(2)

where A_p and v are the piston’s sectional area and relative velocity of piston, respectively.

A_p can be calculated by equation (3)

$$A_p=\pi {r_2}^2$$

(3)

Suppose that the input current in the magnetic coil is I, so air-gap field strength H can be calculated through equation (4)

$$H=\frac{\mathrm{NI}}{x\prime }$$

(4)

where N and x′ are coil turn number and moving displacement of valve plug, respectively.

Suppose that the magnetic field energy in air-gap field is E, the reduced amount dE in the air-gap field can be expressed by equation (5) according to the working principle of electromagnets

$$\mathrm{dE}=\frac{1}{2}{H}^2{\mu }_0\pi {r_5}^2\mathrm{dx}\prime$$

(5)

where μ₀ is magnetic permeability, and r₅ is the cross-sectional radius of the iron core.

According to energy conservation law, the reduced magnetic energy turns into elastic potential energy of magnetic valve. Thus, equation (6) can be achieved

$$\frac{1}{2}{H}^2{\mu }_0\pi {r_5}^2{L}_{\mathrm{v}}=\frac{1}{2}{k}_a(L_{\mathrm{v}}{)}^2$$

(6)

where k_a is spring stiffness of magnetic valve; L_v is maximum displacement of valve plug.

Equation (7) can be achieved by substituting equation (4) into equation (6)

$$x\prime =\sqrt{\frac{{\mu }_0\pi N^2{I}^2{r_5}^2}{k_a{L}_{\mathrm{v}}}}$$

(7)

Because the right end of the valve plug is a cone, suppose that the height and radius of the cone are $h\prime$ and $r\prime$ , respectively. Suppose that the radius of the bigger pipe on the right side of magnetic valve is $r″$ . Therefore, the section area $S″$ of flowing liquid with magnetic valve control can be calculated through equation (8) when valve plug moves to the left for a displacement of $x\prime$

$$S″=\pi {(r″-\frac{x\prime r\prime }{h\prime })}^2$$

(8)

Thus, according to the related theory of fluid mechanics, the flow rate Q₁ can be calculated through equation (9) when the input current in the magnetic valve is I

$$Q_1=\pi v{(r″-\frac{\sqrt{\frac{{\mu }_0\pi N^2{I}^2{r_5}^2}{k_a{L}_{\mathrm{v}}}}r\prime }{h\prime })}^2$$

(9)

Bingham model is employed to calculate damping force of PMMVMRD in this article. Therefore, magnetic induction intensity in damping channel must be calculated first. Magnetic induction intensity in the damping channel is only activated by permanent magnet in the studied MRD. In order to calculate the magnetic induction intensity accurately, equivalent surface current model is employed to conduct the calculation. Because the structure of permanent magnet is a hollow cylinder, rectangular coordinate system XOYZ is built by taking center and axis of cylinder as origin and Z axis, respectively. Because the cylinder is of symmetric structure about Y axis, only XOZ plane needs to be taken into consideration, and XOZ coordinate system is shown in Figure 3(a). As is shown in Figure 3(a), equivalent surface current of the permanent magnet is surface current 1, 2, 3, and 4, which can be regarded as superposition of several endless line-current elements. Suppose that any point coordinate in damping channel is P(x, y, z), then, point P(x, y, z) magnetic field activated by permanent magnet is the value that P(x, y, z) magnetic field activated by equivalent surface current 1 and 2 minus that activated by 3 and 4, according to model of equivalent surface current. Point P(x, y, z) magnetic field activated by line-current element of surface current 1 is shown in Figure 3(b). Because the permanent magnet is of symmetric structure about Y axis, make y = 0 and suppose that coordinate of any point on line-current element of surface current 1 is P(x″,y″,z″). Likewise, make y″ = 0 and point P(x, y, z) magnetic field created by the line-current element can be expressed by equation (10)

$$\{\begin{array}{c}d{B}_{1X}=d\left(\frac{{\mu }_{0}I\prime }{2\pi r_1\prime }\sin \theta \right)\hfill \\ d{B}_{1Y}=d(-\frac{{\mu }_{0}I\prime }{2\pi r_1\prime }\cos \theta )\hfill \end{array}$$

(10)

where $I\prime$ is surface current strength, $r_1\prime$ is length of radius vector from electricity current element to point P(x, y, z) and $\theta$ is an angle indicated in Figure 3.

Figure 3.

Diagram of permanent magnetic field: (a) diagram of permanent magnet surface current; (b) diagram of point P(x, y, z) magnetic field activated by line-current element of surface current 1.

Equation (11) can be achieved according to geometrical relationship in Figure 3(b)

$$\{\begin{array}{c}\theta =\arcsin \left(\frac{z-z″}{r_1\prime }\right)\hfill \\ r_1\prime =\sqrt{{(x″-x)}^2+{(y″-y)}^2}\hfill \end{array}$$

(11)

Suppose that $k_v$ is line-current density, so the magnetic field created by surface current $I\prime$ can be achieved through equation (12)

$$\{\begin{array}{c}{B}_{1X}=\underset{{\theta }_1}{\overset{{\theta }_2}{\int }}\frac{{\mu }_0{k}_v}{2\pi r_1\prime }\sin \theta \mathrm{dz}\prime \hfill \\ B_{1Y}=-\underset{{\theta }_1}{\overset{{\theta }_2}{\int }}\frac{{\mu }_0{k}_v}{2\pi r_1\prime }\cos \theta \mathrm{dz}\prime \hfill \end{array}$$

(12)

where $B_{1X}$ and $B_{1Y}$ are X-direction component and Y-direction component of magnetic induction intensity, respectively.

Equation (13) can be achieved according to equations (11) and (12)

$$\{\begin{array}{c}{B}_{1X}=\underset{{\theta }_1}{\overset{{\theta }_2}{\int }}\frac{{\mu }_0{k}_v}{2\pi \sqrt{{(x″-x)}^2+{(y″-y)}^2}}\sin \left(\arcsin \left(\frac{z-z″}{r_1\prime }\right)\right)\mathrm{dz}\prime \hfill \\ B_{1Y}=-\underset{{\theta }_1}{\overset{{\theta }_2}{\int }}\frac{{\mu }_0{k}_v}{2\pi \sqrt{{(x″-x)}^2+{(y″-y)}^2}}\cos \left(\arcsin \left(\frac{z-z″}{r_1\prime }\right)\right)\mathrm{dz}\prime \hfill \end{array}$$

(13)

Equation (14) can be achieved according to geometrical relationship

$$\{\begin{array}{c}{r}_1\prime =\frac{x-r_2/2}{\cos \theta }\hfill \\ {\theta }_1=\arctan \left(\frac{z+L/2}{x-r_2/2}\right)\hfill \\ {\theta }_2=\arctan \left(\frac{z-L/2}{x-r_2/2}\right)\hfill \\ \mathrm{dz}\prime =-\frac{x-r_2/2}{{\cos }^2\theta }d\theta \hfill \end{array}$$

(14)

In damping channel, point P(x, y, z) magnetic flux density created by surface current 1 can be calculated through equation (15)

$$\{\begin{array}{c}{B}_{1X}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\ln \frac{{(z+L/2)}^2+{(x-r_2/2)}^2}{{(z-L/2)}^2+{(x-r_2/2)}^2}\hfill \\ B_{1Z}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }[\arctan \left(\frac{z-L/2}{x-r_2/2}\right)-\arctan \left(\frac{z+L/2}{x+r_2/2}\right)]\hfill \end{array}$$

(15)

Equation (16) can be achieved through equation (15)

$$\{\begin{array}{c}{B}_{1X}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\ln [(z+L/2{)}^2+(x-r_2/2{)}^2]-\frac{{\mu }_0{k}_v}{4\pi }\ln [(z-L/2{)}^2+(x-r_2/2{)}^2]\hfill \\ B_{1Z}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z-L/2}{x-r_2/2}\right)-\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z+L/2}{x+r_2/2}\right)\hfill \end{array}$$

(16)

Likewise, in damping channel, point P(x, y, z) magnetic induction density created by surface current 2, 3, and 4 can be calculated through equations (17)–(19)

$$\{\begin{array}{c}{B}_{2X}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\ln [(z+L/2{)}^2+(x+r_2/2{)}^2]-\frac{{\mu }_0{k}_v}{4\pi }\ln [(z-L/2{)}^2+(x+r_2/2{)}^2]\hfill \\ B_{2Z}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z+L/2}{x+r_2/2}\right)-\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z-L/2}{x+r_2/2}\right)\hfill \end{array}$$

(17)

$$\{\begin{array}{c}{B}_{3X}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\ln [(z+L/2{)}^2+(x-r_1/2{)}^2]-\frac{{\mu }_0{k}_v}{4\pi }\ln [(z-L/2{)}^2+(x-r_1/2{)}^2]\hfill \\ B_{3Z}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z-L/2}{x-r_1/2}\right)-\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z+L/2}{x+r_1/2}\right)\hfill \end{array}$$

(18)

$$\{\begin{array}{c}{B}_{4X}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\ln [(z+L/2{)}^2+(x+r_1/2{)}^2]-\frac{{\mu }_0{k}_v}{4\pi }\ln [(z-L/2{)}^2+(x+r_1/2{)}^2]\hfill \\ B_{4Z}(x,z)=\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z+L/2}{x+r_1/2}\right)-\frac{{\mu }_0{k}_v}{4\pi }\arctan \left(\frac{z-L/2}{x+r_1/2}\right)\hfill \end{array}$$

(19)

where $B_{2X}$ , $B_{3X}$ , and $B_{4X}$ are X-direction component of point P(x, y, z) magnetic induction intensity created by surface current 2, 3, and 4, respectively; $B_{2Z}$ , $B_{3Z}$ , and $B_{4Z}$ are Z-direction component of point P(x, y, z) magnetic induction intensity created by surface current 2, 3, and 4, respectively

$$B=\sqrt{{(B_{1X}+B_{2X}+B_{3X}+B_{4X})}^2+{(B_{1Z}+B_{2Z}+B_{3Z}+B_{4Z})}^2}$$

(20)

Critical shear yield stress ${\tau }_y$ can be calculated by equation (19) according to selected magnetorheological liquid feature

$${\tau }_y=\alpha \prime B^{k\prime }$$

(21)

where $\alpha \prime$ and $k\prime$ are factors of selected magnetorheological liquid.

Damping force model of normal MRD can be achieved as follows according to Bingham model ²⁵

$$\begin{array}{c}\multicolumn{1}{c}{F_d=(\frac{24\eta {A_p}^{2}l}{b{h}^3}+\frac{2\eta \mathrm{bl}}{h})v}\\ +(\frac{4l{A}_p}{\eta }+2\mathrm{bl}){\tau }_y·\mathrm{sgn}(v)+f\hfill \end{array}$$

(22)

where $F_d$ is the damping force of MRD ; $b$ is the width of flat plate; $h$ is the height of damping channel ; $l$ is the length of damping channel; $\eta$ is the apparent viscosity of magnetorheological liquid; and $f$ is the offset force.

Equation (23) can be achieved through transforming equation (22)

$$F_d=(\frac{24\eta {Q_2}^{2}l}{\mathrm{vb}{h}^3}+\frac{2\eta \mathrm{blv}}{h})+(\frac{4l{Q}_2}{\eta v}+2\mathrm{bl}){\tau }_y·\mathrm{sgn}(v)+f$$

(23)

Damping force F_d can be achieved through putting the parameters of designed PMMVMRD into equation (23)

$$\begin{array}{c}{F}_d=(\frac{24\eta {(Q-Q_1)}^2(L+2L_1)}{\pi (r_3+r_2)v{(r_3-r_2)}^3}+\frac{2\eta \pi (r_3+r_2)(L+2L_1)v}{(r_3-r_2)})\hfill \\ +(\frac{4(L+2L_1)(Q-Q_1)}{\eta v}+2\pi (r_3+r_2)(L+2L_1)){\tau }_y·\mathrm{sgn}(v)+f\hfill \end{array}$$

(24)

Experiments of PMMVMRD

In order to verify working performance of the studied PMMVMRD, a damper is manufactured on the basis of an automotive suspension technique requirement of Chery Automobile Co, Ltd., and the commercial magnetorheological liquid provided by Ningbo Shangon Center of Structrual Monitoring and Control Engineering Co, Ltd. The experiments on PMMVMRD are carried out through PA-20-Z 20 kN electrohydraulic servo suspension system produced by Shanghai Beiyue Machine Test CO, Ltd. In the experiments, on the basis of QCT 545-1999 (bench experiment methods of automotive telescopic shock absorbers) and by taking road excitation of coming real vehicle tests into consideration, sine wave is employed to be input signal, the maximum exciting velocity is 0.3 m/s, amplitude is 30 mm, frequency is 1.95 Hz, and input electric current is 0–1.2 A. The test rig is shown in Figure 4, and the experimental results are shown in Figures 5 and 6, in which, Figure 5(a)–(d) shows damping force–displacement curves and Figure 6 shows damping force–velocity curves when the input electric current is 0, 0.4, 0.8, and 1.2 A, respectively. Figure 5 indicates when input current in the coil is 0 A, the damping force reaches maximum, and the maximum value can be 1 kN. The reason lies in that: when input current is 0 A, the magnetic valve shuts off, and all the magnetorheological liquid flows through damping channel so that large damping force is created. When input current in the coil is 1.2 A, the damping force reaches minimum, nearly equal to 0. The reason lies in that: when the input current in the coil is 1.2 A, the magnetic valve is opened fully, the major magnetorheological liquid flows through magnetic valve and few liquid flows through damping channel so that small damping force is created. With the increase of input current in the coil, the open degree of the magnetic valve increases accordingly and so does flowing rate. Thus, liquid flowing through damping channel decreases, and the damping force decreases accordingly. The conclusion can be drawn from the above experiment analysis that damping force of the designed MRD can be adjusted through changing the input current in the magnetic coil so that PMMVMRD can meet technique requirement of automotive suspensions. Moreover, when the input electric current is 0, the PMMVMRD can output bigger damping force so that the proposed PMMVMRD can solve the fail-safe problem of normal dampers. Meanwhile, because permanent magnets are introduced into PMMVMRD, magnetorheological fluid sedimentation problem of normal MRDs can be solved with magnetic field of permanent magnets. Therefore, the proposed PMMVMRD has better working performance than normal servo-valve dampers and MRDs. It can be used not only in automobiles but also in vibration damping fields with high requirements on dampers, such as buildings and mechanic systems.

Figure 4.

Test rig of PMMVMRD.

Figure 5.

Damping force varying with input current at 30 mm amplitude and 1.95 Hz frequency: (a) damping force versus displacement at 0 A current; (b) damping force versus displacement at 0.4 A current; (c) damping force versus displacement at 0.8 A current; and (d) damping force versus displacement at 1.2 A current.

Figure 6.

Experimental damping force versus velocity.

Figure 7 is damping force–velocity curve simulated by the designed model in this article. For the convenience of analyzing and comparing simulation results of damping force–velocity and tested results of actual MRD, characteristic physical dimensions of PMMVMRD models are similar to those of actual tested MRD. It can be indicated that simulated curves agree with those of experiments under different current and velocity. The comparison of Figures 6 and 7 shows the error between simulation and experimental results is very small. In compression travel, maximum damping force of simulations and actual experiments are 1084 and 1007 N, respectively, with the error being less than7.6%. In rebounding travel, maximum damping force of simulations and actual experiments are 1479 and 1405 N, respectively with the error being less than 5.2%. The above analysis verifies that the damping force model designed in this article is correct.

Figure 7.

Simulated damping force versus velocity.

Design of GSO-PID

Development of GSO

PID controller has three controlling units: proportion, integration, and differentiation. The difference between desired controlling values and actual controlling values is the controller input. Design of PID controller is realized through determining proportionality factor, integration time factor, and differentiating time factor. PID controllers have many advantages, such as rapid responding speed, high controlling accuracy, and good robustness, which can enhance working performance of controlled systems greatly. However, the factors of PID controllers are difficult to determine, so the application of PID controllers is limited. GSO algorithm with many advantages such as greater searching ability and higher optimizing accuracy is a new intelligent optimization algorithm designed by Indian scholars Krishnanand and Ghose in 2005. Therefore, GSO algorithm has been used in many fields, such as assembly sequence planning, mobile robot localization, traveling salesman problems, and wavelength variable selection method in near infrared spectroscopy. It is proved that the application effects are good. If PID controllers and GSO algorithm can be combined to design GSO-PID controlling algorithm, the new controlling algorithm hopes to enhance the working performance of automotive electronic systems greatly. Thus, in order to enhance controlling performance of PMMVMRD air suspension, a novel GSO-PID algorithm is designed and employed in this article to carry out control experiments of suspension system.

The basic working theory is as follows. At the beginning, n glowworms are distributed randomly into a D-dimensional space, and every glowworm has its own fluorescein determined by its position. The better the position a glowworm has, the more the fluorescein it has and the more intensive the light it can radiate. Hence, the optimized objective value is better. Each glowworm has its own decision domain, in which glowworms select others with more fluorescein than themselves to construct neighborhoods. In the neighborhoods, each glowworm selects randomly another glowworm with higher fluorescein and move one step toward the selected one in each iteration. In the end, glowworms will stay around several positions with good objective values through continuous moving so that searching on the whole space can be completed.

Literature ²⁶ shows that many scholars across the world have carried out deep research on GSO algorithm. It has been widely used in controlling collective robots, solving 0-1 Knapsack Problem, optimizing multi-modal functions, and detecting multi signal sources, with desirable applying effects. ²⁷ However, nowadays, GSO algorithm still has some disadvantages, such as low optimization velocity, being easy to fall into local optimization and the like. ²⁸ Thus, in view of the above problems, optional moving operation, moving variation operation, and operation of step size changing with fluorescein proportionality are introduced into normal GSO algorithm. Thus, a novel GSO algorithm is developed and its specific theory is as follows:

1. Optional moving operation. In order to avoid slow optimizing velocity, caused by blind selection of glowworms with more fluorescein to move in optimization, the following operations are conducted: order all the glowworms in neighbors by their fluorescein, select k glowworms with more fluorescein; then, select and move to a glowworm in the k glowworms randomly.

2. Moving variation operation. To avoid falling in local optimum, operation of moving with certain-probability variation is introduced in glowworm movement; that is, in the sixth step, select and move to jth glowworm in its neighbor and renew position according to equation (25)

$$\begin{array}{c}{x}_i(t+1)=x_i(t)+s_t\left(\frac{x_j(t)-x_i(t)}{\Vert x_j(t)-x_i(t)\Vert }\right)\cos \gamma \prime \hfill \\ (\gamma \prime =0^{\circ }\mathrm{or}{180}^{\circ })\hfill \end{array}$$

(25)

where x_i and x_j are ith and jth glowworm position, respectively; s_t is step size; $\cos \gamma \prime$ is the moving factor. When ith glowworm moves to jth glowworm, ith one is permitted to move in opposite direction of jth one with specified probability. $\gamma \prime$ in equation (25) is determined by equation (26)

$$\{\begin{array}{cc}0\le \mathrm{ran}\mathrm{dom}(100)<98\hfill & \gamma \prime =0^{\circ }\hfill \\ \multicolumn{1}{c}{98\le \mathrm{ran}\mathrm{dom}(100)<100}& \gamma \prime ={180}^{\circ }\hfill \end{array}$$

(26)

• 3. Operation of step size changing with fluorescein proportionality. Step size s_t is a constant value in the optimization of normal GSO algorithm, so good objective domain will be missed if the step size is too big, and slow velocity will be caused if the step size is too small. Therefore, in the optimization of developed GSO algorithm in this article, the operation, in which step size changes with ratio between the moving glowworm fluorescein and its selected glowworm fluorescein, is employed. More specifically, big step size is employed for big ratio, while small step size for small ratio. So s_t can be described by equation (27)

$$s_t=s_1(1-\frac{l_i}{l_j})$$

(27)

where s₁ is step size parameter; l_i and l_j are fluorescein of ith and jth glowworm.

In order to verify the optimizing effect of GSO algorithm with the proposed three operations, simulation experiments are conducted by using six standard testing functions $(f_1~f_6)$ , which are from literature. ²⁹ GSO algorithm is employed to search the minimum values of the six functions, with their names, features, and number n of glowworm being displayed in Table 1. In the simulation, the parameters of GSO algorithm are: neighbor threshold n_t = 5, original fluorescein l₀ = 5, changing ratio of neighbor β = 0.08, volatilizing factor of fluorescein ρ = 0.5, and extracting ratio of fitness γ = 0.7. In order to avoid accidental error, experiments on normal GSO algorithm and developed GSO algorithm are conducted for 20 times. The results are shown in Figure 8(a)–(f), which indicate that convergence velocity of the developed GSO algorithm is faster than that of normal GSO algorithm. Similarly, they also indicate the developed GSO algorithm has higher optimization accuracy. Thus, the simulation has proved that the optimization ability of the developed GSO algorithm is stronger than that of normal one.

Table 1.

Tested functions and their minimum values.

	Function names	Searching scope	Minimum value	n
f1	20 + (x2–10cos(2πx) + y2–10cos(2πy))	[−5.12,5.12]	0	50
f2	x2 + y2	[−10,10]	0	100
f3	$0.5+((\sin (\sqrt{x_1^2+x_2^2}){)}^2-0.5)/(1+0.001(x_1^2+x_2^2){)}^2$	[−10,10]	0	100
f4	$\sum _{i=1}^n{\mathrm{ix}}_i^2$	[−5.12,5.12]	0	50
f5	$\sum _{i=1}^n(0.2x_i^2+0.1x_i^2\sin 2x_i)$	[−10,10]	0	100
f6	$\sum _{i=1}^n{(x_i+0.5)}^2$	[−10, 10]	0	100

Figure 8.

Algorithm comparison of optimizing f₁–f₆: (a) algorithm comparison of optimizing tested function f₁; (b) algorithm comparison of optimizing tested function f₂; (c) algorithm comparison of optimizing tested function f₃; (d) algorithm comparison of optimizing tested function f₄; (e) algorithm comparison of optimizing tested function f₅, and (f) algorithm comparison of optimizing tested function f₆.

Design of GSO-PID

In sky-hook algorithm, a damper is supposed to be mounted between inertial coordinate system and sprung weight, and its force direction is opposite to motion direction of the sprung weight. Sky-hook damper can control absolute moving velocity of sprung weight so as to enhance automotive riding comfort. At present, sky-hook algorithm is widely used in automotive suspensions. Thus, PMMVMRD is controlled by sky-hook algorithm, while air spring of semi-active suspension is controlled by GSO-PID algorithm in this article. Before determining the sky-hook factors, reasonable stiffness coefficient and dampers’ damping factors of front and rear suspensions are given by designers according to researched automobile system features. In order to enhance handling and riding performance of automobiles, transfer function H_T of tire dynamic load and transfer function H_S of automotive riding performance are deduced based on suspensions system models. A group of amplitude–frequency characteristic curves can be achieved through changing the value of sky-hook factors and simulations. Hence, the value of satisfying both handling and riding performance of automobiles is selected as the sky-hook factor. Thus, in the process of optimizing sky-hook factors, stiffness coefficients of suspension springs are regarded as constant value, and PID controller does not work in this process. Based on transfer functions of PMMVMRD air suspension systems, simulations are conducted by using different factors of sky-hook algorithm. Finally, the better sky-hook factors of front suspension and rear suspension are determined as 2100 Ns/m and 2300 Ns/m, respectively. In order to enhance robustness of air spring, GSO algorithm is employed to optimize three factors of PID algorithm, namely, proportionality factor k_p, integration time factor k_i, and differentiating time factor k_d. The working principle is as follows. First, every glowworm in a group is distributed in a three-dimensional (3D) space $[{\alpha }_1,{\alpha }_2,{\alpha }_3]$ , determined by the possible value domain of three factors of PID controller, that is, $k_p\in {\alpha }_1,k_i\in {\alpha }_2,k_d\in {\alpha }_3$ . Thus, every glowworm corresponds to a specific position in 3D space, representing specific values of PID factors. Meanwhile, every glowworm radiates fluorescence light determined by objective value of its specific position. For the studied GSO-PID algorithm, automotive dynamic performance indexes of applying PID controller (corresponding to glowworm position) in controlling air suspensions is employed to express objective value of glowworm position. Performance indexes will be explicitly defined by the following equation (28). The glowworm fluorescein can be calculated through equation (28) and will be explicitly defined by the following equation (29).

Meanwhile, the more shining a glowworm is, the better the position it has, which means that the glowworm has better objective value, and corresponding PID controller has better factors. In the space, every glowworm has its own visual domain, in which, every glowworm moves to the more shining one with a specified probability. Every time, the glowworm moves a step size, that is, corresponding PID controller of glowworm position moves to a better one for a step size. Finally, most of the glowworms will gather around several ones with higher fitness through the continuous movement and optimized values, that is, optimized PID controller, are found out.

In order to speed up optimization velocity of GSO algorithm and avoid local optimum, the three operations, such as optional moving operation, moving variation operation, and operation of step size changing with fluorescein proportionality, which are proposed and introduced in GSO algorithm in previous section (sections “Development of GSO” and “Design of GSO-PID”). On the basis of PMMVMRD air suspension models, the performance index J of PID controller is expressed by equation (28)

$$J={\beta }_1{J}_1+{\beta }_2{J}_2+{\beta }_2\prime J_2\prime +{\beta }_3{J}_3+{\beta }_3\prime J_3\prime$$

(28)

where J₁, J₂, J₂′, J₃, and J₃′ are root mean square (RMS) values of acceleration of automobile mass center, front-suspension working space, rear-suspension working space, front-tire dynamic load, and rear-tire dynamic load, respectively; β₁, β₂, β₂′, β₃, and β₃′ are weighted values.

The control scheme of GSO-PID is shown in Figure 9, and the working flow of the developed GSO algorithm is as follows:

• Step 1: Initiate parameters.

• Step 2: Distribute ith glowworm into 3D space $[{\alpha }_1,{\alpha }_2,{\alpha }_3]$ .

• Step 3: Transfer corresponding objective function J(x_i(t)), of x_i(t) position of ith glowworm in tth iteration, into corresponding fluorescein l_i(t)

$$l_i(t)=(1-\rho )l_i(t-1)+\sigma J(x_i(t))$$

(29)

where $\rho$ is the volatilizing factor of fluorescein, and $\sigma$ is the extracting factor of fitness ratio.

• Step 4: Every glowworm selects glowworms with more fluorescein in domain within dynamic decision radius to make up its neighbor N_i(t)

$$N_i(t)=\{j:\Vert x_j(t)-x_i(t)\Vert <r_d^i(t);l_i(t)<l_j(t)\}$$

(30)

where $r_d^i(t)$ is the dynamic decision radius of glowworm; $r_s$ is the perception radius of glowworm; $r_d^i(t)$ and $r_s$ meet the relation 0 < r_dⁱ(t) < r_s.

• Step 5: According to “optional moving operation” proposed in sections “Development of GSO” and “Design of GSO-PID,” order all the glowworms in neighbors and select k glowworms with more fluorescein; then, calculate the probability p_ij(t) of ith glowworm moving to jth glowworm in its k glowworms through equation (31)

$$p_{\mathrm{ij}}(t)=\frac{l_j(t)-l_i(t)}{\sum _{k\in N_i(t)}{l}_k(t)-l_i(t)}$$

(31)

• Step 6: Select jth glowworm through roulette and then move to it. Position is renewed according to equation (25). Thus, “moving variation operation” and “operation of step size changing with fluorescein proportionality” proposed in section “Development of GSO” and “Design of GSO-PID” are achieved through moving factor $\cos \gamma$ and the changing step size s_t.

• Step 7: Renew radius of dynamic decision domain through equation (32)

$$r_d^i(t+1)=\min \{r_s,\max \{0,r_d^i(t)+\beta (n_t-|N_i(t)|)\}\}$$

(32)

where β is the variation ratio of neighbor, and n_t is the threshold of neighbor.

• Step 8: If GSO algorithm reaches maximum iteration, output three factors k_p, k_i, k_d of PID controller represented by optimum position of the glowworm group. Otherwise, calculate dynamic decision domain $r_d^i(t+1)$ of next time for every glowworm; let t = t + 1 and then go to step 3.

Figure 9.

Control scheme of GSO-PID.

First, in the process of optimization simulation, specific numerical value of objective performance index is determined by the designers according to the features of controlled suspension systems and controlling requirement. Then, real-time signals, such as acceleration of automobile mass center, front-suspension working space, rear-suspension working space, front-tire dynamic load, and rear-tire dynamic load, are collected so that the RMS values of these real-time signals can be calculated. On this basis, actual performance index is calculated through equation (28). Actual performance index is the feedback of PID controllers of suspension systems, and difference between actual performance index and objective performance index is the deviation of the PID controller.

Modeling of PMMVMRD air suspension

Modeling of airbag

The structure diagram is shown in Figure 10. ³⁰ Suppose that the pressure in the air gasbag is p and pressure atmosphere is p_a, then output force F_a of air spring can be expressed by equation (33)

$$F_{\mathrm{a}}=(p-p_a)A_e$$

(33)

where A_e is effective area of airbag.

Figure 10.

Structure diagram of airbag.

According to Hooke’s law, stiffness k_b of air spring can be calculated by equation (34)

$$k_b=\frac{d{F}_a}{\mathrm{ds}}=\frac{\mathrm{dp}}{\mathrm{ds}}{A}_e+(p-p_a)\frac{d{A}_e}{\mathrm{ds}}$$

(34)

where s is the vertical displacement of air spring; in this article, s is the difference between automotive sprung mass and unsprung mass.

Air in the bag meets ideal gas state equation (35)

$$p{V}^n=k_0$$

(35)

where k_o is constant value; V is volume of airbag, n is air polytropic exponent.

Equation (36) can be achieved through working out derivation on both sides of equation (35) about s

$$\frac{\mathrm{dp}}{\mathrm{ds}}{V}^n+\mathrm{pn}{V}^{n-1}\frac{\mathrm{dv}}{\mathrm{ds}}=0$$

(36)

Equation (37) can be achieved through transforming equation (36)

$$\frac{\mathrm{dp}}{\mathrm{ds}}=-\frac{\mathrm{pn}}{V}\frac{\mathrm{dV}}{\mathrm{ds}}$$

(37)

Equation (38) can be achieved according to airbag feature

$$\frac{\mathrm{dV}}{\mathrm{ds}}=A_e$$

(38)

Equation (39) can be achieved through putting equations (37) and (38) into equation (34)

$$k_b=-\frac{\mathrm{pn}}{V}{A_e}^2+(p-p_a)\frac{d{A}_e}{\mathrm{ds}}$$

(39)

Because the airbag is used by air suspension, equation (40) can be achieved

$$\frac{d{A}_e}{\mathrm{ds}}\approx 0$$

(40)

Therefore, stiffness of air spring can be calculated through equation (41)

$$k_b=-\frac{\mathrm{pn}}{V}{A_e}^2$$

(41)

Thus, elastic force of air spring can be achieved through equation (41)

$$F_{\mathrm{a}}=k_{b}s$$

(42)

Equation (43) can be achieved through equation (42)

$$F_{\mathrm{a}}=-\frac{\mathrm{pn}}{V}{A_e}^{2}s$$

(43)

It can be concluded that stiffness and elastic force of air spring can be adjusted through controlling the pressure in the airbag.

According to the GSO-PID control scheme shown in Figure 9, the electronic control unit (ECU) of PMMVMRD air suspension collects the signals such as automotive mass center acceleration, front-suspension working space, rear-suspension working space, front-tire dynamic load, and rear-tire dynamic load. The actual performance index is calculated based on the collected signals and compared with the target index in ECU so that pressure in the airbags is calculated by GSO-PID controllers. Inflating or outgassing for airbags are determined through comparison between the calculated and actual pressure in airbags. Finally, pressure in airbags is adjusted through controlling magnetic valves so that working performance of PMMVMRD air suspension is optimized.

Modeling of 4 degree-of-freedom suspension

The automobile with 4 degree-of-freedom is researched, and sky-hook algorithm is employed to control the PMMVMRD air suspension. First, suppose that there is a sky-hook damper on each automobile body corresponding to front and rear axis. The physical model of the PMMVMRD air suspension is shown in Figure 11. Based on the physical model, mathematical models are built.

Figure 11.

Physical model of 4 degree-of-freedom air suspension.

Kinematic equation of automotive body can be expressed as equation (44)

$$\{\begin{array}{c}{m}_s{x}_c″=F_{a1}+F_{a2}+C_1(x_{s1}\prime -x_{t1}\prime )+C_2(x_{s2}\prime -x_{t2}\prime )+C_{\mathrm{skf}1}{x}_{s1}\prime +C_{\mathrm{skf}2}{x}_{s2}\prime \hfill \\ \multicolumn{1}{c}{J_{\varphi }\varphi ″=(F_{a2}+C_2(x_{s2}\prime -x_{t2}\prime )+C_{\mathrm{skf}2}{x}_{s2}\prime )(L_x-b_x)-(F_{a1}+C_1(x_{s1}\prime -x_{t1}\prime )+C_{\mathrm{skf}1}{x}_{s1}\prime )b_x}\end{array}$$

(44)

where m_s is the automotive sprung mass; x_c is the automotive center displacement; L_x is the automotive axle base; b_x is the distance from automotive mass center to front axis; Φ is the pitching angular of automobile; J_Φ is the rotary inertia; F_a1, F_a2 are the elastic force of air springs of front and rear suspensions; C ₁ , C₂ are the dampers’ damping factors of front and rear suspensions; x_t1, x_t2 are the unsprung weight displacement of front and rear suspensions; x_s1, x_s2 are the sprung weight displacement of front and rear suspensions; k_t1, k_t2 are the stiffness coefficient of front and rear tires; and x_r1, x_r2 are the disturbance input from the road at front and rear tires.

Kinematic equation of unsprung weight is expressed as equation (45)

$$\{\begin{array}{c}{m}_1{x}_{t1}″=-F_{a1}-C_1(x_{s1}\prime -x_{t1}\prime )-k_{t1}(x_{t1}-x_{r1})\hfill \\ \multicolumn{1}{c}{m_2{x}_{t2}″=-F_{a2}-C_2(x_{s2}\prime -x_{t2}\prime )-k_{t2}(x_{t2}-x_{r2})}\end{array}$$

(45)

where m_t1, m_t2 are unsprung weight of front and rear suspensions.

Because sky-hook dampers are imaginary, the damping force of sky-hook force is generated by PMMVMRDs. Therefore, damping force of PMMVMRDs of front and rear air suspensions is the resultant force between sky-hook damping force and damping force generated by original damping factors, that is

$$\{\begin{array}{c}{F}_{d1}=C_1(x_{s1}\prime -x_{t1}\prime )+C_{\mathrm{skf}1}{x}_{s1}\prime \hfill \\ \multicolumn{1}{c}{F_{d2}=C_2(x_{s2}\prime -x_{t2}\prime )+C_{\mathrm{skf}2}{x}_{s2}\prime }\end{array}$$

(46)

where F_d1, F_d2 are the actual output force of PMMVMRDs of front and rear suspensions.

In the actual control of PMMVMRD air suspensions, signals such as x_t1, x_t2, x_s1, and x_s2 are collected by the controlling systems, and actual force of PMMVMRDs is calculated according to equation (46). Then, PMMVMRD controlling electric current of generating F_d1 and F_d2 is calculated through equation (24). Hence, the proper electric current is applied to PMMVMRDs through the controlling systems. Proper damping force is generated so that performance of PMMVMRD air suspensions can be optimized.

Simulation experiments

In order to verify the working performance of studied PMMVMRD air suspension controlled by GSO-PID algorithm, a testing automobile from Chery Automobile Co, Ltd is employed to conduct the simulation experiments, and main parameters of the testing automobile are: m_s = 730 kg, J_ϕ = 1230 kg m², m_t1 = 40 kg, m_t2 = 36 kg, k_t1 = 175 KN/m, k_t2 = 175 KN/m, b_x = 1 m, and L_x = 2.8 m. By referring to the method in literature ³¹ and in view of automotive safety, riding comfort, and handling performance comprehensively, weighted values are set as follows: β₁ = 1, β₂ = 2865.7, β₂′ = 2865.7, β₃ = 355.2, and β₃ = 413.6. The simulation experiments are conducted at three working conditions, including input of impact road surface, sine wave road surface, and white noise road surface.

Simulation experiments of impact road surface

Simulation experiments of impact road surface are conducted according to Chinese Standard GB5902-86. In this simulation, the automobile passes a single bump, and the vibration of the automobile is employed to evaluate the experiments. Suppose that the driving speed of automobile is u′, then time domain expression of road pulse impact on tires can be expressed by equation (47)

$$x_r(t)=\frac{A_{\mathrm{m}}}{2}(1-\cos \frac{2\pi u\prime }{L\prime }t)0\le t\le \frac{L\prime }{u\prime }$$

(47)

where u′ is A_m = 0.1 m and L′ = 5 m.

In order to investigate working performance of PMMVMRD air suspension controlled by GSO-PID algorithm, especially the effect of applying PMMVMRD in air suspension, experimental results are compared with those of normal MRD suspension controlled by GSO-PID and passive suspension. During the GSO-PID optimization, impact road surface is employed as road-induced excitation. The normal MRD is from the testing automobile, and the main parameters are as follows: inner diameter of working cylinder is 30 mm; length of damping channel is 20 mm; radius of piston is 13 mm; working gap is 1.0 mm; diameter of piston rod is 10 mm; turn number of coil is 100.

The results are shown in Figure 12(a)–(d), which indicate that peak value of mass center acceleration, body pitching angle acceleration, suspension working space and tire dynamic load of PMMVMRD air suspension controlled by GSO-PID, normal MRD air suspension controlled by GSO-PID, and passive suspension increases in sequence. Meanwhile, adjusting time of the above four values increases in sequence, too. Therefore, better control results and shorter adjusting time can be achieved by the proposed PMMVMRD air suspension. These demonstrate that both of the semi-active suspension can enhance suspension performance. PMMVMRD air suspension controlled by GSO-PID has better working performance than normal MRD air suspension controlled by GSO-PID, which indicates that the effect of applying PMMVMRD in air suspension is better than that of applying normal MRD in air suspension. In order to conduct further comparison between the two semi-active suspensions, peak values of mass center acceleration, pitching angle acceleration, suspension working space, and tire dynamic load are calculated and shown in Table 2. Table 2 indicates that compared with normal MRD air suspension, mass center acceleration, pitching angle acceleration, suspension working space, and tire dynamic load of PMMVMRD suspension controlled by GSO-PID decreases by 20%, 10.9%, 33.3%, and 11.8%, respectively. This demonstrates that riding and handling performance of PMMVMRD air suspension controlled by GSO-PID are better than those of normal MRD air suspension controlled by GSO-PID. Therefore, both quality and quantity analysis show that in the simulation experiment, PMMVMRD air suspension controlled by GSO-PID has good performance and effect of applying PMMVMRD in air suspension is good.

Figure 12.

Time-domain simulation results of impact road surface: (a) mass center acceleration; (b) body pitching angle acceleration; (c) suspension working space; and (d) tire dynamic load.

Table 2.

Simulation result comparison of different suspensions.

		Passive suspension	Normal MRD air suspension	PMMVMRD air suspension
Mass center acceleration	Maximum amplitude/(m s−2)	4.821	4.655	3.726
Pitching angle	Maximum amplitude/(rad s−2)	0.477	0.421	0.375
Suspension working space	Maximum amplitude/m	0.0146	0.012	0.008
Tire dynamic load	Maximum amplitude/mm	1257	1165	1027

Simulation experiments of sine wave road surface

In order to investigate the controlling effect of GSO-PID algorithm used in PMMVMRD air suspension, sine wave road surface is employed to conduct simulation experiments, in which driving speed of automobile is 40 km/h. The expression of sine wave road surface is expressed as equation (45). In order to compare GSO-PID with other control algorithms, particle swarm optimization (PSO) algorithm and genetic algorithm (GA) are employed to optimize PID controller so that PSO-PID and GA-PID controlling algorithms are designed. When PSO algorithm is employed to optimize PID controller, the main parameters are set as follows: population scale, maximum iteration, and learning factor are 50, 400, and 2, respectively. When GA algorithm is employed to optimize PID controller, main parameters are set as follows: population scale, maximum iteration, cross probability, and variation probability are 100, 50, 0.8, and 0.1, respectively. Simulation results are compared with those of PMMVMRD air suspension controlled by PSO-PID, GA-PID, and passive suspension and the results are shown in Figure 13(a)–(f)

$$z_0=0.05\times \sin (4\pi t)$$

(48)

Figure 13(a) indicates that vibration acceleration amplitude of passive suspension, PMMVMRD air suspension controlled by GA-PID, PSO-PID, and GSO-PID decreases in sequence, which shows that under the control of above algorithms, riding comfort increases in sequence. Among them, PMMVMRD air suspension controlled by GSO-PID has the best riding comfort. Figure 13(b) indicates that pitching angle acceleration amplitude of suspension controlled by three algorithms is smaller than that of passive suspension, which means that the three algorithms can control automotive body posture effectively. Among them, suspension controlled by GSO-PID has the smallest pitching angle acceleration and the best performance. Figure 13(c) and (d) indicates that working space of passive suspension and PMMVMRD air suspension controlled by GA-PID, PSO-PID, and GSO-PID decreases in sequence, which means that the probability of knocking at stop block decreases in sequence, too and PMMVMRD air suspension controlled by GSO-PID has the best performance. Figure 13(e) and (f) indicates that under the control of three algorithms, tire dynamic load is smaller than that of passive suspension, which means that three algorithms can suppress tire jumping, and GSO-PID has the best controlling effect.

Figure 13.

Time-domain simulation results of sine wave road surface: (a) mass center acceleration; (b) body pitching angle acceleration; (c) front-suspension working space; (d) rear-suspension working space; (e) front-tire dynamic load; (f) rear-tire dynamic load.

In order to verify the effect of applying GSO-PID algorithm in PMMVMRD air suspension, RMS values of six performance indexes are calculated and shown in Table 3. Table 3 indicates that for PMMVMRD air suspension, controlling effect of PSO-PID is better than that of GA-PID algorithm. It can also indicate that by GSO-PID algorithm, RMS values of mass center acceleration, pitching angle acceleration, front-suspension working space, rear-suspension working space, front-tire dynamic load, rear-tire dynamic load of PMMVMRD air suspension decrease by 23.2%, 28.3%, 35.1%, 25.7%, 19.9%, and 18.1%, respectively than those of PMMVMRD air suspension controlled by PSO-PID algorithm. Hence, a conclusion can be reached that controlling effect of GSO-PID is better than that of PSO-PID with riding comfort and handling performance being enhanced further.

Table 3.

Simulation result comparison of different control algorithms.

Suspension performance index	Passive suspension	GA-PID	PSO-PID	GSO-PID
Mass center acceleration RMS/m s−2	5.117	4.814	3.319	2.547
Pitching angle acceleration RMS/rad s−2	1.294	0.717	0.674	0.653
Front-suspension working space RMS/m	0.054	0.047	0.037	0.024
Rear-suspension working space RMS/m	0.065	0.044	0.035	0.026
Front-tire dynamic load/N	1164	957	846	677
Rear-tire dynamic load/N	1095	935	875	716

Simulation experiments of white noise road surface

In most cases, automobile speed is set as a constant value in suspension simulations. However, the actual speed is changing, which has a great impact on suspension performance. In order to fully investigate the performance of novel PMMVMRD air suspension, road surface of white noise is employed to conduct the simulation, and the results are compared with those of automobile with normal MRD suspensions at different speed. The road surface is taken as the simulation input, and the irregularity coefficient is $G_{z0}(n_0)=64\times {10}^{-6}{\mathrm{m}}^3$ . In the simulation, automotive speed is set as 30, 50, and 80 km/h, and the results are shown in Figure 14(a)–(c). Figure 14(a)–(c) shows the simulation results of mass center acceleration of automobile with normal MRD suspensions, simulation results of mass center acceleration of automobile with PMMVMRD suspensions, and simulation result comparison between automobiles with PMMVMRD and normal MRD suspensions. Through the comparison between Figure 14(a) and (b), it can be indicated: under the working condition at different speed, the automobile with PMMVMRD air suspensions controlled by GSO-PID have better vibration damping performance than the one with normal MRD suspensions. When the automotive speed is 30, 50, and 80 km/h, mass-center-acceleration peak value of automobile with normal MRDs suspension is 5.411, 4.375, and 4.024 m s⁻², respectively. Meanwhile, mass center–acceleration peak value of automobile with PMMVMRD suspension is 4.625, 3.522, and 3.179 m s⁻², respectively. Therefore, under the working condition at three different speed, mass center–acceleration peak value of automobile with PMMVMRD suspension decreases by 10%, 11%, and 13%, respectively, compared with those of automobile with MRD suspensions. Figure 14(c) is about the comparison between mass center–acceleration RMS value of automobiles with PMMVMRD air suspension and normal MRD suspensions at different speed. It indicates that riding comfort of the automobile with PMMVMRD air suspensions increases by 17%, 13%, and 12%, respectively at the speed of 30, 50, and 80 km/h.

Figure 14.

Time-domain simulation results of white noise road surface: (a) simulation results of mass center acceleration of automobile with normal MRD suspensions; (b) simulation results of mass center acceleration of automobile with PMMVMRD suspensions; and (c) simulation result comparison between automobiles with PMMVMRD and normal MRD suspensions.

Conclusion

In order to enhance working performance of air suspension, a PMMVMRD and GSO-PID controlling algorithm are designed. On the basis, three simulation experiments are concluded, and the work in this article is summarized as follows.

1. To enhance comprehensive performance of air suspension, a novel MRD with a permanent magnet and a magnet valve is designed. First, based on equivalent surface current model, magnetic field model in damping channel is built and then mathematical model of PMMVMRD is built by Bingham model and theories about magnetic valve. Experiments on the novel damper are conducted, and experimental results agree with those simulated by PMMVMRD model. On the basis of newly designed damper model, 4 degree-of-freedom automotive suspension model is built.

2. In order to enhance optimization velocity of GSO algorithm and avoid local optimum, normal GSO algorithm is developed. On this basis, GSO-PID suitable for PMMVMRD air suspension is designed. The optimized PID controllers by GSO have better factors so that the controlling systems by GSO-PID algorithm have higher controlling accuracy and responding speed. Thus, GSO-PID controlling algorithm can control the pressure in airbags with higher accuracy so that the PMMVMRD air suspensions have better working performance.

3. Simulation experiments on PMMVMRD air suspension controlled by GSO-PID algorithm are conducted with impact road surface, sine wave road surface and white noise road surface being the input. The experimental results indicate that PMMVMRD air suspension controlling GSO-PID algorithm has good working performance.