1 Introduction

In the last decade a remarkable amount of researchers have concentrated on the analysis of nonlinear quarter-car suspension systems. These complicated physical systems include various distinguished components and are influenced by some external road excitation forces. A magnetorheological (MR) damper constructs the main structure of a suspension system, and the principle of MR-fluid expands its basic physical behaviors. In fact, an MR-fluid fills up the MR-damper and is controlled by an electromagnetic field. In [1] the dynamical behavior of a vehicle with an MR-damper was studied. In [2] the chaotic behavior of a half-vehicle model with semiactive suspension system was investigated, and a discrete optimal controller was designed for chaotic oscillations. In [3] a sensitivity control was proposed for an MR-damper semiactive suspension system, which was modeled by a sixth-order polynomial. In all these studies the considered chaotic systems were modeled in the framework of ordinary differential equations (ODEs). However, ODEs may not be able to specify the complicated dynamical behaviors and many hidden physical features of nonlinear chaotic oscillators.

Fractional calculus (FC), as an extension of integer-order domain, is an efficient tool for the modeling of real-world phenomena with complex physical dynamics [4–8]. Nowadays, a considerable number of researches illustrated that fractional-order differential equations (FDEs) can reveal complex dynamical features more accurately than ODEs [9–11]. This is due to the fact that fractional-order derivatives and integrals are capable to characterize the properties of memory effects as an essential aspect in many real-world phenomena [12–16]. Latterly, FC has been distinguished as an interesting field of modeling and analysis of nonlinear dynamical systems with chaotic features [17–19]. In many researches, it was shown that FC is efficient to illustrate hereditary characteristics and memory effects [20–22]; however, due to the existence of singular kernel in the classic framework of fractional operators, some nonlinear complicated dynamics may not be identified correctly by the classical descriptions. Therefore, to analyze nonlocal dynamics more accurately, new kinds of fractional derivatives including nonsingular kernels were introduced; two of the most efficient and practicable ones are known as Caputo–Fabrizio (CF) [23] and Atangana–Baleanu–Caputo (ABC) [24] derivatives with exponential and Mittag-Leffler (ML) kernels, respectively. These new fractional-order differential operators demonstrate different asymptotic behaviors, and sometimes they are more accurate than their classic counterparts [25, 26]. Moreover, they play an outstanding role in clarifying the hidden characteristics of real-world nonlocal dynamical processes [27, 28]. Nonetheless, the properties of these operators are needed to be carefully analyzed, and efficacious analytical/numerical methods should be developed to make these operators more practical in real-world systems [29, 30]. Motivated by the above argumentation, the main aim of this paper is to apply the CF fractional operator to model and investigate a chaotic suspension system with a sine function as its road excitation force. Considering various orders of the CF differential operator, we recognize the chaotic behavior of this suspension system in the phase-portrait and time-domain states. Comparing the system behavior in both integer- and fractional-order forms, we realize significant differences between the integer-order model and its fractional-order counterpart. Moreover, some hidden features of the considered system are drawn out from the proposed fractional-order model, whereas these characteristics are indistinguishable when ODEs are used. Also, we develop an efficient quadratic numerical method to solve the suspension system-related FDEs. More to the point, we design a state-feedback chaos controller and a chaos optimal control to suppress the negative effect of chaotic oscillations. In the case of designing the optimal controller, we combine the proposed numerical method with an iterative forward-backward algorithm to solve the related fractional optimal control problem (FOCP). To the best of our knowledge, this is the first study that analyzes and controls the chaotic behavior of a nonlinear suspension system in a nonsingular fractional-order framework. Simulation results also verify the effectiveness of our proposed strategy.

In the rest of this paper, we first present some mathematical definitions and preliminaries in Sect. 2. Then a new fractional-order mathematical model for the suspension system and its principal features are introduced in Sect. 3. A quadratic numerical approach for solving the system of FDEs is suggested in Sect. 4. Next, some numerical simulations are given in Sect. 5. Afterward, a chaos state-feedback controller is suggested in Sect. 6. Then an optimal controller and an iterative algorithm for solving the related FOCP are designed in Sect. 7, and lastly, the paper is concluded in Sect. 8.

2 Mathematical preliminaries

In this section, we briefly introduce preliminaries considering the CF fractional operator and investigate its significant features [31].

Definition 2.1

For \(f\in \mathbb{H}^{1}{(0,T)}\) and \(0< q<1\), the CF fractional derivative and its related fractional integral operator are, respectively, presented by

$$\begin{aligned}& {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f(t)}= \frac{1}{1-q} \int _{0}^{t} {\exp {\bigl[-{\beta }(t- \tau ) \bigr]}}\frac{df(\tau )}{d\tau }d{\tau }, \end{aligned}$$
(1)
$$\begin{aligned}& {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q}{f(t)}=(1-q)f(t)+{q} { \int _{0}^{t}{f(\tau )}d{\tau }}, \end{aligned}$$
(2)

where \(\beta =\frac{q}{1-q}\).

Note that Eq. (2) yields the integer-order integral when \(q=1\), and it leads to \(f(t)\) for \(q=0\). In the following, we explain some beneficial formulas of the CF fractional operators [31, 32]:

  • For \(f_{1},f_{2} \in \mathbb{H}^{1} (0,T)\) and \(c_{1},c_{2} \in \mathbb{R}\), the CF derivative and integral are linear operators:

    $$\begin{aligned}& {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{{ \bigl(c_{1}f_{1}(t)+c_{2}f_{2}(t) \bigr)}}=c_{1}~{{}^{CF}_{\:\: 0}\hspace{-0.1cm} \mathscr{D}_{t}^{q}} {{f_{1}(t) }}+c_{2}~{{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}} {{f_{2}(t) }}, \end{aligned}$$
    (3)
    $$\begin{aligned}& {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q}{f(t)} {{\bigl(c_{1}f_{1}(t)+c_{2}f_{2}(t) \bigr)}}=c_{1}~{{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q}} {{f_{1}(t) }}+c_{2}~{{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q}} {{f_{2}(t) }}. \end{aligned}$$
    (4)

  • For a constant function \(f(t)\equiv f_{c}\), the CF differential operator is zero, that is,

    $$ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{c}=0}. $$
    (5)

  • Since the CF fractional derivative (1) is the convolution integral of \(\frac{df(t)}{dt}\) and \(\exp (-\frac{q}{1-q}t)\), applying the convolution theorem for the Laplace transform, we get

    $$ \mathcal{L} {{}\bigl[ } {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f(t)} {\bigr] }=\frac{sF(s)-f(0)}{s+q(1-s)}, $$
    (6)

    where \(F(s)=\mathcal{L}{{}[ }{f(t)}] \).

  • The antiderivative property among the CF derivative (1) and its corresponding integral (2) is

    $$ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q}{{} \bigl[ } {{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f(t)}} {\bigr] }=f(t)-f(0). $$
    (7)

  • Consider \(\Vert f(t) \Vert =\max_{0 < t < T} \vert f(t) \vert \) for \(f \in \mathbb{H}^{1}(0,T)\); then the CF differential operator (1) satisfies the inequality

    $$ \bigl\Vert {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f(t)} \bigr\Vert \leq \frac{1}{1-q} \bigl\Vert f(t) \bigr\Vert . $$
    (8)

  • For any \(f_{1},f_{2} \in \mathbb{H}^{1} (0,T)\), there exists a constant \(L>0\) such that the Lipschitz condition is satisfied for the CF differential operator as follows:

    $$ \bigl\Vert {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}-{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{2}(t)} \bigr\Vert \leq L \bigl\Vert f_{1}(t)-f_{2}(t) \bigr\Vert . $$
    (9)

3 New model

In this section, we investigate a quarter-car model of a nonlinear suspension system and its dynamical equation of motion. Based on the principle of D’Alembert [33], the suspension system dynamical equation is described by the integer-order ODE

$$\begin{aligned}& m\ddot{y}(t)+c_{1}\bigl(y(t)-y_{0}(t) \bigr) \\& \quad + \underbrace{c_{2}\bigl(y(t)-y_{0}(t) \bigr)^{3}+\alpha _{1}\bigl(\dot{y}(t)- \dot{y_{0}}(t)\bigr)+\alpha _{2}\bigl(\dot{y}(t)- \dot{y_{0}}(t)\bigr)^{3}}_{Y_{h}(t)}=0, \end{aligned}$$
(10)

where \(y(t)\) is the displacement of the mass m in the vertical direction, \(c_{1}\) is the stiffness coefficient of suspension, \(Y_{h}(t)\) is the hysteretic nonlinear effective damping force, and \(\alpha _{1}\), \(c_{2}\), \(\alpha _{2}\) are positive constants. Also, \(y_{0}=A\sin ({\Omega }t)\) is the road kinematical excitation with frequency Ω. Supposing a new coordinate \(f(t)=y(t)-y_{0}(t)\), we get the following equation of motion:

$$ \ddot{f}(t)+k_{1}\dot{f}(t)+k_{2} \dot{f}^{3}(t)+k_{3}f(t)+k_{4}f^{3}(t)=A \Omega ^{2} \sin (\Omega t), $$
(11)

where \(k_{1}=\alpha _{1}/m\), \(k_{2}=\alpha _{2}/m\), \(k_{3}=c_{1}/m\), and \(k_{4}=c_{2}/m\). Now introducing the new states \(f_{1}(t)=f(t)\) and \(f_{2}(t)=\dot{f}(t)\), we rewrite Eq. (11) in the form [34]

$$ \textstyle\begin{cases} \dot{f_{1}}(t)=f_{2}(t), \\ \dot{f_{2}}(t)=-k_{1}f_{2}(t)-k_{2}f_{2}^{3}(t)-k_{3}f_{1}(t)-k_{4}f_{1}^{3}(t)+A \Omega ^{2} \sin (\Omega t). \end{cases} $$
(12)

As mentioned in [35], many principles of real-world phenomena cannot be described and analyzed by the theory of calculus of variations. For example, the presented dynamical model (12) does not include memory effects as the substantial aspects of many physical systems with complex dynamics. Nevertheless, FC has prevailed this restriction as it can specify the complicated behavior of many physical phenomena involving the effects of memory and other hereditary features. Thus, in the following, we give a new dynamical model in fractional sense for the considered suspension system, in which the integer-order derivatives in (12) are substituted by the CF fractional operator (1). In other words, by fractionalizing (12), we have the following new model:

$$\begin{aligned} \textstyle\begin{cases} {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=f_{2}(t), \\ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=-k_{1}f_{2}(t)-k_{2}f_{2}^{3}(t)-k_{3}f_{1}(t)-k_{4}f_{1}^{3}(t)+A \Omega ^{2} \sin (\Omega t). \end{cases}\displaystyle \end{aligned}$$
(13)

It is worth noting that this fractional dynamical system is reduced to its integer-order counterpart (12) as \(q\to {1}\).

4 Quadratic numerical method

In this section, we suggest a quadratic numerical method to efficiently solve the fractional-order model (13). To this end, we first convert the FDEs (13) into their corresponding integral equations. Afterward, based on the trapezoidal approach, we introduce a quadratic numerical scheme to approximately solve the related fractional-order integral equations. To introduce the proposed method in the sense of CF, Eq. (13) is first compacted in a well-set form:

$$ \textstyle\begin{cases} {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{{F(t) }} =\Phi (F(t)), & 0\leq t\leq T< \infty , \\ F(0)=F_{0}, \end{cases} $$
(14)

where \(F(t)=(f_{1}(t),f_{2}(t))\), Φ is the real-valued continuous vector function

Φ ( F ( t ) ) = [ f 2 ( t ) − k 1 f 2 ( t ) − k 2 f 2 3 ( t ) − k 3 f 1 ( t ) − k 4 f 1 3 ( t ) + A Ω 2 sin ( Ω t ) ] ,
(15)

satisfying the Lipschitz condition

$$ \bigl\Vert \Phi \bigl(F_{1}(t)\bigr)-\Phi \bigl(F_{2}(t)\bigr) \bigr\Vert \leq {L} \bigl\Vert {F_{1}(t)-F_{2}(t)} \bigr\Vert , \quad L>0, $$
(16)

and \(F_{0}=(f_{1}(0),f_{2}(0))\) includes the initial points. Applying the CF integral operator (2) to Eq. (14) yields

$$ F(t)=F_{0}+{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q}{\Phi }\bigl(F(t)\bigr),\quad 0\leq t \leq T< \infty . $$
(17)

Now we define a uniform mesh \(0=t_{0}< t_{1}<\cdots <t_{M}=T\) on \([0,T]\) with the time step \(\hbar =\frac{T}{M}=t_{m+1}-t_{m}\), where \(M>0\) is an arbitrary integer, and \(t_{m}=m\hbar \), \(m=0,1,2,\ldots ,M-1\). Moreover, we show the approximation of \(F(t_{m})\) by \(F_{m}\). Then we approximate \(\Phi (F(\tau ))\) on every subinterval \([t_{k},t_{k+1}]\) by a piecewise linear interpolation as follows:

$$ \Phi \bigl(F(\tau )\bigr)|_{[t_{k},t_{k+1}]}\approx \frac{ t_{k+1}-\tau }{t_{k+1}-t_{k}}\Phi (F_{k})+ \frac{\tau -t_{k}}{t_{k+1}-t_{k}}\Phi (F_{k+1}), \quad 0 \leq k \leq m. $$
(18)

To have a numerical solution for Eq. (17), we consider the following discretization of CF integral operator (2):

$$ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{I}_{t}^{q} \Phi \bigl(F(t_{m+1})\bigr)=(1-q)\Phi \bigl(F(t_{m+1})\bigr)+q \int _{0}^{t_{m+1}} \Phi \bigl(F(\tau )\bigr)\,d\tau . $$
(19)

Substituting the approximation of \(\Phi (F(\tau ))\) from (18) into (19), we provide the approximation

$$ F_{m+1}=F_{0}+(1-q)\Phi (F_{m+1})+q \hbar \sum_{k=0}^{m+1}a_{m+1,k} \Phi (F_{k}), \quad m=0,1,\ldots ,M-1, $$
(20)

where the coefficient \(a_{m+1,k}\) is as follows:

$$ \textstyle\begin{cases} a_{m+1,0}=1/2, \\ a_{m+1,k}=1, \quad k=1,2,\ldots,m, \\ a_{m+1,m+1}=1/2. \end{cases} $$
(21)

We analyze the stability and convergence of the numerical approach (20) in the following theorems.

Theorem 4.1

The numerical scheme (20) is conditionally stable.

Proof

Let the perturbation of \(F_{0}\) and \(F_{m}\) (\(m = 0,1,\ldots,M-1\)) be denoted by \(\widetilde{F}_{0}\) and \(\widetilde{F}_{m}\), respectively. Then for the numerical scheme (20), we have

$$ F_{m+1}+\widetilde{F}_{m+1}=F_{0}+ \widetilde{F}_{0}+(1-q)\Phi (F_{m+1}+ \widetilde{F}_{m+1})+q \hbar \sum_{k=0}^{m+1}a_{m+1,k} \Phi (F_{k}+ \widetilde{F}_{k}). $$
(22)

Substituting Eq. (20) into Eq. (22) results in

$$\begin{aligned} \widetilde{F}_{m+1} =&\widetilde{F}_{0}+(1-q) \bigl(\Phi (F_{m+1}+ \widetilde{F}_{m+1})-\Phi (F_{m+1})\bigr) \\ &{}+q\hbar \sum_{k=0}^{m+1}a_{m+1,k} \bigl( \Phi (F_{k}+\widetilde{F}_{k})-\Phi (F_{k})\bigr). \end{aligned}$$
(23)

Using the triangle inequality and Lipschitz condition, from Eq. (23) we get

$$\begin{aligned} \Vert \widetilde{F}_{m+1} \Vert =& \Biggl\Vert \widetilde{F}_{0}+(1-q) \bigl( \Phi (F_{m+1}+ \widetilde{F}_{m+1})-\Phi (F_{m+1})\bigr) \\ &{}+q\hbar \sum_{k=0}^{m+1}a_{m+1,k} \bigl( \Phi (F_{k}+\widetilde{F}_{k})-\Phi (F_{k})\bigr) \Biggr\Vert \\ \leq& \Vert \widetilde{F}_{0} \Vert +(1-q) \bigl\Vert \Phi (F_{m+1}+ \widetilde{F}_{m+1})-\Phi (F_{m+1}) \bigr\Vert \\ &{}+q\hbar \sum_{k=0}^{m+1}a_{m+1,k} \bigl\Vert \Phi (F_{k}+\widetilde{F}_{k})-\Phi (F_{k}) \bigr\Vert \\ \leq& \Vert \widetilde{F}_{0} \Vert +(1-q)L \Vert \widetilde{F}_{m+1} \Vert +q\hbar L\sum_{k=0}^{m+1}a_{m+1,k} \Vert \widetilde{F}_{k} \Vert \\ \leq& \Vert \widetilde{F}_{0} \Vert +q\hbar L\sum _{k=0}^{m}a_{m+1,k} \Vert \widetilde{F}_{k} \Vert +L\bigl((1-q)+q\hbar a_{m+1,m+1} \bigr) \Vert \widetilde{F}_{m+1} \Vert . \end{aligned}$$
(24)

Considering the coefficients from Eq. (21), we have

$$ \Vert \widetilde{F}_{m+1} \Vert \leq \Vert \widetilde{F}_{0} \Vert +q\hbar L\sum_{k=0}^{m} \Vert \widetilde{F}_{k} \Vert +L\biggl(1-q+ \frac{q\hbar }{2}\biggr) \Vert \widetilde{F}_{m+1} \Vert . $$
(25)

For all parameters L, ħ, q such that \(L(1-q +\frac{q\hbar }{2}) < 1\), we have

$$ \Vert \widetilde{F}_{m+1} \Vert \leq \phi (q, \hbar ) \Vert \widetilde{F}_{0} \Vert +\phi (q,\hbar )q\hbar L \sum_{k=0}^{m} \Vert \widetilde{F}_{k} \Vert , $$
(26)

where

$$ \phi (q,\hbar )=\frac{1}{1-L(1-q+\frac{\hbar q}{2})}. $$
(27)

Furthermore, there exists a constant \(C_{\phi }\) such that for a sufficiently small ħ, we have

$$ 1 < \phi (q,\hbar ) < C_{\phi }. $$
(28)

Then

$$ \Vert \widetilde{F}_{m+1} \Vert \leq C_{\phi } \Vert \widetilde{F}_{0} \Vert +C_{\phi }q\hbar L \sum_{k=0}^{m} \Vert \widetilde{F}_{k} \Vert . $$
(29)

Eventually, from Lemma 3.3 in [36] and applying the Grönwall inequality, we get \(\Vert \widetilde{F}_{m+1} \Vert \leq C_{1} \Vert \widetilde{F}_{0} \Vert \), in where \(C_{1}\) is a general constant. □

Theorem 4.2

The numerical scheme (20) is conditionally convergent.

Proof

The difference between the exact solution \(F(t_{m+1})\) and the approximate solution \(F_{m+1}\) (Eq. (20)) is computed by

$$\begin{aligned} F(t_{m+1})-F_{m+1} =&(1-q) \bigl(\Phi \bigl(F(t_{m+1})\bigr)-\Phi (F_{m+1})\bigr) \\ &{}+q \Biggl[ \int _{0}^{t_{m+1}}\Phi \bigl(F(\tau )\bigr)\,d\tau - \hbar \sum_{k=0}^{m+1} a_{m+1,k} \Phi (F_{k}) \Biggr] \\ =&(1-q) \bigl(\Phi \bigl(F(t_{m+1})\bigr)-\Phi (F_{m+1}) \bigr) \\ &{}+q \Biggl[ \int _{0}^{t_{m+1}} \Phi \bigl(F(\tau )\bigr)\,d\tau - \hbar \sum_{k=0}^{m+1}a_{m+1,k} \Phi \bigl(F(t_{k})\bigr) \Biggr] \\ &{} +q \hbar \sum_{k=0}^{m+1}a_{m+1,k} \bigl(\Phi \bigl(F(t_{k})\bigr)- \Phi (F_{k}) \bigr). \end{aligned}$$
(30)

Applying the triangle inequality and Lipschitz condition, together with Lemma 3.1 in [36], we obtain

$$\begin{aligned} \bigl\Vert F(t_{m+1})-F_{m+1} \bigr\Vert \leq& L(1-q) \bigl\Vert F(t_{m+1})-F_{m+1} \bigr\Vert + \frac{1}{2}qCT\hbar ^{2} \\ &{} +q \hbar L\sum_{k=0}^{m+1}a_{m+1,k} \bigl\Vert F(t_{k})- F_{k} \bigr\Vert \\ \leq& L\biggl(1-q+\frac{q\hbar }{2}\biggr) \bigl\Vert F(t_{m+1})-F_{m+1} \bigr\Vert + \frac{1}{2}qCT \hbar ^{2} \\ &{} +q \hbar L\sum_{k=0}^{m}a_{m+1,k} \bigl\Vert F(t_{k})- F_{k} \bigr\Vert , \end{aligned}$$
(31)

where C is a general constant. Considering the coefficients in Eq. (21) and all admissible parameters L, ħ, q satisfying the inequality \(L(1-q +\frac{q \hbar }{2}) < 1\), we have

$$ \bigl\Vert F(t_{m+1})-F_{m+1} \bigr\Vert \leq \phi (q,\hbar ) \frac{1}{2}qCT \hbar ^{2}+ \phi (q,\hbar ) q\hbar L \sum_{k=0}^{m} \bigl\Vert F(t_{k})- F_{k} \bigr\Vert , $$
(32)

where \(\phi (q,\hbar )\) is taken into account from Eq. (27) and satisfies inequality (28). Finally, from Lemma 3.3 in [36] and implementing the Grönwall inequality, we derive

$$ \bigl\Vert F(t_{m+1})-F_{m+1} \bigr\Vert \leq C_{1}\hbar ^{2}, $$
(33)

where \(C_{1}\) is a generic constant. □

5 Numerical simulations

In this section, we apply the proposed numerical method to solve the fractional suspension system (13) by considering various values of q. In these simulations, we investigate the related chaotic and nonchaotic motions with the road excitation force considered as a sinusoidal function with \(A=0.05\) and \(\Omega =8\). Figures 1–3 represent the dynamical states and the phase-portraits for \(q=0.992\), 0.993, 0.994, 0.995, 0.996, 0.997, 0.998, 0.999, 1. From these figures we can observe that the system displays periodic solutions for smaller values of q. However, the noninteger-order model presents chaotic oscillations by increasing the value of q, a fact that is distinguished from the periodic manner appeared in the previous plots. Considering greater values of fractional order, we see disorderly behaviors in fractional sense; indeed, the state variables in the phase plane are extended randomly and never converge to a specific region. Accordingly, the dynamical system in the form of fractional model depicts a chaotic behavior in this case. Finally, we see that the simulation results in fractional mode come close to those of integer-order case when the fractional order goes to 1. In this case the system states are irregular, and the phase-portraits establish numerous loops. By these properties we can conclude that the considered dynamical system modeled in terms of integer-order calculus presents only chaotic oscillations. Concerning the fractional model (13), if the order of the fractional operator goes to 1, the system model (13) behaves chaotically like the integer-order one. However, as the fractional order q diverges from 1, the dynamical behavior is completely changed. Furthermore, some dynamical features achieved by the integer-order model present a greater divergence from the real system behavior than the fractional-order model [20]. These results confirm that the fractional-order models can describe nonlinear systems with complicated dynamics more accurately than the integer-order counterparts.

Figure 1
figure 1

The state variable \(f_{1}(t)\) of the fractional suspension system (13) for \(q=0.992,0.993,0.994,0.995, 0.996,0.997,0.998,0.999,1\)

Full size image
Figure 2
figure 2

The state variable \(f_{2}(t)\) of the fractional suspension system (13) for \(q=0.992,0.993,0.994,0.995,0.996, 0.997,0.998,0.999,1\)

Full size image
Figure 3
figure 3

The phase-portraits of the fractional suspension system (13) for \(q=0.992,0.993,0.994,0.995,0.996, 0.997,0.998,0.999,1\)

Full size image

6 State-feedback chaos control

In this section, we design a state-feedback controller to suppress the disorderly undesirable behaviors of the nonlinear suspension system (13). To this aim, we apply the state-feedback inputs \(u_{1}\), \(u_{2}\) to the fractional-order model (13) as follows:

$$\begin{aligned} \textstyle\begin{cases} {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=f_{2}(t)-u_{1}f_{1}(t), \\ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=-k_{1}f_{2}(t)-k_{2}f_{2}^{3}(t)-k_{3}f_{1}(t)-k_{4}f_{1}^{3}(t) \\ \hphantom{{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=}{}+A \Omega ^{2} \sin (\Omega t)-u_{2}f_{2}(t), \end{cases}\displaystyle \end{aligned}$$
(34)

where \(u_{1},u_{2}>0\) are the state-feedback gain constants. We apply the numerical method suggested in Sect. 4 to solve the FDEs (34) in the presence of the state-feedback controller with \(u_{1}=u_{2}=5\). Figures 4–6 show the state trajectories and the phase-portraits. These figures indicate that the controlled suspension system (34) displays nonchaotic behaviors for all values of q. Indeed, the phase-portraits in Fig. 6 converge to a specific region, and the state variables \(f_{1}(t)\) and \(f_{2}(t)\) in the time domain (Figs. 4–5) depict a stable periodic solution. These results confirm that the state-feedback controller overcomes the irregular behavior of the chaotic suspension system under investigation.

Figure 4
figure 4

The state variable \(f_{1}(t)\) of the controlled system (34) for \(q=0.992,0.993,0.994,0.995,0.996, 0.997,0.998,0.999,1\)

Full size image
Figure 5
figure 5

The state variable \(f_{2}(t)\) of the controlled system (34) for \(q=0.992,0.993,0.994,0.995,0.996,0.997, 0.998,0.999,1\)

Full size image
Figure 6
figure 6

The phase-portraits of the controlled system (34) for \(q=0.992,0.993,0.994,0.995,0.996,0.997, 0.998,0.999,1\)

Full size image

7 Optimal chaos control

In this section, we design an optimal chaos controller to suppress the chaotic behavior of the considered suspension system. To do this, we take into account the following generalized model of Eq. (13) including the control inputs \(u_{1}(t)\), \(u_{2}(t)\):

$$\begin{aligned} \textstyle\begin{cases} {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=f_{2}(t)+u_{1}(t), \\ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}=-k_{1}f_{2}(t)-k_{2}f_{2}^{3}(t)-k_{3}f_{1}(t)-k_{4}f_{1}^{3}(t)+A \Omega ^{2} \sin (\Omega t)+u_{2}(t). \end{cases}\displaystyle \end{aligned}$$
(35)

The aim is designing the optimal control inputs \(u_{1}^{*}(t)\), \(u_{1}^{*}(t)\) along with minimizing the quadratic objective index

$$ J=\frac{1}{2} \int _{0}^{T}q_{1}f_{1}^{2}( \tau )+q_{2}f_{2}^{2}(\tau )+r_{1}u_{1}^{2}( \tau )+r_{2}u_{2}^{2}( \tau )\,d\tau , $$
(36)

where \(q_{i}\geq \) and \(r_{i}>0\) are constant weights. To solve the suspension system-related FOCP (35)–(36), we should derive the corresponding optimality necessary conditions. To do so, considering the optimal control theory, we express the scalar Hamiltonian function

$$ \mathscr{H}(t)=\mathscr{L}_{0}(t)+\sum _{1}^{2}g_{i}(t) \mathscr{L}_{i}(t), $$
(37)

where \(\mathscr{L}_{0}\) is the integrand function in (36), \(\mathscr{L}_{i}\) is the right-hand side of the ith equation in (35), and \(g_{i}\) is the Lagrange multiplier, also known as the costate variable. Based on the results available in [37, 38], we write the necessary optimality conditions of the FOCP (35)–(36) as follows:

$$\begin{aligned}& \textstyle\begin{cases} {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{1}(t)}= \frac{\partial \mathscr{H}}{\partial {g_{1}}}(t)=f_{2}(t)+u_{1}^{*}(t), \\ {}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{2}(t)}= \frac{\partial \mathscr{H}}{\partial {g_{2}}}(t) \\ \hphantom{{}^{CF}_{\:\: 0} \hspace{-0.1cm}\mathscr{D}_{t}^{q}{f_{2}(t)}}=-k_{1}f_{2}(t)-k_{2}f_{2}^{3}(t)-k_{3}f_{1}(t)-k_{4}f_{1}^{3}(t)+A \Omega ^{2} \sin (\Omega t)+u_{2}^{*}(t), \end{cases}\displaystyle \end{aligned}$$
(38)
$$\begin{aligned}& \textstyle\begin{cases} {}^{CF}_{\:\:\: t} \hspace{-0.1cm}\mathscr{D}_{T}^{q}{g_{1}(t)}= \frac{\partial \mathscr{H}}{\partial {f_{1}}}(t)=q_{1}f_{1}(t)-k_{3}g_{2}(t)-3k_{4}g_{2}(t)f_{1}^{2}(t), \\ {}^{CF}_{\:\:\: t} \hspace{-0.1cm}\mathscr{D}_{T}^{q}{g_{2}(t)}= \frac{\partial \mathscr{H}}{\partial {f_{2}}}(t)=q_{2}f_{2}(t)+g_{1}(t)-k_{1}g_{2}(t)-3k_{2}g_{2}(t)f_{2}^{2}(t), \end{cases}\displaystyle \end{aligned}$$
(39)
$$\begin{aligned}& \textstyle\begin{cases} \frac{\partial \mathscr{H}}{\partial {u_{1}}}(t)=r_{1}u_{1}(t)+g_{1}(t)=0 \quad \Rightarrow\quad u_{1}^{*}(t)=-\frac{1}{r_{1}}g_{1}(t), \\ \frac{\partial \mathscr{H}}{\partial {u_{2}}}(t)=r_{2}u_{2}(t)+g_{2}(t)=0 \quad \Rightarrow\quad u_{2}^{*}(t)=-\frac{1}{r_{2}}g_{2}(t), \end{cases}\displaystyle \end{aligned}$$
(40)

with initial conditions \(f_{1}(0)\) and \(f_{2}(0)\) and final conditions \(g_{1}(T)=g_{2}(T)=0\). To solve the optimality conditions (38)–(40), several methods have been introduced into the literature [39–41]; here we employ the quadratic numerical method proposed in Sect. 4 in combination with the forward–backward sweep iterative algorithm suggested in [42]. The convergence and stability of this scheme were proved in [43]. Considering \(F(t)=(f_{1}(t),f_{2}(t))\), \(G(t)=(g_{1}(t),g_{2}(t))\), and \(F_{i}(t)\), \(G_{i}(t)\) as the approximations of \(F(t)\), \(G(t)\) at the ith iteration, respectively, we summarize the algorithm as follows.

Algorithm

  • Step 1. Initiate the costate variable \(G_{0}(t)\) by making an initial guess and set \(J_{0}=\infty \).

  • Step 2. Take into account a positive integer M and set \(i=1\).

  • Step 3. Employ the value of \(G_{i-1}(t)\) and calculate \(F_{m+1}\), \(0\leq {m}\leq {M-1}\), forward in time by the numerical procedure designed in Sect. 4.

  • Step 4. Apply a linear interpolation and compute \(F_{i}(t)\) from \(F_{0},F_{1},\ldots ,F_{M}\).

  • Step 5. From Step 4 utilize the value of \(F_{i}(t)\) and obtain \(G_{M-m}\), \(1\leq {m}\leq {M}\), backward in time from the numerical procedure proposed in Sect. 4.

  • Step 6. Apply a linear interpolation and compute \(G_{i}(t)\) from \(G_{0},G_{1},\ldots ,G_{M}\).

  • Step 7. Stop the algorithm if the error \(e_{i}=\lvert {\frac{J_{i}-J_{i-1}}{J_{i}}}\rvert \) satisfies \(e_{i}<\epsilon \) in which \(\epsilon >0\) is a predefined error bound; otherwise, replace i by \(i+1\) and go to Step 3.

To illustrate the efficacy of the above iterative algorithm, we consider \(T=20\), \(q_{1}=q_{2}=1\), \(r_{1}=r_{2}=1\), and solve the FOCP (36)–(37) along with the state and costate optimality equations (38)–(40). Figures 7–9 display the simulation results, which verify that the chaotic behaviors of the fractional dynamical system (13) are compensated by the proposed fractional optimal controller, and the controlled system now exhibits stable periodic solutions.

Figure 7
figure 7

The state variable \(f_{1}(t)\) of the controlled system (35) for \(q=0.992,0.993,0.994,0.995,0.996,0.997, 0.998,0.999,1\)

Full size image
Figure 8
figure 8

The state variable \(f_{2}(t)\) of the controlled system (35) for \(q=0.992,0.993,0.994,0.995,0.996,0.997, 0.998,0.999,1\)

Full size image
Figure 9
figure 9

The phase-portraits of the controlled system (35) for \(q=0.992,0.993,0.994,0.995,0.996,0.997, 0.998,0.999,1\)

Full size image

8 Conclusions

In this paper, the complicated behaviors of a nonlinear suspension system were analyzed in the framework of fractional-order calculus. First, we considered the fractional-order derivative operator in the sense of CF with exponential kernel. Then we proposed a quadratic numerical method to solve the related FDEs. The system states and the phase-portraits shown in Figs. 1–3 demonstrated that both chaotic and nonchaotic behaviors can be drawn out by a fractional-order mathematical model. Indeed, if the fractional order q tends to 1, the fractional-order model reveals chaotic behaviors; however, as q becomes smaller, a twice-periodic behavior is emerged. In addition, a state-feedback controller and also a chaos optimal compensator were suggested, which suppressed the system chaotic motion. To solve the related FOCP, we used the proposed numerical method and designed an iterative forward–backward algorithm. Simulation results confirmed the effectiveness of the proposed FC approach in the modeling and control of nonlinear dynamical systems with both chaotic and nonchaotic oscillations.