# Bursting Oscillations and Sliding Motion in Permanent Magnet Synchronous Motor Systems with Friction Factor

Authors: Rui Qu; Shaolong Li
Journal: Shock and Vibration
 (2021)

---

## Abstract

The main purpose of this paper is to investigate the qualitative effects of external excitation and friction factor on the response of permanent magnet synchronous motor (PMSM) system. Three different modes of bursting oscillations are found. In particular, the introduction of friction function changes the governing equations from a smooth type to a nonsmooth (Filippov) type in which the special sliding motion is observed. The mechanism, attractor structure, vector field structure, and analytic bifurcation conditions of bursting oscillation and sliding motion are discussed in detail. The validity of theoretical results obtained is verified by numerical simulations and analysis.


---

## Body

## 1. Introduction

Due to the high efficiency and small size, permanent magnet synchronous motor (PMSM) has a wide range of engineering applications [1, 2] such as vehicle drive, medical machinery, and servo robot. As a typical nonlinear, multiple-variable, and strongly coupled system, PMSM system has complicated nonlinear behaviors such as bifurcations, chaos, and bursting oscillations [3, 4] with the variation of external excitation and internal parameters. These special dynamical behaviors may affect the stability of system and cause equipment failure.In recent years, the hot research topics of nonlinear dynamic systems mainly involve high-dimensional systems, high codimension bifurcations [5], nonlinear systems with special structures [6, 7], and others. For the high-dimensional systems, their dynamical behaviors are bound to be complicated due to the greater number of state variables. The observation of hidden hyperchaotic attractors [8] in different types of high-dimensional systems [9] is a good proof of their complexity. Although most practical systems are of high dimension, the complexity caused by extra state variables will make it difficult to analyze the coupling effects between nonlinear factors such as excitation and friction. Therefore, this paper employs a three-phase PMSM system as the object to investigate the coupling of external excitation and friction factors on the system which has been rarely reported before. To be specific, the studies involving different types of oscillations resulting from external factors in PMSM system are still relatively lacking. Besides, the great majority of current researches on oscillations are still limited in smooth system [10, 11], which indicates that the possible effects resulting from the introduction of nonsmooth factor [12, 13] are still not clear. In engineering practice, the vibration with friction corresponding to the oscillation with nonsmooth bifurcation may aggravate the equipment damage. Therefore, it is of great theoretical and practical significance to further study and control the complex dynamical behaviors of PMSM system.In this paper, the stability and bifurcation conditions of equilibrium points are derived in detail. The two-parameter bifurcation diagram, Poincaré map, and Lyapunov exponent spectrum are employed to analyze the inherent properties of autonomous PMSM system. Meanwhile, the nonsmooth bursting oscillations of nonautonomous PMSM system with friction factor are studied through equilibrium curves and transformation phase portraits (TPP) [14, 15]. In addition, the mechanism of sliding motion is discussed separately, and the theoretical results obtained from attractor and vector field structure analysis [16, 17] are also verified by the calculation of frequency and sliding time length. The physical meanings of different modes of vibrations are also explained.


## 2. Mathematical Model

The structural sketch and the physical photo of PMSM system are shown in Figure1, and the three-dimensional nonautonomous mathematical model of PMSM system considering external excitation and friction factor can be expressed as(1)diddt=ud−R1id+ωLqiqLd,diqdt=ud−R1iq−ωLdid+ωψfLq,dωdt=npψfiq+npLd−Lqidiq−TL−βω+fωJ.Figure 1
(a) Structural sketch and (b) physical photo of PMSM.
(a)(b)The physical meanings of state variables and parameters in the governing equation (1) are as follows:id is stator current of d-axis; iq is stator current of q-axis; ω is mechanical angular velocity of rotorud is voltage of d-axis; uq is voltage of q-axis; Ld is inductance of d-axis; Lq is inductance of q-axisTL is external torque; ψf is permanent magnetic flux; R1 is stator resistance; np is polar logarithmβ is viscous damping coefficient; J is moment of inertia; fω is friction functionThe control equations of PMSM system (1) can be derived from two parts.


### 2.1. Electromagnetic Part

The voltage equations for the two-phase rotating coordinate model of thed-axis and q-axis can be expressed as(2)ud=Rid+ddtψd−ωψq,uq=Riq+ddtψq+ωψd.where ψd and ψq represent the d-axis and q-axis components of the stator magnetic chain, respectively. Besides, the equations of the stator magnetic chain can be expressed as(3)ψd=Ldid+ψf,ψq=Lqiq.Combining the voltage and the stator magnetic chain equations, it can be concluded that(4)diddt=ud−R1id+ωLqiqLd,diqdt=ud−R1iq−ωLdid+ωψfLq.


### 2.2. Mechanical Motion Part

The mechanical motion equation of the motor can be obtained by Newton’s second law.(5)Jα=Jdωdt,=Te−TL−βω+fω,where the electromagnetic torque Te can be expressed as Te=npiqψf+Ld−Lqid based on the electromechanical energy conversion principle (the electric angle is equal to the product of the mechanical angle and the polar logarithm). Combining the electromagnetic part (4) and mechanical motion part (5), the control equations of three-phase PMSM system can be formed as system (1).In order to simplify the governing equations, the time scale transformation and linear affine transformation are introduced as(6)t=τ⋅t^,x=λx^,where τ=Lq/R1,x=id,iq,ωT,x^=i^d,i^q,ω^T, and(7)λ=LqβLdnpτψf000βnpτψf0001τ.Through the assumption of uniform air gapLd=Lq and dimensionless process,(8)γ=npψf2R1β,σ=LqβR1J,u^q=npψfLquqR12β,u^d=npψfLqudR12β,ε=Lqβ2Ld−LqLdJnpψf2,T^L=Lq2TLR12J,f^ω^=Lq2fω^R12J, and the governing equations can be written as(9)di^ddt^=−LqLdi^d+ω^iq+u^d,di^qdt^=−i^q−ω^i^d+γω^+u^q,dω^dt^=σi^q−ω^+εi^di^q−T^L+f^ω^,⇒Lq=Lddi^ddt^=−i^d+ω^iq+u^d,di^qdt^=−i^q−ω^i^d+γω^+u^q,dω^dt^=σi^q−ω^−T^L+f^ω^,where σ>0,γ>0.


### 2.3. Stability and Bifurcations of the Autonomous System

First of all, it is necessary to briefly investigate the inherent properties of PMSM system, so it can be assumed thatu^d=u^q=0,T^L=f^ω^=0 to further simplify the control equations as(10)di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^.The equilibrium point of the system satisfies the following conditions:(11)−i^d0+ω^0i^q0=0,−i^q0−ω^0i^d0+γω^0=0,σi^q0−ω^0=0,⇒i^d0=ω^02,i^q0=ω^0,ω^0ω^02−γ−1γ−1=0.which indicates that there are three equilibrium points in the system with γ>1:(12)E00,0,0,E+γ−1,γ−1,γ−1,E−γ−1,−γ−1,−γ−1,while E0 is the sole equilibrium with γ≤1.Routh-Hurwitz criterion can be employed to analyze the stability of equilibrium point. The Jacobi matrix of system (10) can be expressed as(13)J=−1ω^0i^q0−ω^0−1−i^d0+γ0σ−σ.It can be seen that all equilibrium points of the system can be expressed in a unified formω^02,ω^0,ω^0. Therefore, the characteristic equation at the equilibrium point can be formed as(14)λ3+σ+2σ2+1+σω^02−σγ+ω^02+2σλ+3σω^02−σγ+σ=0.Based on the inference of Routh–Hurwitz stability criterion, for a cubic polynomial(15)P3z=z3+a1z2+a2z3+a3,its Routh-Hurwitz conditions can be formed as(16)a1>0,a11a3a2>0,a3a1a2−a3>0,⇕a1>0,a3>0,a1a2>a3.Therefore, the stability conditions of the equilibrium point can be expressed as(17)σ>−2,σω^02−γ+1>0,σ+21+σω^0−σγ+ω^02+2σ−σω^02−γ+1>0.When conditiona3>0 degenerates to a3=0, the characteristic equation has a zero root corresponding to the Fold bifurcation conditions.(18)FB:σ>−2,γ=1,ω^0=0,σ+21+σ>0,⇒σ>0γ=1,ω^0=0.When conditiona1a2>a3 degenerates to a1a2=a3, a pair of pure imaginary roots appear in the solution of characteristic polynomial (14) indicating the occurrence of Hopf bifurcation. Therefore, the conditions of Hopf bifurcation are as follows:(19)HB:σ>−2,σω^02−γ+1>0,σ+21+σω^0−σγ+ω^02+2σ−σω^02+γ−1=0,⟹σ>0,γ>1,ω^02=γ−1σ2−σγ+2γ+4σ=0.Based on above analysis, the two-parameter bifurcation diagram onγ,σ plane is shown in Figure 2 in which FB represents the Fold bifurcation curve, while HB describes the subcritical Hopf bifurcation curve. It can be seen that the parameter plane is divided into three parts by the two bifurcation curves motioned. In region 1, the system has a sole stable equilibrium point. With the variation of parameter combinations, it traverses the Fold curve and enters region 2 leading to the appearance of two new stable equilibrium points. Then, the equilibrium points continue to cross the Hopf curve and then enter region 3, resulting in the existence of three unstable equilibrium points.Figure 2
Two-parameter bifurcation diagram onσ,γ plane.Taking the fixed parameterσ=4.0 and variable parameter γ, the bifurcation diagram of equilibrium points is shown in Figure 3 in which the solid line represents the stable equilibrium point, while the dotted line corresponds to the unstable ones. The specific parameters of Fold and Hopf bifurcations are calculated as γ=1.0 and γ=16.0, respectively. Two unstable limit cycles generated by the subcritical Hopf bifurcation points HB± are denoted by LPC1,2. It can be seen that the property changes of equilibrium points are consistent with Figure 2.Figure 3
Bifurcation diagram of autonomous PMSM system withσ=4.0.


### 2.4. Dynamical Behavior Simulation of the Autonomous System

#### 2.4.1. Case A

Taking a set of fixed parameters in region 1 (Figure2) as σ=4.0,γ=0.1, the system has the sole equilibrium point E0′0,0,0 as shown in Figure 4(a), while, in Figure 4(b), there are two stable focal equilibrium points E±4,±2,±2 and an unstable saddle equilibrium point E0′0,0,0 with another set of fixed parameters σ=4.0,γ=5.0 in region 2.Figure 4
Phase portraits onω^,i^d plane with (a) σ=4.0,γ=0.1 and (b) σ=4.0,γ=5.0.
(a)(b)


#### 2.4.2. Case B

Taking another set of parameters in region 2 asσ=4.0,γ=15.0, the number and stability of equilibrium points remain unchanged, but the convergence modes of trajectory to the stable focus are obviously different due to the variation of initial values. In Figure 5(a), the envelope of the trajectory is similar to the chaotic attractor which can be also confirmed by the quasi-chaotic process in time history (Figure 6(a)). Then, the trajectory suddenly approaches the stable focus by the means of focal-type convergence after the quasi-chaotic process. However, with another set of initial values, the trajectory directly converges to the stable focus as shown in Figures 5(b) and 6(b).Figure 5
Phase portraits with initial valuesi^d,i^q,ω^ of (a) 20,3.74,3.74 and (b) 15,3.74,3.74.
(a)(b)Figure 6
Time histories with initial valuesi^d,i^q,ω^ of (a) 20,3.74,3.74 and (b) 15,3.74,3.74.
(a)(b)


#### 2.4.3. Case C

In the two cases mentioned above, the system can operate stably. Taking the last set of parameters in region 3 asσ=4.0,γ=17.0, it can be found from the phase portrait (Figure 7(a)) and time history (Figure 7(c)) that the chaotic attractor appears after the Hopf bifurcation (σ=4.0,γ=16.0). In addition, the dynamical evolution of the autonomous PMSM system (10) and the attraction domain transition of stable attractors can be also proved by Poincaré map (Figure 7(b)) and the maximum Lyapunov exponent diagram (Figure 7(d)). To be specific, maximum Lyapunov exponent λmax=0.3726>0 at the Hopf bifurcation point, which indicates that the two adjacent phase trajectories are going to separate gradually, making the system enter the chaotic state finally.Figure 7
Chaos of the autonomous system: (a) phase portrait; (b) Poincaré map; (c) time history; (d) the maximum Lyapunov exponent diagram.
(a)(b)(c)(d)


#### 2.4.4. Discussion on Physical Meaning

The equilibrium point in the mathematical model corresponds to the stable working state of the machine in engineering practice. The chaotic solution or quasi-chaotic convergence process of the system represents the irregular oscillation of the equipment, which usually leads to the aging and damage of motor and thus is ought to be avoided. In this section, the self-excited operation of the PMSM system without external load is mainly discussed. The results show that the risk of irregular oscillation can be effectively reduced by proper parameter combination and initial value selection. However, due to the existence of multiple excitations (external excitation, parametric excitation, etc.) in engineering practice, the PMSM system may produce multiple-scale effect inducing the bursting oscillation which is more complex. In addition to the oscillation, the existence of friction factors may further aggravate the equipment damage and even change the inherent properties of the system. These problems will be discussed in the next section.


### 2.5. Stability and Bifurcations of the Nonsmooth PMSM System

#### 2.5.1. Bifurcation Analysis of Subsystems

In order to explore the typical dynamical behaviors of nonsmooth PMSM system, the symbolic functionf^ω^=c⋅signω^=−c,ω^>0+c,ω^<0,c≥0 is employed to describe the dry friction, while the excitation function T^L=AsinΩt^ is used to represent the external torque. It should be noted that the definition of symbolic function value with ω^=0 is omitted for convenience because it depends on the equilibrium equation (similar to static friction). Therefore, the control equations of Filippov-type PMSM can be written as(20)di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L+f^ω^.It should be pointed out that the effect of two-time scales [18] may appear when the excitation frequency Ω is far less than the natural frequency Ω0 which can be estimated by imaginary part of equilibrium point eigenvalue.There will be magnitude differences between external excitation frequency and natural frequency if0<Ω≪1. In any natural frequency period TΩ0=2π/Ω0,t^∈t^0,t^0+TΩ0, the variation range of external excitation is WA,WB, where WA=AsinΩt^0 and WB=AsinΩt0+2πΩ/ω. Based on the assumption Ω/ω≪1⇒WA≈WB, it can be concluded that, in any natural frequency period, the external excitation almost remains a constant. Therefore, the whole external excitation term can be regarded as a slow-varying parameter; thus, two generalized autonomous subsystems can be derived as follows:(21)Forω^>0,S1:di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L−c,Forω^<0,S1:di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L+c.The nonsmooth boundary can be denoted byΣΣ:i^d,i^q,ω^|ω^=0, which divides the state phase space into two regions represented by D± corresponding to ω^>0 and ω^<0, respectively. When the trajectory is located in D+, it is governed by subsystem S1. Once the trajectory reaches D− traversing the nonsmooth boundary, the governing equations will convert to S2 immediately.Compared with smooth system, Filippov system (20) may produce more abundant and complex dynamical behaviors. In order to further explore the mechanism of oscillations in nonsmooth PMSM system, the bifurcations of subsystems are ought to be investigated.The equilibrium points of the two subsystems can be uniformly expressed as(22)γω^021+ω^02,γω^01+ω^02,ω^0,where ω^0 satisfies the condition(23)T^L=σω^0γ−1−ω^021+ω^02+csignω^0.The properties of equilibrium points are determined by the characteristic equation of nonsmooth system, which is calculated as(24)λ3+σ+2λ2+1+σγω^021+ω^02−σγ+ω^02+2σλ+σω^02+2γω^021+ω^02−γ+1=0.Two typical bifurcations can also be observed in nonsmooth cases, which can be expressed as(25)FB:σω^02+2γω^021+ω^02−γ+1=0,σ+21+σγω^021+ω^02−σγ+ω^02+2σ−σω^02+2γω^021+ω^02−γ+1>0,T^L=σω^0γ−1−ω^021+ω^02+csignω^0,(26)HB:σω^02+2γω^021+ω^02−γ+1>0,σ+21+σγω^021+ω^02−σγ+ω^02+2σ−σω^02+2γω^021+ω^02−γ+1=0,T^L=σω^0γ−1−ω^021+ω^02+csignω^0.Now we turn to the oscillations of the PSMS systems with external periodic excitation. It should be noted that the bifurcation conditions are related to the slow-varying excitationT^L. Therefore, it is necessary to introduce the conception of transformed phase portrait (TPP) to show the influence of the bifurcations on the dynamics.


#### 2.5.2. TPP for the Trajectory

The traditional phase portraitΠ:i^dt^,i^qt^,ω^t^,t^∈R reflects the relationship between state variables in the phase space or the projection plane. If the slow-varying parameter is regarded as a generalized state variable, the generalized phase trajectory can be defined as Π:i^dt^,i^qt^,ω^t^,T^Lt^,t^∈R, where T^t^L=AsinΩt^. The generalized phase trajectory located in the generalized phase space i^d,i^q,ω^,T^L or its projected phase plane can be named as transformed phase portrait, which is usually employed to reveal the relationship between state variables and slow-varying parameters.


### 2.6. Dynamical Behavior Simulation of the Nonsmooth PMSM System

#### 2.6.1. Case A

It should be noted that the locations of equilibrium points and the bifurcation conditions can be both affected by the variation of the parameter. The equilibrium branches and corresponding bifurcations are plotted in Figure8(a), while the superposition of equilibrium curves and transformed phase portrait is shown in Figure 8(b) with the fixed parameters σ=4,γ=5,A=12,Ω=0.01,andc=0.Figure 8
(a) The equilibrium branches and (b) transformed phase portrait withσ=4,γ=5,A=12,Ω=0.01,andc=0.
(a)(b)It can be seen from Figure8(a) that there are three kinds of equilibrium curves E1±,E0,E2±, where E1± represents stable focal branches, E2± represents unstable saddle branches, and E0 describes unstable focal branch. Based on conditions (25) and (26), the coordinates of the Fold bifurcation points can be calculated as FB±±6.6044,±0.7288, while the sub-Hopf bifurcation points are sub-HB±±5.0497,±1.1975.As shown in Figure8(b), the trajectory moves strictly along the stable branch E1+ until sub-HB+ is reached. At the Hopf bifurcation point sub-HB+, the system is supposed to oscillate due to the instability of the equilibrium point. However, because of the bifurcation lag phenomenon [18], the trajectory keeps moving along the saddle branch E2+ instead of oscillating until arriving FB+. Then, due to the Fold bifurcation, the attraction of sole stable focal branch results in the transition and convergence of trajectory to E1− with a sharp oscillation. The above process describes the dynamical behaviors of the system in half a complete period. The other half is unnecessary to be elaborated because of the symmetry.Generally speaking, the motion state of the system along the stable equilibrium curve is called quiescent state (QS). Hopf bifurcation may result in the instability of equilibrium points and the generation of limit cycle, making the system oscillate locally, which can be named as spiking state (SP). Besides, Fold bifurcation can lead to the transition between equilibrium points, thus connecting different quiescent states and spiking states. The combination of quiescent state and spiking state is named as bursting oscillation. The bursting oscillation of the system can be also verified by the phase portrait and time history shown in Figure9.Figure 9
(a) The phase portrait and (b) time history withσ=4,γ=5,A=12,Ω=0.01,andc=0.
(a)(b)


#### 2.6.2. Case B

Another set of fixed parametersσ=4,γ=10,A=12,Ω=0.01,c=0 is taken to observe the change of dynamical behaviors. The superposition of equilibrium curves and TPP is shown in Figure 10, from which we find that the qualitative property of this periodic solution is different from the previous case because Fold bifurcation points do not participate in the oscillation structure at all.Figure 10
The superposition of equilibrium curves and TPP withσ=4,γ=10,A=12,Ω=0.01,andc=0.To be specific, the lag effect still occurs after the Hopf bifurcation, but the displacement and time length of bifurcation lag are not enough for the trajectory to reach the Fold bifurcation point due to the limited value range of external excitation [−12, +12] and the longer distance betweensub-HB−−3.93,−2.48 and FB−−16.34,−0.84. In other words, the spiking state SP1− (oscillation resulted from Hopf bifurcation) directly connects to the other spiking state SP2+ (convergence to stable focal branch E1+) without any type of quiescent state, which is definitely separate from the previous bursting mode. This interesting phenomenon can be also observed from the phase portrait and time history shown in Figure 11.Figure 11
(a) The phase portrait and (b) time history withσ=4,γ=10,A=12,Ω=0.01,andc=0.
(a)(b)


#### 2.6.3. Case C

The above two cases have described different modes of bursting oscillations in smooth nonautonomous PMSM system. However, the introduction of friction factor (c≠0) is going to destroy the smoothness of the system and make it become Filippov-type, resulting in more complex mixed oscillation mode.Taking a set of fixed parameters asσ=4,γ=5,A=12,Ω=0.01,c=12, the phase portrait and time history are plotted in Figure 12, which are totally different from smooth cases. It can be found that the structure of attractors keeps the symmetry of original system. The trajectory visits the two regions D± in turn, which are divided by the nonsmooth boundary ω^=0. The bursting oscillation combining the relatively large amplitude oscillation corresponding to the spiking state and the small amplitude oscillation related to the quiescent state can be still observed in this case. Furthermore, the special sliding motion on the nonsmooth boundary can be also found, which functions as the connection of two segments of the trajectory located in different regions D±. In order to reveal the mechanism of the sliding-bursting oscillation, the equilibrium curves and the transformed phase portrait of the two subsystems on T^L,ω^ plane are presented in Figure 13.Figure 12
(a) The phase portrait and (b) time history withσ=4,γ=5,A=12,Ω=0.01,andc=12.
(a)(b)Figure 13
(a) The equilibrium branches and (b) transformed phase portrait withσ=4,γ=5,A=12,Ω=0.01,andc=12.
(a)(b)


#### 2.6.4. Mechanism of Bursting Oscillations

Assume that the trajectory starts from the initial pointA1 in Figure 13, at which the slow-varying parameter T^L takes its minimum value T^L=−12. At this moment, due to the control from subsystem (21) and the attraction from the stable equilibrium point, the trajectory moves almost strictly along the stable focal branch E1+. Then, instead of the Hopf bifurcation, the trajectory still moves along the unstable saddle branch appearing in quiescent state until it reaches the Fold point FB1+−5.3956,0.7288 resulting in the transition to the nonsmooth boundary at A2. It should be noted that although there exists the Hopf point sub−HB1+−6.9503,1.1975, the slow passage effect [19] still leads to the lag of bifurcation. It can be found that when the trajectory reaches the nonsmooth boundary, it moves almost strictly along the boundary appearing in a special kind of quiescent state named sliding motion. Finally, the trajectory jumps off the nonsmooth boundary at A3 and converges to the stable equilibrium branch E1− indicating that half a period of bursting oscillation is completed. Due to the symmetry of the system, the other half period is unnecessary to be elaborated.


#### 2.6.5. Mechanism of Sliding Motion

The mechanism of sliding bifurcation can be explained from two aspects. On the one hand, the attractor structure on each side of the sliding trajectory is special. To be specific, when the trajectory locates in regionD+ near the nonsmooth boundary, it is governed by subsystem (21) and attracted by the sole stable focal branch E4+, which is inaccessible but functional, thus tending to cross the boundary. However, once the trajctory locates in region D−, it will be controlled by subsystem (22) and attracted by the sole stable focal branch E4−, thus returning to regionE4+. It should be pointed out that the inaccessible equilibrium E4± denoted by pink-dotted line refers to the equilibrium state that does not exist in the corresponding phase space of the control subsystem. Although the inaccessible equilibrium exists in the phase space of another subsystem, it can still affect the phase trajectory of the current control subsystem. Generally, the alternative control (attraction) of two subsystems (branches) results in the sliding motion along the nonsmooth boundary. In short, the sliding phenomenon may occur when there exists a sole stable attractor on each side of the nonsmooth boundary.On the other hand, the vector structure on each side of sliding trajectory is also interesting. In order to better understand the mechanism of sliding bifurcation, the sliding processA2⟶A3 in Figure 13(b) is numerically analyzed in Table 1.Table 1
Simulation data near points of special behaviors.
T^Lω^dω^dt^A2 (sliding)−5.2203>0−1.0090 < 0−5.2192<022.9106>0A3 (jumping)11.9756>0−23.9536 < 011.9757<0&ap;0,<0It is noteworthy that, in sliding regionA2⟶A3, where ω^≈0, the sign of vector field dω^/dt^ on each side of nonsmooth boundary ω^=0 remains opposite and the difference is approximately equal to 2c=2×12.0=24.0. Once the trajectory reaches the jumping point A3, the direction of vector field on each side of nonsmooth boundary is going to change from the opposite to the same. This slow-varying process can be also described by the sketch of vector field structure shown in Figure 14. It can be concluded that the vector field direction on each side of nonsmooth boundary is opposite in the sliding region, while the direction remains the same in the traversing region. The analytic expressions of the sliding region and the traversing region which were reported in our previous work [13] are omitted in this paper.Figure 14
The sketch of vector field structure in the sliding region.


#### 2.6.6. Time Length of the Sliding Motion

It can be found that the trajectory jumps quickly to the nonsmooth boundary and then moves along with it from the Fold bifurcation pointFB1+. Therefore, the time length of the sliding motion can be approximated by the time needed for slow-varying parameter T^L changing from FB1+ to A3. For FB1+, T^L=AsinΩt^1=−5.3956, while for A3, T^L=AsinΩt^2=12. The theoretical value of time length in Figure 15 can be calculated as Λ=t^1−t^2=203.7152, which coincides with the numerical result Λ1≈197.5186 from the time history (Figure 12(b)).Figure 15
The time history of slow-varying parameterT^L.


#### 2.6.7. Calculation of Natural and Spiking Frequency

The natural frequencyΩ0 of the system can be estimated by the imaginary part of equilibrium point eigenvalue with the amplitude of external excitation A=0. The equilibrium points of two subsystems in D± can be approximated as E±4.7248,±1.1403,±4.1436. The corresponding eigenvalues of Jacobi matrix are calculated as λ1=−4.7106,λ±=−0.6447±4.3104I implying that the approximate value of natural frequency Ω0=4.3104≫0.01=Ω.The spiking oscillations in this case are caused by the rapid convergence of stable focusEf±4,±2,±2 corresponding to the maximum amplitude of external excitation A=±12. The related eigenvalues can be computed as λ1=−5.4280,λ±=−0.2860±2.4111I indicating that the theoretical value of spiking frequency is ΩT=2.4111, which agrees well with the numerical result ΩSP=2π/TSP≈2.6059≈2.4111 obtained from Figure 16.Figure 16
The time history of spiking oscillation.


#### 2.6.8. Discussion on Physical Meaning

The spiking state in bursting corresponds to the high frequency oscillation in engineering practice which is undoubtedly harmful to the stability and life span of the equipment, but it can not be concluded that quiescent state is always beneficial. For instance, in engineering practice, the quiescent state caused by sliding bifurcation (friction factor) implies the stagnation or jamming of machines, which can result in the fatigue and abrasion of the internals. Therefore, it is urgent to study the mechanism of bursting oscillations with sliding bifurcation and propose a reasonable control scheme to eliminate these harmful vibrations.


## 2.1. Electromagnetic Part

The voltage equations for the two-phase rotating coordinate model of thed-axis and q-axis can be expressed as(2)ud=Rid+ddtψd−ωψq,uq=Riq+ddtψq+ωψd.where ψd and ψq represent the d-axis and q-axis components of the stator magnetic chain, respectively. Besides, the equations of the stator magnetic chain can be expressed as(3)ψd=Ldid+ψf,ψq=Lqiq.Combining the voltage and the stator magnetic chain equations, it can be concluded that(4)diddt=ud−R1id+ωLqiqLd,diqdt=ud−R1iq−ωLdid+ωψfLq.


## 2.2. Mechanical Motion Part

The mechanical motion equation of the motor can be obtained by Newton’s second law.(5)Jα=Jdωdt,=Te−TL−βω+fω,where the electromagnetic torque Te can be expressed as Te=npiqψf+Ld−Lqid based on the electromechanical energy conversion principle (the electric angle is equal to the product of the mechanical angle and the polar logarithm). Combining the electromagnetic part (4) and mechanical motion part (5), the control equations of three-phase PMSM system can be formed as system (1).In order to simplify the governing equations, the time scale transformation and linear affine transformation are introduced as(6)t=τ⋅t^,x=λx^,where τ=Lq/R1,x=id,iq,ωT,x^=i^d,i^q,ω^T, and(7)λ=LqβLdnpτψf000βnpτψf0001τ.Through the assumption of uniform air gapLd=Lq and dimensionless process,(8)γ=npψf2R1β,σ=LqβR1J,u^q=npψfLquqR12β,u^d=npψfLqudR12β,ε=Lqβ2Ld−LqLdJnpψf2,T^L=Lq2TLR12J,f^ω^=Lq2fω^R12J, and the governing equations can be written as(9)di^ddt^=−LqLdi^d+ω^iq+u^d,di^qdt^=−i^q−ω^i^d+γω^+u^q,dω^dt^=σi^q−ω^+εi^di^q−T^L+f^ω^,⇒Lq=Lddi^ddt^=−i^d+ω^iq+u^d,di^qdt^=−i^q−ω^i^d+γω^+u^q,dω^dt^=σi^q−ω^−T^L+f^ω^,where σ>0,γ>0.


## 2.3. Stability and Bifurcations of the Autonomous System

First of all, it is necessary to briefly investigate the inherent properties of PMSM system, so it can be assumed thatu^d=u^q=0,T^L=f^ω^=0 to further simplify the control equations as(10)di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^.The equilibrium point of the system satisfies the following conditions:(11)−i^d0+ω^0i^q0=0,−i^q0−ω^0i^d0+γω^0=0,σi^q0−ω^0=0,⇒i^d0=ω^02,i^q0=ω^0,ω^0ω^02−γ−1γ−1=0.which indicates that there are three equilibrium points in the system with γ>1:(12)E00,0,0,E+γ−1,γ−1,γ−1,E−γ−1,−γ−1,−γ−1,while E0 is the sole equilibrium with γ≤1.Routh-Hurwitz criterion can be employed to analyze the stability of equilibrium point. The Jacobi matrix of system (10) can be expressed as(13)J=−1ω^0i^q0−ω^0−1−i^d0+γ0σ−σ.It can be seen that all equilibrium points of the system can be expressed in a unified formω^02,ω^0,ω^0. Therefore, the characteristic equation at the equilibrium point can be formed as(14)λ3+σ+2σ2+1+σω^02−σγ+ω^02+2σλ+3σω^02−σγ+σ=0.Based on the inference of Routh–Hurwitz stability criterion, for a cubic polynomial(15)P3z=z3+a1z2+a2z3+a3,its Routh-Hurwitz conditions can be formed as(16)a1>0,a11a3a2>0,a3a1a2−a3>0,⇕a1>0,a3>0,a1a2>a3.Therefore, the stability conditions of the equilibrium point can be expressed as(17)σ>−2,σω^02−γ+1>0,σ+21+σω^0−σγ+ω^02+2σ−σω^02−γ+1>0.When conditiona3>0 degenerates to a3=0, the characteristic equation has a zero root corresponding to the Fold bifurcation conditions.(18)FB:σ>−2,γ=1,ω^0=0,σ+21+σ>0,⇒σ>0γ=1,ω^0=0.When conditiona1a2>a3 degenerates to a1a2=a3, a pair of pure imaginary roots appear in the solution of characteristic polynomial (14) indicating the occurrence of Hopf bifurcation. Therefore, the conditions of Hopf bifurcation are as follows:(19)HB:σ>−2,σω^02−γ+1>0,σ+21+σω^0−σγ+ω^02+2σ−σω^02+γ−1=0,⟹σ>0,γ>1,ω^02=γ−1σ2−σγ+2γ+4σ=0.Based on above analysis, the two-parameter bifurcation diagram onγ,σ plane is shown in Figure 2 in which FB represents the Fold bifurcation curve, while HB describes the subcritical Hopf bifurcation curve. It can be seen that the parameter plane is divided into three parts by the two bifurcation curves motioned. In region 1, the system has a sole stable equilibrium point. With the variation of parameter combinations, it traverses the Fold curve and enters region 2 leading to the appearance of two new stable equilibrium points. Then, the equilibrium points continue to cross the Hopf curve and then enter region 3, resulting in the existence of three unstable equilibrium points.Figure 2
Two-parameter bifurcation diagram onσ,γ plane.Taking the fixed parameterσ=4.0 and variable parameter γ, the bifurcation diagram of equilibrium points is shown in Figure 3 in which the solid line represents the stable equilibrium point, while the dotted line corresponds to the unstable ones. The specific parameters of Fold and Hopf bifurcations are calculated as γ=1.0 and γ=16.0, respectively. Two unstable limit cycles generated by the subcritical Hopf bifurcation points HB± are denoted by LPC1,2. It can be seen that the property changes of equilibrium points are consistent with Figure 2.Figure 3
Bifurcation diagram of autonomous PMSM system withσ=4.0.


## 2.4. Dynamical Behavior Simulation of the Autonomous System

### 2.4.1. Case A

Taking a set of fixed parameters in region 1 (Figure2) as σ=4.0,γ=0.1, the system has the sole equilibrium point E0′0,0,0 as shown in Figure 4(a), while, in Figure 4(b), there are two stable focal equilibrium points E±4,±2,±2 and an unstable saddle equilibrium point E0′0,0,0 with another set of fixed parameters σ=4.0,γ=5.0 in region 2.Figure 4
Phase portraits onω^,i^d plane with (a) σ=4.0,γ=0.1 and (b) σ=4.0,γ=5.0.
(a)(b)


### 2.4.2. Case B

Taking another set of parameters in region 2 asσ=4.0,γ=15.0, the number and stability of equilibrium points remain unchanged, but the convergence modes of trajectory to the stable focus are obviously different due to the variation of initial values. In Figure 5(a), the envelope of the trajectory is similar to the chaotic attractor which can be also confirmed by the quasi-chaotic process in time history (Figure 6(a)). Then, the trajectory suddenly approaches the stable focus by the means of focal-type convergence after the quasi-chaotic process. However, with another set of initial values, the trajectory directly converges to the stable focus as shown in Figures 5(b) and 6(b).Figure 5
Phase portraits with initial valuesi^d,i^q,ω^ of (a) 20,3.74,3.74 and (b) 15,3.74,3.74.
(a)(b)Figure 6
Time histories with initial valuesi^d,i^q,ω^ of (a) 20,3.74,3.74 and (b) 15,3.74,3.74.
(a)(b)


### 2.4.3. Case C

In the two cases mentioned above, the system can operate stably. Taking the last set of parameters in region 3 asσ=4.0,γ=17.0, it can be found from the phase portrait (Figure 7(a)) and time history (Figure 7(c)) that the chaotic attractor appears after the Hopf bifurcation (σ=4.0,γ=16.0). In addition, the dynamical evolution of the autonomous PMSM system (10) and the attraction domain transition of stable attractors can be also proved by Poincaré map (Figure 7(b)) and the maximum Lyapunov exponent diagram (Figure 7(d)). To be specific, maximum Lyapunov exponent λmax=0.3726>0 at the Hopf bifurcation point, which indicates that the two adjacent phase trajectories are going to separate gradually, making the system enter the chaotic state finally.Figure 7
Chaos of the autonomous system: (a) phase portrait; (b) Poincaré map; (c) time history; (d) the maximum Lyapunov exponent diagram.
(a)(b)(c)(d)


### 2.4.4. Discussion on Physical Meaning

The equilibrium point in the mathematical model corresponds to the stable working state of the machine in engineering practice. The chaotic solution or quasi-chaotic convergence process of the system represents the irregular oscillation of the equipment, which usually leads to the aging and damage of motor and thus is ought to be avoided. In this section, the self-excited operation of the PMSM system without external load is mainly discussed. The results show that the risk of irregular oscillation can be effectively reduced by proper parameter combination and initial value selection. However, due to the existence of multiple excitations (external excitation, parametric excitation, etc.) in engineering practice, the PMSM system may produce multiple-scale effect inducing the bursting oscillation which is more complex. In addition to the oscillation, the existence of friction factors may further aggravate the equipment damage and even change the inherent properties of the system. These problems will be discussed in the next section.


## 2.4.1. Case A

Taking a set of fixed parameters in region 1 (Figure2) as σ=4.0,γ=0.1, the system has the sole equilibrium point E0′0,0,0 as shown in Figure 4(a), while, in Figure 4(b), there are two stable focal equilibrium points E±4,±2,±2 and an unstable saddle equilibrium point E0′0,0,0 with another set of fixed parameters σ=4.0,γ=5.0 in region 2.Figure 4
Phase portraits onω^,i^d plane with (a) σ=4.0,γ=0.1 and (b) σ=4.0,γ=5.0.
(a)(b)


## 2.4.2. Case B

Taking another set of parameters in region 2 asσ=4.0,γ=15.0, the number and stability of equilibrium points remain unchanged, but the convergence modes of trajectory to the stable focus are obviously different due to the variation of initial values. In Figure 5(a), the envelope of the trajectory is similar to the chaotic attractor which can be also confirmed by the quasi-chaotic process in time history (Figure 6(a)). Then, the trajectory suddenly approaches the stable focus by the means of focal-type convergence after the quasi-chaotic process. However, with another set of initial values, the trajectory directly converges to the stable focus as shown in Figures 5(b) and 6(b).Figure 5
Phase portraits with initial valuesi^d,i^q,ω^ of (a) 20,3.74,3.74 and (b) 15,3.74,3.74.
(a)(b)Figure 6
Time histories with initial valuesi^d,i^q,ω^ of (a) 20,3.74,3.74 and (b) 15,3.74,3.74.
(a)(b)


## 2.4.3. Case C

In the two cases mentioned above, the system can operate stably. Taking the last set of parameters in region 3 asσ=4.0,γ=17.0, it can be found from the phase portrait (Figure 7(a)) and time history (Figure 7(c)) that the chaotic attractor appears after the Hopf bifurcation (σ=4.0,γ=16.0). In addition, the dynamical evolution of the autonomous PMSM system (10) and the attraction domain transition of stable attractors can be also proved by Poincaré map (Figure 7(b)) and the maximum Lyapunov exponent diagram (Figure 7(d)). To be specific, maximum Lyapunov exponent λmax=0.3726>0 at the Hopf bifurcation point, which indicates that the two adjacent phase trajectories are going to separate gradually, making the system enter the chaotic state finally.Figure 7
Chaos of the autonomous system: (a) phase portrait; (b) Poincaré map; (c) time history; (d) the maximum Lyapunov exponent diagram.
(a)(b)(c)(d)


## 2.4.4. Discussion on Physical Meaning

The equilibrium point in the mathematical model corresponds to the stable working state of the machine in engineering practice. The chaotic solution or quasi-chaotic convergence process of the system represents the irregular oscillation of the equipment, which usually leads to the aging and damage of motor and thus is ought to be avoided. In this section, the self-excited operation of the PMSM system without external load is mainly discussed. The results show that the risk of irregular oscillation can be effectively reduced by proper parameter combination and initial value selection. However, due to the existence of multiple excitations (external excitation, parametric excitation, etc.) in engineering practice, the PMSM system may produce multiple-scale effect inducing the bursting oscillation which is more complex. In addition to the oscillation, the existence of friction factors may further aggravate the equipment damage and even change the inherent properties of the system. These problems will be discussed in the next section.


## 2.5. Stability and Bifurcations of the Nonsmooth PMSM System

### 2.5.1. Bifurcation Analysis of Subsystems

In order to explore the typical dynamical behaviors of nonsmooth PMSM system, the symbolic functionf^ω^=c⋅signω^=−c,ω^>0+c,ω^<0,c≥0 is employed to describe the dry friction, while the excitation function T^L=AsinΩt^ is used to represent the external torque. It should be noted that the definition of symbolic function value with ω^=0 is omitted for convenience because it depends on the equilibrium equation (similar to static friction). Therefore, the control equations of Filippov-type PMSM can be written as(20)di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L+f^ω^.It should be pointed out that the effect of two-time scales [18] may appear when the excitation frequency Ω is far less than the natural frequency Ω0 which can be estimated by imaginary part of equilibrium point eigenvalue.There will be magnitude differences between external excitation frequency and natural frequency if0<Ω≪1. In any natural frequency period TΩ0=2π/Ω0,t^∈t^0,t^0+TΩ0, the variation range of external excitation is WA,WB, where WA=AsinΩt^0 and WB=AsinΩt0+2πΩ/ω. Based on the assumption Ω/ω≪1⇒WA≈WB, it can be concluded that, in any natural frequency period, the external excitation almost remains a constant. Therefore, the whole external excitation term can be regarded as a slow-varying parameter; thus, two generalized autonomous subsystems can be derived as follows:(21)Forω^>0,S1:di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L−c,Forω^<0,S1:di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L+c.The nonsmooth boundary can be denoted byΣΣ:i^d,i^q,ω^|ω^=0, which divides the state phase space into two regions represented by D± corresponding to ω^>0 and ω^<0, respectively. When the trajectory is located in D+, it is governed by subsystem S1. Once the trajectory reaches D− traversing the nonsmooth boundary, the governing equations will convert to S2 immediately.Compared with smooth system, Filippov system (20) may produce more abundant and complex dynamical behaviors. In order to further explore the mechanism of oscillations in nonsmooth PMSM system, the bifurcations of subsystems are ought to be investigated.The equilibrium points of the two subsystems can be uniformly expressed as(22)γω^021+ω^02,γω^01+ω^02,ω^0,where ω^0 satisfies the condition(23)T^L=σω^0γ−1−ω^021+ω^02+csignω^0.The properties of equilibrium points are determined by the characteristic equation of nonsmooth system, which is calculated as(24)λ3+σ+2λ2+1+σγω^021+ω^02−σγ+ω^02+2σλ+σω^02+2γω^021+ω^02−γ+1=0.Two typical bifurcations can also be observed in nonsmooth cases, which can be expressed as(25)FB:σω^02+2γω^021+ω^02−γ+1=0,σ+21+σγω^021+ω^02−σγ+ω^02+2σ−σω^02+2γω^021+ω^02−γ+1>0,T^L=σω^0γ−1−ω^021+ω^02+csignω^0,(26)HB:σω^02+2γω^021+ω^02−γ+1>0,σ+21+σγω^021+ω^02−σγ+ω^02+2σ−σω^02+2γω^021+ω^02−γ+1=0,T^L=σω^0γ−1−ω^021+ω^02+csignω^0.Now we turn to the oscillations of the PSMS systems with external periodic excitation. It should be noted that the bifurcation conditions are related to the slow-varying excitationT^L. Therefore, it is necessary to introduce the conception of transformed phase portrait (TPP) to show the influence of the bifurcations on the dynamics.


### 2.5.2. TPP for the Trajectory

The traditional phase portraitΠ:i^dt^,i^qt^,ω^t^,t^∈R reflects the relationship between state variables in the phase space or the projection plane. If the slow-varying parameter is regarded as a generalized state variable, the generalized phase trajectory can be defined as Π:i^dt^,i^qt^,ω^t^,T^Lt^,t^∈R, where T^t^L=AsinΩt^. The generalized phase trajectory located in the generalized phase space i^d,i^q,ω^,T^L or its projected phase plane can be named as transformed phase portrait, which is usually employed to reveal the relationship between state variables and slow-varying parameters.


## 2.5.1. Bifurcation Analysis of Subsystems

In order to explore the typical dynamical behaviors of nonsmooth PMSM system, the symbolic functionf^ω^=c⋅signω^=−c,ω^>0+c,ω^<0,c≥0 is employed to describe the dry friction, while the excitation function T^L=AsinΩt^ is used to represent the external torque. It should be noted that the definition of symbolic function value with ω^=0 is omitted for convenience because it depends on the equilibrium equation (similar to static friction). Therefore, the control equations of Filippov-type PMSM can be written as(20)di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L+f^ω^.It should be pointed out that the effect of two-time scales [18] may appear when the excitation frequency Ω is far less than the natural frequency Ω0 which can be estimated by imaginary part of equilibrium point eigenvalue.There will be magnitude differences between external excitation frequency and natural frequency if0<Ω≪1. In any natural frequency period TΩ0=2π/Ω0,t^∈t^0,t^0+TΩ0, the variation range of external excitation is WA,WB, where WA=AsinΩt^0 and WB=AsinΩt0+2πΩ/ω. Based on the assumption Ω/ω≪1⇒WA≈WB, it can be concluded that, in any natural frequency period, the external excitation almost remains a constant. Therefore, the whole external excitation term can be regarded as a slow-varying parameter; thus, two generalized autonomous subsystems can be derived as follows:(21)Forω^>0,S1:di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L−c,Forω^<0,S1:di^ddt^=−i^d+ω^i^q,di^qdt^=−i^q−ω^i^d+γω^,dω^dt^=σi^q−ω^−T^L+c.The nonsmooth boundary can be denoted byΣΣ:i^d,i^q,ω^|ω^=0, which divides the state phase space into two regions represented by D± corresponding to ω^>0 and ω^<0, respectively. When the trajectory is located in D+, it is governed by subsystem S1. Once the trajectory reaches D− traversing the nonsmooth boundary, the governing equations will convert to S2 immediately.Compared with smooth system, Filippov system (20) may produce more abundant and complex dynamical behaviors. In order to further explore the mechanism of oscillations in nonsmooth PMSM system, the bifurcations of subsystems are ought to be investigated.The equilibrium points of the two subsystems can be uniformly expressed as(22)γω^021+ω^02,γω^01+ω^02,ω^0,where ω^0 satisfies the condition(23)T^L=σω^0γ−1−ω^021+ω^02+csignω^0.The properties of equilibrium points are determined by the characteristic equation of nonsmooth system, which is calculated as(24)λ3+σ+2λ2+1+σγω^021+ω^02−σγ+ω^02+2σλ+σω^02+2γω^021+ω^02−γ+1=0.Two typical bifurcations can also be observed in nonsmooth cases, which can be expressed as(25)FB:σω^02+2γω^021+ω^02−γ+1=0,σ+21+σγω^021+ω^02−σγ+ω^02+2σ−σω^02+2γω^021+ω^02−γ+1>0,T^L=σω^0γ−1−ω^021+ω^02+csignω^0,(26)HB:σω^02+2γω^021+ω^02−γ+1>0,σ+21+σγω^021+ω^02−σγ+ω^02+2σ−σω^02+2γω^021+ω^02−γ+1=0,T^L=σω^0γ−1−ω^021+ω^02+csignω^0.Now we turn to the oscillations of the PSMS systems with external periodic excitation. It should be noted that the bifurcation conditions are related to the slow-varying excitationT^L. Therefore, it is necessary to introduce the conception of transformed phase portrait (TPP) to show the influence of the bifurcations on the dynamics.


## 2.5.2. TPP for the Trajectory

The traditional phase portraitΠ:i^dt^,i^qt^,ω^t^,t^∈R reflects the relationship between state variables in the phase space or the projection plane. If the slow-varying parameter is regarded as a generalized state variable, the generalized phase trajectory can be defined as Π:i^dt^,i^qt^,ω^t^,T^Lt^,t^∈R, where T^t^L=AsinΩt^. The generalized phase trajectory located in the generalized phase space i^d,i^q,ω^,T^L or its projected phase plane can be named as transformed phase portrait, which is usually employed to reveal the relationship between state variables and slow-varying parameters.


## 2.6. Dynamical Behavior Simulation of the Nonsmooth PMSM System

### 2.6.1. Case A

It should be noted that the locations of equilibrium points and the bifurcation conditions can be both affected by the variation of the parameter. The equilibrium branches and corresponding bifurcations are plotted in Figure8(a), while the superposition of equilibrium curves and transformed phase portrait is shown in Figure 8(b) with the fixed parameters σ=4,γ=5,A=12,Ω=0.01,andc=0.Figure 8
(a) The equilibrium branches and (b) transformed phase portrait withσ=4,γ=5,A=12,Ω=0.01,andc=0.
(a)(b)It can be seen from Figure8(a) that there are three kinds of equilibrium curves E1±,E0,E2±, where E1± represents stable focal branches, E2± represents unstable saddle branches, and E0 describes unstable focal branch. Based on conditions (25) and (26), the coordinates of the Fold bifurcation points can be calculated as FB±±6.6044,±0.7288, while the sub-Hopf bifurcation points are sub-HB±±5.0497,±1.1975.As shown in Figure8(b), the trajectory moves strictly along the stable branch E1+ until sub-HB+ is reached. At the Hopf bifurcation point sub-HB+, the system is supposed to oscillate due to the instability of the equilibrium point. However, because of the bifurcation lag phenomenon [18], the trajectory keeps moving along the saddle branch E2+ instead of oscillating until arriving FB+. Then, due to the Fold bifurcation, the attraction of sole stable focal branch results in the transition and convergence of trajectory to E1− with a sharp oscillation. The above process describes the dynamical behaviors of the system in half a complete period. The other half is unnecessary to be elaborated because of the symmetry.Generally speaking, the motion state of the system along the stable equilibrium curve is called quiescent state (QS). Hopf bifurcation may result in the instability of equilibrium points and the generation of limit cycle, making the system oscillate locally, which can be named as spiking state (SP). Besides, Fold bifurcation can lead to the transition between equilibrium points, thus connecting different quiescent states and spiking states. The combination of quiescent state and spiking state is named as bursting oscillation. The bursting oscillation of the system can be also verified by the phase portrait and time history shown in Figure9.Figure 9
(a) The phase portrait and (b) time history withσ=4,γ=5,A=12,Ω=0.01,andc=0.
(a)(b)


### 2.6.2. Case B

Another set of fixed parametersσ=4,γ=10,A=12,Ω=0.01,c=0 is taken to observe the change of dynamical behaviors. The superposition of equilibrium curves and TPP is shown in Figure 10, from which we find that the qualitative property of this periodic solution is different from the previous case because Fold bifurcation points do not participate in the oscillation structure at all.Figure 10
The superposition of equilibrium curves and TPP withσ=4,γ=10,A=12,Ω=0.01,andc=0.To be specific, the lag effect still occurs after the Hopf bifurcation, but the displacement and time length of bifurcation lag are not enough for the trajectory to reach the Fold bifurcation point due to the limited value range of external excitation [−12, +12] and the longer distance betweensub-HB−−3.93,−2.48 and FB−−16.34,−0.84. In other words, the spiking state SP1− (oscillation resulted from Hopf bifurcation) directly connects to the other spiking state SP2+ (convergence to stable focal branch E1+) without any type of quiescent state, which is definitely separate from the previous bursting mode. This interesting phenomenon can be also observed from the phase portrait and time history shown in Figure 11.Figure 11
(a) The phase portrait and (b) time history withσ=4,γ=10,A=12,Ω=0.01,andc=0.
(a)(b)


### 2.6.3. Case C

The above two cases have described different modes of bursting oscillations in smooth nonautonomous PMSM system. However, the introduction of friction factor (c≠0) is going to destroy the smoothness of the system and make it become Filippov-type, resulting in more complex mixed oscillation mode.Taking a set of fixed parameters asσ=4,γ=5,A=12,Ω=0.01,c=12, the phase portrait and time history are plotted in Figure 12, which are totally different from smooth cases. It can be found that the structure of attractors keeps the symmetry of original system. The trajectory visits the two regions D± in turn, which are divided by the nonsmooth boundary ω^=0. The bursting oscillation combining the relatively large amplitude oscillation corresponding to the spiking state and the small amplitude oscillation related to the quiescent state can be still observed in this case. Furthermore, the special sliding motion on the nonsmooth boundary can be also found, which functions as the connection of two segments of the trajectory located in different regions D±. In order to reveal the mechanism of the sliding-bursting oscillation, the equilibrium curves and the transformed phase portrait of the two subsystems on T^L,ω^ plane are presented in Figure 13.Figure 12
(a) The phase portrait and (b) time history withσ=4,γ=5,A=12,Ω=0.01,andc=12.
(a)(b)Figure 13
(a) The equilibrium branches and (b) transformed phase portrait withσ=4,γ=5,A=12,Ω=0.01,andc=12.
(a)(b)


### 2.6.4. Mechanism of Bursting Oscillations

Assume that the trajectory starts from the initial pointA1 in Figure 13, at which the slow-varying parameter T^L takes its minimum value T^L=−12. At this moment, due to the control from subsystem (21) and the attraction from the stable equilibrium point, the trajectory moves almost strictly along the stable focal branch E1+. Then, instead of the Hopf bifurcation, the trajectory still moves along the unstable saddle branch appearing in quiescent state until it reaches the Fold point FB1+−5.3956,0.7288 resulting in the transition to the nonsmooth boundary at A2. It should be noted that although there exists the Hopf point sub−HB1+−6.9503,1.1975, the slow passage effect [19] still leads to the lag of bifurcation. It can be found that when the trajectory reaches the nonsmooth boundary, it moves almost strictly along the boundary appearing in a special kind of quiescent state named sliding motion. Finally, the trajectory jumps off the nonsmooth boundary at A3 and converges to the stable equilibrium branch E1− indicating that half a period of bursting oscillation is completed. Due to the symmetry of the system, the other half period is unnecessary to be elaborated.


### 2.6.5. Mechanism of Sliding Motion

The mechanism of sliding bifurcation can be explained from two aspects. On the one hand, the attractor structure on each side of the sliding trajectory is special. To be specific, when the trajectory locates in regionD+ near the nonsmooth boundary, it is governed by subsystem (21) and attracted by the sole stable focal branch E4+, which is inaccessible but functional, thus tending to cross the boundary. However, once the trajctory locates in region D−, it will be controlled by subsystem (22) and attracted by the sole stable focal branch E4−, thus returning to regionE4+. It should be pointed out that the inaccessible equilibrium E4± denoted by pink-dotted line refers to the equilibrium state that does not exist in the corresponding phase space of the control subsystem. Although the inaccessible equilibrium exists in the phase space of another subsystem, it can still affect the phase trajectory of the current control subsystem. Generally, the alternative control (attraction) of two subsystems (branches) results in the sliding motion along the nonsmooth boundary. In short, the sliding phenomenon may occur when there exists a sole stable attractor on each side of the nonsmooth boundary.On the other hand, the vector structure on each side of sliding trajectory is also interesting. In order to better understand the mechanism of sliding bifurcation, the sliding processA2⟶A3 in Figure 13(b) is numerically analyzed in Table 1.Table 1
Simulation data near points of special behaviors.
T^Lω^dω^dt^A2 (sliding)−5.2203>0−1.0090 < 0−5.2192<022.9106>0A3 (jumping)11.9756>0−23.9536 < 011.9757<0&ap;0,<0It is noteworthy that, in sliding regionA2⟶A3, where ω^≈0, the sign of vector field dω^/dt^ on each side of nonsmooth boundary ω^=0 remains opposite and the difference is approximately equal to 2c=2×12.0=24.0. Once the trajectory reaches the jumping point A3, the direction of vector field on each side of nonsmooth boundary is going to change from the opposite to the same. This slow-varying process can be also described by the sketch of vector field structure shown in Figure 14. It can be concluded that the vector field direction on each side of nonsmooth boundary is opposite in the sliding region, while the direction remains the same in the traversing region. The analytic expressions of the sliding region and the traversing region which were reported in our previous work [13] are omitted in this paper.Figure 14
The sketch of vector field structure in the sliding region.


### 2.6.6. Time Length of the Sliding Motion

It can be found that the trajectory jumps quickly to the nonsmooth boundary and then moves along with it from the Fold bifurcation pointFB1+. Therefore, the time length of the sliding motion can be approximated by the time needed for slow-varying parameter T^L changing from FB1+ to A3. For FB1+, T^L=AsinΩt^1=−5.3956, while for A3, T^L=AsinΩt^2=12. The theoretical value of time length in Figure 15 can be calculated as Λ=t^1−t^2=203.7152, which coincides with the numerical result Λ1≈197.5186 from the time history (Figure 12(b)).Figure 15
The time history of slow-varying parameterT^L.


### 2.6.7. Calculation of Natural and Spiking Frequency

The natural frequencyΩ0 of the system can be estimated by the imaginary part of equilibrium point eigenvalue with the amplitude of external excitation A=0. The equilibrium points of two subsystems in D± can be approximated as E±4.7248,±1.1403,±4.1436. The corresponding eigenvalues of Jacobi matrix are calculated as λ1=−4.7106,λ±=−0.6447±4.3104I implying that the approximate value of natural frequency Ω0=4.3104≫0.01=Ω.The spiking oscillations in this case are caused by the rapid convergence of stable focusEf±4,±2,±2 corresponding to the maximum amplitude of external excitation A=±12. The related eigenvalues can be computed as λ1=−5.4280,λ±=−0.2860±2.4111I indicating that the theoretical value of spiking frequency is ΩT=2.4111, which agrees well with the numerical result ΩSP=2π/TSP≈2.6059≈2.4111 obtained from Figure 16.Figure 16
The time history of spiking oscillation.


### 2.6.8. Discussion on Physical Meaning

The spiking state in bursting corresponds to the high frequency oscillation in engineering practice which is undoubtedly harmful to the stability and life span of the equipment, but it can not be concluded that quiescent state is always beneficial. For instance, in engineering practice, the quiescent state caused by sliding bifurcation (friction factor) implies the stagnation or jamming of machines, which can result in the fatigue and abrasion of the internals. Therefore, it is urgent to study the mechanism of bursting oscillations with sliding bifurcation and propose a reasonable control scheme to eliminate these harmful vibrations.


## 2.6.1. Case A

It should be noted that the locations of equilibrium points and the bifurcation conditions can be both affected by the variation of the parameter. The equilibrium branches and corresponding bifurcations are plotted in Figure8(a), while the superposition of equilibrium curves and transformed phase portrait is shown in Figure 8(b) with the fixed parameters σ=4,γ=5,A=12,Ω=0.01,andc=0.Figure 8
(a) The equilibrium branches and (b) transformed phase portrait withσ=4,γ=5,A=12,Ω=0.01,andc=0.
(a)(b)It can be seen from Figure8(a) that there are three kinds of equilibrium curves E1±,E0,E2±, where E1± represents stable focal branches, E2± represents unstable saddle branches, and E0 describes unstable focal branch. Based on conditions (25) and (26), the coordinates of the Fold bifurcation points can be calculated as FB±±6.6044,±0.7288, while the sub-Hopf bifurcation points are sub-HB±±5.0497,±1.1975.As shown in Figure8(b), the trajectory moves strictly along the stable branch E1+ until sub-HB+ is reached. At the Hopf bifurcation point sub-HB+, the system is supposed to oscillate due to the instability of the equilibrium point. However, because of the bifurcation lag phenomenon [18], the trajectory keeps moving along the saddle branch E2+ instead of oscillating until arriving FB+. Then, due to the Fold bifurcation, the attraction of sole stable focal branch results in the transition and convergence of trajectory to E1− with a sharp oscillation. The above process describes the dynamical behaviors of the system in half a complete period. The other half is unnecessary to be elaborated because of the symmetry.Generally speaking, the motion state of the system along the stable equilibrium curve is called quiescent state (QS). Hopf bifurcation may result in the instability of equilibrium points and the generation of limit cycle, making the system oscillate locally, which can be named as spiking state (SP). Besides, Fold bifurcation can lead to the transition between equilibrium points, thus connecting different quiescent states and spiking states. The combination of quiescent state and spiking state is named as bursting oscillation. The bursting oscillation of the system can be also verified by the phase portrait and time history shown in Figure9.Figure 9
(a) The phase portrait and (b) time history withσ=4,γ=5,A=12,Ω=0.01,andc=0.
(a)(b)


## 2.6.2. Case B

Another set of fixed parametersσ=4,γ=10,A=12,Ω=0.01,c=0 is taken to observe the change of dynamical behaviors. The superposition of equilibrium curves and TPP is shown in Figure 10, from which we find that the qualitative property of this periodic solution is different from the previous case because Fold bifurcation points do not participate in the oscillation structure at all.Figure 10
The superposition of equilibrium curves and TPP withσ=4,γ=10,A=12,Ω=0.01,andc=0.To be specific, the lag effect still occurs after the Hopf bifurcation, but the displacement and time length of bifurcation lag are not enough for the trajectory to reach the Fold bifurcation point due to the limited value range of external excitation [−12, +12] and the longer distance betweensub-HB−−3.93,−2.48 and FB−−16.34,−0.84. In other words, the spiking state SP1− (oscillation resulted from Hopf bifurcation) directly connects to the other spiking state SP2+ (convergence to stable focal branch E1+) without any type of quiescent state, which is definitely separate from the previous bursting mode. This interesting phenomenon can be also observed from the phase portrait and time history shown in Figure 11.Figure 11
(a) The phase portrait and (b) time history withσ=4,γ=10,A=12,Ω=0.01,andc=0.
(a)(b)


## 2.6.3. Case C

The above two cases have described different modes of bursting oscillations in smooth nonautonomous PMSM system. However, the introduction of friction factor (c≠0) is going to destroy the smoothness of the system and make it become Filippov-type, resulting in more complex mixed oscillation mode.Taking a set of fixed parameters asσ=4,γ=5,A=12,Ω=0.01,c=12, the phase portrait and time history are plotted in Figure 12, which are totally different from smooth cases. It can be found that the structure of attractors keeps the symmetry of original system. The trajectory visits the two regions D± in turn, which are divided by the nonsmooth boundary ω^=0. The bursting oscillation combining the relatively large amplitude oscillation corresponding to the spiking state and the small amplitude oscillation related to the quiescent state can be still observed in this case. Furthermore, the special sliding motion on the nonsmooth boundary can be also found, which functions as the connection of two segments of the trajectory located in different regions D±. In order to reveal the mechanism of the sliding-bursting oscillation, the equilibrium curves and the transformed phase portrait of the two subsystems on T^L,ω^ plane are presented in Figure 13.Figure 12
(a) The phase portrait and (b) time history withσ=4,γ=5,A=12,Ω=0.01,andc=12.
(a)(b)Figure 13
(a) The equilibrium branches and (b) transformed phase portrait withσ=4,γ=5,A=12,Ω=0.01,andc=12.
(a)(b)


## 2.6.4. Mechanism of Bursting Oscillations

Assume that the trajectory starts from the initial pointA1 in Figure 13, at which the slow-varying parameter T^L takes its minimum value T^L=−12. At this moment, due to the control from subsystem (21) and the attraction from the stable equilibrium point, the trajectory moves almost strictly along the stable focal branch E1+. Then, instead of the Hopf bifurcation, the trajectory still moves along the unstable saddle branch appearing in quiescent state until it reaches the Fold point FB1+−5.3956,0.7288 resulting in the transition to the nonsmooth boundary at A2. It should be noted that although there exists the Hopf point sub−HB1+−6.9503,1.1975, the slow passage effect [19] still leads to the lag of bifurcation. It can be found that when the trajectory reaches the nonsmooth boundary, it moves almost strictly along the boundary appearing in a special kind of quiescent state named sliding motion. Finally, the trajectory jumps off the nonsmooth boundary at A3 and converges to the stable equilibrium branch E1− indicating that half a period of bursting oscillation is completed. Due to the symmetry of the system, the other half period is unnecessary to be elaborated.


## 2.6.5. Mechanism of Sliding Motion

The mechanism of sliding bifurcation can be explained from two aspects. On the one hand, the attractor structure on each side of the sliding trajectory is special. To be specific, when the trajectory locates in regionD+ near the nonsmooth boundary, it is governed by subsystem (21) and attracted by the sole stable focal branch E4+, which is inaccessible but functional, thus tending to cross the boundary. However, once the trajctory locates in region D−, it will be controlled by subsystem (22) and attracted by the sole stable focal branch E4−, thus returning to regionE4+. It should be pointed out that the inaccessible equilibrium E4± denoted by pink-dotted line refers to the equilibrium state that does not exist in the corresponding phase space of the control subsystem. Although the inaccessible equilibrium exists in the phase space of another subsystem, it can still affect the phase trajectory of the current control subsystem. Generally, the alternative control (attraction) of two subsystems (branches) results in the sliding motion along the nonsmooth boundary. In short, the sliding phenomenon may occur when there exists a sole stable attractor on each side of the nonsmooth boundary.On the other hand, the vector structure on each side of sliding trajectory is also interesting. In order to better understand the mechanism of sliding bifurcation, the sliding processA2⟶A3 in Figure 13(b) is numerically analyzed in Table 1.Table 1
Simulation data near points of special behaviors.
T^Lω^dω^dt^A2 (sliding)−5.2203>0−1.0090 < 0−5.2192<022.9106>0A3 (jumping)11.9756>0−23.9536 < 011.9757<0&ap;0,<0It is noteworthy that, in sliding regionA2⟶A3, where ω^≈0, the sign of vector field dω^/dt^ on each side of nonsmooth boundary ω^=0 remains opposite and the difference is approximately equal to 2c=2×12.0=24.0. Once the trajectory reaches the jumping point A3, the direction of vector field on each side of nonsmooth boundary is going to change from the opposite to the same. This slow-varying process can be also described by the sketch of vector field structure shown in Figure 14. It can be concluded that the vector field direction on each side of nonsmooth boundary is opposite in the sliding region, while the direction remains the same in the traversing region. The analytic expressions of the sliding region and the traversing region which were reported in our previous work [13] are omitted in this paper.Figure 14
The sketch of vector field structure in the sliding region.


## 2.6.6. Time Length of the Sliding Motion

It can be found that the trajectory jumps quickly to the nonsmooth boundary and then moves along with it from the Fold bifurcation pointFB1+. Therefore, the time length of the sliding motion can be approximated by the time needed for slow-varying parameter T^L changing from FB1+ to A3. For FB1+, T^L=AsinΩt^1=−5.3956, while for A3, T^L=AsinΩt^2=12. The theoretical value of time length in Figure 15 can be calculated as Λ=t^1−t^2=203.7152, which coincides with the numerical result Λ1≈197.5186 from the time history (Figure 12(b)).Figure 15
The time history of slow-varying parameterT^L.


## 2.6.7. Calculation of Natural and Spiking Frequency

The natural frequencyΩ0 of the system can be estimated by the imaginary part of equilibrium point eigenvalue with the amplitude of external excitation A=0. The equilibrium points of two subsystems in D± can be approximated as E±4.7248,±1.1403,±4.1436. The corresponding eigenvalues of Jacobi matrix are calculated as λ1=−4.7106,λ±=−0.6447±4.3104I implying that the approximate value of natural frequency Ω0=4.3104≫0.01=Ω.The spiking oscillations in this case are caused by the rapid convergence of stable focusEf±4,±2,±2 corresponding to the maximum amplitude of external excitation A=±12. The related eigenvalues can be computed as λ1=−5.4280,λ±=−0.2860±2.4111I indicating that the theoretical value of spiking frequency is ΩT=2.4111, which agrees well with the numerical result ΩSP=2π/TSP≈2.6059≈2.4111 obtained from Figure 16.Figure 16
The time history of spiking oscillation.


## 2.6.8. Discussion on Physical Meaning

The spiking state in bursting corresponds to the high frequency oscillation in engineering practice which is undoubtedly harmful to the stability and life span of the equipment, but it can not be concluded that quiescent state is always beneficial. For instance, in engineering practice, the quiescent state caused by sliding bifurcation (friction factor) implies the stagnation or jamming of machines, which can result in the fatigue and abrasion of the internals. Therefore, it is urgent to study the mechanism of bursting oscillations with sliding bifurcation and propose a reasonable control scheme to eliminate these harmful vibrations.


## 3. Conclusions

The dynamics of autonomous PMSM system correspond to the operation mode of motor without external load, which can reflect its inherent properties. The results show that the appropriate parameter combination and initial state selection can effectively avoid the irregular oscillations caused by chaos and quasi-chaotic convergence, thus reducing the risk of equipment damage.Due to the order gap between the natural frequency and the excitation frequency, the bursting oscillation may occur in nonautonomous PMSM system resulting in severe local vibration of the internal components. Three modes of bursting are found with the mechanism analysis, which indicates that internal parameters and the external excitation can both influence the attractor structure qualitatively. Besides, the amplitude of external excitation can directly determine the range of equilibrium branch involved in bursting oscillations. In particular, the introduction of friction factor may destroy the smoothness of the vector field and cause sliding bifurcation implying the intermittent pause of the motor which is ought be eliminated. The results show that the sliding motion has special structures of attractor and vector field.It is ought to be supplemented that the nonsmooth bursting mode in high-dimensional case is not proposed in this paper due to current research progress. In addition, this paper only involves the external excitation without the parametric excitation, which is also of importance. For instance, friction can lead to the shaft looseness and damage, resulting in the periodical change of rotor’s moment of inertia, which can be expressed as the form of parametric excitation. These problems requiring further study will be discussed separately in other literature.



---

*Source: 9944286-2021-07-12.xml*