Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Periodic solutions of flexible systems under discontinuous control
(USC Thesis Other)
Periodic solutions of flexible systems under discontinuous control
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
PERIODIC SOLUTIONS OF FLEXIBLE SYSTEMS
UNDER DISCONTINUOUS CONTROL
by
Michael K. Borre
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MECHANICAL ENGINEERING)
May 2011
Copyright 2011 Michael K. Borre
ii
Dedication
To my parents
iii
Acknowledgements
I wish to express my deepest gratitude to my research advisor, Professor Henryk Flasher.
His patience and continued interest were instrumental in guiding me throughout my
research efforts at USC. Dr. Flashner also provided countless hours of consultation
outside of normal office hours, which helped make this work possible. I also wish to
thank Professor Ben Yang for serving on my dissertation committee, and for his help in
deepening my understanding of vibrations in continuous systems. I would also like to
thank Professor Petros Ioannou for serving on my dissertation committee.
Above all I would like to thank my parents Patricia and Jim, for their constant
encouragement and unwavering support during this extended period of time.
iv
Table of Contents
Dedication ii
Acknowledgements iii
List of Tables vi
List of Figures vii
Abstract xii
Chapter 1: Introduction 1
Chapter 2: Problem Description 6
2.1 Relay Feedback Control of Flexible Structures...........................................6
2.2 Relay Feedback Control System (RFCS) ....................................................7
2.3 Relay Types...............................................................................................11
Chapter 3: Background 16
3.1 Limit Cycle Analysis in Relay Feedback Control Systems.......................16
3.2 Frequency Domain Techniques.................................................................17
3.2.1 Describing Function Method.....................................................................20
3.2.2 Tsypkin’s Method......................................................................................25
3.3 Time Domain Methods ..............................................................................30
3.4 Methods to Calculate Region of Attraction ...............................................35
3.4.1 Simple Cell Mapping .................................................................................35
3.4.2 Poincaré-Like Simple Cell Mapping .........................................................39
Chapter 4: Switching Surface Cell Mapping Method (SSCM) 43
4.1 Introduction................................................................................................43
4.2 Formulation of the SSCM..........................................................................44
4.2.1 General Plant Dynamics ............................................................................44
4.2.2 Flexible Plant Dynamics............................................................................46
4.2.3 Poincaré Mapping Surface Definition .......................................................48
4.2.4 Poincaré Mapping Surface Intersections ...................................................49
4.2.5 Fixed Point Stability on a General Switching Surface...............................53
4.2.6 Poincaré Map Integration Stopping Criteria..............................................61
4.3 SSCM Procedure........................................................................................63
4.4 SSCM Summary........................................................................................67
v
Chapter 5: Applications 69
5.1 Two-Degree-of-Freedom System Under Relay Control with Hysteresis..70
5.1.1 Plant Description........................................................................................70
5.1.2 Limit Cycle Frequencies and Fixed Points ................................................72
5.1.3 Limit Cycle Fixed Point Region of Attraction Mapping ...........................78
5.1.4 Discussion of Results.................................................................................84
5.2 Multilevel Relay with Dead-Zone and Hysteresis.....................................85
5.2.1 Plant Description........................................................................................85
5.2.2 Limit Cycle Frequencies and Fixed Points ................................................87
5.2.3 Limit Cycle Fixed Point Region of Attraction Mapping ...........................94
5.2.4 Discussion of Results...............................................................................100
5.3 Slewing Beam with Nonlinear Feedback.................................................101
5.3.1 Plant Description......................................................................................101
5.3.2 Limit Cycle Frequencies and Fixed Points ..............................................106
5.3.3 Limit Cycle Fixed Point Region of Attraction Mapping .........................112
5.3.4 Stability Jacobian for Nonlinear Switching Surface................................117
5.3.5 Discussion of Results...............................................................................121
5.4 Nonlinear Plant Under Relay Feedback Control .....................................122
5.4.1 Plant Description......................................................................................122
5.4.2 Limit Cycle Frequencies and Fixed Points ..............................................124
5.4.3 Limit Cycle Fixed Point Region of Attraction Mapping .........................127
5.4.4 Discussion of Results...............................................................................132
Chapter 6: Concluding Remarks 133
References 135
Appendices
Appendix A: Block-Diagonal Forms of Switching Constraint Equation Matrices.....139
Appendix B: Switching Constraint Equations for Other Relay Types........................157
B.1 Relay with Hysteresis and Pure Time Delay ...........................................157
B.2 Relay with Linear-Zone (Saturation Nonlinearity)..................................160
vi
List of Tables
Table 5.1: Relay Reset Time Parameters (Figure 5.17).................................................... 93
Table 5.2: Relay Reset Time Parameters (Figure 5.27).................................................. 111
Table 5.3: Linear Plant Limit Cycle Periods and Fixed Points (Section 5.1.2).............. 125
Table 5.4: Nonlinear Plant Limit Cycle Periods and Fixed Points................................. 127
vii
List of Figures
Figure 2.1: General SISO Feedback Control with Piecewise Constant Nonlinearity......... 8
Figure 2.2: Piecewise Constant Nonlinearity...................................................................... 9
Figure 2.3: General Feedback Control with Relay Nonlinearity and Flexible Modes ..... 10
Figure 2.4: Linear Feedback Control with Relay Nonlinearity and Flexible Modes
Transfer Function Form .............................................................................................. 10
Figure 2.5: Ideal Relay (D=1)........................................................................................... 11
Figure 2.6: Relay with Hysteresis (D=1, δ=.5)................................................................. 12
Figure 2.7: Relay with Dead-Zone (D=1, δ=.5)................................................................ 13
Figure 2.8: Relay with Dead-Zone and Hysteresis (D=1, ε=2, δ=.25) ............................. 14
Figure 2.9: Pulse-Width Pulse-Frequency (PWPF) Modulator ........................................ 14
Figure 2.10: Multilevel Relay with Dead-Zone and Hysteresis ....................................... 14
Figure 2.11: Relay with Linear-Zone (Saturation Nonlinearity) ...................................... 15
Figure 3.1: Relay Feedback Control of Linear Flexible System with r = 0 ..................... 16
Figure 3.2: Equivalent System Used for Limit Cycle Analysis........................................ 17
Figure 3.3: Equivalent System in Describing Function Form .......................................... 20
Figure 3.4: Input/Output for Relay with Hysteresis Under Limit Cycle .......................... 22
Figure 3.5: Single Term, 2 Term, and 10 Term Fourier Series Representation of
Square-Wave Relay Output Under Limit Cycle ......................................................... 22
Figure 3.6: Modified Nyquist Plot of Linear System Frequency Response (red) and
Negative Inverse of Describing Function (blue) for Relay with Hysteresis............... 24
Figure 3.7: Verification of Limit Cycle Predicted by DFM ............................................. 24
Figure 3.8: Sinusoidal Relay Input and Square-Wave Output Under Limit Cycle .......... 26
Figure 3.9: Linear System Frequency Response (red) and Single Term Tsypkin Locus
(green) and Negative Inverse of Describing Function (blue) for Relay with
Hysteresis.................................................................................................................... 27
viii
Figure 3.10: Imaginary part of Frequency Response Im(G) (red) and Single Term
Tsypkin Locus Im(J) (green) and -1/N(a) (blue) for Relay with Hysteresis vs
Frequency.................................................................................................................... 28
Figure 3.11: Imaginary part of Frequency Response Im(G) (red) and Two Term
Tsypkin Locus Im(J) (green) and -1/N(a) (blue) for Relay with Hysteresis vs
Frequency.................................................................................................................... 29
Figure 3.12: Equivalent State-Space Time Domain Model.............................................. 30
Figure 3.13: Symmetric Unimodal Limit Cycle Waveform in RFS................................. 31
Figure 3.14: Frequency response using State-Space Locus f( ω) (red) and Two Term
Tsypkin Locus Im(J) (blue) and – Δ (green) for Relay with Hysteresis vs
Frequency.................................................................................................................... 33
Figure 3.15: 2D Cell Map Centers in State-Space Coordinates........................................ 36
Figure 3.16: 2D Cell Map Centers in Cell-Space Integer Coordinates............................. 37
Figure 3.17: System Trajectory For Integration Period “T” in State-Space Coordinates. 38
Figure 3.18: 2-State System Trajectory Between Poincaré Surface Intersections............ 41
Figure 4.1: SISO Regulator Control (r = 0) of Nonlinear Plant with General Feedback
Control Signal to Relay Nonlinearity ......................................................................... 45
Figure 4.2: SISO Regulator Control (r = 0) of Linear Plant with General Feedback
Control Signal to Relay Nonlinearity ......................................................................... 46
Figure 4.3: Mapping Surface Defined as System Output Causing Relay Switch in
Negative Direction...................................................................................................... 48
Figure 4.4: Mapping Surface Defined as Relay Input Causing Switch in Positive
Direction ..................................................................................................................... 52
Figure 4.5: Stable Limit Cycle Fixed Point Designations for Relay with Dead-Zone..... 53
Figure 4.6: Phase Plane View of Stable Limit Cycle Perturbations on Switching
Surfaces for Relay with Dead-Zone............................................................................ 54
Figure 4.7: Phase Plane View of Stable Limit Cycle Dynamics Between Switching
Surfaces for Relay with Dead-Zone............................................................................ 55
Figure 4.8: Convergence of Switching Surface Intersections to a Limit Cycle Fixed
Point ............................................................................................................................ 62
ix
Figure 4.9: Limit Cycle Region of Attraction Mapping Procedure.................................. 66
Figure 5.1: 4-State System with 1 Rigid-Body Mode and 1 Flexible Mode Under
Relay Control with Hysteresis and PD Feedback....................................................... 70
Figure 5.2: Rearranged Block Diagram with System Input r = 0..................................... 71
Figure 5.3: Symmetric Unimodal Limit Cycle Waveform in Relay Feedback System.... 72
Figure 5.4: Nyquist Plot Used to Estimate Limit Cycle Candidates................................. 73
Figure 5.5: System Frequency Locus Used to Determine Limit Cycle Candidates ......... 74
Figure 5.6: System Frequency Locus Intersections for Various Hysteresis Values......... 75
Figure 5.7: Relay Hysteresis vs Normalized Frequency Solutions ( ω
c
/ ω
n
)...................... 75
Figure 5.8: System Frequency Locus Used to Determine Limit Cycle Candidates ......... 76
Figure 5.9: Four Stable Limit Cycle Solutions Started From Initial Conditions.............. 78
Figure 5.10: End View of Point Mapping Surface and Zero State Startup Transient
Trajectory at 1st Intersection (location “A” in Figure 5.12)....................................... 79
Figure 5.11: Region of Attraction Map Showing Trajectory Intersection “a” (magenta
point) in Figure 5.12 Located in LC4 ROA (redXs): ROA for LC1=blueX,
LC2=greenX, LC3=cyanX, LC4=redX ...................................................................... 80
Figure 5.12: Simulation Showing Different Limit Cycles Excited by Pulse
Disturbance at Relay Output at 1000sec Intervals Starting at 500sec ........................ 81
Figure 5.13: Relay Input Trajectory Transient and Mapped State at First + to –Switch
After Pulse Disturbance: a) zero-state startup, b) pulse at 500 sec, c) pulse at 1500
sec ............................................................................................................................... 82
Figure 5.14: Map Showing Intersection “b” (magenta circle) in Figure 5.13 Located
on Boundary of LC2 ROA (greenXs): LC1=blueX, LC2=greenX, LC3=cyanX,
LC4=redX ................................................................................................................... 83
Figure 5.15: Map Showing Intersection “c” (magenta circle) in Figure 5.13 Located
on Boundary of LC1 ROA (blueXs): LC1=blueX, LC2=greenX, LC3=cyanX,
LC4=redX ................................................................................................................... 83
Figure 5.16: 4-State System with 1 Rigid-Body Mode and 1 Flexible Mode Under 2-
level Relay Control with Dead-Zone and Hysteresis.................................................. 85
x
Figure 5.17: Symmetric Unimodal Limit Cycle Waveform in 2-level Relay Feedback
System......................................................................................................................... 87
Figure 5.18: Describing Function Plot to Estimate Limit Cycle Frequency Range ......... 92
Figure 5.19: Five Stable Limit Cycle Solutions Started From Initial Conditions ............ 94
Figure 5.20: Region of Attraction Map at x
1
=0.2636 and Limit Cycle State
Trajectories for LC1=blueX, LC2=greenX, LC3=magentaX, LC4=cyanX,
LC5=redX ................................................................................................................... 95
Figure 5.21: Region of Attraction Map at x1=0.2636 ROA for LC1=blueX,
LC2=greenX, LC3=magentaX, LC4=cyanX, LC5=redX .......................................... 96
Figure 5.22: Simulation Showing LC3 Excited from Initial State Noted on ROA Map
in Figure 5.21.............................................................................................................. 96
Figure 5.23: Rigid-Body Mode Region of Attraction Map: ROA for LC1=blueX,
LC2=greenX, LC3=magentaX, LC4=cyanX, LC5=redX .......................................... 98
Figure 5.24: Simulation Showing Limit Cycles Excited from Rigid-Body Initial
Conditions Noted on ROA Map in Figure 5.23.......................................................... 99
Figure 5.25: Pinned-free Beam with Rigid-Body Mode and Flexible Modes Under
Relay Control with Dead-Zone and Hysteresis and Nonlinear Feedback ................ 101
Figure 5.26: Rearranged Block Diagram with System Input r = 0................................. 102
Figure 5.27: Symmetric Unimodal Limit Cycle Waveform in Relay Feedback System 107
Figure 5.28: Solutions of (5.51) with 0 to Estimate Limit Cycle Period Range ..... 110
Figure 5.29: Three Stable Limit Cycle Solutions Started from Fixed Point Initial
Conditions, System Output y(t) (blue), Relay Output u(y) (red) .............................. 111
Figure 5.30: ROA Map at Slice Through 3 Fixed Point Locations and Limit Cycle
State Trajectories for LC1=blueX, LC2=greenX, LC3 =redX ................................. 113
Figure 5.31: x2-x4 View of 2D Slice in Figure 5.30. ROA for LC1=blueX,
LC2=greenX, LC3 =redX ......................................................................................... 113
Figure 5.32: Simulation Showing LC2 Excited from Initial State Noted on ROA Map
in Figure 5.31............................................................................................................ 114
Figure 5.33: Rigid-Body Mode Region of Attraction Map: ROA for LC1=blueX,
LC2=greenX, LC3=redX .......................................................................................... 115
xi
Figure 5.34: Simulation Showing Limit Cycles Excited from Rigid-Body Initial
Conditions Noted on ROA Map in Figure 5.33........................................................ 116
Figure 5.35: 4-State System with 1 Rigid-Body Mode and 1 Nonlinear Flexible Mode
Under Relay Control with Hysteresis and PD Feedback.......................................... 122
Figure 5.36: Nonlinear Softening Spring Force Characteristics..................................... 123
Figure 5.37: Linear Plant Limit Cycle Waveforms (Section 5.1.2)................................ 126
Figure 5.38: Nonlinear Plant Limit Cycle Waveforms................................................... 128
Figure 5.39: Nonlinear Plant Rigid-Body Mode Region of Attraction Map: ROA for
LC1=red, LC2=cyan, LC3=magenta, LC4=green, LC5=blue, Divergence=black .. 129
Figure 5.40: Linear Plant Rigid-Body Mode Region of Attraction Map: ROA for
LC1=red, LC2=cyan, LC3=magenta, LC5=blue...................................................... 130
Figure 5.41: Nonlinear Plant Limit Cycles Excited by Initial Conditions in Figure
5.39............................................................................................................................ 131
Figure B.1: Limit Cycle Waveform with Pure Time Delay in Relay Feedback............. 157
Figure B.2: Relay with Linear-Zone (Saturation Nonlinearity)...................................... 160
Figure B.3: Limit Cycle Waveform Switching Criterion for a Relay with Linear-Zone 161
xii
Abstract
A method for calculating all periodic solutions and their domains of attraction for flexible
systems under nonlinear feedback control is presented. The systems considered consist of
mechanical systems with multiple flexible modes and a relay type controller coupled with
a linear or nonlinear control law operating in a feedback configuration. The proposed
approach includes three steps. First, limit cycle frequencies and periodic fixed points are
computed exactly, using a block-diagonal state-space modal representation of the plant
dynamics. Then the relay switching surface is chosen as the Poincaré mapping surface
and is discretized using the cell mapping method. Finally, the region of attraction for each
limit cycle is computed using the cell mapping algorithm and employing an error based
convergence criterion. The method is demonstrated on four application examples, each
with a different type of relay and/or control law to demonstrate the versatility and
accuracy of the method.
1
Chapter 1: Introduction
In many important engineering applications, the objective of mechanical control systems
is to generate periodic behavior (limit cycle) of given amplitude and frequency. For
example, in the case of spacecraft attitude control systems utilizing thruster actuators, the
goal is to generate a periodic motion of given amplitude and frequency, in order to
minimize fuel consumption, while meeting spacecraft pointing requirements. To achieve
the control system specifications for this class of systems, pulse–modulated control
algorithms, which are a special type of discontinuous control laws, are used (Agrawal et
al. [1], Sidi [39], Wie [50]).
Discontinuous control laws also occur in control design for mechanical systems due
to robustness requirements, actuator limitations, and time-optimum specifications. These
requirements lead to the development of control approaches, such as Sliding Mode
method (Edwards and Spurgeon [16], Khalil [31], Slotine and Li [40], Spong [42], Utkin
[48]) and Lyapunov Redesign approach (Corless [13], Corless and Leitmann [14]), which
call for switching between different control levels, leading to discontinuous control
algorithms. Due to switching delays and actuator dynamics, only ultimate boundedness of
the system is guaranteed, and the closed-loop system exhibits periodic behavior. Analysis
of discontinuous control laws of mechanical systems has been mainly performed
assuming a rigid body. However, in many mechanical systems, strict performance
specifications in the form of high bandwidth and the system’s weight limitations cause a
large number of flexible modes to be within the control system’s bandwidth, thus
requiring their inclusion in the system model. The presence of many lightly damped
2
modes in the model may lead to a highly complex periodic behavior that requires careful
analysis. Understanding the relation between the control law, system flexible dynamics
and periodic motion characteristics, such as amplitude and frequency of periodic motion,
are of paramount importance for achieving closed-loop performance specifications.
The class of systems investigated in this dissertation consists of a linear mechanical
system that includes flexible structural modes, a feedback control law that includes a
linear or nonlinear controller, and a discontinuous control law. An important application
of this type of system is controller design for attitude control of a flexible spacecraft
(Agrawal et al. [1], Bryson [12], Junkins [30], Likins [33], Sidi [39], Song and Agrawal
[41], Wie [50]). For example, station-keeping pointing operations require a prescribed
limit cycle with specified frequency and amplitude. When the control torques are
generated using thrusters, strict performance specifications on the limit cycle amplitude
may require including flexible modes in the design model. Flexible modes may alter the
frequency and amplitude of the limit cycle and may generate additional stable periodic
modes of behavior.
In many other applications a limit cycle of a closed-loop system with discontinuous
controller must be avoided. Nonlinear robust control methods such as Sliding Mode and
Lyapunov Redesign lead to discontinuous control laws which may excite flexible modes
and generate a number of stable limit cycles that limit closed-loop accuracy. Finding the
various limit cycles and their domains of attraction are of the utmost importance for
understanding the behavior of the closed-loop system and improve performance.
3
Current solution techniques to solve for the multimodal limit cycle frequencies in
Relay Feedback Control Systems (RFCS) can be grouped into three categories –
Frequency Based, Time Based, and Cell Mapping. Frequency based methods include the
Describing Function Method (DFM) (Gelb and Vander Velde [19], Graham and Duane
[25], Khalil [31], Slotine and Li [40]) and the more accurate Tsypkin’s Method (Tsypkin
[47]). The DFM assumes that the response of the system during a limit cycle is a pure
sinusoid, whereas Tsypkin’s Method allows for more complex limit cycle waveforms,
using a truncated Fourier series approximation. Extensions of the Tsypkin method also
using the Fourier series expansion are done by Atherton (Amrani and Atherton [2],
Atherton [5], Moeini [36]) and Boiko (Boiko [9], [10], [11], Nikolay, Parnferov, and
Boiko [37]). A limitation to the frequency based methods is that they provide no
information on periodic fixed points of the limit cycle. Also, they provide no information
on global stability of the limit cycle, which is especially important when multiple limit
cycles are possible for the system.
Time based methods analyze the relay switching instances using a state-space
formulation (Astrom [4], Goncalves [22], [23], [24]). This allows the determination of a
limit cycle fixed point location in state-space coordinates. Also, the waveform solution is
exact, and is equivalent to the infinite term Fourier series solution used in Tsypkin’s
Method. However, like the frequency methods, these methods provide no information on
global stability for multiple limit cycles.
In order to determine global stability for multiple limit cycles, the cell mapping
method developed by Hsu [28] can be used to map the region of attraction about each
4
periodic fixed point. The cell mapping method discretizes the state space into individual
cells, with each cell center serving as an initial condition for the mapping process. After
the specified integration time, the mapped destination cell of each initial condition is
noted and sorted into different periodic solutions, using a sorting routine developed by
Hsu and Guttalu [27]. Cells with like periodic solutions can then be plotted to show
regions of attraction for each periodic solution. Because accuracy of the map depends on
the cell density, these methods are computationally expensive, and are not attractive for
studying systems with more than four states. An improvement on Hsu’s method was
proposed by Levitas (Levitas et al. [34], Levitas and Weller [35]) where a hyper-surface
is used as a system constraint to reduce the dimension of the cell space by one, compared
to Hsu’s method. This leads to a substantial reduction in the total number of cells,
especially for higher-order systems.
An alternate point mapping method that was investigated finds the region of attraction
of a known periodic solution by backward integrating from a spherical set of seed points
located at each fixed point (Flashner and Hsu [17], Flashner and Guttalu [18], Genesio et
al. [20], Stacey and Stonier [43]). It was found that this method cannot be used with
systems that contain memory, because future knowledge would be necessary for the
correct relay switches to occur during the backward integration process. A feedback
system that contains a relay nonlinearity with memory (hysteresis), therefore is limited to
methods that use forward integration.
A new computational method for accurately computing limit cycles and their regions
of attraction for mechanical systems with multiple flexible modes is presented in this
5
dissertation. The method employs a block-diagonal state-space modal formulation,
allowing for a closed form solution of the relay switching constraint equations used in
determining limit cycle frequencies and fixed points. The method also presents a new cell
mapping technique, which significantly reduces memory requirements with increased
solution accuracy, so that more flexible modes may be considered in the plant model.
In Chapter 2 the problem of relay feedback control of flexible structures is described,
as well as some performance issues that can occur when controller bandwidth approaches
structural resonance frequencies. Also described are the various system configurations
that are considered in this work. In Chapter 3 some theoretical background on current
solution techniques for determining limit cycle frequencies, periodic fixed points, and
their local stability in closed-loop systems is presented. Also presented is the cell
mapping technique, which can be used for determining global stability of limit cycle
fixed points, and their regions of attraction. Chapter 4 presents the formulation of the
Switching Surface Cell Map (SSCM) method, which accurately determines multiple limit
cycle characteristics that can exist in lightly damped flexible structures under relay
feedback control. Four application examples are presented in Chapter 5, which
demonstrate the method on system configurations with different relay types and feedback
control laws. Finally, some concluding remarks and suggestions for future research are
discussed in Chapter 6.
6
Chapter 2: Problem Description
2.1 Relay Feedback Control of Flexible Structures
In many mechanical systems the objective of the control system is to generate predictable
periodic behavior (limit cycle). For example, in the case of spacecraft attitude control
using thrusters, the goal is to generate a periodic motion of given amplitude and
frequency. The amplitude is constrained by positional error requirements, and the
frequency is constrained by minimizing fuel consumption and excitation levels of
structural resonances. The discontinuous nature of the control actuators (on-off) can
excite structural vibrations that can reduce control effectiveness, and may even cause
structural damage. Therefore, it is imperative to evaluate the effects of structural
flexibility on control system performance at the beginning of the design cycle.
The analysis of mechanical systems controlled by discontinuous relay type control
laws has mainly been performed using rigid-body models or models that include at most
one or two flexible modes. These low order models are acceptable for design purposes if
the controlled response is slow enough to avoid exciting the higher frequency vibrational
modes that are ignored in the model. However, performance specifications calling for
lighter weight structures with faster and more accurate positioning requirements results in
the need to consider more flexible modes in the system model. When additional modes
are included in the model, analyzing the effects of flexibility on the control robustness
becomes more complex. The effects are typically analyzed using multiple simulations
within a specified design space. This can become more and more costly, with regards to
7
computational resources and design schedule, as more structural modes are added to the
model.
The research for this dissertation focuses on developing a design methodology and
simulation tool that aid in analyzing the effects of higher-order modes on control system
performance. The ability to analyze the effects of control parameter variation and
structural mode uncertainty needs to be evaluated in a relatively quick timeframe. The
design methodology combines two methods of limit cycle analysis in relay feedback
systems in a new way, and improves on the efficiency and accuracy of the result. The
first part of the method allows for accurate calculation of multiple periodic solutions,
which may occur when higher-order structural modes are considered. The second part of
the method provides a mapping of the region of attraction about the periodic solutions,
which aids in determining robustness of design parameters. The mapping technique has
been improved for quicker solution time, while using less computer memory resources.
What follows is a brief description of the relay feedback system under consideration,
as well as the various relay nonlinearities that are investigated.
2.2 Relay Feedback Control System (RFCS)
The type of feedback control system that this research focuses on is Single-Input-
Single-Output (SISO) position control of a system with a control law that can be
nonlinear in the form of a relay type nonlinearity. In some cases it is assumed that the
velocity output is also available, in order to avoid any signal noise associated with
differentiating the position output. The relay nonlinearity may model the physical on-off
actuator output, or it may be part of a robust control law. Figure 2.1 shows the block
8
diagram of a general relay feedback system, where the position error from the plant
output is used by the controller, to provide input to the relay nonlinearity. In this general
case the plant may be linear or nonlinear, and the control law may also be linear or
nonlinear. However, the actuating force is always a piecewise constant magnitude,
dictated by the input signal.
Figure 2.1: General SISO Feedback Control with Piecewise Constant Nonlinearity
Descriptions for each block in Figure 2.1 are as follows
Plant Dynamics
State Equations
() x fx u B (2.1)
where ( ,) f xt is an n x 1 vector of functions for the n states, and B is an n x 1 vector
of 0’s and 1’s applying the input force to the proper equation.
Output Equations
() y gx (2.2)
where ( ) gx is an m x 1 vector of functions for the m states that are measured or
estimated states, which are feedback to the controller.
() e () N
()
()
x fx u
ygx
B
e
u
y
r
–
Plant Dynamics
Piecewise Constant
Nonlinearity Controller
9
Controller
Tracking problem
error between desired state and measured output
( ) relay input signal
er y
e
(2.3)
where r is an m x 1 vector of functions for the m desired states, and ( ) e is a scalar
function that is a general combination of the error states.
Regulator problem (used for limit cycle analysis)
0 error between 0 and measured output
( ) ( ) relay input signal
ey r
yx
(2.4)
where the desired set point is r = 0, and () x is a scalar function that is a general
combination of the measured or estimated states.
Piecewise Constant Nonlinearity
Figure 2.2: Piecewise Constant Nonlinearity
()0,1,2,...,
i
uN M i n (2.5)
where M
i
is a constant magnitude, and N( ε) is a nonlinear function that determines relay
output u based on relay input ε and dead-zone locations δ
i
.
M
n
δ
1
δ
2
δ
n
M
1
M
2
u
10
In order to investigate the effect of structural modes on limit cycle performance, a
linear plant is used as in Figure 2.3, which includes both a rigid-body mode and
additional flexible modes.
Figure 2.3: General Feedback Control with Relay Nonlinearity and Flexible Modes
When the controller is also linear, e.g. a PID type controller, the control system takes
the form in Figure 2.4. Nonlinear analysis tools using the open-loop frequency response
can then be applied as outlined in Section 3.2.
Figure 2.4: Linear Feedback Control with Relay Nonlinearity and Flexible Modes
Transfer Function Form
The preceding Figure 2.1 shows feedback control of a general nonlinear system,
which contains a relay type nonlinearity. Stable limit cycle frequencies and periodic fixed
points and their region of attraction in such a general system can be found using the new
method presented in this work. However, since a closed form analytical solution is not
generally available, the limit cycle frequencies and fixed points cannot be calculated
() e () N
x xu
yx
AB
C
e
u
y r
–
Linear Flexible
System
Relay
Nonlinearity
General
Controller
u
() Ks () N
()
p
Gs
e y r
–
Linear Flexible
System
Relay
Nonlinearity
Linear
Controller
11
before the region of attraction mapping process is carried out. It is the mapping process
itself that leads to the solution of the stable periodic fixed points and limit cycle
frequencies. This presents some computational disadvantages that will be discussed
further in Section 4.2.1.
Limit cycle characteristics and fixed points for the linear plant versions shown in
Figure 2.3 and Figure 2.4 can be solved for using closed form equations shown in Section
4.2.2 and Appendix A. The closed form equation used depends on the controller and
relay type. The relay may be any of the types shown in Section 2.3.
2.3 Relay Types
Figure 2.5 shows an ideal relay (bang-bang) nonlinearity. This is the simplest of the
relay type nonlinearities, and when combined with a PD controller provides the most
accurate position control when employed in a feedback configuration. However this
nonlinearity has the shortcoming of chattering (high frequency cycling), when the relay
input (usually position error) is near zero. In the case of spacecraft attitude control using
thrusters, these high frequency on-off commands to each thruster results in excessive fuel
consumption.
Figure 2.5: Ideal Relay (D=1)
0 0.5 1 1.5 2
2
1
0
1
2
xt ()
yt ()
t
12
0 0.5 1 1.5 2
2
1
0
1
2
xt ()
yt ()
t
Figure 2.6 shows a relay with hysteresis nonlinearity. This relay is not single valued
within the hysteresis band, and has output depending on the immediately preceding
output, when the magnitude of the input signal is in this area. This is commonly
described as a relay with memory. Note the presence of a phase lag between the input and
output signals, when compared to the above ideal relay. Since there is no dead-zone at
zero input, this relay may cause a relatively low frequency limit cycle at zero input, as
compared to the ideal relay. In the case of spacecraft attitude control with thrusters, this
constant command results in the thruster valve being fully open, even when there is a
zero error input, resulting in wasted fuel.
Figure 2.6: Relay with Hysteresis (D=1, δ=.5)
The next relay shown in Figure 2.7 is classified as on-off, since there is a dead-zone
providing zero command signal at zero input signal. This can be especially desirable in
the case of regulator problems, where the system spends most of its time in a quasi
steady-state condition around the desired set-point, resulting in a relatively constant near
zero input to the relay. The drawback to the dead-zone is that it results in a steady-state
error, since no correcting command is available within the dead-zone. However, for the
spacecraft control problem, this trade-off on positional accuracy is countered by the
13
advantage of fuel efficiency, since the thruster valves remain closed within the dead-
zone.
Figure 2.7: Relay with Dead-Zone (D=1, δ=.5)
Another on-off type relay is shown in Figure 2.8. This relay with dead-zone and
hysteresis, also known as a Schmidt Trigger, has the same advantages and disadvantages
as the above relay with dead-zone. Note that a phase lag between input and output occurs
due to the hysteresis band, when compared to relay with dead-zone. This hysteresis band
results in the same memory characteristic as discussed in the above relay with hysteresis.
This type of nonlinearity is useful for simulating electro-mechanical actuators that
require a higher current to actuate (generate electromagnetic force across an air gap to
engage), than to de-actuate (reduce electromagnetic force below restoring spring load).
However, for the spacecraft control problem, the Schmidt Trigger is commonly used in a
sublevel feedback loop in the controller to produce pulse-width pulse-frequency (PWPF)
modulated signals that power on-off actuators, such as thruster valves (see Figure 2.9).
This control scheme, results in a quasi-linear average output torque applied to the
spacecraft, resulting in better fuel economy.
0 0.5 1 1.5 2
2
1
0
1
2
xt ()
yt ()
t
14
Figure 2.8: Relay with Dead-Zone and Hysteresis (D=1, ε=2, δ=.25)
Figure 2.9: Pulse-Width Pulse-Frequency (PWPF) Modulator
A more general version of the relay with dead-zone and hysteresis is shown in Figure
2.10. In this case there are multiple levels, each with a possibly different dead-zone and
hysteresis.
Figure 2.10: Multilevel Relay with Dead-Zone and Hysteresis
M
n
δ
1
δ
2
δ
n
M
1
M
2
u
Δ
n
Δ
2
Δ
1
0 0.5 1 1.5 2
2
1
0
1
2
xt ()
yt ()
t
15
Also considered in this work is a relay with linear-zone or saturation nonlinearity as
shown in Figure 2.11 (see Appendix B).
Figure 2.11: Relay with Linear-Zone (Saturation Nonlinearity)
In the sections that follow, current methods for analyzing the existence of limit cycles
and their characteristics in linear feedback systems will be outlined. The systems
investigated will be limited to single-input-single-output (SISO) systems with one relay
nonlinearity. In the state-space formulation it is assumed that all of the states are
available for feedback. Multiple-input-multiple-output (MIMO) systems and systems
containing more than one relay are beyond the scope of this research.
M
δ
u
M
k
16
Chapter 3: Background
3.1 Limit Cycle Analysis in Relay Feedback Control Systems
Here we describe existing methods used for limit cycle analysis for linear systems with
linear control laws. The methods are presented in such a way as to highlight their relative
accuracy in calculating limit cycle characteristics. The Time Domain Method presented
in Section 3.3 is typically used to investigate limit cycles in linear process plants with
linear feedback and a single limit cycle (Astrom [4] and Goncalves [22]). A new
contribution of this work is that it extends this method to structural vibration problems
that contain multiple limit cycles (Section 4.2.2) and nonlinear feedback (Sections 4.2.5
and 5.3.4).
Figure 3.1 shows a regulator configuration (r = 0) of a linear flexible system with a
linear controller. Figure 3.2 shows this same system rearranged in order to separate the
relay nonlinearity from the linear portions of the closed loop. Separating the nonlinearity
from the linear portions in the control loop is the first step in the limit cycle analysis
procedure.
Figure 3.1: Relay Feedback Control of Linear Flexible System with r = 0
0 r
() Ke () N
x xu
yx
AB
C
e
u
y
–
Linear Flexible
System
Relay
Nonlinearity
Linear
Controller
17
Figure 3.2: Equivalent System Used for Limit Cycle Analysis
3.2 Frequency Domain Techniques
For the frequency domain techniques it is convenient to put the system equations into
transfer function modal summation form. The modal analysis method is valid for lightly
damped structures, where no physical damping components (e.g. dashpots) are present.
This form can be developed using the undamped n x n system of equations in
physical coordinates
qt qt t MK F (3.1)
where q is an n x 1 vector of physical coordinates, M and K are symmetric n x n mass
and stiffness matrices, and F is and n x 1 force vector. Since the systems investigated are
lightly damped structures, damping is accounted for by adding modal damping to the
transformed set of equations, resulting in the decoupled modal equations with damping
2
(2 ) ( )
T
jj
xtdiag xtdiag xt t UF (3.2)
where the variable substitution q = Ux has been made, and U is the n x n modal matrix
normalized with respect to the mass matrix M (U
T
MU = I). The variables and
j
are
the jth structural natural frequency and assumed modal damping ratio (typically .001-.05
depending on application). The above equation (3.2) is the form that Finite Element
1
() Ne
x xu
y x
AB
KC
e u
y
Linear Controller and
Flexible System
Relay
Nonlinearity
18
Analysis results are reported in, and these results (natural frequency, modal damping and
mode shape) can be applied to the following formulation by substituting the appropriate
values.
Since the equations in (3.2) are uncoupled, they can be solved for individually using
the Laplace Transform method as follows
22
1
1
2 ( is the th row of
or the th column of )
TT T
jjj jj j j
jnj
n
sx s sx s x s s j
j
F
uu
F
UF U U
U
222 2
1
1
()
22
T
n
j
jkjk
k jj jj
s
xsuFs
ss ss
UF
(3.3)
The form of (3.3) assumes zero initial conditions, since we are ultimately interested in the
frequency response, which is a steady-state response of the system. Since the input/output
response is needed in physical coordinates for output feedback control purposes, the
uncoupled solutions in modal coordinates (3.3) need to be back-transformed into the
original physical coordinate system using the modal matrix as follows
111 1 1
1
11
1
() ()
() () () ()
n
nn nn n
n
ii inn ijj
j
qs x s
qu u x
qu u x
qs u x s u x s ux s
U
(3.4)
Substituting (3.3) into (3.4) yields the response in physical coordinates using the modal
form
19
22
11
22
11
1
() ( )
2
() ()
2
nn
iij kjk
jk jj
nn
ij kj
ik
jk jj
qs u u F s
ss
uu
qs F s
ss
(3.5)
Equation (3.5) represents the position response measured at physical coordinate location i
due to force inputs at multiple physical coordinate locations k = 1,2,…,n. For example, if
the flexible structure is simplified as a lumped parameter mass-spring model, equation
(3.5) represents a transfer function for the output measured at mass i due to an input force
on mass k. For the collocated SISO systems under investigation, the single force input
(index k in (3.5)) is at the same physical coordinate as the single measured output (index i
in (3.5)). This results in the following summation form transfer function for a plant with n
flexible modes and collocated input/output at coordinate i.
22
1
2
22
1
( ) ( ) single input at single output at
2
( ) collocated input/output at
2
n
ij kj
ik
j jj
n
ij
i
j jj
uu
qs F s k i
ss
u
Fs k i
ss
2
22
1
()
()
() 2
n
ij
i
i
jijj
u
qs
Gs
Fs s s
(3.6)
where
ij
u is the component of the modal matrix U, located in the ith row and jth column.
The magnitude of this component represents the contribution of the jth structural mode to
the overall response seen at coordinate location i.
The Open-Loop transfer function at location i of the linear part of the system shown
in Figure 2.4 including a PID controller is
20
22
1
() () ()
2
n
j
I
ppd
j j j
a
K
Gs K s G s K K s
ss s
(3.7)
where
2
j ij
au . The parameter ζ is the modal damping ratio, and ω
j
is the undamped
natural frequency of the jth structural mode. The parameters K
p
, K
I
, and K
d
represent the
proportional, integral, and derivative gains of a PID controller for the closed-loop system.
3.2.1 Describing Function Method
The following two frequency domain techniques make use of the frequency response
characteristics of the linear system, while undergoing self-sustained steady-state
oscillations or limit cycles. When the linear system’s response during a limit cycle is
matched with the relay switching frequency, solutions provide possible stable limit cycle
trajectories of the system, subject to stability and waveform criteria.
The Describing Function Method (DFM) for limit cycle analysis can be applied to
systems that contain a single dominant nonlinearity that can be separated from the linear
parts of the system by rearranging the block diagram in Figure 2.4 to that shown in
Figure 3.3.
Figure 3.3: Equivalent System in Describing Function Form
1
() Ne
() ()
p
KsG s
e u
y
Linear Controller and
Flexible System
Relay
Describing Function
21
Under the self-sustained oscillations of a limit cycle, the system regulation about zero
input is analyzed using the following characteristic equation from the block diagram in
Figure 3.3.
1()()()0
p
Ne KsG s (3.8)
Equation (3.8) is called the “harmonic balance equation” (Khalil [31]), and represents the
condition for a limit cycle to exist. This equation can be solved graphically (most
common) or analytically.
The equation is solved graphically using a modified Nyquist plot with the frequency
response of the linear part, plotted against the negative inverse of the Describing Function
(DF). Intersections of these two curves on the plot indicate a solution to the following
alternate form of the harmonic balance equation, and possible solutions for limit cycle
frequency and amplitude.
1
() ()
p
Ks G s
Ne
(3.9)
The DFM assumes a steady-state sinusoidal waveform ( ) sin( ) et A t input to the
relay. This input causes a periodic switching of control magnitude to occur at the
hysteresis settings above and below the zero crossing of the error signal as shown in
Figure 3.4. The relay hysteresis causes a phase lag of the output with respect to the input
of
N
.
22
Figure 3.4: Input/Output for Relay with Hysteresis Under Limit Cycle
It is also assumed that the system plant will act as a low-pass filter, so that the system
will only be excited by the 1
st
harmonic of the square-wave relay control signal (see
Figure 3.5), and that this filtered response will be passed through to the controller via
feedback.
Figure 3.5: Single Term, 2 Term, and 10 Term Fourier Series Representation of Square-
Wave Relay Output Under Limit Cycle
This 1
st
harmonic term is used as a quasi-linear approximation of the nonlinear relay
in limit cycle analysis. With this approximation, the system’s open-loop frequency
response can now be used to determine the existence and characteristics of any limit
cycles, using standard linear system tools such as the Nyquist plot.
The Describing Function (, ) NA acts as an equivalent gain in the control loop
shown in Figure 3.3. Note that in general this equivalent gain can be a function of both
Relay Input
() sin( ) et A t
Relay Output = +/-M
Phase Lag
N
Hysteresis
Only frequency that
excites system
23
limit cycle amplitude and frequency, however for the relay type nonlinearities
investigated, the DF will be a function of amplitude only. Evaluation of the first Fourier
series term of the relay square-wave output yields the DF gain and phase shift that occurs
when a sinusoidal signal is passed through it. Details on how to calculate Describing
Functions for various nonlinearities can be found in Atherton [5], Gelb and Vander Velde
[19], Graham and Duane [25], Khalil [31], and Slotine and Li [40].
For a relay with hysteresis, the Describing Function is
2
4
() 1
M
NA j
AA A
(3.10)
2
1
() 1
() 4 4
A
Aj
NA M A M
(3.11)
where is the hysteresis band, M is the relay magnitude, and A is the amplitude of the
sinusoidal output of the first-order approximation of the relay. The amplitude A is varied
for plotting purposes.
Figure 3.6 shows a sample modified Nyquist plot showing the frequency response of
a single lightly damped flexible mode .01, 10
n
with proportional control (K
p
=
1), plotted against the negative inverse plot of the DF for a relay with hysteresis (3.11)
.1, 1 M .
24
Figure 3.6: Modified Nyquist Plot of Linear System Frequency Response (red) and
Negative Inverse of Describing Function (blue) for Relay with Hysteresis
The above intersection corresponds to a limit cycle of 10.23rad/s (period=0.61s) with
amplitude of .25. The following Simulink simulation using a system block diagram per
Figure 3.3 with an initial output of .25 verifies the results.
Figure 3.7: Verification of Limit Cycle Predicted by DFM
Limit cycle frequency candidates are tested for stability using Loeb’s criteria (3.14) as
described in Atherton [5].
Let the terms in (3.9) be represented by
() () () ()
p
Kj G j U jV (3.12)
0.3 0.2 0.1 0 0.1 0.2 0.3
0.5
0.4
0.3
0.2
0.1
Im G () ()
Im a () ()
Re G () ()Re a () ()
Intersection is a
solution to equation
(3.9), yielding possible
limit cycle frequency
and amplitude.
25
1
() ()
()
PA jQ A
NA
(3.13)
The condition for stability of each limit cycle ( ω
c
, A) satisfying (3.9) is (Atherton [5])
0
c
UQ V P
AA
(3.14)
3.2.2 Tsypkin’s Method
Tsypkin’s method (Tsypkin [47]) is described as an exact method of determining limit
cycle conditions, as opposed to the DFM, which is an approximate method. The
approximation of the DFM is due to the low-pass filtering assumption placed on the
linear plant. Tsypkin’s method matches the relay output per limit cycle period, with the
linear system response. For the regulator problem being considered, the negative open-
loop linear response serves as the input to the relay, causing sustainable self-oscillations
of the system, through negative feedback (see Figure 3.3). Constraints are placed on the
linear system’s response at the relay half-period switching times (refer to Figure 3.8).
These requirements on the relay input cause the relay to switch at the correct limit cycle
period. The advantage over the DFM is that the waveform of the relay input can now be
more complex than a simple sinusoid, thus relaxing the low-pass filtering assumption for
the linear plant.
26
Figure 3.8: Sinusoidal Relay Input and Square-Wave Output Under Limit Cycle
The following Tsypkin frequency response locus J( ω) is used to estimate limit cycle
frequencies in a manner similar to the G( ω) locus for the DFM. Note that the following
summation form is composed of variations on the U( ω) and V( ω) terms used for the DFM
in (3.12). The variation consists of evaluating U(m ω) and V(m ω) for an odd frequency
multiple m ω, at each successive summation term.
1,3,5...
4()
( ) ( ) , where is the relay magnitude
m
MVm
JUmj M
m
(3.15)
Note that the additional m terms in the summation are Fourier Series terms of the true
response of the system to a square-wave input. Including more terms will provide a more
accurate representation of the actual frequency response. In using this method, the
imaginary part of the frequency response ( ) J is equated to the hysteresis setting of the
relay.
1,3,5...
4()
Im ( ) = , where is the relay magnitude
and is the relay hysteresis
m
MVm
JM
m
(3.16)
Relay Input is
negative of linear
system response
Relay Output = +/-M
Phase Lag
N
Hysteresis
Constraint set1
()
() 0
o
o
et
d
et
dt
Constraint set2
(2)
(2) 0
eT
d
eT
dt
ey
27
The frequency that satisfies this constraint is a potential limit cycle frequency. This is
similar to the DFM, where the system frequency response is equated to a first order
approximation of the relay square-wave output (i.e. the relay’s Describing Function). The
improvement over the DFM solution comes from the extra terms now accounted for in
the Fourier Series of the relay response in (3.16). Typically only the first 3 to 5 terms are
necessary to get a reasonable answer on the limit cycle frequency for the problems
investigated.
For purposes of comparing the two methods, note that the first summation term in
(3.15) is related to (3.12) in the following way
44
() () () ( ), for 1
MM
JU jV Gj m
(3.17)
The following Figure shows a scaled version of the Tsypkin locus ()
4
J
M
frequency response of (3.17) to (3.12). Note that the curves lie on top of one another.
Figure 3.9: Linear System Frequency Response (red) and Single Term Tsypkin Locus
(green) and Negative Inverse of Describing Function (blue) for Relay with Hysteresis
0.3 0.2 0.1 0 0.1 0.2 0.3
0.5
0.4
0.3
0.2
0.1
Im G () ()
4M
Im J () ()
Im a () ()
Re G () ()
4M
Re J () () Re a () ()
Intersection is a
solution to equation
(3.9), yielding possible
limit cycle frequency
and amplitude.
28
Since the negative inverse of the DF has an imaginary component that is constant
(refer to (3.11)), an alternate representation plotting the Imaginary part of the frequency
response versus frequency is more convenient for finding frequencies that yield solutions
(see Figure 3.10).
Figure 3.10: Imaginary part of Frequency Response Im(G) (red) and Single Term
Tsypkin Locus Im(J) (green) and -1/N(a) (blue) for Relay with Hysteresis vs Frequency
The same criterion for the DFM can be applied to Tsypkin’s method by rearranging
(3.16) as follows
1,3,5...
1,3,5...
4()
Im ( ) =
()
=
4
m
m
MVm
J
m
Vm
mM
(3.18)
0 5 10 15
0.4
0.2
0
0.2
0.4
Re G () ()
4M
Re J () ()
0 5 10 15
0.1
0.08
0.06
0.04
0.02
0
Im G () ()
4M
Im J () ()
Im a () ()
a
Intersection is a
solution to equation
(3.9), yielding
possible limit cycle
frequency and
amplitude.
The other intersection
is not a solution to
equation (3.9), since it
corresponds to
Re[ ( )] G in the RHP.
29
and noting that the negative inverse of the imaginary part of the Describing Function for
a relay with hysteresis (3.10) is
1
Im[ ]
() 4 NA M
(3.19)
Note that higher-order harmonics are introduced by the Im[J( ω)] locus in (3.18), due
to each additional summation term. Compared to the DFM, which only uses the 1
st
harmonic, this will introduce additional peaks with potential for intersection with the
relay hysteresis line at – Δ (see Figure 3.11).
Figure 3.11: Imaginary part of Frequency Response Im(G) (red) and Two Term Tsypkin
Locus Im(J) (green) and -1/N(a) (blue) for Relay with Hysteresis vs Frequency
These additional intersections present additional potential limit cycle frequencies,
which must be tested with an additional constraint. The waveform constraint ensures that
these new solutions will not violate the initial assumption, and cause more than two
switches per period of the relay (see Figure 3.8). Such solutions are invalid, and will not
cause a limit cycle. Note that this waveform criterion is automatically satisfied when
using the DFM, since it is initially assumed that the limit cycle has the shape of a simple
sinusoid.
0 5 10 15
0.1
0.08
0.06
0.04
0.02
0
Im G () ()
4M
Im J () ()
Im a () ()
a
Additional
possible solution
caused by
additional
summation term in
Tyspkin locus
30
3.3 Time Domain Methods
Frequency domain methods such as the Describing Function Method can provide
estimates of limit cycle frequency and amplitude, but not periodic fixed point locations
(or initial conditions). The limit cycle fixed point locations are necessary for defining the
Poincaré mapping region of interest. In order to find the fixed points, it is necessary to
cast the problem in the following state-space form.
Figure 3.12: Equivalent State-Space Time Domain Model
if or and ( )
()
if or and ( )
output matrix includes linear controller gains
x xu
yx
M yt y ut M
ut
Myt y ut M
AB
C
C
(3.20)
The system input u(t) represents a relay nonlinearity with hysteresis , and
magnitude M, based on the filtered position feedback (see Figure 3.13).
In order to determine the limit cycle frequency using the state-space equations, it is
necessary to calculate the time at which the relay switches under limit cycling conditions.
The formulation of the current problem requires that the limit cycle be symmetric, with
1
() Ne
x xu
yx
AB
C
e u
y
Linear Controller and
Flexible System
Relay
Nonlinearity
31
only two switches per cycle. These are the characteristics of a unimodal periodic motion,
with the waveform and input/output relations shown in Figure 3.13.
Figure 3.13: Symmetric Unimodal Limit Cycle Waveform in RFS
As shown in Astrom [4] and Tsypkin [47], in order for the symmetric unimodal limit
cycle in Figure 3.13 to exist in system (3.20), the following must hold
( ) ( ) the motion is periodic with period
( ) ( 2) the motion has odd symmetry
ytyt T T
yt yt T
(3.21)
This provides constraints on the system motion for establishing the periodic conditions
required for the limit cycle to exist. Using the relations in (3.21), the initial conditions of
the state must satisfy
** *
*
(0) (0)
() ()where 2
() (0)
y x
ytxt t T
xt x
C
C (3.22)
Also, the system is piecewise linear between relay switches, so the solution to (3.20)
during these time periods is
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
1
0
1
1
0
1
Limit Cycle Waveform (Relay Input u(t), Linear Plant Output y(t))
t
yt ()
ut ()
t ()
t ()
t
t = t* = T/2
t = 2t* = T
t = 0
32
() ()
() ( ) ( )
o
o
t
tt t
o
t
xte xt e u d
A A
B (3.23)
For a relay with hysteresis, the input is constant during this time (+/-M), resulting in
the following state dynamics
*
*
for 0
for
x xM t t
x xM t t T
AB
AB
(3.24)
At relay switching time t = t* the output y(t) is just crossing the lower hysteresis
boundary of the relay (see Figure 3.13). The state solution (3.23) at this instant in time
and (3.22) can be used to determine the initial conditions that will start the limit cycle
*
*
*
*
*
0
*
0
() e (0) ( )
() (0) , e ,
(0) (0)
t
t
t
t
xt x e d M
xtx M ed
xx M
AA
AA
B
ΦΓ Φ Γ B
ΦΓ
1
*
(0) x IM x
ΦΓ (3.25)
This yields the following constraint equation that must be satisfied during the limit cycle.
1
(0) (0) y x
M
C
C Φ I Γ
*
*
1
0
(0) e
t
t
yedM
AA
CI B (3.26)
Solutions of t* which satisfy (3.26) must be obtained numerically, and determine possible
limit cycle frequencies
**
t . There are typically multiple solutions, but not all yield
a limit cycle. Since the state matrices in (3.20) can be expressed in block-diagonal modal
33
form, the matrices in (3.26) can be expressed in a closed form (see Section 4.2.2 and
Appendix A for details). This negates the need for numerical techniques when inverting
or integrating the matrix terms in (3.26), thus increasing the accuracy of solutions for t*,
even for large matrices containing many flexible modes.
The following additional waveform criteria must be satisfied for valid solutions (see
Figure 3.13).
*
*
) (0) 0, ) ( ) 0
) ( ) for 0 (no more than two relay
switches per period)
iy ii yt
iii y t t t
(3.27)
In order to compare the solutions of (3.26) with those obtained using frequency
domain methods, substitute
*
t into (3.26) to yield the “State-Space Locus”
1
0
() efedM
A A
CI B (3.28)
and plot the magnitude versus frequency as shown in Figure 3.14.
Figure 3.14: Frequency response using State-Space Locus f( ω) (red) and Two Term
Tsypkin Locus Im(J) (blue) and – Δ (green) for Relay with Hysteresis vs Frequency
0 5 10 15 20
0.15
0.125
0.1
0.075
0.05
0.025
0
f ()
Im J () ()
Additional possible
solutions of (3.28)
correspond to each
additional summation
term in Tyspkin locus
(3.18)
34
Note in Figure 3.14 that the State-Space Locus created by equation (3.28) lies on top
of the imaginary part of the Tsypkin Locus (3.18), for the number of Tsypkin Locus
summation terms. Thus the Tsypkin Locus approaches the State-Space Locus as more
and more terms are included in the Tsypkin Locus summation. Therefore, the most
accurate locus for finding limit cycle frequencies is the State-Space Locus (3.28). In
order to determine valid solutions that will yield possible limit cycles, frequencies
satisfying equation (3.28) are tested using the waveform criteria stated above in (3.27)
and noting that
**
t . Initial conditions x(0) = x* for these valid solutions are fixed
points of the periodic motion such that
**
**
(0)()(2)
() (0)
x xT x t x
xt x x
(3.29)
Constraint equations that yield initial conditions for other relay types can be found in the
application problems in Chapter Chapter 5: and Appendix B.
For fixed point stability, eigenvalues of W must lie within the unit circle (Astrom [4])
v
v
C
WI Φ
C
(3.30)
*
* *
*
*
where initial condition fixed point on switching surface ( )
, limit cycle half-period
state velocity just before switch from to M at (0 )
t
xyx
et
vx M M x
A
Φ
AB
Stability for a relay which includes a dead-zone and more general switching surface is
developed in Section 4.2.5. Also, Section 5.3.4 contains an example of how to determine
stability on a nonlinear switching surface.
35
3.4 Methods to Calculate Region of Attraction
The eigenvalue test on the linearized Poincaré map W presented in Section 3.3 can only
establish local stability of a limit cycle periodic fixed point. This test does not provide
information about how large the perturbation of a fixed point can be, before the state
trajectory is attracted to a different fixed point location, leading to different limit cycle
characteristics. In order to determine the effects of adjacent fixed points, the global
stability of each fixed point must be determined. Point mapping methods provide a tool
for mapping the region of attraction about each fixed point, and can thus be used to
establish the global stability characteristics of the individual limit cycle fixed points.
3.4.1 Simple Cell Mapping
There are two main types of Cell Mapping Methods that are typically used for region of
attraction analysis. The first method, known as Simple Cell Mapping (SCM), was
developed by Hsu [28] and is essentially a time based Poincaré Map. It is time based in
that each point mapping consists of integrating the system equations of motion over a
predetermined fixed integration time period T.
The region of interest in an n dimensional state space is discretized into N cells of
uniform size h
i
along each dimension x
i
, i = 1,2,…,n, defining the cell center location in
state space
() d
j
x for a total of
n
j N cells. Each of these points
() d
j
x serves as the cell
center for cell Z
j
and initial condition during the cell mapping process. Figure 3.15 shows
the discretization of a 2D state space (n = 2) divided into 3 cells (N = 3) along each
36
dimension x
i
, i = 1,2. Therefore, the total number of cells , 1,2,...,
n
j
j N Z , in the cell
space is
2
39.
n
N
Figure 3.15: 2D Cell Map Centers in State-Space Coordinates
By definition the components of each cell center Z
j
in cell-space coordinates are
integers, and are related to the cell center state-space coordinates by
()
1
, 1,2, ,
2
d
i
i
i
x
zInt i n
h
(3.31)
where Int(y), for any real number y, represents the largest integer which is less than or
equal to y, i.e. ( ) Int y y . Figure 3.16 shows the cell-space coordinates of the cell
centers depicted in Figure 3.15.
x
1
(d)
(0,0)
Z
4
Z
7
Z
2
(h
1
,0)
Z
5
Z
8
Z
3
(2h
1
,0)
Z
6
(2h
1
,h
2
)
Z
9
(2h
1
,2h
2
)
x
1
x
2
(h
1
,2h
2
)
(h
1
,h
2
) (0,h
2
)
(0,2h
2
)
h
1
h
2
x
2
(d)
Z
1
x
4
(d)
x
7
(d)
x
3
(d)
x
5
(d)
x
6
(d)
x
8
(d)
x
9
(d)
37
Figure 3.16: 2D Cell Map Centers in Cell-Space Integer Coordinates
The center point of the destination cell that the system trajectory ends up in after one
integration period T is the mapping of the originating cell.
()(())
jp j j p
tT t ZCZ (3.32)
where t
p
is the integration start time for the current mapping. The cell-space coordinates
of the destination cell center ( )
jp
tT Z are determined from the state-space coordinates
of the destination state
()
()
d
jp
tT x by
()
()
1
() ( ( )) , 1,2, ,
2
d
ip
ip i i p
i
xt T
zt T C zt Int i n
h
(3.33)
where
()
()
d
ip
x tT is the ith component of the destination state started from originating
cell center
() d
j
x . For example, referring to Figure 3.17, Cell Z
1
maps to Z
15
( ((1)) (15) zz C ), Z
15
maps to Z
8
( ((15)) (8) z z C ), and Z
8
maps back to Z
15
( ((8)) (15) zz C ). For a periodic solution to exist, trajectories from the originating cell
eventually return to this cell after k integration periods of duration T. This is called a
“P-K” periodic solution. Figure 3.17 shows a P-2 periodic solution, initiating from Cell
Z
1
(0,0)
Z
4
(0,1)
Z
7
(0,2)
Z
2
(1,0)
Z
5
(1,1)
Z
8
(1,2)
Z
3
(2,0)
Z
6
(2,1)
Z
9
(2,2)
z
1
z
2
38
Z
15
and returning to Cell Z
15
, after 2 integration periods. Since the initial trajectory from
Cell Z
1
landed in Cell Z
15
, Cell Z
1
also eventually has this same P-2 periodic solution. In
fact any trajectory that lands in Cell Z
1
or Cell Z
8
or Cell Z
15
has this same P-2 solution,
since it has already been determined that Cell Z
15
maps back to itself after 2 integration
periods. Cells that have trajectories that land outside of the cell-space boundary are
mapped to the Sink Cell. Note that the period of integration T is typically not the same as
the limit cycle period. An efficient sorting algorithm developed by Hsu and Guttalu [26]
sorts and characterizes each individual cell with what Hsu defines in (Hsu [28]) as a P
Number and Group Number, which differentiates the state space into regions of attraction
for each different P-K solution that exists.
Figure 3.17: System Trajectory For Integration Period “T” in State-Space Coordinates
The cell mapping process is essentially defining the mapping of every initial state that
lies within the originating cell boundary to every final state that lies within the boundary
of the destination cell. Thus there is a discretization error that occurs when the cell-space
Z
16
Z
8
Z
12
Z
4
Z
13 Z
14
Z
15
x(0)=x
1
(d)
Z
5
Z
9
Z
2
Z
6
Z
10
Z
3
Z
7
Z
11
x
1
x
2
Z
1
Discretization
Error
2
() T
Discretization
Error
1
() T
x(T)
39
discretization is relatively coarse (see Figure 3.17). To improve the mapping results, the
discretization is refined. As noted above, the cell space consists of
n
N cells, so the cell
space and thus computer memory requirements grow exponentially with refinement.
While this method is attractive in its simplicity, the computer memory requirements and
resulting solution times makes this method impractical for systems with more than 2
degrees-of-freedom (4 states).
3.4.2 Poincaré-Like Simple Cell Mapping
The second Cell Mapping Method by Levitas [34] is a spatially based Poincaré Map.
It is spatially based in that each mapping consists of integrating the system equations of
motion from an originating cell, until a predetermined spatial barrier (defined using state
variables) is crossed, independent of the integration time to reach this barrier. The center
point of the destination cell that the system trajectory ends up in is the mapping of the
originating cell. Note that the number of intersections with the barrier can be assigned as
a parameter for stopping the integration routine, and defining the destination cell (map of
the originating cell). The longer that the system equations are integrated, the closer the
system trajectory will come to its true steady-state solution, and the more accurate will be
the map from the originating cell to the destination cell. Thus, along with cell
discretization size, the number of barrier intersections can be set to obtain a more
accurate solution. The solution for a fixed number of cells can be improved by increasing
the number intersections. This allows for less memory storage requirements than the
above Hsu method. The trade-off is that more intersections with the barrier requires a
longer integration time, which may lead to a longer solution time. The main improvement
40
to Hsu’s method is that the spatial barrier can be defined as a hyper-plane of dimension
1 n , so that the resulting cell space is reduced to
(1) n
N
cells. This makes the method
more attractive for higher-order systems. Levitas’ method takes full advantage of the
sorting algorithm mentioned above, by defining the cell space in a similar manner as Hsu.
Figure 3.18 shows an example of the Poincaré surface for a 2D system (n = 2) that
would be discretized into a 1D cell space (n – 1). Trajectories are integrated forward from
each cell center for a prescribed number of surface intersections (Int
max
), defining the
final location or “map” of each cell (Levitas [34]). In Figure 3.18, the trajectory started at
the center of cell Z
6
is integrated until two intersections with the Poincaré surface are
encountered (Int
max
=2). The cell at this location (Z
5
) is the map of the originating cell.
The center of this cell then becomes the starting point for integration until the specified
number of intersections is reached. This is repeated until a previously processed cell is
reached (Z
5
in this case).
41
Figure 3.18: 2-State System Trajectory Between Poincaré Surface Intersections
Hsu’s method defines the system integration stopping criterion for each mapping as
the integration period T. This predetermined period is not totally arbitrary, and depends
on the system time scale and state-space discretization. The integration time must be long
enough for the system trajectory to exit the originating Cell. The accuracy of the resulting
map which defines the region of attraction is improved by refining the discretization of
the state space.
Levitas’ method defines the system integration stopping criterion as the number of
intersections with the predetermined hypersurface. The number of intersections is not
totally arbitrary, and like the Hsu method depends on the system time scale and state-
space discretization. If the discretization is course and the number of intersections is too
few, originating Cells may incorrectly map to themselves, since the trajectory has not had
42
enough time to evolve (map) to an adjacent cell. The number of intersections needed can
be determined by increasing the number of intersections, until the overall map does not
change, in other words converges to the true cell map with acceptable error.
43
Chapter 4: Switching Surface Cell Mapping Method (SSCM)
4.1 Introduction
The proposed new method called SSCM combines the time domain and cell mapping
methods in a novel way, providing increased accuracy, while using less computer
memory. This allows for the detailed study of limit cycle characteristics and their regions
of attraction for flexible systems with many structural modes. The control law for these
systems can be either linear or nonlinear, and applications of the method for both types
are presented in Chapter Chapter 5:. First, the time based state-space method for limit
cycle determination described in Section 3.3, which is typically applied to plants with
real poles and a single limit cycle (see Astrom [4], Goncalves [22]), is now applied to
structural vibration problems, where the plant contains a rigid-body mode and many
structural modes. These structural systems may exhibit multiple stable limit cycles. Once
the limit cycle frequencies and periodic fixed points are known, the region of interest can
be mapped to determine the regions of attraction and global stability for each limit cycle.
Second, the new method includes an improved version of the spatially based
Poincaré-Like Simple Cell Mapping Method (Levitas [34]), in order to significantly
reduce the memory requirements described above for Simple Cell Mapping (Hsu [28]).
Since the periodic fixed points lie on the relay switching surface, the switching surface is
selected as the Poincaré mapping surface. The improvement is in the integration stopping
criterion. The required number of crossings with the switching surface is automatically
determined by measuring the distance from the previous crossing point on the switching
surface. When this distance becomes less than the predetermined error criterion, the
44
mapping for that point is considered complete, and the algorithm moves on to the next
point in the cell space. Note that a true periodic solution will cross the switching surface
from the same direction at the same location (fixed point). Using this method allows for
the minimum number of switching surface crossings required to obtain the specified
accuracy. It does so automatically, so specified convergence is reached without the need
for multiple trial-and-error runs, as with Levitas’ method. This also means that a
specified accuracy can be achieved with a relatively coarse cell space, so that memory
requirements are much less than with Hsu’s method.
4.2 Formulation of the SSCM
4.2.1 General Plant Dynamics
Figure 4.1 shows the block diagram for a nonlinear plant controlled by a general
control law in the rearranged regulator form (r = 0), where the input vector B applies the
control force u(y) to the appropriate state equation, and ( ) y x is the controller output
that determines the relay switching. For the nonlinear plant, a closed form of the state
solution is generally not available, so it is necessary to solve the system of state equations
numerically. This means that the region of interest where the limit cycle periodic fixed
points lie cannot be obtained using methods from Section 3.3. However, the new cell
mapping process can still be implemented on such a system, and will yield the fixed point
locations, as well as the limit cycle periods that lie within the map boundaries. Therefore,
the only limitation caused by the nonlinear plant is that the boundaries of the cell map
must be determined by other methods. If the nonlinearities are weak, then a linearized
45
version of the plant may provide an acceptable first order estimate of the region of
interest. This is demonstrated in an application example in Section 5.4, where a nonlinear
softening spring is modeled.
An alternative for estimating cell map boundaries is to perform numerical simulations
from initial states at the outer boundaries of expected performance, where intersections
with the switching surface will provide outer boundaries for the cell map. Note that initial
conditions on the map that are attracted to fixed point locations that lie outside of the cell
mapping boundaries will show up as cells that map to the Sink Cell. Subsequent
simulations using these cells as initial conditions will provide locations on the switching
surface for modifying the cell map boundary.
Figure 4.1: SISO Regulator Control (r = 0) of Nonlinear Plant with General Feedback
Control Signal to Relay Nonlinearity
Since the trajectory for a nonlinear system can only be evaluated using numerical
integration, the total solution time is slower than for the linear plant case discussed in the
next section. This is because the block-diagonal modal form used in the linear case
allows for a closed form of the solution, and therefore the exact value of the state
trajectory is known at each time step. Therefore, an accurate trajectory can be produced
with relatively large time steps compared to the numerical integration scheme, so the total
1
() Ne
() ( )
()
x fx u y
yx
B
e u
y
Plant Dynamics with
General Controller
Relay
Nonlinearity
46
solution time is less. However, this is only true between switching instances. When a
relay switch occurs, the state solution changes due to its piecewise linear nature.
Therefore, numerical techniques are required in order to solve for the state on the
switching surface as discussed below in Sections 4.2.3 and 4.2.4.
4.2.2 Flexible Plant Dynamics
In the case where the plant represents a linear flexible structure, the system can be
represented as in Figure 4.2
Figure 4.2: SISO Regulator Control (r = 0) of Linear Plant with General Feedback
Control Signal to Relay Nonlinearity
The A and B matrices represent the state matrix and input vector respectively. For a
linear controller ( ) x x C , where C represents the output vector. Recall the open-loop
transfer function (3.7) repeated below
2
22
0
() () ()
2
n
ij
I
ppd
j jjj
u
K
Gs K s G s K K s
ss s
(4.1)
where K(s) represents a linear PID controller and G
p
(s) represents a flexible structure
with “n” flexible modes and a rigid-body mode (j=0). The parameters from the frequency
domain model in (4.1) can be used to assemble the A, B, and C matrices for the time
domain model by inspection as follows.
1
() Ne
()
()
x xuy
yx
AB
e u
y
Linear Plant with
General Controller
Relay
Nonlinearity
47
Note that the structure of (4.1) allows the state space to be formulated into a block-
diagonal matrix form. This block-diagonal form will allow the state solution to be
obtained in closed form, which will eliminate the need for numerical techniques for
matrix inversion and evaluating the matrix exponential, thus improving solution
accuracy.
For a P or PD controller
00
11
01
,
nn
n
AB
AB
AB
CC C
AB
C
(4.2)
For rigid-body modes
00 00
0
0 01
,, 0
00
i
i
AB Cu
u
(4.3)
For flexible modes
2
01 0
,, 0
2
ii iij
jjj ij
AB Cu
u
(4.4)
Substituting (4.3) and (4.4) into (4.2) yields for a state-space representation with a P or
PD controller the following
48
0
2
111 1
2
00 1 1
01 0
00
01 0
2
,
0
01 0
2
i
i
nnn in
pi d i pi d i pin d in
u
u
u
Ku K u Ku K u Ku K u
AB
C
(4.5)
where K
p
and K
d
are position and velocity feedback PD controller gains. For a PI or PID
type controller, the above equations are augmented with an extra state to calculate the
integral of position term (see Appendix A).
4.2.3 Poincaré Mapping Surface Definition
The following Figure 4.3 shows how switching instances of the relay in the same
direction, are used to define the Poincaré mapping surface. When the system output
reaches the positive hysteresis boundary with positive slope causing a switch, the state at
this instance is noted (mapped) as an intersection with the Relay Switching Surface.
Figure 4.3: Mapping Surface Defined as System Output Causing Relay Switch in
Negative Direction
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
1
0
1
1
0
1
yt ()
ut ()
t
Intersection 1
(on Poincaré surface)
Intersection 2
(on Poincaré surface)
y(t) system output u(t) = ±M relay output
+ Δ
– Δ
u(t)
M
+ Δ
y(t)
– Δ
Relay with Hysteresis
49
The piecewise constant relay output makes a step change at the relay switching time, and
the state at this instance becomes the inital condition for the trajectory until the next
switch is reached. It is therefore critical to accurately predict the system state at the time
of each switch. The numerical procedure used to evaluate the state on the switching
surface in this work is shown in Figure 4.4. The trajectory is started from an initial state
on the switching surface, and is integrated forward with large time steps until a change in
sign of the switching surface constraint equation ( ( ) 0 Sx in Figure 4.4) signals
that the switching surface has been crossed. The program then backs up to the last state
just before the crossing, and advances forward with a much smaller time step until the
surface is crossed again. This will result in a point on the trajectory that is close to the
switching surface at the time of a relay switch.
4.2.4 Poincaré Mapping Surface Intersections
A one step method developed by Henon [26] is used to accurately calculate both the state
on the switching surface, and the time to reach the surface. It was found that this method
is far more accurate than using linear interpolation between two points on opposite sides
of the surface. Also, since the solution is arrived at in one step, this method is much faster
than using a numerical solver (such as fsolve in Matlab) at each switching instance, and
accruacy is of the same order. The method is outlined in the following. First consider the
original nth order system state equations
50
1
11
1
(, , )
(, , )
n
n
nn
dx
fxx
dt
dx
fxx
dt
(4.6)
with a Poincaré surface defined by the relay switching surface S. For the switching
surface depicted in Figure 4.4
1
(, , ) ( ) 0
n
Sx x x (4.7)
Define a new state variable as the relay switching surface
11
(, , )
nn
x Sx x
(4.8)
1
11 1
1
(, , ),
n
n
nnn i
i i
dx S
fx x f f
dt x
(4.9)
Adding this to the original state results in a new set of augmented state equations
1
11
1
1
11
(, , )
(, , )
(, , )
n
n
nn
n
nn
dx
fx x
dt
dx
fx x
dt
dx
f xx
dt
(4.10)
Next make time t a dependent variable by dividing the first n equations in (4.10) through
by (4.9), and inverting the last equation in (4.10). The new system of state equations
(4.11) are now dependent on x
n+1
, with the augmented state vector
1
[, , ,]
T
n
x xxt .
51
1
11
11
1
1
11
1
11
1
(, , )
(, , )
(, , )
(, , )
1
(, , )
n
nn
n
n
nn
nn
n
nn
n
dx
fx x
f xx
dx
dx
fx x
f xx
dx
dt
fx x
dx
(4.11)
Equations in (4.11) are now solved for one integration step using a one step method, such
as a fourth-order Runge-Kutta, with an integration step size of
1 n
x S
(4.12)
using the state and time just before the switch as intial conditions
010 00
[, , , ]
T
n
x xxt .
The solution yields the state on the switching surface S, along with the time to reach
the surface. The only numerical error introduced in this one integration step is from the
fourth-order Runge-Kutta routine. Therefore, the closer to the surface the initial state is,
the more accurate the solution. Note that for the linear plant case, where the state is
known exactly between switches, this one integration step to each switching surface is the
only source of numerical error aside from round-off errors. For nonlinear plants, where
numerical integration in time is necessary for the entire trajectroy, the error during this
one step is on the order of that encoutered during one time step (Henon [26]).
Once the state on the switching surface is known, it is then used as the initial
condition for integrating forward until a change in sign of the next constraint equation.
Depending on the relay type, there will be a different number of constraint equations
where this procedure is carried out. However, only relay switches in the prescribed
52
direction, caused by relay input with the prescribed slope will be mapped (Figure 4.3) for
the subsequent error calculation (Section 4.2.6).
Figure 4.4: Mapping Surface Defined as Relay Input Causing Switch in Positive
Direction
Switching surface constraints used to develop Poincaré switching surfaces for a relay
with hysteresis are found in Section 3.3. Other relay types developed during this research
can be found in Chapter Chapter 5: and Appendix B.
state trajectory ( )
from plant dynamics
()
x t
xfx u B
() 0
switching surface
(Poincare surface)
Sx
State x just before
switch (large time step)
State x just after switch
(large time step)
Solve for x on
switching surface
using small time
step crossing as
initial guess
State x just before
switch (small time step)
used as initial guess
53
4.2.5 Fixed Point Stability on a General Switching Surface
If the system plant is linear, the fixed point locations can be calculated using methods
found in Section 3.3. The local stability of each potential fixed point is then checked to
eliminate any unstable candidates. The stability method found in Astrom [4] pertains to a
linear switching surface only. As an extension to this method, the following includes
local stability of fixed points located on a general switching surface.
The following two figures describe the physical meaning of the fixed point locations
that will be used in the analysis. Figure 4.5 shows the relation between the system output
y(t), and the relay output u(t) for a relay with dead-zone at fixed point locations
** * * *
01 0
(0) , ( ) and ( ) x xx t x xt x .
Figure 4.5: Stable Limit Cycle Fixed Point Designations for Relay with Dead-Zone
Figure 4.6 shows a phase-plane view of the state trajectory as it intersects the switching
surfaces, under a stable limit cycle. During stable limit cycling, the state trajectory
intersects each fixed point at periodic time intervals, equal to the limit cycle period. Also
*
**
0
2,
() –
tt T
xtx
*
*
0
2,
()
tt T
xTx
0 0.5 1 1.5 2
2
1
0
1
2
+ δ
– δ
*
0
0, (0) tx x
*
**
1
,
()
tt
xtx
() ( ( )) ytxt u(t) = 0, ±M
54
shown is a perturbation
*
0
x of the fixed point
*
0
(0) x x at time zero, and resulting
perturbations of fixed points located on the subsequent switching surfaces. In a stable
limit cycle, the initial perturbations at each fixed point will diminish with each passing of
the state trajectory. Note that this property is used as a convergence criterion in the
mapping process described in Section 4.2.6.
Figure 4.6: Phase Plane View of Stable Limit Cycle Perturbations on Switching Surfaces
for Relay with Dead-Zone
In order to determine the Jacobian matrix W for stability analysis purposes, first
consider the limit cycle trajectory between the switching instances shown in Figure 4.7.
** *
20
() xtx x
x
2
(()) xt
x
1
*
0
() (0) xTx x
*
0
x
*
2
x
(()) xt
x xM AB
x xM AB
Stable
Limit Cycle
Trajectory
*
1
x
**
1
() xtx
*
1
x
x x A
x x A
55
Figure 4.7: Phase Plane View of Stable Limit Cycle Dynamics Between Switching
Surfaces for Relay with Dead-Zone
Since the system is piecewise linear with constant input d
i
between relay switches, the
state equation between switching surfaces i and i-1 and its solution on switching surface i
are given by
ii
i
x xd AB (4.13)
() *
1
0
() e
i
ii
i
ii i
x xe d d
AA
B (4.14)
1,2, , , number of switching surface intersections
relay magnitude between switching surfaces and 1
time duration between switching instance and 1
i
i
iss
dii
ii
These equations can represent Poincaré maps from one switching surface to the next. For
example, (4.14) represents the map from the initial condition at
*
0
(0) x x on the
originating switching surface, to the fixed point
1 *
11
() x x on the destination switching
surface, after time duration τ
1
(see Figure 4.7). It also represents the subsequent mapping
from
*
1
x to
*
2
x .
x
2
(()) xt
x
1
**
0
(0)
s
x xx
(()) xt
11
1
x xd AB
11
1
ss
s
x xd
AB
Stable
Limit Cycle
Trajectory
*
11
() x x
*
1 s
x
22
2
x xd A
ss
s
x xd A
2 **
22 0
() x xx
56
The Jacobian used to determine the limit cycle stability is found from the Poincaré
maps, linearized about their respective fixed points
*
i
x . To linearize each map, take
Taylor series expansions in both location and time of (4.14) (ignore higher-order terms)
** *
11 1
() () (,) (,)
ii i i
i i i x ii i t ii i
xx fx xfx
(4.15)
where
** 1
11
**
*
1
( , ) ( ) system state between switches and -1
and
(, ) and (, )
initial state, unperturbed time between switches
ii
n
ii i i
ii nn nn
xi i t i i
ii
f xt x t i i
xx x R t
fx R f x R
x
In order to ensure that the resulting perturbation ()
i
ii
x of the fixed point ()
i
i
x
in (4.15) also lies on the switching surface, a Taylor series expansion of the system
output ( ()) x t is taken about the fixed point ()
i
i
x
** *
11 1
(( )) (()) ( , ) ( , )
ii ii i i
i i i x ii i t ii i
xxgxxgx
(4.16)
where
** 1
11
**1
11
*
1
( , ) ( ( )) system output between switches and -1
and
(, ) and (, )
initial state, unperturbed time between switches
iii
n
ii i i
ii nn
xi i t i i
ii
gxt x t i i
xx x R t
gx R g x R
x
Noting that the fixed point ()
i
i
x and the perturbation point ()
i
ii
x in (4.16) both
lie on the same switching surface requires (()) (( ))
ii i i
iii
xx , yielding
** *
11 1
(, ) ( , ) 0
ii
xi i i t i i i
gx x g x
(4.17)
57
Equation (4.17) can be rearranged to yield the perturbation in time to reach the next
switching surface, caused by the perturbation on the initial switching surface
*
* 1
1
*
1
(, )
(, )
i
xi i
ii
i
ti i
gx
x
gx
(4.18)
Rearranging (4.15) and using (4.18), the incremental Jacobians W
i
between switching
instances are
** *
11 1
*
** * * 1
11 1 1
*
1
() () (,) (,)
(, )
() ( , ) ( , )
(, )
ii i i
i i i x ii i t ii i
i
ii i
xi i
ix i i t i i i ii
i
ti i
xx fx xfx
gx
x fx f x x x
gx
W
(4.19)
**
1 iii
x x
W (4.20)
where
*
()
i
ii
xx and
*
** 1
11
*
1
(, )
(, ) ( , )
(, )
i
ii
xi i
ix i i t i i
i
ti i
gx
fx f x
gx
W (4.21)
Equation (4.20) provides a linearized Poincaré Map for a fixed point perturbation on
the initial surface, to the resulting perturbation of the fixed point on the destination
surface. These maps can be combined to consider the effect of an initial perturbation
*
0
x
on the originating switching surface, to the perturbation of the fixed point on the final
switching surface as follows
**
1210
*
0
sss
x x
x
WW W W
W
(4.22)
where the Jacobian of the combined successive Poincaré maps is
58
*
** 1
11
*
11
1
(, )
(, ) ( , )
(, )
i
ss
ii
xi i
ixii tii
i
ii
ti i
gx
fx f x
gx
WW (4.23)
and s = the number of switching surfaces intersected by the state trajectory.
From the state response at switching instances given in (4.14), calculate the partial
derivative terms in (4.23), evaluated at
*
1
(, )
ii
x
. Note that x is substituted for
*
1 i
x
, and t
is substituted for τ
i
in (4.14), in order to restore dependence of the equations on x and t.
()
0
(, )
(, ) ( , )
t
i tt
i
ii
f xt e x e d d
gxt xt
AA
B
Derivatives with respect to x –
*
1
(, )
i
i
x ii i
fx e
A
Φ (4.24)
**
11
(, ) ( , )
ii
x ii x i i
gx x
(4.25)
Derivatives with respect to t –
*
1
**
(, ) ()
() using ()
ii
ti i i
ii
ii i i i i
fx x
x dx d x x
ABAB
*
1
(, )
i
ti i i
fxv
(4.26)
**
11
(, ) ( , )
ii
ti i t i i
gx x
(4.27)
where
*
ii i
vx d AB is the state velocity just before the relay switches magnitude from
d
i
to d
i+1
at ()
i
i
x
.
59
Thus for periodic solutions of piecewise linear plant dynamics described by (4.13)
and (4.14), and general output equation ( ) y x , the combined Jacobian (4.23) is
*
1
*
1
1
(, )
(, )
i
s
xi i
ii
i
i
ti i
x
v
x
W Φ (4.28)
where
i
i
e
A
Φ and
*
ii i
vx d AB . Equation (4.28) gives the stability Jacobian W for a
linear plant with s general form switching surfaces.
For a linear switching surface, such as a PD controller, the partial derivative terms in
(4.28) for the ith switching surface are evaluated at
*
1
(, )
ii
x
as follows.
()
0
(()) ()
ii i
t
tt
i
xt xt
ex e d d
AA
C
CB
(4.29)
Derivatives with respect to x –
*
1
(, )
i
i
x ii i
xe
A
CC Φ (4.30)
Derivatives with respect to t –
*
1
**
(, ) ()
()
using ( )
ii
ti i i
i
ii
i
ii i i
xx
xd
xdx x
C
CA B
CA B
*
1
(, )
i
ti i i
x v
C (4.31)
where
*
ii i
vx d AB is the state velocity just before the relay switches magnitude from
d
i
to d
i+1
at ()
i
i
x
.
60
Using (4.30) and (4.31) in (4.28) yields
1
1
s
i
ii
i i
s
i
i
i i
v
v
v
v
C Φ
W Φ
C
C
I Φ
C
(4.32)
where
i
i
e
A
Φ and
*
ii i
vx d AB .
Equation (4.32) gives the stability Jacobian W for a linear plant with s linear
switching surfaces. For a relay with hysteresis only, there is only one switching surface
intersection during the half-period occurring at t = t*. Due to the limit cycle symmetry,
only stability of the half-cycle needs to be evaluated. Therefore, for one switching surface
s = 1 and τ
1
= t*, so (4.32) reduces to
1
1
1
v
v
C
WI Φ
C
(4.33)
where
*
1
t
e
A
Φ and
** *
11 1 0 0
() ( ) ( ) vx d x M x M AB A B AB , using the limit
cycle symmetry constraint
**
10
x x and
1
dM . This matches the corresponding result
found in Astrom [4].
Referring to Figure 4.5, a relay with dead-zone has two switches that occur during the
half-period of the limit cycle. One switch at the reset time γt* and one at t*. Assuming a
linear PD controller and using (4.32), the stability Jacobian for this relay system is
21
21
21
vv
vv
CC
WI Φ I Φ
CC
(4.34)
61
where
***
12
()
12
1 *
11 1 1
2 **
22 2 2 0
**
12 2 0
, ,
() ,
() ,
using , 0, and (symmetry)
ttt
ee e e
vx d x M
vx d x x
dMd x x
AA AA
ΦΦ
ABAB
ABA A
Equation (4.34) was derived assuming a relay with dead-zone and no hysteresis as
depicted in Figure 4.5. However, note that the dead-zone and hysteresis parameters do
not appear explicitly in W. The effect of adding hysteresis is that it shifts switching
surface 1 from
*
(( )) yxt to
*
(( )) ( ) xt , and switching surface 2
from
*
(( )) yxt to
*
(( )) ( ) xt (see Section 5.2). This effect shows up in
the calculations of switching times γt* and t* and the resulting fixed points, but does not
change the form of the Jacobian W in (4.34). Therefore, both (4.34) and the more general
Jacobian for a linear plant with general form controller and multiple switching surfaces
given by (4.28) are valid for a relay with or without dead-zone and hysteresis.
4.2.6 Poincaré Map Integration Stopping Criteria
Note that for a periodic limit cycle to exist, the following mapping holds for consecutive
relay switches.
((1))(()) 0,1,2... xn T GxnT n (4.35)
In the general case, the mapping ( ( )) Gx nT is not known, however, for the relay feedback
case at hand, it is determined numerically as the system is integrated from one relay
switch to the next as shown in Figure 4.8. Since the system trajectory is mapped at each
62
switch, convergence to a periodic solution can be determined by an asymptotic
contraction of the distance between each consecutive intersection.
Figure 4.8: Convergence of Switching Surface Intersections to a Limit Cycle Fixed Point
The following error criterion is checked by the program at each relay switch, and a
determination is made whether to keep integrating for one more switch, or to stop and set
the period for the current cell. This period is defined as the switching time for the last
switch, or set to zero if the trajectory has diverged into the Sink Cell (as is the case for the
unstable solutions). An additional divergence criterion can be set to stop integrating when
the measured output exceeds the small deformation assumption used in constructing the
*
state trajectory ( )
from plant dynamics
()
x t
xfx u B
()
switching surface (Poincare surface)
yx
Error between
successive
switches
Fixed point x*
x
1
(t)
x
2
(t)
x
k
(t)
Trajectory between
successive switches
x
1
, x
2
, x
3
,…, x
k
(convergence)
63
structural model, and maps this cell to the Sink Cell. One more additional setting is the
maximum number of switches. This stops the integration for this cell if the minimum
error criterion is not reached within the maximum number of switches, and maps it to the
Sink Cell. The error convergence criterion is given by
1
() ()
0,1,2...
()
( ) state variable coordinates on the switching plane
number of relay switching instance in proper direction
0 for initial condition at center of originating Cell
nn
x
n
n
xt xt
en
xt
xt
n
n
(4.36)
Figure 4.8 demonstrates convergence of successive trajectory intersections to a cell in
a stable Region of Attraction, based on the above error criteria. The program calculates
an average error using the last 3 intersection points in order to smooth the results, and
prevent halting of the integration routine at the first crossing of the error threshold. This
setting can be increased if needed. Note that (4.36) provides a convergence criterion for a
general plant with a general switching surface, and is not limited to linear plant dynamics
with a linear controller.
4.3 SSCM Procedure
The SSCM procedure is outlined in the following. Four example applications
demonstrating the method in detail are included in Chapter Chapter 5:. For linear plants,
limit cycle frequencies and fixed points can be determined by methods in Step I, prior to
the cell mapping process in Step II. The advantage of knowing the fixed point locations
ahead of time is that the cell map boundaries can be established on a region enclosing the
fixed point locations.
64
For nonlinear plants, the limit cycle fixed points cannot in general be determined
before hand. Therefore, a reasonable estimate of the region of interest to be mapped must
be made using other methods. A good starting point is to use a coarse mapping, which
encompasses the system trajectory bounds for the specific design problem. This will
point out any clustering of the fixed points on the switching surface, which can be used to
focus the mapping to a smaller area. Even with a relatively coarse cell map, the unknown
fixed points will be accurately mapped using the integration stopping criterion, and
corresponding limit cycle periods will be calculated using the time between the last relay
switch.
I. Find Limit Cycle Characteristics and Periodic Fixed Points
a) Obtain undamped natural frequencies and mode shapes of discretized
flexible structure. These can be derived from lumped mass models,
assumed modes models, or Finite Element software modal files.
b) Assemble state equations in block-diagonal form shown in Section 4.2.2,
using assumed modal damping ratio and undamped natural frequencies
and modal matrix coefficients obtained in a)
c) Find limit cycle frequencies and corresponding periodic fixed points using
half-period t* locus of constraint equation(s) of the form (3.26). Number
of constraint equations depends on relay type.
d) Check stability of fixed points using eigenvalue test (3.30)
e) Check half-period waveform per (3.27) via simulation from fixed point
initial conditions
65
II. Map Region of Attraction for Fixed Points Found in Step I
a) Determine relay switching surface ( ) ( ) yx to be mapped,
using system output equation ( ) y x and relay dead-zone δ and/or
hysteresis Δ.
b) Establish 1 n independent state coordinates x
i
, 1, ,, inik to be
mapped by the 1 n dimensional cell space, using the switching surface
constraint in a). Note that the resulting dependent state x
k
should be chosen
as an explicit function of the other states.
c) Determine number of divisions N along each independent state dimension
x
i
, using boundaries enclosing fixed point locations from Step I.
d) Create
1 n
N
dimensional cell space
e) Correlate independent state coordinates x
i
of initial condition points with
cell center locations from d). Dependent state coordinate x
k
for initial
condition points is found using switching surface constraint equation in b).
f) Perform cell mapping process by integrating system equations from each
initial condition in e), stopping when
i. Error criterion (4.36) is reached
ii. Maximum number of switches is reached
iii. Maximum spatial criterion is reached (e.g. max spring deflection)
g) For each originating cell center in f), note period between last relay switch
and compare to known limit cycle periods found in Step I. Assign this
66
period to cell. Cell period is assigned as zero (Sink Cell) for integrations
stopped using criteria f)ii or f)iii.
h) Plot initial condition points displaying limit cycle regions of attraction
differentiated by point colors assigned by limit cycle period
Figure 4.9 shows a flowchart summarizing the major tasks performed when
implementing the method. Once a cell map is obtained, it can then be used to investigate
the region of attraction for each stable limit cycle, determining robustness of the design
towards desirable or undesirable limit cycles. Model parameters can be adjusted, and the
process repeated, to determine the effect of parameter variation on limit cycle robustness.
Figure 4.9: Limit Cycle Region of Attraction Mapping Procedure
Assemble State Equations and
Feedback Equation
Find LC Fixed Points
for Region of Interest
Discretize Switching Surface
in Region of Interest
Integrate from selected IC until max
location error or max intersections
Note mapped Cell location and time
between last intersection (LC period)
Pick next
unprocessed
Cell
Assemble State Equations and
Feedback Equation
Find LC Fixed Points
for Region of Interest
Discretize Switching Surface
in Region of Interest
Integrate from selected IC until max
location error or max intersections
Note mapped Cell location and time
between last intersection (LC period)
Pick next
unprocessed
Cell
Step II
Step I
67
4.4 SSCM Summary
The purpose of the SSCM is to provide a tool for evaluating the effect of structural
flexibility on resulting limit cycle performance. If the design requirement is a limit cycle
of prescribed frequency and amplitude, the methods in Step I can be used to study how
changes in plant parameters (mass, stiffness, modal content), feedback gains, and relay
parameters affect the resulting limit cycle(s). If multiple limit cycles exist, these
parameters can be adjusted to provide a dominant limit cycle with the desired limit cycle
characteristics. For linear plants the following limit cycle characteristics are calculated
limit cycle frequencies, amplitudes, and waveforms
periodic fixed points on switching surfaces
switching surfaces may be linear or nonlinear
For nonlinear plants, use cell mapping results. Trajectories originating from cells that
settle on the same steady-state limit cycle will converge to the same destination cell. All
originating cells will possess the same Period characteristic measured at the time of the
last switch before convergence. This is the limit cycle period for this group of cells, and
the fixed point location will be the destination cell shared by all of the originating cells in
this group.
If a specified limit cycle is desired and multiple limit cycles exist, the degree of
robustness of the dominant limit cycle can be qualitatively determined by the region of
attraction map provided by Step II. A large region of attraction that covers the area of
expected state-space performance, provides assurance that bounded perturbations will not
cause undesired limit cycles to be inadvertently excited.
68
If the design requirement is that no limit cycles be excited, methods in Step I provide
analytical techniques for determining if a limit cycle exists that are superior to the
Describing Function Method. Additionally, the cell mapping procedure in Step II can be
run on the final design to verify that no periodic motions exist within the region of
interest. All cells should map to the Sink Cell.
Chapter Chapter 5: contains four examples where application of the method is
described in detail. The examples include a linear plant with linear switching surface, an
example with a 2-level relay, an example with a nonlinear switching surface, and finally a
nonlinear plant with a linear switching surface.
69
Chapter 5: Applications
The following sections present four example problems that demonstrate the method
presented in Chapter Chapter 4:. First, a system of 2 masses connected by a spring is
controlled by a constant force input from a relay with hysteresis, which is triggered by a
linear feedback signal. Example 2 consists of a 2 mass-spring system controlled by a 2-
level relay with dead-zone and hysteresis, triggered by a linear feedback signal. Example
3 consists of a slewing beam with constant level torque input from a relay with dead-zone
and hysteresis, which is triggered by a nonlinear feedback signal. The nonlinear feedback
signal contains a signed velocity squared term. A new stability Jacobian is developed in
this section, since the switching surface is no longer a linear combination of the states.
Lastly, a nonlinear plant is considered, consisting of a 2 mass-spring-damper plant, with a
nonlinear spring. The nonlinear spring is a softening type spring, which creates an
unstable plant, when the spring exceeds a certain length.
70
5.1 Two-Degree-of-Freedom System Under Relay Control with
Hysteresis
5.1.1 Plant Description
Figure 5.1 shows an unconstrained 2 mass-spring plant with one rigid-body and one
flexible mode. Control force is applied to mass 1 according to (5.1).
Figure 5.1: 4-State System with 1 Rigid-Body Mode and 1 Flexible Mode Under Relay
Control with Hysteresis and PD Feedback
The equations of motion are first put into modal form per (3.2), and then arranged in
block-diagonal state variable form per (4.5), with the controller included in the output
equation as shown in Figure 5.2. This is done to define the input to the relay as a function
of the system output y(t). Note that the system input r = 0 in the limit cycle analysis under
consideration.
k
1
() qt
2
() qt
11 1 1
22 2
0
0 0
mq q kk F
mq q kk
12
1, 2, 10 mm k
11
()
pd
ekq kq
Collocated PD
feedback at m1
1, 0.1
pd
kk
Force on m1
1, 0.1 D
e
u
PD Controller with On-Off Actuator with
hysteresis
71
Figure 5.2: Rearranged Block Diagram with System Input r = 0
The resulting state equations for the system are
if or and ( )
()
if or and ( )
()
M yt y ut M
uy
Myt y ut M
x xuy
yx
AB
C
(5.1)
The system input u(y) represents a relay nonlinearity with hysteresis Δ, and magnitude M,
based on mass 1 position and velocity feedback signal
11 pd
y ekq kq .
The resulting eigenvalues and normalized eigenvectors used to assemble the A, B,
and C matrices defined in (4.5) for this plant are as follows, where modal damping of 1%
( 0.01) has been assumed.
0.5774 0.8165
0, 3.8730,
0.5774 0.4082
n
U (5.2)
01 0 0 0
0 0 0 0 0.5774
,
00 0 1 0
0 0 15.0000 0.0775 0.8165
0.5774 0.0577 0.8165 0.0186
AB
C
(5.3)
1
() Ne
() x xuy
yx
AB
C
e u
y
Linear Controller and
Flexible System
Relay
Nonlinearity
72
5.1.2 Limit Cycle Frequencies and Fixed Points
In order to determine the limit cycle frequency using the state-space equations, it is
necessary to calculate the time at which the relay switches under limit cycling conditions.
It is assumed that the limit cycle is symmetric, with only two switches per cycle. These
are the characteristics of a unimodal periodic motion, and the waveform and input/output
relations are shown in Figure 5.3.
+ Δ
– Δ
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
1
0
1
1
0
1
Limit Cycle Waveform (Relay Input u(t), Linear Plant O
t
yt ()
ut ()
t ()
t ()
t
t = t* = T/2
t = 2t* = T
t = 0
Figure 5.3: Symmetric Unimodal Limit Cycle Waveform in Relay Feedback System
As shown in the preceding derivation of equation (3.26) repeated below, at t = 0 the
relay input crosses the positive hysteresis boundary yielding
*
*
1
0
(0) e
t
t
yIedM
AA
CB
(5.4)
Closed form solutions of the matrices in (5.4) are found in Appendix A. This allows for a
more accurate determination of t* and x*, especially for higher-order systems. Solutions
to equation (5.4) can be found in the frequency domain by substituting t* = π/ ω and
rearranging to yield the “State-Space Locus”
73
1
0
() e 0 fIedM
A A
CB (5.5)
A plot of the magnitude versus frequency can then be constructed as shown in Figure 5.5.
Solutions occur where the locus crosses the frequency axis with a positive slope. The
positive slope solutions correspond to positive velocity (dy/dt) at the switching instance at
y(0), which is a criterion for proper waveform defined in (3.27).
The solutions to (5.5) shown in Figure 5.5 for this example system are as follows
0.5074, 0.6132, 0.8344, 1.3932, 4.4959
c
(5.6)
Note that the DFM only predicts 2 solutions at 0.415 and 4.45 rad/sec as shown in Figure
5.4, and by the noted curve in Figure 5.5.
Figure 5.4: Nyquist Plot Used to Estimate Limit Cycle Candidates
()
p
KG
1() Na
ω
a
ω
1 ω
2
74
Figure 5.5: System Frequency Locus Used to Determine Limit Cycle Candidates
From a design standpoint, there are two ways to change the limit cycle frequencies as
well as the number of low frequency solutions. Notice in the right hand plot in Figure 5.5,
that there are no additional downward spikes in the area of the high frequency solution,
however, in the low frequency area there are several. If the locus is shifted up or down,
the number of intersections with the time axis in the low frequency region will increase if
the curve is shifted up, and decrease if shifted down. Referring to Equation (5.5), an
equivalent up or down shifting of the locus can be achieved by increasing or decreasing
the relay hysteresis as shown in Figure 5.6. Note that the stability of the limit cycles may
also change as noted in Figure 5.7, which shows the solutions and their stability for
various hysteresis values depicted in Figure 5.6.
High freq
solution
4 Low freq
solutions
1 DFM Low freq soln
75
Figure 5.6: System Frequency Locus Intersections for Various Hysteresis Values
Figure 5.7: Relay Hysteresis vs Normalized Frequency Solutions ( ω
c
/ ω
n
)
Another way to change the limit cycle frequencies and number of low frequency
solutions is to vary the controller gains. Figure 5.8 shows frequency loci at the nominal
gains as well as ½ nominal gains and 2X nominal gains. The hysteresis is held constant.
+ = stable
x = unstable
Δ =
0.050
0.100
0.150
0.200
Frequency increases and number of low
freq solutions decreases with decreasing Δ.
Nominal
Δ = 0.100
76
Figure 5.8: System Frequency Locus Used to Determine Limit Cycle Candidates
Local stability for each limit cycle at these frequencies is tested using the Jacobian of
the Poincaré map given by equation (3.30) repeated below
*
t
v
e
v
A
C
WI
C
(5.7)
where
*
(0) v x Ax BM . The limit cycles are locally stable if the eigenvalues of
(5.7) lie within the unit circle. The largest eigenvalue for each frequency for the nominal
system is
max
0.8138, 0.8545, 0.9983, 1.0129, 0.9751 (5.8)
.1
for all curves
Increasing K
High frequency solutions
Low frequency solutions
Increasing K
Nominal
k
p
=1, k
d
=.1
(blue curve)
f( ω) vs ω
77
Note that the magnitude of the eigenvalue corresponding to the fourth frequency is larger
than unity, so this limit cycle is unstable. Therefore, the stable limit cycle frequencies and
their corresponding half-period times t* are
0.5074, 0.6132, 0.8344, 4.4959 rad/sec
LC
(5.9)
*
= 6.1915, 5.1231, 3.7650, 0.6988 sec
LC
t
(5.10)
Using (3.25) repeated below, the solutions for t* (5.10) yield the following fixed
points of the stable limit cycles.
*
*
1
*
0
(0)
where e and
t
t
xIMx
ed
AA
ΦΓ
ΦΓ B
(5.11)
*
0.0000 0.0000 0.0000 0.0000
1.7873 1.4789 1.0869 0.2017
, , ,
0.0095 0.0066 0.0129 0.0133
-0.1345 0.1125 0.3276 0.9493
x
(5.12)
Figure 5.9 shows the 4 limit cycle waveforms started from the above fixed point initial
conditions.
78
Figure 5.9: Four Stable Limit Cycle Solutions Started From Initial Conditions
System Output y(t) (blue), Relay Output u(y) (red)
5.1.3 Limit Cycle Fixed Point Region of Attraction Mapping
As shown in Section 3.4.2, if the cell space lies on a Poincaré surface, the dimension of
the cell space is n – 1 for an n dimensional system. Thus for the 4 state system at hand,
the cell space is 3D. For visualization purposes, it is helpful to place other constraints on
the cells in the area of interest, so that the resulting map is 2D. The following 2D slices
through the 3D cell space are taken where the state trajectory intersected the switching
surface, during a simulation of the system. The simulation produces various limit cycles,
which are excited with pulse force inputs to the system as shown in Figure 5.12. The idea
is to see if the map can be used to predict a steady-state limit cycle, based on the first
switching surface intersection in the proper direction of an arbitrary state trajectory.
79
Figure 5.10 shows the system state trajectory started from rest as it intersects the
switching surface in the proper direction for the first time at the location noted in Figure
5.12a. Figure 5.12a shows a simulation of the system starting from rest. The 2D map is
centered at this location. Figure 5.11 shows the x1-x3 projection of the map in Figure
5.10, where the first intersection lands in the high frequency limit cycle region of
attraction (red Xs). The high frequency limit cycle is initially excited as predicted by the
map. Pulse disturbances at the system input at 1000 second intervals starting at 500
seconds show that other stable limit cycles can be excited also (Figure 5.12b-d). These
limit cycles can also be predicted using the map, given the system state on the switching
surface immediately following the disturbance (see Figure 5.13 – Figure 5.15).
State Trajectory
from Zero State
Switching Surface for first
“+ to –“ Relay switch of
state trajectory from Zero
Figure 5.10: End View of Point Mapping Surface and Zero State Startup Transient
Trajectory at 1st Intersection (location “A” in Figure 5.12)
80
Note that 1
st
intersection
with Cell Map
surface is in the
High Freq LC
ROA (red Xs),
which was
verified by
simulation
Figure 5.11: Region of Attraction Map Showing Trajectory Intersection “a” (magenta
point) in Figure 5.12 Located in LC4 ROA (redXs): ROA for LC1=blueX, LC2=greenX,
LC3=cyanX, LC4=redX
81
b)
LC4
Tp4=1.40 sec
Amp4=.198
c)
LC2
Tp2=10.25 sec
Amp2=1.03
d)
LC1
Tp1=12.38 sec
Amp2=1.61
a)
1
st
switch in
proper direction
after startup
.068506
.166001
.122508
.366553
x
LC4 LC1 LC2 LC1 LC2
b)
c)
d)
a)
Figure 5.12: Simulation Showing Different Limit Cycles Excited by Pulse Disturbance at
Relay Output at 1000sec Intervals Starting at 500sec
82
Figure 5.13b-c shows intersection points right after the pulse disturbances that are
mapped in Figure 5.14 and Figure 5.15. Note that the correct ROA is intersected.
a)
b)
c)
b)
1
st
switch after
disturbance
.015708
.658610
.039436
1.042322
x
a)
1
st
switch after
startup
.068506
.166001
.122508
.366553
x
c)
1
st
switch after
disturbance
.170644
1.861206
.081152
.303784
x
LC4 LC1 LC2
Figure 5.13: Relay Input Trajectory Transient and Mapped State at First + to –Switch
After Pulse Disturbance: a) zero-state startup, b) pulse at 500 sec, c) pulse at 1500 sec
83
1
st
intersection
with mapping
surface after
Disturbance at
500 sec lies in
LC2 ROA
(green Xs)
Figure 5.14: Map Showing Intersection “b” (magenta circle) in Figure 5.13 Located on
Boundary of LC2 ROA (greenXs): LC1=blueX, LC2=greenX, LC3=cyanX, LC4=redX
1
st
intersection
with mapping
surface after
Disturbance at
1500 sec lies
in LC1 ROA
(blue Xs)
Figure 5.15: Map Showing Intersection “c” (magenta circle) in Figure 5.13 Located on
Boundary of LC1 ROA (blueXs): LC1=blueX, LC2=greenX, LC3=cyanX, LC4=redX
84
5.1.4 Discussion of Results
This example demonstrated a step-by-step application of the SSCM, and how it can be
used to design limit cycle frequencies by varying relay parameters and/or control gains.
The accuracy of the calculated fixed points was demonstrated by simulations showing the
waveform of each limit cycle, started from its fixed point initial conditions.
The cell map boundaries in this example were chosen based on simulation data, in
order to verify that the map could predict the correct steady-state limit cycle, given an
arbitrary crossing of the switching surface in the proper direction. It correctly predicted
each limit cycle, even at locations that were on the border of adjacent ROAs. This
validates the accuracy of the map even for a relatively coarse cell map (41
3
cells). This
accuracy is due to the SSCM error stopping criterion, which continues integration from
each initial condition in the cell map, until the periodic fixed point has been located
within a specified accuracy. This accurately maps the initial point to the correct limit
cycle fixed point.
85
5.2 Multilevel Relay with Dead-Zone and Hysteresis
5.2.1 Plant Description
Figure 5.16 shows an unconstrained 2 mass-spring plant with one rigid-body and one
flexible mode. Control force is applied to mass 1 according to (5.13).
Figure 5.16: 4-State System with 1 Rigid-Body Mode and 1 Flexible Mode Under 2-level
Relay Control with Dead-Zone and Hysteresis
The resulting state equations for the system are
() x xuy
yx
AB
C
(5.13)
k
1
() qt
2
() qt
11 1 1
22 2
0
0 0
mq q kk F
mq q kk
12
2, 4, 100 mm k
11
()
pd
ekq kq
Collocated PD
feedback at m1
1, 0.1
pd
kk
PD Controller with 2-level Relay
with dead-zone and hysteresis
Force on m1
11 1
22 2
.1, .02, 1
.2, .05, 3
M
M
δ
1
δ
2
M
1
M
2
u
Δ
2
Δ
1
e
86
where
11
11
11
11
if ( ) ( )
or ( ) ( )
and ( )
() if ( ) ( )
or ( ) ( )
and ( )
0otherwise
1,2,...
ii i i i
ii i i
i
ii i ii
ii i i
i
Myt
yt
ut M
uy M y t
yt
ut M
i
, (number of magnitude levels )
i
pM
The system input u(y) represents a relay nonlinearity with multilevel dead-zone δ
i
,
hysteresis Δ
i
, and magnitude M
i
, based on mass 1 position and velocity feedback signal
11 pd
ye kq kq .
The resulting eigenvalues and normalized eigenvectors used to assemble the A, B,
and C matrices defined in (4.5) for this plant are as follows, where modal damping of 1%
(0.01) has been assumed.
0.4082 0.5774
0, 8.6603,
0.4082 0.2887
n
U (5.14)
01 0 0 0
0 0 0 0 0.4082
,
00 0 1 0
0 0 75.0000 0.1732 0.5774
0.4082 0.0408 0.5774 0.0577
AB
C
(5.15)
87
5.2.2 Limit Cycle Frequencies and Fixed Points
In order to determine the limit cycle frequency using the state-space equations, it is
necessary to calculate the time at which the relay switches under limit cycling conditions.
It is assumed that the limit cycle is symmetric, with 4p switches per cycle, where p is the
number of relay levels. These are the characteristics of a unimodal periodic motion, and
the waveform and input/output relations are shown in Figure 5.17.
Figure 5.17: Symmetric Unimodal Limit Cycle Waveform in 2-level Relay Feedback
System
For a symmetric unimodal limit cycle to exist in the system given in (5.13), the
following must hold.
( ) ( ) the motion is periodic with period
( ) ( 2) the motion has odd symmetry
yt yt T T
yt yt T
(5.16)
2
() uy M –M
1
M=0 () yt
22
22
11
11
t = 2t* = T
t = γ
1
t*
t = γ
2
t*
t = γ
3
t*
t = t* = T/2
88
Equation (5.16) provides constraints on the system motion for establishing the periodic
conditions required for the limit cycle to exist. Using the relations in (5.16), the initial
conditions of the state must satisfy
*
() (0) xtx (5.17)
where t* is the limit cycle half-period shown in Figure 5.17.
For a multilevel relay with dead-zone and hysteresis, the method presented in Astrom
[4] is extended to include multiple levels of switching as follows. Since the system is
piecewise linear with constant input between relay switches, the state equation and its
solution are given by
i
x xd AB (5.18)
*
1
*
1
() *
1
0
() e ( )
i
i
tt
tt
ii
xtxt edd
A A
B (5.19)
0
**
1
1, 2, , , number of switches, 0
relay magnitude while
i ii
iss
dttt
For the 2-level relay shown in Figure 5.17, equation (5.18) yields the following state
dynamics between switches
*
21
**
11 2
**
23
**
13
for 0
for
for
for
x xM t t
x xM t t t
xxttt
x xM t t t
AB
AB
A
AB
(5.20)
Since the state is known at each relay switching time t = γ
i
t*, the state solution (5.19) at
these instants in time and (5.17) can be used to determine the initial conditions that result
in a limit cycle
89
*
1
*
1
**
21
**
21
**
32
**
3
**
3
*
12
0
() **
21 1
0
() **
32
() **
31
0
() e (0) ( )
() e ( ) ( )
() e ( ) 0
() e ( ) ( ) (0)
t
t
tt
tt
tt
tt
tt
xt x e d M
xt x t e d M
xt x t
xtxt edMx
A A
A A
A
A A
B
B
B
(5.21)
Combining the equations in (5.21) yields the limit cycle fixed point initial condition
** *
12
1
(1 ) (1 )
21
(0)
tt t
x eI e M e M
AA A
g1 g2g1 τg3
ΓΓΓ
(5.22)
** ***
3 121
00 0
where , ,
tt ttt
ed ed ed
AA A
g1 g2g1 τg3
Γ B Γ B Γ B
Also, referring to Figure 5.17, at t = 0 and subsequent switching instants, the relay input
crosses the hysteresis bands as follows:
22
**
11 22
**
22 11
**
33 11
(0) (0) ( )
() () ( )
() () ( )
() () ( )
yx
yt x t
yt x t
yt x t
C
C
C
C
(5.23)
Solutions of γ
1
, γ
2
, γ
3
, and t* which satisfy (5.23) are obtained numerically, yielding limit
cycle frequencies of ω* = π/t*. Note that an initial guess is required for the parameters γ
i
.
For the assumed symmetric unimodal limit cycle waveform, these parameters obey the
following constraints
12 3
01
(5.24)
90
The following additional waveform criterion must be satisfied in order for the solution to
be a limit cycle (see Figure 5.17).
**
12
**
33
To ensure proper velocity direction at switches
( ) (0) 0, ( ) ( ) 0, ( ) ( ) 0,
( ) () 0, ( ) () 0,
i y ii y t iii y t
iv y t v y t
(5.25)
*
22 1
**
22 1 1 1 2
**
11 11 2 3
**
11 2 2 3
To ensure proper waveform between switches
() ( ) ( ), 0
() ( ) () ( ),
()() () (),
() ( ) () ( ),
iyt tt
ii y t t t t
iii y t t t t
iv y t t t t
(5.26)
Limit cycle stability is verified by extending the Jacobian developed in Astrom [4] to
a multilevel relay with dead-zone and hysteresis. Stability of the limit cycle is satisfied
using the Jacobian of successive Poincaré maps given by
*
121
, 1,2, ,
i
ss
t
i
i
i
v
Ie i s
v
A
WWW WW
C
W
C
(5.27)
where s is the number of switching surfaces,
**
1
()
ii i
tt
is the time between switch i
and i-1, and v
i
is the velocity at switch i at
*
i
tt given by (5.20). The limit cycles are
locally stable if the eigenvalues of W lie within the unit circle. A detailed derivation of
(5.27) is found in Section 4.2.5.
Using (5.27), the Jacobian for a linear switching surface and 2-level relay with dead-
zone and hysteresis is
43 2 1
43 2 1
43 2 1
vv vv
vv vv
CCC C
WI Φ I Φ I Φ I Φ
CC CC
(5.28)
91
22
22
11
where initial condition fixed point on switching surface ( ) ( )
relay reset fixed point on switching surface ( ) ( )
relay reset fixed point on switching surface ( ) ( )
aa
bb
cc
d
** **
32 3 121
11
() (1 ) ()
12 3 4
12 2 1 1
21
relay reset fixed point on switching surface ( ) ( )
, , ,
state velocity just before reset from to at ( )
stat
tt tt
d
ee e e
vb M M M x
vc M
AA AA
Φ ΦΦΦ
AB
AB
12
3 13
41 1 2
e velocity just before reset from to 0 at ( )
state velocity just before switch from 0 to + at ( )
state velocity just before reset from + to + at (0 )
Mx
vd M x
va M M M x
A
AB
The Describing Function Method (DFM) can provide a preliminary estimate of the
frequency range for existing limit cycles. One caveat is that the DFM is approximate, and
may predict no limit cycles, when they do exist. Note that the 2-level relay being
considered can be represented as two 1-level relays in parallel. For nonlinearities in
parallel, an equivalent DF can be assembled by summing the individual DFs (Atherton
[5]). For a relay with dead-zone and hysteresis, the DF is
22
22
22
() 0,
24
() ,
Na a
MM
Na a a j a
aa
(5.29)
where δ is the relay dead-zone, Δ is the hysteresis band, M is the relay magnitude, and a
is the amplitude of the sinusoidal output of the first-order approximation of the relay.
For the 2-level relay shown in Figure 5.16, the combined DF is then
12
22
22 1 11
11 1 1 2 2
22
21 2 2 21 2
22 2 2 2 2
() ()
24
24
Na N a
MM
aa j
aa
MM MM
aa j
aa
(5.30)
92
where δ
i
, Δ
i
, and M
i
are per the relay sketch shown in Figure 5.16. Note that due to the
dead-zone constraint in (5.29), both terms in (5.30) are zero when
11
a and the
second term remains zero until
22
a .
Figure 5.18 shows a modified Nyquist plot of the frequency response of the plant
(5.15), and the amplitude sweep of -1/DF for the 2-level relay (5.30). The DFM predicts
two stable limit cycles (intersections of blue curves and red curve in Figure 5.18). One at
0.872 rad/s (Tp=7.21 sec), which only triggers the first level of the relay, and one at 1.20
rad/s (Tp=5.24 sec), which triggers both levels. This preliminary analysis provides a time
span for the half-period time sweep of t*, used in the numerical analysis to find solutions
to the constraint equations (5.23).
Figure 5.18: Describing Function Plot to Estimate Limit Cycle Frequency Range
The results of the time sweep analysis show 5 stable limit cycle periods of Tp = 6.23,
7.12, 4.82, 4.89, and 5.22 seconds. Note that two of these periods are relatively close to
the DFM predictions. Table 5.1 shows the relay reset time parameters for each limit
cycle. These parameters were solved for by sweeping the half-period time parameter t*
ω
1
=.872
rad/s
amp
1
amp
2
ω
ω
2
=1.20
rad/s
93
from 0 to 5 seconds, and using the Matlab command fsolve on (5.23) at each trial time for
t*, along with an initial guess for each γ
i
. For LC1 and LC2, only the first level of the
relay is active, due to the amplitude of the limit cycle. The guess used was γ
2
=0.5. For
LC3-LC5, the guess used was γ
3
=0.8, γ
2
=( γ
3
)
2
, and γ
1
=( γ
3
)
3
. The range for the time sweep
was arrived at using the limit cycle period estimates (t*=Tp/2) found by the DFM results
shown in Figure 5.18. The resulting solutions for the reset time parameters are then used
in (5.21) to determine the fixed points on each switching surface. The waveform criteria
(5.26) and stability criterion (5.28) for each limit cycle are then verified to determine the
stable fixed points. The fixed points for each limit cycle are x
1
=[0.2461, -0.3403, -0.0025,
-0.0891]
T
, x
2
=[0.2546, -0.5308, -0.0072, 0.0323]
T
, x
3
=[0.3085, -0.7739, -0.0303,
0.0545]
T
, x
4
=[0.2371, -0.8688, 0.0095, 0.2278]
T
, and x
5
=[0.2718, -1.0959, 0.0041,
0.1973]
T
. Figure 5.19 shows the limit cycle waveforms started from fixed point initial
conditions. Note that only Relay Level 1 is active during LC1 and LC2.
Table 5.1: Relay Reset Time Parameters (Figure 5.17)
LC γ
1
γ
2
γ
3
t*(sec)
1 -- 0.5355 -- 3.1136
2 -- 0.7305 -- 3.5597
3 0.4122 0.5053 0.7575 2.4116
4 0.4636 0.6041 0.7906 2.4451
5 0.6023 0.6907 0.8375 2.6088
The above fixed points are located on the switching surface
*
211
() ( ) yt (see
Figure 5.17). This surface was selected for mapping purposes in this example, because all
94
5 limit cycle trajectories cross this surface. As stated above, any switching surface can be
used for mapping purposes.
Figure 5.19: Five Stable Limit Cycle Solutions Started From Initial Conditions
System Output y(t) (blue), Relay Output u(y) (red)
5.2.3 Limit Cycle Fixed Point Region of Attraction Mapping
Figure 5.20 shows the limit cycle trajectories and a 2D slice of the 3D switching surface
cell map in x1-x2-x3 space. The 3D cell map is centered at the cluster of limit cycle fixed
points found in Section 5.2.2. The 2D slice is for visualization purposes, and is taken at
the x
1
coordinate of the center of the fixed point cluster. Figure 5.21 shows the x2-x3
projection of the map in Figure 5.20, as well as a test point for verification purposes. The
95
test point is located in the ROA of LC3 (magenta Xs), and was chosen from a small
region of attraction area on the map, in order to verify the map’s accuracy at predicting
the resulting limit cycle. Figure 5.22 shows a simulation of the system starting from the
test point, and verifies that LC3 is excited from this location on the switching surface.
Figure 5.20: Region of Attraction Map at x
1
=0.2636 and Limit Cycle State Trajectories
for LC1=blueX, LC2=greenX, LC3=magentaX, LC4=cyanX, LC5=redX
LC5
LC4
LC3
LC2
LC1
Full switching surface map is 3D (n-1). Slice taken
at x
1
= 0.2636 (center of fixed point cluster).
96
Figure 5.21: Region of Attraction Map at x1=0.2636 ROA for LC1=blueX, LC2=greenX,
LC3=magentaX, LC4=cyanX, LC5=redX
Figure 5.22: Simulation Showing LC3 Excited from Initial State Noted on ROA Map in
Figure 5.21
Initial state at x = [0.2636 -0.6767 0 0]
T
(rigid-body motion only) will result in
LC3 (magenta Xs), which was verified by
simulation (see Figure 5.22).
97
Limit cycles started from initial conditions that are not on the switching surface can
be predicted by the ROA that is first intersected by the state trajectory in the proper
direction. Since the steady-state convergence has already been determined for all cells on
the switching surface, only the first intersection of the state trajectory in the proper
direction needs to be considered for limit cycle determination.
As noted in Section 3.4.2, if the cell map is constructed with initial conditions that do
not lie on the switching surface, the dimension of the cell space is no longer less than the
dimension of the state space. However, the increase in total number of cells can be
reduced, if only specific modal initial conditions are investigated. In the following case,
only the rigid-body mode initial conditions are mapped, in order to investigate which
limit cycle will be excited with various combinations of rigid-body initial conditions of
displacement and velocity. Note that for this case both masses are initially assumed to
have the same position and velocity. All other modal initial conditions are set to zero
(x(0) = [x
1
(0) x
2
(0) 0 0]
T
). Thus the cell map of initial conditions is only 2D. This would
be the case if the controller is turned on when the flexible modes are not initially excited
(no initial spring deflection or relative velocity between the masses). Note that all modes
are still included in the subsequent dynamics.
The following Figure 5.23 shows a phase plane representation of the regions of
attraction for the limit cycles, when initiated from rigid-body mode initial conditions.
One initial condition in each ROA is selected, and simulation results are shown in Figure
5.24. Compare these results to those in Figure 5.19.
98
Figure 5.23: Rigid-Body Mode Region of Attraction Map: ROA for LC1=blueX,
LC2=greenX, LC3=magentaX, LC4=cyanX, LC5=redX
Black Xs represent maps to Sink Cell (zero velocity ICs
started within relay dead-zone resulting in no motion).
Initial state at a) x = [0 0.2817 0 0]
T
(LC1);
b) [0 0.4557 0 0]
T
(LC2); c) [0 0.6690 0 0]
T
(LC3);
d) [0 0.7746 0 0]
T
(LC4); e) [0 0.9507 0 0]
T
(LC5),
verified by simulation (see Figure 5.24).
99
a)
LC1
Tp1=6.23 sec
Amp1=.164
b)
LC2
Tp2=7.12 sec
Amp2=.245
c)
LC3
Tp3=4.82 sec
Amp3=.271
d)
LC4
Tp4=4.89 sec
Amp4=.287
LC5
Tp5=5.22 sec
Amp5=.365
e)
Figure 5.24: Simulation Showing Limit Cycles Excited from Rigid-Body Initial
Conditions Noted on ROA Map in Figure 5.23
100
5.2.4 Discussion of Results
This example extends the previous example to include both a 2-level relay and a dead-
zone with hysteresis. These result in additional switching surface constraint equations,
which were derived for the general case of multiple switching surfaces. These additional
constraints also require a new stability Jacobian, which was derived for the general case
of multiple switching surfaces.
The cell map boundaries in this example were chosen based on the center point of the
cluster of the calculated fixed point locations. A 2D slice was taken at this location and a
test point location was selected from a small ROA area on the 2D map. A simulation
from this test point verified the accuracy of the map at predicting the steady-state limit
cycle from this ROA. Additionally, an alternate 2D map was constructed with points not
lying on the switching surface, to demonstrate how to handle such a case. Initial
conditions representing the rigid-body mode (x(0) = [x
1
(0) x
2
(0) 0 0]
T
) were chosen, since
only a 2D surface of initial condition points needs to be constructed for this case. The
resulting map correctly predicted the 5 limit cycles from locations on the x2 axis (initial
velocities).
101
5.3 Slewing Beam with Nonlinear Feedback
5.3.1 Plant Description
The following example demonstrates the method using a simplified model of a
spacecraft, represented by a slewing beam system with colocated actuator and sensor at
the base. System damping of 1% is applied as modal damping. The controller has a
nonlinear signed velocity squared term, used in Time-Optimal attitude control of
spacecraft (Sidi [39]). The beam length is assumed to be large compared to dimensions of
base and tip inertias, so the center of mass can be approximated as the attachment point.
Torque input at the base is an on-off relay type with dead-zone and hysteresis.
Figure 5.25: Pinned-free Beam with Rigid-Body Mode and Flexible Modes Under Relay
Control with Dead-Zone and Hysteresis and Nonlinear Feedback
J
h
θ
2
EI, ρ, L
m
t
, I
t
θ
1
τ
x=0
x=L
111
() () () ()
pd
et k t k t t
2, .1
1, .1, =.02
40, 1, 1
1, 2, 0
pd
ht t
kk
EI L
Jm I
e
u
τ
δ
Δ
Δ
102
The following Figure shows the nonlinear nature of the phase-plane switching surface
using the Time-Optimal control law for a rigid-body slewing maneuver to determine the
timing of relay switches.
dot ()
1
2
dot dot
1 0.5 0 0.5 1
1
0.5
0
0.5
1
dot
dot ()
dot ()
1
2
dot dot .1
1 0.5 0 0.5 1
1
0.5
0
0.5
1
dot
dot
dot () dot ()
Figure 8: Nonlinear Switching Surface for Ideal Relay (left)
and Relay with Dead-Zone (right)
The equations of motion are first put into modal form per (3.2), and then arranged in
state variable form (5.31), with the nonlinear controller included in the output equation as
shown in Figure 5.26. This is done to define the input to the relay as a function of the
system output y(t). This form assumes r = 0 in the feedback system depicted in Figure
2.3, which is true for the limit cycle analysis considered here.
Figure 5.26: Rearranged Block Diagram with System Input r = 0
1
() Ne
()
()
x xuy
yx
AB
e u
y
Nonlinear Controller
and Flexible System
Relay
Nonlinearity
103
()
()
x xuy
yx
AB
(5.31)
if ( )
or ( ) and ( )
() if ( )
or ( ) and ( )
0otherwise
Myt
yt ut M
uy M y t
yt ut M
The system input u(y) represents a relay nonlinearity with dead-zone δ, hysteresis Δ, and
magnitude M, based on the feedback signal y(t). The matrices A and B are defined in
(4.5).
The resulting eigenvalues and normalized eigenvectors used to assemble the A and B
matrices defined in (4.5) for this plant are as follows, where modal damping of 1%
(0.01) has been assumed.
0.5477 0.8252
0, 3.9823,
0.5477 1.0046
n
U (5.32)
01 0 0 0
0 0 0 0 0.5477
,
00 0 1 0
0 0 15.8587 0.0796 0.8252
AB (5.33)
The controller in this example is nonlinear, so the C output matrix in (4.5) is not used.
The gains are applied in the nonlinear output equation ( ) y x in (5.31).
Since the mode shapes of the pinned-free beam with hub and tip mass are derived
using the assumed modes method, components of the modal matrix U in (5.32) are
derived from eigenfunctions of the beam at a location x along the beam
104
() ()
T
j beam beam
j
ux x V Φ (5.34)
where for the two mode system at hand j = 1,2 (rigid-body mode, 1st flexible mode). The
vector ()
beam
x Φ in (5.34) is an 1 N vector of admissible functions (eigenfunctions)
for a nominal pinned-free beam without end masses, where N is the number of nominal
beam modes used to represent the real beam (N = 20 in this example). The vector
j
V is
the 1 N normalized eigenvector for the jth mode of the eigenequation for the
discretized real beam system (see Yang [51] for details)
2
MVKV (5.35)
where [M] is the mass matrix and [K] is the stiffness matrix of the discretized system, and
V is normalized using
T
VMV I . Inertial effects of the hub and end mass on the
real beam are accounted for in the mass matrix [M].
For this example problem torque is applied at the base, and angle at the base is
measured. Therefore, the slope of the beam mode shapes is of interest, not the transverse
displacement, so (5.34) is modified as follows
() ()
T
j beam beam
j
d
ux x V
dx
Φ (5.36)
where ( )
j beam
ux represents the slope of the beam mode shape at location x along the
beam. Therefore, the resulting modal matrix U is constructed as follows
111 12 1
221 22 2
x
uu x
uu x
U
(5.37)
105
where θ
1
and θ
2
are angles at the hub and tip mass respectively. Equation (5.37) can also
be stated as follows
111 12 1
221 22 2
1
mode1 mode2
2
, , sensor_location,mode#
angle sensor at location 1
angle sensor at location 2
( )
ij j beami
uu x
uu x
x
uu
x
uux u
(5.38)
so that the i,jth components of U in (5.32) are evaluated using (5.36) at each sensor
location i for each mode j. The resulting components of U for this example are
, , sensor_location,mode#
11,1
1
12,1
1
21,2
2
() () , , 1,2
for the rigid body mode 1
(0) (0) 0.5477
( ) ( ) 0.5477
for the 1st flexible mode 2
(0) (0) 0
T
jbeami beami
j
T
T
T
d
ux x V u i j
dx
d
uVu
dx
d
uL L V u
dx
d
uVu
dx
Φ
Φ
Φ
Φ
22,2
2
.8252
( ) ( ) 1.0046
T d
uL L V u
dx
Φ
(5.39)
The system in the current example has 4 states in modal coordinate form, with
angular position and angular velocity measured at physical coordinate location 1, the hub
at the base of the beam. The relation between physical coordinates and modal coordinates
is
11 2 2
T
q
(5.40)
106
111 12 1
211 122
321 22 3
421 224
00
00
ˆ
00
00
qu u x
qu ux
x
qu u x
qu ux
U (5.41)
122
() ( ( )) () () ()
pd
ytxt KqtKqtqt (5.42)
where
11 1 1
22 2 2
ˆ
( ) ( ) state transformation to physical coordinates
ˆ
( ) ( ) 1 0 0 0 isolates postion term
ˆ
( ) ( ) 0 1 0 0 isolates velocity term
qt x t
qt xt q
qt xt q
U
CU C
CU C
(5.43)
and
ˆ
U is composed of components from the modal matrix U in (5.32).
Using (5.43) in (5.42), the nonlinear output equation becomes
122
12 2
( ( )) () () ()
ˆˆ ˆ
() () ()
pd
pd
xt Kq t K qt qt
Kxt K xt xt
CU CU CU
(5.44)
5.3.2 Limit Cycle Frequencies and Fixed Points
In order to determine the limit cycle frequency using the state-space equations, it is
necessary to calculate the time at which the relay switches under limit cycling conditions.
It is assumed that the limit cycle is symmetric, with only two switches and two resets per
cycle. These are the characteristics of a unimodal periodic motion, and the waveform and
input/output relations are shown in Figure 5.27.
107
0 0.5 1 1.5 2
2
1
0
1
2
+ δ
– δ
t = t* = T/2
t = 2t* = T
t = 0
t = γt*
y(t) u(t) = 0, ±M
Δ
Figure 5.27: Symmetric Unimodal Limit Cycle Waveform in Relay Feedback System
In order for the symmetric unimodal limit cycle in Figure 5.27 to exist in the system
given in (5.31), the following must hold.
( ) ( ) the motion is periodic with period
( ) ( 2) the motion has odd symmetry
y t ytT T
yt y t T
(5.45)
Equation (5.45) provides constraints on the system motion for establishing the periodic
conditions required for the limit cycle to exist. Using the relations in (5.45), the initial
conditions of the state must satisfy
*
() (0) xtx (5.46)
where t* is the limit cycle half-period shown in Figure 5.27. Since the system is
piecewise linear, we have that between relay switches the state is given by
() ()
() ( ) ( )
o
o
t
tt t
o
t
xte xt e u d
A A
B (5.47)
For a relay with dead-zone the input is constant or zero, resulting in the following
state dynamics
*
**
for 0
for
x xM t t
xxttt
AB
A
(5.48)
108
At relay reset time t = γt* the output y(t) is crossing the upper dead-zone minus hysteresis
boundary () of the relay with negative slope (() 0) yt . At relay switching time
t = t* the output y(t) is crossing the lower dead-zone plus hysteresis boundary(( ))
of the relay with negative slope. The state solution (5.47) at these instants in time and the
symmetry constraint (5.46) can be used to determine the initial conditions that result in a
limit cycle
*
**
**
*()
0
*( ) *
() e (0) ( )
() e ( ) (0)
t
tt
tt
xtx e d M
xt x t x
AA
A
B
(5.49)
*
*
*
1
*
yielding (0) e
t
t
tt
x ed M
AA
IB (5.50)
Also, referring to Figure 5.27, at t = 0 and subsequent switching instant, the relay input
crosses the hysteresis bands as follows:
**
(0) ( (0)) ( )
() (( )) ( )
yx
yt x t
(5.51)
Solutions of γ, and
t* which satisfy (5.51) are obtained numerically, yielding limit cycle
frequencies of ω* = π/t*. Note that an initial guess is required for the parameter γ. For the
assumed symmetric unimodal limit cycle waveform, this parameter obeys the following
constraint
01 (5.52)
The following additional waveform criterion must be satisfied in order for the solution to
be a limit cycle (see Figure 5.27).
109
*
To ensure proper velocity direction at switches
( ) (0) 0, ( ) ( ) 0 iy ii y t
(5.53)
*
**
To ensure proper waveform between switches
() ( ) ( ), 0
()() () (),
iyt tt
ii y t t t t
(5.54)
Limit cycle stability is verified by the Jacobian developed in Section 5.3.4 for a relay
with dead-zone and hysteresis and a nonlinear switching surface.
The Describing Function Method (DFM) can often provide a preliminary estimate of
the frequency range of interest for existing limit cycles in linear systems. However, since
the controller in this case is nonlinear, the DFM cannot be used.
Note that when the relay contains a dead-zone, the solution to the limit cycle half-
period t* must satisfy the two constraints shown in (5.51). This requires numerical
methods to solve the two simultaneous equations, and an initial guess for both parameters
γ and t*. If it is initially assumed that the relay does not contain a dead-zone, only one
equation is needed to solve for t*. Note that when there is no dead-zone present γ = 1, so
γt* = t* (see Figure 5.27). Figure 5.28 shows the solution locus generated by rearranging
the first constraint equation in (5.51) as ((0)) 0 x and varying t*, where solutions
occur at intersections with the x axis. These solutions are then verified with the stability
and waveform tests. This preliminary analysis provides a time span for the half-period
time sweep of t*, used in the numerical analysis that finds both γ and t* from (5.51).
110
Figure 5.28: Solutions of (5.51) with 0 to Estimate Limit Cycle Period Range
The results of the time sweep analysis when the dead-zone is included are shown in Table
5.2. There are 3 stable limit cycle periods of Tp = 4.09, 5.03, and 6.39 seconds. Note that
these periods are relatively close to the three larger solutions for the zero dead-zone
reference case shown in Figure 5.28. The high frequency solution does not exist for a
relay with the given dead-zone. Table 5.2 shows the relay reset time parameters for each
limit cycle. These parameters were solved for by sweeping the half-period time parameter
t* from 0 to 5 seconds, and using the Matlab command fsolve on (5.51) at each trial time
for t*, along with an initial guess for γ. The initial guess used was γ = 0.8. The range for
the time sweep was arrived at using the range of t* solutions for the zero dead-zone case
shown in Figure 5.28. The resulting solutions for t* and the reset time parameter γ are
then used in (5.49) to determine the fixed points on each switching surface. The
t* = 0.6711, 2.097, 2.4546, 2.8230 sec
T
p
= 1.3422, 4.0594, 4.9093, 5.6460 sec (T
p
= 2t*)
4 Stable LC periods using eigenvalue test
Solutions to
((0)) 0 x
cross the x axis
111
waveform criteria (5.54) and stability criterion (5.67) for each limit cycle are then
verified to determine the stable fixed points. The fixed points for each limit cycle are
x
1
=[0.0951, 0.4429, 0.0012, 0.1566]
T
, x
2
=[0.1007, 0.5967, -0.0052, 0.1216]
T
, and
x
3
=[0.0979, 0.8090, -0.0043, 0.0038]
T
. Figure 5.29 shows the limit cycle waveforms
started from fixed point initial conditions.
Table 5.2: Relay Reset Time Parameters (Figure 5.27)
LC γ
t*(sec)
Tp(sec)
ω(rad/s)
1 0.7901 2.0469 4.0939 1.5348
2 0.8659 2.5163 5.0325 1.2485
3 0.9243 3.1959 6.3917 0.9830
Figure 5.29: Three Stable Limit Cycle Solutions Started from Fixed Point Initial
Conditions, System Output y(t) (blue), Relay Output u(y) (red)
112
5.3.3 Limit Cycle Fixed Point Region of Attraction Mapping
The switching surface constraint (0) ( (0)) ( ) yx in (5.51) removes one degree of
freedom from the system, by making one state a function of the other remaining states.
This means that only n-1 states need to be included in the cell space, since the remaining
state is no longer independent. If the output is a linear combination of the states, such as a
PD controller, there is freedom to choose any state for elimination. However, due to the
nonlinear nature of the output equation in this example, one of the position states (x
3
) is
chosen as the dependent state, since the velocity states (x
2
and x
4
) are related by the
nonlinear term, and cannot be isolated. Thus the 3D cell space defines the x
1
, x
2
, and x
4
coordinates of the state-space initial condition points, while the x
3
coordinate for these
points is defined by the x
1
, x
2
, and x
4
coordinates, and the switching constraint equation.
Figure 5.30 shows the limit cycle trajectories and a 2D slice of the 3D switching
surface cell map in x1-x2-x4 space. The 3D cell map is centered at the cluster of limit
cycle fixed points found in Section 5.3.2. The 2D slice is for visualization purposes, and
is defined by a plane that passes through the 3 limit cycle fixed points. Figure 5.31 shows
the x2-x4 projection of the map in Figure 5.30, as well as a test point for verification
purposes. The test point is located in the ROA of LC2 (green Xs), and was chosen from a
small region of attraction area on the map, in order to verify the map’s accuracy at
predicting the resulting limit cycle.
Figure 5.32 shows a simulation of the system starting from the test point, and verifies
that LC2 is excited from this location on the switching surface.
113
Figure 5.30: ROA Map at Slice Through 3 Fixed Point Locations and Limit Cycle State
Trajectories for LC1=blueX, LC2=greenX, LC3 =redX
Figure 5.31: x2-x4 View of 2D Slice in Figure 5.30. ROA for LC1=blueX, LC2=greenX,
LC3 =redX
Initial state at x = [0.0945 0.4005 0.0020 0.1724]
T
will result in LC2 (green Xs), which was verified
by simulation (see Figure 5.32).
LC3
LC2
LC1
Full switching surface map is 3D (n-1). Slice is
defined by plane that contains all 3 fixed points.
114
Figure 5.32: Simulation Showing LC2 Excited from Initial State Noted on ROA Map in
Figure 5.31
As noted above in Section 5.2.3, if the cell map is constructed with initial conditions
that do not lie on the switching surface, the dimension of the cell space is no longer less
than the dimension of the state space. However, the increase in total number of cells can
be reduced, if only specific modal initial conditions are investigated. In the following
case, only the rigid-body mode initial conditions are mapped, in order to investigate
which limit cycle will be excited with various combinations of rigid-body initial
conditions of displacement and velocity. Note that for this case both the hub and the tip
mass are initially assumed to have the same angular position and velocity. All other
modal initial conditions are set to zero (x(0) = [x
1
(0) x
2
(0) 0 0]
T
). Thus the cell map of
initial conditions is only 2D. This would be the case if the controller is turned on when
the flexible modes are not initially excited (no initial beam deflection or relative angular
velocity between the hub and tip mass). Note that all modes are still included in the
subsequent dynamics.
115
The following Figure 5.33 shows a phase plane representation of the regions of
attraction for the limit cycles, when initiated from rigid-body mode initial conditions.
One initial condition in each ROA is selected, and simulation results are shown in Figure
5.34. Compare these results to those in Figure 5.29.
Figure 5.33: Rigid-Body Mode Region of Attraction Map: ROA for LC1=blueX,
LC2=greenX, LC3=redX
Black Xs represent maps to Sink Cell (zero velocity ICs
started within Relay dead-zone resulting in no motion).
Initial state at a) x = [0 0.3824 0 0]
T
(LC1);
b) [0 0.5294 0 0]
T
(LC2); c) [0 0.6471 0 0]
T
(LC3),
verified by simulation (see Figure 5.34).
116
c)
LC3
Tp3=6.39 sec
Amp3=.755
b)
LC2
Tp2=5.03 sec
Amp2=.472
a)
LC1
Tp1=4.09 sec
Amp1=.301
Figure 5.34: Simulation Showing Limit Cycles Excited from Rigid-Body Initial
Conditions Noted on ROA Map in Figure 5.33
117
5.3.4 Stability Jacobian for Nonlinear Switching Surface
The following shows how to calculate the Jacobian for a nonlinear switching surface
using (4.28) repeated below. The nonlinearity arises from a signed velocity squared term
in the system output equation
111
(())
pd
yxt k k
.
*
1
*
1
1
(, )
(, )
i
s
xi i
ii
i
i
ti i
x
v
x
W Φ (5.55)
where
i
i
e
A
Φ and
*
ii i
vx d AB .
The system in the current example has 4 states in modal coordinate form, with
angular position and angular velocity measured at physical coordinate location 1, the hub
at the base of the beam. The relation between physical coordinates and modal coordinates
is
11 2 2
T
q
(5.56)
10 11 11
10 11 22
20 21 33
20 22 44
00
00
ˆ
00
00
uu qx
uu qx
x
uu qx
uu qx
U (5.57)
122
() ( ( )) () () ()
pd
ytxt KqtKqtqt (5.58)
where
11 1 1
22 2 2
ˆ
( ) ( ) state transformation to physical coordinates
ˆ
( ) ( ) 1 0 0 0 isolates postion term
ˆ
( ) ( ) 0 1 0 0 isolates velocity term
qt x t
qt xt q
qt xt q
U
CU C
CU C
(5.59)
118
and
ˆ
U is composed of components from the modal matrix U.
Using (5.59) in (5.58), the nonlinear output equation becomes
122
12 2
(()) ( , ) () () ()
ˆˆ ˆ
() () ()
pd
pd
xt xt Kq t K qt qt
K x t K xt xt
CU CU CU
(5.60)
Recall the piecewise linear state response equations from above. The state equation
between switching surfaces i and i-1 and its solution on switching surface i are given by
ii
i
x xd AB (5.61)
() *
1
0
() e
i
ii
i
ii i
x xe d d
AA
B (5.62)
1,2, , , number of switching surfaces
relay magnitude between switching surfaces and 1
time duration between switching instance and 1
i
i
iss
dii
ii
Next calculate the derivative terms in (5.55)
Derivatives with respect to x –
122
12
122
12
1
11
1
2
22 2 2
2
(, ) () () ()
(, ) () () ()
1 ( ,) 2( ,)
1( , ) () ()
2 ( , ) () () 2 ( ) ( )
pd
xp d
xx
xp p
xd d
xt Kq t K qt qt
qq
xt K q t K qt qt
qx q x
xt xt
q
xt K q t K q t
qx x
q
xt K qt qt K q t q t
qx x
12 2
( , ) ( ) 2 () ()
xp d
xtK qt Kqt qt
xx
(5.63)
119
For trajectory started from fixed point
*
1 i
x
on originating switching surface i-1 and
switch on subsequent switching surface i at time
i
t
*
1
*
1
*()
11 1
0
,
1
*
12 2
()
22
0
,
2
ˆ
1( , ) ()
ˆ
2( , ) 2 () ()
ˆ
2()
2()
ii
i
ii
t
ii tt
xi i p p i
xx t
p
ii i
xi i d
t
i tt
di
xx t
i
di
xK qt K exe dd
xx
Ke
x Kqt qt
x
Kq t e x e d d
x
Kq
AA
A
AA
CU B
CU
CU B
2
**
22 2 2 2
ˆ
ˆˆ ˆ ˆ
2using () ()
i
i
ii
di i i i
e
Kx e q x x
A
A
CU
C U CU CU C U
**
11 2 2
*
12 2
ˆˆˆ
(, ) 2
ˆˆ ˆ
2,
ii
i
i
xi i p d i
pd i i i
xK e K x e
KK x e
AA
A
CU CU CU
CU C U C UΦΦ
(5.64)
Derivatives with respect to t –
Using a similar derivation to (5.63)
12 2
122
(, ) () 2 () ()
() 2 ( ) ( )
tp d
pd
xtK qt Kqt qt
tt
Kq t K qt qt
(5.65)
For trajectory started from fixed point
*
1 i
x
on originating switching surface i-1 and
switch on subsequent switching surface i at time
i
t
122
** *
11 1
( , ) ( ) 2 () ()
( ,)1( ,)2( ,)
ii ii
tp d
ii i
ti i t i i ti i
xtKq t Kq tq t
xx x
120
where
*
1
*
1
*
11 1
,
*
1
**
*
12 2
,
22
ˆ
1( ,) () ( )
ˆ
using velocity just before switch
() , ( )
2( , ) 2 () ()
ˆˆ
2()
ii
ii
ii i
ti i p p i
xx t
pi i
ii
ii i i i
iii
ti i d
xx t
ii
di
xKqt K x
Kx d
xxdx x
x K qt qt
KC x C x
CU
CU A B
AB
UU
**
22
()
ˆˆ
2 ( ( ), ( ) same as above)
i
ii
di i i i i
Kx x d x x
CU C U A B
**
11 22
122
*
ˆˆˆ
(, ) 2
ˆˆˆ
2
using ( ) (velocity just before switch)
i
ti i p d i i
pd i
i
ii i i
xK K b x d
KK b v
xd v x
CU C U C U A B
CU C U C U
AB
(5.66)
Substituting (5.64) and (5.66) into (5.55) yields
21 **
12 0 1
22 1 1
21 **
12 0 1
*
12222
22
*
12222
*
12121
11
*
121
(, ) ( , )
(, ) ( , )
ˆˆ ˆ
2
ˆˆ ˆ
2
ˆˆ ˆ
2
ˆˆ
2
xx
tt
pd
pd
pd
pd
xx
vv
xx
KK x
v
KK x v
KK x
v
KK x
WΦΦ
CU CU CU Φ
Φ
CU CU CU
CU CU CU Φ
Φ
CU C U
21
ˆ
v
CU
*
21 202
2
*
12022
*
11 212
1
*
12121
ˆˆˆ
2
ˆˆˆ
2
ˆˆ ˆ
2
ˆˆ ˆ
2
pd
pd
pd
pd
vK K x
KK x v
vK K x
KK x v
CU CU CU
WI Φ
CU CU CU
CU C U C U
I Φ
CU C U C U
(5.67)
121
where
*
12 1
* *
0 0
* *
1 1
() *
12 1
initial condition fixed point on switching surface ( ) ( )
relay reset fixed point on switching surface ( ) ( )
, , relay reset time, limit cycle half
t
xx
xx
ee e t
AA A
ΦΦ
12
*
11 1
* *
20
*
2
-period
, 0 relay magnitude between switching surfaces
state velocity just before reset from to 0 at ( )
( ) state velocity just before switch from 0 to at ( )
dMd
vx M M x
vx Mxt
x
AB
A
*
0
limit cycle symmetry constraint x
5.3.5 Discussion of Results
This example extends the previous example to a discretized beam model, using the
assumed modes method to construct the plant model. Eigenfunctions are now used to
determine components of the modal matrix. Additionally, the control law is nonlinear,
resulting in a nonlinear switching surface. As a result of the nonlinear switching surface,
a new stability Jacobian is developed from the general Jacobian in Section 4.2.5.
For cell mapping purposes, the dependent state variable defined by the nonlinear
switching surface had to be chosen with more care in this case, since the output velocity
term in the control law is a squared term. The output velocity is made up of a
combination of two modal states x
2
and x
4
, so neither of these could be separated and
made an explicit function of the remaining states. Therefore x
3
was chosen as the
dependent state, and x
1
, x
2
, and x
4
were the independent states used for mapping purposes.
A 2D slice was taken at the center of the fixed point cluster, and a limit cycle was
correctly predicted from a small ROA. Additionally, a 2D map of rigid-body mode initial
conditions, was verified by simulation using 3 test points from 3 different ROAs.
122
5.4 Nonlinear Plant Under Relay Feedback Control
The following example demonstrates how to apply the method to a nonlinear plant. The
nonlinearity in this example is due to nonlinear stiffness. The flexible mode is driven by a
nonlinear spring modeled as a softening spring, which causes instability to occur when
the relative deflection between the masses reaches a certain point. This instability is
caused by a reversing of the spring restoring force as shown in Figure 5.36. The
linearized version of the plant ( ε = 0) has the same characteristics as the example in
Section 5.1.
5.4.1 Plant Description
Figure 5.35 shows an unconstrained 2 mass-spring-damper plant with a nonlinear
softening spring and one rigid-body and one flexible mode. Control force is applied to
mass 1 according to (5.69).
Figure 5.35: 4-State System with 1 Rigid-Body Mode and 1 Nonlinear Flexible Mode
Under Relay Control with Hysteresis and PD Feedback
12
1, 2, 10, 0.0516, 3.5 mm k c
11
()
pd
ekq kq
Collocated PD
feedback at m1
1, 0.1
pd
kk
Force on m1
1, 0.1 D
e
u
PD Controller with On-Off Actuator with
hysteresis
eff
k
1
() qt
2
() qt
c
2
21
1( )
eff
kk q q
Effective spring rate
123
Figure 5.36: Nonlinear Softening Spring Force Characteristics
(k = 10, ε = –3.5, z = q
2
– q
1
)
Defining the state variables as
11 2 2
T
x qqq q , the equations of motion in state
variable form are
12
2
23131 42
111
34
2
43131 42
22
12
1
1( ) ( ) ( ) ()
1( ) ( ) ( )
linear spring rate
nonlinear spring coefficient
proportional damping coefficient
()
pd
xx
kc
x xx xx x x uy
mmm
xx
kc
xxxxx xx
mm
ykx kx
k
c
uy
relay forcing term
(5.68)
System becomes unstable
when spring force reverses
direction at
1
1
(3.5)
0.5345
z
stable region
Linear spring
Nonlinear spring
21
zq q
124
The resulting state equations in vector form are
1
() ( )
if or and ( )
()
if or and ( )
01 0 0 , 0 0
T
pd
x fx u y
yx
M yt y ut M
uy
Myt y ut M
Bm Ckk
B
C
(5.69)
The system input u(y) represents a relay nonlinearity with hysteresis Δ, and magnitude M,
based on mass 1 position and velocity feedback signal
11 pd
ye kq kq .
To compare results with the linear plant in Section 5.1, the linear plant modal
damping of 1% is simulated in the nonlinear plant using the following damping
coefficient
12
12
22
neff n
mm
cm
mm
(5.70)
where ζ is the modal damping ratio, ω
n
is the flexible mode frequency of the linear plant,
and m
eff
is the effective mass of the system with regards to natural frequency and the
linear spring stiffness
neff
km .
5.4.2 Limit Cycle Frequencies and Fixed Points
For a nonlinear plant, the limit cycle frequencies and fixed points cannot be determined
ahead of time, since a closed form solution of the state is not available. However, the
mapping step of the SSCM method can accurately predict both the limit cycle frequencies
and their fixed points through the mapping process. Accuracy is achieved through the
integration stopping criterion (4.36). The state trajectory is numerically integrated
125
forward from each initial condition in the cell map, until successive crossings of the
switching surface are crossed within this error. This locates the limit cycle periodic fixed
point on the switching surface. The limit cycle period is the time duration between the
last two switches.
For reference, the limit cycle frequencies and corresponding fixed points from the
linear plant example in Section 5.1 are repeated in Table 5.3, and the limit cycle
waveforms are repeated in Figure 5.37. The natural frequency for the linear plant is
3.8730
n
rad/sec.
Table 5.3: Linear Plant Limit Cycle Periods and Fixed Points (Section 5.1.2)
LC1
LC2 LC3
LC4
(unstable)
LC5
LC
(rad/sec)
0.5074
0.6132
0.8344 1.3932 4.4959
LC
T (sec)
12.3830
10.2462
7.5300 4.5099 1.3975
*
x
0.0000
1.7873
0.0095
-0.1345
0.0000
1.4789
0.0066
0.1125
0.0000
1.0869
0.0129
0.3276
0.0000
0.6510
0.0189
0.5756
0.0000
0.2017
0.0133
0.9493
126
Figure 5.37: Linear Plant Limit Cycle Waveforms (Section 5.1.2)
LC4 not shown because it is unstable (see Table 5.3)
System Output y(t) (blue), Relay Output u(y) (red)
These results will be compared to the limit cycle frequency and fixed point results for the
nonlinear plant, obtained from the cell mapping process in the next section. The nonlinear
plant and relay parameters are the same as in the linear plant example, except for the
additional nonlinear spring coefficient ε. The damping coefficient c has been tuned to
mimic the modal damping of 1% assigned in the linear plant example (see (5.70)).
LC5
127
5.4.3 Limit Cycle Fixed Point Region of Attraction Mapping
The limit cycle periods and fixed points for the nonlinear plant are shown in Table 5.4.
These were found from average period results for cells that were mapped to the same
destination cell. For improved accuracy of the fixed point locations, the actual end state
of the trajectory is reported below, not the center of the destination cell. The waveform
plots in Figure 5.38 show that even with a relatively coarse cell density (31
3
cells),
accurate results for the fixed points can be obtained. This is due to the integration
stopping criterion used in the SSCM method.
Table 5.4: Nonlinear Plant Limit Cycle Periods and Fixed Points
LC1
LC2 LC3
LC4
LC5
LC
(rad/sec)
0.4951
0.5997
0.8025 1.2961 4.1973
LC
T (sec)
12.6901
10.4770
7.8291 4.8478 1.4969
*
x
0.0117
0.8827
-0.0057
1.1448
0.0057
0.9433
-0.0019
0.8386
0.0099
0.9006
-0.0041
0.5288
0.0147
0.8531
-0.0051
0.1789
0.0108
0.8922
-0.0023
-0.2575
Figure 5.38 shows the limit cycle waveforms, when started from fixed point initial
conditions in Table 5.4. Note that LC4 is an unstable limit cycle waveform for the linear
plant model. However, the calculated period of 4.51 seconds for the unstable LC4 shown
in Table 5.3 is close to the 4.85 seconds for the stable LC4 of the nonlinear plant found
by the SSCM procedure. A simulation of the unstable LC4 for the linear plant showed the
same waveform of the nonlinear LC4, however it eventually diverges to the high
128
frequency LC5. The period and amplitude of the nonlinear plant limit cycles are
comparable to those of the linear plant, but the periods are slightly longer (cf. Figure
5.37). This is due to the lower restoring force exerted by the spring for comparable spring
deflections (see Figure 5.36).
Figure 5.38: Nonlinear Plant Limit Cycle Waveforms
System Output y(t) (blue), Relay Output u(y) (red) (relay output scaled for clarity)
The limit cycle waveforms in Figure 5.38 verify that the frequencies and periodic
fixed points in Table 5.4, calculated using the SSCM cell mapping procedure, accurately
predict limit cycle performance for the nonlinear plant model.
129
Next a comparison of 2D 31cell x 31cell ROA maps for both the linear and nonlinear
plants is done. Figure 5.39 shows the ROA for the nonlinear plant from rigid-body initial
conditions
10 30 20 40
(, ) x xx x . The same Region of Interest is used in the ROA map for
the linear plant shown in Figure 5.40. Note that LC4 does not appear in Figure 5.40, since
it is an unstable limit cycle for the linear plant. A simulation of mass1 position q
1
from
the noted points is shown in Figure 5.41, in order to verify the accuracy of each ROA.
Figure 5.39: Nonlinear Plant Rigid-Body Mode Region of Attraction Map: ROA for
LC1=red, LC2=cyan, LC3=magenta, LC4=green, LC5=blue, Divergence=black
10 10 10 10 10
) 0.0, ) 0.5, ) 1.0, ) 1.25, ) 1.75
LC5 LC3 Divergence LC2 LC1
ax b x c x d x e x
Limit Cycle ROA Map for Nonlinear Plant Starting From Rigid-
Body Mode Initial Conditions
10 30 20 40
(, ) x xx x
130
Figure 5.40: Linear Plant Rigid-Body Mode Region of Attraction Map: ROA for
LC1=red, LC2=cyan, LC3=magenta, LC5=blue
Note that the ROA for each limit cycle for the linear plant shown in Figure 5.40 are
very similar to those for the nonlinear plant shown in Figure 5.39. The exceptions are that
LC4 does not exist in the linear plant map, since it is an unstable limit cycle for the linear
plant case. Also, there are no divergent cells in the linear plant map, since the spring
restoring force does not reverse direction, as with the nonlinear spring.
Limit Cycle ROA Map for Linear Plant Starting From Rigid-Body
Mode Initial Conditions
131
Figure 5.41: Nonlinear Plant Limit Cycles Excited by Initial Conditions in Figure 5.39
Mass1 Position q
1
(t) (blue), Relay Output u(y) (red) (relay output scaled for clarity)
a)
LC5
Tp5=1.50 sec
Amp5=0.190
b)
LC3
Tp3=7.83 sec
Amp3=0.673
c)
Divergence (max
spring deflection
exceeded)
d)
LC2
Tp2=10.48 sec
Amp2=1.074
LC1
Tp1=12.69 sec
Amp1=1.696
e)
132
5.4.4 Discussion of Results
This example demonstrated how the SSCM can be applied to a nonlinear plant with a
nonlinear softening spring. Since a closed form of the state solution is not available in
general for nonlinear systems, the limit cycles and their fixed points cannot be
determined before the mapping procedure is implemented. The advantage of knowing
these beforehand is that the cell map boundaries can be selected to enclose the fixed
points, thus reducing the cell map size. Also, the number of expected limit cycles and
their periods can be used for grouping cells with like attributes, during the mapping
process.
Although the limit cycle frequencies and fixed points cannot be determined ahead of
time, they can be found as a byproduct of the mapping process. This is due to the SSCM
error convergence criterion, which stops integrating only after successive crossings of the
switching surface in the proper direction are close. This defines the fixed point for a limit
cycle with the prescribed waveform. As shown in the example, these values can be found
accurately, with a relatively coarse cell map (31
3
in this case).
133
Chapter 6: Concluding Remarks
In this dissertation a new numerical-analytical method allowing the analysis of multiple
limit cycles and their regions of attraction was formulated. This method, called Switching
Surface Cell Mapping (SSCM), provides a tool for analyzing limit cycles occurring in
flexible structures under relay feedback control. When flexible modes are considered in
the system model, multiple limit cycles may exist, compared to the rigid-body case. The
method allows the study of structural vibrations with multiple flexible modes, and
includes the ability to handle nonlinear relay switching surfaces. Local stability analysis
for limit cycles in relay controlled systems in the literature is limited to linear switching
surfaces. The method modifies the cell mapping process by introducing a convergence
criterion that allows for faster convergence.
The following summarizes the method’s capabilities in analyzing limit cycle
characteristics for design purposes, and mapping the regions of attraction of multiple
limit cycles to study the effects of adjacent limit cycles.
Analysis of limit cycle characteristics
a) allows computation of all frequency, amplitude, and local stability characteristics
b) multiple flexible modes included in plant model
c) it can handle linear or nonlinear switching surfaces
d) multiple relay switching surfaces
Computation of limit cycle regions of attraction
a) multiple fixed point regions of attraction
b) use of Poincaré mapping surface increases memory efficiency
134
c) error convergence criterion increases accuracy
d) nonlinear plant capability
For future research two main areas come to mind. One area is extending the method
to limit cycles with more complex switching patterns. The focus of this dissertation was
on symmetric unimodal limit cycles. More complex switching patterns could possibly be
analyzed, if more switching surfaces were used in the mapping process. However, this
would increase the cell-space dimensions, and the focus here was to reduce the cell space,
so that more flexible modes could be included in the model. A second area would be to
apply the method to systems with more complex control laws. Since controller design
was not the focus of this research, a linear PD controller was assumed throughout most of
this paper, with the exception of a nonlinear example that included a signed velocity
squared term.
135
References
[1] Agrawal, B. N., McClelland, R., and Song, G. (1997) “Attitude Control of Flexible
Spacecraft Using Pulse-Width Pulse Frequency Modulated Thrusters”, Journal
of Industrial and Commercial Application in Space, Vol. 17, No. 1, 1997, pp
15-34.
[2] Amrani, D., and Atherton, D. (1989) “Designing Autonomous Relay Systems with
Chaotic Motion”, Proc of the 28
th
Conference on Decision and Control, Tampa,
FL December 1989, 512-517
[3] Asada, H., and Slotine, J., (1986) Robot Analysis and Control, John Wiley and
Sons, New York
[4] Astrom, K. J., (1995) “Oscillations in Systems with Relay Feedback”, The IMA
Volumes in Mathematics and its Applications: Adaptive Control, Filtering, and
Signal Processing, Vol. 74:1-25, Springer-Verlag, New York.
[5] Atherton, D. P. (1975) Nonlinear Control Engineering, Van Nostrand Reinhold Co
Ltd, New York.
[6] Atherton, D. P. (2000) “Relay Autotuning: A Use of Old Ideas in a New Setting”,
Trans of the Inst of Meas and Cont, Vol. 22, No. 1, 2000 pp 103-122
[7] Balas, M. J. (1978) “Feedback Control of Flexible Systems”, IEEE Trans on
Automatic Control, Vol. AC-23, No. 4, August 1978, 673–679.
[8] Balas, M. J. (1982) “Trends in Large Space Structure Control Theory: Fondest
Hopes, Wildest Dreams”, IEEE Trans on Automatic Control, Vol. AC-27, No.
3, June 1982, 522–535.
[9] Boiko, I. (1999) “Input-output Analysis of Limit Cycling Relay Feedback Control
Systems”, Proc of the American Control Council, June 1999 542–546.
[10] Boiko, I. (2006) “Refinement of Periodic Solution Obtained via Describing
Function Method”, Proc of the American Control Council, June 2006 1517–
1522.
[11] Boiko, I. (2009) Discontinuous Control Systems; Frequency-Domain Analysis and
Design, Birkhäuser, Boston.
[12] Bryson, A. E., (1994) Control of Spacecraft and Aircraft, Princeton University
Press, NJ.
136
[13] Corless, M. (1993) “Control of Uncertain Nonlinear Systems”, ASME Journal of
Dynamic Systems and Control, Vol. 115, June 1993, 362-372.
[14] Corless, M. and G. Leitmann (1993) “Bounded Controller for Robust Exponential
Convergence,” J. Optim. Theory Appl., vol. 76, 1993, 1-12,
[15] Ding, Q., Cooper, J. and Leung, A. (2005) “Application of an Improved Cell
Mapping Method to Bilinear Stiffness Aeroelastic Systems”, Journal of Fluids
and Structures 20 35-49
[16] Edwards, C. and Spurgeon, S. (1999) Sliding mode control: theory and
applications, Taylor and Francis, London.
[17] Flashner, H., and Hsu, C. S. (1983) “A Study of Nonlinear Periodic Systems Via
the Point Mapping Method”, International Journal for Numerical Methods in
Engineering, Vol. 19, 185-215
[18] Flashner, H., and Guttalu, R. S. (1988) “A Computational Approach for Studying
Domains of Attraction for Non-linear Systems”, Int. J. Non-linear Mechanics,
Vol. 23, No. 4, 1988, 279–295.
[19] Gelb, A. and Vander Velde, W. (1968) Multiple-Input Describing Functions and
Nonlinear System Design, McGraw-Hill, New York.
[20] Genesio, R., Tartaglia, M., and Vicino, A. (1985) “On the Estimation of Asmptotic
Stability Regions: State of the Art and New Proposals”, IEEE Trans on
Automatic Control, Vol. A-C30, No. 8, August, 1985, 747–755.
[21] Golat, M. C. (2002) “A Global Study of Nonlinear Dynamical Systems by a
Combined Numerical-Analytical Approach”, PhD Thesis, Department of
Mechanical Engineering, USC.
[22] Goncalves, J. M. (2000) “Constructive Global Analysis of Hybrid Systems”, PhD
Thesis, Department of Electrical Engineering and Computer Science, MIT.
[23] Goncalves, J. M. et al. (2001) “Global Stability of Relay Feedback Systems”, IEEE
[24] Goncalves, J. M. (2005) “Regions of Stability for Limit Cycle Oscillations in
Piecewise Linear Systems”, IEEE Trans on Automatic Control, Vol. 50, No. 11,
November 2005, 1877–1882.
[25] Graham, D. and Duane, M. (1961) Analysis of Nonlinear Control Systems, Dover
Publications, New York.
[26] Henon, M. (1982) “On the Numerical Computation of Poincaré Maps”, Physica 5D,
North-Holland Publishing, 412-414.
137
[27] Hsu, C. S., and Guttalu, R. S. (1980) “An Unraveling Algorithm for Global
Analysis of Dynamical Systems: An Application of Cell-to-Cell Mappings”,
Trans of the ASME, Vol. 47, December 1980, 940–948.
[28] Hsu, C. S. (1987) Cell-to Cell Mapping: A Method of Global Analysis for Nonlinear
Systems, Applied Mathematical Sciences 64, Springer-Verlag, New York.
[29] Johansson, K., Barabanov, A., and Astrom, K. (2002) “Limit Cycles With
Chattering in Relay Feedback Systems”, IEEE Transactions on Automatic
Control, Vol. 47, No. 9, Sept 2002, 1414-1423
[30] Junkins, L. J., and Kim, Y., (1993) Introduction to Dynamics and Control of
Flexible Systems, AIAA Education Series, American Institute of Aeronautics
and Astronautics, Inc., Washington D.C.
[31] Khalil, H. (2002) Nonlinear Systems, 3
rd
ed, Prentice Hall, New Jersey.
[32] Li, X., DeCarlo, R., and Corless, M. (2003) “Sliding Manifold Design for a Class of
Uncertain Time Delay Systems”, Proc of the American Control Conference,
Denver, CO, June 4-6, 2003, 1518-1523
[33] Likins, P. W., (1970) “Dynamics and Control of Flexible Space Vehicles”, TR 32-
1329, Rev. 1, 1970, Jet Propulsion Lab., Pasadena, CA.
[34] Levitas, J., Weller, T., and Singer, J. (1994) “Poincaré-Like Simple Cell Mapping
For Non-linear Dynamical Systems”, Journal of Sound and Vibration, 176(5),
641-662
[35] Levitas, J., and Weller, T. (1995) “Poincaré Linear Interpolated Cell Mapping:
Method for Global Analysis of Oscillating Systems”, Journal of Applied
Mechanics, Vol. 62, June 1995, pp 489-495.
[36] Moeini, A. (1996) “Computer Aided Design of Nonlinear Phenomenon in Relay
Control Systems”, PhD Thesis, School of Engineering, University of Sussex.
[37] Nikolay, N., Parnferov, N., and Boiko, I. (2001) “Input-output Analysis of Dead-
Zone Relay Control Systems”, Proc of the American Control Conference,
Arlington, VA, June 25-27, 2001, 1507-1512
[38] Shorten, R., Corless, M. et al. (2009) “A Quadratic Stability Result for Singular
Switched Systems with Application to Anti-windup Control”, 2009 American
Control Conference, St. Louis, MO, June 10-12, 2009, 1917-1922
[39] Sidi, M. J. (2002) Spacecraft Dynamics & Control, University Press, Cambridge.
138
[40] Slotine, J. E., and Li, W. (1991) Applied Nonlinear Control, Prentice Hall, New
Jersey.
[41] Song G. and Agrawal B. (2001) “Vibration Suppression of Flexible Spacecraft
During Attitude Control”, Acta Astronautics, Vol. 49, No. 2, 2001, pp 73-83
[42] Spong, M., and Vidyasagar, M., (1989) Robot Dynamics and Control, John Wiley
and Sons, New York
[43] Stacey, A. and Stonier, R., (1998) “Analytic Estimates for the Boundary of the
Region of Asymptotic Attraction”, Dynamics and Control, 8, 177-189.
[44] Tongue, B. H., and Gue, K., (1988) “Interpolated Cell Mapping of Dynamical
Systems”, Journal of Applied Mechanics, Vol. 55, June 1988, pp 461-468.
[45] Tongue, B. H., and White, M. T., (1995) “Application of Interpolated Cell Mapping
to an Analysis of the Lorenz Equations”, Journal of Sound and Vibration, Vol.
188(2), 1995, pp 209-226.
[46] Tsiotras, P., Corless, M., and Longuski, J. (1993) “Invariant Manifold Techniques
for Attitude Control of Symmetric Spacecraft”, Proc of the 32
nd
Conference on
Decision and Control, San Antonio, TX, December 1993.
[47] Tsypkin, Y. Z. (1984) “Relay Control Systems”, University Press, Cambridge.
[48] Utkin, V., Guldner, J., and Shi, J. (2009) Sliding Mode Control in Electro-
Mechanical Systems, 2
nd
ed, CRC Press, New York
[49] Wie, B. and Plescia, C. T. (1984) “Attitude Stabilization of Flexible Spacecraft
During Stationkeeping Maneuvers”, Journal of Guidance, Control, and
Dynamics, 1984, Vol. 7, No. 4, 430-436
[50] Wie, B. (1998) Space Vehicle Dynamics and Control, AIAA Education Series, VA.
[51] Yang, B. (2005) Stress, Strain, and Structural Dynamics, Elsevier Academic Press
Inc, Boston.
[52] Zak, S., Brehove, J., and Corless, M. (1989) “Control of Uncertain Systems with
Unmodeled Actuator and Sensor Dynamics and Incomplete State Information”,
IEEE Trans on Systems, Man, and Cybernetics, Vol. 19, No. 2, March 1989
139
Appendix A: Block-Diagonal Forms of Switching Constraint
Equation Matrices
The following Appendix shows the block-diagonal form of the constraint equation
matrices found in Section 3.3 that were used in the Matlab program developed for this
research. This form allows accurate determination of limit cycle frequency candidates,
since it eliminates the need for numerical techniques for matrix inversion and solution of
the exponential matrix. This increased accuracy is important for evaluating limit cycles in
flexible systems that contain many lightly damped modes.
The following three equations show the context for the matrices derived in this
appendix.
Output equation: Find limit cycle half-period solutions t = t*
*
*
1
0
(0) e
t
t
yedM
AA
CI B (A.1)
Limit cycle fixed points: Solutions that satisfy (A.1) plus additional waveform criterion
(3.27) yield limit cycle fixed points x*
1
*
(0) x xI M
ΦΓ (A.2)
Limit cycle stability test: Eigenvalues of W are inside the unit circle
v
v
C
WI Φ
C
(A.3)
The form of
1
0
, () , and ()
t
tt
et ed t e I
AA A
Γ B Ψ used in the above equations
depends on the type of controller used. For a PI or PID controller an additional state is
140
needed to account for the integral of position term. Therefore, the derivations are divided
into two sections one P or PD controllers and the second for PI or PID controllers.
P or PD controller
00
11
01
,
nn
n
AB
AB
AB
CC C
AB
C
(A.4)
PI or PID controller (augmented state, input, and output matrices)
22 2 1 21
12
12
12
ˆ ˆ
,
0 0
ˆ
0
ˆ
1 (isolates integral term)
nn n n
n
n
In
A0 B
AB
C
CC
C0
(A.5)
where the block matrix components for A, B, and C are as follows
For rigid-body modes
00 00
0
0 01
,, ()0
() 00
s
a
AB Cvx
vx
(A.6)
For flexible modes
2
01 0
,, ()0
2()
ii iis
iii ia
AB Cvx
vx
(A.7)
141
For P or PD controller substituting (A.6) and (A.7) into (A.4) yields
0
2
111 1
2
01
01 0
00 ( )
01 0
2()
,
0
01 0
2()
( ) 0 ( )0 0 ( )0
a
a
nnn na
ss ns
pd
vx
vx
vx
vx v x v x
KK
AB
C
CC CA
(A.8)
For PI or PID controller substituting (A.6) and (A.7) into (A.5) yields the augmented
state, input, and output matrices
0
2
111 1
2
01
01
01 0 0
00 0 ()
01 0
2 ()
ˆ ˆ
,
01 0 0
20 ()
() 0 ( ) 0 ( ) 0 0 0
ˆ
() 0 ( ) 0 () 0 0
ˆ
0 0 0 0 0 0 1 (isolates integral term)
a
a
nnn na
ss ns
ss ns
I
vx
vx
vx
vx v x v x
vx v x v x
AB
C
C
C
ˆˆˆ ˆ
pd II
KK K CCA C
(A.9)
142
For P or PD controller : evaluate
1
0
, () , and ()
t
tt
et ed t e I
AA A
Γ B Ψ
Rigid-body modes
0
1
1
0
1
1
1
1
2
1
01
00
1
0
11
1
0
At
eL sI A
LsI
s
L
s
ss
L
s
(A.10)
0
1
01
At
t
e
(A.11)
Since
0
A is singular (ref (A.6)), the integral form of ( ) t Γ is used to calculate
0
() t
0
2
00 0
0
00
0 1
() ( )
2
() 01
tt
A
a
a
t
tedB d vx
vx
t
(A.12)
If the relay contains a dead-zone, the above integral term is evaluated in a more general
form as follows.
22
0
11
22
21
00 0
0
21
1
0 1
() ( ) 2
() 01
tt
A
a
a
tt
tt
tedB d vx
vx
tt
(A.13)
143
Flexible modes
1
1
1
1
2
1
1
2
222 2
1
2
222 2
01
2
1
2
2 1
22
22
i
At
i
iii
iii
ii
ii i i i i
i
ii i i i i
eL sI A
LsI
s
L
s
s
ss ss
L
s
ss ss
(A.14)
2
1
cos( ) sin( ) sin( )
sin( ) cos( ) sin( )
i i
ii i
ii
i
ii
ii i
i i
tt i
dd d
dd
At
tt ii
dd d
dd
et t e t
e
et e t t
(A.15)
where
2
, and 1
i
iii d i i
Since
i
A is nonsingular (ref (A.7)), the matrix form of ( ) t Γ is used to calculate ( )
i
t
1
2
1
11 12
2
21 22
1
12
2
22
22
()
01 0 1
2()
1
01
() (substituting )
2
1
1
sin( )
2 1
10 cos(
i
i
i
i
i
At
ii i
ii
ii
iii ia
i
ia i i i
i
ii
t
d
i
d
ii
t
d
tA e I B
ee
vx
ee
e
vx
e
et
e
()
)sin( )1
ii
i
ia
i
d
d
vx
tt
144
2
1cos() sin()
() ( )
1
sin( )
i
ii
i
i
i
i
t i
dd
d
iia i
t
d
d
et t
tvx
et
(A.16)
where
11 12 21 22
, , , and
iii i
eee e are from the exponential matrix (A.15) for the ith mode.
The resulting vector ( ) t Γ is
0
1
()
()
()
()
n
t
t
t
t
(A.17)
If the relay contains a dead-zone the above integral terms for () t Γ are evaluated in a
more general form as follows.
22
11
2
1
2
1
11 12
21 22
1
2
0
()
()
1
sin( )
() ()
cos( ) sin( )
i
i
i
i
i
ii
i
ii tt
A
ii
ii
ia
tt
t
d
d t
i
ia ia
t
i
i
dd
d t
ee
tedB d
vx
ee
ed
vx vx
ed
(A.18)
where
11 12 21 22
, , , and
iii i
eee e are from the exponential matrix (A.15) for the ith mode.
145
2
1
2
1
2
1
1
22
22
1
22
11
and where
1
sin( )
sin( ) cos( )
1
sin( ) cos( )
1
sin( ) cos( )
and
i
i
i
ii i
i
ii
i
ii i
i
ii
ii i
t
id
d t
t
id d d
did
t
t
id d d
i
t
di d
id d d
ed
e
et t
et t
2
1
2
1
2
1
22 2
2
22
22
22
11
cos( )
cos( ) sin( )
cos( ) sin( )
1
cos( ) sin( )
i
i
ii i
i
i
i
ii i
i
i
ii i
ab
ii i
t
a
id
t
t
id d d
id
t
t
id d d
t
id
id d d
i
ed
e
et t
et t
2
1
22
2
22
11
sin( ) cos( )
1
sin( ) cos( )
i
ii i
i
i
ii i
t
dd i d
a
t
id
dd i d
et t
et t
2
1
2
1
2
1
2
22
22
22
11
and
sin( )
sin( ) cos( )
sin( ) cos( )
sin( ) cos( )
i
i
i
ii i
i
ii
i
ii i
i
ii
ii i
t
b i
id
d t
t
id d d
i
did
t
t
id d d
i
t
di d
id d d
i
ed
e
et t
et t
2
1
22
2
22
11
sin( ) cos( )
sin( ) cos( )
i
ii i
i
ii
ii i
t
id d d
b i
t
di d
id d d
et t
et t
146
The inverse term
1
At
eI
is block diagonal also, and can thus be evaluated as
follows.
0
1
n
At
At
t
A t
e
e
e
e
A
(A.19)
0
1
0
1
0
1
1
2
1
2
2
1
2
2
2
1
2
1
2
1
2
()
n
n
n
At
At
At
At
At
At
At
At
At
At
I e
I e
te I
I e
eI
eI
eI
eI
eI
eI
0
1
()
n
t
(A.20)
0
where , and are as follows
i
147
Rigid-body modes
0
1
1
02
1
110
()
01 0 1
11
2
24
02 1
0
2
At
t
te I
t
t
(A.21)
Flexible modes
1
1
11 12
2
21 22
1
11 12
21 22
22 12
11 22 12 21 11 22 12 21
11 12
21 22 21
11 22
10
()
01
1
1
1
11 11
()
1
i
ii
At
i
ii
ii
ii
ii
ii ii ii ii
ii
i
ii i
ii
ee
te I
ee
ee
ee
ee
ee ee ee ee
t
e
ee
11
12 21 11 22 12 21
1
111
i
ii i i ii
e
ee e e ee
(A.22)
where
11 12 21 22
, , , and
iii i
eee e are from the exponential matrix (A.15) for the ith mode.
148
For PI or PID controller: evaluate
1
ˆˆ ˆ
0
ˆ ˆˆ
, () , and ()
t
tt
et ed t e I
AA A
Γ B Ψ
The above results for a PD controller can be used with the augmented state (A.5) repeated
below to form the exponential matrix, gamma, and psi matrices
ˆ
ˆ ˆ
, (), and ()
t
et t
A
for
the augmented system.
22 2 1 21
1
12
12
ˆ ˆ
,
0 0
ˆ
0
ˆ
1 (isolates integral term)
nn n n
n
n
In
A0 B
AB
C
CC
C0
(A.23)
1
ˆ
1
21
1
22 2 1 1
21
12
1
22 2 1 2
1
12
ˆ
0
t
n
nn n
n
n
nn n n
n
eL sI
LsI
sI
L
s
A
A
A0
C
A0
C
(A.24)
The inverse of the above partitioned matrix (A.24) can be solved for by using a 2x2 block
matrix
1
11 11 1
21 33 21 33
11
21 33
11 11
21 33 21 33
0
0
0
0
mn mm m n
mm
nm nn
nm n n
mm m n
mm mm
nm n n
nm nn n m nn
I XY
I ZU
XY I
XZ Y U I
a0 a 0
AA I
aa aa
a0
aa
aa
aa a a
(A.25)
Solving for the four unknown block components X, Y, Z, and U yields
149
1
11 11
11
11
21 33 33 21 11
1
21 33 33
00
0
mm m m
mm mm
mn m n mn
mm
nm n m
nm nn n n nm m m
nn nn
nm nn nn
XI X
YY
XZ Z
YUI U
aa
a
aa a a a
aa a
(A.26)
Substituting matrix components from (A.24), the resulting inverse is
1
11
1
11 1
33 21 11 33
1
22 2 1 2
1
11
12 2 2 2
mn
mm
nn nm mm nn
nn n n
nnn n
XY
ZU
sI
ssI s
a0
A
aa a a
A0
CA
(A.27)
Using (A.27) in (A.24) yields
ˆ
11
1
1
22 2 1 2
1
11 1
12 2 2 2
ˆ 11 2 1 21
22
21 22 21
12 12
1
ˆ
ˆˆ ˆ 1
t
nn n n
nnn n
t
n n t nn
n n
eL
LsI
Ls sI L
s
e e
e
ee e
A
A
A
A
A0
CA
0 0
(A.28)
where
11
22
ˆ
nn
e
the exponential matrix (A.19) from above and
22
ˆ 1 e . The individual
terms of the off-diagonal component
21
12
ˆ
n
e
can be found as follows.
150
1
11
21 1 2 2 2 2 12
01
1
0
1
1
1
21 21 21 21
01 2
ˆ
1
ˆˆ ˆ ˆ
nnn n n
n
i
n
n
eLs sI
CC C
s
sI A
L
sI A
sI A
ee e e
CA
(A.29)
where terms for the output matrix C are taken from (A.6) and (A.7). The terms for
1
22 2 nn n
sI
A are found in (A.10) and (A.14). The individual terms of (A.29) are
Rigid-body modes (ref (A.10))
2
1
11
21 0 0 0
0
11
00 23
12
2
21 21 21 0 0
00 0
11
11
ˆ () 0
1
0
11
() ( )
1
ˆˆ ˆ () ( )
2
s
ss
ss
ss
eL CsIA L vx
ss
s
Lv x L v x
ss
ee e vxtvxt
(A.30)
Flexible modes (ref (A.14))
1
1
21
222 2
1
2
222 2
11
22 22
1
ˆ
2 1
22
1
() 0
22
2 1
() ()
22
ii
i
ii
ii i i i i
is
i
ii i i i i
ii
is is
ii i i i i
eL CsIA
s
s
ss ss
Lvx
s s
ss ss
s
Lvx L vx
ss s s s s
12
21 21 21
ˆˆ ˆ
ii i
ee e
(A.31)
151
where the components of (A.31) are
1
21 2
2
21 2
2 1
ˆ () sin( ) 1 cos( ) sin( )
1
ˆ () 1 cos( ) sin( )
ii
ii i
i i
i
ii
i
ttii
is d d d
i
di d
t i
is d d
i
id
evxe t e t t
evx e t t
(A.32)
The matrix exponential and gamma term from the previous results can be used to
calculate the augmented gamma term as follows
ˆ
21
21
00 12
()
ˆˆ
()
ˆˆ 1 0
t tt
n
n
e t
ted d
e
A
A
0 B Γ
Γ B
γB
(A.33)
where the term ˆ γB is evaluated using (A.31) as follows
0
1
21 0 1
12
0
00 1 1
ˆˆˆˆ ˆ
ˆˆ ˆ
t
n
n
n
nn
B
B
ed
B
BB B
γBB
(A.34)
where for
Rigid-body modes (ref (A.6) and (A.30))
12
0 0 21 21
00
0
00
2
00
0
3
00
0
ˆ ˆˆ
()
1
() ( )
2
1
() ( )
6
tt
a
t
sa
sa
Be d e d
vx
vx td vx
vx vx t
(A.35)
152
Flexible modes
12 2
21 21 21
00 0
2
0
0
ˆ ˆˆ ˆ ()
()
1
() ( ) 1 cos( ) sin( )
ˆˆˆˆ () ( )
i
ii
i
tt t
ii ia
ii i
ia
t
t i
is i a d d
id
ab c
ii i s i a i i i
B e d ed edvx
vx
vx vx e t t d
Bvxvx
(A.36)
2
0
2 0
2
2
0
222
0
22 2
22
where
1
ˆ
1
1
ˆ
and
1
ˆ cos( )
cos( ) sin( )
1
11
cos( ) sin( )
1
ˆ
i
i
ii i i
i
i
ii i
i
t
a
i
i
t
i
a
i
i
t
b
id
i
t
id d d
iid
t
id d d i
ii d
b
i
ii
d
t
ed
e
et t
2
2
0
222
0
2 22
sin( ) cos( )
and
1
ˆ sin( )
sin( ) cos( )
1
1
sin( ) cos( )
i
ii i
i
i
i
i
ii i i
ii
i
ii i i
ii
t
dd i d i
d
t
c i
id
id
t
id d d
i
id i d
t i
id d d d
i di d
et t
ed
e
et t
22 2
ˆ sin( ) cos( )
i
ii i i
ii
t c i
iidddd
id i d
et t
where
2
, and 1
i
iii d i i
.
153
If the relay contains a dead-zone the above integral terms for (A.34) are evaluated in a
more general form as follows.
Rigid-body modes (ref (A.6) and (A.30))
22
11
2
1
12
0 0 21 21
00
0
2
00
33
00 2 1
0
ˆ ˆˆ
()
1
() ( )
2
1
() ( )
6
tt
a tt
t
sa
t
sa
Be d e d
vx
vx td vx
vx vx t t
(A.37)
Flexible modes
22 2
11 1
2
1
12 2
21 21 21
2
0
ˆ ˆˆ ˆ ()
()
1
() ( ) 1 cos( ) sin( )
ˆˆˆˆ () ( )
i
ii
i
tt t
ii ia
ii i
ia tt t
t
t i
is i a d d
id t
ab c
ii i s i a i i i
Be d e d e dvx
vx
vx vx e t t d
Bvxvx
(A.38)
where
2
1
2
1
2
2
21 2
1
ˆ
1
1
ˆ
t
a
i
i t
t
t
i
a
i
i
d
tt
154
2
1
2
1
2
1
2
222
22
22 2
11
and
1
ˆ cos( )
cos( ) sin( )
1
cos( ) sin( )
11
cos( ) sin( )
ˆ
i
i
ii i i
i
i
ii i
i
i
ii i
t
b
id
i t
t
id d d
iid
t
t
id d d
t
ii d
id d d
i
ed
e
et t
et t
2
1
2
1
22
22 2
11
2
222
sin( ) cos( )
1
sin( ) cos( )
and
1
ˆ sin( )
sin( ) cos( )
1
i
ii i
i
i
ii i
i
i
i
ii i i
ii
t
dd i d
b
t
ii d
dd i d
t
c i
id
id t
id d d
i
id i d
et t
et t
ed
e
2
1
2
1
2
1
22
2 22
11
22
22 2
11
sin( ) cos( )
1
sin( ) cos( )
sin( ) cos( )
ˆ
sin( ) cos(
i
ii i
i
ii
ii i
i
ii i
i
ii
ii i
t
t
t
id d d
i
t
i di d
id d d
t
id d d
c i
i
t
id i d
id d d
et t
et t
et t
et t
)
where
2
, and 1
i
iii d i i
.
The inverse term
1
ˆ
t
eI
A
for the augmented system can be evaluated using (A.28)
as follows.
1
1
ˆ
21 2 2 21
21 1 2
12
1
21
21
12
ˆ
()
ˆ 11
ˆ 2
t
nnn n t
n
n
t
n
n
e
te I
e
e
e
A
A
A
0I 0
Ψ
0
I0
(A.39)
155
The block matrix inverse of (A.39) can be evaluated using (A.27) which is repeated
below
1
11
1
11 1
33 21 11 33
mn
mm
nn nm m m nn
a0
A
aa a a
(A.40)
Substituting block matrix components from (A.39) yields
1
21
22
1
1
21
12
22
ˆ
()
11
ˆ
22
t
n
nn
t
n
nn
e
t
ee
A
A
I0
Ψ A
I
(A.41)
Note that the upper left component of (A.41) is the same as the above (A.20).
Substituting this results in
22 2 1
11 2 1 1 22
21 22
21 2 2 12
12
()
ˆ
ˆ
()
11
ˆˆ
ˆ ()
22
nn n
n
nn
nn n
n
t
t
et
Ψ 0
0
Ψ A
Ψ
(A.42)
The off-diagonal term of (A.42) can be evaluated using
21
12
ˆ
n
e
from (A.29) and
22
()
nn
t
Ψ from (A.20) as follows
21 21 2 2
12 12
0
1
21 21 21
01
21 21 0 21 1 21
12 0 1
1
ˆ ˆ ()
2
1
ˆˆ ˆ
2
11 1
ˆ ˆˆ ˆ
22 2
nn
nn
n
n
n
nn
et
ee e
ee e
Ψ
(A.43)
156
Rigid-body modes (ref (A.30) and (A.21))
2
21 0 0 0
0
0
11
11 1 1
24
ˆ () ( )
1 22 2 2
0
2
1
() 0
4
ss
s
t
evxt vxt
vx t
(A.44)
Flexible modes (ref (A.31) and (A.22))
12
11 12
21 21 21
21 22
11 1
ˆˆ ˆ
22 2
ii
i
ii i
ii
ee e
21
12 1 2
21 11 21 21 21 12 21 22
1
ˆ
2
11
ˆˆ ˆˆ
22
i
i
ii ii
ii ii
e
ee ee
(A.45)
where from (A.22)
22 12
11 22 12 21 11 22 12 21
11 12
21 22 21 11
11 22 12 21 11 22 12 21
1
11 11
1
11 11
ii
ii ii ii ii
ii
ii ii
ii ii ii ii
ee
ee ee e e ee
ee
ee ee e e ee
where
11 12 21 22
, , , and
ii i i
eee e are the matrix elements of (A.15) for the ith mode and from
(A.32)
1
21 2
2
21 2
2 1
ˆ ( ) sin( ) 1 cos( ) sin( )
1
ˆ () 1 cos( ) sin( )
ii
ii i
i i
i
ii
i
ttii
is d d d
i
di d
t i
is d d
i
id
evxe t e t t
evx e t t
where
2
, and 1
i
iii d i i
.
157
Appendix B: Switching Constraint Equations for Other Relay
Types
The following Appendix extends the derivations of switching constraint equations shown
in Section 3.3 to other relay types.
B.1 Relay with Hysteresis and Pure Time Delay
Limit cycles in systems described in Section 3.3 that contain a pure time delay in the
feedback loop can be evaluated using a modified version of (3.26) as follows. In the
system shown in Figure 3.12, the relay input is the negative of the system output y(t).
However, when a time delay is present in the feedback signal, the resulting relay
switching will be delayed by a time delay constant τ
d
, as shown in Figure B.1.
Figure B.1: Limit Cycle Waveform with Pure Time Delay in Relay Feedback
The resulting state dynamics during the limit cycle half-period t* for a relay with
hysteresis and delayed input are then
t = 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
1
0
1
1
0
1
Limit Cycle Waveform (Relay Input u(t), Linear Plant Output y(t))
t
yt ()
ut ()
Δ t ()
Δ t ()
t
t = t* = T/2
t = 2t* = T
τ
d
y(t) u(t – τ
d
) = ±M
158
*
for 0
for
d
d
xx M t
x xM t t
AB
AB
(B.1)
The state solution (3.23) at the switching instances, and (3.22) are used to determine the
initial conditions that will start the limit cycle.
At
d
t
0
() e (0) ( )
d
d
d
x xed M
A A
B (B.2)
At t = t*
*
*
() *
0
() e ( ) ( )
d
d
t
t
d
xtx edM
A A
B (B.3)
Using (B.2) in (B.3) and the symmetry constraint
*
() (0) xtx yields the following
initial condition for the limit cycle
*
*
*
* *
*
* *
*
* *
() *
00
()
00
()
00
1
(
0
() e e (0) ( ) ( )
(0) e (0) e ( ) ( )
e(0) e
(0) e
dd
dd
dd
d
dd
d
d
t
t
t
t t
t
t t
t
t t
xtxedMedM
xx ed M ed M
xedM edM
xede
AA AA
A AA A
A AA A
A AA
BB
BB
IB B
IB
)
0
d
d
dM
B
* *
* *
*
1
*
()
00
(0)
where e , , 2
dd
d
d
t t
t t
t
xMx
ed e d ed
A AA A
Φ I Γ1 Γ2
ΦΓ1B Γ BB
(B.4)
159
The above was derived assuming
*
d
t . For the case when
*
d
t , equation (B.4) is
modified as follows (Astrom [4])
* *
1
*
*
1
1
*
0
**
1
ˆˆ
(0) ( 1)
ˆˆ
where e , , 2 ,
interger part of , remainder of
m
t t
t
t
dd
xMx
ed ed
mt t
AA A
Φ I Γ1 Γ2
ΦΓ1B Γ B (B.5)
The above initial condition yields the following constraint equation that must be satisfied
during the limit cycle.
1
(0) (0) y x
M
C
C Φ I Γ1 Γ2
* *
*
*
1
0
(0) e
d
d
t t
t
t
yededM
AA A
CI B B (B.6)
Solutions of t* which satisfy (B.6) must be obtained numerically, and determine possible
limit cycle frequencies
**
t . The waveform must also satisfy (3.27).
160
B.2 Relay with Linear-Zone (Saturation Nonlinearity)
Figure B.2: Relay with Linear-Zone (Saturation Nonlinearity)
For the case when the relay contains a linear zone (Figure B.2), the following state-space
form is used.
if
() if
( ) otherwise,
() () ()
() ( ( ))
Myt
ut M y t
kx t k M
xt xt ut
yt x t
AB
(B.7)
In order for a symmetric unimodal limit cycle to exist in system (B.7), the following
waveform and switching conditions must hold.
M
δ
ey
u
M
k
161
Figure B.3: Limit Cycle Waveform Switching Criterion for a Relay with Linear-Zone
Note in the above Figure B.3 that the output y(t) crosses the positive saturation boundary
at
*
tt (linear-zone reset time) with negative slope, and the negative saturation
boundary at
*
tt (limit cycle half-period) with negative slope. Therefore, this type of
relay requires that two switching conditions be satisfied during the half-cycle, in order for
a symmetric limit cycle to exist.
For
*
0 tt , u(t) = –M (constant saturation magnitude)
()
0
(0)
t
tt
x xM
x ex e Md
AA
AB
B
and at
*
tt
*
*
*
0
() (0)
t
t
xte x e Md
AA
B (B.8)
t*
*
t
T = 2t*
y(t) u(t)
2
162
Note that as the linear-zone width is decreased to zero, the slope parameter k ,
the reset time parameter
**
1, and tt (ref Figure B.3). Equation (B.8) then
matches the switching requirement for an ideal relay as expected.
For
**
tt t , () () () ut ky t k xt C (controller in linear range)
*
1
**
1
11
() *
*
() **
where =
()
and at
() ( )
tt
tt
xxk x k x
xk
xe x t
tt
xt e x t
A
A
ABC A BC
AAABC
(B.9)
Substituting (B.8) into (B.9) yields the required state at the half-period switching time t*.
*
**
**
1
*
****
1
*()
0
() **
()*()
0
() (0)
() ( )
() (0)
t
tt
tt
t
tttt
xt e x e Md
xt e x t
xte e x e Md
AA
A
A AA
B
B
*
** ** *
11
() ( ) *
0
() (0)
t
tt tt t
xte e x e e Md
AA AA
B (B.10)
Again note that as the linear-zone width is decreased to zero, the slope parameter
k , the reset time parameter
** * *
1, , and 0 tt t t . Equation (B.10)
then matches the switching requirement for an ideal relay as expected.
The initial state which results in a limit cycle can be solved for using (B.10) and
noting that due to symmetry
*
() (0) xtx .
163
*
** * * *
11
() ( ) *
0
() (0) (0)
t
tt tt t
xtx e e x e e Md
AA AA
B
rearranging yields the initial state x*
*
** * * *
11
1
() ( ) *
0
(0)
t
tt t t t
x xe e Ie e Md
AA AA
B (B.11)
The initial state (B.11) can now be used in (B.8) to solve for the second switching
criterion at the relay reset time
*
tt .
*
*
**
** ** **
11
*
0
1
() ( )
00
() (0)
t
t
tt
tt tt tt
xt e x e Md
ee e I e e Md e Md
AA
AA AA A A
B
BB
(B.12)
Using (B.11) and (B.12), solve the following two simultaneous equations. Valid
solutions must also satisfy the waveform criteria at each switch
*
**
(0)
() ()
yx
yt x t
C
C
(B.13)
Waveform Criteria (refer to Figure B.3).
*
*
**
*
**
(0)
(0) 0 ensures proper velocity direction at switch 1
( ) 0 ensures proper waveform
() ()
( ) 0 ensures proper velocity direction at switch 2
( ) ensures pro
yx
y
yt t t
yt x t
yt
yt t t t
C
C
per waveform
(B.14)
The above simultaneous equations (B.13) must be solved using numerical methods to get
values of t* and γt*.
Abstract (if available)
Abstract
A method for calculating all periodic solutions and their domains of attraction for flexible systems under nonlinear feedback control is presented. The systems considered consist of mechanical systems with multiple flexible modes and a relay type controller coupled with a linear or nonlinear control law operating in a feedback configuration. The proposed approach includes three steps. First, limit cycle frequencies and periodic fixed points are computed exactly, using a block-diagonal state-space modal representation of the plant dynamics. Then the relay switching surface is chosen as the Poincaré mapping surface and is discretized using the cell mapping method. Finally, the region of attraction for each limit cycle is computed using the cell mapping algorithm and employing an error based convergence criterion. The method is demonstrated on four application examples, each with a different type of relay and/or control law to demonstrate the versatility and accuracy of the method.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Robust control of periodically time-varying systems
PDF
Nonlinear control of flexible rotating system with varying velocity
PDF
Active delay output feedback control for high performance flexible servo systems
PDF
Practical adaptive control for systems with flexible modes, disturbances, and time delays
PDF
Transient modeling, dynamic analysis, and feedback control of the Inductrack Maglev system
PDF
Modeling and dynamic analysis of coupled structure-moving subsystem problem
PDF
Adaptive control with aerospace applications
PDF
Medium and high-frequency vibration analysis of flexible distributed parameter systems
PDF
Control of spacecraft with flexible structures using pulse-modulated thrusters
PDF
An analytical dynamics approach to the control of mechanical systems
PDF
Modeling, analysis and experimental validation of flexible rotor systems with water-lubricated rubber bearings
PDF
Distributed adaptive control with application to heating, ventilation and air-conditioning systems
PDF
Periodic motion analysis in spacecraft attitude dynamics
PDF
Coordinated freeway and arterial traffic flow control
PDF
Analytical and experimental studies in system identification and modeling for structural control and health monitoring
PDF
Integrated control of traffic flow
PDF
Control of mainstream traffic flow: variable speed limit and lane change
PDF
Optimal design, nonlinear analysis and shape control of deployable mesh reflectors
PDF
Flexible formation configuration for terrain following flight: formation keeping constraints
PDF
High-accuracy adaptive vibrational control for uncertain systems
Asset Metadata
Creator
Borre, Michael K.
(author)
Core Title
Periodic solutions of flexible systems under discontinuous control
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Mechanical Engineering
Publication Date
04/28/2011
Defense Date
03/29/2011
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
flexible structures,limit cycles,numerical-analytical methods,OAI-PMH Harvest,region of attraction,relay feedback control
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Flashner, Henryk (
committee chair
), Ioannou, Petros A. (
committee member
), Yang, Bingen (
committee member
)
Creator Email
borre@usc.edu,mborre@mbengineering.com
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m3803
Unique identifier
UC149345
Identifier
etd-Borre-4560 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-464519 (legacy record id),usctheses-m3803 (legacy record id)
Legacy Identifier
etd-Borre-4560.pdf
Dmrecord
464519
Document Type
Dissertation
Rights
Borre, Michael K.
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
flexible structures
limit cycles
numerical-analytical methods
region of attraction
relay feedback control