Analysis of torsional stability and the dynamical characteristics of electromechanical systems

Introduction

Permanent-magnet synchronous motor (PMSM) is of significant interest, particularly for industrial applications such as in the hybrid electric vehicles (HEVs) in the medium-power range owing to its superior features such as compact size, high torque/weight ratio, high torque/inertia ratio, and the absence of rotor losses. The electromechanical coupling phenomenon has been observed in most of the systems. Therefore, recently, coupling has been reported to play a crucial role in determining important evolutions of nonlinear systems.^1–4 The nonlinear dynamic characteristics of the nonlinear electromechanical coupling systems have been extensively investigated and reported in the literature.^5–9 However, in these reports, only the bifurcation and the chaos of the three-dimensional (3D) or four-dimensional (4D) electromechanical coupling systems were analyzed. Few literature reports investigated the nonlinear characteristics of the high-dimensional electromechanical coupling systems. Xue ¹⁰ proposed a methodology called SPTR for stability-preserving trajectory reduction. The key is to get the time response curves at first, and divide the curves into complementary pairs respectively, and then map them into a series of spaces with 2 degrees of freedom (DOF) using suitable linear transformation with stability preserving. Based on this methodology, Shao ¹¹ proposed the coordinate plane projection (CPP) method to study structure stability of general dynamic system. It divides the system into a pair of complementary subsets, namely, subset X₁ of R² and its complementary subset X₂ of Rⁿ⁻².

Over the past few years, the secure and stable operation of the PMSM, which is an essential requirement of industrial automation manufacturing, has received considerable attention. For instance, Li et al. ¹² studied the dynamic characteristics of the PMSM based on the bifurcation and chaos theory. With changing operating parameters, the PMSM exhibits complex behaviors such as limiting cycles and chaos oscillation. Jing et al. ¹³ thoroughly investigated the complex dynamics in the PMSM with a nonsmooth air gap, extending the work on the smooth cases. Their theoretical analysis and numerical simulations show some more new results in the PMSM including period-doubling bifurcation, cyclic-fold bifurcation, single- and double-scroll chaotic attractors, ribbon chaotic attractors, as well as intermittent chaos. A linear transformation was used to develop a more explicit form of the PMSM model, and the stability and chaotic behavior of the resulting system were analyzed. ¹⁴ Zaher ¹⁵ analyzed the stability of the PMSM in the uncontrolled and controlled cases. Hemati and Kwatny ¹⁴ presented a compact representation of the machines’ dynamics, considering the nonlinear dynamic characteristics of variable-speed permanent-magnet machines. The systems under investigation were shown to possess multiple equilibria, which profoundly affected their global stability and dynamic characteristics.

The abovementioned literature analyzed the nonlinear characteristics of the PMSMs. In recent years, based on the theoretical analysis of the electromechanical coupling system, a lot of research achievements about bifurcation and chaos of an electromechanical coupling system have been reported. Liu et al. ¹⁶ deduced the dynamic equation of a nonlinear electromechanical coupling relative rotation system using the dissipation Lagrange equation. The necessary and sufficient conditions of Hopf bifurcation are given by choosing the electromagnetic stiffness as a bifurcation parameter. Wei et al. ¹⁷ investigated the collective dynamics of two coupled PMSMs, indicating that with positive and negative couplings, the coupled PMSMs display rich phenomena such as complete synchronization, amplitude death, antiphase synchronization, and phase-flip transition. R Darula and S Sorokin ³ used the multiple scales method to find a steady-state response of the electromechanical system exposed to a harmonic close-resonance mechanical excitation. Taking into account the nonlinearity of the generator shaft and the interaction of mechanics and electrics in the generator set, Niu and Qiu ¹⁸ obtained a transient model by combining the Park equations^19,20 of the synchronous generator and mechanical equations. The Hopf bifurcation and the chaos were investigated by the nonlinear mode and Floquet theory.

Above all, there are researchers studying dynamic behaviors of high-dimension system, but that is not the electromechanical coupling system. And there are researchers studying nonlinear dynamic behavior in the electrical motor system, but that is a simplified 3-DOF model. In this study, the nonlinear electromechanical coupling model was established considering the coupling effects of both electrical and mechanical parameters, leading to a five-dimensional (5D) autonomous equation with only two quadratic terms. Because of the large dimensional system, the dynamic characteristics cannot be directly solved and deduced. Based on the CPP method, the dynamic characteristics of the system were analyzed. By calculating the weight of all the 10 projection planes of the electromechanical coupling system, four projection planes with less weight values were chosen and analyzed. By analyzing the four projection planes of the system, the parameters of the system determining desired dynamic characteristics such as stable state were investigated. Finally, computer simulations were performed to verify the presented dynamic behaviors in the electromechanical coupling system.

Model and equations

The transmission system scheme in HEV is shown in Figure 1(a), in which k₁ is the characteristic value of the first planet gear, k₂ is the characteristic value of the second planet gear, i indicates the input shaft, and o the output shaft.

OPEN IN VIEWER

Figure 1.

Transmission and dynamic system: (a) transmission scheme and (b) electromechanical coupling system.

Figure 1(a) shows the variable state combination of brake Z, and clutch C could provide variable transmission mode. When the system is operated in one mode, the motor A or B could be coupled to the mechanical system. Not considering the transmission error and speed ratio, the engine excitation and the road excitation are the load excitation for the simplified electromechanical model.

The simplified electromechanical model of an electromechanical coupling system is shown in Figure 1(b). The electrical part is a PMSM, and the mechanical part has the load excitation indicated by T_l. The two parts are coupled using a shaft, whose stiffness and damp are k and β, respectively.

Let φ₁ be the angle of rotation in PMSM and φ₂ be the angle of rotation in the load system. Then, the torsional angle of the coupling shaft is α = φ₁ − φ₂.

In this study, the transient model of the coupling system was obtained by combining the Park equations^19,20 with mechanical equations

{ i · d = u d L d + L q L d n p ω 1 i q − 1 L d R 1 i d i · q = u q L q − L d L q n p ω 1 i d − n p ω 1 ψ f L q − 1 L q R 1 i q ω · 1 = 1 I 1 { 1 . 5 n p [ ψ f i q + ( L d − L q ) i d i q ] − k α − β ( ω 1 − ω 2 ) } ω · 2 = 1 I 2 ( k α + β ( ω 1 − ω 2 ) − T l ) α · = ω 1 − ω 2

(1)

where i_d, i_q, ω₁, ω₂, and α are the state variables, representing the d-axis stator current, q-axis stator current, motor rotary speed, load rotary speed, and torsional angle, respectively; u_d and u_q are the d-axis and q-axis stator voltage, respectively; I₁ is the rotor inertia of the PMSM; I₂ is the rotor inertia of the load; k is the stiffness of the connecting shaft; β is the damp of the connecting shaft; R₁ is the stator winding resistance; T_l is the external load torque; L_d and L_q are the d-axis and q-axis stator inductors, respectively; ψ_f is the permanent-magnet flux; and n_p is the number of pole-pairs.

Electrical vehicles mainly use surface-mounted PMSMs. Then, L_d = L_q = L. Ignoring the damp of the connecting shaft, namely β = 0, the model equation can be written as

{ i · d = u d L + n p ω 1 i q − 1 L R 1 i d i · q = u q L − n p ω 1 i d − n p ω 1 ψ f L − 1 L R 1 i q ω · 1 = 1 I 1 ( 1 . 5 n p ψ f i q − k α ) ω · 2 = 1 I 2 ( k α − T l ) α · = ω 1 − ω 2

(2)

The equilibrium point can be calculated by comparing each equation in system of equation (2) to zero

{ 0 = u d L + n p ω 1 i q − 1 L R 1 i d 0 = u q L − n p ω 1 i d − n p ω 1 ψ f L − 1 L R 1 i q 0 = 1 I 1 ( 1 . 5 n p ψ f i q − k α ) 0 = 1 I 2 ( k α − T l ) 0 = ω 1 − ω 2

(3)

Then, the following equation can be obtained

{ i d = u d R 1 + n p L ω 1 R 1 T l 1 . 5 n p ψ f i q = T l 1 . 5 n p ψ f 0 = n p 2 L 2 T l 1 . 5 n p ψ f ω 1 2 + ( n p L u d + R 1 n p ψ f ) ω 1 + R 1 2 T l 1 . 5 n p ψ f − R 1 u q ω 2 = ω 1 α = T l k

(4)

The Jacobian matrix can be obtained as

D x = [ − R 1 L n p ω 1 n p i q 0 0 − n p ω 1 − R 1 L − n p i d − n p ψ f L 0 0 0 1 . 5 n p ψ f I 1 0 0 − k I 1 0 0 0 0 k I 2 0 0 1 − 1 0 ]

(5)

The equilibrium point can be obtained from equation (4). Substituting equation (4) into equation (5) and solving the eigenvalue of the Jacobian matrix, the dynamic stability of the system can be analyzed.

Dynamic behaviors of the electromechanical coupling system

First solving the characteristic equation of the corresponding Jacobian matrix (5), the solution of x₁, x₂, …, x_n could be obtained. If at least one of the solutions has a positive part, then the bifurcation and chaos occur.

Solving equation (3), assuming a = n_pL²T_l/1.5ψ_f, b = n_pLu_d + R₁n_pψ_f, and c = R 1 2 T l / 1 . 5 n p ψ f − R 1 u q , leads to the equilibrium point as represented by equation (6)

x 1 , 2 = { u d R 1 + n p L ω 1 R 1 T l 1 . 5 n p ψ f T l 1 . 5 n p ψ f − b ± b 2 − 4 ac 2 a − b ± b 2 − 4 ac 2 a T l k

(6)

The Jacobian matrix is the same as represented by equation (5), and equation (7) represents its characteristic equation

[ − R 1 L − x n p ω 1 n p i q 0 0 − n p ω 1 − R 1 L − x − n p i d − n p ψ f L 0 0 0 1 . 5 n p ψ f I 1 − x 0 − k I 1 0 0 0 − x k I 2 0 0 1 − 1 − x ] = 0

(7)

Because of the complicated solution of the characteristic equation, it is very difficult to solve equation (7) directly. Therefore, the stability of the system cannot be analyzed directly, and a new method must be proposed.

A novel methodology, for trajectory stability-preserving dimension reduction (TSPDR), was introduced. The key of the TSPDR is to get the time response curves at first and then divide the curves into complementary pairs and map them into a series of spaces with 2 DOF using suitable linear transformation with preserved stability. Based on the TSPDR, the CPP^10,11 method was applied for studying the structural stability of the electromechanical coupling system. It divides the system into a pair of complementary subsets, namely, subset X₁ of R² and its complementary subset X₂ of Rⁿ⁻². On the coordinate plane X₁, the equation of X₁ can be transformed to a pseudo-linear system with one equilibrium point, as X₂ is handled as a time-varying parameter. At each time t, the characteristics of the equilibrium point can be known by the eigenvalue technology, and the corresponding phase trajectory can be sorted into five modes: stable node, stable focus, unstable focus, unstable node, and saddle. The critical condition for the mode to switch from one to another, namely the bifurcation set B, can be obtained analytically. Bifurcation occurs whenever X intersects B, thereby the point corresponding to the bifurcation is identified. Moreover, the route to chaos via a sequence of bifurcations can also be thoroughly understood.

Shao ¹¹ analyzed the 5-DOF Lorenz system by the CPP method. From the 10 subsets, 4 subsets were chosen to analyze the chaos of the system. The analysis of these four subsets, assuming |x₁(t)| < sqrt(15)/3, leads all the variables to be in stable status.

With these advantages of the CPP method, this study used the CPP method^10,11 to analyze the dynamic characteristics of the electromechanical coupling system.

For the system of n variables, the arbitrarily two variables could form a plane’s trajectory, and there are a total n(n + 1)/2 planes. If the characteristics of the whole n variables need to be analyzed in a coordinate plane set, then the following conditions for the coordinate plane set must be satisfied: (1) the coordinate plane set contains all the n variables and (2) the coordinate plane set should contain as few coordinate planes as possible. Namely, the arbitrarily n − 1 mutual independent coordinate planes contain the whole dynamic information.

In this study, n = 5, then n − 1 = 4. Therefore, four 2D systems must be chosen to reflect the original 5D system.

Analysis of the electromechanical 5-DOF system using CPP method

The electromechanical coupling system is shown by equation (2), and the number of all the projection planes is n(n − 1)/2 = 10. Therefore, the first 10 projection planes were analyzed to obtain the bifurcation set and the motion pattern. ¹¹

The motion pattern partition on the projection plane i_d-i_q

The projection of equation (2) on the i_d and i_q plane can be written as equation (13) in Appendix 1.

From this equation, the projection equilibrium point can be obtained

( u d R 1 + n p ω 1 u q R 1 − ( n p ω 1 ) 2 ( u d L R 1 2 + ψ f R 1 ) ( n p ω 1 ) 2 L R 1 + R 1 L , u q L − n p ω 1 ( u d R 1 + ψ f L ) ( n p ω 1 ) 2 L R 1 + R 1 L )

and the equilibrium point moves along the curve: n_pω₁Li_q + u_d − R₁i_d = 0. The characteristic equation is (−R₁/L−λ)² + (n_pω₁)² = 0, and the eigenvalues are λ_1,2 = −R₁/L ± in_pω₁.

The bifurcation set and the motion pattern of the projection phase trajectory on the projection planes i_d-i_q are as follows

{ ω = 0 stable node ω ≠ 0 stable focus

(8)

2. The motion pattern partition on the projection plane i_d-ω₁

The projection of equation (2) on the i_d-ω₁ plane can be written as equation (14) in Appendix 1.

The projection equilibrium point (i_d, ω₁) is on the curve u_d/L + n_pω₁i_q−R₁i_d/L = 0, when i_q = kα/1.5n_pψ_f, and there is no equilibrium point when i_q ≠ kα/1.5n_pψ_f. The characteristic equation is −λ(−R₁/L − λ) = 0, and the corresponding eigenvalues are λ₁ = −R₁/L and λ₂ = 0.

Above all, the motion pattern of the projection phase trajectory on the projection plane i_d-ω₁ is non-hyperbolic equilibrium position.

3. The motion pattern partition on the projection plane i_d-ω₂

The projection of equation (2) on the plane i_d-ω₂ can be written as equation (15) in Appendix 1.

When the right-hand side of the equations are ordered to be zero, solutions could be obtained that i_d = (u_d + Ln_pω₁i_q)/R₁. And the variable ω₂ cannot be solved. Thus, the projection equilibrium point is {(u_d + Ln_pω₁i_q)/R₁, cannot be calculated}. The characteristic equation is −λ(−R₁/L−λ) = 0, and the eigenvalues are λ₁ = −R₁/L and λ₂ = 0.

Above all, the motion pattern of the projection phase trajectory on the projection plane i_d-ω₂ is in non-hyperbolic equilibrium position.

4. The motion pattern partition on the projection plane i_d-α

The projection of equation (2) on the plane i_d-ω₂ can be written as equation (16) in Appendix 1.

When the right-hand side of the equations are ordered to be zero, solutions could be obtained that i_d = (u_d + Ln_pω₁i_q)/R₁. But the variable α cannot be solved. The projection equilibrium point is {(u_d + Ln_pω₁i_q)/R₁, cannot be calculated}, and the characteristic equation is −λ(−R₁/L−λ) = 0, and the eigenvalues are λ₁ = −R₁/L and λ₂ = 0.

Above all, the motion pattern of the projection phase trajectory on the projection plane i_d-α is in non-hyperbolic equilibrium position.

5. The motion pattern partition on the projection plane i_q-ω₁

The projection of equation (2) on the plane i_q-ω₁ can be written as equation (17) in Appendix 1.

The projection equilibrium point is {kα/1.5n_pψ_f, (u_q−R₁kα/1.5n_pψ_f)/(n_pψ_f + Ln_pi_d)}, and the projection equilibrium point (i_q, ω₁) is on the curve u_q/L − n_pω₁i_d − n_pω₁ψ_f/L − R₁i_q/L = 0; the characteristic equation is λ² + R₁λ/L + 1.5n_pψ_f (n_pi_d + n_pψ_f/L)/I₁ = 0, and the eigenvalues are λ_1,2 = −R₁/2L ± sqrt[(R₁/L)²−6(n_pi_d + n_pψ_f/L)n_pψ_f/I₁]/2.

The bifurcation set and the motion pattern of the projection phase trajectory on the projection plane i_q-ω₁ is defined by the following equation.

6. The motion pattern partition on the projection plane i_q-ω₂

The projection of equation (2) on the plane i_q-ω₂ can be written as equation (19) in Appendix 1.

When the right-hand side of the equations are ordered to be zero, solutions could be obtained that i_q = (u_q − Ln_pω₁i_d − n_pω₁ψ_f)/R₁. But the variable ω₂ cannot be solved. The projection equilibrium point is {(u_q − Ln_pω₁i_d − n_pω₁ψ_f)/R₁}, which cannot be calculated}, and the projection equilibrium point exists when α = T_l/k and doesn’t exist when α ≠ T_l/k. The characteristic equation is −λ(−R₁/L−λ) = 0, and the eigenvalues are λ₁ = −R₁/L, λ₂ = 0.

Above all, the motion pattern of the projection phase trajectory on the projection plane i_q-ω₂ is in the non-hyperbolic equilibrium position.

7. The motion pattern partition on the projection plane i_q-α

The projection of equation (2) on the plane i_q-α can be written as equation (20) in Appendix 1.

When the right-hand side of the equations are ordered to be zero, solutions could be obtained that i_q = (u_q − Ln_pω₁i_d − n_pω₁ψ_f)/R₁. But the variable α cannot be solved. The projection equilibrium point is {(u_q − Ln_pω₁i_d − n_pω₁ψ_f)/R₁}, which cannot be calculated, and the projection equilibrium point exists when ω₁ = ω₂ and doesn’t exist when ω₁ ≠ ω₂. The characteristic equation is −λ(−R₁/L−λ) = 0, and the eigenvalues are λ₁ = −R₁/L, λ₂ = 0.

Above all, the motion pattern of the projection phase trajectory on the projection plane i_q-α is in non-hyperbolic equilibrium position.

8. The motion pattern partition on the projection plane ω₁-ω₂

The projection of equation (2) on the plane ω₁-ω₂ can be written as equation (21) in Appendix 1.

In equation (21), the variables i_q and α are the known numbers. Therefore, the projection equilibrium point cannot be calculated, and the projection equilibrium point exists only when α = T_l/k = 1.5n_pψ_fi_q/k; otherwise, it does not exist. The characteristic equation is λ² = 0, and the eigenvalues are λ_1,2 = 0.

Above all, the motion pattern of the projection phase trajectory on the projection plane ω₁-ω₂ is in non-hyperbolic equilibrium position.

9. The motion pattern partition on the projection plane ω₁-α

The projection of equation (2) on the plane ω₁-α can be written as equation (22) in Appendix 1.

The projection equilibrium point is {ω₂, 1.5n_pψ_fi_q/k}. The characteristic equation is λ² + k/I₁ = 0, and the eigenvalues are λ_1,2 = ±√(k/I₁).

Above all, the motion pattern of the projection phase trajectory on the projection plane ω₁-α follows the center point motion pattern.

10. The motion pattern partition on the projection plane ω₂-α

The projection of equation (2) on the plane ω₂-α can be written as equation (23) in Appendix 1.

The projection equilibrium point is (ω₁,T_l/k). The characteristic equation is λ² + k/I₂ = 0, and the eigenvalues are λ_1,2 = ±√(k/I₂).

Above all, the motion pattern of the projection phase trajectory on the projection plane ω₁-α follows center point motion pattern.

Projection plane set of the electromechanical coupling system

According to the principle in the literature, ¹¹ the complex degree, the weight value of the projection equilibrium point, and the bifurcation set on every projection plane are listed in Table 1.

OPEN IN VIEWER

Table 1.

The CPP analysis result of the electromechanical coupling system.

Projection plane	Equilibrium point	ωij, e	Bifurcation curve	ωij, b	ωij
id-iq	{ u d R 1 + n p ω 1 u q R 1 − ( n p ω 1 ) 2 ( u d L R 1 2 + ψ f R 1 ) ( n p ω 1 ) 2 L R 1 + R 1 L , u q L − n p ω 1 ( u d R 1 + ψ f L ) ( n p ω 1 ) 2 L R 1 + R 1 L }	2	{ ω 1 = 0 stable node ω 1 ≠ 0 stable focus	3	2.2
id-ω1	On the curve ud/L + npω1iq−R1id/L = 0	3.2	Non-hyperbolic equilibrium position	7	3.96
id-ω2	{(ud + Lnpω1iq)/R1, can’t be calculated}	3.4	Non-hyperbolic equilibrium position	7	4.12
id-α	{(ud + Lnpω1iq)/R1, can’t be calculated}	3.4	Non-hyperbolic equilibrium position	7	4.12
iq-ω1	{kα/1.5npψf, (uq−R1kα/1.5npψf)/(npψf + Lnpid)}	3.2	{ i d < − ψ f L saddle point i d > − ψ f L stable node	6	3.76
iq-ω2	{(uq−Lnpω1id−npω1ψf)/R1, can’t be calculated}	3.4	Non-hyperbolic equilibrium position	7	4.12
iq-α	{(uq−Lnpω1id−npω1ψf)/R1, can’t be calculated}	3.4	Non-hyperbolic equilibrium position	7	4.12
ω1-ω2	{can’t be calculated, can’t be calculated}	4	Non-hyperbolic equilibrium position	7	4.6
ω1-α	{ω2, 1.5npψfiq/k}	3.4	Center point	6	3.92
ω2-α	{ω1, Tl/k)	3.2	Center point	6	3.76

The variable ω_ij, _e is defined as following in the literature ¹¹

ω ij , e = { 1 The projection equilibrium point is fixed 2 The projection equilibrium point moves along with a certain curve 3 + m / n The projection equilibrium point moves along with an uncertain curve

where m is the number of variables, which are in the projection equilibrium equation.

ω_ij, _b indicates the dimension of the bifurcation set of the projection systems. And finally, ω_ij = qω_ij, _e  + (1 − q)ω_ij, _b , where the weight coefficient q indicates the extent of the effect on the projection subsystems of the equilibrium point movement and the bifurcation set complexity. The weight coefficient q was set as 0.8. The final chosen result of the projection plane set is shown in Figure 2.

OPEN IN VIEWER

Figure 2.

The optimized CPP analysis plane set of the system.

Principle of the chaotic phenomenon in the electromechanical coupling system

Chaos occurs in the dynamic trajectory of the electromechanical coupling system, and the system parameters are listed in Table 2. The principle of the chaotic trajectory was analyzed in the determined projection set P = {i_d-i_q, i_q-ω₁, ω₁-α, and ω₂-α}, then the dynamic characteristics can be revealed.

OPEN IN VIEWER

Table 2.

System parameters.

R1 (Ω)	0.1	ψf (Wb)	0.63
Ld = Lq (H)	0.091	I1 (kg m2)	0.041
I2 (kg m2)	0.001	k (N/rad)	3e4
Tl (N m)	50	np	3
ud (V)	380	uq (V)	50

Substituting the value of variables in Table 2 into equation (7), characteristic roots are obtained: λ = (0.0006 + 5.5467i, 0.0006 − 5.5467i, −0.5488 + 1.0484i, −0.5488 − 1.0484i, 1.0989) × 10³. According to the five values of λ, there are three positive values. Therefore, the system is unstable. Thus, it must be analyzed by the CPP method.

1. The principle on the projection plane i_d-i_q

According to equation (8), when ω ≠ 0, the projection phase trajectory is at the stable focus state. Therefore, the phase trajectory on the projection plane i_d-i_q moves around the center point, which is the equilibrium. The equilibrium point moves along with the curve n_pω₁Li_q + u_d − R₁i_d = 0; therefore, the center point moves along with the curve n_pω₁Li_q +u_d − R₁i_d = 0.

According to the CPP analysis on the projection plane i_d-i_q, the equilibrium point moves along the curve n_pω₁Li_q + u_d − R₁i_d = 0. The equilibrium point is the stable node when ω₁ = 0, and the corresponding phase trajectory follows the stable node motion pattern. The equilibrium point is the stable focus point when ω ≠ 0, and the corresponding phase trajectory follows the stable focus point motion pattern. Figure 3 shows that during the whole time range, ω₁ ≠ 0, and the phase trajectory follows the stable focus point motion pattern. The phase trajectory’s center point moves with the equilibrium point. Figure 4(a) shows that the equilibrium point moves along with a line when t∈ (42,42.014), and therewith the phase trajectory’s center point. When the machine rotary speed ω₁ becomes larger, namely t∈ (82,82.004), ω₁ fluctuates, leading to the fluctuation of i_d and i_q. The phase trajectory’s center point did not move, since the equilibrium point fluctuates only slightly, as shown in Figure 4(b).

OPEN IN VIEWER

Figure 3.

t–ω₁ curve.

OPEN IN VIEWER

Figure 4.

i_d-i_q phase trajectory: (a) t∈ (42,42.014) phase trajectory and (b) t∈ (82,82.004) phase trajectory.

2. The principle on the projection plane i_q-ω₁

On the projection plane i_q-ω₁, after a period of time, the phase trajectory moves following a circle, as shown in Figure 5.

OPEN IN VIEWER

Figure 5.

Phase trajectory on the projection plane i_q-ω_1.

According to the CPP analysis on the projection plane i_q-ω₁, the phase trajectory follows the saddle point motion pattern when i_d < −ψ_f/L and the stable node motion pattern when i_d > −ψ_f/L. Figure 6(a) shows that when t∈ (10,70), the phase trajectory moves away from the equilibrium point because i_d < −ψ_f/L during the period. Figure 6(b) shows that when t∈ (85,85.001), the phase trajectory switches between the saddle point motion pattern and the stable node motion pattern because i_q fluctuates on the point −ψ_f/L.

OPEN IN VIEWER

Figure 6.

Phase trajectory: (a) phase trajectory t∈ (10,70) and (b) phase trajectory t∈ (85,85.004).

3. The principle on the projection plane ω₁-α

Figure 7 shows that the phase trajectory of the projection plane ω₁-α follows a rotational motion.

OPEN IN VIEWER

Figure 7.

Phase trajectory on the projection plane ω₁-α.

According to the CPP analysis on the projection plane ω₁-α, the phase trajectory follows the center point motion pattern. Figure 8(a) shows the phase trajectory during t∈ (32,32.005). The phase trajectory’s rotary center moves, because the equilibrium point is on a horizontal line and then appears spiral as the equilibrium moves, and finally appears as a limiting circle, as shown in Figure 8(b).

OPEN IN VIEWER

Figure 8.

Phase trajectory: (a) phase trajectory during t∈ (32,32.005) and (b) phase trajectory during t∈ (85,85.005).

4. The principle on the projection plane ω₂-α

The phase trajectory is on the rotational motion on the projection plane ω₂-α, as shown in Figure 9.

OPEN IN VIEWER

Figure 9.

Phase trajectory on the projection plane ω₂-α.

According to the CPP analysis on the projection plane ω₁-α, the phase trajectory follows the center point motion pattern. Figure 10(a) shows the phase trajectory during t∈ (32,32.005). The phase trajectory’s rotary center moves because the equilibrium point is on a horizontal line and then appears spiral. During t∈ (85,85.005), the phase trajectory’s center moves as the equilibrium moves and then appears to be a limiting circle, as shown in Figure 10(b).

OPEN IN VIEWER

Figure 10.

Phase trajectory: (a) phase trajectory during t∈ (32,32.005) and (b) phase trajectory during t∈ (85,85.005).

Dynamic behavior analysis of the electromechanical coupling system

Using the CPP method, not only the dynamic behavior but also the special trajectory in a nonlinear system could be analyzed. The analysis on the projection plane set P = {i_d-i_q, i_q-ω₁, and ω₁-α, ω₂-α} shows that the bifurcation set structure and the moving rule of equilibrium on an arbitrary projection plane is relative to other state variables. The interaction between the five state variables forms the dynamic behavior of the electromechanical coupling system. The complexity of the projection subsystem is relatively small on the projection planes i_d-i_q and ω₂-α. Since the bifurcation set on the two projection planes has the state variable ω₁, the whole system’s dynamic behavior could be deduced by monitoring the range of ω₁ from the numerical integration.

Variable ω₁ converging at a constant C

Since ω₁ = constant C, the phase trajectory on the projection plane i_d-i_q follows the stable node or stable focus motion pattern mode and finally converges at the equilibrium

{ u d R 1 + n p ω 1 u q R 1 − ( n p ω 1 ) 2 ( u d L R 1 2 + ψ f R 1 ) ( n p ω 1 ) 2 L R 1 + R 1 L , u q L − n p ω 1 ( u d R 1 + ψ f L ) ( n p ω 1 ) 2 L R 1 + R 1 L } lim t → ∞ i d ( t ) = u d R 1 + n p ω 1 u q R 1 − ( n p ω 1 ) 2 ( u d L R 1 2 + ψ f R 1 ) ( n p ω 1 ) 2 L R 1 + R 1 L , lim t → ∞ i q ( t ) = u q L − n p ω 1 ( u d R 1 + ψ f L ) ( n p ω 1 ) 2 L R 1 + R 1 L

From the equation ω · 1 = ( 1 / I 1 ) ( 1 . 5 n p ψ f i q − k α ) , the following equation can be obtained: lim t → ∞ α ( t ) = ( 1 . 5 n p ψ f i q ) / k . Then, combining with the equation α · = ω 1 − ω 2 , the following equation can be obtained: lim t → ∞ ω 2 ( t ) = ω 1 .

Above all, as long as the dynamic behavior of electromechanical coupling system satisfies ω₁ = constant C, the other state variables meet the following equations

lim t → ∞ i d ( t ) = u d R 1 + n p ω 1 u q R 1 − ( n p ω 1 ) 2 ( u d L R 1 2 + ψ f R 1 ) ( n p ω 1 ) 2 L R 1 + R 1 L

(9)

lim t → ∞ i q ( t ) = u q L − n p ω 1 ( u d R 1 + ψ f L ) ( n p ω 1 ) 2 L R 1 + R 1 L

(10)

lim t → ∞ α ( t ) = 1 . 5 n p ψ f i q k

(11)

lim t → ∞ ω 2 ( t ) = ω 1 = constant C

(12)

Figure 11 illustrates the simulation result with the system parameters. The data in Table 3 indicate that the system is stable.

OPEN IN VIEWER

Figure 11.

Simulation result of the system variables.

OPEN IN VIEWER

Table 3.

System parameters.

R1 (Ω)	6.1	ψf (Wb)	0.63
Ld = Lq (H)	0.011	I1 (kg m2)	0.001
I2 (kg m2)	0.081	K (N/rad)	6e4
Tl (N m)	50	np	3
ud (V)	191	uq (V)	50

Substituting the value of variables in Table 3 to equation (7), characteristic roots were obtained: λ = (−6.7028 + 0.0677i, −6.7028 − 0.0677i, −0.0003 + 0.0442i, −0.0003 − 0.0442i, −0.0003) × 10⁴. The real parts of λ are all negative values. Therefore, the system is stable, and it is also proved, as shown in Figure 11.

Variable ω₁ fluctuating the constant C

If the value of ω₁ fluctuates, the phase trajectory on the projection plane follows the center point motion pattern. Moreover, the trajectory center point fluctuates with ω₁. Thereby, the 5-DOF system could form the helical trajectory, finally leading to chaos.

Results and discussion

In this section, the CPP method is presented to reduce the 5-DOF system to the 2-DOF system. In order to complete the system, the dynamic characteristics of four projection planes were analyzed, when the other three variables change. Moreover, the complexity of the projection subsystem is relatively small on the projection planes i_d-i_q and ω₂-α. Since the bifurcation set on the two projection planes has the state variable ω₁, the dynamic behavior of whole system could be deduced by monitoring the range of ω₁ from the numerical integration. Therefore, the dynamic stability of the two projective planes was analyzed when the value of ω₁ tends to be stable and fluctuate. Moreover, the results of this study show that the electromechanical system will be stable if the ω₁ value is controlled in a stable area.

Conclusion

In conclusion, the electromechanical coupling system was investigated using a 5-DOF dynamic model. The stability of the model was studied, considering the effect of electromechanical coupling. Because the model which contains the nonlinear characteristics is large, the CPP method was used to reduce the system’s number of 5-DOF to 2-DOF. The theory of nonlinear stability was applied to analyze the system. The following conclusions were drawn:

This study used the CPP method, which is only used in the power system up to now. Through this CPP method, the 5-DOF system was successfully reduced to a 2-DOF system. Moreover, the conditions under which the coupling system exhibits specific behavior are presented. Extending this study to the control area has a large scope. Some control strategies for the electromechanical coupling system could be built to avoid or utilize the chaos of the system.