Parametric Characteristics and Bifurcation Analysis of Multi-Degree-of-Freedom Micro Gyroscope with Drive Stiffness Nonlinearity

The dynamic equations of a four-degree-of-freedom micro gyroscope system were developed considering the nonlinearity of driving stiffness, the primary resonance, and the 1:1 internal resonance. Then, the perturbation analysis was carried out using the method of multiple scales. The influence of stiffness nonlinearity and system parameters on micro-gyro dynamic characteristics, output sensitivity, detection bandwidth, and working stability were discussed based on the analytic and numerical solutions of the dynamic equations. Through the singularity theory, the influence of system parameters on bifurcation behavior was analyzed. The results show that the amplitude jump and multi-stable solutions caused by the nonlinear hardening characteristics of the high robust two-degree-of-freedom drive-mode occur outside the detection bandwidth. In addition, the influence on the bandwidth was weak and the sensitivity of the bandwidth area was slightly reduced. Moreover, saturation existed in the response amplitude of the second drive-mode in spite of the primary resonance being completely tuned or detuned. As a result, although the electrostatic force amplitude was out of the unstable region and even took a larger value, the micro gyroscope obtained a larger stable output. Besides, nonlinearity will lead to energy transfer between various modes of multi-degree-of-freedom micro gyroscopes. That means the response amplitudes could change greatly due to the variation of the external environment even the system is under a constant excitation frequency. Therefore, increasing the stiffness coefficient of the micro beam and the electrostatic force amplitude can maintain the robustness of the system to environmental changes and avoid the occurrence of bifurcation.


Introduction
The micro-electro-mechanical gyroscope is a kind of inertial device used to measure angular velocity and angular displacement, which is developed on the basis of micro-electro-mechanical system (MEMS) processing technology and measurement and control technology. It is widely used in many fields, such as aerospace, weapon guidance, automotive navigation, and robotics, because of its small size, low power consumption, reliability and stability, and batch production [1][2][3][4]. The drive and sense mode of the traditional micro-mechanical gyros mainly utilized the single-degree-of-freedom (1-DOF) structure, which increases the mechanical sensitivity by matching the natural frequencies of the drive and sense mode [5,6]. However, mode matching also leads to a significant reduction in bandwidth and discussed. Wang et al. [28] proposed a method using the stiffness nonlinearity of the 1-DOF micro gyroscope in the driving direction. The inherent hardening characteristics of the stiffness nonlinearity were used to match the resonant frequency to achieve a higher amplitude than the linear design. Using experiments and simulations, Xu et al. [29] verified that the amplitude-frequency curves of the gyroscope have a wide flat region and a higher amplitude than the linear design when the driving direction is nonlinear. A new driving method was proposed to excite the large amplitude nonlinear vibration in the downward scanning characteristic curve. Lajimi et al. [30] presented the nonlinear dynamical features of a gyroscopic system manifesting in a rotation rate sensor. A computational shooting method and Floquet multipliers were used to characterize the response. The study showed that larger bandwidth and higher sensitivity appeared when the system operating in the nonlinear regime. Wen et al. [31] investigated the design principle of the detection bandwidth in micro gyros when nonlinear stiffness existed in driving microbeams and sensing microbeams at the same time. Tatar et al. [20] successfully linearized the electrostatic nonlinearity at the driving comb using a formed comb with a tuned cubic hardening compensation in a triple symmetric silicon-on-insulator (SOI)-MEMS micro gyroscope. Ding et al. [32] proposed an improved digital phase-locked driving method based on the comparison of two methods for controlling stiffness nonlinearity under large amplitude. The stability of the control loop was significantly improved at the expense of the larger driving force. Shang et al. [33] studied a 1-DOF micro-gyroscope model driven by parametric excitation with stiffness nonlinearity using analytical and numerical methods. The influence mechanism of system parameters on amplitude-frequency characteristics and bifurcation behavior of driving and sensing modes was revealed. It was found that the variation of excitation frequency has the potential to cause complex dynamic behaviors when the micro gyro vibration system was under 1:1 internal resonance or large carrier angular velocity, such as multi-stable solution, amplitude jump, and almost periodic response.
Parametric excitation for MEMS gyroscopes can provide resonance in both the drive and the sense modes, even with mismatched natural frequencies. Pakniyat et al. [34] studied the requirements for such a condition by analyzing the effect of each factor on the steady state amplitudes of the two modes. By comparing the gyroscope with matched natural frequencies with a gyroscope with mismatching modes, the result showed that parametric excitation was able to provide high accuracy and robustness for MEMS gyroscopes. Min et al. [35] investigated the partial and full chaotic synchronizations of two nonlinear gyroscope systems with/without noise. The study showed that the simple feedback control can make the noised gyroscope system synchronizing with chaotic behaviors of the expected gyroscope system, and a novel synchronization methodology was presented.
To date, the traditional 2-DOF micro gyroscope, which utilizes the 1-DOF drive-mode and sense-mode, is the mainstream of micro gyroscopes nonlinear dynamics. In fact, nonlinear phenomena are more common in micro gyros with MDOF drive-mode and sense-mode, and the nonlinear dynamic behavior is more complicated. With the gradual popularization of MDOF micro gyro applications [36][37][38][39][40][41][42], facing the possible nonlinear dynamic problems in the MDOF micro gyro system is unavoidable. In the MDOF micro gyroscope system, since the displacement of the sense module is much smaller than that of the drive module, the stiffness nonlinearity resulting from the large deformation of the elastic beam in the drive module could become the main nonlinearity source. Therefore, a new type of micro gyroscope system with 2-DOF drive mode and 2-DOF sense mode was proposed in this study, which takes the stiffness nonlinearity of two drive modes into account simultaneously.
In this paper, the gyroscope's nonlinear dynamic model was developed, which considered the stiffness nonlinearity of two drive modes. The method of multiple scales (MMS) was applied to determine the response of the system under the conditions of primary resonance and 1:1 internal resonance. Then, the Runge-Kutta method was applied to verify the theoretical results. The influence of system parameters on the dynamic characteristics of drive and sense modes were discussed, such as sensitivity, bandwidth, stability, etc. The local bifurcation analysis of the system was carried out.
Through the study of the bifurcation behavior in different parameter spaces, the influence of system parameters on the environmental robustness was discussed, which provides theoretical guidance for the design of the MDOF micro gyroscope structure.

Operational Principle of the Gyroscope
In this paper, a typical 4-DOF micromachinery gyroscope with 2-DOF drive-mode and sense-mode [13] was considered. The schematic diagram of this gyroscope is shown in Figure 1. This gyroscope is mainly composed of driving mass, decoupled frame, proof mass, detection mass, elastic microbeams, and comb electrodes. As shown in Figure 1, the drive direction is along x axis, the sense direction is along y axis, and Ω z is the input angular velocity perpendicular to the x-y plane. The decoupled mass m f and the proof mass m 2 form a two-stage decoupling structure, which can isolate the drive-mode and the sense-mode. The driving mass m 1 vibrates in x direction under the effect of the driving comb electrodes, and the decoupled mass m f starts to vibrate in x direction due to the effect of the microbeam k 2 when the gyroscope works. At the same time, the proof mass m f vibrates along the x direction with the decoupled mass under the effect of the microbeam k 4 . Because of the Coriolis effect, the vibration in x direction causes the resonance in y direction when the system has angular velocity Ω z input in the vertical direction of the x-y plane. Then, the proof mass m 2 and the detection mass m 3 vibrate along the y direction under the constraint of the microbeam k 4 , k 5 , and k 6 . The displacement of detection mass m 3 in y direction is the detection output of the gyroscope, which increases as the angular velocity increases. The detection output amplitude is proportional to the input angular velocity Ω z during structural resonance, so the input angular velocity Ω z of the carrier can be obtained by measuring the output amplitude.

Operational Principle of the Gyroscope
In this paper, a typical 4-DOF micromachinery gyroscope with 2-DOF drive-mode and sensemode [13] was considered. The schematic diagram of this gyroscope is shown in Figure 1. This gyroscope is mainly composed of driving mass, decoupled frame, proof mass, detection mass, elastic microbeams, and comb electrodes. As shown in Figure 1, the drive direction is along x axis, the sense direction is along y axis, and Ωz is the input angular velocity perpendicular to the x-y plane. The decoupled mass mf and the proof mass m2 form a two-stage decoupling structure, which can isolate the drive-mode and the sense-mode. The driving mass m1 vibrates in x direction under the effect of the driving comb electrodes, and the decoupled mass mf starts to vibrate in x direction due to the effect of the microbeam k2 when the gyroscope works. At the same time, the proof mass mf vibrates along the x direction with the decoupled mass under the effect of the microbeam k4. Because of the Coriolis effect, the vibration in x direction causes the resonance in y direction when the system has angular velocity Ωz input in the vertical direction of the x-y plane. Then, the proof mass m2 and the detection mass m3 vibrate along the y direction under the constraint of the microbeam k4, k5, and k6. The displacement of detection mass m3 in y direction is the detection output of the gyroscope, which increases as the angular velocity increases. The detection output amplitude is proportional to the input angular velocity Ωz during structural resonance, so the input angular velocity Ωz of the carrier can be obtained by measuring the output amplitude. Because the mass of the elastic beam is far less than the vibrating mass, it can be neglected when the micro gyroscope rotates at a constant angular velocity in the x-y plane. Therefore, the 4-DOF lumped parameter model can be used to describe the vibration of the micro gyroscope in the x-y plane, which is shown in Figure 2.  Because the mass of the elastic beam is far less than the vibrating mass, it can be neglected when the micro gyroscope rotates at a constant angular velocity in the x-y plane. Therefore, the 4-DOF lumped parameter model can be used to describe the vibration of the micro gyroscope in the x-y plane, which is shown in Figure 2.
In addition, the working environment of such a gyroscope is consistent with [13], the air damping is relatively small, and the nonlinear factor of damping can be ignored. Thus, it can be assumed that the damping of the drive and sense directions in the system are linear damping. The dynamic equations of the drive and sense directions of the gyroscope system are established by Figure 2, respectively. The dynamic equations are written as: Drive direction: .. ..
where m i (i = 1, 2, 3, f ) represents the mass of each mass, c i (i = 1, 2, 3, 4, 5, 6) represents each damping coefficient, k i (i = 1, 2, 3, 4, 5, 6) represents the stiffness coefficient of each elastic microbeam, x i (i = 1, 2) and y i (i = 1, 2) represent the displacement of the i-th degree-of-freedom in x direction and y direction respectively, and F d and F c are electrostatic driving force and Coriolis force, respectively.
x 2 . F is the amplitude of electrostatic driving force, ω 0 is the frequency of electrostatic driving force, and Ω z is input angular velocity of the gyroscope. The values of these physical parameters [13] are shown in Table 1. Because the mass of the elastic beam is far less than the vibrating mass, it can be neglected when the micro gyroscope rotates at a constant angular velocity in the x-y plane. Therefore, the 4-DOF lumped parameter model can be used to describe the vibration of the micro gyroscope in the x-y plane, which is shown in Figure 2.

Perturbation Analysis
In order to consider the case where the model drive beams have nonlinear stiffness, a cubic stiffness term was added to the original drive direction dynamic equations. The nonlinear dynamic equations can be obtained: ..
where K i (i = 1, 2) represents the stiffness nonlinear coefficient of the i-th degree of freedom in the drive direction. Simplify Equation (3), and the following equation can be obtained: .. where Equation (4) is the forced vibration of Duffing system with damping under harmonic excitation. The MMS is used to investigate the approximate periodic response of Equation (4). In order to study the dynamic behavior of the system near the internal resonance, the primary resonance and 1:1 internal resonance of the system were considered. To describe the nearness of the resonance, two detuning parameters σ 1 and σ 2 were introduced and defined by: where ω 1 and ω 2 are the resonance frequencies of the first-order mode and the second-order mode in the drive direction, respectively. ε is introduced as a small nondimensional bookkeeping parameter.
Since the primary resonance of ω 0 close to ω 1 was studied, the electrostatic force term F = O(ε 3 ) and the damping term c i = O(ε 2 ) were considered here for making the effects of the damping term and the forcing term appear in the same perturbation equation as the nonlinear effect. Substituting Equation (5) into Equation (4) and scaling the dissipative terms, Equation (4) can be modified as: ..
The approximate solution of Equation (6) can be written in the following form: where T n = ε n t, (n = 0, 1, 2). Substituting Equation (7) into Equation (6) and equating the coefficients of like powers of ε, the following partial differential equations can be obtained: where D n = ∂ ∂T n , n = (0, 1, 2). The general solution of Equation (8) can be written as: Here, it is convenient to express A 11 and A 21 in the polar form: where a 1 and a 2 indicate the amplitudes of the first-order vibration mode and the second-order vibration mode in the drive direction, respectively. θ 1 and θ 2 indicate the initial phases of the first-order vibration mode and the second-order vibration mode in the drive direction, respectively.
Next, we substituted Equations (10) and (11) into Equation (9). According to the solvability condition without secular terms, the average equations of amplitudes and phases can be obtained: where φ = θ 2 − θ 1 .
To determine the stability of the periodic solution, the Jacobian matrix of Equation (12) at (a 10 , θ 10 , a 20 , φ 0 ) is given: where  The system is stable when all the matrix eigenvalues are negative. Otherwise, the system is unstable.
. φ = 0, the frequency response equations about a 1 and a 2 can be derived as follows: Using the Newton iteration method and pseudo-arc length method to solve nonlinear coupled Equations (14) and (15), and amplitude frequency response of a 1 and a 2 can be obtained. By discussing the Equations (14) and (15), the influence of the system parameters on the amplitude of the drive response can be obtained.
To obtain the detection output response under the driving nonlinear stiffness, the complex exponential method (CEM) was used to solve the dynamic equation in the sense direction. Equation (2) can be expressed as: M ..
where M is the mass matrix of micro gyro in sense direction, C is the damping matrix, K is the stiffness matrix, and F c (t) is the column vector of the excitation force. According to condition K − ω 2 M = 0, the resonant frequency of first and second modes in sense direction ω 3 and ω 4 can be sought out. Here, From the previous calculation, the approximate solution of x 2 is a 2 cos(ω 0 t−θ 2 ). Therefore, the Coriolis force can be expressed as 2m 2 Ω z a 2 ω 0 sin(ω 0 t−θ 2 ). f c is introduced as the amplitude of Coriolis force, so f c = 2m 2 Ω z a 2 ω 0 . Obviously, the amplitude of Coriolis force is related to the drive excitation frequency. The CEM was used to solve and determine the steady-state response of the system. The steady-state response of Equation (2) can be written as: where b 1 , b 2 are complex amplitudes. Substituting Equation (16) into Equation (2), the following can be obtained: Therefore, the complex amplitude can be obtained as: where b 1 , b 2 and ϑ 1 , ϑ 2 are the amplitudes and initial phases of the sense modal steady-state response, respectively. Finally, the steady-state amplitudes in sense direction can be derived as:

Stiffness Nonlinear Analysis
The calculation parameters were selected as follows: Excitation force amplitude F = 3 × 10 −5 N, ε = 1. From the parameters in Table 1, the natural frequencies of the first-order and second-order drive Micromachines 2019, 10, 578 9 of 22 modes can be obtained as ω 1 = 30904.11 rad/s, ω 2 = 31880.86 rad/s. In the case of σ 2 = 5.91 × 10 7 , the effect of internal resonance parameter changes on micro gyroscope response were not considered in this section. The amplitude-frequency response curves of the first-order and second-order drive modes with different nonlinear stiffness coefficients are shown in Figures 3 and 4. In order to verify the theoretical analytical solution obtained by MMS, the Runge-Kutta numerical method was used to solve Equation (1). Comparison to the theoretical solution, there is a good agreement.

Stiffness Nonlinear Analysis
The calculation parameters were selected as follows: Excitation force amplitude F = 3 × 10 −5 N, ε = 1. From the parameters in Table 1, the natural frequencies of the first-order and second-order drive modes can be obtained as ω1 = 30904.11 rad/s, ω2 = 31880.86 rad/s. In the case of σ2 = 5.91 × 10 7 , the effect of internal resonance parameter changes on micro gyroscope response were not considered in this section. The amplitude-frequency response curves of the first-order and second-order drive modes with different nonlinear stiffness coefficients are shown in Figure 3 and Figure 4. In order to verify the theoretical analytical solution obtained by MMS, the Runge-Kutta numerical method was used to solve Equation (1). Comparison to the theoretical solution, there is a good agreement. As shown in Figure 3, the response of the first-order drive mode was basically the same as that of the linear system, when the nonlinear stiffness coefficients K1 and K2 are less than 10 11 N/m 3 . In this case, the influence of stiffness nonlinearity was negligible. The obvious nonlinear hardening characteristics began to appear at the second peak value, when the nonlinear stiffness coefficients K1 and K2 reached 10 11.5 N/m 3 . The sensitivity of the first peak decreased slightly, but the second peak increased considerably. The nonlinear elastic force was only 2.24% of the linear elastic force in this case. Typical nonlinear characteristics, such as multiple steady state solution, amplitude jump, and frequency offset appeared at the second peak, when the nonlinear stiffness coefficients K1 and K2 reached 10 12 N/m 3 . In comparison to the case of the linear stiffness, the sensitivity at the first peak was reduced by 25.5%, and the first-order resonance frequency was slightly shifted to the right. The sensitivity at the second peak was increased by 77.8% and the second-order resonance frequency was obviously shifted to the right, but the sensitivity was reduced by 40% at the original resonance frequency. At this time, the nonlinear elastic force was 5.42% of the linear elastic force. The sensitivity stability near the second resonance frequency was destroyed, and there was another stable solution  The response of second-order drive mode is shown in Figure 4. The amplitude-frequency response was basically the same as that under linear stiffness, when the nonlinear stiffness coefficients K1 and K2 were less than 10 11 N/m 3 . The influence of stiffness nonlinearity was negligible in this case. The amplitude-frequency response began to appear hardened and the sensitivity was basically unchanged. The second-order resonant frequency was slightly offset to the right when the nonlinear stiffness coefficients K1 and K2 reached 10 11.5 N/m 3 . The response curves showed obvious nonlinear characteristics when the nonlinear stiffness coefficients K1 and K2 reached 10 12 N/m 3 . Nonlinear behaviors, such as multiple steady state solutions and amplitude jump, appeared near the second formant. The BW area in Figure 4 represents the 3dB bandwidth of the drive output response. With the increase of nonlinear stiffness, the output response bandwidth remained basically unchanged, but the frequency corresponding to the bandwidth shifted slightly to the right due to the nonlinear frequency offset effect. Within the bandwidth range, the sensitivity of the output response decreased slightly with nonlinear enhancement. However, it is important that the multiple steady state solution region caused by the hardening of the stiffness nonlinearity was outside the bandwidth range, as the system instability caused by the amplitude jump did not occur in the bandwidth range. This is because that the 2-DOF drive-mode structure of the micro-gyro model had a relatively flat and high sensitivity between the two peaks in the amplitude-frequency curve, which made the drive As shown in Figure 3, the response of the first-order drive mode was basically the same as that of the linear system, when the nonlinear stiffness coefficients K 1 and K 2 are less than 10 11 N/m 3 . In this case, the influence of stiffness nonlinearity was negligible. The obvious nonlinear hardening characteristics began to appear at the second peak value, when the nonlinear stiffness coefficients K 1 and K 2 reached 10 11.5 N/m 3 . The sensitivity of the first peak decreased slightly, but the second peak increased considerably. The nonlinear elastic force was only 2.24% of the linear elastic force in this case. Typical nonlinear characteristics, such as multiple steady state solution, amplitude jump, and frequency offset appeared at the second peak, when the nonlinear stiffness coefficients K 1 and K 2 reached 10 12 N/m 3 . In comparison to the case of the linear stiffness, the sensitivity at the first peak was reduced by 25.5%, and the first-order resonance frequency was slightly shifted to the right. The sensitivity at the second peak was increased by 77.8% and the second-order resonance frequency was obviously shifted to the right, but the sensitivity was reduced by 40% at the original resonance frequency. At this time, the nonlinear elastic force was 5.42% of the linear elastic force. The sensitivity stability near the second resonance frequency was destroyed, and there was another stable solution far below the peak value at the instantaneous natural frequency corresponding to the second peak value, because of the amplitude jumping behavior and the appearance of multiple steady state solutions. This phenomenon was due to the dependence of the nonlinear system on the initial conditions. The system was a periodic motion of smaller amplitude when the operating frequency was lowered from the high frequency to the resonant frequency, while the system was a periodic motion of larger amplitude when the operating frequency rose from the low frequency to the resonant frequency. As shown in Figure 3d, the dashed line was an unstable region of the approximate periodic solution and the intermediate solution branch of the multiple steady state solution. The motion corresponding to this area is unlikely to occur in a real micro gyroscope system.
The response of second-order drive mode is shown in Figure 4. The amplitude-frequency response was basically the same as that under linear stiffness, when the nonlinear stiffness coefficients K 1 and K 2 were less than 10 11 N/m 3 . The influence of stiffness nonlinearity was negligible in this case. The amplitude-frequency response began to appear hardened and the sensitivity was basically unchanged. The second-order resonant frequency was slightly offset to the right when the nonlinear stiffness coefficients K 1 and K 2 reached 10 11.5 N/m 3 . The response curves showed obvious nonlinear characteristics when the nonlinear stiffness coefficients K 1 and K 2 reached 10 12 N/m 3 . Nonlinear behaviors, such as multiple steady state solutions and amplitude jump, appeared near the second formant. The BW area in Figure 4 represents the 3dB bandwidth of the drive output response. With the increase of nonlinear stiffness, the output response bandwidth remained basically unchanged, but the frequency corresponding to the bandwidth shifted slightly to the right due to the nonlinear frequency offset effect. Within the bandwidth range, the sensitivity of the output response decreased slightly with nonlinear enhancement. However, it is important that the multiple steady state solution region caused by the hardening of the stiffness nonlinearity was outside the bandwidth range, as the system instability caused by the amplitude jump did not occur in the bandwidth range. This is because that the 2-DOF drive-mode structure of the micro-gyro model had a relatively flat and high sensitivity between the two peaks in the amplitude-frequency curve, which made the drive module have high robustness in the working process. In addition, the drive module was not easily affected in the bandwidth region.
Combining Equations (14) and (15), and (20), the relationship between the steady-state amplitude of the sense mode and the driving frequency can be obtained, that is, the amplitude frequency response relationship of the sense mode. The mathematical analysis software Mathematica was used to plot the amplitude-frequency response curve of the primary resonance of the first-order and second-order sense modes, as shown in Figures 5 and 6. With Runge-Kutta method used to solve Equations (1) and (2), a series of numerical solutions which coincide with the theoretical analytical solution were obtained, and the correctness of theoretical solution was verified.
As illustrated in Figures 5 and 6, the amplitude frequency response of the first-order and second-order sense modes were similar to that of the drive modes. The response curve was basically the same as that of linear stiffness case when the nonlinear stiffness coefficients are less than 10 11 N/m 3 . The amplitude-frequency response began to show nonlinear characteristics, and a slight increase would cause a significant change in the shape of the response curve when the stiffness nonlinearity coefficients exceeded 10 11 N/m 3 . As shown in Figure 5c,d and Figure 6c,d, the steady-state response had typical nonlinear behaviors, such as multiple steady state solution, amplitude jump, and frequency offset. The nonlinear influence was stronger, which can be seen from the bending degree of the peak. The difference is that when the stiffness nonlinearity coefficients increased from 10 11.5 N/m 3 to 10 12 N/m 3 , the frequency band range corresponding to the unstable solution (dashed line) of the first-order sense mode increased, while second-order sense mode decreased. The reduction of the corresponding frequency band range of the unstable solution means that the probability of amplitude jump in steady-state response decreased, that is to say, the possibility of an abrupt change of output signal in the micro-gyroscope system was reduced. As shown in Figure 6, the area marked by BW is the bandwidth for sensing the output response, and it is also the detection bandwidth of the micro gyroscope. It can be seen from Figure 6 that the bandwidth was slightly shifted to the right when the nonlinearity was strengthened. The bandwidth was narrowed and the detection sensitivity decreased in the bandwidth range when the nonlinearity was relatively strong. On the other hand, the region of the multi-stable solution was outside the detection bandwidth and did not affect the stability in the bandwidth range because of the nonlinear stiffness hardening. It is similar to the drive output response. Therefore, the instability behavior, such as amplitude jump, will not occur when the micro gyroscope works in the detection bandwidth, even if there is a strong nonlinear stiffness. This is also due to the high robustness of the sense module. The amplitude in the bandwidth is relatively flat and the stiffness nonlinearity has less effect on the inner part of the bandwidth. In addition, the frequency band region of the drive output response bandwidth is highly matched with the sense response bandwidth region, so that the micro gyroscope has a strong resistance to nonlinear factors when operating in the working bandwidth. module have high robustness in the working process. In addition, the drive module was not easily affected in the bandwidth region. Combining Equations (14), (15), and (20), the relationship between the steady-state amplitude of the sense mode and the driving frequency can be obtained, that is, the amplitude frequency response relationship of the sense mode. The mathematical analysis software Mathematica was used to plot the amplitude-frequency response curve of the primary resonance of the first-order and second-order sense modes, as shown in Figure 5 and Figure 6. With Runge-Kutta method used to solve Equations (1) and (2), a series of numerical solutions which coincide with the theoretical analytical solution were obtained, and the correctness of theoretical solution was verified.  As illustrated in Figure 5 and Figure 6, the amplitude frequency response of the first-order and second-order sense modes were similar to that of the drive modes. The response curve was basically the same as that of linear stiffness case when the nonlinear stiffness coefficients are less than 10 11 N/m 3 . The amplitude-frequency response began to show nonlinear characteristics, and a slight increase would cause a significant change in the shape of the response curve when the stiffness nonlinearity coefficients exceeded 10 11 N/m 3 . As shown in Figure 5c, 5d, 6c, and 6d, the steady-state response had typical nonlinear behaviors, such as multiple steady state solution, amplitude jump, and frequency offset. The nonlinear influence was stronger, which can be seen from the bending degree of the peak. The difference is that when the stiffness nonlinearity coefficients increased from 10 11.5 N/m 3 to 10 12 N/m 3 , the frequency band range corresponding to the unstable solution (dashed line) of the first-order sense mode increased, while second-order sense mode decreased. The reduction of the corresponding frequency band range of the unstable solution means that the probability of amplitude jump in steady-state response decreased, that is to say, the possibility of an abrupt change of output signal in the micro-gyroscope system was reduced. As shown in Figure 6, the area marked by BW is the bandwidth for sensing the output response, and it is also the detection bandwidth of the micro gyroscope. It can be seen from Figure 6 that the bandwidth was slightly shifted to the right when the nonlinearity was strengthened. The bandwidth was narrowed and the detection sensitivity decreased in the bandwidth range when the nonlinearity was relatively strong.

System Parameters Analysis
The parameter σ 2 is an internal resonance tuning parameter, which characterizes the proximity of the resonant frequencies of the first-order and second-order drive modes. σ 2 = ω 2 2 − ω 2 1 , when ε = 1. Therefore, both the internal resonance tuning parameters σ 2 and the steady-state response of the system will change when a certain order drive mode resonance frequency shifts under the influence of external conditions. The variation curves of the response amplitude of each degree of freedom of the system are shown in Figure 7 with internal resonance parameter σ 2 at different operating frequency ω 0 .
As shown in Figure 7, the operating frequency ω 0 takes three different values (30940.11 rad/s, 31500 rad/s, 32000 rad/s), corresponding to three different values of σ 1 (0, 3.5 × 10 7 , 6.67 × 10 7 ), which represent the complete tuning, small detuning, and detuning of the primary resonance, respectively. In Figure 7a, with σ 2 increases the response amplitude of the first-order drive mode decreased from the stable value to the minimum value, and then slowly increased to another stability value, when σ 1 takes different values. However, in Figure 7b, the response amplitude of second-order drive mode increased first, and then decreased gradually after reaching the maximum value with the increase of σ 2 . This is because that ω 1 was closer to ω 2 , and the internal resonance was closer to the tuning state when σ 2 was relatively small. The partial energy generated by the first-order drive mode primary resonance was transferred to the second-order drive mode, and the response amplitude of the second-order drive mode increased gradually, while the amplitude of the first-order drive mode decreased. The degree of internal resonance detuning became larger, and energy transfer was relatively less with the larger σ 2 . The response amplitude of the second-order drive mode gradually decreased, but the first-order drive mode slowly increased.
bandwidth is relatively flat and the stiffness nonlinearity has less effect on the inner part of the bandwidth. In addition, the frequency band region of the drive output response bandwidth is highly matched with the sense response bandwidth region, so that the micro gyroscope has a strong resistance to nonlinear factors when operating in the working bandwidth.

System Parameters Analysis
The parameter σ2 is an internal resonance tuning parameter, which characterizes the proximity of the resonant frequencies of the first-order and second-order drive modes.    , when ε = 1. Therefore, both the internal resonance tuning parameters σ2 and the steady-state response of the system will change when a certain order drive mode resonance frequency shifts under the influence of external conditions. The variation curves of the response amplitude of each degree of freedom of the system are shown in Figure 7 with internal resonance parameter σ2 at different operating frequency ω0.  As shown in Figure 7c,d, in the case of three different σ 1 , the response amplitudes of the first-order and second-order sense modes increased first and then decreased with the increase of σ 2 . Its change trend is similar to that of the second-order drive mode in Figure 7b. This is because the response amplitude of the second-order drive mode in this system was the output response of the drive direction, which determined the amplitude of the Coriolis force in the sense direction. There is no nonlinear stiffness in the sense direction. The steady-state response amplitudes of the first-order and second-order sense modes are proportional to the Coriolis force amplitude, and the change trend is consistent with that of the second-order drive mode. During the operation of the micro gyroscope, the stiffness of elastic micro-beams decreased because of the softening effect of the elastic modulus when the external temperature rose. Then, the modal resonance frequency was reduced and the drift was generated. Frequency drift results in the change of internal resonance parameters, even at constant excitation frequency. The response amplitude will change or even change greatly. In addition, the stiffness of the elastic micro-beams will be affected when there are errors(deviations) in the processing of them. It will cause the resonance frequency to drift, as well as have a greater impact on the amplitude of the output response. Therefore, the external working conditions and processing errors may have a large impact on the output of the gyroscope system with 1:1 internal resonance. In this case, the error factors should be fully considered, and corresponding error compensation should be made.
In order to study the influence of external input energy on the driving response of gyroscope system more clearly, AC voltage with different amplitudes was applied to both ends of driving electrodes to subject the driving mass to the electrostatic driving force of different amplitudes. In addition, the value of parameter σ 2 was selected in the region where the amplitude existed in the multiple steady-state solutions in Figure 7. Here, σ 2 = 10 × 10 7 , K 1 = K 2 = 10 12 N/m 3 were selected. The relationship between the steady-state response amplitude of the first-order and second-order drive modes and electrostatic driving force amplitude were studied under the condition of primary resonance complete tuning (σ 1 = 0) and detuning (σ 1 = 6.67 × 10 7 ), as shown in Figure 8. It can be seen from Figure 8a that both the response amplitudes of first-order and second-order drive modes were uniquely determined. The driving mass and the decoupled frame were both stable periodic motions, under the condition of primary resonance complete tuning, as the electrostatic driving force amplitude increased. Both of the amplitudes of the first-order and second-order drive modes increased with the increase of the electrostatic force amplitude, and the growth relationship was nonlinear. The growth rate decreased gradually when the electrostatic force amplitude was large. The amplitude of the second-order drive mode reached saturation when the electrostatic force amplitude reached a certain threshold. In Figure 8b, the primary resonance was detuned, and the steady-state amplitudes of the first-order and second-order drive modes simultaneously exhibited a multiple stable solution phenomenon in the range of AB. Two of the three solutions were stable, and the initial conditions determined which solution was the real response. There was only one stable solution in other regions, which ws confirmed by numerical calculation. The amplitudes of the firstorder and second-order drive modes were approximately linear with the electrostatic force amplitude when the electrostatic force amplitude was below the corresponding amplitude of point B. The amplitudes of the first-order and second-order drive modes jumped simultaneously when the electrostatic force amplitude reached point B. With jumping from the current amplitude to another higher stable amplitude, the amplitude of the second-order drive mode gradually reached saturation.
The first-order and second-order drive modes will produce vibrations regardless of the It can be seen from Figure 8a that both the response amplitudes of first-order and second-order drive modes were uniquely determined. The driving mass and the decoupled frame were both stable periodic motions, under the condition of primary resonance complete tuning, as the electrostatic driving force amplitude increased. Both of the amplitudes of the first-order and second-order drive modes increased with the increase of the electrostatic force amplitude, and the growth relationship was nonlinear. The growth rate decreased gradually when the electrostatic force amplitude was large. The amplitude of the second-order drive mode reached saturation when the electrostatic force amplitude reached a certain threshold. In Figure 8b, the primary resonance was detuned, and the steady-state amplitudes of the first-order and second-order drive modes simultaneously exhibited a multiple stable solution phenomenon in the range of AB. Two of the three solutions were stable, and the initial conditions determined which solution was the real response. There was only one stable solution in other regions, which ws confirmed by numerical calculation. The amplitudes of the first-order and second-order drive modes were approximately linear with the electrostatic force amplitude when the electrostatic force amplitude was below the corresponding amplitude of point B. The amplitudes of the first-order and second-order drive modes jumped simultaneously when the electrostatic force amplitude reached point B. With jumping from the current amplitude to another higher stable amplitude, the amplitude of the second-order drive mode gradually reached saturation.
The first-order and second-order drive modes will produce vibrations regardless of the magnitude of the electrostatic driving force because of the stiffness of the first-order and second-order drive modes coupled with each other in structure. At the beginning, amplitudes grows linearly with the increase of electrostatic force amplitude. Then, the increase of driving amplitudes becomes nonlinear, and the phenomenon of the amplitude jump occurs when the primary resonance detuning. The amplitudes of the first-order and second-order drive modes increase slower with the increase in electrostatic force amplitude when it is higher than that of jumping region. Then, the amplitude of the second-order drive mode gradually reaches the saturation state. The response of the second-order drive mode is the output response of the driving direction, which determines the amplitudes of Coriolis force and sense mode responses. Therefore, the electrostatic force amplitude should be as large as possible when it is far away from the unstable region in order to make the output response of the gyroscope as large and stable as possible. At the same time, it should be noted that the driving voltage should not exceed the pull-in voltage of the comb electrodes.

Local Bifurcation Analysis
The local bifurcation response of the drive module was discussed for the further study of the influence of external disturbance on the dynamic performance of nonlinear gyroscope system. Here, the transformation of a 2 2 = Z, σ 2 = µ was done, and the recognition condition of the theory of singularity was used to simplify Equations (14) and (15). The local bifurcation equation of the drive direction can be obtained: where . The specific parameters are given in Appendix A.
According to the theory of singularity, the transition sets in the unfolding parameter space is composed of Σ = B ∪ H ∪ D. Here, B is the bifurcation set, H is the hysteresis set, and D is the double limit point set. The specific expressions of the transition set are given in Appendix B.
Since the parameter space is a seven-dimensional space and the transition set cannot be visually expressed in space, the bifurcation behavior in the parameter space on the projection plane was mainly discussed. Here, g 6 , g 7 were selected as the unfolding parameters. The transition sets and bifurcation diagrams on the g 6 -g 7 plane were obtained, as shown in Figure 9. The local bifurcation response of the drive module was discussed for the further study of the influence of external disturbance on the dynamic performance of nonlinear gyroscope system. Here, the transformation of 2 2 aZ  , σ2 = μ was done, and the recognition condition of the theory of singularity was used to simplify Equations (14) and (15). The local bifurcation equation of the drive direction can be obtained: According to the theory of singularity, the transition sets in the unfolding parameter space is composed of B H D     . Here, B is the bifurcation set, H is the hysteresis set, and D is the double limit point set. The specific expressions of the transition set are given in Appendix B.
Since the parameter space is a seven-dimensional space and the transition set cannot be visually expressed in space, the bifurcation behavior in the parameter space on the projection plane was mainly discussed. Here, g6, g7 were selected as the unfolding parameters. The transition sets and bifurcation diagrams on the g6-g7 plane were obtained, as shown in Figure 9. On the projection plane g6-g7, the curves represented by the bifurcation set B and the hysteresis set H divide the plane into five sub-intervals. The bifurcation diagrams corresponding to two points on each sub-interval are topological equivalent, that is, persistent. Any two points that are not in the On the projection plane g 6 -g 7 , the curves represented by the bifurcation set B and the hysteresis set H divide the plane into five sub-intervals. The bifurcation diagrams corresponding to two points on each sub-interval are topological equivalent, that is, persistent. Any two points that are not in the same sub-interval are topologically inequivalent. Taking a random point in each sub-interval, the retentive bifurcation diagram in each region was obtained, as shown in Figure 9 I to V. The system presents a more complicated bifurcation behavior with the change of the unfolding parameters g 6 , g 7 .
There is a unique parameter Z corresponding to the determined parameter µ, and the system is stable in this condition in the interval I and II. There are many values of Z corresponding to the determined parameter µ, and the phenomenon of multiple solutions occurs within intervals III, IV, and V. The system exhibits bifurcation behavior and periodic motion become instable in this moment. In interval III, parameter Z will change from the single solution region to the multi-solution region, when parameter µ changes from the value far from 0 to the vicinity of 0. In this case, the motion behavior of the system will change greatly, that is, the instability region is near 0. In intervals IV and V, parameter Z starts with only a unique solution, then reaches the first multi-solution region, becomes the only solution, passes through the second multi-solution region, and finally becomes the stable solution, with the µ becoming larger gradually. At present, there are two unstable regions in the system. In summary, the dynamic response of the system changes from the initial stable state to the unstable state, in which unstable regions exist, when the unfolding parameters g 6 , g 7 caused by external disturbance change from interval I or II to interval III, IV, or V. The movement of the gyroscope changes essentially, which results in the system not working properly.
It can be seen from Appendix B that the parameters h 7 , h 8 can be converted into the following form: It can be seen from Equation (22) that almost every item in parameter h 7 contains factorβ 4 1 (β 2 −ω 2 2 ). Therefore, the value of h 7 is mainly affected byβ 1 ,β 2 , and the stiffness coefficients k 2 and k 3 play major roles in the disturbance of parameter g 7 . The common factor in parameter h 8 isβ 4 1 f 2 , so it can be seen that the stiffness coefficient k 2 and the electrostatic driving force amplitude F play major roles in the disturbance of parameter g 8 . To keep the system as robust as possible to the environmental changes and avoid the occurrence of multiple solutions, the parameters g 6 , g 7 should be taken in interval I or II. For this reason, the value of parameters g 6 , g 7 could be increased as much as possible and kept at a positive value. That is to say, the stiffness coefficients k 2 and k 3 of the micro-beams should be increased to ensure the value of g 6 . To ensure the value of g 7 , the electrostatic driving force amplitude F and the damping coefficient c 2 should be increased.
The unfolding parameters were used here to characterize the disturbance of the environment. Here, the environmental disturbance is a general concept which contains many factors, such as temperature, electrical, external force, and noise. Temperature noise, electrical noise, and phase noise [43,44] in noise mixing will affect the stability of the system from the system parameters or excitation signals. By adjusting the system parameters k 2 , k 3 , c 2 , F, the system could be placed in a relatively stable area as far as possible. From a qualitative point of view, it can suppress the influence of noise mixing.

Conclusions
In this paper, a 4-DOF micro gyroscope model including double drive-mode and double sense-mode was investigated. Particularly, the complex dynamic characteristics of the micro gyroscope were studied by considering the cubic stiffness nonlinearity, which may occur in both degree-of-freedom in the drive direction. The complex dynamic behaviors of the highly robust micro gyroscope system under nonlinear influence were obtained. The key conclusions are as follows: (1) The dynamic characteristics of the drive and sense modes are basically unaffected when the driving nonlinearity is relatively weak, which is similar to that of the linear design. However, the dynamic characteristics of the drive and sense modes will be obviously affected when the driving nonlinearity is strong. There are obvious amplitude jumps, multiple steady state solutions, and resonance frequency offsets in the amplitude-frequency response curves. Because the bandwidth regions of drive and sense modes match each other and have high robustness, the nonlinear hardening behavior in second-order drive mode and sense modes is no longer typical right bending but is instead right-down bending. Amplitude jumps and unstable solutions occur outside the response bandwidth of drive and sense modes and have only a slight effect on the sensitivity within the bandwidth. (2) When the internal resonance is close to the tuning state, part of the energy generated by the primary resonance of the first-order drive mode is transferred to the second-order drive mode. The response amplitude of the second-order drive mode increases gradually, while that of the first-order drive mode decreases. When the degree of internal resonance detuning increases, the energy transfer is relatively small, the response amplitude of the second-order drive mode decreases gradually, and the first-order drive mode increases slowly. The frequency drift caused by the change of the external environment will change the internal resonance parameter of the system. Even at a constant excitation frequency, the response amplitude will change or even change greatly. Therefore, external working conditions and processing errors should be fully considered, and corresponding error compensation should be made. (3) When the primary resonance complete tuning, the motions of 2-DOF drive mode are stable periodic motions with different electrostatic force amplitudes. When the main resonance is detuned, the response amplitudes of the drive modes will jump, and the multi-stable solution will appear with the increase of the electrostatic force amplitude and reach saturation when the electrostatic force amplitude increases to a certain value. In order to satisfy the stable output response of the gyroscope as much as possible, the electrostatic force amplitude should be as large as possible when it is far away from the unstable region. At the same time, it should be noted that the driving voltage must not exceed the pull-in voltage of the comb electrodes. (4) Through the study of the local bifurcation of the nonlinear gyroscope system, it is found that the stiffness coefficient of the micro-beams connected with the decoupling frame and the electrostatic driving force amplitude play major roles in the parameter perturbation. In order to keep the system as robust as possible to the environment and avoid the occurrence of multiple solutions phenomena, the stiffness coefficients of the micro-beams on both sides of the decoupling frame should be increased, and the damping coefficient between the driving mass and the decoupled frame and the electrostatic driving force amplitude should be appropriately increased.