LIMIT CYCLES FOR SYSTEMS WITH NONLINEAR FRICTION
By
MICHAEL R. JAMES
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
1989
ACKNOWLEDGEMENTS
I am indebted to many people for their cooperation and
assistance over the long period of my doctoral studies.
First, my adviser, Dr. Thomas Bullock, and the other
members of my graduate committee, for their suggestions,
help and cooperation. Dr. V. M. Popov provided a very
helpful discussion on symmetry. Dean Eugene Chenette and
Art Zirger were very helpful during my coursework and
residence, as were many other UF personnel too numerous to
mention.
I received help from many fellow employees of Harris
Corporation: Pete Pitard, for his hard work in obtaining
financial assistance, Lee Almond, for his support, Kevin
Arter, for his review and suggestions, and the
Training/FEEDS and Human Resources personnel. Special
thanks go to Rich and Donna Phelan, whose enthusiastic
support solved many problems during the first years.
Most of all, I want to thank my wife, parents, and family
for their understanding and support during these six long
years. Without their sacrifices I could not have succeeded,
nor would it be meaningful without them.
All of these people share in this achievement, and I
deeply thank them for their support.
TABLE OF CONTENTS
ACKNOWLEDGMENTS............................ .. ....... ii
LIST OF FIGURES ................................. ... v
ABSTRACT............................................ vi
CHAPTERS
I INTRODUCTION............................... 1
Problem Statement .......................... 1
The Piecewise Linear Method................ 4
Illustrative Examples...................... 8
Organization of Dissertation............... 21
II REVIEW OF EXISTING RESULTS................. 23
An Unpopular Area of Study.................. 23
Phase Plane (2D) Analysis of Friction...... 24
Literature on Piecewise Linear Method...... 25
Extensions to MultiDimensional Case
(Approximate Analyses by Describing
Function) ................................ 26
General Results on Periodic Orbits.......... 27
Summary of Literature Review............... 28
III DEVELOPMENT OF STANDARD MODEL.............. 30
System Model with Physical State Variables. 31
Modelling Friction By a Nonlinear,
TwoInput Function .................... 32
The Physical State Model in Piecewise
Linear Regions..... .................... 34
Conversion to Control Canonical Form........ 38
Simplification of DE Solution Via
Coordinate Translation................ 41
Summary of Assumptions and Discussion of
Generality of Model ................... 43
iii
IV EXISTENCE CONDITIONS FOR LIMIT CYCLES...... 49
Necessary Conditions for a Simple Limit
Cycle ................................. 49
Illustrative Examples...................... 59
Exact (Necessary and Sufficient) Conditions 68
V LIMIT CYCLE SYMMETRY AND STABILITY......... 71
Stability of Predicted Limit Cycles......... 71
Stability Calculations for Example Problems 91
Symmetry of Limit Cycles for Odd Systems... 100
VI CONCLUSIONS AND RECOMMENDATIONS FOR FURTHER
RESEARCH................................. 116
Results and Conclusions of Current Research 116
Recommendations for Further Research........ 120
APPENDIX ANALYSIS COMPUTER PROGRAMS ............. 123
LIST OF REFERENCES ................................. 137
BIOGRAPHICAL SKETCH ................................ 139
LIST OF FIGURES
Figure pa
Ii BLOCK DIAGRAM OF SYSTEM UNDER STUDY... 3
I2 TWOINPUT FRICTION MODEL.............. 3
I3 PIECEWISELINEAR REGIONS FOR 2D CASE.. 6
I4 LIMIT CYCLE AMPLITUDES FOR 2D CASE,
(B=l, J=l, LC=1)................... 10
I5 2D SLIDING LIMIT CYCLE, PHASE PLANE
PLOT (B=l, K=l). .................... 12
I6 BLOCK DIAGRAM OF 3D SYSTEM (EXAMPLE 2) 13
I7 3D STICKING LIMIT CYCLE, VIEWED
ALONG VELOCITY AXIS (B=l).......... 15
I8 3D STICKING LIMIT CYCLE, VIEWED
ALONG AXIS OF INTEGRATOR STATE
(B=1)...... ......................... 16
I9 3D STICKING LIMIT CYCLE, VIEWED
ALONG POSITION AXIS (B=l).......... 17
I10 3D SLIDING LIMIT CYCLE, VIEWED
ALONG AXIS OF INTEGRATOR STATE
(B=.9)..... .......................... 19
I11 3D SLIDING LIMIT CYCLE, VIEWED
ALONG VELOCITY AXIS (B=.9)......... 20
IV1 3D LIMIT CYCLE EXAMPLE STICKY CASE,
DETERMINATION OF PERIOD Tl (ZEROS
OF F(T1))............................ 65
V1 3D LIMIT CYCLE EXAMPLE SLIDING CASE,
LOCAL STABILITY NEAR XO............ 94
V2 3D LIMIT CYCLE EXAMPLE STICKY CASE,
LOCAL STABILITY NEAR X3............ 96
V3 3D EXAMPLE SYSTEM PHASE PORTRAIT,
VIEWED ALONG AXIS OF INTEGRATOR
STATE ............................. 97
V4 3D EXAMPLE SYSTEM PHASE PORTRAIT,
VIEWED ALONG AXIS OF INTEGRATOR
STATE (EXPANDED) ....... ............ 98
V5 3D EXAMPLE SYSTEM PHASE PORTRAIT,
VIEWED ALONG AXIS OF INTEGRATOR
STATE (EXPANDED) ........ ........... 99
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
LIMIT CYCLES FOR SYSTEMS WITH NONLINEAR FRICTION
By
MICHAEL R. JAMES
December, 1989
Chairman: Dr. Thomas E. Bullock
Major Department: Electrical Engineering
Exact conditions for the existence of limit cycles due
to friction in a servomechanism have been determined for the
general (ndimensional) case. These conditions are in the
form of nonlinear algebraic equations, for which a solution
consistent with certain assumptions must be found.
A twoinput model for friction is required for accurate
results, including both dry friction and static friction
components. Therefore, standard nonlinear analysis methods
such as absolute stability and describing functions have
limited applicability. A piecewise linear method is used,
which constructs the limit cycle trajectory with pieces from
each region of the state space.
The system model is first converted to a control
canonical form, which simplifies the analysis. Assumptions
made in the formation of the model are specified exactly;
the model represents a rotational servomechanism, but is
shown to be applicable to other systems.
vi
Next, the piecewise linear method is applied to an
assumed limit cycle. The resulting algebraic equations can
be solved for the exact limit cycle period, amplitude, and
trajectory, if any exist. The limit cycles found are of
simple type; that is, complex orbits such as those which
traverse a region more than once would not be found.
Necessity and sufficiency of the existence conditions is
proved.
An analysis of the local orbital stability of limit
cycles found by this method is also presented. Eigenvalues
of a given matrix determine the asymptotic stability of the
periodic orbit.
The property of symmetry of limit cycles for odd
systems in general and friction in particular is examined.
Special cases are presented in which symmetry holds, and
counterexamples for the general case of odd systems are
given. The question of symmetry for friction limit cycles
remains open.
Examples are included to illustrate the method. One
example is a 3D system that exhibits two limit cycles in its
phase portrait, each a different type.
Listings of useful analysis computer programs are
given, in addition to suggestions for further research in
this area.
vii
CHAPTER I
INTRODUCTION
Problem Statement
The objective of this research is to define exact
conditions on the existence and behavior of limit cycles due
to friction. The analysis method to be used is a piecewise
linear approach which, unlike approximate methods such as
describing functions or absolute stability, has the
potential for providing exact (i.e., necessary and
sufficient) conditions for existence.
This problem is of interest in the design of
servomechanisms, which must deal with the nonlinearities in
typical electromechanical devices, such as friction.
Experience has shown that the effect of friction is
generally benign, since it tends to provide increased
damping to the system. However, this cannot be guaranteed
for all systems. Furthermore, it is a standard practice to
precompensate for this effect by deliberately tuning the
system to be underdamped without friction, so that it is not
excessively damped in actual operation with friction. It is
postulated that this technique, when pushed to extremes,
results in largesignal instability, since the damping
effect of friction varies inversely with amplitude. It is,
therefore, a significant research problem to determine the
point where limit cycles can exist.
1
An approach has been developed during this research
that uses the piecewise linear method, that provides an
exact condition for the existence of a limit cycle in terms
of a nonlinear algebraic system of equations. When a
solution to the system is found (and thus a limit cycle
exists), the method automatically gives exact frequency,
amplitude, and state trajectories for the limit cycle. On
the other hand, exact conditions under which the nonlinear
system of equations is solvable are more difficult because
of the nonlinear nature of the system.
The specific problem under study is shown in block
diagram form in Figure I1. The plant in the standard
system is represented by a Laplacetransform model of a
motor/load combination, which includes the effect of load
inertia and damping. The model does not include compliance
in the load, although, as described in Chapter III, this may
be included. The inertia (J) and damping (B) represent the
lumped quantities for the load and motor, with any gear
ratio taken into account. The servo controller is
represented as a linear system with motor drive torque as
the output quantity. Note that despite the specific
description of the problem model, it will be shown in
Chapter III (Development of Standard Model) that the model
is actually quite general.
The model used for nonlinear friction includes both
static and sliding coulomb friction components (see Figure
12). Therefore, it is a twoinput nonlinearity and is
FIGURE I1. BLOCK DIAGRAM OF SYSTEM UNDER STUDY
L F(TORQUE)
L
C
SLIDING MODE Y
(FOR NONZERO VELOCITY, Yn O) VELO
F_ n(VELOCITY)
L
LF= LCsgn[Yn] c
L (TORQUE)
STATIC OR STICTION MODE
(FOR ZERO VELOCITY, Y =0) L
n S
LF= LS sat[LIN/LS ] L
(INPUT
L TORQUE)
/ S
FIGURE 12. TWOINPUT FRICTION MODEL
thus difficult to handle by standard techniques (such as
describing functions or absolute stability criteria). The
extra complexity caused by the use of the twoinput model is
necessary; examples are known where modelling sliding
friction only indicates a limit cycle exists, when, in fact,
the "sticky" effect of the static friction prevents it. On
the other hand, the sliding friction can provide much of the
damping of many servomechanisms which have little or no
viscous friction damping in the load. These comments
indicate that a model including both effects is necessary.
Note that the friction model used is odd, and therefore the
system derivative vector will be an odd function of the
state.
The Piecewise Linear Method
The piecewise linear method is an approach where
systems with piecewise linear nonlinearities can be analyzed
using the tools of linear system theory. Since the
nonlinearity is linear over portions of its domain, the
state space can be divided into regions in which the system
behaves in a purely linear fashion. The linear system
theory then applies, although there will be a different
linear system in each region. Solution trajectories must
then be pasted together from different regions through the
use of boundary conditions (called commutation equations in
relay servo or onoff system analysis such as in Weischedel,
1973).
The nonlinear friction to be analyzed in this study is
piecewise linear, even though it requires two inputs. Five
regions of the state space are required (Figure I3 shows
the case for 2dimensional state space):
(I) velocity < 0,
(II) velocity > 0,
(III) velocity = 0 and magnitude of accelerating input
torque < LS,
(IV) velocity = 0 and torque > LS, and
(V) velocity = 0 and torque < Ls,
where LS is the maximum static friction torque (breakaway
torque). The first two regions are always open halfspaces,
while the last three divide the (nl)dimensional subspace
where velocity equals zero into three regions. For example,
for the 3D case (n=3), the last three regions are portions
of the plane y3 = velocity = 0. As can be seen by
examination of the friction model, the nonlinearity either
has a constant output throughout each region (I, II, IV, and
V), or its output is proportional to a linear combination of
states (region III). Even a simple limit cycle must
traverse at least three of these regions, so the linear
systems of each must be considered in this limit cycle
analysis.
There are several advantages in the use of the
piecewise linear technique. Once a limit cycle is found
(i.e., solutions from four regions that can be pasted
together to form a closed trajectory), the exact trajectory
VELOCITY
Y2 OR DY1/DT
LINEOF 0F
ZERO N
INPUTTORQUE \
Y2 = (K/B) Y1
(FOR 2D CASE) '
REGION IV:
VELOCITY = 0
AND TORQUE > LS
REGION II:
VELOCITY > 0
SAMPLE PERIODIC
ORBIT
REGION III:
UNESEGMENT
OF EQUIL. POINTS
LS / K
Y1
REGION V:
\ / VELOCITY= 0
AND TORQUE <LS
1 \ I
 *~~~~~~^
REGION I:
VELOCITY < 0
FIGURE 13. PIECEWISELINEAR REGIONS FOR 2D CASE
POSITION
~
__
is known. Therefore, the limit cycle amplitude, period, and
characteristics are known exactly, rather than approximately
as in describing functions. Furthermore, the method leads
to a more exact analysis, since the commutation conditions
require the solution trajectories to match exactly at the
region boundaries. Finally, the technique has a geometric
approach that is intuitively appealing and lends itself to
graphical solution.
The method is applicable to the majority of servo
nonlinearities. Backlash, friction, and most other
nonlinearities found in practice are piecewise linear.
The disadvantage of the method is in the difficulty of
solving the commutation conditions for a solutionfinding
the trajectories that will "paste" together to form a cycle.
One way to surmount this problem applies the method as a
secondary step in the analysis in order to refine an
approximate solution found by other methods. For example, a
limit cycle period and amplitude are determined
approximately by describing functions; then the piecewise
linear method is used (with the approximate solution as an
initial solution trajectory) in an iterative way to converge
to the true solution. A composite computer analysis program
could integrate these functions for the convenience of the
servo designer. Therefore, the difficulty of solution does
not really limit the applicability of the method.
Illustrative Examples
Before going into the detailed development of the
method and other results, it is useful to present simple
example systems exhibiting friction limit cycles. The limit
cycles were analyzed by the methods to be presented later,
so we can return to these examples at the appropriate time
to illustrate the method; here only the results are shown.
The second example system is particularly interesting
since it is only 3D (three state variables describe the
system), yet it exhibits two distinct limit cycles of
differing types (when the appropriate parameter values are
selected). The fact that such a simple system can have such
complex behavior indicates the richness of the study of
piecewise linear systems.
Example I1: TwoDimensional (2D) System with Sliding
Friction Limit Cycle
The question of friction limit cycles for the case n=2
can be completely solved; in fact, the piecewise linear
method is the same as standard phase plane methods for this
case. For n=2, the controller can be a constant gain only
(since there are two states in the plant); the state
variable model of the system is then
(1.1) dyl/dt = y2
dy2/dt = (K/J) yl (B/J) y2 LF/J
where LF is the friction torque from the model in Figure
12, and any feedback gain on y2 is lumped in with B. It
can be shown (see Thaler and Pastel, 1962, pp. 96104) that
for B>0 and for K>0, the system is globally asymptotically
stable (i.e., no limit cycles). Thus a stable system
remains stable (in fact it is more stable) with friction (2D
case only).
The case B=0 is also stable; the linear system has
imaginary eigenvalues, the friction gives enough damping to
give asymptotic stability (trajectories spiral to an
equilibrium point that is not necessarily the origin).
Interestingly, a limit cycle exists for every case with
negative damping (B<0). If the position feedback gain K is
sufficiently large to give the linear system a complex pair
of poles, an unstable, unique limit cycle exists, whose
amplitude depends on K, B and the inertia J and sliding
friction level LC. This "sliding" type of friction limit
cycle is one in which the motor or other object being moved
never "sticks"; the reversing torque is sufficient when it
comes to rest to immediately breakaway in the opposite
direction. Therefore, the orbit consists of two parts, each
representing motion, with two switching per cycle. The
initial conditions can be expressed in terms of the system
parameters as
(1.2) X10 = (Lc/K) [(L+1)/(Xl)], x20 = 0
where
(1.3) p = exp[Br/(2JP)], 3 = [4K/J B2/J2]k
Figure I4 shows some values of amplitude (which equals x10)
as a function of gain K for the case J = LC = 1, B = 1; it
0
0
0
0
co
0
S Q
N
o
oo II
II
u E II
HO
p4
U
H
__________________g 0
*
O a
appears that a limit cycle of any amplitude can be found by
adjusting gain properly.
The unstable limit cycle bounds the region of
asymptotic stability: trajectories inside spiral in to
equilibrium near the origin, while those outside go to
infinity (Figure I5 shows phase plane plot). This example
shows why it is dangerous to design underdamped systems, and
depend on the sliding friction to stabilize the system; the
friction damping is amplitudedependent, and at some
critical amplitude, is insufficient for stability.
Example 12: ThreeDimensional (3D) System with Both
Sliding and Sticking Limit Cycles
The second example system to be considered is defined
by the differential equations
(1.4) dy/dt = A + b LF
where
(1.5) A = 0 1 0
0 0 1
K1 K2 B
with b' = [ 0 0 1 ], LF = friction torque, and where the
state variables are y2 = position, y3 = velocity, and
yl = compensator integrator. Figure I6 shows a block
diagram of this system.
Case I: One Stable Pole and One Imaginary Pole Pair (Only
Sticking Limit Cycle Exists)
For specific results, let us first set B = K1 = K2 = 1,
and set the sliding friction torque to 1.0 and sticky
o o 0 0 0
S0 C
> O OU u
C4
0
3
>4l
U
II
""
II
H
c I
II
1
H
oi a
3 &<
c?;~
1
FIGURE 16. BLOCK DIAGRAM OF 3D SYSTEM (EXAMPLE 2)
friction (breakaway torque) to 1.2. Simulation of this case
demonstrates a sticking limit cycle with an initial
condition of y0' = [ 1.0 0.2 0.0 ]. This type of friction
limit cycle repeats a fourpart cycle of sticking (until
torque integrates up to breakaway), sliding to a new
position, sticking again, and sliding back to the start
position. A visual display of this behavior can be seen in
Figures I7 through 19, which show projected views of the
closely similar limit cycle found using the Case II
parameters.
The solution of the piecewiselinear equations results
in a valid limit cycle solution where y = [ 1 0.2 0 ],
T1 (the sliding period) = r, and the sticking period
T2 = 10, yielding an overall limit cycle period of approx.
26.3 seconds, matching the simulation results.
Case II: One Stable Pole and One Unstable Pole Pair (Both
Sticking and Sliding Limit Cycles Exist)
If the same example system is used with B = 0.9, and
K = K2 = 1 still, the closedloop eigenvalues (of the
linear portion of the system) move into the righthalf
plane, to .026 j 1.024, while the other pole is at .9524.
As might be expected, there is still a sticking limit cycle
solution close to that of Case I (Figures I7 through I9).
Numerical methods give a solution at approximately
O = [0.98 0.22 0 ], and approximately the same period
(confirmed by simulation).
1 + +
OM E4OZ
I  I
o
0
oo
o
U<
Ei
o
o
I
0D
U
U II
E4
H
lX]X
uu
HE
i< H
UU
O
aH
0>
,H
Ct<
i
ll
O
0
0
O
0
0
0
0
______ _________ __________ _______ 
0
0
0
0
0
0
I
01
> c OUOU E
II
U
C E
1 0
uO
U E
HZ
U)
EH
co
z
o0
Z
H
HO
3M
CC 3
___ __ __ _ _ _ _ _ _ _4_____ I__ _ _
t t + +
0
CD
o o
o
E4
CD
o
I
0
4
The interesting point about this case, however, is the
similarity to Example I1, the 2D case; that is, an unstable
complex pole pair exists. Since the nonlinear damping
effect of the sliding friction is largest at small
amplitudes, we might expect trajectories to converge for
small amplitudes and diverge at large amplitudes where the
unstable linear poles overcome the nonlinear damping (as in
the 2D case). Therefore, this case would have both a
sticking and a pure sliding friction limit cycle!
Simulation having shown a limit cycle with a
halfperiod T1 approximately equal to 3.1 seconds, numerical
solution of the equations for the symmetric, sliding case
was performed, yielding a solution at T1 = 3.1455 seconds,
and y0 = [0.283 12.4 0.]. Figures 110 and I11 show
some views of the two 3D limit cycles for this case.
It is expected from physical intuition and simulation
results that the sticking orbit would be stable, while the
sliding orbit would be unstable. This is in fact the case,
with the sliding orbit exhibiting saddle point behaviorone
stable and one unstable mode (real eigenvalue / eigenvector
pair), while the sticking orbit has stable node behavior
two real stable modes. The analysis of the orbital
stability will be examined in Chapter V, where Example 12,
case 2 will be used to demonstrate the method used for limit
cycle orbital stability calculations.
O
O
o
r(
S 1
 O
 o
 1 I    \ I   I  I I I    I \  
> OU)HEq >
(4
o
o
4
0
 4
0
o
..I
 0
0 0 0 0
4 4 ( )
U.
0
U*
II
H'
H
*Z
rt>
L M
.x
U
11
Q &>
"Z=EHuC < E,OW OE4
Organization of Dissertation
The main body of this dissertation is organized as
follows. After a review of existing results on this problem
and related areas, the standard model introduced above is
developed more fully. The objective in the model
development (in Chapter III) is to provide more convenient
forms for analysis (a control canonical form and a normal
form), while maintaining the generality and traceability to
actual servomechanisms of the original model.
Once the model is in the proper form, Chapter IV
develops the nonlinear system of equations that represent
the necessary and sufficient conditions for the existence of
a friction limit cycle. Several special cases are examined
that are somewhat simplified from the general problem and
easier to solve. For example, if it can be assumed that the
limit cycle is symmetric about the origin (that is,
halfwave symmetric), the equation system can be simplified
considerably.
The benefits of symmetry in the solution of the limit
cycle conditions led to an examination of the conditions
under which this symmetry existed. Chapter V examines this
problem and describes special cases for which symmetry
exists. Unfortunately, it is not known if friction limit
cycles are symmetric in general, although a search for
counterexamples was unsuccessful.
Using the results of the previous chapters, a complete
stability analysis of friction limit cycles is also included
in Chapter V. Using the eigenvalues of the local stability
matrix developed in this analysis, the stability of limit
cycles previously found is determined. This allows the
global phase portrait of the system to be pieced together
from a definition of all the equilibrium points, limit
cycles, and their stability characteristics.
The results presented here are considered only a start
in an essentially new area: use of piecewise linear
techniques for limit cycle analysis. Many possibilities for
further research exist, both in the completion of the
solution for friction, and in application of the method to
other piecewise linear nonlinearities. As pointed out
above, the majority of the plant nonlinearities encountered
in servo design are piecewise linear, including gear
backlash, saturation, quantization, and dead zone. The last
chapters provide a summary of results and suggestions for
further research in these areas.
CHAPTER II
REVIEW OF EXISTING RESULTS
An Unpopular Area of Study
Surprisingly little work has been done in this area,
for various reasons. First, the nonlinearity is not "nice,"
i.e., continuous or smooth, which prevents the application
of many mathematical tools. Second, when modelled with two
inputs it really becomes a system with multiple
nonlinearities, for which standard approaches cannot be
used, and extremely few results are available.
Finally, and, I believe, most importantly, it is
perceived as generally benign in effect. It can be shown
fairly easily (see Thaler and Pastel, 1962, for example) in
the 2D case (i.e., two state variables), that if the linear
portion of the system is already stable, the system with
pure sliding friction added is also stable. .It seems
probable that this result could be extended to arbitrary
dimensioned systems, based on the following energy argument:
the linear system is stable (hence contracting), and because
the (sliding) friction opposes motion, it damps the system
further, so the map is more contracting than without it.
Therefore, friction is generally held to be benign in effect
(it increases the stability of a system by increasing the
damping).
For these and possibly other reasons, the references
available in this area were few, and no analysis of the
general multidimensional friction limit cycle problem was
found.
In spite of these comments, the study of friction limit
cycles has practical value for control system design.
Although the effect of sliding friction is to increase
damping, it has become fairly standard to precompensate for
this effect in a design. The system is deliberately
designed to be underdamped, so that it is not too well
damped when friction is added. It is conjectured that this
technique, when pushed to extremes, results in largesignal
instability, since the damping effect varies inversely with
amplitude.
In addition, the stiction component can cause
difficulties also. An example was given in Chapter I of a
stable (although not asymptotically stable) thirdorder
linear system which has a limit cycle when nonlinear
friction (including stiction effects) is added. A slight
variation in parameters results in a stable linear system
which destabilizes with the addition of nonlinear friction!
Phase Plane (2D) Analysis of Friction
As mentioned above, Thaler and Pastel (1962), in their
classic text on nonlinear systems, completely solved the
secondorder case for friction, including both sliding and
static components. They also give an exact criterion for
the existence of limit cycles when an input ramp is present.
They showed that the system had no limit cycles (for a
system with positive damping) for the zeroinput case.
Earlier work by the same authors, (Pastel and Thaler, 1960)
actually shows a stability boundary for backlash (as a
function of system damping), then demonstrates the
stabilizing effect of coulomb friction.
Analyses of the 2D friction problem were presented in
several other references from the period when phase plane
analysis was an active area of research.
Literature on Piecewise Linear Method
The approach used in this dissertation to give exact
solutions for limit cycles (if they exist) has been known
and used for many years. In the relay servomechanism
literature, it is known as the piecewise linear method.
The method was used primarily for relay systems, though
it is generally applicable to any piecewise linear
nonlinearity. Note that the piecewise linear method is
essentially the same as phase plane analysis for 2D systems,
so there is overlap between the phase plane references, such
as those cited previously, and those for the piecewise
linear method. Unfortunately, the graphical methods of
phase plane analysis cannot be easily extended to higher
dimensions (see comment below on the work of Ku).
The piecewise linear method was used at least as early
as 1963 (by Kovatch); additional studies can be found in
references such as Kovatch (1964) for two nonlinearities,
O'Donnell (1964) for timeoptimal switching, Marstrander
(1969) for backlash, Negoescu and Sebastion (1971), Langill
(1965), and Urabe (1967). Weischedel (1973) applied the
method to onoff control systems. Ku (1958) made
significant and early contributions, including attempts to
extend the phase plane method to higher dimensions by
various graphical projections.
The Russian controls literature calls this the Andronov
point transformation method; a reference was found in a
translation of material developed in 1956 by E.P. Popov
(translated 1962). Although significant development of this
area has been performed, especially in the Russian journal
Automatika, an application to the friction problem has not
been found.
Extensions to MultiDimensional Case
(Approximate Analyses by Describing Function)
Analysis using approximate techniques such as the
describing function go back at least as far as the piecewise
linear method, or even farther in the Russian journals
(Popov, 1956). Four papers are listed by Gibson (1963),
that he refers to as independent developments of the
describing function method; all are from the late 1940s!
However, the two input difficulty that limits the
applicability of describing functions to friction is usually
dealt with by ignoring the static friction. No references
have been found that discuss the general case for an exact
friction model. No other references were found that
attempted the multidimensional friction problem with other
methods.
General Results on Periodic Orbits
In the mathematics literature on the theory of periodic
solutions of ordinary differential equations, the Poincare
map (Hirsch and Smale, 1974) is frequently used; this map
represents the change in system state after one cycle or
period. By searching for fixed points of the map (which
represent a return to the original state after one cycle), a
periodic solution can be found. This theory was used in
generating the stability results to be stated in Chapter V,
where the orbital stability of friction limit cycles is
examined. The Chapter V results are, therefore, original
only in their application to friction limit cycles, not in
the theory of orbital stability.
There are many results in this literature (although the
Poincare theory itself is restricted to planar systems).
Pliss (1966) stated some results on the existence and
uniqueness of these fixed points, although they could not
apply to friction as they were restricted to continuous
functions, or nonautonomous systems. Hale (1969) discussed
fixedpoint theorems for 2D systems and ndimensional
systems on a torus; in addition he stated the orbital
stability result (that all eigenvalues of the map except one
must be less than one in magnitude). Hayashi (1964) refers
to this approach as the topological approach to periodic
solutions, since the theorist is concerned with the
existence and topological properties of integral curves in
statespace. Finally, Hirsch and Smale (1974) provide an
excellent discussion on an introductory level of the theory
of flow maps; this reference was invaluable in the stability
analysis of Chapter V.
My survey of this literature was thorough, but not
exhaustive; no applications to this nonlinearity were
encountered.
Summary of Literature Review
In the 20 years since the heyday of phase plane and
piecewise linear methods, it appears that almost no new
results have been discovered in this area. However, the
study of friction limit cycles has significant practical
value for designers of servomechanisms, since a servo
designed using standard techniques may exhibit instabilities
in the presence of friction.
The literature survey was extensive in the area of
controls, and moderate in the mathematics field. In the
former, IEEE indices were searched, including 25 years of
the Transactions on Automatic Control. Twentyfive years of
the Russian journal Automation and Remote Control were
examined. The journal Automatica was also reviewed.
Although no journals in the field of pure mathematics were
checked, many texts in the areas of stability, limit cycles,
and dynamical systems were examined.
The main contribution of this research is the
application of known (although relatively obscure in
controls literature) techniques such as the piecewise linear
method and Poincare flow maps to the problem of nonlinear
friction. No reference was found that considers the multi
dimensional friction limit cycle problem using an exact
method such as the piecewise linear technique. All the
references found either analyze the 2D case for friction, or
derive approximate results by using a simpler friction
model. Practical experience with servos indicates that it
is unacceptable to ignore either friction component.
The development presented here, therefore, appears to
be the first to attack the multidimensional friction limit
cycle problem with sufficient accuracy for practical
applications. In addition, the results on limit cycle
symmetry for piecewise linear systems, and the derivation of
limit cycle orbital stability by the same method, are
original, to the author's knowledge.
CHAPTER III
DEVELOPMENT OF STANDARD MODEL
This chapter is concerned with setting up the model of
the system to be studied, and converting it to a convenient
form for analysis. The original system model is in the form
of physical state variables, where the states represent
actual physical quantities in the servomechanism. The model
form to be used for most results is a control canonical
representation, which is developed from the physical states
by a linear transformation. This form has advantages in
clarifying subsequent developments of the theory, including
the role of system zeroes and characterizing limit cycle
symmetry. This chapter also demonstrates how the piecewise
linear friction torque input can be removed by a
translation, which leads to a simplification of the
existence formulas of Chapter IV.
In addition to the derivation of the state forms, this
chapter provides a description of the various assumptions
made in the analysis. In other words, the generality of the
analysis and the solutions are defined by examining the
assumptions made. For example, although the original
motor/load model included no compliance or resonance
effects, this restriction can be removed, allowing the
results to be more generally applicable.
System Model with Physical State Variables
The two primary physical variables are angular velocity
(required as input to the sliding mode friction) and its
integral, angular position. The source of this model is a
rotating servomechanism. The theory is not restricted to
rotational friction problems, however, since angular
variables (torque, angular velocity, etc.) could be replaced
by linear ones (force, velocity, etc.) without affecting the
nature of the problem (this point is discussed further at
the end of this chapter). The model was illustrated in
block diagram form in Figure I1; the controller transfer
function G(s) is replaced by state equations in the model
below.
The remaining n2 state variables, where n is the order
of the system, are controller states. That is, they form
the model of the servo loop compensation, amplifiers, motor
armature effects, and other dynamic portions of the system.
Thus we may set up a model of the following form:
(3.1) = Ac y + a y + a c2yn
c c cl nl c2 n
Yn1 = Yn
J Yn = Kyn1 Byn + b'yc + f(y)
where the yi variables are the physical states and make up
the vector y, the dot (') indicates the time derivative
d/dt, prime (') indicates transpose, f(.) is the friction
torque, J is the lumped system inertia, and B is the
damping. The n2 dimensional linear subsystem in the n2
states in vector yc is formed by the system matrix Ac
c c
(n2 x n2), the input vectors a cand a (n2 x 1), and
the output vector b (n2 x 1) that feeds the output of the
controller into the torque equation.
The only assumptions made in the construction of this
model are that yn1 and yn are angular position and
velocity, respectively. These assumptions completely define
the first equation, and require that the nonlinear friction
torque term f/J appears in the y differential equation
(representing angular acceleration due to the sum of torques
applied). The rest of the system is linear; the friction
term does not appear in any other derivative.
Note that there is an assumption made here that
restricts the generality of the model: that the friction
appears in only one equation. Although examples can be
constructed that do not obey this restriction (such as those
involving differentiators, full state feedback, torque
measurement, or a foundation model), this model is felt to
be quite general, and covers essentially all servo problems
where friction limit cycle information is of.interest. This
point will be discussed in more detail later in this
chapter.
Modelling Friction By a Nonlinear, TwoInput Function
The nonlinearity investigated in this analysis is a
model for mechanical friction that depends on both the
relative velocity of the surfaces and the driving torque of
the moving element. The level of friction opposing the
motion depends on velocity (or more precisely, the direction
of the velocity) when the element is moving. In the case
where the element is stationary, on the other hand, the
static friction torque holding the system motionless must be
set equal to the input driving torque that is attempting to
move the element. Each of these cases is based on the
models of Coulomb friction described in any basic statics
textbook. Figure I2 in Chapter I illustrates the friction
model.
Although an analysis can be performed with only one
element modelled, thereby using a singleinput nonlinearity
and simplifying the analysis, this would lead to results
that would not apply in any physical system, and could be
misleading or erroneous. Cases can be constructed where the
modelling of sliding friction alone predicts a limit cycle,
due to its highly nonlinear, discontinuous characteristic.
In a practical system, however, no limit cycle would exist.
For example, the twodimensional friction limit cycle
is completely solved in Thaler and Pastel (1962), as
described in Chapter I. Neglecting the static friction
component in this analysis leads to the conclusion that a
limit cycle of very small amplitude can be obtained.
However, in reality, such a small initial condition will
stick at the initial point due to the static friction
torque.
On the other hand, the sliding component of friction
cannot be ignored either. As a qualitative example of a
case demonstrating this, consider a system with no damping,
which would limit cycle if static friction alone were
present. The system starts at an initial position with
nonzero error, integrates up in the controller until it
breaks free from the static friction, and overshoots the
desired position to start the cycle all over again. The
inclusion of sliding friction in the model, however, might
provide enough damping to the system to cause it to spiral
in to the equilibrium point at the origin. Another
possibility indicated by the behavior of typical friction
limit cycles found is that a system that is actually stable
for small signals (bounded region of stability) might appear
unstable if analyzed without sliding friction.
These cases indicate that accurately predicting limit
cycles for physical servomechanisms requires a complete
model of friction, as is used in the analysis contained in
this paper.
The Physical State Model in Piecewise Linear Regions
As described in the introduction, application of the
piecewise linear method to the system model results in five
regions of state space, with a linear system for each
region. The form of the (now linear) model in each region
is described, so that the trajectories in each region can be
described.
Region I: Velocity < 0 (v <0). This ndimensional region is
described by the same set of differential equations as in
(3.1), except that the nonlinear friction term f(y) is
replaced by a constant input Lc. This follows from the fact
that the friction in sliding mode is constant (independent
of velocity) as long as the system is in motion in one
direction. From basic mechanics, this constant friction is
a function of the normal force, with the proportionality
constant being the coefficient of sliding friction. The
torque Le is positive, since physical friction always acts
to oppose the direction of motion; velocity is negative, so
the acceleration torque due to friction is positive.
Therefore, the model in region I is a linear system driven
by a constant input.
Region II: Velocity > 0 (y >0). This region model is
exactly the same as in region I, except that f(y) is
replaced by Lc
Region III: Velocity = 0 (y =0) and Acceleration Torque 5
Breakaway Torque. In this case, the system is in a static
condition, and the model for the static component of
friction applies (also known as stiction). The friction
opposes the applied acceleration torque with an equal torque
so that the net acceleration torque is zero,up to a maximum
amount of static friction. The maximum amount of friction
(the breakaway torque) is determined by the coefficient of
static friction and the normal force. Once the acceleration
torque becomes greater than the maximum stiction, the net
torque becomes nonzero, and the system begins to move
(enters region I or II).
Note that the acceleration torque used to define this
region, and which is matched by an opposing friction force,
consists of every term in the torque equation (the
derivative of y ) except the friction f. Therefore, the
linear terms in the torque equation define the boundaries of
this region. The region is contained in the
(nl)dimensional hyperplane y n=0.
Since the net torque is zero and the velocity is zero,
the first two differential equations in (3.1) drop out (the
state yni remains constant and y remains zero over
trajectories in this region), leaving the n2 dimensional
controller system as the only equations required to define
the trajectory in region III. Of course, input vector ac2
drops out (since yn=O), while acl and b are used to model
the effect of the constant angular position on the
controller and the controller output torque (to determine
the breakaway condition).
Region IV: Velocity = 0 (y =0) and Acceleration Torque >
(+)Breakaway Torque. This region is also contained in the
(nl)dimensional hyperplane Yn=0, but the driving
acceleration is sufficient at any point to overcome the
maximum static friction opposing the impending motion. Thus
a system that comes to rest with sufficient torque
(trajectory enters region IV rather than III) immediately
breaks into motion again, entering region II (since
acceleration torque is positive, velocity increases from
zero and trajectory must enter region II).
For this reason, a trajectory can be said to cross this
region (in the sense of crossing a plane, for example), but
is only in the region for an instant. Since the integral of
a finite quantity over a set of measure zero is zero, the
applied torque and the differential equations are irrelevant
to the system behavior. The trajectory leaves the region at
the exact point it entered. It is therefore unnecessary to
examine the form of the system equations in this region.
Region V: Velocity = 0 (y =0) and Acceleration Torque <
() Breakaway Torque. The region V model is identical to
the region IV model, except for the breakaway conditions.
The five regions cover the state space, including the
equilibrium point (actually, line segment) at the origin.
The trajectories in each of the individual regions can be
described by the solution to the linear system of
differential equations:
(3.2) y(t) = F(t,t0) y(t0) + G(t,t0) L
F(t,t0) = exp[A(tt0)]
where L is the (constant) input in that region, A is the
system matrix for that region, and the F and G matrices vary
from region to region. The initial conditions y(t0 ) are
determined by the commutation conditions at the boundary of
the region at which the trajectory enters.
An algebraic trick, to be presented below, can be used
with this form to remove the input term (GL) in equation
(3.2). The system then behaves as an autonomous system,
with considerable simplification of the defining equations.
Conversion to Control Canonical Form
The most convenient form for analysis is the control
canonical form, because the role of system poles and zeroes
in the existence conditions in Chapter IV will be clarified
through the use of this representation. In addition, a
result in Chapter V on symmetry is proved using this form.
Although this section constitutes a digression from the main
line of development, it is justified by the usefulness of
the control form.
The conversion to control form is accomplished by
representing the controller (transfer function G(s) in
Figure I1) in control form and combining this
(n2)dimensional system with the equations for position and
velocity. When completed, the system model will have the
following form:
(3.3) x1 = x2
X2 = x3
Xn1 = xn
xn = a0X1 alx2 .. axn + f()
The coefficients in the differential equation for x are the
n
coefficients of the closed loop polynomial of the linear
system.
The procedure used to convert the system to control
canonical form is as follows. The physical state model is
set up so that the controller (the (n2)dimensional
subsystem represented by state variables yc) is in control
canonical form already. This can be done for any controller
representation that is controllable (if the original
representation for the controller was a transfer function, a
controllable representation can always be found).
Next, the model is transformed so the velocity feedback
into the controller subsystem is zero (i.e., vector a in
c2
(3.1) is zero). This can be done without loss of generality
in the model, as shown by the following argument: Suppose
ac2 is not zero in the original system. We can apply a
similarity transformation
(3.4) y = Ta (1)
where y(1) is the new state variable vector, and
T = 1 0 0 . 0 a 0
0 1 0 Oa 0
c2
n2
O0 0 O l a 2 0
c2
0 0 0 1 0
0 0 1
Note this is equivalent to replacing all (n2) controller
states yc by (Yc(1) + ac2 Ynl) A straight forward
calculation (plug (3.4) into (3.1)) shows that the resulting
system has no feedback of velocity into the controller.
Since a similarity transformation (or the equivalent
change of variables) does not affect the basic system
behavior, but merely its representational model, this
alteration of the model is valid. Physically, this
operation has removed a redundant parameter in the original
model, since the velocity feedback in the original model can
be achieved by adjusting the damping parameter, B, or
adjusting the position feedback gain, K1. Since the models
are equivalent, the superscript on the y state variables can
then be dropped in subsequent equations.
With the controller in canonical form, and the velocity
feedback to the controller eliminated, the following
intermediate form is obtained:
(3.5) yc = Ac + a yn
Yn1 = Yn
J n = K1n1 By + b'c + f()
In order to complete the conversion of the model to
control form, two similarity transformations are then
applied to this intermediate model in succession. As
previously explained, similarity tranformations do not
change the basic behavior, hence are valid alterations of
the model equation. The transformations are:
1
(3.6a) A, = T1 A T1
1
(3.6b) A2 = T21 A T2
where the nxn matrix A is obtained from the model in (3.5),
and
(3.7a) T1 = 1 0 0 0
0 1 0 0
0 0 1 0 0
0 1 n3 1 0
0 .0 0 1
(3.7b) T2 = 1 0 0 0
0 1 0 0
0 0 1 0
0 60 1 n3 1
and where the 's are the coefficients of the
characteristic polynomial of the controller alone (i.e., the
characteristic polynomial of matrix Ac in equation 3.5).
The transformed system is then in control canonical form.
Note that the friction input is unchanged by these
transformations, and is still applied in the yn differential
equation. The composite transformation is required for
later derivations, as it defines the relationship between
the physical states and the control form states:
(3.8) y = T1T2x
where
T = TT = 1 0. 0 0
0 1 0 0
0 1 0 0
P0 16 n3 1 0
0 60 91 n3 1
Simplification of DE Solution Via Coordinate Translation
Once the system is in control form, the following
translation of the state coordinates can be used to convert
it to an autonomous system (within each region). Since the
friction is a constant in region I or II (L or L ), and
C C
since none of the differential equations except that for x
n
involves xl, the xl coordinate can be translated by
(3.9) xla = x1 f/a0
where the subscript "a" indicates autonomous coordinates,
and f is the (constant) friction level in that region (I or
II). Substitution of this translation shows that the x
n
equation still gives the same result, the system matrix is
unchanged, and the input has disappeared. The autonomous
system model is then
(3.10) xla = x2
k2 = 3
n1 = n
x = ax x nx
n 0 la 1 2 n n
Therefore, by performing a translation on the initial
condition and on the final state of the trajectory in the
region, the trajectory can now be described by the simpler
form
(3.11) x (t) = F(t,t0) a (t0)
F(t,t0) = exp[A(tt0)]
instead of equation (3.2), where the matrix A is in control
form (i.e., the matrix A2 produced by transformation 3.6).
Of course, the control form representation is not
required to perform this translation. The equivalent
translation in terms of the physical state variables can be
defined using the transformation matrix as
(3.12) ly = T2 1
Ya = Y ly
where 1' = [ f/a0 0 . 0 ] is the translation of the
control form states (equation (3.9)) and Ya is the
autonomous physical state variables. The resulting
differential equation for y would be of the form
(3.13) ya = A Ya
(with A unchanged). A slightly different form of the
existence conditions of Chapter IV would then result. The
control form has certain advantages, however, as will be
seen.
Summary of Assumptions and Discussion of Generality of Model
The following assumptions and conditions on the system
are made during the analysis of friction limit cycles in
this dissertation:
(Al) The system can be modelled by a model of the form
in Figure I1, and equations (3.1).
(A2) The velocity feedback into the controller is zero
(without loss of generality, as shown by argument given
previously in this chapter).
(A3) The controller portion has a controllable
representation (this is assured if the controller is
representable by a strictly proper transfer function; if
transfer function is proper the direct feedthrough term can
be divided out and lumped with the direct position feedback,
so this case is also allowed).
(A4) Friction appears only in the differential
equation for velocity (the acceleration torque equation).
(A5) The closedloop, linear portion of the system
should not have a pole at the origin.
(A6) The limit cycles to be analyzed are all simple,
that is, they traverse regions I and II once before
returning to the initial point, and therefore have only four
switching per cycle at most.
Note that many of the assumptions were made for
convenience in the subsequent development, so that a
specific model could be used for concreteness. Most can be
removed by a transformation of the original problem as
discussed below, or by minor modifications of the theory,
and so involve no loss of generality.
In fact, the only restriction to the system model that
cannot be removed is that the limit cycle is simple
(assumption A6). The analysis could be extended, as
discussed below, to remove even this restriction. Thus the
theory is quite general.
Several of these assumptions require further
discussion. Assumption (A4) is required to perform the
analysis using the piecewise linear method in Chapter IV,
since the analysis uses the control canonical form heavily.
If the friction input appears in more than one differential
equation, the trick used there to eliminate the friction by
a coordinate translation is invalid.
Examples could be generated that violate this
assumption. If acceleration torque was directly measured,
the friction torque could be inferred and used as an input
to the controller. The same comment holds for systems with
differentiators, since the derivative of the velocity
involves the friction. Friction also causes a reaction
torque on the mounting of the motor or other physical
device, so if a foundation or mount model is included, the
assumption would be violated.
This condition does not restrict the generality of the
method, however. Systems which violate this assumption can
still be transformed to control form, the method applied,
and existence conditions checked. A minor modification of
the Chapter IV conditions (involving the P vector only)
would be required.
Assumption (A5) is not necessary, except for steps
involving taking the inverse of certain matrix quantities in
Chapter IV. However, the assumption results in no loss of
generality, for the following reason: a system with a zero
eigenvalue has a free integrator so that, with the proper
selection of state variables, one state does not feed back
into the system. This can be seen by examining the control
form of the system model; when an eigenvalue is zero, the
coefficient a0 of the characteristic polynomial is zero, and
the first column of the system matrix is zero (see equation
3.3). Therefore the magnitude of the state xl has no effect
on the behavior of the system. Any limit cycle existing in
the system will thus also exist in the (nl)dimensional
system formed by deleting the state xl (although not
necessarily vice versa). Therefore, any system that fails
this assumption can be analyzed by deleting free integrator
states, then looking for limit cycles; any limit cycles
existing in the original system would be found by this
method. Note, however, that solutions may be introduced
that do not exist in the original. For example, an
asymmetric limit cycle in the reduced system might not exist
in the original since the free integrator output might not
be periodic (see Chapter V).
Assumption (A6) is also necessary to allow the
derivation of Chapter IV to proceed. The piecewise linear
method requires the overall shape of the limit cycle to be
known, including which regions are traversed and the order
of traversal. This could be considered a limitation of the
method, since the analysis cannot be performed without
knowing the shape of limit cycles that are to be found by
the analysis!
Note, however, that an even more restrictive assumption
is made in the classical describing function analysis, a
tool that has proven quite useful for decades. That
analysis assumes a sinusoidal input to the nonlinearity,
therefore implicitly assuming a halfwave symmetric type of
solution with two switching per cycle. It is therefore
felt that the assumption is reasonable. Future research
could remove this assumption by listing all possible
solution shapes (and order and number of traversals of
regions during each cycle), and applying the method to each
possibility. The difficulty lies in the large number of
potential solutions to be checked.
If it could be proved that all limit cycles are simple,
the difficulty would disappear. Without such a proof, the
results of this thesis apply only to simple limit cycles.
For example, the necessary conditions for the existence of a
limit cycle found in the next chapter are necessary only for
simple limit cycles. Therefore, the lack of a solution for
these conditions for a particular system only guarantees
that a simple limit cycle cannot exist for that system.
More complicated behaviors are not ruled out.
Finally, some comments on the generality of the system
model assumed are in order (assumption (Al)). The model was
inspired by the case of a rotational servomechanism, such as
a DC motor. Although this is an important application, the
theory developed here can be applied to many other systems,
as long as they can be put in the form specified. For
example, a linear drive system has the same characteristics
of inertia, damping, friction, and acceleration force
(instead of torque), so the theory can be directly applied.
Systems with geartrains and loads can be modelled by
lumping the load and gear inertia into the motor, as long as
the system is rigid and closely coupled (gear backlash can
be ignored). Although the model is for a continuoustime
system, the use of discretetime (i.e., digital) controllers
can be accomodated by using their continoustime
equivalents.
Many systems of practical interest have mechanical
flexure, modelled by resonances. If they can be lumped into
the model of the controller, the theory can still be applied
as is. However, the case of two masses, separated by a
spring, where both masses experience nonlinear friction,
does not fit into the present model. Additional piecewise
linear regions must be specified, with a great increase in
the complexity of the analysis, unless the limit cycle is
known to be low in frequency so the system could be
considered to behave as a rigid body.
In summary, certain assumptions are made to allow a
specific, concrete model to be used in the development.
Examination of the assumptions shows that most are for
convenience only, and the theory applies to a quite general
class of problems.
Given these restrictions and assumptions on the problem
to be studied, and the development of a convenient
representation, the piecewise linear method can be applied
to the problem of friction limit cycles. The next chapter
derives exact conditions for the existence of these limit
cycles, and, once a solution is found, provides exact
information on the limit cycle trajectory.
CHAPTER IV
EXISTENCE CONDITIONS FOR LIMIT CYCLES
This chapter applies the piecewise linear method to the
problem of friction limit cycles, using the model developed
in the previous chapter. The method provides a set of
nonlinear algebraic equations for which a solution must
exist in order to have a limit cycle. These represent
necessary conditions for limit cycle existence, since every
(simple) friction limit cycle for this system model meets
these conditions. The issues of sufficient conditions and
minimal sets of necessary and sufficient conditions are
discussed. Finally, some examples applying these conditions
are presented.
Necessary Conditions for a Simple Limit Cycle
In order to simplify the analysis, the assumption is
made that the limit cycle which the system exhibits is
simple, in the sense that it traverses each of the four
statespace regions exactly once during each limit cycle
period. This is assumption (A6) discussed in Chapter III.
This assumption makes sense from a physical viewpoint, since
it is hard to imagine a system where the velocity changes
sign more than twice before returning to the initial state.
It is clear that this assumption holds in the 2D case,
since the trajectory must go to the left in region I and to
the right in region II (see Figure I3 and I5) and the
49
trajectory cannot cross over itself. However, no claim can
be made that the following analysis is completely general,
due to this assumption. Therefore, although this is a
restriction of the generality of the results, the solution
of the equations requires the limit cycle meet the following
definition:
Simple Limit Cycle: A periodic orbit which traverses each
piecewise linear region once at most (at most twice for
region III).
As a consequence of this assumption, there will be
exactly four distinct switching instants during the period.
The system will be in a particular region, and, therefore,
its behavior will be governed by that region's system model
equations, for T. seconds, for i = 1, 2, 3, 4. As will be
seen below, some of these switching periods may be zero.
Thus the total limit cycle period T is the sum of the four
T.'s:
1
(4.1) T = T1 + T2 + T3 + T4
Note that the subscript on the periods indicate order, and
not the region being traversed during period i.
Since the limit cycle is periodic, any initial point
may be chosen to start the analysis. Define x to be the
system state at t=0. Assuming a limit cycle exists, this
initial state can be chosen to be the point in region III
(or V, see Figure 13, p. 6) where motion is impending. In
other words, initial velocity is zero, and the initial
driving torque is negative, and sufficient to overcome
static friction and begin motion. These requirements on the
initial point can be expressed as (recalling that y
represents the physical state variables, and x the control
canonical states):
(4.2a) Yn(0) = [ 0 60 n3 1 ] x
= 1' 0X = 0
(4.2b) yn(0) = H' A x" < LS (or = LS )
where the notation A' is used to represent the nth row of
the x to y transformation in (3.8) (used to calculate the
velocity yn from the state x), the parameter LS is the
maximum static friction before breakaway, and the matrix A
is the control canonical form of the state matrix, in region
I. Note that the A matrix from region I (velocity < 0) is
used since the rate of change of yn is identically zero in
region III, by design.
Having characterized the requirements on the initial
point, the limit cycle trajectory can be followed in region
I (which the system transits after breaking free) using the
region I model. The autonomous state variables can be
formed as discussed in Chapter III, by making the
translation:
(4.3) xla = x (L/an)
or
x = x 
a
with
1' = [ (Lc/an) 0 0 .
0 ]
where a is the value in the nth row and first column of A,
and LC > 0 is the sliding coulomb friction torque. This
translation is selected so that the friction torque is
positive, and opposes the motion in region I in the negative
direction. Note that the parameter anl cannot be zero,
since this would give a zero column in A, which signals the
presence of a free integrator.
The trajectory in region I is then
At At
(4.4) x (t) = e x = e ( 1)
a aO (
so that the trajectory exits region I at the point where y
(velocity) again goes to zero
(4.5) x(T1) = xa(TI) + 1
AT
= e (x 1) +1
with the requirement that n (T ) = 0:
(4.6) n(T1) = l' X(T1) = 0
At this point, the limit cycle will display two types
of behavior, depending on the driving torque at t=T If
the trajectory transits from region I into region V, this
torque is already sufficient to break loose the system and
drive back in the other direction, then T2 = 0 and the
system spends no time in region V (in other words, system
does not stick, but is only momentarily motionless). As
discussed in Chapter III, the region V equations do not have
to be applied, since the entry and exit points for the
trajectory in this region are the same. If this also occurs
at the other switch point (i.e., T4 = 0), the limit cycle is
driven purely by sliding friction.
The other case is more complex to analyze, since the
trajectory must be calculated through all four switching
periods. If the driving torque is insufficient to overcome
the static friction at t = T1, the mechanical load stops for
a period T2 until the controller can integrate up and break
it loose. The system has come to rest in region III after
leaving region I, so system behavior is governed by the
model for region III during this period. Both static and
sliding friction effects will affect this limit cycle. For
this case, the system state evolves as
(4.7) x(t) = es (tT1 (T) T1 < t < (T1 + T2)
from the initial point when the velocity went to zero at
t = T Note that no translation is necessary, since the
friction term is absorbed into the system model in the
static condition, represented by state matrix A
In either the pure sliding friction case (T2 = 0) or
the sticky case, the requirements on the system state at
t = T1 + T2 is that the driving torque is sufficient to
break loose:
(4.8) Yn(T1 + T2) = V'A x(T1 + T2) > LS (or = LS)
Since the sticky case involved the system evolving until
sufficient torque was developed to break away, the system
will move as soon as torque equals the static friction, so
(4.8) becomes an equality; note that the same comment
applies to condition (4.2b).
There is also the implicit requirement that the
velocity is zero when motion is impending (as in condition
(4.2)), but this requirement is not placed explicitly, since
it is taken care of by the model equations in region III
(which force acceleration to be zero) plus the zero velocity
at t = T1.
The system goes through the same behavior as above, in
the opposite sense, while traversing region II and returning
to III (or IV). The equations for the trajectory in these
regions, along with the conditions on the state at the
switching instants, are found by first performing the
translation to autonomous coordinates for region II, as in
(4.3):
(4.9) xla = x + (Lc/an)
or
x = x + 1
a 
with
i' = [ (Lc/a ) 0 0 . 0 ]
as before. Note the friction torque is negative in this
region for LC > 0, since the velocity is positive. The
trajectory in region II is then
(4.10) x (t) = eA(tT1T2) x (T1 + T2)
= eA(tT1T2 (x(T1 + T2) + i),
for T1 + T2 < t < T1 + T2 + T3
and the trajectory exits region II at the point where y
(velocity) again goes to zero
(4.11) x(T1 + T2 + T3) = x (T1 + T2 + T3) 1
AT
= e 3 ( x(T1 + T2) + 1) 1
with the requirement that y (T + T2 + T) = 0:
(4.12) Yn(T1 + T2 + T3) = x(T + T2 + T3) = 0
At this point, the system trajectory either immediately
breaks free (trajectory in region IV) and leaves the region,
or it sticks (trajectory in region III) and follows the
trajectory defined by
(4.13) x(t) = e (tT1T2T3) x(T1 + T2 + T3),
for T1 + T2 + T < t < T + T2 + T3 + T
from the point when the velocity went to zero. Note that no
translation is necessary, since the friction term is
absorbed into the system model in the static condition as in
the case of region III.
Collecting the conditions (4.5, 4.7, 4.11, 4.13) for
the trajectory through all four regions, and placing the
requirement that the endpoint match the initial point, we
have the following nonlinear algebraic matrix equation for
the initial point:
(4.14) AT AT AT AT
(4.14) I es 4 e 3 eAs 2 e 1) x = 1
with
AT AT AT AT
1 = es 4 (e 3 [e s 2 (I e A} 1 + 1] 1)
AT AT AT AT AT AT AT
= es4 e 3 eAsT2 e 1 1 + e s 4 e 3 e s 2 1
AT AT AT
+ e s 4 e 3 1 es 4 1
In addition, there are four scalar nonlinear equations
from the four conditions on switching points (4.2a/b, 4.6,
4.8, and 4.12, but the last one is redundant). Altogether,
there are n+4 (nonlinear) equations in the n+4 unknowns x0'
T1, T2, T3, and T4. In the pure sliding friction case, of
course, T2 and T4 are both zero, but conditions (4.2b) and
(4.8) become inequalities and thus drop out, and (4.14)
simplifies to
AT AT
(4.15) {I e 3 e 1) 0 = 11
with
AT AT
1 = e 3 [ (I e 1 1 + 1 ] 1
1
= eA3 eAT 1 1 + 2 eAT3 1 1
so there are n+2 equations in the n+2 unknowns xg, T1, and
T3 plus the two inequality conditions on the switch
points.
The results presented so far in this chapter can be
summarized as
Theorem 41: A simple limit cycle exists for the system
model under study only if the following algebraic conditions
hold on the variables xO, T1, T2, T3, and T4:
Sticking Limit Cycle Case
AT AT AT AT
(4.14) (I e s 4 e 3 e s 2 e 1) 0 = 1
where 1 = es 4 (e AT3 [es 2 (I e AT1) 1+ 1] 1)
AT AT AT AT A.T AT AT
= e s 4 e 3 e s 2 e 1 1 + e s 4 e 3 e s 2 1
AT AT AT
+ esT4 eA3 1 es4 1
(4.12) Yn(T1 + T2 + T3)
=[ 0 0 Pi n3 1 ] x(T1 + T2 + T3)
= (T1 + T + T) = 0
(4.2b) Yn(0) = x' A x0 = LS
(4.6) Yn(T1) = X(T1) = 0
(4.8) Yn(T1 + T2) = 2 'A x(T1 + T2) = LS
Pure Sliding Limit Cycle Case (T2 = T3 = 0)
(4.15) {I e 3 e 1) x0 = 1
0 1
AT AT
where 1 = e 3 [ {I e 1) 1 + 1 ] 1
1l
= eA3 eAT 1 1 + 2 e AT3 1 1
(4.2a) y (0) = [ 0 P . n3 1 ] X
= X0 = 0
(4.6) Yn(T1) = x(T) = 0
In other words, the conditions above are necessary for the
existence of a simple limit cycle. Condition (4.12) was
included in the sticking case instead of the equivalent
condition (4.2a) to obtain a more symmetric form of the
conditions.
Note that the conditions requiring velocity = 0 could
be applied at either t = T1 or t = T1 + T2 (equation (4.6)),
and at either t = 0 or t = T1 + T2 + T3 (equation (4.2a)).
This is a consequence of the fact that velocity remains zero
between these times.
Proof: The derivation presented above is a stepbystep
solution, through each region, for the limit'cycle
trajectory that was assumed to exist. Since the conditions
stated in the theorem were derived directly from the
trajectory solutions as seen above, and since a simple limit
cycle by definition must traverse those regions as assumed,
the conditions are necessary, and the proof is complete.
Note that a nonsimple limit cycle, should any exist,
would not have to meet the necessary conditions. Thus, if a
solution to the necessary conditions exists for a system, a
limit cycle may exist; if no solution exists, no simple
limit cycle exists, but a nonsimple one is not ruled out!
A complete theory requires an examination of the possibility
of more complex limit cycles.
It is interesting to note that the same development
could be performed using physical state variables. This
would result in replacing the control canonical A matrix in
the equations by the original system matrix, eliminating the
need for the vector P, and replacing 1 by the translation in
equation (3.12). The form used here clarifies the role of
system poles and zeroes, however, as seen by the following
discussion.
The definition of the initial state X0, equation
(4.14), is completely determined by the poles (for fixed
switching periods), since the A matrix contains the
characteristic polynomial information only. Now, observing
that the transfer function from the friction input to the
velocity as an output is
Yn(s) / l(s) = [0 0 . 0 1] T12 (sIA)1 b
(this holds since y(s)=T12x(s), and x(s)=(sIA)lbl(s) ),
with b' = [0 0 . 1], it is clear that the transfer
function zeroes are defined by the bottom row of T12, that
is, the vector P. Thus, varying the zeroes only would
change the P vector only (in such a way that poles of the
closed loop system are unaffected), X0 would be unchanged,
and the four (or two) switching conditions alone would be
affected. The control form representation clarifies this
relationship, which suggests an interesting line for further
research.
Illustrative Examples
Example IV1: TwoDimensional (2D) System with Sliding
Friction Limit Cycle
This is the same as Example I1 of Chapter I; the
piecewise linear method is the same as standard phase plane
methods for this 2D case. The state variable model of the
system is
(4.16) dyl/dt = y2
dy2/dt = (K/J) yl (B/J) Y2 Lf/J
where Lf is the friction torque, yl is position, and y2 is
velocity. The analysis performed below assumes K, J > 0,
although B can be negative (the other cases can be solved by
similar methods, but are left out for brevity). Chapter I
presents more information about the system, including
results of a previous solution by phase plane methods.
A secondorder system cannot possibly have a limit
cycle with sticking, since once the system sticks, it is at
an equilibrium point (there is insufficient torque to break
free, and no controller integrator to ramp up). Therefore,
we need only examine the necessary conditions for the
sliding case. In addition, it will be assumed for this
example that the limit cycle is symmetric (this property is
proved for 2D systems in Chapter V). Applying the equations
developed in this chapter, the initial condition can be
found from equations (4.5) (which (4.15) simplifies to using
the symmetry assumption) and (4.6):
AT
(4.17) x(T1) = eAT1 (x 1) + 1 = 
(4.18) y2(T1) = 1' x(T1) = x2(T1) = 0
where use is made of the fact that 0' = [ 0 1 ] (physical
state model already in canonical form), and x(T1) = X0 (by
symmetry).
Making the assumption that K is large enough (or B is
small enough) that the system poles are complex, the matrix
exponential for this 2D system can be evaluated as
(4.19) eA = et cospt+(a/p)sinpt (1/p)sinpt
(K/JP)sinpt cos3t(a/p)sin3t
where a = B/(2J), 3 = (1/2) { (4K/J) (B/J)2 ).
Letting 0 = [a 0]' (since x2=0), and noting 1 = [LC/K 0]',
equations (4.17) and (4.18) become:
aYT
(4.20) x2(T1) = 0 = [a (LC/K)] eT 1 (K/JP)sin3T1
which implies T1 = 7/3 (note a > (LC/K) since otherwise x0
is an equilibrium point).
Theoretically, T1 can also be any integer multiple of
at
7/3. However, noting that x2(t) = velocity = C e sin/t,
any multiple greater than one would entail sign changes in
the velocity for t
This is an example of a solution to the necessary conditions
that is not a valid limit cycle; see discussion on
sufficient conditions below.
The other component of x0 can be evaluated from (4.17)
as
(4.21) x1(T1) = a
arT
= [a (Lc/K)] eT l(cospT1+(a/P)sinPT1) + (LC/K)
= [a (LC/K)] ea7/8 cosr + (LC/K)
Therefore
(4.22) a + (LC/K) = [a (LC/K)] ear/
The initial condition can now be evaluated as
(4.23) a = (Lc/K) ((/+1)/(41)), P = ear/
which is the result presented in Chapter I. Since x0 is in
region IV, a must be greater than zero. This implies j>1,
requiring a<0, or B<0.
As the solution from Chapter I predicted, a sliding
type of limit cycle exists for every case with negative
damping (B<0) (if the position feedback gain K is
sufficiently large to give the linear system a complex pair
of poles). The initial conditions can be expressed in terms
of the system parameters by the equations from Chapter I:
(1.2) X10 = (LC/K) [(x+1)/(g1)], x20 = 0
where
(1.3) p = exp[Br/(2JP) ], = P [4K/J B2/J2]
Refer to Figure I4 for a plot of amplitude (which equals
x10) as a function of gain K for the case J=LC=1, B=l.
Example IV2: ThreeDimensional (3D) System with Both
Sliding and Sticking Limit Cycles
The second example system to be considered (same as
Example I2) is defined by the differential equations
(1.4) dy/dt = A y + b L
where
(1.5) A = 0 1 0
0 0 1
K1 K2 B
with b' = [0 0 1], Lf = friction torque, and where the
state variables are y2 = position, y3 = velocity, and
yl = compensator integrator. Figure I6 shows a block
diagram of this system. The results presented in Chapter I
are now to be derived, based on the equations developed in
this chapter.
Case I: One Stable Pole and One Imaginary Pole Pair (Only
Sticking Limit Cycle Exists)
As in the example in Chapter I, let us first set
B = K = K = 1, and set the sliding friction torque to 1.0
and sticky friction (breakaway torque) to 1.2.
The equations (4.2a/b, 4.6, 4.8) defining the switching
periods and initial state (4.14) can be set up and solved
for this specific case, in order to demonstrate the
analytical calculation of the limit cycle trajectory. The
eigenvalues of this system are 1 and jl. The matrix
exponential of A is required for the equations, and can be
found (by, for example, diagonalizing A) as
AT T T
(4.24) 2eAT e +sinT+cosT 2sinT e +sinTcosT
T T
e sinT+cosT 2cosT e +sinT+cosT
T sT
e sinTcosT 2sinT e sinT+cosT
The equations to be applied are
AT
(4.25) x = e AT1 ( 1) + 1
AT
(4.26) x2 = es 2 x x0
where a symmetric limit cycle was assumed (breakaway point
at t = T1 + T2 is at x0), and therefore the simpler
equations were used in place of (4.14). Note that A is the
s
A matrix with the bottom row set to zero, hence the matrix
exponential in the sticking region is
(4.27) eAsT2 = 1 T2 T22/2
0 1 T2
0 0 1
In addition, the breakaway condition that the acceleration
torque equal the static friction is required:
(4.28) 1'Ax0 = dx3(0)/dt = LS
Noting that the translation vector 1' = [ 1 0 0 ], letting
0 = [ a b 0 ] (since x30 must be zero), and letting
xI = [ c d 0 ], we obtain
(4.29) c 1 = (1/2)(al)(eT + sinT + cosT) + bsinT
d = (1/2)(a1)(eT sinT + cosT) + bcosT
0 = (1/2)(a1)(e sinT cosT) bsinT
from application of equation (4.25) and the definition of
AT
the matrix exponential e,
(4.30) c + d T2 = a
d = b
from application of equation (4.26) and the definition of
A T
e s 2, and
(4.31) a + b = L
S
from application of equation (4.28). By eliminating c and d
from these equations and performing additional algebra, they
simplify to
(4.32) b = LS a
T
T2 = [(al)eT1 + (a+l)] / (LS a)
T
a [e 1 + sinT cosT ]
T
= eT1 + (21sf l)sinT1 cosT1
T
a [e 1 sinT1 cosT1 2]
= e T1 (2 LS l)cosT1 sinT1 2LS
The first two equations define b and T2 in terms of a and
TI, while b and T2 were eliminated from the last two, which
can be used to simultaneously solve for a and T1.
Eliminating a and after some algebra, the last two equations
form a nonlinear algebraic equation, whose zeroes are
potential limit cycle solutions:
(4.33) f(T1) = (eT 1 l)(l+cosT1) (eT1 + l)sinT1 = 0
The plot of this function in Figure (IV1) and analysis
shows that there are zeroes at odd multiples of r, and also
near 37/2 (+ 2n7, n = 0, 1, ...). The solution at r results
in a valid limit cycle solution where x = [ 1 0.2 0 ],
0
T1 = r, and T2 = 10, yielding an overall limit cycle period
of approx. 26.3 seconds, matching the simulation results.
The second potential solution is at T1 = 4.73 (approx. 37/2,
obtained by numerical solution of f(T1)=0), and is invalid,
since it results in a negative value for T2 (of 12.214).
Additional solutions are also invalid, for the same
reason as in Example IV1. Examination of the matrix
'
__ _  ,^.
o
i
I
oo
000
O rr m
> 0 o 0 0
u
I0
U
HYA
o4
o
41
, 
I
M U
*1
Cn C
11C ti
exponential and the resulting formula for the velocity state
shows that the velocity changes sign approximately every 7
seconds. Therefore, the trajectory using these longer
values for T1 leaves region II before the switching time, so
it is invalid (see next section).
Therefore, the solution of the piecewiselinear
equations results in a valid limit cycle solution that
matches simulation results. Chapter I contains plots of the
limit cycle trajectory.
Case II: One Stable Pole and One Unstable Pole Pair (Both
Sticking and Sliding Limit Cycles Exist)
If the same example system is used with B = 0.9, and
K1 = K2 = 1 still (case II in Chapter I), the closedloop
eigenvalues (of the linear portion of the system) move into
the righthalf plane, to .026 j 1.024, while the other
pole is at .9524. As stated in Chapter I, there is still a
sticking limit cycle solution close to that of Case I
(Figures I7 through 19, confirmed by simulation). A
similar analysis to that presented above can'be used for
this case, so the details are omitted.
Note, however, that numerical methods may be required
to evaluate the matrix exponential (due to the tediousness
of the calculations; case I was a particularly simple form);
this forces the use of iterative methods to calculate the
limit cycle parameters (see example of this method below).
The second limit cycle solution, as discussed in
Chapter I, is of the sliding type. Simulation having shown
one with a halfperiod T1 approximately equal to 3.1
seconds, an iterative numerical solution of the equations
for the symmetric, sliding case was performed, as follows.
The equation for this case is obtained from (4.5),
where again the assumption is made that the limit cycle is
symmetric:
AT
(4.34) xl = eAT ( 0 1) + 1
= xO
which can be solved as
(4.35) xO = (I + eAT1)1 (I eAT) 1
A computer program calculated values for the initial
condition xg, given trial values of T1, and the process was
repeated until the velocity initial condition XO(3) = 0.
This was not difficult, since an approximate starting
condition was available from the simulation results. Given
a T1 that resulted in a zero initial velocity, the entire
initial condition could be defined from (4.35).
This yielded a solution at T1 = 3.1455 seconds, and
20 = yO = [0.283 12.4 0.]. The solution is valid since
the breakaway condition is also satisfied (torque at zero is
greater than LS). Figures 110 and I11 show some views of
the limit cycle for this case. An appendix contains
computer code that was used to perform this iterative
calculation. The code applies to the sliding case only, and
would have to be modified for the general case.
Exact (Necessary and Sufficient) Conditions
Theorem 41 presented some conditions that were
necessary for the existence of a simple limit cycle. That
they are not sufficient is seen by the results in the
Examples, where solutions were found that met all of the
necessary conditions, yet were not valid limit cycles (the
switching period T2 was less than zero, or velocity changed
sign during a switching period).
In order to expand the previous conditions into a set
that is also sufficient, the concept of a consistent
solution is useful:
Definition: A consistent solution is a solution of the
equations in Theorem 41 that meets the assumptions about
the region containing the trajectory during each switching
period. These assumptions are
(Sticking Case)
(1) Trajectory in region I, 0 <
(2) Trajectory in region III, T
(3) Trajectory in region II,
T1 + T2 < t < T1 + T2
(4) Trajectory in region III,
T1 + T2 + T3 5 t T1
(Pure Sliding Case)
(1) Trajectory in region I, 0 <
(2) Trajectory in region V, t =
(3) Trajectory in region II, T,
t < T1
1 t 2 T1 + T2
+ T3
+ T2 + T3 + T4
< t < T1 + T3
(4) Trajectory in region IV, t = T1 + T3
A check of the consistency of a solution is
straightforward, based on the definitions of the various
regions presented in previous chapters: the sign of the
velocity is examined for assumptions (1) and (3), and the
magnitude of the torque for (2) and (4). Given this
definition, a set of necessary and sufficient existence
conditions can now be presented:
Theorem 42: A simple limit cycle exists in the system
under study if and only if a consistent solution of the
Theorem 41 equations exists. In addition, this limit cycle
has the exact properties (switching periods, initial
condition, and trajectory) defined by the corresponding
solution to the equations.
Proof: (Necessity) Suppose a simple limit cycle exists
with the given properties. Theorem 41 (more exactly, the
derivation in the first part of this chapter) demonstrated
by analysis of this limit cycle that the listed equations in
fact apply to the trajectory.
In addition, this derivation, along with the assumption
that the limit cycle is simple, shows that the solution must
meet the consistency conditions, i.e., it must be in the
appropriate region of the state space during each switching
period. Therefore this condition is also necessary.
(Sufficiency) Now suppose a consistent solution of the
equations exists. A limit cycle can be generated, with the
same properties as defined in the solution, by starting at
x0 and following around to the initial point again. At each
step of the process, the consistency conditions require the
trajectory to be within a given region, so that the system
equations for that region apply. The equations listed in
Theorem 41, derived from the application of those system
equations, then show that the trajectory reaches the next
switching point at the appropriate time and state.
Therefore, a limit cycle must exist with those exact
properties, and the proof is complete.
Although this theorem defines a set of exact
conditions, i.e., they are equivalent (necessary and
sufficient) to the existence of the limit cycle, it is
appropriate to ask if these conditions are minimal. In
other words, are there redundancies in the set of
conditions? Is there a simpler set of conditions that are
still exact? For example, it may be sufficient to merely
require switch periods greater than zero, thus eliminating
cases such as found in Example IV2. These questions are
open, and may be subjects of further research.
CHAPTER V
LIMIT CYCLE SYMMETRY AND STABILITY
This chapter uses the results of chapter IV and other
known facts about the friction nonlinearity to examine the
behavior and characteristics of friction limit cycles. The
model of the system developed previously is used to analyze
the stability of limit cycles predicted by the equations in
chapter IV. In addition, the relationships between the
oddness of the system differential equations, symmetry, and
uniqueness of limit cycles is explored.
Stability of Predicted Limit Cycles
Once a limit cycle has been predicted by the solution
of the equations derived in Chapter IV, it is natural to ask
about its relationship to the global phase portrait. Since
the piecewise linear model used here is exact, the exact
limit cycle period, amplitude, and trajectory can be found
immediately once an initial point on the limit cycle is
known. It is more difficult to determine stability and
other information about the phase portrait, however, since
it can require extensive simulation.
However, local stability of the limit cycle can be
determined by the standard method of linearization of the
periodic flow map at the initial point previously found. In
other words, an approximate model of the trajectories near
the limit cycle initial point can be determined by
71
linearization. This analysis provides information about the
asymptotic stability of a closed orbit solution as discussed
in Hirsch and Smale (1974), chapter 13 (if more information
on orbital stability analysis is of interest, this reference
provides a good description, and was quite useful in the
development of the results in this chapter).
Note that the model provided by this analysis is not a
continuoustime representation of trajectories near the
cycle; instead, this discretetime model shows the deviation
from the initial point and how it evolves after each cycle
period. Specifically, an initial sufficiently small
deviation from the initial point x0 is assumed (so the
trajectory starts at x +6x ) and the trajectory traced once
around the limit cycle. The new deviation from the initial
point when the trajectory completes the cycle is related to
the initial deviation by the equation
(5.1) 6xT = Y 6Sx + higher order terms
where 6xT is the new deviation from the initial point
(0+6xT is the point on the trajectory after one cycle), and
Y is the transition matrix, representing the firstorder
evolution of the deviations from the initial point over
time. Eigenvalues of the transition matrix determine
stability, while eigenvectors determine the phase portrait
near the limit cycle.
This analysis is attempted below; however, because of
the discontinuities in the system nonlinearity the analysis
is difficult. Therefore, the analysis is accomplished, as
in the derivation of Chapter IV, by solving one region at a
time, and then pasting the trajectories together. The first
case analyzed is that where the predicted limit cycle has no
sticking (sliding friction only). The case with sticking
during the cycle is then analyzed by building on the first
case.
Case 1: Limit Cycle with no Sticking
Assume an initial point x0 that meets the equations in
Chapter IV with no sticking. In this case, the first
switching period is equal to the constant T1 in chapter IV,
the second switching period is T3, and the constants T2 and
T4, representing the periods of stickiness, are zero.
A deviation 6x0 is assumed from the initial point of
the limit cycle. We wish to limit deviations, however, so
that the perturbed initial point is still contained in
region IV. This is necessary in order to apply the Chapter
IV equations, since the initial point to which the equations
apply is assumed to be in region IV. Limiting the deviation
in this fashion results in no loss of information about the
limit cycle stability, however; (nl) of the eigenvalues of
the periodic flow map can still be determined, while the nt
eigenvalue (that applies to deviations out of region IV, and
along the limit cycle) would be zero, as will be seen in the
examples below. Note that the derivative of the flow after
one orbital period has an eigenvalue of unity (Hirsch and
Smale, 1974, p. 277), but the Poincare map has a zero
eigenvalue (Hirsch and Smale, chapter 13, section 3).
In order to limit the deviation 6x such that the
initial point remains in region IV, 5x can be constructed
as
1 '
(5.2) 6x = T12 Ir 6y0r
where T12 is the composite linear transformation from
Chapter III used to convert between canonical and physical
state variables (y = T12 x), I is an n x (nl) matrix (Ir
is (n1) x n) used to force the nth element of the deviation
vector Sy0 to be zero:
(5.3) I = 1 0 0 ... 0
0 1 0 ... 0
0 ... 0 1 0
0 ... 0 0 1
0 ... 0 0 0
and, finally, Sy0r is the (nl)vector of deviations in
physical state variables. This construction forces the 6x0
vector to remain in region IV, since the change in initial
velocity (6Sy(n), the nth element of the physical state
vector) is zero. Although the derivation below will be done
in canonical variables, this construction (5.2) will be used
at the end to complete the calculation for physical
variables, and force the deviation to have zero initial
velocity.
Starting at the initial point x +6x the trajectory is
determined by the piecewise linear methods as in chapter IV.
The trajectory in region I is found using autonomous state
variables:
(5.4) x (t) = x(t) 1
a
so the trajectory is
At At
(5.5) x (t) = e x (0) = e [x0 + x
a a
At
x(t) = e [X0 + SX0 1] + 1
For a sufficiently small deviation 6x0, the trajectory
reaches region V (velocity goes to zero preparatory to
reversing direction), but not necessarily in exactly T1
seconds. Denoting the change in this period is ST, the
velocity goes to zero when the trajectory reaches
(5.6) x(T1 + 5T) = eA(T1+6T) [x + 5x 1 1
The quantity of interest is the deviation of this point from
the point where the original limit cycle reaches region V
AT
(5.7) x e 1 [X 1] + 1
(5.8) Sx = x(T1+6T) x
AST AT A(T +aT)
=(e I} e 1 [ ] + e(T1T) 6x
Note that the first term shows the dependence on ST, since
it is zero if ST is zero, while the second term shows the
dependence on S6x0
Q
Now that the exact value for Sx is determined, this
nonlinear formula must be linearized to determine local
behavior. One factor that causes difficulties is that 6T
depends on 6x0, which complicates the calculation for
(5.9) Y = (d(Sx )/d(5x )}Ix=
where Y is the desired transition matrix defined above, and
the derivative is evaluated for 6x = 0. The dependence of
Sx on both ST and 6x0 explicitly is handled using the chain
rule.
Note that since x is not a function of 6x or ST,
(5.10) [d/d(6x0)] (6x ) = [d/d(6x )] {x(T1+6T) x)
= [d/d(Sx0)] (x(T1+6T))
= [d/d(6x0)] (eA(T1+6T) [X0 + S 1])
where the term in x(T1+6T) equal to 1 was also dropped since
it does not depend on x 0. Using the chain rule, this
derivative is
(5.11) d(6x1)/d(6x0) = D(6X1)/D(6x0)
+ [D(6X1)/D(6T)] [d(6T)/d(6x0)]
where the symbol "D" indicates a partial derivative. The
first term is
(5.12) D(6x )/D(6x0) = eA(T1+6T)
AT
= eAT1
when evaluated at 6x0=0 (hence 6T=0). The partial
derivative in the second term can be evaluated as follows.
First define
AT
m = e 1 [X0 + &x 1]
Then, by using (5.10), the desired partial derivative is
(5.13) D(6x1)/D(6T) = [D/D(6T)] {eAST M)
= AST
A e m
AT
= Ae 1 [ 1]
when the equation for m is substituted back in, and the
derivative is evaluated at 6x0=0 (6T=0). The
differentiation step can be verified by expanding the
exponential into a power series, and evaluating term by
term.
The evaluation of the [d(6T)/d(6x0)] factor in the
second term requires the definition of a functional
relationship between 6T and 6x0. We know from the
derivation of x(T1 + 6T) above that the velocity goes to
zero at this point (trajectory reaches region V). As in the
analysis in Chapter IV, this can be stated as
(5.14) P' x(T1 + ST) = 0
(eA(T +6T)
= {e 1 [x + 6x0 1])
since P'1 = 0. This formula implicitly defines the
relationship between 6T and 6x in that we have a function
of the form f(6T,6x0) = 0. The implicit function theorem
(Rudin (1976), pp. 2238) can be used to evaluate the
desired derivative as
i
(5.15) d(6T)/d(6x0) = [Df/D(6T)]1 [Df/D(6x0)]
Evaluating these factors,
(5.16) Df/D(ST) = [D/D(6T)] {('(eA(T1+T) [0+6x011))
A6T
= [D/D(6T)] {(' eAST m)
A6T
= A e A
with m defined as above (to see this, again expand
exponential into series and evaluate term by term). The
second factor is
(5.17) Df/D(6x0) = [D/D(6x0)] {('(eA(T1+6T) [x0+60]))
= [D/D(6x0)] {(' eA(T1+6T) 6x0)
= eA(T1+6T)
The derivative from the implicit function theorem, (5.15),
can now be evaluated as
AT 1 AT
(5.18) d(6T)/d(6x0) = ({'A eAT1 [OI]} ) ('eAT 1
where the factors are again evaluated at 6x =0 (6T=0).
Finally, combining (5.12,13,18) into (5.11), we get the
total derivative representing the linearized periodic map
(5.19) Y = (d(6xl)/d(6x0)) 16x=o
AT AT AT 1 AT
= eA1 (A eAT1 ['X1]) ('A eAT1 [O i{e A1)
Note that the factor to be inverted is a scalar. This
expression can be simplified by noting that
AT
(5.20) i = eAT1 [x 1] + 1
so substituting into (5.19) gives
AT 1 AT
(5.21) Y = eAT1 (A [X11]) (P'A [X1 i)1 ('eAT1)
= eAT1 (A [xl1] 'e AT1) / ({'A [x1])
The linearized flow map for the entire cycle must now
be constructed, using this equation (5.19) for the
linearized flow map for a halfcycle. Note, first of all,
that the linear transformation representing the linearized
flow map for the cycle is the product of the'linear
transformations of the two halves. This is a consequence of
the chain rule, since
(5.22) 6xT = f (6x ) = (f2(6x0))
d(6xT)/d(6xo) = (d(6)/d(6x1)) (d(6x1)/d(6x ))
The second of these two linear transformations has already
been derived; the first must be found and the linear map for
the entire cycle found by multiplication.
Assume first of all that the limit cycle is symmetric,
so that the point at which the limit cycle velocity goes to
zero is exactly onehalf period after the initial point, at
X0. In this case, the halfperiod is equal to the constant
T1 and T3 equals T l
Therefore, the linearization in which we are interested
is completely determined by the first half of the limit
cycle trajectory, because of the odd symmetry. To show
this, assume that the same deviation from the point xl is
taken, except with a change of sign, and trace the
trajectory from this point x 6x Noting that this initial
point equals (x06x0), and using the odd symmetry of the
phase portrait, it is clear that the trajectory returns to
region IV at the point (x 6x ) (since the trajectory from
x0+6x0 goes to x1+6x ). Thus the resulting deviation on the
second half of the cycle is
(5.23) 6x = (x26x ) X
= 6X1
since xl = x0 Linearizing this second leg, using the
initial deviation of 6x ,
(5.24) 6x = Z (6x ) + higher order terms
2
= 6x
= Y 6x + h.o.t.
by the analysis of the first halfcycle. Since this holds
true for any deviation 6x0, we must have Y = Z. Therefore,
the transition matrix for the whole cycle is the square of
the transition matrix for the first halfcycle, Y2
Now remove the assumption of symmetry. Assume that the
limit cycle is not symmetric (but still involves only
sliding, i.e. case 1, no sticking); in this case the
linearized flow map must be examined on the second half of
the cycle also. The same derivation as above (equations
(5.4) through (5.8)) can be performed, but starting at the
initial point x1+6x1, in order to determine SxT in terms of
6x1. The trajectory in region II is found using autonomous
state variables:
(5.25) x (t) = x(t) + 1
a
so the trajectory is
(5.26) x(t) = eA(tT) [ + x1 + 1] 1
For a sufficiently small deviation 6x1, the velocity goes to
zero when the trajectory reaches
(5.27) x(T + 6T3) = eA(T3+6T3) [x + 6x + 1 
The deviation of this point from the point where the
original limit cycle reaches region IV is
(5.28) 6xT = x(T+ST3) xT
= (eAT3 I) eA3 [x + 1] + eA(T3+T3x
Note that by comparing this formula to (5.8), and stating
(5.8) as a function with parameters
(5.29) 6x= (eAST I) eAT1 [x 1] + eA(T1 T)
= f(6x0, 6T; T1, x01)
then we see
(5.30) 6x = f(6xl, ST3; T, x +1)
The linearization for this function has already been
derived, so the linear map for the second half of the cycle
is found by substituting the correct parameters from (5.30)
into (5.19), so that
(5.31) 6ST = Y(T3, x1+1) 6x, + h. o. t.
(5.32) Y(T3, l+1) = d(6XT)/d(6X1) 6x=0
= eAT3 (A eAT3 [+]) (~'A eAT3 [+)1 ('eAT3)
By the argument given previously in (5.22), the total
linearized flow map is the product
(5.33) Ycompos = Y(T3' +1) Y(T1' 0)
6T = Y cx + h. o. t.
T compos 0
Note that this reduces to the symmetric case when T3 = T1
and x = 0:
(5.34) Ycompos Y(T x+1) Y(T, xO1)
= Y2(T1, 1)
since Y is an even function of its second variable.
The final step in the derivation for this case is to
convert back to physical variables, while at the same time
forcing the initial deviation to the desired region. By
following the reasoning given at the start of the analysis
for this case (equation (5.2)), the desired matrix for the
linearized flow map of the entire cycle is
(5.35) Y compos = Y(T' +1) Y )
rcompos r 3' l1 O ( 0
1
Yr(T, x 1) = Ir T 2 Y(T x01) T Ir
(5.36) 6YTr = Yrcompos 60r + higher order terms
where 6yTr is the (nl)vector representing the deviation in
physical variables after going around one cycle (with
velocity = 6y left out since it is zero), 6y0r is the
initial deviation in the same variables (also has (nl)
elements, no deviation in initial velocity allowed), and Yr'
Y are the matrices representing the linearized flow
rcompos
maps derived above, but in physical variables and reduced to
order (n1). Note that the latter matrix is (nl)x(n1), so
its (n1) eigenvalues are those desired for the local
stability analysis (as discussed above, the nth eigenvalue
is always 0).
Case 2: Limit Cycle with Sticking
To analyze this case, it is necessary to derive the
linearized flow map for the portions of the limit cycle in
the sticky region, i.e., region III. These linear maps can
then be combined with the maps previously derived for the
sliding regions to form the composite map for the cycle.
Assume an initial point x3 that meets the equations in
Chapter IV, so that this is the point at which the system
comes to rest in region III (this is the fourth switch point
in the cycle). In this case, the period the system sticks
in region III is equal to the constant T4 in Chapter IV.
This initial point must be used, instead of the point x at
the point of impending motion, in order to derive the
behavior within region III.
A deviation 6x3 is assumed from the initial point of
the limit cycle. This deviation is of course constrained as
before, so that the resulting initial point remains inside
region III (velocity = 0). Starting at the initial point
x3+6x3, the trajectory is determined by the piecewise linear
methods as in chapter IV.
The trajectory in region III is
(5.37) x(t) = eAt x(0) = eAt [X3 + 6x 3
For a sufficiently small deviation 6x3, the trajectory
reaches the condition of impending motion (acceleration
torque sufficient to break loose), but not necessarily in
exactly T4 seconds. Denoting the change in this period is
ST, impending motion is reached when
(5.38) x(T4 + ST) = eAs(T4+6T) [x3 + 6x3]
The quantity of interest is the deviation of this point from
the point where the original limit cycle breaks free
AT
(5.39) x es 4 x
(5.40) 6x = x(T4+6T) 4
A 6T AT A (T +6T)
= {es I) es 4 x + esT4T)
The differentiation of this function to determine the
firstorder dependence of 6x on 6x3 is performed as before;
many of the steps are similar and thus omitted. The desired
derivative is
(5.41) Ys = d(6x4)/d(6x3) 1=0
where Y is the desired linearized map in the stiction
region. Note that since x4 is not a function of 6x3 or 6T,
(5.42) [d/d(6x3)] (6x4} = [d/d(6x3)] {x(T4+6T) x4)
= [d/d(6x3)] (x(T4+6T)}
= [d/d(6x3)] {eAs(T4+6T) [x3 + 6X3])
The dependence of 6x on both ST and 6x3 explicitly is
handled using the chain rule, as before:
(5.43) d(6x4)/d(6x3) = D(6x4)/D(6x3)
+ [D(6x4)/D(6T)] [d(6T)/d(6x3)]
where "D" is again the symbol for partial derivative. The
first term is
(5.44) D(6x4)/D(6x3) = e s
when evaluated at 6x = 6T = 0. The partial derivative in
the second term is
A 6T
(5.45) D(6x4)/D(6T) = [D/D(6T)] {e s m)
A T
= A e s 4 [x3]
where m is a temporary vector independent of 6T, as in the
previous derivation.
The evaluation of the [d(6T)/d(6x3)] factor in the
second term is handled as before by relating the two
variables with an implicit function, and using the implicit
function theorem. The implicit relation in this case is
determined by the breakaway condition (i.e., impending
motion); as in the analysis in Chapter IV, this can be
stated as
(5.46) P' A x(T4+6T) = LS
or, equivalently,
0' A eAs(T4+6T) [x3 + 6x3] + LS = 0
Thus we have a function of the form f(6T, 6x ) = 0. The
derivative is then
(5.47) d(6T)/d(63) = [Df/D(6T)]1 [Df/D(6X3)]
Evaluating these factors,
(5.48) Df/D(6T) = [D/D(6T)] ({' A e s6T m)
A 6T
= e' AA es m
s
The second factor is
(5.49) Df/D(6x3) = [D/D(6x3)] ({' A eAs(T4+6T) 6 x
= A eAs(T4+6T)
so that
AT 1 AT
(5.50) d(6T)/d(6x ) = {('AA eAs 4 x3) {('A eAs 4)
where the factors are again evaluated at 6x =6T=0. Finally,
combining (5.44,45,50) into (5.43), we get the total
derivative representing the linearized periodic map
(5.51) Ys(T4, x3) = d(6x4)/d(6x3) 6x=0
AT AT AT 1 AT
= e s 4 (As e s 4 x3) ('AAs e s 4 x ) ('A e s 4)
Note that the factor to be inverted is a scalar. This
expression can be simplified by noting that
AT
(5.52) x4 = e s x3
so substituting into (5.51) gives
(5.53) Ys(T4x3) = eAsT4 (As x4)('AAs 4)('A e s4)
AT AT
= e s 4 (As x4 'A e s 4) / {('A A x4)
In order to obtain the linear map for the entire cycle
for the case with sticking during the cycle, we can use the
chain rule as before:
(5.54) Ycompos = d( )d(/d(6 2)) (d(6x2)/d(6x ))
x (d(6x )/d(6x0 ) (d(6x4)/d(x3 )
where the last factor was just derived, the first and third
are obtained from the sliding case previously analyzed, and
the second factor can be found from the fourth by symmetry,
as shown below. Note that 6x0 = 6x, since both are the
deviations at the point of impending motion in region III;
this equivalence justifies the juxtaposition of the last two
factors when the chain rule is applied.
Note also that in (5.54) the order of the factors is
different from that used in the sliding friction case; the
initial deviation for the case with sticking is taken as 6x3
at the point x where the system comes to rest, rather than
at the point x where motion is impending. This order is
required since any (small) deviations Sx3 are allowed that
keep the initial velocity zero (hence 6x3 can be constructed
as in (5.2)). On the other hand, certain initial
perturbations would not be allowed at the point x ; any
perturbations at this point of impending motion that altered
the forcing torque would either cause motion (so initial
point no longer in region III), or would mean motion is no
longer impending. This problem does not arise in the
sliding friction case, but in the sticking case there are
two constraints on x0: zero velocity and forcing torque at
breakaway. It is therefore much more convenient to take the
initial point to be x3 for the stability analysis; all
constraints are then taken care of in the derivation
automatically.
The second factor in (5.54) can be found by symmetry,
as follows: taking a deviation from the point x~ at which
the system comes to rest in region III of Sxl, we can derive
by the same steps as above that
(5.55) 6x2 2=(T2+6T) x
=(es I) eAsT2 x + eAs(T2+T)65
1 1
The similarity of this equation to (5.40) can be used to
show that the desired linearization is
(5.56) Ys(T2, 1) = d(6x2)/d(6 ) 6X=
AT AT AT 1 AT
= es 2 (As es 2 xl} ('AAs e s 2 xl)1 'A eAsT2}
Of course, as in the sliding case, this matrix must be
modified to transform to physical variables, and reduced to
order (nl).
The linear transformation for the entire, fourpart
cycle is then
(5.57) Yrcompos = Y(T3' 2+1) Y (T x )
rcompos r 2 rs 2'1
x Yr (T1,Xi) Yrs(T4'x3)
1
Y (T1' 0x1) = I T2 Y(T 1) T2 I
12 r Y1' 0 12 I
1
Y(T2' l) =Ir T 2 Ys(T2' l) T12 I
(5.58) STr = Yrcompos XOr + higher order terms
where, as before, the Syr vectors are the (nl) elements of
the initial and final deviation, respectively, and Y
rcompos
is the desired linear transformation that determines local
stability of the limit cycle.
The results of the stability analysis of this chapter
can be summarized as follows:
Theorem 51: Given the existence of a friction limit cycle,
as predicted by the solution of the nonlinear algebraic
equations derived in Chapter IV (Theorem 41), the local
asymptotic stability of the periodic orbit is determined by
the eigenvalues of the matrix:
(Case 1: Sliding only)
(5.35) Y compos = Y (T +) Yr(T )
rcompos r 3' 1 r(1' 0
(Case 2: Sliding plus sticking)
(5.57) Yrcompos = Y(T3' x2+1) Y rs(T2, 1)
rcompos r 3' 2 rs2' x1
x Y (T1, xO1) Yrs(T, x3)
where
1
Yr(T XO1) = I1 T12 Y(TI, X1) T12 I1
1i
Yrs(T x) = I T2 Y (T, x ) T2 Ir
the x 's are the switching points, the T 's are the
corresponding switching periods, Ir is the (nl) x n matrix
defined in (5.3), T12 is the transformation from physical to
canonical variables defined in Chapter III, and the matrix
functions Y and Y are
s
(5.19) Y(T1, x01) = d(6Xl)/d(6x0) 6x=0
SeAT1 (A eAT1 [X0i]) ('A eAT1 [x 1]){('eAT1)
(5.56) Ys(T2, X1) = d(6x2)/d(6x1) 6x=0
AT AT AT 1 AT
= es 2 (As e s 2 X1) ({'AA eAs 2 X}) ({'A es 2
Before this theorem is proved, some terms must be
defined more rigorously.
Definition: The flow of the differential equation system is
the map (t,x0):'(t,x0), where 9(t,x0) = x(t)'when the
initial condition x(0) = x .
0'
Definition: The Poincare map (Pmap) is the map 6x :g(6x ),
where g(6x0) = x(t) x0, 6x0 = x(0) x, and tl is the
first time at which the trajectory again crosses the local
section around the periodic orbit at x .
0
Proof: The proof follows from the theory in Chapter 13 of
Hirsch and Smale (1974) and the calculations above, except
for one difficulty to be overcome. The reference theory is
stated for the case where f(.) is continuously
differentiable (so that the flow < is continuous). However,
as will be shown below, all that is required to obtain the
desired result is that the flow be a continuous function of
the deviation about the initial point, x + Sx so that
certain limits can be taken, and, in addition, that the
Pmap is differentiable, so that the derivative may be used
to check asymptotic stability.
As a preliminary step, then, these properties must be
shown. The flow is the integral of a realvalued function
of time (see equation 5.69 below), f(x(.)), so it is a
continuous function of time (Rudin, (1976), Theorem 6.20).
It is also clear that the flow is continuous in x within
each piecewise linear region, since the function f(.) is at
least C1 there (and Hirsch and Smale, Chapter 15, shows this
is sufficient for continuity of the flow). Now, taking the
initial condition as a sufficiently small deviation from the
initial point x0' the flow at any finite time is the
composition of a finite number of continuous functions,
formed from the pieces of the trajectory in each region it
traverses. This composition can be formed since the
trajectory is continuous across the discontinuities in f(.);
that is, x(t+) = x(t), where t is a switching time.
Therefore the flow is a continuous function of this initial
condition (Rudin, (1976), Theorem 4.7, the composition of
continuous functions is continuous), for sufficiently small
deviations.
The same approach is used to demonstrate the
differentiability of the Pmap. This map is simply the
composition of the flows in each of the regions, so is
continuous by the previous argument. But, in fact, the
flows in each region are continuously differentiable, so
their composition is also differentiable.
Using these properties, the rest of the proof proceeds
as in Hirsch and Smale (1974). In chapter 13, section 3, it
is shown that if the Pmap has x0 as a sink, then the orbit
is asymptotically stable (using the property of continuity
of the flow to show that the flow converges uniformly to the
periodic orbit if the deviations after each cycle also
converge). But the Pmap is a discretetime dynamical
system, so it is a sink if the derivative of this map has
all eigenvalues less than one in magnitude (which requires
the property of differentiability of the Pmap). The
calculations of this derivative have been presented above,
and show that the matrix Y is in fact the derivative
rcompos
of the Pmap at x Therefore the (nl) eigenvalues of this
matrix determine if the orbit is asymptotically stable, and
the proof is complete.
To summarize the above stability results, it has been
shown that the perturbations from the periodic orbit evolve
according to the equation
(5.58) yTr = Yrcompos 65Or + higher order terms
where 6YTr is the (nl)vector representing the deviation in
physical variables after going around one cycle (with
velocity = Syn left out since it is zero), and Sy0r is the
initial deviation in the same variables (also has (nl)
elements, no deviation in initial velocity allowed). The
eigenvalues of the matrix Yrcompos' therefore, indicate the
stability of the periodic orbit (i.e., if they are all less
than one in magnitude, the orbit is asymptotically stable,
and in fact a periodic attractor).
In addition, it has been shown that the results on
orbital stability in Hirsch and Smale (1974), chapter 13,
can be extended to the case where f(.) is only piecewise
linear, rather than C1.
Stability Calculations for Example Problems
The limit cycle orbital stability equations derived
above can be applied to the example system discussed in
Chapters I and IV (example 12, case 2, a 3D system with
both a sliding and stiction limit cycle). The example
system is defined, as in the example, by
(5.59) dy/dt = A y + b Lf
where
A = 0 1 0
0 0 1
K1 K2 B
with b' = [0 0 1], Lf = friction torque, and where the
state variables are y2 = position, y3 = velocity, and
y = compensator integrator. The transformation matrices T1
and T2 are identity matrices, and the row vector 9' equals
[ 0 1]. The particular case has B = 0.9, K = K2 = 1, with
closedloop poles at .026 j 1.024 and .9524. The
stiction limit cycle solution is at approximately
x = [0.98 0.22 0], with T1 at approximately 7. There is
also a solution of the equations for the symmetric, sliding
case, at T1 = 3.1455 seconds, and xO = [0.283 12.4 0].
Example Vl, Part 1: Stability of Sliding Limit Cycle
For this case, equations (5.35, 5.19) from Theorem 51
apply. Using the parameters defined above, we find that
AT
(5.60) Ae 1 [x 1] = [12.41 .001 13.13]'
so that
P' AeAT1 [x0
The matrix exponential
algorithm was used)
 1] = 13.13
is (approximately,
AT K
(5.61) eAT = .5314
.5227
.5830
Therefore, the matrix Y required
flow map is
.1126
1.054
.0603
to define
a numerical
.5227
.5830
.5294
the linearized
(5.62) Y = .0197 .0556 .0223
.5227 1.054 .583
.0002 .0000 .0002
Although there are unavoidable numerical errors in this
calculation, the true Y matrix apparently does have a zero
eigenvalue, as expected. The upper left block submatrix is
equal to Yr' and the composite stability matrix is
(5.63) Y = .02944 .05749
.5407 1.14
This matrix has eigenvalues of .002125 (with corresponding
eigenvector at [2.104 1.]' ), and 1.168 (with eigenvector
[.05 1.]' ). Thus it has one stable and one unstable mode
(note this is a discretetime system), and is a saddle
point. A simulation of this system was used to verify this,
where the xl and x2 states were sampled whenever the
velocity x3 went through zero (near x0). The plot in Figure
V1 shows the point x0, the eigenvectors of the stable and
unstable modes, and typical trajectories near x The
stable mode is very fast so that simulated trajectories
quickly jump onto the unstable eigenvector, then move along
it. The simulation experimentally verifies the stability
predictions.
Example Vl, Part 2: Stability of Sticking Limit Cycle
Similar calculations can be performed for the other
limit cycle that exists for this system, except that
equation (5.57) from Theorem 51 governs this case. Using
the initial condition defined above, with T1 = v and T 9
(as approximate values), and xl = x(T1) = [1. 0.22 0.]',
and recalling that the symmetry defines the other switching
points and periods, we find after substitution in (5.56)
(5.64) Y = 0 1 9.9
0 1 9
0 0 1
with Y being the upper left block submatrix. The
stability matrix for the sliding portions of the limit cycle
is also required, and the upper left block of Y is
