# Stability, Bifurcation, and Quenching Chaos of a Vehicle Suspension System

Authors: Chun-Cheng Chen; Shun-Chang Chang
Journal: Journal of Advanced Transportation
 (2018)

---

## Abstract

This study investigated the dynamics and control of a nonlinear suspension system using a quarter-car model that is forced by the road profile. Bifurcation analysis used to characterize nonlinear dynamic behavior revealed codimension-two bifurcation and homoclinic orbits. The nonlinear dynamics were determined using bifurcation diagrams, phase portraits, Poincaré maps, frequency spectra, and Lyapunov exponents. The Lyapunov exponent was used to identify the onset of chaotic motion. Finally, state feedback control was used to prevent chaotic motion. The effectiveness of the proposed control method was determined via numerical simulations.


---

## Body

## 1. Introduction

The suspension system of vehicles is meant to isolate the vehicle body and passengers from oscillations imposed by road irregularities and maintain continuous contact between the tires and surface of the road. Suspension systems provide a smoother ride and help the operator to maintain control of the vehicle over rough terrain or during sudden stops [1]. Hysteresis is an important aspect of suspension design, and the development of smart materials for vehicle suspension systems with hysteretic damping has been the subject of considerable studies [2–5]. The fact that the damping force in a vehicle’s suspension system is hysteretic means that road excitation can be modeled as a multifrequency harmonic function. Therefore, automotive suspension systems may be considered a nonlinear hysteretic system with multifrequency excitation.Nonlinear phenomena of codimension bifurcations have been studied extensively, particularly in engineering systems [6–8]. Tseng and Tung [9] examined the bifurcation behavior of a nonlinear structure with a magnetic actuator. Yagasaki [10] investigated nonlinear dynamics and the bifurcation behavior of a pendulum with feedback control. Tian et al. [11] studied the complex nonlinear dynamics of smooth, discontinuous oscillators. Numerous methods have been proposed for nonlinear control. Yin et al. [12] applied an approximated adaptive fuzzy control method to a nonlinear strict-feedback system with unmodeled dynamics. And, Yin et al. [13] developed a fault-tolerant controller based on the adaptive neural network backstepping technique for time-varying delay systems with unmodeled dynamics.The fact that hysteresis is not a smooth process renders the characteristics of a suspension system with hysteretic damping inherently nonlinear. This can greatly complicate the design of a controller. Accurately predicting and/or controlling the performance of a system require that the effects of nonlinear phenomena be taken into account. The problem of rough road profiles and unwanted vehicle vibrations produced by kinematic excitation make it necessary to form accurate predictions of dynamic behavior over a range of operating conditions in the design of suspension systems. Bifurcation and chaos are crucial concepts in the study of stability in nonlinear systems, and numerous researchers have addressed issues pertaining to nonlinear dynamics in suspension systems [14–17] as well as chaotic responses in nonlinear air suspension systems [18, 19]. However, the issue of codimension bifurcations has largely been overlooked.In this study, we examined a suspension system presenting codimension-two bifurcation and homoclinic orbits, both of which determine the nature of the dynamic behavior and the existence of strange attractors. We sought to calculate the strange attractor using numerical simulation in which external forces from a road profile are added to the system as a parameter with a specific range of values. We employed bifurcation diagrams, phase portraits, Poincaré maps, frequency spectra, and Lyapunov exponents to explain periodic and chaotic motions in a vehicle suspension system that exerts hysteretic nonlinear damping forces. Algorithms that computes Lyapunov exponents for smooth dynamic systems [20] were used to determine whether the system exhibits chaotic motion.Chaotic behavior is normally considered undesirable due to the restrictions it imposes on the operating range of electrical and mechanical devices. Many engineering problems can only be resolved by driving a chaotic attractor in a periodic orbit. Numerous methods have been proposed for the control of chaos [21–23], including suspension systems [15, 24, 25]. In this study, we adopted the simple control method proposed by Cai et al. [26], in which state feedback is used to convert chaotic motion into stable motion. Numerical simulations demonstrate the feasibility of the proposed approach.This paper is organized as follows. Section2 describes the proposed vehicle suspension system using a quarter-car model under kinematic excitation by a road surface profile. Section 3 outlines the local bifurcation analysis used to identify codimension-two bifurcation and homoclinic orbits. Section 4 elucidates the complex dynamic behavior in the suspension system using numerical analysis, such as a bifurcation diagrams, phase portraits, Poincaré maps, and frequency spectra. Section 5 presents the Lyapunov exponent used to determine whether the system exhibits chaotic motion. A state feedback control technique for controlling chaos in the suspension system is presented in Section 6. Conclusions are drawn in Section 7.


## 2. 2. Formulation of Problem

Figure1 presents quarter-car model with a single degree-of-freedom and hysteretic nonlinear damping. The equation of motion for the system was formulated by Li et al. [14]:(1)mx¨+k1x-x0+Fh=0,where Fh is an additional nonlinear hysteretic suspension damping and stiffness force, which is dependent on the relative displacement and velocity:(2)Fh=k2x-x03+C1x˙-x˙0+C2x˙-x˙03,x0=a~sin⁡Ωt.In defining the relative vertical displacement as y=x-x0, (1) can be written as follows:(3)y¨+ω2y+B1y3+B2y˙+B3y˙3=Asin⁡Ωt,where ω2=k1/m, B1=k2/m, B2=C1/m, B3=C2/m, A=a~Ω2.Figure 1
Schematic diagram of the quarter-car model.The states are chosen such thaty1 =y, y2 = y˙. Thus, the state-space model for the suspension system can be written as(4a)y˙1=y2,(4b)y˙2=-ω2y1-B1y13-B2y2-B3y23.The amplitude of the road profile excitation isa~=0. Table 1 lists the values of the parameters in the above equations, which were obtained using the methods outlined by Li et al. [14].Table 1
Nomenclature and parameters.
Notation Value Unit Description m 240 kg The body mass k 1 16× 104 N/m The suspension stiffness coefficient k 2 30× 104 N/m3 The nonlinear hysteretic suspension stiffness coefficient C 1 250 Ns/m The linear damping coefficient C 2 −25 Ns3/m3 The nonlinear hysteretic suspension viscous damping coefficient a ~ 0.05 m The amplitude of road excitation x 0  m The road excitation x  m The body’s vertical displacement Ω  rad/s The frequency of road excitation


## 3. Local Codimension-Two Bifurcation Analysis

Our main objective in this analysis is to observe the qualitative behavior of the system and how it changes as system parameters are varied. Changes occurring near the equilibrium point of the system are referred to as local bifurcation. The codimension of bifurcation refers to the minimum number of parameters that must be varied in order to observe a specific type of bifurcation [8].Nonlinear analysis must be based on linear analysis. Bifurcation analysis is based on the Jacobian matrix in (4a) and (4b). When (k1, c1) = (0, 0) and the Jacobian matrix has double-zero eigenvalues of the form of 0100, we have the most degenerate situation (eigenvalues with zero real part), that is, a double-zero codimension-two bifurcation problem. Bifurcation equations can be written as follows:(5a)y˙1=y2,(5b)y˙2=-B1y13+higherorderterms.Bifurcation equations (5a) and (5b) take the normal form because they are the simplest forms of the equations of motion [8, 27–30].The blow-up technique is used to check whether the cubic terms in (5a) and (5b) determine local behavior. The blow-up technique involves a singular change to the coordinates that expands degenerate fixed points into circles containing a finite number of fixed points. If these fixed points are hyperbolic following the first change in the coordinates, then the local flow near the circle (and consequently flow near the original fixed point) remains stable with respect to higher-order terms. If these fixed points are all degenerate, then the coordinates must be changed again. If all of the fixed points become hyperbolic within a finite number of coordinate changes, then a system with only lower-order terms will exhibit the same qualitative behavior as the overall system [8]. The steps used in the blow-up method are as follows:First Blow-Up (6a) y 1 = r cos ⁡ θ , (6b) y 2 = r sin ⁡ θ .Equations (5a), (5b) and (6a), (6b) become(7a)r˙=rsin⁡θ1-B1r2cos2θ+h.o.t.,(7b)θ˙=-sin2θ-B1r2cos4θ+h.o.t..Stability analysis reveals that the system has two degenerate equilibrium points (0, 0) and (0,π), which means that the blow-up procedure must be repeated. In the case of (0, 0), the power series in normal form is(8a)r˙=rθ+h.o.t.,(8b)θ˙=-θ2-B1r2+h.o.t..Second Blow-Up. Let(9a)r=ρsin⁡ϕ,(9b)θ=ρcos⁡ϕ,so that (8a), (8b) and (9a), (9b) become(10a)ρ˙=-ρ2cos⁡ϕB1sin2ϕ+cos2ϕ+h.o.t.,(10b)ϕ˙=ρsin⁡ϕ2cos2ϕ-B1sin2ϕ+h.o.t..After dividing by the unaffected common component,ρ, (10a) and (10b) have two hyperbolic fixed points, (0, 0) and (0, π). Thus, no further blow-up is required.For the other equilibrium point, (0,π), the power series in the normal form is the same as the power series used for (8a) and (8b), such that the same results are obtained. The degenerate field contains cubic terms of the lowest order and unfolding is performed using the two-parameter family, as follows:(11a)y˙1=y2,(11b)y˙2=K1y1+K2y2-B1y13,where K1=-k1/m, K2=-c1/m. The higher-order terms in (11a) and (11b) are disregarded when applied as an alternative for (4a) and (4b).The term “bifurcation set” refers to a set in parameter space on which the system bifurcates. It is obtained from unfolded equations in the normal form.There are three equilibrium points in this case:1.
(0, 0) for allK1, K2;2.
(±K1/B1) for K1>0.By checking the stability of these equilibrium points, we obtain the following:(1  ) For K 1 = 0, a subcritical pitchfork bifurcation occurs if K2 > 0 and a supercritical pitchfork bifurcation occurs if K2<0.(2  ) For K 1 > 0 and K 2 = 0, a rescaling technique and Melnikov’s perturbation method are used to obtain the complete bifurcation set associated with homoclinic orbits in the unfolding bifurcation equations (11a) and (11b) [8].Lety1=ϵu, y2=ϵ2v, K1=ϵ2v1, K2=ϵ3v2, t→ϵt.Thus, after rescaling, (11a) and (11b) are rewritten as(12a)u˙=v,(12b)v˙=v1u+ϵ2vv2-B1u3.Melnikov’s method is then used to characterize the relationship between K1 and K2.After a series of operations, the following expression is obtained:(13a)K2=0,(13b)K1=B1ϵ3.Equations (13a) and (13b) yield the homoclinic saddle connection, referred to as homoclinic bifurcation.(3  ) For K 1 < 0 and K 2 = 0, the Jacobian matrix in (11a) and (11b) is written as(14)DF0,0=01K10,The eigenvalues for DF(0,0) are(15)λ=±K1i.For the following change in coordinates,(16)y1y2=00-K10UV. Equations (12a) and (12b) become(17)U˙V˙=K2-K1K10UV+0B1V3K1.LetgU,V=B1V3/K1 and f(U,V) = 0.Clearly, it can be seen that(18)fvvv=0,fvuu=0gvvu=0,guuu=0fuu=0,fuv=0,fvv=0guu=0,guv=0,gvv=0.According to the Hopf bifurcation theorem [8], coefficienta must not equal zero. In this case, the coefficienta equals zero; therefore, there is no Hopf bifurcation.Numerical simulations using (4a) and (4b) were performed to check the validity of this bifurcation analysis. The parameter space was divided into several regions using bifurcation sets, as shown in Figure 2. The commercial package DIVPRK of IMSL in the FORTRAN subroutines for mathematics applications was used to solve these ordinary differential equations [31]. The resulting numerical simulations are presented in Figures 3–7. The results in each parametric region are discussed below.Figure 2
The bifurcation sets that are associated with typical phase portraits in each region, whereP+ represents a subcritical pitchfork, P- represents a supercritical pitchfork, and SC represents a homoclinic bifurcation.Figure 3
Numerical simulation for region (I), for (K1, K2) = (0.02, 0.01).Figure 4
Numerical simulation for SC, for (K1, K2) = (0.02, 0.0).Figure 5
Numerical simulation for region (II), for (K1, K2) = (0.02, −0.01).Figure 6
Numerical simulation for region (III), for (K1, K2) = (−0.02, −0.01).Figure 7
Numerical simulation for region (V), for (K1, K2) = (−0.02, 0.01).Region (I) includes sources and a saddle. If the initial conditions are in the neighborhood of the two sources, then the mass is distant from the equilibrium points (Figure3). At the boundary between regions (I) and (II) (i.e., in the bifurcation set SC (K2=0, K1>0)), the system is an integrable system with homoclinic orbits (Figure 4). Therefore, if the damping term is disregarded, then the system is conservative. In region (II), the system is not integrable with two sinks and one saddle. If the initial conditions are in the neighborhood of the saddle, then the mass returns to the two sinks, which are stable points, as time approaches infinity (Figure 5). At the boundary between regions (II) and (III) (i.e., in bifurcation set P-), the two sinks move closer to the saddle as the parameters approach bifurcation set, P-, and the stability region shrinks to zero. In region (III), the two sinks disappear and the saddle (0, 0) becomes a sink. This is a supercritical pitchfork bifurcation with an attracting sink region (III) (Figure 6). From regions (III) to (V), the eigenvalues for the Jacobian matrix cross the imaginary axis; however, there is no Hopf bifurcation. This means that the dynamic behavior is linear, which is a similar situation to that when a stable equilibrium point becomes an unstable equilibrium point (Figure 7). Between region (V) and region (I), boundary P+ must be crossed. This is called a subcritical pitchfork bifurcation.In terms of the equilibria and Figure2, K1 determines the number of equilibrium points and K2 determines the local behavior (e.g., stability or instability), and whether there are homoclinic orbits. Clearly, only a positive damping coefficient results in a stable system.


## 4. Analytical Tools for Observing Nonlinear Dynamics in a Vehicle Suspension System

In (3), the amplitude of the road excitation, a~, equals 0.05 m. We performed a series of numerical simulations using (3) to determine the characteristics of this system. Dynamic behaviors can be more clearly observed over a range of parameter values in bifurcation diagrams, which are widely used to describe transitions from periodic motion to chaotic motion in dynamic systems. The resulting bifurcation diagram is presented in Figure 8. This figure clearly shows chaotic motion in regionA. RegionC presents period-1 motion and there is a stable equilibrium point. In this region, each response was characterized using a phase portrait, a Poincaré map (velocity versus phase angle), and a frequency spectrum. Figure 9 shows that the equilibrium point of the system is stable if parameter, Ω, exists in regionC. It is well-known that road excitation covers a wide range frequencies, and higher vehicle velocities increase excitation from the road surface. The range of human sensitivity to frequency is from 4 to 8 Hz (about 25.12 rad/s to 50.24 rad/s); therefore, we limited this investigation to the case of Ω = 40.0 rad/s. Figure 10 shows the periodic motion in the system in cases where the parameter exists in regionC (i.e., Ω = 8.2 rad/s and Ω = 40 rad/s). The amplitude of road excitation is defined as 0.05 m, while the response of the system (even in periodic motion) has an amplitude of near 0.35 m (Figure 9). Clearly, this response is not acceptable. When parameter Ω exists in regionB, quasiperiodic motion occurs in this region. As shown in Figures 11–13, the quasiperiodic trajectory is restricted to an annular-like region of the state space and differs considerably from periodic trajectories. When Ω = 8.0, 8.05, or 8.15, the system exhibits a quasiperiodic motion and the map forms a continuous closed orbit in the Poincaré section, as shown in Figures 11–13. As the forcing frequency (Ω) continues decreasing into regionA (Figure 14), chaotic motion produces chatter vibration. The transition sequence for the system in this paper is from periodic motion to quasiperiodic motion and finally to chaotic motion. The particular features of these two descriptors characterize the chaotic behavior. The Poincaré map contains an infinite set of points, which are referred to as a strange attractor. The frequency spectrum of the chaotic motion is a continuous broad spectrum. The strange attractor and continuous Fourier spectrum are strong indicators of chaos. Figure 14 shows clear indications of chaotic behavior in regionA. Figures 14(a)–14(c), respectively, show the phase portraits, the Poincaré maps, and frequency spectra.Figure 8
Bifurcation diagram for the system when the amplitude of the road excitation,a~, is 0.05 m: (a) Ω = 7.96 rad/s to 40 rad/s, (b) Ω = 7.96 rad/s to 8.25 rad/s.
(a) (b)Figure 9
The period-one orbit atΩ = 8.2 rad/s: (a) the phase portrait, (b) the Poincaré map, and (c) the frequency spectrum.
(a) (b) (c)Figure 10
Phase portrait: (a)Ω = 8.2 rad/s, (b) Ω = 40 rad/s.
(a) (b)Figure 11
Quasiperiodic orbit atΩ = 8.0 rad/s: (a) the phase portrait, (b) the Poincaré map, and (c) the frequency spectrum.
(a) (b) (c)Figure 12
Quasiperiodic orbit atΩ = 8.05 rad/s: (a) the phase portrait, (b) the Poincaré map, and (c) the frequency spectrum.
(a) (b) (c)Figure 13
Quasiperiodic orbit atΩ = 8.15 rad/s: (a) the phase portrait, (b) the Poincaré map, and (c) frequency spectrum.
(a) (b) (c)Figure 14
Chaotic motion atΩ = 7.964 rad/s: (a) the phase portrait, (b) Poincaré map, and (c) the frequency spectrum.
(a) (b) (c)


## 5. Analysis of Chaotic Phenomena: Lyapunov Exponent and Lyapunov Dimension

The analysis presented in Section4 is insufficient to confidently identify chaotic motion in a suspension system. In this section, we adopt Lyapunov exponents to verify the existence of chaotic motion. In every dynamic system, the spectrum of Lyapunov exponents (λ) [20] indicates variations in the lengths, areas, and volumes of the phase space. Verifying the existence of chaos requires calculation of only the largest exponent to indicate whether nearby trajectories diverge (λ>0) or converge (λ<0), on average. Any bounded motion within a system that has at least one positive Lyapunov exponent is defined as chaotic. The Lyapunov exponents for periodic motion are not positive.Figure15 plots the evolution of the largest Lyapunov exponent in a suspension system. This value is computed using the algorithm presented by Wolf et al. [20]. Figure 14 shows that chaotic motion began at approximately Ω=7.965 rad/s, because at point (P1), the sign of the largest Lyapunov exponent changed from negative to positive as the parameter, Ω, slowly decreased. Between points P1 and P2, the largest Lyapunov exponent approached zero. The existence of a forcing frequency, Ω, in these regions denotes the motion in the system as quasiperiodic. When the value of Ω is greater than P1 (such as when Ω=8.2 rad/s), the Lyapunov exponents obtained using (3) are λ1=-0.0010375 and λ2=-0.0013232. The sum of these values is λ1+λ2=-0.0023607; this negative value demonstrates that the motion of mass at these values eventually converges to a stable limit cycle. Using the Lyapunov exponents for a dynamical system where λ1≥⋯≥λn, Kaplan and Yorke [32] provided an estimation of the Lyapunov dimension, dL, as follows:(19)dL=j+1λj+1∑i=1jλi, where j is the largest integer that satisfies ∑i=1jλi>0. This technique yields a Lyapunov dimension of dL=1 for Ω = 8.2 rad/s. The fact that the Lyapunov dimension is an integer means that the system has periodic motion. When parameter Ω increases beyond the bifurcation point (such as when Ω = 7.964 rad/s), the Lyapunov exponents are λ1=0.0054396 and λ2=-0.092153, and the Lyapunov dimension is dL=1.059. In this case, the Lyapunov dimension is not an integer, indicating that the system exhibits chaotic motion.Figure 15
Evolution of the largest Lyapunov exponent.


## 6. Controlling Chaos

Improving the performance of a dynamic system and/or avoiding chaotic behavior require periodic motion, which is more important when working under specific conditions. Cai et al. [26] suggested a simple, effective control method that converts chaotic motions into periodic motion using the linear state feedback as a system variable. This approach to achieving ann-dimensional dynamic system is explained briefly as follows:(20)x˙=fx,t,where x(t)∈Rn is the state vector and f=(f1,…,fi,…,fn), where fi is a linear or a nonlinear function and f includes at least one nonlinear function. If fk(x,t) is a nonlinear function that gives chaotic motion in (20), then only one term of the state feedback for the available system variable, xm, is added to the equation that includes fk, as follows:(21)x˙k=fkx,t+Gxm,k,m∈1,2,…,n, where G is the feedback gain. Other functions remain in their original forms.System (3) with state feedback control can be written as follows:(22)dy1dt=y2dy2dt=-ω2y1-B1y13-B2y2-B3y23+a~Ω2sin⁡Ωt+Gy2. In the absence of state feedback control,G = 0 and (3) describes chaotic motion when Ω = 7.964 rad/s. Figure 16 shows the resulting bifurcation diagram for a system with state feedback control. This figure clearly shows that chaotic motion occurs at approximately G>-0.14 and disappears at approximately G≤-0.14. When feedback gain, G, falls below -0.14, stable periodic motion occurs. As shown in Figure 17, a control signal is applied to the system at 35 s to govern the chaotic oscillation. In this manner, state feedback control can be used to suppress chaotic motion and improve the performance of a vehicle’s suspension system. The form of the control gain in the original system equation (1) is as follows:(23)G~=G·m. The physical meaning of gain, G~, is control gain and the unit of the gain is N·s/kg·m.Figure 16
The bifurcation diagram for a system with state feedback control, whereG denotes the feedback gain.Figure 17
Transforming chaotic motion to a desired period-one orbit forG=-0.3 and Ω=7.964 rad/s: (a) the phase portrait for the controlled system when the state feedback control signal is introduced after 35 s and (b) the time responses of the velocity.
(a) (b)In many general systems that involve chaotic dynamics, this state feedback technique does not require preliminary or on-line analysis of the system dynamics and can be implemented by adding the feedback of a suitable variable to the original system with sufficient control gain to suppress the development of chaos. In this regard, the state feedback control technique is simple and effective approach to chaos suppression.


## 7. Conclusions

This study examined the nonlinear dynamic behavior in a vehicle suspension system with hysteretic damping. Our objective was to develop a method to control chaotic vibration. Dynamic behavior was observed over the entire range of parameter values of the bifurcation sets and bifurcation diagram. Using bifurcation analysis, we identified a codimension-two bifurcation structure with double-zero degeneracy in the system. The bifurcation diagram indicates the existence of quasiperiodic orbit and chaotic motion, thereby demonstrating the effectiveness of using Lyapunov exponents as a tool for identifying chaotic motion. Finally, we proposed a state feedback control method to suppress chaotic motions and improve the ride provided by an automotive suspension system. Simulation results confirm the feasibility of the proposed method. Our results provide insight into the operating conditions under which undesirable dynamic motion takes place in vehicle suspension systems. Our results also provide a useful reference for engineers involved in the design of control systems for vehicle suspensions. It should be noted that in practical engineering, phenomena such as actuator faults, time-varying delays, unmodeled dynamics, and dynamics disturbances often occur simultaneously, and these problems are beyond the scope of the current study. In future work, we will seek to control nonlinear systems with time-varying delays.



---

*Source: 4218175-2018-03-27.xml*