---
abstract: 'Attitude control systems naturally evolve on nonlinear configuration spaces, such as ${\ensuremath{\mathsf{S}}}^2$ and ${\ensuremath{\mathsf{SO(3)}}}$. The nontrivial topological properties of these configuration spaces result in interesting and complicated nonlinear dynamics when studying the corresponding closed loop attitude control systems. In this paper, we review some global analysis and simulation techniques that allow us to describe the global nonlinear stable manifolds of the hyperbolic equilibria of these closed loop systems. A deeper understanding of these invariant manifold structures are critical to understanding the global stabilization properties of closed loop control systems on nonlinear spaces, and these global analysis techniques are applicable to a broad range of problems on nonlinear configuration manifolds.'
author:
- 'Taeyoung Lee, Melvin Leok, and N. Harris McClamroch[^1][^2][^3][^4] [^5]'
bibliography:
- 'CDC11.bib'
title: |
    **Stable Manifolds of Saddle Points\
    for Pendulum Dynamics on ${\ensuremath{\mathsf{S}}}^2 $ and ${\ensuremath{\mathsf{SO(3)}}}$**
---

Introduction
============

Global nonlinear dynamics of various classes of closed loop attitude control systems have been studied in recent years. An overview of results on attitude control of a rotating rigid body is given in  [@ChSaMcCSM11]. Closely related results on attitude control of a spherical pendulum (with attitude an element of the two-sphere ${\ensuremath{\mathsf{S}}}^2$) and of a 3D pendulum (with attitude an element of the special orthogonal group ${\ensuremath{\mathsf{SO(3)}}}$) are given in  [@ChMcIJRNC07; @ChMcBeAut08; @ChaMcCITAC09]. These and other similar publications address the global closed dynamics of smooth vector fields. Assuming that the closed loop vector field has an asymptotically stable equilibrium, as desired in attitude stabilization problems, additional hyperbolic equilibria necessarily exist. The domain of attraction of the asymptotically stable equilibrium is contained in the complement of the union of the stable manifolds of the hyperbolic equilibria. These geometric factors motivate the current paper, in which new analytical and computational results on the stable manifolds of the hyperbolic equilibria are obtained.

To make the development concrete, the presentation is built around two specific closed loop vector fields: one for the attitude dynamics of a spherical pendulum and one for the attitude dynamics of a 3D pendulum. In analyzing these two cases, we introduce new analytical and computational tools that are broadly applicable to studying the geometry of more general attitude control systems.

Spherical Pendulum
==================

A spherical pendulum is composed of a mass $m$ connected to a frictionless pivot by a massless link of length $l$. It is acts under uniform gravity, and it is subject to a control moment $u$. The configuration of a spherical pendulum is described by a unit-vector $q\in{\ensuremath{\mathbb{R}}}^3$, representing the direction of the link with respect to a reference frame.

Therefore, the configuration space is the two-sphere ${\ensuremath{\mathsf{S}}}^2=\{q\in{\ensuremath{\mathbb{R}}}^3\,|\, q\cdot q =1\}$. The tangent space of the two-sphere at $q$, namely ${\ensuremath{\mathsf{T}}}_q{\ensuremath{\mathsf{S}}}^2$, is the two-dimensional plane tangent to the unit sphere at $q$, and it is identified with ${\ensuremath{\mathsf{T}}}_q{\ensuremath{\mathsf{S}}}^2\simeq\{\omega\in{\ensuremath{\mathbb{R}}}^3\,|\, q\cdot \omega =0\}$, using the following kinematics equation: $$\begin{aligned}
\dot q= \omega\times q,\end{aligned}$$ where the vector $\omega\in{\ensuremath{\mathbb{R}}}^3$ represents the angular velocity of the link. The equation of motion is given by $$\begin{gathered}
\dot\omega = \frac{g}{l}q\times e_3 + \frac{1}{ml^2} u,\end{gathered}$$ where the constant $g$ is the gravitational acceleration, and the vector $e_3=[0,0,1]\in{\ensuremath{\mathbb{R}}}^3$ denotes the unit vector along the direction of gravity. The control moment at the pivot is denoted by $u\in{\ensuremath{\mathbb{R}}}^3$.

Control System
--------------

Several proportional-derivative (PD) type control systems have been developed on ${\ensuremath{\mathsf{S}}}^2$ in a coordinate-free fashion [@BulMurN95; @BulLew05]. Here, we summarize a control system that stabilizes a spherical pendulum to a fixed desired direction $q_d\in{\ensuremath{\mathsf{S}}}^2$.

Consider an error function on ${\ensuremath{\mathsf{S}}}^2$, representing the projected distance from the direction $q$ to the desired direction $q_d$, given by $$\begin{aligned}
\Psi(q,q_d)=1-q\cdot q_d.\end{aligned}$$ The derivative of $\Psi$ with respect to $q$ along the direction $\delta q =\xi\times q$, where $\xi\in{\ensuremath{\mathbb{R}}}^3$ and $\xi\cdot q=0$, is given by $$\begin{aligned}
{\ensuremath{\mathbf{D}}}_q \Psi(q,q_d)\cdot\delta q = -(\xi\times q)\cdot q_d = (q_d\times q)\cdot \xi.\end{aligned}$$ For positive constants $k_q,k_\omega$, the control input is chosen as: $$\begin{aligned}
u = ml^2(- k_\omega \omega - k_q q_d\times q -\frac{g}{l}q\times e_3).\end{aligned}$$ The corresponding closed loop dynamics are given by $$\begin{gathered}
\dot \omega = -k_\omega\omega - k_q q_d\times q,\label{eqn:dotw}\\
\dot q= \omega\times q.\label{eqn:dotq}\end{gathered}$$

This yields two equilibrium solutions: (i) the desired equilibrium $(q,\omega)=(q_d,0)$; (ii) additionally, there exists another equilibrium $(-q_d,0)$ at the antipodal point on the two-sphere.

It can be shown that the desired equilibrium is asymptotically stable by using the following Lyapunov function: $$\begin{aligned}
\mathcal{V} = \frac{1}{2} \omega\cdot\omega + k_q \Psi(q,q_d).\end{aligned}$$ In this paper, we analyze the local stability of each equilibrium by linearizing the closed loop dynamics to study the equilibrium structures more explicitly. In particular, we develop a coordinate-free form of the linearized dynamics of [(\[eqn:dotw\])]{}, [(\[eqn:dotq\])]{}, in the following section.

Linearization
-------------

A variation of a curve $q(t)$ on ${\ensuremath{\mathsf{S}}}^2$ is a family of curves $q^\epsilon(t)$ parameterized by $\epsilon\in{\ensuremath{\mathbb{R}}}$, satisfying several properties [@BulLew05]. It cannot be simply written as $q^\epsilon(t)=q(t)+\epsilon\delta q(t)$ for $\delta q(t)$ in ${\ensuremath{\mathbb{R}}}^3$, since in general, this does not guarantee that $q^\epsilon(t)$ lies in ${\ensuremath{\mathsf{S}}}^2$. In [@LeeLeoIJNME08], an expression for a variation on ${\ensuremath{\mathsf{S}}}^2$ is given in terms of the exponential map as follows: $$\begin{aligned}
q^\epsilon(t) = \exp(\epsilon \hat\xi(t))q(t),\label{eqn:qe}\end{aligned}$$ for a curve $\xi(t)$ in ${\ensuremath{\mathbb{R}}}^3$ satisfying $\xi(t)\cdot q(t)=0$ for all $t$. The *hat map* $\hat\cdot:{\ensuremath{\mathbb{R}}}^3\rightarrow{\ensuremath{\mathfrak{so}(3)}}$ is defined by the condition that $\hat x y =x\times y$ for any $x,y\in{\ensuremath{\mathbb{R}}}^3$. The resulting infinitesimal variation is given by $$\begin{aligned}
\delta q (t) = \frac{d}{d\epsilon}\bigg|_{\epsilon=0} q^\epsilon(t) = \xi(t)\times q(t). \label{eqn:delq}\end{aligned}$$ The variation of the angular velocity can be written as $$\begin{aligned}
\omega^\epsilon(t) = \omega(t) + \epsilon \delta \omega(t),\label{eqn:we}\end{aligned}$$ for a curve $\delta w(t)$ in ${\ensuremath{\mathbb{R}}}^3$ satisfying $q(t)\cdot w(t)=0$ for all $t$. Hereafter, we do not write the dependency on time $t$ explicitly.

The time-derivative of $\delta q$ can be obtained either from [(\[eqn:delq\])]{} or by substituting [(\[eqn:qe\])]{}, [(\[eqn:we\])]{} into [(\[eqn:dotq\])]{}, and considering the first order terms of $\epsilon$. In either case, we have $$\begin{aligned}
\delta\dot q= \dot\xi\times q + \xi\times(\omega\times q) = \delta\omega\times q + \omega\times(\xi\times q).\end{aligned}$$ Using the vector cross product identity $a\times (b\times c)=(a\cdot c)b-(a\cdot b) c$ for any $a,b,c\in{\ensuremath{\mathbb{R}}}^3$, this can be written as $$\begin{gathered}
\dot\xi\times q + (\xi\cdot q)w -(\xi\cdot\omega) q = \delta\omega\times q +(\omega\cdot q)\xi -(\omega\cdot\xi)q.\end{gathered}$$ Since $\xi\cdot q=0$, $\omega\cdot q=0$, this reduces to $$\begin{gathered}
\dot\xi\times q = \delta\omega\times q.\end{gathered}$$ Since both sides of the above equation are perpendicular to $q$, this is equivalent to $q\times(\dot\xi\times q) = q\times(\delta\omega\times q)$, which yields $$\begin{gathered}
\dot \xi - (q\cdot\dot\xi) q = q\times(\delta\omega\times q).\end{gathered}$$ Since $\xi\cdot q =0$, we have $\dot\xi\cdot q +\xi\cdot\dot q=0$. Using this, the above equation can be rewritten as $$\begin{aligned}
\dot \xi & = -(\xi\cdot(\omega\times q))q+ q\times(\delta\omega\times q)\nonumber\\
& = (qq^T\hat\omega) \xi +(I-qq^T)\delta\omega.\label{eqn:dotxi}\end{aligned}$$ This corresponds to the linearized equation of motion for [(\[eqn:dotq\])]{}. Similarly, by substituting [(\[eqn:delq\])]{}, [(\[eqn:we\])]{} into [(\[eqn:dotw\])]{}, we obtain $$\begin{aligned}
\delta\dot\omega &= -k_\omega\delta\omega -k_qq_d\times(\xi\times q)\nonumber\\
&=-k_\omega\omega +k_q\hat q_d \hat q\,\xi,\label{eqn:dotdelw}\end{aligned}$$ which is the linearized equation for [(\[eqn:dotw\])]{}.

Equations [(\[eqn:dotxi\])]{}, [(\[eqn:dotdelw\])]{} can be written in a matrix form as $$\begin{aligned}
\dot x =
\begin{bmatrix}\dot\xi \\ \delta\dot\omega \end{bmatrix}
=\begin{bmatrix} qq^T\hat\omega & I-qq^T\\k_q\hat q_d\hat q & -k_wI\end{bmatrix}
\begin{bmatrix}\xi \\ \delta\omega \end{bmatrix}=Ax,\label{eqn:xdot}\end{aligned}$$ where the state vector of the linearized controlled system is $x=[\xi;\delta\omega]\in{\ensuremath{\mathbb{R}}}^6$. A spherical pendulum has two degrees of freedom, but this linearized equation of motion evolves in ${\ensuremath{\mathbb{R}}}^6$ instead of ${\ensuremath{\mathbb{R}}}^4$. Since $q\cdot\omega=0$ and $q\cdot\xi=0$, we have the following two additional constraints on $\xi,\delta\omega$: $$\begin{aligned}
Cx=
\begin{bmatrix} q^T & 0 \\ -\omega^T\hat q & q^T\end{bmatrix}
\begin{bmatrix}\xi \\ \delta\omega \end{bmatrix}
=\begin{bmatrix} 0 \\ 0 \end{bmatrix}.\label{eqn:con}\end{aligned}$$ Therefore, the state vector $x$ should lie in the null space of the matrix $C\in{\ensuremath{\mathbb{R}}}^{2\times 4}$. However, this is not an extra constraint that should be imposed when solving [(\[eqn:xdot\])]{}. As long as the initial condition $x(0)$ satisfies [(\[eqn:con\])]{}, the structure of [(\[eqn:dotw\])]{}, [(\[eqn:dotq\])]{}, and [(\[eqn:xdot\])]{}, guarantees that the state vector $x(t)$ satisfies [(\[eqn:con\])]{} for all $t$, i.e. $\frac{d}{dt}C(t)x(t) =0$ for all $t\geq 0$ when $C(0)x(0)=0$. This means that the null space of $C$ is a flow-invariant subspace.

Equilibrium Solutions
---------------------

We choose the desired direction as $q_d=e_3$. The equilibrium solution $(q_d,0)=(e_3,0)$ is referred to as the hanging equilibrium, and the additional equilibrium solution $(-q_d,0)=(-e_3,0)$ is referred to as the inverted equilibrium. We study the eigen-structure of each equilibrium using the linearized equation [(\[eqn:xdot\])]{}. To illustrate the ideas, the controller gains are selected as $k_q=k_\omega=1$.

### Hanging Equilibrium

The eigenvalues $\lambda_i$, and the eigenvectors $v_i$ of the matrix $A$ at the hanging equilibrium $(e_3,0)$ are given by $$\begin{gathered}
\lambda_{1,2}=(-1\pm\sqrt{3}i)/2,\;
\lambda_{3,4}=\lambda_{1,2},\;\lambda_5=0,\;\lambda_6=-1,\\
v_{1,2}= e_1 + (-1\pm\sqrt{3}i)e_4/2,\;
v_{3,4}= e_2 + (-1\pm\sqrt{3}i)e_5/2,\\
v_5=e_3,\quad v_6=e_6,\end{gathered}$$ where $e_i\in{\ensuremath{\mathbb{R}}}^6$ denotes the unit-vector whose $i$-th element is one, and other elements are zeros. Note that there are repeated eigenvalues, but we obtain six linearly independent eigenvectors, i.e., the geometric multiplicities are equal to the algebraic multiplicities.

The basis of the null space of the matrix $C$, namely $\mathcal{N}(C)$ is $\{e_1,e_2,e_4,e_5\}$. The solution of the linearized equation can be written as $x(t)=\sum_{i=1}^6 c_i \exp(\lambda_it) v_i$ for constants $c_i$ that are determined by the initial condition: $x(0)=\sum_{i=1}^6 c_i v_i$. But, the eigenvectors $v_5,v_6$ do not satisfy the constraint given by [(\[eqn:con\])]{}, since they do not lie in $\mathcal{N}(C)$. Therefore, the constants $c_5,c_6$ are zero for initial conditions that are compatible with [(\[eqn:con\])]{}. We have $\mathrm{Re}[\lambda_i]<0$ for $1\leq i\leq 4$. Therefore, the equilibrium $(q,\omega)=(e_3,0)$ is asymptotically stable.

### Inverted Equilibrium

The eigenvalues $\lambda_i$, and the eigenvectors $v_i$ of the matrix $A$ at the inverted equilibrium $(-e_3,0)$ are given by $$\begin{gathered}
\lambda_{1,2}=-(\sqrt{5}+1)/2,\lambda_{3,4}=(\sqrt{5}-1)/2,\lambda_5=0,\lambda_6=-1,\nonumber\\
v_{1}= e_1 -(\sqrt{5}+1)e_4/2,\, v_2=e_2-(\sqrt{5}+1)e_5/2,\label{eqn:v1v2}\\
v_3=(\sqrt{5}+1)e_1/2 +e_4,\, v_4=(\sqrt{5}+1)e_2/2+e_5,\nonumber\\
v_5=e_3,\, v_6=e_6.\nonumber\end{gathered}$$ The basis of $\mathcal{N}(C)$ is $\{e_1,e_2,e_4,e_5\}$. Hence, the eigenvectors $v_5,v_6$ do not lie in $\mathcal{N}(C)$. Therefore, the solution can be written as $x(t)=\sum_{i=1}^4 c_i \exp(\lambda_it) v_i$ for constants $c_i$ that are determined by the initial condition.

We have $\mathrm{Re}[\lambda_{1,2}]<0$, and $\mathrm{Re}[\lambda_{3,4}]>0$. Therefore, the inverted equilibrium $(q,\omega)=(-e_3,0)$ is a hyperbolic equilibrium, and in particular, a saddle point.

Stable Manifold for the Inverted Equilibrium {#sec:SM}
--------------------------------------------

### Stable Manifold

The saddle point $(-e_3,0)$ has a stable manifold $W^s$, which is defined to be $$\begin{aligned}
W^s(-e_3,0) &= \{ (q,\omega)\in {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2\,|\, \lim_{t\rightarrow\infty} \mathcal{F}^t(q,\omega) = (-e_3,0)\},\end{aligned}$$ where $\mathcal{F}^t:(q(0),\omega(0))\rightarrow(q(t),\omega(t))$ denotes the flow map along the solution of [(\[eqn:dotw\])]{}, [(\[eqn:dotq\])]{}. The existence of $W^s(-e_3,0)$ has nontrivial effects on the overall dynamics of the controlled system. Trajectories in $W^s(-e_3,0)$ converge to the antipodal point of the desired equilibrium $(e_3,0)$, and it takes a long time period for any trajectory near $W^s(-e_3,0)$ to asymptotically converge to the desired equilibrium $(e_3,0)$.

According to the stable and unstable manifold theorem [@Kuz98], a local stable manifold $W^s_{loc}(-e_3,0)$ exists in the neighborhood of $(-e_3,0)$, and it is tangent to the stable eigenspace $E^s(-e_3,0)$ spanned by the eigenvectors $v_1$ and $v_2$ of the stable eigenvalues $\lambda_{1,2}$. The (global) stable manifold can be written as $$\begin{aligned}
W^s(-e_3,0) & = \bigcup_{t>0} \mathcal{F}^{-t} ( W^s_{loc}(-e_3,0)),\label{eqn:Ws}\end{aligned}$$ which states that the stable manifold $W^s$ can be obtained by globalizing the local stable manifold $W^s_{loc}$ by the backward flow map.

This yields a method to compute $W^s(-e_3,0)$ [@KraOsiIJBC05]. We choose a small ball $B_\delta\subset W^s_{loc}(-e_3,0)$ with a radius $\delta$ around $(-e_3,0)$, and we grow the manifold $W^s(-e_3,0)$ by evolving $B_\delta$ under the flow $\mathcal{F}^{-t}$. More explicitly, the stable manifold can be parameterized by $t$ as follows: $$\begin{aligned}
W^s(-e_3,0) =\{ \mathcal{F}^{-t} (B_\delta)\}_{t>0}.\label{eqn:Wc}\end{aligned}$$

We construct a ball in the stable eigenspace of $(-e_3,0)$ with sufficiently small radius $\delta$, i.e. $B_\delta\subset E^s_{loc}(-e_3,0)$. From the stable eigenvectors $v_1,v_2$ at [(\[eqn:v1v2\])]{}, $E^s_{loc}(-e_3,0)$ can be written as $$\begin{aligned}
E^s_{loc}& (-e_3,0) = \{  (q,\omega)\in {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2\,|\,
  q=\exp(\alpha_1\hat e_1+\alpha_2\hat e_2)(-e_3),\nonumber\\
& \omega=-\hat q^2(-(\sqrt{5}+1)/2)(\alpha_1 e_1+\alpha_2e_2)\text{ for $\alpha_1,\alpha_2\in{\ensuremath{\mathbb{R}}}$}\},\label{eqn:Esloc}\end{aligned}$$ where $-\hat q^2$ in the expression for $\omega$ corresponds to the orthogonal projection onto the plane normal to $q$, as required due to the constraint $q\cdot\omega=0$.

We define a distance on ${\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2$ as follows: $$\begin{aligned}
d_{{\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2} ((q_1,\omega_1),(q_2,\omega_2)) = \sqrt{\Psi(q_1,q_2)} + \|\omega_1-\omega_2\|.\label{eqn:dis}\end{aligned}$$ For $\delta>0$, the subset $B_\delta$ of $E^s_{loc}(-e_3,0)$ is parameterized by $\theta\in{\ensuremath{\mathsf{S}}}^1$ as $$\begin{aligned}
B_\delta & = \{  (q,\omega)\in {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2\,|\,
  q=\exp(\alpha_1\hat e_1+\alpha_2\hat e_2)(-e_3),\nonumber\\
& \omega=-\hat q^2(-(\sqrt{5}+1)/2)(\alpha_1 e_1+\alpha_2e_2),\text{ where}\nonumber\\
& \text{$\alpha_1=\frac{\delta}{1/\sqrt{2}+(\sqrt{5}+1)/2}\cos\theta,$\;}\nonumber\\
& \text{$\alpha_2=\frac{\delta}{1/\sqrt{2}+(\sqrt{5}+1)/2}\sin\theta$,} \text{ for $\theta\in{\ensuremath{\mathsf{S}}}^1$} \}.\label{eqn:Bdelta}\end{aligned}$$ The given choice of the constants $\alpha_1,\alpha_2$ guarantees that any point in $B_\delta$ has a distance $\delta$ to $(-e_3,0)$ according to the distance metric [(\[eqn:dis\])]{}.

### Variational Integrators

The parameterization of the stable manifold $W_s$ in [(\[eqn:Wc\])]{} requires the computation of the backward flow map $\mathcal{F}^{-t}$. However, general purpose numerical integrators may not preserve the structure of the two-sphere or the underlying dynamic characteristics, such as energy dissipation rate, accurately, and they may yield qualitatively incorrect numerical results in simulating a complex trajectory over a long-time period [@HaiLub00].

Geometric numerical integration is concerned with developing numerical integrators that preserve geometric features of a system, such as invariants, symmetry, and reversibility. In particular, variational integrators are geometric numerical integrators for Lagrangian or Hamiltonian systems, constructed according to Hamilton’s principle. They have desirable computational properties of preserving symplecticity and momentum maps, and they exhibit good energy behavior [@MarWesAN01]. A variational integrator is developed for Lagrangian or Hamiltonian systems evolving on the two-sphere in [@LeeLeoIJNME08]. It preserves both the underlying symplectic properties and the structures of the two-sphere concurrently.

A variational integrator on ${\ensuremath{\mathsf{S}}}^2$ for the controlled dynamics of a spherical pendulum can be written in a backward-time integration form as follows: $$\begin{aligned}
    q_{{k}} & = -{\ensuremath{\left( h\omega_{k+1} - \frac{h^2}{2ml^2} M_{k+1} \right)}}\times q_{k+1}\nonumber\\
    &\quad + {\ensuremath{\left( 1-{\ensuremath{\left\| h\omega_{k+1} - \frac{h^2}{2ml^2} M_{k+1} \right\|}}^2 \right)}}^{1/2} q_{k+1},\label{eqn:qk}\\
    \omega_{{k}} & = \omega_{k+1} - \frac{h}{2ml^2} M_k -    \frac{h}{2ml^2} M_{k+1},\label{eqn:wk}\end{aligned}$$ where the constant $h>0$ is time step, the subscript $k$ denotes the value of a variable at the time $t_k=kh$, and $M_k = ml^2(-k_\omega \omega_k - k_q q_d\times q_k)$. For given $(q_{k+1},\omega_{k+1})$, we first compute $M_{k+1}$. Then, $q_k$ is obtained by [(\[eqn:qk\])]{}, followed by $M_k$, and $\omega_k$ is computed by [(\[eqn:wk\])]{}. This yields an explicit, discrete inverse flow map $\mathcal{F}_d^{-h}((q_{k+1},\omega_{k+1}))=(q_{k},\omega_{k})$.

### Visualization

We choose 100 points on the surface of $B_\delta$ with $\delta=10^{-6}$, and each point is integrated backward using [(\[eqn:qk\])]{}, [(\[eqn:wk\])]{} with timestep $h=0.002$. The resulting trajectories are illustrated in [Fig. \[fig:Ws\]]{} for several values of $t$. Each colored curve on the sphere represents a trajectory on ${\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2$, since at any point $q$ on the curve, the direction of $\dot q=\omega\times q$ is tangent to the curve at $q$, and the magnitude of $\dot q$ is indirectly represented by color shading.

We observe the following characteristics of the stable manifold $W_s(-e_3,0)$ of the inverted equilibrium:

-   The boundary of the stable manifold $W_s(-e_3,0)\subset {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{S}}}^2$ parameterized by $t$ is circular when projected onto ${\ensuremath{\mathsf{S}}}^2$.

-   Each trajectory in $W_s(-e_3,0)$ is on a great circle, when projected onto ${\ensuremath{\mathsf{S}}}^2$. According to the closed loop dynamics [(\[eqn:dotw\])]{}, and the given initial condition at the surface of $B_\delta$, the direction of $\dot\omega$ is always parallel to $\omega$. Therefore, the direction of $\omega$ is fixed, and the resulting trajectory of $q$ is on a great circle. This also corresponds to the fact that the eigenvalue $\lambda_1$ for the first mode representing the rotations about the first axis is equal to the eigenvalue $\lambda_2$ for the second mode representing the rotations about the second axis at [(\[eqn:v1v2\])]{}, i.e. the convergence rates of these two rotations are identical.

-   The angular velocity decreases to zero as the direction of the pendulum $q$ converges to $-e_3$.

-   The stable manifold $W_s(-e_3,0)$ may cover ${\ensuremath{\mathsf{S}}}^2$ multiple times if $t$ is sufficiently large, as illustrated at [Fig. \[fig:Ws9\]]{}. Therefore, at any point $q\in{\ensuremath{\mathsf{S}}}^2$, we can choose $\omega$ such that $(q,\omega)$ lies in the stable manifold $W^s(-e_3,0)$ (the corresponding value of $\omega$ is not unique, since if it is sufficiently large, $q$ can traverse the sphere several times before converging to $-e_3$). This is similar to *kicking* a damped spherical pendulum carefully such that it converges to the inverted equilibrium.

3D Pendulum
===========

A 3D pendulum is a rigid body supported by a frictionless pivot acting under a gravitational potential. This is a generalization of a planar pendulum or a spherical pendulum, as it has three rotational degrees of freedom. It has been shown that a 3D pendulum may exhibit irregular maneuvers [@ChaLeeJNS11].

We choose a reference frame, and a body-fixed frame. The origin of the body-fixed frame is located at the pivot point. The attitude of a 3D pendulum is the orientation of the body-fixed frame with respect to the reference frame, and it is described by a rotation matrix representing the linear transformation from the body-fixed frame to the reference frame. The configuration manifold of a 3D pendulum is the special orthogonal group, ${\ensuremath{\mathsf{SO(3)}}}=\{R\in{\ensuremath{\mathbb{R}}}^{3\times 3}\,|\, R^T R=I,\mathrm{det}[R]=1\}$.

The equations of motion for a 3D pendulum are given by $$\begin{gathered}
J\dot\Omega + \Omega\times J\Omega = mg \rho\times R^T e_3 + u,\label{eqn:Wdot2}\\
\dot R = R\hat\Omega,\label{eqn:Rdot}\end{gathered}$$ where the matrix $J\in{\ensuremath{\mathbb{R}}}^{3\times 3}$ is the inertia matrix of the pendulum about the pivot, and $\rho\in{\ensuremath{\mathbb{R}}}^3$ is the vector from the pivot to the center of mass of the pendulum represented in the body-fixed frame. The control moment at the pivot is denoted by $u\in{\ensuremath{\mathbb{R}}}^3$.

Control System
--------------

Several control systems have been developed on ${\ensuremath{\mathsf{SO(3)}}}$ [@BulLew05; @ChaMcCITAC09; @LeePACC11]. Here, we summarize a control system to stabilize a 3D pendulum to a fixed desired attitude $R_d\in{\ensuremath{\mathsf{SO(3)}}}$. Consider an attitude error function given by $$\begin{aligned}
\Psi(R,R_d)=\frac{1}{2}{\mbox{tr}\ensuremath{\negthickspace{\ensuremath{\left[ (I-R_d^TR)G \right]}}}},\end{aligned}$$ for a diagonal matrix $G=\mathrm{diag}[g_1,g_2,g_3]\in{\ensuremath{\mathbb{R}}}^{3\times 3}$ with $g_1,g_2,g_3>0$. The derivative of this attitude error function with respect to $R$ along the direction of $\delta R= R\hat\eta$ for $\eta\in{\ensuremath{\mathbb{R}}}^3$ is given by $$\begin{aligned}
{\ensuremath{\mathbf{D}}}_R & \Psi(R,R_d)\cdot\delta R = -\frac{1}{2}{\mbox{tr}\ensuremath{\negthickspace{\ensuremath{\left[ R_d^TR\hat\eta G \right]}}}} \\
& = \frac{1}{2}(GR_d^TR -R^TR_d G)^\vee \cdot\eta\equiv e_R\cdot \eta,\end{aligned}$$ where we use the property that $\mathrm{tr}[\hat x A]=-x\cdot(A-A^T)^\vee$ for any $x\in{\ensuremath{\mathbb{R}}}^3,A\in{\ensuremath{\mathbb{R}}}^{3\times 3}$. The *vee map*, $\vee:{\ensuremath{\mathfrak{so}(3)}}\rightarrow{\ensuremath{\mathbb{R}}}^3$, denotes the inverse of the hat map. An attitude error vector is defined as $e_R = \frac{1}{2}(GR_d^TR -R^TR_d G)\in{\ensuremath{\mathbb{R}}}^3$. For positive constants $k_\Omega,k_R$, we choose the following control input: $$\begin{aligned}
u = -k_R e_R -k_\Omega \Omega -mg \rho\times R^T e_3.\end{aligned}$$ The corresponding closed loop dynamics are given by $$\begin{gathered}
J\dot\Omega =- \Omega\times J\Omega  -k_R e_R -k_\Omega \Omega,\label{eqn:Wdot2}\\
\dot R = R\hat\Omega.\label{eqn:Rdot}\end{gathered}$$ This system has four equilibria: in addition to the desired equilibrium $(R_d,0)$, there exist three other equilibria at $(R_d\exp (\pi\hat e_i,0),0)$ for $i\in\{1,2,3\}$, which correspond to the rotation of the desired attitude by $180^\circ$ about each body-fixed axis.

The existence of additional, undesirable equilibria is due to the nonlinear topological structure of ${\ensuremath{\mathsf{SO(3)}}}$, and it cannot be avoided by constructing a different control system (as long as it is continuous). It has been shown that it is not possible to design a continuous feedback control stabilizing an attitude globally on [$\mathsf{SO(3)}$]{} [@BhaBerSCL00; @KodPICDC98].

The stability of the desired equilibrium can be studied by using the following Lyapunov function, $$\begin{aligned}
\mathcal{V} =\frac{1}{2}\Omega\cdot J\Omega + k_R \Psi(R,R_d).\end{aligned}$$ In this paper, we analyze the stability of each equilibrium by linearizing the closed loop dynamics to study the equilibrium structures more explicitly.

Linearization
-------------

A variation in ${\ensuremath{\mathsf{SO(3)}}}$ can be expressed as  [@LeeLeoPICCA05]: $$\begin{aligned}
R^\epsilon=R\exp(\epsilon\hat\eta),\quad \Omega^\epsilon=\Omega +\epsilon\delta\Omega,\label{eqn:delRdelW}\end{aligned}$$ for $\eta,\delta\Omega\in{\ensuremath{\mathbb{R}}}^3$. The corresponding infinitesimal variation of $R$ is given by $\delta R = R\hat\eta$. Substituting this into [(\[eqn:Rdot\])]{}, $$\begin{aligned}
R\hat\Omega\hat\eta + R\hat{\dot\eta}= R\hat\eta\hat\Omega + R\delta\hat\Omega.\end{aligned}$$ Using the property $\hat x \hat y -\hat y\hat x=\widehat{x\times y}$ for any $x,y,\in{\ensuremath{\mathbb{R}}}^3$, this can be rewritten as $$\begin{aligned}
\dot\eta = \delta\Omega -\hat\Omega\eta.\label{eqn:etadot}\end{aligned}$$ Similarly, by substituting [(\[eqn:delRdelW\])]{} into [(\[eqn:Wdot2\])]{}, we obtain $$\begin{aligned}
J\delta\dot\Omega & = -\delta\Omega\times J\Omega -\Omega\times J\delta\Omega\nonumber\\
&\quad -\frac{1}{2}k_R (GR_d^T R\hat\eta +\hat\eta R^T R_d G) -k_\Omega\delta\Omega,\nonumber\\
& = (\widehat{J\Omega}-\hat\Omega J -k_\Omega I) \delta\Omega -\frac{1}{2}k_RH\eta,\label{eqn:delWdot}\end{aligned}$$ where $H=\mathrm{tr}[R^T R_d G]I-R^T R_d G\in{\ensuremath{\mathbb{R}}}^{3\times 3}$, and we used the property, $\hat x A + A^T\hat x = \mathrm{tr}[A]I-A$ for any $x\in{\ensuremath{\mathbb{R}}}^3, A\in{\ensuremath{\mathbb{R}}}^{3\times 3}$. Equations [(\[eqn:etadot\])]{},[(\[eqn:delWdot\])]{} can be written in matrix form as $$\begin{aligned}
\dot x & = \begin{bmatrix}\dot\eta \\ \delta\dot\Omega \end{bmatrix}
=\begin{bmatrix}
-\hat\Omega & I\\
-\frac{1}{2}k_R J^{-1}H & J^{-1}(\widehat{J\Omega}-\hat\Omega J -k_\Omega I)
\end{bmatrix}
\begin{bmatrix}\eta \\ \delta\Omega \end{bmatrix}\nonumber\\
& = Ax.\label{eqn:xdotSO3}\end{aligned}$$ This corresponds to the linearized equation of motion of [(\[eqn:Wdot2\])]{}, [(\[eqn:Rdot\])]{}.

Equilibrium Solutions
---------------------

We choose the desired attitude as $R_d=I$. In addition to the desired equilibrium $(I,0)$, there are three additional equilibria, namely $(\exp(\pi\hat e_1),0)$, $(\exp(\pi\hat e_2),0)$, $(\exp(\pi\hat e_3),0)$. We study the eigen-structure of each equilibrium using the linearized equation [(\[eqn:xdotSO3\])]{}. We assume that $$\begin{aligned}
J=\mathrm{diag}[3,2,1]\,\mathrm{kgm^2},\; G=\mathrm{diag}[0.9,1,1.1],\; k_R=k_\Omega=1.\end{aligned}$$

### Equilibrium $(I,0)$

The eigenvalues of the matrix $A$ at the desired equilibrium $(I,0)$ are given by $$\begin{gathered}
\lambda_{1,2}=-0.1667\pm0.5676i,\\
\lambda_{3,4}=-0.25\pm 0.6614i,\\
\lambda_{5,6}=-0.5\pm 0.8367i.\end{gathered}$$ This equilibrium is an asymptotically stable focus.

### Equilibrium $(\exp(\pi\hat e_1),0)$

At this equilibrium, the eigenvalues and the eigenvectors of $A$ are given by $$\begin{gathered}
\lambda_1=-0.7813,\quad v_1= e_1-0.7813 e_4,\nonumber\\
\lambda_2=-0.5854,\quad v_2=e_2-0.5854e_5,\nonumber\\
\lambda_3=-1.0477,\quad v_3=e_3-1.0477e_6,\label{eqn:v3_SO31}\\
\lambda_4=0.4480,\quad v_4=e_1+0.4480e_4,\nonumber\\
\lambda_5=0.0854,\quad v_5=e_2+0.0854e_5,\nonumber\\
\lambda_6=0.0477,\quad v_6=e_3+0.0477e_6.\nonumber\end{gathered}$$ Therefore, this equilibrium is a saddle point, where three modes are stable, and three modes are unstable.

### Equilibrium $(\exp(\pi\hat e_2),0)$

At this equilibrium, the eigenvalues and the eigenvectors of $A$ are given by $$\begin{gathered}
\lambda_1=-0.3775,\quad v_1= e_1-0.3775 e_4,\nonumber\\
\lambda_2=-1,\quad v_2=e_2-e_5,\label{eqn:v2_SO32}\\
\lambda_3=-0.9472,\quad v_3=e_3-0.9472e_6,\nonumber\\
\lambda_4=-0.0528,\quad v_4=e_3-0.0528e_6,\nonumber\\
\lambda_5=0.0442,\quad v_5=e_1+0.0442e_4,\nonumber\\
\lambda_6=0.5,\quad v_6=e_2+5e_5.\nonumber\end{gathered}$$ Therefore, this equilibrium is a saddle point, where four modes are stable, and two modes are unstable.

### Equilibrium $(\exp(\pi\hat e_3),0)$

At this equilibrium, the eigenvalues and the eigenvectors of $A$ are given by $$\begin{gathered}
\lambda_1=-0.0613,\quad v_1= e_1-0.0613 e_4,\nonumber\\
\lambda_2=-0.2721,\quad v_2= e_1-0.2721 e_4,\nonumber\\
\lambda_3=-0.1382,\quad v_3= e_2-0.1382 e_5,\nonumber\\
\lambda_4=-0.3618,\quad v_4= e_2-0.3618 e_5,\nonumber\\
\lambda_5=-1.5954,\quad v_5= e_3-1.5954 e_6,\label{eqn:v5_SO33}\\
\lambda_6= 0.5954,\quad v_6= e_2+0.5954 e_6.\nonumber\end{gathered}$$ Therefore, this equilibrium is a saddle point, where five modes are stable, and one mode is unstable.

Stable Manifolds for the Saddle Points
--------------------------------------

The eigen-structure analysis shows that there exist multi-dimensional stable manifolds for each saddle point. They have zero measure as the dimension of stable manifold is less than the dimension of ${\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{SO(3)}}}$. But, the existence of these stable manifolds may have nontrivial effects on the attitude dynamics.

We numerically characterize these stable manifolds using backward time integration, as discussed in Section \[sec:SM\].

The stable eigenspace for each saddle point can be written as $$\begin{aligned}
&E^s_{loc} (\exp(\pi\hat e_1),0) = \{ (R,\Omega)\in {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{SO(3)}}}\,|\,\nonumber\\
&\; R=\exp(\pi\hat e_1)\exp(\alpha_1\hat e_1+\alpha_2\hat e_2+\alpha_3\hat e_3),\\
&\; \Omega=-0.7813\alpha_1e_1-0.5854\alpha_2e_2-1.0477\alpha_3e_3\text{ for $\alpha_i\in{\ensuremath{\mathbb{R}}}$}\},\\
&E^s_{loc} (\exp(\pi\hat e_2),0) = \{ (R,\Omega)\in {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{SO(3)}}}\,|\,\nonumber\\
&\; R=\exp(\pi\hat e_2)\exp(\alpha_1\hat e_1+\alpha_2\hat e_2+(\alpha_3+\alpha_4)\hat e_3),\\
&\; \Omega=-0.37\alpha_1e_1-\alpha_2e_2-(0.94\alpha_3+0.05\alpha_4)e_3\text{ for $\alpha_i\in{\ensuremath{\mathbb{R}}}$}\},\\
&E^s_{loc} (\exp(\pi\hat e_3),0) = \{ (R,\Omega)\in {\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{SO(3)}}}\,|\,\nonumber\\
&\; R=\exp(\pi\hat e_3)\exp((\alpha_1+\alpha_2)\hat e_1+(\alpha_3+\alpha_4)\hat e_2+\alpha_5\hat e_3),\\
&\; \Omega=-(0.06\alpha_1+0.27\alpha_2)e_1-(0.13\alpha_3+0.36\alpha_4)e_2\\
&\;\quad -1.59\alpha_5e_3\text{ for $\alpha_i\in{\ensuremath{\mathbb{R}}}$}\},\end{aligned}$$

We define a distance on ${\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{SO(3)}}}$ as follows: $$\begin{aligned}
d_{{\ensuremath{\mathsf{T}}}{\ensuremath{\mathsf{SO(3)}}}} ((R_1,\Omega_1),(R_2,\Omega_2)) = \sqrt{\Psi(R_1,R_2)} + \|\Omega_1-\Omega_2\|.\end{aligned}$$

A variational integrator for the attitude dynamics of a rigid body on ${\ensuremath{\mathsf{SO(3)}}}$ is developed in [@LeeLeoPICCA05; @LeeLeoCMDA07]. It can be rewritten in a backward integration form as follows: $$\begin{gathered}
h (\Pi_{k+1}-\frac{h}{2}M_{k+1})^\wedge = J_dF_k - F_k^TJ_d,\label{eqn:Fk}\\
R_k = R_{k+1}F_k^T,\label{eqn:Rk}\\
\Pi_k = F_k\Pi_{k+1}-\frac{h}{2}F_kM_{k+1}-\frac{h}{2} M_k,\label{eqn:Pik}\end{gathered}$$ where $M_{k}=u_k+mg\rho\times R^T e_3\in{\ensuremath{\mathbb{R}}}^3$ is the external moment, $\Pi_k=J\Omega_k\in{\ensuremath{\mathbb{R}}}^3$ is the angular momentum. The matrix $J_d\in{\ensuremath{\mathbb{R}}}^{3\times 3}$ denotes a non-standard inertia matrix given by $J_d = \frac{1}{2}\mathrm{tr}[J]I-J$, and the rotation matrix $F_k\in{\ensuremath{\mathsf{SO(3)}}}$ represent the relative attitude update between two integration time steps. For given $(R_{k+1},\Pi_{k+1})$, we first compute $M_{k+1}$, and solve [(\[eqn:Fk\])]{} for $F_k$. Then, $R_k$ is obtained by [(\[eqn:Rk\])]{}, and $\Pi_k$ is computed by [(\[eqn:Pik\])]{}. This yields a discrete inverse flow map, $\mathcal{F}^{-h}_d(R_{k+1},\Pi_{k+1})\rightarrow(R_{k},\Pi_{k})$.

### Visualization of $W_s(\exp(\pi\hat e_1),0)$

In [@LeeLeoPICDC08], a method to visualize a function or a trajectory on ${\ensuremath{\mathsf{SO(3)}}}$ is proposed. Each column of a rotation matrix represents the direction of a body-fixed axis, and it evolves on ${\ensuremath{\mathsf{S}}}^2$. Therefore, a trajectory on ${\ensuremath{\mathsf{SO(3)}}}$ can be visualized by three curves on a sphere, representing the trajectory of three columns of a rotation matrix. The direction of the angular velocity should be chosen such that the corresponding time-derivative of the rotation matrix is tangent to the curve, and the magnitude of angular velocity can be illustrated by color shading. An example of visualizing a rotation about a single axis is illustrated in [Fig. \[fig:visSO3\_demo\]]{}.

We choose 112 points on the surface of $B_\delta\subset E^s_{loc} (\exp(\pi\hat e_1),0)$ with $\delta=10^{-6}$, and each point is integrated backward using [(\[eqn:qk\])]{}, [(\[eqn:wk\])]{} with timestep $h=0.002$. The resulting trajectories are illustrated in [Fig. \[fig:SM1\]]{} for several values of $t$.

In each figure, three body-fixed axes of the desired attitude $R_d=[e_1,e_2,e_3]$, and three body-fixed axes of the additional equilibrium attitude $\exp(\pi\hat e_1)=[e_1,-e_2,-e_3]$ are shown. From these computational results, we observe the following characteristics on the stable manifold $W_s(\exp(\pi\hat e_1),0)$:

-   When $t\leq 15$, the trajectories in $W_s(\exp(\pi\hat e_1),0)$ are close to rotations about the third body-fixed axis $e_3$ to $\exp(\pi\hat e_1)$. This is consistent with the linearized dynamics, where the eigenvalue of the third mode, corresponding to the rotations about $e_3$, has the fastest convergence rate, as seen in [(\[eqn:v3\_SO31\])]{}.

-   When $t\geq 15$, the first mode representing the rotations about $e_1$ starts to appear, followed by the second mode representing the rotation about $e_2$. This corresponds to the fact that the first mode has a faster convergence rate than the second mode, i.e. $|\lambda_1|>|\lambda_2|$.

-   As $t$ is increased further, the third body-fixed axis leaves the neighborhood of $-e_3$, and it exhibit the following pattern:

    ![image](SO3RA_k9001_a){width="0.32\columnwidth"}

-   The stable manifold $W_s(\exp(\pi\hat e_1),0)$ covers a certain part of ${\ensuremath{\mathsf{SO(3)}}}$, when projected on to it. So, when an initial attitude is chosen such that its third body-fixed axis is sufficiently close to $-e_3$, there possibly exist multiple initial angular velocities such that the corresponding solution converges to $\exp(\pi\hat e_1)$ instead of the desired attitude $R_d=I$.

### Visualization of $W_s(\exp(\pi\hat e_2),0)$

We choose 544 points on the surface of $B_\delta\subset E^s_{loc} (\exp(\pi\hat e_2),0)$ with $\delta=10^{-6}$, and each point is integrated backward using [(\[eqn:qk\])]{}, [(\[eqn:wk\])]{} with timestep $h=0.002$. The resulting trajectories are illustrated in [Fig. \[fig:SM2\]]{} for several values of $t$.

In each figure, three body-fixed axes of the desired attitude $R_d=[e_1,e_2,e_3]$, and three body-fixed axes of the additional equilibrium attitude $\exp(\pi\hat e_2)=[-e_1,e_2,-e_3]$ are shown. From these computational results, we observe the following characteristics on the stable manifold $W_s(\exp(\pi\hat e_2),0)$:

-   When $t\leq 12$, the trajectories in $W_s(\exp(\pi\hat e_2),0)$ is close to the rotations about the second body-fixed axis $e_2$. As $t$ increases, rotations about $e_3$ starts to appear. This corresponds to the linearized dynamics where the second mode representing rotations about $e_2$ has the fastest convergence rate, followed by the third mode at [(\[eqn:v2\_SO32\])]{}.

-   As $t$ is increased further, nonlinear modes become dominant. The trajectories in $W_s(\exp(\pi\hat e_2),0)$ almost cover [$\mathsf{SO(3)}$]{}. This suggests that for any initial attitude, we can choose several initial angular velocities such that the corresponding solutions converges to $\exp(\pi\hat e_2)$.

### Visualization of $W_s(\exp(\pi\hat e_3),0)$

Similarly, we choose 976 points on the surface of $B_\delta\subset E^s_{loc} (\exp(\pi\hat e_3),0)$ with $\delta=10^{-6}$, and each point is integrated backward using [(\[eqn:qk\])]{}, [(\[eqn:wk\])]{} with timestep $h=0.002$. The resulting trajectories are illustrated in [Fig. \[fig:SM3\]]{} for several values of $t$.

At each figure, three body-fixed axes of the desired attitude $R_d=[e_1,e_2,e_3]$, and three body-fixed axes of the additional equilibrium attitude $\exp(\pi\hat e_3)=[-e_1,-e_2,e_3]$ are shown. From these computational results, we observe the following characteristics on the stable manifold $W_s(\exp(\pi\hat e_3),0)$:

-   When $t\leq 8$, the trajectories in $W_s(\exp(\pi\hat e_3),0)$ are close to the rotations about the third body-fixed axis $e_3$. This corresponds to the linearized dynamics where the fifth mode representing rotations about $e_3$ has the fastest convergence rate given in [(\[eqn:v5\_SO33\])]{}.

-   The rotations about $e_3$ are still dominant, even as $t$ is increased further. For the given simulation times, all trajectories in $W_s(\exp(\pi\hat e_3),0)$ are close to rotations about $e_3$.

Conclusions
===========

Stable manifolds of saddle points that arise in the closed-loop dynamics of two pendulum models are characterized numerically, and several properties are observed. Although the analytical and computational results have been presented for a spherical pendulum and a 3D pendulum, the methods presented naturally extend to any closed loop attitude control system with configurations in either ${\ensuremath{\mathsf{S}}}^2$ or ${\ensuremath{\mathsf{SO(3)}}}$.

[^1]: Taeyoung Lee, Mechanical and Aerospace Engineering, Florida Institute of Technology, Melbourne, FL 39201 [taeyoung@fit.edu]{}

[^2]: Melvin Leok, Mathematics, University of California at San Diego, La Jolla, CA 92093 [mleok@math.ucsd.edu]{}

[^3]: N. Harris McClamroch, Aerospace Engineering, University of Michigan, Ann Arbor, MI 48109 [nhm@umich.edu]{}

[^4]: ^$*$^This research has been supported in part by NSF under grants CMMI-1029551.

[^5]: ^$\dagger$^This research has been supported in part by NSF under grants DMS-0726263, DMS-1001521, DMS-1010687, and CMMI-1029445.
