\documentclass[twoside]{article} \usepackage{amsfonts,amsmath} % font used for R in Real numbers \usepackage{graphicx} % for including postscript files \pagestyle{myheadings} \markboth{\hfil Dynamics of delayed piecewise linear systems \hfil EJDE/Conf/10} {EJDE/Conf/10 \hfil L. E. Koll\'ar, G. St\'ep\'an, \& J. Turi \hfil} \begin{document} \setcounter{page}{163} \title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent Fifth Mississippi State Conference on Differential Equations and Computational Simulations, \newline Electronic Journal of Differential Equations, Conference 10, 2003, pp 163--185. \newline http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu (login: ftp)} \vspace{\bigskipamount} \\ % Dynamics of delayed piecewise linear systems % \thanks{ {\em Mathematics Subject Classifications:} 34K35, 37C75. \hfil\break\indent {\em Key words:} Stability analysis, backlash, digitally controlled pendulum. \hfil\break\indent \copyright 2003 Southwest Texas State University. \hfil\break\indent Published February 28, 2003. } } \date{} \author{L\'aszl\'o E. Koll\'ar, G\'abor St\'ep\'an, \& J\'anos Turi} \maketitle \begin{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. \end{abstract} \newtheorem{theorem}{Theorem}[section] \numberwithin{equation}{section} \numberwithin{figure}{section} \section{Introduction} Dynamical systems with piecewise linear components in the equations of motion occur frequently in practice. Gear pairs with backlash \cite{Theo}, railway wheelsets with clearance \cite{Lorant}, mechanical oscillators with clearance \cite{Kahraman,Mahfouz1,Mahfouz2, ShawHolmes} or amplitude constraints \cite{Shaw1,Shaw2}, impact dampers \cite{Nats90,Nats93,Nigm}, moving parts with dry friction \cite{vanCampen} and adjacent structures during earthquake \cite{Hogan89,Hogan90} 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 \cite{Enikovmu,Kawazoe}. 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. 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. \section{The inverted pendulum on a cart} The model of the inverted pendulum on a cart can be seen in Figure \ref{pencart}, \cite{Gep98,COC2000,Lyngby,PerPol}. 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 \( \varphi \) of the pendulum and the angle \( \psi \) of the motor axle, respectively. It is assumed that the angle \( \varphi \) 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 \begin{equation}\label{pd} U_{m} = P_{1} x+D_{1} \dot x+P_{2} \varphi +D_{2} \dot \varphi \,, \end{equation} where the constants \( P_1, D_1, P_2, D_2 \) are the so-called PD controller parameters. \( U_m \) together with \( \dot \psi \), the angular velocity of the motor determine \( M_m \), the driving torque of the motor by \begin{equation}\label{momment} M_{m} = L U_{m}-K \dot \psi \,, \end{equation} where \( L \) and \( K \) are given motor parameters. \begin{figure}[t] \begin{center} \includegraphics[width=0.9\textwidth]{cart.eps} \end{center} \caption{\sl The inverted pendulum on a cart and its stability map} \label{pencart} \end{figure} The nonlinear equations of motion are obtained (see also Figure \ref{pencart}) by means of the Lagrange equations of the second kind. The kinetic energy \( \mathcal{T} \), the potential \( \mathcal{U} \), the dissipation \( \mathcal{D} \) and the general force \( Q \) have the following form \begin{equation} \label{gf} \begin{gathered} \mathcal{T}=\frac{1}{2}m \left( \left( \dot x+\frac{l}{2} \dot \varphi \cos \varphi \right)^2+\left( \frac{l}{2} \dot \varphi \sin \varphi \right)^2 \right)+\frac{1}{24}ml^2 \dot \varphi^2+ \frac{1}{2}M \dot x^2+\frac{1}{2}J_m \dot \psi^2,\\ \mathcal{U}=-\frac{1}{2}mgl \left( 1-\cos \varphi \right)+U_s\,,\\ \mathcal{D}=0,\quad Q_1^e=0, \quad Q_2^e=0, \quad Q_3^e=M_m\,, \end{gathered} \end{equation} where the index \( e \) refers to elastic, \( U_s=0 \) if the belt is ideally rigid, \( U_s=s \Delta^2/2 \) if the belt is elastic and \( \Delta \) is the elongation of the spring, \( M=m_m+m_c+4 m_w+4 J_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 \( 4 m_w+ 4 J_w/R_w^2 \). Applying the Lagrange equations of the second kind \begin{equation}\label{lagrange} \frac{\mathrm{d}}{\mathrm{d}t} \frac{\partial \mathcal{T}}{\partial \dot q_k} -\frac{\partial \mathcal{T}}{\partial q_k}+\frac{\partial \mathcal{U}}{\partial q_k}+\frac{\partial \mathcal{D}}{\partial \dot q_k}=Q_k\,, \end{equation} \( k=1,2,3 \), the nonlinear equations of motion for the controlled pendulum are obtained \begin{equation}\label{elbeltnonlin} \begin{aligned} \begin{pmatrix} m+M & \frac{1}{2} ml \cos \varphi & 0 \\ \frac{1}{2} ml \cos \varphi & \frac{1}{3} ml ^2 & 0 \\ 0 & 0 & \frac{1}{2} m_m r_m^2 \end{pmatrix} \begin{pmatrix} \ddot x \\ \ddot \varphi \\ \ddot \psi \end{pmatrix} & \\ +\begin{pmatrix} s \frac{r_w^2}{R_w^2} & 0 & -s r_m \frac{r_w}{R_w} \\ 0 & 0 & 0 \\ -s r_m \frac{r_w}{R_w} & 0 & s r_m^2 \end{pmatrix} \begin{pmatrix} x \\ \varphi \\ \psi \end{pmatrix}+ \begin{pmatrix} -\frac{1}{2}ml \dot \varphi^2 \sin \varphi \\ -\frac{1}{2}mgl \sin \varphi \\ 0 \end{pmatrix}&= \begin{pmatrix} 0 \\ 0 \\ Q_3^e \end{pmatrix}. \end{aligned} \end{equation} \subsection*{The case when the belt is ideally rigid} If the belt is ideally rigid \cite{Gep98}, then \( x \) and \( \psi \) are related by the equation \( \psi=x r_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 \( \psi \) can be eliminated from equation (\ref{elbeltnonlin}). Multiplying the third equation of (\ref{elbeltnonlin}) by \( r_m r_w/R_w \), adding it to the first equation and substituting \( \psi=x r_w/(r_m R_w) \) yields the nonlinear equations of motion \begin{equation}\label{eqrigid} \begin{pmatrix} m_r & \frac{1}{2} ml \cos \varphi\\ \frac{1}{2} ml \cos \varphi & \frac{1}{3} ml ^2 \end{pmatrix} \begin{pmatrix} \ddot x \\ \ddot \varphi \end{pmatrix}- \begin{pmatrix} \frac{1}{2}ml \dot \varphi^2 \sin \varphi \\ \frac{1}{2}mgl \sin \varphi \end{pmatrix}= \begin{pmatrix} Q_1^r \\ 0 \end{pmatrix}, \end{equation} where \( m_r=m+M+m_m r_w^2/(2 R_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 \( \varphi=0 \) we obtain the variational system \begin{equation}\label{eqrigidlin} \begin{pmatrix} m_r & \frac{1}{2} ml \\ \frac{1}{2} ml & \frac{1}{3} ml ^2 \end{pmatrix} \begin{pmatrix} \ddot x \\ \ddot \varphi \end{pmatrix}- \begin{pmatrix} 0 & 0 \\ 0 & \frac{1}{2}mgl \end{pmatrix} \begin{pmatrix} x \\ \varphi \end{pmatrix}= \begin{pmatrix} Q_1^r \\ 0 \end{pmatrix}. \end{equation} Using equations (\ref{pd}), (\ref{momment}) and (\ref{gf}), the control force has the form \begin{equation}\label{q1} Q_1^r= \frac{r_{w}}{r_{m} R_{w}} L \left( P_1 x+D_1 \dot x+ P_2 \varphi +D_2 \dot \varphi \right)- K \frac{r_{w} ^2}{r_{m} ^2 R_{w} ^2} \dot x \,. \end{equation} It is easy to see from (\ref{q1}) that if \begin{equation}\label{pxdx} P_1=0 \quad {\rm and} \quad D_1=\frac{1}{L}\frac{r_w}{r_m R_w} K \,, \end{equation} then \( Q_1^r \) can be written as \begin{equation}\label{cfrigid} Q_1^r= \frac{r_{w}}{r_{m} R_{w}} \left( P \varphi + D \dot \varphi \right) \,, \end{equation} where \( P=L P_2 \) and \( D=L D_2 \). Thus, the control force \( Q_1^r \) is only a function of \( \varphi \) and \( \dot \varphi \) (i.e., \( x \) and \( \dot x \) 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 (\ref{eqrigidlin}) for \( \ddot x \) and substitute the given expression into the second equation of (\ref{eqrigidlin}). A one degree of freedom system is obtained with the following single second order governing equation \begin{equation}\label{onedof} \ddot \varphi-\frac{6g}{l} \varphi=-\frac{6}{m_r l} Q_1^r \,. \end{equation} In our setting stabilization means that \( \varphi \) and \( \dot \varphi \) are both zero, but \( x \) and \( \dot x \) could be nonzero. Due to the relationship between \( \ddot x \), \( \varphi, \dot \varphi \) and \( \ddot \varphi \) as described by equation (\ref{eqrigid}), we get that \( \dot x \) is constant when \( \varphi=\dot \varphi=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 (\ref{cfrigid}). \begin{theorem} \label{thm2.1} The trivial solution of (\ref{cfrigid})--(\ref{onedof}) is asymptotically stable if and only if \begin{equation}\label{stcondrigid} P>P_{0}=m_r g r_m \frac{R_w}{r_w} \quad and \quad D>0 \,. \end{equation} \end{theorem} \paragraph{Proof:} The characteristic equation of (\ref{onedof}) has the following form \begin{equation}\label{chareqrigid} \lambda^2+\frac{6 r_w}{m_r l r_m R_w} D \lambda+ \left( \frac{6 r_w}{m_r l r_m R_w} P-\frac{6g}{l} \right)=0\,. \end{equation} The coefficients in (\ref{chareqrigid}) are clearly positive if the conditions (\ref{stcondrigid}) are satisfied. The statement of the theorem follows. \hfill\( \square \) \subsection*{The case when the belt is elastic} If the belt is elastic, then the nonlinear equations of motion assume the form of (\ref{elbeltnonlin}). This system can be reduced to a two degrees of freedom system if the assumptions (\ref{pxdx}) are applied again and a new general coordinate, \( \Delta \), the elongation of the spring is introduced by \begin{equation}\label{delta} \Delta=r_{m} \psi- \frac{r_{w}}{R_{w}} x \,. \end{equation} Multiplying the first and third equation of (\ref{elbeltnonlin}) by \( \left( -m_m r_m r_w \right)/(2 R_w) \) and \( \left( m+M \right) \), respectively, summing them and using (\ref{delta}) gives the first equation of motion of the reduced system. Solving the first equation of (\ref{elbeltnonlin}) for \( \ddot x \), using (\ref{delta}) and substituting it into the second equation of (\ref{elbeltnonlin}) gives the second equation of motion of the reduced system and altogether we have \begin{equation} \label{nonlinpencart} \begin{aligned} &\begin{pmatrix} \frac{1}{2}\left( m+M \right) m_{m} r_{m} & -\frac{1}{4}m m_{m} l r_m \frac{r_{w}}{R_{w}} \cos \varphi \\ 0 & \frac{1}{3}ml^2-\frac{1}{4}\frac{m}{\left( m+M \right)}ml^2 \cos^2 \varphi \end{pmatrix} \begin{pmatrix} \ddot \Delta \\ \ddot \varphi \end{pmatrix}\\ &+\begin{pmatrix} (m+M) \frac{K}{r_{m}} & 0 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \dot \Delta \\ \dot \varphi \end{pmatrix} + \begin{pmatrix} \frac{1}{4}m m_{m} l r_m \frac{r_{w}}{R_{w}} \dot \varphi^2 \sin \varphi \\ \frac{1}{4}\frac{m}{\left( m+M \right)} ml^2 \dot \varphi^2 \sin \varphi \cos \varphi-\frac{1}{2}mgl \sin \varphi \end{pmatrix}\\ &+ \begin{pmatrix} \left( m+M \right) r_{m}+ \frac{1}{2}m_{m} r_{m} \frac{r_{w}^2}{R_{w} ^2}\\ \frac{1}{2}\frac{m}{\left( m+M \right)}l \frac{r_w}{R_w} \end{pmatrix} R_{s}= \begin{pmatrix} \left( m+M \right) Q \\ 0 \end{pmatrix}\,, \end{aligned} \end{equation} where \begin{equation}\label{forcespring} R_{s}=\begin{pmatrix} -s\frac{r_{w}}{R_{w}} & 0 & sr_{m} \end{pmatrix} \begin{pmatrix} x \\ \varphi \\ \psi \end{pmatrix}=s \Delta \end{equation} is the force in the driving belt and \begin{equation}\label{cfpencart} Q=P \varphi+D \dot \varphi \end{equation} is the control force. Linearizing these equations around the upper equilibrium \( \varphi=0 \) yields \begin{equation}\label{linpencart} \begin{aligned} \begin{pmatrix} \frac{1}{2}\left( m+M \right) m_{m} r_{m} & -\frac{1}{4}m m_{m} l r_m \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^2-\frac{1}{4}\frac{m}{\left( m+M \right)}ml^2 \end{pmatrix} \begin{pmatrix} \ddot \Delta \\ \ddot \varphi \end{pmatrix}&\\ +\begin{pmatrix} (m+M) \frac{K}{r_{m}} & -\left( m+M \right) D \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \dot \Delta \\ \dot \varphi \end{pmatrix}&\\ +\begin{pmatrix} \left( m+M \right) r_{m} s+ \frac{1}{2}m_{m} r_{m} \frac{r_{w}^2}{R_{w} ^2} s & -\left( m+M \right) P \\ \frac{1}{2}\frac{m}{\left( m+M \right)}l \frac{r_w}{R_w} s & -\frac{1}{2}mgl \end{pmatrix} \begin{pmatrix} \Delta \\ \varphi \end{pmatrix}&= \begin{pmatrix} 0 \\ 0 \end{pmatrix}\,. \end{aligned} \end{equation} The subsequent theorem establishes the stability properties of system (\ref{linpencart}). \begin{theorem}\label{pencartstthm} The trivial solution of system (\ref{linpencart}) is asymptotically stable if and only if \begin{equation}\label{pencartstconds} P>P_{0}=m_r g r_m \frac{R_w}{r_w}\,, \quad s>s_{cr}=\frac{3 \left( m+M \right) m_{m} g}{\left( m+4M+2 m_{m} r_{w} ^2/R_{w} ^2 \right) l}\,,\quad q \left( P,D,s \right)<0 \,, \end{equation} where $q \left( P,D,s \right)=a_{22}D^2+a_1 P+a_2 D+a_0$ is a parabola, with the coefficients \begin{align*} &a_{22} = 6m_m r_w R_w r_m^3 s \,, \\ &a_1 = 2 \left( m+4M \right) l R_w^2 K^2 \,, \\ &a_2 = -6 \left( m+M \right) m_m gR_w^2 r_m^2 K -2 \left( m+4M \right) l R_w^2 r_m^2 Ks -4 m_m l r_w^2 r_m^2 Ks \,, \\ &a_0 = 3mm_m glr_w R_w r_m K^2 \,. \end{align*} \end{theorem} \paragraph{Proof:} The characteristic equation of (\ref{linpencart}) can be written as \begin{equation}\label{chareq} b_4 \lambda^4+b_3 \lambda^3+b_2 \lambda^2+b_1 \lambda+b_0=0\,, \end{equation} where \begin{align*} b_4=&\frac{1}{24} \left( m+4M \right) m_m l r_m^2, \\ b_3=&\frac{1}{12} \left( m+4M \right) lK, \\ b_2=&\frac{1}{12} \left( m+4M+2m_m \frac{r_w^2}{R_w^2} \right) l r_m^2 s-\frac{1}{4} \left( m+M \right) m_m g r_m^2, \\ b_1=&\frac{1}{2} r_m \frac{r_w}{R_w} sD-\frac{1}{2} \left( m+M \right) gK, \\ b_0=&\frac{1}{2} r_m \frac{r_w}{R_w} sP-\frac{1}{2} \left( m+M+ \frac{1}{2} m_m \frac{r_w^2}{R_w^2} \right) gr_m^2 s\,. \end{align*} The Hurwitz matrix assumes the form $$ \begin{pmatrix} 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{pmatrix}. $$ According to the Routh-Hurwitz criterion \cite{Farkas}, 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 \left( P,D,s \right)<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 \). \( \square \) 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 \ref{pencartstthm}, and then we present a numerical bifurcation analysis for the original nonlinear system (\ref{nonlinpencart}). The stability charts in the plane of the control parameters are shown in Figure \ref{pencart} 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 \left( P,D,s \right)=0 \). The angular frequency, \( \omega, \) 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 \( \lambda=\mathrm{i} \omega \) in the characteristic equation (\ref{chareq}) and solving either its real or imaginary part for \( \omega \). Using the fact that \( q \left( P,D,s \right)=0 \) we can obtain the frequency of the oscillation as a function of \( P \) or \( D \) only. For example, \( \omega \) as a function of \( P \) can be obtained as the positive solution of the equation \begin{align*} \frac{1}{24} \left( m+4M \right) m_m l r_m^2 \omega^4&\\ -\left( \frac{1}{12} \left( m+4M+2m_m \frac{r_w^2}{R_w^2} \right) l r_m^2 s-\frac{1}{4} \left( m+M \right) m_m g r_m^2 \right) \omega^2&\\ +\frac{1}{2} r_m \frac{r_w}{R_w} sP-\frac{1}{2} \left( m+M+ \frac{1}{2} m_m \frac{r_w^2}{R_w^2} \right) g r_m^2 s&=0\,. \end{align*} The frequency \( f=\omega/\left( 2 \pi \right) \) of the oscillation at the loss of stability and the stability domain are drawn in Figure \ref{biflin}(a) for the given values of parameters: \( m=0.169 \) [kg], \( M=1.136 \) [kg], \( m_m=0.2 \) [kg], \( g=9.81 \; [{\rm m}/{\rm 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 \ref{biflin}(a) shows that there is a unique frequency for every pair \( P,D \) on the parabola. \begin{figure}[t] \begin{center} \includegraphics[width=0.45\textwidth]{pcartf.eps} \includegraphics[width=0.45\textwidth]{linp.eps}\\ \includegraphics[width=0.45\textwidth]{pcartst.eps} \includegraphics[width=0.45\textwidth]{lind.eps} \end{center} \caption{\sl (a) The frequency of the oscillation and the stability domain \enskip (b) Bifurcation diagrams for \( D=2 \) {\rm [Nms]} (upper) and for \( P=20 \) {\rm [Nm]} (lower)} \label{biflin} \end{figure} To study the relationship between (\ref{linpencart}) and (\ref{nonlinpencart}) 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 \ref{biflin}(b) for the values of parameters given in the previous paragraph. The bifurcation parameter is \( P \) in the upper half of Figure \ref{biflin}(b), and the control parameter \( D \) is kept at 2 [Nms]. A pitchfork bifurcation is occurred at \( P_{0}=0.1986 \) [Nm] obtained from (\ref{pencartstconds}), where the upper equilibrium of the pendulum becomes stable. Figure \ref{biflin}(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 \left( P,2,10000 \right)=0 \). An unstable stationary solution appears at the pitchfork bifurcation and \( \varphi \) tends to \( \pi/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 \ref{biflin}(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 \left( 20,D,10000 \right)=0 \). Periodic solution exists at the border of the stability domain where two complex conjugate characteristic roots are crossing the imaginary axis. If \( \varphi \) 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 \( \varphi \) is not small. Thus, the linear variational system (\ref{linpencart}) can be considered a good approximation of the original nonlinear system (\ref{nonlinpencart}) in the stability region given by (\ref{pencartstconds}). \section{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 \begin{equation}\label{forcespringbl} R_{s}=\begin{cases} s\left( \Delta+r_{0}\right) & \mbox{if } \Delta\leq -r_{0}\\ 0 & \mbox{if }\left| \Delta\right| 0, \) and processing, \( \tau_p>0, \) delays are included in the model of the inverted pendulum on a cart \cite{COC2000}, \cite{PerPol}. First, we consider the case of no backlash. Then the linearized equations of motion are given by (\ref{linpencart}) as before, but now delayed arguments will appear due to the fact that the control force, \( Q \left( t \right) \), takes the form $$ Q \left( t \right)=P \varphi \left( \left( j-1 \right) \tau_s \right)+ D \dot \varphi \left( \left( j-1 \right) \tau_s \right), \quad t \in \left[ \left( j-1 \right) \tau_s+\tau_p, \, j \tau_s+\tau_p \right), $$ for $j=1,2,\dots$ We shall assume that the sampling and the processing delays are the same, i.e., \( \tau_s=\tau_p=\tau \). The more general case, \( \tau_s \not = \tau_p, \) could be considered similarly. With that, the control force can be written as \begin{equation}\label{cfblsd} Q \left( t \right)=P \varphi \left( \left( j-1 \right) \tau \right)+ D \dot \varphi \left( \left( j-1 \right) \tau \right), \quad t \in \left[ j \tau, \left( j+1 \right) \tau \right), \quad j=1,2,\dots \end{equation} Without time delay we had a system (\ref{linpencart}) of ordinary differential equations. With control force (\ref{cfblsd}) we obtain the following system of delay differential equations \begin{equation}\label{linpencartsd} \begin{aligned} &\begin{pmatrix} \frac{1}{2}\left( m+M \right) m_{m} r_{m} & -\frac{1}{4}m m_{m} l r_m \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^2-\frac{1}{4}\frac{m}{\left( m+M \right)}ml^2 \end{pmatrix} \begin{pmatrix} \ddot \Delta \left( t \right) \\ \ddot \varphi \left( t \right) \end{pmatrix}\\ &+\begin{pmatrix} (m+M) \frac{K}{r_{m}} & 0 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \dot \Delta \left( t \right) \\ \dot \varphi \left( t \right) \end{pmatrix}\\ &+\begin{pmatrix} \left( m+M \right) r_{m} s+ \frac{1}{2}m_{m} r_{m} \frac{r_{w}^2}{R_{w} ^2} s & 0 \\ \frac{1}{2}\frac{m}{\left( m+M \right)}l \frac{r_w}{R_w} s & -\frac{1}{2}mgl \end{pmatrix} \begin{pmatrix} \Delta \left( t \right) \\ \varphi \left( t \right) \end{pmatrix}\\ &=\begin{pmatrix} \left( m+M \right) \left( P \varphi \left( \left( j-1 \right) \tau \right)+ D \dot \varphi \left( \left( j-1 \right) \tau \right) \right) \\ 0 \end{pmatrix}\,, \quad t \in \left[ j \tau, \left( j+1 \right) \tau \right), \end{aligned} \end{equation} for $j=1,2,\dots$. Transform the second order system (\ref{linpencartsd}) into a system of first order equations, i.e., multiply (\ref{linpencartsd}) by the inverse of the coefficient matrix of the second order term and introduce a new state vector \( {\bf y} \). We obtain \begin{equation}\label{cauchy} {\bf \dot y} \left( t \right)={\bf Ay} \left( t \right)+ {\bf b}Q \left( t \right), \end{equation} where \begin{equation}\label{cauchyelbelt} \begin{gathered} {\bf A}=\begin{pmatrix} 0 & 1 & 0 & 0 \\ -\frac{4 r_w^2 s}{\left( m+4M \right) R_w^2}-\frac{2s}{m_m} & -\frac{2K}{m_m r_m^2} & \frac{3mgr_w}{\left( m+4M \right) R_w} & 0 \\ 0 & 0 & 0 & 1 \\ -\frac{6r_w s}{\left( m+4M \right) lR_w} & 0 & \frac{6 \left( m+M \right) g}{\left( m+4M \right) l} & 0 \\ \end{pmatrix}, \\ {\bf b}=\begin{pmatrix} 0 \\ \frac{2}{m_m r_m} \\ 0 \\ 0 \end{pmatrix},\quad {\bf y}=\begin{pmatrix} \Delta \\ \dot \Delta \\ \varphi \\ \dot \varphi \end{pmatrix}. \end{gathered} \end{equation} Note that \( Q \left( t \right) \) can be expressed as $$ Q \left( t \right)={\bf Ky} \left( \left( j-1 \right) \tau \right), \quad t \in \left[ j \tau, \left( j+1 \right) \tau \right), \quad j=1,2,\dots \nonumber $$ where \begin{equation}\label{cauchyk} {\bf K}=\begin{pmatrix} 0 & 0 & P & D \end{pmatrix}. \end{equation} The general solution of (\ref{cauchy}) for \( t>t_0 \) can be written as \begin{equation}\label{gensol} {\bf y} \left( t \right)=\mathrm{e}^{{\bf A} \left( t-t_0 \right)} {\bf y} \left( t_0 \right)+\int_{t_0}^t \mathrm{e}^ {{\bf A} \left( t-s \right)} {\bf b} Q \left( s \right) \mathrm{d}s. \end{equation} If in (\ref{gensol}) we select \( t_0=j \tau, \ t=\left( j+1 \right) \tau, \) and introduce the notations \( {\bf y}_j={\bf y} \left( j \tau \right), \ Q_j=Q \left( j \tau \right), \ j=1,2,\dots, \) we have \begin{equation}\label{disceq} {\bf y}_{j+1}=\mathrm{e}^{{\bf A} \tau} {\bf y}_j+ \int_{j \tau}^{\left( j+1 \right) \tau} \mathrm{e}^{{\bf A} \left( \left( j+1 \right) \tau-s \right)} {\bf b} \mathrm{d}s Q_j=\mathrm{e}^{{\bf A} \tau} {\bf y}_j+ \int_0^{\tau} \mathrm{e}^{{\bf A} \left( \tau-s \right)} {\bf b} \mathrm{d}s Q_j, \quad j=1,2,\dots \end{equation} Note that in (\ref{disceq}) we have used that \( Q \left( t \right) \) is a piecewise constant function. Let \begin{equation}\label{coefblsd} {\bf A_d}=\mathrm{e}^{{\bf A} \tau}, \quad {\bf b_d}=\int_0^{\tau} \mathrm{e}^{{\bf A} \left(\tau-s \right)} \mathrm{d}s{\bf b}. \end{equation} The discrete-time model, i.e., \( t=j \tau, \ j=1,2,\dots, \) consisting the state and feedback equations assumes the form \begin{align*} {\bf y}_{j+1}=&{\bf A_d y}_j+{\bf b_d} Q_j \\ Q_{j+1}=&{\bf K y}_j \end{align*} for $j=1,2,\dots$, or equivalently \begin{equation}\label{itsd} {\bf z}_{j+1}={\bf S} {\bf z}_j, \quad j=1,2,\dots \end{equation} where \begin{equation}\label{zks} {\bf z}_j=\begin{pmatrix} {\bf y}_j \\ Q_j \end{pmatrix}, \quad {\bf S}=\begin{pmatrix} {\bf A_d} & {\bf b_d} \\ {\bf K} & 0 \end{pmatrix}. \end{equation} Recall that the discrete system (\ref{itsd}) is asymptotically stable if and only if all its characteristic multipliers \( \mu_i, \ i=1, \dots, 5 \) are in modulus less than one, i.e., \( \left| \mu_i \right|<1, \ i=1, \dots, 5 \). The characteristic equation of (\ref{itsd}) has the following form \begin{equation}\label{charblsd} g_5 \mu^5+g_4 \mu^4+g_3 \mu^3+g_2 \mu^2+g_1 \mu+g_0=0\,, \end{equation} where the coefficients \( g_i, \ i=0, \dots, 5 \) are determined by the system parameters. Applying the Moebius-Zukovski transformation \cite{Farkas}, i.e., substituting \( \mu=\left( \nu+1 \right)/\left( \nu-1 \right) \) into (\ref{charblsd}) and multypling it by \( \left( \nu-1 \right)^5 \), we obtain a fifth order equation for \( \nu \) in the form of (\ref{charblsd}). Since \( \left| \mu_i \right|<1, \ i=1, \dots, 5 \) if and only if \( \Re \nu_i<0, \ i=1, \dots, 5 \), the Routh-Hurwitz criterion can be used to determine the stability conditions following the same procedure as in Theorem \ref{pencartstthm}. \begin{figure}[t] \begin{center} \includegraphics[width=0.45\textwidth]{stdomtau.eps} \includegraphics[width=0.45\textwidth]{stcurves.eps} \end{center} \caption{\sl (a) The stability charts for \( \tau=0 \) {\rm [ms]} and \( \tau=5 \) {\rm [ms]} \enskip (b) The domain of 2 kinds of stable motions for \( \tau=5 \) {\rm [ms]}} \label{stabtaupencart} \end{figure} 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 \ref{stabtaupencart}(a) for \( \tau=0 \) [ms] and \( \tau=5 \) [ms]. \begin{figure}[t] \begin{center} \includegraphics[width=0.4\textwidth]{tau-l.eps} \includegraphics[width=0.4\textwidth]{ofrp.eps}\\ \includegraphics[width=0.4\textwidth]{tau-s.eps} \includegraphics[width=0.4\textwidth]{stcurv2.eps} \end{center} \caption{\sl (a) The critical time delay vs. the length of the pendulum and the spring stiffness \enskip (b) The frequency of the oscillation at the loss of stability and the stability domain} \label{tauls} \end{figure} 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 \ref{tauls}(a). 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 \ref{tauls}(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 \ref{tauls}(a), \( l=0.5 \) [m] and \( s_{min}=3.14 \) [N/m]. Now we study the behavior of system (\ref{itsd}) on the boundary of the stability region. A real characteristic multiplier crosses at the point \( \left( 1,0 \right) \) 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, \( \mu, \) of a delay equation is defined as \( \mu=\mathrm{e}^{\lambda \tau} \), where \( \lambda \) is a root of the characteristic equation. The angular frequency \( \omega \) of the oscillation at the loss of stability is the imaginary part of the characteristic root with zero real part. Since $$ \lambda=\frac{1}{\tau} \ln \mu=\frac{1}{\tau} \left( \ln \left| \mu \right|+i \left( \Theta+2 k \pi \right) \right), $$ where \( \Theta \) is the angle of \( \mu \) and \( k \) is an integer, then it follows that \( \Re \lambda=0 \) if \( \left| \mu \right|=1 \) and $$ \omega=\Im \lambda=\frac{1}{\tau} \Theta=\frac{1}{\tau} \arctan \frac{\Im \mu}{\Re \mu} $$ as we choose \( k=0 \). Note that \( \left| \mu \right|=1 \) corresponds to \( \Re \nu=0 \), so \( \nu=i \kappa \) and \( \mu=\left( i \kappa+1 \right)/\left( i \kappa-1 \right)=\left( 1-\kappa^2 \right)/\left( -1-\kappa^2 \right)+2 \kappa i/\left( -1-\kappa^2 \right) \). We have \( \Im \mu/\Re \mu=2 \kappa/ \left( 1-\kappa^2 \right) \) and therefore \begin{equation}\label{ofrsd} \omega=\frac{1}{\tau} \arctan \frac{2 \kappa}{1-\kappa^2}. \end{equation} Choosing control parameters from the boundary of the stability domain, substituting \( \nu=i \kappa \) in the characteristic equation, and solving either its real or imaginary part for \( \kappa \), the frequency of the oscillation is obtained dividing equation (\ref{ofrsd}) by \( 2 \pi \). Since \( \arctan 2 \kappa/\left( 1-\kappa^2 \right)<\pi/2 \), we have the estimate \begin{equation}\label{ofreqsd} f<\frac{1}{4 \tau} \end{equation} 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/ \tau \). The frequency \( f \) of the oscillation at the loss of stability is drawn in Figure \ref{tauls}(b). If we have backlash in the system we use (\ref{eqblout})-(\ref{eqblin}) instead of (\ref{linpencart}). The delayed control force (\ref{cfblsd}) leads to the following systems of delay differential equations \begin{equation}\label{eqblsdout} \begin{aligned} &\begin{pmatrix} \frac{1}{2}\left( m+M \right) m_{m} r_{m} & -\frac{1}{4}m m_{m} l r_m \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^2-\frac{1}{4}\frac{m}{\left( m+M \right)}ml^2 \end{pmatrix} \begin{pmatrix} \ddot \Delta \left( t \right) \\ \ddot \varphi \left( t \right) \end{pmatrix}\\ &+\begin{pmatrix} \left( m+M \right) \frac{K}{r_{m}} & 0 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \dot \Delta \left( t \right) \\ \dot \varphi \left( t \right) \end{pmatrix}\\ &+\begin{pmatrix} \left( m+M \right) 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{pmatrix} \begin{pmatrix} \Delta \left( t \right) \\ \varphi \left( t \right) \end{pmatrix}\\ &= \begin{pmatrix} ( m+M )( P \varphi ( ( j-1 ) \tau) +D \dot \varphi (( j-1) \tau ))\\[-1pt] +( m+M+\frac{1}{2}m_m \frac{r_w^2}{R_w^2}) r_m s r_0 \mathop{\rm sgn} \Delta (t) \\[6pt] \frac{1}{2}\frac{m}{m+M}l\frac{r_w}{R_w}sr_0 \mathop{\rm sgn} \Delta \left( t \right) \end{pmatrix}, \end{aligned} \end{equation} for $t \in [ j \tau, ( j+1 ) \tau )$ and $j=1,2,\dots$, if \( \left| \Delta \right| \geq r_0 \); and \begin{equation}\label{eqblsdin} \begin{aligned} &\begin{pmatrix} \frac{1}{2}\left( m+M \right) m_{m} r_{m} & -\frac{1}{4}m m_{m} l r_m \frac{r_{w}}{R_{w}} \\ 0 & \frac{1}{3}ml^2-\frac{1}{4}\frac{m}{\left( m+M \right)}ml^2 \end{pmatrix} \begin{pmatrix} \ddot \Delta \left( t \right) \\ \ddot \varphi \left( t \right) \end{pmatrix}\\ &+ \begin{pmatrix} \left( m+M \right) \frac{K}{r_{m}} & 0 \\ 0 & 0 \end{pmatrix} \begin{pmatrix} \dot \Delta \left( t \right) \\ \dot \varphi \left( t \right) \end{pmatrix}+ \begin{pmatrix} 0 & 0 \\ 0 & -\frac{1}{2}mgl \end{pmatrix} \begin{pmatrix} \Delta \left( t \right) \\ \varphi \left( t \right) \end{pmatrix}\\ &= \begin{pmatrix} \left( m+M \right) \left( P \varphi \left( \left( j-1 \right) \tau \right)+D \dot \varphi \left( \left( j-1 \right) \tau \right) \right) \\ 0 \end{pmatrix}, \quad t \in \left[ j \tau, \left( j+1 \right) \tau \right) \end{aligned} \end{equation} for $j=1,2,\dots$, if \( \left| \Delta \right|