# Formulation and Simulations of the Conserving Algorithm for Feedback Stabilization on Rigid Body Rotations

Authors: Yong-Ren Pu; Thomas A. Posbergh
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The problem of stabilization of rigid bodies has received a great deal of attention for many years. People have developed a variety of feedback control laws to meet their design requirements and have formulated various but mostly open loop numerical algorithms for the dynamics of the corresponding closed loop systems. Since the conserved quantities such as energy, momentum, and symmetry play an important role in the dynamics, we investigate the conserved quantities for the closed loop control systems which formally or asymptotically stabilize rigid body rotation and modify the open loop numerical algorithms so that they preserve these important properties. Using several examples, the authors first use the open loop algorithm to simulate the tumbling rigid body actions and then use the resulting closed loop one to stabilize them.


---

## Body

## 1. Introduction

The ideas of modern geometry may always find a role in the applications of the engineering sciences. Through geometry people can explore the structure of nonlinear problems in a natural and intuitive way. Feedback control of a rotating rigid body has become one of such research subjects. The model remains of fundamental interest from both a theoretical and applications point of view. From a theoretical point of view the nonlinear rigid body dynamics provides an important vehicle for the development of nonlinear control theory [1, 2]. At the same time the rigid body equations of motion appear in important applications ranging from spacecraft attitude control [3] to industrial robots [4].For many systems it is essential that the algorithms used to simulate the dynamics produce physically meaningful results. To achieve this end algorithms should be based on physically sound models and capture without approximation the important physical characteristics of the actual system. Such quantities include symmetry, energy, momentum, Poisson, and symplectic structures. Several algorithms which preserve such invariants have been studied in detail. Early work in this area includes that of [5] who looked at energy and momentum conserving algorithms, [6] who focused on the integration of Poisson structure, [7] who considered algorithms which preserved symplectic structure, and [8] who used an imposed Lie-Poisson structure to determine the rigid body trajectory in Euclidean space. The conservation of such quantities is crucial in meaningful numerical simulation of systems possessing such structure. For example, consider the dynamics of a spin stabilized satellite. Such a satellite may have a rotational period on the order of several rotations per minute for the mission duration on the order of years. Failure to conserve angular momentum when simulating such a system will lead to significant errors and inaccurate conclusions about system performance.In the area of stability analysis of uniform rigid body rotation the geometric framework has seen some development of methods. The energy-Casimir method was used in the context of Lie-Poisson systems to find stabilizing feedback control laws [9, 10]. A class of feedback control laws for the stabilization of uniform rigid body rotation was introduced in [11]. The energy and momentum conserving algorithms described in [12] were used as the basic methodology for the numerical validation of these stabilizing control laws. These algorithms were suitably modified to incorporate feedback controls based on the methodology introduced and developed in these papers.In this paper we will discuss the modifications which are necessary to apply the algorithm introduced in [12] to rigid body motion with feedback. While the modifications are straightforward they raise important issues with regard to numerical simulation of such systems. The results represent special cases of algorithms which conserve energy and momentum for general Hamiltonian systems [13, 14]. The remainder of the paper is then as follows. In the next section we outline the formulation of rigid body dynamics from a geometric point of view. We then show how the dynamics are modified when feedback is applied to the system. The following section then discusses the energy and momentum conserving algorithm. We then introduce and justify the modifications to the algorithm in the presence of feedback. We conclude this paper with several examples.


## 2. Rigid Body Dynamics

The unforced dynamics of a rigid body free to rotate about a fixed point possess two conserved quantities, the kinetic energy and the (spatially referenced) angular momentum. These quantities are constant with respect to unforced motion of the rigid body. If we allow this motion to be modified by state feedback, it is a remarkable fact that the closed loop systems also possess two quantities analogous to the energy and momentum of the original system [9].


### 2.1. The Open Loop Rigid Body Dynamics

The dynamics of the rigid body free to move about a fixed point benefit from a geometric description [15, 16]. The configuration of the rigid body at any time can be described by an element of the special orthogonal group, that is, the set of all orthogonal matrices with determinant equal to one. We denote this group as S
O
(
3
). This group is a Lie group.If we consider a motion of this rigid body this is represented by a continuous curve of elements ofS
O
(
3
) parameterized by time. The tangent space to any point Λ in S
O
(
3
) consists of the tangent vectors to all the possible curves passing through such a point; we denote it by  T
Λ
S
O
(
3
). Since S
O
(
3
) is a Lie group, it has an associated Lie algebra which we denote by s
o
(
3
). We can identify any vector in the tangent space  T
Λ
S
O
(
3
) with two elements in the Lie algebra of S
O
(
3
) by translating this element to the identity by left or right translation. These two elements in the Lie algebra correspond to the convected angular velocity  W and the spatial angular velocity  w. In addition to the tangent space at the current configuration we can define a cotangent space T
Λ
*
S
O
(
3
) through the duality pairing 〈
π
,
w
〉
=
(
1
/
2
)
t
r
w
T
π. As with an element of the tangent space, an element of the cotangent space can be translated to two elements in the dual to the Lie algebra s
o
*
(
3
). For the natural pairing these two elements correspond to the convected angular momentum Π and the spatial angular momentum π.The rigid body is a Hamiltonian system which has a phase space, a Hamiltonian, and a symplectic form. The phase space is denoted byP and defined as the cotangent bundle of S
O
(
3
)
(1)
P
=
T
*
S
O
(
3
)
=
⋃
A
T
A
*
S
O
(
3
)
.
Locally we coordinatize elements of this phase space as (
Λ
,
π
) with Λ a rotation matrix and π a three-element angular momentum vector. We define a Hamiltonian function H
(
Λ
,
π
) on this phase space. In spatial coordinates
(2)
H
(
Λ
,
π
)
=
1
2
π
·
I
-
1
π
.
Here I
=
Λ
J
Λ
T denotes the spatially referenced inertia matrix and J is the body fixed inertia. Hamiltonian’s equations are then
(3)
Λ
˙
=
w
^
Λ
,
π
˙
=
0
.
Application of the Legendre transformation yields
(4)
Λ
˙
=
w
^
Λ
,
I
w
˙
=
I
w
×
w
,
where w
∈
ℜ
3 lies in the space of Lie algebra with the standard isomorphism ^
:
ℜ
3
→
s
o
(
3
). Note that ℜ represents the set of real numbers. By standard arguments for Hamiltonian systems one then shows that the energy is conserved [2].For the rigid body the symmetry induced by the action ofS
O
(
3
) on the phase space provides an additional conserved quantity, the angular momentum. This symmetry is used to define a momentum map J
:
P
→
s
o
*
(
3
) corresponding to the invariance under rigid body rotations. It is easily shown that this corresponds to the classical notion of angular momentum; that is
(5)
J
=
π
.


### 2.2. The Closed Loop System

Application of state feedback can be used to stabilize the uniform rotation of the rigid body about an arbitrary axis in space. Such state feedback modifies the vector field corresponding to rigid body motion. In what follows we consider feedback of the velocity components of the rigid body. Equation (3) is now
(6)
Λ
˙
=
w
^
Λ
,
I
w
˙
=
I
w
×
w
+
u
,
where u is an external moment applied to the system.Consider a general form of velocity feedback for the rigid body given by(7)
u
=
a
w
×
w
+
b
×
w
,
where  a
=
Λ
A
Λ
T, b
=
Λ
B and A
∈
ℜ
3
×
3,  B
∈
ℜ
3 are constant matrices. Using (7) in (6) we obtain the closed loop dynamics
(8)
Λ
˙
=
w
^
Λ
,
I
w
˙
=
[
(
I
+
a
)
w
+
b
]
×
w
.Remarkably, with a more general definition of the phase space this system possesses conserved quantities analogous to the energy and momentum of the free rigid body. In what follows we will consider a special case of (8) by setting a
=
0 in (7). Such feedback arises, for example, in the case of a spacecraft carrying momentum wheels. The second equation of (8) then takes the form
(9)
I
w
˙
=
(
I
w
+
b
)
×
w
.The phase space of this system is the cotangent bundle ofS
O
(
3
). We define the comomentum vector as π
*
=
I
w
+
b. The vector b has the effect of “shifting” the angular momentum due to the rotation of the rigid body. It is straightforward to verify that π
*
=
0 along the trajectory. The Hamiltonian
(10)
H
(
Λ
,
π
*
)
=
1
2
(
π
*
-
b
)
·
I
-
1
(
π
*
-
b
)
has the property that H
˙
=
0 along the flow given by (8) from which we conclude that the Hamiltonian is constant.As before there exists a momentum mapJ which is conserved by the flow. This momentum map reflects the invariance of the dynamics with respect to rigid body rotations. It is shown that
(11)
J
=
I
w
+
b
,
and that J
˙
=
0 is then a direct consequence of Noether’s theorem [2].


## 2.1. The Open Loop Rigid Body Dynamics

The dynamics of the rigid body free to move about a fixed point benefit from a geometric description [15, 16]. The configuration of the rigid body at any time can be described by an element of the special orthogonal group, that is, the set of all orthogonal matrices with determinant equal to one. We denote this group as S
O
(
3
). This group is a Lie group.If we consider a motion of this rigid body this is represented by a continuous curve of elements ofS
O
(
3
) parameterized by time. The tangent space to any point Λ in S
O
(
3
) consists of the tangent vectors to all the possible curves passing through such a point; we denote it by  T
Λ
S
O
(
3
). Since S
O
(
3
) is a Lie group, it has an associated Lie algebra which we denote by s
o
(
3
). We can identify any vector in the tangent space  T
Λ
S
O
(
3
) with two elements in the Lie algebra of S
O
(
3
) by translating this element to the identity by left or right translation. These two elements in the Lie algebra correspond to the convected angular velocity  W and the spatial angular velocity  w. In addition to the tangent space at the current configuration we can define a cotangent space T
Λ
*
S
O
(
3
) through the duality pairing 〈
π
,
w
〉
=
(
1
/
2
)
t
r
w
T
π. As with an element of the tangent space, an element of the cotangent space can be translated to two elements in the dual to the Lie algebra s
o
*
(
3
). For the natural pairing these two elements correspond to the convected angular momentum Π and the spatial angular momentum π.The rigid body is a Hamiltonian system which has a phase space, a Hamiltonian, and a symplectic form. The phase space is denoted byP and defined as the cotangent bundle of S
O
(
3
)
(1)
P
=
T
*
S
O
(
3
)
=
⋃
A
T
A
*
S
O
(
3
)
.
Locally we coordinatize elements of this phase space as (
Λ
,
π
) with Λ a rotation matrix and π a three-element angular momentum vector. We define a Hamiltonian function H
(
Λ
,
π
) on this phase space. In spatial coordinates
(2)
H
(
Λ
,
π
)
=
1
2
π
·
I
-
1
π
.
Here I
=
Λ
J
Λ
T denotes the spatially referenced inertia matrix and J is the body fixed inertia. Hamiltonian’s equations are then
(3)
Λ
˙
=
w
^
Λ
,
π
˙
=
0
.
Application of the Legendre transformation yields
(4)
Λ
˙
=
w
^
Λ
,
I
w
˙
=
I
w
×
w
,
where w
∈
ℜ
3 lies in the space of Lie algebra with the standard isomorphism ^
:
ℜ
3
→
s
o
(
3
). Note that ℜ represents the set of real numbers. By standard arguments for Hamiltonian systems one then shows that the energy is conserved [2].For the rigid body the symmetry induced by the action ofS
O
(
3
) on the phase space provides an additional conserved quantity, the angular momentum. This symmetry is used to define a momentum map J
:
P
→
s
o
*
(
3
) corresponding to the invariance under rigid body rotations. It is easily shown that this corresponds to the classical notion of angular momentum; that is
(5)
J
=
π
.


## 2.2. The Closed Loop System

Application of state feedback can be used to stabilize the uniform rotation of the rigid body about an arbitrary axis in space. Such state feedback modifies the vector field corresponding to rigid body motion. In what follows we consider feedback of the velocity components of the rigid body. Equation (3) is now
(6)
Λ
˙
=
w
^
Λ
,
I
w
˙
=
I
w
×
w
+
u
,
where u is an external moment applied to the system.Consider a general form of velocity feedback for the rigid body given by(7)
u
=
a
w
×
w
+
b
×
w
,
where  a
=
Λ
A
Λ
T, b
=
Λ
B and A
∈
ℜ
3
×
3,  B
∈
ℜ
3 are constant matrices. Using (7) in (6) we obtain the closed loop dynamics
(8)
Λ
˙
=
w
^
Λ
,
I
w
˙
=
[
(
I
+
a
)
w
+
b
]
×
w
.Remarkably, with a more general definition of the phase space this system possesses conserved quantities analogous to the energy and momentum of the free rigid body. In what follows we will consider a special case of (8) by setting a
=
0 in (7). Such feedback arises, for example, in the case of a spacecraft carrying momentum wheels. The second equation of (8) then takes the form
(9)
I
w
˙
=
(
I
w
+
b
)
×
w
.The phase space of this system is the cotangent bundle ofS
O
(
3
). We define the comomentum vector as π
*
=
I
w
+
b. The vector b has the effect of “shifting” the angular momentum due to the rotation of the rigid body. It is straightforward to verify that π
*
=
0 along the trajectory. The Hamiltonian
(10)
H
(
Λ
,
π
*
)
=
1
2
(
π
*
-
b
)
·
I
-
1
(
π
*
-
b
)
has the property that H
˙
=
0 along the flow given by (8) from which we conclude that the Hamiltonian is constant.As before there exists a momentum mapJ which is conserved by the flow. This momentum map reflects the invariance of the dynamics with respect to rigid body rotations. It is shown that
(11)
J
=
I
w
+
b
,
and that J
˙
=
0 is then a direct consequence of Noether’s theorem [2].


## 3. Conserving Algorithms

Discretization of the closed loop dynamics (4) will not in general preserve the conserved quantities of energy and momentum. In fact, standard numerical methods such as Runge-Kutta or direct integration will in general remain valid for only short periods of time. For long term dynamics to retain the conserved quantities we must use algorithms which have been designed for this purpose.Reference [17] developed a Lie-Poisson Hamilton-Jacobi theory leading to a symplectic integrator that preserves exactly the Lie-Poisson structure but does not typically conserve the energy. Actually it showed that the energy conserving algorithms cannot be symplectic at the same time in the nonintegrable or reduced space. Another almost Poisson symplectic integrator introduced by [18] with mid-point rule conserves the first integrals containing only linear and quadratic terms, especially the energy, but only conserves the Poisson structure up to the second order. Reference [12] constructed energy-momentum algorithms by balancing the spatial angular momentum which also has the second order accuracy in the Poisson structure that is compared with the symplectic integrator in [13].In contrast with the previous algorithms for the open loop systems, we are interested in the closed loop systems based on the energy-momentum algorithm with a family of stabilizing feedback introduced in the work of [19]. The idea can be found in [15, 16] that for the Hamiltonian system (rigid body rotation in this case) with symmetries induced by the action of a Lie group (S
O
(
3
) here) on the phase space, there exist conserved energy and momentum map along the integral line of the dynamics. The reason we choose energy-momentum conserving algorithm is that the spatial angular momentum plays the most important role in the rigid body rotation and is dealt with explicitly in this algorithm.


### 3.1. Energy-Momentum Conserving Algorithm for Open Loop System

We first review the open loop energy-momentum conserving algorithm of [12] in this section and then discuss the modifications for the closed loop system in the following section.The equations of motion of a rigid body in spatial coordinate are expressed in (4). They can be interpreted in body coordinate with applied torque as
(12)
Λ
˙
=
Λ
W
^
,
J
W
˙
=
J
W
×
W
+
U
,
with initial conditions,
(13)
Λ
(
t
)
|
t
=
0
=
Λ
0
,
W
(
t
)
|
t
=
0
=
W
0
,
where Λ
∈
S
O
(
3
) is an element of the configuration space, W
∈
s
o
(
3
) is an element of the velocity field, U
=
Λ
T
u is the body fixed applied torque, and J is the body fixed moment of inertia.The update of the configuration here is the exponential map exp:s
o
(
3
)
→
S
O
(
3
) from the Lie algebra to the group, which has the closed form called Rodrigues formula
(14)
Λ
=
exp
⁡
⁡
(
Θ
^
)
=
1
+
sin
⁡
∥
Θ
∥
∥
Θ
∥
Θ
^
+
1
2
sin
2
⁡
(
(
1
/
2
)
∥
Θ
∥
)
(
(
1
/
2
)
∥
Θ
∥
)
2
Θ
^
2
,
where Θ
∈
ℜ
3 is the convected relative rotation vector lying in the linear vector space of Lie algebra.Therefore the open loop energy-momentum conserving algorithm for the initial value problem (12) is based on the Newmark algorithm [20] outlined as follows.For a typical time interval[
t
n
,
t
n
+
1
]
⊂
[
0
,
T
] with the known variables Λ
n, W
n, A
n, u
(
t
), and h
=
t
n
+
1
-
t
n where A
∈
ℜ
3 is the convected angular acceleration, we let A
n
+
1 to be zero initially.Step 1.
(15)
Θ
=
W
n
h
+
h
2
[
(
1
2
-
β
)
A
n
+
β
A
n
+
1
]
,
(16)
W
n
+
1
=
γ
β
h
Θ
-
W
n
+
(
2
-
γ
β
)
(
W
n
+
1
2
h
A
n
)
.Step 2.
(17)
Λ
n
+
1
=
Λ
n
exp
⁡
⁡
(
Θ
^
)
.Step 3.
(18)
R
∶
=
Λ
n
+
1
J
W
n
+
1
-
Λ
n
J
W
n
-
h
u
n
+
α
,
0
≤
α
≤
1
.Step 4.
If∥
R
∥  < Tolerance, then
(19)
A
n
+
1
=
W
n
+
1
-
W
n
γ
h
+
(
1
-
1
γ
)
A
n
.
Letn
←
n
+
1, go to Step 1.
ElseStep 5.
(20)
K
∶
=
Λ
n
+
1
[
γ
β
h
J
T
(
Θ
)
-
J
W
^
n
+
1
]
.Step 6.
(21)
Δ
Θ
=
K
-
1
R
,
Θ
←
Θ
+
Δ
Θ, go to Step 2.To accomplish the iteration by the Newton-Raphson method we need to linearize (18) about (
Λ
,
Θ
). Equation (18) is nonlinear in the finite rotation vector Θ
.By [21] the linearization needs the following operations of the directional derivative in the direction ΔΘ
(22)
D
exp
⁡
⁡
(
Θ
^
)
·
Δ
Θ
=
exp
⁡
⁡
(
Θ
^
)
Δ
Θ
,
D
Θ
·
Δ
Θ
=
T
(
Θ
)
Δ
Θ
,
where D
(
·
)
·
Δ
Θ is the directional derivative of (
·
) in the direction Δ
Θ and T
:
ℜ
3
→
s
o
(
3
) is a linear map having the expression
(23)
T
(
Θ
)
:
=
e
⊗
e
+
(
1
/
2
)
∥
Θ
∥
tan
⁡
(
(
1
/
2
)
∥
Θ
∥
)
(
1
-
e
⊗
e
)
+
1
2
Θ
^
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
e
=
Θ
∥
Θ
∥
.
Thus the tangent matrix K calculated by linearization is derived. We will use these equations again to derive new tangent matrix K
* for the closed loop systems later.It can be shown that for the choice ofγ
/
β
=
2 this algorithm conserves energy exactly and thus is unconditionally stable [12]. Moreover, if we choose γ
=
1 and β
=
1
/
2, the updated angular acceleration A
n
+
1 does not depend on the initial A
n
+
1 during each time interval and we may avoid the accumulated errors from A
n. Moreover, it can also be proved that if, in addition, we use the trapezoidal rule for the applied torque, that is, α
=
1
/
2, then the algorithm is second order accurate both in configuration and velocities.


### 3.2. New Result for Closed Loop System

When we close the loop with the linear stabilizing feedback ofU
=
B
×
W by energy-momentum method mentioned in [20], (12) become
(24)
Λ
˙
=
Λ
W
^
,
J
W
˙
=
(
J
W
+
B
)
×
W
,
which are the equations of motion for the closed loop system in body coordinate.As in Section2.2 we define the convected shifted angular momentum Π
*
=
J
W
+
B. Clearly the norm of  Π
* is constant along the flow. We can also define the Hamiltonian H
=
(
1
/
2
)
(
Π
*
-
B
)
·
J
-
1
(
Π
*
-
B
) in the dual Lie algebra, or as normally H
=
(
1
/
2
)
J
W
·
W in the Lie algebra, and both can be easily seen to be conserved.Since we have closed the system with the constant feedback gain, there is no need to trace the control input and we do not need the parameterα anymore in our algorithm. We only have to set the constant vector B beforehand.In the closed loop energy-momentum conserving algorithm we hope to have the same nice properties as those in the open loop one. We can achieve that by modifying the algorithm in the previous section. First, with the open loop algorithm at hand (24) can be discretized in the form
(25)
J
W
n
+
1
-
J
W
n
=
h
(
J
W
n
+
β
/
γ
+
B
)
×
W
n
+
β
/
γ
between any typical time interval [
t
n
,
t
n
+
1
]
⊂
[
0
,
T
] with the following:
(26)
Π
n
*
=
J
W
n
+
B
,
Π
n
+
1
*
=
J
W
n
+
1
+
B
,
π
n
*
=
Λ
n
(
J
W
n
+
B
)
,
π
n
+
1
*
=
Λ
n
+
1
(
J
W
n
+
1
+
B
)
.
It can be proved that the norm of the shifted angular momentum is conserved in the closed loop algorithm when γ
/
β
=
2.In addition to conservation of norm of the shifted momenta the direction of the shifted spatial angular momentum is also conserved if we enforce thatπ
n
+
1
+
b
n
+
1
=
π
n
+
b
n. To do this we just replace Step 3 in the previous section by
(27)
R
=
Λ
n
+
1
J
W
n
+
1
-
Λ
n
J
W
n
+
Λ
n
+
1
B
-
Λ
n
B
to modify the algorithm. Note that the linear feedback needs not to stabilize the rigid body rotation. In general, any feedback in the form U
=
B
×
W would conserve the shifted spatial angular momentum which is actually the momentum map for the Hamiltonian system. We will see this in some examples later.The conservation of energy depends on properly choosing the parametersα and β as well. From the kinematic equations (16), (17) and the fact that exp
⁡
⁡
(
Θ
^
)
Θ
=
Θ we update
(28)
w
n
+
1
=
w
n
+
Λ
n
W
n
+
1
-
Λ
n
+
1
W
n
+
(
2
-
γ
β
)
(
Λ
n
+
1
-
Λ
n
)
W
n
*
,
where W
n
*
=
W
n
+
(
1
/
2
)
h
A
n. Besides, (16) yields W
n
+
1
+
W
n
=
(
γ
/
β
h
)
Θ
+
(
2
-
γ
/
β
)
W
n
*. By eliminating W
n
+
1 in the previous equations and substituting into the balance equation
(29)
Λ
n
+
1
J
W
n
+
1
+
Λ
n
+
1
B
=
Λ
n
J
W
n
+
Λ
n
B
it turns out to be the following:
(30)
J
W
n
-
B
-
γ
β
h
J
Θ
+
exp
⁡
(
-
Θ
^
)
J
W
n
+
exp
⁡
(
-
Θ
^
)
B
=
(
2
-
γ
β
)
J
W
n
*
.It can also be proved that the energy is conserved in the closed loop algorithm whenγ
/
β
=
2. Again, if we choose the parameters γ
=
1 and β
=
1
/
2, we are free from the error propagation of angular acceleration A
n and the closed loop algorithm would be second order accurate. Moreover, it is unconditionally stable since the energy is conserved.The linearization process of (18) starts from the balance equation (29) with the orientation update (17)
(31)
Λ
n
[
J
W
n
exp
⁡
⁡
(
Θ
^
)
J
(
γ
β
h
Θ
-
W
n
)
+
B
-
exp
⁡
⁡
(
Θ
^
)
B
]
=
0
,
which is a nonlinear equation in Θ. Define
(32)
G
(
Θ
)
∶
=
Λ
n
[
J
W
n
exp
⁡
⁡
(
Θ
^
)
J
(
γ
β
h
Θ
-
W
n
)
k
k
k
k
k
k
+
B
-
exp
⁡
⁡
(
Θ
^
)
B
(
γ
β
h
Θ
-
W
n
)
]
;
then the linearization ofG(Θ) with respect to Θ becomes
(33)
L
[
G
(
Θ
)
]
=
G
(
Θ
0
)
+
D
G
(
Θ
)
·
Δ
Θ
,
where Θ
0 is the fixed point at which the linearization is exerted. By (22) and (23) with the definition (32), we can compute the linearization of G
(
Θ
) along Δ
Θ in body coordinate as follows:
(34)
D
G
(
Θ
)
·
Δ
Θ
=
-
Λ
n
+
1
[
γ
β
h
J
T
(
Θ
)
-
J
W
^
n
+
1
-
B
^
]
Δ
Θ
=
-
K
*
Δ
Θ
.Note that there is a term-
B
^ in the tangent matrix K
* that is different from the result (20) of open loop system. This means we need to put one more term in dealing with the linearization of the balance equation for the closed loop system.Based on the analysis above we treat the constant vectorB as a controller and put it inside the algorithm. Modify the linearization process in the previous section and replace the momenta Π and π by Π
* and π
* to construct the new closed loop algorithm.


## 3.1. Energy-Momentum Conserving Algorithm for Open Loop System

We first review the open loop energy-momentum conserving algorithm of [12] in this section and then discuss the modifications for the closed loop system in the following section.The equations of motion of a rigid body in spatial coordinate are expressed in (4). They can be interpreted in body coordinate with applied torque as
(12)
Λ
˙
=
Λ
W
^
,
J
W
˙
=
J
W
×
W
+
U
,
with initial conditions,
(13)
Λ
(
t
)
|
t
=
0
=
Λ
0
,
W
(
t
)
|
t
=
0
=
W
0
,
where Λ
∈
S
O
(
3
) is an element of the configuration space, W
∈
s
o
(
3
) is an element of the velocity field, U
=
Λ
T
u is the body fixed applied torque, and J is the body fixed moment of inertia.The update of the configuration here is the exponential map exp:s
o
(
3
)
→
S
O
(
3
) from the Lie algebra to the group, which has the closed form called Rodrigues formula
(14)
Λ
=
exp
⁡
⁡
(
Θ
^
)
=
1
+
sin
⁡
∥
Θ
∥
∥
Θ
∥
Θ
^
+
1
2
sin
2
⁡
(
(
1
/
2
)
∥
Θ
∥
)
(
(
1
/
2
)
∥
Θ
∥
)
2
Θ
^
2
,
where Θ
∈
ℜ
3 is the convected relative rotation vector lying in the linear vector space of Lie algebra.Therefore the open loop energy-momentum conserving algorithm for the initial value problem (12) is based on the Newmark algorithm [20] outlined as follows.For a typical time interval[
t
n
,
t
n
+
1
]
⊂
[
0
,
T
] with the known variables Λ
n, W
n, A
n, u
(
t
), and h
=
t
n
+
1
-
t
n where A
∈
ℜ
3 is the convected angular acceleration, we let A
n
+
1 to be zero initially.Step 1.
(15)
Θ
=
W
n
h
+
h
2
[
(
1
2
-
β
)
A
n
+
β
A
n
+
1
]
,
(16)
W
n
+
1
=
γ
β
h
Θ
-
W
n
+
(
2
-
γ
β
)
(
W
n
+
1
2
h
A
n
)
.Step 2.
(17)
Λ
n
+
1
=
Λ
n
exp
⁡
⁡
(
Θ
^
)
.Step 3.
(18)
R
∶
=
Λ
n
+
1
J
W
n
+
1
-
Λ
n
J
W
n
-
h
u
n
+
α
,
0
≤
α
≤
1
.Step 4.
If∥
R
∥  < Tolerance, then
(19)
A
n
+
1
=
W
n
+
1
-
W
n
γ
h
+
(
1
-
1
γ
)
A
n
.
Letn
←
n
+
1, go to Step 1.
ElseStep 5.
(20)
K
∶
=
Λ
n
+
1
[
γ
β
h
J
T
(
Θ
)
-
J
W
^
n
+
1
]
.Step 6.
(21)
Δ
Θ
=
K
-
1
R
,
Θ
←
Θ
+
Δ
Θ, go to Step 2.To accomplish the iteration by the Newton-Raphson method we need to linearize (18) about (
Λ
,
Θ
). Equation (18) is nonlinear in the finite rotation vector Θ
.By [21] the linearization needs the following operations of the directional derivative in the direction ΔΘ
(22)
D
exp
⁡
⁡
(
Θ
^
)
·
Δ
Θ
=
exp
⁡
⁡
(
Θ
^
)
Δ
Θ
,
D
Θ
·
Δ
Θ
=
T
(
Θ
)
Δ
Θ
,
where D
(
·
)
·
Δ
Θ is the directional derivative of (
·
) in the direction Δ
Θ and T
:
ℜ
3
→
s
o
(
3
) is a linear map having the expression
(23)
T
(
Θ
)
:
=
e
⊗
e
+
(
1
/
2
)
∥
Θ
∥
tan
⁡
(
(
1
/
2
)
∥
Θ
∥
)
(
1
-
e
⊗
e
)
+
1
2
Θ
^
,
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
k
e
=
Θ
∥
Θ
∥
.
Thus the tangent matrix K calculated by linearization is derived. We will use these equations again to derive new tangent matrix K
* for the closed loop systems later.It can be shown that for the choice ofγ
/
β
=
2 this algorithm conserves energy exactly and thus is unconditionally stable [12]. Moreover, if we choose γ
=
1 and β
=
1
/
2, the updated angular acceleration A
n
+
1 does not depend on the initial A
n
+
1 during each time interval and we may avoid the accumulated errors from A
n. Moreover, it can also be proved that if, in addition, we use the trapezoidal rule for the applied torque, that is, α
=
1
/
2, then the algorithm is second order accurate both in configuration and velocities.


## 3.2. New Result for Closed Loop System

When we close the loop with the linear stabilizing feedback ofU
=
B
×
W by energy-momentum method mentioned in [20], (12) become
(24)
Λ
˙
=
Λ
W
^
,
J
W
˙
=
(
J
W
+
B
)
×
W
,
which are the equations of motion for the closed loop system in body coordinate.As in Section2.2 we define the convected shifted angular momentum Π
*
=
J
W
+
B. Clearly the norm of  Π
* is constant along the flow. We can also define the Hamiltonian H
=
(
1
/
2
)
(
Π
*
-
B
)
·
J
-
1
(
Π
*
-
B
) in the dual Lie algebra, or as normally H
=
(
1
/
2
)
J
W
·
W in the Lie algebra, and both can be easily seen to be conserved.Since we have closed the system with the constant feedback gain, there is no need to trace the control input and we do not need the parameterα anymore in our algorithm. We only have to set the constant vector B beforehand.In the closed loop energy-momentum conserving algorithm we hope to have the same nice properties as those in the open loop one. We can achieve that by modifying the algorithm in the previous section. First, with the open loop algorithm at hand (24) can be discretized in the form
(25)
J
W
n
+
1
-
J
W
n
=
h
(
J
W
n
+
β
/
γ
+
B
)
×
W
n
+
β
/
γ
between any typical time interval [
t
n
,
t
n
+
1
]
⊂
[
0
,
T
] with the following:
(26)
Π
n
*
=
J
W
n
+
B
,
Π
n
+
1
*
=
J
W
n
+
1
+
B
,
π
n
*
=
Λ
n
(
J
W
n
+
B
)
,
π
n
+
1
*
=
Λ
n
+
1
(
J
W
n
+
1
+
B
)
.
It can be proved that the norm of the shifted angular momentum is conserved in the closed loop algorithm when γ
/
β
=
2.In addition to conservation of norm of the shifted momenta the direction of the shifted spatial angular momentum is also conserved if we enforce thatπ
n
+
1
+
b
n
+
1
=
π
n
+
b
n. To do this we just replace Step 3 in the previous section by
(27)
R
=
Λ
n
+
1
J
W
n
+
1
-
Λ
n
J
W
n
+
Λ
n
+
1
B
-
Λ
n
B
to modify the algorithm. Note that the linear feedback needs not to stabilize the rigid body rotation. In general, any feedback in the form U
=
B
×
W would conserve the shifted spatial angular momentum which is actually the momentum map for the Hamiltonian system. We will see this in some examples later.The conservation of energy depends on properly choosing the parametersα and β as well. From the kinematic equations (16), (17) and the fact that exp
⁡
⁡
(
Θ
^
)
Θ
=
Θ we update
(28)
w
n
+
1
=
w
n
+
Λ
n
W
n
+
1
-
Λ
n
+
1
W
n
+
(
2
-
γ
β
)
(
Λ
n
+
1
-
Λ
n
)
W
n
*
,
where W
n
*
=
W
n
+
(
1
/
2
)
h
A
n. Besides, (16) yields W
n
+
1
+
W
n
=
(
γ
/
β
h
)
Θ
+
(
2
-
γ
/
β
)
W
n
*. By eliminating W
n
+
1 in the previous equations and substituting into the balance equation
(29)
Λ
n
+
1
J
W
n
+
1
+
Λ
n
+
1
B
=
Λ
n
J
W
n
+
Λ
n
B
it turns out to be the following:
(30)
J
W
n
-
B
-
γ
β
h
J
Θ
+
exp
⁡
(
-
Θ
^
)
J
W
n
+
exp
⁡
(
-
Θ
^
)
B
=
(
2
-
γ
β
)
J
W
n
*
.It can also be proved that the energy is conserved in the closed loop algorithm whenγ
/
β
=
2. Again, if we choose the parameters γ
=
1 and β
=
1
/
2, we are free from the error propagation of angular acceleration A
n and the closed loop algorithm would be second order accurate. Moreover, it is unconditionally stable since the energy is conserved.The linearization process of (18) starts from the balance equation (29) with the orientation update (17)
(31)
Λ
n
[
J
W
n
exp
⁡
⁡
(
Θ
^
)
J
(
γ
β
h
Θ
-
W
n
)
+
B
-
exp
⁡
⁡
(
Θ
^
)
B
]
=
0
,
which is a nonlinear equation in Θ. Define
(32)
G
(
Θ
)
∶
=
Λ
n
[
J
W
n
exp
⁡
⁡
(
Θ
^
)
J
(
γ
β
h
Θ
-
W
n
)
k
k
k
k
k
k
+
B
-
exp
⁡
⁡
(
Θ
^
)
B
(
γ
β
h
Θ
-
W
n
)
]
;
then the linearization ofG(Θ) with respect to Θ becomes
(33)
L
[
G
(
Θ
)
]
=
G
(
Θ
0
)
+
D
G
(
Θ
)
·
Δ
Θ
,
where Θ
0 is the fixed point at which the linearization is exerted. By (22) and (23) with the definition (32), we can compute the linearization of G
(
Θ
) along Δ
Θ in body coordinate as follows:
(34)
D
G
(
Θ
)
·
Δ
Θ
=
-
Λ
n
+
1
[
γ
β
h
J
T
(
Θ
)
-
J
W
^
n
+
1
-
B
^
]
Δ
Θ
=
-
K
*
Δ
Θ
.Note that there is a term-
B
^ in the tangent matrix K
* that is different from the result (20) of open loop system. This means we need to put one more term in dealing with the linearization of the balance equation for the closed loop system.Based on the analysis above we treat the constant vectorB as a controller and put it inside the algorithm. Modify the linearization process in the previous section and replace the momenta Π and π by Π
* and π
* to construct the new closed loop algorithm.


## 4. Numerical Examples

There are three examples that we are going to discuss in this section. The first two examples are rigid body rotations about the intermediate and arbitrary axes. The last example is a rotation with the asymptotic stabilizing control.


### 4.1. Intermediate Axis Stabilizing Control

The linear feedbackU
=
B
×
W introduced here is based on the energy-momentum method and is a robustly stable controller in the face of uncertainties of moment of inertia. We will use the open loop algorithm to simulate the tumbling action and then use the closed loop one to stabilize it.Consider a rigid body with moment of inertiaJ
=
diag
⁡
(
5.0
,
4.0
,
3.0
) rotating about the intermediate axis with the convected angular velocity W
0
=
(
0.0
,
20.0
,
0.0
)
T which corresponds to the convected angular momentum Π
0
=
(
0.0
,
80.0
,
0.0
)
T.Initially we perturb the angular velocity a little bit toW
0
=
(
0.01
,
20.01
,
0.01
)
T so that the numerical instability occurs. See Figure 1. Without applied torque the rigid body tumbles and its convected angular momentum Π is shown in Figure 2(a). The energyH and momentum map π are conserved in Figures 2(b) and 2(c), respectively.Figure 1
Tumbling rigid body about intermediate axis without control.Rigid body rotation about intermediate axis without control. (a) Convected angular momentum. (b) Energy. (c) Spatial angular momentum.
(a)
(b)
(c)With the stabilizing vectorB
=
(
0.0
,
40.0
,
0.0
)
T corresponding to λ
=
2 (see [19]) the closed loop algorithm generates the constant shifted momenta Π
* and π
* in Figure 3, while the original π loses its stability. Actually π
* is our momentum map in this case. The energy in Figure 3(c) shows that the closed loop system is Hamiltonian.Rigid body rotation about intermediate axis with controlB
=
(
0.0
,
40.0
,
0.0
)
T. (a) Shifted convected angular momentum. (b) Shifted spatial angular momentum. (c) Energy.
(a)
(b)
(c)Strictly speaking,B
=
(
0.0
,
40.0
,
0.0
)
T is the stabilizing control vector for W
0
=
(
0.0
,
20.0
,
0.0
)
T instead of the perturbed W
0
P; therefore the shifted body angular momentum actually oscillates nearby the axis W
0 due to the initial W
0
P and the oscillation cannot be eliminated by our control vector B. We will resolve this problem by the variable structure control later in Section 4.3.


### 4.2. Arbitrary Axis Stabilizing Control

Consider the same rigid body rotating about an arbitrary axis withW
0
=
(
2.0
,
6.0
,
1.0
)
T and  Π
0
=
(
10.0
,
24.0
,
3.0
)
T.Without applied torque the system’s convected momentumΠ tumbles as shown in Figure 4 and the conserved energyH and momentum map π are shown in Figure 5.Figure 4
Tumbling rigid body about arbitrary axisW
0
=
(
2.0
,
6.0
,
1.0
)
T without control.Rigid body rotation about arbitrary axisW
0
=
(
2.0
,
6.0
,
1.0
)
T without control. (a) Convected angular momentum. (b) Energy. (c) Spatial angular momentum.
(a)
(b)
(c)Applying the stabilizing vectorB
=
(
6.0
,
12.0
,
1.0
)
T, corresponding to λ
=
2 again, in the closed loop algorithm the energyH, momentum map π
*, and stabilized convected momentum Π
* are stable in Figure 6.Rigid body rotation about arbitrary axisW
0
=
(
2.0
,
6.0
,
1.0
)
T with control B
=
(
6.0
,
12.0
,
1.0
)
T. (a) Energy. (b) Shifted spatial angular momentum. (c) Shifted convected angular momentum.
(a)
(b)
(c)In addition, we apply the same vectorB
=
(
6.0
,
12.0
,
1.0
)
T to the rotation of the same system about another arbitrary axis W
0
=
(
6.0
,
2.0
,
0.0
)
T. As expected before, we can observe in Figure 7 that the energyH and momentum map π
* are still conserved, just like the unforced rigid body rotations.Rigid body rotation about arbitrary axisW
0
=
(
6.0
,
2.0
,
0.0
)
T with control B
=
(
6.0
,
12.0
,
1.0
)
T. (a) Energy. (c) Shifted spatial angular momentum.
(a)
(b)


### 4.3. Asymptotic Stabilizing Control

In (7) the quadratic gyroscopic control a
w
×
w modifies the phase flow in order to make some designated axes become equilibria, and the linear part b
×
w shifts the angular momentum to do the same job. In what follows, we consider the special case u
=
a
w
×
w (or U
=
A
W
×
W in body coordinates) by letting b
=
0.If the initial body angular momentum of a closed loop rigid body system is away from any of the stable equilibria, the angular momentum will follow a periodic orbit about that equilibrium on the momentum ball. We say this equilibrium is Lyapunov stable in contrast with asymptotically stable.For different feedback controls there are different shapes of the periodic orbits, especially the orbit of the body angular momentum projected on the plane perpendicular to the stabilized axis. Based on this fact, we can develop an asymptotic stabilizing control to asymptotically pull the initially perturbed angular momentum back to the designated equilibrium axis by using the variable structure control method. We next give an example to clarify this idea.Consider the same rigid body with the initial body angular velocityW
0
=
(
0.0
,
20.0
,
0.0
)
T (or angular momentum Π
0
=
(
0.0
,
80.0
,
0.0
)
T) as in Section 4.1. This system rotates about the intermediate axis and thus is unstable. The projected phase flow around the intermediate axis on Π
1
-
Π
3 plane is plotted in Figure 8. It shows that the intermediate axis is a saddle point.Figure 8
Phase flow near intermediate axis onΠ
1
-
Π
3 plane for a rigid body rotation about intermediate axis without control.With the feedback control lawU
=
A
W
×
W the closed loop equation becomes
(35)
J
W
˙
=
(
J
+
A
)
W
×
W
=
J
c
W
×
W
,
where J
c is the modified moment of inertia.The first quadratic stabilizing control matrix in this example is chosen to be(36)
A
(
1
)
=
diag
⁡
(
-
2
,
0
,
-
1
)
,
and the associate modified moment of inertia and control law are
(37)
J
c
(
1
)
=
diag
⁡
(
3
,
4
,
2
)
,
U
(
1
)
=
[
W
2
W
3
,
W
3
W
1
,
-
2
W
1
W
2
]
T
,respectively. Figure9 shows its phase flow on Π
1
-
Π
3 plane. We can see the rotation about the intermediate axis (the center) has been stabilized and any of other rotations deviating from the center undergoes a periodic motion. Observe that in the first and third quadrants the flow goes away from the origin and in the second and fourth quadrants the flow goes near the center.Figure 9
Phase flow near intermediate axis onΠ
1-Π
3 plane for a rigid body rotation about intermediate axis with control A
(
1
).The second stabilizing control matrix we use is(38)
A
(
2
)
=
diag
⁡
(
-
4
,
0
,
0
)
,
and its modified moment of inertia and control law are
(39)
J
c
(
2
)
=
diag
⁡
(
1
,
4
,
3
)
,
U
(
2
)
=
[
0
,
4
W
3
W
1
,
-
4
W
1
W
2
]
T
,
respectively. We plot the phase flow in Figure 10. Again, we observe that the flow in the first and third quadrants approaches the center and the flow in the second and fourth quadrants goes in the other way.Figure 10
Phase flow near intermediate axis onΠ
1-Π
3 plane for a rigid body rotation about intermediate axis with control A
(
2
).In order to asymptotically stabilize the rotation near the intermediate axis we combine these two control laws and take advantage of both feedbacks in alternative quadrants by switching; that is,(40)
U
=
U
(
1
)
if
Π
1
Π
3
<
0
,
U
=
U
(
2
)
if
Π
1
Π
3
>
0
.Figure11 demonstrates the simulation of the trajectory on Π
1
-
Π
3 plane when starting at Π
0
=
(
20.0
,
80.0
,
0.0
)
T. As we expect the angular momentum asymptotically approaches the intermediate axis. Once the approach achieves some small level, the switching control can be turned off and any of stabilizing control family can be used to keep the intermediate axis rotation stable.Figure 11
Phase flow near intermediate axis onΠ
1
-
Π
3 plane for a rigid body rotation about intermediate axis with switching control consisting of A
(
1
) and A
(
2
).It should be noted that when the angular momentum is away from the equilibrium both control laws do work on the system. Therefore, the energy is changing but converges to a constant associated with the intermediate axis rotation, and the spatial angular momentum has an attitude drift in space. As long as the initial perturbation is not large the deviations in both properties would not be significant.


## 4.1. Intermediate Axis Stabilizing Control

The linear feedbackU
=
B
×
W introduced here is based on the energy-momentum method and is a robustly stable controller in the face of uncertainties of moment of inertia. We will use the open loop algorithm to simulate the tumbling action and then use the closed loop one to stabilize it.Consider a rigid body with moment of inertiaJ
=
diag
⁡
(
5.0
,
4.0
,
3.0
) rotating about the intermediate axis with the convected angular velocity W
0
=
(
0.0
,
20.0
,
0.0
)
T which corresponds to the convected angular momentum Π
0
=
(
0.0
,
80.0
,
0.0
)
T.Initially we perturb the angular velocity a little bit toW
0
=
(
0.01
,
20.01
,
0.01
)
T so that the numerical instability occurs. See Figure 1. Without applied torque the rigid body tumbles and its convected angular momentum Π is shown in Figure 2(a). The energyH and momentum map π are conserved in Figures 2(b) and 2(c), respectively.Figure 1
Tumbling rigid body about intermediate axis without control.Rigid body rotation about intermediate axis without control. (a) Convected angular momentum. (b) Energy. (c) Spatial angular momentum.
(a)
(b)
(c)With the stabilizing vectorB
=
(
0.0
,
40.0
,
0.0
)
T corresponding to λ
=
2 (see [19]) the closed loop algorithm generates the constant shifted momenta Π
* and π
* in Figure 3, while the original π loses its stability. Actually π
* is our momentum map in this case. The energy in Figure 3(c) shows that the closed loop system is Hamiltonian.Rigid body rotation about intermediate axis with controlB
=
(
0.0
,
40.0
,
0.0
)
T. (a) Shifted convected angular momentum. (b) Shifted spatial angular momentum. (c) Energy.
(a)
(b)
(c)Strictly speaking,B
=
(
0.0
,
40.0
,
0.0
)
T is the stabilizing control vector for W
0
=
(
0.0
,
20.0
,
0.0
)
T instead of the perturbed W
0
P; therefore the shifted body angular momentum actually oscillates nearby the axis W
0 due to the initial W
0
P and the oscillation cannot be eliminated by our control vector B. We will resolve this problem by the variable structure control later in Section 4.3.


## 4.2. Arbitrary Axis Stabilizing Control

Consider the same rigid body rotating about an arbitrary axis withW
0
=
(
2.0
,
6.0
,
1.0
)
T and  Π
0
=
(
10.0
,
24.0
,
3.0
)
T.Without applied torque the system’s convected momentumΠ tumbles as shown in Figure 4 and the conserved energyH and momentum map π are shown in Figure 5.Figure 4
Tumbling rigid body about arbitrary axisW
0
=
(
2.0
,
6.0
,
1.0
)
T without control.Rigid body rotation about arbitrary axisW
0
=
(
2.0
,
6.0
,
1.0
)
T without control. (a) Convected angular momentum. (b) Energy. (c) Spatial angular momentum.
(a)
(b)
(c)Applying the stabilizing vectorB
=
(
6.0
,
12.0
,
1.0
)
T, corresponding to λ
=
2 again, in the closed loop algorithm the energyH, momentum map π
*, and stabilized convected momentum Π
* are stable in Figure 6.Rigid body rotation about arbitrary axisW
0
=
(
2.0
,
6.0
,
1.0
)
T with control B
=
(
6.0
,
12.0
,
1.0
)
T. (a) Energy. (b) Shifted spatial angular momentum. (c) Shifted convected angular momentum.
(a)
(b)
(c)In addition, we apply the same vectorB
=
(
6.0
,
12.0
,
1.0
)
T to the rotation of the same system about another arbitrary axis W
0
=
(
6.0
,
2.0
,
0.0
)
T. As expected before, we can observe in Figure 7 that the energyH and momentum map π
* are still conserved, just like the unforced rigid body rotations.Rigid body rotation about arbitrary axisW
0
=
(
6.0
,
2.0
,
0.0
)
T with control B
=
(
6.0
,
12.0
,
1.0
)
T. (a) Energy. (c) Shifted spatial angular momentum.
(a)
(b)


## 4.3. Asymptotic Stabilizing Control

In (7) the quadratic gyroscopic control a
w
×
w modifies the phase flow in order to make some designated axes become equilibria, and the linear part b
×
w shifts the angular momentum to do the same job. In what follows, we consider the special case u
=
a
w
×
w (or U
=
A
W
×
W in body coordinates) by letting b
=
0.If the initial body angular momentum of a closed loop rigid body system is away from any of the stable equilibria, the angular momentum will follow a periodic orbit about that equilibrium on the momentum ball. We say this equilibrium is Lyapunov stable in contrast with asymptotically stable.For different feedback controls there are different shapes of the periodic orbits, especially the orbit of the body angular momentum projected on the plane perpendicular to the stabilized axis. Based on this fact, we can develop an asymptotic stabilizing control to asymptotically pull the initially perturbed angular momentum back to the designated equilibrium axis by using the variable structure control method. We next give an example to clarify this idea.Consider the same rigid body with the initial body angular velocityW
0
=
(
0.0
,
20.0
,
0.0
)
T (or angular momentum Π
0
=
(
0.0
,
80.0
,
0.0
)
T) as in Section 4.1. This system rotates about the intermediate axis and thus is unstable. The projected phase flow around the intermediate axis on Π
1
-
Π
3 plane is plotted in Figure 8. It shows that the intermediate axis is a saddle point.Figure 8
Phase flow near intermediate axis onΠ
1
-
Π
3 plane for a rigid body rotation about intermediate axis without control.With the feedback control lawU
=
A
W
×
W the closed loop equation becomes
(35)
J
W
˙
=
(
J
+
A
)
W
×
W
=
J
c
W
×
W
,
where J
c is the modified moment of inertia.The first quadratic stabilizing control matrix in this example is chosen to be(36)
A
(
1
)
=
diag
⁡
(
-
2
,
0
,
-
1
)
,
and the associate modified moment of inertia and control law are
(37)
J
c
(
1
)
=
diag
⁡
(
3
,
4
,
2
)
,
U
(
1
)
=
[
W
2
W
3
,
W
3
W
1
,
-
2
W
1
W
2
]
T
,respectively. Figure9 shows its phase flow on Π
1
-
Π
3 plane. We can see the rotation about the intermediate axis (the center) has been stabilized and any of other rotations deviating from the center undergoes a periodic motion. Observe that in the first and third quadrants the flow goes away from the origin and in the second and fourth quadrants the flow goes near the center.Figure 9
Phase flow near intermediate axis onΠ
1-Π
3 plane for a rigid body rotation about intermediate axis with control A
(
1
).The second stabilizing control matrix we use is(38)
A
(
2
)
=
diag
⁡
(
-
4
,
0
,
0
)
,
and its modified moment of inertia and control law are
(39)
J
c
(
2
)
=
diag
⁡
(
1
,
4
,
3
)
,
U
(
2
)
=
[
0
,
4
W
3
W
1
,
-
4
W
1
W
2
]
T
,
respectively. We plot the phase flow in Figure 10. Again, we observe that the flow in the first and third quadrants approaches the center and the flow in the second and fourth quadrants goes in the other way.Figure 10
Phase flow near intermediate axis onΠ
1-Π
3 plane for a rigid body rotation about intermediate axis with control A
(
2
).In order to asymptotically stabilize the rotation near the intermediate axis we combine these two control laws and take advantage of both feedbacks in alternative quadrants by switching; that is,(40)
U
=
U
(
1
)
if
Π
1
Π
3
<
0
,
U
=
U
(
2
)
if
Π
1
Π
3
>
0
.Figure11 demonstrates the simulation of the trajectory on Π
1
-
Π
3 plane when starting at Π
0
=
(
20.0
,
80.0
,
0.0
)
T. As we expect the angular momentum asymptotically approaches the intermediate axis. Once the approach achieves some small level, the switching control can be turned off and any of stabilizing control family can be used to keep the intermediate axis rotation stable.Figure 11
Phase flow near intermediate axis onΠ
1
-
Π
3 plane for a rigid body rotation about intermediate axis with switching control consisting of A
(
1
) and A
(
2
).It should be noted that when the angular momentum is away from the equilibrium both control laws do work on the system. Therefore, the energy is changing but converges to a constant associated with the intermediate axis rotation, and the spatial angular momentum has an attitude drift in space. As long as the initial perturbation is not large the deviations in both properties would not be significant.


## 5. Conclusion

In this paper we have outlined the formulations of rigid body dynamics from a geometric point of view for both open loop systems and closed loop ones with a family of the stabilizing linear feedback. It turns out that the momentum map of the closed loop system needs to be modified such that it preserves the analogous conserved quantities as well as the open loop system. To be more specific, we have to translate the spatial angular momentum by an amount of the constant feedback gains B and calculate all properties in this shifted dual to the Lie algebra. We also reviewed the energy-momentum conserving algorithm and made some modifications for the closed loop systems. We derived the conditions of choosing the parameters of Newmark algorithm to make the modified conserving algorithm unconditionally stable and second order accurate and conserve the magnitude and direction of momentum and the energy. A correction term also needs to be placed into the tangent matrix in the linearization process. Several numerical examples have been used to demonstrate the validity of the closed loop algorithm. The energy and shifted momentum map are conserved in these cases, and even if the feedback vectorB does not stabilize the system.



---

*Source: 674910-2014-04-27.xml*