- Email: [email protected]
Regular self-oscillating and chaotic behaviour of a PID controlled gimbal suspension gyro
Chaos, Solitons and Fractals 21 (2004) 1057–1074 www.elsevier.com/locate/chaos
Regular self-oscillating and chaotic behaviour of a PID controlled gimbal suspension gyro Manuel F. Perez Polo
a,*
, Manuel Perez Molina
b,1
a
b
Department of Fısica, Ingenierıa de Sistemas y Teorıa de la Se~nal, Universidad de Alicante, Escuela Politecnica Superior, Campus de San Vicente, Alicante 03071, Spain Facultad de Ciencias Matematicas, Universidad Nacional de Educacion a Distancia, UNED, C/Boyero 12-1A, Alicante 03007, Spain Accepted 6 January 2004
Abstract The dynamics of a gyro in gimbal with a PID controller to obtain steady state, self-oscillating and chaotic motion is considered in this paper. The mathematical model of the whole system is deduced from the gyroscope nutation theory and from a feedback control system formed by a PID controller with constrained integral action. The paper shows that the gyro and the associated PID feedback control system have multiple equilibrium points, and from the analysis of a Poincare–Andronov–Hopf bifurcation at the equilibrium points, it is possible to deduce the conditions, which give regular and self-oscillating behaviour. The calculation of the first Lyapunov value is used to predict the motion of the gyro in order to obtain a desired equilibrium point or self-oscillating behaviour. The mechanism of the stability loss of the gyro under small vibrations of the gyro platform and the appearance of chaotic motion is also presented. Numerical simulations are performed to verify the analytical results. 2004 Elsevier Ltd. All rights reserved.
1. Introduction Throughout the years, the dynamics of gyroscopic devices has been considered as one of the most exciting problems of mechanics. From the classical handbook by Klein and Sommerfield [1] at the beginning of the last century until the present time, many books and papers concerning the stability, nonlinear phenomena and technical applications of gyroscopic systems have been be written [2–6]. Gyros for sensing angular motion are used in airplane automatic pilots, rocket-vehicle launch-guidance, space-vehicle attitude systems, ship’s gyrocompasses and submarine inertial autonavigators. Today, fibber optic gyro based has many applications in avionics, platform stabilization and automotive testing. In balanced gimbal suspension gyro systems, it is well known that the exact equations which describe the nutation motion are nonlinear [7,8]. However, for purpose of simplified analysis, a linear model has been deduced in order to obtain a transfer function that can be used in the design of a controller. Exact equations of motion from an ideal twoaxis gyro, the linear approximation for small motions, and the response to a constant torque to give precession and nutation motion can be found in [9]. Singh [10], Ge and Chen [11], taking into account harmonic external disturbances have already analyzed the nonlinear dynamics of a single-axis gyro mounted on a space-vehicle. These papers show that from averaging method [12,13] and the harmonic balance it is possible to obtain chaotic behaviour. A similar system, but with feedback
*
Corresponding author. Tel.: +34-9-659-09673; fax: +34-9-65-909750. E-mail addresses: manolo@dfists.ua.es (M.F. Perez Polo), [email protected] (M.P. Molina). 1 Tel.: +34-9-651-01988; fax: +34-9-659-09750. 0960-0779/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.01.011
1058
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
proportional control was considered in [14], using normal forms and Melnikov method [12,13]. This paper shows that regular and chaotic motions are both possible. In the paper of Agafonov [15], the stability of stationary motion, the periodic motions from a Poincare–Andronov–Hopf bifurcation and the chaotic behaviour of a gimbal suspension gyro is researched. However, in some previously cited papers, the gyro mathematical model is very simplified without taking into account the control problem, and in others papers a very simple proportional control is only used to obtain a desired set point. So, the whole effect of a PID linear controller dynamics and the exact nonlinear mathematical model of the gyro (nutation motion) are completely unknown. This paper considers the motion of a heavy symmetrical gyro in gimbal. The mathematical model is deduced from Euler–Lagrange equations, considering that the gyro moment of inertia is constant. A PID feedback control system is added, taking the angle of the internal gimbal ring as feedback signal to produce the error signal equal to set point minus feedback signal. The error signal is the input to an amplifier–PID controller to produce the output control signal, which is sent to a DC motor. The shaft of the DC motor is connected by means of an ideal gear train to the external gimbal. The external gimbal ring is also connected to an external body, which is subject to a constant external disturbance torque. The objective of the control system is that the external gimbal reaches a desired set point by considering the external disturbance torque. The whole system is mounted on a platform which is subject to periodic vibrations, so when the motion gyro is self-oscillating, chaotic motion can appear. The problem of the constrained integral action in the controller is also considered in this paper. It is well known that the dynamic range of actuators is limited, so when the input reaches the saturation limit, the integrator in the controller will continue integrating while the input is constrained. Consequently, the state of the integrator can reach an unacceptable high value that can produce a poor transient response. In order to avoid this undesirable effect, an anti-windup integrator is considered, so the integral action remains constant from the saturation of the actuator as soon as the actuator saturates [16–18]. When the integral action is constrained, the system has three equilibrium points P1 , P2 and P3 . Point P1 is identified as the set point, and the analysis of equilibrium points P2 and P3 reveals that two interesting bifurcation phenomena can appear. Considering a constrained value of the integral action, it is possible to prove that in the equilibrium point P2 the system has double-zero eigenvalues, and in the equilibrium point P3 the system has a pair of pure imaginary eigenvalues. Consequently, a Bogdanov–Takens bifurcation at point P2 , and a Poincare–Andronov–Hopf bifurcation with a weak focus at point P3 are the responsible of a complicated behaviour of the system [12,13,19,20]. In this paper, we show from the analysis of integral action of the PID controller and the Poincare–Andronov–Hopf bifurcation, that it is possible to predict the motion of the gyro from the first Lyapunov value. The center manifold theorem and the normal forms have been used in order to obtain the first Lyapunov value [19,20]. From a fixed value of the integral action, it is possible to prove that the system only has one reachable equilibrium point. In this case the first Lyapunov value is used to determine if this equilibrium point is a stable or unstable weak focus. When all the equilibrium points are unstable, the system reaches a self-oscillation motion. In this situation an external vibration of the platform can give chaotic behaviour, which has been researched from the sensitive dependence, Lyapunov exponents and the power spectral density.
2. Mathematical model Consider a gyroscope in gimbal schematically shown in Fig. 1. The following notation has been used: OX1 Y1 Z1 and OX2 Y2 Z2 are respectively the internal and external gimbal ring reference systems, and OXYZ is the gyro reference system. OZ1 is the gyro spin axis, OX2 is the external gimbal rotation axis, c is the gyro spin angle in the internal gimbal, b is the internal gimbal rotation angle, and a is the external gimbal rotation angle. Taking into account that the principal axis of the rotor are: A ¼ B, C and Ai , Bi , Ci (i ¼ 1; 2) are the moment of inertia of the internal and external gimbal respectively, and denoting pi , qi , ri the angular velocity respect to the OXi , OYi , OZi axis respectively, the following kinetic energies can be written: • Kinetic energy of the external gimbal: T2 ¼ 12ðA2 p22 þ B2 q22 þ C2 r22 Þ;
q2 ¼ r2 ¼ 0;
p2 ¼ a_ ) T2 ¼ 12A2 a_ 2
ð1Þ
• Kinetic energy of the internal gimbal: T1 ¼ 12ðA1 p12 þ B1 q21 þ C1 r12 Þ;
p1 ¼ a_ cos b;
T1 ¼ 12ðA1 a_ 2 cos2 b þ B1 b_ 2 þ C1 a_ 2 sin2 bÞ
_ q1 ¼ b;
r1 ¼ a_ sin b
ð2Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1059
Fig. 1. Physical model of the gimbal suspension gyro and reference frames.
• Kinetic energy of the rotor: T0 ¼ 12ðAp2 þ Bq2 þ Cr2 Þ; p ¼ a_ cos b; 2 2 2 1 _ _ 2Þ T0 ¼ ðAa_ cos b þ Bb þ Cð_c þ a_ sin bÞ
_ q ¼ b;
r ¼ c_ þ a_ sin b
ð3Þ
2
where c_ is the spin velocity of the rotor, and the dot denotes differentiation respect to time. Considering the projections of the angular velocities in the inertial reference system On0 s0 f0 and considering negligible friction, the whole system kinetic energy has the form: T ¼ 12½J ðbÞa_ 2 þ Hr b_ 2 þ Cð_c þ a_ sin bÞ2
ð4Þ
where Hr ¼ A þ B1 , and J ðbÞ is defined by the equation: J ðbÞ ¼ A2 þ ðA1 þ AÞ cos2 b þ C1 sin2 b
ð5Þ
Eq. (5) has the following physical interpretation: Without the term A cos2 b, the value of J ðbÞ represents the total moment of inertia of the external and internal gimbal rings respect to the rotation axis OX2 . Considering that a, b, c are independent and holonomic variables, and applying Euler–Lagrange equations [5,7,8] the gyro motion equations without control are the following: 9 d da dc da > > J ðbÞ þ C þ sin b sin b ¼ 0 > > > dt dt dt dt > > > = 2 2 d b 1 dJ ðbÞ da dc da da ð6Þ Hr 2 C þ sin b cos b ¼ 0 > dt 2 db dt dt dt dt > > > > > d dc da > > C þ sin b ¼0 ; dt dt dt Eqs. (6) represent the exact gyroscope nutation theory [5,7]. Setting A1 þ A ¼ C1 in Eq. (5), Eq. (6) can be simplified as follows: 9 d da > J þ H sin b ¼ 0 > = dt dt ð7Þ 2 > db da > Hr 2 H cos b ¼ 0 ; dt dt
1060
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 2. Gyroscope and PID controller feedback control system.
where H ¼C
dc da þ sin b ; dt dt
J ¼ A þ A1 þ A2
ð8Þ
The value of H is the gyro kinetic moment, which is a constant obtained from the initial conditions of the motion, and J ðbÞ J is constant. In order to use the gyro as a control device, a feedback control system is added as shown in Fig. 2. The internal gimbal ring angle is sensed, and the output signal of the set performed by the potentiometer, operational amplifier and PID controller is applied to a DC motor. R, L, Kb and Km are the motor resistance, inductance, back-electromotive constant and the motor torque constant respectively. Considering the potentiometer and the circuit of the DC motor the following equations can be deduced: V2 ðtÞ ¼ Kpo bðtÞ; L
V1 ¼ Kpo br ;
Vs ðtÞ ¼ V2 ðtÞ V1 ðtÞ ¼ Kpo ½bðtÞ br ð9Þ
diðtÞ du þ RiðtÞ ¼ V ðtÞ Kb dt dt
where Kpo is the potentiometer constant, u is the giro angle of the shaft motor and V ðtÞ is the outlet voltage of the amplifier–PID controller. Taking into account that the gear train is considered as ideal, from the balance of the torque we get the following equations: MðtÞ ¼ Im
d2 u Ma ; dt2 j
MðtÞ ¼ Km iðtÞ;
d2 a Ma ¼ j Jm 2 jKm iðtÞ; dt 2
N2 j¼ N1
du da ¼ j dt dt
ð10Þ
where Im is the moment of inertia of motor shaft, Ma is the torque due to the feedback control and j is the tooth ratio of the gear train. Considering that equilibrium torque M acts in the OX2 , Eqs. (7) and (8) are the model of the gyro when: M þ Ma ¼ 0
ð11Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1061
In the general case, Eq. (11) is not verified, so the mathematical model of the gyro can be written as follows: J
d2 a db ¼ Ma þ M þ H cos b dt2 dt
d2 b da H cos b ¼ 0 2 dt dt dc da þ sin b ; H ¼C dt dt
ð12Þ
Hr
J ¼ A þ A1 þ A2
Eqs. (12) represent the gyroscope model with the torque Ma due to the feedback control system and the equilibrium torque M . In Eqs. (12) we assume that the platform of the whole system is static. However if the platform has an external oscillation hðtÞ, Eqs. (10) are not valid, because another torque appear in the On0 axis. In this case the relative velocity of the shaft motor is: dðu hÞ dða hÞ ¼ j dt dt
ð13Þ
where the h angle is measured respect to an horizontal static plane. Taking into account Eq. (13), the value of Ma can be rewritten as follows: Ma ¼ j2 Jm
d2 a d2 h þ jðj þ 1ÞJ jKm iðtÞ m dt2 dt2
ð14Þ
and consequently, from Eqs. (12) and (14) we get: 9 d2 a db d2 a d2 h 2 > > = ¼ j þ H cos b J þ jðj þ 1ÞJ þ M m m dt2 dt dt2 dt2 > d2 b da > ; Hr 2 H cos b ¼ 0 dt dt
J
ð15Þ
The amplifier–PID controller equation can be written as follows: Z 1 t dbðtÞ ðbðsÞ br Þ ds þ sd V ðtÞ ¼ Ka Kp bðtÞ br þ si 0 dt
ð16Þ
where Ka is the amplifier constant, Kp is the proportional constant (including the potentiometer constant Kpo ), si is the reset time, sd is the derivative time, and br is the reference angle of the internal gimbal ring. Note that if si ! 1 and sd ! 0, Eq. (13) represents the most simple strategy of control, i.e. control P [17,18]. From Eqs. (15) and (16) the equations of the system with PID feedback control and external vibration of the platform are the following: 9 d2 aðtÞ dbðtÞ d2 hðtÞ > > ¼ jK þ H cos b iðtÞ þ M þ jðj þ 1ÞJ > m m > dt2 dt dt2 > > > > 2 > > d bðtÞ daðtÞ > > ¼ 0 Hr H cos bðtÞ = dt2 dt Z t > 1 dbðtÞ > > > V ðtÞ ¼ Ka Kp bðtÞ br þ ðbðsÞ br Þ ds þ sd > > si 0 dt > > > > > diðtÞ daðtÞ dhðtÞ > ; þ RiðtÞ ¼ V ðtÞ þ jKb jKb L dt dt dt JT
ð17Þ
where JT ¼ J þ j2 Jm . The equations of the mathematical model show that the a variable is cyclic, so it is suitable to remove it by using Eq. (17). Introducing the state variables: x1 ðtÞ ¼ bðtÞ;
x2 ðtÞ ¼
dbðtÞ ; dt
x3 ðtÞ ¼ iðtÞ;
x4 ðtÞ ¼
diðtÞ ; dt
x5 ðtÞ ¼
Z o
t
ðbðsÞ br Þ ds
ð18Þ
1062
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
the following equations can be written:
9 x_ 1 ðtÞ ¼ x2 ðtÞ > > > > > H Ka Kp > > cos x1 ðtÞ Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ x5 ðtÞ x_ 2 ðtÞ ¼ > > > jKb Hr si > > > > H > _ > þ cos x1 ðtÞhðtÞ > > > Hr > > > > x_ 3 ðtÞ ¼ x4 ðtÞ = 2 j Km Kb R jKb H sd Ka Kp H > cos x1 ðtÞ x3 ðtÞ x4 ðtÞ x2 ðtÞ cos x1 ðtÞ þ x_ 4 ðtÞ ¼ > > L LJT jKb Hr L LJ > > T > > Ka Kp Ka Kp > > Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ x5 ðtÞ þ ½x1 ðtÞ x1r > > > si Lsi > > > > Ka Kp jKb jKb s HK K > d a p > € _ > x2 ðtÞ þ þ M þ ½jðj þ 1ÞIm JT hðtÞ þ cos x1 ðtÞhðtÞ > > L LJT LJT LHr > > ; x_ 5 ðtÞ ¼ x1 ðtÞ x1r
ð19Þ
It is interesting to point out that depending on the values of Kp , sd , si , it is possible to research different control strategies using the exact mathematical model of the gyro.
3. Analysis of the equilibrium points when the integral action is constrained If the integral action is constrained, Eqs. (19) can be simplified taking into account that it is assumed x5 ¼ Im , where Im is the maximum value that the PID controller integral action reaches. Assuming that the platform of the system is static, Eqs. (19) can be rewritten as follows: 9 x_ 1 ðtÞ ¼ x2 ðtÞ > > > > > Ka Kp > > Im x_ 2 ðtÞ ¼ ac cos x1 ðtÞ Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ > > si > > > > > x_ 3 ðtÞ ¼ x4 ðtÞ = ð20Þ x_ 4 ðtÞ ¼ bc x3 ðtÞ cc x2 ðtÞ cos x1 ðtÞ þ dc cos x1 ðtÞ > > > Ka Kp Ka Kp > > Im þ ½x1 ðtÞ x1r > Rx3 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ > > si Lsi > > > > Ka Kp jKb > > ; M x2 ðtÞ þ þ L LJT where ac bc cc and dc are parameters defined by the following equations: ac ¼
H ; jKb Hr
bc ¼
j2 Km Kb ; LJT
cc ¼
jKb H ; LJT
dc ¼
sd Ka Kp H sd Ka Kp ¼ ac jKb Hr L L
ð21Þ
It is interesting to point out that while the integral action is smaller than the maximum value Im , it is determined by Eq. (19), but when the value Im is reached, the change of the state variables is defined by Eq. (20). So, Eq. (19) can be considered as a transitory state for the variable x5 ðtÞ. Let us consider Eq. (20) when the integral action has reached the steady state value Im . It is clear that in the equilibrium state, the state variable x2 and x4 are always zero, consequently three different equilibrium points are possible: • Point P1 (x1e ; x2e ¼ 0, x3e ; x4e ¼ 0) In this case, considering fixed values of x1r , M , Kp and si , the value of x1e and x3e can be obtained from the following equations: 9 Ka Kp > Im ¼ 0 > = si Ka Kp jKb > ; ðx1e x1r Þ þ M ¼ 0> bc x3e þ Lsi LJT
Rx3e Ka Kp ðx1e x1r Þ
ð22Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
• Point P2 (x1e ¼ p=2, x2e ¼ 0, x3e ; x4e ¼ 0) Where the value of x3e can be found from the following equation: jK Ka Kp p b p=2 þ M ¼ 0 x1r þ bc x3e Lsi 2 LJT
1063
ð23Þ
• Point P3 (x1e ¼ 3p=2, x2e ¼ 0, x3e ; x4e ¼ 0) In this case, the value of x3e can be found from the following equation: Ka Kp 3p jKb 3p=2 M ¼0 bc x3e þ x1r þ Rsi RJT 2
ð24Þ
Note that for n integer, and x1e ¼ ð4n þ 1Þp=2, x1e ¼ ð4n þ 1Þ3p=2 we also obtain the equilibrium points P2 and P3 respectively. On the other hand, when the integral action is constrained to an arbitrary fixed value Im the following inequalities at points P2 and P3 are verified: 9 p KK a p p=2 > x1r Ka Kp Im 6¼ 0 > P2 : Rx3e = 2 si ð25Þ 3p Ka Kp > > x P3 : Rx3p=2 K K I ¼ 6 0 ; a p 1r m 3e 2 si Nevertheless, it is possible to find the value of Im at points P2 and P3 . Taking into account Eqs. (22), at point P2 (p=2) the following equations are verified: 9 p KK a p p=2 > > x Rxp=2 K K I ¼ 0 = a p 1r 3e 2 si ð26Þ K K p jK > a p b ; x1r þ bc xp=2 M ¼ 0 > 3e þ Lsi 2 LJT Eliminating xp=2 3e by using Eqs. (26) one obtains: si RM RJT p x1r si 2 I p=2 ¼ 2 Ka Kp jKm j Km Kb
ð27Þ
The same reasoning allows us to obtain the integral action which belongs to equilibrium point P3 : si RM RJT 3p x1r si 2 I 3p=2 ¼ 2 Ka Kp jKm j Km Kb
ð28Þ
Note that from Eqs. (27) and (28) t is deduced that if x1r < p=2 then I 3p=2 < I p=2 . Once the equilibrium points are determined, by considering the parameter values of the PID controller and from certain initial conditions, the system will evolve to point P1 if the eigenvalues of the Jacobian matrix of the linearized system have negative real part, and points P2 and P3 are unreachable. Consequently it is necessary to research the behaviour of the system at points P1 , P2 and P3 . The analysis of the equilibrium points is carried out transforming the fixed points to the origin by the translation: x01 ðtÞ ¼ x1 ðtÞ x1e ;
x02 ðtÞ ¼ x2 ðtÞ x2e ;
x03 ðtÞ ¼ x3 ðtÞ x3e ;
x04 ðtÞ ¼ x4 ðtÞ x4e
ð29Þ
3.1. Analysis of point P1 Substituting Eq. (29) into Eq. (20) and taking into account the equilibrium equations (22) we get the following equations: 9 x_ 01 ðtÞ ¼ x02 ðtÞ > 0 > > 0 0 0 0 0 > x_ 2 ðtÞ ¼ ac cosðx1 ðtÞ þ x1e Þ Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp x1 ðtÞ sd Ka Kp x2 ðtÞ > > > 0 0 > x_ 3 ðtÞ ¼ x4 ðtÞ = ð30Þ c x_ 04 ðtÞ ¼ dc cosðx01 ðtÞ þ x1e Þ Rx03 ðtÞ þ Lx04 ðtÞ Ka Kp x01 ðtÞ sd Ka Kp þ c x02 ðtÞ > > > dc > > > > Ka Kp 0 Ka Kp 0 R > ; x2 ðtÞ bc x03 ðtÞ x04 ðtÞ þ x1 ðtÞ þ Lsi L L
1064
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
The equilibrium point of Eq. (30) is the origin, and the eigenvalues of the Jacobian matrix are defined by the following equations: k4 þ b3 k3 þ b2 k2 þ b1 k þ b0 ¼ 0 R b2 ¼ bc þ Lac cc cos2 x1e ; b3 ¼ ; L R cos x1e b0 ¼ ac Ka Kp bc Lsi
b1 ¼
Rac cc cos x1e þ sd Ka Kp ac cc ac
Ka Kp si
cos x1e
ð31Þ
From Routh–Hurwitz criterion [17] it is deduced that point P1 will be stable if the roots of the characteristic equation have negative real part, i.e. when p=2 < x1e < p=2, b2 b3 > b1 and b1 ðb2 b3 b1 Þ > b0 b23 . 3.2. Analysis of point P2 Considering Eq. (30), substituting x1e ¼ p=2, and taking Im ¼ I p=2 , Im > 0 in Eq. (27), the Jacobian matrix of the system (30) has two zero eigenvalues, so a Takens–Bogdanov bifurcation can appear. Assume that the integral action is constrained to Im ¼ I p=2 e with e > 0, Im > 0, then the system (30) has a pure imaginary pair of eigenvalues and so point P2 will be a weak focus. In contrast to a stable focus in a weak focus the convergence of the trajectory is not exponential, i.e. the trajectory is not a logarithmic spiral and its length tends to infinity as the radius of the trajectory tends to zero [19]. However, if the integral action is constrained to Im ¼ I p=2 þ e with e > 0, Im > 0 it is easy to show that the equilibrium point P2 is a saddle, and consequently, this point will be unreachable. From the previous reasoning, it is deduced that while the integral action (x5 ðtÞ) is smaller than I p=2 the system will evolve to point P2 if this point is a stable weak focus. On the other hand, if P2 is an unstable weak focus the system will jump to another equilibrium point, P1 or P3 . Finally, when the integral action reaches a value greater than I p=2 the point P2 will be unreachable and the system will change to point P1 or P3 . Consequently point P2 will be reachable when the following conditions are satisfied: • The integral action is constrained to Im ¼ I p=2 e with e > 0, Im > 0 and a weak focus with a pure imaginary pair of eigenvalues appear. • The weak focus is stable. 3.3. Analysis of point P3 Assuming that the integral action is constrained to Im ¼ I p=2 þ e, with e > 0, Im > 0 and substituting x1e by 3p=2 in Eq. (30), the equations of the system are the following: 9 x_ 01 ðtÞ ¼ x02 ðtÞ > > > > x_ 02 ðtÞ ¼ ac sin x01 ðtÞ½Rx03 ðtÞ þ Lx04 ðtÞ Ka Kp x01 ðtÞ sd Ka Kp x02 ðtÞ þ p > > > > x_ 03 ðtÞ ¼ x04 ðtÞ = c ð32Þ x_ 04 ðtÞ ¼ dc sin x01 ðtÞ Rx03 ðtÞ þ Lx04 ðtÞ Ka Kp x01 ðtÞ sd Ka Kp þ c x02 ðtÞ þ p > > > dc > > > Ka Kp 0 Ka Kp 0 R > > ; x2 ðtÞ bc x03 ðtÞ x04 ðtÞ þ x1 ðtÞ þ L Lsi L where 3p=2
p ¼ Rx3e
Ka Kp
3p x1r 2
Ka Kp p=2 ðI þ eÞ si
ð33Þ
and I p=2 between (24), (27) and (33), the following equation is deduced: Eliminating x3p=2 3e JT R Ka Kp p ¼ Ka Kp p 1 2 e ð34Þ si j Km Kb si pffiffiffiffiffiffiffiffiffiffi It is easy to show that if p < 0, Eq. (32) have eigenvalues k1;2 ¼ i ac jpj and consequently it is necessary to research the equilibrium state of point P3 carefully. 4. Analysis of the stability of point P3 From the conclusions of the previous section it can be shown that the equilibrium point P3 is a weak focus. The stability of a weak focus for systems of high dimension can be difficult to determine [12,13,19–21]. In the case of point
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1065
P3 the system has four dimensions, so it will be necessary to carry out a set of previous steps in order to put Eq. (32) in an appropriate form, to determine the first Lyapunov value [19,20]. These steps are summarized as follows: Step 1 Make a Taylor expand of the function sin x01 ðtÞ up to order 3. Substituting in Eq. (32), we next split off the linear and nonlinear part of the vector field to get 2 3 2 32 0 3 2 3 0 0 x1 0 1 0 0 6 x_ 10 7 6 0 7 6 0 0 0 0 7 6 x_ 2 7 6 c 0 0 0 7 76 x20 7 þ 6 f2 ðx1 ; x2 ; x3 ; x4 Þ 7 þ Oðkx0 k4 Þ 6 0 7¼6 ð35Þ 54 x 5 4 5 6 x_ 7 4 0 0 0 0 1 3 4 30 5 0 0 0 0 0 f4 ðx1 ; x2 ; x3 ; x4 Þ x_ 4 x4 a b bc n where a ¼ dc p þ
c ¼ ac p;
Ka Kp ; Lsi
b¼
Ka Kp ; L
n¼
R L
ð36Þ
9 = f2 ¼ ac x01 ðRx03 þ Lx04 Ka Kp x01 sd Ka Kp x02 Þ a6c p x033 h i cc dc p 03 0 0 0 0 0 f4 ¼ dc x1 Rx3 þ Lx4 Ka Kp x1 sd Ka Kp þ dc x2 6 x3 ;
ð37Þ
The eigenvalues of linear part of (35) are: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi k2;3 ¼ n=2 i bc ðn=2Þ2 k1;2 ¼ ix ¼ i ac jpj;
ð38Þ
and consequently an Andronov–Poincare–Hopf bifurcation appears [21,22]. Step 2 Let P be the matrix that transforms transformation: 2 3 2 x01 0 1 6 07 6 6x 7 6 x 0 6 27 x0 ¼ P x ) 6 0 7 ¼ 6 6 x3 7 6 4 5 4 g f x04
f ¼
aðbc x2 Þ þ bnx2 ðbc
x 2 Þ2
Eq. (35) becomes 2 3 2 x_ 0 6 y_ 7 6 x 6 7¼6 4 z_ 5 4 0 w_ a
2
þn
x 0 0 b
x2
0 0 a1 b1
xf
xg
;
g¼
the linear part of (35) into Jordan canonical form. Then, under the linear 0 0 0 b1
32
3 x 76 7 0 76 y 7 76 7; 6 7 1 7 54 z 5 0
a1
R a1 ¼ ; 2L
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 R b1 ¼ bc 2L
ð39Þ
w
x½bðbc x2 Þ an ðbc x2 Þ2 þ n2 x2
2 3 32 3 0 0 x 6 7 6 7 0 7 76 y 7 þ P 1 6 f2 ðx; y; z; wÞ 7 þ Oðkxk4 Þ 4 5 b1 54 z 5 0 a1 f4 ðx; y; z; wÞ w
ð40Þ
where f2 and f3 are the same functions of Eq. (37) but substituting coordinates x0 by x. Step 3 We now simplify Eq. (35) as much as possible using the center manifold theorem [12,13,19,20] in order to reduce the dimensionality to two in a neighbourhood of the equilibrium point P3 . In the center manifold the variables z and w are approximated by: 2 z h ðx; yÞ a x þ a2 y 2 þ a3 xy ¼ 1 2 ð41Þ ¼ 1 h2 ðx; yÞ w b1 x þ b2 y 2 þ b3 xy Coefficients ai , bi (i ¼ 1; 2; 3) are determined by substituting Eq. (41) into the equation which defines the center manifold [12,13,20]: x þ f ðx; y; h1 ðx; yÞ; h2 ðx; yÞÞ ¼ Bhðx; yÞ þ gðx; y; h1 ðx; yÞ; h2 ðx; yÞÞ Dx hðx; yÞ½A
ð42Þ
1066
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
where " # x x ¼ ; y " B¼
" ¼ A
a1
b1
b1
a1
# ;
" # 1 f2 ðx; y; hðx; yÞÞ ; f ðx; y; hðx; yÞÞ ¼ x x 0 0 2 3 ðga1 xf Þf2 f4 þ 7 6 oh oh b1 x b1 7 6 g¼6 Dx h ¼ ; 7; 4 5 ox oy gf 0
x
#
ð43Þ
2
x We then equate powers of x2 , xy, y 2 in the left and right sides of Eq. (42), so coefficients ai , bi can be computed. The simplified system at the center manifold in normal form can be written as follows: 1 f2 ðx; yÞ 0 x x x_ ¼ þ ð44Þ y_ x 0 y 0 x where f2 ðx; yÞ ¼ a02 y 2 þ a11 xy þ a03 y 3 þ ac ðR þ La1 Þyðb1 x2 þ b2 y 2 þ b3 xyÞ þ ac Lb1 yða1 x2 þ a2 y 2 þ a3 xyÞ ac p a02 ¼ ac ðRf Lxf Ka Kp Þ; a11 ¼ ac ðRg þ Lxf sd Ka Kp xÞ; a03 ¼ 6 Step 4 Considering a system of the form: 1 0 x x x_ f ðx; yÞ þ 2 ¼ y x 0 y_ f ðx; yÞ
ð45Þ
ð46Þ
there is an explicit equation that gives the first Lyapunov value L1 [12,13,19–22]: i 1 1 1 h 1 1 1 2 2 fxy ðfxx þ fyy1 Þ fxy2 ðfxx2 þ fyy2 Þ fxx1 fxx2 þ fyy1 fyy2 L1 ¼ ½fxxx þ fxyy þ fxxy þ fyyy þ 16 16x
ð47Þ
Applying Eq. (47) to Eqs. (45), (46), and considering that all partial derivatives are evaluated at the bifurcation point P3 (x ¼ 0, y ¼ 0), the value of L1 is: L1 ¼
a11 a02 ac ðR þ La1 Þb3 ac Lb1 a3 þ þ 8x3 8x 8x
ð48Þ
so, only values of a3 and b3 are necessary to know the first Lyapunov value. The value of parameters a3 and b3 obtained from the center manifold theorem (Step 3) are the following: g a3 ¼ F1 ðp1 c11 Þ þ F2 p2 a11 x g b3 ¼ F2 ðp1 c11 Þ þ F1 p2 a11 x
ð49Þ
where F1 ¼
a1 þ a1 f1 2
ða1 þ a1 f1 Þ þ ðb1 b1 f1 Þ
2
;
F2 ¼
b1 b1 f1 2
ða1 þ a1 f1 Þ þ ðb1 b1 f1 Þ
9 ðga1 xf Þa02 b02 ga02 > > þ b a þ p1 ¼ 2x > 1 1 bc > b1 x b1 x > > > > > > ðga xf Þa b ga 1 02 02 02 > > þ þ a b p2 ¼ 2x 1 > 1 bc > b1 x b1 x > = b02 ¼ dc ðRf Lxg Ka Kp Þ > > > ðga1 xf Þa11 b11 > > þ c11 ¼ > > > b1 x b1 > > > > > a11 ¼ ac ðRg þ Lxf sd Ka Kp Þ > > ; b11 ¼ dc ðRg þ Lxf sd Ka Kp Þ cc x
2
;
f1 ¼
4x2 þ b21
a21
ð50Þ
ð51Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1067
If L1 > 0 the equilibrium point P3 will be unstable. When L1 < 0 the equilibrium point will be stable, and if L1 ¼ 0 it is necessary to reduce the system (35) on the center manifold to the formal norm up to the terms of fifth order [19–22]. Consequently Eqs. (45), (46) and (48)–(51) give us a criterion to know the stability proprieties of point P3 . 5. Regular and self-oscillating behaviour The analysis of the equilibrium points reveals that when the integral action is constrained to a value given by Im ¼p=2 þe the equilibrium point P2 will be unreachable and if the first Lyapunov value defined by Eqs. (48)–(51) is positive, the equilibrium point P3 will be an unstable weak focus. Consequently, the system changes to equilibrium point P1 if all roots of the characteristic equation (31) have negative real part and there is no external vibration of the platform. Therefore, if one chooses adequate values for the parameter Kp , si , sd of the PID controller it is possible to reach a desired set point x1e . In this case we say that the behaviour of the gyro is regular. Nevertheless, it is possible to obtain other movements when the gyro changes to equilibrium points P2 , P3 or reaches a self-oscillating state. Note that it is possible to choose different combinations of the parameters of the PID controller in order to obtain a desired set point. In this section several examples are carried out to examine different behaviours of the gyro by numerical simulation. The Runge–Kutta–Felhberg method has been used in numerical simulations with integration steps between 0.001 and 0.0005 s. The numerical values and the parameters of the gyro are shown in Table 1. We first consider the dynamics of the gyro when the set point x1e and the reference x1r are different. The simulation results are shown in Figs. 3 and 4. The parameters of the PID controller Kp , si and sd are chosen from the standard procedures [17,18,23] taking into account that for fixed values of si , sd the proportional constant Kp must give a positive first Lyapunov value, so the equilibrium point P3 (3p=2) will be an unstable weak focus as shown in Fig. 3. The value e ¼ 0:5 > 0 gives an equilibrium point P1 of 1.297 rad and the integral action is constrained to 8.79 s, as shown in Fig. 4. Note that the eigenvalues of point P1 have negative real part, and consequently this point will be reachable. The moment of inertia of the motor shaft must be Jm . Description and values are correct. The table is the following: Table 1 Gyroscope and control system parameter values Variable
Description
Value
A1 B1 C1 A2 A C Hr J Jm j JT nr H Km Kb R L Ka Kp si sd I Im ac bc cc dc M
Internal gimbal moment of inertia respect to OX1 (kg m2 ) Internal gimbal moment of inertia respect to OY1 (kg m2 ) Internal gimbal moment of inertia respect to OZ1 (kg m2 ) External gimbal moment of inertia respect to OX2 (kg m2 ) Rotor equator moment of inertia (kg m2 ) Rotor polar moment of inertia (kg m2 ) A þ B1 (kg m2 ) A þ A1 þ A2 (kg m2 ) Moment of inertia of the motor shaft (kg m2 ) N2 =N1 Tooh ratio J þ j2 Jm (kg m2 ) Rotor angular velocity (rad/s) C nr Gyrokinetic moment (kg m2 rad/s) Motor torque constant (N m/A) Back-electromotive constant (V s) Motor resistance (X) Motor inductance (H) Amplifier constant Proportional constant of the PID controller (V) Reset time of the PID controller (s) Derivative time of the PID controller (s) Integral action (s) Maximum value of the integral action (s) Parameter H =jKb Hr (rad/V s) Parameter j2 Km Kb =LJT (s1 ) Parameter jKb =LJT (A) Parameter sd Ka Kp =jKb Hr L (A V/s) Disturbance torque (N m)
0.045 0.045 0.05625 0.04875 0.01125 0.0225 0.05625 0.105 0.001 10 0.205 104.72–523.59 2.3562–11.781 5 0.04 0.5–1 0.01 10 0.01–0.05 1–2 0.2–0.5
9756.097 195.122 10–200
1068
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 3. Simulation results for transient and set point steady state response from the following values: nr ¼ 1500 r.p.m.; M ¼ 200 N m; x1r ¼ p=8 rad; e ¼ 0:5; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; Im ¼ 8:7937. Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0, x5 ð0Þ ¼ 0. Eigenvalues at the set point: [)95.8 ± 48.5j, )3.656 ± 1.357j]. First Lyapunov value at P3 : L1 ¼ 0:02569 > 0.
Fig. 4. Steady state of the PID controller integral action. Simulation values from Fig. 3.
Similar results are obtained when the set point or equilibrium point P1 coincides with the reference angle br ¼ x1r . It is important to remark that if the value of the reference angle br ¼ x1r is greater than a certain value, from the last equation of (19), it is possible that the integral action reaches a negative value. Consequently, the steady state of the integral action will be Im . In this case, applying the same procedure as before, the equilibrium point P2 is a stable weak focus and the gyro has not regular behaviour, i.e. is not possible to reach the set point. Fig. 5 shows this case. Fig. 6 shows a case in which the first Lyapunov value is negative, so point P3 is a stable weak focus. In this case it is not possible to reach the equilibrium point. Note that while the integral action x5 ðtÞ is lower than Im point P2 is a stable weak focus, but when the limit of the integral action is reached, the gyro jumps to point P3 . The self-oscillating behaviour is researched from the previous consideration of Section 3. From Eqs. (22) and (26), if the integral action is constrained to I p=2 þ e with a small positive value of e, then the equilibrium point P1 will be almost coincident with the equilibrium point P2 . The basins of attraction of P1 and P2 collide and by considering the discussion of Section 3, these points will be unreachable. So, the only possibility of the gyro is to evolve to point P3 . However, point P3 is an unreachable weak focus if the first Lyapunov value is positive. So, if L1 > 0 a self-oscillatory behaviour
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1069
Fig. 5. Simulation results for transient and steady state response when the PID controller integral action is constrained to a negative value. Point P2 is a stable weak focus. Simulation values: nr ¼ 1500 r.p.m.; M ¼ 100 N m; x1r ¼ 2 rad; e ¼ 105 ; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; I p=2 ¼ 7:8088; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0, x5 ð0Þ ¼ 0.
Fig. 6. Simulation results when the first Lyapunov at point P3 is negative L1 ¼ 0:5 and P1 P2 . Simulation values: nr ¼ 2000 r.p.m.; M ¼ 100 N m; x1r ¼ p=8; e ¼ 105 ; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; I p=2 ¼ 4:88; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0; x5 ð0Þ ¼ 0. Eigenvalues at the set point P1 x1e p=2 [)25 ± 95.56j; )6.73 · 105 ± 172.4j].
appears. This behaviour has been corroborated by simulation as shown in Fig. 7. Note that if the integral action is constrained to an adequate value it will always possible to chose parameters of the PID controller in order to obtain L1 > 0.
6. Chaotic behaviour Nowadays, the theory of chaotic oscillations is considered as an interesting field of research. The qualitative features of many models of oscillation have been previously applied to systematic experimental and theoretical developments of mechanical oscillators, both in numerical and experimental studies [24–27]. In dissipative electromechanical systems, the presence of a strange attractor, which appears when the system has chaotic behaviour, is of paramount importance in order to understand new nonlinear proprieties of the electromechanical device. Nowadays, in order to carry out theoretical and experimental studies of a strange attractor, the following steps are recognized as necessary [28,29]:
1070
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 7. Self-oscillating behaviour when the first Lyapunov value at point P3 is positive L1 ¼ 0:01 and P1 P2 . Simulation values: nr ¼ 1000 r.p.m.; M ¼ 100 N m; x1r ¼ p=8; e ¼ 105 ; Kp ¼ 0:02 V; si ¼ 2 s; sd ¼ 0:5 s; I p=2 ¼ 7:65; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0; x5 ð0Þ ¼ 0. Eigenvalues at the set point P1 x1e p=2 [)25 ± 95.56j; )2.59 · 105 ± 1.02 · 102 j].
ii(i) The development of mathematical models based on plausible subsystems of the electromechanical system. i(ii) The analysis and design of the whole system in order to understand strange behaviours. (iii) The possibility of simulating the mathematical model with a clear criterion to elucidate between regular, transient chaotic and chaotic behaviour. With these ideas the study of chaotic behaviour of the gyro is carried out taking into account the model given by Eq. (19). When the platform of the gyro has an external vibration movement, the relationship of the platform external disturbance can be written as follows: _ ¼ Ap xp cos xp t; hðtÞ ¼ Ap sin xp t ) hðtÞ
€hðtÞ ¼ Ap x2 sin xp t p
ð52Þ
The values of the amplitude Ap and the angular frequency xp must be properly chosen in order to obtain chaotic behaviour. Previously, the system is simulated with Ap ¼ 0, xp ¼ 0 and the integral action constrained in order to obtain a first Lyapunov value L1 > 0. In this case, the behaviour of the system will be self-oscillating. Fig. 8 shows the value of L1 as
Fig. 8. Variation of the first Lyapunov value vs. Kp and self-oscillating behaviour with L1 > 0 and P1 P2 . Simulation values: nr ¼ 1200 r.p.m.; M ¼ 100 N m; x1r ¼ p=8; e ¼ 105 ; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; I p=2 ¼ 4:88; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0; x5 ð0Þ ¼ 0. Eigenvalues at the set point P1 x1e p=2 [)25 ± 95.56j; )4.04 · 105 ± 1.335 · 105 j].
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1071
Fig. 9. Strange attractor when the platform has an external harmonic disturbance. Simulation values of Fig. 8. Amplitude Av ¼ 0:065 m, frequency xv ¼ 8:6 rad/s.
function of the proportional constant Kp of the PID controller. Note that for Kp ¼ 0:026 the value of L1 is positive, and the point P3 is an unstable weak focus. The period of the self-oscillation is close to 2.65 s. When the system has self-oscillating behaviour, an external disturbance in the platform of the gyro as shown in Eq. _ (52), can produce chaotic behaviour. Fig. 9 shows the strange attractor drawn in the bðtÞ–bðtÞ plane. The zones with strong stretching and folding as well as zones with Cantor’s sets are shown in Fig. 9. On the other hand, it is well known that it is very difficult to show rigorously that a strange attractor is chaotic. It can be said that a nonlinear system has chaotic behaviour if (a) it has sensitive dependence on initial conditions, i.e. two very close initial conditions diverge exponentially with time, (b) the orbits are dense in a state region i.e. the orbits fill the phase space zone of the strange attractor, and (c) the orbits are topologically transitive in i.e. for any two open sets X1 ; X2 X there is a time from which any orbit starting at X1 ends at X2 [13]. While the conditions (a) and (b) can be approximately verified by simulation, proving the condition (c) is generally very difficult. Nevertheless, the fulfilment of (a) and (b) can be enough to assure the long time chaotic behaviour. Fig. 10 shows that the system has sensitive dependence to initial conditions, and Fig. 11 shows Lyapunov’s exponents of the system. The algorithm that has been used can be found in [26,27,30,31]. Eqs. (19) are transformed in an autonomous dynamical system considering a new variable x6 ðtÞ defined by: x6 ðtÞ ¼ xp t )
dx6 ðtÞ ¼ xp dt
Fig. 10. Sensitive dependence and error vs. time with two initial conditions very close. Simulation values of Fig. 8.
ð53Þ
1072
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 11. Lyapunov’s exponents for the chaotic strange attractor of Fig. 9.
Note the presence of a positive Lyapunov exponent, which is an indicator of chaos, and a Lyapunov exponent nil, corresponding to x6 ðtÞ periodic variable. In order to verify the computations, the divergence of vector field from Eq. (19) is calculated as follows: div~ f ¼ sd Ka Kp
H R sd Ka Kp H cos x1 ðtÞ cos x1 ðtÞ þ jKb Hr L jKb H
ð54Þ
then, the sum of Lyapunov exponents are determined, and the following relation: div~ f ¼
6 X
ð55Þ
NLi
i¼1
is verified by simulation. The results obtained with an integration step of 0.004 s are the following: div~ f ¼ 50:13;
6 X
NLi ¼ 50:01
ð56Þ
i¼1
therefore the computations can be considered correct. Fig. 12 shows another indicator of chaos, i.e. the Fourier power spectrum for the variable ðtÞ ¼ x2 ðtÞ shown at Fig. 10. The power spectrum is a representation of the relative abundance of different frequencies in a given time series. Note that the power spectrum is broadband and contains substantial at low frequencies.
Fig. 12. Power spectrum for the bðtÞ variable of Fig. 9.
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1073
Finally, it is important to point out that although both a positive Lyapunov exponent and broad spectrum do not guarantee sensitivity to initial conditions it is, in practice, an indicator of chaos. The straight line AB has been drawn through the harmonics peaks to remark the character of chaotic dynamic, probably due to a route of period-doubling sequences [12,13,26–29]. 7. Conclusions From the results of the research presented in this paper, it has been corroborated that a gyroscope in gimbal with a feedback PID controller, can be operated by using a strategy that allows to reach a desired set point, self-oscillation dynamics as well as a more complex behaviour, characterized by chaotic dynamics. By using the exact nutation theory of the gyroscope and the PID controller with constrained integral action, we show that the equations of the gyroscope and the feedback control system have multiple equilibrium points. The analysis of these equilibrium points reveals that depending on a fixed value of the PID controller integral action three different equilibrium points may appear. One of them, point P1 , with positive or negative real part of eigenvalues, another point P2 (p=2), with a pure imaginary pair or two zero eigenvalues, and a third point P3 (3p=2), with a pure imaginary pair of eigenvalues. For a certain integral action value, we show that point P2 is unreachable, and point P3 is a weak focus. Then, the stability of point P3 , by using the manifold center theorem, the normal form theory and the first Lyapunov value is researched. The sign of the first Lyapunov value gives us a clear criterion in order to elucidate if point P3 is a stable or unstable weak focus. The determination of normal form reveals that the algebraic machinery can be formidable, nevertheless it follows straightforward rules. Advances in both digital and symbolic computational methods can be investigated in connection with the analysis of complex system with several equilibrium points, as it is shown in this paper. When points P1 and P2 are unreachable, if the first Lyapunov value at point P3 is negative, this point will be a stable weak focus, as corresponds to an Andronov–Poincare–Hopf bifurcation with a pure imaginary pair of eigenvalues. On the other hand, if the first Lyapunov value is positive, point P3 is unreachable too and the gyro has a self-oscillating motion. The self-oscillating behaviour allows us to find chaotic oscillations when the platform is submitted to an external harmonic disturbance so, a new family of strange attractors has been found. The chaotic behaviour has been researched from the sensitive dependence to initial conditions, Lyapunov’s exponents and the power spectral density, from which it is possible to verify that the mechanism of chaotic oscillations is due to a period-doubling phenomenon. Another interesting aspect of this paper is the possibility of designing an algorithm, which gives the PID parameters values i.e. the proportional constant Kp , the reset time si , and the derivative time sd , in order to obtain a previous desired set point. This algorithm will be based in the eigenvalues of point P1 with negative real part, point P2 with a saddle and point P3 with a positive first Lyapunov value. The previously analyzed phenomena should be known by the mechanical designer in order to either avoid or properly use them, depending on gyro type. The paper shows a research line, where more elaborated models of gimbal suspension gyro with structural flexibility, and more sophisticated control strategies can be used to plan both theoretical and experimental studies. References € Klein F, Sommerfield A. Uber die theorie des kreisels. H. 1–4. Leipzig-Berlin: Teubner, 1897/1898–1903/1910. Timoshenko S, Young DH. Advanced dynamics. New York: McGraw-Hill; 1948. p. 392–448. Pars LA. A treatise on analytical dynamics. New York: John Wiley & Sons; 1966. p. 159–87. Cabannes H. Cours de Mecanique Generale. Paris: Dunod; 1966. p. 248–69. Goldstein H. Classical mechanics. New York: Addison Wesley; 1980. Marsden JE, Ratiu SI. Introduction to mechanics and symmetry. New York: Springer Verlag; 1999. p. 481–512. Chetaev NG. Theoretical mechanics. Moscow: Mir; 1989. p. 221–5. Arnold VI. Mathematical methods of classical mechanics. New York: Springer Verlag; 1974. p. 75–90, 157–95. Cannon RH. Dynamics of physical systems. New York: McGraw-Hill; 1967. p. 159–63, 617–26. Singh SN. Stability of gyro with harmonic nonlinearity in spinning vehicle. IEEE Trans Aerospace Electron Syst 1983;19:182–9. Ge ZM, Chen HH. Bifurcations and chaos in a rate gyro with harmonic excitation. J Sound Vib 1996;194(1):107–17. Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems and bifurcations of vector fields. New York: Springer Verlag; 1983. p. 117–80. [13] Wiggins S. Introduction to applied nonlinear dynamical systems and chaos. New York: Springer Verlag; 1990. p. 193–202, 211–29, 253–84. [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
1074
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
[14] Ge ZM, Chen HH. Double degeneration and chaos in a rate gyro with feedback control. J Sound Vib 1998;209(5):753–69. [15] Agafonov SA. Gimbal suspension gyro: stability, bifurcation, and chaos. J Dyn Control Syst 2001;7(3):339–51. [16] Kothare M, Campo P, Morari M, Nett C. A unified framework for the study of anti wind-up designs. Automatica 1994;30(12):1869–83. [17] Franklin GF, Powell JD, Emani-Naeni A. Feedback control of dynamics systems. New York: Prentice-Hall; 2000. p. 95–100. [18] Goodwin GC, Graebe SF, Salgado ME. Control system design. Saddle River, New Jersey: Prentice-Hall; 2001. p. 159–76. [19] Shilnikov LP, Shilnikov AL, Turaev DV, Chua LO. Methods of qualitative theory in nonlinear dynamics. Part II. New Jersey: World Scientific; 2001. p. 451–65, 598–610. [20] Kuznetsov YA. Elements of applied bifurcation theory. New York: Springer Verlag; 1998. p. 86–109, 151–71. [21] Glendinning P. Stability, instability and chaos: an introduction to the theory of nonlinear differential equations. New York: Cambridge University Press; 1996. p. 224–43. [22] Marsden JE, McCraken M. The Hopf bifurcation and its applications. New York: Springer Verlag; 1976. [23] Ogata K. Modern engineering control. New York: Prentice-Hall; 2000. [24] Rudiger S. Practical bifurcation and stability analysis. From equilibrium to chaos. New York: Springer Verlag; 1994. p. 327–57. [25] Wggins S. Global bifurcations and chaos. Analytical methods. New York: Springer Verlag; 1988. p. 75–108. [26] Lichtenberg AJ, Lieberman MA. Regular and chaotic dynamics. New York: Springer Verlag; 1992. p. 62–5, 312–20, 457–78. [27] Holden AV, editor. Chaos. Princeton, New Jersey: Princeton University Press; 1986. [28] Hilborn RC. Chaos and nonlinear dynamics. Oxford: Oxford University Press; 2000. p. 121–48, 323–41. [29] Ott E. Chaos in dynamical systems. Cambridge: Cambridge University Press; 2002. p. 115–37, 212–36, 305–10. [30] Benettin G, Galgani L, Giorgilly A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them, Part I: Theory. Meccanica 1980;15:9–20. [31] Benettin G, Galgani L, Giorgilly A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them, Part II: Numerical applications. Meccanica 1980;15:21–30.
Regular self-oscillating and chaotic behaviour of a PID controlled gimbal suspension gyro Manuel F. Perez Polo
a,*
, Manuel Perez Molina
b,1
a
b
Department of Fısica, Ingenierıa de Sistemas y Teorıa de la Se~nal, Universidad de Alicante, Escuela Politecnica Superior, Campus de San Vicente, Alicante 03071, Spain Facultad de Ciencias Matematicas, Universidad Nacional de Educacion a Distancia, UNED, C/Boyero 12-1A, Alicante 03007, Spain Accepted 6 January 2004
Abstract The dynamics of a gyro in gimbal with a PID controller to obtain steady state, self-oscillating and chaotic motion is considered in this paper. The mathematical model of the whole system is deduced from the gyroscope nutation theory and from a feedback control system formed by a PID controller with constrained integral action. The paper shows that the gyro and the associated PID feedback control system have multiple equilibrium points, and from the analysis of a Poincare–Andronov–Hopf bifurcation at the equilibrium points, it is possible to deduce the conditions, which give regular and self-oscillating behaviour. The calculation of the first Lyapunov value is used to predict the motion of the gyro in order to obtain a desired equilibrium point or self-oscillating behaviour. The mechanism of the stability loss of the gyro under small vibrations of the gyro platform and the appearance of chaotic motion is also presented. Numerical simulations are performed to verify the analytical results. 2004 Elsevier Ltd. All rights reserved.
1. Introduction Throughout the years, the dynamics of gyroscopic devices has been considered as one of the most exciting problems of mechanics. From the classical handbook by Klein and Sommerfield [1] at the beginning of the last century until the present time, many books and papers concerning the stability, nonlinear phenomena and technical applications of gyroscopic systems have been be written [2–6]. Gyros for sensing angular motion are used in airplane automatic pilots, rocket-vehicle launch-guidance, space-vehicle attitude systems, ship’s gyrocompasses and submarine inertial autonavigators. Today, fibber optic gyro based has many applications in avionics, platform stabilization and automotive testing. In balanced gimbal suspension gyro systems, it is well known that the exact equations which describe the nutation motion are nonlinear [7,8]. However, for purpose of simplified analysis, a linear model has been deduced in order to obtain a transfer function that can be used in the design of a controller. Exact equations of motion from an ideal twoaxis gyro, the linear approximation for small motions, and the response to a constant torque to give precession and nutation motion can be found in [9]. Singh [10], Ge and Chen [11], taking into account harmonic external disturbances have already analyzed the nonlinear dynamics of a single-axis gyro mounted on a space-vehicle. These papers show that from averaging method [12,13] and the harmonic balance it is possible to obtain chaotic behaviour. A similar system, but with feedback
*
Corresponding author. Tel.: +34-9-659-09673; fax: +34-9-65-909750. E-mail addresses: manolo@dfists.ua.es (M.F. Perez Polo), [email protected] (M.P. Molina). 1 Tel.: +34-9-651-01988; fax: +34-9-659-09750. 0960-0779/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.chaos.2004.01.011
1058
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
proportional control was considered in [14], using normal forms and Melnikov method [12,13]. This paper shows that regular and chaotic motions are both possible. In the paper of Agafonov [15], the stability of stationary motion, the periodic motions from a Poincare–Andronov–Hopf bifurcation and the chaotic behaviour of a gimbal suspension gyro is researched. However, in some previously cited papers, the gyro mathematical model is very simplified without taking into account the control problem, and in others papers a very simple proportional control is only used to obtain a desired set point. So, the whole effect of a PID linear controller dynamics and the exact nonlinear mathematical model of the gyro (nutation motion) are completely unknown. This paper considers the motion of a heavy symmetrical gyro in gimbal. The mathematical model is deduced from Euler–Lagrange equations, considering that the gyro moment of inertia is constant. A PID feedback control system is added, taking the angle of the internal gimbal ring as feedback signal to produce the error signal equal to set point minus feedback signal. The error signal is the input to an amplifier–PID controller to produce the output control signal, which is sent to a DC motor. The shaft of the DC motor is connected by means of an ideal gear train to the external gimbal. The external gimbal ring is also connected to an external body, which is subject to a constant external disturbance torque. The objective of the control system is that the external gimbal reaches a desired set point by considering the external disturbance torque. The whole system is mounted on a platform which is subject to periodic vibrations, so when the motion gyro is self-oscillating, chaotic motion can appear. The problem of the constrained integral action in the controller is also considered in this paper. It is well known that the dynamic range of actuators is limited, so when the input reaches the saturation limit, the integrator in the controller will continue integrating while the input is constrained. Consequently, the state of the integrator can reach an unacceptable high value that can produce a poor transient response. In order to avoid this undesirable effect, an anti-windup integrator is considered, so the integral action remains constant from the saturation of the actuator as soon as the actuator saturates [16–18]. When the integral action is constrained, the system has three equilibrium points P1 , P2 and P3 . Point P1 is identified as the set point, and the analysis of equilibrium points P2 and P3 reveals that two interesting bifurcation phenomena can appear. Considering a constrained value of the integral action, it is possible to prove that in the equilibrium point P2 the system has double-zero eigenvalues, and in the equilibrium point P3 the system has a pair of pure imaginary eigenvalues. Consequently, a Bogdanov–Takens bifurcation at point P2 , and a Poincare–Andronov–Hopf bifurcation with a weak focus at point P3 are the responsible of a complicated behaviour of the system [12,13,19,20]. In this paper, we show from the analysis of integral action of the PID controller and the Poincare–Andronov–Hopf bifurcation, that it is possible to predict the motion of the gyro from the first Lyapunov value. The center manifold theorem and the normal forms have been used in order to obtain the first Lyapunov value [19,20]. From a fixed value of the integral action, it is possible to prove that the system only has one reachable equilibrium point. In this case the first Lyapunov value is used to determine if this equilibrium point is a stable or unstable weak focus. When all the equilibrium points are unstable, the system reaches a self-oscillation motion. In this situation an external vibration of the platform can give chaotic behaviour, which has been researched from the sensitive dependence, Lyapunov exponents and the power spectral density.
2. Mathematical model Consider a gyroscope in gimbal schematically shown in Fig. 1. The following notation has been used: OX1 Y1 Z1 and OX2 Y2 Z2 are respectively the internal and external gimbal ring reference systems, and OXYZ is the gyro reference system. OZ1 is the gyro spin axis, OX2 is the external gimbal rotation axis, c is the gyro spin angle in the internal gimbal, b is the internal gimbal rotation angle, and a is the external gimbal rotation angle. Taking into account that the principal axis of the rotor are: A ¼ B, C and Ai , Bi , Ci (i ¼ 1; 2) are the moment of inertia of the internal and external gimbal respectively, and denoting pi , qi , ri the angular velocity respect to the OXi , OYi , OZi axis respectively, the following kinetic energies can be written: • Kinetic energy of the external gimbal: T2 ¼ 12ðA2 p22 þ B2 q22 þ C2 r22 Þ;
q2 ¼ r2 ¼ 0;
p2 ¼ a_ ) T2 ¼ 12A2 a_ 2
ð1Þ
• Kinetic energy of the internal gimbal: T1 ¼ 12ðA1 p12 þ B1 q21 þ C1 r12 Þ;
p1 ¼ a_ cos b;
T1 ¼ 12ðA1 a_ 2 cos2 b þ B1 b_ 2 þ C1 a_ 2 sin2 bÞ
_ q1 ¼ b;
r1 ¼ a_ sin b
ð2Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1059
Fig. 1. Physical model of the gimbal suspension gyro and reference frames.
• Kinetic energy of the rotor: T0 ¼ 12ðAp2 þ Bq2 þ Cr2 Þ; p ¼ a_ cos b; 2 2 2 1 _ _ 2Þ T0 ¼ ðAa_ cos b þ Bb þ Cð_c þ a_ sin bÞ
_ q ¼ b;
r ¼ c_ þ a_ sin b
ð3Þ
2
where c_ is the spin velocity of the rotor, and the dot denotes differentiation respect to time. Considering the projections of the angular velocities in the inertial reference system On0 s0 f0 and considering negligible friction, the whole system kinetic energy has the form: T ¼ 12½J ðbÞa_ 2 þ Hr b_ 2 þ Cð_c þ a_ sin bÞ2
ð4Þ
where Hr ¼ A þ B1 , and J ðbÞ is defined by the equation: J ðbÞ ¼ A2 þ ðA1 þ AÞ cos2 b þ C1 sin2 b
ð5Þ
Eq. (5) has the following physical interpretation: Without the term A cos2 b, the value of J ðbÞ represents the total moment of inertia of the external and internal gimbal rings respect to the rotation axis OX2 . Considering that a, b, c are independent and holonomic variables, and applying Euler–Lagrange equations [5,7,8] the gyro motion equations without control are the following: 9 d da dc da > > J ðbÞ þ C þ sin b sin b ¼ 0 > > > dt dt dt dt > > > = 2 2 d b 1 dJ ðbÞ da dc da da ð6Þ Hr 2 C þ sin b cos b ¼ 0 > dt 2 db dt dt dt dt > > > > > d dc da > > C þ sin b ¼0 ; dt dt dt Eqs. (6) represent the exact gyroscope nutation theory [5,7]. Setting A1 þ A ¼ C1 in Eq. (5), Eq. (6) can be simplified as follows: 9 d da > J þ H sin b ¼ 0 > = dt dt ð7Þ 2 > db da > Hr 2 H cos b ¼ 0 ; dt dt
1060
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 2. Gyroscope and PID controller feedback control system.
where H ¼C
dc da þ sin b ; dt dt
J ¼ A þ A1 þ A2
ð8Þ
The value of H is the gyro kinetic moment, which is a constant obtained from the initial conditions of the motion, and J ðbÞ J is constant. In order to use the gyro as a control device, a feedback control system is added as shown in Fig. 2. The internal gimbal ring angle is sensed, and the output signal of the set performed by the potentiometer, operational amplifier and PID controller is applied to a DC motor. R, L, Kb and Km are the motor resistance, inductance, back-electromotive constant and the motor torque constant respectively. Considering the potentiometer and the circuit of the DC motor the following equations can be deduced: V2 ðtÞ ¼ Kpo bðtÞ; L
V1 ¼ Kpo br ;
Vs ðtÞ ¼ V2 ðtÞ V1 ðtÞ ¼ Kpo ½bðtÞ br ð9Þ
diðtÞ du þ RiðtÞ ¼ V ðtÞ Kb dt dt
where Kpo is the potentiometer constant, u is the giro angle of the shaft motor and V ðtÞ is the outlet voltage of the amplifier–PID controller. Taking into account that the gear train is considered as ideal, from the balance of the torque we get the following equations: MðtÞ ¼ Im
d2 u Ma ; dt2 j
MðtÞ ¼ Km iðtÞ;
d2 a Ma ¼ j Jm 2 jKm iðtÞ; dt 2
N2 j¼ N1
du da ¼ j dt dt
ð10Þ
where Im is the moment of inertia of motor shaft, Ma is the torque due to the feedback control and j is the tooth ratio of the gear train. Considering that equilibrium torque M acts in the OX2 , Eqs. (7) and (8) are the model of the gyro when: M þ Ma ¼ 0
ð11Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1061
In the general case, Eq. (11) is not verified, so the mathematical model of the gyro can be written as follows: J
d2 a db ¼ Ma þ M þ H cos b dt2 dt
d2 b da H cos b ¼ 0 2 dt dt dc da þ sin b ; H ¼C dt dt
ð12Þ
Hr
J ¼ A þ A1 þ A2
Eqs. (12) represent the gyroscope model with the torque Ma due to the feedback control system and the equilibrium torque M . In Eqs. (12) we assume that the platform of the whole system is static. However if the platform has an external oscillation hðtÞ, Eqs. (10) are not valid, because another torque appear in the On0 axis. In this case the relative velocity of the shaft motor is: dðu hÞ dða hÞ ¼ j dt dt
ð13Þ
where the h angle is measured respect to an horizontal static plane. Taking into account Eq. (13), the value of Ma can be rewritten as follows: Ma ¼ j2 Jm
d2 a d2 h þ jðj þ 1ÞJ jKm iðtÞ m dt2 dt2
ð14Þ
and consequently, from Eqs. (12) and (14) we get: 9 d2 a db d2 a d2 h 2 > > = ¼ j þ H cos b J þ jðj þ 1ÞJ þ M m m dt2 dt dt2 dt2 > d2 b da > ; Hr 2 H cos b ¼ 0 dt dt
J
ð15Þ
The amplifier–PID controller equation can be written as follows: Z 1 t dbðtÞ ðbðsÞ br Þ ds þ sd V ðtÞ ¼ Ka Kp bðtÞ br þ si 0 dt
ð16Þ
where Ka is the amplifier constant, Kp is the proportional constant (including the potentiometer constant Kpo ), si is the reset time, sd is the derivative time, and br is the reference angle of the internal gimbal ring. Note that if si ! 1 and sd ! 0, Eq. (13) represents the most simple strategy of control, i.e. control P [17,18]. From Eqs. (15) and (16) the equations of the system with PID feedback control and external vibration of the platform are the following: 9 d2 aðtÞ dbðtÞ d2 hðtÞ > > ¼ jK þ H cos b iðtÞ þ M þ jðj þ 1ÞJ > m m > dt2 dt dt2 > > > > 2 > > d bðtÞ daðtÞ > > ¼ 0 Hr H cos bðtÞ = dt2 dt Z t > 1 dbðtÞ > > > V ðtÞ ¼ Ka Kp bðtÞ br þ ðbðsÞ br Þ ds þ sd > > si 0 dt > > > > > diðtÞ daðtÞ dhðtÞ > ; þ RiðtÞ ¼ V ðtÞ þ jKb jKb L dt dt dt JT
ð17Þ
where JT ¼ J þ j2 Jm . The equations of the mathematical model show that the a variable is cyclic, so it is suitable to remove it by using Eq. (17). Introducing the state variables: x1 ðtÞ ¼ bðtÞ;
x2 ðtÞ ¼
dbðtÞ ; dt
x3 ðtÞ ¼ iðtÞ;
x4 ðtÞ ¼
diðtÞ ; dt
x5 ðtÞ ¼
Z o
t
ðbðsÞ br Þ ds
ð18Þ
1062
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
the following equations can be written:
9 x_ 1 ðtÞ ¼ x2 ðtÞ > > > > > H Ka Kp > > cos x1 ðtÞ Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ x5 ðtÞ x_ 2 ðtÞ ¼ > > > jKb Hr si > > > > H > _ > þ cos x1 ðtÞhðtÞ > > > Hr > > > > x_ 3 ðtÞ ¼ x4 ðtÞ = 2 j Km Kb R jKb H sd Ka Kp H > cos x1 ðtÞ x3 ðtÞ x4 ðtÞ x2 ðtÞ cos x1 ðtÞ þ x_ 4 ðtÞ ¼ > > L LJT jKb Hr L LJ > > T > > Ka Kp Ka Kp > > Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ x5 ðtÞ þ ½x1 ðtÞ x1r > > > si Lsi > > > > Ka Kp jKb jKb s HK K > d a p > € _ > x2 ðtÞ þ þ M þ ½jðj þ 1ÞIm JT hðtÞ þ cos x1 ðtÞhðtÞ > > L LJT LJT LHr > > ; x_ 5 ðtÞ ¼ x1 ðtÞ x1r
ð19Þ
It is interesting to point out that depending on the values of Kp , sd , si , it is possible to research different control strategies using the exact mathematical model of the gyro.
3. Analysis of the equilibrium points when the integral action is constrained If the integral action is constrained, Eqs. (19) can be simplified taking into account that it is assumed x5 ¼ Im , where Im is the maximum value that the PID controller integral action reaches. Assuming that the platform of the system is static, Eqs. (19) can be rewritten as follows: 9 x_ 1 ðtÞ ¼ x2 ðtÞ > > > > > Ka Kp > > Im x_ 2 ðtÞ ¼ ac cos x1 ðtÞ Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ > > si > > > > > x_ 3 ðtÞ ¼ x4 ðtÞ = ð20Þ x_ 4 ðtÞ ¼ bc x3 ðtÞ cc x2 ðtÞ cos x1 ðtÞ þ dc cos x1 ðtÞ > > > Ka Kp Ka Kp > > Im þ ½x1 ðtÞ x1r > Rx3 ðtÞ Ka Kp ðx1 ðtÞ x1r Þ sd Ka Kp x2 ðtÞ > > si Lsi > > > > Ka Kp jKb > > ; M x2 ðtÞ þ þ L LJT where ac bc cc and dc are parameters defined by the following equations: ac ¼
H ; jKb Hr
bc ¼
j2 Km Kb ; LJT
cc ¼
jKb H ; LJT
dc ¼
sd Ka Kp H sd Ka Kp ¼ ac jKb Hr L L
ð21Þ
It is interesting to point out that while the integral action is smaller than the maximum value Im , it is determined by Eq. (19), but when the value Im is reached, the change of the state variables is defined by Eq. (20). So, Eq. (19) can be considered as a transitory state for the variable x5 ðtÞ. Let us consider Eq. (20) when the integral action has reached the steady state value Im . It is clear that in the equilibrium state, the state variable x2 and x4 are always zero, consequently three different equilibrium points are possible: • Point P1 (x1e ; x2e ¼ 0, x3e ; x4e ¼ 0) In this case, considering fixed values of x1r , M , Kp and si , the value of x1e and x3e can be obtained from the following equations: 9 Ka Kp > Im ¼ 0 > = si Ka Kp jKb > ; ðx1e x1r Þ þ M ¼ 0> bc x3e þ Lsi LJT
Rx3e Ka Kp ðx1e x1r Þ
ð22Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
• Point P2 (x1e ¼ p=2, x2e ¼ 0, x3e ; x4e ¼ 0) Where the value of x3e can be found from the following equation: jK Ka Kp p b p=2 þ M ¼ 0 x1r þ bc x3e Lsi 2 LJT
1063
ð23Þ
• Point P3 (x1e ¼ 3p=2, x2e ¼ 0, x3e ; x4e ¼ 0) In this case, the value of x3e can be found from the following equation: Ka Kp 3p jKb 3p=2 M ¼0 bc x3e þ x1r þ Rsi RJT 2
ð24Þ
Note that for n integer, and x1e ¼ ð4n þ 1Þp=2, x1e ¼ ð4n þ 1Þ3p=2 we also obtain the equilibrium points P2 and P3 respectively. On the other hand, when the integral action is constrained to an arbitrary fixed value Im the following inequalities at points P2 and P3 are verified: 9 p KK a p p=2 > x1r Ka Kp Im 6¼ 0 > P2 : Rx3e = 2 si ð25Þ 3p Ka Kp > > x P3 : Rx3p=2 K K I ¼ 6 0 ; a p 1r m 3e 2 si Nevertheless, it is possible to find the value of Im at points P2 and P3 . Taking into account Eqs. (22), at point P2 (p=2) the following equations are verified: 9 p KK a p p=2 > > x Rxp=2 K K I ¼ 0 = a p 1r 3e 2 si ð26Þ K K p jK > a p b ; x1r þ bc xp=2 M ¼ 0 > 3e þ Lsi 2 LJT Eliminating xp=2 3e by using Eqs. (26) one obtains: si RM RJT p x1r si 2 I p=2 ¼ 2 Ka Kp jKm j Km Kb
ð27Þ
The same reasoning allows us to obtain the integral action which belongs to equilibrium point P3 : si RM RJT 3p x1r si 2 I 3p=2 ¼ 2 Ka Kp jKm j Km Kb
ð28Þ
Note that from Eqs. (27) and (28) t is deduced that if x1r < p=2 then I 3p=2 < I p=2 . Once the equilibrium points are determined, by considering the parameter values of the PID controller and from certain initial conditions, the system will evolve to point P1 if the eigenvalues of the Jacobian matrix of the linearized system have negative real part, and points P2 and P3 are unreachable. Consequently it is necessary to research the behaviour of the system at points P1 , P2 and P3 . The analysis of the equilibrium points is carried out transforming the fixed points to the origin by the translation: x01 ðtÞ ¼ x1 ðtÞ x1e ;
x02 ðtÞ ¼ x2 ðtÞ x2e ;
x03 ðtÞ ¼ x3 ðtÞ x3e ;
x04 ðtÞ ¼ x4 ðtÞ x4e
ð29Þ
3.1. Analysis of point P1 Substituting Eq. (29) into Eq. (20) and taking into account the equilibrium equations (22) we get the following equations: 9 x_ 01 ðtÞ ¼ x02 ðtÞ > 0 > > 0 0 0 0 0 > x_ 2 ðtÞ ¼ ac cosðx1 ðtÞ þ x1e Þ Rx3 ðtÞ þ Lx4 ðtÞ Ka Kp x1 ðtÞ sd Ka Kp x2 ðtÞ > > > 0 0 > x_ 3 ðtÞ ¼ x4 ðtÞ = ð30Þ c x_ 04 ðtÞ ¼ dc cosðx01 ðtÞ þ x1e Þ Rx03 ðtÞ þ Lx04 ðtÞ Ka Kp x01 ðtÞ sd Ka Kp þ c x02 ðtÞ > > > dc > > > > Ka Kp 0 Ka Kp 0 R > ; x2 ðtÞ bc x03 ðtÞ x04 ðtÞ þ x1 ðtÞ þ Lsi L L
1064
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
The equilibrium point of Eq. (30) is the origin, and the eigenvalues of the Jacobian matrix are defined by the following equations: k4 þ b3 k3 þ b2 k2 þ b1 k þ b0 ¼ 0 R b2 ¼ bc þ Lac cc cos2 x1e ; b3 ¼ ; L R cos x1e b0 ¼ ac Ka Kp bc Lsi
b1 ¼
Rac cc cos x1e þ sd Ka Kp ac cc ac
Ka Kp si
cos x1e
ð31Þ
From Routh–Hurwitz criterion [17] it is deduced that point P1 will be stable if the roots of the characteristic equation have negative real part, i.e. when p=2 < x1e < p=2, b2 b3 > b1 and b1 ðb2 b3 b1 Þ > b0 b23 . 3.2. Analysis of point P2 Considering Eq. (30), substituting x1e ¼ p=2, and taking Im ¼ I p=2 , Im > 0 in Eq. (27), the Jacobian matrix of the system (30) has two zero eigenvalues, so a Takens–Bogdanov bifurcation can appear. Assume that the integral action is constrained to Im ¼ I p=2 e with e > 0, Im > 0, then the system (30) has a pure imaginary pair of eigenvalues and so point P2 will be a weak focus. In contrast to a stable focus in a weak focus the convergence of the trajectory is not exponential, i.e. the trajectory is not a logarithmic spiral and its length tends to infinity as the radius of the trajectory tends to zero [19]. However, if the integral action is constrained to Im ¼ I p=2 þ e with e > 0, Im > 0 it is easy to show that the equilibrium point P2 is a saddle, and consequently, this point will be unreachable. From the previous reasoning, it is deduced that while the integral action (x5 ðtÞ) is smaller than I p=2 the system will evolve to point P2 if this point is a stable weak focus. On the other hand, if P2 is an unstable weak focus the system will jump to another equilibrium point, P1 or P3 . Finally, when the integral action reaches a value greater than I p=2 the point P2 will be unreachable and the system will change to point P1 or P3 . Consequently point P2 will be reachable when the following conditions are satisfied: • The integral action is constrained to Im ¼ I p=2 e with e > 0, Im > 0 and a weak focus with a pure imaginary pair of eigenvalues appear. • The weak focus is stable. 3.3. Analysis of point P3 Assuming that the integral action is constrained to Im ¼ I p=2 þ e, with e > 0, Im > 0 and substituting x1e by 3p=2 in Eq. (30), the equations of the system are the following: 9 x_ 01 ðtÞ ¼ x02 ðtÞ > > > > x_ 02 ðtÞ ¼ ac sin x01 ðtÞ½Rx03 ðtÞ þ Lx04 ðtÞ Ka Kp x01 ðtÞ sd Ka Kp x02 ðtÞ þ p > > > > x_ 03 ðtÞ ¼ x04 ðtÞ = c ð32Þ x_ 04 ðtÞ ¼ dc sin x01 ðtÞ Rx03 ðtÞ þ Lx04 ðtÞ Ka Kp x01 ðtÞ sd Ka Kp þ c x02 ðtÞ þ p > > > dc > > > Ka Kp 0 Ka Kp 0 R > > ; x2 ðtÞ bc x03 ðtÞ x04 ðtÞ þ x1 ðtÞ þ L Lsi L where 3p=2
p ¼ Rx3e
Ka Kp
3p x1r 2
Ka Kp p=2 ðI þ eÞ si
ð33Þ
and I p=2 between (24), (27) and (33), the following equation is deduced: Eliminating x3p=2 3e JT R Ka Kp p ¼ Ka Kp p 1 2 e ð34Þ si j Km Kb si pffiffiffiffiffiffiffiffiffiffi It is easy to show that if p < 0, Eq. (32) have eigenvalues k1;2 ¼ i ac jpj and consequently it is necessary to research the equilibrium state of point P3 carefully. 4. Analysis of the stability of point P3 From the conclusions of the previous section it can be shown that the equilibrium point P3 is a weak focus. The stability of a weak focus for systems of high dimension can be difficult to determine [12,13,19–21]. In the case of point
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1065
P3 the system has four dimensions, so it will be necessary to carry out a set of previous steps in order to put Eq. (32) in an appropriate form, to determine the first Lyapunov value [19,20]. These steps are summarized as follows: Step 1 Make a Taylor expand of the function sin x01 ðtÞ up to order 3. Substituting in Eq. (32), we next split off the linear and nonlinear part of the vector field to get 2 3 2 32 0 3 2 3 0 0 x1 0 1 0 0 6 x_ 10 7 6 0 7 6 0 0 0 0 7 6 x_ 2 7 6 c 0 0 0 7 76 x20 7 þ 6 f2 ðx1 ; x2 ; x3 ; x4 Þ 7 þ Oðkx0 k4 Þ 6 0 7¼6 ð35Þ 54 x 5 4 5 6 x_ 7 4 0 0 0 0 1 3 4 30 5 0 0 0 0 0 f4 ðx1 ; x2 ; x3 ; x4 Þ x_ 4 x4 a b bc n where a ¼ dc p þ
c ¼ ac p;
Ka Kp ; Lsi
b¼
Ka Kp ; L
n¼
R L
ð36Þ
9 = f2 ¼ ac x01 ðRx03 þ Lx04 Ka Kp x01 sd Ka Kp x02 Þ a6c p x033 h i cc dc p 03 0 0 0 0 0 f4 ¼ dc x1 Rx3 þ Lx4 Ka Kp x1 sd Ka Kp þ dc x2 6 x3 ;
ð37Þ
The eigenvalues of linear part of (35) are: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffi k2;3 ¼ n=2 i bc ðn=2Þ2 k1;2 ¼ ix ¼ i ac jpj;
ð38Þ
and consequently an Andronov–Poincare–Hopf bifurcation appears [21,22]. Step 2 Let P be the matrix that transforms transformation: 2 3 2 x01 0 1 6 07 6 6x 7 6 x 0 6 27 x0 ¼ P x ) 6 0 7 ¼ 6 6 x3 7 6 4 5 4 g f x04
f ¼
aðbc x2 Þ þ bnx2 ðbc
x 2 Þ2
Eq. (35) becomes 2 3 2 x_ 0 6 y_ 7 6 x 6 7¼6 4 z_ 5 4 0 w_ a
2
þn
x 0 0 b
x2
0 0 a1 b1
xf
xg
;
g¼
the linear part of (35) into Jordan canonical form. Then, under the linear 0 0 0 b1
32
3 x 76 7 0 76 y 7 76 7; 6 7 1 7 54 z 5 0
a1
R a1 ¼ ; 2L
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 R b1 ¼ bc 2L
ð39Þ
w
x½bðbc x2 Þ an ðbc x2 Þ2 þ n2 x2
2 3 32 3 0 0 x 6 7 6 7 0 7 76 y 7 þ P 1 6 f2 ðx; y; z; wÞ 7 þ Oðkxk4 Þ 4 5 b1 54 z 5 0 a1 f4 ðx; y; z; wÞ w
ð40Þ
where f2 and f3 are the same functions of Eq. (37) but substituting coordinates x0 by x. Step 3 We now simplify Eq. (35) as much as possible using the center manifold theorem [12,13,19,20] in order to reduce the dimensionality to two in a neighbourhood of the equilibrium point P3 . In the center manifold the variables z and w are approximated by: 2 z h ðx; yÞ a x þ a2 y 2 þ a3 xy ¼ 1 2 ð41Þ ¼ 1 h2 ðx; yÞ w b1 x þ b2 y 2 þ b3 xy Coefficients ai , bi (i ¼ 1; 2; 3) are determined by substituting Eq. (41) into the equation which defines the center manifold [12,13,20]: x þ f ðx; y; h1 ðx; yÞ; h2 ðx; yÞÞ ¼ Bhðx; yÞ þ gðx; y; h1 ðx; yÞ; h2 ðx; yÞÞ Dx hðx; yÞ½A
ð42Þ
1066
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
where " # x x ¼ ; y " B¼
" ¼ A
a1
b1
b1
a1
# ;
" # 1 f2 ðx; y; hðx; yÞÞ ; f ðx; y; hðx; yÞÞ ¼ x x 0 0 2 3 ðga1 xf Þf2 f4 þ 7 6 oh oh b1 x b1 7 6 g¼6 Dx h ¼ ; 7; 4 5 ox oy gf 0
x
#
ð43Þ
2
x We then equate powers of x2 , xy, y 2 in the left and right sides of Eq. (42), so coefficients ai , bi can be computed. The simplified system at the center manifold in normal form can be written as follows: 1 f2 ðx; yÞ 0 x x x_ ¼ þ ð44Þ y_ x 0 y 0 x where f2 ðx; yÞ ¼ a02 y 2 þ a11 xy þ a03 y 3 þ ac ðR þ La1 Þyðb1 x2 þ b2 y 2 þ b3 xyÞ þ ac Lb1 yða1 x2 þ a2 y 2 þ a3 xyÞ ac p a02 ¼ ac ðRf Lxf Ka Kp Þ; a11 ¼ ac ðRg þ Lxf sd Ka Kp xÞ; a03 ¼ 6 Step 4 Considering a system of the form: 1 0 x x x_ f ðx; yÞ þ 2 ¼ y x 0 y_ f ðx; yÞ
ð45Þ
ð46Þ
there is an explicit equation that gives the first Lyapunov value L1 [12,13,19–22]: i 1 1 1 h 1 1 1 2 2 fxy ðfxx þ fyy1 Þ fxy2 ðfxx2 þ fyy2 Þ fxx1 fxx2 þ fyy1 fyy2 L1 ¼ ½fxxx þ fxyy þ fxxy þ fyyy þ 16 16x
ð47Þ
Applying Eq. (47) to Eqs. (45), (46), and considering that all partial derivatives are evaluated at the bifurcation point P3 (x ¼ 0, y ¼ 0), the value of L1 is: L1 ¼
a11 a02 ac ðR þ La1 Þb3 ac Lb1 a3 þ þ 8x3 8x 8x
ð48Þ
so, only values of a3 and b3 are necessary to know the first Lyapunov value. The value of parameters a3 and b3 obtained from the center manifold theorem (Step 3) are the following: g a3 ¼ F1 ðp1 c11 Þ þ F2 p2 a11 x g b3 ¼ F2 ðp1 c11 Þ þ F1 p2 a11 x
ð49Þ
where F1 ¼
a1 þ a1 f1 2
ða1 þ a1 f1 Þ þ ðb1 b1 f1 Þ
2
;
F2 ¼
b1 b1 f1 2
ða1 þ a1 f1 Þ þ ðb1 b1 f1 Þ
9 ðga1 xf Þa02 b02 ga02 > > þ b a þ p1 ¼ 2x > 1 1 bc > b1 x b1 x > > > > > > ðga xf Þa b ga 1 02 02 02 > > þ þ a b p2 ¼ 2x 1 > 1 bc > b1 x b1 x > = b02 ¼ dc ðRf Lxg Ka Kp Þ > > > ðga1 xf Þa11 b11 > > þ c11 ¼ > > > b1 x b1 > > > > > a11 ¼ ac ðRg þ Lxf sd Ka Kp Þ > > ; b11 ¼ dc ðRg þ Lxf sd Ka Kp Þ cc x
2
;
f1 ¼
4x2 þ b21
a21
ð50Þ
ð51Þ
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1067
If L1 > 0 the equilibrium point P3 will be unstable. When L1 < 0 the equilibrium point will be stable, and if L1 ¼ 0 it is necessary to reduce the system (35) on the center manifold to the formal norm up to the terms of fifth order [19–22]. Consequently Eqs. (45), (46) and (48)–(51) give us a criterion to know the stability proprieties of point P3 . 5. Regular and self-oscillating behaviour The analysis of the equilibrium points reveals that when the integral action is constrained to a value given by Im ¼p=2 þe the equilibrium point P2 will be unreachable and if the first Lyapunov value defined by Eqs. (48)–(51) is positive, the equilibrium point P3 will be an unstable weak focus. Consequently, the system changes to equilibrium point P1 if all roots of the characteristic equation (31) have negative real part and there is no external vibration of the platform. Therefore, if one chooses adequate values for the parameter Kp , si , sd of the PID controller it is possible to reach a desired set point x1e . In this case we say that the behaviour of the gyro is regular. Nevertheless, it is possible to obtain other movements when the gyro changes to equilibrium points P2 , P3 or reaches a self-oscillating state. Note that it is possible to choose different combinations of the parameters of the PID controller in order to obtain a desired set point. In this section several examples are carried out to examine different behaviours of the gyro by numerical simulation. The Runge–Kutta–Felhberg method has been used in numerical simulations with integration steps between 0.001 and 0.0005 s. The numerical values and the parameters of the gyro are shown in Table 1. We first consider the dynamics of the gyro when the set point x1e and the reference x1r are different. The simulation results are shown in Figs. 3 and 4. The parameters of the PID controller Kp , si and sd are chosen from the standard procedures [17,18,23] taking into account that for fixed values of si , sd the proportional constant Kp must give a positive first Lyapunov value, so the equilibrium point P3 (3p=2) will be an unstable weak focus as shown in Fig. 3. The value e ¼ 0:5 > 0 gives an equilibrium point P1 of 1.297 rad and the integral action is constrained to 8.79 s, as shown in Fig. 4. Note that the eigenvalues of point P1 have negative real part, and consequently this point will be reachable. The moment of inertia of the motor shaft must be Jm . Description and values are correct. The table is the following: Table 1 Gyroscope and control system parameter values Variable
Description
Value
A1 B1 C1 A2 A C Hr J Jm j JT nr H Km Kb R L Ka Kp si sd I Im ac bc cc dc M
Internal gimbal moment of inertia respect to OX1 (kg m2 ) Internal gimbal moment of inertia respect to OY1 (kg m2 ) Internal gimbal moment of inertia respect to OZ1 (kg m2 ) External gimbal moment of inertia respect to OX2 (kg m2 ) Rotor equator moment of inertia (kg m2 ) Rotor polar moment of inertia (kg m2 ) A þ B1 (kg m2 ) A þ A1 þ A2 (kg m2 ) Moment of inertia of the motor shaft (kg m2 ) N2 =N1 Tooh ratio J þ j2 Jm (kg m2 ) Rotor angular velocity (rad/s) C nr Gyrokinetic moment (kg m2 rad/s) Motor torque constant (N m/A) Back-electromotive constant (V s) Motor resistance (X) Motor inductance (H) Amplifier constant Proportional constant of the PID controller (V) Reset time of the PID controller (s) Derivative time of the PID controller (s) Integral action (s) Maximum value of the integral action (s) Parameter H =jKb Hr (rad/V s) Parameter j2 Km Kb =LJT (s1 ) Parameter jKb =LJT (A) Parameter sd Ka Kp =jKb Hr L (A V/s) Disturbance torque (N m)
0.045 0.045 0.05625 0.04875 0.01125 0.0225 0.05625 0.105 0.001 10 0.205 104.72–523.59 2.3562–11.781 5 0.04 0.5–1 0.01 10 0.01–0.05 1–2 0.2–0.5
9756.097 195.122 10–200
1068
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 3. Simulation results for transient and set point steady state response from the following values: nr ¼ 1500 r.p.m.; M ¼ 200 N m; x1r ¼ p=8 rad; e ¼ 0:5; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; Im ¼ 8:7937. Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0, x5 ð0Þ ¼ 0. Eigenvalues at the set point: [)95.8 ± 48.5j, )3.656 ± 1.357j]. First Lyapunov value at P3 : L1 ¼ 0:02569 > 0.
Fig. 4. Steady state of the PID controller integral action. Simulation values from Fig. 3.
Similar results are obtained when the set point or equilibrium point P1 coincides with the reference angle br ¼ x1r . It is important to remark that if the value of the reference angle br ¼ x1r is greater than a certain value, from the last equation of (19), it is possible that the integral action reaches a negative value. Consequently, the steady state of the integral action will be Im . In this case, applying the same procedure as before, the equilibrium point P2 is a stable weak focus and the gyro has not regular behaviour, i.e. is not possible to reach the set point. Fig. 5 shows this case. Fig. 6 shows a case in which the first Lyapunov value is negative, so point P3 is a stable weak focus. In this case it is not possible to reach the equilibrium point. Note that while the integral action x5 ðtÞ is lower than Im point P2 is a stable weak focus, but when the limit of the integral action is reached, the gyro jumps to point P3 . The self-oscillating behaviour is researched from the previous consideration of Section 3. From Eqs. (22) and (26), if the integral action is constrained to I p=2 þ e with a small positive value of e, then the equilibrium point P1 will be almost coincident with the equilibrium point P2 . The basins of attraction of P1 and P2 collide and by considering the discussion of Section 3, these points will be unreachable. So, the only possibility of the gyro is to evolve to point P3 . However, point P3 is an unreachable weak focus if the first Lyapunov value is positive. So, if L1 > 0 a self-oscillatory behaviour
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1069
Fig. 5. Simulation results for transient and steady state response when the PID controller integral action is constrained to a negative value. Point P2 is a stable weak focus. Simulation values: nr ¼ 1500 r.p.m.; M ¼ 100 N m; x1r ¼ 2 rad; e ¼ 105 ; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; I p=2 ¼ 7:8088; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0, x5 ð0Þ ¼ 0.
Fig. 6. Simulation results when the first Lyapunov at point P3 is negative L1 ¼ 0:5 and P1 P2 . Simulation values: nr ¼ 2000 r.p.m.; M ¼ 100 N m; x1r ¼ p=8; e ¼ 105 ; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; I p=2 ¼ 4:88; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0; x5 ð0Þ ¼ 0. Eigenvalues at the set point P1 x1e p=2 [)25 ± 95.56j; )6.73 · 105 ± 172.4j].
appears. This behaviour has been corroborated by simulation as shown in Fig. 7. Note that if the integral action is constrained to an adequate value it will always possible to chose parameters of the PID controller in order to obtain L1 > 0.
6. Chaotic behaviour Nowadays, the theory of chaotic oscillations is considered as an interesting field of research. The qualitative features of many models of oscillation have been previously applied to systematic experimental and theoretical developments of mechanical oscillators, both in numerical and experimental studies [24–27]. In dissipative electromechanical systems, the presence of a strange attractor, which appears when the system has chaotic behaviour, is of paramount importance in order to understand new nonlinear proprieties of the electromechanical device. Nowadays, in order to carry out theoretical and experimental studies of a strange attractor, the following steps are recognized as necessary [28,29]:
1070
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 7. Self-oscillating behaviour when the first Lyapunov value at point P3 is positive L1 ¼ 0:01 and P1 P2 . Simulation values: nr ¼ 1000 r.p.m.; M ¼ 100 N m; x1r ¼ p=8; e ¼ 105 ; Kp ¼ 0:02 V; si ¼ 2 s; sd ¼ 0:5 s; I p=2 ¼ 7:65; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0; x5 ð0Þ ¼ 0. Eigenvalues at the set point P1 x1e p=2 [)25 ± 95.56j; )2.59 · 105 ± 1.02 · 102 j].
ii(i) The development of mathematical models based on plausible subsystems of the electromechanical system. i(ii) The analysis and design of the whole system in order to understand strange behaviours. (iii) The possibility of simulating the mathematical model with a clear criterion to elucidate between regular, transient chaotic and chaotic behaviour. With these ideas the study of chaotic behaviour of the gyro is carried out taking into account the model given by Eq. (19). When the platform of the gyro has an external vibration movement, the relationship of the platform external disturbance can be written as follows: _ ¼ Ap xp cos xp t; hðtÞ ¼ Ap sin xp t ) hðtÞ
€hðtÞ ¼ Ap x2 sin xp t p
ð52Þ
The values of the amplitude Ap and the angular frequency xp must be properly chosen in order to obtain chaotic behaviour. Previously, the system is simulated with Ap ¼ 0, xp ¼ 0 and the integral action constrained in order to obtain a first Lyapunov value L1 > 0. In this case, the behaviour of the system will be self-oscillating. Fig. 8 shows the value of L1 as
Fig. 8. Variation of the first Lyapunov value vs. Kp and self-oscillating behaviour with L1 > 0 and P1 P2 . Simulation values: nr ¼ 1200 r.p.m.; M ¼ 100 N m; x1r ¼ p=8; e ¼ 105 ; Kp ¼ 0:026 V; si ¼ 1:827 s; sd ¼ 0:4568 s; I p=2 ¼ 4:88; Im I p=2 . Initial conditions: bð0Þ ¼ 0, dbð0Þ=dt ¼ 0, ið0Þ ¼ 0:01 A, x4 ð0Þ ¼ 0; x5 ð0Þ ¼ 0. Eigenvalues at the set point P1 x1e p=2 [)25 ± 95.56j; )4.04 · 105 ± 1.335 · 105 j].
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1071
Fig. 9. Strange attractor when the platform has an external harmonic disturbance. Simulation values of Fig. 8. Amplitude Av ¼ 0:065 m, frequency xv ¼ 8:6 rad/s.
function of the proportional constant Kp of the PID controller. Note that for Kp ¼ 0:026 the value of L1 is positive, and the point P3 is an unstable weak focus. The period of the self-oscillation is close to 2.65 s. When the system has self-oscillating behaviour, an external disturbance in the platform of the gyro as shown in Eq. _ (52), can produce chaotic behaviour. Fig. 9 shows the strange attractor drawn in the bðtÞ–bðtÞ plane. The zones with strong stretching and folding as well as zones with Cantor’s sets are shown in Fig. 9. On the other hand, it is well known that it is very difficult to show rigorously that a strange attractor is chaotic. It can be said that a nonlinear system has chaotic behaviour if (a) it has sensitive dependence on initial conditions, i.e. two very close initial conditions diverge exponentially with time, (b) the orbits are dense in a state region i.e. the orbits fill the phase space zone of the strange attractor, and (c) the orbits are topologically transitive in i.e. for any two open sets X1 ; X2 X there is a time from which any orbit starting at X1 ends at X2 [13]. While the conditions (a) and (b) can be approximately verified by simulation, proving the condition (c) is generally very difficult. Nevertheless, the fulfilment of (a) and (b) can be enough to assure the long time chaotic behaviour. Fig. 10 shows that the system has sensitive dependence to initial conditions, and Fig. 11 shows Lyapunov’s exponents of the system. The algorithm that has been used can be found in [26,27,30,31]. Eqs. (19) are transformed in an autonomous dynamical system considering a new variable x6 ðtÞ defined by: x6 ðtÞ ¼ xp t )
dx6 ðtÞ ¼ xp dt
Fig. 10. Sensitive dependence and error vs. time with two initial conditions very close. Simulation values of Fig. 8.
ð53Þ
1072
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
Fig. 11. Lyapunov’s exponents for the chaotic strange attractor of Fig. 9.
Note the presence of a positive Lyapunov exponent, which is an indicator of chaos, and a Lyapunov exponent nil, corresponding to x6 ðtÞ periodic variable. In order to verify the computations, the divergence of vector field from Eq. (19) is calculated as follows: div~ f ¼ sd Ka Kp
H R sd Ka Kp H cos x1 ðtÞ cos x1 ðtÞ þ jKb Hr L jKb H
ð54Þ
then, the sum of Lyapunov exponents are determined, and the following relation: div~ f ¼
6 X
ð55Þ
NLi
i¼1
is verified by simulation. The results obtained with an integration step of 0.004 s are the following: div~ f ¼ 50:13;
6 X
NLi ¼ 50:01
ð56Þ
i¼1
therefore the computations can be considered correct. Fig. 12 shows another indicator of chaos, i.e. the Fourier power spectrum for the variable ðtÞ ¼ x2 ðtÞ shown at Fig. 10. The power spectrum is a representation of the relative abundance of different frequencies in a given time series. Note that the power spectrum is broadband and contains substantial at low frequencies.
Fig. 12. Power spectrum for the bðtÞ variable of Fig. 9.
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
1073
Finally, it is important to point out that although both a positive Lyapunov exponent and broad spectrum do not guarantee sensitivity to initial conditions it is, in practice, an indicator of chaos. The straight line AB has been drawn through the harmonics peaks to remark the character of chaotic dynamic, probably due to a route of period-doubling sequences [12,13,26–29]. 7. Conclusions From the results of the research presented in this paper, it has been corroborated that a gyroscope in gimbal with a feedback PID controller, can be operated by using a strategy that allows to reach a desired set point, self-oscillation dynamics as well as a more complex behaviour, characterized by chaotic dynamics. By using the exact nutation theory of the gyroscope and the PID controller with constrained integral action, we show that the equations of the gyroscope and the feedback control system have multiple equilibrium points. The analysis of these equilibrium points reveals that depending on a fixed value of the PID controller integral action three different equilibrium points may appear. One of them, point P1 , with positive or negative real part of eigenvalues, another point P2 (p=2), with a pure imaginary pair or two zero eigenvalues, and a third point P3 (3p=2), with a pure imaginary pair of eigenvalues. For a certain integral action value, we show that point P2 is unreachable, and point P3 is a weak focus. Then, the stability of point P3 , by using the manifold center theorem, the normal form theory and the first Lyapunov value is researched. The sign of the first Lyapunov value gives us a clear criterion in order to elucidate if point P3 is a stable or unstable weak focus. The determination of normal form reveals that the algebraic machinery can be formidable, nevertheless it follows straightforward rules. Advances in both digital and symbolic computational methods can be investigated in connection with the analysis of complex system with several equilibrium points, as it is shown in this paper. When points P1 and P2 are unreachable, if the first Lyapunov value at point P3 is negative, this point will be a stable weak focus, as corresponds to an Andronov–Poincare–Hopf bifurcation with a pure imaginary pair of eigenvalues. On the other hand, if the first Lyapunov value is positive, point P3 is unreachable too and the gyro has a self-oscillating motion. The self-oscillating behaviour allows us to find chaotic oscillations when the platform is submitted to an external harmonic disturbance so, a new family of strange attractors has been found. The chaotic behaviour has been researched from the sensitive dependence to initial conditions, Lyapunov’s exponents and the power spectral density, from which it is possible to verify that the mechanism of chaotic oscillations is due to a period-doubling phenomenon. Another interesting aspect of this paper is the possibility of designing an algorithm, which gives the PID parameters values i.e. the proportional constant Kp , the reset time si , and the derivative time sd , in order to obtain a previous desired set point. This algorithm will be based in the eigenvalues of point P1 with negative real part, point P2 with a saddle and point P3 with a positive first Lyapunov value. The previously analyzed phenomena should be known by the mechanical designer in order to either avoid or properly use them, depending on gyro type. The paper shows a research line, where more elaborated models of gimbal suspension gyro with structural flexibility, and more sophisticated control strategies can be used to plan both theoretical and experimental studies. References € Klein F, Sommerfield A. Uber die theorie des kreisels. H. 1–4. Leipzig-Berlin: Teubner, 1897/1898–1903/1910. Timoshenko S, Young DH. Advanced dynamics. New York: McGraw-Hill; 1948. p. 392–448. Pars LA. A treatise on analytical dynamics. New York: John Wiley & Sons; 1966. p. 159–87. Cabannes H. Cours de Mecanique Generale. Paris: Dunod; 1966. p. 248–69. Goldstein H. Classical mechanics. New York: Addison Wesley; 1980. Marsden JE, Ratiu SI. Introduction to mechanics and symmetry. New York: Springer Verlag; 1999. p. 481–512. Chetaev NG. Theoretical mechanics. Moscow: Mir; 1989. p. 221–5. Arnold VI. Mathematical methods of classical mechanics. New York: Springer Verlag; 1974. p. 75–90, 157–95. Cannon RH. Dynamics of physical systems. New York: McGraw-Hill; 1967. p. 159–63, 617–26. Singh SN. Stability of gyro with harmonic nonlinearity in spinning vehicle. IEEE Trans Aerospace Electron Syst 1983;19:182–9. Ge ZM, Chen HH. Bifurcations and chaos in a rate gyro with harmonic excitation. J Sound Vib 1996;194(1):107–17. Guckenheimer J, Holmes P. Nonlinear oscillations, dynamical systems and bifurcations of vector fields. New York: Springer Verlag; 1983. p. 117–80. [13] Wiggins S. Introduction to applied nonlinear dynamical systems and chaos. New York: Springer Verlag; 1990. p. 193–202, 211–29, 253–84. [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12]
1074
M.F. Perez Polo, M.P. Molina / Chaos, Solitons and Fractals 21 (2004) 1057–1074
[14] Ge ZM, Chen HH. Double degeneration and chaos in a rate gyro with feedback control. J Sound Vib 1998;209(5):753–69. [15] Agafonov SA. Gimbal suspension gyro: stability, bifurcation, and chaos. J Dyn Control Syst 2001;7(3):339–51. [16] Kothare M, Campo P, Morari M, Nett C. A unified framework for the study of anti wind-up designs. Automatica 1994;30(12):1869–83. [17] Franklin GF, Powell JD, Emani-Naeni A. Feedback control of dynamics systems. New York: Prentice-Hall; 2000. p. 95–100. [18] Goodwin GC, Graebe SF, Salgado ME. Control system design. Saddle River, New Jersey: Prentice-Hall; 2001. p. 159–76. [19] Shilnikov LP, Shilnikov AL, Turaev DV, Chua LO. Methods of qualitative theory in nonlinear dynamics. Part II. New Jersey: World Scientific; 2001. p. 451–65, 598–610. [20] Kuznetsov YA. Elements of applied bifurcation theory. New York: Springer Verlag; 1998. p. 86–109, 151–71. [21] Glendinning P. Stability, instability and chaos: an introduction to the theory of nonlinear differential equations. New York: Cambridge University Press; 1996. p. 224–43. [22] Marsden JE, McCraken M. The Hopf bifurcation and its applications. New York: Springer Verlag; 1976. [23] Ogata K. Modern engineering control. New York: Prentice-Hall; 2000. [24] Rudiger S. Practical bifurcation and stability analysis. From equilibrium to chaos. New York: Springer Verlag; 1994. p. 327–57. [25] Wggins S. Global bifurcations and chaos. Analytical methods. New York: Springer Verlag; 1988. p. 75–108. [26] Lichtenberg AJ, Lieberman MA. Regular and chaotic dynamics. New York: Springer Verlag; 1992. p. 62–5, 312–20, 457–78. [27] Holden AV, editor. Chaos. Princeton, New Jersey: Princeton University Press; 1986. [28] Hilborn RC. Chaos and nonlinear dynamics. Oxford: Oxford University Press; 2000. p. 121–48, 323–41. [29] Ott E. Chaos in dynamical systems. Cambridge: Cambridge University Press; 2002. p. 115–37, 212–36, 305–10. [30] Benettin G, Galgani L, Giorgilly A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them, Part I: Theory. Meccanica 1980;15:9–20. [31] Benettin G, Galgani L, Giorgilly A, Strelcyn JM. Lyapunov characteristic exponents for smooth dynamical systems and for Hamiltonian systems; a method for computing all of them, Part II: Numerical applications. Meccanica 1980;15:21–30.