Bifurcation Control of an Electrostatically-Actuated MEMS Actuator with Time-Delay Feedback

The parametric excitation system consisting of a flexible beam and shuttle mass widely exists in microelectromechanical systems (MEMS), which can exhibit rich nonlinear dynamic behaviors. This article aims to theoretically investigate the nonlinear jumping phenomena and bifurcation conditions of a class of electrostatically-driven MEMS actuators with a time-delay feedback controller. Considering the comb structure consisting of a flexible beam and shuttle mass, the partial differential governing equation is obtained with both the linear and cubic nonlinear parametric excitation. Then, the method of multiple scales is introduced to obtain a slow flow that is analyzed for stability and bifurcation. Results show that time-delay feedback can improve resonance frequency and stability of the system. What is more, through a detailed mathematical analysis, the discriminant of Hopf bifurcation is theoretically derived, and appropriate time-delay feedback force can make the branch from the Hopf bifurcation point stable under any driving voltage value. Meanwhile, through global bifurcation analysis and saddle node bifurcation analysis, theoretical expressions about the system parameter space and maximum amplitude of monostable vibration are deduced. It is found that the disappearance of the global bifurcation point means the emergence of monostable vibration. Finally, detailed numerical results confirm the analytical prediction.


Introduction
Microelectromechanical systems (MEMS) have been widely applied in gyroscopes [1,2], filter [3][4][5] and so on. Parametric resonance in MEMS was first proposed for amplification of harmonically-excited oscillators [6], and since then parametric excitation has been investigated for increasing sensitivity in scanning probe microscopy [7], mass sensing [8], and tuning [9,10]. With the existence of structure nonlinearity and nonlinear electrostatic forces, they can exhibit rich static and dynamic behaviors [11], such as nonlinear jump phenomena [12] and chaos [13]. Rhoads et al. [14] studied a single degree of freedom parametric excitation equation with both the linear and cubic terms and provided a complete description of the dynamic response. Welte et al. [15] studied parametric excitation in a two degree of freedom MEMS. Design parameters were included in the model by lumping them into non-dimensional parameters, thereby allowing for an easier understanding of their effects and the interaction between the mechanical and electrical forces. A variety of nonlinear dynamic behaviors exist in the nonlinear parametric excitation system. However, only a small part of the dynamic behaviors are desired, which requires us to control the bifurcation behavior of the system. To control the dynamic behavior of the system, researchers proposed various parameter optimization methods; for instance, choosing the correct geometry and appropriate voltage. However, it is not enough to improve the performance of the system. To enhance the actuator performance, tip tracking control [16], time-delayed feedback control [17], and pole placement control [18] were introduced. Here, we mainly care about the performance of a time delay feedback controller. Generally, the delayed signal can be displacement, velocity [19], and acceleration [20]. With a properly designed time delay, a delayed feedback controller has been proven to stabilize systems, including atomic force microscopes (AFM) [21] and magneto-elastic beam systems [22]. Shehrin et al. [23] proposed the active control of effective stiffness, damping, and mass of MEMS by applying feedback forces that are proportional to displacement, velocity, and acceleration of its proof mass. Mehta et al. [24] used position feedback to control the effective stiffness of a micro-cantilever to improve the quality factor for biological sensing applications. Morrison et al. [25] investigated the dynamic behaviors of a delayed nonlinear Mathieu equation, and the method of averaging (valid for small ε) was used to obtain a slow flow that was analyzed for stability and bifurcation. Alsaleem et al. [26] presented a study for the stabilization of a MEMS resonator by using a delayed feedback controller. The controller showed a good performance in rejecting disturbances. Warminski [27] analyzed vibrations of a parametrically-excited MEMS device driven by external excitation and time delay inputs. Alsaleem et al. [28] investigated the stability and integrity of parallel-plate MEMS resonators by using a delayed feedback controller. The perturbation method plays an important role in the nonlinear dynamic analysis of MEMS resonators [14]. Kaminski et al. [29,30] studied stochastic nonlinear dynamic behaviors of a MEMS device using the generalized stochastic perturbation technique.
It can be concluded from the above analysis that nonlinear dynamic behaviors and parameter optimization are both important in the design of parametric excitation MEMS and should be taken into account [31][32][33]. Meanwhile, a time delay feedback controller is gradually applied to the design of the MEMS. However, to the best of our knowledge, there are fewer quantitative studies about a general analysis of parametric excitation comb systems with time delay feedback controllers. Additionally, the parametric excitation system consisting of flexible beam and shuttle mass widely exists in MEMS. Early studies mainly focus on a single degree of freedom, which cannot accurately describe the nonlinear dynamic behaviors. In this paper, a new method is used to solve partial differential equations and the detailed mathematical derivation is proposed to quantitatively make a complete description of the transition mechanism of nonlinear jumping phenomena. It is noteworthy that we are mainly concerned with the nonlinear dynamic behavior of MEMS actuator. Here, the influence of parasitic capacitance is not considered.
The structure of this paper is as follows: in Section 2, the partial differential governing equation with parametric excitation and time-delay feedback is obtained; in Section 3, we apply the method of multiple scales directly to the partial differential equation to produce an approximate solution; in Section 4, we analyze the stability and bifurcation near the origin and the discriminant of Hopf bifurcation is theoretically derived; in Section 5, the global bifurcation and saddle note bifurcation are studied. With a time delay feedback controller, monostable vibration is realized. In Section 6, the numerical simulation is given; and Concluding remarks are given in Section 7.

Problem Formulation
Parametrically-excited MEMS was proposed for use in a number of sensing, actuating, resonator, and filtering applications. A conceptual model of an electrostatic comb-finger actuator is shown in Figure 1a, which has been investigated in many studies [5,14]. Such actuators consist of a shuttle mass, namely the actuator's backbone, connected to anchors via beam springs, and excited by a pair of non-interdigitated electrostatic comb drives, which are powered by a voltage source. The actuator's motion is assumed to be described by the movement of the shuttle mass in one direction in the plane. In this paper, multiple electrodes are introduced to optimize the dynamic behavior [34], as shown in Figure 1b. Here, the electrostatic comb-finger actuator is driven by a sensing electrode, a driving electrode, and a controlling electrode. Morrison et al. [25] proved that time delay feedback force can control stability and bifurcation in parametrically-excited nonlinear differential equations. In order to improve the global dynamic behavior of the system, the time delay feedback is introduced as shown in Figure 2, where τ is the time delay and G is the amplitude of the displacement feedback controller with unit V/m. shown in Figure 1b. Here, the electrostatic comb-finger actuator is driven by a sensing electrode, a driving electrode, and a controlling electrode. Morrison et al. [25] proved that time delay feedback force can control stability and bifurcation in parametrically-excited nonlinear differential equations. In order to improve the global dynamic behavior of the system, the time delay feedback is introduced as shown in Figure 2, where τ is the time delay and G is the amplitude of the displacement feedback controller with unit V/m.  The electrostatic attractive force is derived from the very well-known equation: where E is the electric energy stored in the capacitor, i.e., E = CV 2 /2. Here, C represents capacitance. The electrostatic force arising from the driving electrode can be modeled for small displacements [14]: where r1 and r3 are electrostatic coefficients, which can provide harmonic excitation to the device. Here, w is displacement of shuttle mass and V is driving voltage. Then, the electrostatic forces arising from the sensing electrode and controlling electrode can be derived using Equation (1) as [12]: shown in Figure 1b. Here, the electrostatic comb-finger actuator is driven by a sensing electrode, a driving electrode, and a controlling electrode. Morrison et al. [25] proved that time delay feedback force can control stability and bifurcation in parametrically-excited nonlinear differential equations. In order to improve the global dynamic behavior of the system, the time delay feedback is introduced as shown in Figure 2, where τ is the time delay and G is the amplitude of the displacement feedback controller with unit V/m.  The electrostatic attractive force is derived from the very well-known equation: where E is the electric energy stored in the capacitor, i.e., E = CV 2 /2. Here, C represents capacitance. The electrostatic force arising from the driving electrode can be modeled for small displacements [14]: where r1 and r3 are electrostatic coefficients, which can provide harmonic excitation to the device.
Here, w is displacement of shuttle mass and V is driving voltage. Then, the electrostatic forces arising from the sensing electrode and controlling electrode can be derived using Equation (1) as [12]: The electrostatic attractive force is derived from the very well-known equation: where E is the electric energy stored in the capacitor, i.e., E = CV 2 /2. Here, C represents capacitance. The electrostatic force arising from the driving electrode can be modeled for small displacements [14]: where r 1 and r 3 are electrostatic coefficients, which can provide harmonic excitation to the device.
Here, w is displacement of shuttle mass and V is driving voltage. Then, the electrostatic forces arising from the sensing electrode and controlling electrode can be derived using Equation (1) as [12]: where d is the gap width, h is thickness of the microbeam and ε 0 is the dielectric constant of the gap medium. Here, we consider V dc >> Gw(t − τ). Then F control = − ε 0 h d (V 2 dc + 2V dc Gw(t − τ)) is obtained. The electrostatic comb-finger actuator can be described by two clamped-clamped beams and a shuttle mass. An equivalent mechanical model is shown in Figure 3. In this paper, we consider a clamped-clamped Euler microbeam actuated by a concentrated load and subject to a viscous damping c per unit length. By using Hamilton's principle, the equation of motion that governs the transverse deflection w(x,t) is written as [35]: with the boundary conditions: where . w = ∂w/∂t and w = ∂w/∂x, x is the position along the plate length, A and I are the area and moment of inertia of the microbeam, t is time, E is Young's modulus, is the material density, L is the microbeam length, b is the microbeam width, M is mass of shuttle, and δ represents impulse function. where d is the gap width, h is thickness of the microbeam and ε0 is the dielectric constant of the gap medium. Here, we consider Vdc ＞＞ Gw(t − τ).
The electrostatic comb-finger actuator can be described by two clamped-clamped beams and a shuttle mass. An equivalent mechanical model is shown in Figure 3. In this paper, we consider a clamped-clamped Euler microbeam actuated by a concentrated load and subject to a viscous damping c per unit length. By using Hamilton's principle, the equation of motion that governs the transverse deflection w(x,t) is written as [35]: with the boundary conditions: where w w t = ∂ ∂  and w w x ′ = ∂ ∂ , x is the position along the plate length, A and I are the area and moment of inertia of the microbeam, t is time, E is Young's modulus, ρ is the material density, L is the microbeam length, b is the microbeam width, M is mass of shuttle, and δ represents impulse function. For convenience, we introduce the non-dimensional variables: where x0 is a characteristic length of the system. x0 = 2 μm is taken.
Substituting the non-dimensional variables into Equations (5) and (6), yields the following non-dimensional equation of motion: n n n n n n n n n n n n n n n n n n n Ax w w c w w dx w I with boundary conditions: Here, β represents the time delay feedback gain of the system.

Perturbation Analysis
In this section, the method of multiple scales [36] are directly used to investigate the response of the MEMS actuator with small vibration amplitude around an equilibrium position. To indicate the significance of each term in the equation of motion, ε is introduced as a small non-dimensional bookkeeping parameter. Considering that the inertia force of the flexible beam is much smaller than that of the shuttle mass, is given. Then, scaling the damping, nonlinear electrostatic force, periodic excitation, and time delay feedback force, we obtain: For convenience, we introduce the non-dimensional variables: where x 0 is a characteristic length of the system. x 0 = 2 µm is taken. Substituting the non-dimensional variables into Equations (5) and (6), yields the following non-dimensional equation of motion: ..
with boundary conditions: Here, β represents the time delay feedback gain of the system.

Perturbation Analysis
In this section, the method of multiple scales [36] are directly used to investigate the response of the MEMS actuator with small vibration amplitude around an equilibrium position. To indicate the significance of each term in the equation of motion, ε is introduced as a small non-dimensional bookkeeping parameter. Considering that the inertia force of the flexible beam is much smaller than that of the shuttle mass, .. w n = O(ε) is given. Then, scaling the damping, nonlinear electrostatic force, periodic excitation, and time delay feedback force, we obtain: w n + r 1n w n V 2 + εr 1n w n V 2 cosΩ n t n +εr 3n w 3 n V 2 (1 + cosΩ n t n ) − εβ n w n (t n − τ n )] (9) with the boundary conditions: where α = Ax 2 0 /2I. In order to investigate the primary parametric resonance when the driving frequency is close to two times of the natural frequency, a detuning parameter σ is introduced and defined by: where ω is the natural frequency of the system. We seek the approximate solution of Equation (9) in the form: where T 0 = t and T 1 = εt. Substituting Equations (11) and (10) into Equation (9) and equating coefficients of like powers of ε 0 and ε 1 , yield ε 0 : Considering a single mode vibration, the general solution of Equation (12) can be written as: where cc indicates the complex conjugate of the preceding terms. Substituting Equation (16) into Equations (12) and (13), yields: where ω can be defined by Equation (12). We obtain ω = (192 + r 1n V 2 )/M n . The general solution of Equation (14) can be written as: Substituting Equation (18) into Equation (14), multiplying by ϕ 0 , and integrating the outcome from x = 0 to 1, we obtain [36,37]: Substituting Equations (10) and (17) into Equation (19) and eliminating the secular term, we obtain: where the overbar indicates the complex conjugate and |A| represents modulus of A. At this point, it is convenient to express A in the polar form: where a is the amplitude. Substituting Equation (21) into Equation (20), and separating the imaginary and real parts, yields: where ψ = θ − σ/2 is the phase of oscillator's response, and several variables are defined as: The steady-state response can be obtained by imposing the condition a' = ψ = 0. The first-order approximate solution is obtained:

Bifurcation Control near the Origin
For determining the steady states of this system, it turns out to be advantageous to introduce the new unknown variables. We obtain an alternate form of the equation by transforming from polar coordinates a and ψ to rectangular coordinates u and v, where: It is noted that u represents the displacement signal and v represents the velocity signal when t = 0.

Stability of the Origin
The trace and determinant of the Jacobian matrix evaluated at an equilibrium point contain the local stability information. From Equations (22) and (23), 8ξ + 4gsinωτ n represents the effective damping of the system. When it is negative, all of the solutions are not stable. Here, 8ξ + 4gsinωτ n > 0 is considered. For stability of the origin, a critical point generically occurs when Det(J) = 0. From Equation (28), the Jacobian matrix value of the origin is obtained: Through Equation (29), it is found that the stability of the origin varies with g and τ n . Here, according to the bifurcation behaviors of the origin, two different situations are divided. Case one: when (8ξ + 4gsinωτ n ) 2 − 4λ 1 2 ≥ 0, the origin is always stable. Case two: when (8ξ + 4gsinωτ n ) 2 − 4λ 1 2 < 0, there exists resonance in the system. When Det(J) = 0, the bifurcation points can be obtained: The origin is unstable from σ 1 to σ 2 , which means resonance occurs. Then, the size of the unstable region is obtained: From Equations (30)-(32), it is found that a negative feedback gain can increase the resonance frequency and widens the resonance band when ωτ n < π/2. However, the maximum of the resonance band is λ 1 , which is decided by the linear electrostatic force.
In this paper, the model is improved from the traditional one. Part of the system parameters are defined as stated in Table 1 [5]. Figure 4 shows the stability of origin in the case of V = 30 V. It is found that the positive feedback gain can increase the resonance frequency when ωτ n = π and decrease the frequency band when ωτ n = π/2. Through choosing the right time delay, we can control the resonance frequency and frequency band individually. Figure 4c,d study the resonance frequency and frequency band variance with β = 0.225 and β = 0.36. When feedback gain is more than the critical value λ 1 /2 − 2ξ, the saddle-node bifurcation occurs as shown in Figure 4d. In practical engineering, we can choose the right time delay and gain to meet the engineering demand.  Figure 4 shows the stability of origin in the case of V = 30 V. It is found that the positive feedback gain can increase the resonance frequency when ωτn = π and decrease the frequency band when ωτn = π/2. Through choosing the right time delay, we can control the resonance frequency and frequency band individually. Figure 4c,d study the resonance frequency and frequency band variance with β = 0.225 and β = 0.36. When feedback gain is more than the critical value λ1/2 − 2ξ, the saddle-node bifurcation occurs as shown in Figure 4d. In practical engineering, we can choose the right time delay and gain to meet the engineering demand.

Stability of the Untrivial Solution
Stability of the origin was studied in the previous section. Through Equations (26) and (27), it is found that Hopf bifurcation will occur at σ1 and σ2, which can lead to an untrivial solution and change the stability of the origin. Then, stability of the untrivial solution is investigated. The supercritical Hopf bifurcation of σ1 and σ2 can lead to stable branches between two critical points. On the contrary, subcritical Hopf bifurcation of σ1 and σ2 can lead to unstable branches. In this section, we study Hopf bifurcation of critical points to determine the stability of periodic vibration.

Stability of the Untrivial Solution
Stability of the origin was studied in the previous section. Through Equations (26) and (27), it is found that Hopf bifurcation will occur at σ 1 and σ 2 , which can lead to an untrivial solution and change the stability of the origin. Then, stability of the untrivial solution is investigated. The supercritical Hopf bifurcation of σ 1 and σ 2 can lead to stable branches between two critical points. On the contrary, subcritical Hopf bifurcation of σ 1 and σ 2 can lead to unstable branches. In this section, we study Hopf bifurcation of critical points to determine the stability of periodic vibration. A series of calculations is made in Appendix A, and the discriminants of Hopf bifurcation points are obtained: where ∆ 1 stands for the discriminant of σ 1 , and ∆ 2 stands for the discriminant of σ 2 . The case ∆ 1 < 0 and ∆ 2 < 0 results in the supercritical Hopf bifurcation of σ 1 and σ 2 . Likewise, the case ∆ 1 > 0 and ∆ 2 > 0 results in the subcritical Hopf bifurcation of σ 1 and σ 2 . However, due to the nature of this equation of motion, two mixed cases also exist, namely, ∆ 1 > 0 and ∆ 2 < 0, or ∆ 1 < 0 and ∆ 2 > 0, which correspond to the two branches from the critical point bending toward each other. When ∆ 1 = 0 or ∆ 2 = 0, the type of bifurcation cannot be captured without the inclusion of high-order nonlinear terms, which is not considered here. Figure 5 shows different types of bifurcation diagrams. The stability of solution can be obtained by Equation (28), as mentioned in [14]. Meanwhile, the values of the discriminant are obtained: Among them, ∆ 1 < 0, supercritical Hopf bifurcation occurs at σ 1 . Similarly, with ∆ 2 (20) > 0, subcritical Hopf bifurcation occurs at σ 2 . With ∆ 2 (30) < 0 and ∆ 2 (40) < 0, supercritical Hopf bifurcation occurs at σ 2 , which coincides with the numerical results of Equations (22) and (23). From V = 20 V to V = 40 V, the resonance frequency increases. Additionally, feedback force can also change the resonance frequency and transform the types of bifurcation points, as shown in Figure 5d. ∆ 2 (20) < 0 is obtained in the case of β = 0.354 and ωτ n = π/4. Subcritical Hopf bifurcation of the point σ 2 is transformed to supercritical Hopf bifurcation.
where 1 Δ stands for the discriminant of σ1, and 2 Δ stands for the discriminant of σ2. The case 1 Δ < 0 and 2 Δ < 0 results in the supercritical Hopf bifurcation of σ1 and σ2. Likewise, the case 1 Δ > 0 and 2 Δ > 0 results in the subcritical Hopf bifurcation of σ1 and σ2. However, due to the nature of this equation of motion, two mixed cases also exist, namely, 1 Δ > 0 and 2 Δ < 0, or 1 Δ < 0 and 2 Δ > 0, which correspond to the two branches from the critical point bending toward each other. When 1 Δ = 0 or 2 Δ = 0, the type of bifurcation cannot be captured without the inclusion of high-order nonlinear terms, which is not considered here. Figure 5 shows different types of bifurcation diagrams. The stability of solution can be obtained by Equation (28), as mentioned in [14]. Meanwhile, the values of the discriminant are obtained:   Through Equations (33) and (34), we can obtain that supercritical Hopf bifurcation always occurs at σ 1 and σ 2 when 9k 2 3 − 4λ 2 3 < 0. However, when 9k 2 3 − 4λ 2 3 > 0, subcritical Hopf bifurcation may occur at σ 1 in case of k 3 < 0 and occur at σ 2 in the case of k 3 > 0. Meanwhile, when gsinωτ n is close to λ 1 /2 − 2ξ, ∆ 1 and ∆ 2 are negative. Thus, with an appropriate feedback force, supercritical Hopf bifurcation always occurs under any driving voltage. The detailed derivation process is seen in Appendix B. Figure 6 shows time delay feedback control on Hopf bifurcation. As shown in Figure 6a, with a voltage value less than 64.96 V, the branch from σ 1 is always stable. However, the branch becomes unstable with the increase of the voltage. On the contrary, with the voltage value is greater than 29.04 V, the branch from σ 2 is always stable. However, the branch becomes unstable with the decrease of the voltage, as shown in Figure 6b. With the introduction of the time-delay feedback force, the unstable branches become stable. The curves shown in Figure 6 represent the minimum value of gsinωτ n to make the branches stable.
Through Equations (33) and (34), we can obtain that supercritical Hopf bifurcation always occurs at σ1 and σ2 when supercritical Hopf bifurcation always occurs under any driving voltage. The detailed derivation process is seen in Appendix B. Figure 6 shows time delay feedback control on Hopf bifurcation. As shown in Figure 6a, with a voltage value less than 64.96 V, the branch from σ1 is always stable. However, the branch becomes unstable with the increase of the voltage. On the contrary, with the voltage value is greater than 29.04 V, the branch from σ2 is always stable. However, the branch becomes unstable with the decrease of the voltage, as shown in Figure 6b. With the introduction of the time-delay feedback force, the unstable branches become stable. The curves shown in Figure 6 represent the minimum value of sin ωτ n g to make the branches stable.

Global Bifurcation
The monostable vibration, which can eliminate jump phenomena between multiple stable solutions and improve system stability, has been widely applied in engineering, such as resonators. In Section 4, there is one or more periodic solutions under the single frequency. To obtain the parameter space of the monostable vibration, global dynamics analysis is introduced. This is almost unaffected by damping, since it is small, and zero damping is assumed to simplify the analysis.
First, we consider the case of V = 30 V without time-delay feedback force. The bifurcation diagram of the no-damping system is shown in Figure 7. Here, four frequency points 0, 0.085, 1, and 2 are taken to analyze the global dynamics. There exists one periodic solution at 0, 2 and two periodic solutions at 0.085, 1, as shown in Figure 7.
Representative phase planes corresponding to each of the frequency response points are delineated in Figure 8. Figure 8a describes an unstable origin and a center, which corresponds to a in Figure 7. Both Figure 8b and Figure 8c depict unstable origins and two centers. It is worth noting that different initial disturbance can lead to different movement tracks. The system easily arrives at b, e by velocity disturbance. Likewise, the system easily arrives at c, d by displacement disturbance. When σ = 0.2, the origin becomes a center, as shown in Figure 8d. From Figure 8b,c, the global bifurcation occurs at the intersection, as shown in Figure 7. Meanwhile, the phenomena of multistability exist in the system.

Global Bifurcation
The monostable vibration, which can eliminate jump phenomena between multiple stable solutions and improve system stability, has been widely applied in engineering, such as resonators. In Section 4, there is one or more periodic solutions under the single frequency. To obtain the parameter space of the monostable vibration, global dynamics analysis is introduced. This is almost unaffected by damping, since it is small, and zero damping is assumed to simplify the analysis.
First, we consider the case of V = 30 V without time-delay feedback force. The bifurcation diagram of the no-damping system is shown in Figure 7. Here, four frequency points 0, 0.085, 1, and 2 are taken to analyze the global dynamics. There exists one periodic solution at 0, 2 and two periodic solutions at 0.085, 1, as shown in Figure 7.  Representative phase planes corresponding to each of the frequency response points are delineated in Figure 8. Figure 8a describes an unstable origin and a center, which corresponds to a in Figure 7. Both Figures 8b and 8c depict unstable origins and two centers. It is worth noting that different initial disturbance can lead to different movement tracks. The system easily arrives at b, e by velocity disturbance. Likewise, the system easily arrives at c, d by displacement disturbance. When σ = 0.2, the origin becomes a center, as shown in Figure 8d. From Figure 8b,c, the global bifurcation occurs at the intersection, as shown in Figure 7. Meanwhile, the phenomena of multistability exist in the system.  The above analysis results are based on an undamped system, the existence of the damping brings the center into focus and shifts each of their locations at the same time.

Saddle Node Bifurcation
Through Equation (22), when the amplitude is in a certain range, there is no periodic vibration, as shown in Figure 5. The existence of periodic vibration requires: The above analysis results are based on an undamped system, the existence of the damping brings the center into focus and shifts each of their locations at the same time.

Saddle Node Bifurcation
Through Equation (22), when the amplitude is in a certain range, there is no periodic vibration, as shown in Figure 5. The existence of periodic vibration requires: When equal, saddle node bifurcation occurs, as shown in Figure 9. From Equation (23), the perturbation frequency of the node bifurcation point is obtained: where: a 2 = 4gsinωτ n + 8ξ − 2λ 1 λ 3 (37) or: The branch from the bifurcation point σ is unstable, which can be proved by Equations (22) and (23).
perturbation frequency of the node bifurcation point is obtained: 3 cos 4 σ = + − ωτ a k k g (36) where:  (37) or: The branch from the bifurcation point σ is unstable, which can be proved by Equations (22) and (23).

Monostable Parameter Vibration
Now, time-delay feedback force is introduced to eliminate the global bifurcation point and make the system translate to a monostable one from a bistable one. When the saddle node bifurcation point and global bifurcation point coincide, the bistable phenomenon disappears, and there is only one stable periodic solution under the single frequency.
From Equation (23), we obtain: It is found that there is only one amplitude corresponding to one frequency when global bifurcation occurs, which becomes true when 2λ1 + 2a 2 λ3 = 0. In other words, 1 3 −λ λ is the maximum amplitude under monostable parameter vibration, which is only related to the structure parameter. When the saddle node bifurcation point and global bifurcation point coincide, the following expression is obtained with Equations (36) and (39): From Equation (40), it is noted that the maximum amplitude of the system is less than

Monostable Parameter Vibration
Now, time-delay feedback force is introduced to eliminate the global bifurcation point and make the system translate to a monostable one from a bistable one. When the saddle node bifurcation point and global bifurcation point coincide, the bistable phenomenon disappears, and there is only one stable periodic solution under the single frequency.
From Equation (23), we obtain: It is found that there is only one amplitude corresponding to one frequency when global bifurcation occurs, which becomes true when 2λ 1 + 2a 2 λ 3 = 0. In other words, √ −λ 1 /λ 3 is the maximum amplitude under monostable parameter vibration, which is only related to the structure parameter.
When the saddle node bifurcation point and global bifurcation point coincide, the following expression is obtained with Equations (36) and (39): From Equation (40), it is noted that the maximum amplitude of the system is less than √ −λ 1 /λ 3 in the case of λ 1 /2 − 2ξ > gsinωτ n > λ 1 /4 − 2ξ. Thus, there is only one periodic solution under this parameter space.
Reasonable design of time-delay feedback in engineering practice can make the system reach a monostable motion state, as shown in Figure 10. This shows the motion state of the system in the different voltage and gain for ωτ n = π/2. Figure 11 shows the bifurcation diagram of the monostable vibration with V = 40 V and β = 0.55.
From Equation (40), although delay feedback force can eliminate the global bifurcation point and make the system translate to a monostable one from a bistable one, it cannot change the maximum amplitude of vibration, which is only related to the structure parameters of the comb. Additionally, in the Section 4, when the voltage is too large or too small, the branches from σ 1 and σ 2 become unstable. A large enough force feedback is introduced to make them stable and reduce the maximum amplitude at the same time. In the Figure 12, the relationship between the maximum amplitude of the monostable vibration and the voltage is given. It shows that the maximum amplitude is constant when the voltage is greater than 28.96 V. On the contrary, the maximum amplitude decreases with the decrease of voltage when the voltage is less than 28.96 V.
Reasonable design of time-delay feedback in engineering practice can make the system reach a monostable motion state, as shown in Figure 10. This shows the motion state of the system in the different voltage and gain for ωτn = π/2. Figure 11 shows the bifurcation diagram of the monostable vibration with V = 40 V and β = 0.55.
From Equation (40), although delay feedback force can eliminate the global bifurcation point and make the system translate to a monostable one from a bistable one, it cannot change the maximum amplitude of vibration, which is only related to the structure parameters of the comb. Additionally, in the Section 4, when the voltage is too large or too small, the branches from σ1 and σ2 become unstable. A large enough force feedback is introduced to make them stable and reduce the maximum amplitude at the same time. In the Figure 12, the relationship between the maximum amplitude of the monostable vibration and the voltage is given. It shows that the maximum amplitude is constant when the voltage is greater than 28.96 V. On the contrary, the maximum amplitude decreases with the decrease of voltage when the voltage is less than 28.96 V.    monostable motion state, as shown in Figure 10. This shows the motion state of the system in the different voltage and gain for ωτn = π/2. Figure 11 shows the bifurcation diagram of the monostable vibration with V = 40 V and β = 0.55.
From Equation (40), although delay feedback force can eliminate the global bifurcation point and make the system translate to a monostable one from a bistable one, it cannot change the maximum amplitude of vibration, which is only related to the structure parameters of the comb. Additionally, in the Section 4, when the voltage is too large or too small, the branches from σ1 and σ2 become unstable. A large enough force feedback is introduced to make them stable and reduce the maximum amplitude at the same time. In the Figure 12, the relationship between the maximum amplitude of the monostable vibration and the voltage is given. It shows that the maximum amplitude is constant when the voltage is greater than 28.96 V. On the contrary, the maximum amplitude decreases with the decrease of voltage when the voltage is less than 28.96 V.    monostable motion state, as shown in Figure 10. This shows the motion state of the system in the different voltage and gain for ωτn = π/2. Figure 11 shows the bifurcation diagram of the monostable vibration with V = 40 V and β = 0.55.
From Equation (40), although delay feedback force can eliminate the global bifurcation point and make the system translate to a monostable one from a bistable one, it cannot change the maximum amplitude of vibration, which is only related to the structure parameters of the comb. Additionally, in the Section 4, when the voltage is too large or too small, the branches from σ1 and σ2 become unstable. A large enough force feedback is introduced to make them stable and reduce the maximum amplitude at the same time. In the Figure 12, the relationship between the maximum amplitude of the monostable vibration and the voltage is given. It shows that the maximum amplitude is constant when the voltage is greater than 28.96 V. On the contrary, the maximum amplitude decreases with the decrease of voltage when the voltage is less than 28.96 V.

Numerical Simulation
The above analysis is carried out based on the perturbation theory. This section gives the numerical results of the partial differential equation obtained by the finite difference method and long-time integration to verify the validity of the perturbation theory. Here, the foregoing analytical results based on the slow flow Equations (22) and (23) are compared with direct numerical integration. The numerical integration is completed by using fourth-order Runge-Kutta method with fixedstep. Figure 13 shows the frequency response curve. Here, we choose a small displacement disturbance as the initial value. The numerical method has yielded results that are in good agreement with the analytical solution, especially near the origin.

Numerical Simulation
The above analysis is carried out based on the perturbation theory. This section gives the numerical results of the partial differential equation obtained by the finite difference method and long-time integration to verify the validity of the perturbation theory. Here, the foregoing analytical results based on the slow flow Equations (22) and (23) are compared with direct numerical integration. The numerical integration is completed by using fourth-order Runge-Kutta method with fixedstep. Figure 13 shows the frequency response curve. Here, we choose a small displacement disturbance as the initial value. The numerical method has yielded results that are in good agreement with the analytical solution, especially near the origin.

Conclusions
This article theoretically investigates the nonlinear jumping phenomena and bifurcation conditions of a class of electrostatically-driven MEMS actuators with a time-delay feedback controller. A detailed analysis method is proposed to deal with the parametric excitation system consisting of a flexible beam and shuttle mass. It is found that our model can improve the resonance frequency and stability of the system. Meanwhile, an appropriate time-delay feedback force can make the branch from the Hopf bifurcation point stable under any driving voltage value.
Additionally, monostable vibration can eliminate dynamic bifurcation and improve system stability, which is desired in many MEMS applications. Theoretical expressions about the system parameter space and maximum amplitude of monostable vibration are deduced. It is found that the disappearance of the global bifurcation point means the emergence of monostable vibration. A method is proposed to translate the bistable state to a monostable state with time-delay feedback.
Finally, detailed numerical results confirm the analytical prediction. The analysis presented here describes a complete picture of this behavior, and provides designers of these devices with useful predictive tools.

Conclusions
This article theoretically investigates the nonlinear jumping phenomena and bifurcation conditions of a class of electrostatically-driven MEMS actuators with a time-delay feedback controller. A detailed analysis method is proposed to deal with the parametric excitation system consisting of a flexible beam and shuttle mass. It is found that our model can improve the resonance frequency and stability of the system. Meanwhile, an appropriate time-delay feedback force can make the branch from the Hopf bifurcation point stable under any driving voltage value.
Additionally, monostable vibration can eliminate dynamic bifurcation and improve system stability, which is desired in many MEMS applications. Theoretical expressions about the system parameter space and maximum amplitude of monostable vibration are deduced. It is found that the disappearance of the global bifurcation point means the emergence of monostable vibration. A method is proposed to translate the bistable state to a monostable state with time-delay feedback.
Finally, detailed numerical results confirm the analytical prediction. The analysis presented here describes a complete picture of this behavior, and provides designers of these devices with useful predictive tools. Author Contributions: L.L. and Q.Z. conceived and designed the model; L.L. and W.W. contributed theoretical analysis; L.L., W.W., and J.H. analyzed the data; L.L. wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.