EXPLOITING IMPULSIVE INPUTS FOR STABILIZATION OF UNDERACTUATED

ROBOTIC SYSTEMS: THEORY AND EXPERIMENTS

By

Nilay Kant

A DISSERTATION

Michigan State University

in partial fulfillment of the requirements

Submitted to

for the degree of

Mechanical Engineering – Doctor of Philosophy

2020

ABSTRACT

EXPLOITING IMPULSIVE INPUTS FOR STABILIZATION OF UNDERACTUATED

ROBOTIC SYSTEMS: THEORY AND EXPERIMENTS

By

Nilay Kant

Robots have become increasingly popular due to their ability to perform complex tasks and operate
in unknown and hazardous environments. Many robotic systems are underactuated i.e., they have
fewer control inputs than their degrees-of-freedom (DOF). Common examples of underactuated
robotic systems are legged robots such as bipeds, flying robots such as quadrotors, and swimming
robots. Due to limited control authority, underactuated systems are prone to instability. This work
includes impulsive inputs in the set of admissible controls to address several challenging control
problems. It has already been shown that continuous-time approximation of impulsive inputs can
be realized in physical hardware using high-gain feedback.

Stabilization of an equilibrium point is an important control problem for underactuated systems.
The ability of the system to remain stable in the presence of disturbances depends on the size of
the region of attraction of the stabilized equilibrium. The sum of squares and trajectory reversing
methods are combined to generate a large estimate of the region of attraction. This estimate is then
effectively enlarged by applying the impulse manifold method to stabilize equilibria from points
lying outside the estimated region of attraction. Simulation results are provided for a three-DOF
tiptoebot and experimental validation is carried out on a two-DOF pendubot.
Impulsive inputs
are also utilized to control the underactuated inertia-wheel pendulum (IWP). When subjected to
impulsive inputs, the dynamics of the IWP can be described by algebraic equations. Optimal
sequences of inputs are designed to achieve rest-to-rest maneuvers and the results are applied to
the swing-up control problem. The novel problem of juggling a devil-stick using impulsive inputs
is also investigated. Impulsive forces are applied to the stick intermittently and the impulse of the
force and its point of application are modeled as inputs to the system. A dead-beat design for one of
the inputs simplifies the control problem and results in a discrete linear time invariant system. To

achieve symmetric juggling, linear quadratic regulator (LQR) and model predictive control (MPC)
based designs are proposed and validated through simulations.

A repetitive motion is described by closed orbits and therefore, stabilization of closed orbits is
important for many applications such as bipedal walking and steady swimming. We first investigate
the problem of energy-based orbital stabilization using continuous inputs and intermittent impulsive
braking. The orbit is a manifold where the active generalized coordinates are fixed and the total
mechanical energy of the system is equal to some desired value. Simulation and experimental
results are provided for the tiptoebot and the rotary pendulum, respectively. The problem of orbital
stabilization using virtual holonomic constraints (VHC) is also investigated. VHCs are enforced
using a continuous controller which guarantees existence of closed orbits. A Poincaré section is
constructed on the desired orbit and the orbit is stabilized using impulsive inputs which are applied
intermittently when the system trajectory crosses the Poincaré section. This approach to orbital
stabilization is general, and has lower complexity and computational cost than control designs
proposed earlier.

This dissertation is lovingly dedicated to parents; my father Shri. Rajni Kant, and to my mother
Smt. Kumari Jyotsna Kant as without their sacrifices, I would not have been able to accomplish
this. I would also like to dedicate this work to Swami Vivekananda, an Indian yogi whose life,

works and teachings have greatly transformed my life.

iv

ACKNOWLEDGEMENTS

I would like to express my utmost gratitude to my PhD advisor, Dr. Ranjan Mukherjee, as his
influence is as profound on this dissertation as it is in my life. In the ancient Indian civilization,
education was not very formal as it is today and a student spent several years at a teacher’s house
to attain knowledge. The close association with the teacher resulted in a strong human relationship
and the student mastered the most subtle aspect of a subject along with the desired human values.
I am extremely lucky to have experienced such a process in the last five years, and this has been
possible due to my association with Dr. Mukherjee. His office was always open during the day and
I could simply walk inside anytime if I was stuck in a proof. There were countless occasions where
both of us brainstormed together on a topic and hours just passed by. I vividly remember that he
once showed me how to derive the Lagrange’s equation on a notebook with minute details such as
defining suitable coordinate system and drawing clear diagrams using a pencil. These instances
may appear insignificant, but for me, they were life changing. It was only due to such involved
and patient training, I became capable enough to accomplish this work. Despite over twenty five
years of being a professor, Dr. Mukherjee unfailingly comes to office every weekend. Such hard
work coupled with a genius mind inspired me throughout to push my boundaries and would forever
continue to do so. I have experienced a transformation in the last five years - admirable traits of
mine, if any, are due to Dr. Mukherjee, and the rest are my fault. Additionally, I am sincerely
thankful to my PhD committee members Dr. Hassan Khalil, Dr. Vaibhav Srivastava and Dr.
Gouming Zhu for their continued guidance, motivation and support. I am also very grateful to my
lab mates; Sheryl, Mahmoud, Sanders, Connor, Krissy and Amer as their friendship did not let
me miss my family back in India. I would also like to thank Mr. Roy Bailiff of the mechanical
engineering department whose help in the machine shop ensured me to build my experimental
setup.

I would like to express my gratitude to my parents are they have made several sacrifices to
ensure that I get the best education possible. If it was not their conditional love and support, I

v

would not have come this far. I am also very thankful to my brother Dr. Kislay Kant as he has
always set a high standard due to his achievements and conduct which has greatly inspired me
to dream big and sincerely work towards it. I would also like to pay my salutations to my late
grandparents Shri. Ramaeshward Prasad, Shri Hari Narayan Prasad, Smt. Rajkishori Devi and
Smt. Gyneshwari Devi as their morals and life struggles were the source of inspiration in difficult
times. Lastly, I would like to express my deepest gratitude to my wife Priyanka, as she was my
strongest support throughout, and her love made me finish this journey with ease.

vi

TABLE OF CONTENTS

.
.

xi

.
.
.
.

.
.
.
.

.
.
.
.

LIST OF TABLES .
.
.
LIST OF FIGURES .
.
.
LIST OF ALGORITHMS .
KEY TO SYMBOLS .
.
.
CHAPTER 1

.
.
.
.
INTRODUCTION .
.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvi
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii
1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1
1.1 Motivation .
.
2
1.2 Estimating and Enlarging the Region of Attraction of Equilibria . . . . . . . . . . .
1.3 Rest-to-Rest Maneuver of the Inertia Wheel Pendulum . . . . . . . . . . . . . . . .
4
5
.
1.4 Devil-Stick Juggling .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
1.5 Energy-Based Orbital Stabilization . . . . . . . . . . . . . . . . . . . . . . . . .
.
1.6 Orbital Stabilization using Virtual Holonomic Constraints . . . . . . . . . . . . . .
8

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

.

.

.

.

.

.

.

.

CHAPTER 2 ESTIMATION OF THE REGION OF ATTRACTION OF EQUILIBRIA

.
.

.
.

.
.

.
.

. .
.
.

2.1
Introduction .
2.2 Background .

AND ITS ENLARGEMENT USING IMPULSIVE INPUTS . . . . . . . . . 10
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 The Impulse Manifold Method . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.2 The Sum of Squares Method . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Enlarging the Region of Attraction Using SOS and IMM . . . . . . . . . . . . . . 15
2.3.1 An Illustrative Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.3.2
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

Simulation Results

2.4 Algorithm for Computing ReA Using the Method of

.
.
.

Trajectory Reversing .
. .
2.4.1 Definitions
2.4.2 Algorithm .
.
.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.5 Enlarging the Region of Attraction Using SOS, CHART and IMM . . . . . . . . . 25
2.6 Experimental Verification of Enlargement of ReA using SOS, CHART and IMM . . 29
2.6.1 Hardware Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
2.6.2 Design of Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.6.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

CHAPTER 3

Introduction .

IMPULSIVE DYNAMICS AND CONTROL OF THE INERTIA-WHEEL
PENDULUM .
.
.

3.1
3.2 Mathematical Model and Problem Statement . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
. 37
3.2.1 Equations of Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.2.2 Effect of Impulsive Input and Constants of Free Motion . . . . . . . . . . . 37
3.2.3
Problem Statement: Rest-to-Rest Maneuvers . . . . . . . . . . . . . . . . . 38

. .

.
.

.
.

.
.

vii

3.3 Rest-to-Rest Maneuvers: Case of One and Two Impulsive Inputs

. . . . . . . . . . 39
3.3.1 Case of One Impulsive Input (N “ 1) . . . . . . . . . . . . . . . . . . . . . 39
3.3.2 Case of Two Impulsive Inputs (N “ 2) . . . . . . . . . . . . . . . . . . . . 40
3.4 Rest-to-Rest Maneuvers: Generalization to N Inputs . . . . . . . . . . . . . . . . . 44
3.4.1 Revisiting the Problem Statement
. . . . . . . . . . . . . . . . . . . . . . 44
3.4.2 Rest-to-Rest Maneuvers with Even Number of Impulsive Inputs
. . . . . . 45
. . . . . . . . . . . . 48
3.4.3 Rest-to-Rest Maneuvers with Odd Number of Inputs
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.5.1 Optimal Swing-Up Using Even Impulse Sequences . . . . . . . . . . . . . 50
3.5.2
Implementation Using High-Gain Feedback . . . . . . . . . . . . . . . . . 50
3.5.3 Discussion of Results in the Literature . . . . . . . . . . . . . . . . . . . . 51
. . . . . . . . . . . . . . . . . . 51
Energy Based Controller . . . . . . . . . . . . . . . . . . . . . . 55

3.5.3.1 Globally Stabilizing Controller
3.5.3.2

3.5 The Swing-Up Problem .

.

CHAPTER 4

.
.
.

.
.
.

. .

Introduction .

Impulsive Dynamics

IMPULSIVE CONTROL OF A DEVIL-STICK: PLANAR SYMMETRIC
JUGGLING .
.
.

. 58
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.1
.
4.2 Problem Description .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.3 Dynamics of the Devil-Stick . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3.1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.3.2 Continuous-time Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.3.3
Poincaré Sections and Half-Return Maps . . . . . . . . . . . . . . . . . . . 61
4.3.4 Coordinate Transformation . . . . . . . . . . . . . . . . . . . . . . . . . . 62
4.3.5
. . . . . . . . . . . . . . . . 64
4.4 State Feedback Control Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
Steady-State Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.4.1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.4.2 Error Dynamics .
4.4.3
. . . . . . . . . . . . . . . . . 68
4.4.4 Residual Control Design . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.5.1
. . . . . . . . . . . . . . . . . . 72
4.5.2 Results for the LQR-based Design . . . . . . . . . . . . . . . . . . . . . . 72
4.5.3 Results for the MPC-based Design . . . . . . . . . . . . . . . . . . . . . . 74

Single Return Map and Discrete-Time Model

Partial Control Design: Dead-Beat Control

System Parameters and Initial Conditions

4.5 Simulation Results

.

.

.

.

.

CHAPTER 5 ENERGY-BASED ORBITAL STABILIZATION USING CONTINUOUS

.

.
.

.
.

.
.

.
.

Introduction .

INPUTS AND IMPULSIVE BRAKING . . . . . . . . . . . . . . . . . . . 77
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.1
. .
5.2 Problem Statement
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
5.3 Modeling of System Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.1 Lagrangian Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3.2 Effect of Impulsive Inputs
. . . . . . . . . . . . . . . . . . . . . . . . . . 80
5.3.3 Hybrid Dynamical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.4.1 Main Result .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
5.4.2 Choice of Controller Gains . . . . . . . . . . . . . . . . . . . . . . . . . . 89

5.4 Hybrid Control Design .
.

.

viii

5.5

Illustrative Example - The Tiptoebot
5.5.1
5.5.2

Selection of Controller Gains . . . . . . . . . . . . . . . . . . . . . . . .
Simulation Results
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

System Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.6.1
5.6.2
Selection of Controller Gains . . . . . . . . . . . . . . . . . . . . . . . .
5.6.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . . . 90
. 91
. 92
5.6 Experimental Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
. 94
. 95
. 97
5.7 Proofs and Additional Discussions . . . . . . . . . . . . . . . . . . . . . . . . . . 99
. 99
5.7.1
. 100
5.7.2
. 101
5.7.3
5.7.4 Nonoccurrence of Zeno Phenomenon . . . . . . . . . . . . . . . . . . . . 102
. . . . . . . . . . . . . . . . . . . . . 103
5.7.5 Well-posedness of Switching Times
. 103
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.7.6 Quasi-Continuous Dependence Property . . . . . . . . . . . . . . . . . . . 104
. 104
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.7.7 Verification of Assumption 2 for Tiptoebot and Rotary Pendulum . . . . . . 105

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.7.5.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.7.5.2

5.7.6.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.7.6.2

Proof of Lemma 1 .
Proof of Lemma 2 .
Proof of Lemma 3 .

Proof .

Proof .

CHAPTER 6 ORBITAL STABILIZATION USING VIRTUAL HOLONOMIC CON-

. .

.

.

.
.

.

.
.

6.3.1
6.3.2
6.3.3

Introduction .

Problem Statement

6.1
.
6.2 Problem Formulation .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.2.1
6.2.2
6.2.3 Zero Dynamics and Periodic Orbits
6.2.4

STRAINTS AND IMPULSE CONTROLLED POINCARÉ MAPS . . . . . 106
. 106
. 107
System Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
. . . . . . . . . . . . . . 108
Imposing Virtual Holonomic Constraints (VHC)
. . . . . . . . . . . . . . . . . . . . . 109
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110
6.3 Main Result: Stabilization of Od . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Poincaré Map .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
Impulse Controlled Poincaré Map (ICPM) . . . . . . . . . . . . . . . . . . 113
Implementation of Control Design . . . . . . . . . . . . . . . . . . . . . . 116
6.3.3.1 Numerical Computation of A and B matrices
. . . . . . . . . . 116
Impulsive Input using High-Gain Feedback . . . . . . . . . . . . 117
6.3.3.2
. 117
System Dynamics and VHC . . . . . . . . . . . . . . . . . . . . . . . . . 117
Stabilization of VHC and Od . . . . . . . . . . . . . . . . . . . . . . . . . 119
Simulation Results
. 120
. . . . . . . . . . . . . . . . . . . . . . . . . 122
. 122
. . . . . . . . . . . . . . . . . . . . . 123
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

Illustrative Example: Cart-Pendulum . . . . . . . . . . . . . . . . . . . . . . . .
6.4.1
6.4.2
6.4.3
Illustrative Example - The Tiptoebot
6.5.1
6.5.2
6.5.3
6.5.4

System Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
Imposing VHC and Selection of Od
Stabilization of Od
Simulation Results

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6.4

6.5

CHAPTER 7 CONCLUSION AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 127

ix

BIBLIOGRAPHY .

.

.

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132

x

LIST OF TABLES

Table 2.1: Kinematic and Dynamic Parameters for Tiptoebot in SI units . . . . . . . . . . . 16

Table 2.2: Lumped Parameters for Pendubot in SI units . . . . . . . . . . . . . . . . . . . . 30

Table 3.1: Swing-up time and maximum magnitude of high-gain torque required for

different values of N .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

Table 5.1: Tiptoebot lumped parameters in SI units . . . . . . . . . . . . . . . . . . . . . . 91

Table 6.1: Tiptoebot lumped parameters in SI units . . . . . . . . . . . . . . . . . . . . . . 122

xi

LIST OF FIGURES

Figure 1.1: The inertia wheel pendulum.

(a) RA and a point pq˚, 9q˚q that lies outside RA, (b) rRApq˚q and the point
9q˚, (c) IMp 9q˚q, rRApq˚q and pRApq˚, 9q˚q. . . . . . . . . . . . . . . . . . . . . . 12

. . . . . . . . . . . . . . . . . . . . . . . . . .

4

.

(a) The Tiptoebot - a three-link underactuated system with two active joints
at the hip and knee, and one passive joint at the toe (b) a simple model of a
human balancing on the tip of the toe.

. . . . . . . . . . . . . . . . . . . . . . 15

Figure 2.1:

Figure 2.2:

Figure 2.3: Slices of sos and lin, and their intersection with the impulse manifold, IM.

A transformed coordinate system is used for better visualization.

. . . . . . . . 19

coordinates; the intersection set pReA is an ellipse. . . . . . . . . . . . . . . . . . 20

Intersection of the slice of sos with the impulse manifold IM in the original

Figure 2.4:

Figure 2.5: Effective enlargement of Re
A

using SOS-IMM: Plot of joint angles, their
velocities, and control inputs with time for stabilization of the Tiptoebot from
the initial configuration given by (2.3.1). The effect of high-gain feedback
is shown using a dilated time scale; the effect of the stabilizing controller is
shown using a normal time scale.

. . . . . . . . . . . . . . . . . . . . . . . . . 21

Figure 2.6:

A “ 0 at t “ t0 yields
(a) Discretization of the boundary of conservative Re
S0, (b) Convex hull Conv(S1) constructed at time t “ t1, (c) Enlarged Re
obtained as Conv(S1) at time t “ t1`k.

. . . . . . . . . . . . . . . . . . . . . . 25

A

Figure 2.7: Slices of chart, sos and lin, and their intersection with the impulse man-
ifold, IM. For better visualization and for ease of comparison, the same
coordinate system of Fig.2.3 is used.

. . . . . . . . . . . . . . . . . . . . . . . 26

Intersection of the slices of chart and sos with the impulse manifold IM in

the original coordinates; the intersection set pReA of chart is a convex region
and encloses the pReA of sos.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 27

Figure 2.8:

Figure 2.9: Effective enlargement of Re
A

using SOS-CHART-IMM: Plot of joint angles,
their velocities, and control inputs with time for stabilization of the Tiptoebot
from the initial configuration given by (2.3.1).

. . . . . . . . . . . . . . . . . . 28

xii

Figure 2.11: An arbitrary configuration of the two-link pendubot.

of chart with the impulse manifold IM, shown in the original coordinates.

Figure 2.10: (a) Initial configuration of the Tiptoebot, given by (2.5.2), lies outside chart,

(b) rReA of chart (shown in transformed coorrdinates for better visualization)
is a convex hull whereas rReA of sos is a null set, (c) intersection of the slice
Figure 2.12: Intersection of rReA of chart with IM; the intersection set pReA of chart is
configuration p 9q2, 9q1q “ p2.00,0.05q was chosen to lie on pReA of chart. A
zoomed-out scale has been used to show the entire rReA of chart; due to this

the dark line segment. The impulse manifold IM was determined by the
initial velocity configuration p 9q2, 9q1q “ p0.00,0.25q. The desired velocity

scaling the initial and desired velocity configurations appear to be close.

. . . . 32

. . 29

. . . . . . . . . . . . . . . 30

Figure 2.13: Steps of experimental verification using the pendubot shown in Fig.2.11: (a)
1 s - high-gain feedback, (c)

rt0,t´
1 q - feedback linearizing control, (b) rt´
pt`
1 ,t fs - SOS stabilizing control

. . . . . . . . . . . . . . . . . . . . . . . . . 33

1 ,t`

Figure 2.14: Experimental verification of enlargement of Re
A

of pendubot: Plot of joint
angles, their velocities, and control input with time. The effect of high-gain
feedback is shown using a dilated time scale; the effect of the stabilizing
controller is shown using a normal time scale.

. . . . . . . . . . . . . . . . . . 34

Figure 2.15: Simulation results of ideal pendubot behavior for comparison with experi-

mental results presented in Fig.2.14.

. . . . . . . . . . . . . . . . . . . . . . . 34

Figure 3.1: The inertia wheel pendulum is shown in an arbitrary configuration.

. . . . . . . 36

Figure 3.2: An example showing the initial, intermediate, and final configuration of the

IWP for the case with two impulsive inputs (N “ 2).

. . . . . . . . . . . . . . . 42

Figure 3.3: Simulation results for the globally stabilizing controller [1] with controller

parameter values: c0 “ ´π{10, c1 “ 13, c2 “ 16 and c3 “ 8.0.

. . . . . . . . . 51

Figure 3.4: Simulation results for the globally stabilizing controller [1] with controller

parameter values: c0 “ ´π{10, c1 “ 13, c2 “ 16 and c3 “ 4.5.

. . . . . . . . . 52

Figure 3.5: High-gain feedback implementation of two impulsive inputs (N “ 2) for
swing-up of the IWP. For the purpose of comparison with the results in
Figs.3.3 and 3.4, the controller is designed to keep θ in the domain p´π{2,3π{2s. 53

Figure 3.6: High-gain feedback implementation of the sequence of two optimal impulsive
inputs (N “ 2) for swing-up of the IWP. The controller is designed to keep θ
in the domain p´3π{2,π{2s.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 54

xiii

Figure 3.7: High-gain feedback implementation of the sequence of eight optimal impul-
sive inputs (N “ 8) for swing-up of the IWP. The high-gain controller was
implemented with ε “ 0.02. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Figure 3.8: Simulation using the PFBLC + AL controller in [2]; the controller parameters
were chosen as ke “ 3.1ˆ107 and kv “ 0.1. This choice of parameters ensured
that the time required for swing-up and the magnitude of the maximum
control torque in simulations matched those of the experiments. The initial
configuration was chosen to be slightly different from θ0 “ ´π{2 since the
controller is unable to swing-up from this configuration. . . . . . . . . . . . . . 56

Figure 4.1: A three degree-of-freedom of a devil-stick. . . . . . . . . . . . . . . . . . . . . 58

Figure 4.2: Symmetric configurations of the devil-stick in Fig. 4.1.

. . . . . . . . . . . . . 59

Figure 4.3:

(a) Ambidexterous juggler standing at P and applying control actions with
both hands, (b) right-handed juggler standing at P and applying control action
with right hand, (c) right-handed juggler standing at Q and applying control
action with right hand. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. 63

Figure 4.4: State variables and total energy E of the devil-stick at sampling instants k,

k “ 1,2,¨¨¨ ,10, for the LQR design.

. . . . . . . . . . . . . . . . . . . . . . . 73

Figure 4.5: Control inputs for the devil-stick at sampling instants k, k “ 1,2,¨¨¨ ,10, for
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

the LQR design.

. .

. 73

Figure 4.6: State variables and control inputs of the devil-stick at sampling instants k,

k “ 1,2,¨¨¨ ,15, for the MPC design.
(a) Trajectory of the center-of-mass from the initial configuration to steady-
state and (b) symmetric configurations and seven intermediate configurations
of the devil-stick in steady state for the MPC design. . . . . . . . . . . . . . . . 75

. . . . . . . . . . . . . . . . . . . . . . . 75

Figure 4.7:

Figure 5.1: Plots of the joint angles θ1, θ2, θ3, error in the desired energy pE ´ Edesq,

and the active joint velocities 9θ2, 9θ3 of the Tiptoebot.

. . . . . . . . . . . . . . 93

Figure 5.2: Plots showing (a) phase portrait of passive joint angle θ1, and (b) variation of
the Lyapunov-like function V. The desired orbit is shown using dashed green
line in (a).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

.

.

.

.

.

.

Figure 5.3: Schematic of a rotary pendulum.

. . . . . . . . . . . . . . . . . . . . . . . . . 94

xiv

Figure 5.4: Rotary pendulum experimental results: (a)-(d) are plots of joint angles and
joint velocities, (e) control torque, and (f) derivative of Lyapunov-like func-
tion. The brake pulses are shown within plots (e) and (f), the peaks represent
time intervals when the brakes were engaged. . . . . . . . . . . . . . . . . . . . 97

Figure 5.5: Rotary pendulum simulation results.

. . . . . . . . . . . . . . . . . . . . . . . 98

Figure 6.1: Schematic of ICPM approach to orbital stabilization. . . . . . . . . . . . . . . . 115

Figure 6.2:

Inverted pendulum on a cart.

. . . . . . . . . . . . . . . . . . . . . . . . . .

. 118

Figure 6.3: Orbital stabilization for the cart-pendulum system for the initial conditions
in [3]: (a) plot of ρ with respect to time, (b) phase portrait of the pendulum,
(c) plot of the norm of the error of the discrete-time system.

. . . . . . . . . . . 121

Figure 6.4: Orbital stabilization for the cart-pendulum system for the initial conditions in
(6.4.10): (a) plot of ρ with respect to time, (b) phase portrait of the pendulum;
Od is shown in red, (c) plot of the norm of the error of the discrete-time system. 122
Figure 6.5: Phase portrait of tiptoebot zero dynamics. The desired orbit Od is shown in red. 124
Figure 6.6: Orbital stabilization for the tiptoebot using ICPM: (a) and (b) provide plots
of ρ1 and ρ2, (c) and (d) provide plots of 9ρ1 and 9ρ2, (e) and (f) provide the
phase portrait of the passive joint, (g) provides the norm of the error for the
discrete-time system.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126

.

xv

LIST OF ALGORITHMS

CHART
(Convex Hull Algorithm Using Reversal of Trajectories)

. . . . . . . . . . . . . . 24

xvi

KEY TO SYMBOLS

g

(cid:96)1

acceleration due to gravity, (m/s2)

distance of center-of-mass of pendulum from the passive joint, (m)

length of the pendulum, (m)

(cid:96)2
m1 mass of the pendulum, (kg)

m2

t0

ti

t f

xy

I1

I2

N

θ

φ

9θ

9φ

i

θi
9θ´
9θ`
9φ´

i

i

combined mass of the motor and the wheel, (kg)

initial time, (s)
instant of time at which the i-th impulsive torque, is applied, i “ 1,2,¨¨¨ ,N, (s)
final time, (s)

inertially fixed Cartesian coordinate frame

mass moment of inertia of the pendulum about its center of mass, (kg.m2)

mass moment of inertia of the wheel about its center of mass, (kg.m2)

number of impulsive torques applied

angular position of the pendulum, measured CCW with respect to the x axis, (rad)

angular position of the wheel, measured CCW with respect to the pendulum, (rad)

angular velocity of the pendulum, (rad/s)

angular velocity of the wheel, (rad/s)

value of θ at time ti
value of 9θ immediately before application of impulsive torque at time ti
value of 9θ immediately after application of impulsive torque at time ti
value of 9φ immediately before application of impulsive torque at time ti

xvii

9φ`

i

τ

τi

τhg

value of 9φ immediately after application of impulsive torque at time ti
torque applied by the motor on the wheel, (Nm)

impulsive torque applied at time ti, (Nm)
continuous-time implementation of impulsive torque using high-gain feedback, (Nm)

angular impulse of impulsive torque τi applied at time ti, (Nm-s)

Ii
I N the N-dimensional vector rI1 I2 ¨¨¨ INs

xviii

CHAPTER 1

INTRODUCTION

1.1 Motivation

A few decades ago, robots were primarily used for doing repetitive tasks in manufacturing lines

and were supposed to work in risk-free and familiar environments. In recent times, there has been

an increasing interest among engineers and scientists to create robotic systems which can operate

in non-familiar environments. The applications of such robotic systems are immense as they can

be used in application where human health is at risk such as in extra-terrestrial space explorations,

rescue operations in the aftermath of disasters, construction in hostile environments such as deep

sea etc. The new-age robots must be capable of locomotion and as a result, design, control and

stability of robotic systems such as legged robots have gained popularity.

Legged robots fall under the category of underactuated systems as there is no actuation at

the ground-foot contact point. While underactuation is inherent in legged robots, it is purposely

introduced in bio-inspired models such as in the pendubot [4], acrobot [5], tiptoebot [6] etc. Several

control objectives have been studied for underactuated robotic systems such as stabilization about

an equilibrium configuration, swing-up to an equilibrium configuration corresponding to highest

potential energy and generating stable repetitive motion. Since the dynamics of such systems are

typically nonlinear, and due to the presence of fewer control inputs than their degree-of -freedom,

their control design is often very challenging. Due to this reason, it is difficult to design general

control methodologies for underactuated systems.

Impulsive inputs are theoretically modeled as Dirac-delta functions which causes discontinuous

change in the velocity of a dynamical system. Impulsive forces occur naturally in several mechanical

systems that experience impacts; for example walking [7] and hopping robots [8], juggling systems

[9] etc. Modeling of such systems have been widely investigated in the framework of hybrid

dynamical systems [10]. However, the use of impulsive inputs for control of systems, especially

1

underactuated mechanical systems have been rather limited. Few works that have previously

proposed the use of impulsive inputs in control of some underactuated systems can be found

in [11,12]. Impulsive inputs can be extremely useful where sudden change in system’s configuration

is required to attain certain control objective. This work presents several control problem where

impulsive inputs are exploited to simplify the control of underactuated robotic systems; these

problems are discussed next.

1.2 Estimating and Enlarging the Region of Attraction of Equilibria

Underactuated mechanical systems typically have multiple isolated equilibria and often the

control design objective is to stabilize one of the equilibrium points that is unstable. Several methods

have been proposed for stabilizing an equilibrium point for underactuated systems. These include

the controlled Lagrangian method [13, 14], the interconnection and damping assignment passivity-

based control (IDA-PBC) method [15], and the λ-method [16]. Other than these general methods,

control designs have been proposed for stabilizing equilibrium points of specific underactuated

systems such as the pendubot [17–19], acrobot [5,20], inverted pendulum on a cart [21], the reaction-

wheel pendulum [1, 2], and the ball on beam system [22]. Irrespective of the control method used,

a stabilized equilibrium of an underactuated system will have a finite region of attraction since the

majority of all control designs cannot guarantee global stabilization. It is not trivial to obtain a non-

conservative estimate of the region of attraction; and with a conservative estimate, it is only possible

to guarantee stable behavior for small disturbances. To enable underactuated systems remain stable

for large disturbances, we first present a method for obtaining a non-conservative estimate of the

region of attraction. Then, the impulse manifold method (IMM) [23] is used to identify the region

that lies outside this estimate and from where the equilibrium can be stabilized. A portion of this

region may lie outside the region of attraction implying that the IMM can effectively enlarge the

estimate of the region of attraction.

For underactuated systems, the IMM [23] uses impulsive inputs to move the configuration

of the system from outside the estimated region of attraction to inside the region.

Impulsive

2

inputs have been used for control of underactuated systems1 [11, 12, 19, 24] and experimental

investigations [23, 25] have relied on high-gain feedback for continuous approximation of the

impulsive inputs. In an early work [26], high-gain control was used to force the trajectory of an

acrobot into the region of attraction of a stabilizing controller. The limitations of this work are

that the stabilizing controller must be pseudolinear [27] and the system cannot have more than one

input. In [23], impulsive inputs were used to move the configuration of the system from a point

lying immediately outside the region of attraction to inside the region; the boundary of the region

was found iteratively by trial and error. Since it is not possible to determine the region of attraction

for the general case, the IMM is used here in conjunction with an estimate of the region.

Both analytical and numerical methods have been used for estimating regions of attraction,

[28–30], for example. The sum of squares (SOS) method uses a polynomial approximation of the

system dynamics to maximize the estimate of the region of attraction [31–34]. It can also be used

to enlarge the region of attraction by optimizing the controller in polynomial form [35,36]. Starting

from an initial estimate, the trajectory reversing method enlarges an estimate using time-reversed

trajectories of the system, [37–40], for example. The points obtained by time-reversal of trajectories

lie on the surface of a larger estimate but it is challenging to obtain an analytical description of that

surface in the general case.

Both the SOS method and the method of trajectory reversing have been used with the IMM in

our earlier work [41, 42]. We combine the two methods in this work to enhance the applicability

of the IMM. For a given order of the controller polynomial, the SOS method is first used to obtain

the largest estimate of the region of attraction. Using this estimate as the initial estimate, a larger

estimate is obtained using trajectory reversing. The generality of the combined method is illustrated

using the example of the Tiptoebot - a three-link underactuated system with one passive joint. The

result is notable since the trajectory reversing method has been illustrated in the literature using

systems with two states only; the system considered here has six states. In addition to simulations,

experimental validation is provided by stabilizing an equilibrium of a pendubot from a point lying

1Impulsive control of underactuated systems should be differentiated from control of underac-

tuated systems subjected to impulsive disturbances.

3

outside the region of attraction estimated using the SOS and trajectory reversing methods. The

impulse manifold passing through this point does not intersect the region of attraction estimated

using the SOS method alone; this illustrates the usefulness of combining the SOS and trajectory

reversing methods for application of the IMM.

1.3 Rest-to-Rest Maneuver of the Inertia Wheel Pendulum

y

g

c.m.

τ

ℓ2

ℓ1

φ

motor

θ

pendulum

inertia wheel

x

Figure 1.1: The inertia wheel pendulum.

The IWP is comprised of a simple pendulum and a motor mounted at its distal end; the motor

drives a wheel - see Fig.1.1. The torque produced by the motor can be used to accelerate the

wheel in both CW and CCW directions and the resulting reaction torque can be used to control the

pendulum angle. The IWP bears similarity to the underactuated Acrobot [5] where in place of the

the wheel, a link is driven by the motor. Since the wheel is symmetric, the angular displacement of

the wheel is not included in the state-space representation and unlike the Acrobot, the only source

of nonlinearity in the equation of motion is the gravity term. The IWP is an ideal candidate for

impulsive control since the dynamics of the system can be described by simple algebraic equations:

the effect of impulsive forces can be described by changes in the velocities of the system along the

impulse manifold [23], and conservation of energy and conservation of wheel momentum describe

the dynamics when no torque is applied by the motor.

4

For the IWP, stabilization of its upright posture and swing-up to this configuration have been

investigated in the literature, [1, 2, 15, 43–46], for example. Some of the early work on the IWP

can be found in [2] where a hybrid controller comprised of separate swing-up and balancing

controllers was designed using energy methods. Swing-up and stabilization was achieved using

a single controller using the IDA-PBC method [15, 43] and a global change of coordinates that

converted the dynamics of the IWP into strict feedback form [1]. The IDA-PBC method based

controllers [15, 43] and the globally stabilizing controller [1] result in high wheel velocities and

actuator torques.

The swing-up problem of the IWP is a rest-to-rest maneuver between its lowest and highest

potential energy configurations. In this work, we address the general problem of designing rest-to-

rest maneuvers using a sequence of impulsive inputs. The use of impulsive inputs simplifies the

dynamics of the IWP and permits the design of optimal sequences under some general assumptions.

It is found that a sequence with an odd number of inputs is less optimal than the two adjacent

sequences with even number of inputs. The high wheel velocities and/or actuator torques associated

with some results in the literature [1,15,43], is explained using the analysis based on two impulsive

inputs. An optimal sequence of impulsive inputs, implemented using high-gain feedback [23],

is used for swing-up of the IWP; the simulations results show similarities with the energy-based

methods in [2]. It is shown that an optimal input sequence can be designed to accommodate torque

constraint of the actuator; this has not been discussed in earlier works.

1.4 Devil-Stick Juggling

A devil-stick is typically juggled using two hand sticks and several tricks can be performed

depending on the proficiency of the juggler. Some of the common tricks are: standard-idle,

flip-idle, airplane-spin or propeller, top-only idle, and helicopter [47]. The top-only idle is one

of the simplest tricks and is the focus of this investigation; a video tutorial for learning this trick

can be found in [48]. In top-only idle, intermittent forces are applied to the devil-stick. Since the

devil-stick is never grasped, the juggling problem can be viewed as a non-prehensile manipulation

5

problem.

If robots are to perform this trick, the motion of the end-effectors would have to be

coordinated and controlled to apply the correct magnitude of impulsive forces to the devil-stick at

appropriate locations. We do not focus on the end-effector motion control problem (see [49, 50]

for application to ball juggling); instead, we investigate the magnitude and location of the forces

needed to perform the top-only idle trick.

Many juggling tasks, including the top-only idle trick, involve intermittent application of

impulsive forces and several researchers [51–55] have studied the controllability and stability of

such systems. Although impulsive control of the devil-stick has not been investigated, the control

problem associated with juggling of balls and air-hockey pucks has seen several solutions [9,56–58].

In all of these solutions, the object being juggled has been modeled as a point mass and its orientation

is excluded from the dynamic model. In contrast, for devil-stick tricks such as top-only idle, the

stick is shuffled between two symmetric configurations about the vertical; therefore, the orientation

of the stick must be included in the dynamic model.

In earlier works on the devil-stick [59, 60], controllers have been designed for airplane-spin or

propeller-type motion; a single hand-stick is used to rotate the devil-stick about a virtual horizontal

axis using continuous-time inputs. The dynamics model and control design of top-only idle motion

of the devil-stick has not appeared in the literature; to the best of our knowledge, it is presented here

for the first time. It is assumed that impulsive forces are applied intermittently to the devil-stick and

the control inputs are the impulse of the force and its point of application on the stick. The control

inputs are designed to juggle the stick between two symmetric configurations about the vertical,

starting from an arbitrary initial configuration.

1.5 Energy-Based Orbital Stabilization

For underactuated mechanical systems, the problem of stabilization of an equilibrium has been

investigated widely, see [2, 13–16, 18, 21, 61], for example, and the references therein. In many

applications, such as in legged locomotion, the system is required to undergo periodic motion;

therefore, the control objective is to stabilize an orbit rather than an equilibrium. To achieve

6

repetitive motion, the virtual holonomic constraint (VHC) approach [3, 62, 63] has been proposed.

The VHC imposes a set of geometric constraints on the generalized coordinates of the system

and has been demonstrated in bipedal locomotion [7, 64]. Experimental validation of orbital

stabilization using VHC can be found in [65, 66].

Orbital stabilization has also been used for swing-up control of underactuated systems with

one passive degree-of-freedom (DOF). Some examples include two-DOF systems such as the

pendubot [17], the acrobot [67], the reaction-wheel pendulum [68], inverted pendulum on a cart

[69, 70], the rotary pendulum [71], and the three-DOF gymnast robot [72]. These controllers

stabilize an orbit that include the equilibrium, which is typically unstable. The controllers are

designed to gradually pump energy in and out of the system and converge the active DOFs to their

desired configuration. Such control designs are typically based on a Lyapunov-like function that is

comprised of terms involving positions and velocities of the active DOFs and the total mechanical

energy of the system. The passivity property of the system is used to make the time derivative

of the Lyapunov-like function negative semidefinite. A stability analysis is then carried out using

LaSalle’s invariance principle [73]; this typically results in certain restrictions on controller gains

and initial conditions. In all of these prior works, the control strategies for orbital stabilization

have been designed separately for individual systems. Although the structure of the Lyapunov-like

function is identical, the stability analysis is different for each system due to the difference in the

nature of their nonlinear dynamics. Despite the effectiveness of the individual controllers, a general

methodology for energy-based orbital stabilization does not exist. We present a control strategy

for n-DOF underactuated systems with one passive DOF based on continuous time inputs and

impulsive braking.

A combination of continuous and impulsive inputs have been recently used for stabilization

of a homoclinic orbits of two-DOF underactuated systems [74]. The current work provides a

generalization of the theory to n-DOF systems and experimental validation. The main contributions

of this work are as follows:

1. The control design is applicable to a class of underactuated systems which include a majority

7

of the commonly investigated underactuated systems.

2. The stability analysis is presented for the general case and it results in sufficient conditions that

are not restrictive and can be checked for a given system.

3. Experimental validation is provided.

4. Impulsive braking is accomplished using a friction brake; this eliminates the need for high-gain

feedback [23] which may result in actuator saturation.

1.6 Orbital Stabilization using Virtual Holonomic Constraints

For underactuated systems, Virtual Holonomic Constraint (VHC) based control designs have

gained popularity due to their conceptual simplicity and applicability to control of repetitive

motion. They have been used for stabilization of gaits in bipeds [64, 75–80], which undergo

impacts, and trajectory control for systems with open kinematic chains that are not subjected to

impacts [3, 62, 63, 65, 66, 81–84]. VHCs parameterize the active joint variables in terms of the

passive joint variables and confine system trajectories to a constraint manifold [62]. To enforce the

VHC, the constraint manifold has to be stabilized using feedback. Typically, a constraint manifold

contains a dense set of periodic orbits and the choice of repetitive motion determines the specific

orbit that has to be stabilized. To stabilize gaits in bipeds, for example, Grizzle et.al [64, 76, 79]

enforced the VHC and periodic loss of energy due to ground-foot interaction was exploited for

orbital stabilization.

A special class of underactuated systems are those with open kinematic chains and one passive

DOF. For such systems, Shiriaev and collaborators [3, 65, 66] used VHC to select the desired orbit.

For an n-DOF system, the 2n dimensional dynamics is linearized about the desired orbit; this results

in a 2n´1 dimensional system. A periodic Ricatti equation is then solved to design a time-varying
controller that stabilizes the orbit. It should be noted that the control designs in [3, 65, 66] stabilize

the orbit but do not enforce the VHC. A control scheme that enforces the VHC and simultaneously

stabilizes the orbit was recently proposed in [63]. The key idea is that the VHC is made time-

8

varying using a scalar parameter which is controlled via feedback. The stabilization problem

involves solving a periodic Ricatti equation; however, unlike [3,65,66], where the dimension of the

system is 2n´1, the dimension of the system in [63] is always three. For systems with more than
two DOF, the method in [63] reduces the computational complexity of control implementation.

Also, by enforcing the VHC, it improves control over transient characteristics of the trajectory [85].

Similar to [63, 85], we design a continuous controller with the objective of enforcing the VHC and

stabilizing the constraint manifold. We then design impulsive inputs that work in tandem with the

continuous inputs to stabilize a desired periodic orbit on the constraint manifold. Compared to the

methods proposed earlier [3, 63, 65, 66], where periodic Ricatti equations have to be solved, our

method has lower computational cost and complexity. Since impulsive inputs are used to control the

Poincaré map, the dynamics of the closed-loop system can be described by the Impulse Controlled

Poincaré Map (ICPM). The simplicity and generality of the ICPM approach to orbital stabilization

is demonstrated using the examples of the 2-DOF cart-pendulum and the 3-DOF tiptoebot.

9

CHAPTER 2

ESTIMATION OF THE REGION OF ATTRACTION OF EQUILIBRIA AND ITS

ENLARGEMENT USING IMPULSIVE INPUTS

2.1 Introduction

Stabilization of an equilibrium point is an important control problem for underactuated systems.

For a given stabilizing controller, the ability of the system to remain stable in the presence of

disturbances depends on the size of the region of attraction of the stabilized equilibrium. The

sum of squares (SOS) and trajectory reversing methods are combined together to generate a large

estimate of the region of attraction. Then, this estimate is effectively enlarged by applying the

impulse manifold method (IMM), which can stabilize equilibria from points lying outside the

estimated region of attraction. The IMM and the SOS method are reviewed in section 2.2. The

Tiptoebot is introduced in section 2.3 and the IMM is used to stabilize its equilibrium from a point

lying outside the region of attraction estimated by the SOS method. A convex hull algorithm for

estimating the region of attraction using trajectory reversing is presented in section 2.4. The benefit

of using the SOS method, the convex hull algorithm, and IMM together is illustrated using the

example of the Tiptoebot in section 2.5. Experimental results are presented in section 2.6.

2.2 Background

2.2.1 The Impulse Manifold Method

For an n degree-of-freedom dynamical system, Lagrange’s equations have the form:

Mpqq :q` Hpq, 9qq “ T

(2.2.1)

where q P Rn and T P Rn represent the vectors of generalized coordinates and generalized forces;
Mpqq P Rnˆn denotes the generalized mass matrix; and Hpq, 9qq P Rn denotes the vector of centrifu-
gal, Coriolis, gravitational and friction forces. If the system is underactuated with m active and

10

pn´ mq passive degrees-of-freedom, m ă n, (2.2.1) can be rewritten in the form:

»—–

M11pqq M12pqq
MT
12pqq M22pqq

fiffifl

»—–

:q1

:q2

fiffifl`»—–

h1pq, 9qq
h2pq, 9qq

fiffifl “»—–

0

u

fiffifl

(2.2.2)

1 , qT

2 qT , and q1 P Rpn´mq and q2 P Rm denote the generalized coordinates correspond-
where q “ pqT
ing to the passive and active degrees-of-freedom, respectively, and u P Rm is the vector of control
inputs.

a stabilizing controller.

We assume that pq, 9qq “ p0, 0q is an equilibrium point of the unforced system and u “ us is
If RA is the region of attraction of the equilibrium with u “ us, the
equilibrium cannot be stabilized from configurations that lie outside RA. However, it may be

possible to stabilize the equilibrium after moving the configuration of the system from outside RA

to inside RA using the Impulse Manifold Method (IMM) [23]. In the IMM, impulsive inputs are

applied to the active degrees-of-freedom. Ideal impulsive inputs cause no change in the generalized

coordinates but result in discontinuous changes in both active and passive generalized velocities.

These changes can be computed from the following equation which is obtained by integration of

(2.2.2) over the infinitesimal interval of time Δt, during which the impulsive inputs act [86]:

»—–

M11 M12

MT

12 M22

fiffifl

»—–

Δ 9q1

Δ 9q2

fiffifl “»—–

0

I

fiffifl ,

I fiż Δt

0

u dt

(2.2.3)

In the above equation, I P Rm is the vector of impulses, and the velocity jumps, Δ 9q1 and Δ 9q2, are
defined as

Δ 9q1 fi p 9q`

1 ´ 9q´
1 q,

Δ 9q2 fi p 9q`

2 ´ 9q´
2 q

where 9q´ and 9q` are the generalized velocities immediately before and after application of the

impulsive inputs. Due to the underactuated nature of the system, the discontinuous changes in

the velocities of the passive coordinates are dependent on those of the active coordinates; this

dependence is obtained from (2.2.3) as follows:

p 9q`
1 ´ 9q´

1 q “ ´M´1

11 M12p 9q`

2 ´ 9q´
2 q

(2.2.4)

11

In the n-dimensional velocity space, the impulse manifold is the m-dimensional manifold repre-

sented by (2.2.4). It was stated earlier that the generalized coordinates do not change under the

application of an impulse. Hence, for a fixed value of q, the impulse manifold passing through the
point p 9q´

2 q is defined by the set

1 , 9q´

IMp 9q´

1 , 9q´

2 q “ t 9q`

1 , 9q`

2 | p 9q`

1 ´ 9q´

1 q “ ´M´1

11 M12p 9q`

2 ´ 9q´
2 qu

(2.2.5)

The IMM is explained with the help of Fig.2.1. The region of attraction RA is shown in Fig.2.1
(a). An arbitrary point outside RA is denoted by pq˚, 9q˚q. For q “ q˚, the slice of RA in the
n-dimensional velocity space is denoted by rRApq˚q and is shown in Fig.2.1 (b). Since pq˚, 9q˚q lies
outside RA, the point 9q˚ lies outside rRApq˚q. The m-dimensional impulse manifold IMp 9q˚q is
depicted by the plane in Fig.2.1 (c); its intersection with rRApq˚q is shown by the hatched region

and is mathematically defined by the set

Impulsive inputs can be applied to change the velocity of the system from 9q˚ to any other point

pRApq˚, 9q˚q “ IMp 9q˚qX rRApq˚q

(2.2.6)

(desired velocity) on the impulse manifold. However, to move the configuration of the system to

a point inside RA, the desired velocity should be chosen to lie in pRApq˚, 9q˚q. For the IMM to be
applicable, pRApq˚, 9q˚q should not be a null set and the desired velocity should be chosen such that

the magnitude of the impulse I does not violate actuator constraints.

IMp 9q˚q

RA Ă R

2n

pq˚, 9q˚q

9q˚

9q˚

rRApq˚q Ă Rn

(a)

(b)

pRApq˚, 9q˚q

(c)

Figure 2.1: (a) RA and a point pq˚, 9q˚q that lies outside RA, (b) rRApq˚q and the point 9q˚, (c)
IMp 9q˚q, rRApq˚q and pRApq˚, 9q˚q.

12

Ideal impulsive inputs, which are Dirac-delta functions, cannot be implemented in real physical

systems. For an underactuated system of the form as in (2.2.2), impulsive inputs can be implemented

using high-gain feedback [23]; the control law is given by1

u “ uhg “ rCT M´1

Cs´1„CT M´1
where C P Rnˆm is the matrix C “ r0 EsT , E P Rmˆm is the identity matrix, 9qdes
desired velocities of the active joints, Λ fi diag„λ1 λ2

2 ´ 9q2¯
Λ´ 9qdes
¨¨¨ λm, where λi, i “ 1, 2,¨¨¨ , m are

is the vector of

(2.2.7)

H`

2

1

µ

positive numbers, and µ is a small positive number.

In general, an exact determination of RA is not possible and an estimate of the region of

attraction Re

A, Re

A Ă RA, has to be used. As we replace RA with Re

same manner that rRA and pRA were defined in relation to RA.

A, we define rRe

A and pRe

A in the

From the discussion in this section it is obvious that the IMM will be more effective for larger

Re

A. When applicable, the IMM will effectively enlarge the Re

A of the equilibrium.

2.2.2 The Sum of Squares Method

The simplest way to obtain an Re

A is through linearization [73]: the dynamics of the closed-loop

underactuated system in (2.2.2) is first represented in state-space form

9x “ frx, u “ uspxqs “ gpxq

(2.2.8)

where x “ pqT , 9qTqT and uspxq is a stabilizing controller. The Re
A of the asymptotically stable
equilibrium x “ 0 can be defined with the help of the quadratic Lyapunov function V0pxq as follows
(2.2.9)

Re
A “ tx | V0pxq ď cu fi Ωlin,

V0pxq fi xT Px

where P is the positive definite matrix obtained by solving the Lyapunov equation

PA` AT P “ ´Q,

A fi„Bg

Bxx“0

, Q “ QT ą 0

1The high-gain feedback in [23] has been modified here for its applicability to systems with
more than one active joint:
the diagonal matrix Λ has been introduced to provide flexibility in
choosing different gains for the different active joints. It can be easily shown that this modification
does not affect the continuous approximation of the impulsive inputs.

13

and c is the largest number that satisfies

tV0pxq ď cu Ď t 9V0pxq ď 0u

The Sum of Squares (SOS) method [31] automates the process of finding a Lyapunov function

using numerical tools from convex optimization. Subsequently, an Re

A can be obtained by solving

the following optimization problem [36]:

subject to

maximize ρ

Vpxq
Lpxq
´ 9Vpxq` LpxqpVpxq´ ρq

is SOS

is SOS

is SOS

$’’’’’’&’’’’’’%

where Vpxq is the Lyapunov function and Lpxq is a polynomial multiplier. The decision variables
in this optimization problem are Lpxq, Vpxq and ρ. The third constraint guarantees

tVpxq ď ρu Ď t 9Vpxq ď 0u

and therefore an estimate of the region of attraction is

Re
A “ tx | Vpxq ď ρu fi Ωsos

(2.2.10)

In the above optimization problem, Re

The controller can also be optimized to further enlarge Re

A is maximized for a pre-determined stabilizing controller.
A. The controller upxq is introduced as

another decision variable in the third constraint of the SOS optimization problem:

9Vpxq “ BV
Bx

frx, upxqs

This non-convex optimization problem was studied in [35, 36], where the following steps are

performed sequentially and repeatedly till no further enlargement of Re

A is observed:

• Vpxq is fixed; Lpxq and upxq are varied to maximize ρ

• Lpxq and upxq are fixed; Vpxq is varied to maximize ρ

14

The SOS algorithm requires an initial guess of the Lyapunov function Vpxq. If this initial guess

is chosen to be V0pxq in (2.2.9), we can guarantee

Ωsos Ě Ωlin

which implies that the SOS method can improve upon the Re

A obtained through linearization.

2.3 Enlarging the Region of Attraction Using SOS and IMM

2.3.1 An Illustrative Example

From our discussion in section 2.2 it is clear that an Re

A can be obtained using the SOS method
and that this estimate can be effectively enlarged using the IMM. In this section we demonstrate the

effective enlargement of the Re

A for a three-link underactuated system with two active joints and
one passive joint. A majority of articulated underactuated systems that have been considered in the

literature have two degrees-of-freedom with one passive joint; examples include the pendubot, the

g

ℓ3

θ3

τ3

pm2, J2q

d2

y

ℓ1ℓ1

d1

pm3, J3q

ℓ3

upper body

upper leg

ℓ2

d3

ℓ2

θ2

τ2

pm1, J1q

ℓ1

lower leg

θ1

x

(a)

(b)

Figure 2.2: (a) The Tiptoebot - a three-link underactuated system with two active joints at the hip
and knee, and one passive joint at the toe (b) a simple model of a human balancing on the tip of the
toe.

15

acrobot, the rotary pendulum, inverted pendulum on a cart, etc. The three-link system is purposely

chosen to demonstrate the applicability of the method to more complex systems.

Consider the three-link underactuated system2 in Fig 2.2 (a). We have christened it the Tiptoebot

since it models a human balancing on the tip of its toe - see Fig 2.2 (b). The three links are analogous

to the lower leg, the upper leg, and the upper body comprised of the torso and head. The knee

joint connecting the upper and lower legs, and the hip joint connecting the upper body and upper

leg are actuated; the torques applied by the actuators in these joints are assumed to be positive in

the counter-clockwise direction and are denoted by τ2 and τ3. The toe provides a simple point

of support and is modeled as a passive joint. The lower leg, upper leg, and upper body have

link lengths ℓ1, ℓ2 and ℓ3; their center-of-masses are located at distances d1, d2 and d3 as shown

in Fig.2.2 (a); their mass and mass moments of inertia about their individual center-of-mass are

denoted by pm1, J1q, pm2, J2q and pm3, J3q. The joint angles of the links are measured positive in
the counter-clockwise direction and are denoted by θ1, θ2 and θ3; θ1 is measured relative to the

x-axis whereas θ2 and θ3 are measured relative to the second and third links.

Table 2.1: Kinematic and Dynamic Parameters for Tiptoebot in SI units

m1

0.20

J1

m2

0.20

J2

m3

ℓ1

ℓ2

ℓ3

0.30

0.20

0.20

0.30

J3

d1

d2

d3

6.66ˆ 10´4

6.66ˆ 10´4

15.0ˆ 10´4

0.10

0.10

0.15

The kinematic and dynamic parameters of the Tiptoebot in Fig.2.2 are provided in Table 6.1.

For these set of parameters, the equations of motion are of the form given by (2.2.2), where n “ 3,
m “ 2 and u “ pτ2 τ3qT . The vertically upright posture or configuration pθ1,θ2,θ3, 9θ1, 9θ2, 9θ3q “
pπ{2, 0, 0, 0, 0, 0q is the desired equilibrium and hence q1 “ pθ1 ´ π{2q, q2 “ pθ2 θ3qT . The
matrices M11, M12, M22, h1 and h2, which can be described as functions of θk, 9θk, k “ 1, 2, 3, are

2This system is a special case of the N-link system considered in [87].

16

provided below:

M11 “ 0.018C23 ` 0.032C2 ` 0.018C3 ` 0.046

M12 “»—–
M22 “»—–

0.009C23 ` 0.016C2 ` 0.018C3 ` 0.023

0.009C23 ` 0.018C3 ` 0.009

0.018C3 ` 0.023
0.009C3 ` 0.009

0.009C3 ` 0.009

0.009

T

fiffifl

fiffifl

h1 “ 0.100 9θ1 ` 0.441C123 ` 0.784C12 ` 1.177C1 ´ 0.016 9θ2
´ 0.009p 9θ2 ` 9θ3q2
h2 fi„h21 h22T
h21 “ 0.100 9θ2 ` 0.441C123 ` 0.784C12 ` 0.016 9θ2

1 S2 ´ 0.009 9θ2

2 S2 ´ 0.009 9θ2

3 S3

S23 ´ 0.018 9θ1p 9θ2 ` 9θ3q S23 ´ 0.018p 9θ1 ` 9θ2q 9θ3 S3 ´ 0.032 9θ1

9θ2 S2

3 S3 ` 0.009 9θ2

1 S23

´ 0.018p 9θ1 ` 9θ2q 9θ3 S3

h22 “ 0.100 9θ3 ` 0.441C123 ` 0.009 9θ2

1 S23 ` 0.009p 9θ1 ` 9θ2q2

S3

where Si and Ci denote sinθi and cosθi, Si j and Ci j denote sinpθi`θjq and cospθi`θjq, and Ci jk
denotes cospθi`θj`θkq. We assumed the presence of viscous friction in the joints of the Tiptoebot.
This accounts for the terms linearly proportional to the angular velocities in the expressions for h1,

h21 and h22.

The stabilizing controller us in (2.2.8) was obtained through linearization by placing the poles
of the closed-loop system at -8, -3, -2, -1, -0.7, -0.4. The controller has the form us “ rτ2 τ3sT
where

τ2 “ ´10.204pθ1 ´ π{2q´ 5.789θ2 ´ 2.140θ3 ´ 2.114 9θ1 ´ 1.144 9θ2 ´ 0.496 9θ3
τ3 “ ´7.996pθ1 ´ π{2q´ 4.284θ2 ´ 1.878θ3 ´ 1.632 9θ1 ´ 0.974 9θ2 ´ 0.295 9θ3

17

»——————————————–

fiffiffiffiffiffiffiffiffiffiffiffiffiffiffifl

The Lyapunov function V0 was obtained using (2.2.9); the matrices P and Q are

5271.9 2668.5 973.5 1054.4 626.9 241.7

2668.5 1350.7 492.7 533.75 317.3 122.3

P “

973.55 492.80 179.8 194.73 115.7 44.65

1054.4 533.75 194.7 210.91 125.3 48.35

, Q “ 0.01ˆ diag„3 3 3 1 1 1

626.90 317.33 115.7 125.39 74.55 28.75

241.75 122.37 44.65 48.359 28.75 11.09

The estimated region of attraction Ωlin in (2.2.9) was defined by the value of c “ 0.020; this value
was determined by following the steps [73] given below:

1. Find a set of points αi, i “ 1, 2,¨¨¨ , M, where 9Vpαiq “ 0.

2. Choose c ă ¯c “ min

i

Vpαiq.

The SOS method was used to improve upon the estimate Ωlin. By fixing V “ V0, L and u were
varied to maximize ρ and determine Ωsos in (2.2.10). The input u was chosen to be a first-order

polynomial of the states and the maximum value of ρ was found to be 0.044. The SOS controller

was obtained as follows:

τ2 “ ´3.866pθ1 ´ π{2q´ 2.603θ2 ´ 1.025θ3 ´ 0.898 9θ1 ´ 0.417 9θ2 ´ 0.259 9θ3
τ3 “ 10.157pθ1 ´ π{2q` 4.887θ2 ` 1.410θ3 ` 1.942 9θ1 ` 1.156 9θ2 ` 0.476 9θ3

2.3.2 Simulation Results

To demonstrate the enlargement of the Re

A “ Ωsos using the IMM, we consider the following initial

configuration that lies outside Ωsos:

x “„q˚T 9q˚TT

“„0.0 0.0 0.0 0.00 0.00 0.21T

(2.3.1)

For this configuration, it can be verified that V0pxq “ 0.48 ą ρ “ 0.044. Following the discussion
Apq˚ “ 0q corresponding to Ωlin and Ωsos can be expressed as

in section 2.2.1, the slices rRe

18

o f Ω

lin

A

rRe

IM

o f Ωsos

A

rRe

Figure 2.3: Slices of Ωsos and Ωlin, and their intersection with the impulse manifold, IM. A
transformed coordinate system is used for better visualization.

V0 |q“q˚“0ď c and V0 |q“q˚“0ď ρ, respectively. Both these slices are ellipsoids and are shown in
Fig.2.3; the coordinates were scaled and transformed for better visualization. The impulse manifold

IMp0, 0, 0.21q is a plane (two-dimensional manifold) in the three-dimensional velocity space; it is
shown in Fig.2.3 in the transformed coordinates and mathematically described by the relation:

9q1 ` 0.583 9q2 ` 0.236 9q3 ´ 0.049 “ 0

(2.3.2)

Apq˚ “ 0q “ t 9q |
Apq˚ “ 0q “ t 9q | V0p0, 9qq ă cu
A in (2.2.6) obtained using the SOS method
is non-empty whereas that obtained using linearization is a null set; this illustrates the benefit of

It can be seen from Fig.2.3 that the impulse manifold intersects the slice rRe
V0p0, 9qq ă ρu obtained using the SOS method but not the slice rRe
obtained using linearization. Consequently, the set pRe

using a larger Re

A in terms of applicability of the IMM.

A obtained using the SOS method is an ellipse in the three-dimensional velocity space;
it is shown together with the impulse manifold IM in Fig.2.4 in the original coordinates. To take

The pRe

the configuration of the Tiptoebot from outside Re

A, we choose the following desired

velocity configuration that lies in the interior of pRe

„ 9q1 9qT

2 “„2.36 ´4.04 0.16

A to inside Re
A:

(2.3.3)

19

3.0

2.0

1.0

9q3

0.0

-1.0

A

pRe

p2.36, ´4.04, 0.16q

2.0

0.0

9q1

p0, 0, 0.21q

IM

0.0

-2.0

9q2

-4.0

Figure 2.4: Intersection of the slice of Ωsos with the impulse manifold IM in the original coordi-

nates; the intersection set pRe

A is an ellipse.

The desired effect of the impulsive inputs is shown by the dotted line joining the initial and desired

gain feedback in (2.2.7) was used with 9qdes

velocity configurations, both of which lie on IM - see Fig.2.4. To implement the IMM, the high-
2 “ 9q2 in (2.3.3), µ “ 0.0004, λ1 “ 10 and λ2 “ 1; the
results are shown in Fig.2.5. Figure 2.5 plots the joint angles, their velocities, and the control inputs

with time. The high-gain controller is active over t P r0, 0.004s; its effect is shown using a dilated
time scale. It can be seen that the change in the joint angles are negligible but there are significant
changes in the joint velocities. At t “ 0.004 when 9q2 « 9qdes
off and the stabilizing controller is switched on; at this time, the system configuration is inside Re
A

2 , the high-gain controller is switched

and therefore the equilibrium can be stabilized. The stabilization of the equilibrium after t ą 0.004
is shown using a normal time scale.

The high-gain feedback was based on µ “ 0.0004. The change in the joint angles were
negligible, of the order of 0.01 rad. A larger value of µ could be used to reduce the magnitude of

the control inputs but it would result in larger changes in the joint angles. The computation of pRe

is based on the assumption that there is no change in the joint angles and therefore a small value of

A

µ is necessary to ensure the applicability of the IMM.

For the general case, if δ denotes the small change in the joint angles due to high-gain feedback,

20

θ2

θ1

9θ3

2.0

1.0

π{2

0.0

3.0

0.0

-4.0

0.0

9θ1

9θ2

-1.5

τ2

-3.0

0.0

-2.0

-4.0

τ3

dilated time scale

normal time scale

(rad)

θ3

(rad/s)

(Nm)

(Nm)

0.0

0.004

time (s)

1.5

3.0

Figure 2.5: Effective enlargement of Re
A using SOS-IMM: Plot of joint angles, their velocities,
and control inputs with time for stabilization of the Tiptoebot from the initial configuration given
by (2.3.1). The effect of high-gain feedback is shown using a dilated time scale; the effect of the
stabilizing controller is shown using a normal time scale.

we should ensure that the desired velocity 9qdes

2 , which was chosen to lie in the interior of pRApq˚, 9q˚q,

should be chosen

also lies in the interior of pRApq˚ ` δ, 9q˚q. To this end, the desired velocity 9qdes
sufficiently far from the boundary of pRApq˚, 9q˚q.

There is significant flexibility in the choice of the desired velocity 9qdes

2

2 . This choice can be

used judiciously to minimize the overall control effort or some other cost function.

It was mentioned in section 2.2.1 that the IMM will be more effective for larger Re

end, we present an algorithm in the next section for obtaining large estimates of Re

A. To this
A; the algorithm

21

is based on the method of trajectory reversing [37].

2.4 Algorithm for Computing Re

A Using the Method of

Trajectory Reversing

2.4.1 Definitions

For the closed-loop system in (2.2.8), an Re

(2.2.10). The method of trajectory reversing [37] can be used to obtain a larger Re

A of x “ 0 is given by Ωlin in (2.2.9) and Ωsos in
A starting from a
A is a convergent sequence of simply connected
domains generated by backward integration from a discrete set of points on the boundary of the

conservative Re

A such as Ωlin or Ωsos. The larger Re

conservative estimate. In the literature, the method of trajectory reversing has been used to obtain

Re

A for example systems described by two state variables. Here, we present an algorithm [42] that
has been shown to provide larger estimates for the two-state systems considered in the literature;

additionally, it will be used to obtain an Re

A for the Tiptoebot, whose dynamics is described six state
variables. To the best of our knowledge, estimation of the region of attraction using the method of

trajectory reversing has not been reported in the literature for systems with more than two states.

Before we present the algorithm, we present the following definitions:

Definition 1. Convex Hull: A convex hull of a set of points S “ tx1, x2,¨¨¨ , xNu in Rn, denoted by

Conv(S), is the smallest convex set that contains all the points.

A convex hull can be mathematically represented by a set of equations and inequalities. In the

two-dimensional plane, a convex hull is described by line segments which can be represented by

equations of the lines and inequalities based on their lengths. For higher dimensions, the convex

hull is represented by planes/hyper-planes and inequalities. The Quick Hull Algorithm [88] in

Matlab can be used to compute the convex hull of a set of points efficiently. For a given convex

hull, the Matlab function inhull [89] can be used to determine whether an arbitrary point lies in the

interior of the hull or not.

Definition 2. Slice of a Convex Hull: The slice of a convex hull is the set obtained by intersecting

22

the convex hull with a hyper-plane.

The slice of a convex hull is also convex.

It can be computed using the Matlab function

intersectionhull [90].

Definition 3. Cluster: For a given set of points, a cluster is a subset where any two points in the

subset is density-reachable for some pre-defined distance metric ε.

Simply stated, two points are said to be density-reachable if there is a chain of intermediate

points and the distance between neighboring points does not exceed ε. A more formal definition

of cluster and a cluster identification algorithm can be found in [91]. A cluster of points can be

identified using the Matlab function DBSCAN [92].

2.4.2 Algorithm

Nu on the boundary of the conservative Re

We choose a discrete set of points S0 “ tx0
time t0, where N is a sufficiently large number. The conservative Re

A at
A, denoted by Ω0, can be Ωlin
or Ωsos, for example. Now consider the time sequence ti, i “ 0, 1, 2,¨¨¨ , where tk “ t0 ` kΔt, and
Δt is some pre-defined time interval. Define the set Si “ txi
Nu, i “ 0, 1, 2,¨¨¨ , which is
obtained by reversing the trajectories of the points in S0 for iΔt interval of time. Conv(Si) is an

2,¨¨¨ , xi

1, x0

2,¨¨¨ , x0

1, xi

enlarged Re

A if:

1. none of the points in Si`k, where k is some pre-defined positive integer, lies in the interior of

Conv(Si), and

2. the vertices of Conv(Si) form a single cluster around the origin with distance metric ε [42].

For i, j “ 0, 1, 2,¨¨¨ , let m
a time ti` j and NC(Conv(Si)) represent the number of clusters formed by the vertices of Conv(Si).

denote the number of points which remain in the interior of Conv(Si) at

i` j
i

The complete algorithm for obtaining a large estimate of Re

A, which we refer to as CHART (Convex

Hull Algorithm Using Reversal of Trajectories), is presented next.

23

Algorithm CHART
(Convex Hull Algorithm Using Reversal of Trajectories)

Initialization: Choose a sufficiently large integer N which represents the number of points
located uniformly on the boundary of Ω0 and a sufficiently small distance metric ε.
Step 1: Select a suitable time interval Δt and sufficiently large positive integers k1 and k2,
k1 ă k2.
Step 2:
Let ℓ “ 0
For i “ 1 to k1
If ℓ “ i´ 1
If m

For j “ 1 to k2

i` j

i “ 0 and NC(Conv(Si)) “ 1
ℓ “ i
Break (exit from the inner For loop which

checks the condition on j)

End If
End For

End If
End For
Step 3: If ℓ “ 0, repeat Steps 1 to 2 with a lower value of Δt than chosen in Step 1
Result: ConvpSℓq is the enlarged Re
A.
Parameter Selection: The choice of k1 and k2 is a trade-off between the enlargement of the Re
A
and the computation time. A few examples showing the choice of N, ε, Δt, k1, k2 are provided
in [42].

For the ease of understanding, the CHART is explained for a system with two state variables

2,¨¨¨ , x0

with the help of Fig.2.6. At time t “ t0, the boundary of the conservative estimate Ω0 is discretized
into the set of points S0 “ tx0
Nu, where N is a large number - see Fig.2.6 (a). At time
t “ t1, it is found that some points in S1 lie in the interior of Conv(S1), i.e. m1
1 ­“ 0 - see Fig.2.6
(b). However, after k intervals of time Δt, it is found that none of the points in S1`k lie in the
1 “ 0. Additionally, the vertices of Conv(S1) form a single cluster, i.e.
A “ ConvpS1q - see Fig.2.6 (c). The single cluster requirement
ensures that the vertices of the convex hull are sufficiently close to one another and hence the

NC(Conv(S1)) = 1. Therefore, Re

interior of Conv(S1), i.e. m

1, x0

1`k

boundary of the convex hull is a positively invariant set.

24

Conv(S1)

1
x
2

1
x
1

1
x
N

1`k
x
2

1`k
x
3

1`k
x
1

1`k
x
N

Re
A

0
x
1

0
x
N

1
x
3

p0, 0q

p0, 0q

p0, 0q

0
x
2

0
x
3

0
x
i

Ω0

t “ t0

S0 “ tx

0

N u

0
1, ¨ ¨ ¨ , x
(a)

Conv(S1`k)

t “ t1`k “ pt1 ` kΔtq
1`k
N u

1`k
S1`k “ tx
1

, ¨ ¨ ¨ , x

t “ t1 “ pt0 ` Δtq
N u

S1 “ tx

1

1
1, ¨ ¨ ¨ , x
(b)

Figure 2.6: (a) Discretization of the boundary of conservative Re
Convex hull Conv(S1) constructed at time t “ t1, (c) Enlarged Re
t “ t1`k.

(c)

A “ Ω0 at t “ t0 yields S0, (b)
A obtained as Conv(S1) at time

2.5 Enlarging the Region of Attraction Using SOS, CHART and IMM

Re

In the last section, we presented CHART for obtaining a larger Re

A starting from a conservative
A for the
Tiptoebot starting from the conservative estimate Ωsos. It should be noted that the SOS method

A. To improve the effectiveness of the IMM, we now apply CHART to obtain an Re

optimizes the stabilizing controller to obtain the largest Re

A and the CHART improves upon this

estimate further without changing the controller.

The CHART was implemented for the Tiptoebot equilibrium configurationpθ1,θ2,θ3, 9θ1, 9θ2, 9θ3q“

pπ{2, 0, 0, 0, 0, 0q using the following parameter values, which were selected by trial and error:

N « 3ˆ 10

6, ε “ 0.01, Δt “ 0.0026,

k1 “ 10,

k2 “ 65

The enlarged estimate Re

A, which is a convex hull in a six-dimensional space, was obtained as

Ωchart “ ConvpS4q.

25

o f Ωchart

A

rRe

o f Ωsos

A

rRe

A of Ωlin

pRe

IM

Figure 2.7: Slices of Ωchart, Ωsos and Ωlin, and their intersection with the impulse manifold, IM.
For better visualization and for ease of comparison, the same coordinate system of Fig.2.3 is used.

given by (2.3.1), the slice rRe

We now reconsider the simulation in section 2.3.2. For the initial configuration of the Tiptoebot
Apq˚ “ 0q of Ωchart is shown in Fig.2.7 together with the slices
of Ωlin and Ωsos that were previously shown in Fig.2.3.
A of Ωchart is
significantly larger than that of Ωsos. The impulse manifold IMp0, 0, 0.21q described by (2.3.2),

It can be seen that rRe

A of both Ωchart and Ωsos in Fig.2.7.

shown previously in Fig.2.3, can be seen to intersect rRe
The pRe
original coordinates. The pRe
in the interior of pRe

A of Ωchart is a convex region on the impulse manifold and is shown in Fig.2.8 in the
A of Ωsos, shown previously in Fig.2.4 and now shown in Fig.2.8, lies
A to
A of

A, we now choose the following velocity configuration that lies in the interior of pRe

A of Ωchart. To take the configuration of the Tiptoebot from outside Re

inside Re

Ωchart:

„ 9q1 9qT

2 “„1.57 ´2.69 0.18

(2.5.1)

The desired effect of the impulsive inputs is shown by the dotted line joining the initial and desired

velocity configurations, both of which lie on IM - see Fig.2.8. To implement the IMM, the high-
2 “ 9q2 in (2.5.1), µ “ 0.0004, λ1 “ 10 and λ2 “ 1; the

gain feedback in (2.2.7) was used with 9qdes

26

A of Ωchart

pRe

3.0

2.0

1.0

9q3

0.0

-1.0

A of Ωsos

pRe

p2.36, ´4.04, 0.16q

p1.57, ´2.69, 0.18q

2.0

0.0

9q1

p0, 0, 0.21q

IM

0.0

-2.0

9q2

-4.0

Figure 2.8: Intersection of the slices of Ωchart and Ωsos with the impulse manifold IM in the
A of

A of Ωchart is a convex region and encloses the pRe

original coordinates; the intersection set pRe

Ωsos.

results are shown in Fig.2.9. The high-gain controller is active over t P r0, 0.004s; its effect is shown
using a dilated time scale. The high-gain controller is switched off and the stabilizing controller
is switched on at t “ 0.004 when the system configuration is inside Re
equilibrium after t ą 0.004 is shown using a normal time scale.

A. The stabilization of the

The desired velocity configuration in (2.5.1) was chosen to lie on the line joining the initial

velocity configuration p0.0, 0.0, 0.21q and the desired velocity configuration chosen for the SOS-
IMM method in section 2.3, given by (2.3.3). It can be seen from Fig.2.8 that the initial velocity

configuration is closer to the desired velocity configuration in (2.5.1) than the desired velocity

configuration in (2.3.3). Consequently, the magnitudes of the high-gain feedback for the SOS-

CHART-IMM method are less than that of the SOS-IMM method - see plots of τ2 and τ3 in

Figs.2.5 and 2.9 for t P r0.00, 0.004s.

We complete this section with a second example that further illustrates the benefit of using

SOS-CHART-IMM for enlarging the Re

A over SOS-IMM. In this example, the initial configuration

27

2.0

1.0

π{2

0.0

3.0

0.0

-4.0

0.0

-1.5

-3.0

0.0

-2.0

-4.0

9θ1

9θ2

τ2

τ3

θ2

θ1

9θ3

(rad)

θ3

(rad/s)

(Nm)

(Nm)

dilated time scale

normal time scale

0.0

0.004

time (s)

1.5

3.0

Figure 2.9: Effective enlargement of Re
A using SOS-CHART-IMM: Plot of joint angles, their ve-
locities, and control inputs with time for stabilization of the Tiptoebot from the initial configuration
given by (2.3.1).

of the Tiptoebot is chosen to be

x “„q˚ 9q˚T

“„´0.2 0.3 0.3 0.0 0.0 0.0T

(2.5.2)

This initial configuration, shown in Fig.2.10 (a), lies outside Ωchart, and hence outside Ωsos. The
Apq˚q of Ωsos is also a null set.
Therefore, the IMM cannot be used to enlarge Ωsos. However, it can be seen from Fig.2.10 (b)

rRe
Apq˚q of Ωsos is found to be a null set and consequently pRe
Apq˚q of Ωchart is not a null set and the IM intersects it. The pRe
that rRe

A of Ωchart is a convex
region on IM and is shown in Fig.2.10 (c). To take the configuration of the Tiptoebot from outside
Re
A “ Ωchart to inside Ωchart, we choose the following velocity configuration that lies in the interior

28

A of Ωchart

rRe

9q3

chart

A of Ω

pRe

p0, 0, 0q

p2.0, ´3.3, ´0.3q

9q1

IM

IM

(c)

9q2

(a)

(b)

A
of Ωsos is a null set, (c) intersection of the slice of Ωchart with the impulse manifold IM, shown in
the original coordinates.

Figure 2.10: (a) Initial configuration of the Tiptoebot, given by (2.5.2), lies outside Ωchart, (b) rRe
of Ωchart (shown in transformed coorrdinates for better visualization) is a convex hull whereas rRe
of pRe

A of Ωchart:

„ 9q1 9qT

2 “„2.0 ´3.3 ´0.3

(2.5.3)

A

The desired effect of the impulsive inputs is shown by the dotted line joining the initial and desired

velocity configurations on the IM - see Fig.2.10 (c). To avoid repetition, the simulation results of

high-gain feedback followed by stabilization are not presented.

2.6 Experimental Verification of Enlargement of Re

A using SOS, CHART

and IMM

2.6.1 Hardware Description

Experiments were done with a two-link Pendubot system. A schematic of the pendubot is shown

in Fig.2.11. The equations of motion of the pendubot are in the form of (2.2.2) where n “ 2, m “ 1
and

u “ τ

(2.6.1)

q “»—–

q1

q2

fiffifl “»—–

fiffifl ,

θ2 ´ π
θ1

29

The equations of motion in terms of the variables q1 and q2 can be found in [23]; the lumped

parameters of the pendubot that appear therein are given in Table 2.2.

The torque τ in the shoulder joint (joint with angle q2 “ θ1) of the pendubot is applied by
a 90-volt direct-drive permanent magnet brushed DC motor3 The motor is driven by a power

amplifier4 with a maximum current limit of 25 A, operating in current mode. The motor torque

constant is 0.4 Nm/A and the amplifier gain is 11.5 A/volt. The torque constant and the current

limit together imply that the maximum torque that can be applied by the motor is 10 Nm. The

pendubot was interfaced with a dSpace DS1104 board and the controller was implemented in the

Matlab/Simulink environment for real-time control.

τ

d1

ℓ1

Y

g

θ1

c

.

m

.

m1, I1

X

m2, I2

c.m.

θ2

d2

ℓ2

Figure 2.11: An arbitrary configuration of the two-link pendubot.

Table 2.2: Lumped Parameters for Pendubot in SI units

β1

β2

β3

β4

β5

0.0147

0.0051

0.0046

0.1001

0.0302

3The motor manufacturer is Minnesota Electric Technology.
4The amplifier is a product of Advanced Motion Control.

30

2.6.2 Design of Experiment

The equilibrium configuration chosen for stabilization was x “ pqT , 9qTqT “ 0 ñ pθ1,θ2, 9θ1, 9θ2q “
p0,π, 0, 0q. The SOS method was used to design a stabilizing second-order polynomial controller
uspxq; the controller is given below

2

τ “´ 14.926θ1 ´ 1.951θ1
´ 0.105θ1

` 0.003 9θ1
´ 0.766pθ2 ´ πq 9θ2 ´ 0.018 9θ1

2 ´ 17.728pθ2 ´ πq´ 2.514pθ2 ´ πq2 ´ 4.469θ1pθ2 ´ πq´ 1.224 9θ1
9θ1 ´ 0.121pθ2 ´ πq 9θ1 ´ 2.813 9θ2 ´ 0.055 9θ2

´ 0.671θ1

9θ2

2

9θ2

(2.6.2)

The corresponding estimated region of attraction, Ωsos in (2.2.10), was obtained as Vpxq ď ρ “
0.083 where

Vpq, 9qq “ 1.24 q2

2 ` 3.01 q1q2 ` 0.15 q2 9q2 ` 0.43 q2 9q1 ` 2.03 q1

2 ` 0.18 q1 9q2 ` 0.54 q1 9q1

` 0.005 9q2

2 ` 0.026 9q1 9q2 ` 0.040 9q1

2

(2.6.3)

Subsequently, the enlarged estimated region of attraction was obtained as Ωchart “ ConvpS8q using
the following parameter values:

N “ 729000, ε “ 0.01, Δt “ 0.0035,

k1 “ 10,

k2 “ 18

To demonstrate the enlargement of the Re

A “ Ωchart using the IMM, we consider the following

initial configuration that lies outside Ωchart:

x “„q˚T 9q˚TT

“„0.00 ´0.72 0.25 0.00T
For the initial configuration of the pendubot given by (2.6.4), the slice rRe
Apq˚
1 “ 0.00, q˚
of Ωchart is shown in Fig.2.12 together with the slice of Ωsos. It can be seen that rRe

2 “ ´0.72q
A of Ωchart is
significantly larger than that of Ωsos. The impulse manifold IMp0.25, 0.00q, which is line, is given
by the relation

(2.6.4)

9q1 `p1´ β3{β2q 9q2 ´ 0.25 “ 0

(2.6.5)

and is found to intersect rRe
segment of the impulse manifold (shown in Fig.2.12 using a thicker line) whereas pRe

A of Ωchart but not that of Ωsos. Therefore, pRe

A of Ωchart is a line
A of Ωsos is a

31

8

4

0

)
s
/
d
a
r
(

2
9θ
“
1
9q

A of Ωchart

rRe

IM

-4

-20

0.3

0.0

p2.00, 0.05q

p0.00, 0.25q

0.0

1.0

2.0

A of Ωsos

rRe

20

A ofΩ

pRe

chart

0
9q2 “ 9θ1 (rad/s)

Figure 2.12: Intersection of rRe
lie on pRe

A of Ωchart is the dark
line segment. The impulse manifold IM was determined by the initial velocity configuration
p 9q2, 9q1q “ p0.00, 0.25q. The desired velocity configuration p 9q2, 9q1q “ p2.00, 0.05q was chosen to
A of Ωchart; due to
this scaling the initial and desired velocity configurations appear to be close.

A of Ωchart with IM; the intersection set pRe
A of Ωchart. A zoomed-out scale has been used to show the entire rRe

null set. To take the configuration of the pendubot from outside Re

A to inside Re

A, we now choose

the following velocity configuration that lies in the interior of pRe
9θ1 “„0.05 2.00

„ 9q1 9q2 “„ 9θ2

A of Ωchart:

(2.6.6)

Our experimental verification involved the following steps:

1. Feedback Linearization: A controller was designed using feedback linearization [73] to
fix the configuration of the first link at pq2, 9q2q “ pθ1, 9θ1q “ p´0.72, 0.00q. With the first
link fixed, the second link behaves as a simple pendulum. To achieve the configuration
pq1, 9q1q “ pθ2 ´ π, 9θ2q “ p0.00, 0.25q, the second link was manually taken to an angle (at
time t0) and released from rest; this angle was determined through trial and error. At time

t´
1 , the configuration in (2.6.4) is achieved - see Fig.2.13 (a), and the feedback linearizing

controller was switched off.

2. High-Gain Control: To implement the IMM, the high-gain controller in (2.2.7) was switched

on at time t´
1 . For the high-gain controller, we used 9qdes
λ1 “ 1. The high-gain controller was switched off at time t`
a small neighborhood of (2.6.6) - see Fig.2.13 (b).

2 “ 9q2 in (2.6.6), µ “ 0.001, and
1 when the joint velocities reached

32

t0

Y

X

Y

X

Y

t´
1

rt´
1 ,t`
1 s

t`
1

t f

(a)

(b)

(c)

Figure 2.13: Steps of experimental verification using the pendubot shown in Fig.2.11: (a) rt0,t´
1 q -
feedback linearizing control, (b) rt´
1 ,t fs - SOS stabilizing control

1 s - high-gain feedback, (c) pt`

1 ,t`

3. Stabilizing Control: The SOS controller in (2.6.2) was invoked at time t`

1 . The pendubot

configuration reaches the desired equilibrium at time t f - see Fig.2.13 (c).

2.6.3 Experimental Results

The experimental results are shown in Fig.2.14. Simulation results are presented in Fig.2.15 to

compare the experimental results with the ideal behavior. The total time scale is divided into

three intervals corresponding to the three steps discussed in the last section. In the experiment,

the purpose of first step was to create the initial configuration of the pendubot in (2.6.4).

In

the simulation, the first step was not necessary as the pendubot was simply assumed to have this

configuration. In the experiment, at time t0, the second link was released from θ2 “ ´2.28 rad
« ´131 deg with the feedback linearizing controller active. At time t´
1 , the feedback linearizing
controller was switched off and the high-gain controller switched on since the configuration of the

pendubot was close to the desired initial configuration in (2.6.4); the pendubot configuration was
pθ1,θ2, 9θ1, 9θ2q “ p´0.72,π, 0.00, 0.42q - see Fig.2.14.

For both the experiment and the simulation, the high gain controller was switched on at time

t´
1 . For better visualization, the effect of high-gain control is illustrated using a dilated time scale.

Due to the initial error in the actual and desired velocities of the actuated joint, the high-gain

controller invokes a large torque at time t´

1 . Through our simulation, this torque was found to be

33

feedback

linearization

high-gain

control

-2.28

1.92

0.42

3.14

0
-2
4

0
-2

3

0

10

0

stabilizing

control

θ2 (rad)

θ1 (rad)

9θ1 (rad/s)

9θ2 (rad/s)

τ (Nm)

t0
0.00

t´
1
0.995

t`
1
1.03

time (s)

t f
4.85

Figure 2.14: Experimental verification of enlargement of Re
A of pendubot: Plot of joint angles,
their velocities, and control input with time. The effect of high-gain feedback is shown using a
dilated time scale; the effect of the stabilizing controller is shown using a normal time scale.

3.14

-0.72

1.92

0.42

0
-2

10

0
-5
2
0

-4
20

10

0

high-gain

control

stabilizing

control

θ2 (rad)

θ1 (rad)

9θ1 (rad/s)

9θ2 (rad/s)

τ (Nm)

t´
1
0.0

t`
1
0.003

time (s)

t f
5.0

Figure 2.15: Simulation results of ideal pendubot behavior for comparison with experimental
results presented in Fig.2.14.

34

of the order of 20 Nm. In the ideal case, the high-gain controller remains active for a very short

interval of time and the control input decays exponentially during this interval [23]. This can be

verified from the simulation where the time interval is equal to 0.003 s. In our experiment, due to

the current limit (25 A) of the power amplifier (see section 2.6.1), the torque applied by the motor

was saturated at 10 Nm. Therefore, the high-gain controller remained active for a longer duration

of time 0.035 s and the control torque decays after remaining saturated for a significant fraction

of that interval. The high-gain torque was switched off at time t`
1 when the joint velocities were
close to the desired values in (2.6.6); these values are p 9θ1, 9θ2q “ p1.92, 0.00q for our experiment and
p 9θ1, 9θ2q “ p1.96, 0.05q for our simulation. The changes in the joint angles for both the simulation
and the experiment were negligible, as expected. Although small, the difference in the joint

velocities between the experiment and the simulation at time t`

1 can be attributed to several factors.

These include: (a) error between the actual and desired configuration at time t´

1 , (b) saturation

of the high-gain torque, and (c) sensitivity of the slope of the impulse manifold (line) IM to the

system parameters - see (2.6.5).

The SOS stabilizing controller was switched on at time t`

1 and the pendubot reached its

equilibrium configuration shortly thereafter. The time required for stabilization and the transient

behavior of the system are similar for both the experiment and the simulation. A video of the

experiment can be found on the weblink:

https://www.egr.msu.edu/~mukherji/RofAEnlargeStableUnstable.mp4

35

IMPULSIVE DYNAMICS AND CONTROL OF THE INERTIA-WHEEL PENDULUM

CHAPTER 3

3.1 Introduction

The problem of designing rest-to-rest maneuvers for the inertia-wheel pendulum using a se-

quence of impulsive inputs is addressed ˙The use of impulsive inputs simplifies the dynamics of the

IWP and permits the design of optimal sequences under some general assumptions1. This chapter

is organized as follows:

the mathematical model of the inertia-wheel pendulum and the formal

problem statement is presented in section 3.2. In section 3.3, rest-to-rest maneuver is designed us-

ing a sequence of two impulsive inputs and a generalization of this result which utilizes a sequence

of N impulsive inputs is presented in section 3.4. In section 3.4, it is also proven that a sequence

with an odd number of inputs is less optimal than the two adjacent sequences with even number

of inputs. Finally simulation results are presented in section 3.5 which provides insights into high

wheel velocity for the inertia-wheel pendulum obtained in prior works.

y

g

c.m.

τ

ℓ2

ℓ1

φ

motor

θ

pendulum

inertia wheel

x

Figure 3.1: The inertia wheel pendulum is shown in an arbitrary configuration.

1The frequently used symbols in this chapter are defined earlier - see key to symbols before

chapter 1.

36

3.2 Mathematical Model and Problem Statement

3.2.1 Equations of Motion

The Inertia-Wheel Pendulum (IWP), shown in Fig.3.1, is an underactuated system comprised of

a simple pendulum and an inertia wheel. The inertia wheel, also referred to as a reaction wheel,

is driven by a motor mounted at the distal end of the pendulum. The kinetic energy T and the

potential energy V of the IWP are given by the relations:

2

,

9θ2 ` m22p 9θ` 9φq
2 ` I1¯{2,

1 ` m2l

2

V “ β sinθ

m22 “ I2{2

T “ m11

m11 “´m1l

2

β “ pm1l1 ` m2l2qg

Using (3.2.1), the Euler-Lagrange equation of motion of the IWP can be written as:

2pm11 ` m22q :θ` 2 m22

:φ` β cosθ “ 0
:θ` 2 m22
:φ “ τ

2 m22

(3.2.1)

(3.2.2a)

(3.2.2b)

Note that the angle φ does not appear in (3.2.2a) or (3.2.2b) - this is because the wheel is symmetric

about its axis of rotation.

3.2.2 Effect of Impulsive Input and Constants of Free Motion

The application of an ideal impulsive torque τi by the motor at time ti results in discontinuous

jumps in the velocities of both the pendulum and the wheel while their angular positions remain
unchanged. By integrating (3.2.2a) with respect to time over the infinitesimal interval rt´
can be shown [23] that these jumps satisfy the relation:

i s, it

i ,t`

p 9θ`
i ´ 9θ´

i q “ ´Cp 9φ`

i ´ 9φ´
i q,

C fi

m22

pm11 ` m22q

(3.2.3)

37

Using (3.2.2b), the jump in the pendulum velocity can be shown to be related to the angular impulse

as follows:

Ii fiż t`

i
t´
i

τi dt “ 2 m22”p 9θ`

i qı
i ´ 9φ´

i ´ 9θ´

i q`p 9φ`

ñ Ii “ ´2 m11p 9θ`

i ´ 9θ´
i q

where (3.2.4b) was obtained from (3.2.4a) using (3.2.3).

(3.2.4a)

(3.2.4b)

If τ “ 0, the dynamics of the IWP described by (3.2.2a) or (3.2.2b) simplifies to the form

2 m11

:θ` β cosθ “ 0
:θ` :φ “ 0

(3.2.5a)

(3.2.5b)

The above equations are the differential forms for conservation of energy of the IWP and conser-

vation of angular momentum of the wheel, namely

m11

9θ2 ` β sinθ “ constant
9θ` 9φ “ constant

(3.2.6a)

(3.2.6b)

It can be seen from (3.2.6a) and (3.2.6b) that the dynamics of the pendulum and wheel are decoupled

and simplified when τ “ 0. For controlling the IWP, this motivates the use of impulsive torques at
discrete instants of time.

3.2.3 Problem Statement: Rest-to-Rest Maneuvers

We consider the problem of designing a set of discrete impulsive inputs over the time interval rt0,t fs
for rest-to-rest maneuvers of the IWP between two configurations where the final configuration has

higher potential energy than the initial configuration, i.e.

9θ0 “ 9φ0 “ 0
9θf “ 9φf “ 0
βpsinθf ´ sinθ0q ą 0

38

(3.2.7a)

(3.2.7b)

(3.2.7c)

and θ lies in the range p´3π{2,π{2s. For a given value of N, the task is to design the vector
I N “ rI1 I2 ¨¨¨ INs where

signpIkq “ ´signpIk´1q,

k “ 2, 3,¨¨¨ , N

(3.2.8)

and t0 ă t1 ă t2 ă ¨¨¨ ă tN ă t f . The rationale for imposing the constraint in (3.2.8) will be provided
later, at the beginning of section 3.4. For the sake of simplicity, it is assumed that the time instants

ti are such that the angular velocity of the pendulum is momentarily zero, i.e.,

9θ´
i

fi 9θpt´

i q “ 0,

i “ 1, 2,¨¨¨ , N

(3.2.9)

Without loss of generality, it is assumed that the first impulsive torque is applied at time t1, where

t´
1 “ t0. Therefore,

9θ´
1 “ 9θ0 “ 0,

9φ´
1 “ 9φ0 “ 0

(3.2.10)

Since the angular positions of the pendulum and the wheel do not change during application of an

impulsive torque, we have

θ1 “ θ0,

φ1 “ φ0

(3.2.11)

3.3 Rest-to-Rest Maneuvers: Case of One and Two Impulsive Inputs

3.3.1 Case of One Impulsive Input (N “ 1)
Using (3.2.3) and (3.2.10), we can write

9θ`
1 “ ´C 9φ`

1

(3.3.1)

The main result can now be stated as follows:

Result 1: (One Impulsive Input) The rest-to-rest maneuver described by (3.2.7) cannot be accom-
plished using I 1 “ rI1s.

39

Discussion: Since τ “ 0 for rt`

1 ,t fs, we get from (3.2.6), and (3.2.11):

2

9
θ`
1

m11

2

` β sinθ0 “ m11
9θ`
1 ` 9φ`

9θf
1 “ 9θf ` 9φf

` β sinθf

(3.3.2a)

(3.3.2b)

Using (3.2.7b), (3.3.2b), and (3.3.1), it can be shown that

9θ`
1 “ 0. From (3.3.2a) we now get

βpsinθf ´ sinθ0q “ 0, which violates (3.2.7c). This establishes Result 1 by contradiction.

A single impulsive torque can always be chosen to impart sufficient angular momentum to

the pendulum such that it reaches its desired configuration with zero angular velocity. However,

this impulsive torque will also cause the wheel to have a nonzero angular velocity in the inertial
reference frame. This implies that 9φ ­“ 0 when 9θ “ 0.

3.3.2 Case of Two Impulsive Inputs (N “ 2)
It is assumed that the first impulse is applied at time t1 and therefore (3.3.1) is still valid. Two

results are presented next. In the first result (Result 2), we relax the assumption in (3.2.9) and design

a more general sequence of impulses that satisfies (3.2.7). We will show that the second result

(Result 3), which is a special case of the first, automatically satisfies the assumption in (3.2.9).

Result 2: (Two Impulsive Inputs) The rest-to-rest maneuver described by (3.2.7) can be accom-
plished using I 2 “ rI1 I2s, where I2 “ ´I1.
Discussion: From (3.2.4a) and (3.2.7a) we have

and since (3.2.6b) holds good for rt`

1 ` 9φ`
1 q

I1 “ 2 m22p 9θ`
1 ,t´
2 s, we can write

I1 “ 2 m22p 9θ`

1 ` 9φ`

1 q “ 2 m22p 9θ´

2 ` 9φ´
2 q

(3.3.3)

From (3.2.7b) we have p 9θf ` 9φfq “ 0 and since (3.2.6b) holds good for rt`
p 9θ`
2 ` 9φ`

2 q “ 0. Using (3.2.4a), we can now write

2 ,t fs, we can write

I2 “ ´2 m22p 9θ´

2 ` 9φ´
2 q

(3.3.4)

40

It is clear from (3.3.3) and (3.3.4) that the conditions in (3.2.7b) require I2 “ ´I1.

For rt`

1 ,t´

2 s and rt`

2 ,t fs, the conservation laws in (3.2.6), together with (3.2.7b) and (3.2.11),

give

2

9
θ`
1

m11

2

9
θ`
2

m11

` β sinθ2

2

9
θ´
` β sinθ0 “ m11
2
1 ` 9φ`
2 “ 9θ`
9θ´
2 ` 9φ´
` β sinθ2 “ β sinθf
9θ`
2 ` 9φ`

2 “ 9θf ` 9φf “ 0

1

(3.3.5a)

(3.3.5b)

(3.3.5c)

(3.3.5d)

For the second impulse, the relationship between the velocity jumps can be obtained from (3.2.3)

as:

2 ´ 9θ´
p 9θ`
2 ´ 9φ´
2 q
2 from (3.3.5d) and 9φ´
2 “ ´ 9θ`
2 “ 9θ`
the right-hand side of (3.3.6) and simplifying using (3.3.1), we get

By substituting the relations 9φ`

2 q “ ´Cp 9φ`

(3.3.6)

1 ` 9φ`

1 ´ 9θ´

2 from (3.3.5b) into

p 9θ`
2 ´ 9θ´

2 q “ ´ 9θ`

1

(3.3.7)

By combining (3.3.5a) and (3.3.5c) and substituting (3.3.7), we get

β`sinθf ` sinθ0 ´ 2 sinθ2˘ “ m11„ 9

θ`

2

“ 2m11

9
θ´
2

2

9
θ´
2

9
θ`
1

´

2

2

`
9
θ`
2

Substituting the expressions for

9
2 and 9
θ´
θ`
2 from (3.3.5a) and (3.3.5c) in the above equation and

simplifying, we get

| 9θ`

1 |“

1

2d β

m11

`sinθf ´ sinθ0˘
asinθf ´ sinθ2

From (3.2.7c) we know that psinθf ´ sinθ0q ą 0 and it can be seen from (3.3.5c) that psinθf ´
sinθ2q ą 0; therefore, the above equation is well-defined. For a given pair tθ0,θfu, (3.3.8) provides
a functional relationship between the initial angular velocity 9θ`
1 (resulting from application of the
first impulse I1) and the configuration θ2 where the second impulse (I2 “ ´I1) is applied. There

(3.3.8)

41

Final Configuration

Intermediate Configuration

y

θf

θ0

x

rt`
2 ,t f s

θ2

rt`
1 ,t´
2 s

Initial Configuration

Figure 3.2: An example showing the initial, intermediate, and final configuration of the IWP for
the case with two impulsive inputs (N “ 2).
are infinite solutions given by the pair t 9θ`
1 can be used to
compute the impulses I1 and I2 (I2 “ ´I1), using (3.3.1), (3.3.3), and (3.3.4). An example
showing the initial, intermediate, and final configuration of the IWP for the case with two impulsive

1 ,θ2u; for each solution, the value of 9θ`

inputs is shown in Fig.3.2. The next result pertains to the particular solution that minimizes the

magnitude of the impulses.

Result 3: (Optimal Input I 2) The minimum magnitude of the impulsive inputs required for the

rest-to-rest maneuver described by (3.2.7) is

| I1|“| I2|“b2m11β`sinθf ´ sinθ0˘

Discussion: Using (3.3.1) and (3.3.3), I1 can be expressed as

I1 “ 2 m22p 9θ`

1 ` 9φ`

1 q “ 2 m22p1´Cq 9θ`

1

(3.3.9)

(3.3.10)

Therefore, the magnitude of I1 can be minimized by minimizing the magnitude of 9θ`

1 . From

42

(3.3.5a) it can be seen that

2

9
θ`
1

m11

2

9
θ´
2

“ m11

` βpsinθ2 ´ sinθ0q

9θ´
2 ­“ 0 but the minimum value of 9θ`

If psinθ2 ´ sinθ0q ď 0,
1 is equal to zero. This implies
psinθf ´ sinθ0q “ 0 from (3.3.8), which contradicts (3.2.7c). Since psinθ2 ´ sinθ0q must be
positive, the minimum magnitude of 9θ`
2 “ 02; this magnitude is

1 can be obtained by choosing 9θ´

equal to

1 |“d β
| 9θ`

m11 psinθ2 ´ sinθ0q

By equating (3.3.8) and (3.3.11), we get

sinθ2 “

1

2psinθ0 ` sinθfq

(3.3.11)

(3.3.12)

where θ2 is the angle at which the second impulse is applied3. Substituting (3.3.12) into (3.3.10)

and (3.3.11), and comparing (3.3.3) and (3.3.4) we get

1 | “d β
| 9θ`

2m11`sinθf ´ sinθ0˘

ñ | I1 | “| I2 |“b2m11β`sinθf ´ sinθ0˘

(3.3.13a)

(3.3.13b)

This establishes Result 3.

From Result 3 it can be seen that the time instant t2 is automatically known when the magnitudes

of the individual impulses are minimized. This is different from Result 2, where the choice of t2 is

not unique and each feasible choice of t2 (alternatively θ2) uniquely determines the magnitudes of

the impulses.

2This choice automatically satisfies the assumption in (3.2.9).
3It is clear that (3.3.12) can have multiple solutions for θ2. The procedure for computing the
correct solution will be discussed in section 3.4.2. Knowing the value of θ2, it will be possible to
determine the time instant t2.

43

3.4 Rest-to-Rest Maneuvers: Generalization to N Inputs

3.4.1 Revisiting the Problem Statement

We start this section with the result that justifies the rationale for imposing the constraint in (3.2.8).

Result 4: (Sum of Two Consecutive Impulses) Consider two impulses Ik and Ik`1 applied at

times tk and tk`1. The net effect of these two impulses in terms of change in the velocities of the
pendulum and wheel over the interval rt´
¯I at time tk
over the interval rt´

k`1s can be achieved by a single impulse

k s where

k ,t`

k ,t`

Discussion: From (3.2.4a) we have

¯I “ Ik ` Ik`1

(3.4.1)

k`1s, (3.2.6b) can be used to rewrite the above equations as follows

Since τ “ 0 for rt`

k ,t´

k ´ 9θ´
k`1 ´ 9θ´

Ik “ 2 m22”p 9θ`
Ik`1 “ 2 m22”p 9θ`

k qı
k q`p 9φ`
k ´ 9φ´
k`1qı
k`1 ´ 9φ´
k`1q`p 9φ`
k qı
k`1 ´ 9φ´
k`1qı
k`1 ´ 9φ´
k qı
k`1 ´ 9φ´
The change in the velocities of the pendulum and wheel over the interval rt´
k q, respectively. To achieve the same change over rt´
k ,t`
and p 9φ`
k s,
k qı
k`1 ´ 9φ´

Ik “ 2 m22”p 9θ´
Ik`1 “ 2 m22”p 9θ`
ñ Ik ` Ik`1 “ 2 m22”p 9θ`

k q`p 9φ´
k`1q`p 9φ`
k q`p 9φ`

k`1 ´ 9θ´
k`1 ´ 9θ´
k`1 ´ 9θ´

k`1 ´ 9θ´

k q`p 9φ`

k`1 ´ 9φ´

¯I “ 2 m22”p 9θ`
ñ ¯I “ Ik ` Ik`1

k`1s are p 9θ`
k ,t`
¯I must satisfy

k`1´ 9θ´
k q

This establishes Result 4.

Result 4 clearly indicates that two consecutive impulses of the same sign can be replaced by a

single impulse of the same sign. This justifies the constraint imposed in our problem statement that

consecutive impulses must have opposite sign.

44

An extension of Result 4 is now considered. For a rest-to-rest maneuver using two impulsive
¯I “ I1 ` I2 “ 0. This is true since (3.2.6b) and (3.2.7b)
f q “ 0, and (3.2.7a) and (3.2.10) implies p 9θ´
k`1q “ p 9θ`
k `
1 q “ 0. This is consistent with Result 2, where it was shown that I2 “ ´I1. A

inputs pN “ 2q, Result 4 implies that
implies p 9θ`
2 ` 9φ`
9φ´
k q “ p 9θ´
generalization of this result in stated next.

k`1 ` 9φ`
1 ` 9φ´

2 q “ p 9θ`

f ` 9φ`

Result 5: (Zero Sum of Impulses) For a rest-to-rest maneuver involving N impulsive inputs, N ě 2,
the following equation must hold.

Nÿi“1

Ii “ 0

(3.4.2)

Discussion: A sequence of N impulses, N ě 2, can be replaced by two impulses by applying Result
4 iteratively. For a rest-to-rest maneuver, the sum of these two impulses is zero - this follows from

our discussion above.

It is clear from the discussion above that both Result 4 and Result 5 are quite general and they

do not require the assumption in (3.2.9) to be satisfied.

With the motivation of investigating the minimum values of the magnitudes of the impulsive

torques, we investigate rest-to-rest maneuver of the IWP with I 3 and I 4. As in the cases with

I 1 and I 2, (3.3.1) holds good.

3.4.2 Rest-to-Rest Maneuvers with Even Number of Impulsive Inputs

Theorem 1. (Optimality of Even Impulse Sequence)

For a rest-to-rest maneuver of the IWP that satisfies (3.2.7) and (3.2.9) and uses 2n impulsive

inputs, n “ 1, 2,¨¨¨ , I 2n

8 is minimized by the following choice of inputs:

| Ii|“

@ i “ 1, 2,¨¨¨ , 2n

(3.4.3)

The angles where the impulsive inputs are applied satisfy the following relation

1

?nb2m11β`sinθf ´ sinθ0˘,
2n  sinθf
 sinθ0 `„ i´ 1

sinθi “„ 2n´ i` 1

2n

45

i “ 1, 2,¨¨¨ , 2n

(3.4.4)

Proof: We use induction to first prove (3.4.3). Assuming that (3.4.3) is satisfied for 2m impulsive

inputs, i.e., n “ m. We express the magnitudes of the impulses for the case with p2m` 2q inputs
using the relation

| Ii|“

ki?mb2m11β`sinθf ´ sinθ0˘,

i “ 1, 2,¨¨¨ ,p2m` 2q

(3.4.5)

where ki, i “ 1, 2,¨¨¨ ,p2m` 2q, are arbitrary positive numbers. Using (3.2.4b) and (3.2.9) we can
show

ki?md β

| 9θ`

i

| “

2m11`sinθf ´ sinθ0˘,

i “ 1, 2,¨¨¨ ,p2m` 2q

(3.4.6)

Using (3.2.9), the conservation law in (3.2.6a) for the time intervals rt`
¨¨¨ ,rt`

2m`2,t fs can be written as:

1 ,t´

2 s, rt`

2 ,t´

3 s, ¨¨¨rt`

j ,t´

j`1s,

2

2

9
θ`
1

9
θ`
2

m11

m11

m11

2

9
θ`
j

` β sinθ0 “ β sinθ2
` β sinθ2 “ β sinθ3

...

` β sinθj “ β sinθj`1

...

m11

9θ`2
2m`2 ` β sinθ2m`2 “ β sinθf

(3.4.7)

Addition of equations in (3.4.7) and substitution of | 9θ`
gives the following

i

| from (3.4.6) in the resulting equation

i “ rk1 k2 ¨¨¨ k2m`2s2
Using (3.4.8) and the property of norms it can be shown that

k

2

2 “ 2m

2m`2ÿi“1

(3.4.8)

(3.4.9)

?2m` 2 rk1 k2 ¨¨¨ k2m`2s8 ě rk1 k2 ¨¨¨ k2m`2s2
ñ rk1 k2 ¨¨¨ k2m`2s8 ěc m

m` 1

It can be shown that ki “am{pm` 1q, i “ 1, 2,¨¨¨ ,p2m` 2q, satisfy (3.4.8) and minimize

rk1 k2 ¨¨¨ k2m`2s8. Substitution of these values of ki in (3.4.5) shows that (3.4.3) is satisfied for
p2m` 2q impulsive inputs, i.e. n “ pm` 1q. It has been shown earlier in (3.3.9) that (3.4.3) is

46

| 9θ`

i

| “

1

?nd β

2m11`sinθf ´ sinθ0˘

(3.4.10)

satisfied for N “ 2 (n “ 1). By induction we can now claim that (3.4.3) will be satisfied for any
even number of impulsive inputs.

The values of | 9θ`

i

|, i “ 1, 2,¨¨¨ , 2n, can be obtained from (3.4.3), (3.2.4b), and (3.2.9), namely

By substituting (3.4.10) in the energy conservation laws (similar to (3.4.7)) for 2n impulses and

solving sequentially as shown below, we get the relations between the angles where the impulsive

inputs should be applied

S1 “ S0
S2 “

1

S3 “
...

Si “

2n

2n pS f ´ S0q` S0 “„ 2n´ 2` 1
2n pS f ´ S0q` S2 “„ 2n´ 3` 1
2n pS f ´ S0q` Si´1 “„ 2n´ i` 1

 S0 `„ 2´ 1
2n  S f
 S0 `„ 3´ 1
2n  S f
 S0 `„ i´ 1
2n  S f

2n

2n

1

1

(3.4.11)

where Sp “ sinθp, p “ 1, 2,¨¨¨ f . It can be seen that these angles satisfy 3.4.4.

It follows from (3.4.3) in Theorem 1 and (3.2.8) that all the 2n impulses have the same magnitude

and consecutive impulses have opposite signs. Since the pendulum and wheel are both at rest at

the initial time, it follows that each pair of consecutive impulses (starting with I1 and I2) result

in a rest-to-rest maneuver.

It follows from (3.2.4b), (3.2.8), and (3.2.9) that the velocity of the pendulum immediately

k`1q “ ´signp 9θ`

after application of an impulsive input will have opposite sign for two consecutive impulses, i.e.,
signp 9θ`
constructed to determine the unique value of θi, i “ 1, 2,¨¨¨ , 2n, from (3.4.4).
If θf belongs to quadrant I (includes θ “ π{2) or IV, then

k q, k “ 1, 2,¨¨¨ ,p2n´ 1q. Using this fact, the following algorithm can be

For m “ 1, 2,¨¨¨ , n, compute sinθ2m using (3.4.4)

47

If sinθ2m ą 0, θ2m belongs to quadrant II
Elseif sinθm ă 0, θ2m belongs to quadrant III
Else θ2m “ ´π.
For m “ 0, 1,¨¨¨ ,pn´ 1q, compute sinθ2m`1 using (3.4.4)
If sinθ2m`1 ą 0, θ2m`1 belongs to quadrant I
If sinθ2m`1 ă 0, θ2m`1 belongs to quadrant IV
Else θ2m`1 “ 0.

Else θf belongs to quadrant II or III, then

For m “ 1, 2,¨¨¨ , n, compute sinθ2m using (3.4.4)
If sinθ2m ą 0, θ2m belongs to quadrant I
Elseif sinθm ă 0, θ2m belongs to quadrant IV
Else θm “ 0.
For m “ 0, 1,¨¨¨ ,pn´ 1q, compute sinθ2m`1 using (3.4.4)
If sinθ2m`1 ą 0, θ2m`1 belongs to quadrant II
If sinθ2m`1 ă 0, θ2m`1 belongs to quadrant III
Else θ2m`1 “ ´π.

Endif

3.4.3 Rest-to-Rest Maneuvers with Odd Number of Inputs

We generalize the result presented in Remark 4.

Theorem 2. (Lack of Optimality of Odd Impulse Sequence)

It is not possible to design an odd impulse sequence for which the magnitudes of all the impulsive

inputs are less than the optimal magnitude for the preceding and succeeding even impulse sequence.

In other words, the following inequality holds for n “ 1, 2,¨¨¨ .

I 2n`1

8 ą min I 2n

8 ą min I 2n`2

8

(3.4.12)

48

Proof:

It is clear from (3.4.3) that min I 2n

8 ą min I 2n`2

8. We proceed to prove the left

inequality by contradiction and therefore assume

| Ii| “ ki min I 2n
ki?nd β

| “

ñ | 9θ`

i

8 |

ki P p0, 1s

2m11`sinθf ´ sinθ0˘

From (3.2.8), (3.4.2), and (3.4.13) we can show

k1 ` k3 `¨¨¨ k2n`1 “ k2 ` k4 `¨¨¨ k2n

Substituting the expression for 9θ`
i

for the p2n` 1q impulses, we get

from (3.4.13) in the energy conservation laws (similar to (3.4.7))

k

2

1 ` k

2

2 `¨¨¨ k

2

2n`1 “ 2n

Since ki P p0, 1s, the following inequality holds true

k

2

1 ` k

2

2 `¨¨¨ k

2

2n`1 ď k1 ` k2 ` k3¨¨¨ k2n`1

Using (3.4.14), and (3.4.15) we get

(3.4.13)

(3.4.14)

(3.4.15)

(3.4.16)

(3.4.17)

ě n

k2 ` k4 `¨¨¨ k2n

l               jh               n

n terms

The left hand side of the inequality above is equal its right hand side @ kq “ 1, q “ 2, 4,¨¨¨ 2n.
However, this choice of kq violates (3.4.14) and (3.4.15) to be satisfied simultaneously. Therefore,

(3.4.17) can be modified as

ą n

(3.4.18)

k2 ` k4 `¨¨¨ k2n

l               jh               n

n terms

There exist no choice of kq such that the above inequality holds true @kq P p0, 1s and thus, completes
the proof by contradiction.

49

3.5 The Swing-Up Problem

3.5.1 Optimal Swing-Up Using Even Impulse Sequences

The swing-up problem is a rest-to-rest maneuver where the final pendulum angle is θf “ π{2.
We consider special case where the initial angle of the pendulum is in the vertically downward

configuration, i.e., θ0 “ ´π{2. For swing-up using an even number of impulsive inputs, the optimal
solution can be obtained from (3.4.3). For N “ 2n, n “ 1, 2,¨¨¨ , the optimal solution is given by a
sequence of equal and opposite impulses of the following magnitude:

(3.5.1)
Since the magnitude of the impulses is inversely proportional to ?n, it is clear that if the number of

@ i “ 1, 2,¨¨¨ , 2n,

| Ii|“

?n

I,

I fi 2am11β

1

impulsive inputs are increased by an even number, the magnitude of each impulse in the sequence

is reduced. This information will be useful for designing an impulse sequence that takes into

consideration actuator saturation.

3.5.2

Implementation Using High-Gain Feedback

Ideal impulsive inputs are Dirac-delta functions and cannot be generated by actuators.

In real

physical systems, continuous-time implementation of impulsive inputs has be achieved using high-

gain feedback in both theory and experiments [19,23,25]. The high-gain feedback [93] for the IWP

can be obtained as

τhg “ rKT M´1

Ks´1rKT M´1
pm11 ` m22q m22
m22

m22

M “ 2»—–

i ´ 9φ´
i qs
β cosθ

0

H `

1

εp 9φ`

fiffifl , H “»—–

fiffifl

(3.5.2)

where the matrices M and H above were reconstructed from (3.2.2), K fi„0 1T

, and ε ą 0 is
a small number. Implementation of impulsive inputs using high-gain feedback also enables us to

compare our results with those published in the literature. A discussion of select results in the

literature is presented next.

50

3.5.3 Discussion of Results in the Literature

3.5.3.1 Globally Stabilizing Controller

We implemented the globally stabilizing controller in [1] using their kinematic and dynamic

parameter values, which are given below

m11 “ 4.83ˆ 10´3, m22 “ 32ˆ 10´6, β “ 37.9ˆ 10´2

(3.5.3)

Due to brevity of space, we omit complete expression of the controller used in [1]; rather, we recall

the controller’s structure which is in the form

c0 c1 µ0 ` c1 µ1 ` c2 µ2 ` c3 µ3

where µ0, µ1, µ2 and µ3 are nonlinear functions of state variables. The controller parameters were

chosen as: c0 “ ´π{10, c1 “ 13, c2 “ 16 and c3 “ 8.0; the results are shown in Fig.3.3. The
plots are slightly different from those presented in [1] but the overall trends are similar. It is clear

from the plot of θ that the control input drives the pendulum directly towards the desired value of

θf “ π{2. There is a small overshoot beyond π{2 but the wheel velocity is extremely high, of the

2.0

0.0

-2.0

0.0

0

-8000

0.0

π{2

θ (rad)

10

0

-5
0.0

10.0

0.2

0.0

9θ (rad/s)

10.0

9φ (rad/s)

τ (Nm)

time (s)

-0.4

0.0

10.0

time (s)

10.0

Figure 3.3: Simulation results for the globally stabilizing controller [1] with controller parameter
values: c0 “ ´π{10, c1 “ 13, c2 “ 16 and c3 “ 8.0.

51

2.0

8

-2.0

0.0

0

θ (rad)

9φ (rad/s)

0
0.0

10.0

0.2

0.0

-5000

0.0

time (s)

-0.4

10.0

0.0

9θ (rad/s)

1.0

10.0

τ (Nm)

time (s)

10.0

Figure 3.4: Simulation results for the globally stabilizing controller [1] with controller parameter
values: c0 “ ´π{10, c1 “ 13, c2 “ 16 and c3 “ 4.5.
order of 8000 rad/s. A change in the controller parameter c3 from 8.0 to 4.5 increases the overshoot

slightly but reduces the maximum wheel velocity by almost 50% - see Figs.3.3 and 3.4.

The above observation can be explained by the analysis presented in section 3.3 even though the

nature of the inputs are completely different (continuous inputs in [1] vs a pair of impulsive inputs,

N “ 2). By changing the domain of θ from p´3π{2,π{2s to p´π{2, 3π{2s4 and using (3.3.1) and
(3.3.8), we get for θf “ π{2 and θ0 “ ´π{2:

1

Cd

| 9φ`

1 |“

β

m11p1´ sinθ2q

(3.5.4)

It is clear from (3.5.4) that the wheel velocity immediately after application of the first impulse

depends only on the angle where the second impulse is applied, namely θ2, and tends to infinity
when θ2 “ π{2`, i.e., when the overshoot approaches zero. While it is clear from (3.5.4) that
θ2 “ π{2` is not a good choice for application of the second impulse, the value of θ2 that minimizes
the magnitude of the wheel velocity | 9φ`
1 | can be obtained using the energy conservation law in
4This change in the domain is necessary to ensure that the trajectory of θ is similar to that in [1]

but it does not change the analysis whatsoever.

52

(3.2.6a). For the IWP to cross the upright configuration, the following inequality must be satisfied:

2

9
θ`
1

m11

` β sinθ0 ą β ñ | 9φ`

1 |ą

1

Ca2β{m11

(3.5.5)

where θ0 “ ´π{2 and (3.3.1) were used. Comparing (3.5.4) and (3.5.5), we can show that
θ2 “ p5π{6q´ minimizes | 9φ`
1 |.

A simulation was performed using the high-gain feedback law in (3.5.2) with ε “ 0.01, the
parameter values in (3.5.3), and θ2 « 3π{4 (slightly less that 5π{6); the results are shown in
Fig.3.5. After the IWP reached a neighborhood of θf “ π{2, a linear controller was invoked for
stabilization. The linear controller was designed to place the poles of the closed loop system at

´4˘ 2i and ´8. The simulation results indicate that swing-up is achieved in less than 1.0 s, which
is much faster than that achieved in [1]. The maximum velocity of the wheel is still quite high

(3000 rad/s) but it is significantly lower than that in Fig.3.3. The torque required is quite high (« 13
Nm) but this can be reduced significantly by simply changing the domain of θ2, as we will show in

the next simulation.

The simulation results presented in Fig.3.5 were obtained by assuming θ P p´π{2, 3π{2s; this
was motivated by the need to generate trajectories of the IWP similar to those generated in [1], for

2.0

0.0

-2.0

0.0

0

θ (rad)

0.5

1.0

-3000

0.0

0.5

9φ (rad/s)

1.0

time (s)

20

0

-10

1.5

0.0

10

0

-10

9θ (rad/s)

0.5

1.0

1.5

τ (Nm)

1.5

0.0

0.5

1.0
time (s)

1.5

Figure 3.5: High-gain feedback implementation of two impulsive inputs (N “ 2) for swing-up of the
IWP. For the purpose of comparison with the results in Figs.3.3 and 3.4, the controller is designed
to keep θ in the domain p´π{2, 3π{2s.

53

comparison. If we switch the domain of θ back to p´3π{2,π{2s, as defined in section 3.2.3, the
maximum wheel speed and the magnitude of the maximum torque can both be reduced from their

values in Fig.3.5. Simulation results of high-gain feedback implementation of the optimal impulse

sequence based on two inputs, described by (3.4.3) and (3.4.4) with n “ 1, is shown in Fig.3.6.
The high-gain controller was implemented using ε “ 0.02 and stabilization of the equilibrium was
achieved by the same linear controller that was used in the last simulation. It can be seen from

Fig.3.6 that the second impulse is applied when θ2 « ´π rad. Similar to the results in Fig.3.5,
swing-up is achieved in less than 1.0 s, but the maximum wheel velocity is now reduced from 3000

rad/s to 2000 rad/s and magnitude of the maximum torque is reduced from « 13 Nm to « 3 Nm.
The maximum torque of « 3 Nm in Fig.3.6, although larger than those reported in the literature,
is not a significant concern because it is applied for a very short duration of time. Motors can

apply substantially larger torques5 than their maximum continuous torque over short time intervals.

The maximum torque of « 3 Nm also corresponds to the continuous-time implementation of the
optimal I 2. The magnitude of this torque, as well as the maximum velocity of the wheel, can

π{2

0.0

θ (rad)

20

0

9θ (rad/s)

-3.0

0.0

2000

0

0.0

1.0

-15

2.0

0.0

1.0

2.0

9φ (rad/s)

3

0

-3

τ (Nm)

1.0
time (s)

2.0

0.0

1.0
time (s)

2.0

Figure 3.6: High-gain feedback implementation of the sequence of two optimal impulsive inputs
(N “ 2) for swing-up of the IWP. The controller is designed to keep θ in the domain p´3π{2,π{2s.
5This is referred to a peak torque [94]; for different motors, the peak torque can be twice to ten

times larger than the maximum continuous torque.

54

be easily reduced if we consider continuous-time implementation of I 2n, n “ 2, 3,¨¨¨ . This is
discussed in the next section.

Our results are not compared with those obtained using the IDA-PBC method, [15, 43], for

example, due to space constraints. Similar to the globally stabilizing controller [1], the IDA-PBC

method also takes the pendulum directly to the desired upright configuration and results in a large

continuous torque [15] or large wheel velocity [43] during swing-up.

3.5.3.2 Energy Based Controller

When the number of impulses are increased from N “ 2 to N “ 8, for example, the magnitude
of the impulsive torques are reduced by a factor of ?4 “ 2; consequently, the magnitude of the

maximum high-gain torque and the wheel velocity are reduced proportionately - see Fig.3.7. The

trajectories of the state variables in Fig.3.7 resemble those of the energy based controllers [2, 68]

during swing-up phase of the IWP; the PFBLC + AL energy-based controller presented in [2] is

simulated here to show the similarities in the trajectories - see Fig.3.8. It can be seen from Figs.3.7

and 3.8, that, unlike the globally stabilizing controller [1] (see Fig.3.3) where the pendulum is

aggressively driven towards its desired configuration, both controllers (presented here and in [2])

0.0

π{2

θ (rad)

-4.0

0.0

950

0
0.0

20

0

-20

4.5

0.0

1.5

0

-1.5

9θ (rad/s)

4.5

τ (Nm)

9φ (rad/s)

time (s)

4.5

0.0

time (s)

4.5

Figure 3.7: High-gain feedback implementation of the sequence of eight optimal impulsive inputs
(N “ 8) for swing-up of the IWP. The high-gain controller was implemented with ε “ 0.02.

55

0.0

π{2

θ (rad)

-4.0

0.0

600

0

-700

0.0

20

0

-20

5.0

0.0

0.15

0

-0.15

5.0

0.0

9φ (rad/s)

time (s)

9θ (rad/s)

5.0

τ (Nm)

time (s)

5.0

Figure 3.8: Simulation using the PFBLC + AL controller in [2]; the controller parameters were
chosen as ke “ 3.1ˆ 107 and kv “ 0.1. This choice of parameters ensured that the time required for
swing-up and the magnitude of the maximum control torque in simulations matched those of the
experiments. The initial configuration was chosen to be slightly different from θ0 “ ´π{2 since
the controller is unable to swing-up from this configuration.

gradually add energy to the pendulum over several cycles of oscillation.

Table 3.1: Swing-up time and maximum magnitude of high-gain torque required for different values
of N

N

ts (s)

max

t0ďtďt f |τhg| (Nm)

2

4

6

8

10

12

0.6

1.36

1.96

2.90

3.67

4.53

3.0

2.12

1.73

1.50

1.35

1.20

A comparison of Figs.3.7 and 3.8 indicates that the magnitude of the maximum torque required

by our method is larger than that required by the approach proposed in [2]. However, since the

torques are applied intermittently over very short intervals of time, feasibility of our approach is

determined by the peak torque rating of the actuator as opposed to the maximum continuous torque

rating, which is always lower [94]. A salient feature of our approach is that, given any actuator,

an optimal impulse sequence can be designed such that the peak torque rating of the motor is not

exceeded. A higher value of N reduces the peak torque requirement of the motor but increases the

56

time time required for swing-up. The swing-up time and the magnitude of the maximum torque for

several different values of N are presented in Table 6.1.

57

IMPULSIVE CONTROL OF A DEVIL-STICK: PLANAR SYMMETRIC JUGGLING

CHAPTER 4

4.1 Introduction

The problem of juggling a devil-stick is investigated. Assuming that the stick remains confined

to the vertical plane, the task is to juggle between two symmetric configurations. Impulsive forces

are applied to the stick intermittently and the impulse of the force and its point of application are

modeled as inputs to the system. The juggling problem is formally described in section 4.2. The

dynamics of the devil-stick is presented in section 4.3; it is comprised of impulsive dynamics due to

the control inputs and continuous dynamics due to torque-free motion under gravity. A coordinate

transformation is used to simplify the control problem and the dynamics is described by a nonlinear

discrete-time system. The control design is provided in section 4.4. By choosing one of the control

inputs to be dead-beat, the nonlinear system is simplified to a linear discrete-time system. For

stable juggling, the linear system is controlled using linear quadratic regulator (LQR) and model

predictive control (MPC) techniques. Simulation results are presented in section 4.5.

y

g

m, J

G ” phx, hyq

ℓ

θ

x

Figure 4.1: A three degree-of-freedom of a devil-stick.

58

4.2 Problem Description

Consider the three degree-of-freedom devil-stick shown in Fig. 4.1, which can move freely

in the xy vertical plane. The stick has length ℓ, mass m, and mass moment of inertia J about its

center-of-mass G. The configuration of the stick is described by the three generalized coordinates:

pθ, hx, hyq, where θ is the orientation of the stick with respect to the positive x axis, measured
counter-clockwise, and phx, hyq are the Cartesian coordinates of G. The objective is to juggle
the stick between two configurations that are symmetric with respect to the vertical axis. The
coordinates of the stick in these two configurations are pθ˚, h˚
yq, where
θ˚ P p0,π{2q - see Fig. 4.2. It is assumed that juggling is achieved by applying impulsive forces
perpendicular to the stick; they are applied only when the orientation of the stick is θ “ θ˚ or
θ “ π´ θ˚. Therefore, the time of application of the impulsive force is not a part of the control
design. The control inputs are the pair pI, rq, where I, I ě 0, is the impulse of the impulsive force
and r is the distance of the point of application of the force from G. The value of r is considered to

yq and pπ´ θ˚,´h˚

x, h˚

x, h˚

be positive if the angular impulse of the impulsive force about G is in the positive z direction when
θ “ θ˚, and is in the negative z direction when θ “ π´θ˚. The control inputs that juggle the stick
between the symmetric configurations are denoted by the pair pI˚, r˚q.

y

g

I “ I˚

r

“

r˚

p´h˚

x , h˚
yq

π´ θ˚

˚

r “ r
ph˚

x , h˚
yq

θ˚

I “ I˚

x

Figure 4.2: Symmetric configurations of the devil-stick in Fig. 4.1.

59

4.3 Dynamics of the Devil-Stick

4.3.1

Impulsive Dynamics

The dynamics of the three-DOF devil-stick is described by the six-dimensional state vector X,

where

X “„θ ω hx vx hy vyT

, ω fi 9θ, vx fi 9hx, vy fi 9hy

Let tk, k “ 1, 2, 3,¨¨¨ , denote the instants of time when the impulsive inputs are applied. Further-
more, without loss of generality, let k “ p2n´ 1q, n “ 1, 2,¨¨¨ denote the instants of time when the
impulsive inputs are applied at θ “ θ˚, and k “ 2n, n “ 1, 2,¨¨¨ denote the instants of time when the
impulsive inputs are applied at θ “ π´ θ˚. If t´
k denote the instants of time immediately
before and after application of the impulsive inputs, the linear and angular impulse-momentum

k and t`

relationships can be used to describe the impulsive dynamics1 as follows, for k “ 1, 3, 5,¨¨¨

Xpt`

k q “ Xpt´

k q`

and for k “ 2, 4, 6,¨¨¨

Xpt`

k q “ Xpt´

k q`

(4.3.1)

(4.3.2)

0

pIk rk{Jq

0

´pIk{mq sinθ˚

0

pIk{mq cosθ˚

0

´pIk rk{Jq

0

pIk{mq sinθ˚

0

pIk{mq cosθ˚

fiffiffiffiffiffiffiffiffiffiffiffiffiffiffifl
fiffiffiffiffiffiffiffiffiffiffiffiffiffiffifl

»——————————————–
»——————————————–

60

1Impulsive inputs cause discontinuous jumps in the velocity coordinates but no change in the
position coordinates. The dynamics of underactuated systems subjected to impulsive inputs is
discussed in [6, 25, 74, 95].

where pIk, rkq denote the control inputs at time tk. Between two consecutive impulsive inputs, the
devil-stick undergoes torque-free motion under gravity; this is discussed next.

4.3.2 Continuous-time Dynamics

Over the interval t P rt`
k`1s, the devil-stick will be in flight; its center-of-mass G will undergo
projectile motion and its angular momentum will remain conserved. This dynamics is described

k ,t´

by the differential equation:

where the initial condition Xpt`
k is odd or even.

9X “„ω 0 vx 0 vy ´gT

(4.3.3)

k q can be obtained from ( 4.3.1) or ( 4.3.2), depending on whether

4.3.3 Poincaré Sections and Half-Return Maps

For the hybrid system, described by impulsive dynamics of section 4.3.1 and continuous dynamics

of section 4.3.2, we define two Poincaré sections2, 3 [96] Sr and Sl as follows:

Sr : tX P R
Sl : tX P R

6 | θ “ θ˚u
6 | θ “ π´ θ˚u

(4.3.4)

These Poincaré sections are chosen since the impulsive inputs are applied only when θ is equal to
θ˚ or pπ´ θ˚q. Any point on Sr and Sl can be described by the vector Y , Y Ă X, where

Y “„ω hx vx hy vyT

(4.3.5)

(4.3.6)

The map Pr : Sr Ñ Sl can be determined from ( 4.3.1) and ( 4.3.3) as follows:

Ypt´

k`1q “ AYpt´

k q` Br

2Poincaré sections have been previously used for design of gaits for bipedal robots [76, 80].
3It is assumed that the initial conditions of the devil-stick are such that its trajectory intersects

one of the two Poincaré sections before the first impulsive control input is applied.

61

A fi

0 0

0

, Br fi

»——————————–

1 0

0

0

0 1 δk 0

1

0

0

0 0

0 0

1 δk

0

1

0

0

0

fiffiffiffiffiffiffiffiffiffiffifl

»——————————–

pIk rk{Jq

´pIk{mq sinθ˚δk
´pIk{mq sinθ˚

pIk{mq cosθ˚δk´p1{2qgδ2

k

pIk{mq cosθ˚´gδk

fiffiffiffiffiffiffiffiffiffiffifl

where δk

fi pt´

k`1 ´ t´

k q and k “ p2n´ 1q, n “ 1, 2,¨¨¨ . Similarly, the map Pl : Sl Ñ Sr can be

determined from ( 4.3.2) and ( 4.3.3) as follows

Ypt´

k`1q “ AYpt´

k q` Bl

(4.3.7)

fi

Bl

»——————————–

´pIk rk{Jq
pIk{mq sinθ˚δk
pIk{mq sinθ˚

pIk{mq cosθ˚δk´p1{2qgδ2

k

pIk{mq cosθ˚´gδk

fiffiffiffiffiffiffiffiffiffiffifl

where k “ 2n, n “ 1, 2,¨¨¨ . Both Pr and Pl in ( 4.3.6) and ( 4.3.7), respectively, can be viewed as
half-return maps4 since the composition of these maps are the return maps Pr ˝ Pl : Sl Ñ Sl and
Pl ˝ Pr : Sr Ñ Sr. In the next section we introduce a coordinate transformation to show that the map
Pl, in the transformed coordinates, is identical to Pr. This simplifies the analysis of the problem.

4.3.4 Coordinate Transformation

Consider Fig. 4.3, where z “ 0 denotes the xy plane in which the devil-stick is juggled. Typically,
the juggler will stand at a point on the positive z axis, denoted by P in Fig. 4.3 (a), and face the
z “ 0 plane. The juggler will apply a control action with the right hand when θ “ θ˚, and with the
left hand when θ “ π´ θ˚, i.e., the juggler is ambidextrous. Instead of alternating between the
right and left hands, the juggler can choose to apply all control actions using the right hand only.

4Half-return maps have been used to analyze the behavior of dynamical systems such as the van

der Pol oscillator [97, 98].

62

(a)

y

g

p´h˚

x,h˚
yq

π´θ˚

ph˚

x,h˚
yq

θ˚

Q

x

y

z

ph˚

x,h˚
yq

θ˚

(b)

x

(c)

z

ph˚

x,h˚
yq

θ˚

x

P

y

Q

P

z

Figure 4.3: (a) Ambidexterous juggler standing at P and applying control actions with both hands,
(b) right-handed juggler standing at P and applying control action with right hand, (c) right-handed
juggler standing at Q and applying control action with right hand.

This juggler, whom we will now refer to as the right-handed juggler, To this end, the juggler will
apply the control action standing at P when θ “ θ˚ - see Fig. 4.3 (b), and apply the next control
action after changing location to Q (mirror image of P) and facing the z “ 0 plane when θ “ π´θ˚
- see Fig. 4.3 (c). When the devil stick has the orientation θ “ π´ θ˚, as seen by an observer at
P, it will have the orientation θ “ θ˚ for the right-handed juggler. After applying control action at
Q, the right-handed juggler will return back to P. If xyz denotes the rotating coordinate frame of

the right-handed juggler, the change in position of this juggler can be described by the coordinate

63

transformation:

where

»————–

x

y

z

fiffiffiffiffiflQ

“ Ry,π

»————–

x

y

z

fiffiffiffiffiflP

, »————–

x

y

z

fiffiffiffiffiflP

“ Ry,π

»————–

x

y

z

fiffiffiffiffiflQ

Ry,π fi diagr´1, 1, ´1s

Since Ry,π changes the sign of the x and z coordinates and leaves the y coordinate unchanged, we

can show

where

YQ “ RYP, YP “ RYQ

R “ R´1 fi diagr´1, ´1, ´1, 1, 1s

and YP and YQ denote the vector Y as seen by the right-handed juggler standing at points P and Q,

respectively.

4.3.5 Single Return Map and Discrete-Time Model

In the reference frame of the right-handed juggler, who alternates between positions P and Q, the

two Poincaré sections Sl and Sr are identical, and equal to

S : tX P R

6 | θ “ θ˚u

(4.3.8)

This follows from our discussion in section 4.3.4 as well as Figs. 4.3 (b) and (c). The half-return

maps Pr and Pl in ( 4.3.6) and ( 4.3.7) can be rewritten as follows:

YPpt´
YPpt´

k`1q “ AYPpt´
k`1q “ AYPpt´

k q` Br,
k q` Bl,

k “ 1, 3, 5,¨¨¨
k “ 2, 4, 6,¨¨¨

(4.3.9a)

(4.3.9b)

64

to explicitly indicate the reference frame of Y . Since the right-handed juggler alternates between

positions P and Q, the half-return map Pl in ( 4.3.9b) can be transformed as follows:

RYPpt´
ñ YQpt´
ñ YQpt´

k`1q “ RAYPpt´
k`1q “ ARYPpt´
k`1q “ AYQpt´

k q` RBl
k q` RBl
k q` Br,

k “ 2, 4, 6,¨¨¨

(4.3.10)

where we used the relations RA “ AR and RBl “ Br. It is clear from ( 4.3.9a) and ( 4.3.10) that
the half-return maps Pr and Pl in ( 4.3.6) and ( 4.3.7) are identical in the reference frame of the

right-handed juggler. This implies that the hybrid dynamics of the devil-stick between any two

control actions can be described by a single return map if the change in reference frame of the

right-handed juggler is incorporated in the dynamic model. This map, P : S Ñ S, can be obtained
by first rewriting the left-hand-sides of ( 4.3.6), ( 4.3.9a) and ( 4.3.10) as follows:

k q` Br,

RYQpt´
ñ YQpt´
RYPpt´
ñ YPpt´

k`1q “ AYPpt´
k`1q “ R“AYPpt´
k`1q “ AYQpt´
k`1q “ R“AYQpt´

k “ 1, 3, 5,¨¨¨
k “ 1, 3, 5,¨¨¨
k “ 2, 4, 6,¨¨¨
k “ 2, 4, 6,¨¨¨

(4.3.11a)

(4.3.11b)

k q` Br,

k q` Br‰ ,
k q` Br‰ ,

Then, by accounting for the change in reference frame of the right-handed juggler after each control

action, ( 4.3.11a) and ( 4.3.11b) can be combined into the following single equation which represents

the return map P:

k`1q “ R“A ¯Ypt´
¯Ypt´

k q` Br‰ ,

k “ 1, 2, 3,¨¨¨

65

where ¯Y denotes the state vector Y in the reference frame of the right-handed juggler. The above

equation results in the following discrete-time equations:

k q´pIk rk{Jq

k`1q “ ´ωpt´
k`1q “ ´hxpt´
k`1q “ ´vxpt´
k`1q “ hypt´
k`1q “ vypt´

k q´“vxpt´
k q´pIk{mq sinθ˚‰δk
k q`pIk{mq sinθ˚
k `“vypt´
k q´p1{2qgδ2
k q`pIk{mq cosθ˚´gδk

k q`pIk{mq cosθ˚‰δk

ωpt´
hxpt´
vxpt´
hypt´
vypt´
k`1 ´ t´

fi pt´

where δk
k q, k “ 1, 2,¨¨¨ , is the time of flight between two consecutive control actions.
During this time duration, the devil-stick rotates by a net angle π´ 2θ˚. Since the angular velocity
of the stick remains constant in the interval rt`

k ,t´

k`1s, δk is given as follows
, Δθ fi pπ´ 2θ˚q

(4.3.12a)

(4.3.12b)

(4.3.12c)

(4.3.12d)

(4.3.12e)

(4.3.13)

k q`pIk rk{Jq
The control design for juggling is presented next.

ωpt´

δk “

Δθ

4.4 State Feedback Control Design

4.4.1 Steady-State Dynamics

From the discussion in section 4.3.5 it becomes clear that when the change of reference frame

of the juggler is taken into account, the problem of juggling between the two distinct configura-
tions pθ˚, h˚
configurations pθ˚, h˚
by

yq is converted to the problem of juggling between identical
x, h˚
yq. If the state variables at this configuration are denoted

x, h˚
yq and pθ˚, h˚

x, h˚

x, h˚

yq and pπ´ θ˚,´h˚

¯Y ˚ fi„ω˚ h˚

x

v˚
x h˚
y

v˚

yT

(4.4.1)

66

where I˚

k , r˚

value of the time of flight. Since h˚

k denote the steady-state values of the control inputs and δ˚ denote the steady-state
y is eliminated from ( 4.4.2d), ( 4.4.2) represents six equations
x, v˚
y, I˚, r˚, and δ˚. By choosing δ˚, the remaining six

x, v˚

in seven unknowns, namely, ω˚, h˚

unknowns are obtained as follows:

tanθ˚{4

ω˚ “ ´Δθ{δ˚,
h˚
x “ gδ˚2
v˚
y “ ´gδ˚{2
x “ g tanθ˚δ˚{2, v˚
I˚ “ mgδ˚{ cosθ˚,r˚ “ 2J cosθ˚Δθ{pmgδ˚2q

Since the point of application of the impulsive force must lie on the stick, r˚ in ( 4.4.3) must satisfy
0 ă r˚ ă ℓ{2. This imposes the following constraint of the value of δ˚:

mgℓ
It should be noted that for a given value of δ˚, the value of h˚

y is not unique.

δ˚ ą ¯δ,

¯δ fi 2d J cosθ˚Δθ

(4.4.2a)

(4.4.2b)

(4.4.2c)

(4.4.2d)

(4.4.2e)

(4.4.2f)

(4.4.3)

(4.4.4)

(4.4.5)

then ¯Y ˚ “ Pp ¯Y ˚q is a fixed point and ( 4.3.12) and ( 4.3.13) give

k{mq sinθ˚‰δ˚
y `pI˚

ω˚ “ ´ω˚ ´pI˚
k r˚
k{Jq
x´“v˚
x´pI˚
x “ ´h˚
h˚
x `pI˚
v˚
x “ ´v˚
k{mq sinθ˚
y ´p1{2qgδ˚2 `”v˚
y “ h˚
h˚
y `pI˚
v˚
y “ v˚
k{mq cosθ˚´gδ˚
Δθ
δ˚ “
k r˚
ω˚ `pI˚

k{Jq

k{mq cosθ˚ıδ˚

4.4.2 Error Dynamics

To converge the states to their desired values, we first define the discrete error variables:

k q´ ω˚
k q´ h˚
x,
k q´ h˚
y,

rωpkq fi ωpt´
rhxpkq fi hxpt´
k q´ v˚
rhypkq fi hypt´
k q´ v˚
ru1pkq fi pIkrk ´ I˚r˚q{J, ru2pkq fi pIk ´ I˚q{m

rvxpkq fi vxpt´
rvypkq fi vypt´

x

y

67

Using ( 4.3.12) and ( 4.4.2a)-( 4.4.2e), the error dynamics can now be written as

rωpk` 1q “´rωpkq´ru1pkq
rhxpk` 1q “´rhxpkq´δkrvxpkq` δk sinθ˚ru2pkq
rvxpk` 1q “´rvxpkq` sinθ˚ru2pkq
rhypk` 1q “rhypkq` δkrvypkq` δk cosθ˚ru2pkq`pg{2q”δkδ˚
rvypk` 1q “rvypkq` cosθ˚ru1pkq´grδk ´ δ˚

k s

kı
k ´ δ2

where δk, defined in ( 4.3.13), can be written in terms of the error variables as follows:

(4.4.6a)

(4.4.6b)

(4.4.6c)

(4.4.6d)

(4.4.6e)

(4.4.7)

δk “

Δθδ˚

rrωpkq`ru1pkqsδ˚ ` Δθ

It is clear from ( 4.4.6) and ( 4.4.7) that the error dynamics is nonlinear. In the next section we

present a partial control design that converts the nonlinear system into a linear system and simplifies

the remaining control design.

4.4.3 Partial Control Design: Dead-Beat Control

The error dynamics in ( 4.4.6) involves two control inputs, namely, ru1pkq and ru2pkq. The input
ru1pkq appears only in ( 4.4.6a). To this end, we first designru1pkq as follows:

(4.4.8)

to guarantee dead-beat convergence of the error staterωpkq. Substitution of ( 4.4.8) in ( 4.4.7) yields
δk “ δ˚. Since, δ˚ is user-defined and is a constant, the choice of control in ( 4.4.8) is special as it

ru1pkq “ ´rωpkq

68

transforms the remaining dynamics in ( 4.4.6b)-( 4.4.6e) into the linear system:

zpk` 1q “ A zpkq` Bru2pkq
zpkq fi„rhxpkq rvxpkq rhypkq rvypkqT
»———————–
»———————–
´1 ´δ˚ 0
0 ´1

fiffiffiffiffiffiffiffifl

, B fi

1 δ˚

0

0

0

0

0

1

0

0

0

δ˚ sinθ˚

sinθ˚

δ˚ cosθ˚

cosθ˚

fiffiffiffiffiffiffiffifl

A fi

(4.4.9)

It can be verified that the pair pA , Bq is controllable since θ˚ P p0,π{2q and δ˚ ą 0.

4.4.4 Residual Control Design

The error dynamics in ( 4.4.9) is linear and therefore the states can be converged to zero by simply

designing a linear controller. However, it should be noted that the control inputru2pkq determines
the value of Ik which also appears in the dead-beat control designru1pkq - see ( 4.4.5). By using the
values ofru2pkq from ( 4.4.5) in ( 4.4.8), we get:

(4.4.10)

rk “ rI˚r˚ ´ Jrωpkqs{Ik

Since the point of application of impulsive force must lie of the stick, rk must satisfy´ℓ{2ă rk ă ℓ{2.
By imposing this condition on the value of rk in ( 4.4.10), we get the following constraints on the

inputru2pkq:

(4.4.11)

ru2pkq ą r 2I˚r˚ ´ 2Jrωpkq´ I˚ℓs{pmℓq
ru2pkq ą r´2I˚r˚ ` 2Jrωpkq´ I˚ℓs{pmℓq

Since I˚ and r˚ are both positive, as it can be seen from ( 4.4.3), the inequalities in ( 4.4.11) can be

combined into the single inequality:

ru2pkq ą ¯a` ¯b |rωpkq|
¯a fi p2r˚ ´ ℓqI˚{pmℓq,

¯b fi 2J{pmℓq

69

(4.4.12)

Sinceru1pkq is dead-beat, rωpkq “ 0, k “ 2, 3¨¨¨ . Thus, ( 4.4.12) can also be written as

The input ru2pkq is designed using Linear Quadratic Regulator (LQR) and Model Predictive

Control (MPC) methods. For an LQR design, the control minimizes the cost function

ru2pkq ą ¯a` ¯b |rωpkq|,
ru2pkq ą ¯a,

k “ 1
k “ 2, 3,¨¨¨

J “

2pkqı
8ÿk“1”zpkqT Q zpkq` Rru

2

(4.4.13)

(4.4.14)

(4.4.15)

where, Q and R are constant weighting matrices that can be chosen by trial and error to satisfy the

constraints in ( 4.4.13). The closed-form solution of the control input ru2pkq can be obtained by

solving the Ricatti equation [99].

For a receding horizon MPC design, the constraint in ( 4.4.13) can be explicitly included in

the optimization problem. In the MPC design5, it is necessary to calculate the predicted output

with future control input as the adjusted variable. Since the current control input cannot affect the

output at the same time for receding horizon control, the system dynamics must be represented in

terms of the difference between the current and the predicted control input. To this end, we define

the following variables based on the augmented state-space model6 in [100]:

Δupkiq firu2pkiq´ru2pki ´ 1q

Δupkiq
Δupki ` 1q

...

ΔUi fi

, Zi fi

zpki ` 1 | kiq
zpki ` 2 | kiq

...

Δupki ` Nc ´ 1q

zpki ` Np | kiq

fiffiffiffiffiffiffiffifl

»———————–

»———————–

fiffiffiffiffiffiffiffifl

where ki is the current sampling instant, zpkiq is the state vector in ( 4.4.9) measured at ki, Nc is
the control horizon, Np is the prediction horizon, and zpki` m | kiq is the predicted state variable at
ki ` m with state measurements zpkiq.

5A detailed discussion of MPC design for discrete-time systems can be found in Chapters 1-3

in [100].

6The augmented state-space model is controllable; this was verified using Theorem 1.2 in [100].

70

We now construct the following N-step receding horizon optimal control problem:

Nÿi“1”ZT

i Zi ` ΔU T

i

ΔUiı

minimize J “

subject to

zpki ` 1q “ A zpkiq` Bru2pkiq
ru2pkiq ą ¯a` ¯b |rωp1q|,
ru2pkiq ą ¯a,

i “ 1
i “ 2, 3¨¨¨ , N

(4.4.16)

(4.4.17)

In every sampling period, the optimization problem determines the best control parameter ΔUi that

attempts to converge the sequence of states in Zi to zero. Although ΔUi contains Nc number of future

control inputs, only the first entry is implemented as the actual control input. This optimization

process is repeated using a more recent measurement of the states. It should be emphasized that the

window.

input constraint in ( 4.4.17), namely,ru2pkiq ą ¯a` ¯b |rωp1q| is imposed only in the first optimization
In subsequent optimization windows, the constraint is relaxed to ru2pkiq ą ¯a. This is
necessary forru2 to converge to zero since ¯a is negative - see ( 4.4.12), whereas ¯a` ¯b |rωp1q| can
assume positive values based on the initial value of rω.
Remark 1. The control inputru2pkq is obtained as the numerical solution of the optimal control

problem in ( 4.4.16) and ( 4.4.17). These inputs are applied at discrete time instants and the

optimization solver is required to compute these inputs within the sampling time interval, which is

equal to the time of flight δ˚. Since δ˚ is relatively large, there is sufficient time for the optimization

solver to generate the solution. This, along with the fact that the input constraint can be explicitly

considered in the problem formulation, makes MPC well-suited for this problem.

71

4.5 Simulation Results

4.5.1 System Parameters and Initial Conditions

We present simulation results of both LQR- and MPC-based control designs. The physical param-

eters of the devil-stick are provided below in SI units:

m “ 0.1,

ℓ “ 0.5,

J “ 0.0021

(4.5.1)

Using these physical parameters and by choosing the values of θ˚ “ π{6 rad and δ˚ “ 0.5 sec, the
steady-state values of state variables and control inputs are obtained from ( 4.4.3) as

ω˚ “ ´4.18 rad{s h˚
v˚
y “ ´2.45 m{s

x “ 0.353 m v˚
I˚ “ 0.565 Ns

x “ 1.414 m{s
r˚ “ 0.030 m

Since h˚

y can be chosen arbitrarily, we chose

h˚
y “ 3.0 m

At the initial time, we assume θ “ θ˚ “ π{6 rad and the states variables (in SI units) are

ωp0q “ 0,

hxp0q “ 0.53,

vxp0q “ 2.0,

hyp0q“ 1.0,

vyp0q“ ´2.0

(4.5.2)

(4.5.3)

(4.5.4)

For the physical parameters in ( 4.5.1), steady-state values of the states in ( 4.5.2) and ( 4.5.3), and

initial conditions in ( 4.5.4), the controlru1pkq was chosen according to ( 4.4.8). The control input
ru2pkq was designed using LQR and MPC methods and simulation results are presented next.

4.5.2 Results for the LQR-based Design

For the LQR design, the weight matrix Q for the states was chosen to be the identity matrix and the

control weight R was chosen as 0.2. The control was obtained as

ru2pkq “ Fzpkq, F “„´0.43 ´0.77 0.43 0.66

72

0.0

-5.0

0.6

0.1

3.0

1.0

1

(a)

(c)

(e)

5
k

ωpkq (rad/s)

3.5

1.0

2.2

(b)

(d)

Epkq (J)

hxpkq (m)

vxpkq (m/s)

hypkq (m)

1.0

-0.9

-2.6

1

10

(f)

5
k

vypkq (m/s)

10

Figure 4.4: State variables and total energy E of the devil-stick at sampling instants k, k “
1, 2,¨¨¨ , 10, for the LQR design.
(a)

(b)

0.7

0.5

1

Ik (Ns)

5
k

0.04

0.01

10

1

rk (m)

5
k

10

Figure 4.5: Control inputs for the devil-stick at sampling instants k, k “ 1, 2,¨¨¨ , 10, for the LQR
design.

The simulation results are shown in Figs. 4.47 and 4.5. It can be seen from Fig. 4.4 (a) that

the dead-beat control ru1pkq converges ωpkq to ω˚ in one sampling interval. The control ru2pkq

converges the remaining states to their steady-state values given in ( 4.5.2) in approximately k “ 10
steps - see Figs. 4.4 (c)-(f). The control inputs Ik and rk are shown in Figs. 4.5 (a) and (b). It can

7It should be noted that the state variables are shown in the reference frame of the right-handed

juggler.

73

be seen that both control inputs converge to their steady-state values defined in ( 4.5.2); also the

control input rk remains well inside the constraint boundary | rk |ă ℓ{2. The convergence of both
the states and control inputs to their desired values imply that the devil-stick is juggled between

two symmetric configurations. Since the magnitudes of vx, vy, ωx, and hy are the same in the two

symmetric configurations, the total energy E (kinetic plus potential) reaches a constant value at

steady state - see Fig. 4.4(b).

Remark 2. The total energy of the devil-stick is the same at the symmetric configurations. Also,

it is conserved during the flight phase. Therefore, in steady-state, the control inputs I˚ and r˚ do

zero work on the system.

4.5.3 Results for the MPC-based Design

The control horizon, prediction horizon, and the number of steps were taken as

Nc “ 5, Np “ 10, N “ 15

The MPC problem, defined by ( 4.4.16) and ( 4.4.17) were solved using quadratic programming in

Matlab8. The state variables hxpkq, hypkq, vxpkq and vypkq and the control inputs Ik and rk are shown
in Fig. 4.6. The state variable ωpkq is not shown as it converged to its desired value in one sampling
interval by the dead-beat controller. Similar to the LQR design, the control input rk remains well

inside the constraint boundary. The trajectory of the center-of-mass of the devil stick is shown in

Fig. 4.7 (a); it starts from the initial configuration phx, hyq “ p0.53, 1.00q and is eventually juggled
between the symmetric coordinates ph˚
yq “ p´0.353, 3.00q in
steady state. Typically, N is chosen to be large to guarantee convergence. For our system, the

yq “ p0.353, 3.00q and p´h˚

x, h˚

x, h˚

states rapidly converged to zero with N “ 15. In Fig. 4.7 (b), the devil-stick is shown at the two
symmetric configurations where θ˚ “ π{6 and several intermediate configurations that are equal
time intervals apart.

8The quadprog Matlab function was used.

74

1.0

0.0

3.0

1.0

1.5

(a)

(c)

(e)

hxpkq (m)

hypkq (m)

Ik (Ns)

2.0

1.0

-1.5

-2.6

0.04

(b)

(d)

(f)

vxpkq (m/s)

vypkq (m/s)

rk (m)

0.1
1

5

10

k

0.01

1

15

5

10

15

k

Figure 4.6: State variables and control inputs of the devil-stick at sampling instants k, k “
1, 2,¨¨¨ , 15, for the MPC design.
(a)

(b)

3.0

)

m

(

y
h

1.0

0.5

hx (m)

3.50

)

m

(

y
h

θ˚

2.75

-0.6

1.5

θ˚

0.6

0.0
hx (m)

Figure 4.7: (a) Trajectory of the center-of-mass from the initial configuration to steady-state and (b)
symmetric configurations and seven intermediate configurations of the devil-stick in steady state
for the MPC design.

Remark 3. In both simulations, the stick rotates by an angle pπ´ 2θ˚q between two consecutive

control inputs. This corresponds to “top-only idle" juggling [47]. The controller is quite general

and the stick can be controlled to rotate by pqπ´ 2θ˚q, q “ 2, 3,¨¨¨ , by simply changing the

75

definition of Δθ in ( 4.3.13) from Δθ “ pπ´ 2θ˚q to Δθ “ pqπ´ 2θ˚q. In other words, the stick

can be made to flip multiple times in the flight phase, if desired. The “flip-idle" in [47] corresponds

to the case where q “ 2.

76

CHAPTER 5

ENERGY-BASED ORBITAL STABILIZATION USING CONTINUOUS INPUTS AND

IMPULSIVE BRAKING

5.1 Introduction

We present a controller for energy-based orbital stabilization of the class of Euler-Lagrange

systems that have one passive revolute joint. Examples of systems in this class include the pendubot,

acrobot, cart-pendulum system, rotary pendulum etc. The orbit is defined by fixed positions of

the active coordinates and a desired energy of the overall system. The controller is comprised of

continuous-time inputs and impulsive inputs for braking. At first, a general result for underactuated

systems is presented which shows that an impulsive input causes an instantaneous jump in the

mechanical energy of the system. This jump is shown to be explicitly dependent on the change in

the active velocities due to an impulsive input. This result is then used to show that impulsive braking

causes a negative jump in the mechanical energy of the system as well as the Lyapunov-like function.

Finally, using a state-dependent impulsive dynamical system model [101], we present sufficient

conditions for stabilization. To demonstrate the generality of our approach, we demonstrate orbital

stabilization for the three-DOF Tiptoebot [6] through simulations. Experimental validation of the

control design is carried out on a rotary pendulum to show the applicability of our approach in real

hardware.

5.2 Problem Statement

Consider an n degree-of-freedom underactuated system with one passive degree-of-freedom.

The generalized coordinates of the system are denoted by q, q fi“qT

1 q2‰T , where q1 P Rn´1 and

q2 P R are the coordinates associated with the active and passive degrees-of-freedom. Our control
objective is to stabilize the orbit defined by

pq1, 9q1, Eq “ p0, 0, Edesq

(5.2.1)

77

where E is the total mechanical energy of the system and is given by the relation

Epq, 9qq “

1

2

9qT Mpqq 9q` Fpqq

(5.2.2)

and Edes is the desired value of E. In (5.2.2), Mpqq P Rnˆn is the mass matrix, assumed to be
positive definite, and Fpqq is the potential energy, assumed to be a smooth function of q. The mass
matrix is partitioned as

Mpqq “»—–

M11pqq M12pqq
MT
12pqq M22pqq

fiffifl

(5.2.3)

where M11 P Rpn´1qˆpn´1q and M22 P R.
Assumption 1. The mechanical energy of the system is assumed to be periodic in the passive

coordinate q2, such that Epq2 ` 2πq “ Epq2q.

Remark 4. Assumption 4 is easily satisfied if the passive degree-of-freedom is a revolute joint.

Assumption 2. The elements of the mass matrix Mpqq are bounded and the potential energy Fpqq
is lower bounded. Furthermore, if q1 “ q˚
dq2 rkP1pq˚

1, q2q` P2pq˚

1 is constant, then

1, q2qs ­“ 0

d

where k is any scalar constant, and

P1pq˚

1, q2q “

2

M22˜„BM12
Bq2 ´

P2pq˚

1, q2q “ ´P1pq˚

1, q2qF ´

1

´

M12

Bq2 ¸
Bq1 T
2 M22„BM22
2„BM22
Bq2 `„BF
Bq1T

BF

M12

M22

Remark 5. The boundedness property of Mpqq and Fpqq is satisfied for systems that have no

prismatic joints. For a given underactuated system with one passive DOF, assumption 5 can be

easily verified.

78

5.3 Modeling of System Dynamics

5.3.1 Lagrangian Dynamics

For our system described in section 5.2, the equations of motion can be written as:

d

dtˆBL
dtˆBL

B 9q1˙´ˆBL
B 9q2˙´ˆBL

Bq1˙ “ u
Bq2˙ “ 0

d

(5.3.1)

where L pq, 9qq is the Lagrangian and u P Rn´1 is the vector of independent control inputs. Each
element of the vector u is a continuous function of time for all t ě 0 except at discrete instants
t “ τk, k “ 1, 2,¨¨¨ , where it is impulsive in nature. At t “ τk, the impulsive input vector has the
form upτkq “ Ik δpt ´ τkq, where δpt ´ τkq is the Dirac measure at time τk and Ik P Rn´1 is the
impulse of the impulsive input. The Lagrangian has the form

L pq, 9qq “

1

2

9qT Mpqq 9q´ Fpqq

By substituting (5.2.3) in (5.3.2), the Lagrangian is written as

L pq, 9qq “

1

2

9qT
1 M11 9q1 `

1

2

2
M22 9q

2 ` 9qT

1 M12 9q2 ´ F

and by substituting (5.3.3) in (5.3.1), the equations of motion are written in the form:

M11 :q1 ` M12 :q2 ` h1pq, 9qq “ u
MT
:q1 ` M22 :q2 ` h2pq, 9qq “ 0

12

(5.3.2)

(5.3.3)

(5.3.4a)

(5.3.4b)

where

h1 “ 9M11 9q1 ` 9M12 9q2 ´

h2 “ 9M22 9q2 ` 9qT

1

9M12 ´

1

Bq1pM11 9q1q 9q1 ´„BpM12 9q2q
2„ B
Bq2  9q
2„BM22
1„BpM11 9q1q

´

Bq1

Bq2

9qT

1

1

2

1

2
9q

Bq1 T
 9q1 ´
2„BM22
1„BpM12 9q2q

Bq1T
2 `„BF
` BF

Bq1

Bq2

2

2 ´ 9qT

(5.3.5a)

(5.3.5b)

79

Equations (5.3.4a) and (5.3.4b) can be rewritten in the form

:q1 “ Apq, 9qq` Bpqqu

:q2 “ ´p1{M22q”MT

12tApq, 9qq` Bpqquu` h2ı

where

Bpqq “”M11 ´p1{M22qM12 MT

12ı´1

Apq, 9qq “ p1{M22qBpqqrM12 h2 ´ h1M22s

(5.3.6a)

(5.3.6b)

(5.3.7)

(5.3.8)

Using the properties of the mass matrix Mpqq and the Schur complement theorem [102], it can be
shown that Bpqq is well-defined, symmetric and positive-definite, i.e., Bpqq “ BTpqq ą 0.

5.3.2 Effect of Impulsive Inputs

When the control input u in (5.3.4a) is impulsive, it causes discontinuous jumps in the velocities

p 9q1, 9q2q, while the coordinates pq1, q2q remain unchanged. For the impulsive input at t “ τk, the
jump in the velocities is computed by integrating (5.3.4) as follows [86]:

»—–

M11 M12

MT

12 M22

In the above equation, Δ 9q1 and Δ 9q2 are defined as

fiffifl

»—–

Δ 9q1

Δ 9q2

fiffifl “»—–

Ik

0

fiffifl , Ik

fiż t`

k
t´
k

uptkq dt

Δ 9q1 fi p 9q`

1 ´ 9q´
1 q,

Δ 9q2 fi p 9q`

2 ´ 9q´
2 q

where 9q´ fi 9qpτ´
k q denote the generalized velocities immediately before and after
application of the impulsive inputs. Since the system is underactuated, the jump in 9q2 is dependent

k q and 9q` fi 9qpτ`

on the jumps in 9q1; this dependence is described by the one-dimensional impulse manifold [23] or

impulse line, obtained from the equation above:

9q`
2 “ 9q´

2 ´p1{M22qMT

12p 9q`

1 ´ 9q´
1 q

(5.3.9)

80

The kinetic energy undergoes an instantaneous change due to jumps in the generalized velocities.

This change is also equal to the change in the total mechanical energy of the system since the

potential energy is only a function of the generalized coordinates. A formal statement of this result

is provided next.

Lemma 1. For the dynamical system in (5.3.4), the jump in the total mechanical energy due to

application of an impulsive input is given by

ΔE fi pE` ´ E´q “

1

2

9q`T
1 B´1pqq 9q`
1 ´

1

2

9q´T
1 B´1pqq 9q´

1

(5.3.10)

where E´ and E` are the energies immediately before and after application of the impulsive input.

Proof: See section 5.7.1.

Remark 6. The proof of Lemma 1 is provided for the general case where the number of active and

passive degrees-of-freedom are pn´ mq and m, respectively. This general result indicates that the
change in mechanical energy due to an impulsive input depends only on the velocities of the active

degrees-of-freedom immediately before and after application of the input. The result is analogous

to the passivity property for the continuous-time case [103], where the power input to the system

is the inner product of the velocities of the active degrees-of-freedom and control inputs.

It is

important to note that results similar to Lemma 1 appeared earlier in [104].

We define impulsive braking as the process of applying impulsive inputs that result in 9q`

1 “ 0.
It follows from Lemma 1 that impulsive braking results in a loss of mechanical energy, given by

the expression

ΔE “ ´

1

2

9q´T
1 B´1pqq 9q´

1

(5.3.11)

The discontinuous jumps in the generalized velocities and kinetic energy of the system also

occurs when mechanical systems experience forces due to impact. In this regard, results similar to

the ones presented above can be found in [104–106].

We now state an important result related to impulsive braking.

81

Lemma 2. Consider the scalar function

V “

1

2”qT
1 Kp q1 ` 9qT

1 Kd

9q1 ` KepE ´ Edesq2ı

(5.3.12)

where Kp and Kd are diagonal positive definite constant matrices and Ke is a positive constant.

Impulsive braking results in a discontinuous jump in the function given by

ΔV fi pV ` ´V ´q “ ´

1

2

9q´T

1 „ 1

4"Ke 9q´T

1 B´1pqq 9q´

1* B´1pqq` Kd ` KepE` ´ EdesqB´1pqq 9q´

1

(5.3.13)

where V ´ and V ` are values of the function immediately before and after impulsive braking.

Furthermore, if”Kd ` KepE` ´ EdesqB´1pqqı is positive definite, then ΔV ď 0, and ΔV “ 0 if and

only if 9q´

1 “ 0.

Proof: See section 5.7.2.

5.3.3 Hybrid Dynamical Model

For orbital stabilization of our system in (5.3.4), we propose a control strategy comprised of

continuous and impulsive inputs. The impulsive inputs will be used for impulsive braking of the

active coordinates, i.e.,

9q`
1 “ 0. As a result, the change in the velocities can be obtained using

(5.3.9) as follows:

Δ 9q1 “ 0´ 9q´
2 ´ 9q´
Δ 9q2 “ 9q`

1 “ ´ 9q´
2 “ p1{M22qMT

1

12

9q´

1

(5.3.14)

In addition to the impulsive braking inputs, we will reset the passive coordinate q2 periodically to

confine it to the compact set r´3π{2,π{2s. To describe the dynamics of our system, we adopt the
state-dependent impulsive dynamical model in [101, pg.20]:

9xptq “ fcrxptqs,
Δxptq “ fdrxptqs,

xp0q “ x0,

xptq R Z
xptq P Z

(5.3.15a)

(5.3.15b)

82

where Z defines the set where the impulsive inputs are applied and/or periodic resetting occurs.

For our system,

1

q2

9qT
1

xptq fi rqT
Δxptq fi xpt`q´ xpt´q

9q2sT P D Ď R

2n

In the above expression, D is the open set where q2 P pa, bq, a ă ´3π{2, b ą π{2, and xpt´q, xpt`q
are the values of the state variables immediately before and after application of impulsive inputs or

coordinate resetting. Using (5.3.6), (5.3.9) and (5.3.14), it can be shown

fiffiffiffiffiffiffiffifl

: xptq P Z1

: xptq P Z2

: xptq P Z3

(5.3.16)

(5.3.17)

fc “

fd “

»———————–
$’’’’’’’&’’’’’’’%

9q1

9q2

Apq, 9qq` Bpqqu

12tApq, 9qq` Bpqquu` h2‰
1T
p1{M22qMT

9q´

12

1

´p1{M22q“MT
„0 0 ´ 9q´
„0 2π 0 0T
„0 ´2π 0 0T

Z “ Z1 Y Z2 Y Z3, Z1 is the set where impulsive braking inputs are applied (to be defined in
Theorem 2 where the control design will be presented), and Z2 fi tq2 “ ´3π{2, 9q2 ă 0u and
Z3 fi tq2 “ π{2, 9q2 ą 0u are the sets where coordinate resetting occurs.

Note that the solution xptq of (5.3.15) is left-continuous, i.e., it is continuous everywhere except
at time instants tk, k “ 1, 2,¨¨¨ , where tk’s denote the time instants when impulsive inputs are
applied or coordinate resetting occurs. Since impulsive inputs are applied at τk, k “ 1, 2,¨¨¨ ,
tτ1,τ2,¨¨¨u Ă tt1,t2,¨¨¨u. Thus, xpt`q and xpt´q, in the equation after (5.3.15), were defined

83

without the time stamp, i.e.

xpt´q fi xpt´
xpt`q fi xpt´

k q “ lim
εÑ0`
k q` fdrxpt´

xptk ´ εq

k qs “ lim
εÑ0`

xptk ` εq

The well-posedness of the tk’s and the quasi-continuous dependence of the solution to (5.3.15)

with respect to initial conditions [101, 104, 107] will be established later in the proof of Theorem 4

where Z1 will be defined.

5.4 Hybrid Control Design

5.4.1 Main Result

For the control objective in (5.2.1), we propose a hybrid control strategy comprised of a continuous

controller and impulsive braking1. Theorem 4 provides the design of the continuous controller

and defines the set Z1, where impulsive braking is applied. The proof of Theorem 2 is based on

a Lyapunov-like function. The continuous controller is invoked as long as the derivative of the

Lyapunov-like function is negative semi-definite; when this condition is not satisfied, impulsive

braking is applied to produce negative jumps in the Lyapunov-like function. Before stating Theorem

4, we present Lemma 3 and an invariant set theorem [101, pg.38] that will be used in the proof of

Theorem 4.

Lemma 3. For the system in (5.3.4) subjected to continuous control, q2 is constant if 9q1 and 9u are

identically zero.

Proof: See section 5.7.3.

Remark 7. Lemma 3 implies that the active and passive generalized coordinates are dynamically

coupled. Due to this coupling, the active generalized coordinates cannot be held stationary by

1It is assumed that the active degrees-of-freedom will have a friction brake such that they can

be stopped instantaneously.

84

constant generalized forces when the passive generalized coordinate is non-stationary. The proof

of Lemma 3 is based on assumption 5, which is satisfied when dynamical coupling exists. The

existence of such coupling has been verified for an inverted pendulum on a cart [74], pendubot,

acrobot, and reaction-wheel pendulum; in this work it is shown for the tiptoebot and the rotary

pendulum.

Theorem 3. [101, pg.38] Consider the impulsive dynamical system given by (5.3.15), assume that

Dc Ă D is a compact positively invariant set with respect (5.3.15), and assume that there exists a
continuously differentiable function W : D Ñ R such that

rBWpxq{Bxs fcpxq ď 0,
Wpx` fdpxqq

ď Wpxq,

x P Dc,
x P Dc,

x R Z
x P Z

(5.4.1a)

(5.4.1b)

Let R fi tx P Dc : x R Z , rBWpxq{Bxs fcpxq “ 0uYtx P Dc : x P Z ,Wpx` fdpxqq “ Wpxqu and let
M denote the largest invariant set contained in R. If x0 P Dc, then xptq Ñ M as t Ñ 8.

Remark 8. Although the invariance principle in [101] will be used (Theorem 3), the reader should

be aware of the invariance principles in [55, 106, 108] that can be used for analysis of non-smooth

and hybrid systems.

Theorem 4. For the impulsive dynamical system defined by (5.3.15), (5.3.16) and (5.3.17), and

x0 P D such that

´3π{2 ă q2p0q ă π{2

the following choice of control design:

u “ ´rpKd ` Kcq Bpqq` KepE ´ Edesq Is´1“Kp q1 `pKd ` KcqApq, 9qq‰
Z1 “ txptq | rApq, 9qq` BpqqusT Kc 9q1 ď 0, 9q1 ­“ 0u

(5.4.2)

(5.4.3a)

(5.4.3b)

where I is the identity matrix and Kc is a diagonal positive-definite matrix, guarantees asymptotic

stability of the orbit in (5.2.1) if the gain matrices Kp, Kd and Ke satisfy the following conditions:

85

(i) ”Kd ` KepE ´ EdesqB´1pqqı is positive definite for all q and 9q,

(ii) If q˚

1 and q˚

2 are constant values of q1 and q2, then the following system of equations:

Kp q˚

1

KerFpq˚q´ Edess

„BF
Bq1T
q“q˚ “ ´
Bq2q“q˚ “ 0
„BF

yields a finite number of solutions with q˚

1 “ 0, and

(iii) For all possible solutions of q˚

2 obtained from (ii) and for the function V in (5.3.12), the

following inequality is satisfied

Vpt “ 0q ă mintV | q1 “ 0, 9q1 “ 0, E P SEztEdesuu

where SE is the set of values of E evaluated at q1 “ 0, q2 “ q˚

2, 9q “ 0.

Proof:

Consider the Lyapunov-like function V defined in (5.3.12); V is zero on the orbit defined in

(5.2.1) and positive everywhere else. The time derivative of V is

1 Kd

9V “ qT
1 Kp 9q1 ` :qT
“”qT
1 Kp ` :qT

9q1 ` KepE ´ Edesq 9E
1 Kd ` KepE ´ Edesq uTı 9q1

(5.4.4)

where 9E “ uT 9q1 follows from the passivity property of underactuated Euler-Lagrange systems -
see [103]2 and proposition 2.5 of [109]. By substituting :q1 from (5.3.6a) in (5.4.4) and using the

symmetry of Bpqq, we get

9V “”qT

1 Kp ` AT Kd ` uT B!Kd ` KepE ´ EdesqB´1)ı 9q1

(5.4.5)

The following choice of u

uT “ ´”qT

1 Kp ` AT Kd ` :qT

1 Kcı”B!Kd ` KepE ´ EdesqB´1)ı´1

,

(5.4.6)

2The proof of the passivity property follows from the fact that the matrix [ 9M ´ 2C] is skew-

symmetric for our choice of generalized coordinates.

86

which is well defined based on condition (i), results in

9V “ ´ :qT

1 Kc 9q1

(5.4.7)

Substitution of (5.3.6a) in (5.4.6) followed by algebraic manipulation gives the expression for u in

(5.4.3a). Substitution of (5.3.6a) in (5.4.7) gives

9V “ ´rApq, 9qq` BpqqusT Kc 9q1

(5.4.8)

Based on the expression of 9V , three cases may arise:

case (a): if rA` BusT Kc 9q1 ą 0, then 9V ă 0,
case (b): if rA` BusT Kc 9q1 ď 0, 9q1 ­“ 0, then x P Z1 and impulsive braking is applied - see (5.4.3b).
Since condition (i) is satisfied, Lemma 2 indicates that V undergoes a discontinuous

change ΔV , where ΔV ă 0, and

case (c): if 9q1 “ 0, then 9V “ 0.
For case (b), impulsive braking results in 9q1 “ 0 at t` and the trajectories of the system leave Z1.
If 9q1 ” 0 for all t ą t`, the trajectories of the system remain outside Z1 and 9V ” 0. If 9q1 ı 0 for
t ą t`, V decreases and the trajectories of the system move away from Z1 since (5.4.7) implies

9Vpt`q “ 0,

:Vpt`q “ ´ :qT

1 pt`qKc :q1pt`q ă 0

ñ 9Vpτq ă 0, τ P pt`,t` ` εq

for some ε ą 0 since Kc is positive-definite and :q1 ­“ 0. The trajectories may re-enter Z1, but
not arbitrarily quickly. Hence Zeno phenomenon will not be observed; this is discussed in section
5.7.4. Case (c) implies that either 9q1 ” 0 ñ 9V ” 0, or 9q1 ı 0 and V continues to decrease again;
this follows from our discussion of the nature of trajectories after impulsive braking. Cases (a), (b)

and (c) imply that for t ą 0, Vptq ď Vp0q fi c and therefore the set

Dc fi tV ď cuXt´3π{2 ď q2 ď π{2u

87

is positively invariant.

Cases (a), (b) and (c) together satisfy the conditions in Theorem 33 with Dc defined above and
Wpxq “ Vpxq. Since (b) implies ΔV ă 0, tx P Dc : x P Z , ΔV “ 0u is an empty set. Therefore,
xptq Ñ M Ă R “ tx P Dc : x R Z , 9V “ 0u as t Ñ 8. From case (c), 9V “ 0 implies 9q1 “ 0 and thus
R “ tx P Dc : 9q1 ” 0u. In R, :q1 “ 0. Substitution of :q1 “ 0 in (5.3.6a) and (5.4.6) yields

uT “ ´AT B´1

uT BKd “ ´KepE ´ EdesquT ´ qT

1 Kp ´ AT Kd

Substitution of (5.4.9a) into (5.4.9b) gives

uT KepE ´ Edesq` qT

1 Kp “ 0

(5.4.9a)

(5.4.9b)

(5.4.10)

The definition of R in Theorem 3 implies V is constant in R. Also, q1 is constant and 9q1 “ 0 in
R. Therefore, from the definition of V in (5.3.12), we can claim that E is constant in R. Let q˚

1

and E˚ be the constant values of q1 and E. We now discuss two cases that can arise:

case 1: If E˚ “ Edes, we have q˚
case 2: if E˚ ­“ Edes, we get from (5.4.10)

1 “ 0 from (5.4.10). This implies that M is the orbit in (5.2.1).

u fi u˚ “ ´

Kp q˚

1

KepE˚ ´ Edesq

(5.4.11)

where u˚ is the constant value of the continuous control in R.

For case 2, both q1 and u are constants. Therefore, based on Lemma 3, we claim q2 “ q˚
2 is a
constant. It follows from (5.2.2) that E˚ “ Fpq˚q. Using (5.3.4) and (5.3.5), we can show that the
trajectories in R satisfy

Bq1T
„BF

q“q˚ “ u˚,

„BF
Bq2q“q˚ “ 0

3Before we can apply Theorem 3, we must ensure that the time instants tk, k “ 1, 2,¨¨¨ , are
well-posed and the hybrid dynamical system in (5.3.15) satisfies the quasi-continuous dependence
property. The conditions for well-posedness and its proof appear in section 5.7.5. A sufficient
condition for satisfying the quasi-continuous dependence property and its proof appear in section
5.7.6.

88

Substituting the expression for u˚ from (5.4.11) in the above equation along with E˚ “ Fpq˚q, we
can use condition (ii) to claim q˚
1 “ 0. Using (5.3.12) and cases (a) and (b), we can claim that as
t Ñ 8, V Ñ V ˚, where

V ˚ “

1

2

KepE˚ ´ Edesq2 ď Vpt “ 0q

where E˚ P SE. Since V ˚ ď Vpt “ 0q, we can claim using condition (iii) that E˚ “ Edes, i.e., V ˚ “ 0.
Thus the largest invariant set M is the orbit defined in (5.2.1). This concludes the proof.

(cid:3)

5.4.2 Choice of Controller Gains

It can be easily shown that condition (i) in Theorem 4 is satisfied if

p1{KeqλminpKdq ą rEdes ´ minpFqsλmaxrB´1pqqs

where λminpKdq and λmaxrB´1pqqs are the minimum and maximum eigenvalues of Kd andrB´1pqqs.
Assumption 5 implies λmaxrB´1pqqs and minpFq exist and therefore Kd and Ke can always be
chosen to satisfy condition (i).

For the choice of Ke satisfying condition (i), Kp has to be chosen to satisfy condition (ii).

Although we do not prove that condition (ii) can be simultaneously satisfied for the general case,

several combinations of gains pKp, Kd, Keq were found to exist for the inverted pendulum on
a cart [74]. The authors have independently verified that condition (ii) can be easily satisfied

for several other underactuated mechanical systems, namely, the pendubot, the acrobot, and the

reaction-wheel pendulum. It is shown that conditions (i) and (ii) can be simultaneously satisfied

for the three-DOF Tiptoebot and the rotary pendulum. These examples indicate that condition (ii)

is not restrictive.

Once the controller gains Kp, Kd and Ke have been chosen to satisfy conditions (i) and (ii) in

Theorem 4, condition (iii) imposes no additional restrictions on the gains but simply provides an

estimate of the region of attraction of the orbit. Since Kc does not appear in conditions (i)-(iii), it

can be chosen without restriction.

89

5.5 Illustrative Example - The Tiptoebot

We consider the tiptoebot example as presented in 2.3.1. Using the following definition for the

joint angles

qT
1 “ r θ2 θ3 sT ,

q2 “ θ1

(5.5.1)

the dynamics of the tiptoebot can be expressed in the form of (5.3.4); the components of mass

matrix in (5.2.3) are

M11 “»—–
M12 “»—–

fiffifl

α2`α3`2α5 cosθ3 α3`α5 cosθ3

α3`α5 cosθ3

α3

α2`α3`α4 cosθ2`2α5 cosθ3`α6 cospθ2`θ3q

α3`α5 cosθ3`α6 cospθ2`θ3q

fiffifl

(5.5.2)

(5.5.3)

M22 “ α1`α2`α3 ` 2rα4 cosθ2 ` α5 cosθ3`α6 cospθ2`θ3qs

where αi, i “ 1, 2,¨¨¨ , 6, are lumped parameters, defined as follows:

α1 fi m1pℓ2

2 ` ℓ2
α4 fi m2ℓ1ℓ2 ` m3ℓ1ℓ2

1 ` ℓ2

3q, α2 fi pm2 ` m3qℓ2

2, α3 fi m3ℓ2
3,

α5 fi m3ℓ2ℓ3,

α6 fi m3ℓ1ℓ3

The sum of Coriolis, centrifugal and gravitational force terms, h1 and h2, can be obtained using

(5.3.5), where Fpqq has the expression

F “ β1 sinθ1 ` β2 sinpθ1 ` θ2q` β3 sinpθ1 ` θ2 ` θ3q
β1 fi pm1 ` m2 ` m3qℓ1 g, β2 fi pm2 ` m3qℓ2 g, β3 fi m3ℓ3 g

(5.5.4)

The control input is defined as u “ rτ2 τ3sT . In the compact set θ1 P r´3π{2,π{2s, as defined in
section 5.3.3, the upright equilibrium configuration of the tiptoebot is defined by

θ1 “ ´3π{2 or π{2,„ θ2 θ3

9θ1

9θ2

9θ3  “„ 0 0 0 0 0 

is unstable, but can be stabilized, by a linear controller, for example. The stabilized equilibrium

will typically have a finite region of attraction; therefore, to stabilize from an arbitrary initial

90

configurations, we first use the controller in section 5.3 to stabilize an orbit that intersects the region

of attraction. The obvious choice for such an orbit is the one where Edes equals the potential energy
of the system at the equilibrium. Substitution of θ1 “ ´3π{2 or π{2 and θ2 “ θ3 “ 0 in (5.5.4)
yields Edes “ β1 ` β2 ` β3. The control objective in (5.2.1) can therefore be written as

θ2 “ θ3 “ 0,

9θ2 “ 9θ3 “ 0, Edes “ pβ1 ` β2 ` β3q

(5.5.5)

The feasibility of our control design is discussed next.

5.5.1 Selection of Controller Gains

The initial configuration of the tiptoebot is taken as

r θ1 θ2 θ3

9θ1

9θ2

9θ3 s “ r 0 π π 0 0 0 s

(5.5.6)

In this configuration, the tiptoebot is coiled up: the first link is horizontal, the second link folds

back on the first link, and the third link folds back on the second link. The links were chosen to

have the same mass m1 “ m2 “ m3 “ 0.1 kg and the same length ℓ1 “ ℓ2 “ ℓ3 “ 0.6 m. For this
choice of mass and length parameters, the lumped parameters of the tiptoebot, defined in (5.5.3)

and (5.5.4), are provided in Table 6.1 below:

Table 5.1: Tiptoebot lumped parameters in SI units

α1

0.108 α4

0.072 β1

1.764

α2

0.072 α5

0.036 β2

1.176

α3

0.036 α6

0.036 β3

0.588

The passive joint of the tiptoebot is revolute and therefore assumption 1 holds good. Assumption

2 also holds good - this is discussed in 5.7.7.

The following choice of gains satisfy condition (i) and (ii):

Kp “»—–

70

0

0

70

fiffifl , Kd “»—–

91

2.8

0

0

2.8

fiffifl , Ke “ 2.2

(5.5.7)

Condition (ii) results in θ˚

2 “ θ˚

3 “ 0, which upon substitution in (5.3.4b) and (5.3.5b) yields

BF
Bq2 “ 0 ñ cosθ˚

1 “ 0

(5.5.8)

From section 5.3.3 we know that q2 lies in the compact set r´3π{2,π{2s. Thus θ1 lies in the same
compact set - see (5.5.1). In this set, the possible solutions of (5.5.8) are θ˚
1 “ t´3π{2,´π{2,π{2u.
3 “ 0, we know that E “ Edes. Therefore, to satisfy condition

1 “ ´3π{2 or π{2, and θ˚

2 “ θ˚

For θ˚

(iii), we use θ˚

1 “ ´π{2; this results in the following inequality
Vpt “ 0q ă 2KerEpq˚

1 “ 0, q˚

2 “ ´π{2q´ Edess2 “ 2Kepβ1 ` β2 ` β3q2

For the initial configuration in (5.5.6), Ke in (5.5.7) satisfies the inequality above. The matrix Kc

was chosen as

Kc “»—–

1.2

0

0

1.2

fiffifl

(5.5.9)

5.5.2 Simulation Results

For the initial configuration in (5.5.6) and controller gains in (5.5.7) and (5.5.9), the simulation

results are shown in Figs.5.1 and 5.2. The effect of impulsive braking can be seen in Figs.5.1 (d)

and (f), where 9θ2 and 9θ3 (the velocities of the active joints) jump to zero on multiple occasions.

Each impulsive braking also results in a negative jump in the mechanical energy (follows from

Lemma 1) which can be seen in Fig.5.1 (b). Since impulsive inputs cause no jumps in the joint

angles, there is no change in θ1, θ2 and θ3 at the time of impulsive braking - see Figs.5.1 (a), (c)

and (e). In Fig.5.1 (a), θ1 never leaves the set r´3π{2,π{2s and therefore virtual impulsive inputs
are not applied.

While impulsive brakings cause negative jumps in the total energy E, the continuous-time

controller in (5.4.3a) adds energy to the system; together, they converge the energy to the desired

values Edes - see Fig.5.1 (b). The phase portrait of the passive joint is shown in Fig.5.2 (a). The

jumps in the phase portrait (vertical drops in 9θ1, twice) is due to impulsive braking. The variation

92

1.57

-1.57

-4.71

3.0

0.0

3.0

0.0

0

(a)

(c)

(e)

t (s)

θ1 (rad)

θ2 (rad)

10.0

0.0

-3.5

0.0

-14.0

0.0

(b)

pE ´ Edesq (J)

(d)

(f)

9θ2 (rad/s)

θ3 (rad)

9θ3 (rad/s)

-10.0

6

0

t (s)

6

Figure 5.1: Plots of the joint angles θ1, θ2, θ3, error in the desired energy pE ´ Edesq, and the
active joint velocities 9θ2, 9θ3 of the Tiptoebot.

(a)

(b)

700

9θ1 vs θ1

V

9

0

-5

-4.71

-1.57

1.57

0

0

t (s)

6

Figure 5.2: Plots showing (a) phase portrait of passive joint angle θ1, and (b) variation of the
Lyapunov-like function V . The desired orbit is shown using dashed green line in (a).

of the Lyapunov-like function V with time is shown in Fig.5.2 (b) - it can be seen that V decreases

monotonically due to the action of the continuous-time controller and undergoes negative jumps

intermittently due to impulsive brakings. The continuous controller and impulsive brakings work

together to converge V to zero.

93

The gain matrices in (5.5.7) and (5.5.9) were chosen such that convergence to the desired orbit

is fast. The simulation results indicate that the system trajectories reach a close neighborhood of

the desired orbit very quickly, at approximately 3 s. For stabilization of the equilibrium in (5.5.6),

a linear controller was designed using LQR. The matrices Q and R of the algebraic Ricatti equation

were chosen to be I6ˆ6 and 2I2ˆ2, where Ikˆk is the identity matrix of size k. The linear controller
was invoked when V ď 0.05 and | θ1 ´ π{2 |ď 0.05.
5.6 Experimental Validation

5.6.1 System Description

Experiments were done with a rotary pendulum. As shown in Fig.5.3, the system is comprised of a

horizontal arm OA of mass ma and length ℓa, which rotates about point O, and a pendulum of mass

mp and length ℓp, that rotates about point A. The center-of-mass of the horizontal arm is located

at a distance da from O and the center-of-mass of the pendulum is located at a distance dp from

A. The horizontal arm is actively controlled by an external torque τ and its angular displacement

about the z axis is denoted by φ. The pendulum is passive and its angular displacement about the εr

axis is denoted by θ. The accleration due to gravity is denoted by g. With the following definition:

rq1 q2sT “ rφ θsT

(5.6.1)

z

O
φ

g

τ

x

y

ℓa

B

ℓp

θ

εθ

A

εr

Figure 5.3: Schematic of a rotary pendulum.

94

the dynamics of the system can be expressed in the from given by (5.3.4), where u “ τ, and

2 θ, M12 “ γ3 sinθ, M22 “ γ2

M11 “ γ1 ` γ2 cos
h1 “ γ3 cosθ 9θ2 ´ 9φ 9θγ2 sin 2θ,
a, γ2 fi mpd
γ1 fi mad

2

a ` mp ℓ2

h2

2
p, γ3

9φ2

“ γ2
sinθ cosθ` γ4 cosθ
fi ´mp ℓadp, γ4 fi mpgdp

(5.6.2)

The physical parameters of the experimental setup are

γ1 “ 0.0120, γ2 “ 0.0042, γ3 “ ´0.0038, γ4 “ 0.1190

(5.6.3)

The control torque was applied by a 24-Volt permanent magnet brushed DC motor4.The motor is

driven by a power amplifier5 operating in current mode. The motor torque constant is 37.7 mNm/A

and the amplifier gain is 4.4 A/volt. An electromagnetic friction brake6 was integrated to the shaft

of the DC motor. In the OFF state, the brake engages a friction pad to the shaft of the motor which

prevents the shaft from turning; in the ON state, the brake is disengaged and the motor shaft rotates

freely. For impulsive braking, the brake was kept engaged till the active velocity 9φ reached a close

neighborhood of zero. The brake was powered ON/OFF by sending command voltage signals

through an n-channel mosfet. The rotary pendulum was interfaced with a dSpace DS1104 board

and the Matlab/Simulink environment was used for real-time data acquisition and control with a

sampling rate of 1 Khz. The angular positions of the links were measured using incremental optical

encoders; the angular velocities were obtained by differentiating and low-pass filtering the position

signals.

5.6.2 Selection of Controller Gains

The total energy of the system is obtained from (5.2.2) as follows

4The motor manufacturer is Faulhaber Drive Systems. The motor has a gearbox with a reduction

ratio of 3.71 : 1.

5The amplifier is a product of Advanced Motion Control.
6The electromagnetic brake is manufactured by Anaheim Automation, model BRK-20H-480-

024. The brake can withhold torques up to 3.4 Nm.

95

1

2pγ1`γ2 cos

E “
F “ γ4 sinθ

2 θq 9φ2 `

1

2

γ2

9θ2 ` γ3 sinθ 9φ 9θ` F

(5.6.4)

For the control objective in (5.2.1), we choose Edes to be equal to the energy associated with the

homoclinic orbit that contains the upright equilibrium

„ φ θ 9φ 9θ  “„ 0 π{2 0 0  or „ 0 ´3π{2 0 0 

Using (5.6.4), the energy associated with the homoclinic orbit can be written as

Edes “ γ4

(5.6.5)

The passive joint of the rotary pendulum is revolute and thus assumption 1 holds good. Assumption

5 also holds good - this is shown in 5.7.7. From (5.6.4) we know that F is only a function
of θ and therefore condition (ii) is trivially satisfied resulting in the solution φ˚ “ 0.
In the
compact set r´3π{2,π{2s, the possible solutions of θ˚ obtained from condition (ii) are θ˚ “
t´3π{2,´π{2,π{2u. At θ˚ “ π{2 or θ˚ “ ´3π{2 and φ˚ “ 0, E “ Edes. Using condition (iii),
we therefore get θ˚ “ ´π{2; this implies that Ke should be chosen to satisfy

Vpt “ 0q ă 2Keγ2

4

(5.6.6)

At the lower equilibrium configuration where rφ θ 9φ 9θs “ r0 ´π{2 0 0s, we have V “ 2Keγ2
violates the inequality in (5.6.6). This implies that our controller cannot swing-up the pendulum

4 . This

when the system is exactly at the lower equilibrium.Therefore, in experiments, a small external

perturbation was provided such that the system is not at the lower equilibrium at the initial time.

For the experimental results presented herein, the initial configuration of the system after the

perturbation was measured as

„φp0q θp0q

9φp0q

9θp0qT

“„0.01 ´1.42 0.05 0T

(5.6.7)

For the initial conditions in (5.6.7) and physical parameter values in (5.6.3), the following gains

satisfied conditions (i)-(iii):

Kp “ 0.5, Kc “ 0.08, Kd “ 0.3, Ke “ 100

(5.6.8)

96

For the above set of gains, the experimental results are presented next.

5.6.3 Experimental Results

The experimental results are shown in Fig.5.4. The hybrid controller for orbital stabilization

was active for the first 20 s. At the end of this period, the system trajectories reached a close
neighborhood of the upright equilibrium rφ θ 9φ 9θs “ r0 ´ 3π{2 0 0s and the following linear
controller was invoked for stabilization:

τs “ 1.4φ´ 20.23pθ` 3π{2q` 1.14 9φ´ 1.98 9θ

The poles of the closed-loop system were located at ´37.0˘ 20.0 i and ´1.0˘ 1.2 i.

1.0

0.0

-1.0

1.57

-1.57

-4.71

0.3

-0.3

0

φ (rad)

θ (rad)

(a)

(c)

(e)

1.0

0.0

-1.0

1.5

0.0

-1.5

0

9φ (rad/s)

9θ (rad/s)

(b)

(d)

(f)

τ (Nm)

t (s)

-15

0

35

9V

t (s)

35

Figure 5.4: Rotary pendulum experimental results:
(a)-(d) are plots of joint angles and joint
velocities, (e) control torque, and (f) derivative of Lyapunov-like function. The brake pulses are
shown within plots (e) and (f), the peaks represent time intervals when the brakes were engaged.

97

0.5

0.0

-0.5

1.57

-1.57

-4.71

0

φ (rad)

θ (rad)

2.0

0.0

-2.0

7.0

0.0

-7.0

9φ (rad/s)

9θ (rad/s)

t (s)

20

0

t (s)

20

Figure 5.5: Rotary pendulum simulation results.

The pulses shown on the top of Figs.5.4 (e) and (f) correspond to the time intervals when the

brake was engaged (OFF) during orbital stabilization. The brake was disengaged (ON) when the
condition | 9φ |ď µ was satisfied; the value of µ was chosen to be small, equal to 0.1 rad/s. The time
intervals required for braking were very short (« 0.04 s, on average); this implies that the brakings
were impulsive in nature. The effect of impulsive braking can be seen in Fig.5.4 (b) where 9φ jumps

to almost zero value upon engagement of the brake on multiple occasions.

It can be seen from Fig.5.4 (c) that the amplitude of the pendulum gradually increases and

finally reaches a close neighborhood of the upright equilibrium configuration. The derivative of

the Lyapunov-like function is shown in Fig.5.4 (f). It can be seen that 9V never becomes positive;

this is because the brake is engaged every time when 9V is about to become positive7. Since 9V

is always negative, V decreases monotonically and orbital stabilization is achieved. A plot of the

motor torque is shown in Fig.5.4 (e). To minimize wear and tear of the brake, the commanded

motor torque was set to zero when the brake was engaged. A video of this experiment can be found

on the weblink: www.egr.msu.edu/~mukherji/RotaryPendulum.mp4

Simulation results for the same set of initial conditions and controller gains in (5.6.7) and

(5.6.8) are presented in Fig.5.5. A comparison of Figs.5.4 and 5.5 indicate that the joint velocities

7When | 9φ|ď µ « 0, the brake is not engaged since 9V « 0 - see (5.4.7).

98

in experiments are lower than those in simulations - this can be attributed to the presence of friction

and other dissipative forces. The amplitude of the active joint φ is larger in experiments than

simulations - this is due to the fact that the controller has to overcome the dissipative losses and

additional energy is added through larger amplitude of motion. As expected, the time needed for

stabilization is less in simulations than experiments.

Remark 9. For comparison, we considered the rotary pendulum example in [71]. Taking identical

initial conditions and physical parameters of the system therein, we simulated our controller with

the gains

Kp “ 0.20, Kd “ 0.12, Ke “ 50, Kc “ 0.70

The gains were tuned such that the magnitude of the motor torque did not exceed 0.3 Nm. The

system trajectories converged to the desired orbit in approx. 30 s. The controller in [71] took

approx. 100 sec and the magnitude of the maximum torque was 8 Nm. Our controller performed

well, both in terms of motor torque requirement and speed of convergence. This better performance,

however, comes at the cost of additional brake hardware.

5.7 Proofs and Additional Discussions

5.7.1 Proof of Lemma 1

The proof of Lemma 1 is provided here for the general case where the underactuated system has
m passive and n´ m active generalized coordinates, i.e. q1 P Rn´m, q2 P Rm and u P Rn´m. The
Euler-Lagrange equation has the same form as in (5.3.4) with M11 P Rpn´mqˆpn´mq, M22 P Rmˆm,
h1 P Rpn´mq, and h2 P Rm. The change in the total energy due to application of an impulsive input
is equal to the change in the kinetic energy:

1

9q`T

2

1

ΔE “
2„ 9q`T
“
` 9q`T

Mpqq 9q` ´
1 M11 9q`

1

2

Mpqq 9q´
1`
1 M11 9q´

9q´T
1 ´ 9q´T
1 M12 9q´

2

1 M12 9q`

2 ´ 9q´T

1

2„ 9q`T

2 M22 9q`

2 ´ 9q´T

2 M22 9q´

2

(A.1)

99

The impulse manifold, given in (5.3.9) for m “ 1, is
2 ´ M´1

9q`
2 “ 9q´

22 MT

12p 9q`

1 ´ 9q´
1 q

Substitution of 9q`

2 from (A.2) into (A.1) yields

(A.2)

1

2„ 9q`T

9q´T
2 M22 9q´

1

2

ΔE “

´

1 M11 9q`

1 M11 9q´

1 ´ 9q´T
2 ´ 9q´T

1 M12 9q´

22 MT

1

1`
2” 9q´
2 ´ M´1
1 M12” 9q´
2 ` 9q`T

2 ´ M´1

12p 9q`
22 MT

M22” 9q´
1 qıT
1 ´ 9q´
1 qı
1 ´ 9q´
12p 9q`

2 ´ M´1

22 MT

12p 9q`

1 qı
1 ´ 9q´

Expanding, canceling, and regrouping the terms on the right-hand side of the above equation yields

ΔE “

9q`T

1

1

2

”M11 ´ M12M´1

22 MT

12ı 9q`

1 ´

1

2

9q´T

1 ”M11 ´ M12M´1

22 MT

12ı 9q´

1

Similar to (5.3.7), Bpqq is defined for the general case as follows

Bpqq “”M11 ´ M12M´1

22 MT

12ı´1

(A.3)

(A.4)

From the properties of the mass matrix Mpqq, it can be shown that Bpqq is well-defined; also, it
is symmetric and positive-definite, i.e., Bpqq “ BTpqq ą 0. Substitution of (A.4) into (A.3) gives
(5.3.10).

(cid:3)

5.7.2 Proof of Lemma 2

Impulsive inputs result in no change in the generalized coordinates. Additionally, impulsive braking

results in 9q`

1 “ 0. Therefore, from the definition of V in (5.3.12), ΔV for impulsive braking can be

expressed as

ΔV “

“

“

1

1

2„KepE`´ Edesq2 ´ KepE´´ Edesq2 ´ 9q´T
2„KepE`` E´´ 2EdesqΔE ´ 9q´T
1
2„Kep2E`´ ΔE ´ 2EdesqΔE ´ 9q´T
1

1 Kd

1 Kd

9q´

9q´

1

1 Kd

9q´

1

100

where ΔE is defined in (5.3.10). Substitution of ΔE from (5.3.11) in the equation above yields

1

1

2

1

2

1

4

1

4

9q´

9q´T

1 Kd

“ ´

1 B´1 9q´

ΔV “ ´

9q´T
1 B´1 9q´

9q´T
1 B´1 9q´

1 qKetE`´ Edes `

1 u` 9q´T
1 uB´1 ` Kd 9q´

2„p 9q´T
1 „KetE`´ Edes `
1 „ 1

1
1* B´1 ` Kd ` KepE` ´ EdesqB´1 9q´
which is the same as in (5.3.13). Since B, defined in (A.4), is positive-definite, tKe 9q´T
1 uB´1
is positive-definite. Therefore, if”Kd ` KepE` ´ EdesqB´1pqqı is positive-definite, ΔV ď 0 and
ΔV “ 0 iff 9q´

4"Ke 9q´T

1 B´1 9q´

1 B´1 9q´

1 “ 0.

“ ´

9q´T

(cid:3)

1

1

5.7.3 Proof of Lemma 3

Let q˚
1 and u˚ denote the constant values of q1 and u, respectively. From the passivity property of
underactuated Euler-Lagrange systems [103,109], it follows that 9E “ uT 9q1 “ 0. Let E˚ denote the
constant value of E. Since 9q1 ” 0, using (5.3.4) and (5.3.5), we get the following

(A.5a)

(A.5b)

(A.6)

(A.7)

(A.8)

M12 :q2 `„BM12
Bq2  9q

2

2 ´

1

2
9q

Bq1 T
Bq1T
2 `„BF
2„BM22
“ u˚
Bq2  9q
2„BM22
2 ` BF
Bq2 “ 0

1

2

M22 :q2 `

Eliminating :q2 using (A.5a) and (A.5b), we get

1

˜„BM12
Bq2 ´

2 M22„BM22
Since E “ E˚ is constant and 9q1 “ 0, (5.2.2) gives

Bq1 T
2„BM22

Bq2 ¸ 9q

M12

´

2

2 ´

M12

M22

BF

Bq2 `„BF
Bq1T

“ u˚

2
2 “
9q

2pE˚ ´ Fq

M22

Substitution of 9q2

2 from (A.7) in (A.6) yields

E˚ P1pq˚

1, q2q` P2pq˚

1, q2q “ u˚

101

where P1pq˚
with respect to time results in

1, q2q, P2pq˚

1, q2q are defined in assumption 5. Differentiation of the above equation

dq2 rkP1pq˚
Using assumption 5, we claim q2 is constant.

d

1, q2q` P2pq˚

1, q2qs 9q2 “ 0

5.7.4 Nonoccurrence of Zeno Phenomenon

Impulsive braking results in 9q1pt`q “ 0 and the time derivative of p5.4.7) yields:

:Vpt`q “ ´ :qT

1 pt`qKc :q1pt`q ă 0

(cid:3)

(A.9)

since Kc is positive-definite and :q1pt`q ­“ 0. If ε denotes the time for the trajectories to re-enter
Z1, we can write

9Vpt`` εq “ 9Vpt`q`ż t``ε

t`

:Vpτqdτ

(A.10)

9Vpt`` εq “ 9Vpt`q “ 0, and :Vpt`q ­“ 0, :V must change sign in t P pt`,t`` εq. If ε is a

Since,
small number, then the continuity of :V in the interval rt`,t`` εs implies:

Using (A.9) and (A.11), we get

| :Vpt`q|“ Opεq
} :q1pt`q} ďd | :Vpt`q|

λminpKcq “ Opε1{2q

(A.11)

(A.12)

where λminpKcq is the smallest eigenvalue of Kc. Integrating :q1 with respect to time and using
(A.12), we get

} 9q1pt`` εq} ďż t``ε

t`

} :q1pτq}dτ “ Opε3{2q

(A.13)

Let us now assume that Zeno phenomenon occurs. Then, the time interval between consecutive

impulsive inputs, which is equal to ε, converges to zero and the trajectory of the system at time
t “ pt`` εq lies in the set Z1, where 9q1pt` ` εq ­“ 0. Since ε3{2 converges to zero faster than ε
converges to zero, (A.13) implies that 9q1pt`` εq also converges to zero faster than ε converges
to zero. This implies that the trajectory of the system does not lie in Z1 at t “ pt`` εq, and, by
contradiction, proves the nonoccurence of the Zeno phenomenon.

102

5.7.5 Well-posedness of Switching Times

5.7.5.1 Background

Let ψpt, 0, x0q denote the solution of the continuous-time dynamics in (5.3.15a) for 0 ď t ď t´
If and when the trajectory reaches
starting from the initial condition xp0q “ x0 at time t “ 0.
Z at t1, the states are instantaneously transferred to xpt`
resetting in (5.3.17). The trajectory xptq, t`
The resetting times are well-posed if the following conditions are satisfied [101, pg.13]:

1 q` fdrxpt´
2 , is then given by ψpt,t`

1 qs, according to the
1 , xpt`
1 qq, and so on.

1 ď t ď t´

1 q “ xpt´

1

1. If xptq P ¯Z zZ , then there exists ε ą 0 such that, for all 0 ă δ ă ε,

ψpt ` δ,t, xptqq R Z

2. If xpt´

k q P BZ X Z , then there exists ε ą 0 such that, for all 0 ď δ ă ε,

ψpt`

k ` δ,t`

k , xpt´

k q` fdrxpt´

k qsq R Z

Condition (1) ensures that when the system trajectory reaches the closure of the set Z , the trajectory

must be directed away from Z . Condition (2) ensures that when the trajectory intersects Z , it

instantaneously exists Z .

5.7.5.2 Proof

Consider the set Z1 defined in (5.4.3b). If xptq P ¯Z1zZ1, then 9q1ptq “ 0, which implies 9Vptq “ 0.
Differentiating (5.4.7) with respect to time and substituting 9q1ptq“ 0, we get :Vptq“´ :qT
1 ptqKc :q1ptqă
9Vpt`εq ă 0. Using (5.4.8), rA` BusT Kc 9q1pt`
0. Since 9Vptq “ 0 and :Vptq ă 0, it follows lim
εq ą 0, the trajectory of the system is directed away from Z1 after it reaches the closure of Z1.
This implies condition (1) is satisfied for Z1. It can be shown that

εÑ0`

BZ1 X Z1 “ trApq, 9qq` BpqqusT Kc 9q1 “ 0, 9q1 ­“ 0u

103

k q P BZ1 X Z1, impulsive braking is applied. This results in 9q1pt`

k q “ 0, which implies
k q R Z1. Thus the trajectory instantaneously exists Z1 after intersecting Z1 and condition (2)

If xpt´
xpt`
is satisfied for Z1.

Now consider the set Z2 fi tq2 “ ´3π{2, 9q2 ă 0u, which was defined after (5.3.17). If xptq P
¯Z2zZ2 i.e., xptq P tq2 “ ´3π{2, 9q2 “ 0u. Then, if lim
εÑ0` 9q2pt ` εq ě 0, then it follows from
the definition of Z2 that xpt ` εq R Z2. Alternatively, if lim
εÑ0` 9q2pt ` εq ă 0, then coordinate
resetting occurs and q2pt ` εq “ π{2 Ñ xpt ` εq R Z2. Thus condition (1) is satisfied for Z2.
Condition (2) is trivially satisfied for Z2 since BZ2 X Z2 is an empty set.

Similar to the set Z2, it can be shown that conditions (1) and (2) are satisfied for the set

Z3 fi tq2 “ π{2, 9q2 ą 0u. Since conditions (1) and (2) are satisfied individually for Z1, Z2, Z3,
they are satisfied for Z “ Z1 Y Z2 Y Z3.

(cid:3)

5.7.6 Quasi-Continuous Dependence Property

5.7.6.1 Background

Let spt, 0, x0q denote the solution of the hybrid dynamical system in (5.3.15a) and (5.3.15b) for t ě 0
starting from the initial condition xp0q “ x0 at time t “ 0. Let Tx0
fi r0,8qztt1px0q,t2px0q,¨¨¨u
denote all times except the discrete instants at which impulsive inputs are applied or coordinate

resetting occurs. We now state the following assumption [101, Assumption 2.1, pg.27]:

Assumption 3. Consider the hybrid dynamical system in (5.3.15). For every x0 P D and for
every ε ą 0 and t P Tx0, there exists δpε, x0,tq ą 0 such that if }x0 ´ y} ă δpε, x0,tq, y P D , then
}spt, 0, x0q´ spt, 0, yq} ă ε.

The above assumption is a generalization of the standard continuous dependence property for

dynamical systems with continuous flows to dynamical systems with left-continuous flows. We now

state a proposition that provides sufficient conditions for satisfying Assumption 3 [101, Proposition

2.1, pg.32]:

104

Proposition 1. Consider the hybrid dynamical system in (5.3.15) and suppose that conditions (1)

and (2) in section 5.7.5.1 are satisfied; then the system in (5.3.15) satisfies Assumption 3 if for all

x0 P D , 0 ď t1px0q ă 8, t1px0q is continuous, and limkÑ8 tkpx0q Ñ 8.

5.7.6.2 Proof

Our domain D was defined in section 5.3.3 as the open set where q2 P pa, bq, a ă ´3π{2, b ą π{2.
If q2p0q P pa,´3π{2s Y rπ{2, bq, we can reset q2 to satisfy (5.4.2) since q2 corresponds to a
revolute joint. From (5.3.15) we know that x0 R Z . This implies that t1px0q ą 0. Furthermore,
t1px0q is a continuous function of x0 since the system is described by continuous-time dynamics
(ordinary differential equations) over the interval t P r0,t1q. Also, as shown in section 5.7.4, Zeno
phenomenon does not occur. Thus limkÑ8 tkpx0q Ñ 8.

(cid:3)

5.7.7 Verification of Assumption 2 for Tiptoebot and Rotary Pendulum

Tiptoebot: It can be seen from (5.5.2) and (5.5.4) that all elements of the mass matrix are bounded

and the potential energy is lower bounded. It can also be seen from (5.5.2) that M12 and M22 are

only functions of q1. Therefore P1 in (A.8), defined earlier in assumption 5, is not a function of

q2. Since F in (5.5.4) is a function of both q1 and q2, it can be easily shown that pdP2{dq2q ­“ 0.
Thus assumption 5 holds.

Rotary Pendulum: It can be seen from (5.6.2) and (5.6.4) that all elements of the mass matrix are

bounded and the potential energy is lower bounded. It can also be seen from (5.6.2) and (5.6.4)

that the mass matrix and potential energy are independent of q1 “ φ. Therefore

d
dθ rkP1pθq` P2pθqs
d

cosθp2E˚ ´ 3γ4 sinθq ­“ 0

dθ„γ3

γ2

“

Thus assumption 5 holds.

105

CHAPTER 6

ORBITAL STABILIZATION USING VIRTUAL HOLONOMIC CONSTRAINTS AND

IMPULSE CONTROLLED POINCARÉ MAPS

6.1 Introduction

The problem of orbital stabilization of underactuated mechanical systems with one passive

degree-of-freedom (DOF) is revisited in this chapter. The system dynamics is presented in section

6.2 and the results in [62] are utilized to design a continuous controller that can enforce the

VHC such that the resulting zero dynamics is Euler-Lagrange. A periodic orbit is selected on the

constraint manifold and a method for orbital stabilization is presented in section 6.3. To stabilize

the orbit, a Poincaré section is defined at a point on the orbit and the return map is linearized about

the fixed point; this results in a 2n´1 dimensional discrete linear time-invariant (LTI) system. To
control this system and stabilize the orbit, impulsive inputs are applied when the system trajectory

crosses the Poincaré section. The controllability of the orbit can be verified by simply checking

the controllability of the linear system. This is simpler than the approach in [3, 65, 66] where

controllability is verified numerically along the orbit for most systems. Since the system is LTI,

the control design involves constant gains that can be computed off-line. Compared to the methods

proposed earlier [3, 63, 65, 66], where periodic Ricatti equations have to be solved, our method

has lower computational cost and complexity. Since impulsive inputs are used to control the

Poincaré map, the dynamics of the closed-loop system can be described by the Impulse Controlled

Poincaré Map (ICPM). The simplicity and generality of the ICPM approach to orbital stabilization

is demonstrated using the examples of the 2-DOF cart-pendulum in section 6.4 and the 3-DOF

tiptoebot in section 6.5.

106

6.2 Problem Formulation

6.2.1 System Dynamics

Consider an n DOF underactuated system with one passive DOF, where the passive DOF is a

revolute joint. Let q, q fi“qT
1 q2‰T , denote the generalized coordinates, where q1 P Rn´1 and
q2 P S1, S1 “ R modulo 2π, are the coordinates of the active and passive DOFs. The configuration
space of the system is denoted by Qn, Qn P Rn´1ˆ S1. The Lagrangian of the system can be written
as

Lpq, 9qq “

1

2

9qT Mpqq 9q` Fpqq

In the equation above, Mpqq P Rnˆn denotes the symmetric, positive-definite mass matrix, parti-
tioned as

Mpqq “»—–

M11pqq M12pqq
MT
12pqq M22pqq

fiffifl

where M11 P Rpn´1qˆpn´1q, M22 P R and Fpqq is the potential energy of the system. The Euler-
Lagrange equation of motion can be written as follows

M11pqq :q1 ` M12pqq :q2 ` h1pq, 9qq “ u
MT
12pqq :q1 ` M22pqq :q2 ` h2pq, 9qq “ 0

(6.2.1a)

(6.2.1b)

where u P Rn´1 is the control input, and rhT
forces. In compact form, (6.2.1a) and (6.2.1b) can be rewritten as

1 , h2sT is the vector of Coriolis, centrifugal and gravity

:q1 “ Apq, 9qq` Bpqqu
:q2 “ Cpq, 9qq` Dpqqu

where,

Bpqq “”M11 ´p1{M22qM12 MT
Dpqq “ ´p1{M22qMT

12ı´1

, Apq, 9qq “ p1{M22qBpqqrM12 h2 ´ h1M22s

12 Bpqq, Cpq, 9qq “ ´p1{M22q”MT

12 Apq, 9qq` h2ı

Similar to [62], we make the following assumption:

(6.2.2a)

(6.2.2b)

(6.2.3)

107

Assumption 4. For some ¯q fi r ¯qT
are even with respect to ¯q, i.e.,

1 , ¯q2sT P Qn, the mass matrix Mpqq and the potential energy Fpqq

Mp ¯q` qq “ Mp ¯q´ qq, Fp ¯q` qq “ Fp ¯q´ qq

6.2.2

Imposing Virtual Holonomic Constraints (VHC)

A holonomic constraint enforced by feedback is referred to as VHC. The current and the next

subsection summarizes relevant results from [62]. For a wide class of mechanical systems, a

comprehensive discussion on VHC can be found in [62], [85].

A VHC for (6.2.1) is described by the relation ρpqq “ 0 where, ρ : Qn Ñ Rn´1 is smooth and
rankrJqpρqs “ n´1 for all q P ρ´1p0q. Here, Jqpρq is the Jacobian of ρ with respect to q. The
VHC is stabilizable if there exists a smooth feedback ucpq, 9qq that asymptotically stabilizes the set

C “ tpq, 9qq : ρpqq “ 0, Jqpρq 9q “ 0u

(6.2.4)

The set C , which is referred to as the constraint manifold, is controlled invariant [62]. For the
system described by (6.2.1), the set C is the tangent bundle of ρ´1p0q, which is a closed embedded
submanifold of Qn of codimension pn´ 1q [110, Def. 4.2].

An important goal of this work is to generate repetitive motion, which can be described by closed
orbits. Consequently, ρ´1p0q must be a smooth and closed curve without any self-intersection.
The VHC can be described as

ρpqq “ q1 ´ Φpq2q “ 0

(6.2.5)

where Φ : S1 Ñ Rn´1 is a smooth vector-valued function. The constraint manifold C in (6.2.4) can
be expressed as:

C “"pq, 9qq : q1 “ Φpq2q, 9q1 “„ BΦ

Bq2 9q2*

(6.2.6)

It should be noted that since q2 P S1, Φpq2 ` 2πq “ Φpq2q and ρ´1p0q is closed. Following the
notion of odd VHC [62], we make another assumption.

108

Assumption 5. For ¯q which satisfies Assumption 4, Φpq2q is odd with respect to ¯q2, i.e.,

To stabilize C , we investigate the dynamics of ρpqq; differentiating ρpqq twice with respect to

Φp ¯q2 ` q2q “ ´Φp ¯q2 ´ q2q

time, we get

Substitution of :q1 and :q2 from (6.2.2a) and (6.2.2b) in (6.2.7) yields

2
2

2ff 9q
Bq2 :q2 ´«B2Φ
:ρ “ :q1 ´„ BΦ
Bq2
2ff 9q
:ρ “ A´«B2Φ
Bq2C`„B´„ BΦ
2 ´„ BΦ
Bq2
2ff 9q
Bq2 D´1«´A`«B2Φ
2 `„ BΦ
Bq2

Bq2 D u
Bq2C´ kpρ´ kd

2

uc “„B´„ BΦ

2

(6.2.7)

(6.2.8)

(6.2.9)

(6.2.10)

9ρff

The following choice of linearizing control

where kp and kd are positive definite matrices, results in

:ρ` kd

9ρ` kp ρ “ 0

This implies that limtÑ8 ρptq Ñ 0 exponentially and uc in (6.2.9) stabilizes the VHC in (6.2.5). If
the initial conditions are chosen such that ρp0q “ 9ρp0q “ 0, uc in (6.2.9) enforces the VHC and the
constraint manifold C is controlled invariant.

Remark 10. For uc in (6.2.9) to be well-defined, the matrix rB´pBΦ{Bq2qDs must be invertible.
It can be shown that rB´pBΦ{Bq2qDs is invertible iff MT
12pBΦ{Bq2q` M22 ‰ 0. This is also a
necessary and sufficient condition for the VHC in (6.2.5) to be regular; since the VHC in (6.2.5) is

in parametric form, this is also a necessary and sufficient condition for C to be stabilizable - see

proposition 3.2 of [62].

6.2.3 Zero Dynamics and Periodic Orbits

On the constraint manifold C , the dynamics of the system satisfies ρpqq ” 0; this implies

q1 “ Φpq2q, 9q1 “„ BΦ

2ff 9q
Bq2 9q2, :q1 “«B2Φ
Bq2

2

2 `„ BΦ
Bq2 :q2

(6.2.11)

109

Substitution of q1, 9q1 and :q1 from (6.2.11) in (6.2.1b) provides the zero dynamics, which can be

expressed in the following form

It was shown in [3, 62, 83, 84] that the equation above has an integral of motion of the form

:q2 “ α1pq2q` α2pq2q 9q

2
2

(6.2.12)

Epq2, 9q2q “ p1{2qMpq2q 9q

Mpq2q “expˆ´2ż q2

0

α2pτqdτ˙ , Ppq2q “ ´ż q2

0

2

2 ` Ppq2q

α1pτqMpτq dτ

(6.2.13)

where Mpq2q represents the mass and Ppq2q represents the potential energy of the reduced
system in (6.2.12). Since both Assumption 4 and 5 are satisfied, the zero dynamics represents an
Euler-Lagrange system with the Lagrangian1 equal to p1{2qMpq2q 9q2

2 ´ Ppq2q.

Since Assumptions 1 and 2 hold, the zero dynamics in (6.2.12) is similar to the dynamics of

a simple pendulum [62] and its qualitative properties can be described by the potential energy

Ppq2q. Let Pmin and Pmax denote the minimum and maximum values of P. If an energy level
set is denoted by Epq2, 9q2q “ c, then c P pPmin, Pmaxq corresponds to a periodic orbit where the
sign of 9q2 changes periodically and c ą Pmax corresponds to an orbit where the sign of 9q2 does
not change [62].

6.2.4 Problem Statement

Since the zero dynamics in (6.2.12) has an Euler-Lagrange structure, there cannot exist any non-

trivial isolated periodic orbit - this follows from the Poincaré-Lyapunov-Liouville-Arnol’d theorem

[62, 85, 111]. A direct implication of this theorem is that the reduced dynamics possesses a dense

set of closed orbits that are stable, but not asymptotically stable. For a desired repetitive motion, the

corresponding orbit must be stabilized. Consider the desired closed orbit Od, defined as follows:

Od “ tq, 9q P C : Epq2, 9q2q “ cdu,

cd ą Pmin

(6.2.14)

1Assumptions 4 and 5 provide necessary and sufficient conditions for the reduced system to be

Euler-Lagrange - the proof of this result can be found in [62, 85].

110

Let x, x fi rqT , 9qTsT , denote the states of the system in (6.2.1). We define an ε-neighborhood of
Od by

Uε “ tx P Qn ˆ Rn : distpx, Odq ă εu
distpx, Odq fi inf

d }x´ y}

yPO

We now define stability of the orbit Od from [112].

Definition 1. The orbit Od in (6.2.14) is

• stable, if for every ε ą 0, there is a δ ą 0 such that xp0q P Uδ ùñ xptq P Uε, @t ě 0.

• asymptotically stable if it is stable and δ can be chosen such that limtÑ8 distpxptq, Odq “ 0.

The control uc in (6.2.9) asymptotically stabilizes C but does not asymptotically stabilize Od.
If q, 9q P Od, uc enforces the VHC and trajectories stay on Od; however, a perturbation of the states
will cause the trajectories to converge to a different orbit on C . The objective of this work is to

utilize the control uc in (6.2.9), which was designed to stabilize C , together with impulsive inputs

to stabilize Od directly, that lies on C .

6.3 Main Result: Stabilization of Od

6.3.1 Poincaré Map

The system in (6.2.1) with u “ uc defined in (6.2.9), has the state-space representation

9x “ fpxq

(6.3.1)

The stability characteristics of periodic orbits can be studied using Poincaré maps [113]. To this

end, we define the Poincaré section Σ of Od as follows2:

Σ “ tx P Qn ˆ Rn : q2 “ q˚
2Since the periodic orbit Od can pass through q2 “ q˚

2, 9q2 ě 0u
2 with both positive and negative velocity
9q2, the condition 9q2 ě 0 implicitly assumes that Σ intersects the periodic orbit Od at a single point.
In the definition of Σ in (6.3.2), 9q2 ě 0 can be replaced with 9q2 ď 0 without any loss of generality.

(6.3.2)

111

where q˚

2 is a constant. Let z, z fi rqT

1 , 9qTsT P Rp2n´1q, denote the states of the system on Σ. The
Poincaré map P : Σ Ñ Σ is obtained by following trajectories of z from one intersection with Σ to
the next. Let tk, k “ 1, 2,¨¨¨ denote the time of the k-th intersection and zpkq “ zptkq. Then, zpk` 1q
can be described with the help of the map P

zpk` 1q “ Przpkqs

(6.3.3)

The point of intersection of Σ and Od is the fixed point of P denoted by z˚; it satisfies the following

relation

z˚ “ Ppz˚q

(6.3.4)

The stability characteristics of the orbit Od can be studied by investigating the stability properties
of z˚, which is an equilibrium point of the discrete-time system in (6.3.3); this can be done by
linearizing the map P about z˚. For zpkq “ z˚ ` ν, where }ν} is a small number, we can write

zpk` 1q “ Ppz˚ ` νq “ Ppz˚q`r∇zPpzqsz“z˚ rzpkq´ z˚s` Op}ν}2q

(6.3.5)

Using Ppz˚q “ z˚ from (6.3.4) and neglecting higher-order terms in }ν}, the above equation can be
written as

epk` 1q “ A epkq

epkq fi zpkq´ z˚, A fi r∇zPpzqsz“z˚

(6.3.6)

The stability properties of z˚ is governed by the eigenvalues of A , which are referred to as

the Floquet multipliers of Od [113].

If the Floquet multipliers lie inside the unit circle, Od is

exponentially stable - see Theorem 7.3 of [112]. From our discussion in section 6.2.4 we know

that the desired orbit Od is not asymptotically stable, i.e., not all eigenvalues of A lie inside the
unit circle. To asymptotically stabilize the orbit, i.e., to asymptotically stabilize z˚, we design an

impulsive controller in the next subsection.

112

6.3.2

Impulse Controlled Poincaré Map (ICPM)

To asymptotically stabilize the desired orbit Od, our controller in (6.2.9) is modified as follows

u “ uc ` uI

(6.3.7)

where uI is an impulsive input which is applied only when xptq P Σ. The dynamics of the system,
with uI as the new input, can be written as

M11pqq :q1 ` M12pqq :q2 ` ¯h1pq, 9qq “ uI
MT
12pqq :q1 ` M22pqq :q2 ` h2pq, 9qq “ 0

(6.3.8a)

(6.3.8b)

where ¯h1 fi ph1 ´ ucq. Impulsive inputs cause discontinuous changes in the generalized velocities
while there is no change in the generalized coordinates. On the Poincaré section Σ, the jump in

velocities can be computed by integrating (6.3.8) as follows [86]:

uI dt

(6.3.9)

»—–

M11 M12

MT

12 M22

fiffifl

»—–

Δ 9q1

Δ 9q2

fiffifl “»—–

I

0

fiffifl , I fiż Δt

0

In the above equation, Δt is the infinitesimal interval of time for which uI is active, I P Rn´1 is
the impulse of the impulsive input, and Δ 9q1 and Δ 9q2 are defined as

Δ 9q1 fi p 9q`

1 ´ 9q´
1 q,

Δ 9q2 fi p 9q`

2 ´ 9q´
2 q

(6.3.10)

where 9q´ and 9q` are the velocities immediately before and after application of uI . Since the

system is underactuated, the jump in the passive velocity 9q2 is dependent on the jumps in the active

velocity 9q1; this relationship is described by the pn´1q dimensional impulse manifold [6,23], which
can be obtained from (6.3.9):

IM “ t 9q`

1 , 9q`

2 | Δ 9q2 “ ´p1{M22qMT

12

Δ 9q1u

(6.3.11)

Since impulsive inputs can cause the system states to move on Σ, we exploit this property to design

a feedback law that asymptotically stabilizes z˚, i.e., asymptotically stabilizes Od. The control

113

input applied at tk is denoted by I pkq3. The dynamics of the impulse controlled system in (6.3.3)
can be described by the map

zpk` 1q “ Przpkq, I pkqs

(6.3.12)

where I pkq “ 0 if zpkq “ z˚ 4. By linearizing the above map about the fixed point z “ z˚ and I “ 0,
we get

epk` 1q “ A epkq` B I pkq

A fi”∇zPpz, I qız“z˚, I “0

, B fi”∇I Ppz, I qız“z˚, I “0

(6.3.13)

where A P Rp2n´1qˆp2n´1q and B P Rp2n´1qˆpn´1q can be obtained numerically. Since A is not
Hurwitz (see discussion in the last sub-section), we make the following proposition to exponentially

stabilize Od:

Proposition 2. If the pair tA , Bu is stabilizable, the orbit Od can be exponentially stabilized

locally using the discrete impulsive feedback

I pkq “ K epkq

(6.3.14)

where the matrix K is chosen such that pA ` BK q is Hurwitz.

Remark 11. The matrices A and B depend on the continuous controller in (6.2.9), and hence on

the choice of controller gains kp and kd. Thus, the stabilizability of the pair tA , Bu depends on
the choice of kp and kd.

The above approach to stabilization, which we refer to as the impulse controlled Poincaré map

(ICPM) approach, is explained with the help of the schematic in Fig.6.1. The desired orbit Od is
shown in red and it intersects the Poincaré section Σ at the fixed point z˚. A trajectory starting from

an arbitrary initial condition, shown by the point 1 , intersects Σ at 2 . The impulsive input in

3As long as Δt is sufficiently small, the effect of the impulsive input uI depends solely on the

value of I - see (6.3.9). Thus I can be viewed as the control input.

4If the system trajectory is passing through z˚, the trajectory is evolving on Od, which lies on C .
Since I “ 0, the invariance of Od is guaranteed solely by the continuous controller uc in (6.2.9).

114

Σ

3

4

z˚

2

IM

1

O

d

Figure 6.1: Schematic of ICPM approach to orbital stabilization.

(6.3.14) moves the configuration of the system from 2 to 3 along the impulse manifold IM. The
point 3 lies on Σ since impulsive inputs cause no change in the position coordinates and q2 “ q˚
is maintained during the transition from 2 to 3 along IM. Furthermore, the condition 9q`
2 ě 0
is satisfied since the impulsive controller in (6.3.14) is designed using linearization and 3 lies in

2

some small neighborhood of z˚ on Σ. Hereafter, the altered system trajectory evolves under the

continuous control uc and 4 denotes its next intersection with Σ. A series of ICPMs, similar to
the map 2 Ñ 4 exponentially converge the intersection point of the trajectory on Σ to z˚.
Remark 12. By stabilizing the constraint manifold C , the continuous controller keeps the system

trajectory close to C . By stabilizing the fixed point, the impulsive controller works in tandem with

the continuous controller to converge the system trajectory to Od, which lies on C . Although the

impulsive control inputs perturb the system trajectory intermittently, the trajectory converges to

Od on C due to the vanishing nature of the perturbations. While global exponential stability of

C guarantees that C remains exponentially stable despite the perturbations, the magnitudes of

the perturbations exponentially converge to zero as the intersection point of the trajectory with the

Poincaré section converges to the fixed point.

115

Remark 13. In the ICPM approach, impulsive inputs create discontinuous jumps in the states of the

system. Although the subsequent continuous-time dynamics remains unchanged, the discontinuous

jumps in the states result in a change in the Poincaré map. This is different from the well-known

OGY method of chaos control [114] where the continuous-time dynamics is altered by discretely

changing system parameters on a Poincaré section. The OGY method has been utilized in control

of dynamical systems described by Poincaré maps; examples include bipeds [115] and hopping

robots [8].

6.3.3

Implementation of Control Design

6.3.3.1 Numerical Computation of A and B matrices

Let δi, i “ 1, 2,¨¨¨ , 2n´1, denote the i-th column of ε1Ip2n´1q, where ε1 is a small number and
Ip2n´1q is the identity matrix of size p2n´1q. If Ai denotes the i-th column of A , then Ai can be
numerically computed as follows:

Ai “

1

ε1 rPpz˚ ` δiq´ z˚s

(6.3.15)

Let Q P Rnˆpn´1q and S P Rpn´1q be defined as follows:

Q fi»—–

Ipn´1q

01ˆpn´1q

fiffifl ,

S fi»—–

0pn´1qˆ1
Mpqq´1ηi

fiffifl

where 0iˆ j is a matrix of zeros of dimension iˆ j, and ηi, i “ 1, 2,¨¨¨ , n´1, denote the i-th column
of ε2Q, where ε2 is a small number. If Bi denotes the i-th column of B, then Bi can be numerically

computed as follows:

Bi “

1

ε2 rPpz` Sqsz“z˚ ´ z˚(

The above expression has been obtained using (6.3.9).

(6.3.16)

Remark 14. Numerical computation of A and B matrices is sensitive to the choice of ε1 and

ε2. While ε1 and ε2 should be small, excessively small values can lead to numerical errors. For

116

example, the eigenvalues of A may be found to lie inside the unit circle, which cannot be the case

since orbits on C are not asymptotically stable.

6.3.3.2

Impulsive Input using High-Gain Feedback

Impulsive inputs are Dirac-delta functions and cannot be realized in real physical systems. Using

singular perturbation theory [116], it was shown that continuous-time approximation of impulsive

inputs can be carried out using high-gain feedback [23]; this has allowed implementation of

impulsive control in physical systems using standard hardware [6, 25]. To obtain the expression for

the high-gain feedback, we substitute (6.3.14) in (6.3.9) to get

Δ 9q1pkq “ B I pkq “ BK epkq

(6.3.17)

where B is defined in (6.2.3) and is evaluated at tk; Δ 9q1pkq is the jump in the active velocities
generated by the input I pkq. From (6.3.10) and (6.3.17), the desired active joint velocities at tk is

9qdes
1 pkq “ 9q1pkq` BK epkq

(6.3.18)

where 9q1pkq “ 9q1ptkq. To reach the desired velocities in a very short period of time, we use the
high-gain feedback [6]

µ

uhg “ B´1„ 1
which remains active for as along as } 9qdes
9qdes
1
h1 with ¯h1. Furthermore, Λ fi diagrλ1 λ2
numbers, and µ ą 0 is a small number.
6.4 Illustrative Example: Cart-Pendulum

1 pkq´ 9q1} ě ε3, where ε3 is a small number. In (6.3.19),
is obtained from (6.3.18) and ¯A is obtained from the expression for A in (6.2.3) by replacing
¨¨¨ λn´1s, where λi, i “ 1, 2,¨¨¨ , n´1 are positive

1 pkq´ 9q1¯´ ¯A
Λ´ 9qdes

(6.3.19)

6.4.1 System Dynamics and VHC

Consider the frictionless cart-pendulum system in Fig.6.2. The masses of the cart and pendulum

are denoted by mc and mp, ℓ denotes the length of the pendulum, and g is the acceleration due to

117

mp

θ

ℓ

g

u

x

mc

Figure 6.2: Inverted pendulum on a cart.

gravity. The control input u is the horizontal force applied on the cart. The cart position is denoted

by x, x P R, and the angular displacement of the pendulum, measured clock-clockwise with respect
to the vertical, is denoted by θ, θ P S1. We consider physical parameters of the system to be the
same as those in [3]: mp “ mc “ ℓ “ 1. With the following definition

q “ rq1 q2sT “ rx θsT

and the potential energy of the system, given by

the equations of motion can be obtained as

F “ cosθ

»—–

sinθ 9θ2

g sinθ

fiffifl “»—–

u

0

which is of the form in (6.2.1). The VHC in (6.2.5) is chosen as

2

cosθ

cosθ

1

fiffifl

»—–

:x

:θ

fiffifl´»—–

ρ “ x` 1.5 sinθ “ 0

(6.4.1)

(6.4.2)

(6.4.3)

(6.4.4)

fiffifl

which is identical to that considered in [3]. It can be verified that the mass matrix in (6.4.3) and

the choice of VHC in (6.4.4) satisfy Assumptions 4 and 5 for ¯q “ p0, 0q. For the VHC in (6.4.4) to
be stabilizable, Remark 10 provides the following condition that needs to be satisfied:

1´ 1.5 cos

2 θ ‰ 0 ñ θ ­“ ˘0.61 rad

(6.4.5)

118

Clearly, the VHC in (6.4.4) is not regular and the corresponding constraint manifold is therefore

not stabilizable. Since our control design requires the VHC to be regular, we cannot use the VHC
in (6.4.4) with θ P S1. However, by restricting the domain of θ to be p´0.61, 0.61q, it is possible to
compare the performance of our control design with that presented in [3]. Through trial and error,

the controller gains are chosen such that θ lies in the interval p´0.61, 0.61q.

6.4.2 Stabilization of VHC and Od

The ICPM approach relies on stabilization of both the constraint manifold C , and the orbit Od

on C . This is a distinctive difference between our approach and the approach in [3] where Od is

stabilized without stabilizing C , i.e., without enforcing the VHC. With the objective of enforcing

the VHC, the gains kp and kd in (6.2.9) are chosen as follows:

We choose the desired orbit Od to pass through the point:

kp “ 2,

kd “ 1

px, θ, 9x,

9θq “ p0.0, 0.0, ´0.675, 0.450q

(6.4.6)

(6.4.7)

which is approximately the desired orbit in [3] - see Fig.2 therein. For exponential stabilization of

Od, we define the Poinacaré section

Σ “ tx P Q

2 ˆ R

2

: θ “ 0,

9θ ě 0u

(6.4.8)

The states of the system on Σ are

Since z˚ lies on Od, using (6.4.7) and (6.4.8) we get

z “ r x

9x

9θsT

z˚ “ r0.0 ´0.675 0.450sT

119

Using the values of ε1 “ 0.02 and ε2 “ 0.01, the matrices A and B in (6.3.15) and (6.3.16) are
obtained as

0.115

0.435

0.600

´0.510 ´0.640 ´2.465
´0.145
1.325

0.215

, B “

fiffiffiffiffifl

»————–

´0.06
1.80

´1.09

fiffiffiffiffifl

A “

»————–

It can be verified that all eigenvalues of A do not lie inside the unit circle but the pair tA , Bu
is controllable and satisfy Proposition 2. Using LQR design, the gain matrix K in (6.3.14) was

obtained as

K “ r 0.163 0.288 1.198 s

(6.4.9)

The eigenvalues of pA ` BK q are located at 0.13 and ´0.06˘ 0.48i; thus, the impulsive feedback
in (6.3.14) exponentially stabilizes the desired orbit Od.

6.4.3 Simulation Results

The initial configuration of the system is taken from [3]:

rx θ 9x

9θs “ r0.1 0.4 ´0.1 ´0.2s

For the controller gains in (6.4.6) and (6.4.9), simulation results for the ICPM are shown in Fig.6.3;

ρ is plotted with time in Fig.6.3 (a) and the phase portrait of the pendulum is shown in Fig.6.3

(b). To show the convergence of system trajectories to Od, }epkq}2 is plotted with respect to k in
Fig.6.3 (c). It can be seen from Fig.6.3 (a) that the system trajectories converge to the constraint

manifold. To stabilize Od on C , the impulsive controller in (6.3.14) is implemented using the
high-gain feedback in (6.3.19) with Λ “ 1 and µ “ 0.005. It can be seen from the phase portrait in
Fig.6.3 (b) that the pendulum trajectory converges exponentially to Od, shown in red. The effect of

discrete impulsive feedback can be seen in Fig.6.3 (b) where 9θ jumps when trajectories cross the

Poincaré section Σ defined in (6.4.8). The system trajectories reach a close neighborhood of Od in

approximately 10 sec; this is comparable to the results in [3]. We now consider the following initial

120

(b)

9θ vs θ

(a)

ρ

2.0

0.0

-2.0

time (s)

10

-0.5

0.0

0.5

}epkq}2

(c)

k

5

0.6

0.0

-0.6

0

3.5

0.0

0

Figure 6.3: Orbital stabilization for the cart-pendulum system for the initial conditions in [3]: (a)
plot of ρ with respect to time, (b) phase portrait of the pendulum, (c) plot of the norm of the error
of the discrete-time system.

condition that lies far away from Od:

rx θ 9x

9θs “ r0.0 0.0 0.0 0.0s

(6.4.10)

We used the same controller gains as that used in the previous simulation. It can be seen from the

results in Fig.6.4 that Od is asymptotically stabilized. For the same initial conditions in (6.4.10),

the control design in [3] fails to converge the pendulum trajectory to Od. Several other initial

conditions yielded similar results. While this may be indicative of a larger region of attraction

of Od with the ICPM approach than with the approach in [3], the region of attraction has to be

estimated for both designs and compared before any conclusion can be drawn. To demonstrate the

generality of the ICPM approach, we consider the three DOF tiptoebot, which is presented next.

121

(a)

ρ

time (s)

}epkq}2

0.09

0.00

0

0.1

0.0

0

0.5

0.0

(b)

9θ vs θ

-0.5

-0.15

15

(c)

0.0

0.15

k

5

Figure 6.4: Orbital stabilization for the cart-pendulum system for the initial conditions in (6.4.10):
(a) plot of ρ with respect to time, (b) phase portrait of the pendulum; Od is shown in red, (c) plot
of the norm of the error of the discrete-time system.

6.5 Illustrative Example - The Tiptoebot

6.5.1 System Description

Consider the three DOF tiptoebot [6]. Using the following definition for the joint angles and control

inputs

1 “ r θ2 θ3 sT ,
qT

q2 “ θ1,

u “ rτ2 τ3sT

(6.5.1)

Table 6.1: Tiptoebot lumped parameters in SI units

α1

0.386 α4

0.065 β1

4.307

α2

0.217 α5

0.054 β2

1.102

α3

0.247 α6

0.104 β3

1.764

122

where θ1 P S1, θ2,θ3 P R, the dynamics of the tiptoebot can be expressed in the form given in
(2.2.1), where the components of the mass matrix and the potential energy are:

α2`α3`2α5 cosθ3 α3`α5 cosθ3

α3`α5 cosθ3

α3

fiffifl

M11 “»—–
M12 “»—–

α2`α3`α4 cosθ2`2α5 cosθ3`α6 cospθ2`θ3q

(6.5.2)

α3`α5 cosθ3`α6 cospθ2`θ3q

fiffifl

M22 “ α1`α2`α3 ` 2rα4 cosθ2 ` α5 cosθ3`α6 cospθ2`θ3qs
F “ β1 cosθ1 ` β2 cospθ1 ` θ2q` β3 cospθ1 ` θ2 ` θ3q

where αi, i “ 1, 2,¨¨¨ , 6, and βi, i “ 1, 2, 3 are lumped physical parameters; their values are given
in Table 6.1. It can be verified that Assumption 4 is satisfied for ¯q “ p0 0 0qT .

6.5.2

Imposing VHC and Selection of Od

The VHC in (6.2.5) is chosen as

ρ “»—–

ρ1

ρ2

fiffifl “»—–

θ2 ` 2.0θ1
θ3 ´ 0.1θ1

fiffifl “»—–

0

0

fiffifl

kp “»—–

1.0 0.0

0.0 1.0

fiffifl ,

kd “»—–

0.1 0.0

0.0 0.1

fiffifl

(6.5.3)

(6.5.4)

which satisfies Assumption 5 for ¯q “ p0 0 0qT . Also, it is stabilizable as it satisfies the condition in
Remark 10. With the objective of enforcing the VHC, the gain matrices in (6.2.9) were chosen as

The phase portrait of the zero dynamics in (6.2.12) is shown in Fig.6.5. It can be seen that the
equilibrium pθ1, 9θ1q “ p0, 0q is a center, surrounded by a dense set of closed orbits. We choose the
desired orbit Od to be the one that passes through pθ1, 9θ1q “ p0.0, 3.0q.

123

5

)
s
/
d
a
r
(

1
9θ

0

-5

-3

O

d

0

θ1 (rad)

3

Figure 6.5: Phase portrait of tiptoebot zero dynamics. The desired orbit Od is shown in red.

6.5.3 Stabilization of Od

We define the Poincaré section of Od as follows

Σ “ tx P Q

3 ˆ R

3

: θ1 “ 0,

9θ1 ě 0u

(6.5.5)

The states on Σ are

z “ rθ2 θ3

9θ1

9θ2

9θ3sT

Substituting pθ1, 9θ1q “ p0.0, 3.0q in (6.5.3) and its derivative gives
z˚ “ r0.0 0.0 3.0 ´6.0 0.3sT

(6.5.6)

Using the values of ε1 “ 0.01 and ε2 “ 0.004, the matrices A and B in (6.3.15) and (6.3.16) were
obtained as

0.800

0.050

1.530

´0.380 ´0.080
0.000 ´0.460 ´0.080 ´0.003
2.770
1.230

1.890

6.120

0.730

4.050

´3.210 ´3.770 ´13.360 ´6.090 ´8.100
0.120 ´0.560
0.100

0.280

0.670

1.525 ´3.700 ´17.700
4.875 ´8.650

22.650 ´43.850 ´0.325

34.325

0.875

124

A “

»——————————–
B “»—–

fiffiffiffiffiffiffiffiffiffiffifl
fiffifl

T

All eigenvalues of A do not lie inside the unit circle but the pair tA , Bu is controllable and satisfy
Proposition 2. Using LQR, the gain matrix K in (6.3.14) is obtained as

K “»—–

fiffifl

0.028

0.024 0.197

0.094

0.138

´0.034 ´0.051 0.116 ´0.049 ´0.055

(6.5.7)

The eigenvalues of pA ` BK q are located at 0.14, ´0.47˘ 0.73i and ´0.12˘ 0.56i; thus Od is
exponentially stable under the impulsive feedback.

6.5.4 Simulation Results

The initial configuration of the tiptoebot is taken as

rθ1 θ2 θ3

9θ1

9θ2

9θ3s “ r´0.1 0.2 0.05 3.3 ´6.0 0.4s

For the controller gains in (6.5.4) and (6.5.7), simulation results of the ICPM approach are shown

in Fig.6.6. The plots of ρ1, ρ2, 9ρ1 and 9ρ2 with time are shown in Figs.6.6 (a)-(d); it can be seen

that the continuous controller uc in (6.2.9) stabilizes the constraint manifold C . The phase portrait
of the passive joint θ1 is shown in Fig.6.6 (e) for 0 ď t ď 20 s and Fig.6.6 (f) for t ą 20 s. To
asymptotically stabilize the desired orbit Od, shown in red in both Figs.6.6 (e) and (f), the impulsive

controller in (6.3.14) is implemented using the high-gain feedback in (6.3.19); Λ was chosen to be

an identity matrix and µ was chosen as 0.0001. To show the convergence of system trajectories

to Od, }epkq}2 is plotted with respect to k in Fig.6.6 (e). It can be seen that for large values of k,
}epkq}2 Ñ 0; this implies that Od is exponentially stable.

125

0.6

0.0

-0.6

0.5

0.0

-0.5

0

(a)

(c)

ρ1

9ρ1

0.45

0.00

-0.45

0.5

0.0

(b)

(d)

ρ2

9ρ2

-0.5

0

100

time (s)

(e)

time (s)

(f)

100

4

t ď 20

t ą 20

1
9θ

-4

-2

3.5

0.0

0

θ1

2

-2

(g)

θ1

2

}epkq}2

k

40

Figure 6.6: Orbital stabilization for the tiptoebot using ICPM: (a) and (b) provide plots of ρ1 and
ρ2, (c) and (d) provide plots of 9ρ1 and 9ρ2, (e) and (f) provide the phase portrait of the passive joint,
(g) provides the norm of the error for the discrete-time system.

126

CHAPTER 7

CONCLUSION AND FUTURE WORK

In this work we investigated several problems where impulsive inputs were exploited to stabilize

equlibria and orbits of underactuated robotic systems. In chapter 2, using both simulations and

experiments, we demonstrated a method for stabilizing the equilibria of underactuated systems

from configurations lying outside the estimate of their region of attraction. The method, known

as the IMM, uses impulsive inputs to force the configuration of the system to move inside the

region of attraction along a manifold of dimension equal to the number of active degrees-of-

freedom. Justified by singular perturbation theory [23], the impulsive inputs were implemented

using high-gain feedback.

In the general case, the region of attraction is not known and the

IMM requires an estimate of the region, which is computed off-line. If the initial configuration

of the system is outside this estimate, the IMM can be used if the impulse manifold intersects

this region. The likelihood of intersection improves if the size of the estimated region is large,

and therefore, a procedure for obtaining large estimates was developed by combining the SOS

and trajectory reversing methods. For a fixed order of the stabilizing controller, the SOS method

maximizes the estimate of the region of attraction and this estimate is further enlarged using an

algorithm (CHART) based on the method of trajectory reversing. The advantage of combining the

two methods was demonstrated in simulations using a three-link underactuated system (Tiptoebot)

and in experiments with a pendubot. For both systems, it was shown that the impulse manifold

may intersect the enlarged estimate but not the estimate obtained using the SOS method, implying

that the CHART enhances the utility of the IMM. In another simulation with the Tiptoebot, it was

shown that the combination of SOS, CHART and IMM requires a lower magnitude of the impulsive

input compared to the combination of SOS and IMM. In addition to providing a large estimate, the

CHART makes it possible to generalize the results to higher-dimensional systems. In the literature,

several algorithms have been proposed for estimating the region of attraction using trajectory

reversing, but the CHART is the only algorithm that has been demonstrated for a system with as

127

many as six states. Earlier experimental results [23] relied on finding the boundary of the region of

attraction through trial and error and used the IMM to move the configuration of the system from a

point lying immediately outside the region to inside the region. The experimental results presented

here are more meaningful and practical as they are based on an estimate of the region of attraction

and points chosen outside this region are not very close to the boundary. Experimental results of

stabilization from two different initial conditions using the SOS, CHART and IMM were presented.

In one of the experiments, through trial and error, the initial condition was purposely chosen to

lie outside the region of attraction; it was shown that the equilibrium could not be stabilized with

the SOS controller alone but was stabilized by the combination of SOS, CHART and IMM. Future

work will investigate the possibility of extending the approach to multiple impulsive inputs. This

will be useful in situations where the impulse manifold does not intersect the estimated region of

attraction: a first impulsive input may be used to take the configuration closer to the boundary of

the estimate and additional impulses may be applied to move the system configuration inside the

region in future time. Of course, this would require real-time computation of the intersection of

the impulse manifold with the boundary of the estimated region. Compared to the SOS method,

which provides an analytical representation of the estimate of the region of attraction, the trajectory

reversing method provides a numerical representation and will involve higher computational cost.

Considering the fact that the trajectory reversing method provides the opportunity to enlarge the

estimate, it is necessary to investigate the trade-off between computational cost and size of the

estimate for real-time applications.

In chapter 3, rest-to-rest maneuvers of the inertia-wheel pendulum was studied in the framework

of impulsive control. Assuming a set of discrete impulsive inputs, optimal sequences were designed

to minimize their infinity-norm. It was shown analytically that a sequence with an odd number

of inputs is less optimal than the two adjacent sequences with even number of inputs. Analytical

and simulation results with two inputs were used to explain the high wheel velocities and large

continuous torques associated with methods that attempt to take the pendulum directly to its

desired configuration and minimize overshoot. Implementation of impulsive control using high-

128

gain feedback also results in large torques but they act over short intervals of time; therefore,

feasibility of impulsive control is determined by the peak torque rating of the actuator, which is

always larger than the continuous torque rating. It was shown that the number of impulsive inputs

can be increased to not exceed the peak torque rating of the actuator; this, of-course, increases

the time required for swing-up. Simulation results for swing-up showed similarities between the

optimal trajectories and the trajectories obtained using the energy-based controllers. Future work

will investigate the possibility of extending the method presented here to other underactuated

systems.

Using impulsive forces only, the problem of juggling a devil-stick was presented in chapter

4.

Impulsive forces were applied intermittently for juggling the stick between two symmetric

configurations. A dynamic model of the devil-stick and a control design for the juggling task was

presented for the first time. The control inputs are the impulse of the impulsive force and its point of

application on the stick. The control action is event-based and the inputs are applied only when the

stick has the orientation of one of the two symmetric configurations. The dynamics of the devilstick

due to the control action and torque-free motion under gravity is described by two Poincaré sections;

the symmetric configurations are fixed points of these sections. A coordinate transformation is

used to exploit the symmetry and convert the problem into that of stabilization of a single fixed

point. A dead-beat controller is designed to convert the nonlinear system into a controllable linear

discrete-time system with input constraints. LQR and MPC methods are used to design the control

inputs and achieve symmetric juggling. The LQR method has a closed-form solution and is easier

to implement but requires trial and error to satisfy the input constraints. The MPC method has no

closed-form solution as it is obtained by solving an optimization problem online. However, the

optimization problem directly takes into account the input constraint. The computational cost of the

MPC method, which can be a concern for many problems, is not a concern for the juggling problem

since the time between consecutive control actions is relatively large.Simulation results validate

both control designs and demonstrate non-prehensile manipulation solely using impulsive forces.

Our future work will focus on robotic juggling; this includes design of experimental hardware,

129

feedback compensation of energy losses due to inelastic collisions between the devil-stick and hand

sticks, and motion planning and control of the robot end-effector for generating the impulsive forces

designed by the control algorithms.

In chapter 5, a hybrid control strategy was presented for orbital stabilization of underactuated

systems with one passive DOF. The orbit is defined with the help of a Lyapunov-like function

that has been commonly used in the literature. Unlike existing energy-based methods, that have

relied on continuous control inputs alone, our control strategy uses continuous control inputs and

intermittent impulsive brakings. The continuous control is designed to make the time derivative

of the Lyapunov-like function negative semi-definite. When this condition cannot be enforced, the

impulsive inputs are invoked. This results in negative jumps in the Lyapunov-like function and

guarantees its negative semi-definiteness under continuous control for some finite time interval.

Thus, a combination of continuous and impulsive inputs guarantees monotonic convergence of

the system trajectories to the desired orbit, which can be periodic, or non-periodic as in the case

of homoclinic orbits, depending on the choice of desired energy. More importantly, it allows

us to develop a general framework for energy-based orbital stabilization, which is an important

contribution of this work. A set of conditions, that impose constraints on the choice of controller

gains, have to be satisfied for applicability of the control strategy. These conditions are easily

satisfied by systems commonly studied in the literature such as the pendubot, acrobot, inertia-wheel

pendulum, and pendulum on a cart. In this work, the hybrid control strategy was demonstrated

in a three-DOF underactuated system using simulations and the two-DOF rotary pendulum using

experiments.

In experiments, impulsive brakings were not applied by the motor; instead, they

were applied by a friction brake mounted co-axially with the motor shaft. This requires additional

hardware but there are two important advantages of using the brake. In physical systems, impulsive

inputs are implemented using high-gain feedback, which can result in actuator saturation. Since our

impulsive control strategy requires the active velocities to be reduced to zero, a brake is a natural

choice and it eliminates the possibility of motor torque saturation. The advantage of using a brake

is also manifested in the time required for orbital stabilization. A comparison of our approach with

130

an approach in the literature shows significant reduction in the time for convergence for the same

set of initial conditions.

Finally, the control design for stabilization of VHC based orbits was presented in chapter 6.

Repetitive motion in underactuated systems are typically designed using VHCs. A VHC results

in a family of periodic orbits and stabilization of a desired orbit is an important problem. A

hybrid control design is presented to stabilize a VHC-generated periodic orbit for underactuated

systems with one passive DOF. A continuous controller is first designed to enforce the VHC and

stabilize the corresponding constraint manifold. The desired orbit on the constraint manifold is

exponentially stabilized by applying the continuous inputs together with impulsive inputs that are

periodically applied on a Poincaré section. The impulsive inputs alter the Poincaré map and this

impulse controlled Poincaré map (ICPM) is described by a discrete LTI system. The problem of

orbital stabilization is thus simplified to exponential stabilization of the fixed point of the ICPM. The

controllability of the system can be easily verified and the control design can be easily carried out

using standard techniques such as pole-placement and LQR. The identification of the linear system

and computation of the controller gains are performed off-line. The complexity and computational

cost of the ICPM approach is less than existing methods as it eliminates the need for on-line solution

of a periodic Ricatti equation. The ICPM approach is demonstrated using standard cart-pendulum

system; its applicability to higher-dimensional systems is demonstrated using the tiptoebot. Future

work will focus on gait stabilization of legged robots and experimental validation.

131

BIBLIOGRAPHY

132

BIBLIOGRAPHY

[1]

R. Olfati-Saber, Global stabilization of a flat underactuated system: the inertia wheel pen-
dulum, in: Proc. IEEE Conference on Decision and Control, Vol. 4, 2001, pp. 3764–3765.

[2] M. W. Spong, P. Corke, R. Lozano, Nonlinear control of the reaction wheel pendulum,

Automatica 37 (11) (2001) 1845–1851.

[3]

A. Shiriaev, J. W. Perram, C. Canudas-de Wit, Constructive tool for orbital stabilization
of underactuated nonlinear systems: Virtual constraints approach, IEEE Transactions on
Automatic Control 50 (8) (2005) 1164–1176.

[4] M. W. Spong, D. J. Block, The pendubot: A mechatronic system for control research and
education, in: Decision and Control, 1995., Proceedings of the 34th IEEE Conference on,
Vol. 1, IEEE, 1995, pp. 555–556.

[5] M. W. Spong, The swing up control problem for the acrobot, IEEE Control Systems 15 (1)

(1995) 49–55.

[6]

[7]

[8]

[9]

N. Kant, R. Mukherjee, D. Chowdhury, H. K. Khalil, Estimation of the region of attraction
of underactuated systems and its enlargement using impulsive inputs, IEEE Transactions on
Robotics 35 (3) (2019) 618–632.

F. Plestan, J. W. Grizzle, E. R. Westervelt, G. Abba, Stable walking of a 7-dof biped robot,
IEEE Transactions on Robotics and Automation 19 (4) (2003) 653–668.

F. B. Mathis, R. Mukherjee, Apex height control of a two-mass robot hopping on a rigid
foundation, Mechanism and Machine Theory 105 (2016) 44–57.

R. Ronsse, P. Lefevre, R. Sepulchre, Rhythmic feedback control of a blind planar juggler,
IEEE Transactions on Robotics 23 (4) (2007) 790–802.

[10] R. Goebel, R. G. Sanfelice, A. R. Teel, Hybrid dynamical systems, IEEE Control Systems

29 (2) (2009) 28–93.

[11] Y. Aoustin, D. T. Romero, C. Chevallereau, S. Aubin, Impulsive control for a thirteen-
link biped, in: 9th IEEE International Workshop on Advanced Motion Control, 2006, pp.
439–444.

[12] S. Weibel, G. W. Howell, J. Baillieul, Control of single-degree-of-freedom Hamiltonian
systems with impulsive inputs, in: Proc. 35th IEEE Conference on Decision and Control,
Vol. 4, 1996, pp. 4661–4666.

[13] A. M. Bloch, N. E. Leonard, J. E. Marsden, Controlled Lagrangians and the stabilization of
mechanical systems. i. the first matching theorem, IEEE Transactions on Automatic Control
45 (12) (2000) 2253–2270.

133

[14] A. M. Bloch, D. E. Chang, N. E. Leonard, J. E. Marsden, Controlled Lagrangians and the
stabilization of mechanical systems. ii. potential shaping, IEEE Transactions on Automatic
Control 46 (10) (2001) 1556–1571.

[15] R. Ortega, M. W. Spong, F. Gómez-Estern, G. Blankenstein, Stabilization of a class of un-
deractuated mechanical systems via interconnection and damping assignment, IEEE Trans-
actions on Automatic Control 47 (8) (2002) 1218–1233.

[16] D. Auckly, L. Kapitanski, W. White, Control of nonlinear underactuated systems, Commu-

nications on Pure and Applied Mathematics 53 (3) (2000) 354–369.

[17]

I. Fantoni, R. Lozano, M. W. Spong, Energy based control of the pendubot, IEEE Transactions
on Automatic Control 45 (4) (2000) 725–729.

[18] M. Zhang, T.-J. Tarn, Hybrid control of the pendubot, IEEE/ASME Transactions on Mecha-

tronics 7 (1) (2002) 79–86.

[19] T. Albahkali, R. Mukherjee, T. Das, Swing-up control of the pendubot: An impulse-

momentum approach, IEEE Transactions on Robotics 25 (4) (2009) 975–982.

[20] R. Jafari, F. B. Mathis, R. Mukherjee, Swing-up control of the acrobot: An impulse-

momentum approach, in: Proc. American Control Conference, 2011, pp. 262–267.

[21]

J. Lee, R. Mukherjee, H. K. Khalil, Output feedback stabilization of inverted pendulum on
a cart in the presence of uncertainties, Automatica 54 (2015) 146–157.

[22] F. Andreev, D. Auckly, S. Gosavi, L. Kapitanski, A. Kelkar, W. White, Matching, linear

systems, and the ball and beam, Automatica 38 (12) (2002) 2147 – 2152.

[23] R. Jafari, F. B. Mathis, R. Mukherjee, H. Khalil, Enlarging the region of attraction of
equilibria of underactuated systems using impulsive inputs, IEEE Transactions on Control
Systems Technology 24 (1) (2016) 334–340.

[24] F. B. Mathis, R. Jafari, R. Mukherjee, Efficient swing-up of the acrobot using continuous
torque and impulsive braking, in: Proc. American Control Conference, 2011, pp. 268–273.

[25] F. B. Mathis, R. Jafari, R. Mukherjee, Impulsive actuation in robot manipulators: Experi-
mental verification of pendubot swing-up, IEEE/ASME Transactions on Mechatronics 19 (4)
(2014) 1469–1474.

[26] D. Davidson, S. A. Bortoff, Enlarge your region of attraction using high-gain feedback, in:

Proc. IEEE Conference on Decision and Control, 1994, pp. 634–639.

[27] C. Reboulet, C. Champetier, A new method for linearizing non-linear systems: the pseudo-

linearization, International Journal of Control 40 (4) (1984) 631–638.

[28] T. Hu, Z. Lin, Composite quadratic Lyapunov functions for constrained control systems,

IEEE Transactions on Automatic Control 48 (3) (2003) 440–450.

134

[29] R. Bakhtiari, M. Yazdanpanah, Designing a linear controller for polynomial systems with the
largest domain of attraction via LMIs, in: Proc. IEEE International Conference on Control
and Automation, Vol. 1, 2005, pp. 449–453.

[30] A. I. Zečević, D. D. Šiljak, Estimating the region of attraction for large-scale systems with

uncertainties, Automatica 46 (2) (2010) 445–451.

[31] P. A. Parrilo, Structured semidefinite programs and semialgebraic geometry methods in

robustness and optimization, Ph.D. thesis, California Institute of Technology (2000).

[32] W. Tan, A. Packard, Stability region analysis using polynomial and composite polynomial
Lyapunov functions and sum-of-squares programming, IEEE Transactions on Automatic
Control 53 (2) (2008) 565–571.

[33] U. Topcu, A. Packard, P. Seiler, Local stability analysis using simulations and sum-of-squares

programming, Automatica 44 (10) (2008) 2669–2675.

[34] U. Topcu, A. K. Packard, P. Seiler, G. J. Balas, Robust region-of-attraction estimation, IEEE

Transactions on Automatic Control 55 (1) (2010) 137–142.

[35] A. Majumdar, A. A. Ahmadi, R. Tedrake, Control design along trajectories with sums of
squares programming, in: Proc. IEEE International Conference on Robotics and Automation,
2013, pp. 4054–4061.

[36] A. Majumdar, R. Tedrake, Funnel libraries for real-time robust feedback motion planning,

The International Journal of Robotics Research 36 (8) (2017).

[37] R. Genesio, A. Vicino, New techniques for constructing asymptotic stability regions for

nonlinear systems, IEEE Transactions on Circuits and Systems 31 (6) (1984) 574–581.

[38] R. Genesio, M. Tartaglia, A. Vicino, On the estimation of asymptotic stability regions:
State of the art and new proposals, IEEE Transactions on Automatic Control 30 (8) (1985)
747–755.

[39] M. Loccufier, E. Noldus, A new trajectory reversing method for estimating stability regions

of autonomous nonlinear systems, Nonlinear dynamics 21 (3) (2000) 265–288.

[40] F. Hamidi, H. Jerbi, W. Aggoune, M. Djemai, M. N. Abdkrim, Enlarging region of attraction
via LMI-based approach and genetic algorithm, in: Proc. IEEE International Conference on
Communications, Computing and Control Applications, 2011, pp. 1–6.

[41] D. Chowdhury, N. Kant, R. Mukherjee, H. K. Khalil, Enlarging the region of attraction of
equilibria of underactuated systems using sum of squares and impulse manifold method, in:
Proc. American Control Conference, 2017, pp. 893–898.

[42] N. Kant, D. Chowdhury, R. Mukherjee, H. K. Khalil, An algorithm for enlarging the region
of attraction using trajectory reversing, in: Proc. American Control Conference, 2017, pp.
4171–4176.

135

[43] M. Ryalat, D. S. Laila, A simplified ida-pbc design for underactuated mechanical systems

with applications, European Journal of Control 27 (2016) 1–16.

[44] N. Qaiser, N. Iqbal, A. Hussain, Stabilization of non-linear inertia wheel pendulum system
using a new dynamic surface control based technique, in: IEEE International Conference on
Engineering of Intelligent Systems, 2006, pp. 1–6.

[45] A. Khoroshun, Stabilization of the upper equilibrium position of a pendulum by spinning an

inertial flywheel, International applied mechanics 52 (5) (2016) 547–556.

[46] A. Zhang, C. Yang, S. Gong, J. Qiu, Nonlinear stabilizing control of underactuated inertia
wheel pendulum based on coordinate transformation and time-reverse strategy, Nonlinear
Dynamics 84 (4) (2016) 246702476.

[47] HNWpodcasts,

Five

easy

devil

stick

https://www.youtube.com/watch?v=Z7Pv-p-nEYo,
2020] (2013).

tricks
[Online;

with

instructions,
accessed 01-May-

[48] N. Duinker,

Devilstick

tutorial

-

basic

tricks

to

https://www.youtube.com/watch?v=MAuAtwZ7BF4,
2020] (2015).

[Online;

you

get
started,
accessed 01-May-

[49] H. Hirai, F. Miyazaki, Dynamic coordination between robots: Self-organized timing selec-
tion in a juggling-like ball-passing task, IEEE Transactions on Systems, Man, and Cyber-
netics, Part B (Cybernetics) 36 (4) (2006) 738–754.

[50]

J. Kober, M. Glisson, M. Mistry, Playing catch and juggling with a humanoid robot, in: 12th
IEEE-RAS International Conference on Humanoid Robots (Humanoids 2012), IEEE, 2012,
pp. 875–881.

[51] A. Tornambe, Modeling and control of impact in mechanical systems: Theory and experi-

mental results, IEEE Transactions on Automatic Control 44 (2) (1999) 294–309.

[52] K. M. Lynch, C. K. Black, Recurrence, controllability, and stabilization of juggling, IEEE

Transactions on Robotics and Automation 17 (2) (2001) 113–124.

[53] A. Zavala-Rio, B. Brogliato, Direct adaptive control design for one-degree-of-freedom

complementary-slackness jugglers, Automatica 37 (7) (2001) 1117–1123.

[54] B. Brogliato, M. Mabrouk, A. Z. Rio, On the controllability of linear juggling mechanical

systems, Systems & control letters 55 (4) (2006) 350–367.

[55] R. Goebel, R. G. Sanfelice, A. R. Teel, Hybrid Dynamical Systems: modeling, stability, and

robustness, Princeton University Press, 2012.

[56] M. W. Spong, Impact controllability of an air hockey puck, Systems & Control Letters 42 (5)

(2001) 333–345.

[57] S. Schaal, C. G. Atkeson, Open loop stable control strategies for robot juggling, in: Proc.

IEEE International Conference on Robotics and Automation, IEEE, 1993, pp. 913–918.

136

[58] R. G. Sanfelice, A. R. Teel, R. Sepulchre, A hybrid systems approach to trajectory tracking
control for juggling systems, in: Decision and Control, 2007 46th IEEE Conference on,
IEEE, 2007, pp. 5282–5287.

[59] S. Nakaura, Y. Kawaida, T. Matsumoto, M. Sampei, Enduring rotary motion control of devil

stick, IFAC Proceedings Volumes 37 (13) (2004) 805–810.

[60] A. Shiriaev, A. Robertsson, L. Freidovich, R. Johansson, Generating stable propeller mo-
tions for devil stick, in: 3rd IFAC Workshop on Lagrangian and Hamiltonian Methods for
Nonlinear Control, 2006, 2006.

[61]

J. Zhao, M. W. Spong, Hybrid control for global stabilization of the cart–pendulum system,
Automatica 37 (12) (2001) 1941–1951.

[62] M. Maggiore, L. Consolini, Virtual holonomic constraints for Euler- Lagrange systems,

IEEE Transactions on Automatic Control 58 (4) (2013) 1001–1008.

[63] A. Mohammadi, M. Maggiore, L. Consolini, Dynamic virtual holonomic constraints for
stabilization of closed orbits in underactuated mechanical systems, Automatica 94 (2018)
112–124.

[64] E. R. Westervelt, C. Chevallereau, J. H. Choi, B. Morris, J. W. Grizzle, Feedback Control of

Dynamic Bipedal Robot Locomotion, CRC press, 2007.

[65] A. S. Shiriaev, L. B. Freidovich, A. Robertsson, R. Johansson, A. Sandberg, Virtual-
holonomic-constraints-based design of stable oscillations of Furuta pendulum: Theory and
experiments, IEEE Transactions on Robotics 23 (4) (2007) 827–832.

[66] L. Freidovich, A. Robertsson, A. Shiriaev, R. Johansson, Periodic motions of the pendubot
via virtual holonomic constraints: Theory and experiments, Automatica 44 (3) (2008) 785–
791.

[67] X. Xin, T. Yamasaki, Energy-based swing-up control for a remotely driven acrobot: Theo-
retical and experimental results, IEEE Transactions on Control Systems Technology 20 (4)
(2012) 1048–1056.

[68]

I. Fantoni, R. Lozano, M. Spong, Stabilization of the reaction wheel pendulum using an
energy approach, in: European Control Conference (ECC), IEEE, 2001, pp. 2552–2557.

[69] R. Lozano, I. Fantoni, D. J. Block, Stabilization of the inverted pendulum around its homo-

clinic orbit, Systems & Control Letters (2000) 197–204.

[70] H. Oka, Y. Maruki, H. Suemitsu, T. Matsuo, Nonlinear control for rotational movement of
cart-pendulum system using homoclinic orbit, International Journal of Control, Automation
and Systems 14 (5) (2016) 1270–1279.

[71]

I. Fantoni, R. Lozano, Stabilization of the Furuta pendulum around its homoclinic orbit,
International Journal of Control 75 (6) (2002) 390–398.

137

[72] X. Xin, M. Kaneda, Swing-up control for a 3-dof gymnastic robot with passive first joint:

design and analysis, IEEE Transactions on Robotics 23 (6) (2007) 1277–1285.

[73] H. K. Khalil, Nonlinear control, Prentice Hall, 2014.

[74] N. Kant, R. Mukherjee, H. K. Khalil, Stabilization of homoclinic orbits of two degree-of-

freedom underactuated systems, in: 2019 American Control Conference (ACC), 2019.

[75]

[76]

J. W. Grizzle, C. Chevallereau, R. W. Sinnet, A. D. Ames, Models, feedback control, and
open problems of 3d bipedal robotic walking, Automatica 50 (8) (2014) 1955–1988.

J. W. Grizzle, G. Abba, F. Plestan, Asymptotically stable walking for biped robots: Analysis
via systems with impulse effects, IEEE Transactions on automatic control 46 (1) (2001)
51–64.

[77] C. Canudas-de Wit, On the concept of virtual constraints as a tool for walking robot control

and balancing, Annual Reviews in Control 28 (2) (2004) 157–166.

[78] C. Chevallereau, G. Abba, Y. Aoustin, F. Plestan, E. R. Westervelt, C. Canudas-De-Wit, J. W.
Grizzle, RABBIT: A testbed for advanced control theory, IEEE Control Systems Magazine
23 (5) (2003) 57–79.

[79] C. Chevallereau, J. W. Grizzle, C.-L. Shih, Asymptotically stable walking of a five-link

underactuated 3-D bipedal robot, IEEE Trans. on Robotics 25 (1) (2009) 37–50.

[80] E. R. Westervelt, J. W. Grizzle, D. E. Koditschek, Hybrid zero dynamics of planar biped

walkers, IEEE Transactions on Automatic Control 48 (1) (2003) 42–56.

[81] C. Canudas-de Wit, B. Espiau, C. Urrea, Orbital stabilization of underactuated mechanical

systems, in: Proc. 15th IFAC World Congress, 2002.

[82] A. Mohammadi, E. Rezapour, M. Maggiore, K. Y. Pettersen, Maneuvering control of planar
snake robots using virtual holonomic constraints, IEEE Trans. on Control Systems Technol-
ogy 24 (3) (2015) 884–899.

[83] A. Shiriaev, A. Robertsson, J. Perram, A. Sandberg, Periodic motion planning for virtually

constrained Euler-Lagrange systems, Systems & Control Letters 55 (11) (2006) 900–907.

[84] A. S. Shiriaev, L. B. Freidovich, S. V. Gusev, Transverse linearization for controlled mechan-
ical systems with several passive degrees of freedom, IEEE Trans. on Automatic Control
55 (4) (2010) 893–906.

[85] A. Mohammadi, Virtual holonomic constraints for Euler-Lagrange control systems, Ph.D.

thesis, University of Toronto (2016).

[86] L. L. Flynn, R. Jafari, R. Mukherjee, Active synthetic-wheel biped with torso, IEEE Trans-

actions on Robotics 26 (5) (2010) 816–826.

[87] A. Majumdar, A. A. Ahmadi, R. Tedrake, Control and verification of high-dimensional
systems with DSOS and SDSOS programming, in: Proc. 53rd IEEE Conference on Decision
and Control, 2014, pp. 394–401.

138

[88] C. B. Barber, D. P. Dobkin, H. Huhdanpaa, The quickhull algorithm for convex hulls, ACM

Transactions on Mathematical Software 22 (4) (1996) 469–483.

[89]

J. D’Errico, Inhull, mathworks.com Updated on 06 Sep 2012 (March 2006).

[90] M. J, Analyze n-dimensional polyhedra in terms of vertices or (In)Equalities, mathworks.com

Updated on 16 Sep 2015 (March 2011).

[91] M. Ester, H.-P. Kriegel, J. Sander, X. Xu, et al., A density-based algorithm for discovering

clusters in large spatial databases with noise, in: Kdd, Vol. 96, 1996, pp. 226–231.

[92] Yarpiz, Dbscan clustering algorithm, mathworks.com 06 Sep 2015 (March 2011).

[93] N. Kant, R. Mukherjee, H. K. Khalil, Swing-up of the inertia wheel pendulum using im-
pulsive torques, in: 2017 IEEE 56th Annual Conference on Decision and Control (CDC),
IEEE, 2017, pp. 5833–5838.

[94] H. A. Toliyat, G. B. Kliman, Handbook of electric motors, Vol. 120, CRC press, 2004.

[95] N. Kant, R. Mukherjee, Impulsive dynamics and control of the inertia-wheel pendulum,

IEEE Robotics and Automation Letters 3 (4) (2018) 3208–3215.

[96] S. Wiggins, Introduction to applied nonlinear dynamical systems and chaos, Vol. 2, Springer

Science & Business Media, 2003.

[97]

J. Guckenheimer, K. Hoffman, W. Weckesser, The forced Van der Pol equation i: The slow
flow and its bifurcations, SIAM Journal on Applied Dynamical Systems 2 (1) (2003) 1–35.

[98] K. Bold, C. Edwards, J. Guckenheimer, S. Guharay, K. Hoffman, J. Hubbard, R. Oliva,
W. Weckesser, The forced Van der Pol equation ii: Canards in the reduced system, SIAM
Journal on Applied Dynamical Systems 2 (4) (2003) 570–608.

[99] P. J. Antsaklis, A. N. Michel, A Linear Systems Primer, Birkhauser, 2007.

[100] L. Wang, Model predictive control system design and implementation using MATLAB®,

Springer Science & Business Media, 2009.

[101] W. M. Haddad, V. Chellaboina, S. G. Nersesov, Impulsive and hybrid dynamical systems,

Princeton Series in Applied Mathematics (2006).

[102] B. Brogliato, R. Lozano, B. Maschke, O. Egeland, Dissipative Systems Analysis and Control,

Springer, 2020.

[103] D. Liu, W. Guo, J. Yi, D. Zhao, Passivity-based-control for a class of underactuated mechani-
cal systems, in: IEEE International Conference on Intelligent Mechatronics and Automation,
2004, pp. 50–54.

[104] B. Brogliato, Nonsmooth Mechanics. Models, Dynamics and Control, Springer Nature

Switzerland AG, 2020.

139

[105] B. Brogliato, S. . Niculescu, P. Orhant, On the control of finite-dimensional mechanical
systems with unilateral constraints, IEEE Transactions on Automatic Control 42 (2) (1997)
200–215.

[106] R. I. Leine, N. Van de Wouw, Stability and Convergence of Mechanical Systems with

Unilateral Constraints, Vol. 36, Springer Science & Business Media, 2007.

[107] D. D. Bainov, P. S. Simeonov, Systems with Impulse Effect: Stability, Theory, and Applica-

tions, Ellis Horwood, 1989.

[108] R. G. Sanfelice, R. Goebel, A. R. Teel, Invariance principles for hybrid systems with con-
nections to detectability and asymptotic stability, IEEE Transactions on Automatic Control
52 (12) (2007) 2282–2297.

[109] R. Ortega, J. A. L. Perez, P. J. Nicklasson, H. J. Sira-Ramirez, Passivity-based control
of Euler-Lagrange systems: Mechanical, Electrical and Electromechanical Applications,
Springer Science & Business Media, 2013.

[110] L. Consolini, A. Costalunga, M. Maggiore, A coordinate-free theory of virtual holonomic

constraints, Journal of Geometric Mechanics 10 (4) (2018) 467–502.

[111] N. Nekhoroshev, The Poincaré-Lyapunov-Liouville-Arnol’d theorem, Functional Analysis

and Its Applications 28 (1994) 128–129.

[112] H. K. Khalil, Nonlinear Systems, 2nd Edition, Prentice Hall, Upper Saddle River, NJ, 1996.

[113] S. H. Strogatz, Nonlinear dynamics and chaos: With applications to physics, biology,

chemistry, and engineering, CRC Press, 2018.

[114] E. Ott, C. Grebogi, J. A. Yorke, Controlling chaos, Physical Review Letters 64 (11) (1990)

1196.

[115] J. Grizzle, E. Westervelt, C. Canudas-de Wit, Event-based PI control of an underactuated
biped walker, in: 42nd Int. Conference on Decision and Control, Vol. 3, IEEE, 2003, pp.
3091–3096.

[116] P. Kokotović, H. K. Khalil, J. O’reilly, Singular perturbation methods in control: analysis

and design, SIAM, 1999.

140