Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electronic Journal of Differential Equations, Conference 10, 2003, pp 163–185. http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Dynamics of delayed piecewise linear systems ^{∗}

László E. Kollár, Gábor Stépán, & János Turi

Abstract

In this paper the dynamics of the controlled pendulum is investigated assuming backlash and time delays. The upper equilibrium of the pendulum is stabilized by a piecewise constant control force which is the linear combination of the sampled values of the angle and the angular velocity of the pendulum. The control force is provided by a motor which drives one of the wheels of the cart through an elastic teeth belt. The contact between the teeth of the gear (rigid) and the belt (elastic) introduces a nonlinearity known as “backlash” and causes the oscillation of the controlled pendulum around its upper equilibrium. The processing and sampling delays in the determination of the control force tend to destabilize the controlled system as well. We obtain conditions guaranteeing that the pendulum remains in the neighborhood of the upper equilibrium. Experimental findings obtained on a computer controlled inverted pendulum cart structure are also presented showing good agreement with the simulation results.

1 Introduction

Dynamical systems with piecewise linear components in the equations of motion occur frequently in practice. Gear pairs with backlash [23], railway wheelsets with clearance [13], mechanical oscillators with clearance [6, 14, 15, 19] or amplitude constraints [20, 21], impact dampers [16, 17, 18], moving parts with dry friction [12] and adjacent structures during earthquake [4, 5] are modelled by systems with piecewise linear stiffness, damping or restoring force.

Control is often added to such systems to stabilize unstable equilibria. Digital balancing of the inverted pendulum is examined in the following sections. The pendulum is placed on a cart and controlled by a motor which drives one of the wheel-axle of the cart by an elastic teeth belt [2, 7]. Backlash occurring at the driving wheel of the motor makes the system piecewise linear and creates an unstable zone around the otherwise stable equilibrium of the controlled system. Consequently, the best possible outcome of the stabilization process is to keep the pendulum in a small neighborhood of its upper equilibrium.

∗Mathematics Subject Classifications: 34K35, 37C75.

Key words: Stability analysis, backlash, digitally controlled pendulum.

\begin{matrix}c \\ ◯ \end{matrix}2003 Southwest Texas State University.

Published February 28, 2003.

Figure 2.1: The inverted pendulum on a cart and its stability map

In Section 2, a model of the inverted pendulum on a cart is derived. Stability analysis is performed using the assumption that the belt is perfectly elastic, (i.e. no backlash). Stability chart for control parameters is also obtained and parameterized by the stiffness of the driving belt. In Section 3, backlash at the contact of the driving belt and the axle of the motor is considered. The effect of backlash and time delay together is investigated in Section 4. Experimental observations are presented in Section 5.

2 The inverted pendulum on a cart

The model of the inverted pendulum on a cart can be seen in Figure 2.1, [8, 9, 10, 11]. The cart can move only in the horizontal direction. The motor drives one of the wheels of this cart through a driving-belt with stiffness s. The system has three degrees of freedom described by the general coordinates q_{1}, q_{2}, q_{3}, denoting the displacement x of the cart, the angle φ of the pendulum and the angle ψ of the motor axle, respectively. It is assumed that the angle φ of the pendulum and the displacement x of the cart are observed together with their derivatives and the observed values provide the motor voltage U_{m} according to the equation

U_{m} = P_{1}x + D_{1}ẋ + P_{2}φ + D_{2}φ̇ ,  (2.1)

where the constants P_{1}, D_{1}, P_{2}, D_{2} are the so-called PD controller parameters. U_{m} together with ψ̇, the angular velocity of the motor determine M_{m}, the driving torque of the motor by

M_{m} = LU_{m} − Kψ̇ ,  (2.2)

where L and K are given motor parameters.

The nonlinear equations of motion are obtained (see also Figure 2.1) by means of the Lagrange equations of the second kind. The kinetic energy T , the potential U, the dissipation D and the general force Q have the following form

T = \frac{1}{2}m ((ẋ + \frac{l}{2}φ̇ cos φ)^{2} + (\frac{l}{2}φ̇ sin φ)^{2}) + \frac{1}{24}ml^{2}φ̇^{2} + \frac{1}{2}Mẋ^{2} + \frac{1}{2}J_{m}ψ̇^{2},

U = −\frac{1}{2}mgl (1 − cos φ) + U_{s} ,

D = 0, Q_{1}^{e} = 0, Q_{2}^{e} = 0, Q_{3}^{e} = M_{m} ,

(2.3)

where the index e refers to elastic, U_{s} = 0 if the belt is ideally rigid, U_{s} = sΔ^{2}/2 if the belt is elastic and Δ is the elongation of the spring, M = m_{m} + m_{c} + 4m_{w} + 4J_{w}/R_{w}^{2}, the sum of the mass of the motor m_{m}, the cart m_{c} and the reduced mass of inertia of the wheels 4m_{w} + 4J_{w}/R_{w}^{2}. Applying the Lagrange equations of the second kind

\frac{d}{dt} \frac{∂T}{∂q̇_{k}} − \frac{∂T}{∂q_{k}} + \frac{∂U}{∂q_{k}} + \frac{∂D}{∂q̇_{k}} = Q_{k} ,  (2.4)

k = 1, 2, 3, the nonlinear equations of motion for the controlled pendulum are obtained

(\begin{matrix}m + M & \frac{1}{2}ml cos φ & 0 \\ \frac{1}{2}ml cos φ & \frac{1}{3}ml^{2} & 0 \\ 0 & 0 & \frac{1}{2}m_{m}r_{m}^{2} \end{matrix}) (ẍ ; φ̈ ; ψ̈)

+ (\begin{matrix}s\frac{r_{w}^{2}}{R_{w}^{2}} & 0 & −sr_{m} \frac{r_{w}}{R_{w}} \\ 0 & 0 & 0 \\ −sr_{m} \frac{r_{w}}{R_{w}} & 0 & sr_{m}^{2} \end{matrix} ) (\begin{matrix}x \\ φ \\ ψ \end{matrix}) + (\begin{matrix}−\frac{1}{2}mlφ̇^{2} sin φ \\ −\frac{1}{2}mgl sin φ \\ 0 \end{matrix}) = (\begin{matrix}0 \\ 0 \\ Q_{3}^{e} \end{matrix}) .

(2.5)

The case when the belt is ideally rigid

If the belt is ideally rigid [8], then x and ψ are related by the equation ψ = xr_{w}/(r_{m}R_{w}), where r_{w} and R_{w} are radii of the wheel and r_{m} is the radius of the motor axle and then ψ can be eliminated from equation (2.5). Multiplying the third equation of (2.5) by r_{m}r_{w}/R_{w}, adding it to the first equation and substituting ψ = xr_{w}/(r_{m}R_{w}) yields the nonlinear equations of motion

(\begin{matrix}m_{r} & \frac{1}{2}ml cos φ \\ \frac{1}{2}ml cos φ & \frac{1}{3}ml^{2} \end{matrix}) (ẍ ; φ̈) − (\begin{matrix}\frac{1}{2}mlφ̇^{2} sin φ \\ \frac{1}{2}mgl sin φ \end{matrix}) = (\begin{matrix}Q_{1}^{r} \\ 0 \end{matrix}) ,  (2.6)

where m_{r} = m + M + m_{m}r_{w}^{2}/(2R_{w}^{2}), m and l are the mass and the length of the pendulum and Q_{1}^{r} = Q_{3}^{e}r_{w}/(r_{m}R_{w}) with index r referring to rigid. Linearizing these equations around the upper equilibrium φ = 0 we obtain the variational system

(\begin{matrix}m_{r} & \frac{1}{2}ml \\ \frac{1}{2}ml & \frac{1}{3}ml^{2} \end{matrix}) (ẍ ; φ̈) − (\begin{matrix}0 & 0 \\ 0 & \frac{1}{2}mgl \end{matrix}) (\begin{matrix}x \\ φ \end{matrix}) = (\begin{matrix}Q_{1}^{r} \\ 0 \end{matrix}) .  (2.7)

Using equations (2.1), (2.2) and (2.3), the control force has the form

Q_{1}^{r} = \frac{r_{w}}{r_{m}R_{w}} L (P_{1}x + D_{1}ẋ + P_{2}φ + D_{2}φ̇) − K \frac{r_{w}^{2}}{r_{m}^{2}R_{w}^{2}} ẋ .  (2.8)

It is easy to see from (2.8) that if

P_{1} = 0 and D_{1} = \frac{1}{L} \frac{r_{w}}{r_{m}R_{w}} K ,  (2.9)

then Q_{1}^{r} can be written as

Q_{1}^{r} = \frac{r_{w}}{r_{m}R_{w}} (P φ + Dφ̇) ,  (2.10)

where P = LP_{2} and D = LD_{2}. Thus, the control force Q_{1}^{r} is only a function of φ and φ̇ (i.e., x and ẋ do not appear in Q_{1}^{r}). The displacement x of the cart can then be eliminated from the equations of motion as follows. Solve the first equation of (2.7) for ẍ and substitute the given expression into the second equation of (2.7). A one degree of freedom system is obtained with the following single second order governing equation

φ̈ − \frac{6g}{l} φ = −\frac{6}{m_{r}l} Q_{1}^{r} .  (2.11)

In our setting stabilization means that φ and φ̇ are both zero, but x and ẋ could be nonzero. Due to the relationship between ẍ, φ, φ̇ and φ̈ as described by equation (2.6), we get that ẋ is constant when φ = φ̇ = 0. (I.e., the cart is moving with constant velocity or possibly standing.)

We have the following result concerning finding suitable control parameters P, D in (2.10).

Theorem 2.1 The trivial solution of (2.10)–(2.11) is asymptotically stable if and only if

P > P_{0} = m_{r}gr_{m} \frac{R_{w}}{r_{w}} and D > 0 .  (2.12)

Proof: The characteristic equation of (2.11) has the following form

λ^{2} + \frac{6r_{w}}{m_{r}lr_{m}R_{w}} Dλ + (\frac{6r_{w}}{m_{r}lr_{m}R_{w}} P − \frac{6g}{l} ) = 0 .  (2.13)

The coefficients in (2.13) are clearly positive if the conditions (2.12) are satisfied. The statement of the theorem follows. □

The case when the belt is elastic

If the belt is elastic, then the nonlinear equations of motion assume the form of (2.5). This system can be reduced to a two degrees of freedom system if the assumptions (2.9) are applied again and a new general coordinate, Δ, the elongation of the spring is introduced by

Δ = r_{m}ψ − \frac{r_{w}}{R_{w}} x .  (2.14)

Multiplying the first and third equation of (2.5) by (−m_{m}r_{m}r_{w}) /(2R_{w}) and (m + M), respectively, summing them and using (2.14) gives the first equation of motion of the reduced system. Solving the first equation of (2.5) for ẍ, using (2.14) and substituting it into the second equation of (2.5) gives the second equation of motion of the reduced system and altogether we have

(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} cos φ \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} cos^{2} φ \end{matrix}) (Δ̈ ; φ̈)

+ (\begin{matrix}(m + M)\frac{K}{r_{m}} & 0 \\ 0 & 0 \end{matrix}) (Δ̇ ; φ̇) + (\begin{matrix}\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} φ̇^{2} sin φ \\ \frac{1}{4} \frac{m}{(m+M)}ml^{2}φ̇^{2} sin φ cos φ − \frac{1}{2}mgl sin φ \end{matrix})

+ (\begin{matrix}(m + M) r_{m} + \frac{1}{2}m_{m}r_{m} \frac{r_{w}^{2}}{R_{w}^{2}} \\ \frac{1}{2} \frac{m}{(m+M)}l \frac{r_{w}}{R_{w}} \end{matrix} ) R_{s} = (\begin{matrix}(m + M) Q \\ 0 \end{matrix}) ,

(2.15)

where

R_{s} = (−s\frac{r_{w}}{R_{w}} 0 sr_{m}) (\begin{matrix}x \\ φ \\ ψ \end{matrix}) = sΔ  (2.16)

is the force in the driving belt and

Q = P φ + Dφ̇  (2.17)

is the control force. Linearizing these equations around the upper equilibrium φ = 0 yields

\begin{matrix}(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} \end{matrix}) & (Δ̈ ; φ̈) \\ + (\begin{matrix}(m + M)\frac{K}{r_{m}} & − (m + M) D \\ 0 & 0 \end{matrix}) (Δ̇ ; φ̇) \\ + (\begin{matrix}(m + M) r_{m}s + \frac{1}{2}m_{m}r_{m} \frac{r_{w}^{2}}{R_{w}^{2}} s & − (m + M) P \\ \frac{1}{2} \frac{m}{(m+M)}l \frac{r_{w}}{R_{w}} s & −\frac{1}{2}mgl \end{matrix} ) (\begin{matrix}Δ \\ φ \end{matrix}) & = (\begin{matrix}0 \\ 0 \end{matrix}) . \end{matrix} (2.18)

The subsequent theorem establishes the stability properties of system (2.18).

Theorem 2.2 The trivial solution of system (2.18) is asymptotically stable if and only if

P > P_{0} = m_{r}gr_{m} \frac{R_{w}}{r_{w}} , s > s_{cr} = \frac{3 (m + M) m_{m}g}{(m + 4M + 2m_{m}r_{w}^{2}/R_{w}^{2}) l} , q (P, D, s) < 0 ,

(2.19)

where q (P, D, s) = a_{22}D^{2} + a_{1}P + a_{2}D + a_{0} is a parabola, with the coefficients

a_{22} = 6m_{m}r_{w}R_{w}r_{m}^{3}s ,

a_{1} = 2 (m + 4M) lR_{w}^{2}K^{2} ,

a_{2} = −6 (m + M) m_{m}gR_{w}^{2}r_{m}^{2}K − 2 (m + 4M) lR_{w}^{2}r_{m}^{2}Ks − 4m_{m}lr_{w}^{2}r_{m}^{2}Ks ,

a_{0} = 3mm_{m}glr_{w}R_{w}r_{m}K^{2} .

Proof: The characteristic equation of (2.18) can be written as

b_{4}λ^{4} + b_{3}λ^{3} + b_{2}λ^{2} + b_{1}λ + b_{0} = 0 ,  (2.20)

where

b_{4} =\frac{1}{24} (m + 4M) m_{m}lr_{m}^{2},
b_{3} =\frac{1}{12} (m + 4M) lK,
b_{2} =\frac{1}{12} (m + 4M + 2m_{m} \frac{r_{w}^{2}}{R_{w}^{2}} ) lr_{m}^{2}s − \frac{1}{4} (m + M) m_{m}gr_{m}^{2},
b_{1} =\frac{1}{2}r_{m} \frac{r_{w}}{R_{w}} sD − \frac{1}{2} (m + M) gK,
b_{0} =\frac{1}{2}r_{m} \frac{r_{w}}{R_{w}} sP − \frac{1}{2} (m + M + \frac{1}{2}m_{m} \frac{r_{w}^{2}}{R_{w}^{2}} ) gr_{m}^{2}s .

The Hurwitz matrix assumes the form

(\begin{matrix}b_{1} & b_{0} & 0 & 0 \\ b_{3} & b_{2} & b_{1} & b_{0} \\ 0 & b_{4} & b_{3} & b_{2} \\ 0 & 0 & 0 & b_{4} \end{matrix}) .

According to the Routh-Hurwitz criterion [3], the statement of the theorem holds if all the coefficients of the characteristic equation and all the leading principal minors of the Hurwitz matrix, i.e. the Hurwitz determinants, are positive. The conditions b_{3} > 0 and b_{4} > 0 hold because b_{3} and b_{4} depend only on positive physical parameters. The conditions b_{0} > 0, b_{2} > 0 and H_{3} = b_{1}b_{2}b_{3} − b_{1}^{2}b_{4} − b_{0}b_{3}^{2} > 0 are equivalent to P > P_{0}, s > s_{cr} and q (P, D, s) < 0, respectively. The positivity of the Hurwitz determinants H_{1} = b_{1} and H_{2} = b_{1}b_{2} − b_{0}b_{3} follow from the positivity of the coefficients b_{0}, b_{2}, b_{3}, b_{4} and the Hurwitz determinant H_{3}, because these conditions imply that b_{1}b_{2} − b_{0}b_{3} > b_{1}^{2}b_{4}/b_{3} > 0 and b_{1} > b_{0}b_{3}/b_{2} > 0, i.e., H_{2} > 0 and H_{1} > 0. □

In the remainder of this section first we study the dynamic behavior of the linearized model of the controlled pendulum on the boundary of the stability region determined in Theorem 2.2, and then we present a numerical bifurcation analysis for the original nonlinear system (2.15).

The stability charts in the plane of the control parameters are shown in Figure 2.1 for rigid and elastic belts. The stability domain for rigid belt is a half plane and for the elastic belt it is bordered by a line and a parabola. A real characteristic root changes its sign crossing the line, while a complex conjugate pair turns up in the right half of the complex plane crossing the parabola q (P, D, s) = 0.

The angular frequency, ω, of the oscillation at the loss of stability is the imaginary part of the characteristic roots with zero real parts, thus it can be derived by substituting λ = iω in the characteristic equation (2.20) and solving Figure 2.2: (a) The frequency of the oscillation and the stability domain (b) Bifurcation diagrams for D = 2 [Nms] (upper) and for P = 20 [Nm] (lower)

either its real or imaginary part for ω. Using the fact that q (P, D, s) = 0 we can obtain the frequency of the oscillation as a function of P or D only. For example, ω as a function of P can be obtained as the positive solution of the equation

\frac{1}{24} (m + 4M) m_{m}lr_{m}^{2}ω^{4}

− (\frac{1}{12} (m + 4M + 2m_{m} \frac{r_{w}^{2}}{R_{w}^{2}} ) lr_{m}^{2}s − \frac{1}{4} (m + M) m_{m}gr_{m}^{2}) ω^{2}

+\frac{1}{2}r_{m} \frac{r_{w}}{R_{w}} sP − \frac{1}{2} (m + M + \frac{1}{2}m_{m} \frac{r_{w}^{2}}{R_{w}^{2}} ) gr_{m}^{2}s = 0 .

The frequency f = ω/ (2π) of the oscillation at the loss of stability and the stability domain are drawn in Figure 2.2(a) for the given values of parameters: m = 0.169 [kg], M = 1.136 [kg], m_{m} = 0.2 [kg], g = 9.81 [m/s^{2}], l = 0.5 [m], r_{w} = 0.02 [m], R_{w} = 0.03 [m], r_{m} = 0.01 [m], K = 0.01 [Nms] and s = 10000 [N/m]. Figure 2.2(a) shows that there is a unique frequency for every pair P, D on the parabola.

To study the relationship between (2.18) and (2.15) we performed a numerical bifurcation analysis using the software package AUTO (available at ftp://ftp.cs.concordia.ca/pub/doedel/auto/). Results are presented in Figure 2.2(b) for the values of parameters given in the previous paragraph. The bifurcation parameter is P in the upper half of Figure 2.2(b), and the control parameter D is kept at 2 [Nms]. A pitchfork bifurcation is occurred at P_{0} = 0.1986 [Nm] obtained from (2.19), where the upper equilibrium of the pendulum becomes stable. Figure 2.2(b) indicates that the upper equilibrium maintains its stability till a supercritical Hopf-bifurcation occurs at P_{1} = 139.7 [Nm] obtained from the condition q (P, 2, 10000) = 0. An unstable stationary solution appears at the pitchfork bifurcation and φ tends to π/2 along this solution as P increases. A stable periodic solution appears at the supercritical Hopf-bifurcation.

The bifurcation parameter is D in the lower half of Figure 2.2(b), and the control parameter P is kept at 20 [Nm]. The upper equilibrium is stable between the Hopf-bifurcation points. They occur at D_{1} = 0.1991 [Nms] and D_{2} = 5.916 [Nms] obtained from the condition q (20, D, 10000) = 0. Periodic solution exists at the border of the stability domain where two complex conjugate characteristic roots are crossing the imaginary axis.

If φ is less than its value at the unstable stationary solution then trajectories tend to the stable equilibrium. It follows that if P is chosen from the stability domain and not very close to P_{0} then the upper equilibrium of the pendulum is stabilized even if the initial value of the angle φ is not small. Thus, the linear variational system (2.18) can be considered a good approximation of the original nonlinear system (2.15) in the stability region given by (2.19).

3 Backlash at the gears

Backlash appears in the system as a nonlinearity in the belt-drive. In particular, the force in the belt, R_{s}, is given by

R_{s} = \begin{cases}s (Δ + r_{0}) & if Δ ≤ −r_{0} \\ 0 & if |Δ| < r_{0} \\ s (Δ − r_{0}) & if Δ ≥ r_{0}, \end{cases}  (3.1)

where s and Δ is defined as before and r_{0} is half of the backlash.

Consequently, the equations of motion (2.15) with (3.1) and control force (2.17), after linearization with respect to φ, yield

(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} \end{matrix}) (Δ̈ ; φ̈)

+ (\begin{matrix}(m + M) \frac{K}{r_{m}} & − (m + M) D \\ 0 & 0 \end{matrix}) (Δ̇ ; φ̇)

+ (\begin{matrix}(m + M) r_{m}s + \frac{1}{2}m_{m}r_{m} \frac{r_{w}^{2}}{R_{w}^{2}} s & − (m + M) P \\ \frac{1}{2} \frac{m}{m+M} l \frac{r_{w}}{R_{w}} s & −\frac{1}{2}mgl \end{matrix} ) (\begin{matrix}Δ \\ φ \end{matrix})

= (\begin{matrix}((m + M) r_{m}sr_{0} + \frac{1}{2}m_{m}r_{m} \frac{r_{w}^{2}}{R_{w}^{2}} sr_{0}) sgn Δ \\ \frac{1}{2} \frac{m}{m+M} l \frac{r_{w}}{R_{w}} sr_{0} sgn Δ \end{matrix}) ,

(3.2)

Figure 3.1: (a), (b) Bifurcation diagrams, the bifurcation parameter is P and D, respectively (c) The bifurcation diagram near the point where the computation is interrupted (d) The stability chart with the bifurcation curves

if |Δ| ≥ r_{0}, and

\begin{matrix}(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} \end{matrix}) & (Δ̈ ; φ̈) \\ + (\begin{matrix}(m + M) \frac{K}{r_{m}} & − (m + M) D \\ 0 & 0 \end{matrix}) (Δ̇ ; φ̇) \\ + (\begin{matrix}0 & − (m + M) P \\ 0 & −\frac{1}{2}mgl \end{matrix} ) (\begin{matrix}Δ \\ φ \end{matrix}) & = (\begin{matrix}0 \\ 0 \end{matrix}) , \end{matrix}  (3.3)

if |Δ| < r_{0}.

We have calculated exact (Mathematica) and approximate (MATLAB) solutions for (3.2)–(3.3) using the same parameter values as in Section 2, and r_{0} = 0.001 [m] for the backlash in the system. Our simulations were performed with control parameters chosen from the region given by (2.19). We also performed numerical bifurcation analysis using AUTO. The bifurcation diagrams are shown in Figure 3.1(a) and 3.1(b), where the bifurcation parameters are the control parameters, P and D, respectively. D = 2 [Nms] in Figure 3.1(a) and 3.1(c) (which is an enlargement of a section of Figure 3.1(a)). There is a branch point at one of the borders, P_{0}, of the stability domain (i.e. the domain given by (2.19)), where stable fix points appear. This is the same value which is determined in Theorem 2.2. A homoclinic orbit occurs at P = P_{cr}, so there is a homoclinic bifurcation there, where a limit cycle appears. P_{cr} is approximately 4 according to the AUTO-computations and it is 3.46 according to the MATLAB-simulations. Only the fix points are stable between P_{0} and P_{cr}, while all the fix points and the limit cycle are stable for values of P between P_{cr} and P_{2}. All of them have a domain of attraction, so trajectories spiral to one of them depending on the initial conditions. There is another bifurcation point at P_{2}, where the fix points become unstable. The limit cycle retains its stability with increasing amplitude till the parameter P reaches the other border P_{1} of the stability domain.

P = 40 [Nm] in Figure 3.1(b). There is a stable limit cycle between the borders of the stability domain but the fix points are stable only along the continuous line between the bifurcation points indicated with squares in the figure.

According to the results summarized above, the stability domain in the plane of the control parameters is constructed and sketched in Figure 3.1(d). It has the same boundary as in the case of no backlash (i.e. (2.19)). Only the fix points are stable in a small domain near the line. Stable limit cycle appears at the homoclinic bifurcation points indicated with the dotted line. Fix points lose their stability at the other bifurcation points indicated with the smashed line, so all the fix points and the limit cycle are stable between the dotted and the smashed line, and only the limit cycle is stable in the remaining part of the stability domain.

We note that a periodic solution corresponds to the oscillation of the stick around its vertical equilibrium. The physical meaning of the stable equilibria is that the control force does not push the stick further than the vertical line and it oscillates with decreasing amplitude on either side of the vertical position.

4 The inverted pendulum on a cart with time delay

In this section, sampling, τ_{s} > 0, and processing, τ_{p} > 0, delays are included in the model of the inverted pendulum on a cart [9], [11]. First, we consider the case of no backlash. Then the linearized equations of motion are given by (2.18) as before, but now delayed arguments will appear due to the fact that the control force, Q (t), takes the form

Q (t) = P φ ((j − 1) τ_{s}) + Dφ̇ ((j − 1) τ_{s}) , t ∈ [(j − 1) τ_{s} + τ_{p}, jτ_{s} + τ_{p}) ,

for j = 1, 2, … We shall assume that the sampling and the processing delays are the same, i.e., τ_{s} = τ_{p} = τ. The more general case, τ_{s} ≠ τ_{p}, could be considered similarly. With that, the control force can be written as

Q (t) = P φ ((j − 1) τ) + Dφ̇ ((j − 1) τ) , t ∈ [jτ, (j + 1) τ) , j = 1, 2, …

(4.1)

Without time delay we had a system (2.18) of ordinary differential equations. With control force (4.1) we obtain the following system of delay differential equations

\begin{matrix}(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} \end{matrix}) & (\begin{matrix}Δ̈ & (t) \\ φ̈ & (t) \end{matrix}) \\ + (\begin{matrix}(m + M)\frac{K}{r_{m}} & 0 \\ 0 & 0 \end{matrix}) (\begin{matrix}Δ̇ & (t) \\ φ̇ & (t) \end{matrix}) \\ + (\begin{matrix}(m + M) r_{m}s + \frac{1}{2}m_{m}r_{m} \frac{r_{w}^{2}}{R_{w}^{2}} s & 0 \\ \frac{1}{2} \frac{m}{(m+M)}l \frac{r_{w}}{R_{w}} s & −\frac{1}{2}mgl \end{matrix}) (\begin{matrix}Δ (t) \\ φ (t) \end{matrix}) \\ = (\begin{matrix}(m + M) & (P φ ((j − 1) τ) + Dφ̇ ((j − 1) τ)) \\ & 0 \end{matrix}) , t ∈ [jτ, (j + 1) τ) , \end{matrix} (4.2)

for j = 1, 2, . . .. Transform the second order system (4.2) into a system of first order equations, i.e., multiply (4.2) by the inverse of the coefficient matrix of the second order term and introduce a new state vector y. We obtain

ẏ (t) = Ay (t) + bQ (t) ,  (4.3)

where

\begin{matrix}A = (\begin{matrix}0 & 1 & 0 & 0 \\ −\frac{4r_{w}^{2}s}{(m+4M)R_{w}^{2}} − \frac{2s}{m_{m}} & −\frac{2K}{m_{m}r_{m}^{2}} & \frac{3mgr_{w}}{(m+4M)R_{w}} & 0 \\ 0 & 0 & 0 & 1 \\ −\frac{6r_{w}s}{(m+4M)lR_{w}} & 0 & \frac{6(m+M)g}{(m+4M)l} & 0 \end{matrix}) , \\ b = ( 0 ; \frac{2}{m_{m}r_{m}} ; 0 ; 0 ) , y = (Δ ; Δ̇ ; φ ; φ̇) . \end{matrix}  (4.4)

Note that Q (t) can be expressed as

Q (t) = Ky ((j − 1) τ) , t ∈ [jτ, (j + 1) τ) , j = 1, 2, …

where

K = (0 0 P D) .  (4.5)

The general solution of (4.3) for t > t_{0} can be written as

y (t) = e^{A(t−t_{0})}y (t_{0}) + ∫_{t_{0}}^{t} e^{A(t−s)}bQ (s) ds.  (4.6)

If in (4.6) we select t_{0} = jτ, t = (j + 1) τ, and introduce the notations y_{j} = y (jτ) , Q_{j} = Q (jτ) , j = 1, 2, …, we have

y_{j+1} = e^{Aτ} y_{j}+∫_{jτ}^{(j+1)τ} e^{A((j+1)τ−s)}bdsQ_{j} = e^{Aτ} y_{j}+∫_{0}^{τ} e^{A(τ−s)}bdsQ_{j}, j = 1, 2, …

(4.7)

Note that in (4.7) we have used that Q (t) is a piecewise constant function. Let

A_{d} = e^{Aτ} , b_{d} = ∫_{0}^{τ} e^{A(τ−s)}dsb.  (4.8)

The discrete-time model, i.e., t = jτ, j = 1, 2, …, consisting the state and feedback equations assumes the form

y_{j+1} =A_{d}y_{j} + b_{d}Q_{j}
Q_{j+1} =Ky_{j}

for j = 1, 2, …, or equivalently

z_{j+1} = Sz_{j}, j = 1, 2, …  (4.9)

where

z_{j} = (\begin{matrix}y_{j} \\ Q_{j} \end{matrix}) , S = (\begin{matrix}A_{d} & b_{d} \\ K & 0 \end{matrix}) .  (4.10)

Recall that the discrete system (4.9) is asymptotically stable if and only if all its characteristic multipliers μ_{i}, i = 1, …, 5 are in modulus less than one, i.e., |μ_{i}| < 1, i = 1, …, 5. The characteristic equation of (4.9) has the following form

g_{5}μ^{5} + g_{4}μ^{4} + g_{3}μ^{3} + g_{2}μ^{2} + g_{1}μ + g_{0} = 0 ,  (4.11)

where the coefficients g_{i}, i = 0, …, 5 are determined by the system parameters. Applying the Moebius-Zukovski transformation [3], i.e., substituting μ = (ν + 1) / (ν − 1) into (4.11) and multypling it by (ν − 1)^{5}, we obtain a fifth order equation for ν in the form of (4.11). Since |μ_{i}| < 1, i = 1, …, 5 if and only if ℜν_{i} < 0, i = 1, …, 5, the Routh-Hurwitz criterion can be used to determine the stability conditions following the same procedure as in Theorem 2.2.

The stability domain is bordered by the same line as in case of the system without time delay and a parabola which depends on the time delay, resulting in a shrinking stability domain for increasing time delays. The stability chart is given in Figure 4.1(a) for τ = 0 [ms] and τ = 5 [ms].

The stability domain disappears at a certain critical value of the time delay, so balancing is impossible above that value. The critical delay depends on the parameters describing the system. The most interesting is the dependence on the length of the pendulum and the stiffness of the belt. The connections between the critical time delay and the length of the pendulum and between the critical time delay and the stiffness of the belt are shown in Figure 4.2(a).

Figure 4.1: (a) The stability charts for τ = 0 [ms] and τ = 5 [ms] (b) The domain of 2 kinds of stable motions for τ = 5 [ms]

In the latter case we can see an asymptote which corresponds to the result gained for the system with rigid belt. As the length of the pendulum decreases, the critical time delay also decreases. There is a minimal length l_{min} of the pendulum, where the critical time delay becomes 0. It is proportional to 1/s. If the belt were ideally rigid, then l_{min} = 0, but since the belt is elastic, l_{min} > 0. Our result is shown in Figure 4.2(a), where s = 10 [kN/m] and l_{min} = 0.16 [mm]. Similarly, as the stiffness of the belt decreases, the critical time delay also decreases and at a certain minimal value s_{min}, it becomes 0. The minimal stiffness of the belt is proportional to 1/l. In case of Figure 4.2(a), l = 0.5 [m] and s_{min} = 3.14 [N/m].

Now we study the behavior of system (4.9) on the boundary of the stability region. A real characteristic multiplier crosses at the point (1, 0) of the complex plane as P approaches the linear segment of the boundary of the stability region, while a complex conjugate pair crosses the unit circle of the complex plane as P approaches the parabolic part of the boundary of the stability region. The characteristic multiplier, μ, of a delay equation is defined as μ = e^{λτ}, where λ is a root of the characteristic equation. The angular frequency ω of the oscillation at the loss of stability is the imaginary part of the characteristic root with zero real part. Since

λ = \frac{1}{τ} ln μ = \frac{1}{τ} (ln |μ| + i (Θ + 2kπ)) ,

where Θ is the angle of μ and k is an integer, then it follows that ℜλ = 0 if |μ| = 1 and

ω = ℑλ = \frac{1}{τ} Θ = \frac{1}{τ} arctan \frac{ℑμ}{ℜμ}

as we choose k = 0. Note that |μ| = 1 corresponds to ℜν = 0, so ν = iκ and μ = (iκ + 1) / (iκ − 1) = (1 − κ^{2}) / (−1 − κ^{2}) + 2κi/ (−1 − κ^{2}). We have ℑμ/ℜμ = 2κ/ (1 − κ^{2}) and therefore

ω = \frac{1}{τ} arctan \frac{2κ}{1 − κ^{2}} .  (4.12)

Figure 4.2: (a) The critical time delay vs. the length of the pendulum and the spring stiffness (b) The frequency of the oscillation at the loss of stability and the stability domain

Choosing control parameters from the boundary of the stability domain, substituting ν = iκ in the characteristic equation, and solving either its real or imaginary part for κ, the frequency of the oscillation is obtained dividing equation (4.12) by 2π. Since arctan 2κ/ (1 − κ^{2}) < π/2, we have the estimate

f < \frac{1}{4τ}  (4.13)

for the frequency of oscillation. This means that the frequency of the self-excited oscillation arising after the loss of stability can vary between zero and 25% of the sampling frequency 1/τ. The frequency f of the oscillation at the loss of stability is drawn in Figure 4.2(b).

If we have backlash in the system we use (3.2)-(3.3) instead of (2.18). The delayed control force (4.1) leads to the following systems of delay differential equations

(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} \end{matrix}) (\begin{matrix}Δ̈ & (t) \\ φ̈ & (t) \end{matrix})

+ (\begin{matrix}(m + M) \frac{K}{r_{m}} & 0 \\ 0 & 0 \end{matrix}) (\begin{matrix}Δ̇ & (t) \\ φ̇ & (t) \end{matrix})

+ (\begin{matrix}(m + M) r_{m}s + \frac{1}{2}m_{m}r_{m} \frac{r_{w}^{2}}{R_{w}^{2}} s & 0 \\ \frac{1}{2} \frac{m}{m+M} l \frac{r_{w}}{R_{w}} s & −\frac{1}{2}mgl \end{matrix}) (\begin{matrix}Δ (t) \\ φ (t) \end{matrix})

= (\begin{matrix}(m + M)(P φ((j − 1)τ) + Dφ̇((j − 1)τ)) \\ +(m + M + \frac{1}{2}m_{m} \frac{r_{w}^{2}}{R_{w}^{2}} )r_{m}sr_{0} sgn Δ(t) \\ \frac{1}{2} \frac{m}{m+M} l \frac{r_{w}}{R_{w}} sr_{0} sgn Δ (t) \end{matrix}) ,

(4.14)

for t ∈ [jτ, (j + 1)τ) and j = 1, 2, …, if |Δ| ≥ r_{0}; and

\begin{matrix}(\begin{matrix}\frac{1}{2} (m + M) m_{m}r_{m} & −\frac{1}{4}mm_{m}lr_{m} \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^{2} − \frac{1}{4} \frac{m}{(m+M)}ml^{2} \end{matrix}) & (\begin{matrix}Δ̈ & (t) \\ φ̈ & (t) \end{matrix}) \\ + (\begin{matrix}(m + M) \frac{K}{r_{m}} & 0 \\ 0 & 0 \end{matrix}) (\begin{matrix}Δ̇ & (t) \\ φ̇ & (t) \end{matrix}) + (\begin{matrix}0 & 0 \\ 0 & −\frac{1}{2}mgl \end{matrix}) (\begin{matrix}Δ (t) \\ φ (t) \end{matrix}) \\ = (\begin{matrix}(m + M) & (P φ ((j − 1) τ) + Dφ̇ ((j − 1) τ)) \\ & 0 \end{matrix}) , t ∈ [jτ, (j + 1) τ) \end{matrix} (4.15)

for j = 1, 2, …, if |Δ| < r_{0}. Transform the second order system (4.14) and (4.15) into a system of first order equations. We obtain

ẏ (t) =Ay (t) + bQ (t) + f
Q (t) =Ky ((j − 1) τ) ,

for t ∈ [jτ, (j + 1) τ), and j = 1, 2, …; and

ẏ (t) =A^{∗}y (t) + bQ (t)
Q (t) =Ky ((j − 1) τ) ,

for t ∈ [jτ, (j + 1) τ) and j = 1, 2, …, respectively, where y, A, b are given in (4.4) and K is given in (4.5), while f and A^{∗} assume the form

f = (0 (\frac{4r_{w}^{2}sr_{0}}{(m+4M)R_{w}^{2}} + \frac{2sr_{0}}{m_{m}} ) sgn Δ 0 \frac{6r_{w}sr_{0}}{(m+4M)lR_{w}} sgn Δ)^{T} ,

A^{∗} = (\begin{matrix}0 & 1 & 0 & 0 \\ 0 & −\frac{2K}{m_{m}r_{m}^{2}} & \frac{3mgr_{w}}{(m+4M)R_{w}} & 0 \\ 0 & 0 & 0 & 1 \\ 0 & 0 & \frac{6(m+M)g}{(m+4M)l} & 0 \end{matrix}) .

The discrete-time model, associated with (4.14) and (4.15), is introduced in the form

y_{j+1} =A_{d}y_{j} + b_{d}Q_{j} + f_{d}
Q_{j+1} =Ky_{j} (4.16)

Figure 4.3: Δ̇ − Δ phase-plane (a) P = 1 [Nm], D = 0.2 [Nms] (b) P = 10 [Nm], D = 0.5 [Nms]

for j = 1, 2, …, if |Δ| ≥ r_{0}; and

y_{j+1} =A_{d}^{∗}y_{j} + b_{d}Q_{j}
Q_{j+1} =Ky_{j} (4.17)

for j = 1, 2, …, if |Δ| < r_{0}, respectively, where A_{d} and b_{d} are given in (4.8), while the constant term f_{d} and the matrix A_{d}^{∗} is obtained as follows

f_{d} = ∫_{0}^{τ} e^{A(τ−s)}dsf, A_{d}^{∗} = e^{A^{∗}τ} .

System (4.16) and (4.17) can be rewritten in the form

z_{j+1} =Sz_{j} + f_{j} |Δ| ≥ r_{0}
z_{j+1} =S^{∗}z_{j} |Δ| < r_{0}, (4.18)

for j = 1, 2, …, where z_{j} and S is given in (4.10), while

S^{∗} = (\begin{matrix}A_{d}^{∗} & b_{d} \\ K & 0 \end{matrix}) , f_{j} = (\begin{matrix}f_{d} \\ 0 \end{matrix}) , j = 1, 2, …

Considering system (4.18), two kinds of stable motions are obtained using MATLAB simulations. The equilibria (Δ, Δ̇, φ, φ̇) = (±r_{0}, 0, 0, 0) are stable

left to the dotted line in Figure 4.1(b). A trajectory in the Δ̇ − Δ phaseplane is drawn in Figure 4.3(a) for P = 1 [Nm], D = 0.2 [Nms]. Since the phase-space is symmetric, the same kind of motion exists around the other equilibrium. A different type of stable motion occurs right to the dotted line in Figure 4.1(b). P = 10 [Nm], D = 0.5 [Nms] in Figure 4.3(b), where this type of motion is shown. Trajectories approach and penetrate a stable set, but since there exists no limit cycle inside the stable set, the motions appear irregular (i.e. quasiperiodic or even chaotic). Peaks can be seen along the trajectories of both motion. They appear at each sampling, when the control force changes.

Figure 4.4: Poincarè maps at the enter to the zone of backlash at Δ = r_{0} (a) The 3 dimensional Poincarè map (b), (c), (d) Series of 1 dimensional Poincarè maps

Note that the presented phase-diagrams are results of simulations, thus further studies could provide additional information about the type of motions obtained. In particular, Poincarè maps, Fourier spectra, behavior of nearby trajectories and Lyapunov exponents are investigated numerically with MATLAB in order to learn more about these motions. For the remainder of this section, we used fixed values of the control parameters, i.e., P = 10 [Nm] and D = 0.5 [Nms].

Poincarè maps are shown in Figure 4.4. The three dimensional space Δ = r_{0} = 1 [mm] intersects the four dimensional phase-space. For chaotic motions the Poincarè map shows irregularly spread points inside the attractor in the cross section. Figure 4.4(a) shows this three dimensional cross section. Points are spread, but some regularity can be explored. One dimensional Poincarè maps are also drawn, in Figure 4.4(b), 4.4(c) and 4.4(d). Points are situated in an attractor, but the longer we let the simulations run, the denser the filling of the attractor become.

The Fourier spectrum plotted in Figure 4.5(a) is obtained by applying the Fast Fourier Transformation. It shows the highest peak clearly, but many other peaks also appear. It implies that although the motion has a dominating frequency, many other frequencies arise, so the motion is more complicated than periodic.

Nearby trajectories of chaotic motions stretch each other. It means that Figure 4.5: (a) Fourier spectrum (b) The deviation between nearby trajectories

two trajectories with arbitrarily close initial conditions may move arbitrarily far from each other inside the attractor. Figure 4.5(b) draws the difference between the Δ̇ components of two trajectories in time. The initial difference is 0.001 [m/s]. We can state that trajectories neither approach nor move away from each other, their distance changes between 0 and an upper limit.

The Lyapunov exponents describe the stretching rate of trajectories in a direction. Negative Lyapunov exponent shows that trajectories approach each other in that direction, while positive Lyapunov exponent indicates that trajectories stretch each other in that direction. The Lyapunov exponents are determined by applying the following algorithm [1]. We choose an orthogonal basis {w_{1}^{0}, …, w_{5}^{0}} and compute the vectors v_{1}^{1}, …, v_{5}^{1} by multiplying the chosen basis by the matrix S or S^{∗} of (4.18) at the initial condition. Use the Gram-Schmidt orthogonalization procedure to orthogonalize these vectors and get w^{∗}_{1}^{1}, …, w^{∗}_{5}^{1}. To eliminate the problem of extremely large and small numbers, we normalize these vectors to obtain the orthogonal basis of the next step {w_{1}^{1}, …, w_{5}^{1}}. Repeating this step n times, we get the following estimate for the ith Lyapunov exponent

λ_{i} = \frac{ln ∥w^{∗}_{i}^{n}∥ + ⋯ + ln ∥w^{∗}_{i}^{1}∥}{n}, i = 1, …, 5.

The result of this algorithm shows that three of the Lyapunov exponents are negative, but two of them tend to 0, as n increases. This result predicts that the largest Lyapunov exponent is 0 which means that trajectories neither approach nor move away from each other. This shows good agreement with the discussion in the previous paragraph.

5 Experimental observations

A pendulum-cart structure was constructed [22] and experiments to validate our simulation results were carried out. Measurements were taken with different values of the control parameters P and D and with different time delays τ. In theoretical/simulation studies above we used the assumptions that position Figure 5.1: The theoretical stability domain and the measured points for τ = 5 [ms]

of the cart was arbitrary, but we had to impose restrictions on the position of the cart in the experimental setting. The proportional and the differential gain of the cart were different from 0, but they were fixed, and those of the pendulum were varied. Therefore, the stability domain in the plane of the control parameters P and D was different from that gained from the system investigated in this paper.

The time delay is τ = 5 [ms], the proportional and the differential gain of the cart are P_{x} = 7.7 [N] and D_{x} = 15.4 [Ns], respectively, in Figure 5.1. The continuous curve surrounds the analitically determined stability domain, while the small circles, squares, crosses and triangle represent measurement points. The circles indicate stable controllers. Since in reality the backlash is always positive and all the disturbances cannot be eliminated, the pendulum oscillates around its upper equilibrium with small amplitude. Figure 5.2(a) shows the time evolution of the angle of the pendulum for P = 100 [Nm], D = 6 [Nms]. The oscillation frequency is about 0.45 [Hz]. The squares in Figure 5.1 indicate controllers at the boundary of the stability domain, where oscillations cause loss of stability. We took measurements for P = 700 [Nm], D = 6 [Nms], the time evolution of the angle of the pendulum is depicted in Figure 5.2(b). The pendulum swings with higher frequency around its upper equilibrium. The oscillation frequency is about 12 [Hz] which is still less than 25 % of the sampling frequency 1/τ. It corresponds to the estimate given in (4.13). The crosses in Figure 5.1 represent unstable controllers. Figure 5.2(c) shows the time evolution of the voltage of the motor providing the control force for P = 100 [Nm], D = 2 [Nms]. This figure does not show increasing values of the voltage, because it cannot exceed its maximum value. The motor voltage reaches its maximum in a few tenth of a second since the controller is unstable and the pendulum loses its stability. The triangle in Figure 5.1 indicates a controller which works with quite large control force, since both the proportional and the differential gain is large, i.e., P = 700 [Nm], D = 23 [Nms]. Figure 5.2(d) shows that the Figure 5.2: Time evolutions for different controllers (a) P = 100 [Nm], D = 6 [Nms] (b) P = 700 [Nm], D = 6 [Nms] (c) P = 100 [Nm], D = 2 [Nms] (d) P = 700 [Nm], D = 23 [Nms]

voltage of the motor reaches its maximum value quite frequently which leads to oscillations and loss of stability. The experimental observations do not coincide exactly with the theoretical results, but they show good agreement. If the control force too small, then the upper equilibrium of the pendulum is unstable. There is a domain in the plane of the control parameters where the equilibrium is stable. Increasing the control force causes oscillations and loss of stability. The oscillation frequency arising at the loss of stability is less than 25 % of the sampling frequency.

The theoretical critical time delay is τ_{cr} = 88 [ms]. Stable motions is obtained for τ = 20 [ms], but controllers with delay of 30 [ms] do not stabilize the upper equilibrium of the pendulum. The stability domain is quite small if the time delay is greater than 20 [ms], therefore it is quite difficult to find it experimentally or it does not exist due to some unmodeled effects and disturbances.

The frequency of the motion presented in Figure 5.2(a) and 5.2(b) is about 0.45 [Hz] and 12 [Hz]. The largest peaks occur in the Fourier spectra at these values (see Figure 5.3(a) and 5.3(b)).

Conclusions

Dynamics of digitally controlled inverted pendulum on a cart is examined. It is easy to stabilize the upper equilibrium of the simplest model of the pendulumcart system. A stability domain has been found in the plane of the control parameters. Backlash occuring at the driving wheel introduces small ampli-Figure 5.3: Fourier spectra (a) P = 100 [Nm], D = 6 [Nms] (b) P = 700 [Nm], D = 6 [Nms]

tude periodic oscillations of the controlled pendulum around its upper equilibrium. The combination of backlash and time delays lead to irregular oscillations. (Our studies so far indicates the appearance of motions which are ”marginally” chaotic.) We plan to give a more detailed characterization of these irregularities in future projects.

References

[1] Alligood, K. T., Sauer, T. D., Yorke, J. A., Chaos (An introduction to dynamical systems), Springer-Verlag, New York, 1996.

[2] Enikov, E., Stépán, G., Micro-Chaotic Motion of Digitally Controlled Machines, J. of Vibration and Control, 4 pp. 427-443, 1998.

[3] Farkas, M., Periodic Motions, Springer-Verlag, New York, 1994.

[4] Hogan, S. J., On the dynamics of rigid block motion under harmonic forcing, Proc. R. Soc. Lond. A 425 pp. 441-476, 1989.

[5] Hogan, S. J., The many stady state responses of a rigid block under harmonic forcing, Earthquake Engineering and Structural Dynamics, Vol. 19, pp. 1057-1071, 1990.

[6] Kahraman, A., On the Response of a Preloaded Mechanical Oscillator with a Clearance: Period-Doubling and Chaos, Nonlinear Dynamics 3, pp. 183-198, 1992.

[7] Kawazoe, Y., Manual control and computer control of an inverted pendulum on a cart, Proc. 1st Int. Conf. on Motion and Vibration Control, pp. 930-935, Yokohama, 1992.

[8] Kollár, L. E., Backlash in Machines Stabilized by Control Force, Proc. of 1st Conference on Mechanical Engineering Vol. 1, pp. 147-151, Budapest, 1998.

[9] Kollár, L. E., Stépán, G., Digital Controlling of Piecewise Linear Systems, Proc. of 2nd Conference on Control of Oscillations and Chaos, Vol. 2, pp. 327-330, St.Petersburg, 2000.

[10] Kollár, L. E., Stépán, G., Hogan, S. J., Backlash in Balancing Systems Using Approximate Spring Characteristics, Proc. of 3rd European Nonlinear Oscillations Conference, 2000.

http://serv1.imm.dtu.dk/documents/users/mps/ENOC/proceedings/Kollar/),

[11] Kollár, L. E., Stépán, G., Hogan, S. J., Sampling Delay and Backlash in Balancing Systems, Periodica Polytechnica Ser. Mech. Eng., Vol. 44, No. 1, pp. 77-84, 2000.

[12] Leine, R. I., van Campen, D. H., Fold Bifurcations in Discontinuous Systems, Proc. of DETC’99, 1999 ASME Design Engineering Technical Conferences, Las Vegas, 1999.

[13] Lóránt, G., Stépán, G., The Role of Non-Linearities in the Dynamics of a Single Railway Wheelset, Machine Vibration 5, pp. 18-26, 1996.

[14] Mahfouz, I. A., Badrakhan, F., Chaotic Behaviour of Some Piecewise-Linear Systems, Part 1: Systems with Set-up Spring or with Unsymmetric Elasticity, Journal of Sound and Vibration, 143(2), pp. 255-288, 1990.

[15] Mahfouz, I. A., Badrakhan, F., Chaotic Behaviour of Some Piecewise-Linear Systems, Part 2: Systems with Clearance, Journal of Sound and Vibration, 143(2), pp. 289-328, 1990.

[16] Natsiavas, S., On the dynamics of oscillators with bi-linear damping and stiffness, Int. J. Non-Linear Mechanics, 25 pp. 535-554, 1990.

[17] Natsiavas, S., Dynamics of multiple-degree-of-freedom oscillators with colliding components, Journal of Sound and Vibration, 165(3), pp. 439-453, 1993.

[18] Nigm, M. M., Shabana, A. A., Effect of an Impact Damper on a Multidegree of Freedom System, Journal of Sound and Vibration, 89(4), pp. 541-557, 1983.

[19] Shaw, S. W., Holmes, P. J., A Periodically Forced Piecewise Linear Oscillator, Journal of Sound and Vibration, 90(1), pp. 129-155, 1983.

[20] Shaw, S. W., The Dynamics of a Harmonically Excited System Having Rigid Amplitude Constraints, Part 1: Subharmonic Motions and Local Bifurcations, ASME Journal of Applied Mechanics, 52, pp. 453-458, 1985.

[21] Shaw, S. W., The Dynamics of a Harmonically Excited System Having Rigid Amplitude Constraints, Part 2: Chaotic Motions and Global Bifurcations, ASME Journal of Applied Mechanics 52, pp. 459-464, 1985.

[22] Szászi, I., Gáspár, P., Robust servo control design using the H_{∞}/μ method, Periodica Polytechnica Ser. Mech. Eng., Vol. 27, No. 1-2, pp. 3-16, 1999.

[23] Theodossiades, S., Natsiavas, S., Nonlinear Dynamics of Gear-Pair Systems with Periodic Stiffness and Backlash, Journal of Sound and Vibration, Vol. 229, pp. 287-310, 2000.

László E. Kollár Department of Applied Sciences University of Québec at Chicoutimi 555, Boul. de l’Université, Chicoutimi G7H 2B1 Québec, Canada E-mail address: laszlo kollar@uqac.ca Gábor Stépán Department of Applied Mechanics Budapest University of Technology and Economics H-1521 Budapest, Hungary E-mail address: stepan@mm.bme.hu János Turi Department of Mathematical Sciences The University of Texas at Dallas P.O. Box 830688, MS EC 35, Richardson, TX 75083-0688, USA E-mail address: turi@utdallas.edu