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

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 rst-order approximate periodic solution for the Mathieu equation are derived. Caputo’s denition is adopted to process the fractional-order delayed feedback term.e general formulas of this system are suitable for not only the weakly but also the strongly nonlinear fractional-order system.rough 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. e 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 eects of the time delay, feedback coecient, 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 inuence caused by time delay in fractional-order delayed feedback control. ose new system characteristics will provide theoretical guidance to the design and the control of this kind system.


Introduction
e de nition of fractional calculus has existed for more than 300 years. Fractional calculus has many advantages. For example, it could better re ect 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 di erent research elds are attracted by fractional calculus [1][2][3][4][5].
At present, fractional calculus has played a very important role in di erent elds, such as electrochemistry, signal processing, bioengineering, mechanics, and control. e 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][7][8][9][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. en, the stabilities of the various fractional-order systems are studied, and an e ective 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 di erential equations based on the Simpson method and put forward some numerical calculation methods of fractional-order partial di erential 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. ose 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][23][24][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][27][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. at means that IHBM can be applied into those different nonlinear systems and also is widely used in the engineering field. e 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. e 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. e dynamical phenomena of the system under the combined action of parametric excitation and multifrequency excitation are more complicated [32][33][34]. Accordingly, the general formulas of the Mathieu equation with multifrequency excitation under fractionalorder delayed feedback control based on IHBM are obtained in the following. e dynamical analysis of the Mathieu equation with multifrequency excitation under fractionalorder 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 Section 2. 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. en, 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.

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 where x(t) is the system displacement, μ is the damping coefficient, 2ε cos ωt is the periodic time-varying stiffing coefficient, and F 0 and L k�1 F k cos kωt are the constant excitation and multifrequency excitation, respectively. U(t) is fractional-order delayed feedback control. Here, When p � 0, fractional-order delayed feedback becomes displacement delayed feedback. When p � 1, fractional-order delayed feedback becomes velocity delayed feedback.
(4b) (3), using Taylor series expression, and omitting the higher-order terms of the small increment Δx, one can obtain the differential equation as Defining the following vectors X � [1, cos τ, cos 2 τ, . . . cos nτ, sin τ, one can get 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 � ∞. e 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: Rearrange equation (7), and it yields en, obtaining the 2N + 1 linearized equations about ΔA 0 where and expanding to matrix forms Shock and Vibration 3 one could obtain the explicit expression of M 1 and R 1 as follows: where where δ ij is the Kronecker symbol, and the value of θ is 2π. en, the integration process of fractional-order delayed feedback is presented. M p 2 and R p 2 are the Jacobi matrix and the corrective vector, respectively, and the forms are as follows: According to equations (2) and (15a), it can be written as Letting s � τ − u and du � −ds, and substituting into equation (16), one can get

Shock and Vibration
From equation (17), it could be found the expression for equation (17) has two integrals. Defining the first integral as A 1 and the second integral as A 2 , and integrating it by parts, respectively, one have sin jτ τ p dτ.

(18b)
To simplify the process, we define the first part of equation (18a) as A 11 , the second part of equation (18a) as A 12 , the first part of equation (18b) as A 21 , and the second part of equation (18b) as A 22 . Here, two basic formulas [17] are introduced According to equation (19b), we can get

Shock and Vibration
According to equation (19a), we can get Combining equations (17) Substitute the system original parameters, and equation (22) becomes Similarly, the explicit expressions of the other matrices in equation (14) where δ ij is Kronecker's notation.
Based on equations (12a)-(12g) and (24), the higherorder 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 A 0 , ΔA can be obtained through iteration. en the initial value A 0 of the next iteration would be get, and so on. Repeat iteration until the accuracy of ΔA meets the requirement.

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). e numerical calculation form of fractional-order derivative [1,2] is where t l � lh, h is the sample step of calculation, and C p j is the binomial coefficient which satisfies the following iterative relationship: 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 According to equations (26)- (28), the numerical iterative algorithm of equation (1) can be expressed as where . Select L � 3, and the system external excitation parameters are as follows: F 0 � 1, F 1 � 0.2, F 2 � 0.05, and F 3 � 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. e solution by NIM is shown by the blue circles based on equation (29). e 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. erefore, 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 amplitudefrequency curves in Figures 1 and 2, it could be found that there exits two main peaks of the amplitude-frequency curve. e 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.

Shock and Vibration
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. erefore, for a complex parametric system under multifrequency excitation, the dynamical characteristics of the system in different resonance regions can be simultaneously re ected based on IHBM proposed in this paper.

e E ect 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, F 0 1, F 1 0.2, F 2 0.05, F 3 0.02, K p 0.2, and τ 1 0. e amplitude-frequency curve near ω ≈ 1 is the primary resonance response, where the circles represent the maximum amplitude in the primary resonance region with di erent fractional order, respectively. e 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 rst decrease and then increase as the increasing of the fractional order (p ∈ 0 2 ), the largest maximum amplitude appears at p 2, and the smallest maximum amplitude appears at p 1. Simultaneously, one could nd 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. is is because the equivalent linear sti ness 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 rst 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 ∈ 0 1 ). When p 1, the fractional-order delayed feedback is equivalent to the velocity feedback, which brings the maximum damping e ect, and at this time the maximum amplitude of the system is minimal.

e E ect 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, F 0 1, F 1 0.2, F 2 0.05, F 3 0.02, K p 0.2, and p 0.5. It could be found from Figure 5 that maximum amplitudes rstly 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 rstly result into the rightwards moving and then leftwards moving of the maximum amplitude position. e dynamical phenomena shows that the change of time delay will change the sti ness e ect and the damping e ect from fractional-order delayed feedback.
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 Figure 6. From Figure 6, with the increase of time delay, maximum amplitudes increase rst and then decrease, then increase and decrease, and so on in the resonance region.
ose 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). rough the analysis of equation (24), we could nd that time delay a ects the dynamic characteristics of the system in the form of trigonometric function. erefore, 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 di erence 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 di erent resonance regions can be reduced by choosing appropriate time delay. ose results will help us to design the parameters of fractional-order delayed feedback in the Mathieu equation with multifrequency excitation.
From the analysis of equation (24), it could be found that the fractional order a ects the dynamical characteristics of the system in the form of trigonometric function too. A reasonable choice of the fractional order can o set the complex dynamical problem caused by time delay. e equation of the fractional order and time delay relation is Select μ 0.1, ε 0.2, F 0 1, F 1 0.2, F 2 0.05, F 3 0.02, and K p 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. erefore, it can be inferred that the appropriate fractional order can o set the e ect of time delay. When time delay is disadvantageous, the complex dynamical problem of the system caused by time delay can be eliminated through selecting di erent viscoelastic components.
is will provide a new design idea for the control strategy of the system.

e E ect of the Feedback Gain Coe cient.
When the feedback gain coe cient F p is changed, the amplitudefrequency curves are obtained as are shown in Figure 8, where μ 0.1, ε 0.2, F 0 1, F 1 0.2, F 2 0.05, F 3 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 coe cient F p . e amplitude-frequency curve near ω ≈ 0.3 represents the superharmonic resonance response, which will disappear with the larger of the feedback gain coe cient in this region. It can be concluded that the feedback gain coe cient can suppress the maximum vibration amplitude of the system. e resonance frequency in each resonance region is moved to the right with the increase of the feedback gain coe cient. It is shown that the increase of the feedback gain coe cient can increase both the equivalent linear damping and the equivalent linear sti ness of the system.

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. e more accurate solution could be obtained through using the obtained general formulas. en, the numerical result of this system is also studied. e 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. e appropriate fractional order can offset the disadvantageous dynamical problem caused by time delay. e feedback gain coefficient can suppress the maximum vibration amplitude of the system, and so on.
ose 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.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest.