Micro-chaotic Behaviour in PD-controlled Systems

Micro-chaotic Behaviour in PD-controlled Systems G´ abor Csern´ ak G´ abor St´ ep´ an HAS-BUTE Research Group on Dynamics of Machines and Vehicles, Budapest, Hungary (e-mail: [email protected]) † Department of Applied Mechanics, Budapest University of Technology and Economics, Budapest, Hungary (e-mail: [email protected]) ∗

Abstract: The nonlinearity caused by the application of digital control may lead to chaotic behaviour. Our goal is to analyse the problem of computer-controlled stabilization of unstable equilibria, with the application of the PD control law. We consider a linear equation of motion. However, as a consequence of the digital effects, i.e., the sampling and the round-off error, the solutions of the system can be described by a two dimensional piecewise linear map. We show that this system may perform chaotic behaviour. The amplitude of the evolving oscillations is usually very small, this is why these are called micro-chaotic vibrations. Keywords: digital control, micro-chaos, proof of chaos, two-dimensional map 1. INTRODUCTION

Parallel with the fast development of computer technology, more and more commercial products are operated using digital control loops. Unfortunately, several engineers are not aware of the effects of the discretization in time (sampling), the discretization of the measured data (round-off error), and the processing delay, however, these digital effects may lead to chaotic behaviour, as it was shown by Kuo (1977), Chen and Dong (1998), Delchamps (1990), and Haller and St´ep´ an (1996). Our goal is to analyse the problem of computer-controlled stabilization of unstable equilibria, with the application of the PD control law. One of the most frequently studied models in this field is the digitally controlled inverted pendulum on a cart. Several researchers, e.g. Landry et. al (2005), Sieber and Krauskopf (2004), and Kelly et. al (2005) dealt with this, or similar problems. Enikov and St´ep´an (1998) constructed an experimental device for the examination of the system and found some stochastic vibrations around the desired upper equilibrium position of the pendulum. PD control was applied, and the goal of the control was to stabilize the pendulum, while the position of the cart could be arbitrary. They took into account all the three digital effects, and they found that the solutions of the equation of motion could be described by a three-dimensional piecewise linear map. Great efforts were made in order to prove that the map is chaotic, but the proof is not complete yet. In the present contribution, we neglect the processing delay but consider the effects of sampling and round-off error. The resulting 2D map is more manageable than the original 3D one, this is why we can show that the corresponding map is chaotic.

2. MATHEMATICAL MODEL As a first step, we take into account the effect of temporal discretization (sampling) only. We assume that the position and velocity are measured periodically, with period τ , and the obtained values are used to determine the control force, without delay. In this case, the control force is kept constant between two successive sampling instants. Thus, the equation of motion of the examined system assumes the following form: Mx ¨(t) − kx(t) = −P xj − Dx˙ j , {z } |

t ∈ [jτ, (j + 1)τ ). (1)

uj

Here M denotes mass, −k is the stiffness, τ is the sampling time, xj = x(jτ ), and x˙ j = x(jτ ˙ ), while P and D are the proportional and differential gains, respectively. Dividing the equation by the coefficient of the acceleration, we obtain x ¨(t) − β 2 x(t) = −γ(P xj + Dx˙ j ),

(2)

where β 2 = k/M and γ = 1/M . The general solution of equation (2) assumes the following form in t ∈ [jτ, (j + 1)τ ): γ x(t) = C1 sinh(βt) + C2 cosh(βt) + 2 (P xj + Dx˙ j ). (3) β Since the coefficients C1 and C2 can be determined according to the initial and matching conditions, one can construct a linear two-dimensional map xj+1 = Axj , where xj = col[xj , x˙ j ] and  ch − 1 ch − γP  β2 A= sh βsh − γP β

sh ch − 1  − γD β β2  sh  . ch − γD β

(4)

(5)

Here we used the notations ch ≡ cosh(βτ ), and sh ≡ sinh(βτ ). The originally unstable fixed point of the system becomes asymptotically stable if the eigenvalues of A are inside the unit circle. Applying Jury’s criterions, i.e., the Moebiustransformed Routh-Hurwitz criterions, one finds that the conditions of asymptotic stability are the following: β2 (6) γ β(1 + cosh(βτ )) D < a0cond ≡ (7) γ sinh(βτ ) cosh(βτ ) − 1 D > a1cond ≡ P. (8) β sinh(βτ ) As it is shown in Fig. 1, the domain of stability formes a triangle on the P -D parameter plane. In the following, we will restrict ourselves to this parameter domain. P > a2cond ≡

D

ˆ ˆ ˆ j + sinh(β) y˙ ′ − cosh(β) − 1 m (11) yj+1 = cosh(β)y j βˆ βˆ2 ′ ˆ ˆ j + cosh(β)y ˆ ′ − 1 sinh(β)m. (12) yj+1 = βˆ sinh(β)y j βˆ To have a more concise form, let us introduce new variables and notations:  ˆ    ˆ sinh(β) yj cosh( β) , (13) yj = ′ , B =  βˆ yj ˆ cosh(β) ˆ βˆ sinh(β)   ˆ 1 − cosh(β)     Pˆ   βˆ2 b= , c = (14)  ˆ ˆ .  D sinh(β)  − βˆ Using these notations, map (11,12) can be rewritten as a vector-vector map f : f : yj+1 = Byj + bInt (c · yj ) . (15) It can be checked easily that det B = 1 and the eigenvalues ˆ while the eigenvectors are of B are λ1,2 = exp(±β), ˆ s1,2 = col[1, ±β].

a0cond S

The fixed points of the map (15) are given by the following formula   m −1 m r = (I − B) bm = col ,0 , (16) βˆ2

a1cond a2cond

where I denotes the unit matrix. P

Fig. 1. The asymptotic stability of the continuous system. The stable domain disappears if the three border lines have a common intersection point. This situation arises if βτ = 0. Since the parameter β is strictly positive, we find that – with appropriate gains P and D–, the equilibrium can be asymptotically stabilized at every positive value of the sampling time τ . The output signal has a finite resolution h in the digitally controlled system. Taking into account this quantization, too, the control force cannot be given as in (2), but   P xj + Dx˙ j x ¨(t) − β 2 x(t) = −γhInt ≡ −γhm, (9) h where m is an integer number. Let us introduce the nondimensional time T = t/τ , and the ′ new variables y = τ 2xγh and y ′ = τ 2xγh , where []′ denotes differentiation with respect to T . Using these notations, Equation (9) can be rewritten in nondimensional form:   ˆ ′ ≡ −m, y ′′ (T ) − βˆ2 y(T ) = −Int Pˆ yj + Dy (10) j ˆ = Dγτ . where βˆ = βτ , Pˆ = P γτ , and D 2

The solutions of this equation can be expressed by the following 2D map, instead of (4):

Thus, the domain of definition of the map (15) can be divided into parallel to the value of the   bands, according ˆ ′ . These bands are separated number m ≡ Int Pˆ yj + Dy j by the lines m − Pˆ yj yj′ = . (17) ˆ D Fig. 2 shows these lines and bands, together with a numerically obtained solution representing the attractor of the system. The fixed point rm must be in band m, thus |m| (|m| + 1) |m| ≤ < . Pˆ Pˆ βˆ2

(18)

The first inequality implies Pˆ ≥ βˆ2 , that corresponds to (6), one of the conditions of asymptotic stability of the single fixed point of the continuous system (2). The second inequality implies that the discrete map (15) has 2mmax +1 fixed points if βˆ2 mmax < . (19) Pˆ − βˆ2 In the example shown in Fig. 2, βˆ = 1 and Pˆ = 1.5, thus, mmax < 2 must be fulfilled. Consequently, mmax = 1 and there are 3 fixed points, at yj = −1, yj = 0, and yj = 1. The eigenvectors, corresponding to the stable and unstable directions of the fixed point at the origin, are also shown

2

m=0

m=1

m=2

m=3

m=4

∞ X k |y∞ | = max − A bχk .

m=5

1.5

To give an approximation for the sum in (23), we transformed the matrix A into diagonal form:   µ1 0 ˜ A≡ = T−1 AT, (24) 0 µ2

1 0.5

yj′

0

−0.5 −1 −1.5 −2 −3

(23)

k=0

m = −5

m = −4

−2

m = −3

m = −2

−1

m = −1

m=0

0

1

yj

2

3

ˆ = 1 and Fig. 2. The attractor of the system at βˆ = D ˆ P = 1.5. in Fig. 2. It is clearly seen in the figure, that the points of the attractor are distributed along threads, parallel with the unstable direction. According to (19), there was an infinity of fixed points in the case Pˆ = βˆ2 – all of them at the borders between the bands –, while P ≥ 2βˆ2 would mean that the origin was the single fixed point. 3. PROOF OF CHAOS 3.1 Existence of an Attractor Map (15) can be rewritten as yj+1 = (B + b ◦ c)yj − bχj = (B + b ◦ c)j y0 −

j−1 X

(B + b ◦ c)k bχk , (20)

k=0

where −1 < χj < 1. It can be observed that B + b ◦ c ≡ A.

(21)

According to numerical evidence, the numbers χk vary irregularly and the sequence (20) is neither convergent, nor divergent for almost all initial vectors y0 . However, for sufficiently large index j, the vectors yj stay in the neighbourhood of the origin. The maximal possible norm |y∞ | of the vectors yj in the limit j → ∞ provides an estimation for the size of the attracting domain at the origin.

˜ ≡ [b1 , b2 ]T = T−1 b, and transformed the vector b into b where matrix T consists of the eigenvectors of matrix A, and µi , i = 1, 2 are the eigenvalues. The maximal possible norm can be given as ∞ ∞ X X k k ˜ k ˜ bχ |y∞ | ≡ max A bχk = max T A k=0 k=0 ∞  X µk1 b1 χk    k=0   . = max T  (25) ∞ X   k  µ b χ 2 2 k k=0

Note that – since both µi and bi are complex numbers, (i = 1, 2) –, the products µki bi correspond to vectors on the complex plane, with different orientations in the cases of the different powers k. Thus, the choice χk ≡ 1, k = 0, 1, 2, . . . would not lead to the maximal norm. We use the following estimation, instead:  ∞  X k |µ1 | |b1 |    k=0   . |y∞ | ≤ kTk  (26) ∞ X   k  |µ2 | |b2 | k=0 ˆ = 1, and In the example shown above, with βˆ = 1, D Pˆ = 1.5, the transformed matrix and vector assume the forms   0.5482 + 0.5822i 0 ˜ A≈ (27) 0 0.5482 − 0.5822i   ˜ ≈ −0.456 + 0.5876i , b (28) −0.456 − 0.5876i and |y∞ | is estimated as |y∞ | ≤ 8.61712.

(29)

In our numerical simulations – as it can also be seen in Fig. 2 – the maximal norm of yj was approximately 2.78, that corresponds to this result.

Since A is not a normal matrix, its Euclidean norm is given by the greatest singular value σ1 , which is not equal to the spectral radius ρ. Thus, the assumption ||A|| = ρ is not valid and cannot be exploited to approximate the limit |y∞ |, as it was done in Enikov and St´ep´ an (1998).

Since the moduli of eigenvalues of A are less than one in the considered cases, |y∞ | is necessarily finite. Thus, we can claim that there is a finite attracting domain – an absorbing sphere – in the neighbourhood of the origin.

Although the greatest singular value – the norm of A – can be greater than one,

3.2 Sensitivity to Initial Conditions

lim ||Aj || = lim ρj = 0

j→∞

j→∞

is fulfilled, since ρ < 1, according to (6)-(8). Thus,

(22)

ˆ thus, the The eigenvalues of B are λ1,2 = exp(±β), distance of two solutions with neighbouring initial points taken from the same band increases exponentially. Since the unstable manifolds are transversal to the border lines

1.5 1 0.5

yj′ 0 −0.5 −1 −1.5

11111 00000 00000 11111 00 11 00000 11111 00000 11111 00 111111111111 11 00000 11111 000000000000 0000000000000 1111111111111 0000000 1111111 000000000000 111111111111 0000000000000 1111111111111 0000000 10 0 11111111 0 000000000000 111111111111 00 11 1 0000000000000 1111111111111 000000000000 111111111111 00 11 000 111 0000000000000 0001111111111111 111 0000000000000 1111111111111 V0

W0s

V1

W1u

F0

H1

H0

V0

W1s

F1

V1

Switching line

W0u

0

0.2

0.4

yj

0.6

0.8

1

ˆ = 1 and Pˆ = 1.5. Fig. 3. Parallelograms H0,1 and their images V0,1 form a Smale horseshoe at βˆ = D between the bands, we can assume without loss of generality that after n iterations the two solutions arrive at different bands. Once the points are separated in this manner, their orbits can no longer be considered close, since they behave essentially independently.

3.3 Topological Transitivity and Smale’s Horseshoes ˆ = 1, Pˆ = 1.5, there In our example parameter set βˆ = D are three fixed points of saddle type on the (yj , yj′ ) plane: F−1 (−1, 0), F0 (0, 0), and F1 (1, 0). These fixed points are found in three different bands, according to the integer ˆ ′ . Fixed points F0 and F1 are shown part of Pˆ yj + Dy j s in Fig. 3, together with their stable (W0,1 ) and unstable u (W0,1 ) manifolds, and the switching line, forming the border between the bands m = 0 and m = 1. As it can be seen, the examined map is similar to the well-known baker map, examined by Farmer et. al (1983), and conjugate to a modified version of the Smale horseshoe map, introduced by Smale (1967). We placed two parallelograms H0 and H1 along the stable manifolds of the fixed points. The sides are parallel with either the stable, or the unstable manifolds. The coordinates of the rightmost vertex of H1 are (1.04,-0.03), while the lengths of the sides are defined as 0.5 times and 0.08 times the length of the eigenvector ˆ T , respectively. The leftmost vortex of H0 is s1 = [1, β] found at (0.15,0.1), and the lengths of the sides are 0.52|s1 | and 0.1|s1 |, respectively. We constructed the second images of H0 and H1 : V0 = f 2 (H0 ), and V1 = f 2 (H1 ). Due to the existence of the switching line, both images consist of two parts. As it is shown in the figure, V0 and V1 fully intersect H1 , and V1 fully intersects H0 . The intersection between H0 and V0 is only partial, but the obtained structure is still topologically conjugate to a Smale horseshoe. The third images of H0,1 already intersect both parallelograms, several times.

The existence of a Smale horseshoe implies sensitivity on initial conditions, the map is topologically transitive and there exists a countable infinity of periodic orbits, an uncountable infinity of nonperiodic orbits and a dense orbit. Thus, the examined map is chaotic in a finite parameter domain in the neighbourhood of the examined parameters. Note, that this proof could be performed in more details using interval arithmetic, as it was done by Mischaikov and Mrozek (1995), and Zgliczinski (1996) in the case of other maps. 4. CONCLUSION We showed in a simple example that the digital effects due to computer control may lead to chaotic vibrations of a PD-controlled system. The used parameters cannot be considered to be typical in engineering applications. Thus, in the future we would like to perform a similar calculation for a realistic control system, as well. ACKNOWLEDGEMENTS This research was supported by the Hungarian National Science Foundation under grant no. OTKA F049242, and by the J´anos Bolyai Research Scholarship of the Hungarian Academy of Sciences. REFERENCES B.C. Kuo, Digital Control Systems, SRL Publishing, Champaign, IL, USA, 1977. G. Chen, X. Dong, From Chaos to Order: Methodologies, Perspectives and Applications, World Scientific, Singapore, 1998 F.D. Delchamps, Stabilizing a linear system with quantized state feedback, IEEE Trans. Autom. Contr., 35 pp. 916-924, 1990. Haller, G. – St´ep´an, G.: Micro-Chaos in Digital Control, J. Nonlinear Sci., Vol. 6, 1996, 415-448

M. Landry, S.A. Campell, K. Morris and C.O. Aguilar, Dynamics of an inverted pendulum with delayed feedback control, SIAM Journal on Applied Dynamical Systems 4(2) pp. 333-351, 2005. J. Sieber, B. Krauskopf, Complex Balancing Motions of an Inverted Pendulum Subject to Delayed Feedback Control, Physica D 197(3-4), pp. 332-345, 2004. R. Kelly, V. Santibanez, A. Loria, Control of Robot Manipulators in Joint Space, Springer, London, 2005. E. Enikov, G. St´ep´an, Micro-Chaotic Motion of Digitally Controlled Machines, J. of Vibration and Control, 4, pp. 427-443, 1998. Mischaikov K and Mrozek M. Chaos in the Lorenz equations: a computer-assisted proof. Bull. Amer. Math. Soc., 32, pp. 66-72, 1995. Zgliczinski P. Fixed point index for iterations of maps, topological horseshoe and chaos. Topol. Methods. Nonlin. Anal., 8, pp. 169-177, 1996. J.D. Farmer, E. Ott, J.A. Yorke, The Dimension of Chaotic Attractors, Physica, 7D pp. 153-180, 1983. S. Smale, Differentiable Dynamical Systems, Bulletin of AMS, 73, 1967