1 Introduction

Before 1971, it was generally believed that basic circuit elements comprised three well-known components: resistance, inductance, and capacitance. However, after the publication of the missing circuit element by Chua^[1], the discovery of memristor changed people's traditional perception of electronic circuits. Memristive systems with integer-order have been well discussed, but fractional-order memristive systems have rarely been studied.

Fractional calculus plays a very important role in various fields of science and engineering. In 1974, Oldham & Spanier ^[2] gave it a complete mathematical presentation. A large number of studies show that the fractional-order model can depict the characteristics of the actual system more accurately. One year later^[3], described a fractional-order memristor-based classical Chua's oscillator, in which the memristor is a flux-controlled memristor. Since then, a lot of research has been focused on fractional-order memristor models and fractional-order memristor-based systems ^[4-8].

The nonlinear dynamic of the fractional-order delayed chaotic system has been analyzed ^[9-14]. In Ref.[13], the researcher analyzed the bifurcation oscillators with delay. In addition, the crucial focus chaos behaviors in the fractional time-delay chaotic system were also investigated^[15-16]. However, most of the current researches are limited to mathematical models. No circuit experiments have been conducted^[17-18]. In this paper, we propose a nonlinear time-delay chaotic circuit based on fractional-order memristive system. Results of the system circuit simulation accomplished by Multisim are consistent with theoretical analysis and numerical simulation. The structure of this article is as follows: In Section 2, we propose a fractional-order memristor and analyze the electrical characteristics of the fractional-order memristor. In Section 3, we propose a fractional-order memristive delayed chaotic circuit and present the phase portraits of the proposed fractional-order memristive delayed chaotic circuit. Section 4 provides a detailed analysis on the dynamic behaviors of the memristive circuit mentioned. Remarkably, the results of the system circuit simulation accomplished by Multisim are perfectly consistent with theoretical analyses and numerical simulation in Section 5. Finally, the conclusions are given in Section 6.

2 Modeling and Analysis of Fractional-Order Memristor 2.1 Fractional Calculus

Several existing concept of fractional calculus were proposed in Ref.[2]. The fundamental operator _aD_t^q is defined as follows:

$_aD_t^q = \left\{ {\begin{array}{*{20}{l}} {\frac{{{{\rm{d}}^q}}}{{{\rm{d}}{t^q}}}, }&{q > 0}\\ {1, }&{q = 0}\\ {\int_a^t {{{({\rm{d}}\tau )}^{ - q}}, } }&{q < 0} \end{array}} \right. $

(1)

where a, t are the upper and lower limits of the operation, respectively; q is the order of the calculus operator, whose value can be any complex number. If q > 0, _aD_t^q means fractional derivative; if q < 0, _aD_t^q is expressed as fractional integral ^[2].

In the development of fractional calculus theory, the initial conditions corresponding to the Riemann-Liouville derivatives have no physical meaning, while the initial conditions of the systems described by the Caputo derivatives have clear physical meanings and the forms are similar to the integer-order differential equations. The fractional derivative of f(t) is as follows:

$ _aD_t^q\cdot f\left( t \right) = \frac{1}{{\Gamma \left( {n - q} \right)}}\int_a^t {\frac{{{f^n}\left( \tau \right)}}{{{{\left( {t - \tau } \right)}^{q - n + 1}}}}{\rm{d}}\tau } $

(2)

where 0≤n-1≤q≤n, n=[q]+1, [q] is the largest integer not exceeding q, Γ(n) is the Euler Gamma function, x⁽ⁿ⁾(t) is the n-order derivative of x(t), Γ(·) is the Gamma function, which is a function of a variable x defined via a convergent improper integral

$ \Gamma \left( x \right) = \int_0^\infty {{t^{x - 1}}{e^{ - t}}{\rm{d}}t} $

(3)

n∈N is the first integer which is not less than α. It is a fractional derivative of a continuous function of time, namely f(t), which is given in terms of a time integral. Under natural conditions on the function f(t), the Laplace transform satisfies

$ \begin{array}{l} L\left\{ {_aD_t^q\cdot f\left( t \right)} \right\} = {s^q}F\left( s \right) - \sum\limits_{k = 0}^{n - 1} {{s^{q - k - 1}}{f^{(k)}}\left( 0 \right)} , \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {n - 1 < q \le n} \right) \end{array} $

(4)

2.2 Fractional-Order Flux-Controlled Memristor

Muthuswamy et al.^[19] studied the flux-controlled memristor model equivalent realization circuit shown in Fig. 1. It consists of three remote arithmetic amplifiers, five resistors, two multipliers, and one capacitor. The first stage op amp U₁ circuit can avoid the load effect. The second stage op amp U₂ circuit connected to resistor R₀ and capacitor C₀ is an integrator.

The mathematical model describing the memristive behavior is established as follows:

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}{x_2}}}{{{\rm{d}}t}} = \frac{1}{{{R_0}{C_0}}}{x_1}\left( t \right)\\ i\left( t \right) = f\left( {{x_1}, {x_2}} \right) = \left( {\frac{1}{{{R_1}}} - \frac{1}{{{R_2}}} + \frac{{{g_1}{g_2}}}{{{R_2}}}{{\left[ {{x_2}\left( t \right)} \right]}^2}} \right){x_1}\left( t \right) \end{array} \right. $

(5)

where x₁ is the input of the memristor, x₂ is internal state of the memristor, f(x₁, x₂) is the output of the memristor. g₁ is the variable scale factor in multiplier M₁, g₂ is the variable scaling factor in multiplier M₂. Assume α=R₁^-1-R₂^-1, β=g₁g₂/R₂ in the active magnetron memristive equivalent circuit described by the quadratic nonlinear characteristic curve in Fig. 1. Thus, we set R₀=8.2 k Ω, C₀=47 nF, R₁=750 Ω, R₂=1.5 k Ω, R₃=R₄=8.2 k Ω, g₁=0.1, g₂=1.3. The equivalent parameters for the proposed active magnetron memristor can be calculated as α=0.666 7 ms and β=0.008 7 S/Wb.

In practical calculations, there are a variety of methods to solve fractional calculus. By obtaining the complex frequency domain transfer function 1/s^q, the expanded form of the complex frequency domain is obtained and the complex frequency domain form is converted into the time domain form for numerical solution. Ref.[20] introduced a Bode graph approximation method. We determined the expansion of F=1/s^q. Ref.[21] derived the 1/s^q expansion of q from 0.1 to 0.9. In order to implement the equivalent circuit of fractional-order memristor shown in Eq. (8), we first considered applying the approximation technique proposed by Oustaloup to get the approximate transfer function of the fractional-order integrator operator ^[22].

We designed a unit circuit that implements a 1/s^q expansion of q from 0.1 to 0.9 based on the 1/s^q expansion of the approximate error of 2 dB, as shown in Fig. 3. The complex frequency domain expression of the equivalent circuit is as follows:

$ H\left( s \right) = {C^{ - 1}}{s^{ - q}} = \sum\limits_{k - 1}^n {\frac{{C_k^{ - 1}}}{{\left( {s + {{({R_k}{C_k})}^{ - 1}}} \right)}}} $

(6)

where C is the unit parameter. Let C=1 μF, H(s)=1/s^q.

Table 1 gives the approximate transfer functions of fractional integrator F(s) when fractional order q is selected as 0.6, 0.7, 0.8, 0.9. The corresponding resistance and capacitance values are shown in Table 2.

Then, we used electronic circuits to achieve the fractional-order flux-controlled memristor in the complex frequency domain as shown in Fig. 2, with a unit circuit (Fig. 3) for implementing the approximate transfer functions of fractional-order integrator operator depicted in Table 2.

According to the concept of fractional calculus, the real capacitor can be extended to fractional-order to describe the real electricity characteristics ^[23]. In this paper, the integer-order capacitor C in the memristor is actually fractional-order capacitor in reality. The fractional-order unified model of memristor is derived as follows:

$\begin{array}{l} \frac{{{\rm{d}}x_2^q}}{{{\rm{d}}t}} = \frac{1}{{{R_0}C_0^q}}{x_1}\\ f({x_{1, }}{x_2}) = \left( {\frac{1}{{{R_1}}} - \frac{1}{{{R_2}}} + \frac{{{g_1}{g_2}}}{{{R_2}}}{{[{x_2}]}^2}} \right){x_1} \end{array} $

(7)

where q is the order of the fractional-order.

To prove the validity of the theory mentioned above, circuit simulations were carried out to realize the memristor. The driven voltage was chosen as x₁(t)=X₁sin(2πft)V. When X₁=4 V, the current-voltage curve is described in Fig. 4.

Fig. 4(a) exhibits the pinched hysteresis loop when f=100 Hz and the fractional q is 0.6, 0.8, and 0.9, respectively. Fig. 4(b) exhibits the pinched hysteresis loop when q=0.9 and f equals to 100, 200, and 500 Hz, respectively. Through simulation experiments, we concluded that in the same frequency, as the order q increases, the area of the pinched hysteresis loop increases monotonically; in the same order, as the frequency f increases, the area of the pinched hysteresis loop decreases monotonically.

2.3 Delay Circuit

Ref.[24] compared the APF with the LCL, Bessel delay circuit. Other delay circuits have their own shortcomings. APF is the most efficient compared with other delay blocks. Therefore, APF experimental system was chosen in this paper.

This delay unit is included by a series of All Pass Filter (APF) as shown in Fig. 5 ^[25]. The same delay unit is used in Refs.[26-27]. Variable delay was studied as delay blocks. The transfer function of the APF is as follows:

$ T(s) = - {a_1}\frac{{s - {w_0}}}{{s + {w_0}}} $

(8)

The flat gain q₁=1 is determined by R₈ and R_g. TD≈RC is the time delay provided by each APF block. Thus this delay unit needs i sections to realize a delay TD≈iRC, (i=1, 2, 3).

Because TD changes with the resistance change of the sliding rheostat, in order to avoid emergency, we used two time-delay modules in series as shown in Fig. 5.

3 Fractional-Order Memristive Time-Delay Chaotic Circuit

In this module, we present a circuital realization of the fractional-order memristive delayed Chaotic system (9) to prove the feasibility and correctness of the theoretical model. The designed circuit is shown in Fig. 6. Fractional-order memristive device and delayed unit was analyzed in Ref.[28]. The equation of the circuit in Fig. 6 is derived as follows:

$ \left\{ \begin{array}{l} \frac{{{\rm{d}}x_2^q}}{{{\rm{d}}{t^q}}} = l{x_1}\\ \frac{{{{\rm{d}}^q}{x_1}}}{{{\rm{d}}{t^q}}} = a{x_1} + b{x_{1\tau }} + c{x_1}{\left( {{x_2}} \right)^2} \end{array} \right. $

(9)

where a=-α/C₁^q, b=-(RC₁^q)^-1, c=-β/C₁^q, l=(R₀C₀^q)^-1, α=0.666 7 ms, β=0.008 7 S/Wb, and τ is the time-delay.

The fractional-order q is set to be 0.9, and initial values are set to be x₁=0.1 V, x₂=0.1 V. In this paper, we studied delayed memristive chaotic system in Fig. 6 ^[29].

Subsequently, we performed numerical simulation on the decomposed fractional mathematical model. The simulation was carried out by MATLAB, via utilizing the predictor-corrector method to simulate the fractional differential Eq. (9).Other parameters include fractional order q=0.9, l=2.5, b=-2, τ=1 and initial values (0.1, 0.1). In order to easily see the influence of the parameter a on the time-delayed chaotic circuit of the fractional memristor, we drew the two-dimensional phase diagrams and time-series plots of different parameter a, as is shown in Fig. 7. The phase portrait is a stable limit cycle with parameter a=1.62 as shown in Fig. 7(a). The corresponding time-series plots is shown in Fig. 7(b). When parameter a > 1.64, the system goes into a chaotic state. When parameter a > 1.64, the real time and phase plane plots are shown in Fig. 7(c) and Fig. 7(d), respectively. When parameter a≥1.7, the system enters double-scroll chaotic attractor in Fig. 7(e)-(f).

4 Analysis of Dynamical Behaviors 4.1 Stability of Equilibrium Points

We presume that the system (9) has an equilibrium point X^*=(x₁^*+x₂^*). Then, the linearized system near the equilibrium point X^* is as follows:

$ D_t^q\mathit{\boldsymbol{X}} = \mathit{\boldsymbol{JX}} + {\mathit{\boldsymbol{J}}_\tau }{\mathit{\boldsymbol{X}}_\tau } $

(10)

where X and X_τ are the column vectors (x₁, x₂) at t and t-τ,

$ \mathit{\boldsymbol{J}} = \left[ {\begin{array}{*{20}{c}} {a + c{{(x_2^*)}^2}}&{2cx_1^*x_2^*}\\ l&0 \end{array}} \right] $

(11)

and

$ {\mathit{\boldsymbol{J}}_\tau } = \left[ {\begin{array}{*{20}{c}} b&0\\ 0&0 \end{array}} \right] $

(12)

Theorem 1  We can get all the roots of the characteristic equation^[30]

$ {\rm{det}}\left( {{s^q}\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{J}} - {\rm{exp}}\left( { - s\tau } \right){\mathit{\boldsymbol{J}}_\tau }} \right) = 0 $

(13)

where I is the identity matrix.

The necessary factor for the system (9) to generate chaos is that Eq.(13) has the root of the real part.

Thus we have the following equation:

$ {\rm{det}}\left( {\begin{array}{*{20}{c}} {{s^q} - \left( {a + c{{\left( {x_2^*} \right)}^2} + b{e^{ - s\tau }}} \right)}&{2cx_1^*x_2^*}\\ { - l}&{{s^q}} \end{array}} \right) = 0 $

(14)

Based on Eq.(14), the characteristic equation is as below:

$ {s^{2q}} + {A_1}{s^q} + {A_2}{s^q}{e^{ - s\tau }} + {A_3} = 0 $

(15)

where

$ \left\{ \begin{array}{l} {A_1} = - \left( {a + c{{\left( {x_2^*} \right)}^2}} \right)\\ {A_2} = - b\\ {A_3} = - 2lcx_1^*x_2^* \end{array} \right. $

(16)

Assuming that a root of Eq.(15) is s=iw= |w|(cos(π/2)+isin(±(π/2))), and w is a real number, when w > 0, we obtain isin(π/2); when w < 0, we obtain -isin(π/2). Real and imaginary parts are separated, then we get

$ \left\{ \begin{array}{l} {\omega ^{2q}}{\rm{cos}}\;q{\rm{ \mathsf{ π} }} + {A_1}{\omega ^q}{\rm{cos}}\frac{{q{\rm{ \mathsf{ π} }}}}{2} + {A_3} + \\ \;\;\;\;\;\;{A_2}{\omega ^q}{\rm{cos}}\frac{{q{\rm{ \mathsf{ π} }}}}{2}{\rm{cos}}\;\omega \tau + {A_2}{\omega ^q}{\rm{sin}}\frac{{q{\rm{ \mathsf{ π} }}}}{2}{\rm{sin}}\;\omega \tau = 0\\ {\omega ^{2q}}{\rm{sin}}\;q{\rm{ \mathsf{ π} }} + {A_1}{\omega ^q}{\rm{sin}}\frac{{q{\rm{ \mathsf{ π} }}}}{2} - \\ \;\;\;\;\;\;\;{A_2}{\omega ^q}{\rm{cos}}\;\frac{{q{\rm{ \mathsf{ π} }}}}{2}{\rm{sin}}\;\omega \tau + {A_2}{\omega ^q}{\rm{sin}}\frac{{q{\rm{ \mathsf{ π} }}}}{2}{\rm{cos}}\;\omega \tau = 0 \end{array} \right. $

(17)

We assume

$ \left\{ \begin{array}{l} {B_1} + C{\rm{cos}}\;\omega \tau + D{\rm{sin}}\;\omega \tau = 0\\ {B_2} + C{\rm{sin}}\;\omega \tau + D{\rm{sin}}\;\omega \tau = 0 \end{array} \right. $

(18)

and

$ \left\{ \begin{array}{l} {B_1} = {\omega ^{2q}}{\rm{cos}}\;q{\rm{ \mathsf{ π} }} + {A_1}{\omega ^q}{\rm{cos}}\frac{{q{\rm{ \mathsf{ π} }}}}{2} + {A_3}\\ {B_2} = {\omega ^{2q}}{\rm{sin}}\;q{\rm{ \mathsf{ π} }} + {A_1}{\omega ^q}{\rm{sin}}\frac{{q{\rm{ \mathsf{ π} }}}}{2}\\ C = {A_2}{\omega ^q}{\rm{cos}}\frac{{q{\rm{ \mathsf{ π} }}}}{2}\\ D = {A_2}{\omega ^q}{\rm{sin}}\frac{{q{\rm{ \mathsf{ π} }}}}{2} \end{array} \right. $

(19)

We set

$ B_1^2 + B_2^2 = {C^2} + {D^2} $

(20)

Then, we get

$ {\left| \omega \right|^{4q}} + {\psi _1}{\left| \omega \right|^{3q}} + {\psi _2}{\left| \omega \right|^{2q}} + {\psi _3}{\left| \omega \right|^q} + {\psi _4} = 0 $

(21)

where

$ \left\{ \begin{array}{l} {\psi _1} = 2{A_1}\left( {{\rm{cos}}\;\frac{{q{\rm{ \mathsf{ π} }}}}{2}{\rm{cos}}\;q{\rm{ \mathsf{ π} }} + {\rm{sin}}\;q{\rm{ \mathsf{ π} }}\;{\rm{sin}}\;\frac{{q{\rm{ \mathsf{ π} }}}}{2}} \right)\\ {\psi _2} = 2\left( {{A_3}{\rm{cos}}\;q{\rm{ \mathsf{ π} }} + A_1^2 - A_2^2} \right)\\ {\psi _3} = 2{\rm{cos}}\frac{{q{\rm{ \mathsf{ π} }}}}{2}{A_1}{A_3}\\ {\psi _4} = A_3^2 \end{array} \right. $

(22)

We can calculate the root of Eq.(21).

We set

$ y\left( \omega \right) = {\left| \omega \right|^{4q}} + {\psi _1}{\left| \omega \right|^{3q}} + {\psi _2}{\left| \omega \right|^{2q}} + {\psi _3}{\left| \omega \right|^q} + {\psi _4} $

(23)

Lemma 1 For Eq.(15), the following results hold:

(a) If Ψ_k > 0(k=1, 2, 3, 4) and A₃+A₄0, Eq.(15) has no root with zero real parts for all τ≥0.

(b) If Ψ_k > 0(k=1, 2, 3, 4) and Ψ_k < 0, Eq.(15) has a pair of purely imaginary roots ±w when τ=τ_j, j=0, 1, 2, …, where

$ \tau _j^{(k)} = \frac{1}{{{\omega _k}}}\left[ {{\rm{arccos}}\left( { - \frac{{{B_1}C + {B_2}D}}{{{C^2} + {D^2}}}} \right) + 2j{\rm{ \mathsf{ π} }}} \right] $

(24)

where j=0, 1, 2, …; k=1, 2, 3, …, and w₀ is the unique positive zero of the function y(w).

Proof

(a) From Ψ_K > 0(k=1, 2, 3, 4), we obtain y=(0)=Ψ₄ > 0,

$ \begin{array}{l} {y^\prime }\left( \omega \right) = 4q{\left| \omega \right|^{4q - 1}} + 3{\psi _1}q{\left| \omega \right|^{3q - 1}} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;2{\psi _2}q{\left| \omega \right|^{2q - 1}} + {\psi _3}q{\left| \omega \right|^{q - 1}} \end{array} $

(25)

Combine q > 0 and Ψ_k > 0(k=1, 2, 3, 4) we get that Eq.(15) has no real root. Hence, Eq.(15) has no purely imaginary root, provided that A₃+A₄≠0, Ψ₄ < 0 is not a root of Eq.(15).

(b) If Ψ₄ < 0, it is simple to get y(0)=Ψ₄ < 0. Then, by $ \mathop {{\rm{lim}}}\limits_{w \to \infty } y\left( w \right) = + \infty $ and y′(w) > 0 for w > 0, there exists a unique positive number w₀ such that s(τ)=ξ(τ)+iw(τ), where ω₀ is a root of Eq.(21), τ_j is defined in Eq.(18), and s(τ)=ξ(τ)+iw(τ) is a root of Eq.(18). It can be proved that ±w₀ is a pair of purely imaginary roots of Eq.(18) when τ=τ_j, j=0, 1, 2, …

4.2 The Largest Lyapunov Exponent and Bifurcation Diagram

As the system parameter a changes, balance point of the fractional-order delayed chaotic system, stability of the system, and state of the system will alter. The largest Lyapunov exponent and bifurcation diagram of the fractional-order delayed chaotic system are considered in order to demonstrate the chaos. Under the initial condition (0.1, 0.1), when a changes in the range [1.35, 2], the maximum Lyapunov exponent and bifurcation diagram of the fractional-order delayed chaotic system is shown in Fig. 8. For better observation, the maximum Lyapunov exponent of the fractional-order delayed chaotic system is basically consistent with the bifurcation diagram.

It is obvious that the orbits of the memristive circuit start from the periodic state, then when a passes through the threshold, Hopf bifurcation appears, and finally turn into chaotic state with the increase of a. When parameter a increases to 1.62, the phase portrait of the system (9) turns to be a limit cycle. Besides, with decrease of parameter a, the orbit degenerates into a single-scroll attractor as shown in Fig. 7(c).

4.3 Coexistence Attractors and Coexisting Bifurcation Mode

When we set q=0.9, l=2.5, b=02, τ=1, parameter a changes in the range [1.35, 2], and the state initial values are (0.1, 0.1) and (-0.1, -0.1) by blue and red, respectively, it is clear from Fig. 9 that there is a coexistence bifurcation mode phenomenon. In the case of different initial values, the system's bifurcation diagram is inconsistent.

According to dynamics given in Fig. 10, when parameter a > 1.64, the system (9) manifests coexisting attractors. When parameter a=1.68, the state initial values are (0.1, 0.1) and (-0.1, 0.1) by blue and red, respectively. The coexistence of single-scroll chaotic attractors are shown in Fig. 10(a), and the state initial values are (0.1, 0.1) and (-0.1, -0.1) by blue and red, respectively, as is shown in Fig. 10(b). When parameter a=1.73, the state initial values are (0.1, 0.1) and (-0.1, -0.1). The coexistence of double-scroll chaotic attractors are expressed in Fig. 10(c). From this phenomenon, we know that the fractional-order memristive time-delay chaotic system has wealth dynamic characteristics.

5 Experimental Research on Circuit

To confirm the theoretical investigations, circuit simulations are carried out by Multisim. The structure of the circuit was designed as Fig. 11 when q is 0.9. The proposed circuit was designed on Multisim software using op-amp AD711KN and multiplier AD633JN with ±15 volts and R=11.24 k Ω. The phase plane and time-series plots for variable a is as shown in Fig. 12. The phase portrait is a stable limit cycle when a=1.62, a single-scroll chaotic attractor when a=1.68, and a double-scroll attractor when a=1.73. Simulation results demonstrate the influence of parameter a on chaotic behaviors.

6 Conclusions

In this paper, we presented a fractional-order memristive delayed chaotic system. Some dynamical properties, including equilibrium points, stability, bifurcation, and Lyapunov exponent of the oscillator, were investigated in detail by theoretical analyses and simulations. Moreover, we also found the nonlinear phenomena of coexisting bifurcation and coexisting attractor. The phenomenon shows that the state of this oscilator was highly sensitive to its initial value, which is called coexistent oscillation in this paper. Besides, the results of the circuit simulation accomplished by Multisim were perfectly consistent with theoretical analyses.