# Dynamical Behavior of Fractional-Order Delayed Feedback Control on the Mathieu Equation by Incremental Harmonic Balance Method

Authors: Shaofang Wen; Yongjun Shen; Jiangchuan Niu; Yunfei Liu
Journal: Shock and Vibration
 (2022)

---

## Abstract

In this study, the dynamical analysis of the Mathieu equation with multifrequency excitation under fractional-order delayed feedback control is investigated by the incremental harmonic balance method (IHBM). IHBM is applied to the fractional-order delayed feedback control system, and the general formulas of the first-order approximate periodic solution for the Mathieu equation are derived. Caputo’s definition is adopted to process the fractional-order delayed feedback term. The general formulas of this system are suitable for not only the weakly but also the strongly nonlinear fractional-order system. Through the analysis of the general formulas of this system, it shows that fractional-order delayed feedback control has two functions, which are velocity delayed feedback control and displacement delayed feedback control. Next, the numerical simulation of the system is carried out. The comparison between the approximate analytical solution and the numerical iterative result is made, and the accuracy of the approximate analytical result by IHBM is proved to be high. At last, the effects of the time delay, feedback coefficient, and fractional order are investigated, respectively. It is generally known that time delay is common and inevitable in the control system. But the fractional order can be used to adjust the influence caused by time delay in fractional-order delayed feedback control. Those new system characteristics will provide theoretical guidance to the design and the control of this kind system.


---

## Body

## 1. Introduction

The definition of fractional calculus has existed for more than 300 years. Fractional calculus has many advantages. For example, it could better reflect the viscoelastic and the memory properties of materials, and also provide better control performance. With the development of science and the increase in the demand of complex engineering application, a lot of scholars in different research fields are attracted by fractional calculus [1–5].At present, fractional calculus has played a very important role in different fields, such as electrochemistry, signal processing, bioengineering, mechanics, and control. The research of basic theory has been greatly improved on fractional calculus, and the achievement of the basic theoretical research mainly focuses on the following aspects: qualitative analysis, numerical calculation, and analytical research [6–10]. For example, Wang et al. [11, 12] studied the dynamics of the linear fractional-order system under external excitation and proved that the fractional-order feedback can not only adjust the damping force but also adjust the elastic force. Then, the stabilities of the various fractional-order systems are studied, and an effective judging method for the stability of the fractional-order system is proposed. Trigeassou et al. [13] investigated the stability of the fractional-order system by Lyapunov theory. Li et al. [14, 15] presented some higher-order iterative algorithms of fractional-order differential equations based on the Simpson method and put forward some numerical calculation methods of fractional-order partial differential equations. Ebaid et al. [16] developed a new approach to solve a class of first-order fractional initial value problems based on Riemann–Liouville fractional derivative. Failla et al. [17] proposed a numerical method to investigate the dynamical response of the fractional-order system with random excitation. Shen et al. [18, 19] studied the approximate analytical solutions of several fractional-order systems by the improved average method and the dynamics of the strongly nonlinear fractional-order system by the improved incremental harmonic balance method.In the engineering application, the system model is generally complex, and it is often affected by multifrequency excitation and strong nonlinearity. At this moment, the system cannot be regarded as a perturbed system for its derived linear system because the strong nonlinearity will seriously affect the stability of the periodic solution. Those classical perturbation methods for the weakly nonlinear system are used to study the strongly nonlinear problems will cause enormous error [20, 21].In the present survey, some efficient analytical methods and numerical integration methods for the strongly nonlinear system have been presented now. For the numerical integration method, because of its slow convergence speed, it is very difficult to analyze the effects of those parameters in complex dynamical systems. However, some efficient analytical methods, such as the generalized multiple-scale method with its improved version, the perturbation-incremental method and the homotopy analysis method [22–25] could provide satisfactory results for those strongly nonlinear oscillators in some complex engineering. IHBM is a classical, effective, and semianalytical method with many advantages. For example, it could be used to deal with strongly and weakly nonlinear systems at the same time [26–28], and its convergence speed is very fast. Moreover, the solution with higher accuracy can be obtained through the increasing of the solution order based on IHBM. That means that IHBM can be applied into those different nonlinear systems and also is widely used in the engineering field.The application of fractional calculus in the control system is also unique. Chen et al. [29] proposed a new fractional-order control method for the four-wheel steering system, which improved the transient response of the vehicle in the steering process. Chen et al. [30] studied the nonlinear dynamical characteristics of the fractional-order Duffing system with random excitation, and presented an efficient bifurcation control method based on the fractional-order PIλDμ feedback controller. Avci et al. [31] proposed spectral formulation for a fractional optimal control problem defined in spherical coordinates. The dynamical response of the Mathieu–Duffing oscillator with fractional-order delayed feedback control is studied by the average method, and it can be concluded that fractional-order delayed feedback control can better control the dynamical characteristics of the system [19]. As we all know, time delay is inevitable in control engineering. The dynamical phenomena of the system under the combined action of parametric excitation and multifrequency excitation are more complicated [32–34]. Accordingly, the general formulas of the Mathieu equation with multifrequency excitation under fractional-order delayed feedback control based on IHBM are obtained in the following. The dynamical analysis of the Mathieu equation with multifrequency excitation under fractional-order delayed feedback control is investigated by the obtained general formulas. Also, the general formulas can be used to study a variety of complex strongly nonlinear systems.In this paper, the Mathieu equation with multifrequency excitation is presented in Section2. In addition, the general formulas of the Mathieu equation with multifrequency excitation under fractional-order delayed feedback control are derived by IHBM in Section 2. Then, in Section 3, the numerical result is studied and the correctness of the result by IHBM is verified by numerical simulation. In Section 4, the influence of the fractional order, the feedback gain coefficient, and time delay is investigated, respectively. From the analysis of the general formulas and the simulation results, the correlation between the fractional order and time delay is studied. Fractional order and time delay can be adjusted with each other, so fractional-order delayed feedback control will provide a wider design and adjustment space in the control system. At last, the main conclusions are drawn.


## 2. The General Formulas by IHBM

In this paper, the Mathieu equation with multifrequency excitation is considered. At first, this system under fractional-order delayed feedback control is(1)d2xtdt2+μdxtdt+1+2εcosωtxt+Ut=F0−∑k=1LFkcoskωt,where xt is the system displacement, μ is the damping coefficient, 2εcosωt is the periodic time-varying stiffing coefficient, and F0 and ∑k=1LFkcoskωt are the constant excitation and multifrequency excitation, respectively. Ut is fractional-order delayed feedback control. Here, Ut=FpDpxt−τ1, τ1 is the time delay, Fp is the feedback gain coefficient, and Dpxt−τ1 represents p-order (0≤p≤2) derivative of xt−τ1 to t. When p  = 0, fractional-order delayed feedback becomes displacement delayed feedback. When p  = 1, fractional-order delayed feedback becomes velocity delayed feedback.At present, there are several most frequently used definitions for fractional calculus, such as Grünwald–Letnikov, Riemann–Liouville, and Caputo definitions, which are equivalent under some conditions [1, 2]. Here, Caputo’s definition is adopted as(2)Dpxt=1Γ1−p∫0tx′ut−updu,where Γy is Gamma function satisfying Γy+1=yΓy.Letτ=ωt, τ0=ωτ1, and equation (1) becomes(3)ω2x¨τ+ωμx˙τ+1+2εcosτxτ+FpωpDτpxτ−τ0=F0+∑k=1LFkcoskτ.The symbol “⋅” represents the derivative with respect to τ. The periodic solution with N-order harmonic terms of equation (3) is expressed as follows:(4a)x0=a0+∑n=1Nancosnτ+bnsinnτ,and(4b)Δx0=Δa0+∑n=1NΔancosnτ+Δbnsinnτ.Substitutingx=x0+Δx0 into equation (3), using Taylor series expression, and omitting the higher-order terms of the small increment Δx, one can obtain the differential equation as(5)ω2Δx¨0+ωμΔx˙0+1+2εcosτΔx0+FpωpDτpΔx0τ−τ0=F0−∑k=1LFkcoskτ−FpωpDτpΔx0τ−τ0−f0,where f0=ω2x¨0+ωμx˙0+1+2εcosτx0.Defining the following vectors(6a)X=1,cosτ,cos2τ,…cosnτ,sinτ,sin2τ,…sinnτ,(6b)Xτ−τ0=1,cosτ−τ0,cos2τ−τ0,…cosnτ−τ0,sinτ−τ0,sin2τ−τ0,…sinnτ−τ0,(6c)A0=a0,a1,…an,b1,b2,…bnT,(6d)ΔA0=Δa0,Δa1,…Δan,Δb1,Δb2,…ΔbnT,one can get(6e)x0=XA0,Δx0=XΔA0,x0τ−τ0=Xτ−τ0A0,Δx0τ−τ0=Xτ−τ0ΔA0.Based on Galerkin’s procedure [26, 27], the integral for equation (5) is deduced in the following. In this procedure, fractional-order delayed feedback item is an aperiodic function, so the time terminal is selected as T=∞. The other items of equation (5) are periodic functions with a period of 2π, so the time terminal is selected as T=2π. Accordingly, one can get the following equations:(7)12π∫02πδΔx0ω2Δx¨0+ωμΔx˙0+1+2εcosτΔx0dτ+1T∫0TδΔx0FpωpDτpΔx0τ−τ0dτ=12π∫02πδΔx0F0−∑k=1LFkcoskτ−f0dτ−1T∫0TδΔx0FpωpDτpx0τ−τ0dτ.Rearrange equation (7), and it yields(8)δΔA0T12π∫02πXTω2X¨+ωμX˙+1+2εcosτXdτΔA0+δΔA0T1T∫0TXTFpωpDτpXτ−τ0dτΔA0=δΔA0T12π∫02πXTF0−∑k=1LFkcoskτ−f0dτ−δΔA0T1T∫0TXTFpωpDτpx0τ−τ0dτ.Then, obtaining the 2N + 1 linearized equations about ΔA0(9)MΔA0=R,where(10a)M=M1+M2p,(10b)M1=12π∫02πXTω2X¨+ωμX˙+1+2εcosτXdτ,(10c)M2p=1T∫0TXTFpωpDτpXτ−τ0dτ,(10d)R=R1+R2p,(10e)R1=12π∫02πXTF0−∑k=1LFkcoskτ−f0dτ,(10f)R2p=1T∫0TXTFpωpDτpx0τ−τ0dτ,and expanding to matrix forms(11)M1=M11M12M21M22,R1=R10R11R12,one could obtain the explicit expression ofM1 and R1 as follows:(12a)M11=δij+εsin1+i−jθ1+i−jθ+sin1−i+jθ1−i+jθ,i=0,j=0,1,…,N,j=0,i=1,…,N,−12δijj2ω2+δij2+ε2sin1+i−jθ1+i−jθ+sin1−i+jθ1−i+jθ,i=1,…,N,j=1,…,N,(12b)M12=12μδijjω,i=0,1,…,N,j=1,…,N,(12c)M21=−12μδijjω,i=1,…,N,j=0,1,…,N,(12d)M22=−12δijj2ω2+δij2+ε2sin1+i−jθ1+i−jθ+sin1−i+jθ1−i+jθ,i=1,…,N,j=1,…,N,(12e)R10=F0−a0+R10iCJ,i=0,(12f)R11=12fi+i2ω2ai−μiωbi−ai2+R11iCJ,i≤L,12i2ω2ai−μiωbi−ai2+R11iCJ,i>L,,i=1,…,N,(12g)R12=12i2ω2bi+μiωai−bi2+R12iCJ,i=1,…,N,where(13)R10iCJ=R11iCJ=−∑j=0Nε2ajsin1+i+jθ1+i+jθ+sin1+i−jθ1+i−jθ+sin1−i+jθ1−i+jθ+sin1−i−jθ1−i−jθ,R12iCJ=−∑j=1Nε2bjsin1+i+jθ1+i+jθ+sin1+i−jθ1+i−jθ+sin1−i+jθ1−i+jθ+sin1−i−jθ1−i−jθ,where δij is the Kronecker symbol, and the value of θ is 2π.Then, the integration process of fractional-order delayed feedback is presented.M2p and R2p are the Jacobi matrix and the corrective vector, respectively, and the forms are as follows:(14)M2p=M11pM12pM21pM22p,R2p=R200pR21ipR22ip,where(15a)M11ijp=limT⟶∞1T∫0TFpωpcosiτDτpcosjτ−τ0dτ,i=0,1,…,N,j=0,1,…,N,(15b)M12ijp=limT⟶∞1T∫0TFpωpcosiτDτpsinjτ−τ0dτ,i=0,1,…,N,j=1,…,N,(15c)M21ijp=limT⟶∞1T∫0TFpωpsiniτDτpcosjτ−τ0dτ,i=1,…,N,j=0,1,…,N,(15d)M22ijp=limT⟶∞1T∫0TFpωpsiniτDτpsinjτ−τ0dτ,i=1,…,N,j=1,…,N,(15e)R200p=limT⟶∞−1T∫0TFpωpDτpx0τ−τ0dτ,(15f)R21ip=limT⟶∞−1T∫0TFpωpDτpx0τ−τ0cosiτdτ,i=1,…,N,(15g)R22ip=limT⟶∞−1T∫0TFpωpDτpx0τ−τ0siniτdτ,i=1,…,N.According to equations (2) and (15a), it can be written as(16)M11ijp=limT⟶∞1T∫0TFpωpcosiτDτpcosjτ−τ0dτ=limT⟶∞1T∫0TFpωpcosiτ1Γ1−p∫0τ−jsinju−τ0τ−updudτ.Lettings=τ−u and du=−ds, and substituting into equation (16), one can get(17)M11ijp=−jFpωpΓ1−plimT⟶∞1T∫0Tcosiτ∫0τsinjτ−τ0−jsspdsdτ=−jFpωpΓ1−plimT⟶∞1T∫0Tcosiτsinjτ−τ0∫0τcosjsspdsdτ−−jFpωpΓ1−plimT⟶∞1T∫0Tcosiτcosjτ−τ0∫0τsinjsspdsdτ.From equation (17), it could be found the expression for equation (17) has two integrals. Defining the first integral as A1 and the second integral as A2, and integrating it by parts, respectively, one have(18a)A1=−jFpωpΓ1−plimT⟶∞1T∫0T12siniτ+jτ−τ0−siniτ−jτ+τ0∫0τcosjsspdsdτ=jFpωp2Γ1−plimT⟶∞1Tcosiτ+jτ−τ0i+j−cosiτ−jτ+τ0i−j∫0τcosjsspds|0T−jFpωp2Γ1−plimT⟶∞1T∫0Tcosiτ+jτ−τ0i+j−cosiτ−jτ+τ0τi−jcosjττpdτ,(18b)A2=−−jFpωpΓ1−plimT⟶∞1T∫0T12cosiτ+jτ−τ0+cosiτ−jτ+τ0∫0τsinjsspdsdτ=jFpωp2Γ1−plimT⟶∞1Tsiniτ+jτ−τ0i+j+siniτ−jτ+τ0i−j∫0τsinjsspds|0T−jFpωp2Γ1−plimT⟶∞1T∫0Tsiniτ+jτ−τ0i+j+siniτ−jτ+τ0i−jsinjττpdτ.To simplify the process, we define the first part of equation (18a) as A11, the second part of equation (18a) as A12, the first part of equation (18b) as A21, and the second part of equation (18b) as A22. Here, two basic formulas [17] are introduced(19a)limT⟶∞∫0Tsinjttpdt=jp−1Γ1−pcospπ2,(19b)limT⟶∞∫0Tcosjttpdt=jp−1Γ1−psinpπ2.According to equation (19b), we can get(20a)A11=jFpωp2Γ1−plimT⟶∞1Tcosiτ+jτ−τ0i+j−cosiτ−jτ+τ0i−j∫0τcosjsspds|0T=Fpωpjpsinpπ/22limT⟶∞1TcosiT+jT−τ0i+j−cosiT−jT+τ0i−j=Fpωpjpsinpπ/22limT⟶∞1Tcosi+jTcosτ0i+j+sini+jTsinτ0i+j−cosi−jTcosτ0i+j+sini−jTsinτ0i−j=0,i≠j,Fpωpjpsinpπ/2sinτ02,i=j,(20b)A12=−jFpωp2Γ1−plimT⟶∞1T∫0Tcosiτ+jτ−τ0i+j−cosiτ−jτ+τ0τi−jcosjττpdτ=0.According to equation (19a), we can get(21a)A21=jFpωp2Γ1−plimT⟶∞1Tsiniτ+jτ−τ0i+j+siniτ−jτ+τ0i−j∫0τsinjsspds|0T=Fpωpjpcospπ/22limT⟶∞1TsiniT+jT−τ0i+j+siniT−jT+τ0i−j=Fpωpjpcospπ/22limT⟶∞1Tsini+jTcosτ0i+j−cosi+jTsinτ0i+j+sini−jTcosτ0i−j+cosi−jTsinτ0i−j=0,i≠j,Fpωpjpcospπ/2cosτ02,i=j,(21b)A22=−jFpωp2Γ1−plimT⟶∞1T∫0Tsiniτ+jτ−τ0i+j+siniτ−jτ+τ0i−jsinjττpdτ=0.Combining equations (17), (18a), (18b), (20a), (20b), (21a), and (21b) we obtain the explicit expression of M11P(22)M11ijp=0,i≠j,Fpωpjpsinpπ/2sinτ0+cospπ/2cosτ02,i=j.Substitute the system original parameters, and equation (22) becomes(23)M11ijp=0,i≠j,Fpωpjpcospπ/2−ωτ12,i=j.Similarly, the explicit expressions of the other matrices in equation (14) also could be derived. M2p and R2p are as follows:(24)M11ijp=δijFpωpip2cospπ2−ωτ1,i=0,1,…,N,j=0,1,…,N,M12ijp=δijFpωpip2sinpπ2−ωτ1,i=0,1,…,N,j=1,…,N,M22ijp=δijFpωpip2cospπ2−ωτ1,i=1,…,N,j=1,…,N,R200p=0,R21ip=−Fpωpaiip2cospπ2−ωτ1+biip2sinpπ2−ωτ1,i=1,…,N,R22ip=−Fpωp−aiip2sinpπ2−ωτ1+biip2cospπ2−ωτ1,i=1,…,N,where δij is Kronecker’s notation.(25)δij=0,ifi≠j,1,ifi=j.Based on equations (12a)–(12g) and (24), the higher-order approximate analytical result of the strongly nonlinear system with parametric excitation and multifrequency excitation under fractional-order delayed feedback could be obtained. When time delay τ1 is 0, the explicit expression forms of equations (12a)–(12g) and (24) are the same as the explicit forms of the fractional-order derivative in Reference [28], so the correctness of the general formulas in this paper is proved indirectly. By providing the initial value of A0, ΔA can be obtained through iteration. Then the initial value A0 of the next iteration would be get, and so on. Repeat iteration until the accuracy of ΔA meets the requirement.


## 3. Comparison between IHBM and Numerical Integration

In order to analyze the general formulas more correctly, the numerical integration method (NIM) is introduced to simulate equation (1). The numerical calculation form of fractional-order derivative [1, 2] is(26)Dpxtl≈h−p∑j=0lCjpxtl−j,where tl=lh, h is the sample step of calculation, and Cjp is the binomial coefficient which satisfies the following iterative relationship:(27)C0p=1,Cjp=1−1+pjCj−1p.Fractional-order delayed feedback involves time delayτ1 in equation (1). Let τ1=i×h, i is the natural number and then the fractional-order delayed derivative can be written as(28)Dpxt−τ1=Dpxt−ih.According to equations (26)–(28), the numerical iterative algorithm of equation (1) can be expressed as(29)Z1tl=Z2tl−1h−∑j=1lCj1Z1tl−j,Z2tl=F0−∑k=1LFkcoskωtl−1−KpZ3tl−1−i−μZ2tl−1−1+2εcosωtl−1Z1tl−1h−∑j=1lCj1Z2tl−j,Z3tl=Z2tl−1h1−p−∑j=1lCj1−pZ3tl−j,where Z1=xt, Z2=x˙t, and Z3=Dpxt.SelectL = 3, and the system external excitation parameters are as follows: F0=1, F1=0.2, F2=0.05, and F3=0.02. In the following analysis, the five-order approximately analytical solution is considered, that is to say, the value of N in those general formulas of equations (12a)–(12g) and (24) is selected as 5. Here, we can choose a larger value of N to get a higher precision solution. The solution by NIM is shown by the blue circles based on equation (29). The comparisons between the IHBM results in this paper and the NIM results of equation (1) are shown in Figures 1 and 2. It could be found that the amplitude-frequency curves not only with small parameters but also with large parameters are almost coincided by two methods. Therefore, it could be concluded that the precision of the solution by IHBM in this paper is very high, and it is suitable for both the weakly and the strongly nonlinear fractional-order system. From the observation of the amplitude-frequency curves in Figures 1 and 2, it could be found that there exits two main peaks of the amplitude-frequency curve. The amplitude-frequency curve near ω≈1 is the primary resonance response caused by the forcing excitation, and the amplitude-frequency curve near ω≈0.5 is the main parametric resonance response caused by the parametric excitation.Figure 1
The amplitude-frequency curves with small parameters. (μ = 0.1, ε = 0.05, Fp = 0.1, p = 0.5, and τ1 = 0.1).Figure 2
The amplitude-frequency curves with larger parameters. (μ = 0.3, ε = 0.2, Fp = 0.2, p = 0.5, and τ1 = 2).According to the general formulas of IHBM in this paper, the corresponding phase diagrams ofω=1, ω=1/2, and ω=1/3 are shown in Figures 3(a)–3(c), respectively, by the red solid line. From the observation of Figure 3, whether in the primary resonance region (ω≈ω0), the main parametric resonance region (ω≈ω0/2), or the superharmonic resonance region (ω≈ω0/3), it could be found that the solution of equation (1) obtained by IHBM is in good agreement with the solution by NIM. Therefore, for a complex parametric system under multifrequency excitation, the dynamical characteristics of the system in different resonance regions can be simultaneously reflected based on IHBM proposed in this paper.Figure 3
(a) Phase portraitω=1 (μ = 0.1, ε = 0.05, Fp = 0.1, p = 0.5, and τ1 = 0.1). (b) Phase portrait ω  = 1/2 (μ = 0.1, ε = 0.05, Fp = 0.1, p = 0.5, and τ1 = 0.1). (c) Phase portrait ω  = 1/3 (μ = 0.1, ε = 0.05, Fp = 0.1, p = 0.5, and τ1 = 0.1).
(a)(b)(c)


## 4. The Effects of the Feedback Control Parameters

### 4.1. The Effect of the Fractional Order

We have the fractional order from 0 to 2 due to the causality and physical realizability [35], and the amplitude-frequency curves according to the general formulas are obtained, as are shown in Figure 4, where μ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, Kp=0.2, and τ1=0. The amplitude-frequency curve near ω≈1 is the primary resonance response, where the circles represent the maximum amplitude in the primary resonance region with different fractional order, respectively. The amplitude-frequency curve near ω≈0.5 is the main parametric resonance response, where the asterisks represent the maximum amplitude in the main parametric resonance region. From the observation of Figure 4, it could be found that the maximum amplitudes show a trend of first decrease and then increase as the increasing of the fractional order (p∈02), the largest maximum amplitude appears at p=2, and the smallest maximum amplitude appears at p=1. Simultaneously, one could find that whether in the primary resonance region or in the main parametric resonance region, the increasing of the fractional order p would result into the leftwards moving of the maximum amplitude position. This is because the equivalent linear stiffness from fractional-order delayed feedback will become smaller as the fractional order increase, and the equivalent linear damping from fractional-order delayed feedback will become larger at first and then smaller in this procedure. It can be obtained that the vibration behavior of the system will be restrained with the increasing of the fractional order (p∈01). When p=1, the fractional-order delayed feedback is equivalent to the velocity feedback, which brings the maximum damping effect, and at this time the maximum amplitude of the system is minimal.Figure 4
The amplitude-frequency curves.


### 4.2. The Effect of Time Delay

With the change of time delay (τ1∈ [0 4]), the amplitude-frequency curves are shown in Figure 5, where μ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, Kp=0.2, and p = 0.5. It could be found from Figure 5 that maximum amplitudes firstly increase and then decrease with the increase of time delay in the primary resonance region. But maximum amplitudes constantly increase with the increase of time delay in the main parametric resonance region in this procedure. Whether in the primary resonance region or in the main parametric resonance region, the increase of time delay would firstly result into the rightwards moving and then leftwards moving of the maximum amplitude position. The dynamical phenomena shows that the change of time delay will change the stiffness effect and the damping effect from fractional-order delayed feedback.Figure 5
The amplitude-frequency curves.When time delay is increased from 0 to 10, the change of the amplitude-frequency curve with the increase of time delay is shown in Figure6. From Figure 6, with the increase of time delay, maximum amplitudes increase first and then decrease, then increase and decrease, and so on in the resonance region. Those dynamical phenomena is cyclical with the increase of time delay, which could also be obtained from the analysis of the general formulas of equation (24). Through the analysis of equation (24), we could find that time delay affects the dynamic characteristics of the system in the form of trigonometric function. Therefore, it can be inferred that the dynamical phenomena is also cyclical with the increase of time delay in the parametric resonance region. Because of the difference of the resonance frequency between the primary resonance and the parametric resonance, the period of cycle of the primary resonance motion response is faster than the parametric resonance motion response. Vibration amplitudes in different resonance regions can be reduced by choosing appropriate time delay. Those results will help us to design the parameters of fractional-order delayed feedback in the Mathieu equation with multifrequency excitation.Figure 6
The amplitude-frequency-delay graph.From the analysis of equation (24), it could be found that the fractional order affects the dynamical characteristics of the system in the form of trigonometric function too. A reasonable choice of the fractional order can offset the complex dynamical problem caused by time delay. The equation of the fractional order and time delay relation is(30)Δp=2ωτ0π.Selectμ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, and Kp=0.2. When p = 0.2 and τ0=0, the amplitude-frequency curve is shown in Figure 7 by the red solid line. From the comparison of the two curves, it could be found that maximum amplitudes become larger both in the primary resonance region and in the parametric resonance region when time delay is 0.5. Substituting τ0=0.5 into equation (30), we can get Δp≈0.3392. Let the fractional order p = 0.2 + 0.3392 = 0.5392, and the amplitude-frequency curve is shown in Figure 7 by circles. It could be found that those circles almost coincided with the red solid line. Therefore, it can be inferred that the appropriate fractional order can offset the effect of time delay. When time delay is disadvantageous, the complex dynamical problem of the system caused by time delay can be eliminated through selecting different viscoelastic components. This will provide a new design idea for the control strategy of the system.Figure 7
The amplitude-frequency curves.


### 4.3. The Effect of the Feedback Gain Coefficient

When the feedback gain coefficientFp is changed, the amplitude-frequency curves are obtained as are shown in Figure 8, where μ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, p=0.5, and τ1=0.2. From the observation of Figure 8, it could be found that maximum amplitudes decrease both in the primary resonance region and in the parametric resonance region with the increase of the feedback gain coefficient Fp. The amplitude-frequency curve near ω≈0.3 represents the superharmonic resonance response, which will disappear with the larger of the feedback gain coefficient in this region. It can be concluded that the feedback gain coefficient can suppress the maximum vibration amplitude of the system. The resonance frequency in each resonance region is moved to the right with the increase of the feedback gain coefficient. It is shown that the increase of the feedback gain coefficient can increase both the equivalent linear damping and the equivalent linear stiffness of the system.Figure 8
The amplitude-frequency curves.


## 4.1. The Effect of the Fractional Order

We have the fractional order from 0 to 2 due to the causality and physical realizability [35], and the amplitude-frequency curves according to the general formulas are obtained, as are shown in Figure 4, where μ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, Kp=0.2, and τ1=0. The amplitude-frequency curve near ω≈1 is the primary resonance response, where the circles represent the maximum amplitude in the primary resonance region with different fractional order, respectively. The amplitude-frequency curve near ω≈0.5 is the main parametric resonance response, where the asterisks represent the maximum amplitude in the main parametric resonance region. From the observation of Figure 4, it could be found that the maximum amplitudes show a trend of first decrease and then increase as the increasing of the fractional order (p∈02), the largest maximum amplitude appears at p=2, and the smallest maximum amplitude appears at p=1. Simultaneously, one could find that whether in the primary resonance region or in the main parametric resonance region, the increasing of the fractional order p would result into the leftwards moving of the maximum amplitude position. This is because the equivalent linear stiffness from fractional-order delayed feedback will become smaller as the fractional order increase, and the equivalent linear damping from fractional-order delayed feedback will become larger at first and then smaller in this procedure. It can be obtained that the vibration behavior of the system will be restrained with the increasing of the fractional order (p∈01). When p=1, the fractional-order delayed feedback is equivalent to the velocity feedback, which brings the maximum damping effect, and at this time the maximum amplitude of the system is minimal.Figure 4
The amplitude-frequency curves.


## 4.2. The Effect of Time Delay

With the change of time delay (τ1∈ [0 4]), the amplitude-frequency curves are shown in Figure 5, where μ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, Kp=0.2, and p = 0.5. It could be found from Figure 5 that maximum amplitudes firstly increase and then decrease with the increase of time delay in the primary resonance region. But maximum amplitudes constantly increase with the increase of time delay in the main parametric resonance region in this procedure. Whether in the primary resonance region or in the main parametric resonance region, the increase of time delay would firstly result into the rightwards moving and then leftwards moving of the maximum amplitude position. The dynamical phenomena shows that the change of time delay will change the stiffness effect and the damping effect from fractional-order delayed feedback.Figure 5
The amplitude-frequency curves.When time delay is increased from 0 to 10, the change of the amplitude-frequency curve with the increase of time delay is shown in Figure6. From Figure 6, with the increase of time delay, maximum amplitudes increase first and then decrease, then increase and decrease, and so on in the resonance region. Those dynamical phenomena is cyclical with the increase of time delay, which could also be obtained from the analysis of the general formulas of equation (24). Through the analysis of equation (24), we could find that time delay affects the dynamic characteristics of the system in the form of trigonometric function. Therefore, it can be inferred that the dynamical phenomena is also cyclical with the increase of time delay in the parametric resonance region. Because of the difference of the resonance frequency between the primary resonance and the parametric resonance, the period of cycle of the primary resonance motion response is faster than the parametric resonance motion response. Vibration amplitudes in different resonance regions can be reduced by choosing appropriate time delay. Those results will help us to design the parameters of fractional-order delayed feedback in the Mathieu equation with multifrequency excitation.Figure 6
The amplitude-frequency-delay graph.From the analysis of equation (24), it could be found that the fractional order affects the dynamical characteristics of the system in the form of trigonometric function too. A reasonable choice of the fractional order can offset the complex dynamical problem caused by time delay. The equation of the fractional order and time delay relation is(30)Δp=2ωτ0π.Selectμ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, and Kp=0.2. When p = 0.2 and τ0=0, the amplitude-frequency curve is shown in Figure 7 by the red solid line. From the comparison of the two curves, it could be found that maximum amplitudes become larger both in the primary resonance region and in the parametric resonance region when time delay is 0.5. Substituting τ0=0.5 into equation (30), we can get Δp≈0.3392. Let the fractional order p = 0.2 + 0.3392 = 0.5392, and the amplitude-frequency curve is shown in Figure 7 by circles. It could be found that those circles almost coincided with the red solid line. Therefore, it can be inferred that the appropriate fractional order can offset the effect of time delay. When time delay is disadvantageous, the complex dynamical problem of the system caused by time delay can be eliminated through selecting different viscoelastic components. This will provide a new design idea for the control strategy of the system.Figure 7
The amplitude-frequency curves.


## 4.3. The Effect of the Feedback Gain Coefficient

When the feedback gain coefficientFp is changed, the amplitude-frequency curves are obtained as are shown in Figure 8, where μ=0.1, ε=0.2, F0=1, F1=0.2, F2=0.05, F3=0.02, p=0.5, and τ1=0.2. From the observation of Figure 8, it could be found that maximum amplitudes decrease both in the primary resonance region and in the parametric resonance region with the increase of the feedback gain coefficient Fp. The amplitude-frequency curve near ω≈0.3 represents the superharmonic resonance response, which will disappear with the larger of the feedback gain coefficient in this region. It can be concluded that the feedback gain coefficient can suppress the maximum vibration amplitude of the system. The resonance frequency in each resonance region is moved to the right with the increase of the feedback gain coefficient. It is shown that the increase of the feedback gain coefficient can increase both the equivalent linear damping and the equivalent linear stiffness of the system.Figure 8
The amplitude-frequency curves.


## 5. Conclusion

In this paper, the complex dynamical characteristics of the Mathieu equation with multifrequency excitation under fractional-order delayed feedback control are investigated by IHBM. At first, the general formulas of the approximate analytical solution for this system is obtained. The more accurate solution could be obtained through using the obtained general formulas. Then, the numerical result of this system is also studied. The comparison between the result by NIM and the approximate analytical solution by IHBM is given, and the correctness precision of the general forms of the approximate analytical solution is verified. At last, the effects of fractional-order delayed feedback are investigated. Both the fractional order and time delay of fractional-order delayed feedback will affect the dynamical response of the system in the form of trigonometric function. The appropriate fractional order can offset the disadvantageous dynamical problem caused by time delay. The feedback gain coefficient can suppress the maximum vibration amplitude of the system, and so on. Those results could help us to design the control parameters of fractional-order delayed feedback in the Mathieu equation with multifrequency excitation. It can also provide beneficial reference for other fractional-order control systems, even for the strongly nonlinear systems.



---

*Source: 7515080-2022-07-19.xml*