# A Systematic Methodology for Modeling and Attitude Control of Multi-body
Space Telescopes

Yaguang Yang, William Bentz, Lia Lewis Aerospace Engineers, Mission
Engineering and Systems Analysis Division, Goddard Space Flight Center, NASA,
8800 Greenbelt Rd, Greenbelt, 20771 MD.

###### Abstract

This paper derives a symbolic multi-body rigid nonlinear model for a space
telescope using Stoneking’s implementation of Kane’s method. This symbolic
nonlinear model is linearized using Matlab symbolic functions diff and inv
because the analytic linearization is intractable for manual derivation. The
linearized system model is then used to design the controllers using both
linear quadratic regulator (LQR) and robust pole assignment methods. The
closed-loop systems for the two designs are simulated using both the rigid
model as well as a second model containing flexible modes. The performances of
the two designs are compared based on the simulation testing results. Our
conclusion is that the robust pole assignment design offers better performance
than that of the LQR system in terms of actuator usage and pointing accuracy.
However, the LQR approach remains an effective first design step that can
inform the selection of real eigenvalues for robust pole assignment. The
proposed method may be used for the modeling and controller designs for
various multi-body systems.

keywords: Modeling, space telescope control, multibody dynamics, LQR, robust
pole assignment.

##  1 Introduction

The Large UV Optical Infrared Surveyor (LUVOIR) (see Figure 1), to be placed
to Sun-Earth L2 point, is a concept proposed for the key science goal of
characterizing a wide range of exoplanets some of which are potentially
habitable. Although the telescope is still in the concept phase, NASA has
engaged multiple engineering disciplines to conduct preliminary design studies
[1]. The telescope is a typical multi-body dynamical system.

Multi-body dynamical systems can be found in many applications including
machine design, spacecraft dynamics, and robotics. Modern modeling techniques
for multi-body dynamics are based on d’Alembert’s principle in which dynamical
systems were essentially converted into static ones through the introduction
of inertial forces. In 1788, Lagrange formalized this approach by combining
the fundamental ideas of d’Alembert’s principle with explicit descriptions of
virtual work and generalized coordinates [2]. An extension of d’Alembert’s
principle valid for holonomic systems was presented in 1909 by Jourdain [3].
As many engineering systems are nonholonomic, Kane extended d’Alembert’s
principle to this general case in 1961 [4]. Kane’s method has many
applications particularly in robotics for systems of rigid bodies linked by
rotational joints that have an arbitrary number of degrees of freedom (see [5,
6, 7]). Therefore, it is now included in a few engineering handbooks, such as
[5, 8]. Although Kane’s method has become popular, controversy exists
surrounding its originality and efficiency when compared with the Gibbs-Appell
equations [9, 10]. Recently, Piedboeuf indicated that Kane’s equations are
consistent to the Jourdain principle [11]. In [12], Stoneking demonstrated
that Kane’s method can be particularly useful in modeling the case of multiple
rigid bodies connected via rotary joints, e.g., space telescopes. As this was
a conference paper, it was limited in both scope as well as exposure. The
implementation, however, is available as open source software [13].

Using Kane’s method as described by Stoneking [12], Bentz and Lewis derived a
two-body rigid dynamics model for the LUVOIR telescope and simulated the
initial condition response of a LQR design [14]. This work was followed by
similar testing of a higher fidelity three-body rigid dynamics model in [15].
However, the derivation is brief and typos exist (which may mislead the
readers). In addition, the full derivation for the three-body model is left to
but can be difficult for the readers. As linked multi-body systems are widely
seen in robotics and space applications, we extend the preliminary research of
[14, 15] and provide all detailed derivations, which, to our best knowledge,
are not available in any published work. Through the derivation of our three-
body model, we aim to generate further exposure within the aerospace community
of Stoneking’s implementation of Kane’s dynamics and analysis technique that
can efficiently model rigid multiple bodies, connected via rotary joints
having arbitrary degrees of freedom, arranged in tree topologies.

Although plenty of flexible system modeling methods exist (for example [16,
17, 18]), we are particularly interested in rigid model because the rigid
model size is much smaller and its states are normally measurable. Therefore,
the rigid model is more suitable for the control system design than flexible
models, and using a rigid model for the controller design is widely used in
practice. Our ultimate goal is to design a controller for the LUVOIR telescope
in compliance with some arbitrary pointing requirement. As Kane’s multi-body
dynamics are nonlinear, and many powerful control techniques such as LQR and
robust pole assignment are based upon linear models, we have chosen to
linearize the symbolic model for the purpose of controller design. Two
controllers are designed based on the linearized model and their performances
are compared for both rigid and flexible models to give us some confidence
that the designed controller will work for the real system.

There are other multibody modeling methods in the literature. For example, Li
et al [19] discussed a flexible multibody spacecraft modeling which has a
center service module, supporting trusses, and a mirror module. It is assumed
that the center service module and the mirror module are rigid but the trusses
are flexible. In addition, the rigid center service module’s translational
motion is not considered, and connection of the rigid center service module
and the trusses are fixed. Therefore, their model is more specific than the
one discussed in this paper because we do consider translational motion for
all bodies and all connections are not fixed. Hu et al [20] derived a more
general flexible multi-body system modeling method, which has much more
states. Therefore, the model is more suitable for validating the controller
design but is not practical for controller design.

The rest of the paper is organized as follows: Section 2 presents some useful
background and formulas to be used in the rest of the paper. Section 3
provides the detailed derivations of the three-body rigid nonlinear time
invariant dynamics models. Section 4 discusses linearization of the dynamics
models, and two different controller designs, LQR and robust pole assignment,
based on the rigid linearized time invariant dynamics model. Since the LUVOIR
telescope is a flexible structure, we need to validate the designs using both
the rigid linearized time invariant model as well as the higher fidelity
flexible model. Redesigns are necessary if a design works for simulation
testing on the rigid linearized time invariant model but fails the simulation
testing on the high fidelity flexible model. The final LQR and robust pole
assignment designs are tested and their performances for the rigid linearized
time invariant model and flexible model are compared in Section 4 which
highlights the advantages of robust pole assignment.

Figure 1: The concept of LUVOIR telescope.

##  2 Preliminary

This section provides important concepts and formulas in dynamics theory to be
used in this paper and a brief discussion of Kane’s method.

###  2.1 Basic concepts and important formulas

Before we proceed, we present some basic concepts and important formulas which
can be found in [7] and will be used repeatedly in the remainder of the paper.
Let ${\bf B}_{\mathcal{F}}=[{\bf b}_{1},{\bf b}_{2},{\bf b}_{3}]$ be a set of
bases of the frame ${\mathcal{F}}$, then a general vector $\vec{{\bf v}}$
resolved in frame ${\mathcal{F}}$ can be written as:

$\vec{{\bf v}}=v_{1}{\bf b}_{1}+v_{2}{\bf b}_{2}+v_{3}{\bf b}_{3}=[{\bf
b}_{1},{\bf b}_{2},{\bf b}_{3}][v_{1},v_{2},v_{3}]^{{\rm T}}={\bf
B}_{\mathcal{F}}{{\bf v}},$ (1)

where ${\bf v}=[v_{1},v_{2},v_{3}]^{{\rm T}}$. Let $\vec{\omega}_{B/A}$ be the
angular rate of frame B relative to frame A resolved in B. Invoking [7,
(2.3.1)], for any moving vector $\vec{{\bf x}}$ resolved in B, its derivative
in frame A and frame B can be related as:

$\frac{d\vec{{\bf x}}}{dt}\Big{|}_{A}=\frac{d\vec{{\bf
x}}}{dt}\Big{|}_{B}+\vec{\omega}_{B/A}\times\vec{{\bf x}},$ (2)

where $\times$ denotes the cross multiplication of two vectors. If $\vec{{\bf
x}}$ is fixed in frame B, then $\frac{d\vec{{\bf x}}}{dt}\Big{|}_{B}={\bf 0}$.
Therefore, we obtain (see [7, (2.1.2)]),

$\frac{d\vec{{\bf x}}}{dt}\Big{|}_{A}=\vec{\omega}_{B/A}\times\vec{{\bf x}}.$
(3)

The angular velocity of a rigid body B relative to a reference frame A can be
expressed in the following form involving $n$ auxiliary references
$A_{1},\ldots,A_{n}$ [7, (2.4.1)]:

$\vec{\omega}_{B/A}=\vec{\omega}_{B/A_{1}}+\vec{\omega}_{A_{1}/A_{2}}+\cdots+\vec{\omega}_{A_{n}/A}.$
(4)

The angular acceleration of a rigid body B relative to a reference frame A is
defined as the first time-derivative in A of the angular velocity of
$\vec{\omega}_{B/A}$ as [7, (2.5.1)]:

$\vec{\alpha}_{B/A}=\frac{d\vec{\omega}_{B/A}}{dt}.$ (5)

If $P$ and $Q$ are two points fixed on a rigid body B having an angular
velocity $\vec{\omega}_{B/A}$ relative to a reference frame A, then the
velocity of P in A, denoted as $\vec{{\bf v}}_{P/A}$, and the velocity of Q in
A, denoted as $\vec{{\bf v}}_{Q/A}$, are related to each other as [7,
(2.7.1)]:

$\vec{{\bf v}}_{P/A}=\vec{{\bf v}}_{Q/A}+\vec{\omega}_{B/A}\times\vec{{\bf
r}},$ (6)

where $\vec{{\bf r}}$ is the position vector from Q to P. The relationship
between the acceleration of P in A, denoted as $\vec{{\bf a}}_{P/A}$, and the
acceleration of Q in A, denoted as $\vec{{\bf a}}_{Q/A}$, is given as [7,
(2.7.2)]:

$\vec{{\bf a}}_{P/A}=\vec{{\bf
a}}_{Q/A}+\vec{\omega}_{B/A}\times(\vec{\omega}_{B/A}\times\vec{{\bf
r}})+\vec{\alpha}_{B/A}\times\vec{{\bf r}}.$ (7)

###  2.2 Kane’s method

We will derive the three-body rigid nonlinear model for LUVOIR telescope by
using Kane’s method [12]. The notations in this section are defined in [7, 12]
and will become clear for the readers who follow the derivation to the end of
the next section. At that time, readers will see the beauty of Stoneking’s
form of Kane’s method [12]. Let $\\{\tau\\}$ be the general torque vector of
the system, $[{\bf J}]$ be the general inertia matrix of the system,
$\\{\alpha\\}$ be the general angular acceleration vector of the system,
$\\{\omega\\}$ be the general angular rate vector of the system, $\\{{\bf
h}\\}$ be general angular momentum vector of the system, $\\{{\bf f}\\}$ be
the general force vector of the system, $[{\bf M}]$ be the general mass matrix
of the system, $\\{{\bf a}\\}$ be the general linear acceleration vector of
the system, $\Omega$ be the partial angular velocity dyad, and ${\bf V}$ be
the partial velocity dyad ($\Omega$ and ${\bf V}$ will be defined in (24) and
(34)). The Kane’s equation in matrix form can be expressed as

$\Omega^{{\rm T}}\left(\\{\tau\\}-[{\bf J}]\\{\alpha\\}-\\{\omega\times{\bf
h}\\}\right)+{\bf V}^{{\rm T}}\left(\\{{\bf f}\\}-[{\bf M}]\\{{\bf
a}\\}\right)={\bf 0},$ (8)

where the expression in the first parentheses is Euler’s equation, and the
expression in the second parentheses is Newton’s second low of motion.
Therefore, formula (8) appears at first glance to be trivial; however, there
are some significant merits to using Kane’s equation for multi-body models as
discussed in [12]. We will see that the following relations hold in the rest
development.

$\displaystyle\\{\alpha\\}=\Omega\dot{{\bf x}}_{g}+\\{\alpha_{r}\\},$ (9a)
$\displaystyle\\{{\bf a}\\}={\bf V}\dot{{\bf x}}_{g}+\\{{\bf a}_{r}\\},$ (9b)

where ${\bf x}_{g}$ is the generalized speeds of the multi-body system,
$\\{\alpha_{r}\\}$ and $\\{{\bf a}_{r}\\}$ are items that do not include
$\dot{{\bf x}}_{g}$. Substituting equations (9a) and (9b) into (8), and then
grouping on $\dot{{\bf x}}_{g}$ yields Stoneking’s form of Kane’s equation

$\displaystyle\left(\Omega^{{\rm T}}[{\bf J}]\Omega+{\bf V}^{{\rm T}}[{\bf
M}]{\bf V}\right)\dot{{\bf x}}_{g}$ (10) $\displaystyle=$
$\displaystyle\Omega^{{\rm T}}\left(\\{\tau\\}-[{\bf
J}]\\{\alpha_{r}\\}-\\{\omega\times{\bf h}\\}\right)$ $\displaystyle+{\bf
V}^{{\rm T}}\left(\\{{\bf f}\\}-[{\bf M}]\\{{\bf a}_{r}\\}\right),$

which is the rigid multi-body system model. A similar idea was proposed and a
similar formula to (10) is obtained by Hu et al [20] in 2012 for flexible
multi-body system modeling. The advantages of using Kane’s method for
multibody system modeling with tree structure was discussed in [21]. In the
next section, we will provide details of using (10) for rigid multi-body
system modeling.

##  3 Three-body rigid model for LUVOIR telescope

The LUVOIR-A telescope model is assumed to be composed of three rigid bodies
connected in serial by two rotary joints as illustrated in Figure 2.

The three bodies are the spacecraft bus, the boom (tower, or payload
articulation system), and the payload. Spacecraft bus includes many subsystems
such as electrical power system, propulsion, attitude control system,
avionics, command and data handling, thermal management system, mechanical and
structure. The boom can repoint the payload to any position on sky. The
payload includes optical telescope assembly, the high definition imager, the
extreme coronagraph for living planetary system, and ultraviolet multi-object
spectrograph [1]. Several frames of the LUVOIR telescope will be considered.
Let the spacecraft body frame be denoted as $\mathcal{F}_{s}=[{\bf x}_{s},{\bf
y}_{s},{\bf z}_{s}]$, the inertial frame be denoted as $\mathcal{F}_{I}=[{\bf
x}_{I},{\bf y}_{I},{\bf z}_{I}]$, the boom body frame be denoted as
$\mathcal{F}_{b}=[{\bf x}_{b},{\bf y}_{b},{\bf z}_{b}]$, the payload body
frame be denoted as $\mathcal{F}_{p}=[{\bf x}_{p},{\bf y}_{p},{\bf z}_{p}]$.
The spacecraft frame may be defined relative to the inertial frame by

$\mathcal{F}_{s}^{{\rm T}}=\mathcal{O}_{s/I}\mathcal{F}_{I}^{{\rm T}}.$ (11)

where $\mathcal{O}_{s/I}$ is the orientation matrix whose subscript ${s/I}$
represents that the orientation of $\mathcal{F}_{s}$ is relative to
$\mathcal{F}_{I}$. Using standard $3-2-1$ sequence of the intrinsic Euler
angle rotations by yaw angle $\psi$, pitch angle $\theta$, and roll angle
$\phi$, the orientation matrix $\mathcal{O}_{s/I}$ can be expressed as

$\mathcal{O}_{s/I}=\left[\begin{array}[]{ccc}1&0&0\\\
0&\cos(\phi)&\sin(\phi)\\\
0&-\sin(\phi)&\cos(\phi)\end{array}\right]\left[\begin{array}[]{ccc}\cos(\theta)&0&-\sin(\theta)\\\
0&1&0\\\
\sin(\theta)&0&\cos(\theta)\end{array}\right]\left[\begin{array}[]{ccc}\cos(\psi)&\sin(\psi)&0\\\
-\sin(\psi)&\cos(\psi)&0\\\ 0&0&1\end{array}\right].$ (12)

Since orientation matrix $\mathcal{O}_{s/I}$ is an orthogonal matrix, we have

$\mathcal{O}_{I/s}=\mathcal{O}_{s/I}^{{\rm T}}$ (13)

Let the boom gimbal angle be $\gamma$ and the payload gimbal angle be
$\lambda$. The boom body frame to spacecraft body frame orientation matrix can
be expressed as

$\mathcal{O}_{b/s}=\left[\begin{array}[]{ccc}\cos(\gamma)&0&-\sin(\gamma)\\\
0&1&0\\\ \sin(\gamma)&0&\cos(\gamma)\end{array}\right].$ (14)

The payload body frame to boom body frame orientation matrix can be expressed
as

$\mathcal{O}_{p/b}=\left[\begin{array}[]{ccc}\cos(\lambda)&0&-\sin(\lambda)\\\
0&1&0\\\ \sin(\lambda)&0&\cos(\lambda)\end{array}\right].$ (15)

The payload body frame to spacecraft body frame orientation matrix can be
expressed as

$\mathcal{O}_{p/s}=\left[\begin{array}[]{ccc}\cos(\gamma+\lambda)&0&-\sin(\gamma+\lambda)\\\
0&1&0\\\ \sin(\gamma+\lambda)&0&\cos(\gamma+\lambda)\end{array}\right].$ (16)

Figure 2: The description of the three bodies of LUVOIR telescope.

Let the angular velocity of the boom relative to the inertial frame be denoted
as $\vec{\omega}_{b/I}$. We will use similar notations in the remainder of
this paper, for example, $\vec{\omega}_{b/s}$, $\vec{\omega}_{p/b}$, and
$\vec{\omega}_{s/I}$. Let $\Gamma_{1}=[0,1,0]^{{\rm T}}$ and
$\Gamma_{2}=[0,1,0]^{{\rm T}}$, $\sigma_{1}=\dot{\gamma}$ and
$\sigma_{2}=\dot{\lambda}$ be the generalized speeds of the rotary joints of
$G1$ and $G2$. Then the angular rate of the rotary joint $G1$ represented in
the boom frame and the angular rate of the rotary joint $G2$ resolved in the
payload frame can be written as 111For LUVOIR-B where the connection between
payload and boom has two degrees of freedom, the following equations may be
replaced by (24) in [12], but the rest derivation remains essentially the
same.

$\vec{\omega}_{b/s}=\vec{\Gamma}_{1}\sigma_{1},\hskip
7.22743pt\vec{\omega}_{p/b}=\vec{\Gamma}_{2}\sigma_{2}.$ (17)

Using these notations and (4) (as consistent with [7, (2.4.1)]), we have

$\vec{\omega}_{b/I}=\vec{\omega}_{b/s}+\vec{\omega}_{s/I}=\vec{\Gamma}_{1}\sigma_{1}+\vec{\omega}_{s/I}.$
(18)

Let ${\bf B}_{I}$ be the bases of inertial frame, ${\bf B}_{s}$ be the bases
of the spacecraft frame, ${\bf B}_{b}$ be the bases of the boom frame, ${\bf
B}_{p}$ be the bases of the payload frame. Then, we may write (18) as

${\bf B}_{b}{\omega}_{b/I}={\bf B}_{b}\Gamma_{1}\sigma_{1}+{\bf
B}_{s}{\omega}_{s/I}.$ (19)

By premultiplying ${\bf B}_{b}^{{\rm T}}$, we may clear the base dyads and
obtain

${\omega}_{b/I}=\Gamma_{1}\sigma_{1}+\mathcal{O}_{b/s}{\omega}_{s/I}.$ (20)

Similarly,

$\displaystyle\vec{\omega}_{p/I}=\vec{\omega}_{p/b}+\vec{\omega}_{b/s}+\vec{\omega}_{s/I}=\vec{\Gamma}_{2}\sigma_{2}+\vec{\Gamma}_{1}\sigma_{1}+\vec{\omega}_{s/I}.$
(21)

We may write (21) as

$\displaystyle{\bf B}_{p}{\omega}_{p/I}$ $\displaystyle=$ $\displaystyle{\bf
B}_{p}{\omega}_{p/b}+{\bf B}_{b}{\omega}_{b/s}+{\bf B}_{s}{\omega}_{s/I}$ (22)
$\displaystyle=$ $\displaystyle{\bf B}_{p}\Gamma_{2}\sigma_{2}+{\bf
B}_{b}\Gamma_{1}\sigma_{1}+{\bf B}_{s}{\omega}_{s/I},$

and we may clear the base dyads by pre-multiplying ${\bf B}_{p}^{{\rm T}}$ and
obtain

$\displaystyle{\omega}_{p/I}=\Gamma_{2}\sigma_{2}+\mathcal{O}_{p/b}\Gamma_{1}\sigma_{1}+\mathcal{O}_{p/s}{\omega}_{s/I}.$
(23)

Combining (20) and (23) yields

$\left[\begin{array}[]{c}\omega_{s/I}\\\ \omega_{b/I}\\\
\omega_{p/I}\end{array}\right]=\underbrace{\left[\begin{array}[]{cccc}{\bf
I}_{3}&{\bf 0}_{31}&{\bf 0}_{31}&{\bf 0}_{33}\\\
\mathcal{O}_{b/s}&\Gamma_{1}&{\bf 0}_{31}&{\bf 0}_{33}\\\
\mathcal{O}_{p/s}&\mathcal{O}_{p/b}\Gamma_{1}&\Gamma_{2}&{\bf
0}_{33}\end{array}\right]}_{\Omega}\left[\begin{array}[]{c}\omega_{s/I}\\\
\sigma_{1}\\\ \sigma_{2}\\\ {\bf v}_{s/I}\end{array}\right],$ (24)

where ${\bf I}_{3}$ is the three-dimensional identity matrix, ${\bf 0}_{31}$
is the $3\times 1$ all zero matrix, ${\bf 0}_{33}$ is the $3\times 3$ all zero
matrix, and ${\bf v}_{s/I}$ is the velocity vector of the center of mass of
the spacecraft relative to the inertial frame.

Now, we consider the linear velocity of the center of the mass for the boom
and the linear velocity of the center of the mass for the payload. First, we
introduce a notation. For any vector ${\bf a}=[a_{1},a_{2},a_{3}]^{{\rm T}}$,
let a skew symmetric matrix related to ${\bf a}$ be defined as

${\bf a}^{\times}=\left[\begin{array}[]{ccc}0&-a_{3}&a_{2}\\\
a_{3}&0&-a_{1}\\\ -a_{2}&a_{1}&0\end{array}\right].$ (25)

The cross product of two vectors ${\bf a}\times{\bf b}$ can be written as a
multiplication of the matrix ${\bf a}^{\times}$ and the vector ${\bf b}$,
i.e., ${\bf a}^{\times}{\bf b}$. Let $\vec{{\bf v}}_{s/I}$ be the velocity of
the center of mass of the spacecraft in the inertial frame, $\vec{{\bf
v}}_{b/I}$ be the velocity of the center of mass of the boom in the inertial
frame, $\vec{{\bf v}}_{G1/I}$ be the velocity of $G1$ in the inertial frame
(see Figure 2), $\vec{{\bf r}}_{G1/s}$ be the position vector from the center
of mass of the spacecraft to the joint $G1$, $\vec{{\bf r}}_{G1/b}$ be the
position vector from the center of mass of the boom to the joint $G1$,
$\vec{{\bf v}}_{p/I}$ be the velocity of the center of mass for the payload in
the inertial frame, $\vec{{\bf v}}_{G2/I}$ be the velocity of $G2$ in the
inertial frame, $\vec{{\bf r}}_{G2/p}$ be the position vector from the center
of mass of the payload to the joint $G2$, $\vec{{\bf r}}_{G2/b}$ be the
position vector from the center of mass of the boom to the joint $G2$. Note
that all these $\vec{{\bf v}}$ and $\vec{{\bf r}}$ vectors are in the inertial
frame (see Figure 2). Since $G1$ is a point on both the spacecraft and the
boom, from (6) (see [7, (2.7.1)]), we have

$\displaystyle\vec{{\bf v}}_{G1/I}=\vec{{\bf
v}}_{s/I}+\vec{\omega}_{s/I}\times\vec{{\bf r}}_{G1/s},$ (26a)
$\displaystyle\vec{{\bf v}}_{G1/I}=\vec{{\bf
v}}_{b/I}+\vec{\omega}_{b/I}\times\vec{{\bf r}}_{G1/b}.$ (26b)

Substituting (26a) into (26b) and invoking (18) yield

$\displaystyle\vec{{\bf v}}_{b/I}$ $\displaystyle=$ $\displaystyle\vec{{\bf
v}}_{s/I}+\vec{\omega}_{s/I}\times\vec{{\bf
r}}_{G1/s}-\vec{\omega}_{b/I}\times\vec{{\bf r}}_{G1/b}$ $\displaystyle=$
$\displaystyle\vec{{\bf v}}_{s/I}+\vec{\omega}_{s/I}\times\vec{{\bf
r}}_{G1/s}-(\vec{\omega}_{b/s}+\vec{\omega}_{s/I})\times\vec{{\bf r}}_{G1/b}$

We may represent each vector in an appropriate basis and write (LABEL:G1) as

$\displaystyle{\bf B}_{I}{{\bf v}}_{b/I}$ $\displaystyle=$ $\displaystyle{\bf
B}_{I}{{\bf v}}_{s/I}+{\bf B}_{s}{\omega}_{s/I}\times{\bf B}_{I}{{\bf
r}}_{G1/s}$ (28) $\displaystyle-({\bf B}_{b}\Gamma_{1}\sigma_{1}+{\bf
B}_{s}{\omega}_{s/I})\times{\bf B}_{I}{{\bf r}}_{G1/b}$

Using the notations that ${{\bf r}}_{G1/s}|_{I}={\bf B}_{I}{{\bf r}}_{G1/s}$
(where ${{\bf r}}_{G1/s}|_{I}$ means that the vector ${{\bf r}}_{G1/s}$ is
expressed in the inertial frame) and ${{\bf r}}_{G1/b}|_{I}={\bf B}_{I}{{\bf
r}}_{G1/b}$, we may clear the base dyads by pre-multiplying ${\bf B}_{I}^{{\rm
T}}$ and obtain

$\displaystyle{{\bf v}}_{b/I}$ $\displaystyle=$ $\displaystyle{{\bf
v}}_{s/I}-{{\bf r}}_{G1/s}|_{I}\times\mathcal{O}_{I/s}{\omega}_{s/I}+{{\bf
r}}_{G1/b}|_{I}\times\mathcal{O}_{I/b}\Gamma_{1}\sigma_{1}$ (29)
$\displaystyle+{{\bf r}}_{G1/b}|_{I}\times\mathcal{O}_{I/s}{\omega}_{s/I}$
$\displaystyle=$ $\displaystyle{\bf v}_{s/I}|_{I}+\underbrace{\left[{\bf
r}_{G1/b}|_{I}^{\times}-{\bf
r}_{G1/s}|_{I}^{\times}\right]\mathcal{O}_{I/s}}_{{\bf V}_{21}\in{\bf
R}^{3\times 3}}\omega_{s/I}$ $\displaystyle+\underbrace{{\bf
r}_{G1/b}|_{I}^{\times}\mathcal{O}_{I/b}\Gamma_{1}}_{{\bf v}_{22}\in{\bf
R}^{3\times 1}}\sigma_{1}.$

Applying the same idea to the joint $G2$ and invoking (6) (see [7, (2.7.1)]),
we have

$\displaystyle\vec{{\bf v}}_{G2/I}=\vec{{\bf
v}}_{b/I}+\vec{\omega}_{b/I}\times\vec{{\bf r}}_{G2/b},$ (30a)
$\displaystyle\vec{{\bf v}}_{G2/I}=\vec{{\bf
v}}_{p/I}+\vec{\omega}_{p/I}\times\vec{{\bf r}}_{G2/p}.$ (30b)

Substituting (30a) into (30b) and invoking (21) yield

$\displaystyle\vec{{\bf v}}_{p/I}$ $\displaystyle=$ $\displaystyle\vec{{\bf
v}}_{b/I}+\vec{\omega}_{b/I}\times\vec{{\bf
r}}_{G2/b}-\vec{\omega}_{p/I}\times\vec{{\bf r}}_{G2/p}$ (31) $\displaystyle=$
$\displaystyle\vec{{\bf
v}}_{b/I}+(\vec{\omega}_{b/s}+\vec{\omega}_{s/I})\times\vec{{\bf r}}_{G2/b}$
$\displaystyle-(\vec{\omega}_{p/b}+\vec{\omega}_{b/s}+\vec{\omega}_{s/I})\times\vec{{\bf
r}}_{G2/p}$ $\displaystyle=$ $\displaystyle\vec{{\bf
v}}_{b/I}-\vec{\omega}_{p/b}\times\vec{{\bf
r}}_{G2/p}+\vec{\omega}_{b/s}\times(\vec{{\bf r}}_{G2/b}-\vec{{\bf
r}}_{G2/p})$ $\displaystyle+\vec{\omega}_{s/I}\times(\vec{{\bf
r}}_{G2/b}-\vec{{\bf r}}_{G2/p}).$

We may represent each vector in an appropriate basis and write (31) as

$\displaystyle{\bf B}_{I}{{\bf v}}_{p/I}$ $\displaystyle=$ $\displaystyle{\bf
B}_{I}{{\bf v}}_{b/I}-{\bf B}_{p}\Gamma_{2}\sigma_{2}\times{\bf B}_{I}{{\bf
r}}_{G2/p}$ (32) $\displaystyle+{\bf B}_{b}\Gamma_{1}\sigma_{1}\times({\bf
B}_{I}{{\bf r}}_{G2/b}-{\bf B}_{I}{{\bf r}}_{G2/p})$ $\displaystyle+{\bf
B}_{s}{\omega}_{s/I}\times({\bf B}_{I}{{\bf r}}_{G2/b}-{\bf B}_{I}{{\bf
r}}_{G2/p})$

Using the notations that ${{\bf r}}_{G2/b}|_{I}={\bf B}_{I}{{\bf r}}_{G2/b}$
and ${{\bf r}}_{G2/p}|_{I}={\bf B}_{I}{{\bf r}}_{G2/p}$, and invoking (29), we
may clear the base dyads by pre-multiplying ${\bf B}_{I}^{{\rm T}}$ in (32)
and obtain

$\displaystyle{{\bf v}}_{p/I}$ $\displaystyle=$ $\displaystyle{{\bf
v}}_{b/I}+{{\bf
r}}_{G2/p}|_{I}\times\mathcal{O}_{I/p}\Gamma_{2}\sigma_{2}+({{\bf
r}}_{G2/p}|_{I}-{{\bf
r}}_{G2/b}|_{I})\times\mathcal{O}_{I/b}\Gamma_{1}\sigma_{1}$ (33)
$\displaystyle+({{\bf r}}_{G2/p}|_{I}-{{\bf
r}}_{G2/b}|_{I})\times\mathcal{O}_{I/s}{\omega}_{s/I}$ $\displaystyle=$
$\displaystyle{{\bf v}}_{s/I}+{\bf V}_{21}{\omega}_{s/I}+{\bf
v}_{22}\sigma_{1}+{{\bf
r}}_{G2/p}|_{I}\times\mathcal{O}_{I/p}\Gamma_{2}\sigma_{2}$
$\displaystyle+({{\bf r}}_{G2/p}|_{I}-{{\bf
r}}_{G2/b}|_{I})\times\mathcal{O}_{I/b}\Gamma_{1}\sigma_{1}+({{\bf
r}}_{G2/p}|_{I}-{{\bf r}}_{G2/b}|_{I})\times\mathcal{O}_{I/s}{\omega}_{s/I}$
$\displaystyle=$ $\displaystyle{\bf v}_{s/I}+\underbrace{\left[{\bf
r}_{G1/b}|_{I}^{\times}-{\bf r}_{G1/s}|_{I}^{\times}+{\bf
r}_{G2/p}|_{I}^{\times}-{\bf
r}_{G2/b}|_{I}^{\times}\right]\mathcal{O}_{I/s}}_{{\bf V}_{31}\in{\bf
R}^{3\times 3}}(\omega_{s/I})$ $\displaystyle+\underbrace{\left[{\bf
r}_{G1/b}|_{I}^{\times}+{\bf r}_{G2/p}|_{I}^{\times}-{\bf
r}_{G2/b}|_{I}^{\times}\right]\mathcal{O}_{I/b}\Gamma_{1}}_{{\bf
v}_{32}\in{\bf R}^{3\times 1}}\sigma_{1}+\underbrace{{\bf
r}_{G2/p}|_{I}^{\times}\mathcal{O}_{I/p}\Gamma_{2}}_{{\bf v}_{33}\in{\bf
R}^{3\times 1}}\sigma_{2}.$

Combining (29) and (33) yields

$\left[\begin{array}[]{c}{\bf v}_{s/I}\\\ {\bf v}_{b/I}\\\ {\bf
v}_{p/I}\end{array}\right]=\underbrace{\left[\begin{array}[]{cccc}{\bf
0}_{33}&{\bf 0}_{31}&{\bf 0}_{31}&{\bf I}_{3}\\\ {\bf V}_{21}&{\bf
v}_{22}&{\bf 0}_{31}&{\bf I}_{3}\\\ {\bf V}_{31}&{\bf v}_{32}&{\bf
v}_{33}&{\bf I}_{3}\end{array}\right]}_{{\bf
V}}\left[\begin{array}[]{c}\omega_{s/I}\\\ \sigma_{1}\\\ \sigma_{2}\\\ {\bf
v}_{s/I}\end{array}\right].$ (34)

In the sequel, we show that $\Omega$, defined in (24), and ${\bf V}$, defined
in (34), are the same ones defined in the differential equations (9) which
will be used to obtain the multi-body system model (10). Let
$\vec{\alpha}_{s/I}$ be the angular acceleration of the center of the mass of
the spacecraft relative to the inertial frame, $\vec{\alpha}_{b/I}$ be the
angular acceleration of the center of the mass of the boom relative to the
inertial frame, and $\vec{\alpha}_{p/I}$ be the angular acceleration of the
center of the mass of the payload relative to the inertial frame,
respectively. Taking the derivative for (18) and invoking (2), we have

$\displaystyle\vec{\alpha}_{b/I}$ $\displaystyle=$
$\displaystyle\frac{d\vec{\omega}_{b/I}}{dt}=\frac{d(\vec{\omega}_{b/s}+\vec{\omega}_{s/I})}{dt}$
(35) $\displaystyle=$
$\displaystyle\frac{d\vec{\omega}_{b/s}}{dt}+\frac{d\vec{\omega}_{s/I}}{dt}$
$\displaystyle=$
$\displaystyle\vec{\Gamma}_{1}\dot{\sigma_{1}}+\vec{\omega}_{b/I}\times\vec{\Gamma}_{1}\sigma_{1}+\vec{\alpha}_{s/I}.$

Representing each vector in an appropriate base yields

$\displaystyle{\bf B}_{b}{\alpha}_{b/I}={\bf
B}_{b}\Gamma_{1}\dot{\sigma_{1}}+{\bf B}_{b}{\omega}_{b/I}\times{\bf
B}_{b}\Gamma_{1}\sigma_{1}+{\bf B}_{s}{\alpha}_{s/I}.$ (36)

Pre-multiplying ${\bf B}_{b}^{{\rm T}}$ on both sides of (36) clears the base
dyads and yields

$\displaystyle{\alpha}_{b/I}$ $\displaystyle=$
$\displaystyle\Gamma_{1}\dot{\sigma_{1}}+{\omega}_{b/I}\times\Gamma_{1}\sigma_{1}+\mathcal{O}_{b/s}{\alpha}_{s/I}$
(37) $\displaystyle=$
$\displaystyle\Gamma_{1}\dot{\sigma_{1}}+\mathcal{O}_{b/s}{\alpha}_{s/I}+{\alpha}_{b/I}^{r},$

where

${\alpha}_{b/I}^{r}={\omega}_{b/I}\times\Gamma_{1}\sigma_{1},\hskip
7.22743pt{\omega}_{b/I}=\mathcal{O}_{b/s}\omega_{s/I}.$ (38)

Taking the derivative for (21) and invoking (2) yields

$\displaystyle\vec{\alpha}_{p/I}$ $\displaystyle=$
$\displaystyle\frac{d\vec{\omega}_{p/I}}{dt}=\frac{d(\vec{\omega}_{p/b}+\vec{\omega}_{b/s}+\vec{\omega}_{s/I})}{dt}$
$\displaystyle=$
$\displaystyle\frac{d\vec{\omega}_{p/b}}{dt}+\frac{d\vec{\omega}_{b/s}}{dt}+\frac{d\vec{\omega}_{s/I}}{dt}$
$\displaystyle=$
$\displaystyle\vec{\Gamma}_{2}\dot{\sigma}_{2}+\vec{\omega}_{p/I}\times\vec{\Gamma}_{2}\sigma_{2}+\vec{\Gamma}_{1}\dot{\sigma}_{1}$
$\displaystyle+\vec{\omega}_{b/I}\times\vec{\Gamma}_{1}\sigma_{1}+\vec{\alpha}_{s/I}.$

Representing each vector in an appropriate base yields

$\displaystyle{\bf B}_{p}{\alpha}_{p/I}$ $\displaystyle=$ $\displaystyle{\bf
B}_{p}\Gamma_{2}\dot{\sigma}_{2}+{\bf
B}_{p}{\omega}_{p/I}\times\Gamma_{2}\sigma_{2}+{\bf
B}_{b}\Gamma_{1}\dot{\sigma}_{1}$ (40) $\displaystyle+{\bf
B}_{b}{\omega}_{b/I}\times\Gamma_{1}\sigma_{1}+{\bf B}_{s}{\alpha}_{s/I}.$

Pre-multiplying ${\bf B}_{p}^{{\rm T}}$ on both sides of (40) clears the base
dyads and yields

$\displaystyle{\alpha}_{p/I}$ $\displaystyle=$
$\displaystyle\Gamma_{2}\dot{\sigma}_{2}+{\omega}_{p/I}\times\Gamma_{2}\sigma_{2}+\mathcal{O}_{p/b}\Gamma_{1}\dot{\sigma}_{1}$
(41)
$\displaystyle+\mathcal{O}_{p/b}{\omega}_{b/I}\times\Gamma_{1}\sigma_{1}+\mathcal{O}_{p/s}\alpha_{s/I}$
$\displaystyle=$
$\displaystyle\Gamma_{2}\dot{\sigma}_{2}+\mathcal{O}_{p/b}\Gamma_{1}\dot{\sigma}_{1}+\mathcal{O}_{p/s}\alpha_{s/I}+{\alpha}_{p/I}^{r}.$

where

$\displaystyle{\alpha}_{p/I}^{r}$ $\displaystyle=$
$\displaystyle{\omega}_{p/I}\times\Gamma_{2}\sigma_{2}+\mathcal{O}_{p/b}{\omega}_{b/I}\times\Gamma_{1}\sigma_{1}$
$\displaystyle=$
$\displaystyle{\omega}_{p/I}\times\Gamma_{2}\sigma_{2}+\mathcal{O}_{p/b}{\alpha}_{b/I}^{r},\hskip
7.22743pt{\omega}_{p/I}=\mathcal{O}_{p/b}{\omega}_{b/I}.$

Denote $\alpha_{1}=\dot{\sigma}_{1}$, $\alpha_{2}=\dot{\sigma}_{2}$, and

$\dot{{\bf
x}}_{g}=[\dot{\omega}_{s/I},~{}\dot{\sigma}_{1},~{}\dot{\sigma}_{2},~{}\dot{{\bf
v}}_{s/I}]^{{\rm T}}.$ (43)

Combining (37) and (41) yields

$\displaystyle\\{\alpha\\}:=\left[\begin{array}[]{c}\alpha_{s/I}\\\
\alpha_{b/I}\\\ \alpha_{p/I}\end{array}\right]$ (47) $\displaystyle=$
$\displaystyle\underbrace{\left[\begin{array}[]{cccc}{\bf I}_{3}&{\bf
0}_{31}&{\bf 0}_{31}&{\bf 0}_{33}\\\ \mathcal{O}_{b/s}&\Gamma_{1}&{\bf
0}_{31}&{\bf 0}_{33}\\\
\mathcal{O}_{p/s}&\mathcal{O}_{p/b}\Gamma_{1}&\Gamma_{2}&{\bf
0}_{33}\end{array}\right]}_{\Omega}\underbrace{\left[\begin{array}[]{c}\dot{\omega}_{s/I}\\\
\dot{\sigma}_{1}\\\ \dot{\sigma}_{2}\\\ \dot{{\bf
v}}_{s/I}\end{array}\right]}_{\dot{{\bf
x}}_{g}}+\underbrace{\left[\begin{array}[]{c}{\bf 0}_{31}\\\
{\alpha}_{b/I}^{r}\\\
{\alpha}_{p/I}^{r}\end{array}\right]}_{\\{\alpha_{r}\\}}$ (58)
$\displaystyle=$ $\displaystyle\Omega\dot{{\bf x}}_{g}+\\{\alpha_{r}\\},$ (59)

which is equivalent to (9a). We also showed that $\Omega$ defined in (24) is
the same as the one defined in (59) or in (9a). Next, we derive equation (9b).
Let $\vec{{\bf a}}_{s/I}$ be the acceleration of the center of the mass of the
spacecraft relative to the inertial frame, $\vec{{\bf a}}_{b/I}$ be the
acceleration of the center of the mass of the boom relative to the inertial
frame, $\vec{{\bf a}}_{p/I}$ be the acceleration of the center of the mass of
the payload relative to the inertial frame, $\vec{{\bf a}}_{G1/I}$ be the
acceleration of the joint $G1$ relative to the inertial frame, and $\vec{{\bf
a}}_{G2/I}$ be the acceleration of the joint $G2$ relative to the inertial
frame, respectively. Since $G1$ is a point on both the spacecraft and the
boom, applying (7) to the joint $G1$, we have

$\displaystyle\vec{{\bf a}}_{G1/I}=\vec{{\bf
a}}_{s/I}+\vec{\omega}_{s/I}\times(\vec{\omega}_{s/I}\times\vec{{\bf
r}}_{G1/s})+\vec{\alpha}_{s/I}\times\vec{{\bf r}}_{G1/s},$ (60a)
$\displaystyle\vec{{\bf a}}_{G1/I}=\vec{{\bf
a}}_{b/I}+\vec{\omega}_{b/I}\times(\vec{\omega}_{b/I}\times\vec{{\bf
r}}_{G1/b})+\vec{\alpha}_{b/I}\times\vec{{\bf r}}_{G1/b}.$ (60b)

Substituting (60a) into (60b) yields

$\displaystyle\vec{{\bf a}}_{b/I}$ $\displaystyle=$ $\displaystyle\vec{{\bf
a}}_{s/I}+\vec{\omega}_{s/I}\times(\vec{\omega}_{s/I}\times\vec{{\bf
r}}_{G1/s})+\vec{\alpha}_{s/I}\times\vec{{\bf r}}_{G1/s}$ (61)
$\displaystyle-\vec{\omega}_{b/I}\times(\vec{\omega}_{b/I}\times\vec{{\bf
r}}_{G1/b})-\vec{\alpha}_{b/I}\times\vec{{\bf r}}_{G1/b}.$

Representing each vector in an appropriate base and invoking (36) yields

$\displaystyle{\bf B}_{I}{{\bf a}}_{b/I}$ $\displaystyle=$ $\displaystyle{\bf
B}_{I}{{\bf a}}_{s/I}+{\bf B}_{s}{\omega}_{s/I}\times({\bf
B}_{s}{\omega}_{s/I}\times{\bf B}_{I}{{\bf r}}_{G1/s})$ $\displaystyle+{\bf
B}_{s}{\alpha}_{s/I}\times{\bf B}_{I}{{\bf r}}_{G1/s}$ $\displaystyle-{\bf
B}_{b}{\omega}_{b/I}\times({\bf B}_{b}{\omega}_{b/I}\times{\bf B}_{I}{{\bf
r}}_{G1/b})$ $\displaystyle-{\bf B}_{b}{\alpha}_{b/I}\times{\bf B}_{I}{{\bf
r}}_{G1/b}$ $\displaystyle=$ $\displaystyle{\bf B}_{I}{{\bf a}}_{s/I}+{\bf
B}_{s}{\omega}_{s/I}\times({\bf B}_{s}{\omega}_{s/I}\times{\bf B}_{I}{{\bf
r}}_{G1/s})$ $\displaystyle+{\bf B}_{s}{\alpha}_{s/I}\times{\bf B}_{I}{{\bf
r}}_{G1/s}$ $\displaystyle-{\bf B}_{b}{\omega}_{b/I}\times({\bf
B}_{b}{\omega}_{b/I}\times{\bf B}_{I}{{\bf r}}_{G1/b})$ $\displaystyle-({\bf
B}_{b}\Gamma_{1}\dot{\sigma_{1}}+{\bf
B}_{b}{\omega}_{b/I}\times\Gamma_{1}\sigma_{1}+{\bf
B}_{s}{\alpha}_{s/I})\times{\bf B}_{I}{{\bf r}}_{G1/b}.$

Pre-multiplying ${\bf B}_{I}^{{\rm T}}$ on both sides of
(LABEL:accelerationB1) clears the base dyads and yields

$\displaystyle{{\bf a}}_{b/I}$ (63) $\displaystyle=$ $\displaystyle{{\bf
a}}_{s/I}+\mathcal{O}_{I/s}{\omega}_{s/I}\times({\omega}_{s/I}\times{{\bf
r}}_{G1/s}|_{I})$ $\displaystyle+\mathcal{O}_{I/s}{\alpha}_{s/I}\times{{\bf
r}}_{G1/s}|_{I}$
$\displaystyle-\mathcal{O}_{I/b}{\omega}_{b/I}\times({\omega}_{b/I}\times{{\bf
r}}_{G1/b}|_{I})$
$\displaystyle-(\mathcal{O}_{I/b}\Gamma_{1}\dot{\sigma_{1}}+\mathcal{O}_{I/b}{\omega}_{b/I}\times\Gamma_{1}\sigma_{1}+\mathcal{O}_{I/s}{\alpha}_{s/I})\times{{\bf
r}}_{G1/b}|_{I}$ $\displaystyle=$ $\displaystyle{{\bf
a}}_{s/I}+\underbrace{({{\bf r}}_{G1/b}|_{I}-{{\bf
r}}_{G1/s}|_{I})\times\mathcal{O}_{I/s}}_{{\bf V}_{21}\in{\bf R}^{3\times
3}}{\alpha}_{s/I}$ $\displaystyle+\underbrace{{{\bf
r}}_{G1/b}|_{I}\times\mathcal{O}_{I/b}\Gamma_{1}}_{{\bf v}_{22}\in{\bf
R}^{3\times 1}}\dot{\sigma_{1}}+{{\bf a}}_{b/I}^{r},$

where

$\displaystyle{{\bf a}}_{b/I}^{r}$ $\displaystyle=$
$\displaystyle\mathcal{O}_{I/s}{\omega}_{I/s}\times({\omega}_{s/I}\times{{\bf
r}}_{G1/s}|_{I})$ (64)
$\displaystyle-\mathcal{O}_{I/b}{\omega}_{b/I}\times({\omega}_{b/I}\times{{\bf
r}}_{G1/b}|_{I})$ $\displaystyle+{{\bf
r}}_{G1/b}|_{I}\times(\mathcal{O}_{I/b}{\omega}_{b/I}\times\Gamma_{1}\sigma_{1})$
$\displaystyle=$
$\displaystyle\mathcal{O}_{I/s}{\omega}_{I/s}^{\times}({\omega}_{s/I}^{\times}{{\bf
r}}_{G1/s}|_{I})-\mathcal{O}_{I/b}{\omega}_{b/I}^{\times}({\omega}_{b/I}^{\times}{{\bf
r}}_{G1/b}|_{I})$ $\displaystyle+{{\bf
r}}_{G1/b}|_{I}^{\times}\mathcal{O}_{I/b}\alpha_{b/I}^{r}.$

Applying (7) to the joint $G2$, we have

$\displaystyle\vec{{\bf a}}_{G2/I}=\vec{{\bf
a}}_{b/I}+\vec{\omega}_{b/I}\times(\vec{\omega}_{b/I}\times\vec{{\bf
r}}_{G2/b})+\vec{\alpha}_{b/I}\times\vec{{\bf r}}_{G2/b},$ (65a)
$\displaystyle\vec{{\bf a}}_{G2/I}=\vec{{\bf
a}}_{p/I}+\vec{\omega}_{p/I}\times(\vec{\omega}_{p/I}\times\vec{{\bf
r}}_{G2/p})+\vec{\alpha}_{p/I}\times\vec{{\bf r}}_{G2/p}.$ (65b)

Substituting (65a) into (65b) yields

$\displaystyle\vec{{\bf a}}_{p/I}$ $\displaystyle=$ $\displaystyle\vec{{\bf
a}}_{b/I}+\vec{\omega}_{b/I}\times(\vec{\omega}_{b/I}\times\vec{{\bf
r}}_{G2/b})+\vec{\alpha}_{b/I}\times\vec{{\bf r}}_{G2/b}$ (66)
$\displaystyle-\vec{\omega}_{p/I}\times(\vec{\omega}_{p/I}\times\vec{{\bf
r}}_{G2/p})-\vec{\alpha}_{p/I}\times\vec{{\bf r}}_{G2/p}.$

Representing each vector in an appropriate base yields

$\displaystyle{\bf B}_{I}{{\bf a}}_{p/I}$ $\displaystyle=$ $\displaystyle{\bf
B}_{I}{{\bf a}}_{b/I}+{\bf B}_{b}{\omega}_{b/I}\times({\bf
B}_{b}{\omega}_{b/I}\times{\bf B}_{I}{{\bf r}}_{G2/b})$ (67)
$\displaystyle+{\bf B}_{b}{\alpha}_{b/I}\times{\bf B}_{I}{{\bf r}}_{G2/b}$
$\displaystyle-{\bf B}_{p}{\omega}_{p/I}\times({\bf
B}_{p}{\omega}_{p/I}\times{\bf B}_{I}{{\bf r}}_{G2/p})$ $\displaystyle-{\bf
B}_{p}{\alpha}_{p/I}\times{\bf B}_{I}{{\bf r}}_{G2/p}.$

Pre-multiplying ${\bf B}_{I}^{{\rm T}}$ on both sides of (67) to clear the
base dyads, and substituting (37), (41), and (63) into the formula yields

$\displaystyle{{\bf a}}_{p/I}$ $\displaystyle=$ $\displaystyle{{\bf
a}}_{b/I}+\mathcal{O}_{I/b}{\omega}_{b/I}\times({\omega}_{b/I}\times{{\bf
r}}_{G2/b}|_{I})+\mathcal{O}_{I/b}{\alpha}_{b/I}\times{{\bf r}}_{G2/b}|_{I}$
(68)
$\displaystyle-\mathcal{O}_{I/p}{\omega}_{p/I}\times({\omega}_{p/I}\times{{\bf
r}}_{G2/p}|_{I})-\mathcal{O}_{I/p}{\alpha}_{p/I}\times{{\bf r}}_{G2/p}|_{I}$
$\displaystyle=$ $\displaystyle{{\bf a}}_{s/I}+({{\bf r}}_{G1/b}|_{I}-{{\bf
r}}_{G1/s}|_{I})\times\mathcal{O}_{I/s}{\alpha}_{s/I}+{{\bf
r}}_{G1/b}|_{I}\times\mathcal{O}_{I/b}\Gamma_{1}\dot{\sigma_{1}}+{{\bf
a}}_{b/I}^{r}$
$\displaystyle+\mathcal{O}_{I/b}{\omega}_{b/I}\times({\omega}_{b/I}\times{{\bf
r}}_{G2/b}|_{I})+\mathcal{O}_{I/b}(\Gamma_{1}\dot{\sigma_{1}}+\mathcal{O}_{b/s}{\alpha}_{s/I}+{\alpha}_{b/I}^{r})\times{{\bf
r}}_{G2/b}|_{I}$
$\displaystyle-\mathcal{O}_{I/p}{\omega}_{p/I}\times({\omega}_{p/I}\times{{\bf
r}}_{G2/p}|_{I})-\mathcal{O}_{I/p}(\Gamma_{2}\dot{\sigma}_{2}+\mathcal{O}_{p/b}\Gamma_{1}\dot{\sigma}_{1}+\mathcal{O}_{p/s}\alpha_{s/I}+{\alpha}_{p/I}^{r})\times{{\bf
r}}_{G2/p}|_{I}$ $\displaystyle=$ $\displaystyle{{\bf
a}}_{s/I}+\underbrace{({{\bf r}}_{G1/b}|_{I}-{{\bf r}}_{G1/s}|_{I}+{{\bf
r}}_{G2/p}|_{I}-{{\bf r}}_{G2/b}|_{I})\times\mathcal{O}_{I/s}}_{{\bf
V}_{31}\in{\bf R}^{3\times 3}}{\alpha}_{s/I}$
$\displaystyle+\underbrace{({{\bf r}}_{G1/b}|_{I}+{{\bf r}}_{G2/p}|_{I}-{{\bf
r}}_{G2/b}|_{I})\times\mathcal{O}_{I/b}\Gamma_{1}}_{{\bf v}_{32}\in{\bf
R}^{3\times 1}}\dot{\sigma}_{1}+\underbrace{{{\bf
r}}_{G2/p}|_{I}\times\mathcal{O}_{I/p}\Gamma_{2}}_{{\bf v}_{33}\in{\bf
R}^{3\times 1}}\dot{\sigma}_{2}+{{\bf a}}_{p/I}^{r},$

where

$\displaystyle{{\bf a}}_{p/I}^{r}$ $\displaystyle=$ $\displaystyle{{\bf
a}}_{b/I}^{r}-{{\bf r}}_{G2/b}|_{I}\times\mathcal{O}_{I/b}{\alpha}_{b/I}^{r}$
(69) $\displaystyle+{{\bf
r}}_{G2/p}|_{I}\times\mathcal{O}_{I/p}{\alpha}_{p/I}^{r}$
$\displaystyle+\mathcal{O}_{I/b}{\omega}_{b/I}\times({\omega}_{b/I}\times{{\bf
r}}_{G2/b}|_{I})$
$\displaystyle-\mathcal{O}_{I/p}{\omega}_{p/I}\times({\omega}_{p/I}\times{{\bf
r}}_{G2/p}|_{I}).$

Combining (63) and (68) yields

$\displaystyle\\{{\bf a}\\}:=\left[\begin{array}[]{c}{\bf a}_{s/I}\\\ {\bf
a}_{b/I}\\\ {\bf a}_{p/I}\end{array}\right]$ (73) $\displaystyle=$
$\displaystyle\underbrace{\left[\begin{array}[]{cccc}{\bf 0}_{33}&{\bf
0}_{31}&{\bf 0}_{31}&{\bf I}_{3}\\\ {\bf V}_{21}&{\bf v}_{22}&{\bf
0}_{31}&{\bf I}_{3}\\\ {\bf V}_{31}&{\bf v}_{32}&{\bf v}_{33}&{\bf
I}_{3}\end{array}\right]}_{{\bf
V}}\underbrace{\left[\begin{array}[]{c}\dot{\omega_{s/I}}\\\
\dot{\sigma}_{1}\\\ \dot{\sigma}_{2}\\\ \dot{{\bf
v}}_{s/I}\end{array}\right]}_{\dot{{\bf
x}}_{g}}+\underbrace{\left[\begin{array}[]{c}{\bf 0}_{31}\\\ {{\bf
a}}_{b/I}^{r}\\\ {{\bf a}}_{p/I}^{r}\end{array}\right]}_{\\{{\bf a}_{r}\\}}$
(84) $\displaystyle:=$ $\displaystyle{\bf V}\dot{{\bf x}}_{g}+\\{{\bf
a}_{r}\\},$ (85)

which is (9b). We also showed that ${\bf V}$ defined in (85) is the same as
the one defined in (34) or in (9b).

Model (10) is a very general multi-body rigid system model. For a three-body
rigid system like LUVOIR, assume that the $3\times 3$ inertia matrices for the
spacecraft, the boom, and the payload are given as ${\bf J}_{s}$, ${\bf
J}_{b}$, and ${\bf J}_{p}$, then we have

$[{\bf J}]=\left[\begin{array}[]{ccc}{\bf J}_{s}&{\bf 0}_{33}&{\bf 0}_{33}\\\
{\bf 0}_{33}&{\bf J}_{b}&{\bf 0}_{33}\\\ {\bf 0}_{33}&{\bf 0}_{33}&{\bf
J}_{p}\end{array}\right].$ (86)

Assume that the masses of the spacecraft, the boom, and the payload are given
as $m_{s}$, $m_{b}$, and $m_{p}$, then we have

$[{\bf M}]=\left[\begin{array}[]{ccc}m_{s}{\bf I}_{3}&{\bf 0}_{33}&{\bf
0}_{33}\\\ {\bf 0}_{33}&m_{b}{\bf I}_{3}&{\bf 0}_{33}\\\ {\bf 0}_{33}&{\bf
0}_{33}&m_{p}{\bf I}_{3}\end{array}\right].$ (87)

Assume there are no external forces acting on the rigid bodies, then we have
$\\{{\bf f}\\}={\bf 0}$. Assume that the control torques ${\bf u}$ on the
spacecraft, the boom, and the payload are $\tau_{s}$, $\tau_{b}$, and
$\tau_{p}$, i.e.,

${\bf u}=[\tau_{s}^{{\rm T}},~{}\tau_{b}^{{\rm T}},~{}\tau_{p}^{{\rm
T}}]^{{\rm T}},$ (88)

then, we have

$\\{\tau\\}=\left[\begin{array}[]{c}\tau_{s}-\tau_{b}\\\ \tau_{b}-\tau_{p}\\\
\tau_{p}\end{array}\right].$ (89)

Finally, we can express $\\{\omega\times{\bf h}\\}$ in terms of the angular
rates $\omega_{s/I}$, $\omega_{b/s}$, and $\omega_{p/b}$ of the spacecraft,
the boom, and the payload as follows:

$\\{\omega\times{\bf h}\\}=\left[\begin{array}[]{c}\omega_{s/I}\times{\bf
J}_{s}\omega_{s/I}\\\ \omega_{b/s}\times{\bf J}_{b}\omega_{b/s}\\\
\omega_{p/b}\times{\bf J}_{p}\omega_{p/b}\end{array}\right].$ (90)

Let

${\bf L}=\left(\Omega^{{\rm T}}[{\bf J}]\Omega+{\bf V}^{{\rm T}}[{\bf M}]{\bf
V}\right),$ (91) ${\bf r}_{1}=\Omega^{{\rm T}}\left(\\{\tau\\}-[{\bf
J}]\\{\alpha_{r}\\}-\\{\omega\times{\bf h}\\}\right),$ (92)

and

${\bf r}_{2}={\bf V}^{{\rm T}}\left(\\{{\bf f}\\}-[{\bf M}]\\{{\bf
a}_{r}\\}\right).$ (93)

Substituting (86), (87), (89), (90), (91), (93), and $\\{{\bf f}\\}={\bf 0}$
into (10), we have the three-body rigid system model:

${\bf L}\dot{{\bf x}}_{g}={\bf r}_{1}+{\bf r}_{2},$ (94)

or

$\dot{{\bf x}}_{g}={\bf L}^{-1}({\bf r}_{1}+{\bf r}_{2}).$ (95)

Equation (95) looks very simple, but it is a nonlinear system because ${\bf
L}$, ${\bf r}_{1}$ and ${\bf r}_{2}$ have nonlinear components of the state
and control variables. We need a linear system model so that we can apply LQR
or robust pole assignment designs. First, we must rescope the model for the
purpose of pure attitude control. We take our generalized speeds ${\bf
x}_{g}$, discard the ${\bf v}_{s/I}$ component which decouples from the
attitude states, and add the spacecraft’s Euler angles $\phi$, $\theta$, and
$\psi$ in order to define our state vector ${\bf
x}=[\phi,\theta,\psi,\gamma,\lambda,\omega_{s/I},\sigma_{1},\sigma_{2}]^{{\rm
T}}$. The kinematical differential equations associated with these Euler
angles in the reference inertial frame are given as [22, Page 429, Space-three
1-2-3]:

$\displaystyle\left[\begin{array}[]{c}\dot{\phi}\\\ \dot{\theta}\\\
\dot{\psi}\end{array}\right]=\left[\begin{array}[]{c}\left<\omega_{s/I},[1,~{}\sin(\phi)\tan(\theta),~{}\cos(\phi)\tan(\theta)]^{{\rm
T}}\right>\\\ \left<\omega_{s/I},[0,~{}\cos(\phi),~{}-\sin(\phi)]^{{\rm
T}}\right>\\\
\left<\omega_{s/I},[0,~{}\sin(\phi)\sec{\theta},~{}\cos(\phi)\sec{\theta}]^{{\rm
T}}\right>\end{array}\right],$ (102)

where $\left<{\bf a},{\bf b}\right>$ denotes the inner product of two vectors
of ${\bf a}$ and ${\bf b}$. Therefore, the revised state space nonlinear
system is given as:

$\displaystyle\dot{{\bf x}}:=\left[\begin{array}[]{c}\dot{\phi}\\\
\dot{\theta}\\\ \dot{\psi}\\\ \dot{\gamma}\\\ \dot{\lambda}\\\
\dot{\omega}_{s/I}\\\ \dot{\sigma}_{1}\\\ \dot{\sigma}_{2}\\\
\end{array}\right]=\left[\begin{array}[]{c}\left<\omega_{s/I},[1,~{}\sin(\phi)\tan(\theta),~{}\cos(\phi)\tan(\theta)]^{{\rm
T}}\right>\\\ \left<\omega_{s/I},[0,~{}\cos(\phi),~{}-\sin(\phi)]^{{\rm
T}}\right>\\\
\left<\omega_{s/I},[0,~{}\sin(\phi)\sec{\theta},~{}\cos(\phi)\sec{\theta}]^{{\rm
T}}\right>\\\ \sigma_{1}\\\ \sigma_{2}\\\ \left[{\bf I}_{3},~{}~{}~{}{\bf
0}_{35}\right]({\bf L}^{-1}({\bf r}_{1}+{\bf r}_{2}))\\\
\left[0,~{}0,~{}0,~{}1,~{}0,~{}0,~{}0,~{}0\right]({\bf L}^{-1}({\bf
r}_{1}+{\bf r}_{2}))\\\ \left[0,~{}0,~{}0,~{}0,~{}1,~{}0,~{}0,~{}0\right]({\bf
L}^{-1}({\bf r}_{1}+{\bf r}_{2}))\end{array}\right].$ (119)

###### Remark 3.1

The procedure of the 3-body modeling can easily be applied to any multibody
system with tree structure, and the modeled system has the structure described
in [12]. It is also worthwhile to note that the final state space model
discarded some states in ${\bf x}_{g}$ and added some states into ${\bf x}$,
therefore, the dimensions of ${\bf x}_{g}$ and ${\bf x}$ are different.
Finally, (119) involves an analytic inverse matrix ${\bf L}^{-1}$ (its
computation will be discussed in the next section), and is slightly different
from Stoneking’s implementation (10).

##  4 Linearization and controller design

To use popular controller design methods, we need to have a linearized rigid
dynamics model.

###  4.1 Linearization

Now, we linearize the nonlinear system (119) about a desired new equilibrium
state (when this equilibrium state is attained, ${\bf u}={\bf 0}$) so that we
will have a symbolic linear system. Assume that this equilibrium state is at
${\bf
x}_{d}=[\phi_{d},~{}\theta_{d},~{}\psi_{d},~{}\gamma_{d},~{}\lambda_{d},~{}0,~{}0,~{}0,~{}0,~{}0]^{{\rm
T}}$ and the control torques are zeros, i.e., ${\bf u}={\bf 0}$. Therefore,

$\dot{{\bf x}}=\frac{\partial{\bf f}({\bf x},{\bf u})}{\partial{\bf
x}}\Bigg{|}_{\begin{array}[]{c}{\bf x}={\bf x}_{d}\\\ {\bf u}={\bf
0}\end{array}}({\bf x}-{\bf x}_{d})+\frac{\partial{\bf f}({\bf x},{\bf
u})}{\partial{\bf u}}\Bigg{|}_{\begin{array}[]{c}{\bf x}={\bf x}_{d}\\\ {\bf
u}={\bf 0}\end{array}}{\bf u}.$ (120)

where $\phi_{d}=\pi/2$, and $\theta_{d}=\psi_{d}=\gamma_{d}=\lambda_{d}=0$.

###### Remark 4.1

In this case, the target equilibrium state is a $90^{o}$ rotation of the
spacecraft in roll axis from the current state. Our simulation in the next
section will show that the designed controller works in such a large
rotational maneuver. In the next section, we will discuss a method to obtain
the analytic formula for (120).

###  4.2 Symbolic inverse for linearization

Clearly, it will be very tedious, if it is not impossible, to find the
analytic partial derivatives for (120), which involves the calculation of the
analytic partial derivatives of ${\bf L}^{-1}$. Bentz and Lewis suggested in
[14] using Matlab symbolic function ‘diff’ and the symbolic inverse function
‘inv’ for matrix ${\bf L}$. For this $8\times 8$ matrix ${\bf L}$, even using
Matlab symbolic inverse, the computation is still too complex to handle.
Fortunately, we are only interested in the first five states in (95), (see
(43)), which are the last five states in (119), We can use the method proposed
in [15]. Let

${\bf L}=\left[\begin{array}[]{cc}{\bf L}_{1}&{\bf L}_{2}\\\ {\bf L}_{2}^{{\rm
T}}&{\bf L}_{3}\end{array}\right],\hskip 7.22743pt{\bf p}={\bf r}_{1}+{\bf
r}_{2}=\left[\begin{array}[]{c}{\bf p}_{1}\\\ {\bf
p}_{2}\end{array}\right],\hskip 7.22743pt\dot{{\bf
x}}_{g}=\left[\begin{array}[]{c}\dot{{\bf x}}_{g,1}\\\ \dot{{\bf
x}}_{g,2}\end{array}\right],$

then, we have

$\left[\begin{array}[]{cc}{\bf L}_{1}&{\bf L}_{2}\\\ {\bf L}_{2}^{{\rm
T}}&{\bf L}_{3}\end{array}\right]\left[\begin{array}[]{c}\dot{{\bf
x}}_{g,1}\\\ \dot{{\bf
x}}_{g,2}\end{array}\right]=\left[\begin{array}[]{c}{\bf p}_{1}\\\ {\bf
p}_{2}\end{array}\right].$ (121)

Solving the second equation of (121) for $\dot{{\bf x}}_{g,2}$ gives

$\dot{{\bf x}}_{g,2}={\bf L}_{3}^{-1}({\bf p}_{2}-{\bf L}_{2}^{{\rm
T}}\dot{{\bf x}}_{g,1}).$ (122)

Substituting (122) into the first equation of (121) gives

$\dot{{\bf x}}_{g,1}=\left({\bf L}_{1}-{\bf L}_{2}{\bf L}_{3}^{-1}{\bf
L}_{2}^{{\rm T}}\right)^{-1}\left({\bf p}_{1}-{\bf L}_{2}{\bf L}_{3}^{-1}{\bf
p}_{2}\right),$ (123)

which involves symbolic inverses of a $3\times 3$ matrix ${\bf L}_{3}^{-1}$
and a $5\times 5$ matrix $\left({\bf L}_{1}-{\bf L}_{2}{\bf L}_{3}^{-1}{\bf
L}_{2}^{{\rm T}}\right)^{-1}$.

###  4.3 Representation of vectors in inertial frame

All constants in the multi-body model are provided by mechanical engineers
according to the spacecraft designs. Some constants in the multi-body model
are independent to the frames, such as mass of spacecraft, mass of payload,
etc., but some constants are dependent to the frames. Most likely, the
distance vectors in a rigid body are given in that rigid body frame, but we
need to represent these distance vectors in the inertial frame in the model
(123) as discussed in the previous section.

Let ${\bf r}_{G1/s}|_{s}=a_{1}{\bf x}_{s}+a_{2}{\bf y}_{s}+a_{3}{\bf z}_{s}$
be the position vector from the center of mass of the spacecraft pointing to
the joint $G_{1}$ represented in the spacecraft frame, ${\bf
r}_{G1/b}|_{b}=b_{1}{\bf x}_{b}+b_{2}{\bf y}_{b}+b_{3}{\bf z}_{b}$ be the
position vector from the center of mass of the boom pointing to the joint
$G_{1}$ represented in the boom frame, ${\bf r}_{G2/b}|_{b}=c_{1}{\bf
x}_{b}+c_{2}{\bf y}_{b}+c_{3}{\bf z}_{b}$ be the position vector from the
center of mass of the boom pointing to the joint $G_{2}$ represented in the
boom frame, and ${\bf r}_{G2/p}|_{p}=d_{1}{\bf x}_{p}+d_{2}{\bf
y}_{p}+d_{3}{\bf z}_{p}$ be the position vector from the center of mass of the
payload pointing to the joint $G_{2}$ represented in the payload frame. Denote
${\bf r}_{G1/s}|_{I}$ be the the position vector from the center of mass of
the spacecraft pointing to the joint $G_{1}$ represented in the inertial
frame, ${\bf r}_{G1/b}|_{I}$ be the the position vector from the center of
mass of the boom pointing to the joint $G_{1}$ represented in the inertial
frame, ${\bf r}_{G2/b}|_{I}$ be the the position vector from the center of
mass of the boom pointing to the joint $G_{2}$ represented in the inertial
frame, and ${\bf r}_{G2/p}|_{I}$ be the the position vector from the center of
mass of the payload pointing to the joint $G_{2}$ represented in the inertial
frame, then using (12), (13), (14), (15), and (16), we have

$\displaystyle{\bf r}_{G1/s}|_{I}$ $\displaystyle=$
$\displaystyle\mathcal{O}_{I/s}{\bf
r}_{G1/s}|_{s}=\mathcal{O}_{I/s}[a_{1},a_{2},a_{3}]^{{\rm T}}$
$\displaystyle=$ $\displaystyle\mathcal{O}_{s/I}^{{\rm
T}}[a_{1},a_{2},a_{3}]^{{\rm T}},$ (124a) $\displaystyle{\bf r}_{G1/b}|_{I}$
$\displaystyle=$ $\displaystyle\mathcal{O}_{I/b}[b_{1},b_{2},b_{3}]^{{\rm
T}}=\mathcal{O}_{b/I}^{{\rm T}}[b_{1},b_{2},b_{3}]^{{\rm T}}$ $\displaystyle=$
$\displaystyle\left(\mathcal{O}_{b/s}\mathcal{O}_{s/I}\right)^{{\rm
T}}[b_{1},b_{2},b_{3}]^{{\rm T}},$ (124b) $\displaystyle{\bf r}_{G2/b}|_{I}$
$\displaystyle=$ $\displaystyle\mathcal{O}_{I/b}[c_{1},c_{2},c_{3}]^{{\rm
T}}=\mathcal{O}_{b/I}^{{\rm T}}[c_{1},c_{2},c_{3}]^{{\rm T}}$ $\displaystyle=$
$\displaystyle\left(\mathcal{O}_{b/s}\mathcal{O}_{s/I}\right)^{{\rm
T}}[c_{1},c_{2},c_{3}]^{{\rm T}},$ (124c) $\displaystyle{\bf r}_{G2/p}|_{I}$
$\displaystyle=$ $\displaystyle\mathcal{O}_{I/p}[d_{1},d_{2},d_{3}]^{{\rm
T}}=\mathcal{O}_{p/I}^{{\rm T}}[d_{1},d_{2},d_{3}]^{{\rm T}}$ $\displaystyle=$
$\displaystyle\left(\mathcal{O}_{p/s}\mathcal{O}_{s/I}\right)^{{\rm
T}}[d_{1},d_{2},d_{3}]^{{\rm T}},$ (124d) $\displaystyle{\bf r}_{G1/p}|_{I}$
$\displaystyle=$ $\displaystyle{\bf r}_{G2/p}|_{I}-{\bf r}_{G2/b}|_{I}+{\bf
r}_{G1/b}|_{I},$ (124e) $\displaystyle{\bf r}_{s/b}|_{I}$ $\displaystyle=$
$\displaystyle{\bf r}_{G1/b}|_{I}-{\bf r}_{G1/s}|_{I},$ (124f)
$\displaystyle{\bf r}_{s/p}|_{I}$ $\displaystyle=$ $\displaystyle{\bf
r}_{G1/p}|_{I}-{\bf r}_{G1/s}|_{I}.$ (124g)

Using the definition of (25), we can write

${\bf
r}_{s/b}|_{I}^{\times}=\left[\begin{array}[]{ccc}0&-r_{{s/b}_{3}}&r_{{s/b}_{2}}\\\
r_{{s/b}_{3}}&0&-r_{{s/b}_{1}}\\\
-r_{{s/b}_{2}}&r_{{s/b}_{1}}&0\end{array}\right],$ (125) ${\bf
r}_{s/p}|_{I}^{\times}=\left[\begin{array}[]{ccc}0&-r_{{s/p}_{3}}&r_{{s/p}_{2}}\\\
r_{{s/p}_{3}}&0&-r_{{s/p}_{1}}\\\
-r_{{s/p}_{2}}&r_{{s/p}_{1}}&0\end{array}\right],$ (126)

and similarly, we can define ${\bf r}_{G1/b}|_{I}^{\times}$, ${\bf
r}_{G1/p}|_{I}^{\times}$, and ${\bf r}_{G2/p}|_{I}^{\times}$.

Using the parameters of the LUVOIR telescope, we use a Matlab code (which is
provided in [15]) to generate the rigid linearized time invariant model

$\dot{{\bf x}}={\bf A}{\bf x}+{\bf B}{\bf u}$ (127)

with ${\bf A}$ and ${\bf B}$ given as follows:

${\bf A}=\left[\begin{array}[]{rrrrrrrrrr}0&0&0&0&0&1&0&0&0&0\\\
0&0&0&0&0&0&0&-1&0&0\\\ 0&0&0&0&0&0&1&0&0&0\\\ 0&0&0&0&0&0&0&0&1&0\\\
0&0&0&0&0&0&0&0&0&1\\\ 0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0&0&0\end{array}\right],$ ${\bf
B}=10^{-5}\left[\begin{array}[]{rrrrr}0&0&0&0&0\\\ 0&0&0&0&0\\\ 0&0&0&0&0\\\
0&0&0&0&0\\\ 0&0&0&0&0\\\ 0.027838&0.000983&0.011040&-0.001203&-0.000149\\\
0.000983&0.337771&0.007301&-0.473521&0.142849\\\
0.011040&0.007301&0.082504&-0.009908&0.002876\\\
-0.001203&-0.473521&-0.009908&0.835383&-0.391328\\\
-0.000149&0.142849&0.002876&-0.391328&0.331638\end{array}\right].$

###### Remark 4.2

The correctness of the rigid linearized model is indirectly validated when the
controller designed by this rigid model stabilizes a separately developed
flexible telescope model.

###  4.4 LQR and robust pole assignment designs

Using the linearized model obtained from Stoneking’s form of Kane’s method, we
consider two well-known state-space controller design methods: LQR and robust
pole assignment. For the linearized control system, the LQR design is to find
an optimal feedback gain matrix ${\bf K}_{LQR}$ to minimize the following cost
function [23]:

$J=\frac{1}{2}\int_{0}^{\infty}({\bf x}^{{\rm T}}{\bf Q}{\bf x}+{\bf u}^{{\rm
T}}{\bf R}{\bf u})dt$ (128)

under the state space system constraint $\dot{{\bf x}}={\bf A}{\bf x}+{\bf
B}{\bf u}$ with ${\bf A}$ and ${\bf B}$ being given in the previous section;
while the robust pole assignment design is to find an optimal feedback gain
matrix ${\bf K}_{rpa}$ such that (a) the close-loop eigenvalues of $({\bf
A}-{\bf B}{\bf K}_{rpa})$ are in the desired locations for the rigid model,
and (b) the sensitivity to the modeling uncertainty (because of using less
accurate rigid model in controller design) of the close-loop eigenvalues of
$({\bf A}-{\bf B}{\bf K}_{rpa})$ is minimized [24, 25]. Let $\Lambda={\rm
diag}(\lambda_{1},\ldots,\lambda_{n})$ and ${\bf X}=[{\bf x}_{1},{\bf
x}_{2},\ldots,{\bf x}_{n}]$ be the closed-loop diagonal eigenvalue matrix and
the corresponding eigenvector matrix of $({\bf A}-{\bf B}{\bf K}_{rpa})$, the
object of the robust pole assignment design is to solve the following
optimization problem [24, 25]:

$\displaystyle\min$ $\displaystyle\frac{1}{2}\det({\bf X}^{H}{\bf X})$
$\displaystyle s.t.$ $\displaystyle({\bf A}-{\bf B}{\bf K}_{rpa}){\bf X}={\bf
X}\Lambda$ $\displaystyle{\bf x}_{i}^{H}{\bf x}_{i}=1,\hskip
7.22743pti=1,\ldots,n,$

where the superscript $H$ is used for complex-conjugate transpose, it reduces
to a transpose if all elements of $\Lambda$ are real. A very efficient
algorithm is developed to solve this problem in [25]. A Matlab code that
implements the algorithm of [25] is available on the website of [26]. For more
details on robust pole assignment design discussed in this paper, the readers
are referred to [27, 25]. A concise description of the robust pole assignment
is available in [28, Appendix C].

The LQR design has been widely used in aerospace applications because it is
considered a good choice when energy consumption is a major consideration.
Pole assignment design is not as popular as the LQR design in this case
because it is not clear whether the design will consume more energy than the
LQR design. However, users have noticed that robust pole assignment design
[25] normally generates a small feedback gain matrix which is a good sign of
efficient use of energy. There are two other merits associated with the robust
pole assignment approach. First, the performance of the closed-loop system is
robust to modeling errors. This is important because the high fidelity model
of LUVOIR used in testing will include flexible modes that are ignored during
control system design due to modeling complexity. Second, robust pole
assignment can predict approximations of closed-loop system performance
characteristics such as settling time, oscillation frequency, etc., by
assigning the closed-loop poles in desired areas [29]. For example, to avoid
the oscillations for a second order system, all poles should be assigned to be
real according to [29, Chapter 5]. For higher order systems, the performance
is determined by dominate poles which are closer to the imaginary axis,
therefore, the dominate poles should be assigned to be real. LQR cannot do
this. Our strategy is to use LQR approach as an effective first design step
that informs the selection of the real eigenvalues for robust pole assignment
such that these poles are close to the real parts of the closed-loop
eigenvalues of LQR.

###  4.5 Simulation testing on rigid model

Figure 3: LQR and robust pole assignment design comparison for rigid model:
(a) $x_{1}$ initial state response (b) $x_{2}$ initial state response.

Figure 4: LQR and robust pole assignment design comparison for rigid model:
(a) $x_{3}$ initial state response (b) $x_{4}$ initial state response.

Figure 5: LQR and robust pole assignment design comparison for rigid model:
(a) $x_{5}$ initial state response (b) $x_{6}$ initial state response.

We calculated state feedback gain matrices for both LQR and robust pole
assignment designs. Then, we compared the system performances by analyzing
initial state responses and energy consumption. The ${\bf Q}$ and ${\bf R}$
matrices we used in LQR design are exactly the same as the ones used in [15].
The target of robust pole assignment is to have a similar settling time to the
LQR design but with fewer oscillations—this is of particular importance to
space telescopes with precision pointing requirements. This can be achieved by
choosing the prescribed closed-loop eigenvalues of the robust pole assignment
design to have similar real parts to that of the LQR design, i.e., placing all
closed-loop eigenvalues on the real axis of the complex plane.

For LQR design, the ${\bf Q}$ and ${\bf R}$ matrices are selected exactly the
same as the ones in [15]:

${\bf Q}=\left[\begin{array}[]{cc}1000{\bf I}_{8}&{\bf 0}_{82}\\\ {\bf
0}_{28}&2000{\bf I}_{2}\end{array}\right],\hskip 7.22743pt{\bf R}={\bf
I}_{5}.$ (130)

For robust pole assignment design, the prescribed close-loop eigenvalues are
selected as

$(-0.0141,-0.0135,-0.0059,-0.0058,-0.0037,-0.0036,-0.0029,-0.0028,-0.00265,-0.00262)$

which are close to the real parts of the eigenvalues of $({\bf A}-{\bf B}{\bf
K}_{LQR})$. The feedback gain matrices of the LQR and the robust pole
assignment are obtained by a Matlab functions lqr and robpole respectively.

$\displaystyle{\bf
K}_{LQR}=10^{4}\left[\begin{array}[]{rrrrrrrrrr}0.0031&0.0000&0.0000&-0.0000&0.0000&1.5433&0.0030&-0.1319&0.0045&0.0053\\\
-0.0000&0.0000&0.0031&-0.0000&0.0000&0.0030&1.2432&-0.0113&0.7782&0.4049\\\
0.0000&-0.0031&0.0000&-0.0000&0.0000&-0.1319&-0.0113&0.8907&0.0019&0.0000\\\
0.0000&-0.0000&0.0000&0.0031&-0.0000&0.0045&0.7782&0.0019&0.9157&0.5169\\\
-0.0000&0.0000&-0.0000&0.0000&0.0031&0.0053&0.4049&0.0000&0.5169&0.8071\end{array}\right]$
(136) $\displaystyle{\bf
K}_{rpa}=10^{4}\left[\begin{array}[]{rrrrrrrrrr}-0.0088&-0.0013&0.0067&-0.0000&-0.0025&-3.8298&1.8600&0.5231&0.0202&-0.6649\\\
0.0052&0.0005&-0.0112&-0.0025&-0.0042&1.4826&-4.5333&-0.0981&-1.8008&-2.2176\\\
0.0013&0.0022&-0.0010&-0.0002&0.0008&0.5224&-0.2667&-1.1273&-0.0748&0.2379\\\
0.0040&0.0004&-0.0089&-0.0023&-0.0044&1.1578&-3.6498&-0.0756&-1.6557&-2.2215\\\
0.0023&0.0001&-0.0056&-0.0015&-0.0045&0.6759&-2.3400&-0.0008&-1.1043&-2.1248\end{array}\right]$
(142)

We demonstrated that the controllers stabilize both the rigid linearized time
invariant model and the flexible LUVOIR telescope model via simulation using
Matlab and Simulink.

####  4.5.1 Oscillation comparison of the two designs

Figures $\ref{figure2}-\ref{figure10}$ compare the three-body rigid linearized
model initial state responses of the LUVOIR telescope for LQR and robust pole
assignment designs. It is clear that the initial state responses of robust
pole assignment design have fewer oscillations in general for all $10$ states
(indicating more stable pointing). If we amplify the figures, for the LQR
design, the oscillation still can be seen after $5000$ seconds, but robust
pole assignment initial state responses do not have this kind of long term
oscillations. This meets our expectation as discussed earlier. Since both
controller designs are based on rigid model, when the controllers are applied
to the rigid model, we don’t see the jitter that will be seen when the
controllers are applied to flexible model. Figure 8 depicts the long term
oscillation effects of $x_{9}$ and $x_{10}$. The osculations will adversely
affect the telescope pointing because disturbances can occur at any time for
many different reasons.

####  4.5.2 Energy consumption comparison of the two designs

Using the least energy consumption to achieve the desired performance is
always an important design consideration in space missions. Therefore, we
compare the energy consumption of the two designs. For both LQR and robust
pole assignment designs, the energy consumption can be measured by

$\int_{0}^{\infty}\|{\bf u}(t)\|dt,$ (143)

where $\|\cdot\|$ denotes the Euclidean norm. Using formula (143), we get
$\int_{0}^{10000}\|{\bf u}_{LQR}\|dt=1.45506e+05$ and $\int_{0}^{10000}\|{\bf
u}_{rpa}\|dt=1.28432e+05$, which shows that robust pole assignment consumes
noticeably less energy.

Figure 6: LQR and robust pole assignment design comparison for rigid model:
(a) $x_{7}$ initial state response (b) $x_{8}$ initial state response.

Figure 7: LQR and robust pole assignment design comparison for rigid model:
(a) $x_{9}$ initial state response (b) $x_{10}$ initial state response.

Figure 8: LQR and robust pole assignment design comparison for rigid model:
(a) $x_{9}$ initial state response (b) $x_{10}$ initial state response.

###  4.6 Simulation testing on the flexible model

Figure 9: LQR and robust pole assignment design comparison for flexible model:
(a) Roll angle initial state response (b) Pitch angle initial state response.

Figure 10: LQR and robust pole assignment design comparison for flexible
model: (a) Yaw angle initial state response (b) Roll angular rate initial
state response.

Figure 11: LQR and robust pole assignment design comparison for flexible
model: (a) Pitch angular rate initial state response (b) Yaw angular rate
initial state response.

Figure 12: LQR and robust pole assignment design comparison for flexible
model: (a) Pitch augular rate initial state response (b) Yaw angular rate
initial state response.

Figure 13: LQR and robust pole assignment design comparison for flexible
model: (a) Gimbal 1 torque initial state response (b) Gimbal 2 torque initial
state response.

The LUVOIR flexible model has $20$ modes on the spacecraft, $9$ on the boom,
and $100$ on the payload [15], the system has about $260$ states and most of
these states are not measurable. To have an attainable design and a seamless
implementation, we designed the controllers based on the coarse rigid model.
Since the rigid model approximates the flexible model, we need to validate the
designs by using the simulation for the high fidelity flexible model and
examining the performances. Modeling a flexible mechanical system has been
discussed in [30] and a Simulink implementation for the LUVOIR telescope was
described in [15, 31]. As mentioned in the introduction section, if LQR and/or
robust pole assignment designs stabilize the rigid model but cannot stabilize
the high fidelity flexible model, then, redesigns are necessary. As a matter
of fact, the feedback gain matrices given in the previous section are the ones
of the final design which are obtained after a few iterations.

For the designs given in the previous section, the total energy consumption is
$1.5448e+03$ for the LQR design and is $1.2683e+03$ for the robust pole
assignment design. Again, the energy consumption for the robust pole
assignment design is slightly less than the one for the LQR design.

The performances of these two designs are compared and displayed in Figures
9-13. The LQR design settles the spacecraft faster than the robust pole
assignment design, which is good. It can also be seen, from Figures 13, that
torques requested for LQR design have oscillations at about 0.45 Hz, which is
not good because not only does it consume more energy, but it may also
introduce jitters to the telescope. Figures in 12 (a) and (b) are the
amplified pitch and yaw angular rate of the spacecraft body, which shows the
0.45 HZ oscillations after 200 minutes.

We summarize the systematic space telescope design methodology in the
following procedure:

  * 1.

Develop a rigid symbolic nonlinear multibody model using Stoneking’s form of
Kane’s equation (10).

  * 2.

Take symbolic inverse for Kane’s model to obtain the symbolic nonlinear state
space model $\dot{{\bf x}}_{g}={\bf f}_{g}({\bf x},{\bf u})$.

  * 3.

Determine spacecraft kinematical differential equations associated with the
spacecraft Euler angle using the method provided in [22].

  * 4.

Determine rotary angular dynamics.

  * 5.

Extract relavent states from the symbolic nonlinear state space model
$\dot{{\bf x}}_{g}={\bf f}_{g}({\bf x},{\bf u})$.

  * 6.

Combine states obtained in Steps 3, 4, and 5 to form a rigid nonlinear
symbolic state space model $\dot{{\bf x}}={\bf f}({\bf x},{\bf u})$.

  * 7.

Symbolically linearize the nonlinear system about the desired equilibium point
to get a symbolic rigid linear system model.

  * 8.

Using the spacecraft parameters to populate the symbolic model to get the
spacecraft specific rigid linear system model.

  * 9.

Design a LQR controller which stabilizes both the rigid linear system model
and the flexible system model (developed separately from the rigid linear
system model).

  * 10.

Design a robust pole assignment controller, whose desired closed-loop poles
are all real and the value of the real poles are close to the real part of the
closed loop poles of LQR design, such that the robust pole assignment design
stabilizes both the rigid linear system model and the flexible system model.

##  5 Conclusion

In this paper, we presented a modeling method for multi-body system using
Kane’s method. A rigid model for LUVOIR telescope is established as a result.
LQR and robust pole assignment methods are used to design the controllers
using the linearized rigid model. Simulation test of the closed loop system
using both rigid and flexible models are performed. The test result based on
the rigid linearized time invariant model shows that robust pole assignment
has better performance in terms of energy consumption and pointing accuracy
(measured by the low frequency oscillation around the equilibrium point). For
the test on the flexible Simulink model (which was developed by Roger Chen
[31]), LQR design has a shorter (better) settling time but the gimbal commands
have oscillations at about 0.45 Hz which may cause the jitters problem and
affect the image quality of the telescope; while robust pole assignment design
has a similar gimbal command oscillation problem at the beginning, it
attenuates fast to zero. Overall, we recommend the robust pole assignment
technique for this application with the caveat that the LQR approach is
effective as a method of jump-starting the robust pole assignment design. That
is, the LQR approach allows the designer to make an initial pass in which
setting time and other performance characteristics are tuned through the
traditional cost function weight matrices. This provides a set of desired real
eigenvalue components that can be targeted through robust pole assignment in
order to refine the performance, e.g., in order to improve damping.

##  6 Acknowledgments

The authors are grateful to the anonymous reviewers for their constructive
comments which are invaluable in the preparation of the final version. The
first author thanks Eric Stoneking and Gary Mosier of Goddard Space Flight
Center at NASA for introducing him into this interesting project.

## References

  * [1] D. Fischer, and et. al., LUVOIR final report, NASA, Washington D.C., 2021, Chap. 11, pp. 1–24.
  * [2] P. Eberhard and W. Schiehlen, “Computational dynamics of multibody system: history, formulism, and application,” Journal of Computational and Nonlinear Dynamics, Vol. 1, pp. 3–12, 2006.
  * [3] P. E. B. Jourdain, “Note on an analogue of Gauss’ principle of least constraint,” The Quarterly Journal of Pure and Applied Mathematics, Vol. XL, pp. 153–157, 1909.
  * [4] T. R. Kane, “Dynamics of nonholonomic systems,” Transactions of the ASME Journal of Applied Mechanics, Vol. 28, pp. 574–578, 1961.
  * [5] K. W. Buffinton, “Kane’s Method in Robotics,” in Robotics and Automation Handbook, CRC Press, New York, 2005, Chap. 6, pp. 89–119.
  * [6] M. A. Chace and Y. O. Bayazitoglu, “Development and application of a generalized D’Alembert force for multifreedom mechanical system”, J. Manuf. Sci. Eng., Vol. 93, pp. 317–326, 1971.
  * [7] T. R. Kane and D. A. Levinson, Dynamics Theory and Applications, McGraw-Hill Book Company, Boston, 1985, Chap. 2, pp. 15–57.
  * [8] W. Schiehlen, Multibody Systems Handbook, Springer, Berlin, 1990, Chap. 1, pp. 1–9.
  * [9] E. A. Desloge, “Development and application of a generalized D’Alembert force for multifreedom mechanical system,” Journal of Guidance, Control, and Dynamics, Vol. 10, pp. 120–122, 1987.
  * [10] D. A. Levinson and A. K. Banerjee, “Comment on “Efficiency of the Gibbs-Appell equations”,” Journal of Guidance, Control, and Dynamics, Vol. 13, pp. 381–382, 1990.
  * [11] J.-C. Piedboeuf, “Kane’s equations or Jourdain’s principle?” in Proceedings of 36th Midwest Symposium on Circuits and Systems, IEEE, Detroit, MI, USA, 1993, pp. 1471–1474.
  * [12] E. Stoneking, Implementation of Kane’s method for a spacecraft composed of multiple rigid bodies, in AIAA Guidance, Navigation, and Control (GNC) Conference, Detroit, MI, USA, 2013, pp. 46–49.
  * [13] E. Stoneking, “42: A General-Purpose Spacecraft Simulation,” NASA Software Designation GSC-16720-1, 2010-, [Online]. Available: https://sourceforge.net/projects/fortytwospacecraftsimulation, 2010.
  * [14] W. Bentz and L. Sacks, “Derivation of equations of motion and model-based control for 3D rigid body approximation of LUVIOR,” NASA technical note, 2018. Available on https://ntrs.nasa.gov/citations/20240000697
  * [15] W. Bentz and L. Sacks, “Summary of model-based attitude control of LUVOIR in modular dynamics analysis (MDA) simulink environment,” NASA technical report, 2018. Available on https://ntrs.nasa.gov/citations/20240000711.
  * [16] J. Chae and T. Park, “Dynamic Modeling and Control of Flexible Space Structures,” KSME International Journal, vol. 17 No. 12, 2003, pp. 1912–1921.
  * [17] M. Liu, D. Cao, and D. Zhu, “Equivalent Dynamic Model of the Space Antenna Truss with Initial Stress,” AIAA Journal, Vol. 58, No. 4, 2020, pp. 1851–1863.
  * [18] W. He, and S. S. Ge, “Dynamic Modeling and Vibration Control of a Flexible Satellite,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 51, No. 2, 2015, pp. 1422–1431.
  * [19] Y. Li, S. Li, and M. Xin, “Dynamic Modeling and Attitude Control of Large-Scale Flexible Parallel Multibody Spacecraft,” Journal of Guidance, Control, and Dynamics, Volume 45, Issue 12, 2022, pp. 2304-2317.
  * [20] Q. Hu, Y. H. Jia, and S. J. Xu, “A new computer-oriented approach with efficient variables for multibody dynamics with motion constraints,” Acta Astronautica, Vol. 81, p. 380–389, 2012.
  * [21] R. P. Singh, R. J. V. Voort, and P.W. Likins, “Dynamics of flexible bBodies in tree topology—a computer-oriented approach,” Journal of Guidance, Vol. 8, pp. 584–590, 1985.
  * [22] T. R. Kane, P. W. Likins, and D. A. Levinson, Spacecraft Dynamics, McGraw-Hill Book Company, Boston, 1993, Chap. Appendix, p. 429.
  * [23] F. L. Lewis, D. Vrabie, and V. L. Syrmos, Optimal Control, John Wiley & Sons, Inc., New Jersey, 2012, Chap. 2, pp. 19–53.
  * [24] Y. Yang, “A globally convergent algorithm for robust pole assignment,” Science in China (Series A), Vol. 23, pp. 1126–1131, 1993.
  * [25] A. L. Tits, and Y. Yang, “Globally convergent algorithms for robust pole assignment by state feedback,” IEEE transactions on Automatic Control, Vol. 41, pp. 1432–1452, 1996.
  * [26] Y. Yang, Robpole, Matlab file exchange, [Online]. Available: https://www.mathworks.com/matlabcentral/fileexchange/53969-robpole, 2022.
  * [27] J. Kautsky, N. K. Nichols, and P. Van Dooren, Robust pole assignment in linear state feedback, International Journal of Control, Vol. 41, pp. 1129-1155, 1985.
  * [28] Y. Yang, Spacecraft Modeling, Attitude Determination, and Control:: Quaternion-based Approach, CRC Press, 2019.
  * [29] R. C. Dorf, and R. H. Bishop, Modern Control Systems, Pearson Printice Hall, New Jersey, 2008, Chaps. 5.4,5.5, pp. 287–295.
  * [30] L. Meirovitch, Fundamentals of Vibrations, McGraw-Hill, Boston, 2001, Chap. 1, pp. 1–79.
  * [31] R. Chen, “Modular dynamic analysis,” NASA unpublished technical note, 2016.