ESE 547 F’22
Problem Set 2
Computing and Visualizing Attractors and Basins via
Lyapunov Functions
D. E. Koditschek
February 13, 2022
Released Tue., Jan. 18
Corrected, Wed., Feb. 2
Due Tue., Feb. 22 at 11:59pm
This problem set has three regular questions, each of which include some Extra Credit
components. These will add to the numerator of your score without adding to the denomina-
tor. Recall from the course description document that the five scheduled problem sets count
cumulatively toward 60% of your final grade. The breakdown of points for this specific prob-
lem set are posted in the rubric details up at the Canvas portal where you will be submitting
your work.
Some partial tables of notation are included at the end of this document for your con-
venience. Please note that the hyperref links in this pdf only work for items and equations
defined within this document. For live links back to the lecture notes (e.g., equation and
section references in the tables of notation), use the version of this problem set included in
the article formatted Lecture notes.
Please remember all work must be uploaded to the canvas portal as a typeset pdf docu-
ment in order to be graded.
1. Here, we reconsider the damped harmonic oscillator VF , fDHO(x) (2.1.7), characterized
by the parameter values m = k = b = 1 as introduced in PS 1 Exercise 1. Since this
problem set explores related properties in two different coordinate representations of the
same (two-dimensional) state space, it will be convenient to name two different copies
of the plane to distinguish physical (position and velocity) coordinates, x ∈ X := R2
from the abstract (real canonical) coordinates y ∈ Y := R2.
(a) Regenerate the ODE-solver from PS 1 Exercise 1-1a and compare the computed
change in total energy along traejctories over time, ηHO (2.1.15), with the power
function η˙HO (2.1.16) as follows.
Plot all output graphs for solutions generated using the eight different initial con-
ditions x0 ∈ { [cos θisin θi
]
: θi = iπ/4, i ∈ {0, . . . 7}} ⊂ X .
1
i. Present a “three dimensional plot” depicting the energy values, ηHO◦f tDHO
(
x0
)
,
as a curve of (nonnegative) “heights” over the statespace plane, X , on which
you superimpose the forward time orbits, fR+DHO
(
x0
)
. Your output will look
something along the lines of Lecture Notes Figure 2.7
ii. Superimpose on that previous output an additional set of eight “three dimen-
sional curves” that plot the values of the power function, η˙DHO ◦ f tDHO
(
x0
)
as
a curve of (nonpositive) “depths” below the statespace plane, X . Now your
output will juxtapose a graphic along the lines of Lecture Notes Figure 3.5 ly-
ing “below” the X -plane that appears to be supporting “above” it the previous
graphic (Lecture Notes Figure 2.7).
iii. Compute a numerical estimate of the rate of energy change
̂˙ηHO (tk) := 1
∆Tk
[
ηHO ◦ f tkDHO
(
x0
)− ηHO ◦ f tk−1DHO(x0)] , (1.1)
where k ranges over the numerical values of time that your simulation was
allowed to run for each IC , x0, and ∆Tk denotes the time step the solver used
to advance from the (k − 1)th to the kth time step of its computed solution.
Substitute ̂˙ηHO for η˙HO in the previous “three dimensional plot”.
iv. Plot the error between the closed form power function, η˙HO, and its numeri-
cally computed estimate, ̂˙ηHO,
ϵ (tk) := η˙HO ◦ f tkDHO
(
x0
)− ̂˙ηHO (tk)
in a new “three dimensional plot” where each error value is plotted as a “height”
over or under the statespace plane, X on which you superimpose the forward
time orbits.
v. Extra Credit
Explain why ηHO cannot be used in conjunction with Lyapunov’s theorem
(Lecture Notes Section 3.1.2.2) to mathematically prove that the origin of X
is asymptotically stable under the flow generated by fDHO.
Instead, use ηHO in conjunction with Lasalle’s theorem (Lecture Notes Section
3.1.2.2) to prove mathematically that the origin is asymptotically stable.
(b) Although Lord Kelvin “knew” [1] that the PD property of ηHO coupled with the
NSD property of η˙HO was sufficient to conclude the asymptotic stability of the
origin, it fell to Lasalle [2] to give a formal mathematical justification for this
conclusion.
Here, we will use Lyapunov’s theorem by appeal to his generalized concept of
energy [3] in the form of a different PD function that turns out to have a ND
derivative.
Specifically, we will introduce a PD (generalized energy) function
ηRC : Y → R : y 7→ 1
2
yTy = ∥y∥2/2 (1.2)
whose associated power function, η˙RC is ND along the motions of fRC (2.1.13),
and use the CC hRC to “pull it back” into the x-coordinate plane where it can be
compared to the physical total energy, ηHO.
2
i. Compute the closed form generalized power function, η˙RC and simplify terms
to show that it is a multiple of the squared Euclidean norm in the y-plane.
Now use the code of the previous question, Exercise 1-1a, to redo those same
steps (i. through iv.) to visualize ηRC & η˙RC in plots over the y-coordinate
plane, Y. Your plots should reveal that η˙RC lies strictly “beneath” the y-plane,
Y, except at the origin as required to be ND.
...
v. Now compute the conjugate of ηRC in physical coordinates x ∈ X = h−1RC(Y)
as
η˜RC : X → R : x 7→ ηRC ◦ hRC(x)
and use fDHO to compute its associated (generalized) power function
˙˜ηRC (x) = Dη˜RC (x) · fDHO(x)
Once again, repeat the same steps of Exercise 1-1a (i. through iv.) to visualize
η˜RC and ˙˜ηRC in physical coordinates, x ∈ X .
...
ix. Explain how your plots corroborate the claim that ηRC is a strict LF for
fDHO
x. Extra Credit
Verify mathematically that ηRC is really a strict LF for fDHO and use Lya-
punov’s Theorem to prove that the origin of X is asymptotically stable.
2. This problem aims to provide a computational example of a dynamical (i.e. non-
equilibrium) steady state behavior called a limit cycle that will become our archetype
for understanding and commanding gaits.
We start with a nonlinear modification of fRC in the y-coordinate plane, Y,
fnRC(y) := [σ(∥y∥)I2 + ωJ2] y (1.3)
where ω is a constant as in (2.1.13) but σ is the scalar function
σ(ρ) := 1− ρ2.
and we will use reuse ηRC (1.2), but this time in a composition that accords it a role
in the construction of a Lasalle function.
(a) Use the ODE solver of your choice to simulate the solutions of the dynamical
system (1.3) generated by fnRC from the eight different initial conditions
y0 ∈ {2γ
[
(−1)βα, (−1)β(1− α)]T : α, β ∈ {0, 1} & γ ∈ {−1, 1}} ⊂ Y
i. Plot the resulting trajectories as functions of time so that your outputs have
the general appearance of the pair of left hand side traces in the panels of
Figure 2.3.
3
We don’t yet have good physical interpretation of this system, nevertheless,
write a sentence or two that describes the empirical observation you would
make about the state trajectories in contrast to the state trajectories of the
flow generated by fRC.
ii. Plot the orbits of the flow through these eight initial conditions so that your
output has the general appearance of Figure 2.4.
Describe qualitatively in a sentence or two the nature of the forward limit set
of the entire state space plane, Y.
iii. Extra Credit
Plot the trajectories through these eight initial conditions as spatial curves of
two-dimensional solution vectors parametrized by time. The output should
look like a set of eight corkscrew curves in three-space that start on the state
space plane at time t = 0 and wind around the time axis toward some desti-
nation
(b) Given the generalized energy,
ηnRC(y) :=
1
2
[σ ◦ 2ηRC(y)]2 = 1
2
[
σ
(
yTy
)]2
, (1.4)
compute the associated generalized power function along the motions generated
by fnRC,
η˙nRC (y) = DηnRC (y) · fnRC(y).
Once again, use the code of the earlier question, Exercise 1-1a, to redo those same
steps (i. through iv.) to visualize ηnRC & η˙nRC in the real canonical coordinate
plane, y ∈ Y.
(c) Interpret the numerical evidence from your plots as suggesting that ηnRC is a
Lasalle function and use the associated theorem to conjecture a formula for the
forward limit set, f∞nRC
(Y).
(d) Extra Credit
Prove mathematically that ηnRC is a Lasalle function and use the associated the-
orem to furnish a formula for the forward limit set f∞nRC
(Y).
3. This problem aims to provide an advance computational familiarity with what we shall
come to call the “active damping stance mode controller” for vertical hopping following
a method presented in [4].
Specifically in Lecture Notes Section 5.3.2 (still in preparation for upload to the Canvas
site) we will introduce a controlled version of the damped harmonic oscillator (5.3.23)
fcDHO(x, τ) =
[
0 ω
−ω −2ζ¯ω
]
x+
[
0
τ/ω
]
where the parameter assignments, ω := 1 & ζ¯ := 12 put the unforced system (i.e.,
where τ ≡ 0) into correspondence with the instance of fDHO we have been exploring
numerically. Let the actuator input take the form (5.3.24) (again, to appear in Lecture
Notes Section 5.3.2)
τ(x) := ktx2∥x∥+ ϵ
4
where we will let kt := 1 (it can be an arbitrary positive gain) and ϵ := 1/100 (it can
be any “sufficiently small” positive offset). We will denote the resulting closed loop VF
fAD(x) := fcDHO(x, τ(x)).
Finally, define a Lasalle function candidate for fAD,
ηAD(x) :=
1
2
[
−ζ¯ + kt/2ω
2
∥x∥+ ϵ
]2
(a) Repeat the steps of Exercise 2-2a for fAD where the set of eight IC for this numerical
study is given as
x0 ∈ {ρ
[
(−1)βα, (−1)β(1− α)]T : α, β ∈ {0, 1} & ρ ∈ {ρ−, ρ+}} ⊂ X
and
ρ+ ∈ {ρ ∈ R+ : ζ¯ > kt/2ω
2
ρ+ ϵ
}; ρ− ∈ {ρ ∈ R+ : ζ¯ < kt/2ω
2
ρ+ ϵ
}.
(b) Given that these numerical simulations use physical coordinates, x ∈ X , describe
in a sentence or two the intuitive steady state behavior of the mass when subject
to the control τ(x) that yields the closed loop fAD.
(c) Repeat the steps of Exercise 2-2a for the unforced system, fcDHO(x, 0) with ζ¯ := 0
(which coincides exactly with fDHO for b := 0) for the same initial conditions, x0
as just above, substituting ηHO, the total mechanical energy, in place of ηAD.
Describe in a sentence or two the intuitive steady state behavior of the mass under
these conditions, making sure to explain how it is different from that produced by
fAD.
(d) Extra Credit
Show that ηAD is a Lasalle function for fAD and use Lasalle’s Theorem to give a
mathematical characterization of the forward limit set f∞AD
(X ).
5
2 Appendix - Tables of Notation
Abbreviated Terminology
Symbol Segment Equation Description
VF Section 2.1.1 (2.1.3) vector field
IC Section 2.1.1 (2.1.3) initial condition
CC Section 2.1.1 (2.1.11) general change of coordinates
CCL Section 2.1.1 (2.1.11) linear change of coordinates
FP Section 2.2.2 fixed point of a VF
PD Section 3.1.2 positive definite
Generic Symbols and Spaces
Symbol Segment Equation Description
:= Section 2.1.1 (2.1.1) term := definition of term
=: Section 2.1.1 (2.1.2) definition of term =: term
· Section 2.1.2 (2.1.9) matrix multiplication
fDHO
hdiag∼ fdiag Section 2.1.2 (2.1.11) fdiag conjugate to fDHO, via CC hdiag
Rn Section 2.1.2 space of n× 1 real vectors
Rn×m Section 2.1.2 space of n×m real matrices
Cn Section 2.1.2 space of n× 1 complex vectors
Cn×m Section 2.1.2 space of n×m complex matrices
ℜ Section 2.1.2 (2.1.13) real component of a complex number
ℑ Section 2.1.2 (2.1.13) real component of a complex number
◦ Section 2.1.3 (2.1.15) function composition
S1 Section 2.2.1 (2.2.1) Unit circle
TS Section 2.2.1 (2.2.1) Unit circle tangent space
Constants
Symbol Segment Equation Description
In Section 2.1.2 (2.1.13) n× n identity matrix
J2 Section 2.1.2 (2.1.13) 2× 2 skew symmetric matrix
e1, e2 Section 2.2.1 planar unit vectors
g Section 2.2.1 (2.2.4) gravitational constant
Variables & Functions
6
Symbol Segment Equation Description
x Section 2.1.1 prismatic joint variable
ΦD Section 2.1.1 (2.1.1) scalar damper force function
fD Section 2.1.1 (2.1.3) scalar damper vector field
f tD Section 2.1.1 (2.1.4) scalar damper flow
φHS Section 2.1.1 (2.1.5) Hooke’s spring potential
x Section 2.1.1 (2.1.6) prismatic tangent vector
ΞDHO Section 2.1.1 (2.1.6) damped harmonic oscillator acceleration
fDHO Section 2.1.1 (2.1.7) damped harmonic oscillator vector field
f tDHO Section 2.1.1 (2.1.10) damped harmonic oscillator flow
hdiag Section 2.1.1 (2.1.11) diagonalizing linear coordinate transformation
f˜DHO Section 2.1.1 (2.1.11) conjugate VF to fDHO
z Section 2.1.1 (2.1.11) diagonalized (complex) coordinates for x
y Section 2.1.1 (2.1.13) rational canonical (RC) coordinates for x
ηHO Section 2.1.1 (2.1.15) harmonic oscillator total energy
η˙HO Section 2.1.1 (2.1.16) damped harmonic oscillator power function
Variables & Functions, ctd.
Symbol Segment Equation Description
q Section 2.2.1 (2.2.1) revolute tangent vector
gR Section 2.2.1 (2.2.2) 1 DoF Revolute Kinematics
b Section 2.2.1 (2.2.2) Cartesian coordinates for pendulum point mass body
κR Section 2.2.1 (2.2.3) Kinetic Energy for 1 DoF Revolute Kinematics
φR Section 2.2.1 (2.2.4) Potential Energy for 1 DoF Revolute Kinematics
λR Section 2.2.1 (2.2.5) Lagrangian for 1 DoF Revolute Kinematics
References
[1] Sir William Thomson and Peter Guthrie Tait. Treatise on natural philosophy. University
of Cambridge, 1888. 2
[2] J. P. LaSalle. The Stability of Dynamical Systems. Society for Industrial Mathematics,
1976. 2
[3] Aleksandr Mikha˘ı Lyapunov. Problème général de la stabilité du mouvement. Princeton
University Press, 1949. 2
[4] A. De and D.E. Koditschek. Parallel composition of templates for tail-energized planar
hopping. In 2015 IEEE International Conference on Robotics and Automation (ICRA),
page 4562–4569. IEEE, May 2015. 4
7